熱および電気弾性連成の例 (単純な非線形連成問題，モデルオブジェクト，汎用アセンブリ，解析および可視化)¶
この例では，変位フィールド，温度フィールド，および電位フィールドの非線形連成のマルチフィジックス問題の簡単な例を紹介します．また，C++ ライブラリと各インターフェイスの使用法を比較することも目的としています．対応するデモファイルは， GetFEM のテストディレクトリにあります( tests/ ， interface/tests/python ， interface/scr/scilab/demos ， interface/tests/matlab )．
問題の設定¶
\(\Omega \subset \rm I\hspace{-0.15em}R^2\) は，2次元板の基準配置です( ここ のジオメトリを参照してください)，ここで厚さ \(\varepsilon\) は外部力，電位および加熱に寄与します．ここで， \(\theta : \Omega \rightarrow \rm I\hspace{-0.15em}R\) は( °C 内の)温度フィールド， \(V : \Omega \rightarrow \rm I\hspace{-0.15em}R\) は電位フィールドです．また， \(u : \Omega \rightarrow \rm I\hspace{-0.15em}R^2\) は膜変位場です．
熱問題¶
プレートの表裏面は，熱伝達係数 \(D\) の (20 °C の) 空気と熱交換しており，プレートの側面は断熱であると想定されています．
熱方程式 \(\theta\) と境界条件は次のように書くことができます．
\(\kappa\) は熱電導率， \(T_0\) は空気の温度， \(\partial \Omega\) は，領域 \(\Omega\) の境界であり \(n\) は \(\partial \Omega\) 上の \(\Omega\) への外向きの単位法線ベクトルです．
項 \(\sigma|\nabla V|^2\) は，Joule加熱項に対応する非線形連成項であり， \(\sigma\) は電気伝導率です．
電位問題¶
プレートの左右の側面の間には \(0.1V\) の電位差があるとします．他の面は電気的に絶縁されていると考えます．電位の方程式は次の通りです．
ここでも \(\sigma\) は電気伝導率です．さらに， \(\sigma\) は次のように温度に依存すると考えます．
ここで \(T_0\) は参照温度 (ここでは空気温度) であり， \(\rho_0\) は \(T_0\) 上の抵抗温度係数であり，\(\alpha\) は第2抵抗温度係数です．
変形問題¶
右側面に加わる力の下でプレートの微小変形を考慮し, プレートの加熱によって影響を受けるとします．変位 \(u\) は，次の(線形化された弾性)問題の解になります．
ここで \(F\) は右横方向の境界に加えられる力密度であり， \(\bar{\sigma}(u)\) は次のように定義されたCauchy応力テンソルです．
\(\bar{\varepsilon}(u) = (\nabla u + (\nabla u)^T)/2\) は，線形化されたひずみテンソルであり， \(I\) は2次単位テンソルで， \(\lambda^*, \mu\) は Lamé 係数であり次のように定義されます．
\(E\) は材料のヤング係数であり \(\nu\) はポアソン比です．
項 \(\beta(T_0-\theta) I\) は，熱膨張の項に対応しています．ここで \(\beta = \alpha_{th} E/(1-2\nu)\) で， \(\alpha_{th}\) は熱膨張係数です．
弱定式化¶
方程式の連成系方程式の弱定式化は重要なステップです．有限要素の定式化は弱定式化 (Galerkin近似) に基づいており，弱定式化は追加される項を表現する唯一の方法であるため，これはきわめて重要なステップです．
各偏微分方程式の弱定式化は，未知数に対応する試験関数を式に掛けることにより得られます．未知数はDirichlet条件を満たし条件は均質であるとします．その後 \(\Omega\) 領域上の積分と(Greenの式を使用した)部分のいくつかの積分を実行します．偏微分方程式系の弱定式化は次の通りです．
ここで \(\delta_{\theta}, \delta_V, \delta_u\) は，それぞれ， \(\theta, V, u\) に対応する試験関数です， \(\Gamma_N\) は力密度 \(F\) がかかる右の境界を示します．そして \(\bar{\sigma}:\bar{\varepsilon}\) は2次テンソル同士のFrobeniusスカラー積です．
C++ での実装とインターフェイス¶
次に，問題を近似する GetFEM の使用法についての詳細なプレゼンテーションを行います．C++，Python，Scilab と Matlab プログラムを同時に構築します．Matlab および Scilab プログラムでは，オブジェクト指向コマンドを使用しません (使用する方法は GetFEM OOコマンド を参照してください)
初期化¶
まず，C++ では，モデルオブジェクト，ジェネリックアセンブリ，線形代数インターフェイス (Gmm++) ，実験用メッシュ，およびエクスポート機能に対して，特定の数のヘッダーをインクルードしなければなりません．Python の場合，これは簡単です． GetFEM をグローバルにインポートするだけです(numpy もインポートする必要があります)．Scilab の場合，最初にライブラリを Scilab コンソールにロードします (ここでは説明しません)，Matlab の場合，すべての GetFEM 変数をクリアする gf_workspace('clear all') 以外，何も必要ありません．
C++ |
#include "getfem/getfem_model_solvers.h"
#include "getfem/getfem_export.h"
#include "gmm/gmm.h"
#include "getfem/getfem_mesher.h"
#include "getfem/getfem_generic_assembly.h"
using bgeot::size_type;
using bgeot::base_node;
using bgeot::base_small_vector;
typedef getfem::model_real_plain_vector plain_vector;
int main(void) {
|
Python |
import getfem as gf
import numpy as np
|
Scilab |
gf_workspace('clear all'); // In case of multiple runs
|
Matlab |
gf_workspace('clear all'); % In case of multiple runs
|
モデルのパラメータ¶
ここで，問題のさまざまな物理パラメータおよび数値パラメータを定義しましょう．スクリプト言語 (Python，Scilab と Matlab) には違いはありません．
C++ |
double epsilon = 1.; // Thickness of the plate (cm)
double E = 21E6; // Young Modulus (N/cm^2)
double nu = 0.3; // Poisson ratio
double clambda = E*nu/((1+nu)*(1-2*nu));
double cmu = E/(2*(1+nu));
double clambdastar = 2*clambda*cmu/(clambda+2*cmu);
double F = 100E2; // Force density at the right boundary (N/cm^2)
double kappa = 4.; // Thermal conductivity (W/(cm K))
double D = 10.; // Heat transfer coefficient (W/(K cm^2))
double air_temp = 20;// Temperature of the air in oC.
double alpha_th = 16.6E-6; // Thermal expansion coefficient (/K)
double T0 = 20.; // Reference temperature in oC
double rho_0 = 1.754E-8; // Resistance temperature coefficient at T0
double alpha = 0.0039; // Second resistance temperature coefficient
double h = 2. // Approximate mesh size
bgeot::dim_type elements_degree = 2; // Degree of the finite element methods
|
Scripts |
epsilon = 1.; E = 21E6; nu = 0.3;
clambda = E*nu/((1+nu)*(1-2*nu));
cmu = E/(2*(1+nu));
clambdastar = 2*clambda*cmu/(clambda+2*cmu);
F = 100E2; kappa = 4.; D = 10;
air_temp = 20; alpha_th = 16.6E-6;
T0 = 20; rho_0 = 1.754E-8;
alpha = 0.0039;
h = 2; elements_degree = 2;
|
メッシュ生成¶
GetFEM には，ここで説明するいくつかの制約のあるメッシュ機能があります．ここではそれらを使用するつもりです．しかし，得られたメッシュの品質と適合性は保証していません．そのため， GetFEM のメッシュ作成機能を使用する場合はメッシュを確認することをお勧めします．また，外部メッシャ (例えばGiD または Gmsh) を使用してインポートすることもできます( メッシュの保存とロード をご覧ください)．
領域のジオメトリは，3つの円形の穴がある長方形であるとします ( 得られたメッシュ． を参照)．ジオメトリは，いくつかの初等幾何と union/setminus 操作によって記述します ( src/getfem/getfem_mesher.h ファイルを参照してください)．以下では， h はメッシュサイズを表し， 2 はメッシュの次数を表します (これは変換が次数2であり，曲線エッジを使用することを意味します)．
C++ |
getfem::mesh mesh;
getfem::pmesher_signed_distance
mo1 = getfem::new_mesher_rectangle(base_node(0., 0.), base_node(100., 25.)),
mo2 = getfem::new_mesher_ball(base_node(25., 12.5), 8.),
mo3 = getfem::new_mesher_ball(base_node(50., 12.5), 8.),
mo4 = getfem::new_mesher_ball(base_node(75., 12.5), 8.),
mo5 = getfem::new_mesher_union(mo2, mo3, mo4),
mo = getfem::new_mesher_setminus(mo1, mo5);
std::vector<getfem::base_node> fixed;
getfem::build_mesh(mesh, mo, h, fixed, 2, -2);
|
Python |
mo1 = gf.MesherObject('rectangle', [0., 0.], [100., 25.])
mo2 = gf.MesherObject('ball', [25., 12.5], 8.)
mo3 = gf.MesherObject('ball', [50., 12.5], 8.)
mo4 = gf.MesherObject('ball', [75., 12.5], 8.)
mo5 = gf.MesherObject('union', mo2, mo3, mo4)
mo = gf.MesherObject('set minus', mo1, mo5)
mesh = gf.Mesh('generate', mo, h, 2)
|
Scilab |
mo1 = gf_mesher_object('rectangle', [0 0], [100 25]);
mo2 = gf_mesher_object('ball', [25 12.5], 8);
mo3 = gf_mesher_object('ball', [50 12.5], 8);
mo4 = gf_mesher_object('ball', [75 12.5], 8);
mo5 = gf_mesher_object('union', mo2, mo3, mo4);
mo = gf_mesher_object('set minus', mo1, mo5);
mesh = gf_mesh('generate', mo, h, 2);
|
Matlab |
mo1 = gf_mesher_object('rectangle', [0 0], [100 25]);
mo2 = gf_mesher_object('ball', [25 12.5], 8);
mo3 = gf_mesher_object('ball', [50 12.5], 8);
mo4 = gf_mesher_object('ball', [75 12.5], 8);
mo5 = gf_mesher_object('union', mo2, mo3, mo4);
mo = gf_mesher_object('set minus', mo1, mo5);
mesh = gf_mesh('generate', mo, h, 2);
|
得られたメッシュ．¶
境界の選択¶
境界のそれぞれの部分には異なる境界条件を設定するため，境界のさまざまな部分には番号を付けます (穴には，熱と電気の絶縁と，応力自由境界条件が想定されます)．したがって，メッシュ上の要素面を選択し，メッシュ領域を定義する必要があります ( メッシュ領域 を参照)，1，2，3，4はそれぞれ右境界，左境界，上境界，下境界です．これらの境界番号は，モデルのブリックで使用されます．
C++ |
getfem::mesh_region border_faces;
getfem::outer_faces_of_mesh(mesh, border_faces);
getfem::mesh_region fb1
= getfem::select_faces_in_box(mesh, border_faces, base_node(1., 1.),
base_node(99., 24.));
getfem::mesh_region fb2
= getfem::select_faces_of_normal(mesh, border_faces,
base_small_vector( 1., 0.), 0.01);
getfem::mesh_region fb3
= getfem::select_faces_of_normal(mesh, border_faces,
base_small_vector(-1., 0.), 0.01);
getfem::mesh_region fb4
= getfem::select_faces_of_normal(mesh, border_faces,
base_small_vector(0., 1.), 0.01);
getfem::mesh_region fb5
= getfem::select_faces_of_normal(mesh, border_faces,
base_small_vector(0., -1.), 0.01);
size_type RIGHT_BOUND=1, LEFT_BOUND=2, TOP_BOUND=3, BOTTOM_BOUND=4;
mesh.region( RIGHT_BOUND) = getfem::mesh_region::subtract(fb2, fb1);
mesh.region( LEFT_BOUND) = getfem::mesh_region::subtract(fb3, fb1);
mesh.region( TOP_BOUND) = getfem::mesh_region::subtract(fb4, fb1);
mesh.region(BOTTOM_BOUND) = getfem::mesh_region::subtract(fb5, fb1);
|
Python |
fb1 = mesh.outer_faces_in_box([1., 1.], [99., 24.])
fb2 = mesh.outer_faces_with_direction([ 1., 0.], 0.01)
fb3 = mesh.outer_faces_with_direction([-1., 0.], 0.01)
fb4 = mesh.outer_faces_with_direction([0., 1.], 0.01)
fb5 = mesh.outer_faces_with_direction([0., -1.], 0.01)
RIGHT_BOUND=1; LEFT_BOUND=2; TOP_BOUND=3; BOTTOM_BOUND=4; HOLE_BOUND=5;
mesh.set_region( RIGHT_BOUND, fb2)
mesh.set_region( LEFT_BOUND, fb3)
mesh.set_region( TOP_BOUND, fb4)
mesh.set_region(BOTTOM_BOUND, fb5)
mesh.set_region( HOLE_BOUND, fb1)
mesh.region_subtract( RIGHT_BOUND, HOLE_BOUND)
mesh.region_subtract( LEFT_BOUND, HOLE_BOUND)
mesh.region_subtract( TOP_BOUND, HOLE_BOUND)
mesh.region_subtract(BOTTOM_BOUND, HOLE_BOUND)
|
Scilab |
fb1 = gf_mesh_get(mesh, 'outer faces in box', [1 1], [99 24]);
fb2 = gf_mesh_get(mesh, 'outer faces with direction', [ 1 0], 0.01);
fb3 = gf_mesh_get(mesh, 'outer faces with direction', [-1 0], 0.01);
fb4 = gf_mesh_get(mesh, 'outer faces with direction', [0 1], 0.01);
fb5 = gf_mesh_get(mesh, 'outer faces with direction', [0 -1], 0.01);
RIGHT_BOUND=1; LEFT_BOUND=2; TOP_BOUND=3; BOTTOM_BOUND=4; HOLE_BOUND=5;
gf_mesh_set(mesh, 'region', RIGHT_BOUND, fb2);
gf_mesh_set(mesh, 'region', LEFT_BOUND, fb3);
gf_mesh_set(mesh, 'region', TOP_BOUND, fb4);
gf_mesh_set(mesh, 'region', BOTTOM_BOUND, fb5);
gf_mesh_set(mesh, 'region', HOLE_BOUND, fb1);
gf_mesh_set(mesh, 'region subtract', RIGHT_BOUND, HOLE_BOUND);
gf_mesh_set(mesh, 'region subtract', LEFT_BOUND, HOLE_BOUND);
gf_mesh_set(mesh, 'region subtract', TOP_BOUND, HOLE_BOUND);
gf_mesh_set(mesh, 'region subtract', BOTTOM_BOUND, HOLE_BOUND);
|
Matlab |
fb1 = gf_mesh_get(mesh, 'outer faces in box', [1 1], [99 24]);
fb2 = gf_mesh_get(mesh, 'outer faces with direction', [ 1 0], 0.01);
fb3 = gf_mesh_get(mesh, 'outer faces with direction', [-1 0], 0.01);
fb4 = gf_mesh_get(mesh, 'outer faces with direction', [0 1], 0.01);
fb5 = gf_mesh_get(mesh, 'outer faces with direction', [0 -1], 0.01);
RIGHT_BOUND=1; LEFT_BOUND=2; TOP_BOUND=3; BOTTOM_BOUND=4; HOLE_BOUND=5;
gf_mesh_set(mesh, 'region', RIGHT_BOUND, fb2);
gf_mesh_set(mesh, 'region', LEFT_BOUND, fb3);
gf_mesh_set(mesh, 'region', TOP_BOUND, fb4);
gf_mesh_set(mesh, 'region', BOTTOM_BOUND, fb5);
gf_mesh_set(mesh, 'region', HOLE_BOUND, fb1);
gf_mesh_set(mesh, 'region subtract', RIGHT_BOUND, HOLE_BOUND);
gf_mesh_set(mesh, 'region subtract', LEFT_BOUND, HOLE_BOUND);
gf_mesh_set(mesh, 'region subtract', TOP_BOUND, HOLE_BOUND);
gf_mesh_set(mesh, 'region subtract', BOTTOM_BOUND, HOLE_BOUND);
|
メッシュの描画¶
メッシュをプレビューし，その妥当性を制御するために，次の手順を使用します．
C++ |
getfem::vtk_export exp("mesh.vtk", false);
exp.exporting(mesh);
exp.write_mesh();
// You can view the mesh for instance with
// mayavi2 -d mesh.vtk -f ExtractEdges -m Surface
|
Python |
mesh.export_to_vtk('mesh.vtk');
# You can view the mesh for instance with
# mayavi2 -d mesh.vtk -f ExtractEdges -m Surface
|
Scilab |
scf(1);
gf_plot_mesh(mesh, 'refine', 8, 'curved', 'on', 'regions', ...
[RIGHT_BOUND LEFT_BOUND TOP_BOUND BOTTOM_BOUND]);
title('Mesh');
sleep(1000);
|
Matlab |
gf_plot_mesh(mesh, 'refine', 8, 'curved', 'on', 'regions', ...
[RIGHT_BOUND LEFT_BOUND TOP_BOUND BOTTOM_BOUND]);
title('Mesh');
pause(1);
|
C++ および Python インターフェイスでは，外部グラフィカルポストプロセッサを使用する必要があります (たとえば，gmsh，Mayavi2，PyVista，または Paraview)．Scilab および Matlab インタフェースでは，内部プロット機能を使用することができます (結果は 得られたメッシュ． を参照)．
有限要素法と積分法の定義¶
3つの有限要素法を定義します．最初の1つは，変位フィールドを近似する mfu です．これはベクトルフィールドで C++ では次のように定義されます．
getfem::mesh_fem mfu(mesh, 2);
mfu.set_classical_finite_element(elements_degree);
ここで， 2 はベクトル場の次元を表します．2行目は，使用する有限要素を設定します． classical_finite_element は，連続したLagrange要素を意味し，elements_degree は 2 に設定されています．これは2次の (アイソパラメトリック) 要素を使用することを意味します．
GetFEM では，既存の有限要素法を幅広く選択肢できます． 付録A.有限要素法リスト を参照してください．しかし，実際には Lagrange 有限要素法が最も使用されています．
第2の有限要素法はスカラーで，温度場と電位場の両方を近似する mft です．単一の有限要素法で任意の数の有限要素変数を近似すると便利です．
3番目の有限要素法は不連続なスカラーLagrangeで，1変数の導関数を補間することができます (たとえば，Von Mises応力の補間)．
最後に定義するのは，積分法 mim です．GetFEM にデフォルトの積分法はありません．したがって，これは積分法を定義するためには必須です．もちろん，積分法の次数は，選択された有限要素法に好都合な積分を行うため，十分に選定しなければなりません．ここでは， elements_degree の2乗で十分です．
C++ |
getfem::mesh_fem mfu(mesh, 2);
mfu.set_classical_finite_element(elements_degree);
getfem::mesh_fem mft(mesh, 1);
mft.set_classical_finite_element(elements_degree);
getfem::mesh_fem mfvm(mesh, 1);
mfvm.set_classical_discontinuous_finite_element(elements_degree);
getfem::mesh_im mim(mesh);
mim.set_integration_method(2*elements_degree);
|
Python |
mfu = gf.MeshFem(mesh, 2)
mfu.set_classical_fem(elements_degree)
mft = gf.MeshFem(mesh, 1)
mft.set_classical_fem(elements_degree)
mfvm = gf.MeshFem(mesh, 1)
mfvm.set_classical_discontinuous_fem(elements_degree)
mim = gf.MeshIm(mesh, elements_degree*2)
|
Scilab |
mfu = gf_mesh_fem(mesh, 2);
gf_mesh_fem_set(mfu, 'classical fem', elements_degree);
mft = gf_mesh_fem(mesh, 1);
gf_mesh_fem_set(mft, 'classical fem', elements_degree);
mfvm = gf_mesh_fem(mesh, 1);
gf_mesh_fem_set(mfvm, 'classical discontinuous fem', elements_degree-1);
mim = gf_mesh_im(mesh, elements_degree*2);
|
Matlab |
mfu = gf_mesh_fem(mesh, 2);
gf_mesh_fem_set(mfu, 'classical fem', elements_degree);
mft = gf_mesh_fem(mesh, 1);
gf_mesh_fem_set(mft, 'classical fem', elements_degree);
mfvm = gf_mesh_fem(mesh, 1);
gf_mesh_fem_set(mfvm, 'classical discontinuous fem', elements_degree-1);
mim = gf_mesh_im(mesh, elements_degree*2);
|
モデルの定義¶
GetFEM のモデルオブジェクトは(未知)変数，データ，およびモデルブリックと呼ばれるものの集まりです．モデルブリックは，単一の変数か，複数の変数をリンクしているモデルの一部 (線形または非線形項) です．これらは，(接線) 線形システムの構築のために使用されます (詳細については modelオブジェクト を参照してください)．
直接構築を行い(接線) 線形システムを作成することも可能なため，モデルオブジェクトを使用することは厳密には必須ではありません．しかし，モデルオブジェクトは，モデルのほとんどの既往な部分が事前にプログラムされているので，モデルの素早い構築が可能です． 標準境界条件，標準偏微分方程式，制約を処理するための乗数の使用などが整備されています．さらに，いくつかのブリックは，標準ブリック (ジェネリックアセンブリブリック，陽なマトリックスブリックなど) の拡張が可能なように設計されています．そのため，モデルオブジェクトのフレームワークを使用することをお勧めします．
モデルには実数と複素数の2つのバージョンがあります．複素数モデルは複雑な線形システムを解くのに有利な特殊なアプリケーション (例えば，いくつかの電磁気学の問題など) のために予約されています．
計算される3つのフィールドに対応する3つの変数を持つ実際のモデルを宣言してみましょう ．
C++ |
getfem::model md;
md.add_fem_variable("u", mfu);
md.add_fem_variable("theta", mft);
md.add_fem_variable("V", mft);
|
Python |
md=gf.Model('real');
md.add_fem_variable('u', mfu)
md.add_fem_variable('theta', mft)
md.add_fem_variable('V', mft)
|
Scilab |
md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mfu);
gf_model_set(md, 'add fem variable', 'theta', mft);
gf_model_set(md, 'add fem variable', 'V', mft);
|
Matlab |
md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mfu);
gf_model_set(md, 'add fem variable', 'theta', mft);
gf_model_set(md, 'add fem variable', 'V', mft);
|
弾性膜変形問題¶
ここでは，弾性変形問題から始めましょう．以下の add_isotropic_linearized_elasticity_brick によって追加されている定義済みのブリックを使用します．対応する項は以下の通りです．
この追加を接線線形システムに対して行います．このモデルブリックを使用するために， Lamé 係数に対応するデータは，最初にモデルに追加する必要があります．ここでは， Lamé 係数は，領域に対して一定です．ただし，一部の非定数データを定義することもできます．また，この定義済みのブリックを使用する代わりに，GWFL，汎用弱形式言語項 add_linear_term(md mim, "lambda*(Div_u*Div_Test_u) + mu*((Grad_u + Grad_u'):Grad_Test_u)" を使用することもできます．
連成項
は定義済みのブリックはなく，GWFL項 add_linear_term(md mim, "beta*theta*Div_Test_u)" を直接使用します．GWFLの詳細については， 任意の項を計算する - 高水準の汎用的な構築手順 - 弱形式言語 (GWFL) を参照してください．基本的に，アセンブリ文字列は，各Gauss点で実行される最適化されたアセンブリ命令のリストでコンパイルされます．
以下のプログラムは，全体の弾性変形方程式を考慮に入れることができます．左側の境界に Dirichlet 条件を規定するために，既定のブリックを使用します．Dirichlet条件を定義するいくつかのオプションがあります( Dirichlet条件ブリック要素 を参照)．
C++ |
md.add_initialized_scalar_data("cmu", cmu);
md.add_initialized_scalar_data("clambdastar", clambdastar);
md.add_initialized_scalar_data("T0", T0);
getfem::add_isotropic_linearized_elasticity_brick
(md, mim, "u", "clambdastar", "cmu");
getfem::add_Dirichlet_condition_with_multipliers
(md, mim, "u", bgeot::dim_type(elements_degree-1), LEFT_BOUND);
md.add_initialized_fixed_size_data("Fdata", base_small_vector(F*epsilon,0.));
getfem::add_source_term_brick(md, mim, "u", "Fdata", RIGHT_BOUND);
md.add_initialized_scalar_data("beta", alpha_th*E/(1-2*nu));
getfem::add_linear_term(md, mim, "beta*(T0-theta)*Div_Test_u");
|
Python |
md.add_initialized_data('cmu', [cmu])
md.add_initialized_data('clambdastar', [clambdastar])
md.add_initialized_data('T0', [T0])
md.add_isotropic_linearized_elasticity_brick(mim, 'u', 'clambdastar', 'cmu')
md.add_Dirichlet_condition_with_multipliers(mim, 'u', elements_degree-1, LEFT_BOUND)
md.add_initialized_data('Fdata', [F*epsilon, 0])
md.add_source_term_brick(mim, 'u', 'Fdata', RIGHT_BOUND)
md.add_initialized_data('beta', [alpha_th*E/(1-2*nu)])
md.add_linear_term(mim, 'beta*(T0-theta)*Div_Test_u')
|
Scilab |
gf_model_set(md, 'add initialized data', 'cmu', [cmu]);
gf_model_set(md, 'add initialized data', 'clambdastar', [clambdastar]);
gf_model_set(md, 'add initialized data', 'T0', [T0]);
gf_model_set(md, 'add isotropic linearized elasticity brick', mim, 'u', 'clambdastar', 'cmu');
gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', elements_degree-1, LEFT_BOUND);
gf_model_set(md, 'add initialized data', 'Fdata', [F*epsilon, 0]);
gf_model_set(md, 'add source term brick', mim, 'u', 'Fdata', RIGHT_BOUND);
gf_model_set(md, 'add initialized data', 'beta', [alpha_th*E/(1-2*nu)]);
gf_model_set(md, 'add linear term', mim, 'beta*(T0-theta)*Div_Test_u');
|
Matlab |
gf_model_set(md, 'add initialized data', 'cmu', [cmu]);
gf_model_set(md, 'add initialized data', 'clambdastar', [clambdastar]);
gf_model_set(md, 'add initialized data', 'T0', [T0]);
gf_model_set(md, 'add isotropic linearized elasticity brick', mim, 'u', 'clambdastar', 'cmu');
gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', elements_degree-1, LEFT_BOUND);
gf_model_set(md, 'add initialized data', 'Fdata', [F*epsilon, 0]);
gf_model_set(md, 'add source term brick', mim, 'u', 'Fdata', RIGHT_BOUND);
gf_model_set(md, 'add initialized data', 'beta', [alpha_th*E/(1-2*nu)]);
gf_model_set(md, 'add linear term', mim, 'beta*(T0-theta)*Div_Test_u');
|
電位問題¶
同様に，以下のプログラムは，電位方程式を記述しています．電気伝導率 \(\sigma\) とGWFLの項の使用方法に注意してください．
C++ |
std::string sigmaeps = "(eps/(rho_0*(1+alpha*(theta-T0))))";
md.add_initialized_scalar_data("eps", epsilon);
md.add_initialized_scalar_data("rho_0", rho_0);
md.add_initialized_scalar_data("alpha", alpha);
getfem::add_nonlinear_term
(md, mim, sigmaeps+"*(Grad_V.Grad_Test_V)");
getfem::add_Dirichlet_condition_with_multipliers
(md, mim, "V", bgeot::dim_type(elements_degree-1), RIGHT_BOUND);
md.add_initialized_scalar_data("DdataV", 0.1);
getfem::add_Dirichlet_condition_with_multipliers
(md, mim, "V", bgeot::dim_type(elements_degree-1), LEFT_BOUND, "DdataV");
|
Python |
sigmaeps = '(eps/(rho_0*(1+alpha*(theta-T0))))'
md.add_initialized_data('eps', [epsilon])
md.add_initialized_data('rho_0', [rho_0])
md.add_initialized_data('alpha', [alpha])
md.add_nonlinear_term(mim, sigmaeps+'*(Grad_V.Grad_Test_V)')
md.add_Dirichlet_condition_with_multipliers(mim, 'V', elements_degree-1, RIGHT_BOUND)
md.add_initialized_data('DdataV', [0.1])
md.add_Dirichlet_condition_with_multipliers(mim, 'V', elements_degree-1, LEFT_BOUND, 'DdataV')
|
Scilab |
sigmaps = '(eps/(rho_0*(1+alpha*(theta-T0))))';
gf_model_set(md, 'add initialized data', 'eps', [epsilon]);
gf_model_set(md, 'add initialized data', 'rho_0', [rho_0]);
gf_model_set(md, 'add initialized data', 'alpha', [alpha]);
gf_model_set(md, 'add nonlinear term', mim, sigmaeps+'*(Grad_V.Grad_Test_V)');
gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'V', elements_degree-1, RIGHT_BOUND);
gf_model_set(md, 'add initialized data', 'DdataV', [0.1]);
gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'V', elements_degree-1, LEFT_BOUND, 'DdataV');
|
Matlab |
sigmaps = '(eps/(rho_0*(1+alpha*(theta-T0))))';
gf_model_set(md, 'add initialized data', 'eps', [epsilon]);
gf_model_set(md, 'add initialized data', 'rho_0', [rho_0]);
gf_model_set(md, 'add initialized data', 'alpha', [alpha]);
gf_model_set(md, 'add nonlinear term', mim, [sigmaeps '*(Grad_V.Grad_Test_V)']);
gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'V', elements_degree-1, RIGHT_BOUND);
gf_model_set(md, 'add initialized data', 'DdataV', [0.1]);
gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'V', elements_degree-1, LEFT_BOUND, 'DdataV');
|
熱問題¶
ここで，熱問題を記述するプログラムは次の通りです．
C++ |
md.add_initialized_scalar_data("kaeps", kappa*epsilon);
getfem::add_generic_elliptic_brick(md, mim, "theta", "kaeps");
md.add_initialized_scalar_data("D2", D*2);
md.add_initialized_scalar_data("D2airt", air_temp*D*2);
getfem::add_mass_brick(md, mim, "theta", "D2");
getfem::add_source_term_brick(md, mim, "theta", "D2airt");
getfem::add_nonlinear_term
(md, mim, "-"+sigmaeps+"*Norm_sqr(Grad_V)*Test_theta");
|
Python |
md.add_initialized_data('kaeps', [kappa*epsilon])
md.add_generic_elliptic_brick(mim, 'theta', 'kaeps')
md.add_initialized_data('D2', [D*2])
md.add_initialized_data('D2airt', [air_temp*D*2])
md.add_mass_brick(mim, 'theta', 'D2')
md.add_source_term_brick(mim, 'theta', 'D2airt')
md.add_nonlinear_term(mim, '-'+sigmaeps+'*Norm_sqr(Grad_V)*Test_theta')
|
Scilab |
gf_model_set(md, 'add initialized data', 'kaeps', [kappa*epsilon]);
gf_model_set(md, 'add generic elliptic brick', mim, 'theta', 'kaeps');
gf_model_set(md, 'add initialized data', 'D2', [D*2]);
gf_model_set(md, 'add initialized data', 'D2airt', [air_temp*D*2]);
gf_model_set(md, 'add mass brick', mim, 'theta', 'D2');
gf_model_set(md, 'add source term brick', mim, 'theta', 'D2airt');
gf_model_set(md, 'add nonlinear term', mim, '-'+sigmaeps+'*Norm_sqr(Grad_V)*Test_theta');
|
Matlab |
gf_model_set(md, 'add initialized data', 'kaeps', [kappa*epsilon]);
gf_model_set(md, 'add generic elliptic brick', mim, 'theta', 'kaeps');
gf_model_set(md, 'add initialized data', 'D2', [D*2]);
gf_model_set(md, 'add initialized data', 'D2airt', [air_temp*D*2]);
gf_model_set(md, 'add mass brick', mim, 'theta', 'D2');
gf_model_set(md, 'add source term brick', mim, 'theta', 'D2airt');
gf_model_set(md, 'add nonlinear term', mim, ['-' sigmaeps '*Norm_sqr(Grad_V)*Test_theta']);
|
モデルの求解¶
モデルを正しく定義したら，次のようにして簡単に解くことができます．
C++ |
gmm::iteration iter(1E-9, 1, 100);
getfem::standard_solve(md, iter);
|
Python |
md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
|
Scilab |
gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 'noisy');
|
Matlab |
gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 'noisy');
|
問題は大域的に非線形であるため，Newton法を用いて問題を反復的に解きます．数回のイテレーションが必要です (この場合は約4です)．
2ステップでのモデルの求解¶
別の解き方として，最初に熱と電位の問題を解くこともできます．今回のモデルでは，熱および電位は変形に依存しません．熱と電位の問題を解いてから変形の問題を解くには，次のように実行します．
C++ |
gmm::iteration iter(1E-9, 1, 100);
md.disable_variable("u");
getfem::standard_solve(md, iter);
md.enable_variable("u");
md.disable_variable("theta");
md.disable_variable("V");
iter.init();
getfem::standard_solve(md, iter);
|
Python |
md.disable_variable('u')
md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
md.enable_variable('u')
md.disable_variable('theta')
md.disable_variable('V')
md.solve('max_res', 1E-9, 'max_iter', 100, 'noisy')
|
Scilab |
gf_model_set(md, 'disable variable', 'u');
gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 'noisy');
gf_model_set(md, 'enable variable', 'u');
gf_model_set(md, 'disable variable', 'theta');
gf_model_set(md, 'disable variable', 'V');
gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 'noisy');
|
Matlab |
gf_model_set(md, 'disable variable', 'u');
gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 'noisy');
gf_model_set(md, 'enable variable', 'u');
gf_model_set(md, 'disable variable', 'theta');
gf_model_set(md, 'disable variable', 'V');
gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', 100, 'noisy');
|
解のエクスポート/可視化¶
以上で有限要素問題が解けました．図のように解をプロットすることができます．C++ および Python プログラムでは，外部のグラフィカルポストプロセッサを使用する必要があることに注意してください．また，汎用補間を使用して任意の数量を後処理できることにも注意してください(後述の ga_interpolation_Lagrange_fem を参照)．また，複雑なエクスポートやスライスの作成も可能です ( 解の出力と表示 をご覧ください)．
C++ |
plain_vector U(mfu.nb_dof()); gmm::copy(md.real_variable("u"), U);
plain_vector V(mft.nb_dof()); gmm::copy(md.real_variable("V"), V);
plain_vector THETA(mft.nb_dof()); gmm::copy(md.real_variable("theta"),THETA);
plain_vector VM(mfvm.nb_dof());
getfem::compute_isotropic_linearized_Von_Mises_or_Tresca
(md, "u", "clambdastar", "cmu", mfvm, VM, false);
plain_vector CO(mfvm.nb_dof() * 2);
getfem::ga_interpolation_Lagrange_fem(md, "-"+sigmaeps+"*Grad_V", mfvm, CO);
getfem::vtk_export exp("displacement_with_von_mises.vtk", false);
exp.exporting(mfu);
exp.write_point_data(mfu, U, "elastostatic displacement");
exp.write_point_data(mfvm, VM, "Von Mises stress");
cout << "\nYou can view solutions with for instance:\n\nmayavi2 "
"-d displacement_with_von_mises.vtk -f WarpVector -m Surface\n" << endl;
getfem::vtk_export exp2("temperature.vtk", false);
exp2.exporting(mft);
exp2.write_point_data(mft, THETA, "Temperature");
cout << "mayavi2 -d temperature.vtk -f WarpScalar -m Surface\n" << endl;
getfem::vtk_export exp3("electric_potential.vtk", false);
exp3.exporting(mft);
exp3.write_point_data(mft, V, "Electric potential");
cout << "mayavi2 -d electric_potential.vtk -f WarpScalar -m Surface\n"
<< endl;
}
|
Python |
U = md.variable('u')
V = md.variable('V')
THETA = md.variable('theta')
VM = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u', 'clambdastar', 'cmu', mfvm)
CO = np.reshape(md.interpolation('-'+sigmaeps+'*Grad_V', mfvm), (2, mfvm.nbdof()), 'F')
mfvm.export_to_vtk('displacement_with_von_mises.vtk', mfvm,
VM, 'Von Mises Stresses', mfu, U, 'Displacements')
print ('You can view solutions with for instance:')
print ('mayavi2 -d displacement_with_von_mises.vtk -f WarpVector -m Surface')
mft.export_to_vtk('temperature.vtk', mft, THETA, 'Temperature')
print ('mayavi2 -d temperature.vtk -f WarpScalar -m Surface')
mft.export_to_vtk('electric_potential.vtk', mft, V, 'Electric potential')
print ('mayavi2 -d electric_potential.vtk -f WarpScalar -m Surface')
|
Scilab |
U = gf_model_get(md, 'variable', 'u');
V = gf_model_get(md, 'variable', 'V');
THETA = gf_model_get(md, 'variable', 'theta');
VM = gf_model_get(md, 'compute_isotropic_linearized_Von_Mises_or_Tresca', 'u', 'clambdastar', 'cmu', mfvm);
CO = matrix(gf_model_get(md, 'interpolation', '-'+sigmaeps+'*Grad_V', mfvm), [2 gf_mesh_fem_get(mfvm, 'nbdof')]);
hh = scf(2);
hh.color_map = jetcolormap(255);
subplot(3,1,1);
gf_plot(mfvm, VM, 'mesh', 'off', 'deformed_mesh','off', 'deformation', U, 'deformation_mf', mfu, 'deformation_scale', 100, 'refine', 8);
colorbar(min(VM),max(VM));
title('Von Mises stress in N/cm^2 (on the deformed configuration, scale factor x100)');
subplot(3,1,2);
drawlater;
gf_plot(mft, V, 'mesh', 'off', 'deformed_mesh','off', 'deformation', U, 'deformation_mf', mfu, 'deformation_scale', 100, 'refine', 8);
colorbar(min(V),max(V));
gf_plot(mfvm, CO, 'quiver', 'on', 'quiver_density', 0.1, 'mesh', 'off', 'deformed_mesh','off', 'deformation_mf', mfu, ...
'deformation', U, 'deformation_scale', 100, 'refine', 8);
title('Electric potential in Volt (on the deformed configuration, scale factor x100)');
drawnow;
subplot(3,1,3);
gf_plot(mft, THETA, 'mesh', 'off', 'deformed_mesh','off', 'deformation', U, 'deformation_mf', mfu, 'deformation_scale', 100, 'refine', 8);
colorbar(min(THETA),max(THETA));
title('Temperature in °C (on the deformed configuration, scale factor x100)');
|
Matlab |
U = gf_model_get(md, 'variable', 'u');
V = gf_model_get(md, 'variable', 'V');
THETA = gf_model_get(md, 'variable', 'theta');
VM = gf_model_get(md, 'compute_isotropic_linearized_Von_Mises_or_Tresca', 'u', 'clambdastar', 'cmu', mfvm);
CO = reshape(gf_model_get(md, 'interpolation', ['-' sigmaeps '*Grad_V'], mfvm), [2 gf_mesh_fem_get(mfvm, 'nbdof')]);
figure(2);
subplot(3,1,1);
gf_plot(mfvm, VM, 'mesh', 'off', 'deformed_mesh','off', 'deformation', U, 'deformation_mf', mfu, 'deformation_scale', 100, 'refine', 8);
colorbar;
title('Von Mises stress in N/cm^2 (on the deformed configuration, scale factor x100)');
subplot(3,1,2);
hold on;
gf_plot(mft, V, 'mesh', 'off', 'deformed_mesh','off', 'deformation', U, 'deformation_mf', mfu, 'deformation_scale', 100, 'refine', 8);
colorbar;
gf_plot(mfvm, CO, 'quiver', 'on', 'quiver_density', 0.1, 'mesh', 'off', 'deformed_mesh','off', 'deformation', U, 'deformation_mf', ...
mfu, 'deformation_scale', 100, 'refine', 8);
colorbar;
title('Electric potential in Volt (on the deformed configuration, scale factor x100)');
hold off;
subplot(3,1,3);
gf_plot(mft, THETA, 'mesh', 'off', 'deformed_mesh','off', 'deformation', U, 'deformation_mf', mfu, 'deformation_scale', 100, 'refine', 8);
colorbar;
title('Temperature in ?C (on the deformed configuration, scale factor x100)');
|
解のプロット．¶