B16 MLPG 法による熱弾塑性解析法の開発と溶接問題への応用
堀友(指導教員：正岡孝治､柴原正和)

則Development of Thermal Elastic Plastic Analysis using Meshless Local Petrov-Galerkin Method and the Application for Welding Problem
byTomonori HoriAbstract
A meshless local Petrov-Galerkin method(MLPG) for the analysis of the thermal elastic-plastic problem is presented in thispaper. The advantages of meshless method are that no element connectivity and no remeshing are required in pre and postprocessing. The MLPG method uses the moving least squares (MLS) to approximate the field functions and the Heaviside stepfunction as the test function of the weighted residual method in each sub-domain. The advantage of this method is to avoid the needfor a domain integral and only needs to integrate the local boundary, if body force is ignored. In the formulation of plasticity theory,flow theory is adopted. It is suitable to solve problems under complex loading conditions involving unloading. As a numericalexample, the present method is applied to simulate the bead on plate welding in order to consider the influence of specificparameters related to MLPG parameters. The results with those parameters show that the welding deformations and the residual
stress distributions along welding line have a good agreement with those of FEM. These observations show that the proposed
method is valid and useful for the solution of the welding problem.

1. 緒

言

現在，数値計算による構造解析手法として，最も発展し，用いられている解析手法に有限要素法(Finite Element
Method 以下 FEM)がある．FEM を用いて解析を行う場合，解析対象を要素に分割する必要があるが，解析対象の形状および力学現象等により，最適な要素分割は異なり，経験に依存するところが大きい．
そこで，本研究ではメッシュレス法のひとつであるMLPG(Meshless Local Petrove-Galerkin)法を用いた構造解析手法に注目した．MLPG 法では解析に要素分割を全く使用せず，要素・節点間におけるコネクティビィティ情報が不要となるため，プリプロセッシング時間の大幅な縮小が期待できる．さらに，再解析のために要素再分割を行う必要がないということも利点として挙げられる．
これまでに，MLPG 法を用いた，構造解析では，大変形解析 1)，熱弾性解析 2)，弾塑性解析 3)などが実施されている．本研究では，新たに MLPG 法を用いた溶接構造解析手法を確立することを目的とし，MLPG 法による熱弾塑性解析手法の開発を行い，本解析手法の基礎的特性を示すために溶接継ぎ手問題に対して本手法を適用・検討した．

ある領域 Ω 内における任意の評価点の位置ベクトルをx ，近似関数を u h (x ) とすると次式のように表される．
u h (x ) =m

∑ p (x)a (x) = p (x)a(x) x ∈ ΩJ(1)

T

JJ =1

ここで， p T (x ) = {p1 (x ), p2 (x ),…, pm (x )} は基底関数であり 2次元解析においては以下のようなものを用いる．

p T (x ) = {1 , x , y} 1 次基底：m=3
p T (x ) = {1 , x , y , x 2 , y 2 , xy} 2 次基底：m=6(2)(3)

a(x ) は未定係数ベクトルであり評価関数J (a(x )) =N∑
I =1[[w(x − x I ) p T (x I )a(x ) − uˆ Iuˆ T = uˆ 1 , uˆ 2 ,…, uˆ N2](4)](5)を最小化することによって決定する． û I は座標 x I における見かけの節点物理量， Ω SI は各節点がもつ影響領域，
ρ I は影響領域の半径， w x − x I は節点 I を含む影響領域 Ω SI 内における重み関数である．本解析において影響領域に Fig.1 のような円を用いている．N は w x − x I ≠ 0 となる影響領域内に含まれる節点の個数を表す．

()

()

ΩSI

ρI

2. MLPG 法による熱弾塑性解析理論
2. 1 内挿関数の作成
メッシュレス法では，点ごとに内挿関数を作成するので，注目する局所領域における点と点の間の物理量を補間する必要がある．このとき，良く用いられる手法として，移動最小二乗近似（Moving Least Square Approximation
以下 MLS）がある．以下にその概要について示す．

ΩSxrIxIΩx

Fig.1 Support domains which are related to the points ofinterest.

B16-1

大阪府立大学 工学部 海洋システム工学科 平成 20 年度 卒論・修論概要集(2009 年 2 月)

