SENAC Vol. 54, No. （2 2021. 4）― 1 ―［共同研究成果］

数値解の安定条件を考慮した y 形溶接構造の冷却過程における水素凝集挙動解析および本解析プログラムによる数値解析
高速化へ向けての考察
尾関郷 1，横堀壽光 1，大見敏仁 2，小川道夫 1，糟谷正 3
1
帝京大学 戦略的イノベーション研究センター2
湘南工科大学 機械工学科
3
東京大学大学院 工学系研究科
溶接部における強度安全性維持における重要な問題として、低温割れが挙げられる。低温割れは溶接部冷却後に発生する残留応力と、溶接金属中に存在する水素の凝集により発現する水素脆化が原因とされている。特に金属中の水素拡散・凝集挙動の時系列変化を把握することが重要であるが、水素拡散の様子を実験的に観察することは非常に困難であることから、数値計算による予測が工学的に重要である。本研究では、著者らが新たに開発した熱伝導－熱応力－局所応力誘起水素拡散連成解析を用いて、予熱を含んだ y 形溶接部の冷却過程における水素拡散・凝集挙動を解析するとともに、数値解析の安定条件についての考察を行った。また、長時間解析に適応するため、本解析プログラムの高速化についても検討した。
記号説明
Length of analytical model
l0C
Hydrogen concentration
Hydrogenconcentration
inC0atmosphere
DDiffusion coefficient(= D0 exp(-Q/RT))
Diffusion constant independent ofD0temperature
T
Absolute temperature
Tout
Temperature in atmosphere
QActivation energyR
Gas constant
VolumechangeduetoVaccommodation of a hydrogen
Hydrostatic stresspE
Young's modulus
Poisson’s ratioYield stressys1.

ysEquivalent stress
Yield stressHp

Work hardening coefficientp

Equivalent plastic strainc1, , m1Constants

htcth

Hydrogen transfer coefficientTime of hydrogen diffusion
Time of heat transferDensity

cSpecific heatkHhTTHeat conductivity
Heat transfer coefficient=H/kThermal expansion coefficient

緒言
構造物の安全性を確保する上で、溶接部の信頼性確保は非常に重要である。溶接部における強度安全性維持における重要な問題として、低温割れが挙げられる。低温割れは溶接部冷却後に発生する残留応力と、溶接金属中に存在する水素の凝集により発現する水素脆化が原因とされている [1]。低温割れを防ぐ為に、母材を予め加熱してから溶接を行う予熱処理が施される。これは鋼中の水素拡散を促進し、溶接部の水素量を減らすことが狙いである。しかし、溶接部冷却過程における水素凝集挙動には、鋼種、応力集中、予熱温度など様々な要因が複雑に影響するため適切な予熱条件を設定するのが困難である。特に金属中の水素拡散・凝集挙動の時系列変化を把握することが重要であるが、水素拡散の様子を実験的に観察することは非常に困難であることか

― 2 ―

SENAC Vol. 54, No. （2 2021. 4）

ら、数値計算による予測が工学的に重要である。
著者らはこれまで、マルティプリケーション法[2,3] を組み入れた熱伝導－熱応力－局所応力
水素拡散連成解析プログラムを開発し、y 形溶接構造を対象として溶接部の冷却時における水素拡散・凝集挙動の解析を行ってきた[4~6,8]。しかし、溶接部の冷却のような急激な温度変化に起因する応力変動を伴う解析の場合、拡散方程式の一次微分項において振動解を誘起する可能性があるため、安定な数値解を得られる条件についても考慮する必要がある[6, 7]。著者の一人は、この影響を除くために Δt / (ΔL)2 < λc (t:時間、L:位置、λc:収束度限界パラメータ)なる条件を導入した[2,3]。
本研究では、これまで開発してきた熱伝導－熱応力－局所応力水素拡散連成解析プログラムを用いて、様々な予熱温度における y 形溶接構造の水素輸送挙動を解析するとともに、数値解析安定条件のために、Δt / (ΔL)2 < λc の効果を検討した。また、長時間解析に適応するため、本解析プログラムの高速化についても検討した。

2.

熱伝導-熱熱応力-水
水素拡散連成解析 [4~6,8]
溶接部における冷却時の水素拡散挙動を解析するため、著者らは熱伝導、熱応力、応力誘起水素拡散を連成した解析プログラムを開発してきた。本解析は、応力解析を有限要素法 (FEM)、熱伝導と拡散解析を有限差分法 (FDM)で行う、FEM-FDM 法[2,9]を用いている。以下に解析手法を示す。まず熱伝導を FDM で計算し、その後、FEM の各節点へ温度を補間する。次に、補間された温度を用いて各節点の熱応力を計算する。その後、得られた熱応力を差分格子点へ補間し、応力誘起水素拡散解析を行う。これらの計算を連続的に繰返すことで冷却時における水素拡散解析を実施した。本解析の流れを Fig. 1 に示す[4,5]。

