一种基于CFD预测消浪过程的方法及其生态修复应用与流程

本发明涉及生态工程技术领域，特别涉及一种基于cfd预测消浪过程的方法及其生态修复应用。

背景技术：

湖泊生态环境修复对于改善生态环境、提高居民生活质量具有重要的意义，传统的湖泊生态修复方法可分为物理方法、化学方法和生物方法，其中化学方法难以治本，尤其对湖泊中常见的水生微生物集中爆发灾害，例如蓝藻，沿湖居民深受其害。由于湖水波浪较为湍急，对湖岸用于生态修复的人工植被等具有较强的摧毁力，人工栽培的水生植物与湖岸床底较难固定，存活时间难以满足消纳水生微生物的作用。为防止洪峰期间出现洪涝灾害，很多湖畔构筑过人工石坝，湖浪撞击石坝时对湖床会形成“马蹄涡效应”，新沉积的湖床土层难以稳定，对湖畔人工栽培的水生植物的稳固和存活极为不利。因此，通过修筑布局合理的消浪设施(例如消波坝、群桩等)，对涌进湖畔的风生浪进行扰动，撞击过程可吸收波浪动能，使湖水浪流由湍急形态转换到平缓状态，减少浪爬高，使得湖畔用于消纳有害微生物的水生植物根茎周围可以逐渐形成稳固的湖床，从而创造湖泊生态修复的基本条件。

现有湖水生态修复工程中用到的消浪方案未进行准确的定量分析，消浪设施的消浪效果无法得到准确预测，同时消浪设施在湖床的稳固性和抗浪冲击疲劳性也无法得到精确评估，最终直接导致消浪设施被湖浪击垮，人工栽培的水生植物无法长期存活，修复效果不明显或者修复失败。

技术实现要素：

针对上述技术问题，本发明提供一种基于cfd预测消浪过程的方法，包括如下步骤：

s100、拟合湖床地形数字图像：基于实测的湖床测绘图，拟合生成地形数字图像；

s200、拟合地形曲面：采集地形数字图像上的点坐标数据，拟合生成地形曲面；

s300、利用消浪设施与地形曲面构造流体域几何模型：基于地形曲面的边界线，利用几何拓扑学创建辅助曲面，封闭地形曲面、辅助曲面和消浪设施围成的区域，形成流体域几何模型；

s400、建立消浪过程流动数值计算模型：根据流体力学基本方程建立包含体积分数函数输运计算格式和水气交界自由表面重构算法的数值计算模型；

s500、对所述流体域几何模型进行网格划分，得到流体域网格模型，用于离散流体力学基本方程和体积分数函数输运方程；

s600、基于流体域网格模型，设定所述的初始条件、边界条件(含入射波)，进行数值求解，得到消浪过程的定量结果；

进一步的，所述步骤s100拟合湖床地形数字图像的步骤为：

s1001、编制拟合计算程序，读取实测湖床地形测绘图中的等高线数据；

s1002、设定拟合步长，基于所述等高线数据得到地形数字图像。

进一步的，所述边界条件包括壁面边界、流体域入口边界、流体域出口与侧面边界；

其中：

所述壁面边界：采用无滑移壁面边界条件；

所述流体域出口与侧面边界：采用静水压力边界；

所述流体域入口边界：采用人工波浪边界条件。

进一步的，所述入射波采用二阶stokes波。

进一步的，所述步骤s400中的数值计算模型包括n-s方程离散格式、体积分数函数输运模型以及自由表面重构算法模型。

进一步的，所述步骤s500中进行网格划分时采用结构化网格划分，将流体域几何模型剖分成若干个具备扫掠特征的子特征，设置网格密度控制数据，基于扫掠特征生成流体域的六面体网格。

优选的，所述六面体网格包括棱柱体网格、正方体网格、长方体网格。

进一步的，所述初始条件包括基于流体域网格模型的初始水位分布、初始流速。

进一步的，所述数值求解能够采用的方法包括有限差分法和有限体积法。

进一步的，所述方法还包括：

s700、利用步骤s600求解得到的结果，指导消浪设施本体设计。

进一步的，所述消浪设施包括消波坝、消浪墙和群桩。

优选的，所述消波坝、消浪墙和群桩的边界条件采用无滑移壁面边界条件。

优选的，所述消浪墙采用段状分布。

进一步的，所述消浪设施本体设计包括壁面压力分布模式统计、壁面压力时间历程分析、消浪墙配筋计算。

进一步的，所述方法还包括：

s800、将上述计算结果及设计结果作为生态修复的量化依据，应用于生态修复工程。

