［境界方程式に関する留意点（ラプラス型）］

［境界方程式に関する留意点（ラプラス型）］

 

　前回でやっと、境界方程式の具体的な離散化方程式にたどり着きました(^^;)。

  

　しかし式(1)を具体的に解くためには、なおいくつかの手順が必要です。

（境界要素法って、なんて面倒(^^;)）

 

［1］境界条件の考慮

　式(1)の係数行列B，Hはn×nの正方行列で、既知ベクトルwの次元はn，未知ベクトルψとqの次元もnです。よってこのままでは、条件数n＜未知数の数2n、となり不定解しか出せません。これは偏微分方程式一般に言える事で、偏微分方程式の解は境界条件を与えて初めて一意に決まります。逆に言えば、境界条件を与えない限り、式(1)は不定解以外は出していけないのです。これもある意味、数学ってなんて正直・・・(^^;)。

 

　ラプラス型の境界条件では、解析領域境界C上の全ての点で、ψ(c)かq(c)の値を与える必要があります。よって、式(1)の未知数2n個のうちn個の値が与えられる事になり、条件数n＝未知数の数nとなって解けるようになります。

　最も簡単に境界条件を考慮する方法は、境界条件を表す式を補助方程式として追加してやる事です。

  

　ここでj＝j₁，j₂，・・・とk＝k₁，k₂，・・・は、ψ_jとq_kの値が与えられる節点番号で、最も一般的には規則性はないと思った方が無難です。Ψ_jとQ_kは境界条件により与える値です。

　式(2)の左辺全体を行列D₁，D₂で表します。また対応する右辺をベクトルd₁，d₂で表します。補助方程式を追加すると式(1)は、

  

の形になります。ψとqを合わせた未知ベクトルの次元は2n。(D₁ D₂)はn×2n行列なので、n×2nの(B H)と合わせて、全体は2n×2nの正方な係数行列になります。右辺の既知ベクトルはもちろん、d₁とd₂を合わせた次元はn，wの次元もnなので、全体で2nです。

　(D₁ D₂)を(B H)の上に置くのには理由があります。(D₁ D₂)に適当な行の並べ替えをかけると対角行列になり、(3)のψ，qに関する掃き出し計算が高速化されるからです。この性質は行の並べ替えを行わなくても保たれ、普通にやればよろしい！という事になります(^^)。

 

［2］ポテンシャル値問題

　境界条件として最もすっきりするのは、境界上の全ての点でψ(c)の値が与えられるか、q(c)の値が与えられた場合です。これらのケースでは式(3)の形を使うまでもなく、

  

のどちらかを解けばOKです。

　ところが式(5)のBは特異行列になります。つまり不定解しか出ません。これの理由は物理的事情です。ラプラス（ポアソン）方程式の解関数ψ(x，y)は、必ず何かの物理量のポテンシャル関数になります。そしてポテンシャル自体は物理量に対応せず、その微分だけが意味を持ち、ポテンシャル関数には定数分の不定性があるからです。すなわち式(5)は、不定解を出さなければならないのです。・・・数学ってなんて正直 (^^;)。

 

　物理的理由がはっきりしているので、対処は簡単です。ポテンシャル関数の不定性を固定してやります。例えば節点1に対してψ₁＝0という補助方程式を追加してやり、節点1に関する境界方程式を削除します。節点1の境界方程式とは、

 

  

の事です。実際には、ケースバイケースに最も都合の良い条件を選びます(^^)。

 

［3］滑らかな境界点と本当の角点の区別

　前回の結果では、境界方程式(1)の作成にあたって解析的積分式を使用すれば自然に、

  

の自由項k/2/π×ψが出てくるのでした。理論的には滑らかな境界の場合はk/2/π＝1/2で、角点の場合は角の内角をk_iとしてk_i/2/πです。

　という事は滑らかな境界に折れ線近似を行った場合、解析的積分式からは、近似折れ線間の角度k_iが現れてしまう事になります。近似が十分なら、k_iはπに近い数値にはなりますが。

　そこで通常は、強引にk_iをπにおきかえます。近似値より正確な真の値を使うんだから精度はより上がるはずだ！、という方針ですが、本当かなぁ～？(^^;)。実際にはBマトリックスを作った後に、滑らかな境界点に対応する対角成分iに(1/2－k_i/2/π)を足し込んでやるのが実用的です。またこれは計算プログラムが、滑らかな境界点と角点を区別する情報を持たねばならない事をも意味します。

 

［4］2重流速問題

　これは本当の角点で発生し、角点問題とも言われます。[1]で述べたように、各境界節点iには未知数への条件として境界方程式、

  

が一つだけ存在します。ψ_iかq_iが境界条件で与えられるので、全体として(ψ_j)や(q_j)を決定できる訳です。

　ところが角点には2個のq_iがあります。何故なら角点iの前後で外法線方向が不連続にジャンプするからです。q(c)はψ(x，y)の境界上での外法線方向微分値でした。角点iで2方向の外法線方向β⁽¹⁾とβ⁽²⁾があるなら（図-1）、当然外法線微分値q_iもq_i⁽¹⁾，q_i⁽²⁾と二つある事になります。これは、Hマトリックスとqベクトルの構成を、多少面倒にします。同一節点について、2つのq_iをみなければならないからです。

 

　もし境界条件でq_i⁽¹⁾とq_i⁽²⁾の両方が与えられるなら普通に解けますが、ψ_iとq_iの一方だったり、ψ_iしか与えられない場合には、節点iで未知数2個に対し境界方程式一本は条件不足です。これに対する対処法は、1980年代には色々考えられました。

 

