![\"\"](\"http:\/\/81.70.49.155\/wp-content\/uploads\/2021\/12\/untitled.jpg\")
<\/p>\n <\/p>\n
\n\tFrom: \n<\/p>\n
https:\/\/palabos.unige.ch\/get-started\/lattice-boltzmann\/lattice-boltzmann-sample-codes-various-other-programming-languages\/
\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\n% thermalLB.m: Rayleigh Benard Convection, using a LB method,
\n% based on [Z.Guo, e.a., http:\/\/dx.doi.org\/10.1002\/fld.337].
\n% Boussinesq approximation is used for the buoyancy term:
\n% - Fluid is approximated with incompressible Navier-Stokes
\n% equations including a body force term, and simulated
\n% with a BGK model
\n% - Temperature is approximated with advection-diffusion
\n% equation and simulated with a BGK model
\n%
\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\n% Lattice Boltzmann sample, written in matlab
\n% Copyright (C) 2008 Andrea Parmigiani, Orestis Malaspinas, Jonas Latt
\n% Address: Rue General Dufour 24, 1211 Geneva 4, Switzerland
\n% E-mail: andrea.parmigiani@terre.unige.ch
\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\n% This program is free software; you can redistribute it and\/or
\n% modify it under the terms of the GNU General Public License
\n% as published by the Free Software Foundation; either version 2
\n% of the License, or (at your option) any later version.
\n% This program is distributed in the hope that it will be useful,
\n% but WITHOUT ANY WARRANTY; without even the implied warranty of
\n% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
\n% GNU General Public License for more details.
\n% You should have received a copy of the GNU General Public
\n% License along with this program; if not, write to the Free
\n% Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
\n% Boston, MA 02110-1301, USA.
\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<\/p>\n
clear<\/p>\n
% GENERAL FLOW CONSTANTS<\/p>\n
ly = 51;
\naspect_ratio = 2;
\nlx = aspect_ratio*ly;
\ndelta_x = 1.\/(ly-2);
\nPr = 1.;
\nRa = 20000.; % Rayleigh number
\ngr = 0.001; % Gravity
\nbuoyancy = [0,gr];<\/p>\n
Thot = 1; % Heating on bottom wall
\nTcold = 0; % Cooling on top wall
\nT0 = (Thot+Tcold)\/2;<\/p>\n
delta_t = sqrt(gr*delta_x);
\n% nu: kinematic viscosity in lattice units
\nnu = sqrt(Pr\/Ra)*delta_t\/(delta_x*delta_x);
\n% k: thermal diffusivity
\nk = sqrt(1.\/(Pr*Ra))*delta_t\/(delta_x*delta_x);
\nomegaNS = 1.\/(3*nu+0.5); % Relaxation parameter for fluid
\nomegaT = 1.\/(3.*k+0.5); % Relaxation parameter for temperature<\/p>\n
maxT = 80000; % total number of iterations
\ntPlot = 100; % iterations between successive graphical outputs
\ntStatistics = 10; % iterations between successive file accesses<\/p>\n
% D2Q9 LATTICE CONSTANTS
\ntNS = [4\/9, 1\/9,1\/9,1\/9,1\/9, 1\/36,1\/36,1\/36,1\/36];
\ncxNS = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
\ncyNS = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
\noppNS = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];<\/p>\n
% D2Q5 LATTICE CONSTANTS
\ntT = [1\/3, 1\/6, 1\/6, 1\/6, 1\/6];
\ncxT = [ 0, 1, 0, -1, 0];
\ncyT = [ 0, 0, 1, 0, -1];
\noppT = [ 1, 4, 5, 2, 3];<\/p>\n
[y,x] = meshgrid(1:ly,1:lx);<\/p>\n
% INITIAL CONDITION FOR FLUID: (rho=1, u=0) ==> fIn(i) = t(i)
\nfIn = reshape( tNS' * ones(1,lx*ly), 9, lx, ly);<\/p>\n
% INITIAL CONDITION FOR TEMPERATURE: (T=0) ==> TIn(i) = t(i)
\ntIn = reshape( tT' *Tcold *ones(1,lx*ly), 5, lx, ly);
\n% Except for bottom wall, where T=1
\ntIn(:,:,ly)=Thot*tT'*ones(1,lx);
\n% We need a small trigger, to break symmetry
\ntIn(:,lx\/2,ly-1)= tT*(Thot + (Thot\/10.));<\/p>\n
% Open file for statistics
\nfid = fopen('thermal_statistics.dat','w');
\nfprintf(fid,'Thermal Statistics: time-step --- uy[nx\/2,ny\/2] --- Nu\\n\\n\\n');<\/p>\n
% MAIN LOOP (TIME CYCLES)
\nfor cycle = 1:maxT
\n % MACROSCOPIC VARIABLES
\n rho = sum(fIn);
\n T = sum(tIn); %temperature
\n ux = reshape ( (cxNS * reshape(fIn,9,lx*ly)), 1,lx,ly) .\/rho;
\n uy = reshape ( (cyNS * reshape(fIn,9,lx*ly)), 1,lx,ly) .\/rho;<\/p>\n
% MACROSCOPIC BOUNDARY CONDITIONS
\n % NO-SLIP for fluid and CONSTANT at lower and upper
\n % boundary... periodicity wrt. left-right
\n % COLLISION STEP FLUID
\n for i=1:9
\n cuNS = 3*(cxNS(i)*ux+cyNS(i)*uy);
\n fEq(i,:,:) = rho .* tNS(i) .* ...
\n ( 1 + cuNS + 1\/2*(cuNS.*cuNS) - 3\/2*(ux.^2+uy.^2) );
\n force(i,:,:) = 3.*tNS(i) .*rho .* (T-T0) .* ...
\n (cxNS(i)*buoyancy(1)+cyNS(i)*buoyancy(2))\/(Thot-Tcold);
\n fOut(i,:,:) = fIn(i,:,:) - omegaNS .* (fIn(i,:,:)-fEq(i,:,:)) + force(i,:,:);
\n end<\/p>\n
% COLLISION STEP TEMPERATURE
\n for i=1:5
\n cu = 3*(cxT(i)*ux+cyT(i)*uy);
\n tEq(i,:,:) = T .* tT(i) .* ( 1 + cu );
\n tOut(i,:,:) = tIn(i,:,:) - omegaT .* (tIn(i,:,:)-tEq(i,:,:));
\n end<\/p>\n
% MICROSCOPIC BOUNDARY CONDITIONS FOR FLUID
\n for i=1:9
\n fOut(i,:,1) = fIn(oppNS(i),:,1);
\n fOut(i,:,ly) = fIn(oppNS(i),:,ly);
\n end<\/p>\n
% STREAMING STEP FLUID
\n for i=1:9
\n fIn(i,:,:) = circshift(fOut(i,:,:), [0,cxNS(i),cyNS(i)]);
\n end<\/p>\n
% STREAMING STEP FLUID
\n for i=1:5
\n tIn(i,:,:) = circshift(tOut(i,:,:), [0,cxT(i),cyT(i)]);
\n end<\/p>\n
% MICROSCOPIC BOUNDARY CONDITIONS FOR TEMEPERATURE
\n %
\n tIn(5,:,ly) = Tcold-tIn(1,:,ly)-tIn(2,:,ly)-tIn(3,:,ly)-tIn(4,:,ly);
\n tIn(3,:,1) = Thot-tIn(1,:,1) -tIn(2,:,1) -tIn(4,:,1) -tIn(5,:,1);<\/p>\n
% VISUALIZATION
\n if (mod(cycle,tStatistics)==0)
\n u = reshape(sqrt(ux.^2+uy.^2),lx,ly);
\n uy_Nu = reshape(uy,lx,ly); % vertical velocity
\n T = reshape(T,lx,ly);
\n Nu = 1. + sum(sum(uy_Nu.*T))\/(lx*k*(Thot-Tcold));
\n fprintf(fid,'%8.0f %12.8f %12.8f\\n',cycle,u(int8(lx\/2),int8(ly\/2))^2, Nu);
\n if(mod(cycle,tPlot)==0)
\n subplot(2,1,1);
\n imagesc(u(:,ly:-1:1)');
\n title('Fluid velocity');
\n axis off; drawnow
\n subplot(2,1,2);
\n imagesc(T(:,ly:-1:1)')
\n title(['Temperature (Nusselt number is ' num2str(Nu) ')']);
\n axis off; drawnow
\n end
\n end
\nend<\/p>\n
fclose(fid);<\/p>\n","protected":false},"excerpt":{"rendered":"