CN107992715B - 一种基于虚拟边界力法的弧板式防波堤受力计算方法 - Google Patents

一种基于虚拟边界力法的弧板式防波堤受力计算方法 Download PDF

Info

Publication number
CN107992715B
CN107992715B CN201810070914.3A CN201810070914A CN107992715B CN 107992715 B CN107992715 B CN 107992715B CN 201810070914 A CN201810070914 A CN 201810070914A CN 107992715 B CN107992715 B CN 107992715B
Authority
CN
China
Prior art keywords
equation
grid
term
boundary
virtual boundary
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
CN201810070914.3A
Other languages
English (en)
Other versions
CN107992715A (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.)
Yantai Junfeng Software Technology Co.,Ltd.
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 CN201810070914.3A priority Critical patent/CN107992715B/zh
Publication of CN107992715A publication Critical patent/CN107992715A/zh
Application granted granted Critical
Publication of CN107992715B publication Critical patent/CN107992715B/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]

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)
  • Revetment (AREA)
  • Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)

Abstract

本发明提供一种基于虚拟边界力法的弧板式防波堤受力计算方法，该方法是基于雷诺时均N‑S方程组模拟考虑粘性的波浪运动，利用湍流模型封闭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根据速度、压强初始值，求出新时刻流场近似解，通过SIMPLE算法迭代调整压力，使得内部流体单元满足连续方程，自由表面单元满足自由表面动力边界条件，同时调整速度场；
g利用上一时刻的流体体积函数F值和已知的速度场，求出新时刻的F值，据此确定流体自由表面的位置；
h利用虚拟边界力法模拟波浪与弧板式防波堤之间的相互作用；
i判断计算过程是否满足数值稳定条件和收敛条件，若满足，则输出压力场和速度场结果；
j重复上述e至i步骤所述过程，直至计算时间达到程序所设定的总时间。
进一步的，所述的雷诺时均N-S方程组包括增加虚拟边界力项的水平方向时均动量方程1、竖直方向时均动量方程2和连续方程3；所述的湍流模型为K-ε模型，由K方程4和ε方程5组成；
Figure GDA0002759564310000031
Figure GDA0002759564310000032
Figure GDA0002759564310000033
Figure GDA0002759564310000034
Figure GDA0002759564310000035
其中，u为x方向的速度分量，v为y方向的速度分量，gx为水平方向重力加速度，取值为零，gy为垂直方向的重力加速度，取值为9.81N/kg，p为流体压力，ρ为流体密度，ν为流体运动粘滞系数，
Figure GDA0002759564310000036
是紊动粘性系数；fxvbf和fyvbf分别为虚拟边界力在x和y方向的分量，θ为部分单元体参数，即结构物在整个网格单元中所占的面积与网格单元总面积的比值，范围在0～1之间；其它参数Cu＝0.09，Cε1＝1.43，σk＝1.0，σε＝0.1643，Cε2＝1.92。
进一步的，所述的同位网格是所有参变量均定义在网格单元的中心点，包括压力pi,j、流体体积函数Fi,j、紊动动能Ki,j、紊动耗散率εi,j、水平方向速度ui,j、竖直方向速度vi,j。
进一步的，所述的利用有限差分法对控制方程进行离散，水平方向时均动量方程的差分格式，见方程6：
Figure GDA0002759564310000041
其中，FUX代表水平方向对流项，FUY代表竖直方向对流项，VISX和TUBX分别代表运动粘性项和紊动粘性项，fxvbf为x方向的虚拟边界力项；
水平方向对流项边界网格点采用一阶迎风格式和二阶中心格式线性组合的偏心差分格式，见方程7和方程8；
Figure GDA0002759564310000042
其中
Figure GDA0002759564310000043
α是控制迎风差分量的参数，当α取值为零时，上述差分方程式为二阶中心差分；当α取值为1时，上述差分方程退化成一阶迎风格式；式中sign是符号函数的记号：
Figure GDA0002759564310000044
Figure GDA0002759564310000045
其中，
Figure GDA0002759564310000046
Figure GDA0002759564310000051
水平方向对流项内部网格点采用三阶迎风差分格式；
当网格右侧边界中心点的水平速度大于零时，对流项差分格式，见方程9：
Figure GDA0002759564310000052
当网格右侧边界中心点的水平速度小于零时，对流项差分格式，见方程10：
Figure GDA0002759564310000053
其中,
Figure GDA0002759564310000054
Figure GDA0002759564310000055
Figure GDA0002759564310000056
Figure GDA0002759564310000057
当网格右侧边界中心点的竖直速度v*＞0时，具体计算式见方程11；
Figure GDA0002759564310000058
当网格右侧边界中心点的竖直速度v*＜0时，具体计算式见方程12；
Figure GDA0002759564310000059
网格右侧边界中心点的竖直速度v*，可由相邻网格单元上下边界中心点的竖直速度取平均值获得；具体计算式见方程13；
Figure GDA0002759564310000061
其中，所采用各参数的具体表达式如下：
Figure GDA0002759564310000062
运动粘性项采用二阶中心差分格式，具体表达式见方程14：
Figure GDA0002759564310000063
紊动粘性项，采用二阶中心差分格式，具体表达式见方程15：
Figure GDA0002759564310000064
竖直方向时均动量方程的差分格式见方程16：
Figure GDA0002759564310000071
方程中FVX，FVY，VISY，TUBY可同理写出；
K方程和ε方程采用隐式线性化处理，以保证紊动动能K和紊动耗散率ε恒为正值；
隐式线性化处理后的K方程表达式详见方程17：
Figure GDA0002759564310000072
K方程中的水平方向对流项的离散格式，见方程18：
Figure GDA0002759564310000073
其中，
Figure GDA0002759564310000074
同理可写出FKYi,j，运动粘性项和紊动粘性项见方程19和20；
Figure GDA0002759564310000075
Figure GDA0002759564310000076
隐式线性化处理后的ε方程表达式详见方程21：
Figure GDA0002759564310000077
其中，水平与垂直方向对流项，运动粘性项和紊动粘性见方程22和23；
Figure GDA0002759564310000081
Figure GDA0002759564310000082
连续方程采用中心差分格式，具体离散形式详见方程24：
Figure GDA0002759564310000083
其中，AR与AT为网格单元右侧边界与上边界可通过流体部分的面积系数；AC是网格单元的体积系数；上述连续方程的离散形式只能作为流场是否收敛的判定条件；计算过程中，为满足连续方程24，需同时对速度和压力进行调整，反复迭代。
进一步的，所述的数值水槽上下边界条件均设为自由可滑移边界条件；上边界条件设置详见方程25，下边界条件设置详见方程26；
Figure GDA0002759564310000084
Figure GDA0002759564310000085
进一步的，所述的SIMPLE算法用于流体单元，迭代过程中的压力修正方程表达式为方程27：
δp＝-sβω (27)
其中，s为上一时刻连续方程右边的一个不为零的源项，迭代须进行到所有网格上的s不大于0.001为止；ω是为提高计算精度而在压力修正项的右端引入的压力松弛因子，取值1.7；β是与网格参数和时间步长有关的量，表达式见方程28：
Figure GDA0002759564310000091
其中
Figure GDA0002759564310000092
求解迭代过程中，流体单元采用的速度修正方程29如下：
Figure GDA0002759564310000093
进一步的，所述的虚拟边界力法，指的是无需布置物面边界条件，通过一组离散的边界力模拟波浪与不规则结构物之间的相互作用，具体是在时均动量方程1和方程2的右端项中添加一个附件力项来反映；
根据水平方向时均动量方程6的表达式，推导出水平方向虚拟边界力的计算表达式如方程30所示：
Figure GDA0002759564310000101
同理，推导出竖直方向虚拟边界力的表达式如方程31所示：
Figure GDA0002759564310000102
波浪在物理实际中与不规则结构物的作用力分布于其外表面，多数情况下不与离散网格单元重合；因此，计算虚拟边界力需要用到的速度值，需利用邻近结构物表面的网格单元的速度插值计算得到；以竖直方向为例，虚拟边界力的计算式如方程32所示；若结构物边界与网格单元中心重合，可直接利用网格单元中心点速度求虚拟边界力，具体见方程式33；同理，可推导出不同情况下水平方向虚拟边界力的计算式；
Figure GDA0002759564310000103
其中，
Figure GDA0002759564310000104
为物面上B点处的竖直速度，D点的竖直速度vD由上一迭代步计算获取；
Figure GDA0002759564310000105
本发明的一种基于虚拟边界力法的弧板式防波堤受力计算方法与已有计算方法相比具有突出的实质性特点和显著进步：1、所采用的控制方程考虑流体粘性，数值计算结果更加接近物理实际；2、采用虚拟边界力法模拟波浪与弧板式防波堤之间的相互作用，直接在矩形网格上差分求解控制方程，可有效提高计算效率。
附图说明：
图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 GDA0002759564310000111
Figure GDA0002759564310000121
Figure GDA0002759564310000122
Figure GDA0002759564310000123
Figure GDA0002759564310000124
其中，u为x方向的速度分量，v为y方向的速度分量，t为计算时间，gx为水平方向重力加速度，取值为零，gy为垂直方向的重力加速度，取值为9.81N/kg，p为流体压力，ρ为流体密度，ν为流体运动粘滞系数，
Figure GDA0002759564310000125
是紊动粘性系数，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、流体体积函数Fi,j、紊动动能Ki,j、紊动耗散率εi,j、水平方向速度ui,j、竖直方向速度vi,j；
利用有限差分法对控制方程进行离散，水平方向时均动量方程的差分格式，见方程6：
Figure GDA0002759564310000131
其中，FUX代表水平方向对流项，FUY代表竖直方向对流项，VISX和TUBX分别代表运动粘性项和紊动粘性项，fxvbf为x方向的虚拟边界力项；
水平方向对流项边界网格点采用一阶迎风格式和二阶中心格式线性组合的偏心差分格式，见方程7和方程8；
Figure GDA0002759564310000132
其中
Figure GDA0002759564310000133
α是控制迎风差分量的参数，当α取值为零时，上述差分方程式为二阶中心差分；当α取值为1时，上述差分方程退化成一阶迎风格式；式中sign是符号函数的记号：
Figure GDA0002759564310000134
Figure GDA0002759564310000135
其中，
Figure GDA0002759564310000136
Figure GDA0002759564310000137
水平方向对流项内部网格点采用三阶迎风差分格式；
当网格右侧边界中心点的水平速度大于零时，对流项差分格式，见方程9：
Figure GDA0002759564310000141
当网格右侧边界中心点的水平速度小于零时，对流项差分格式，见方程10：
Figure GDA0002759564310000142
其中，
Figure GDA0002759564310000143
Figure GDA0002759564310000144
当网格右侧边界中心点的竖直速度v*＞0时，具体计算式见方程11；
Figure GDA0002759564310000145
当网格右侧边界中心点的竖直速度v*＜0时，具体计算式见方程12；
Figure GDA0002759564310000146
网格右侧边界中心点的竖直速度v*，可由相邻网格单元上下边界中心点的竖直速度取平均值获得；具体计算式见方程13；
Figure GDA0002759564310000151
其中，所采用各参数的具体表达式如下：
Figure GDA0002759564310000152
运动粘性项采用二阶中心差分格式，具体表达式见方程14：
Figure GDA0002759564310000153
紊动粘性项，采用二阶中心差分格式，具体表达式见方程15：
Figure GDA0002759564310000154
竖直方向时均动量方程的差分格式见方程16：
Figure GDA0002759564310000155
方程中FVX，FVY，VISY，TUBY可同理写出；
K方程和ε方程采用隐式线性化处理，以保证紊动动能K和紊动耗散率ε恒为正值；
隐式线性化处理后的K方程表达式详见方程17：
Figure GDA0002759564310000161
K方程中的水平方向对流项的离散格式，见方程18：
Figure GDA0002759564310000162
其中，
Figure GDA0002759564310000163
同理可写出FKYi,j，运动粘性项和紊动粘性项见方程19和20；
Figure GDA0002759564310000164
Figure GDA0002759564310000165
隐式线性化处理后的ε方程表达式详见方程21：
Figure GDA0002759564310000166
其中，水平与垂直方向对流项，运动粘性项和紊动粘性见方程22和23；
Figure GDA0002759564310000171
Figure GDA0002759564310000172
连续方程采用中心差分格式，具体离散形式详见方程24：
Figure GDA0002759564310000173
其中，AR与AT为网格单元右侧边界与上边界可通过流体部分的面积系数；AC是网格单元的体积系数；上述连续方程的离散形式只能作为流场是否收敛的判定条件；计算过程中，为满足连续方程24，需同时对速度和压力进行调整，反复迭代；
第五步，设置计算区域的造波边界条件、开边界条件、数值水槽上下边界条件；
造波边界条件由虚拟的可吸收式造波机产生，包括行进波和抵消造波板二次反射波的附加波；造波板的运动速度Um(t)的计算表达式见下述方程：
Figure GDA0002759564310000174
其中，η0为所需要的余弦波面，η为造波板前的实际波面，ω为行进波和附加波的角频率，W和L为传递函数，表达式如下：
Figure GDA0002759564310000175
造波边界条件的差分形式如下所示：
Figure GDA0002759564310000181
开边界条件应用波浪线性辐射条件，方程如下所示：
Figure GDA0002759564310000182
其中，Φ为速度势函数；
数值水槽上下边界条件均设为自由可滑移边界条件；上边界条件设置见方程25，下边界条件设置，见方程26：
Figure GDA0002759564310000183
Figure GDA0002759564310000184
第六步，根据速度、压强初始值，求出新时刻流场近似解，通过SIMPLE算法迭代调整压力，使得内部流体单元满足连续方程，自由表面单元满足自由表面动力边界条件，同时调整速度场；
SIMPLE算法用于流体单元，迭代过程中的压力修正方程表达式为方程27：
δp＝-sβω (27)
其中，s为上一时刻连续方程右边的一个不为零的源项，迭代须进行到所有网格上的s不大于0.001为止；ω是为提高计算精度而在压力修正项的右端引入的压力松弛因子，取值1.7；β是与网格参数和时间步长有关的量，表达式见方程28：
Figure GDA0002759564310000191
其中
Figure GDA0002759564310000192
求解迭代过程中，流体单元采用的速度修正方程29如下：
Figure GDA0002759564310000193
第七步，利用上一时刻的流体体积函数F值和已知的速度场，求出新时刻的F值，据此确定流体自由表面的位置；
采用流体体积方法，应用施主与受主单元计算流体体积函数F值；控制方程如下所示：
Figure GDA0002759564310000194
施主与受主单元根据其交界面上的速度方向来确定，上游单元为施主单元，用FD表示，下游单元为受主单元，用FA表示，FAD表示或是施主单元或是受主单元，由自由表面方向和流体流动方向综合确定；在时间δt内，通过边界面的F变化量如下所示：
ΔF＝MIN(FAD(Uδt)+CF,FDδx)
其中，
Figure GDA0002759564310000201
自由表面方向通过在VOF方法中引入函数NF来表示，具体见表1所示；
表1自由表面方向与函数NF的关系
Figure GDA0002759564310000202
第八步，利用虚拟边界力法模拟波浪与弧板式防波堤之间的相互作用；
虚拟边界力法，指的是无需布置物面边界条件，通过一组离散的边界力模拟波浪与不规则结构物之间的相互作用，具体是在时均动量方程1和方程2的右端项中添加一个附件力项来反映；
根据水平方向时均动量方程6的表达式，推导出水平方向虚拟边界力的计算表达式如方程30所示：
Figure GDA0002759564310000203
同理，推导出竖直方向虚拟边界力的表达式如方程31所示：
Figure GDA0002759564310000211
波浪在物理实际中与不规则结构物的作用力分布于其外表面，多数情况下不与离散网格单元重合；因此，计算虚拟边界力需要用到的速度值，需利用邻近结构物表面的网格单元的速度插值计算得到；以竖直方向为例，虚拟边界力的计算式如方程32所示；若结构物边界与网格单元中心重合，可直接利用网格单元中心点速度求虚拟边界力，具体见方程式33；同理，可推导出不同情况下水平方向虚拟边界力的计算式；
Figure GDA0002759564310000212
其中，
Figure GDA0002759564310000213
为物面上B点处的竖直速度，D点的竖直速度vD由上一迭代步计算获取；
Figure GDA0002759564310000214
第九步，判断计算过程是否数值稳定条件和收敛条件，若满足，则输出压力场和速度场结果；
数值稳定条件包括Courant条件、扩散稳定条件和迎风系数稳定条件分别如下所示：
δt＜min(δxi/|ui+1,j|,δyj/|vi,j+1|)
Figure GDA0002759564310000221
第十步，重复上述第五步至第九步步骤所述过程，直至计算时间达到程序所设定的总时间。

