一种风帆助航船舶航线节能评估方法

本发明涉及船舶航行节能技术领域，尤其是一种风帆助航船舶航线节能评估方法。

背景技术：

根据国际原子能机构发布的相关数据，船舶运输业作为大能源消耗产业之一，其排放的sox、nox和co2分别占全球排放总量的10％、25％和2.7％，所耗能源约占世界总能耗的2.5％。中国提出在2030年前二氧化碳排放达到峰值，力争2060年前实现碳中和。如何节能减排，提高船舶营运的经济性已成为航运业越来越关注的问题。

自2005年imo公布的《防止船舶造成大气污染规则》生效后，海运己被纳入对节能环保具有积极意义的产业。船舶节能与替代能源技术，尤其是船舶可再生能源技术，是具有前瞻性的船舶工业核心技术，这事关国家船舶工业可持续发展能力。交通运输部颁布的《综合运输服务“十三五”发展规划》指出，风能作为地球上取之不尽用之不竭的清洁能源，具有非常大的发展潜力。

对于现代风帆船，由于计算机控制技术的高度发展，风帆在不同风况下智能控制已成为可能，同时由于卫星等气象监控技术的发展，根据船舶航线上天气预测调整风帆应用也变得可行。总之现代风帆船，随着科学技术的发展正日趋成熟，将成为未来新型环保船舶的重要发展方向。

研究全球范围内海洋风力资源的分布特征，是评估其对风帆助航船舶的可用性以及推进效率的基础。但是，全球海上风能资源的分布并不均匀:从大气环流的角度看，全球海上风能资源的分布与全球气压带和风带的分布密切相关，一般越往赤道风越小，越往两极风越大，同时有区域气候和特殊地形形成的区域性风资源。并且，在同一条航线上，不同时间，不同航向下风力资源分布也不同，因此建立一种针对风帆助航船舶的航线节能评估方法对实现风帆助航船舶的推广及其技术的完善是很有必要的。

技术实现要素：

本发明针对以上问题提出了一种风帆助航船舶航线节能评估方法。

本发明采用的技术手段如下：

一种风帆助航船舶航线节能评估方法，包括以下步骤：

步骤1、获取目标航线和全球风场数据库，从所述全球风场数据库中提取所述目标航线上的风场数据；

步骤2、获得船舶航速和船舶航向，根据所述船舶航速和船舶航向计算目标航线上的风场数据与船舶的相对风向和相对风速；

步骤3、建立包括相对风向和相对风速的风场分布矩阵；

步骤4、统计所述目标航线上的风场数据中落入每个风场分布矩阵的矩阵元素范围内风场的个数；

步骤5、根据所述落入每个矩阵元素范围内风场的个数与目标航线上的风场数据中风场的总个数计算每个矩阵元素范围内风场概率获得风场概率分布矩阵；

步骤6、根据所述风场概率分布矩阵计算风帆载荷概率矩阵；

步骤7、根据所述风帆载荷概率矩阵计算风帆总推力，并通过所述风帆总推力计算风帆的节能效率。

进一步地，所述风帆载荷概率矩阵中每个矩阵元素通过公式(1)进行计算：

fij＝fij×fsk(1)

其中：fij为风帆载荷概率矩阵中每个矩阵元素值，fij为风场概率分布矩阵中每个矩阵元素值，fsk为fij对应的风速风向下风帆能够提供的推力。

进一步地，所述风帆总推力fm通过公式(2)进行计算：

其中，i,j为风载荷概率矩阵fij索引值，n为风载荷概率矩阵fij的数据个数；

所述风帆的节能效率η通过公式(3)进行计算：

其中：sv为航速，p为船舶的有效功率。

与现有技术比较，本发明公开的风帆助航船舶航线节能评估方法具有以下有益效果：将风帆空气动力学特性、风帆助航船舶及其航线相结合来有效评估节能效果,本方法中通过获取全球风场数据库，并根据全球风场数据库建立目标航线的风场概率分布矩阵，进而可以根据目标航线的风场概率分布矩阵计算风帆载荷概率矩阵，从而可以计算风帆总推力获得风帆节能效率。

附图说明

图1为本发明公开的风帆助航船舶航线节能评估方法的流程图；

图2为航向确定的示意图；

