CN115495919B - 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 - Google Patents

一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 Download PDF

Info

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
Application number
CN202211205938.8A
Other languages
English (en)
Other versions
CN115495919A (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.)
702th Research Institute of CSIC
Original Assignee
702th Research Institute of CSIC
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 702th Research Institute of CSIC filed Critical 702th Research Institute of CSIC
Priority to CN202211205938.8A priority Critical patent/CN115495919B/zh
Publication of CN115495919A publication Critical patent/CN115495919A/zh
Application granted granted Critical
Publication of CN115495919B publication Critical patent/CN115495919B/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
    • 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

一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法
技术领域
本申请涉及海洋科学技术领域，尤其是一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法。
背景技术
在近岸工程、海洋工程、近岛礁工程和海洋科学等技术领域中经常同时涉及波浪折射和绕射，这时使用单纯折射或绕射的数学模型均不适用，解决这些问题需要采用综合考虑折射和绕射的波浪数学模型。时域缓冲方程能很好地描述缓慢变化地形上的复合波的绕射、折射、浅水变形、波浪破碎、海底摩擦等过程，且更能用于描述有限带宽不规则波的传播，成为了目前最常用的波浪模型之一，具有重要的海洋工程应用价值。
目前常用的时域缓坡方程的数值求解方法主要包括有限体积法、边界元法、有限元法等，但是这些求解方法普遍都存在算法实现难度大、边界不易处理、计算量大而耗时较长等问题，影响了时域缓坡方程的推广应用。
发明内容
本申请人针对上述问题及技术需求，提出了一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法，本申请的技术方案如下：
一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法，该方法包括：
对时域缓坡方程的空间计算域进行网格划分；
对于任意第j行、第i列的网格celli,j，初始化第0个时间步t0的波面升高
Figure GDA0004168657850000011
任意第r个分布函数/>
Figure GDA0004168657850000012
任意第r个平衡函数/>
Figure GDA0004168657850000013
任意第r个源项函数
Figure GDA0004168657850000014
R为参数；
根据入射波边界在时间步tm+1的速度势
Figure GDA0004168657850000015
以及入射波边界相邻网格在时间步tm的分布函数，更新得到入射波边界处的各个网格在时间步tm+1的分布函数，m为起始值为0的参数；
基于各个网格的分布函数计算得到对应网格在时间步tm+1的波面升高
Figure GDA0004168657850000016
基于任意网格celli,j在时间步tm+1的波面升高
Figure GDA0004168657850000017
更新得到入射波边界以外的任意网格celli,j在时间步tm+1的平衡函数/>
Figure GDA0004168657850000018
以及源项函数/>
Figure GDA0004168657850000019
并计算得到入射波边界以外的各个网格在时间步tm+1的分布函数/>
Figure GDA0004168657850000021
更新反射边界处各个网格的分布函数，令m＝m+1并再次执行根据入射波边界在时间步tm+1的速度势
Figure GDA0004168657850000022
以及入射波边界相邻网格在时间步tm的分布函数的步骤，直到时间步tm+1为最后一个时间步时，输出波浪参数。/>
其进一步的技术方案为，对时域缓坡方程的空间计算域进行网格划分的方法包括：
对时域缓坡方程的空间计算域在x轴方向上划分为M列网格，对时域缓坡方程的空间计算域在y轴方向上划分为N行网格，其中，x轴方向平行于入射波边界，y轴方向垂直于x轴方向，i∈[1,M]，j∈[1,N]；
则入射波边界处的网格包括j＝1时的各个网格，入射波边界相邻网格包括j＝2时的各个网格，入射波边界以外的各个网格包括j≠1时的各个网格。
其进一步的技术方案为，该方法还包括：初始化任意网格celli,j在时间步t0的速度势
Figure GDA0004168657850000023
更新得到入射波边界以外的任意网格celli,j在时间步tm+1的平衡函数的方法包括：
根据任意网格celli,j在时间步tm+1的波面升高
Figure GDA0004168657850000024
结合网格celli,j在时间步tm的速度势/>
Figure GDA0004168657850000025
更新得到任意网格celli,j在时间步tm+1的速度势/>
Figure GDA0004168657850000026
根据入射波边界以外的任意网格celli,j对应的波面升高
Figure GDA0004168657850000027
和速度势/>
Figure GDA0004168657850000028
更新得到网格celli,j在时间步tm+1的平衡函数/>
Figure GDA0004168657850000029
其进一步的技术方案为，更新得到入射波边界处的各个网格在时间步tm+1的分布函数，包括将入射波边界处的任意网格celli,1在时间步tm+1的任意第k个分布函数更新为：
Figure GDA00041686578500000210
其中，
Figure GDA00041686578500000211
是网格celli,1对应的入射波边界相邻网格celli,2在时间步tm的第r个分布函数，/>
Figure GDA00041686578500000212
Figure GDA00041686578500000213
表示网格celli,2在时间步tm的波面升高，/>
Figure GDA00041686578500000214
表示网格celli,2在时间步tm的速度势，r＝0的分布函数位于网格的中心，r≠0的分布函数分布在网格的其他区域。
其进一步的技术方案为，计算得到对应网格在时间步tm+1的波面升高
Figure GDA0004168657850000031
包括确定/>
Figure GDA0004168657850000032
其进一步的技术方案为，更新得到任意网格celli,j在时间步tm+1的速度势
Figure GDA0004168657850000033
的方法包括：
根据任意网格celli,j在时间步tm+1的波面升高
Figure GDA0004168657850000034
计算得到网格celli,j在时间步tm+1的速度势的源项/>
Figure GDA0004168657850000035
g是重力加速度；/>
更新得到网格celli,j在时间步tm+1的速度势
Figure GDA0004168657850000036
dt是相邻两个时间步之间的步长。
其进一步的技术方案为，更新得到网格celli,j在时间步tm+1的平衡函数
Figure GDA0004168657850000037
的方法，包括确定入射波边界以外的任意网格celli,j的/>
Figure GDA0004168657850000038
为
Figure GDA0004168657850000039
的分布函数位于网格的中心，r≠0的分布函数分布在网格的其他区域。
其进一步的技术方案为，更新得到入射波边界以外的任意网格celli,j在时间步tm+1的源项函数
Figure GDA00041686578500000310
的方法包括：
根据任意网格celli,j在历史多个时间步的波面升高的使用上跨零算法计算得到网格celli,j在时间步tm+1的波高
Figure GDA00041686578500000311
基于地形数据确定任意网格celli,j的水深hi,j，并结合网格celli,j在时间步tm+1的波高
Figure GDA00041686578500000312
根据色散关系计算得到网格celli,j在时间步tm+1的波数/>
Figure GDA00041686578500000313
相速度/>
Figure GDA00041686578500000314
和群速度
Figure GDA00041686578500000315
更新入射波边界以外的任意网格celli,j在时间步tm+1的源项函数
Figure GDA00041686578500000316
其中，ωIN是入射波圆频率，g是重力加速度，
Figure GDA00041686578500000317
其进一步的技术方案为，计算得到入射波边界以外的各个网格在时间步tm+1的分布函数的方法，包括确定入射波边界以外的任意网格celli,j在时间步tm+1的任意第r个分布函数
Figure GDA0004168657850000041
为：
Figure GDA0004168657850000042
其中，
Figure GDA0004168657850000043
是网格celli,j在时间步tm的第r个分布函数，τ为参数，dt是相邻两个时间步之间的步长。
其进一步的技术方案为，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行网格的一侧；
更新反射边界处各个网格的分布函数包括：
将位于第M列的各个网格cellM,j在时间步tm+1的r＝3的分布函数
Figure GDA0004168657850000044
更新为与网格cellM,j在时间步tm+1的r＝1的分布函数/>
Figure GDA0004168657850000045
相同；
以及，将第1列的各个网格cell1,j在时间步tm+1的r＝1的分布函数
Figure GDA0004168657850000046
更新为与网格cell1,j在时间步tm+1的r＝3的分布函数/>
Figure GDA0004168657850000047
相同；
以及，将第N行的各个网格celli,N在时间步tm+1的r＝4的分布函数
Figure GDA0004168657850000048
更新为与网格celli,N在时间步tm+1的r＝2的分布函数/>
Figure GDA0004168657850000049
相同。
本申请的有益技术效果是：
本申请公开了一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法，该方法引入格子玻尔兹曼方法用于完成时域缓坡方程的数值求解，边界处理方便，易于实现，相比于现有的各种数值求解方法来说，降低了实现难度，提高了求解的效率，提出了一种全新的高效求解方法，有利于时域缓坡方程的推广应用。
另外得益于格子玻尔兹曼的特性，该方法适用于GPU并行处理，可以进一步加快计算速率，而且格子玻尔兹曼方法在处理类似的扩散方程中有大量的应用，因此使用LBM方法求解时域缓坡方程能为实现基于缓坡方程的大范围波浪数值演化打下基础。
附图说明
图1是本申请一个实施例中的数值求解方法的求解数据流示意图。
图2是本申请一个实施例中的对时域缓坡方程的空间计算域的网格划分示意图以及单个网格内的5个分布函数的示意图。
具体实施方式
下面结合附图对本申请的具体实施方式做进一步说明。
本申请公开了一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法，常见的时域缓坡方程表示为
Figure GDA0004168657850000051
其中，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的波面升高
Figure GDA0004168657850000052
初始化网格celli,j在时间步t0的任意第r个分布函数/>
Figure GDA0004168657850000053
任意第r个平衡函数/>
Figure GDA0004168657850000054
任意第r个源项函数/>
Figure GDA0004168657850000055
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的速度势/>
Figure GDA0004168657850000061
步骤130，更新得到入射波边界处的各个网格在时间步tm+1的分布函数。如上所述，入射波边界处的各个网格即为j＝1的各个网格celli,1，在更新得到各个网格celli,1在时间步tm+1的任意第r个分布函数
Figure GDA0004168657850000062
时，需要使用到入射波边界在时间步tm+1的速度势/>
Figure GDA0004168657850000063
以及入射波边界相邻网格在时间步tm的分布函数。
由于入射波边界处的各个网格即为j＝1的各个网格celli,1，因此入射波边界相邻网格包括j＝2时的各个网格celli,2。结合
Figure GDA0004168657850000064
以及网格celli,1对应的入射波边界相邻网格celli,2的网格参数可以更新得到/>
Figure GDA0004168657850000065
为：
Figure GDA0004168657850000066
其中，
Figure GDA0004168657850000067
是网格celli,1对应的入射波边界相邻网格celli,2在时间步tm的第r个分布函数。另外：
Figure GDA0004168657850000068
Figure GDA0004168657850000069
在上式中，
Figure GDA00041686578500000610
表示网格celli,2在时间步tm的波面升高，/>
Figure GDA00041686578500000611
表示网格celli,2在时间步tm的速度势。入射波边界在时间步tm+1的速度势/>
Figure GDA00041686578500000612
其中，ωIN是入射波圆频率，HIN是入射波高，g是重力加速度。/>
步骤140，基于各个网格的分布函数计算得到对应网格celli,j在时间步tm+1的波面升高
Figure GDA00041686578500000613
包括：基于入射波边界处的各个j＝1的网格更新后的时间步tm+1的/>
Figure GDA00041686578500000614
确定对应的网格的波面升高/>
Figure GDA00041686578500000615
以及根据入射波边界以外的各个j≠1的网格的还未更新的/>
Figure GDA00041686578500000616
确定对应的网格的波面升高/>
Figure GDA00041686578500000617
也即有
Figure GDA0004168657850000071
步骤150，更新得到入射波边界以外的各个网格在时间步tm+1的分布函数
Figure GDA0004168657850000072
包括：基于任意网格celli,j在时间步tm+1的波面升高/>
Figure GDA0004168657850000073
更新得到入射波边界以外的任意网格celli,j在时间步tm+1的平衡函数/>
Figure GDA0004168657850000074
以及源项函数/>
Figure GDA0004168657850000075
并计算得到入射波边界以外的各个网格在时间步tm+1的分布函数/>
Figure GDA0004168657850000076
分如下几部分分别介绍：
1、计算网格celli,j在时间步tm+1的平衡函数
Figure GDA0004168657850000077
包括：
(1)根据任意网格celli,j在时间步tm+1的波面升高
Figure GDA0004168657850000078
结合网格celli,j在时间步tm的速度势/>
Figure GDA0004168657850000079
更新得到任意网格celli,j在时间步tm+1的速度势/>
Figure GDA00041686578500000710
包括：
首先根据任意网格celli,j在时间步tm+1的波面升高
Figure GDA00041686578500000711
计算得到网格celli,j在时间步tm+1的速度势的源项/>
Figure GDA00041686578500000712
然后可以更新得到任意网格celli,j在时间步tm+1的速度势/>
Figure GDA00041686578500000713
dt是相邻两个时间步之间的步长。
(2)根据入射波边界以外的任意网格celli,j对应的波面升高
Figure GDA00041686578500000714
和速度势/>
Figure GDA00041686578500000715
更新得到网格celli,j在时间步tm+1的平衡函数/>
Figure GDA00041686578500000716
为：
Figure GDA00041686578500000717
2、计算网格celli,j在时间步tm+1的源项函数
Figure GDA00041686578500000718
包括：
(1)根据任意网格celli,j在历史多个时间步的波面升高的使用上跨零算法计算得到网格celli,j在时间步tm+1的波高
Figure GDA00041686578500000719
(2)如上所述，基于地形数据可以确定任意网格celli,j的水深hi,j，则结合网格celli,j的水深hi,j以及网格celli,j在时间步tm+1的波高
Figure GDA00041686578500000720
根据色散关系计算得到网格celli,j在时间步tm+1的波数/>
Figure GDA00041686578500000721
相速度/>
Figure GDA00041686578500000722
和群速度/>
Figure GDA00041686578500000723
/>
(3)基于网格celli,j在时间步tm+1的波数
Figure GDA00041686578500000724
相速度/>
Figure GDA00041686578500000725
和群速度/>
Figure GDA00041686578500000726
更新入射波边界以外的任意网格celli,j在时间步tm+1的源项函数/>
Figure GDA00041686578500000727
为：
Figure GDA0004168657850000081
其中，
Figure GDA0004168657850000082
3、对于入射波边界以外的任意网格celli,j，根据网格celli,j在时间步tm+1的平衡函数
Figure GDA0004168657850000083
以及源项函数/>
Figure GDA0004168657850000084
计算得到该网格celli,j在时间步tm+1的分布函数
Figure GDA0004168657850000085
计算方法为：
Figure GDA0004168657850000086
其中，
Figure GDA0004168657850000087
是网格celli,j在时间步tm的第r个分布函数，τ为参数，dt是相邻两个时间步之间的步长。
步骤160，更新反射边界处各个网格的分布函数。反射边界包括第M列的网格构成的边界、第1列的网格构成的边界、第N行的网格构成的边界，则反射边界处的各个网格包括i＝M的网格、i＝1的网格以及j＝N的网格。在更新反射边界处的各个网格的分布函数时，按照分布函数的对称关系，对于i＝M的各个网格cellM,j，令r＝3的分布函数
Figure GDA0004168657850000088
与网格cellM,j的r＝1的分布函数/>
Figure GDA0004168657850000089
相等。对于i＝1的各个网格cell1,j，令网格cell1,j的r＝1的分布函数/>
Figure GDA00041686578500000810
与r＝3的分布函数/>
Figure GDA00041686578500000811
相等。对于j＝N的各个网格celli,N，令网格celli,N的r＝4的分布函数/>
Figure GDA00041686578500000812
与r＝2的分布函数/>
Figure GDA00041686578500000813
相等。可以表示为/>
Figure GDA00041686578500000814
步骤170，如果当前时间步tm+1为最后一个时间步，则输出所需的波浪参数，完成数值求解过程。如果当前时间步tm+1还不是最后一个时间步，则令m＝m+1并再次执行步骤130～170执行下一个时间步的循环，直至完成最后一个时间步的计算输出所需的波浪参数。
以上所述的仅是本申请的优选实施方式，本申请不限于以上实施例。可以理解，本领域技术人员在不脱离本申请的精神和构思的前提下直接导出或联想到的其他改进和变化，均应认为包含在本申请的保护范围之内。

