BGK / MRT turbulence in Matlab

citadel

用的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)=[0 1 0 -1 0 1 -1 -1 1];
ey(1,1,1:9)=[0 0 1 0 -1 1 1 -1 -1];
c0=[0,0];c1=[1,0];c2=[0,1];c3=[-1,0];c4=[0,-1];c5=[1,1];c6=[-1,1];c7=[-1,-1];c8=[1,-1];
w(1,1,1:9)=[4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36];

%% setup solid boundary
[i,j]=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,[1,1,9]);

%% 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)=[0 1 0 1 0 1 1 1 1];
eyy(1,1,1:9)=[0 0 1 0 1 1 1 1 1];
exy(1,1,1:9)=[0 0 0 0 0 1 -1 1 -1];

%% 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=[0,1.63,1.14,0,1.92,0,1.92];
% S0=[0,1.64,1.54,0,1.7,0,1.7];
S0=ones(x*y,1)*S0;

%% UI control
restart = uicontrol('position',[20 20 80 20],'style','toggle','string','restart');
quit = uicontrol('position',[120 20 80 20],'style','toggle','string','close');
pausebutton = uicontrol('position',[220 20 80 20],'style','toggle','string','pause');
bgkmrt = uicontrol('position',[320 20 80 20],'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,[x*y,9])*M;    ux_=ux(;uy_=uy(:);usq_=usq(:);
    moment_eq=bsxfun(@times,rho(:),[ones(x*y,1),-2+3*usq_,1-3*usq_,ux_,-ux_,uy_,-uy_,ux_.*ux_-uy_.*uy_,ux_.*uy_]);
   
    %% 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=[S0,1./tau(:),1./tau(:)];
    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,[x,y,9]);
%     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([1 y 1 x])
        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([1 y 1 x])
        drawnow
    end

end

end
close(gcf)

citadel

= ' :'+')'

citadel

更正： S0=[1,1.63,1.14,1,1.92,0,1.92];
外力才起作用

luo@odu.edu

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 编辑 ]

citadel

没有仔细研究......为什么呢？用宏观量算是否没这个问题？感谢罗老师回答

为什么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)=[0 1 -1 0  0 0  0 1 -1  1 -1 1 -1 -1  1 0  0  0  0];
ey(1,1,1,1:19)=[0 0  0 1 -1 0  0 1 -1 -1  1 0  0  0  0 1 -1  1 -1];
ez(1,1,1,1:19)=[0 0  0 0  0 1 -1 0  0  0  0 1 -1  1 -1 1 -1 -1  1];
ee=[squeeze(ex),squeeze(ey),squeeze(ez)];
w0=1/3;w1=1/18;w2=1/36;
w(1,1,1,1:19)=[w0, w1,w1,w1,w1,w1,w1, w2,w2,w2,w2,w2,w2,w2,w2,w2,w2,w2,w2];

%% setup solid boundary
[j,i,k]=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,[1,1,1,19]);

%% 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=[ones(19,1),19*e2-30,(21*e4-53*e2+24)/2,ex_,(5*e2-9).*ex_,ey_,(5*e2-9).*ey_,ez_,(5*e2-9).*ez_,3*ex2-e2,(3*e2-5).*(3*ex2-e2),ey2-ez2,(3*e2-5).*(ey2-ez2),ex_.*ey_,ey_.*ez_,ex_.*ez_,(ey2-ez2).*ex_,(ez2-ex2).*ey_,(ex2-ey2).*ez_];
Mt=M';
M0=diag(1./sum(M.*M))*Mt;
S0=[0,1.19,1.4,1,1.2,1,1.2,1,1.2,1,1.4,1,1.4,1,1,1,1.98,1.98,1.98];
S=ones(x*y*z,1)*S0;
St=1./tau(:);

%% UI control
restart = uicontrol('position',[20 20 80 20],'style','toggle','string','restart');
quit = uicontrol('position',[120 20 80 20],'style','toggle','string','close');
pausebutton = uicontrol('position',[220 20 80 20],'style','toggle','string','pause');
bgkmrt = uicontrol('position',[320 20 80 20],'style','toggle','string','switch to bgk');
fastslow = uicontrol('position',[420 20 80 20],'style','toggle','string','fast simulation');
cbar = uicontrol('position',[520 20 80 20],'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,[x*y*z,19])*M;    rho_=rho(:);jx_=ux(:).*rho_;jy_=uy(:).*rho_;jz_=uz(:).*rho_;usq_=jx_.*jx_+jy_.*jy_+jz_.*jz_;
        moment_eq=[rho_,-11*rho_+19*usq_,3*rho_-5.5*usq_,jx_,-2/3*jx_,jy_,-2/3*jy_,jz_,-2/3*jz_,3*jx_.*jx_-usq_,-1.5*(jx_.*jx_-usq_/3),jy_.*jy_-jz_.*jz_,-.5*(jy_.*jy_-jz_.*jz_),jx_.*jy_,jy_.*jz_,jx_.*jz_,zeros(x*y*z,3)];
        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)),[x,y,z]);
    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(:,[9+1,11+1,13+1,14+1,15+1])=repmat(St,[1,5]);
        moment=moment+bsxfun(@times,moment_eq-moment,S);
        density_temp=reshape(moment*M0,[x,y,z,19]);
    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([1 y 1 x])
        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([1 y 1 x])
        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([1 y 1 x])
        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([1 y 1 x])
        drawnow
    end

end

end
close(gcf)

luo@odu.edu

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.

bitxyq

非常好的程序，非常好的楼主，谢谢！2D的已经用起来了

erhu10

感谢楼主分享，还有Luo老师的讲解，：）

bitxyq

~~

lwkobe

帮助很大 请问有没有C++代码实现的呢 新手学习