Claims (6)

1.一种基于虚拟边界力法的弧板式防波堤受力计算方法，其特征在于，包括以下步骤：
a根据所建立的数值水槽尺寸和弧板式防波堤结构所在的位置，采用笛卡尔坐标系统生成计算区域的网格；沿水槽水平和垂直方向均匀划分网格，水深方向的网格总长度大于水深与最大波高之和，且留有一定数量的空网格，以确保波动水面不受上边界限制；
b初始化计算区域的数值参数，计算区域所有网格点上的速度均设置为零，全域的压强值按静水压强分布设置；
c基于雷诺时均N-S方程组模拟考虑流体粘性的波浪运动，利用湍流模型封闭N-S方程组，作为计算方法的控制方程；所述的雷诺时均N-S方程组包括动量方程和连续方程；
d利用有限差分法对控制方程进行离散；所述的有限差分法的差分网格采用同位网格；
e设置计算区域的造波边界条件、开边界条件、数值水槽上下边界条件；
f根据速度、压强初始值，求出新时刻流场近似解，通过SIMPLE算法迭代调整压力，使得内部流体单元满足连续方程，自由表面单元满足自由表面动力边界条件，同时调整速度场；
g利用上一时刻的流体体积函数F值和已知的速度场，求出新时刻的F值，据此确定流体自由表面的位置；
h利用虚拟边界力法模拟波浪与弧板式防波堤之间的相互作用；
i判断计算过程是否满足数值稳定条件和收敛条件，若满足，则输出压力场和速度场结果；
j重复上述e至i步骤，直至计算时间达到程序所设定的总时间；
其中，c步骤所述的雷诺时均N-S方程组包括增加虚拟边界力项的水平方向时均动量方程(1)、竖直方向时均动量方程(2)和连续方程(3)，
Figure FDA0002906714980000021
Figure FDA0002906714980000022
Figure FDA0002906714980000023
u为x方向的速度分量，v为y方向的速度分量，t为计算时间，gx为水平方向重力加速度，取值为零，gy为垂直方向的重力加速度，取值为9.81N/kg，p为流体压力，ρ为流体密度，ν为流体运动粘滞系数，
Figure FDA0002906714980000024
是紊动粘性系数，k为紊动动能，ε紊动耗散率；fxvbf和fyvbf分别为虚拟边界力在x和y方向的分量，θ为部分单元体参数，即结构物在整个网格单元中所占的面积与网格单元总面积的比值，范围在0～1之间；
其中，d步骤所述的利用有限差分法对控制方程进行离散，水平方向时均动量方程的差分格式详见方程(6)：
Figure FDA0002906714980000025
其中，δt为时间步长，FUX代表水平方向对流项，FUY代表竖直方向对流项，VISX和TUBX分别代表运动粘性项和紊动粘性项，fxvbf为x方向的虚拟边界力项；
其中，h步骤所述的虚拟边界力法，指的是无需布置物面边界条件，通过一组离散的边界力模拟波浪与不规则结构物之间的相互作用具体是在时均动量方程(1)和方程(2)的右端项中添加一个附件力项来反映；
根据水平方向时均动量方程(6)的表达式，推导出水平方向虚拟边界力的计算表达式如方程(30)所示：
Figure FDA0002906714980000031
同理，推导出竖直方向虚拟边界力的表达式如方程(31)所示：
Figure FDA0002906714980000032
波浪在物理实际中与不规则结构物的作用力分布于其外表面，多数情况下不与离散网格单元重合；因此，计算虚拟边界力需要用到的速度值，需利用邻近结构物表面的网格单元的速度插值计算得到；以竖直方向为例，虚拟边界力的计算式如方程(32)所示；若结构物边界与网格单元中心重合，可直接利用网格单元中心点速度求虚拟边界力，具体见方程式(33)；同理，可推导出不同情况下水平方向虚拟边界力的计算式；
Figure FDA0002906714980000041
其中，
Figure FDA0002906714980000042
Figure FDA0002906714980000043
为网格中心点距离，
Figure FDA0002906714980000044
为物面上的点距最近网格下边界中心点距离，
Figure FDA0002906714980000045
为物面上的点距最近网格上边界中心点距离，vB为物面上B点处的竖直速度，D点的竖直速度vD由上一迭代步计算获取；
Figure FDA0002906714980000046
2.根据权利要求1所述的一种基于虚拟边界力法的弧板式防波堤受力计算方法，其特征在于，所述的湍流模型为k-ε模型，由K方程4和ε方程(5)组成；
Figure FDA0002906714980000047
Figure FDA0002906714980000048
其中，其它参数Cu＝0.09，Cε1＝1.43，σk＝1.0，σε＝0.1643，Cε2＝1.92。
3.根据权利要求1所述的一种基于虚拟边界力法的弧板式防波堤受力计算方法，其特征在于，所述的同位网格是所有参变量均定义在网格单元的中心点，包括压力pi,j、流体体积函数Fi,j、紊动动能Ki,j、紊动耗散率εi,j、水平方向速度ui,j、竖直方向速度vi,j。
4.根据权利要求1所述的一种基于虚拟边界力法的弧板式防波堤受力计算方法，其特征在于，所述的利用有限差分法对控制方程进行离散，水平方向时均动量方程的差分格式详见方程(6)：
Figure FDA0002906714980000051
其中，δt为时间步长，FUX代表水平方向对流项，FUY代表竖直方向对流项，VISX和TUBX分别代表运动粘性项和紊动粘性项，fxvbf为x方向的虚拟边界力项；
水平方向对流项边界网格点采用一阶迎风格式和二阶中心格式线性组合的偏心差分格式，见方程(7)和方程(8)：
Figure FDA0002906714980000052
其中
Figure FDA0002906714980000053
α是控制迎风差分量的参数，当α取值为零时，上述差分方程式为二阶中心差分；当α取值为1时，上述差分方程退化成一阶迎风格式；式中sign是符号函数的记号：
Figure FDA0002906714980000054
Figure FDA0002906714980000061
其中，
Figure FDA0002906714980000062
Figure FDA0002906714980000063
水平方向对流项内部网格点采用三阶迎风差分格式；
当网格右侧边界中心点的水平速度大于零时，对流项差分格式详见方程(9)：
Figure FDA0002906714980000064
当网格右侧边界中心点的水平速度小于零时，对流项差分格式详见方程(10)：
Figure FDA0002906714980000065
其中，
Figure FDA0002906714980000066
Figure FDA0002906714980000067
Figure FDA0002906714980000068
Figure FDA0002906714980000069
当网格右侧边界中心点的竖直速度v*＞0时，具体计算式见方程(11)；
Figure FDA0002906714980000071
当网格右侧边界中心点的竖直速度v*＜0时，具体计算式见方程(12)；
Figure FDA0002906714980000072
网格右侧边界中心点的竖直速度v*，可由相邻网格单元上下边界中心点的竖直速度取平均值获得；具体计算式见方程(13)；
Figure FDA0002906714980000073
其中，所采用各参数的具体表达式如下：
Figure FDA0002906714980000074
Figure FDA0002906714980000075
Figure FDA0002906714980000076
Figure FDA0002906714980000077
运动粘性项采用二阶中心差分格式，具体表达式，见方程(14)：
Figure FDA0002906714980000081
紊动粘性项，采用二阶中心差分格式，具体表达式，见方程(15)：
Figure FDA0002906714980000082
竖直方向时均动量方程的差分格式，见方程(16)：
Figure FDA0002906714980000083
方程中FVX，FVY，VISY，TUBY可同理写出；
K方程和ε方程采用隐式线性化处理，以保证紊动动能K和紊动耗散率ε恒为正值；
隐式线性化处理后的K方程表达式，见方程(17)：
Figure FDA0002906714980000084
K方程中的水平方向对流项的离散格式，见方程(18)：
Figure FDA0002906714980000085
其中，
Figure FDA0002906714980000086
同理可写出FKYi,j；运动粘性项和紊动粘性项见方程(19)和(20)；
Figure FDA0002906714980000091
Figure FDA0002906714980000092
隐式线性化处理后的ε方程表达式详见方程(21)：
Figure FDA0002906714980000093
其中，水平与垂直方向对流项，运动粘性项和紊动粘性见方程(22)和(23)；
Figure FDA0002906714980000094
Figure FDA0002906714980000095
连续方程采用中心差分格式，具体离散形式，见方程(24)：
Figure FDA0002906714980000096
其中，AR与AT为网格单元右侧边界与上边界可通过流体部分的面积系数；AC是网格单元的体积系数；上述连续方程的离散形式只能作为流场是否收敛的判定条件；计算过程中，为满足连续方程(24)，需同时对速度和压力进行调整，反复迭代。
5.根据权利要求1所述的一种基于虚拟边界力法的弧板式防波堤受力计算方法，其特征在于，所述的数值水槽上下边界条件均设为自由可滑移边界条件；上边界条件设置详见方程(25)，下边界条件设置详见方程(26)；
Figure FDA0002906714980000101
Figure FDA0002906714980000102
6.根据权利要求1所述的一种基于虚拟边界力法的弧板式防波堤受力计算方法，其特征在于，所述的SIMPLE算法用于流体单元，迭代过程中的压力修正方程为：
δp＝-sβω (27)
其中，s为上一时刻连续方程右边的一个不为零的源项，
迭代须进行到所有网格上的s不大于0.001为止；ω是为提高计算精度而在压力修正项的右端引入的压力松弛因子，取值1.7；β是与网格参数和时间步长有关的量，表达式见方程(28)：
Figure FDA0002906714980000103
其中，
Figure FDA0002906714980000104
Figure FDA0002906714980000111
求解迭代过程中，流体单元采用的速度修正方程(29)如下：
Figure FDA0002906714980000112
CN201810070914.3A 2018-01-25 2018-01-25 一种基于虚拟边界力法的弧板式防波堤受力计算方法 Active CN107992715B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810070914.3A CN107992715B (zh) 2018-01-25 2018-01-25 一种基于虚拟边界力法的弧板式防波堤受力计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810070914.3A CN107992715B (zh) 2018-01-25 2018-01-25 一种基于虚拟边界力法的弧板式防波堤受力计算方法