图3为曲面差值的原理图；

图4为航向定义的示意图；

图5为船舶运动坐标系的示意图；

图6为真实风向的计算示意图；

图7为相对风向的计算示意图；

图8为相对风速的计算示意图；

图9为风场概率矩阵计算的流程图；

图10为风帆布置的示意图；

图11为机翼原理示意图；

图12为风帆坐标系及受力分析示意图；

图13为风帆受力图。

具体实施方式

如图1所示为本发明公开的风帆助航船舶航线节能评估方法，包括以下步骤：

步骤1、获取目标航线和全球风场数据库，从所述全球风场数据库中提取所述目标航线上的风场数据；

具体地，目标航线的获取过程如下：确定已知的航线经纬度信息，各航行节点记为pi，pi(longitude，latitude)，其中longitude表示该航行节点经度，latitude表示为该航行节点纬度，航线航程记为l(km)，航速记为sv(kn)，以及到达每个航行节点的时间t。根据相邻两个航行节点可确定航向sd(0-360°，以正北为0°，顺时针为正)。航向确定方式为如图2中，p1、p2为相邻两个航行节点，将p2、p1经纬度相减，取反正切函数可得航向sd。

本申请中全球风场数据库是从海洋气象数据中获得，海洋气象数据为来自欧洲中期天气预报中心(ecmwf)第五代再分析资料era5。era5提供了大量大气、陆地和海洋气候变量的每小时估计数。这些数据覆盖地球30公里的网格，并使用137个从地面到80公里高度的水平解析大气。era5包括关于所有变量在空间和时间分辨率降低时的不确定性的信息。era5在3个月内实时发布数据集，初步更新可在实时的5天内提供给用户，其组合大量的历史观察使用先进的建模和数据同化系统进行全球估计。该数据集的时间范围为1979年至今，包含着大气、海洋波浪以及陆地表面的气象数据，其时间分辨率为6h(每天00：00，06:00，12:00，18:00utc)，平面分辨率最高可为0.125°×0.125°。

本申请中的海洋气象数据的时间分辨率为6h(每天00：00，06:00，12:00,18:00utc)，平面分辨率为0.25°×0.25°，调取的海洋气象数据为：10m高度u、v方向的风场数据、平均压力以及海洋表面温度。依据历史航行数据，以时间为单位，提取相应航行位置点的风场数据。

从所述全球风场数据库中提取所述目标航线上的风场数据具体为：以航线所处经纬度和航行时间为依据，采用时空插值的方法，从全球风场数据库中提取该航线上所有风速风向数据。

具体地，由于地球是球体，船舶航行时实际上是在球面上航行，本申请中获得的从ecmwf上下载包含各个航线的风场数据(nc格式文件)为球面数据，因此插值时需要用到曲面插值。

确定曲面插值方法是基于信息点之间的相似程度或者整个曲面的光滑性来创建一个拟合曲面，其原理如图3所示：图中gridpoints(网格点)是指已知的节点的值，interpolatedpoints(插值点)是指曲线上的插值，f(x，y)是指曲面上的插值。

船舶所处风场信息提取具体方法如下：

(1)从ecmwf上分别下载包含各个航线的风场数据：海平面10m处的u方向上的风速和v方向的风速，网格精度为0.25×0.25，每天的采样时间点为00:00:00、06:00:00、12:00:00、18:00:00，格式为nc格式。

(2)读取nc文件中的time(时间)，longitude(经度)，latitude(纬度)数据。

(3)读取nc文件中所需要提取的风场信息数据(u10，v10)。

(4)按utc时间输入航线坐标时间，并且将其转换为nc文件格式。

(5)根据航线坐标经纬度及时间，使用时空插值的方法提取风场数据。

(6)首先判断坐标时间是否为6小时整数倍，如果在整数时间点上，则进行二维曲面插值，提取当前坐标点的风场数据。

(7)如不在整数时间节点上，则读取前后两个时间节点上的数据，然后进行插值。

通过以上步骤提取数据结果为一个一维列向量，其每个元素为航线对应坐标的天气信息。利用上述方法提取所需的所有船舶所处风场的风场信息。

步骤2、获得船舶航速和船舶航向，根据所述船舶航速和船舶航向计算目标航线上的风场数据与船舶的相对风向和相对风速；