1) 定数要素を用いる方法

　今まで線形近似でやって来ましたが、もちろん定数近似だって可能です。境界要素上で未知量ψ(c)やq(c)は変化しない定数だ、という近似です。境界積分の式からはcの（tの)一次の項が消え、計算はより簡単になります。その際の節点配置は通常、境界要素の中央に置かれ図-2のようになります。

　図から明らかなように定数要素を使えば、2重流速問題は避けて通れます。積分計算が簡単な事もあり実用的にもまずまずで、境界要素法の入門本ではけっこう推奨されてますが、線形近似くらいしたいよねぇ～、ってのが個人的本音です(^^;)。

 

2)陰解法を併用する方法

　陰解法の利点は、基本解の特異点が積分領域外にあるので、未知数を増やす事なくいくらでも境界方程式を立てられる事です。だったら目標角点での直接法の境界方程式をやめて、陰解法の境界方程式を角点で必要な数だけ追加すれば良いじゃない、という考えです（図-4）。

　この方法を使えば、今までやってきた積分式をそのまま流用でき、そのうえ精度的にも直接法と遜色ないんですよ。実用的には、とってもお奨めです。

　ところがやっぱり外部特異点の最適配置と定量的誤差評価は？、と始まりまして、理論家達には不評でした。特異点配置が自由過ぎて、評価できなかっただけなんだろ？と思っちゃう訳ですが、そういう訳で自分はときどき好んで使います(^^)。

 

3)方向微分の公式を用いる方法

　たぶん現在の主流はこれです。かつ入門本では、あまり見た事ありません(^^;)。

 

　2重流速問題に対処するためには、未知数を増やす事なく補助方程式を追加できれば良い訳です。この意味で3)陰解法併用法はすじが通っています。

　境界要素法では、未知関数ψとその微分値qを独立な近似関数として扱います。そこが微分値と関数値の精度が同等になる境界要素法の利点の一つです。有限要素法などでは微分値の近似関数は、関数値の近似関数を微分した関数で与えるので、微分値の近似次数は関数値より一段落ちるのが普通です。qはψの微分値なので理論的には本当は、qとψの間に一般的な制約条件が存在します。その意味では有限要素法の方が妥当なんです。

　まず節点j+1における要素kと要素k+1でのψ(c)の接線方向微分値p^(k)とp^(k+1)を、ψ_j，ψ_j+1，ψ_j+2で近似します（図-4）。長さはそれぞれL_kとL_k+1です。

(9)

　　

 

　(8)，(9)の接線微分の方向は図-4より、β^(k)＋π/2，β^(k+1) ＋π/2になります。従って節点j+1におけるx方向，y方向微分値ψ_x，ψ_yとは方向微分の公式から、以下の関係があります。

　　

 

　一方、図-4のβ^(k)，β^(k+1)方向の法線微分値q_j+1⁽¹⁾とq_j+1⁽²⁾についても同様に、

　　

 

が成り立ちます。(10)，(11)から(ψ_x，ψ_y)を消去すると（補足参照）、

　　

が得られます。さらに(8)，(9)を考慮して、

　　

です。これが節点j+1近傍で、関数値と外法線微分値が満たすべき制約条件です。当然の事ながらβ^(k+1)＝β^(k)の時は（滑らかな境界点の時は）、条件(12)はq_j⁽¹⁾＝q_j⁽²⁾＝q_jを与えるので不要です。これの意味するところは、あまり扁平な角点は頑張らずに、滑らかな境界点とみなした方が安全よ(^^)、という事です。

　境界条件で節点j+1両側のq_j+1が与えられた時、および片側のq_j+1およびψ_j+1がわかる時は、普通に境界方程式を立てます。片側のq_j+1だけ、またはψ_j+1しか既知でない時は、境界方程式をやめて条件(12)を使います。

 

　以上で計算プログラムに突入する準備は整ったぁ～！、と思ったら、まだ領域積分項の扱いがありましたね(^^;)。もぉ～、境界要素法って面倒くさいんだから・・・。

 

（執筆：ddt³さん）

 

補足

［境界方程式の組み立て（ラプラス型）］

［境界方程式の組み立て（ラプラス型）］

　ここまでで境界方程式（内点方程式）を作成する計算公式は全て出揃ったので、いよいよ境界方程式を組み立てます。

 

　境界方程式と内点方程式は以下でした。

　　

 

　内点方程式(2)の左辺の境界積分を、図-1の一要素kだけ取り出して書いてやると、境界上でψ(c)とその外法線微分値q(c)を線形近似した場合、

　　

の形になります。

　ψ_jやq_jは、境界要素kの両端点にある境界節点jとj+1でのψ(c)とq(c)の値を表し、これらが未知数です。b_j^(k)やh_j^(k)は、境界要素kの配置と特異点(ξ，η)の位置だけから計算できるのでした。具体的な形は前回にあります。

 

　図-1に示したように、要素kの節点j+1は、隣の要素k+1と共有されるので、そこに注意して(2)の左辺を(3)，(4)を使って書くと、

　　

となります。いま境界要素はn個あり、境界要素と節点は左回りに順序付けられているとします。ψ₁とq₁の係数にb₁⁽ⁿ⁾とh₁⁽ⁿ⁾が現れるのは、要素nと要素1が節点1を共有するからです（図-1）。また境界要素と節点が同数あるのは、図-1から明らかです。式(5)で各ψ_j，q_jの係数を、B_j，H_jと書く事にします。

　　

