DefineConstant
DefineConstant[ Flag_AnalysisType = {2, Choices{ 0="Steady-state", 2="Transient"}, Name "Parameters/0Type of analysis"} ];
今回は、Gmshのメニューで静解析、過渡解析を選択できるようにする。そのためのメニューはこのようにDefineConstantオブジェクトで指定する。 このメニューを利用することにより、Flag_AnalysisTypeという変数に0:Steady-state、2:Transientが代入され、後でそれを使って解析を選択できる。
Constraint
ここでは、境界条件を設定する。
Constraint { { Name Temperature ; Case { { Region l_rod ; Type Assign; Value kleft ; } { Region r_rod ; Type Assign; Value kright ; } { Region Surf_Domain ; Type Assign; Value kdomain ; } If(Flag_AnalysisType != 0) { Region Domain ; Type Init; Value kinit ; } EndIf } } }
ここで、l_rod, r_rod, Surf_Domainは境界領域で、それぞれ左、右の定温棒、領域外部壁を示す。 また、時間発展解の場合(Flag_AnalysisType=2)、空間内の温度初期値がkinitに設定される。 ここで、kleft,kright,kdomain,kinitはそれぞれ定数が既に代入されている(Function)。
FunctionSpace
スカラー値なので、ノード毎に物理量が存在するので、FunctionとしてBF_Nodeを指定する。 またConstraintとして、前節で設定したTemperatureを指定する。
FunctionSpace { { Name Hgrad_T; Type Form0; BasisFunction { { Name sn; NameOfCoef Tn; Function BF_Node; Support Domain; Entity NodesOf[All]; } } Constraint { { NameOfCoef Tn; EntityType NodesOf ; NameOfConstraint Temperature; } } } }
Formulation
ここでは、既に得られている熱伝導方程式の弱形式をそのままGetDPスクリプトで記述する。 また、一応熱源の項、qVolも含めておいた。 この関数に適当な値を代入することによって、熱源が存在する問題も解くことが出来る。 この例では、qVol[] = 0となっている。
Formulation { { Name The_T ; Type FemEquation; Quantity { { Name T; Type Local; NameOfSpace Hgrad_T; } } Equation { Galerkin { [ k[] * Dof{d T} , {d T} ]; In Domain; Integration I1; Jacobian JVol; } Galerkin { DtDof [ rhoc[] * Dof{T} , {T} ]; In Domain; Integration I1; Jacobian JVol; } Galerkin { [ -qVol[] , {T} ]; In Domain; Integration I1; Jacobian JVol; } } } }
Resolution
今回のモデルは、静解析と過渡解析をメニューで選択できるようにした。 解法の指定は、このResolution Objectで行う。 If ... EndIf文がGetDPスクリプトの中で使うことが出来る。 Formulationの中で、DtDofという時間微分項が存在するが、その項は、解き方の指定をここですることによって、静解析にも過渡解析にも対応できるようになっている。
Resolution { { Name analysis; System { { Name T; NameOfFormulation The_T; } } Operation { If(Flag_AnalysisType == 0) // steady state Generate[T] ; Solve T ; SaveSolution T; EndIf If(Flag_AnalysisType == 2) // transient linear fast InitSolution[T] ; SaveSolution[T] ; GenerateSeparate[T] ; TimeLoopTheta [t0, t1, dt, 1.0] { Update[T] ; SolveAgain[T] ; } EndIf } } }
過渡解析の場合、初期状態を設定しなければならない。そのため
InitSolution
が呼ばれている。 また、時間微分項に対して分離できる行列を、
GenerateSeparete[T]
を使うことによって、保持し、後で
Update[T]
を呼ぶことによって行列の再生成を省略できる。
TimeLoopTheta[ t0,t1,dt,$\theta$]
は時間積分の台形公式[Trapezoidal Scheme]を指定している。 ここでは、$\theta=1.0$、つまりCrank-Nickolson公式の陰解法を指定している。 これらの命令の詳しくはGetDPマニュアル[2]を参照して欲しい。ただ難解である。