■境界要素法でTorsionのPoisson’s equationを解く■
境界要素法で今回の問題を解いてみました。 と言いてもPoisson’s equationのGreen関数は無いのでLaplace EquationのGreen関数を使い、2Gθは領域積分します。 通常、境界要素法では境界のみを要素分割しますが、今回は領域も要素分割する必要があります。

境界要素法は有限要素法の重み付け残差法と同様に、下式に示す様に残差×重み関数の領域積分からスタートします。

上式のG(ξ,x)がGreen関数で、Laplace Equationの場合Green関数は以下です。

ここにξはSource point、xをObservation point と言います。 これらの2つの点間の距離が上式のrになります。 下図を見て下さい。ξは領域内にxが境界に有る場合を示しています。 今回の場合、2Gθの領域積分があるのでxは領域内にも存在します。 境界要素法では記述を簡略化するために、点を表すときに(x,y)とせずにベクトルxで表します。

有限要素法の場合は、変数の値を計算したい節点がSource pointです。その点を節点Aとすると、節点Aを共有している要素の節点がObservation pointになります。節点Aを共有していない要素はObservation pointになりません。境界要素法では、上のGreen関数で分かる通り、全ての境界要素上の節点がObservation pointになります。 有限要素法と境界要素法の相違点ですが、もう1つあります。境界要素法では部分積分を2回使います。有限要素法では1回でした。境界要素法の詳細についてはサイトを参考にしてください。 まず、1回目の部分積分を上の式に施します。結果は以下です。

\begin{eqnarray} \oint{\frac{\partial\phi\left(\boldsymbol{x}\right)}{\partial x_i}G\left(\boldsymbol{\xi},\boldsymbol{x}\right)n_ids}-\iint{\frac{\partial\phi\left(\boldsymbol{x}\right)}{\partial x_i}\frac{\partial G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i}dA+2G\theta\iint{G\left(\boldsymbol{\xi},\boldsymbol{x}\right)dA=0}} \end{eqnarray}

\begin{eqnarray} \oint{\frac{\partial\phi\left(\boldsymbol{x}\right)}{\partial x_i}G\left(\boldsymbol{\xi},\boldsymbol{x}\right)n_ids}-\oint{\phi\left(\boldsymbol{x}\right)\frac{\partial G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i}n_ids}+\iint{\phi\left(\boldsymbol{x}\right)\frac{\partial^2G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i\partial x_i}dA+2G\theta\iint{G\left(\boldsymbol{\xi},\boldsymbol{x}\right)dA=0}} \end{eqnarray}

\begin{eqnarray} \frac{\partial G}{\partial r}=-\frac{1}{2\pi}\left(\frac{1}{r}\right) \end{eqnarray}

\begin{eqnarray} \frac{\partial^2G}{\partial r\partial r}=-\frac{1}{2\pi}\left(-\frac{1}{r^2}\right) \end{eqnarray}

\begin{eqnarray} \frac{\partial^2G}{\partial\theta\partial\theta}=0 \end{eqnarray}

\begin{eqnarray} \iint{\phi\left(\boldsymbol{x}\right)\frac{\partial^2G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i\partial x_i}dA}=\phi\left(\boldsymbol{\xi}\right)\iint{\frac{\partial^2G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i\partial x_i}dA}=\phi\left(\boldsymbol{\xi}\right)\oint{\frac{\partial G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i}n_ids} \end{eqnarray}

\begin{eqnarray} \phi\left(\boldsymbol{\xi}\right)\oint{\frac{\partial G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i}n_ids}=\phi\left(\boldsymbol{\xi}\right)\int_{0}^{2\pi}{\frac{dG\left(\boldsymbol{r}\right)}{dr}\varepsilon}d\theta=\phi\left(\boldsymbol{\xi}\right)2\pi\left(-\frac{1}{2\pi}\frac{1}{\varepsilon}\right)\varepsilon=-\phi\left(\boldsymbol{\xi}\right) \end{eqnarray}

\begin{eqnarray} \phi\left(\boldsymbol{\xi}\right)=\oint{\frac{\partial\phi\left(\boldsymbol{x}\right)}{\partial x_i}G\left(\boldsymbol{\xi},\boldsymbol{x}\right)n_ids}-\oint{\phi\left(\boldsymbol{x}\right)\frac{\partial G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i}n_ids}+2G\theta\iint G\left(\boldsymbol{\xi},\boldsymbol{x}\right)dA \end{eqnarray}

\begin{eqnarray} q_n\left(\boldsymbol{x}\right)=\frac{\partial\phi\left(\boldsymbol{x}\right)}{\partial x_i}n_i \end{eqnarray}

\begin{eqnarray} F\left(\boldsymbol{\xi},\boldsymbol{x}\right)=\frac{\partial G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i}n_i \end{eqnarray}

