今回は，境界形状に摂動が加わったときの解析方法を紹介する。例として，一様流中に置かれた気泡の変形問題を考える。簡単のため流れをポテンシャル流れとし，微小変形する場合を扱う。なお，流体力学を勉強すると，ポテンシャル流れは，非常に特殊な流れで，現実に存在するほとんどの流れは，ポテンシャル流れとは大きく異なる速度分布を持つことを知る。しかし，球形気泡および球形気泡からわずかに変形した気泡周りの流れは，固体粒子と異なり気泡表面での渦度の生成量が少ないためにポテンシャル流れが良い近似となる。実際，非線形のNavier-Stokes方程式を数値計算により解いた球形気泡周りの速度分布は，レイノルズ数が高くなるとポテンシャル流れの速度分布に漸近するという特異な性質を持つ。

領域摂動法（Domain Perturbation Method）

領域摂動法とは，ある形状に対して解が与えられるとき，その形状からわずかにずれた形状に対しての近似解を求めるための方法である。ラプラス方程式やポアソン方程式などの境界値問題に対して，解析解が与えられる単純な形状からずれた場合の解を計算するのに用いられることが多い。

ここでは，球体周りのポテンシャル流れの解を利用して，球形からわずかに変形した物体周りのポテンシャル流れの解を求める方法を説明する。ポテンシャル流れとは，流れの中に渦度の存在しない流れで，このとき流れの速度ベクトル${\bf u}$は，スカラー関数$\varphi$の勾配を用いて，

$${\bf u} = \nabla \varphi$$

(80)

で与えられる。この$\varphi$を速度ポテンシャルと呼び，さらに流れが非圧縮であるときには，非圧縮流れの連続の式（$\nabla \cdot {\bf u} = 0$）を用いると，スカラー関数$\varphi$は，ラプラス方程式

$$\nabla^2 \varphi = 0$$

(81)

を満たす。このラプラス方程式の解は，与えられた境界条件のもと，変数分離法あるいは，本連載記事の前半で説明した微分演算子に対する固有値問題として解く方法により解析解を得ることができる。ここでは，まずは一様流中に存在する球形物体周りの解を取得し，その結果を利用して変形物体の解を求めるため，球形物体周りの解を求めるのに適した球座標系による解析解を取得する。詳細を参考資料(5)に示すが，球座標系$(r,\theta,\phi)$による式(81)の解は，一様流の方向を球座標の軸方向（$\theta = 0,\pi$の方向）と一致させると，球体周りの流れは軸対称になり，直交関数系であるルジャンドル多項式（$P_n(\cos \theta)$）を用いて，次式で与えられる。

$$\varphi (r,\theta) = \sum_{k=0}^\infty \left( A_k r^k + \frac{B_k}{r^{k+1}} \right) P_k(\cos \theta)$$

(82)

具体的に書くと，$P_1(\cos \theta) = \cos \theta$，$P_2(\cos \theta) = \dfrac{1}{2}(3\cos^2\theta – 1)$，$P_3(\cos \theta) = \dfrac{1}{2}(5\cos^3 \theta – 3\cos \theta) = \dfrac{5}{8} \cos 3 \theta + \dfrac{3}{8} \cos \theta$，となる。

境界条件は，遠方と粒子表面で与えられ，以下の様になる。

i）遠方では，一様流$U_\infty$より，無限遠で

$${\bf u} = U_\infty {\bf e}_x$$

(83)

ここで，${\bf e}_x$は，一様流の方向の単位ベクトル。

ii）物体表面（$r = R_0$）では，表面を貫く流れがない条件より

$${\bf u} \cdot {\bf n} = 0$$

(84)

ここで，${\bf n}$は物体表面の法線ベクトル。

速度ベクトルは，

$${\bf u} = \nabla \varphi = \frac{\partial \varphi}{\partial r} {\bf e}_r + \frac{1}{r} \frac{\partial \varphi}{\partial \theta} {\bf e}_\theta$$

(85)

で与えられるので，式(83)，(84)の境界条件を用いると，

$$\text{無限遠で，}\quad\qquad \varphi = U_\infty r \cos \theta$$

(86)

$$\text{物体表面で，}\qquad u_r = \left. \frac{\partial \varphi}{\partial r} \right|_{r = R_0} = 0$$

(87)

となる。この境界条件を用いて，式(82)の係数を決めると，一様流中に存在する球形物体に対するポテンシャル流れの解として次の解を得る。

$$\varphi = U_\infty \left( r + \frac{R_0^3}{2r^2} \right) P_1(\cos \theta) = U_\infty \left( r + \frac{R_0^3}{2r^2} \right) \cos \theta$$

(88)

これで，摂動する前のリーディングオーダーの解が得られたので，この解を用いて形が球形からずれたときの解を求める。

