D. FreeFem++ のマニュアル§3.9 から: Stokes方程式の cavity flow


// stokes.edp (in the manual, chapter 3)
int n=3;
mesh Th=square(10*n,10*n);
plot(Th, wait=1, ps="stokes-p1b-mesh.eps");
fespace Uh(Th,P1b); Uh u,v,uu,vv;
fespace Ph(Th,P1); Ph p,pp;
solve stokes([u,v,p],[uu,vv,pp]) =
int2d(Th)(dx(u)*dx(uu)+dy(u)*dy(uu) + dx(v)*dx(vv)+ dy(v)*dy(vv)
+ dx(p)*uu + dy(p)*vv + pp*(dx(u)+dy(v)))
+ on(1,2,4,u=0,v=0) + on(3,u=1,v=0);
plot([u,v],p,wait=1,ps="stokes-p1b-solution.eps");

$\displaystyle \Vector{u}_h=\begin{pmatrix}u_1  u_2 \end{pmatrix} \leftarrow
...
... \texttt{pp},\quad
x_1 \leftarrow \texttt{x},\quad
x_2 \leftarrow \texttt{y}
$

とすると、

      $\displaystyle \dint_{\Omega_h} \left( \frac{\rd u_1}{\rd x_1}\frac{\rd v_1}{\rd...
...\rd v_2}{\rd x_1} +\frac{\rd u_2}{\rd x_2}\frac{\rd v_2}{\rd x_2} \right) \DxDy$
      $\displaystyle + \dint_{\Omega_h} \left( \frac{\rd p_h}{\rd x_1}v_1+\frac{\rd p_...
...h} q_h \left( \frac{\rd u_1}{\rd x_1}+\frac{\rd u_2}{\rd x_2} \right) \DxDy =0,$

すなわち

$\displaystyle a(\Vector{u}_h,\Vector{v}_h)+(\nabla p_h,\Vector{v}_h)
+(q_h,\nabla\cdot\Vector{u_h})=0$   $\displaystyle \mbox{($\forall(\Vector{v}_h,q_h)$)}$$\displaystyle .$

これは

$\displaystyle \left\{
\begin{array}{l}
a(\Vector{u}_h,\Vector{v}_h)+(\nabla p...
...,\nabla\cdot\Vector{u_h})=0
\quad\mbox{($\forall q_h$)}
\end{array} \right.
$

と同値である。

P1b/P1 を採用した結果は、図 7 となる。

図 7: stokes-p1b.edp の結果
\includegraphics[width=10cm]{stokes-p1b/stokes-p1b-solution.eps}

P1/P1 にすると上手く計算できないことは常識とされているが、 自分でプログラムを書く場合、なかなかそこまで試す気にはなれない。 しかし FreeFem++ でそれをするのは簡単である (図 8)。

図 8: もしも P1/P1 とすると…なんだ、これは?!
\includegraphics[width=10cm]{stokes-p1b/stokes-p1p1-solution.eps}

P2/P1 は、良く知られているようにきちんと解ける (図 9)。

図 9: P2/P1 (分割は粗くした)
\includegraphics[width=10cm]{stokes-p1b/stokes-p2p1-solution.eps}

他にも、悪い冗談のようだが P2/P2 などが気軽に試せる。

当り前かもしれないが、弱形式を少しいじっただけで結果がコロコロ変わる。

桂田 祐史
2017-12-05