基于响应面分析的行波式能量吸收装置的参数优化

Optimization of Traveling-Wave Capability Absorption Parameters Based on Response Surface Analysis

摘要: 针对行波式能量吸收装置，其能量吸收效率受多种因素影响。为优化其性能，本研究对单个行波板进行了建模与仿真分析。运用响应面法，以无量纲波幅、无量纲波速及波数为变量，能量吸收效率为响应目标，进行了参数优化设计。结果显示，各因素对能量吸收效率的影响程度依次为：无量纲波幅、无量纲波速、波数。经过响应面优化，得到最优参数组合为无量纲波幅0.111、无量纲波速0.601、波数4.518，此时能量吸收效率可达84.27%。为进一步验证优化参数的有效性，进行了数值模拟对比，结果显示响应面预测值与数值计算值之间的误差仅为0.821%，从而验证了响应面模型在行波式能量吸收装置优化设计中的有效性。

Abstract: For the traveling wave energy absorption device, there are many factors that affect its energy absorption efficiency. To optimize its performance, this study conducted modeling and simulation analysis on a single traveling wave plate. Using the response surface method, the dimensionless wave amplitude, dimensionless wave speed, and wave number were studied as variables, with energy absorption efficiency as the response target for parameter optimization design. The results showed that the impact of each factor on energy absorption efficiency was ranked as follows: dimensionless wave amplitude, dimensionless wave speed, and wave number. After response surface optimization, the optimal parameter combination was obtained with a dimensionless wave amplitude of 0.111, a dimensionless wave speed of 0.601, and a wave number of 4.518, at which point the energy absorption efficiency could reach 84.27%. To further verify the validity of the optimized parameters, a numerical simulation comparison was conducted, which showed that the error between the predicted value of the response surface and the numerical calculation value was only 0.821%, thus verifying the effectiveness of the response surface model in the optimization design of the traveling wave energy absorption device.

1. 引言

风能作为一种无尽且可再生的清洁能源，对于降低人类对化石燃料的依赖具有显著意义 [1] 。然而，传统的流体机械存在灵活性不足和效率不高的问题，这凸显了探索新型能量吸收方法的紧迫性。通过研究和开发新型技术，我们可以更有效地利用风能，为可持续发展作出贡献。

自然界的动物，为适应多变且复杂的环境，经过漫长的自然选择过程，逐渐演化出高效且机动性强的各种运动模式 [2] 。这些模式不仅体现了动物的生存智慧，也给人类探索更高效的机械系统和运动方式提供了灵感。Lighthill [3] 将用于分析飞艇船体空气动力学的“长体理论”(Sender-Body Theory, SBT)引入鱼类波状摆动推进的分析中，提出了“细长体理论”(Elongated-Body Theory, EBT)。在同一时间，Wu [4] 将鱼体简化为二维平板，采用线性无粘流动以及振动可变形翼型的一般理论求解了任意非定常前进速度下游动的二维柔性板的一般问题，提出了二维波动板理论(2 Dimensional Wave Panel Theory, 2DWPT)。Lighthill [5] 将“细长体理论”进行了拓展，使之可以预测任意振幅、规则或不规则的鱼体与水之间的反作用力，提出了“大振幅细长体理论”(Large-amplitude Elongated-Body Theory, LEBT)。基于这个理论，童秉纲 [6] 和Cheng [7] 将扁平鱼类的波状游动简化为一块具有任意形状的无限薄柔性平板作波状摆动推进，从鱼类形态与其生存环境的适应性角度出发，建立了半解析–半数值的三维波动板理论(3 Dimensional Wave Panel Theory, 3DWPT)。Li [8] 等人采用网格变形技术，将鱼的行波简化为二维变形板的运动，以此进行水动力性能的数值模拟。在研究中，他们分析了波长、线性波幅、波频等主要参数对水动力性能的影响，并得出了平均推力系数和推进效率，以及随不同参数变化的压力曲线。

