FreeFEM 実例集

有限要素法による汎用の偏微分方程式解法プログラムとしては、 FreeFEMは極めて実用的なツールですが、構文を説明した文法書だけでは理解しにくいと思いますので、いくつかの実例を集めました。これらのほとんどは FreeFEMの配布ソースに説明抜きで入っていますし、一部は新しい FreeFEM++ の構文で、FreeFem++ のマニュアルに収録されています。

1. 弾性膜

弾性膜 Ω が平面の硬い枠 Γ に固定され、f(x)*dx の圧力を膜の全体に与えるものとすれば、膜の垂直方向の変移 ψ(x) は、ラプラス方程式

 (1.1) -△ψ = f in Ω

と、ディリクレ境界条件

 (1.2) ψ = 0 at Γ

を満たします。２次元問題だけを考えますので、ラプラス演算子 △ は

 (1.3) △ψ = (∂/∂x)^2(ψ) + (∂/∂y)^2(ψ)

として定義されています。

1.1. 弾性薄膜の歪み

長半径 a = 2、短半径 b = 1の楕円膜を考えることにして、表面に加わる力を 1 とすれば、下記のようなプログラムになります。


border(1,0,2*pi,40)	/* The number 1 boundary */
 begin
  x:=2*cos(t);
  y:=sin(t);
 end;
buildmesh(800);

solve(v) begin
 onbdy(1) v=0;
 pde(v) -laplace(v)=1
end;
plot(v);

三角分割はキーワード buildmesh で行われます。これはデローニー(Delaunay)テストに基づいたもので、まず境界点のみで三角形に分割し、その後、それらの辺を再分割することで内部の点を追加しますが、指示された最大値(このプログラムでは 800)を上回ることはありません。三角分解の細かさは辺の細かさで制御できます。

偏微分方程式はこれらの三角形を使って１次の有限要素に離散化され、頂点数に等しい連立方程式に変換された後、ガウスの LU 分解で解かれます。

2. 静電場

2.1. キャパシタ

容器 C0 内部に導体 Ci (i=1,2,..N) が存在し、個々の導体の電位を ψi、容器 C0 の電位を 0 とすれば、

 (2.1) △ψ = 0 in Ω,  ψ|Γ = g

が成り立ちますから、円筒状シールドケースに封入した、板状電極によるキャパシタの電位分布なら、次のようになります。


/* a circle of radius 5 centered at (0,0) */
border(1,0,2*pi,60) { x:=5*cos(t); y:=5*sin(t) };
/* The rectangle on the right */
border(2,0,1,4) { x:=1+t; y:=3 };
border(2,0,1,24) { x:=2; y:=3-6*t };
border(2,0,1,4) { x:=2-t; y:=-3 };
border(2,0,1,24) { x:=1; y:=-3+6*t };
/* The rectangle on the left */
border(3,0,1,4) { x:=-2+t; y:=3 };
border(3,0,1,24) { x:=-1; y:=3-6*t };
border(3,0,1,4) { x:=-1-t; y:=-3 };
border(3,0,1,24) { x:=-2; y:=-3+6*t };
buildmesh(800);

/* Boundary conditions and PDE */
solve(v) {
  onbdy(1) v = 0;
  onbdy(2) v = 1;
  onbdy(3) v = -1;
  pde(v) -laplace(v) = 0;
};
plot(v);

3. 熱伝導

3.1. 長方形断面のロッド

周囲温度 te (t0) に固定したときの、ロッド内部の温度分布を求めます。

板の中央を通る垂直な対称面で考えると２次元問題になりますが、対流と放射による損失が t0 + t1 - 2*te に比例し、その比例定数を α とすれば、

 (3.1) -△u +α*u = 0 in Ω
 (3.2) u|Γ1 = t0, u|Γ2 = t1, ∂u/∂n|Γ3 = 0

ですから、ロッドの長さを 10 とすれば、下記のようになります。



border(2,0,1,4)   begin x:=10;    y:=t   end;
border(4,0,1,4)   begin x:=0;     y:=1-t end;
border(1,0,10,40) begin x:=t;     y:=0   end;
border(3,0,10,40) begin x:=10-t;  y:=1   end;
buildmesh(800);

