FreeFem++で電流分布を求める

複雑な形状の電磁界解析は解析的に解けないものですから有限要素解析により数値計算で求めることが一般的です．
フリーな有限要素解析ソフトウェアとしてFreeFem++が有名です．ドキュメントには有限要素解析のやり方が詳細記述されているのですが，電流解析のサンプルがなかったので，使い方習得を兼ねて方法をまとめてみます．やっていることはFreeFem++のドキュメントにあるGetting startedそのままです．

お題

左図の電流分布および電界分布を求めます．灰色で示す金属電極が黄色で示す絶縁物の板の上に乗っているものが，導電率ρの一様な導電体の中に埋まっています．境界条件は，金属電極が電圧1V，無限遠が0Vとします．

基本の数式

電磁気の基本式は上図左の式で与えられます．ベクトルはシンボルの上にバーを引いて表しています．まず電流ベクトルJは電界ベクトルEと導電率ρで与えられます．また電圧ベクトルEは電圧Vのgradientで与えられます．ここで電荷保存則はJとρの式で表されますが，ここでは定常電流を扱うので，dρ/dtは0，つまりdivJ=0となります．これらから解くべき式はとラプラスの式を解けばよいです．

FreeFem++のシミュレーションファイルを書く

どんなシミュレーションをするかをファイルに書かねばなりません．以下の順序で記述していきます:

領域の定義
メッシュの切り方を書く
解きたい問題を式で表す
シミュレーションの実行と結果表示

領域の定義

まずは有限要素解析を行う領域を定義します．サンプルとして半径2の円の中で辺々2の正方形以外の領域を指定してみます．C++ライクな書式で，こんな感じで書きます:

border C1(t=0, 2 *pi) {x=2*cos(t), y=2*sin(t)}; //従属変数tを使って円を表す．指定する領域は，ベクトル進行方向左手の領域．
 int C2 = 0; // ラベル付けのための適当な変数．値は何でもいい．
border C21(t=-1, 1){ x=t; y = 1; label=C2}; // 正方形の辺々をそれぞれ指定する．4つの辺をラベルC2で結びつけて1つの正方形を表す．
...略...
border C24(t=-1, 1){ x=-1; y =t; label=C2}; // 最後の辺

境界線を従属変数で表して領域を分けていきます．複数の境界線を結合するにはラベルを使います．ここのポイントは，境界線のどちらが指定する領域になるかです．境界線を指定する従属変数が増加する方向に向かったとき，左手の領域がシミュレーションを行う領域になります．例えば，サンプルの円C1だと，従属変数t = 0で座標(1, 0)から始まりt=0.5*piで座標(0,1)に移動しますから，指定する領域はその左手方向，原点を含む円形領域，になります．

メッシュを切る

有限要素解析は，シミュレーションをする領域をメッシュに切り，解くべき式をそのメッシュの頂点に対する連立方程式で表して数値解析するものです．このためメッシュの切り方で結果の精度が違ってしまいます．このメッシュを切るテクニックというかノウハウは知らないのですが，FreeFem++は自動的にメッシュを切ってくれるのが特徴みたいです．では，メッシュを切ってみます．

mesh Th = buildmesh(C0(50) + C21(10) + C22(10) ...略... C24(10));

変数 mesh Th は関数buildmeshで作ります．それぞれの境界にいくつ頂点を置くかを指定するだけで，あとは自動的にメッシュが作製されるみたいです．

問題を式に表す

解くべき式はのラプラスの式です．FreeFem++には積分形で書かないといけないので，そのあたりを復習します．
ラプラスの式はポアソンの式の特別な形です．ポアソンの式に対して，バウンダリで0になる任意の関数vをかけて領域Ωで積分しても両辺は同じです．
.この式をグリーンの公式で変形すると
となるそうです...む?
直感的にいまいち分からないので，数学的にむちゃくちゃだと思うのですが，こんなものかなと求めてみると．スカラ関数fとgに対して，∇(f･g)=f･∇g+∇f･g．ここでfをv, gを∇uとすれば，∇･(v･∇u)=∇u･∇v+v･Δu．この両辺を体積積分すると，左辺はベクトル関数(v･∇u)の発散の体積積分．発散の体積積分は，ガウスの発散定理から(v･∇u)の境界上での線積分になる．ところがvは境界上で0だから左辺はゼロ．したがって∇u･∇v+v･Δu=0よりv･Δu=-∇u･∇v．
さて問題に戻ります．最初の式J=ρ(x,y)E，E=-∇V,∇･J=0から，∇･J=∇･(ρE)=-∇･∇ρV=-ΔρV=0とラプラスの式が求まります．均一な導電体であればρは定数ですから，見慣れたΔV=0になりますが，今は絶縁板があったりと導電率が不均質ですからρが消えません．
これを書くとこうなります．

solve a(u, v) = int2d(Th)(kappa*(dx(u)*dx(v)+dy(u)*dy(v)))
+on(C0, u= 0) + on(C1, u= 1);

変数 kappa は導電率を表します．kappa は以下のように表しました．読みやすさのために，表記中に"絶縁板の厚さ"および"絶縁板の幅"という変数名で書いています．シミュレーションファイル中では，これらの変数名は適当な英数文字にしています．

Vh kappa = 0.001
-0.001 * (x < 0.5 * 絶縁板の幅) * (x > -0.5 * 絶縁板の幅) * (y < -絶縁板の厚さ*0.5) * (y > 絶縁板の厚さ*0.5);

