LBE.zip_LBE _sun

function [ux,uy,press,time,vel,velDifference,i] = anb(tmax,density,accel,omega,rrey,obst) % % ANB was originally written in FORTRAN % LBGK Lattice Boltzmann, Bhatnagar-Gross-Krook % % Flow between two plates with obstacles % The file par.m contains the parameters for the obstacles definitions % % Usage: % [ux,uy,press,time,vel,velDifference,i] = anb(tmax,density,accel,omega,rrey,obst) % where: % ux = matrix of x-velocity components % uy = matrix of y-velocity components % press = matrix of pressure values % i = number of time steps used lx=400; ly=200; tmax=3000; density=0.1; accel=0.005; omega=1/1.85; rrey=10; obst=zeros(lx,ly); obst(:,[1,ly])=1; % Node matrix (distribution matrix) node = zeros(lx,ly,9); nhlp = zeros(lx,ly,9); vel = 0; function [time, vel] = writeVelocity(lx,ly,time,obst,node,vel) % obtains the macroscopic velocity x = lx/2; nfree = 0; ux = 0; for y = 1:ly if( obst(x,y) == 0 ) dloc = 0; for i = 1:9 dloc = dloc + node(x,y,i); end ux(x,y) = ( node(x,y,1) + node(x,y,5) + node(x,y,8)... -( node(x,y,3) + node(x,y,6) + node(x,y,7) ) )/dloc; nfree = nfree + 1; end end vel = ux/nfree; function nhlp = propagate(lx,ly,node,nhlp) % nhlp = output matrix % moves the fluid particles to the next lattice node for x = 1:lx for y = 1:ly % Define north, south, east and west yn = mod(y,ly) + 1; xe = mod(x,lx) + 1; ys = ly - mod(ly+1-y, ly); xw = lx - mod(lx+1-x, lx); % density propagation nhlp(x ,y ,9) = node(x,y,9); nhlp(xe ,y ,1) = node(x,y,1); nhlp(x ,yn ,2) = node(x,y,2); nhlp(xw ,y ,3) = node(x,y,3); nhlp(x ,ys ,4) = node(x,y,4); nhlp(xe ,yn ,5) = node(x,y,5); nhlp(xw ,yn ,6) = node(x,y,6); nhlp(xw ,ys ,7) = node(x,y,7); nhlp(xe ,ys ,8) = node(x,y,8); end end function node = redistribute(lx,ly,obst,node,accel,density) % redistributes the density in the first lattice column % drives the motion by modifying the distribution in the first lattice % column % weighting factors t1 = density*accel/9; t2 = density*accel/36; for y = 1:ly if (obst(1,y) == 0 & node(1,y,3)-t1 >= 0 & node(1,y,6)-t2 >= 0 & ... node(1,y,7)-t2 >= 0 ) node(1,y,1) = node(1,y,1) + t1; node(1,y,3) = node(1,y,3) - t1; node(1,y,5) = node(1,y,5) + t2; node(1,y,6) = node(1,y,6) - t2; node(1,y,7) = node(1,y,7) - t2; node(1,y,8) = node(1,y,8) + t2; end end function node = relaxation(density, omega, lx,ly,node,nhlp,obst) % implements the BGK approximation. % Finds local equilibrium distributions and relaxes the distributions % towards equilibrium % Weighting factors in 9-bit LBGK t0 = 4/9; t1 = 1/9; t2 = 1/36; % c = speed of sound % csqu = c squared csqu = 1/3; % Main relaxation loop for x = 1:lx for y = 1:ly if (obst(x,y)==0) % initialize variable dloc dloc = 0; for i = 1:9 dloc = dloc + nhlp(x,y,i); end % ux and uy velocity components ux = (nhlp(x,y,1) + nhlp(x,y,5) + nhlp(x,y,8) -... (nhlp(x,y,3) + nhlp(x,y,6) + nhlp(x,y,7)))/ dloc; uy = (nhlp(x,y,2) + nhlp(x,y,5) + nhlp(x,y,6) -... (nhlp(x,y,4) + nhlp(x,y,7) + nhlp(x,y,8)))/ dloc; % usqu = magnitude of the velocity usqu = ux^2 + uy^2; % n-velocity components (n = lattice node connection vectors) un(1) = ux; un(2) = uy; un(3) = -ux; un(4) = -uy; un(5) = ux + uy; un(6) = -ux + uy; un(7) = -ux - uy; un(8) = ux - uy; % Equilibrium densities % index 0 in the fortran code is index 9 here nequ(9) = t0*dloc*(1 - usqu/(2*csqu)); nequ(1) = t1*dloc*(1 + un(1)/csqu + un(1)^2/(2*csqu^2) - usqu/(2*csqu)); nequ(2) = t1*dloc*(1 + un(2)/csqu + un(2)^2/(2*csqu^2) - usqu/(2*csqu)); nequ(3) = t1*dloc*(1 + un(3)/csqu + un(3)^2/(2*csqu^2) - usqu/(2*csqu)); nequ(4) = t1*dloc*(1 + un(4)/csqu + un(4)^2/(2*csqu^2) - usqu/(2*csqu)); nequ(5) = t2*dloc*(1 + un(5)/csqu + un(5)^2/(2*csqu^2) - usqu/(2*csqu)); nequ(6) = t2*dloc*(1 + un(6)/csqu + un(6)^2/(2*csqu^2) - usqu/(2*csqu)); nequ(7) = t2*dloc*(1 + un(7)/csqu + un(7)^2/(2*csqu^2) - usqu/(2*csqu)); nequ(8) = t2*dloc*(1 + un(8)/csqu + un(8)^2/(2*csqu^2) - usqu/(2*csqu)); % % relaxation towards the equilibrium for i= 1:9 node(x,y) = nhlp(x,y,i) + omega*(nequ(i) - nhlp(x,y,i)); end end end end function [ux,uy,press] = writeResults(lx,ly,obst,node,density) % obtains the macroscopic velocity components ux and uy and the pressure % press csqu = 1/3; % in the original fortran code the matrices ux, uy, press are written to a % file. Here the matrices ux, uy and press are passed to the function anb.m ux = zeros(lx,ly); uy = zeros(lx,ly); press = zeros(lx,ly); for y = 1:ly for x = 1:lx if (obst(x,y) == 1) ux(x,y) = 0; uy(x,y) = 0; press(x,y) = density*csqu; else dloc = 0; for i = 1:9 dloc = dloc + node(x,y,i); end % velocity components ux(x,y) = ( node(x,y,1) + node(x,y,5) + node(x,y,8)... -( node(x,y,3) + node(x,y,6) + node(x,y,7) ) )/dloc; uy(x,y) = ( node(x,y,2) + node(x,y,5) + node(x,y,6)... -( node(x,y,4) + node(x,y,7) + node(x,y,8) ) )/dloc; press(x,y) = dloc*csqu; end end end function node = initdensity(lx,ly,density) % initdensity(lx,ly,density,node) finds the initial density and a node-matrix % the function is run before any time steps and sets the initial % distributions t0 = density*4/9; t1 = density/9; t2 = density/36; for x = 1:lx for y = 1:ly node(x,y,9) = t0; node(x,y,1) = t1; node(x,y,2) = t1; node(x,y,3) = t1; node(x,y,4) = t1; node(x,y,5) = t2; node(x,y,6) = t2; node(x,y,7) = t2; node(x,y,8) = t2; end end % Initial densities node = initdensity(lx,ly,density); % % An iterative time step follows % write interim results [ux,uy,press] = writeResults(lx,ly,obst,node,density); % i = 1; velDifference(1) = 1; while i<tmax & velDifference(i) > 1/200000 % while i<tmax i = i + 1; if i >= 10 & mod(i, tmax/10) == 0 totalMass = checkDensity(lx,ly,node,i); disp(sprintf('Time step: %d Total density: %0.5g', i, totalMass)) % disp(sprintf('Node difference: %0.5g', nodeDiff(node, nhlp))) end %disp('Step'); disp(i); % A Control printout to find how long it takes to reach stationary flow % Compare the old velocity matrices to the next velocity matrices oldUx = ux; oldUy = uy; % end of control printout node = redistribute(lx,ly,obst,node,accel,density); nhlp = propagate(lx,ly,node,nhlp); for i=1:9 nhlp(i,obst)=node(opp(i),bbregion); end node = bounceback(lx,ly,obst,node,nhlp); node = relaxation(density, omega, lx,ly,node,nhlp,obst); % The difference in velocities between the two time steps [ux,uy,press] = writeResults(lx,ly,obst,node,density); velDifference(i) = sum(abs(ux - oldUx) + abs(uy - oldUy)) / sum(abs(ux) + abs(uy)); [time, vel] = writeVelocity(lx,ly,i,obst,node,vel); end %while [ux,uy,press] =