重み関数としては，4 次のスプライン関数を用いている．
式(4)の停留値から次の関係が得られる．
(7)
Α(x )a(x ) = B(x )uˆ

ここで A(x ) ， B(x ) は以下のように表される．
A(x ) =N∑wI p(x I )p T (x I )

(10)

I

∑

[]p J (x) A −1 (x )B(x )J =1LS−(11)ttt i + dt i = t i + d t itiiti

on Γu(13)on Γt(14)α∫ ( u + du − u − du )v dΓ = 0iit(20)iiitiiΓst∫ ( b + db )dΩ + α ∫ ( u + du )dΓ − ∫ σ n dΓ − α ∫ u dΓtitiiti

Γsutij j(22)iΓsuLSiΓsu∫t + dtt + dttdΓ

ΓstbdΩ + αΩs∫∫t +dtu dΓ − n σdΓ − αΓsut∫ udΓ(24)tΓsuLS

熱弾塑性解析理論は材料の温度依存性と非線形性を考慮するために，増分型の構成式を用いる．本研究においては，金属の溶接問題を取り扱うので，降伏条件に VonMises の降伏条件を用いた．また，塑性ひずみ増分は，降伏関数の応力成分の微係数に比例する流れ理論による．
(25)
dσ = Ddε − gdT
ただし， dε は全ひずみ増分， dT は温度増分を表し， Dと g は弾性域においてdE 1 e −1 ⎞⎛(26)
D = De , g = Dp ⎜ a −D σ⎟dT E⎠⎝
塑性域においては以下のものを用いる．
dE 1 e −1 ⎞⎛D = Dp , g = Dp ⎜a −D σ ⎟ − 2σ Y H T D e φ/S (27)dT E⎠⎝
変位増分ベクトルは MLS により得られる内挿関数を用いて以下のように表される．

)

tii

Γsu

−tij , j + dσ ij , j + b j + db j vi dΩ +Ωsi

t

Γsu

LS

