# Bug #2-1 — flag=1（直交異方性）で c_hat の副対称性が欠落 **対象バージョン**: FrontISTR 5.9 — commit `6d3a1684` (2026-05-19) **修正コミット**: `fdf8efbe`（ブランチ `debug/bug2-c-hat-symmetry`） --- ## 1. 現象 直交異方性積層材（flag=1, `!ITEM=1, SUBITEM = 1+8×n_layers`）を MITC3 三角形要素 (731) に 使用すると、等方性相当の物性を入力しても等方性材（flag=0）より固有振動数が高くなる （剛性過大）。MITC4 四辺形要素や規則格子では差がほとんど出ず、**不規則（デローニ）格子の MITC3 で顕著**に現れる。 本フォルダは、同一の不規則格子・同一物性で flag=0 と flag=1 を比較し、 両者が一致すべきこと（=flag=1 のずれがバグであること）を確認する。 --- ## 2. 原因 `fistr1/src/lib/physics/ElasticLinear.f90`、`LinearElastic_Shell` サブルーチンの flag=1 ブランチ。4 階構成テンソル `c_hat(i,j,k,l)` は副対称性 ``` c(i,j,k,l) = c(j,i,k,l) = c(i,j,l,k) = c(j,i,l,k) ``` を満たす必要がある。flag=0 ブランチはせん断成分を 4 つすべて明示設定しているが、 flag=1 ブランチは D_hat(5×5) → c_hat への変換に用いるインデックス配列に逆順成分 `(2,1) (3,2) (1,3)` を含まないため、対称化ループがスキップされ、以下 9 成分が 0 のまま残る。 | 欠落成分 | 正しい値 | |------|------| | `c_hat(2,1,1,2)`, `c_hat(1,2,2,1)`, `c_hat(2,1,2,1)` | G12 | | `c_hat(3,2,2,3)`, `c_hat(2,3,3,2)`, `c_hat(3,2,3,2)` | G23 | | `c_hat(1,3,3,1)`, `c_hat(3,1,1,3)`, `c_hat(1,3,1,3)` | G31 | ### MITC3 デローニ格子にのみ顕著な理由 剛性積分は局所デカルト基底 `e_hat` と反変基底の内積 `e_hat_dot_cg` を介して座標変換される。 規則格子（軸整列）では `e_hat_dot_cg ≈ I` のため欠落成分が積に現れないが、 デローニ等の一般配置では `e_hat_dot_cg ≠ I` となり、欠落成分が剛性に寄与して過大評価になる。 --- ## 3. 修正 `ElasticLinear.f90` の変換ループ直後（line 419〜）に副対称 9 成分を明示補完する。 ```fortran end do end do ! minor symmetry 補完: flag=1 ブランチはインデックス配列に逆順成分がないため欠落 c_hat(2,1,1,2) = c_hat(1,2,1,2) c_hat(1,2,2,1) = c_hat(1,2,1,2) c_hat(2,1,2,1) = c_hat(1,2,1,2) c_hat(3,2,2,3) = c_hat(2,3,2,3) c_hat(2,3,3,2) = c_hat(2,3,2,3) c_hat(3,2,3,2) = c_hat(2,3,2,3) c_hat(1,3,3,1) = c_hat(3,1,3,1) c_hat(3,1,1,3) = c_hat(3,1,3,1) c_hat(1,3,1,3) = c_hat(3,1,3,1) ``` > **注記（中間的な修正）**: この 9 行補完は最終コードには残らない **中間的な変更** である。 後続の Bug #2-2（`ce5ba6de`）が、対角限定のこの補完を **D_hat ループ内での全副対称成分の補完に 置き換え**（9 行は削除）、θ≠0 で生じる結合成分まで含めて統合した。出荷される最終形はループ方式 （`ce5ba6de`）であり、本 `fdf8efbe` はその前段にあたる。 --- ## 4. テストモデル アルミ片持ち梁を MITC3 デローニ格子でモデル化し、5 モードの固有値解析を行う。 | 項目 | 値 | |----|---| | 形状 | 500 × 100 mm | | 要素 | MITC3 (type 731)、**デローニ不規則格子** | | 規模 | 164 節点・267 要素 | | 板厚 | 10 mm（1 層、積分点 2） | | 境界条件 | x=0 端 全 DOF 固定 (`Constraint_1`, 1–6) | | 解析 | `!SOLUTION, TYPE=EIGEN`、5 モード | | ヤング率 E | 71000 MPa | | ポアソン比 ν | 0.33 | | 密度 ρ | 2.770E-9 t/mm³ | 1 層モデルのため Bug #1（質量行列の層数倍）の影響はなく、密度補正も不要。 --- ## 5. ケース構成 同一格子・同一物性で flag=0 と flag=1 を比較する。flag=1 はせん断弾性率を 等方相当 `G = E / (2(1+ν)) = 71000 / 2.66 ≈ 26692 MPa` とし、繊維角 angle=0 とすることで 物理的に flag=0 と等価にしてある。理論上は flag=0 と flag=1 が一致すべき。 | ケース | 要素 | 材料 | flag | 材料定義（`!ITEM=1`） | |-----|----|----|:----:|-----------------| | `case_C_731_flag0` | 731 / デローニ | `AL_ISO`（等方） | 0 | `SUBITEM=4`: `0, 71000.0, 0.33, 10.0` | | `case_D_731_flag1` | 731 / デローニ | `AL_ORTH`（直交異方・等方相当） | 1 | `SUBITEM=9`: `1, 71000.0, 0.33, 71000.0, 26692.0, 26692.0, 26692.0, 0.0, 10.0` | 各シナリオを before（未修正）/ after（修正版）の両方で実行。フラット構成（01,02=before / 03,04=after）: `01_case_C_731_flag0_before` / `02_case_D_731_flag1_before` / `03_case_C_731_flag0_after` / `04_case_D_731_flag1_after`。 - `case_C`（flag=0）: 修正の影響を受けず before = after。比較の基準。 - `case_D`（flag=1）: 修正前は基準よりずれる（剛性過大）。修正後は基準にほぼ一致する。 各ケースディレクトリは入力一式（`Job_1.msh` / `Job_1.cnt` / `hecmw_ctrl.dat`）と 実行ログ（`run.log` ほか）・結果（`Job_1.res.0.*`）を含む。 --- ## 6. 検証結果 固有振動数 \[Hz\]（5 モード）。`before` = 未修正バイナリ（git_hash `6d3a1684`）、 `after` = 修正版（git_hash `fdf8efbe`）。TOTAL MASS は全ケース `1.3850E-03` で一致 （本バグは剛性側のため質量は不変）。 | 実行 | flag | mode1 | mode2 | mode3 | mode4 | mode5 | |----|:----:|------:|------:|------:|------:|------:| | `01_case_C_731_flag0_before` | 0（基準） | 33.20 | 207.2 | 310.4 | 328.8 | 579.8 | | `02_case_D_731_flag1_before` | 1 | 33.78 | 211.9 | 308.1 | 327.5 | 593.8 | | `03_case_C_731_flag0_after` | 0（基準） | 33.20 | 207.2 | 310.4 | 328.8 | 579.8 | | `04_case_D_731_flag1_after` | 1 | 33.21 | 207.4 | 311.2 | 328.8 | 580.4 | **flag=1 の flag=0（基準）からのずれ** | モード | before（flag1 vs flag0） | after（flag1 vs flag0） | |-----|-----------------------:|----------------------:| | mode1 | +1.7% | +0.03% | | mode2 | +2.3% | +0.1% | | mode3 | −0.7% | +0.3% | | mode4 | −0.4% | ±0.0% | | mode5 | +2.4% | +0.1% | **確認ポイント** - flag=0（`03_C`）は修正の影響を受けず before = after（基準として不変）。 - flag=1（`04_D`）は修正前に基準から最大 **+2.4%**（mode5）ずれていたが、 修正後は全モードで **+0.3% 以下**に収束し、flag=0 とほぼ一致 :white_check_mark:。 - 同一物性（等方相当）の flag=0 / flag=1 が一致したことで、c_hat 副対称性の補完が 正しく効いていることを確認。 [03_case_C_731_flag0_after.7z](/uploads/70032cff3ac44ac2659a529e54bdb907/03_case_C_731_flag0_after.7z) [02_case_D_731_flag1_before.7z](/uploads/f20c5de04b0dc8946b1a9adf961f1437/02_case_D_731_flag1_before.7z) [01_case_C_731_flag0_before.7z](/uploads/681684d123d2a693acd49789f1a40f3c/01_case_C_731_flag0_before.7z) [0002-fix-c_hat-minor-symmetry-Bug-2.patch](/uploads/c8a89165563c49861df4b9a36d01030b/0002-fix-c_hat-minor-symmetry-Bug-2.patch)