です（j＝1の場合は、適当に変更して下さい(^^;)）。Σ記号を導入して、

　　

と書けます。一方(2)の右辺はψ^*もgも既知関数なので、なんとかすれば具体的値がわかるだろう、という事で(^^;)、たんにwと書きます。よって、

　　

　ところで内点方程式(2)の目的は、基本解の特異点(ξ，η)を任意に動かして、解析領域R内の任意の位置における未知関数ψの値を、ψ(ξ，η)の形で得る事でした。(ξ，η)の位置が変われば、B_j，H_j，wの値も当然変わります。そこで、(ξ，η)＝(ξ_i，η_i)に取った時の値をB_ij，H_ij，w_iと書く事にします。

　i＝1，2，・・・，[何個でも良い]です。

　そのような意味で式(8)は、

　　

と書くべきだぁ～という事になります。

 

　次に前回の結果から、特異点(ξ_i，η_i)を節点jに近付けて行けば、式(9)左辺のψ_jに関する和、

　　

から自然に、

　　

が導かれ、内点方程式(2)は連続的に境界方程式(1)に移行できるのでした。k_jは、節点jにおける境界の内角です。

　(ξ_i，η_i)→節点jの時のB_ij，H_ijの具体的形も前回の結果で与えられます。一般に式(11)に相当する項は境界要素法では、自由項（free Term）と呼ばれます（←あまり役に立たない蘊蓄(^^;)）。

　(10)，(11)が起きるのは、(ξ_i，η_i)が節点jに一致する時だけだという事に注意し、式(9)を境界方程式として書き直すと、

　　

になりますが、今度は(ξ_i，η_i)がどれかの節点jと一致するので、i＝1，2，・・・，nです。

　なおd_ijは、クロネッカーのデルタです。普通それはδ_ijと書かれますが、今までδはデルタ関数や変分の意味にさんざん使ってきたので、ここはdにしました(^^;)。

 

　式(12)をさらに整理するために突然ですが、行列とベクトルの積を思い出して下さい。行列A＝(a_ij)とベクトルx＝(x_i)との積は、

　　

って書きますよね？。これを成分で書くと、

　　

ですよね？。・・・式(12)と同じじゃないですか！(^^)。

　　

です。

　さらに、

　　

と以後略記します。(14)で行列はn×nの正方行列、ベクトルの次元は必ずnです。

　　

　形式的には式(15)が、境界節点で離散化して組み立てられた境界方程式の全てです。ψとqの係数行列BとHは既知であり、wも既知ベクトルです。後は(15)を、ψとqに関する連立一次方程式とみなして解けば良い訳ですが、なお留意点がいくつかあります。

 

［内点方程式から境界方程式へ（ラプラス型）］

［内点方程式から境界方程式へ（ラプラス型）］

 

　前回までで内点方程式の境界積分公式を与える事ができました。同様に境界方程式の境界積分公式も作れるはずです。境界方程式では、基本解の特異点(ξ，η)は境界節点j＝1，2，・・・のどれかと一致します。そこで、(ξ，η)が最初から節点jやj+1と一致した図-2を考え、最初から(ξ，η)を境界上において積分計算を実行すれば良い訳です。それは実行可能ですが、初めから積分計算をやり直す事にもなります(^^;)。

　それならばという事で、今までの結果を再利用するという手があります。つまり前回までの結果で(ξ，η)が節点jやj+1に近づく極限を試す、という方法です。これは上手く行く上に、境界要素法に関わるある種の曖昧さの一つを解消できます。

 

　図-2の(ξ，η)が節点j+1に近づく極限を考えます。これはr₂→0という事なのでr₂＝εと表し、またγ₂がεの方向である事を明示するためにγ^(j+1)と書く事にすると、容易に（図をにらむだけ(^^)）、

　　

を導けます。これらの条件を考慮すると、前回の積分公式の極限は、

　　

を使って、

　　　
　　

 　　

 　　

　同様に(ξ，η)→節点jでは、r₁＝ε，γ₁＝γ^(j)と表す事にして、

　　

なので、

　　

を得ます。

　これらを、［境界積分の積分部品（ラプラス型）］の最後のボトムアップ公式、

　　

 に代入すれば要素k上で、境界方程式の境界積分は、

　(ξ，η)→節点j+1の時：

　　

　(ξ，η)→節点jの時：

　　

となります。

 

　注目は次です。以上の計算は、基本解の特異点を図-1の解析領域Rの内点として、(ξ，η)→節点j，j+1とした状況なので、内点方程式、

　　

の(ξ，η)→節点j，j+1の極限を取った状況です。式(17)の境界積分では要素kだけでなく、全ての要素の積分値を足すので、要素k上の積分値には必ずその隣の要素k+1での値も足されます。

　式(17)左辺において、節点j+1を共有する要素kと要素k+1に注目します（図-3）。図-3で(ξ，η)が節点j+1に近づいた時のψ_j+1の項をまとめます。

 

　式(12)，(15)と図-3の記号を使って、

　　

 

　ここにβ^(k+1)とβ^(k)は、要素k+1と要素kの外法線方向である事の明示です。

　(ξ，η)→節点j+1でした。従って、ψ(ξ，η)→ψ_j+1です。

　　

になります。

　ところが図-4に示すように、要素k+1は節点j+1からj+2にいたる線分で、その方向はβ^(k+1)＋π/2，要素kは節点jからj+1にいたる線分で方向はβ^(k)＋π/2。よって(β^(k+1)－β^(k))は、要素k+1と要素kの方向差です。図-4からπ－(β^(k+1)－β^(k))は、節点j+1の内角k_j+1となり、式(18)が得られます。