式(12)に対して，重み付残差法を用いて弱形式化する．また，ペナルティー法により基本境界条件を付与すると評価点近傍の積分領域 Ω S において式(16)となる．
t)ij , j + dσ ij , j vi , j dΩ∫ ( b + db )v dΩ + α ∫ ( u + du − u − du )v dΓ = 0∫ ndσdΓ + α ∫ dudΓ = − ∫ここで， σ ij は応力テンソル， b j は体積力， Γu は基本境界条件を与える境界面であり， ui は節点 i の変位， Γt は自然境界条件を与える境界面である．また， t i は以下で示すような境界面の法線方向の応力である．
(15)
t i = t σ ij n j , dt i = dσ ij n j∫(σti

ここで，式(22)を行列・ベクトル形式に書き直す．

2. 2 MLPG 法熱弾塑性解析理論の定式化
本節においては，本研究において開発した，熱弾塑性
MLPG 法理論の概要について示す．
時刻 t における解析対象の全体領域 Ω ，境界を Γ としたとき応力の釣り合い方程式および境界条件はt(12)σ ij , j + dσ ij , j + t b j + db j = 0 in Ωu i + du i = t u i + du iiΩsΩsJIti

Γstij j以上の式から得られた，以下の式が本 MLPG 法で用いる内挿関数となる．
mti∫ dσ n dΓ + α ∫ du dΓ = − ∫ ( t + dt )dΓ

I=1φ I ( x) =i

本研究において test function にヘビサイド関数を用いる．
⎧1 at x ∈ Ω S(21)v( x) = ⎨⎩0 at x ∉ Ω SΩ S 内において，ヘビサイド関数は一定なので，vi , j = 0となる．これを式(20)に適用し，式(15)を用いると(9)れる．
ILSi

Ωs

式(9)を a(x ) について解くと，以下で表す近似関数が得ら
Nt+I=1∑φ (x )û∫ ( t + dt )v dΓ + ∫ ( t + dt )v dΓ − ∫ ( σ(8)B (x ) = [w1p (x1 ), w2 p (x 2 ),… , wN p (x N )]u h (x ) =式(17)と自然境界条件を式(16)に適用するとJJ⎧⎪du h ⎫⎪⎡φ J0 ⎤ ⎧ duˆ1 ⎫JJ(28)du (x ) = ⎨ 1h ⎬ = ⎢⎬ = Φ duˆJ ⎥⎨⎪⎩du 2 ⎪⎭⎣ 0 φ ⎦ ⎩duˆ 2 ⎭式(24)に式(25)を代入する．また，式(28)の関係を用いて，(16)

hJi iΓsuvi は重み関数(test function)， α はペナルティー係数を表し α >> 1 である．次に式(16)に対して，次の関係式を用いる．
(17)
σ ij , j vi = (σ ij , j vi ) , j − σ ij vi , j∫ n DB duˆ dΓ + α ∫ SΦ duˆ dΓ = − ∫τJJJΓsuLS−J∫t + dtbdΩ + αΩs∫t + dtt + dtt dΓ(29)Γst∫u d Γ − n t σdΓ − αΓsu

∫

tΓsu

LS

∫ud Γ + nτ gdTdΓ
LS

τ は収束計算により，残差力が十分に小さいときの値である．
N∑ τ K duˆ = df + df + df
IJ0IJ1I(30)2IJ =1

ただし，
τ∫∫

K IJ = nτ DB J dΓ + α SΦ J dΓΓsuLS∫df I0 = −t + dtt dΓ −Γst∫∫df I2 = nτ gdTdΓ

Fig.2 Schematic illustration of MLPG method

LSB16-2t + dt

bdΩ + αΩsdf I1 = − n t σdΓ − αLS∫Γsu∫ ud ΓtΓsu∫

t + dtu dΓ平成 20 年度 卒論・修論概要集(2009 年 2 月)
3.5

Material Properties

B=200[mm]σ y (× 10 2 MPa) ν=0.3=const

3.0Ayα =1.2=const(× 10 -5 ℃ -1 )2.5

E(× 10 2 GPa)

2.0

TS

1.5

SS

1.0d0.5x0.0

L=200[mm]Fig.3 Bead on plate model

0500
1000
Tem parature [℃ ]1500

Fig.4 Temperature dependent material property

Fig.5 Geometrical Parameters of MLPG

3.MLPG 法による溶接変形・残留応力解析
3.1 解析モデル
Fig.3 に本章に用いた解析モデルを示す．板厚方向に平面応力状態を仮定し，境界条件については図に示した通りで，対称線上の全節点の y 方向変位を拘束および A 点に x 方向変位を一点拘束する．今回用いたモデルに対する節点分布として，節点数は 861 で，板長方向に 41 点，板幅方向に 21 点で均一な節点分布となっている．
溶接条件は，入熱量 105[J/mm]，溶接線速度 v=5[mm/s]，熱効率 0.7 とし，平板表面における熱伝達の影響を考慮する．また，Fig.4 に材料定数の温度依存性を示す．

Fig.6 Two types of in-plate deformation mode
横収縮については，SS の大きさにかかわらず，FEM による解析結果と良好に一致することが確認できた．一方で，縦収縮は SS の大きさが小さい範囲においては，FEM による解析結果と比べてずれがある．しかし，SS が 2.1d 以上の値を用いるとき，FEM 解析結果とほぼ完全に一致することが確認できる．溶接線近傍における溶接線方向の残留応力 σ x は，MLPG，FEM ともに降伏応力に近い値となっており，SS の大きさにかかわらず，良好に一致していることが確認できる． σ y についても FEM 解析結果と非常に合っていることがわかる．以上結果から，本提案手法を用いて，解析を行うとき，溶接変形，残留応力が共に，FEM 解析結果と良好に一致するのは，SS が 2.1d以上を用いた場合である．

3.2MLS 影響領域の最適値検証
Fig.5 に示す，MLPG 特有のパラメータである影響領域半径(Support size 以下 SS)と積分領域半径(Test functionsize 以下 TS)がある．ｄは節点間の距離を表す．SS の大きさが異なると，各積分領域において得られる内挿関数の性質が異なるため，積分精度に大きく影響を与える．
TS は，その大きさが小さすぎると，Collocation methodのような強形式で解析することに近い状態になることが知られている 4）．これは領域というより点で解析する状態に近くなるためである．本解析においては，節点間の距離 d は均一な節点分布のため一定で TS は 0.6d 付近の値を用いて解析を行っている．以上の条件のもとで，内挿関数作成に用いる基底ベクトルの 1 次基底(MONO)および 2 次基底(POLY)それぞれについて，SS と Fig.6 に示すような，横収縮，平均縦収縮について，FEM による解析結果と比較・検討を行った結果を Fig.7 に，残留応力分布について，その結果を Fig.8 に示す．

0-0.2-0.4
-0.6
1.5

1.7

1.9
2.1
Support size SS/d2.3

2.5-0.02

Proposed Method (MONO)
Proposed Method (POLY)
FEM-0.03
-0.04
-0.05
-0.06
-0.07
1.51.71.9
2.1
Support size SS/d

2.32.5250
200FEMMONO SS=1.5d
MONO SS=2.0d
MONO SS=2.1d
MONO SS=2.3d
POLY SS=2.1d
POLY SS=2.2d
POLY SS=2.3d150
100500-50050100150
X coordinate along welding line [mm](a)Stress distribution σ x200100T300

0FEM
MONO SS=1.5d
MONO SS=2.0d
MONO SS=2.1d
MONO SS=2.3d
POLY SS=2.1d
POLY SS=2.2d
POLY SS=2.3d-100-200-300

200 -400

Transverse shrinkage (S /2) [mm]350Residual stress distributions [MPa]Residual stress distributions [MPa](a) Transverse shrinkage(b) Longitudinal shrinkage
Fig.7 The Influence of Support size (SS)050100150
X coordinate along welding line [mm]2000
Proposed Method (MONO)-0.1

Propsed Method (POLY)
FEM

-0.2
-0.3
-0.4
-0.5
-0.6
-0.7

(b) Stress distribution σ y0

Proposed Method (MONO)
Proposed Method (POLY)
FEM-0.01

L0.2

0-0.01Longitudinal shrinkage (S /2) [mm]Longitudinal shrinkage SL/2 [mm]

Propsed Method (MONO)
Propsed Method (POLY)
FEM

0.4

T

Transverse shrinkage S /2 [mm]0.6

3.3 入熱量の影響
本節では 3.2 節において，FEM 解析と良好に一致した
MLPG パラメータ SS=2.1～2.3d を用いて，各種溶接条件における解析結果の検証を行う．
溶接線速度 3[mm/s]で一定とし，単位長さの溶接線に与え ら れ る 入 熱 量 98[J/mm] ， 140[J/mm] ， 168[J/mm] ，210[J/mm]，280[J/mm]のそれぞれについて解析を行い，横収縮，平均縦収縮を求め，FEM による解析結果と精度比較を行った．また，溶接線方向残留応力について，入熱量 280[J/mm]のときの分布を示している．
Fig.9，Fig.10 に，それぞれの結果を示す．
入熱量が大きくなるにつれて，溶接変形が大きくなるので，一般に解析は難しくなるが，Fig.9 からわかるように，2 次基底を用いた解析においては，FEM による解析結果と良好に一致していることがわかる．一方，1 次基底関数を用いた解析結果に関しては，横収縮量については定性的にも定量的にも良好に一致している．しかし，縦050100150200
Heat input [J/mm]

250300-0.02-0.03-0.04-0.05050

100150200
Heat input [J/mm]

Fig.9 Influence of heat inputs on welding deformation

Fig.8 Residual stress distributions on welding line
B16-3250300

4000300

Present Method (MONO)

Residual stress distributions [MPa]大阪府立大学 工学部 海洋システム工学科 平成 20 年度 卒論・修論概要集(2009 年 2 月)

Q =280[J/mm]

200
100
0

Proposed Method (MONO) σ x
Proposed Method (MONO) σ yProposed Method (MONO) σ xy
Proposed Method (POLY) σ x
Proposed Method (POLY) σ y
Proposed Method (POLY) σ xy
FEM σ x
FEM σ yFEM σ xy

-100-200-300-400050100
150
X coordinate along w elding line [mm ]

-0.2
-0.3
-0.4
-0.5
-0.6200

Proposed M ethod (M ONO)
Proposed M ethod (POLY)
FEM-0.102468v[mm/s]

101214Fig.12 Transverse shrinkage in case of influence of welding speeds

Fig.10 Residual stress distributions on welding line[×10MPa]

Present Method (MONO) σx400(b) FEM(a ) Propsed method (MLPG)
Fig.11 The accuracy of stress distribution σ x300v=13[m m /s]

2001000

Proposed Method (MONO) σ x
Proposed Method (MONO) σ y
Proposed Method (MONO) σ xy
Proposed Method (POLY) σ x
Proposed Method (POLY) σ y
Proposed Method (POLY) σ xy
FEM σ x
FEM σ yFEM σ xy

-100-200-300-400050100150X coordinate along w elding line [mm ]200

Fig.13 Residual stress distributions in case of influence of welding
speed on welding line[mm](b) FEM(a) Proposed method (MLPG)
Fig.12 The accuracy of deformation distribution in y directionがわかる．一方で，1 次基底関数に関しては，v=13[mm/s]のときにおける横収縮の解析結果が他の解析結果よりずれが大きい，しかし，他の溶接速度においては，定量的に，ほぼ一致している．残留応力分布については，1 次基底，2 次基底共に，溶接始端部，および終端部においてもFEM 解析結果と良好に一致していることが確認できる．

4.結言
収縮に関して，入熱量の増大と共に増加する傾向は両者で一致するがやや大きめな値となっている．残留応力分布については，Fig.10 から，2 次基底を用いた解析においては，FEM 解析結果と良好に一致している．一方で，1次基底による解析で x=125[mm]付近の y 方向応力の分布がやや小さく表れているが，その分布傾向は一致している．また，x 方向応力，せん断応力については FEM 解析結果と良好に一致していることがわかる．
次に，MLPG 法による解析が解析領域全体においても十分な解析が行えていることを示すために，解析領域全Fig.12体のコンター図を示す．Fig.11 は x 方向応力分布図，は，y 方向変位を示している．FEM 解析と MLPG 法による解析結果を比較すると，変位，応力共に，両者の分布形状および，定量的な分布についても，ほぼ一致していることが確認できる．

本研究において，新たに，MLPG 法による熱弾塑性解析法の開発を行い，基本的な溶接継手問題に適用し，FEM解析結果と比較して以下の知見を得た．
1) 本提案手法を用いることにより，溶接変形，残留応力の解析を実用上十分な精度で解析することが可能であることを確認した．
2) MLPG 法による解析において，重要なパラメータの一つである，影響領域 SS について最適値を求めた結果，本解析条件においては，2.1d 以上を用いるときFEM 解析と良好に一致することを確認した．
3) 入熱量や溶接速度といった溶接条件の違いにかかわらず，最適な影響領域 SS，および，MLS 基定関数を用いることにより，FEM 解析とほぼ同等の結果を得ることができることを確認した．

