［境界要素法プログラムを設計する（境界積分ｻﾌﾞﾙｰﾃｨﾝ：具体調査編）］

［境界要素法プログラムを設計する（境界積分ｻﾌﾞﾙｰﾃｨﾝ：具体調査編）］

 

１．最後はボトムアップになる

 

 

　例は1個で十分です。

　　

 

　式(1)は、［内点方程式の積分公式（ラプラス型）］であげた、境界要素kの節点jのψ_jの係数を与える積分結果です（ずいぶん前のような気がしてますが(^^;)）。特異点が要素節点に一致しない場合、境界方程式でもこの形の計算は必ずします。というかこういう場合の方が、圧倒的に多いです。

　式中のL_k，s，γ₂，γ₁，h，r₂，r₁が図-1から定められる積分パラメータです。このうちL_k，s，h，r₂，r₁は単純に計算できます。例えばr₁なら明らかに、

 

　　

 

ですし、s，hはちょっと手がかかりますが、基本は(x0，y0)，(x1，y1)，(x2，y2)で係数が決まる2行2列の連立方程式を解くだけです。問題はγ₂，γ₁の与え方です。これは一見簡単そうに見えます。例えばγ₁なら、

 

　　

でいいんでないの？と。まず汎用言語の数学関数のAtn（tan^-1の事）は、値域がふつう－π/2～π/2です。なんか鬱陶しそうだから、x1－x0とy1－y0の符号に注意して場合分けすれば、容易に0～2πをカバーできます。これでγ2－γ1＞0は容易に計算できそうです。でも図-2のような場合はどうしますか？。

 

 

　Atnを0～2πをカバーできるように改造したとしても、図-2では明らかにγ2＜γ1でγ2－γ1＜0です。しかし線積分（境界積分）は常に左回りなので、図中の角度矢印で示したように（左回り）、γ2－γ1＞0でなければならず、しかも0～2πで表した|γ2－γ1|は本当の(γ2－γ1)と値が違います（2π大きい）。

　では境界要素がx軸方向を横切る時だけ、2π引いて符号を反転させればOKなのか？というと、一般的には図-2への対処も必要です。

 

 

　この場合は角度矢印が示すように（右回り）、γ2－γ1＜0が正解です。図-3では0～2πの角度定義でたまたまγ2＜γ1なので、そのままやって正解になりますが、図-2と3の複合ケースは？。・・・もう鬱陶しくてやってられません(^^;)。

 

　つまりγ1とγ2は（βkも）単独で計算しては駄目なんですよ。あくまでγ1とγ2の角度差分を局所情報に基づいて計算する必要があります。こういう事は、普通の積分では起こりません。普通は図-1のように一周したら値は元に戻って問題ないはずです。これは積分結果がtan^-1になるような正則でない関数を積分せざる得ない、特異積分を使っている境界要素法の特徴の反映です（複素積分を知ってる人なら、イメージできるはず）。そういう事情から、前回メインで大雑把に決めた関数の引数並びも、角度差分を与えるように変更する必要があります。

　これはボトムアップです。このように最後はボトムアップになる事は多いです。トップダウンだけでは問題の個別状況の全ては把握できないからです。でも最後までボトムアップは我慢した方が、たいてい上手く行きます。

 

 

２．角度差分を計算するアルゴリズムを考える

　一つだけ確かなのは、|γ2－γ1|はπ以下だ、という点です。|γ2－γ1|がπに成り得るのは、特異点が境界要素の中点上などにある場合なので、πを越える事はありません。そうすると格好の数学関数があります。Acos（cos^-1）です。Acosの値域は、汎用言語においてたいてい0～πです。内積の定義から、

 

　　

 

として、

　　

なので

　　

と求められます。|γ2－γ1|の符号は要するに、図-2，3の角度矢印の向きがわかれば良いので、ベクトル(x01，y01)と(x02，y02)の外積を利用する手があります。(x01，y01)はγ1方向，(x02，y02)はγ2方向のベクトルなので。(x01，y01)と(x02，y02)を3次元化して外積を取ると、

 

　　

 

