A parallel simulation method for ocean acoustic field based on DR-SETD

Abstract

The application discloses a kind of ocean sound field parallel simulation methods based on DR-SETD, belong to the technical field of simulation calculation.This application includes: three-dimensional fluid-structure interaction model is simplified to axisymmetric model to reduce the scale of calculation;Adopt high-order spectral element method based on orthogonal polynomial to carry out spatial dispersion, and numerical singularity is handled for axisymmetric property;The weak form of time-domain wave equation of fluid domain and solid domain is respectively constructed, and is solved by interface coupling condition;Efficient parallelization of solving process is realized using explicit time integration format and shared memory parallel computing architecture.This application realizes the high-precision, high-efficiency full-wave simulation of fluid-structure interaction sound field containing complex terrain.

Description

Technical Field

[0001] This invention belongs to the field of computational simulation technology, and in particular relates to a parallel simulation method for ocean acoustic fields based on DR-SETD. Background Technology

[0002] Ocean acoustic wave propagation simulation is a key technology for exploring acoustic behavior in complex underwater media, with applications spanning seabed resource exploration, underwater communication, sonar detection, and marine environmental monitoring. With the increasing demands of marine engineering and military applications, more precise modeling of the propagation characteristics of sound waves in complex marine environments is required. The acoustic propagation characteristics in complex marine environments are influenced by various factors, especially the coupling effects of seabed boundary undulations and elastic media. Changes in seabed topography significantly alter the modal distribution of sound waves, leading to energy leakage or diffusion, thus affecting long-range propagation efficiency. Furthermore, the interaction between the elastic seabed medium and the seawater interface excites interface wave modes, which have relatively small attenuation and significantly influence the energy distribution and directionality of very low frequency (VLF) sound propagation. The complexity of these physical phenomena makes the study of VLF sound propagation one of the core challenges in the field of marine acoustics.

[0003] In underwater acoustic propagation research, ray acoustic models and parabolic equation models are two classic methods for handling horizontally varying waveguide problems. Ray acoustic models, based on high-frequency approximations of the wave equation, have a clear physical mechanism and are easy to interpret intuitively, but they have significant limitations. On the one hand, this model cannot describe the phase transition phenomenon at the seawater-seabed interface and the excitation process of interface waves; on the other hand, its high-frequency assumptions result in insufficient computational power for low-frequency sound fields, making it difficult to meet the modeling requirements of very low-frequency sound waves. In contrast, parabolic equation models achieve wavefield propagation simulation through spatial recursion algorithms, achieving some improvement in computational accuracy. However, this method still faces stability challenges when dealing with continuously undulating seabed boundaries. When the seabed topography changes drastically, numerical simulations are prone to divergence, thus limiting its adaptability to complex three-dimensional spatial structures. Furthermore, based on the unidirectional propagation assumption, this method cannot handle the significant backscattering caused by complex topography. The inherent limitations of the two types of models mentioned above—the inability to strike a balance between computational efficiency, adaptability to complex terrain, and the ability to capture full-wave physical phenomena—highlight the urgent need to develop new numerical methods to accurately simulate the propagation process of the full-wave sound field, including backscattering and interface waves. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention proposes a parallel simulation method for ocean acoustic fields based on DR-SETD, thereby resolving the issues present in the prior art.

[0005] Firstly, to achieve the above objectives, this invention provides a parallel simulation method for ocean acoustic fields based on DR-SETD, comprising the following steps:

[0006] The physical parameters and topographic data of the target sea area are obtained. Based on the axisymmetric assumption, the three-dimensional fluid-structure coupled acoustic field model is reduced in dimension to obtain an axisymmetric calculation model.

[0007] The higher-order spectral element method based on orthogonal polynomials is adopted, and two different types of integration nodes are configured on the unit of the axisymmetric computation model to complete the spatial discretization of the computation domain.

[0008] We construct and derive the weak forms of the time-domain wave equations that satisfy the stress and displacement continuity conditions for fluid and solid media respectively, and assemble them to obtain the global coupled system equations.

[0009] The global coupled system equations are solved by time-progression using an explicit time integration scheme. Within each time step, the computational tasks of each discrete unit are processed in parallel and data is synchronized using a shared memory parallel computing architecture.

[0010] Optionally, the process of reducing the dimensionality of the three-dimensional fluid-structure interaction acoustic field model to obtain an axisymmetric computational model includes:

[0011] Input two-dimensional topographic profile data of the target sea area, which includes changes in seawater depth and seabed interface morphology;

[0012] The two-dimensional terrain profile is rotated about one of its vertical sides to generate a three-dimensional axisymmetric computational domain, thereby simplifying the original three-dimensional problem into a computational problem within the two-dimensional profile.

[0013] Optionally, the spatial discretization process can be completed using a higher-order spectral element method based on orthogonal polynomials, including:

[0014] First-type Gauss-Lobart integration nodes are configured in the axial direction of the element;

[0015] Second-type Gaussian-Lobart integral nodes are configured in the radial direction of the element. The second-type nodes are constructed based on modified Legendre polynomials.

[0016] Higher-order interpolation of field variables within the element is performed using two types of nodes, and the corresponding local stiffness matrix and diagonal mass matrix are calculated.

[0017] Optionally, the process of constructing and deriving the weak form of the time-domain wave equation includes:

[0018] For the fluid domain, a scalar velocity potential is introduced, and the displacement field is expressed as the gradient of this potential function, thus deriving the weak form of the wave equation with the potential function as the variable.

[0019] For the solid domain, the weak form of the elastic wave equation is derived directly using the displacement vector as the variable, based on the constitutive relations of elasticity.

[0020] At the shared interface between the fluid and solid regions, a boundary integral term is introduced to achieve forced coupling of the continuity conditions of normal stress and normal displacement.

[0021] Optionally, the process of parallel processing and data synchronization using a shared-memory parallel computing architecture includes:

[0022] The spatial domain decomposition method is used to divide all discretized spectral units into sub-regions whose number matches the number of parallel threads;

[0023] By using OpenMP parallel compilation guidelines, the cyclic computation tasks of each cell within a sub-region are dynamically scheduled and allocated to multiple threads for simultaneous execution.

[0024] At the coupling interface, a thread synchronization mechanism is used to ensure data consistency between adjacent units computed by different threads.

[0025] Optionally, the process of using an explicit time integration scheme for time advancement includes:

[0026] The field variables at each time step are updated using a second-order accurate Newmark-β explicit integral scheme;

[0027] Within their assigned sub-regions, each thread independently and rapidly solves the motion equations of the unit nodes using the local diagonal mass matrix, thus advancing the local time step.

[0028] Secondly, the present invention also provides a parallel simulation system for ocean acoustic fields based on DR-SETD, for implementing a parallel simulation method for ocean acoustic fields based on DR-SETD, the system comprising:

[0029] The axisymmetric modeling module is used to acquire the physical parameters and topographic data of the target sea area, and to simplify the three-dimensional fluid-structure coupled acoustic field model based on the axisymmetric assumption to obtain an axisymmetric calculation model.

[0030] The spectral element discretization module, connected to the axisymmetric modeling module, is used to configure integration nodes on the units of the axisymmetric computation model using a higher-order spectral element method based on orthogonal polynomials to complete the spatial discretization of the computational domain.