\begin{eqnarray} \phi\left(\boldsymbol{\xi}\right)=\oint{q_n\left(\boldsymbol{x}\right)G\left(\boldsymbol{\xi},\boldsymbol{x}\right)ds}-\oint\phi\left(\boldsymbol{x}\right)F\left(\boldsymbol{\xi},\boldsymbol{x}\right)ds+2G\theta\iint G\left(\boldsymbol{\xi},\boldsymbol{x}\right)dA \end{eqnarray}

\begin{eqnarray} -1=\oint{\frac{\partial G\left(\boldsymbol{\xi},\boldsymbol{x}\right)}{\partial x_i}n_ids}=\oint F\left(\boldsymbol{\xi},\boldsymbol{x}\right)ds=\int_{0}^{2\pi}{\frac{dG\left(\boldsymbol{r}\right)}{dr}\varepsilon}d\theta \end{eqnarray}

\begin{eqnarray} \oint F\left(\boldsymbol{\xi},\boldsymbol{x}\right)ds=\int_{0}^{\pi}{\frac{dG\left(\boldsymbol{r}\right)}{dr}\varepsilon}d\theta=-0.5 \end{eqnarray}

■SET4-TORSION-QUARTER.FOR■
とりあえず有限要素法のtorsion4q.forまたはtorsion9q.forで取り上げたのと同じ問題を解いてみましょう。 下はNEX=10とNEY=10を用いた境界の分割図です。境界は1次要素で分割してあります。計算領域はfemと同じ1/4モデルです。

■モーメントの計算■
有限要素法のときと同様にここでもモーメントを計算できるようにプログラムを書きました。 モーメントを計算するために下式の領域積分を行う必要があります。 今のプログラムでは2つの方法でモーメントを計算するプログラムを書きました、１つはφの描画の値を使って計算しています。 φの値は領域の要素分割の節点で境界要素式を用い計算されています。 それらの値と形状関数を使ってモーメントを計算します。結果はファイルTORQUEBEM.OUTに入ります。

■角問題(Corner Problems)■
有限要素法では発生しなくて境界要素法で発生する現象が角問題(Corner Problems)です。 理由は、温度、密度、Stress Functionφ等のスカラーを微分するとベクトルになります。 境界でスカラー(φ)を指定するのがDirichlet境界条件で、 スカラーの法線微分(q_n=∂φ/∂n)を指定するのがNeumann境界条件です。 境界要素法では、各境界要素においてDirichlet境界条件かNeumann境界条件かを指定しなければなりません。 要素を線積分すると結果は節点に集約されます。 直線状の境界(スムース境界と言う)にこの節点(内角180度)があれば問題ないです。 問題になるのが、Dirichlet境界条件を持つ節点がスムース境界に無い(内角180度以外)ときです。 下図はStress Functionφが既知のスムース境界での節点と要素を表しています。 つまり、節点での内角が180度になっています。

この場合、要素i-1での q_n2と要素iでのq_n1の方向がそれぞれ同じですよね。 ですので、節点iでStress Functionφの法線方向の微分が1つですので問題無しです。

角に2つの節点を置き要素iを設けます。この要素の長さは出来るだけ短くします。 すると問題になっていた角のqnは要素i-1と要素i+1に分けることが出来ます。節点Iの内角は135度になり、90度から少しスムースな境界になっています。 境界要素法では一定要素があります。 要素の中央で境界値を設定します。 この場合、スムースな境界で未知数が計算されますので角問題は発生しません。 BEMサイト では、もう１つ角問題フリーのMIXED要素も紹介してあります。 この要素ではDirichlet境界を節点で指定し、Neumann境界は要素の中央で指定します。

■例題の計算■
角問題についてお話しましたが、角問題の処理を行わず計算を行ってみました。 NEX=NEY=10で境界および領域を要素分割しました。 Shear Modulus G=26538461.5384615[Pa]、θ= 0.0000670055862315403を使いました。
下は同じ条件のFEMの結果をBEMの結果に重ねて表示してあります。

■Warpingの解析■
WarpingについてはTorsion理論のところ示した通り、Neumann境界は文献に示されていましたが、解析に必要なDirichlet境界条件については不明でした。しかし、 Torsionの3次元弾性解析でシャフトが矩形の場合、 断面の4隅のシャフト長手方向の変位がゼロで有ることが判明しました。この事実を使ってWarpingの解析を行ってみます。 まず解かなければ式は下のLaplace Equationでした。

■Warping計算の実行■
境界条件が明白になったので、計算を行ってみます。まず計算に必要なデータを作成する必要があります。ここでは、 SET4-WARPING-FULL.FOR が作成してくれます。そしてLaplace Equationは BEM8LINQ-WARPING.FOR が解いてくれます。以下にこれらのプログラムが入力および出力するファイルを以下に示します。