214INSS JOURNAL,Vol.17(2010),pp214-221

弾性波伝播有限要素解析における要素分割選定方法の考案
A Proposal for a Determination Method of Element Division on an Analytical Modelfor Finite Element Elastic Waves Propagation Analysis
石田仁志（Hitoshi Ishida)＊ 1，

飯井俊行（Toshiyuki Meshii)＊ 2

要約 弾性波の伝播現象の有限要素解析を行う解析モデルに対して最適な要素分割を選定する手法を考案した．この手法は，次のような特徴を有する．（1）要素分割により変化する解析モデル全体のひずみエネルギーに着目し，その収束性に基づいて要素分割を選定する．
（2）時間とともに変化する加振入力変位のうち，離散化した時間ステップにおいて最大に変化する変位とこれにより生じる慣性力を境界条件として入力する静的解析（時間 1 ステップのみの動的解析に相当）を行う．時間 1 ステップの間の変位と慣性力が衝撃的に変化するという特徴に着目したことから IM（ImpactMeshing）法と名づけた．
3 次元弾性体における超音波伝播問題を例にとって，IM 法によって選定した要素分割された解析モデルを使っても，有限要素解析による解が収束することを示した．さらに，時間 100 ステップの解析の場合，従来の有限要素解析の解の収束判定によるという従来の要素分割選定方法に比べて，IM 法では要素分割選定のための計算時間が約 1/6 に削減されることも示した．
キーワード

弾性波，音波，有限要素解析，要素寸法，ひずみエネルギー

Abstract
This study proposes an element size selection method named the ³Impact-Meshing (IM) method· for a finite element waves propagation analysis model, which is characterizedby (1) determination of element division of the model with strain energy in the whole model, (2)static analysis (dynamic analysis in a single time step) with boundary conditions which gives amaximum change of displacement in the time increment and inertial (impact) force caused by thedisplacement change.
In this paper, an example of application of the IM method to 3D ultrasonic wave propagation
problem in an elastic solid is described. These examples showed an analysis result with a model
determined by the IM method was convergence and calculation time for determination of element
subdivision was reduced to about 1/6 by the IM Method which did not need determination ofelement subdivision by a dynamic transient analysis with 100 time steps.
Keywords

1.

elastic waves, acoustic waves, finite element analysis, element size, strain energy

はじめに問題や求める解の扱いに応じて要素寸法と時間増分が個別に選定されたり，または同様の先例が適用さ鋼材内の超音波や地盤内の地震波などの伝播現象れたりしていることが多い．の解明のために，弾性波伝播解析が行われ，その手すなわち，時間増分については，劉(1)による入射法の一つとして，有限要素解析が用いられる．有限波形を滑らかに表現するために入力波の周期の分割

要素解析においては，適切な解を得ることだけでな数による選定方法，および大槻(2)，Jaleel(3)，伊藤(4)，く，計算効率の点からも，要素寸法と時間増分の選

林 (5)，Ju (6)，Zhang (7) らによる弾性波の伝播速度に定が重要である．

基づく（要素寸法を伝播する時間による）選定方法

しかし，解くべき弾性波伝播に係る方程式が同じがよく使用されている．

であるにもかかわらず，要素寸法と時間増分の選定

要素寸法については，Harumi(8)，三木(9)，劉(1)，

方法に必ずしも統一されたものはなく，対象とする

Guan (10)，大槻 (2) らは，入力波の波形を滑らかに表＊1＊2(株)原子力安全システム研究所
福井大学大学院工学研究科

技術システム研究所

215

現できるような波長の分割数によって選定している．

題に対し解の収束を判定するパラメータとして，入一方，Galdos(11)，Moser(12)，Mirkhani(13)，吉川(14)，

力された変位と慣性力から求められ，解析対象領域

Nair(15)らは，要素細分化に伴って変化する解析結果全体を代表する変化量，すなわち解析モデル全体にの収束によって選定している．前者の選定基準は定蓄えられるひずみエネルギーを考えることとする．

