高性能要素定式化¶

3 次元ソリッド要素の定式化で示した変位ベースの標準定式化は，非圧縮に近い材料や曲げ支配の薄肉構造に適用するとロッキング（体積ロッキング・せん断ロッキング）と呼ばれる不当な過剛性を示す．これを回避するため FrontISTR は，B 行列または変形勾配の体積部のみを差し替える B-bar 法・F-bar 法，内部自由度を導入する非適合要素，圧力を独立な未知場とする u-p 混合要素，および板・梁構造に特化した MITC シェル要素・梁要素を備える．

本章ではこれら高性能要素・構造要素の定式化を要素単位で整理する．

B-bar 法¶

8 節点線形六面体要素を準非圧縮材料に用いると，1 要素内のひずみが体積一定の制約と矛盾し，体積ロッキングと呼ばれる過剛性が現れる．B-bar 法は B 行列のうち体積膨張に寄与する成分を要素中心で評価した値で置き換え，制約過多を緩和する [Hughes1980]．

要素中心 \(\boldsymbol{r} = \boldsymbol{0}\) で形状関数の空間微分から計算した B 行列を \(\bar{\boldsymbol{B}}\)，積分点 \(\boldsymbol{r}\) で計算した通常の B 行列を \(\boldsymbol{B}(\boldsymbol{r})\) とする．節点 \(\alpha\) の自由度 \(i\) に対する変位－ひずみ関係において，体積ひずみ成分 \((\varepsilon_{11}, \varepsilon_{22}, \varepsilon_{33})\) には

\[ \Delta B_{i\alpha} = \tfrac{1}{3}\bigl(\bar{B}_{i\alpha}(\boldsymbol{0}) - B_{i\alpha}(\boldsymbol{r})\bigr) \]

を加え，せん断成分 \((\varepsilon_{12}, \varepsilon_{23}, \varepsilon_{31})\) には通常の \(\boldsymbol{B}\) を用いる．こうして得た B-bar 行列で要素剛性と内力ベクトルを組み立てる．

FrontISTR ではこの定式化を 8 節点線形六面体（要素 ID 361，要素番号体系）専用に提供し，微小変形・Total Lagrange・Updated Lagrange のいずれにも適用できる．

F-bar 法¶

有限変形下では体積変化が変形勾配 \(\boldsymbol{F}\) を通じて非線形に効くため，B-bar 法と同等の体積ロッキング対策を変形勾配の段階で行うのが F-bar 法 [deSouzaNeto1996] である．

要素中心 \(\boldsymbol{r} = \boldsymbol{0}\) で評価した変形勾配の体積比を \(J_0 = \det \boldsymbol{F}(\boldsymbol{0})\)，積分点での変形勾配の体積比を \(J = \det \boldsymbol{F}(\boldsymbol{r})\) とし，積分点の変形勾配を

\[ \bar{\boldsymbol{F}} = \left(\frac{J_0}{J}\right)^{1/3} \boldsymbol{F} \]

で差し替える．これにより \(\det \bar{\boldsymbol{F}} = J_0\) となり，要素全域で体積比が要素中心値に揃う．応力評価とひずみ－変位行列の構築には差し替え後の \(\bar{\boldsymbol{F}}\) を用い，要素剛性ではこの差し替えに伴う追加項を含めて接線行列を組み立てる．

FrontISTR では F-bar 法を 8 節点線形六面体専用に実装し，微小変形・Total Lagrange / Updated Lagrange のいずれの非線形幾何にも適用できる．

非適合要素¶

8 節点線形六面体要素は曲げモードに対して必要なひずみ自由度を持たず，曲げ支配の問題で曲げロッキングを示す．非適合要素 [Taylor1976] はこの欠陥を補うため，要素内に追加の変位モードを導入する．

要素節点変位 \(\boldsymbol{u}^e\) に加えて，要素内部にのみ存在する非適合モード自由度 \(\boldsymbol{\alpha} \in \mathbb{R}^{9}\) を 1 要素あたり 3 方向 × 3 モード分導入し，変位場を

