注意

(1)

$v_n:=\bm{v} \cdot\bm{n}$ は $\dsp\int_{\Gamma}v_n\;\D\sigma=0$ を満たしている必要がある。実際、Gaussの発散定理と非圧縮性の仮定から

$$\int_{\Gamma}v_n\;\D\sigma =\int_{\Gamma}\bm{v}\cdot\bm{n}\;\D\sigma =\int_{\Omega}\Div\bm{v}\;\D\bm{x} =\int_{\Omega}0\;\D\bm{x}=0.$$

サンプルプログラムでは、 円盤領域 $\Omega=\left\{(x,y)\in\mathbb{R}^2\relmiddle\vert x^2+y^2<1\right\}$, 一様流 $\bm{v}=\begin{pmatrix}1\\ 2\end{pmatrix}$ であったので、

$$\bm{n}=\begin{pmatrix}x\\ y\end{pmatrix},\quad v_n=\bm{v}\cdot\b... ... \begin{pmatrix}1\\ 2\end{pmatrix}\cdot \begin{pmatrix}x\\ y\end{pmatrix}=x+2y$$

としてある。 $\Div\bm{v}=0$ であるから、 当然 $\dsp\int_{\Gamma}v_n\;\D\sigma=0$ も成り立つ。

(2)

湧き出しや吸い込み、点渦など、特異点が $\Omega$ 内にあるような問題は、 この方法では解くことが出来ない。

potential2d-v0.edp

   1 // potential2d-v0.edp
   2 //   http://nalab.mind.meiji.ac.jp/~mk/complex2/potential2d-v0.edp
   3 //   2次元非圧縮ポテンシャル流
   4 //     速度ポテンシャル，速度を求め、等ポテンシャル線, 速度場を描く
   5 
   6 border Gamma(t=0,2*pi) { x = cos(t); y = sin(t); } // 円盤領域
   7 int m=40;
   8 mesh Th=buildmesh(Gamma(m));
   9 plot(Th, wait=1, ps="Th.eps");
  10 // 次の2行は区分1次多項式を使うという意味
  11 fespace Vh(Th,P1);
  12 Vh phi, v, v1, v2;
  13 // 境界条件の設定
  14 func Vn=x+2*y; // Ωが単位円で, V=(1,2) のとき V・n=x+2y
  15 
  16 // 速度ポテンシャルφを求め、その等高線 (等ポテンシャル線) を描く
  17 solve Laplace(phi,v) =
  18   int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma)(Vn*v);
  19 plot(phi,ps="contourpotential.eps",wait=1);
  20 
  21 // ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
  22 v1=dx(phi); v2=dy(phi);
  23 plot([v1,v2],ps="vectorfield.eps",wait=1);
  24 
  25 // 等ポテンシャル線とベクトル場を同時に描く
  26 plot([v1,v2],phi,ps="both.eps", wait=1);

• 6行目で領域$\Omega$の境界$\Gamma$を指定している。
• 8行目で、$\Gamma$を$m$分割して、$\Gamma$の囲む範囲を三角形分割して、 それを mesh 型の変数 Th に代入している。
• 11行目、有限要素空間 $V_h$ を区分的1次関数の空間と定義している。 この辺についてはこの講義では説明を省略する (知りたい人は有限要素法のテキストを読んで下さい)。
• 12行目、$\phi$ と試験関数 $v$, 流速ベクトル場 $\bm{v}$ の成分 $v_1,v_2$ を、 Vh の要素とする。
• 14行目、境界条件の設定をしている。 ここでは $\bm{v}\cdot\bm{n}=x+2y$ とする理由は上に書いた。
• 17〜18行目、弱形式の定義。ここを修正する必要が生じる可能性は低い。