flag=1（直交異方性）で c_hat の副対称性が欠落
Bug #2 (closed)-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 成分を明示補完する。
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 (closed)-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 (closed)（質量行列の層数倍）の影響はなく、密度補正も不要。
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 とほぼ一致✅ 。 - 同一物性（等方相当）の flag=0 / flag=1 が一致したことで、c_hat 副対称性の補完が 正しく効いていることを確認。
01_case_C_731_flag0_before.7z 0002-fix-c_hat-minor-symmetry-Bug-2.patch