Main Content

시뮬레이션된 시스템과 풍력 터빈 블레이드의 모드 해석

이 예제에서는 실험 데이터에서 주파수 응답 함수(FRF)와 모드 파라미터를 추정하는 방법을 보여줍니다. 첫 번째 섹션에서는 일련의 해머 충격으로 3-자유도(3DOF) 시스템에 진동을 가하고 그 결과로 나타나는 변위를 기록하는 시뮬레이션된 실험을 설명합니다. 주파수 응답 함수, 고유 주파수, 감쇠비, 모드 형상 벡터가 구조물의 세 가지 모드에 대해 추정됩니다. 두 번째 섹션에서는 풍력 터빈 블레이드 실험에서 얻은 주파수 응답 함수 추정값에서 모드 형상 벡터를 추정합니다. 터빈 블레이드 측정값 구성과 그에 따른 모드 형상이 시각화됩니다. 이 예제를 실행하려면 System Identification Toolbox(TM)가 필요합니다.

시뮬레이션된 빔에 대한 고유 주파수와 감쇠

SISO(단일 입력 단일 출력) 해머 가진(Excitation)

일련의 해머 타격으로 3자유도 시스템에 진동을 가하면 센서가 결과로 생성되는 변위를 기록합니다. 이 시스템은 비례적으로 감쇠되며, 감쇠 행렬은 질량 행렬과 경직성(Stiff) 행렬의 선형 결합입니다.

가진(Excitation) 신호, 응답 신호, 시간 신호, ground truth 주파수-응답 함수를 포함하여 두 측정값 세트에 대한 데이터를 가져옵니다. 첫 번째 응답 신호 세트 Y1은 첫 번째 질량의 변위를 측정하고, 두 번째 응답 신호 세트 Y2는 두 번째 질량의 변위를 측정합니다. 각각의 가진 신호는 10개의 이어지는 해머 충격으로 구성되며, 각각의 응답 신호는 이에 대응되는 변위를 포함합니다. 각 충격 신호의 지속 시간은 2.53초입니다. 가진 신호와 응답 신호에 가산성 잡음이 존재합니다. 첫 번째 가진(Excitation)과 첫 번째 측정값의 응답 채널을 시각화합니다.

[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 samples
modalfrf(X0,Y0,fs,winLen,'Sensor','dis')

디폴트 'H1' 추정량을 사용하여 추정된 FRF는 측정된 주파수 대역에서 세 개의 두드러진 피크를 포함합니다. 이 피크는 세 가지 유연한 진동 모드에 해당합니다. 코히어런스는 이러한 피크 근처에서는 1에 가까우며, 응답 측정값의 신호 대 잡음비가 낮은 반공진 영역에서는 낮습니다. 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);

여기서 FRF3은 피드스루 항을 포함하고 범위 1:10의 최적 차수를 갖는 상태공간 모델을 피팅하여 추정됩니다. 'H1', 'H2', 'subspace' 방법을 사용하여, 추정된 FRF와 이론적 FRF를 비교합니다.

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

모든 추정량이 응답 피크 근처에서 상당히 유사한 결과를 나타내지만, 'H2' 추정량은 반공진 영역에서 응답을 과대 추정합니다. 어떤 추정량을 선택해도 코히어런스는 영향을 받지 않습니다.

다음으로, 피크 추출 알고리즘(Peak-Picking Algorithm)을 사용하여 각 모드의 고유 주파수를 추정합니다. 피크 추출 알고리즘은 FRF에서 피크를 식별하는 데 사용할 수 있는 간단하면서도 빠른 절차입니다. 각 추정값이 단일 주파수-응답 함수에서 생성되므로 이 알고리즘은 국소적 방법입니다. 또한, 각 모드의 피크가 서로 독립인 것으로 간주되기 때문에 단일 자유도(SDOF) 방법입니다. 결과적으로, 각 FRF에 대해 모드 파라미터 세트가 생성됩니다. 이전 플롯에 기반하여 주파수 범위를 200Hz ~ 1600Hz로 지정합니다. 이 범위에는 세 개의 피크가 있습니다.

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

   1.0e+03 *

    0.3727
    0.8525
    1.3707

고유 주파수는 대략적으로 373Hz, 853Hz, 1371Hz입니다. 재생성된 FRF를 플로팅하고 modalfit을 사용하여 이 FRF를 측정된 데이터와 비교합니다. FRF는 주파수-응답 함수 행렬 FRF1에서 추정된 모드 파라미터를 사용하여 재생성됩니다. 출력 인수 없이 modalfit를 다시 호출하여 재구성된 FRF를 포함하는 플롯을 생성합니다.

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

재구성된 FRF는 첫 번째 가진과 응답 채널에 대해 측정된 FRF와 일치합니다. 다음 섹션에서는 두 개의 추가 가진 위치가 고려됩니다.

이동 해머 가진

디폴트 'H1' 추정량을 사용하여 세 개 센서 모두의 응답에 대해 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% 미만으로 변동하는 경우 감쇠에서도 안정적인 것으로 분류됩니다. 두 기준 모두 다른 값으로 조정할 수 있습니다. 안정적인 극점의 위치를 기반으로 373Hz, 852.5Hz, 1371Hz의 고유 주파수를 선택합니다. 이러한 주파수는 주파수 상에서 안정적인 다른 극점의 고유 주파수와 함께 modalsd의 출력값 fn에 포함되어 있습니다. 물리적으로 존재하는 모드 개수보다 더 높은 모델 차수는 일반적으로 LSCE 알고리즘을 사용하여 양호한 모드 파라미터 추정값을 생성하는 데 필요합니다. 이 경우, 네 개 모드의 모델 차수는 세 개의 안정적인 극점을 나타냅니다. 원하는 주파수는 fn의 네 번째 행에 있는 처음 세 개 열에 나타납니다.

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

