Matlab气体模拟扩散编程代码

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

%陈金刚3100104568工程力学1001 修正版
clear all；
N=50;%球数，可变 
rad=rand（N1）*2+4;%球半径，[4，6] 
pos=[rand（N1）*90+5 rand（N1）*190+5];%初始位置：左半边区域 
vel=rand（N2）*36-18;%各球初始速度==可变
color=rand（N3）;%各球颜色，随机产生
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（i1）-pos（j1）pos（i2）-pos（j2）];
            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<1000  
 %=========以下各球两两碰撞最小时间计算
     k=1;
    for i=1:N
        for j=i:N
          dis=[pos（j1）-pos（i1）pos（j2）-pos（i2）];
          vv=[vel（j1）-vel（i1）vel（j2）-vel（i2）];
          dist=norm（dis）;
          rr=rad（i）+rad（j）;
          cosAlpha =abs（ sqrt（1-（rr/dist）^2））;
          cosTheta=（dot（disvv）/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（i1）>0
      tx（i）=（Right-pos（i1）+rad（i））/vel（i1）;%易错点，少了半径
   else 
      tx（i）=（Left-pos（i1）+rad（i））/vel（i1）;
    end
   if vel（i2） >0
      ty（i）=（Up-pos（i2）+rad（i））/vel（i2）;
   else 
      ty（i）=（Down-pos（i2）+rad（i））/vel（i2）;
   end
   end
   twall=min（txty）;
   %====
   tBallWall=min（tballtwall）;
  dt=dt+tBallWall;
  %====与墙相撞改变速度
for i=1:N
   if pos（i1）-rad（i）<=Left 
       pos（i1）=rad（i）;%修正边
     couleft=couleft+mass（i）*abs（vel（i1））;
       vel（i1）=-vel（i1）;
     elseif pos（i1）+rad（i）>=Right 
         pos（i1）=200-rad（i）;%修正边
        couright=couright+mass（i）*abs（vel（i1））;
       vel（i1）=-vel（i1）;
       end
   if  pos（i2）-rad（i）<=Down 
       pos（i2）=rad（i）;
     coudown=coudown+mass（i）*abs（vel（i2））;
       vel（i2）=-vel（i2）;
   elseif pos（i2）+rad（i）>=Up 
       pos（i2）=200-rad（i）;
       couup=couup+mass（i）*abs（vel（i2））;
       vel（i2）=-vel（i2）;
   end
end
%==========两球碰撞改变速度
for i=1:N
  for j=i+1:N
      twoBall=[pos（i1）-pos（j1）pos（i2）-pos（j2）];
      D=norm（twoBall）;
  if D-rad（i）-rad（j）<=0;
       if D          DRt=（rad（i）+rad（j）-D）/norm（[vel（i1）-vel（j1） vel（i2）-vel（j2）]）;
          pos（i1）=pos（i1）-DRt*vel（i1）;
          pos（i2）=pos（i2）-DRt*vel（i2）;
        

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件      422291  2012-09-06 13:45  气体扩散模拟实验报告.pdf
     文件        4157  2012-09-07 01:24  a.m