　この意味するところは、内点方程式(17)で(ξ，η)を境界に近づけたら、内点方程式は境界方程式、

　　

に化けたという事です。

　境界要素法の普通の定式化では、境界上でのデルタ関数の再評価で式(18)を導くので、(18)は(17)と無関係に不連続に成り立つように見えますが、ちゃんと積分してやれば、内点方程式から境界方程式へちゃんと連続的に移行できるのがわかります。数学ってなんて正直！(^^)。

　もっとも内点方程式から境界方程式への移行には式(1)が必要で、この計算はコンピュータには無理ですから、計算プログラム上は人間が式を使い分けてやる必要はあります(^^;)。

 

　これで計算準備はほとんど整いました。残ったのは式(17)，(18)右辺の領域積分項ですが、これは比較的簡単に扱えるので後回しにします。次にやる事は、数値計算用の境界方程式を組立てです。

 

※

　ここまでの内容を日本で最初に論文にしたのは、静電場解析にラプラス型の境界要素法を応用した、北大電気工学科の先生です（1980年代）。先生の名前は忘れましたが、2次元境界要素法，解析的積分あたりでググれば（かつ幸運なら）、その論文なんかにヒットするかも知れません。この先生はさらに、3次元領域表面を3角形メッシュで離散化した場合の、ラプラス型境界要素法の境界積分の解析的積分公式も与えています。いずれにしろ、「岩波数学公式集ゅう～～！」・・・です(^^)。

 

［内点方程式の積分公式（ラプラス型）］

［内点方程式の積分公式（ラプラス型）］

 

 

 

　今までの計算では、基本解の特異点(ξ，η)は図-1の解析領域Rの内点であるとして全てやってきました。すなわち今までの結果は全て、内点方程式、

　　

の境界積分に対応するものです。

　前回の結果から、前々回（［境界積分の積分部品（ラプラス型）］）のボトムアップ公式の内、式(27)～(34)は既に片付いています。

　　

　　

　ここではボトムアップ公式の残りも再掲し、一つ一つ計算して行きます。

　　

　　

　　

　　

 

　　

 

　単純に代入してくだけですよ(^^)。

　　

 　　

 

　そういう訳で・・・、

　　

　　

　という訳で後は、境界要素ごとに式(14)，(15)を計算し、式(1)に従い足し込んで行けばψ(ξ，η)の値が得られます。ここで積分パラメータs，h，r₁，r₂，γ₁，γ₂は、基本解の特異点位置(ξ，η)と境界要素kの配置および長さL_kから事前に計算できるのでした。また内点方程式(1)を使用する段階では、境界未知量ψ_jとq_jは全て求まった後と想定して良いのでした(^^)。

［境界積分の積分公式（ラプラス型）］

［境界積分の積分公式（ラプラス型）］

　・・・という訳で、「岩波数学公式集」です(^^)。一つの境界要素k上でのラプラス型の境界積分は、次の形の積分計算が出来れば良いのでした。t＝c－sとして（cとsは図-2参照）、

　　

ここでn≧0は整数です。

　(1)に部分積分を使えば、

　　

になるので結局、式(2)のu_n+2(t)がわかればOKです。

 

　そこで突然ですが、形式的に次の等比数列の和、

　　

 

を考えます。

　pは何でも良いので、p＝－t²/h²とすれば、

　　

となり、順次変形して行けば、

　　

が得られます。

　両辺にtをかければ、

　　

です。

　n＝0，1，・・・だったので、

　　2(n＋1)　 ＝2，4，6，・・・

　　2(n＋1)＋1＝3，5，7，・・・

となり、(4)，(5)を式(2)に代入した姿を想像すれば、u₀(t)とu₁(t)さえ計算できればOKとわかります。

　　

なので、

　　

です。

 

　これらを式(3)に代入して、

　　

　　

を得ます。

 

　v₀(t)とv₁(t)については、式(3)，(8)，(9)から計算すれば、

　　

となります。

 

　ところで境界上の未知関数ψ^*(c)，q^*(c)を線形近似した場合、v₀(t)，v₁(t)とu₀(t)，u₁(t)があれば良いのでした(^^)。
これらは既に計算してあるので、後はこれらを、c＝0→L_kについて定積分化するだけです。

 

　t＝c－sでした。またhとsは、境界要素kの配置と特異点(ξ，η)の位置だけで決まり積分に対して定数なので、積極的理由のない限り出来るだけそのままの形で使用した方が便利です。そうすると図-2より、例えば、

　　

が得られます。ここにγ₁，γ₂は、図-2に示したc＝0とc＝L_kの時のr＝r₁，r₂の方向（角度）です。同様に、

　　

　　

が得られます。

（執筆　ddt³さん）

これは、ネムネコのひとりごとです。

　　

の積分は、次のように置換積分を使ったほうがよいでしょう。

　　

したがって、

　　

 また、

　　

 

［境界積分の積分部品（ラプラス型）］

［境界積分の積分部品（ラプラス型）］

前回の最後で、境界積分を解析的に実行すると決心したのでした(^^;)。そこで境界積分に必要な積分部品をトップダウンで特定し、一つ一つ解析的に積分してボトムアップする事にします。そういう訳でここでは、記述を系統立てる記号の定義に終始します。

　図-1に示した一つの境界要素kでの境界積分は、境界上で解関数ψとその外法線微分qを線形近似した場合、

　　