고유 주파수와 감쇠를 추정하고, 재구성된 후 측정한 FRF를 플로팅합니다. 안정성 다이어그램 'PhysFreq'에서 확인한 네 가지 모드와 물리적 주파수를 지정합니다. 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를 계산하고 다른 위치의 센서를 사용하여 두 번째 해머 충격 세트에 대한 안정화 다이어그램을 플로팅합니다. 안정성 기준을 주파수의 경우 0.1%로 변경하고 감쇠의 경우 2.5%로 변경합니다.

[FRF,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(최소제곱 유리 함수) 접근 방식은 공유된 분모 전달 함수를 MIMO(다중 입력, 다중 출력) 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 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  

마지막으로, 모수적 방법인 FRF 추정 방법('subspace')과 모드 파라미터 추정 방법('lsrf')은 System Identification Toolbox에서 동적 모델을 시간-영역 신호나 주파수-응답 함수에 피팅할 때 사용되는 방법과 비슷합니다. 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 validation
ValidationData = 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.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

모델 sys를 사용하여 모드 파라미터를 계산합니다.

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

풍력 터빈 블레이드의 모드 형상 벡터(Mode Shape Vector)

풍력 터빈 블레이드의 역학적 동작을 이해하는 것은 운영 효율성을 최적화하고 블레이드 고장을 예측하는 데 중요합니다. 이 섹션에서는 풍력 터빈 플레이드에 대한 실험적 모드 해석 데이터를 분석하고 블레이드의 모드 형상을 시각화합니다. 해머로 20개 위치에서 터빈 블레이드에 진동을 가하고, 기준 가속도계로 위치 18에서 응답을 측정합니다. 알루미늄 블록이 블레이드 베이스에 장착되어 있으며 블레이드에는 블레이드의 편평한 부분에 수직인 플랩 방향으로 진동이 가해집니다. 각 위치에서 FRF가 수집됩니다. FRF 데이터는 메사추세츠 대학교 로웰 캠퍼스의 Structural Dynamics and Acoustics Systems Laboratory에서 제공했습니다. 먼저, 측정 위치의 공간적 배치 상태를 시각화합니다.

위치 18 및 위치 20에 대한 풍력 터빈 블레이드의 FRF 추정값을 불러오고 플로팅합니다. 처음 몇 개의 피크를 확대합니다.

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

처음 두 개 모드는 37Hz와 111Hz 근방에서 피크로 나타납니다. 고유 주파수를 추정하기 위해 안정화 다이어그램을 플로팅합니다. 모델 차수 14에 대해 반환된 처음 두 개 값은 주파수와 감쇠비에서 안정적입니다.

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

다음으로, modalfit을 사용하여 처음 두 개 모드에 대한 모드 형상을 추출합니다. 이전 플롯을 기반으로 하여 주파수 피팅 범위를 0Hz ~ 250Hz로 제한합니다.

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

모드 형상은 각 위치에서 구조물의 각 모드에 대한 운동 진폭을 수치화합니다. 모드 형상 벡터를 추정하려면 주파수-응답 함수 행렬의 한 행 또는 한 열이 필요합니다. 실제로, 이는 구조물의 모든 측정 위치에서 가진이 필요하거나(이 경우 이동 해머) 모든 위치에서 응답 측정이 필요하다는 것을 의미합니다. 모드 형상은 FRF의 허수부를 검토하여 시각화할 수 있습니다. 블레이드 한 면의 위치에 대해 FRF 행렬의 허수부에 대한 폭포 다이어그램(Waterfall Diagram)을 플로팅합니다. 주파수 범위를 최대 150Hz로 제한하여 처음 두 개 모드를 검토합니다. 플롯의 피크는 모드 형상을 나타냅니다.

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

플롯에서 피크의 등고선으로 표시되는 형상은 블레이드의 1차 및 2차 굽힘 모멘트를 나타냅니다. 다음으로, 동일한 측정 위치에 대한 모드 형상 벡터의 크기를 플로팅합니다.

helperPlotModalAnalysisExample(ms,measlocs);

진폭이 다르게 스케일링되었지만(모드 형상 벡터는 단위 모드 A로 스케일링됨), 모드 형상 등고선들의 모양이 일치합니다. 첫 번째 모드 형상에는 큰 선단 변위와 두 개의 노드(진동 진폭이 0임)가 있습니다. 두 번째 모드에도 큰 선단 변위와 세 개의 노드가 있습니다.

요약

이 예제에서는 이동 해머로 진동을 가하는 3자유도(3DOF) 시스템에 대해 시뮬레이션된 모드 해석 데이터셋을 분석하고 비교했습니다. 또한 안정화 다이어그램과 LSCE/LSRF 알고리즘을 사용하여 고유 주파수와 감쇠를 추정했습니다. 모드 파라미터는 두 개의 측정값 세트에서 일관적이었습니다. 개별 사용 사례에서, 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.

참고 항목

| |