熱弾塑性解析の基礎⑥：1次元棒要素FEMと熱応力解析【MATLABコード付き】

はじめに

　これまでの記事では、1つの材料点、または棒全体が同じ温度を受ける問題として、熱弾塑性解析を見てきました。 その扱いでは、温度、応力、塑性ひずみは棒のどこでも同じ値として考えます。 この仮定は、材料モデルの意味を理解するにはとても便利です。 一方で、実際の構造物では、温度は場所によって変わります。 加熱されている部分、まだ冷たい部分、周囲から拘束されている部分が同時に存在します。

　温度や応力が場所によって変わる問題へ進むには、棒を空間方向に分けて考える必要があります。 そこでこの記事では、熱弾塑性へ進む前の土台として、1次元棒要素の有限要素法を整理します。 ここでは塑性はまだ入れず、線形弾性の範囲で、

• 棒の力の釣合いから支配方程式を定めること
• 支配方程式を弱形式へ変換すること
• 棒を節点と要素に分けること
• 試験関数、基底関数、形状関数を対応付けること
• 節点変位から要素ひずみを求めること
• 弱形式から要素剛性行列と等価熱荷重を作ること
• 要素剛性を全体剛性行列へ組み立てること
• 熱ひずみを等価節点力として入れること

　を確認します。

　今回の記事の目的は、いきなり複雑なFEMプログラムを作ることではありません。 材料点で行っていた「ひずみから応力を計算する」という操作が、FEMの中ではどこに現れるのかを見えるようにすることです。 この入口を作っておくと、次の記事で各積分点へ熱弾塑性の応力更新を入れたときにも、式の意味を追いやすくなります。

熱弾塑性解析の基礎⑤：拘束条件が及ぼす影響【MATLABコード付き】

はじめに　これまでの記事では、両端を完全に固定した1次元棒を主な対象としてきました。完全拘束では、温度上昇によって材料が自由に伸びようとしても、棒全体の伸びが許されません。そのため、熱ひずみが弾性ひずみとして打ち消され、熱応力が発生します。...

taku-naka.site

2026.07.21

解析対象と支配方程式

　長さ \(L\)、断面積 \(A\) の1次元棒を考えます。棒の軸方向を \(x\)、軸方向変位を \(u(x)\)、軸方向応力を \(\sigma(x)\) とします。単位長さ当たりの分布荷重を \(p(x)\) とし、右向きを正とします。

　支配方程式を導くため、図1のように、棒から \(x\) と \(x+\Delta x\) に挟まれた微小区間を取り出します。断面に作用する軸力を \(N(x)=A\sigma(x)\) と置き、引張軸力を正とします。正の引張軸力は微小区間を両側へ引っ張るため、左断面では左向きに \(N(x)\)、右断面では右向きに \(N(x+\Delta x)\) が作用します。一方、分布荷重 \(p(x)\) が微小区間内でほぼ一定なら、その合力は \(p(x)\Delta x\) です。

　微小区間に作用する力を右向き正で足し合わせると、左断面の軸力は負、右断面の軸力と分布荷重の合力は正です。したがって、有限の長さ \(\Delta x\) に対する力の釣合いから始め、\(\Delta x\) で割って極限を取ると、

となります。1行目は微小区間全体の力の釣合い、2行目は軸力の変化を区間長で割った差分商、3行目と4行目はその極限を微分として表したものです。最後の行では \(N=A\sigma\) を代入しています。つまり式\((1)\)は、「単位長さ当たりの軸力の変化」と「単位長さ当たりの外力」が釣り合うことを表しています。

　たとえば分布荷重がなければ \(p(x)=0\) です。このとき式\((1)\)から \(dN/dx=0\) となるため、軸力 \(N=A\sigma\) は棒の中で一定です。逆に \(p(x)\neq 0\) なら、分布荷重を受け持つために軸力は位置 \(x\) とともに変化します。図1で左右の軸力を最初から同じ大きさにしていないのは、この変化を許しているためです。

　支配方程式だけでは変位は一意に決まらないため、棒の端には境界条件を与えます。変位を指定する境界を \(\Gamma_u\)、端部力を指定する境界を \(\Gamma_t\) とすると、

です。ここで \(n\) は棒端の外向き法線で、左端では \(n=-1\)、右端では \(n=1\) です。\(\bar{u}\) は既知変位、\(\bar{t}\) は既知の端部力です。式\((1)\)と式\((2)\)を合わせたものが、今回FEMで解く境界値問題です。

支配方程式を弱形式へ変換する

　有限要素法では、式\((1)\)をそのまま節点へ当てはめるのではなく、まず領域全体の積分式へ変換します。変位境界 \(\Gamma_u\) で0となる任意の試験関数 \(v(x)\) を用意し、支配方程式へ掛けて積分すると、

となります。ここで \(v\) は、支配方程式の残差を領域内のどのような分布で調べるかを指定する重みです。

　式\((3)\)のままでは、未知量である軸力 \(A\sigma\) を \(x\) で微分しなければなりません。そこで、軸力の微分を含む第1項へ部分積分を適用し、微分を \(A\sigma\) から試験関数 \(v\) へ移します。

　まず、式\((3)\)の角括弧を外し、軸力の微分を含む項と分布荷重の項に分けます。

　式\((4)\)の第1項だけに部分積分を適用します。 一般の部分積分公式 \(\int f g’ dx=[fg]-\int f’g dx\) において、\(f=v\)、\(g=A\sigma\) と置きます。 ただし、第1項には積分全体の前に負号があるため、部分積分後の各項の符号も反転します。

　式\((5)\)では、微分が未知の軸力 \(A\sigma\) から既知の試験関数 \(v\) へ移っています。 この結果を式\((4)\)の第1項へ代入すると、

となります。 式\((6)\)の境界項を右辺へ移し、分布荷重の項も右辺へ移すと、

と整理できます。 式\((7)\)の \([vA\sigma]_0^L\) は、左右の棒端で評価される境界項です。 変位を指定する境界 \(\Gamma_u\) では試験関数が \(v=0\) であり、端部力を指定する境界 \(\Gamma_t\) では式\((2)\)の \(nA\sigma=\bar{t}\) を使えます。 したがって、境界項を既知の端部力で書き換えると、

となります。式\((8)\)が今回の弱形式です。左辺は棒内部の応力による作用、右辺は分布荷重と端部力による作用を表します。強形式では応力から作る軸力の微分が必要でしたが、弱形式では試験関数の1階微分と応力そのものを積分すればよくなります。

　この段階では、まだ棒を要素へ分けていません。ここから試験関数と変位を有限個の基底関数で近似し、積分を要素ごとに分解することで、有限要素方程式を作ります。

棒を節点と要素に分ける

　前節で定めた連続体の棒を、有限個の要素に分けます。 1次元棒では、もっとも基本的な要素として、両端に2つの節点を持つ線形棒要素を使えます。

