熱弾塑性解析の基礎②：線形硬化【MATLABコード付き】

はじめに

　前回は、両端固定された1次元棒を題材に、完全弾塑性材料の熱弾塑性応答を整理しました。完全塑性モデルでは、応力が降伏応力に達すると、それ以上は応力が増えず、塑性ひずみだけが更新されます。これは考え方を理解するには分かりやすい一方で、実際の金属材料では、塑性変形が進むにつれてさらに変形しにくくなることがあります。

　この資料では、その効果を最も単純に表すために、等方線形硬化を導入します。先に理論を整理し、その後でMATLABによる最小限の数値実験を使って、完全塑性との違いを確認します。狙いは、材料点レベルで

線形硬化とは何を表すのか
降伏条件が完全塑性とどう変わるのか
累積塑性ひずみ \(\kappa\) がなぜ必要なのか
加熱・冷却後の残留応力が完全塑性とどう違うのか

を見通せるようにすることです。

モデル化する問題

　対象は前回と同じく、両端固定された長さ \(L_0\) の1次元棒です。棒は軸方向に伸び縮みできないため、全ひずみは常に

\begin{align}
\varepsilon = 0 \tag{1}
\end{align}

です。この式は、棒全体として見た軸方向の伸びがゼロであることを表しています。

　材料そのものが伸びたくないという意味ではなく、熱膨張しようとする変形が外部拘束によって打ち消される、という条件です。温度変化 \(\Delta T\) によって材料は自由膨張しようとしますが、その変形が拘束されるため、棒の中に熱応力が発生します。

今回の仮定は次の通りです。

1次元棒として扱います。
温度は棒全体で一様です。
全ひずみは両端固定により常にゼロです。
材料は線形弾性と等方線形硬化塑性で表します。
ヤング率、線膨張係数、初期降伏応力、硬化係数は温度に依存しないとします。

ひずみの分解

　熱弾塑性解析では、全ひずみを弾性ひずみ、塑性ひずみ、熱ひずみに分けます。

\begin{align}
\varepsilon=\varepsilon_e+\varepsilon_p+\varepsilon_{\mathrm{th}} \tag{2}
\end{align}

この式は、観測される全ひずみを、応力を生む弾性ひずみ、戻らない塑性ひずみ、温度変化だけで生じる熱ひずみに分けて考えるという意味です。

　熱弾塑性解析では、どの成分が応力を生み、どの成分が履歴として残るのかを分けることが重要になります。熱ひずみは

\begin{align}
\varepsilon_{\mathrm{th}} = \alpha \Delta T \tag{3}
\end{align}

です。この式は、拘束がなければ温度上昇に比例して棒が自由に伸びる、という熱膨張の最も単純な表現です。この段階では応力や塑性変形は関係しておらず、材料が本来取りたい変形量だけを表しています。

　今回の棒は両端固定なので、式(1)より

\begin{align}
\varepsilon_e=-\varepsilon_p-\varepsilon_{\mathrm{th}} \tag{4}
\end{align}

となります。この式は、固定された棒の中で弾性ひずみが調整役になることを表しています。

　熱ひずみや塑性ひずみが変わるたびに、全ひずみをゼロに保つため、弾性ひずみが反対向きに決まります。したがって、温度が上がると熱ひずみが正になり、それを打ち消すように圧縮側の弾性ひずみと応力が生じます。

等方線形硬化モデル

　弾性範囲では、応力はフックの法則で与えられます。

\begin{align}
\sigma = E \varepsilon_e \tag{5}
\end{align}

この式で応力に直接効いているのは弾性ひずみだけです。熱ひずみや塑性ひずみそのものが直接応力になるのではなく、それらによって必要になった弾性ひずみが応力を生みます。

　完全塑性では、降伏応力は一定値 \(\sigma_Y\) でした。一方、等方線形硬化では、塑性変形が進むほど降伏応力が大きくなると考えます。

　物理的には、一度塑性変形した材料の内部で、転位などの微視的な変形機構がさらに動きにくくなり、次の塑性変形を進めるためにより大きな応力が必要になる、というイメージです。ここではその複雑な内部構造の変化を細かく追跡せず、塑性変形の蓄積量に比例して降伏応力が増える、という単純な近似で表します。

　また「等方」と呼ぶのは、引張側だけ、あるいは圧縮側だけが強くなるのではなく、応力空間で降伏できる範囲が正負両側に同じように広がると考えるためです。そのため、現在の降伏応力

\begin{align}
\sigma_{\mathrm{y}}^{\mathrm{current}}
=\sigma_{Y0}+H \kappa \tag{6}
\end{align}

と置きます。ここで \(\sigma_{Y0}\) は初期降伏応力、\(H\) は等方線形硬化係数、\(\kappa\) は累積塑性ひずみです。この式は、材料が塑性変形を経験するほど、次に降伏させるために必要な応力が大きくなることを表しています。\(\kappa\) は向きを持つ塑性ひずみではなく、これまでにどれだけ塑性変形したかを数える履歴量です。

　降伏条件は

\begin{align}
f(\sigma,\kappa)
=|\sigma|-\left(\sigma_{Y0}+H\kappa\right)\le 0 \tag{7}
\end{align}

となります。この式は、現在の応力の大きさが、現在の降伏応力の内側に収まっているかを判定する条件です。\(f<0\) なら弾性範囲、\(f=0\) なら降伏面上、数値計算で試行的に \(f>0\) となった場合は塑性補正が必要です。完全塑性との違いは、右辺の降伏応力が一定ではなく、\(\kappa\) によって広がる点です。

