波浪能是海洋能中品位最高、分布最广的能源形式之一^[1]，作为清洁的可再生能源，对缓解能源危机、减少对传统化石燃料的依赖、改善环境污染等问题具有重要的现实意义. 离岸式振荡水柱（off-shore oscillating water column，off-shore OWC）波能转换装置由于结构简单、易于安装维护、转换效率高等优点而被广泛应用，其基本工作原理是利用海浪上下起伏运动造成气室内外空气压强差，推动气室顶部的空气透平旋转，从而把波浪能转换成机械能.

国内外学者已对OWC波能转换装置进行了不少研究. Evans^[2]和Falnes^[3]在线性波理论框架下，提出OWC波能转换理论. Evans^[4]对早期理论进行改进，将OWC气室内的压强简化为均匀分布的空气压强，并允许气室内水体的自由液面在空间上发生变形，不需要一定保持为完整的平面，得到OWC装置最大和最优效率的解析表达式. Evans等^[5]应用特征函数展开法，研究有限水深情况下OWC简化模型的水动力学特性，通过近似表征流体流动在结构大变形处的奇异行为，有效提升数值计算的计算效率和数值精度. 在物理模型试验方面，Ambli等^[6]提出多共振振荡水柱装置，利用港口效应聚集波浪能，提高能量转换效率. 梁贤光等^[7]研究在不同入射波参数和管口情况下，OWC气室内的水面运动情况. 史宏达等^[8-9]对沉箱式防波堤兼作岸式波力发电装置在非规则波下的能量转换效率进行试验研究. 在数值方面，You^[10]运用数值手段考察地形对近岸OWC装置工作效率的影响. 史宏达等^[11]利用雷诺时均方程和VOF模型，计算波浪的传播和水柱在气室内的振荡情况. 宁德志等^[12-13]基于域内源造波技术的时域高阶边界元方法，研究岸式OWC波能转换装置，得到共振频率与前墙吃水深度、墙体厚度及气室宽度之间的关系. Kamath 等^[14]依据达西方程，模拟不同输出阻尼系数下的OWC水动力学响应特性. Luo等^[15]考虑入射波波高变化对波能捕获效率的影响. Wang等^[16]分析波浪的非线性和黏性效应对OWC装置工作效率的影响.

从理论求解和数值模拟2个角度出发，首先基于势流理论，采用匹配特征函数展开法求解微幅波与离岸式振荡水柱波能转换装置相互作用的边界值问题，其次利用商业软件FLUENT建立二维完全非线性波-OWC数值模型，深入研究一系列参数对装置能量转换效率的影响，参数包括OWC装置的吃水深度、气室宽度以及墙体厚度等. 此外，基于经过验证的OWC数值模型，进一步研究入射波波高变化对OWC装置转换效率的影响情况.

如图1所示，OWC波能转换装置安装于静水深为h的波浪水槽中，气室尺寸为2w₁c. 装置吃水深度为d，墙体厚度为w₂−w₁. 将直角坐标系的原点置于气室内静水面的中点，x轴与静水面重合，z轴垂直向上，波浪沿x轴正向传播. 根据边界条件的差异，将整个流场划分为5个区域（Ω₁、Ω₂、Ω₃、Ω₄、Ω₅）. 记入射波频率为ω，速度势在复数域中可表示为 ${\rm{\varPhi }}\left( {x,\,z,\,t} \right) = {\rm{Re\,[}}\phi {\left( {x,\,z} \right)^{ - {\rm i}\omega t}}]$. 其中，t为时间变量；i为虚部；Re为取变量实部的运算；ϕ(x, z)为分离出时间变量的速度势，在整个流场中满足拉普拉斯方程：

(1) $\frac{{{\partial ^2}\phi (x,z)}}{{\partial {x^2}}} + \frac{{{\partial ^2}\phi (x,z)}}{{\partial {z^2}}} = 0.$

相应的自由表面边界条件为