鱼类在水中游动，将能量传递给水，产生反卡门涡街，以此获得向前的推力。童秉纲 [9] 等人在研究三维波动板的推进性能时发现，只有当无量纲相速度大于1时，波动板才会获得向前的推力。近年来，Huang [10] 等经过研究发现，当行波运动的相速度大于来流速度时可以产生推力，此时，做行波运动的柔性体处于推进模式下；当行波运动的相速度小于来流速度时，做行波运动的柔性体可以从来流中吸收能量，此时行波运动处于能量吸收模式下。Ji等以NACA0012翼型为研究对象，通过数值模拟的方法，对这种新型能量吸收方式的影响因素，如Re [11] 和模型几何参数 [12] 进行了研究并与推进模式下的行波运动进行了对比。研究发现在Re一定时，行波运动的相速度与来流速度的比值是区分行波式推进模式和行波式能量吸收模式的关键参数，并且在这两种模式下模型几何参数的影响效果是完全相反的。林燕 [13] 等研究了NACA4、NACA63与NACA65系列不同最大厚度的翼型的获能特性，研究发现同一系列翼型的获能效率随最大厚度的增大而减小，且同一厚度，不同系列翼型随波长、振幅和无量纲波速的变化对获能效率影响不大。Zhang [14] [15] 则以柔性波动板为研究对象系统研究了波形几何参数，包括振幅(等振幅及不等振幅)和波长对其能量吸收特性的影响，发现在一定的参数范围内，这种新型能量吸收方式的能量吸收效率可以维持在较高的水平，值得进一步研究。Qi [16] 计算研究了钻石阵列中NACA0012翼型在纵向间距和相位差变化下的能量吸收效率。发现在钻石阵列中，随着间距增加，中游翼型效率先降后稳，下游翼型效率呈简谐波动。最佳排列时，阵列效率比单翼型提高35%。

然而，以往的研究主要侧重于单一因素对行波装置能量吸收效率的影响，实际上，行波装置的能量吸收效率是由多种因素共同决定的。因此，单因素研究存在一定的局限性，无法全面揭示行波装置能量吸收效率的复杂机制。

本文旨在探讨多种因素共同作用下，行波装置的能量吸收效率与相关参数之间的关系。为此，本文采用了响应面分析法(Response Surface Methods, RSM)，通过构建目标函数并研究其在多种参数下的敏感性，来深入分析各参数对能量吸收效率的影响。

通过这种方法，不仅能够更全面地理解行波装置能量吸收效率的影响因素，还能为优化行波装置的设计提供科学依据。这对于提高风能利用效率、推动清洁能源的发展具有重要意义。

2. 数值计算方法

2.1. 数学模型

本文采用的行波运动的数学模型如图所示，二维柔性波动板的运动方程为：

$y\left(x,t\right)=A\left(x\right)\cos \left[\frac{2\pi }{\lambda \left(x\right)}x-\frac{2\pi }{T}t+{\theta }_0\right]$ (1)

其中， $A\left(x\right)$ 为沿x方向变化的行波振幅， $\lambda \left(x\right)$ 为沿x方向变化的行波波长，T为行波运动的周期，t为流动时间， ${\theta }_0$ 为行波运动的初始相位。本文中柔性波动板的行波运动为如图1所示的等振幅、等波长行波运动。另外，通过对行波运动方程的线积分获取柔性行波板的模型长度：

$L={\displaystyle {\int }_0^{n\lambda }y\left(x,t\right)\text{d}s}$ (2)

式中L为柔性行波板长度，单位m，n为行波运动的波形数量，简称波数。

2.1.1. 模型的无量纲化和参数化定义

在对于模型的无量纲描述上，选取波长作为波形的特征长度，据此定义了无量纲波幅a并以来流速度对行波运动的相速度进行无量纲化定义了无量纲波速c，二者的定义式如下：

$a=\frac{A}{\lambda }$ (3)

$c=\frac{\lambda }{UT}$ (4)

Figure 1. Diagram of traveling wave motion