具体地，船舶的航速可以根据实际需要进行确定。船舶航向的获得过程如下：船舶航向的定义，本申请中所指的航向均为真航向，其定义为船舶航行时，真北线与船首线之间的夹角。真航向以真北线为基准，顺时针量至船首线，如图4所示，船首朝向正北为0°，顺时针为正向，范围为0°-360°。

船舶运动坐标系是固结在船体上的，随船一起运动，如图5所示。

图5中，o-x1y1为大地坐标系，o-x1正方向为正北，o-y1正方向为正东；g-xyz构成一右手法则的固结在船体上的坐标系；g是动坐标系的原点，通常可选取船舶重心或者船中剖面处；gx轴为船纵轴，其方向指向船首为正；gy轴与纵剖面垂直，以指向右舷为正；gz轴垂直于水线面，以指向龙骨为正。

船舶实际航速vs与船艏中线间的夹角称为漂角γ,规定当f转向gx正方向为顺时计时为正，而船舶偏转方向与船舶航速和装载等因素有关；vs在x和y轴方向的速度分量分别记为u和v；船体绕gz轴的角速度为r；真风风速用vt表示，真风来风角度记为βt；船上直接测量、作用于风帆及船体上风能均为相对风，相对风速记为va，相对风来风角度为β,规定左舷来风为正，右舷为负；相对风向与风帆弦向夹角即为风帆攻角，记为αs，规定由相对风至风帆弦向顺时针为正；风帆转角为φ，规定当gx至风帆弦向为顺时针时为正。

相对风向和相对风速的计算方法如下：

1、计算真实风向

由于本专利获取的nc文件中风速的存储方式是按uv坐标存储的，uv坐标是指所有的图象文件都是二维的一个平面。水平方向是u，垂直方向是v，通过这个平面的，二维的uv坐标系，因此我们首先需要利用uv方向风的分量计算真实风向，如图6所示。

真实风向计算方法为，首先根据提取的航向节点上的uv风向，按不同正负号分为四个象限来处理，一象限，uv均为正号；二象限u为正号，v为负号；三象限，uv均为负号；四象限，u为负号，v为正号。然后在每个象限上分别利用反正切函数，可算得该节点的真实风向。

2、计算相对风向

本专利中的相对风向定义为真实风向与船舶航向的矢量叠加方向，如图7所示。

相对风向计算方法为用绝对风向(真实风向)与航向的差取绝对值，即

β＝|θ-sd|(4)

式中θ为绝对风向，sd为航向。

3、计算相对风速

本专利相对风速计算方式为首先在航线节点处将船舶航速sv按uv坐标方向进行分解，得到船舶航速在uv方向上的分量，然后分别在uv方向上与提取的nc文件中的真实风速在uv方向上的分量进行矢量叠加，叠加后得到相对风速在uv方向上的分量，最后利用三角函数可由相对风速在uv方向上的分量求得相对风速，如图8所示。

图中sdv是航速在v方向上的分量，sdu是航速在u方向上的分量；vtv是真实风速在v方向上的分量，vtu是真实风速在u方向上的分量；由sdv和vtv可计算出相对风速在v方向上的分量vav，由sdu和vtu可计算出相对风速在u方向上的分量vau；由vav和vau可计算出该航线节点处的相对风速va。

步骤3、建立包括相对风向和相对风速的风场分布矩阵；

步骤4、统计所述目标航线上的风场数据中落入每个风场分布矩阵的矩阵元素范围内风场的个数；

步骤5、根据所述落入每个矩阵元素范围内风场的个数与目标航线上的风场数据中风场的总个数计算每个矩阵元素范围内风场概率获得风场概率分布矩阵；

具体地，如图9所示，依据所下载的风场数据，提取航线上各节点的风场数据，并将该数据按照imo815通函给出的风场矩阵格式保存，将相对风向角(0°-360°)以5°为单位长度分为72段，将相对风速以1m/s为单位分为26段(最后一段为相对风速>25m/s)，则每段相对风向角和每段相对风速确定一个矩阵元素获得风场分布矩阵。然后提取并统计出每个矩阵元素位置处的风场个数，计算出每个矩阵元素处的风场个数之后除以总风场个数可得到每个矩阵元素处的风场概率。

