VBAで粒子法を勉強中

メッシュレスで数値解析できるという、粒子法に興味を持ち、勉強中です。 プログラムのスキルが無い為、ExcelとVBAでのプログラム作成を目指しています。

【SPH粒子法のプログラムをVBAで作ってみる】 8.一連の計算のまとめと粒子変位の表示

一連の計算のまとめと粒子変位の表示

 粒子の圧力計算、力の計算、変位の計算までの計算モジュールを記述したので、これをルーチンとしてまとめます。さらに、エクセル上で粒子の位置を表示するモジュールを作成します。

 計算モジュールと、表示モジュールを1つのルーチンとし、これを繰り返します。





 上記をプログラムにすると、以下の様になります。



'===========================
'粒子の圧力、力、変位の計算までのルーチン化:
'===========================
Private Sub process(Particles() As Particle)

  Call calculate_density_and_pressure(Particles())
  Call calculate_force(Particles())
  Call calculate_position(Particles())

End Sub





'===========================
'粒子の変位の表示:
'===========================
Private Sub output_particles(Particles() As Particle)

  Dim num As Integer  'カウンタ
  Dim iP As Integer   '粒子No


'---計算ループの7回に1度実行---------
If LOOPs Mod 7 = 1 Then
        
  '---カウンタをリセット---------
    num = 0
 
  
'---キャンバス代わりに白い枠無し長方形を描きます---------
'  (アクティブなシートに、400x200の大きさです)

    sN = ActiveSheet.Name
    ActiveSheet.Shapes.AddShape(msoShapeRectangle, 1, 1, 400, 200).Select
    With Selection.ShapeRange.Fill
        .Visible = msoTrue
        .ForeColor.ObjectThemeColor = msoThemeColorBackground1
        .ForeColor.TintAndShade = 0
        .ForeColor.Brightness = 0
        .Transparency = 0
        .Solid
    End With
 
    Selection.ShapeRange.Line.Visible = msoFalse
        

    '---全ての粒子に対して実行---------
      
      For iP = 0 To nP - 1
       
    
    '--粒子を直径4の楕円のシェイプで示します---
    '  X座標は20+5X 、 Y座標は190-5X
        
        ActiveSheet.Shapes.AddShape(msoShapeOval _
        , 20 + 5 * Particles(iP).r.X - 2 _
        , 190 - 5 * Particles(iP).r.Y - 2 _
        , 4 _
        , 4 _
       ).Select
        With Selection.ShapeRange.Fill
            .Visible = msoTrue
           .ForeColor.RGB = RGB(0, 0, 255)
            .Transparency = 0
            .Solid
        End With
        Selection.ShapeRange.Line.Visible = msoFalse
          
    
      
    Next iP


'---シェイプをカットしてピクチャ貼付け---
    '---"描画"というシートに貼付けます---
    ActiveSheet.Shapes.SelectAll
    Selection.Cut
    Sheets("描画").Select
    Range("A1").Select
    
    '---ピクチャを最背面に移動します---
    ActiveSheet.PasteSpecial Format:="図 (拡張メタファイル)", Link:=False, _
        DisplayAsIcon:=False
    Selection.ShapeRange.ZOrder msoSendToBack
        
    '---元のアクティブシートに復帰---
    Sheets(sN).Select


''======================================
End If
  LOOPs = LOOPs + 1


End Sub



'===========================
'ルーチンを再帰呼び出しで回す:
'===========================
Private Sub smain(Particles() As Particle, iLoop As Integer)

  If iLoop > 0 Then

    '---初回のルーチン---
    Call process(Particles)
    
    '---描画ルーチン---
    Call output_particles(Particles)
    
    '---ここが再帰部分---
    Call smain(Particles(), iLoop - 1)

  Else
     Exit Sub
  End If
End Sub


 表示ルーチンは、エクセルのシェイプを利用して、長方形の上に粒子とみなした楕円を描き、これをピクチャとして貼り付けています。また、ルーチンの繰り返しは、再帰呼び出しで実現しています。

 一連のプログラムができましたので、次はメインモジュールでこれを動かしてみます。



にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

【SPH粒子法のプログラムをVBAで作ってみる】 7.各粒子の加速度、速度、変位の計算

【SPH粒子法のプログラムをVBAで作ってみる】 7.各粒子の加速度、速度、変位の計算