Fig. 1 Flowchart of coupled analysis of heat conduction-thermal stress-hydrogen diffusion [4,5]2.1 解析モデルおよび境界条件
解析モデルは JIS 規格[10]に基づいて、溶金 (WM)、母材(BM)、熱影響部（HAZ: Heat AffectedZone）から成る y 形溶接構造を作成した。解析モデルを Fig. 2 に示す。本解析は 2 次元解析で行い、FEM の総節点数は 6,144、総要素数は 11,713 であり、FDM の総格子点数は 1,500 である。応力解析に用いた物性値を Table 1 に示す。ヤング率および降伏応力は Fig. 3、Fig. 4 に示すように、温度に対応するようにしている [11]。高温からの冷却を考え、降伏応力はマルテンサイトの物性値を用いた。また、降伏応力は式(1)の条件にて解析を行った [12]。

数値解の安定条件を考慮した y 形溶接構造の冷却過程における水素凝集挙動解析および本解析プログラムによる数値解析高速化へ向けての考察 ys (WM , BM )  0 .8 ys ( HAZ )― 3 ―(1).

線膨張係数は、一般的な鉄の値を用いた。ヤング率、ポアソン比、線膨張係数は WM や BM といった組織による違いはないものとし、応力解析は平面ひずみ状態にて行った。水素拡散および熱伝導解析の境界条件は Fig. 5 に示すとおりである。

Table 1 Material properties for FEM analysis.
E[GPa]



ref. [Fig. 3]

0.3

 ys

p

[MPa]

[GPa]

T
[1/K]

ref. [Fig. 4]

0.01×E

1.2E-05

WM
HAZ
BMFig. 2 Analytical model.

Fig. 3 Temperature dependence ofYoung's modulus [11].

Fig. 4 Temperature dependence ofyield stress[11].

Fig. 5 Boundary conditions of heat conduction and hydrogen diffusion analysis.

― 4 ―SENAC Vol. 54, No. （2 2021. 4）

2.2 熱伝導解析
熱伝導解析には式(2)用いた。

  2T  2T 
T k  2  2 cthy  x(2).
式(2)を Crank-Nicolson 法で離散化した後に、SOR (Successive-Over-Relaxation)法で解いた。
初期温度は、WM を 1500oC とし、予熱無しの場合、BM および HAZ は室温(20oC)とした。予熱条件では BM および HAZ を 50、75、100、125、150oC とした。本解析では WM が 1500oC からの冷却のみを模擬した解析となっている。熱伝導解析に用いた物性値を Table2 に示す[13]。
Table 2 Material properties used for heat transfer analysis [13].
oka0.054
12.0･c

4.53×10-52.09×10H-3J/( C･mm･sec)
2mm /seco3J/( C･mm )o2J/( C･mm ･sec)

2.3 応力解析
応力解析には 2 次元弾塑性有限要素法解析プログラム EPIC-I に熱応力解析のルーチンを組み入れた[4,14]。応力ひずみ則は熱伝導に伴う温度変化に対応するようにし、降伏・塑性変形は温度に対応した応力ひずみ則に従って進行するようにした。弾性域ではフックの法則に従い、塑性域では式(3)に示す加工硬化則を用いた。本解析では m1=1 とする線形加工硬化則を用いている。

Hp m1 1d m1c1  *  pd p.

(3)

2.4 水素拡散解析
水素拡散の基礎式 [15]を式(4)に示す。

VC  DC  DC
 p 
t CRT.

(4)

第 1 項は濃度勾配による拡散項、第 2 項は応力勾配による拡散項をそれぞれ示す。ここで、右辺第 2 項の値は第 1 項よりも低い値となるため、濃度勾配が存在する場合は第 1 項の影響のみが発現し、第 2 項の影響が発現しないことが報告されている[2,3,9]。この場合、各項の影響が同程度となるように重み係数を導入することで、応力誘起の拡散･凝集効果が発現することが明らかとなっている。本手法はマルティプリケーション法として提案されている[2,3,9]。式(4)を展開し、各項に重み係数を導入すると式(5)のようになる。

DV
CDV 1 D 2C   2C p   3
C  2 p .
RT
t C
RT(5)