*図１等方線形硬化を持つ材料の単調応力–ひずみ線図。降伏後も硬化係数に応じて応力が増えます。*

応力更新アルゴリズム

　数値解析では、温度履歴を小さなステップに分け、各ステップで応力、塑性ひずみ、累積塑性ひずみを更新します。まず、現在の塑性ひずみを固定したまま、弾性で進んだと仮定した試行応力を計算します。

\begin{align}
\varepsilon_e^{\mathrm{trial}}=\varepsilon-\varepsilon_p-\varepsilon_{\mathrm{th}}\tag{8}
\end{align}

これは、「もし今回の温度ステップを弾性だけで進めたら、弾性ひずみはいくつになるか」を仮に計算したものです。まだ塑性ひずみは更新せず、まず弾性予測として応力状態を見積もります。

\begin{align}
\sigma^{\mathrm{trial}}=E \varepsilon_e^{\mathrm{trial}} \tag{9}
\end{align}

この試行応力は、弾性予測だけで得られる仮の応力です。この値が降伏面の内側ならそのまま採用し、外側なら降伏面へ戻す補正を行います。

　試行降伏関数は

\begin{align}
f^{\mathrm{trial}}=|\sigma^{\mathrm{trial}}|-\left(\sigma_{Y0}+H\kappa\right)\tag{10}
\end{align}

です。この式は、試行応力が現在の降伏面をどれだけ超えているかを測る量です。\(f^{\mathrm{trial}}\) が正であるほど、弾性予測が許容範囲から大きく外れていることになります。

弾性の場合

\begin{align}
f^{\mathrm{trial}} \le 0 \tag{11}
\end{align}

なら、降伏していないので塑性ひずみと累積塑性ひずみは更新しません。

\begin{align}
\sigma &= \sigma^{\mathrm{trial}} \\
\varepsilon_p^{\mathrm{new}} &= \varepsilon_p \\
\kappa^{\mathrm{new}} &= \kappa \\ \tag{12}
\end{align}

塑性の場合

　\(f^{\mathrm{trial}}>0\) の場合、弾性のまま進めると応力が現在の降伏面の外側に出てしまいます。実際の弾塑性モデルでは、そのような応力状態は許されないため、塑性ひずみを増やして応力を降伏面上へ戻します。

　まず決めたいのは、「どれだけ塑性変形を進めれば、応力が降伏面上に戻るか」です。この補正の大きさを表す量を、塑性乗数増分 \(\Delta\gamma\) と呼びます。1次元の線形硬化モデルでは、試行応力が降伏面を超えた量 \(f^{\mathrm{trial}}\) を、弾性剛性 \(E\) と硬化係数 \(H\) の和で割ることで、\(\Delta\gamma\) を求めます。

\begin{align}
\Delta\gamma=\frac{f^{\mathrm{trial}}}{E+H} \tag{13}
\end{align}

分子の \(f^{\mathrm{trial}}\) は、試行応力が降伏面をどれだけ超えたかを表します。分母の \(E+H\) は、その超過量を弾性応力の低下と降伏応力の増加で受け持つことを表しています。そのため、硬化係数 \(H\) が大きいほど、同じ超過量に対して必要な塑性ひずみ増分は小さくなります。

　次に、塑性ひずみをどちらの向きに進めるかを決めます。1次元では、試行応力の符号

\begin{align}
s = \mathrm{sign}(\sigma^{\mathrm{trial}}) \tag{14}
\end{align}

を使えば十分です。この \(s\) を用いて、塑性ひずみを

\begin{align}
\varepsilon_p^{\mathrm{new}}
=\varepsilon_p+\Delta\gamma s \tag{15}
\end{align}

と更新します。圧縮側で降伏していれば負方向へ、引張側で降伏していれば正方向へ、塑性ひずみが進みます。

　塑性ひずみの向きが決まったら、次に硬化の程度を表す履歴変数を更新します。等方線形硬化では、降伏応力の増加量を累積塑性ひずみ \(\kappa\) で管理します。\(\kappa\) は塑性変形の「向き」ではなく「量」を蓄積する変数なので、今回生じた塑性増分 \(\Delta\gamma\) をそのまま足します。

\begin{align}
\kappa^{\mathrm{new}}
=\kappa+\Delta\gamma \tag{16}
\end{align}

　塑性ひずみ \(\varepsilon_p\) は符号を持ちますが、\(\kappa\) は塑性変形量の累積なので、塑性変形が起きるたびに増加します。

　最後に応力を更新します。塑性ひずみが増えると、その分だけ弾性ひずみは小さくなります。応力は弾性ひずみから決まるため、試行応力から \(E\Delta\gamma s\) を差し引きます。

\begin{align}
\sigma=\sigma^{\mathrm{trial}}-E\Delta\gamma s \tag{17}
\end{align}

　この更新後の応力と、更新後の \(\kappa\) から決まる降伏応力は、降伏条件を満たすように整合します。つまり塑性の場合の応力更新は、弾性予測で外へ出た応力を、塑性ひずみと硬化変数を更新しながら降伏面上へ戻す手続きだと考えればよいです。

数値例の設定

　材料定数と解析条件は次の通りです。

