衝撃応答スペクトル

SRS解析への考察B
衝撃試験解析手法
by ke-ta
SRS解析とは？

Shock Response Spectrum (SRS)
訳：衝撃応答スペクトル
ある供試体が受けた衝撃に対するダメージポテンシャルを評価解析す
る手法の１つ。
試験する供試体が、あらゆる固有振動数をもった物体から構成されて
いると仮定し、これを数学モデルに置き換え衝撃に対する応答を固有振
動数ごとに計算し、その最大値を周波数軸上にプロットしたもの。
SRSの概念

SRSは、データを数学モデルに当てはめる
入力
G
Aデータ
出力
t
m
G
G
Bデータ
t
G
k
c
Aデータ
Bデータ
Cデータ
f
Cデータ
各固有振動数ごとに最大加速度応
答が得られる。
t
別々の試験にて取得されたデータA,B,C
（供試体に加えられた衝撃結果）
このままでは、比較も出来無いし、どういっ
た結果を意味しているのかよくわからない。
そこで、ある数学モデルを仮定する。この
数学モデルに同等の衝撃を印加する（計
算上）。
等価比較が可能となる。
数学モデル上でのピーク等が分か
る。
供試体のダメージ評価ができる
FFTとの違い？

FFTは周波数の成分分解
入力
出力
入力した波形の周波数成分
に分解する。逆FFTをかけれ
ば、元に戻る。
FFT
t
f

SRSは数学モデルを仮定した周波数解析
入力
出力
入力した波形にて、数学モデ
ル上に衝撃を印加し、最大応
答を見る。不可逆計算となる。
SRS
t
f
数学モデル

バネ・マス・ダンパモデル
運動方程式
m
x
u
d 2x
m 2  k ( x  u )  c( x  u )
dt
y  x u
k
c
固有角振動数
n 
衝撃印加
k
m
減衰比
1
 
2Q
Qは、ユーザが規定するパラメータ
H2Aユーザマニュアル等では、
Q=10となっている。
計算の流れ
1.Qを任意に決定する
2.解析したい数学モデルのωnを決定する
3.衝撃データを入力する
4.各時間ごとの加速度応答を求める
5.加速度応答の絶対最大値を求める
入力
Q=10
ωn=100
出力
この点（絶対最大値）が、ωn=100
時の加速度応答となる
t
t
周波数グラフを作るには、各ωnについ
て繰り返し解いていく。
計算例1
422K
Half Sin波の750GをSRSにかけたもの
422K@750g's HalfSin/ 0.5msec（Q=10）
1.E+04
1.E+03
1.E+02
Gsrs[g]

1.E+01
1.E+00
10
100
1000
1.E-01
Freq.[Hz]
10000
計算例2
とある試験との比較
422K@750g's HalfSin/ 0.5msec（Q=10）
422K規定
JPOD
1.E+05
1.E+04
1.E+03
Gsrs[g]

1.E+02
1.E+01
1.E+00
10
100
1000
10000
100000
1.E-01
Freq.[Hz]
スペックアウトしているような…？（＾＾;
どの衝撃が厳しいのか？1/2

Half Sin波
80G/3ms
50G/11ms
50G/6ms
80G
50G
3ms
50G
11ms
どの衝撃が厳しいのか？→SRS解析
6ms
どちらの衝撃が厳しいのか？2/2
Half Sin波
SRS衝撃ポテンシャル比較
80G/3ms
50G/11ms
50G/6ms
1.E+03
100Hz以上であれば、80G/3ms
の方が厳しい条件となる
100Hz以下に固有振動数をもつ
供試体であれば、50G/11msの
方が厳しい条件となる
Gsrs[g]

1.E+02
50G/11msと比べ立ち上がりが遅
くなる。
1.E+01
1
10
100
Freq.[Hz]
1000
10000
ソースコード

