M01 GPU 並列化を用いた理想化陽解法 FEM による大規模溶接力学解析システムの構築
生島一樹(指導教員 柴原・深沢・池田・馬場)

Development of the large scale analysis system for welding mechanics using GPU parallelized idealized explicit FEM
byKazuki Ikushima

Abstract
In ship building, joining process is necessary and welding technique is widely used because of its ease in utilizing. However,residual stress and deformation occurs with the practice of welding. And these residual stress and deformation raise the problem not
only in construction, but also under the operation. In order to solve these problems, quantitative prediction method is necessary andFE analysis is the most powerful tool in present. But even using FE analysis, the prediction is very difficult due to the large
computing time and memory consumption. This is because of the strong non-linearity of the welding phenomenon. In this research,
the authors proposed a new numerical method for welding mechanics based on the dynamic explicit FEM, named Idealized Explicit
FEM. In the proposed method, the displacements are computed using dynamic explicit FEM until the whole system reaches thestatic equilibrium state. And, to achieve the static equilibrium state faster, modified mass and damping matrix is introduced. The
proposed method and static implicit FEM are compared at the final path of multilayer welding of thick bead-on-plat. And it isverified that the transient and residual deformation and stress distribution of the proposed method show good agreement with thoseof static implicit FEM. In addition, the computing time and memory consumption of the proposed method are 1/12 and 1/40 times
shorter than those of static implicit FEM, respectively, in 243,243 degree of freedom model. In addition, parallelization using
Graphic Processing Unit (GPU) are applied to idealized explicit FEM to achieve more fast computation. The developed method isapplied to large scale welding problem that has more than 1 million degree of freedom. And as a result, it is found that thedeveloped method is very useful to analyze the large scale welding mechanics problem that has almost the same simulation scale as
the practical structures.

1. 緒

言

時間ステップに分割するので，溶接問題の解析においては，時間ステップが増大し，解析が困難になる．そのため，動的陽解法 FEM をそのまま溶接力学解析に適用することは困難である．
そこで，著者らは，動的陽解法 FEM を基に，溶接過渡応力変形解析手法である理想化陽解法 FEM の開発を進めてきた．その結果，理想化陽解法 FEM は良好な解析精度を有し，高速かつ省メモリに解析できることが確認した．
加えて，本研究では，近年普及しつつある計算手法である Graphic Processing Unit (GPU) による並列化を用いたシステムを構築することで，計算の一層の高速化を図り，実構造物クラスの解析規模を有する問題の解析を行い，開発した手法の有効性を確認した．

船舶に代表される大型構造物の製造において，部材を組上げる際には，部材同士を接合する工程が必ず存在するため，接合は必要不可欠な技術である．接合技術の中で，簡便さなどから，広く用いられているのが溶接である．
しかし，溶接の施工に伴い，残留応力や収縮・曲げなどの溶接変形が必ず生じ，構造物の製作のみならず，その使用中においても様々な問題を引き起こすことがある．
溶接変形や残留応力の問題を克服するためには，溶接変形・残留応力を定量的に予測することのできる手法が必要である．溶接変形・残留応力を予測する最も有力な手法は，計算機によるシミュレーションである．しかし，溶接現象は，材料の局所的な溶融を伴う非常に複雑な現象であるため，計算機を用いた解析では解析時間・メモリ消費量の両面から解析が困難である．このため，現在の解析技術では，適用範囲が溶接継手レベルに限られ，大規模問題への適用例が少ないのが現状である 1)．
一方，現在の解析技術において，大規模な問題を非常に短時間かつ，少ないメモリ消費で解析が可能な手法として，動的陽解法 FEM が挙げられる．動的陽解法 FEM は連立方程式の求解過程を含まない．そのため，非常に高速かつ省メモリな手法である．しかしながら，動的陽解法 FEM は，解析を行う際に，現象の時間を非常に小さい

2. 理想化陽解法 FEM の開発
2. 1 動的陽解法 FEM
理想化陽解法 FEM は動的陽解法 FEM を基にした手法であるが，動的陽解法 FEM は以下のような手順で定式化される．まず，節点変位ベクトル uを用いて，中心差分により節点速度ベクトル u および節点加速度ベクトルuを離散化する．
U t  U t  t  2U t  U t  t  t 2(1)U  U   U   2 t (2)tM01-1t  tt  t

平成 22 年度 修論・卒論概要集(2011 年 2 月)

Load step 3

Displacement

Load step 2

Static equilibrium
state

Load step
Load step 1tcr

Tstatic

Time step progression

Time step

