一、振动微分方程的建立

无阻尼多自由度线性系统的运动方程可由达朗贝尔原理、影响系数方法或拉格朗日方程建立。对于保守系统，也可利用机械能守恒关系。

考虑不计水平方向运动、无外载荷作用的刚硬翼段。线弹簧刚度为 $K_h$，扭转弹簧刚度为 $K_\alpha$；翼段质量为 $m$，绕质心的转动惯量为 $J_0$。取弯心偏离平衡位置的位移 $h$ 和翼段绕弯心的转角 $\alpha$ 为广义坐标，正方向如图所示。质心与弯心的距离为 $b$，在线性小转角近似下，质心位移为 $h+b\alpha$。

1、达朗贝尔原理

质心加速度为 $\ddot h+b\ddot\alpha$。沿位移正方向列平动方程，有

$$ -K_hh=m(\ddot h+b\ddot\alpha) $$

绕弯心列达朗贝尔力矩平衡方程，平动惯性力对弯心的力矩为 $-mb(\ddot h+b\ddot\alpha)$，绕质心的惯性力偶矩为 $-J_0\ddot\alpha$，故

$$ -K_\alpha\alpha-mb(\ddot h+b\ddot\alpha)-J_0\ddot\alpha=0 $$

整理得

$$ \begin{cases} m\ddot h+mb\ddot\alpha+K_hh=0\\ mb\ddot h+(J_0+mb^2)\ddot\alpha+K_\alpha\alpha=0 \end{cases} $$

记质量对弯心的静矩为 $S_\alpha=mb$，绕弯心的转动惯量为 $J_\alpha=J_0+mb^2$，则

$$ \begin{cases} m\ddot h+S_\alpha\ddot\alpha+K_hh=0\\ S_\alpha\ddot h+J_\alpha\ddot\alpha+K_\alpha\alpha=0 \end{cases} $$

定义位移列向量、质量矩阵和刚度矩阵

$$ \boldsymbol{x}=\begin{bmatrix}h\\\alpha\end{bmatrix},\qquad \boldsymbol{M}=\begin{bmatrix}m&S_\alpha\\S_\alpha&J_\alpha\end{bmatrix},\qquad \boldsymbol{K}=\begin{bmatrix}K_h&0\\0&K_\alpha\end{bmatrix} $$

运动方程的矩阵形式为

$$ \boldsymbol{M}\ddot{\boldsymbol{x}}+\boldsymbol{K}\boldsymbol{x}=\boldsymbol{0} $$

质量矩阵也称惯性矩阵，其元素 $m_{ij}$ 为质量影响系数：仅令坐标 $j$ 产生单位加速度，其余坐标的加速度为零时，在坐标 $i$ 上所需施加的广义力。刚度矩阵也称弹性矩阵，其元素 $k_{ij}$ 为刚度影响系数：仅令坐标 $j$ 产生单位位移，其余坐标的位移为零时，在坐标 $i$ 上所需施加的广义力。坐标排列一旦确定，矩阵行列的排列也随之确定。

对于线性保守系统，$\boldsymbol{M}$ 和 $\boldsymbol{K}$ 为实对称矩阵，即 $m_{ij}=m_{ji}$、$k_{ij}=k_{ji}$。由二次型动能和势能可写出

$$ m_{ij}=\frac{\partial^2T}{\partial\dot x_i\partial\dot x_j},\qquad k_{ij}=\frac{\partial^2V}{\partial x_i\partial x_j} $$

二阶混合偏导数的对称性给出了影响系数的互等关系。

矩阵的对角元素称为主项，非对角元素称为耦合项。当前坐标下，$S_\alpha\ne0$ 时存在惯性耦合；刚度矩阵的非对角元素为零，不存在弹性耦合。

若改取质心位移为平动坐标，并仍记为 $h$，则弯心位移为 $h-b\alpha$。此时

$$ T=\frac12m\dot h^2+\frac12J_0\dot\alpha^2,\qquad V=\frac12K_h(h-b\alpha)^2+\frac12K_\alpha\alpha^2 $$

由此得到

$$ \begin{cases} m\ddot h+K_h(h-b\alpha)=0\\ J_0\ddot\alpha+K_\alpha\alpha+K_h(b\alpha-h)b=0 \end{cases} $$

即

$$ \begin{bmatrix}m&0\\0&J_0\end{bmatrix} \begin{bmatrix}\ddot h\\\ddot\alpha\end{bmatrix} +\begin{bmatrix}K_h&-K_hb\\-K_hb&K_\alpha+K_hb^2\end{bmatrix} \begin{bmatrix}h\\\alpha\end{bmatrix} =\begin{bmatrix}0\\0\end{bmatrix} $$

此时只有弹性耦合。可见，惯性耦合和弹性耦合取决于坐标选择。对质量矩阵正定、刚度矩阵对称的系统，可选择使二者同时对角化的主坐标，也称模态坐标。

2、影响系数方法

柔度影响系数 $F_{ij}$ 是仅在广义坐标 $j$ 上施加单位广义力时，广义坐标 $i$ 产生的静位移。对于线性保守系统，有 $F_{ij}=F_{ji}$。

令 $Q_1$、$Q_2$ 分别为与 $h$、$\alpha$ 对应的广义力，则静力关系为

$$ \begin{cases} h=F_{11}Q_1+F_{12}Q_2\\ \alpha=F_{21}Q_1+F_{22}Q_2 \end{cases} $$

弯心坐标下，两个弹簧独立提供恢复力和恢复力矩，故

$$ F_{11}=\frac1{K_h},\qquad F_{22}=\frac1{K_\alpha},\qquad F_{12}=F_{21}=0 $$

自由振动时，将惯性广义力作为瞬时等效载荷，有

$$ Q_1=-m(\ddot h+b\ddot\alpha),\qquad Q_2=-mb\ddot h-(J_0+mb^2)\ddot\alpha $$

代入静力关系得

$$ h=-\frac1{K_h}(m\ddot h+mb\ddot\alpha) $$

$$ \alpha=-\frac1{K_\alpha}\left[mb\ddot h+(J_0+mb^2)\ddot\alpha\right] $$

定义柔度矩阵

$$ \boldsymbol{F}=\begin{bmatrix}1/K_h&0\\0&1/K_\alpha\end{bmatrix} $$

则运动方程为

$$ \boldsymbol{F}\boldsymbol{M}\ddot{\boldsymbol{x}}+\boldsymbol{x}=\boldsymbol{0} $$

记动力矩阵 $\boldsymbol{D}=\boldsymbol{F}\boldsymbol{M}$，也可写为

$$ \boldsymbol{D}\ddot{\boldsymbol{x}}+\boldsymbol{x}=\boldsymbol{0} $$

$\boldsymbol{F}$、$\boldsymbol{M}$ 虽然各自对称，但

$$ \boldsymbol{D}^{\mathrm T}=\boldsymbol{M}\boldsymbol{F} $$

只有两矩阵可交换时，$\boldsymbol{D}$ 才对称。

也可利用刚度影响系数写出

$$ \begin{cases} Q_1=k_{11}h+k_{12}\alpha\\ Q_2=k_{21}h+k_{22}\alpha \end{cases} $$

由 $k_{11}=K_h$、$k_{22}=K_\alpha$、$k_{12}=k_{21}=0$，代入惯性广义力得

$$ \begin{cases} K_hh=-m(\ddot h+b\ddot\alpha)\\ K_\alpha\alpha=-mb\ddot h-(J_0+mb^2)\ddot\alpha \end{cases} $$

两种影响系数方法所得方程相同。柔度系数可通过单位载荷下的静位移确定；刚度系数则要求约束其余坐标，仅使一个坐标产生单位位移。

当 $\boldsymbol{K}$ 可逆时，由 $\boldsymbol{Q}=\boldsymbol{K}\boldsymbol{x}$ 得

$$ \boldsymbol{x}=\boldsymbol{K}^{-1}\boldsymbol{Q},\qquad \boldsymbol{F}=\boldsymbol{K}^{-1} $$

