matlab-基于matlab的三维气体扩散模拟,包括速度分布、理想气体定律、真实气体、热流和扩散等-源码

clc; clear; close all; warning off; addpath(genpath(pwd)); diary dat.txt; NA = 6.022e23; % Avogadro's number kbJoules = 1.380658e-23; % Boltzmann's constant [Joules/K] Dalton = 1.6605402e-27; % [kg] Angstrom = 1e-10; % [meters] ps = 1e-12; % [seconds] kjpermol = 1000 / NA; % [1.66e-21 J] eV = 1.602e-19; % [Joules] Mewton = Dalton * Angstrom / ps^2; % [Newtons] Moule = Dalton * Angstrom^2 / ps^2; % [Joules] Mascal = Mewton / Angstrom^2; % [Pascals] Molvo = Angstrom^3; % [m^3] % Conversion Factors NewtonsPerMewton = Mewton; % [1.66054e-13 Newtons/Mewton] JoulesPerMoule = Moule; % [1.66054e-23 Joules/Moule] PascalsPerMascal = Mascal; % [1.66054e7 Pascal/Mascal] MewtonsPerNewton = 1/Mewton; % [6.022137377e12 Mewtons/Newton] MoulesPerJoule = 1/Moule; % [6.022137e22 Moules/Joule] MascalsPerPascal = 1/Mascal; % [6.022137e-8 Mascal/Pascal] AngstromsPerMeter = 1/Angstrom; KilogramsPerDalton = Dalton; CubicMetersPerMolvo = Angstrom^3; eVperMoule = Moule / eV; % [eV/Moule] eVperJoule = 1/eV; % [eV/Joule] AtmospheresPerPascal = 1/101325; % [Atm/Pa] AtmospheresPerMascal = AtmospheresPerPascal * PascalsPerMascal; % [Atm/Mascal] kb = kbJoules * MoulesPerJoule; % Boltzmann's constant [0.83145116 Moules/K] g = 9.8 * AngstromsPerMeter / (1e12)^2; % Acceleration of gravity [9.8e-14 A/ps^2] SpeedOfLight = 3e8 * 1e10/1e12; % [3e6 A/ps] nd = 3; % number of dimensions b3dDisplay = 1; % turn 3d display on / off % Define different atomic species Species(1).Name = 'Argon'; Species(1).Mass = 40; % [D] Species(1).Radius = 0.95; % [A] Species(1).Color = 'b'; % blue Species(1).Quantity = 0; Species(2) = Species(1); % just a copy of species 1 Species(2).Color = 'g'; % green % Diffusion2 %-------------------------------------------------- sExperiment='Diffusion'; bGasCollisions = 1; InitTemperature = 300; % [K] sPotentialType = 'none'; % hard spheres / ideal gas nPotentialType = 0; sParticleDistribution = 'diff2'; sVelocityDistribution = 'gaussian2'; DisplaySteps = 2; %round(1/TimeStep); % [steps/display] %L = 40; % [A] BoxSize = [20 20 8]; % flat box bMembrane = 0; np = 220; % Number of particles bVelocityHistogram = 0; TimeRange = 10; % [ps] TimeStep = 1/40; % [ps] nSpecies = 2; Species(1).Quantity = np*0.9; Species(2).Quantity = np*0.1; bShiftParticles = 0; Species(1).Temperature = InitTemperature; Species(2).Temperature = InitTemperature; SampleSteps = round(1/TimeStep); % number of time steps between samples. gives 1 sample per ps nx=10; Species(2).Name = 'dye'; Species(2).Mass = 800; % [D] Species(2).Radius = 4; % [A] % Container %----------------------------------------------------------------------------------- BoxMatrix = [BoxSize(1) 0 0; 0 BoxSize(2) 0; 0 0 BoxSize(3)]; BoxArea = 2 * (BoxSize(1) * BoxSize(2) + BoxSize(2) * BoxSize(3) + ... BoxSize(1) * BoxSize(3)); % [A^2] ie 2LH + 2HW + 2LW InitVolume = prod(BoxSize); % [A^3] % Calculated Parameters %----------------------------------------------------------------------------------- TotalTimeSteps = TimeRange / TimeStep; % total number of time steps during run SampleTime = TimeStep * SampleSteps; % amount of time elapsed between samples [ps] TotalSamples = TotalTimeSteps / SampleSteps; % total number of samples to take during run xrange = 0:(BoxSize(1)/nx):BoxSize(1); trange = 0:SampleTime:TimeRange; xrange = xrange(1:length(xrange)-1); trange = trange(1:length(trange)-1); % Assign data to arrays %----------------------------------------------------------------------------------- Radius = zeros(np,1); Mass = zeros(np,1); Color = zeros(np,1); Position = zeros(np,nd); Velocity = zeros(np,nd); Acceleration = zeros(np,nd); iStart = 1; for i = 1:nSpecies AtomMass = Species(i).Mass; AtomRadius = Species(i).Radius; AtomVolume = 4/3*pi*AtomRadius^3; % [A^3] kT = kb * Species(i).Temperature; % Thermal energy [Moules] VxRms = sqrt(kT / AtomMass); % std deviation of one velocity component [A/ps] Quantity = Species(i).Quantity; % number of atoms of this species iEnd = iStart + Quantity - 1; Range = iStart:iEnd; % range of indices for this species (eg 20:39) Radius(Range) = Species(i).Radius; Mass(Range) = Species(i).Mass; Color(Range) = Species(i).Color; Species(i).Range = Range; r = Species(i).Radius; BoxMatrixInner = BoxMatrix - 2 * r * eye(3); switch sParticleDistribution case 'half' Position(Range,:) = rand(Quantity,nd) * BoxMatrixInner + r; Position(Range,1) = Position(Range,1) * 0.5 + BoxSize(1)*(i-1)/2; % shift to left or right side of the box case 'test' % two particles symmetrically located in box Position(1,:) = [BoxSize(1)/4, BoxSize(2)/2, BoxSize(3)/3]; Position(2,:) = [3 * BoxSize(1)/4, BoxSize(2)/2, BoxSize(3)/3]; case 'fcc' % arrange particles evenly throughout box volume for ix=1:npod for iy=1:npod for iz=1:npod n = iz + npod*(iy-1) + npod*npod*(ix-1); Position(n,1) = ix-1; Position(n,2) = iy-1; Position(n,3) = iz-1; end end end Position = Position / npod * BoxMatrix + BoxSize(1)/npod/2; case 'diff2' if i==1 Position(Range,:) = rand(Quantity,nd) * BoxMatrixInner + r; Position(Range,3) = ones(Quantity,1) * BoxSize(3)/2; % z's in middle else theta = rand(Quantity,1) * 2 * pi; rpos = rand(Quantity,1) * 3; for j = 1:Quantity n = j+Range(1) - 1; Position(n,:) = [BoxSize(1)/2 + rpos(j)*cos(theta(j)), BoxSize(2)/2 + rpos(j)*sin(theta(j)), BoxSize(3)/2]; % Position(Range,:) = randn(Quantity,nd) * 5 + BoxSize(1)/2 end Position(Range,3) = ones(Quantity,1) * BoxSize(3)/2; % z's in middle end otherwise Position(Range,:) = rand(Quantity,nd) * BoxMatrixInner + r; end % Velocity Distribution switch sVelocityDistribution case 'gaussian' Velocity(Range,:) = VxRms * randn(Quantity,nd); for j=1:nd ActualRms = sqrt(mean(Velocity(Range,j) .* Velocity(Range,j))); % eg 1.96 should be 2.03 Scale = VxRms / ActualRms; Velocity(Range,j) = Velocity(Range,j) * Scale; end case 'gaussian2' Velocity(Range,:) = VxRms * randn(Quantity,nd); for j=1:nd ActualRms = sqrt(mean(Velocity(Range,j) .* Velocity(Range,j))); % eg 1.96 should be 2.03 Scale = VxRms / ActualRms; Velocity(Range,j) = Velocity(Range,j) * Scale; end Velocity(Range,3) = zeros(Quantity,1); % no z velocity case 'flat' Velocity(Range,:) = VxRms * (rand(Quantity,nd)-0.5); % flat distribution case 'same' % all have same speed but random directions Velocity(Range,:) = VxRms * sign(rand(Quantity,nd)-0.5); case 'none' Velocity(Range,:) = zeros(Quantity,nd); end iStart = iEnd + 1; end AverageIntermolecularSpacing = (InitVolume/np)^(1/nd); % [A] SpeedRms = sqrt(3 * kT / AtomMass); % [A/ps] ExpectedPressure = np * kT / InitVolume * AtmospheresPerMascal; % [atm] ExpectedCollisionRate = 4 * pi * sqrt(2) * AtomRadius^2 * np * SpeedRms / InitVolume * 1000; % [collisions/mlc/fs] ExpectedMeanFreeTime = 1 / ExpectedCollisionRate; % [ps] ExpectedMeanFreePath = SpeedRms * ExpectedMeanFreeTime; % [A] % Maxwell-Boltzmann probability distribution if bVelocityHistogram Speed = sqrt(kT / AtomMass); % [A/ps] s = Speed * (0:1/50:5); % Speed range ExpectedSpeedDistribution = np/2 * 4 * pi * (AtomMass/(2*pi*kT))^(3/2) * (s.^2) .* exp(-0.5 * AtomMass * s.^2 / kT); s2 = Speed * (-2:1/50:2); % 1 dim speed range ExpectedSpeedDistribution1d = np/2 * sqrt (AtomMass/(2*pi*kT)) * exp(-0.5 * AtomMass * s2.^2 / kT); end ExpectedMostProbableSpeed = sqrt(2 * kT / AtomMass); ExpectedAverageSpeed = sqrt(8 * kT / pi / AtomMass); ExpectedRmsSpeed = sqrt(3 * kT / AtomMass); VolumeFraction = AtomVolume * np / InitVolume; MassDensity = np * AtomMass / InitVolume * KilogramsPerDalton / CubicMetersPerMolvo; % [kg/m^3] fprintf('---------------------------------------------