航线内风场概率矩阵格式，依据imo815通函给出的格式统计航线上近10年间的所有风场数据综合考虑之后的总矩阵如表1所示；每个小单元表示概率水平，而不是总的样本点个数。

表1imo815通函给出风场矩阵格式

第一列为相对风速区间：m/s,[含下限，不含上限)

<1，表示[0，1)

<2，表示[1，2)

<3，表示[2，3)

…

<24，表示[23，24)

<25，表示[24，25)

≥25，表示[25，∞)

第一行为相对风向区间：°,[含下限，不含上限)

0，表示[0，5)

5，表示[5，10)

10，表示[10，15)

15，表示[15，20)

…

340，表示[340，345)

345，表示[345，350)

350，表示[350，355)

355，表示[355，360)

表格中fi,j计算公式：

例如f1,1为航线上相对风向角为[0，5)，相对风速为[0，1)的概率水平，其计算公式为：

步骤6、根据所述风场概率分布矩阵计算风帆载荷概率矩阵；

具体地，根据风帆空气动力学特性以及考虑风帆尺寸等因素，通过分析计算得到不同相对风向角下相对风速为1m/s时的风帆推力，其具体计算过程如下：

1、确定风帆信息

对于风帆助航船舶，风帆是提供辅助推进力的主要设备，能够实现主机减额输出，从而降低船舶能耗，减少碳排放，选取的风帆应具有较好的空气动力性能，对船舶运营影响较小，结构较为简单，控制方便，风帆可收放，对船舶相关设计影响较小，本专利中提到的风帆指船用硬质翼面风帆，风帆安装时需要确定的条件包括风帆的布置位置，数量，尺寸以及空气动力学特性。

1.1风帆布置位置

根据中国船级社《船用硬质翼面帆评估与检验指南2020》要求，安装用于节能目的风帆装置后，应保证船舶和船上人员的安全性，将风帆装置对船舶安全航行和人员正常操作的风险降至最低；风帆装置在正常工作和放置状态时，均应将其对船舶布置和其它设备的影响控制在尽量小的范围内；应采用公认的风险评估方法，对风帆装置引起的船舶的附加风险进行评估，并通过降低风险的措施或优化设计，将安装风帆装置船舶的总体安全风险，控制在与未安装风帆装置的同年代建造的同类型船舶相当的水平；在分析或评价风帆装置对船舶的总体性影响时，应考虑风帆装置可能导致的附加风险；发生故障时，风帆均可以在可控情况下恢复至所设定的安全状态。

总的来说，风帆布置主要考虑以下几方面：1)保持稳性的基础上尽可能增加风帆面积以实现更大的推力；2)充分利用布置空间，但尽量减少甚至不破坏船上的重要舱室；3)避开船内管路设备；4)风帆桅杆依靠船舶舱壁生长，增加刚度。如图10所示为本申请中风帆布置效果图，本发明中提到的风帆满足上述要求。

1.2风帆数量

风帆的有效使用面积受船舶宽度的制约比较明显，所以采用双帆或多帆是一种必然的有效措施，本专利可根据布置要求和风帆尺寸，确定船舶假装的风帆数量，用符号n(大于等于1)来表示加装的风帆数量。

1.3风帆的尺寸

风帆尺寸参数包括装置高度、宽度、帆体高度、桅杆高度、基座高度、风帆节数等，本专利只涉及能够提供助推效果的帆体尺寸，即风帆宽度b(m)，风帆高度h(m)，最大装帆面积a(m²)。对于可伸缩的风帆来说，还应包括风帆节数(nw)。

1.3.1风帆宽度

一般来说，现代风帆船风帆宽度基本和船宽一致。

1.3.2风帆高度

展弦比确定后，风帆帆高主要取决于帆宽和其他限制。

1.3.3最大装帆面积

船舶最大装帆面积主要由风况、船舶稳性富余及船舶具体布置情况决定。

1.4风帆的空气动力学特性

