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!
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]
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]
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
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:
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!