[0031] The coupled solution module, connected to the spectral element discretization module, is used to construct weakly coupled forms of the time-domain wave equation for both fluid and solid media, and assemble them to obtain the global system equation for solution.

[0032] The parallel computing module, connected to the coupled solution module, is used to perform parallel processing and data synchronization of the computational tasks of each discrete unit based on shared memory during the solution process using an explicit time integration scheme.

[0033] Thirdly, the present invention also provides a computer terminal device, comprising:

[0034] One or more processors;

[0035] A memory, coupled to the processor, for storing one or more programs;

[0036] When the one or more programs are executed by the one or more processors, the one or more processors implement the steps of the parallel simulation method for ocean acoustic fields based on DR-SETD in the first aspect described above.

[0037] Fourthly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the computer program is executed by a processor, it implements the steps of the parallel simulation method for ocean acoustic fields based on DR-SETD in the first aspect described above.

[0038] Fifthly, the present invention also provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the parallel simulation method for ocean acoustic fields based on DR-SETD described in the first aspect.

[0039] Compared with the prior art, the present invention has the following advantages and technical effects:

[0040] This invention provides a parallel simulation method for ocean acoustic fields based on DR-SETD, effectively overcoming the shortcomings of existing technologies. By simplifying the three-dimensional model to axisymmetric characteristics, this invention significantly reduces the computational scale and, combined with high-order spectral element discretization and shared memory parallel computing, achieves a substantial improvement in computational efficiency. This method exhibits exponential numerical convergence, enabling high-precision simulation of sound wave propagation and fully capturing all wavefield physical phenomena, including backscattering and fluid-solid interface waves. Furthermore, its unstructured mesh based on the spectral element method can flexibly adapt to complex terrains such as inclined seabeds. Overall, this invention achieves synergistic improvements in computational accuracy, efficiency, physical completeness, and adaptability to complex environments, providing a general solution for high-fidelity, high-efficiency simulation of complex ocean acoustic fields. Attached Figure Description

[0041] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings:

[0042] Figure 1 This is a two-dimensional cross-sectional view of a homogeneous flat-bottomed FSI according to an embodiment of the present invention;

[0043] Figure 2The image shows a comparison of the displacement components COMSOL and DR-SETD at receiving point 2 in this embodiment of the invention, where (a) is the horizontal component and (b) is the vertical component.

[0044] Figure 3 The image shows a comparison of the displacement components COMSOL and DR-SETD at receiving point 3 in this embodiment of the invention, where (a) is the horizontal component and (b) is the vertical component.

[0045] Figure 4 These are instantaneous snapshots of the wave field at different times in a homogeneous flat-bottom case according to an embodiment of the present invention. (a) is a snapshot of the direct longitudinal wave propagation at 0.3s, (b) is a snapshot of the wavefront extending into the solid region and the early morphology of the Scholte wave at 0.9s, (c) is a snapshot of the Scholte interface wave formation at 1.08s, and (d) is a snapshot of the mature morphology of the Scholte wave and the boundary absorption effect at 1.5s.

[0046] Figure 5 This is a two-dimensional cross-sectional view of an inclined seabed FSI case according to an embodiment of the present invention;

[0047] Figure 6 The above is a comparison of acoustic propagation transmission loss curves in the case of inclined seabed according to the present invention, where (a) represents the transmission loss of 1-4km and (b) represents the transmission loss of 5-9km.

[0048] Figure 7 These are instantaneous snapshots of the wavefield at different times in the case of the inclined seabed in this embodiment of the invention. (a) is a snapshot of the initial stage of wavefield propagation at 0.49s, (b) is a snapshot of the wavefront interacting with the inclined topography at 1.46s, (c) is a snapshot of the wavefield propagating along the inclined interface and the backscattering effect at 2.32s, and (d) is a snapshot of the long-term wavefield attenuation and numerical stability at 7.32s. Detailed Implementation

[0049] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other. The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0050] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0051] Example 1

[0052] This embodiment provides a parallel simulation method for ocean acoustic fields based on DR-SETD, including:

[0053] The physical parameters and topographic data of the target sea area are obtained. Based on the axisymmetric assumption, the three-dimensional fluid-structure coupled acoustic field model is reduced in dimension to obtain an axisymmetric calculation model.

[0054] The higher-order spectral element method based on orthogonal polynomials is adopted, and two different types of integration nodes are configured on the unit of the axisymmetric computation model to complete the spatial discretization of the computation domain.

[0055] We construct and derive the weak forms of the time-domain wave equations that satisfy the stress and displacement continuity conditions for fluid and solid media respectively, and assemble them to obtain the global coupled system equations.

[0056] The global coupled system equations are solved by time-progression using an explicit time integration scheme. Within each time step, the computational tasks of each discrete unit are processed in parallel and data is synchronized using a shared memory parallel computing architecture.

[0057] As one implementation method in this embodiment, the process of reducing the dimension of the three-dimensional fluid-structure coupled acoustic field model to obtain an axisymmetric computational model includes:

[0058] Input two-dimensional topographic profile data of the target sea area, which includes changes in seawater depth and seabed interface morphology;

[0059] The two-dimensional terrain profile is rotated about one of its vertical sides to generate a three-dimensional axisymmetric computational domain, thereby simplifying the original three-dimensional problem into a computational problem within the two-dimensional profile.

[0060] As one implementation method in this embodiment, the spatial discretization process using the high-order spectral element method based on orthogonal polynomials includes:

[0061] First-type Gauss-Lobart integration nodes are configured in the axial direction of the element;

[0062] Second-type Gaussian-Lobart integral nodes are configured in the radial direction of the element. The second-type nodes are constructed based on modified Legendre polynomials.

[0063] Higher-order interpolation of field variables within the element is performed using two types of nodes, and the corresponding local stiffness matrix and diagonal mass matrix are calculated.

[0064] As one implementation method in this embodiment, the process of constructing and deriving the weak form of the time-domain wave equation includes:

[0065] For the fluid domain, a scalar velocity potential is introduced, and the displacement field is expressed as the gradient of this potential function, thus deriving the weak form of the wave equation with the potential function as the variable.

[0066] For the solid domain, the weak form of the elastic wave equation is derived directly using the displacement vector as the variable, based on the constitutive relations of elasticity.

[0067] At the shared interface between the fluid and solid regions, a boundary integral term is introduced to achieve forced coupling of the continuity conditions of normal stress and normal displacement.

[0068] As one implementation method in this embodiment, the process of parallel processing and data synchronization using a shared memory parallel computing architecture includes:

[0069] The spatial domain decomposition method is used to divide all discretized spectral units into sub-regions whose number matches the number of parallel threads;

[0070] By using OpenMP parallel compilation guidelines, the cyclic computation tasks of each cell within a sub-region are dynamically scheduled and allocated to multiple threads for simultaneous execution.

[0071] At the coupling interface, a thread synchronization mechanism is used to ensure data consistency between adjacent units computed by different threads.

[0072] As one implementation method in this embodiment, the process of using an explicit time integration scheme for time advancement includes:

[0073] The field variables at each time step are updated using a second-order precision Newmark-β explicit integral scheme;

