SPH 法による傾斜機能材料の動的応答解析*

304
日本機械学会論文集(A 編)
79 巻 799 号 (2013-3)
原著論文 No.2012-JAR-0624
SPH 法による傾斜機能材料の動的応答解析*
戎 圭明
∗1
, 木須 博行
∗1
Investigation of the Dynamic Response of Functionally Graded Materials Using Smoothed
Particle Hydrodynamics
∗1
Guiming RONG
∗1
and Hiroyuki KISU
Department of Mechanical Engineering, School of Engineering, Nagasaki University
1-14 Bunkyo-machi, Nagasaki 852, Japan
In the present study, the problem of functionally graded materials (FGMs) under a stress pulse is analyzed based on
smoothed particle hydrodynamics (SPH) using the formulation in which the deviatoric stress and the continuity equation
are used. First, the formulation of SPH for this problem is described and several benchmark calculations for the primary
wave are performed compared to one-dimensional analytical solutions to show that the present formulation can be used
to analyze the dynamic response with considerable accuracy at the initial stage. The behavior of FGMs subjected to a
stress pulse is then investigated for several cases, including individual influence from distribution of density or Young’s
modulus to the stress induced in the material, etc.. In the two-dimensional case where the system permits free deformation
perpendicular to the direction of the external force, the primary wave as well as the secondary wave are both checked. At
that time, induced stress waves due to the free surface are observed.
Key Words : Smoothed Particle Hydrodynamics (SPH), Functionally Graded Materials (FGMs), Impact Loading,
Dynamic Response, Primary Wave, Secondary Wave
1. は
じ
め
に
(1)
傾斜機能材料(FGM)は元々,超高温断熱材として開発された が,近年,運動部品の表面の衝撃緩衝層とし
ても応用されるようになった.そのため,熱応力解析だけではなく,衝撃に対する性能も研究されている.さら
(2)
に衝撃低減のため FGM の製造について傾斜発泡材料(FGFM) の製造も提案されている.これらの開発にあたっ
ての課題は,FGM についてはたとえば衝撃に対して母材と FGM との接合面においてき裂や破壊等が発生しやす
いこと,FGFM については,材料の分布をどのようにすれば最も衝撃を吸収できるようにできるか,などがある.
このように,FGM に対しても FGFM 対しても衝撃解析が不可欠のものとなり,多くの研究者によって研究が始め
(3)
られている.例えば,Chiu と Erdogan はラプラス変換を使って FGM 板を1次元的に解析し,両端の境界条件に
(4)
ついて自由-自由と固定-自由の二種類の組合せの下での衝撃波の理論解を導いた.Bruck は,FGM の成分をうま
(5)
く分布させ,応力波を制御できるようなデザインを研究した.Jabbari ら は厚い FGM 円筒の熱応力と衝撃応力に
(6)
対する理論解を与えた.Yas らは layer-wise finite element formulation を使って FGM 材料の弾性波の伝播を解析
(7)
した.Gupta らは FGFM に関して変動荷重条件下での応答の様子を解析した.しかし,これらの研究は未だ一次
元や軸対称問題に留まっており,一般的な2次元問題としてはまだ扱われていない.さらに3次元解析も必要にな
ると思われる.
(8)
このような観点から,著者らはこの種の問題に対する一般的な数値解析法の確立が必要と考え,SPH 法 に着目
した.SPH 法は要素分割が不要なので,特に大変形など有限要素法などで処理しにくい場合は非常に効率的であ
*
原稿受付 2012 年 9 月 10 日
正員,長崎大学大学院工学研究科(〒852-8521 長崎県長崎市文教町 1-14)
E-mail: [email protected]
*1
― 91 ―
© 2013 The Japan Society of Mechanical Engineers
305
SPH 法による傾斜機能材料の動的応答解析
る.また,非線形問題などにおいてもモデルの作成が簡単であるためスムーズに対応できる.SPH 法による弾性
(9)
衝撃問題解析には酒井ら が疎密波(P 波)に関して解析を行っている.ただし,重み関数や影響半径などの基礎
的な検討の段階であり,解析時間も反射波が生じない範囲にとどまっている.さらに,Blanc と Pastor
(11)
法を土壌の疎密波による衝撃問題に適用している.著者らは傾斜機能材料の熱衝撃問題
(10)
は SPH
を SPH 法で解析し,そ
の有効性を確認している.本研究ではこれをさらに FGM の衝撃問題解析に拡張して定式化を行い,二次元問題の
解析プログラムを開発し,疎密波とせん断波(S波)とも観察した.ここでは,2つの興味ある課題,すなわち,
動的応答に及ぼす FGM の密度分布とヤング率分布の個別的な影響,および境界条件の違いによる影響,主に自由
表面による弾性応力波が発生して応力波が重畳する現象について報告する.
2. SPH 法による固体材料動的応答問題の定式化
2·1 弾性体の動的応答問題の支配方程式
弾性体の動的応答問題の支配方程式として以下を適用する.
∂ui 1 ∂σi j
=
+ gi ,
∂t
ρ ∂x j
(1)
ここで ui は単位質量部分の速度の i 番目成分, t は時間,ρ は密度,σi j は応力テンソル, x j は単位質量部分の位
置ベクトル, gi は単位体積当たりの体積力ベクトルである.応力テンソルとひずみ速度 ε˙ i j の関係は
(
)
1 ∂ui ∂u j
ij
e ij
ij
σ =D ε ,
ε˙ =
+
2 ∂x j ∂xi
で,De は弾性マトリックスである.Blanc と Pastor
アルゴリズムと応力ノード
(10)
(2)
の研究によると,以上のような定式化では Runge-Kutta-Taylor
(12)(13)
を用いなければ張力不安定問題(tensile instability)が生じて計算不能に陥ってし
まうようである.そこで,本研究では σi j を直接に計算に用いるのではなく,以下の偏差応力を用いた
σi j = −pδi j + S i j ,
(14)
.
(3)
上式において p は圧力,S i j は偏差応力である.Hooke’s 定理を適用すると,ヤング率を E(x) として偏差応力速度
S˙ i j は以下の式で表わされる.
(
)
E(x) i j 1 ik k j
ij
˙
S =
ε˙ − δ ε˙ ,
(4)
1+ν
3
ここに ν はポアソン比である.また,質量保存則と状態方程式
(15)(16)
により,以下の式が成り立つ.
ρ˙ = −ρ∇ · u,
p = c20 (ρ − ρ0 ),
(5)
c20 = K/ρ0 ,
K=
E
.
3(1 − 2ν)
(6)
ここに,ρ0 は基準密度,K は体積弾性率である.ポアソン比 ν = 一定とする.
以上の定式化をこの種の問題に適用するのは本研究が初めてと思われるが,これによれば上述した張力不安定
問題は生じない.一番の大きな違いは,連続の式 (5) を用いて体積のわずかな変化まで計算に組み込んでいること
であり,詳細はさらに調べたいが,これが大きな寄与を為しているのではないかと考えている.
2·2 SPH による離散化
本研究では式 (1) を SPH で離散化した式として以下を用いた.


