0 引　言

随着对电力的大量需求和化石燃料的快速消耗，迫切需要发展可再生能源，世界范围内，新一轮的能源革命正在进行。2021年，全球可再生能源新增装机容量占全部新增装机容量的81%，传统的化石能源正加速向可再生能源过度^[1]。太阳能是一种清洁的可再生能源，经过多年的发展技术最为成熟，且具备大规模开发的条件。光伏系统是将太阳能转化为电能的最常见方式。

海域面积宽广，日照充足且无遮挡，布局光伏有天然的优势，海上漂浮式光伏成为新能源发展的新战场。Yan等^[2]介绍了一种专门适用于海洋环境的 FPV 平台的新型模块化设计，研究了多体 FPV 系统在极端海况和运行海况下的流体动力学和耦合效应，研究分析了铰接和固定连接边界下的海上光伏系统平台的连接器强度和系泊系统受力。Magkouris等^[3]提出一种基于边界元法的双船体结构的漂浮式光伏，并研究了其在波浪作用下，该浮式光伏结构在不同水深的水动力特性。

近年来，离散元方法得到广泛关注，尤其是在对水动力分析领域发挥着重要作用。Zhu等^[4]使用球形离散元方法构建浮式结构的水平冰和系泊缆，模拟水平冰与半潜式结构相互作用过程中的冰荷载和系泊力。翟必垚等^[5]将离散元方法与二维水动力学相耦合，建立了河冰的动力数值模型，以模拟河冰在水面的运动姿态以及冰坝形成的动力过程。刘璐等^[6]采用离散元方法分析海冰与海洋结构作用中的冰荷载，分别计算了船舶结构和海洋平台结构的冰荷载特点。王祥等^[7]采用扩展多面体离散元方法建立冰船动力作用模型，在此基础上给出了对冰船的六自由度运动响应和冰阻力结果。Ji等^{[8 − 9]}利用GPU并行技术开发了一种粘结球形离散元方法，以模拟水平冰的破裂和运动，他们还通过单轴压缩和三点弯曲试验的数值模拟，定量研究了DEM计算参数与海冰宏观力学性能之间的关系。这些研究表明，DEM适用于模拟浮体的水动力特性分析。

本文使用球形离散元法来构建海上漂浮式光伏电站，对不同波高、不同周期和不同结构直径的海上漂浮式光伏电站系统的水动力特性进行数值计算。通过对计算结果的比较及分析，探讨了波浪要素和结构直径对海上漂浮式光伏电站系统的运动响应和受力特性的影响。

1 研究对象及其简化

本文对水深h=30 m，波高H=4～7 m，周期T=5.4～8.9 s的波浪作用下，海上漂浮式薄膜光伏电站系泊系统的水动力特性进行数值模拟分析，如图1所示。这一范围的波浪具有非线性，可以采用Stokes二阶波浪理论^[10]进行计算。

采用球形离散元方法来构建海上漂浮式薄膜光伏电站系统的数值模型，将海上漂浮式薄膜光伏电站系统视为刚体，球形元件之间不计算接触力，通过不同的单元配置和计算参数来定义球形元件的不同的力学性能，外力作用在球形元件上，根据数值模型的结构求出质心坐标，然后求出每个球形元件对质心的外力矩和外力，最后对海上漂浮式薄膜光伏电站系统进行运动求解。如图2(a)所示。薄膜式光伏电站由环形浮力环和一层薄膜组成，光伏组件通过导轨安装在薄膜上，为了简化计算过程，将薄膜和光伏板的重量计入导轨上，并建立数值模型，将环形浮力环和导轨视为细长体，如图2(b)所示。

2 计算理论 2.1 波浪场理论

根据Stokes二阶波浪理论，波面升高和水质点运动速度为：

