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) + 2*PIxy(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