3.2 FreeFem++ プログラム (その2) 流体の速度ポテンシャルを求める

速度ポテンシャルを求める境界値問題

(3.4)   $\displaystyle \Laplacian\phi=0$   (in $ \Omega $)
(3.5)   $\displaystyle \frac{\rd\phi}{\rd\bm{n}}=\bm{v}\cdot\bm{n}$   (on $ \rd\Omega$)

の場合は、 $ \Gamma_1=\emptyset$, $ \Gamma_2=\rd\Omega$, $ f=0$, $ g_2=\bm{v}\cdot\bm{n}$.

特に $ \Omega=\left\{(x,y)\in\mathbb{R}^2\mid x^2+y^2<1\right\}$, $ \bm{v}\equiv \begin{pmatrix}
1  2
\end{pmatrix}$ のときは、 $ (x,y)\in\rd\Omega$ において $ \bm{n}=\begin{pmatrix}x  y \end{pmatrix}$ であるから、

$\displaystyle g_2=\bm{v}\cdot\bm{n}=
\begin{pmatrix}
1  2
\end{pmatrix} \cdot
\begin{pmatrix}
x  y
\end{pmatrix} =x+2y.
$

速度ポテンシャルを求める potential2d-v0.edp ← これを保存する

// potential2d-v0.edp
//   http://nalab.mind.meiji.ac.jp/~mk/complex2/potential2d-v0.edp
//   2次元非圧縮ポテンシャル流
//     速度ポテンシャル,速度を求め、等ポテンシャル線, 速度場を描く

border Gamma(t=0,2*pi) { x = cos(t); y = sin(t); } // 円盤領域
int m=40;
mesh Th=buildmesh(Gamma(m));
plot(Th, wait=1, ps="Th.eps");
// 次の2行は区分1次多項式を使うという意味
fespace Vh(Th,P1);
Vh phi, v, v1, v2;
// 境界条件の設定
func Vn=x+2*y; // Ωが単位円で, V=(1,2) のとき V・n=x+2y

// 速度ポテンシャルφを求め、その等高線 (等ポテンシャル線) を描く
solve Laplace(phi,v) =
  int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma)(Vn*v);
plot(phi,ps="contourpotential.eps",wait=1);

// ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
v1=dx(phi); v2=dy(phi);
plot([v1,v2],ps="vectorfield.eps",wait=1);

// 等ポテンシャル線とベクトル場を同時に描く
plot([v1,v2],phi,ps="both.eps", wait=1);
ターミナルでこうして入手できる
curl -O http://nalab.mind.meiji.ac.jp/~mk/complex2/potential2d-v0.edp

プログラムはテキスト・エディター (Cot Editor, テキスト・エディット5, Visual Studio Code, Vim, emacs など) で作成し、ターミナルから、
こんなふうにして実行
FreeFem++ potential2d-v0.edp
とタイプして実行できる。リターン・キーを打つごとに次の図に進み、 最後 (ベクトル場と等ポテンシャル線を描いてお終い) は Escape キーを打って終了する。

図 2: $ \Omega $ の三角形分割
\includegraphics[width=5cm]{Th.eps}

図 3: 一様流の等ポテンシャル線, 流線, ベクトル場
\includegraphics[width=5cm]{potential-graph/contourpotential.eps} \includegraphics[width=5cm]{potential-graph/streamline.eps} \includegraphics[width=5cm]{potential-graph/vectorfield.eps}

このプログラムは、 図の PostScript データ “counterpotential.eps”, “vectorfield.eps”, “both.eps” を出力している。


(実は、上の弱形式は解の一意性がないので ($ \rd\Omega$ 全体で Neumann 境界条件だから)、 このプログラムは少し危ういところがある。)



桂田 祐史