Claims (9)

1.一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法，其特征在于，所述方法包括：
对时域缓坡方程的空间计算域进行网格划分；
对于任意第j行、第i列的网格celli,j，初始化第0个时间步t0的波面升高
Figure FDA0004168657840000011
任意第r个分布函数/>
Figure FDA0004168657840000012
任意第r个平衡函数/>
Figure FDA0004168657840000013
任意第r个源项函数
Figure FDA0004168657840000014
r∈[0,R]，R为参数；
根据入射波边界在时间步tm+1的速度势
Figure FDA0004168657840000015
以及入射波边界相邻网格在时间步tm的分布函数，更新得到入射波边界处的各个网格在时间步tm+1的分布函数，m为起始值为0的参数；
基于各个网格的分布函数计算得到对应网格在时间步tm+1的波面升高
Figure FDA0004168657840000016
基于任意网格celli,j在时间步tm+1的波面升高
Figure FDA0004168657840000017
更新得到入射波边界以外的任意网格celli,j在时间步tm+1的平衡函数/>
Figure FDA0004168657840000018
以及源项函数/>
Figure FDA0004168657840000019
并计算得到入射波边界以外的各个网格在时间步tm+1的分布函数/>
Figure FDA00041686578400000110
更新反射边界处各个网格的分布函数，令m＝m+1并再次执行所述根据入射波边界在时间步tm+1的速度势
Figure FDA00041686578400000111
以及入射波边界相邻网格在时间步tm的分布函数的步骤，直到时间步tm+1为最后一个时间步时，输出波浪参数；
其中，更新得到入射波边界以外的任意网格celli,j在时间步tm+1的源项函数
Figure FDA00041686578400000112
的方法包括：
根据任意网格celli,j在历史多个时间步的波面升高的使用上跨零算法计算得到网格celli,j在时间步tm+1的波高
Figure FDA00041686578400000113
基于地形数据确定任意网格celli,j的水深hi,j，并结合网格celli,j在时间步tm+1的波高
Figure FDA00041686578400000114
根据色散关系计算得到网格celli,j在时间步tm+1的波数/>
Figure FDA00041686578400000115
相速度/>
Figure FDA00041686578400000116
和群速度
Figure FDA00041686578400000117
更新入射波边界以外的任意网格celli,j在时间步tm+1的源项函数
Figure FDA0004168657840000021
其中，ωIN是入射波圆频率，g是重力加速度，
Figure FDA0004168657840000022
2.根据权利要求1所述的数值求解方法，其特征在于，所述对时域缓坡方程的空间计算域进行网格划分的方法包括：
对时域缓坡方程的空间计算域在x轴方向上划分为M列网格，对所述时域缓坡方程的空间计算域在y轴方向上划分为N行网格，其中，x轴方向平行于入射波边界，y轴方向垂直于x轴方向，i∈[1,M]，j∈[1,N]；
则入射波边界处的网格包括j＝1时的各个网格，入射波边界相邻网格包括j＝2时的各个网格，入射波边界以外的各个网格包括j≠1时的各个网格。
3.根据权利要求2所述的数值求解方法，其特征在于，所述方法还包括：初始化任意网格celli,j在时间步t0的速度势
Figure FDA0004168657840000023
更新得到入射波边界以外的任意网格celli,j在时间步tm+1的平衡函数的方法包括：
根据任意网格celli,j在时间步tm+1的波面升高
Figure FDA0004168657840000024
结合网格celli,j在时间步tm的速度势/>
Figure FDA0004168657840000025
更新得到任意网格celli,j在时间步tm+1的速度势/>
Figure FDA0004168657840000026
根据入射波边界以外的任意网格celli,j对应的波面升高
Figure FDA0004168657840000027
和速度势/>
Figure FDA0004168657840000028
更新得到网格celli,j在时间步tm+1的平衡函数/>
Figure FDA0004168657840000029
4.根据权利要求3所述的数值求解方法，其特征在于，所述更新得到入射波边界处的各个网格在时间步tm+1的分布函数，包括将入射波边界处的任意网格celli,1在时间步tm+1的任意第k个分布函数更新为：
Figure FDA00041686578400000210
其中，
Figure FDA00041686578400000211
是网格celli,1对应的入射波边界相邻网格celli,2在时间步tm的第r个分布函数，/>
Figure FDA00041686578400000212
Figure FDA00041686578400000213
表示网格celli,2在时间步tm的波面升高，/>
Figure FDA00041686578400000214
表示网格celli,2在时间步tm的速度势，r＝0的分布函数位于网格的中心，r≠0的分布函数分布在网格的其他区域。
5.根据权利要求4所述的数值求解方法，其特征在于，所述计算得到对应网格在时间步tm+1的波面升高
Figure FDA0004168657840000031
包括确定/>
Figure FDA0004168657840000032
6.根据权利要求3所述的数值求解方法，其特征在于，更新得到任意网格celli,j在时间步tm+1的速度势
Figure FDA0004168657840000033
的方法包括：
根据任意网格celli,j在时间步tm+1的波面升高
Figure FDA0004168657840000034
计算得到网格celli,j在时间步tm+1的速度势的源项/>
Figure FDA0004168657840000035
g是重力加速度；
更新得到网格celli,j在时间步tm+1的速度势
Figure FDA0004168657840000036
dt是相邻两个时间步之间的步长。/>
7.根据权利要求3所述的数值求解方法，其特征在于，更新得到网格celli,j在时间步tm+1的平衡函数
Figure FDA0004168657840000037
的方法，包括确定入射波边界以外的任意网格celli,j的/>
Figure FDA0004168657840000038
为
Figure FDA0004168657840000039
r＝0的分布函数位于网格的中心，r≠0的分布函数分布在网格的其他区域。
8.根据权利要求1所述的数值求解方法，其特征在于，计算得到入射波边界以外的各个网格在时间步tm+1的分布函数的方法，包括确定入射波边界以外的任意网格celli,j在时间步tm+1的任意第r个分布函数
Figure FDA00041686578400000310
为：
Figure FDA00041686578400000311
其中，
Figure FDA00041686578400000312
是网格celli,j在时间步tm的第r个分布函数，τ为参数，dt是相邻两个时间步之间的步长。
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行网格的一侧；
所述更新反射边界处各个网格的分布函数包括：
将位于第M列的各个网格cellM,j在时间步tm+1的r＝3的分布函数
Figure FDA0004168657840000041
更新为与网格cellM,j在时间步tm+1的r＝1的分布函数/>
Figure FDA0004168657840000042
相同；
以及，将第1列的各个网格cell1,j在时间步tm+1的r＝1的分布函数
Figure FDA0004168657840000043
更新为与网格cell1,j在时间步tm+1的r＝3的分布函数/>
Figure FDA0004168657840000044
相同；
以及，将第N行的各个网格celli,N在时间步tm+1的r＝4的分布函数
Figure FDA0004168657840000045
更新为与网格celli,N在时间步tm+1的r＝2的分布函数/>
Figure FDA0004168657840000046
相同。/>
CN202211205938.8A 2022-09-30 2022-09-30 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 Active CN115495919B (zh)

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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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