[0074] Within the allocated sub-region of a cell, each thread independently and rapidly solves the motion equations of the cell nodes using the local diagonal mass matrix, thus advancing the local time step.

[0075] This invention discloses a large-scale parallel method for acoustic fields in fluid-structure interaction (FSI) models based on axisymmetric time-domain spectral elements (DR-SETD). It includes: constructing the mathematical theoretical framework for DR-SETD to handle FSI problems, and deriving the weak formal equations and numerical implementation steps under axisymmetric conditions in detail. The reliability and applicability of DR-SETD are verified through two typical cases: firstly, the computational accuracy of DR-SETD is verified based on a homogeneous flat-bottom FSI case; secondly, for a case involving inclined seabed FSI with significant backscattering effects, the advantages of DR-SETD in capturing complex wavefield behavior (such as Stoneley-Scholte interface wave excitation) are demonstrated, illustrating the necessity of full-wave time-domain simulation compared to approximate methods such as parabolic equations. Finally, the performance of DR-SETD in large-scale parallel computing is evaluated from the perspective of computational efficiency. This invention provides a general-purpose simulation tool that combines high computational accuracy with superior large-scale parallel scalability, providing strong technical support for scientific research and engineering design in fields such as acoustic propagation, target detection, and identification in complex marine environments. It can be used for acoustic detection and inversion research of seabed geological structures, and by comparing simulation results with actual observation data, it can deepen the understanding of the wave mechanism of different seabed interfaces.

[0076] Methods and solutions:

[0077] The DR-SETD system constructed in this invention focuses on the actual needs of marine acoustic exploration. Through "high-order precision modeling + efficient parallel computing", it achieves three major goals: "adaptation to complex terrain, full wave field fidelity, and real-time simulation support", providing quantitative basis for sonar system optimization, underwater target identification, and seabed resource positioning.

[0078] The system consists of 7 core modules, which work together to achieve the entire process of "accurate modeling - efficient calculation - reliable output". Its key technological advantages are as follows:

[0079] 1. Marine Topography Input Interface: Input satellite mapping / sonar scanning data (such as seabed topographic maps and water depth profiles) into the exploration target area: Provide an accurate modeling basis for complex terrains such as sloping seabeds / flat bottoms, and avoid sound field distortion caused by terrain simplification in traditional methods.

[0080] 2. Axisymmetric Modeling Engine: Rotates 2D terrain profiles around an axis to generate an axisymmetric 3D model. Dimensionality reduction and compression: Compresses the full 3D computation to 2.5D, reducing the degrees of freedom by 90% (e.g., 21600 elements → 2160 elements), making the cost of a single exploration simulation decrease from 500 to 50.

[0081] 3. High-order spectral element discretizer: The output GLL / GLJ node distribution and Jacobian matrix have an exponentially improved accuracy: the error is <0.003 when N=4, which is 3 times more accurate than COMSOL, avoiding resource misjudgment caused by wave field distortion on exploration vessels.

[0082] Fluid-structure interaction solver: Input material parameters (density / velocity of sound), output pressure / displacement field with full wave field fidelity: accurately simulate Scholte interface waves and backscattering, solve the problem that traditional frequency domain methods cannot capture interface waves, and improve the underwater target recognition rate.

[0083] Strong form of elastic wave equation

[0084]

[0085] in: For density, For stress tensor, For external force source, This represents displacement acceleration.

[0086] For fluid regions, the stress tensor in the fluid contains only hydrostatic pressure. There is no shear component, meaning the stress tensor can be simplified to: Tensor divergence Substituting into (1), we get:

[0087]

[0088] To simplify calculations and eliminate the complexity of displacement vector fields, a scalar potential is introduced. Defined as: .

[0089] vector displacement field The solution is transformed into a scalar potential. The solution brings the following advantages: (1) significantly reduces the degree of freedom; (2) the scalar equation is easier to discretize and can be directly adapted to the Newmark explicit time integration algorithm; (3) the scalar potential form facilitates the coupling of the displacement field between the fluid and the solid region.

[0090] Substituting into equation (2) yields

[0091]

[0092] This indicates It is an irrotational field. Scalar potentials are permissible. A complete description of the displacement field. Also includes the constitutive relations of the fluid. ,in Let be the bulk modulus of the fluid. Combining equation (3), we get:

[0093]

[0094] use In the substitution formula (4) This leads to the partial wave equation of the fluid.

[0095]

[0096] exist Add point source

[0097]

[0098] Multiply both sides of the wave equation by the test function. and in the fluid domain By integrating, we can obtain the weak form.

[0099]

[0100] Applying the integral by parts formula to the divergence term on the right:

[0101]

[0102] Where: the first term is the stiffness integral, and the second term is the boundary integral. It is the unit vector normal to the boundary.

[0103] Physical constraints and mathematical implementation of boundary conditions

[0104] (1) The upper surface of the fluid is a free surface with zero pressure. The boundary integral term disappears here.

[0105] (2) Fluid-circular body interface Following the principle of continuity of normal stress, it is mathematically represented as: ,in Let be the displacement velocity on the solid side. The boundary integral is retained as:

[0106] To accurately inject excitation at a single point within the computational domain, a point source term is constructed using the filtering property of the Dirac delta function: for any continuous field f(x), the integral relation is satisfied.

[0107]

[0108] Combined equations (7)-(9),

[0109]

[0110] The momentum implications of each term in equation (10) are analyzed as follows:

[0111] Quality Items : Corresponds to the time rate of change of the system's kinetic energy, characterizing the transient response features dominated by inertial effects;

[0112] Stiffness term : Describes the gradient characteristics of potential energy, reflecting the constraint law of the elastic constitutive properties of the medium on the field distribution;

[0113] Coupling terms Characterizing momentum exchange at the fluid-solid interface, it is the core carrier of cross-medium dynamic coupling;

[0114] Source Item Point source energy injection provides initial excitation or external load drive for the wave field.

[0115] Weak-form spatial discretization: the weak-form integral (10) of the fluid region is extracted from physical space. Transition to reference unit The process of "coordinate transformation - integral term decomposition - boundary mapping" achieves the universality and precision control of numerical discretization.

[0116] Coordinate mapping: physical unit Through isoparametric mapping Convert to reference unit It relies on shape functions to achieve unified interpolation of geometry and field quantities.

[0117] During the mapping process, the Jacobian determinant Describe the area scaling relationship between the physical unit and the reference unit:

[0118]

[0119] Integral reduction correction in axisymmetric scenarios: In cylindrical coordinates, the integral volume element is... However, the weak form already implicitly contains axisymmetry (integral around...). (Symmetrical), therefore simplified to .

[0120] The weak form (10) of the fluid region is decomposed and transformed according to mass, stiffness, and boundary, following the logic of "integrand interpolation + Jacobian scaling + axisymmetric correction":

[0121] Conversion of mass integral terms:

[0122]

[0123] Transformation of stiffness integral terms:

[0124]

[0125] Transformation of boundary integral terms:

[0126]

