So-net無料ブログ作成

[境界要素法プログラムを設計する(境界積分サブルーティン:実装編)] [ネコ騙し数学]

[境界要素法プログラムを設計する(境界積分サブルーティン:実装編)]

1.積分パラメータ関数の実装

 

bem14-fig-001.png

 

 前々回、メインの中に与えた境界積分サブルーティン(関数)は、以下でした。

 

    REM  (x0,y0),(x1,y1),(x2,y2)から、要素kの境界積分値を計算
    REM  ψj1に関する結果をB1,qj1に関する結果をH1で表す
    REM  ψj2に関する結果をB2,qj2に関する結果をH2で表す

    Lk = Length(x1,y1,x2,y2)
    Beta_k = Angle(x2 - x1, y2 - y1, Lk)
    r1 = Length(x0,y0,x1,y1)
    r2 = Length(x0,y0,x2,y2)
    Gamm_1 = Angle(x1 - x0, y1 - y0, r1)
    Gamm_2 = Angle(x2 - x0, y2 - y0, r2)
    h = Normal_Length(x2 - x1, y2 - y1, x0,y0)
    s = Normal_Foot(x2 - x1, y2 - y1, x0,y0, h)

  B1 = B_Inte_B1(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
  B2 = B_Inte_B2(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
  H1 = B_Inte_H1(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
  H2 = B_Inte_H2(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)

 B1B2H1H2の計算前に必要な積分パラメータを全て揃えます。そこには、B_Inte_?の中身として[内点方程式の積分公式(ラプラス型)]と[内点方程式から境界方程式へ(ラプラス型)]で導出した式を単純に並べたいという意図があります。なのでまず、積分パラメータを与える関数を検討します。

 

 最初にLength関数ですが、明らかに次でOKです。

 

    External Function Length(x1,y1,x2,y2)
      Length = Sqr((x2 - x1)^2 + (y2 - y1)^2)
    End Function

 前々回に予告したように、関数とサブルーティンはローカル変数を持てる、外部手続きで統一します。Normal_LengthNormal_Footの計算に疑問の余地はないのですが、ちょっと長そうなので後にします。

 次にBeta_kの扱いですが、前回の検討からGamm_1Gamm_2の特殊ケースとして扱って良いとわかりました。という訳でGamm_1Gamm_2の計算です。これらは単独に計算しては駄目で、Gamm_2Gamm_1の形で計算する必要がありました。それをGam_12で表します。

 

 まず特異点(x0y0)が要素節点(x1y1)(x2y2)に一致しないケースから考え始めると、わかりやすいと思います。以下、前回を参照しつつ・・・。

 

 計算開始の前提として、図-1r1r2および(x0y0)(x1y1)(x2y2)は引数渡しされるとします。また特異点と要素節点の一致/不一致判定のために、0距離判定値Epsも必要です。

 

    External Function Angle(x0, y0, x1,y1,x2,y2, r1, r2, Eps)

      If r1 < Eps Then
        x01 = 1
        y01 = 0
        x02 = x2 - x1
        y02 = y2 - y1
        r1 = 1
      ElseIf r2 < Eps Then
        x01 = x1 - x2
        y01 = y1 - y2
        x02 = 1
        y02 = 0
        r2 = 1
      Else
        x01 = x1 - x0
        y01 = y1 - y0
        x02 = x2 - x0
        y02 = y2 - y0
      End If

      Gum_12 = Acos((x01 * x02 + y01 * y02) / r1 / r2)
      Sign = Sgn(x01 * y02 - x02 * y01)
      Gum_12 = Gum_12 * Sign

      Angle = Gum_12
    End Function

 ElseEnd Ifの設定が基本です。

 r1 < Epsr10)の時は、(x01y01)(10)とするのでした。この時(x0y0)(x1y1)に一致するとみなせるので、(x02y02)(x2 - x1y2 - y1)です。そしてr11とします。

 r2 < Epsr20)の時は、(x02y02)(10)とするのでした。この時(x0y0)(x2y2)に一致するとみなせるので、(x01y01)(x1 - x2y1 - y2)です。そしてr21とします。

 

 以上より、積分パラメータを与える関数呼び出しは、

 

    Lk = Length(x1,y1,x2,y2)
    r1 = Length(x0,y0,x1,y1)
    r2 = Length(x0,y0,x2,y2)
    Gam_12 = Angle(x0, y0, x1,y1,x2,y2, r1, r2, Eps)
    h = Normal_Length(x0, y0, x1,y1,x2,y2)
    s = Normal_Foot(x0, y0, x1,y1,x2,y2)

になります(Normal_LengthNormal_Footの引数が前回のままでは不味いと気付いた(^^;))。

 

 Normal_Lengthの実装は以下です。Normal_Lengthは図-1の垂線距離hを求めます。

a = x2 - x1

b = y2 - y1

とすると、(x1y1)(x2y2)を通る直線は、yb/a*xに平行です。a0の場合も考慮して、これをbxay0の形にし、(x1y1)を通る事を考慮すると要素kの直線の方程式は、b(xx1)a(yy1)0になります(なんか受験を思い出すなぁ~(^^))。直線の式からその垂線ベクトルの方位は、明らかに(b,-a)です。従って点(x0y0)からの垂線は、tを任意の実数として、

  

と書けます。(1)b(xx1)a(yy1)0に載る条件は、(1)を代入し、

  

 (2)tについて解くと、

  

 従ってhは、t(3)の値の時、

  

 

 tを利用すれば垂線の足もわかるので、図-1sも出せます。垂線の足の座標は、t(3)の値にした(1)です。よってsは、

  

から、

  

です。

 こうなるとNormal_LengthNormal_Footで同じパラメータtを利用したいですよね。こうしましょう。

 


    t = Normal_Param(x0, y0, x1,y1,x2,y2)
    h = Normal_Length(x1,y1,x2,y2, t)
    s = Normal_Foot(x0, y0, x1,y1,x2,y2, t)

    External Function Normal_Param(x0, y0, x1,y1,x2,y2)
      a = x2 - x1
      b = y2 - y1
      t = (b * (x1 - x0) - a * (y1 - y0)) / (a^2 + b^2)

      Normal_Param = t
    End Function

    External Function Normal_Length(x1,y1,x2,y2, t)
      a = x2 - x1
      b = y2 - y1
      h = Abs(t) * Sqr(a^2 + b^2)

      Normal_Length = h
    End Function

    External Function Normal_Foot(x0, y0, x1,y1,x2,y2, t)
      a = x2 - x1
      b = y2 - y1
      s = Sqr((t * b + x0 - x1)^2 + (-t * a + y0 - y1)^2)

      Normal_Foot = s
    End Function


2.境界積分関数の実装

 以上で境界積分に必要な積分パラメータは全て揃いました。

 

  B1 = B_Inte_B1(Lk, r1, r2, Gam_12, h, s, Eps)
  B2 = B_Inte_B2(Lk, r1, r2, Gam_12, h, s, Eps)
  H1 = B_Inte_H1(Lk, r1, r2, Gam_12, h, s, Eps)
  H2 = B_Inte_H2(Lk, r1, r2, Gam_12, h, s, Eps)


で良いはずです。引数の最後にEpsが入るのは、これらでも特異点と要素節点の一致/不一致判定が必要だからです。

 

 後はB_Inte_?の中に、[内点方程式の積分公式(ラプラス型)]と[内点方程式から境界方程式へ(ラプラス型)]で導出した式を並べるだけですが、関数の内部構造はExternal Function Angleと同一構造になるとわかります。なのでIfブロックの形を決め、「関数名 = 値」の行を書いてから、計算式を並べる事をお奨めします。



    External Function B_Inte_B1(Lk, r1, r2, Gam_12, h, s, Eps)

      If r1 < Eps Then
        B1 = Gam_12 / 2 /Pai
      ElseIf r2 < Eps Then
        B1 = 0
      Else
        B1 = Gum_12 - 1 / Lk * (s * Gum_12 + h * (Log(r2) - Log(r1)))
        B1 = B1 / 2 /Pai
      End If

      B_Inte_B1 = B1
    End Function

    External Function B_Inte_B2(Lk, r1, r2, Gam_12, h, s, Eps)

      If r1 < Eps Then
        B2 = 0
      ElseIf r2 < Eps Then
        B2 = Gam_12 / 2 /Pai
      Else
        B2 = 1 / Lk * (s * Gum_12 + h * (Log(r2) - Log(r1)))
        B2 = B1 / 2 /Pai
      End If

      B_Inte_B2 = B2
    End Function

    External Function B_Inte_H1(Lk, r1, r2, Gam_12, h, s, Eps)

      If r1 < Eps Then
        H1 = Lk / 4 / Pai * (Log(Lk) - 3 / 2)
      ElseIf r2 < Eps Then
        H1 = Lk / 4 / Pai * (Log(Lk) - 1 / 2)
      Else
        H1 = -1 / 4 / Pai / Lk * (r2^2 * Log(r2) - r1^2 * Log(r1) - 1 / 2 *((Lk - s)^2 - s^2))
        H1 = H1 + 1 / 2 / Pai * (1 - s / Lk) * (h * Gum_12 + (Lk - s) * Log(r2) + s * Log(r1) - Lk)
      End If

      B_Inte_H1 = H1
    End Function

    External Function B_Inte_H2(Lk, r1, r2, Gam_12, h, s, Eps)

      If r1 < Eps Then
        H2 = Lk / 4 / Pai * (Log(Lk) - 1 / 2)
      ElseIf r2 < Eps Then
        H2 = Lk / 4 / Pai * (Log(Lk) - 3 / 2)
      Else
        H2 = 1 / 4 / Pai / Lk * (r2^2 * Log(r2) - r1^2 * Log(r1) - 1 / 2 *((Lk - s)^2 - s^2))
        H2 = H2 + 1 / 2 / Pai * s / Lk * (h * Gum_12 + (Lk - s) * Log(r2) + s * Log(r1) - Lk)
      End If

      B_Inte_H2 = H2
    End Function

 以上の結果は、[計算部][外部手続き]のところへは、まだはめ込みません。[計算部]では、まだやる事があるからです(もぉ~、境界要素法って面倒くさいなぁ~(^^;))。


nice!(0)  コメント(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。