\begin{align}
\begin{array}{c}
\text{表: 数値例で用いる材料定数と解析条件。} \\
\begin{array}{ccl}
\hline
\text{記号} & \text{値} & \text{意味} \\
\hline
E & 210000\,\mathrm{MPa} & \text{ヤング率} \\
\alpha & 12\times10^{-6}\,\mathrm{K}^{-1} & \text{線膨張係数} \\
\sigma_{Y0} & 250\,\mathrm{MPa} & \text{初期降伏応力} \\
H & 10000\,\mathrm{MPa} & \text{等方線形硬化係数} \\
\Delta T_{\max} & 400\,\mathrm{K} & \text{最大温度上昇} \\
\hline
\end{array}
\end{array}
\end{align}

　塑性ひずみがまだゼロの初期状態では、降伏開始温度は完全塑性の場合と同じです。

\begin{align}
\Delta T_Y=\frac{\sigma_{Y0}}{E\alpha}=\frac{250}{210000\times12\times10^{-6}}\simeq 99.2 \ \mathrm{K} \tag{18}
\end{align}

数値実験による確認

　数値実験では、温度を \(0 \rightarrow 400\) K まで上げ、その後 \(0\) K まで戻します。各温度ステップで、上の応力更新を適用します。ここで確認したいのは、温度を一往復させたあとに、応力と履歴変数がどのような状態で残るかです。次の値はこの温度履歴を与えて計算した結果です。

Initial yield temperature = 99.21 K
Final stress = 304.413 MPa
Final plastic strain = -1.449587e-03
Final accumulated plastic strain = 5.441322e-03
Final yield stress = 304.413 MPa

　最初の `Initial yield temperature` は、初期状態で圧縮降伏が始まる温度上昇量です。両端固定棒では、温度上昇による自由膨張が拘束されるため、温度を上げるほど圧縮応力が増えます。その圧縮応力が初期降伏応力 \(250\) MPa に達する温度が、約 \(99.21\) K です。

`Final stress` は、温度を \(0\) K まで戻した最終ステップに残っている応力です。温度は元に戻っていますが、途中で塑性変形が発生しているため、応力はゼロには戻りません。この値は、熱サイクル後に残る残留応力を表しています。

`Final plastic strain` は、最終ステップに残った符号付きの塑性ひずみです。負の値になっているのは、加熱中に圧縮側の塑性変形が大きく進み、その一部が冷却後も残っているためです。

　一方、`Final accumulated plastic strain` は塑性変形量の累積値です。こちらは向きを持たないため、加熱時と冷却時に塑性変形が起きるたびに増えていきます。

　最後の `Final yield stress` は、線形硬化によって広がった最終的な降伏応力です。初期降伏応力は \(250\) MPa でしたが、塑性変形が蓄積したため、最終的には約 \(304\) MPa まで大きくなっています。この計算では最終応力と最終降伏応力がほぼ同じ値になっており、冷却後の最終状態が引張側の降伏面上にあることを示しています。

計算結果の解釈

　まず、計算結果の全体像を次の図に示します。

図３線形硬化モデルで得られた温度履歴、応力履歴、応力–温度線図、塑性ひずみと累積塑性ひずみの履歴。

　左上は温度履歴、右上は応力履歴と現在の降伏応力、左下は応力–温度線図、右下は塑性ひずみと累積塑性ひずみの履歴です。完全塑性モデルと比べて重要なのは、降伏応力が一定ではなく、塑性変形の蓄積に応じて赤い破線が広がっていく点です。

降伏後も応力が増える

　完全塑性モデルでは、圧縮降伏後の応力はほぼ \(-\sigma_Y\) に張り付きました。一方、線形硬化モデルでは、塑性ひずみが進むにつれて \(\kappa\) が増え、現在の降伏応力

\begin{align}
\sigma_{Y0}+H\kappa \tag{19}
\end{align}

が大きくなります。そのため、加熱中の圧縮応力は \(-250\) MPa で止まらず、さらに少し大きな圧縮応力へ進みます。

図４線形硬化材料の応力–ひずみ線図上に、熱負荷による応力経路を重ねた図。

ひずみ成分と温度変化の見方

　次に、熱ひずみ、塑性ひずみ、弾性ひずみを温度変化 \(\Delta T\) に対して整理します。次の図では、実線が加熱過程、破線が冷却過程です。熱ひずみは温度だけで決まるため、加熱時も冷却時も同じ直線上にあります。

　一方、塑性ひずみと弾性ひずみは履歴の影響を受けるため、同じ温度でも加熱時と冷却時で異なる値を取ります。

図５温度変化に対する熱ひずみ、塑性ひずみ、弾性ひずみの変化。実線は加熱、破線は冷却を表します。

累積塑性ひずみ \(\kappa\) は戻らない

　塑性ひずみ \(\varepsilon_p\) は符号を持つ量です。加熱中は圧縮側の塑性ひずみとして負方向に進み、冷却中に引張側へ反転すると値が戻る方向へ動きます。

　一方、累積塑性ひずみ \(\kappa\) は塑性変形量の蓄積です。そのため、塑性変形が起きるたびに増え、減りません。図３の右下を見ると、\(\varepsilon_p\) は途中で向きが変わりますが、\(\kappa\) は単調に増えていることが分かります。

冷却後の残留応力が大きくなる

　今回の計算では、冷却後の応力は約 \(304\) MPa です。これは初期降伏応力 \(250\) MPa より大きい値です。理由は、加熱・冷却の途中で塑性変形が蓄積し、降伏応力そのものが

\begin{align}
\sigma_{Y0}+H\kappa \tag{20}
\end{align}

