PSDデータから時刻歴の無相関ランダム入力を作成する方法
Altair EmployeeOverview
Altair MotionSolveなどの時刻歴の解析にてランダム入力解析を行いたい場合に、PSDデータから時刻歴の無相関ランダム入力を作成するスクリプトです。
Pre-Requisite
スクリプト、サンプルファイルのダウンロードはこちらです。
Usage/Installation Instructions
PSDデータの準備
入力ファイル：load_profile.csv
CSVファイルで横軸周波数、縦軸PSDで入力を作成してください。両対数グラフで補間されます。
今回は1Hzまではフラット、1Hz以上で振幅が低下する特性としました。
変換の設定
Altair Composeで添付のスクリプトstc_make_random_input_rev1.0.omlを開きます。
ファイルブラウザで開いているフォルダがカレントディレクトリとなります。csvファイルがあるフォルダを開いてください。
先に作成した入力ファイルを処理して、入力ファイルの特性を満足する時刻歴の無相関ランダム波形を作成します。
7行目で入力数を指定します。今回は自動車の4輪に無相関の入力を与えると仮定してNC=4とします。相関の無い4つの時刻歴の入力波形が作成されます。振幅は入力ファイルによりPSDで規定されていますので、位相を乱数rand関数を用いて設定し、逆FFTにより時刻歴波形を算出します。乱数を用いていますので、時刻歴波形はスクリプト実行のたびに変化します。4輪への入力は振幅は同じですが、位相が異なる状態となります。
11行目~で作成する時刻歴の時間を指定します。時間刻みはdt=0.005sec（200Hz）、時間は100secとしました。
16行目~で作成した入力に対する確認計算のため、オートパワースペクトル、クロスパワースペクトルを算出するためのサンプリング時間、オーバーラップ比を指定します。
開始をクリックすると、ランダム時刻歴波形を格納したload_profile_his.h3dと、そのランダム波形から算出したオートパワースペクトル、クロスパワースペクトル結果を格納したload_profile_fft.h3dが出力されます。
作成データの確認（時刻歴）
HyperGraphでload_profile_his.h3dを開くと、時刻歴波形を確認できます。
4種類の振動波形が格納されています。
作成データの確認（オートパワースペクトル、クロスパワースペクトル、伝達関数）
HyperGraphのcomplex plotでload_profile_fft.h3dを開くと、作成した時刻歴のランダム波形を再度周波数分析した結果を確認できます。
1)FFT
Loadcase 1にFFT処理した結果を格納しています。
0~100secのデータを一度にFFT処理しています。
Displacement : FFT処理した結果
PS Displacement : 複素共役をかけて、パワースペクトルに変換したもの
PSD Displacement : さらに周波数刻みで除算してパワースペクトル密度（PSD）に変換したもの
PSD Displacementを確認すると、PSDを定義したload_profile.csvと一致していることがわかります。時刻歴波形にDC成分は含まれませんので、0Hzの値はゼロです。
2)オートパワースペクトル密度、クロスパワースペクトル密度、伝達関数
Loadcase 2~4にオートパワースペクトル密度、クロスパワースペクトル密度、伝達関数の算出結果を格納しています。
2種類の信号AfとAを考えた場合、オートパワースペクトル密度APSD、クロスパワースペクトル密度CPSD、伝達関数TFはそれぞれ
となります。2種類の信号が異なる場合がクロスパワースペクトル密度、同じ場合はオートパワースペクトル密度でクロスパワーをオートパワーで割ったものが伝達関数です。
クロスパワースペクトル密度の処理パラメータはスクリプトの16行目～で設定しています。
切り出し時間を2sec = 周波数刻み0.5Hz
オーバーラップ率50%
ハニングウィンドウ
下記のように各切り出し時間に関してFFT処理を行い、2種類の信号で相関を計算します。その和を取ることで平均化処理が行われ、2種類の信号で相関のある信号のみが残り、相関の無い信号は相殺されて無くなります。
(a) オートパワースペクトル密度
Loadcase 2にオートパワースペクトル密度の結果が格納されています。pwelch関数で求めています。
こちらの結果もload_profile.csvと一致していることがわかります。
(b) クロスパワースペクトル密度
Loadcase 3にクロスパワースペクトル密度の結果が格納されています。cpsd関数で求めています。
クロスパワーは2つの信号の組み合わせですので、4 x 4の結果が格納されています。
下記に波形1を基準とした場合の波形1~4のクロスパワースペクトル密度を示します。
波形1の結果はオートパワースペクトル密度の結果と一致します。残りの波形2~4の結果は波形1にくらべ小さな値となります。これは、波形2~4は波形1と相関性が小さいことを表しています。無相関を仮定してランダムで位相を設定しましたので、妥当な結果と言えます。
波形2～4を基準とした場合も同様です。
(c) 伝達関数
Loadcase 4に伝達関数の結果が格納されています。クロスパワースペクトル密度をオートパワースペクトル密度で除算したものです。
伝達関数も2つの信号の組み合わせですので、4 x 4の結果が格納されています。
下記に波形1を基準とした場合の波形1~4の伝達関数を示します。
波形1の結果は自分自身の伝達関数ですので結果は1です。残りの波形2~4の結果は1に比べ小さな値となります。これは、波形2~4は波形1と相関性が小さいことを表しています。無相関を仮定してランダムで位相を設定しましたので、妥当な結果と言えます。理想的な無相関では0となります。
波形2～4を基準とした場合も同様です。
Post-Requisite
この無相関ランダム入力を用いて、MotionSolveの時刻歴解析を行い、計算結果を周波数分析したサンプルは下記で公開しております。あわせてご利用ください。
https://community.altair.com/community/?id=kb_article_view&sysparm_article=KB0121838
使用製品：Altair Compose、Altair MotionSolve/MotionView
よくあるエンジニアからの質問はこちら
Comments
-
この記事を参考に同じようにやってみましたが、
line 32のdata=csvread(filename_in,1,0);
のところで「Unknown function」が出てきてしまい波形を作成できません。
composeのバージョンによってはcsvreadが使えないのでしょうか？
composeのバージョンは2019.4を使用しております。
0