# Smagorinsky subgrid model

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

if your flow is laminar, then it’s to be expected that Q is negligible. you must increase the reynolds number to get turbulence; the purpose of the smagorinsky model is to stabilize underresolved turbulent simulation, not to produce turbulence.

It’s clear Adam; when I write “in order to simulate a turbulent flow” I mean that my Reynolds number will be set to 1000, or 10000. And Q is yet low when Reynolds is 10^8!

Hallo gianluca,

at the moment I’m trying the exact same thing like you two years ago. And this with the same problem!
Could you solve the implementation of the Smagorinsky subgrid model in Matlab? If yes, could you please tell or show me how you did it?