若系统存在刚体自由度，则 $\boldsymbol{K}$ 奇异，整体柔度矩阵不能由上述逆矩阵定义。刚度矩阵形式仍可用于描述系统运动。

3、拉格朗日方程方法

取 $q_1=h$、$q_2=\alpha$，系统动能与势能为

$$ \begin{aligned} T&=\frac12m(\dot h+b\dot\alpha)^2+\frac12J_0\dot\alpha^2\\ &=\frac12m\dot h^2+mb\dot h\dot\alpha+\frac12(J_0+mb^2)\dot\alpha^2 =\frac12\dot{\boldsymbol{x}}^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}} \end{aligned} $$

$$ V=\frac12K_hh^2+\frac12K_\alpha\alpha^2 =\frac12\boldsymbol{x}^{\mathrm T}\boldsymbol{K}\boldsymbol{x} $$

无外加广义力时，拉格朗日方程为

$$ \frac{\mathrm d}{\mathrm dt}\frac{\partial T}{\partial\dot q_i} -\frac{\partial T}{\partial q_i} +\frac{\partial V}{\partial q_i}=0 $$

对 $h$，有

$$ \frac{\partial T}{\partial\dot h}=m\dot h+mb\dot\alpha,\qquad \frac{\partial T}{\partial h}=0,\qquad \frac{\partial V}{\partial h}=K_hh $$

因而

$$ m\ddot h+mb\ddot\alpha+K_hh=0 $$

对 $\alpha$，有

$$ \frac{\partial T}{\partial\dot\alpha}=mb\dot h+(J_0+mb^2)\dot\alpha,\qquad \frac{\partial T}{\partial\alpha}=0,\qquad \frac{\partial V}{\partial\alpha}=K_\alpha\alpha $$

从而

$$ mb\ddot h+(J_0+mb^2)\ddot\alpha+K_\alpha\alpha=0 $$

对于自由度较多的系统，直接写出动能、势能后应用拉格朗日方程，通常比逐一分析惯性力及力矩更方便。

二、振动微分方程的本征解

考虑 $n$ 自由度无阻尼线性保守系统

$$ \boldsymbol{M}\ddot{\boldsymbol{x}}+\boldsymbol{K}\boldsymbol{x}=\boldsymbol{0},\qquad \boldsymbol{x}^{\mathrm T}=[x_1\ x_2\ \cdots\ x_n] $$

其中 $\boldsymbol{M}$ 为实对称正定矩阵，$\boldsymbol{K}$ 为实对称半正定矩阵；没有零频模态时，$\boldsymbol{K}$ 为正定矩阵。

1、广义本征方程

取简谐形式的特解

$$ \boldsymbol{x}=a\boldsymbol{\varphi}\sin(\omega t+\theta),\qquad \boldsymbol{\varphi}^{\mathrm T}=[\varphi_1\ \varphi_2\ \cdots\ \varphi_n] $$

$\boldsymbol{\varphi}$ 表示各坐标振幅的比例，$a$ 为幅值系数，$\theta$ 为初相位。求导得

$$ \ddot{\boldsymbol{x}}=-\omega^2a\boldsymbol{\varphi}\sin(\omega t+\theta) $$

代入运动方程，要求其对任意时刻成立，得到

$$ (\boldsymbol{K}-\omega^2\boldsymbol{M})\boldsymbol{\varphi}=\boldsymbol{0} $$

令 $\lambda=\omega^2$，则

$$ \boldsymbol{K}\boldsymbol{\varphi}=\lambda\boldsymbol{M}\boldsymbol{\varphi} $$

称为广义本征方程。$\lambda$ 为本征值，非零向量 $\boldsymbol{\varphi}$ 为本征向量，二者组成一个本征对。

由于 $\boldsymbol{M}$ 可逆，也可写成标准本征方程

$$ \boldsymbol{M}^{-1}\boldsymbol{K}\boldsymbol{\varphi}=\lambda\boldsymbol{\varphi} $$

当 $\boldsymbol{K}$ 可逆时，还可写成

$$ \boldsymbol{D}\boldsymbol{\varphi}=\lambda^{-1}\boldsymbol{\varphi},\qquad \boldsymbol{D}=\boldsymbol{K}^{-1}\boldsymbol{M} $$

为保持对称性，对 $\boldsymbol{M}$ 作 Cholesky 分解

$$ \boldsymbol{M}=\boldsymbol{L}\boldsymbol{L}^{\mathrm T} $$

其中 $\boldsymbol{L}$ 是对角元素为正的可逆下三角矩阵。定义

$$ \widetilde{\boldsymbol{\varphi}}=\boldsymbol{L}^{\mathrm T}\boldsymbol{\varphi},\qquad \boldsymbol{\varphi}=\boldsymbol{L}^{-\mathrm T}\widetilde{\boldsymbol{\varphi}} $$

代入广义本征方程，有

$$ \boldsymbol{K}\boldsymbol{L}^{-\mathrm T}\widetilde{\boldsymbol{\varphi}} =\lambda\boldsymbol{L}\widetilde{\boldsymbol{\varphi}} $$

左乘 $\boldsymbol{L}^{-1}$，得到

$$ \widetilde{\boldsymbol{K}}\widetilde{\boldsymbol{\varphi}} =\lambda\widetilde{\boldsymbol{\varphi}},\qquad \widetilde{\boldsymbol{K}}=\boldsymbol{L}^{-1}\boldsymbol{K}\boldsymbol{L}^{-\mathrm T} $$

因为

$$ \widetilde{\boldsymbol{K}}^{\mathrm T} =\boldsymbol{L}^{-1}\boldsymbol{K}^{\mathrm T}\boldsymbol{L}^{-\mathrm T} =\widetilde{\boldsymbol{K}} $$

所以广义本征问题化为实对称矩阵的标准本征问题。并且，对任意非零向量 $\boldsymbol{y}$，有

$$ \boldsymbol{y}^{\mathrm T}\widetilde{\boldsymbol{K}}\boldsymbol{y} =(\boldsymbol{L}^{-\mathrm T}\boldsymbol{y})^{\mathrm T} \boldsymbol{K}(\boldsymbol{L}^{-\mathrm T}\boldsymbol{y})\ge0 $$

故 $\widetilde{\boldsymbol{K}}$ 也为半正定矩阵，其本征值均为非负实数，并具有一组完备的实正交本征向量。

2、固有模态和主振动

齐次线性方程组

$$ (\boldsymbol{K}-\lambda\boldsymbol{M})\boldsymbol{\varphi}=\boldsymbol{0} $$

存在非零解的充要条件是

$$ p(\lambda)=\det(\boldsymbol{K}-\lambda\boldsymbol{M})=0 $$

这是关于 $\lambda$ 的 $n$ 次方程，称为本征方程；以 $\lambda=\omega^2$ 表示时，也称频率方程。固有频率按从小到大排列为

$$ 0\le\omega_1\le\omega_2\le\cdots\le\omega_n,\qquad \omega_j=\sqrt{\lambda_j} $$

$\omega_1$ 称为第一阶固有频率或基频。刚度矩阵的半正定性指

$$ \boldsymbol{z}^{\mathrm T}\boldsymbol{K}\boldsymbol{z}\ge0 \qquad(\boldsymbol{z}\in\mathbb{R}^n) $$

若存在零频，则 $\det\boldsymbol{K}=0$；若 $\boldsymbol{K}$ 正定，则所有固有频率均为正。由本征方程左乘 $\boldsymbol{\varphi}_j^{\mathrm T}$ 还可得到

$$ \lambda_j= \frac{\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{K}\boldsymbol{\varphi}_j} {\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_j}\ge0 $$

将 $\lambda_j$ 代回齐次方程，得到对应的固有振型 $\boldsymbol{\varphi}_j$，也称主振型或模态向量。固有频率与固有振型共同组成第 $j$ 阶固有模态。

当 $\omega_j>0$ 时，第 $j$ 阶主振动为

$$ \boldsymbol{x}_j=a_j\boldsymbol{\varphi}_j\sin(\omega_jt+\theta_j) $$