(2) $\frac{{\partial \phi }}{{\partial z}} - \frac{{{\omega ^2}}}{g}\phi = \left\{ \begin{gathered} {{{\rm i}\omega }}{p_{\rm{c}}}/({{\rho g}}),\quad \left| x \right| < {w_1},\,z = 0; \\ 0,\quad\quad\,\,\;\,\quad\quad \left| x \right| > {w_2},\,z = 0. \\ \end{gathered} \right.$

式中：p_c为拟封闭气室内的压强，g为重力加速度，ρ为海水密度.

由于结构物处于固定状态，物面及海底边界条件均为

(3) ${{\partial \phi }}/{{\partial n}} = 0.$

式中：n为物面法线方向，指向流体区域外. 此外，ϕ(x, z)在远场还需满足Sommerfeld辐射边界条件.

波浪与OWC装置相互作用的过程可分为波浪衍射及辐射问题，即

(4) $\phi = {\phi ^{\rm{S}}} + {\phi ^{\rm{R}}}.$

式中：ϕ^S为衍射势，在气室内外气压一致时，由波浪与OWC装置前后2个墙体相互作用产生；ϕ^R为在没有入射波的情况下，气室内压力振荡作用于液面而产生的辐射速度势.

在区域Ω₁、Ω₃、Ω₅中，上部为自由表面边界（在式（2）中令p_c=0），下部为不可穿透的海底边界. 因此，可将速度势展开为特征函数基 ${Z_m}\left( {{k_m}z} \right)(m = 0,1,2, \cdots )$ 线性叠加的形式. ${Z_m}\left( {{k_m}z} \right) = N_m^{ - 1/2}\cos\;\left[{k_m}\left( {z + h} \right)\right]\left(m = 0,1,2, \cdots \right)$ ，特征函数在区间 $z \in ( - h,\,0)$ 上相互正交，其中：

(5) $N_m=\left\{ \begin{aligned} &\left[2kh + \sinh\; (2kh)\right]/(4k),\quad\;\, m={\rm 0}; \\ &\left[2{k_m}h + \sin\; (2{k_m}h)\right]/(4{k_m}),\,m=1,2,3 \cdots . \end{aligned} \right.$

式中：k为行进波的波数. k、k_m可分别由以下方程求出：

(6) $\left. \begin{gathered} {\omega ^2}{\rm{ = }}gk\tanh\; (kh), \\ {\omega ^2} + g{k_m}\tan\; ({k_m}h){\rm{ = 0}}. \\ \end{gathered} \right\}$

取k₀=−ik，将相应区域的速度势分别展开为

(7) $\begin{split} \phi _1^{\rm{S}}(x,z) =& - \frac{{{\rm i}gA}}{\omega }\left[\displaystyle\frac{{N_0^{1/2}}}{{\cosh\;(kh)}}{{\rm {exp}}\;\left[{{\rm i}k(x + {w_2})}\right]}{Z_0}({k_0}z) +\right. \\ &\left. \sum\limits_{m = 0}^\infty {A_m^{\rm{S}}{{\rm {exp}}\;\left[{{k_m}(x + {w_2})}\right]}{Z_m}({k_m}z)} \right], \end{split} $

(8) $\phi _3^{\rm{S}}(x,z) = - \frac{{{\rm i}gA}}{\omega }\sum\limits_{l = 0}^\infty {\left[D_l^{\rm{S}}{{\rm {exp}}\;({{k_l}x}}) + E_l^{\rm{S}}{{\rm {exp}}\;({ - {k_l}x}})\right]{Z_l}({k_l}z)}, $

(9) $\phi _5^{\rm{S}}(x,z) = - \frac{{{\rm i}gA}}{\omega }\sum\limits_{m = 0}^\infty {H_m^{\rm{S}}{{\rm {exp}}\;\left[{ - {k_m}(x - {w_2})}\right]}{Z_m}({k_m}z)} .$

式中：A为入射波振幅， $A_m^{\rm S}$、 $D_l^{\rm S}$、 $E_l^{\rm S}$、 $H_m^{\rm S}$均为未知系数.

在区域Ω₂、Ω₄中，上下为不可穿透的固壁边界，相应的速度势展开为

(10) $ \begin{split} \phi _2^{\rm{S}}(x,z) =& - \displaystyle\frac{{{\rm i}gA}}{\omega }\Bigg\{ \left[B_0^{\rm{S}} + C_0^{\rm{S}}(x + {w_1})\right]{Y_0}({\lambda _0}z){\rm{ + }}\Bigg.\\ & \displaystyle\sum\limits_{n = 1}^\infty {\left[B_n^{\rm{S}}{{\rm {exp}}\;\left[{{\lambda _n}(x + {w_1})}\right]} \right.}+\\& \Bigg.\left.C_n^{\rm{S}}{{\rm {exp}}\;\left[{ - {\lambda _n}(x + {w_1})}\right]}\right]{Y_n}({\lambda _n}z) \Bigg\} , \end{split} $

(11) $ \begin{split} \phi _4^{\text{S}}(x,z) =& - \frac{{{\rm i}gA}}{\omega }\Bigg\{ {\left[ {F_0^{\text{S}} + G_0^{\text{S}}(x - {w_1})} \right]{Y_0}({\lambda _0}z) + } \Bigg. \\[-3pt] & \sum\limits_{n = 1}^\infty \left[ {F_n^{\text{S}}{{\rm {exp}}\;\left[{{\lambda _n}(x - {w_1})}\right]}} \right.+\\[-3pt] & \Bigg.\left. G_n^{\text{S}}{{\rm {exp}}}\;\left[{ - {\lambda _n}(x - {w_1})}\right] \right]{Y_n}({\lambda _n}z) \Bigg\}. \end{split} $

式中： $B_n^{\rm S}$、 $C_n^{\rm S}$、 $F_n^{\rm S}$、 $G_n^{\rm S}$均为未知系数. ${Y_n}({\lambda _n}z) $的表达式为

(12) $ {Y_n}({\lambda _n}z) = M_n^{ - 1/2}\cos\;\left[ {\lambda _n}(z + h)\right],\,n = 0,1,2, \cdots , $

(13) $ \left.\begin{aligned} {M_n} = \left\{ \begin{aligned} &h - d,\quad\quad\;\;\; n = 0; \\ &(h - d)/2,\,\quad n = 1,2,3 \cdots , \\ \end{aligned} \right. \\ {\lambda _n} ={n}/({{h - d}});\quad\;\;\; n = 0,1,2, \cdots . \end{aligned} \right\}$

对于结构物尖角附近流体运动存在奇异性的问题，采用Deng等^[17]的处理方法，在区域间的公共界面上引入正交多项式将速度分布函数近似展开. 速度分布函数表达式为

(14) $U_{{\rm{ - }}{w_2}}^{\rm{S}}(z) = - \frac{{{\rm i}gA}}{\omega }\sum\limits_{r = 0}^\infty {\alpha _r^{\rm{S}}{u_{1r}}(z)} ,$

(15) $U_{{\rm{ - }}{w_1}}^{\rm{S}}(z) = - \frac{{{\rm i}gA}}{\omega }\sum\limits_{s = 0}^\infty {\beta _s^{\rm{S}}{u_{1s}}(z)} ,$

(16) $U_{{w_1}}^{\rm{S}}(z) = - \frac{{{\rm i}gA}}{\omega }\sum\limits_{q = 0}^\infty {\gamma _q^{\rm{S}}{u_{2q}}(z)} ,$

(17) $U_{{w_2}}^{\rm{S}}(z) = - \frac{{{\rm i}gA}}{\omega }\sum\limits_{p = 0}^\infty {\kappa _p^{\rm{S}}{u_{2p}}(z)} .$

式中： $U_{ - {w_2}}^{\rm{S}}(z)$、 $U_{ - {w_1}}^{\rm{S}}(z)$、 $U_{{w_1}}^{\rm{S}}(z)$、 $U_{{w_2}}^{\rm{S}}(z)$分别为相应界面上的水平速度分布函数， $\alpha _q^{\rm{S}}$、 $\beta _s^{\rm{S}}$、 $\gamma _r^{\rm{S}}$、 $\kappa _p^{\rm{S}}$ 均为未知系数. ${u}_{1{r}}(z){\text {、}} {u}_{2{r}}(z)$的表达式分别为

(18) $\begin{split} {u_{1r}}(z) =& \displaystyle\frac{{{{( - 1)}^r}(2r)!\varGamma (1/6)}}{{\pi {{(h - d)}^{{1}/{3}}}\varGamma (2r + 1/3)}}\times \\& {[{(h - d)^2} - {(z + h)^2}]^{ - 1/3}}{C}_{2r}^{1/6}\left(\displaystyle\frac{{z + h}}{{h - d}}\right), \end{split} $

(19) $\begin{split} {u_{2r}}(z) =& \displaystyle\frac{{{{( - 1)}^r}(2r)!\varGamma (1/6)}}{{\pi {{(h - d)}^{{1}/{3}}}\varGamma (2r + 1/3)}}\times \\& {[{(h - d)^2} - {(z + h)^2}]^{ - 1/3}}{ C}_{2r}^{1/6}\left(\displaystyle\frac{{z + h}}{{h - d}}\right). \end{split} $

式中：Γ为Gamma函数，C为盖根堡多项式.

在获得各区域速度势展开式及公共界面上速度分布函数的基础上，应用区域间的速度连续条件及压力连续条件，构建能求解出所有未知系数的线性方程组.

在辐射问题中，气室内液面上存在气压波动，因此区域Ω₃的上边界为非齐次边界条件，即式（2）中 ${p_{\rm c}} \ne 0$. 在区域Ω₃中，除需求解出满足齐次边界条件的通解之外，还需得到1个特解. 与求解衍射问题类似，先求得各区域速度势的展开表达式,然后再将区域间公共界面上的速度分布函数进行近似展开. 最后,应用速度与压力连续条件，构建线性方程组进行数值求解. 其具体求解过程可参阅文献[17].

假定气室内空气可压，Sarmento 等^[18]提出关联线性能量输出(power take-off, PTO)装置与气室内空气质量变化的方程：

(20) $M = - \frac{{{\rm{d}}({\rho _{\rm{c}}}V)}}{{{\rm{d}}t}} = - \rho _{\rm{c}}^0\frac{{{\rm{d}}V}}{{{\rm{d}}t}} - {V_0}\frac{{{\rm{d}}{\rho _{\rm{c}}}}}{{{\rm{d}}t}} = \frac{{\hat KD}}{{\hat N}}{p_{\rm{c}}}.$

式中：M为空气质量，ρ_c为空气密度，V为气室体积，V₀、 $\rho _{\rm{c}}^0 $分别为静止时气室内空气体积和密度， $\hat K$ 为与透平几何形状相关的阻尼系数， $\hat N$ 为涡轮叶片转速，D为涡轮转子直径。为了获得最优转换效率，涡轮转速、转子直径和阻尼系数的影响由可调控的透平参数χ统一表征：

(21) $\chi = \frac{{\hat KD}}{{h\hat N}}{\left( {\frac{g}{h}} \right)^{1/2}}.$

取

(22) ${p_{\rm{c}}} = {\hat p_{\rm{c}}}\rho gA,$

(23) ${\displaystyle\int_{ - {w_1}}^{{w_1}} {\left. {\frac{{\partial \phi _3^{\rm{S}}}}{{\partial z}}} \right|} _{z = 0}}{\rm d}x = {Q^{\rm{S}}}A,$

(24) ${\displaystyle\int_{ - {w_1}}^{{w_1}} {\left. {\frac{{\partial \phi _3^{\rm{R}}}}{{\partial z}}} \right|} _{z = 0}}{\rm d}x = {Q^{\rm{R}}}{p_{\rm{c}}}.$

式中： ${\hat p}_{\rm c} $为无量纲气室压强，Q^S、Q^R分别为绕射和辐射问题下的体积通量.

采用Deng等^[17]的方式，式（20）可以改写为

(25) $\left({Q^{\rm{R}}} + \frac{{{\rm i}\omega {V_0}}}{{\rho _{\rm{c}}^0c_{\rm a}^2}} - \frac{{\hat KD}}{{\rho _{\rm{c}}^{\rm{0}}\hat N}}\right){\hat p_{\rm{c}}} = - \frac{{{Q^{\rm{S}}}}}{{\rho g}}.$

式中： $c_{\rm a} $为声速.

根据Deng等^[17]的研究，1个周期内被转换的波浪能为

(26) ${P_{{\rm{out}}}} = \frac{{\hat KD}}{{2\rho _{\rm{c}}^{\rm{0}}\hat N}}{\left| {{p_{\rm{c}}}} \right|^2} = \frac{{{\rho ^2}{g^2}{A^2}\hat KD}}{{2\rho _{\rm{c}}^{\rm{0}}\hat N}}{\left| {{{\hat p}_{\rm{c}}}} \right|^2}.$

将其除以入射波的波能流 ${P_{{\rm{in}}}} = \rho g{A^2}{C_{\rm g}}$，其中C_g为波群速度，可得无量纲的波能转换效率：

(27) $E = \frac{{\rho h{{\left( {gh} \right)}^{1/2}}}}{{\rho _{\text{c}}^{\text{0}}{C_{\rm g}}}}\frac{{\hat KD}}{{h\hat N}}{\left( {\frac{g}{h}} \right)^{1/2}}{\left| {{{\hat p}_{\text{c}}}} \right|^2}.$

考虑二维不可压缩流体，其流动控制方程包括连续方程和N-S方程：

(28) $\frac{{\partial u}}{{\partial x}} + \frac{{\partial w}}{{\partial z}} = 0,$

(29) $\frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} + w\frac{{\partial u}}{{\partial z}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial x}} + \nu \left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right),$