になります。(7)右辺のz成分が＋なら左回りで、γ2－γ1＝|γ2－γ1|。(7)右辺のz成分が－なら右回りで、γ2－γ1＝－|γ2－γ1|と判定できます。ただし式(6)と(7)は小さすぎるr1やr2に対しては、結果が不安定になるのは目に見えています。本当にr1やr2が0なら、division by zeroエラーになります。なりますが、特異点が要素節点に一致する場合も、もともと扱うのでした。だとすればそういう特殊ケースは、特異点が要素節点に一致するケースに含めれば十分です。

 

３．特異点が要素節点に一致するケース

 

 

　ここでもγ2－γ1の扱いが問題になります。［内点方程式から境界方程式へ（ラプラス型）］でFree Termの項は、次のように計算されました。

 

　　

 

　式(8)上段の1/2π×()内の1項目が要素k+1のγ2－γ1すなわち、図-4のγ₃－γ^(j+1)です。2項目は要素kのγ2－γ1すなわち、図-4のγ^(j+1) －γ₁です。ただしε→0の極限で考えてます。

　でも待って下さい。ε→0の場合（特異点が要素節点に一致する場合）、どうやってγ^(j+1) を与えれば良いんでしょう？。だって数値計算上はε＝0なんですから。そしてすぐに気づくのは、特異点が要素節点に一致するケースでは、「個々の要素の境界積分値は、特異点の要素節点への近づき方に依存する！」という驚愕の事態が発覚します。もろな特異積分！。境界積分は結局、可積分でない？(^^;)。

　・・・とけっこう猛烈に焦っちゃう訳ですが、式(8)下段では不定なγ^(j+1)が引き算でキャンセルされて、「境界積分全体では近づき方に依存せず、値は一意に定まる」事になります。境界要素法は、こういう危ういところを持っています。ブレビアさんがなりふり構わず、コーシーの主値積分を持ち込んだ気持ちもわかりますよ(^^;)。

 

　つまり足せばγ^(j+1)は消えるのですから、逆に言えばなんでも良い訳です。例えば要素kで(x02，y02)＝(1，0)とし、要素k+1でも(x01，y01)＝(1，0)にするとか。共有節点でγの方法が揃ってればOKなはず。どうしてかと言うと特異積分値は、局所情報に基づいて計算する必要があるから。

　特異点が要素節点に一致するケースでは、その点でγ1，γ2は(1，0)方向に取ると決定します。後で後悔するかも知れませんが(^^;)。

 

　あと考えるべき事は、r1またはr2がどれくらい0に近かったら、0とみなすかです。これは【桁落ち誤差と丸め誤差はプログラマーの敵！】のところで書いたように、

　　http://nekodamashi-math.blog.so-net.ne.jp/archive/c2306081239-2

 

「問題のスケール規模」で決めるのが妥当と考えられます。r1，r2は特異点と要素節点の距離です。これに対応するスケール規模は明らかに、解析領域の大きさです。0判定Epsはどこで与えれば良いでしょう？。

　手がかりはたいてい「実行者」と「妥当な実行タイミング」にあります。「実行者」の第一候補は[入力部]です。[入力部]のRoll（役割）は、解析領域という外部情報をセットし[計算部]へ渡す事です。Epsは解析領域の大きさで決まるので、これは外部情報のセットです（そして渡す）。

　次に実行タイミングですが、Epsのセットは[計算部]の実行前に、事前に一回だけで十分です。[計算部]の実行前には、[入力部]の実行があります。「実行者」と「妥当な実行タイミング」が理想的に一致したので、[入力部]の中でEpsのセットを行います。セットの場所は、全ての入力が完了した後の方がなんか安心じゃないですか？。

 