さて，ここで重要となるのは，非球形の変形をどのように表現するかである。注目する点は，ラプラス方程式の解が，式(82)の様に$\theta$方向に関して，ルジャンドル多項式で与えられることである。したがって，変形に関してもルジャンドル多項式を用いて変形モードを表現すると，以下に見るようにルジャンドル多項式の直交関係を用いて高次の解を求めて行くことができる。すなわち，球形からのずれを表すのに，球の中心から物体表面までの距離を$r$とおき，物体表面の形状をルジャンドル多項式（$P_n(\cos \theta)$）を用いて，

$$r = R_0(1 + \varepsilon_n P_n (\cos \theta)) \qquad (\varepsilon_n \ll 1)$$

(89)

と記述する。ここで与えた変形は，方位角$\phi$の成分を含んでおらず，軸対称な変形モードを表している。実際に$n = 2 {\sim} 4$の変形モードを，球形モードに重ねたものを図1に示す。$n = 2$が回転楕円体モードで，$n$が増えるに従って，複雑な変形モードとなっているのがわかる。

図1　変形モードの様子（$n = 2 {\sim} 4$）

さて，これらの変形モードに対する近似解を，式(88)で与えられる球形周りの解をもとに求めていく。この様に形状に対して摂動を与え，近似解を求めていく手法を領域摂動法と呼ぶ。式(88)で与えられる$\varphi$を$\varphi_0$と置き，微小変形量$\varepsilon_n$を用いて，求める速度ポテンシャルの解を次のように表現する。

$$\varphi \equiv \varphi_0 + \varepsilon_n \varphi_1 + O(\varepsilon_n^2)$$

(90)

したがって，速度ポテンシャルのラプラス方程式は，

$$\nabla^2 \varphi \cong \nabla^2 \varphi_0 + \varepsilon_n \nabla^2 \varphi_1 = 0$$

(91)

と近似でき，$\varepsilon_n \ll 1$の下でスケールの分離が可能となり，

$$O(1) \text{の解として，} \nabla^2 \varphi_0 = 0$$

(92)

$$O(\varepsilon_n) \text{の解として，} \nabla^2 \varphi_1 = 0$$

(93)

を得る。

物体表面では，${\bf u} \cdot {\bf n} = 0$の境界条件を満たすが，この式を用いるためには，非球形の物体表面における法線ベクトル${\bf n}$を求める必要がある。物体表面は，$r = R_0(1 + \varepsilon_n P_n(\cos \theta ))$と表わされることより，

$$F(r,\theta ) = r – R_0(1 + \varepsilon_n P_n(\cos \theta))$$

(94)

とすると，物体表面は，$F(r,\theta) = 0$と表わされる。したがって，物体表面における法線ベクトル${\bf n}$は，

$$\begin{split} {\bf n} &{} = \left. \frac{\nabla F}{\left| \nabla F \right|} \right|_{F=0} = \left. \frac{\left( \dfrac{\partial F}{\partial r}, \dfrac{1}{r} \dfrac{\partial F}{\partial \theta} \right)}{\sqrt{\left( \dfrac{\partial F}{\partial r} \right)^2 + \left( \dfrac{1}{r} \dfrac{\partial F}{\partial \theta} \right)^2}} \right|_{F=0} \\ &{} = {\bf e}_r – \varepsilon_n \frac{\partial P_n(\cos \theta)}{\partial \theta} {\bf e}_\theta + O(\varepsilon_n^2) \end{split}$$

(95)

と近似することができる。したがって，物体表面で，

$${\bf u} \cdot {\bf n} \cong \left. \left[ \frac{\partial \varphi_0}{\partial r} + \varepsilon_n \frac{\partial \varphi_1}{\partial r} – \varepsilon_n \frac{1}{r} \frac{\partial \varphi_0}{\partial \theta} \frac{{\rm d} P_n(\cos \theta)}{{\rm d} \theta} \right] \right|_{F=0} = 0$$

(96)

さらに，

$$\begin{split} \left. \frac{\partial \varphi_0}{\partial r} \right|_{F=0} &{} = \left. \frac{\partial \varphi_0}{\partial r} \right|_{r = R_0(1 + \varepsilon_n P_n(\cos \theta))} \\ &{} = \left. \frac{\partial \varphi_0}{\partial r} \right|_{r = R_0} + \varepsilon_n R_0 P_n \left. \frac{\partial^2 \varphi_0}{\partial r^2} \right|_{r = R_0} + O(\varepsilon_n^2) \end{split}$$

(97)

の関係を用いて，式(96)を書き直すと，$O(\varepsilon_n)$の境界条件として次式を得る。

$$R_0 P_n \left. \frac{\partial^2 \varphi_0}{\partial r^2} \right|_{r = R_0} + \left. \frac{\partial \varphi_1}{\partial r} \right|_{r = R_0} – \left. \frac{1}{R_0} \frac{\partial \varphi_0}{\partial \theta} \right|_{r = R_0} \frac{{\rm d} P_n}{{\rm d} \theta} = 0$$

