シミュレーションされるビームの固有振動数および減衰

単入力/単出力のハンマー励起

一連のハンマー打撃により 3DOF システムを励起し、センサーで結果の変位を記録します。システムは比例的に減衰され、減衰行列が質量および剛性行列の線形結合となります。

励起信号、応答信号、時間信号およびグラウンド トゥルース周波数応答関数を含む 2 セットの測定のデータをインポートします。応答信号の最初のセットY1は最初の質量の変位を測定し、2 番目のY2は 2 番目の質量を測定します。各励起信号は 10 個の連結されたハンマー打撃から構成され、各応答信号には対応する変位が含まれます。各打撃信号の持続時間は 2.53 秒です。励起および応答信号には加法性ノイズが存在します。最初の測定における最初の励起および応答チャネルを可視化します。

[t,fs,X1,X2,Y1,Y2,f0,H0] = helperImportModalData(); X0 = X1(:,1); Y0 = Y1(:,1); helperPlotModalAnalysisExample([t' X0 Y0]);

最初の励起および応答チャネルの FRF を動的柔軟性に関して計算してプロットしますが、これは力に対する変位の測定です [1]。既定の設定では、FRF はウィンドウが適用されたセグメントのスペクトルを平均して計算します。各ハンマー励起が実質的に減衰してから次の励起が開始するため、箱型ウィンドウを使用できます。センサーを変位として指定します。

winLen = 2.5275*fs;% window length in samplesmodalfrf(X0,Y0,fs,winLen,'Sensor','dis')

既定の'H1'推定器を使用して推定された FRF には、測定された周波数帯域に 3 つの顕著なピークが含まれ、振動の 3 つの柔軟なモードに対応します。コヒーレンスはこれらのピーク付近では 1 に近く、応答測定の S/N 比が低い反共振領域では低くなります。1 に近いコヒーレンスは推定の質が高いことを示します。'H1'の推定は、ノイズが出力測定にのみ存在する場合に最適ですが、'H2'推定器は、入力にのみ加法性ノイズが存在する場合に最適です [2]。この FRF の'H1'および'H2'推定値を計算します。

[FRF1,f1] = modalfrf(X0,Y0,fs,winLen,'Sensor','dis');% Calculate FRF (H1)[FRF2,f2] = modalfrf(X0,Y0,fs,winLen,'Sensor','dis','Estimator','H2');

著しい測定ノイズがある場合、または励起が不十分な場合、パラメトリック法を使用すると、データから FRF を正確に抽出する追加オプションが得られます。「部分空間」法は最初にデータ [3] に対して状態空間モデルを適用してから、周波数応答関数を計算します。状態空間モデルの次数 (極の数に等しい) および直達の有無を指定して、状態空間推定を設定できます。

[FRF3,f3] = modalfrf(X0,Y0,fs,winLen,'Sensor','dis','Estimator','subspace','Feedthrough',true);

ここで、直達項および 1:10 の範囲で最適な次数を含む状態空間モデルを適用することで、FRF3 が推定されます。'H1'、'H2'、および'subspace'の手法を使用して、推定された FRF と理論的な FRF を比較します。

helperPlotModalAnalysisExample(f1,FRF1,f2,FRF2,f3,FRF3,f0,H0);

推定器は応答ピーク付近では同等に機能しますが、'H2'推定器は反共振では応答を過大に推定します。コヒーレンスは推定器の選択の影響を受けません。

次に、ピーク選択アルゴリズムを使用して各モードの固有振動数を推定します。ピーク選択アルゴリズムは、シンプルで高速な手続きで FRF のピークを識別します。これはローカル メソッドであるため、各推定が単一の周波数応答関数から生成されます。また、単一自由度 (SDOF) 法でもあるため、各モードのピークが個別に考慮されます。その結果、モーダル パラメーターのセットが各 FRF に対して生成されます。前のプロットに基づいて、3 つのピークが含まれる 200 ～ 1600 Hz の周波数範囲を指定します。

fn = modalfit(FRF1,f1,fs,3,'FitMethod','PP','FreqRange',[200 1600])

fn = 1.0e+03 * 0.3727 0.8525 1.3707

固有周波数はおよそ 373、853 および 1371 Hz です。再構成された FRF をプロットし、modalfitを使用して測定したデータと比較します。FRF は、周波数応答関数の行列FRF1から推定されたモーダル パラメーターを使用して再構成されます。modalfitを出力引数を指定せずに再度呼び出し、再構成された FRF を含むプロットを生成します。

modalfit(FRF1,f1,fs,3,'FitMethod','PP','FreqRange',[200 1600])

再構成された FRF は、最初の励起および応答チャネルの測定された FRF と一致します。次の節では、2 つの追加の励起位置を検討します。

ロービング ハンマー励起

既定の'H1'推定器を使用し、3 つすべてのセンサーの応答について FRF を計算してプロットします。ロービング ハンマー励起であるため、測定タイプを'rovinginput'として指定します。

modalfrf(X1,Y1,fs,winLen,'Sensor','dis','Measurement','rovinginput')

前のセクションでは、単一の FRF から単一セットのモーダル パラメーターを計算しました。今度は、最小二乗複素指数 (LSCE) アルゴリズムを使用してモーダル パラメーターを推定します。LSCE および LSRF アルゴリズムは、複数の応答信号を同時に解析して、単一セットのモーダル パラメーターを生成します。これらはグローバルな多自由度 (MDOF) 法であるため、複数の周波数応答関数からすべてのモードのパラメーターが同時に推定されます。

LSCE アルゴリズムは計算モードを生成しますが、これは物理的に構造には存在しません。安定化ダイアグラムを使用し、モード数の増加に伴う極の安定性を調べて物理モードを特定します。物理モードの固有周波数および減衰比は、同じ位置にとどまる (すなわち安定している) 傾向があります。安定化ダイアグラムを作成し、周波数が安定している極の固有周波数を出力します。

[FRF,f] = modalfrf(X1,Y1,fs,winLen,'Sensor','dis','Measurement','rovinginput'); fn = modalsd(FRF,f,fs,'MaxModes',20,'FitMethod','lsce');% Identify physical modes

既定では、極の固有振動数の変化が 1% 未満の場合、極は周波数が安定しているとして分類されます。周波数が安定している極は、減衰比の変化が 5% 未満の場合、さらに減衰が安定しているとして分類されます。両方の基準を異なる値に調整できます。安定した極の位置に基づいて、固有周波数として 373、852.5 および 1371 Hz を選択します。これらの周波数は、他の周波数が安定した極の固有周波数と共に、modalsdの出力fnに含まれます。物理的に存在するモード数より高いモデル次数は、一般に、LSCE アルゴリズムを使用して良好なモーダル パラメーターの推定を生成するために必要です。この場合、4 つのモードのモデル次数で 3 つの安定した極を示します。目的の振動数はfnの4 行目の最初の 3 列に出現します。

physFreq = fn(4,[1 2 3]);

固有周波数と減衰を推定し、再構成および測定された FRF をプロットします。安定性ダイアグラム'PhysFreq'から決定される 4 つのモードおよび物理周波数を指定します。modalfitは、指定したモードのモーダル パラメーターのみを返します。

modalfit(FRF,f,fs,4,'PhysFreq',physFreq)

[fn1,dr1] = modalfit(FRF,f,fs,4,'PhysFreq',physFreq)

fn1 = 1.0e+03 * 0.3727 0.8525 1.3706 dr1 = 0.0008 0.0018 0.0029

次に、FRF を計算し、センサー位置が異なる 2 番目のセットのハンマー打撃の安定化ダイアグラムをプロットします。安定化基準を周波数の場合は 0.1%、減衰の場合は 2.5% に変更します。

润扬悬索桥[f] = modalfrf (X2, Y2、fs winLen,'Sensor','dis','Measurement','rovinginput'); fn = modalsd(FRF,f,fs,'MaxModes',20,'SCriteria',[0.001 0.025]);

より厳密な基準にすると、大部分の極が周波数が不安定として分類されます。周波数および減衰が安定している極は平均された FRF とかなりの程度で一致しており、これらが測定データに存在することが示唆されます。

physFreq = fn(4,[1 2 3]);

この測定のセットにおけるモーダル パラメーターを抽出し、測定の最初のセットにおけるモーダル パラメーターと比較します。励起および応答測定が一致する位置に対応する駆動点 FRF のインデックスを指定します。固有周波数の差は 1% 以内、減衰比の差は 4% 未満のため、モーダル パラメーターは測定と測定の間で一致することが示されます。

[fn2,dr2] = modalfit(FRF,f,fs,4,'PhysFreq',physFreq,'DriveIndex',[1 2])

fn2 = 1.0e+03 * 0.3727 0.8525 1.3705 dr2 = 0.0008 0.0018 0.0029

Tdiff2 = table((fn1-fn2)./fn1,(dr1-dr2)./dr1,'VariableNames',{'diffFrequency','diffDamping'})

Tdiff2 = 3x2 table diffFrequency diffDamping _____________ ___________ 2.9972e-06 -0.031648 -5.9335e-06 -0.0099076 1.965e-05 0.0001186

モーダル パラメーター推定に対してパラメトリック法を使用すると、FRF に測定ノイズがある場合や FRF が高いモーダル密度を示している場合にピーク選択および LSCE の手法として役立つ代替方法が得られます。最小二乗有理関数 (LSRF) アプローチは、分母を共有する伝達関数を多入力、多出力 FRF に適用し、それによって単一のモーダル パラメーターのグローバルな推定 [4] を取得します。LSRF のアプローチを使用した手順は、LSCE と類似しています。安定化ダイアグラムを使用して静止モードを識別し、識別した物理周波数に対応するモーダル パラメーターを抽出できます。

[FRF,f] = modalfrf(X1,Y1,fs,winLen,'Sensor','dis','Measurement','rovinginput'); fn = modalsd(FRF,f,fs,'MaxModes',20,'FitMethod','lsrf');% Identify physical modes using lsfrphysFreq = fn(4,[1 2 3]); [fn3,dr3] = modalfit(FRF,f,fs,4,'PhysFreq',physFreq,'DriveIndex',[1 2],'FitMethod','lsrf')

fn3 = 372.6832 372.9275 852.4986 dr3 = 0.0008 0.0003 0.0018

Tdiff3 = table((fn1-fn3)./fn1,(dr1-dr3)./dr1,'VariableNames',{'diffFrequency','diffDamping'})

Tdiff3 = 3x2 table diffFrequency diffDamping _____________ ___________ -7.8599e-06 0.007982 0.56255 0.83086 0.37799 0.37626

パラメトリック法についての最後のノートは、次のとおりです。FRF 推定法 ('subspace') およびモーダル パラメーター推定法 ('lsrf') は、動的モデルを時間領域信号や周波数応答関数に対して適用するために System Identification Toolbox で使用するものと類似しています。このツールボックスを使用できる場合、tfestやssestなどのコマンドを使用して、データを適用するためのモデルを識別できます。compareおよびresidコマンドを使用して、識別されたモデルの品質を評価できます。モデルの品質を検証すると、それらをモデル パラメーターの抽出に使用できます。これは、状態空間推定器を使用して、簡潔に示されます。

Ts = 1/fs;% sample time% Create a data object to be used for model estimation.EstimationData = iddata(Y0(1:1000), X0(1:1000), 1/fs);% Create a data object to be used for model validationValidationData = iddata(Y0(1001:2000), X0(1001:2000), 1/fs);

直達項を含む 6 次の連続時間状態空間モデルを識別します。

sys = ssest(EstimationData, 6,'Feedthrough', true)

sys = Continuous-time identified state-space model: dx/dt = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x4 x5 x6 x1 4.05 -1765 149.8 -1880 -49.64 -358 x2 1764 -0.3332 2197 -232.5 -438.3 -128.4 x3 -152.4 -2198 2.85 4715 255.9 547.5 x4 1879 228.2 -4713 -15.9 -1216 -28.79 x5 59.42 440.9 -275.5 1217 35.05 -8508 x6 363.7 120.2 -545.4 -44.02 8508 -92.45 B = u1 x1 -0.1513 x2 -1.911 x3 4.439 x4 -3.118 x5 -0.9416 x6 -8.039 C = x1 x2 x3 x4 x5 x6 y1 3.135e-05 2.511e-06 8.634e-06 -1.416e-05 2.218e-06 -6.271e-06 D = u1 y1 7.564e-09 K = y1 x1 3.513e+07 x2 -3.244e+06 x3 -3.598e+07 x4 -1.059e+07 x5 1.724e+08 x6 7.521e+06 Parameterization: FREE form (all coefficients in A, B, C free). Feedthrough: yes Disturbance component: estimate Number of free coefficients: 55 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using SSEST on time domain data "EstimationData". Fit to estimation data: 99.3% (prediction focus) FPE: 1.235e-16, MSE: 1.189e-16

検証データに対してどれほど適切に適用されるかをチェックすることで、モデルの品質を評価します。

clf compare(ValidationData, sys)% plot shows good fit

モデルsysを使用してモーダル パラメーターを計算します。

[fn4, dr4] = modalfit(sys, f, 3);

風力タービン ブレードのモード形状ベクトル

風力タービン ブレードの動的挙動を理解することは、運転効率を最適化し、ブレードの故障を予測するために重要です。このセクションでは、風力タービン ブレードの実験的なモード解析データを解析し、ブレードのモード形状を可視化します。ハンマーでタービン ブレードの 20 個の位置を励起し、基準加速度計で位置 18 の応答を測定します。アルミニウム ブロックはブレードの土台に取り付けられ、ブレードはブレードのフラップ部分に対して垂直なフラップ方向に励起されます。各位置の FRF を収集します。FRF データはマサチューセッツ大学ローウェル校の Structural Dynamics and Acoustics Systems Laboratory により提供されました。まず、測定位置の空間的配置を可視化します。

位置 18 および 20 の風力タービン ブレードの FRF 推定を読み込んでプロットします。最初のいくつかのピークを拡大します。

[FRF,f,fs] = helperImportModalData(); helperPlotModalAnalysisExample(FRF,f,[18 20]);

最初の 2 つのモードは 37 Hz および 111 Hz 付近のピークとして現れます。安定化ダイアグラムをプロットして、固有周波数を推定します。モデル次数 14 について返される最初の 2 つの値は、周波数および減衰比が安定しています。

fn = modalsd(FRF,f,fs,'MaxModes',20); physFreq = fn(14,[1 2]);

次に、modalfitを使用して、最初の 2 つのモードのモード形状を抽出します。前のプロットに基づいて、近似を周波数範囲 0 ～ 250 Hz に制限します。

[~,~,ms] = modalfit(FRF,f,fs,14,'PhysFreq',physFreq,'FreqRange',[0 250]);

モード形状は、各位置における構造の各モードの運動振幅を定量化します。モード形状ベクトルを推定するには、周波数応答関数の行列の 1 つの行または列が必要です。実際には、これは構造 (この場合はロービング ハンマー) のすべての測定位置に励起が必要なこと、またはすべての位置の応答測定が必要なことを意味します。モード形状を可視化するには、FRF の虚数部を調べます。ブレードの片側にある位置の FRF 行列における虚数部のウォーターフォール ダイアグラムをプロットします。周波数範囲を最大 150 Hz に制限して、最初の 2 つのモードを調べます。プロットのピークがモード形状を表します。

measlocs = [3 6 9 11 13 15 17 19 20];%在叶片边缘测量位置helperPlotModalAnalysisExample(FRF,f,measlocs,150);

ピークの等高線によりプロットに示された形状は、ブレードの最初と 2 番目の曲げモーメントを表します。次に、同じ測定位置について、モード形状ベクトルの大きさをプロットします。

helperPlotModalAnalysisExample(ms,measlocs);

振幅はスケーリングが異なりますが (モード形状ベクトルはモーダル 1 の A に対してスケーリングされる)、モード形状の等高線は形状が一致します。最初のモードの形状には大きな先端の変位および 2 つのノードがあり、振動振幅はゼロです。2 番目のモードにも大きな先端の変位および 3 つのノードがあります。

まとめ

この例では、ロービング ハンマーによって励起される 3DOF システムに関してシミュレートされたモード解析データセットを解析して比較しました。LSCE と LSRF アルゴリズムおよび安定化ダイアグラムを使用して、固有振動数および減衰を推定しました。モーダル パラメーターは 2 セットの測定において一致しました。個別のケースで、FRF 行列の虚数部およびモード形状ベクトルを使用して、風力タービン ブレードのモード形状を可視化しました。

謝辞

風力タービン ブレードの実験データをご提供いただいたマサチューセッツ大学ローウェル校 Structural Dynamics and Acoustics Systems Laboratory の Peter Avitabile 博士に謝意を申し上げます。

参考文献

[1] Brandt, Anders.Noise and Vibration Analysis: Signal Analysis and Experimental Procedures. Chichester, UK: John Wiley and Sons, 2011.

[2] Vold, Håvard, John Crowley, and G. Thomas Rocklin."New Ways of Estimating Frequency Response Functions."Sound and Vibration.Vol. 18, November 1984, pp. 34-38.

[3] Peter Van Overschee and Bart De Moor."N4SID: Subspace Algorithms for the Identification of Combined Deterministic-Stochastic Systems."Automatica.Vol. 30, January 1994, pp. 75-93.

[4] Ozdemir, A. A., and S. Gumussoy."Transfer Function Estimation in System Identification Toolbox via Vector Fitting."Proceedings of the 20th World Congress of the International Federation of Automatic Control.Toulouse, France, July 2017.