同一阶主振动中，所有非零分量具有相同频率。振型分量同号的坐标同相运动，异号的坐标反相运动，相位差分别为 $0$ 和 $\pi$；振型分量为零的坐标在该阶主振动中不动。

没有零频时，系统自由响应可表示为各阶主振动之和

$$ \begin{aligned} \boldsymbol{x} &=\sum_{j=1}^n\boldsymbol{x}_j =\sum_{j=1}^na_j\boldsymbol{\varphi}_j\sin(\omega_jt+\theta_j)\\ &=\sum_{j=1}^n\boldsymbol{\varphi}_j \left(A_j\cos\omega_jt+B_j\sin\omega_jt\right) \end{aligned} $$

其中 $A_j=a_j\sin\theta_j$、$B_j=a_j\cos\theta_j$，由初始条件

$$ \boldsymbol{x}(0)=\boldsymbol{x}_0,\qquad \dot{\boldsymbol{x}}(0)=\dot{\boldsymbol{x}}_0 $$

确定。零频对应的解为时间的一次函数，需另行计入。

3、固有振型正交性和主坐标

（1）不同固有频率对应的振型正交性

定理 设 $\boldsymbol{M}$ 实对称正定、$\boldsymbol{K}$ 实对称。若 $\omega_i^2\ne\omega_j^2$，则对应振型满足

$$ \boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_j=0,\qquad \boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{K}\boldsymbol{\varphi}_j=0 $$

证明
两个本征方程分别为

$$ \boldsymbol{K}\boldsymbol{\varphi}_j=\omega_j^2\boldsymbol{M}\boldsymbol{\varphi}_j,\qquad \boldsymbol{K}\boldsymbol{\varphi}_i=\omega_i^2\boldsymbol{M}\boldsymbol{\varphi}_i $$

第一式左乘 $\boldsymbol{\varphi}_i^{\mathrm T}$，得到

$$ \boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{K}\boldsymbol{\varphi}_j =\omega_j^2\boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_j $$

第二式左乘 $\boldsymbol{\varphi}_j^{\mathrm T}$ 后转置，利用两矩阵的对称性，得到

$$ \boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{K}\boldsymbol{\varphi}_j =\omega_i^2\boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_j $$

两式相减，有

$$ (\omega_j^2-\omega_i^2)\boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_j=0 $$

由 $\omega_i^2\ne\omega_j^2$，得质量正交关系。再代入第一式，得到刚度正交关系。该结论同样适用于零频振型与非零频振型之间。

同一重频对应的振型可以选择为质量正交的，具体构造见第四节。

定义第 $j$ 阶模态质量和模态刚度

$$ M_{pj}=\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_j,\qquad K_{pj}=\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{K}\boldsymbol{\varphi}_j $$

其中下标 $p$ 表示主坐标。由于 $\boldsymbol{M}$ 正定，有 $M_{pj}>0$；并且

$$ K_{pj}=\omega_j^2M_{pj},\qquad \omega_j=\sqrt{\frac{K_{pj}}{M_{pj}}} $$

非零频模态有 $K_{pj}>0$，零频模态有 $K_{pj}=0$。

（2）振型矩阵与坐标解耦

将已选为质量正交的 $n$ 个固有振型按列排列，得到振型矩阵

$$ \boldsymbol{\Phi}= \begin{bmatrix}\boldsymbol{\varphi}_1&\boldsymbol{\varphi}_2&\cdots&\boldsymbol{\varphi}_n\end{bmatrix} =\begin{bmatrix} \varphi_{11}&\varphi_{12}&\cdots&\varphi_{1n}\\ \varphi_{21}&\varphi_{22}&\cdots&\varphi_{2n}\\ \vdots&\vdots&&\vdots\\ \varphi_{n1}&\varphi_{n2}&\cdots&\varphi_{nn} \end{bmatrix} $$

$\varphi_{ij}$ 表示第 $j$ 阶振型在第 $i$ 个物理坐标上的分量。若

$$ \sum_{j=1}^nc_j\boldsymbol{\varphi}_j=\boldsymbol{0} $$

左乘 $\boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{M}$，利用正交性得

$$ c_iM_{pi}=0 $$

由 $M_{pi}>0$ 知所有 $c_i=0$。因此振型线性无关，$\boldsymbol{\Phi}$ 可逆，构成 $\mathbb{R}^n$ 的一组基。

任意位移可唯一展开为

$$ \boldsymbol{x}=\sum_{j=1}^n\boldsymbol{\varphi}_jq_j=\boldsymbol{\Phi}\boldsymbol{q} $$

其中 $\boldsymbol{q}^{\mathrm T}=[q_1\ q_2\ \cdots\ q_n]$ 为主坐标列向量，$q_j$ 也称模态坐标。代入运动方程并左乘 $\boldsymbol{\Phi}^{\mathrm T}$，得

$$ \boldsymbol{M}_p\ddot{\boldsymbol{q}}+\boldsymbol{K}_p\boldsymbol{q}=\boldsymbol{0} $$

其中

$$ \boldsymbol{M}_p=\boldsymbol{\Phi}^{\mathrm T}\boldsymbol{M}\boldsymbol{\Phi} =\operatorname{diag}(M_{p1},\ldots,M_{pn}) $$

$$ \boldsymbol{K}_p=\boldsymbol{\Phi}^{\mathrm T}\boldsymbol{K}\boldsymbol{\Phi} =\operatorname{diag}(K_{p1},\ldots,K_{pn}) $$

二者分别称为模态质量矩阵和模态刚度矩阵，也称主质量矩阵和主刚度矩阵。原来的耦合方程化为 $n$ 个相互独立的方程

$$ \ddot q_j+\omega_j^2q_j=0,\qquad j=1,2,\ldots,n $$

当 $\omega_j>0$ 时，其解为

$$ q_j=A_j\cos\omega_jt+B_j\sin\omega_jt $$

（3）正交性的能量意义

将 $\boldsymbol{x}=\boldsymbol{\Phi}\boldsymbol{q}$ 代入动能和势能，有

$$ \begin{aligned} T&=\frac12\dot{\boldsymbol{x}}^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}} =\frac12\dot{\boldsymbol{q}}^{\mathrm T}\boldsymbol{M}_p\dot{\boldsymbol{q}} =\frac12\sum_{j=1}^nM_{pj}\dot q_j^2\\ V&=\frac12\boldsymbol{x}^{\mathrm T}\boldsymbol{K}\boldsymbol{x} =\frac12\boldsymbol{q}^{\mathrm T}\boldsymbol{K}_p\boldsymbol{q} =\frac12\sum_{j=1}^nK_{pj}q_j^2 \end{aligned} $$

正交性使不同模态之间的交叉项全部消失。第 $j$ 阶模态能量为

$$ E_j=\frac12M_{pj}\dot q_j^2+\frac12K_{pj}q_j^2 $$

沿系统运动对时间求导，得

$$ \frac{\mathrm dE_j}{\mathrm dt} =\dot q_j(M_{pj}\ddot q_j+K_{pj}q_j)=0 $$

因此，无外力、无阻尼时，各阶模态分别守恒能量，彼此不交换能量。

（4）振型的定标和模态质量归一化

固有振型乘以任意非零常数仍是同一本征值的振型。若改用 $c\boldsymbol{\varphi}_j$，则模态质量和模态刚度均乘以 $c^2$，二者之比不变。

一种定标方式是将振型的某个非零分量设为 $1$，例如 $\varphi_{1j}=1$。另一种方式是模态质量归一化，定义

$$ \boldsymbol{\varphi}_{Nj}=\frac{\boldsymbol{\varphi}_j}{\sqrt{M_{pj}}} $$

于是

$$ \boldsymbol{\varphi}_{Nj}^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_{Nj} =\frac{\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_j}{M_{pj}}=1 $$

$$ \boldsymbol{\varphi}_{Nj}^{\mathrm T}\boldsymbol{K}\boldsymbol{\varphi}_{Nj} =\frac{K_{pj}}{M_{pj}}=\omega_j^2 $$

