当船舶和海洋结构物产生大幅参数共振时，会使结构遭受破坏和失效而导致人员和财产损失. 1998年10月，一艘超巴拿马型C11级集装箱船在迎浪海况中航行时发生了幅值达35^o~40^o的参数横摇运动，导致甲板上的集装箱大量损坏^[1]. 参数横摇可发生在随浪^[2]、迎浪^[3-4]和斜浪^[5]航行的船舶上. Yu 等^[6]通过检测纵摇和横摇的频率差来判断是否会发生参数横摇. Liu等^[7]利用基于雷诺平均Navier-Stokes方程的计算流体力学（computational fluid dynamics, CFD）程序模拟预测了内倾舰艇在迎浪规则波中的参数横摇. Ghamari等^[8]采用数值模拟和物理实验的方法，研究了渔船在无航速和有航速2种情况下的参数横摇. 其中数值模拟中Froude-Krylov（F-K）力和静水回复力的非线性影响通过积分到船体瞬时湿表面实现. 海洋浮式平台同样存在参数共振问题^[9-10]. Yang等^[11]通过数值模拟分析了锚泊缆绳、立管和阻尼对Spar平台参数稳定性的影响. Yang等^[12]分析了不规则波中Spar平台几何形状对参数纵摇的影响. Jang等^[13]采用与文献[8]类似的方法对极地Spar平台进行数值模拟，发现当波浪频率等于2倍纵摇/横摇固有频率时有参数共振产生，当波浪频率为垂荡固有频率的2倍时，观察到垂荡参数共振.

波浪能作为可再生能源，对其进行开发利用，对保障能源安全、实现可持续发展具有重要意义. 在已有的波浪能利用技术中，固定式振荡水柱（oscillating water column, OWC）^[14]或浮式OWC^[15-16]是重要的波浪能转换装置，某些小型浮式OWC型装置已成功为航标提供电力供给^[17]. Gomes等^[18]对松弛锚泊的Spar型OWC波浪能转换装置进行时域数值模拟和模型试验，发现当波浪周期为一半横摇/纵摇固有周期时会产生参数共振. 为了预测参数共振，Giorgi等^[19]利用非线性F-K模型模拟文献[18]中的装置，结果表明，参数共振产生时，波能俘获效率下降53%. 另一类受研究者关注的是浮式点吸收波浪能转换装置^[20-22]，Tarrant等^[23]采用F-K模型对双体式垂荡点吸收装置Wavebob进行数值模拟，假设动力输出（power take-off, PTO）机构阻尼力与2个物体间的相对垂荡速度成正比. 本研究认为，在发生大幅参数横摇/纵摇情况下，Wavebob装置PTO阻尼力应修改为与2个物体间的轴向相对速度成正比. 垂荡式点吸收装置只能俘获垂荡方向的波浪能，而多自由度波能转换装置能够比垂荡式点吸收装置俘获更多的能量^[24-25].

本研究 1）提出由浮体和承载结构组成的多自由度波浪能装置，承载结构与固定平台铰接，垂向PTO、纵摇PTO及横摇PTO分别用于俘获浮体沿铰接承载结构轴向滑动、浮体纵荡和纵摇及浮体横荡和横摇方向的能量. 2）利用势流理论建立该装置的频域和时域仿真模型. 3）探讨装置发生参数激励运动的原因，分析波高、周期、各PTO阻尼系数等对参数共振及波能俘获效率的影响.

如图1所示，多自由度波浪能转换装置由浮体组合结构（包括浮体、卡箍和撑杆）和承载结构（包括管状壳体、丝杆、推杆及内部构件）组成. 将垂向PTO机构安装在承载壳体内，并简称此管状承载结构为PTO管. 图中，oxyz为空间固定坐标系，原点o位于静水面，z轴垂直向上；G_bx_by_bz_b为固连于浮体的坐标系，原点G_b位于浮体重心处；G_px_py_pz_p为固连于PTO管的坐标系，原点G_p位于PTO管重心处；o′x′y′z′为固连于铰接点的坐标系，L为铰接点离水面高度. 在PTO管外侧设置滑轨，在通过撑杆与浮体固定相连卡箍的上下端设置滑轮，使卡箍能够沿PTO管上的滑轨滑动. 推杆的作用是将浮体的运动传递到滚珠丝杆机构的螺母上，使浮体组合结构能够沿PTO管往复运动，进而带动与发电机相连的丝杆旋转以捕获浮体沿PTO管滑动方向的能量. PTO管上端与固定平台铰接，并分别与纵摇和横摇方向的PTO机构连接，俘获浮体纵荡与纵摇方向、横荡与横摇方向的能量. 当装置的纵摇和横摇运动幅值较小时，垂向PTO用于俘获浮体垂荡运动方向的能量.

