はじめに
これまでの記事では、両端を完全に固定した1次元棒を主な対象としてきました。完全拘束では、温度上昇によって材料が自由に伸びようとしても、棒全体の伸びが許されません。そのため、熱ひずみが弾性ひずみとして打ち消され、熱応力が発生します。
しかし、実際の構造物がいつも完全拘束されているとは限りません。自由に伸びられる場合もあれば、周囲の部材や治具によって部分的に拘束される場合もあります。このとき、同じ材料、同じ温度履歴であっても、発生する熱応力や塑性ひずみは大きく変わります。
この記事では、1次元棒を対象として、
- 自由膨張
- 完全拘束
- ばね拘束
を比較します。目的は、材料モデルそのものを複雑にすることではなく、熱応力は材料定数だけでなく、拘束条件によって決まることを確認することです。
熱応力と拘束条件
温度が \(\Delta T\) だけ上がると、棒は自由であれば熱ひずみによって伸びます。線膨張係数を \(\alpha\) とすると、自由熱ひずみは
\varepsilon_{\mathrm{th}}=\alpha \Delta T \tag{1}
\end{align}
です。この式は、応力がなくても温度変化だけで生じるひずみを表しています。長さ \(L_0\) の棒であれば、自由に伸びたい量は
u_{\mathrm{free}}=L_0 \alpha \Delta T \tag{2}
\end{align}
です。つまり、温度が上がるほど、棒は自然に長くなろうとします。
一方、実際に許される伸びを \(u\) とすると、棒の全ひずみは
\varepsilon=\frac{u}{L_0} \tag{3}
\end{align}
です。この全ひずみは、弾性ひずみ、塑性ひずみ、熱ひずみに分けて
\varepsilon=\varepsilon_e+\varepsilon_p+\varepsilon_{\mathrm{th}} \tag{4}
\end{align}
と書けます。この式は、実際に観測される伸びが、応力を生む弾性ひずみ、履歴として残る塑性ひずみ、温度だけで生じる熱ひずみの和であることを表します。したがって、弾性ひずみは
\varepsilon_e =\varepsilon-\varepsilon_p-\varepsilon_{\mathrm{th}} \tag{5}
\end{align}
です。応力はフックの法則により
\sigma=E\varepsilon_e=E\left(\varepsilon- \varepsilon_p- \varepsilon_{\mathrm{th}}\right) \tag{6}
\end{align}
となります。この式を見ると、熱応力を決めているのは熱ひずみだけではないことが分かります。
実際の全ひずみ \(\varepsilon\)、つまり棒がどれだけ伸びることを許されるかも同じくらい重要です。
図2では、黄色の棒を自由膨張したい長さ、色付きの棒を実際に許された長さとして表示します。自由膨張では両者が一致します。完全拘束では、温度が上がって黄色の棒が伸びても、実際の棒の長さは固定されたままです。ばね拘束では、完全固定よりは伸びますが、自由膨張ほどは伸びません。この「伸びたい量」と「実際に伸びた量」の差が、熱応力を生む原因になります。
3つの拘束条件
自由膨張
自由膨張では、棒は熱によって伸びたいだけ伸びることができます。塑性ひずみがまだない状態であれば、
\varepsilon=\varepsilon_{\mathrm{th}} \tag{7}
\end{align}
となります。このとき式(5)より、弾性ひずみはゼロです。したがって応力もゼロになります。熱ひずみは発生していますが、それが拘束されていないため、応力にはなりません。
完全拘束
完全拘束では、棒の両端が固定されているため、伸びは
u = 0,\qquad\varepsilon = 0 \tag{8}
\end{align}
です。この場合、温度が上がっても棒全体は伸びられません。塑性ひずみがまだなければ、式(6)から
\sigma=-E\alpha\Delta T \tag{9}
\end{align}
となります。符号が負になるのは、棒が伸びたい方向と逆向きに圧縮応力が生じるためです。つまり完全拘束は、自由熱膨張をすべて応力に変える極端な条件です。
ばね拘束
ばね拘束では、棒の端部が外部ばねにつながっていると考えます。棒が伸びると、ばねはそれに抵抗する反力を生じます。棒の断面積を \(A\)、外部ばね剛性を \(k_s\) とすると、棒の内力とばね反力のつり合いは
\sigma A +k_s u=0 \tag{10}
\end{align}
です。この式は、棒が圧縮応力を持つとき、その内力とばねの引き戻す力がつり合うことを表しています。\(u=L_0\varepsilon\) を代入すると、
\sigma A+k_s L_0 \varepsilon=0 \tag{11}
\end{align}
となります。完全拘束では \(\varepsilon=0\) を先に決めましたが、ばね拘束では式(11)を満たすように \(\varepsilon\) が決まります。つまり、応力更新だけでなく、棒がどれだけ伸びるかも同時に決める必要があります。
拘束剛性の見方
ばね拘束の強さは、棒自身の軸剛性 \(EA/L_0\) と比べると分かりやすくなります。ここで出てくる \(EA/L_0\) は、棒を軸方向のばねとして見たときの剛性です。なぜこの形になるのかを確認しておきます。
棒に軸力 \(N\) が作用し、長さ \(L_0\) の棒が \(u\) だけ伸びるとします。このとき、棒の平均ひずみは
\varepsilon=\frac{u}{L_0} \tag{12}
\end{align}
です。線形弾性であれば、応力は
\sigma=E\varepsilon=E\frac{u}{L_0} \tag{13}
\end{align}
となります。軸力は応力に断面積を掛けたものなので、
N=\sigma A=\frac{EA}{L_0}u \tag{14}
\end{align}
です。
一方、ばねの力と変位の関係は
N =k_{\mathrm{bar}} u \tag{15}
\end{align}
と書けます。式(14)と式(15)を比べると、棒自身の軸方向ばね剛性は
k_{\mathrm{bar}}=\frac{EA}{L_0} \tag{16}
\end{align}
になります。つまり、\(EA/L_0\) は、棒を端部から引っ張ったときに、単位変位あたりどれだけ軸力が発生するかを表す量です。断面積 \(A\) やヤング率 \(E\) が大きいほど棒は硬く、長さ \(L_0\) が長いほど棒は柔らかくなります。そこで、無次元の拘束剛性比を
r=\frac{k_s}{EA/L_0}=\frac{k_s L_0}{EA} \tag{17}
\end{align}
と定義します。\(r=0\) なら外部ばねがないので自由膨張に近づきます。\(r\) が大きくなるほど外部ばねが硬くなり、完全拘束に近づきます。
今回の数値例では、中間的な拘束として \(r=1\) を使います。これは、外部ばねの剛性が棒自身の軸剛性と同程度であることを意味します。
弾塑性応力更新
拘束条件を変えても、材料点で行う応力更新の考え方はこれまでと同じです。違うのは、各温度ステップで使う全ひずみ \(\varepsilon\) の決め方です。ここを急いでしまうと、同じ応力更新式を使っているのに、なぜ拘束条件だけで結果が変わるのかが見えにくくなります。そこでこの節では、まず拘束条件ごとに全ひずみの決まり方を整理し、その後で共通の弾塑性応力更新を説明します。ばね拘束の具体的な数値例は、応力更新式を定義した後の小節で述べます。
応力更新に入る前に決めるもの
数値計算では、温度履歴を小さなステップに分け、前ステップを \(n\)、今回更新するステップを \(n+1\) と書きます。応力更新へ入る前に必要なのは、今回ステップの全ひずみ \(\varepsilon^{n+1}\) です。
全ひずみが決まると、前ステップまでの塑性ひずみ \(\varepsilon_p^n\) と今回の熱ひずみ \(\varepsilon_{\mathrm{th}}^{n+1}\) を使って、塑性補正前の試行応力を計算できます。
\sigma^{\mathrm{trial}}=E\left(\varepsilon^{n+1}-\varepsilon_p^n-\varepsilon_{\mathrm{th}}^{n+1}\right) \tag{18}
\end{align}
試行応力とは、今回ステップの全ひずみと熱ひずみは使う一方で、塑性ひずみだけは前ステップ値のまま固定して計算した応力です。この式は、「今回の拘束条件で決まった全ひずみに対して、まだ塑性補正を入れなければ応力はいくつになるか」を予測するものです。
したがって、拘束条件ごとの差は、式(18)そのものではなく、ここへ代入する \(\varepsilon^{n+1}\) の違いとして現れます。
自由膨張の場合
自由膨張では、外部から棒の伸びを妨げる力がありません。したがって、棒は熱ひずみと塑性ひずみを含めて、応力がゼロになるように伸びます。応力がゼロであることは、弾性ひずみがゼロであることと同じです。
\sigma^{n+1}=E\varepsilon_e^{n+1}= 0\qquad \Longrightarrow\qquad\varepsilon_e^{n+1}=0 \tag{19}
\end{align}
全ひずみの分解式
\varepsilon^{n+1}=\varepsilon_e^{n+1}+\varepsilon_p^{n+1}+\varepsilon_{\mathrm{th}}^{n+1} \tag{20}
\end{align}
に式(19)を入れると、
\varepsilon^{n+1}=\varepsilon_p^{n+1}+\varepsilon_{\mathrm{th}}^{n+1} \tag{21}
\end{align}
となります。ただし、応力更新を始める時点では、まだ今回ステップの塑性ひずみ \(\varepsilon_p^{n+1}\) は分かっていません。そこで数値計算では、前ステップの塑性ひずみ \(\varepsilon_p^n\) を使って、まず
\varepsilon_{\mathrm{free}}^{n+1}=\varepsilon_p^n+\varepsilon_{\mathrm{th}}^{n+1} \tag{22}
\end{align}
と置きます。この全ひずみを式(18)へ入れると、
\sigma_{\mathrm{free}}^{\mathrm{trial}}=E\left( \varepsilon_{\mathrm{free}}^{n+1}- \varepsilon_p^n-\varepsilon_{\mathrm{th}}^{n+1}\right)=0 \tag{23}
\end{align}
となります。つまり、自由膨張では試行応力の時点でゼロになり、降伏判定に進んでも通常は塑性化しません。熱ひずみは確かに発生していますが、棒がその分だけ伸びることを許されるため、弾性ひずみに変換されないのです。
完全拘束の場合
完全拘束では、棒の両端が固定されているため、実際の伸びは常にゼロです。したがって、今回ステップの全ひずみは材料状態に関係なく
\varepsilon_{\mathrm{fixed}}^{n+1}=0 \tag{24}
\end{align}
です。この条件を式(18)へ入れると、
\sigma_{\mathrm{fixed}}^{\mathrm{trial}}=E\left(0 -\varepsilon_p^n-\varepsilon_{\mathrm{th}}^{n+1}\right)=-E\left(\varepsilon_p^n+\varepsilon_{\mathrm{th}}^{n+1}\right) \tag{25}
\end{align}
となります。最初の加熱では \(\varepsilon_p^n=0\) なので、
\sigma_{\mathrm{fixed}}^{\mathrm{trial}}=-E\alpha\Delta T_{n+1} \tag{26}
\end{align}
です。これは、自由熱膨張したい分をすべて弾性ひずみで打ち消すことを意味します。
温度が上がるほど \(\varepsilon_{\mathrm{th}}^{n+1}\) が大きくなり、試行応力は圧縮側へ大きくなります。この試行応力が降伏応力を超えると、塑性補正によって応力が降伏面上へ戻され、塑性ひずみが更新されます。
ばね拘束の場合
ばね拘束では、自由膨張のように応力をゼロにすることも、完全拘束のように伸びをゼロにすることも、どちらも先に決められません。棒が少し伸びればばね反力が発生し、その反力と棒の内力がつり合う位置で伸びが決まります。したがって、今回ステップでは
\sigma^{n+1} A +k_s L_0 \varepsilon^{n+1}=0 \tag{27}
\end{align}
を満たす全ひずみ \(\varepsilon^{n+1}\) を求めます。ここで \(\sigma^{n+1}\) は、仮に選んだ \(\varepsilon^{n+1}\) から応力更新を通して決まる量です。そのため、ばね拘束では「全ひずみを決める処理」と「応力を更新する処理」が互いに結びつきます。この結びつきの具体的な解き方は、共通の応力更新式を定義した後、図3を使って説明します。
3つの拘束条件と応力更新の関係
ここまでをまとめると、応力更新に入る前の全ひずみは次のように扱います。
\varepsilon^{n+1}=\begin{cases}\varepsilon_p^n+\varepsilon_{\mathrm{th}}^{n+1},
& \text{自由膨張} \\
0,& \text{完全拘束} \\
g(\varepsilon^{n+1})=0\text{を満たす値},\qquad& \text{ばね拘束}
\end{cases} \tag{28}
\end{align}
この式は、材料モデルの違いではなく、境界条件の違いを表しています。ばね拘束の行に出てくる \(g(\varepsilon)\) は、後で式(40)として定義するつり合い残差です。この全ひずみが決まった後に、同じ応力更新アルゴリズムを適用します。
ここまでの流れを、計算全体の手順として見ると図3のようになります。3つの拘束条件で異なるのは、今回ステップの全ひずみ \(\varepsilon^{n+1}\) を決める部分です。その後の試行応力、降伏判定、塑性補正は共通の応力更新として扱えます。
全ひずみが決まった後の共通応力更新
全ひずみが決まったら、式(18)で試行応力を計算し、その大きさが現在の降伏応力を超えるかどうかを判定します。等方線形硬化では、現在の降伏応力を
\sigma_Y^n=\sigma_{Y0}+H\kappa_n \tag{29}
\end{align}
とします。\(\sigma_{Y0}\) は初期降伏応力、\(H\) は等方線形硬化係数、\(\kappa_n\) は前ステップまでの蓄積塑性ひずみです。試行降伏関数は
f^{\mathrm{trial}}=\left|\sigma^{\mathrm{trial}}\right|-\left(\sigma_{Y0}+H\kappa_n\right) \tag{30}
\end{align}
です。\(f^{\mathrm{trial}}\le0\) なら弾性ステップとして扱います。\(f^{\mathrm{trial}}>0\) なら塑性補正が必要です。このとき、塑性乗数増分は
\Delta\gamma=\frac{f^{\mathrm{trial}}}{E+H} \tag{31}
\end{align}
です。分母に \(E\) と \(H\) が入るのは、塑性補正によって応力が弾性的に戻る効果と、塑性変形の蓄積によって降伏応力が上がる効果が同時に働くためです。
応力、塑性ひずみ、蓄積塑性ひずみは
s=\mathrm{sign}\left(\sigma^{\mathrm{trial}}\right) \tag{32}
\end{align}
を用いて
\sigma^{n+1}=\sigma^{\mathrm{trial}}-E\Delta\gamma s \tag{33}
\end{align}
\varepsilon_p^{n+1}=\varepsilon_p^n +\Delta\gamma s \tag{34}
\end{align}
\kappa_{n+1}=\kappa_n+\Delta\gamma \tag{35}
\end{align}
と更新します。
ばね拘束の解き方を詳しく見る
ここからは、ばね拘束における全ひずみの求め方をもう少し詳しく説明します。まず、材料が弾性のままであれば、ばねとのつり合いから全ひずみを手で求められます。この弾性ばね解を基準にしておくと、その後に塑性が入ったとき、なぜ二分法のような反復計算が必要になるのかが見えやすくなります。
最初に、弾性ばね解の意味を確認します。ばね拘束を弾性範囲だけで考えると、つり合い式は手で解けます。応力を
\sigma=E\left(\varepsilon-\varepsilon_p^n-\varepsilon_{\mathrm{th}}^{n+1}\right) \tag{36}
\end{align}
とおいて式(27)へ代入すると、
EA\left(\varepsilon-\varepsilon_p^n-\varepsilon_{\mathrm{th}}^{n+1}\right)+k_sL_0\varepsilon=0 \tag{37}
\end{align}
です。これを \(\varepsilon\) について整理すると、
\varepsilon_{\mathrm{spring,elastic}}^{n+1}=\frac{EA}{EA+k_sL_0}\left(\varepsilon_p^n+\varepsilon_{\mathrm{th}}^{n+1}\right) \tag{38}
\end{align}
となります。さらに式(17)の \(r=k_sL_0/(EA)\) を使えば、
\varepsilon_{\mathrm{spring,elastic}}^{n+1}=\frac{1}{1+r}\left(\varepsilon_p^n+\varepsilon_{\mathrm{th}}^{n+1}\right) \tag{39}
\end{align}
です。\(r=0\) なら自由膨張と同じになり、\(r\) が大きくなるほど全ひずみはゼロに近づきます。つまり、ばね拘束は自由膨張と完全拘束の中間を連続的につなぐ条件です。
ただし、弾塑性解析では、式(39)だけで終わりではありません。この式は、材料が弾性のままである場合の閉じた解です。実際には、その全ひずみを使って式(18)以降の応力更新を行い、降伏が起これば塑性補正後の応力で改めて式(41)を満たす必要があります。そのため、塑性が入る場合は、次に示すように全ひずみを反復的に探します。
塑性が入る場合の二分法
弾性範囲では式(39)のように全ひずみを手で求められました。しかし塑性が入ると、応力は単純な弾性式だけでは決まりません。仮に全ひずみを \(\varepsilon\) と選び、その全ひずみに対する応力更新後の応力を \(\sigma_{\mathrm{update}}(\varepsilon)\) と書きます。このとき、棒の内力とばね反力のつり合い残差を
g(\varepsilon)=\sigma_{\mathrm{update}}(\varepsilon)A+k_sL_0\varepsilon \tag{40}
\end{align}
と定義します。ここで \(\sigma_{\mathrm{update}}(\varepsilon)\) は、式(18)から式(35)までの応力更新を、仮の全ひずみ \(\varepsilon\) に対して実行した結果です。
求めたい条件は
g(\varepsilon^{n+1})=0 \tag{41}
\end{align}
です。つまり、ばね拘束では、式(41)を満たす全ひずみを外側の反復計算で探し、その全ひずみに対して応力、塑性ひずみ、蓄積塑性ひずみを確定します。この手順を図3に示します。
図3の緑色の箱が、仮の全ひずみに対する応力更新です。この中身は、試行応力を計算する式(18)と、塑性補正を行う式(31)から式(35)に対応します。ここでは「仮の \(\varepsilon\) から \(\sigma_{\mathrm{update}}(\varepsilon)\) を返す処理」として見ておけば十分です。
この記事の計算では、全ひずみの探索に二分法を使います。二分法は、解をはさむ区間を少しずつ半分にしていく方法です。ニュートン法のように残差の微分を使わないため、ここで扱うような簡易モデルでは実装しやすく、安定に動かしやすいという利点があります。
まず残差の符号が変わる2つの全ひずみ \(\varepsilon_L\) と \(\varepsilon_R\) を用意します。
g(\varepsilon_L)g(\varepsilon_R)\le0 \tag{42}
\end{align}
この条件は、区間の左端と右端の間に \(g(\varepsilon)=0\) となる点があることを意味します。
次に、中点
\varepsilon_M=\frac{\varepsilon_L+\varepsilon_R}{2} \tag{43}
\end{align}
を計算し、\(g(\varepsilon_M)\) の符号を調べます。左端と中点で符号が変わる場合は、解は左半分にあります。そのため、
g(\varepsilon_L)g(\varepsilon_M)\le0\quad\Rightarrow\quad\varepsilon_R\leftarrow\varepsilon_M \tag{44}
\end{align}
とします。そうでなければ、解は右半分にあるので、
\varepsilon_L\leftarrow\varepsilon_M \tag{45}
\end{align}
とします。この操作を繰り返し、区間幅または残差が十分小さくなったところで
\varepsilon^{n+1}\simeq\frac{\varepsilon_L+\varepsilon_R}{2} \tag{46}
\end{align}
を採用します。この \(\varepsilon^{n+1}\) が、ばね拘束で実際に許される全ひずみです。
次の小節で、弾性ばね解から出発して二分法が必要になる様子を具体的な数値で確認します。
ばね拘束の小さな数値例
ここでは、図3の手順を、具体的な数値で一度だけなぞります。表1の同じ材料定数を使い、ある温度ステップで \(\Delta T_{n+1}=300 \mathrm{K}\)、\(\varepsilon_p^n=0\)、\(\kappa_n=0\)、\(r=1\) であったとします。
\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{初期降伏応力} \\
H & 10000\,\mathrm{MPa} & \text{等方線形硬化係数} \\
\Delta T_{\max} & 400\,\mathrm{K} & \text{最大温度上昇} \\
L_0 & 100\,\mathrm{mm} & \text{可視化用の棒長さ} \\
r & 1.0 & \text{ばね拘束の剛性比} \\
\hline
\end{array}
\end{array}
\end{align}
このとき熱ひずみは
\varepsilon_{\mathrm{th}}^{n+1}=\alpha\Delta T_{n+1}=12\times10^{-6}\times300=3.6\times10^{-3} \tag{47}
\end{align}
です。もし材料が弾性のままだと仮定すれば、式(39)から
\varepsilon_{\mathrm{spring,elastic}}^{n+1}=\frac{1}{1+1}\left(0+3.6\times10^{-3}\right)=1.8\times10^{-3} \tag{48}
\end{align}
となります。ここまでは、図3の「仮の全ひずみ」を弾性ばね解で与えたことに対応します。棒の自由膨張量の半分だけ実際に伸びる、という弾性ばね問題の答えです。しかし、この全ひずみを式(18)へ入れると、
\sigma^{\mathrm{trial}}=210000\left(1.8\times10^{-3}-3.6\times10^{-3}\right)=-378\,\mathrm{MPa} \tag{49}
\end{align}
となります。この時点、つまり塑性補正を行う前の降伏応力は \(250 \mathrm{MPa}\) なので、この仮定ではすでに圧縮降伏しています。
つまり、式(39)で得た全ひずみは、弾性問題としては正しくても、弾塑性問題の最終解ではありません。この全ひずみに対して式(31)の塑性補正を行うと、圧縮側なので
\Delta\gamma=\frac{\left|\sigma^{\mathrm{trial}}\right|-\sigma_{Y0}-H\kappa_n}{E+H}=\frac{378-250}{210000+10000}=5.82\times10^{-4} \tag{50}
\end{align}
となり、更新後の応力は
\sigma^{n+1}=-378+210000\times5.82\times10^{-4}=-255.8\,\mathrm{MPa} \tag{51}
\end{align}
となります。ここで、式(51)の応力の大きさが \(250 \mathrm{MPa}\) を少し超えている点に注意が必要です。これは計算ミスではありません。今回の材料モデルは線形硬化を含むため、塑性補正後の応力は初期降伏応力 \(\sigma_{Y0}\) ではなく、塑性変形によって上昇した現在の降伏応力に一致します。この数値例では
\sigma_Y^{n+1} =\sigma_{Y0}+ H \left(\kappa_n+\Delta\gamma\right)=250+10000\times5.82\times10^{-4}=255.8\,\mathrm{MPa} \tag{52}
\end{align}
となります。したがって、式(51)は、圧縮側の硬化後の降伏面\(\sigma^{n+1}=-\sigma_Y^{n+1}\)に戻った結果です。
ところが、ばね拘束では棒の応力だけでなく、ばね反力とのつり合いも満たす必要があります。式(27)を断面積 \(A\) で割ると、今回の \(r=1\) では
\sigma^{n+1}+rE\varepsilon^{n+1}=0 \tag{53}
\end{align}
と書けます。いま仮に用いた \(\varepsilon^{n+1}=1.8\times10^{-3}\) では
\sigma^{n+1}+rE\varepsilon^{n+1}=-255.8+210000\times1.8\times10^{-3}=122.2\,\mathrm{MPa}\neq0 \tag{54}
\end{align}
となり、つり合いが満たされていません。弾性と仮定して求めた伸びでは、塑性補正後の棒の応力とばね反力が釣り合わない、ということです。
そこで実際の数値計算では、仮の \(\varepsilon\) を少しずつ変えながら、式(41)を満たす値を探します。この数値例では、つり合いを満たす値はおよそ
\varepsilon^{n+1}=1.24\times10^{-3} \tag{55}
\end{align}
です。このとき
\sigma^{n+1}\simeq -261.0\,\mathrm{MPa},\qquad rE\varepsilon^{n+1}\simeq261.0\,\mathrm{MPa} \tag{56}
\end{align}
となり、棒の圧縮応力とばね反力が打ち消し合います。このように、ばね拘束では、全ひずみを決める問題と応力更新の問題が互いに結びついています。この結びつきが、完全拘束や自由膨張よりもばね拘束が一段だけ複雑になる理由です。
数値例の条件
ここでは、表 1 の材料定数を使い、これまでと同じ温度履歴として、温度を 0 K から 400 K まで上げ、その後 0 K まで戻します。材料モデルは等方線形硬化です。前節の小さな数値例では、ある 1ステップだけを取り出して、ばね拘束で全ひずみを探す考え方を確認しました。ここからは、同じ計算を温度履歴全体に適用し、自由膨張、完全拘束、ばね拘束の結果を比較します。
図4では、温度を上げてから戻す過程で、応力、伸び、塑性ひずみが同時にどのように変化するかを確認できます。静止画だけを見ると最終的な履歴線に目が行きますが、アニメーションでは「どの温度域で応力が増え、どの時点で塑性ひずみが進むか」を順番に追うことができます。
この後の各図では、同じ内容を応力、伸び、塑性ひずみに分けて読み直します。この数値例で特に見たいのは、材料定数や温度履歴を変えなくても、拘束条件だけで結果が大きく変わる点です。自由膨張、完全拘束、ばね拘束は、材料モデルとしては同じ応力更新式を使っています。
それにもかかわらず応力や塑性ひずみが変わるのは、応力更新に入る前の全ひずみ \(\varepsilon^{n+1}\) が異なるためです。つまり、ここから先の計算結果は、材料モデルの違いではなく、「どれだけ伸びることを許したか」の違いを見ていることになります。
計算結果：応力–温度線図
図4に、拘束条件ごとの応力–温度線図を示します。自由膨張では、棒が温度上昇に合わせて伸びられるため、応力はゼロのままです。完全拘束では、温度上昇とともに圧縮応力が増え、降伏後は塑性変形を伴いながら履歴を描きます。ばね拘束では、自由膨張と完全拘束の中間的な応答になります。
ここで重要なのは、温度履歴も材料定数も同じなのに、応力履歴が変わることです。つまり、熱応力は材料だけで決まる量ではありません。どれだけ自由膨張が妨げられるかが、応力の大きさを直接決めます。
完全拘束の線では、加熱初期に応力が直線的に圧縮側へ増えます。これは、温度上昇による自由熱ひずみ \(\alpha\Delta T\) が、そのまま弾性ひずみとして打ち消されるためです。しかし、応力が降伏応力に達すると、そこから先は温度を上げても応力だけが比例的に増え続けるわけではありません。塑性ひずみが進むことで、増えた熱ひずみの一部を塑性変形が受け持つようになるためです。
冷却時には、単純に加熱時の線を逆向きになぞるわけではありません。加熱中に生じた塑性ひずみが履歴として残っているため、同じ温度でも、加熱時と冷却時で弾性ひずみの内訳が変わります。このため、応力–温度線図には履歴が現れます。
熱弾塑性解析で残留応力が生じる本質は、この「温度は戻っても、塑性ひずみは完全には戻らない」という点にあります。ばね拘束の線は、自由膨張と完全拘束の中間にあります。これは、ばね拘束では自由熱膨張の一部が実際の伸びとして許され、残りだけが弾性ひずみ、つまり応力に変わるためです。
拘束剛性比 \(r\) を大きくすれば完全拘束に近づき、小さくすれば自由膨張に近づきます。したがって、ばね拘束の結果は、実構造物で周囲の剛性や治具の硬さが変わると残留応力も変わることを示す簡易モデルとして読めます。
計算結果：自由膨張したい長さと実際の伸び
図5は、自由膨張したい伸びと、実際に許された伸びを比較したものです。自由膨張では、実際の伸びが自由熱膨張量と一致します。完全拘束では、実際の伸びは常にゼロです。ばね拘束では、完全固定ほど拘束されず、自由膨張ほどは伸びない中間的な値になります。
この図を見ると、熱応力を直感的に理解できます。自由膨張したい伸びと実際の伸びの差が大きいほど、棒の内部には大きな弾性ひずみが生じます。その弾性ひずみが応力になります。
自由膨張の場合、黄色で示す自由膨張したい伸びと、実際の伸びが重なります。このとき、熱ひずみは棒全体の伸びとしてそのまま現れるため、弾性ひずみに変換されません。したがって、応力も塑性ひずみもほとんど発生しません。これは「温度が上がったから必ず熱応力が出る」のではなく、「伸びたいのに伸びられないときに熱応力が出る」ことを表しています。
完全拘束の場合、実際の伸びはゼロに固定されています。一方で、自由膨張したい伸びは温度とともに増えます。この差がすべて棒の内部変形として処理されるため、加熱初期には大きな圧縮応力が生じます。降伏後は、差の一部が塑性ひずみとして残るため、冷却して温度が戻っても内部状態は初期状態に戻りません。
ばね拘束の場合は、実際の伸びがゼロでも自由膨張量そのものでもなく、その間の値になります。これは、外部ばねが棒の伸びを完全には止めず、しかし自由にもさせないためです。
ばね拘束の図を見ると、熱応力を小さくするには「材料を柔らかくする」だけでなく、「構造として伸びを逃がす」ことも有効であることが分かります。
計算結果：塑性ひずみ
図６に、塑性ひずみと温度の関係を示します。自由膨張では応力が発生しないため、塑性ひずみも発生しません。完全拘束では、加熱中に圧縮側の塑性ひずみが発生し、冷却中には応力が反転して引張側へ戻ります。その結果、温度を戻しても塑性ひずみは履歴として残ります。
ばね拘束では、完全拘束より応力が小さくなるため、塑性ひずみの発生量も変わります。ただし、完全に自由ではないため、温度履歴によっては塑性変形が残ります。このように、残留応力や残留塑性ひずみを考えるときには、材料モデルだけでなく拘束条件を一緒に見る必要があります。
塑性ひずみの図で見るべき点は、塑性ひずみが「応力が降伏応力を超えた後にだけ」進むことです。自由膨張では応力がゼロなので、塑性ひずみは発生しません。完全拘束では、自由膨張が許されないため、比較的早い段階で降伏に到達し、加熱中に圧縮側の塑性ひずみが増えます。つまり、拘束が強いほど、同じ温度上昇でも塑性化しやすくなります。
冷却時に注目すると、温度は元に戻っていきますが、塑性ひずみは弾性ひずみのようには自動的に消えません。塑性ひずみは不可逆な履歴量として残るため、最終的に温度が初期値へ戻っても、棒の内部状態は初期状態とは異なります。この残った塑性ひずみと、拘束条件によって許される全ひずみとの整合を取るために、残留応力が現れます。
ばね拘束では、完全拘束よりも塑性ひずみの進み方が緩やかになります。これは、ばねが一部の伸びを許すことで、同じ温度上昇に対して発生する応力が小さくなるためです。ただし、ばね拘束でも応力が降伏応力に達すれば塑性ひずみは残ります。したがって、ばね拘束は「残留応力を完全になくす条件」ではなく、「拘束の強さに応じて残留応力と塑性ひずみを変える条件」と理解するのが自然です。
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 は下の設定を true にしたものだけ保存します。
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
%% 2. 出力設定
% ブログに貼るコードでは、まず静止画だけを作る設定にしておきます。
% MP4 は時間がかかるため、必要なものだけ true にしてください。
saveFigures = true;
saveResponseAnimation = false; % 温度・応力・流動応力・kappa の4分割動画
saveStressStrainAnimation = false; % SSカーブ上の応力経路動画
saveStrainAnimation = false; % ひずみ成分と温度の関係動画
saveReturnMappingAnimation = false; % Delta gamma の決まり方動画
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";
% 図と動画で共通して使う温度レベルと色です。
temperatureLevels = [0, 100, 200, 400];
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];
% 温度履歴: 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;
%% 3. 材料点解析
% 応力更新アルゴリズムを全温度ステップに適用します。
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));
%% 4. 静止画の作成
% Figure 1: 温度ごとのフローカーブ
if saveFigures
kappaAxis = linspace(0, 0.010, 600);
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);
end
%% 5. アニメーションの作成
% Animation 1: 温度履歴、応力、流動応力、蓄積塑性ひずみの全体像
if saveResponseAnimation
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 saveStressStrainAnimation
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 saveStrainAnimation
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 saveReturnMappingAnimation
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
%% 6. 解析用の関数
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 における単調負荷の応力-ひずみ線図を作ります。
% 横軸 epsAxis は機械ひずみ epsilon_m です。
% 各 epsilon_m に対して、
%
% 1. 弾性範囲かどうかを判定する
% 2. 塑性範囲なら kappa を求める
% 3. sigma = sign(epsilon_m) * sigma_flow(kappa, dT) とする
%
% という順番で、SSカーブ上の点を計算します。
sigmaCurve = zeros(size(epsAxis));
for i = 1:numel(epsAxis)
epsM = epsAxis(i);
s = sign(epsM);
absEpsM = abs(epsM);
% 塑性ひずみがまだない状態の流動応力です。
sigmaFlow0 = flow_stress(0, dT, sigmaY0, Q, b, beta, lowerBound);
if E * absEpsM <= sigmaFlow0
% 弾性範囲: epsilon_m がそのまま弾性ひずみです。
sigmaCurve(i) = E * epsM;
else
% 塑性範囲:
% 単調負荷では、求める kappa は塑性ひずみの大きさに対応します。
% E*(|epsilon_m|-kappa) = sigma_flow(kappa,dT) を解いています。
kappaPoint = solve_plastic_multiplier(E * absEpsM, 0, dT, E, sigmaY0, Q, b, beta, lowerBound);
sigmaCurve(i) = s * flow_stress(kappaPoint, dT, sigmaY0, Q, b, beta, lowerBound);
end
end
end
%% 7. 描画用の補助関数
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次元棒、同じ材料定数、同じ温度履歴に対して、拘束条件だけを変えたときの熱弾塑性応答を比較しました。自由膨張では、熱ひずみは発生しても、それが拘束されないため応力は生じません。完全拘束では、自由膨張がすべて妨げられるため、大きな熱応力が発生します。ばね拘束では、その中間として、外部ばねの剛性に応じた応力と伸びが決まります。
この比較から分かるように、熱応力は材料定数だけで決まるものではありません。材料がどれだけ伸びることを許されるか、つまり拘束条件が非常に重要です。完全拘束は熱応力を理解するための分かりやすい理想化ですが、実際の構造物では、周囲の部材や治具の剛性によって拘束の強さが変わります。
したがって、残留応力や塑性ひずみを考えるときには、材料モデルと同じくらい、拘束条件の設定にも注意する必要があります。
コメント