# 仿真系统和风力涡轮叶片的模态分析

### 仿真波束的固有频率和阻尼

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

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

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

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

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

```fn = modalfit(FRF1,f1,fs,3,'FitMethod','PP','FreqRange',[200 1600]) ```
```fn = 1.0e+03 * 0.3727 0.8525 1.3707 ```

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

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

LSCE 算法生成的是计算模态，这些模态在结构中并不实际存在。可以使用稳定图通过检查模态数增加时极点的稳定性来识别物理模态。物理模态的固有频率和阻尼比倾向于保持在同一位置，即“保持稳定”。创建一个稳定图，并输出频率稳定的极点的固有频率。

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

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

```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,f] = modalfrf(X2,Y2,fs,winLen,'Sensor','dis','Measurement','rovinginput'); fn = modalsd(FRF,f,fs,'MaxModes',20,'SCriteria',[0.001 0.025]); ```

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

```[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,f] = modalfrf(X1,Y1,fs,winLen,'Sensor','dis','Measurement','rovinginput'); fn = modalsd(FRF,f,fs,'MaxModes',20, 'FitMethod', 'lsrf'); % Identify physical modes using lsfr physFreq = 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 ```

```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 validation ValidationData = iddata(Y0(1001:2000), X0(1001:2000), 1/fs); ```

```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.041 -1765 149.8 -1880 -49.64 -358 x2 1764 -0.3434 2197 -232.5 -438.3 -128.4 x3 -152.4 -2198 2.839 4715 255.9 547.5 x4 1880 228.2 -4714 -15.91 -1216 -28.79 x5 59.42 440.9 -275.5 1217 35.03 -8508 x6 363.7 120.2 -545.4 -44.02 8508 -92.47 B = u1 x1 0.2777 x2 -0.6085 x3 0.07123 x4 -3.658 x5 0.04771 x6 -7.642 C = x1 x2 x3 x4 x5 x6 y1 4.46e-05 -5.315e-06 -7.46e-06 -1.641e-05 2.964e-06 -5.871e-06 D = u1 y1 7.997e-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.26% (prediction focus) FPE: 1.355e-16, MSE: 1.304e-16 ```

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

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

### 风力涡轮叶片的模态振形向量

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

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

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

```measlocs = [3 6 9 11 13 15 17 19 20]; % Measurement locations on blade edge helperPlotModalAnalysisExample(FRF,f,measlocs,150); ```

```helperPlotModalAnalysisExample(ms,measlocs); ```

### 参考资料

[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.