放射音測定と平板の曲げ理論を用いた衝撃力の同定
Identification of the Impact Force by using Acoustic Response and Bending Theory of Plate
精密工学専攻
6
号 井坂隆之 Takayuki Isaka1.
緒言緒言緒言緒言物体の破損に関係する衝撃力を測定することは重要であ り,これまでに衝撃力を同定する研究が多数報告されている.
例えば,井上らは,衝撃力が加えられた弾性体に生じるひず み応答を用いて衝撃力の同定を行っている(1)(2).
一方,衝撃力により発生する放射音にはその衝撃力の情報 が含まれているはずである.放射音の測定による衝撃力の同 定が出来れば,衝撃体や被衝撃体から離れた位置に置いた測 定素子による非接触測定が可能であり,また,ひずみゲージ による測定のように測定素子の貼り付け等の手間がかから ない.辻らにより,衝撃力と放射音の伝達関数を用いて,測 定した放射音より衝撃力を同定する研究が報告されている(3). この方法では,平板の曲げ理論を用いることで平板の変位か ら音圧を計算し,衝撃力と放射音の関係を求めている.さら に,逆問題解析手法を用いて,実測した放射音から衝撃力の 時間変化を同定している.そこでは中実円柱が平板に衝突し た際の衝撃力を,測定した放射音を用いて同定している.し かし,実用性を考慮すれば衝撃力が一点では無く,分布する ような場合も考える必要がある.
そこで本研究では,先の方法の改良を行い,衝撃力が分布 する問題としてリング状領域の荷重が平板に与えられた際 の衝撃力の時間変化を同定する手法を提案する.そして,実 際の衝撃実験において,提案している手法により,衝撃力の 時間変化の同定が可能であることを示し,本手法の有効性を 明らかにする.
2.
同同同同定原理定原理定原理定原理本研究では音圧を測定するマイクロフォンを設置し,測定 される音圧の時間変化のデータを用いて衝撃力の時間変化 の同定を行う.そのために,平板の曲げ理論を用いて音圧を 計算し,衝撃力と放射音の関係を求めている.このとき,与 えた衝撃力の時間変化をP(t) ,求められた音圧をS(t)をとす る.
任意の衝撃力の時間変化Pcal (t)は,P(t)の重ね合わせにより 式(1)のような未知係数Anの級数表示で表すことができる.
∑
=−
= N
n
n n
cal t A P t t
P
1
) ( )
( (1)
ここで,tn= ∆t×(n-1),∆tは時間刻み,Nは任意の重ね合わ せ回数である.
P(t)が加わったときの,音圧の時間変化を平板の曲げ理論 を用いて計算し求める.衝撃力と放射音圧には線形関係があ る.すなわち,衝撃力がAn倍となったとき,放射音圧もAn 倍となる.したがって,式(1)で表される衝撃力Pcal (t)に対す る放射音圧の時間変化Scal (t)は,S(t)を用いて以下のように表 される.
∑
=−
= N
n
n n
cal t A S t t
S
1
) ( )
( (2)
このとき実際に測定されたマイクロフォンの音圧の時間 変化Sexp(t)と式(2)の計算値Scal (t)の差の二乗和Eerrorは,次の ように表される.
2
1
exp 1
max
) ( )
∑ ∑
(= =
− −
=N
i
i n
i N
n n
error A S t t S t
E (3)
ここで,ti (i=1,2,...Nmax)は音圧を測定した時刻,Nmaxは測定時 刻の総数である.最小二乗法によれば,Anは次の連立一次方 程式の解として得られる.
) , , 2 , 1 ( ) ( ) (
) ( ) (
max max
1 exp
1 1
N m
t t S t S
t t S t t S A
N
i
m i i N
i
m i n i N
n n
= K
−
=
−
−
∑
∑∑
=
= =
(4)
式(4)よりAnが求まれば,式(1)より衝撃力の時間変化Pcal (t) が同定できる.Fig.1 に衝撃力の時間変化の同定までのフロ ーチャートを示す.
Fig.1 The flow chart to identify the impact force
3.
実験方法実験方法実験方法実験方法3.1 3.1 3.1
3.1 実験装置実験装置実験装置実験装置
Fig.2に測定システムの概略を示す.衝撃による平板からの 放射音圧は,マイクロフォンで測定され,デジタルオシロス コープで記録される.使用したマイクロフォンはBruel&Kjaer 製のType4190(直径12.7mm,長さ17.6mm)で,測定可能な 音圧変化の周波数帯域は3〜20kHzであり,各周波数に対して 更正された音圧値が電圧として出力される.辻らの報告によ るとマイクロフォンは,前面がキャップで保護された構造と なっているため,実際の音圧の測定素子はマイクロフォンの 表面から5mm内側に入っていると推定されている(4).マイク ロフォンの測定位置は,この位置を実際の位置としている.
またアクリル中空円柱側面に貼り付けたひずみゲージによ り,ひずみデータも同時に記録する.このデータより衝撃力 を算出し,同定により得られた結果と比較し,衝撃力の同定 結果の有効性を検討する.全ての記録されたデータはPCによ り処理を行う.
Fig.2 The Measuring System of the sound and the strain data oscilloscope
conditioning amp
microphone
signal conditioner bridge box
hollow cylinder strain gauge
plate aluminum bar acryl Impact
An : Eq.(4)
Pcal(t) : Eq.(1) Sexp(t)
Measure sound pressure Calculate relationship P(t)
S(t)
Fig.3,4に示すように,アルミニウムの円板(直径500mm,
厚さ10mm)の中心にアクリルの中空円柱(半径rm[mm],高 さ50mm,肉厚10mm)を設置する.このアクリルの中空円柱 にアルミニウムの丸棒(直径30mm,長さ300mm)を高さ30mm から自由落下させることで,円板に半径rm のリング状領域の 荷重を与える.これは,リング状領域の荷重を衝突して与え ることが困難であるため,設置した状態で丸棒を落とすこと により荷重を与えている.アクリルの中空円柱の先端から 15mmの所にひずみゲージを貼り付け,中空円柱の側面に生 じるひずみ変化から衝撃力を測定する.この測定値は衝撃力 の同定の精度の確認のために用いる.アルミニウム丸棒はリ ニアブシュをガイドして垂直に落ちるようにしている.円板 の下に厚さ5mm,一辺20mmの正方形のシリコンラバーシー トを設置することにより支持している.これは,円板に外乱 の力を作用させず,なるべく自由に振動させるためである.
円板の下面から下方10mmの位置にマイクロフォンを設置し,
音圧を測定する.
Fig.5にrm=20, 40mmの場合の音圧の時間変化Sexp(t)を示す.
aluminum plate
silicone rubber
sheet aluminum bar
30mm
Φ500mm 300mm Φ30
10mm
10mm microphone
50mm
guide
strain gauge
acryl hollow cylinder
Fig.3 The Experiment equipment
50 30
15 10
rm
dout
din
Fig.4 The configuration of the acryl hollow cylinder
-100 -50 0 50 100 150 200
0 0.2 0.4 0.6 0.8 1
t [ms]
Sound pressure [Pa]
rm=40mm 20mm
Fig.5 The time dependence of the measured sound pressure with rm = 20, 40mm
3.23.2
3.23.2 衝撃力の測定衝撃力の測定衝撃力の測定衝撃力の測定
アクリル中空円柱の断面内の応力分布が一様であると仮 定すれば,アクリル断面に加わる力の時間変化Pexp(t)は,ア クリルに貼付けた半導体ゲージより得られたひずみの時間
変化εhollow(t)より次のように与えられる.
) 4(
)
( 2 2
exp t Ehollow hollow dout din
P = ε π − (5)
ここで,Ehollow ,dout ,din はそれぞれ,アクリル中空円柱の 縦弾弾性係数と外径,内径である.本研究では,この値を円 板上面に加わる衝撃力の時間変化の参考値として用いてい る.
4.
平板の曲げ理論平板の曲げ理論平板の曲げ理論平板の曲げ理論を用いた放射音圧の計算を用いた放射音圧の計算を用いた放射音圧の計算 を用いた放射音圧の計算平板の曲げ理論を用いて,平板の変形から放射音圧を計算 する.まず,平板の中心に衝撃力が加わった際の任意の位置 での放射音圧を計算する.次に,得られた放射音圧を積分す ることでリング状領域に衝撃力が加わった際の放射音圧を 求める.
impact distribution h
r z
f(r,t)=P0H(t) a2 e
r2 a2
Fig.6 The plate with impact force
Fig.6に示すような厚さhの平板の上表面に軸対称の衝撃分 布荷重 f(r, t)が加わる問題を平板の曲げ理論を用いて計算す る.分布荷重 f(r, t)を以下のように置く.
f(r,t)=P0
H(t) πa2 e
−r2
a2 (6)
ここで,aは衝撃力の分布半径,H(t)はヘビサイドのステップ
関数.P0は分布荷重の合計,すなわち衝撃力である.r = 0に おいて,分布荷重は最大でr = aにおいて,1/eとなり,またa
→0において,f(r ,t)→P0H(t)δ(r)である.平板のz方向の時間変 化w(r, t)は以下の微分方程式を満足する(5).
∇4w+ 1 cb2
∂2w
∂t2 = f(r,t)
D (7)
∇2≡ ∂2
∂r2 +1 r
∂
∂r, 1 cb2 =
ρ
bhD , D= Eh3 12 (1−ν2)
ここで,D,E,ρはそれぞれ,平板の曲げ剛性,縦弾性係数,
密度を表す.
ハンケル,ラプラス変換を以下のように示す.
∫
∞= 0 ( , ) 0( ) )
,
~(
dr r rJ t r f t
f ξ ξ (8)
∫
∞ −= 0 ( , ) )
,
~(
dt e t r f p
r
f pt (9)
式(6)はハンケルラプラス変換領域において以下のように 表される.
0 4
2 2
1 ) 2 ,
~( ξ
ξ π
a
pe p P
f = − (10)
式(10)を用いて式(7)をw~について解く.
4 2 4
2 2 0
2 2
) (
1 ) 2
,
~( ξ
π ξ ξ
a
b
b e
c p D p c p P
w −
= + (11)
さらに,式(11)に逆ハンケル,ラプラス変換を行う.平板のz 方向の変位w(r, t)は次のようになる.
w(r,t)= P0 2πD
1−cos(cbξ2t) ξ3 e
−1 4a2ξ2
J0(rξ)dξ
0
∫
∞ (12)式(12)を時刻tで偏微分することで,平板のz方向の加速度 が次式のように得られる.
&&
w(r,t)= P0cb2
2πD ξcos(cbξ2t)e−
1 4a2ξ2
J0(rξ)dξ
0
∫
∞ (13)さらに右辺の積分を実行すれば,平板表面の加速度は次式に 示すような閉じた形で求めることができる.
&&
w(r,t)2πDa2 P0cb2 =2e
−(r*)2 1+(t*)2
1+(t*)2 t*sin(t* (r*)2 1+(t*)2)
+cos(t* (r*)2
1+(t*)2)
(t*=4c
b
a2 t , r*= r
a , ξ*=aξ) (14)
y
x
z
) , (rt w&&
dA R )
, (rt f
θ1
r r1 impact distribution
X(x,y,z)
sound pressureS(t,X) Fig.7 The relation between location of the sound pressure and
acceleration of the region on the plate
Fig.7に示すように,平板上の微小部分dAがz軸方向に加速
度w&&(r,t)で運動する時,平板下の位置X(x,y,z)での音圧dSは次
式に表される(6).
c dA t R r w R X
t
dS ( , )
) 2 ,
( = && −
π
ρ (15)
R= r12+z2,r= (x+r1cosθ1)2+(y+r1sinθ1)2 ここで,ρ,c はそれぞれ空気の密度,音速を表す.式(15) を平板の表面で積分することで位置Xでの放射音圧S(t, X)は,
∫ ∫
− −−
= π θ
π ρ
0 0 1 1
1 2 2 2
) , ( )
( 2 ) ,
( tc z drd
c t R r Rw z r
tc H X t
S && (16)
R= r12+z2,r= (x+r1cosθ1)2+(y+r1sinθ1)2 となり,式(14)を代入して計算することができる.
y’
x’
z
X’(x’,y’,z)
θ
r
myx
sound pressureS'(t,X')
Fig.8 The relationship between location of the sound pressure and impact force applied to circular region
Fig.8に示すように,得られた放射音圧S(t, X)を数値積分す
ることで,リング状領域に衝撃力が加わった際の放射音圧 S’(t, X’)は,
∫
= 2π θ
0 ( , )
) ' , (
' t X S t X d
S (17)
(
X'(x',y',z)=X(x'−rmcosθ,y'−rmsinθ,z))
と表すことができる.
実 験 条 件 に 合 わ せ , 各 定 数 はcb =2202s/m2, E=70GPa, ρ=1.3kg/m3 ,c = 340m/sとしてS'(t, X’)を求める.
Fig.8にX’(0,0,10)mmとし,rm を変化させた場合の音圧の時 間変化を示す.rm が大きくなるにつれ音圧の最大値は小さく なり,生じる時刻も遅くなっている.また,rm = 20, 40, 80mm の波形は立ち上がる際に振動していることが分かる.ただし,
これらの波形は全て衝撃力を一定にした際の比較である.
Fig.9にrm = 40mmとし,音圧の位置をx'方向にを移動させた
場合の音圧の時間変化を示す.中心から離れるにつれ音圧の 最大値は小さくなることが分かる.また,衝撃位置に最も近
い位置x' = 40mmで最も早く立ち上がり,x' = 0mmでは負に立
ち上がることが分かる.
-0.05 0 0.05 0.1 0.15
0 0.1 0.2 0.3 0.4 0.5
t [ms]
Sound pressure [Pa]
rm=0mm 20mm 40mm
80mm
Fig.9 The time dependence of the sound pressure with x' = 0mm, y' = 0mm, z = 10mm, rm = 0, 20, 40, 80mm,dθ=2π/64
-0.05 0 0.05 0.1 0.15
0 0.1 0.2 0.3 0.4 0.5
t [ms]
Sound pressure [Pa]
x’=0mm
40mm 20mm
60mm
Fig.10 The time dependence of the sound pressure with x' = 0,20, 40, 60mm, y' = 0mm, z = 10mm, rm = 40mm,dθ=2π/64
5.
同定同定同定同定結果結果結果結果Fig.11にrm=20mmの衝撃実験で測定した放射音圧の時間変
化Sexp(t)と,式(2)から求めた音圧の時間変化Scal (t)を示す.
Scal(t)はt = 0~1msを200分割,すなわちNmax = 200, ∆t =0.005ms とし,式(4)の最小二乗法により未知係数Anを求めている.
Scal(t)は,N = 10, 50, 200と変化させてある.N が大きくなる ほど測定値に近くなることが分かる.
Fig.12に測定した衝撃力の時間変化Pexp(t)と式(1)から求め
た衝撃力の時間変化の同定結果Pcal(t)を示す.式(5)より,
Pexp(t)はEhollow = 3.14GPa, dout = 90mm, din = 70mmとして計算 している.Pcal (t)は,N = 200の値を示してある.衝撃力の測 定値と同定値の立ち上がり,最大値は共にほぼ一致している ことが確認できる.しかし,測定値が0.04ms早く立ち上がっ ている.また,最大値については同定値が15%大きくなって いる.
同様にFig.13にrm=40mmの衝撃実験で測定した放射音圧の
時間変化Sexp(t)と,式(2)から求めた音圧の時間変化Scal (t)を示 す.N が大きくなるほど測定値に近くなることが分かる.
Fig.14に測定した衝撃力の時間変化Pexp(t)と式(1)から求め
た衝撃力の時間変化の同定結果Pcal (t)を示す.式(5)より,
Pexp(t)はEhollow =3.14GPa, dout = 50mm, din= 30mmとして計算し ている.Pcal (t)は,N = 200の値を示してある.衝撃力の測定 値と同定値の立ち上がり,最大値は共にほぼ一致しているこ とが確認できる.しかし,測定値が0.01ms早く立ち上がって いる.また,最大値については同定値が14%大きくなってい る.
これらの衝撃力の測定値と同定値の誤差については,測定 値の測定位置が衝撃点でなく,側面のひずみから衝撃力を計 算していることが原因だと考えられる.また,0.5ms 以降,
同定値は大きく振動してしまっている.これは,実験と平板 の曲げ理論での平板の境界条件が異なっていることが原因 であり,音と力の関連性が異なっているためであると考えら れる.
しかし,rm の大きさを変化させた場合でも,同定値と実測 値の波形の立ち上がり,最大値共にほぼ一致している.この 結果から,リング状領域に加わった衝撃力の同定は可能であ ることが確認できた.
6.
結論結論結論結論本研究では,平板の曲げ理論を用いてリング状領域に加わ る衝撃力と放射音の関連性を利用して,測定した放射音より 非接触に衝撃力を求める手法を提案した.実際に衝撃実験を 行い本手法による衝撃力の同定が可能であることを確認し た.以下に,本研究により得られた結果を示す.
(1) 従来の研究では,一点で加えられた衝撃力の同定であっ たが,リング状領域に加えられた衝撃力の場合でも同定 は可能であることが分かった.
(2) 平板の曲げ理論を用いることで,平板の変位から音圧を 計算することができる.このことから,リング上領域に 加わる衝撃力と放射音の関係を得ることができる.
参考文献 参考文献 参考文献 参考文献
(1)井上 裕嗣,岸本 喜久雄,渋谷 寿一,小泉 尭,逆解析 による衝撃荷重の推定,逆解析のための最適伝達関数, 日本機械学会論文集.A編,57-543 (1991) pp.2727-2734.
(2)井上 裕嗣,石田 宏之,岸本 喜久雄,渋谷 寿一,逆解 析手法を用いた衝撃荷重の測定,伝達関数の推定法の比 較検討と計装化シャルピー衝撃試験への応用,日本機械 学会論文集 A編,57-534 (1991) pp.424-429.
(3) Tomoaki Tsuji,Yuya Tominaga,Takayuki Isaka,Yuki Tsuchiya,Determination of the time dependence of the impact force by measuring the radiated sound from the impacted thin plate,ACMFMS2012,(2012).
(4)辻 知章,栗本 貴文,渋谷 寿一,FEM解析を併用した放 射音圧測定による平板の衝撃力の同定,日本機械学会論 文集.A編,74-742 (2008) pp.858-863.
(5)中原 一郎,渋谷 寿一,土田 栄一郎,笠野 英秋,辻 知 章,井上 裕嗣,弾性学ハンドブック,朝倉書店,東京
(2001),pp.516-517.
(6)Philip, M., Morse, K., and Uno, Ingard,. Theoretical acoustics, Princeton University Press (1986), P.387.
-100 -50 0 50 100 150 200
0 0.2 0.4 0.6 0.8 1
t [ms]
Sound pressure [Pa]
measured calculated (N=10) calculated (N=50) calculated (N=200)
Fig.11 The time dependence of the sound pressure (rm =20mm)
-500 0 500 1000 1500
0 0.2 0.4 0.6 0.8 1
t [ms]
Impact force [N]
measured calculated (N=200)
Fig.12 The time dependence of the impact force (rm =20mm)
-100 -50 0 50 100 150 200
0 0.2 0.4 0.6 0.8 1
t [ms]
Sound pressure [Pa]
measured calculated (N=10) calculated (N=50) calculated (N=200)
Fig.13 The time dependence of the sound pressure (rm =40mm)
-500 0 500 1000 1500
0 0.2 0.4 0.6 0.8 1
t [ms]
Impact force [N]
measured calculated (N=200)
Fig.14 The time dependence of the impact force (rm =40mm)