由归一化振型组成矩阵 $\boldsymbol{\Phi}_N$，则

$$ \boldsymbol{\Phi}_N^{\mathrm T}\boldsymbol{M}\boldsymbol{\Phi}_N=\boldsymbol{I},\qquad \boldsymbol{\Phi}_N^{\mathrm T}\boldsymbol{K}\boldsymbol{\Phi}_N=\boldsymbol{\Omega} $$

其中单位矩阵 $\boldsymbol{I}$ 与 $\boldsymbol{\Omega}=\operatorname{diag}(\omega_1^2,\ldots,\omega_n^2)$ 均为 $n$ 阶矩阵。

例1
三个质量均为 $m$ 的质点由四根刚度均为 $k$ 的弹簧串联，两端固定。取三个质点向右的位移 $x_1,x_2,x_3$ 为坐标。求固有频率、固有振型，并作模态质量归一化。

解
系统动能和势能为

$$ T=\frac12m(\dot x_1^2+\dot x_2^2+\dot x_3^2) $$

$$ V=\frac12kx_1^2+\frac12k(x_2-x_1)^2+\frac12k(x_3-x_2)^2+\frac12kx_3^2 $$

逐项求偏导，得

$$ \frac{\partial V}{\partial x_1}=k(2x_1-x_2),\qquad \frac{\partial V}{\partial x_2}=k(-x_1+2x_2-x_3),\qquad \frac{\partial V}{\partial x_3}=k(-x_2+2x_3) $$

代入拉格朗日方程，得到

$$ \boldsymbol{M}\ddot{\boldsymbol{x}}+\boldsymbol{K}\boldsymbol{x}=\boldsymbol{0},\qquad \boldsymbol{M}=m\boldsymbol{I},\qquad \boldsymbol{K}=k\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{bmatrix} $$

本征方程为

$$ p(\lambda)=\begin{vmatrix}2k-\lambda m&-k&0\\-k&2k-\lambda m&-k\\0&-k&2k-\lambda m\end{vmatrix}=0 $$

沿第一行展开，得

$$ \begin{aligned} p(\lambda) &=(2k-\lambda m)\left[(2k-\lambda m)^2-k^2\right]-k^2(2k-\lambda m)\\ &=(2k-\lambda m)\left[(2k-\lambda m)^2-2k^2\right] \end{aligned} $$

因此

$$ \lambda_1=(2-\sqrt2)\frac{k}{m},\qquad \lambda_2=\frac{2k}{m},\qquad \lambda_3=(2+\sqrt2)\frac{k}{m} $$

固有频率为

$$ \omega_1=\sqrt{(2-\sqrt2)\frac{k}{m}},\qquad \omega_2=\sqrt{\frac{2k}{m}},\qquad \omega_3=\sqrt{(2+\sqrt2)\frac{k}{m}} $$

对第一阶振型，代入 $\lambda_1$ 并令 $\varphi_{11}=1$，得到

$$ \begin{cases} \sqrt2-\varphi_{21}=0\\ -1+\sqrt2\varphi_{21}-\varphi_{31}=0\\ -\varphi_{21}+\sqrt2\varphi_{31}=0 \end{cases} $$

第一式给出 $\varphi_{21}=\sqrt2$，第二式给出 $\varphi_{31}=1$，并满足第三式，故

$$ \boldsymbol{\varphi}_1=\begin{bmatrix}1\\\sqrt2\\1\end{bmatrix} $$

对第二阶，$2k-\lambda_2m=0$，方程给出 $\varphi_{22}=0$、$\varphi_{32}=-\varphi_{12}$。令 $\varphi_{12}=1$，得

$$ \boldsymbol{\varphi}_2=\begin{bmatrix}1\\0\\-1\end{bmatrix} $$

对第三阶，$2k-\lambda_3m=-\sqrt2k$。令 $\varphi_{13}=1$，由第一行得 $\varphi_{23}=-\sqrt2$，由第二行得 $\varphi_{33}=1$，故

$$ \boldsymbol{\varphi}_3=\begin{bmatrix}1\\-\sqrt2\\1\end{bmatrix} $$

模态质量分别为

$$ M_{p1}=m(1^2+(\sqrt2)^2+1^2)=4m $$

$$ M_{p2}=m(1^2+0^2+(-1)^2)=2m,\qquad M_{p3}=m(1^2+(-\sqrt2)^2+1^2)=4m $$

归一化振型为

$$ \boldsymbol{\varphi}_{N1}=\frac1{\sqrt m}\begin{bmatrix}1/2\\1/\sqrt2\\1/2\end{bmatrix},\qquad \boldsymbol{\varphi}_{N2}=\frac1{\sqrt m}\begin{bmatrix}1/\sqrt2\\0\\-1/\sqrt2\end{bmatrix},\qquad \boldsymbol{\varphi}_{N3}=\frac1{\sqrt m}\begin{bmatrix}1/2\\-1/\sqrt2\\1/2\end{bmatrix} $$

三阶主振动分别为

$$ \boldsymbol{x}_1=a_1\begin{bmatrix}1\\\sqrt2\\1\end{bmatrix}\sin(\omega_1t+\theta_1),\qquad \boldsymbol{x}_2=a_2\begin{bmatrix}1\\0\\-1\end{bmatrix}\sin(\omega_2t+\theta_2) $$

$$ \boldsymbol{x}_3=a_3\begin{bmatrix}1\\-\sqrt2\\1\end{bmatrix}\sin(\omega_3t+\theta_3) $$

第一阶中三个质点同相运动，内部无节点；第二阶中质点 $1$ 与质点 $3$ 反相运动，质点 $2$ 不动，内部有一个节点；第三阶中两端质点同相，中央质点与两端反相，振幅连线在内部有两个节点。

三、自由振动的振型叠加分析方法

1、坐标变换与固有模态的不变性

设可逆坐标变换为 $\boldsymbol{x}=\boldsymbol{A}\boldsymbol{y}$。将其代入运动方程并左乘 $\boldsymbol{A}^{\mathrm T}$，得

$$ \boldsymbol{M}'\ddot{\boldsymbol{y}}+\boldsymbol{K}'\boldsymbol{y}=\boldsymbol{0},\qquad \boldsymbol{M}'=\boldsymbol{A}^{\mathrm T}\boldsymbol{M}\boldsymbol{A},\qquad \boldsymbol{K}'=\boldsymbol{A}^{\mathrm T}\boldsymbol{K}\boldsymbol{A} $$

因此

$$ \begin{aligned} \det(\boldsymbol{K}^\prime-\lambda\boldsymbol{M}^\prime) &=\det\left[\boldsymbol{A}^{\mathrm T}(\boldsymbol{K}-\lambda\boldsymbol{M})\boldsymbol{A}\right]\\ &=(\det\boldsymbol{A})^2\det(\boldsymbol{K}-\lambda\boldsymbol{M}) \end{aligned} $$

由于 $\det\boldsymbol{A}\ne0$，变换前后的本征值相同。若新坐标下的振型为 $\boldsymbol{\psi}_j$，则物理坐标下的振型为 $\boldsymbol{\varphi}_j=\boldsymbol{A}\boldsymbol{\psi}_j$。坐标变换改变振型的分量表示，但不改变固有频率及其所代表的物理运动形式。

取 $\boldsymbol{A}=\boldsymbol{\Phi}$，即可将系统化为相互独立的主振动方程。分别求解各阶主坐标，再将其叠加恢复物理位移，称为振型叠加法或模态叠加法。

2、由初始条件确定主坐标

对 $\omega_j>0$，有

$$ \boldsymbol{x}(t)=\sum_{j=1}^n\boldsymbol{\varphi}_j \left(A_j\cos\omega_jt+B_j\sin\omega_jt\right) $$

代入初始条件，得到

$$ \boldsymbol{x}_0=\sum_{j=1}^n\boldsymbol{\varphi}_jA_j,\qquad \dot{\boldsymbol{x}}_0=\sum_{j=1}^n\boldsymbol{\varphi}_j\omega_jB_j $$

