CN113221392B - 一种非均匀粘声波在无限域内传播模型的构建方法 - Google Patents
一种非均匀粘声波在无限域内传播模型的构建方法 Download PDFInfo
- Publication number
- CN113221392B CN113221392B CN202110103651.3A CN202110103651A CN113221392B CN 113221392 B CN113221392 B CN 113221392B CN 202110103651 A CN202110103651 A CN 202110103651A CN 113221392 B CN113221392 B CN 113221392B
- Authority
- CN
- China
- Prior art keywords
- equation
- domain
- viscous
- acoustic wave
- complex
- 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
- 238000010276 construction Methods 0.000 title claims abstract description 14
- 238000000034 method Methods 0.000 claims abstract description 75
- 238000004364 calculation method Methods 0.000 claims abstract description 27
- 238000004088 simulation Methods 0.000 claims abstract description 14
- 238000012892 rational function Methods 0.000 claims abstract description 12
- 238000003860 storage Methods 0.000 claims abstract description 10
- 238000001514 detection method Methods 0.000 claims abstract description 9
- 238000012876 topography Methods 0.000 claims abstract description 9
- 229910052500 inorganic mineral Inorganic materials 0.000 claims abstract description 7
- 239000011707 mineral Substances 0.000 claims abstract description 7
- 238000006073 displacement reaction Methods 0.000 claims description 23
- 238000005457 optimization Methods 0.000 claims description 15
- 230000001419 dependent effect Effects 0.000 claims description 14
- 230000010354 integration Effects 0.000 claims description 14
- 238000012360 testing method Methods 0.000 claims description 9
- 239000011159 matrix material Substances 0.000 claims description 6
- 238000003780 insertion Methods 0.000 claims description 5
- 230000037431 insertion Effects 0.000 claims description 5
- 238000013507 mapping Methods 0.000 claims description 5
- 238000004422 calculation algorithm Methods 0.000 claims description 4
- 238000004590 computer program Methods 0.000 claims description 4
- 239000013598 vector Substances 0.000 claims description 4
- 239000002245 particle Substances 0.000 claims description 3
- 102100029469 WD repeat and HMG-box DNA-binding protein 1 Human genes 0.000 claims description 2
- 101710097421 WD repeat and HMG-box DNA-binding protein 1 Proteins 0.000 claims description 2
- 230000004069 differentiation Effects 0.000 claims description 2
- 238000000518 rheometry Methods 0.000 claims description 2
- 235000019333 sodium laurylsulphate Nutrition 0.000 claims 3
- 238000001370 static light scattering Methods 0.000 claims 3
- 230000003595 spectral effect Effects 0.000 abstract description 9
- 230000008569 process Effects 0.000 abstract description 8
- 230000006870 function Effects 0.000 description 14
- 238000013016 damping Methods 0.000 description 11
- 238000001228 spectrum Methods 0.000 description 11
- 239000006185 dispersion Substances 0.000 description 8
- 230000009466 transformation Effects 0.000 description 7
- 238000013459 approach Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 230000008901 benefit Effects 0.000 description 5
- 238000009795 derivation Methods 0.000 description 5
- 230000001052 transient effect Effects 0.000 description 5
- 230000007774 longterm Effects 0.000 description 4
- 238000012937 correction Methods 0.000 description 3
- 230000008878 coupling Effects 0.000 description 3
- 238000010168 coupling process Methods 0.000 description 3
- 238000005859 coupling reaction Methods 0.000 description 3
- 230000005923 long-lasting effect Effects 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 239000000203 mixture Substances 0.000 description 2
- 239000003208 petroleum Substances 0.000 description 2
- 238000004904 shortening Methods 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 241001344923 Aulorhynchidae Species 0.000 description 1
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000001364 causal effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 229910052956 cinnabar Inorganic materials 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 230000008034 disappearance Effects 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000005055 memory storage Effects 0.000 description 1
- 230000008092 positive effect Effects 0.000 description 1
- 238000005381 potential energy Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 230000001105 regulatory effect Effects 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
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
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Geometry (AREA)
- Operations Research (AREA)
- Evolutionary Computation (AREA)
- Algebra (AREA)
- Computer Hardware Design (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明属于应用基础分析领域，公开了一种非均匀粘声波在无限域内传播模型的构建方法，涉及一种适用于时域伽辽金连续有限元波动模拟的粘声性波动方程，同时基于一种高精度时域伽辽金连续有限元—勒让德谱元建立了这一方程的显式时空解耦离散方案。本发明创新引入了串联的广义标准线性体以建立了方程中复柔度的有理函数逼近式，有效降低了所构建时域粘声性波动方程有限元离散过程中存在的高计算存储，解决了以往时域粘声性波动方程难以适用于时域伽辽金连续有限元波动模拟的困难，对分析高精度的海底地形、地貌探测、海底矿产位置与储量的确定手段具有重要价值，可为我国海洋强国战略的实施提供坚实的科技支撑。
Description
技术领域
本发明属于数据处理技术领域，尤其涉及一种非均匀粘声波在无限域内传播模型的构建方法。
背景技术
目前，在频域，粘声波动方程和声波动方程的形式一致，其中粘声波动方程的体积模量为复数，而声波动方程的为实数，因此，在频域内更容易实现粘声波动的模拟，然而频域方程需采用直接或迭代求解，因而需要更多的计算存储和计算量，对于大尺度的三维问题，其难以有效应用。
对于时域方程而言，目前的工作主要有：
现有技术一基于Kjartansson的常Q本构模型建立了空间分数阶导数的粘声方程，其中空间分数阶导数采用伪谱法计算，若采用有限差分或有限元等离散，需对空间分数阶导数进行近。但相比于基于广义标准线性体建立的粘声波动方程，其所需计算量更大。
现有技术二基于广义标准线性体本构建立了粘声波动方程，方程中体积应变为唯一未知量，其方程难以直接采用有限元或谱元离散，因其弱形式方程中边界积分项非自然边界条件，难以直接处理。
现有技术三采用时间二阶空间八阶精度的有限差分格式离散压强形式粘声波动方程，但该方程对于需分析分层非均匀介质时，难以采用有限元或谱元离散。
由于其简单性和实用性，用于模拟声衰减介质中波传播的波动方程的标量粘声近似被广泛用于海洋声学，地球物理学和医学成像。当使用这种近似时，粘声介质的经典特征是其密度和频率相关的体积模量。由介质的粘性性质引起的衰减效应表现为波耗散及其因果关系，即波散。通常，衰减的幅度是使用品质因数Q来衡量的，品质因数Q是复数体积模量的实部和虚部之间的比率。
在频域中，根据对应原理，除了本构关系的体积模量很复杂之外，粘声波方程类似于声波方程。它的解决方案取得了一些丰硕的成果。然而，在某些情况下，可能会出现困难。例如，对于一个大型三维问题，使用直接或迭代求解器的关联线性系统的解决方案需要存储大量的内存，因此需要大量的计算工作。此外，当需要宽带解决方案时，频域解决方案的效率降低。使用时域技术可以部分缓解这些问题。这就是为什么它们通常优于频域技术的原因。
但是，对于一般的复数模量，时域技术的发展是困难的。主要原因是存在卷积积分，其计算需要存储压力和体积应变的整个时间历史。然而，使用适当的Q模型可以避免这些积分的计算。一种常用的模型是Kjartansson的常数Q模型(CQM)。在该模型中，假设Q的值恒定且与频率无关。本构关系包含分数时间导数，其直接计算仍需要存储压力和体积应变的整个历史。朱和哈里斯提出了在空间中使用两个分裂分数拉普拉斯算子的分数时间导数的近似值。然而，分数拉普拉斯算子的数值实现特别复杂，尤其是当分析在具有非平滑空间分布Q的非均质介质中波传播的仿真时。各种方法，例如局部均质逼近，低秩为了解决这个问题，已经提出了近似或混合域近似。然后可以使用有限差分法，频谱/伪谱法或最近开发的有限元方法对近似分数阶粘声波方程进行时空离散。但是，处理分数导数的计算成本通常很高，这就是为什么通常首选基于广义标准线性实体(GSLS)的时域技术的原因。
Liu提出了GSLS模型，该模型由几个平行连接的标准线性实体(SLS)建立，用于近似粘性介质。对于给定的频带，通过增加SLS的数量，可以使用线性或非线性全局优化算法来计算其特征，可以使基于GSLS的Q模型任意接近于参考Q(取决于频率或独立于参考Q)。使用Carcione等人的方法，可以使用GSLS的本构关系导出两种类型的粘声波方程，将体积应变或压力作为标量未知变量。由于体积应变是位移的发散，因此选择体积应变作为未知变量使得难以明确地分析介质界面处的边界条件，即压力和正常粒子位移(速度或加速度)的连续性。使用压力公式，Bai等。开发了一种时空解耦有限差分方案，该方案在时间上精确到二阶，在空间上精确到八阶，用于求解粘声波方程。按照他们的方法，也可以开发高阶有限差分或频谱方法。但是，由于压力公式的推导基于涉及体积模量的本构关系，因此所产生的波动方程包含随介质性质而变化的核与压力发散的卷积，因此阻碍了使用连续有限或频谱元素方法。据作者所知，仍然没有时空解耦方案来模拟粘声波在任意异质介质中的传播，并明确地包含了边界和界面条件。
综上所述，现有技术存在的问题是：
(1)现有技术无法解决复柔度形式的弹簧阻尼模型在单自由度问题中出现的瞬态增长问题。
(2)现有粘声波动方程计算量大，难以直接处理，难以采用有限元或谱元离散。
(3)现有使用有限差分法，频谱/伪谱法或最近开发的有限元方法对近似分数阶粘声波方程进行时空离散。处理分数导数的计算成本通常很高，准确性差。
发明内容
针对现有技术存在的问题，本发明提供了一种非均匀粘声波在无限域内传播模型的构建方法。
本发明是这样实现的，一种非均匀粘声波在无限域内传播模型的构建方法，包括：
步骤一，基于衰减本构关系与其对应的弹性本构关系形式一致，将位移势标量声波动方程中柔度由实数代替为复数，得到粘声波动方程；
步骤二，利用若干串联链接的标准线性体给出波动方程的本构模型；通过非线性优化算法根据任一目标品质因子Q，并求解其各个弹簧阻尼系数，准确计入介质衰减信息；
步骤三，引入辅助变量，通过傅里叶反变换由柔度形式本构给出的粘声频域波动方程，得到对应的可直接用于有限元或谱元离散的时域方程。
进一步，所述非均匀粘声波在无限域内传播模型的构建方法海包括在内域方程的基础上，利用复坐标延拓其弱形式粘声波动方程，建立与之相匹配的完美匹配层。
具体的，所述非均匀粘声波在无限域内传播模型的构建方法包括：
步骤一，基于衰减本构关系与对应的弹性本构关系形式的一致性，将位移势标量声波动方程中的柔度由实数代替为复数，构建复柔度形式的频域粘声性波动方程；
步骤二，基于串联的广义标准线性体建立复柔度的有理函数逼近式，并对逼近式中的参数建立非线性优化计算方法；
步骤三，通过引入辅助变量，通过傅里叶反变换由柔度形式本构给出的粘声频域波动方程，得到适用于时域伽辽金连续有限元模拟的粘声性波动方程。
本发明的另一目在于提供一种执行所述非均匀粘声波在无限域内传播模型的构建方法的非均匀粘声波在无限域内传播模型的构建系统。
本发明的另一目在于提供一种接收用户输入程序存储介质，所存储的计算机程序使电子设备执行所述非均匀粘声波在无限域内传播模型的构建方法。
本发明的另一目在于提供一种存储在计算机可读介质上的计算机程序产品，包括计算机可读程序，供于电子装置上执行时，提供用户输入接口以实施所述非均匀粘声波在无限域内传播模型的构建方法所述的方法。
本发明的另一目在于提供一种执行所述非均匀粘声波在无限域内传播模型的构建方法的海底地形、地貌探测、海底矿产位置与储量探测设备。
综上所述，本发明的优点及积极效果为：本发明构建了基于复柔度形式本构的标量粘性波动方程，其中复柔度本构由串联连接的标准线性体给出，通过引入辅助变量规避了时域方程中所需的卷积计算。本发明新构建的公式同时解决了复柔度形式的弹簧阻尼模型在单自由度问题中出现的瞬态增长。相比于以往有限元离散公式，本发明可分析非均匀介质分界面，同时更节省计算存储。此外，在这一内域方程的基础上推导了与之相匹配的完美匹配层。数值结果表明，新构建的粘声波动内域以及完美匹配层具有极高的精度和数值稳定性，可进一步推广至粘弹性介质及流固耦合介质完美匹配层，同时保证数值计算的长持时稳定。
本发明公开了一种适用于时域伽辽金连续有限元波动模拟的粘声性波动方程，同时基于一种高精度时域伽辽金连续有限元—勒让德谱元建立了这一方程的显式时空解耦离散方案。具体工作包括：构建了复柔度形式的频域粘声性波动方程；基于串联的广义标准线性体建立了复柔度的有理函数逼近式；对逼近式中的参数建立了一种非线性优化计算方法；同时为规避频域粘声性波动方程变换至时域后存在卷积而导致的计算困难，通过引入辅助变量，得到了适用于时域伽辽金连续有限元模拟的粘声性波动方程。本发明创新引入了串联的广义标准线性体以建立了方程中复柔度的有理函数逼近式，有效降低了所构建时域粘声性波动方程有限元离散过程中存在的高计算存储，解决了以往时域粘声性波动方程难以适用于时域伽辽金连续有限元波动模拟的困难，提供了一种新的辅助变量引入方式，解决了新构建时域粘声性波动方程有限元离散过程中存在的数值失稳问题。本发明对分析高精度的海底地形、地貌探测、海底矿产(如可燃冰、石油)位置与储量的确定手段具有重要价值，可为我国海洋强国战略的实施提供坚实的科技支撑。
相比于现有技术，本发明的优点进一步包括：
本发明构建了基于复柔度形式本构的标量粘性波动方程，其中复柔度本构由串联连接的标准线性体给出，通过引入辅助变量规避了时域方程中所需的卷积计算。本发明新构建的公式同时解决了复柔度形式的弹簧阻尼模型在单自由度问题中出现的瞬态增长。相比于以往有限元离散公式，本发明可分析非均匀介质分界面，同时更节省计算存储。此外，在这一内域方程的基础上推导了与之相匹配的完美匹配层。数值结果表明，新构建的粘声波动内域以及完美匹配层具有极高的精度和数值稳定性，可进一步推广至粘弹性介质及流固耦合介质完美匹配层，同时保证数值计算的长持时稳定。
本发明创新引入了串联的广义标准线性体以建立了方程中复柔度的有理函数逼近式，有效降低了所构建时域粘声性波动方程有限元离散过程中存在的高计算存储，解决了以往时域粘声性波动方程难以适用于时域伽辽金连续有限元波动模拟的困难，提供了一种新的辅助变量引入方式，解决了新构建时域粘声性波动方程有限元离散过程中存在的数值失稳问题。本发明对分析高精度的海底地形、地貌探测、海底矿产(如可燃冰、石油)位置与储量的确定手段具有重要价值，可为我国海洋强国战略的实施提供坚实的科技支撑。
本发明通过其频率相关的复数体积柔量来代替常规使用的复数体积模量来表征粘声介质；使用串行连接的标准线性实体(SSLS)建立机械模型，以获取复杂体积顺应性的合理近似值，其参数是通过自适应非线性优化方法计算；利用获得的基于体积柔量的本构关系，在频域中导出了一个新的二阶粘声波方程，其中的弱公式可以物理上解释为虚拟功方程，因此可以使用连续谱元方法离散化空间。此外，引入了一种新方法来解决傅里叶逆变换中涉及的卷积项，然后可以使用显式时间方案实现精确的时间积分，从而避免了经典方法中存在的瞬态增长。最终的完整时空解耦方案可以处理任意异构介质中的波传播。此外，为了处理波在无限域中的传播，推导了弱公式中的完美匹配层，以截断虚拟功方程的通孔复数坐标拉伸。只需稍作修改，就可以使用与截断域内的波动方程式相同的时间方案来实现最终的完美匹配层。通过数值算例证明了该新方案的准确性，数值稳定性和通用性。
附图说明
图1是本发明实施例提供的非均匀粘声波在无限域内传播模型的构建方法流程图。
图2是本发明实施提供的左：包含L个SLS的GSLS并联连接，右图：SSLS包含L个SLS的串联连接图。
图3是本发明实施提供的使用GSLS(蓝线)获得Qref＝100(左)和Qref＝5(右)的频率相关品质因数的倒数和SSLS(黑线)使用非线性优化中的L＝2(上部)，L＝4(中心)和L＝6(底部)SLS的数量频段[0.1Hz，10Hz]。红线是参考模型的倒数Q曲线。水平轴具有对数刻度图。
图4是本发明实施提供的Qref(ω)＝100(ω/ω0)-0.2(左)和Qref(ω)＝5(ω/ω0)-0.2(右)ω0＝2πf0，f0＝1Hz使用GSLS(蓝线)和SSLS(黑线)获得的频率相关品质因数的倒数，其中L＝2(上部)，L＝4(中心)和L＝6(底部)来自[0.1Hz，10Hz]频带中非线性优化的SLS。红线是参考的倒数Q曲线模型。水平轴具有对数刻度图。
图5是本发明实施提供的常数Q模型的频率相关的相速度曲线，其中Qref＝100(左)和Qref＝5(右)使用GSLS(蓝线)和SSLS(黑线)，使用来自非线性的L＝2(上)，L＝4(中心)和L＝6(下)的SLS[0.1Hz，10Hz]频段的优化。红线是参考模型的相速度。横轴有一个对数刻度图。
图6是本发明实施提供的频率相关的Q模型Qref(ω)＝100(ω/ω0)-0.2(左)和Qref(ω)＝5(ω/ω0)-0.2 ω0＝2πf0，f0＝1Hz(右)的频率相关的相速度曲线，使用GSLS(蓝线)和SSLS(黑线)，使用多个L＝2获得(上部)，L＝4(中心)和L＝6(底部)SLS源自[0.1Hz，10Hz]频带中的非线性优化。横轴具有对数刻度图。
图7是本发明实施提供的计算域Ω划分为弯曲元素(二维的四边形和三维的六面体)对于SEM的形状适合体积和媒体界面的边缘图。
图8是本发明实施提供的PML截断的无限域的示意图。
图9是本发明实施提供的单自由度阻尼系统的位移时间历史。半解析解可以通过直接获得解方程(80)通过傅立叶变换，通过用Newmark时间积分求解方程(83)，获得了新的数值解方法，并用Newmark时间积分方法结合方程式(81)求解方程(81)，得到经典的数值解。卷积的二阶离散化方案。粘声介质通过SSLS近似，L＝2(上)，L＝4(中心)，L＝6(底部)。由于该方案的高精度，半解析曲线该解与新颖的数值解的解一致。
图10是本发明实施提供的高斯源时间情况下(3000m，0m)(左列)和(5000m，0m)(右列)的压力时间序列在(0m，0m)处施加功能压力源图。蓝线显示了本发明方案的数值解，与半解析相比Q＝100时的解(红线)。绿线表示Q＝5时本发明的方案的数值解。无论是数值解还是半解析解，通过SSLS对粘声介质进行近似处理，L＝2(上部)，L＝4(中心)，和L＝6(底部)。由于该方案的高精度，半解析解的曲线与数值解是一致的图。
图11是本发明实施提供的500个接收器的等压压力时间序列在(200，2933.33)和(6200，2933.33)之间等距。压力表示为Pa。左轴是以秒为单位的时间图。
图12是本发明实施提供的计算域中的压力快照。时间步长从左上角的600到右下角的2400，每600个时间步的间隔。压力色条与图10右侧所示的色条相同(以Pa为单位)。两个值的单位水平和垂直轴均为米图。
图13是本发明实施提供的计算域内总能量和动能的时间演化，用于模拟a的传播半无限域中的非均质粘声波图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
针对现有技术存在的问题，本发明提供了一种非均匀粘声波在无限域内传播模型的构建方法，下面结合附图对本发明作详细的描述。
如图1所示，本发明实施例提供的非均匀粘声波在无限域内传播模型的构建方法包括：
S101，基于衰减本构关系与对应的弹性本构关系形式的一致性，将位移势标量声波动方程中的柔度由实数代替为复数，构建复柔度形式的频域粘声性波动方程；
S102，基于串联的广义标准线性体建立复柔度的有理函数逼近式，并对逼近式中的参数建立非线性优化计算方法；
S103，通过引入辅助变量，通过傅里叶反变换由柔度形式本构给出的粘声频域波动方程，得到适用于时域伽辽金连续有限元模拟的粘声性波动方程。
下面结合具体实施例对本发明的技术方案作进一步说明。
1、基于顺应性的位移电位粘声波方程
首先简单回顾根据位移电势对声波方程的推导。从频域中涉及顺应性(体积模量的倒数)的本构关系开始，并使用新的SSLS模型对其进行近似，从而获得了新的粘声波方程。
与标准情况一样，本发明证明了其弱形式可以解释为一个虚拟的工作方程。另外，对于质量因子Q取决于频率的情况，可以采用非线性优化的方法来调整SSLS的参数。
1.1声波方程和相应的虚功方程
对于非衰减介质，根据位移推导声波方程潜力和它的弱公式，本发明在这里作了小幅修改。
本发明证明，这样的弱公式确实是一个虚拟的工作方程式。对于以虚拟的工作方程为特征的声学介质密度ρ和体积顺应性即体积模量的倒数κ，压力为：
其中，p表示压力，是梯度算子，u是粒子位移矢量，f是外部压力身体负荷。
动量守恒方程为：
为了求解方程(1)和(2)，本发明引入位移势χ：
其中的点表示相对于时间的微分。假如说：
插入式(3)进入(2)，并进行时间积分，本发明得到：
代入等式(5)和(3)式(1)，本发明根据位移势得到声波方程：
公式(6)通过将其乘以标量检验函数χ得到有限体积Ω：
其中，Γ表示体积的边界，n表示垂直于边界的单位。解释φ可以将虚拟压力服从公式(2)，视为虚拟惯性力的体积密度。和共享体积应变的物理尺寸，即等式中的第一项和第二项。
因此可以看作是内部的虚拟工作。第三项可以理解为通过评估的边界牵引力完成的外部虚拟功整个边界，以及第四项作为由人体载荷完成的外部虚拟功。
对于异质介质中的波传播，内部的法向位移和压力是连续的两个域之间的接口：
另外，沿着外边界(即流体的自由表面)，电位遵循自由表面条件：
此外，根据工作原理，虚拟压力是与边界或界面边界的约束是：
φ＝0. (11)
插入式(8)-(11)至式(7)中，第三项沿不同域之间的接口抵消。因此在一般的非均质介质中使用连续光谱元素法是直截了当。
1.2频域中的粘声波方程和相应的虚功方程
在频域中，粘声介质中基于依从性的本构关系可以写为：
其中符号上的帽子代表其傅立叶变换，而代表频率相关的体积柔量，这是复数模的倒数M(ω)：
为了简化表示法，稍后将省略对ω的依赖。代入式的傅立叶变换。(5)和(3)式(11)，本发明在频域中获得了粘声波方程：
通过引入与上述相同的虚拟压力测试函数，本发明获得了粘声中的虚拟功方程介质：
出于与第2.1节末尾讨论的相同的原因，等式的实现。(15)一般使用连续光谱元素方法的非均质介质很简单(请参见第3节)。用于模拟异质介质中的电磁波，公式为(15)优于基于体积的经典公式模，其波动方程为：
给定连续的频域测试函数方程式的弱公式(16)可以写成：
通过部分集成读取：
分析等式中的第一个方程(8)和的连续性表明不连续的体积模量M阻碍了标准有限元或谱元中界面边界条件的经典实现方法；即等式中的第三项(18)不会沿不同域之间的接口抵消。而且，本发明无法为测试函数和等式中的术语找到清晰的物理类比。
1.3SSLS模型用于复杂总体合规性的合理近似
将复杂的，频域的，批量合规性表示为一种合理的功能，对于以后的开发式的时域对应(15)允许时空离散化。GSLS的直接应用(图2，左)近似给出不合理的函数：
其中，MR是松弛模量，L是GSLS中SLS的数量。τεl和τσl(l＝1，2，…，L)分别是应变(或声学情况下的体积应变)松弛时间和应力(或声学情况下的压力)松弛第lth SLS的时间。式中M的详细推导。为了获得复杂的合规性作为有理函数的形式，本发明在这里采用由串联连接的SLS组成的线性流变模型(图2，右)：SSLS。
与GSLS情况相反，SSLS的频域一致性可以表示为有理函数(请参见附录A)：
其中，
τσl＝ηl/δMl，τεl＝ηl[1/δMl+1/(LMR)]， (21)
MRl＝L×MR和l＝1，2…，L。根据O'Connell，SSLS的有效Q定义为：
插入式(20)式(22)，获得：
其中，
看来，等式(23)具有与等式相同的表达式(28)。因此，本发明调整了他们的非线性优化方法以适合SSLS的参数。同样，本发明在kl和θl保证总能量SSLS的衰减。在本发明的例子中，约束是：
kl≥0 and θl≥0. (25)
然后本发明重写式(23)为：
在理想情况下，有效Q-1(ω)SSLS的数量应等于目标引用用于反转在[ωmin，ωmax]中定义了一组K个角频率最小最大它们以对数尺度线性分布为：
然后本发明的目标误差函数定义为：
此外，为了遵守等式中给出的约束。(22)，本发明介绍k′l与θ′l和kl＝k′l 2等式θl＝θ′l 2(28)。经过这些操作后，可以采用与[14]中相同的方法来应用SolvOpt算法，以获取所有所需的信息。参数。图3显示了有效在角频带Q-1中包括N＝2、4和6个SLS的SSLS中的[0.1，10Hz]，给定和并将其与有效值进行比较Q-1组成的GSLS相同数量的SLS。使用(通常>2)具有约束的非线性最优化给出的曲线几乎是恒定的，并且拟合精确值很好。使用GSLS或SSLS的近似质量几乎相同。此外，随着τεl和τσl的价值MR可通过拟合给定频率下的波相速度来计算。
图5是常数Q模型的频率相关的相速度曲线，其中Qref＝100(左)和Qref＝5(右)使用GSLS(蓝线)和SSLS(黑线)，使用来自非线性的L＝2(上)，L＝4(中心)和L＝6(下)的SLS[0.1Hz，10Hz]频段的优化。红线是参考模型的相速度。横轴有一个对数刻度。
图6是频率相关的Q模型Qref(ω)＝100(ω/ω0)-0.2(左)和Qref(ω)＝5(ω/ω0)-0.2ω0＝2πf0，f0＝1Hz(右)的频率相关的相速度曲线，使用GSLS(蓝线)和SSLS(黑线)，使用多个L＝2获得(上部)，L＝4(中心)和L＝6(底部)SLS源自[0.1Hz，10Hz]频带中的非线性优化。横轴具有对数刻度。
图7是计算域Ω划分为弯曲元素(二维的四边形和三维的六面体)。对于SEM的形状适合体积和媒体界面的边缘。
2、粘声波方程的时空解耦方案
在本节中，本发明首先重点介绍等式离散化的关键要素。在频域中使用连续光谱元素法(SEM)，提出一种新的方法来解决在时域对角线上产生的卷积。频域SEM方程。推导该方法时，此方法比经典方法更准确时域PML。
2.1使用SEM在频域中对粘声波方程进行空间离散
根据Komatitsch和Tromp，本发明通过将体积Ω细分为非重叠来生成网格六面体(在三维情况下)或四边形(在二维情况下)元素Ωe，其中e＝1，2，…，ne这样(图4)。每个六面体体积元素Ωe可以映射到参考多维数据集。
此参考立方体内的点由坐标ξ＝(ξ，η,ζ)的向量表示，其中-1≤ξ≤1，-1≤η≤1和-1≤ζ≤1。六面体元素内的点与参考立方体之间的映射可以写为：
其中，na表示用于定义元素几何的锚节点总数(至少8个角节点)是六面体元素所必需的)。的元素给定元素内的体积dxdydzΩe与参考多维数据集中的体积dξdηdζ的元素相关：
其中，J表示体积雅可比矩阵J的行列式。然后，在元素水平上，等式中的第一项可以写成：
频域位移势可以在一个元素上扩展为：
其中lσ(ξ)表示度数的拉格朗日多项式nl，在nl+1选择的高斯-洛巴托-莱根特雷上定义沿ξ方向指向。为了方便起见，多项式稍后将nl省略为上和符号。本发明还将虚拟压力和扩展为：
代入等式(33)至(35)(31)和(32)并介绍了高斯-洛巴托-莱根特雷数值整合，本发明有：
插入式(20)成等式(36)和(37)，本发明有
因此，由于及其任意性，在由Ne元素，本发明有
其中，
在这里，LTGR是指元素的局部索引与网格尺度下的全局索引之间的关系。细节从等式左侧的第二项计算(15)。
2.2、在时域求解SEM方程的显式时间积分方法
等式的时域对应。(40)是使用傅立叶逆变换以与时域PML源自。本发明先写
其中星号是卷积算符，H(t)是Heaviside函数。插入式(44)和(45)式(40)阅读
其中，
通过结合常用的预测-校正，可以及时离散(46)二阶Newmark方案和以下二阶离散卷积方案用于计算和
其中，和Δt是离散时间步长。但是，如第5节所述，离散时间卷积的方式会在幅度和相位上引入误差。而且，这样的错误不能通过缩短时间步长删除。为了解决这个问题，本发明建议引入一个内存变量。让本发明先将(40)重写为：
然后，通过逆傅立叶变换：
因为fi是已知的，并且相对容易计算，所以(50)的右部分被视为已知值Fi。本发明将(52)改写为：
其中，等式的时间离散化。(54)可以通过采用广泛使用的预测校正二阶Newmark方案：
等式(60)建议，对于每个节点，必须将局部矩阵求逆，这将很耗时。但是这个通过存储节点位置矩阵Mi+ΔtCi/2的逆，可以避免明显的问题这是
A的倒数可以通过解析方式计算为：
值得一提的是，如果是相同的，则可以大大减少内存变量的数量、等式的时间离散化，也可以使用最近开发的高阶低耗散进行低色散Runge-Kutta方案。
3、推导无穷域的弱形式PML
由于PML被公式化为波动方程从实际空间域到波场的解析延续。复杂的对应物，拉伸复杂坐标的技术已成为推导的主要方法PML。
传统上，PML是从波动方程的坐标的扩展形式以强形式得出的，但是它们也可以从波动方程以虚拟形式导出。在经典方法中，PML内部的波场方程及其边界和/或界面条件的推导是独立完成的两者可能无法正确匹配，从而导致数值不稳定和数值精度下降。但是，在虚拟形式方法中，自然可以避免这种不匹配，因为自由和接口同时以一致的方式扩展条件，因此获得的PML呈弱形式拉伸变量的物理含义。但是，仍可以使用连续SEM直接离散化它。
在以下各节中，本发明首先从虚拟中的频域粘声波方程式推导出PML形成。其频谱元素实现的显着点(包括明确的时间方案离散化)将然后突出显示。请注意，只有转角PML会被处理，因为它很容易退化为边缘和表面PML。本发明假设坐标的原点位于角PML和计算域。
在虚拟形式的波动方程经过复数坐标拉伸后，实际中的弱形式PML通过逆坐标变换获得坐标系。傅立叶逆变换最终产生弱时域中的PML形式。
因此，本发明首先将复杂的拉伸坐标引入：
其中表示拉伸函数，通常表示为：
在此，di(xi)是阻尼因子。αi(xi)表示频率偏移因子，可确保PML内部的波取决于频率，因此有效地提供了Butterworth滤波器，κi(xi)是比例因子可以用来改善e逝波的衰减。最后，βi(xi)＝αi(xi)+di(xi)/κi(xi)(不加总)。
映射方程(15)进入复杂空间，本发明得到：
其中，是新坐标系中的体积元素。之所以消失，是因为由于PML的完美匹配特性，它在计算域中沿着PML的中间接口和内部接口抵消了。术语也被忽略，因为PML中的然后，
根据等式。公式(70)借助连锁规则，获得：
在SSLS和等式中插入将(72)变成(70)：
在极点(1/τσ1+iω)sxsysz/(1/τεl+iω)，szsy/sx，sxsz/sy的情况下sxsy/sz不相等，本发明可以将它们分解为反复使用：
在存在多个极点的情况下，本发明可以调整PML中的参数，因为已知PML的精度对拉伸函数的参数的微小变化不敏感。常用的参数定义如下：
αi(xi)＝aπf0(1-|xi|/L)， (77)
κi(xi)＝1+(κmax-1)(|xi|/L)2， (79)
其中，f0是源项的主导频率，VP是PML附近物理域中P波的速度，xi是从界面沿法线方向测量的距离，L是PML的厚度，并且R＝0.001。的值介于1.0到2.0之间，而最佳间隔为最大值κmax用于不同型号之间1.0和13.0。本发明将极点分隔为aπf0Δx0/(8L)值，其中Δx0是PML衰减方向上相邻节点之间的最小距离。式的时空离散。(74)与等式相同。(15)除了卷积需要根据
通过引入本发明有：
因此，对于的时间演变PML E t，本发明可以使用等式中的离散卷积方案，(49)或套用低耗散和低耗散Runge-Kutta方案(81)。
4、数值测试
在本节中，将给出三个数值示例，以验证本发明的精度和数值稳定性。无限域中模拟粘声波的公式及其实现。第一个例子说明了本发明新方法在解决卷积方面的优势。在第二个示例中，本发明验证了在简单无限域中传播粘声波的时空解耦方案存在半分析结果。第三个示例说明了本发明的方案在以下情况下的长期稳定性和准确性：应用于异质介质中的波传播。
4.1单自由度阻尼系统
由于整体刚度矩阵的对称性和正性，单元刚度矩阵内部的组装内部等式中的(40)，本发明的方法在解决卷积上的优势可以用跟随单自由度阻尼系统。
其中，m是质量，1/k是松弛的柔度，以及F(ω)的傅立叶变换源函数f(t)。遵循卷积项的经典方法，并注意逆的傅立叶变换，等式的时域对应。公式(82)可以被写为：
相比之下，使用本发明的方法，Eq的时域对应公式(82)，是：
等式中的γl＝m(1-τσl/τεl)/τεl(86)。El(t)的控制方程式是从Fourier换导出的。
方程的解决方案。公式(81)可以使用预测校正Newmark时间积分方法与卷积(方程式(49))的二阶离散化方案，如2.2节所述，而方程式的解(85)只能使用预测校正Newmark时间积分方法获得。给出了SSLS，表示系统的品质因数等于5。变量的初始状态及其变量一阶时间导数设置为零。
f0＝1Hz和t0＝2s。图9显示了使用等式获得的位移。(83)和(85)同时步长为0.01s。为了比较，对方程式的解决方案。(80)是通过傅立叶变换直接计算的，这个解是称为半解析解。如图9清楚所示，经典方法在两个方面都引入了错误幅度和相位是固有的，无法通过缩短时间步长来消除。总体反应是与具有较大阻尼率的半解析解决方案相比，它最初更大。相反，获得的结果使用式(85)与半解析解一致。
图9是单自由度阻尼系统的位移时间历史。半解析解可以通过直接获得解方程(80)通过傅立叶变换，通过用Newmark时间积分求解方程(83)，获得了新的数值解方法，并用Newmark时间积分方法结合方程式(81)求解方程(81)，得到经典的数值解。
卷积的二阶离散化方案。粘声介质通过SSLS近似，L＝2(上)，L＝4(中心)，L＝6(底部)。由于该方案的高精度，半解析曲线该解与新颖的数值解的解一致。
4.2均匀无限域中的粘声波传播
本发明的时空解耦方案的精度是在无限大范围内对粘声波的传播进行测试的域。附录中给出了二维无限均质粘声介质的解析解。本发明将密度设置为ρ＝1200kg/m3，而参考频率为1Hz时的声波速度为cr＝2000m/s。为了测试低衰减介质和高衰减介质，粘声介质的品质因子Q为选择等于100或5。介质使用SSLS近似。
在网孔的中心(0m，0m)处施加声压载荷。源时间函数是高斯函数由...描述
其中，A是信号源的幅度因子t0＝0s，f0＝1Hz，并且频带是0.1–10Hz。网格包含600X600个正方形元素，长度和宽度为Δx×Δz＝200m×200m。根据以前结果，本发明将离散时间步长设置为Δt＝0.005s。每个元素包括25个高斯-洛巴托-莱热德雷整合点。本发明沿着网格的边缘将压力设置为零。网格的长度和宽度为足够大，可以避免波浪在边缘反射。图10显示了在两个接收器处记录的压力时间序列，定位在(5000m，0m)和(3000m，0m)处，可以通过本发明的方案或分析解决方案获得。该方案的数值解与解析结果几乎相同，验证了该方案的准确性。在Q＝5的情况下，信号的非常早到达是由于本发明平均设置了声波低衰减介质和高衰减介质在参考频率下的速度。但是，图10清楚地表明在高衰减介质中，振幅耗散和波色散更为明显。
图10中，高斯源时间情况下(3000m，0m)(左列)和(5000m，0m)(右列)的压力时间序列在(0m，0m)处施加功能压力源。蓝线显示了本发明方案的数值解，与半解析相比Q＝100时的解(红线)。绿线表示Q＝5时本发明的方案的数值解。无论是数值解还是半解析解，通过SSLS对粘声介质进行近似处理，L＝2(上部)，L＝4(中心)，和L＝6(底部)。由于该方案的高精度，半解析解的曲线与数值解是一致的。
4.3半无限域中的异构波传播
第三个例子显示了在传播情况下本发明的方案的准确性和长期数值稳定性半无限域中的非均质粘声波的特性。本发明分析两层吸收介质，即内部其界面由正弦曲线描述，z(x)＝2400+180sin(12πx/6400)。上层的顶表面是平坦的，并且那里的压力被设置为零。准时压力源是应用于(2908.33m，3100m)，其源时间函数设置为Ricker小波(高斯的二阶导数)由等式表示。(86)与t0＝0s and f0＝10Hz，通常在海洋声学模拟中使用。本发明近似具有SSLS的粘性介质，在[1，100Hz]频带中具有四个SLS。为了与常规做法保持一致，参考波速是在无限频率下给出的，即在非松弛波下给出的速度。出于实际原因，源函数要乘以1.0×109，以及数值的开始时间模拟选择为-0.12s。PML用于无限域的截断。截断的网格包括nx×nz＝144×108＝15552元素，其总水平长度和深度等于6400m×4800m。四阶多项式用于每个光谱元素内的Lagrange插值。截断的网格包括沿左侧，右侧和底部边缘的三元素厚PML。本发明设置离散时间升至Δt＝0.0008s。在(200，2933.33)和(6200，2933.33)之间等距的500个接收器上记录该字段。图11显示了在所有接收器处模拟的压力时间序列。本发明清楚地看到直接到达和衍射在弯曲的界面波动来自自由表面的反射以及以后的反射和衍射，相比较而言很小到直接到达，但仍值得注意。图12显示了在时间步长600处获得的压力波场的一些快照，1200、1800和2400，它们清楚地表明PML以高精度吸收了入射波。几乎没有沿截短域的PML接口观察到反射。图13显示了能量随时间的衰减计算域以检查本发明方案的长期稳定性。总能量是动能的总和和势能。总共计算了200,000个时间步长，并且解决方案保持稳定。因此，本发明对异质粘声波传播的时空解耦方案是准确，高效的，并且数值稳定的。
5、通过根据其频率相关的复杂整体顺应性来表征粘声介质，而不是通过传统上使用的复数体积模量。本发明推导了频域粘声波方程，用于分析波在异质介质中的传播。它的弱形式对应于一个虚拟的工作方程，因此可以是在空间中使用连续SEM离散化。使用以下公式建立了复杂散装法规遵从性的合理近似值：由几个串联的SLS组成的机械模型。在时域，本发明提出了一种新的方法解决在频域SEM方程的时域对立部分出现的卷积项。这个这种方法避免了不稳定性，并允许推导新的PML公式以截断计算通过虚拟工作方程的复杂坐标拉伸来确定域。生成的PML公式可以是使用与截断域内的波动方程几乎相同的方案来实现。通过几个示例，本发明证明了时空解耦的准确性，长期数值稳定性和多功能性无限域中非均质粘声波传播的方案，现在正在努力扩展当前在非均质粘声/粘弹性耦合介质中传播的方案。
下面结合具体实验数据对本发明的技术方案作进一步说明。
本发明利用衰减本构关系与其对应的弹性本构关系形式一致这一性质，将位移势标量声波动方程中柔度由实数代替为复数，即可得到粘声波动方程，为了便于波动方程的数值离散，其本构模型由若干串联链接的标准线性体给出，通过非线性优化算法根据任一目标品质因子Q求解其各个弹簧阻尼系数，从而准确的计入介质衰减信息。通过引入辅助变量，通过傅里叶反变换由柔度形式本构给出的粘声频域波动方程，得到其对应的时域方程，从而可直接用于有限元或谱元离散。此外，在这一内域方程的基础上，利用复坐标延拓其弱形式粘声波动方程，建立了与之相匹配的完美匹配层。
针对新构建的基于柔度本构的粘声内域以及完美匹配层，进行了以下数值实验：
通过对比全空间弱衰减及强衰减介质半解析解，验证了新构建的内域离散格式的高精度。
分别将新格式与旧格式应用至单自由度体系，其结果表明旧格式存在瞬态增长现象，新格式与半解析解几乎重合。
在分析较复杂介质分界面的两层粘声介质中，采用新构建的完美匹配层，其数值结果验证了新构建的完美匹配层可高精度无反射的应用于分层介质，同时保证数值计算的长持时稳定。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (7)
1.一种非均匀粘声波在无限域内传播模型的构建方法，其特征在于，所述非均匀粘声波在无限域内传播模型的构建方法包括：
构建复柔度形式的频域粘声性波动方程；基于串联的广义标准线性体建立复柔度的有理函数逼近式；对逼近式中的参数建立非线性优化计算方法；同时通过引入辅助变量，得到适用于时域伽辽金连续有限元模拟的粘声性波动方程；
基于获取的粘声性波动方程对海底地形、地貌探测、海底矿产位置与储量进行分析；
所述非均匀粘声波在无限域内传播模型的构建方法进一步包括：
步骤一，基于衰减本构关系与对应的弹性本构关系形式的一致性，将位移势标量声波动方程中的柔度由实数代替为复数，构建复柔度形式的频域粘声性波动方程；
步骤二，基于串联的广义标准线性体建立复柔度的有理函数逼近式，并对逼近式中的参数建立非线性优化计算方法；
步骤三，通过引入辅助变量，通过傅里叶反变换由柔度形式本构给出的粘声频域波动方程，得到适用于时域伽辽金连续有限元模拟的粘声性波动方程；基于获取的粘声性波动方程对海底地形、地貌探测、海底矿产位置与储量进行分析；
步骤一中，所述复柔度形式的频域粘声性波动方程的构建方法为：
(1)声波方程和相应的虚功方程
对于以虚拟的工作方程为特征的声学介质密度ρ和体积顺应即体积模量的倒数κ，压力为：
其中，p表示压力，是梯度算子，u是粒子位移矢量，f是外部压力身体负荷；动量守恒方程为：
其中，为求解方程(1)和(2)，引入位移势χ：
其中的点表示相对于时间的微分；假如：
插入式(3)进入(2)，并进行时间积分，得到：
代入等式(5)和(3)式(1)，根据位移势得到声波方程：
公式(6)通过将其乘以标量检验函数φ得到有限体积Ω：
其中，Γ表示体积的边界，n表示垂直于边界的单位；解释φ可以将虚拟压力服从公式(2)，视为虚拟惯性力的体积密度；和共享体积应变的物理尺寸；
对于异质介质中的波传播，内部的法向位移和压力是连续的两个域之间的接口：
另外，沿着外边界，电位遵循自由表面条件：
虚拟压力与边界或界面边界的约束：
φ＝0(11)
插入式(8)-(11)至式(7)中；
(2)频域中的粘声波方程和相应的虚功方程
在频域中，粘声介质中基于依从性的本构关系写为：
其中符号上的帽子代表其傅立叶变换，而代表频率相关的体积柔量，复数模的倒数为M(ω)：
通过傅立叶变换将式(5)和式(3)代入式(11)，在频域中获得粘声波方程：
通过引入虚拟压力测试函数，获得粘声中的虚拟功方程介质：
其中，用于模拟异质介质中的电磁波，公式为(15)优于基于体积的经典公式模型，波动方程为：
给定连续的频域测试函数方程式的弱公式(16)写成：
通过部分集成读取：
2.如权利要求1所述的非均匀粘声波在无限域内传播模型的构建方法，其特征在于，GSLS的直接应用近似给出不合理的函数：
其中，MR是松弛模量，L是GSLS中SLS的数量；τεl和τσl(l＝1，2，…，L)分别是应变松弛时间和应力松弛第lth SLS的时间；式中M的详细推导，获得复杂的合规性作为有理函数的形式；
采用由串联连接的SLS组成的线性流变模型SSLS，SSLS的频域一致性可以表示为有理函数：
τol＝ηl/δMl， τεl＝ηl[1/δM l+1/(LMR)] (21)
MRl＝L×MR和l＝1，2…，L；根据O′Connell，SSLS的有效Q定义为：
插入式(20)式(22)，获得：
其中，
等式(23)具有与等式(28)相同的表达式，调整非线性优化方法以适合SSLS的参数；在kl和θl保证总能量SSLS的衰减，约束是：
kl≥0andθl≥0(25)
重写式(23)为：
在理想情况下，有效Q-1(ω)SSLS的数量应等于目标引用用于反转在[ωmin，ωmax]中定义一组K个角频率最小最大以对数尺度线性分布为：
目标误差函数定义为：
此外，应用SolvOpt算法，获取所有所需的信息。
3.如权利要求1所述的非均匀粘声波在无限域内传播模型的构建方法，其特征在于，所述粘声波方程的时空解耦方案为：
(1)使用SEM在频域中对粘声波方程进行空间离散
通过将体积Ω细分为非重叠来生成网格六面体或四边形元素Ωe，其中e＝1，2，…，ne这样每个六面体体积元素Ωe映射到参考多维数据集；
参考立方体内的点由坐标ξx＝(ξ，η，ζ)的向量表示，其中-1≤ξ≤1，-1≤η≤1和-1≤ζ≤1；六面体元素内的点与参考立方体之间的映射写为：
其中，na表示用于定义元素几何的锚节点总数是六面体元素所必需的，
给定元素内的体积dxdydzΩe与参考多维数据集中的体积dξdηdζ的元素相关：
其中，J表示体积雅可比矩阵J的行列式；
频域位移势在一个元素上扩展为：
其中，lσ(ξ)表示度数的拉格朗日多项式nl，在nl+1选择的高斯-洛巴托-莱根特雷上定义沿ξx方向指向；多项式稍后将nl省略为上和符号；
将虚拟压力和扩展为：
将等式(33)、(31)和(32)代入至(35)，有：
插入式(20)至等式(36)和(37)，有：
由于 及其任意性，在由Ne元素，有：
LTGR是指元素的局部索引与网格尺度下的全局索引之间的关系；(2)在时域求解SEM方程的显式时间积分方法
时域对应是使用傅立叶逆变换以与时域PML得到的：
其中星号是卷积算符，H(t)是Heaviside函数；插入得到：
4.如权利要求1所述的非均匀粘声波在无限域内传播模型的构建方法，其特征在于，推导无穷域的弱形式PML的方法，包括：
首先将复杂的拉伸坐标引入：
其中表示拉伸函数，
映射方程(15)进入复杂空间，得到：
其中，是新坐标系中的体积元素，消失，在计算域中沿着PML的中间接口和内部接口抵消；
借助连锁规则，获得：
插入变成
在极点(1/τσl+iω)sxsy sz/(1/τεl+iω)，szsy/sx，sxsz/sy的情况下sxsy/sz不相等，分解为
在存在多个极点的情况下，调整PML中的参数，将极点分隔为aπf0Δx0/(8L)值，其中Δx0是PML衰减方向上相邻节点之间的最小距离；根据：
通过引入有：
对于的时间演变PML Et，使用离散卷积，或套用低耗散和低耗散Runge-Kutta方法。
5.一种执行权利要求1-4任意一项所述非均匀粘声波在无限域内传播模型的构建方法的非均匀粘声波在无限域内传播模型的构建系统。
6.一种接收用户输入程序存储介质，所存储的计算机程序使电子设备执行权利要求1-4任意一项所述非均匀粘声波在无限域内传播模型的构建方法。
7.一种执行权利要求1-4任意一项所述非均匀粘声波在无限域内传播模型的构建方法的海底地形、地貌探测、海底矿产位置与储量探测设备。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110103651.3A CN113221392B (zh) | 2021-01-26 | 2021-01-26 | 一种非均匀粘声波在无限域内传播模型的构建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110103651.3A CN113221392B (zh) | 2021-01-26 | 2021-01-26 | 一种非均匀粘声波在无限域内传播模型的构建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113221392A CN113221392A (zh) | 2021-08-06 |
| CN113221392B true CN113221392B (zh) | 2023-12-19 |
Family
ID=77084571
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110103651.3A Active CN113221392B (zh) | 2021-01-26 | 2021-01-26 | 一种非均匀粘声波在无限域内传播模型的构建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113221392B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113821948B (zh) * | 2021-08-19 | 2024-01-19 | 东南大学 | 一种埋入或浸入式储液管道中导波模态的建模方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2010125405A2 (en) * | 2009-05-01 | 2010-11-04 | Dynamic Dinosaurs Bv | Method and apparatus for applying vibrations during borehole operations |
| CN104965223A (zh) * | 2015-05-29 | 2015-10-07 | 中国石油天然气股份有限公司 | 粘声波全波形反演方法及装置 |
| CN107798156A (zh) * | 2016-09-02 | 2018-03-13 | 赵建国 | 一种频率域2.5维粘弹性波数值模拟方法及装置 |
| CN109190169A (zh) * | 2018-08-02 | 2019-01-11 | 电子科技大学 | 一种三维时域电磁学杂交时域间断伽辽金数值方法 |
-
2021
- 2021-01-26 CN CN202110103651.3A patent/CN113221392B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2010125405A2 (en) * | 2009-05-01 | 2010-11-04 | Dynamic Dinosaurs Bv | Method and apparatus for applying vibrations during borehole operations |
| CN104965223A (zh) * | 2015-05-29 | 2015-10-07 | 中国石油天然气股份有限公司 | 粘声波全波形反演方法及装置 |
| CN107798156A (zh) * | 2016-09-02 | 2018-03-13 | 赵建国 | 一种频率域2.5维粘弹性波数值模拟方法及装置 |
| CN109190169A (zh) * | 2018-08-02 | 2019-01-11 | 电子科技大学 | 一种三维时域电磁学杂交时域间断伽辽金数值方法 |
Non-Patent Citations (3)
| Title |
|---|
| IMPLEMENTATION AND EXPERIMENTAL VALIDATION OF A NEW VISCOTHERMAL ACOUSTIC FINITE ELEMENT FOR ACOUSTO-ELASTIC PROBLEMS;W.M.Beltman等;《Journal of Sound and Vibration》;19980910;第216卷(第1期);159-185 * |
| Improved forward wave propagation and adjoint-based sensitivity kernel calculations using a numerically stable finite-element PML;Zhinan xie等;《Geophysical Journal International》;20140721;第198卷(第3期);1714-1747 * |
| 海域地震动模拟关键技术研究;郑永路;《中国优秀硕士学位论文全文数据库基础科学辑》;20200115(第01期);A010-141 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113221392A (zh) | 2021-08-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Birk et al. | An improved continued‐fraction‐based high‐order transmitting boundary for time‐domain analyses in unbounded domains | |
| Ford et al. | A nonhydrostatic finite-element model for three-dimensional stratified oceanic flows. Part I: Model formulation | |
| De La Puente et al. | Discontinuous Galerkin methods for wave propagation in poroelastic media | |
| Bueno-Orovio et al. | Continuous adjoint approach for the Spalart-Allmaras model in aerodynamic optimization | |
| Askan et al. | Full waveform inversion for seismic velocity and anelastic losses in heterogeneous structures | |
| Chen et al. | A high‐order approach for modelling transient wave propagation problems using the scaled boundary finite element method | |
| Seriani et al. | Spectral element method for acoustic wave simulation in heterogeneous media | |
| Bazyar et al. | Scaled boundary finite-element method for solving non-homogeneous anisotropic heat conduction problems | |
| Cohen et al. | Mixed spectral finite elements for the linear elasticity system in unbounded domains | |
| Schneider et al. | Decoupling simulation accuracy from mesh quality | |
| Proot et al. | Least-squares spectral elements applied to the Stokes problem | |
| Prempramote et al. | High‐order doubly asymptotic open boundaries for scalar wave equation | |
| Vandenhoeck et al. | Implicit high-order flux reconstruction solver for high-speed compressible flows | |
| Daxini et al. | Parametric shape optimization techniques based on Meshless methods: A review | |
| CN101369024A (zh) | 一种地震波动方程生成方法及系统 | |
| CN113221392B (zh) | 一种非均匀粘声波在无限域内传播模型的构建方法 | |
| Lee et al. | Absorbing boundary condition for scalar-wave propagation problems in infinite media based on a root-finding algorithm | |
| Williams et al. | Nodal points and the nonlinear stability of high-order methods for unsteady flow problems on tetrahedral meshes | |
| Hanyga et al. | Wave field simulation for heterogeneous transversely isotropic porous media with the JKD dynamic permeability | |
| Lin et al. | Isogeometric least-squares collocation method with consistency and convergence analysis | |
| Whiteley et al. | Error estimation and adaptivity for incompressible hyperelasticity | |
| AlSalem et al. | Embedded boundary methods for modeling 3D finite-difference Laplace-Fourier domain acoustic-wave equation with free-surface topography | |
| Shevchenko et al. | Boundary and contact conditions of higher order of accuracy for grid-characteristic schemes in acoustic problems | |
| Fu et al. | An efficient high-order multiscale finite element method for frequency-domain elastic wave modeling | |
| Zhang et al. | A reduced high-order compact finite difference scheme based on proper orthogonal decomposition for the generalized Kuramoto-Sivashinsky equation |
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 |