性的であり，後者は適切な解を得ることにおいて有ここで，与えられた Δt に対して，変位の絶対値効であると考えられるが，いつの時刻でどのパラメが最大に変化する時刻 t=nΔt~n+1Δt の間の変位ータの収束を見ればよいのかが，必ずしも明確に示

ベクトル変化 Δu  を特定し，この Δu  がステされてはいない．

ップ状に与えられる結果発生する慣性力 F  を境

先例のある解析対象を扱うのであれば，その解析界条件として与えることにする．絶対値が最大に変例を参考として時間増分や要素寸法を選定すれば実化する時刻の変位変化によって蓄えられるひずみエ用上問題はないと考えられるが，これまでに実績のネルギーが収束すれば，最大に変化する変位絶対値ない解析対象に対し弾性波伝播有限要素解析を行う

以下の他の各時刻においては，ひずみエネルギーは場合にも適用できる，汎用的な時間増分および要素当然収束すると考えられる．Δu  という変位入力寸法選定方法があれば，実用上有益であると考えらとこれによって衝撃的に発生する慣性力 F  とをれる．

解析モデルに加えることによって，解析モデル全体

そこで，本研究では，時間増分がすでに選定されに蓄えられるひずみエネルギーを評価し，その収束ているとの前提の下で，対象とする問題によらない，性をもとに要素分割を選定する方法であることから，統一的，かつ判断根拠が明確であるという意味におこの方法を IM（Impact Meshing）法と呼ぶことといてより合理的な要素分割選定法を提案する．以下，する．

提案する方法の着眼点と手順を示し，弾性波伝播解析例によってその有用性を示す．

2. 要素分割選定方法

2.2

要素分割選定手順

図 2 に，要素分割選定の手順を示す．具体的な手順は以下のとおりである．なお，IM 法では，時間2.1 着眼点増分は既に与えられていることを前提とする．
Step1

変位の絶対値が最大に変化する時刻正弦波状に変化する変位を入力し，弾性波としてt=nΔt~n+1Δt の間の変位ベクトルの変伝播させるという動的解析問題について考える．図化 Δu  を選定し，これを解析モデルに1（1）のように時々刻々連続的に変位が変化する正与える変位とする．

弦波は，時刻歴変化を扱う有限要素解析においては，

Step2

Δu  とこれにより生じる慣性力 F 

変位の変化を図 1（2）のように時間増分 Δt によりが，解析モデルに与えられる力となる．こ離散化し，ステップ状に変化する変位として入力さのとき，F  は Step3 において動的解析れる．このステップ状に変化する変位 Δu によって方法によって慣性力として扱うことが可能

速度が，さらにその速度とそれ以前の時刻の速度とである．の変化によって加速度が生じ，その結果図 1（3）のStep3

解析領域をモデル化し，要素分割する．

ように慣性力 F が生じる．したがって，時間ステッStep4

境界条件は，実際に解くべき時刻歴解析

プごとに離散化された動的解析問題は，その時間 1問題における加振変位 ut を Δu  に置

ステップにおける変位変化および慣性力が境界条件き換える他は，実際の問題と同一とする．

として入力される静的解析問題と同等であり，図 1

静止状態を初期条件とし，動的解析を時刻（4）のように，動的解析問題は，1 コマずつの静的Δt の時間 1 ステップだけ行い，解析の結解析問題を連ねたものに置き換えられる．

果得られる応力 σ とひずみ ε の分布から，蓄えられたひずみエネルギー U  を式

このことから，弾性波伝播という動的解析問題に（1）により算出する．

おける要素分割の選定は，各時間ステップにおける静的解析問題に対する適切な要素分割を選ぶ問題に帰着できると考えられる．そして，この静的解析問

U =1σ  ε+σ  ε+σ  ε+τ  γ  +τ  γ  +τ  γ  dV (1)
2 

