CN115495919B - 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 - Google Patents
一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 Download PDFInfo
- Publication number
- CN115495919B CN115495919B CN202211205938.8A CN202211205938A CN115495919B CN 115495919 B CN115495919 B CN 115495919B CN 202211205938 A CN202211205938 A CN 202211205938A CN 115495919 B CN115495919 B CN 115495919B
- Authority
- CN
- China
- Prior art keywords
- grid
- time step
- boundary
- distribution function
- incident 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 50
- 238000005315 distribution function Methods 0.000 claims abstract description 59
- 238000004364 calculation method Methods 0.000 claims abstract description 21
- 230000001133 acceleration Effects 0.000 claims description 5
- 239000006185 dispersion Substances 0.000 claims description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000013178 mathematical model Methods 0.000 description 2
- 239000002131 composite material Substances 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
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
本申请公开了一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法,涉及海洋科学技术领域,该方法对时域缓坡方程的空间计算域进行网格划分后,在每个时间步根据入射波边界的速度势更新入射波边界处的各个网格的分布函数,继而根据各个网格的分布函数求解各个网格的波面升高后,可以更新入射波边界以外的各个网格的分布函数,然后更新反射边界处各个网格的分布函数进入下一个时间步的迭代求解直至最终输出需要的波浪参数。本申请引入格子玻尔兹曼方法,提出了一种全新的高效求解方法,边界处理方便,易于实现,提高了数值求解的效率,为实现基于时域缓坡方程的大范围波浪数值演化打下基础。
Description
技术领域
本申请涉及海洋科学技术领域,尤其是一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法。
背景技术
在近岸工程、海洋工程、近岛礁工程和海洋科学等技术领域中经常同时涉及波浪折射和绕射,这时使用单纯折射或绕射的数学模型均不适用,解决这些问题需要采用综合考虑折射和绕射的波浪数学模型。时域缓冲方程能很好地描述缓慢变化地形上的复合波的绕射、折射、浅水变形、波浪破碎、海底摩擦等过程,且更能用于描述有限带宽不规则波的传播,成为了目前最常用的波浪模型之一,具有重要的海洋工程应用价值。
目前常用的时域缓坡方程的数值求解方法主要包括有限体积法、边界元法、有限元法等,但是这些求解方法普遍都存在算法实现难度大、边界不易处理、计算量大而耗时较长等问题,影响了时域缓坡方程的推广应用。
发明内容
本申请人针对上述问题及技术需求,提出了一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法,本申请的技术方案如下:
一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法,该方法包括:
对时域缓坡方程的空间计算域进行网格划分;
其进一步的技术方案为,对时域缓坡方程的空间计算域进行网格划分的方法包括:
对时域缓坡方程的空间计算域在x轴方向上划分为M列网格,对时域缓坡方程的空间计算域在y轴方向上划分为N行网格,其中,x轴方向平行于入射波边界,y轴方向垂直于x轴方向,i∈[1,M],j∈[1,N];
则入射波边界处的网格包括j=1时的各个网格,入射波边界相邻网格包括j=2时的各个网格,入射波边界以外的各个网格包括j≠1时的各个网格。
更新得到入射波边界以外的任意网格celli,j在时间步tm+1的平衡函数的方法包括:
其进一步的技术方案为,更新得到入射波边界处的各个网格在时间步tm+1的分布函数,包括将入射波边界处的任意网格celli,1在时间步tm+1的任意第k个分布函数更新为:
其中,是网格celli,1对应的入射波边界相邻网格celli,2在时间步tm的第r个分布函数, 表示网格celli,2在时间步tm的波面升高,表示网格celli,2在时间步tm的速度势,r=0的分布函数位于网格的中心,r≠0的分布函数分布在网格的其他区域。
其进一步的技术方案为,r∈[0,4],r=0的分布函数位于网格的中心,r=1和r=3分别位于网格的中心沿着x轴方向的两侧且沿着y轴相对于网格中心对称,r=2和r=4分别位于网格的中心沿着y轴方向的两侧且沿着x轴相对于网格中心对称,r=1位于靠近第1列网格的一侧、r=3位于靠近第M列网格的一侧,r=2位于靠近第1行网格的一侧、r=4位于靠近第N行网格的一侧;
更新反射边界处各个网格的分布函数包括:
本申请的有益技术效果是:
本申请公开了一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法,该方法引入格子玻尔兹曼方法用于完成时域缓坡方程的数值求解,边界处理方便,易于实现,相比于现有的各种数值求解方法来说,降低了实现难度,提高了求解的效率,提出了一种全新的高效求解方法,有利于时域缓坡方程的推广应用。
另外得益于格子玻尔兹曼的特性,该方法适用于GPU并行处理,可以进一步加快计算速率,而且格子玻尔兹曼方法在处理类似的扩散方程中有大量的应用,因此使用LBM方法求解时域缓坡方程能为实现基于缓坡方程的大范围波浪数值演化打下基础。
附图说明
图1是本申请一个实施例中的数值求解方法的求解数据流示意图。
图2是本申请一个实施例中的对时域缓坡方程的空间计算域的网格划分示意图以及单个网格内的5个分布函数的示意图。
具体实施方式
下面结合附图对本申请的具体实施方式做进一步说明。
本申请公开了一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法,常见的时域缓坡方程表示为其中,c表示相速度,cg表示群速度,k表示波数,φ表示势函数。进一步的,该时域缓坡方程还有其他的诸如抛物型缓坡方程、双曲型缓坡方程等变式,本申请不一一列举。请参考图1所示的计算流程图,该方法包括如下步骤:
步骤110,对时域缓坡方程的空间计算域进行网格划分。
在一个实施例中,请参考图2,以x轴方向平行于入射波边界,y轴方向垂直于x轴方向,且y=0作为入射波边界,对时域缓坡方程的空间计算域在x轴方向上划分为M列网格,对时域缓坡方程的空间计算域在y轴方向上划分为N行网格。
以网格celli,j表示任意第j行、第i列的网格,i∈[1,M],j∈[1,N]。则可以确定,入射波边界处的网格包括j=1时的各个网格celli,1,入射波边界以外的各个网格包括j≠1时的各个网格celli,j。
本申请除了对空间计算域进行网格划分之外,实际还对计算时段进行时间步的划分,也即将计算时段按照步长dt划分为若干个时间步,任意第m个时间步表示为时间步tm,m为起始值为0的参数,本申请的方法按照各个时间步依次迭代处理。
在对空间计算域进行网格划分后,就可以基于空间计算域的地形数据确定任意网格celli,j的水深hi,j。
步骤120,对网格进行初始化参数设置,也即初始化任意一个网格celli,j在时间步t0的参数,包括:初始化网格celli,j在时间步t0的波面升高初始化网格celli,j在时间步t0的任意第r个分布函数任意第r个平衡函数任意第r个源项函数R为参数。在一个实施例中,每个网格包括位于网格中心处的分布函数以及位于网格内其他位置处的多个分布函数,一种实现方式是,请参考图2,R=4,每个网格都有5个分布函数,r=0的分布函数位于网格的中心,r≠0的分布函数分布在网格的其他区域,其中,r=1和r=3分别位于网格的中心沿着x轴方向的两侧且沿着y轴相对于网格中心对称,r=2和r=4分别位于网格的中心沿着y轴方向的两侧且沿着x轴相对于网格中心对称,r=1位于靠近第1列网格的一侧、r=3位于靠近第M列网格的一侧,r=2位于靠近第1行网格的一侧、r=4位于靠近第N行网格的一侧。除此之外,还初始化任意网格celli,j在时间步t0的速度势
步骤130,更新得到入射波边界处的各个网格在时间步tm+1的分布函数。如上所述,入射波边界处的各个网格即为j=1的各个网格celli,1,在更新得到各个网格celli,1在时间步tm+1的任意第r个分布函数时,需要使用到入射波边界在时间步tm+1的速度势以及入射波边界相邻网格在时间步tm的分布函数。
由于入射波边界处的各个网格即为j=1的各个网格celli,1,因此入射波边界相邻网格包括j=2时的各个网格celli,2。结合以及网格celli,1对应的入射波边界相邻网格celli,2的网格参数可以更新得到为:
步骤140,基于各个网格的分布函数计算得到对应网格celli,j在时间步tm+1的波面升高包括:基于入射波边界处的各个j=1的网格更新后的时间步tm+1的确定对应的网格的波面升高以及根据入射波边界以外的各个j≠1的网格的还未更新的确定对应的网格的波面升高也即有
步骤150,更新得到入射波边界以外的各个网格在时间步tm+1的分布函数包括:基于任意网格celli,j在时间步tm+1的波面升高更新得到入射波边界以外的任意网格celli,j在时间步tm+1的平衡函数以及源项函数并计算得到入射波边界以外的各个网格在时间步tm+1的分布函数分如下几部分分别介绍:
首先根据任意网格celli,j在时间步tm+1的波面升高计算得到网格celli,j在时间步tm+1的速度势的源项然后可以更新得到任意网格celli,j在时间步tm+1的速度势dt是相邻两个时间步之间的步长。
(2)如上所述,基于地形数据可以确定任意网格celli,j的水深hi,j,则结合网格celli,j的水深hi,j以及网格celli,j在时间步tm+1的波高根据色散关系计算得到网格celli,j在时间步tm+1的波数相速度和群速度
步骤160,更新反射边界处各个网格的分布函数。反射边界包括第M列的网格构成的边界、第1列的网格构成的边界、第N行的网格构成的边界,则反射边界处的各个网格包括i=M的网格、i=1的网格以及j=N的网格。在更新反射边界处的各个网格的分布函数时,按照分布函数的对称关系,对于i=M的各个网格cellM,j,令r=3的分布函数与网格cellM,j的r=1的分布函数相等。对于i=1的各个网格cell1,j,令网格cell1,j的r=1的分布函数与r=3的分布函数相等。对于j=N的各个网格celli,N,令网格celli,N的r=4的分布函数与r=2的分布函数相等。可以表示为
步骤170,如果当前时间步tm+1为最后一个时间步,则输出所需的波浪参数,完成数值求解过程。如果当前时间步tm+1还不是最后一个时间步,则令m=m+1并再次执行步骤130~170执行下一个时间步的循环,直至完成最后一个时间步的计算输出所需的波浪参数。
以上所述的仅是本申请的优选实施方式,本申请不限于以上实施例。可以理解,本领域技术人员在不脱离本申请的精神和构思的前提下直接导出或联想到的其他改进和变化,均应认为包含在本申请的保护范围之内。
Claims (9)
1.一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法,其特征在于,所述方法包括:
对时域缓坡方程的空间计算域进行网格划分;
更新反射边界处各个网格的分布函数,令m=m+1并再次执行所述根据入射波边界在时间步tm+1的速度势以及入射波边界相邻网格在时间步tm的分布函数的步骤,直到时间步tm+1为最后一个时间步时,输出波浪参数;
2.根据权利要求1所述的数值求解方法,其特征在于,所述对时域缓坡方程的空间计算域进行网格划分的方法包括:
对时域缓坡方程的空间计算域在x轴方向上划分为M列网格,对所述时域缓坡方程的空间计算域在y轴方向上划分为N行网格,其中,x轴方向平行于入射波边界,y轴方向垂直于x轴方向,i∈[1,M],j∈[1,N];
则入射波边界处的网格包括j=1时的各个网格,入射波边界相邻网格包括j=2时的各个网格,入射波边界以外的各个网格包括j≠1时的各个网格。
9.根据权利要求2所述的数值求解方法,其特征在于,r∈[0,4],r=0的分布函数位于网格的中心,r=1和r=3分别位于网格的中心沿着x轴方向的两侧且沿着y轴相对于网格中心对称,r=2和r=4分别位于网格的中心沿着y轴方向的两侧且沿着x轴相对于网格中心对称,r=1位于靠近第1列网格的一侧、r=3位于靠近第M列网格的一侧,r=2位于靠近第1行网格的一侧、r=4位于靠近第N行网格的一侧;
所述更新反射边界处各个网格的分布函数包括:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211205938.8A CN115495919B (zh) | 2022-09-30 | 2022-09-30 | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211205938.8A CN115495919B (zh) | 2022-09-30 | 2022-09-30 | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN115495919A CN115495919A (zh) | 2022-12-20 |
CN115495919B true CN115495919B (zh) | 2023-05-26 |
Family
ID=84472923
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202211205938.8A Active CN115495919B (zh) | 2022-09-30 | 2022-09-30 | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115495919B (zh) |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107515987A (zh) * | 2017-08-25 | 2017-12-26 | 中国地质大学(北京) | 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 |
CN109325309A (zh) * | 2018-10-23 | 2019-02-12 | 哈尔滨工程大学 | 船舶大幅横摇运动的三维数值模拟方法 |
CN109446634A (zh) * | 2018-10-23 | 2019-03-08 | 哈尔滨工程大学 | 基于泰勒展开边界元方法的船舶运动预报方法 |
CN110275733A (zh) * | 2019-06-27 | 2019-09-24 | 上海交通大学 | 基于有限体积法求解声子玻尔兹曼方程的gpu并行加速方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US11010509B2 (en) * | 2018-05-23 | 2021-05-18 | Nvidia Corporation | Systems and methods for computer simulation of detailed waves for large-scale water simulation |
-
2022
- 2022-09-30 CN CN202211205938.8A patent/CN115495919B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107515987A (zh) * | 2017-08-25 | 2017-12-26 | 中国地质大学(北京) | 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 |
CN109325309A (zh) * | 2018-10-23 | 2019-02-12 | 哈尔滨工程大学 | 船舶大幅横摇运动的三维数值模拟方法 |
CN109446634A (zh) * | 2018-10-23 | 2019-03-08 | 哈尔滨工程大学 | 基于泰勒展开边界元方法的船舶运动预报方法 |
CN110275733A (zh) * | 2019-06-27 | 2019-09-24 | 上海交通大学 | 基于有限体积法求解声子玻尔兹曼方程的gpu并行加速方法 |
Non-Patent Citations (2)
Title |
---|
基于SVS-603传感器的波浪反演技术及海上比测数据分析;王志等;山东科学;第34卷(第3期);1-8 * |
格子玻尔兹曼方法中的特征无反射边界条件研究;邵卫东;李军;;工程热物理学报(02);277-284 * |
Also Published As
Publication number | Publication date |
---|---|
CN115495919A (zh) | 2022-12-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106842306B (zh) | 一种全局优化的交错网格有限差分正演模拟方法和装置 | |
CN110058307B (zh) | 一种基于快速拟牛顿法的全波形反演方法 | |
CN113259034B (zh) | 一种并行的耦合海洋声学预报系统及运行方法 | |
CN111858465B (zh) | 大规模矩阵qr分解并行计算系统 | |
CN116774292B (zh) | 一种地震波走时确定方法、系统、电子设备及存储介质 | |
CN113158492B (zh) | 一种时变电磁场的全隐式双时间步计算方法 | |
Kaliukh et al. | Constructing the adaptive algorithms for solving multi-wave problems | |
CN107290793A (zh) | 一种基于加权多策略蛙跳算法的超高密度电法并行反演方法 | |
CN115495919B (zh) | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 | |
CN115034360A (zh) | 三维卷积神经网络卷积层的处理方法和处理装置 | |
CN109343006B (zh) | 基于增广拉格朗日遗传算法的nflm信号优化方法及装置 | |
CN113792440A (zh) | 计算浮式结构物在非稳态载荷作用下结构力响应的方法 | |
CN105844018A (zh) | 一种大型反射面天线反射体俯仰关键模态的选取方法 | |
CN104346488B (zh) | 电大复杂外形金属目标混合建模及电磁散射快速仿真方法 | |
CN110826197B (zh) | 一种基于改进Cholesky分解闭合解的风速场模拟方法 | |
CN106662665A (zh) | 用于更快速的交错网格处理的重新排序的插值和卷积 | |
CN108983290A (zh) | 一种三维横向各向同性介质中旅行时确定方法及系统 | |
CN109298451A (zh) | 一种改进偏斜度的自动拾取s波震相方法 | |
CN114154111A (zh) | 一种频率域连分式展开的位场数据向下延拓方法 | |
CN110649912B (zh) | 空间滤波器的建模方法 | |
CN113866827A (zh) | 一种解释性速度建模地震成像方法、系统、介质和设备 | |
CN115563838B (zh) | 一种基于comsol的地震超材料能带结构和频域分析方法 | |
CN104391824A (zh) | 一种基于a=lr三角分解法快速求取电力系统节点阻抗矩阵的方法 | |
CN112578431B (zh) | 一种限存状态全波形反演波场最优化存储方法及系统 | |
CN101866381A (zh) | 一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法 |
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 | ||
GR01 | Patent grant | ||
GR01 | Patent grant |