等方性線形弾性ブリック要素

このブリック要素は項

\[-div(\sigma) = \ldots\]

\[\begin{split}\sigma &= \lambda\mbox{tr}(\varepsilon(u))I + 2\mu\varepsilon(u) \\ \varepsilon(u) &= (\nabla u + \nabla u^T)/2\end{split}\]

を表します. \(\varepsilon(u)\) は微小ひずみテンソルです. \(\sigma\) は応力テンソルです. \(\lambda\)\(\mu\) は Lamé 係数です.これは線形等方性弾性の系を表します.Stokes問題を構築するために線形非圧縮要素と一緒に \(\lambda = 0\) で使うこともできます.

Lamé 係数とYoung係数 \(E\) とPoisson比 \(\nu\)

\[\lambda = \dfrac{E\nu}{(1+\nu)(1-2\nu)}, ~~~ \mu = \dfrac{E}{2(1+\nu)},\]

平面応力近似(2次元モデル)を除き,

\[\lambda^* = \dfrac{E\nu}{(1-\nu^2)}, ~~~ \mu = \dfrac{E}{2(1+\nu)},\]

このブリック要素をモデルに追加し Lamé 係数をパラメータとするモデルの関数は次の通りです.

ind_brick = getfem::add_isotropic_linearized_elasticity_brick
            (md, mim, varname, data_lambda, data_mu,
             region = size_type(-1));

ここで, dataname_lambdadataname_mu は, Lamé 係数を表すモデルのデータです.

このブリック要素をモデルに追加し,Young率とPoisson比でパラメータ化した関数は次の通りです.

ind_brick = getfem::add_isotropic_linearized_elasticity_pstrain_brick
            (md, mim, varname, data_E, data_nu, region = size_type(-1));

この要素は,2次元メッシュ(および3次元メッシュ上の標準モデル)に適用されるときの平面歪み近似を表します. 次の関数により2次元メッシュの平面応力近似を得ることができます.

ind_brick = getfem::add_isotropic_linearized_elasticity_pstress_brick
            (md, mim, varname, data_E, data_nu, region = size_type(-1));

3次元メッシュの場合,前の2つの要素で同じ結果が得られます.

関数は次の通りです.

getfem::compute_isotropic_linearized_Von_Mises_or_Tresca
  (md, varname, dataname_lambda, dataname_mu, mf_vm, VM, tresca_flag = false);

varname に格納されている変位フィールド上のVon Mises基準( tresca_flag がtrueに設定されている場合はTresca)を計算します.応力は mesh_fem mf_vm を呼び出して VM というベクトルに格納します. 2次元平面応力近似には有効ではなく, Lamé 係数でパラメータ化されています.関数は次の通りです.

getfem::compute_isotropic_linearized_Von_Mises
  (md, varname, data_E, data_nu, mf_vm, VM);

getfem::compute_isotropic_linearized_Von_Mises
  (md, varname, data_E, data_nu, mf_vm, VM);

Young率とPoisson比でパラメータ化されたVon Mises応力を計算し,2番目の方程式が2次元メッシュに適用されたときに2次元平面応力近似は有効です(2つの関数は3次元問題に対して同じ結果を与えます).

プログラム tests/elastostatic.cc は線形化された等方性弾性要素の使用モデルを得ることができます.

線形非圧縮性(またはほぼ非圧縮性)ブリック要素

このブリック要素は次のような種類の問題で線形非圧縮性条件(またはほぼ非圧縮性の条件)を追加します.

\[\mbox{div}(u) = 0,\quad (\mbox{ or } \mbox{div}(u) = \varepsilon p)\]

この制約は混合式で導入された圧力を表すLagrange乗数によって強制されます.

この非圧縮性条件を追加する関数は次の通りです.

ind_brick = getfem::add_linear_incompressibility
            (md, mim, varname, multname_pressure, region = size_type(-1),
             dataexpr_penal_coeff = std::string());

ここで varname は非圧縮性条件が規定されている変数であり, multname_pressure は乗数(圧力)を表すスカラー有限要素法で記述される変数であり, dataexpr_penal_coeff はほぼ非圧縮状態のオプションのペナルティ係数です.

ほぼ非圧縮性の均質線形弾性では, \(\varepsilon = 1 / \lambda\) のようになります.ここで, \(\lambda\) はLamé係数の1つであり,\(\varepsilon\) はペナルティ係数です.

例えば,次のプログラムは,境界項0上のソース項と均質なDirichlet条件を持つStokes問題を定義しています. mf_umf_datamf_p は同じメッシュ上の有効な有限要素法でなければなりません. mim は同じメッシュ上で有効な積分法でなければなりません.

typedef std::vector<getfem::scalar_type> plain_vector;
size_type N = mf_u.linked_mesh().dim();

getfem::model Stokes_model;

laplacian_model.add_fem_variable("u", mf_u);

getfem::scalar_type mu = 1.0;
Stokes_model.add_initialized_data("lambda", plain_vector(1, 0.0));
Stokes_model.add_initialized_data("mu", plain_vector(1, mu));

getfem::add_isotropic_linearized_elasticity_brick(Stokes_model, mim,
                                                  "u", "lambda", "mu");

laplacian_model.add_fem_variable("p", mf_p);
getfem::add_linear_incompressibility(Stokes_model, mim, "u", "p");

plain_vector F(mf_data.nb_dof()*N);
for (int i = 0; i < mf_data.nb_dof()*N; ++i) F(i) = ...;
Stokes_model.add_initialized_fem_data("VolumicData", mf_data, F);
getfem::add_source_term_brick(Stokes_model, mim, "u", "VolumicData");

getfem::add_Dirichlet_condition_with_multipliers(Stokes_model, mim,
                                                 "u", mf_u, 1);

gmm::iteration iter(residual, 1, 40000);
getfem::standard_solve(Stokes_model, iter);

plain_vector U(mf_u.nb_dof());
gmm::copy(Stokes_model.real_variable("u"), U);

ほぼ非圧縮性条件の例は,プログラム tests/elastostatic.cc にあります.