数理解析研究所講究録 1441 巻 2005 年 144-158

144

応用力学における多重極法について-

周期境界値問題への適用

京都大学・学術情報メディアセンター-

西村 直志 (Naoshi Nishimura)Academic Center for Computing and Media Studies,
Kyoto University

1

はじめに

境界積分方程式法は、 有限要素法や差分法と並んで、 初期値・境界値問題を数値的に解く代表的な計算力学の手法の一つである。 考える領域の境界のみに着目して解析を実行できるという点で、 外部問題、 或いはクラック問題に適している。 しかし、 最終的に帰着される線形連立一次方程式の係数行列が密となるため、 それを例えば反復法を用いて解くと要素数 $N$ に対して

$N^{2}$のオーダーの計算量が必要となり、 大規模問題に適用する際の難点となっていた。 それを克服したのが高速多重極法 [$1|$に代表される高速解法である。 境界積分方

程式法に高速多重工法を適用すれば、 計算量や必要なメモリの量を $O(N)$ にまで下げることができる。
一方、 近年、 カーボンナノチューブに代表される新しい素材を用いた複合材料が開発されており、 その力学的特性の解明が工学上の興味となっている。複合材料のような微視構造を有する物質の巨視的な力学特性を導く方法として、 微視構造に比べてある程度大きい解析領域をとって通常の力学解析を適用する事も行なわれているが$[2]_{\text{、}}$

精度や効率を考えると、 均質化法 $[3, 4]$ による方法が有望である。 均質化法では、 周期的な微視構造をもつ物質が巨視的な構造を成していると仮定し、 マクロとミクロのスケールで漸近展開を行なうことによって、 微視的な構造を反映させた巨視的な解析を行なうが、 その際に、 微視的な構造に関する周期境界値問

題を解くことが必要になる。

微視構造がクラックなどの鋭い形状変化部を含む場合、 その解析手法として境界積分方程式法が有利である。 しかし、 微視構造が複雑な形状を有している時、 従来法を用いたのでは最初に述べた $O(N^{2})$ の計算量がネックとなって、 十分に大きな解析がおこなえない可能性がある。 従って、 周期境界値問題に対する高速解法を研究することは有意義であると考えられる。 これまで均質化法に境界積分方程式法を適用した論文はいくつかあるが (Koro and Abe [5] の文献参照)、 特に高速解法に関する研究として Koro and Abe[5] を挙げることが出来る。 この研究では wavelet 境界積分方程式法によって速度向上を図っている。 ただし、 ユニットセルの外周に周期境界条件を科すタイプの手法であるので、 通常の境界値問題に比べて行列の圧縮が十分に行なえない傾向がある。 そこで本稿では、 外周の境界条件を考慮する必要のない高速多重極境界積分方程式法を定式化する。 これまで、 周期境界条件の下で高速多重極法を用いた研究として、 Laplace 方程式における Greengardと Rokhtin[6] の初期の仕事や、 弾性学における Greengard and Helsing[7] を上げることが出来るが、 その数

学的取り扱いには不明瞭な部分が残っている。 そこで本稿では 2 次元 Laplace 方程式$[8]_{\text{、}}$

及び弾性学 [9] のクラック問題を取り上げ、 その周期境界値問題における多重極法を、 数学的明瞭さを保ちつつ定式化する。 さらに、 均質化法を用いて多数のクラックから成る微視構造を有する物質の巨視的な弾性定数を求める。

145

2 定式化 (Laplace の場合)
2.1

境界積分方程式法

2 次元 Laplace 方程式におけるクラック周期境界値問題は、 以下のように表される。
基礎方程式

(in

$\Delta u=0$$\Omega\backslash S$(1))

クラック上での境界条件(on$\frac{c98_{-}’}{\partial n}=g$$S$)$(’\mathit{2})$周期境界条件$u(x^{1})$$\frac{\partial u(x^{1})}{\partial n}$

ここに、$\Omega$はユニットセル、$\Gamma$は$\Omega$$=$$=$$u(x^{2})$(on$\Gamma$)(3)$- \frac{\partial u(x^{2})}{\partial n}$(on$\Gamma$)(4)の外側の境界、$S$

はクラック、 9C よ $S$ 上の既知関数、は外向き、 クラックについては任意に定めた正の向きの法線方向微分であり、$x^{1}$と x2\sim よ$\frac{\partial}{\partial n}$$\Gamma$\iota よ境界につ U‘て-E で互 U‘に反対側の点である。 この問題の解は加算的な定数を除いて一意である。
$\mathrm{x}^{\mathrm{t}}$$\llcorner_{\mathrm{x}\prime}\mathrm{X}2$

図 1

次に、 この問題の$\text{解表}$現.
$\text{を}\backslash \ovalbox{\tt\small REJECT}.$く。 以下

2 次元クラック周期境界値問題$\mathrm{a}\mathrm{e}_{\mathrm{B}}\mathrm{g}^{1}\text{のために適当に}\ \mathrm{f}\wedge,\frac{\Phi}{\tau}\backslash \text{をとって_{、}}$

が 1 の正方形領域であるとする。

まず、 無限領域中にクラック $S$ がある場合を考える。 このとき、 解$u$$\phi$

はクラックの開口変位、$G$

は原点を中心とする一辺の長さの表現 G よ次のようになる。

$u(x)= \int_{S}\frac{\partial G(x-y)}{\partial n_{y}}\phi(y)dS_{y}$ここに、$\Omega$(5)

は 2 次元 Laplace 方程式の基本解で ‘$G(x-y)= \frac{1}{2\pi}\log\frac{1}{|x-y|}$

(6)

146

である。
$z=x_{1}+\mathrm{i}x_{2}$,$\xi=y_{1}+\mathrm{i}y_{2}$(7)とおくと、 式 (5) は、$u(x)=-{\rm Re}[ \frac{1}{2\pi \mathrm{i}}\int_{S}\frac{\phi(\xi)}{z-\xi}d\xi]$と書ける。 これは、 Weierstrass の$\zeta$関数$\zeta(_{\sim}’)=\frac{1}{z}+\sum_{w\in N}(\frac{1}{z-w}+\frac{1}{w}+\frac{z}{w^{2}})$を用いて疑周期化することができる$[10]_{\text{。}}$$N=$

なる複素数の集合である。

$\zeta(z)$

(8)

(9)

ここに、 $N$ は{ $l=m+\mathrm{i}n|m,$$n$: 整数, $l\neq 0$ }が$\zeta(z+1)=\pi-\tau\zeta \mathrm{l}(z)$$\zeta(z+\mathrm{i})=-\pi i+\zeta(z)$(10)

を満たすことに注意すると、 周期境界値問題の解表現は、$u(x)={\rm Re}[- \frac{1}{2\pi i}\int_{S}((z-\xi)\phi(\xi)d\xi+\frac{\overline{z}}{2i}\int_{S}\phi(\xi)d\xi]$(11)となることが分かる。 よって、$\frac{\partial u(x)}{\partial x_{1}}-\dot{\iota}\frac{\partial u(x)}{\partial x_{2}}=\frac{1}{2\pi\dot{x}}(\int_{S}\wp(z-\xi)\phi(\xi)d\xi-\pi\overline{M}_{0}(0))$(12)

が得られる。 ここに、$\wp(z)=-\zeta’$ $(z)= \frac{1}{z^{2}}+\sum_{w\in N}(\frac{1}{(z-w)^{2}}-\frac{1}{w^{2}})$は Weierstrass の

$\wp$

(13)

関数であり、 $M_{0}(0)$ はユニットセルにおける 0 次の多重極モーメント$M_{0}(0)= \int_{S}\phi d\xi$,

および、 は複素共役をあらわす。
-

22 高速多重極法
式 (12) の右辺の積分を計算するために、 高速多重極法を用いる。 Greengard と Rokhlin の有名な最初の論
文 [6] は、逸早く周期境界条件を取り扱っているが、発散項を物理的考察により有限にする「くりこみ」を使っているために数学的に意味不明な点を残している。 以下では、 基本的には Greengard-Rokhlin のアイデアに従いながら、 式 (12) に基づいて多重極法を定式化し、 数学的不明瞭さを取り除くこととする。 また、$I_{p}(z)= \frac{z^{p}}{p!}$$O_{p}(z)= \frac{p!}{z^{1[perp]_{\mathrm{P}}}}$

(14)

147

という記号を用いる。

まず、$\Omega$

をレベル 0 のセルとし、 多重極モーメントる。 ここに、$\xi_{0}$はレベル 0, 1 を含む全レベルで計算することにす$M_{p}$

をセル $C$ の中心とすると、(15)$M_{q}( \xi_{0})=\int_{C}I_{q}(\xi-\xi_{0})\phi(\xi)d\xi$である。 多重極モーメントの計算法は、 通常の多重調法のものと全く同じである。

この際 レベ次に、 レベル 0 のセルから始まって下向きに、全てのセルの中心での局所展開係数を求める。
関数の寄与のうち、 式 (13) の和の中のル 0 のセルの局所展開係数は、 式 (12) への $M0$ の寄与と、$\wp$$N’=N\backslash N’’$からの寄与を考える。 ここに$N”=\{l=m+in| m, n\in\{-1,0,1\}_{?}l\neq 0\}$である。 従って、 レベル 0 のセルの局所展開係数は、 式 (12) の右辺中括弧内の該当部分を$\sum_{q=0}^{\infty}L_{q+1}(0)I_{q}\langle z)$の形に局所展開したとすると、(16)$\sum(-1)^{p-q}O_{p}(w)M_{p-q}(0)-\delta_{q1}\pi\overline{M}_{0}(0)$$L_{q}(0)= \sum_{p=\max(2,q)w\in N’}^{\infty}$となる。 レベ J 口以下のセルの局所展開係数は、 レベル 0 の局所展開係数で考慮されていない項の式(12) へ

(同図での寄与を考慮して計算する。 これは図 2 において、 neighbour level 0cell と書かれている位置にの極法での局所と記されている) のレプリカが置かれてし ると考えた上で$\Omega$は$\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{n}^{r}\mathrm{a}$llevel0cell$\Phi\text{、}- \mathrm{r}\mathrm{a}\mathrm{e}$$\backslash$

係数の算法を適用して求められるものに外ならない。
式 (16) の評価には次式で表される lattice sum$\sum_{z\in N’}\frac{1}{z^{r}}$の値が$\prime A^{\backslash }\backslash \text{要}$である。 これらの和は、

$r$

が 4 の

$\mathrm{P}_{\mathrm{D}}^{\cdot}\text{数}$,$\tau=3,$ $\ldots$

であるとき以外は対より 0 となることが容易に分力 る。
$\mathrm{a}$$\Psi\prime \mathrm{f}^{\mathrm{J}}\mathrm{a}1\not\simeq$

数値計算においてはこれらの lattice sum の値はあらかじめ計算しておき、それを用いれば良い。

23 均質化法
え ae\not\in めるために勇ここで以下の考察は Laplace \hslash \Xi 式で支配される{fJ=\Xi ‘の物 fE\Re 象に当てはまるが、, ( \Re i g“’|4 定数 $=1$ ) を想定し、 は 方向の変位を表すものとする。 その上で、 微 ] なスケー断$\mathrm{F}\mathrm{h}\text{考}$$7\mathrm{E}\ovalbox{\tt\small REJECT}_{\text{、}^{}\prime}$$\text{弾^{}\iota}\prime \mathbb{E}_{\acute{\mathrm{A}}}^{7}\Gamma\text{、}\pi’/$$u$$\Rightarrow \mathrm{R}^{|}$

ルにおいて、$x_{3}$

多数のクラックを含む物質を考える。各クラックは表面力を受けないものとする。

この時、 この物質の巨視的な弾‘を数は、 均質化法の fgp により次のように求められる [11](詳細 \iota fX$=\mathrm{r}\prime \mathrm{f}\mathrm{f}$$\mathrm{E}\text{定}$参照されたい)。
$\phi^{j}n_{i}dy$$\delta_{\dot{\mathrm{z}}j}+\frac{1}{|\Omega|}\int_{S_{c}}$

.

(17)

148

}