図2 1次元棒を節点と要素に分けた模式図。橙色の矢印は各節点の変位自由度 \(u_i\) とその正方向を表します。要素の中では節点変位を補間して変位場を表します。

　図2のように棒を分割すると、未知量は各節点の変位になります。 つまり、連続関数 \(u(x)\) をそのまま未知にするのではなく、節点変位

　式\((9)\)の節点変位ベクトルを未知量として解きます。 この式は、FEMでは変位場全体を節点値の集合として近似することを表しています。 未知量が有限個になるため、連立一次方程式として解けるようになります。

基底関数と形状関数で変位を近似する

　式\((9)\)で未知量を有限個の節点変位 \(u_i\) に置き換えましたが、ひずみや応力を求めるには、節点間を含む棒全体の変位分布も必要です。そこで、節点 \(i\) に対応する全体基底関数 \(g_i(x)\) を使い、節点変位から連続した変位場を再構成します。

　ここで用いる \(\tilde{u}(x)\) の上の波線は、厳密な変位場 \(u(x)\) そのものではなく、有限個の節点変位から作った近似変位場であることを示します。節点上では \(\tilde{u}(x_i)=u_i\) となりますが、節点間の分布は選んだ基底関数によって近似されます。この近似変位場と試験関数を

と表します。式\((10)\)の左側では、各節点変位 \(u_i\) に、その節点の影響範囲を表す \(g_i(x)\) を掛け、すべての節点について足し合わせています。右側では、ガラーキン法として試験関数 \(v_i\) を同じ全体基底関数 \(g_i\) に一致させています。したがって、弱形式へ導入した試験関数は消えたのではなく、各節点の方程式を作る基底関数として使われます。

　次に、1つの要素だけを取り出して考えます。 要素長さを \(L_e\)、左節点変位を \(u_1\)、右節点変位を \(u_2\) とします。 2節点線形要素では、要素内の変位を節点変位の線形補間で表します。

　式\((11)\)において、\(N_1\) と \(N_2\) は形状関数です。 要素左端を \(x=0\)、右端を \(x=L_e\) と置けば、

　式\((12)\)が、今回の2節点棒要素で使う形状関数です。 記号 \(g_i|_e\) は、棒全体で定義された基底関数 \(g_i\) を、要素 \(e\) の内部だけで見ることを表します。つまり、全体では \(g_i\) と呼んでいた関数を、要素内では形状関数 \(N_i\) として扱います。 この式は、左端では \(N_1=1, N_2=0\)、右端では \(N_1=0, N_2=1\) になるように作られています。 つまり、要素内の任意の位置の変位を、左右の節点変位からなめらかにつないでいます。

　図3を見ると、\(N_1\) と \(N_2\) が足して常に1になることが分かります。 これは、要素内の変位が左右の節点変位の重み付き平均として表されることを意味します。

　試験関数についても同じ形状関数を使います。2節点要素では、形状関数を並べた行ベクトルと任意係数ベクトルを

と定義します。ここで \(c_1,c_2\) は、弱形式がどのような試験関数に対しても成り立つことを表す任意の係数です。また、\(\boldsymbol{B}=d\boldsymbol{N}/dx\) は形状関数を \(x\) で微分した行ベクトルであり、次節で具体的な成分を求めます。式\((13)\)によって、弱形式の試験関数 \(v\) が、後で要素行列に現れる \(\boldsymbol{N}^{\mathrm T}\) と \(\boldsymbol{B}^{\mathrm T}\) へつながります。

節点変位からひずみを求める

　1次元棒の全ひずみは、変位を \(x\) で微分したものです。

　式\((14)\)へ式\((11)\)を代入して微分すると、

となります。 式\((15)\)は、要素両端の変位差を要素長さで割ると、要素ひずみが得られることを表しています。 行列で書けば、

　式\((16)\)において、\(\boldsymbol{B}\) は節点変位をひずみに変換する行列です。 FEMの中で材料モデルが使われるのは、基本的にこのあとです。 節点変位からひずみを求め、そのひずみから応力を計算します。

熱ひずみを含む応力

　温度変化がある場合、棒は自由に伸びようとします。 線膨張係数を \(\alpha\)、温度上昇を \(\Delta T\) とすると、熱ひずみは

　式\((17)\)が、温度上昇から自由膨張量をひずみとして求める関係です。 ここまでと同じく、棒に実際に生じる全ひずみを \(\varepsilon\) と書き、応力を生む弾性ひずみを新たに \(\varepsilon_{\mathrm{ela}}\) と書きます。 全ひずみは、弾性ひずみと熱ひずみの和として、

　と分解します。 式\((18)\)の1行目がひずみの加算分解で、2行目は弾性ひずみについて解き直した形です。 添字 \(\mathrm{ela}\) は弾性ひずみ、\(\mathrm{th}\) は熱ひずみを表します。 全ひずみには添字を付けず、これまでどおり \(\varepsilon\) と表すことで、前節までの式とのつながりを保ちます。 また、ひずみの種類を表す \(\mathrm{ela}\) と、後で要素番号として使う \(e\) を区別します。

　線形弾性では弾性ひずみにヤング率を掛けて応力を求めるため、式\((18)\)を使うと、

となります。 この式は、棒が熱で伸びたい量と、実際に許された変形量の差が応力になることを表しています。 自由に伸びられるなら \(\varepsilon=\varepsilon_{\mathrm{th}}\) となり、弾性ひずみ \(\varepsilon_{\mathrm{ela}}\) と応力はゼロです。

　FEMでは、この一般式を要素ごとに使います。 要素 \(e\) の全ひずみは式\((16)\)より \(\varepsilon_e=\boldsymbol{B}\boldsymbol{u}_e\) です。 また、要素の線膨張係数を \(\alpha_e\)、温度上昇を \(\Delta T_e\) とすると、熱ひずみは式\((17)\)より \(\varepsilon_{\mathrm{th},e}=\alpha_e\Delta T_e\) です。 したがって、要素の弾性ひずみは \(\varepsilon_{\mathrm{ela},e}=\varepsilon_e-\varepsilon_{\mathrm{th},e}\) です。 ここで、コンマの後の \(e\) は要素番号を表し、\(\mathrm{ela}\) は弾性ひずみという種類を表します。 これらを式\((19)\)へ代入すると、

となります。 式\((20)\)が、後で一要素の弱形式へ代入する要素応力です。 一般式と比べて新しい材料則を加えたのではなく、全ひずみを節点変位で表しただけです。

　完全に伸びを拘束されるなら \(\varepsilon=0\) となり、圧縮応力

が発生します。 この関係は、以前の記事で扱った完全拘束棒の弾性応答と同じです。 FEMでは、この考えを各要素で使います。

弱形式から熱を含む一要素方程式を作る