t0:=10; t1:=100; te:=25;
alpha:=0.05*(t0+t1-2*te);

solve(v) {
 onbdy(4) v=t0;
 onbdy(2) v=t1;
 onbdy(1,3) dnu(v)=0;	/* the 2 horizontal edges */
 pde(v) -laplace(v)+id(v)*alpha=0;
};
plot(v);

3.2. 円形断面のロッド

ロッドの断面が円形の場合は円筒座標を使って、

 (3.3) -∂r(r∂ru) - ∂z(r∂zu) + rαu = 0 in Ω
 (3.4) u|Γ1 = t0, u|Γ2 = t1, ∂u/∂n|Γ3 = 0

ですから、


border(2,0,1,4)   begin x:=10;    y:=t   end;
border(4,0,1,4)   begin x:=0;     y:=1-t end;
border(1,0,10,40) begin x:=t;     y:=0   end;
border(3,0,10,40) begin x:=10-t;  y:=1   end;
buildmesh(800);

t0:=10; t1:=100; te:=25;
alpha=0.1*((t0+t1)/2-te);

solve(v) {
 onbdy(4) v=t0;
 onbdy(2) v=t1;
 onbdy(1,3) dnu(v)=0;
 pde(v) -laplace(v)*y+id(v)*alpha*y=0;
};
plot(v);

3.3. 対流と放射

前記の対流と放射のモデルは極度に単純化したものですが、対流を

 (3.5) ∂u/∂n = -b*(u - te)

で近似し、放射を黒体放射の４乗則

 (3.6) ∂u/∂n - c*((u + 273)^4 - (te + 273)^4) = 0

で近似して、これらの合計を

 (3.7) ∂u/∂n + b*(u - te) - c*((u + 273)^4 - (te + 273)^4) = 0

とすれば、はるかによい精度が得られます。

問題は非線形ですから、反復解法を使うことになって、

           m+1        m+1          m+1        m
  (3.8) ∂u/∂n + b*(u - te) - c*(u   - te)*(u + te + 546)
         m
     *((u + 273)^2 - (te + 273)^2) = 0

を解くことになります。a^4 - b^4 = (a - b)*(a + b)*(a^2 + b^2) に注意してください。


border(2,0,1,4)   begin x:=10;    y:=t   end;
border(4,0,1,4)   begin x:=0;     y:=1-t end;
border(1,0,10,40) begin x:=t;     y:=0   end;
border(3,0,10,40) begin x:=10-t;  y:=1   end;
buildmesh(800);

t0:=10; t1:=100; te:=25; tk:=273;
b=0.1; c=-5.0e-9;
u=(t0*(10-x)+t1*x)/10;
iter(10) {
 w=b-c*(u+te-2*tk)*((u-tk)^2+(te-tk)^2);
 solve(v) {
  onbdy(4) v=t0-te;
  onbdy(2) v=t1-te;
  onbdy(1,3) id(v)*w+dnu(v)=0;
  pde(v) -laplace(v)=0;
 };
 u=v+te;
 plot(u);
};

3.4. 一様でない板

長方形断面の問題に戻って、サーモスタットのように、２枚の板の張り合わせを考えます。熱伝導率 k が場所によって違いますから、

 (3.9) α*u - ▽・(k▽u) = 0

を解くことになって、下記のようになります。


h0:=1.0; h1:=1.0;
border(2,0,1,4)   begin x:=10;   y:=h0*t        end;
border(2,0,1,4)   begin x:=10;   y:=h0+h1*t     end;
border(4,0,1,4)   begin x:=0;    y:=h0+h1*(1-t) end;
border(4,0,1,4)   begin x:=0;    y:=h0*(1-t)    end;
border(1,0,10,40) begin x:=t;    y:=0           end;
border(3,0,10,40) begin x:=10-t; y:=h0+h1       end;
border(0,0,10,40) begin x:=t;    y:=h0          end;
buildmesh(800);

