積層シェル要素の質量行列が層数倍になる
Bug #1 (closed) — 積層シェル要素の質量行列が層数倍になる
対象バージョン: FrontISTR 5.9 — commit 6d3a1684 (2026-05-19) 修正コミット: df42b03c（ブランチ debug/bug1-mass-matrix-fix）
1. 現象
積層シェル材料（!ITEM=1, SUBITEM ≥ 7）を用いた解析で、質量行列が層数 n_layers 倍に 過大評価される。固有値解析では固有振動数が正しい値の 1/√n_layers 倍 に過小評価される。
等価な単層モデルと比較することで確認できる。本フォルダの 3 ケースがその比較を構成している。
2. 原因
fistr1/src/lib/static_LIB_shell.f90 の mass_Shell サブルーチン、質量行列の層積分ループ。
積層シェルでは各層のガウス積分に層厚ヤコビアン shell_var(n_layer)%weight （正規化層厚、Σ = 1.0）を乗じる必要がある。剛性行列ループ（同ファイル line 1064）では 正しく乗じられているが、質量行列ループではこの因子が欠落していた。
各層のガウス積分は層区間 [-1+Σ_{m<n}2w_m, -1+Σ_{m≤n}2w_m] に対する積分で、 変数変換ヤコビアンが層厚分率 w_n になる。この因子が無いと各層が全厚 [-1,1] 相当の 寄与を返し、n_layers 層を足し合わせて質量が n_layers 倍になる。
3. 修正
static_LIB_shell.f90 line 3548 の 1 行修正。w_w_w_det に層厚ヤコビアンを畳み込むことで、 直後の mass 行と totalmass 行の両方が同時に正しくなる。
!--------------------------------------------
- w_w_w_det = w_w_lx*w_ly*det
+ w_w_w_det = w_w_lx*w_ly*det*gausses(1)%pMaterial%shell_var(n_layer)%weight
mass(1:nsize,1:nsize) = mass(1:nsize,1:nsize)+ matmul( transpose(N), N )*w_w_w_det*rho
totalmass = totalmass + w_w_w_det*rho
!--------------------------------------------n_layers = 1 のとき weight は正規化後 1.0 になるため、単層モデルの結果は不変。
4. テストモデル
アルミ片持ち梁。x=0 端を全 DOF 拘束し、5 モードの固有値解析（!SOLUTION, TYPE=EIGEN）。
| 項目 | 値 |
|---|---|
| 形状 | 500 × 100 mm |
| 総厚 | 30 mm |
| 要素 | MITC4 (type 741)、10 × 2 = 20 要素 |
| 節点 | 33 節点（規則格子） |
| 境界条件 | x=0 端 全 DOF 固定 (Constraint_1, 1–6) |
| ヤング率 E | 71000 MPa |
| ポアソン比 ν | 0.33 |
| 密度 ρ | 2.770E-9 t/mm³ |
5. ケース構成
| ディレクトリ | モデル | 材料定義 | 層構成 |
|---|---|---|---|
01_reference_1layer/ |
等価単層（基準） | 等方 SUBITEM=2 |
30 mm × 1 |
02_before_fix/ |
3 層積層・修正前バイナリ | 等厚積層 SUBITEM=10, flag=0 |
10 mm × 3 |
03_after_fix/ |
3 層積層・修正後バイナリ | 等厚積層 SUBITEM=10, flag=0 |
10 mm × 3 |
3 層は同一物性（E=71000, ν=0.33）を 3 枚重ねたもので、物理的には単層 30 mm と等価。 したがって正しい実装なら 01 と 03 は一致し、02（修正前）だけがずれる。
各ケースは入力一式（Job_1.msh / Job_1.cnt / hecmw_ctrl.dat）と 実行ログ（run.log ほか）・結果（Job_1.res.0.*）を含む。
6. 検証結果
TOTAL MASS と低次 3 モードの固有振動数：
| ケース | TOTAL MASS | mode1 [Hz] | mode2 [Hz] | mode3 [Hz] | 判定 |
|---|---|---|---|---|---|
01_reference_1layer（基準） |
4.1550E-03 | 98.81 | 336.3 | 608.6 | — |
02_before_fix（3層・修正前） |
1.2465E-02 | 57.05 | 194.2 | 351.4 | バグ |
03_after_fix（3層・修正後） |
4.1550E-03 | 98.81 | 336.3 | 608.6 | 修正 |
確認ポイント
- 質量: 修正前は
1.2465E-02 = 4.1550E-03 × 3（ちょうど 3 倍）。修正後は基準と完全一致。 - 振動数: 修正前は全モードで基準の
1/√3倍に過小評価。- mode1: 98.81 / 57.05 = 1.732 = √3
- mode2: 336.3 / 194.2 = 1.732 = √3
- mode3: 608.6 / 351.4 = 1.732 = √3
- 修正後は
03が基準01と全モード一致✅