　ここから、式\((8)\)の弱形式を要素ごとの行列へ変換します。領域全体の積分を要素区間の積分の和へ分け、そのうち一つの要素 \(e\) を取り出します。

　まず変位について、式\((11)\)の \(u^e=N_1u_1+N_2u_2\) へ、式\((13)\)で定義した \(\boldsymbol{N}=[N_1\ N_2]\) と、要素節点変位ベクトル \(\boldsymbol{u}_e=[u_1\ u_2]^{\mathrm T}\) を当てはめます。すると、スカラーの補間式は行列積 \(\boldsymbol{N}\boldsymbol{u}_e\) と同じです。これを \(x\) で微分すると、節点変位は \(x\) に依存しないため、形状関数だけが微分されて \(\boldsymbol{B}\boldsymbol{u}_e\) になります。

　試験関数についても、式\((13)\)で同じ操作をすでに行っています。以上をスカラー表示から行列表示まで並べると、

となります。式\((22)\)は新しい仮定ではなく、式\((11)\)と式\((13)\)を、後の積分へ代入しやすい行列形に書き直したものです。

　次に、弱形式\((8)\)から要素 \(e\) に対応する部分だけを取り出します。

　式\((23)\)は、棒全体の弱形式を一要素の積分区間 \(x_1^e\leq x\leq x_2^e\) に制限したものです。

　ここへ式\((22)\)の \(v^e=\boldsymbol{N}\boldsymbol{c}_e\) と \(dv^e/dx=\boldsymbol{B}\boldsymbol{c}_e\) を代入します。 また、応力 \(\sigma_e\) には7節で導いた要素応力の式\((20)\)を代入します。 スカラー \(\boldsymbol{B}\boldsymbol{c}_e\) は \((\boldsymbol{c}_e)^{\mathrm T}\boldsymbol{B}^{\mathrm T}\) と書き直せるため、

となります。 式\((24)\)の右辺の角括弧は、分布荷重と既知端部力を要素の2節点へ振り分けた量です。 ここで、この角括弧全体を等価節点外力ベクトルとして

と置きます。 式\((25)\)の第1項は分布荷重を左右の節点へ配分した量、第2項は要素端が外力境界に一致する場合の既知端部力です。 機械的な分布荷重も端部力もなければ、\(\boldsymbol{f}^{\mathrm{ext}}_e=\boldsymbol{0}\) です。

　式\((24)\)の右辺を \((\boldsymbol{c}_e)^{\mathrm T}\boldsymbol{f}^{\mathrm{ext}}_e\) と書き換えます。 そのうえで左辺を、節点変位に比例する項と熱ひずみによる項へ分けます。 さらに、\(\boldsymbol{c}_e\) と \(\boldsymbol{u}_e\) は積分変数 \(x\) に依存しないため、積分の外へ出すと、

を得ます。 式\((26)\)の1行目では、左辺第1項が節点変位による弾性抵抗、左辺第2項が自由熱膨張による応力の減少、右辺が機械的外力を表します。 2行目では、1行目の右辺を左辺へ移し、すべての項に共通する \((\boldsymbol{c}_e)^{\mathrm T}\) をくくり出しています。

　ここで、\(\boldsymbol{c}_e\) は任意です。したがって、式\((26)\)の2行目がすべての \(\boldsymbol{c}_e\) に対して成り立つには、中括弧の中がゼロでなければなりません。中括弧内の第1項を \(\boldsymbol{k}_e\boldsymbol{u}_e\)、第3項を \(\boldsymbol{f}^{\mathrm{th}}_e\) と定義し、それぞれを整理すると、

を得ます。式\((27)\)の1行目が、支配方程式から弱形式を経て得られた一要素の熱弾性方程式です。2行目は1行目に現れる要素剛性行列 \(\boldsymbol{k}_e\)、3行目は等価熱荷重 \(\boldsymbol{f}^{\mathrm{th}}_e\) の定義です。試験関数 \(v\) は消えたのではなく、その微分が \(\boldsymbol{B}^{\mathrm T}\)、その値が機械荷重側の \(\boldsymbol{N}^{\mathrm T}\) として要素積分に残っています。

　式\((27)\)は一要素の方程式ですが、\(\boldsymbol{k}_e\)、\(\boldsymbol{f}^{\mathrm{th}}_e\)、\(\boldsymbol{f}^{\mathrm{ext}}_e\) の具体的な成分はまだ計算していません。ここからは別の考え方へ移らず、式\((27)\)の3項を2節点線形棒要素について順番に計算し、その結果を全体方程式へ組み立てます。

　要素 \(e\) の左端を \(x=0\)、右端を \(x=L_e\) とし、要素内では \(E_e\)、\(A_e\)、\(\alpha_e\)、\(\Delta T_e\) が一定であるとします。式\((12)\)と式\((16)\)から、計算に使う形状関数、ひずみ–変位行列、節点変位ベクトルは

です。

　まず、式\((27)\)の要素剛性行列へ式\((28)\)の \(\boldsymbol{B}\) を代入します。

　列ベクトルと行ベクトルを掛けると、\(2\times2\) 行列になります。

　式\((30)\)を式\((29)\)へ入れます。被積分関数は要素内で一定なので、積分によって要素長さ \(L_e\) が一つ掛かります。

　次に、式\((27)\)の等価熱荷重を計算します。式\((28)\)の \(\boldsymbol{B}^{\mathrm T}\) を代入すると、

となります。材料定数と温度上昇を積分の外へ出し、\(\int_0^{L_e}dx=L_e\) を使うと、

を得ます。\(q_e\) は、熱ひずみを節点力の形へ置き換えたときの大きさです。左節点成分が \(-q_e\)、右節点成分が \(+q_e\) となるため、要素を自由に伸ばす向きのベクトルになっています。

　残る機械的外力は、式\((25)\)で定義した \(\boldsymbol{f}^{\mathrm{ext}}_e\) です。一般には分布荷重と既知端部力から計算し、その2成分を \(F^{\mathrm{ext}}_{1,e}\)、\(F^{\mathrm{ext}}_{2,e}\) と書きます。たとえば要素内で一定の分布荷重 \(p_e\) だけが作用する場合には、

となります。一様な分布荷重の合力 \(p_eL_e\) が、左右の節点へ半分ずつ配分されています。端部力がある場合は、その端部に対応する成分へ \(\bar{t}\) を加えます。

　以上で式\((27)\)の3項がすべて成分表示できました。式\((31)\)、式\((33)\)、式\((34)\)を式\((27)\)へ入れると、

という一要素の方程式になります。行ごとに書けば、

です。式\((36)\)は、要素左端と右端における節点力の釣合いをそれぞれ表しています。

　最後に、この一要素の成分を全体行列へ足します。要素 \(e\) が全体節点 \(i=e\) と \(j=e+1\) を結ぶ場合、式\((35)\)から次の加算規則が得られます。

