海上浮动式风力机设计的前提是要对其进行准确的气动-水动-伺服-弹性时域全耦合分析, 即通过数值积分方法求解其全耦合运动微分方程从而获得风机系统响应结果的时间历程数据。在这组微分方程中存在一个关于水动力辐射载荷求解的时域卷积项, 作者在前期研究的基础上[1-3]发现该卷积项的计算效率并不高, 为了保证计算精度, 每一时刻的卷积计算均需要获取先前足够长时间段内的风机浮动平台的运动信息[4-5]。卷积项的计算规模与风机浮动平台运动自由度数的平方成正比关系, 数值积分求解将会非常耗时且需要占用大量的计算存储空间。此外, 该卷积项使系统运动方程难以转换为状态空间的表述形式, 从而增加浮动风机控制系统分析和设计的困难[6]。
为了克服上述时域卷积计算方法的缺点, 利用状态空间模型进行水动力辐射载荷的计算则可作为一种近似替代方法, 只需根据当前时刻平台的速度信息即可演化未来时刻的水动力辐射载荷。本质上说, 水动力辐射载荷的记忆效应已包含于系统状态矢量中, 因此可避免卷积计算需要大量过去时刻的速度信息及对应的脉冲响应函数值。目前, 利用状态空间模型等效计算水动力辐射载荷的方法主要应用于船舶[7-10]和海洋波能转换装置[11-13]等领域研究中, 而在海上风机领域的研究报道仅见到一篇[14]。
水动力辐射载荷状态空间模型法的关键是状态空间模型的系统识别问题, Taghipour等利用时域最小二乘法, 数据实现理论和频域回归分析等3种系统辨识方法识别现代大型货船的水动力辐射载荷的状态空间模型, 仿真所得水动力辐射载荷均具有较高的精度[8-9]。Rogne等利用频域最小二乘法成功识别了海洋波能转换装置的水动力辐射传递函数模型, 进而将其转换为状态空间模型, 较好的模拟了该设备的水动力辐射载荷[13]。Duarte等利用最小二乘法识别了TLP和OC4浮动风机的水动力辐射状态空间模型, 发现模型精度随状态阶数的提高而增加, TLP的计算精度高于OC4[14]。然而Duarte并未将状态空间模型与风机全耦合运动方程相耦合, 风机运动的累积误差会影响水动力辐射载荷的计算精度, 因此, Duarte识别的状态空间模型的精度有待重新评估。
本文提出海上浮动式风力机水动力辐射载荷的状态空间实现法。利用一组参数化状态空间模型替代传统时域卷积法来实现海上浮动式风力机水动力辐射载荷的等效计算。分别从频域和时域2个方面选取4种系统识别方法进行状态空间模型的参数识别; 将状态空间模型与风机动力学分析程序FAST及其子程序HydroDyn进行耦合。以时域卷积法仿真结果为标准, 从仿真精度和时间两方面综合评估状态空间法水动力辐射载荷的计算效率。本文工作为海上浮动式风力机动力学仿真奠定了基础。
1 水动力辐射载荷时域卷积模型 1.1 基于脉冲响应函数的时域卷积法海上浮动式风机水动力辐射载荷是由于平台运动引起的出流波反作用于平台结构而产生的脉动压力作用, 包括水动力附加质量和附加阻尼两部分效应[15], 其表达式为
|
(1) |
式中, 等号右端第一项为由平台运动加速度引起的水动力附加质量力, 可知其与平台运动加速度同相位; 第二项表示水动力附加阻尼力即水动力辐射载荷的记忆效应。目前几乎所有的仿真程序均对积分时间进行了截断, 即
|
(2) |
式中, tmemory表示截断时间, 即只存储tmemory时间段内的平台速度信息与有限长度的脉冲响应函数, 目前一般取为60 s。卷积法的精度与取决于截断时间长度, 与计算时间、数据存储空间和脉冲响应函数求解规模相矛盾。
1.2 脉冲响应函数的时频域关系(2) 式中的卷积项在频域中被表示为
|
(3) |
式中, A(ω)、B(ω)分别表示频域水动力附加质量矩阵和附加阻尼矩阵。
利用三角函数加权形式的傅里叶变换方法可推得频域水动力系数矩阵和时域脉冲响应函数矩阵的关系为
|
(4) |
(4) 式两式均可反推时域脉冲响应函数, 但为了避免对脉冲响应函数在初始时刻进行误差修正, 学术界一般通过(4)式第二式来反推时域脉冲响应函数, 即
|
(5) |
为了更加高效的计算水动力辐射载荷的记忆效应, 本文提出求解波浪辐射载荷的状态空间实现法, 假定系统具有因果性和时不变性, 则(2)式中的卷积项可由状态空间模型近似替代:
|
(6) |
式中, Fdamping表示水动力辐射载荷中的附加阻尼力; x为状态向量; A、B、C分别表示系统矩阵、输入矩阵和输出矩阵, 此3个矩阵为待识别的参数矩阵。
关于该状态空间模型的识别问题, 本文分别从频域和时域2个方面进行系统参数识别, 识别过程如图 1所示。
|
| 图 1 状态空间模型参数识别过程 |
提出如(7)式所示的评估指标, 该评估指标反映了参数识别模型的变异性。该指标越接近1, 表示参数模型识别的精度越高。
|
(7) |
式中, a, b=1, 2, …, 6;Kab为参考脉冲响应函数平均值; 对于频域识别方法, l=ω, 即按照各频率点值求和; 对于时域识别方法, l=t, 即按照各时间点值求和。
2.3 频域参数模型识别频域脉冲响应函数矩阵中的每一项均可被拟合为一个适当阶数的传递函数模型:
|
(8) |
式中, u=[pm, pm-1, …, p0, qn-1, qn-2 …, q0], 为传递函数分子P(s, u)、分母Q(s, u)各项的系数组成的待识别参数矢量; s=jω。
根据水动力辐射势流理论, 脉冲响应函数应当满足一些时频域上的约束条件[16], 为了节省篇幅, 此处仅列出结论, 如表 1所示。根据这些约束条件可确定参数矢量u中的部分系数。
通过(3)式得到参考频域脉冲响应函数矩阵, 由此可对(8)式传递函数模型进行拟合, 此问题可转化为(9)式所示的最小二乘法优化问题, 采用Alves等提出的考虑频率范围权重因子的频域最小二乘法[17](Freq-LS)和Perez等提出的频域迭代最小二乘法[18](Freq-ILS)进行求解, 方法详细内容见相关文献, 此处不再赘述。
|
(9) |
时域脉冲相应函数识别方法同样可以得到状态空间模型, 时域识别方法需要利用参考时域脉冲响应函数, 可以通过以下2种方法获得:
1) 逆傅里叶变换(IFFT):对参考频域脉冲响应函数作Kab(ω)傅里叶逆变换即可得到参考时域脉冲响应函数Kab(t)。但这种方法受到奈奎斯特抽样频率的限制, 会带来附加误差, 且经过逆傅里叶变换得到的Kab(t)是在[0, t]内等间隔均匀取值的, 结合时域脉冲响应函数在初始时间段内变换剧烈的特点, 当模拟时间t较大时, 该法有可能导致不能精确描述初始时间段内的脉冲响应函数值。
2) 余弦变换:根据时域脉冲响应函数和附加阻尼矩阵之间的余弦转换积分关系式(5), 本文利用梯形积分法对其求解, 表达式为
|
(10) |
式中, kmax为频率离散点数, Δω为频率间隔; 该法比IFFT计算更方便, 且具有更高的计算精度, 因此本文选用此法。选取频率上限为ωu=5 rad/s, kmax=256, 则Δω≈0.02 rad/s; 时间上限tu=100 s, 时间间隔Δt=0.1 s。
采用基于Prony算法[19]的时域最小二乘法(Time-LS)和基于Hankel矩阵奇异值分解法的时域数据实现理论[20](Time-RT), 同样方法详细内容见相关文献, 此处不再赘述。
3 模型识别结果对比 3.1 模型及工况选取本文选取Barge式、Spar式和TLP式3种典型的浮动风机为研究对象, 其浮动平台的详细参数见文献[3]。
为了比较本文状态空间模型和传统卷积模型在真实海洋环境工况下的仿真结果, 选取和文献[21]一致的不同有效波高和谱峰周期组合的8种参考海洋工况, 如表 2所示, 其中工况1为微幅波况, 工况8为极端波况。各工况下的仿真时间取为1 h, 步长为0.012 5 s。
时域脉冲响应函数矩阵和频域脉冲响应函数矩阵均有10个非零元素, 即主对角线的6个元素K11、K22、…、K66, 非主对角线的4个元素K15、K51、K24和K42。根据脉冲响应函数矩阵的对称性可知K15=K51、K24=K42, 只需识别K15和K24即可。因此, 每个脉冲响应函数矩阵有8个待识别元素。分别预设评估指标精度为0.97和0.99, 利用上述4种方法进行识别。为了节省篇幅, 仅展示Barge式风机平台利用时频域方法识别的脉冲响应函数矩阵各元素的识别精度与所需状态阶数, 如图 2所示。
|
| 图 2 Barge式平台时频域脉冲响应函数精度及所需状态阶数 |
从图 2a)可知, 4种方法在预设精度为0.97的条件下识别所得状态空间模型均具有良好的稳定性, 并未出现响应发散的情况。Freq-ILS法和Time-RT法所需的状态阶数较少, 而Freq-ILS法识别的精度略低于Freq-LS法和Time-RT法, 某些元素的识别精度甚至低于预设精度0.97, 这是因为Freq-ILS法的直接识别对象是附加质量和附加阻尼矩阵, 而脉冲响应函数是间接所得。Time-RT法兼有所需状态阶数少和时频域响应精度均较好的优点。Freq-LS法识别精度虽然也较好, 但所需状态阶数较Freq-ILS法和Time-RT法多。Time-LS法所需状态阶数最多, 但识别精度却不理想, 其中频域响应函数K33的识别精度甚至低于0.75。
从图 2b)可知, 当预设识别精度提高至0.99时, Freq-LS法和Time-LS法识别的状态空间模型具有不稳定性, 导致响应结果发散。结合其他2种风机识别结果来看, 当预设识别精度提高后, 各方法(除发散情况外)所得识别结果的精度也相应得到了提高, 但同时也会增加所需的状态阶数。
3.3 不同预设精度下模型所需状态阶数对比状态空间模型阶数越高仿真所需时间越长, 为了保证状态空间模型相比传统时域卷积模型具有仿真时间的优势, 在精度满足的条件下, 应尽量减少模型的状态阶数。从3.2节可知, 在相同的预设精度条件下, 4种方法识别所需的状态阶数有较大不同。
图 3为不同预设精度条件下4种方法识别的3类平台状态空间模型所需的状态阶数, 可知4种方法识别的状态空间模型阶数随预设识别精度的增加而增加。当预设识别精度小于0.97时, 各识别方法所需状态阶数的增加幅度较为平缓; 当预设识别精度大于0.97后, 状态阶数的增加幅度较大, 该规律对于Freq-ILS法和Time-RT法更为明显。因此, 利用较低阶的状态空间模型来达到较好识别精度的思路是可行的。对比4种识别方法, 可以更加清楚地看到Freq-ILS法和Time-RT法识别所需状态阶数最少, 其次为Freq-LS法, 而Time-LS法所需状态阶数最多。对比3类平台, 可知各预设识别精度下, Spar平台和TLP平台所需状态阶数相当, 而Barge平台则需要更高的状态阶数, 说明Barge平台的立方体外形对其水动力辐射特性的影响比Spar和TLP平台的圆柱状外形复杂, 因此, 可进一步推论平台外形越复杂, 其水动力辐射特性越复杂, 识别所需的状态阶数越高。
|
| 图 3 不同预设精度下4种方法识别的状态空间模型阶数对比 |
本节将利用状态空间模型方法计算水动力辐射载荷, 并和传统时域卷积模型法的计算结果进行对比。为了对状态空间模型法计算水动力辐射效应进行全面评估, 考虑状态模型法的非耦合和耦合2种计算方式, 如图 4所示。
|
| 图 4 状态空间模型法的非耦合和耦合计算方式 |
图中下半部分表示传统时域卷积模型法的计算过程, 上半部分为本文状态空间法的计算过程。
4.1 状态空间模型法求解的水动力辐射载荷的精度分析以海况6为例, 各识别方法在不同预设识别精度下的状态空间模型计算的3种典型平台的纵向、横向水动力辐射载荷与卷积模型方法所得结果的相对误差如图 5、图 6所示。4种识别方法中, Time-LS法状态空间模型的水动力辐射载荷计算精度最差, Time-LS法和Freq-LS法状态空间模型对于Barge平台横向水动力载荷的计算结果发散, 说明这2种方法识别的状态空间模型的稳定性低于另外两种方法。Spar平台和TLP平台的计算精度均很高, 其中Spar平台的纵向和横向水动力辐射载荷计算精度可达到0.999 996和0.998, TLP平台则为0.999 97和0.999 92, 且两者的计算精度与预设识别精度关系不大, 因此, 对于这2种平台, 建议采用较低识别精度下的状态空间模型进行水动力辐射载荷计算。而Barge平台的计算精度则比Spar和TLP要差, 当预设识别精度小于0.95时, 其水动力辐射载荷的计算精度与预设识别精度的关系不大; 当预设识别精度大于0.95时, 其水动力辐射载荷的计算精度随预设识别精度的增大而提高, 当预设识别精度从0.95增大到0.995时, 其纵向水动力辐射载荷计算精度从0.991提高至0.999 4, 而横向水动力辐射载荷精度则从0.854提高至0.978。
|
| 图 5 不同方法计算的纵向水动力辐射载荷精度比较 |
|
| 图 6 不同方法计算的横向水动力辐射载荷精度比较 |
图 7给出了不同海况下利用状态空间模型法计算的纵、横向水动力辐射载荷与卷积模型法所得结果的相对误差, 其中Barge平台采用高预设精度0.995, 而Spar和TLP平台则采用低预设精度0.7。图中所示3种方法识别的状态空间模型在各海况下计算的水动力辐射载荷的精度均很高, 其中横向水动力辐射载荷的计算精度略低于纵向, 这是因为波浪方向为纵向, 导致横向水动力辐射载荷比纵向水动力辐射载荷小约2个数量级, 致使横向水动力辐射载荷的相对误差大于纵向水动力辐射载荷。
|
| 图 7 不同海况下状态空间模型水动力辐射载荷相对于卷积模型的精度 |
定义仿真计算效率指标Tratio=Tsim/Tcpu, 其中, Tsim为水动力仿真模拟时长, Tcpu为计算机仿真计算用时。该指标说明在水动力仿真模拟时长一定的条件下, 计算机仿真计算用时越短, 仿真计算的效率越高。
状态空间模型在解耦和耦合2种计算方式下的仿真时间效率分别如图 8、图 9所示。由图可知, 3种平台的仿真时间效率由高到低依次为Barge、TLP和Spar。非耦合计算方式下时域卷积法的仿真时间效率与卷积截断时间有关, 随截断时间的增加呈反比例下降趋势, 而耦合计算方式下则呈线性下降趋势。各辨识方法所得状态空间模型的仿真时间效率大致相同。除Barge平台在非耦合计算方式下的仿真时间效率随状态阶数的增加而略有降低外, 其余所有情况的仿真时间效率与状态阶数的关系并不大。状态空间模型法的仿真时间效率远高于传统时域卷积法, 以时域卷积模型常用的60 s截断时间下的仿真时间效率为标准, Barge、TLP和Spar 3种平台在非耦合计算方式下状态空间模型法的平均仿真时间效率为时域卷积模型法的11.1、3.8和3.3倍, 对应仿真用时可减少90%、74%和71%;耦合计算方式下仿真状态空间模型法的平均仿真时间效率为时域卷积模型法的1.8、1.5和1.4倍, 对应仿真用时可减少45%、33%和31%。
|
| 图 8 状态空间模型法与时域卷积模型法非耦合计算仿真时间效率对比 |
|
| 图 9 状态空间模型法与时域卷积模型法耦合计算仿真时间效率对比 |
本文提出海上浮动式风力机水动力辐射载荷仿真的状态空间实现法, 利用4种系统识别方法识别了状态空间模型的参数, 比较了各方法识别的状态空间模型在浮动风机水动力辐射载荷仿真方面的优缺点, 主要结论如下:
1) 频域最小二乘法、迭代最小二乘法和时域实现理论均可用较少的状态阶数较好地识别系统时频域脉冲响应函数, 而时域最小二乘法需要最高的状态阶数, 且识别结果不甚理想。
2) 频域迭代最小二乘法状态空间模型和时域实现理论状态空间模型所得3种典型风机浮动平台水动力辐射载荷仿真结果具有很高的计算精度。对于Spar和TLP平台, 计算精度与预设识别精度关系不大, 建议采用低预设识别精度的状态空间模型; 而对于Barge平台, 当预设识别精度小于0.95时, 计算精度与预设识别精度的关系不大, 当预设识别精度大于0.95时, 计算精度随预设识别精度的增大而提高较明显, 建议采用高预设识别精度的状态空间模型。
3) 浮动风机水动力辐射载荷的状态空间实现法比传统时域卷积法具有明显的时间仿真优势, 其中非耦合计算方式下仿真用时可减少70%以上, 耦合计算方式下可减少仿真用时30%以上。
| [1] |
贺尔铭, 胡亚琪, 张扬. 基于TMD的海上浮动风力机结构振动控制研究[J]. 西北工业大学学报, 2014, 32(1): 55-61.
He Erming, Hu Yaqi, Zhang Yang. Structural Vibration Control of Offshore Floating Wind Turbine Based on TMD[J]. Journal of Northwestern Polytechnical University, 2014, 32(1): 55-61. (in Chinese) |
| [2] | He E M, Hu Y Q, Zhang Y. Optimization Design of Tuned Mass Damper For Vibration Suppression of a Barge-type Offshore Floating Wind Turbine[J]. Journal of Engineering for the Maritime Environment, 2017, 231(1): 302-315. |
| [3] |
贺尔铭, 张扬, 胡亚琪. 3种典型海上浮动式风力机动力学特性比较分析[J]. 太阳能学报, 2015, 36(12): 2874-2881.
He Erming, Zhang Yang, Hu Yaqi. Comparison And Analysis of Dynamic Characteristics of Three Typical Floating Wind Turbine[J]. Acta Energiae Solaris Sinica, 2015, 36(12): 2874-2881. DOI:10.3969/j.issn.0254-0096.2015.12.006 (in Chinese) |
| [4] | Jonkman J M, Sclavounos P D. Development of Fully Coupled Aeroelastic And Hydrodynamic Models for Offshore Wind Turbines[R]. National Renewable Energy Laboratory, Golden, CO, 2006 |
| [5] | Jonkman J M. Dynamics of Offshore Floating Wind Turbines——Model Development And Verification[J]. Wind Energy, 2009, 12(5): 459-492. DOI:10.1002/we.v12:5 |
| [6] | Sheng W, Alcom R and Lewis A. A New Method for Radiation Forces for Floating Platforms in Waves[J]. Ocean Engineering, 2015, 105: 43-53. DOI:10.1016/j.oceaneng.2015.06.023 |
| [7] | Unneland K. Identification And Order Reduction of Radiation Force Models of Marine Structures[D]. Trondheim, Norwegian University of Science and Technology, 2007 |
| [8] | Taghipour R, Perez T, Moan T. Hybrid Frequency-Time Domain Models for Dynamic Response Analysis of Marine Structures[J]. Ocean Engineering, 2008, 35(7): 685-705. DOI:10.1016/j.oceaneng.2007.11.002 |
| [9] | Taghipour R, Perez T, Moan T. Time-Domain Hydroelastic Analysis of a Flexible Marine Structure Using State-Space Models[J]. Journal of Offshore Mechanics and Arctic Engineering, 2009, 13(1): 011603. |
| [10] | Hatecke H. The Impulse Response Fitting And Ship Motions[J]. Ship Technology Research, 2015, 62(2): 97-106. DOI:10.1179/2056711115Y.0000000001 |
| [11] | Alves M. Numerical Simulation of The Dynamics of Point Absorber Wave Energy Converters Using Frequency And Time Domain Approaches[D]. Lisbon, Technical University of Lisbon, 2012 |
| [12] | Armesto J A, Guanche R, Iturrioz A, et al. Identification of State-Space Coefficients for Oscillating Water Columns Using Temporal Series[J]. Ocean Engineering, 2014, 79: 43-49. DOI:10.1016/j.oceaneng.2014.01.013 |
| [13] | Rogne Ø Y, Moan T, Ersdal S. Identification of Passive State-Space Models of Strongly Frequency Dependent Wave Radiation Forces[J]. Ocean Engineering, 2014, 92: 114-128. DOI:10.1016/j.oceaneng.2014.09.032 |
| [14] | Duarte T, Alves M, Jonkman J, et al. State-Space Realization of The Wave-Radiation Force within FAST[C]//Proceedings of the ASME 201332nd International Conference on Ocean, Offshore and Arctic Engineering, Nantes, France, 2013 |
| [15] | Jonkman J M. Dynamics Modeling And Loads Analysis of An Offshore Floating Wind Turbine[D]. Boulder, University of Colorado, 2007 |
| [16] | Perez T, Fossen T. Time-vs. Frequency-Domain Identification of Parametric Radiation Force Models for Marine Structures at Zero Speed[J]. Modeling, Identification and Control, 2008, 29(1): 1-19. DOI:10.4173/mic.2008.1.1 |
| [17] | Alves M, Vicente M, Sarmento A, et al. Implementation and Verification of a Time Domain Model to Simulate the Dynamics of OWCs[C]//9th European Wave and Tidal Energy Conference, Southampton, UK, 2011 |
| [18] | Perez T, Fossen T. Practical Aspects of Frequency-Domain Identification of Dynamic Models of Marine Structures from Hydrodynamic Data[J]. Ocean Engineering, 2011, 38(2): 426-435. |
| [19] | Zygarlicki J, Mroczka J. Short Time Algorithm of Power Waveforms Fundamental Harmonic Estimation with Use of Prony's Methods[J]. Metrology and Measurement Systems, 2011, 18(3): 371-378. |
| [20] | Alvin K F, Robertson A N, Reich G W, et al. Structural System Identification:from Reality to Models[J]. Computers & Structures, 2003, 81(12): 1149-1176. |
| [21] | Jonkman J. Definition of the Floating System for Phase IV of OC3[R]. National Renewable Energy Laboratory(NREL), Golden, CO, 2010 |