各粒子の加速度、速度、変位の計算モジュール

 各粒子に作用する力と粒子質量から粒子の加速度を求め、さらに重力加速度を加える事で、最終的な粒子加速度を求めます。

 但し粒子が壁にめりこんでいる場合は跳ね返りの処理をする必要があります。これは、求めた加速度の符号を入替えて減衰係数等を乗じる事で表現しています。なお、壁は凝った形ではなく単純にX軸(小さい方、大きい方)Y軸(小さい方、大きい方)で表される矩形のものとして扱っております。この場合、めり込み判定は単に粒子のXY各座標を壁の座標と比較する事で可能です。この部分はやや冗長にコードを記述した為、本モジュールの大半は、壁に対する跳ね返りの計算に使われています。

 加速度が決まったら、これに刻み時間を乗じて速度とし、この速度に刻み時間を乗じて粒子の変位としています。(いわゆるオイラー法での近似計算をしています)





 上記をプログラムにすると、以下の様になります。



'===========================
'粒子の加速度、速度、変位の計算:
'===========================

Private Sub calculate_position(Particles() As Particle)


  Dim accel As Vec2D  '加速度
  Dim g As Vec2D      '重力加速度
  Dim Norm As Vec2D   '法線ベクトル
  
  Dim speed As Double  '速度
  Dim diff As Double   '壁へのめりこみ
  Dim adj As Double    '壁からの反力
  
  Dim iP As Integer  'カウンタ
  
  '---重力加速度の代入---
  g.X = 0
  g.Y = -9.8
   
   
  '---全粒子で計算---

   For iP = 0 To nP - 1


  '---加速度を計算---
  ' 力/粒子質量 の筈ですが、プログラムは力*粒子質量です。
  ' 参考としたプログラムでも // why "*"? とあり
  '   理由は確認中です
    
    accel.X = Particles(iP).F.X * SPH_PMASS
    accel.Y = Particles(iP).F.Y * SPH_PMASS
    
  '---速度の変数に加速度^2を代入---
    speed = Norm2(accel)
     

    '---加速度^2が下限速度^2以上なら計算実施---
     If speed > SPH_LIMIT * SPH_LIMIT Then
     
        '---下記の式:下限速度/|速度| の意味は確認中---
        accel.X = accel.X * (SPH_LIMIT / Sqr(speed))
        accel.Y = accel.Y * (SPH_LIMIT / Sqr(speed))
     End If


   

    '---X軸壁の跳ね返り計算---

    '---座標の小さい側のX軸壁を計算---
    
        '---壁粒子の直径と粒子X座標の差(めりこみ量)を計算---
        diff = 2 * RADIUS - (-Particles(iP).r.X + MIN.X) * SPH_SIMSCALE
    
        '---粒子のめりこみが計算下限値以上なら計算実施---
         If diff > EPSILON Then
    
            '---壁への法線(単位Xベクトル)の値を代入---
            Norm.X = 1
            Norm.Y = 0
            
            '---反力値を計算:1項目は硬さ・めりこみ、2項目は減衰係数・(法線ベクトル・速度ベクトル内積)---
            adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP _
            * (Norm.X * Particles(iP).V.X + Norm.Y * Particles(iP).V.Y)
            
            '---加速度に反力ベクトル(反力値・法線ベクトル)を加算---
            accel.X = accel.X + adj * Norm.X
            accel.Y = accel.Y + adj * Norm.Y
            
        End If
    
    '---座標の大きい側のX軸壁を計算---
        
        '---壁粒子の直径と粒子X座標の差(めりこみ量)を計算---
        diff = 2 * RADIUS - (-MAX.X + Particles(iP).r.X) * SPH_SIMSCALE

        '---粒子のめりこみが計算下限値以上なら計算実施---
        If diff > EPSILON Then


            '---壁への法線(単位Xベクトル)の値を代入---
            Norm.X = -1
            Norm.Y = 0
            
            '---反力値を計算:1項目は硬さ・めりこみ、2項目は減衰係数・(法線ベクトル・速度ベクトル内積)---
            adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP _
            * (Norm.X * Particles(iP).V.X + Norm.Y * Particles(iP).V.Y)
            
            '---加速度に反力ベクトル(反力値・法線ベクトル)を加算---
            accel.X = accel.X + adj * Norm.X
            accel.Y = accel.Y + adj * Norm.Y
        End If
    

    '---Y軸壁の跳ね返り計算---

        '---座標の小さい側のY軸壁を計算---
        
        '---壁粒子の直径と粒子X座標の差(めりこみ量)を計算---
        diff = 2 * RADIUS - (-Particles(iP).r.Y + MIN.Y) * SPH_SIMSCALE
    
        '---粒子のめりこみが計算下限値以上なら計算実施---
        If diff > EPSILON Then
    
           '---壁への法線(単位Xベクトル)の値を代入---
            Norm.X = 0
            Norm.Y = 1
            
            '---反力値を計算:1項目は硬さ・めりこみ、2項目は減衰係数・(法線ベクトル・速度ベクトル内積)---
            adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP _
            * (Norm.X * Particles(iP).V.X + Norm.Y * Particles(iP).V.Y)
            
           '---加速度に反力ベクトル(反力値・法線ベクトル)を加算---
            accel.X = accel.X + adj * Norm.X
            accel.Y = accel.Y + adj * Norm.Y
        
        End If

        '---座標の大きい側のY軸壁を計算---
    
        diff = 2 * RADIUS - (-MAX.X + Particles(iP).r.X) * SPH_SIMSCALE
    
        '---粒子のめりこみが計算下限値以上なら計算実施---
        If diff > EPSILON Then
    
    
            '---壁への法線(単位Xベクトル)の値を代入---
            Norm.X = 0
            Norm.Y = -1
            
            '---反力値を計算:1項目は硬さ・めりこみ、2項目は減衰係数・(法線ベクトル・速度ベクトル内積)---
            adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP _
            * (Norm.X * Particles(iP).V.X + Norm.Y * Particles(iP).V.Y)
            
           '---加速度に反力ベクトル(反力値・法線ベクトル)を加算---
            accel.X = accel.X + adj * Norm.X
            accel.Y = accel.Y + adj * Norm.Y
            
        
        End If

   
    '---加速度に重力を加算---
     accel.X = accel.X + g.X
     accel.Y = accel.Y + g.Y
     
    
    '---加速度に刻み時間を乗じて速度を計算---
    Particles(iP).V.X = Particles(iP).V.X + accel.X * DT
    Particles(iP).V.Y = Particles(iP).V.Y + accel.Y * DT

    '---速度に刻み時間を乗じて変位を計算---
    'p->r += p->v * DT / SPH_SIMSCALE;
    Particles(iP).r.X = Particles(iP).r.X + Particles(iP).V.X * DT / SPH_SIMSCALE
    Particles(iP).r.Y = Particles(iP).r.Y + Particles(iP).V.Y * DT / SPH_SIMSCALE
    

  Next iP

