# Message Boards

Answer
(Unmark)

GROUPS:

6

# [WSS22] Fluid flow analysis through squared arrangements of pipes

Posted 23 days ago

Fluid flow analysis through squared arrangements of pipes by Giuseppe Parasiliti Rantone Institut Jean Le Rond d’Alembert, Sorbonne Université and the company Gaztransport & Technigaz (GTT) The main goal is to study fluid flow through a closed arrangement of pipes using the finite elements method. The equations we take into consideration are Navier-Stokes equations eventually coupled with the energy conservation equation. We study the behaviour of laminar flows on different types of meshes, we start from a simple two dimensional mesh, we then extend the analysis to the Sierpinski carpet and the Menger Sponge. Sierpinski carpet is a two dimensional mesh that iteration by iteration assume fractal patterns, the Menger Sponge is its three dimensional extension. We deal with different boundary conditions, we impose both pressure and temperature gradients between the walls. We study both stationary and unsteady problems. We find several flow patterns and we understand how laminar flow behaves on fractal meshes and how boundary conditions are crucial to this behaviour. Next steps will be to extend the analysis to turbulent flows and to make an unsteady analysis of the stationary flows we considered.
Introduction Fluid flows are natural processes we experiment in our everyday life and despite that, their features are still partly unknown. Fluid motion is the result of viscous forces, pressure variations and gravitational forces. Given some initial condition and some boundary constraints, the equations fluids are governed by are called Navier-Stokes equations. They consist in the mass and the momentum conservation laws : ∂ t ∂ t u p g σ ρ C p ∂ t ∂ t C p α D(u) k
Pressure-controlled laminar flows Let us construct a symmetric squared arrangement of pipes. Notice that the domain is constructed as a Region and that the function we decided to use for refinement is MaxBoundaryCellMeasure. L=3; x 1 y 1 β 1 x 1 y 1 x 1 y 1 β 2 x 1 y 1 x 1 y 1 β 3 x 1 y 1 x 1 y 1 β 4 x 1 y 1 x 1 y 1 β 1 β 2 β 3 β 4 Let us assume that the density of the flow is constant and that gravity could be neglected. The Navier-Stokes equations for two dimensional domains reduce to the following: ρ( ∂ t ∂ x ∂ y ∂ x ∂ x,x ∂ y,y ∂ t ∂ x ∂ y ∂ y ∂ x,x ∂ y,y ∂ x ∂ y μ ρ=1 μ= -4 10 ρ n+1 u n u Δt n u ∂ x n+1 u n v ∂ y n+1 u ∂ x n+1 p ∂ x,x n+1 u ∂ y,y n+1 u n+1 v n v Δt n u ∂ x n+1 v n v ∂ y n+1 v ∂ y n+1 p ∂ x,x n+1 v ∂ y,y n+1 v ∂ x n+1 u ∂ y n+1 v In[]:= dt=0.02;op1[i_]={ρ((u[x,y]-xVel[i-1][x,y])/dt+xVel[i-1][x,y] ∂ x ∂ y ∂ x,x ∂ y,y ∂ x ∂ x ∂ y ∂ x,x ∂ y,y ∂ y ∂ x ∂ y -4 10 Let us give time independent boundary conditions for closed domains by imposing a pressure gradient between the walls and no-slip conditions for velocity on all boundary edges: boundaryconditions1={DirichletCondition[{u[x,y]==0.,v[x,y]==0.},True],DirichletCondition[p[x,y]==0,x==L],DirichletCondition[p[x,y]==1,x==0]}; Let us give some initial conditions and create a time loop in which at each time step we solve Navier-Stokes equations in space on the squared arrangement of pipes. We will use a finite elements method with a quadratic order of interpolation for velocity and a linear order for pressure. xVel[0][x_,y_]:=0;yVel[0][x_,y_]:=0;pressure[0][x_,y_]:=0;Table[{xVel[i],yVel[i],pressure[i]}=NDSolveValue[{op1[i],boundaryconditions1},{u,v,p},{x,y}∈mesh,Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,p->1}}],{i,1,15}]; Let us show the flow: In[]:= StreamPlot[{xVel[currenttime][x,y],yVel[currenttime][x,y]},{x,y}∈xVel[currenttime]["ElementMesh"],AspectRatio->Automatic,PlotRange->All,PlotLegendsAutomatic,ImageSize->Medium] From the StreamPlot we see that the horizontal velocity is practically zero in the vertical pipes and that the flow is completely horizontal. Moreover the fluid flows in opposite direction of the gradient of pressure accordingly to what could be predictable a priori. Let us show the values of Reynolds number in both x and y directions to confirm that we are in laminar flow. {minX,maxX}=MinMax[xVel[currenttime]["ValuesOnGrid"]];{minY,maxY}=MinMax[yVel[currenttime]["ValuesOnGrid"]];characteristicUx=(Abs[minX]+Abs[maxX])/2;characteristicUy=(Abs[minY]+Abs[maxY])/2;ReynoldsX={ρUL/μ}/.{μ-> -4 10 -4 10 Out[]= {2981.81} Out[]= {5711.32} Let us now deal with more complex boundary conditions, we will impose the pressure to be 1 0 boundaryconditionsdiagonalflow={DirichletCondition[{u[x,y]==0.,v[x,y]==0.},True],DirichletCondition[p[x,y]==0,x==L&&y<L],DirichletCondition[p[x,y]==1,x==0],DirichletCondition[p[x,y]==0,y==0&&x>0],DirichletCondition[p[x,y]==1,y==L]}; In this case a StreamPlot of component of velocities at the last time gives the following: Out[]= An animation of the ContourPlot of the module of velocity for each time step follows: Out[]= We notice that the flow is diagonal as we could expect. Other boundary conditions we can use consist in imposing no-slip conditions on all boundary but the left and right walls, the left wall will be an inflow boundary in which we impose parabolic horizontal velocity and null vertical velocity while the right wall will be an outflow boundary with pressure imposed to 0 boundaryconditionsopen={DirichletCondition[u[x,y]==0.,x>0&&x<L],DirichletCondition[v[x,y]==0.,x>0&&x<L],DirichletCondition[p[x,y]==0,x==L],DirichletCondition[v[x,y]==0.,x==0],DirichletCondition[u[x,y]==4*0.3*y*(L-y)/L^2,x==0]}; The StreamPlot of components of velocity gives: In[]:= StreamPlot[{xVel[15][x,y],yVel[15][x,y]},{x,y}∈xVel[15]["ElementMesh"],AspectRatio->Automatic,PlotRange->All,PlotLegendsAutomatic,ImageSize->Medium] Out[]= Notice that we obtain a result very similar to the first simulation, the flow is still horizontal, but here we have that velocity has higher values, especially at the left border where now is no more null.
Temperature-controlled laminar flows Let us maintain for now the same mesh we have used in the last section. Let us still assume that the density of the flow is constant and that gravity could be neglected. The Navier-Stokes equations coupled with the energy equation for two dimensional domains reduce to the following: ρ( ∂ t ∂ x ∂ y ∂ x ∂ 2 x ∂ 2 y ∂ t ∂ x ∂ y ∂ y ∂ 2 x ∂ 2 y C p ∂ t ∂ x ∂ y ∂ 2 x ∂ 2 y ∂ x ∂ y ρ=1 μ= Pr Ra C p k= 1 √(PrRa) 7.1 2* 5 10 op=ρ (1,0,0) u ∇ {x,y} ∇ {x,y} ∇ {x,y} (0,1,0) p (1,0,0) v ∇ {x,y} ∇ {x,y} ∇ {x,y} (0,0,1) p (0,1,0) u (0,0,1) v ρC p (1,0,0) T ∇ {x,y} ∇ {x,y} ρC p ∇ {x,y} C p Let us give boundary conditions for closed domains by imposing a temperature gradient between the walls, no-slip conditions for velocity on all boundary edges and a value for pressure in a point of the boundary: boundaryconditions2={DirichletCondition[{u[x,y]==0.,v[x,y]==0.},True],DirichletCondition[T[x,y]==1,x==0],DirichletCondition[T[x,y]==0,x==L],DirichletCondition[p[x,y]==0,x==0&&y==0]}; We will solve full Navier-Stokes equations through NDSolve. Let us use the Monitor function to control that times go on without problems during the computation. The structure of NDSolve is similar to the previous cases, we just enrich the Method Option. In[]:= Monitor[AbsoluteTiming[{xVel,yVel,pressure,temperature}=NDSolveValue[{op{0,0,0,0},bcs,ic},{u,v,p,T},{x,y}∈ω,{t,0,300},Method{"PDEDiscretization"{"MethodOfLines","SpatialDiscretization"{"FiniteElement","MeshOptions"{"MaxCellMeasure"maxCellMeasure},"InterpolationOrder"{u2,v2,p1,T2}}}},EvaluationMonitor(currentTime=Row[{"t = ",CForm[t]}])];],currentTime] The ContourPlots of pressure and temperature give the following: In[]:= The StreamPlot of component of velocities gives: Out[]= Notice that a circular flow on the more external pipes generates while the flow in the interior pipes is negligible. The values of Re x 1400 Re y 1700 We could make an animation of the StreamPlot of velocity coupled with the StreamPlot of the temperature: numberOfFrames=30;frames=TableShowbmesh["Wireframe"],ContourPlottemperature[t,x,y],{x,y}∈temperature["ElementMesh"],
Out[]=
The Sierpinski carpet Let us now construct a different type of meshes: the Sierpinski Carpet. They are created by iteration: ◼ at iteration 0 ◼ at iteration 1 the square is divided in a grid of nine squares and the central square is removed; ◼ at iteration 2 the same procedure of the previous iteration is applied to every one of the eight remaining square and so on. In Mathematica there is a function MengerMesh that takes in input the number of iterations and the dimension, that in two dimensions creates the Sierpinski carpet. In[]:= mengermesh[n_]=ToElementMesh[MengerMesh[n,2]];mengermesh1=mengermesh[1];mengermesh2=mengermesh[2];mengermesh3=mengermesh[3];{mengermesh1["Wireframe"],mengermesh2["Wireframe"],mengermesh3["Wireframe"]} Let us take the pressure-controlled flow’s operator with the same initial and boundary conditions of the first simulation. The StreamPlots for the first three orders of iteration at final time give: Notice that also in this case the flow is mainly horizontal. Moreover the flow is more and more difficult due to the exponential increasing of obstacles and this is testified by the fact that the bigger is the iteration the bigger is the Reynolds number. For example in the case of Re x If we take the temperature-controlled flow’s operator with the same initial and boundary conditions of the stationary simulation previously performed, the VectorPlots and StreamPlots for the first three orders of iteration give: Notice that also in this case we obtain a flow pattern similar to the previous, the flow is circular, but there is more. It seems that the increasing number of obstacles helps the flow to stabilize. Indeed the flow becomes more and more laminar, forexample Re x
The Menger sponge The Menger Sponge is the three dimensional version of the Sierpinski carpet. The way it is constructed is the extension to one more dimension of the one we have explained in the previous section. We have created the mesh with the same function as before but in dimension 3. generalmengermesh[n_,d_]=ToElementMesh[MengerMesh[n,d]];mengersponge1=generalmengermesh[1,3];mengersponge2=generalmengermesh[2,3];mengersponge3=generalmengermesh[3,3];{mengersponge1["Wireframe"],mengersponge2["Wireframe"],mengersponge3["Wireframe"]} Out[]= , , Let us still assume that the density of the flow is constant and that gravity could be neglected. The Navier-Stokes equations for three dimensional domains by components appear as: ρ( ∂ t ∂ x ∂ y ∂ z ∂ x ∂ x,x ∂ y,y ∂ z,z ∂ t ∂ x ∂ y ∂ z ∂ y ∂ x,x ∂ y,y ∂ z,z ∂ t ∂ x ∂ y ∂ z ∂ z ∂ x,x ∂ y,y ∂ z,z ∂ x ∂ y ∂ z ρ=1 μ= -4 10 ρ n+1 u n u Δt n u ∂ x n+1 u n v ∂ y n+1 u n w ∂ z n+1 u ∂ x n+1 p ∂ x,x n+1 u ∂ y,y n+1 u ∂ z,z n+1 u n+1 v n v Δt n u ∂ x n+1 v n v ∂ y n+1 v n w ∂ z n+1 v ∂ y n+1 p ∂ x,x n+1 v ∂ y,y n+1 v ∂ z,z n+1 v n+1 w n w Δt n u ∂ x n+1 w n v ∂ y n+1 w n w ∂ z n+1 w ∂ z n+1 p ∂ x,x n+1 w ∂ y,y n+1 w ∂ z,z n+1 w ∂ x n+1 u ∂ y n+1 v ∂ z n+1 w In[]:= dt=0.02;op5[i_]={ρ((u[x,y,z]-xVel[i-1][x,y,z])/dt+xVel[i-1][x,y,z] ∂ x ∂ y ∂ z ∂ x,x ∂ y,y ∂ z,z ∂ x ∂ x ∂ y ∂ z ∂ x,x ∂ y,y ∂ z,z ∂ y ∂ x ∂ y ∂ z ∂ x,x ∂ y,y ∂ z,z ∂ z ∂ x ∂ y ∂ z -4 10 Let us give the same boundary conditions of the first simulation extended to three dimensions, so no-slip boundary conditions for the three components of velocity and gradient of pressure between two walls (that now are surfaces and no more lines): In[]:= boundaryconditions5={DirichletCondition[{u[x,y,z]==0.,v[x,y,z]==0.,w[x,y,z]==0.},True],DirichletCondition[p[x,y,z]==0,x==1],DirichletCondition[p[x,y,z]==1,x==0]}; The solver also is an extension of the previous: xVel[0][x_,y_,z_]:=0;yVel[0][x_,y_,z_]:=0;zVel[0][x_,y_,z_]:=0;pressure[0][x_,y_,z_]:=0;Table[{xVel[i],yVel[i],zVel[i],pressure[i]}=NDSolveValue[{op5[i],boundaryconditions5},{u,v,w,p},{x,y,z}∈mengersponge,Method->{"FiniteElement","InterpolationOrder"->{u- |