模態分析主要研究頻率域內系統動態特性。

通過模態分析方法搞清楚了結構物在某一易受影響的頻率范圍內的各階主要模態的特性，就可以預言結構在此頻段內在外部或內部各種振源作用下產生的實際振動響應。

模態分析主要分為解析模態分析和試驗模態分析

解析模態分析

通常我們針對一個線性定常系統進行動力學描述可以得到方程組：

其中，[M] 是質量矩陣，[C] 是阻尼矩陣，[K] 是剛度矩陣，{x(t)} 是位移向量, {F(t)} 是力矩陣。 我們的目標是求解這個線性定常系統振動微分方程組得到 {x(t)}，也就是系統上各點隨時間的位移。 對于單自由度的系統是方便求解的，所以對于這種多自由度系統我們首先想到的是通過坐標變換，用一組新的正交基（模態空間里的基）去描述 {x(t)}，例如 [P?1]x(t)，來實現對方程組 (1) 解耦，從而將問題轉化為互相獨立的子方程（組），用更少自由度甚至單自由度的方程進行求解。 其中[P?1]就是模態矩陣，其每列為模態振型，它描述的是在新的解耦后的坐標系中的各維坐標與原來坐標系中各個維度的線性關系。

例如對于一個簡化的多自由度的彈簧系統，這個線性定常系統由 3 個相同重量的模塊組成，m?=m?=m?=m，他們中間用彈簧鏈接, 為了簡化問題，我們設彈簧的勁度系數 k?=k?=k?=k?=k，阻尼系數 c?=c?=c?=c?=0。 定義 x?,x?,x??為每個質量塊的位移，另外 k? ,k?邊緣兩端的位移。由于兩端固定，都為0。每個質量塊的運動方程分別為：

將上述方程寫為矩陣形式，上述運動學方程組可以簡化為：

其中

對于方程組 (2)，如何進行坐標解耦呢？ 計算時運動學方程組(2) 右側項可以忽略, 只保留質量矩陣項 [M] 和剛度矩陣 [K] 項，即

對于方程組(3) 通常的做法是轉換為：

本質對方程組(4) 解耦實際上就是求解質量矩陣關于剛度矩陣的廣義特征值的問題。也就是計算得到特征值，

和特征向量，

使得[M][P]=[K][P][D] (5) 其中特征向量 [P] 即為模態向量（空間基向量），為對應的特征值對角陣對應解耦后固有頻率的平方，即

具體此處不做推導，但可以簡單的從必要性上進行解釋： 假設已經通過空間變換矩陣得到新的坐標

我們計算一下特征向量矩陣是否將原始方程組坐標由 {X} 變換為后 {Q} 得到單自由度的二階常微分方程組。 我們將{X}=[P]{Q} 帶入方程(3)得到

同時我們將(5) 帶入(6) 可以得到

在 [K][P] 均可逆的條件下，我們得到方程

即：

也就是完全解耦的單自由度的二階常微分方程，接下來可以單獨求解 q?(t), q?(t), q?(t), 最后只需要再做一次 [P] 變換將模態空間坐標變換到物理坐標即可。

?????????

% 'M' 是質量矩陣,'K' 是剛度矩陣. 'm' 質量塊質量,單位 Kgs

% 'k' 剛度系數,單位N/m. 'P' 和'D' 是特征向量和特征值.

m = 0.003;

k = 180000;

M = m*eye(3);

K = k*[2 -1 0 ;

-1 2 -1 ;

0 -1 2 ];

% P為特征向量：變換矩陣，D為特征值：固有頻率的平方，w_nat 包含自然響應頻率.

[P,D] = eig(K,M)

w_nat = sqrt(D)

% 我們查看低階模態振型，也就是低頻振型，可以嘗試設置number

% 查看三階模態振型'EIGS' 函數可以用來計算特征值和特征向量

number = 3;

[P,D]=eigs(K,M,number,'smallestabs');

% 因為系統兩端固定，模態振型坐標在這兩端為0

vect_aug1 = [0 0 0;P;0 0 0];

c = ['m','b','r'];