End Sub


 ようやく、粒子の位置までたどり着きました。次は一連のモジュールをまとめるモジュールを記します。


にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

【SPH粒子法のプログラムをVBAで作ってみる】 6.粒子に作用する力の計算

粒子に作用する力の計算モジュール

 各粒子に作用する力(重力以外)を求めます。力には、圧力によるものと粘性によるものがあります。

 圧力は以下の式がベースになっています。

f:id:particlemethod:20170721053705p:plain

 ここで、カーネル関数の勾配である∇Wは、以下の様な式です。

f:id:particlemethod:20170721053945p:plain

 この関数の形は、以下のグラフの緑の所を見てください。粒子の距離が近いほど大きな値になっています。

f:id:particlemethod:20170721054055p:plain

 一方、粘性による力は以下の式がベースになっています。ポイントは速度差に比例している事です。(勾配を用いた式ですが、本プログラムではラプラシアンを用いた式を使用しています。これについては確認中です)

f:id:particlemethod:20170721054233p:plain

 ここで、カーネル関数のラプラシアンである∇2Wは、以下の様な式です。

f:id:particlemethod:20170721054337p:plain

 この関数の形は、以下のグラフの青の所を見てください。粒子の距離のマイナスに比例した値になっています。

f:id:particlemethod:20170721054646p:plain

 上記をプログラムにすると、以下の様になります。


'===========================
'粒子に作用する力の計算:
'===========================