(30) $\frac{{\partial w}}{{\partial t}} + u\frac{{\partial w}}{{\partial x}} + w\frac{{\partial w}}{{\partial z}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial z}} + \nu \left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}} + \frac{{{\partial ^2}w}}{{\partial {z^2}}}} \right) - g.$

式中：u、w分别为x轴和z轴2个方向上的速度分量，p为流体压强，ν为流体运动黏滞系数.

采用VOF方法追踪水气交界面的波动，a_q为体积分数函数，其定义为单元内第q相流体所占体积与该单元的体积之比. a_q=1表示第q相流体占据整个流体单元；a_q=0表示在整个流体单元内不含第q相流体. 若a_q $\in $ (0, 1.0)，则表示该单元为自由液面单元. 每个单元内所有相的体积分数之和为1，即

(31) $\sum\limits_{q = 0}^1 {{a_q}} = 1.$

式中：q=0为空气相，q=1为水相.

此外，a_q还应满足体积分数连续方程：

(32) $\frac{{\partial {a_q}}}{{\partial t}} + \frac{{\partial \left( {u{a_q}} \right)}}{{\partial x}} + \frac{{\partial \left( {v{a_q}} \right)}}{{\partial y}} = 0.$

借助FLUENT软件的用户自定义函数（user-defined function, UDF），在入射边界给定线性波的速度表达式及相应波面函数生成入射波. 入射波的波浪要素有自由液面高度η，水平方向速度u，垂直方向速度w，可以分别表示为

