CIP 法による 3 次元弾性体の波動伝播解析 WAVE PROPAGATION

49
法政大学情報メディア教育研究センター研究報告 Vol.29
2015 年
CIP 法による 3 次元弾性体の波動伝播解析
WAVE PROPAGATION ANALYSIS OF
THE THREE DIMENSIONAL ELASTIC BODY BY CIP METHOD
柳川智隆１）
松下周平１）
吉田長行 2）
Tomotaka Yanagawa, Syuhei Matsushita, Nagayuki Yoshida
１）
法政大学大学院デザイン工学研究科建築学専攻
法政大学デザイン工学部建築学科
２）
When we calculate the earthquake response of the building, the effect of Dynamic
Soil-Structure Interaction is considered. Generally the finite element method (FEM) is used to
analyze the dynamic behavior of the soil. This method is well-known to be effective in dealing
with the nonlinear behavior. The analytical domain of the FEM has to be closed with the
boundaries. But the dissipated wave is propagated over the boundaries and not reflected.
Therefore some numerical technique is required to realize the open boundary. The advection
equation can decompose the wave of the field into the incoming wave and the out -going one.
This property enables us to make the open boundary processing for the FEM analysis. The
CIP(Constrained Interpolation Profile Scheme) method is a typical procedure with high precision for
the advection equation. This paper investigates the basic characteristics of the CIP method for
the transmission of the out-going wave against the three dimensional elastic field.
Keywords : Elastic wave , CIP method, Three dimensional field, Open boundary
1．はじめに
近年、地盤の非線形な動的挙動が活発に研究され
ている。非線形問題を扱う場合、有限要素法が有効
かつ柔軟な手法であることはよく知られている [1]。
有限要素法によって無限あるいは半無限地盤を離散
化するためには開境界処理を必要とする。境界処理
法として、境界にダッシュポットを設ける粘性境界
が代表的であるが完全な波動透過は実現されていな
い[2],[3]。
そこで、本論では CIP 法(Constrained Interpolation
Profile Scheme)を用いた新しい境界処理法の可能性
を検討する。CIP 法は矢部ら[4],[5]によって考案された
原稿受付 2015 年 3 月 9 日
発行
2015 年 4 月 1 日
Copyright © 2015 Hosei University
移流方程式を解く高精度差分法である。
本論の目的は新しい境界処理法のための基礎研究
として CIP 法を用いた 3 次元弾性体波動伝播解析を
行い、その特性を把握することにある。
2．弾性体の移流方程式とその CIP 解法
弾性体の移流方程式は波動方程式から導くことが
できる。
2.1 弾性体の波動方程式
次のように𝑛(= 1,2,3)次元弾性体の振動方程式は表
せる。
50

1

un   n  n , un  un
t

t
(1)
 n   Dn  n   Dn n un
(2)
T
ここで、変位を{𝑢}、応力を{𝜎}、ひずみを{𝜀}、偏微
分演算子をマトリクス表記したものを [] 、密度
を  、時間を t、ヤング率 E とポアソン比 からな
る応力ひずみマトリクスを[𝐷]とする。
式(1)に式(2)を代入したものが波動方程式(3)となる。
2
1
T
u   n  D n  n un
2  n

t
(3)
2.2 弾性体の移流方程式
式(3)は[0]を任意のゼロマトリクス、[𝐼]を任意の単
位マトリクスとおくと、まとめて次のように表せる。

 u 
 [0]
  
t  n
[ D]n

[ I ]n  []nT

  [0]
[0]  
1
[0]  u 
   (4)
[]n   n
簡易に表すと、

F   AQF
t



  q y    qz 
x   y
z
(5)
(6)
これを用いて式(5)を書き換えると、




F   Ax  F   Ay  F   Az  F  (7)
t
x
y
z
ここで、
 As    Aqs  , s  x, y, z
fi n1  ai 3  bi 2  gin  fi n
gin1  3ai 2  2bi  gin
ai 
bi 
ここで、
Q    qx 


