[009] 空気の流れをシミュレートするプログラム

 投稿者：mike  投稿日：2017年 3月12日(日)09時38分43秒
  　２次元で圧縮性流体の流れをシミュレートする実数型格子ガス法プログラム(009expansionRLG2D.bas)を公開します。
気体分子は非常に多いため、分子運動を追跡することでマクロな気体の流れをシミュレートするのには
多大な計算が必要になります。そこで、代表粒子を運動させ、それを統計処理することで、マクロな流れを
追跡することを考えます。
　代表粒子は位置rと速度vをもち、時間dt後には、r=v*dtに移動すると考えます。
空間は小さな格子状のセルに分割され、同じセル内の代表粒子は運動量とエネルギーを保存しながら
運動量とエネルギーを交換するとします。この際、交換のやり方は一意に決まらなくて、セル内の代表粒子を
同じ角度だけ回転することに相当する自由度が存在します。この自由度は各セルごとに乱数で決めます。
　本プログラムでは16万個の代表粒子と200x200の格子を使っています。つまり、各セルに平均4個の
代表粒子が存在することになります。統計的に意味のある平均の流れや密度、温度、圧力を計算するため、
5x5のセル（平均100個の代表粒子）について統計をとります。
　本プログラムでは箱の一部を真空にし、そこへの気体の膨張をシミュレートします。
この種の極端な問題は通常のNavier?Stokesを差分法で解く方法では結構難しいです。

表示の説明：
平均の流れは線分で、密度は円の大きさで、温度あるいは圧力は円の色で表現しています。
色は温度あるいは圧力が赤色に近いほど高く、青色に近いほど低いことを現します。（青色＜水色＜緑色＜黄色＜赤色）

試験環境：
本プログラムは十進BASIC 6.6.1 / macOS 10.7.5, 十進BASIC Ver 7.7.8 / windows 10　でテストしました。


------------------

!
! ========= real-coded lattice gas model 2D ==========
!
! 009expansionRLG2D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1   2017.03.12  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB rlg2d.setInitialCondition, rlg2d.moveParticles, rlg2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$

LET dispMode = 3 !3:mean flow+dens+temperature
CALL setInitialCondition

DO
   CALL moveParticles
   CALL drawParticles(dispMode)
   LET S$=INKEY$
   IF S$="." THEN EXIT DO
   IF S$="0" THEN LET dispMode = 0 !0:500 sample particles
   IF s$="1" THEN LET dispMode = 1 !1:mean flow
   IF s$="2" THEN LET dispMode = 2 !2:mean flow+dens
   IF s$="3" THEN LET dispMode = 3 !3:mean flow+dens+temperature
   IF s$="4" THEN LET dispMode = 4 !4:mean flow+dens+pressure
LOOP
END

EXTERNAL FUNCTION INKEY$
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- rlg2d module (real-coded lattice gas model 2D) ----------
!
! (1) move particles
!    r = r + v*dt
! (2) collision - all particles collide each other in the same cell
!    vm = mean v in the cell
!    vi = vm + rotate(theta)(vi0-vm),  (theta = PI/2 or -PI/2 : minmum viscosity)
!
MODULE rlg2d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, moveParticles, drawParticles
PUBLIC FUNCTION systemTime,systemTemp,totalKineticEnergy
SHARE NUMERIC Nsx,Nsy,NNs,NinCell,NNp,sysTime,t,dt,dx,dy,xMax,yMax,gravity
SHARE NUMERIC xx(400000),yy(400000) ! position of particle
SHARE NUMERIC vx(400000),vy(400000) ! velocity of particle
SHARE NUMERIC section(300,300,30)   ! section(Nsx,Nsy,NNs)
SHARE NUMERIC cellAttribute(300,300)! cellAttribute(Nsx,Nsy)
SHARE NUMERIC refTemp(4)
LET Nsx = 200    ! x-max number of section(x,y,n)
LET Nsy = 200    ! y-max number of section(x,y,n)
LET NNs = 30     ! n-max number of section(x,y,n)
LET NinCell = 4  !
LET NNp = Nsx*Nsy*NinCell
LET sysTime = 0.0
LET t = 0.0
LET dt = 1.0
LET dx = 1.0
LET dy = 1.0
LET xMax = Nsx*dx
LET yMax = Nsy*dy
LET gravity = 0.0 !e.g. -0.001 (gravitational force direction is -y direction)