本解析では、各項の影響が同程度となるように1:2:3=1:1000:5 とした。本解析は、温度変化を考慮した解析であるため、拡散係数や熱応力も位置(x, y)の関数となる。これを考慮し、式(5)について、式(6a)-(6e)を用いて無次元化すると式(7)のようになる。これを Crank-Nicolson 法で離散化した後に、SUR (Successive-Under-Relaxation)法を用いて数値解析を行った[2,3]。
x xl0,(6a)数値解の安定条件を考慮した y 形溶接構造の冷却過程における水素凝集挙動解析および本解析プログラムによる数値解析高速化へ向けての考察― 5 ―y yl0 ,(6b)C CC0,(6c)D DD0,(6d)t DtC  C 2 0l0 .

(6e)   2 C   2 C    D  C  D  C   C   1  D  
22 
t Cy  y   y    x x  x  V   C   D   p D   p  D   C    p C    p    2 y y   T  x  x y y  R   T  x  x 22D  C   T   p T   p  D  C     p   p  3y  y   T 2  x  x T  x  2y  2  (7)係数の物理的意味については、濃度勾配および応力勾配を駆動力とする拡散過程において、その駆動ポテンシャルの相違に起因する両者の拡散係数のエントロピーの比を表すものとされている[16]。を考慮しない（=1 とする）ことは、異なる駆動ポテンシャルの下での拡散過程における原子配列の相違に起因するエントロピー変化を同じと仮定していることとなる。はこのようなエントロピー変化を考慮した係数である。
また初期条件として、溶金内の水素濃度を C/C0=3.0 とし、それ以外を C/C0=1.0 とした。本解析で用いた物性値を Table 3 に示す[17,18]。
Table 3 Material properties used for hydrogen diffusion analysis[17,18].

D0-6QR

5.54×10326.81×108.314m/sec
J/mol
J/K molV2.0×10-6m /mol

33. 解析結果
3.1 予熱を伴わない場合の冷却過程における水素凝集挙動
熱伝導-熱応力-水素拡散連成解析を用いた先行研究の結果を Fig. 6 に示す [4,5]。Fig. 6 に示すように本解析手法から得られる水素分布は、溶接底鈍角側に水素が凝集した。この結果は、Fig. 7 に示す実際の y 形溶接割れ試験で得られるき裂発生位置と良い一致を示す[19]。
次に、Fig. 8 に示す溶接底鈍角側を含む x=62、y=6 座標における静水圧応力分布の時系列変化
を Figs. 9-10 にそれぞれ示す。これらの結果から、静水圧応力は時間経過とともに増加し、溶接底鈍角側において極大値を示した。このことから、水素は静水圧応力の極大値領域に凝集することが示された。

― 6 ―

SENAC Vol. 54, No. （2 2021. 4）

Fig. 6 Hydrogen concentration distribution after coolingobtained by analysis[4,5].

Fig. 7 Result of y-grooved weld cracking
test [19].

Fig. 8 Collection position of hydrogen and hydrostatic stress .

Fig. 9 Time sequential change of hydrostatic stress
distributions at y=6 [4,5].

Fig. 10 Time sequential change of hydrostatic
stress distribution at x=62 [4,5].

3.2 冷却時の静水圧応力の時系列変化に及ぼす予熱の効果
2.2 節で示したように溶接部冷却過程において、水素は溶接底部の鈍角側に凝集した。そこで、それぞれの予熱条件での溶接底部鈍角側における静水圧応力の時系列変化を Fig. 11 に示す。Fig. 11 より、予熱温度の上昇に伴って静水圧応力の値は低下した。また、時間経過に伴う静水圧応力の低下率は予熱温度が上昇すると大きくなることがわかった。

数値解の安定条件を考慮した y 形溶接構造の冷却過程における水素凝集挙動解析および本解析プログラムによる数値解析高速化へ向けての考察

― 7 ―

3500

Hydrostatic Stress p, MPa

3000
2500
2000
1500
1000

NoPreheatPre 50oC
Pre 75oC

500
0400800Time, sec

Pre 100oC
Pre 125oC
Pre 150oC
12001600

Fig. 11 Time sequential change of hydrostatic stress at the bottom of weld metal

3.3 冷却時の水素濃度の時系列変化に及ぼす予熱の効果と数値解の安定性
次に、溶接底部鈍角側におけるそれぞれの予熱温度での水素濃度の時系列変化を Fig. 12 に示す。
Fig. 12 より、予熱温度の上昇に伴い、水素凝集濃度のピーク発生時間が早くなることが分かった。
冷却後の水素濃度は予熱温度 100oC において最も低くなり、100oC 以上では逆に増加することがわかった。また、予熱 100oC 以上では、水素濃度の時系列変化は若干不安定な挙動を示した。この原因について、横堀らは拡散方程式の数値解の安定性（振動解の制御）について次のように述べている [2,3, 20]。Clank-Nicolson 法は、位置と時間の刻み幅に関する制限はないとされているが、微分方程式の性質によっては収束するものの、Fig. 13 に示すように振動を示す場合がある。そこで、解が振動せず安定に収束する条件として、式(3) の条件を満たすことが必要であると指摘されている。ここで t+は無次元化時間、r+は無次元化距離である。