Fig. 1 Concept of the proposed method.
ここで， t は時刻， t は時間増分を示し，ut  t ，ut ，ut  t はそれぞれ時刻 t  t ，t ，t  t における変位を示す．
M  U  C  U  BT  dV  F (3) tt1
 1
C U t t  2 M  2t tTF t  B   dV  2 2 M U tt11C U t t  2 M  2t
 t


t

ことができ，静的陰解法 FEM と同等の解を得ることが可能になると考えられる．本研究においては，この収束条件に関して，慣性力ベクトル，減衰力ベクトルのノルムと拘束節点反力ベクトルのノルムの比により動的効果の評価を行い，残差力ベクトルのノルムと拘束節点反力のノルムの比により静的平衡状態の評価を行うものとした．

K U   F   G  F  f nn(5)以上の手法をとることで，溶接問題，特に数万以上の自由度を有する大規模溶接構造物の過渡変形・応力問題に対して，従来の静的陰解法 FEM と同等の解析精度を実現し，計算時間の短縮，並びに省メモリ化を実現できるものと考えられる．
2.2.2 質量・減衰マトリックスの導出
理想化陽解法 FEM では，必要な計算ステップ数を少なくするために，以下のような質量マトリックス M  および減衰マトリックス C  の導出方法を提案する．
提案手法では，密度を要素の各自由度方向の最小サイズ
lmin_ i i : x, y, z  とヤング率 E を用いて以下のように表す．
2 i  a 2 t cr2 E l min_ i(6)
ここで， a は 1 以上の値で，今回は，要素内の応力波の伝播を十分に評価するため，10 とした．
式(6)により求めた要素の各自由度毎の密度を用い，次式(7)に示す積分を実行することにより質量マトリックスを求める．

(4)

M     dVVe

ここで，質量マトリックス M  と減衰マトリックス C が節点集中型，すなわち，対角成分のみ非零であるとす
ると，式(4)は，ut および ut  t を代入することで，ut  tを求めることができる．このことにより，計算規模が大きくなればなるほど陰解法と比べて計算時間が小さくなり，また，使用メモリも大幅に小さくなる．
しかしながら，クーラン条件と呼ばれる解析における時間増分値の制限により，計算時間が非常に長くなるという問題が生じる 7)．特に，溶接問題においては，加熱開始時から完全冷却時までを解く場合，約 1.9  10 10 もの膨大な時間ステップが必要となる．
2. 2 理想化陽解法 FEM の開発
2.2.1 解析手法の提案
本研究で提案する理想化陽解法 FEM では，通常の静的解析と同様に，熱伝導解析により得られた温度増分により発生する荷重を外力として用いる．そして，Fig. 1 に示すように，この温度増分による荷重を式(4)に与え，次節で示す収束性を高めるために調整を行った質量マトリッ
クス M  ，減衰マトリックス C  を用いて同式を解くことで各時間ステップにおける変位を求め，各荷重ステップにおいて変位が収束するまで式(4)の計算を繰り返す．このとき，各時間ステップにおいては，静的陰解法 FEM における反復計算と同様に，応力，塑性ひずみなどの物理量の更新は行わず，慣性項と減衰項の影響が無視できる程度に小さくなり，かつ，解析対象が式(5)で示される静的平衡状態を満たす場合にのみ解が収束したものとみなし，物理量の更新を行う．このような手法をとることで，慣性項と減衰項を付加したことによる動的効果を抑えるi(7)ただし，式(7)の積分は，行成分を対角項に集中させ，節点集中型のマトリックスとする．
減衰マトリックス C  は，静的解への収束の速さを重視し，1 次元振動理論を基に決定する．すなわち，1 次元振動理論において，質量を m ，減衰係数を c ，バネ定数を kとするとき， c  2 mk において臨界減衰となり，最も速く安定的に変位が静的解に収束することから，これを基に，減衰マトリックスの対角成分を次式のように決定する．

cii  2 mii kii(8)ここで， cii ， m ii ， k ii はそれぞれ減衰マトリックス，質量マトリックス，要素剛性マトリックスの対角成分であり，減衰マトリックス C  も，質量マトリックスと同様に対角マトリックスであると仮定する．
以上の方法により，クーラン条件を基に，要素サイズや材料定数の温度依存性をも考慮に入れた質量マトリック
ス，減衰マトリックスを導出することで，収束に要する時間ステップを削減することが可能になると考えられる．
2.2.3 理想化陽解法 FEM の性能開発した理想化陽解法 FEM を基礎的な溶接問題に適用し，手法の性能を検証する．解析に使用したモデルは Fig.
2 に示すビードオンプレート溶接問題である．また，解析に使用した材料の温度依存性データは，Fig. 3 に示す一般的な軟鋼材のものとした．
Fig. 4 に，理想化陽解法 FEM と従来手法である静的陰解法 FEM の y 方向残留応力と z 方向の残留変位の分布を示す．Fig. 4 より，理想化陽解法 FEM は静的陰解法 FEMの解析結果が良好に一致していることが分かる．