Private Sub calculate_force(Particles() As Particle)


  Dim pterm As Double  '圧力項
  Dim vterm As Double  '速度項
  Dim r As Double      '位置
  Dim c As Double      '有効半径-粒子間距離
 
  Dim dr As Vec2D      '相手粒子との距離
  Dim force As Vec2D   '力
  Dim fcurr As Vec2D   '力(計算途中)
  
  Dim iP As Integer    'カウンタ
  Dim jP As Integer    'カウンタ
  
    
 '---全ての粒子に対して計算---
  For iP = 0 To nP - 1

    '---力をリセット---
    force.X = 0
    force.Y = 0


    '---全ての相手粒子に対して計算---
     For jP = 0 To nP - 1

      '---相手粒子が自身の場合除外---
       If iP = jP Then Exit For


      '---相手粒子との距離を計算---
       dr.X = (Particles(iP).r.X - Particles(jP).r.X) * SPH_SIMSCALE
       dr.Y = (Particles(iP).r.Y - Particles(jP).r.Y) * SPH_SIMSCALE

      '---距離ベクトルの大きさを計算---
      r = Norm(dr)
      
      '---有効半径内の場合に計算---
      If H > r Then

         
        '---圧力項を計算---
        c = H - r
        pterm = -0.5 * c * SpikyKern * (Particles(iP).p + Particles(jP).p) / r
        
        '---速度項を計算---
        vterm = LapKern * SPH_VISC
        
        '---圧力項x粒子間距離 と 速度項x粒子速度差 和を計算---
        fcurr.X = pterm * dr.X + vterm * (Particles(jP).V.X - Particles(iP).V.X)
        fcurr.Y = pterm * dr.Y + vterm * (Particles(jP).V.Y - Particles(iP).V.Y)
        
        
        '---上記値に密度を掛合せ力を算出する(密度Rhoは逆数1/ρになっている)---
        fcurr.X = fcurr.X * c * Particles(iP).Rho * Particles(jP).Rho
        fcurr.Y = fcurr.Y * c * Particles(iP).Rho * Particles(jP).Rho
        
        '---力を計算---
         force.X = force.X + fcurr.X
         force.Y = force.Y + fcurr.Y

     End If

    Next jP

   '---力の計算結果を粒子変数に代入---
    Particles(iP).F.X = force.X
    Particles(iP).F.Y = force.Y
  
  Next iP

End Sub


 まだ勉強中で、先に記述した式と上記のプログラムが完全にうまくリンクしていません。後で記しますが、それぞれ別の出典を参考にしている為です。いずれ整理して完全にリンクできればと考えています。ここでは理論については、雰囲気程度にご理解下さい。


にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

【SPH粒子法のプログラムをVBAで作ってみる】 5.密度と圧力の計算

密度と圧力の計算モジュール

 各粒子の粒子密度を求めます。粒子密度は、他の粒子との距離rに応じて重み係数Wを乗じたものになります。重み係数はカーネル係数とも呼ばれますが、SPH粒子法でよく使用されるのは6次関数で、以下の形をしています。

f:id:particlemethod:20170720035043p:plain

 変わった形の式ですが、この値を関数グラフにした場合、その勾配(グラフの傾き)やラプラシアン(グラフの曲率の様なもの)がシンプルであり各種計算に都合がよいのでSPH粒子法ではよく利用されます。グラフにすると、以下の様です。

f:id:particlemethod:20170720034934p:plain

 グラフ赤線が重み係数Wです。そして、SPH粒子法では、近傍粒子の距離rに応じた重み係数Wの値を合計し、これに粒子質量を乗じたものを粒子密度ρとします。


 求めた粒子密度ρと定常密度との間に差異がある場合、差異に応じて圧力Pが発生するものとします。

 上記をプログラムにすると、以下の様になります。


'===========================
'密度と圧力の計算:
'===========================

Private Sub calculate_density_and_pressure(Particles() As Particle)

    
Dim dr  As Vec2D    '相手粒子との距離ベクトル
Dim R2 As Double    '粒子間距離^2
Dim H2 As Double    '有効半径^2
Dim c As Double     '有効半径^2-粒子間距離^2

Dim Summ As Double  '計算和を保存する変数
Dim iP As Integer   '粒子計算に対するカウンタ
Dim jP As Integer   '粒子計算に対するカウンタ


'---有効半径の2乗---
H2 = H * H