本研究所提装置总共有3个PTO机构，每个PTO机构均有独立的发电机与之相连. 无论装置的横摇运动和纵摇运动如何变化，与垂向PTO连接的发电机均能正常发电，但参数横摇、参数纵摇会引起垂向发电机输出功率的减小，因此有必要采取措施抑制参数激励运动. 遭遇随浪或迎浪海况时，与垂向PTO连接的发电机和安装在纵摇铰接轴上的发电机可以正常发电. 遭遇横浪海况时，与垂向PTO连接的发电机和安装在横摇铰接轴上的发电机可以正常发电. 遭遇斜浪海况时3个发电机均可以正常发电. 因此本研究所提装置能够同时俘获浮体多个运动自由度方向的能量.

由于PTO管铰接约束下的装置无首摇运动，在线性规则波作用下浮体相对于自身重心位置做频率为ω的五自由度简谐运动

(1) $ {\xi }_{j}(t)=\mathrm{Re}\left({\xi }_{j0}{\text{e}}^{-\text{i}\omega \text{ }t}\right)\text{，} j=1, 2, …, 5. $

式中：ξ_j(t)为t时刻浮体重心位置j自由度的运动位移，ξ_j0为复振幅. PTO管绕铰接点o′作横摇（j=4）和纵摇(j=5) 2种自由度运动，且浮体与PTO管具有相同的纵/横摇运动. 装置的运动幅值通过拉格朗日运动方程求解：

(2) $ \frac{{\text{d}}}{{{\text{d}}t}}\left( {\frac{{\partial {E_{\text{k}}}}}{{\partial {{\dot \xi }_j}}}} \right) - \frac{{\partial {E_{\text{k}}}}}{{\partial {\xi _j}}}+\frac{{\partial {E_{\text{P}}}}}{{\partial {\xi _j}}}+\frac{{\partial {E_{\text{D}}}}}{{\partial {{\dot \xi }_j}}} = {F_{\xi _{_j}}}. $

式中：E_k为装置动能，E_P为势能，E_D为负载阻尼吸收的能量； ${ {F_\xi }_ {_j}}$为与ξ_j有关的广义作用力，包括作用在浮体上的静水力，波浪扰动力和辐射流体动力载荷. 在小幅值运动假设下，浮体重心位置的纵荡和纵摇、横荡及横摇运动间存在如下关系：

(3) $ {\xi }_{1}=-(L-{z}_{G_\text{b}}){\xi }_{5}\text{，} {\xi _2} = (L - {z_{G{_\text{b}}}}){\xi _4}. $

式中： ${z_{G_{\text{b}}}}$为浮体重心在空间固定坐标系下的垂向坐标值，ξ₁、ξ₂分别为浮体重心位置的纵荡和横荡运动，ξ₄、ξ₅分别为装置的横摇和纵摇. 结合式（1）、（3）可知，在线性规则波作用下，轴对称型波浪能装置可用浮体垂荡(j=3)、装置相对于铰接点o′的横摇(j=4)和纵摇(j=5) 3个单自由度方程来描述其运动响应. 经推导，将式（2）中各自由度的运动复幅值表示为