t0:=10; t1:=100; alpha:=1;
kappa=0.01+one(y>h0);
solve(v) {
 onbdy(4) v=t0;
 onbdy(2) v=t1;
 onbdy(1,3) dnu(v)=0;
 pde(v) -laplace(v)*kappa+id(v)*alpha=0;
};
plot(v);

4. 音響

静止状態の空気圧の振動は波動方程式

 (4.1) (∂/∂t)^2(u) - c^2*△u = f

に従います。

解が単音なら、u は u(x,t) = v(x)*exp(i*k*t) 型の解を持ち、v はヘルムホルツ方程式(Helmholtz's equation)

 (4.2) (∂/∂t)^2(u) - c^2*△ = f in Ω
 ∂v/∂n|Γ = 0

の解になります。ラプラシアンの前の符号が + であることに注意してください。プラス符号の場合は特定の係数で共振を生じ、病的なふるまいをします。

音速を c、周波数を k、音源の分布を f とすれば、次の例はコンサートホールのステージ中央から指数的に減衰する音源からの垂直断面に於ける音波の伝搬状態を表します。


border(1,0,1,20) begin x:=5;     y:=1+2*t end;
border(1,0,1,20) begin x:=5-2*t; y:=3     end;
border(1,0,1,20) begin x:=3-2*t; y:=3-2*t end;
border(1,0,1,10) begin x:=1-t;   y:=1     end;
border(1,0,1,10) begin x:=0;     y:=1-t   end;
border(1,0,1,10) begin x:=t;     y:=0     end;
border(1,0,1,20) begin x:=1+4*t; y:=t     end;

buildmesh(1000);

f=exp(-10*((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)));
solve(v) {
 onbdy(1) dnu(v)=0;
 pde(v) laplace(v)+id(v)*0.8=f;
};
plot(v);

4. 渦なし流れ

4.1. 非圧縮流れ

粘性のない非圧縮性流体の渦なし流れの流速は

 (4.1) u = (∂ψ/∂y,-∂ψ/∂x)
       ψ は △ψ = 0 の解

になりますが、壁面に垂直方向の流速は 0 ですから、ψ = 0 が境界条件になります。

今、x-y 軸の両方に対称なノズル

  y = +-b*(x/a)^2*(2-(x/a)^2)+1  (-a<x<a)

内部の流れを考えることにすれば、下記のようになります。


a:=6; b:=1;
border(1,0,1,8) begin x:=-a; y:=1+b-2*(1+b)*t end;
border(2,0,1,16) {
 x:=-a+2*a*t;
 y:=-1-b*(x/a)*(x/a)*(2-(x/a)*(x/a))
};
border(3,0,1,8) begin x:=a; y:=-1-b+2*(1+b)*t end;
border(4,0,1,16) {
 x:=a-2*a*t;
 y:=1+b*(x/a)*(x/a)*(2-(x/a)*(x/a))
};
buildmesh(800);
solve(psi) {
 onbdy(1) psi=y;
 onbdy(2) psi=-1-b;
 onbdy(3) dnu(psi)=0;	/* free outflow */
 onbdy(4) psi=1+b;
 pde(psi) -laplace(psi)=0;
};
plot(psi);

4.2. 翼と揚力

翼の場合は境界に於ける ψ の値が未知数になります。一様な流れに於かれた翼 S を考えることにして、無限遠を大きな円 Γ∞ で表現します。前記と同様に

 (4.2) △ψ = 0 in Ω
       ψ|s = c, ψ|Γ∞ = u∞1*y - u∞2*x

を解くことになりますが、Ω は流体が占める空間、u∞ は無限遠の流速、 c は ∂nψ が翼の後端で連続になるように決めなければならない定数 (Kutta-Joukowski 条件) で、揚力は c に比例します。

c を決めるために、重ね合わせ法を使うことにします。すべての方程式は線形ですから、 ψ0 を c = 0 に対する解、ψ1 を c = 1 で無限遠の速度が 0 に対する解とすれば、

 (4.3) ψc = ψ0 + c*ψ1

と分解することができて、この２つの解から c を決めることができます。

古典的な翼形 NACA0012 の上表面を表す方程式は

 (4.4) y = 0.17735*sqrt(x)-0.075597*x-0.212836*x^2+0.17363*x^3-0.06254*x^4

ですが、迎角 α を tan(α)=1 とすれば、 Γ2 を翼形、Γ1 を無限遠近似として、

 (4.5) -△ψ = 0 in Ω, ψ|Γ1 = y - 0.1*x, ψ|Γ2 = c

を解くことになりますから、

 (4.6) -△ψ0 = 0 in Ω, ψ0|Γ1 = y - 0.1*x, ψ0|Γ2 = c
 (4.7) -△ψ1 = 0 in Ω, ψ0|Γ1 = 0, ψ0|Γ2 = 1

の２つを解くことで c が得られます。

ψ = ψ0 + c*ψ1 から c を求めるには、翼の後端 P = (1,0) に於いて ∂nψ が連続になるという条件を利用して、P+ を P からわずか δ だけ法線方向に上、 P- をわずか下とすれば、

 ∂ψ ≒ -(ψ(P+) - ψ(P))/δ
  n

から、

∂nψ は (ψ0|P+ + c*(ψ1|P+ - 1)) + (ψ0|P- + c*(ψ1|P- - 1)) を δ で割った値になって、

 (4.8) c = -(ψ0|P+ + ψ0|P-)/(ψ1|P+ + ψ1|P- - 2)

になりますが、プログラムでは、下記のように計算しています。

 (4.9) c = (ψ0(0.99,0.01)+ψ0(0.99,-0.01)/(ψ1(0.99,0.01)+ψ1(0.99,-0.01)-2)

border(1,0,2*pi,40) { x:=5*cos(t); y:=5*sin(t) };
border(2,0,2,71) {
 if(t<=1)then {	/* upper surface */
  x:=t; y:=0.17735*sqrt(t)-0.075597*t
	-0.212836*(t^2)+0.17363*(t^3)-0.06254*(t^4);
 }
 else {	/* lower surface */
  x:=2-t; y:=-(0.17735*sqrt(2-t)-0.075597*(2-t)-0.212836*((2-t)^2)
	+0.17363*((2-t)^3)-0.06254*(2-t)^4)
 };
};
buildmesh(1800);

solve(psi0) {
 onbdy(1) psi0=y-0.1*x;
 onbdy(2) psi0=0;
 pde(psi0) -laplace(psi0)=0;
};

solve(psi1,-1) {	/* -1 => matrix already factored */
 onbdy(1) psi1=0;
 onbdy(2) psi1=1;
 pde(psi1) -laplace(psi1)=0;
};
beta:=-(psi0(0.99,0.01)+psi0(0.99,-0.01))
	/(psi1(0.99,0.01)+psi1(0.99,-0.01)-2);
psi=beta*psi1+psi0;
cp=-dx(psi)^2-dy(psi)^2;	/* Presure */
plot(cp);

5. 輸送拡散問題

5.1. 数値解法

u(x) を速度場として、汚染物質 c(x,t) が初期値 c0(x) ないし内部供給源 f や境界供給源 g として存在するケースを考えます。c の拡散方程式は、υ を拡散係数、流体を非圧縮性とすれば、

 (5.1) ∂c/∂t + u・▽c - υ△c = f in Ω×]0,T[

圧縮性流体の場合は、u・▽c を ▽・(u*c) に書き換えなければなりません。

下記の初期条件と境界条件で、c が決まります。

 (5.2) c(x,0) = c0(x) ∀x∈Ω,
 (5.3) c|Γ1 = cΓ, ∂c/∂n|Γ2 = g with Γ1∪Γ2 = Γ, Γ1∩Γ2 = 0

(5.1) を時間領域で離散化する手法はいくつかありますが、Crank-Nicolson 解法が精度と安定性を備えています。

                 m+1  m               m+1  m            m+1  m
 (5.4) (1/δt)*(c   -c ]+(1/2)*u・▽(c   +c )-(υ/2)△(c   +c ) = f

補助変数

       m+1  m
 w = (c   +c )/2

を使えば、(5.4) は次のように簡略化でき、

                 1   m                          m+1         m
 (5.5) (2/δt)*(w - c ] + u・▽w - υ△w = f,  c   = 2*w - c

下記のようなプログラムが書けます。

 iter(...)
 begin...
  solve(w) {
   ...
   pde(w) id(w)*2/dt+dx(w)*u1+dy(w)*u2
             -laplace(w)*nu=f+c*2/dt;
  };
  c=2*w-c;
 end

υ が極めて小さい場合、この方法は発振しますが、その理由は次のとおりです。まず、υ=0 のケースも物理的に可能で(例えば、c が混合不能な水中の油滴の場合)、この場合、輸送拡散方程式は純粋は輸送方程式

 (5.6) ∂c/∂t + u・▽c = f in Ω×]0,T[

になりますが、(5.3) との組合せのは適していませんから、次の条件を使います。

 (5.7) c(x,0) = c0(x) ∀x∈Ω,
 (5.8) c|Γ1 = cΓ, ∂c/∂n|Γ2 = g with Γ1∪Γ2 = Γ-, Γ1∩Γ2 = 0

Γ- は流体が Ω に流れ込む境界の一部です。

 Γ- = {x∈Γ:u(x).n(x)<0}

n(x) は点 x に於ける外向き法線です。

これは、後続情報の必要性を示唆しますから、中心差分でなく、例えば次のような手法を使うことになります。

                                m+1     m                          m+1
 (5.9) ∂c+u・▽c = (Dc/Dt)(x,t) ≒ (c   (x)-c *(x-u(x)δt))/δt at t = t
        t

このアイデアを実現するために Gfem で用意されているのが convect 関数で、

 convect(c,u1,u2,dt)

の呼出で、c から c(x-u(x)*dt) を作るようになっています。u1, u2 は速度の２つの成分、dt は時間間隔です。

かくて、この問題のもう一つの解法は次のようになります。

 iter(...)
 begin...
  solve(w){...
   pde(w) id(w)*2/dt-laplace(w)*nu
     =f+c/dt+convect(c,u1,u2,dt)/dt
  };
  c=2*w-c;
 end

以上を理解するために、渦なしで粘性のない流体がながれるノズルを考えます。ノズルは２つに分かれていて、最初に (0.8*a,0) に存在した汚染物質がどう拡散するかを知りたいとすれば、次のようになります。


a:=6; b:=1; c:=0.5;
border(1,0,1,8) { x:=-a; y:=1+b-2*(1+b)*t };
border(2,0,1,46) { x:=-a+2*a*t; y:=-1-b*(x/a)*(x/a)*(3-2*abs(x)/a) };
border(3,0,1,8) { x:=a; y:=-1-b+(1+b)*t };
border(4,0,1,20) { x:=a-a*t; y:=0 };
border(4,0,pi,8) { x:=-c*sin(t)/2; y:=c/2-c*cos(t)/2 };
border(4,0,1,30) { x:=a*t; y:=c };
border(3,0,1,8) { x:=a; y:=c+(1+b-c)*t };
border(5,0,1,55) { x:=a-2*a*t; y:=1+b*(x/a)*(x/a)*(3-2*abs(x)/a) };
buildmesh(800);

nowait;
solve(psi) {
 onbdy(1) psi=y;
 onbdy(2) psi=-1-b;
 onbdy(3) dnu(psi)=0;
 onbdy(4) psi=0;
 onbdy(5) psi=1+b;
 pde(psi) -laplace(psi)=0;
};

u1=dy(psi); u2=-dx(psi);
dt=0.4; plot(psi); nu=0.01; j:=1;
yy=(x+a*0.8)^2+y^2;
phi=exp(max((-0.4*yy),-10));

iter(10) {
 f=convect(phi,u1,u2,dt)/dt;
 solve(w,j) {
  onbdy(1,2,4,5) w=0;
  onbdy(3) dnu(w)=0;
  pde(w) id(w)*2/dt-laplace(w)*nu=phi/dt+f;
 };
 phi=2*w-phi; plot(phi); j:=-1;
};

5.2. 回転するベル

古典的な輸送問題の１つが「回転するベル」(rotating bell)です。Ω が中心 0 の１枚の円盤で、速度

 (5.10) u1 = y, u2 = -x

で回転しているものとして、下記の問題を考えます。

 (5.11) ∂c + u・▽c = 0 in Ω
         t
 (5.12) c(t=0) = c0

正確な解は c(x,y,t) = c0(X(t),Y(t)) で、(X,Y) は (x,y) を原点を中心として θ=t だけ回転したものと等しく、３次元的に見てベル型に見えるような c0 を考えると、c も全く同じ形で、単に回転したかどうかの違いになります。

このゲームは T=2*π つまり完全に１回転するまで方程式を解いて、最終解を初期解と比較し、その一致を確かめることです。


twopi:=2*pi;
border(1,0,twopi,100) { x:=cos(t); y:=sin(t) };
buildmesh(800);
nowait;
r2=(x-0.3)^2+(y-0.3)^2;
v=exp(-10*r2);	/* initial condition */
plot(v);
dt:=0.3; u1=y; u2=-x; i:=1;

iter(20) {
 solve(v,i) {
  pde(v) id(v)=convect(v,u1,u2,dt);
 };
 i:=-1; plot(v);
};

注: この例では、プログラムの最初に「precise」命令を追加して補間関数の精度を高めることで、より良い結果が得られます。

6. ナビエ・ストークス(NAVIER-STOKES)方程式

非圧縮性の粘性流体は下記の方程式を満足します。

 (6.1) ∂u+u・▽u+▽p-υ△u = 0, ▽・u = 0 in Ω×]0,T[
        t
 (6.2) u|t=0 = u0, u|Γ = uΓ

考えられるアルゴリズムとしては、

                 m+1  m  m     m    m
 (6.3) (1/δt)*(u   -u0*X )+▽p-υ△u = 0, u|Γ=uΓ
          m+1        m  m      m+1
 (6.4) -△   = -▽・u0*X ,  ∂p   = 0
                             n

(5.9) の手法で

 ∂u + u・▽u
  t

を近似すれば、u0*X(x) = u(x - u(x)*δt になります。

注: 無限遠を円で近似する外部問題では、u,∂np を必要とする内向きの境界と、 ∂nu, p を必要とする外向き境界を区別しなければなりません。(Projection algorithm が p に対する境界条件を許します。)

6.1. 階段状水路の例

流入路が流出路より狭い、階段を後ろ向きに見たような水路を考えます。この形状は循環ゾーンを生みますから、そこを正確に把握しなければなりませんが、これは三角分割が十分に細かいか、良く流れに沿ったものであるときのみ可能です。


/* NAVIER STOKES EQUATIONS */
/* Poor but better than none algorithm*/
border(1,0,1,6)  begin x:=0;          y:=1-t    end;
border(2,0,1,15) begin x:=2*t;        y:=0      end;
border(2,0,1,10) begin x:=2;          y:=-t     end;
border(2,0,1,20) begin x:=2+3*t;      y:=-1     end;
border(2,0,1,35) begin x:=5+15*t;     y:=-1     end;
border(3,0,1,10) begin x:=20;         y:=-1+2*t end;
border(4,0,1,35) begin x:=5+15*(1-t); y:=1      end;
border(4,0,1,40) begin x:=5*(1-t);    y:=1      end;
buildmesh(1500);

nowait;
nu=0.001; dt:=0.4;	/* initial value */
u=y*(1-y)*one(y>0); v=0; p=2*nu*x*one(y>0);
un=u; vn=v; i:=1; j:=2; k:=3;
iter(180) {
 f=convect(un,u,v,dt); g=convect(vn,u,v,dt);
 solve(u,i) {	/* Horzontal velocity */
  onbdy(1) u=y*(1-y);
  onbdy(2,4) u=0;
  onbdy(3) dnu(u)=0;
  pde(u) id(u)/dt-laplace(u)*nu=f/dt-dx(p);
 }; plot(u);

 solve(v,j) {	/* Vertical velocity */
  onbdy(1,2,3,4) v=0;
  pde(v) id(v)/dt-laplace(v)*nu=g/dt-dy(p);
 };

 solve(p,k) {	/* Presure */
  onbdy(1,2,4) dnu(p)=0;
  onbdy(3) p=0;
  pde(p) -laplace(p)=-(dx(f)+dy(g))/dt;
 };
 un=u; vn=v; i:=-1; j:=-2; k:=-3;
};

7. 連立方程式の例

7.1. 弾性

固体に力を加えると変形し、当初 (x,y,z) にあった点はある時間の後に (X,Y,Z) に移動しますが、ベクトル U = (u,v,w) = (X-x,Y-y,Z-z) は歪み(displacement)と呼ばれています。歪みが小さくて、応力テンソル σ と歪みテンソル ε の間にフックの法則(Hooke's law)が成立すれば、下記のようになります。

 (7.1) σ = λ*δ*▽・U + μ
        ij      ij         ij
       δ  = 1  (i = j)
         ij
       ε  = (1/2)*(∂U /∂x + ∂U /∂x )
                       i    j     j    i

λ と μ は固体の機械的特性を表す定数ですが、よく知られた定数、ヤング率 (Young's modulus) E とポアソン比(Poisson.s ratio) υ との間に次の関係があります。

 (7.2) μ = E/(1 + υ), λ = E*υ/(1 + υ)/(1 - 2*υ)

7.2. ラメのシステム

長方形断面を持ち z 軸方向に於かれた梁を考えます。断面 Ω の歪みの x, y 成分 u(x,y) が軸に直交する応力 f を受けているとすると、下記が成り立ます。

 (7.3) -λ*△u - μ*▽(▽・u) = f in Ω

λ と μ は前記のラメ定数(Lame coefficients)です。

中心に割け目のある長方形断面の梁の水平面に垂直な引張力が加えられた場合を考えます。垂直面は自由に動けます。対称性を考慮すれば全体の 1/4 を解析すれば済みますので、割け目、下面の残りの面、右側の垂直面、上側の水平面、左側の垂直面をそれぞれ、 Γ1,Γ2,Γ3,Γ4,Γ5 とすれば、境界条件は下記のようになります。

 (7.4) λ*∂u/∂n+μ*(▽・u)*n = 0 on Γ1
       v = 0, λ*∂u/∂n+μ*▽・u = 0 on Γ2
       u = 0, λ*∂v/∂n+μ*▽・u = 0 on Γ3,Γ5
       u = 0 on 左の垂直面
                                     T
       λ*∂u/∂n+μ*(▽・u)*n = (0,1) on Γ4

u は u(u,v) の２つの成分を持つベクトルで、ラメの方程式(Lame's equation)は

 (7.5) -λ*△u - μ*(∂/∂x)^2(u) - μ*(∂^2/∂x/∂y)(v) = f1
       -λ*△v - μ*(∂/∂y)^2(v) - μ*(∂^2/∂x/∂y)(u) = f2

になります。

この２つの方程式は混合した微分によって密接に関連していますから、別々に反復解法を適用するのは危険で、Gfem が用意した連立方程式の解法を使うべきです。


/* Lame's system */
border(1,0,1,20) { x:=1-(t-1)^2; y:=0 };
border(2,1,2,20) { x:=(t-1)^2+1; y:=0 };
border(3,0,1,5) { x:=2; y:=t };
border(4,0,2,10) { x:=2-t; y:=1 };
border(5,0,1,5) { x:=0; y:=1-t };
buildmesh(500);

E:=2.15; sigma:=0.29; lambda:=E/(2*(1+sigma));
nu:=E*sigma/((1+sigma)*(1-2*sigma));
mu:=nu+sigma;

solve(u,v) {
 onbdy(3,5) u=0;
 onbdy(2) v=0;
 onbdy(4) dnu(v)=1;
 pde(u) -laplace(u)*mu-dxx(u)*nu-dxy(v)*nu=0;
 pde(v) -laplace(v)*mu-dyx(u)*nu-dyy(v)*nu=0;
};
plotmesh;
x=x+0.1*u;
y=y+0.1*v;
plotmesh;

7.3. ツィエンキビッチの問題

古典的な例[Zienkiewicz(1971)]として、穴のある無限遠まで広がる板に水平方向の一定な力 σx = -1 が加わる場合を考えます。

計算に際しては大きな正方形の板を考えることにして、さらに対称性を利用して 1/4 の領域だけを考えることにします。対称条件を穴を除く座標軸についてのみ適用すれば、下記のプログラムが得られます。


border(1,0,1.8,10) { x:=0; y:=2-t };
border(2,pi/2,0,10) { x:=0.2*cos(t); y:=0.2*sin(t) };
border(3,0,1.8,20) { x:=0.2+t; y:=0 };
border(4,0,2,10) { x:=2; y:=t };
border(5,0,2,20) { x:=2-t; y:=2 };
buildmesh(800);

E:=1;
nu:=0.29;
mu:=E/(1+nu);
lambda:=E*nu/((1+nu)*(1-2*nu));

solve(u,v) {
 onbdy(1) u=0;
 onbdy(1) v=0;
 pde(u) -laplace(u)*lambda-dxx(u)*mu-dxy(v)*mu=-1;
 pde(v) -laplace(v)*lambda-dyx(u)*mu-dyy(v)*mu=0;
};
plot(u); plot(v);

7.4. ストークスの問題

３番目の座標軸について不変の２次元流れでレイノルズ数(Reynolds number)が小さいとき、流れ関数 Ψ、速度 u = (u1,u2) と渦度 ω との間に下記の関係が成り立ちます。

 (7.6) u1 = ∂Ψ, u2 = -∂Ψ
             y           x
 (7.7) -△Ψ = ω, -△ω = 0

境界条件は

 (7.8) Ψ = ΨΓ, ∂ψ/∂n = g on Γ

ですが、ΨΓ と g は縁に於ける速度から決まります。Ψ に対して２つの条件があり、 ω に対する条件が存在しないことに注意してください。このシステムは実際には４次ですが、(7.8) を

  Ψ + ε*∂ψ/∂n = Ψ  ,  ∂Ψ/∂n = g on Γ
                      Γ
  ここに、ε<<1

で置換すると２次のシステムにすることができます。

下記の例では上面が speed(1,0) で移動するキャビティになっています。


border(1,0,1,20) { x:=t; y:=0 };
border(1,0,1,20) { x:=1; y:=t };
border(2,0,1,20) { x:=1-t; y:=1 };
border(1,0,1,20) { x:=0; y:=1-t };
buildmesh(1000);

eps:=0.00001;

solve(p,o) {
 onbdy(1,2) id(p)/eps + dnu(o) = 0;
 onbdy(1) dnu(p) = 0;
 onbdy(2) dnu(p) = 1;
 pde(o) - laplace(o) = 0;
 pde(p) id(o) + laplace(p) = 0;
};

plot(o); plot(p);

8. 複素数の例

電子レンジの熱源は電磁波による分子の加速です。単一周波数の平面波はヘルムホルツ方程式(Helmholtz's equation)

 (8.1) βυ + △υ = 0

で表現できます。

正方形のオーブンを考えることにして、電磁波は上面の一部から放射されることにすれば、 v = 0 の境界 Γ1 と υ = sin((π*(y-c)/(c-d)) の境界 Γ2 = [c,d] の２つになります。

料理の対象を b とすれば、熱源は υ^2 に比例し、平衡条件から、

 (8.2) -△θ = υ^2*Ib,  θ = 0
                          Γ

Ib は料理される物体内部で 1、それ以外では 0 になります。

下記のプログラムでは空中で β = 1/(1-I/2)、料理される物体内部で β = 2/(1-I/2) としています。(I = sqrt(-1))


border(1,0,1,20) { x:=t; y:=0 };
border(1,0,1,20) { x:=1; y:=t };
border(2,0,1,20) { x:=1-t; y:=1 };
border(1,0,1,20) { x:=0; y:=1-t };
buildmesh(1000);

eps:=0.00001;

solve(p,o) {
 onbdy(1,2) id(p)/eps + dnu(o) = 0;
 onbdy(1) dnu(p) = 0;
 onbdy(2) dnu(p) = 1;
 pde(o) - laplace(o) = 0;
 pde(p) id(o) + laplace(p) = 0;
};

plot(o); plot(p);