Run of turbulent channel flow produces a laminar flow (AIAA journal case))

Dear PYFR developers,

We are trying to reproduce the case of turbulent channel flow at Re_tau = 180 reported in the journal paper by Giangaspero, Witherden and Vincent (Synthetic Turbulence Generation for High-Order Scale-Resolving Simulations on Unstructured Grids).

The case is periodic one with the periodic boundary conditions imposed on the streamwise and spanwise boundaries of the computational domain. The grid mesh and the ini file are obtained from the paper. Following the paper, the case is run with 2nd order scheme for some time and then switched to the 4th order scheme. However, the resolved flow is always laminar no matter how long the case is run. The Q-criterion map for this case and ini file are attached. It will be greatly appreciated if you can have a look at the result and the ini file and give some suggestion on the how to set up this case in the ini file. The key to this case is the ini file, but we do not know how to correctly specify the setting in the ini file.

Thank you very much

The ini file:

[backend]
precision = double
rank-allocator = linear

[backend-cuda]
gimmik-max-nnz = 512
device-id = local-rank
;device-id = round-robin

[backend-openmp]
cc = gcc
cblas = /usr/lib/x86_64-linux-gnu/blas/libblas.so.3
;cblas = /libxsmm/lib/libxsmm.so
cblas-type = parallel

[constants]
gamma = 1.4
mu = 5.05967e-06
Pr = 0.72
pC = 1.0
Um = 1.0
Uc = 1.15833
Ul = 1.5
rhom = 0.014
tauw = 5.92464e-05
delta = 1
nu = 0.000361405
utau = 0.06397952655150352 ; =Um/15.63 (Kim Moin Moser 1987)
ac-zeta = 5.0

[solver]
system = ac-navier-stokes
order = 4
anti-alias = none

[solver-time-integrator]
formulation = dual
scheme = bdf2
pseudo-scheme = rk45
controller = none
pseudo-controller = local-pi
tstart = 0.0
tend = 1221.0
dt = 0.02 ;working for all
;pseudo-dt = 0.003 ;working p2
pseudo-dt = 0.0007 ;working p4 ; all hexes
pseudo-niters-min = 1
pseudo-niters-max = 25 ;15
pseudo-resid-norm = l2
pseudo-resid-tol = 2e-5 ;1e-5 ;5e-5
pseudo-resid-tol-p = 1.0
atol = 1e-6
pseudo-dt-max-mult = 2.5
safety-fact =0.9 
min-fact = 0.98
max-fact= 1.01

[solver-dual-time-integrator-multip]
pseudo-dt-fact = 1.9 ; working p2 and p4 all hexes
;cycle = [(2, 1), (1, 1), (0, 4), (1, 1), (2, 2)] ;working p2
cycle = [(4, 1), (3, 1), (2, 1), (1, 1), (0, 10), (1, 2), (2, 2), (3, 2), (4, 4)] ;working p4

[solver-interfaces]
riemann-solver = rusanov
ldg-beta = 0.5
ldg-tau = 0.1

[solver-interfaces-line]
flux-pts = gauss-legendre
quad-deg = 11
quad-pts = gauss-legendre

[solver-interfaces-tri]
flux-pts = williams-shunn
quad-deg = 11
quad-pts = williams-shunn

[solver-interfaces-quad]
flux-pts = gauss-legendre
quad-deg = 11
quad-pts = gauss-legendre

[solver-elements-tri]
soln-pts = williams-shunn
quad-deg = 11
quad-pts = williams-shunn

[solver-elements-quad]
soln-pts = gauss-legendre
quad-deg = 11
quad-pts = gauss-legendre

[solver-elements-hex]
soln-pts = gauss-legendre
quad-deg = 11
quad-pts = gauss-legendre

[solver-elements-tet]
soln-pts = shunn-ham
quad-deg = 7
quad-pts = shunn-ham

[solver-elements-pri]
soln-pts = williams-shunn~gauss-legendre
quad-deg = 11
quad-pts = williams-shunn~gauss-legendre

[solver-elements-pyr]
soln-pts = gauss-legendre
quad-deg = 11
quad-pts = witherden-vincent

[solver-interfaces-line-mg-p3]
flux-pts = gauss-legendre
quad-deg = 9
quad-pts = gauss-legendre

[solver-interfaces-tri-mg-p3]
flux-pts = williams-shunn
quad-deg = 9
quad-pts = williams-shunn

[solver-interfaces-quad-mg-p3]
flux-pts = gauss-legendre
quad-deg = 9
quad-pts = gauss-legendre

[solver-elements-tri-mg-p3]
soln-pts = williams-shunn
quad-deg = 9
quad-pts = williams-shunn

[solver-elements-quad-mg-p3]
soln-pts = gauss-legendre
quad-deg = 9
quad-pts = gauss-legendre

[solver-elements-hex-mg-p3]
soln-pts = gauss-legendre
quad-deg = 9
quad-pts = gauss-legendre

[solver-elements-tet-mg-p3]
soln-pts = shunn-ham
quad-deg = 7
quad-pts = shunn-ham