まで広がったためです。最終的に引張側の降伏面に到達しているので、最終応力は現在の降伏応力とほぼ一致します。

完全塑性との違い

　完全塑性と線形硬化の違いをまとめると、次のようになります。

完全塑性では降伏応力が一定です。
線形硬化では塑性変形の蓄積に応じて降伏応力が大きくなります。
完全塑性では降伏後の応力が頭打ちになります。
線形硬化では降伏後も応力が少しずつ増えます。
線形硬化では内部変数 \(\kappa\) を保存する必要があります。

この簡易モデルから分かること

　線形硬化を入れると、熱弾塑性解析で保存すべき履歴変数が増えます。完全塑性では塑性ひずみ \(\varepsilon_p\) が主な履歴変数でしたが、等方線形硬化では累積塑性ひずみ \(\kappa\) も必要になります。この違いは小さく見えますが、有限要素法では重要です。なぜなら各積分点で、応力だけでなく、塑性ひずみ、累積塑性ひずみ、場合によっては硬化変数を保存し、次の荷重ステップへ引き継ぐ必要があるからです。

MATLABコード

　今回用いたコードを以下に記載します。

MATLABコード

%% thermal_elastoplastic_linear_hardening_demo.m
% 両端固定された1次元棒の熱弾塑性応答を、等方線形硬化モデルで計算します。
%
% 完全塑性モデルとの違い:
% - 降伏後も応力は一定値に張り付かず、塑性ひずみの蓄積に応じて少しずつ増えます。
% - 降伏応力は sigmaY + H * kappa として広がります。
% - kappa は累積塑性ひずみで、塑性変形の履歴を表す内部変数です。

clear; clc; close all;

%% 材料定数と解析条件
E = 210e3;             % ヤング率 [MPa]
alpha = 12e-6;         % 線膨張係数 [1/K]
sigmaY = 250;          % 初期降伏応力 [MPa]
H = 10e3;              % 等方線形硬化係数 [MPa]

Tmax = 400;            % 最大温度上昇 [K]
nUp = 160;             % 加熱ステップ数
nDown = 160;           % 冷却ステップ数
L0 = 100;              % 可視化用の棒長さ [mm]

saveFigures = true;
saveMp4 = true;
animationStride = 4;
videoFrameRate = 12;
responseMp4FileName = "linear_hardening_response_animation.mp4";
stressStrainMp4FileName = "linear_hardening_stress_strain_path.mp4";
strainTemperatureMp4FileName = "linear_hardening_strain_temperature_path.mp4";

% 温度履歴: 0 -> Tmax -> 0
Tup = linspace(0, Tmax, nUp);
Tdown = linspace(Tmax, 0, nDown);
dT = [Tup, Tdown(2:end)];
nStep = numel(dT);

%% 履歴配列
sigmaHist = zeros(1, nStep);     % 応力履歴 [MPa]
epsPHist = zeros(1, nStep);      % 塑性ひずみ履歴 [-]
epsThHist = zeros(1, nStep);     % 熱ひずみ履歴 [-]
epsEHist = zeros(1, nStep);      % 弾性ひずみ履歴 [-]
kappaHist = zeros(1, nStep);     % 累積塑性ひずみ履歴 [-]
yieldHist = zeros(1, nStep);     % 現在の降伏応力 [MPa]
epsMechHist = zeros(1, nStep);   % 機械ひずみ eps - eps_th [-]
stateHist = strings(1, nStep);   % Elastic / Plastic

epsTotal = 0;  % 両端固定なので全ひずみはゼロです。
epsP = 0;      % 塑性ひずみ
kappa = 0;     % 累積塑性ひずみ

%% 応力更新: 1次元の弾性予測・塑性補正
for i = 1:nStep
    epsTh = alpha * dT(i);
    epsETrial = epsTotal - epsP - epsTh;
    sigmaTrial = E * epsETrial;

    currentYield = sigmaY + H * kappa;
    fTrial = abs(sigmaTrial) - currentYield;

    if fTrial <= 0
        sigma = sigmaTrial;
        epsE = sigma / E;
        stateHist(i) = "Elastic";
    else
        direction = sign(sigmaTrial);
        dGamma = fTrial / (E + H);

        epsP = epsP + dGamma * direction;
        kappa = kappa + dGamma;

        sigma = sigmaTrial - E * dGamma * direction;
        epsE = sigma / E;
        stateHist(i) = "Plastic";
    end

    sigmaHist(i) = sigma;
    epsPHist(i) = epsP;
    epsThHist(i) = epsTh;
    epsEHist(i) = epsE;
    kappaHist(i) = kappa;
    yieldHist(i) = sigmaY + H * kappa;
    epsMechHist(i) = epsTotal - epsTh;
end