! ---------- set initial condition

EXTERNAL SUB setInitialCondition
   DECLARE EXTERNAL FUNCTION cellAttributeAt
   DECLARE EXTERNAL SUB setColorMix
   LET t = 0.0
   FOR ic=0 TO Nsx-1
      FOR jc=0 TO Nsy-1
         LET cellAttribute(ic,jc) = 0
         IF (ic>=Nsx-Nsx/5 AND ic<=Nsx-Nsx/5+3 AND jc<Nsy/2) THEN
            LET cellAttribute(ic,jc) = 1
         END IF
      NEXT jc
   NEXT ic
   LET i=0
   DO WHILE (i<NNp)
      LET x = xMax*RND*0.799
      LET y = yMax*RND
      IF (cellAttributeAt(x,y)=0) THEN
         LET xx(i) = x
         LET yy(i) = y
         LET vx(i) = 2.0*(RND-0.5)
         LET vy(i) = 2.0*(RND-0.5)
         LET i = i+1
      END IF
   LOOP
   CALL setColorMix
   SET WINDOW 0,500,0,500
END SUB

EXTERNAL SUB setColorMix !modified decimal BASIC library
!  color pallet 100-180 blue(100) - green(140)  - red(180)
   FOR i=0 TO 80
      LET x = 8 - 0.1*i
      SELECT CASE x
      CASE 0 TO 2
         LET red=1
         LET green=x/2
         LET blue=0
      CASE 2 TO 4
         LET red=2-x/2
         LET green=1
         LET blue=0
      CASE 4 TO 6
         LET red=0
         LET green=1
         LET blue=x/2-2
      CASE 6 TO 8
         LET red=0
         LET green=4-x/2
         LET blue=1
      END SELECT
      SET COLOR MIX(100+i) red,green,blue
   NEXT i
END SUB

! ---------- move particles

EXTERNAL SUB moveParticles
   DECLARE EXTERNAL SUB movement,collision
   CALL movement
   CALL collision
   LET t = t + dt
END SUB

EXTERNAL SUB movement
   DECLARE EXTERNAL FUNCTION cellAttributeAt
   FOR i=0 TO NNp-1
      LET vy(i) = vy(i)+gravity*dt
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
      LET ca = cellAttributeAt(xx(i),yy(i))
      IF (ca>0) THEN
         LET xx(i) = xx(i)-vx(i)*dt
         LET yy(i) = yy(i)-vy(i)*dt
         LET r = 1
         IF (refTemp(ca)>0.0) THEN
            LET temp = 0.5*(vx(i)*vx(i)+vy(i)*vy(i))
            LET r = SQR(refTemp(ca)/temp)*(1.0-(RND-0.5))
         ELSE
            LET r = 1
         END IF
         LET vx(i) = -r*vx(i)
         LET vy(i) = -r*vy(i)
      END IF
   NEXT i
   !--- wall reflection
   LET rr = 1
   FOR i=0 TO NNp-1
      IF (xx(i)<0.0) THEN
         LET xx(i) = 0.0
         LET vx(i) = -rr*vx(i)
         LET vy(i) = rr*vy(i)
      END IF
      IF (xx(i) > xMax) THEN
         LET xx(i) = xMax
         LET vx(i) = -rr*vx(i)
         LET vy(i) = rr*vy(i)
      END IF
      IF (yy(i) < 0.0) THEN
         LET yy(i) = 0.0
         LET vx(i) = rr*vx(i)
         LET vy(i) = -rr*vy(i)
      END IF
      IF (yy(i) > yMax) THEN
         LET yy(i) = yMax
         LET vx(i) = rr*vx(i)
         LET vy(i) = -rr*vy(i)
      END IF
   NEXT i
END SUB

EXTERNAL FUNCTION cellAttributeAt(x,y)
   LET ix = INT(x/dx)
   IF ix>=Nsx THEN LET ix = Nsx-1
   IF ix<0 THEN LET ix = 0
   LET iy = INT(y/dy)
   IF iy>=Nsy THEN LET iy = Nsy-1
   IF iy<0 THEN LET iy = 0
   LET cellAttributeAt = cellAttribute(ix,iy)
END FUNCTION