　[メモ欄]と[宣言部]はできれば一対一に対応させたいので、[メモ欄]の最後に「REM Eps：0判定，Eps0：0ｵｰﾀﾞｰ」を追加し、[宣言部]の最後には「REM Eps：0判定」と「Eps0 = 10^(-6) REM 0ｵｰﾀﾞｰ」のコードを追加します。コードの心は、解析領域の大きさの1/100万より小さい距離は0とみなすです。そして[入力部]の最後に、

 

    REM 全ての入力完了
    REM 解析領域のx方向幅とy方向幅を取得
    x_Min = 10^10    REM x座標最小値初期化
    x_Max = -10^10   REM x座標最大値初期化
    y_Min = 10^10    REM y座標最小値初期化
    y_Max = -10^10   REM y座標最大値初期化

    For j = 1 to BN_Count
      h = BN(j, 1)  REM j節点のx座標

      If h < x_Min Then
        x_Min = h
      End If

      If x_Max < h Then
        x_Max = h
      End If


      h = BN(j, 2)  REM j節点のy座標

      If h < y_Min Then
        y_Min = h
      End If

      If y_Max < h Then
        y_Max = h
      End If

    Next j

    x_Width = x_Max - x_Min  REM 解析領域のx方向幅
    y_Width = y_Max - y_Min  REM 解析領域のy方向幅

    REM 0判定距離をセット
    If x_Width < y_Width Then
      Eps = y_Width * Eps0
    Else
      Eps = x_Width * Eps0
    End If

 

　このEpsの与え方は、解析領域の最大幅をスケール規模とするやり方です。これでは駄目な場合ももちろんあり得ますが、とりあえず(^^;)。今回は[メモ欄][宣言部][入力部]のみあげます。[計算部]は、引数並びから検討し直す必要があるとわかったので。

 

　ところでもう一個忘れてました。式(1)を見ると、定数πも必要です。これも[メモ欄]と[宣言部]に追加しておきます。

※ 後で調べてみたら、MAX関数とMIN関数と定数PIが組み込みになってました。まぁ～、好きにして下さい(^^;)。

REM [メモ欄]

REM B_Count：境界要素数
REM EN_B(k，2)：境界要素－節点対応表，k=1～B_Count
REM (k，1)：要素kの始端の節点No.，(k，2)：終端の節点No.

REM R_Count：領域要素数
REM EN_R(k，4)：領域要素－節点対応表，k=1～R_Count
REM (k，j)，j=1～4：jに従って要素kの要素節点No.を左回りに格納

REM BN_Count：境界節点数
REM BN(j，2)：境界節点座標，j=1～BN_Count，(j，1)：節点jのx座標，(j，2)：y座標
REM RN_Count：領域節点数
REM RN(j，2)：領域節点座標，j=1～RN_Count，(j，1)：節点jのx座標，(j，2)：y座標

REM BN_BP(j，2)：節点の境界条件の有無，j=1～BN_Count
REM (j，1)：節点jのψjが既知なら1、そうでないなら0
REM (j，2)：節点jのqjが既知なら1、そうでないなら0

REM BN_BV(j，2)：境界条件の値，j=1～BN_Count
REM (j，1)：節点jのψjに与える値
REM (j，2)：節点jのqjに与える値

REM z(i):未知ベクトル，i=1～2*BN_Count
REM BN_Z(j，2)：節点－未知数対応表，j=1～BN_Count
REM (j，1)：ψjのzの中の位置，(j，2)：qjのzの中の位置
REM 節点j→ψj→z(j)，節点j→qj→z(BN_Count+j)

REM A(i, j)：全体係数マトリックス、i,j=1～2*BN_Count

REM Eps：0判定
REM Eps0：0ｵｰﾀﾞｰ

REM Pai：π


REM [宣言部]
REM 配列寸法の上限は1000

REM B_Count：境界要素数
REM R_Count：領域要素数
Dim EN_B(1000，2)  REM 境界要素－節点対応表，k=1～B_Count
Dim EN_R(1000，4)  REM 領域要素－節点対応表，k=1～R_Count

REM BN_Count：境界節点数
REM RN_Count：領域節点数
Dim BN(1000，2)  REM 境界節点座標，j=1～BN_Count
Dim RN(1000，2)  REM 領域節点座標，j=1～RN_Count

Dim BN_BP(1000，2) REM 節点の境界条件の有無，j=1～BN_Count
Dim BN_BV(1000，2) REM 境界条件の値，j=1～BN_Count