$ \eta = {{A}}\left[ { - \frac{{Ak}}{{2{\rm sh}2kh}} + {\rm cos}\theta + \frac{{Ak}}{4}\frac{{{\rm ch}kh\left( {2{\rm c{h}^2}kh + 1} \right)}}{{{\rm s{h}^4}kh}}{\rm cos}2\theta } \right] ，$

(1)

$ {V_x} = A\omega \left[ {\frac{{{\rm ch}k\left( {z + h} \right)}}{{{\rm sh}kh}}{\rm cos}\theta + \frac{3}{4}Ak\frac{{{\rm ch}2k\left( {z + h} \right)}}{{{\rm s{h}^4}kh}}{\rm cos}2\theta } \right]，$

(2)

$ {V_z} = A\omega \left[\frac{{{\rm sh}k(z + h)}}{{{\rm sh}kh}}\sin \theta + \frac{3}{4}Ak\frac{{{\rm sh}2k(z + h)}}{{{\rm s{h}^4}kh}}{\rm sin}2\theta \right] 。$

(3)

波浪水质点的加速度通过对式（2）和式（3）求导得到：

$ {a_x} = A{\omega ^2}\left[ {\frac{{{\rm ch}k\left( {z + h} \right)}}{{{\rm sh}kh}}{\rm sin}\theta + \frac{3}{2}Ak\frac{{{\rm ch}2k\left( {z + h} \right)}}{{{\rm s{h}^4}kh}}{\rm sin}2\theta } \right] ，$

(4)

$ {a_z} = - A{\omega ^2}\left[ {\frac{{{\rm sh}k\left( {z + h} \right)}}{{{\rm sh}kh}}{\rm cos}\theta + \frac{3}{2}Ak\frac{{{\rm sh}2k\left( {z + h} \right)}}{{{\rm s{h}^4}kh}}{\rm cos}2\theta } \right] 。$

(5)

式中：$\theta = kx - \omega t$；${V_x}$和$ {V_z} $分别为波浪水质点的水平和竖直方向速度；${a_x}$和${a_z}$分别为水平和竖直方向的加速度；$A$为波幅；$\omega $为波浪频率；$k$为波数；$h$为水深；$t$为时间；$x$、$z$为波浪水质点的位置。

2.2 构建运动方程

海上漂浮式光伏电站系统在波浪作用下受到的外力有重力、浮力、波浪力以及系泊力，波浪力主要由两部分组成：一是水质点运动速度引起的拖曳力$ {\stackrel{\rightharpoonup }{F}}_{D} $，另一个是水质点加速度引起的惯性力$ {\stackrel{\rightharpoonup }{F}}_{I} $。当流体质点和浮体存在相对运动时，波浪力$ {\stackrel{\rightharpoonup }{F}}_{W} $采用Morison^[11]方程计算为：

$ {\stackrel{\rightharpoonup }{F}}_{D}=\rho {C}_{D}S\frac{{\stackrel{\rightharpoonup }{V}}_{R}\left|{\stackrel{\rightharpoonup }{V}}_{R}\right|}{2} ，$

(6)

$ {\stackrel{\rightharpoonup }{F}}_{I}=\rho {K}_{m}{V}_{f}\frac{\partial \left(V-\dot{\xi }\right)}{\partial t}+\rho {V}_{f}\frac{\partial V}{\partial t} ，$

(7)

$ {\stackrel{\rightharpoonup }{F}}_{W}={\stackrel{\rightharpoonup }{F}}_{D}+{\stackrel{\rightharpoonup }{F}}_{I} 。$

(8)

式中：$ {\stackrel{\rightharpoonup }{F}}_{W} $为波浪力；$ {\stackrel{\rightharpoonup }{F}}_{D} $为拖曳力；$ {\stackrel{\rightharpoonup }{F}}_{I} $为惯性力；$ \rho $为流体密度，$ \rho $取值1035 kg/m³；$ {C_D} $为阻力系数取值1.2；$ S $为颗粒在流速方向的投影面积；$ {\stackrel{\rightharpoonup }{V}}_{R}=\stackrel{\rightharpoonup }{V}-\dot{\xi } $，$ \stackrel{\rightharpoonup }{V} $为波浪场速度，$ \dot \xi $为颗粒的运动速度，$ {\stackrel{\rightharpoonup }{V}}_{R} $为颗粒与水质点的相对速度；$ {K_m} $为质量系数，取值1.0；$ {V_f} $为浸水体积。