本发明方法以完备的理论模型和数值模型作为基础，综合考虑湖床地形及消浪设施干扰、湖水波浪模式等复杂因素，通过数值计算方法实现消浪过程模拟，可对消浪过程进行准确的定量预测，给出准确的消浪效果，与实际环境相符，有效避免以往生态修复工程中因对消浪设施仅做定性分析带来的弊端。通过消浪，创造湖水挟带泥沙沉淀形成稳定湖床的基本条件，同时消浪过程的计算结果可用于指导并优化消浪设施的设计，合理化的降低消浪设施的建造成本，最后将计算结果及设计结果作为生态修复的量化依据，应用于生态修复工程，使得修复效果明显，为改善生态环境提供了重要功效。

附图说明

图1为流体域几何模型示意图；

图2为第9步重构的水气交界自由表面示意图；

图3为流体域几何模型网格示意图；

图4为自由表面区域网格加密示意图；

图5为第21步流速矢量示意图；

图6为第50步流速矢量示意图；

图7为第9步消浪墙前后波浪高度分布示意图；

图8为第9步消浪墙前后波浪高度对比图；

图9为中段消浪墙在水深中点处承受的波浪压力时间历程示意图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细描述，但不作为对本发明的限定。

本发明公开了一种基于cfd预测消浪过程的方法及其生态修复应用，通过数值计算方法对消浪过程进行准确的定量计算，同时根据消浪过程的计算结果指导并优化消浪设施的设计，并将计算结果及设计结果作为生态修复的量化依据，应用于生态修复工程，使得修复效果明显，为改善生态环境提供重要功效。

本发明包括如下步骤：

s100、拟合湖床地形数字图像：基于实测的湖床测绘图，拟合生成地形数字图像。具体步骤为：

s1001、编制拟合计算程序，读取实测湖床地形测绘图中的等高线数据；

s1002、设定拟合步长，基于所述等高线数据得到地形数字图像。

s200、拟合地形曲面：采集地形数字图像上的点坐标数据，拟合生成地形曲面。

s300、利用消浪设施与地形曲面构造流体域几何模型：基于地形曲面的边界线，利用几何拓扑学创建辅助曲面，封闭地形曲面、辅助曲面和消浪设施围成的区域，形成流体域几何模型。

s400、建立消浪过程流动数值计算模型：根据流体力学基本方程建立包含体积分数函数输运计算格式和水气交界自由表面重构算法的数值计算模型。

进一步的，数值计算模型包括n-s方程离散格式、体积分数函数输运模型以及自由表面重构算法模型。

s500、对所述流体域几何模型进行网格划分，得到流体域网格模型，用于离散流体力学基本方程。

进一步的，进行网格划分时采用结构化网格划分，将流体域几何模型剖分成若干个具备扫掠特征的子特征，设置网格密度控制数据，基于扫掠特征生成流体域的六面体网格。

优选的，六面体网格包括棱柱体网格、正方体网格、长方体网格。

s600、设定所述流体域网格模型的初始条件、边界条件以及入射波，进行数值求解，得到消浪过程的定量计算。

进一步的，初始条件包括流体域网格模型的初始水位分布、初始流速。

进一步的，边界条件包括壁面边界、流体域入口边界、流体域出口与侧面边界；

其中，所述壁面边界：采用无滑移壁面边界条件；所述流体域出口与侧面边界：采用静水压力边界；所述流体域入口：采用人工波浪边界条件。

进一步的，入射波采用二阶stokes波。

进一步的，数值求解能够采用的方法包括有限差分法和有限体积法。

s700：指导消浪设施本体设计：利用步骤s600求解得到的结果，指导消浪设施本体设计。

进一步的，消浪设施包括消波坝、消浪墙和群桩。

优选的，消波坝、消浪墙和群桩的边界条件采用无滑移壁面边界条件。

优选的，消浪墙采用段状分布。

进一步的，消浪设施本体设计包括壁面压力分布模式统计、壁面压力时间历程分析、消浪墙配筋计算。

s800、将上述计算结果及设计结果作为生态修复的依据，应用于生态修复工程。

在一个实施例中，公开了某内陆湖预测消浪过程的方法及其生态修复的应用：

1、地形数字图像拟合：

编制拟合计算程序，具体步骤是：

s1001、编制拟合计算程序，读取实测湖床地形测绘图中的等高线数据；

s1002、设定拟合步长，基于所述等高线数据得到地形数字图像。

2、拟合地形曲面：

根据设定的分辨率采集地形数字图像上的点坐标数据，拟合生成地形曲面。

3、构造流体域几何模型：