Patent Citations (4)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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
CN112925012B (zh) 地震全波形反演方法及装置
CN106842306B (zh) 一种全局优化的交错网格有限差分正演模拟方法和装置
CN110058307B (zh) 一种基于快速拟牛顿法的全波形反演方法
CN113259034B (zh) 一种并行的耦合海洋声学预报系统及运行方法
CN116774292B (zh) 一种地震波走时确定方法、系统、电子设备及存储介质
CN101937425B (zh) 基于gpu众核平台的矩阵并行转置方法
CN113158492B (zh) 一种时变电磁场的全隐式双时间步计算方法
CN109657794B (zh) 一种基于指令队列的分布式深度神经网络性能建模方法
CN107290793A (zh) 一种基于加权多策略蛙跳算法的超高密度电法并行反演方法
CN115495919B (zh) 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法
CN106662665B (zh) 用于更快速的交错网格处理的重新排序的插值和卷积
CN110826197B (zh) 一种基于改进Cholesky分解闭合解的风速场模拟方法
US6301192B1 (en) Method for generating 2 and 3-dimensional fluid meshes for structural/acoustic finite element analysis in infinite medium
CN115293978A (zh) 卷积运算电路和方法、图像处理设备
CN109725346A (zh) 一种时间-空间高阶vti矩形网格有限差分方法及装置
CN109298451A (zh) 一种改进偏斜度的自动拾取s波震相方法
CN115034360A (zh) 三维卷积神经网络卷积层的处理方法和处理装置
CN104849747A (zh) 一种优化气枪阵列的方法和装置
CN114154111A (zh) 一种频率域连分式展开的位场数据向下延拓方法
CN117010260A (zh) 一种裂缝性油藏自动历史拟合模型预测方法、系统及设备
CN113866827A (zh) 一种解释性速度建模地震成像方法、系统、介质和设备
CN115563838B (zh) 一种基于comsol的地震超材料能带结构和频域分析方法
CN104391824A (zh) 一种基于a＝lr三角分解法快速求取电力系统节点阻抗矩阵的方法
CN103837894A (zh) 获取剩余静校正量的方法
JP7379801B2 (ja) ノイズ除去方法

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