EXTERNAL SUB collision
   DECLARE EXTERNAL SUB dividingSection,collisionInTheCell
   CALL dividingSection
   FOR ic=0 TO Nsx-1
      FOR jc=0 TO Nsy-1
         LET nn = section(ic,jc,0)
         IF nn>1 THEN
            IF cellAttribute(ic,jc)=0 THEN
               CALL collisionInTheCell(ic,jc)
            END IF
         END IF
      NEXT jc
   NEXT ic
END SUB

EXTERNAL SUB collisionInTheCell(ic,jc)
   LET cTh = 0
   IF RND<0.5 THEN LET sTh = -1 ELSE LET sTh = 1
   LET n = section(ic,jc,0)
   LET vxm = 0.0
   LET vym = 0.0
   FOR i=1 TO n
      LET k = section(ic,jc,i)
      LET vxm = vxm+vx(k)
      LET vym = vym+vy(k)
   NEXT i
   LET vxm = vxm/n
   LET vym = vym/n
   FOR i=1 TO n
      LET k = section(ic,jc,i)
      LET vxs = vx(k)-vxm
      LET vys = vy(k)-vym
      LET vx(k) = vxm +cTh*vxs +sTh*vys
      LET vy(k) = vym -sTh*vxs +cTh*vys
   NEXT i
END SUB

EXTERNAL SUB dividingSection
   FOR ic=0 TO Nsx-1
      FOR jc=0 TO Nsy-1
         LET section(ic,jc,0) = 0
      NEXT jc
   NEXT ic
   FOR ipp=0 TO NNp-1
      LET ic = INT(Nsx*xx(ipp)/xMax)
      IF ic>=Nsx THEN LET ic = Nsx-1
      LET jc = INT(Nsy*yy(ipp)/yMax)
      IF jc>=Nsy THEN LET jc = Nsy-1
      LET iq = section(ic,jc,0) + 1
      IF iq<NNs THEN
         LET section(ic,jc,0) = iq
         LET section(ic,jc,iq) = ipp
      END IF
   NEXT ipp
END SUB

EXTERNAL FUNCTION maxSection
   LET m = 0
   FOR ic=0 TO Nsx-1
      FOR jc=0 TO Nsy-1
         IF section(ic,jc,0)>m THEN LET m = section(ic,jc,0)
      NEXT jc
   NEXT ic
   LET maxSection = m
END FUNCTION

! ---------- utility

EXTERNAL FUNCTION systemTime
   LET systemTime = t
END FUNCTION

EXTERNAL FUNCTION totalKineticEnergy
   LET en = 0
   FOR i=0 TO NNp-1
      LET en = en + 0.5*(vx(i)*vx(i)+vy(i)*vy(i))
   NEXT i
   LET totalKineticEnergy = en
END FUNCTION

EXTERNAL FUNCTION systemTemp
   DECLARE EXTERNAL FUNCTION totalKineticEnergy
   LET systemTemp = totalKineticEnergy/NNp
END FUNCTION

! ---------- draw particles

EXTERNAL SUB drawParticles(drawMode)
   DECLARE EXTERNAL FUNCTION maxSection
   DECLARE EXTERNAL PICTURE boxWall,sampleParticles,meanFlow
   SET DRAW MODE HIDDEN
   CLEAR
   LET sc = 300/Nsx  !300 = boxsize in graphic window
   LET xp = 100
   LET yp = 100
   IF drawMode=0 THEN
      CALL boxWall(xp,yp,sc)
      CALL sampleParticles(500,xp,yp,sc)
   ELSEIF drawMode>=1 AND drawMode<=4 THEN
      CALL boxWall(xp,yp,sc)
      CALL meanFlow(drawMode-1,xp,yp,sc)
   END IF
   !draw caption
   SET TEXT HEIGHT 8
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50, 70 ,USING "time =####.#,   N =########":t,NNp
   PLOT TEXT, AT 50, 55 ,USING  "box size =####.# x ####.# ":xMax,yMax
   PLOT TEXT, AT 50, 40 ,USING "max section =####": maxSection
   PLOT TEXT, AT 50, 25 ,USING "total energy(=kinetic energy) =######.###########": totalKineticEnergy
   PLOT TEXT, AT 50, 10 :"real-coded lattice gas model 2D"
   SET TEXT HEIGHT 10
   PLOT TEXT, AT 30,480 :"'.':exit   '0':500-sample particles "
   PLOT TEXT, AT 30,460 :"'1':mean v         '2':mean v + dens"
   PLOT TEXT, AT 30,440 :"'3':mean v + Temp  '4':mean v + Pres"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL SUB boxWall(x0,y0,mag)
   DECLARE EXTERNAL FUNCTION cellAttributeAt
   SET LINE COLOR 1 ! black : wall
   PLOT LINES: x0,y0; x0+Nsx*mag,y0; x0+Nsx*mag,y0+Nsy*mag; x0,y0+Nsy*mag; x0,y0
   SET LINE COLOR 8 ! gray : wall
   FOR ic=1 TO Nsx
      FOR jc=1 TO Nsy
         IF cellAttributeAt(ic,jc)=1 THEN
            DRAW circle WITH SCALE(0.5)*SHIFT(x0+ic*mag,y0+jc*mag)
         END IF
      NEXT jc
   NEXT ic
