同軸上に配置された２つの円筒を均一にある温度まで加熱された時の円筒の伸びを計算します。
この問題は、不静定問題といわれており、材料力学の教科書や技術士一次試験（機械問題）の問題にも出題されています。
＜動作確認環境＞
Salome-Meca 2019 (CAELinux2020Lite)
１．解析対象
下図に解析モデルを示します。
１．１ モデル寸法・物性値
・円筒1：内径7.5mm x 外径12.5ｍｍ x 高さ40ｍｍ
・円筒2：内径17.5mm x 外径22.5ｍｍ x 高さ40ｍｍ
・円盤：外径30、高さ2.5ｍｍの剛体
・円柱1 、円柱2 、円盤の温度：20℃（初期）→100℃
表 物性値一覧
| モデル | ヤング率(MPa) | ポアソン比 | 熱膨張係数(1/K) |
| 円筒１ | 200000 | 0.3 | 1.1x10-5 |
| 円筒２ | 70000 | 0.3 | 2.3x10-5 |
| 円盤(*) | 200000000 | 0 | 0 |
（*)剛体指定が不明だったため、ヤング率を円筒1に対して1000倍大きくして、事実上の変形による影響をなくした
１．２ 解析モデル
・軸対称モデル、Y方向（高さ方向）も対称
・メッシュ：四角形１次要素
・拘束条件：円柱１、円柱２の底面にX,Y方向固定、円盤をX方向固定
・解析：熱応力、拘束解析
・出力：円盤と円柱との接触位置の伸び（Y方向変形）
１．３ 解析手順
１．３．１ 連成解析
解析ファイル（commファイル）を伝熱解析用、構造解析用の2種類を作り解析を行います。
１）伝熱解析にて、解析モデル全体を100℃に設定する。
２）構造解析にて、20℃→100℃になったときの物体の変形分だけ熱応力が発生する。
１．３．２ 熱解析ウィザードによる解析
SalomeMecaによるのthermo-mechanical analysisを使用しました。
thermo-mechanical analysis は、ウィザード方式の設定であり、メッシュの選択、初期温度の入力、時間の入力、FEMモデルの設定（軸対称、3Dモデルなど）、荷重などの拘束条件入力の順番で計算条件を入力します。
引用文献：Thermal analysys, code_aster,Salome_meca course material P.19以降を参照ください。
URL:https://www.code-aster.org/UPLOAD/DOC/Formations/06-thermal_analysis.pdf
２．解析結果
どちらの手法でも同様の解析結果でした。
円筒１、円筒２の変位をみるとほぼ同値ですので、 円盤のヤング率を極端に大きくすることで剛体を近似できていました。
２．１ 連成解析による結果
Y方向変位は、それぞれ円筒１：0.269mm、円筒2：0.268mmです。
本解析では対称の考え方から変位は半分の伸びをあらわしていますので、伸びは、変位を2倍して、円筒１：0.538mm、円筒2：0.536mmとなります。
２．２ 熱解析ウィザードによる解析結果
伸びは、それぞれ円筒１：0.269mm、円筒2：0.268mmで、連成解析と同値でした。
こちらも 連成解析と同様に対象の考え方から半分の伸びをあらわしていますので2倍して、 伸びは、変位を2倍して、円筒１：0.538mm、円筒2：0.536mmとなります。
３．検算
３．１ 材料力学の公式との比較
今回の熱応力の場合、円筒の伸びλは、材料力学の公式より、0.0504mmです。
E1 : 円筒1のヤング率(N/mm2)、α1 : 円筒1の熱膨張係数(1/K) 、A1 : 円筒1の断面積(mm2)
E2 : 円筒2のヤング率(N/mm2) 、 α2 : 円筒2の熱膨張係数(1/K) 、A2 : 円筒2の断面積(mm2)
Δt : 温度差(K) 、 L:円筒の長さ (mm)
参照資料：材料力学（改訂版）P.26,P.27,P.28より 黒木剛司郎著（森北出版）
材料力学との比較をすると、5%程度のずれがありました。比較結果を下表に示します。
表 SalomeMecaの計算値と材料力学との比較結果
| SalomeMecaの計算値 | 計算結果 | 材料力学との差異（％） | |
| 円筒１ | 0.0537 | 0.0504 | 5.4％ |
| 円筒２ | 0.0538 | 0.0504 | 5.2% |
３．２ 材料力学の公式との比較結果を考察する
SalomeMecaと材料力学の理論式との乖離が大きいので、円筒１と円筒２の材料をアルミに変更してY方向 の変位（軸方向の伸び）を比較しました。 伸びλは 熱膨張の公式より以下で表します。
伸びλ ＝ α2 x Δt x L ＝2.3x10-5x （100-20） x 40 = 0.736 mm(半分 0.368mm)
α2 : 円筒2の熱膨張係数(1/K) 、Δt : 温度差(K) 、 L:円筒の長さ (mm)
結果を見ていくと5%程度大きく、剛体の拘束をうけたものと考えられます。
４．まとめ
熱変形解析は、連成解析でしかできなかったのですが、熱応力解析ウィザードによる解析も可能です。使いこなせると便利だと考えます。
実務で解析する場合、物体には温度分布があることが多いため、理論値がモデルで計算を行いFEMモデルと理論値の差異を確認してから、実機の誤差を予測しておくのが望ましいです。
本内容を5月上に開催された第21回オープンCAE勉強会＠関東（構造など）で発表しました。指摘があったのは、LIEAON_UNIFというコマンドを使うと、剛体拘束を再現できるとのことでした。
５．補足
(1) 連成解析
伝熱解析
円柱１、円柱２、円盤 を100℃に設定します。
DEBUT(identifier='0:1',
LANG='EN')
mesh = LIRE_MAILLAGE(identifier='1:1',
FORMAT='MED',
UNITE=3)
model = AFFE_MODELE(identifier='2:1',
AFFE=_F(MODELISATION=('AXIS', ),
PHENOMENE='THERMIQUE',
TOUT='OUI'),
MAILLAGE=mesh)
mater0 = DEFI_MATERIAU(identifier='3:1',
THER=_F(LAMBDA=0.016))
fieldma0 = AFFE_MATERIAU(identifier='4:1',
AFFE=_F(MATER=(mater0, ),
TOUT='OUI'),
MODELE=model)
temp = AFFE_CHAR_THER(identifier='5:1',
MODELE=model,
TEMP_IMPO=_F(TEMP=100.0,
TOUT='OUI'))
thload = THER_LINEAIRE(identifier='6:1',
CHAM_MATER=fieldma0,
EXCIT=_F(CHARGE=temp),
MODELE=model)
IMPR_RESU(identifier='7:1',
FORMAT='MED',
RESU=_F(RESULTAT=thload),
UNITE=4)
FIN(identifier='8:1',
)
構造解析
POORSUITEを用いることで、伝熱解析の結果（円柱１、円柱２、円盤全体が100℃に加熱された結果)を入力します。
POURSUITE(identifier='0:2',
LANG='EN')
model1 = AFFE_MODELE(identifier='1:2',
AFFE=_F(MODELISATION=('AXIS', ),
PHENOMENE='MECANIQUE',
TOUT='OUI'),
MAILLAGE=mesh)
mater1 = DEFI_MATERIAU(identifier='2:2',
ELAS=_F(ALPHA=1.1e-05,
E=200000.0,
NU=0.3))
mater3 = DEFI_MATERIAU(identifier='3:2',
ELAS=_F(ALPHA=0.0,
E=200000000.0,
NU=0.3))
mater2 = DEFI_MATERIAU(identifier='4:2',
ELAS=_F(ALPHA=2.3e-05,
E=70000.0,
NU=0.3))
fieldma1 = AFFE_MATERIAU(identifier='5:2',
AFFE=(_F(GROUP_MA=('innerring_Faces', ),
MATER=(mater1, )),
_F(GROUP_MA=('top_Faces', ),
MATER=(mater3, )),
_F(GROUP_MA=('outering_Faces', ),
MATER=(mater2, ))),
AFFE_VARC=_F(EVOL=thload,
NOM_VARC='TEMP',
TOUT='OUI',
VALE_REF=20.0),
MODELE=model1)
fix = AFFE_CHAR_MECA(identifier='6:2',
DDL_IMPO=(_F(DY=0.0,
GROUP_NO=('outer_base', 'inner_base')),
_F(DX=0.0,
GROUP_NO=('top_yslide', ))),
MODELE=model1)
reslin = MECA_STATIQUE(identifier='7:2',
CHAM_MATER=fieldma1,
EXCIT=_F(CHARGE=fix),
MODELE=model1)
reslin = CALC_CHAMP(identifier='8:2',
reuse=reslin,
CONTRAINTE=('SIGM_NOEU', 'SIEF_NOEU'),
CRITERES=('SIEQ_NOEU', ),
DEFORMATION=('EPSI_NOEU', ),
FORCE=('FORC_NODA', 'REAC_NODA'),
RESULTAT=reslin)
IMPR_RESU(identifier='9:2',
FORMAT='MED',
RESU=_F(NOM_CHAM=('EPSI_NOEU', 'SIEQ_NOEU', 'SIGM_NOEU', 'DEPL', 'SIEF_NOEU'),
RESULTAT=reslin),
UNITE=80)
FIN(identifier='10:2',
)
(2) 熱応力ウィザードによる解析
1つの解析ファイルで解析を行っていきます。 LAMBDA,CHO_CPは計算には使わないので、任意の値を入力しています。こちらの計算は、一部計算に不具合があったので、線膨張係数は固定値だったため、温度差を80℃を反映すべく、基準の温度を0℃(VALE=0)、温度上昇分を80℃（TEMP=80）にて計算をしています。
DEBUT(LANG='EN')
mesh = LIRE_MAILLAGE(identifier='0:1',
FORMAT='MED',
UNITE=20)
modelth = AFFE_MODELE(identifier='1:1',
AFFE=_F(MODELISATION=('AXIS', ),
PHENOMENE='THERMIQUE',
TOUT='OUI'),
MAILLAGE=mesh)
modmeca = AFFE_MODELE(identifier='2:1',
AFFE=_F(MODELISATION=('AXIS', ),
PHENOMENE='MECANIQUE',
TOUT='OUI'),
MAILLAGE=mesh)
mater1 = DEFI_MATERIAU(identifier='3:1',
ELAS=_F(ALPHA=2.3e-05,
E=70000.0,
NU=0.3),
THER=_F(LAMBDA=0.5,
RHO_CP=1.0))
mater2 = DEFI_MATERIAU(identifier='4:1',
ELAS=_F(ALPHA=2.3e-05,
E=70000.0,
NU=0.3),
THER=_F(LAMBDA=0.5,
RHO_CP=1.0))
mater3 = DEFI_MATERIAU(identifier='5:1',
ELAS=_F(ALPHA=2.3e-05,
E=70000.0,
NU=0.3),
THER=_F(LAMBDA=0.5,
RHO_CP=1.0))
matth = AFFE_MATERIAU(identifier='6:1',
AFFE=(_F(GROUP_MA=('innerring_Faces', ),
MATER=(mater1, )),
_F(GROUP_MA=('outering_Faces', ),
MATER=(mater2, )),
_F(GROUP_MA=('top_Faces', ),
MATER=(mater3, ))),
MAILLAGE=mesh)
linst = DEFI_LIST_REEL(identifier='7:1',
VALE=(0.0, 1.0))
loadst = AFFE_CHAR_THER(identifier='8:1',
MODELE=modelth,
TEMP_IMPO=_F(GROUP_MA=('innerring_Faces', 'top_Faces', 'outering_Faces'),
TEMP=80.0))
mecabc = AFFE_CHAR_MECA(identifier='9:1',
DDL_IMPO=(_F(DX=0.0,
DY=0.0,
GROUP_NO=('outer_base', 'inner_base')),
_F(DX=0.0,
GROUP_NO=('top_yslide', ))),
MODELE=modmeca)
resuth = THER_LINEAIRE(identifier='10:1',
CHAM_MATER=matth,
ETAT_INIT=_F(VALE=0.0),
EXCIT=_F(CHARGE=loadst),
INCREMENT=_F(LIST_INST=linst),
MODELE=modelth)
matmeca = AFFE_MATERIAU(identifier='11:1',
AFFE=(_F(GROUP_MA=('innerring_Faces', ),
MATER=(mater1, )),
_F(GROUP_MA=('outering_Faces', ),
MATER=(mater2, )),
_F(GROUP_MA=('top_Faces', ),
MATER=(mater3, ))),
AFFE_VARC=_F(EVOL=resuth,
NOM_VARC='TEMP',
TOUT='OUI',
VALE_REF=0.0),
MODELE=modmeca)
resumeca = MECA_STATIQUE(identifier='12:1',
CHAM_MATER=matmeca,
EXCIT=_F(CHARGE=mecabc),
LIST_INST=linst,
MODELE=modmeca)
IMPR_RESU(identifier='13:1',
FORMAT='MED',
RESU=(_F(RESULTAT=resuth),
_F(RESULTAT=resumeca)),
UNITE=80)
FIN()