BGK / MRT turbulence in Matlab
用的Smagorinsky model,稍微简化了一下, 不过现有几个model似乎差别不大clc
clear all
close all
x=100
y=200
tau=0.5*ones(x,y);
step=100
u0=0.05;
%% D2Q9 setup
ex(1,1,1:9)=;
ey(1,1,1:9)=;
c0=;c1=;c2=;c3=[-1,0];c4=;c5=;c6=[-1,1];c7=[-1,-1];c8=;
w(1,1,1:9)=;
%% setup solid boundary
=meshgrid(1:y,1:x);
Solid=false(x,y);
Solid=Solid | (j-x/3).^2+(i-y/5).^2<(x/20)^2;
Solid=Solid | (j-x/2).^2+(i-y/2).^2<(x/20)^2;
Solid=Solid | (j-x/1.5).^2+(i-y/1.25).^2<(x/20)^2;
Solid(1,:)=true;
Solid(x,:)=true;
Fluid=~Solid;
Inlet=Fluid & i==1;
Outlet=i==y;
FluidAll=repmat(Fluid,);
%% used for streaming indexing
Ind=reshape(1:x*y*9,x,y,9);
Ind(:,:,1+1)=circshift(Ind(:,:,1+1),c1);
Ind(:,:,2+1)=circshift(Ind(:,:,2+1),c2);
Ind(:,:,3+1)=circshift(Ind(:,:,3+1),c3);
Ind(:,:,4+1)=circshift(Ind(:,:,4+1),c4);
Ind(:,:,5+1)=circshift(Ind(:,:,5+1),c5);
Ind(:,:,6+1)=circshift(Ind(:,:,6+1),c6);
Ind(:,:,7+1)=circshift(Ind(:,:,7+1),c7);
Ind(:,:,8+1)=circshift(Ind(:,:,8+1),c8);
%% used for bounceback indexing
Blank=false(x,y);
SolidBB1=cat(3,Blank,Solid,Solid,Blank,Blank,Solid,Solid,Blank,Blank);
SolidBB2=cat(3,Blank,Blank,Blank,Solid,Solid,Blank,Blank,Solid,Solid);
temp=Ind(SolidBB1);Ind(SolidBB1)=Ind(SolidBB2);Ind(SolidBB2)=temp;
%% viscosity
Re=Inf;
Cs=5e-3
viscosity=u0*y/Re
tau0=3*viscosity+.5;
%% LES setup
exx(1,1,1:9)=;
eyy(1,1,1:9)=;
exy(1,1,1:9)=;
%% MRT setup
M=[1 -4 4 0 0 0 0 0 0
1 -1 -2 1 -2 0 0 1 0
1 -1 -2 0 0 1 -2 -1 0
1 -1 -2 -1 2 0 0 1 0
1 -1 -2 0 0 -1 2 -1 0
1 2 1 1 1 1 1 0 1
1 2 1 -1 -1 1 1 0 -1
1 2 1 -1 -1 -1 -1 0 1
1 2 1 1 1 -1 -1 0 -1];
Mt=M';
M0=diag(1./sum(M.*M))*Mt;
S0=;
% S0=;
S0=ones(x*y,1)*S0;
%% UI control
restart = uicontrol('position',,'style','toggle','string','restart');
quit = uicontrol('position',,'style','toggle','string','close');
pausebutton = uicontrol('position',,'style','toggle','string','pause');
bgkmrt = uicontrol('position',,'style','toggle','string','switch to bgk');
%%
while get(quit,'value') == 0
set(restart,'value',0);
set(pausebutton,'value',0);
%% initial condition
nTimestep=0;
ux=zeros(x,y);
uy=u0*ones(x,y);
forcex=zeros(x,y);
forcey(1:x,1:y)=2e-5;
rho=ones(x,y);
usq=ux.*ux+uy.*uy;
c=(bsxfun(@times,ex,ux)+bsxfun(@times,ey,uy));
density=bsxfun(@times,rho,w) .* bsxfun(@plus,1+3*c+4.5*c.*c,-1.5*usq);
%% loop
while get(restart,'value')==0 && get(quit,'value')==0
while get(pausebutton,'value')==1 && get(restart,'value')==0 && get(quit,'value')==0
set(pausebutton,'string','continue');pause(.01);end;set(pausebutton,'string','pause');%% restart and quit are volatile, so extra test is required
%% calc macro
rho=sum(density,3);
ux=sum(bsxfun(@times,density,ex),3)./rho;
uy=sum(bsxfun(@times,density,ey),3)./rho;
%% forcing
ux=ux+forcex./rho;
uy=uy+forcey./rho;
%% collision
usq=ux.*ux+uy.*uy;
c=(bsxfun(@times,ex,ux)+bsxfun(@times,ey,uy));
density_eq=bsxfun(@times,rho,w) .* bsxfun(@plus,1+3*c+4.5*c.*c,-1.5*usq);
density_delta=density_eq-density;
moment=reshape(density,)*M; ux_=ux(:);uy_=uy(:);usq_=usq(:);
moment_eq=bsxfun(@times,rho(:),);
%% eddy viscosity with Smagorinsky model
qxx=sum(bsxfun(@times,density_delta,exx),3);
qyy=sum(bsxfun(@times,density_delta,eyy),3);
qxy=sum(bsxfun(@times,density_delta,exy),3);
q=sqrt(qxx.*qxx+qyy.*qyy+2*qxy.*qxy);
tau=tau0+.5*(sqrt(tau0.*tau0+36*Cs*q)-tau0);
%% collision cont'd
if(get(bgkmrt,'value')==0)
set(bgkmrt,'string','switch to bgk')
S=;
else
set(bgkmrt,'string','switch to mrt')
S=1./tau(:)*ones(1,9);
end
moment=moment+bsxfun(@times,moment_eq-moment,S);
density_temp=reshape(moment*M0,);
% density_temp=density+bsxfun(@rdivide,density_delta,tau);
density(FluidAll)=density_temp(FluidAll);
%% stream & bounce back
density=density(Ind);
%%
nTimestep=nTimestep+1;
if(~mod(nTimestep,step))
uxplot=uy.*Fluid;
uyplot=ux.*Fluid;
umag=sqrt(uxplot.*uxplot+uyplot.*uyplot);
vmax=max(umag(:));
vel=(256/vmax)*(umag);
vor=curl(uxplot,uyplot)*20000+128;
subplot(2,1,1)
pcolor(vel)
xlabel('Velocity')
shading interp
axis equal
axis()
title(sprintf('%d, tau=%f, u_m_a_x=%f, Vorticity_m_a_x=%f, Mass=%f',nTimestep,tau(x/2,y/5),vmax,max(vor(:)),sum(rho(:))));
subplot(2,1,2)
pcolor(vor)
xlabel('Vorticity')
shading interp
axis equal
axis()
drawnow
end
end
end
close(gcf) :) = ' :'+')' 更正： S0=;
外力才起作用 Note: in LES, you need to use the special relationship between s_q and s_\nu so that the boundary locations is fixed.
[ 本帖最后由 luo@odu.edu 于 2013-5-30 16:51 编辑 ]
回复 4# luo@odu.edu 的帖子
没有仔细研究......为什么呢？用宏观量算是否没这个问题？感谢罗老师回答:D为什么3D下如果用 Lallemand & Luo 2000建议的we = 0, wej = 475/63, wxx = 0会出现流动方向上的抖动纹？
3D MRT/LBGK LES:
clc
clear all
close all
x=60
y=100
z=60
tau=0.5*ones(x,y,z);
step=1
u0=.1;
%% D3Q19 setup
ex(1,1,1,1:19)=;
ey(1,1,1,1:19)=;
ez(1,1,1,1:19)=;
ee=;
w0=1/3;w1=1/18;w2=1/36;
w(1,1,1,1:19)=;
%% setup solid boundary
=meshgrid(1:y,1:x,1:z);
Solid=false(x,y,z);
Solid=Solid | (i-x/2).^2+(j-y/5).^2+(k-z/2).^2<(x/10)^2;
Solid(1,:,:)=true;
Solid(x,:,:)=true;
Solid(:,:,1)=true;
Solid(:,:,z)=true;
Fluid=~Solid;
Inlet=Fluid & j==1;
Outlet=j==y;
FluidAll=repmat(Fluid,);
%% used for streaming indexing
Ind=reshape(1:x*y*z*19,x,y,z,19);
for n=1:18
Ind(:,:,:,n+1)=circshift(Ind(:,:,:,n+1),ee(n+1,:));
end
%% used for bounceback indexing
Blank=false(x,y,z);
SolidBB1=cat(4,Blank, Solid,Blank,Solid,Blank,Solid,Blank, Solid,Blank,Solid,Blank,Solid,Blank,Solid,Blank,Solid,Blank,Solid,Blank);
SolidBB2=cat(4,Blank, Blank,Solid,Blank,Solid,Blank,Solid, Blank,Solid,Blank,Solid,Blank,Solid,Blank,Solid,Blank,Solid,Blank,Solid);
temp=Ind(SolidBB1);Ind(SolidBB1)=Ind(SolidBB2);Ind(SolidBB2)=temp;
%% viscosity
Re=2e5;
Cs=0.02
viscosity=u0*y/Re
tau0=3*viscosity+.5;
%% LES setup
exx=ex.*ex;eyy=ey.*ey;ezz=ez.*ez;exy=ex.*ey;eyz=ey.*ez;exz=ex.*ez;
%% MRT setup
ex_=squeeze(ex);ey_=squeeze(ey);ez_=squeeze(ez);
ex2=ex_.*ex_;ey2=ey_.*ey_;ez2=ez_.*ez_;e2=ex2+ey2+ez2;e4=e2.*e2;
M=;
Mt=M';
M0=diag(1./sum(M.*M))*Mt;
S0=;
S=ones(x*y*z,1)*S0;
St=1./tau(:);
%% UI control
restart = uicontrol('position',,'style','toggle','string','restart');
quit = uicontrol('position',,'style','toggle','string','close');
pausebutton = uicontrol('position',,'style','toggle','string','pause');
bgkmrt = uicontrol('position',,'style','toggle','string','switch to bgk');
fastslow = uicontrol('position',,'style','toggle','string','fast simulation');
cbar = uicontrol('position',,'style','toggle','string','show colorbar');
%%
while get(quit,'value') == 0
set(restart,'value',0);
set(pausebutton,'value',0);
%% initial condition
nTimestep=0;
ux=zeros(x,y,z);
uy=u0*ones(x,y,z).*Fluid;
uz=zeros(x,y,z);
forcex=zeros(x,y,z);
forcey(1:x,1:y,1:z)=2e-5 *Fluid;
forcez=zeros(x,y,z);
rho=ones(x,y,z);
usq=ux.*ux+uy.*uy+uz.*uz;
c=bsxfun(@times,ex,ux)+bsxfun(@times,ey,uy)+bsxfun(@times,ez,uz);
density=bsxfun(@times,rho,w) .* bsxfun(@plus,1+3*c+4.5*c.*c,-1.5*usq);
%% loop
while get(restart,'value')==0 && get(quit,'value')==0
while get(pausebutton,'value')==1 && get(restart,'value')==0 && get(quit,'value')==0
set(pausebutton,'string','continue');pause(.01);end;set(pausebutton,'string','pause');%% restart and quit are volatile, so extra test is required
%% calc macro
rho=sum(density,4);
ux=sum(bsxfun(@times,density,ex),4)./rho;
uy=sum(bsxfun(@times,density,ey),4)./rho;
uz=sum(bsxfun(@times,density,ez),4)./rho;
%% forcing
ux=ux+forcex./rho;
uy=uy+forcey./rho;
uz=uz+forcez./rho;
%% collision
usq=ux.*ux+uy.*uy+uz.*uz;
c=bsxfun(@times,ex,ux)+bsxfun(@times,ey,uy)+bsxfun(@times,ez,uz);
density_eq=bsxfun(@times,rho,w) .* bsxfun(@plus,3*c+4.5*c.*c,1-1.5*usq);
density_delta=density_eq-density;
if(get(bgkmrt,'value')==0)
set(bgkmrt,'string','switch to bgk')
moment=reshape(density,)*M; rho_=rho(:);jx_=ux(:).*rho_;jy_=uy(:).*rho_;jz_=uz(:).*rho_;usq_=jx_.*jx_+jy_.*jy_+jz_.*jz_;
moment_eq=;
moment_delta_1 =moment_eq(:,1 +1)-moment(:,1 +1);
moment_delta_9 =moment_eq(:,9 +1)-moment(:,9 +1);
moment_delta_11=moment_eq(:,11+1)-moment(:,11+1);
moment_delta_13=moment_eq(:,13+1)-moment(:,13+1);
moment_delta_14=moment_eq(:,14+1)-moment(:,14+1);
moment_delta_15=moment_eq(:,15+1)-moment(:,15+1);
else
set(bgkmrt,'string','switch to mrt')
end
%% eddy viscosity with Smagorinsky model
if(get(bgkmrt,'value')==0)
set(bgkmrt,'string','switch to bgk')
S1=S(:,1+1);
qxx=-(S1.*moment_delta_1/38+St.* moment_delta_9/2);
qyy=-(S1.*moment_delta_1/38-St.*(moment_delta_9-3*moment_delta_11)/4);
qzz=-(S1.*moment_delta_1/38-St.*(moment_delta_9+3*moment_delta_11)/4);
qxy=-1.5*St.*moment_delta_13;
qyz=-1.5*St.*moment_delta_14;
qxz=-1.5*St.*moment_delta_15;
q=reshape(sqrt(qxx.*qxx+qyy.*qyy+qzz.*qzz+2*(qxy.*qxy+qyz.*qyz+qxz.*qxz)),);
else
set(bgkmrt,'string','switch to mrt')
qxx=sum(bsxfun(@times,density_delta,exx),4);
qyy=sum(bsxfun(@times,density_delta,eyy),4);
qzz=sum(bsxfun(@times,density_delta,ezz),4);
qxy=sum(bsxfun(@times,density_delta,exy),4);
qyz=sum(bsxfun(@times,density_delta,eyz),4);
qxz=sum(bsxfun(@times,density_delta,exz),4);
q=sqrt(qxx.*qxx+qyy.*qyy+qzz.*qzz+2*(qxy.*qxy+qyz.*qyz+qxz.*qxz));
end
tau=tau0+.5*(sqrt(tau0.*tau0+36*Cs*q)-tau0);
%% collision cont'd
if(get(bgkmrt,'value')==0)
set(bgkmrt,'string','switch to bgk')
St=1./tau(:);
S(:,)=repmat(St,);
moment=moment+bsxfun(@times,moment_eq-moment,S);
density_temp=reshape(moment*M0,);
else
set(bgkmrt,'string','switch to mrt')
density_temp=density+bsxfun(@rdivide,density_delta,tau);
end
density(FluidAll)=density_temp(FluidAll);
%% stream & bounce back
density=density(Ind);
%%
nTimestep=nTimestep+1;
if( get(fastslow,'value')==1 )
set(fastslow,'string','slow simulation');
step=10;
else
set(fastslow,'string','fast simulation');
step=1;
end
if(~mod(nTimestep,step))
uxplot=uy(:,:,z/2);
uyplot=ux(:,:,z/2);
rho_plot=rho(:,:,z/2);
q_plot=q(:,:,z/2);
tau_plot=tau(x/2,y/2,z/2);
Mass=sum(rho(:));
umag=sqrt(uxplot.*uxplot+uyplot.*uyplot);
vmax=max(umag(:));
vel=umag;
vor=curl(uxplot,uyplot);
subplot(2,2,1)
pcolor(vel)
if(get(cbar,'value')==1)
set(cbar,'string','hide colorbar');
colorbar
else
set(cbar,'string','show colorbar')
end
xlabel('Velocity')
shading interp
axis equal
axis()
title(sprintf('%d, tau=%f, u_m_a_x=%f, Vorticity_m_a_x=%f, Mass=%f',nTimestep,tau_plot,vmax,max(vor(:)),Mass));
subplot(2,2,3)
pcolor(vor)
if(get(cbar,'value')==1)
set(cbar,'string','hide colorbar');
colorbar
else
set(cbar,'string','show colorbar')
end
xlabel('Vorticity')
shading interp
axis equal
axis()
subplot(2,2,2)
pcolor(rho_plot);
if(get(cbar,'value')==1)
set(cbar,'string','hide colorbar');
colorbar
else
set(cbar,'string','show colorbar')
end
xlabel('Rho')
shading interp
axis equal
axis()
subplot(2,2,4)
pcolor(q_plot)
if(get(cbar,'value')==1)
set(cbar,'string','hide colorbar');
colorbar
else
set(cbar,'string','show colorbar')
end
xlabel('Local strain rate')
shading interp
axis equal
axis()
drawnow
end
end
end
close(gcf) First, the parameters we suggested was based on LINEAR analysis, and the flow you are simulating is NONLINEAR. The linear analysis is only valid for linear STABILITY and nothing else.
Second, boundary conditions are a separate issue. Perhaps you can read one of my recent paper:
L.-S. Luo, W. Liao, X. Chen, Y. Peng, and W. Zhang.
Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations.
Physical Review E 83(5):056710 (May 2011).
Which you can also download from my website. 非常好的程序，非常好的楼主，谢谢！2D的已经用起来了:lol
谢谢楼主分享
感谢楼主分享，还有Luo老师的讲解，：） 本帖最后由 bitxyq 于 2015-5-13 02:51 编辑:handshake ~~
帮助很大 请问有没有C++代码实现的呢 新手学习
页:
[1]