Axisymmetric Flow Driven by a Pressure Gradient in a Constricted Pipe

Compute and visualize transient axisymmetric viscous flow in a pipe with a narrowed middle section under a time-dependent pressure gradient

A viscous, incompressible fluid flows through a long cylindrical pipe that becomes narrower in the middle. The motion is caused by a pressure difference applied between the inlet and the outlet. The goal is to compute how the velocity and pressure evolve in time and to examine how the narrowed section affects the velocity profile at the outlet. Since both the pipe geometry and the loading are rotationally symmetric about the central axis, the problem can be described using an axisymmetric formulation.

Define the final simulation time (tEnd), the maximum pipe radius (R), and the pipe length (L):

In[1]:=

tEnd=5;R=3/2;L=10;

Let the pipe cross-section Ω be a rectangle with a circular constriction removed:

In[2]:=

Ω=BooleanRegion[#1 && !#2 & , {Rectangle[{0, 0}, {3/2, 10}], Disk[{3/2, 5}, 3/4]}]

Out[2]=

The flow is modeled by the incompressible Navier-Stokes equations ρ(

∂

+ (v · ∇)) v) = -∇p + μ

∇

v with ∇ · v = 0, where v = (u, w) denotes the radial and axial velocity components and p is the pressure. The system is solved in an axisymmetric (r, z) cross-section of the pipe region Ω, meaning the solution is independent of the angular coordinate:

In[3]:=

pdes=FluidFlowPDEComponent[{{u[t,r,z],0,w[t,r,z],p[t,r,z]},t,{r,θ,z}},"DynamicViscosity"1/10,"MassDensity"1,"RegionSymmetry""Axisymmetric"]

Out[3]=

{u[t,r,z],w[t,r,z]}.

∇

{r,z}

u[t,r,z]+

(0,1,0)

[t,r,z]+

u[t,r,z]

(0,1,0)

[t,r,z]

-

(0,0,2)

[t,r,z]-

(0,1,1)

[t,r,z]+



(0,0,2)

[t,r,z]+

(0,1,1)

[t,r,z]-

(0,2,0)

[t,r,z]+

(1,0,0)

[t,r,z],{u[t,r,z],w[t,r,z]}.

∇

{r,z}

w[t,r,z]+

(0,0,1)

[t,r,z]-

(0,0,2)

[t,r,z]+

(0,0,1)

[t,r,z]-

(0,1,0)

[t,r,z]

10r

-

(0,1,1)

[t,r,z]-

(0,2,0)

[t,r,z]+



(0,1,1)

[t,r,z]+

(0,2,0)

[t,r,z]+

(1,0,0)

[t,r,z],

u[t,r,z]

(0,0,1)

[t,r,z]+

(0,1,0)

[t,r,z]

No-slip conditions enforce zero velocity on the pipe wall, symmetry sets the radial velocity to zero on the axis, inlet pressure is ramped in time using

Min

[t,1]

to avoid an abrupt start, outlet pressure is fixed to zero as a reference, and the fluid is initially at rest:

In[4]:=

AbsoluteTiming[{rVel,zVel,pressure}=NDSolveValue[{pdes{0,0,0},{DirichletCondition[u[t,r,z]0,z0],DirichletCondition[u[t,r,z]0,r0],DirichletCondition[{u[t,r,z]0,w[t,r,z]0},r>0&&z≠0&&z≠L],DirichletCondition[p[t,r,z]Min[t,1]*3/2*L,z0],DirichletCondition[{p[t,r,z]0,u[t,r,z]0},zL]},{u[0,r,z]0,w[0,r,z]0,p[0,r,z]0}},{u,w,p},{r,z}∈Ω,{t,0,tEnd},Method{"PDEDiscretization"{"MethodOfLines","SpatialDiscretization"{"FiniteElement","InterpolationOrder"{u2,w2,p1}}}}];]

Out[4]=

{6.37792,Null}

Visualize the velocity for various times:

In[5]:=

Ωfull=DiscretizeRegion@RegionUnion[RegionDifference[Rectangle[{0,0},{R,L}],Disk[{R,L/2},R/2]],TransformedRegion[RegionDifference[Rectangle[{0,0},{R,L}],Disk[{R,L/2},R/2]],ScalingTransform[{-1,1}]]];

In[6]:=

frames=VectorPlot[{Sign[r]rVel[#,Abs[r],z],zVel[#,Abs[r],z]},{r,z}∈Ωfull,AspectRatioAutomatic,PlotLegendsBarLegend[{Automatic,{0,6}}],ColorFunctionColorData[{"WL12DefaultVectorDensityGradient",{0,6}}],ColorFunctionScalingFalse,PlotLabel("Time: "<>ToString[#])<>" s",FrameFalse,VectorPointsFine,VectorScalingAutomatic]&/@Range[1,tEnd,0.2];frames=Rasterize[#,"Image",RasterSize400,ImageSize300]&/@frames;ListAnimate[frames,ControlPlacementTop,AnimationRepetitions2]

Out[6]=