Natural Convection & Unit Conversion

Hey guys,

I would be nice if you take a look at my calculations. This is my attempt to simulate 2D natural convection in an enclosure. I tried my best to follow Dr. Jonas Latt (2008) unit conversion paper, but my LBM results tend to be unstable (see image below). My calculations worked fine on arbitrary parameters based on Mohamad A.A. (2011) Pg. 97-98, so I believe I might have made a mistake in unit conversion. Please give me a tip!

Physical basics:


r_l   = 6093; %density [kg/m^3]
vis   = 1.81e-3/r_l; %kinematic viscosity [m^2/sec] 
gbeta = 9.81*1.20e-4; %gravity * thermal expansion coefficient [m/(sec^2.K)]
al    = 1.38e-5; %thermal diffusivity [m^2/sec]

domainHeight = 6.4e-2; %domainHeight [m]
n            = 64; %spatial resolution in vertical direction
dx           = domainHeight/n; %space step size [m]
totalTime    = 1020; %total time [sec]
iter         = 1020000; %total time iterations
dt           = totalTime/iter; %time step size [sec]

Thot        = 311; %hot wall temperature [K]
Tcold       = 301; %cold wall temperature [K]

Ra = gbeta*((Thot-Tcold)*domainHeight^3)/(vis*al); %Rayleigh number [phys]
Pr = vis/al; %Prandtl number [phys]

Unit conversions:


delta_x = 1; %space step size [LB unit]
L0      = dx/delta_x; %conversion factor – length [m]
delta_t = 1; %time step size [LB unit]
T0      = dt/delta_t; %conversion factor – time [sec]
rhoo    = 1; %density [LB unit]
M0      = r_l/rhoo*L0^3; %conversion factor – mass [kg]
K0      = (Thot-Tcold); %conversion factor – temperature [K]

LB units:


vis_LB = vis*T0/L0^2; %kinematic viscosity [LB unit]
al_LB  = al*T0/L0^2; %thermal diffusivity [LB unit]

gbeta_LB    = gbeta*K0*T0^2/L0; %gravity * thermal expansion coefficient
Thot_LB     = (Thot-Tcold)/K0;
Tcold_LB    = (Tcold-Tcold)/K0;
Tref_LB     = 0.5*(Thot_LB-Tcold_LB); %Tref for Boussinesq approximation

Ra_LB  = gbeta_LB*((Thot_LB-Tcold_LB)*(n*delta_x)^3)/(vis_LB*al_LB); %Rayleigh number [LB unit]
Pr_LB  = vis_LB/al_LB; %Prandtl number [LB unit]
Ma     = sqrt(gbeta_LB*(Thot_LB-Tcold_LB)*(n*delta_x))/(delta_x/delta_t/sqrt(3)); %Mach number [LB unit]

omega   = 1./(3*vis_LB  +0.5); %density relaxation parameter
omegat  = 1./(3*al_LB +0.5); %temperature relaxation parameter

Collision step:


cx  = [0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0]; %x-lattice velocity for D2Q9
cy  = [0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0]; %y-lattice velocity for D2Q9
w   = [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]; %weightage for D2Q9
cxg = [0.0,1.0,0.0,-1.0,0.0]; %x-lattice velocity for D2Q5
cyg = [0.0,0.0,1.0,0.0,-1.0]; %y-lattice velocity for D2Q5
wg  = [1/3, 1/6, 1/6, 1/6, 1/6]; %weightage for D2Q5

for k=1:9
term1        = u.*u+v.*v;
term2        = u.*cx(k)+v.*cy(k);
force(k,:,:) = 3.*w(k).*gbeta_LB.*(th-Tref_LB).*cy(k).*rho/(Thot_LB-Tcold_LB);
feq(k,:,:)   = rho.*w(k).*(1.0+3.0.*term2+4.50.*term2.*term2-1.50.*term1);
f(k,:,:)     = omega.*feq(k,:,:)+(1.-omega).*f(k,:,:)-delta_t*force(k,:,:);
end

for k=1:5
geq(k,:,:)  = th.*wg(k).*(1.0+3.0*(u.*cxg(k)+v.*cyg(k)));
g(k,:,:)    = omegat.*geq(k,:,:)+(1.-omegat).g(k,:,:);
end

Image - unstable LBM results:
https://1drv.ms/i/s!Aqf8s2sy-YOLbNNhEAZ2YQxD3MY
https://1drv.ms/i/s!Aqf8s2sy-YOLa220ATbLGDwpepY

I check my unit conversion by equating both Ra and Ra_LB as well as Pr and Pr_LB. It seems that current both values are equal, but instability still incur. On top of that, Ma=0.0468<0.1 is well within the incompressible limit. Please advise. I hope my query can help others with similar problems too!

Hi Victoria
I am also getting the same problem . I am even following A,A. Mohamad book but i am not able to get correct answer . Did you get it ?