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");

$$\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}$$

とすると、

$$\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$$

$$+ \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,$$

すなわち

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

これは

$$\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.$$

と同値である。

第2項が $-(p_h,\nabla\cdot\Vector{v}_h)$ でないところが目につく。 そうするためには

int2d(Th)(dx(u)*dx(uu)+dy(u)*dy(uu) + dx(v)*dx(vv)+ dy(v)*dy(vv)
-p*(dx(uu)+dy(vv))+ dy(p)*vv + pp*(dx(u)+dy(v)))
+ on(1,2,4,u=0,v=0) + on(3,u=1,v=0);

とするわけだろう。

とにかく P1b/P1 を採用した結果は、 図 6 となる。

図 6: stokes-p1b.edp の結果

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

図 7: もしも P1/P1 とすると…なんだ、これは？！

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

図 8: P2/P1 (分割は粗くした)

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

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