如图11所示，现代翼型风帆产生的推力与机翼十分相似，当帆翼与来风成一个夹角叫做攻角(迎角)时，帆翼背风面拱面的外面空气的流速大于迎风面拱面的内面的空气流速，帆翼背风面的流线较密，迎风面流线较疏。根据伯努利定律：流体流速大的地方流体产生的压强小流体流速小的地方流体产生的压强大。因此帆翼上面空气流速大，空气产生的压强小，帆翼下面空气流速小，空气产生的压强大，这样机翼的上下面就产生了一个压强差，此压强差产生了空气作用于机翼上的由下而上的压力。

风帆的使用对于天气海况、船舶航行姿态、自身的操纵控制等都有一定的要求。故在不适宜风帆使用的条件下，收起风帆，保证船舶的正常航行，减少风帆的额外损耗。在适宜风帆使用的情况下，需要制定一套风帆的最佳操纵策略，保证风帆效率达到最佳。

1.4.1风帆受力分析

风帆的坐标系定义和受力特性如图12所示。

坐标系原点o选在风帆轴线上，ox轴平行地面指向风帆尾缘，oz轴垂直地面指向下，oy轴根据右手法则确定。定义风帆迎流前缘向左偏转为风帆攻角α正方向；定义风帆风阻力正方向顺风指向下游，升力垂直阻力方向指向下为正，转帆力矩顺时针为正。

海面上的绝对风与船舶自身的速度合成后，形成相对风，相对风速度为v，风帆以功角α迎风时，产生垂直于来风方向的升力l和沿着风速方向的阻力d，二者沿航向的分力的合力为风帆提供的推进力，二者沿垂直航向分力的合力为风帆产生的横向力。

1.4.2风洞试验

确定风帆的展弦比、拱度比等参数后，需进行模型风洞试验，以确定风帆的气动力特性，主要是确定风帆的升力系数、阻力系数、转帆系数和合力系数。各系数可用以下无因次量系数表示：

式中cl、cd、cmz、cr为升力、阻力、转帆力矩、合力的无因次系数，ρ为空气密度，kg/m³；v为相对风速，m/s；s为风帆投影面积，m²。

通过不断改变风帆功角，分别测量风帆受力情况，最终计算出各个攻角下风帆气动力系数。功角变化范围为：-90°～+90°，间隔10°，升力极值点附近加密。根据实验可以得到cl、cd、cmz、cr与相对风向角的关系。

1.4.3风帆最佳回转角度分析

辅助船舶前进的风帆助推力是升力与阻力合成的结果，风帆的使用会增加舵角，利用横向的舵力来抵消风帆横向力的影响，维持船舶平衡，风帆的净推力是风帆助推力与保向舵力的差值。随着相对风向角变化，需要调整风帆转角，以求得最大的风帆净推力。

如图13所示，图中横坐标为风帆阻力cd，其方向和相对风向平行，纵坐标为风帆升力cl，曲线上的每一个点代表了一个风帆攻角。过该坐标原点并且沿船舶航向的为船舶所受风帆助推力，过原点垂直于船舶航向的是船舶所受风帆侧推力。由图13中几何关系可以得出风帆助推力系数ct与侧推力系数ch和风帆升力系数cl与阻力系数cd之间的关系。

ct＝clsinβ-cdcosβ(7)

ch＝clcosβ+cdsinβ(8)

其中，β为相对风向角。

因此，风帆助推力tw和风帆侧推力th也可以表示为：

可知，风帆助推力取决于相对风速和推力系数，而推力系数是由相对风向角和风帆的升力、阻力系数决定的，风帆的升力、阻力系数和风帆攻角有关。在相对风向角一定时，风帆助推力取决于风帆攻角。如图13所示，船舶漂角γ、风帆攻角α、相对风向角β和风帆回转角度之间的几何关系为：

不同的风帆攻角对应不同的风帆助推力和横向力，在计算净推力时要将风帆的横向力考虑在内，即计算保向舵力。

yr＝th(12)

yr为舵横向力，由之前的舵力模型可知：

yr＝(1+αh)fncosδ(13)

可计算当风帆处于1°至180°任意攻角时的保向舵力所对应的保向舵角。根据保向舵角即可求得舵纵向力，即保向舵力，将风帆助推力减去保向舵力即为风帆的净助推力：

tw'＝xr-rx'(15)

tw'为风帆净助推力，rx'为保向舵力。

