FreeFEM++の使い方

FreeFEM++は、フランスのパリ第6大学のJ.L.Lions研究所やINRIAで開発が進められている有限要素法のためのインタープリター環境とでも言うべきものである。日本では、広島国際大学の大塚厚二教授がプロジェクトに参加されていて、啓蒙活動を進められている。筆者も学会で隣に座られた教授に直接啓蒙されて、マニュアルまでいただいたという様な経緯がある。

1次元から3次元まで解くことが出来る環境であるが、2次元のマイクロストリップ線路の準TEMモード解析を例にして説明する。

このモデルを表現するのに、FreeFEM++では、まず領域の枠（border）を指定し、それを用いて2次元領域を定義してメッシュを作成させる。ここまでの計算メッシュを作成するまでのスクリプトが、以下である。

real b=0.65e-3, w=0.2e-3, t=0.035e-3; /* the hight, width and thickness of line */
real er=4.7, e0=8.854e-12, c0=2.998e8;
real lx=20.0e-3, ly=10.0e-3;
real p1, p2, p3, p4, C, C0, L, Z0;
real ax0, ax1, ay0, ay1; /* boundary */
real bx0, bx1, by0, by1; /* boundary of line */
ax0=0.0; ay0=0.0; ax1=lx; ay1=ly;
bx0=(lx-w)*0.5; by0=b; bx1=(lx+w)*0.5; by1=b+t;

int n = 200; /* initial mesh gen */
mesh Th = buildmesh( a1(n) + a2(n) + a3(n) + a4(n)
+ b1(-5) + b2(-3) + b3(-5) + b4(-3)
+ d1(5) + d2(5) ); plot(Th, wait=1, value=false);

buildmesh関数で、ボーダーから領域を指定している。この時、border変数に引数が割り当てられているが、これはそのボーダーで作られるメッシュの数を指定している。また、2次元領域を指定する際、ボーダーの方向に対し左側が領域に対応するようになっている。逆向きである場合、b1などであるように数にマイナス符号をつけることで表現できる。

最後にplot命令で図示される。この時のメッシュはかなり偏ったひどいものである事に注意していただきたい(図２)。後で改良する。

中心にマイクロストリップがあり、その下側は基板（誘電体）である。そして計算領域の外周として外側の枠がある。この内部で、ラプラス方程式で表されるポテンシャル方程式を解く問題に帰結することは、ものの本などに書いてあるおなじみのものである。

$$

\nabla \varepsilon \nabla \phi = 0

$$

境界条件として、$\phi=1\ \mbox{on b1,b2,b3,b4}$, $phi=0\ \mbox{on a1,a2,a3,a4)}である。これを表現したスクリプトが以下である。

fespace Ph(Th, P0);
fespace Vh(Th, P1);
Vh uh, vh, vx, vy;
Ph reg=region; /* the boundary of dielectrics */

plot(Th, reg, wait=1 );

int r1=reg(ay0, by0);
cout << "Region is " << r1 << endl;
Ph erc = er*(region==r1) + 1.0 * (region!=r1);
problem laplace(uh, vh) = int2d(Th) ( -erc*(dx(uh)*dx(vh)+dy(uh)*dy(vh)) )
+ on(1, 2, 3, 4, uh=0)
+ on(5, 6, 7, 8, uh=1); /* definition of governing eq */
laplace; /* solve */

ここで、on(1,2,3,4)とあるのは境界条件の指定で、数字は境界のラベルである。border命令順に番号が割り振られる。

先に表示されたメッシュは偏ったものなので、それによって得られた解も甘い良いものでは無い。そこで、adaptmesh命令を用いて解からそれに適したメッシュを作成するツールを使う。その際、残渣を指定することによって精度を上げることが出来る。使い方は以下の通り。

real error=0.1;
int i = 0;
for(i=0; i<5; i++) {
Th = adaptmesh(Th, uh, inquire=0, err=error); /* adaptive mesh */
error = error * 0.5;
laplace;
};

 

 

 最終的に得られた計算メッシュは図３にあるようにかなりなめらかなものとなっている。

個のメッシュ状で得られたポテンシャル分布のコンター表示が図４である。

最終的に線路のインピーダンスを求めるために、後処理も行うことが出来る。

ここで、電力$P$は、

$$P = VI = V^2/Z = E_p \cdot c_{eff},$$

で与えられる。ここで、$E_p, c_{eff}$はマイクロストリップラインの断面上でのエネルギーと、光速度を示す。また誘電体が無い場合のエネルギーと光速度、$E_0,c_0$は、

$$P_0 = E_0 \cdot c_0$$

なる。また、

$$c_0 = \frac{1}{\sqrt{\varepsilon_0 \mu_0}},$$

$$c_{eff} = \frac{1}{\sqrt{\varepsilon_0\varepsilon_r \mu_0}}$$

と表すとする。また、

$$E = \int_\Omega \varepsilon_{eff} \varepsilon_0 |\nabla \phi |^2 \ d\Omega,$$

でエネルギーは計算できる。これを踏まえ、特性インピーダンスと、有効誘電率$\varepsilon_{eff}$を後処理で求めることが出来る。

$$\varepsilon_{eff} = E / E_0, $$

$$Z = \frac{V}{E_{eff} c_{eff}}, $$

である。V=1であるので以下の様な手続きを得る。

vx=dx(uh); vy=dy(uh); /* calc. for charge */
p1=int1d(Th, b1)( erc*vy); p2=int1d(Th, b2)(-erc*vx);
p3=int1d(Th, b3)(-erc*vy); p4=int1d(Th, b4)( erc*vx);
C=e0*(p1+p2+p3+p4);
E = e0* int2d(Th)(erc*(dx(uh))*(dx(uh)) + erc*(dy(uh))*(dy(uh)));

erc = 1.0;
laplace;
vx=dx(uh); vy=dy(uh); /* calc. for charge */
p1=int1d(Th, b1)( vy); p2=int1d(Th, b2)(-vx);
p3=int1d(Th, b3)(-vy); p4=int1d(Th, b4)( vx);

C0 = e0*(p1+p2+p3+p4); L = 1/(c0^2*C0);
E0 = e0* int2d(Th)(erc*(dx(uh))*(dx(uh)) + erc*(dy(uh))*(dy(uh)));
Ep = E / E0; /* effective permittivity ratio */
cr = c0 / sqrt(Ep);
Z = 1.0/(E*cr);

cout << "C = " << C << " L = " << L << endl;
cout << "Impedance = " << Z << endl;
cout << "Effective Permittivity = " << Ep << endl;