Hi dear all,
I have a problem with the implementation of the Smagorinsky turbulence model.
I wrote a Matlab code able to perform the simulation of a laminar flow in a cylinder using a d2q9 scheme. I imposed a constant velocity profile at the inlet (u=0.01 in lattice Boltzmann units) and a costant density at the outlet, according to Zou and He boundary conditions. I use bounce back conditions for the walls. At low Reynolds I obtain the Poiseulle velocity profile after about 8000 iterations, so I suppose the code to be correct.
In order to simulate a turbulent flow, I referred to the Smagorinsky turbulence model; in this case I calculate the total local relaxation parameter as follows:
omegatot(x,y)=(0.5*(((1/omega(x,y))^2 + (CSD^2)9(8*Q(x,y))^(0.5)/rho(x,y) )^0.5+(1/omega(x,y))))^(-1);
where CSD could approximately be 0.1 as founded into literature and Q is related to the magnitude of the strain-rate tensor. In order to determine the value of Q I wrote the code as follows:
[b]for x=1:nx
for y=1:ny
for i=1:9
Feq(i,x,y)=rho(x,y)*w(i)*(1 + 3*(cx(i)*ux(x,y)+cy(i)*uy(x,y))/c^2 + (9/(2*c^4))*(cx(i)*ux(x,y)+cy(i)*uy(x,y))^2 - (3/(2*c^2))*(ux(x,y)^2+uy(x,y)^2));
Fneq(i,x,y)=F(i,x,y)-Feq(i,x,y);
end
end
end
%Calcolo del tensore degli sforzi
for x=1:nx
for y=1:ny
PIxx(x,y)=cx(1)*cx(1)*Fneq(1,x,y)+cx(2)*cx(2)*Fneq(2,x,y)+cx(3)*cx(3)*Fneq(3,x,y)+cx(4)*cx(4)*Fneq(4,x,y)+cx(5)*cx(5)*Fneq(5,x,y)+cx(6)*cx(6)*Fneq(6,x,y)+cx(7)*cx(7)*Fneq(7,x,y)+cx(8)*cx(8)*Fneq(8,x,y)+cx(9)*cx(9)*Fneq(9,x,y);
PIyy(x,y)=cy(1)*cy(1)*Fneq(1,x,y)+cy(2)*cy(2)*Fneq(2,x,y)+cy(3)*cy(3)*Fneq(3,x,y)+cy(4)*cy(4)*Fneq(4,x,y)+cy(5)*cy(5)*Fneq(5,x,y)+cy(6)*cy(6)*Fneq(6,x,y)+cy(7)*cy(7)*Fneq(7,x,y)+cy(8)*cy(8)*Fneq(8,x,y)+cy(9)*cy(9)*Fneq(9,x,y);
PIxy(x,y)=cx(1)*cy(1)*Fneq(1,x,y)+cx(2)*cy(2)*Fneq(2,x,y)+cx(3)*cy(3)*Fneq(3,x,y)+cx(4)*cy(4)*Fneq(4,x,y)+cx(5)*cy(5)*Fneq(5,x,y)+cx(6)*cy(6)*Fneq(6,x,y)+cx(7)*cy(7)*Fneq(7,x,y)+cx(8)*cy(8)*Fneq(8,x,y)+cx(9)*cy(9)*Fneq(9,x,y);
Q(x,y)=PIxx(x,y)PIxx(x,y) + 2PIxy(x,y)*PIxy(x,y) + PIyy(x,y)*PIyy(x,y);
CSD = 0.1;
omegatot(x,y)=(0.5*(((1/omega(x,y))^2 + (CSD^2)9(8*Q(x,y))^(0.5)/rho(x,y) )^0.5+(1/omega(x,y))))^(-1);
nut(x,y)=((1/omegatot(x,y))-0.5)/3-nu;
end
end[/b]
where:
cx = [0 1 0 -1 0 1 -1 -1 1];
cy = [0 0 1 0 -1 1 1 -1 -1];
The problem is that the magnitude of Q is so low (10^-7) that the term CSD^2)9(8*Q(x,y))^(0.5)/rho(x,y) is so small to become uninfluential on the omegatot definition, although its profile looks to be correct (lower in the centre than near the wall).
Is my definition of Q wrong?
Thank you very much, I’m quite inexpert as you can see and I don’t know how to solve this problem.
Gianluca