216

1.01.0∆u i |t =(n +1) ∆t㞫ᩋ໩∆u i |t =n ∆t

Fiu i /A 0u i /A 00.00.0

F i |t =n ∆t
F i |t =(n +1)∆t∆t-1.0-1.0

n n+1t/∆t(2) 離散化された入力変位t =n∆t
㻧u|t =n∆tt/∆t(3) 慣性力

1.0F|t =n∆tui/A0ንనථງ ut/∆t∆u i |t =n ∆tFi(1) 弾性波となる入力波∆tn n +10.0

F i |t =n ∆tt =(n +1)∆tF|t =(n +1) ∆t∆u|t =(n +1)∆tnui/A01.0∆u i |t =(n +1) ∆t0.0-1.0nt/ ∆tt/∆tFi-1.0

F i |t =(n +1)∆tn +1t/∆t

n +1t/∆t(4) 時間ステップ毎の変位変化と慣性力による静的解析
図1

弾性波伝播問題

え，8 節点六面体要素を用いて，x ，y ，z 方向の寸法
Step5

要素細分化前の U  と今回の要素細分をそれぞれ n ，n ，n に等分割した．

化後の U  を比較し，工学的に U 

時間増分 Δt は，図 5 に示すように，正弦波を離が収束していないと判断すれば Step3 以散化したとき，正弦波の振幅値が誤差 5％以内で再降を繰り返し，収束していると判断すれば現可能な分割数となるように，加振入力変位の周期要素細分化を打ち切る．

T の 1/10 とした．
この問題に対し，IM 法を適用し，適切な要素分3.

適用例割を選定する．図 2 の Step4 は，具体的には初期条件を弾性波の一つである鋼材中の超音波の伝播問題への IM 法の適用を検討した．

u0=0(2)図 3 に，検討した平鋼板表面上の探触子から超音波を入射する 3 次元超音波伝播問題を示す．
寸法 L=2.0×1010

m ，W =2.0×10m ，H =2.0×次の境界条件

m の静止状態にある鋼板において，その表面上に 配 置 す る 探 触 子 寸 法 L=5.0×10 m ，W = 5.0×10とし，x ，y ，z 方向の節点変位 u ，v ，w により，mmの範囲が一様に時刻 t>0 に変位 ut を与u<L/4,<W /8,0;Δt=Δu
uL,W /2,0;Δt=0,uL,W /2,H ;Δt=0えることにより加振されるとした．加振入力条件を

vx,0,z;Δt=0

表 1 に材料定数を表 2 に示す．

w<L/4,<W /8,0;Δt=Δu

図 3 の問題を図 4 のようにモデル化し，有限要素解析を行う．対称性を考慮し，y ←0 の 1/2 領域を考w0,W /2,H ;Δt=0,wL,W /2,H ;Δt=0(3)217Step 2F i impact∆ui impact

0.0

Fiui/A0ථງንన㻃∆uimpact 䛴㐽ᏽ

Step 1

1.0

Step 1㛜ጙ∆t-1.0

ៈᛮງ Fimpact

Step 2nn+10t/∆t∆u impact

Step 3
Step 3ゆᯊ䝦䝋䝯䛴こ⣪ฦ๪

Step 4ゆᯊ( U impact 䛴⟤ฝ)z

U impact 䛴཭ᮨYes཭ᮨ

U impact ＝こ⣪ᩐ

図2H

㗨ᯀx

㉰㡚ἴ

(0,0,0)(0,0,H )W(L ,0,0)䝿䝿 䝿xzy

Ly1
σ ij ε ij d V
2 V∫

IM 法による要素分割選定の手順

᥀よᏄWt

ゆᯊ䝦䝋䝯

Step 4 䛸䛠䜅䜬䝑䝯䜲䞀⤂஡Lt

FimpactxU impact

Not/∆ty

Step 5
Step 5

1(0,W /2,0)z

