はじめに
本記事では、境界要素法（BEM: Boundary Element Method）を用いたポテンシャル流れの数値解析について、理論的な背景から実装までを解説します。
流体解析において、ポテンシャル流れは最も基本的かつ重要な概念の一つです。これは非圧縮・非粘性・非回転という理想化された条件下での流れを指します。実際の流体現象を完全に表現することはできませんが、航空力学や水理学をはじめとする多くの工学分野で有用な近似解を提供します。特に、物体周りの流れや複雑な形状を持つ領域内の流れを解析する際、ポテンシャル流れの理論は強力なツールとなります。
BEMは境界のみを離散化することで計算コストを大幅に削減できる数値解析手法です。有限要素法（FEM）や有限差分法（FDM）といった手法では領域全体をメッシュ化する必要がありますが、BEMでは解析対象の次元を1つ減らすことができます。例えば、2次元問題では境界が1次元となるため、計算量が大幅に削減されます。この特性は、ポテンシャル流れのように支配方程式がラプラス方程式となる問題に対して特に有効です。
本記事では、迷路のような複雑な環境において、入口から出口へ至る流れ場をBEMによって計算する手法を扱います。以下の図は、迷路状領域におけるポテンシャル流れの例であり、BEMを用いて内部領域におけるポテンシャル値と速度ベクトルを計算した結果です。
図: 迷路状領域におけるポテンシャル流れの例
本記事は以下のように構成されています。本記事が、BEMやポテンシャル流れ解析に興味を持つ読者の皆様にとって、理論と実装の両面から理解を深める一助となれば幸いです。
- 第1章：シミュレーションの目的と解析フロー、ポテンシャル流れとBEMの基礎概念
- 第2章：ポテンシャル流れの支配方程式から境界積分方程式を導出する定式化
- 第3章：境界積分方程式の離散化と数値計算の具体的手法
本記事で解説する手法の実装コードは、以下のGitHubリポジトリで公開しています。
目次
1. シミュレーションの概要
1.1 本シミュレーションの目的と解析フロー
本シミュレーションでは、障害物が存在する複雑な環境において、始点から終点へ至る流れ場を計算します。具体的には、迷路のような環境における入口から出口への流れを、ポテンシャル流れと境界要素法（BEM: Boundary Element Method）を用いて解析します。
解析は以下の4つの手順で進めます。
- 形状定義：解析領域の境界（壁、入口、出口）を閉曲線として定義し、線分要素に分割します
- 境界条件の設定：入口と出口にはポテンシャル値を与え、壁面には法線流速ゼロ（壁を突き抜けない）という条件を設定します
- 境界値の決定：BEMを用いて境界積分方程式を解き、境界上の未知量（ポテンシャルまたは流速）を決定します
- 内部流速の計算：得られた境界値を用いて、領域内部の任意点における速度ベクトルを計算します
1.2 ポテンシャル流れとは
ポテンシャル流れとは、流体の速度場 $\boldsymbol{v}$ がスカラー関数である速度ポテンシャル $\phi$ の勾配として表される流れです。数式で表すと、
$$\boldsymbol{v} = \nabla \phi \tag{1.1}$$
を満たす流れを指します。
ポテンシャル流れが成立するためには、以下の3つの条件を満たす必要があります。
- 非回転性（渦なし）：流体の回転成分である渦度 $\boldsymbol{\omega} = \nabla \times \boldsymbol{v}$ がゼロであること。渦なし流れであれば、必ず速度ポテンシャルが存在し、ポテンシャル流れとなります
- 非圧縮性：流体の密度 $\rho$ が一定であること。連続の式 $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \boldsymbol{v}) = 0$ より、非圧縮性の場合は $\nabla \cdot \boldsymbol{v} = 0$ が成り立ちます
- 非粘性：流体の粘性が無視できること
また、本シミュレーションでは、さらに以下の条件も仮定します。
- 定常性：時間に依存しない流れであること ($\frac{\partial \boldsymbol{v}}{\partial t} = 0$)
- 外力の不存在：重力などの外力が存在しないこと
注：条件3の非粘性は、式 $(1.1)$ を満たすための必要条件ではありません。しかしながら、流体が存在する領域全体で非回転性を保つためには必要となります。詳細については補遺B.1をご参照ください。
ポテンシャル流れには、計算コストが低く理論的な解析が容易であるという利点があります。一方で、粘性の影響が支配的な壁面近傍の流れや、渦が発生する流れなど、多くの実際の流体現象を正確に表現することはできません。それでも、航空力学や水理学などの分野では、ポテンシャル流れを用いた近似解析が広く利用されています。特に、次節で説明する境界要素法（BEM）を用いた数値解析が盛んに行われています。
1.3 BEMとは
境界要素法（BEM: Boundary Element Method）は、境界上の情報のみを用いて領域内部の物理量を計算する数値解析手法です。
有限要素法（FEM）や有限差分法（FDM）といった他の主要な数値解析手法では、領域全体を離散化（メッシュ化）する必要があります。これに対して、BEMは境界のみを離散化する点に大きな特徴があります。
図 1.1 : FEM, FDM, BEMの比較（イメージ図）
BEMの主な利点は以下の通りです。
- 計算コストの削減：解析対象の次元を1つ減らすことができます。例えば、2次元領域の解析では境界が1次元となるため、計算コストが大幅に削減されます。
- 無限領域への適用：領域が無限に広がる問題であっても、境界上の情報のみを扱うため、容易に解析できます。これは、無限遠まで広がる流れ場の解析に特に適しています。
一方、BEMの欠点として、離散化によって得られる連立方程式の係数行列が密行列かつ非対称となるため、要素数が増加すると計算コストが急激に増大するというものがあります。この問題は、高速多重極法（FMM: Fast Multipole Method）などの手法を組み合わせることで、ある程度緩和することができます。本記事では、基本的なBEMの理論と実装に焦点を当てるため、これらの高度な手法については扱いません。
2. 境界積分方程式の定式化
2.1 支配方程式
本章では、前節で紹介したBEMを用いて、領域内部の任意の点における速度ベクトルを計算する方法について説明します。まず、ポテンシャル流れの支配方程式を確認しましょう。
節1.2で述べたように、ポテンシャル流れでは、点 $\boldsymbol{x}$ における速度ベクトル $\boldsymbol{v} = (u, v)$ が、スカラー関数である速度ポテンシャル $\phi(\boldsymbol{x})$ を用いて式 $(1.1)$ のように表されます。さらに、非圧縮性の条件 $\nabla \cdot \boldsymbol{v} = 0$ を適用すると、速度ポテンシャル $\phi$ は以下のラプラス方程式を満たすことがわかります。
$$\nabla^2 \phi = \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} = 0 \tag{2.1}$$
したがって、ポテンシャル流れの解析は、与えられた領域内でラプラス方程式 $(2.1)$ と境界条件を満たす速度ポテンシャル $\phi$ を求める問題に帰着されます。
2.2 境界条件の設定
次に、閉曲線境界 $\Gamma$ 上での条件を設定します。閉曲線によって区切られた、流体が存在する側を内部、その補集合を外部として、境界上の外向き法線ベクトルを $\boldsymbol{n}$ とします。
境界は、以下の2種類に分類されます。
- 開口部：流体の出入りがある部分
- 壁面：流体が壁を貫通しない部分
開口部の境界条件
開口部では、法線方向の速度成分を既知の値として指定します。図 2.1 (a), (b) に示すように、流体領域へ流入する部分では $V_{\text{in/out}} < 0$、流体領域から流出する部分では $V_{\text{in/out}} > 0$ となります。
$$\frac{\partial \phi}{\partial n} = \boldsymbol{v} \cdot \boldsymbol{n} = V_{\text{in/out}} \tag{2.2}$$
ここで、$V_{\text{in/out}}$ は既知の値であり、法線の向きの定義に依存して流入なら負、流出なら正となります。
壁面の境界条件
壁面では、流体が壁を突き抜けないため、法線方向の速度成分はゼロになります。
$$\frac{\partial \phi}{\partial n} = \boldsymbol{v} \cdot \boldsymbol{n} = 0 \tag{2.3}$$
このとき、図 2.1 (c) に示すように、流れは壁に沿って「滑る」形となります。節1.2で述べたように、ポテンシャル流れでは粘性の影響が無視されるため、壁面での滑りなし条件（壁面において速度ベクトルがゼロとなる条件）は課されません。
図 2.1 : 点 $\boldsymbol{\xi}$ における境界条件と流れのイメージ
質量保存則による制約
流体では連続の式（質量保存則）により、閉じた領域への総流入量と総流出量は一致しなければなりません。ポテンシャル流れでは流体は非圧縮であるため、全境界での流入体積流量と流出体積流量の和はゼロとなります。
$$\int_{\Gamma} \frac{\partial \phi}{\partial n} ds = 0 \tag{2.4}$$
内部境界がある場合
内部境界 $\Gamma_{\text{inner}}$（障害物など）を設ける場合、その法線ベクトルはその境界が囲う領域の内側を向くように定義します。このとき、流体領域は外部境界 $\Gamma_{\text{outer}}$ に内包され、内部境界の外側に位置する領域となります。
図 2.2 : 内部境界と外部境界における法線ベクトルの向き
プログラム上での境界条件の設定について
本プログラムでは、以下の2種類の境界条件を使い分けます。
- 開口部：速度ポテンシャル $\phi$ の値を指定（ディリクレ条件）
- 壁面：法線速度成分 $\partial \phi / \partial n$ を指定（ノイマン条件）
開口部における $V_{\text{in/out}}$ は、指定されたポテンシャル値から計算されるため、プログラム上では明示的に指定する必要はありません。
設定例：$(0, 0)-(1, 0)-(1, 1)-(0, 1)$ の四角形領域において、左辺を入口、右辺を出口、上辺と下辺を壁面とする場合
- 入口（左辺）: $\phi = 0$ （ディリクレ条件）
- 出口（右辺）: $\phi = 1$ （ディリクレ条件）
- 上辺および下辺（壁面）: $\frac{\partial \phi}{\partial n} = 0$ （ノイマン条件）
ポテンシャルの値 $\phi$ は任意の値を設定可能ですが、一般には入口を $0$、出口を正の値に設定することが多いです。
なお、1つでもディリクレ条件を設定すれば式 $(2.4)$ は自動的に満たされるため、特に考慮する必要はありません。内部境界を全て壁面として扱う場合も同様です。すべてノイマン条件で境界条件を設定する場合のみ、各線分における（長さ × 指定した値）の和がゼロとなるように注意する必要があります。
2.3 境界積分方程式による解析
任意形状の境界を持つ問題において、ラプラス方程式の解を求めるために、グリーンの第2恒等式を利用します。まず、必要な記号を定義しましょう。
- 観測点 $\boldsymbol{x}$：ポテンシャルを求めたい点
- 積分点 $\boldsymbol{\xi}$：境界 $\Gamma$ 上の点
- 相対位置ベクトル $\boldsymbol{r} = \boldsymbol{\xi} - \boldsymbol{x}$、その大きさ $r = |\boldsymbol{r}|$
基本解の導入
点 $\boldsymbol{x}$ に単位強度の点源を置いたときに生じるポテンシャル場である、2次元ラプラス方程式の基本解は次式で与えられます。
$$G(\boldsymbol{x}, \boldsymbol{\xi}) = -\frac{1}{2\pi} \ln r \tag{2.5}$$
境界積分方程式
グリーンの第2恒等式を適用すると（参考文献[1]ではGreen-Gaussの定理から導出）、領域内部の点 $\boldsymbol{x}$ における速度ポテンシャル $\phi(\boldsymbol{x})$ は、境界 $\Gamma$ 上の積分で次のように表されます。
$$c(\boldsymbol{x}) \phi(\boldsymbol{x}) = \oint_{\Gamma} \left(G(\boldsymbol{x}, \boldsymbol{\xi}) \frac{\partial \phi(\boldsymbol{\xi})}{\partial n} - \phi(\boldsymbol{\xi}) \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n} \right) ds(\boldsymbol{\xi}) \tag{2.6}$$
ここで、$\frac{\partial \phi(\boldsymbol{\xi})}{\partial n}$ はポテンシャルの法線方向微分（法線方向流速）です。また、係数 $c(\boldsymbol{x})$ は、観測点の位置に依存する幾何学的係数であり、式 $(2.7)$ のように定義されます[1, Eq. (29)]。
c(\boldsymbol{x}) = \begin{cases}
1 & (\boldsymbol{x} \text{ が領域内部にある場合}) \\
\frac{1}{2} & (\boldsymbol{x} \in \Gamma \text{ かつ境界が滑らかな場合}) \\
1 - \frac{\theta}{2\pi} & (\boldsymbol{x} \in \Gamma \text{ かつ境界が角を持つ場合}) \\
0 & (\boldsymbol{x} \text{ が領域外にある場合})
\end{cases} \tag{2.7}
本シミュレーションでは、境界を線分要素で近似しているため、境界上の点 $\boldsymbol{x} \in \Gamma$ では $c(\boldsymbol{x}) = 1/2$ とします。
境界条件との関係
式 $(2.6)$ の右辺には、境界上の2つの量が現れます。
- $\boldsymbol{q}(\boldsymbol{\xi}) = \frac{\partial \phi(\boldsymbol{\xi})}{\partial n}$：ノイマン条件が指定された場合は既知、それ以外では未知量
- $\phi(\boldsymbol{\xi})$：ディリクレ条件が指定された場合は既知、それ以外では未知量
したがって、式 $(2.6)$ は境界上の未知量 $\boldsymbol{q}(\boldsymbol{\xi}), \phi(\boldsymbol{\xi})$ に関する積分方程式となります。この積分方程式を離散化して連立方程式として解くことで、境界上の全ての点における $\boldsymbol{q}(\boldsymbol{\xi}), \phi(\boldsymbol{\xi})$ を求めることができます。
内部点での速度ベクトルの計算
境界上の $\phi$ と $\frac{\partial \phi}{\partial n}$ が求まれば、内部点 $\boldsymbol{x}$ における速度ベクトル $\boldsymbol{v}(\boldsymbol{x}) = \nabla \phi(\boldsymbol{x})$ を計算できます。内部点 $\boldsymbol{x}$ では $c(\boldsymbol{x}) = 1$ であることから、式 $(2.6)$ の両辺を $\boldsymbol{x}$ で微分すると、以下の式 $(2.8)$ が得られます。ここで、最後の式は、境界を $N$ 個の要素 $\Gamma_j$ に分割した場合の離散化表現です。
\begin{aligned}
\boldsymbol{v}(\boldsymbol{x}) &= \nabla_{\boldsymbol{x}} \phi(\boldsymbol{x}) \\\
&= \oint_{\Gamma} \left( \nabla_{\boldsymbol{x}} G(\boldsymbol{x}, \boldsymbol{\xi}) \frac{\partial \phi(\boldsymbol{\xi})}{\partial n} - \phi(\boldsymbol{\xi}) \nabla_{\boldsymbol{x}} \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n} \right) ds(\boldsymbol{\xi}) \\\
&= \sum_{j=1}^{N} \int_{\Gamma_j} \left( \nabla_{\boldsymbol{x}} G(\boldsymbol{x}, \boldsymbol{\xi}) \frac{\partial \phi(\boldsymbol{\xi})}{\partial n} - \phi(\boldsymbol{\xi}) \nabla_{\boldsymbol{x}} \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n} \right) ds(\boldsymbol{\xi})
\end{aligned} \tag{2.8}
以降の流れ
式 $(2.8)$ を用いて速度ベクトルを計算するためには、以下の4つの量を求める必要があります。本稿では、前者2つ（基本解の微分）は数学的に導出し、後者2つ（境界上の物理量）はBEMによる数値計算で求めます。
| 量 | 説明 | 解説箇所 |
|---|---|---|
| $\nabla_{\boldsymbol{x}} G(\boldsymbol{x}, \boldsymbol{\xi})$ | 基本解の勾配 | 節3.1 |
| $\nabla_{\boldsymbol{x}} \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n}$ | 基本解の法線微分の勾配 | 節3.1 |
| $\phi(\boldsymbol{\xi})$ | 境界上のポテンシャル | 節3.2 |
| $\frac{\partial \phi(\boldsymbol{\xi})}{\partial n}$ | 境界上の法線方向流速 | 節3.2 |
3. BEMによる計算
3.1 基本解とその微分の導出
前述の通り、2次元ラプラス方程式の基本解は式 $(2.5)$ で与えられます。式 $(2.8)$ から、速度ベクトル $\boldsymbol{v}(\boldsymbol{x})$ を計算するためには、基本解の $\boldsymbol{x}$ に関する勾配 $\nabla_{\boldsymbol{x}} G(\boldsymbol{x}, \boldsymbol{\xi})$ と、基本解の法線方向微分の $\boldsymbol{x}$ に関する勾配 $\nabla_{\boldsymbol{x}} \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n}$ を求める必要があります。
まず、$r$ の $\boldsymbol{x}, \boldsymbol{\xi}$ に関する勾配は以下のように展開されます。計算の詳細は補遺Aをご参照ください。
$$\nabla_{\boldsymbol{x}} r = -\frac{\boldsymbol{r}}{r}, \quad \nabla_{\boldsymbol{\xi}} r = \frac{\boldsymbol{r}}{r} \tag{3.1}$$
したがって、基本解の $\boldsymbol{x}$ に関する勾配 $\nabla_{\boldsymbol{x}} G(\boldsymbol{x}, \boldsymbol{\xi})$、および基本解の法線方向微分 $\frac{\partial G}{\partial n}$ は以下のようになります。
\begin{aligned}
\nabla_{\boldsymbol{x}} G(\boldsymbol{x}, \boldsymbol{\xi}) &= -\frac{1}{2\pi} \nabla_{\boldsymbol{x}} \ln r \\\
&= -\frac{1}{2\pi} \cdot \frac{1}{r} \nabla_{\boldsymbol{x}} r \\\
&= -\frac{1}{2\pi} \cdot \frac{1}{r} \left( -\frac{\boldsymbol{r}}{r} \right) \\\
&= \frac{1}{2\pi r^2} \boldsymbol{r}
\end{aligned} \tag{3.2}
\begin{aligned}
\frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n} &= \nabla_{\boldsymbol{\xi}} G(\boldsymbol{x}, \boldsymbol{\xi}) \cdot \boldsymbol{n} \\\
&= \left(-\frac{1}{2\pi} \cdot \frac{1}{r} \nabla_{\boldsymbol{\xi}} r\right) \cdot \boldsymbol{n} \\\
&= \left(-\frac{1}{2\pi} \cdot \frac{1}{r} \cdot \frac{\boldsymbol{r}}{r}\right) \cdot \boldsymbol{n} \\\
&= -\frac{\boldsymbol{r} \cdot \boldsymbol{n}}{2\pi r^2}
\end{aligned} \tag{3.3}
次に、このスカラー量を $\boldsymbol{x}$ で微分します。商の微分公式を用いると、
$$\nabla_{\boldsymbol{x}} \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n} = -\frac{1}{2\pi} \frac{r^2 \nabla_{\boldsymbol{x}} (\boldsymbol{r} \cdot \boldsymbol{n}) - (\boldsymbol{r} \cdot \boldsymbol{n}) \nabla_{\boldsymbol{x}} (r^2)}{r^4} \tag{3.4}$$
となります。ここで、
$$\nabla_{\boldsymbol{x}} (\boldsymbol{r} \cdot \boldsymbol{n}) = -\boldsymbol{n}, \quad \nabla_{\boldsymbol{x}} (r^2) = -2\boldsymbol{r} \quad\text{(3.5)}$$
です（補遺A）。これらを式 $(3.4)$ に代入することで、式 $(3.6)$ が得られます。
\begin{aligned}
\nabla_{\boldsymbol{x}} \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n} &= -\frac{1}{2\pi} \frac{r^2 (-\boldsymbol{n}) - (\boldsymbol{r} \cdot \boldsymbol{n}) (-2\boldsymbol{r})}{r^4} \\\
&= \frac{1}{2\pi r^2} \left( \boldsymbol{n} - \frac{2 (\boldsymbol{r} \cdot \boldsymbol{n})}{r^2} \boldsymbol{r} \right)
\end{aligned} \tag{3.6}
このとき、式 $(2.5)$, 式 $(3.3)$ を利用すると、内部点 $\boldsymbol{x}$ におけるポテンシャル値 $\phi(\boldsymbol{x})$ の式 $(2.6)$ は以下のように展開できます。詳細な計算は補遺Aを参照してください。
$$c(\boldsymbol{x})\phi(\boldsymbol{x}) = \frac{1}{2\pi}\sum_{j=1}^{N}\left(-q_j \left( \frac{s_2}{2} \ln r_2^2 - \frac{s_1}{2} \ln r_1^2 - (s_2 - s_1) \right) + (\phi_j - q_j d) \Delta\theta_j \right) \tag{3.7}$$
ここで、各記号は以下のように定義されます。
- $\boldsymbol{r}_i, r_i \ (i = 1, 2)$：点 $\boldsymbol{x}$ から線分要素 $\Gamma_j$ の各端点へ向かうベクトルとその大きさ
- $\boldsymbol{t}$：線分要素 $\Gamma_j$ における単位接線ベクトル
- $s_i = \boldsymbol{r}_i \cdot \boldsymbol{t} \ (i = 1, 2)$：各端点への相対ベクトルの接線方向成分
- $d = \boldsymbol{r}_i \cdot \boldsymbol{n} \ (i = 1, 2)$：観測点から線分要素への符号付最短距離
- $\Delta \theta_j = \arctan \frac{s_2}{d} - \arctan \frac{s_1}{d}$：点 $\boldsymbol{x}$ から見た線分要素 $\Gamma_j$ の両端点に対する角度差
プログラム上では、角度差 $\Delta \theta_j$ は
arctan2関数を用いて計算され、$(-\pi, \pi]$ の範囲に正規化されます。
3.2 境界積分方程式の離散化
本節では、境界積分方程式 $(2.6)$ を離散化し、連立一次方程式として解く手順を説明します。
境界の離散化
図3.1のように、閉曲線 $\Gamma$ を $N$ 個の線分要素 $\Gamma_j \ (j = 1, 2, \ldots, N)$ に分割します。各要素 $\Gamma_j$ 上では、ポテンシャル $\phi$ と法線方向流速 $q = \partial \phi / \partial n$ を一定と仮定します。この仮定により、境界上の未知量を有限個の値で近似できます。
$$\phi(\boldsymbol{\xi}) \approx \phi_j, \quad \frac{\partial \phi}{\partial n}(\boldsymbol{\xi}) = q(\boldsymbol{\xi}) \approx q_j, \quad \left( \boldsymbol{\xi} \in \Gamma_j \right) \tag{3.8}$$
図 3.1 : 境界の離散化. 外部の線分要素を $N_1$ 個、内部の線分要素を $N_2$ 個に分割している例 ($N = N_1 + N_2$)
離散化された境界積分方程式
各要素の中点を $\boldsymbol{\xi}_j$、長さを $l_j$ と表します。中点 $\boldsymbol{\xi}_i$ における境界積分方程式 $(2.6)$ を離散化すると、以下のようになります。ここで、$\phi$ の項を左辺に、$q$ の項を右辺に移項しています。
\begin{aligned}
&\quad c(\boldsymbol{\xi}_i) \phi(\boldsymbol{\xi}_i) + \oint_{\Gamma} \phi(\boldsymbol{\xi}) \frac{\partial G(\boldsymbol{\xi}_i, \boldsymbol{\xi})}{\partial n} ds(\boldsymbol{\xi}) = \oint_{\Gamma} G(\boldsymbol{\xi}_i, \boldsymbol{\xi}) q(\boldsymbol{\xi}) ds(\boldsymbol{\xi}) \\\
&\Rightarrow c(\boldsymbol{\xi}_i) \phi_i + \sum_{j=1}^{N} \phi_j \int_{\Gamma_j} \frac{\partial G(\boldsymbol{\xi}_i, \boldsymbol{\xi})}{\partial n} ds(\boldsymbol{\xi}) = \sum_{j=1}^{N} q_j \int_{\Gamma_j} G(\boldsymbol{\xi}_i, \boldsymbol{\xi}) ds(\boldsymbol{\xi}) \\\
&\Rightarrow \sum_{j=1}^{N} H_{ij} \phi_j = \sum_{j=1}^{N} G_{ij} q_j
\end{aligned} \tag{3.9}
ここで、行列要素 $H_{ij}$ と $G_{ij}$ は、クロネッカーのデルタ $\delta_{ij}$ を用いて以下のように定義されます。計算の詳細は補遺Aをご参照ください。
\begin{aligned}
H_{ij} &= \int_{\Gamma_j} \frac{\partial G(\boldsymbol{\xi}_i, \boldsymbol{\xi})}{\partial n} ds(\boldsymbol{\xi}) + \delta_{ij} c(\boldsymbol{\xi}_i) \\\
&= \int_{\Gamma_j} \left( -\frac{\boldsymbol{r} \cdot \boldsymbol{n}}{2\pi r^2} \right) ds(\boldsymbol{\xi}) + \delta_{ij} c(\boldsymbol{\xi}_i) \\\
G_{ij} &= \int_{\Gamma_j} G(\boldsymbol{\xi}_i, \boldsymbol{\xi}) ds(\boldsymbol{\xi}) \\\
&= \int_{\Gamma_j} \left( -\frac{1}{2\pi} \ln r \right) ds(\boldsymbol{\xi})
\end{aligned} \tag{3.10}
対角項 $H_{ii}$ について：境界を線分で近似しているため、境界上の点 $\boldsymbol{\xi} _j$ において位置ベクトル $\boldsymbol{r}$ と法線ベクトル $\boldsymbol{n}$ は直交します。したがって $\boldsymbol{r} \cdot \boldsymbol{n} = 0$ となり、$H _{ii} = c(\boldsymbol{\xi} _i) = 1/2$ となります。
行列方程式への変換
式 $(3.9)$ を $i = 1, 2, \ldots, N$ について連立させると、以下の行列方程式が得られます。
\begin{bmatrix}
H_{11} & H_{12} & \cdots & H_{1N} \\\
H_{21} & H_{22} & \cdots & H_{2N} \\\
\vdots & \vdots & \ddots & \vdots \\\
H_{N1} & H_{N2} & \cdots & H_{NN}
\end{bmatrix}
\begin{bmatrix}
\phi_1 \\\ \phi_2 \\\ \vdots \\\ \phi_N
\end{bmatrix} =
\begin{bmatrix}
G_{11} & G_{12} & \cdots & G_{1N} \\\
G_{21} & G_{22} & \cdots & G_{2N} \\\
\vdots & \vdots & \ddots & \vdots \\\
G_{N1} & G_{N2} & \cdots & G_{NN}
\end{bmatrix}
\begin{bmatrix}
q_1 \\\ q_2 \\\ \vdots \\\ q_N
\end{bmatrix} \tag{3.11}
連立一次方程式の構築
節2.2で説明したように、各要素 $\Gamma_j$ において $\phi_j$ または $q_j$ のどちらかが境界条件として指定されます。したがって、式 $(3.11)$ は未知量に関する連立一次方程式として解くことができます。
具体的には、既知の値を右辺ベクトル $\mathbf{b}$ にまとめ、未知量を左辺ベクトル $\mathbf{x}$ としてまとめることで、以下の形に変換できます。
$$A \mathbf{x} = \mathbf{b} \tag{3.12}$$
具体例：四角形領域での設定
境界 $\Gamma$ が4つの要素 $\Gamma_1, \Gamma_2, \Gamma_3, \Gamma_4$ に分割されており、$\Gamma_1$ と $\Gamma_3$ でポテンシャル $\phi$ が指定され（ディリクレ条件）、$\Gamma_2$ と $\Gamma_4$ で流速 $q$ が指定されている（ノイマン条件）場合を考えます。これは、四角形領域の2辺を開口部、残りの2辺を壁面として設定した場合に対応します（図 3.2）。
図 3.2 : 台形領域における各線分要素の境界条件の設定例
このとき、式 $(3.11)$ は以下のように変形できます。ここで、既知の変数には添え字 "k" を付けています。
\begin{bmatrix}
H_{11} & H_{12} & H_{13} & H_{14} \\\
H_{21} & H_{22} & H_{23} & H_{24} \\\
H_{31} & H_{32} & H_{33} & H_{34} \\\
H_{41} & H_{42} & H_{43} & H_{44}
\end{bmatrix}
\begin{bmatrix}
\phi_{1,k} \\\ \phi_2 \\\ \phi_{3,k} \\\ \phi_4
\end{bmatrix} =
\begin{bmatrix}
G_{11} & G_{12} & G_{13} & G_{14} \\\
G_{21} & G_{22} & G_{23} & G_{24} \\\
G_{31} & G_{32} & G_{33} & G_{34} \\\
G_{41} & G_{42} & G_{43} & G_{44}
\end{bmatrix}
\begin{bmatrix}
q_1 \\\ q_{2,k} \\\ q_3 \\\ q_{4,k}
\end{bmatrix}
既知の値を右辺に移項し、未知量を左辺にまとめると、式 $(3.12)$ の形になります。
\begin{bmatrix}
-G_{11} & H_{12} & -G_{13} & H_{14} \\\
-G_{21} & H_{22} & -G_{23} & H_{24} \\\
-G_{31} & H_{32} & -G_{33} & H_{34} \\\
-G_{41} & H_{42} & -G_{43} & H_{44}
\end{bmatrix} \begin{bmatrix}
q_1 \\\ \phi_2 \\\ q_3 \\\ \phi_4
\end{bmatrix} = \begin{bmatrix}
- H_{11} \phi_{1,k} + G_{12} q_{2,k} - H_{13} \phi_{3,k} + G_{14} q_{4,k} \\\
- H_{21} \phi_{1,k} + G_{22} q_{2,k} - H_{23} \phi_{3,k} + G_{24} q_{4,k} \\\
- H_{31} \phi_{1,k} + G_{32} q_{2,k} - H_{33} \phi_{3,k} + G_{34} q_{4,k} \\\
- H_{41} \phi_{1,k} + G_{42} q_{2,k} - H_{43} \phi_{3,k} + G_{44} q_{4,k}
\end{bmatrix} \tag{3.13}
壁面では $q_{2,k} = q_{4,k} = 0$ とし、開口部では例えば $\phi_{1,k} = 0$（入口）、$\phi_{3,k} > 0$（出口）と設定します。この連立方程式を解くことで、$\Gamma_1$ から $\Gamma_3$ へ向かう流れが計算されます。
要素分割について：上の例のように1辺を1つの線分要素として設定すると、辺が長い場合に計算精度が低下することがあります。より高い精度を得るためには、各辺を複数の線分要素に分割することが望ましいです。本プログラムでは、境界条件設定時に最大の線分長さ (BoundaryConditionsコンストラクタのmax_length引数) に基づき、自動的に要素分割を行っています。そのため、形状定義で与えた線分の数よりも、実際に使用される線分要素の数が多くなる場合があります。
3.3 内部点での速度ベクトルの計算
前節で境界上の未知量 $\phi_j$ と $q_j$ が求まりました。本節では、これらの値を用いて、任意の内部点 $\boldsymbol{x}$ における速度ベクトル $\boldsymbol{v}(\boldsymbol{x})$ を計算する方法を説明します。
速度ベクトルの導出
式 $(2.8)$ に、基本解の勾配 $(3.2)$ と基本解の法線微分の勾配 $(3.6)$ を代入すると、速度ベクトルは以下のように表されます。
\begin{aligned}
\boldsymbol{v}(\boldsymbol{x}) &= \sum_{j=1}^{N} \int_{\Gamma_j} \left( \nabla_{\boldsymbol{x}} G(\boldsymbol{x}, \boldsymbol{\xi}) \frac{\partial \phi(\boldsymbol{\xi})}{\partial n} - \phi(\boldsymbol{\xi}) \nabla_{\boldsymbol{x}} \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n} \right) ds(\boldsymbol{\xi}) \\\
&= \sum_{j=1}^{N} \int_{\Gamma_j} \left( \frac{1}{2\pi r^2} \boldsymbol{r} \cdot q_j - \phi_j \cdot \frac{1}{2\pi r^2} \left( \boldsymbol{n} - \frac{2 (\boldsymbol{r} \cdot \boldsymbol{n})}{r^2} \boldsymbol{r} \right) \right) ds(\boldsymbol{\xi}) \\\
&= \frac{1}{2\pi} \sum_{j=1}^{N} \left( q_j \int_{\Gamma_j} \frac{\boldsymbol{r}}{r^2} ds(\boldsymbol{\xi}) - \phi_j \int_{\Gamma_j} \frac{1}{r^2} \left( \boldsymbol{n} - \frac{2 (\boldsymbol{r} \cdot \boldsymbol{n})}{r^2} \boldsymbol{r} \right) ds(\boldsymbol{\xi}) \right)
\end{aligned} \tag{3.13}
積分の解析的評価
式 $(3.13)$ に現れる2つの積分は、解析的に計算することができます。したがって、速度ベクトルは以下のように表されます。導出の詳細については、補遺Aをご参照ください。
\boldsymbol{v}(\boldsymbol{x}) = \frac{1}{2\pi} \sum_{j=1}^{N} \left(
\quad \left( q_j \ln \frac{r_2}{r_1} - \phi_j \left( \frac{d}{r_2^2} - \frac{d}{r_1^2} \right) \right) \boldsymbol{t} \\\
+ \left( q_j (\theta_2 - \theta_1) - \phi_j \left( -\frac{s_2}{r_2^2} + \frac{s_1}{r_1^2} \right) \right) \boldsymbol{n} \quad \right) \tag{3.15}
ここで、$s_1 = \boldsymbol{r}_1 \cdot \boldsymbol{t}$、$s_2 = \boldsymbol{r}_2 \cdot \boldsymbol{t}$, および $d = \boldsymbol{r}_1 \cdot \boldsymbol{n} = \boldsymbol{r}_2 \cdot \boldsymbol{n}$ である。
角度差 $(\theta_2 - \theta_1)$ は、arctan2 関数を用いて計算します。要素の両端点の座標を $(x_{j1}, y_{j1})$, $(x_{j2}, y_{j2})$、観測点の座標を $(x_i, y_i)$ とすると、以下のように求められます。
$$\begin{aligned}\theta_1 &= \text{atan2}(y_{j1} - y_i, x_{j1} - x_i) \\ \theta_2 &= \text{atan2}(y_{j2} - y_i, x_{j2} - x_i)\end{aligned}$$
この差分 $(\theta_2 - \theta_1)$ を $[-\pi, \pi]$ の範囲に正規化し、式 $(3.15)$ に代入することで速度ベクトルを求めます。
参考文献
[1] Tara LaForce, "PE281 Boundary Element Method Course Notes", Stanford University, 2006, https://web.stanford.edu/class/energy281/BoundaryElementMethod.pdf
Appendix
A. 計算
本節では、本文で述べた一部の式について、これを整理・導出する過程を示します。
A-Eq3.1
観測点 $\boldsymbol{x}$ と境界上の点 $\boldsymbol{\xi}$ に関する相対ベクトルは $\boldsymbol{r} = \boldsymbol{\xi} - \boldsymbol{x}$ です。本節では、$r = |\boldsymbol{r}|$ の $\boldsymbol{x}$ および $\boldsymbol{\xi}$ に関する勾配を求めます。
2次元の場合
$\boldsymbol{x} = (x, y)$、$\boldsymbol{\xi} = (\xi_x, \xi_y)$ とすれば $\boldsymbol{r} = (\xi_x - x, \xi_y - y)$、$r = \sqrt{(\xi_x - x)^2 + (\xi_y - y)^2}$ となるので、$r$ の $\boldsymbol{x}, \boldsymbol{\xi}$ に関する勾配は以下のように求められます。
\begin{aligned}
\nabla_{\boldsymbol{x}} r &= \left( \frac{\partial r}{\partial x}, \frac{\partial r}{\partial y} \right) \\\
&= \left( \frac{-(\xi_x - x)}{r}, \frac{-(\xi_y - y)}{r} \right) \\\
&= -\frac{\boldsymbol{r}}{r}
\end{aligned} \tag{A31-1}
\begin{aligned}
\nabla_{\boldsymbol{\xi}} r &= \left( \frac{\partial r}{\partial \xi_x}, \frac{\partial r}{\partial \xi_y} \right) \\\
&= \left( \frac{\xi_x - x}{r}, \frac{\xi_y - y}{r} \right) \\\
&= \frac{\boldsymbol{r}}{r}
\end{aligned} \tag{A31-2}
3次元の場合
$\boldsymbol{x} = (x, y, z)$、$\boldsymbol{\xi} = (\xi_x, \xi_y, \xi_z)$ とすれば $\boldsymbol{r} = (\xi_x - x, \xi_y - y, \xi_z - z)$、$r = \sqrt{(\xi_x - x)^2 + (\xi_y - y)^2 + (\xi_z - z)^2}$ となるので、$r$ の $\boldsymbol{x}, \boldsymbol{\xi}$ に関する勾配は以下のように求められます。
\begin{aligned}
\nabla_{\boldsymbol{x}} r &= \left( \frac{\partial r}{\partial x}, \frac{\partial r}{\partial y}, \frac{\partial r}{\partial z} \right) \\\
&= \left( \frac{-(\xi_x - x)}{r}, \frac{-(\xi_y - y)}{r}, \frac{-(\xi_z - z)}{r} \right) \\\
&= -\frac{\boldsymbol{r}}{r}
\end{aligned} \tag{A31-3}
\begin{aligned}
\nabla_{\boldsymbol{\xi}} r &= \left( \frac{\partial r}{\partial \xi_x}, \frac{\partial r}{\partial \xi_y}, \frac{\partial r}{\partial \xi_z} \right) \\\
&= \left( \frac{\xi_x - x}{r}, \frac{\xi_y - y}{r}, \frac{\xi_z - z}{r} \right) \\\
&= \frac{\boldsymbol{r}}{r}
\end{aligned} \tag{A31-4}
A-Eq3.5
式 $(3.6)$ で使用する、$\nabla_{\boldsymbol{x}} (\boldsymbol{r} \cdot \boldsymbol{n})$ と $\nabla_{\boldsymbol{x}} (r^2)$ を求めます。
2次元の場合
\begin{aligned}
\nabla_{\boldsymbol{x}} (\boldsymbol{r} \cdot \boldsymbol{n}) &= \nabla_{\boldsymbol{x}} ((\xi_x - x) n_x + (\xi_y - y) n_y) = -\boldsymbol{n} \\
\nabla_{\boldsymbol{x}} (r^2) &= \nabla_{\boldsymbol{x}} ((\xi_x - x)^2 + (\xi_y - y)^2) = -2\boldsymbol{r}
\end{aligned} \tag{A35-1}
3次元の場合
\begin{aligned}
\nabla_{\boldsymbol{x}} (\boldsymbol{r} \cdot \boldsymbol{n}) &= \nabla_{\boldsymbol{x}} ((\xi_x - x) n_x + (\xi_y - y) n_y + (\xi_z - z) n_z) = -\boldsymbol{n} \\
\nabla_{\boldsymbol{x}} (r^2) &= \nabla_{\boldsymbol{x}} ((\xi_x - x)^2 + (\xi_y - y)^2 + (\xi_z - z)^2) = -2\boldsymbol{r}
\end{aligned} \tag{A35-2}
A-Eq3.7
式 $(2.6)$ を展開します。式 $(2.5)$, 式 $(3.3)$ を代入すると、以下のようになります。
\begin{aligned}
c(\boldsymbol{x}) \phi(\boldsymbol{x}) &= \oint_{\Gamma} \left(G(\boldsymbol{x}, \boldsymbol{\xi}) \frac{\partial \phi(\boldsymbol{\xi})}{\partial n} - \phi(\boldsymbol{\xi}) \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n} \right) ds(\boldsymbol{\xi}) \\\
&= \sum_{j=1}^{N} \int_{\Gamma_j} \left( G(\boldsymbol{x}, \boldsymbol{\xi}) q_j - \phi_j \frac{\partial G(\boldsymbol{x}, \boldsymbol{\xi})}{\partial n} \right) ds(\boldsymbol{\xi}) \\\
&= \sum_{j=1}^{N} \int_{\Gamma_j} \left( -\frac{q_j}{2\pi} \ln r + \phi_j \frac{\boldsymbol{r} \cdot \boldsymbol{n}}{2\pi r^2} \right) ds(\boldsymbol{\xi})
\end{aligned} \tag{A26-1}
この計算のため、相対ベクトル $\boldsymbol{r}$ を分解し、局所座標系 $s$ を導入します。線分要素 $\Gamma_j$ の単位接線ベクトルを $\boldsymbol{t}$、単位法線ベクトルを $\boldsymbol{n}$ とすれば、$\boldsymbol{r}$ は以下のように分解できます。
\begin{aligned}
&\boldsymbol{r} = s \boldsymbol{t} + d \boldsymbol{n} \\\
& \text{where} \quad \left\lbrace \begin{aligned}
s &= \boldsymbol{r} \cdot \boldsymbol{t} \\\
d &= \boldsymbol{r} \cdot \boldsymbol{n}
\end{aligned} \right.
\end{aligned} \tag{A26-2}
ここで、$d$ は観測点 $\boldsymbol{x}$ から線分要素 $\Gamma_j$ までの符号付最短距離であり、$s$ は線分要素上の局所座標です。したがって、$r = \sqrt{s^2 + d^2}$ となります。以降、積分変数を $ds(\boldsymbol{\xi}) = ds$ とし、線分の始点を $s = s_1$、終点を $s = s_2$ と表します。
これらを式 $\text{(A26-1)}$ に代入すると、要素 $\Gamma_j$ に関する積分は以下のように書き換えられます。
\begin{aligned}
I_j &= \int_{\Gamma_j} \left( -\frac{q_j}{2\pi} \ln r + \phi_j \frac{\boldsymbol{r} \cdot \boldsymbol{n}}{2\pi r^2} \right) ds(\boldsymbol{\xi}) \\\
&= -\frac{q_j}{2\pi} \int_{s_1}^{s_2} \ln \sqrt{s^2 + d^2} ds + \frac{\phi_j}{2\pi} \int_{s_1}^{s_2} \frac{d}{s^2 + d^2} ds \\\
\end{aligned} \tag{A26-3}
これらの積分は、以下のように計算できます。
\begin{aligned}
\int_{s_1}^{s_2} \ln \sqrt{s^2 + d^2} ds &= \frac{1}{2} \int_{s_1}^{s_2} \ln (s^2 + d^2) ds \\\
&= \frac{1}{2} \left[ s \ln (s^2 + d^2) - 2s + 2d \arctan \left( \frac{s}{d} \right) \right]_{s_1}^{s_2} \\\
&= \frac{s_2}{2} \ln r_2^2 - \frac{s_1}{2} \ln r_1^2 - (s_2 - s_1) + d \left( \arctan \frac{s_2}{d} - \arctan \frac{s_1}{d} \right) \\\
\int_{s_1}^{s_2} \frac{d}{s^2 + d^2} ds &= \left[ \arctan \left( \frac{s}{d} \right) \right]_{s_1}^{s_2} \\\
&= \arctan \frac{s_2}{d} - \arctan \frac{s_1}{d}
\end{aligned}
以降、$\Delta\theta_j = \arctan \frac{s_2}{d} - \arctan \frac{s_1}{d}$ と表します。これらを式 $\text{(A26-3)}$ に代入すると、以下が得られ、これをを全要素について合計することで、式 $(2.6)$ の右辺が計算されます。
\begin{aligned}
&I_j = -\frac{q_j}{2\pi} \left( \frac{s_2}{2} \ln r_2^2 - \frac{s_1}{2} \ln r_1^2 - (s_2 - s_1) \right) + \frac{\phi_j - q_j d}{2\pi} \Delta\theta_j \\\
\end{aligned} \tag{A26-4}
実際のプログラムでは、
arctan2関数を用いて $\Delta \theta_j$ を計算しています。
A-Eq3.10
式 $(3.10)$ の積分を計算します。
\begin{aligned}
H_{ij} &= \int_{\Gamma_j} \left( -\frac{\boldsymbol{r} \cdot \boldsymbol{n}}{2\pi r^2} \right) ds(\boldsymbol{\xi}) + \delta_{ij} c(\boldsymbol{\xi}_i) \\\
G_{ij} &= \int_{\Gamma_j} \left( -\frac{1}{2\pi} \ln r \right) ds(\boldsymbol{\xi})
\end{aligned} \tag{3.10}
$H_{ij}$ の計算：
まず、対角項 ($i = j$) について考えます。このとき、前述の通り $\boldsymbol{r} \cdot \boldsymbol{n} = 0$ となるため、
$$H_{ii} = c(\boldsymbol{\xi}_i) = \frac{1}{2} \tag{A310-1}$$
とできます。次に、非対角項 ($i \neq j$) について考えます。このとき、$\boldsymbol{r}$ と $\boldsymbol{n}$ のなす角を $\theta$ とすると、以下のように変形できます。
\begin{aligned}
H_{ij} &= \int_{\Gamma_j} \left( -\frac{\boldsymbol{r} \cdot \boldsymbol{n}}{2\pi r^2} \right) ds(\boldsymbol{\xi}) \\\
&= -\frac{1}{2\pi} \int_{\Gamma_j} \frac{r \cos \theta}{r^2} ds(\boldsymbol{\xi}) \\\
&= -\frac{1}{2\pi} \int_{\Gamma_j} \frac{\cos \theta}{r} ds(\boldsymbol{\xi})
\end{aligned}
ここで、観測点 $i$ から見て、微小要素 $ds$ が張る角度を $d\alpha$ とすると、演習の接線方向に対する要素 $ds$ の射影成分は $ds \cos \theta$ であり、これが弧長 $r d\alpha$ に等しいことから、$ds \cos \theta = r d\alpha$ が成り立ちます。したがって、上の式は以下のように変形できます。
\begin{aligned}
H_{ij} &= -\frac{1}{2\pi} \int_{\Gamma_j} \frac{\cos \theta}{r} ds(\boldsymbol{\xi}) \\\
&= -\frac{1}{2\pi} \int_{\alpha_1}^{\alpha_2} d\alpha \\\
&= -\frac{1}{2\pi} (\alpha_2 - \alpha_1)
\end{aligned} \tag{A310-2}
すなわち、線分要素 $\Gamma_j$ の両端点 $\boldsymbol{\xi} _{j1}, \boldsymbol{\xi} _{j2}$ が観測点 $\boldsymbol{\xi} _i$ となす角度 $(\alpha_2 - \alpha_1)$ を用いて、$H _{ij}$ を表現できます。
プログラムでは、要素の両端点の座標を $(x_{j1}, y_{j1}), (x_{j2}, y_{j2})$、観測点の座標を $(x_i, y_i)$ として、以下のように計算できます。この差分 $(\alpha_2 - \alpha_1)$ を $[-\pi, \pi]$ の範囲に正規化し、上の式に代入することで $H_{ij}$ を求めます。
$$\begin{aligned}\alpha_1 &= \text{atan2}(y_{j1} - y_i, x_{j1} - x_i) \\ \alpha_2 &= \text{atan2}(y_{j2} - y_i, x_{j2} - x_i)\end{aligned}$$
$G_{ij}$ の計算：
まず、対角項 ($i = j$) について考えます。線分要素 $\Gamma_j$ の長さを $l_j$ とすると、中点 $\boldsymbol{\xi}_j$ を原点に取った局所座標 $s \in \left[-\frac{l_j}{2}, \frac{l_j}{2}\right]$ を導入できます。このとき $r = |s|$ となるから、以下のように計算できます。
\begin{aligned}
G_{ii} &= -\frac{1}{2\pi} \int_{-\frac{l_j}{2}}^{\frac{l_j}{2}} \ln |s| ds \\\
&= -\frac{1}{\pi} \int_{0}^{\frac{l_j}{2}} \ln s \ ds \\\
&= -\frac{1}{\pi} \left[ s \ln s - s \right]_{0}^{\frac{l_j}{2}} \\\
&= -\frac{l_j}{2\pi} \left( 1 - \ln\frac{l_j}{2} \right)
\end{aligned} \tag{A310-3}
次に、非対角項 ($i \neq j$) について考えます。観測点 $\boldsymbol{\xi}_i$ から線分 $\Gamma_j$ までの最短距離を $d$ とすると、$r = \sqrt{t^2 + d^2}$ となります。両端点における局所座標 $t$ の値をそれぞれ $t_1, t_2$ とすると、以下のように計算できます。
\begin{aligned}
G_{ij} &= -\frac{1}{2\pi} \int_{t_1}^{t_2} \ln \sqrt{t^2 + d^2} \ dt \\\
&= -\frac{1}{4\pi} \int_{t_1}^{t_2} \ln (t^2 + d^2) \ dt
\end{aligned}
ここで、積分部分を $I$ と置いて部分積分を行うと、
\begin{aligned}
I &= \int \ln (t^2 + d^2) \ dt \\\
&= [t \ln (t^2 + d^2)] - \int t \cdot \frac{2t}{t^2 + d^2} \ dt \\\
&= t\ln(t^2 + d^2) - 2 \int \frac{t^2}{t^2 + d^2} \ dt \\\
&= t\ln(t^2 + d^2) - 2 \left( t - d^2 \int \frac{1}{t^2 + d^2} \ dt \right) \\\
&= t\ln(t^2 + d^2) - 2t + 2d\ \text{arctan} \frac{t}{d} \\\
&= 2\left(t \ln r - t + d\ \text{arctan} \frac{t}{d} \right)
\end{aligned}
したがって、$G_{ij}$ は以下のように求められる。
\begin{aligned}
G_{ij} &= -\frac{1}{4\pi} [I]_{t_1}^{t_2} \\\
&= -\frac{1}{2\pi} \left[ t \ln r - t + d\ \text{arctan} \frac{t}{d} \right]_{t_1}^{t_2} \\\
&= \frac{1}{2\pi} \left( t_1 \ln r_1 - t_2 \ln r_2 + (t_2 - t_1) - d(\theta_2 - \theta_1) \right)
\end{aligned} \tag{A310-4}
A-Eq3.15
式 $(3.13)$ の積分を計算します。
$$\boldsymbol{v}(\boldsymbol{x}) = \frac{1}{2\pi} \sum_{j=1}^{N} \left( q_j \int_{\Gamma_j} \frac{\boldsymbol{r}}{r^2} ds(\boldsymbol{\xi}) - \phi_j \int_{\Gamma_j} \frac{1}{r^2} \left( \boldsymbol{n} - \frac{2 (\boldsymbol{r} \cdot \boldsymbol{n})}{r^2} \boldsymbol{r} \right) ds(\boldsymbol{\xi}) \right) \tag{3.13}$$
この計算のため、A-Eq3.7と同様に、相対ベクトル $\boldsymbol{r}$ を分解し、局所座標系 $s$ を導入します。以降、積分変数を $ds(\boldsymbol{\xi}) = ds$ とし、線分の始点を $s = s_1$、終点を $s = s_2$ と表します。
(1) $\int_{\Gamma_j} \frac{\boldsymbol{r}}{r^2} ds(\boldsymbol{\xi})$ の計算：
\begin{aligned}
\int_{\Gamma_j} \frac{\boldsymbol{r}}{r^2} ds &= \int_{s_1}^{s_2} \frac{s \boldsymbol{t} + d \boldsymbol{n}}{s^2 + d^2} ds \\\
&= \boldsymbol{t} \int_{s_1}^{s_2} \frac{s}{s^2 + d^2} ds + \boldsymbol{n} \int_{s_1}^{s_2} \frac{d}{s^2 + d^2} ds \\\
&= \boldsymbol{t} \left[ \ln(r) \right]_{s_1}^{s_2} + \boldsymbol{n} \left[ \text{arctan} \frac{s}{d} \right]_{s_1}^{s_2} \\\
&= \left( \ln \frac{r_2}{r_1} \right) \boldsymbol{t} + \left( \theta_2 - \theta_1 \right) \boldsymbol{n}
\end{aligned} \tag{A313-2}
(2) $\int_{\Gamma_j} \frac{1}{r^2} \left( \boldsymbol{n} - \frac{2 (\boldsymbol{r} \cdot \boldsymbol{n})}{r^2} \boldsymbol{r} \right) ds(\boldsymbol{\xi})$ の計算：
\begin{aligned}
&\quad \int_{\Gamma_j} \frac{1}{r^2} \left( \boldsymbol{n} - \frac{2 (\boldsymbol{r} \cdot \boldsymbol{n})}{r^2} \boldsymbol{r} \right) ds \\\
&= \int_{s_1}^{s_2} \left( \frac{\boldsymbol{n}}{s^2 + d^2} - \frac{2d(s \boldsymbol{t} + d^2 \boldsymbol{n})}{(s^2 + d^2)^2} \right) ds \\\
&= \boldsymbol{t} \int_{s_1}^{s_2} \left( -\frac{2ds}{(s^2 + d^2)^2} \right) ds + \boldsymbol{n} \int_{s_1}^{s_2} \left( \frac{1}{s^2 + d^2} - \frac{2d^2}{(s^2 + d^2)^2} \right) ds \\\
&= \left[ \frac{d}{r^2} \right]_{s_1}^{s_2} \boldsymbol{t} + \left[ -\frac{s}{r^2} \right]_{s_1}^{s_2} \boldsymbol{n} \\\
&= \left( \frac{d}{r_2^2} - \frac{d}{r_1^2} \right) \boldsymbol{t} + \left( -\frac{s_2}{r_2^2} + \frac{s_1}{r_1^2} \right) \boldsymbol{n}
\end{aligned} \tag{A313-3}
以上より、式 $(3.13)$ は以下のように表されます。
\begin{aligned}
\boldsymbol{v}(\boldsymbol{x}) &= \frac{1}{2\pi} \sum_{j=1}^{N} \left( q_j \left( \left( \ln \frac{r_2}{r_1} \right) \boldsymbol{t} + \left( \theta_2 - \theta_1 \right) \boldsymbol{n} \right) - \phi_j \left( \left( \frac{d}{r_2^2} - \frac{d}{r_1^2} \right) \boldsymbol{t} + \left( -\frac{s_2}{r_2^2} + \frac{s_1}{r_1^2} \right) \boldsymbol{n} \right) \right) \\\
&= \frac{1}{2\pi} \sum_{j=1}^{N} \left( \left( q_j \ln \frac{r_2}{r_1} - \phi_j \left( \frac{d}{r_2^2} - \frac{d}{r_1^2} \right) \right) \boldsymbol{t} + \left( q_j (\theta_2 - \theta_1) - \phi_j \left( -\frac{s_2}{r_2^2} + \frac{s_1}{r_1^2} \right) \right) \boldsymbol{n} \right)
\end{aligned} \tag{A313-4}
ここで、$s_1 = \boldsymbol{r}_1 \cdot \boldsymbol{t}$、$s_2 = \boldsymbol{r}_2 \cdot \boldsymbol{t}$, および $d = \boldsymbol{r}_1 \cdot \boldsymbol{n} = \boldsymbol{r}_2 \cdot \boldsymbol{n}$ です。
B. 理論的な背景
B.1 ポテンシャル流れにおけるナビエ–ストークス方程式の簡略化
流体の速度を $\boldsymbol{v}(\boldsymbol{x}, t)$, 流体にかかる圧力を $p(\boldsymbol{x}, t)$ とします。また、流体の密度を $\rho$, 粘性係数を $\mu$, 単位体積あたりの外力を $\boldsymbol{f}(\boldsymbol{x}, t)$ とします。このとき、ナビエ–ストークス方程式は以下のように表されます。
$$\rho \left( \frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \cdot \nabla) \boldsymbol{v} \right) = -\nabla p + \mu \nabla^2 \boldsymbol{v} + \rho \boldsymbol{f} \tag{A11-1}$$
また、連続の式は次のように表されます。
$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \boldsymbol{v}) = 0 \tag{A11-2}$$
非圧縮性を仮定すると、$\rho = \text{const.}$ となるため、式 $\text{(A11-2)}$ は以下のように簡略化されます。
$$\nabla \cdot \boldsymbol{v} = 0 \tag{A11-3}$$
さらに、非回転性 ($\nabla \times \boldsymbol{v} = 0$) を仮定すると、式 $\text{(A11-1)}$ は式 $\text{(A11-4)}$ のように簡略化されます。この時点で、粘性項 $\mu \nabla^2 \boldsymbol{v}$ はゼロになり、粘性の影響は無視されます。
\begin{aligned}
\text{LHS} &= \rho \left( \frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \cdot \nabla) \boldsymbol{v} \right) \\\
&= \rho \left( \frac{\partial \boldsymbol{v}}{\partial t} + \nabla \left( \frac{1}{2} |\boldsymbol{v}|^2 \right) - \boldsymbol{v} \times (\nabla \times \boldsymbol{v}) \right) \\\
&= \rho \left( \frac{\partial \boldsymbol{v}}{\partial t} + \nabla \left( \frac{1}{2} |\boldsymbol{v}|^2 \right) \right) \\\
\text{RHS} &= -\nabla p + \mu \nabla^2 \boldsymbol{v} + \rho \boldsymbol{f} \\\
&= -\nabla p + (\nabla (\nabla \cdot \boldsymbol{v}) - \nabla \times (\nabla \times \boldsymbol{v})) + \rho \boldsymbol{f} \\\
&= -\nabla p + \rho \boldsymbol{f} \\\
\end{aligned}
$$\therefore \quad \rho \left( \frac{\partial \boldsymbol{v}}{\partial t} + \nabla \left( \frac{1}{2} |\boldsymbol{v}|^2 \right) \right) = -\nabla p + \rho \boldsymbol{f} \tag{A11-4}$$
粘性項について：粘性を考慮する場合、壁面では滑りなし ($\boldsymbol{v} = 0$) の境界条件を課す必要があります。この場合、壁面近傍では流体に速度差が生じるため、非回転性の仮定が破れます。したがって、ポテンシャル流れでは、非粘性流体を仮定し、壁面では流体が壁面に対して滑る（すなわち、壁面に垂直な速度成分のみゼロ）と考えることで、非回転性を保ちます。
本シミュレーションは、迷路を解くような（終点領域が存在する場合に、そこへ向かう流れが発生する）流れを計算することを目的としています。この場合、壁面で滑りなしの境界条件を課す必要はなく、ポテンシャル流れの仮定が妥当と考えます。
さらに、流れの定常性 ($\frac{\partial \boldsymbol{v}}{\partial t} = 0$) と外力の不存在 ($\boldsymbol{f} = 0$) を仮定すると、式 $\text{(A11-4)}$ は以下のように簡略化されます（ベルヌーイの定理の微分形）。
$$\frac{\rho}{2} \nabla |\boldsymbol{v}|^2 = -\nabla p \tag{A11-5}$$