注意

(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 
  11 fespace Vh(Th,P1);
  12 Vh phi, v, v1, v2;
  13 func Vn=x+2*y; // Ωが単位円で, V=(1,2) のとき V・n=x+2y
  14 
  15 // 速度ポテンシャルφを求め、その等高線 (等ポテンシャル線) を描く
  16 solve Laplace(phi,v) =
  17   int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v))
  18   -int1d(Th,Gamma)(Vn*v);
  19 plot(phi,ps="contourpotential.eps",wait=1);
  20 
  21 // ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
  22 v1=dx(phi);
  23 v2=dy(phi);
  24 plot([v1,v2],ps="vectorfield.eps",wait=1);
  25 
  26 // 等ポテンシャル線とベクトル場を同時に描く
  27 plot([v1,v2],phi,ps="both.eps", wait=1);

• 6行目で領域$\Omega$ の境界$\Gamma$ を指定している。
• 8行目で、$\Gamma$ を$m$ 分割して、$\Gamma$ の囲む範囲を三角形分割して、 それを Th に代入している。
• 11行目、有限要素空間 $V_h$ を区分的1次関数の空間と定義している。 この辺についてはこの講義では説明を省略する (知りたい人は有限要素法のテキストを読んで下さい)。
• 12行目、$\phi$ と対応する試験関数 $v$ を Vh の要素とする。
• 16〜18行目、弱形式の定義。ここを修正する必要が生じる可能性は低い。