(12)
 f s   [ s ]  f s  , s  x, y, z
t
s
2.3 移流方程式の CIP 解法
CIP 法は、矢部らによって考案された高精度差分
法あり、数値拡散が少ない高精度な移流方程式解法
である。移流方程式


f  x, t   u
f  x, t   0
(13)
t
x
は移流関数 f が x 軸を速度 u で伝播していることを
表す。この式を CIP 法で解くことは次式を用いるこ
ととなる。𝑛ステップ時における求める i 点とその隣
接点 iup 点の移流値 f とその空間微分値 g の情報か
ら∆𝑡秒後の𝑛 + 1ステップでの移流値と空間微分値
を導くことができる。
(8)
gi  giup
2

D
3( fiup  fi )
2( fi  fiup )
D3
2 gi  giup
(14)
(15)
(16)
(17)

D
D2
ここで、
  ut である。また、∆𝑠を節点間距離
とすると、𝑢 ≥ 0のとき𝐷 = −∆𝑠，𝑖𝑢𝑝 = 𝑖 − 1、𝑢 ≤ 0
のとき𝐷 = +∆𝑠，𝑖𝑢𝑝 = 𝑖 + 1である。
2.4 分離解法
分離解法は、多くの分野において使用される多次
元解法であり、1 次元問題に落とし込むことで、容
易に解くことが可能な手法である。
CIP 法の分離解法には M 型 CIP と C 型 CIP の 2
種類がある。
方向分離すると、


(9)
F   As  F , s  x, y, z
t
s
このまま移流方程式を解くことは困難であるため、
[𝐴𝑠 ]の対角化を行い独立した移流方程式を導く。し
たがって次の固有値問題を解く。
 As F  s F , s  x, y, z
(10)
上式より得られる固有値𝜆𝑠 を対角にならべたマト
リクス[𝛬s ]と、固有ベクトルを並べた固有マトリク
ス[𝜑𝑠 ]を用いると次式が成立する。
[s ]1[ As ][s ]  [s ], s  x, y, z
(11)
以上を利用して、{𝐹𝑠 } = [𝜑𝑠 ]{𝑓𝑠 }とすると、独立し
た移流方程式を導くことができる。
Copyright © 2015 Hosei University
2.4.1 M 型 CIP
方向分離解法を用いて、2 次元 CIP 法を解く場合
の手順を示す。
CIP
CIP
f * , g n 
g*
STEP1： f n 
時刻 n において、まず x 方向の移流から始めたと
する。移流値𝑓 𝑛 とその x 方向微分値 g n が CIP 法に
よって中間値となる値が求まる。
CIP
CIP
f n 1 , h* 
hn 1
STEP2： f * 
中間値から y 方向の移流を行い次ステップの移流値
𝑓 𝑛+1 とその y 方向微分値ℎ𝑛+1 が得られる。しかし
STEP2 で必要な𝑓 ∗ は STEP1 で求められていないため、
何らかの方法で求めなければならない。
法政大学情報メディア教育研究センター研究報告 Vol.29
51
この移流方向と垂直方向の微分値を求める手段と
して 2 つの手法があり、その 1 つが M 型 CIP と呼ば
れる。これは 1 次線形補間で中間微分値を求める手
法で、
 y fi *,j   y fi ,nj  ut
 y f
n
i, j
 u t
 y fi n1, j   y fi n1, j
2x
 y fi n1, j   y fi n1, j