(33) $\eta = A\cos\; (kx - \omega t),$

(34) $u = A\omega \frac{{\cosh\; [k\left( {z + h} \right)]}}{{\sinh\; (kh)}}\cos\; \left( {kx - \omega t} \right),$

(35) $w = A\omega \frac{{\sinh\; [k\left( {z + h} \right)]}}{{\sinh\; (kh)}}\sin\; \left( {kx - \omega t} \right).$

由于压力出口边界条件不能保证波浪完全透射，为避免出口边界处反射波的形成，在水槽末端设置阻尼消波区，其实现方法是在动量方程中添加1个消波阻尼源项，表达式为

(36) ${S_i} = \beta \frac{{ {x - {x_1}} }}{{\left| {{x_1} - {x_2}} \right|}}{v_i}.$

式中：x₁为消波区的起点，x₂为消波区的终点，β为消波系数，v_i为速度分量.

边界条件的设置如下：

1）上边界与右边界设为压力出口边界条件；左边界设为速度入口边界条件；下边界与OWC装置表面设为固壁边界条件.

2）控制方程中的扩散项采用中心差分格式离散，动量方程和压力方程的对流项分别采用二阶迎风格式和Body Force Weighted格式，采用PISO算法求解压力与速度耦合方程.

OWC气室上方孔口处存在1个线性空气透平，气室内的压强与孔口处的空气流速呈线性关系. 通过引入多孔阶跃层来代替该线性透平，阶跃层满足达西定律：