颗粒在流速方向的投影面积$ S $和浸水体积$ {V_f} $随着运动不断改变，图3为颗粒截面图。其在各个方向的投影面积$ S $和浸水体积$ {V_f} $为：

入水深度$\Delta h$为：

$ \Delta h=\eta (x,z,t)-(z\left(t\right)-R)，$

(9)

$\alpha $为Z轴的轴线至水面夹角的二倍：

$ \alpha = 2{\rm{cos}}{^{ - 1}}\left( {\frac{{\left| {D/2 - \Delta h} \right|}}{{D/2}}} \right) ，$

(10)

$ {S_x} = {S_y}\left\{ {\begin{array}{l} {{\text{π}} {R^2} - 1/2{R^2}\left( {\theta - {\rm sin}\theta } \right)}，{0 < \Delta h \leqslant R}，\\ {1/2{R^2}\left( {\theta - {\rm sin}\theta } \right)}，{R < \Delta h \leqslant 2R}，\\ {{\text{π}} {R^2}}，{2R < \Delta h}，\end{array}} \right. $

(11)

$ {S}_{\text{z}}\left\{\begin{array}{l}{\text{π}} \left[{R}^{2}-{\left(R-\Delta h\right)}^{2}\right]，0 < \Delta h < R，\\ {\text{π}} {R}^{2} ，R\leqslant\Delta h\leqslant 2R，\end{array}\right. $

(12)

$\left\{\begin{array}{l} {V}_{f} = \displaystyle\frac{{\text{π}} \Delta h}{6} \left(3({R}^{2} - {\left(R - \Delta h\right)}^{2}\right) + \Delta {h}^{2})，0 < \Delta h \leqslant 2R，\\ {V}_{f}=\displaystyle\frac{4}{3}{\text{π}} {R}^{3}，\Delta h > 2R。\end{array}\right. $

(13)

系泊力是由缆绳的伸长产生的，与缆绳的固有参数有关，表示为：

$ {F_T} = {C_p}{d^2}{\left( {\frac{{l - {l_0}}}{{{l_0}}}} \right)^n} 。$

(14)

式中：${F_T}$为张力；${l_0}$为缆绳原始长度；$l$为缆绳变形后长度；$d$为缆绳直径；${C_p}$为缆绳的弹性系数；$n$为指数。

将海上漂浮式薄膜光伏电站系统所在的空间坐标系分为整体坐标系、以质心为坐标原点的整体坐标系和局部坐标系。如图4所示。整体坐标系为整个光伏电站系统的坐标系，用$ G $表示。以质心为坐标原点的整体坐标系的坐标原点固定在单元质心上，随着质心的运动而运动，3个坐标轴与整体坐标系平行，简称为$ GL $坐标系。局部坐标系的坐标原点同样固定在质心上，随着质心运动而运动，其3个坐标轴不一定与整体坐标系平行，光伏电站的转动惯量即是在该坐标系下求得，采用$ B $表示。

对于光伏电站系统的运动求解，视为局部坐标系在整体坐标系下的运动。该运动分解为平动和转动。平动方程在光伏电站的质心上进行求解，局部坐标系随之发生平动；转动视为局部坐标系绕质心的转动，本质为局部坐标系和$ GL $坐标系之间的坐标转换过程。对于局部坐标系下的坐标$ {e^{{B}}} $和$ GL $坐标系下的坐标$ {e^{GL}} $，其转换关系为：

$ {e^{{B}}} = \bar {\boldsymbol A}{e^{{{GL}}}}。$

(15)

式中：$ B $为固定在单元上的局部坐标系；$ \bar {\boldsymbol A} $为转换矩阵，其$ {\bar {\boldsymbol A}^{ - 1}} = {\bar {\boldsymbol A}^{\text T}} $。

光伏电站的平动方程为：

$ m\frac{{{{\rm d}^2}x}}{{{\rm d}{t^2}}} = F 。$

(16)

式中：$m$为质量；$F$为所受合力。

光伏电站的转动方程为：

$ \left\{ {\begin{array}{*{20}{c}} {{I_{xx}}\dot \omega _x^B = M_x^B + \left( {{I_{yy}} - {I_{zz}}} \right)\omega _y^B\omega _z^B}，\\ {{I_{yy}}\dot \omega _y^B = M_y^B + \left( {{I_{zz}} - {I_{xx}}} \right)\omega _z^B\omega _x^B} ，\\ {{I_{zz}}\dot \omega _z^B = M_z^B + \left( {{I_{xx}} - {I_{yy}}} \right)\omega _x^B\omega _y^B} 。\end{array}} \right. $

(17)

式中：$ B $为固定在光伏电站上的局部坐标系；$ \Bigr\{ M_x^B, M_y^B,M_z^B \Bigr\}^{\text T} = M_{}^B $$ {\left\{{ {M}}_{ {x}}^{ {B}}，{ {M}}_{ {y}}^{ {B}}，{ {M}}_{ {z}}^{ {B}}\right\}}^{ {\text T}}={ {M}}^{ {B}} $是局部坐标系下所受力矩；$ {\left\{ {\omega _x^B,\omega _y^B,\omega _z^B} \right\}^{\text T}} = \omega _{}^B $$ {\left\{{ {\omega }}_{ {x}}^{ {B}}，{ {\omega }}_{ {y}}^{ {B}}，{ {\omega }}_{ {z}}^{ {B}}\right\}}^{ {\text T}}= { {\omega }}^{ {B}} $是局部坐标系下的角速度。

采用四元数的方法对光伏电站的转动方程进行求解。

3 对比验证

为了验证本文采用的方法对海上漂浮式薄膜光伏电站系统水动力特性分析的准确性，分别对浮体的运动和受力进行验证。对于浮体运动的验证采用文献[12]，计算工况为水深h=20 m，波高H=4 m，周期T=7.16 s，波浪理论为Stokes二阶波，验证结果对比图如图5所示。对于系泊力的验证采用文献[13]方法，计算工况为水深15 m，波高为4.2～7 m，周期为7.2 s，波浪理论为线性规则波，验证结果对比如图6所示。

从2个对比验证的结果图可知，本文的数值结果与所对比文献的数值结果基本吻合。因此可以确定本文所使用的数值方法适用于海上漂浮式薄膜光伏电站的水动力特性分析。

4 数值计算

为探讨海上漂浮式薄膜光伏电站系统在波浪作用下的质心位移、波浪力和系泊力的变化，对其分别在不同波高、不同周期和不同结构直径的3种情况下进行了水动力特性数值模拟计算。

4.1 波高对海上漂浮式薄膜光伏电站水动力特性的影响

计算工况如表1所示。从图7～图9中可以看出，在相同周期不同波高情况下，可以得出以下结论：

1） 随着波高的增大，海上漂浮式薄膜光伏电站系统的质心最大水平位移和最大垂直位移呈现出逐渐增大的趋势。这是由于较高的波高代表波浪的振幅越大，因此位移的幅度也会增大，这意味着波的能量和强度的增大。

2）海上漂浮式薄膜光伏电站系统所受的水平波浪力大于垂直波浪力，随着波高的增大，其所受的水平波浪力和垂直波浪力呈现出逐渐增大的趋势，波高越大，其所受波浪力越大。这是因为较高的波高会导致更大的水平振幅，从而增大水平波浪力对结构的影响，而垂直波浪力的变化较小，因为波浪在垂直方向上的振动幅度相对较小。

3）海上漂浮式薄膜光伏电站系统所受的迎浪侧系泊力大于背浪侧系泊力，随着波高的增大，其所受的迎浪侧系泊力呈现逐渐增大的趋势，背浪侧系泊力呈现出先增大后减小再增大的趋势。这是因为对于迎浪侧系泊力来说，更大的波高会引起更强的波浪冲击力，从而使得迎浪侧系泊力更大，但对于背浪侧系泊力，波浪力在结构的背面会受到波浪的遮挡，从而减小，但随着波浪的继续传递，背浪侧系泊力会再次增大。

4.2 周期对海上漂浮式薄膜光伏电站水动力特性的影响

计算工况如表2所示。从图10～图12中可知，在相同波高不同波浪周期周期的情况下的以下结论：

1） 随着波浪周期的增大，海上漂浮式薄膜光伏电站系统的质心最大水平位移呈现出先减小后增大再减小再增大的趋势，而质心最大垂直位移呈现出一直增大的趋势，周期越大，最大水平位移越小，最大垂直位移越大。造成这一现象的原因为较长周期的波浪会导致结构在波浪中的频率减小，从而减小水平位移的幅度。然而，当周期继续增大时，质心的水平位移再次增大，因为它开始与波浪的周期发生共振。

2）海上漂浮式薄膜光伏电站系统所受水平波浪力大于垂直波浪力，随着波浪周期的不断增大，其所受水平波浪力呈现出先增大再减小的趋势，当波高H=4.0 m和波高H=5.0 m时，其所受垂直波浪力随着周期的增大逐渐减小，而当波高H=5.0 m和H=5.5 m时，其所受垂直波浪力随着周期的增大逐渐增大。这是因为波浪与海上漂浮式薄膜光伏电站之间的相互作用在一定周期范围内引起共振现象，初始阶段，波浪的周期与结构的自然周期相匹配，导致水平波浪力增大。然而，当波浪周期与结构的自然周期不再匹配时，水平波浪力会减小。对于所受垂直波浪力呈现出2种不同趋势，这是因为较小波高的波浪在较短周期下施加的力较大，但随着周期的增大，垂直波浪力会减小。较大波高的波浪在较长周期下产生更大的冲击力，因此垂直波浪力会逐渐增大。

3） 海上漂浮式薄膜光伏电站系统所受迎浪测系泊力大于背浪侧系泊力，随着波浪周期的不断增大，其所受迎浪侧系泊力呈现出先增大再减小的趋势，其所受背浪侧系泊力呈现出先增大后减小再增大的趋势。造成这一现象的原因为海上漂浮式薄膜光伏电站与波浪产生了共振效应，导致迎浪侧系泊力呈现出这样的趋势。而背浪侧所受系泊力随着周期的增大而增大，因为较长周期的波浪在单位时间内对系统施加更多力量，这导致背浪侧的系泊力增大。然而，当周期进一步增大时，系泊系统内的阻尼效应和非线性效应开始显著起作用，以抵消较大周期波浪的影响，因此系泊力会减小。这种减小导致背浪侧系泊力呈现先增大后减小再增大的趋势。

4.3 结构直径对海上漂浮式薄膜光伏电站水动力特性的影响

计算工况如表3所示。从图13～图15可知，在相同波高和波浪周期，不同结构直径的情况下的以下结论：

1） 在波高H＜5.5 m，海上漂浮式薄膜光伏电站系统在结构直径D=50 m时，质心最大水平位移小于其结构直径D=30 m时的质心最大水平位移，在波高H > 5.5 m，其结构直径D=50 m时，质心的最大水平位移大于结构直径D=30 m时的质心最大水平位移，在其结构直径D=50 m时，质心的最大垂直位移都小于结构直径D=30 m时的质心最大垂直位移。这是因为较小波高的波浪施加在较大的直径结构上的力的响应小，而较小的直径结构受到的力的响应大。而较大波高的波浪会施加更大的力，尤其是在较大直径的结构上，从而导致水平位移增大。垂直位移受到垂直波浪力的影响，结构直径越大，结构所受的垂直波浪力越大，位移也就越大。

2）海上漂浮式薄膜光伏电站系统所受水平波浪力大于垂直波浪力，在结构直径D=50 m时，海上漂浮式薄膜光伏电站系统所受的水平波浪力和垂直波浪力都大于结构直径D=30 m时系统所受的波浪力。造成这一现象的原因为较大直径的结构具有更大的表面积，导致结构在水中的横截面积更大，这增加了结构所暴露于波浪力下的区域，从而受到更大的水平波浪力；垂直波浪力通常与波浪振幅相关，而波浪振幅受到波浪的波高影响。当波高相同的情况下，较大直径的结构通常会受到更大的垂直波浪力，因为更大的结构会暴露在波浪的上下运动中，导致更大的垂直力。

3）海上漂浮式薄膜光伏电站系统所受迎浪测系泊力大于背浪侧系泊力，其在结构直径D=50 m时所受迎浪侧系泊力大于结构直径D=30 m时所受迎浪侧系泊力，其在结构直径D=50 m时所受背浪侧系泊力总体趋势为大于结构直径D=30 m所受背浪侧系泊力，但在波高H=5.0 m和H=5.5 m时，其在结构直径D=50 m时背浪侧系泊力会小于结构直径D=30 m时所受背浪侧系泊力。造成这一现象的原因为较大直径的结构在波浪的迎浪侧受到更大的波浪力，因为更大的结构暴露在波浪面上的面积更大，这导致迎浪侧系泊力更强大。同样，较大直径的结构在背浪侧受到更大的波浪力，因为它们也具有更大的阻力，更多的结构暴露在波浪的背浪侧，从而受到更大的背浪侧系泊力。在波高H=5.0 m和H=5.5 m时，海上漂浮式薄膜光伏电站系统在结构直径D=50 m时所受背浪侧系泊力小于结构直径D=30 m时所受背浪侧系泊力，是因为在结构直径D=50 m时结构的倾斜度小于结构直径D=30 m时结构的倾斜度。

5 结　语

1） 与文献[12 − 13]数值结果对比，验证结果吻合，确定了本文所采用的数值方法在对海上漂浮式薄膜光伏电站在波浪作用下的水动力分析的适用性。

2） 随着波高增大，海上漂浮式薄膜光伏电站系统的水平和垂直位移逐渐增大，其水平和垂直波浪力逐渐增大，其迎浪侧系泊力逐渐增大，背浪侧系泊力呈现先增大后减小再增大的趋势，总体趋势为三者随着波高增大逐渐增大。

3） 随着周期增大，海上漂浮式薄膜光伏电站系统的水平位移呈现先减小后增大再减小的趋势，垂直位移逐渐增大，其水平波浪力呈现先增大后减小的趋势，在波高H=4.0 m和H=5.0 m时，其垂直波浪力逐渐减小，在波高H=5.0 m和H=5.5 m时，其垂直波浪力逐渐增大，其迎浪侧系泊力呈现先增大后减小的趋势，背浪侧系泊力呈现先增大后减小再增大的趋势。

4） 不同结构直径条件下，当海上漂浮式薄膜光伏电站系统结构直径D=50 m时，在波高小于5.5 m时其位移小，在波高大于5.5 m时其位移大，其所受水平和垂直波浪力大，其所受迎浪侧和背浪侧系泊力大，总体趋势为三者随着结构直径增大逐渐增大。

5）在本文所给的各种工况下，海上漂浮式光伏电站系统所受的水平波浪力大于垂直波浪力，迎浪测系泊力大于背浪侧系泊力。