(98)

式(88)より$\varphi_0 = U_\infty \left( r + \dfrac{R_0^3}{2r^2} \right) \cos \theta$となるので，式(98)は，$\varphi_1$に対する境界条件になっている。この式より，境界条件を介して，$\varphi_0$が$\varphi_1$に影響を与えるのがわかる。

式(93)より，$O(\varepsilon_n)$の解は，ラプラス方程式を満たし，かつ無限遠で，$\varphi_1 = 0$を満たすことより，

$$\varphi_1 = \sum_{k=0}^\infty \frac{B_k}{r^{k+1}} P_k(\cos \theta)$$

(99)

$$\therefore \quad \left. \frac{\partial \varphi_1}{\partial r} \right|_{r = R_0} = \sum_{k=0}^\infty \frac{-(k + 1) B_k}{R_0^{k+2}} P_k(\cos \theta)$$

(100)

と与えられ，さらに物体表面で，式(98)を満たすことより

$$\begin{split} \left. \frac{\partial \varphi_1}{\partial r} \right|_{r = R_0} &{} = \left. \left( -R_0 \frac{\partial^2 \varphi_0}{\partial r^2} P_n + \frac{1}{R_0} \frac{\partial \varphi_0}{\partial \theta} \frac{{\rm d} P_n}{{\rm d} \theta} \right) \right|_{r = R_0} \\ &{} = -3U_\infty \cos \theta P_n(\cos \theta) – \frac{3}{2} U_\infty \sin \theta \frac{{\rm d} P_n}{{\rm d} \theta} \end{split}$$

(101)

すなわち，式(100)，(101)より，次式を得る。

$$\sum_{k=0}^\infty \frac{-(k + 1) B_k}{R_0^{k+2}} P_k(\cos \theta) = -3U_\infty \cos \theta P_n(\cos \theta) – \frac{3}{2} U_\infty \sin \theta \frac{{\rm d} P_n}{{\rm d} \theta}$$

(102)

詳細は参考資料(5)に記載するが，この式にルジャンドル多項式の直交関係式，

$$\int_0^\pi P_n(\cos \theta) P_m(\cos \theta) \sin \theta \, d\theta = \frac{2}{2n + 1} \delta_{n\,m}$$

(103)

および，ルジャンドル多項式が満たす漸化式（資料(5)参照）を用いて，整理すると，次式を得る。

$$\begin{split} \varphi_1 &{} = \left[ \frac{B_{n-1}}{r^n} P_{n-1} (\cos \theta) + \frac{B_{n+1}}{r^{n+2}} P_{n+1} (\cos \theta) \right] \\ &{} = \frac{3}{2} U_\infty \left[ -\frac{R_0^{n+1}}{r^n} \frac{n – 1}{2n + 1} P_{n – 1}(\cos \theta) + \frac{R_0^{n+3}}{r^{n+2}} \frac{n + 1}{2n + 1} P_{n + 1}(\cos \theta) \right] \end{split}$$

(104)

すなわち，球形のときの解に対して，図1に示すような$n$次の形状モードと，式(88)で与えられる速度ポテンシャルの$P_1(\cos \theta)$成分との相互作用により，式(90)の$\varphi_1$の項には，速度ポテンシャルの$n – 1$と$n + 1$の項が現れることになる。具体的に$n = 2$のときを与えておくと，

$$\varepsilon_2 \varphi_1 = \varepsilon_2 \frac{3}{2} U_\infty \left[ -\frac{1}{5} \frac{R_0^3}{r^2} P_1(\cos \theta) + \frac{3}{5} \frac{R_0^5}{r^4} P_3(\cos \theta) \right]$$

(105)

と，一様流である$P_1$の流れと，回転楕円体形状を表現する$P_2$の形状モードの相互作用で，$P_1$と$P_3$の流れが生じるのがわかる。

ここで示した領域摂動法は，解析解が得られる問題からスタートして，境界形状がそれからずれたときの解を摂動法により求める方法であるが，同様の考え方は，数値計算でも利用できる。たとえば，球形周りの数値解を球座標系で計算した後に，変形の効果を直交関数系の級数により表現し，各変形モードに対する計算を実施することも可能である。そのような手法は，数値解を直交関数系で表現するスペクトル法において特に有効である。数値計算の場合には，境界形状を微小変形とする必要はなく，式(89)に示した様な直交関数による形状表現を，大変形に対しても適用することができる。

 

参考

(5) https://www.jsme.or.jp/jsme/uploads/2018/11/KikaiyaMath1811.pdf

---

＜フェロー＞

高木 周

◎東京大学・大学院工学系研究科・機械工学専攻 教授

---