VBAで粒子法を勉強中

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

【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


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


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