第一式左乘 $\boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{M}$，得到

$$ \boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0 =\sum_{j=1}^n\boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_jA_j =M_{pi}A_i $$

同理，第二式给出

$$ \boldsymbol{\varphi}_i^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}}_0 =\omega_iM_{pi}B_i $$

于是

$$ A_j=\frac{\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0}{M_{pj}},\qquad B_j=\frac{\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}}_0}{\omega_jM_{pj}} $$

主坐标位移为

$$ q_j(t)=\frac1{M_{pj}}\left( \boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0\cos\omega_jt +\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}}_0 \frac{\sin\omega_jt}{\omega_j}\right) $$

依次求导，得

$$ \dot q_j(t)=\frac1{M_{pj}}\left( -\omega_j\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0\sin\omega_jt +\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}}_0\cos\omega_jt \right) $$

$$ \ddot q_j(t)=-\omega_j^2q_j(t) $$

对于零频模态，由

$$ \lim_{\omega\to0}\cos\omega t=1,\qquad \lim_{\omega\to0}\frac{\sin\omega t}{\omega}=t $$

得到

$$ q_j(t)=\frac1{M_{pj}}\left( \boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0 +\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}}_0\,t\right) \qquad(\omega_j=0) $$

该式也可由 $\ddot q_j=0$ 直接积分得到。它描述恒定位移或匀速运动，而非往复振动。

最终物理位移为

$$ \boldsymbol{x}(t)=\sum_{j=1}^n\boldsymbol{\varphi}_jq_j(t) $$

对正频模态，也可合写为

$$ \boldsymbol{x}(t)=\sum_{j=1}^n\frac{\boldsymbol{\varphi}_j}{M_{pj}} \left(\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0\cos\omega_jt +\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}}_0\frac{\sin\omega_jt}{\omega_j}\right) $$

含零频时，将相应求和项按上述极限解释即可。实际计算依次求出固有模态、投影初始条件、叠加各阶响应。

3、单阶主振动、零频运动和周期性

若初始位移仅沿第 $i$ 阶振型，即

$$ \boldsymbol{x}_0=a\boldsymbol{\varphi}_i,\qquad \dot{\boldsymbol{x}}_0=\boldsymbol{0} $$

则由质量正交性，$A_i=a$，其余 $A_j=0$，所有 $B_j=0$，故

$$ \boldsymbol{x}(t)=a\boldsymbol{\varphi}_i\cos\omega_it $$

系统只作第 $i$ 阶主振动。

零频模态满足 $K_{pj}=0$，其主坐标为 $q_j=a_j+b_jt$。相应势能恒为零，动能为常量

$$ V_j=0,\qquad T_j=\frac12M_{pj}b_j^2 $$

多个正频模态叠加时，响应存在共同周期 $T>0$ 的条件是所有实际参与响应的频率满足

$$ \omega_jT=2\pi N_j,\qquad N_j\in\mathbb{N} $$

即各参与频率之比均为有理数，也可表述为存在非零有理数 $c_j$，使

$$ c_1\omega_1=c_2\omega_2=\cdots=c_s\omega_s $$

这里 $s$ 为实际参与响应的正频模态数。若存在共同周期，由位移和速度的唯一模态展开可知，各非零模态均须在 $T$ 后回到同一状态，故满足上述条件；反之，这些条件使每个参与模态在 $T$ 后重复，叠加响应也随之重复。只有一阶正频模态参与时，响应为简谐振动；多阶不可通约频率共同参与时，响应不具有严格周期。

例2
对例1的三质点弹簧系统，给定初始条件

$$ \boldsymbol{x}_0=\begin{bmatrix}0\\1\\0\end{bmatrix},\qquad \dot{\boldsymbol{x}}_0=\begin{bmatrix}0\\0\\0\end{bmatrix} $$

用振型叠加法求自由振动响应。

解
由例1，三个振型及其模态质量为

$$ \boldsymbol{\varphi}_1=\begin{bmatrix}1\\\sqrt2\\1\end{bmatrix},\qquad \boldsymbol{\varphi}_2=\begin{bmatrix}1\\0\\-1\end{bmatrix},\qquad \boldsymbol{\varphi}_3=\begin{bmatrix}1\\-\sqrt2\\1\end{bmatrix} $$

$$ M_{p1}=4m,\qquad M_{p2}=2m,\qquad M_{p3}=4m $$

初始位移在三阶振型上的质量投影分别为

$$ \boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0=m\sqrt2,\qquad \boldsymbol{\varphi}_2^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0=0,\qquad \boldsymbol{\varphi}_3^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0=-m\sqrt2 $$

所以

$$ A_1=\frac{\sqrt2}{4},\qquad A_2=0,\qquad A_3=-\frac{\sqrt2}{4} $$

由于初速度为零，$B_1=B_2=B_3=0$，因而

$$ q_1=\frac{\sqrt2}{4}\cos\omega_1t,\qquad q_2=0,\qquad q_3=-\frac{\sqrt2}{4}\cos\omega_3t $$

代回 $\boldsymbol{x}=\boldsymbol{\Phi}\boldsymbol{q}$，得

$$ \boldsymbol{x}(t)=\frac{\sqrt2}{4} \begin{bmatrix}1\\\sqrt2\\1\end{bmatrix}\cos\omega_1t -\frac{\sqrt2}{4}\begin{bmatrix}1\\-\sqrt2\\1\end{bmatrix}\cos\omega_3t $$

即

$$ \begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix} =\begin{bmatrix} \dfrac{\sqrt2}{4}(\cos\omega_1t-\cos\omega_3t)\\ \dfrac12(\cos\omega_1t+\cos\omega_3t)\\ \dfrac{\sqrt2}{4}(\cos\omega_1t-\cos\omega_3t) \end{bmatrix} $$

取 $t=0$ 得 $[0,1,0]^{\mathrm T}$，对时间求导后取 $t=0$ 得零向量，满足全部初始条件。

由

$$ \frac{\omega_3}{\omega_1} =\sqrt{\frac{2+\sqrt2}{2-\sqrt2}} =\sqrt{3+2\sqrt2}=1+\sqrt2 $$

可知两阶参与频率之比为无理数，因此该响应没有严格周期。

对任意响应，左乘 $\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}$ 可得

$$ q_j(t)=\frac{\boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{x}(t)}{M_{pj}} $$

主坐标表示在给定振型定标下，第 $j$ 阶模态对位移响应的贡献系数。

四、重频、零频和高频

1、与重频对应的振型向量的正交性

设 $\omega_1=\omega_2=\cdots=\omega_r$，取该本征子空间的一组线性无关振型 $\boldsymbol{\varphi}_1,\ldots,\boldsymbol{\varphi}_r$。这些振型不一定彼此质量正交，但均与不同本征值的振型质量正交。

实对称矩阵 $\widetilde{\boldsymbol{K}}$ 可以正交对角化，故重本征值对应的本征子空间具有与其代数重数相同的维数。经 $\boldsymbol{L}^{-\mathrm T}$ 变换后，广义本征问题也具有足够多的线性无关振型。可在此子空间内按质量内积进行正交化。

第一步，取

$$ \widetilde{\boldsymbol{\varphi}}_1=\boldsymbol{\varphi}_1 $$

第二步，令

$$ \widetilde{\boldsymbol{\varphi}}_2=c_1\widetilde{\boldsymbol{\varphi}}_1+\boldsymbol{\varphi}_2 $$

要求 $\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}\widetilde{\boldsymbol{\varphi}}_2=0$，即

$$ c_1\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}\widetilde{\boldsymbol{\varphi}}_1 +\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_2=0 $$

因此

$$ c_1=-\frac{\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_2} {\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}\widetilde{\boldsymbol{\varphi}}_1} $$

第三步，令

$$ \widetilde{\boldsymbol{\varphi}}_3=c_2\widetilde{\boldsymbol{\varphi}}_1+c_3\widetilde{\boldsymbol{\varphi}}_2+\boldsymbol{\varphi}_3 $$

