VBAで粒子法を勉強中

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

【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


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


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