図33 次元超音波伝播問題(0,W /8,0)䝿䝿(0,0,0) 䝿䝿䝿䝿(L /4,0,0)䝿䝿 䝿䝿䝿 䝿x(L ,0,0)

表 1 加振入力条件
波形
周波数f振幅A

正弦波
1.0 × 106Hz
1.0 × 10−11m

図4解析モデル（3 次元超音波伝播解析）

表 3 に示す要素数の異なる 5 種類の解析モデルにつ
表2

材料定数

ヤング率E
ポアソン比 v
密度ρ206GPa0.37700kg/m3

いて求めた．表 3 に示す 5 種類の解析モデルによるIM 解析結果から，ひずみエネルギー U  の収束状況を図 6 に示す．図中，U   は表 3 のモデルⅤに対する値であり，モデルⅢが 0.9772，モデルⅣが0.9996 であった．モデルⅣからⅤに要素数を増加さを与え，t=0~Δt のみの動的解析を行った．解析は，せたときの，ひずみエネルギーの変化は 0.04％であ汎 用 構 造 解 析 コ ー ド MSC Marc を 用 い て 陰 解 法

り，収束していると判断した．

（Houbolt 法）により行った．

次に，本来解くべき図 4 に示す正弦波を入力する

この t=Δt における解析結果を，以下 IM 解析結

時刻歴解析を解析モデルⅠ〜Ⅴの 5 ケースに対して果と称する．その時のひずみエネルギー U  を，行った．表 3 に示した 5 種類の要素数の異なる解析

218

表31.0Ⅰ8Ⅱ20Ⅲ40Ⅳ80Ⅴ120nn

411222044086012

要素数
要素寸法(mm)322.5

4801.03 20025 600 86 400

0.50.25

u i /A 00.50.0

解析モデル（3 次元超音波伝播解析）

モデル
n

0.050.17-0.5-1.0246810∆t図5時間増分の選定

U / U impact V01.0䕝 䕱䕵 䕹model 䊗 䊘 䊙 䊚 䊛

モデルによる解析結果として，図 7 にひずみエネルギーの時刻歴変化を，図 8 に観測点（ L ，0，0）に0.5おける各種計算値についての時刻歴変化を示す．図

107 および図 8 において，横軸は入力変位の周期 T に1102より規格化した時刻である．また，Δt 毎のひずみエネルギーおよび各種計算値を直線で結んで示してい図6る．
図 7 に示した時間 100 ステップまでのひずみエネによって評価した時刻 Δt だけでなく，それ以降の2.0U x 10

122.5model 䊗 䊘 䊙 䊚 䊛1.5

うち，加振入力の変位方向である x および z 方向に0.5関係するものとして，変位 u ，w ，体積ひずみ ε ，0.00.0

せん断ひずみ γ ，回転 ω を示している．
これらの計算値について，IM 法によりひずみエデルⅢおよびモデルⅣの誤差を評価した．表 4 に，51.0

図 8 には，観測点（ L ，0，0）における計算値のネルギーが収束したと判断したモデルⅤに対するモ

10

ひずみエネルギーの収束（3 次元超音波伝播解析）のひずみエネルギーの差は 2％以内であった．IM 法

さく，モデルⅤにおいて収束していると言える．

43.0

ルギーの時刻歴変化において，モデルⅣとモデルⅤ時刻においてもモデルⅣとモデルⅤの結果の差は小31010
こ⣪ᩐ2.04.0t/T

6.0

8.010.0

図 7 ひずみエネルギーの時系列解析結果（3 次元超音波伝播解析）

観測点（ L ，0，0）において弾性波の伝播に相当す
る変位 u ，w ，回転 ω のピーク値における誤差を示す．IM 法による評価においてひずみエネルギーが4.

考察収束しているモデルⅣとモデルⅤの値はよく一致した結果が得られていると言える．