Publications (2)

Publication Number Publication Date
CN107992715A CN107992715A (zh) 2018-05-04
CN107992715B true CN107992715B (zh) 2021-03-23

Family

ID=62039789

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810070914.3A Active CN107992715B (zh) 2018-01-25 2018-01-25 一种基于虚拟边界力法的弧板式防波堤受力计算方法

Country Status (1)

Country Link
CN (1) CN107992715B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109408863B (zh) * 2018-09-11 2020-08-21 中山大学 波浪与防波堤相互作用的预报方法、装置和计算机设备
CN110232229B (zh) * 2019-05-29 2021-08-27 大连理工大学 一种护岸防浪墙断面的参数测定优化方法
CN112182995B (zh) * 2020-10-27 2021-06-04 中国海洋大学 一种粘性势流理论分析方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103542999A (zh) * 2013-09-30 2014-01-29 上海交通大学 核电站防浪堤抗波浪冲击性能数值检测方法
CN105200957A (zh) * 2015-09-17 2015-12-30 河海大学 带圆弧板的π型桩基透空式防波堤及其设计方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5865806B2 (ja) * 2012-09-05 2016-02-17 株式会社東芝 半導体装置の製造方法及び半導体製造装置

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103542999A (zh) * 2013-09-30 2014-01-29 上海交通大学 核电站防浪堤抗波浪冲击性能数值检测方法
CN105200957A (zh) * 2015-09-17 2015-12-30 河海大学 带圆弧板的π型桩基透空式防波堤及其设计方法

