境界要素法

境界要素法（きょうかいようそほう、英: boundary element method、BEM）とは、汎用性の高い離散化解析手法の1つで、有限差分法、有限体積法、有限要素法と並び、汎用離散化解析手法の主要3解法の1つとして理工学の分野で受け入れられている。電子計算機の発明・発展以前から進められてきた、応用数学における積分方程式の研究に端を発していることもあり、境界積分方程式法（Boundary Integral Equation Method、略してBIEM）と呼ばれることもある。 電磁気学では、この境界要素法をもちいた電磁界解析をモーメント法（Method of Moments、MOM）と呼んでいる。 解析手法は、積分方程式の定式化と離散化の2段階を経て構成される。 境界要素法では、まず対象とする問題の支配（微分）方程式から境界積分方程式を導出する。定式化には直接法と間接法の2種類がある。今日では、支配方程式の未知量をそのまま積分方程式の未知量として取り扱うことのできる、直接法定式化を採用する場合が多い。ここでは、2次元ラプラス問題を例に、直接法定式化による境界積分方程式の導出方法を説明する。 ラプラス問題は、支配方程式： ∇ 2 u = ∂ 2 u ∂ x 1 2 ( x ) + ∂ 2 u ∂ x 2 2 ( x ) = 0 , ( x ∈ Ω ) {\displaystyle \nabla ^{2}u={\frac {\partial ^{2}u}{\partial x_{1}^{2}}}({\boldsymbol {x}})+{\frac {\partial ^{2}u}{\partial x_{2}^{2}}}({\boldsymbol {x}})=0,\quad \quad ({\boldsymbol {x}}\in \Omega )} と、境界条件： u ( x ) = u ¯ ( x ) , ( x   o n   Γ u ) , q ( x ) = ∂ u ∂ n ( x ) = q ¯ ( x ) , ( x   o n   Γ q ) , {\displaystyle {\begin{aligned}u({\boldsymbol {x}})&={\bar {u}}({\boldsymbol {x}}),\quad ({\boldsymbol {x}}\ \mathrm {on} \ \Gamma _{u}),\\q({\boldsymbol {x}})&={\frac {\partial u}{\partial n}}({\boldsymbol {x}})={\bar {q}}({\boldsymbol {x}}),\quad ({\boldsymbol {x}}\ \mathrm {on} \ \Gamma _{q}),\end{aligned}}} とを同時に満たす解（ポテンシャル）u を求める問題である。ここで、Ωは領域であり、領域の境界Γは、ポテンシャルu が規定されている境界Γu と、フラックス q = ∂ u / ∂ n {\displaystyle q=\partial u/\partial n} が規定されている境界Γq からなり、 Γ u ∪ Γ q = Γ , Γ u ∩ Γ q = ∅ {\displaystyle \Gamma _{u}\cup \Gamma _{q}=\Gamma ,\Gamma _{u}\cap \Gamma _{q}=\emptyset } であるものとする。また、n は境界での外向き法線方向を示す。 上で示した支配方程式と関数u* とをかけ合わせてΩに関する領域積分を考えると、u が真の解であれば支配方程式を満足するため、これを含む項を積分しても0となる。 ∫ Ω ( ∂ 2 u ∂ x 1 2 ( x ) + ∂ 2 u ∂ x 2 2 ( x ) ) u ∗ ( x ) d Ω = 0 {\displaystyle \int _{\Omega }\left({\frac {\partial ^{2}u}{\partial x_{1}^{2}}}({\boldsymbol {x}})+{\frac {\partial ^{2}u}{\partial x_{2}^{2}}}({\boldsymbol {x}})\right)u^{*}({\boldsymbol {x}})d\Omega =0} この恒等式を2回部分積分すると、 q ∗ = ∂ u ∗ / ∂ n {\displaystyle q^{*}=\partial u^{*}/\partial n} として、 ∫ Γ q ( x ) u ∗ ( x ) d Γ x − ∫ Γ u ( x ) q ∗ ( x ) d Γ x + ∫ Ω u ( x ) ( ∂ 2 u ∗ ∂ x 1 2 ( x ) + ∂ 2 u ∗ ∂ x 2 2 ( x ) ) d Ω x = 0 , {\displaystyle \int _{\Gamma }q({\boldsymbol {x}})u^{*}({\boldsymbol {x}})d\Gamma _{x}-\int _{\Gamma }u({\boldsymbol {x}})q^{*}({\boldsymbol {x}})d\Gamma _{x}+\int _{\Omega }u({\boldsymbol {x}})\left({\frac {\partial ^{2}u^{*}}{\partial x_{1}^{2}}}({\boldsymbol {x}})+{\frac {\partial ^{2}u^{*}}{\partial x_{2}^{2}}}({\boldsymbol {x}})\right)d\Omega _{x}=0,} を得る。なお、下添字 x は点x に関する積分であることを表している。 得られた式を見ると、領域積分が1つだけ残っている。これを消し去るために、関数u* は次式を満足するように与える。 ∂ 2 u ∗ ∂ x 1 2 ( x ; ξ ) + ∂ 2 u ∗ ∂ x 2 2 ( x ; ξ ) + δ ( x − ξ ) = 0 , {\displaystyle {\frac {\partial ^{2}u^{*}}{\partial x_{1}^{2}}}({\boldsymbol {x}};{\boldsymbol {\xi }})+{\frac {\partial ^{2}u^{*}}{\partial x_{2}^{2}}}({\boldsymbol {x}};{\boldsymbol {\xi }})+\delta ({\boldsymbol {x}}-{\boldsymbol {\xi }})=0,} この関数u* は基本解と呼ばれており、ラプラス問題では、 r = | x − ξ | {\displaystyle r=|{\boldsymbol {x}}-{\boldsymbol {\xi }}|} （2点間の距離）として、 u ∗ ( x ; ξ ) = − 1 2 π ln ⁡ r , {\displaystyle u^{*}({\boldsymbol {x}};{\boldsymbol {\xi }})=-{\frac {1}{2\pi }}\ln r,} で与えられる。u* は、空間内の点ξに大きさ 1 の単位わきだしがあったときの、点x で観測されるポテンシャル値を与える関数、と解釈することができる。そのため、 ξ {\displaystyle {\boldsymbol {\xi }}} はソース点、 x {\displaystyle {\boldsymbol {x}}} は観測点と呼ばれる。 この定義式を上の積分方程式に代入すると、 u ( ξ ) = ∫ Γ u ∗ ( x ; ξ ) q ( x ) d Γ x − ∫ Γ q ∗ ( x ; ξ ) u ( x ) d Γ x , ( ξ ∈ Ω ) {\displaystyle u({\boldsymbol {\xi }})=\int _{\Gamma }u^{*}({\boldsymbol {x}};{\boldsymbol {\xi }})q({\boldsymbol {x}})d\Gamma _{x}-\int _{\Gamma }q^{*}({\boldsymbol {x}};{\boldsymbol {\xi }})u({\boldsymbol {x}})d\Gamma _{x},\quad ({\boldsymbol {\xi }}\in \Omega )} を得る。この式は、境界上のポテンシャルとフラックスの分布が得られている時に、領域内部の点x におけるポテンシャル値を計算する際に用いることができる。ラプラス問題や静弾性問題などでは、観測されるポテンシャル値に及ぼす境界上の解の変動の影響は距離が離れると共に小さくなるため、境界要素法によると内部の点でのポテンシャル値は精度よく計算できると考えられている。 次に、この積分方程式において、ポテンシャル値を評価する点ξを領域内部から境界上の点に移動させる。基本解u* ，q* は r = 0 で関数値が無限大に発散するため、点ξの境界上への移動は、境界積分が有限確定値となるように注意しながら、極限の意味で考える必要がある。その結果、先に示した積分方程式は、極限操作によって次のようになる。 c ( ξ ) u ( ξ ) + ∫ Γ q ∗ ( x ; ξ ) u ( x ) d Γ x − ∫ Γ u ∗ ( x ; ξ ) q ( x ) d Γ x = 0 , ( ξ , x ∈ Γ ) {\displaystyle c({\boldsymbol {\xi }})u({\boldsymbol {\xi }})+\int _{\Gamma }q^{*}({\boldsymbol {x}};{\boldsymbol {\xi }})u({\boldsymbol {x}})d\Gamma _{x}-\int _{\Gamma }u^{*}({\boldsymbol {x}};{\boldsymbol {\xi }})q({\boldsymbol {x}})d\Gamma _{x}=0,\quad \quad ({\boldsymbol {\xi }},{\boldsymbol {x}}\in \Gamma )} ここで、c (ξ) は、点ξの境界形状から決まる定数で、境界が滑らかであれば 1/2、かど点であれば当該点での内角の大きさから与えられる。この式が境界積分方程式 (Boundary Integral Equation) であり，境界要素法の離散化の出発点となる重要な方程式である。 境界要素法は、先に示した境界積分方程式を離散化し、近似解を得るための方法である。離散化においては、 境界上の未知量（ラプラス問題であればポテンシャルu とフラックスq ）の近似 近似関数を代入した後に得られる（積分方程式の）残差方程式の取り扱い 境界形状の近似 境界上の積分計算 が必要となる。以下では、先に取り上げた2次元ラプラス問題を例に、順を追ってその内容を説明する。 先に示したように、2次元ラプラス問題ではポテンシャルu とフラックスq とが変数（未知量）である境界積分方程式の離散化が必要になる。そこでまず、u とq とをN 個の補間関数を用いて近似する。 u ( x ) ≈ u ~ ( x ) = ∑ j = 1 N ϕ j ( x ) U j , q ( x ) ≈ q ~ ( x ) = ∑ j = 1 N ϕ j ( x ) Q j , {\displaystyle u({\boldsymbol {x}})\approx {\tilde {u}}({\boldsymbol {x}})=\sum _{j=1}^{N}\phi _{j}({\boldsymbol {x}})U_{j},\quad \quad q({\boldsymbol {x}})\approx {\tilde {q}}({\boldsymbol {x}})=\sum _{j=1}^{N}\phi _{j}({\boldsymbol {x}})Q_{j},} ここで、φj は補間関数であり、有限要素法で用いられている曲線要素や三角形要素、四辺形要素などをそのまま利用できる。なお、境界要素近似においては、定式化の上で特段の制約がない限り、区間一定近似の導入が可能である。その簡易さと境界積分の計算のしやすさから、多くの場合で区間一定近似が用いられる。 上で示した近似関数を境界積分方程式に代入すると、次の残差方程式を得る。 R ( ξ ) := c ( ξ ) u ~ ( ξ ) + ∫ Γ q ∗ ( x ; ξ ) u ~ ( x ) d Γ x − ∫ Γ u ∗ ( x ; ξ ) q ~ ( x ) d Γ x ≠ 0. {\displaystyle R({\boldsymbol {\xi }}):=c({\boldsymbol {\xi }}){\tilde {u}}({\boldsymbol {\xi }})+\int _{\Gamma }q^{*}({\boldsymbol {x}};{\boldsymbol {\xi }}){\tilde {u}}({\boldsymbol {x}})d\Gamma _{x}-\int _{\Gamma }u^{*}({\boldsymbol {x}};{\boldsymbol {\xi }}){\tilde {q}}({\boldsymbol {x}})d\Gamma _{x}\neq 0.} 近似関数をN 個の補間関数を用いて定義したことに注意して、残差方程式に対して次のいずれかの条件を課し、N 個の（積分）方程式を導出する。 境界上にN 個の代表点（選点）ξi (i = 1, 2, ... , N ) を置き、この各点で残差について R ( ξ i ) = 0 {\displaystyle R({\boldsymbol {\xi }}_{i})=0} であることを求める（選点法）。 N 個の補間関数φi (i = 1, 2, ... , N ) と残差方程式との境界積分を考え、各々が全て 0 となることを求める（ガラーキン法（英語版））。 境界要素法では、前者の選点法を採用して離散化を進めるのが一般的である。その結果、 R ( ξ i ) = ∑ j = 1 N [ c ( ξ i ) ϕ j ( ξ i ) + ∫ Γ q ∗ ( x ; ξ i ) ϕ j ( x ) d Γ x ] U j − ∑ j = 1 N [ ∫ Γ u ∗ ( x ; ξ i ) ϕ j ( x ) d Γ x ] Q j = 0 , {\displaystyle R({\boldsymbol {\xi }}_{i})=\sum _{j=1}^{N}\left[c({\boldsymbol {\xi }}_{i})\phi _{j}({\boldsymbol {\xi }}_{i})+\int _{\Gamma }q^{*}({\boldsymbol {x}};{\boldsymbol {\xi }}_{i})\phi _{j}({\boldsymbol {x}})d\Gamma _{x}\right]U_{j}-\sum _{j=1}^{N}\left[\int _{\Gamma }u^{*}({\boldsymbol {x}};{\boldsymbol {\xi }}_{i})\phi _{j}({\boldsymbol {x}})d\Gamma _{x}\right]Q_{j}=0,} ここで、 G i j = ∫ Γ u ∗ ( x ; ξ i ) ϕ j ( x ) d Γ , H i j = c ( ξ i ) ϕ j ( ξ i ) + ∫ Γ q ∗ ( x ; ξ i ) ϕ j ( x ) d Γ x {\displaystyle G_{ij}=\int _{\Gamma }u^{*}({\boldsymbol {x}};{\boldsymbol {\xi }}_{i})\phi _{j}({\boldsymbol {x}})d\Gamma ,\quad H_{ij}=c({\boldsymbol {\xi }}_{i})\phi _{j}({\boldsymbol {\xi }}_{i})+\int _{\Gamma }q^{*}({\boldsymbol {x}};{\boldsymbol {\xi }}_{i})\phi _{j}({\boldsymbol {x}})d\Gamma _{x}} とおくと、次のN 元の代数方程式を得る。 ∑ j = 1 N H i j U j = ∑ j = 1 N G i j Q j , ( i = 1 , 2 , … , N ) {\displaystyle \sum _{j=1}^{N}H_{ij}U_{j}=\sum _{j=1}^{N}G_{ij}Q_{j},\quad \quad (i=1,2,\ldots ,N)} なお、解の一意性が保証される場合では、境界値Uj , Qj はどちらか一方が未知で、もう一方が既知である。そのため、未知境界値をまとめてXj 、未知境界値に乗じられている係数をAij 、既知境界値と係数成分との乗算結果をまとめてbi で表すと、次の連立一次方程式を得る。 ∑ j = 1 N A i j X j = b i , ( i = 1 , 2 , … , N ) {\displaystyle \sum _{j=1}^{N}A_{ij}X_{j}=b_{i},\quad \quad (i=1,2,\ldots ,N)} この式を解くことで、境界上のポテンシャルとフラックスが近似的に得られることになる。 前節で示したように、境界積分方程式の離散化においては境界値の近似と境界積分の計算が必要となる。その際，当然のことながら物体形状も定義しておく必要がある。領域形状の近似表現においても、有限要素法で用いられる曲線要素や平面・曲面要素がそのまま利用できる。ただし、境界値の近似では区間一定近似が適用可能であったが、領域形状の近似においては区間一定近似を用いることはできない。 先に述べた選点法で積分方程式を代数方程式に置き換える場合、次の境界積分の計算が必要となる。 ∫ Γ u ∗ ( x ; ξ i ) ϕ j ( x ) d Γ x , ∫ Γ q ∗ ( x ; ξ i ) ϕ j ( x ) d Γ x . {\displaystyle \int _{\Gamma }u^{*}({\boldsymbol {x}};{\boldsymbol {\xi }}_{i})\phi _{j}({\boldsymbol {x}})d\Gamma _{x},\quad \int _{\Gamma }q^{*}({\boldsymbol {x}};{\boldsymbol {\xi }}_{i})\phi _{j}({\boldsymbol {x}})d\Gamma _{x}.} これらの積分は、被積分関数u* , q* が x = ξi で無限大となる特異性がある。境界要素解析において満足のいく結果を得るためには、この「特異性」を示す関数の積分をいかに精度よく、効率よく処理するかが重要である。この特異積分は、可能ならば解析的に（手計算で）処理し、不可能ならば特異性を除去した上で数値的に処理するか、または剛体移動条件や一定ポテンシャル条件などの物理的に満たさねばならない条件を用いて間接的に計算することになる。 なお、境界積分は、選点が境界上にない場合でもその取り扱いには注意を要する。選点と積分領域との距離が積分領域の代表長さに比べて小さい場合には、被積分関数が積分領域内で大きく変動し、ガウス求積などの数値積分公式を用いて積分計算を実行する場合に積分精度が大幅に低下することがある。境界積分の計算は係数行列の作成において必要となり、その積分誤差が大きくなると近似解の誤差も増大する。改善のためには、積分領域を細分割して積分を計算する方法が最も簡単である。 主要な境界値問題・初期値境界値問題における基本解は、次の通りである。 ラプラス問題（ポテンシャル問題） r = | x − ξ | {\displaystyle r=|{\boldsymbol {x}}-{\boldsymbol {\xi }}|} ：ソース点と積分点との距離、 u ∗ ( x ; ξ ) = − 1 2 π ln ⁡ r , (for 2D) , u ∗ ( x ; ξ ) = 1 4 π r , (for 3D) {\displaystyle u^{*}({\boldsymbol {x}};{\boldsymbol {\xi }})=-{\frac {1}{2\pi }}\ln r,\quad {\text{(for 2D)}},\quad \quad u^{*}({\boldsymbol {x}};{\boldsymbol {\xi }})={\frac {1}{4\pi r}},\quad {\text{(for 3D)}}} 静弾性問題（等方均質の場合, Kelvin解） r = | x − ξ | {\displaystyle r=|{\boldsymbol {x}}-{\boldsymbol {\xi }}|} ：ソース点と積分点との距離，E ：ヤング率、ν：ポアソン比、 G = E 2 ( 1 + ν ) {\displaystyle G={\frac {E}{2(1+\nu )}}} ：せん断弾性係数、δij ：クロネッカーのデルタ、 u i j ∗ ( x ; ξ ) = − 1 8 π ( 1 − ν ) G [ ( 3 − 4 ν ) ln ⁡ r δ i j − ∂ r ∂ x i ∂ r ∂ x j ] , (for 2D) {\displaystyle u_{ij}^{*}({\boldsymbol {x}};{\boldsymbol {\xi }})=-{\frac {1}{8\pi (1-\nu )G}}\left[(3-4\nu )\ln r\delta _{ij}-{\frac {\partial r}{\partial x_{i}}}{\frac {\partial r}{\partial x_{j}}}\right],\quad {\text{(for 2D)}}} ただし、この式は平面ひずみ問題の基本解である。平面応力問題の場合には、νを ν ′ = ν / ( 1 + ν ) {\displaystyle \nu '=\nu /(1+\nu )} に置き換えて基本解を構成すればよい。 u i j ∗ ( x ; ξ ) = 1 16 π ( 1 − ν ) G r [ ( 3 − 4 ν ) δ i j + ∂ r ∂ x i ∂ r ∂ x j ] , (for 3D) {\displaystyle u_{ij}^{*}({\boldsymbol {x}};{\boldsymbol {\xi }})={\frac {1}{16\pi (1-\nu )Gr}}\left[(3-4\nu )\delta _{ij}+{\frac {\partial r}{\partial x_{i}}}{\frac {\partial r}{\partial x_{j}}}\right],\quad {\text{(for 3D)}}} 定常スカラー波動問題 境界要素法には、以下のような特徴、および利点・欠点がある。 境界要素法の最大の特徴は、対象とする問題によっては「境界上の離散化のみで近似解が得られる」ことにある。境界上の離散化は、3次元問題ならば曲面上、2次元問題ならば曲線上で行われる。そのため、有限要素法のように領域内の離散近似が必要な解法と比べ、離散化に必要な要素や節点の数が少なくて済む。 境界上の離散化だけで問題が解ける場合としては、静的問題・定常問題では ラプラス問題、線形弾性問題、定常波動問題などのように、線形問題で離散化の際に用いられる基本解が解析的に厳密に得られ、かつ内部ソースや物体力のような支配方程式の非同次項が存在しない場合である。ただし、支配方程式が非同次項を含んでいても常にこの特徴が失われる訳ではなく、非同次項の種類によっては非同次項を含む領域積分を境界積分に変換できる場合もある（例：線形弾性問題における重力の作用）。 時間発展型問題において境界上の離散化のみで近似解を得るためには、線形問題の際に課された条件の他に、時間に関する離散化方法にも注意が必要である。具体的には、与えられた問題に対応する時間と空間に関する積分方程式（時間域積分方程式）を定式化の出発点とし、空間・時間双方を離散化した上で、当該の定式化の下での基本解と初期条件との領域積分（定式化の結果として残る積分項）が消滅するか、または境界積分に置換可能な場合に限り、時間発展問題の境界要素解析でも境界上の離散化だけで近似解が得られることになる。有限要素解析や差分計算の場合のように、時間方向の離散化を時間積分法で近似的に処理すると、解析における各時刻において現時刻での場の値と基本解とを含む領域積分が生じ、上述の特徴は失われてしまうことになる。 なお、境界要素法は、幾何学的非線形問題や材料非線形問題のように、領域内部で満たすことを求められる支配方程式や構成方程式そのものに非線形性がある場合でも近似解を得ることが可能ではある。しかし、定式化の取り扱いの中で領域積分が副次的に生じ、境界要素法の最大の利点である「境界上の離散化だけで近似解が得られる」点が失われてしまい、現在ではあまり用いられない。 離散化に用いる要素や節点の数、場の変数の評価点の数が小さくなれば、最終的に得られる代数方程式（多元連立（一次）方程式）の規模（元数、未知量の総数）も小さくなる。線形問題・非線形問題を問わず、汎用の離散化解析手法では支配方程式を最終的に連立一次方程式に帰着させ、この方程式の解から近似解を構成するため、連立方程式の元数の大小は解析時の計算負荷（使用メモリ、計算時間）に直結する。当然のことながら、問題の規模を小さくすれば、計算負荷はより小さく抑えられることになる。 上述のように、境界要素法では規模の小さい連立一次方程式を取り扱うことができるものの、方程式の係数行列の成分はほぼ全て 0 でないものとなる。そのため、係数行列の保存に要する記憶量は方程式の元数N に比例する。また、連立方程式の解を得るためには、ガウスの消去法に代表される直接法を用いれば N3 に、反復法 (数値計算)を用いてもN2 に比例する計算量が必要となる。 領域内の離散化が必要となる有限要素法や有限差分法では、係数行列の成分のほとんどが 0 である疎行列となるため、多少問題の規模が大きくなっても使用メモリや計算量は境界要素法と比べて少なくて済む場合が少なくない。そのため、この点は境界要素法の最大の欠点の一つとして考えられている。解決策としては、多体問題の解析の高速化に用いられていた高速多重極展開法（英語版）の適用や、ウェーブレットの利用が提案されている。 上述のような当該解法の特徴，および利点・欠点を考慮して，今日において境界要素法が得意とする問題としては以下のものがある． 波動伝搬問題とは、対象とする領域内で物理量の擾乱が「波動」として有限な速さで伝播していく問題であり、その多くは開領域の問題（領域に無限遠を含む問題）または半無限領域の問題として定義されることが多い。境界要素法では、開領域の問題をそのまま取り扱うことができ、特に波動問題では、無限遠での波動の放射が近似処理なしに表現できる。有限差分法や有限要素法では動的応答の観測点から十分離れたところに仮想的に境界を設け、そこでは波動の放射を表現するような近似的な取り扱いが必要となる。そのような点から、境界要素法は地盤振動解析や地震波の伝播解析、音響問題の解析、電磁場解析などで用いられることが多い。ただし、閉領域を対象とした動的問題（振動問題など）においては、有限要素解析の場合のようなモード解析ができない上、有限要素法と比べて計算時間を要することから、あまり多用されないようである。 境界要素法の利点の1つに、境界上の離散化だけで問題を解くことができる点があった。形状最適化問題とは、工学分野の構造部材の形状を、所定の目的関数と制約条件の下で自動的に最適化する問題である。部材の供用を弾性限界内に考えた場合、弾性応答は境界積分方程式を解くことで把握でき、設計感度の計算も同様となる。感度計算は形状の変更のたびに必要となるが、境界上の離散化のみでよいため、要素分割等の作業の手間を大幅に削減することができる。 数値解析 CAE 有限要素法 代用電荷法