となりました。L_kは要素kの長さ、ψ_jやq_jは図-2に示した節点jとj+1でのψ(c)とq(c)の値です。

　(1)，(2)の境界未知数ψ_j+1，ψ_j，q_j+1，q_jの係数を、要素k上の未知量ψ_jやq_jに関する係数という意味で、

　　

　　

と書きます。(3)～(6)を使うと(1)，(2)は、

　　

と書けます。図-2に示した幾何学的な積分パラメータによって、b_j^(k)やh_j^(k)の具体的形を与える事が、当面の目標です。

　これも前回の結果から、

　　

です。(9)，(10)を(3)～(6)に代入します。

　　
　　
　式(11)～(14)を眺めると、

　　

 　　

という定積分を求積出来れば良いとわかります。ところで積分計算の一般的傾向として、「分母は出来るだけ簡単に、logの中身も出来るだけ簡単に」した方が、「計算としてまだマシ！」ってのがありますよね？。

　よって上記は、c－s＝tと置換して、

　　

　　

とやるのが安全です。式(19)～(22)の積分部品を拾うと、

　　

　　

　さらに定積分は不定積分がわかれば良いので、式(23)～(26)の不定積分を(t)付きで表し、以上を全部のまとめて、ボトムアップ公式として書き出します。

　　

　　

　　

　　

　　

　　

　　

　　

　大変そうに見えますが、式(27)～(30)が計算できたとすれば、その結果を、式(31)～(34) → 式(35)～(38) → 式(39)～(42) → 式(43)，(44)と順番に代入して行けば良いだけです。こういう事はコンピュータの最も得意とするところです。
 　よってあと人間のやるべき事は、式(27)～(30)の不定積分を決定する事だけですよ(^^)。

（執筆 ddt³さん）

［ラプラス型境界要素法の境界積分の詳細］

［ラプラス型境界要素法の境界積分の詳細］

　ラプラス型境界要素法の基礎式は以下でした。

 

　内点方程式：(ξ，η)はRの内点。

　　

　境界方程式：(ξ，η)はRの境界C上。

　　

　ここで求めたい未知関数ψ(x，y)はポアソン方程式、

　　

を満たします。g(x，y)は既知関数です。(ξ，η)を特異点と呼びます。

 

　図-1に示した領域Rが、解きたい偏微分方程式(3)の解析領域でCはRの境界とします。式(1)，(2)の∫_C dcは境界C上の線積分，∫_R dxdyはRでの領域積分です。式(2)のkはk(ξ，η)で、その点での境界Cの内角を表します。

　ψ^*はラプラス方程式の基本解で、

 

　　

を取れます。qとq^*は、図-1に示した境界C上での外法線の方向βへの、ψとψ^*の外法線微分値になります。また式(4)のrの方向をγとします。境界上のψとqが具体的に定めたい未知数になります。何故ならそれらが求まると式(1)から、領域内任意点のψの値を計算できるからです。

 

　最初に式(1)の境界積分項を検討します。右辺の領域積分は、比較的簡単に処理できます。

　ψとqを具体的に近似するために、領域Rを折れ線近似し、折れ線の各線分を境界要素と呼んで要素番号kを与えます（図-1）。線積分の定義に従い、kの番号付けはCを左回りに一周が便利です（k＝1，2，・・・n）。

　図-1の境界要素を一個取り出して、外法線方向βがちょうど上向きになるように描いたのが、図-2です。

　∫_C dcは境界積分なので、図-2の境界要素上の点(x，y)が積分点です。線積分の定義より積分点(x，y)は、右から左へ走る事になります。それに伴って積分パラメータcは、要素kの長さをL_kとすれば、c＝0→L_kと増加し、これが要素k上での∫_C dcの積分区間になります。

 

　要素kの長さL_kが十分小さければ、要素k上でψ(c)とq(c)の変化は線形近似くらいで十分です。ψ(c)とq(c)を線形近似するために境界要素の両端に節点を配置し、節点番号jを与えます。線積分の定義に従い、jの番号付けもCを左回りに一周させるのが、一番便利です。明らかにj＝1，2，・・・nで、境界要素数と同数の節点が配置されます。節点jとj+1でのψ(c)およびq(c)の値を、(ψ_j，ψ_j+1)および(q_j，q_j+1)で表します。

 

　線形近似だったので要素k上では、

　　

です。これを(1)の境界積分項に代入すると、要素k上で、

　　　

 となります。

　かなりごちゃごちゃしてますが（数値計算なんてそんなもの(^^;)）、式(1)，(2)が基礎式である限り、具体的に積分できなければお話にもなりません。そういう目で(8)，(9)を見直すと、基本解に関するψ^*(c)とq^*(c)の具体的形が必要なのがわかります。

 

　まずψ^*(c)については式(4)よりrのみの関数なので、rをcで表せば良い事になります。特異点(ξ，η)から境界要素kへ下した垂線の長さをh，垂線の足の位置を要素局所座標系cで表した値をsとします。c＝0の時とc＝L_kの時のrをr₁，r₂として、rの方向をγで表すのにならって、r₁とr₂の方向をγ₁，γ₂とします。

 

　図-2から、

　　

ですが、r₁，β，γ₁は境界要素kの配置と特異点の位置だけで決まるので、これらはcに対して定数です。

　従って、

　　

に関する積分をすれば良い事になります。