M01-2

平成 22 年度 修論・卒論概要集(2011 年 2 月)
Proposed method

(MPa)
250.0

torch
120.0 mm

Proposed method

(mm)
0.120

140.00.082

30.00.044-80.00.006-190.0-0.032・・・

40 divisionz60.0 mmy30 division
・・・

x20.0 mm-300.010 division-0.07
Static Implicit FEM

1000000

Physical constantsY (×102 MPa)

6.0c (×10 J/g/℃)
(×10-5 J/mm2/sec/℃)

(×10-2 J/mm/sec/℃)4.02.0(×10-1)
(×10-5 oC-1)

Computing time (s) (×10-3 g/mm3)
-1

18000

Proposed method
Static Implicit FEM

10.0

800000600000400000112200000

E (×105 MPa)

0.0
0.0

500.0

1000.0Temperature (℃)

Fig. 3 Temperature dependent
physical constants.

1500.000 20000

40000

60000

Number of Nodes

80000

100000

Memory consumption (MB)Fig. 2 Computational model forverification.

8.0

Static Implicit FEM(a) Stress distribution  y .
(b) Displacement in z direction.
Fig. 4 Comparison of computing results between proposed method and staticimplicit FEM.
Proposed method
Static Implicit FEM

15000
12000
9000

1406000
3000
00 20000

40000

60000

Number of Nodes

80000

100000(a) Stress distribution .
(b) Displacement in z direction.
Fig. 5 Comparison of computing time and memory consumption between proposed
method and static implicit FEM.

また，モデルの要素分割数を変更し，計算時間，及び，メモリ消費量について両手法を比較したものをそれぞれFig. 5 に示す．Fig. 5 より，理想化陽解法 FEM は静的陰解法 FEM に対して，高速かつ，省メモリに解析可能であり，8 万節点のモデルにおいては，計算時間，メモリ消費量をそれぞれ 1/12，1/40 まで削減できている．
このように，理想化陽解法 FEM は，従来手法とほぼ同等の解析精度を高速かつ省メモリに実現することが可能である．

3. GPU による並列化を用いた高速化
3. 1 理想化陽解法 FEM への GPU 並列化の導入
理想化陽解法 FEM は，動的陽解法 FEM に基づき変位を求め，得られた変位を用いて応力を求めるため，連立方程式を解く過程が含まれず，要素の応力積分を行う箇所の計算負荷が最も高くなる．
そこで，本研究では，GPU 上で要素の応力積分を行う手法を採った．GPU は，コンピュータグラフィックス(CG)処理に用いられてきたプロセッサであるが，近年の CG技術の発展を受け，非常に高速かつ汎用化されており，数値計算に応用可能であることが示されている 2)．
本研究で採用した手法では，GPU 上の 1 つの計算ユニ
ットで 1 つの要素の応力を算出し，その応力を積分することで，要素等価節点力を算出する．算出された要素等価節点力は，CPU 上ですべての要素等価節点力を足し合わせることで，解析領域全体の等価節点力ベクトルを作成する．また，変位の算出に関しても，CPU 上で計算した結果を GPU に転送する際の転送時間を考慮して，GPUの 1 プロセッサ上で 1 自由度の処理を行うような手法で計算するものとした．以上の解析手順をフロー形式で表M01-3

したものを Fig. 6 に示す．
このような手法を採ることで，理想化陽解法 FEM において，計算コストがもっとも高い，要素毎の計算，並びに自由度毎の計算を並列化することが可能となり，高速化が達成されるものと考えられる．
3. 2 GPU 並列化の効果
先ほど用いたモデルを用いて，計算時間の比較を行った結果を Fig. 7 に示す．計算に用いた GPU は，NVIDIA
GeForce 285GTX で，計算ユニット数は 240 である．計算時間の定義は，解析開始から 100 温度ステップまでの解析時間とした．
Fig. 7 より，静的陰解法 FEM は，計算時間が解析自由度の 2 乗に比例しているのに対して，GPU による並列化を適用した理想化陽解法 FEM と通常の理想化陽解法
FEM の計算時間は，解析自由度に対してほぼ線形となっていることが分かる．さらに，GPU を用いた並列化理想化陽解法 FEM は通常の理想化陽解法 FEM に対して，15倍以上高速に計算できている．また，GPU 並列化を適用した理想化陽解法 FEM は，静的陰解法 FEM に対しては，24 万自由度程度の問題において，180 倍以上高速に計算できていることを確認した．
3. 3 多層溶接問題への適用
以上のように，GPU による並列化を適用した理想化陽解法 FEM は，大規模な問題をさらに高速に解析できることが分かった．本節では，Fig. 8 に示す多層溶接問題をGPU 並列化を適用した理想化陽解法 FEM にて解析し，手法の有用性を示す．
本解析モデルの節点数は 427735，要素数は 412400 で，自由度数は 120 万以上である．溶接条件は溶接電圧を120V，溶接電流を 17A，溶接速度を 6.66mm/s とし，熱効平成 22 年度 修論・卒論概要集(2011 年 2 月)

