This commit is contained in:
yize 2025-06-03 11:00:42 +08:00
parent 8359d58723
commit 2ce0ee8043
2 changed files with 560 additions and 0 deletions

Binary file not shown.

After

Width:  |  Height:  |  Size: 88 KiB

View File

@ -0,0 +1,560 @@
<style>
p {
text-indent: 2em;
text-align: justify;
font-family: "SimSun", "Times New Roman", serif;
font-size: 12pt;
line-height: 1.5;
}
</style>
[TOC]
<div style="page-break-after: always;"></div>
# 1. 基于模态叠加法的叶片变形计算
## 1.1. 叶片模态计算原理
在本节会出现时间导数和距离导数函数,约定时间导数用变量正上方加点,距离导数变量右上方加撇,如$u(z,t)$为时间t和距离的函数$\dot{u}$表示时间导数,$u^{\prime}$表示距离导数。
### 1.1.1. 叶片模态计算
把叶片看成是一端固定于轮毂另一端自由的柔性悬臂梁, 基于模态叠加法,叶片侧向变形$u(z,t)$用前N阶模态表示为
$$
u(z,t)=\sum_{a=1}^{N}\phi_{a}(z)q_{a}(t) \tag{1.1}
$$
其中,$\phi_{a}(z)$为第a阶模态的振型函数是只关于的函数是柔性叶片径向距离$=$表示悬臂固定端, $=z$表示悬臂自由端。$q_{a}(t)$为第a阶模态广义坐标 是只关于时间的函数。对于悬臂梁,正交模态的广义坐标习惯上取自由端的变形量,每个模态可用柔性叶片长度无量纲化和归一化,故自由端的模态值为1。
Rayleigh-Ritz法$\phi_{a}(z)$可表示为
$$
\phi_a(z)=\sum_{b=P}^{N+p-1} D_{a b} \varphi_b(z) \quad(a=1,2, \cdots, N ; b=P, P+1, \cdots, N+P-1)\tag{1.2}
$$
式中$\varphi_b(z)$为形函数, $D_{a b}$为与第a个模态第b个形函数相关的系数。$u(z,t)$可表示为
$$
\begin{gathered}
u(z, t)=\sum_{a=1}^N \phi_a(z) q_a(t)=\sum_{a=1}^N\left[\sum_{b=P}^{N+P-1} D_{a b} \varphi_b(\mathrm{z})\right] q_a(t) \\
=\sum_{b=P}^{N+P-1} \varphi_b(\mathrm{z}) \sum_{a=1}^N D_{a b} q_a(t)=\sum_{b=P}^{N+P-1} \varphi_b(z) d_b(t)
\end{gathered}\tag{1.3}
$$
$$
d_b(t)=\sum_{a=1}^N D_{a b} q_a(t)\tag{1.4}
$$
式中,$d_b(t)$为形函数$\varphi_b(z)$的广义坐标。
假定悬臂梁的每个模态可用一个多项式表示,$\varphi_b(z)$可定义为
$$
\varphi_b(z)=\left(\frac{z}{Z}\right)^b \tag{1.5}
$$
在悬臂梁的固定端形函数的导数必为0 故有$P \geq 2$ 。本例在柔性叶片计算中取$P = 2$$N = 5$。叶片模态可表示为
$$
\phi_a(z)=\sum_{b=P}^{N+P-1} D_{a b} \left(\frac{z}{Z}\right)^b \tag{1.6}
$$
### 1.1.2. 相关系数$D_{ab}$计算
对于有N个自由度的保守定常系统运动方程为
$$
\sum_{b=P}^{N+P-1} m_{a b} \ddot{d}_b(t)+\sum_{b=P}^{N+P-1} k_{a b} d_b(t)=0 \tag{1.7}
$$
其中,$m_{a b}$、$k_{a b}$与广义坐标$d_b(t)$对应的广义质量、广义刚度。
当柔性梁以第$m
$阶模态振动,有
$$
q_a(t)= \begin{cases}Q_a \sin \left(\omega_a t+\beta_a\right) & (a=m) \\ 0 & (a \neq m)
\end{cases} \tag{1.8}
$$
所以
$$
d_b(t)=\sum_{a=1}^N D_{a, b} q_a(t)=D_{m, b} q_m(t)=D_{m, b} Q_m \sin \left(\omega_m t+\beta_m\right) \tag{1.9}
$$
1.9代入1.7),有
$$
\sum_{b=P}^{N+P-1}\left(-\omega^2 m_{m b}+k_{m b}\right) D_{m, b}=0 \tag{1.10}
$$
写出矩阵形式为
$$
\left(-\omega^2 M_d+K_d\right) D=0 \tag{1.11}
$$
其中,$M_d$$K_d$为与广义坐标$d_b(t)$对应的广义质量矩阵、广义刚度矩阵,$D$为模态系数矩阵$N*1$列向量。把1.11)改写为
$$
\left(\frac{1}{\omega^2} \boldsymbol{I}-\boldsymbol{K}_d^{-1} \boldsymbol{M}_d\right) \boldsymbol{D}=\mathbf{0} \tag{1.12}
$$
由方程式(4.12)可知有N个特征值$\frac{1}{\omega_a^2}$对应的N个特征向量$D_a$,即
$$
\boldsymbol{D}_a=\left(D_{a, p} D_{a, p+1} \cdots D_{a, N+p-1}\right)^T \tag{1.13}
$$
### 1.1.3. 广义质量矩阵$M_d$和广义刚度矩阵$K_d$计算
假设轮毂的旋转角速度为$\Omega$,且柔性叶片在挥舞方向和摆振方向的变形是相互独立的。不考虑叶片的锥角、叶尖质量。基于一维弹性体系统的模态叠加法,质量矩阵元素$m_{ij}$计算式为
$$
m_{i j}=\int_0^{R-R_H} \mu_B(z) \varphi_i(z) \varphi_j(z) d z \quad(i=1,2, \cdots, N ; j=1,2, \cdots, N) \tag{1.14}
$$
式中$\mu_B(z)$为叶片分布线密度,$R$表示风轮半径,$R_H$为轮毂半径。式(4.14)对挥舞和摆振方向都是有效的。
刚度矩阵元素$k_{ij}$计算式为
$$
k_{i j}=\int_0^{R-R_H} E I_B(z){\varphi}_i^{\prime \prime}(z) \varphi_j^{\prime \prime}(z) d z+\Omega^2 \int_0^{R-R_H}\left(\int_r^{R-R_H} \mu_B(\gamma)\left(R_H+\gamma\right) d \gamma\right) {\varphi}_i^{\prime}(z) \varphi_j^{\prime}(z) d z \tag{1.15}
$$
其中$E I_B(z)$为叶片的分布刚度,$r$为叶片长度的虚变量。式(4.15)中的第二项为离心刚度项。
至此, 通过式(4.14-4.15)求得广义质量矩阵$M_d$和广义刚度矩阵$K_d$后, 再通过式(4.12)求出模态系数$D_{ab}$,再根据式(4.6)可求出叶片的各阶模态。
## 1.2. 叶片弹性恢复力计算
柔性或弹性体的弯曲会产生弹性恢复力, 使其恢复到未变形位置。对于柔性叶片,与广义坐标$d_b(t)$对应的广义质量矩阵$M_q$和广义刚度矩阵$K_q$的计算式为
$$
m_{i j}=\int_0^{R-R_H} \mu_B(z) \phi_i(z) \phi_j(z) d z \tag{1.16}
$$
$$
k_{i j}=\int_0^{R-R_H} E I_B(z) \phi_i(z) \phi_j(z) d z+\Omega^2 \int_0^{R-R_H}\left(\int_r^{R-R_H} \mu_B(\gamma)\left(R_H+\gamma\right) d \gamma\right) \phi_i(z) \phi_j(z) d z \tag{1.17}
$$
式中,$i=1,2, \cdots, N$$j=1,2, \cdots, N$N等于某方向表示叶片变形的模态个数。在接下来的式(4.12)(4.20)中ij也如此。
对于大型风机叶片,由于转速相对较小, 通常可不考虑离心钢化作用。
叶片固有频率$f_j$公式为
$$
f_j=\frac{1}{2 \pi} \sqrt{\frac{k_{i j}}{m_{i j}}} \tag{1.18}
$$
其中下标j表示第j阶模态频率。
叶片阻尼矩阵$C_q$的元素$C_{ij}$计算为
$$
c_{i j}=\frac{\zeta_j k_{i j}}{\pi f_j} \tag{1.19}
$$
其中$\zeta_j$为第j阶模态阻尼
柔性叶片在变形方向的弹性恢复力计算表达式为
$$
{ }^{\text {Ela }} F_i^B=-\sum_j k_{i j} q_j(t)-\sum_j c_{i j} \dot{q}_j(t) \tag{1.20}
$$
## 1.3. 有结构扭角时叶片变形和转角计算
叶片在挥舞方向的变形用两个振动模态计算, 而叶片在摆振方向的柔性较小,故其变形用一个模态计算。设$\phi_{1F}(z)$、$\phi_{2F}(z)$为无扭角叶片的一
、二阶挥舞模态,$\phi_{1E}(z)$为无扭角叶片的一阶摆振模态,$q_{1}(t)$、$q_{2}(t)$为叶尖一、二阶挥舞变形, $q_{3}(t)$为叶尖一阶摆振变形, 为表示方便, 在后面的表达中省略模态函数中的z 和叶尖变形中的t。在叶片局部坐标系中 局部挥舞曲率可表示为$q_1 {\phi_{1\mathrm{~F}}}^{\prime \prime}+q_2 {\phi_{2 \mathrm{~F}}}^{\prime \prime}$,局部摆振曲率可表示为$q_3 {\phi_{1{\mathrm{E}}}}^{\prime \prime}$。假定局部扭角为$\theta_{\mathrm{tw}}(z)$ 为叶片结构预扭角和叶片桨矩角的和,则面外方向曲率为
$$
{u}^{\prime \prime}=\cos \theta_{\mathrm{tw}}\left[q_1 {\phi_{1\mathrm{~F}}}^{\prime \prime}+q_2 {\phi_{2 \mathrm{~F}}}^{\prime \prime}\right]+\sin \theta_{\mathrm{tw}} q_3 {\phi_{1{\mathrm{E}}}}^{\prime \prime} \tag{1.21}
$$
面内方向的曲率为
$$
v^{\prime \prime}=-\sin \theta_{\mathrm{tw}}\left[q_1 {\phi_{1\mathrm{~F}}}^{\prime \prime}+q_2 {\phi_{2 \mathrm{~F}}}^{\prime \prime}\right]+\cos \theta_{\mathrm{tw}} q_3 {\phi_{1{\mathrm{E}}}}^{\prime \prime} \tag{1.22}
$$
将式(1.21)和(1.22)表示为
$$
{u}^{\prime \prime}=q_1 {\varphi_1}^{\prime \prime}+q_2 {\varphi_2}^{\prime \prime}+q_3 {\varphi_3}^{\prime \prime} \tag{1.23}
$$
$$
v^{\prime \prime}=q_1 {\varphi_1}^{\prime \prime}+q_2 {\varphi_2}^{\prime \prime}+q_3
{\varphi_3}^{\prime \prime} \tag{1.24}
$$
定义${\varphi_1}^{\prime \prime} = \cos \theta_{\mathrm{tw}}{\phi_{1\mathrm{~F}}}^{\prime \prime}$${\varphi_2}^{\prime \prime} = \cos \theta_{\mathrm{tw}}{\phi_{2\mathrm{~F}}}^{\prime \prime}$${\varphi_3}^{\prime \prime} = \sin \theta_{\mathrm{tw}}{\phi_{1\mathrm{~E}}}^{\prime \prime}$为面外扭转形函数,$\psi_1^{\prime \prime}=-\sin \theta_{\mathrm{tw}}{\phi_{1\mathrm{~F}}}^{\prime \prime}$$\psi_2^{\prime \prime}=-\sin \theta_{\mathrm{tw}}{\phi_{2\mathrm{~F}}}^{\prime \prime}$$\psi_3^{\prime \prime}=\cos \theta_{\mathrm{tw}}{\phi_{1\mathrm{~E}}}^{\prime \prime}$为面内扭转形函数。对面外面内扭转形函数关于z进行两次积分得面外模态$\varphi_1$$\varphi_2$$\varphi_3$,和面内模态$\psi_1$、$\psi_2$、$\psi_3$。叶片变形统一在叶根坐标系$b_{m 0}\left(x_{b m 0} y_{b m 0} z_{b m 0}\right)$中表示, 在$x_{b m 0}$方向的变形(挥舞变形)$u$和$y_{b m 0}$方向的变形(摆振变形)$v$可表示为
$$
u=\sum_{i=1}^3 q_i \varphi_i \tag{1.25}
$$
$$
v=\sum_{i=1}^3 q_i \psi_i \tag{1.26}
$$
面外转角$\theta_{u}$和面内转角$\theta_{v}$可表示为
$$
\theta_u=\frac{\partial u}{\partial z}=\sum_{i=1}^3 q_i \varphi_i^{\prime} \tag{1.27}
$$
$$
\theta_v=\frac{\partial v}{\partial z}=\sum_{i=1}^3 q_i \psi_i^{\prime} \tag{1.28}
$$
$\theta_{u}$沿$y_{b m 0}$轴正向为正,$\theta_{v}$沿$x_{b m 0}$轴负向为正。
假定叶片弯曲变形后总的叶片长度不变,则在叶片锥坐标系$z_{mc}$轴方向的投影必减少, 减少量表达式为
$$
w=-0.5 \sum_{i=1}^3 \sum_{j=1}^3 s_{i j} q_i q_j \tag{1.29}
$$
式中有$s_{i j}=\int_0^z\left(\dot{\varphi}_i \dot{\varphi}_j^{\prime}+\psi_i^{\prime} \psi_j^{\prime}\right) d z$,为轴向缩减因子。
根据式(1.25)、(1.26)、(1.29),变形后的叶片节点当前位置向量为
$$
\boldsymbol{u}(z, t)=u i_{m c}+v j_{m c}+[r(z)+w] k_{m c} \tag{1.30}
$$
式中$r(z)$为叶片变形前从叶根到叶素节点的局部半径;$i_{m c}$、$j_{m c}$、$k_{m c}$为第个叶片锥坐标系各轴的基矢量。
# 2. 基于模态叠加法的塔架变形计算
## 2.1. 塔架变形计算
把塔架看成是倒立的悬臂梁, 且在自由端存在一个点质量$m_{Top}$,这个点质量即为塔顶质量,等于机舱、轮毂及风轮的质量和。假定塔架在前后向(fore-aft)和侧向(side-side)的变形是独立的,且塔架在这两个方向的刚度相等,故在两方向的固有频率和模态也是相同的。只要计算塔架一个方向的模态即可。
塔架与叶片的模态计算原理及弹性恢复力计算公式是一样的,只是塔架为自由端有质量的倒立悬臂梁,其广义质量矩阵和阻尼矩阵的元素计算稍有区别,为
$$
m_{i j}=m_{\mathrm{Top}}+\int_0^H \mu_{\mathrm{T}}(z) \varphi_i(z) \varphi_j(z) d z \tag{2.1}
$$
$$
k_{i j}=\int_0^H E I_{\mathrm{T}}(z) \varphi_i^{\prime \prime}(z) \varphi_j^{\prime \prime}(z) d z-g \int_0^H\left[m_{\mathrm{Top}}+\int_z^H \mu_{\mathrm{T}}(\gamma) d \gamma\right] \varphi_i^{\prime}(z) \varphi_j^{\prime}(z) d z \tag{2.2}
$$
上两式中的$\mu_{\mathrm{T}}(z)$为塔架的分布线密度,$E I_{\mathrm{T}}(z)$为塔架的分布刚度。从式(2.2)可知塔顶质量和塔架自身的重量减少了塔架的刚度。
塔架在前后方向和侧向的变形各用两个振动模态计算。设$\phi_{\mathrm{1FA}}(z)$、$\phi_{\mathrm{2FA}}(z)$为塔架前后方向的一、二阶模态,$\phi_{\mathrm{1SS}}(z)$、$\phi_{\mathrm{2SS}}(z)$为塔架在侧向的一、二阶模态,$q_{\mathrm{1FA}}(t)$、$q_{\mathrm{2FA}}(t)$为塔顶一、二阶前后向变形,$q_{\mathrm{1SS}}(t)$、$q_{\mathrm{2SS}}(t)$为塔顶一、二阶侧向变形,为表示方便, 在后面的表达中省略模态函数中的z和塔顶变形中的t。在塔基坐标系$t_o(x_{t0}, y_{t0}, z_{t0})$中,塔架在任意高度位置的前后向变形$u_{\mathrm{FA}}$和侧向变形$v_{\mathrm{SS}}$表示为
$$
\left\{\begin{array}{c}
u_{\mathrm{FA}}=\phi_{\mathrm{1FA}} q_{\mathrm{1FA}}+\phi_{ \mathrm{2FA}} q_{ \mathrm{2FA}} \\
v_{\mathrm{SS}}=\phi_{\mathrm{1SS}} q_{\mathrm{1SS}}+\phi_{\mathrm{2SS}} q_{\mathrm{2SS}}
\end{array}\right. \tag{2.3}
$$
假定塔架在$z_{t0}$轴向无转角,有$y_{t0}$轴向转角$\theta_{\mathrm{FA}}$和$x_{t0}$轴向转角$\theta_{\mathrm{SS}}$为
$$
\left\{\begin{array}{l}
\theta_{\mathrm{FA}}=\frac{\partial u_{\mathrm{FA}}}{\partial z} \\
\theta_{\mathrm{SS}}=\frac{\partial v_{\mathrm{SS}}}{\partial z}
\end{array}\right. \tag{2.4}
$$
塔架的轴向缩减为
$$
\left\{\begin{array}{l}
w_{\mathrm{T}}=-0.5\left(\sum_{i=1}^2 \sum_{j=1}^2 s_{i j \mathrm{FA}} q_{\mathrm{iFA}} q_{j \mathrm{FA}}+\sum_{i=1}^2 \sum_{j=1}^2 s_{i j \mathrm{SS}} q_{i \mathrm{SS}} q_{j \mathrm{SS}}\right) \\
s_{i j \mathrm{FA}}=\int_0^z \phi_{\mathrm{iFA}}^{\prime} \phi_{j \mathrm{FA}}^{\prime} d z \quad s_{i j \mathrm{SS}}=\int_0^z \phi_{i \mathrm{SS}}^{\prime} \phi_{j \mathrm{SS}}^{\prime} d z
\end{array}\right. \tag{2.5}
$$
根据(2.3)和(2.5),变形后塔架任意节点当前位置向量为
$$
\boldsymbol{u}_{\mathrm{T}}(z, t)=u_{\mathrm{FA}} \boldsymbol{i}_{t 0}+v_{\mathrm{SS}} \boldsymbol{j}_{t 0}+\left[Z(z)+w_T\right] \boldsymbol{k}_{t 0} \tag{2.6}
$$
式中$Z(z)$为塔架变形前从塔基到塔架节点的距离;$i_{t0}$、$j_{t0}$、$k_{t0}$为塔基坐标系各轴的基矢量。
# 3. 整机刚柔混合多体动力学建模
在多体动力学建模过程中,把水平轴三叶片上风向海上浮式风电机抽象成由浮式平台、塔架、机舱、低速轴、轮毂、叶片、高速轴及电机构成的多体系统。其中浮式平台、机舱、高速轴及电机看成是刚性的,塔架、低速轴及叶片看成柔性的。塔架和叶片的变形用模态叠加法表示,低速轴的柔性用旋转惯量- 刚度阻尼器表示。用Kane方法建立整机结构动力学模型时把叶片离散成17个叶素单元体 把塔架离散成11个单元体。
## 3.1. 定义参考系
对于结构复杂的多体系统为方便描述建立各组成件的连体参考系并确定各连体参考系间的转换关系。海上浮式风电机坐标系如图3.1所示, 图中$O(XYZ)$为惯性参考系,$p(x_{p}, y_{p}, z_{p})$、$t_i(x_{ti}, y_{ti}, z_{ti})$、$t_p(x_{top}, y_{top}, z_{top})$、$n(x_{n}, y_{n}, z_{n})$、$l(x_{l}, y_{l}, z_{l})$、$h(x_{h}, y_{h}, z_{h})$、$c(x_{mc}, y_{mc}, z_{mc})$、$b_{m0}(x_{bm0}, y_{bm0}, z_{bm0})$及$b_{mi}(x_{bmi}, y_{bmi}, z_{bmi})$分别表示浮式平台参考点$p$、塔架单元节点$t_i$、塔架顶点$t_p$、机舱参考点$n$、低速轴参考点$l$、轮毂参考点$h$、风轮锥顶点$c$、第$m$个叶片叶根节点$b_{m0}$及叶片离散单元(叶素)节点$b_{mi}$处点连体参考系。图中的$d$、$f$、$e$点分别表示平台、机舱及轮毂的质心。$\theta_l$为风轮轴倾斜角,$\theta_c$为叶片锥角。假定惯性参考系坐标轴的基矢量分别为$X=(1,0,0)^T$、$Y=(0,1,0)^T$、$Z=(0,0,1)^T$,根据各参考系间的变换矩阵可求得任意时刻各连体系各坐标轴的单位矢量。
<figure style="text-align: center;">
<img src="/坐标系示意图.png" alt="图片描述">
<figcaption style="text-align: center;">图 1: 坐标系示意图</figcaption>
</figure>
假定浮式平台相对固定参考系的转动为小角转动,塔架和叶片因变形产生的转动也为小角转动,故从惯性系到平台连体系、从塔基系到塔架单元局部系及从叶根系到叶素局部系间的变换可以不考虑转动角顺序,采用小角变换矩阵:
$$
T_{S A T}\left(\theta_1, \theta_2, \theta_3\right)=\left(\begin{array}{ccc}
\frac{\theta_1^2 b+\theta_2^2+\theta_3^2}{a b} & \frac{\theta_3 a+\theta_1 \theta_2(b-1)}{a b} & \frac{-\theta_2 a+\theta_1 \theta_3(b-1)}{a b} \\
\frac{-\theta_3 a+\theta_1 \theta_2(b-1)}{a b} & \frac{\theta_1^2+\theta_2^2 b+\theta_3^2}{a b} & \frac{\theta_1 a+\theta_2 \theta_3(b-1)}{\left(\theta_1^2+\theta_2^2+\theta_3^2\right) \sqrt{1+\theta_1^2+\theta_2^2+\theta_3^2}} \\
\frac{\theta_2 a+\theta_1 \theta_3(b-1)}{a b} & \frac{\theta_1 a+\theta_2 \theta_3(b-1)}{a b} & \frac{\theta_1^2+\theta_2^2+\theta_3^2 b}{a b}
\end{array}\right) \tag{3.1}
$$
式中$a = \theta_1^2+\theta_2^2+\theta_3^2$$b = \sqrt{1+\theta_1^2+\theta_2^2+\theta_3^2}$$\theta_1$、$\theta_2$、$\theta_3$为绕三个坐标轴的转动角。
根据惯性系到平台参考系的小角变换矩阵有
$$
\left(\begin{array}{l}
x_{\mathrm{p}} \\
y_{\mathrm{p}} \\
z_{\mathrm{p}}
\end{array}\right)=T_{\mathrm{SAT}}\left(\theta_X, \theta_Y, \theta_Z\right)\left(\begin{array}{l}
X\\
Y \\
Z
\end{array}\right) \tag{3.2}
$$
式中$\theta_X, \theta_Y, \theta_Z$为浮式平台横摇角、纵摇角及艏摇角。$x_p, y_p, z_p$为平台系三个坐标轴的基矢量。
从平台参考系$p(x_{p}, y_{p}, z_{p})$通过小角变换到塔顶系$t_p(x_{top}, y_{top}, z_{top})$矢量变换式为
$$
\left(\begin{array}{l}
x_{\text {top }} \\
y_{\text {top }} \\
z_{\text {top }}
\end{array}\right)=T_{\mathrm{SAT}}\left(\theta_{x_{\mathrm{p}}}, \theta_{y_{\mathrm{p}}}, \theta_{z_{\mathrm{p}}}\right)\left(\begin{array}{l}
x_{\mathrm{p}} \\
y_{\mathrm{p}} \\
z_{\mathrm{p}}
\end{array}\right) \tag{3.3}
$$
式中$\theta_{x_{\mathrm{p}}} = -\theta_{SS}$$\theta_{y_{\mathrm{p}}} = \theta_{FA}$$\theta_{z_{\mathrm{p}}} = 0$。
从塔顶系$t_p(x_{top}, y_{top}, z_{top})$的$z_{top}$正向旋转偏航角$\theta_y$到机舱系$n(x_{n}, y_{n}, z_{n})$,机舱系的基矢量
$$
\left(\begin{array}{l}
x_{\mathrm{n}} \\
y_{\mathrm{n}} \\
z_{\mathrm{n}}
\end{array}\right)=\left(\begin{array}{ccc}
\cos \theta_y & \sin \theta_y & 0 \\
-\sin \theta_y & \cos \theta_y & 0 \\
0 & 0 & 1
\end{array}\right)\left(\begin{array}{l}
x_{\text {top }} \\
y_{\text {top }} \\
z_{\text {top }}
\end{array}\right) \tag{3.4}
$$
从机舱系$n(x_{n}, y_{n}, z_{n})$的$y_{n}$轴正向旋转轴倾斜角$\theta_l$到低速轴系$l(x_{l}, y_{l}, z_{l})$,低速轴系的基矢量为
$$
\left(\begin{array}{l}
x_1 \\
y_1 \\
z_l
\end{array}\right)=\left(\begin{array}{ccc}
\cos \theta_l & 0 & -\sin \theta_l \\
0 & 1 & 0 \\
\sin \theta_l & 0 & \cos \theta_l
\end{array}\right)\left(\begin{array}{l}
x_n \\
y_{\mathrm{n}} \\
z_{\mathrm{n}}
\end{array}\right)\tag{3.5}
$$
从低速轴系$l(x_{l}, y_{l}, z_{l})$的$x_l$轴正向旋转风轮方位角$\theta_a$到轮毂系$h(x_{h}, y_{h}, z_{h})$,轮毂系随风轮转动,其$z_{h}$轴与叶片1方位角一致
$$
\left(\begin{array}{l}
x_{\mathrm{h}} \\
y_{\mathrm{h}} \\
z_{\mathrm{h}}
\end{array}\right)=\left(\begin{array}{ccc}
1 & 0 & 0 \\
\cos \theta_a & \sin \theta_a & 0 \\
-\sin \theta_a & \cos \theta_a & 0
\end{array}\right)\left(\begin{array}{l}
x_1 \\
y_1 \\
z_{\mathrm{l}}
\end{array}\right) \tag{3.6}
$$
绕轮毂系$h(x_{h}, y_{h}, z_{h})$的$x_{h}$轴正向旋转角度$2 \pi(m-1) / 3$后,再绕$y_{h}$轴负向旋转叶片锥角$\theta_c$得到第$m$个叶片的锥坐标系$c(x_{mc}, y_{mc}, z_{mc})$,通常三个叶片的锥顶点$c$、与轮毂参考点重合,有
$$
\left(\begin{array}{l}
x_{\mathrm{cm}} \\
\boldsymbol{y}_{\mathrm{cm}} \\
z_{\mathrm{cm}}
\end{array}\right)=\left(\begin{array}{ccc}
\cos \theta_c & 0 & \sin \theta_c \\
0 & 1 & 0 \\
-\sin \theta_c & 0 & \cos \theta_c
\end{array}\right)\left(\begin{array}{ccc}
1 & 0 & 0 \\
\cos [2 \pi(m-1) / 3] & \sin [2 \pi(m-1) / 3] & 0 \\
-\sin [2 \pi(m-1) / 3] & \cos [2 \pi(m-1) / 3] & 0
\end{array}\right)\left(\begin{array}{l}
x_{\mathrm{h}} \\
y_{\mathrm{h}} \\
z_{\mathrm{h}}
\end{array}\right) \tag{3.7}
$$
式中的$m=1,2,3$,指代三个叶片。绕锥坐标系$c(x_{mc}, y_{mc}, z_{mc})$的$z_{mc}$轴负向旋转桨矩角$\theta_p$得到叶根系$b_{m0}(x_{m0}, y_{m0}, z_{m0})$,有
$$
\left(\begin{array}{l}
x_{\mathrm{bm} 0} \\
y_{\mathrm{bm} 0} \\
z_{\mathrm{bm} 0}
\end{array}\right)=\left(\begin{array}{ccc}
\cos \theta_p & -\sin \theta_p & 0 \\
\sin \theta_p & \cos \theta_p & 0 \\
0 & 0 & 1
\end{array}\right)\left(\begin{array}{l}
x_{\mathrm{cm}} \\
y_{\mathrm{cm}} \\
z_{\mathrm{cm}}
\end{array}\right) \tag{3.8}
$$
从叶根系$b_{m0}(x_{m0}, y_{m0}, z_{m0})$的$z_{m0}$轴负向旋转叶素局部扭角$\theta_t$得到叶素局部系$b_{mi}(x_{mbi}, y_{mbi}, z_{mbi})$,有
$$
\left(\begin{array}{c}
x_{\mathrm{bmi}} \\
y_{\mathrm{bmi}} \\
z_{\mathrm{bmi}}
\end{array}\right)=\left(\begin{array}{ccc}
\cos \theta_t & -\sin \theta_t & 0 \\
\sin \theta_t & \cos \theta_t & 0 \\
0 & 0 & 1
\end{array}\right)\left(\begin{array}{l}
x_{\mathrm{bm} 0} \\
y_{\mathrm{bm} 0} \\
z_{\mathrm{bm} 0}
\end{array}\right) \tag{3.9}
$$
## 3.2. 定义整机运动自由度
对三叶片水平轴上风向海上浮式风电机共设立22个自由度浮式平台设置6个自由度三个平移自由度$q_1$$q_2$$q_3$,三个旋转自由度$q_4$$q_5$$q_6$。塔架4个自由度分别是塔顶一阶、二阶前后向(fore-aft)变形,用$q_7$$q_9$表示,塔顶一阶、二阶侧向(side-side)变形,用$q_8$$q_{10}$表示。一个机舱偏航自由度$q_{11}$,一个低速轴转动自由度$q_{12}$,一个低速轴扭转柔性自由度$q_{13}$。每一叶片三个自由度,叶尖一阶、二阶挥舞(flapwise)变形$q_{14+3(m-1)}$、$q_{16+3(m-1)}$和一阶摆振(edgewise)变形$q_{15+3(m-1)}$下标中,${m =1,2,3}$表示不同的叶片。
## 3.3. 基于Kane法的整机多体结构动力学建模
对同一多体系统广义速率有不同的选取方案, 系统中同一质量和同一刚体也将有不同的偏速度和偏角速度。故在选取广义速率时应尽量使对应偏速度和偏角速度的表达式简单,以便于求导运算,易于建立动力学方程和使方程间的耦合程度减弱。选取广义速率$u_r=\dot{q}_r, r=1,2, \cdots, 22$。$\dot{q}_r$为各自由度速度。在论文中用$\omega_X^{(r)}$表示刚体$X$对应第$r$个广义速率的偏角速度,在海上浮式风电机中刚体$X$指代平台$P$、塔架单元$T_i$,塔顶基板$T_{tp}$、机舱$N$、低速轴$L$,电机$G$、轮毂$H$,叶素单元$B_{mi}$。用$v_X^{(r)}$表示质点或参考点$x$的第$r$个广义速率的偏线性速度,同样$x$表示平台参考点$p$;平台质心点$d$、塔架单元节点$t_i$、塔顶点$t_p$、机舱参考点$n$、机舱质心点$f$、轮毂参考点$h$、轮毂质心点$e$及叶素单元节点$b_{mi}$。符号$\overrightarrow{r x y}$指在惯性系中,从点$x$到点$y$的空间位置向量。
### 3.3.1. 各刚体偏角速度
浮式$p$的偏角速度
$$
\boldsymbol{\omega}_{\mathrm{P}}^{(1)}=0 \quad \boldsymbol{\omega}_{\mathrm{P}}^{(2)}=0 \quad \boldsymbol{\omega}_{\mathrm{P}}^{(3)}=0 \quad \boldsymbol{\omega}_{\mathrm{P}}^{(4)}=\boldsymbol{X} \quad \boldsymbol{\omega}_{\mathrm{P}}^{(5)}=\boldsymbol{Y} \quad \boldsymbol{\omega}_{\mathrm{P}}^{(6)}=\boldsymbol{Z} \quad \boldsymbol{\omega}_{\mathrm{P}}^{(r)}=0(r=7,8, \cdots, 22) \tag{3.10}
$$
塔顶$T_{tp}$的偏角速度
$$
\begin{aligned}
& \boldsymbol{\omega}_{\mathrm{Ttp}}^{(1)}=0 \quad \boldsymbol{\omega}_{\mathrm{Ttp}}^{(2)}=0 \quad \omega_{\mathrm{Ttp}}^{(3)}=0 \quad \omega_{\mathrm{Ttp}}^{(4)}=X \quad \omega_{\mathrm{Ttp}}^{(5)}=Y \quad \omega_{\mathrm{Ttp}}^{(6)}=Z \quad \omega_{\mathrm{Ttp}}^{(7)}=\phi_{1 F A}^{tt p} y_p \\
& \boldsymbol{\omega}_{\mathrm{Ttp}}^{(8)}=-\phi_{1 S S}^{ttp} \boldsymbol{x}_p \quad \boldsymbol{\omega}_{\mathrm{Ttp}}^{(9)}=\phi_{2 \mathrm{FA}}^{\mathrm{ttp}} \boldsymbol{y}_p \quad \boldsymbol{\omega}_{\mathrm{Ttp}}^{(10)}=-\phi_{2 S S}^{\mathrm{ttp}} \boldsymbol{x}_p \quad \boldsymbol{\omega}_{\mathrm{Ttp}}^{(r)}=0(r=11,12, \cdots, 22)
\end{aligned} \tag{3.11}
$$
其中$\phi_{1 F A}^{tt p}$、$\phi_{2 F A}^{tt p}$为塔顶前后向一、二阶模态,$\phi_{1 SS}^{tt p}$、$\phi_{2SS}^{tt p}$为塔顶侧向一、二阶模态。
机舱$N$的偏角速度
$$
\omega_{\mathrm{N}}^{(r)}=\omega_{\mathrm{T}_{\mathrm{tp}}}^{(r)}(r=1,2, \cdots, 10) \quad \omega_{\mathrm{N}}^{(11)}=z_n \quad \omega_{\mathrm{N}}^{(r)}=0(r=11,12, \cdots, 22) \tag{3.12}
$$
低速轴$L$的偏角速度
$$
\omega_{\mathrm{N}}^{(r)}=\omega_{\mathrm{N}}^{(r)}(r=1,2, \cdots, 11) \quad \omega_{\mathrm{H}}^{(12)}=\boldsymbol{x}_l \quad \omega_{\mathrm{H}}^{(13)}=\boldsymbol{x}_l \quad \omega_{\mathrm{H}}^{(r)}=0(r=14,15, \cdots, 22) \tag{3.13}
$$
电机$G$的偏角速度
$$
\omega_{\mathrm{G}}^{(r)}=\omega_{\mathrm{N}}^{(r)}(r=1,2, \cdots, 11) \quad \omega_{\mathrm{G}}^{(12)}=\lambda x_l \quad \omega_{\mathrm{G}}^{(r)}=0(r=13,14, \cdots, 22) \tag{3.14}
$$
式中$\lambda$为传动比。
叶素单元${B_{mi}}$的偏角速度
$$
\begin{aligned}
& \boldsymbol{\omega}_{\mathrm{B} m i}^{(r)}=\omega_{\mathrm{H}}^{(r)}(r=1,2, \cdots, 13) \quad \omega_{\mathrm{B} m i}^{(14+3 m-3)}=\left(\varphi_l^{b m i}\right)^{\prime} y_{\mathrm{b} m 0}-\left(\psi_l^{b m i}\right)^{\prime} x_{\mathrm{b} m 0} \\
& \boldsymbol{\omega}_{\mathrm{B} m i}^{(15+3 m-3)}=\left(\varphi_2^{\mathrm{bmi}}\right)^{\prime} y_{\mathrm{b} m 0}-\left(\psi_2^{\mathrm{bmi}}\right)^{\prime} x_{\mathrm{b} m 0} \quad \boldsymbol{\sigma}_{\mathrm{B} m i}^{(16+3 m-3)}=\left(\varphi_3^{\mathrm{bmi}}\right)^{\prime} \boldsymbol{y}_{\mathrm{b} m 0}-\left(\psi_3^{\mathrm{b} m i}\right)^{\prime} x_{\mathrm{b} m 0}
\end{aligned} \tag{3.15}
$$
### 3.3.2. 质心或参考点偏线速度
平台参考点$p$的偏线速度为
$$
\boldsymbol{v}_{\mathrm{p}}^{(1)}=\boldsymbol{X} \quad \boldsymbol{v}_{\mathrm{p}}^{(2)}=\boldsymbol{Y} \quad \boldsymbol{v}_{\mathrm{p}}^{(3)}=\boldsymbol{Z} \quad \boldsymbol{v}_{\mathrm{p}}^{(r)}=0(r=4,5, \cdots, 22) \tag{3.16}
$$
平台质心$d$的偏线速度为
$$
\begin{gathered}
v_{\mathrm{d}}^{(1)}=v_{\mathrm{p}}^{(1)}(r=1,2,3) \quad v_{\mathrm{d}}^{(4)}=\omega_{\mathrm{p}}^{(4)} \times \overrightarrow{r p d} \quad v_{\mathrm{d}}^{(5)}=\omega_{\mathrm{P}}^{(5)} \times \overrightarrow{r p d} \\
v_{\mathrm{d}}^{(6)}=\omega_{\mathrm{P}}^{(6)} \times \overrightarrow{r p d} \quad v_{\mathrm{d}}^{(r)}=0(r=7,8, \cdots, 22) ;
\end{gathered} \tag{3.17}
$$
式中$\overrightarrow{r p d}$为平台质心$d$到参考点$p$的向量。
塔架顶点$t_{tp}$的偏线速度为
$$
\begin{aligned}
& v_{\mathrm{ttp}}^{(r)}=v_{\mathrm{p}}^{(r)}(r=1,2,3) \quad v_{\mathrm{tpp}}^{(4)}=\omega_{\mathrm{P}}^{(4)} \times \overrightarrow{r p t_{t p}} \quad v_{\mathrm{ttp}}^{(5)}=\omega_{\mathrm{P}}^{(5)} \times \overrightarrow{r p t_{t p}} \quad v_{\mathrm{tpp}}^{(6)}=\omega_{\mathrm{P}}^{(6)} \times \overrightarrow{r p t_{t p}} \\
& v_{\mathrm{ttp}}^{(7)}=x_{\mathrm{p}}-\left(S_{\mathrm{FA} 11}^{\mathrm{ttp}} q_7+S_{\mathrm{FA} 12}^{\mathrm{ttp}} q_9\right) z_{\mathrm{p}} \quad v_{\mathrm{ttp}}^{(8)}=x_{\mathrm{p}}-\left(S_{\mathrm{SS} 11}^{\mathrm{ttp}} q_8+S_{\mathrm{SSt} 12}^{\mathrm{ttp}} q_9\right) z_{\mathrm{p}} \\
& v_{\mathrm{ttp}}^{(9)}=x_{\mathrm{p}}-\left(S_{\mathrm{FA} 12}^{\mathrm{ttp}} q_7+S_{\mathrm{FA} 22}^{\mathrm{ttp}} q_9\right) z_{\mathrm{p}} \quad v_{\mathrm{ttp}}^{(10)}=x_{\mathrm{p}}-\left(S_{\mathrm{SS} 12}^{\operatorname{ttp}} q_8+S_{\mathrm{SS} 22}^{\operatorname{ttp}}q_{10}\right) z_{\mathrm{p}} \\
& v_{\text {ttp }}^{(r)}=0(r=11,12, \cdots, 22)
\end{aligned} \tag{3.18}
$$
式中$S_{\mathrm{FA} 11}^{\mathrm{ttp}}$、$S_{\mathrm{FA} 12}^{\mathrm{ttp}}$、$S_{\mathrm{FA} 22}^{\mathrm{ttp}}$为塔顶前后向缩减因子,$S_{\mathrm{SS} 11}^{\mathrm{ttp}}$、$S_{\mathrm{SS} 12}^{\mathrm{ttp}}$、$S_{\mathrm{SS} 22}^{\mathrm{ttp}}$为塔顶侧向缩减因子;$\overrightarrow{r p t_{t p}}$为参考点$p$到塔架顶点$t_{tp}$的位置向量。塔架单元节点$t_i$的偏线速度计算与塔顶点的相似,只是位置向量和缩减因子不同。
机舱质心$f$的偏线速度
$$
\begin{aligned}
& v_f^{(r)}=v_{\mathrm{tpp}}^{(r)}(r=1,2,3) \quad v_f^{(r)}=v_{\mathrm{ttp}}^{(r)}+\omega_{\mathrm{N}}^{(r)} \times \overrightarrow{r o f}(r=4,5, \cdots, 10) \\
& \boldsymbol{v}_f^{(11)}=\omega_{\mathrm{N}}^{(11)} \times \overrightarrow{r o f} \quad v_f^{(r)}=0(r=12,13, \cdots, 22)
\end{aligned} \tag{3.19}
$$
式中$\overrightarrow{r o f}$为参考点$o$到机舱质心$f$的位置向量。
轮毂参考点$h$的偏线速度
$$
\begin{aligned}
v_{\mathrm{h}}^{(r)}=v_{\text {ttp }}^{(r)}(r & =1,2,3) \quad v_{\mathrm{h}}^{(r)}=v_{\mathrm{ttp}}^{(r)}+\omega_{\mathrm{N}}^{(r)} \times \overrightarrow{r o h}(r=4,5, \cdots, 10) \\
v_{\mathrm{h}}^{(11)} & =\boldsymbol{\omega}_{\mathrm{N}}^{(11)} \times \overrightarrow{r o h} \quad v_{\mathrm{h}}^{(r)}=0(r=12,13, \cdots, 22)
\end{aligned} \tag{3.20}
$$
式中$\overrightarrow{r o h}$为参考点$o$到轮毂参考点$h$的位置向量。
轮毂质心$e$的偏线速度
$$
\begin{aligned}
v_{\mathrm{e}}^{(r)}=v_{\mathrm{h}}^{(r)}(r=1,2,3) \quad \nu_{\mathrm{e}}^{(r)}=\boldsymbol{v}_{\mathrm{h}}^{(r)}+\omega_{\mathrm{H}}^{(r)} \times \overline{r h e}(r=4,5, \cdots, 11) \\
\boldsymbol{v}_{\mathrm{e}}^{(12)}=\boldsymbol{\omega}_{\mathrm{H}}^{(12)} \times \overrightarrow{r h e} \quad \boldsymbol{v}_{\mathrm{e}}^{(13)}=\boldsymbol{\omega}_{\mathrm{H}}^{(13)} \times \overrightarrow{r h e} \quad \boldsymbol{v}_{\mathrm{h}}^{(r)}=0(r=14,15, \cdots, 22)
\end{aligned} \tag{3.21}
$$
式中$\overrightarrow{r h e}$为参考点$h$到轮毂质心$e$的位置向量。
叶素节点$b_{mi}$的偏线速度
$$
\begin{aligned}
& v_{\mathrm{b} m i}^{(r)}=v_{\mathrm{h}}^{(r)}(r=1,2,3) \quad v_{\mathrm{bmi}}^{(r)}=v_{\mathrm{h}}^{(r)}+\omega_{\mathrm{H}}^{(r)} \times \overrightarrow{r h b_{m i}}(r=4,5, \cdots, 11) \\
& \boldsymbol{v}_{\mathrm{bmi}}^{(\mathrm{I2)}}=\boldsymbol{\omega}_{\mathrm{H}}^{(12)} \times \overrightarrow{r h b_{m i}} \quad \boldsymbol{v}_{\mathrm{b} m i}^{(13)}=\boldsymbol{\omega}_{\mathrm{H}}^{(13)} \times \overline{r h b_{m i}} \\
& \boldsymbol{v}_{\mathrm{b} m i}^{(14+3 m-3)}=\varphi_1^{\mathrm{b} m i} \boldsymbol{x}_{\mathrm{b} m 0}+\psi_1^{\mathrm{b} m i} \boldsymbol{y}_{\mathrm{b} m 0}-\left(S_{11}^{\mathrm{bmi}} q_{14+3 m-3}+S_{12}^{\mathrm{b} m i} q_{15+3 m-3}+S_{13}^{\mathrm{b} m i} q_{16+3 m-3}\right) z_{\mathrm{b} m 0} \\
& \boldsymbol{v}_{\mathrm{b} m i}^{(15+3 m-3)}=\varphi_3^{\mathrm{b} m i} \boldsymbol{x}_{\mathrm{b} m 0}+\psi_3^{\mathrm{bmi}} \boldsymbol{y}_{\mathrm{b} m 0}-\left(S_{31}^{\mathrm{bmi}} q_{14+3 m-3}+S_{32}^{\mathrm{bmi}} q_{15+3 m-3}+S_{33}^{\mathrm{bmi}} q_{16+3 m-3}\right) z_{\mathrm{b} m 0} \\
& v_{\mathrm{b} m i}^{(16+3 m-3)}=\varphi_2^{\mathrm{b} m i} x_{\mathrm{b} m 0}+\psi_2^{\mathrm{bmi}} y_{\mathrm{b} m 0}-\left(S_{21}^{\mathrm{b} m i} q_{14+3 m-3}+S_{22}^{\mathrm{bmi}} q_{15+3 m-3}+S_{23}^{\mathrm{b} m i} q_{16+3 m-3}\right) z_{\mathrm{b} m 0}
\end{aligned} \tag{3.22}
$$
其中$\varphi_1^{\mathrm{b} m i} $、$\varphi_2^{\mathrm{b} m i} $、$\varphi_3^{\mathrm{b} m i} $表示第$m$个叶片第$i$个叶素节点处挥舞方向的扭转形函数,$\psi_1^{\mathrm{bmi}}$、$\psi_2^{\mathrm{bmi}}$、$\psi_3^{\mathrm{bmi}}$表示摆振方向的扭转形函数,$S_{11}^{\mathrm{b} m i}$、$S_{12}^{\mathrm{b} m i}$ 等为轴向缩减因子
己知偏角速度后,可用它们表示刚体的角速度和角加速度为
$$
\omega_{\mathrm{X}}=\sum_{r=1}^{22} \omega_{\mathrm{X}}^{(r)} \dot{q}_r \tag{3.23}
$$
$$
\boldsymbol{\beta}_{\mathrm{X}}=\dot{\boldsymbol{\omega}}_{\mathrm{X}}=\sum_{r=1}^{22} \boldsymbol{\omega}_{\mathrm{X}}^{(r)} \ddot{q}_r+\sum_{r=1}^{22} \dot{\omega}_{\mathrm{X}}^{(r)} \dot{q}_r \tag{3.24}
$$
己知偏线速度后,各质心或参考点的线速度和线加速度为
$$
v_x=\sum_{r=1}^{22} v_x^{(r)} \dot{q}_r \tag{3.25}
$$
$$
\alpha_x=\dot{v}_x=\sum_{r=1}^{22} v_x^{(r)} \ddot{q}_r+\sum_{r=1}^{22} \dot{v}_x^{(r)} \dot{q}_r \tag{3.26}
$$
由式(3.24)和(3.26)可知角加速度和线加速度由两部分组成,一部分与自由度的加速度有关,另一部分与自由度的加速度无关。
### 3.3.3. 整机动力学方程
海上浮式风电机各组成件受到的惯性力和力矩主要包括:浮式平台受到的惯性力$F_d$和惯性力矩$L_p$,塔架离散单元惯性力$F_{ti}$;机舱平动惯性力$F_f$和和偏航惯性力矩$L_{N}$,电机转动惯性力矩$L_{G}$;轮毂惯性力$F_H$和惯性力矩$L_{H}$,叶片离散单元惯性力$F_{Bmi}$。各刚体的惯性力统一表示为
$$
F_x=m_x a_x=m_x\left(\sum_{r=1}^{22} v_x^{(r)} \ddot{q}_r+\sum_{r=1}^{22} \dot{v}_x^{(r)} \dot{q}_r\right) \tag{3.27}
$$
式中的下标$x$指代浮式平台质心$d$、塔架单元节点$t_i$及机舱质心$f$、轮毂质心$e$及叶素节点$b_{mi}$$m_x$为质心或节点处的质量。根据动量矩定理,各惯性力矩统一表示为
$$
\begin{gathered}
\boldsymbol{L}_{\mathrm{Xx}}=\boldsymbol{\omega}_{\mathrm{X}} \times\left[I_{\mathrm{Xx}} \boldsymbol{x}\left(\boldsymbol{x} \cdot \boldsymbol{\omega}_{\mathrm{X}}\right)\right]+I_{\mathrm{Xx}} x\left(x \cdot \boldsymbol{\beta}_{\mathrm{X}}\right) \\
=\boldsymbol{\omega}_{\mathrm{X}} \times\left[I_{\mathrm{Xx}} x\left(x \cdot \omega_{\mathrm{X}}\right)\right]+I_{\mathrm{Xx}} x\left[x \cdot\left(\sum_{r=1}^{22} \boldsymbol{\omega}_{\mathrm{X}}^{(r)} \ddot{q}_r+\sum_{r=1}^{22} \dot{\omega}_{\mathrm{X}}^{(r)} \dot{q}_r\right)\right]
\end{gathered} \tag{3.28}
$$
式中的$X$指代平台$P$、机舱$N$、轮毂$H$及电机$G$$x$指代过质心的转动轴的单位矢量,$I_{\mathrm{Xx}}$刚体$X$在$x$方向的转动惯性根据Kane方法 广义惯性力$F_r^*$为
$$
F_r^*=\sum_x F_x \cdot v_x^{(r)}+\sum_{\mathrm{X}} L_{\mathrm{X}} \cdot \boldsymbol{\omega}_{\mathrm{X}}^{(r)}(r=1,2, \cdots, 22) \tag{3.29}
$$
多体系统的主动力包括系统外部对系统内物体的作用力,也包括系统内物体的相互作用力。对于海上浮式风电机,广义主动力包括:广义重力 $ F^g_r $ 为
$$
F^g_r = m_d g Z \cdot v_d^{(r)} + \sum_i m_{ti} g Z \cdot v_{ti}^{(r)} + m_f g Z \cdot v_f^{(r)} + m_e g Z \cdot v_e^{(r)} + \sum_m \sum_i m_{bmi} g Z \cdot v_{bmi}^{(r)} \tag{3.30}
$$
叶素节点处的气动力 $ F_{bmi}^A $ 和气动力矩 $ L_{Bmi}^A $,广义气动力为
$$
F_r^{Aero} = \sum_m \sum_i F_{bmi}^A \cdot v_{bmi}^{(r)} + \sum_m \sum_i L_{Bmi}^A \cdot \omega_{Bmi}^{(r)} \tag{3.31}
$$
电机传动反力矩 $ L^G $,由变速变桨控制策略确定,对应广义力为
$$
\begin{cases}
F_{13}^G = -\lambda L^G \\
F_r^G = 0 \, (r \neq 13)
\end{cases} \tag{3.32}
$$
偏航刚度和阻尼产生的广义偏航力矩 $ L^{Yaw} $ 为
$$
\begin{cases}
L_{11}^{Yaw} = -K_{Yaw} q_{11} - C_{Yaw} \dot{q}_{11} \\
L_r^{Yaw} = 0 \, (r \neq 11)
\end{cases} \tag{3.33}
$$
柔性传动轴刚度和阻尼产生的广义力矩 $ L^{Dri} $ 为
$$
\begin{cases}
L_{13}^{Dri} = -K_{Dri} q_{13} - C_{Dri} \dot{q}_{13} \\
L_r^{Dri} = 0 \, (r \neq 13)
\end{cases} \tag{3.34}
$$
柔性塔架广义刚度和阻尼产生的广义恢复力 $ { }^{\mathrm{Ela}} \boldsymbol{F}_r^{\mathrm{T}} $ ($ r = 7, 8, 9, 10 $),和柔性叶片广义刚度和阻尼产生的广义恢复力 $ { }^{\mathrm{Ela}} \boldsymbol{B}_r^{\mathrm{T}} $ ($ r = 14, 15, \dots, 22 $) 都采用式 (1.20) 计算,$ r $ 为其它自由度时,广义恢复力为 0。
浮式平台参考点处受到的载荷分为两部分,即与平台运动加速度有关的力和无关的力,表示为
$$
F_i^{Platform} = F_i^{Platform1} + F_i^{Platform2} \\
\left\{\begin{array}{l}
F_i^{\text {Platform1 }}=-A_{i j} \ddot{q}_j \\
F_i^{\text {Platform2 }=}=F_i^{\text {Hydrostatic }}+F_i^{\text {Wave }}-\int_0^t K_{i j}(t-\tau) \dot{q}_j(\tau) d \tau+F_i^{\mathrm{Lines}}+{ }^{(2)} F_i+F_i^{\text {Monison }}
\end{array}\right. \tag{3.35}
$$
浮式平台参考点处受到的与平台运动加速度无关的力及力矩 $F_p^{Hydro}$、$L_p^{Hydro}$,与平台运动加速度有关的附加惯性力及力矩 ${ }^{\mathrm{A}} F_{\mathrm{p}}^{\text {Hydro }}$、${ }^{\mathrm{A}} L_p^{Hydro}$ 可表示为
$$
\left\{\begin{array}{l}
\boldsymbol{F}_p^{\text {Hydro }}=F_1^{\text {Plafform2 }} \boldsymbol{X}+F_2^{\text {Platform2 } } \boldsymbol{Y}+F_3^{\text {Platform2 }} \boldsymbol{Z} \\
\boldsymbol{L}_{\mathrm{p}}^{\text {Hydro }}=F_4^{\text {Platform2 } } \boldsymbol{X}+F_5^{\text {Platform2 }} \boldsymbol{Y}+F_6^{\text {platform2 }} \boldsymbol{Z}
\end{array}\right. \tag{3.36}
$$
$$
\left\{\begin{array}{l}
{ }^{\mathrm{A}} \boldsymbol{F}_{\mathrm{p}}^{\text {Hydro }}=\sum_{i=1}^6\left(-A_{1 i} X-A_{2 i} \boldsymbol{Y}-A_{3 i} Z\right) \ddot{q}_i \\
{ }^{\mathrm{A}} \boldsymbol{L}_{\mathrm{p}}^{\text {Hydro }}=\sum_{i=1}^6\left(-A_{3 i} X-A_{4 i} \boldsymbol{Y}-A_{5 i} Z\right) \ddot{q}_i
\end{array}\right. \tag{3.37}
$$
广义水动力 $F_r^{Hydro}$ 为
$$
F_t^{Hydro} = F_p^{Hydro}\cdot v_p^{(r)} + L_p^{Hydro}\cdot \omega_p^{(r)} + A_p^{Hydro}\cdot v_p^{(r)} + A_p^{Hydro}\cdot \omega_p^{(r)} \tag{3.38}
$$
总广义主动力 $F_r$ 为
$$
F_r=F_r^{\mathrm{g}}+F_r^{\text {Aero }}+F_r^{\mathrm{G}}+F_r^{\mathrm{Yaw}}+F_r^{\mathrm{Dri}}+{ }^{\text {Ela }} F_r^{\mathrm{T}}+{ }^{\text {Ela }} F_r^{\mathrm{B}}+F_r^{\text {Hydro }} \tag{3.39}
$$
对于有 22 个自由度的海上浮式风电机系统,凯恩动力学方程表示为
$$
F_r + F_r^* = 0 \, (r=1,2,...,22) \tag{3.40}
$$
通过整理动力学方程表示为
$$
M_{22\times 22}({q}_r, t)\ddot{q}_{22\times 1} = f_{22\times 1}(\dot{q}_r,q_r,t) \tag{3.41}
$$
式中 $M_{22\times 22}$为对称质量矩阵: $\ddot{q}_{22\times 1}$ 为 22 个自由度加速度构成的列向量, 有 $\ddot{q}_{22\times 1} = (\ddot{q}_1 \ddot{q}_2 ... \ddot{q}_{22})^T$$f_{22\times 1}(\dot{q}_r,q_r,t)$ 为与自由度加速度$\ddot{q}_r$无关的力向量。