Transition to turbulence

Question

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

https://community.wolfram.com/groups/-/m/t/1433064 Solver for unsteady flow with the use of the Navier-Stokes and Mathematica FEM

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}];