CN108287965B - 一种计算波浪与不规则结构物相互作用的方法 - Google Patents

一种计算波浪与不规则结构物相互作用的方法 Download PDF

Info

Publication number
CN108287965B
CN108287965B CN201810070890.1A CN201810070890A CN108287965B CN 108287965 B CN108287965 B CN 108287965B CN 201810070890 A CN201810070890 A CN 201810070890A CN 108287965 B CN108287965 B CN 108287965B
Authority
CN
China
Prior art keywords
equation
grid
boundary
fluid
wave
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810070890.1A
Other languages
English (en)
Other versions
CN108287965A (zh
Inventor
李雪艳
王庆
谢天
朱小松
战超
张金芝
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ludong University
Original Assignee
Ludong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ludong University filed Critical Ludong University
Priority to CN201810070890.1A priority Critical patent/CN108287965B/zh
Publication of CN108287965A publication Critical patent/CN108287965A/zh
Application granted granted Critical
Publication of CN108287965B publication Critical patent/CN108287965B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06—Power analysis or power optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种计算波浪与不规则结构物相互作用的方法，其特点是：该方法是基于雷诺时均N‑S方程组模拟粘性波动流体，利用湍流模型封闭N‑S方程组，作为计算方法的控制方程；通过SIMPLE算法迭代调整压力，使得内部流体单元满足连续方程，自由表面单元满足自由表面动力边界条件，同时调整速度场；同传统的计算波浪与不规则结构物相互作用的方法相比较，本发明基于考虑流体粘性的N‑S方程、采用矩形网格划分计算区域，具有数值结果更符合物理实际的优点。

Description

一种计算波浪与不规则结构物相互作用的方法
技术领域：
本发明涉及防波堤技术领域，具体地讲是一种计算波浪与不规则结构物相互作用的方法；该发明可用于波浪与圆柱体、弧形防波墙、海洋平台、半圆型防波堤等海洋工程不规则结构物相互作用问题的研究。
背景技术：
随着海洋经济的快速发展和国家大力推动“海上丝绸之路”经济带建设，许多大型港口、码头、海洋平台、跨海大桥、浮式飞机场和人工岛等正逐步向海况恶劣、地质条件复杂的较深海域发展。以斜坡式和沉箱重力式为主的传统防波堤结构型式存在工程造价高、施工困难、不利于掩护海域水体交换等缺点，已无法满足海洋工程建设对防波堤结构型式及其功能的要求。近年来，专家学者们陆续开展新型防波堤结构研究，并提出各种型式的透空式和浮式防波堤。为满足不同结构与功能要求，各种新型防波堤结构的几何外形呈现不规则化。关于波浪与此类不规则结构物相互作用的研究，部分学者采用基于势流理论的边界元法(Liu et al.,2009；Ning et al.,2014；Zhou et al.,2015；王双强等,2016)开展。此方法将流体概化为无粘、无旋的理想流体，与实际流体运动的有旋性、漩涡的扩散性和能量的耗散性等存在较大不同，忽略流体粘性作用的势流模型，应用受到一定限制。另有学者采用贴体网格法(Liu,et al.,2009；Li,et al.,2011；李雪艳等,2013)开展；此方法中贴体网格生成较为复杂，且需对控制方程进行坐标变换，影响模型的计算效率。
发明内容：
本发明的目的是克服上述已有计算方法的不足，提供一种计算波浪与不规则结构物相互作用的方法；主要解决现有的方法应用受限及影响模型的计算效率等问题。
本发明的技术方案是：一种计算波浪与不规则结构物相互作用的方法，其特殊之处在于，包括以下步骤：
a根据设计的数值水槽尺寸和不规则结构物所在位置，生成计算区域的网格；所述的网格由网格单元组成，包括空网格单元、内部流体网格单元和近壁区流体网格单元；
b对计算区域的流场，设置初始条件；所述的初始条件包括速度和压强的初始值；初始条件指的是计算域内所有网格点上的速度均设置为零，全域的压强值按静水压强分布设置；
c基于雷诺时均N-S方程组模拟粘性波动流体，利用湍流模型封闭N-S方程组，作为计算方法的控制方程；所述的雷诺时均N-S方程组包括动量方程和连续方程。
d利用有限差分法对控制方程进行离散；所述的有限差分法的差分网格采用交错网格；
e设置计算区域的边界条件，所述的边界条件包括造波边界条件、开边界条件、数值水槽上下边界条件；
f将设定速度和压强的初始值代入动量方程，求出计算区域新时刻速度场的近似解；
g通过SIMPLE算法迭代调整压力，使得内部流体网格单元满足连续方程，自由表面单元满足自由表面动力边界条件，同时调整速度场；
h利用壁函数技术计算近壁区流体网格单元的紊动动能k和紊动耗散率ε；
i利用上一时刻的流体体积函数F值和已知的速度场，求出新时刻的F值，据此确定流体自由表面的位置；
j利用虚拟边界力法模拟波浪与不规则结构物之间的相互作用；
k进行数值稳定条件和收敛条件的判断，满足条件则输出压力场和速度场的结果；
l重复上述e至k步骤所述过程，直至计算时间达到程序所设定的总时间。
进一步的，所述的雷诺时均N-S方程组包括增加虚拟边界力项的水平方向时均动量方程1、竖直方向时均动量方程2和连续方程3；所述的湍流模型为k-ε模型，由K方程4和ε方程5组成；
Figure GDA0002906715430000031
Figure GDA0002906715430000032
Figure GDA0002906715430000033
Figure GDA0002906715430000034
Figure GDA0002906715430000035
其中，u为x方向的速度分量，v为y方向的速度分量，t为计算时间，gx为水平方向重力加速度，取值为零，gy为垂直方向的重力加速度，取值为9.81N/kg，p为流体压力，ρ为流体密度，ν为流体运动粘滞系数，
Figure GDA0002906715430000036
是紊动粘性系数，k为紊动动能，ε紊动耗散率；fxvbf和fyvbf分别为虚拟边界力在x和y方向的分量，θ为部分单元体参数，即结构物在整个网格单元中所占的面积与网格单元总面积的比值，范围在0～1之间；其它参数Cu＝0.09，Cε1＝1.43，σk＝1.0，σε＝0.1643，Cε2＝1.92。
进一步的，所述的交错网格是指压力与速度交错布置，压力pi,j位于网格单元中心点；
Figure GDA0002906715430000041
和
Figure GDA0002906715430000042
分别位于网格单元右边界和上边界中心点。
进一步的，所述的利用有限差分法对控制方程进行离散，水平方向时均动量方程的差分格式详见方程6：
Figure GDA0002906715430000043
其中，δt为时间步长，FUX代表水平方向对流项，FUY代表竖直方向对流项，VISX和TUBX分别代表运动粘性项和紊动粘性项，fxvbf为x方向的虚拟边界力项；
水平方向对流项边界网格点采用一阶迎风格式和二阶中心格式线性组合的偏心差分格式，见方程7和方程8：
Figure GDA0002906715430000044
其中
Figure GDA0002906715430000045
α是控制迎风差分量的参数，当α取值为零时，上述差分方程式为二阶中心差分；当α取值为1时，上述差分方程退化成一阶迎风格式；式中sign是符号函数的记号：
Figure GDA0002906715430000046
Figure GDA0002906715430000051
其中，
Figure GDA0002906715430000052
Figure GDA0002906715430000053
水平方向对流项内部网格点采用三阶迎风差分格式；
当网格右侧边界中心点的水平速度大于零时，对流项差分格式，见方程9：
Figure GDA0002906715430000054
当网格右侧边界中心点的水平速度小于零时，对流项差分格式，见方程10：
Figure GDA0002906715430000055
其中，
Figure GDA0002906715430000056
Figure GDA0002906715430000057
Figure GDA0002906715430000058
Figure GDA0002906715430000061
当网格右侧边界中心点的竖直速度v*＞0时，具体计算式见方程11；
Figure GDA0002906715430000062
当网格右侧边界中心点的竖直速度v*＜0时，具体计算式见方程12；
Figure GDA0002906715430000063
网格右侧边界中心点的竖直速度v*，可由相邻网格单元上下边界中心点的竖直速度取平均值获得；具体计算式见方程13；
Figure GDA0002906715430000064
其中，所采用各参数的具体表达式如下：
Figure GDA0002906715430000065
Figure GDA0002906715430000066
Figure GDA0002906715430000067
Figure GDA0002906715430000071
运动粘性项采用二阶中心差分格式，具体表达式，见方程14：
Figure GDA0002906715430000072
紊动粘性项，采用二阶中心差分格式，具体表达式，详见方程15：
Figure GDA0002906715430000073
进一步的，所述的利用有限差分法对控制方程进行离散，竖直方向时均动量方程的差分格式，见方程16：
Figure GDA0002906715430000074
方程中FVX，FVY，VISY，TUBY可同理写出；
K方程和ε方程采用隐式线性化处理，以保证紊动动能K和紊动耗散率ε恒为正值；
隐式线性化处理后的K方程表达式，见方程17：
Figure GDA0002906715430000075
K方程中的水平方向对流项的离散格式，见方程18：
Figure GDA0002906715430000081
其中，
Figure GDA0002906715430000082
同理可写出FKYi,j，运动粘性项和紊动粘性项见方程19和20；
Figure GDA0002906715430000083
Figure GDA0002906715430000084
隐式线性化处理后的ε方程表达式，见方程21：
Figure GDA0002906715430000085
其中，水平与垂直方向对流项，运动粘性项和紊动粘性见方程22和23；
Figure GDA0002906715430000086
Figure GDA0002906715430000087
连续方程采用中心差分格式，具体离散形式，见方程24：
Figure GDA0002906715430000088
其中，AR与AT为网格单元右侧边界与上边界可通过流体部分的面积系数；AC是网格单元的体积系数；上述连续方程的离散形式只能作为流场是否收敛的判定条件；计算过程中，为满足连续方程24，需同时对速度和压力进行调整，反复迭代。
进一步的，所述的造波边界条件由虚拟的可吸收式造波机产生，包括行进波和抵消造波板二次反射波的附加波；造波板的运动速度Um(t)的计算表达式详见方程25：
Figure GDA0002906715430000091
其中，η0为所需要的余弦波面，η为造波板前的实际波面，ω为行进波和附加波的角频率，W和L为传递函数，表达式见方程26：
Figure GDA0002906715430000092
造波边界条件的差分形式，如方程27所示：
Figure GDA0002906715430000093
进一步的，所述的数值水槽上下边界条件均设为自由可滑移边界条件；上边界条件设置见方程28，下边界条件设置，见方程29：
Figure GDA0002906715430000094
Figure GDA0002906715430000101
进一步的，所述的SIMPLE算法用于内部流体网格单元，迭代过程中的压力修正，见方程30：
δp＝-sβω (30)
其中，s为上一时刻连续方程右边的一个不为零的源项，迭代须进行到所有网格上的s不大于0.001为止；ω是为提高计算精度而在压力修正项的右端引入的压力松弛因子，取值1.7；β是与网格参数和时间步长有关的量，表达式见方程31：
Figure GDA0002906715430000102
其中
Figure GDA0002906715430000103
Figure GDA0002906715430000104
求解迭代过程中，流体单元采用的速度修正方程32如下：
Figure GDA0002906715430000105
进一步的，所述的确定流体自由表面的位置，采用流体体积方法，应用施主与受主单元计算流体体积函数F值；控制方程如33所示：
Figure GDA0002906715430000111
施主与受主单元根据其交界面上的速度方向来确定，上游单元为施主单元，用FD表示，下游单元为受主单元，用FA表示，FAD表示或是施主单元或是受主单元，由自由表面方向和流体流动方向综合确定；在时间δt内，通过边界面的F量见方程34：
ΔF＝MIN(FAD(Uδt)+CF,FDδx) (34)
其中
CF＝max[(<F>-FAD)|Uδt|-(<F>-FD)δx,0.0]
<F>＝max{FD,FAD,1.0}
自由表面方向通过在VOF方法中引入函数NF来表示，具体见下表1所示。
表1自由表面方向与函数NF的关系
Figure GDA0002906715430000112
进一步的，所述的虚拟边界力法，指的是无需布置物面边界条件，通过一组离散的边界力模拟波浪与不规则结构物之间的相互作用，具体是在时均动量方程1和2的右端项中添加一个附件力项来反映；
根据水平方向时均动量方程6的表达式，推导出水平方向虚拟边界力的计算表达式，如方程35所示：
Figure GDA0002906715430000121
同理，推导出竖直方向虚拟边界力的表达式，如方程36所示：
Figure GDA0002906715430000122
波浪在物理实际中与不规则结构物的作用力分布于其外表面，多数情况下不与离散网格单元重合；因此，计算虚拟边界力需要用到的速度值，需利用邻近结构物表面的网格单元的速度插值计算得到；以竖直方向为例，虚拟边界力的计算式如方程37所示；若结构物边界与网格单元边界中心重合，可直接利用网格单元边界中心点速度求虚拟边界力，具体见方程38；同理，可推导出不同情况下水平方向虚拟边界力的计算式；
Figure GDA0002906715430000123
其中，
Figure GDA0002906715430000124
为网格中心点距离，
Figure GDA0002906715430000125
为物面上的点距最近网格下边界中心点距离，
Figure GDA0002906715430000126
为物面上的点距最近网格上边界中心点距离，vB为物面上B点处的竖直速度，D点的竖直速度vD由上一迭代步计算获取；
Figure GDA0002906715430000131
进一步的，所述的数值稳定条件包括Courant条件见方程39、扩散稳定条件见方程40和迎风系数稳定条件见方程41；
δt＜min(δxi /|ui+1/2,j|,δyj/|vi,j+1/2|) (39)
Figure GDA0002906715430000132
Figure GDA0002906715430000133
本发明的一种计算波浪与不规则结构物相互作用的方法与已有计算方法相比具有突出的实质性特点和显著进步：1、采用基于考虑流体粘性的N-S方程描述流体运动，利用湍流模型封闭N-S方程组，作为本发明计算方法的控制方程，使计算结果更加接近物理实际；2、在所述的控制方程中增加虚拟边界力项，用来模拟波浪与不规则结构物之间的相互作用，直接在矩形网格上差分求解控制方程，且无需布置物面边界条件，使计算效率大大提高；3、基于虚拟边界力法的波浪与不规则结构物相互作用方法，可推广应用到浮式采油平台、浮式飞机场等大型不规则结构物的计算中，具有重要的工程应用价值。
附图说明：
图1是本发明的自由表面单元压力计算图；
图2是本发明的虚拟边界力法示意图；
图3是本发明的速度插值示意图；
图4是本发明的计算流程图。
具体实施方式：
为了更好的理解与实施，下面结合附图给出具体实施例详细说明本发明一种计算波浪与不规则结构物相互作用的方法；所举实施例仅用于解释本发明，并非用于限制本发明的范围。
实施例1，参见图1、2、3、4，第一步，根据设计的数值水槽尺寸和不规则结构物所在位置，生成计算区域的网格；所述的网格由网格单元组成，包括空网格单元、内部流体网格单元和近壁区流体网格单元；
第二步，对计算区域的流场，设置初始条件；所述的初始条件包括速度和压强的初始值；初始条件指的是计算域内所有网格点上的速度均设置为零，全域的压强值按静水压强分布设置；
第三步，基于雷诺时均N-S方程组模拟粘性波动流体，利用湍流模型封闭N-S方程组，作为计算方法的控制方程；所述的雷诺时均N-S方程组包括动量方程和连续方程；雷诺时均N-S方程组包括增加虚拟边界力项的水平方向时均动量方程1、竖直方向时均动量方程2和连续方程3；所述的湍流模型为k-ε模型，由K方程4和ε方程5组成；
Figure GDA0002906715430000141
Figure GDA0002906715430000142
Figure GDA0002906715430000143
Figure GDA0002906715430000144
Figure GDA0002906715430000145
其中，u为x方向的速度分量，v为y方向的速度分量，t为计算时间，gx为水平方向重力加速度，取值为零，gy为垂直方向的重力加速度，取值为9.81N/kg，p为流体压力，ρ为流体密度，ν为流体运动粘滞系数，
Figure GDA0002906715430000151
是紊动粘性系数，k为紊动动能，ε紊动耗散率；fxvbf和fyvbf分别为虚拟边界力在x和y方向的分量，θ为部分单元体参数，即结构物在整个网格单元中所占的面积与网格单元总面积的比值，范围在0～1之间；其它参数Cu＝0.09，Cε1＝1.43，σk＝1.0，σε＝0.1643，Cε2＝1.92；
第四步，利用有限差分法对控制方程进行离散；所述的有限差分法的差分网格采用交错网格；交错网格是指压力与速度交错布置，压力pi,j位于网格单元中心点；
Figure GDA0002906715430000152
和
Figure GDA0002906715430000153
分别位于网格单元右边界和上边界中心点；
利用有限差分法对控制方程进行离散，水平方向时均动量方程的差分格式详见方程6：
Figure GDA0002906715430000154
其中，δt为时间步长，FUX代表水平方向对流项，FUY代表竖直方向对流项，VISX和TUBX分别代表运动粘性项和紊动粘性项，fxvbf为x方向的虚拟边界力项；
水平方向对流项边界网格点采用一阶迎风格式和二阶中心格式线性组合的偏心差分格式，见方程7和方程8：
Figure GDA0002906715430000161
其中
Figure GDA0002906715430000162
α是控制迎风差分量的参数，当α取值为零时，上述差分方程式为二阶中心差分；当α取值为1时，上述差分方程退化成一阶迎风格式；式中sign是符号函数的记号：
Figure GDA0002906715430000163
Figure GDA0002906715430000164
其中，
Figure GDA0002906715430000165
Figure GDA0002906715430000166
水平方向对流项内部网格点采用三阶迎风差分格式；
当网格右侧边界中心点的水平速度大于零时，对流项差分格式，见方程9：
Figure GDA0002906715430000167
当网格右侧边界中心点的水平速度小于零时，对流项差分格式，见方程10：
Figure GDA0002906715430000171
其中，
Figure GDA0002906715430000172
Figure GDA0002906715430000173
Figure GDA0002906715430000174
Figure GDA0002906715430000175
当网格右侧边界中心点的竖直速度v*＞0时，具体计算式见方程11；
Figure GDA0002906715430000176
当网格右侧边界中心点的竖直速度v*＜0时，具体计算式见方程12；
Figure GDA0002906715430000177
网格右侧边界中心点的竖直速度v*，可由相邻网格单元上下边界中心点的竖直速度取平均值获得；具体计算式见方程13；
Figure GDA0002906715430000178
其中，所采用各参数的具体表达式如下：
Figure GDA0002906715430000181
Figure GDA0002906715430000182
Figure GDA0002906715430000183
Figure GDA0002906715430000184
运动粘性项采用二阶中心差分格式，具体表达式，见方程14：
Figure GDA0002906715430000185
紊动粘性项，采用二阶中心差分格式，具体表达式详见方程15：
Figure GDA0002906715430000186
利用有限差分法对控制方程进行离散，竖直方向时均动量方程的差分格式，见方程16：
Figure GDA0002906715430000187
方程中FVX，FVY，VISY，TUBY可同理写出；
K方程和ε方程采用隐式线性化处理，以保证紊动动能K和紊动耗散率ε恒为正值；
隐式线性化处理后的K方程表达式，见方程17：
Figure GDA0002906715430000191
K方程中的水平方向对流项的离散格式，见方程18：
Figure GDA0002906715430000192
其中，
Figure GDA0002906715430000193
同理可写出FKYi,j，运动粘性项和紊动粘性项见方程19和20；
Figure GDA0002906715430000194
Figure GDA0002906715430000195
隐式线性化处理后的ε方程表达式，见方程21：
Figure GDA0002906715430000196
其中，水平与垂直方向对流项，运动粘性项和紊动粘性见方程22和23；
Figure GDA0002906715430000197
Figure GDA0002906715430000201
连续方程采用中心差分格式，具体离散形式，见方程24：
Figure GDA0002906715430000202
其中，AR与AT为网格单元右侧边界与上边界可通过流体部分的面积系数；AC是网格单元的体积系数；需要指出，上述连续方程的离散形式只能作为流场是否收敛的判定条件；计算过程中，为满足连续方程24，需同时对速度和压力进行调整，反复迭代；
第五步，设置计算区域的边界条件，所述的边界条件包括造波边界条件、开边界条件、数值水槽上下边界条件；
造波边界条件由虚拟的可吸收式造波机产生，包括行进波和抵消造波板二次反射波的附加波；造波板的运动速度Um(t)为方程25：
Figure GDA0002906715430000203
其中，η0为所需要的余弦波面，η为造波板前的实际波面，ω为行进波和附加波的角频率，W和L为传递函数，表达式见方程26：
Figure GDA0002906715430000204
造波边界条件的差分形式，如方程27所示：
Figure GDA0002906715430000205
开边界条件应用波浪线性辐射条件，方程如下所示：
Figure GDA0002906715430000211
其中，Φ为速度势函数；
数值水槽上下边界条件均设为自由可滑移边界条件；上边界条件设置见方程28，下边界条件设置，见方程29：
Figure GDA0002906715430000212
Figure GDA0002906715430000213
第六步，将设定速度和压强的初始值代入动量方程，求出计算区域新时刻速度场的近似解；
第七步，通过SIMPLE算法迭代调整压力，使得内部流体网格单元满足连续方程，自由表面单元满足自由表面动力边界条件，同时调整速度场；
SIMPLE算法用于内部流体网格单元，迭代过程中的压力修正，见方程30：
δp＝-sβω (30)
其中，s为上一时刻连续方程右边的一个不为零的源项，迭代须进行到所有网格上的s不大于0.001为止；ω是为提高计算精度而在压力修正项的右端引入的压力松弛因子，取值1.7；β是与网格参数和时间步长有关的量，表达式见方程31：
Figure GDA0002906715430000221
其中
Figure GDA0002906715430000222
Figure GDA0002906715430000223
求解迭代过程中，流体单元采用的速度修正方程32如下：
Figure GDA0002906715430000224
自由表面单元其压力可由自由表面处的压力与最接近于自由表面法向方向的相邻流体单元，即插值单元的压力线性插值得到；
Figure GDA0002906715430000225
其中，
Figure GDA0002906715430000226
为插值单元的压力，
Figure GDA0002906715430000227
为自由表面处的压力，取值为零，h＝dc/d为插值因子，dc为自由表面单元中心点到插值单元中心的距离，d为自由表面到插值单元中心的距离；
第八步，利用壁函数技术计算近壁区流体网格单元的紊动动能k和紊动耗散率ε；
壁函数技术为张长高的混合长度壁函数技术；
Figure GDA0002906715430000231
Figure GDA0002906715430000232
其中，β0是与壁面粗糙度有关的常数，取值0.0005，L是特征长度，取为近壁区网格中心到壁面的距离；
第九步，利用上一时刻的流体体积函数F值和已知的速度场，求出新时刻的F值，据此确定流体自由表面的位置；
确定流体自由表面的位置，采用流体体积方法，应用施主与受主单元计算流体体积函数F值；控制方程如33所示：
Figure GDA0002906715430000233
施主与受主单元根据其交界面上的速度方向来确定，上游单元为施主单元，用FD表示，下游单元为受主单元，用FA表示，FAD表示或是施主单元或是受主单元，由自由表面方向和流体流动方向综合确定；在时间δt内，通过边界面的F变化量为方程34：
ΔF＝MIN(FAD(Uδt)+CF,FDδx) (34)
其中
CF＝max[(<F>-FAD)|Uδt|-(<F>-FD)δx,0.0]
<F>＝max{FD,FAD,1.0}
自由表面方向通过在VOF方法中引入函数NF来表示，具体见下表1：
表1自由表面方向与函数NF的关系
Figure GDA0002906715430000234
Figure GDA0002906715430000241
第十步，利用虚拟边界力法模拟波浪与不规则结构物之间的相互作用；
虚拟边界力法，指的是无需布置物面边界条件，通过一组离散的边界力模拟波浪与不规则结构物之间的相互作用，具体是在时均动量方程1和2的右端项中添加一个附件力项来反映；
根据水平方向时均动量方程6的表达式，推导出水平方向虚拟边界力的计算表达式，如方程35所示：
Figure GDA0002906715430000242
同理，推导出竖直方向虚拟边界力的表达式，如方程36所示：
Figure GDA0002906715430000243
波浪在物理实际中与不规则结构物的作用力分布于其外表面，多数情况下不与离散网格单元重合；因此，计算虚拟边界力需要用到的速度值，需利用邻近结构物表面的网格单元的速度插值计算得到；以竖直方向为例，虚拟边界力的计算式如方程37所示；若结构物边界与网格单元边界中心重合，可直接利用网格单元边界中心点速度求虚拟边界力，具体见方程38；同理，可推导出不同情况下水平方向虚拟边界力的计算式；
Figure GDA0002906715430000251
其中，
Figure GDA0002906715430000252
为网格中心点距离，
Figure GDA0002906715430000253
为物面上的点距最近网格下边界中心点距离，
Figure GDA0002906715430000254
为物面上的点距最近网格上边界中心点距离，vB为物面上B点处的竖直速度，D点的竖直速度vD由上一迭代步计算获取；
Figure GDA0002906715430000255
第十一步，进行数值稳定条件和收敛条件的判断，满足条件则输出压力场和速度场的结果；
数值稳定条件包括Courant条件见方程39、扩散稳定条件见方程40和迎风系数稳定条件见方程41；
δt＜min(δxi/|ui+1/2,j|,δyj/|vi,j+1/2|) (39)
Figure GDA0002906715430000256
Figure GDA0002906715430000257
收敛条件指的是采用SIMPLE迭代所得到的每一个网格单元的速度代入连续方程后，方程右端的源项值均小于0.001；
第十二步，重复上述第五步至第十一步，直至计算时间达到程序所设定的总时间。

Claims (9)

1.一种计算波浪与不规则结构物相互作用的方法，其特征在于，包括以下步骤：
a根据设计的数值水槽尺寸和不规则结构物所在位置，生成计算区域的网格；所述的网格由网格单元组成，包括空网格单元、内部流体网格单元和近壁区流体网格单元；
b对计算区域的流场，设置初始条件；所述的初始条件包括速度和压强的初始值；初始条件指的是计算域内所有网格点上的速度均设置为零，全域的压强值按静水压强分布设置；
c基于雷诺时均N-S方程组模拟粘性波动流体，利用湍流模型封闭N-S方程组，作为计算方法的控制方程；所述的雷诺时均N-S方程组包括动量方程和连续方程；
d利用有限差分法对控制方程进行离散；所述的有限差分法的差分网格采用交错网格；
e设置计算区域的边界条件，所述的边界条件包括造波边界条件、开边界条件、数值水槽上下边界条件；
f将设定速度和压强的初始值代入动量方程，求出计算区域新时刻速度场的近似解；
g通过SIMPLE算法迭代调整压力，使得内部流体网格单元满足连续方程，自由表面单元满足自由表面动力边界条件，同时调整速度场；
所述的SIMPLE算法用于内部流体网格单元，迭代过程中的压力修正，见方程30：
δp＝-sβω 30
其中，s为上一时刻连续方程右边的一个不为零的源项，迭代须进行到所有网格上的s不大于0.001为止；ω是为提高计算精度而在压力修正项的右端引入的压力松弛因子，取值1.7；β是与网格参数和时间步长有关的量，表达式见方程31：
Figure FDA0002935041830000021
其中
Figure FDA0002935041830000022
Figure FDA0002935041830000023
求解迭代过程中，流体单元采用的速度修正方程32如下：
Figure FDA0002935041830000024
h利用壁函数技术计算近壁区流体网格单元的紊动动能k和紊动耗散率ε；
i利用上一时刻的流体体积函数F值和已知的速度场，求出新时刻的F值，据此确定流体自由表面的位置；
j利用虚拟边界力法模拟波浪与不规则结构物之间的相互作用；
所述的虚拟边界力法，指的是无需布置物面边界条件，通过一组离散的边界力模拟波浪与不规则结构物之间的相互作用，具体是通过一个附件力项来反映；
水平方向虚拟边界力的计算表达式，如方程35所示：
Figure FDA0002935041830000031
竖直方向虚拟边界力的表达式，如方程36所示：
Figure FDA0002935041830000032
波浪在物理实际中与不规则结构物的作用力分布于其外表面，多数情况下不与离散网格单元重合；因此，计算虚拟边界力需要用到的速度值，需利用邻近结构物表面的网格单元的速度插值计算得到；以竖直方向为例，虚拟边界力的计算式如方程37所示；若结构物边界与网格单元边界中心重合，可直接利用网格单元边界中心点速度求虚拟边界力，具体见方程38；同理，可推导出不同情况下水平方向虚拟边界力的计算式；
Figure FDA0002935041830000033
其中，
Figure FDA0002935041830000034
Figure FDA0002935041830000035
为网格中心点距离，
Figure FDA0002935041830000036
为物面上的点距最近网格下边界中心点距离，
Figure FDA0002935041830000037
为物面上的点距最近网格上边界中心点距离，vB为物面上B点处的竖直速度，D点的竖直速度vD由上一迭代步计算获取；
Figure FDA0002935041830000041
k进行数值稳定条件和收敛条件的判断，满足条件则输出压力场和速度场的结果；
l重复上述e至k步骤，直至计算时间达到程序所设定的总时间。
2.根据权利要求1所述的一种计算波浪与不规则结构物相互作用的方法，其特征在于，所述的雷诺时均N-S方程组包括增加虚拟边界力项的水平方向时均动量方程1、竖直方向时均动量方程2和连续方程3；所述的湍流模型为K-ε模型，由K方程4和ε方程5组成；
Figure FDA0002935041830000042
Figure FDA0002935041830000043
Figure FDA0002935041830000044
Figure FDA0002935041830000045
Figure FDA0002935041830000046
其中，u为x方向的速度分量，v为y方向的速度分量，t为计算时间，gx为水平方向重力加速度，取值为零，gy为垂直方向的重力加速度，取值为9.81N/kg，p为流体压力，ρ为流体密度，ν为流体运动粘滞系数，
Figure FDA0002935041830000051
是紊动粘性系数，k为紊动动能，ε紊动耗散率；fxvbf和fyvbf分别为虚拟边界力在x和y方向的分量，θ为部分单元体参数，即结构物在整个网格单元中所占的面积与网格单元总面积的比值，范围在0～1之间；其它参数Cu＝0.09，Cε1＝1.43，σk＝1.0，σε＝0.1643，Cε2＝1.92。
3.根据权利要求1所述的一种计算波浪与不规则结构物相互作用的方法，其特征在于，所述的交错网格是指压力与速度交错布置，压力pi,j位于网格单元中心点；
Figure FDA0002935041830000052
和
Figure FDA0002935041830000053
分别位于网格单元右边界和上边界中心点。
4.根据权利要求1所述的一种计算波浪与不规则结构物相互作用的方法，其特征在于，所述的利用有限差分法对控制方程进行离散，水平方向时均动量方程的差分格式，见方程6：
Figure FDA0002935041830000054
其中，δt为时间步长，FUX代表水平方向对流项，FUY代表竖直方向对流项，VISX和TUBX分别代表运动粘性项和紊动粘性项，fxvbf为x方向的虚拟边界力项；
水平方向对流项边界网格点采用一阶迎风格式和二阶中心格式线性组合的偏心差分格式，见方程7和方程8：
Figure FDA0002935041830000061
其中
Figure FDA0002935041830000062
α是控制迎风差分量的参数，当α取值为零时，上述差分方程式为二阶中心差分；当α取值为1时，上述差分方程退化成一阶迎风格式；式中sign是符号函数的记号：
Figure FDA0002935041830000063
Figure FDA0002935041830000064
其中，
Figure FDA0002935041830000065
Figure FDA0002935041830000066
水平方向对流项内部网格点采用三阶迎风差分格式；
当网格右侧边界中心点的水平速度大于零时，对流项差分格式，见方程9：
Figure FDA0002935041830000067
当网格右侧边界中心点的水平速度小于零时，对流项差分格式，见方程10：
Figure FDA0002935041830000071
其中，
Figure FDA0002935041830000072
Figure FDA0002935041830000073
Figure FDA0002935041830000074
Figure FDA0002935041830000075
当网格右侧边界中心点的竖直速度v*>0时，具体计算式见方程11；
Figure FDA0002935041830000076
当网格右侧边界中心点的竖直速度v*<0时，具体计算式见方程12；
Figure FDA0002935041830000077
网格右侧边界中心点的竖直速度v*，可由相邻网格单元上下边界中心点的竖直速度取平均值获得；具体计算式见方程13；
Figure FDA0002935041830000078
其中，所采用各参数的具体表达式如下：
Figure FDA0002935041830000081
Figure FDA0002935041830000082
Figure FDA0002935041830000083
Figure FDA0002935041830000084
运动粘性项采用二阶中心差分格式，具体表达式，见方程14：
Figure FDA0002935041830000085
紊动粘性项，采用二阶中心差分格式，具体表达式，见方程15：
Figure FDA0002935041830000086
5.根据权利要求1所述的一种计算波浪与不规则结构物相互作用的方法，其特征在于，所述的利用有限差分法对控制方程进行离散，竖直方向时均动量方程的差分格式，见方程16：
Figure FDA0002935041830000087
方程中FVX，FVY，VISY，TUBY可同理写出；
K方程和ε方程采用隐式线性化处理，以保证紊动动能K和紊动耗散率ε恒为正值；
隐式线性化处理后的K方程表达式，见方程17：
Figure FDA0002935041830000091
FKX代表水平方向对流项，FKY代表竖直方向对流项，VISK和SOUK分别代表运动粘性项和紊动粘性项；
K方程中的水平方向对流项的离散格式，见方程18：
Figure FDA0002935041830000092
其中，
Figure FDA0002935041830000093
同理可写出FKYi,j；运动粘性项和紊动粘性项见方程19和20；
Figure FDA0002935041830000094
Figure FDA0002935041830000095
隐式线性化处理后的ε方程表达式，见方程21：
Figure FDA0002935041830000101
其中，水平与垂直方向对流项，运动粘性项和紊动粘性见方程22和23；
Figure FDA0002935041830000102
Figure FDA0002935041830000103
连续方程采用中心差分格式，具体离散形式，见方程24：
Figure FDA0002935041830000104
其中，AR与AT为网格单元右侧边界与上边界可通过流体部分的面积系数；AC是网格单元的体积系数；上述连续方程的离散形式只能作为流场是否收敛的判定条件；计算过程中，为满足连续方程24，需同时对速度和压力进行调整，反复迭代。
6.根据权利要求1所述的一种计算波浪与不规则结构物相互作用的方法，其特征在于，所述的造波边界条件由虚拟的可吸收式造波机产生，包括行进波和抵消造波板二次反射波的附加波；造波板的运动速度Um(t)为方程25：
Figure FDA0002935041830000105
其中，η0为所需要的余弦波面，η为造波板前的实际波面，ω为行进波和附加波的角频率，W和L为传递函数，表达式见方程26：
Figure FDA0002935041830000106
造波边界条件的差分形式，如方程27所示：
Figure FDA0002935041830000111
7.根据权利要求1所述的一种计算波浪与不规则结构物相互作用的方法，其特征在于，所述的数值水槽上下边界条件均设为自由可滑移边界条件；上边界条件设置见方程28，下边界条件设置，见方程29：
Figure FDA0002935041830000112
Figure FDA0002935041830000113
8.根据权利要求1所述的一种计算波浪与不规则结构物相互作用的方法，其特征在于，所述的确定流体自由表面的位置，采用流体体积VOF方法，应用施主与受主单元计算流体体积函数F值；控制方程如方程33所示：
Figure FDA0002935041830000114
施主与受主单元根据其交界面上的速度方向来确定，上游单元为施主单元，用FD表示，下游单元为受主单元，用FA表示，FAD表示或是施主单元或是受主单元，由自由表面方向和流体流动方向综合确定；在时间δt内，通过边界面的F变化量为方程34：
△F＝MIN(FAD(Uδt)+CF,FDδx) 34
其中
CF＝max[(<F>-FAD)|Uδt|-(〈F>-FD)δx,0.0]
〈F>＝max{FD,FAD,1.0}
自由表面方向通过在VOF方法中引入函数NF来表示。
9.根据权利要求1所述的一种计算波浪与不规则结构物相互作用的方法，其特征在于，所述的数值稳定条件包括Courant条件见方程39、扩散稳定条件见方程40和迎风系数稳定条件见方程41；
δt<min(δxi/|ui+1/2,j|,δyj/|vi,j+1/2|) 39
Figure FDA0002935041830000121
Figure FDA0002935041830000122
CN201810070890.1A 2018-01-25 2018-01-25 一种计算波浪与不规则结构物相互作用的方法 Active CN108287965B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810070890.1A CN108287965B (zh) 2018-01-25 2018-01-25 一种计算波浪与不规则结构物相互作用的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810070890.1A CN108287965B (zh) 2018-01-25 2018-01-25 一种计算波浪与不规则结构物相互作用的方法

Publications (2)

Publication Number Publication Date
CN108287965A CN108287965A (zh) 2018-07-17
CN108287965B true CN108287965B (zh) 2021-04-06

Family

ID=62835837

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810070890.1A Active CN108287965B (zh) 2018-01-25 2018-01-25 一种计算波浪与不规则结构物相互作用的方法

Country Status (1)

Country Link
CN (1) CN108287965B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109446581B (zh) * 2018-09-29 2022-12-02 中国船舶重工集团公司第七一九研究所 一种波浪作用下浮体的水动力响应的测量方法及系统
CN109373982A (zh) * 2018-10-31 2019-02-22 重庆交通大学 一种近海浅水区海床波压力的计算方法
CN110348171B (zh) * 2019-07-31 2021-01-08 山东大学 基于OpenFOAM的三维海浪与负载交互的数值仿真模型建立方法
CN111241756A (zh) * 2020-01-08 2020-06-05 大连理工大学 基于Fluent UDF的一种数值波浪水槽造波方法
CN111238770B (zh) * 2020-03-04 2021-11-12 宁波大学科学技术学院 造波机
CN111353253B (zh) * 2020-03-31 2020-12-01 温州科技职业学院 基于Flow-3D软件仿物理造波手段实现多列不同参数波浪的方法
CN111967082B (zh) * 2020-08-21 2023-12-26 中国海洋大学 一种波浪对结构物作用的模拟分析方法
CN112182995B (zh) * 2020-10-27 2021-06-04 中国海洋大学 一种粘性势流理论分析方法
CN115525978B (zh) * 2022-10-31 2023-04-14 南京航空航天大学 用于飞行器波浪情况水动力分析的速度入口数值造波方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102495968A (zh) * 2011-12-14 2012-06-13 中国科学院工程热物理研究所 重力条件下高温蓄热容器内固/液相变的数值模拟方法
CN103235878A (zh) * 2013-04-15 2013-08-07 大连理工大学 一种柔性网衣对波浪传播影响的模拟方法
WO2013147511A1 (ko) * 2012-03-27 2013-10-03 Kim Sug-Moon 3차원 방파제 시뮬레이션 시스템 및 그 시뮬레이션 방법

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102495968A (zh) * 2011-12-14 2012-06-13 中国科学院工程热物理研究所 重力条件下高温蓄热容器内固/液相变的数值模拟方法
WO2013147511A1 (ko) * 2012-03-27 2013-10-03 Kim Sug-Moon 3차원 방파제 시뮬레이션 시스템 및 그 시뮬레이션 방법
CN103235878A (zh) * 2013-04-15 2013-08-07 大连理工大学 一种柔性网衣对波浪传播影响的模拟方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于BFC_VOF方法的弧形防浪墙水动力数值模拟;李雪艳等;《水道港口》;20110831;第32卷(第4期);第270-275页 *

Also Published As

Publication number Publication date
CN108287965A (zh) 2018-07-17

Similar Documents

Publication Publication Date Title
CN108287965B (zh) 一种计算波浪与不规则结构物相互作用的方法
CN110008509B (zh) 一种考虑背景流场下的内孤立波作用力特性分析方法
Yuan et al. An implicit three‐dimensional fully non‐hydrostatic model for free‐surface flows
CN107977542B (zh) 一种波浪与弧板式防波堤相互作用的计算方法
CN107992715B (zh) 一种基于虚拟边界力法的弧板式防波堤受力计算方法
Schiller et al. Modeling river plume dynamics with the HYbrid Coordinate Ocean Model
Ravensbergen et al. A variational multiscale framework for atmospheric turbulent flows over complex environmental terrains
CN103136436B (zh) 考虑泥沙吸附的核素输移扩散分析方法
Zhang et al. Three-dimensional simulation of horseshoe vortex and local scour around a vertical cylinder using an unstructured finite-volume technique
Dippner A frontal-resolving model for the German Bight
CN111241756A (zh) 基于Fluent UDF的一种数值波浪水槽造波方法
CN108170993B (zh) 一种弧板式防波堤结构受力的计算方法
Xu et al. Flow pattern and anti-silt measures of straight-edge forebay in large pump stations.
CN107256312A (zh) 一种基于潮流环境下海湾纳潮变化量计算方法
Jifu et al. Influences of the waterway project in the Yangtze River Estuary on sediment transport
Chen et al. On rotational flows with discontinuous vorticity beneath steady water waves near stagnation
Zou Coastal sediment transport simulation by smoothed particle hydrodynamics
Kim et al. Direct numerical study of air layer drag reduction phenomenon over a backward-facing step
Li et al. Numerical simulation of hydrodynamic characteristics on an arc crown wall using volume of fluid method based on BFC
Andersson et al. Effect of spatial resolution of rough surfaces on numerically computed flow fields with application to hydraulic engineering
Kim Three dimensional sediment transport model.
Guo et al. A new depth-integrated non-hydrostatic model for free surface flows
Xiangju et al. Progress in numerical simulation of high entrained air-water two-phase flow
Liu et al. Numerical simulation of atomized flow diffusion in deep and narrow goeges
LU et al. Numerical simulation and analysis for combinational effects of two bridges on the tidal bore in the Qiantang River

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB03 Change of inventor or designer information

Inventor after: Li Xueyan

Inventor after: Wang Qing

Inventor after: Xie Tian

Inventor after: Zhu Xiaosong

Inventor after: Chao chao

Inventor after: Zhang Jinzhi

Inventor before: Li Xueyan

Inventor before: Wang Qing

Inventor before: Zhu Xiaosong

Inventor before: Chao chao

Inventor before: Zhang Jinzhi

CB03 Change of inventor or designer information
GR01 Patent grant
GR01 Patent grant
PE01 Entry into force of the registration of the contract for pledge of patent right

Denomination of invention: A method for calculating the interaction between waves and irregular structures

Effective date of registration: 20211216

Granted publication date: 20210406

Pledgee: Yantai financing guarantee Group Co.,Ltd.

Pledgor: LUDONG University

Registration number: Y2021980015152

PE01 Entry into force of the registration of the contract for pledge of patent right
PC01 Cancellation of the registration of the contract for pledge of patent right

Date of cancellation: 20220317

Granted publication date: 20210406

Pledgee: Yantai financing guarantee Group Co.,Ltd.

Pledgor: LUDONG University

Registration number: Y2021980015152

PC01 Cancellation of the registration of the contract for pledge of patent right
EE01 Entry into force of recordation of patent licensing contract

Application publication date: 20180717

Assignee: Blue coast environmental protection technology (Dalian) Co.,Ltd.

Assignor: LUDONG University

Contract record no.: X2022370000008

Denomination of invention: A method for calculating the interaction between waves and irregular structures

Granted publication date: 20210406

License type: Common License

Record date: 20220610

EE01 Entry into force of recordation of patent licensing contract