IM 法の最も優れていると考えられる点は，従来以上から，IM 法による解析モデルの要素分割選のような解析モデルの要素細分化に際し，要素分割定が，実際の時刻歴解析においても適切であることの異なる解析モデルの時系列解析を繰り返し実施し，が示された．

その結果を比較して解の収束を判定するという手順が不要となり，要素分割の異なる解析モデルについて時間 1 ステップのみの解析で済むことである．

4.00.5

2.0

31.0

γzx µε×10u/A0219

0.0-0.50.0-2.0model 䊗 䊘 䊙 䊚 䊛model 䊗 䊘 䊙 䊚 䊛-1.0-4.00.0

2.0

4.0

6.0

8.0

10.0

0.0

2.04.0

t/T

4.0

0.5

2.0

31.0ωy µε×10w/A010.0

8.0

10.0(4) せん断ひずみγ 0.0-0.50.0-2.0model 䊗 䊘 䊙 䊚 䊛model 䊗 䊘 䊙 䊚 䊛-1.0-4.00.0

2.0

4.0

6.0

8.0

10.0

0.0

2.0

4.0t/T

6.0t/T(2) 変位w(5) 回転ω
図84.0
38.0t/T(1) 変位uεv µε×106.0

弾性波に係る計算値の時系列解析結果（3 次元超音波伝播解析）

2.0

表40.0-2.0model 䊗 䊘 䊙 䊚 䊛-4.00.0

2.0

4.0

6.0

8.0

10.0t/T(3) 体積ひずみε

表 5 に，要素分割を選定するために，時系列解析による方法と IM 法とによる計算時間の比較を示す．
3 項で適用例として示した 3 次元超音波伝播解析に

モデル
変位u変位w
回転ω
表5

解析モデルによる計算値の誤差（3 次元超音波伝播解析）
Ⅲ0.503-0.758
0.615-0.775
0.118-0.619

Ⅳ
1.034-1.042
0.993-1.020
0.936-1.026

Ⅴ1.0001.000
1.000

要素分割選定のための CPU 計算時間
（Unit：s）

モデル
従来の方法
IM 法

Ⅰ
1.1
0.1

Ⅱ
9.1
0.3

Ⅲ
61.4
2.3

Ⅳ
503.1
82.5

Ⅴ
1 928.4
368.4

おける，表 3 の 5 つの解析モデルについての，時間100 ステップの時系列解析の計算時間と IM 法すな

また，図 6 で示したひずみエネルギーのモデルⅣ

わち時間 1 ステップの計算時間を示している．従来からモデルⅤへの変化は，0.04％に対して，表 4 にの時系列解析結果による選定方法の約 1/6 の計算時示した計算値のモデルⅣからモデルⅤへの変化は，間で選定が可能であった．

0.02-0.07％であった．ひずみエネルギーの変化からさらに大きな領域を対象とし，超音波の入射位置計算値の誤差の推定が可能なことを示唆するものでから伝播させる目的とする位置までの距離が長くなあるが，今後さらなる解析事例によるデータの蓄積るような伝播問題を解析する場合には，超音波の伝が望ましい．

播時間が長くなることによって時間ステップ数も増大する．このような場合には，IM 法の適用がさらに有効であると考えられる．

220

5.

おわりに³線集束超音波探触子により励起された漏洩表
面波伝播の有限要素法解析，· 日本機械学会論

弾性波伝播有限要素解析において適正な解を与え文集．A 編，65［631］，589-596（1999）．

るための，適切に要素分割された解析モデルを選定(５) 林山，伊藤智啓，川嶋紘一郎，³薄肉円管周方するための方法である IM 法を考案し，鋼材におけ向の波動伝搬の FEM 解析と軸方向溝状欠陥のる超音波伝播問題への適用性を示した．

寸法同定，· 日本機械学会論文集．A 編，67

IM 法の特徴は以下のとおりである．

［653］，36-42（2001）．

(1)要素分割により変化する解析モデル全体のひずみ(６) S. H. Ju, Y. M. Wang, ³Time-dependent