[0127] When the computational unit approximates the unit along the axis of symmetry, the integral will exhibit... The singularity of the wavefield is a concern. Traditional Gauss-Lobatto-Legendre (GLL) integrals are constrained by the weighting function and node distribution characteristics, easily exhibiting oscillations and divergences in the paraxial region, leading to instability in wavefield simulations. To address this issue, the Gauss-Lobatto-Jacobi (GLJ) integral rule is introduced. By constructing modified Legendre polynomials and adapted basis functions, robust suppression of paraxial singularities can be achieved.

[0128] Revise the definition of the Legendre polynomial:

[0129]

[0130] In the formula and yes order and Legendre polynomial of order 1. Factoring the denominator. Ensure that When both the numerator and denominator approach zero, singularity can be avoided using L'Hôpital's rule. The modified polynomial... It exhibits smoother behavior in axisymmetric coordinate systems, making it suitable for integration of paraxial elements.

[0131] Based on the modified Legendre polynomial, the GLJ basis functions must satisfy a match between orthogonality and integral accuracy. hour

[0132]

[0133] The basis functions for interior points are adjusted by modifying the derivatives of the Legendre polynomials and the nodal spacing to ensure interpolation accuracy and orthogonality. In the denominator... Ensure that the basis functions are on the nodes At one node, the value is 1; at other nodes, the value is 0. (The molecule contains...) Terms and modified polynomials By combining derivatives, the integral stability in the paraxial region is optimized.

[0134] when The time base function degenerates to the standard value 1, which is compatible with the traditional GLL integral and ensures continuity between units. Finally, the GLJ integral shown in equation (17) is obtained.

[0135]

[0136] Introducing the GLJ integral brings numerical stability and compatibility: the GLJ integral avoids numerical divergence in the paraxial region by correcting the basis functions and integration points; in the off-axis region, the GLJ integral degenerates into the standard GLL integral, ensuring the high accuracy of the overall method.

[0137] In the reference unit The above uses GLJ and GLL nodes to access arbitrary functions. Expand:

[0138]

[0139] in: It is a function At GLJ and GLL nodes The value at that location, These are GLJ basis functions (used in the radial direction). ,deal with (the singularity) These are standard GLL basis functions (for axial directions). ).

[0140] Within the spectral element method discretization framework of the wave equation (9), after completing the intra-element GLL / GLJ integration and field interpolation, the displacement field continuity of the multi-element system needs to be achieved through cross-element coupling assembly. The core technical step lies in establishing the node mapping between the local and global meshes and superimposing the mechanical contributions of shared nodes, similar to the "system assembly" process of traditional finite element methods, thereby obtaining the global algebraic equation system:

[0141]

[0142] in The value of the collection potential, Collect all internal forces at the GLL / GLJ points. The term represents the mass matrix. The time discretization uses the explicit second-order precision Newmark scheme, which is conditionally stable; its Courant-Friedrichs-Lewy (CFL) stability condition depends on the ratio of the maximum P-wave velocity to the minimum grid spacing. The main computational cost stems from small-scale local matrix-vector multiplications between the element-level local field variables and the local stiffness matrix.

[0143] From equation (1), assuming there is no external force, then Select any vector to test the function Dot product with the equation and in the solid region Points accumulation:

[0144]

[0145] Using identities Substituting into the integral, we get:

[0146]

[0147] According to the divergence theorem, the first term can be transformed into a boundary integral:

[0148]

[0149] in It is the boundary of the solid region. It is the outward normal unit vector;

[0150] The solid boundary is divided into two parts: the free surface. :satisfy Therefore, the integral contribution is zero; fluid-solid coupling interface Following the normal stress continuity condition, the boundary integral term must be retained to describe the transmedium momentum transfer. Therefore, equation (22) degenerates into:

[0151]

[0152] Combining the volume inertia term and the constitutive term, the weak form of the solid domain is:

[0153]

[0154] In an axisymmetric coordinate system Under, stress tensor and strain tensor It exhibits symmetrical characteristics:

[0155]

[0156] Expanded in cylindrical coordinates:

[0157]

[0158] The partial expansion of strain energy in an elastic body must take into account both geometric nonlinearity and constitutive anisotropy. Based on the principles of continuum mechanics, the total strain energy can be decoupled into four-dimensional energy components: radial, axial, shear, and circumferential.

[0159] radial strain energy It characterizes the energy interaction between radial displacement gradient and radial normal stress, reflecting the energy storage of the medium during radial compression and tensile deformation.

[0160] Axial strain energy It describes the work density integral of axial displacement gradient and axial normal stress, and dominates the longitudinal bearing characteristics of columnar structures.

[0161] Shear strain energy It originates from the product of radial-axial shear strain and shear stress, reflecting the in-plane shear energy storage caused by non-uniform deformation.

[0162] Circular strain energy A term unique to axisymmetric systems, reflecting circumferential stress. For radial displacement The contribution comes from the definition of circumferential strain in cylindrical coordinates. This section passes the test function. The coupling with circumferential stress reflects the influence of geometric curvature.

[0163] Following similar steps to those used for the fluid component of the model, a similar set of global equations is obtained after assembling the system.

[0164]

[0165] in Represents all values ​​of displacement. This represents all internal forces at each GLL / GLJ point.

[0166] OpenMP Parallel Computation Unit: Calls 4 threads in parallel (omp_set_num_threads(4))

[0167] Significant improvement in computational efficiency: With 4 threads, CPU time decreased from 345s to 74.9s, a 4.6x speedup for a single simulation, supporting real-time decision-making for exploration vessels. While MPI parallel computing can achieve cross-node process-level parallelism, developers need to manually perform task partitioning, data distribution, and synchronization operations, significantly increasing the complexity of program development and debugging, and requiring high refactoring capabilities for serial programs. In contrast, OpenMP single-threaded parallel mode has less intrusiveness to the source code; therefore, this invention adopts the OpenMP parallel architecture. Its core design includes:

[0168] (1) Task decomposition and data distribution strategy:

