www.pudn.com > cavity.rar > cavity.m, change:2011-04-20,size:3684b

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % cavity2d.m: 2D cavity flow, simulated by a LB method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Lattice Boltzmann sample, Matlab script % Copyright (C) 2006-2008 Jonas Latt % Address: Rue General Dufour 24, 1211 Geneva 4, Switzerland % E-mail: Jonas.Latt@cui.unige.ch % % Implementation of 2d cavity geometry and Zou/He boundary % condition by Adriano Sciacovelli %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This program is free software; you can redistribute it and/or % modify it under the terms of the GNU General Public License % as published by the Free Software Foundation; either version 2 % of the License, or (at your option) any later version. % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % You should have received a copy of the GNU General Public % License along with this program; if not, write to the Free % Software Foundation, Inc., 51 Franklin Street, Fifth Floor, % Boston, MA 02110-1301, USA. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear % GENERAL FLOW CONSTANTS lx = 128; ly = 128; uLid = 0.05; % horizontal lid velocity vLid = 0; % vertical lid velocity Re = 100; % Reynolds number nu = uLid *lx / Re; % kinematic viscosity omega = 1. / (3*nu+1./2.); % relaxation parameter maxT = 40000; % total number of iterations tPlot = 10; % cycles for graphical output % D2Q9 LATTICE CONSTANTS t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36]; cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1]; cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1]; opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7]; lid = [2: (lx-1)]; [y,x] = meshgrid(1:ly,1:lx); obst = ones(lx,ly); obst(lid,2:ly) = 0; bbRegion = find(obst); % INITIAL CONDITION: (rho=0, u=0) ==> fIn(i) = t(i) fIn = reshape( t' * ones(1,lx*ly), 9, lx, ly); % MAIN LOOP (TIME CYCLES) for cycle = 1:maxT % MACROSCOPIC VARIABLES rho = sum(fIn); ux = reshape ( (cx * reshape(fIn,9,lx*ly)), 1,lx,ly ) ./rho; uy = reshape ( (cy * reshape(fIn,9,lx*ly)), 1,lx,ly ) ./rho; % MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS ux(:,lid,ly) = uLid; %lid x - velocity uy(:,lid,ly) = vLid; %lid y - velocity rho(:,lid,ly) = 1 ./ (1+uy(:,lid,ly)) .* ( ... sum(fIn([1,2,4],lid,ly)) + 2*sum(fIn([3,6,7],lid,ly)) ); % MICROSCOPIC BOUNDARY CONDITIONS: LID (Zou/He BC) fIn(5,lid,ly) = fIn(3,lid,ly) - 2/3*rho(:,lid,ly).*uy(:,lid,ly); fIn(9,lid,ly) = fIn(7,lid,ly) + 1/2*(fIn(4,lid,ly)-fIn(2,lid,ly))+ ... 1/2*rho(:,lid,ly).*ux(:,lid,ly) - 1/6*rho(:,lid,ly).*uy(:,lid,ly); fIn(8,lid,ly) = fIn(6,lid,ly) + 1/2*(fIn(2,lid,ly)-fIn(4,lid,ly))- ... 1/2*rho(:,lid,ly).*ux(:,lid,ly) - 1/6*rho(:,lid,ly).*uy(:,lid,ly); % COLLISION STEP for i=1:9 cu = 3*(cx(i)*ux+cy(i)*uy); fEq(i,:,:) = rho .* t(i) .* ... ( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) ); fOut(i,:,:) = fIn(i,:,:) - omega .* (fIn(i,:,:)-fEq(i,:,:)); end % MICROSCOPIC BOUNDARY CONDITIONS: NO-SLIP WALLS (bounce-back) for i=1:9 fOut(i,bbRegion) = fIn(opp(i),bbRegion); end % STREAMING STEP for i=1:9 fIn(i,:,: ) = circshift(fOut(i,:,: ), [0,cx(i),cy(i)]); end % VISUALIZATION if (mod(cycle,tPlot)==0) u = reshape(sqrt(ux.^2+uy.^2),lx,ly); u(bbRegion) = nan; imagesc(u(:,ly:-1:1)'./uLid); colorbar axis equal off; drawnow end end