第七节 提高计算速度（理论篇）

1、耦合状态强迫振动法原理

1.频域识别方法基本原理

驱动断面做二自由度耦合强迫振动，其中竖向和扭转均为稳态简谐振动。即： $$h(t)=h_0e^{i(\omega_ht+\phi_h)}\\ \alpha(t)=\alpha_0e^{i(\omega_\alpha+\phi_\alpha)} \tag{1}$$ 将上式带入Scanlan自激力模型式(2)，可以得到式(3) $$L_{se}=\rho U^2B[K_hH_1^*\frac{\dot h}{U}+K_\alpha H_2^*\frac{B\dot \alpha}{U}+K_\alpha ^2H_3^*\alpha+K_h^2H_4^*\frac{h}{B}] \\ M_{se}=\rho U^2B^2[K_hA_1^*\frac{\dot h}{U}+K_\alpha A_2^*\frac{B\dot \alpha}{U}+K_\alpha ^2A_3^*\alpha+K_h^2A_4^*\frac{h}{B}] \tag{2}$$

$$L_{se}=\rho U^2B[K_h^2H_1^*\frac{h}{B}i+K_\alpha^2H_2^*\alpha i+K_\alpha^2H_3^*\alpha+K_h^2H_4^*\frac{h}{B}]\\ M_{se}=\rho U^2B^2[K_h^2A_1^*\frac{h}{B}i+K_\alpha^2A_2^*\alpha i+K_\alpha^2A_3^*\alpha+K_h^2A_4^*\frac{h}{B}] \tag{3}$$

同《第三节 FLUENT计算颤振导数》中介绍的分析方法，自激气动力可以描述为： $$L_{se}(t)=L_1e^{i(\omega_ht+\phi_{L1})}+L_2e^{i(\omega_\alpha t+\phi_{L2})}\\ M_{se}(t)=M_1e^{i(\omega_ht+\phi_{M_1})}+M_2e^{i(\omega_\alpha t+\phi_{M_2})}\\ \tag{4}$$ 将式(4)代入式(3)，并对两端分别进行傅里叶变换。两端的复数力在强迫振动频率处相等，有： $$H_1^*(K_h)=\frac{1}{\rho B^2h_0\omega _h^2}\text{Im}[L_1e^{i(\phi_{L1}-\phi_h)}]\\ H_4^*(K_h)=\frac{1}{\rho B^2h_0\omega _h^2}\text{Re}[L_1e^{i(\phi_{L1}-\phi_h)}]\\ A_1^*(K_h)=\frac{1}{\rho B^3h_0\omega _h^2}\text{Im}[M_1e^{i(\phi_{M1}-\phi_h)}]\\ H_4^*(K_h)=\frac{1}{\rho B^3h_0\omega _h^2}\text{Re}[M_1e^{i(\phi_{M1}-\phi_h)}]\\ H_2^*(K_\alpha)=\frac{1}{\rho B^3\alpha_0\omega _\alpha^2}\text{Im}[L_2e^{i(\phi_{L2}-\phi_h)}]\\ H_3^*(K_\alpha)=\frac{1}{\rho B^3\alpha_0\omega _\alpha^2}\text{Re}[L_2e^{i(\phi_{L2}-\phi_h)}]\\ A_2^*(K_\alpha)=\frac{1}{\rho B^4\alpha_0\omega _\alpha^2}\text{Im}[M_2e^{i(\phi_{M2}-\phi_h)}]\\ H_3^*(K_\alpha)=\frac{1}{\rho B^4\alpha_0\omega _\alpha^2}\text{Re}[M_2e^{i(\phi_{M2}-\phi_h)}]\\$$ 关于自激力的提取