　次にq^*(c)については、ψ^*(x，y，ξ，η)の境界上での外法線方向微分値なので、最初に、(x，y)におけるψ^*(x，y，ξ，η)のx方向偏微分とy方向偏微分を計算します。

　　

ここでrは式(4)で表したものであり、cosθ＝(x－ξ)/r，sinθ＝(y－η)/rです。(x，y)が境界上にある場合は、rは式(12)で表され、図-2から明らかにθ＝γになりますが、

　　

も図-2から得られます。従って外法線β方向への方向微分の公式に(14)～(17)を使うと、

　　

が得られます。

 

　よって式(8)，(9)の積分は、

　　

の積分計算に帰着できるのがわかります。

　(21)と(22)は、有理関数の積分になってます。なので、「やりゃ～必ず求積できる」事にはなります。(19)と(20)は、logの前にかかった1とcを先に積分する部分積分を実行すれば、やはり有理関数の積分に帰着でき、「やりゃ～必ず求積できる」事にはなりますが・・・。あえてやりたい計算ではないですよね？(^^;)。

 

　しかしやるんです。やるべきなんです！。どうしてかと言うと、大学以上の数学において「計算すりゃぁ～答えが得られる」なんていう事態は、ほとんど万に一つの幸運だからです！。もっと賢いカッコ良い方法は、計算結果を得た後でいくらでも考え付けます。そういう事を出来るようになるためにも、地道に計算するんですよ。それが大人の数学だと思います。

 

　ただし「岩波公式集」なんかは使って良いよ、という事にもなります。これも大学以上の数学の大人の対応です(^^)。

 （執筆　ddt³さん）

 

境界要素法入門５

境界要素法入門５

 

［境界要素法-3］

　境界要素法の内点方程式。

　　

　は解析領域R上の積分，はRの境界C上の線積分，φはΔφ＝g(x，y)を満たす未知関数です。qとはφとの外法線微分値。はδをデルタ関数として、

を満たし、

　　

を取れます。

　以上が出発点です。

 

　φ(ξ，η)は、

　　

の積分結果なのでした。(ξ，η)を積分領域Rの内部に置いた場合（内点とした場合）、(1)が得られますが、解けない形でした。(ξ，η)を積分領域Rの外部に置いた場合（外点とした場合）は、

　　

となり離散化すれば解けますが、不評なのでした。とすればもう後できる事は、(ξ，η)をRの境界Cに置く（境界点とする）しかないじゃないですか(^^;)。

　デルタ関数δ(ξ，η)を境界上で考えRで積分する場合、の再評価が必要になります。評価結果は簡単ですが、ここでたいてい怪しげな計算が登場し（コーシーの主値積分）、一部の理論家からは酷評される事になります。デルタ関数の数学的に厳密な定義にまで遡って再評価するなんて、やってられないですもんね(^^;)。そこでここでは、どうせいい加減になるならばという事で、の再評価にデルタ関数の等方性を使います。

 

　デルタ関数δ(ξ，η)って、特異点(ξ，η)を中心に等方的ですよね？（本当は関数じゃないけど）。だって(x，y)≠(ξ，η)ではδ＝0で、(x，y)＝(ξ，η)ではδ＝∞なんですから、これを等方的と言わずして何と言う？です(^^)。

　次に(x，y)≠(ξ，η)ではδ＝0ですから、(ξ，η)を内点として含む限りどんな積分領域を取っても、積分結果は同じです。なので特に積分領域として、(ξ，η)を中心とした半径εの円を取れます。(ξ，η)を中心とした円は、(ξ，η)を中心に等方的です。

　等方的な関数を等方的な領域で積分したら、どうなりますか？。例えばを、(ξ，η)を中心とした円で積分したら。もしになったとしたら、半円での積分値は1/2ですよね？。1/4の扇型に積分領域を制限したら、明らかに積分値は1/4ですよね？。

 

　δ(ξ，η)を境界上に置いた場合、半径εの円が十分小さければ(ξ，η)の近傍で、積分領域Rの境界Cはふつう直線とみなせます。半径εの円は、その直線によって半分だけRに引っかかります。よって、

　　

です。εがいくら小さくても直線とみなせないケースもあります。(ξ，η)がCの角点になるケースです。この時は角の内角をkとすれば明らかに、

　　

です。従って、

　　

です。デルタ関数はこういう事が、普通の関数と同じく実用的に出来るように、非常に注意深く造られたものです(^^)。

 

　前回やったように(1)を離散化します。結果は、

　　

でした。ここでは具体的な数値で与えられます。iは特異点番号，jは節点番号です。

　いま特異点は境界上にあるので、iはどれかのjと一致します。一致するj＝iの内角をとすると、(4)左辺のは、におきかわる事になります。i＝1，2，・・・を考慮して行列記法で書けば、

　　

です。iが角でない場合なので

　　

はクロネッカーのデルタ、すなわち、

　　

です。

　(5)の未知数も境界未知量のみなので、これも境界方程式と呼ばれます。(5)を解いて境界上のφ，qを定め(1)を併用するのが、境界要素法の直接法です。

 

　(5)から明らかなように、直接法では自動的に未知数の数に等しい条件数が得られます。また(5)の理論的誤差は、境界Cを折れ線近似した境界要素の長さと近似関数で決まるので、誤差評価も容易になります。

　こうしてユーザーにはよりわかりやすく（← ホントか？(^^;)）、理論家達も満足する定式化が得られました。

 

　しかし代償もあります。基本解

　　

