著者 佐々木 豊, 柳川 智隆, 松下 周平, 吉田 長行
出版者 法政大学情報メディア教育研究センター
雑誌名 法政大学情報メディア教育研究センター研究報告
巻 29
ページ 49‑51
発行年 2015‑04‑01
URL http://doi.org/10.15002/00012061
法政大学情報メディア教育研究センター研究報告 Vol.29 2015年
原稿受付 2015年3月9日 発行 2015年 4月1 日
CIP 法を用いた 2 次元 FEM 領域の開境界処理法
Open Boundary Processing for the Two Dimensional FEM Domain by CIP Method
佐々木 豊１） 柳川 智隆１） 松下 周平１） 吉田 長行２）
Yutaka Sasaki, Tomotaka Yanagawa, Syuhei Matsushita, Nagayuki Yoshida
１）法政大学大学院デザイン工学研究科建築学専攻
２）法政大学デザイン工学部建築学科
When we analyze wave propagation, we have to model the ground as if it is extended infinitely in the limited analytical region. Now we try to make an infinite ground model by FEM and CIP combined method. CIP method is usually used to analyze sound, noise and electromagnetic field.
In this research, we propose the wave transmitting boundary by CIP method attached to the 2-dimensional FEM domain.
Keywords : CIP Method, FEM, Wave Transmitting Boundary
1. はじめに
地盤の上に建つ構造物は、地盤の影響を受けて振動 している。構造物の周期は地盤の影響を受けて増大 し、構造物がもつ振動エネルギーの一部は地盤に放 出され、構造物の振動エネルギーは減少する。ある いは、地震動がもつ高い振動数の成分は、構造物を 振動させるのに有効な入力とならないこと、隣接す る構造物は地盤を介してお互いに干渉し合うという こともある。これらは地盤と構造物とが影響しあう ことに起因しており、そのような現象を地盤と構造 物 と の 動 的 相 互 作 用 (Dynamic Soil-Structure
Interaction)と呼ばれている[1]。この動的相互作用を
考慮した地震応答を計算する方法の 1 つとして、地 盤と建物を一体の系として計算する方法は直接法ま たは一体解法と呼ばれている。この解析は複雑な地 層構造や地盤の非線形性が容易に解析に取り込める という長所があり、しばしば有限要素法が利用され る。しかし、有限要素法は本来、有限領域を対象と する数値解析手法である。そのため、無限あるいは 半無限弾性体の波動伝播問題に適用する場合には、
Fig.1のように、解析領域の内部と外部に境界を設け、
外部から内部へ伝わる入射波と内部から外部へ逸散 する反射波の双方同時の処理が要求される。有限な 狭領域で波動の完全透過ができる境界処理法は確立 されておらず、実現すれば地盤と建物の解析効率の 向 上 に と っ て 有 益 で あ る[2][3]。 本 研 究 で は 、 CIP(Constrained Interpolation Profile)法を用いた新し い境界処理法の確立を目指している[4][5]。CIP法は移 流方程式を解く有効な解法であり、主に有限要素法 とは異なる分野で用いられている。そこで、如何に して有限要素法とCIP法を組み合わせるかが課題と なる。
Incident wave
Reflected wave Boundary processing
Copyright © 2015 Hosei University 法政大学情報メディア教育研究センター研究報告 Vol.29 図.1 解析領域
Fig.1 Analytical region.
2. 研究過程
移流方程式の定式化 ↓
境界処理法の決定 ↓
入反射場CIP法解析 ↓ 時刻歴応答解析
↓
入反射場解析の検証 ↓
END
図.2 解析プロセス Fig.2 Analytical process.
3 次元への拡張を見据えた地盤の解析を目的とし、
その足がかりを得るために、1 次元棒材モデルでの 入反射場解析[6]を基に、2次元地盤モデルを用いた手 法の提案と検討を行う。
3. 解析
本研究では、FEMとの結合解析を2つの手法を用 いて行う。
3.1 解析モデル
２次元解析におけるモデル図と材料特性を以下に 示す。
●：FEM mesh
○：CIP grid 図.3 解析モデル
Fig.3 Analytical model.
表 1 材料特性 Table 1 Material property.
Length Lx =18 m ,
y 17 L = m Secondary wave velocity cs =150m s/
Primary wave velocity cp =281m s/ Density ρ=1500kg m3
/
Elastic shear modulus G= ×ρ cs2
Poisson ratio υ=0 3. 3.2 移流値算定
解析モデルへ外部から入射する波形データの作成 が必要である。FEM-CIP結合解放では入射値を境界 力として与えるために、速度や応力を移流値に変換 する必要がある。以下では2次元問題を考慮した速 度入力の場合での計算式を示す。
SH波の場合：
0 0
( / )
( )sin ( ) cos
z s
u f t s c
s x x θ y y θ
= −
= − + − (1)
,
,
/ / ,
/ / / sin /
/ / / cos /
z
z x s
z y s
u f t f t f
u f x f t t x f c
u f y f t t y f c
θ θ
∗
∗ ∗
∗ ∗
= ∂ ∂ = ∂ ∂ =
= ∂ ∂ = ∂ ∂ ⋅∂ ∂ = −
= ∂ ∂ = ∂ ∂ ⋅∂ ∂ = −
(2)
, ,
sin / cos /
zx zx z x s
yz yz z y s
G Gu Gf c
G Gu Gf c
τ γ θ
τ γ θ
= = = −
= = = −
(3)
1
2
2
1 ( 1 sin )
2 2
1 ( 1 sin )
2 2
1 ( 1 cos )
2 2
zx z x
s s
zx z x
s s
z yz y
s s
u f
f c G c
u f
f c G c
u f
f c G c
τ θ
τ θ
τ θ
= + + = + −
= − + = − −
= − + = − −
(4)
P・SV波の場合：
0 0
( / ), ( )sin ( ) cos
s s
u = f t−s c s= −x x θ+ −y y θ (5)
cos cos , sin sin
x s y s
u = −u θ= f θ u =u θ= f θ (6)
2
, ,
2
, ,
sin cos / , cos / ,
sin / , sin cos /
x x s x y s
y x s y y s
u f c u f c
u f c u f c
θ θ θ
θ θ θ
= =
= − = −
(7)
y
x
Copyright © 2015 Hosei University 法政大学情報メディア教育研究センター研究報告 Vol.29
1 , 2 , 1
2 , 1 , 1
2 2
, ,
2 sin cos ,
2 sin cos ,
( ) (cos sin )
s
x x x y y
p p
s
y x x y y
p p
xy x y y x
s
f c
D u D u D
c c
f c
D u D u D
c c
G u u G f
c
σ θ θ
σ θ θ
τ θ θ
= + = + ⋅
= + = − ⋅
= + = −
(8)
1
1
2
1
2 2
3
2 4
1 1 2 sin cos
2 2
1 1 2 sin cos
2 2
1 ( sin cos sin )
2 2
1 ( sin cos
2 2
x x s
x
p p p
x x s
x
p p p
y xy
x
s s
y xy
x
s s
u f c
f c D c c
u f c
f c D c c
u f
f c G c
u f
f c G c
σ θ θ
σ θ θ
τ θ θ θ
τ θ θ
= + + = − +
= − + = + +
= + + = + + −
= − + = − + −
2
2
1
2 2
4
sin )
1 1 2 cos sin
2 2
1 ( cos cos sin )
2 2
y y s
y
p p p
x xy y
s s
u f c
f c D c c
u f
f c G c
θ
σ θ θ
τ θ θ θ
= − + = − −
= − + = + + −
(9)
ここで、1−D2/D1=(D1−D2) /D1 =2 /G D1=2cs2/cp2
3.3 FEM-CIP 結合解析
2次元問題を取り扱う場合にも、1次元問題と同様 な手法を使うことが可能である。半無限モデルの場 合は地表面を覆わない形を想定する。つまり Fig.3 のように、コの字型のCIPグリッドで各々開境界処 理をする。
しかしCIP法の特性上、方型メッシュのみ解析が 可能であり、一度にコの字型のグリッドの解析を行 うことができない。そのため、グリッドをいくつか のCIPメッシュに分割する必要がある。
3.3.1 Open Boundary Side scheme(OBS)
境界面に与える外力（以下，境界力）を Fig.4 の ように右面、左面、下面のそれぞれ一組のCIPグリ ッドで計算する手法である。入射用グリッドには入 射波形を記憶し、逸散波用グリッドには FEM デー タを受渡する。FEM領域外のグリッドにはFEMデ ータを外挿している。計算効率は良いものの、メモ リ量はFEM分割に比例して大きなものとなる。
図.4 解析イメージ(OBS) Fig.4 Analytical image(OBS).
3.3.2 Open Boundary Cell scheme(OBC)
Δt秒間での移流で、境界点に影響を与える情報量 は限られる。そこで、境界点ごとに Fig.5 のような 4×3のCIPグリッドを用いて、境界力を計算する手 法を提案する。境界点の数だけ計算が行われるため、
手数は多くなるもののFEM使用メモリを除き、4×3 メッシュだけのメモリで計算が可能であり、大きな モデルを対象とした場合メモリの節約に期待できる。
基本的な計算過程はOBSと同じである。
図.5 解析イメージ(OBC) Fig.5 Analytical image(OBC).
3.4 計算手順
SH 波問題、P・SV 波問題における計算手順をそ れぞれ以下に示す。OBS とOBCとでは計算手数が 違うだけであり、OBSでは境界辺の数だけ、OBCで は、境界点の数だけ下記の境界力算定を行うことと なる。
( , ) fO i j ( , )
f i jI
(i) (i)
(j) (j)
FEM region CIP mesh
Incoming
Evaluation point for open boundary force Outerpolation point
Input point /
Copying point from FEM Outgoing
Advection
direction Sliding scheme
3 2 (j) 1
3 (j) 2 1
1 2 3 4 (i) (i) 1 2 3 4
FEM region CIP mesh
Incoming
Evaluation point for open boundary force Outerpolation point
Input point /
Copying point from FEM Outgoing
Copyright © 2015 Hosei University 法政大学情報メディア教育研究センター研究報告 Vol.29 表 2 SH 波計算手順(OBS)
Table.2 SH wave procedure(OBS)
Input wave f i j iI( , ), =1, 2,3, 4,j=1, 2,,NS Copying
point from FEM
( , ), 2,3, 4, 1, 2, ,
O S
f i j i= j= N
Here, j=1,Ns obtained by outerpolation
( , ) 2 ( , 1) ( , 2),
( ,1) 2 ( , 2) ( ,3)
O s O s O s
O O O
f i N f i N f i N
f i f i f i
= − − −
= −
Dotted point obtained by outerpolation (1, ) (4, ) 3( (3, ) (2, ))
O O O O
f j = f j − f j − f j
Gradient
culculation gO:i=2,3, j=1,2, , NS(CDP) CIP
Advection
Dissipation advection 2, 1,2, , S
i= j= N
Boundary
force τ τzx, yz
表 3 SH 波計算手順(OBC) Table.3 SH wave procedure(OBC) Input wave f i j iI( , ), =1, 2,3, 4,j=1, 2,3
Copying point from
FEM
( , ), 2,3, 4, 1, 2,3 fO i j i= j=
Thinking boundary force of the corner, 1
j= orj=3 obtained by outerpolation
( ,3) 2 ( , 2) ( ,1), ( ,1) 2 ( , 2) ( ,3)
O O O
O O O
f i f i f i
f i f i f i
= −
= −
Dotted point obtained by outerpolation (1, ) (4, ) 3( (3, ) (2, ))
O O O O
f j = f j − f j − f j
Gradient
culculation gO:i=2,3, j=1,2,3(CDP) CIP
Advection Dissipation advection i=2, j=2 Boundary
force τ τzx, yz
表 4 P・SV 波計算手順(OBS) Table.4 P・SV wave procedure(OBS) Input wave f i j iI( , ), =1, 2,3, 4,j=1, 2,,NS
Copying point from
FEM
( , ), 2,3, 4, 1, 2, ,
O S
f i j i= j= N
Here, j=1,Ns obtained by outerpolation
( , ) 2 ( , 1) ( , 2),
( ,1) 2 ( , 2) ( ,3)
O s O s O s
O O O
f i N f i N f i N
f i f i f i
= − − −
= −
Dotted point obtained by outerpolation (1, ) (4, ) 3( (3, ) (2, ))
O O O O
f j = f j − f j − f j
Gradient
culculation gO:i=2,3, j=1,2, , NS(CDP)
: 1,2,3, 2, , 1
O S
h i= j= N − (CDP)
: 1,2,3, 1,
O S
h i= j= N (1DLiner)
CIP Advection
Dissipation advection 2,3, 1,2, , S
i= j= N
T matrix transform Vertical advection
2, 2, , S 1
i= j= N − Boundary
force σ σx, y
表 5 P・SV 波計算手順(OBC) Table.5 P・SV wave procedure(OBC) Input wave f i j iI( , ), =1, 2,3, 4,j=1, 2,3
Copying point from
FEM
( , ), 2,3, 4, 1, 2,3 fO i j i= j=
Thinking boundary force of the corner, 1
j= orj=3 obtained by outerpolation
( ,3) 2 ( , 2) ( ,1), ( ,1) 2 ( , 2) ( ,3)
O O O
O O O
f i f i f i
f i f i f i
= −
= −
Dotted point obtained by outerpolation (1, ) (4, ) 3( (3, ) (2, ))
O O O O
f j = f j − f j − f j
Gradient
culculation gO:i=2,3, j=1,2,3(CDP) : 1,2,3, 2
hO i= j= (CDP) : 1,2,3, 1,3
hO i= j= (1DLiner) CIP
Advection
Dissipation advection i=2,3, j=1,2,3 T matrix transform
Vertical advectioni=2, j=2 Boundary
force σ σx, y
Copyright © 2015 Hosei University 法政大学情報メディア教育研究センター研究報告 Vol.29 4. 解析結果
以下ではsin波1波を入射角30°で下方より斜め入 射した場合について、SH波問題(Fig.6)、P・SV波問
題(Fig.7)それぞれにおける全質点の変位挙動を示す。
OBSと OBCで同じ解析値が得られたため、ここで はOBCで得られた解析結果を掲載する。SH波問題 では面外変位描画を、P・SV波問題では面内変位描 画結果を表す。
SH波問題：
20step 40step
60step 80step
100step 120step
140step 160step 図.6 半無限モデルにおける SH 波解析結果 Fig.6 The half space model subjected to SH wave
with
30-degree of incident angle: Analysis by OBC.
20step 40step
60step 80step
100step 120step
140step 160step 図.7 半無限モデルにおける P・SV 波解析結果 Fig.7 The half space model subjected to P・SV wave
with
30-degree of incident angle: Analysis by OBC.
20step
Copyright © 2015 Hosei University 法政大学情報メディア教育研究センター研究報告 Vol.29 以下ではFig.8に示す質点番号86、176、186、333
のOBSとOBCのP・SV波問題での時刻歴変位結果 を示す。
図.8 半無限モデルの解析点 Fig.8 Observation point in the half space model.
図.9 質点番号 86 の時刻歴変位挙動 Fig.9 Displacement of node 86 in (Incident angle 30°).
図.10 質点番号 176 の時刻歴変位挙動 Fig.10 Displacement of node 176 in (Incident angle 30°).
図.11 質点番号 186 の時刻歴変位挙動 Fig.11 Displacement of node 186 in (Incident angle 30°).
図.12 質点番号 333 の時刻歴変位挙動 Fig.12 Displacement of node 333 in (Incident angle 30°).
5．結論・考察
本研究の結論として以下のものが挙げられる。
OBCとOBSとでSH波、P･SV波伝播が共にほ ぼ同一の結果が得られた。
多次元においてもCIP法による境界処理の有効 性が確認できた．
OBCと OBSのそれぞれの解析結果を見るとほぼ 同一の結果であると判断できる。これにより OBC の解析の有効性が確認できる。CIP法ではΔt後の微 小時間での波動データしか扱いわないため、CIP 法 メッシュを制限してもCIP移流されるこの手法は3 次元問題にでも扱えると考える。また、斜め入射や 鉛直入射等の入射角度を変えても同様な結果を得ら れた。
86
176 186
333
( 0.001425sec) time×
( 0.001425sec) time×
[m] [m] [m]
( 0.001425sec) time×
( 0.001425sec) time×
[m]
Copyright © 2015 Hosei University 法政大学情報メディア教育研究センター研究報告 Vol.29 参考文献
[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] 神崎壮一郎, 吉田長行、" FEM-CIPによる1次元 入反射場解析"、法政大学情報メディア情報教育セ ンター研究報告集Vol.27、pp.68-72、2013