2x
,
u0
(18)
,
u0
を用いて計算される。ここで𝜕𝑦 は、𝜕/𝜕𝑦 である。こ
の式は中央差分で解いたことと同じである。
2.4.2 C 型 CIP
もう 1 つの手法が C 型 CIP である。一次線形補間
で値を導いた M 型に対して、C 型では CIP 法を用い
て値を導く。すなわち、更に 2 階空間微分値も記憶
することによって、微分値を CIP 法で求めている。
M 型と C 型を比較すると精度は C 型が優位だが、
記憶する情報が増える点に短所がある。
2.5 T マトリクス
s 方向に移流したのち、r 方向に移流するために s
方向移流方程式を r 方向移流方程式に変換する必要
がある。変換には T マトリクスが用いられる。
T マトリクスは以下の式で表される。
Trs   r  s 
1
(19)
また、T マトリクスには以下のような性質がある。
一般的な関係則を示す。下添え字の数字は空間次元
を表す。
{
[𝑇𝑥 ]1 = [𝑇]1
[𝑇𝑦𝑥 ]2 = [𝑇𝑥𝑦 ]2 = [𝑇]2
(20)
[𝑇𝑦𝑥 ]3 = [𝑇𝑧𝑦 ]3 = [𝑇𝑥𝑧 ]3 = [𝑇]3
[𝑇]1 = [𝐼]
{ [𝑇]2 ∙ [𝑇]2 = [𝐼]
[𝑇]3 ∙ [𝑇]3 ∙ [𝑇]3 = [𝐼]
(21)
移流値の変換は
 fr   Trs  f s 
で計算される。
3．3 次元弾性体波動伝播解析
(22)
で解くこととする。CTP の精度の方が CDP より上で
ある。しかし、小メモリで計算可能なメリットが
CDP にはあるので両者を比較検討する。
3.1.1 CDP 手順
n
CIP1D
*
Tmat
 f y* ,
STEP1： f x  f x 
f y*
f xn CIP1D
f * Tmat
 x 

,
x
x
x
f y*
f xn CDP
f * Tmat

 x 

,
y
y
y
f y*
f xn CDP
f x* Tmat




,
z
z
z
*
CIP1D
**
Tmat
 f z** ,
STEP2： f y  f y 
f y*
y
f
*
y
x
f y*
z

CIP1D
CDP


CDP


f y**
y
f y**
x
f y**
z
Tmat


f z**
,
y
Tmat


f z**
,
x
Tmat


f z**
,
z
**
CIP1D
***
Tmat
 f xn 1 ,
STEP3： f z  f z 
f n 1
f z** CIP1D
f *** Tmat
 z 
 x ,
z
z
z
f xn 1
f z** CDP
f z*** Tmat




,
x
x
x
f n 1
f z** CDP
f *** Tmat

 z 
 x ,
y
y
y
ここで CIP1D(CIP 1 Dimension Solution)は 1 次元
CIP 解法、 CDP は中央差分補間、 Tmat(T Matrix
Transformation)は T マトリクス変換を表す。
3.1.2 CTP 手順
n
CIP1D
*
Tmat
 f y* ,
STEP1： f x  f x 
f y*
f xn
f * Tmat
CIP1D

 x 

,
x
x
x
f y*
f xn
f * Tmat
CIP1D

 x 

,
y
y
y
f y*
f xn
f * Tmat
CIP1D

 x 

,
xy
xy
xy
f y*
f xn
f * Tmat
CIP1D

 x 

,
z
z
z
3.1 計算手順
f y*
f xn
f x* Tmat
CIP1D
3 次元弾性体の波動伝播解析には先に紹介した、



,
zx
zx
zx
M 型 CIP を中央差分で解いた CDP(Central Difference
f y*
f xn
f * Tmat
CIP1D
Profile)と C 型 CIP(CTP：C Type CIP)の 2 種類の手法

 x 

,
yz
yz
yz
Copyright © 2015 Hosei University
法政大学情報メディア教育研究センター研究報告 Vol.29
52
このモデルの中央節点に速度を z 方向へ sin 波形
で入力した波動伝播解析を行う。計算手法は先に述
べた CDP と CTP の 2 種類である。
f y*
f xn
f x* Tmat
CIP1D




,
xyz
xyz
xyz
CIP1D
Tmat
 f y** 
 f z** ,
STEP2： f y* 
f y*
y
f
*
y
x
f y*
xy
f
*
y
z
f y*
CIP1D


CIP1D


CIP1D


CIP1D


f y**
y
f y*
x
f y**
xy
f y*
z
f y**
f z**
,
y
Tmat


Tmat


f z**
,
x
f z**
,
xy
Tmat


Tmat


LZ
f z**
,
z
LX
f **


 z ,
yz
yz
yz
f y*
zx
f y*
xyz
CIP1D
CIP1D


CIP1D


f y*
Fig.1
Tmat
Tmat


zx
f y**
xyz
f z**
,
zx
Tmat


f z***
z
f ***
CIP1D

 z
x
f ***
CIP1D

 z
zx

f ***
CIP1D

 z
y
CIP1D


f z**
,
xyz
𝐿𝑋 = 𝐿𝑌 = 𝐿𝑍 = 20.0𝑚
𝐷𝑋 = 𝐷𝑌 = 𝐷𝑍 = 1.0𝑚
  0.3
  1.5[103 kg m3 ]
f xn 1
,
z
f n 1
Tmat

 x ,
x
f n 1
Tmat

 x ,
zx

f n 1
Tmat

 x ,
y
Tmat


E  56[MN / m2 ]
G  21.54[ MN / m2 ]
Cs  119.83m / sec
Cp  224.18m / sec
t  1.0m s
f n 1
f z** CIP1D
f *** Tmat
 z 
 x ,
yz
yz
yz
3.3 解析結果
解析結果は初期入力点を含む xy 平面の面外速度
描画結果を示す(Fig.2)。
f n 1
f z** CIP1D
f *** Tmat
 z 
 x ,
xy
xy
xy
f n 1
f z**
f *** Tmat
CIP1D

 z 
 x ,
xyz
xyz
xyz
CDP
3.2 3 次元弾性体無限地盤モデル
Fig.1 に 3 次元弾性体無限地盤波動伝播解析に用い
るモデルを示す。Table 1 は解析に用いた物性値等諸
量である。
ここで、LX、LY、LZ は各方向のモデル長さ、DX、
DY、DZ は各方向の節点距離、 はポアソン比、E は
ヤング率、G はせん断弾性係数、Cs は S 波速度、Cp
は P 波速度、∆𝑡は 1 ステップ秒を表す。
Copyright © 2015 Hosei University
図.1 3 次元解析モデル
Three Dimensional Analytical Model
表 1 三次元モデル物性値等諸量
Table 1 Parameters of Analytical Model
**
CIP1D
***
Tmat
 f xn 1 ,
STEP3： f z  f z 
f z**
z
f z**
x
f z**
zx
f z**
y
LY
CTP
10step
法政大学情報メディア教育研究センター研究報告 Vol.29
53
した伝播を見せているが、CDP はざわめきが大きく、
最後には数値拡散を起こしている。 CTP は描画を
見る限り、対称的に伝播していることがわかる。
3.5 3 次元弾性体半無限地盤モデル
CTP による解析の精度を検討するために、半無限
体の表面波伝播理論解である Lamb の解と比較する。
3 次元弾性体半無限地盤モデルは無限モデルに自由
端を表現する境界条件を入れることによって実現さ
れる。
解析モデルは、無限地盤解析時と同様な諸量を持
つ、21 × 21 × 22のモデルで z 正方向境界面に自由表
面処理を行うための節点面を設ける。地表面の中央
節点へ応力𝜎𝑧 を負の方向へインパルスで与えた波動
伝播解析を行う。
30step
60step
3.6 地表面境界条件
CIP は移流方向に対して独立して計算できること
から、逸散波の移流値を入射移流値として代入する
ことで、逸散波が地表面で反射しそれが入射波とな
って入射する現象を表現することができる。
面𝑍 − 1を地表面と考えた場合、その下の面𝑍 − 2
の逸散移流量𝐹𝑂 を面 z の入射移流量𝐹𝐼 として代入す
る。具体的には、x、y 方向の CIP 計算後
90step
120step
 FO  Z  2  FI  Z  , FO  Z  2  FI  Z 
(23)
とした後、z 方向の移流計算を行う。
3.7 Lamb の解
この解[6]は地表面へインパルス力を加え、その点
から距離 r 離れた地点で観測される z 方向変位を表
す。
150step
3.8 解析結果
Fig.4 は、地表面を表す図(Fig.3)の節点番号 1
～4 における、z 方向時刻歴変位応答と、その節点
に対応する Lamb の解を同じグラフに示したもの
である。Fig.5 は、x、y 方向時刻歴変位応答をグラ
フにしたものである。
180step
210step
2
3
1
4
図.2 3 次元無限モデルにおけるｚ速度結果
Fig.2 z Velocity Result of 3D Infinite Space Model
3.4 考察
両描画結果を比較して分かるように、CTP は安定
Copyright © 2015 Hosei University
法政大学情報メディア教育研究センター研究報告 Vol.29
54
Fig.3
図.3 半無限地表面上の観測点
Observation Point on Surface of Half Space
z displacement × 1011 𝑚
各グラフにある 3 本の縦線は左から順に理論解に
おける P 波、S 波、Rayleigh 波の到達を意味する。
無限地盤モデルでは目立たなかったが、x 方向と y
方向とで大きく対称性がくずれている。 C 型でも方
向分離解法を行っている以上は、完全な対称性を保
持することは難しい。
本解析を通して、CTP による 3 次元無限あるいは
半無限地盤解析のプロファイル精度の基本的な特性
を確認することができた。
4．結論
□M 型 CIP では短時間解析については満足なプロフ
ァイルを得る。
□C 型 CIP は無限、半無限に関わらず安定した解析
値を得ることができる。
□方向分離解法の特性上一部の対称性が崩れて、総
合的には十分な対称性を保持する。
time(0.001sec)
xy displacement × 1011 𝑚
図.4 インパルス応力を受ける 3 次元半無限モデル
の 1-4 点ｚ方向変位解析結果と Lamb 解析
Fig.4 z Displacement Result of 3D Half Space Model
Subjected to Impulse at 1-4 and After Lamb
time(0.001sec)
図.5 インパルス応力を受ける 3 次元半無限モデル
の 1-4 点 xy 方向変位解析結果
Fig.5 xy Displacement Result of 3D Half Space model
Subjected to Impulse at 1-4
参考文献
[1]日本建築学会，
“入門・建物と地盤との動的相互
作用”，日本建築学会，1996
[2]伊野慎二，吉田長行，
“波動透過境界の最適化に
関する研究”，法政大学情報メディア情報教育
センター研究報告集 Vol.21，pp.101-108，2008
[3]古谷忍，吉田長行，
“最適化手法による波動透過
境界処理に関する研究”，法政大学情報メディア
情報教育センター研究報告集 Vol.22，pp.55-61，
2009
[4]矢部，尾形，滝沢，
“CIP 法－原子から宇宙までを
解くマルチスケール解法－”，森北出版，2003
[5]矢部，尾形，滝沢，
“CIP 法と JAVA による CG
シミュレーション”，森北出版，2007
[6]Karl F. Graff，
“Wave Motion in Elastic Solids”，Dover，
1991
3.9 考察
z 方向変位結果では Rayleigh 波の到達がはっきり
と表れ、Lamb の解と同時期に到達していることか
ら、本解析法は妥当だと言える。しかし、波通過後
に残留変位が残ってしまっている。sin 波速度入力な
ど、最終的に変位を 0 に抑え込む解析では残留変位
が残らないことから、インパルス入力特有の問題と
考えられる。
x、y 方向変位結果では P 波と S 波の到達で、ピー
クが現れるがこれも、Lamb と同時期に到達を確認
できる。やはりこの結果でも残留変位がみられる。
Copyright © 2015 Hosei University
法政大学情報メディア教育研究センター研究報告 Vol.29

Download Report