figure(1)

for i=1:size(vect_aug1,2)

plot(vect_aug1(:,i),c(i))

hold on;

end

最終

和

的求解以及繪制都可以用通過 MATLAB 腳本實現。大家可以自己嘗試。 當然對于我們的系統不可能是這種簡單的系統，通常要結合有限元的手段進行計算。 MATLAB 中的 Partial Differential EquationToolbox 對于滿足我們一些基礎的需求可以提供求解方案。 我們看一個基于 MATLAB 有限元計算模態的示例。 本示例是對 KinovaGen3 機械臂肩部關節部件進行模態和頻率響應分析。機械臂通過多個關節鏈接，一端固定。這些鏈接結構強度要夠大以避免電機帶動負載運動時產生振動。 機械臂終端的負載會讓每個鏈接處產生壓力。壓力的方向取決于負載的方向。此示例主要展示如何分析 Kinova Gen3 超輕量級機械臂的肩部關節連接部件在一定壓力下可能的形變。

模態分析MATLAB 中有限元求解流程

model = createpde('structural','modal-solid');

importGeometry(model,'Gen3Shoulder.stl');

generateMesh(model);

structuralProperties(model,'YoungsModulus',E, ...

'PoissonsRatio',nu, ...

'MassDensity',rho);

將 facelabel 可視化出來方便設置邊界約束和負載。

face3 是固定的，face4 是活動的。設置約束

structuralBC(model,'Face',3,'Constraint','fixed');

在關心的頻率范圍進行模型求解。

RF = solve(model,'FrequencyRange',[-1,10000]*2*pi);

modeID = 1:numel(RF.NaturalFrequencies);

通過對結果除以2pi,得到Hz單位的頻率值：

