FreeFEM++ 例題9-6-1 非圧縮NS方程式
uは流速、pは圧力である。まず、キャビティ流れを扱おう。この問題はボックスの三方を固定壁で囲まれ、上側の側面がu=1で移動することで、領域内の流体が対流を起こすという問題である。境界条件として、上側で流速u=1、その他の壁面ではすべりなし条件(no slip条件:u,v=0)を与える。
 
label=1
label=2
label=3
label=4
(u=1,v=0)
ペナルティ法: 非圧縮の条件に微少量εと圧力pを掛けたものを加える。
離散化すると
real  nu=1./100.;
real dt=0.1;
real alpha=1/dt;
Xh up1,up2;
problem  NS ([u1,u2,p],[v1,v2,q],solver=Crout,init=i) =
    int2d(Th)(
             alpha*( u1*v1 + u2*v2 -convect([up1,up2],-dt,up1)*v1 -convect([up1,up2],-dt,up2)*v2)
            + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1) +  dx(u2)*dx(v2) + dy(u2)*dy(v2) )
            - p*q*(0.000001)  - p*dx(v1)- p*dy(v2)
            - dx(u1)*q- dy(u2)*q
           )
   + on(3,u1=1,u2=0) + on(1,2,4,u1=0,u2=0) ;     ////境界条件
for (i=0;i<=20;i++)
 { up1=u1;  up2=u2;
   NS;
   if ( !(i % 10))
    plot(coef=0.2,cmm=" [u1,u2] et p  ",p,[u1,u2]);
 } ;