[0169] A spatial domain-based partitioning method is adopted, dividing the computing domain into several sub-regions and allocating them to various threads. Automatic task allocation is achieved through the OpenMP parallel constructor instruction (#pragma omp parallel for), and load balancing is optimized by combining the schedule clause (dynamic scheduling or static block scheduling).

[0170] (2) Parallel execution model and data localization:

[0171] Based on the Fork-Join model, after the main thread creates a thread pool, all threads collaboratively execute local computational tasks within a parallel region. Efficient memory access management is achieved through data privatization and sharing mechanisms, reducing cross-thread data contention. For global communication needs such as coupled interfaces, critical sections or atomic operations are used to ensure data consistency.

[0172] (3) Dynamic load balancing and performance optimization:

[0173] To address scenarios with uneven computational loads, a runtime dynamic task queue (task + taskwait instructions) is introduced, employing a work-stealing mechanism to achieve adaptive load balancing among threads. Simultaneously, OpenMP's thread binding (OMP_PROC_BIND) feature is used to bind threads to physical cores, maximizing cache locality and computational efficiency.

[0174] (4) Global data aggregation and synchronization mechanism:

[0175] After the parallel computing phase is completed, global synchronization is achieved through the `omp_critical` or `omp_barrier` instructions to ensure that each thread completes its local data updates. For operations requiring global reduction (such as residual calculation or convergence judgment), the `omp_reduce` instruction is used to achieve efficient aggregation, ultimately forming a complete parallel computing process.

[0176] Through the above design, the framework fully leverages the computational potential of multi-core architecture while maintaining code maintainability, providing reliable technical support for high-performance solutions to complex fluid-structure interaction problems.

[0177] The mathematical properties of the spectral element method provide a natural advantage for parallel computing. Through non-overlapping subdomain partitioning, the computational domain is discretized into independent spectral element units. The stiffness matrix calculation of each element depends only on its internal GLL or GLJ sampling points, and the mass matrix within the element has a diagonal form, allowing the inversion of the mass matrix to be performed independently within the element without relying on a global linear system solution. This characteristic fundamentally eliminates the computational bottleneck of global matrix assembly and iterative solution in traditional finite element methods.

[0178] By employing a second-order explicit Newmark time integral, the time step depends only on the displacement and velocity fields of the previous moment. Each process independently advances its local time step. This local update mechanism significantly reduces the global dependency during the time stepping process, providing a natural parallel granularity for parallel computing.

[0179] The continuity constraint of the physical field only needs to be established at the boundaries of adjacent spectral elements. The communication volume between OpenMP threads is linearly proportional to the boundary length, which is far lower than the global communication overhead required by traditional finite difference methods or global explicit methods (such as the time-domain finite element method). Through the mechanism of "local computation + domain decomposition + explicit advancement", the complex problem of three-dimensional full-wave simulation is decomposed into parallel independent unit tasks. Furthermore, at the engineering level, it achieves improved computational efficiency (the diagonal mass matrix eliminates the bottleneck of solving global linear systems), enhanced scalability, and complexity control (avoiding exponential complexity growth), which is especially suitable for long-distance, high-resolution simulation scenarios such as marine acoustics.

[0180] The iterative solution process for fluid-structure interaction based on the OpenMP parallel architecture includes mesh generation for both fluid and solid domains, setting material properties and boundary conditions, initializing the number of OpenMP threads (omp_set_num_threads(N)); fluid solution: parallel calculation of the flow field (pressure, velocity); data transfer: transferring fluid pressure / shear force to the solid domain (coupling interface); solid solution: parallel calculation of structural displacement / stress; convergence check: checking whether the residuals at the coupling interface meet the convergence conditions.

[0181] 4. PML Boundary Manager: Automatically embeds the PML layer to absorb outgoing waves. Infinite domain simulation: Eliminates boundary reflection interference and avoids the 15% detection error caused by reflected waves in traditional methods.

[0182] 5. Sound field output module: Output receiver transmission loss curve (TL) directly supports decision-making: Generate TL curve at 30m underwater, providing quantitative basis for sonar system optimization (error <1dB).

[0183] The system relies on high-performance computing hardware to ensure the needs of large-scale data processing and high-precision computing: The computer CPU of this invention is an Intel Xeon Silver 4314, the memory is 32GB DDR4, the memory required for large-scale data processing and high-precision computing, the storage is 2TB NVMe SSD for fast data read and write, accelerating data loading and result storage, the network is InfiniBand HDR 200G, high-speed node communication, and supports large-scale parallel computing.

[0184] This invention constructs a homogeneous flat-bottomed FSI as a benchmark test case. The two-dimensional profile of this model is shown below. Figure 1 As shown, the structure comprises an upper 2400-meter-deep fluid medium and a lower 2400-meter-thick homogeneous viscoelastic solid medium, with an overall horizontal width of 6400m. Detailed material parameters are listed in Table 1. The Poisson's ratio for the fluid medium is 0.5, and for the solid medium, it is 0.25. Based on this two-dimensional model, a three-dimensional computational model with axisymmetric characteristics is constructed by expanding circumferentially along the right vertical boundary line as the axis of rotational symmetry. This modeling method effectively reduces the degrees of freedom in the three-dimensional simulation while preserving the original structural information, thus improving computational efficiency. The pressure source is located inside the seawater layer, 500 meters above the water surface. The source time function uses a Ricker wavelet with a dominant frequency of 10Hz. Receiver points are arranged horizontally at a distance of 2500 meters from the source, and vertically at depths of 2300 meters and 3300 meters above the water surface, respectively. This arrangement aims to cover the wavefield response at different depths, thereby comprehensively monitoring the energy attenuation and wave characteristics as depth changes during three-dimensional sound wave propagation. A perfectly matched layer (PML) was placed around the perimeter and bottom of the model to absorb the energy of the outwardly propagating wave field. The polynomial degree was N=4, and the total number of spectral elements was 21600. The time step was set to 0.6 ms, with a total of 4500 time steps and a simulation duration of 2.7 s. The effectiveness of the design was verified using a homogeneous flat-bottom case.

[0185] Table 1 Material parameters for homogeneous flat-bottomed FSI case

[0186] fluid 1000 1500 -- solid 2500 3400 1963

[0187] Figure 2 The results show a comparison of displacement component results at receiver point 2 within the range of 1.2s to 2.2s, obtained using the DR-SETD solver and the COMSOL solver. Figure 2 (a) in the figure represents the horizontal component. Figure 2 The vertical component (b) shows that the relative error between the two solvers is less than 0.003. This verification demonstrates that DR-SETD can accurately simulate wave propagation in fluid-solid coupled media.

[0188] The DR-SETD solver has 522,222 DoFs, while the COMSOL solver has 531,746 DoFs. The COMSOL solver's CPU time consumption is 345 seconds, while the four-thread accelerated DR-SETD solver's CPU time consumption is 74.9 seconds, significantly lower than the COMSOL solver, with a computational efficiency approximately 4.6 times that of the former. This indicates that the DR-SETD solver has a significant efficiency advantage, particularly due to its use of OpenMP parallel acceleration, which improves the computational efficiency of the element mass and stiffness matrices.

[0189] Table 2 Speedup and Parallel Efficiency

[0190] Thread count CPU time (s) acceleration ratio Parallel efficiency 1 215.5 1.0 100% 2 127.4 1.69 84.5% 4 74.9 2.88 72% 8 61.9 3.48 43.5%

[0191] The DR-SETD solver exhibits excellent scalability (efficiency >72%) in the 1-4 thread range, demonstrating the high parallelism of its computational core. With 8 threads, the speedup reaches 3.48 (efficiency 43.5%), with the decrease primarily limited by memory bandwidth saturation and non-communication overhead. Doubling the number of threads means an increase in the "boundaries" requiring data exchange or points requiring synchronization. Even if the amount of data per thread decreases, the overall communication frequency or management overhead increases. Current results indicate that DR-SETD reaches its optimal cost-performance ratio with a 4-physical-thread configuration.

[0192] Figure 4 The presentation displays instantaneous wavefield snapshots at four key time points: 0.3s, 0.9s, 1.08s, and 1.5s. Red and blue contour lines characterize the positive and negative distribution of displacement / pressure disturbances, illustrating the propagation process of elastic waves in a fluid (upper blue)-solid (lower gray) coupled system. It covers four core wave types: direct waves, free surface reflected waves, interface reflected / converted volume waves, and Scholte interface waves. Four sub-plots (ad) correspond to the time series of wave evolution, showcasing the physical chain of "direct wave dominance → interface interaction → volume wave separation → interface wave dominance."

[0193] Figure 4 In (a): The direct longitudinal wave (P-wave) of the fluid is centered on the hypocenter. Within the fluid layer, a spherical wavefront (alternating red and blue rings) is formed. The wavefront morphology is pure and unaffected by the free surface (upper boundary of the fluid) or the solid-liquid interface (lower boundary of the fluid). The wavefront propagation speed matches the fluid sound speed (≈1500 m / s), and the amplitude decreases with increasing radial distance, perfectly matching the free propagation law of elastic waves in a homogeneous fluid medium.

[0194] Figure 4In (b, c), the wavefront propagation rate in the solid region is significantly faster than in the fluid region. This is entirely consistent with the theoretical expectation that "solid media, due to their shear modulus, can support transverse waves (S-waves), and longitudinal waves (P-waves) have higher elastic stiffness, resulting in a volume wave velocity far greater than that of fluid media (which have no shear modulus) that can only propagate longitudinal waves," directly confirming the rationality of the medium parameter settings. No reflected echoes were observed at the right boundary of the solid region (the red and blue stripes did not rebound or oscillate at the boundary). This indicates that the perfectly matched layer (PML), as an absorbing boundary condition, has achieved efficient non-reflective absorption of the outgoing wave, successfully simulating the wave radiation characteristics of an "infinite domain," and verifying the stability and accuracy of the numerical algorithm in open boundary problems. Near the solid-liquid coupling interface, a narrow-band energy-concentrated red and blue stripe structure can be observed—this is the early excitation mode of Scholte interface waves.

[0195] Figure 4 In (d): At the solid-liquid coupling interface, the narrow-band continuous red and blue stripes clearly present the mature form of the Scholte interface wave: its energy is strictly confined to the solid-liquid interface region, with no significant energy diffusion above the fluid layer and below the solid layer, which is highly consistent with the theoretical model that "Scholte waves act as interface-constrained guiding waves, and their energy decays exponentially along the interface normal"; by comparing the propagation characteristics of the rapidly propagating volume wave wavefront in the solid layer and the reflected wave in the fluid layer, it can be seen that the propagation rate of the Scholte wave is significantly lower, directly verifying the dynamic law that "Scholte wave speed is lower than that of fluid longitudinal waves and solid volume waves (P / S waves)"; the stripes at the interface maintain continuous and regular phase alternation characteristics, without scattering or energy distortion, fully demonstrating the high-precision characterization capability (high-order accuracy and numerical stability) of the DR-SETD solver for strongly coupled boundary conditions of multiple media. Meanwhile, the non-reflective absorption mechanism of PML was directly verified above the fluid layer: after the upward wave propagating in the fluid layer reached the upper boundary, no wave rebounded, indicating that PML achieved "full absorption" of the outgoing wave and successfully simulated the open boundary conditions of the "infinite fluid domain"; there were no false reflections or numerical oscillations in the PML region, proving that the absorption layer parameters (damping gradient, thickness) were reasonably designed, avoiding the contamination of the core wave field by boundary reflections, and providing a reliable numerical environment for long-term wave simulation.

[0196] To verify the applicability of DR-SETD in handling complex terrain and backscattering scenarios, this invention designs and constructs an axisymmetric inclined seabed FSI model. This model is based on... Figure 5The two-dimensional cross-sectional structure shown comprises an upper fluid medium and a lower homogeneous viscoelastic solid medium, with an overall horizontal width of 10 km and a vertical width of 3 km. Specifically, the fluid depth exhibits a segmented gradient variation with horizontal position. Using the right boundary line as a reference, the water depth is constant at 600 m in the 0-2 km range, while the solid slope in the 2-6 km range slopes towards the sea surface at a ratio of 7.1:1, forming a continuously gradually changing slope topography. The water depth in the 2-6 km range is constant at 100 m. Detailed material parameters are shown in Table 3.

[0197] Based on this two-dimensional model, a three-dimensional computational domain with axisymmetric characteristics is constructed by circumferentially expanding the right vertical boundary line as the axis of rotational symmetry. The acoustic excitation source is set at 590 meters above the water surface, using a Ricker wavelet with a dominant frequency of 5 Hz. The receiving array is fixed vertically at 30 meters underwater and horizontally distributed at equal intervals of 10 meters from the right boundary outwards, covering a range of 0 to 10 km, with a total of 1000 observation points.

[0198] By absorbing boundary wavelength energy through PML, with a total of 31,349 spectral units, a polynomial order of N=4, a time step of 1.22ms, and 20,000 simulation steps totaling 24.4s, this configuration ensures computational efficiency while accurately capturing the interface and backscattering effects caused by slope topography, providing reliable data support for sound field modeling and propagation loss analysis in complex marine environments.

[0199] Table 3 Material parameters for inclined seabed FSI case studies

[0200] fluid 1000 1500 -- solid 2000 3400 1200

[0201] Figure 6 The numerical solution calculated by DR-SETD and the frequency domain solution simulated by COMSOL are compared and verified: the transmission loss curves at a center frequency of 5 Hz are compared and analyzed at a depth of 30m underwater. Figure 6 In the figures (a) and (b), the transmission loss curves of DR-SETD (blue dashed line) and COMSOL (red solid line) in different propagation distance ranges are shown to verify the consistency of their morphology and quantization error characteristics. Their oscillation trends and extreme values ​​(corresponding to the enhancement / attenuation nodes of acoustic interference during propagation) are completely coincident. Figure 6 In (a), multiple sets of resonance peaks and valleys in the 0.5–3.5 km range are strictly aligned in both amplitude and phase dimensions, which intuitively demonstrates the ability of the two methods to effectively characterize nonlinear physical processes such as acoustic multipath interference and terrain backscattering. Figure 6In (b) of the diagram, the low-frequency oscillation characteristics of the curves in the 5–10 km long propagation segment (such as the TL peak at 6 km and the fluctuation trend in the 8–9 km range) are highly consistent, demonstrating the algorithm's accurate capture of long-range wave field attenuation and scattering coupling effects. The error of the TL curves across the entire field is less than 1 dB, verifying the accuracy of DR-SETD in simulating the propagation of complex topographic waves with backscattering, and making it particularly suitable for scenarios where traditional approximation methods fail.

[0202] The degrees of freedom (DoF) of the two solver models are similar: DR-SETD solver has 980,666 DoF, and COMSOL solver has 971,746. The CPU time consumption of the COMSOL solver is 1182s, while the CPU time consumption of the eight-thread accelerated DR-SETD solver is 236.2s, significantly lower than the former, and its computational efficiency is 5 times that of COMSOL.

[0203] The shared-memory parallel performance of the DR-SETD solver was further tested.

[0204] Table 4 Speedup and Parallel Efficiency

[0205] Thread count CPU time (s) acceleration ratio Parallel efficiency 1 758.8 1.0 100% 2 580.1 1.30 65% 4 370.7 2.05 51% 8 236.2 3.21 40%

[0206] The efficiency of dual-threaded and quad-threaded solutions is 65% and 51% respectively, which is typical of memory-intensive tasks. The eight-threaded solution achieves a speedup of 3.21 (40% efficiency), and memory bandwidth becomes the main bottleneck for shared memory expansion. Current results indicate that DR-SETD has reached its optimal cost-performance ratio with an 8-physical-thread configuration. Future work will deploy a hybrid parallel scheme of MPI+OpenMP based on a distributed memory system to address the needs of solving larger-scale problems.

[0207] Figure 7 The (ad) section shows instantaneous wave field snapshots of the inclined seabed FSI case at three key time points: 0.49s, 1.46s, 2.32s, and 7.32s, systematically revealing its propagation process and dynamic behavior.

[0208] Figure 7 In (a) and (b): the absence of significant reflected echoes at the top boundary indicates that the PML absorbing boundary has taken effect, preventing boundary reflections from interfering with the core scattering process and verifying the stability of the algorithm in open boundary scenarios. At this point, the waveform has not yet entered the sloping terrain, and the wavefield analysis is similar to the flat-bottom FSI described earlier.

[0209] Figure 7(c) In this case, the fluid-solid interface exhibits a sloping topography with a higher elevation on the left and a lower elevation on the right. The wavefield excited by the seismic source quickly adapts to the topographic morphology. Within the fluid layer, the wavefront propagates along the sloping solid-liquid interface, and the curvature and density gradient of the red and blue stripes visually reflect the degree of energy concentration. Within the solid layer, the wavefront is constrained by the topography, exhibiting a dual-mode propagation of "spreading along the sloping interface + radiating into the deep solid domain," confirming the geometric adaptation capability of DR-SETD to non-horizontal boundaries. The wavefront morphology adaptively adjusts with the sloping solid-liquid interface (smaller wavefront curvature in the fluid domain and larger wavefront curvature in the solid domain), visually demonstrating DR-SETD's geometric adaptation capability to non-horizontal boundaries (sloping seabed)—high-order spectral element discretization accurately fits the topographic undulations, ensuring the physical authenticity of the wavefront morphology and energy distribution.

[0210] Figure 7 (d) shows the energy decay characteristics after long-term evolution. The red and blue stripes of the core wave field almost disappear, indicating that most of the energy has been dissipated through the PML absorption edge or the natural decay of the wave field. It successfully simulates the radiation characteristics of infinite domain waves and verifies the reliability of the algorithm in long-term open domain simulation. There are no false reflections or numerical oscillations in the solid layer, which proves the numerical stability of DR-SETD in complex terrain + open boundary scenarios. The coupling of high-order time integral and spectral element spatial discretization not only ensures the physical reality of wave field evolution, but also avoids the "numerical dispersion / oscillation" problem that is prone to occur in low-order methods.

[0211] Based on this, the present invention provides a parallel simulation method for ocean acoustic fields based on DR-SETD. Through collaborative innovation in multiple dimensions such as theoretical models, spatial discretization, time integration, and parallel implementation, it constructs a fluid-structure interaction acoustic field simulation tool with high accuracy, high efficiency, and high robustness. Its core advantages are reflected in the following aspects:

[0212] 1. Revolutionary computational efficiency:

[0213] Fundamental dimensionality reduction: Through ingenious processing of the three-dimensional axisymmetric model, a cross-order-of-magnitude compression of the computational scale is achieved, fundamentally solving the bottleneck of the "curse of dimensionality" in traditional three-dimensional simulation.

[0214] High-performance parallel architecture: Employing fine-grained parallel strategies and dynamic load balancing technology, it fully leverages the performance of modern multi-core processors, transforming computing resources into tangible solution speed.

[0215] 2. Excellent numerical accuracy:

[0216] Exponential convergence: Integrating the essence of the spectral method, the accuracy of the solution increases exponentially with the element order, rather than algebraically with the traditional finite element method, thus achieving extremely high numerical accuracy with fewer degrees of freedom.

[0217] Full-wave field physical fidelity: As a time-domain method, it can fully capture all wave phenomena (such as reflection, refraction, scattering and complex interface waves), avoiding the inherent shortcomings of frequency-domain methods or parabolic approximation methods at the physical level.

[0218] 3. Strong engineering applicability and robustness:

[0219] Flexible geometric adaptation capability: Inheriting the advantages of the finite element method in handling complex regions, it can accurately fit various irregular boundaries (such as inclined seabed and curved surface structures), and get rid of the dependence on regular meshes.

[0220] Unconditionally stable time progression: The Newmark implicit integral scheme is adopted to ensure absolute numerical stability in long-term simulations and overcome the step size limitation caused by stability conditions in explicit methods.

[0221] 4. Paradigm transcendence in comprehensive performance:

[0222] Collaborative optimization system engineering: This method is not an improvement of a single technology, but rather an improvement in "accuracy, efficiency and stability" through the linkage of multiple technologies such as "dimensionality reduction + high-order discretization + implicit integration + fine parallelism".

[0223] Maturity benchmarked against industry leaders: After rigorous comparison and verification with authoritative commercial software COMSOL, this method demonstrates superior solution speed and numerical robustness while maintaining equal or even higher accuracy (error <1dB), and has the potential to replace traditional solutions.

[0224] In summary, the DR-SETD method successfully combines cutting-edge theories of computational mechanics with parallel techniques of modern computer science, providing an advanced solution that comprehensively surpasses existing mainstream tools in terms of accuracy, efficiency, and reliability for solving highly complex fluid-structure interaction problems in fields such as marine acoustics and aeroelasticity.

[0225] Example 2

[0226] In this embodiment, a computer terminal device is provided, including:

[0227] One or more processors;

[0228] A memory, coupled to the processor, for storing one or more programs;

[0229] When the one or more programs are executed by the one or more processors, the one or more processors implement the steps of the above-described parallel simulation method for ocean acoustic fields based on DR-SETD.

[0230] In this embodiment, a computer-readable storage medium is also provided, on which a computer program is stored. When the computer program is executed by a processor, it implements the steps of the above-described parallel simulation method for ocean acoustic fields based on DR-SETD.

[0231] In this embodiment, an electronic device is also provided, including a memory and a processor. The memory stores a computer program, and the processor is configured to run the computer program to perform the steps of the above-described parallel simulation method for ocean acoustic fields based on DR-SETD.

[0232] In this embodiment, a computer program product is also provided, including a computer program that, when executed by a processor, implements the steps of the above-described parallel simulation method for ocean acoustic fields based on DR-SETD.

[0233] The aforementioned program can run on a processor or be stored in memory (or a computer-readable medium). Computer-readable media includes both permanent and non-permanent, removable and non-removable media, and information storage can be achieved by any method or technology. Information can be computer-readable instructions, data structures, program modules, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random-access memory (SRAM), dynamic random-access memory (DRAM), other types of random-access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other non-transfer medium that can be used to store information accessible by a computing device.

[0234] These computer programs may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes can be implemented by different modules for different steps.

[0235] This embodiment provides such a device or system. The system, referred to as a DR-SETD-based parallel simulation system for ocean acoustic fields, includes:

[0236] The axisymmetric modeling module is used to acquire the physical parameters and topographic data of the target sea area, and to simplify the three-dimensional fluid-structure coupled acoustic field model based on the axisymmetric assumption to obtain an axisymmetric calculation model.

[0237] The spectral element discretization module, connected to the axisymmetric modeling module, is used to configure integration nodes on the units of the axisymmetric computation model using a higher-order spectral element method based on orthogonal polynomials to complete the spatial discretization of the computational domain.

[0238] The coupled solution module, connected to the spectral element discretization module, is used to construct weakly coupled forms of the time-domain wave equation for both fluid and solid media, and assemble them to obtain the global system equation for solution.

[0239] The parallel computing module, connected to the coupled solution module, is used to perform parallel processing and data synchronization of the computational tasks of each discrete unit based on shared memory during the solution process using an explicit time integration scheme.

[0240] As one implementation method in this embodiment, the axisymmetric modeling module includes:

[0241] The topographic data acquisition unit is used to input two-dimensional topographic profile data of the target sea area, including changes in seawater depth and seabed interface morphology.

[0242] The model generation unit, connected to the terrain data acquisition unit, is used to rotate the two-dimensional terrain profile around one of its vertical sides to generate a three-dimensional axisymmetric computational domain.

[0243] As one implementation method in this embodiment, the spectral element discretization module includes:

[0244] A node configuration unit is used to configure first-kind Gauss-Lobart integral nodes in the axial direction of the unit and second-kind Gauss-Lobart integral nodes based on modified Legendre polynomials in the radial direction.

[0245] The local matrix calculation unit, connected to the node configuration unit, is used to perform high-order interpolation of the element field variables through the two types of nodes, and to calculate the element's local stiffness matrix and diagonal mass matrix.

[0246] As one implementation method in this embodiment, the coupled solution module includes:

[0247] The fluid equation building unit is used to introduce a scalar velocity potential into the fluid domain and derive the weak form of the wave equation with the potential function as the variable.

[0248] Solid equation building unit, used to derive the weak form of linear wave equations for solid domains with displacement vectors as variables;

[0249] An interface coupling unit connects the fluid equation construction unit and the solid equation construction unit, and is used to introduce boundary integral terms at the fluid-solid interface to achieve the continuity condition of stress and displacement.

[0250] As one implementation method in this embodiment, the parallel computing module includes:

[0251] The domain decomposition unit is used to divide all discretized spectral units into multiple sub-regions according to the spatial domain decomposition method.

[0252] The task scheduling unit, connected to the region decomposition unit, is used to dynamically allocate computing tasks within each sub-region to multiple threads for execution using OpenMP parallel instructions.

[0253] As one implementation method in this embodiment, the coupled solution module further includes:

[0254] Explicit time-progression unit for updating field variables using a second-order accurate Newmark-β explicit integration scheme;

[0255] The parallel computing module is configured to enable each thread to independently advance the computation by calling the explicit time advancement unit and using the local diagonal mass matrix within its allocated sub-region.

[0256] The system or apparatus is used to implement the functions of the methods in the above embodiments. Each module in the system or apparatus corresponds to each step in the method, as has been described in the method and will not be repeated here.

[0257] The above implementation method solves the problem of parallel simulation of marine acoustic field based on DR-SETD in related technologies, thereby ensuring that the problems existing in the prior art are resolved.

[0258] The above are merely preferred embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A parallel simulation method for ocean acoustic fields based on DR-SETD, characterized in that, Includes the following steps: The physical parameters and topographic data of the target sea area are obtained. Based on the axisymmetric assumption, the three-dimensional fluid-structure coupled acoustic field model is reduced in dimension to obtain an axisymmetric calculation model. The higher-order spectral element method based on orthogonal polynomials is adopted, and two different types of integration nodes are configured on the unit of the axisymmetric computation model to complete the spatial discretization of the computation domain. We construct and derive the weak forms of the time-domain wave equations that satisfy the stress and displacement continuity conditions for fluid and solid media respectively, and assemble them to obtain the global coupled system equations. The global coupled system equations are solved by time-progression using an explicit time integration scheme. Within each time step, the computational tasks of each discrete unit are processed in parallel and data is synchronized using a shared memory parallel computing architecture.

2. The method according to claim 1, characterized in that, The process of reducing the dimension of the three-dimensional fluid-structure interaction acoustic field model to obtain an axisymmetric computational model includes: Input two-dimensional topographic profile data of the target sea area, which includes changes in seawater depth and seabed interface morphology; The two-dimensional terrain profile is rotated about one of its vertical sides to generate a three-dimensional axisymmetric computational domain, thereby simplifying the original three-dimensional problem into a computational problem within the two-dimensional profile.

3. The method according to claim 1, characterized in that, The process of spatial discretization using the higher-order spectral element method based on orthogonal polynomials includes: First-type Gauss-Lobart integration nodes are configured in the axial direction of the element; Second-type Gaussian-Lobart integral nodes are configured in the radial direction of the element. The second-type nodes are constructed based on modified Legendre polynomials. Higher-order interpolation of field variables within the element is performed using two types of nodes, and the corresponding local stiffness matrix and diagonal mass matrix are calculated.

4. The method according to claim 1, characterized in that, The process of constructing and deriving the weak form of the time-domain wave equation includes: For the fluid domain, a scalar velocity potential is introduced, and the displacement field is expressed as the gradient of this potential function, thus deriving the weak form of the wave equation with the potential function as the variable. For the solid domain, the weak form of the elastic wave equation is derived directly using the displacement vector as the variable, based on the constitutive relations of elasticity. At the shared interface between the fluid and solid regions, a boundary integral term is introduced to achieve forced coupling of the continuity conditions of normal stress and normal displacement.

5. The method according to claim 1, characterized in that, The process of parallel processing and data synchronization using a shared-memory parallel computing architecture includes: The spatial domain decomposition method is used to divide all discretized spectral units into sub-regions whose number matches the number of parallel threads; By using OpenMP parallel compilation guidelines, the cyclic computation tasks of each cell within a sub-region are dynamically scheduled and allocated to multiple threads for simultaneous execution. At the coupling interface, a thread synchronization mechanism is used to ensure data consistency between adjacent units computed by different threads.

6. The method according to claim 5, characterized in that, The process of using an explicit time integration scheme for time advancement includes: The field variables at each time step are updated using a second-order precision Newmark-β explicit integral scheme; Within the allocated sub-region of a cell, each thread independently and rapidly solves the motion equations of the cell nodes using the local diagonal mass matrix, thus advancing the local time step.

7. A parallel simulation system for marine acoustic fields based on DR-SETD, characterized in that, The system for implementing the method of any one of claims 1-6 comprises: The axisymmetric modeling module is used to acquire the physical parameters and topographic data of the target sea area, and to simplify the three-dimensional fluid-structure coupled acoustic field model based on the axisymmetric assumption to obtain an axisymmetric calculation model. The spectral element discretization module, connected to the axisymmetric modeling module, is used to configure integration nodes on the units of the axisymmetric computation model using a higher-order spectral element method based on orthogonal polynomials to complete the spatial discretization of the computational domain. The coupled solution module, connected to the spectral element discretization module, is used to construct weakly coupled forms of the time-domain wave equation for both fluid and solid media, and assemble them to obtain the global system equation for solution. The parallel computing module, connected to the coupled solution module, is used to perform parallel processing and data synchronization of the computational tasks of each discrete unit based on shared memory during the solution process using an explicit time integration scheme.

8. A computer terminal device, characterized in that, include: One or more processors; A memory, coupled to the processor, for storing one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors perform the steps of the method as described in any one of claims 1-6.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method as described in any one of claims 1-6.

10. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1-6.

Images