\[ \boldsymbol{u}(\boldsymbol{r}) = \sum_{\alpha=1}^{8} N_\alpha^e(\boldsymbol{r})\, \boldsymbol{u}^e_\alpha + \sum_{k=1}^{3} M_k(\boldsymbol{r})\, \boldsymbol{\alpha}_k \]

で近似する．非適合形状関数は自然座標 \(\boldsymbol{r} = (\xi, \eta, \zeta)\) に対し \(M_1 = 1 - \xi^2\)，\(M_2 = 1 - \eta^2\)，\(M_3 = 1 - \zeta^2\) ととる．これらは要素境界で連続性を保証しないが，曲げモードのひずみを再現する空間を要素内に追加する．

要素剛性は外部自由度－内部自由度のブロック分解

\[ \begin{bmatrix} \boldsymbol{K}_{dd} & \boldsymbol{K}_{d\alpha} \\ \boldsymbol{K}_{\alpha d} & \boldsymbol{K}_{\alpha\alpha} \end{bmatrix} \begin{bmatrix} d\boldsymbol{u}^e \\ d\boldsymbol{\alpha} \end{bmatrix} = \begin{bmatrix} \boldsymbol{F}^e_{\text{ext}} \\ \boldsymbol{0} \end{bmatrix} \]

の形に組み立てた後，内部自由度を \(d\boldsymbol{\alpha} = -\boldsymbol{K}_{\alpha\alpha}^{-1}\boldsymbol{K}_{\alpha d}\,d\boldsymbol{u}^e\) で消去する静的縮約を行い，外部自由度のみの要素剛性

\[ \boldsymbol{K}^e = \boldsymbol{K}_{dd} - \boldsymbol{K}_{d\alpha}\,\boldsymbol{K}_{\alpha\alpha}^{-1}\,\boldsymbol{K}_{\alpha d} \]

を全体組立に渡す．

FrontISTR では非適合要素を 8 節点線形六面体（C3D8IC）専用に実装し，微小変形・Total Lagrange・Updated Lagrange のいずれにも適用できる．

U-P 混合要素¶

B-bar 法・F-bar 法が変位ベースの枠組みの中で体積成分を補正するのに対し，u-p 混合（U-P）要素は圧力 \(\lambda\) を変位とは独立な未知場として導入する混合定式化 [Bathe1996] である．準非圧縮材料（ポアソン比が 0.5 に極めて近いゴム状材料や塑性変形後の金属など）では，変位場だけで体積一定制約を満たそうとすると体積ロッキングが生じるが，圧力を独立変数とすることでこの制約を緩和する．

応力を偏差成分と圧力に分離し，

\[ \boldsymbol{\sigma} = \boldsymbol{\sigma}_{\mathrm{dev}} + \lambda\,\boldsymbol{I}, \qquad \boldsymbol{\sigma}_{\mathrm{dev}} = \mathbf{D}_{\mathrm{dev}}\,\boldsymbol{\varepsilon} \]

と書く．ここで \(\mathbf{D}_{\mathrm{dev}}\) は弾性行列からバルク係数 \(K\) に比例する体積部を除いた偏差弾性行列である．圧力 \(\lambda\) と体積ひずみ \(g = \mathrm{tr}\,\boldsymbol{\varepsilon}\) は，圧縮率 \(\alpha^{-1} = 1/K\) を介した拘束条件

\[ g - \alpha^{-1}\lambda = 0 \]

で結ばれる．変位 \(\boldsymbol{u}\) と圧力 \(\lambda\) を未知数とする離散化から，要素の連立式

