VBAで粒子法を勉強中

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

【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


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


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