はr＝0に特異性を持ちます。間接法の場合、の特異点(ξ，η)は積分領域Rの外部にあるので、(1)右辺の積分は普通に行えます。直接法の場合(ξ，η)は境界上にあるので、r＝0を含む特異積分を処理する必要に迫られます。

 

　の特異性はlog(r)のオーダーなので可積です。の特異性はアークタンジェントの計算に帰着できます。ただしそのような計算はコンピューターには出来ないので、特異点を持つ境界要素の積分では、人間の指示で場合分けし、解析的積分公式をプログラムに書いてやるか、特殊要素を用いる事になります。

 

　でもやれば、何とかなりますよ(^^)。

 

境界要素法入門４

境界要素法入門４

 

［境界要素法-2］

　境界要素法の内点方程式。

　　

　は解析領域R上の積分，はRの境界C上の線積分，φは

　　

Δφ＝g(x，y)を満たす未知関数です。qとはφとの外法線微分値。

　はδをデルタ関数として、

　　

を満たすものなら、何でもOKでした。そこで(2)を満たすとして出来るだけ簡単なものを考えます。これは実際に計算可能になる事が多いです。例えばとして(ξ，η)を中心に等方的なものを選び、境界条件は考えないとすれば、

　　

が容易に得られます。

　　http://nekodamashi-math.blog.so-net.ne.jp/archive/c2305704962-1　　　(a)

　　の［コーシーの積分公式の周辺-1］

　(1)の入出力関係を整理します。は具体的に決められたので、右辺で *つきのものは全て既知関数です。領域積分項はg(x，y)も既知関数なので、具体的に計算可能です。未知なのは、境界積分項に現れるφとqだけです。そうすると境界C上のφとqさえわかれば、(1)から未知関数φ(x，y)を計算できる事になります。

 

　解析領域Rの境界Cを折れ線近似し、折れ線の角に節点jを配置した絵を想像して下さい。jは節点番号で、左回りにCを一周します。節点jとj＋1を端点として持つ線分を境界要素と呼び、それに要素番号kを与えます（左回りにCを一周）。節点jでのφとqの値は、とで表します。

　境界要素kの要素長が十分に小さければ、k上の関数値：φ(c)，q(c)，は、(と)および(と)で線形近似でもすれば、十分なはずです。，0≦t≦1として、

　　

 

　一方、特異点によるとを、で表します。iは特異点の場所を識別する番号で、i＝1，2，・・・で十分です。

　これらを(1)の境界積分項に代入し、要素k上で考えてやれば、

　　

となります。

　上記右辺の積分において、とは既知関数なので、要素k上で具体的にの形に書けます。tは積分パラメータ（0≦t≦1)です。

　従って右辺の積分は、具体的な数値とになります。に関するについての積分を，に関するについての積分をで表しました。

　境界積分全体の値は、これらの集計です。節点番号jについて和をとる事になります。

　　

ただしは、

　　

は、

　　

です。

 

　ここで行列記法を思い出して下さい。i＝1，2，・・・だから、

　　かつi＝1，2，・・・は、と書ける。

　　かつi＝1，2，・・・は、と書ける。

・・・ですよね（← OKですか？(^^;)）。

 

　　

となります。

　(3)では、具体的に計算可能な既知量でした。よって(3)は、未知量とに関する連立一次方程式の形に、ほぼなっています。係数行列はとです。もし(3)を解いてとを定められれば、線形近似の形で境界上のφとqが得られた事になり、内点方程式(1)から未知関数φ(x，y)の近似解を求められます。左辺のさえ何とかなれば・・・。

　は、

　　

なのでした。でもこれは、「が積分領域Rの内点の場合には」なのでした。つまりは、Rの外にあっても（外点でも）OKです。その時は、デルタ関数の性質から、

　　

です(^^)。そういう訳で、

 

　　

によって、境界未知量を計算できる事になり、(1)を併用します。

　実際には境界条件として与えられた既知量を(4)に代入し、式の数を調整します。典型的には、境界上のφの値全部かqの値全部が与えられます。理屈の上では特異点は、必要な数だけR外部の任意の位置に、自由に設定すれば良い事になります。

 

　(4)には境界未知量しか出てきません。そこで(4)を境界方程式と呼びます。(4)の理論的形はグリーン関数法のいわば特殊解法として、昔から知られていたものでもありましたが、FEMになれた人達が気づいたのは、次の点でした。

 

　においてとすると、は（も）線形独立な関数系を構成します。従って、i＝1，2，・・・によって与えられた(4)は、独立な線形条件を成すはずだと。

　またいかに連続関数であろうと、離散化して考えれば節点自由度の有限の自由度しか持たないと。だとすれば、変分は完全に任意である必要はなく、i＝1，2，・・・で十分に任意化されている、と。

 

　このような発想が手軽に出来るようになったのは、やはり関数空間以後ですが、それだけでは理論止まりだったと思われます。「要素に区切ってやっつける」計算スタイルがFEMにより具体的に普及した事が、最後の決定打になった気がします。

 

　以上のやり方は理論的にも具体的な計算方法としても、最も明解なものです。ところが境界要素法の中で(4)は間接法と言われ、じつはマイナーなんです。というのは特異点の外部配置が余りにも自由過ぎて、数値計算の専門家達にも正確な誤差評価が無理だったからです。

　自分は、実用的な数値解が得られりゃOKよ程度のユーザーなので、(4)も時々使いますが、大概は専門家達に不評な方法なら、趨勢は右にならえです。

 

　という訳で、境界要素法の直接法という話になります。

 

（執筆：ddt³）

境界要素法入門３

境界要素法入門３

 