\[ \begin{bmatrix} \mathbf{K}_{uu} & \mathbf{K}_{up} \\ \mathbf{K}_{up}^{T} & \mathbf{K}_{pp} \end{bmatrix} \begin{bmatrix} d\boldsymbol{u} \\ d\lambda \end{bmatrix} = \begin{bmatrix} \boldsymbol{f}_{u} \\ \boldsymbol{f}_{p} \end{bmatrix} \]

が得られる．ここで \(\mathbf{K}_{uu}\) は偏差弾性成分と（有限変形では）幾何剛性，\(\mathbf{K}_{up}\) は体積ひずみと圧力を結合する行列，\(\mathbf{K}_{pp} = -\int \alpha^{-1}\,\boldsymbol{N}_p \boldsymbol{N}_p^{T}\,dV\) は圧力の安定化項である（\(\boldsymbol{N}_p\) は圧力の形状関数）．圧力自由度は要素内に閉じているため，

\[ \mathbf{K}_{\mathrm{eff}} = \mathbf{K}_{uu} - \mathbf{K}_{up}\,\mathbf{K}_{pp}^{-1}\,\mathbf{K}_{up}^{T} \]

の形に静的縮約して外部（変位）自由度のみの有効剛性を全体組立に渡す．

FrontISTR では U-P 要素を 8 節点線形六面体専用に実装し，圧力自由度は 1 要素あたり 1 個（要素内で一定）とする．微小変形・Total Lagrange・Updated Lagrange のいずれにも適用でき，Updated Lagrange では偏差応力を客観応力速度（Jaumann／Hughes-Winget 型）で更新したうえで，圧力 \(\lambda\,\boldsymbol{I}\) を静的縮約で求めた値から課す．

シェル要素¶

薄肉板殻構造には Reissner-Mindlin 板殻理論に基づくシェル要素を用いる．低次の変位ベース板殻要素は厚さが薄くなるにつれ横せん断ひずみが体積ロッキングと類似の機構（せん断ロッキング）で過大評価され，曲げモードに対する剛性が発散的に増大する．MITC（Mixed Interpolation of Tensorial Components）法 [Dvorkin1984] [Bathe1986] は，せん断ひずみ成分のみを要素内のあらかじめ定めたタイイング点で再標本化し，標本値を要素内に内挿し直すことでこの問題を回避する．

MITC シェル要素の節点は曲面上の中央面節点であり，各節点の自由度は並進 3 成分と中央面法線まわりの回転 3 成分を合わせた 6 自由度を持つ．要素剛性は中央面の自然座標と厚さ方向の合計 3 次元のガウス積分で評価し，板厚 \(h\) は構成則計算時に要素属性として与える．

FrontISTR は中央面 1 層で扱う MITC3（要素 ID 731）・MITC4（741）・MITC9（743）と，板厚を貫く方向に節点を 2 層配置した層状シェル要素 MITC3-shell361（761，3\(\times\)2 節点・節点 3 自由度）・MITC4-shell361（781，4\(\times\)2 節点・節点 3 自由度）を提供する．層状シェル要素では節点自由度は並進 3 成分のみとし，回転自由度に相当する曲げモードを 2 層配置で表現する．

梁要素¶

桁・骨組構造などの線材は梁要素で離散化する．FrontISTR はせん断変形を考慮した Timoshenko 梁定式化を用い，曲げとせん断の両者を変位および回転自由度の関数として記述する．

梁要素の節点自由度は並進 3 成分と軸まわり・横軸まわりの回転 3 成分を合わせた 6 自由度を持ち，要素剛性は梁の軸方向に沿った 1 次元の数値積分で評価する．断面積 \(A\) と曲げ・捩り方向の断面 2 次モーメント \(I\) は梁の断面定数として要素属性に与え，材料の縦弾性係数 \(E\)・せん断弾性係数 \(G\) とともに伸び・曲げ・捩り・せん断の各剛性係数を構成する．

FrontISTR は 2 節点直線梁要素（要素 ID 611）と，3 節点で表現する 4 節点四面体ソリッド－梁ハイブリッド要素（641，自由度混在用）を提供する．