FEM 解析を利用した放射音による衝撃位置と衝撃力の同定

FEM 解析を利用した放射音による衝撃位置と衝撃力の同定

An Identification Method of the Impact Position and the Time Dependence of Impact Force by Using Acoustic Response and FEM Analysis

精密工学専攻 19 号 栗本 貴文

Takafumi Kurimoto

1. 緒言

衝撃問題を考慮するうえで，物体の破損に関係する衝撃力 を測定することは重要であり，これまでに衝撃力を同定する 研究が多数報告されている．例えば井上らは，衝撃力を加え られた弾性体に生じるひずみ応答を用いて衝撃力の同定を 行っている^(1)～(3)．この方法は衝撃力と被衝撃体のひずみ応答 との関係が畳み込み積分で表されることを利用し，衝撃力と ひずみ応答との伝達関数を用いて，ひずみゲージより測定し たひずみ応答から衝撃力の同定を行っている．そのため，予 め実験により衝撃力とひずみ応答との伝達関数を求める必 要がある．しかし，衝撃力の正確な測定には困難を伴う事が 多い．またひずみゲージによる測定には，衝撃力を直接測定 できない，衝突による破損の恐れ，測定対象への貼り付けの 手間がかかる，などの欠点がある．

一方，衝撃力により放射される衝撃音には，その衝撃力の 情報が含まれているはずである．放射音の測定による衝撃力 の同定が出来れば，衝撃体や被衝撃体から離れた位置に置い た測定素子による非接触測定が可能である．筆者らにより，

衝撃力と放射音との伝達関数を用いて，測定した放射音より 衝撃力を同定する研究が報告されている⁽⁴⁾．この方法の場合，

衝撃力と放射音との伝達関数を予め実験により求める必要 がある．しかし，先に記述したように衝撃力の正確な測定に は困難を伴う事が多い．そこで筆者らにより，有限要素法解 析より求めた衝撃力と放射音との関連性を利用して，測定し た放射音より衝撃力を同定する研究が報告されている⁽⁵⁾．こ の方法ではFEM解析により求めた任意の衝撃力に対する音 圧のシミュレーション値を基に，未知係数を用いて未知の衝 撃力と音圧を級数表示する．この未知の音圧と測定した音圧 を用いて，最小二乗法より未知係数を決定し，衝撃力を同定 している．そのため衝撃力と放射音との伝達関数を求める予 備実験を行う必要がない．しかし，

FEM

解析により衝撃力と 放射音との関連性を求める際に正確な衝突した位置の情報 が必要であるが，実際の現場では正確な衝撃位置を特定する のは非常に困難である．

一方，測定した放射音には衝撃位置の情報も含まれている はずである．したがって，放射音の測定位置を増やせば，衝 撃位置を特定することも可能だと考えられる．そこで本研究 では，測定した放射音より衝撃位置と衝撃力を同定する手法 を提案する．

2. 同定原理

シミュレーションで加える任意の衝撃力の時間変化を

P

sim

(t)

，衝突した位置だと予想される任意の衝撃位置を

X

simと し，シミュレーションにおいて衝撃位置X_sim にP_sim

(t)を加え

る．そのとき

FEM

解析で求められる異なる位置

x

1と

x

2での二 つの音圧の時間変化のシミュレーション値をそれぞれp_sim1

(t)

と

p

sim2

(t)

とする．これらの関数を用いて未知の衝撃力

P(t)

とそ れに対する音圧の時間変化p₁

(t)とp

₂

(t)を以下のように級数表

示する．

∑

=

+

=

^N

n

A

n

P

sim

t t

n

t P

1

) ( )