计算在相对风向角为0-180°时，风帆攻角从0°变化至180°过程中，风帆净助推力tw'最大的风帆攻角，再求得风帆最佳操帆角，形成风帆最佳操纵策略。

根据计算所得的推力系数和阻力系数，可得阻力系数与推力系数的极坐标图，然后计算风帆回转角与相对风向角、攻角、漂角之间的关系，实际计算中，如漂角小于3°，其对船舶航行的影响较小，在忽略漂角的情况下，风帆回转角与相对风向角之间的互关系如下：

风帆助航船舶可根据相对风向角大小来旋转风帆，以达到最佳的助推效果，本专利涉及的分析均基于最佳风帆操纵策略。

根据风帆空气动力学特性以及考虑风帆尺寸等因素，通过分析计算得到不同相对风向角下相对风速为1m/s时的风帆推力，单位是kn。因此当确定了相对风向角和相对风速后，即可计算得到风帆能够提供的助推力，如表2所示(注：船舶两侧风向推力效果相同)。相对风向角的定义如图4所示：

表2风帆净推力数据表

在航线上风场概率矩阵的基础上，考虑风帆在不同风速风向下的助推力，计算得到航线上风帆的载荷概率矩阵，计算方式为以表2格式中的实际数据，与风场概率矩阵中的相对风向与相对风速相乘，计算得到风帆载荷概率矩阵。

由目标航线内可用风场概率矩阵和不同风向角下的风帆助推力，可分别算得往返航线内风帆载荷概率矩阵。

首先，将可用风场概率矩阵进行格式转换，对于每个区间上的风帆推力计算，用于插值无量纲矩阵的相对风向角，取每个区间的中间值，如

0，表示[0，5)，无量纲推力矩阵插值，使用2.5°

5，表示[5，10)，无量纲推力矩阵插值，使用7.5°

10，表示[10，15)，无量纲推力矩阵插值，使用12.5°

…

对于每个区间上的风帆推力计算，用于计算风载荷的风速大小，取区间上限，如

<1，表示[0，1)，计算风帆推力，使用1m/s

<2，表示[1，2)，计算风帆推力，使用2m/s

<3，表示[2，3)，计算风帆推力，使用3m/s

…

由风帆不同风向角下的净推力系数以及风帆的尺寸可计算得到不同风向角下的风帆产生的净推力。将其带入到更改后的可用风场概率矩阵中，可计算得到航线内风场载荷概率矩阵fij。

即表1中数据按照风向取区间中间值，风速取区间最大值的原则，将表2中风帆推力数据表中的数据，按照对应的风速风向相乘，得到风帆载荷矩阵，如下式所示：

fij＝fij×fsk(1)

式中，fij为不同风速风向时风帆的推力载荷概率，fij为风场概率，fsk为fij对应的风速风向下风帆能够提供的推力。如表3所示。

表3风翼推力载荷矩阵

同时，由于风帆的实际推力效果与风帆所处的高度有关，因此当船舶在不同装载条件下风帆所能产生的推力也不同，即表1中的推力数据不同，并且船舶在实际营运过程中，一般分为满载和压载状态航行，所以分别按照满载和压载状态分别计算得到相应的风场载荷概率矩阵，进而分别对两种航行状态进行节能评估计算。

步骤7、根据所述风帆载荷概率矩阵计算风帆总推力，并通过所述风帆总推力计算风帆的节能效率。

具体地，所述风帆总推力fm通过公式(2)进行计算：

其中：i,j为风载荷概率矩阵fij索引值，n为风载荷概率矩阵fij的数据个数，即风场载荷概率矩阵元素累加可计算得到所需要分析所有年份的风帆在该航线上的总推力fm；

所述风帆的节能效率η通过公式(3)进行计算：

其中：sv为航速，p为船舶的有效功率，即由于总推力fm与航速sv的乘积即为风帆能够产生的功率，则风帆的效率η为风帆产生的功率比上船舶的有效功率p，从而可以获得评价目标航线风场对风帆助航效果的标准。

本发明通过获取全球风场数据库，并根据全球风场数据库建立目标航线的风场概率分布矩阵，进而可以根据目标航线的风场概率分布矩阵计算风帆载荷概率矩阵，从而可以计算风帆总推力获得风帆节能效率。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。