エネルギーに着目し，その収束によって要素分割

absorbing boundary conditions for elastic waveを選定する．

propagation,·(2)時間とともに変化する加振入力変位のうち，離散化した時間ステップにおいて最大に変化する変位とこれにより生じる慣性力を境界条件として入力

International

Journalfor

Numerical Methods in Engineering, 50 [9] ,2159-2174 (2001).
(７) J. Zhang, D. J. Verschuur, ³Elastic wave

する静的解析（時間 1 ステップのみの動的解析にpropagation相当）を行う．

media

このことから，IM 法の優れている点が，以下のmethod.," Geophysics, 67 [2] , 625-638

とおり挙げられる．

inusingheterogeneousthelumped

anisotropic

finite-element(2002).

(1)解析モデル全体のひずみエネルギーを評価するた(８) K. Harumi, ³Computer simulation of ultraso-

め，解析モデルの大きさ，加振入力点および弾性

nics in a solid,· NDT International, 19 [5],波観測点の位置，弾性波の伝播時間にかかわらず適用可能である．

315-332 (1986).
(９) 三木千寿，高橋政則，舘石和雄，劉銘崇，³超(2)異なる要素分割の解析モデルについての時系列解

音波探傷の高精度定量化に関する研究 その析による解の収束評価が不要であり，これに伴っ1．超音波の数値シミュレーションとモード分て要素分割選定のための解析時間の短縮が可能で離画像化について，· 東京工業大学土木工学科ある．

研究報告，75-100（1992）．

ひずみエネルギーの収束判定時の変化量によって，(10) J. Guan, Z. Shen, X. Ni, J. Lu, J. Wang, B.

解析結果である弾性波に係る計算値の変化量が推定

Xu, ³Numerical simulation of the ultrasonic

できれば，IM 法はより有益なものになると考えら

waves generated by ring-shaped laser illu-

れる，今後さらなる解析事例によるデータの蓄積がmination patterns,· Optics & Laser Technolo-望ましい．

gy, 39 [6], 1281-1287 (2007).
(11) A. Galdos, H. Okuda, G. Yagawa, ³Finite

文献

element simulation of ultrasonic wave propagation in pipe and pressure vessel walls.,·(１) 劉銘崇，三木千尋，³シャープな欠陥からの超音波端部エコーの発生挙動とその特性，· 土木学会論文集，525［I-33］，109-116（1995）．

Finite Element Analysis and design, 7 [1] ,1-13 (1990).
(12) F. Moser, L. J. Jacobs, J. Qu, ³Application(２) 大槻明，春海佳三郎，³地震時の大規模地下空of finite element methods to study transient

洞周辺の波動伝播機構，· 土木学会論文報告

wave propagation in elastic wave guides,·集，339，227-230（1983）．

Review of Progress in Quantitative Nondes-(３) K. M. A. Jaleel, N. N. Kishore, V. Sundar-tructive Evaluation, 17, 161-167 (1998).

arajan, ³Finite-element simulation of elastic(13) K. Mirkhani, C. Chaggares, C. Masterson,

wave propagation in orthotropic composite

M. Jastrzebski, T. Dusatko, A. Sinclair, R.

materials,· Materials Evaluation, 51 [7] ,J. Shapoorabadi, A. Konrad, M. Papini,830-838 (1993).

³Optimal design of EMAT transmitters,·(４) 伊藤智啓，川嶋紘一郎，森貴一郎，表竜二，

NDT & E International, 37 [3] , 181-193

221(2004).
(14) 吉川孝男，前田正広，岸田和人，中村栄隆，前田修平，³数値シミュレーションを用いた超音波伝播挙動の可視化研究，· 西部造船会 110 回例会論文梗概，119-125（2005）．
(15) A. K. Nair, P. R. Heyliger, ³Elastic waves
in combinatorial material libraries,· WaveMotion, 43 [7], 529-543 (2006).