3.4 溶接速度の影響
本節においても 3.3 節と同様に，SS=2.1～2.3d を用いて，次の溶接条件において解析を行う．入熱量 Q=126[J/mm]で一定のもとで，溶接速度 v=3,7,10,13[mm/s]のときの溶接変形・残留応力について評価を行う．溶接速度が速くなるほど，トーチの通過時間が短くなるので，熱伝導しにくくなり，温度増分が大きく解析が難しくなる．Fig.12には，横収縮および，溶接速度 13[mm/s]のときにおける残留応力分布を示している．
横収縮に関して，Fig.12 からわかるように，2 次基定関数を用いた解析において，横収縮は，溶接速度の速さにかかわらず，FEM 解析の結果と良好に一致していること

1)

B16-4

参 考 文 献2)3)

4)Ching H.K, Chen J.K. ： Thermomechanical analysis offunctionally graded composites under laser heating by theMLPG method, CMES, vol13, 2006, pp.199-217
Han Z.D, Rajendran A.M, Atluri S.N.：Meshless Local
Petrov-Galerkin(MLPG) approaches for solving nonlinear
problems with large deformations and rotations, CMES,
vol10, 2005, pp.1-12
Long S.Y, Liu K.Y, Li G.Y. ： An Analysis for theElasto-Plastic Fracture Problem by the Meshless Local
Petrov-Galerkin Method, CMES, vol28, 2008,
pp.203-216
Liu G.R, Gu Y.T.：An Introduction to Meshfree Methods
and Their Programming, 2005, Springer