Hello,
I developped LES LBM code, and I have some problems, for uniforme mesh grid dx=dy=1 with smagroinsky turbulence model, I can’t simulate a high Reynolds number, and the Mach number is high then 0.5. when I use a space step more smaller, the code has a problem, the pressure increase every step, and the flow still in the initial conditions.
Thanks for any advice,
Nouri
Perhaps could you give a preview of the LES code part? It sounds like your stress tensor part of the program that determines the eddy viscosity isn’t working correctly?
hi Ivan ,I use fortran.
c=dx/dt
write(*,*)'vitesse de propagation =',c
do k=1,4
cx(k)=cos((k-1)*pi/2)*c
cy(k)=sin((k-1)*pi/2)*c
ENDDO
do k=5,8
cx(k)=cos((k-5)*pi/2+pi/4)*c*sqrt(2.)
cy(k)=sin((k-5)*pi/2+pi/4)*c*sqrt(2.)
ENDDO
cx(0)=0
cy(0)=0
csound=c/sqrt(3.)
write(*,*)'Vitesse de son =', csound
molecular relaxation time
Re=700
write(*,*)'Nombre de Re' , Re
rhoo=5
obstr=20
write(*,*)'Diamètre Lattice',obstr
viscosl=(obstr*u0)/Re
write(*,*)'Viscosité Lattice', viscosl
alpha=viscosl/pr
tau=3*(dt)*viscosl/(dx**2)+0.5
print *,' tau =',tau
omega=1/(tau)
Equilibre function
do i=0,n
do j=0,m
t1=(u(i,j)*u(i,j)+v(i,j)*v(i,j))/(csound**2)
do k=0,8
t2=(u(i,j)*cx(k)+v(i,j)*cy(k))/(csound**2)
feq(k,i,j)=rho(i,j)*w(k)*(1.0+t2+1./2.*(t2*t2)-1./2.*t1)
end do
end do
end do
share stress tensor
do i=0,n
do j=0,m
do k=0,8
PXX(i,j)=PXX(i,j)+cx(k)*cX(k)*(f(k,i,j)-feq(k,i,j))
PXY(i,j)=Pxy(i,j)+cx(k)*cy(k)*(f(k,i,j)-feq(k,i,j))
PYY(i,j)=Pyy(i,j)+cy(k)*cy(k)*(f(k,i,j)-feq(k,i,j))
enddo
Q(i,j)=Pxx(i,j)*Pxx(i,j) + 2*Pxy(i,j)*Pxy(i,j) + Pyy(i,j)*Pyy(i,j)
Q(i,j)=sqrt(q(i,j))
enddo
enddo
collision step
omegatot(i,j)=1d0/(0.5*(sqrt((tau*tau) + (C_Smagorinsky**2)*9*((8*Q(i,j))**(0.5))/rho(i,j) )+(tau)))
f(k,i,j)=(omegatot(i,j))*feq(k,i,j)+(1.-omegatot(i,j))*f(k,i,j)
Perhaps try the smagorinsky part of the code without dividng by rho. Atleast that is how my code works.
I didn’t go and check that the velocities and speed of sound are correct but generally I think it should work. Otherwise check you’ve got the correct boundary conditions.
[code=“fortran”]
! lattice velocities
c=dx/dt
write(,)‘vitesse de propagation =’,c
do k = 1, 4
cx(k)=cos((k-1)*pi/2)*c
cy(k)=sin((k-1)*pi/2)*c
end do
do k = 5, 8
cx(k)=cos((k-5)pi/2+pi/4)csqrt(2.)
cy(k)=sin((k-5)pi/2+pi/4)csqrt(2.)
end do
cx(0)=0
cy(0)=0
csound=c/sqrt(3.0e0)
write(,)‘Vitesse de son =’, csound
! relaxation time
Re=700
write(,)‘Nombre de Re’ , Re
rhoo=5
obstr=20
write(,)‘Diamètre Lattice’,obstr
viscosl=(obstru0)/Re
write(,)‘Viscosité Lattice’, viscosl
alpha=viscosl/pr
tau=3(dt)*viscosl/(dx**2)+0.5
print *,’ tau =’,tau
omega=1/(tau)
! Equilibre function
do i = 0, n
do j = 0, m
t1 = (u(i,j)u(i,j) + v(i,j)v(i,j))/(csound2)
do k=0,8
t2=(u(i,j)*cx(k) + v(i,j)*cy(k))/(csound2)
feq(k,i,j) = w(k)rho(i,j)(1.0 + t2 +1./2.(t2t2) - 1./2.*t1)
end do
end do
end do
! compute stress tensor
do i = 0, n
do j = 0, m
do k = 1, 8
PXX(i,j) = PXX(i,j) + cx(k)cx(k)(f(k,i,j) - feq(k,i,j))
PXY(i,j) = PXY(i,j) + cx(k)cy(k)(f(k,i,j) - feq(k,i,j))
PYY(i,j) = PYY(i,j) + cy(k)cy(k)(f(k,i,j) - feq(k,i,j))
enddo
Q(i,j) = PXX(i,j)PXX(i,j) + 2.0e0PXY(i,j)*PXY(i,j) + PYY(i,j)*PYY(i,j)
Q(i,j) = sqrt(Q(i,j))
enddo
enddo
! compute omega
tau_t = 0.5e0*(sqrt(tau2+18.0e0*(C_Smagorinsky2)*sqrtVar)-tau)
omegaTot(i,j) = 1.0e0/(tau_t + tau)
! collision step
f(k,i,j)=(omegaTot(i,j))*feq(k,i,j) + (1.0e0 - omegaTot(i,j))*f(k,i,j)
This is a subroutine from my Fortran program that produces the results in the image beneath:
[code="fortran"]
!=======================================================================
! Stream & Collide
!=======================================================================
subroutine StreamCollide(f,rho)
use simParam, only: ny, nx, uMax, tau0, ni, old, new,num
use d2q9Lattice, only: wf,ef
implicit none
real(kind=num), intent(in out):: f (0:nx+1,0:ny+1,0:8,2)
real(kind=num), intent(out):: rho(nx)
real(kind=num) :: ux(nx), uy(nx), uxy(nx),uSqr(nx), ne(nx,0:8)
real(kind=num):: omega(nx), minOmega(nx)
real(kind=num):: C = 0.16, stressTen(0:2), sqrtVar, S, tau_t
integer:: x, y, i
do y = 1, ny
! macroscopic variables
do x = 1, nx
rho(x) = sum(f(x,y,:,old))
ux(x) = (f(x,y,1,old)+f(x,y,5,old)+f(x,y,8,old)-f(x,y,3,old)-f(x,y,6,old)-f(x,y,7,old))/rho(x)
uy(x) = (f(x,y,2,old)+f(x,y,5,old)+f(x,y,6,old)-f(x,y,4,old)-f(x,y,7,old)-f(x,y,8,old))/rho(x)
uSqr(x) = ux(x)*ux(x)+uy(x)*uy(x)
end do
! non-equilibrium parts
do x = 1, nx
ne(x,0) = wf(0)*rho(x)*(1.0 -1.5*uSqr(x))
ne(x,1) = wf(1)*rho(x)*(1.0+3.0*ux(x) +4.5*ux(x)*ux(x) -1.5*uSqr(x))
ne(x,2) = wf(2)*rho(x)*(1.0+3.0*uy(x) +4.5*uy(x)*uy(x) -1.5*uSqr(x))
ne(x,3) = wf(3)*rho(x)*(1.0-3.0*ux(x) +4.5*ux(x)*ux(x) -1.5*uSqr(x))
ne(x,4) = wf(4)*rho(x)*(1.0-3.0*uy(x) +4.5*uy(x)*uy(x) -1.5*uSqr(x))
uxy(x) = ux(x) + uy(x)
ne(x,5) = wf(5)*rho(x)*(1.0+3.0*uxy(x)+4.5*uxy(x)*uxy(x)-1.5*uSqr(x))
uxy(x) = -ux(x) + uy(x)
ne(x,6) = wf(6)*rho(x)*(1.0+3.0*uxy(x)+4.5*uxy(x)*uxy(x)-1.5*uSqr(x))
uxy(x) = -ux(x) - uy(x)
ne(x,7) = wf(7)*rho(x)*(1.0+3.0*uxy(x)+4.5*uxy(x)*uxy(x)-1.5*uSqr(x))
uxy(x) = ux(x) - uy(x)
ne(x,8) = wf(8)*rho(x)*(1.0+3.0*uxy(x)+4.5*uxy(x)*uxy(x)-1.5*uSqr(x))
end do
! compute local strain rate
do x = 1, nx
stressTen = 0.0
do i = 1, 8
stressTen(0) = stressTen(0) + ef(i,0)*ef(i,0)*(f(x,y,i,old)-ne(x,i))
stressTen(1) = stressTen(1) + ef(i,0)*ef(i,1)*(f(x,y,i,old)-ne(x,i))
stressTen(2) = stressTen(2) + ef(i,1)*ef(i,1)*(f(x,y,i,old)-ne(x,i))
end do
! compute magnitude of strain rate
sqrtVar = sqrt(stressTen(0)**2+2.0*stressTen(1)**2+stressTen(2)**2)
tau_t = 0.5*(sqrt(tau0**2+18.0*C**2*sqrtVar)-tau0)
! compute local collision frequency
omega(x) = 1.0/(tau0+tau_t)
minOmega(x) = 1.0 - omega(x)
end do
! stream collide
do x = 1, nx
f(x,y,0,new) = f(x,y,0,old)*minOmega(x) + omega(x)*ne(x,0)
end do
do x = 1, nx
f(x+1,y,1,new) = f(x,y,1,old)*minOmega(x) + omega(x)*ne(x,1)
end do
do x = 1, nx
f(x,y+1,2,new) = f(x,y,2,old)*minOmega(x) + omega(x)*ne(x,2)
end do
do x = 1, nx
f(x-1,y,3,new) = f(x,y,3,old)*minOmega(x) + omega(x)*ne(x,3)
end do
do x = 1, nx
f(x,y-1,4,new) = f(x,y,4,old)*minOmega(x) + omega(x)*ne(x,4)
end do
do x = 1, nx
f(x+1,y+1,5,new) = f(x,y,5,old)*minOmega(x) + omega(x)*ne(x,5)
end do
do x = 1, nx
f(x-1,y+1,6,new) = f(x,y,6,old)*minOmega(x) + omega(x)*ne(x,6)
end do
do x = 1, nx
f(x-1,y-1,7,new) = f(x,y,7,old)*minOmega(x) + omega(x)*ne(x,7)
end do
do x = 1, nx
f(x+1,y-1,8,new) = f(x,y,8,old)*minOmega(x) + omega(x)*ne(x,8)
end do
end do
end subroutine StreamCollide
Thank you, Ivan, can you send me please your code to verify other simulation parameter, nouri.mdafer@gmail.com