［境界要素法-1］

　線形問題に限って言えば、境界要素法は恐らく最強です。精度に優れ計算量も少なく、任意点の結果を取りやすい。計算プログラムの構成も明解でシンプル。

　難点は特異積分（広義積分）を扱う必要のある点ですが、今では解析的積分公式やさまざまな計算テクニックも開発されてるので、やりゃ～何とかなります。

 

　境界要素法（Boundary Element Method, BEM）の原理は有限要素法という計算スタイルより先行し、グリーン関数法です。

　　http://nekodamashi-math.blog.so-net.ne.jp/archive/c2305704962-1　　　(a)

　　の［コーシーの積分公式の周辺-1］

 

　グリーン関数法は条件さえ合えば確かに便利な手段でしたが、簡単な境界条件でしか便利にならないという計算法としては致命的な所がありました。それでいかなる境界条件にも理論上は簡単に適用可能な、変分原理に基づいた有限要素法（Finite Element Method, FEM）が先に発展します。そしてFEMがそろそろ常識化しかけた頃、という事は「解析領域を要素に区切って数値計算」という計算スタイルに大抵の人が慣れてきた頃、一部のFEM人がラプラス方程式などをそういう目で見直し、それらも「要素に区切ってやっつけられる」事に気づきます（恐らくブレビア（※）あたり）。1970～80年代の出来事です。で事後検証してみると、「要素に区切ってやっつける定式化」は、グリーン関数法の現代風アレンジになっていたぁ～、という次第です(^^)。

 

（※）　ブレビア　境界要素法入門（培風館）

https://goo.gl/ZwtTaK

 

　変分原理から導けるのはラプラス方程式だけではありません。ポアソン方程式も導けます。

　　

としてJを最小化すれば（やり方はこの前といっしょ）、ポアソン方程式、

　　

が得られます。

　(2)は湧き出し密度分布g(x，y)を持つ、非圧縮性渦なし完全流体の支配方程式で、g(x，y)は既知関数です。Jの意味は、流体の運動エネルギーと湧き出しのポテンシャルエネルギーも考慮した、流体全体のエネルギーを表します。

　φが静電ポテンシャルなら、g(x，y)は与えられた電荷密度分布。Jは電荷のエネルギーも考慮した静電場エネルギーです。

　グリーンの定理に戻ります。グリーンの定理はグリーンさんが特にグリーン関数法のために、ストークスさん供に導いたものです。

　　

　(3)を＝0とおけばラプラス方程式（やポアソン方程式）が得られるのでした。従って(3)とともにラプラスやポアソン方程式を与えれば、変分を取ったのと同等になります（本当は要証明ですが(^^;)）。

　　

 

　ところで変分δφは実質、連続関数であれば何でもOKでした。それでδφとして、φと同じタイプの方程式を満たすものを考えます。問題としては、偏微分方程式(2)を解きたいとします。

　変分δφは、(2)のタイプとして特に次の条件を満たすとします。

　　

　δ(ξ，η)は点(ξ，η)に特異点を持つデルタ関数で、をラプラス方程式の基本解と呼ぶ事が多いです。ここでが(x，y)だけでなく(ξ，η)の関数にもなってるのは、(4)からはδ(ξ，η)で決定されるので、はきっとその特異点(ξ，η)もパラメータとして含むよね、という意味です。

 

　形式的にを(3)に代入します。

　　

　で、はポアソン方程式を満たすのでした。という事は今度はに対してφを、その変分だと考える事も可能です。次式は直接計算しても出せます。

　　

は明らかです（・は内積）。

(5)と(6)の辺々引くと、

　　

を得ます。

　ここでφは(2)を満たし、は(4)を満たします。

　いま(ξ，η)が積分領域Rの内点なら、デルタ関数の性質から、

　　

です。従って先の関係式から、

　　

という事になります。

　ここで、もう少し表現を簡単にします。

は、Rの境界Cの線素をdcとした時、方向微分の公式から

　　

になるのがわかります。∂φ/∂nは、C上でのφの外法線方向微分値です。これをqで表しますが、についても同様です。は、境界上の流速値とか外法線微分とか呼ばれる事が多いです。

　　

　(ξ，η)は解析領域Rの内点であれば、どこでもOKでした。よってφ(ξ，η)は、Rの任意の内点(ξ，η)における解関数φ(x，y)の値です。それで(7)を、境界要素法の内点方程式と呼びます。

　(7)の関係式自体は、グリーン関数法の開発過程で相反定理の名のもとに知られていました。当時はデルタ関数の方法がなかったので、こんなに簡単には導けませんでしたが。

　じっさい基本解がφと同じ境界条件を満たすとすれば、これはRの境界C上でという事なので、(7)のC上の境界積分項は打ち消しあって0となり、参考URL(a)にあげたオリジナルなグリーン関数法の計算式と同等なものになります（グリーン関数の対称性）。

　しかし(4)を満たし、その上φと同じ境界条件を満たすをみつける計算の方が、(2)を解くよりよっぽど難しいという事態はよく起こります。そこがグリーン関数法の実用上の問題でした。

 

　(7)は、変分関係式(3)に支配方程式(2)を持ちこんだものです。もし変分δφが任意であれば、それは変分を取ったのと同等であり、(7)から陰的に解関数φを計算できる可能性が開けます。しかし(7)では、(4)を満たすという特定の変分しか用いていません。ところが条件(4)があればは十分に任意化されていると、FEMになれた一部の計算技術者や理論家達は気づいたのでした。

（執筆：ddt³）