(37) $\Delta p = -{{\mu \Delta m}}v/{\alpha }.$

式中：Δp为气室内外压强差，Δm为阶跃层厚度，μ为空气动力黏滞系数， $\alpha $为渗透系数. 针对线性透平情况，Sarmento等^[18]引入无量纲参数C，表达式为

(38) $C = \frac{{{\rho _{\rm{w}}}}}{{{\rho _{\rm a}}}}{\left(\frac{g}{{{w_1}}}\right)^{{1}/{2}}}\frac{M}{{{p_{\rm{c}}} - {p_{\rm a}}}}.$

式中：ρ_w为水的密度；ρ_a为空气密度；M=Lvρ_a，L为PTO处孔径；p_a为大气压强. α可以表示为

(39) $\alpha = \frac{{\mu \Delta mC}}{{\rho b{{(g/{w_1})}^{{1}/{2}}}}}.$

式中：b为气室宽度.

此外，理论模型中的χ与C有如下转换关系：

(40) $\chi {\rm{ = }}\frac{{C{{\left( {g/h} \right)}^{{1}/{2}}}{\rho _{\rm a}}}}{{h{{(g/{w_1})}^{{1}/{2}}}\rho }}.$

对于入射波波流能的计算，通常采用基于势流理论得到的理论公式（41）进行求解. 考虑到理论解与实际情况存在一定的偏差，为了准确计算OWC装置的波能捕获效率，本研究采用断面积分法，根据式（42）直接计算入射波传播过程中监控断面处的波能通量E，并将监控得到的波能流结果按式（43）取周期平均，记为P_w.

