So-net無料ブログ作成

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

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

 

1.最後はボトムアップになる

 

bem14-fig-001.png

 

 例は1個で十分です。

  

 

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

 式中のLksγ2γ1hr2r1が図-1から定められる積分パラメータです。このうちLkshr2r1は単純に計算できます。例えばr1なら明らかに、

 

  

 

ですし、shはちょっと手がかかりますが、基本は(x0y0)(x1y1)(x2y2)で係数が決まる22列の連立方程式を解くだけです。問題はγ2γ1の与え方です。これは一見簡単そうに見えます。例えばγ1なら、

 

  

でいいんでないの?と。まず汎用言語の数学関数のAtntan-1の事)は、値域がふつう-π/2π/2です。なんか鬱陶しそうだから、x1x0y1y0の符号に注意して場合分けすれば、容易に0をカバーできます。これでγ2γ10は容易に計算できそうです。でも図-2のような場合はどうしますか?。

 

bem14-fig-002.png

 

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

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

 

bem14-fig-003.png

 

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

 

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

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

 

 

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

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

 

  bem14-001.png

 

として、

  

なので

  

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

 

  bem14-002.png

 

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

 

3.特異点が要素節点に一致するケース

 

bem14-fig-004.png

 

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

 

  

 

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

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

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

 

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

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

 

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

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

 

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

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

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

 

 [メモ欄][宣言部]はできれば一対一に対応させたいので、[メモ欄]の最後に「REM Eps0判定,Eps00オーダー」を追加し、[宣言部]の最後には「REM Eps0判定」と「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

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

nice! 0

コメント 0

コメントを書く

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