# Transition to turbulence

I found that a 3D model of viscous flow in a rectangular channel with jump section implemented on the basis of FEM demonstrates a transition to turbulence. We use an algorithm for unsteady flows described on

Effectively this algorithm works as LES. At one step, the average like values {UX[i], VY[i], WZ[i], P0[i]} are calculated here, and at another step it is used to calculate {u, v, w, p}. Here is the code and visualization of the flow:

H = 1/2; L = 2;
\[CapitalOmega] =
ImplicitRegion[
0 <= z <= H && 0 <= x <= L &&
0 <= y <=
H && ! (x >= 1 && y <= 0.1) && ! (x >= 1 &&
y >= 0.4) && ! (x >= 1 && z <= 0.1) && ! (x >= 1 &&
z >= 0.4), {x, y, z}];
RegionPlot3D[\[CapitalOmega]]
H=1/2;L=2;
\[CapitalOmega]=ImplicitRegion[0<=z<=H&&0<=x<=L&&0<=y<=H&&!(x>=1&&y<=0.1)&&!(x>=1&&y>=0.4)&&!(x>=1&&z<=0.1)&&!(x>=1&&z>=0.4),{x,y,z}];
RegionPlot3D[\[CapitalOmega]]

k=10;Um=45/100;\[Nu]=10^-4;t0=1/30;
Um=45/100;
U0[y_,z_]:=16*Um*y*z*(H-y)*(H-z)/H^4

UX[0][x_,y_,z_]:=0;
VY[0][x_,y_,z_]:=0;WZ[0][x_,y_,z_]:=0;
P0[0][x_,y_,z_]:=0;
Do[
{UX[i],VY[i],WZ[i],P0[i]}=NDSolveValue[{{-\[Nu]*Laplacian[u[x,y,z],{x,y,z}]+(p^(1,0,0))[x,y,z]+(u[x,y,z]-UX[i-1][x,y,z])/t0+UX[i-1][x,y,z]*D[u[x,y,z],x]+VY[i-1][x,y,z]*D[u[x,y,z],y]+WZ[i-1][x,y,z]*D[u[x,y,z],z],-\[Nu]*Laplacian[v[x,y,z],{x,y,z}]+(p^(0,1,0))[x,y,z]+(v[x,y,z]-VY[i-1][x,y,z])/t0+UX[i-1][x,y,z]*D[v[x,y,z],x]+VY[i-1][x,y,z]*D[v[x,y,z],y]+WZ[i-1][x,y,z]*D[v[x,y,z],z],-\[Nu]*Laplacian[w[x,y,z],{x,y,z}]+(p^(0,0,1))[x,y,z]+(w[x,y,z]-WZ[i-1][x,y,z])/t0+UX[i-1][x,y,z]*D[w[x,y,z],x]+VY[i-1][x,y,z]*D[w[x,y,z],y]+WZ[i-1][x,y,z]*D[w[x,y,z],z],(u^(1,0,0))[x,y,z]+(v^(0,1,0))[x,y,z]+(w^(0,0,1))[x,y,z]}=={0,0,0,0},{
DirichletCondition[{u[x,y,z]==U0[y,z],v[x,y,z]==0,w[x,y,z]==0},x==0],DirichletCondition[{u[x,y,z]==0,v[x,y,z]==0,w[x,y,z]==0},0<x<L],
DirichletCondition[p[x,y,z]==0,x==L]}},{u,v,w,p},{x,y,z}\[Element]\[CapitalOmega],Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,w->2,p->1},"MeshOptions"->{AccuracyGoal->5,PrecisionGoal->5,"MaxCellMeasure"->0.0001}}],{i,1,k}];


We see that this is similar to turbulence, although we did not use any model of turbulence. There is a question about the algorithm for calculating the interpolation function in 3D in the case of NDSolveValue[] with FEM. How is the interpolation function calculated in NDSolveValue[] with 3D FEM? Is it possible to make this algorithm faster?

Update 1. This algorithm can be used up to Reynolds number $$Re=4500$$ (it corresponds to $$\nu =5\times 10^{-5}$$ in the code above). Flow visualisation in the cross-section $$x=0.95$$ with k=22

frames = Table[
StreamDensityPlot[{VY[i][.95, y, z], WZ[i][.95, y, z]}, {y, 0,
H}, {z, 0, H}, ColorFunction -> "Rainbow", Frame -> False], {i,
1, 22}];


• I don't understand: This looks like you are solving the steady flow equations. Where does time come in? There is no "turbulence" in steady flows.
– Pirx
Mar 27 '20 at 22:22
• Do you see type (u[x,y,z]-UX[i-1][x,y,z])/t0 terms there? This is approximation of $\frac {\partial u}{\partial t}$ in a numerical scheme that I developed for transient flows. Mar 27 '20 at 23:33
• O.k., got you. We're looking at one long piece of "spaghetti code", so I missed that.
– Pirx
Mar 27 '20 at 23:39
• This does not work in another layout. Mar 27 '20 at 23:42