( (1)

∑

=

+

=

^N

n

A

n

p

sim

t t

n

t p

1 1

1

( ) ( ) (2)

∑

=

+

=

^N

n

A

n

p

sim

t t

n

t p

1 2

2

( ) ( ) (3)

ここで，

t

n

= Δt

×

(n-1)

×

N

max

/N

，

Δt

は時間刻み，

A

nは未知係数，

Nは任意の重ね合わせ回数である．位置x

₁とx₂において測定し

た二つの音圧の時間変化をそれぞれ

p

exp1

(t)

と

p

exp2

(t)

とすれば，

測定値p_exp1

(t)と式(2)の計算値p

₁

(t)の二乗誤差E

_error1は，以下の ように表される．

2 1 exp 1

max

1 1

1

{

_sim

(

_{i n}

) (

_i

)}

N i

N

n n

error

A p t t p t

E = ∑ ∑ + −

= =

(4)

ここで，

t

i

(i = 1, 2, …N

max

)

は測定値から得られた時刻，

N

maxは 測定点の数である．最小二乗法より，

E

error1が最小となる未 知係数

A

nは以下の連立一次方程式の解として求められる．

) , 2 , 1 ( ) ( ) (

) ( ) (

max 1

1 1 exp

1 1

max

1 1

N m

t t p t p

t t p t t p A

N

i i sim i m

m i sim n i sim N

i N n

n

= K +

=

+ +

∑

∑ ∑

=

= =

(5)

式(5)から

A

nが求まれば，式(3)から

p

2

(t)を求めることができ

る．ここで，シミュレーションの衝撃位置が実際の衝撃位 置と同じであれば，

p

2

(t)と p

exp2

(t)は限りなく等しくなると考

えられる．しかし，シミュレーションの衝撃位置が実際の 衝撃位置と異なる位置ならば，その位置ずれの分だけ

p

2

(t)

に違いが表れる．したがって，測定値pexp2

(t)と式(3)のp

2

(t)

の二乗誤差

E

error2が最小となるときのシミュレーションの衝 撃位置Xsimが実際の衝撃位置だと考えられる．

2 2 exp 2

max

1 1

2

{

_sim

(

_{i n}

) (

_i

)}

N i

N

n n

error

A p t t p t

E = ∑ ∑ + −

= =

(6)

式(6)より求めたXsimのときのpsim1

(t)を用いてA

nを求めれば，

式

(1)

より測定値

p

exp1

(t)

に対する衝撃力の時間変化

P(t)

が 同定出来る．

3. 実験方法

3.1 実験装置の概略

Fig.1

に，測定システムの概略を示す．衝撃による，被衝撃

体からの放射音圧はマイクロフォンで測定し，デジタルオシ ロスコープで記録にする．また，衝撃棒側面に貼付けたひず みゲージより，ひずみデータも同時に記録する．記録された データは PCにより処理を行う．

Fig.2

に，衝撃体および被衝撃体の詳細形状を示す．衝撃体

として，アルミニウム丸棒（φ30mm×300mm）を，被衝撃 体としてアルミニウム平板（500mm×200mm×10mm）を用

いている．丸棒の先端側面は丸めてあり，先端は直径

10mm

の円型平面になっている．板の四隅の下に厚さ

5mm，一辺 20mm

の正方形のシリコンラバーシートを設置し支持してい る．これは，板に外乱の力を作用させず，なるべく自由に振 動させるためである．また，これにより板を支持している実 験装置からの影響も無視できると考える．

丸棒を高さ

H

から自由落下させることにより衝撃実験を行 なっている．板上方に内径

31mmのプラスチックの円筒を設

置し，円筒内に丸棒を通過させることにより板に落下するよ うにしてある．衝撃位置はFig.2 のように，板中心を原点，

板の長手方向を

X

方向とし，中心から

X

exp移動させて，実験を 行う．丸棒の先端から

20mmの所にひずみゲージを貼り，衝

撃棒に生じるひずみ変化から衝撃力を測定する．この測定値 は衝撃力の同定の精度の確認のために用いる．

3.2 放射音の測定

測定に使用したマイクロフォンは，Bruel＆Kjaer 社製の

Type4190（直径 12.7mm，長さ 17.6mm）で，測定可能な音圧

変化の周波数帯域は

3〜20kHz

であり，各周波数に対して更 正された音圧値が電圧として出力される．マイクロフォンは，

前面がキャップで保護された構造となっているため，実際の 音圧の測定素子はマイクロフォンの表面から若干奥に入っ た所にある．別途行った実験より，マイク表面から

5mm

内 側に入っていると推定した．以後，マイクロフォンの測定位 置は，この位置を実際の位置としている．

衝突の際の放射音圧を二つのマイクロフォンで測定する．

これらの二つのマイクを便宜上，マイク

1

とマイク

2

に区別 し，Fig.2 のように設置する．板裏面からマイクの音圧測定 部までの距離を

h

，板中心からの

X

方向の距離をそれぞれ

x

1と

x

2とする．今回の実験では，マイクの測定位置はh

= 10mm，

x

1

= -30mm

，

x

2

= 30mm

に固定して測定を行った．

Fig.1 The Measuring System of the sound and the strain.

Fig.2 The configuration of the impact bar and the aluminum plate.

3.3 衝撃力の測定

衝撃棒の断面内の応力分布が一様であると仮定すれば，衝 撃棒断面に加わる力の時間変化

P

_exp

(t)

は，衝撃棒に貼付け た半導体ゲージより得られたひずみの時間変化

ε

bar

(t)

より 次のように与えられる．

2 exp

( t ) E

_{bar bar}

( t ) 4 d

_bar

P = ε π (7)

ここで，Ebar，

d

bar はそれぞれ，衝撃棒の縦弾弾性係数と直 径である．本研究では，この値が板上面に加わる衝撃力の時 間変化とほぼ等しいと考え，衝撃力の直接測定値として用い ている．

4. シミュレーション

4.1 シミュレーションモデルの概要

Fig.3

にシミュレーションモデルを，Table 1と

2

に板と空

気のそれぞれの機械的性質等を示す．モデルは対称性より板 の

1/2

の領域を考え，板の奥行方向に

80mm，下方に 100mm

の空気層を置いて，構造と流体の三次元連性解析を行う．空 気層側面および下面において音波は反射されるので，この部 分の境界条件は実際と異なる条件となっている．しかし，測 定時間 0.5ms のときの音波の到達距離は 173mm であり，

この距離は，板から放射された音波が境界に反射して戻る距 離より短いため，測定時間

0.5msにおいては，音圧の計測位

置における反射音の影響は無いと考えられる．また，板の境 界条件はシリコンラバーシートの上に設置しているので非 常に複雑であるが，板の境界条件の影響が表れるのは

0.35ms

以降であることが確認されている ⁽⁵⁾．よって，衝突が非常に 短い時間内で行われているため，板の境界条件の影響が無視 できると考えられる．したがって，板には拘束条件を与えな いものとする．

4.2 シミュレーション結果

Fig.4

に

FEM

モデル

(

衝撃位置

X

sim

= 0mm)

の要素分割の様子 を示す．衝撃位置の直径d = 10mm の部分に一様に圧力を加 える事により衝撃力を加えた．

Fig.5

に，板に与えた圧力の 時間変化を示す．Fig.5 に示すように，この圧力は滑らかに 立ち上がる単位ステップ状の関数であり，立ち上がり時間は

0.1ms，最大圧力は 1MPaである．シミュレーションの衝撃力

はこの圧力に直径

d

の面積を乗じて求めている．

Fig.6

に板中心(X_sim

= 0mm)に圧力を加えたときの，解析に

より得られた

h = 10mm

，

x

1

= -30mm

と

x

2 ＝

30mm

の位置にお けるそれぞれの音圧の時間変化p_sim1

(t)とp

_sim2

(t)を示す． p

_sim1

(t)

と

p

sim2

(t)

は衝撃開始時間より若干遅れて立ち上がり，その後 上昇した後に緩やかに振動している．また，p_sim1

(t)とp

_sim2

(t)

の測定位置の中心に圧力を加えたので，

p

sim1

(t)

と

p

sim2

(t)

は完 全に一致している．

Fig.3 The configuration of the simulation model of the plate and

the air.

Table 1 The mechanical properties of the aluminum plate for the FEM analysis.

material Aluminum configuration 100×500×10 mm

Young’s modulus 70 GPa Poisson's ratio 0.34

density 2700 kg/m

³

N. of element 3705

Table 2 The mechanical properties of air for the FEM analysis.

material air configuration 180 × 500 × 110 mm

acoustic velocity 346 m/s

density 1.3 kg/m

³

N. of element 24020

Fig.4 FEM mesh for the plate and the air environment.

0 0.2 0.4 0.6 0.8 1 1.2

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time　ms

P re ssur e 　 MP a

Fig.5 The time dependence of the pressure as the input for the FEMsimulation (T = 0.1ms).

5. 同定結果

5.1 X

exp

=0mm での衝撃位置と衝撃力の同定

板中心(衝撃位置X_exp

= 0mm)に，高さH = 30mmから衝撃棒

を落下させ，衝撃実験を行う．その際に測定した放射音圧の 時間変化p_exp1

(t)をFig.7

に示す．t = 0～1msを

200

分割，すな

わち

N

max

= 200, Δt = 0.005ms

とし，式

(5)

の最小二乗法によ

り未知係数

A

_n を求めた．式(2) から求めた音圧の時間変化

p

1

(t)

もグレーの線を用いて同時に示してある．

p

1

(t)

は，

N = 20,

30, 40

と変化させてある．N = 40以上では測定値とほとんど

一致したので，以下の計算では

N = 40

として行っている．

Fig.8

に，衝撃位置X_exp

= 0mmのときのマイク 2

の測定値

p

exp2

(t)

とシミュレーションの衝撃位置

X

sim

= -20, 0, 20mm

のと きの式(3)のp₂

(t)との比較を示す．なお，X

_expは未知として扱 う．

Fig.8

より，

X

sim

= 0mm

のときの計算値

p

2

(t)

が測定値

p

exp2

(t)

に一致するので，実際の衝撃位置X_exp

= 0mmでないかと考え

られる．

より詳細に調べるために，式(6)の測定値p_exp2

(t)と計算値 p

2

(t)との二乗誤差E

error2の平方根をFig.9 に示す．ここで，シ ミュレーションの衝撃位置X_simは-100～100mmの範囲で，

5mm刻みにp

₂

(t)の計算を行う．また，板の境界条件の影響が

0.35ms以降表れるので，0～0.35msの範囲でE

error2の計算を行 う．Fig.9より，シミュレーションの衝撃位置X_simが実際の衝 撃位置Xexp

= 0mmに近付くにつれE

error2が小さくなっていき，

X

sim

= 0mmでは E

error2は最小値をとる．したがって，実際の衝 撃位置はXexp

= 0mmと特定することは十分可能である．

Fig.10

に，衝撃位置

X

exp

= 0mmのときの衝撃力の測定値

P

exp

(t)とシミュレーションの衝撃位置X

sim

= -20, 0, 20mmのと

きの式(1)の

P(t)との比較を示す．Fig.10

より，

X

sim

= 0mmのと

きの計算値P(t)が測定値Pexp

(t)に最も一致することが確認で

きる．この

P(t)の方が測定値 P

exp

(t)に比べて若干大きいが，波

形の立ち上がりは完全に一致し，

0.35msまでは測定値とはほ

ぼ一致しているので，本手法により衝撃力を同定することは 十分可能だと考える．0.35ms以降，P(t)が測定値に比べて大 きく振動しているのは板の境界条件の影響である．

-2 0 2 4 6 8 10 12 14 16 18 20

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time　ms

S o und pr es su re 　 Pa mike1(x1 = -30mm) mike2(x2 = 30mm)

Fig.6 The time dependence of the sound pressure by FEM simulation at the mike 1 and mike 2.

-100 -50 0 50 100 150 200

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time　ms

Sound pressure　Pa

measured calculated(N = 20) calculated(N = 30) calculated(N = 40)

Fig.7 The comparison between the measured and the calculated with various values of N.

-150 -100 -50 0 50 100 150 200 250 300

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time　ms

S o un d pr es su re 　 Pa

measured (Xexp = 0mm) calculated(Xsim = -20mm) calculated(Xsim = 0mm) calculated(Xsim = 20mm)

Fig.8 The comparison between the measured and the calculated at

the microphone 2 with the values of X

sim

= -20, 0, 20mm.

5.2 さまざまな衝撃位置での衝撃位置と衝撃力の同定

それぞれ衝撃位置Xexp

= 0, -5, -10, -20, -30, -40, -60, -80mmに，

高さH = 30mmから衝撃棒を落下させ，放射音を測定する．

Fig.11

に，衝撃位置Xexp

= 0, -5, -10, -20, -30, -40, -60, -80mmの

ときの，式(6)の二乗誤差

E

rror2の平方根を示す．ここで，シミ ュレーションの衝撃位置Xsimは-100～20mmの範囲で，

5mm刻

みに

p

2

(t)の計算を行う．Fig.12

より，

X

exp

= 0, -5, -10, -20, -60, -80mmのときは， X

simが実際の衝撃位置XexpにおいてError2が最 小値をとる．しかし

X

exp

= -30, -40mmでは， X

sim

= -35, -45mm

のときにError2が最小値をとり，実際の位置Xexpとはずれてい る．これは，衝撃位置がマイク

1

の測定位置から距離が離れ るにつれ，シミュレーションの衝撃位置と実際の衝撃位置と のずれが音圧に与える影響が減少し，

E

rror2の最小値を判断す るのが難しくなるからだと考える．また実験で半径

30mmの

衝撃棒を内径

31mmの円筒を通して落下させているので，衝

撃棒と円筒の

1mmの隙間によって衝撃位置にずれが生じて

いる可能性も考えられる．

X

exp

= -20mm以降は衝撃位置を正

確に判断するのは難しいが，±5mmの範囲で衝撃位置を求め ることは十分可能である．したがって，衝撃位置を板中心か らX方向に移動させても衝撃位置を同定することは十分可能 である．

Fig.12

に，衝撃位置Xexp

= 0, -5, -10, -20, -30, -40, -60, -80mm

での衝撃力の測定値

P

exp

(t)と計算値 P(t)との二乗誤差の平方

根を示す．P(t)は衝撃位置Xexpと同じ位置でのXsimのpsim1

(t)を

用いて，式

(5)

から

A

nを決定し，式

(1)

より求めている．

Fig.13

から，誤差の値はXexpごとにかなりばらつきがある．しかし 一定の誤差内にあるので，衝撃位置を板中心から

X

方向に移 動させても衝撃力を同定することは十分可能である．

0 10 20 30 40 50 60 70 80 90 100

-100 -80 -60 -40 -20 0 20 40 60 80 100 X

_sim　

mm

Square error of sound pressure　Pa

Fig.9 The relationship between the expected impact position X

sim

and the square error of the sound presser by the mike 2.

-500 0 500 1000 1500 2000

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time 　 ms

Imp act fo rce 　 N

measured (Xexp = 0mm) calculated(Xsim = -20mm) calculated(Xsim = 0mm) calculated(Xsim = 20mm)

Fig.10 The comparison between the measured and the calculated of the impact force with the values of X

sim

= -20, 0, 20mm.

6. 結論

本研究では，

FEM

解析により求めた衝撃力と放射音との関 連性を利用して，測定した放射音より非接触に衝撃位置と衝 撃力を同定する方法を提案した．実験より，本手法による衝 撃位置と衝撃力の同定が十分可能であることを確認した．以 下に，本研究により得られた結果を示す．

1.

本手法による衝撃位置と衝撃力の時間変化の同定が十分 可能であることを確認した．

2.

放射音の測定位置から離れた位置においても，±5mmの 誤差範囲で衝撃位置を特定することができた．

3.

放射音の測定位置から離れた位置においても，衝撃力の 時間変化の同定が十分可能であることを確認した．

参考文献

(1)

井上 裕嗣，岸本 喜久雄，渋谷 寿一，小泉 尭，逆解析による 衝撃荷重の推定，逆解析のための最適伝達関数, 日本機械学会論 文集.A編，57-543

(1991) pp.2727-2734.

(2)

井上 裕嗣，石田 宏之，岸本 喜久雄，渋谷 寿一，逆解析手法 を用いた衝撃荷重の測定，伝達関数の推定法の比較検討と計装 化シャルピー衝撃試験への応用，日本機械学会論文集.A編，

57-534

(1991) pp.424-429.

(3)

井上 裕嗣，池田 直人，岸本 喜久雄，渋谷 寿一，小泉 尭，衝 撃力の大きさと方向の逆問題解析，日本機械学会論文集.A編，

59-559

(1993) pp.572-579.

(4)

辻 知章，川田 祐嗣，鈴木 喜浩，山口 友康, 野田 直剛, 被衝 撃体からの放射音による衝撃力の同定，逆問題解析を用いた非 接触測定による同定実験

,

日本機械学会論文集

.A

編，65-632

(1999) pp.701-707.

(5)

辻 知章，栗本 貴文，渋谷 寿一，FEM解析を併用した放射音圧 測定による平板の衝撃力の同定，日本機械学会論文集

.A

編，

74-742

(2008) pp.858-863.

0 20 40 60 80 100 120 140

-100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 Xsim　mm

Square error of sound pressure　Pa

Xexp = 0mm Xexp = -5mm Xexp = -10mm Xexp = -20mm Xexp = -30mm Xexp = -40mm Xexp = -60mm Xexp = -80mm

Fig.11 The relationship between the expected impact position X

_sim

and the square error of the sound presser by the mike 2, when the impact position is X

exp

= 0, -5, -10, -20, -30, -40, -60, -80mm.

0 20 40 60 80 100 120 140 160 180

-80 -70 -60 -50 -40 -30 -20 -10 0 Xexp　mm

Square error of impact force　N

Fig.12 The relationship between the impact position X

_exp

and the square error of the impact force, when the impact position is X

exp

=

0, -5, -10, -15, -20, -30, -40, -50, -60, -80mm.