分别左乘 $\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}$、$\widetilde{\boldsymbol{\varphi}}_2^{\mathrm T}\boldsymbol{M}$，并利用已得的正交性，有

$$ 0=c_2\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}\widetilde{\boldsymbol{\varphi}}_1 +\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_3 $$

$$ 0=c_3\widetilde{\boldsymbol{\varphi}}_2^{\mathrm T}\boldsymbol{M}\widetilde{\boldsymbol{\varphi}}_2 +\widetilde{\boldsymbol{\varphi}}_2^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_3 $$

故

$$ c_2=-\frac{\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_3} {\widetilde{\boldsymbol{\varphi}}_1^{\mathrm T}\boldsymbol{M}\widetilde{\boldsymbol{\varphi}}_1},\qquad c_3=-\frac{\widetilde{\boldsymbol{\varphi}}_2^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_3} {\widetilde{\boldsymbol{\varphi}}_2^{\mathrm T}\boldsymbol{M}\widetilde{\boldsymbol{\varphi}}_2} $$

依次进行，得到一般公式

$$ \widetilde{\boldsymbol{\varphi}}_j=\boldsymbol{\varphi}_j -\sum_{i=1}^{j-1} \frac{\widetilde{\boldsymbol{\varphi}}_i^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_j} {\widetilde{\boldsymbol{\varphi}}_i^{\mathrm T}\boldsymbol{M}\widetilde{\boldsymbol{\varphi}}_i} \widetilde{\boldsymbol{\varphi}}_i,\qquad j=2,\ldots,r $$

由于初始向量线性无关，减去前面向量的投影后仍非零；质量正定性保证所有分母为正。

这些线性组合仍为同一重本征值的振型。事实上，若 $\lambda=\omega_1^2$，则

$$ \boldsymbol{K}\left(\sum_{i=1}^ra_i\boldsymbol{\varphi}_i\right) =\lambda\boldsymbol{M}\left(\sum_{i=1}^ra_i\boldsymbol{\varphi}_i\right) $$

因此正交化不改变本征值。同时

$$ \widetilde{\boldsymbol{\varphi}}_i^{\mathrm T}\boldsymbol{K}\widetilde{\boldsymbol{\varphi}}_j =\lambda\widetilde{\boldsymbol{\varphi}}_i^{\mathrm T}\boldsymbol{M}\widetilde{\boldsymbol{\varphi}}_j=0 \qquad(i\ne j) $$

所以质量正交化后也具有刚度正交性。再作模态质量归一化，即可得到满足 $\boldsymbol{\Phi}_N^{\mathrm T}\boldsymbol{M}\boldsymbol{\Phi}_N=\boldsymbol{I}$ 的振型矩阵。数值求解后，可通过本征方程残差及两类正交关系检验结果。

2、零频主振动

（1）零频振型与零频个数

令 $\lambda=\omega^2=0$，广义本征方程化为

$$ \boldsymbol{K}\boldsymbol{\varphi}=\boldsymbol{0} $$

存在非零零频振型的充要条件为 $\det\boldsymbol{K}=0$。由秩与零空间维数的关系，零频模态个数为

$$ r=n-\operatorname{rank}(\boldsymbol{K}) $$

其与零本征值的重数相同，这是因为系统可化为实对称本征问题。对于无内部机构的自由三维弹性体，零频对应三个刚体平移和三个刚体转动；约束可减少刚体模态数，内部机构也可能增加零能运动模态。

若仅有一个零频，记 $\omega_1=0$，则

$$ \ddot q_1=0,\qquad q_1=a_1+b_1t $$

其中

$$ a_1=\frac{\boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0}{M_{p1}},\qquad b_1=\frac{\boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}}_0}{M_{p1}} $$

相应物理位移为

$$ \boldsymbol{x}_1(t)=\boldsymbol{\varphi}_1q_1(t) =\frac{\boldsymbol{\varphi}_1}{M_{p1}} \left[\boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{x}_0 +(\boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\dot{\boldsymbol{x}}_0)t\right] $$

即使这一分量不是往复运动，它仍是系统运动方程的真实解。

（2）通过降阶去除零频

若只研究不含刚体运动的振动响应，可在叠加时略去零频项；也可先施加质量正交约束

$$ \boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{x}=0 $$

因为 $q_1=\boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{x}/M_{p1}$，该约束恰好令零频主坐标为零。若有 $r$ 个零频，则施加

$$ \boldsymbol{\varphi}_j^{\mathrm T}\boldsymbol{M}\boldsymbol{x}=0,\qquad j=1,\ldots,r $$

取这些约束的解空间的一组基，按列组成 $n\times(n-r)$ 矩阵 $\boldsymbol{H}$，令

$$ \boldsymbol{x}=\boldsymbol{H}\overline{\boldsymbol{x}} $$

代入动能、势能并应用拉格朗日方程，得到降阶系统

$$ \overline{\boldsymbol{M}}\ddot{\overline{\boldsymbol{x}}} +\overline{\boldsymbol{K}}\overline{\boldsymbol{x}}=\boldsymbol{0},\qquad \overline{\boldsymbol{M}}=\boldsymbol{H}^{\mathrm T}\boldsymbol{M}\boldsymbol{H},\qquad \overline{\boldsymbol{K}}=\boldsymbol{H}^{\mathrm T}\boldsymbol{K}\boldsymbol{H} $$

定理 若上述约束去除了全部零频模态，则降阶系统的固有频率恰为原系统的全部非零固有频率。

证明
将非零频振型组成矩阵 $\boldsymbol{\Phi}_v$。由质量正交性，这些振型都满足零频约束；又因线性无关且个数为 $n-r$，它们构成该约束空间的一组基。

$\boldsymbol{H}$ 与 $\boldsymbol{\Phi}_v$ 是同一空间的两组基，故存在可逆矩阵 $\boldsymbol{C}$，使

$$ \boldsymbol{H}=\boldsymbol{\Phi}_v\boldsymbol{C} $$

因而

$$ \overline{\boldsymbol{M}}=\boldsymbol{C}^{\mathrm T}\boldsymbol{M}_{pv}\boldsymbol{C},\qquad \overline{\boldsymbol{K}}=\boldsymbol{C}^{\mathrm T}\boldsymbol{K}_{pv}\boldsymbol{C} $$

其中 $\boldsymbol{M}_{pv}$、$\boldsymbol{K}_{pv}$ 是非零频模态的对角质量矩阵和刚度矩阵。因此

$$ \begin{aligned} \det(\overline{\boldsymbol{K}}-\lambda\overline{\boldsymbol{M}}) &=(\det\boldsymbol{C})^2\det(\boldsymbol{K}_{pv}-\lambda\boldsymbol{M}_{pv})\\ &=(\det\boldsymbol{C})^2\prod_{j=r+1}^nM_{pj}(\omega_j^2-\lambda) \end{aligned} $$

由于 $\det\boldsymbol{C}\ne0$、$M_{pj}>0$，降阶本征值恰为 $\omega_{r+1}^2,\ldots,\omega_n^2$。

（3）拉格朗日乘子与约束动力学

仅有一个零频时，构造带约束的拉格朗日函数

$$ L=T-V-\lambda\boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{x} $$

这里 $\lambda=\lambda(t)$ 为拉格朗日乘子，不再表示本征值。对 $\boldsymbol{x}$ 和 $\lambda$ 分别取变分，得到

$$ \begin{cases} \boldsymbol{M}\ddot{\boldsymbol{x}}+\boldsymbol{K}\boldsymbol{x} +\lambda\boldsymbol{M}\boldsymbol{\varphi}_1=\boldsymbol{0}\\ \boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{x}=0 \end{cases} $$

约束对应的雅可比矩阵为 $\boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}$，其转置为 $\boldsymbol{M}\boldsymbol{\varphi}_1$。

令时间步长为 $\Delta t$，采用平均加速度格式