(41) $E = \rho g{A^2}{C_{\rm g}}/2,$

(42) $\begin{gathered} E = \int_{ - h}^\eta {{p_{_{\rm{D}}}}u{\rm d}z} = \int_{ - h}^\eta {(p + \rho gz)u{\rm d}z} =\\ \int_{ - h}^{{Z_0}} { \varphi{_{{\rm{w}}}}(p + \rho gz)u{\rm d}z} . \\ \end{gathered} $

式中：p_D为动压压强，p为静压压强， $ \varphi $_w为网格单元内水的体积分数.

(43) ${P_{\rm{w}}} = \frac{1}{T}\int\limits_t^{t + {{T}}} E {\rm d}t.$

式中：P_w为监控断面处波能通量的周期平均值，T为周期.

如图2所示为在入射波长为12.9 m的情况下，断面处(距离入口64.5 m)单宽波峰线的波能流E的时间历程曲线. 从第16T开始断面处监控得到的曲线趋于稳定，波能流峰值略有起伏. 取16T~30T即若干个周期内的波能流进行平均，得到入射波波能. 如图3所示为水槽中入射波长为12.9 m的情况下，通过断面积分法得到的PTO处的波能捕获值，记为C_p，采取周期平均的方式进行处理.

(44) ${P_{{\rm t}}} = \frac{1}{T}\int\limits_t^{t + {{T}}} {\int\limits_l {{p_{\rm{c}}}{v_{\rm{n}}}{\rm d}l{\rm d}t} } .$

式中：P_t为捕获的能量，p_c为气室内的压强，v_n为孔口处的空气流速. OWC装置的波能捕获效率的表达式为

(45) $\varepsilon {\rm{ = }}{{{P_{\rm{t}}}}}/{{{P_{\rm{w}}}}}.$

中国海波能功率密度呈现南高北低的状态，其中南海波能资源最为丰富，波高为0.4~3.0 m，周期为3.1~23.0 s^[19]. 以南海海域波况为依据，进行OWC装置工作效率的研究. 根据弗罗德相似准则（几何比尺为1∶10），选取波高为0.04 m，周期为2.51~4.76 s的入射波进行数值模型验证，分析波浪作用下OWC装置的波能转换效率，并将所得数值结果与解析解进行对比，波况如表1所示. 表中，λ为波长.