$\Phi$

the口neighbour口

far leve[0cells

図2$\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{n}\mathrm{a}|$$\mathrm{e}\mathrm{v}\mathrm{e}|0$$|\mathrm{e}\mathrm{v}\mathrm{e}|0$

cell

ce[ls

レベル 0 のセルの位置関係

ここに、 がは、 クラック上の境界条件が、$\frac{\partial\chi^{j}}{\partial n}=n_{j}$(onで与えられる Laplace 方程式の周期クラック境界値問題の解$\chi^{J}$$S_{c}$(18))の開口変位である。 以上より、 巨視的な弾性係数を評価するためには、 周期境界条件のもとでクラックの開口変位がを求めなければならないことが分かる。

3 数値計算 ($\mathrm{L}\mathrm{a}\mathrm{p}^{1}\mathrm{a}\mathrm{c}\mathrm{e}$の場合)提案する周期境界値問題の多重極法を用いて、 多数のクラックを含む物質の巨視的弾性定数を求めた。 以下
ではいくつかの数値例を示す。

3.1

問題設定、 記号など

まず、 種々の記号を定義する。以下では取り扱うユニットセルは-

般には複数の、 同一長さ $2a$ の直線クラックを含むものとする。 クラックは基本的に座標方向に整列した場合を取り扱い、$n_{x\text{、}}x_{2}$方向のクラック数を$n_{y}$$x_{1}$

方向のクラック数をとする。 クラックー本当たりの要素分割数は同一’ (n) とする。

148

32 直接法による解との比較
周期境界値問題における高速多重極法の解の妥当性を検証するために、 従来法を用いて求めた開口変位‘位 (FMM) \emptyset lbrik\epsilon 行なった。従来法では、によって求め(Direct) と、 高速の境界積分方程式を離散化し、 上の境界条件と 上の周期境界条件を直接離散化して連立一次方程式に帰の値を陽着させ、 開口変位を求めた。 その際、 解は定数の自由度を持っているため、 境界上の要素の一つで$\Gamma+S$$_{\vee}^{\wedge}\ovalbox{\tt\small REJECT}\# 7\text{口^{}\mathcal{X}\prime}\$$\text{多重極}\backslash \text{法}$

$S$

上の

$\Gamma$$u$に与えた。 条件は、 $a=0.4,$ $n_{x}=1,$ $n_{y}=1,$ $n=200$ でクラックの向きは
条件は$\underline{\partial u}=n_{2}=1.0$$x_{1}$

方向とし、 クラック上の境界とした。 得られた開口変位は図 3 のようになり、 異なる方法から得られた$\varpi i]^{i}$g の結果$\partial_{7l}$が一致したので、 本稿の方法の正しいことが示された。

図 3

直接法及び高速多重極法によって求めた開 [」変位

33 高速多重極法による巨視的な弾性係数の算定
1 辺の長さが 1 のユニットセル内に、 $a=0.04$ のクラックを 10並行$\mathrm{f}\mathrm{f}\mathrm{i}\ovalbox{\tt\small REJECT}$と一$\frac{\pi}{4}$

から

まで ]
$\sigma$$\frac{\pi}{4}$$\mathrm{F}\mathrm{c}$$\mathrm{x}10=100$

本 H 己置し、 これらにランダム b こ

囲でランダムな回転も与えた場合 (図 4) を考える。 本稿の高速多分方程式法を用いて巨視的な弾性定数$E_{ij}$$\sum..\mathrm{t}\underline{\varpi}$

境界積を求めたところ、 各成分は次のようになった。

$E_{ij}=(\begin{array}{ll}0.9140 0243300.02433 0.6486\end{array})$4

クラック問題の定式化 (弾性体の場合)
ここでは$\text{、}$ffl 弾性平面ひずみ問

配方程式を満たす ($\lambda,$$\mu$$\text{題}\#$:考える。

上にクラック $S$ が分布してし るとする。 変 (立$\mathrm{a}$$\mathbb{R}^{2}$は Lan)\’e 定数)。
$\mu u_{i,jj}+(\lambda+\mu)\tau x_{j,ji}=0$in$C_{i_{J}kl}u_{k,l}n_{j}=g_{i}$

ここに、

,

$u_{i}\iota\mathrm{h}^{\text{、_{}i\pi a\mathit{3}\leq}}$$C_{\dot{\tau.}jk}f=\lambda\delta_{ij}\delta_{kl}+\mu(\delta ik\delta jl+\delta il\delta jk)$である。

$\mathbb{R}^{2}\backslash S$

on

$S$

(19)
(20)

150

1– /—- $\backslash /-/$$\backslash$$09$’ $-\backslash /\sim-\backslash \sim\backslash$$0.8’\backslash -\backslash ’\sim/\sim--$$0.7/\backslash \sim/\backslash ’//_{\sim}/$$0.6$$-\backslash \vee-\backslash /---$,$0.\mathrm{s}\backslash /-’/\backslash -\sim//$$0A-\vee\vee-///-\sim-$$\mathfrak{a}.\mathrm{a}$$\backslash \backslash \backslash /\backslash /-\backslash /-$$0.2/–’\backslash -//-$$0.\{0_{00102030A0S0B0\mathrm{J}0\mathrm{B}0\mathrm{B}}\sim.-\sim\backslash \backslash /////$

図4,

cracks

二次元静弾性学の基本解は次のように書ける。
$G \text{り}=\frac{1}{2\pi\mu}[\delta_{ij}\log\frac{1}{r}+\frac{1}{8(1-\nu)}$

ここに、

$l/$

哉

$\partial_{j}$

(

(21)

$r^{2} \log\frac{1}{r}\rangle]$

はポアソン比である。基本解を用いると、 支配方程式と等価な $S$ 上の境界積分方程式は次のよう

に書ける。
.

$g_{\mathrm{r}\iota}=\ovalbox{\tt\small REJECT}_{S}^{c_{abik}c_{\mathrm{c}djl}\frac{\partial}{\partial x_{k}}\frac{\partial}{\partial y_{l}}G_{ij}(x-y)n_{b}(x)n_{c}(y)\phi_{d}(y)dS_{y}}$ここに、$\phi_{i}$はクラックの開口変位である。 また、 変位$u_{i}$(22)の解表示は次式で得られる。

$u_{i}(x)= \int_{S}C_{cdj\mathrm{J}}\frac{\partial}{\partial y_{l}}G_{ij}(x-y)n_{c}(y)\phi_{d}(y)dS_{y}$$x\in \mathbb{R}^{2}\backslash S$(23)

式 (23) の複素表現は次のように書ける。

$V= \frac{1}{2\pi(\kappa+1)}\int_{S}\ovalbox{\tt\small REJECT}\frac{\kappa\nu(\zeta)U(\zeta)}{z-\zeta}+\frac{\overline{l/(\zeta)U(\zeta)}(z-\zeta)}{(\overline{z}-\overline{\zeta})^{2}}+\frac{\iota/(()\overline{U(\zeta)}+\overline{lJ(\zeta)}U(\zeta)}{(\overline{z}-\overline{\zeta})}\ovalbox{\tt\small REJECT} dS_{\zeta}$(24).
ここに、 $V=u_{1}$ 十 $\mathrm{i}u_{2^{\text{、}}}U=\phi_{1}+i\phi_{2^{\text{、}}}z=x_{1}+\mathrm{i}x_{2^{\text{、}}}(=y_{1}+\mathrm{i}y_{2^{\text{、}}}l/=n1(y)+\mathrm{i}n_{2}(y)$ である。 また、$\kappa=3-4\iota/$

41

である。

多重極展開

次に、式 (24) の多重極展開を考える。以下では、 を観測点、$z$$\zeta$を積分点、$z_{0}$

を観測点側のセルの中心、$\zeta_{0}$を積分点側のセルの中心とする。 ここで、 $z-\zeta=$ を $(z-z_{0})+(z_{0}-\zeta_{0})+(\zeta_{0}-\zeta)$ の形に分けて、各項につ151

いて変数分離すると、 次のようになる。
$V= \frac{1}{2\pi(\kappa+1)}\sum_{q=0}^{\infty}(-1)^{q}\ovalbox{\tt\small REJECT}^{riI_{q}(z-z_{0})\sum_{p=q}^{\infty}O_{p}\int_{S}I_{p-q}\nu UdS_{\zeta}}$$-(z-z_{0})I_{q-1}(z-z_{0}) \sum_{p=q}^{\infty}O_{p}\int_{S}I_{p-q}(\zeta-(_{0})\nu UdS_{\zeta}$$+(z_{0}- \zeta_{0})I_{q}(z-z_{0})\sum_{p=q}^{\infty}O_{p+1}\int_{S}I_{p-q}(\zeta-(_{0})\nu UdS_{\zeta}$$-I_{q}(z-z_{0}) \sum_{p=q}^{\infty}O_{p}\int_{S}\overline{(\zeta-\zeta_{0})}I_{p-q-1}(\zeta-\zeta_{0})\iota/Ud\zeta$$+\overline{I_{q}(z-z_{0})\sum_{p=q}^{\infty}O_{p}\int_{S}I_{p-q}(\zeta-\zeta_{0})(\overline{l/}U+\nu\overline{U})dS_{\zeta}}||$

ただし、 $I_{-1}(z)=0$ と定義する。 なお、 $O$ の引数は全て$V= \frac{1}{2\pi(\kappa+1)}$[\kap a 重$z_{0}-\zeta_{0}$(25)である。 上式の形から$(z-z_{0})-l_{\backslash }z-z_{0})\overline{\Psi’(z-z_{0})}+\overline{\Phi(z-z_{\zeta\}})}]$(26)と書けることが分かる。
と書き、 の多重極モーメント、 局所展開係重の多重極モーメント、 局所展開係数をそれぞれと書くことにする。 式 (25) を用いると、 これらの多重極モーメントは次のように数をそれぞれ$M^{\langle 1)}\text{、}L^{(1)}$$\Phi$$M_{\text{、}^{}(2\}}L^{(2)}$計算されることが分かる。
(27)$M_{p}^{(1)}((_{0})= \int_{S}I_{p}(\zeta-\zeta_{0})uUdS_{\zeta}$$M_{\mathrm{p}}^{(2)}( \zeta_{0})=\int_{S}[I_{p}(\zeta-\zeta_{0})(\overline{\nu}U+I/\overline{U})-\overline{((-\zeta_{0})}I_{p-1}(\zeta-\zeta_{0})\nu U]dS_{\zeta}$

また、 多重極モーメントの展開中心の移動$(\mathrm{M}2\mathrm{M})$(28)は次式によって行われる。

(29)$M_{p}^{(1)}( \zeta_{1})=\sum_{q=0}^{p}M_{q}^{(1)}(\zeta_{0})I_{p-q}(\zeta_{0}-\zeta_{1})$$M_{p}^{(2)}( \zeta_{1})=\sum_{q=0}^{p}M_{q}^{(2)}(\zeta 0)I_{p-q}(\zeta_{0}-(_{1})-\overline{(\zeta_{0}-\zeta_{1})}M_{\}2}^{\dot{[}1)}-1((_{1})$多重極モーメントから局所展開係数へは次式によって変換される(30)$(\mathrm{M}2\mathrm{L})_{\text{。}}$$L_{q}^{(1)}(z_{0})=(-1)^{Q} \sum_{p=q}^{\infty}O_{p}(z_{0}-(_{0})M_{p-q}^{(1)}(\zeta_{0})$$L_{q}^{(2)}(z_{0})=(-1)^{q} \sum_{p=q}^{\infty}O_{p}(z_{0}-\zeta_{0})M_{p-q}^{(2)}(\zeta_{0})-\overline{(z_{0}-(_{0})}L_{q+1}^{(1)}(z_{0})$

(31)
(32)

152

局所展開係数の展開中心の移動

局所展開係数を用いて$(\mathrm{L}2\mathrm{L})$は次式によって行われる。

$L_{p}^{(1)}(z_{1})= \sum_{q=\mathrm{p}}^{\infty}L_{q}^{(1)}(z_{0})I_{q-p}(z_{1}-z_{0})$(33)$L_{p}^{(2)}(z_{1})= \sum_{q=\mathrm{p}}^{\infty}L_{q}^{(2)}(z_{0})I_{q-p}(z_{1}-z_{0})-\overline{(z_{1}-z_{0})}L_{p+1}^{(1)}(z_{1})$(34)$\Psi_{\text{、}}\Phi$は次のように計算される。
$\Psi=\sum_{\rho=0}^{\infty}L_{p}^{(1)}(z_{0})I_{\rho}(z-z_{0})$(35)$\Phi=\sum_{\rho=0}^{\infty}L_{p}^{(2)}(z_{0})I_{p}(z-z_{0})$(36)$T= \frac{\mu}{\pi(\kappa+1)}[n(\Psi’+\overline{\Psi’})-\overline{n}((z-z_{0})\overline{\Psi’’}-\overline{\Phi’})]$(37)

式 (26) のトラクションはにより求まる。

42 解の周期化
2 次元静弾性クラック周期境界値問題は次のように表される。
$\mu u_{i,jj}+(\lambda+\mu)uj,ji=0$

in

(38)

$\Omega\backslash S$

on

$S$

(39)

. $u_{i}(x^{2})$(40)$C_{ijkl}.u_{k},\iota(x^{1})nj(x^{1})=-C_{ijkl}u_{k,l}(x^{2})n_{j}(x^{2})$(41)$C_{ijkl}u_{k,l}n_{j}=g_{i}$$u_{\dot{\mathrm{z}}}(x^{1})--$

ここに、$\Omega$はユニットセル、$\Gamma$は$\Omega$の外側の境界、おける外向きの単位法線ベクトルである。
座標をとって、$\Omega$$x^{1}$と$x^{2}$

は原点を中心とする一辺の長さ力

はユニットセル中のクラックである。 また、$S$は$\backslash ^{\backslash }1$$\Gamma$$n_{i}$

は$\Gamma$に上で互いに反対側の点である。 簡単のために適当にの正方形領域であるとする。

次に、 周期境界条件を満たす解を得るため、 式 (24) を周期化する。 基本的なアイデアは、 ユニットセルと同じセルが周期的に無限個並んでいると考えるというものである。 しかし、 式 (24) 中に現れる各項のレプリ
カを単に加えると$\sum_{w\in N^{0}}\frac{1}{z-w}$(42)$\sum_{w\in N^{\mathrm{O}}}\frac{z-w}{(\overline{z}-\overline{w})^{2}}$(43)といった絶対収束しない級数が現れるため、 式の修正を要する。 ここに$N^{0}=${ $l=m+\mathrm{i}n|m,$$n$: 整数 }$N^{0}$

は、

(44)

153

に関して一次以下の項 (当然、なる複素数の集合である。 級数を絶対収束級数にするため、ある) を付け加えて和を取る。 \approx \Xi ‘{*的に 1 ま、 式 (42) (43) はそれぞれ次のように修正される。
$z_{\text{、}}\overline{z}$弾性学の解で(45)$\zeta(z)=\frac{1}{z}+\sum_{w\in N}(\frac{1}{z-w}+\frac{1}{w}+\frac{z}{w^{2}})$(46)$z \overline{\wp(z)}+\overline{\iota(z)}=\frac{z}{\overline{z}z}+\sum_{w\in N}(\frac{z-w}{(\overline{z}-\overline{w})^{2}}+\frac{w}{\overline{w}^{2}}-\frac{z}{\overline{w}^{2}}+\frac{2w\overline{z}}{\overline{w}^{3}})$

ここに、 $N$ は{ $l=m+\mathrm{i}n|m,$ : 整数$N=$である。 関数$\iota(z)$$n$$l\neq \mathrm{t}1$

(47)

}

は

(48)

$\iota(z)=\sum_{1\mathrm{L}\prime\in N}[\frac{-u}{(z-w)^{2}}\overline’+\frac{\overline{w}}{w^{2}}+\frac{2\overline{w}z}{w^{3}}]$と定義した。 は以下の性質を持つ。
$\iota$(49)$\iota(z+1)=\iota(z)-\wp(z)+\alpha$$(_{\iota}50)$$\iota(z+\mathrm{i})=\iota(z)+\mathrm{i}\wp(z)+\mathrm{i}\alpha$

ここに$\alpha$は$\alpha=4+\sum_{\uparrow r’\in N}[$$\frac{2(1-2w)}{(1-2\overline{w})^{2}}+\frac{2(1+2w)}{(1+2\overline{w})^{2}}-\frac{1}{\overline{w}^{2}}+\frac{2w}{\overline{w}^{3}}]$(51)\sim 50153097
なる定数である。 以上の$\zeta_{\text{、}}\iota$

を用いて解は以下のように$1_{\mathrm{D}}\exists$化される (\not\supset :的でなし ‘定数を除く)。

$V= \frac{1}{2\pi(/\sigma+1)}\int_{S}[\mathcal{K}\mathrm{I}/(s)U(s)\zeta(z-s)$$-(z-s)\overline{\nu(s)U(s)\zeta’(z-s)}$$+\overline{(\nu(s)\overline{U}(s)+\overline{\nu}(s)U(s))\zeta(z-s)}+\overline{U(s)\nu(s)\iota(z-.\mathrm{s})}$(52)$-\pi(\nu\overline{U}+\iota^{-}/U)z-(\pi h.\mathcal{U}(s)U(s)+\alpha\overline{\nu}(s)\overline{U}(6’))\overline{z}]dS_{\mathit{8}}$式 (52) の積分の多重極法による算法を次節で述べる。
t こより \Theta 境界値 5Q7
なお、 ユニットセルの外周 \iota c\rightarrow g 界ことも可能であり

rjb 式化 )$\not\in$$\text{、}$

要が生 b“*知数がaaegとして g#W7 する域での境界{直.
に沿って gF‘を E4- する必易であると考えられる。 しかし、 ユニットセル$\text{を取}$$7-7\text{題}$$\beta\S\rceil|\Re$$\sigma \mathit{3}\ovalbox{\tt\small REJECT}$$\overline{\mathrm{A}^{\text{、}}\prime}$\lambda - ることと、 後に示す$.\text{数}$ら、 本手法の方が有利であると考えられる。

{

R 結果からも分力 ‘ るように係数行列が悪条件であること力.

154

43 周期境界値問題における多重極法
周期境界値問題の多重極法は Laplace の場合とほぼ同じであるので、 ここでは、 レベル 0 の局所展開係数の最終結果を記すに留める。
$\sum\infty$$(-l)^{}$$o_{n}(w)M_{\mathfrak{n}.-\mathrm{p}}^{(1)}(0)$$L_{p}^{(1)}(0)= \sum_{w\in N’n=\max(2,\mathrm{p})}$(53)$+ \mathit{5}_{p1}\pi\frac{M_{0}^{(2)}(0)}{1-\kappa}$$(-1)^{n-p+1}O_{n}(w)M_{n-p}^{(2)}(0)|+$$\sum\infty$$L_{p}^{(2)}(0)= \sum_{w\in N’n=\max(2,p)}$$\sum_{w\in \mathrm{A}’’}\sum_{n=\max(2,p)}^{\infty}(-1)^{n-p+1}\overline{w}O_{n+1}(w)M_{n-p}^{(1)}(0)$(54)$-\delta_{p1}(\pi\kappa\overline{M_{0}^{(1)}(0)}+(\alpha-4)M_{0}^{(1)}(0))$

ここで、 式 (53)

式 (54) の評価には(55)$\sum_{v)\in N’}O_{n}(w)$(56)$\sum_{w\in N’}\overline{w}O_{\tau\iota+1}(w)$の値が $n\geq 2$ に対して必要となる。対称性から、 lattice sum が非零であるのは、式 (55) では-1) の時に限り、 式 (56) では$n$が$\zeta 4$$n$が (4 の倍数の倍数 +1) の時に限ることが分かる。 これらの lattice sum の値はあ

らかじめ数値計算により求めておく。

44 均質化法
以上の周期クラック問題の高速多重極法の応用として、 均質化法を取り上げる。 弾性定数 Q お’ を有する均質な線形弾性体を考える。 この物体は、 微視的なスケールにおいて多数のクラックを含み、 各クラックは表面力を受けないものとする。 この時、 この物質の巨視的な弾性定数は、 均質化法の理論により次のように求められる (詳細は Lions[ll] を参照されたい)。

(57)$C_{ijkl}- \frac{c_{ijpq}}{|\Omega|}\int_{S_{c}}\phi_{p}^{kl}n_{q}dy$

ここに、$\phi_{j}^{kl}$は、 クラック上の境界条件が、$c_{ijj}pq^{\frac{\partial\chi_{p}^{kl}}{\partial y_{q}}n=-C_{ijkl}n}j$で与えられる静弾性学の周期クラック境界値問題の解.
$\chi_{j}^{kl}$(on$S_{c}$)(58)の開口変位である。 以上より、 巨視的な弾性係数を評価するためには、周期境界条件のもとでクラックの開口変位$\phi_{j}^{kl}$を求めなければならないことが分かる。 この開口変位を求めるのに、 前節まで議論してきた多重極法を用いることができる。

155

図5

ユニットセル中のクラック配置

表1

反復回数

5 数値計算 (弾性体の場合)提案する周期境界値闇題の多重極法を用いて数値解析を行った。

5.1

直接法との解の比較

周期境界値問題における高速多重極法の解の妥当性を検証するために、 従来法によって求めた開口変位行った。従来法では、 ユニットセルの縁と高速多重極法によって求めた開口変位 (FMM) の比に要素を配置することによって周期境界条件を課した。 その際、 変位が定数の自由度を持つことを考慮しユ(Direct)$\mathrm{r}\mathrm{j}\mathrm{t}\text{を}$ニットセルの縁上の要素の つで変位を陽に与えた。
$-arrow$の正

適当な無次元化を行なってユニットセルは頂点が $(-0.5, -0.5)_{\text{、}}(-0.5,0.5)_{\text{、}}(0.5, -0.5)_{\text{、}}(0.5,0.5)$
数 l ま 200 である。 また。 クラック上の方形とし、半クラック長 0.4 のクラックを図 5 に示すように配$\text{要}\not\equiv_{\theta_{\text{、}^{}\mathrm{i}}}$$\text{置}$$\llcorner$

た、従来法におけるユニットセル縁の要素数は 1 辺あたり 300 とした。式 (39) における界条件) は$g_{i}=\delta_{2i}$(mode I) と得られた開口変 (1}‘‘ k 図$g_{i}=\delta_{1i}$(mode$\mathrm{I}\mathrm{I}$) の 2$1\text{、}\mathrm{f}\underline{\mathrm{f}\mathrm{i}}\text{り}\not\in:\mathrm{f}R$に示す。 両者の 課はほぼ一致しており、 . p’r-\rightarrow ffl 文の手法の妥当性が示されたと言$\text{本}$$\#_{\backslash }$$6_{\text{、}}7$

離散化して得られる連立一次
素の相互作用 g
分かるように(クラック上の境

つた。 Lam\’e 定数は $\lambda=\mu=1$ とした。

えるとした。 この時の反$g_{i}$$\ovalbox{\tt\small REJECT}/\fbox_{\mathfrak{o}}$数$\text{を表}$$\text{方程}$式はリスター$\text{ト}\Phi\tau\backslash \text{、}$

しの GMRES によって解1 に示す。 なお、 多重極法の場合の前$\beta\Delta \mathrm{f}\mathrm{E}$$\iota\backslash$

た。 M

半\not\in-

準は $10^{-5}$はリーフ内 (最大要素数 20) の要す行列を前処理行列とする右前処理である。従来法では前処理を行っていない。表 1 力 ら$\text{多重極法の方}\mathrm{B}^{\backslash }\backslash \acute{\{}j\backslash \xi \text{来法よりも反}\mathfrak{S}\fbox_{\circ}$

られることが分かる。

$\backslash$

数が少なく、 前処理を施すことによりさら$\mathrm{t}_{\check{\mathrm{L}}}$g 好な収束が得158

図6

mode I の場合の

図 7

mode

$\mathrm{I}\mathrm{I}$$x_{2}$

方向開口変位の場合の $x1$ 方向開口変位

5.2 均質化法への応用
周期境界値問題における高速多重極法を用いて、 多数のクラックを含む物体の巨視的な弾性定数を計算する。

いま、微視的構造が図 5 と同様な、$x_{1}$

軸に平行な単一の直線クラックからなり、半クラック長だけが 0.1 となった場合を考える。 ただし、 物性値は $\lambda=\mu=1$ とする。 すると、 式 (57)(58) により巨視的な弾性定数は表 2 のように求められる。 これをコンプライアンス $D_{ijkl}$ に換算すると、 非 0 成分のうちクラックの影響を受けるものは D2222 $=0.3990_{\text{、}}D_{1212}=\mathrm{L}024$ となる。 これは、 Horii and Sahasakmontri[12] による近似
解 (クラック長/ クラック聞隔が小さい時に有効) $D_{2222}=0.3989_{\text{、}}D_{1212}=1.024$ と良く$–\wedge$

致している。 なお、 Horii and Sahasakmontri の $D_{2222}$ の値はグラフを読み取って求めたので、 4 桁目は概数である。

次に、 微視構造のクラック配置が図 8 であるような物体を考える。 半クラック長 005 のクラックが中心閤
隔 0125 で 8 $\mathrm{x}8=64$ 本並んでいる。 クラックー本当たりの要素数は 200 であるので総自由度数は 25600 である。 クラックの向きはランダムに与えた。 この場合、 巨視的な Cijk\beta よ、 表 3 のように求まった。

157

表3

巨視的な弾性定数(random)$C_{ijkl}$

6 結論
本稿では 2 次元 Laplace 及び静式化を$\mathfrak{M}_{\overline{\mathrm{D}}}^{-}\dashv$$\overline{\tau_{\text{、}}\prime}\text{し}$、$\ovalbox{\tt\small REJECT}^{\backslash }\ovalbox{\tt\small REJECT}^{\backslash }’|*$

クラック

によってその妥当$5\mathrm{f}\mathrm{f}\mathrm{i}^{\mathrm{B}}\ovalbox{\tt\small REJECT}$における周$l\mathrm{E}^{i}\mathrm{E}\mathrm{f}\mathrm{f}\mathrm{i}_{\vec{\overline{0}}}^{=}\mathrm{f}$し為 また、–$\backslash$—$\backslash$$0.2-/\sim$$/$$\backslash _{\backslash }$$|$$\backslash$$/$$\xi\Re\Psi \mathrm{g}$$\backslash$

$|$

界値

$\mathfrak{B}’$$7_{\mathrm{R}}5^{\mathrm{a}}\ovalbox{\tt\small REJECT}$

$/$

$/$

$/$

$\backslash$$\backslash \backslash /$$/\backslash /--$

$|$

$/$

$\backslash /$0$\backslash \sim\backslash$$.02//$$\backslash$$\backslash /$

$|$

$|$

$\backslash -/$$\backslash$$\backslash ---/\backslash -$$\backslash$

$/$

/
$/$

$0A$

$/_{04}.\sim/0.2$

図 8

$/0/$

.

$\text{多重}$ $\varpi_{-\mathrm{f}\mathrm{l}}^{\mathrm{p}\text{界要}}$素法

化\not\in ‘k 用し ‘て多数のクラックを含む弾

04$\backslash /$に対する高速$\backslash 0.2\backslash /04$

クラック配置 (ランダム配置)

158

性体の巨視的な弾性定数を求めた。
今後の課題としては、 インクルージョン問題における定式化、 あるいは 3 次元の静弾性問題に対する定式化などが挙げられる。

参考文献[1] N. Nishimura (2002), Fast multipole accelerated boundary integral equation methods, Applied Mechanics Reviews, 55, PP.299-324. Chen and H. Munakata (2005), A fast boundary. Liu, N. Nishimura, Y. Otani, T. Takahashi,[2]element method for the analysis of fiber-reinforced composites based on a rigid-inclusion model,$\mathrm{X}.\mathrm{L}$$\mathrm{Y}.\mathrm{J}$

Journal of Applied Mechanics, 72, PP 115-128[3] 寺田賢二郎・菊池昇 (2003), 均質化法入門, 丸善[4]$\mathrm{h}\mathrm{t}\mathrm{t}\mathrm{p}://\mathrm{w}\mathrm{w}\mathrm{w}.\mathrm{c}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{l}.\mathrm{t}\mathrm{o}\mathrm{h}\mathrm{o}\mathrm{k}\mathrm{u}.\mathrm{a}\mathrm{c}.\mathrm{j}\mathrm{p}/$.htm$\mathrm{t}\mathrm{e}\mathrm{i}/\mathrm{e}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{h}/\mathrm{h}\mathrm{o}\mathrm{m}\mathrm{o}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{i}\mathrm{z}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$Koro and K. Abe(2002), A wavelet method for reducing the computational cost of BE-based
homogenization analysis, Engineering Analysis with Boundary Elements, 27, pp.439-454[6] L. Greengard and V. Rokhlin(1987), A Fast Algorithm for Particle Simulations, Journal of Compu[5]$\mathrm{K}$tational Physics, 73, Pp.325-348[7] L. Greengard and J. Helsing (1998), On the numerical evaluation of elastostatic fields in localty
1441-1462.
Mech. Phys. Solids, 46,
isotropic two-dimensional composites,$\mathrm{J}$$\mathrm{p}\mathrm{p}$[8] 松村知樹, 西村直志 (2004), クラックの周期境界値問題における高速多重極法と均質化法への応用について, 計算数理工学論文集,、 4pp.95-100[9] 大谷佳広, 西村直志 (2004), 2 次元静弾性クラック周期境界値問題における高速多重坐法, 境界要素法論文
集, 21, pp.71-76. Stegun (1972), Handbook of Mathematical Functions, Dover, New York;
M. Abram owitz andChapter 18 Weierstrass Elliptic and Related Functions. Lions(1980), Remarks on Som Asymptotic Problems in Composite and in Perforated Materials,
lrll]
Variational Methods in the Mechanics of Solids, pp.3-20[12] H. Horii and K. Sahasakmontri (1990), Mechanical properties of eraeked solids: validity of theAbe),
Taya and. Weng,self-consistent method, Micromechanics and Inhomogeneity (Eds.
[

$10|$

$\mathrm{I}.\mathrm{A}$$\cdot$$\mathrm{J}.\mathrm{L}$$\mathrm{e}$$\mathrm{G}.\mathrm{J}$Springer, 137-159.

$\mathrm{M}$$\mathrm{H}$