$$ \boldsymbol{x}_{t+\Delta t} =\boldsymbol{x}_t+\Delta t\dot{\boldsymbol{x}}_t +\frac{\Delta t^2}{4}(\ddot{\boldsymbol{x}}_t+\ddot{\boldsymbol{x}}_{t+\Delta t}) $$

$$ \dot{\boldsymbol{x}}_{t+\Delta t} =\dot{\boldsymbol{x}}_t+\frac{\Delta t}{2}(\ddot{\boldsymbol{x}}_t+\ddot{\boldsymbol{x}}_{t+\Delta t}) $$

由第一式解出

$$ \ddot{\boldsymbol{x}}_{t+\Delta t} =\frac4{\Delta t^2}(\boldsymbol{x}_{t+\Delta t}-\boldsymbol{x}_t) -\frac4{\Delta t}\dot{\boldsymbol{x}}_t-\ddot{\boldsymbol{x}}_t $$

将其代入 $t+\Delta t$ 时刻的动力学方程，得

$$ \left(\frac4{\Delta t^2}\boldsymbol{M}+\boldsymbol{K}\right)\boldsymbol{x}_{t+\Delta t} +\lambda_{t+\Delta t}\boldsymbol{M}\boldsymbol{\varphi}_1 =\frac4{\Delta t^2}\boldsymbol{M}\boldsymbol{x}_t +\frac4{\Delta t}\boldsymbol{M}\dot{\boldsymbol{x}}_t +\boldsymbol{M}\ddot{\boldsymbol{x}}_t $$

利用当前时刻的平衡关系

$$ \boldsymbol{M}\ddot{\boldsymbol{x}}_t =-\boldsymbol{K}\boldsymbol{x}_t-\lambda_t\boldsymbol{M}\boldsymbol{\varphi}_1 $$

消去旧时刻加速度，最终得到递推方程

$$ \begin{bmatrix} \dfrac4{\Delta t^2}\boldsymbol{M}+\boldsymbol{K}&\boldsymbol{M}\boldsymbol{\varphi}_1\\ \boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}&0 \end{bmatrix} \begin{bmatrix}\boldsymbol{x}_{t+\Delta t}\\\lambda_{t+\Delta t}\end{bmatrix} =\begin{bmatrix} \left(\dfrac4{\Delta t^2}\boldsymbol{M}-\boldsymbol{K}\right)\boldsymbol{x}_t +\dfrac4{\Delta t}\boldsymbol{M}\dot{\boldsymbol{x}}_t -\lambda_t\boldsymbol{M}\boldsymbol{\varphi}_1\\0 \end{bmatrix} $$

求解该线性方程组后，速度可直接由

$$ \dot{\boldsymbol{x}}_{t+\Delta t} =\frac2{\Delta t}(\boldsymbol{x}_{t+\Delta t}-\boldsymbol{x}_t)-\dot{\boldsymbol{x}}_t $$

更新。初始位移和初始速度须满足对应的位移约束和速度约束。若给定初值包含刚体成分，应先减去其零频模态投影，再计算去除零频后的响应。

3、频率很大情况

对固定的有限维系统，若 $\boldsymbol{M}$ 正定、$\boldsymbol{K}$ 的元素有限，则

$$ \lambda_j=\lambda_j(\widetilde{\boldsymbol{K}}) $$

均为有限数，因而所有固有频率均有限。连续系统则可能具有无穷多个固有频率，并随阶数增大而趋于无穷。

由本征方程可写出

$$ \boldsymbol{M}\boldsymbol{\varphi}_j =\frac1{\omega_j^2}\boldsymbol{K}\boldsymbol{\varphi}_j $$

频率较高时，惯性项与刚度项仍保持平衡；振型长度本身可以任意定标，不能仅由高频推出振型接近零。

自由振动中，第 $j$ 阶主坐标为

$$ q_j(t)=q_j(0)\cos\omega_jt+\dot q_j(0)\frac{\sin\omega_jt}{\omega_j} $$

初速度所贡献的位移幅值含有 $1/\omega_j$，而初位移项没有这一衰减因子。因此，高频模态能否忽略，须结合初始条件在该模态上的投影及所关心的响应量判断。

若只保留前 $s$ 阶模态，记位移误差为

$$ \boldsymbol{e}=\boldsymbol{x}-\sum_{j=1}^s\boldsymbol{\varphi}_jq_j =\sum_{j=s+1}^n\boldsymbol{\varphi}_jq_j $$

则由质量正交性

$$ \boldsymbol{e}^{\mathrm T}\boldsymbol{M}\boldsymbol{e} =\sum_{j=s+1}^nM_{pj}q_j^2 $$

只有被截去模态的贡献足够小时，模态截断才可满足所需精度。

例3
三个圆盘绕共同轴线转动，各盘绕轴转动惯量均为 $J$，相邻圆盘之间两段轴的扭转刚度均为 $k$。三个转角依次为 $\alpha_1,\alpha_2,\alpha_3$。

1. 求固有频率和固有振型，验证振型正交性。
2. 通过降阶求非零固有模态。
3. 给定初值 $\boldsymbol{x}_0=[1,2,3]^{\mathrm T}$、$\dot{\boldsymbol{x}}_0=[1,0,0]^{\mathrm T}$，取 $J=k=1$，求自由响应，并计算去除零频贡献后的数值响应。

解
（1）固有频率、固有振型与正交性。

取 $\boldsymbol{x}=[\alpha_1,\alpha_2,\alpha_3]^{\mathrm T}$，有

$$ T=\frac12J(\dot\alpha_1^2+\dot\alpha_2^2+\dot\alpha_3^2),\qquad V=\frac12k\left[(\alpha_1-\alpha_2)^2+(\alpha_3-\alpha_2)^2\right] $$

拉格朗日方程给出

$$ \begin{cases} J\ddot\alpha_1+k(\alpha_1-\alpha_2)=0\\ J\ddot\alpha_2+k(2\alpha_2-\alpha_1-\alpha_3)=0\\ J\ddot\alpha_3+k(\alpha_3-\alpha_2)=0 \end{cases} $$

因此

$$ \boldsymbol{M}=J\boldsymbol{I},\qquad \boldsymbol{K}=k\begin{bmatrix}1&-1&0\\-1&2&-1\\0&-1&1\end{bmatrix} $$

频率方程为

$$ p(\omega)=\begin{vmatrix} k-J\omega^2&-k&0\\ -k&2k-J\omega^2&-k\\ 0&-k&k-J\omega^2 \end{vmatrix}=0 $$

令 $\mu=J\omega^2/k$，展开得

$$ \begin{aligned} p(\omega) &=k^3\left[(1-\mu)((2-\mu)(1-\mu)-1)-(1-\mu)\right]\\ &=k^3(1-\mu)\left[(2-\mu)(1-\mu)-2\right]\\ &=k^3(1-\mu)\mu(\mu-3) \end{aligned} $$

因而

$$ \omega_1=0,\qquad\omega_2=\sqrt{\frac{k}{J}},\qquad\omega_3=\sqrt{\frac{3k}{J}} $$

对 $\mu=0$，第一、三行分别给出 $\varphi_{11}=\varphi_{21}$、$\varphi_{31}=\varphi_{21}$；令第一分量为 $1$，得

$$ \boldsymbol{\varphi}_1=\begin{bmatrix}1\\1\\1\end{bmatrix} $$

对 $\mu=1$，第一行给出 $\varphi_{22}=0$，第二行给出 $\varphi_{32}=-\varphi_{12}$，故

$$ \boldsymbol{\varphi}_2=\begin{bmatrix}1\\0\\-1\end{bmatrix} $$

对 $\mu=3$，第一行给出 $\varphi_{23}=-2\varphi_{13}$，第三行给出 $\varphi_{33}=\varphi_{13}$，故

$$ \boldsymbol{\varphi}_3=\begin{bmatrix}1\\-2\\1\end{bmatrix} $$

验证质量正交性：

$$ \boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_2=J(1+0-1)=0 $$

$$ \boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_3=J(1-2+1)=0,\qquad \boldsymbol{\varphi}_2^{\mathrm T}\boldsymbol{M}\boldsymbol{\varphi}_3=J(1+0-1)=0 $$