'---全粒子で計算---

   For iP = 0 To nP - 1

    '---計算和をリセット---
      Summ = 0#
    
    '---相手粒子に対し計算---
    For jP = 0 To nP - 1
       
       '---自身の場合は除外---
       If iP = jP Then Exit For
       
       '---粒子間距離を計算---
       dr.X = (Particles(iP).r.X - Particles(jP).r.X) * SPH_SIMSCALE
       dr.Y = (Particles(iP).r.Y - Particles(jP).r.Y) * SPH_SIMSCALE
 
       '---|粒子間距離|^2を計算---
       R2 = Norm2(dr)

       '---有効半径内の場合に計算---
      If H2 > R2 Then

        '---Σ(H^2-R^2)^3を計算---
        c = H2 - R2
        Summ = Summ + c * c * c
      
      End If
      
     Next jP


    '-----------------------
    '重み関数 Wpoly6(r, H)を以下とし、
    
    '  Wpoly6(r, H)= α * |(H^2 - r^2)^3 |_〔if r<H〕
    
    '  ここで α=4/πH^8 〔at 2D〕, α=315/(64πH^9)〔at 2D〕
    
    '粒子密度ρ=Σm・Wpoly6(r, H)を計算。
    '-----------------------
    
    Particles(iP).Rho = Summ * SPH_PMASS * Poly6Kern
    
    
    '-----------------------
    '粒子圧力P=(粒子密度ρ-定常密度)・硬さ係数 を計算。
    '-----------------------
    
    Particles(iP).p = (Particles(iP).Rho - SPH_RESTDENSITY) * SPH_INTSTIFF
    

    
    '-----------------------
    '後の計算の為、粒子密度ρの逆数を保存
    '-----------------------
   
     If Particles(iP).Rho > 0 Then
       Particles(iP).Rho = 1 / Particles(iP).Rho
     End If

  Next iP

End Sub


 近傍粒子との計算に於いては、有効半径内の粒子を効率よく探す技法など、いくつか効率化の方法がありますが、ここではそうした策は講じず、シンプルなコードにしています。


にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

【SPH粒子法のプログラムをVBAで作ってみる】 4.初期設定モジュール

初期設定モジュール

 各パラメーターの読み込みや、配列の初期化を行うモジュールを記します。

 


