Fix Arrhenius shift factor calculation in TRS

inagakirequested to merge
feature/fix-trs-arrhenius into master

close #484 (closed)

Arrhenius型TRSシフトファクターのバグ修正

粘弾性材料の温度依存性を定義する!TRSキーワードのArrhenius型（DEFINITION=ARRHENIUS）において、シフトファクター計算に不正な除算が混入していたバグを修正した。

問題設定と解決策

Viscoelastic.f90のtrs関数およびtrsinc関数において、Arrhenius型のシフトファクター計算でmvar(4)による除算が行われていた。しかし、fstr_ctrl_get_TRSはmatval(1), matval(2), matval(3)の3パラメータのみを読み込むため、mvar(4)にはTRSとは無関係な値（M_THICK=!SECTIONのTHICKパラメータ）が格納されており、計算結果が不正となっていた。

Arrhenius型シフトファクターの定義は以下の通りである。

\ln(A) = \frac{E_0}{R}\left(\frac{1}{\theta - \theta^Z} - \frac{1}{\theta_0 - \theta^Z}\right)

ここで：

  • \theta_0: 参照温度（mvar(1)）
  • C_1 = E_0/R: 活性化エネルギーとガス定数の比（mvar(2)）
  • C_2 = \theta^Z: 絶対零度オフセット（mvar(3)）

この式には4番目のパラメータは存在しないため、/mvar(4)の除算を削除することで修正した。

入出力仕様

解析制御ファイルでArrhenius型TRSを使用するには、!VISCOELASTICの後に以下を記述する。

!TRS, DEFINITION=ARRHENIUS
50.0, 500.0, -273.15

データ行のパラメータ：θ₀（参照温度）, C₁=E₀/R（活性化エネルギー/ガス定数）, C₂=θ^Z（絶対零度オフセット）。入出力の仕様自体に変更はなく、既存の読み込み処理（fstr_ctrl_get_TRSで3パラメータ読み込み）はそのまま機能する。

機能の設計

既存設計

trs関数（シフトファクター計算）とtrsinc関数（シフトファクター増分計算）は、mtypeの値によりWLF型（VISCOELASTIC+1）とArrhenius型（VISCOELASTIC+2）を分岐する。mvar(:)はmaterials%variables配列へのポインタであり、mvar(1)〜mvar(3)にTRSパラメータが格納される。

fstr_ctrl_get_TRSはDEFINITIONパラメータの値（ipt: WLF=1, ARRHENIUS=2）を判別し、mattype = mattype + iptでセットする。TRSパラメータは3つ（matval(1), matval(2), matval(3)）のみ読み込む。

変更内容

Viscoelastic.f90のArrhenius分岐から/mvar(4)を削除した（3行）。WLF分岐のコードは変更なし。読み込み処理（fstr_ctrl_material.f90）への変更もなし。

【Viscoelastic.f90 — trs関数（修正後）】

    if(mtype==VISCOELASTIC+2) then   ! Arrhenius
      hsn = mvar(2)*( 1.d0/(tn-mvar(3))-1.d0/(mvar(1)-mvar(3)) )
    else   ! WLF

【Viscoelastic.f90 — trsinc関数（修正後）】

    if(mtype==VISCOELASTIC+2) then   ! Arrhenius
      hsn = -mvar(2)*( 1.d0/(tn-mvar(3))-1.d0/(mvar(1)-mvar(3)) )
      hsn1 = -mvar(2)*( 1.d0/(tn1-mvar(3))-1.d0/(mvar(1)-mvar(3)) )
    else   ! WLF

機能テスト

テスト仕様

単位立方体（8節点5要素tet4）・一軸ひずみ・定温のテストviscoe_thermal_trs_arrheniusを作成した。小ひずみ粘弾性の解析解と照合し、Arrhenius式の実装を検算する。

メッシュ・境界条件：

  • 単位立方体[0,1]³を5つのtet4要素で分割
  • x=0面: u_1 = 0、x=1面: u_1 = 0.0001（処方変位）
  • 全節点: u_2 = u_3 = 0（一軸ひずみ状態）
  • 熱膨張なし。温度は初期=適用=125°C一定（TRSのみに影響）

材料パラメータ：

  • E = 1000, \nu = 0.3, Prony級数1項: g_1 = 0.5, \tau_1 = 1.0
  • TRS Arrhenius: \theta_0 = 50°C, C_1 = 500, C_2 = -273.15

VISCOステップ: \Delta t = T = 1.0（単一ステップ）

解析解の導出（UpdateViscoelasticの計算手順を追跡）：

  1. シフトファクター: \text{trsinc}(125, 125) = \exp\left(500\left(\frac{1}{323.15} - \frac{1}{398.15}\right)\right) = \exp(0.2915) = 1.339
  2. 有効時間: \Delta t_{\text{eff}} = 1.339 \times 1.0 = 1.339
  3. 緩和関数: \Delta\tau = 1.339, H = (1 - e^{-1.339})/1.339 = 0.5512
  4. \Delta q = 0.5 \times 0.5512 = 0.2756, \mu_0 = 0.5
  5. G = 384.6, K = 833.3, \varepsilon_{11} = 0.0001
  6. \sigma_{11} = 2G(\mu_0 \cdot 2\varepsilon_{11}/3 + \Delta q \cdot 2\varepsilon_{11}/3) + K\varepsilon_{11} = 0.03978 + 0.08333 = 0.1231
  7. \sigma_{22} = \sigma_{33} = 0.0634, \sigma_{\text{Mises}} = 0.0597

理論的バウンド：瞬間弾性応答 \sigma_{11} = 0.1346、完全緩和 \sigma_{11} = 0.1090。期待値はバウンド内にあり、\theta > \theta_0で緩和促進（応力低下）する物理的傾向も正しい。

判定基準：

  • ソルバー出力が解析解と0.1%以内で一致すること
  • CMakeビルドがwarning/errorなく正常終了すること
  • ctest（serial, hybrid）が全件PASSすること

テスト結果

全て正常動作を確認した。

ビルド: CMakeビルドが正常終了。Viscoelastic.f90のコンパイルにwarning/errorなし。

解析解との照合: ソルバーが正常終了（FrontISTR Completed !!）。出力値と解析解の比較は以下の通り。

物理量 ソルバー出力 解析解 相対誤差
\sigma_{11} 1.2313E-01 0.1231 0.02%
\sigma_{22}, \sigma_{33} 6.3443E-02 0.0634 0.07%
\sigma_{\text{Mises}} 5.9684E-02 0.0597 0.02%

いずれも0.1%以内で一致。残差はNLSTATICの大変形定式化（対数ひずみ）と小ひずみ近似の差に起因する（u = 0.0001ではO(u^2)の影響が0.02%程度）。

作業ディレクトリ：work_trs_arrhenius/

リグレッションテスト:

  • serial: 63/63 テストPASS（viscoe_thermal_trs（WLF型）を含む）
  • hybrid: 57/57 テストPASS

妥当性確認テスト

機能テストが妥当性確認を兼ねる。UpdateViscoelasticルーチンの計算手順を追跡した小ひずみ解析解との照合（相対誤差0.1%以内）により、Arrhenius式の修正が正しいことを確認した。また、\theta > \theta_0で緩和が促進される物理的傾向と、応力が理論的バウンド（瞬間弾性〜完全緩和）内にあることも確認した。

Merge request reports