基于地形曲面的边界线，利用几何拓扑学创建辅助曲面，封闭地形曲面、辅助曲面和消浪设施围成的区域，形成流体域几何模型，如图1所示。流体域包括消浪设施影响湖水流动的区域。流体域纵向长度200m，宽度100m，高度20m，消浪墙各段净距为4.0m。

4、建立消浪过程流动数值计算模型：

流动控制方程：

采用雷诺时均法对流动控制方程进行变换，变换后的连续性方程和动量方程如下所示，

式中，ρ——空气密度(kg/m³)，t——时间(s)，ui与uj——速度张量(mm/s)，xi与xj——坐标张量(mm)，p——压力(pa)，f——体积力张量(n/mm³)，τji——空气切应力张量(pa)，——雷诺应力张量(pa)。其中，上标“′”表示脉动值，上标“-”表示取平均值。采用vof方法计算水气交界自由表面体积分数函数满足输运方程，如下所示：

f为体积分数，v为流场速度向量。流体体积分数函数的控制体按下式确定：ω＝(xi-1/2，xi+1/2)×(yj-1/2，yj+1/2)×(zk-1/2，zk+1/2)×(tn，tn+1)(4)

式(3)的展开型式为：

在控制体ω内对(5)进行积分得到：

对(6)按下式展开：

水气交界自由表面由(5)中体积分数函数的界面值进行重构得到，重构采用自然邻近插值算法，插值函数如下式所示：

式中通过相邻网格点的体积分数函数值计算得到，下一时刻体积分数函数的输运通过式(7)与自由面外法线向量确定，外法线向量通过式(8)确定的曲面拟合函数在网格离散点上的值进行计算，第9时间步重构的水气交界自由表面如图2所示。

为保证自由表面的拟合精度，对初判的自由表面分布区域的网格进行细化，细化部分的网格尺寸取0.1m，详见步骤5，湖水自由表面初始高度取为10m，初始流速为0m/s，初始水位分布为平面，(1)与(2)中时间项离散采用向前一次差分，扩散项采用中心差分格式，对流项采用高阶quick格式，计算时长取10s，时间步长取0.1s。

5、进行网格划分：

采用结构化网格对流体域进行分网，将流体域的几何模型剖分成若干个具备扫掠特征的子特征，设置网格密度控制数据，基于扫掠特征生成流体域的六面体(或三棱柱)网格，网格单元数量448835。如图3、图4所示。

6、施加边界条件：

壁面边界：湖床、消浪墙、消波坝体和群桩表面采用无滑移壁面边界条件；

流体域出口与侧面边界：采用静水压力边界；

流体域入口边界：采用人工波浪边界条件；

入射波：采用二阶stokes波，其波浪势函数与波面函数如下式所示：

上式中h为波高，h为水深，l为波长，t为时间，σ为波动圆频率，x为横向坐标。

基于步骤5的分网数据，将初始条件与边界条件代入步骤4中的数值计算模型，进行数值计算。流速矢量如图5、图6所示。

如图7所示，图中右侧为消浪墙入射波侧的波浪高度，平均浪爬高约为1.213m，左侧为波浪经过消浪墙干扰后的波浪高度，平均浪爬高约为0.357m。由此可见，本实例中段状分布的消浪墙(消波坝)的消浪作用显著，波浪振幅缩减比例达到0.294，后续群桩可对消浪后的流动起到滤波作用。

消浪前后的波高对比曲线图如图8所示，wallright表示消浪前的波浪高度沿纵向的分布，wallleft表示消浪后的波浪高度沿纵向的分布，水位基准高度取8500mm。本实例中的消浪墙布置型式对典型波浪具有显著的消波作用，采用本发明中提供的数值计算格式，可对消浪过程进行准确的定量预测，为湖水挟带泥沙沉淀形成稳定湖床创造条件，经过消浪，湖畔栽培的水生植物(用于消纳湖水富营养化形态中的氮磷元素，并可吸附蓝藻、绿藻等巨害生物体)可免受巨浪持续不断的撞击，其在湖床稳固的耐久性得到有效保障。

中段消浪墙在水深中点处承受的波浪压力时间历程如图9所示。按波浪压力峰值1.6x10⁵pa计算得到每延米的消浪墙底弯矩约为8000kn*m，消浪墙墙截面厚度取1500mm，墙体竖向受力钢筋配置的计算值约为17036mm²，钢筋强度等级为hrb400。同时根据消浪墙承受的波浪压力时间历程和消浪墙本体疲劳特性(可表征为广义s-n曲线)，可进一步预测消浪墙的抗浪冲击疲劳寿命。

以上所述实施方式仅表达了本发明的一个实施方式，并不能因此而理解为对本发明范围的限制，应当指出，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。