シェル要素の不具合修正（積層質量・直交異方性副対称性・繊維角度単位・6自由度固有値出力）

inagakirequested to merge
fix/amakawa-shell-bugs into master

Closes #763 (closed), #764 (closed), #765 (closed), #766 (closed), #773 (closed)

Hiroaki Amakawa さん報告のシェル要素関連の不具合5件をまとめて修正します。

修正内容

  • #763 (closed) 積層シェル質量行列が層数倍 — static_LIB_shell.f90 の mass_Shell 層積分ループに層厚ヤコビアン shell_var(n_layer)%weight が欠落しており、質量が n_layers 倍（固有振動数が 1/√n_layers 倍）になっていた。剛性ループと同じ重みを畳み込み。単層（weight=1.0）は不変。
  • #764 (closed) / #765 (closed) flag=1 で c_hat の副対称性欠落 — ElasticLinear.f90 の D_hat→c_hat 変換ループが逆順インデックスを持たず、副対称成分が 0 のまま残っていた。不規則格子で剛性過大（#764 (closed)）、繊維角 ≠ 0 で偽の零エネルギーモード（#765 (closed)）。ループ内で 4 つの副対称置換へ同時代入し、対角せん断＋結合（オフ対角）成分の双方を補完。
  • #766 (closed) flag=1 の繊維角 angle に度→ラジアン変換が欠落 — fstr_get_prop.f90 で度入力をそのままラジアンとして扱っていた（梁は変換済み）。読み込み時に angle/180×π を適用。angle=0 は不変。
  • #773 (closed) 6自由度シェルの固有値解析で変形可視化不可 — fstr_EIG_output.f90 が固有ベクトルを 6 成分 DISPLACEMENT として出力し、VTK/ParaView が対称テンソルと解釈して Warp By Vector が使えなかった。make_result.f90（ndof==6）に倣い DISPLACEMENT（並進 3）と ROTATION（回転 3）に分離。NDOF==3（ソリッド）は従来どおり。

検証

各 issue 添付モデルを修正版バイナリで実行し、issue 記載の「修正後（after）」期待値と一致を確認しました。

  • #763 (closed): 3 層積層が単層基準と質量・全モード完全一致（修正前は ×3、振動数 1/√3）
  • #764 (closed): 等方相当 flag=1 が flag=0 基準と全モード ≤0.3% 一致（修正前 最大 +2.4%）
  • #765 (closed): −45° の偽零モードが消失、±45° とも after 値を厳密再現
  • #766 (closed): 0° と 180° が完全一致（0°≡180° 不変性の回復、修正前は 180°=1.834 Hz と乖離）
  • #773 (closed): 結果ファイルの DISPLACEMENT が 6→3 成分 ＋ ROTATION 3 成分（res ヘッダ 1 0/6→2 0/3 3）

回帰テスト

tests/analysis/eigen/shell_laminate/ に小規模 5 モデルを追加（参照結果は修正版バイナリで生成）。

テスト 要素 ガード対象
K741LAMINATE3 MITC4 #763 (closed) 質量行列層数倍
K731ORTHOISO MITC3 #764 (closed) c_hat 対角副対称
K731FIBERM45 / K731FIBERP45 MITC3 #765 (closed) c_hat 結合副対称（零モード）
K741FIBER180 MITC4 #766 (closed) 角度 度→ラジアン
  • ソルバは既存 analysis/eigen の慣習に合わせ CG（PRECOND=3）に統一。MUMPS 非依存で常時実行でき、当該モデルでは CG と MUMPS で固有値が一致します。
  • #765 (closed) の繊維角は修正後の単位である度（±45.0）で記述。
  • compare_res.pl は eigen では EIGENVALUE のみ比較するため、#763 (closed)–#766 (closed) の数値退行を捕捉します（#773 (closed) の出力構造そのものは本ハーネスの対象外）。

Merge request reports