如图4所示，水槽总长为8L_o，水槽末端设置消波区，长度为L_o，OWC气室宽度为2.40 m。气室高度为1.00 m，前后墙吃水深度为1.75 m，墙体厚度为0.35 m，PTO处孔径为0.35 m. 经网格收敛性验证，单位波长内划分80个网格，自由液面处网格加密，单位波高内划分10个网格，时间步长 ${{{\rm \Delta} t = T}}/{\rm{400}}.$ 如表1所示，10组工况下，气室宽度和波长的比值为0.071~0.245，在装置吃水深度相同的情况下，依据Luo等^[20]的结果可知，该区间内PTO出口处最优阻尼系数C为1.0~5.0，为在所考虑的波浪频率范围内获得的最优转换效率. 在以下计算中C=4.3，为装置共振频率附近的最优值.

考虑到在工况1~10不同周期入射波的作用下，气室内波面的升降规律基本相似，以工况6为例，研究1个完整波周期内OWC气室内的水面波动及附近流场情况. 如图5所示，当波峰到达OWC装置前方时，气室内水位开始抬升，水体流入OWC气室；当波谷到来时，气室内水位下降，水体流出气室. 气室内部水位平整，呈现活塞式的周期性波动，运动周期与入射波周期保持一致. 入射波与离岸式OWC装置相互作用，水动力过程较复杂，入射波波能存在被反射、吸收、透射以及耗散的情况. 如图6所示，在OWC前方3.0 m、后方3.7 m处，以及内部分别设置水位监测点，分别记为监测点1、3、2.

以工况6为例，如图7所示为各测点处的水位波动情况. 由图可以看出，OWC前方水面波动最大，入射波在OWC前墙发生部分反射，形成驻波. OWC内部水面起伏较大，波幅超过0.02 m，出现共振现象. OWC后方水面波动幅值相对于入射波幅较小，透射波波能较小，入射波在经过OWC装置前后能量衰减较大. 如图8所示为工况6下，OWC气室内的压强p_c和PTO出口流速 $ \widetilde v$的时间历程曲线. 在t=4T时波浪传播到OWC装置附近，从t=13T开始气室内压强与流速呈现稳定的正弦型周期波动，且两者相位保持一致. 这是由于在OWC数值模型中考虑了线性透平的情况且忽略了空气的可压缩性.

考虑到入射波长对OWC装置波能吸收效率的影响，如图9所示为在波高为0.04 m，波长分别为12.9、16.5、18.9、33.7 m的情况下，OWC气室内压强p_c，PTO出口流速 $\widetilde v$以及装置波能捕获量C_p的时间历程曲线. 当入射波长分别为12.9、16.5、18.9、33.7 m时，对应的气室内压强分别为22、41、48、34 Pa，PTO处流速分别为0.019、0.036、0.042、0.029 m/s. 当入射波长为19.8 m时，OWC气室内压强与流速的振荡幅值均最大，捕获的波浪能最多.

为了验证OWC数值模型的准确性，将10组波况下得到的波能吸收效率与解析解进行对比. 在数值模拟中选取的相对波高较小，因此水体黏性耗散和波浪非线性效应对结果的影响可以忽略不计. 由图10可以看出，在工况6下，当入射波周期为3.48 s时，OWC装置的能量捕获效率ε最高，达到约46%，其波频率接近共振频率，约为1.8 Hz. 当入射波频率远离结构物共振频率时，波能的吸收效率下降. 由图10可以看出，数值结果与解析解吻合良好.

进一步分析能量的守恒情况. 采用Goda两点法分离入射波与反射波波高，如图11所示为不同波浪频率下得到的反射系数C_r和透射系数C_t，采用B样条曲线进行数据拟合. 可以看出，入射波长的增大导致反射系数减小，大部分入射波能透过前墙传入气室，同时透射系数也逐渐增大，更多的波能透射出去. 采用式（41）计算反射波与透射波波能流. 根据能量守恒原理可知，入射波波能流E_total=反射波波能流E_r+透射波波能流E_t+OWC装置捕获的能量E_ab+能量耗散E_d. 如图12所示为不同波浪频率下各部分能量占入射波波能流的比值P_r，图中的曲线为B样条曲线，其中能量耗散占比较小，这是由于墙体厚度相对于波长较小. 如图5所示也反映出墙体附近没有明显的涡旋存在. 若忽略能量的耗散，各部分能量总和约等于入射波波能流，进一步验证了模型的有效性，所建立的OWC数值模型能够准确模拟波浪与OWC装置的相互作用.

