はじめに
この記事では、「両端固定された1 次元棒の加熱・冷却問題」を題材に、熱弾塑性解析の基本的な定式化を整理します。先に理論を立て、その理論を確かめるための最小限の数値実験として、後半でMATLAB による計算例を示します。狙いは、有限要素法の一般論へ進む前に、材料点レベルで
- 熱ひずみがどのように応力を生むのか
- 弾性応答と塑性応答をどのように判定するのか
- 塑性ひずみをどのように更新するのか
- 加熱・冷却後に残留応力や残留塑性ひずみがなぜ残るのか
を見通せるようにすることです。
モデル化する問題
長さ \(L_0\) の棒を考えます。棒は両端で完全に固定されており、軸方向には伸び縮みできないものとします。温度変化 \(\Delta T\) を与えると、材料は本来であれば熱膨張します。しかし、両端固定により全体の伸びが拘束されるため、棒の中には熱応力が発生します。
今回の簡易モデルでは以下を仮定します。
- 1次元の棒として扱う。
- 温度は棒全体で一様である。
- 全ひずみは両端固定により常にゼロである。
- 材料は線形弾性と完全塑性を組み合わせたモデルで表す。
- 降伏応力は温度に依存しない。
- 硬化は考えない。
この仮定は現実の熱弾塑性解析としては非常に単純です。しかし、熱ひずみ、拘束、降伏、塑性ひずみ、残留応力という重要な要素がすべて入っています。
り実際には長さが変わらない棒を表します。熱膨張差は見やすさのため拡大表示しています。
ひずみの分解
熱弾塑性解析では、全ひずみ \(\varepsilon\) を弾性ひずみ \(\varepsilon_e\)、塑性ひずみ \(\varepsilon_p\)、熱ひずみ \(\varepsilon_{\mathrm{th}}\) に分けて考えます。
\begin{align}
\varepsilon
=\varepsilon_e+\varepsilon_p+\varepsilon_{\mathrm{th}} \tag{1}
\end{align}
\varepsilon
=\varepsilon_e+\varepsilon_p+\varepsilon_{\mathrm{th}} \tag{1}
\end{align}
熱ひずみは、線膨張係数 \(\alpha\) と温度変化 \(\Delta T\) を用いて
\begin{align}
\varepsilon_{\mathrm{th}} = \alpha \Delta T \tag{2}
\end{align}
\varepsilon_{\mathrm{th}} = \alpha \Delta T \tag{2}
\end{align}
と書けます。今回の棒は両端固定なので、全ひずみは
\begin{align}
\varepsilon = 0 \tag{3}
\end{align}
\varepsilon = 0 \tag{3}
\end{align}
です。したがって弾性ひずみは
\begin{align}
\varepsilon_e
=- \varepsilon_p – \varepsilon_{\mathrm{th}} \tag{4}
\end{align}
\varepsilon_e
=- \varepsilon_p – \varepsilon_{\mathrm{th}} \tag{4}
\end{align}
となります。温度が上がると \(\varepsilon_{\mathrm{th}}\) は正になるため、塑性ひずみがまだゼロであれば弾性ひずみは負になります。つまり、棒には圧縮応力が生じます。
応力と降伏条件
弾性範囲では、応力 \(\sigma\) はフックの法則で与えられます。
\begin{align}
\sigma = E \varepsilon_e \tag{5}
\end{align}
\sigma = E \varepsilon_e \tag{5}
\end{align}
ここで \(E\) はヤング率です。完全塑性モデルでは、降伏条件を
\begin{align}
f(\sigma) = |\sigma| – \sigma_Y \le 0 \tag{6}
\end{align}
f(\sigma) = |\sigma| – \sigma_Y \le 0 \tag{6}
\end{align}
とします。\(\sigma_Y\) は降伏応力です。
- \(|\sigma| < \sigma_Y\) なら弾性状態です
- \(|\sigma| = \sigma_Y\) なら降伏状態です。
- 試行応力が \(|\sigma| > \sigma_Y\) になった場合、塑性ひずみを更新して応力を降伏面上へ戻します。
応力–ひずみ線図で見る降伏応力
降伏応力を視覚的に理解するには、応力–ひずみ線図、いわゆるSSカーブを見ると分かりやすいです。完全塑性モデルでは、弾性域では応力がひずみに比例し、降伏応力に達すると応力が頭打ちになります。
図では、\(+\sigma_Y\) と \(-\sigma_Y\) の破線が降伏応力です。応力がこの線に達すると、完全塑性モデルではそれ以上応力が増えず、代わりに塑性ひずみが更新されます。
この図は材料モデルそのものを表しており、温度履歴を与えたときにこの曲線上でどのような経路をたどるかは、後の数値実験で確認します。
応力更新アルゴリズム
数値解析では、温度履歴を小さなステップに分け、各ステップで応力と塑性ひずみを更新します。現在の塑性ひずみ \(\varepsilon_p\) が既知であるとして、まず弾性で進んだと仮定した試行弾性ひずみを計算します。
\begin{align}
\varepsilon_e^{\mathrm{trial}}
=\varepsilon-\varepsilon_p-\varepsilon_{\mathrm{th}} \tag{7}
\end{align}
\varepsilon_e^{\mathrm{trial}}
=\varepsilon-\varepsilon_p-\varepsilon_{\mathrm{th}} \tag{7}
\end{align}
今回の問題では \(\varepsilon = 0\) なので
\begin{align}
\varepsilon_e^{\mathrm{trial}}
=-\varepsilon_p-\alpha \Delta T \tag{8}
\end{align}
\varepsilon_e^{\mathrm{trial}}
=-\varepsilon_p-\alpha \Delta T \tag{8}
\end{align}
です。試行応力は
\begin{align}
\sigma^{\mathrm{trial}}
=E \varepsilon_e^{\mathrm{trial}} \tag{9}
\end{align}
\sigma^{\mathrm{trial}}
=E \varepsilon_e^{\mathrm{trial}} \tag{9}
\end{align}
で求めます。
弾性の場合
試行応力が降伏応力以内であれば、塑性ひずみは更新しません。
\begin{align}
|\sigma^{\mathrm{trial}}| \le \sigma_Y \tag{10}
\end{align}
|\sigma^{\mathrm{trial}}| \le \sigma_Y \tag{10}
\end{align}
なら、
\begin{align}
\sigma &= \sigma^{\mathrm{trial}} \tag{11}\\
\varepsilon_p^{\mathrm{new}} &= \varepsilon_p \tag{12}
\end{align}
\sigma &= \sigma^{\mathrm{trial}} \tag{11}\\
\varepsilon_p^{\mathrm{new}} &= \varepsilon_p \tag{12}
\end{align}
とします。
塑性の場合
試行応力が降伏応力を超えた場合、
\begin{align}
|\sigma^{\mathrm{trial}}| > \sigma_Y \tag{13}
\end{align}
|\sigma^{\mathrm{trial}}| > \sigma_Y \tag{13}
\end{align}
なので、完全塑性の仮定により応力を降伏応力に制限します。
\begin{align}
\sigma=\mathrm{sign}(\sigma^{\mathrm{trial}})\sigma_Y \tag{14}
\end{align}
\sigma=\mathrm{sign}(\sigma^{\mathrm{trial}})\sigma_Y \tag{14}
\end{align}
このとき弾性ひずみは
\begin{align}
\varepsilon_e = \frac{\sigma}{E} \tag{15}
\end{align}
\varepsilon_e = \frac{\sigma}{E} \tag{15}
\end{align}
です。ひずみ分解式
\begin{align}
\varepsilon = \varepsilon_e + \varepsilon_p + \varepsilon_{\mathrm{th}} \tag{16}
\end{align}
\varepsilon = \varepsilon_e + \varepsilon_p + \varepsilon_{\mathrm{th}} \tag{16}
\end{align}
を満たすように塑性ひずみを更新します。
\begin{align}
\varepsilon_p^{\mathrm{new}}
=\varepsilon-\varepsilon_{\mathrm{th}}-\varepsilon_e \tag{17}
\end{align}
\varepsilon_p^{\mathrm{new}}
=\varepsilon-\varepsilon_{\mathrm{th}}-\varepsilon_e \tag{17}
\end{align}
今回の両端固定棒では \(\varepsilon = 0\) なので、
\begin{align}
\varepsilon_p^{\mathrm{new}}
=-\alpha \Delta T-\frac{\sigma}{E} \tag{18}
\end{align}
\varepsilon_p^{\mathrm{new}}
=-\alpha \Delta T-\frac{\sigma}{E} \tag{18}
\end{align}
となります。
数値例の設定
ここまでの定式化がどのような応答を生むかを見るため、代表的な鋼材を意識した数値例を設定します。材料定数と温度履歴は表の通りです。
\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_Y & 250\,\mathrm{MPa} & \text{降伏応力} \\
\Delta T_{\max} & 400\,\mathrm{K} & \text{最大温度上昇} \\
L_0 & 100\,\mathrm{mm} & \text{可視化用の棒長さ} \\
\hline
\end{array}
\end{array}
\end{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_Y & 250\,\mathrm{MPa} & \text{降伏応力} \\
\Delta T_{\max} & 400\,\mathrm{K} & \text{最大温度上昇} \\
L_0 & 100\,\mathrm{mm} & \text{可視化用の棒長さ} \\
\hline
\end{array}
\end{array}
\end{align}
塑性ひずみがゼロの初期状態で、降伏が始まる温度上昇は
\begin{align}
E\alpha \Delta T_Y = \sigma_Y \tag{19}
\end{align}
E\alpha \Delta T_Y = \sigma_Y \tag{19}
\end{align}
より
\begin{align}
\Delta T_Y
=
\frac{\sigma_Y}{E\alpha}
=
\frac{250}{210000 \times 12\times10^{-6}}
\simeq 99.2 \ \mathrm{K} \tag{20}
\end{align}
\Delta T_Y
=
\frac{\sigma_Y}{E\alpha}
=
\frac{250}{210000 \times 12\times10^{-6}}
\simeq 99.2 \ \mathrm{K} \tag{20}
\end{align}
です。つまり、約 \(99\) K までは弾性圧縮、それ以上加熱すると圧縮降伏が始まります。
数値実験による確認
前節までの式に沿って、温度を \(0 \rightarrow 400\) K まで上げ、その後 \(0\) K まで戻す数値実験を行います。実装では、温度ステップごとに上の応力更新を適用し、応力、塑性ひずみ、熱ひずみを履歴として保存します。代表的な出力は以下の通りです。
Yield starts around dT = 99.21 K if eps_p = 0.
Final stress = 250.000 MPa
Final plastic strain = -1.190476e-03
Max thermal strain = 4.800000e-03 図５では、横軸を
\begin{align}
\varepsilon_{\mathrm{mech}}
=\varepsilon – \varepsilon_{\mathrm{th}} \tag{21}
\end{align}
\varepsilon_{\mathrm{mech}}
=\varepsilon – \varepsilon_{\mathrm{th}} \tag{21}
\end{align}
としました。これは、熱ひずみを除いた機械的なひずみに相当します。今回の両端固定棒では \(\varepsilon=0\) なので、
\begin{align}
\varepsilon_{\mathrm{mech}}
=-\varepsilon_{\mathrm{th}} \tag{22}
\end{align}
\varepsilon_{\mathrm{mech}}
=-\varepsilon_{\mathrm{th}} \tag{22}
\end{align}
です。加熱すると、棒は自由膨張したいが固定されているため、数値実験で得られる経路は機械的には圧縮側へ進みます。
加熱初期: 弾性圧縮
加熱初期では、熱ひずみ \(\varepsilon_{\mathrm{th}} = \alpha\Delta T\) が増えます。しかし棒は両端固定されているため、全ひずみはゼロです。この差を打ち消すために弾性ひずみは負になり、応力は圧縮側へ進みます。塑性ひずみがまだゼロなら、
\begin{align}
\sigma = -E\alpha\Delta T \tag{23}
\end{align}
\sigma = -E\alpha\Delta T \tag{23}
\end{align}
であり、応力と温度は直線関係になります。
降伏後: 応力はほぼ一定、塑性ひずみが増える
\(\Delta T \simeq 99.2\) K を超えると、圧縮応力が \(-\sigma_Y\) に達します。完全塑性モデルでは、それ以上応力の絶対値は増えません。代わりに、追加の熱膨張分は塑性ひずみとして蓄積されます。このため、応力履歴を見ると圧縮側の降伏応力で頭打ちになり、塑性ひずみ履歴を見ると加熱中に値が変化します。
最大温度での状態
最大温度上昇は \(400\) K なので、最大熱ひずみは
\begin{align}
\varepsilon_{\mathrm{th,max}}=
12\times10^{-6} \times 400=4.8\times10^{-3} \tag{24}
\end{align}
\varepsilon_{\mathrm{th,max}}=
12\times10^{-6} \times 400=4.8\times10^{-3} \tag{24}
\end{align}
です。可視化用の棒長さ \(L_0=100\) mm に対して、自由膨張できた場合の伸びは
\begin{align}
\Delta L_{\mathrm{free}}
=L_0 \varepsilon_{\mathrm{th,max}}
=100 \times 4.8\times10^{-3}
=0.48 \ \mathrm{mm} \tag{25}
\end{align}
\Delta L_{\mathrm{free}}
=L_0 \varepsilon_{\mathrm{th,max}}
=100 \times 4.8\times10^{-3}
=0.48 \ \mathrm{mm} \tag{25}
\end{align}
です。実際の伸びは固定によりゼロなので、この \(0.48\) mm 分の伸びたい効果が応力と塑性ひずみに変換されていると考えると直感的です。
冷却中: 応力の反転
冷却すると熱ひずみ \(\varepsilon_{\mathrm{th}}\) は小さくなります。しかし、加熱中に生じた塑性ひずみ \(\varepsilon_p\) は残っています。そのため、単純に元の経路を戻るのではなく、弾性ひずみの符号が変わり、応力は圧縮側から引張側へ移ります。
今回の完全塑性モデルでは、最終的に応力が \(+250\) MPa となっています。これは、冷却過程で引張側の降伏応力に到達したことを意味します。
冷却後に残るもの
最終温度では \(\Delta T = 0\) なので、熱ひずみはゼロです。それでも塑性ひずみは
\begin{align}
\varepsilon_p = -1.190476\times10^{-3} \tag{26}
\end{align}
\varepsilon_p = -1.190476\times10^{-3} \tag{26}
\end{align}
だけ残っています。両端固定により全ひずみはゼロなので、
\begin{align}
0 = \varepsilon_e + \varepsilon_p \tag{27}
\end{align}
0 = \varepsilon_e + \varepsilon_p \tag{27}
\end{align}
となり、
\begin{align}
\varepsilon_e = -\varepsilon_p \tag{28}
\end{align}
\varepsilon_e = -\varepsilon_p \tag{28}
\end{align}
です。したがって
\begin{align}
\sigma=E(-\varepsilon_p)=210000 \times 1.190476\times10^{-3}
=250 \ \mathrm{MPa} \tag{29}
\end{align}
\sigma=E(-\varepsilon_p)=210000 \times 1.190476\times10^{-3}
=250 \ \mathrm{MPa} \tag{29}
\end{align}
となります。これが冷却後に残る残留応力です。
なぜ弾性ひずみと塑性ひずみが残るのか
ここで注意したいのは、冷却後に弾性ひずみが残ることです。
ただし、これは弾性ひずみが永久変形として残るという意味ではありません。最後まで両端固定されているため、棒の全ひずみが
\begin{align}
\varepsilon = 0 \tag{29}
\end{align}
\varepsilon = 0 \tag{29}
\end{align}
に拘束されたままです。
一方、冷却後は \(\Delta T=0\) なので、
\begin{align}
\varepsilon_{\mathrm{th}} = 0 \tag{30}
\end{align}
\varepsilon_{\mathrm{th}} = 0 \tag{30}
\end{align}
です。しかし、加熱中に発生した塑性ひずみ \(\varepsilon_p\) は残ります。そのため、ひずみ分解式
\begin{align}
\varepsilon
=\varepsilon_e+\varepsilon_p+\varepsilon_{\mathrm{th}} \tag{31}
\end{align}
\varepsilon
=\varepsilon_e+\varepsilon_p+\varepsilon_{\mathrm{th}} \tag{31}
\end{align}
は、冷却後には
\begin{align}
0 = \varepsilon_e + \varepsilon_p \tag{32}
\end{align}
0 = \varepsilon_e + \varepsilon_p \tag{32}
\end{align}
となります。つまり、残った塑性ひずみを打ち消すように弾性ひずみが残り、その弾性ひずみによって残留応力が生じます。もし最後に固定を外せば、この弾性ひずみは解放され、応力も消えます。そのときに残るのが塑性ひずみ、つまり永久ひずみです。
では、なぜ温度を上げて戻しただけなのに塑性ひずみが残るのでしょうか。理由は、塑性ひずみが現在の温度だけで決まる量ではなく、途中で降伏した履歴によって決まる量だからです。加熱中、棒は自由膨張したいのに両端固定で伸びられないため、圧縮応力が増えます。
そして \(\Delta T \simeq 99.2\) K を超えると、圧縮側で降伏します。降伏すると、材料はそれ以上の熱膨張差を弾性ひずみだけでは受け止められなくなり、余った分を塑性ひずみとして更新します。塑性ひずみはばねの伸び縮みのように自動的には戻りません。
したがって、冷却によって熱ひずみがゼロに戻っても、途中で降伏して作られた塑性ひずみは残ります。この「途中で一度降伏した」という履歴が、最終状態に影響します。
この簡易モデルから分かること
このモデルで特に大事なのは、熱応力が単に「温度が上がったから出る」のではなく、**熱膨張したい変形が拘束されるために出る**という点です。
また、塑性が入ると、温度を元に戻しても状態は元に戻りません。これは塑性ひずみが履歴変数だからです。熱弾塑性解析では、現在の温度と変位だけでなく、過去にどのような温度・応力経路を通ったかが結果に影響します。
MATLABコード
今回用いたコードを以下に記載します。
%% thermal_elastoplastic_bar_demo.m
% 両端固定された1次元棒の熱弾塑性応答を追う簡易デモです。
%
% ポイント:
% - 棒は両端固定なので、全ひずみ eps_total は常に 0 とします。
% - 温度が上がると、本来は alpha * dT だけ自由膨張したくなります。
% - しかし両端固定で伸びられないため、圧縮応力が発生します。
% - 応力が降伏応力を超えると、完全塑性モデルとして塑性ひずみを更新します。
% - 加熱後に冷却すると、応力履歴と塑性ひずみ履歴が残ります。
clear; clc; close all;
%% 材料定数と解析条件
E = 210e3; % ヤング率 [MPa] = [N/mm^2]
alpha = 12e-6; % 線膨張係数 [1/K]
sigmaY = 250; % 降伏応力 [MPa]
Tmax = 400; % 最大温度上昇 [K]
nUp = 160; % 加熱ステップ数
nDown = 160; % 冷却ステップ数
L0 = 100; % 初期棒長さ [mm]。可視化用です。
visualExpansionScale = 25; % 熱膨張差を見やすくするための表示倍率。計算結果には影響しません。
doAnimation = true; % アニメーション表示をする場合は true
animationStride = 4; % 何ステップごとに描画するか
saveMp4 = true; % アニメーションをMP4で保存する場合は true
mp4FileName = "thermal_elastoplastic_bar_animation.mp4";
ssPathMp4FileName = "thermal_elastoplastic_stress_strain_path.mp4";
videoFrameRate = 12; % MP4のフレームレート [fps]
savePdfFigures = true; % PDFやブログに貼る静止画を保存する場合は true
saveStressStrainMp4 = true; % SSカーブ上の履歴アニメーションを保存する場合は true
% 温度履歴: 0 -> Tmax -> 0
Tup = linspace(0, Tmax, nUp);
Tdown = linspace(Tmax, 0, nDown);
dT = [Tup, Tdown(2:end)];
nStep = numel(dT);
% アニメーション中に横軸が動くと伸縮が分かりにくいので、
% 自由膨張を最大表示したときの右端位置を先に決めておきます。
maxFreeExpansion = L0 * alpha * max(dT);
LfreeVisualMax = L0 + visualExpansionScale * maxFreeExpansion;
barSketchXLim = [-16, LfreeVisualMax + 16];
%% 履歴配列
sigmaHist = zeros(1, nStep); % 応力履歴 [MPa]
epsPHist = zeros(1, nStep); % 塑性ひずみ履歴 [-]
epsThHist = zeros(1, nStep); % 熱ひずみ履歴 [-]
epsMechHist = zeros(1, nStep); % 機械ひずみ履歴 [-]。eps - eps_th に相当します。
stateHist = strings(1, nStep); % Elastic / Plastic
% 両端固定なので、棒の全ひずみはゼロです。
epsTotal = 0;
% 塑性ひずみの初期値
epsP = 0;
%% 熱弾塑性の応力更新
for i = 1:nStep
% 熱ひずみ
epsTh = alpha * dT(i);
% trial elastic strain:
% eps_total = eps_e + eps_p + eps_th より
% eps_e_trial = eps_total - eps_p - eps_th
epsETrial = epsTotal - epsP - epsTh;
sigmaTrial = E * epsETrial;
% 完全塑性モデル:
% trial 応力が降伏応力以内なら弾性、超えたら降伏面に戻します。
if abs(sigmaTrial) <= sigmaY
sigma = sigmaTrial;
stateHist(i) = "Elastic";
else
sigma = sign(sigmaTrial) * sigmaY;
epsE = sigma / E;
% eps_total = eps_e + eps_p + eps_th を満たすように塑性ひずみを更新
epsP = epsTotal - epsTh - epsE;
stateHist(i) = "Plastic";
end
sigmaHist(i) = sigma;
epsPHist(i) = epsP;
epsThHist(i) = epsTh;
epsMechHist(i) = epsTotal - epsTh;
end
%% アニメーション: 温度上昇、応力、棒の「自由膨張したい長さ」を同時に見る
if doAnimation
figAnim = figure("Name", "Thermo-elasto-plastic bar animation", ...
"Color", "w", "Position", [80 80 1100 720]);
if saveMp4
videoWriter = VideoWriter(mp4FileName, "MPEG-4");
videoWriter.FrameRate = videoFrameRate;
open(videoWriter);
end
for i = 1:animationStride:nStep
clf(figAnim);
tiledlayout(figAnim, 3, 2, "TileSpacing", "compact", "Padding", "compact");
% 温度履歴
nexttile;
plot(1:nStep, dT, "Color", [0.25 0.25 0.25], "LineWidth", 1.2); hold on;
plot(i, dT(i), "o", "MarkerFaceColor", [0.85 0.25 0.15], ...
"MarkerEdgeColor", "none", "MarkerSize", 8);
grid on;
xlabel("Step");
ylabel("\DeltaT [K]");
title("Temperature history");
xlim([1 nStep]);
ylim([0 Tmax * 1.08]);
% 応力履歴
nexttile;
plot(1:nStep, sigmaHist, "Color", [0.1 0.3 0.7], "LineWidth", 1.4); hold on;
yline(sigmaY, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.0);
yline(-sigmaY, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.0);
plot(i, sigmaHist(i), "o", "MarkerFaceColor", [0.1 0.3 0.7], ...
"MarkerEdgeColor", "none", "MarkerSize", 8);
grid on;
xlabel("Step");
ylabel("\sigma [MPa]");
title("Stress response");
xlim([1 nStep]);
ylim(1.15 * [-sigmaY sigmaY]);
% 応力-温度線図
nexttile;
plot(dT, sigmaHist, "Color", [0.15 0.45 0.25], "LineWidth", 1.8); hold on;
yline(sigmaY, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.0);
yline(-sigmaY, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.0);
plot(dT(i), sigmaHist(i), "o", "MarkerFaceColor", [0.15 0.45 0.25], ...
"MarkerEdgeColor", "none", "MarkerSize", 8);
grid on;
xlabel("\DeltaT [K]");
ylabel("\sigma [MPa]");
title("Stress-temperature path");
xlim([0 Tmax * 1.05]);
ylim(1.15 * [-sigmaY sigmaY]);
% 塑性ひずみ履歴
nexttile;
plot(1:nStep, epsPHist * 1e3, "Color", [0.55 0.2 0.65], "LineWidth", 1.4); hold on;
plot(i, epsPHist(i) * 1e3, "o", "MarkerFaceColor", [0.55 0.2 0.65], ...
"MarkerEdgeColor", "none", "MarkerSize", 8);
grid on;
xlabel("Step");
ylabel("\epsilon_p [x10^{-3}]");
title("Plastic strain history");
xlim([1 nStep]);
% 棒の模式図:
% 実際の熱膨張差は小さいため、模式図では差分だけを拡大して表示します。
nexttile([1 2]);
freeExpansion = L0 * epsThHist(i);
LfreeVisual = L0 + visualExpansionScale * freeExpansion;
drawBarSketch(L0, LfreeVisual, freeExpansion, visualExpansionScale, ...
sigmaHist(i), sigmaY, stateHist(i), barSketchXLim);
title(sprintf("Constrained bar: step %d / %d, \\DeltaT = %.1f K, \\sigma = %.1f MPa, %s", ...
i, nStep, dT(i), sigmaHist(i), stateHist(i)));
drawnow;
if saveMp4
writeVideo(videoWriter, getframe(figAnim));
end
pause(0.02);
end
if saveMp4
close(videoWriter);
fprintf("Saved animation MP4: %s\n", mp4FileName);
end
end
%% 最終まとめ図
figSummary = figure("Name", "Thermo-elasto-plastic bar summary", ...
"Color", "w", "Position", [120 120 1100 760]);
tiledlayout(figSummary, 2, 2, "TileSpacing", "compact", "Padding", "compact");
nexttile;
plot(1:nStep, dT, "k", "LineWidth", 1.5);
grid on;
xlabel("Step");
ylabel("\DeltaT [K]");
title("Temperature history");
xlim([1 nStep]);
nexttile;
plot(1:nStep, sigmaHist, "Color", [0.1 0.3 0.7], "LineWidth", 1.8); hold on;
yline(sigmaY, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.0);
yline(-sigmaY, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.0);
grid on;
xlabel("Step");
ylabel("\sigma [MPa]");
title("Stress history");
xlim([1 nStep]);
ylim(1.15 * [-sigmaY sigmaY]);
legend("\sigma", "+\sigma_Y", "-\sigma_Y", "Location", "best");
nexttile;
plot(dT, sigmaHist, "Color", [0.15 0.45 0.25], "LineWidth", 1.8); hold on;
yline(sigmaY, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.0);
yline(-sigmaY, "--", "Color", [0.7 0.1 0.1], "LineWidth", 1.0);
grid on;
xlabel("\DeltaT [K]");
ylabel("\sigma [MPa]");
title("Stress-temperature path");
xlim([0 Tmax * 1.05]);
ylim(1.15 * [-sigmaY sigmaY]);
nexttile;
plot(1:nStep, epsThHist * 1e3, "Color", [0.85 0.25 0.15], "LineWidth", 1.6); hold on;
plot(1:nStep, epsPHist * 1e3, "Color", [0.55 0.2 0.65], "LineWidth", 1.6);
grid on;
xlabel("Step");
ylabel("Strain [x10^{-3}]");
title("Thermal strain and plastic strain");
xlim([1 nStep]);
legend("\epsilon_{th}", "\epsilon_p", "Location", "best");
if savePdfFigures
exportgraphics(figSummary, "thermal_elastoplastic_summary.png", "Resolution", 200);
fprintf("Saved summary figure: thermal_elastoplastic_summary.png\n");
end
%% PDF用の応力-ひずみ線図を書き出し
if savePdfFigures
figMaterialStressStrain = figure("Name", "Material stress-strain curve", ...
"Color", "w", "Position", [180 180 1100 760]);
epsY = sigmaY / E;
epsRefMaterial = linspace(-4 * epsY, 4 * epsY, 500);
sigmaMaterial = min(max(E * epsRefMaterial, -sigmaY), sigmaY);
plot(epsRefMaterial * 1e3, sigmaMaterial, ...
"Color", [0.15 0.25 0.55], "LineWidth", 2.4); hold on;
lineSigmaYPos = yline(sigmaY, "--", "+\sigma_Y", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
lineSigmaYNeg = yline(-sigmaY, "--", "-\sigma_Y", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
lineEpsYPos = xline(epsY * 1e3, ":", "\epsilon_Y", "Color", [0.25 0.25 0.25], "LineWidth", 1.0);
lineEpsYNeg = xline(-epsY * 1e3, ":", "-\epsilon_Y", "Color", [0.25 0.25 0.25], "LineWidth", 1.0);
set([lineSigmaYPos, lineSigmaYNeg, lineEpsYPos, lineEpsYNeg], "FontSize", 20);
grid on;
set(gca, "FontSize", 18, "LineWidth", 1.1);
xlabel("Strain \epsilon [x10^{-3}]", "FontSize", 22);
ylabel("Stress \sigma [MPa]", "FontSize", 22);
title("Elastic-perfectly plastic stress-strain curve", "FontSize", 24);
ylim(1.25 * [-sigmaY sigmaY]);
exportgraphics(figMaterialStressStrain, "thermal_elastoplastic_material_ss_curve.png", "Resolution", 200);
copyfile("thermal_elastoplastic_material_ss_curve.png", "stress_strain_curve_material.png");
fprintf("Saved material SS curve: thermal_elastoplastic_material_ss_curve.png\n");
figStressStrain = figure("Name", "Stress-strain curve and thermal path", ...
"Color", "w", "Position", [180 180 1200 820]);
epsLimit = max(abs(epsMechHist)) * 1.1;
epsRef = linspace(-epsLimit, epsLimit, 500);
sigmaElasticPerfectPlastic = min(max(E * epsRef, -sigmaY), sigmaY);
plot(epsRef * 1e3, sigmaElasticPerfectPlastic, ...
"Color", [0.55 0.55 0.55], "LineWidth", 2.2); hold on;
scatter(epsMechHist * 1e3, sigmaHist, 26, 1:nStep, "filled");
plot(epsMechHist * 1e3, sigmaHist, "Color", [0.05 0.25 0.55 0.35], "LineWidth", 1.2);
lineSigmaYPos = yline(sigmaY, "--", "+\sigma_Y", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
lineSigmaYNeg = yline(-sigmaY, "--", "-\sigma_Y", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
lineEpsYPos = xline(sigmaY / E * 1e3, ":", "\epsilon_Y", "Color", [0.25 0.25 0.25], "LineWidth", 1.0);
lineEpsYNeg = xline(-sigmaY / E * 1e3, ":", "-\epsilon_Y", "Color", [0.25 0.25 0.25], "LineWidth", 1.0);
set([lineSigmaYPos, lineSigmaYNeg, lineEpsYPos, lineEpsYNeg], "FontSize", 20);
plot(epsMechHist(1) * 1e3, sigmaHist(1), "ko", "MarkerFaceColor", "w", "MarkerSize", 7);
plot(epsMechHist(end) * 1e3, sigmaHist(end), "ks", "MarkerFaceColor", "k", "MarkerSize", 7);
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 view of yielding and thermal loading path", "FontSize", 24);
legend("Elastic-perfectly plastic material curve", ...
"Thermal loading/cooling path", "Path line", ...
"Location", "best", "FontSize", 16);
colorbar("Ticks", [1 nStep], "TickLabels", ["start", "end"], "FontSize", 16);
ylim(1.2 * [-sigmaY sigmaY]);
exportgraphics(figStressStrain, "thermal_elastoplastic_stress_strain.png", "Resolution", 200);
copyfile("thermal_elastoplastic_stress_strain.png", "stress_strain_curve_with_history.png");
fprintf("Saved stress-strain figure: thermal_elastoplastic_stress_strain.png\n");
end
%% SSカーブ上の履歴アニメーションを書き出し
if saveStressStrainMp4
figStressStrainAnim = figure("Name", "Stress-strain path animation", ...
"Color", "w", "Position", [200 200 920 640]);
ssVideoWriter = VideoWriter(ssPathMp4FileName, "MPEG-4");
ssVideoWriter.FrameRate = videoFrameRate;
open(ssVideoWriter);
epsLimit = max(abs(epsMechHist)) * 1.1;
epsRef = linspace(-epsLimit, epsLimit, 500);
sigmaElasticPerfectPlastic = min(max(E * epsRef, -sigmaY), sigmaY);
for i = 1:animationStride:nStep
clf(figStressStrainAnim);
hold on;
% 材料モデルとしての完全弾塑性SSカーブ
plot(epsRef * 1e3, sigmaElasticPerfectPlastic, ...
"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_Y", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
lineSigmaYNeg = yline(-sigmaY, "--", "-\sigma_Y", "Color", [0.7 0.1 0.1], "LineWidth", 1.2);
lineEpsYPos = xline(sigmaY / E * 1e3, ":", "\epsilon_Y", "Color", [0.25 0.25 0.25], "LineWidth", 1.0);
lineEpsYNeg = xline(-sigmaY / E * 1e3, ":", "-\epsilon_Y", "Color", [0.25 0.25 0.25], "LineWidth", 1.0);
set([lineSigmaYPos, lineSigmaYNeg, lineEpsYPos, lineEpsYNeg], "FontSize", 20);
grid on;
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)));
legend("Elastic-perfectly plastic material curve", ...
"Thermal loading/cooling path", "Location", "southoutside");
colorbar("Ticks", [1 nStep], "TickLabels", ["start", "end"]);
xlim(epsLimit * 1e3 * [-1 1]);
ylim(1.2 * [-sigmaY sigmaY]);
drawnow;
writeVideo(ssVideoWriter, getframe(figStressStrainAnim));
end
close(ssVideoWriter);
fprintf("Saved stress-strain path MP4: %s\n", ssPathMp4FileName);
end
%% PDF用の模式図を書き出し
if savePdfFigures
[~, maxTempIndex] = max(dT);
figSketch = figure("Name", "Constrained bar sketch", ...
"Color", "w", "Position", [160 160 1000 320]);
freeExpansion = L0 * epsThHist(maxTempIndex);
LfreeVisual = L0 + visualExpansionScale * freeExpansion;
drawBarSketch(L0, LfreeVisual, freeExpansion, visualExpansionScale, ...
sigmaHist(maxTempIndex), sigmaY, stateHist(maxTempIndex), barSketchXLim);
title(sprintf("Maximum temperature: \\DeltaT = %.1f K, actual free expansion = %.3f mm, shown x%d", ...
dT(maxTempIndex), freeExpansion, visualExpansionScale));
exportgraphics(figSketch, "thermal_elastoplastic_bar_sketch.png", "Resolution", 200);
fprintf("Saved bar sketch figure: thermal_elastoplastic_bar_sketch.png\n");
end
%% 結果の要約をコマンドウィンドウに表示
yieldTemperature = sigmaY / (E * alpha);
fprintf("\n=== Summary ===\n");
fprintf("Yield starts around dT = %.2f K if eps_p = 0.\n", yieldTemperature);
fprintf("Final stress = %.3f MPa\n", sigmaHist(end));
fprintf("Final plastic strain = %.6e\n", epsPHist(end));
fprintf("Max thermal strain = %.6e\n", max(epsThHist));
%% ローカル関数
function drawBarSketch(Lfixed, LfreeVisual, freeExpansion, visualScale, sigma, sigmaY, stateText, xLimits)
% 固定された棒と、自由膨張できた場合の長さを模式的に描きます。
% 実際の熱膨張差はとても小さいため、表示上は差分だけを visualScale 倍しています。
cla;
hold on;
axis equal;
axis off;
wallColor = [0.18 0.18 0.18];
barColor = stressToColor(sigma, sigmaY);
freeColor = [0.95 0.7 0.25];
x0 = 0;
x1 = Lfixed;
y = 0;
h = 7;
% 固定壁。左右の壁は常に同じ位置に描くことで、
% 「棒は固定されているが、本当は右へ伸びたい」という見え方にします。
rectangle("Position", [x0 - 4, y - 12, 3, 24], "FaceColor", wallColor, "EdgeColor", "none");
rectangle("Position", [x1 + 1, y - 12, 3, 24], "FaceColor", wallColor, "EdgeColor", "none");
% 自由膨張できた場合の長さ。実際の棒の背後に薄く表示します。
rectangle("Position", [x0, y - h / 2 - 9, LfreeVisual, h], ...
"FaceColor", freeColor, "EdgeColor", "none", "FaceAlpha", 0.35);
text(min(Lfixed * 0.5, LfreeVisual * 0.5), y - 15, "free thermal expansion", ...
"HorizontalAlignment", "center", "FontSize", 10, "Color", [0.45 0.28 0.05]);
% 拘束により実際には伸びられない差分を、両矢印で強調します。
if LfreeVisual > Lfixed
plot([Lfixed Lfixed], [y - 16 y - 7], "Color", [0.45 0.28 0.05], "LineWidth", 1.2);
plot([LfreeVisual LfreeVisual], [y - 16 y - 7], "Color", [0.45 0.28 0.05], "LineWidth", 1.2);
quiver(Lfixed, y - 11.5, LfreeVisual - Lfixed, 0, 0, ...
"Color", [0.45 0.28 0.05], "LineWidth", 1.6, "MaxHeadSize", 0.8);
quiver(LfreeVisual, y - 11.5, -(LfreeVisual - Lfixed), 0, 0, ...
"Color", [0.45 0.28 0.05], "LineWidth", 1.6, "MaxHeadSize", 0.8);
text((Lfixed + LfreeVisual) / 2, y - 21, ...
sprintf("actual free expansion = %.3f mm (shown x%d)", freeExpansion, visualScale), ...
"HorizontalAlignment", "center", "FontSize", 10, "Color", [0.45 0.28 0.05]);
end
% 実際の固定棒
rectangle("Position", [x0, y - h / 2, Lfixed, h], ...
"FaceColor", barColor, "EdgeColor", [0.1 0.1 0.1], "LineWidth", 1.2);
text(Lfixed / 2, y + 9, "actual length is fixed", ...
"HorizontalAlignment", "center", "FontSize", 11, "Color", [0.1 0.1 0.1]);
% 応力の向きを矢印で表示。圧縮は内向き、引張は外向きです。
if sigma < 0
quiver(-10, y, 8, 0, 0, "Color", [0.1 0.3 0.7], "LineWidth", 2, "MaxHeadSize", 1.2);
quiver(Lfixed + 10, y, -8, 0, 0, "Color", [0.1 0.3 0.7], "LineWidth", 2, "MaxHeadSize", 1.2);
stressLabel = "compression";
elseif sigma > 0
quiver(2, y, -8, 0, 0, "Color", [0.7 0.1 0.1], "LineWidth", 2, "MaxHeadSize", 1.2);
quiver(Lfixed - 2, y, 8, 0, 0, "Color", [0.7 0.1 0.1], "LineWidth", 2, "MaxHeadSize", 1.2);
stressLabel = "tension";
else
stressLabel = "stress free";
end
text(Lfixed / 2, y - 31, sprintf("%s, %s", stressLabel, stateText), ...
"HorizontalAlignment", "center", "FontSize", 11, "FontWeight", "bold");
xlim(xLimits);
ylim([-40, 28]);
end
function c = stressToColor(sigma, sigmaY)
% 応力状態を色に変換します。
% 圧縮側は青、引張側は赤、ゼロ付近は薄いグレーです。
r = min(abs(sigma) / sigmaY, 1);
neutral = [0.75 0.78 0.80];
tension = [0.82 0.20 0.18];
compression = [0.18 0.36 0.78];
if sigma >= 0
c = (1 - r) * neutral + r * tension;
else
c = (1 - r) * neutral + r * compression;
end
endまとめ
両端固定棒の熱弾塑性問題は非常に単純ですが、熱弾塑性解析の基本構造をよく表しています。
- 熱ひずみは \(\varepsilon_{\mathrm{th}}=\alpha\Delta T\) で表される。
- 両端固定では全ひずみがゼロなので、熱ひずみは弾性ひずみや塑性ひずみとして吸収される。
- 弾性範囲では \(\sigma=-E\alpha\Delta T\) となる。
- 降伏後は応力が降伏応力に制限され、塑性ひずみが更新される。
- 冷却しても塑性ひずみは残るため、残留応力が発生する。
この流れを押さえると、より一般的な有限要素法の熱弾塑性解析でも、各積分点で何を更新しているのかが見えやすくなります。
コメント