perl
#!/usr/bin/perl
# 以下は変更する必要なし
#----------------------------------------------------------------------------#
# csv2srs (version $Id: eodv,v 1.4 2007/05/22 02:22:33 fukuda Exp $)
#
# SRSの算出
#
# Usage:
# > csv2srs.pl hogehoge.csv
#
# options:
#
特に無し
#
# example:
# csv2srs.pl hogehoge.csv
#(1)
#
#----------------------------------------------------------------------------# Libraries
my @dt, @dg, @alpa, @beta, @A, @B, @C, @D, @yy;
my $tt,$bs,$bb,$ymax;
$Q
$PI
$g
$delt
$zi
$wn
$fbgn
$fmax
$df
$i
$cnt
$k
$m
$c
$oct
$aa
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
10;
3.1415926;
9.80665;
1.00E-06;
1/(2*$Q);
0.0;
1;
10000;
0;
0;
0;
0;
1;
1;
100;
0;
# [m/s2]
# でるｔt
# [Hz]
# [Hz]
# Δ[Hz]
#----------------------------------------------------------------------------# Code begin
#----------------------------------------------------------------------------($opt) = @ARGV;
if ($opt eq "-h") {
&show_usage();
exit;
}
# オプションがない場合
if ($#ARGV == -1){
}
# オプションが１つある場合
if ($#ARGV == 0){
print "Error:ファイル名を指定してください\n";
# 対象ファイルのオープン
open(IN,"<$ARGV[0]") || die "Error: can't open $ARGV[0]\n";
close(IN);
# SRS算出
$cnt = -1;
foreach (@line){
@line = <IN>;
# ヘッダ情報を飛ばす
if ($cnt == -1){$cnt++; next;}
# 時間[sec]／加速度取得[m/s2]
$dt[$cnt+1] = (split(/\,/,$_))[0];
$dg[$cnt+1] = (split(/\,/,$_))[1];
}
# インクリメント処理
$cnt++;
# ヘッダ出力
$disp = sprintf("f[Hz],G[g],Q=%f\n",$Q);
print $disp;
$aa = &log10(($fmax))/$oct;
for($j=$fbgn;$j<=$oct;$j++)
{
}
sub cal_acc()
{
}
# 計算
$wn = 2*$PI*(10**($j*$aa));
$ymax = &cal_acc();
# 結果出力
$disp = sprintf("%f,%e\n",$wn/(2*$PI),$ymax);
print $disp;
# 初期化
$C[0] = $D[0] = $ymax = 0;
# ヘッダ-debug
#$disp = sprintf("#,time,応答G,αn,βn,Cn,Dn,Cn-1,Dn-1,a'(tn),b'(tn),c'(tn),d'(tn),An,Bn,An-1,Bn-1,Mn(tn),Mn-1(tn)\n");
#print $disp;
for($i=1;$i<$cnt;$i++)
{
# 加速度計算
#$A[$i] = $B[$i] = $C[$i] = 0;
# α,β計算・例外
$alpa[$i] = (($dg[$i+1] - $dg[$i])/($delt))*$g;
$beta[$i] = (($dg[$i]*$dt[$i+1] - $dg[$i+1]*$dt[$i])/($delt))*$g;
# 特解計算
$C[$i] = -$alpa[$i]/($wn*$wn);
$D[$i] = -$beta[$i]/($wn*$wn) + 2*$zi*$alpa[$i]/($wn*$wn*$wn);
# 一般解計算
$tt = $dt[$i];
$bs = &fd($tt)*(&fM($i-1,$tt)-&fM($i,$tt))-&fb($tt)*($C[$i-1] - $C[$i]);
$bb = &fa($tt)*&fd($tt)-&fb($tt)*&fc($tt);
$A[$i] = $bs/$bb/exp(-$zi*$wn*$tt) + $A[$i-1];
$bs = -&fc($tt)*(&fM($i-1,$tt)-&fM($i,$tt))+&fa($tt)*($C[$i-1] - $C[$i]);
$B[$i] = $bs/$bb/exp(-$zi*$wn*$tt) + $B[$i-1];
# yy出力
#$yy[$i] = &fM($tt) + $A[$i]*&fa($tt) + $B[$i]*&fb($tt);
$yy[$i] = -(&fM($i,$tt)+($A[$i]*&fa($tt) + $B[$i]*&fb($tt))*exp(-$zi*$wn*$tt))*$wn*$wn/$g;
# ymaxの算出
if(abs($yy[$i]) > $ymax){$ymax = abs($yy[$i]);}
}
}
# 出力-debug
#$disp = sprintf("%d,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e\n",$i,$dt[$i],$yy[$i],$alpa[$i],$beta[$i],$C[$i],$D[$i],$C[$i-1],$D[$i-1],&fa($tt),&fb($tt),&fc($tt),&fd($tt),$A[$i],$B[$i],$A[$i-1],$B[$i-1],&fM($i,$tt),&fM($i-1,$tt));
#print $disp;
return $ymax;
# オプションが２つある場合
if ($#ARGV == 1){
}
exit;
#----------------------------------------------------------------------------# Sub code begin
#----------------------------------------------------------------------------sub fM()
{
}
sub fa()
{
return $C[$_[0]]*$_[1] + $D[$_[0]];
if($_[0] == 0){
}else{
}
sub fb()
{
}
if($_[0] == 0){
}else{
}
sub fc()
{
}
sub fd()
{
}
sub log10 {
my $n = shift;
return log($n)/log(10);
}
sub show_usage()
{
print "EoDV(放電末期電圧)を抽出・可視化する\n";
print "\n";
print "Usage:\n";
print "> eodv [-h]\n";
print "\n";
print "options:\n";
print "
-h : 利用方法を提示する\n";
print "
-make : Level-2を更新する\n";
print "\n";
print "example:\n";
print "eodv -h
#(1) HELP\n";
}
}
return 1;
#return exp(-$zi*$wn*$_[0])*cos(sqrt(1-$zi*$zi)*$wn*$_[0]);
return cos(sqrt(1-$zi*$zi)*$wn*$_[0]);
return 0;
#return exp(-$zi*$wn*$_[0])*sin(sqrt(1-$zi*$zi)*$wn*$_[0]);
return sin(sqrt(1-$zi*$zi)*$wn*$_[0]);
#return $wn*exp(-$zi*$wn*$_[0])*(-$zi*cos(sqrt(1-$zi*$zi)*$wn*$_[0])-sqrt(1-$zi*$zi)*sin(sqrt(1-$zi*$zi)*$wn*$_[0]));
return $wn*(-$zi*cos(sqrt(1-$zi*$zi)*$wn*$_[0])-sqrt(1-$zi*$zi)*sin(sqrt(1-$zi*$zi)*$wn*$_[0]));
#return $wn*exp(-$zi*$wn*$_[0])*(-$zi*sin(sqrt(1-$zi*$zi)*$wn*$_[0])+sqrt(1-$zi*$zi)*cos(sqrt(1-$zi*$zi)*$wn*$_[0]));
return $wn*(-$zi*sin(sqrt(1-$zi*$zi)*$wn*$_[0])+sqrt(1-$zi*$zi)*cos(sqrt(1-$zi*$zi)*$wn*$_[0]));
計算反省点





Excelにてcsv保存すると表示されている有
効数字に短縮される
参考PDFと参考Excelの計算式が異なる。
ただし、解いているものは一緒
Expを掛け算にしていた（式ミス）。
割と計算の有効桁数がシビアに効いてくる
ため、定数の桁にも注意が必要
Perlには、常用対数log10がない
底変換の式にて代用

Download Report