END SUB

EXTERNAL SUB sampleParticles(nSample,x0,y0,mag)
   FOR i=1 TO nSample
      SET LINE COLOR 2 ! blue : particle
      DRAW circle WITH SCALE(2)*SHIFT(x0+xx(i)*mag,y0+yy(i)*mag)
      SET LINE COLOR 4 ! red : velocity
      PLOT LINES: x0+xx(i)*mag,y0+yy(i)*mag;
      PLOT LINES: x0+xx(i)*mag+vx(i)*mag*10,y0+yy(i)*mag+vy(i)*mag*10
   NEXT i
END SUB

EXTERNAL SUB meanFlow(mode,x0,y0,mag) !mode=0: mvx,mvy  1:0+dens  2:1+temp 3:1+pres
   DECLARE EXTERNAL FUNCTION cellAttributeAt
   LET nStep = 5
   SET LINE COLOR 2 ! blue : velocity
   FOR ic=1 TO Nsx STEP nStep
      FOR jc=1 TO Nsy STEP nStep
         LET n = 0
         LET mvx = 0
         LET mvy = 0
         FOR i=ic TO ic+nStep-1
            FOR j=jc TO jc+nStep-1
               LET nn = section(i,j,0)
               LET n = n + nn
               FOR ipp=1 TO nn-1
                  LET k = section(i,j,ipp)
                  LET mvx = mvx+vx(k)
                  LET mvy = mvy+vy(k)
               NEXT ipp
            NEXT j
         next i
         IF n>0 AND cellAttribute(ic+INT(nStep/2),jc+INT(nStep/2))=0 THEN
            LET mvx = mvx/n
            LET mvy = mvy/n
            LET dens = (n/(1.0*nStep*nStep))/NinCell
            !--- temp,pres
            IF mode>=2 THEN
               LET ke = 0
               FOR i=ic TO ic+nStep-1
                  FOR j=jc TO jc+nStep-1
                     LET nn = section(i,j,0)
                     FOR ipp=1 TO nn-1
                        LET k = section(i,j,ipp)
                        LET ke = ke + 0.5*((vx(k)-mvx)*(vx(k)-mvx)+(vy(k)-mvy)*(vy(k)-mvy))
                     NEXT ipp
                  NEXT j
               NEXT i
               LET temp = ke/n
               LET pres = dens*temp !here, no use
               IF mode=2 THEN
                  LET tempCol = INT(MIN(temp*160,80))
                  SET AREA COLOR 100+tempCol
               ELSEIF mode=3 THEN
                  LET presCol = INT(MIN(pres*160,80))
                  SET AREA COLOR 100+presCol
               END IF
            END IF
            !--- dens+temp plot
            IF mode>=1 THEN
               IF mode=1 THEN SET AREA COLOR 5 ! cyan
               DRAW disk WITH SCALE(dens*1.2*mag)*SHIFT(x0+(ic+0.5)*mag,y0+(jc+0.5)*mag)
            END IF
            !--- mvx,mvy plot
            IF ABS(mvx)>ABS(mvy) THEN
               IF mvx>=0 THEN SET LINE COLOR 2 ELSE SET LINE COLOR 4
            ELSE
               IF mvy>=0 THEN SET LINE COLOR 7 ELSE SET LINE COLOR 10
            END IF
            PLOT LINES: x0+(ic+0.5)*mag,y0+(jc+0.5)*mag;
            PLOT LINES: x0+(ic+0.5)*mag+mvx*mag*50,y0+(jc+0.5)*mag+mvy*mag*50
         END IF
      NEXT jc
   NEXT ic
END SUB

END MODULE
 

戻る