气体扩散模拟实验报告

气体扩散模拟实验报告

工程力学1001 陈金刚3100104568

一．实验背景：

要求用matlab编程模拟分子碰撞，演示气体扩散情况。本实验中的模型采用简化形式，所发生碰撞均为完全弹性碰撞。

壁面压强产生的原因是大量气体分子对容器壁的持续的、无规则撞击产生的。各个壁面所受力和压强为单位时间内的平均值来代替。由Ft=mv-（mv）=2mv，四个壁面在经历相同时间t的情况下，所受的平均压力F与2mv成正比。又四个壁面长度S一样，所以由P=F/S 可知，壁面所受压强P与所受平均压力成正比。

因此，只需统计分析在相同时间内，四个壁面所受碰撞的总的2mv 即可知道各个壁面所受压力和压强的情况。

二．实验基本情况说明：

区域尺寸：2维空间，200*200

粒子数目：N个，N可变化

粒子半径：各不相同，根据实验情况在一定范围内随机产生

粒子质量：各不相同，与粒子半径的平方成正比

初始位置：随机分布在左半区域

初始速度：随机大小

三．实验结果分析：（所有运行时间都一致）

由表格数据分析可知：

对于不同粒子数目，在运行足够长时间下，区域内分子运动会近似稳定，各个壁面所受压力，趋向稳定，在误差允许范围内，各个壁面所受压强可认为近似相同。

四．源代码：

N=50;%球数，可变

rad=rand(N,1)*2+4;%球半径，[4，6]

pos=[rand(N,1)*90+5 rand(N,1)*190+5];%初始位置：左半边区域

vel=rand(N,2)*20-10;%各球初始速度

color=rand(N,3);%各球颜色，随机产生

mass=5*rad.^2;%各球质量，与半径的平方成正比

Left=0; %左边界，下同

Right=200;

Up=200;

Down=0;

figure;

axis manual;

axis equal square; %固定坐标

axis([0 200 0 200]);

hold on;

%===============画不重叠的N个小球====

for i=2:N

index=0;

while index==0

index2=0;

for j=1:i-1

pp(j,:)=[pos(i,1)-pos(j,1),pos(i,2)-pos(j,2)];

rr(j)=rad(i)+rad(j);

ppp(j)=norm(pp(j,:));

if ppp(j)

%修正一开始若球相切时，后面判断可能误认为相撞

index2=1;break;

end

end

if ~index2

index=1;

else pos(i,:)=[rand()*90+5 rand()*190+5];

end

end

end

%==================================================================== couleft=0;%左壁面撞击的总mv，以下类似

couright=0;

couup=0;

coudown=0;

dt=0;%最小时间

qq=0;%用于循环次数控制

while qq<200

%=========以下各球两两碰撞最小时间计算

k=1;

for i=1:N

for j=i:N

dis=[pos(j,1)-pos(i,1),pos(j,2)-pos(i,2)];

vv=[vel(j,1)-vel(i,1),vel(j,2)-vel(i,2)];

dist=norm(dis);

rr=rad(i)+rad(j);

cosAlpha =abs( sqrt(1-(rr/dist)^2));

cosTheta=(dot(dis,vv)/norm(vv)/dist);

if cosTheta>=cosAlpha && cosTheta<1

dd=dist*cosTheta-sqrt(rr^2-(dist*sqrt(1-cosTheta^2))^2);

time(k)=dd/norm(vv);

k=k+1;

end

end

end

tball=min(time);

%================各球碰墙最小时间计算for i=1:N

if vel(i,1)>0

tx(i)=(Right-pos(i,1))/vel(i,1);

else

tx(i)=(Left-pos(i,1))/vel(i,1);

end

if vel(i,2) >0

ty(i)=(Up-pos(i,2))/vel(i,2);

else

ty(i)=(Down-pos(i,2))/vel(i,2);

end

end

twall=min(tx,ty);

%====

tBallWall=min(tball,twall);

dt=dt+tBallWall;

%与墙相撞改变速度

for i=1:N

if pos(i,1)-rad(i)<=Left && vel(i,1)~=0

pos(i,1)=rad(i);%修正边

couleft=couleft+mass(i)*abs(vel(i,1));

vel(i,1)=-vel(i,1);

elseif pos(i,1)+rad(i)>=Right && vel(i,1)~=0 pos(i,1)=200-rad(i);%修正边

couright=couright+mass(i)*abs(vel(i,1));

vel(i,1)=-vel(i,1);

end

if pos(i,2)-rad(i)<=Down && vel(i,2)~=0

pos(i,2)=rad(i);

coudown=coudown+mass(i)*abs(vel(i,2));

vel(i,2)=-vel(i,2);

elseif pos(i,2)+rad(i)>=Up && vel(i,2)~=0

pos(i,2)=200-rad(i);

couup=couup+mass(i)*abs(vel(i,2));

vel(i,2)=-vel(i,2);

end

end

%两球碰撞改变速度