(4) $ \left. {\begin{array}{*{20}{l}} {{\xi _{30}} = \dfrac{{{\zeta _{\text{a}}}{E_{30}}}}{{{K_{33}}+{K_{{\text{p3}}}} - {\omega ^2}({m_{\text{b}}}+{A_{33}}) - {\text{i}}\omega ({B_{33}}+{B_{{\text{p}}3}})}};} \\ {{\xi _{j0}} = \dfrac{{{\zeta _{\text{a}}}E_{j0}^{o'}}}{{K_{jj}^{o'}+{K_{{\text{p}}j}} - {\omega ^2}(I_{jj}^{o'}+A_{jj}^{o'}) - {\text{i}}\omega (B_{jj}^{o'}+{B_{{\text{p}}j}})}},j = 4,5.} \end{array}} \right\} $

(5) $ \left. {\begin{array}{*{20}{l}} {{K_{33}} = \rho g{A_{{\text{WP}}}};} \\ {K_{44}^{o'} = \rho g{V_0}{{\overline {{\rm{GM}}} }_{\text{T}}}+{m_{\text{p}}}g(L - {z_{G{\text{p}}}}),} \\ {K_{{\text{55}}}^{o'} = \rho g{V_0}{{\overline {{\rm{GM}}} }_{\text{L}}}+{m_{\text{p}}}g(L - {z_{G{\text{p}}}}).} \end{array}} \right\} $

(6) $ I_{jj}^{o'} = I_{jj}^{{G_{_\text{b}}}}+I_{jj}^{{G_{\text{p}}}}+{m_{\text{p}}}{(L - {z_{G{_\text{p}}}})^2}+{m_{\text{b}}}{(L - {z_{G_{\text{b}}}})^2},j = 4,5 . $

(7) $ \left. {\begin{array}{*{20}{c}} {A_{44}^{o'} = {A_{44}}+{A_{22}}{{(L - {z_{G{_\text{b}}}})}^2}+({A_{24}}+{A_{42}})(L - {z_{G{_\text{b}}}}),} \\ {A_{55}^{o'} = {A_{55}}+{A_{11}}{{(L - {z_{G{_\text{b}}}})}^2} - ({A_{15}}+{A_{51}})(L - {z_{G{_\text{b}}}}).} \end{array}} \right\} $

(8) $ \left. {\begin{array}{*{20}{c}} {B_{44}^{o'} = {B_{44}}+{B_{22}}{{(L - {z_{G{_\text{b}}}})}^2}+({B_{24}}+{B_{42}})(L - {z_{G{_\text{b}}}}),} \\ {B_{55}^{o'} = {B_{55}}+{B_{11}}{{(L - {z_{G{_\text{b}}}})}^2} - ({B_{15}}+{B_{51}})(L - {z_{G{_\text{b}}}}).} \end{array}} \right\} $

(9) $ \left. {\begin{array}{*{20}{c}} {E_{40}^{o'} = {E_{40}}+{E_{20}}(L - {z_{G{_\text{b}}}}),} \\ {E_{50}^{o'} = {E_{50}} - {E_{10}}(L - {z_{G{_\text{b}}}}).} \end{array}} \right\} $

式中：m_b为包括浮体、撑杆和卡箍在内的组合浮体质量，m_p为PTO管质量， ${z_{G_{\text{p}}}}$为PTO管重心在空间固定坐标系下的垂向坐标， $ I_{jj}^{{G_{\text{b}}}} $(j=4, 5)为组合浮体相对于浮体自身重心G_b的质量惯性矩， $I_{jj}^{{G_{\text{p}}}}$_{(j=4, 5)为PTO管相对于自身重心Gp的质量惯性矩，ρ为水密度，g为重力加速度，AWP为浮体水线面面积，V0为浮体静浮时的排水体积，Kpj、Bpj(j=3, 4, 5)分别为浮体沿PTO管滑动、装置横摇和纵摇运动方向的线性负载刚度系数和负载阻尼系数， $ {\overline {{\rm{GM}}} _{\text{T}}} $、 $ {\overline {{\rm{GM}}} _{\text{L}}} $分别为浮体的横稳性高和纵稳性高，Ej0为单位波幅产生的作用在浮体重心Gb第j自由度方向的波浪扰动力复幅值，ζa为入射波波幅，Aij、Bij分别为浮体相对于自身重心Gb的附加质量和辐射阻尼系数. 规则波作用下垂向PTO(j=3)、横摇PTO（j=4）和纵摇PTO（j=5）俘获的平均功率为}

(10) $ P_{{\text{m}}j}^{} = 0.5{\omega ^2}B_{{\text{p}}j}^{}{\left| {{\xi _{j0}}} \right|^2},\; j=3, 4, 5. $

将式（4）代入式（10），得到各自由度方向的平均输出功率. 为了最大化各PTO的平均输出功率，式（10）对B_pj求导，并令导数等于零，得到各PTO对应的最佳负载阻尼系数B_optj. 令K_pj=0，有

(11) $ \left. {\begin{array}{*{20}{l}} {B_{{\text{opt}}3}^{} = \left( {{{\left( {\dfrac{{{K_{33}} - {\omega ^2}({m_{\text{b}}}+{A_{33}})}}{\omega }} \right)}^2}+{{\left( {{B_{33}}} \right)}^2}} \right)^{0.5};} \\ B_{{\text{opt}}j}^{} = \left( {{{\left( {\dfrac{{K_{jj}^{o'} - {\omega ^2}(I_{jj}^{o'}+A_{jj}^{o'})}}{\omega }} \right)}^2}+{{\left( {B_{jj}^{o'}} \right)}^2}} \right) ^{0.5},\;j = 4,5. \end{array}} \right\} $

对本垂直轴对称波浪能装置而言，利用频域方法求解时，通过式（11）可以很方便地求出各自由度对应的线性最佳负载阻尼系数B_optj（j=3, 4, 5）.

若波浪能装置发生参数激励运动，则大幅值的纵摇或横摇运动将影响装置的输出功率和波浪能俘获效率. 由于线性频域运动方程无法分析参数激励运动，对于可能产生参数共振的非线性系统须采用时域方法求解，利用Cummins方程描述辐射力. 设组合浮体为结构1，PTO管为结构2，则将浮体的时域运动方程表示为

(12) $ \left. {\begin{array}{*{20}{l}} {\left( {{{\boldsymbol{m}}_1}+{{\boldsymbol{A}}_\infty }} \right){{{\boldsymbol{\ddot \zeta }}}_1}(t)+{{\boldsymbol{c}}_1}{{{\boldsymbol{\dot \zeta }}}_1}(t)+\displaystyle\int_0^t {{\boldsymbol{R}}(t - \tau ){{{\boldsymbol{\dot \zeta }}}_1}(\tau ){\text{d}}\tau } } =\\ \quad\quad { {{\boldsymbol{F}}_1}(t)+{\boldsymbol{T}}(t)+{{\boldsymbol{F}}_{{\text{p}}1}}(t);} \\ {{\boldsymbol{R}}(t) = {2 {\text{π}} }^{-1}\displaystyle\int_0^\infty {{\boldsymbol{B}}(\omega )} \cos \omega t{\text{d}}\omega ;} \\ {{{\boldsymbol{A}}_\infty } = {\boldsymbol{A}}(\omega )+{\omega }^{-1}\displaystyle\int_0^\infty {{\boldsymbol{R}}(t)\sin \;} \omega t{\text{d}}\omega .} \end{array}} \right\} $

式中： ${{\boldsymbol{\zeta}} }_{\text{1}}={\left[\begin{array}{ccccc}{\xi }_{1}, \;\; {\xi }_{2}, \;\; {\xi }_{3}, \;\; {\xi }_{4}, \;\; {\xi }_{5}\end{array}\right]}^{\text{T}}$为浮体重心位置的运动位移矩阵， ${{\boldsymbol{\dot \zeta }}_1}(t)$为浮体速度矩阵， ${{\boldsymbol{\ddot \zeta }}_1}(t)$为浮体加速度矩阵， ${{\boldsymbol{m}}_1}$为浮体的质量矩阵，A_∞为频率趋于无穷时浮体的附加质量矩阵， $\tau $为时间差，A(ω)、B(ω)分别为浮体广义附加质量矩阵和辐射阻尼矩阵，c₁为包括辐射阻尼在内的浮体线性阻尼矩阵，R为浮体的速度脉冲函数矩阵，F₁为作用在浮体上的重力、静水力、流体黏性力、波浪力/力矩等外力，F_p1为作用在浮体上的PTO负载力/力矩，T为卡箍和PTO管间的接触力/力矩. PTO管的时域运动方程为

(13) $ {{\boldsymbol{m}}_{\boldsymbol{2}}}{{\boldsymbol{\ddot \zeta }}_2}(t) = {{\boldsymbol{F}}_2}(t) - {\boldsymbol{T}}(t)+{{\boldsymbol{F}}_{{\text{p2}}}}(t). $

式中： ${\zeta }_{2}={\left[\begin{array}{ccccc}{\xi }_{2,1}, \;\; {\xi }_{2,2}, \;\; {\xi }_{2,3}, \;\; {\xi }_{4}, \;\; {\xi }_{5}\end{array}\right]}^{\text{T}}$为PTO管重心G_p处的位移矩阵（PTO管和浮体具有相同的横摇和纵摇运动）， $ {{\boldsymbol{\ddot \zeta }}_{\boldsymbol{2}}} $为PTO管重心位置的加速度矩阵， ${{\boldsymbol{m}}_2}$为PTO管的质量矩阵，F₂为作用在PTO管上的重力/力矩及作用在PTO管铰接点处的约束力/力矩之和，F_p2为作用在PTO管上的PTO装置负载力. 设浮体静浮时PTO管重心和浮体重心间的距离为 ${l_0} = {z_{G_{\text{p}}}} - {z_{G_{\text{b}}}}$. 在波浪作用下，PTO管重心和浮体重心间的瞬时距离为

(14) $ l = \sqrt {{{({\xi _{{\text{2}},1}} - \xi _1^{})}^2}+{{({\xi _{{\text{2}},2}} - {\xi _2})}^2}+{{\left( {{l_0}+{\xi _{{\text{2}},3}} - {\xi _3}} \right)}^2}} . $

相对于平衡位置，定义浮体沿PTO管向上滑动时的位移Δl为正，有 $ \Delta l={-}(l{-}{l}_{0}) $，在规则波作用下，

(15) $ \Delta l = {\text{Re}}\left( {\Delta {l_{\text{a}}}{{\text{e}}^{{{ - {\rm{i}}}}\omega t}}} \right). $

式中： $ \Delta {l_{\text{a}}} $为浮体沿PTO管滑动的位移复幅值. 浮体沿PTO管滑动的瞬时速度

(16) $ \Delta \dot l = \frac{{{\text{d}}(\Delta l)}}{{{\text{d}}t}} = \sum\limits_{j = 1}^3 {} \frac{{\partial (\Delta l)}}{{\partial {\xi _j}}}{\dot \xi _j}+\sum\limits_{j = 1}^3 {} \frac{{\partial (\Delta l)}}{{\partial {\xi _{{\text{2}},j}}}}{\dot \xi _{{\text{2}},j}}. $

垂向PTO作用在结构1的阻尼力F_p1在固定坐标系x、y和z方向的分量

(17) $ {F_{{\text{p1,}}j}} = - {B_{{\text{p3}}}}\Delta \dot l\frac{{\partial (\Delta \dot l)}}{{\partial {{\dot \xi }_j}}},\;j = 1,{\text{ }}2,{\text{ }}3. $

垂向PTO、横摇PTO及纵摇PTO作用在结构2上的阻尼力/力矩F_p2中各分量

(18) $ {F_{{\text{p2,}}j}} = \left\{{\begin{array}{*{20}{l}} { - {B_{{\text{p3}}}}\Delta \dot l\dfrac{{\partial (\Delta \dot l)}}{{\partial {{\dot \xi }_{2,j}}}},\;j = 1,{\text{ }}2,{\text{ }}3;} \\ { - {B_{{\text{p}}j}}\dot \xi _j^{}, \qquad \;\;\; j = 4,{\text{ }}5.} \end{array}} \right. $

装置由2个结构组成，2个结构体间的接触力传递通过设置刚度系数很大的护舷实现^[26]，使得浮体和PTO管具有同步的横摇和纵摇运动. 时域运动方程采用AQWA-NAUT模块求解，其中与发电机相连的PTO负载力和流体黏性力的计算在User Force程序中编写，使用Fortran编译动态链接库文件（.dll）实现外部载荷的定义和调用.

设垂直轴对称浮体水线面处直径为D，规则波中直径D范围内的入射波平均功率为

(19) $ P_{{\text{in}}}^{} = \frac{{\rho g\omega {H^2}D}}{{16k}}\left( {1+\frac{{2kd}}{{\sinh \,(2kd)}}} \right) . $

式中：k为波数；d为水深；H为波高，H=2ζ_a. 采用频域方法求解时，规则波作用下各PTO机构俘获的平均功率P_mj由式（10）计算得到. 采用时域方法求解时，垂向PTO、横摇PTO和纵摇PTO各自俘获的平均功率分别表示为

(20) $ \left. {\begin{array}{*{20}{l}} {P_{{\text{m}}3}^{} = {{(2{T_{\text{w}}}})}^{-1}\displaystyle\int_0^{2{T_{\text{w}}}} {B_{{\text{p}}3}^{}{{\left( {\Delta \dot l} \right)}^2}{\text{d}}t;} } \\ {P_{{\text{m}}j}^{} = {({2{T_{\text{w}}}})}^{-1}\displaystyle\int_0^{2{T_{\text{w}}}} {B_{{\text{p}}j}^{}\dot \xi _j^2(t){\text{d}}t} ,\;j = 4,{\text{ }}5.} \end{array}} \right\} $

式中：T_w为波浪周期. 式（20）表示在2个波浪周期范围内求各PTO的平均输出功率，目的是考虑参数共振的影响，因为参数激励运动主要发生在波浪频率为2倍纵摇或横摇固有频率及其附近频率范围. 将各PTO俘获的输出功率相加得到装置总的输出功率P_T。各PTO对应的波能俘获宽度比

(21) $ {\eta }_{j}^{}={{P}_{\text{m}j}}/{{P}_{\text{in}}}\text{，}j=3\text{，}4\text{，}5 . $

式（20）得到的结果为发电机前端的功率，即式（21）计算的俘获宽度比是没有考虑发电机本身效率时的结果，也称为一次转换效率，且在求解运动方程式（12）、（13）时忽略丝杆摩擦损失的影响. 装置总波能俘获宽度比定义为各PTO的波能俘获宽度比之和。

如图2所示，根据缩尺比为1∶4的多向联动复合PTO机构，建立装置水动力性能数值仿真模型. 取计算水深为1.6 m，模型吃水为0.43 m，静浮时浮体重心在空间固定坐标系下的垂向坐标 ${z_{G_{\text{b}}}}$=−0.041 m，质量m_b=75.5 kg，浮体相对于自身重心的质量惯性矩 $I_{{\text{44}}}^{{G_{\text{b}}}} 、 I_{{\text{55}}}^{{G_{\text{b}}}}$均为6.8 kg·m². PTO管加上推杆的最大高度为1.62 m，最小高度为1.2 m. PTO管上端铰接点离静水面距离L=1.518 5 m，PTO管质量m_p=30 kg，PTO管重心垂向坐标 ${z_{G_{\text{p}}}} $=1.005 3 m，PTO管相对于自身重心的质量惯性矩 $ I_{{\text{44}}}^{{G_{\text{p}}}} 、 I_{{\text{55}}}^{{G_{\text{p}}}} $均为4.8 kg·m². 装置纵摇和横摇的固有周期均为8.7 s.

基于线性的频域方法无法模拟参数共振问题，而时域方法可以处理非线性力的作用. 为了寻找参数激励运动发生的原因，时域方法中的F-K力和静水力均采用线性方法计算，因此所有流体作用力项均是线性的，且假设垂向PTO的阻尼力与浮体沿PTO管的滑行速度成正比. 在线性波作用下的装置也有可能产生大幅值参数共振，由式（17）可知，此时垂向PTO力分解到水平和垂直方向时，须保留非线性项的贡献. 为了探讨垂向PTO的负载阻尼系数对装置运动响应的影响，取随浪工况，波高为0.2 m进行数值模拟，只有垂向PTO连接负载，纵摇和横摇PTO皆无负载. 取B_p3=0、500、800、1 140.4 N·s/m进行数值仿真，其中B_p3=1 140.4 N·s/m为周期2.5 s时垂荡方向的最佳负载阻尼系数.

如图3所示为采用时域计算方法得到的装置运动响应与频域结果的比较。图中， ${\zeta'{_3}} $为浮体的垂荡运动频率响应函数， ${\zeta'{_3}}=\left|\xi_{30}\right| /{ \zeta_{\mathrm{a}}}$； ${\zeta '{_{\rm{a}}}}$为浮体沿PTO管的位移， ${\zeta '{_{\rm{a}}}} = \left| {\Delta {l_{\rm{a}}}} \right|/{\zeta _{\rm{a}}} $； ${\zeta '{_5}} $为装置的纵摇运动频率响应函数， ${\zeta '{_5}} = \left| {{\xi _{50}}} \right|/{\zeta _{\rm{a}}} $。垂向PTO输出的平均功率及其波能俘获宽度比如图4所示. 频域结果对应沿PTO管滑动方向的负载阻尼力，分解到水平和垂向时只保留线性项，即只有垂荡方向具有PTO力作用. 本算例时域计算结果显示，当波浪频率在2倍纵摇固有频率附近，会发生严重的参数纵摇现象. 与参数纵摇相比，参数横摇运动相对较小且随时间不规则变化，无法有效统计，因此本研究未给出横摇运动响应的数据. 在无参数共振发生的情况下，时域结果与频域结果吻合良好. 当发生参数共振时，大幅值纵摇运动导致浮体的垂荡运动幅值增加，但沿PTO管滑动的位移幅值减小，导致垂向PTO的输出功率和波能俘获宽度比明显下降. 在波浪周期为4.5 s，B_p3=500、800、1 140.4 N·s/m时，由于发生参数共振，垂向PTO的波能俘获宽度比与不考虑参数共振的频域结果相比，分别下降11.2%、41.2%和53.6%. 当B_p3=500 N·s/m时，参数共振发生的周期范围较窄，随着B_p3的增大，发生参数共振的周期范围扩大，且参数纵摇的幅值也增大. 本算例结果表明，在垂向PTO阻尼系数B_p3=0时，无参数共振产生，因此本研究建立的装置产生参数共振是垂向PTO阻尼力引起的.

随浪工况，垂向PTO连接负载，纵摇PTO和横摇PTO皆无负载，取模型周期2.5 s时垂荡方向的最佳负载阻尼系数B_p3=B_opt3=1 140.4 N·s/m，并取H=0.1、0.15、0.2 m共 3种波高进行数值模拟. 如图5所示为3种波高对应的浮体沿PTO管滑动位移和装置纵摇运动频率响应函数随波浪周期的变化曲线. 基于线性方程的频域方法所得频率响应函数在B_p3、B_p5一定的情况下不随波高变化，因此只有1条曲线. 在B_p3一定的情况下，由B_p3引起的参数纵摇的频率范围随波高增大而变宽，而纵摇频率响应函数的最大值随波高增大而降低. 如图6所示为3种波高对应的垂向PTO的平均输出功率和波能俘获宽度比随周期变化的曲线. 由图5、6可知，大幅值参数纵摇运动导致浮体沿PTO管滑动的位移变小，造成垂向PTO的输出功率降低，且随着波高增大，发生参数共振的区域向高频（短周期）方向扩展. 当波高为0.1 m时，与频域结果比较，在参数共振区浮体沿PTO管的滑动位移幅值、垂向PTO的输出功率及波能俘获宽度比均下降.

随浪工况，垂向PTO和纵摇PTO皆连接负载用于多自由度波能发电，横摇PTO无负载. 取T_w=2.5 s时垂荡方向的最佳负载阻尼系数B_p3=B_opt3=1 140.4 N·s/m，波高H=0.1 m和纵摇PTO阻尼系数B_p5=0、42.3、344.4、788.0 N·m·s/rad进行数值模拟. B_p5=42.3、344.4、788.0 N·m·s/rad分别对应T_w=8.0、4.5、2.5 s时纵摇方向的最佳负载阻尼系数.

如图7所示为3种B_p5对应的装置运动频率响应函数随波浪周期的变化曲线。图中， ${\zeta '_4} $为装置的横摇运动响应函数， ${\zeta '_4} = \left| {{\xi _{40}}} \right|/{\zeta _{\rm{a}}} $。图8、9分别为垂向PTO和纵摇PTO的输出功率及其对应的波能俘获宽度比. 由式(4)可知，基于线性方程的频域方法所得垂荡运动频率响应函数只与B_p3有关，与B_p5无关，因此在B_p3一定的情况下，浮体垂荡运动响应不随B_p5变化. 由图7可知，在没有发生参数共振的情况下，时域结果与频域结果吻合良好. 时域计算结果表明，在B_p3=1 140.4 N·s/m，B_p5=0 N·m·s/rad时，主要发生参数纵摇运动，参数横摇幅值较小（不规则小幅运动，图中忽略）；取较小的纵摇PTO阻尼系数（B_p5=42.3 N·m·s/rad）时，模型周期T_w=3.8 ~4.7 s，参数纵摇和参数横摇同时发生；取较大的B_p5（B_p5=344.4、788.0 N·m·s/rad）时，参数纵摇受到抑制，但参数横摇运动的幅值较大. 与参数纵摇情况相比，参数纵摇受到抑制而转向参数横摇时，浮体沿PTO管滑动位移的下降幅值相对较大，导致如图8所示垂向PTO的输出功率和波能俘获宽度比的下降幅度增大. 图9结果显示，取B_p3=1 140.4 N·s/m，B_p5=42.3 N·m·s/rad时，与频域结果比较，在参数纵摇和参数横摇同时产生的频率范围，纵摇PTO的输出功率及波能俘获宽度比有所增加，但由于负载阻尼系数太小，在图示频率范围的输出功率过低.

如图10所示为装置总输出功率P_T和总波能俘获宽度比 $\eta^{}_{\rm{T}} $随周期及B_p5变化的曲线. 结果表明，与垂向PTO单独获能的情况相比，纵摇PTO和垂向PTO同时获能发电时的输出功率及波能俘获宽度比均有所提升.

随浪工况，取T_w=2.5 s时垂荡方向的最佳负载阻尼系数B_p3=B_opt3=1 140.4 N·s/m，纵/横摇PTO阻尼系数B_p4=B_p5=B_opt5=788.0 N·m·s/rad，波高H=0.1 、0.2 m进行数值模拟. 随浪工况时的横摇PTO不发电，对横摇施加PTO阻尼力矩的目的是抑制参数横摇. 如图11所示为参数共振得到抑制后总输出功率和总波能俘获宽度比随周期变化的曲线. 可以看出，参数共振得到抑制后，时域与频域曲线吻合良好.

（1）采用线性势流理论建立由浮体和PTO管组成的轴对称型多自由度波浪能装置的频域计算模型. 模型可以简化成3个单自由度方程，具有简单快速的计算优势. 建立时域计算模型，假设垂向PTO阻尼力与浮体沿PTO管滑动的速度成正比，纵/横摇方向PTO阻尼力矩与纵/横摇的角速度成正比，在将垂向PTO力分解到水平和垂直方向时，保留非线性项作用以评估参数激励运动. 研究结果表明，在没有产生参数激励运动的频率范围，时域计算结果与频域结果吻合良好.

（2）随浪工况，纵摇PTO和横摇PTO皆不连负载，只有在垂向PTO施加负载阻尼用于俘获波浪能的过程中，在波浪频率为2倍横摇/纵摇固有频率附近范围，有严重的参数纵摇运动，此时的参数横摇幅值较小，且随时间的变化不规则. 当有参数共振时，浮体沿PTO管滑动的位移减小，导致垂向PTO的输出功率及波能俘获宽度比下降. 发生参数纵摇运动的频率范围随波高和垂向PTO阻尼系数B_p3的增大而变宽；当有参数纵摇时，纵摇运动频率响应函数的最大值随波高增大而降低，随B_p3的增大而增大. 随着波高和B_p3的增大，参数共振的区域向高频（短周期）方向扩展. 随浪工况，横摇PTO不连负载，在垂向PTO和纵摇PTO连接负载用于多自由度获能的过程中，纵摇PTO阻尼系数B_p5较小时，参数纵摇和参数横摇运动同时发生；B_p5较大时，参数纵摇受到抑制，但发生较大幅值的参数横摇运动. 参数纵摇受到抑制情况下的参数横摇运动幅值大于参数纵摇和横摇同时产生时的参数横摇幅值.

（3）随浪工况，对纵摇PTO和横摇PTO均施加一定大小的PTO阻尼力矩系数时，能够有效抑制参数共振. 随浪工况，与垂向PTO单独获能比较，纵摇PTO和垂向PTO同时获能时的总输出功率及总波能俘获宽度比均有所提升.

（4）本研究只分析了规则波中的参数激励运动及对波能俘获能力的影响，后续计划针对随机波中的问题开展进一步研究.