メッシュレス有限要素法による弾性体の解析
日大生産工（院）○大川 功次郎 日大生産工 登坂 宣好
１ まえがき
現在、我々が生活を営む上で逆問題は意識す るしないに関わらず、身近なところに存在して いる。その典型的な例として、原子力発電所の 配管の非破壊検査・地下資源探査・超音波診断 などが有る。
逆問題を解く上で大切なのが、逆問題の対の 考えになる順問題をいかに精度良く解くか、と いうことである。例えば、物性値同定逆問題で は求めたい物性値をまず仮定しFEMなどの数値 計算手法を用い変位を求める。その変位を測定 データとし、フィルタ理論やベイズ推定などの 逆解析手法を用いて、物性値を正解値に収束さ せてゆく。変位を求める部分が順問題であり、
物性値を収束させる部分が逆問題となる。
本論では、物性値同定逆問題の順問題にあた る2次元弾性体の解析を、要素分割を必要とし ないメッシュレス有限要素を適用し、その有効 性を検討する。
２ メッシュレス離散化
2.1 支配方程式
本手法では、2次元弾性体において次式のよ うに表された弱形式（仮想仕事の原理）を支配 方程式として用いる。
{ } { } { } { } { } { }
∫∫
∫∫ ∫
=
−
−
V
T T
V S
T
dV b u δ
dS t δu dV
σ δε
0
1 (１)
2 i
i u on S
u = (２)
ここに{ } { } { }σ 、ε 、u は、それぞれ応力、ひずみ、
変位を表し、 は表面力、体積力を表す。
また、Vは解析領域内部、 は自然境界条件を
与える領域境界、 は基本境界条件を与える 領域境界を表す。
} { } {t 、b
S1
S2
2.2 内挿関数の作成
本手法では、被積分関数を評価するために FEMで用いられているような要素単位での内挿 ではなく移動最小２乗法 （The Moving Least Square Method：以下MLSMと記す）による内挿 方法が用いられる。
) 1 (
MLSMを用いる場合、各節点における近似関数 を次式のように設定する。
φh
∑
= m
j j j
h(x,y) p (x,y)a (x,y)
φ
（３）
)}
, ( { )}
, (
{P x y T a x y
≡
ここで、Pjは空間座標x,yを含む多項式で、
は未定係数である。mは展開に用いた項数
で、例えば
aj
=3
m の場合は次式のように選ばれ る。
] , , 1 [ )}
, (
{P x y T = x y （４）
MLSMでは次式で定義される評価関数J を最
小化させるように未定係数ajを決定する。
∑ − −
= n
I
I h
w I
J ({x} {x })[φ ({x}) φ ]2
∑ − −
= n
I I
T
I P
w({x} {x })[{ ({x})} {a({x})} φ ]2
（５）
ここで{x}=(x,y)Tであり、n、{xI}はそ れぞれ評価点 の近傍に位置する節点数、節 点座標である。
} {x
φIは近似関数φhの節点I での
値であり、w({x}−{xI})は評価点{x}への節 点I の寄与度を重みづける関数で、以下「重み 関数」と記す。
これらの 個の節点での値が、評価点での値 を内挿するために用いられる。その影響の大き さを表現する重み関数は節点上で１、節点から の距離が増加するにつれて零に減衰するよう な滑らかな連続関数であるものとする。
n
式（５）の未定係数{a}についての停留条件よ
Analysis of Elastic Body by Meshless Finite Element Method
Kojiro OKAWA and Nobuyosi TOSAKA
り、 が求まり、さらにそれを式（3）に代 入することによって次式を得る。
} {a
} ){
( })
({ I I
h N φ
φ x = x （６）
ただし、
})]
({
[ })]
({
[ })}
({
{ )
( T I 1 I
NI x = p x A x − B x
（７）
ここで[A]、 はそれぞれm行 列、 行
列のマトリックスで、節点座標 ]
[B m m n
}
{xI から計算
できる。
式（６）によって、任意の評価点 での近 似間数
} {x
}) ({x
φh は、その点の近傍の節点値{φI} を用いて表現されることになる。評価点まわり の影響半径と重み関数を適切に選ぶことによ り、式（７）を内挿関数として用いることが可 能となる。
本手法では、影響半径を節点密度に依存しな い合理的な決定方法を用いる。すなわち、影響 半径 を積分評価点とその点を囲む最小の 多角形を構成する節点との距離の最大値 を基準として次式のように決定する。
dm
dmax
dmax
k
dm = × （８）
ここに、無次元数 をスケールファクタと称す る。 は、1.1から2.0程度の値である。
k k
また、重み関数として次式で示す四次のスプ ライン関数を用いる。
4 3
2 8.0( ) 3.0( ) )
( 0 . 6 0 . 1 ) (
I I
I I
r r
r r
w = − ρ + ρ − ρ
) 0
( ≤r≤ρI （９）
ここで、r、ρIはそれぞれ評価点 からの
距離、重み関数の台の半径である。また、
} {x })
{ }
({ I
w x − x を簡単のためwI(r)と表した。
2.3 数値積分
弱形式表示された基礎方程式を離散化する ためには、領域積分操作が必要である。節点情 報のみからこの操作を実施するためには節点 積分のルールを設定する必要がある。
本論では、図 1 に示すバケットを利用して節 点毎に領域積分のための重みを計算しておく。
ここで、バケットとは節点探索を効率的に行う ために用いられる解析領域被うように与えら れているものである。
具体的には、各バケットの面積をその中心か らの近傍の影響範囲内にある節点に分配する。
面積の分配率はバケットの中心と節点との距 離から計算される重み関数の値から決定する。
すなわち、バケット内のサンプル点iの、そ の近傍に取られた影響範囲内にある節点 へ のサンプル点の面積 の寄与分 を次式で計 算する。
k Ai aik
∑
== n
j ij
i ik
ik w A w
a
1
/ （10）
ここに は影響範囲内にある総節点数、 は 積分点iの座標と節点 との座標から計算さ
れる距離の関数となる重み関数である。 は 内挿関数を作成する際にMLSMで用いられた 重み関数と同様なものを用いることができる。
また、影響半径の大きさは式（8）で示した方 法で決定される。
n wik
k
wik
式（10）で与えられる面積 を、節点ごと に総和をとりこれを と表し、節点での重み とする。この、 を用いて解析領域の面積は 数値積分のための重みとして、節点上に分配さ れることになる。
aik
Wk
Wk
ここで、解析領域内部Vでの被積分関数
の面積分を次式のように近似評価する。
) , (x y F
∫∫ ∑
≅ = V
NP
I F xI yI WI
dxdy y x F
1
) , ( )
,
( （11）
ここに、NPは解析領域内部に分布する節点 の総数、xI、yIは節点I の座標である。
Bucket
Ai aik ｋ
dm:Radius of influence domain
Bucket center Node(reffered)
Node
Distribution of bucket area
図１ 領域積分のための節点重み なお、解析領域の面積 は次式を用いて近 似的に計算されることになる。
total
A
∑
== N
k k
total W
A
1
（12）
もし、別の方法で解析領域の面積が評価可能 な が既知であれば、式（12）を用いて節 点での重み の妥当性を検証することがで きる。
total
A
Wk
2.3 節点剛性評価
本手法では、節点単位での剛性方程式の成分 の寄与を計算する。すなわち、式（1）の左辺 第１項を、式（11）にならい次式のように近似 的に評価する。
∫∫ ∑
= V ≅
NP
I I I
T I
T dV W
1
} { } { }
{ }
{δε σ δε σ （13）
ここに{σI}、{εI}は節点I における応力、ひ
ずみを表すものとする。
節点I における変位 は次式のように表 すことができる。
} {uI
} ]{
[ }
{uI = NI UI （14）
ここに[NI]は式（7）のように節点I で局所的 に定義される内挿関数、{UI}は節点I の近傍 に位置する節点での変位成分を表す。
式（14）の微分を用いて、節点I でのひずみ
}
{εI は次式のように表される。
} ]{
[ }
{εI = BI UI （15）
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
=
x n y n x
y
y n y
x n x
I
N N N
N
N N
N N
, , ,
1 , 1
, ,
1
, ,
1
0 0
0 0
B ・・・ （16）
線形弾性体の応力―ひずみ関係は次式のよ うに表される。
} ]{
[ }
{σI = D εI （17）
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
− −
=
2 ) 1 0 ( 0
0 1
0 1
1 2 ν ν
ν ν
D E （18）
ここで、Eはヤング率、ν はポアソン比を表
す。
式 (14),(15),(17),を式(13)の右辺に代入 して整理して次式を得る。
∫∫ ∑
= V ≅
NP
I I I
T I
T dV U k U
1
} ]{
[ } { }
{ }
{δε σ δ （19）
ここに、 は節点単位で局所的に定義される 剛性成分で次のように表される。
] [kI
I T
I
I B D BW
k ] [ ] [ ][ ]
[ = （20）
すなわち、剛性マトリックスの成分は節点単位 で計算され、これを重ね合わせることによって 系全体の剛性マトリックスが組み立てられる。
式（１）の左辺第2項、第三項で表される表 面力、物体力も同様な方法で離散化することに より、式（１）は次式のように表される。
0 }) { } ]{
([
}
{δU T KG U − F = （21）
∏
== NP
k I
G k
K
1
]
[ （22）
また、 は全体節点変位ベクトル、 は全 体節点荷重ベクトルである。式（22）の
}
{U {F}
∏は 全体化の操作を表す。式（21）から得られる全 体剛性方程式を解くことにより、未知変位を求 めることができる。
3 数値計算例
本手法の計算精度を検討するために、2次元 弾性問題を上記の手法を用いて解く。その結果 を厳密解またはFEMによる解析結果と比較する。
3.1 片もち梁の解析
図2に示す片もち梁の自由端に荷重を与え変 位を求める。節点分布を図3に示す。スケール ファクタ を1.4、内挿関数を決定するための MLSMでの重み関数として4次のスプライン、バ ケットの一辺の長さを節点間距離程度とした。
k
この問題の、厳正解は次式で与えられる。
} ) 3 ( ) 5 4 4( ) 1 ( 3 6 {
4 )}
)( 1 2 ( ) 3 6 6 {(
2 2
2
2 2
x x L x D x
L EI y u P
D y x
x EI L
u Py
y x
− + +
+
−
=
− + +
− −
=
ν ν
ν
（23）
ここで、Iは断面2次モーメントである。各節 点の変位について厳密解と比較して表1,2に示 す。
D=1mm
L=10mm y
x Perfectly fixed
P
20600kgfmm2
E=
kgf P=1.0
3 .
=0 ν
図２ 片もち梁
n u m b er o f n o d es 36 9
図3 節点分布図
表1 変位の比較（x方向）
-5.E-04 -4.E-04 -3.E-04 -2.E-04 -1.E-04 0.E+00 1.E-04 2.E-04 3.E-04 4.E-04 5.E-04
-6.E-04 -4.E-04 -2.E-04 0.E+00 2.E-04 4.E-04 6.E-04
Meshless Method[mm]
exact[mm]
表2 変位の比較（y方向）
0.E+00 5.E-04 1.E-03 2.E-03 2.E-03 3.E-03 3.E-03
0.E+00 5.E-04 1.E-03 2.E-03 2.E-03 3.E-03 3.E-03 Meshless Method[mm]
exact[mm]
3.2 正方形モデルの解析
物性値同定逆問題では、図４に示すモデルを 用いる。スケールファクタ 、重み関数、バケ ットの一辺の長さは3.1と同様とし解析を行う。
表３にモデルの上端における変位をFEMと比 較して示す。
k
number of nodes 36
20600kgfmm2
E
kgf P 1.0
3 . ν 0
P
L 10m
図4 正方形モデル
表2 モデル上端における変位 FEM Meshless Method
8.714E-05 4.001E-05 9.692E-05 3.821E-05 7.057E-05 1.487E-05 7.585E-05 1.646E-05 6.091E-05 2.998E-06 6.829E-05 7.538E-06 5.565E-05 -5.056E-06 5.785E-05 -1.005E-05 5.362E-05 -1.289E-05 5.677E-05 -1.599E-05 5.336E-05 -2.266E-05 5.063E-05 -3.230E-05
uy ux
ux uy
4 おわりに
本論では、メッシュレス有限要素法を2次元 弾性体に適用しその計算精度について比較し た。数値計算例により、本手法を厳密解が得ら れている問題に適用し解がほぼ一致する結果 を得られた。さらに、厳密解が得られない問題 にも、同じ節点配置を用いたFEMによる解析結 果と比較することにより同等の結果が得られ ることを示した。
計算時間に関しては、入力情報として要素―
節点コネクティビティを必要としない本手法 は、各節点での重みを計算するプロセスの分だ け、FEMと比較して多くの計算時間を必要とし た。
今後は、2次元弾性体において静的加重によ る変位を観測値とし、ヤング率やポアソン比な どの弾性定数を同定する物性値同定逆問題の 解析に本手法を援用して適用を検討してゆく。
参考文献
[1]Belytschko,T.,Lu,Y.Y.andGu,L.,Element- Free Galerkin Methods,Int. J. Num.
Methods Eng., 37(1994),PP239-256 [2]奥田洋司・長嶋利夫・矢川元基、エレメン
トフリーガラーキン法に関する基礎的検討 (第1報)常微分方程式への適用、機会学会論 文集 (A編)61―590（1995）PP2302−2308 [3]奥田洋司,長嶋利夫、野口裕久、わかるエレ
メントフリー法、日本機械学会講習会教材、
No. 98-67,（1998） pp1-78
[4]H.Okuda,Y.Satoh and O.Hazama,Digital Analysis Procedure Using Image ProcessIng and Element-Free Galerkin Method,Proc.ICES98,MODELING AND SIMU- LATION BASED ENGINEERING, Vo1Ⅱ, (1998), PP1825-1830