N
∑
 σi j σi j
 ∂W(rab , h)
a
b
i
j
i
mb  2 + 2 + Πab δ 
u˙ =
+ gi ,
j
ρ
ρ
∂xa
a
b
b=1
― 92 ―
(7)
© 2013 The Japan Society of Mechanical Engineers
306
SPH 法による傾斜機能材料の動的応答解析
Table 1 Parameters used in analysis for a FGM slab
ρ0 (Kg/m3 ) E0 (GPa)
ν
l (mm)
a
m
5331
166.198
0.3
5
-0.12354
-1.8866
n
-3.8866
ここに a は注目粒子の番号,b はそれ以外の粒子を表わす.W(rab , h) は核関数であり,半径 h 以内を影響領域とす
る.rab は粒子 a と b の距離で,mb は粒子 b の質量, ρb は粒子 b の密度,N は粒子の総数である.また,Πab は
Monaghan(8)による人工粘性項であり,二つの成分を含んでいる.


−αc0 ϕab + βϕ2ab




, (ua − ub ) · rab < 0
h(ua − ub ) · rab

ρ¯ ab
Πab = 
,
ϕab = 2
,


rab + (0.1h)2


0,
(ua − ub ) · rab ≥ 0
1
ρ¯ ab = (ρa + ρb ).
2
(8)
係数 (α, β) の選択について,普通 (1, 1) の方が多いが,予備計算の末,本研究では (α, β)=(0.1, 0) とした.
各粒子の速度勾配は次のように離散化される.
∂uia
j
∂xa
=−
N
∑
mb
ρb
b=1
(uia − uib )
∂W(rab , h)
j
∂xa
.
(9)
また,質量保存則 (5) を離散化すると以下の式となる.
ρ˙ a = −
N
∑
mb (ua − ub )∇W(rab , h).
(10)
b=1
2·3 理論解との比較検証
(3)
文献 には1次元問題の理論解が示されている.この問題を開発した2次元問題用プログラムで解くことによ
り,以上の定式化とプログラムの検証を行う.以下の解析において計算モデルはすべて z 方向に無限に連続してい
る想定で,2次元平面ひずみ条件を設定した.
(3)
文献 に示された1次元解析モデルの一つは,図 1 に示すように固定境界と自由境界に挟まれた無限に長い FGM
領域層が想定されている.図のように厚さ方向に x 軸を設定すると,一つの FGM の例ではジルコニアとニッケル
からなる成分は次のような分布となっている.通常はヤング率の分布と密度の分布は材料制作上の都合から相互
に連動するため,これらの分布関数は互いに独立ではない.
)m
)n
( x
( x
E ′ (x) = E0′ a + 1 ,
ρ(x) = ρ0 a + 1 ,
l
l
where
E′ =
E(1 − ν)
.
(1 + ν)(1 − 2ν)
(11)
1.8
E/E0
ρ/ρ0
1.6
x
1.4
x
1.2
σ0 f (t)
1
Y
y
0.8
l
0
ZrO2
σ0 f (t)
0.2
0.4
0.6
0.8
x/l
1
Ni
l
Fig. 3 Calculation model
Fig. 1 Model of an FGM slab
Fig. 2
Variation of density and elastic
modulus in model calculation
― 93 ―
© 2013 The Japan Society of Mechanical Engineers
307
SPH 法による傾斜機能材料の動的応答解析
1.5
2.0
Theo.
With art.
No art.
1.0
Theo.
With art.
No art.
1.0
0.0
0
1
2
3
4
5
6
σxx/σ0
σxx/σ0
0.5
0
0
1
2
3
4
5
6
-0.5
-1.0
-1.0
-2.0
-1.5
Time(µs)
Time(µs)
Fig. 4
Variation of σ with time at x = 0.
Results
Fig. 5
obtained by SPH and theoretical values.
Variation of σ with time at x = l/2. Results
obtained by SPH and theoretical values.
ここに,l は傾斜機能層の厚さで,各パラメータの数値は表 1 に示す通りである.境界条件は x = 0(母材と
の接合面)を固定し, x = l の自由境界に時刻 t = 0 において外部から σ0 f (t) のパルス荷重を与える.本解析では
f (t) =H(t)H(−t + t0 ) とすると,t0 = 0.2µs, σ0 = 1Mpa とした.ただし,H(t) は Heaviside 関数である.
(3)
本研究ではこの問題を図 3 のような二次元モデルを使って解析し,破線上(y = Y/2) の解析結果を文献 による
一次元の理論解と比べる.離散化は領域を 160×80 個の粒子で均一に分布させ,時間積分法としては蛙とび積分法
(17)
を用いた.影響半径については文献
(16)
を参考に h/rab = 2 とした.
計算の順序としては次のようになる.まず式 (4) と式 (10) を使って速度ベクトルから偏差応力と密度を評価し,
その密度を使って 式 (6) から圧力を得て,式 (3) で応力を求めることとなる.
解析の結果得られた数値解を図 4 と図 5 に示す.図はそれぞれ x = 0 と x = l/2(いずれも y = Y/2)の点における
σ xx の時間変化であり,時刻 6µs までの結果を示す.図中,実線が文献(3)で示された理論解であり,計算結果は人
工粘性を使った場合と使わなかった場合の両方が比較されている.もし FGM 材料でなく均質材であったら,x = 0
における応力振幅は |σ xx /σ0 | = 2 となるべきものである.図 4 の理論解は |σ xx /σ0 | = 2 を若干下回っており,この
タイプの FGM による応力緩和がわずかとはいえ実現されていることになる.
ところで,2つの図からわかるように,人工粘性を使った場合と使わない場合とでその差はあまり大きくなく,
特に周期についてはほとんど影響されないようである.一方,振幅に関しては人工粘性を使わないデータの方が
やや大きな値を示す.これらのことから,人工粘性は使っても使わなくても大差がないと思われる.しかし,一般
的に言って,人工粘性は非物理的(数値的)な振動を抑制する効果が大きいために導入される.実際,FGM の物
性値が x = 0 と x = l の境界間で差が激しい時は,人工粘性なしの場合に数値が発散するケースもあった.このよ
うなことを総合的に考えて,本研究では人工粘性を導入することとした.
ただし,解析結果には問題点も見られる.すなわち,周期に関してはほぼ十分な精度の結果を示すものの,振
幅については時間が経過すると減衰しているようである.解析対象は弾性体であるから,言うまでもなく減衰は
生じないので,これは何らかの数値解析上の問題だと思われる.しかし,衝撃問題においては初期段階の応答に
興味がある場合が多く,長時間の解析でなければ,十分な実用性を有している.今後,長時間の解析にも堪えられ
るように,計算手法を開発していきたい.
3. SPH による傾斜機能材料の動的応答の挙動の解析
(3)
以上の基礎的検討を経て,本研究では傾斜機能材料の動的応答の挙動の解析を行った.Chiu ら は,傾斜機能
材料におけるヤング率と密度の分布が,弾性衝撃によって誘起される応力の大きさにどのように影響を与えるの
かについて,興味深い知見を得ている.例えば,図 1 において,固定側 x = 0 のヤング率と密度が荷重負荷側 x = l
のヤング率と密度より小さい場合の方が,その逆の場合よりも固定境界における応力の振幅が小さいことを示し
た.また,応力振幅の分布は E(x)ρ(x)/E(l)ρ(l) に依存した関数で表されると結論した.ただし,ヤング率と密度は,
材料形成の都合で互いに独立した分布にはできないものとされている.
― 94 ―
© 2013 The Japan Society of Mechanical Engineers
308
SPH 法による傾斜機能材料の動的応答解析
2.25
2
1.75
1.5
E/E0 or ρ/ρ0
1.25
1
0.75
0
0.2
0.4
0.6
0.8
1
x/l
Fig. 6 Variation of density and elastic modulus
本研究ではこれらの成果を踏まえ,FGFM への応用も考慮して,
(1)ヤング率 E(x) と密度 ρ(x) の独立した影響
を明らかにすること,および(2)固定壁(母材との接合面)以外の境界条件による影響を明らかにすることを
目的とした.
(1)については,ヤング率と密度は材料の構成により決まるので互いに独立した分布にはできないように思わ
れるが,最近,例えば中空微粒子を分散させることにより密度を調節するようなことも可能
(18)
になっている.こ
こではこのような事情を考慮して問題を設定した.
また,
(2)については,これまで傾斜機能層が無限に長い状態を想定して来たのに対し,現実には有限の領域
で形成されることから,境界条件の違いによる動的応答への影響を明らかにするためのものである.
3·1 ヤング率と密度分布の個別の影響
再三述べたように,従来の研究は材料定数の分布に関しては E(x) と ρ(x) の分布は連動しており,例えば図 2 の
分布でも,E(x) と ρ(x) の個別の分布は独立ではない.しかし,本研究では今後のために,同一分布でも E(x) と
ρ(x) の変化を独立して設定し,それらの影響を見ていくことにする.具体的には下記のように設定した.
材料定数の分布曲線を図 6 のように設定し,個別の変化を次のように与える.
(
x)
f (x) = 0.950213 + 0.049787 exp 3.0485
l
(12)
ここに,材料定数の分布は 4 ケースを設定する.すなわち,1. E(x) = E0 f (x) と ρ = ρ0 (= const) のケース,
2. ρ(x) = ρ0 f (x) と E = E0 (= const) のケース,3. E(x) = E0 f (x) と ρ(x) = ρ0 f (x) のケースの3ケースで,
各ケースそれぞれの計算を経てこれらの個別の影響を調べるのである.また,誤差に関する考察のため,4.
E(x) = E0 (= const) と ρ = ρ0 (= const) のケースも計算した.
ただし,E0 = 196 GPa, ρ0 = 7, 850 (kg/m3 ) とした.結局, x = 0 と x = l における Eρ の値の比は,ケース 1 と 2
では2,ケース 3 では4となる.その他,衝撃荷重の負荷時間 t0 = 0.2µs,その大きさ σ0 = 1, 000 MPa とおいた.
計算モデルは図 1 と同じであり,計算時間は 5µs までとした.以上の解析結果を図 7 と図 8 に示す.材料定数が
一定の場合,生じる応力波のピーク値は2倍となるが,本解析によれば 2.11 となり,良い一致を示している.た
だし,一周期後では約 20%程度減衰し,2.3 節で示した数値的な減衰がここにも見られる.また,これらの図を見
ると,ヤング率と密度の分布が独立に変化した場合の影響は決して小さくないことがわかる.まず,ヤング率の
みを変化させたケース 1 の場合,応力波の伝播速度が速くなる.そして生じる応力波のピーク値は,ヤング率一
定のケース 2 の場合に比べて約 10% 大きくなった.密度のみを変化させたケース 2 の場合に応力波の伝播速度が
最も遅くなった.これらのことから,応力緩和に対しては密度を変化させる方がわずかながら効果が大きいので
はないかと思われる.そこでさらに,ヤング率を一定とし,x = 0 での値に対する x = l での密度の比を2から10
まで変化させて計算し,固定壁上での応力緩和程度を比較してみた.その結果を図 9 に示す.その比が10とな
る時,応力緩和はほぼ 50% に達することが分かった.
― 95 ―
© 2013 The Japan Society of Mechanical Engineers
309
SPH 法による傾斜機能材料の動的応答解析
2.0
8.0
1.5
1.0
4.0
σxx/σ0
σxx/σ0
5.0
0.0
-5.0
0
1
2
3
4
5
-1.5
-2.0
-2.5
1
-8.0
-1.2
Time(µs)
Variation of σ with time at x = 0. Results of: 1.
Fig. 7
0
-4.0
Vary both
Vary ρ
Vary E
Const.
-1.0
0.0
Fig. 8
2
3
Vary both
Vary ρ
Vary E
Const.
4
5
Time(µs)
Variation of σ with time at x = l/2. Results
of: 1. both ρ and E vary; 2. ρ varies while
both ρ and E vary; 2. ρ varies while E remains
E remains uniform; 3. E varies while ρ remains
uniform; 3. E varies while ρ remains uniform;
uniform; 4. both ρ and E are uniform.
4. both ρ and E are uniform.
0.6
Reduction rate(%)
0.5
0.4
0.3
0.2
0.1
0
2
Fig. 9
3
4
5
6
7
Ratio of ρ(l)/ρ(0)
8
9
10
Maximum stress reduction rate for cases
Fig. 10
when E=const. while ρl /ρ0 varies from 2
Distribution of σ x at the time of t = 0.6µs, by the
model of Fig. 3.
to 10.
3·2 有限領域効果の解析
これまでは図 1 のような FGM 層の長さを無限と考えた取り扱いであった.従って,計算モデルは図3これまで
の計算結果はすべて対称軸上のデータである.しかし,実際には図3のように2次元モデルとして有限境界部分
が存在する.その影響を見るため解析領域全体の σ x の分布を図 10 に示す.それによると,固定境界付近の狭い
範囲で σ x の値が大きくなる部分(赤い部分)が生じている.これはy方向に固定されている境界といえども,ご
く内部ではポアソン効果によって y 方向の変位を生じるので,それに伴い σ x も副次的に誘起されたものと考えら
れる.しかし,その影響範囲は固定境界付近であるために限定されている.これに対し,自由境界付近ではその
影響範囲は大きくなるものと思われ,次に検討する.
さらに,自由境界の影響を見るため,図 11 の解析モデルを設定し,二次元問題として衝撃問題を解析する.な
お,有限領域について,FGM 材料だけではなく,均質材も同様な効果があることを記したい.
前節までとの本質的違いは y = Y の辺の境界条件であり,ここでは応力フリーの自由境界条件を与える.これ
以外の計算に必要なパラメータは 2.3 節と同じであるが,矩形の縦の長さ Y に関しては,傾斜機能材料は厚さよ
り広さの方がはるかに大きいことを反映させ,領域は l × Y =5 mm×10 mm とした.その他の境界条件は図示の通
りである.また,パルス応力 σ0 = 1, 000 MPa,および作用時間 t0 = 0.2µs も同じとした.SPH ノードについては,
160×320 個を設定した.
図 12 と図 13 はそれぞれ,衝撃を受けた直後の t = 0.6µs の時点における σ x と τ xy の分布図である.これらの図
― 96 ―
© 2013 The Japan Society of Mechanical Engineers
SPH 法による傾斜機能材料の動的応答解析
310
x
Y
y
σ0(t)
l
Fig. 11 A 2-D model
Fig. 12 Distribution of σ x at the time of t = 0.6µs.
Fig. 13 Distribution of τ xy at the time of t = 0.6µs.
Fig. 14 Distribution of σy at the time of t = 0.6µs.
で注目されるのは,矩形領域の右下の角点( x = l, y = Y )を起点とした斜めに進行する応力波が認められることで
ある.すなわち,入力されたパルス状外力によって直接生起された応力波(これを主波と呼ぶ)以外に,自由表
面境界により誘起される斜めに進行する膨張波とせん断応力波が観察される.
これは境界が自由表面であるための二次元波動伝播の効果である.すなわち,衝撃応力によって x 方向に変形
した材料はポアソン効果によって y 方向の変形を誘起する.その y 方向の変形は主波の進行に伴って x 軸方向に
次々と伝播し,未到達部分との境界でせん断応力成分 τ xy を誘起し,すると同時に新たな σ x と σy を二次的に生
起するものと考えられる.実際,同時刻において,σy も図 14 に示すように分布している.
このような二次的応力波は主波の進行方向と垂直な法線を持つ自由境界辺でのみ生じ,主波が到達した時点で
発生するため,斜めに進行する応力波となることが理解できる.この斜めに進行する応力波の角度はポアソン比
に依存するが,当然ながら均質材料ではないことの影響を受ける.
― 97 ―
© 2013 The Japan Society of Mechanical Engineers
SPH 法による傾斜機能材料の動的応答解析
311
Fig. 15 Distribution of σ x at t = 0.95µs. The free-deformation wave overlaps the main stress wave.
このシミュレーションでは t = 0.95µs において斜めに進行する膨張波の起点が固定境界に到達する.その時点に
おける σ x の分布を図 15 に示す.この図の左下の角点 (x = 0, y = Y) のあたりにおいて,応力が内部の固定境界と
比べて約 16% 大きくなったが,特異性のある角点の影響とみられる.
4. 結
言
本研究では SPH 法による傾斜機能材料の動的な応答解析について報告した.主要な結論を列記する.
1. 偏差応力と連続の式を用いた動的問題の定式化を提案した.これにより張力不安定問題を避けることが出来た.
2. 長時間の解析では減衰のしすぎとなって未だ不十分であり,さらに研究の余地があるものの,衝撃初期の問
題の解析には十分な精度の解析が可能となった.また,疎密波だけではなく,せん断波や自由表面により誘
起された膨張波も観察できた.
3. 傾斜機能材料における密度の分布とヤング率の分布の独立した影響を調べた.母材との接合面である固定境
界(本研究では x = 0 の辺)における値より,衝撃荷重負荷境界(本研究では x = l の辺)の値を大きくした
場合,ヤング率一定で密度を傾斜分布させる方が,その逆よりも応力波の伝播速度が遅くなること,応力緩
和も大きくなることがわかった.密度の比が 10 になると応力緩和の程度が 50% にも達する結果が得られた.
文
献
(1) M. Yamanouchi, M. Koizumi, T. Hirai and I. Shiota (editors) FGM ’90, Proceedings of the lst International Symposium on
Functionally Gradient Materials, Tokyo, Japan: FGM Forum., (1990).
(2) El-Hadek M. M. and Tippur H. V., “Dynamic fracture parameters and constraint effects in Functionary Graded Syntactic Epoxy
Foams”, Int. J. of Solids and Structures, Vol. 40 (2003), pp. 1885-1906.
(3) T. C. Chiu and F. Erdogan, “One-dimensional wave propagation in a functionally graded elastic medium”, J. Sound & Vibration,
Vol. 222, No. 3 (1999), pp. 453-487.
(4) Hugh A. Bruck, “A one-dimensional model for designing functionally graded materials to manage stress waves”, Int. J. of Solids
and Structures, Vol. 37 (2000), pp. 6383-6395.
(5) M. Jabbari, S. Sohrabpour and M. R. Eslami, “Mechanical and thermal stresses in a Functionally graded hollow cylinder due to
radially symmetric loads”, Int. J. Pressure Vessels and Piping, Vol. 79 (2002), pp. 493-497.
― 98 ―
© 2013 The Japan Society of Mechanical Engineers
SPH 法による傾斜機能材料の動的応答解析
312
(6) M. H. Yas, M. Shakeri, M. Heshmati and S. Mohammadi, “Layer-wise finite element analysis of functionally graded cylindrical
shell under dynamic load”, Journal of Mechanical Science and Technology, Vol. 25, No. 3 (2011), pp. 597-604.
(7) Gupta N., Gupta S. K. and Mucller B.J., “Analysis of a Functionally Graded Particular Composite under flexural loading
conditions”, Material Science and Engineering, Part A, Vol. 485 (2008), pp. 439-447.
(8) J.J. Monaghan, “Smoothed particle hydrodynamic”, Annu. Rev. Astron. Astrophys., Vol. 30 (1992), pp. 543-574.
(9) 酒井譲,山下彰彦,“SPH理論に基づく粒子法による構造解析の基礎的検討 ” ,日本機械学会論文集 A 編, Vol. 67, No. 659
(2001), pp. 1093-1102.
(10) T. Blanc and M. Pastor, “A stabilized Fractional Step, Runge-Kutta Taylor SPH algorithm for coupled problems in geomechanics”,
Comput. Methods Appl. Mech. Engrg., Vol. 221-222 (2012), pp. 41-53.
(11) 戎圭明,木須博行,
“ SPH 法による傾斜機能材料の非線形熱解析 ”,日本機械学会論文集 A 編, Vol. 76,No. 772 (2010),
pp. 1756-1763.
(12) C. T. Dyka, P. W. Randles and R. P. Ingel, “Stress points for tension instability in SPH”, Int. J. for Numerical methods in
engineering, Vol. 40 (1997), pp. 2325-2341.
(13) P. W. Randles and L. D. Libersky, “Normalized SPH with stress points”, Int. J. for Numerical methods in engineering, Vol. 48
(2000), pp. 1445-1462.
(14) L.D. Libersky, “High strain Lagrangian hydrodynamics”, J. Comput. Phys., Vol. 109 (1993), pp. 67-75.
(15) J. P. Gray, J. J. Monaghan and R. P. Swift, “SPH elastic dynamics”, Comput. Methods Appl. Mech. Eng., Vol. 190, No. 49-50
(2001), pp. 6641-6662.
(16) J.J. Monaghan, “SPH without a tensile instability”, J. Comput. Phys., Vol. 159 (2000), pp. 290-311.
(17) G. R. Liu and M. B. Liu, Smoothed Particle Hydrodynamics–A mesh free particle method (2003), p. 141, World Scientific,.
(18) Pradeep K. Rohatgi, Takuya Matsunaga and Nikhil Gupta, “Compressive and ultrasonic properties of polyester/fly ash
composites”, J Mater Sci, Vol. 44 (2009), pp.1485-1493.
― 99 ―
© 2013 The Japan Society of Mechanical Engineers