不管是试验还是数值模拟，采集得到的气动力中除了自激力外还存在其他周期性荷载如惯性力、周期性机械干扰、模型相对静止空气运动产生的气动力等。因此做出如下假定：强迫振动模型在特定风速下受到的动态力除自激气动力外其他周期性荷载与零风速下保持一致。即： $$L_{se}(t)=L_{dy,U}(t)-L_{dy,0}(t)\\ M_{se}(t)=M_{dy,U}(t)-M_{dy,0}(t)\\ \tag{6}$$ 由于节段模型即使在强迫简谐振动中，其振动时程信号也不是完全相同，因此不能够采用将所测得的动荷载时程信号按照对应相位相减的方法得到自激力，而是应该将足够长的动荷载时程信号转化到频域内，得到对应于振动频率$\omegah,\omega\alpha$的复数形式力（符号表示为$F^{'}$）相减来得到自激力，即： $$\text{Re}[L_{se}^{'}(\omega_h,\omega_\alpha)]=\text{Re}[L_{dy}^{'}(\omega_h,\omega_\alpha)]-\text{Re}[L_{dy,0}^{'}(\omega_h,\omega_\alpha)]\\ \text{Im}[L_{se}^{'}(\omega_h,\omega_\alpha)]=\text{Im}[L_{dy}^{'}(\omega_h,\omega_\alpha)]-\text{Im}[L_{dy,0}^{'}(\omega_h,\omega_\alpha)]\\ \text{Re}[M_{se}^{'}(\omega_h,\omega_\alpha)]=\text{Re}[M_{dy}^{'}(\omega_h,\omega_\alpha)]-\text{Re}[M_{dy,0}^{'}(\omega_h,\omega_\alpha)]\\ \text{Im}[M_{se}^{'}(\omega_h,\omega_\alpha)]=\text{Im}[M_{dy}^{'}(\omega_h,\omega_\alpha)]-\text{Im}[M_{dy,0}^{'}(\omega_h,\omega_\alpha)]\\ \tag{7}$$

2.时域识别方法基本原理

同《第三节 FLUENT计算颤振导数》中介绍的分析方法：基于式(2)，令 $$\begin{bmatrix} H_1\\ H_2\\ H_3\\ H_4 \end{bmatrix} =\rho U^2B \begin{bmatrix} \frac{K_h}{U}\\ & \frac{K_\alpha B}{U}\\ & & K_\alpha^2\\ & & & \frac{K_h^2}{B} \end{bmatrix} \cdot \begin{bmatrix} H_1^*\\ H_2^*\\ H_3^*\\ H_4^* \end{bmatrix} \tag{8}$$

$$\begin{bmatrix} A_1\\ A_2\\ A_3\\ A_4 \end{bmatrix} =\rho U^2B^2 \begin{bmatrix} \frac{K_h}{U}\\ & \frac{K_\alpha B}{U}\\ & & K_\alpha^2\\ & & & \frac{K_h^2}{B} \end{bmatrix} \cdot \begin{bmatrix} A_1^*\\ A_2^*\\ A_3^*\\ A_4^* \end{bmatrix} \tag{9}$$

则，式(2)变为： $$L_{se}=H_1\dot h+H_2 \dot \alpha+H_3\alpha+H_4h\\ M_{se}=A_1\dot h+A_2 \dot \alpha+A_3\alpha+A_4h\\ \tag{10}$$ 在得到了自激力、断面位移以及速度时程信号之后，上式成为两个相互不耦合的线性方程组，求方程组的最小二乘解就可以得到全部气动导数。

上式也可以写成矩阵形式： $$\mathbf{Q_{se}}=\mathbf{\Phi}\mathbf{D} \tag{11}$$ 式中$\mathbf{Q_{se}}$为耦合自由度稳态简谐振动自激力矩阵。 $$\mathbf{Q_{se}}= \begin{bmatrix} L_{se1} & M_{se1}\\ L_{se2} & M_{se2}\\ \vdots & \vdots\\ L_{sen} & M_{sen}\\ \end{bmatrix} \tag{12}$$

$$\mathbf{\Phi}= \begin{bmatrix} \dot h_1 & \dot \alpha_1 & \alpha_1 & h_1\\ \dot h_2 & \dot \alpha_2 & \alpha_2 & h_2\\ \vdots & \vdots & \vdots & \vdots\\ \dot h_n & \dot \alpha_n & \alpha_n & h_n\\ \end{bmatrix} \tag{13}$$

$$\mathbf{D}= \begin{bmatrix} H_1 & A_1\\ H_2 & A_2\\ H_3 & A_3\\ H_4 & A_4 \end{bmatrix} \tag{14}$$

式中$n$为一个采样样本所包含的数点数。

利用最小二乘法可以得到$\mathbf{D}$的最小二乘估计值： $$\mathbf{D}=(\mathbf{\Phi}^T\mathbf{\Phi})^{-1}\mathbf{\Phi^T}\mathbf{Q_se} \tag{15}$$ 利用式(8)和式(9)可以得到无量纲化的颤振导数$H_\text{i}^*$、$A_\text{i}^*(i=1,2,3,4)$。

关于自激力的提取

时域方法与频域方法拥有相同的荷载假定：强迫振动模型在特定风速下受到的动态力除自激气动力外其他周期性荷载与零风速下保持一致。采用时域方法识别颤振导数时，传统的方法是将起始相位相同的特定风速下外力时程信号与零风速下外力时程信号相减得到自激力。即： $$L_{se}(t)=L_{dy,U}(t)-L_{dy,0}(t)\\ M_{se}(t)=M_{dy,U}(t)-M_{dy,0}(t)\\ \tag{16}$$ 对于试验而言，结构的振动并不是理想的稳态简谐振动，强迫振动装置本身的机械振动及外界振动都会作为噪声信号引入力信号及加速度信号中，这样无论是竖向、横向还是扭转振动，即使起始相位相同，各风速下加速度时程信号并非完全相等，惯性力(矩)等也不完全相等。因而采用起始相位相同的外力时程信号相减的方法获取自激力就会产生一定的误差。

而对于数值模拟，特定风速与零风速工况唯一差别即风速参数，其他参数可控制严格不变。因此时域传统的方法式$\eqref{eq:16}$能直接应用，得到自激力信号。

2、分状态多频强迫振动法原理

该方法结构做分状态强迫振动。不失一般性，选取结构做纯竖向运动的情况进行讨论。

1.频域识别方法基本原理

驱动结构做竖向小振幅的稳态多频简谐运动，即 $$\tag{17} \alpha (t)=0,p(t)\equiv0\\ h(t)=\sum_{k=1}^{n}h_k=h_0\sum_{k=1}^{n} e^{i(\omega_kt+\phi_h)}\\$$ 式中$h_0$为竖向振幅，$\omega _k$为竖向振动圆频率，$n$为考虑的频率个数，$\phi_h$为竖向振动初始相位。

在线性小振幅的情况下，自激气动力可以描述为： $$L_{se}(t)=\sum_{k=1}^{n} L_ke^{i(\omega_kt+\phi_{L_k})}\\ M_{se}(t)=\sum_{k=1}^{n}M_ke^{i(\omega_kt+\phi_{M_k})}\\ \tag{18}$$ 同时将式(17)代入式(2)得： $$L_{se}(t)=\rho U^2B\sum_{k=1}^{n}[K_{h_k}^2H_{1k}^*\frac{h_k}{B}i+K_{h_k}^2H_{4k}^*\frac{h_k}{B}]\\ M_{se}(t)=\rho U^2B^2\sum_{k=1}^{n}[K_{h_k}^2A_{1k}^*\frac{h_k}{B}i+K_{h_k}^2A_{4k}^*\frac{h_k}{B}] \tag{19}$$ 将式(18)代入式(19)，并对两端分别进行傅里叶变换。两端的复数力在强迫振动频率处相等，有： $$H_{1k}^*(K_{h_k})=\frac{1}{\rho B^2h_0\omega _k^2}\text{Im}[L_ke^{i(\phi_{L_k}-\phi_h)}]\\ H_{4k}^*(K_{h_k})=\frac{1}{\rho B^2h_0\omega _k^2}\text{Re}[L_ke^{i(\phi_{L_k}-\phi_h)}]\\ A_{1k}^*(K_{h_k})=\frac{1}{\rho B^3h_0\omega _k^2}\text{Im}[M_ke^{i(\phi_{M_k}-\phi_h)}]\\ H_{1k}^*(K_{h_k})=\frac{1}{\rho B^3h_0\omega _k^2}\text{Re}[M_ke^{i(\phi_{M_k}-\phi_h)}]\\$$ 同理，使结构做扭转小振幅的稳态多频简谐运动，可以得到： $$H_{2k}^*(K_{\alpha_k})=\frac{1}{\rho B^3\alpha_0\omega _k^2}\text{Im}[L_ke^{i(\phi_{L_k}-\phi_h)}]\\ H_{3k}^*(K_{\alpha_k})=\frac{1}{\rho B^3\alpha_0\omega _k^2}\text{Re}[L_ke^{i(\phi_{L_k}-\phi_h)}]\\ A_{2k}^*(K_{\alpha_k})=\frac{1}{\rho B^4\alpha_0\omega _k^2}\text{Im}[M_ke^{i(\phi_{M_k}-\phi_h)}]\\ H_{3k}^*(K_{\alpha_k})=\frac{1}{\rho B^4\alpha_0\omega _k^2}\text{Re}[M_ke^{i(\phi_{M_k}-\phi_h)}]\\$$

关于自激力的提取

同前述耦合状态强迫振动法频域识别方法。

注：采用分状态多频强迫振动方法时结构作强迫振动时的振幅不能过大，需要严格满足小振幅前提。一是小振幅情况下才能进行线性自激力的叠加；二是需要减少气动力倍频现象，影响气动力振幅的识别。

2.时域识别方法基本原理

利用式(8)和式(9)，式(2)变为： $$L_{se}=\sum_{k=1}^{m}[H_{1k}\dot h_k+H_{4k}h_k]\\ M_{se}=\sum_{k=1}^{m}[A_{1k}\dot h_k+A_{4k}h_k]\\ \tag{22}$$ 在得到了自激力、断面位移以及速度时程信号之后，上式成为两个相互不耦合的线性方程组，求方程组的最小二乘解就可以得到全部气动导数。

上式也可以写成矩阵形式： $$\mathbf{Q_{se}}=\mathbf{\Phi_h}\mathbf{D_h} \tag{23}$$ 式中$\mathbf{Q_{se}}$为单自由度竖向稳态简谐振动自激力矩阵。 $$\mathbf{Q_{se}}= \begin{bmatrix} L_{se1} & M_{se1}\\ L_{se2} & M_{se2}\\ \vdots & \vdots\\ L_{sen} & M_{sen}\\ \end{bmatrix} \tag{24}$$

$$\mathbf{\Phi_h}= \begin{bmatrix} \dot h_{11} & h_{11} & \dot h_{12} & h_{12} & \dots & \dot h_{1m} & h_{1m}\\ \dot h_{21} & h_{21} & \dot h_{22} & h_{22} & \dots & \dot h_{2m} & h_{2m}\\ \vdots & \vdots & \vdots & \vdots\\ \dot h_{n1} & h_{n1} & \dot h_{n2} & h_{n2} & \dots & \dot h_{nm} & h_{nm}\\ \end{bmatrix}_{n\times m} \tag{25}$$

$$\mathbf{D_h}= \begin{bmatrix} H_{11} & A_{11}\\ H_{41} & A_{41}\\ H_{12} & A_{12}\\ H_{42} & A_{42}\\ \vdots & \vdots\\ H_{1m} & A_{1m}\\ H_{4m} & A_{4m}\\ \end{bmatrix}_{m\times 2} \tag{26}$$

式中$n$为一个采样样本所包含的数点数，$m$为考虑的频率个数。

利用最小二乘法可以得到$\mathbf{D_h}$的最小二乘估计值： $$\mathbf{D_h}=(\mathbf{\Phi}_h^T\mathbf{\Phi_h})^{-1}\mathbf{\Phi_h^T}\mathbf{Q_se} \tag{27}$$ 同理可以得到$\mathbf{D_\alpha}$。利用式(8)和式(9)可以得到无量纲化的颤振导数$H_\text{i}^*$、$A_\text{i}^*(i=1,2,3,4)$。

关于自激力的提取

同前述耦合状态强迫振动法时域识别方法。

至此，提高计算速度的理论介绍完毕。