Cライクな表記のためか，数式中での条件判定はTrueが1，Falseが0に対応するようです．ちょうどユニット関数として使えます．これを使い，絶縁板のある部分ではkappaが0になるように組んでみました．

ソースコード

全ソースコードです．

// Electric field and current flow simulation
// Jul/15/2008 Akihiro Uehara
// Detail of this simulation is described in 
// freefem++ document, chapter "learning by example", "heat change".

// Variable declarations.
real eht = 0.001;  // Eelctrode height [m]
real ewd = 0.0005; // Electrode width  [m]

real sht = 0.000005;  // substrate height [m]
real swd = 0.006;     // substrate width  [m]
real sofs= 0.000010;  // distance between the electrode and the substrate (to prevent boarder overrup)

int C1 = 98, C2 = 97, C3 =96; // could be anything

// simulation resgion 
border C0(t=0, 2*pi) { x=0.05*cos(t); y=0.05*sin(t);}

// electrode region
border C11(t=-0.5, 0.5){x=t    * ewd; y=eht;        label=C1;}
border C12(t= 1,     0){x=0.5  * ewd; y=t    * eht; label=C1;}
border C13(t= 0.5,-0.5){x=t    * ewd; y=0;          label=C1;}
border C14(t=0,      1){x=-0.5 * ewd; y=t    * eht; label=C1;}

// substrate
border C21(t=0.5, -0.5){x=t    * swd; y=-sofs;            label=C2;}
border C22(t=0,     -1){x=-0.5 * swd; y=t    * sht -sofs; label=C2;}
border C23(t=-0.5, 0.5){x=t    * swd; y=-1 *   sht -sofs; label=C2;}
border C24(t=-1,     0){x=0.5  * swd; y=t    * sht -sofs; label=C2;}

// target cell
border RGC(t=0, 2*pi) { x=0.0001*cos(t); y=0.0001*sin(t) + 0.0012;}

// plot border
plot(C0(50)
+C11(10)+C12(10)+C13(10)+C14(10)
+C21(-10)+C22(-10)+C23(-10)+C24(-10)
+RGC(10),
wait=1
);

// Build mesh
mesh Th = buildmesh(C0(50)
+C11(10)+C12(10)+C13(10)+C14(10)
+C21(-10)+C22(-10)+C23(-10)+C24(-10)
+RGC(10)
);
plot(Th, wait=1);

// solve
fespace Vh(Th, P1); 
Vh u, v;
Vh kappa = 0.001
-0.001 * (x < 0.5 * swd) * (x > -0.5 * swd) * (y < -sofs) * (y > -sofs - sht)
;

solve a(u, v) = int2d(Th)(kappa*(dx(u)*dx(v)+dy(u)*dy(v)))
+on(C0, u= 0) + on(C1, u= 1);
plot(u, wait=true, value=true, fill=true);

シミュレーション

では上記のファイルを使いシミュレーションしてみます．まずFreeFem++をインストールしてから，FreeFem++ GUIを立ち上げます．そして上記のファイルを読み込みます．
マウスをクリックするたびに，waitにより表示がいったん停止します．マウスをクリックするたびに，ボーダー，生成したメッシュ，そして以下の図のようなシミュレーション結果が表示されます．

シミュレーション結果のファイルダンプ

さて，設計の基礎データを求めるために，このような数値シミュレーションをしていると思います．そのような場合は，たいていはグラフ表示だけでははなくて，なんらかの設計値を求めるなどの，シミュレーション結果の解析が不可欠になります．
ここでは，手ごろな使い慣れた解析ツールを使うために，FreeFemのシミュレーション結果をテキストダンプしてみます．

まず先ほど求めた電界をファイル"f.txt"にダンプするには，前回のコードに続けて次のコードを書きます．これは，まず頂点Vhで電位の関数fを定義します．次に出力ファイル fileを開いて，関数fの内容を全てダンプします．

Vh f = u;
{
ofstream file("f.txt");
file << f[] << endl;
}

次に電界の大きさをダンプしてみます．次のコードを，前回のコードに追加します．電界のダンプは電位の傾き大きさで与えられます．まず関数efを定義します．次に先ほどと同じようにファイルダンプします．欲しいデータは電極周辺の電界大きさだけですので， if( (Th[i][0].x ^2 + Th[i][0].y ^2) < 0.6e-4 ) という条件式を入れて，画面中心の半径600umの領域のデータだけをダンプしています．ファイル出力はC++の書式と同じです．Th[i][0].x および Th[i][0].y は，それぞれxおよびy座標を表します．次の efh[][Vh(i,0)] が電界の大きさです. ダンプした後に，確認のために電界の大きさをプロットしています．

// electric field
func ef = sqrt(dx(u)^2 + dy(u)^2);
Vh efh = ef; 

// dump data (document frefem++doc.pdf pp.23)
{ 
  ofstream ff("graph.txt");
  for (int i=0;i<Th.nt;i++)
  { 
      if( (Th[i][0].x ^2 + Th[i][0].y ^2) < 0.6e-4 )
      {
//    for(int j=0; j<3; j++)
//      ff<< Th[i][j].x << " " << Th[i][j].y << " "<<efh[][Vh(i,j)] << endl;
//    ff<< Th[i][0].x << " " << Th[i][0].y << " "<<efh[][Vh(i,0)] << "\n\n\n";
      ff<< Th[i][0].x << " " << Th[i][0].y << " "<<efh[][Vh(i,0)] << endl;
      }
  }
}

plot(Th, efh, wait=true, bb=[[0.002,0.002],[-0.002,-0.002]], value=true, fill=false);