　式\((37)\)を要素1、2、3について順に実行すると、次節の式\((38)\)が得られます。つまり、式\((38)\)は突然導入する全体行列ではなく、式\((27)\)を具体的に積分し、その \(2\times2\) の要素成分を共有節点へ足し合わせた結果です。

全体剛性行列への組み立て

　前節では、式\((27)\)の積分を具体的に計算し、式\((35)\)の一要素方程式を得ました。さらに、式\((37)\)として、要素成分を全体節点番号に対応する位置へ足す規則を示しました。この操作を組み立てと呼びます。

　すべての要素を足し込むと、全体のつり合い式は \(\boldsymbol{K}\boldsymbol{u}=\boldsymbol{F}\) になります。 ここで \(\boldsymbol{K}\) は全体剛性行列、\(\boldsymbol{u}\) は全節点変位ベクトル、\(\boldsymbol{F}\) は外力と等価熱荷重を含む全体荷重ベクトルです。

　行列の中身を見えるようにするため、3要素を直列につないだ4節点の棒を例にします。 要素 \(e\) の軸剛性を \(k_e=E_eA_e/L_e\)、熱膨張による等価節点力の大きさを \(q_e=E_eA_e\alpha_e\Delta T_e\) と置きます。 このとき、全体系の式は

になります。 式\((38)\)の左辺を見ると、要素1の剛性は節点1と2に対応する位置、要素2の剛性は節点2と3に対応する位置へ入っています。 節点2は要素1と要素2の両方に共有されるため、対角成分は \(k_1+k_2\) になります。 同じように、節点3の対角成分は \(k_2+k_3\) になります。

　右辺でも同じ組み立てが行われます。 要素1の等価熱荷重は節点1へ \(-q_1\)、節点2へ \(+q_1\) として入り、要素2の等価熱荷重は節点2へ \(-q_2\)、節点3へ \(+q_2\) として入ります。 したがって、共有節点2の熱荷重成分は \(q_1-q_2\) です。 外力が別に作用する場合は、その節点の \(F_i^{\mathrm{ext}}\) を加えます。

　たとえば式\((38)\)の2行目だけを取り出すと、 \(-k_1u_1+(k_1+k_2)u_2-k_2u_3=F^{\mathrm{ext}}_2+q_1-q_2\) です。 これは、節点2に接続する要素1と要素2が生む節点力の合計と、節点2へ作用する外力および等価熱荷重がつり合うことを表しています。 要素数が増えても、この「隣接要素の成分を共有節点へ足し込む」という規則は変わりません。

図4 8要素の1次元棒で組み立てた全体剛性行列。隣り合う節点だけが同じ要素を共有するため、非ゼロ成分は対角付近に集まります。

　図4では、全体剛性行列の非ゼロ成分が対角付近に並んでいます。 これは、1次元棒要素では、ある節点が直接つながるのは隣の節点だけだからです。 FEMの全体行列は、要素ごとの局所的な関係を足し込んで作られるため、このような疎な構造になります。

　図5では、要素を1つずつ追加すると、全体剛性行列のどこに成分が入るかを確認できます。 FEMの組み立ては抽象的に見えますが、実際には各要素の \(2\times2\) 行列を、対応する節点番号の場所へ足すだけです。

境界条件を入れて変位を解く

　式\((38)\)をそのまま解くことはできません。 剛体移動を止めるために、少なくとも1つの変位を固定する必要があります。 たとえば左端節点を固定するなら、

　式\((39)\)を境界条件として与えます。 両端固定なら、

　式\((40)\)の二つの変位を境界条件として与えます。 既知変位を与えた自由度を取り除き、未知変位だけに対する連立方程式を解くと、全節点変位が得られます。 次に、この節点変位から要素のひずみと応力を評価します。

応力とひずみを評価する位置

　ここまでに求めた未知量は、式\((38)\)の節点変位 \(\boldsymbol{u}\) です。 応力とひずみは直接解く未知量ではなく、解けた節点変位を各要素へ戻し、式\((28)\)の \(\boldsymbol{B}\) と構成式\((19)\)を使って評価します。

　この評価位置を積分点と呼ぶのは、式\((27)\)の要素剛性や等価熱荷重が、もともと要素内の積分として定義されているためです。 一般の要素では、これらの積分を次の数値積分で計算します。

　ここで、\(\Omega_e\) は要素領域、\(x_i\) は積分点、\(w_i\) はその重み、\(q(x)\) は積分する量です。 したがって積分点は、結果を見るために後から任意に選ぶ点ではなく、要素方程式の積分を計算するために、要素タイプと積分次数に応じて先に決める点です。

　今回の2節点棒要素では、\(\boldsymbol{B}\) が要素内で一定です。 また、今回は温度上昇 \(\Delta T_e\) も要素内で一定としています。 このため、要素中央の1点 \(x_{\mathrm{gp}}=L_e/2\) で評価すれば十分です。 式\((28)\)の \(\boldsymbol{B}\) をそのまま使い、

と計算します。 ここで、\(\boldsymbol{B}_{\mathrm{gp}}=\boldsymbol{B}\) であり、添字 \(\mathrm{gp}\) は積分点で評価した量を表します。 \(\boldsymbol{B}\) と \(\Delta T_e\) が要素内で一定なので、式\((42)\)から位置 \(x_{\mathrm{gp}}\) が消え、得られるひずみと応力は要素内で一定になります。 したがって、今回の「積分点応力」はそのまま「要素応力」として扱えます。

　高次要素、多次元要素、あるいは材料状態が要素内で変化する解析では、\(\boldsymbol{B}\) や応力が位置によって変わるため、積分点ごとに異なる値を持ちます。 節点応力として表示される値は、通常、積分点の応力を外挿または平均した後処理の値です。 したがって、「節点変位は直接解く量」「応力とひずみは積分点で評価する量」と区別して考えることが大切です。

支配方程式から計算手順までの全体像

　ここまでの定式化を図6にまとめます。最初に支配方程式と境界条件を定め、弱形式へ変換します。その後でガラーキン法として試験関数を基底関数と一致させ、要素内では基底関数を形状関数として扱います。この順序を経て初めて、一要素の剛性行列と等価熱荷重が得られます。

　図6の上半分は、連続体の境界値問題を有限個の要素方程式へ変換する定式化です。下半分は、得られた要素方程式を組み立て、境界条件を反映して節点変位と積分点応力を求める数値計算です。この流れが、今回のMATLABコードの骨格です。 重要なのは、材料モデルは全体方程式を解いたあとにだけ出てくるのではない、という点です。 要素剛性を作るところでも \(E\) が入り、応力を評価するところでも \(E\) と熱ひずみが入ります。 熱弾塑性へ拡張すると、この「応力を評価するところ」が応力更新アルゴリズムに置き換わります。

数値例の解析条件