%% 材料の応力-ひずみ線図
if saveFigures
    epsY = sigmaY / E;
    epsRef = linspace(-5 * epsY, 5 * epsY, 700);
    sigmaMaterial = bilinearStress(epsRef, E, sigmaY, H);

    figMaterial = figure("Name", "Linear hardening material curve", ...
        "Color", "w", "Position", [160 140 1100 760]);
    plot(epsRef * 1e3, sigmaMaterial, "Color", [0.15 0.25 0.55], "LineWidth", 2.4); hold on;
    yline(sigmaY, "--", "+\sigma_{Y0}", "Color", [0.7 0.1 0.1], "LineWidth", 1.2, "FontSize", 18);
    yline(-sigmaY, "--", "-\sigma_{Y0}", "Color", [0.7 0.1 0.1], "LineWidth", 1.2, "FontSize", 18);
    xline(epsY * 1e3, ":", "\epsilon_{Y0}", "Color", [0.25 0.25 0.25], "LineWidth", 1.0, "FontSize", 18);
    xline(-epsY * 1e3, ":", "-\epsilon_{Y0}", "Color", [0.25 0.25 0.25], "LineWidth", 1.0, "FontSize", 18);
    grid on;
    set(gca, "FontSize", 18, "LineWidth", 1.1);
    xlabel("Strain  \epsilon  [x10^{-3}]", "FontSize", 22);
    ylabel("Stress  \sigma  [MPa]", "FontSize", 22);
    title("Elastic-plastic stress-strain curve with linear hardening", "FontSize", 23);
    exportgraphics(figMaterial, "linear_hardening_material_ss_curve.png", "Resolution", 200);
end

