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

/****************************************************************************/
/* This file is part of FreeFEM.                                            */
/*                                                                          */
/* FreeFEM is free software: you can redistribute it and/or modify          */
/* it under the terms of the GNU Lesser General Public License as           */
/* published by the Free Software Foundation, either version 3 of           */
/* the License, or (at your option) any later version.                      */
/*                                                                          */
/* FreeFEM is distributed in the hope that it will be useful,               */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of           */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            */
/* GNU Lesser General Public License for more details.                      */
/*                                                                          */
/* You should have received a copy of the GNU Lesser General Public License */
/* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>.          */
/****************************************************************************/

// Parameters
int n = 3; // mesh quality

// Mesh
mesh Th = square(10*n, 10*n);

// Fespace
fespace Uh(Th, P1b);
Uh u, v;
Uh uu, vv;
fespace Ph(Th, P1);
Ph p, pp;

// Problem
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))
    -1e-10*p*pp
  )
  + on(1, 2, 4, u=0, v=0)
  + on(3, u=1, v=0)
  ;

// Plot
plot([u,v],p,wait=1);

version 4.9 では、 /usr/local/ff++/share/FreeFEM/4.9/examples/examples/stokes.edp に置いてある。

version 4.12 では /Applications/FreeFem++.app/Contents/ff-4.12/share/FreeFEM/4.12/examples/examples/stokes.edp においてある。

cp /usr/local/ff++/share/FreeFEM/4.9/examples/examples/stokes.edp .

あるいは

cp /Applications/FreeFem++.app/Contents/ff-4.12/share/FreeFEM/4.12/example/stokes.edp .

とカレント・ディレクトリにコピーしておいて

FreeFem++ stokes.edp

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

と同値である。

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

図 8: stokes-p1b.edp の結果

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

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

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

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

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