ここでは、支配方程式（２）を元に、GetDPでこの問題をどのように表記して、解いていくのかを紹介する。

GetDPの文法

GetDPでは、弱形式と境界条件で定義された問題を独自のスクリプト言語で記述していく。 この記述スクリプトの構成（オブジェクトと呼ばれている）は大きく分けて、10の項目に分かれる。

  • Group
  • Function
  • Jacobian
  • Integration
  • Constraint
  • FunctionSpace
  • Formulation
  • Resolution
  • PostProcessing
  • PostOperation

Groupでは、計算領域の名称を定義する。Functionでは、計算領域内で使われる物理量の空間・時間分布を定義する。 Constraintでは、境界条件や初期条件を定義する。FunctionSpaceでは要素のタイプや基底関数などを定義する。 Formulationでは、支配方程式を記述する。Resolutionでは、求解のためのシークエンスを定義する。 PostProcessingでは、後処理のための定義が行われる。PostOperationでPostProcessingで定義された後処理を実行し出力する。 Jacobian、Integrationでは、それぞれ幾何学的な変換、要素の積分方式を定義する。

各オブジェクトの詳細はOneLabサイトで提供されている、GetDPマニュアルもしくはソースコード;-)を参照していただきたい。

次節から、electrode.proの各オブジェクトの具体例で説明していく。

Group

Group {
  D = Region[1000];
  FAR  = Region[2000];
  ROD1 = Region[2001];
  ROD2 = Region[2002];
}

electrode.geoでPhysical Regionを定義した番号に対し、変数に定義していく。点数名は自由なのでわかりやすい名称にするのが良い。 Physical Regionでは体積領域には1000番台を、境界領域には2000番台をつけていた。その領域は、計算領域と、境界条件を与えるための領域で、後で利用するために名称が必要なわけである。 ここで、DはDomain、計算領域をPhysical Regionで定義した1000番に指定する。

次に、領域の境界条件の対象である、領域外部の境界をFARとして2000を当てる。また、境界条件として定電圧を印加される 計算領域D内にある二つの電極を、ROD1、ROD2としてそれぞれ2001、2002を当てる。 この問題は単純なので、これで領域の定義、つまりGroup Objectの設定は以上である。

 

Function

ここでは、領域内の物理量などの関数を定義する。

Function {
  eps[] = 1;
  rho[] = 1;
}

関数の[]内にはGroupで定義したRegion変数名を入れて指定することも出来る。

 eps[D] = 1;

省略すると全定義領域となる。

Jacobian

Jacobian {
  { Name V;
    Case {
      { Region All; Jacobian Vol; }
    }
  }
}

内部定義されているJacobian、ここではVol、に対し領域を当て、それに対して名称をつけている。 これは後で、FormulationやPostProcessingなどで使われる。

Integration

Integration {
  { Name I;
    Case {
      { Type Gauss;
        Case {
          { GeoElement Point; NumberOfPoints  1; }
          { GeoElement Line; NumberOfPoints  3; }
          { GeoElement Triangle; NumberOfPoints  3; }
          { GeoElement Quadrangle; NumberOfPoints  4; }
          { GeoElement Tetrahedron; NumberOfPoints  4; }
          { GeoElement Hexahedron; NumberOfPoints  6; }
          { GeoElement Prism; NumberOfPoints  9; }
          { GeoElement Pyramid; NumberOfPoints  8; }
        }
      }
    }
  }
}

要素の形状に対し、使われる積分方式を定義するオブジェクトである。 一般にはこの例にある様にGauss積分を使うので変更する必要は無い。

 

Constraint

Constraint {
  { Name BC; Type Assign;
    Case {
      { Region FAR ; Value 0.; }
      { Region ROD1; Value -1.; }
      { Region ROD2; Value  1.; }
    }
  }
}

Constraintオブジェクトでは、この例にあるように、境界条件を設定することが出来る。 ここでは、電位の境界条件として、領域の外側は０、内部の電極ROD1、ROD2にそれぞれ−１，１を与えている。 そしてこの境界条件に対し、名前BCを当てて、後で使えるようにしてある。

FunctionSpace

FunctionSpace {
  { Name E1; Type Form0;
    BasisFunction {
      { Name sn; NameOfCoef vn; Function BF_Node; Support D; Entity NodesOf[All]; }
    }
    Constraint {
      { NameOfCoef vn; EntityType NodesOf;
        NameOfConstraint BC; }
    }
  }
}

FunctionSpaceオブジェクトでは、物理量（ここでは電位であるvn）に対し、どのような要素タイプで基底関数を何にするかを指定できる。 定義されるFunctionSpaceに対し、この例ではE1という名前をつけている。 この例では、E1において、vnは組み込み基底関数BF_Nodeであり、計算メッシュ領域D内のすべて(All)のNodeに当てはめられている、となっている。 さらに、Constraintとして、先に定義されたBCが設定されている。

 

Formulation

Formulation {
  { Name Poisson; Type FemEquation;
    Quantity {
      { Name v; Type Local; NameOfSpace E1; }
    }
    Equation {
      Galerkin { [ eps[] * Dof{d v}, {d v} ] ; In D; Jacobian V; Integration I; }
      Galerkin { [ rho[], {v} ] ; In D; Jacobian V; Integration I; }
    }
  }
}

この例では、vを既にFunctionSpaceで定義されたE1と宣言し、Equationで支配方程式(2)を記述している。 Galerkin文で弱形式の各項を表記し、各項の和がゼロになるように定義する。 ここで、Galerkin文にある、{d v}は電位vに対する外微分を意味する。 この場合は、Grad vとも書くことが出来る。

 

Resolution

Resolution {
  { Name EleSta_v;
    System {
      { Name Sys_Ele; NameOfFormulation Poisson; }
    }
    Operation {
      Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
    }
  }
}

Resolutionオブジェクトでは、求解のための計算を実行させる命令を記述している。

 

PostProcessing/PostOperation

PostProcessing {
  { Name EleSta_v; NameOfFormulation Poisson;
    PostQuantity {
      { Name v; Value { Term { [ {v} ]; In D; Jacobian V; } } }
      { Name e; Value { Term { [ -{d v} ]; In D; Jacobian V; } } }
    }
  }
}

PostOperation {
  { Name EleSta_v; NameOfPostProcessing EleSta_v;
    Operation {
      Print[ e, OnElementsOf D, File "EleSta_v_e.pos" ];
      Print[ v, OnElementsOf D, File "MagSta_v_v.pos" ];
    }
  }
}

ここでは後処理のための処理の定義（PostProcessing）と実行（PostOperation）を記述している。 ここでは、直接の解である電位vとそれに対する電場eを定義し、それらをEleSta_v_e.posとEleSta_v_v.posという指定されたファイルに出力させている。 電場はResolutionでは求められていないが、電位の勾配を求めることによって得ることが出来る。 PostProcessingでも微分などの演算子を処理させることが出来る。

以上で、２つの電極を内在させる矩形領域での２次元電位問題を解くための記述が説明された。

 

解析結果

Figure 5: 自動的に張られた計算メッシュ
Figure 6: 解析結果
[6a:計算された電位分布]
 [6b:計算された電場分布]

図5に自動的に張られた計算格子と、図6(a)、6(b)に解析された結果である電位と電場の分布を示す。 シンプルさと軽さに留意したモデルであるので、計算メッシュはかなり粗いが、見たところもっともらしい解析結果が得られている。 より細かいメッシュにしたい場合は、Gmshで定義したlcの値を小さくすれば良い。 ２次元なのでlc=0.01程にしてもそれほど計算時間はかからない。

 

FaLang translation system by Faboba