　ここからは、同じ1次元棒を使って計算結果を確認します。 棒は長さ \(100\) mm、断面積 \(10\) mm\(^2\) とし、8個の要素に分割します。 材料は線形弾性とし、塑性変形はまだ考えません。

　比較する計算ケースは3つです。 1つ目は、左端を固定し、右端に引張荷重 \(P\) を与えるケースです。 2つ目は、左端だけを固定し、棒全体を一様に加熱するケースです。 この場合、右端は自由に動けるので、棒は自由膨張できます。 3つ目は、両端を固定し、棒全体を一様に加熱するケースです。 この場合、熱膨張が拘束されるため、圧縮応力が発生します。 この3つの条件を模式的に描くと、図7のようになります。 同じ棒を使っていても、外力で伸ばすのか、熱で自由に伸びるのか、熱で伸びたいのに両端で止められるのかで、応力の意味が変わります。

計算結果

代表値の確認

　表2に、得られた代表値を示します。

　右端荷重のケースでは、棒全体が一様に引張られるため、要素応力は \(P/A=100\) MPa になります。 右端変位は \(PL/(EA)\) に一致します。 これは、連続体の棒の理論解と同じです。

　左端固定の自由熱膨張では、右端変位が \(L\alpha\Delta T=0.120\) mm になります。 棒は熱で伸びたい量だけ伸びられるため、理論上の応力はゼロです。 数値計算結果に極めて小さな値が出る場合は、物理的な応力ではなく丸め誤差として扱います。 一方、両端固定の一様加熱では、全体の伸びがゼロに拘束されます。 そのため、式\((21)\)の通り、応力は \(-E\alpha\Delta T=-252\) MPa になります。

変位、応力、ひずみ成分

　図8は、上段と下段で見ている量が違います。 上段左は3つのケースの節点変位、上段右は3つのケースの要素応力を比較しています。 一方、下段の3つの棒グラフは、各ケースで全ひずみ、熱ひずみ、弾性ひずみがどのように分かれているかを示しています。 したがって、下段左がケース1の右端荷重、下段中央がケース2の自由熱膨張、下段右がケース3の両端固定加熱に対応します。

　ケース1の右端荷重では、温度変化を与えていないので熱ひずみはゼロです。 そのため、外力によって生じた全ひずみがそのまま弾性ひずみになり、一様な引張応力を生みます。 ケース2の自由熱膨張では、棒が熱ひずみの分だけ自由に伸びられます。 この場合、全ひずみは熱ひずみと一致し、弾性ひずみはゼロです。 弾性ひずみがゼロなので、応力もゼロになります。

　ケース3の両端固定加熱では、節点変位が拘束されるため全ひずみはゼロです。 しかし、材料としては熱ひずみが発生しています。 全ひずみをゼロに保つためには、熱ひずみと同じ大きさで符号が反対の弾性ひずみが必要になります。 この負の弾性ひずみが圧縮応力を生みます。 つまり、FEMで計算しても、熱応力の本質はこれまでの記事と同じです。 違うのは、その関係を各要素で評価できる形にしたことです。

　図9では、節点変位を使って棒の変形を表示します。 ここで動いている量は、実時間に沿った変位ではありません。 今回の計算は各ケースの静的なつり合い問題なので、得られるのは最終的な節点変位です。 アニメーションでは、その最終変位に表示倍率をかけ、未変形形状から最終変形形状まで連続的に補間しています。 したがって、膨張したあとに物理的に戻る挙動を表しているわけではありません。 右端荷重と自由熱膨張はどちらも右端が動きますが、応力の意味は異なります。 右端荷重では弾性ひずみが応力を生みます。 自由熱膨張では熱ひずみによって伸びているだけなので、応力は生じません。 両端固定では変位が見た目には動かないにもかかわらず、内部には圧縮応力が発生します。 このように、変位だけでなく、ひずみ分解と応力を同時に見ることが大切です。

MATLABコード

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

% fem_1d_introduction_demo.m
% 1次元2節点棒要素による線形熱弾性FEM
%
% 次の3ケースを同じ要素方程式で解析します。
%   case 1: 左端固定、右端に機械荷重 P
%   case 2: 左端固定、右端自由、一様温度上昇
%   case 3: 両端固定、一様温度上昇
%
% 計算手順は以下のとおりです。
%   1. 節点座標と要素接続を作る
%   2. 要素剛性行列と等価熱荷重を全体系へ組み立てる
%   3. 変位境界条件を適用して節点変位を解く
%   4. 各要素で全ひずみ、熱ひずみ、弾性ひずみ、応力を求める
%
% 単位系は N, mm, MPa (= N/mm^2), K で統一しています。

clear; close all; clc;

%% 1. 材料定数と棒の寸法
E = 210000;       % ヤング率 [MPa]
A = 10.0;         % 断面積 [mm^2]
L = 100.0;        % 棒の長さ [mm]
alpha = 12e-6;    % 線膨張係数 [1/K]

%% 2. 節点と要素を作る
nElem = 8;        % 要素数
nNode = nElem + 1;

x = linspace(0, L, nNode).';                  % 節点座標 [mm]
conn = [(1:nElem).', (2:nElem + 1).'];        % 各行: [左節点, 右節点]

%% 3. 荷重条件と温度条件
P = 1000;         % case 1の右端荷重 [N]
dTUniform = 100;  % case 2, 3の一様温度上昇 [K]

temperatureZero = zeros(nElem, 1);
temperatureUniform = dTUniform * ones(nElem, 1);

%% 4. 出力設定
saveFigures = true;
saveDeformationAnimation = false;  % MP4が必要な場合だけtrueにします
saveAssemblyAnimation = false;     % MP4が必要な場合だけtrueにします
videoFrameRate = 15;

deformationMp4FileName = "fem_1d_deformation_animation.mp4";
assemblyMp4FileName = "fem_1d_assembly_animation.mp4";

%% 5. 3つの境界条件を解析する
% case 1: 左端変位を0に固定し、右端へ荷重Pを与えます。
caseLoad = solve_bar_case( ...
    x, conn, E, A, alpha, temperatureZero, ...
    1, 0, nNode, P);

% case 2: 左端だけを固定します。右端は自由なので熱膨張できます。
caseFreeThermal = solve_bar_case( ...
    x, conn, E, A, alpha, temperatureUniform, ...
    1, 0, [], []);

% case 3: 左右端の変位を0に固定し、熱膨張を拘束します。
caseFixedThermal = solve_bar_case( ...
    x, conn, E, A, alpha, temperatureUniform, ...
    [1, nNode], [0; 0], [], []);