'===========================
'初期化モジュール
'===========================
 
 Private Sub init(Particles() As Particle)

   Dim X As Double  '粒子のX位置
   Dim Y As Double  '粒子のY位置
   Dim d As Double  '粒子間隔
   Dim iP As Integer '粒子カウンタ
  
    MAX_LOOP = Range("最大ループ数").Value
    Pi = Range("円周率").Value
     
    H = Range("有効半径").Value
    SPH_PMASS = Range("粒子質量").Value

  
    SPH_INTSTIFF = Range("圧縮硬さ").Value
    SPH_EXTSTIFF = Range("境界圧縮硬さ").Value
    SPH_EXTDAMP = Range("境界速度減衰率").Value
    
    DT = Range("タイムステップ").Value
    RADIUS = Range("境界部粒子半径").Value
    EPSILON = Range("境界めりこみ判定").Value
      
    INITMIN.X = Range("水塊最小座標_X").Value
    INITMIN.Y = Range("水塊最小座標_Y").Value
    INITMAX.X = Range("水塊最大座標_X").Value
    INITMAX.Y = Range("水塊最大座標_Y").Value
    
    
    MIN.X = Range("境界最小座標_X").Value
    MIN.Y = Range("境界最小座標_Y").Value
    MAX.X = Range("境界最大座標_X").Value
    MAX.Y = Range("境界最大座標_Y").Value
    
    SPH_RESTDENSITY = Range("定常密度").Value
    SPH_PDIST = (SPH_PMASS / SPH_RESTDENSITY) ^ (1 / 3#)
    SPH_SIMSCALE = Range("スケール").Value
    SPH_VISC = Range("粘性係数").Value
    SPH_LIMIT = Range("上限速度").Value
  
    '----各カーネル係数を計算----
    '重み係数
    Poly6Kern = 4 / (Pi * H ^ 8)
    '3次元ではPoly6Kern = 315# / (64# * Pi * H ^ 9)
    
    '圧力係数
    SpikyKern = -30 / (Pi * H ^ 5)
    '3次元ではSpikyKern = -45# / (Pi * H ^ 6)
    
    '粘性係数(のラプラシアン)
    LapKern = 20 / (Pi * H ^ 5)
    '3次元ではLapKern = 45# / (Pi * H ^ 6)
    
    '---粒子間隔を設定---
    d = SPH_PDIST * 0.87 / SPH_SIMSCALE


'---粒子の位置・速度の初期値を設定---

iP = 0
Randomize
   
   For Y = INITMIN.Y To INITMAX.Y Step d
     For X = INITMIN.X To INITMAX.X Step d
         Particles(iP).r.X = X
         Particles(iP).r.Y = Y
         
         Particles(iP).r.X = Particles(iP).r.X - 0.05 + Rnd * 0.1 '揺らぎを与える
         Particles(iP).r.Y = Particles(iP).r.Y - 0.05 + Rnd * 0.1  '揺らぎを与える
         
         Particles(iP).V.X = 0
         Particles(iP).V.Y = 0
         
         iP = iP + 1
     Next X
  Next Y

'---粒子の数を保存---
nP = iP


End Sub


 沢山の変数がありますが、具体的な数値はエクセルのシートに範囲名も設定し、そのセルに記しています。この様にしておくと、各パラメータの設定がとてもやり易いです。エクセルVBAを使うメリットの1つです。

 各カーネル係数については、別途説明を試みたいと思います。

 また、参考にしたWEBページがありますので、これもその際にまとめて記します。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

【SPH粒子法のプログラムをVBAで作ってみる】 3.パブリック変数と関数の設定

パブリック変数と関数を設定する

 次に、モジュールの記述に先立ち、パブリック変数と関数の宣言をします。

先回設定した型を用いています。

 


'===========================
'パブリック変数宣言:
'===========================
Dim nP As Integer               '粒子数
Dim Particles(10000) As Particle '粒子変数(仮に1万粒子)
Dim MAX_LOOP As Integer         '最大ループ数

Dim LOOPs As Integer            'ループのカウンタ
Dim Pi As Double                '円周率

Dim H As Double                 '有効半径
Dim DT As Double                'タイムステップ

Dim SPH_PMASS As Double         '粒子質量
Dim SPH_INTSTIFF As Double      '圧縮硬さ
Dim SPH_EXTSTIFF As Double      '境界圧縮硬さ
Dim SPH_EXTDAMP As Double       '境界速度減衰率

Dim Poly6Kern As Double         '重み係数
Dim SpikyKern As Double         '圧力係数
Dim LapKern As Double           '粘性係数

Dim RADIUS As Double            '境界粒子半径
Dim EPSILON As Double           '境界めりこみ判定

Dim INITMIN As Vec2D            '初期の水塊の範囲(最小)
Dim INITMAX As Vac2D            '初期の水塊の範囲(最大)
Dim MIN As Vac2D                '境界の範囲(最小)
Dim MAX As Vac2D                '境界の範囲(最大)

Dim SPH_PDIST As Double         '初期粒子間隔
Dim SPH_SIMSCALE As Double      '実長/仮長 スケール係数
Dim SPH_RESTDENSITY As Double   '定常密度
Dim SPH_VISC As Double          '粘性係数
Dim SPH_LIMIT As Double         '上限速度


'===========================
'関数宣言:
'===========================

'ベクトルの長さ^2(ノルム^2)を算出する関数
Function Norm2(a As Vec2D) As Double
 Norm2 = a.X * a.X + a.Y * a.Y
End Function

'ベクトルの長さ(ノルム)を算出する関数
Function Norm(a As Vec2D) As Double
 Norm = Sqr(a.X * a.X + a.Y * a.Y)
End Function


 沢山の変数がありますが、内容は追い追い説明します。

 

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

【SPH粒子法のプログラムをVBAで作ってみる】 2.変数の型の設定

Type構文で変数の型を設定する

 まず、モジュールの記述に先立ち、専用の型宣言をします。

簡単の為、2次元で考えたいので、2次元のベクトル Vec2D と、粒子が保有するパラメータ Particle を記載します。

 



'=========================== 
'型宣言: 
'=========================== 
'----2次元ベクトルの型----- 
Type Vec2D 
 X As Double 'X成分 
 Y As Double 'Y成分 
End Type 

'----粒子の型----- 
Type Particle r As Vec2D '変位 
 V As Vec2D '速度 F As Vec2D '力 
 Rho As Double '密度ρ(ロー) 
 p As Double '圧力 
End Type 


 内容は、見て頂ければわかると思いますが、ポイントは、Particle の要素の中で、Vec2D タイプの変数がある所。つまり、VBAの型宣言は、入れ子構造にしても良いのですね。結構便利。

 

 

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村