Error in dissipation rate profile

[Abishek Sarkar]

Dear all,

I am currently working in turbulent channel flows, and trying to generate the variation of the mean dissipation rate flow along the z-direction (direction normal to the walls of the channels). The mean TKE profile looks perfectly fine, but the dissipation rate profiles gives me vague oscillations and unphysical results. I guess, it is primarily related to the definition of the fluctuation velocity field and the corresponding formation of the strain rate tensor to define the dissipation rate. I am attaching the concerned excerpts of the code.

I will be very thankful to your kind suggestions.

Lx, Ly, Lz = 2*np.pi, np.pi, 1
Nx, Ny, Nz = 256, 160, 128

coords = d3.CartesianCoordinates('x', 'y', 'z')
dist = d3.Distributor(coords, dtype=dtype)

xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx), dealias=dealias)
ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(0, Ly), dealias=dealias)
zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz), dealias=dealias)

p = dist.Field(name='p', bases=(xbasis, ybasis, zbasis))
u = dist.VectorField(coords, name='u', bases=(xbasis, ybasis, zbasis))
tau_p  = dist.Field(name='tau_p')
tau_u1 = dist.VectorField(coords, name='tau_u1', bases=(xbasis, ybasis))
tau_u2 = dist.VectorField(coords, name='tau_u2', bases=(xbasis, ybasis))

x, y, z = dist.local_grids(xbasis, ybasis, zbasis)

lift_basis = zbasis.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
integ = lambda A: d3.Integrate(d3.Integrate(d3.Integrate(A, coords[0]), coords[1]), coords[2])
vol = Lx*Ly*Lz
avg = lambda A: integ(A)/vol  
dot = lambda A, B: d3.DotProduct(A, B)
ex, ey, ez = coords.unit_vector_fields(dist)

mean_uz_grad = d3.Average(d3.Average(uz_grad, coords['x']), coords['y'])

# Planar average of the velocity vector field
mean_u = d3.Average(d3.Average(u, coords['x']), coords['y'])

# Velocity fluctuation vector field
U_fluc = u - mean_u
mean_ux_xy = dot(mean_u, ex)

grad_u = d3.grad(u) + ez*lift(tau_u1)
grad_u_fluc = d3.grad(U_fluc)

strain_eps = 0.5 * (grad_u_fluc + d3.transpose(grad_u_fluc))
s2_eps = d3.trace(d3.dot(strain_eps, strain_eps))
epsilon = 2 * nu * s2_eps
epsilon_profile = d3.Average(d3.Average(epsilon, coords['x']), coords['y'])

profiles = solver.evaluator.add_file_handler('profiles', sim_dt=0.05, max_writes=50)

profiles.add_task(epsilon_profile, name='epsilon_profile')

I get rapid fluctuations, when I try to plot the mean (time-averaged) "epsilon_profile" vs z.

Any help would be greatly appreciated

Thanks 

Sincerely

Abishek