t r  2 c .

(3)

そこで、これまでの解析ではc=0.3~0.42,3)として計算していたところを、c=0.1 に変更して再解析を行った。c=0.1 として水素濃度分布の時系列変化を解析したものを Fig. 14 に示す。Fig. 14 より、予熱 100oC 以上の条件にて、t=400~800sec 付近で見られた不安定な挙動が改善されたことがわかる。このことから、式(3) に示すc の条件を適切に考慮することで振動解を抑制し、より安定な解析結果を得ることができた。

― 8 ―

SENAC Vol. 54, No. （2 2021. 4）

Hydrogen concentration C/Co43NoPreheat
Pre 50oC
Pre 75oC

Pre 100oC
Pre 125oC
Pre 150oC

210400800Time, sec

12001600

Fig. 12 Time sequential change of hydrogen concentration at the bottom of weld metal (c=0.3)

Fig. 13 Example of oscillatory solution for stress induced hydrogen diffusion analysis[20]数値解の安定条件を考慮した y 形溶接構造の冷却過程における水素凝集挙動解析および本解析プログラムによる数値解析高速化へ向けての考察― 9 ―Hydrogen concentration C/Co

43NoPreheat
Pre 50oC
Pre 75oC

Pre 100oC
Pre 125oC
Pre 150oC

210400800
Time, sec

12001600

Fig. 14 Time sequential change of hydrogen concentration at the bottom of weld metal (c=0.1)

4. 本プログラムによる数値解析高速化に向けて溶接部冷却時の水素凝集挙動に及ぼす予熱の効果は、実際には 24 時間以上の経過を見る必要があることから、長時間の計算に適応するために数値解析の高速化が望まれる。
本プログラムは、1) 熱伝導解析、2) 局所応力誘起水素拡散解析、3) 有限要素法弾塑性熱応力解析の連成解析を行っている[4~6,8,21]。1) および 2) は、差分法により離散化された多元連立方程式を SUR 法で解いており、3) は共役傾斜法（反復法）で多元連立方程式を解いている。また、
2) はC なる 1 次微分項が含まれており、振動解を生じる可能性があるため、3.3 節で述べたように、<c なる条件と SUR 法を用いている。本手法は差分法と有限要素法の特徴を生かして両者をハイブリッドさせたプログラムであるが、差分格子点と有限要素法の節点との間での温度および応力のデータ相互補間を行うプログラムなど複雑なアルゴリズムも内在している。また、連成解析であることから、1) と 2) および 3) ではそれぞれに同じ冷却時間での値が必要となることから、両者での時間補間も行っている。
そこで、本プログラムの各サブルーチンで消費される CPU 時間を調べた。その結果、本プログ
ラムでの各サブルーチンでの計算時間は、ほぼ非対角行列を共役傾斜法で解く 3) の部分で費やされていることが分かった（90%以上）。有限要素熱応力解析では、行列の非零要素のみを記憶させて解析するという効率化を図って、当初から反復法である共役傾斜法を用いているが、ベクトル化が困難である。そのため、スーパーコンピュータを使用して、計算の高速化をはかり、長時間の水素凝集挙動を再現するためには、この部分を計算機が提供するライブラリの利用も含めた、ベクトル化や並列化に適したソルバーにして 3) の解析プログラムの構築を図っていくことが必要と考える。
そのためには、連立一次方程式の部分を LAPACK/BLAS に置き換え、 メモリ割り当てを静的割当(スタック使用)から動的割当(ヒープ使用) に変更することにより大規模計算に対する体制を整え、高速化・並列化とより多くのメモリを利用できるように改良することが必要と考える。

― 10 ―

SENAC Vol. 54, No. （2 2021. 4）

5.

