線形静解析の熱荷重がsubstepごとに過大になる不具合の修正

Hiroaki Amakawarequested to merge
Hiroaki_Amakawa/FrontISTR:bugfix/linear-thermal-substep into master

線形静解析の熱荷重がsubstepごとに過大になる不具合の修正

Issue #908 (closed) に対するMRです。

線形静解析（!SOLUTION, TYPE=STATIC）で弾性材料に熱荷重をかけ、!VERSION 5 で !STEP, SUBSTEPS を2以上にすると、substepが進むにつれて変位・応力が累積的に過大になる。10 substepでは最終変位が正解の約4.6倍になる。計算はエラーなく正常終了するため、利用者は誤りに気づきにくい。

不具合の内容

線形解析（INFINITESIMAL）の応力更新は増分形 σ = σ_bak + D (dstrain − strain_bak) である。弾性材料では、dstrain が初期温度基準の全熱ひずみを差し引いた機械ひずみであるのに対し、strain_bak は熱ひずみを含む全ひずみである。そのためsubstepごとに、前substepまでの熱ひずみに相当する応力が余分に差し引かれる。

期待される動作と現状

  • 期待される動作：線形解析では荷重係数が時間に比例するので、substep n の変位は1step解の n/10 倍になり、最終substepの結果は SUBSTEPS=1 の解と一致する。
  • 現状：最大変位の比（対 1step解）が次のとおり累積的に大きくなる。
substep 1 2 3 4 5 6 7 8 9 10
修正前 0.100 0.200 0.400 0.700 1.100 1.600 2.200 2.900 3.701 4.601
期待値 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

不具合の再現方法

  • 使用バイナリ：origin/master（ 007551a0）をクリーンビルドしたもの（bin_base/fistr1）。ビルド環境は env_build.sh、構成ログは out_cmake.log、ビルドログは out_build_base.log。
  • 入力：ユーザー提供の /data/a/a878/frontistr/calc/test_thermalstress/ の2ケースを作業フォルダにコピーした。
    • run_base_exF_before_1step/：SUBSTEPS=1（参照解）
    • run_base_exF_before_1step_version5/：!VERSION 5、SUBSTEPS=10
    • run_base_substep10_noversion/：上記から !VERSION 行のみを削除したもの
    • 共通条件：361要素、弾性（E=4000, ν=0.3）、!REFTEMP 20、初期温度20℃、!TEMPERATURE 120℃、!BOUNDARY FIX 1-3
  • 実行：各ディレクトリで OMP_NUM_THREADS=4 ../bin_base/fistr1。比較は python3 compare_steps.py <参照> <対象> で行った。
  • 結果：
    • run_base_exF_before_1step_version5：上表のとおり再現した。最終substepの参照解との最大相対差は DISPLACEMENT 3.60、NodalSTRESS 5.15。
    • run_base_substep10_noversion：!VERSION なし（VERSION 3 相当）では、線形静解析の !STEP, SUBSTEPS が無視され、1 substepで終了する。結果は参照解と完全に一致した。つまり本不具合は、!VERSION 5 でsubstep分割が有効になったときに顕在化する。

不具合修正の試行

原因に対する仮説

  1. 応力更新での熱ひずみの扱いが不整合（主因）
    • fstr_UpdateNewton（fistr1/src/analysis/static/fstr_Update.f90:137-146）では、基準温度 tt0 を材料で切り替えている。
      • 弾塑性/NORTON材料：前substep温度 last_temp → 熱ひずみは増分
      • それ以外（弾性など）：初期温度 → 熱ひずみは全量
    • Update_C3 等（fistr1/src/lib/static_LIB_3d.f90 ほか）では、dstrain = ε(u) − EPSTH(tt0→ttc) を作る一方で、gauss%strain = dstrain + EPSTH（全ひずみ）を保存し、fstr_UpdateState がこれを strain_bak にコピーする。
    • Update_Stress3D の INFINITESIMAL 分岐では σ = σ_bak + D (dstrain − strain_bak) となり、弾性材料では Δε − εth_{n+1} が使われる（正しくは Δε − (εth_{n+1} − εth_n)）。
    • その結果、内力が過小になって次substepの右辺が過大になる。線形解析はNewton反復1回で打ち切る（fstr_solve_NonLinear.f90:166）ため、誤差は補正されずに累積する。
    • 1次元モデル u_{n+1} = u_n − σ_n + a、σ_{n+1} = σ_n + (u_{n+1} − (n+1)a − u_n) で計算した比は 0.1, 0.2, 0.4, 0.7, 1.1, …, 4.6 となり、実測と完全に一致した。
    • 解決策として、次の2案を比較した。
      • 案A：線形時に tt0 = last_temp とする（fstr_Update.f90）。この案では、全量形で応力を計算する2D要素（UPDATE_C2、static_LIB_2d.f90:509）・UP要素（StressUpdate_up）・粘弾性を除外する条件分岐が必要になる。
      • 案B：Update_Stress3D の INFINITESIMAL 分岐で、前状態の熱ひずみ εth_n = EPSTH(tt0→ttn) を戻して機械ひずみ増分を作る。
      • 案Bは、影響範囲がこのルーチンを呼ぶ増分形の3D要素（C3/Bbar/IC/Fbar/SESNS）に限られる。弾塑性材料では tt0 = ttn なので εth_n = 0 となり、挙動は変わらない。ユーザーと協議のうえ、案Bを採用した。

仮説1の検証

  • 修正内容：Update_Stress3D の INFINITESIMAL 分岐で、ttn ≠ tt0 のとき Cal_Thermal_expansion_C3(tt0, ttn, …) により EPSTHN を求める。これを用いて機械ひずみ増分 dmstrain = dstrain − strain_bak + EPSTHN を作り、応力増分とひずみエネルギーの計算に使う。異方性の判定は Update_C3 と同様に MC_ORTHOEXP の有無で行う。
  • 結果（run_fix_exF_before_1step_version5/、バイナリ bin_fix/fistr1）：比が 0.100, 0.200, 0.300, …, 1.000 となり、累積は解消した。最終substepの参照解との最大相対差は DISPLACEMENT 1.66e-4、NodalSTRESS 3.98e-4、NodalMISES 1.58e-4 である。
  • SUBSTEPS=1 のケース（run_fix_exF_before_1step/）は修正前と完全に一致し、影響はない。

テスト結果

  • 解析実行ディレクトリ：
    • run_fix_exF_before_1step_version5/：不具合ケース。可視化出力は BIN_VTK で vis_out/ にある。
    • run_fix_exF_before_1step/：SUBSTEPS=1 ケース。
    • 結果は「仮説1の検証」のとおりである。pvbatch がない環境のため、画像による確認は省略した。
  • serialテスト（ctest --test-dir FrontISTR/build -L ',serial$'、ログ out_ctest_serial.log）：112件中111件がPASSした。
    • test_serial_solver/saamg のみ FAIL した。A361/A761/A641_mpc で SA-AMG が収束しない（indefinite preconditioner / 反復上限）。
    • 修正前のバイナリ（bin_base/fistr1）で同じテストを実行しても、同じ3メッシュで同様に失敗した（out_test_saamg_base.log）。したがって本修正とは無関係の既存事象である。
    • 熱荷重テスト analysis/static/exF と、要素ライブラリのテスト（3dIC, Bbar, Fbar, selectiveESNS, UP）、材料テスト（Elastoplastic, Viscoelastic など）はすべてPASSした。

test計算の可視化結果を添付します熱応力解析.pdf

Merge request reports