3. 熱・力学連成解析を用いた評価
3.1 熱・力学連成解析プログラムの開発
式3-2および式3-3を増分系として整理する。
Δ𝑥 𝑥 𝛿 Δ𝑥 Δ𝑡 3 4
Δ𝑥 1
𝛽Δ𝑡 Δ𝑥 𝑥 Δ𝑡 Δ𝑡
2 𝑥 3 5
式3-4および式3-5を式3-1に代入して応答値を算出する。ここで、𝐾 はLRBの履歴則に基づ き逐次変化するため、Newton法を用いた収斂計算を行う。発生した残差ベクトル 𝛾 は、負担力 ベクトル 𝐸 として次ステップに持ち越している。
式3-1から算出した鉛プラグの吸収エネルギー量を、式3-6を用いて単位体積当たりの発熱エ ネルギーに変換する。
𝑄 𝑊
𝐶 ∙ 𝑉 3 6
ここで、
𝑄：単位体積・単位時間当たりの発熱量（kW/m3） 𝑊：鉛プラグの吸収エネルギー量（kW・m）
𝐶：鉛プラグの熱容量（kW・m/K）
𝑉：発熱範囲の鉛プラグ体積（m3）
次に、式3-6で算出した熱エネルギーを、式3-7に示す熱伝導方程式2)に入力し、LRBの温度 変化を算出する。
𝜌𝑐𝜕𝑇
𝜕𝑡 ∇ 𝜆∇𝑇 𝑄 3 7
ここで、
𝑇：要素内温度（℃）
𝜌：要素密度（ton/m3） 𝑐：要素比熱（kJ/ton・K）
𝜆：要素熱伝導率（kW/m・K）
本プログラムでは、LRBと外部環境の境界条件として、ノイマン境界（式3-8）とディリクレ
境界（式3-9）を組込むこととした。これは、2章で示した繰返し加力試験結果でも確認できる
通り、LRBは鉛プラグの周辺を積層ゴムで被覆されており、装置外周部の温度上昇は殆ど見ら れないことから、対流を考慮せず、外気温一定条件下で熱伝達を再現したノイマン境界を選択し た。また、LRBの上下端には、繰返し加力試験では断熱材が設置されており、実建物では一般
的にコンクリート造のペデスタルが設置されている。よって、断熱材もしくはペデスタルから外 部に伝導する熱量は半径方向に伝導する熱量より小さいと考えられることから、境界上の温度を 予め規定できるディリクレ境界を選択した。
𝑞 𝛼 𝑇 𝑇 3 8
𝑇 𝑇 3 9
ここで、
𝛼 ：外部熱伝導率（kW/m・K）
𝑇：外部温度（℃）
𝑇：規定温度（℃）
式3-7に要素内温度の変分𝛿𝑇を乗じ、要素領域𝑉全体で積分を行う。
𝜌𝑐 𝛿𝑇𝑇𝑑𝑉 𝛿𝑇 ∇ 𝜆∇𝑇 𝑑𝑉 𝛿𝑇 ∙ 𝑄 𝑑𝑉
∇ 𝛿𝑇 𝜆∇𝑇 𝑑𝑉 ∇𝛿𝑇 𝜆∇𝑇 𝑑𝑉 𝛿𝑇 ∙ 𝑄 𝑑𝑉
次に、式3-10の右辺第2項を移項し、右辺第1項に発散定理3)を適用する。なお、𝑑𝑆は境界面 の微小面積を示す。右辺第1項は、要素領域𝑉の境界面法線方向熱流速となるため、フーリエ則
（𝑞 𝑛 ∙ ∇ 𝜆𝑇 ）から式3-11となる。
𝜌𝑐𝛿𝑇𝑇𝑑𝑉 ∇𝛿𝑇 𝜆∇𝑇 𝑑𝑉 𝑛 ∙ 𝛿𝑇 𝜆∇𝑇 𝑑𝑆 𝛿𝑇 ∙ 𝑄 𝑑𝑉
𝛿𝑇 ∙ 𝑄 𝑑𝑉 𝛿𝑇 ∙ 𝑞 𝑑𝑆
既往研究4)では、差分法により熱伝導解析を行う熱・力学連成解析手法が開発されている。し かし、差分法は要素サイズにより解析刻みに制限を受けるため、極小の接触要素をモデル化する 本研究の提案手法には適していない。そこで、本プログラムでは、有限要素法を用いて熱伝導解 析を行う。まず、式3-11を要素内積分の重ね合わせで表現する。なお、右辺第1項は発熱量ベ クトルを表しているため、要素の重ね合わせは鉛プラグの発熱範囲のみである。
𝜌 𝑐𝛿𝑇𝑇𝑑𝑉 ∇𝛿𝑇 𝜆∇𝑇 𝑑𝑉 𝛿𝑇 ∙ 𝑄 𝑑𝑉 𝛿𝑇 ∙ 𝑞 𝑑𝑆 3 12
次に、要素内温度を節点温度と内挿関数𝑁（境界上の内挿関数は𝑁）に置換する。ここで、要素 3 10
3 11
構成節点番号を𝑝,𝑞としている。
𝑇 𝑁 𝑇 3 13
𝛿𝑇 𝑁 𝛿𝑇 3 14
𝛿𝑇 𝑁 𝛿𝑇 3 15
式3-13～式3-15を式3-12に代入して整理する。
𝜌𝑐𝑁 𝛿𝑇 𝑁 𝑇 𝑑𝑉 ∇𝑁 𝛿𝑇 𝜆∇𝑁 𝑇 𝑑𝑉
𝑁 𝛿𝑇 ∙ 𝑄 𝑑𝑉 𝑁 𝛿𝑇 ∙ 𝑞 𝑑𝑆
𝛿𝑇 𝜌𝑐𝑁 𝑁 𝑑𝑉 𝑇 𝛿𝑇 ∇𝑁 λ∇𝑁 𝑑𝑉 𝑇
𝛿𝑇 𝑁 ∙ 𝑄 𝑑𝑉 𝛿𝑇 𝑁 ∙ 𝑞 𝑑𝑆
要素構成節点番号𝑝,𝑞と領域全体での節点番号𝑎,𝑏の対応付けができるものとすると、式3-17は 以下のように整理できる。
𝛿𝑇 𝜌𝑐𝑁 𝑁 𝑑𝑉 𝑇 𝛿𝑇 ∇𝑁 𝜆∇𝑁 𝑑𝑉 𝑇
𝛿𝑇 𝑁 ∙ 𝑄 𝑑𝑉 𝛿𝑇 𝑁 ∙ 𝑞 𝑑𝑆
𝜌𝑐𝑁 𝑁 𝑑𝑉 𝑇 ∇𝑁 𝜆∇𝑁 𝑑𝑉 𝑇
𝑁 ∙ 𝑄 𝑑𝑉 𝑁 ∙ 𝑞 𝑑𝑆
3 16
3 17
3 18
3 19
本研究では、円筒形のLRBを対象にモデル化を行うため、3次元直交座標系である式を円筒 座標系に置き換える5)。直交座標系と円筒座標系の関係は式3-20の通りである。ここで、𝑟
（𝑟 𝑥 𝑦 ）は半径方向、𝜃は円周方向の自由度を示す。
𝑥 𝑦 𝑧
𝑟𝑐𝑜𝑠𝜃 𝑟𝑠𝑖𝑛𝜃
𝑧
3 20
次に、領域積分と境界積分を円筒座標系に置き換える。領域積分を式3-21に示す。
𝑑𝑉 𝑑𝑥𝑑𝑦𝑑𝑧 𝜕 𝑥 𝑦 𝑧
𝜕 𝑟 𝜃 𝑧 𝑑𝑟𝑑𝜃𝑑𝑧 𝑑𝑒𝑡 𝑐𝑜𝑠𝜃 𝑟𝑠𝑖𝑛𝜃 0
𝑠𝑖𝑛𝜃 𝑟𝑐𝑜𝑠𝜃 0
0 0 1
𝑑𝑟𝑑𝜃𝑑𝑧 𝑟 ∙ 𝑑𝑟𝑑𝜃𝑑𝑧
境界積分6)を式3-22に示す。なお、𝑑𝑠は境界面の線素（境界面微小面積における1辺の長さ）
を示す。
𝑑𝑆 𝑟 ∙ 𝑑𝑠𝑑𝜃 1 2𝑑𝑠 𝑑𝜃
≒ 𝑟 ∙ 𝑑𝑠𝑑𝜃
式3-21および式3-22から、式3-19の非定常熱伝導方程式を円筒座標系に変換する。
𝜌𝑐𝑁 𝑁 𝑟 𝑑𝑧𝑑𝜃𝑑𝑟 𝑇 ∇𝑁 λ∇𝑁 𝑟 𝑑𝑧𝑑𝜃𝑑𝑟 𝑇
𝑁 ∙ 𝑄 ∙ 𝑟 𝑑𝑧𝑑𝜃𝑑𝑟 𝑁 ∙ 𝑞 ∙ 𝑟 𝑑𝑠𝑑𝜃
LRBは軸対称の形状であるため、𝜃方向一定条件を加えると式3-24に整理できる。ここで、左 辺第1項は熱容量マトリクス、第2項は熱伝導マトリクス、右辺第1項は発熱量ベクトル、第2 項は熱流束境界条件を示す。
3 21
3 22
3 23
𝜌𝑐𝑁 𝑁 𝑟 𝑑𝑧𝑑𝑟 𝑇 ∇𝑁 λ∇𝑁 𝑟 𝑑𝑧𝑑𝑟 𝑇
𝑁 ∙ 𝑄 ∙ 𝑟 𝑑𝑧𝑑𝑟 𝑁 ∙ 𝑞 ∙ 𝑟 𝑑𝑠
本研究では、非線形熱伝導方程式の解法として後退差分法を用いることとする。ここでは、簡 単のためにマトリクス表示とし、未知量である温度ベクトル 𝑇 を左辺に移行して整理する。
𝑀 𝑇 𝐾 𝑇 𝑄 𝑓 𝑓 𝐾 𝑇 3 25
𝑀 𝑇 𝐾 𝑇 𝐹 3 26
ここで、
𝑀 ：要素熱容量マトリクス 𝐾 ：要素熱伝導マトリクス 𝐾 ：境界要素熱伝導マトリクス
𝐹 ：発熱量と熱流速による境界面からの放熱量ベクトル
時間刻みを𝑑𝑡として後退差分近似する。
𝑀 𝑇 𝑇
𝑑𝑡 𝐾 𝑇 𝐹 3 27
左辺に未知量である 𝑇 を移項すると以下の連立1次方程式を得る。
1
𝑑𝑡 𝑀 𝐾 𝑇 𝐹 1
𝑑𝑡 𝑀 𝑇 3 28
𝑀 𝐾 は実対称行列であるため、LDLT分解することができる。7)そこで、連立1次方程 式の解は、以下に示す前進消去（式3-31）と後退代入（式3-32）により算出する。
𝐴 𝑥 𝑦 3 29
3 24
𝐿 𝐷 𝐿 𝑥 𝑦 3 30
𝐿 𝑏 𝑦 3 31
𝐷 𝐿 𝑥 𝑏 3 32
ここで、
𝐴 1
𝑑𝑡 𝑀 𝐾 𝐿 𝐷 𝐿 3 33
𝑥 𝑇 3 34
𝑦 𝐹 1
𝑑𝑡 𝑀 𝑇 3 35
なお、 𝐷 は対角行列、𝐿 は下三角行列を示す。行列 𝐴が実対称行列であれば、 𝐴 𝐿 𝐷 𝐿 と表すことができ、LDLT分解が可能となる。
式3-28から算出した鉛プラグの温度分布から、発熱範囲の体積平均温度を算出し、式2-5に 示す既往提案式に代入することで、次ステップの地震応答解析で使用する鉛プラグの降伏応力度 を算出する。
𝑇 ∑ 𝑇 ∙ 𝑉
𝑉 3 36
ここで、
𝑇 ：鉛プラグの体積平均温度（℃）
𝑇 ：発熱範囲の節点温度（℃）
𝑉：節点𝑛の支配体積（m3）
本解析プログラムでは、式3-1の運動方程式における残差力が基準値以下となるまで収斂計算 を行い、各積分ステップの地震応答解析および熱伝導解析の解を算出する。熱・力学連成解析の フローチャートを図3-1に示す。
図3-1 熱・力学連成解析のフローチャート 熱・力学連成解析開始
時間ステップ開始
時間ステップ終了 地震応答解析(1ステップ)
発熱量算出（式3-6）
熱伝導解析(1ステップ)
鉛プラグの体積平均温度算出（式3-36）
鉛プラグの降伏応力度算出（式2-5）
収斂計算開始
履歴追跡
収斂計算終了
熱・力学連成解析終了