Also Published As

Publication number Publication date
CN107992715A (zh) 2018-05-04

Similar Documents

Publication Publication Date Title
CN108287965B (zh) 一种计算波浪与不规则结构物相互作用的方法
Yuan et al. An implicit three‐dimensional fully non‐hydrostatic model for free‐surface flows
CN107992715B (zh) 一种基于虚拟边界力法的弧板式防波堤受力计算方法
CN110008509B (zh) 一种考虑背景流场下的内孤立波作用力特性分析方法
CN107977542B (zh) 一种波浪与弧板式防波堤相互作用的计算方法
Cohan et al. Numerical simulation and analysis of the effect of rain and surface property on wind-turbine airfoil performance
CN103235878B (zh) 一种柔性网衣对波浪传播影响的模拟方法
Wang et al. ISPH simulation of scour behind seawall due to continuous tsunami overflow
Zhang et al. Three-dimensional simulation of horseshoe vortex and local scour around a vertical cylinder using an unstructured finite-volume technique
Uddin et al. Application of volume of fluid (VOF) method for prediction of wave generated by flow around cambered hydrofoil
CN108170993B (zh) 一种弧板式防波堤结构受力的计算方法
CN105631101A (zh) 一种地震作用下高耸升船机塔柱结构动态分布系数分析方法
CN107256312B (zh) 一种基于潮流环境下海湾纳潮变化量计算方法
Xu et al. Flow pattern and anti-silt measures of straight-edge forebay in large pump stations.
Chen et al. The effect of vortex induced vibrating cylinders on airfoil aerodynamics
Chen et al. On rotational flows with discontinuous vorticity beneath steady water waves near stagnation
Kim et al. Direct numerical study of air layer drag reduction phenomenon over a backward-facing step
Jifu et al. Influences of the waterway project in the Yangtze River Estuary on sediment transport
Liu et al. Numerical study on a coupled-pitching flexible hydrofoil under the semi-passive mode
Zou Coastal sediment transport simulation by smoothed particle hydrodynamics
Li et al. Numerical simulation of hydrodynamic characteristics on an arc crown wall using volume of fluid method based on BFC
Huang et al. Characteristics of the flow around a circular OWC-type wave energy converter supported by a bottom-sitting C-shaped structure
Masuda et al. Use of the Mps method to estimate the energy conversion efficiency of the owc-wec (first report)
Amiri et al. A review of physical and numerical modeling techniques for horizontal-axis wind turbine wakes
Xiangju et al. Progress in numerical simulation of high entrained air-water two-phase flow

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: Zhang Zhenhua

Inventor after: Wang Qing

Inventor after: Jiang Rui

Inventor after: Li Bo

Inventor after: Qi Hanwen

Inventor after: Zhang Xiaomin

Inventor before: Li Xueyan

Inventor before: Zhang Zhenhua

Inventor before: Jiang Rui

Inventor before: Li Bo

Inventor before: Qi Hanwen

Inventor before: Zhang Xiaomin

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 force calculation method of arc plate Breakwater Based on virtual boundary force method

Effective date of registration: 20211216

Granted publication date: 20210323

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: 20210323

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
TR01 Transfer of patent right

Effective date of registration: 20220822

Address after: No. 10-4, Ping'an Street, Zhifu District, Yantai City, Shandong Province, 264000

Patentee after: Yantai Junfeng Software Technology Co.,Ltd.

Address before: 264000 No. 186 Hongqi Middle Road, Zhifu District, Shandong, Yantai

Patentee before: LUDONG University

TR01 Transfer of patent right