tmodalResults = table(modeID.',RF.NaturalFrequencies/2/pi);

tmodalResults.Properties.VariableNames = {'Mode','Frequency'};

disp(tmodalResults);

同樣我們需要可視化模態振型。最好的方式是以各階模態頻率的簡諧振動動畫顯示出來，此處顯示前六階模態振型：

頻率響應分析模擬機械臂關節在壓力載荷下的動力學，假設附加其上的連桿對各半面分別施加大小相等方向相反的壓力。分析面上某點的頻率響應和形變。

同樣跟上述流程一樣，創建結構，導入幾何形狀，生成網格。

其他過程跟模態分析相同，區別在于加一個力，使用 pressFcnFR 函數在面 4 上施加邊界載荷。 這個函數作用一個推力和一個扭轉壓力信號。推壓分量是均勻的。扭力對左側面施加正壓力，對右側施加負壓力。

structuralBoundaryLoad(fmodel,'Face',4,'Pressure',@(region,state)pressFcnFR(region, state),'Vectorized','on');

這個壓力函數跟頻率無關，作用于所有頻率。

同樣進行求解

R = solve(fmodel,flist,'ModalResults',RF)

我們可以看面 4對應的負向負荷最大的點（0.003; 0.0436; 0.1307）對應的位移

queryPoint= [0.003; 0.0436; 0.1307];

queryPointDisp =R.interpolateDisplacement(queryPoint);

響應的峰值出現在 2662Hz 附近，與二階模態接近。在接近 1947Hz 的一階模態上也會出現小的響應。

通過使用 max 函數找到峰值響應頻率對應的峰值和索引。

[M, I] = max(abs(queryPointDisp.uy))

M = 1.1256e-04

I = 152 繪制峰值響應頻率處的變形。

可以看到所施加的載荷主要激發了部件的開口模態和彎曲模態。 通過上面兩個示例，我們針對計算模態的場景可以在 MATLAB 中通過相應的內置函數做探索研究。

試驗模態分析

試驗模態分析主要是通過實測實驗數據進行頻率響應估計并計算相應模態參數的一種方法。

我們通過 MATLAB 自帶的一個例子來介紹試驗模態分析。

詳見：https://ww2.mathworks.cn/help/releases/R2021a/ident/ug/modal-analysis-of-a-flexible-flying-wing-aircraft.html

這是明尼蘇達大學無人飛行器實驗室的小型柔性飛翼飛機的示例。這個例子展示了柔性機翼飛機彎曲模態的計算過程。

通過沿翼型進行機翼的振動響應的多點采集獲得數據，用于辨識系統的動態模型。

從辨識出的模型中提取模態參數。

將模態參數數據與傳感器位置信息相結合，可視化機翼的各種彎曲模態。

實驗設置

實驗的目的是收集飛機在外部激勵下不同位置的振動響應。

這架飛機懸掛在一個木制框架上，其重心由一根彈簧支撐。該彈簧具有足夠的彈性保證彈簧-質量振蕩的固有頻率不會干擾飛機的基頻。通過一個電動激振器在飛機中心附近施加輸入力。

沿著翼展放置 20 個加速度計來收集振動響應，如下圖所示

激振器輸入命令指定為一個恒定振幅的 chirp 信號 Asin(ω(t)t), chirp 信號的頻率隨時間線性增加，ω(t)=c0+c1t, 頻率范圍為 3 - 35hz。試驗數據由兩個加速度計(在 x 方向的前緣和后緣)一次收集。總共進行了 10 組實驗來收集所有 20 個加速度計的響應。加速度計和力傳感器的測量頻率都是 2000hz。

數據準備

數據由 10 組單輸入/雙輸出信號表示，每組包含 600K 個樣本。

變量 MeasuredData 是一個結構體。結構體的每個域還是一個結構體，包含兩個響應 y 和對應輸入 u。隨機繪制第一次實驗的數據。

fs = 2000; % data sampling frequency

Ts = 1/fs; % sample time

y = MeasuredData.Exp1.y; % 加速度計的輸出值，每個包含兩列

u = MeasuredData.Exp1.u; % input force data

t = (0:length(u)-1)' * Ts;

為了用于系統辨識，將數據封裝到 iddata 對象中，并將 10 次試驗合并。

Data = merge(Data{:})

系統辨識

我們想識別一個頻率響應與實際飛機盡可能接近的動態模型。

動態模型將系統的輸入和輸出之間的數學關系轉化為微分方程或差分方程。例如傳遞函數和狀態空間模型都是動態模型。

動態模型可以通過在時域或頻域對試驗測量數據運行 fest 和 sest 之類的模型估計命令來創建。

這個例子中，我們先用 etfe 命令進行傳遞函數估計，從而將測量數據從時域轉換為頻率響應。然后利用估計的頻響函數用于辨識飛機振動響應的狀態空間模型。

當然直接利用時域數據進行狀態空間模型辨識也是可以的。但頻響函數的形式可以將大數據集壓縮成更少的樣本，并且更方便的在相關的頻率范圍調整估計過程。

針對每組實驗數據（兩輸出/單輸入）進行頻響函數(FRF)估算。使用 24000 個頻率點進行無窗響應計算。

G = cell(1, 10);

N = 24000;

for k = 1:10

% Convert time-domain data into a FRF using ETFE command

Data_k = getexp(Data, k);

G{k} = etfe(Data_k, [], N); % G{k} is an @idfrd object

end

G = cat(1,G{:}); % 將所有的頻響函數合并為一個（20輸出/一個輸入）的頻響

G.OutputName = 'y'; % name outputs 'y1', 'y2', ..., 'y20'

G.InterSample = 'bl';

我們隨便挑選第 1 個和第 15 個輸出幅值進行繪制來看一下頻率響應函數的估計結果。我們關注的頻率范圍(4 - 35hz)。

頻響函數顯示至少 9 個諧振頻率(峰值)。我們最關心飛機的機翼彎曲模態，這些模態主要集中在 6- 35hz 的頻率范圍，因此我們只選擇這個范圍的頻響。

f =G.Frequency/2/pi; % 單位變換

Gs = fselect(G,f>6 & f<=32)??? % 選擇頻響頻率范圍 (6.5 - 35 Hz)

接下來，計算一個狀態空間模型來逼近 Gs。子空間估計器 n4sid 提供了一個快速的非迭代估計。狀態空間模型參數配置會影響精度:

1. 模型階數的選擇。通常要多次嘗試。

2. 模型結構。例如是否包含饋通項(狀態空間模型中的 D 矩陣是否為零)，等等。

3. 設置權重項，確保低振幅和高振幅對結果有相同的影響。

FRF =squeeze(Gs.ResponseData);

Weighting = mean(1./sqrt(abs(FRF))).';

n4Opt =n4sidOptions('EstimateCovariance',false,...

'WeightingFilter',Weighting,...

'OutputWeight',eye(20));

sys1 = n4sid(Gs,24,'Ts',0,'Feedthrough',true,n4Opt);

Fit = sys1.Report.Fit.FitPercent'

通過查看 Fit 的效果，顯示這 20 個響應中最好的和最差的

可以看到 sys1 結果仍然有待提高。通過對模型參數進行非線性最小二乘優化迭代，可以進一步改善模型的擬合效果。這可以使用 sest 命令來實現。這一次也估計了參數協方差。

ssOpt = ssestOptions('EstimateCovariance',true,...

'WeightingFilter',n4Opt.WeightingFilter,...

'SearchMethod','lm',... % use Levenberg-Marquardt search method

'Display','on',...

'OutputWeight',eye(20));

sys2 = ssest(Gs, sys1, ssOpt); % estimate state-space model (this takesseveral minutes)

Fit = sys2.Report.Fit.FitPercent'

我們同樣看看擬合效果：最好的和最差的幅值進行可視化。同時將 1-sd 置信區間可視化出來。

優化后的狀態空間模型 sys2 在 7 - 20hz 區域的頻響擬合效果很好。多數共振位置的不確定性區間都很窄。我們最開始設置的是一個 24 階模型，這意味著在系統 sys2 中最多有 12 個諧振模態。使用 modalfit 命令獲取這些模態的固有頻率。

f = Gs.Frequency/2/pi; % data frequencies (Hz)

fn = modalfit(sys2, f, 12); % naturalfrequencies (Hz)

fn 中的值表明兩個頻率非常接近 7.8 Hz，三個接近 9.4 Hz。查看這些頻率附近的各點頻率響應，峰值位置在輸出中確實發生了一些偏移。

這些差異可以通過更好地控制實驗過程，直接利用頻率為中心的狹窄范圍內的時域數據進行直接辨識，或將這些頻率附近的頻率響應擬合為單個振動模態。本例中不討論這些替代方案。

模態參數計算

現在我們可以使用模型 sys2 來提取模態參數。通過查看頻響結果找到 10 個模態頻率，大約在頻率 [5 7 10 15 17 23 27 30]Hz 附近左右。通過使用 modalsd 函數可讓估計更加準確，該函數嘗試不同模型的階數來檢查模態參數的穩定性。

通過穩定圖可以看到一組更好的自然響應頻率

Freq = [7.8 9.4 15.3 19.3 23.0 27.3 29.231.7];

我們使用 Freq 向量的值作為從模型 sys2 中選擇主要模態的基準。使用 modalfit 函數實現

[fn,dr,ms] = modalfit(sys2,f,length(Freq),'PhysFreq',Freq);

fn 是固有頻率 (Hz)， dr 是相應的阻尼系數，ms 是模態振型向量 (每個固有頻率一列)。

模態振型可視化

我們使用上述模態參數可視化各種彎曲模態。此外，我們需要各測量點位置的信息。

模態振型包含在矩陣 ms 中，其中每一列對應于一個頻率的振型。通過在加速度傳感器位置坐標上疊加模態振型的振幅并以模態固有頻率改變振幅來做動畫展示。

結論

這個例子展示了一種基于參數模型辨識的模態參數估計方法。利用 24 階狀態空間模型，在 6 ~ 32 Hz 頻率范圍內提取了 8 個穩定的振動模態。將模態信息與位置信息相結合，可視化各種彎曲模態。

當然，您也可以對其他設備例如風力發電機的葉片振型進行計算，了解風力機葉片的動態特性從而優化運行效率和預測葉片失效，可以按同樣的思路實現。