Hello guys;
here is a matlab code for simulating rayleigh benard convection;now we want to use this code to simulate natural convection in square cavity;so wat are the changes should we do t o get best results ; Thank u so much ,best regards
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% thermalLB.m: Rayleigh Benard Convection, using a LB method, %
% based on [Z.Guo, e.a., http://dx.doi.org/10.1002/fld.337].
% Boussinesq approximation is used for the buoyancy term: %
% - Fluid is approximated with incompressible Navier-Stokes
% equations including a body force term, and simulated %
% with a BGK model
% - Temperature is approximated with advection-diffusion %
% equation and simulated with a BGK model
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, written in matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
% GENERAL FLOW CONSTANTS
ly = 51;
aspect_ratio = 2;
lx = aspect_ratio*ly;
delta_x = 1./(ly-2);
Pr = 1.;
Ra = 20000.; % Rayleigh number
gr = 0.001; % Gravity
buoyancy = [0,gr];
Thot = 1; % Heating on bottom wall
Tcold = 0; % Cooling on top wall
T0 = (Thot+Tcold)/2;
delta_t = sqrt(grdelta_x);
% nu: kinematic viscosity in lattice units
nu = sqrt(Pr/Ra)delta_t/(delta_xdelta_x);
% k: thermal diffusivity
k = sqrt(1./(PrRa))delta_t/(delta_xdelta_x);
omegaNS = 1./(3*nu+0.5); % Relaxation parameter for fluid
omegaT = 1./(3.*k+0.5); % Relaxation parameter for temperature
maxT = 80000; % total number of iterations
tPlot = 100; % iterations between successive graphical outputs
tStatistics = 10; % iterations between successive file accesses
% D2Q9 LATTICE CONSTANTS
tNS = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cxNS = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
cyNS = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
oppNS = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];
% D2Q5 LATTICE CONSTANTS
tT = [1/3, 1/6, 1/6, 1/6, 1/6];
cxT = [ 0, 1, 0, -1, 0];
cyT = [ 0, 0, 1, 0, -1];
oppT = [ 1, 4, 5, 2, 3];
[y,x] = meshgrid(1:ly,1:lx);
% INITIAL CONDITION FOR FLUID: (rho=1, u=0) ==> fIn(i) = t(i)
fIn = reshape( tNS’ * ones(1,lx*ly), 9, lx, ly);
% INITIAL CONDITION FOR TEMPERATURE: (T=0) ==> TIn(i) = t(i)
tIn = reshape( tT’ Tcold ones(1,lxly), 5, lx, ly);
% Except for bottom wall, where T=1
tIn(:,:,ly)=ThottT’ones(1,lx);
% We need a small trigger, to break symmetry
tIn(:,lx/2,ly-1)= tT(Thot + (Thot/10.));
% Open file for statistics
fid = fopen(‘thermal_statistics.dat’,‘w’);
fprintf(fid,‘Thermal Statistics: time-step — uy[nx/2,ny/2] — Nu\n\n\n’);
% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT
% MACROSCOPIC VARIABLES
rho = sum(fIn);
T = sum(tIn); %temperature
ux = reshape ( (cxNS * reshape(fIn,9,lxly)), 1,lx,ly) ./rho;
uy = reshape ( (cyNS * reshape(fIn,9,lxly)), 1,lx,ly) ./rho;
% MACROSCOPIC BOUNDARY CONDITIONS
% NO-SLIP for fluid and CONSTANT at lower and upper
% boundary… periodicity wrt. left-right
% COLLISION STEP FLUID
for i=1:9
cuNS = 3*(cxNS(i)ux+cyNS(i)uy);
fEq(i,:,
= rho . tNS(i) . …
( 1 + cuNS + 1/2*(cuNS.cuNS) - 3/2(ux.^2+uy.^2) );
force(i,:, = 3.tNS(i) .rho . (T-T0) . …
(cxNS(i)*buoyancy(1)+cyNS(i)buoyancy(2))/(Thot-Tcold);
fOut(i,:, = fIn(i,:, - omegaNS . (fIn(i,:,:)-fEq(i,:,:)) + force(i,:,:);
end
% COLLISION STEP TEMPERATURE
for i=1:5
cu = 3*(cxT(i)ux+cyT(i)uy);
tEq(i,:, = T . tT(i) . ( 1 + cu );
tOut(i,:, = tIn(i,:, - omegaT .* (tIn(i,:,:)-tEq(i,:,:));
end
% MICROSCOPIC BOUNDARY CONDITIONS FOR FLUID
for i=1:9
fOut(i,:,1) = fIn(oppNS(i),:,1);
fOut(i,:,ly) = fIn(oppNS(i),:,ly);
end
% STREAMING STEP FLUID
for i=1:9
fIn(i,:, = circshift(fOut(i,:,:), [0,cxNS(i),cyNS(i)]);
end
% STREAMING STEP FLUID
for i=1:5
tIn(i,:, = circshift(tOut(i,:,:), [0,cxT(i),cyT(i)]);
end
% MICROSCOPIC BOUNDARY CONDITIONS FOR TEMEPERATURE
%
tIn(5,:,ly) = Tcold-tIn(1,:,ly)-tIn(2,:,ly)-tIn(3,:,ly)-tIn(4,:,ly);
tIn(3,:,1) = Thot-tIn(1,:,1) -tIn(2,:,1) -tIn(4,:,1) -tIn(5,:,1);
% VISUALIZATION
if (mod(cycle,tStatistics)==0)
u = reshape(sqrt(ux.^2+uy.^2),lx,ly);
uy_Nu = reshape(uy,lx,ly); % vertical velocity
T = reshape(T,lx,ly);
Nu = 1. + sum(sum(uy_Nu.T))/(lxk*(Thot-Tcold));
fprintf(fid,’%8.0f %12.8f %12.8f\n’,cycle,u(int8(lx/2),int8(ly/2))^2, Nu);
if(mod(cycle,tPlot)==0)
subplot(2,1,1);
imagesc(u(:,ly:-1:1)’);
title(‘Fluid velocity’);
axis off; drawnow
subplot(2,1,2);
imagesc(T(:,ly:-1:1)’)
title(['Temperature (Nusselt number is ’ num2str(Nu) ‘)’]);
axis off; drawnow
end
end
end
fclose(fid);