結言
本研究では、溶接部冷却過程における水素拡散解析を用いて、y 形溶接構造の水素凝集挙動を解析し、以下の結言を得た。
1) 熱伝導解析、熱応力解析および局所応力誘起水素拡散解析を連成させるプログラムを開発して、溶接部の冷却時における水素凝集挙動時系列特性に及ぼす予熱の効果を明らかにした。
2) 水素拡散方程式には、C なる 1 次微分項が含まれており、振動解を生じる可能性があるために、<c なる条件と SUR 法を用いている。本手法は差分法と有限要素法の特徴を生かして両者をハイブリッドさせたプログラムであるが、差分格子点と有限要素法の節点との間での温度および応力のデータ相互補間を行うプログラムなど複雑なアルゴリズムも内在している。計算時間の大半は非対角行列を共役傾斜法で解く有限要素熱応力解析の部分で費やされており、この部分の改善により計算の高速化が期待できる。
謝辞
本研究は、東北大学サイバーサイエンスセンターとの共同研究「マルチマテリアルにおける水素拡散行数挙動に及ぼすポテンシャル誘起駆動力特定解析プログラムの開発」の一環として実施されました。また、本研究は、総合科学技術・イノベーション会議の SIP（戦略的イノベーション創造プログラム）「革新的構造材料」（管理法人：JST）によって実施されました。
参考文献[1] N. Yurioka and T. Kasuya, Quarterly Journal of Japan Welding Society, Vol. 13, No. 3, pp.347-357 (1995).
[2] A. T. Yokobori, Jr., T. Nemoto, K. Satoh and T. Yamada, Engineering Fracture Mechanics, Vol. 55, No. 1,pp. 47-60 (1996).
[3] A. T. Yokobori, Jr., T. Ohmi, T. Murakawa, T. Nemoto, T. Uesugi and R. Sugiura, Strength, Fracture andComplexity, Vol. 7, No. 2, pp. 215-233 (2011).
[4] G. Ozeki, A. T. Yokobori, Jr., T. Ohmi, T. Kasuya, N. Ishikawa, S. Minamoto and M. Enoki, Proceedingsof the ASME 2018 Pressure Vessels and Piping Conference, PVP2018-84178.
[5] G. Ozeki, A. T. Yokobori, Jr., T. Ohmi, T. Kasuya, N. Ishikawa, S. Minamoto and M. Enoki, Advanced
Materials Letters, Vol. 9, No. 10, pp.677-683 (2018).
[6] A. T. Yokobori, Jr., G. Ozeki,. Ohmi, T. Kasuya, N. Ishikawa, S. Minamoto and M. Enoki, Materials
Transactions, 60,2(2019),pp.222-229.
[7] S.V. Patankar, Numerical heat transfer and fluid flow, Trans:Y. Mizutani and M. Katsuki, p. 70, Morikita
Pub. (1985), in Japanese.
[8] 尾関郷, 横堀壽光, 大見敏仁, 糟谷正, 石川信行, 源聡, 榎学, 日本材料強度学会学術講演会
講演論文集, pp.19-28 (2019).
[9] T. Ohmi, A. T. Yokobori, Jr. and K. Takei, Defect and Diffusion Forum, Vol. 326-328 (2012), pp.626-631.
[10] JIS Z 3158, Method of y-groove weld cracking test (2016).
[11] Y. Mikami, N. Kawabe, N. Ishikawa and M. Mochizuki, Quarterly Journal of the Japan Welding Society,Vol. 34, No. 2, pp.67-80 (2016), (in Japanese).
[12] K. Satoh, T. Terasaki and Y. Yamashita, J. of the Japan Welding Society, Vol. 48, No. 7(1979), pp.504509 (1979), (in Japanese).
[13] T. Kasuya, N. Yurioka, Welding Journal, Vol.72, No.2, p.107 (1993).
[14] Y. Yamada:, Sosei - Nendansei, pp. 180-219, Baifu-kan Pub. (1980), (in Japanese).
[15] Oriani, R. A., Stress Corrosion Cracking and Hydrogen Embrittlement of Iron Base Alloys, NACE-5,pp. 351-358 (1973).
[16] A. T. Yokobori Jr. and T. Ohmi, Strength, Fracture and Complexity, Vol. 8, No. 2, pp.117-124 (2014).
[17] T. Kasuya, Y. Hashiba, H. Inoue, S. Nakamura, T. Takai, Welding in the world, Vol.57, No.4, p.581(2013).
[18] Van Leewen, H. P., Corrosion, Vol. 31, No. 2, pp. 42-50 (1975).
[19] N. Ishikawa, SIP meeting at the University of Tokyo, Dec. 2017.
[20] 横堀壽光, 大見敏仁, 根本剛直, 上杉智治, SENAC 東北大学大型計算機センター広報, Vol.42,
No. 4, pp.27-35 (2009).
[21] A. T. Yokobori, Jr., G. Ozeki,. Ohmi, T. Kasuya, N. Ishikawa, S. Minamoto and M. Enoki, Proceedingof the ICF14, CDrom (2017)