[solver-elements-pri-mg-p3]
soln-pts = williams-shunn~gauss-legendre
quad-deg = 9
quad-pts = williams-shunn~gauss-legendre

[solver-elements-pyr-mg-p3]
soln-pts = gauss-legendre
quad-deg = 9
quad-pts = witherden-vincent

[solver-interfaces-line-mg-p2]
flux-pts = gauss-legendre
quad-deg = 7
quad-pts = gauss-legendre

[solver-interfaces-tri-mg-p2]
flux-pts = williams-shunn
quad-deg = 7
quad-pts = williams-shunn

[solver-interfaces-quad-mg-p2]
flux-pts = gauss-legendre
quad-deg = 7
quad-pts = gauss-legendre

[solver-elements-tri-mg-p2]
soln-pts = williams-shunn
quad-deg = 7
quad-pts = williams-shunn

[solver-elements-quad-mg-p2]
soln-pts = gauss-legendre
quad-deg = 7
quad-pts = gauss-legendre

[solver-elements-hex-mg-p2]
soln-pts = gauss-legendre
quad-deg = 7
quad-pts = gauss-legendre

[solver-elements-tet-mg-p2]
soln-pts = shunn-ham
quad-deg = 7
quad-pts = shunn-ham

[solver-elements-pri-mg-p2]
soln-pts = williams-shunn~gauss-legendre
quad-deg = 7
quad-pts = williams-shunn~gauss-legendre

[solver-elements-pyr-mg-p2]
soln-pts = gauss-legendre
quad-deg = 7
quad-pts = witherden-vincent

[solver-interfaces-line-mg-p1]
flux-pts = gauss-legendre
quad-deg = 5
quad-pts = gauss-legendre

[solver-interfaces-tri-mg-p1]
flux-pts = williams-shunn
quad-deg = 5
quad-pts = williams-shunn

[solver-interfaces-quad-mg-p1]
flux-pts = gauss-legendre
quad-deg = 5
quad-pts = gauss-legendre

[solver-elements-tri-mg-p1]
soln-pts = williams-shunn
quad-deg = 5
quad-pts = williams-shunn

[solver-elements-quad-mg-p1]
soln-pts = gauss-legendre
quad-deg = 5
quad-pts = gauss-legendre

[solver-elements-hex-mg-p1]
soln-pts = gauss-legendre
quad-deg = 5
quad-pts = gauss-legendre

[solver-elements-tet-mg-p1]
soln-pts = shunn-ham
quad-deg = 5
quad-pts = shunn-ham

[solver-elements-pri-mg-p1]
soln-pts = williams-shunn~gauss-legendre
quad-deg = 5
quad-pts = williams-shunn~gauss-legendre

[solver-elements-pyr-mg-p1]
soln-pts = gauss-legendre
quad-deg = 5
quad-pts = witherden-vincent

[solver-interfaces-line-mg-p0]
flux-pts = gauss-legendre
quad-deg = 3
quad-pts = gauss-legendre

[solver-interfaces-tri-mg-p0]
flux-pts = williams-shunn
quad-deg = 3
quad-pts = williams-shunn

[solver-interfaces-quad-mg-p0]
flux-pts = gauss-legendre
quad-deg = 3
quad-pts = gauss-legendre

[solver-elements-tri-mg-p0]
soln-pts = williams-shunn
quad-deg = 3
quad-pts = williams-shunn

[solver-elements-quad-mg-p0]
soln-pts = gauss-legendre
quad-deg = 3
quad-pts = gauss-legendre

[solver-elements-hex-mg-p0]
soln-pts = gauss-legendre
quad-deg = 3
quad-pts = gauss-legendre

[solver-elements-tet-mg-p0]
soln-pts = shunn-ham
quad-deg = 3
quad-pts = shunn-ham

[solver-elements-pri-mg-p0]
soln-pts = williams-shunn~gauss-legendre
quad-deg = 3
quad-pts = williams-shunn~gauss-legendre

[solver-elements-pyr-mg-p0]
soln-pts = gauss-legendre
quad-deg = 3
quad-pts = witherden-vincent

[soln-plugin-nancheck]
nsteps = 100

[soln-plugin-pseudostats]
flushsteps = 50
file = pseudo-residual.csv
header = true

[soln-plugin-writer]
dt-out = 100.0
basedir = .
basename = Channel-{t:010.4f}

[soln-bcs-wall]
type = no-slp-wall

[soln-ics]
; parabolic profile
;u = Ul-Ul*y*y/delta/delta
; DNS turb profile
u = 1.36346011e+05*pow(1-abs(y), 14)-1.00157155e+06*pow(1-abs(y), 13)+3.30835371e+06*pow(1-abs(y), 12)-6.49079211e+06*pow(1-abs(y), 11)+8.41556333e+06*pow(1-abs(y), 10)-7.58969284e+06*pow(1-abs(y), 9)+4.87934144e+06*pow(1-abs(y), 8)-2.25289875e+06*pow(1-abs(y), 7)+7.41903706e+05*pow(1-abs(y), 6)-1.70090684e+05*pow(1-abs(y), 5)+2.56967948e+04*pow(1-abs(y), 4)-2.21227945e+03*pow(1-abs(y), 3)+4.33607541e+01*pow(1-abs(y), 2)+1.10361842e+01*pow(1-abs(y), 1)+3.28951994e-04*pow(1-abs(y), 0)
v = 0.0
w = 0.0
p = pC - tauw*x/delta/rhom