图1. 行波运动示意图

2.1.2. 侧向功率和能量吸收效率

柔性体在行波运动的过程中主要受到摩擦力和压差力 [17] 的作用，二者作用在模型表面微元上沿y方向的分量分别为：

$上表面\{\begin{array}{l}{f}_y^f=\mu \left[2\frac{\partial v}{\partial y}-\frac{\text{d}y}{\text{d}x}\left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right)\right]\hfill \\ f_y^p=-p\hfill \end{array}$ (5)

$下表面\{\begin{array}{l}{f}_y^f=\mu \left[-2\frac{\partial v}{\partial y}+\frac{dy}{dx}\left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right)\right]\hfill \\ f_y^p=p\hfill \end{array}$ (6)

其中，μ是流体的动力粘度，u和v分别为速度在x和y方向上的分量，p是压差力。

行波的能量吸收功率定义为：

$P={\displaystyle \oint \left(f_y^f+f_y^p\right)}{v}_y\text{d}s=P^f+P^p$ (7)

式中 $v_y$ 是模型表面微元体沿y方向的速度， $P^f$ 和 $P^p$ 分别为摩擦力和压差力的横向功率，用行波运动单位时间内所能吸收的最大来流动能E对能量吸收功率进行无量纲化，得到能量吸收效率(Energy Absorption Efficiency, EAE)：

$E=\frac{1}{2}\stackrel{\dot{}}{m}{U}^2=\frac{1}{2}\rho H{U}^3$ (8)

${\eta }_E\text{=}\frac{P}{E}$ (9)

式中，H为行波运动的最大扫掠高度，即波峰与波谷在y方向上的距离， $H=2A+d$ ，d为行波板的厚度， $\stackrel{\dot{}}{m}$ 为质量流量(kg/s)，ρ为来流密度。

2.2. 网格及边界条件

数值计算域及其尺寸如图2所示，其中区域B为半径3λ的圆域，区域A为半径40λ的圆域。流域入口设置为速度进口，出口设置为压力出口。采用了重叠网格技术配合动网格技术中的网格光顺保证良好的网格质量。区域A中的计算网格与前景网格均采用结构网格划分，区域B中的计算网格采用结构网格与非结构网格共同填充。网格细节如图3所示，其中灰色部分为背景网格，红色部分为前景网格，前景网格第一层网格高度为1.0e−5 m，严格保证计算网格y⁺小于1。本文采用二阶精度的计算格式，行波板表面采用无滑移边界条件，工质为空气，湍流模型为Transition SST模型。

Figure 2. Grid diagram

图2. 网格示意图

Figure 3. Overlapping mesh detail

图3. 重叠网格细节

2.3. 网格无关性验证和实验验证

为了保证计算网格的独立性以及时间步长的收敛性，分别在网格量为10w，20w和40w，时间步长为T/500，T/1000和T/2000下对网格进行了验证计算。如表1所示，最终选择网格量为20w，时间步长为T/1000进行后续计算。

Table 1. Independence verification

表1. 无关性验证

Nishio [18] [19] 等人对鱼类波状摆动推进的模型进行试验，如图4所示，实验以弦长为0.15 m，展长为0.2 m的NACA0018标准翼型为实验对象，研究了模型在不同巡航速度下头部与尾部相位滞后角φ对于模型“动态阻力系数Cdd”的影响。本文对相同条件下的二维NACA0018翼型进行数值模拟，模拟结果与实验值对比如图5所示。可以看出，实验值与模拟值变化趋势保持一致，可以认为本文所采用的计算方法和计算模型是可信的。

Figure 4. Schematic diagram of the hydrodynamic unit

图4. 水动力装置示意图

Figure 5. Comparison of numerical simulation pre-experiment results

图5. 数值模拟预实验结果对比

2.4. 实验方案设计

选择响应面分析法中的BBD (Box-Behnken Design)实验设计方法作为分析方法，实验设计中的3种参数的范围如表2所示。

Table 2. Horizontal distribution of factors