Dim BN_Z(1000，2) REM 節点－未知数対応表，j=1～BN_Count

Dim z(1000)  REM 未知ベクトル，i=1～2*BN_Count
Dim A(1000, 1000)  REM 全体係数マトリックス、i,j=1～2*BN_Count

REM Eps：0判定
Eps0 = 10^(-6)  REM 0ｵｰﾀﾞｰ

Pai = Acos(-1)  REM π：Fortranになじんだ人なら、懐かしくて泣いちゃう式(^^)


REM [入力部]

REM境界要素数B_Countと領域要素数R_Countは、事前に与えられると仮定する。
REM 境界要素と領域要素データは、要素ごとに節点番号が与えられると仮定する。
REM 与えられる節点番号は、左回りに順序付けられたものが与えられると仮定する。

REM 境界要素を要素番号kで読み込むLoop
REM 　Loop内で要素ごとの節点番号を読み、境界要素－節点対応表EN_B(k，2)をセット

  For k = 1 to B_Count
    EN_B(k，1) = ？
    EN_B(k，2) = ？
  Next k

REM 領域要素を要素番号kで読み込むLoop
REM 　Loop内で要素ごとの節点番号を読み、領域要素－節点対応表EN_R(k，4)をセット

  For k = 1 to R_Count
    For j = 1 to 4
      EN_R(k，j) = ？
    Next j
  Next k

REM境界節点数BN_Countと領域節点数RN_Countは、事前に与えられると仮定する。

REM 境界節点を節点番号jで読み込むLoop
REM 　 Loop内で境界節点座標を読み、BN(j，2)をセット

  For j = 1 to BN_Count
    BN(j，1) = ？
    BN(j，2) = ？
  Next j

REM 領域節点を節点番号jで読み込むLoop
REM 　 Loop内で領域節点座標を読み、RN(j，2)をセット

  For j = 1 to RN_Count
    RN(j，1) = ？
    RN(j，2) = ？
  Next j


REM 境界条件の場所を節点番号jで読み込むLoop
REM 　Loop内で節点jの境界条件の有無を読み、BN_BP(j，2)をセット

  For j = 1 to BN_Count
    BN_BP(j，1) = ？
    BN_BP(j，2) = ？
  Next j

REM 境界条件の値を節点番号jで読み込むLoop
REM 　Loop内で節点jの境界条件の値を読み、BN_BV(j，2)をセット

  For j = 1 to BN_Count
    BN_BV(j，1) = ？
    BN_BV(j，2) = ？
  Next j

REM 節点番号jのLoopで、節点－未知数対応表BN_Z(j，2)をセット
REM (j，1)：ψjのzの中の位置，(j，2)：qjのzの中の位置
REM 対応は、節点j→ψj→z(j)，節点j→qj→z(BN_Count+j)

  For j = 1 to BN_Count
    BN_Z(j，1) = j
    BN_Z(j，2) = BN_Count + j
  Next j


REM 全ての入力完了 **********************************************
REM 解析領域のx方向幅とy方向幅を取得
  x_Min = 10^10    REM x座標最小値初期化
  x_Max = -10^10   REM x座標最大値初期化
  y_Min = 10^10    REM y座標最小値初期化
  y_Max = -10^10   REM y座標最大値初期化

  For j = 1 to BN_Count
    h = BN(j, 1)  REM j節点のx座標

    If h < x_Min Then
      x_Min = h
    End If

    If x_Max < h Then
      x_Max = h
    End If


    h = BN(j, 2)  REM j節点のy座標

    If h < y_Min Then
      y_Min = h
    End If

    If y_Max < h Then
      y_Max = h
    End If

  Next j

  x_Width = x_Max - x_Min  REM 解析領域のx方向幅
  y_Width = y_Max - y_Min  REM 解析領域のy方向幅

REM 0判定距離をセット
  If x_Width < y_Width Then
    Eps = y_Width * Eps0
  Else
    Eps = x_Width * Eps0
  End If

タグ：数値解析 境界要素法