%% 6. 代表値を表示する
fprintf('=== 1D bar FEM introduction ===\n');
fprintf('Mechanical load case: right-end displacement = %.6f mm\n', caseLoad.u(end));
fprintf('Mechanical load case: element stress = %.3f MPa\n', mean(caseLoad.sigma));
fprintf('Free thermal case: right-end displacement = %.6f mm\n', caseFreeThermal.u(end));
fprintf('Free thermal case: max abs stress = %.3e MPa\n', max(abs(caseFreeThermal.sigma)));
fprintf('Fixed thermal case: right-end displacement = %.6f mm\n', caseFixedThermal.u(end));
fprintf('Fixed thermal case: element stress = %.3f MPa\n', mean(caseFixedThermal.sigma));

%% 7. 図を保存する
if saveFigures
    make_mesh_concept_figure(x, conn);
    make_shape_function_figure();
    make_assembly_matrix_figure(x, conn, E, A);
    make_result_summary_figure(x, caseLoad, caseFreeThermal, caseFixedThermal, dTUniform);
end

%% 8. 必要な場合だけMP4を保存する
if saveDeformationAnimation
    make_deformation_animation(x, caseLoad, caseFreeThermal, caseFixedThermal, ...
        deformationMp4FileName, videoFrameRate);
end

if saveAssemblyAnimation
    make_assembly_animation(x, conn, E, A, assemblyMp4FileName, videoFrameRate);
end