表2. 因素水平分布

3. 结果与讨论

响应面的17组实验结果已详细列于表3中。其中，最后五组实验为重复性实验，用以验证实验结果的稳定性和可靠性。为了确保这些重复性实验的结果具有统计显著性，并考虑到实际操作中可能存在的误差，本文人为地引入了一定的误差值。这些误差值的设计使得最后五组实验值数值平均分布在真实值的两侧，从而更全面地评估实验的可重复性。

Table 3. Response surface test and test results

表3. 响应面试验与试验结果

3.1. 显著性检验及模型验证

采用二次多项式方程对表中的实验数据进行拟合获得行波式能量吸收效率为响应值的回归方程：

$\begin{array}{c}{\eta }_E=-3.18429+17.32242A+7.37813B+0.44380C\\ +9.33186AB-0.67222AC+0.27519BC\\ -97.65156A^2-8.27101B^2-0.057659C^2\end{array}$

其中：A——无量纲波幅a；B——无量纲波速c；C——波形数量n。

Table 4. Significance analysis of traveling wave energy absorption efficiency regression model

表4. 行波式能量吸收效率回归模型显著性分析

对拟合得到的二次回归方程进行基于方差分析的显著性检验，结果如表4所示。其中，模型P值用于解释响应值的变化程度以及模型整体的显著性。当P值小于0.05时，表明模型达到极显著水平，即模型能够很好地描述实验数据的变化趋势。另一方面，失拟项P值用于评估模型预测结果与实验结果的不符程度。若失拟项不显著，则意味着模型预测值与实验结果之间不存在较大的误差，模型具有较高的预测精度。

通过表4中的数据分析，可以得出以下结论：首先，模型P值小于0.05，说明该二次回归模型在描述实验数据方面达到了极显著水平，模型具有较高的可信度。其次，失拟项P值不显著，表明模型预测值与实验结果之间的一致性较好，模型预测精度较高。

综上所述，该二次回归方程能够很好地描述实验数据的变化趋势，并且具有较高的预测精度。因此，本文可以基于该模型进行后续的优化设计和参数分析，以提高行波装置的能量吸收效率。

图6展示了预测值与试验值残差的正态概率分布图。从图中可以看出，残差近似分布在一条直线上，这表明残差近似满足正态分布。正态分布是回归分析中一个重要的假设，它保证了回归模型的稳定性和可靠性。因此，残差的正态分布说明了本文所采用的二次回归模型在预测行波装置能量吸收效率方面是合适的。

图7则展示了预测值与试验值的分布关系。从图中可以观察到，预测值主要分布在自左下向右上方的对角线周围。这意味着预测值与试验值之间存在着较高的一致性，二者之间的偏差较小。这种分布情况进一步验证了回归模型的预测效果，说明模型能够较为准确地反映实际情况下行波装置的能量吸收效率。

图8~图10为固定1个参数时剩余2个参数对 ${\eta }_E$ 的影响，观察各参数对于目标函数的影响趋势，对目

Figure 6. Normal probability distribution of residuals

图6. 残差的正态概率分布图

Figure 7. Distribution of predicted and experimental values

图7. 预测值与试验值分布图

标函数影响较大的参数进行分析。由图8和图9的等高线分布可见，等高线在因素A (无量纲波幅a)的坐标轴上呈现密集分布，这意味着无量纲波幅a对行波获能效率具有显著影响。相比之下，无量纲波速和波数的影响则相对较小，等高线在这些因素的坐标轴上分布较为稀疏。这一结果表明，在优化行波装置的设计时，应重点考虑无量纲波幅a的调整，以实现更高的能量吸收效率。

进一步观察图10，可以分析因素B (无量纲波速c)和因素C (波数n)之间的交互效应对获能效率的影响。从图中可以看出，因素B (无量纲波速c)的变化对获能效率的影响更为显著，等高线在该因素方向上呈现出更明显的变化。综上所述，行波获能效率的影响程度按照因素的影响大小排列为A > B > C，这表明改变无量纲波幅a可以更加有效地提高行波获能效率。