又因

$$ \boldsymbol{K}\boldsymbol{\varphi}_1=\boldsymbol{0},\qquad \boldsymbol{K}\boldsymbol{\varphi}_2=k\boldsymbol{\varphi}_2,\qquad \boldsymbol{K}\boldsymbol{\varphi}_3=3k\boldsymbol{\varphi}_3 $$

可知异阶之间的刚度内积也全部为零。模态质量和模态刚度分别为

$$ (M_{p1},M_{p2},M_{p3})=(3J,2J,6J),\qquad (K_{p1},K_{p2},K_{p3})=(0,2k,18k) $$

第一阶为三个圆盘同角度的刚体转动，不引起轴的扭转；其主坐标满足 $\ddot q_1=0$，相应运动为

$$ \boldsymbol{x}_1=\begin{bmatrix}1\\1\\1\end{bmatrix}(a_1+b_1t) $$

（2）利用质量正交约束降阶。

由 $\boldsymbol{\varphi}_1^{\mathrm T}\boldsymbol{M}\boldsymbol{x}=0$，有

$$ \alpha_1+\alpha_2+\alpha_3=0,\qquad\alpha_3=-(\alpha_1+\alpha_2) $$

代回动能、势能，得到

$$ T=\frac12J\left[\dot\alpha_1^2+\dot\alpha_2^2+(\dot\alpha_1+\dot\alpha_2)^2\right] $$

$$ V=\frac12k\left[(\alpha_1-\alpha_2)^2+(\alpha_1+2\alpha_2)^2\right] $$

展开为

$$ T=\frac12J(2\dot\alpha_1^2+2\dot\alpha_1\dot\alpha_2+2\dot\alpha_2^2),\qquad V=\frac12k(2\alpha_1^2+2\alpha_1\alpha_2+5\alpha_2^2) $$

因此，令 $\overline{\boldsymbol{x}}=[\alpha_1,\alpha_2]^{\mathrm T}$，有

$$ \overline{\boldsymbol{M}}=J\begin{bmatrix}2&1\\1&2\end{bmatrix},\qquad \overline{\boldsymbol{K}}=k\begin{bmatrix}2&1\\1&5\end{bmatrix} $$

降阶频率方程为

$$ \begin{aligned} \det(\overline{\boldsymbol{K}}-\overline\omega^2\overline{\boldsymbol{M}}) &=(2k-2J\overline\omega^2)(5k-2J\overline\omega^2)-(k-J\overline\omega^2)^2\\ &=3(k-J\overline\omega^2)(3k-J\overline\omega^2)=0 \end{aligned} $$

故

$$ \overline\omega_1=\sqrt{\frac{k}{J}},\qquad \overline\omega_2=\sqrt{\frac{3k}{J}} $$

对 $\overline\omega_1^2=k/J$，系数矩阵为 $k\begin{bmatrix}0&0\\0&3\end{bmatrix}$，故取

$$ \overline{\boldsymbol{\varphi}}_1=\begin{bmatrix}1\\0\end{bmatrix} $$

对 $\overline\omega_2^2=3k/J$，系数矩阵为 $k\begin{bmatrix}-4&-2\\-2&-1\end{bmatrix}$，故取

$$ \overline{\boldsymbol{\varphi}}_2=\begin{bmatrix}1\\-2\end{bmatrix} $$

利用 $\alpha_3=-(\alpha_1+\alpha_2)$ 恢复第三分量，得到的完整振型分别为 $[1,0,-1]^{\mathrm T}$ 和 $[1,-2,1]^{\mathrm T}$，与原系统非零频振型一致。

（3）解析响应及约束数值响应。

取 $J=k=1$，则 $\omega_2=1$、$\omega_3=\sqrt3$。零频主坐标的初值为

$$ q_1(0)=\frac{1+2+3}{3}=2,\qquad \dot q_1(0)=\frac{1+0+0}{3}=\frac13 $$

因此

$$ q_1(t)=2+\frac t3 $$

第二阶的初值为

$$ q_2(0)=\frac{1-3}{2}=-1,\qquad \dot q_2(0)=\frac{1-0}{2}=\frac12 $$

第三阶的初值为

$$ q_3(0)=\frac{1-4+3}{6}=0,\qquad \dot q_3(0)=\frac{1-0+0}{6}=\frac16 $$

所以

$$ q_2(t)=-\cos t+\frac12\sin t,\qquad q_3(t)=\frac1{6\sqrt3}\sin(\sqrt3t) $$

将三阶响应叠加，得到完整解

$$ \boldsymbol{x}(t)=\begin{bmatrix} 2-\cos t+\dfrac t3+\dfrac12\sin t+\dfrac1{6\sqrt3}\sin(\sqrt3t)\\ 2+\dfrac t3-\dfrac1{3\sqrt3}\sin(\sqrt3t)\\ 2+\cos t+\dfrac t3-\dfrac12\sin t+\dfrac1{6\sqrt3}\sin(\sqrt3t) \end{bmatrix} $$

去除刚体分量 $\boldsymbol{\varphi}_1q_1(t)$ 后，将振动响应记为 $\boldsymbol{x}_v$，则

$$ \boldsymbol{x}_v(t)=\begin{bmatrix} -\cos t+\dfrac12\sin t+\dfrac1{6\sqrt3}\sin(\sqrt3t)\\ -\dfrac1{3\sqrt3}\sin(\sqrt3t)\\ \cos t-\dfrac12\sin t+\dfrac1{6\sqrt3}\sin(\sqrt3t) \end{bmatrix} $$

其初始条件应同步去除刚体分量：

$$ \boldsymbol{x}_{v0}=\boldsymbol{x}_0-2\boldsymbol{\varphi}_1 =\begin{bmatrix}-1\\0\\1\end{bmatrix} $$

$$ \dot{\boldsymbol{x}}_{v0}=\dot{\boldsymbol{x}}_0-\frac13\boldsymbol{\varphi}_1 =\begin{bmatrix}2/3\\-1/3\\-1/3\end{bmatrix} $$

两者均满足三个分量之和为零。令

$$ \Delta t=\frac{2\pi/\omega_3}{20}=\frac{\pi}{10\sqrt3}\approx0.181380,\qquad \lambda_0=0 $$

把 $\boldsymbol{x}_{v0}$、$\dot{\boldsymbol{x}}_{v0}$ 作为约束递推方程的初值，每步求解

$$ \begin{bmatrix} 4\boldsymbol{I}/\Delta t^2+\boldsymbol{K}&\boldsymbol{\varphi}_1\\ \boldsymbol{\varphi}_1^{\mathrm T}&0 \end{bmatrix} \begin{bmatrix}\boldsymbol{x}_{v,t+\Delta t}\\\lambda_{t+\Delta t}\end{bmatrix} =\begin{bmatrix} (4\boldsymbol{I}/\Delta t^2-\boldsymbol{K})\boldsymbol{x}_{v,t} +4\dot{\boldsymbol{x}}_{v,t}/\Delta t-\lambda_t\boldsymbol{\varphi}_1\\0 \end{bmatrix} $$

随后更新速度。圆盘 $1$ 的部分计算结果如下，表中的误差为转角绝对误差。

步数	时间 $t$	解析转角	数值转角	绝对误差
0	0.000000	-1.000000	-1.000000	0.000000
20	3.627599	0.650656	0.654724	0.004068
40	7.255197	-0.150628	-0.182308	0.031680
60	10.882796	-0.384284	-0.367592	0.016692
80	14.510395	0.830201	0.780448	0.049752
100	18.137994	-1.083851	-1.093193	0.009341
200	36.275987	-0.641694	-0.595007	0.046687

在 $0\le t\le200\Delta t$ 内，三个转角之和的最大绝对残差小于 $10^{-14}$，说明数值响应满足去除刚体模态的约束。给定步长下，圆盘 $1$ 的最大转角绝对误差约为 $0.142$，长时间积分可观察到相位误差累积。