Current : 120A
Voltage : 17V
Welding speed : 6.66mm/s
Heat efficiency : 0.5

Construct nodes, elements

σytorch(MPa)450.090mm

Renew temperature field
Compute idealized [M], [C]Converging loop

zGPU Computationx

Compute displacementby explicit FEM

29NoAll temperature
steps computed?

No313032-90.0(b) 20th pass
-270.0σy765-450.0

43Yes

σy90.033272823241920
1715141211926252122181613108Yes(a) 1st passyFig. 8 Mesh division and welding
condition of multi-pass model.

Compute nodal internal force

Static equilibrium state?

270.02z

1endyFig. 6 Computing flow of parallelized
idealized explicit FEM.

Fig. 9 Welding sequence of multi- passwelding(c) final pass (33rd)

Fig. 10 History of stress  y distribution.

4. 結
40000
Parallelized Idealized
Explicit FEM (GPU)Computing Time (s)30000

Idealized Explicit
FEM
Static Implicit
FEM

20000

10000

0Fig.

060000120000 180000
Degree of Freedom

240000

7 Comparison of computing time amongparallelized idealized explicit FEM, idealizedexplicit FEM and static implicit FEM.

率を 0.5，パス間温度を 250℃とした．パス数は 33 パスで，総荷重ステップ数は 56320 となった．また，本解析では，溶接の施工順序を Fig. 9 に示すようなものとした．
Fig. 10 に，溶接線垂直方向応力  y のパスごとの溶接線中央横断面における履歴を示す．Fig. 10 より，パスが進むに従い，肉盛が行われ，応力分布が詳細に解析されている様子がわかり，y 方向応力に関しては，表面で引張りの残留応力が生じ，内部では圧縮の残留応力となっている．これは，今回解析したような厚板の多層溶接においては，すでに溶接されて母材と一体となった下層の溶接部と母材が，各層の溶接金属が横収縮に母材を引寄せようとするのに対して抵抗するためである．以上のように，本事例では，多層溶接において特徴的な現象を十分に解析できている．
このように，本手法は従来手法では，詳細な解析が困難であった多層溶接問題に対して，詳細な検討を行うことができ，非常に有効な手法であるといえる．

M01-4

言

本研究では，溶接力学問題の大規模解析を目的とした理想化陽解法 FEM を開発し，基礎的な溶接問題に適用し，手法の検証を行った．また，理想化陽解法 FEM に対して，GPU による並列化手法を導入し，高速化を図り，多層溶接移動熱源問題に適用した結果，以下の知見が得られた．
1) ビードオンプレート溶接問題に関して，理想化陽解法FEM を適用した結果，残留変形・残留応力において，静的陰解法 FEM と同等の精度での解析可能であることを示した．
2) 理想化陽解法 FEM 及び静的陰解法 FEM における計算時間に関して比較を行った結果，理想化陽解法 FEM
は静的陰解法 FEM に対して，8 万節点のモデルにおいて，12 倍以上高速に解析可能であることを確認した．
3) 理想化陽解法 FEM は静的陰解法 FEM と比べて，非常に少ないメモリ消費量で解析を行えることを示した．
4) 理想化陽解法 FEM に対して，GPU 並列化手法を導入し，計算時間に関して，並列化前の手法と比較した結果，並列化によって 15 倍以上高速化したことを確認した．
5) GPU 並列化を用いた理想化陽解法 FEM を多層溶接移
動熱源問題に対して適用した結果，従来手法では困難であった，溶接順序を考慮した残留応力の解析が可能であることを示した．

参 考 文 献1)

2)

H. Nishikawa, I. Oda, M. Shibahara, H. Serizawa, H.
Murakawa: Three-dimensional Thermal-elastic-plastic
FEM Analysis for Predicting Residual Stress andDeformation Under Multi-pass Welding, Proceedings ofISOPE ‘04, Toulon, France, 2004, CD-ROM.
藤澤 智光，奥田 洋司:計算工学での GPU コンピュー
ティング, 計算工学, Vol. 15, No. 4, 2010, pp.5-6.