%% まとめ図
if saveFigures
    figSummary = figure("Name", "Linear hardening thermal elastoplastic summary", ...
        "Color", "w", "Position", [120 100 1200 880]);
    tiledlayout(figSummary, 2, 2, "TileSpacing", "compact", "Padding", "compact");

    nexttile;
    plot(1:nStep, dT, "k", "LineWidth", 1.6);
    grid on; xlim([1 nStep]);
    xlabel("Step", "FontSize", 16);
    ylabel("\DeltaT [K]", "FontSize", 16);
    title("Temperature history", "FontSize", 18);
    set(gca, "FontSize", 14, "LineWidth", 1.0);

    nexttile;
    plot(1:nStep, sigmaHist, "Color", [0.1 0.3 0.7], "LineWidth", 1.8); hold on;
    plot(1:nStep, yieldHist, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
    plot(1:nStep, -yieldHist, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
    grid on; xlim([1 nStep]);
    xlabel("Step", "FontSize", 16);
    ylabel("\sigma [MPa]", "FontSize", 16);
    title("Stress response and expanding yield stress", "FontSize", 18);
    legend("\sigma", "+(\sigma_{Y0}+H\kappa)", "-(\sigma_{Y0}+H\kappa)", ...
        "Location", "best", "FontSize", 12);
    set(gca, "FontSize", 14, "LineWidth", 1.0);

    nexttile;
    plot(dT, sigmaHist, "Color", [0.15 0.45 0.25], "LineWidth", 1.8); hold on;
    scatter(dT, sigmaHist, 20, 1:nStep, "filled");
    grid on; xlim([0 Tmax]);
    xlabel("\DeltaT [K]", "FontSize", 16);
    ylabel("\sigma [MPa]", "FontSize", 16);
    title("Stress-temperature path", "FontSize", 18);
    colorbar("Ticks", [1 nStep], "TickLabels", ["start", "end"], "FontSize", 12);
    set(gca, "FontSize", 14, "LineWidth", 1.0);

    nexttile;
    plot(1:nStep, epsPHist * 1e3, "Color", [0.55 0.2 0.65], "LineWidth", 1.8); hold on;
    plot(1:nStep, kappaHist * 1e3, "Color", [0.15 0.55 0.55], "LineWidth", 1.8);
    grid on; xlim([1 nStep]);
    xlabel("Step", "FontSize", 16);
    ylabel("Strain [x10^{-3}]", "FontSize", 16);
    title("Plastic strain and accumulated plastic strain", "FontSize", 18);
    legend("\epsilon_p", "\kappa", "Location", "best", "FontSize", 12);
    set(gca, "FontSize", 14, "LineWidth", 1.0);

    exportgraphics(figSummary, "linear_hardening_summary.png", "Resolution", 200);
end

%% ひずみ成分-温度線図
if saveFigures
    figStrainTemperature = figure("Name", "Linear hardening strain components vs temperature", ...
        "Color", "w", "Position", [160 100 1100 900]);
    tiledlayout(figStrainTemperature, 3, 1, "TileSpacing", "compact", "Padding", "compact");

    heatingIndex = 1:nUp;
    coolingIndex = nUp:nStep;
    strainScale = 1e3;

    nexttile;
    plot(dT(heatingIndex), epsThHist(heatingIndex) * strainScale, "-", ...
        "Color", [0.85 0.25 0.15], "LineWidth", 2.0); hold on;
    plot(dT(coolingIndex), epsThHist(coolingIndex) * strainScale, "--", ...
        "Color", [0.85 0.25 0.15], "LineWidth", 2.0);
    grid on; xlim([0 Tmax]);
    ylabel("\epsilon_{th} [x10^{-3}]", "FontSize", 16);
    title("Thermal strain vs temperature change", "FontSize", 18);
    legend("heating", "cooling", "Location", "northwest", "FontSize", 12);
    set(gca, "FontSize", 14, "LineWidth", 1.0);

    nexttile;
    plot(dT(heatingIndex), epsPHist(heatingIndex) * strainScale, "-", ...
        "Color", [0.55 0.2 0.65], "LineWidth", 2.0); hold on;
    plot(dT(coolingIndex), epsPHist(coolingIndex) * strainScale, "--", ...
        "Color", [0.55 0.2 0.65], "LineWidth", 2.0);
    yline(0, ":", "Color", [0.35 0.35 0.35], "LineWidth", 1.0);
    grid on; xlim([0 Tmax]);
    ylabel("\epsilon_p [x10^{-3}]", "FontSize", 16);
    title("Plastic strain vs temperature change", "FontSize", 18);
    set(gca, "FontSize", 14, "LineWidth", 1.0);

    nexttile;
    plot(dT(heatingIndex), epsEHist(heatingIndex) * strainScale, "-", ...
        "Color", [0.1 0.3 0.7], "LineWidth", 2.0); hold on;
    plot(dT(coolingIndex), epsEHist(coolingIndex) * strainScale, "--", ...
        "Color", [0.1 0.3 0.7], "LineWidth", 2.0);
    grid on; xlim([0 Tmax]);
    xlabel("\DeltaT [K]", "FontSize", 16);
    ylabel("\epsilon_e [x10^{-3}]", "FontSize", 16);
    title("Elastic strain vs temperature change", "FontSize", 18);
    set(gca, "FontSize", 14, "LineWidth", 1.0);

    exportgraphics(figStrainTemperature, "linear_hardening_strain_temperature.png", "Resolution", 200);
end

%% 応力-ひずみ線図上の履歴
if saveFigures
    figPath = figure("Name", "Linear hardening stress-strain path", ...
        "Color", "w", "Position", [160 120 1200 820]);
    epsLimit = max(abs(epsMechHist)) * 1.15;
    epsRef = linspace(-epsLimit, epsLimit, 800);
    sigmaMaterial = bilinearStress(epsRef, E, sigmaY, H);

    plot(epsRef * 1e3, sigmaMaterial, "Color", [0.55 0.55 0.55], "LineWidth", 2.2); hold on;
    plot(epsMechHist * 1e3, sigmaHist, "Color", [0.05 0.25 0.55 0.45], "LineWidth", 1.4);
    scatter(epsMechHist * 1e3, sigmaHist, 28, 1:nStep, "filled");
    plot(epsMechHist(1) * 1e3, sigmaHist(1), "ko", "MarkerFaceColor", "w", "MarkerSize", 8);
    plot(epsMechHist(end) * 1e3, sigmaHist(end), "ks", "MarkerFaceColor", "k", "MarkerSize", 8);
    grid on;
    set(gca, "FontSize", 18, "LineWidth", 1.1);
    xlabel("Mechanical strain  \epsilon - \epsilon_{th}  [x10^{-3}]", "FontSize", 22);
    ylabel("Stress  \sigma  [MPa]", "FontSize", 22);
    title("Stress-strain path with linear hardening", "FontSize", 24);
    legend("Monotonic bilinear material curve", "Thermal loading/cooling path", ...
        "Path points", "start", "end", "Location", "best", "FontSize", 14);
    colorbar("Ticks", [1 nStep], "TickLabels", ["start", "end"], "FontSize", 14);
    exportgraphics(figPath, "linear_hardening_stress_strain.png", "Resolution", 200);
end

%% ブログ用アニメーション: 全体応答
if saveMp4
    figResponseAnim = figure("Name", "Linear hardening response animation", ...
        "Color", "w", "Position", [120 100 1200 880]);
    responseVideoWriter = VideoWriter(responseMp4FileName, "MPEG-4");
    responseVideoWriter.FrameRate = videoFrameRate;
    open(responseVideoWriter);

    for i = 1:animationStride:nStep
        clf(figResponseAnim);
        tiledlayout(figResponseAnim, 2, 2, "TileSpacing", "compact", "Padding", "compact");

        nexttile;
        plot(1:nStep, dT, "k", "LineWidth", 1.6); hold on;
        plot(i, dT(i), "o", "MarkerFaceColor", "k", "MarkerEdgeColor", "w", "MarkerSize", 8);
        grid on; xlim([1 nStep]); ylim([0 Tmax * 1.08]);
        xlabel("Step"); ylabel("\DeltaT [K]");
        title(sprintf("Temperature history: step %d / %d", i, nStep));
        set(gca, "FontSize", 14, "LineWidth", 1.0);

        nexttile;
        plot(1:nStep, sigmaHist, "Color", [0.1 0.3 0.7], "LineWidth", 1.8); hold on;
        plot(1:nStep, yieldHist, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
        plot(1:nStep, -yieldHist, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
        plot(i, sigmaHist(i), "o", "MarkerFaceColor", [0.1 0.3 0.7], "MarkerEdgeColor", "w", "MarkerSize", 8);
        grid on; xlim([1 nStep]);
        xlabel("Step"); ylabel("\sigma [MPa]");
        title(sprintf("Stress response: \\sigma = %.1f MPa,  %s", sigmaHist(i), stateHist(i)));
        legend("\sigma", "+(\sigma_{Y0}+H\kappa)", "-(\sigma_{Y0}+H\kappa)", ...
            "Location", "best", "FontSize", 10);
        set(gca, "FontSize", 14, "LineWidth", 1.0);

        nexttile;
        plot(dT, sigmaHist, "Color", [0.15 0.45 0.25], "LineWidth", 1.8); hold on;
        plot(dT(i), sigmaHist(i), "o", "MarkerFaceColor", [0.15 0.45 0.25], ...
            "MarkerEdgeColor", "w", "MarkerSize", 8);
        grid on; xlim([0 Tmax]);
        xlabel("\DeltaT [K]"); ylabel("\sigma [MPa]");
        title("Stress-temperature path");
        set(gca, "FontSize", 14, "LineWidth", 1.0);

        nexttile;
        plot(1:nStep, epsPHist * 1e3, "Color", [0.55 0.2 0.65], "LineWidth", 1.8); hold on;
        plot(1:nStep, kappaHist * 1e3, "Color", [0.15 0.55 0.55], "LineWidth", 1.8);
        plot(i, epsPHist(i) * 1e3, "o", "MarkerFaceColor", [0.55 0.2 0.65], ...
            "MarkerEdgeColor", "w", "MarkerSize", 8);
        plot(i, kappaHist(i) * 1e3, "s", "MarkerFaceColor", [0.15 0.55 0.55], ...
            "MarkerEdgeColor", "w", "MarkerSize", 8);
        grid on; xlim([1 nStep]);
        xlabel("Step"); ylabel("Strain [x10^{-3}]");
        title("Plastic strain and accumulated plastic strain");
        legend("\epsilon_p", "\kappa", "current \epsilon_p", "current \kappa", ...
            "Location", "best", "FontSize", 10);
        set(gca, "FontSize", 14, "LineWidth", 1.0);

        drawnow;
        writeVideo(responseVideoWriter, getframe(figResponseAnim));
    end

    close(responseVideoWriter);
    fprintf("Saved response animation MP4: %s\n", responseMp4FileName);
end

%% ブログ用アニメーション: 応力-ひずみ線図上の履歴
if saveMp4
    figStressStrainAnim = figure("Name", "Linear hardening stress-strain animation", ...
        "Color", "w", "Position", [200 160 980 720]);
    stressStrainVideoWriter = VideoWriter(stressStrainMp4FileName, "MPEG-4");
    stressStrainVideoWriter.FrameRate = videoFrameRate;
    open(stressStrainVideoWriter);

    epsLimit = max(abs(epsMechHist)) * 1.1;
    epsRef = linspace(-epsLimit, epsLimit, 500);
    sigmaMaterial = bilinearStress(epsRef, E, sigmaY, H);
    sigmaLimit = max(abs([sigmaMaterial(:); sigmaHist(:)])) * 1.15;

    for i = 1:animationStride:nStep
        clf(figStressStrainAnim);
        axStressStrainAnim = axes(figStressStrainAnim);
        hold on;

        % 材料モデルとしての線形硬化SSカーブ
        plot(epsRef * 1e3, sigmaMaterial, ...
            "Color", [0.58 0.58 0.58], "LineWidth", 2.4);

        % 現在ステップまでに実際にたどった履歴
        plot(epsMechHist(1:i) * 1e3, sigmaHist(1:i), ...
            "Color", [0.05 0.25 0.55], "LineWidth", 2.0);
        scatter(epsMechHist(1:i) * 1e3, sigmaHist(1:i), ...
            22, 1:i, "filled");

        % 現在位置を大きめに表示
        plot(epsMechHist(i) * 1e3, sigmaHist(i), "o", ...
            "MarkerFaceColor", [0.9 0.25 0.1], ...
            "MarkerEdgeColor", "k", "MarkerSize", 9, "LineWidth", 1.0);

        lineSigmaYPos = yline(sigmaY, "--", "+\sigma_{Y0}", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
        lineSigmaYNeg = yline(-sigmaY, "--", "-\sigma_{Y0}", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
        lineEpsYPos = xline(sigmaY / E * 1e3, ":", "\epsilon_{Y0}", "Color", [0.25 0.25 0.25], "LineWidth", 1.0);
        lineEpsYNeg = xline(-sigmaY / E * 1e3, ":", "-\epsilon_{Y0}", "Color", [0.25 0.25 0.25], "LineWidth", 1.0);
        set([lineSigmaYPos, lineSigmaYNeg, lineEpsYPos, lineEpsYNeg], "FontSize", 20);

        grid on;
        box on;
        set(axStressStrainAnim, "FontSize", 18, "LineWidth", 1.1, ...
            "Layer", "top", "Position", [0.12 0.24 0.66 0.58]);
        xlabel("Mechanical strain  \epsilon - \epsilon_{th}  [x10^{-3}]");
        ylabel("Stress  \sigma  [MPa]");
        title(sprintf("Stress-strain path: step %d / %d,  \\DeltaT = %.1f K,  \\sigma = %.1f MPa,  %s", ...
            i, nStep, dT(i), sigmaHist(i), stateHist(i)), "FontSize", 17);
        legend("Linear hardening material curve", ...
            "Thermal loading/cooling path", "Location", "southoutside");
        clim(axStressStrainAnim, [1 nStep]);
        cbStressStrainAnim = colorbar("Ticks", [1 nStep], "TickLabels", ["start", "end"]);
        cbStressStrainAnim.FontSize = 16;
        cbStressStrainAnim.Position = [0.83 0.28 0.024 0.48];
        set(axStressStrainAnim, "Position", [0.12 0.24 0.66 0.58]);
        xlim(epsLimit * 1e3 * [-1 1]);
        ylim(sigmaLimit * [-1 1]);

        drawnow;
        writeVideo(stressStrainVideoWriter, getframe(figStressStrainAnim));
    end

    close(stressStrainVideoWriter);
    fprintf("Saved stress-strain path MP4: %s\n", stressStrainMp4FileName);
end

%% ブログ用アニメーション: ひずみ成分-温度線図
if saveMp4
    figStrainTemperatureAnim = figure("Name", "Linear hardening strain-temperature animation", ...
        "Color", "w", "Position", [160 100 1100 900]);
    strainTemperatureVideoWriter = VideoWriter(strainTemperatureMp4FileName, "MPEG-4");
    strainTemperatureVideoWriter.FrameRate = videoFrameRate;
    open(strainTemperatureVideoWriter);

    heatingIndex = 1:nUp;
    coolingIndex = nUp:nStep;
    strainScale = 1e3;

    for i = 1:animationStride:nStep
        clf(figStrainTemperatureAnim);
        tiledlayout(figStrainTemperatureAnim, 3, 1, "TileSpacing", "compact", "Padding", "compact");

        nexttile;
        plot(dT(heatingIndex), epsThHist(heatingIndex) * strainScale, "-", ...
            "Color", [0.85 0.25 0.15], "LineWidth", 2.0); hold on;
        plot(dT(coolingIndex), epsThHist(coolingIndex) * strainScale, "--", ...
            "Color", [0.85 0.25 0.15], "LineWidth", 2.0);
        plot(dT(i), epsThHist(i) * strainScale, "o", "MarkerFaceColor", "k", ...
            "MarkerEdgeColor", "w", "MarkerSize", 8);
        grid on; xlim([0 Tmax]);
        ylabel("\epsilon_{th} [x10^{-3}]");
        title(sprintf("Strain components vs temperature: step %d / %d,  \\DeltaT = %.1f K", ...
            i, nStep, dT(i)));
        set(gca, "FontSize", 14, "LineWidth", 1.0);

        nexttile;
        plot(dT(heatingIndex), epsPHist(heatingIndex) * strainScale, "-", ...
            "Color", [0.55 0.2 0.65], "LineWidth", 2.0); hold on;
        plot(dT(coolingIndex), epsPHist(coolingIndex) * strainScale, "--", ...
            "Color", [0.55 0.2 0.65], "LineWidth", 2.0);
        yline(0, ":", "Color", [0.35 0.35 0.35], "LineWidth", 1.0);
        plot(dT(i), epsPHist(i) * strainScale, "o", "MarkerFaceColor", "k", ...
            "MarkerEdgeColor", "w", "MarkerSize", 8);
        grid on; xlim([0 Tmax]);
        ylabel("\epsilon_p [x10^{-3}]");
        title("Plastic strain vs temperature change");
        set(gca, "FontSize", 14, "LineWidth", 1.0);

        nexttile;
        plot(dT(heatingIndex), epsEHist(heatingIndex) * strainScale, "-", ...
            "Color", [0.1 0.3 0.7], "LineWidth", 2.0); hold on;
        plot(dT(coolingIndex), epsEHist(coolingIndex) * strainScale, "--", ...
            "Color", [0.1 0.3 0.7], "LineWidth", 2.0);
        plot(dT(i), epsEHist(i) * strainScale, "o", "MarkerFaceColor", "k", ...
            "MarkerEdgeColor", "w", "MarkerSize", 8);
        grid on; xlim([0 Tmax]);
        xlabel("\DeltaT [K]");
        ylabel("\epsilon_e [x10^{-3}]");
        title("Elastic strain vs temperature change");
        set(gca, "FontSize", 14, "LineWidth", 1.0);

        drawnow;
        writeVideo(strainTemperatureVideoWriter, getframe(figStrainTemperatureAnim));
    end

    close(strainTemperatureVideoWriter);
    fprintf("Saved strain-temperature path MP4: %s\n", strainTemperatureMp4FileName);
end

fprintf("\n=== Linear hardening summary ===\n");
fprintf("Initial yield temperature = %.2f K\n", sigmaY / (E * alpha));
fprintf("Final stress              = %.3f MPa\n", sigmaHist(end));
fprintf("Final plastic strain      = %.6e\n", epsPHist(end));
fprintf("Final accumulated plastic strain = %.6e\n", kappaHist(end));
fprintf("Final yield stress        = %.3f MPa\n", yieldHist(end));

function sigma = bilinearStress(epsRef, E, sigmaY, H)
    epsY = sigmaY / E;
    sigma = zeros(size(epsRef));
    elasticIndex = abs(epsRef) <= epsY;
    plasticIndex = ~elasticIndex;

    sigma(elasticIndex) = E * epsRef(elasticIndex);
    sigma(plasticIndex) = sign(epsRef(plasticIndex)) .* ...
        (sigmaY + (E * H / (E + H)) .* (abs(epsRef(plasticIndex)) - epsY));
end

まとめ

　この記事では、完全塑性モデルに等方線形硬化を加えた1次元熱弾塑性解析を扱いました。完全塑性では、降伏後の応力は一定値に張り付きますが、線形硬化を入れると、塑性変形の蓄積に応じて材料がさらに変形しにくくなる様子を表せます。そのため、加熱中の応力履歴、冷却後の残留応力、SSカーブ上の経路は、完全塑性モデルよりも少し現実の材料応答に近い形になります。

線形硬化では、降伏応力を \(\sigma_{Y0}+H\kappa\) として表します。
\(\kappa\) は累積塑性ひずみであり、塑性変形の履歴を表します。
降伏後も応力は一定値に張り付かず、硬化係数に応じて増えます。
冷却後の残留応力は、現在の降伏応力と関係します。
熱弾塑性解析では、現在の温度だけでなく、過去の塑性履歴を保存することが重要です。

　特に重要なのは、応力を計算するだけでなく、塑性ひずみと累積塑性ひずみを次のステップへ引き継ぐという考え方です。熱弾塑性解析では、同じ温度であっても、そこに至るまでの加熱・冷却履歴によって応力やひずみが変わります。線形硬化モデルは、その「履歴を持つ材料応答」を理解するためのよい入口になります。