はじめに
前回までは、降伏応力を一定値として扱うモデル、塑性ひずみの蓄積によって降伏応力が直線的に大きくなるモデル、さらに降伏応力そのものが温度で低下するモデルを順番に見てきました。
これらは熱弾塑性解析の入口として分かりやすい一方で、実際の材料試験から得られる応力–ひずみ関係は、必ずしも一本の直線硬化で表せるとは限りません。
そこでこの記事では、材料の塑性変形中の抵抗をフローカーブとして与える考え方を扱います。フローカーブとは、塑性変形が進んだときに材料がどの程度の応力で流動するかを表す関係です。ここでは、材料点レベルの1次元モデルとして、
- 塑性変形が進むと流動応力が上がること
- 温度が上がると流動応力が下がること
- 現在の塑性履歴と温度から降伏判定を行うこと
- 塑性補正量 \(\Delta\gamma\) がどのように決まるか
を整理します。今回も、対象は両端固定された1次元棒です。モデルを複雑にしすぎないため、空間分布や有限要素法にはまだ進まず、1つの材料点の応答として理解します。
対象とする1次元棒
対象は、長さ \(L_0\) の1次元棒です。棒は両端で固定されているため、棒全体としての軸方向ひずみは
\varepsilon = 0 \tag{1}
\end{align}
です。この式は、棒が自由に伸び縮みできないことを表しています。温度が上がると材料は自由膨張しようとしますが、両端固定のため、その熱膨張は拘束されます。その結果、棒の内部には熱応力が発生します。
全ひずみは、弾性ひずみ、塑性ひずみ、熱ひずみに分けて考えます。
\varepsilon=\varepsilon_e+\varepsilon_p+\varepsilon_{\mathrm{th}} \tag{2}
\end{align}
この式は、観測される全体のひずみを、応力を生む弾性ひずみ、履歴として残る塑性ひずみ、温度変化だけで生じる熱ひずみに分けるものです。
熱ひずみは
\varepsilon_{\mathrm{th}}=\alpha \Delta T \tag{3}
\end{align}
で表します。ここで \(\alpha\) は線膨張係数、\(\Delta T\) は基準温度からの温度上昇です。両端固定条件である式(1)を式(2)へ代入すると、
\varepsilon_e=-\varepsilon_p-\varepsilon_{\mathrm{th}} \tag{4}
\end{align}
となります。つまり、熱ひずみが正に増えると、それを打ち消すように弾性ひずみは負、すなわち圧縮側になります。
応力はフックの法則
\sigma =E\varepsilon_e \tag{5}
\end{align}
で計算します。この関係により、弾性ひずみが圧縮側なら応力も圧縮側になります。
フローカーブとは何か
完全弾塑性モデルでは、塑性変形が始まる応力を一定の降伏応力 \(\sigma_Y\) として扱いました。線形硬化モデルでは、塑性変形が進むほど降伏応力が直線的に大きくなると考えました。フローカーブは、この考えをさらに一般化したものです。
フローカーブでは、塑性変形中に材料が示す抵抗応力を、蓄積した塑性変形量の関数として表します。
1次元では、蓄積した塑性変形量を \(\kappa\) と書きます。ただし、実際の材料試験で得られるフローカーブは、蓄積塑性ひずみだけで決まるとは限りません。温度だけでなく、塑性変形がどの速さで進むか、つまり塑性ひずみ速度にも依存することがあります。
その場合、本来は流動応力を
\sigma_{\mathrm{flow}}= f\left(\kappa,\dot{\kappa},\Delta T\right) \tag{6}
\end{align}
のように考えます。ここで \(\dot{\kappa}\) は蓄積塑性ひずみの時間変化率であり、塑性変形が速く進むほど材料抵抗が変わる効果を表します。この依存性まで入れると、応力を降伏面へ戻すだけではなく、塑性変形がどの速度で進むかを同時に扱う必要があります。つまり、速度依存塑性や粘塑性の考え方に近づきます。
この記事では、まずフローカーブを熱弾塑性解析へ入れる感覚をつかむため、ひずみ速度依存性には立ち入りません。準静的な計算を想定し、\(\dot{\kappa}\) の効果は固定または無視できるものとして、流動応力を \(\kappa\) と \(\Delta T\) の関数に簡略化します。具体的には、次の形で与えます。
\sigma_{\mathrm{flow}}(\kappa,\Delta T)=\max\left[\sigma_{\min},\left\{\sigma_{Y0}+Q\left(1-\exp(-b\kappa)\right)\right\}\left(1-\beta\Delta T\right)\right] \tag{7}
\end{align}
この式の中で、\(\sigma_{Y0}\) は基準温度での初期降伏応力です。
\(Q(1-\exp(-b\kappa))\) は、塑性変形が進むと流動応力が増える効果を表しています。ただし、増加量はいつまでも直線的に増えるのではなく、徐々に飽和します。これにより、実際の金属材料でよく見られる「最初は硬化しやすく、その後は硬化が緩やかになる」挙動を簡単に表せます。
一方、\((1-\beta\Delta T)\) は温度上昇による軟化を表します。温度が上がるほど、同じ塑性変形量 \(\kappa\) に対する流動応力が低下します。\(\sigma_{\min}\) は、単純な式を高温側へ外挿したときに流動応力が非現実的に小さくなりすぎることを避けるための下限値です。
図1を見ると、同じ \(\kappa\) でも温度が高いほど流動応力が小さいことが分かります。また、同じ温度では、\(\kappa\) が増えるほど流動応力が大きくなります。つまり、このモデルでは、塑性履歴による硬化と温度上昇による軟化が同時に働きます。
フローカーブから見たSSカーブ
フローカーブは、横軸を蓄積塑性ひずみ \(\kappa\)、縦軸を流動応力 \(\sigma_{\mathrm{flow}}\) とした材料データです。一方、一般に見る応力–ひずみ線図、いわゆるSSカーブは、横軸が全ひずみ、または機械ひずみです。両者は同じものではありませんが、単調引張や単調圧縮を考えれば、フローカーブからSSカーブの形を作ることができます。この節では熱ひずみを差し引いた機械ひずみ
\varepsilon_{\mathrm{m}}=\varepsilon-\varepsilon_{\mathrm{th}}=\varepsilon_e+\varepsilon_p \tag{8}
\end{align}
を使ってSSカーブを描きます。この式は、温度による自由膨張分を取り除き、応力と塑性変形に関わるひずみだけを見るためのものです。
では、フローカーブからSSカーブを作るときに、具体的に何をしているのかを見ておきます。ここでは、ある温度 \(\Delta T\) を固定し、その温度で単調に引張または圧縮する状況を考えます。
SSカーブを描くときには、まず横軸の値である機械ひずみ \(\varepsilon_{\mathrm{m}}\) を指定します。その指定した \(\varepsilon_{\mathrm{m}}\) に対して、縦軸の応力 \(\sigma\) を1点ずつ計算していきます。つまり、SSカーブ作成の基本は
\varepsilon_{\mathrm{m}}\quad\longrightarrow\quad\sigma \tag{9}
\end{align}
という対応関係を作ることです。この式は、新しい物理法則ではありません。SSカーブを描くために、横軸を先に決め、そこから縦軸を計算するという作業の流れを表しています。
ただし、塑性域では \(\sigma\) を直接計算できません。なぜなら、応力は流動応力で決まりますが、その流動応力は蓄積塑性ひずみ \(\kappa\) に依存するからです。したがって、塑性域ではまず \(\kappa\) を求め、その後で \(\sigma\) を求めます。単調負荷では、塑性ひずみの向きは途中で反転しないので、機械ひずみの符号を
s=\mathrm{sign}\left(\varepsilon_{\mathrm{m}}\right) \tag{10}
\end{align}
として扱えます。
まず、弾性範囲から考えます。弾性範囲では塑性ひずみがまだ発生していないので、
\kappa=0,\qquad\varepsilon_p=0 \tag{11}
\end{align}
です。したがって、機械ひずみはそのまま弾性ひずみであり、応力は単純に
\sigma=E\varepsilon_{\mathrm{m}} \tag{12}
\end{align}
です。この式は、指定した横軸 \(\varepsilon_{\mathrm{m}}\) にヤング率を掛ければ、そのまま縦軸 \(\sigma\) が得られることを表します。
ただし、この弾性予測が許されるのは、応力の大きさが初期状態の流動応力を超えない場合だけです。その判定は
E\left|\varepsilon_{\mathrm{m}}\right|\le\sigma_{\mathrm{flow}}\left(0,\Delta T\right) \tag{13}
\end{align}
で行います。この式の左辺は、横軸に選んだ機械ひずみをすべて弾性ひずみだと仮定したときの応力の大きさです。右辺は、その温度でまだ塑性ひずみがない状態の流動応力です。左辺が右辺以下なら、まだ降伏していないので、式(12)でSSカーブ上の点が決まります。
一方、この条件を超えると塑性域です。塑性域では、指定した機械ひずみ \(\varepsilon_{\mathrm{m}}\) の一部が塑性ひずみに変わり、残りだけが弾性ひずみとして応力を生みます。単調負荷では、蓄積塑性ひずみ \(\kappa\) を塑性ひずみの大きさとして見られるため、
|\varepsilon_p|=\kappa \tag{14}
\end{align}
と置けます。このとき、指定した機械ひずみの大きさ \(|\varepsilon_{\mathrm{m}}|\) は、弾性ひずみの大きさ \(|\varepsilon_e|\) と塑性ひずみの大きさ \(\kappa\) に分かれます。
|\varepsilon_{\mathrm{m}}| = |\varepsilon_e|+\kappa \tag{15}
\end{align}
したがって、弾性ひずみの大きさは
|\varepsilon_e| =\left(|\varepsilon_{\mathrm{m}}|-\kappa\right) \tag{16}
\end{align}
です。応力は弾性ひずみにヤング率 \(E\) を掛けたものなので、応力の大きさは
|\sigma|=E\left(|\varepsilon_{\mathrm{m}}|-\kappa\right) \tag{17}
\end{align}
と書けます。この式は、\(\kappa\) が大きくなるほど、同じ機械ひずみの中で弾性ひずみとして残る分が小さくなり、その分だけ弾性応力も小さくなることを表しています。
一方、塑性変形が続いている点では、応力の大きさはフローカーブ上の流動応力と一致していなければなりません。
|\sigma|=\sigma_{\mathrm{flow}}\left(\kappa,\Delta T\right) \tag{18}
\end{align}
この式は、塑性域の応力が材料の流動応力に制限されることを表します。式(17)と式(18)を同時に満たす必要があるので、SSカーブ上の塑性点では
E\left(|\varepsilon_{\mathrm{m}}|-\kappa\right)=\sigma_{\mathrm{flow}}\left(\kappa,\Delta T\right) \tag{19}
\end{align}
を満たす \(\kappa\) を求めます。左辺は、全体の機械ひずみから塑性ひずみ分を差し引いた後に残る弾性応力です。右辺は、その塑性履歴 \(\kappa\) における材料の流動応力です。この2つが一致するところが、フローカーブと弾性則が整合する点です。ここでの未知数は \(\sigma\) ではなく、まず \(\kappa\) です。
横軸 \(\varepsilon_{\mathrm{m}}\) と温度 \(\Delta T\) はすでに指定済みなので、式(19)を満たす \(\kappa\) を1次元の方程式として解きます。
\kappa=\kappa\left(|\varepsilon_{\mathrm{m}}|,\Delta T\right) \tag{20}
\end{align}
この式は、機械ひずみと温度を決めると、その点で必要な蓄積塑性ひずみが決まる、という意味です。実際の計算では、ニュートン法などを使って式(19)を解きます。求めた \(\kappa\) をフローカーブへ戻せば、応力の大きさが得られます。
最後に、引張側か圧縮側かを符号 \(s\) で戻して、
\sigma=s\,\sigma_{\mathrm{flow}}\left(\kappa,\Delta T\right) \tag{21}
\end{align}
として描けます。つまり、フローカーブからSSカーブを作る操作は、各ひずみ点ごとに次の順番で行います。
- 横軸の値 \(\varepsilon_{\mathrm{m}}\) と温度 \(\Delta T\) を指定します。
- 式(13)で、弾性範囲か塑性範囲かを判定します。
- 弾性範囲なら、式(12)により \(\sigma=E\varepsilon_{\mathrm{m}}\) とします。
- 塑性範囲なら、式(19)を満たす \(\kappa\) を求めます。
- 求めた \(\kappa\) を式(21)へ代入し、SSカーブ上の応力 \(\sigma\) を求めます。
このように、フローカーブはそのままSSカーブではありません。フローカーブは \(\kappa\) と \(\sigma_{\mathrm{flow}}\) の関係であり、SSカーブは \(\varepsilon_{\mathrm{m}}\) と \(\sigma\) の関係です。両者の間をつなぐ役割をしているのが、機械ひずみの分解式と、式(19)で表される整合条件です。
図２に、ここまでの操作を温度ごとに繰り返して作成したSSカーブを記載します。ただし、これは単調負荷を仮定した説明用のSSカーブです。加熱・冷却を受ける熱拘束問題の履歴そのものは、後の応力更新アルゴリズムで別に計算します。
原点近くの傾きはどの温度でも同じです。これは、ヤング率 \(E\) を温度によらず一定としたためです。一方で、塑性域に入った後の応力レベルは温度によって変わります。温度が高いほど材料は低い応力で流動するため、SSカーブ全体が低い応力側へ移ります。
降伏判定
現在の応力が、現在の温度と塑性履歴で決まる流動応力を超えると、材料は塑性状態に入ります。その判定は、1次元では次のように書けます。
f(\sigma,\kappa,\Delta T)=|\sigma|-\sigma_{\mathrm{flow}}(\kappa,\Delta T)\le 0 \tag{22}
\end{align}
絶対値が入っているのは、引張側でも圧縮側でも、応力の大きさが流動応力に達すれば塑性変形が進むと考えるためです。両端固定棒を加熱すると、最初に問題になるのは圧縮側の降伏です。
ただし冷却過程では、過去に発生した塑性ひずみの影響により、引張側へ応力が戻ることもあります。
応力更新アルゴリズム
数値計算では、温度履歴を小さなステップに分け、各ステップで応力と内部変数を更新します。フローカーブを使う場合も、基本はこれまでと同じです。
ただし今回は、流動応力が蓄積塑性ひずみ \(\kappa\) に依存するため、どの量が前ステップの値で、どの量が今回ステップの値なのかを区別する必要があります。そこで、前ステップを \(n\)、今回更新した後のステップを \(n+1\) と書きます。
まず、今回の温度ステップで塑性ひずみがまだ増えないと仮定します。この仮定のもとでは、塑性ひずみは前ステップの値 \(\varepsilon_p^n\) のままです。一方、温度はすでに今回ステップの温度 \(\Delta T_{n+1}\) まで進めるので、熱ひずみは \(\varepsilon_{\mathrm{th}}^{n+1}\) になります。したがって、試行応力は次のように計算します。
\sigma^{\mathrm{trial}}=E\left(\varepsilon^{n+1}-\varepsilon_p^n-\varepsilon_{\mathrm{th}}^{n+1}\right) \tag{23}
\end{align}
これまでの記事では、読みやすさを優先して \(\sigma^{\mathrm{trial}}=E(\varepsilon-\varepsilon_p-\varepsilon_{\mathrm{th}})\) のようにステップ番号を省いて書いていました。式(23)は、それと同じ関係をステップ計算用に書き直したものです。つまり、\(\varepsilon_p\) を「前ステップから持ち越した塑性ひずみ」、\(\varepsilon_{\mathrm{th}}\) を「今回の温度で決まる熱ひずみ」と明示しただけです。
両端固定棒では全ひずみが常に拘束されるため、数値例では \(\varepsilon^{n+1}=0\) として計算します。この式は、「前ステップまでの塑性ひずみを固定したまま、温度だけを更新したら応力はいくつになるか」を予測するものです。
次に、前ステップの蓄積塑性ひずみ \(\kappa_n\) を使って、試行降伏関数を評価します。
f^{\mathrm{trial}}=\left|\sigma^{\mathrm{trial}}\right|-\sigma_{\mathrm{flow}}\left(\kappa_n,\Delta T_{n+1}\right) \tag{24}
\end{align}
\(f^{\mathrm{trial}}\le0\) であれば、試行応力は現在のフローカーブの内側にあります。したがって、そのステップは弾性として扱い、
\sigma^{n+1}=\sigma^{\mathrm{trial}} \tag{25}
\end{align}
\end{align}
とします。つまり、応力だけが弾性的に変化し、塑性履歴は変わりません。
塑性の場合
\(f^{\mathrm{trial}}>0\) の場合、試行応力はフローカーブの外側に出ています。このままでは降伏条件を満たさないため、塑性ひずみを増やして応力をフローカーブ上へ戻します。ここで使う補正量が、塑性乗数増分 \(\Delta\gamma\) です。
まず、試行応力の符号を
s=\mathrm{sign}\left(\sigma^{\mathrm{trial}}\right) \tag{27}
\end{align}
とします。\(s=1\) なら引張側、\(s=-1\) なら圧縮側です。
塑性補正を \(\Delta\gamma\) だけ入れると、応力は弾性的に戻されます。
\sigma^{n+1}=\sigma^{\mathrm{trial}}-E\Delta\gamma s \tag{28}
\end{align}
この式は、塑性ひずみが増えた分だけ弾性ひずみが小さくなり、その分だけ応力が試行応力から戻ることを表しています。圧縮側の試行応力であれば \(s=-1\) なので、応力は圧縮側から原点方向へ戻ります。
同時に、塑性ひずみと蓄積塑性ひずみは
\varepsilon_p^{n+1}=\varepsilon_p^n+\Delta\gamma s \tag{29}
\end{align}
\kappa_{n+1} =\kappa_n+\Delta\gamma \tag{30}
\end{align}
と更新します。\(\varepsilon_p\) は向きを持つ塑性ひずみなので符号 \(s\) が付きます。
一方、\(\kappa\) は塑性変形がどれだけ蓄積したかを表す量なので、常に増加します。
塑性乗数増分の決め方
塑性補正後の応力は、更新後の流動応力と一致していなければなりません。したがって、
\left|\sigma^{n+1}\right|=\sigma_{\mathrm{flow}}\left(\kappa_{n+1},\Delta T_{n+1}\right) \tag{31}
\end{align}
が成り立ちます。式(28)と式(30)をこの条件へ代入すると、
\left|\sigma^{\mathrm{trial}}\right|-E\Delta\gamma=\sigma_{\mathrm{flow}}\left(\kappa_n+\Delta\gamma,\Delta T_{n+1}\right) \tag{32}
\end{align}
となります。左辺は、試行応力から \(\Delta\gamma\) だけ弾性的に戻した後の応力の大きさです。右辺は、塑性変形が \(\Delta\gamma\) だけ進んだ後の流動応力です。この2つが一致する \(\Delta\gamma\) が、そのステップで必要な塑性補正量です。
線形硬化ではこの式を手で解けました。しかし、今回のフローカーブは指数関数を含むため、一般には数値的に解きます。
残差関数を
R(\Delta\gamma)=\left|\sigma^{\mathrm{trial}}\right|-E\Delta\gamma-\sigma_{\mathrm{flow}}\left(\kappa_n+\Delta\gamma,\Delta T_{n+1}\right) \tag{33}
\end{align}
とおくと、求めたい条件は \(R(\Delta\gamma)=0\) です。この \(R(\Delta\gamma)\) は、式(32)の左辺と右辺の差です。つまり、試行応力から弾性的に戻した応力レベルと、その時点の流動応力がどれだけずれているかを表します。
\(R(\Delta\gamma)>0\) であれば、まだ応力の戻し量が足りず、\(R(\Delta\gamma)<0\) であれば戻しすぎです。したがって、残差がゼロになる位置が、ちょうど整合条件を満たす塑性乗数増分になります。
ニュートン法を使うなら、例えば
\Delta\gamma\leftarrow\Delta\gamma+\frac{R(\Delta\gamma)}{E+\displaystyle\left.\frac{\partial\sigma_{\mathrm{flow}}}{\partial\kappa}\right|_{\kappa=\kappa_n+\Delta\gamma}} \tag{34}
\end{align}
のように更新します。分母に \(E\) だけでなく \(\partial\sigma_{\mathrm{flow}}/\partial\kappa\) が入るのは、塑性補正によって応力が戻るだけでなく、同時に材料側の流動応力も硬化によって上がるためです。
図3では、横軸が \(\Delta\gamma\) です。青線は \(\Delta\gamma\) を大きくするほど下がります。これは塑性ひずみを増やすほど、弾性ひずみが減り、応力が戻るためです。
赤線は \(\Delta\gamma\) を大きくするほど少し上がります。これは、塑性変形が進むほど \(\kappa\) が増え、フローカーブ上の流動応力が上がるためです。
黒い点では、青線で表した応力レベルと赤線で表した流動応力が一致します。この一致点が、そのステップでの塑性補正量です。計算上は、この一致を残差 \(R(\Delta\gamma)=0\) として解いている、と考えればよいです。
最後に、更新後の弾性ひずみは
\varepsilon_e^{n+1}=\frac{\sigma^{n+1}}{E} \tag{35}
\end{align}
として計算します。
弾性ひずみは、塑性ひずみを更新した後に、最終的に得られた応力と整合する量として決まります。
数値例の解析条件
ここまでで、フローカーブを用いるための定式化を整理しました。ここからは、同じ考え方を両端固定棒の加熱・冷却問題に適用します。
温度は \(0\) K から \(400\) K まで上げ、その後 \(0\) K まで戻します。
\begin{array}{c}
\text{表: 数値例で用いる材料定数と解析条件。} \\
\begin{array}{ccl}
\hline
\text{記号} & \text{値} & \text{意味} \\
\hline
E & 210000\,\mathrm{MPa} & \text{ヤング率} \\
\alpha & 12\times10^{-6} 1/K & \text{線膨張係数} \\
\sigma_{Y0} & 250\,\mathrm{MPa} & \text{基準温度での初期降伏応力} \\
Q & 180\,\mathrm{MPa} & \text{飽和硬化量} \\
b & 350 & \text{硬化の進みやすさ} \\
\beta & 8.0\times10^{-4} 1/K & \text{流動応力の温度低下率} \\
\Delta T_{\max} & 400\,\mathrm{K} & \text{最大温度上昇} \\
\hline
\end{array}
\end{array}
\end{align}
今回の例では、塑性変形が進むほど材料は硬化します。一方で、温度が上がるほど材料は軟化します。したがって、加熱中の応力履歴は、熱膨張による圧縮応力、塑性履歴による硬化、温度上昇による軟化の釣り合いとして決まります。
計算結果
代表値の確認
この条件で得られる代表値を表2に示します。降伏開始温度は約 \(93.3\) K です。完全弾塑性で降伏応力を一定にした場合の約 \(99.2\) K より少し低くなります。これは、温度上昇によって初期の流動応力が下がるためです。
\begin{array}{c}
\text{表: フローカーブを用いた数値例の代表値。} \\
\begin{array}{cl}
\hline
\text{項目} & \text{値} \\
\hline
\text{降伏開始温度} & \text{約}93.3\,\mathrm{K} \\
\text{最大温度での応力} & -257.4\,\mathrm{MPa} \\
\text{最大温度での蓄積塑性ひずみ} & 3.57\times10^{-3} \\
\text{最大温度での流動応力} & 257.4\,\mathrm{MPa} \\
\text{冷却後の残留応力} & 401.2\,\mathrm{MPa} \\
\text{冷却後の塑性ひずみ} & -1.91\times10^{-3} \\
\text{冷却後の蓄積塑性ひずみ} & 5.24\times10^{-3} \\
\text{冷却後の流動応力} & 401.2\,\mathrm{MPa} \\
\hline
\end{array}
\end{array}
\end{align}
最大温度では、温度上昇による軟化のために流動応力は下がります。それでも完全弾塑性の温度依存モデルより圧縮応力の絶対値が大きくなるのは、加熱中に塑性変形が進み、\(\kappa\) が増えて材料が硬化するためです。つまり、温度は材料を軟化させますが、塑性履歴は材料を硬化させます。今回の結果は、この2つが同時に働いた結果です。
冷却後の残留応力が \(401\) MPa 程度まで大きくなる点も重要です。完全弾塑性モデルでは、冷却後の応力は基準温度での降伏応力付近に戻りました。しかしフローカーブを使うと、冷却後には温度低下により軟化が消え、さらに加熱・冷却中に蓄積した \(\kappa\) による硬化が残ります。そのため、最終的な流動応力そのものが大きくなり、残留応力も高くなります。
履歴図とアニメーションによる確認
図5では、加熱中に応力が圧縮側へ進み、降伏後は負側の流動応力に沿って変化します。
冷却に入ると熱ひずみが減るため、応力は引張側へ戻ります。このとき、すでに塑性履歴が蓄積しているため、単に来た道を戻るのではありません。流動応力の図が一つの線ではなく履歴を持つように見えるのは、横軸の温度だけでは状態を決めきれないためです。現在の流動応力は、温度 \(\Delta T\) と履歴変数 \(\kappa\) の両方で決まります。
図6では、温度ごとのSSカーブの上に、熱拘束で得られた応力経路を重ねます。背景のSSカーブは材料側の性質を表し、移動する点は熱拘束問題の応答を表します。両者を重ねると、材料データとしてのフローカーブと、境界条件によって生じる応力履歴を分けて見やすくなります。
図7では、熱ひずみ、塑性ひずみ、弾性ひずみを温度に対して分けて表示します。両端固定では全ひずみがゼロなので、これらの成分は互いに打ち消し合います。特に、塑性ひずみが履歴として残るため、同じ温度に戻っても弾性ひずみと応力が元の値に戻らないことが確認できます。
ひずみ成分の確認
図8を見ると、熱ひずみは温度に対して単純な直線です。一方、塑性ひずみと弾性ひずみは加熱時と冷却時で異なる経路を通ります。これは、塑性変形が履歴として残るためです。
冷却後に温度が \(0\) K に戻ると、熱ひずみはゼロになります。しかし塑性ひずみはゼロに戻りません。両端固定では全ひずみがゼロなので、式(2)より、残った塑性ひずみを打ち消す弾性ひずみが必要になります。その弾性ひずみが残留応力を生みます。
このモデルから分かること
フローカーブを使うと、降伏応力を一つの値として扱うのではなく、塑性履歴と温度に応じて変わる材料抵抗として扱えるようになります。今回の例では、温度上昇は流動応力を下げ、塑性変形の蓄積は流動応力を上げます。熱弾塑性解析では、この二つが同時に起こるため、応力履歴は単純な直線にはなりません。
特に重要なのは、状態を温度だけで決められないことです。同じ温度でも、加熱の途中なのか、冷却の途中なのか、どれだけ塑性変形を経験した後なのかによって、現在の流動応力と応力状態は変わります。そのため、フローカーブを使う解析では、応力や塑性ひずみだけでなく、蓄積塑性ひずみ \(\kappa\) のような内部変数を更新し続ける必要があります。
また、非線形のフローカーブを使うと、塑性補正量 \(\Delta\gamma\) は閉じた式で簡単に求まらない場合があります。その場合は、整合条件を満たすように数値的に \(\Delta\gamma\) を解きます。これは実務的な材料モデルへ進むための重要な入口です。
フローカーブのひずみ速度依存性
本来のフローカーブは、式(6)で触れたように、蓄積塑性ひずみ \(\kappa\)、塑性ひずみ速度 \(\dot{\kappa}\)、温度 \(\Delta T\) に依存して与えられることがあります。この場合、たとえば同じ温度、同じ塑性ひずみ量でも、ゆっくり変形させた場合と速く変形させた場合で流動応力が変わります。
\sigma_{\mathrm{flow}}= f(\kappa,\dot{\kappa},\Delta T)\quad\Longrightarrow\quad\text{塑性変形の進む速さも未知量として扱う} \tag{36}
\end{align}
ただし、ひずみ速度依存性を入れると、塑性変形が「降伏条件を超えたら即座に戻す」問題から、「どの速度で塑性変形が進むか」を解く問題へ変わります。時間刻みの取り方や粘塑性の発展式も関係してくるため、定式化は一段複雑になります。
そのため本記事では、そこへは立ち入らず、まずは温度依存フローカーブを使った応力更新に集中しました。この形を理解しておくと、次にひずみ速度依存性を加えるときにも、何が新しく増える未知量なのかを整理しやすくなります。
MATLABコード
今回使用したコードを以下に記載しました。
% flow_curve_material_point_demo.m
%
% 温度依存フローカーブを用いた 1 次元熱弾塑性解析のサンプルコードです。
% 両端固定された棒を 0 K -> 400 K -> 0 K と加熱・冷却し、応力、
% 塑性ひずみ、蓄積塑性ひずみ、流動応力を逐次更新します。
%
% 流動応力は一定の降伏応力ではなく、蓄積塑性ひずみ kappa と温度上昇 dT
% に依存する量として、次式で与えます。
%
% sigma_flow(kappa, dT)
% = [sigmaY0 + Q * (1 - exp(-b * kappa))] * (1 - beta * dT)
%
% この 1 ファイルを実行すると、記事用の静止画 PNG とアニメーション MP4 を
% 保存します。
clear; close all; clc;
%% 1. 解析条件
% 単位系は MPa, mm, K を想定します。
E = 210000; % MPa, Young's modulus
alpha = 12e-6; % 1/K, coefficient of linear thermal expansion
sigmaY0 = 250; % MPa, initial yield stress at reference temperature
Q = 180; % MPa, saturation hardening amount
b = 350; % -, saturation hardening rate with kappa
beta = 8.0e-4; % 1/K, thermal softening coefficient
flowStressLower = 50; % MPa, lower bound used only for numerical robustness
Tmax = 400; % K, maximum temperature rise
nHalf = 240; % number of increments during heating
% アニメーションを作らない場合は false にします。
saveMp4 = true;
videoFrameRate = 15;
animationStride = 4;
responseMp4FileName = "flow_curve_response_animation.mp4";
stressStrainMp4FileName = "flow_curve_stress_strain_path.mp4";
strainTemperatureMp4FileName = "flow_curve_strain_temperature_path.mp4";
returnMappingMp4FileName = "flow_curve_return_mapping_animation.mp4";
% 温度履歴: 0 -> Tmax -> 0
dTUp = linspace(0, Tmax, nHalf + 1);
dTDown = linspace(Tmax, 0, nHalf + 1);
dT = [dTUp, dTDown(2:end)];
step = 0:numel(dT)-1;
%% 2. 材料点解析
% 応力更新アルゴリズムを全温度ステップに適用します。
result = run_flow_curve_case(dT, E, alpha, sigmaY0, Q, b, beta, flowStressLower);
fprintf('=== Flow curve material point ===\n');
fprintf('Initial yield temperature = %.2f K\n', result.initialYieldTemperature);
fprintf('Stress at Tmax = %.3f MPa\n', result.stressAtTmax);
fprintf('Accumulated plastic strain at Tmax = %.6e\n', result.kappaAtTmax);
fprintf('Flow stress at Tmax = %.3f MPa\n', result.flowStressAtTmax);
fprintf('Final stress = %.3f MPa\n', result.sigma(end));
fprintf('Final plastic strain = %.6e\n', result.epsP(end));
fprintf('Final accumulated plastic strain = %.6e\n', result.kappa(end));
fprintf('Final flow stress = %.3f MPa\n', result.sigmaFlow(end));
%% 3. 静止画の作成
% Figure 1: 温度ごとのフローカーブ
temperatureLevels = [0, 100, 200, 400];
kappaAxis = linspace(0, 0.010, 600);
colors = [0.10 0.30 0.75;
0.10 0.55 0.70;
0.90 0.55 0.15;
0.85 0.20 0.15];
figFlow = figure('Color', 'w', 'Position', [120 120 1200 720]);
hold on;
for k = 1:numel(temperatureLevels)
Tlevel = temperatureLevels(k);
sigmaFlowCurve = flow_stress(kappaAxis, Tlevel, sigmaY0, Q, b, beta, flowStressLower);
plot(kappaAxis * 1000, sigmaFlowCurve, 'LineWidth', 3.0, ...
'Color', colors(k, :), ...
'DisplayName', sprintf('\\DeltaT = %d K', Tlevel));
end
grid on; box on;
xlabel('accumulated plastic strain \kappa (10^{-3})');
ylabel('flow stress \sigma_{flow} (MPa)');
title('temperature-dependent flow curves');
legend('Location', 'southeast', 'FontSize', 14);
set(gca, 'FontSize', 18, 'LineWidth', 1.0);
xlim([0 10]);
ylim([130 430]);
ax = gca;
text(1.55, 405, {'larger \kappa', 'hardening'}, ...
'FontSize', 15, 'Color', [0.25 0.25 0.25], ...
'BackgroundColor', 'w', 'Margin', 3, ...
'HorizontalAlignment', 'left', 'VerticalAlignment', 'middle');
draw_annotation_arrow(ax, 2.95, 405, 3.75, 405);
draw_annotation_arrow(ax, 7.2, 405, 7.2, 305);
text(7.45, 355, {'higher \DeltaT', 'lower flow stress'}, ...
'FontSize', 15, 'Color', [0.25 0.25 0.25], ...
'BackgroundColor', 'w', 'Margin', 3, ...
'HorizontalAlignment', 'left', 'VerticalAlignment', 'middle');
exportgraphics(figFlow, 'flow_curve_material_curves.png', 'Resolution', 220);
% Figure 2: フローカーブから作る応力-ひずみ線図
figSS = figure('Color', 'w', 'Position', [140 140 1200 720]);
hold on;
epsAxis = linspace(-0.006, 0.006, 1201);
for k = 1:numel(temperatureLevels)
Tlevel = temperatureLevels(k);
sigmaCurve = monotonic_stress_strain_curve(epsAxis, E, Tlevel, sigmaY0, Q, b, beta, flowStressLower);
plot(epsAxis * 1000, sigmaCurve, 'LineWidth', 2.8, ...
'Color', colors(k, :), ...
'DisplayName', sprintf('\\DeltaT = %d K', Tlevel));
end
plot([0 0], [-450 450], 'k-', 'LineWidth', 0.8, 'HandleVisibility', 'off');
plot(xlim, [0 0], 'k-', 'LineWidth', 0.8, 'HandleVisibility', 'off');
grid on; box on;
xlabel('mechanical strain \epsilon - \epsilon_{th} (10^{-3})');
ylabel('stress \sigma (MPa)');
title('stress-strain curves generated from the flow curves');
legend('Location', 'northwest', 'FontSize', 13);
set(gca, 'FontSize', 18, 'LineWidth', 1.0);
xlim([-6 6]);
ylim([-430 430]);
exportgraphics(figSS, 'flow_curve_stress_strain_curves.png', 'Resolution', 220);
% Figure 3: 解析結果の全体像
figSummary = figure('Color', 'w', 'Position', [100 100 1300 980]);
tiledlayout(2, 2, 'Padding', 'loose', 'TileSpacing', 'loose');
nexttile;
plot(step, dT, 'k-', 'LineWidth', 2.3);
grid on; box on;
xlabel('step');
ylabel('\DeltaT (K)');
title('temperature history');
set(gca, 'FontSize', 16, 'LineWidth', 1.0);
nexttile;
plot(dT, result.sigma, 'Color', [0.9 0.25 0.2], 'LineWidth', 2.4); hold on;
plot(dT, result.sigmaFlow, '--', 'Color', [0.35 0.35 0.35], 'LineWidth', 1.8);
plot(dT, -result.sigmaFlow, '--', 'Color', [0.35 0.35 0.35], 'LineWidth', 1.8);
grid on; box on;
xlabel('\DeltaT (K)');
ylabel('\sigma (MPa)');
title('stress - temperature');
legend({'stress', '\pm flow stress'}, 'Location', 'southoutside', ...
'Orientation', 'horizontal', 'FontSize', 12);
set(gca, 'FontSize', 16, 'LineWidth', 1.0);
ylim([-440 440]);
nexttile;
plot(dT, result.sigmaFlow, 'Color', [0.1 0.45 0.2], 'LineWidth', 2.4);
grid on; box on;
xlabel('\DeltaT (K)');
ylabel('\sigma_{flow} (MPa)');
title('current flow stress');
set(gca, 'FontSize', 16, 'LineWidth', 1.0);
nexttile;
plot(dT, result.kappa * 1000, 'Color', [0.2 0.45 0.9], 'LineWidth', 2.4);
grid on; box on;
xlabel('\DeltaT (K)');
ylabel('\kappa (10^{-3})');
title('accumulated plastic strain');
set(gca, 'FontSize', 16, 'LineWidth', 1.0);
exportgraphics(figSummary, 'flow_curve_summary.png', 'Resolution', 220);
% Figure 4: 熱ひずみ、塑性ひずみ、弾性ひずみ
figStrain = figure('Color', 'w', 'Position', [120 120 1200 620]);
plot(dT, result.epsTh * 1000, 'Color', [0.1 0.45 0.2], 'LineWidth', 2.4); hold on;
plot(dT, result.epsP * 1000, 'Color', [0.9 0.25 0.2], 'LineWidth', 2.4);
plot(dT, result.epsE * 1000, 'Color', [0.2 0.45 0.9], 'LineWidth', 2.4);
grid on; box on;
xlabel('\DeltaT (K)');
ylabel('strain (10^{-3})');
title('strain components');
legend({'thermal strain \epsilon_{th}', 'plastic strain \epsilon_p', 'elastic strain \epsilon_e'}, ...
'Location', 'best', 'FontSize', 14);
set(gca, 'FontSize', 17, 'LineWidth', 1.0);
exportgraphics(figStrain, 'flow_curve_strain_components.png', 'Resolution', 220);
% Figure 5: 塑性乗数増分 Delta gamma の決まり方
plasticIndices = find(result.dGamma > 0);
conceptIndex = plasticIndices(round(0.45 * numel(plasticIndices)));
figReturn = plot_return_mapping_concept(result, conceptIndex, E, sigmaY0, Q, b, beta, flowStressLower);
exportgraphics(figReturn, 'flow_curve_return_mapping_concept.png', 'Resolution', 220);
%% 4. アニメーションの作成
% Animation 1: 温度履歴、応力、流動応力、蓄積塑性ひずみの全体像
if saveMp4
nStep = numel(dT);
figAnim = figure('Name', 'Flow curve response animation', ...
'Color', 'w', 'Position', [100 100 1300 980]);
writer = VideoWriter(responseMp4FileName, 'MPEG-4');
writer.FrameRate = videoFrameRate;
open(writer);
for i = 1:animationStride:nStep
clf(figAnim);
tiledlayout(figAnim, 2, 2, 'Padding', 'loose', 'TileSpacing', 'loose');
nexttile;
plot(step, dT, 'k-', 'LineWidth', 2.3); hold on;
plot(step(i), dT(i), 'o', 'MarkerFaceColor', 'k', ...
'MarkerEdgeColor', 'w', 'MarkerSize', 8);
grid on; box on;
xlabel('step');
ylabel('\DeltaT (K)');
title(sprintf('temperature history: step %d / %d', i, nStep));
xlim([step(1) step(end)]);
ylim([0 Tmax * 1.08]);
set(gca, 'FontSize', 16, 'LineWidth', 1.0);
nexttile;
plot(dT, result.sigma, 'Color', [0.9 0.25 0.2], 'LineWidth', 2.4); hold on;
plot(dT, result.sigmaFlow, '--', 'Color', [0.35 0.35 0.35], 'LineWidth', 1.8);
plot(dT, -result.sigmaFlow, '--', 'Color', [0.35 0.35 0.35], 'LineWidth', 1.8);
plot(dT(i), result.sigma(i), 'o', 'MarkerFaceColor', [0.9 0.25 0.2], ...
'MarkerEdgeColor', 'w', 'MarkerSize', 8);
grid on; box on;
xlabel('\DeltaT (K)');
ylabel('\sigma (MPa)');
title(sprintf('stress - temperature: \\sigma = %.1f MPa', result.sigma(i)));
ylim([-440 440]);
legend({'stress', '\pm flow stress'}, 'Location', 'southoutside', ...
'Orientation', 'horizontal', 'FontSize', 12);
set(gca, 'FontSize', 16, 'LineWidth', 1.0);
nexttile;
plot(dT, result.sigmaFlow, 'Color', [0.1 0.45 0.2], 'LineWidth', 2.4); hold on;
plot(dT(i), result.sigmaFlow(i), 'o', 'MarkerFaceColor', [0.1 0.45 0.2], ...
'MarkerEdgeColor', 'w', 'MarkerSize', 8);
grid on; box on;
xlabel('\DeltaT (K)');
ylabel('\sigma_{flow} (MPa)');
title(sprintf('flow stress: %.1f MPa', result.sigmaFlow(i)));
set(gca, 'FontSize', 16, 'LineWidth', 1.0);
nexttile;
plot(dT, result.kappa * 1000, 'Color', [0.2 0.45 0.9], 'LineWidth', 2.4); hold on;
plot(dT(i), result.kappa(i) * 1000, 'o', 'MarkerFaceColor', [0.2 0.45 0.9], ...
'MarkerEdgeColor', 'w', 'MarkerSize', 8);
grid on; box on;
xlabel('\DeltaT (K)');
ylabel('\kappa (10^{-3})');
title(sprintf('accumulated plastic strain: %.3f x10^{-3}', result.kappa(i) * 1000));
set(gca, 'FontSize', 16, 'LineWidth', 1.0);
drawnow;
writeVideo(writer, getframe(figAnim));
end
close(writer);
fprintf('Saved response animation MP4: %s\n', responseMp4FileName);
end
% Animation 2: 応力-ひずみ線図上の応力経路
if saveMp4
nStep = numel(dT);
epsMech = result.epsE + result.epsP;
xPath = epsMech * 1000;
xLim = [-6.2, 6.2];
yLim = [-450, 450];
epsAxisAnim = linspace(xLim(1), xLim(2), 1000) / 1000;
figSSAnim = figure('Name', 'Flow curve stress-strain path animation', ...
'Color', 'w', 'Position', [140 140 1200 720]);
writer = VideoWriter(stressStrainMp4FileName, 'MPEG-4');
writer.FrameRate = videoFrameRate;
open(writer);
for i = 1:animationStride:nStep
clf(figSSAnim);
hold on;
for k = 1:numel(temperatureLevels)
Tlevel = temperatureLevels(k);
sigmaCurve = monotonic_stress_strain_curve(epsAxisAnim, E, Tlevel, sigmaY0, Q, b, beta, flowStressLower);
plot(epsAxisAnim * 1000, sigmaCurve, 'LineWidth', 2.4, ...
'Color', colors(k, :), ...
'DisplayName', sprintf('\\DeltaT = %d K', Tlevel));
end
plot(xPath, result.sigma, '-', 'Color', [0.65 0.65 0.65], ...
'LineWidth', 1.4, 'DisplayName', 'thermal loading path');
plot(xPath(1:i), result.sigma(1:i), '-', ...
'Color', [0.05 0.25 0.55], 'LineWidth', 2.4, 'HandleVisibility', 'off');
plot(xPath(i), result.sigma(i), 'o', ...
'MarkerFaceColor', [0.9 0.25 0.1], 'MarkerEdgeColor', 'k', ...
'MarkerSize', 9, 'LineWidth', 1.0, 'HandleVisibility', 'off');
plot([0 0], yLim, 'k-', 'LineWidth', 0.8, 'HandleVisibility', 'off');
plot(xLim, [0 0], 'k-', 'LineWidth', 0.8, 'HandleVisibility', 'off');
grid on; box on;
xlabel('mechanical strain \epsilon - \epsilon_{th} (10^{-3})');
ylabel('stress \sigma (MPa)');
title(sprintf('stress-strain path: step %d / %d, \\DeltaT = %.1f K', i, nStep, dT(i)));
xlim(xLim);
ylim(yLim);
legend('Location', 'northwest', 'FontSize', 12);
set(gca, 'FontSize', 18, 'LineWidth', 1.0);
drawnow;
writeVideo(writer, getframe(figSSAnim));
end
close(writer);
fprintf('Saved stress-strain path MP4: %s\n', stressStrainMp4FileName);
end
% Animation 3: 各ひずみ成分と温度の関係
if saveMp4
nStep = numel(dT);
heatingIndex = 1:(nHalf + 1);
coolingIndex = (nHalf + 1):nStep;
figStrainAnim = figure('Name', 'Flow curve strain-temperature animation', ...
'Color', 'w', 'Position', [160 100 1100 900]);
writer = VideoWriter(strainTemperatureMp4FileName, 'MPEG-4');
writer.FrameRate = videoFrameRate;
open(writer);
for i = 1:animationStride:nStep
clf(figStrainAnim);
tiledlayout(figStrainAnim, 3, 1, 'TileSpacing', 'compact', 'Padding', 'compact');
plot_strain_tile(dT, result.epsTh, i, heatingIndex, coolingIndex, ...
'\epsilon_{th} (10^{-3})', 'thermal strain vs temperature', [0.1 0.45 0.2], true);
plot_strain_tile(dT, result.epsP, i, heatingIndex, coolingIndex, ...
'\epsilon_p (10^{-3})', 'plastic strain vs temperature', [0.9 0.25 0.2], false);
plot_strain_tile(dT, result.epsE, i, heatingIndex, coolingIndex, ...
'\epsilon_e (10^{-3})', 'elastic strain vs temperature', [0.2 0.45 0.9], false);
xlabel('\DeltaT (K)');
drawnow;
writeVideo(writer, getframe(figStrainAnim));
end
close(writer);
fprintf('Saved strain-temperature path MP4: %s\n', strainTemperatureMp4FileName);
end
% Animation 4: 塑性乗数増分 Delta gamma の反復的な決まり方
if saveMp4
plasticIndices = find(result.dGamma > 0);
selected = plasticIndices(1:animationStride:end);
figReturnAnim = figure('Name', 'Flow curve return mapping animation', ...
'Color', 'w', 'Position', [140 140 1200 720]);
writer = VideoWriter(returnMappingMp4FileName, 'MPEG-4');
writer.FrameRate = videoFrameRate;
open(writer);
fixedLimits = return_mapping_fixed_limits(result, selected, E, sigmaY0, Q, b, beta, flowStressLower);
for jj = 1:numel(selected)
i = selected(jj);
clf(figReturnAnim);
plot_return_mapping_axes(result, i, E, sigmaY0, Q, b, beta, flowStressLower, false, fixedLimits);
title(sprintf('plastic correction: step %d, \\DeltaT = %.1f K', i, dT(i)));
drawnow;
writeVideo(writer, getframe(figReturnAnim));
end
close(writer);
fprintf('Saved return mapping animation MP4: %s\n', returnMappingMp4FileName);
end
%% 5. 解析用の関数
function result = run_flow_curve_case(dT, E, alpha, sigmaY0, Q, b, beta, flowStressLower)
% 温度履歴 dT に沿って、両端固定棒の応力と内部変数を更新します。
% 全ひずみは 0 とし、熱ひずみ・塑性ひずみ・弾性ひずみの和で拘束条件を満たします。
n = numel(dT);
sigma = zeros(1, n);
epsP = zeros(1, n);
epsE = zeros(1, n);
epsTh = zeros(1, n);
kappa = zeros(1, n);
sigmaFlow = zeros(1, n);
dGamma = zeros(1, n);
yielded = false(1, n);
sigmaTrialHistory = zeros(1, n);
initialYieldTemperature = NaN;
for i = 1:n
% 今回ステップの熱ひずみ
epsTh(i) = alpha * dT(i);
% 前ステップで確定した塑性ひずみと蓄積塑性ひずみ
if i == 1
epsPold = 0;
kappaOld = 0;
else
epsPold = epsP(i - 1);
kappaOld = kappa(i - 1);
end
% 試行応力:
% 塑性ひずみを前ステップ値に固定し、温度だけを今回値へ進めた応力です。
sigmaTrial = E * (0 - epsPold - epsTh(i));
sigmaTrialHistory(i) = sigmaTrial;
% 試行降伏判定
sigmaFlowTrial = flow_stress(kappaOld, dT(i), sigmaY0, Q, b, beta, flowStressLower);
fTrial = abs(sigmaTrial) - sigmaFlowTrial;
if fTrial <= 1e-10
% 弾性ステップ: 塑性履歴は更新しません。
sigma(i) = sigmaTrial;
epsP(i) = epsPold;
kappa(i) = kappaOld;
sigmaFlow(i) = sigmaFlowTrial;
else
% 塑性ステップ: return mapping により Delta gamma を求めます。
yielded(i) = true;
if isnan(initialYieldTemperature)
initialYieldTemperature = dT(i);
end
s = sign(sigmaTrial);
gamma = solve_plastic_multiplier(abs(sigmaTrial), kappaOld, dT(i), ...
E, sigmaY0, Q, b, beta, flowStressLower);
% 応力、塑性ひずみ、蓄積塑性ひずみを確定します。
dGamma(i) = gamma;
epsP(i) = epsPold + gamma * s;
kappa(i) = kappaOld + gamma;
sigma(i) = sigmaTrial - E * gamma * s;
sigmaFlow(i) = flow_stress(kappa(i), dT(i), sigmaY0, Q, b, beta, flowStressLower);
end
epsE(i) = sigma(i) / E;
end
[~, idxTmax] = max(dT);
result.sigma = sigma;
result.epsP = epsP;
result.epsE = epsE;
result.epsTh = epsTh;
result.kappa = kappa;
result.sigmaFlow = sigmaFlow;
result.dGamma = dGamma;
result.yielded = yielded;
result.sigmaTrial = sigmaTrialHistory;
result.initialYieldTemperature = initialYieldTemperature;
result.stressAtTmax = sigma(idxTmax);
result.kappaAtTmax = kappa(idxTmax);
result.flowStressAtTmax = sigmaFlow(idxTmax);
end
function sigmaFlow = flow_stress(kappa, dT, sigmaY0, Q, b, beta, lowerBound)
% 温度依存フローカーブ。
% kappa が増えると硬化し、dT が増えると熱軟化します。
thermalFactor = 1 - beta * dT;
sigmaFlow = (sigmaY0 + Q * (1 - exp(-b * kappa))) .* thermalFactor;
sigmaFlow = max(lowerBound, sigmaFlow);
end
function gamma = solve_plastic_multiplier(absSigmaTrial, kappaOld, dT, E, sigmaY0, Q, b, beta, lowerBound)
% 塑性乗数増分 Delta gamma をニュートン法で解きます。
%
% 解きたい式は
%
% R(gamma)
% = absSigmaTrial - E * gamma
% - sigma_flow(kappaOld + gamma, dT)
% = 0
%
% です。gamma が増えると、応力は E*gamma だけ戻り、同時に
% kappaOld + gamma が増えるので流動応力も変化します。
gamma = 0;
for iter = 1:30
sigmaFlow = flow_stress(kappaOld + gamma, dT, sigmaY0, Q, b, beta, lowerBound);
residual = absSigmaTrial - E * gamma - sigmaFlow;
if abs(residual) < 1e-10
break;
end
hardeningSlope = flow_slope(kappaOld + gamma, dT, Q, b, beta, lowerBound, sigmaFlow);
gamma = gamma + residual / (E + hardeningSlope);
gamma = max(gamma, 0);
end
end
function slope = flow_slope(kappa, dT, Q, b, beta, lowerBound, sigmaFlow)
% d sigma_flow / d kappa。
% lowerBound に張り付いた場合は、下限値を一定値として扱うため傾き 0 とします。
if sigmaFlow <= lowerBound + 1e-12
slope = 0;
else
slope = Q * b * exp(-b * kappa) * (1 - beta * dT);
end
end
function sigmaCurve = monotonic_stress_strain_curve(epsAxis, E, dT, sigmaY0, Q, b, beta, lowerBound)
% 指定した温度 dT における単調負荷の応力-ひずみ線図を作ります。
% フローカーブは kappa の関数なので、各ひずみ点で return mapping と
% 同じ考え方により応力を求めています。
sigmaCurve = zeros(size(epsAxis));
for i = 1:numel(epsAxis)
epsValue = epsAxis(i);
s = sign(epsValue);
absEps = abs(epsValue);
sigmaY = flow_stress(0, dT, sigmaY0, Q, b, beta, lowerBound);
if E * absEps <= sigmaY
sigmaCurve(i) = E * epsValue;
else
gamma = solve_plastic_multiplier(E * absEps, 0, dT, E, sigmaY0, Q, b, beta, lowerBound);
sigmaCurve(i) = s * (E * absEps - E * gamma);
end
end
end
%% 6. 描画用の補助関数
function draw_annotation_arrow(ax, x1, y1, x2, y2)
color = [0.25 0.25 0.25];
[nx1, ny1] = data_to_normalized(ax, x1, y1);
[nx2, ny2] = data_to_normalized(ax, x2, y2);
annotation('arrow', [nx1 nx2], [ny1 ny2], ...
'Color', color, 'LineWidth', 1.8, ...
'HeadLength', 10, 'HeadWidth', 10);
end
function [nx, ny] = data_to_normalized(ax, x, y)
pos = ax.Position;
xl = ax.XLim;
yl = ax.YLim;
nx = pos(1) + (x - xl(1)) / (xl(2) - xl(1)) * pos(3);
ny = pos(2) + (y - yl(1)) / (yl(2) - yl(1)) * pos(4);
end
function fig = plot_return_mapping_concept(result, i, E, sigmaY0, Q, b, beta, flowStressLower)
fig = figure('Color', 'w', 'Position', [140 140 1200 720]);
plot_return_mapping_axes(result, i, E, sigmaY0, Q, b, beta, flowStressLower, true);
title('plastic multiplier increment from the return mapping equation');
end
function fixedLimits = return_mapping_fixed_limits(result, selected, E, sigmaY0, Q, b, beta, flowStressLower)
plasticIndices = find(result.dGamma > 0);
gammaMax = max(result.dGamma(plasticIndices)) * 1.08;
gammaAxis = linspace(0, max(gammaMax, 1e-5), 500);
yValues = [];
for idx = selected(:).'
dTNow = result.epsTh(idx) / 12e-6;
gammaRoot = result.dGamma(idx);
kappaOld = result.kappa(idx) - gammaRoot;
absSigmaTrial = abs(result.sigmaTrial(idx));
leftLine = absSigmaTrial - E * gammaAxis;
flowLine = flow_stress(kappaOld + gammaAxis, dTNow, sigmaY0, Q, b, beta, flowStressLower);
yValues = [yValues; leftLine(:); flowLine(:)]; %#ok<AGROW>
end
fixedLimits.x = [-0.03 * gammaMax * 1000, gammaMax * 1000];
fixedLimits.y = [floor((min(yValues) - 10) / 10) * 10, ceil((max(yValues) + 10) / 10) * 10];
end
function plot_return_mapping_axes(result, i, E, sigmaY0, Q, b, beta, flowStressLower, compactView, fixedLimits)
if nargin < 9
compactView = false;
end
if nargin < 10
fixedLimits = [];
end
dTNow = result.epsTh(i) / 12e-6;
gammaRoot = result.dGamma(i);
kappaOld = result.kappa(i) - gammaRoot;
absSigmaTrial = abs(result.sigmaTrial(i));
plasticIndices = find(result.dGamma > 0);
gammaMaxFull = max(result.dGamma(plasticIndices)) * 1.08;
if compactView
gammaMax = max(1.38 * gammaRoot, 3.2e-5);
else
gammaMax = gammaMaxFull;
end
gammaAxis = linspace(0, max(gammaMax, 1e-5), 500);
leftLine = absSigmaTrial - E * gammaAxis;
flowLine = flow_stress(kappaOld + gammaAxis, dTNow, sigmaY0, Q, b, beta, flowStressLower);
if compactView
yMinFixed = 260;
yMaxFixed = 275;
elseif ~isempty(fixedLimits)
yMinFixed = fixedLimits.y(1);
yMaxFixed = fixedLimits.y(2);
else
yMinFixed = min([leftLine(:); flowLine(:)]) - 12;
yMaxFixed = max([leftLine(:); flowLine(:)]) + 16;
end
hStress = plot(gammaAxis * 1000, leftLine, 'Color', [0.05 0.25 0.55], 'LineWidth', 2.8); hold on;
hFlow = plot(gammaAxis * 1000, flowLine, 'Color', [0.9 0.25 0.2], 'LineWidth', 2.8);
yRoot = flow_stress(kappaOld + gammaRoot, dTNow, sigmaY0, Q, b, beta, flowStressLower);
hSolution = plot(gammaRoot * 1000, yRoot, 'o', 'MarkerFaceColor', 'k', ...
'MarkerEdgeColor', 'w', 'MarkerSize', 9);
ylabel('stress level (MPa)');
ylim([yMinFixed yMaxFixed]);
xline(gammaRoot * 1000, ':', 'Color', [0.2 0.2 0.2], 'LineWidth', 1.4);
yArrow = yMinFixed + 0.20 * (yMaxFixed - yMinFixed);
xLeftTip = 0;
xRightTip = gammaRoot * 1000;
xMaxPlot = max(gammaAxis) * 1000;
headLength = min(0.014 * xMaxPlot, 0.28 * max(xRightTip - xLeftTip, eps));
headHeight = 0.018 * (yMaxFixed - yMinFixed);
plot([xLeftTip + headLength, xRightTip - headLength], [yArrow yArrow], '-', ...
'Color', [0.15 0.15 0.15], 'LineWidth', 1.4, 'HandleVisibility', 'off');
patch([xLeftTip, xLeftTip + headLength, xLeftTip + headLength], ...
[yArrow, yArrow + headHeight, yArrow - headHeight], ...
[0.15 0.15 0.15], 'EdgeColor', [0.15 0.15 0.15], 'HandleVisibility', 'off');
patch([xRightTip, xRightTip - headLength, xRightTip - headLength], ...
[yArrow, yArrow + headHeight, yArrow - headHeight], ...
[0.15 0.15 0.15], 'EdgeColor', [0.15 0.15 0.15], 'HandleVisibility', 'off');
text(0.5 * (xLeftTip + xRightTip), yArrow + 0.055 * (yMaxFixed - yMinFixed), '\Delta\gamma', ...
'FontSize', 18, 'HorizontalAlignment', 'center', 'BackgroundColor', 'w', 'Margin', 3);
blueLabelX = 0.18 * xMaxPlot;
blueLabelY = interp1(gammaAxis * 1000, leftLine, blueLabelX) + 0.055 * (yMaxFixed - yMinFixed);
redLabelX = 0.60 * xMaxPlot;
redLabelY = interp1(gammaAxis * 1000, flowLine, redLabelX) + 0.075 * (yMaxFixed - yMinFixed);
text(blueLabelX, blueLabelY, ...
'|\sigma^{trial}| - E\Delta\gamma', 'FontSize', 16, ...
'Color', [0.05 0.25 0.55], 'BackgroundColor', 'w', 'Margin', 3, ...
'Rotation', -8);
text(redLabelX, redLabelY, ...
'\sigma_{flow}(\kappa_n+\Delta\gamma,\DeltaT)', 'FontSize', 16, ...
'Color', [0.9 0.25 0.2], 'BackgroundColor', 'w', 'Margin', 3, ...
'Rotation', 2);
grid on; box on;
xlabel('plastic multiplier increment \Delta\gamma (10^{-3})');
legend([hStress hFlow hSolution], ...
{'stress after elastic unloading', 'flow curve value', 'solution'}, ...
'Location', 'northeast', 'FontSize', 12);
set(gca, 'FontSize', 17, 'LineWidth', 1.0);
if ~isempty(fixedLimits)
xlim(fixedLimits.x);
else
xlim([-0.03 * xMaxPlot, xMaxPlot]);
end
end
function plot_strain_tile(dT, strain, i, heatingIndex, coolingIndex, yLabelText, titleText, color, isFirst)
nexttile;
plot(dT(heatingIndex), strain(heatingIndex) * 1000, '-', ...
'Color', color, 'LineWidth', 2.2); hold on;
plot(dT(coolingIndex), strain(coolingIndex) * 1000, '--', ...
'Color', color, 'LineWidth', 2.2);
plot(dT(i), strain(i) * 1000, 'o', 'MarkerFaceColor', 'k', ...
'MarkerEdgeColor', 'w', 'MarkerSize', 8);
yline(0, ':', 'Color', [0.35 0.35 0.35], 'LineWidth', 1.0);
grid on; box on;
xlim([0 max(dT)]);
ylabel(yLabelText);
if isFirst
title(sprintf('%s: step %d, \\DeltaT = %.1f K', titleText, i, dT(i)));
else
title(titleText);
end
set(gca, 'FontSize', 14, 'LineWidth', 1.0);
endまとめ
この記事では、1次元の両端固定棒を題材に、温度依存フローカーブを用いた熱弾塑性解析を整理しました。
フローカーブは、塑性変形中の材料抵抗を表す関係であり、単なる初期降伏応力よりも広い情報を持っています。塑性変形が進むと \(\kappa\) が増え、流動応力は上がります。一方、温度が上がると材料は軟化し、流動応力は下がります。
この二つの効果が同時に働くため、加熱・冷却の応力履歴は、完全弾塑性モデルや線形硬化モデルよりも履歴依存性が強くなります。冷却後に温度が元へ戻っても、塑性履歴は残ります。その履歴によって流動応力も変化するため、残留応力の大きさも初期降伏応力だけでは判断できません。
フローカーブを材料点レベルで扱えるようになると、次にFEMへ拡張したときにも、各積分点で何を更新しているのかが見えやすくなります。今回のモデルは簡単ですが、実務的な熱弾塑性解析へ進むための足場になります。
なお、実際の高温変形や高速変形では、フローカーブがひずみ速度にも依存する場合があります。その場合は、流動応力を \(\sigma_{\mathrm{flow}}=f(\kappa,\dot{\kappa},\Delta T)\) として扱い、塑性変形の進む速さも含めて解く必要があります。これは速度依存塑性や粘塑性の領域に入るため、本記事では扱いませんでした。
まずは今回のように \(\kappa\) と \(\Delta T\) への依存を理解しておくと、次の拡張に進みやすくなります。
コメント