DSMC蒙特卡洛算法

function [tauratio,Nbar,ux]=coufull(omega,Kn,Sw,Nc,Nppc,Nsteady) %% omega is viscosity power, vis=visref*(T/Tref)^omega %% if omega is negative, hard sphere collision model is used. %% otherwise vhs, for 0.0 < omega < 1. %% Couette flow, with given wall speed ratio Sw = Vw/sqrt(2RTw) %% and given Knudsen number = Kn = lambda/H, H = wall separation %% Nc cells, Nppc particles/cell (on average) %% therefore, number of paritcles is Nppc*Nc. %% at y = H plate has speed Vw, at y = 0 plate has speed 0 %% fspecref = fraction of specular reflection at walls, remainder diffuse %% dbstop if all error RTw = 1.0, rho1 = 1.0, H = 1.0 %% ARBITRARY UNITS %% these could be set to real values in MKS units %% RTw = ordinary gas constant (R)*wall temperature (Tw) R*Tw lambda_hs = %% Kn*H %% lambda_hs is a nominal mean free path = 32/(5*pi)*vis(rho*cmean) %% a slighty different nominal mean free path was used in J. Comput. Phys. %% v173, 600-619 (2001) %% there lambda_{nom}=2*vis/(rho*cmean). lambda_{hs}is equal to the actual %% mean free path for hard spheres, but not the 'actual mean free path' for %% VHS molecules. lambda_hs (in both cases) is a measure of the viscosity, %% density and temperature in the reference state. global Vw ymax; Vw = Sw*sqrt(2*RTw); ymax = H; % full height dy = ymax/Nc; nden1 = Nppc/vol; % molecules/unit volume global csd cmp; csd = sqrt(RTw); cmp = sqrt(2*RTw); % char. thermal speeds cmean1 = sqrt(8*RTw/pi); % a characteristic (mean) thermal speed if omega < 0.5 % use special case constant hard sphere routines omega = 0.5; HS = 1; VHS = 0; 'hard sphere'; sigma = 1/(sqrt(2)*nden1*lambda_hs); % = pi*d^2, d = molecular diameter % for hard sphere, upsilon = 0, and 6/Gamma(4-0) = 1; % same equations as VHS, see below tau_hs = lambda_hs/cmean1; tau1 = tau_hs; elseif omega <= 1 % if omega = 0.5, vhs routine will give hard sphere result HS = 0; VHS = 1; 'VHS', upsilon = omega - 0.5; omega, upsilon global sigmaref gref twoup sigmaref = (6/gamma(4-upsilon))/(sqrt(2)*nden1*lambda_hs); % OR: mass = rho1/nden1; mu_w = (5*pi/32)*mass*nden1*cmean1*lambda_hs; % sigmaref = (15*mass*sqrt(pi*RT_W))/(8*Gamma(4-upsilon)*mu_w); lambda_vhs = ((3 - upsilon)*(2 - upsilon)/6)*lambda_hs; % lambda_vhs = actual mean free path for VHS, less than nominal tau_vhs = lambda_vhs/cmean1; % expected mean free time for vhs collisions tau1 = tau_vhs; % not equal to tau_hs, used to check actual collision rate gref = sqrt(4*RTw); upsilon = omega - 0.5; twoup = 2*upsilon; else 'omega > 1', return; end % start at wall temperature y = rand(1,Nm)*ymax; vx = y/ymax*Vw + randn(1,Nm)*csd; % initial linear velocity profile vy = randn(1,Nm)*csd; vz = randn(1,Nm)*csd; tsteady = Nsteady*H/cmean1; tlimit = 3*tsteady; % long time required dt = tau1/4; % small time step maxstep = ceil(tlimit/dt) + 1; global nsteps TN_in c_in dtvol sigdtvol gmax gAmax twoup gmin; c_in = zeros(1,Nc); dtvol = dt/vol; if HS; sigdtvol = sigma*dt/vol; end gmax = 4.5*csd*ones(1,Nc); if VHS twoup = 2*upsilon; A = sigmaref*(gref./gmax).^twoup; % cross-section for max. velocity gAmax = gmax.*A; end TN_in = zeros(1,Nc); c_in = zeros(1,Nc); nsteps = 0; % for sampling Nbar global Ntot vxtot vytot vztot vxxtot vyytot vzztot % samples nsamp = 0; Ntot = zeros(1,Nc); vxtot = zeros(1,Nc); vytot = zeros(1,Nc); vztot = zeros(1,Nc); vxxtot = zeros(1,Nc); vyytot = zeros(1,Nc); vzztot = zeros(1,Nc); steady = 0; % starts as unsteady flow for step = 1:maxstep if ~steady & step*dt>tsteady % clear samples steady = 1; 'assumed steady state reached' % assumed steady now nsteps = 0; TN_in = zeros(1,Nc); c_in = zeros(1,Nc); end y = y + vy*dt; for m = 1:Nm % check boundaries, for all molecules if y(m) < 0 % reflect from stationary wall [y(m),vx(m),vy(m),vz(m)] = reflectwall0(y(m),vx(m),vy(m),vz(m)); end if y(m) > ymax % reflect from moving wall [y(m),vx(m),vy(m),vz(m)] = ... reflectwall1(y(m),vx(m),vy(m),vz(m)); end if y(m) < 0 | y(m) > ymax; 'y out of range', y(m), vy(m), dtr, end end pcell = ceil(y/dy); % cell # for this molecule, max(cell),min(cell) [plist,N_in,badd] = indx (pcell,Nc); if steady & mod(step+2,4) == 0 % samples one mean free time apart, before collisions [nsamp] = sample(vx,vy,vz,plist,badd,N_in,nsamp); % add to sample end if HS [vx,vy,vz] = HS(vx,vy,vz,plist,badd,N_in); % hard-sphere collisions else [vx,vy,vz] = VHS(vx,vy,vz,plist,badd,N_in); % VHS collisions end if steady & mod(step,4) == 0 % after collisions, add to sample. final average before and after [nsamp] = sample(vx,vy,vz,plist,badd,N_in,nsamp); if mean(Ntot) > 2e7, 'sample size large enough?', break; end end end cpp = (2*c_in./TN_in); % colls/particle for each cell < 1 figure,plot(cpp),ylabel('colls/particle/time step') crate = cpp/dt; % collision rate colls/particle/time tauratio = (1./crate)/tau1; Nbar = TN_in/nsteps; ux = vxtot./Ntot/Vw; % mean velocity ux/V_w in each cell cond = [Kn Sw Vw ymax/H RTw Nsteady]; save cond.txt cond -ascii; save nsamp.txt nsamp -ascii save Ntot.txt Ntot -ascii; save vxtot.txt vxtot -ascii save vytot.txt vytot -ascii; save vztot.txt vztot -ascii save vxxtot.txt vxxtot -ascii; save vyytot.txt vyytot -ascii save vzztot.txt vzztot -ascii; return