[solver-source-terms]
u = tauw/delta/rhom

[soln-plugin-tavg]
nsteps = 1
dt-out = 100.0
basedir = .
basename = Channel-tavg-{t:010.4f}
avg-p = p
avg-u = u
avg-v = v
avg-w = w
avg-pp = p*p
avg-uu = u*u
avg-uv = u*v
avg-uw = u*w
avg-vv = v*v
avg-vw = v*w
avg-ww = w*w

;[soln-plugin-sampler-1]
;nsteps = 1
;samp-pts = [(3.00, -0.99, 0.00), (3.00, -0.997, 0.00), (3.00, -0.995, 0.00), (3.00, -0.989, 0.00), (3.00, -0.984, 0.00), (3.00, -0.978, 0.00), (3.00, -0.973, 0.00), (3.00, -0.967, 0.00), (3.00, -0.957, 0.00), (3.00, -0.951, 0.00), (3.00, -0.946, 0.00), (3.00, -0.94, 0.00), (3.00, -0.935, 0.00), (3.00, -0.93, 0.00), (3.00, -0.924, 0.00), (3.00, -0.919, 0.00), (3.00, -0.892, 0.00), (3.00, -0.864, 0.00), (3.00, -0.837, 0.00), (3.00, -0.81, 0.00), (3.00, -0.783, 0.00), (3.00, -0.729, 0.00), (3.00, -0.675, 0.00), (3.00, -0.621, 0.00), (3.00, -0.512, 0.00), (3.00, -0.458, 0.00), (3.00, -0.35, 0.00), (3.00, -0.241, 0.00), (3.00, -0.187, 0.00), (3.00, -0.079, 0.00), (3.00, 0.00, 0.00)]
;format = primitive
;file = ./probes/vertical.csv
;header = true

[soln-plugin-sampler-0]
;z = np.linspace(0, np.pi/2., num=24, endpoint=False)
nsteps = 1
samp-pts = [[3.0,0.0,0.0],[3.0,0.0,0.06544984694978735],[3.0,0.0,0.1308996938995747],[3.0,0.0,0.19634954084936207],[3.0,0.0,0.2617993877991494],[3.0,0.0,0.32724923474893675],[3.0,0.0,0.39269908169872414],[3.0,0.0,0.4581489286485115],[3.0,0.0,0.5235987755982988],[3.0,0.0,0.5890486225480862],[3.0,0.0,0.6544984694978735],[3.0,0.0,0.7199483164476609],[3.0,0.0,0.7853981633974483],[3.0,0.0,0.8508480103472356],[3.0,0.0,0.916297857297023],[3.0,0.0,0.9817477042468102],[3.0,0.0,1.0471975511965976],[3.0,0.0,1.112647398146385],[3.0,0.0,1.1780972450961724],[3.0,0.0,1.2435470920459597],[3.0,0.0,1.308996938995747],[3.0,0.0,1.3744467859455345],[3.0,0.0,1.4398966328953218],[3.0,0.0,1.505346479845109]]
format = primitive
;file = ./probes/horizontal_x3y0.csv
file = point-data.csv
header = true

;[soln-plugin-sampler-1]
;nsteps = 1
;samp-pts = [(20.0000,-0.9700,-1.4960),(20.0000,-0.9700,-0.7480),(20.0000,-0.9700,0.0000),(20.0000,-0.9700,0.7480),(20.0000,-0.9700,1.4960)]
;format = primitive
;file = ./probes/horizontal_x20.0000_y-0.9700.csv
;header = true

;[soln-plugin-sampler-2]
;nsteps = 1
;samp-pts = [(20.0000,-0.1710,-1.4960),(20.0000,-0.1710,-0.7480),(20.0000,-0.1710,0.0000),(20.0000,-0.1710,0.7480),(20.0000,-0.1710,1.4960)]
;format = primitive
;file = ./probes/horizontal_x20.0000_y-0.1710.csv
;header = true

Turbulent channel flows can be tricky because of this. They often can take quite a long time to transition if you don’t give them a perturbation to get them going.

Looking at your initial conditions you are just initialising the flow with a boundary layer profile, you should add on some small sine wave perturbations to the velocities. This should get the boundary layer to transition sooner.

Hi. I want to add something.
Another way is to reduce viscosity at the beginning. For example, start with the one tenth of the viscosity. you can then switch to the targeted viscosity as soon as you observe some instability/turbulence. Re=180 is indeed really low. together with the numerical viscosity, the viscosity hampers initial small instabilities.