格子玻尔兹曼MATLAB运用(LBGK_D2Q9_poiseuille_channel2D).docVIP

% Gianni Schena July 2005, schena@units.it % Lattice Boltzmann LBE, geometry: D2Q9, model: BGK % Application to permeability in porous media Restart=false % to restart from an earlier convergence logical(Restart); if Restart==false; close all, clear all % start from scratch and clean ... Restart=false; % type of channel geometry ; % one of the flollowing flags == true Pois_test=true, % no obstacles in the 2D channel % porous systems obs_regolare=false % obs_irregolare=false % tic % IN % |vvvv| + y % |vvvv| ^ % |vvvv| | - + x % OUT % Pores in 2D : Wet and Dry locations (Wet ==1 , Dry ==0 ) wXh_Dry=[3,1];wXh_Wet=[3,4]; if obs_regolare, % with internal obstacles A=repmat([zeros(wXh_Dry),ones(wXh_Wet)],[1,3]);A=[A,zeros(wXh_Dry)]; B=ones(size(A)); C=[A;B] ; D=repmat(C,4,1); D=[B;D] end if obs_irregolare, % with int obstacles A1=repmat([zeros(wXh_Dry),ones(wXh_Wet)],[1,3]); A1=[A1,zeros(wXh_Dry)] ; B=ones(size(A1)); C1=repmat([ones(wXh_Wet),zeros(wXh_Dry)],[1,3]); C1=[C1,ones(wXh_Dry)]; E=[A1;B;C1;B]; D=repmat(E,2,1); D=[B;D] end if ~Pois_test figure,imshow(D,[]) Channel2D=D; Len_Channel_2D=size(Channel2D,1); % Length Width=size(Channel2D,2); % should not be hod Channel_2D_half_Width=Width/2, end % test without obstacles (i.e. 2D channel no obstacles) if Pois_test %over-writes the definition of the pore space clear Channel2D Len_Channel_2D=36, % lunghezza canale 2d Channel_2D_half_Width=8; Width=Channel_2D_half_Width*2; Channel2D=ones(Len_Channel_2D,Width); % define wet area %Channel2D(6:12,6:8)=0; % put fluid obstacle imshow(Channel2D,[]); end [Nr Mc]=size(Channel2D); % Number rows and Munber columns % porosity porosity=nnz(Channel2D==1)/(Nr*Mc) % FLUID PROPERTIES % physical properties cs2=1/3; % cP_visco=0.5; % [cP] 1 CP Dinamic water viscosity 20 C density=1.; % fluid density Lky_visco=cP_visco/density; % lattice kinematic viscosity omega=(Lky_visco/cs2+0.5).^-1; % omega: relaxation frequency %Lky_visco=cs2