在水深一定的情况下，OWC结构尺寸是影响其工作效率的主要因素，现考虑墙体吃水深度、墙体厚度以及气室宽度对波能转换效率的影响. 计算域水深为10 m，入射波波高为0.04 m.

在OWC墙体厚度为0.35 m，气室宽度为2.80 m，装置吃水深度分别为1.50、1.75、2.10 m时，不同频率下OWC装置的能量转换效率如图13所示. 图中，d_a、d_n分别为d的解析解和数值解. 由图可以看出在当前计算范围内，改变吃水深度对低频区的波能转换效率的影响较小，而在高频区，吃水深度的增大导致转换效率明显下降. 这是因为在低频区吃水深度相对波长较小，能量转换效率对吃水深度的变化不敏感，而在高频区吃水深度相对波长较大，能量转换效率对吃水深度的变化较敏感. 同时，可以看出吃水深度的增大，会导致装置高效频率带宽变窄，峰值向低频区移动. 这是由于吃水深度的增加，使得气室内水体质量增大，从而导致OWC的共振频率减小.

在OWC装置吃水深度为1.75 m，气室宽度为2.80 m，墙体厚度分别为0.35、0.60、0.90 m的情况下，不同频率下OWC装置的能量转换效率如图14所示. 图中，(w₂−w₁)^a、(w₂−w₁)ⁿ分别为w₂−w₁的解析解和数值解. 可以看出，墙体厚度对低频区波能转换效率的影响较小，在高频区，墙体厚度的增大导致波能转换效率明显下降. 此外，墙体厚度的增大，导致OWC波能转换效率峰值向低频区移动.

在OWC装置吃水深度为1.75 m，前后墙厚度为0.35 m，气室宽度分别为2.0、2.4、2.8 m时，不同频率下OWC装置的能量转换效率如图15所示. 图中， $2w^{\rm a}_1 $、 $2w_1^{\rm n} $分别为2w的解析解和数值解. 结果表明，气室宽度的增大导致装置高效频率带宽变宽，波能转换曲线的有效频率工作区间变宽，峰值向低频区移动，波能的转换效率在低频区明显提高，在高频区略有下降.

在线性波理论下，波能转换效率与入射波波高无关，而在实际情况下Luo等^[15]、Wang等^[16]以及Elhanafi等^[21]的工作表明，入射波波高增大导致的非线性效应会对OWC装置的工作效率产生重要影响.

如图16所示为4组不同波幅H下波能转换效率随波浪频率的变化，由图可以看出波幅的增大会导致波能转换效率下降，在共振频率附近尤为明显. 当波长为18.9 m，波幅分别为0.04、0.16 m时，OWC装置附近的涡量Ω及流场分布情况如图17所示. 由图可以看出，在H=0.04 m的情况下，墙体底部涡量不明显，涡量值较小，而当H=0.16 m时，墙体底部附近流速较大，涡量明显增强，即波幅的增大导致墙体附近黏性效应增强，能量耗散加大.

此外，通过FFT变换得到不同波幅下的入射波频谱，如图18所示. 由图可以看出，当H=0.16 m时非线性效应不明显，无高阶波的出现，而当H=0.48 m时非线性效应增强，在频率0.574 Hz附近出现二阶波，即入射波波能从低频波向高频波转移. 由于高频波波能难以传入OWC装置内部，OWC装置波能转换效率下降.

从理论和数值上研究OWC装置吃水深度、气室宽度以及墙体厚度对波能转换效率的影响，通过对比发现数值结果与解析解吻合较好. 由结果可知，墙体吃水深度的增加会导致波能转换效率的高效频率带宽变窄，峰值向低频区移动；墙体厚度的增加对低频区波能转换效率的影响较小，但会导致高频区波能转换效率明显下降；气室宽度的增大会导致OWC装置的共振频率向低频区移动，同时在低频区，波能转换效率明显提高，而在高频区略有下降. 此外，波幅的增大带来的黏性及非线性效应会导致OWC装置转换效率下降，在共振频率附近尤为明显. 所建立的OWC数值模型能够准确模拟波与离岸式OWC装置的相互作用，为进一步研究海底地形变化对OWC装置波能转换效率的影响奠定基础.