%% 9. FEM解析関数
function result = solve_bar_case(x, conn, E, A, alpha, elementTemperature, ...
        prescribedDofs, prescribedValues, forceDofs, forceValues)
    % 1次元棒の1ケースを解析します。
    %
    % 入力:
    %   x                  節点座標 [mm]
    %   conn               要素接続行列 [左節点, 右節点]
    %   E, A, alpha        材料定数と断面積
    %   elementTemperature 各要素の温度上昇 [K]
    %   prescribedDofs     変位を指定する節点番号
    %   prescribedValues   指定変位 [mm]
    %   forceDofs          節点荷重を与える節点番号
    %   forceValues        節点荷重 [N]
    %
    % 出力resultには、節点変位、要素ひずみ、要素応力、K、Fを格納します。

    nNode = numel(x);
    nElem = size(conn, 1);
    K = zeros(nNode, nNode);  % 全体剛性行列 [N/mm]
    F = zeros(nNode, 1);      % 全体荷重ベクトル [N]

    % --- 要素行列と等価熱荷重の組み立て ---
    for e = 1:nElem
        nodes = conn(e, :);
        x1 = x(nodes(1));
        x2 = x(nodes(2));
        Le = x2 - x1;

        % 2節点線形棒要素ではB行列が要素内で一定です。
        Be = [-1 / Le, 1 / Le];

        % Bが一定なので、積分は要素長さLeを掛ければ計算できます。
        % ke = integral(B' E A B dx) = B' E A B Le
        ke = Be.' * E * A * Be * Le;

        % fth = integral(B' E alpha*dT A dx) = B' E alpha*dT A Le
        dTe = elementTemperature(e);
        fth = Be.' * E * alpha * dTe * A * Le;

        % 要素eの2×2行列と2×1ベクトルを、対応する全体成分へ加えます。
        K(nodes, nodes) = K(nodes, nodes) + ke;
        F(nodes) = F(nodes) + fth;
    end

    % 機械的な節点荷重を全体荷重ベクトルへ加えます。
    if ~isempty(forceDofs)
        F(forceDofs) = F(forceDofs) + forceValues(:);
    end

    % 既知変位を考慮し、未知節点変位を解きます。
    u = solve_with_prescribed_displacements( ...
        K, F, prescribedDofs, prescribedValues);

    % --- 解いた節点変位から要素量を回収する後処理 ---
    epsTotal = zeros(nElem, 1);
    epsThermal = zeros(nElem, 1);
    epsElastic = zeros(nElem, 1);
    sigma = zeros(nElem, 1);
    xElem = zeros(nElem, 1);

    for e = 1:nElem
        nodes = conn(e, :);
        Le = x(nodes(2)) - x(nodes(1));
        ue = u(nodes);

        % 全ひずみ epsilon_e = B_e u_e
        Be = [-1 / Le, 1 / Le];
        epsTotal(e) = Be * ue;

        % epsilon = epsilon_ela + epsilon_th の加算分解です。
        epsThermal(e) = alpha * elementTemperature(e);
        epsElastic(e) = epsTotal(e) - epsThermal(e);

        % 線形弾性構成式 sigma = E*epsilon_ela
        sigma(e) = E * epsElastic(e);

        % 今回はBが一定なので、要素中央を代表位置として保存します。
        xElem(e) = mean(x(nodes));
    end

    % 結果を一つの構造体にまとめます。
    result.u = u;
    result.epsTotal = epsTotal;
    result.epsThermal = epsThermal;
    result.epsElastic = epsElastic;
    result.sigma = sigma;
    result.xElem = xElem;
    result.K = K;
    result.F = F;
end

function u = solve_with_prescribed_displacements( ...
        K, F, prescribedDofs, prescribedValues)
    % 変位が既知の自由度と未知の自由度に分け、縮約方程式を解きます。
    n = numel(F);
    allDofs = (1:n).';
    prescribedDofs = prescribedDofs(:);
    freeDofs = setdiff(allDofs, prescribedDofs);

    u = zeros(n, 1);
    u(prescribedDofs) = prescribedValues(:);

    % K_ff*u_f = F_f - K_fp*u_p
    Ffree = F(freeDofs) ...
        - K(freeDofs, prescribedDofs) * u(prescribedDofs);
    u(freeDofs) = K(freeDofs, freeDofs) \ Ffree;
end

%% 10. 描画関数
function make_mesh_concept_figure(x, conn)
    % 節点、要素、節点変位自由度の模式図を保存します。
    fig = figure('Color', 'w', 'Position', [120 120 1200 300]);
    hold on; axis off;
    y = 0.5;
    plot([x(1), x(end)], [y, y], 'Color', [0.15 0.35 0.75], 'LineWidth', 12);
    plot(x, y * ones(size(x)), 'ko', 'MarkerFaceColor', 'w', 'MarkerSize', 8, 'LineWidth', 1.6);
    Le = min(diff(x));
    arrowColor = [0.88 0.30 0.10];
    for i = 1:numel(x)
        % quiverは横長の図で矢じりが不安定になるため、線と三角形で描きます。
        arrowStart = x(i) + 0.05 * Le;
        arrowTip = x(i) + 0.40 * Le;
        headLength = 0.075 * Le;
        headHalfHeight = 0.012;
        arrowY = y + 0.035;
        plot([arrowStart, arrowTip - headLength], [arrowY, arrowY], ...
            'Color', arrowColor, 'LineWidth', 1.8);
        patch([arrowTip, arrowTip - headLength, arrowTip - headLength], ...
            [arrowY, arrowY + headHalfHeight, arrowY - headHalfHeight], ...
            arrowColor, 'EdgeColor', arrowColor);
        text(x(i) + 0.225 * Le, y + 0.095, sprintf('u_{%d}', i), ...
            'Interpreter', 'tex', 'Color', [0.70 0.20 0.05], ...
            'HorizontalAlignment', 'center', 'FontSize', 15);
        text(x(i), y - 0.09, sprintf('%d', i), ...
            'HorizontalAlignment', 'center', 'FontSize', 15);
    end
    for e = 1:size(conn, 1)
        xc = mean(x(conn(e, :)));
        text(xc, y + 0.23, sprintf('e_{%d}', e), ...
            'Interpreter', 'tex', 'HorizontalAlignment', 'center', 'FontSize', 14);
    end
    text(x(1) - 1.18 * Le, y + 0.23, 'element', ...
        'HorizontalAlignment', 'left', 'FontSize', 16);
    text(x(1) - 1.18 * Le, y - 0.09, 'node', ...
        'HorizontalAlignment', 'left', 'FontSize', 16);
    xlim([x(1) - 1.32 * Le, x(end) + 0.75 * Le]);
    ylim([y - 0.14, y + 0.27]);
    hide_axes_toolbars(fig);
    exportgraphics(fig, 'fem_1d_mesh_concept.png', 'Resolution', 220);
    close(fig);
end

function make_shape_function_figure()
    % 2節点線形要素の形状関数N1, N2を保存します。
    xi = linspace(0, 1, 200);
    N1 = 1 - xi;
    N2 = xi;
    fig = figure('Color', 'w', 'Position', [120 120 980 620]);
    plot(xi, N1, 'Color', [0.1 0.3 0.75], 'LineWidth', 3.0); hold on;
    plot(xi, N2, 'Color', [0.85 0.25 0.15], 'LineWidth', 3.0);
    grid on; box on;
    xlabel('local coordinate x / L_e');
    ylabel('shape function');
    title('linear shape functions of a 2-node bar element');
    legend({'N_1 = 1 - x/L_e', 'N_2 = x/L_e'}, 'Location', 'east', 'FontSize', 15);
    set(gca, 'FontSize', 18, 'LineWidth', 1.0);
    ylim([-0.05 1.08]);
    hide_axes_toolbars(fig);
    exportgraphics(fig, 'fem_1d_shape_functions.png', 'Resolution', 220);
    close(fig);
end

function make_assembly_matrix_figure(x, conn, E, A)
    % 全要素を組み立てた後の全体剛性行列を可視化します。
    nNode = numel(x);
    K = zeros(nNode, nNode);
    fig = figure('Color', 'w', 'Position', [120 120 760 680]);
    for e = 1:size(conn, 1)
        nodes = conn(e, :);
        Le = x(nodes(2)) - x(nodes(1));
        ke = E * A / Le * [1 -1; -1 1];
        K(nodes, nodes) = K(nodes, nodes) + ke;
    end
    imagesc(K);
    axis equal tight;
    colormap(parula);
    colorbar;
    title('assembled global stiffness matrix');
    xlabel('node number');
    ylabel('node number');
    set(gca, 'FontSize', 17, 'LineWidth', 1.0);
    hide_axes_toolbars(fig);
    exportgraphics(fig, 'fem_1d_global_stiffness.png', 'Resolution', 220);
    close(fig);
end

function make_result_summary_figure(x, loadCase, freeThermal, fixedThermal, dTUniform)
    % 3ケースの変位、応力、ひずみ分解を一枚にまとめます。
    fig = figure('Color', 'w', 'Position', [100 100 1600 980]);
    tiledlayout(2, 3, 'Padding', 'loose', 'TileSpacing', 'loose');

    nexttile([1 2]);
    hold on;
    plot(x, loadCase.u, 'LineWidth', 2.7, 'Color', [0.1 0.3 0.75], 'DisplayName', 'case 1: end load');
    plot(x, freeThermal.u, 'LineWidth', 2.7, 'Color', [0.1 0.55 0.35], 'DisplayName', 'case 2: free thermal');
    plot(x, fixedThermal.u, 'LineWidth', 2.7, 'Color', [0.85 0.25 0.15], 'DisplayName', 'case 3: fixed thermal');
    grid on; box on;
    xlabel('x (mm)');
    ylabel('displacement u (mm)');
    title('nodal displacement: case comparison');
    legend('Location', 'northwest', 'FontSize', 13);
    set(gca, 'FontSize', 16, 'LineWidth', 1.0);

    nexttile;
    hold on;
    plot(loadCase.xElem, loadCase.sigma, '-o', 'LineWidth', 2.4, 'Color', [0.1 0.3 0.75], 'DisplayName', 'case 1: end load');
    plot(freeThermal.xElem, freeThermal.sigma, '-o', 'LineWidth', 2.4, 'Color', [0.1 0.55 0.35], 'DisplayName', 'case 2: free thermal');
    plot(fixedThermal.xElem, fixedThermal.sigma, '-o', 'LineWidth', 2.4, 'Color', [0.85 0.25 0.15], 'DisplayName', 'case 3: fixed thermal');
    yline(0, ':', 'Color', [0.3 0.3 0.3], 'HandleVisibility', 'off');
    grid on; box on;
    xlabel('element center x (mm)');
    ylabel('stress \sigma (MPa)');
    title('element stress: case comparison');
    legend('Location', 'best', 'FontSize', 13);
    set(gca, 'FontSize', 16, 'LineWidth', 1.0);

    nexttile;
    bar(loadCase.xElem, [loadCase.epsTotal, loadCase.epsThermal, loadCase.epsElastic] * 1000, 'grouped');
    grid on; box on;
    xlabel('element center x (mm)');
    ylabel('strain (10^{-3})');
    title('case 1 strain: end load');
    legend({'total', 'thermal', 'elastic'}, 'Location', 'southoutside', ...
        'Orientation', 'horizontal', 'FontSize', 11);
    set(gca, 'FontSize', 14, 'LineWidth', 1.0);

    nexttile;
    bar(freeThermal.xElem, [freeThermal.epsTotal, freeThermal.epsThermal, freeThermal.epsElastic] * 1000, 'grouped');
    grid on; box on;
    xlabel('element center x (mm)');
    ylabel('strain (10^{-3})');
    title('case 2 strain: free thermal');
    legend({'total', 'thermal', 'elastic'}, 'Location', 'southoutside', ...
        'Orientation', 'horizontal', 'FontSize', 11);
    set(gca, 'FontSize', 14, 'LineWidth', 1.0);

    nexttile;
    bar(fixedThermal.xElem, [fixedThermal.epsTotal, fixedThermal.epsThermal, fixedThermal.epsElastic] * 1000, 'grouped');
    grid on; box on;
    xlabel('element center x (mm)');
    ylabel('strain (10^{-3})');
    title(sprintf('case 3 strain: fixed thermal, \\DeltaT = %d K', dTUniform));
    legend({'total', 'thermal', 'elastic'}, 'Location', 'southoutside', ...
        'Orientation', 'horizontal', 'FontSize', 11);
    set(gca, 'FontSize', 14, 'LineWidth', 1.0);

    hide_axes_toolbars(fig);
    exportgraphics(fig, 'fem_1d_result_summary.png', 'Resolution', 220);
    close(fig);
end

function make_deformation_animation(x, loadCase, freeThermal, fixedThermal, fileName, frameRate)
    % 最終変位までを線形補間し、変形の比較アニメーションを保存します。
    fig = figure('Name', '1D FEM deformation animation', 'Color', 'w', 'Position', [120 120 1200 620]);
    writer = VideoWriter(fileName, 'MPEG-4');
    writer.FrameRate = frameRate;
    open(writer);

    cases = {loadCase, freeThermal, fixedThermal};
    names = {'end load', 'free thermal expansion', 'fixed thermal expansion'};
    colors = [0.1 0.3 0.75; 0.1 0.55 0.35; 0.85 0.25 0.15];
    scale = 140;  % 小さな変位を見やすくする表示倍率です。

    nFrame = 75;
    nRamp = 45;
    for frame = 1:nFrame
        ratio = min(1, (frame - 1) / (nRamp - 1));
        clf(fig);
        tiledlayout(fig, 3, 1, 'Padding', 'compact', 'TileSpacing', 'compact');
        for c = 1:3
            nexttile;
            hold on;
            y = zeros(size(x));
            plot(x, y, 'Color', [0.72 0.72 0.72], 'LineWidth', 6);
            xd = x + scale * ratio * cases{c}.u;
            plot(xd, y, '-o', 'Color', colors(c, :), 'MarkerFaceColor', 'w', ...
                'LineWidth', 4, 'MarkerSize', 7);
            title(sprintf('%s: final deformation interpolation, scale x%d', names{c}, scale));
            xlim([-5, max(x + scale * freeThermal.u) + 5]);
            ylim([-0.35 0.35]);
            yticks([]);
            grid on; box on;
            set(gca, 'FontSize', 15, 'LineWidth', 1.0);
        end
        drawnow;
        writeVideo(writer, getframe(fig));
    end

    close(writer);
    close(fig);
end

function make_assembly_animation(x, conn, E, A, fileName, frameRate)
    % 要素剛性行列が全体剛性行列へ加算される過程を保存します。
    nNode = numel(x);
    K = zeros(nNode, nNode);
    fig = figure('Name', '1D FEM stiffness assembly animation', 'Color', 'w', 'Position', [120 120 760 680]);
    writer = VideoWriter(fileName, 'MPEG-4');
    writer.FrameRate = frameRate;
    open(writer);

    for e = 1:size(conn, 1)
        nodes = conn(e, :);
        Le = x(nodes(2)) - x(nodes(1));
        ke = E * A / Le * [1 -1; -1 1];
        K(nodes, nodes) = K(nodes, nodes) + ke;

        clf(fig);
        imagesc(K);
        axis equal tight;
        colormap(parula);
        colorbar;
        title(sprintf('assembling element %d into the global stiffness matrix', e));
        xlabel('node number');
        ylabel('node number');
        set(gca, 'FontSize', 17, 'LineWidth', 1.0);
        drawnow;
        for repeatFrame = 1:8
            writeVideo(writer, getframe(fig));
        end
    end

    close(writer);
    close(fig);
end

function hide_axes_toolbars(fig)
    % exportgraphicsに軸ツールバーが写り込まないようにします。
    axesHandles = findall(fig, 'Type', 'axes');
    for i = 1:numel(axesHandles)
        axesHandles(i).Toolbar.Visible = 'off';
    end
end

　定式化では支配方程式から弱形式を導きましたが、MATLABコードは、すでに得られた要素方程式\((27)\)を数値的に組み立てるところから始まります。 まず節点座標と要素接続を作ります。 その後、各要素について式\((31)\)の要素剛性行列を作り、式\((37)\)に従って全体剛性行列へ足し込みます。 温度がある場合は、式\((33)\)の等価熱荷重 \(\boldsymbol{f}^{\mathrm{th}}_e\) も同じ加算規則で足し込みます。 境界条件を入れて節点変位を解いた後、各要素で式\((16)\)により全ひずみを計算し、式\((19)\)により応力を計算します。

　コードの流れを式で書けば、

　式\((43)\)は、新しい支配方程式ではありません。 FEMコードがどの順番で物理量を計算しているかを表したものです。 熱弾塑性へ進むと、最後の \(\varepsilon_e,\sigma_e\) を求める部分で、各要素または各積分点の塑性ひずみや内部変数も更新することになります。

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

　1次元棒要素FEMでは、未知量は節点変位です。 応力を直接未知量として解くのではなく、まず節点変位を解き、そこから要素ひずみと応力を計算します。 この構造を理解しておくと、熱弾塑性解析で「どこで材料モデルを呼ぶのか」が見えやすくなります。

　また、熱ひずみは単に温度から計算される自由膨張成分です。 応力になるかどうかは、境界条件によって決まります。 自由に伸びられる要素では熱ひずみが全ひずみとして現れ、応力はゼロになります。 拘束される要素では、熱ひずみを打ち消す弾性ひずみが発生し、その結果として熱応力が生じます。 FEMは、この関係を要素ごとに評価し、全体の釣合いを満たす節点変位を求める方法です。

塑性ひずみの考慮

　この記事では、温度を要素内で一様とし、材料も線形弾性のまま扱いました。 次の記事では、まず温度条件は一様なまま、各要素の積分点へ塑性ひずみを導入します。 応力式は、

　式\((44)\)のようになり、各積分点で降伏判定、塑性補正、塑性ひずみの更新を行います。 温度を一様なままにすることで、FEMへの材料モデルの組み込み方だけに集中できます。

まとめ

　この記事では、1次元棒要素FEMの基本を、熱応力解析へつながる形で整理しました。 最初に棒の力の釣合いから支配方程式を定め、試験関数を掛けて部分積分することで弱形式へ変換しました。 次に、ガラーキン法として試験関数を全体基底関数へ一致させ、要素内では基底関数を形状関数として表しました。 そのうえで、弱形式から一要素の剛性行列と等価熱荷重を導き、それらを全体方程式へ組み立てる流れを確認しました。 また、熱ひずみを等価節点力として扱い、自由熱膨張と両端固定加熱の違いを数値例で確認しました。

　重要なのは、FEMでも熱応力の本質は変わらないという点です。 熱ひずみが発生しても、それが自由に実現できれば応力は生じません。 一方で、境界条件や周囲の要素によって熱膨張が妨げられると、弾性ひずみが生じ、応力が発生します。 FEMは、この関係を要素ごとに計算し、全体として釣り合う変位場を求めるための枠組みです。

　次の記事では、この1次元棒要素FEMへ熱弾塑性の応力更新を組み込みます。