Figure 8. Efficiency model diagram of a and c

图8. a和c的效率模型图

Figure 9. Efficiency model diagram of a and n

图9. a和n的效率模型图

Figure 10. Efficiency model diagram of c and n

图10. c和n的效率模型图

3.2. 配合比优化与验证

基于回归模型的目标优化过程，本文旨在最大化行波装置的获能效率，通过调整无量纲波幅a、无量纲波速c和波数n这三个关键参数来实现。经过优化计算，本文得到了最优参数组合：无量纲波幅a = 0.111，无量纲波速c = 0.601，波数n = 4.518。根据回归模型的预测，这一最优组合下的获能效率可达84.27%。

为了验证模型预测的准确性，重新进行了数值模拟验证。数值模拟的结果显示，实际的效率为83.73%。将这一数值与模型预测的84.27%进行比较，发现二者之间的误差仅为0.821%，表明模型预测与实际模拟结果高度一致。

3.3. 流场分析

图11展示了行波板周围的总压云图，该图直观地反映了流体与行波板相互作用过程中的压力分布情况。从图中可以观察到，在一个周期的不同时刻，流体流经行波板后，在行波板的尾部区域出现了明显的低压区。这一现象表明，行波板在工作过程中确实吸收了一部分来流的能量，导致来流在经过行波板后能量有所减少。在行波板前面几个波处，总压减小较小，等到第二个波处，出现了总压的大幅降低。

图12展示了运动周期前半周期的静压云图，它揭示了行波板在运行过程中压力分布的动态变化。从图中可以看到，在行波板的前端，存在显著的压差，这意味着在此区域流体经历了较大的压力变化。这种压差的存在是行波板能够吸收能量的关键，因为它通过产生压力差来驱动流体的运动，进而实现能量的转换。

值得注意的是，尽管在行波板的前面几个波处，出现了少量与运动方向相反的逆压梯度，但随后的波形却全部顺着压力梯度横向运动。这种顺压力梯度的流动模式对于提升行波板的做功能力至关重要。它使得流体能够更加顺畅地通过行波板，减少了流动过程中的能量损失，并提高了整体的能量转换效率。

图13清晰地展示了行波板运动时尾迹的涡街结构。从图中可以看到，上游波峰和波谷处的分离泡在到达尾缘之前就已经脱落，这是由于行波板的周期性运动与来流之间的相互作用导致的。这些脱落的涡在来流中向下游运动的速度介于行波板的相速度和来流速度之间。

Figure 11. Total pressure cloud picture

图11. 总压云图

Figure 12. Static pressure cloud picture

图12. 静压云图

Figure 13. Vortex cloud picture

图13. 涡量云图

值得注意的是，这些上游波峰和波谷处脱落的涡在向下游运动的过程中，会受到下游波峰或波谷的阻拦。当它们接近下游的波峰或波谷时，会与这些位置附近的其他脱落涡融合。这种相互作用和融合导致涡的尺寸逐渐增大，最终在行波板的尾缘处脱落出尺寸更大的旋涡。这些大尺寸的旋涡在行波板的尾迹中形成了卡门涡街。

4. 结论

通过对行波板各参数设计进行响应面分析，得出以下结论：

(1) 通过对行波板的3种参数敏感性分析得出了不同参数对目标函数的影响，并按照由强到弱的顺序进行排序，a > c > n；同时发现，无量纲波幅a对行波板获能效率的影响最突出，其次时无量纲波速c。

(2) 对行波板获能进行仿真分析，明确了不同参数下效率的变化规律，为回归方程建立奠定了基础。

(3) 在研究参数范围内，得出了三种参数的最佳匹配(a = 0.111, c = 0.601, n = 4.518)，使得获能效率达到最高，达到84.27%。

基金项目

国家自然科学基金重点项目(52036005)；国家自然科学基金原创探索计划项目(52350139)。