Acoustic Field Acquisition Method of Two-Dimensional Acoustic Field Model Based on Finite Element Method

Abstract

The invention belongs to the technical field of underwater sound propagation, and in particular, relates to a sound field acquisition method of a two-dimensional sound field model based on a finite element method, comprising: discretizing a target area into a finite number of calculation units according to ocean environment parameters of the acquired target area; constructing a basis function of each calculation unit according to a unit shape and a number of nodes; using a standard Galerkin variation method to derive a weak form of the Helmholtz equation; integrating unit stiffness matrices of all calculation units into a total stiffness matrix according to a one-to-one correspondence between each local number of the node and each overall number; integrating unit external load vectors of all calculation units into a total external load vector according to a one-to-one correspondence between each local number of the node and each overall number; composing the above results into a linear algebraic equation group as a two-dimensional sound field model to obtain a sound field, that is, a sound pressure of each node; and obtaining a sound pressure value of each node in space by selecting a basis function and an interpolation method.

Description

Technical Field

[0001] The invention belongs to the technical field of underwater sound propagation in hydroacoustic physics and hydroacoustic engineering, and in particular relates to a sound field acquisition method of a two-dimensional sound field model based on a finite element method. Background Art

[0002] Sound waves are the only physical form known to date that can propagate over long distances in seawater. They are an important carrier for detecting the marine environment, developing marine resources, and realizing underwater information transmission. The study of the propagation characteristics of sound waves in seawater is one of the important areas of oceanographic research and is also the basis for studying other marine acoustic issues.

[0003] The commonly used sound field calculation models are difficult to accurately calculate the sound field in complex environments. For example, the ray method is only suitable for calculating the propagation of high-frequency sound waves; the parabolic equation method has poor accuracy for calculating the sound field in an environment with drastic level changes; the simple normal wave method and wave number integration method without special treatment cannot calculate the sound field in an environment with level changes; the coupled simple normal wave method cannot accurately calculate the sound field changes caused by asymmetric terrain.

[0004] Therefore, existing methods cannot accurately calculate the sound field in a two-dimensional complex ocean environment, and cannot solve the problem of underwater sound propagation in a two-dimensional complex environment. Summary of the invention

[0005] In order to solve the above defects in the prior art, the present invention proposes a sound field acquisition method of a two-dimensional sound field model based on a finite element method, which is a finite element method for calculating two-dimensional underwater sound propagation problems. The method includes:

[0006] According to the ocean environment parameters of the target area, the target area is discretized; the target area is discretized into a finite number of calculation units, and each calculation unit is connected by a node; according to the unit shape and the number of nodes, the basis function of each calculation unit is constructed;

[0007] Using the standard Galerkin variational method, the weak form of the Helmholtz equation is derived;

[0008] Establish the unit stiffness matrix of each calculation unit, and integrate the unit stiffness matrices of all calculation units into the overall stiffness matrix according to the one-to-one correspondence between the local numbers of the nodes and the overall numbers;

[0009] Establish the unit external load vector of each calculation unit, and integrate the unit external load vectors of all calculation units into the overall external load vector according to the one-to-one correspondence between each local number of the node and each overall number;

[0010] The overall stiffness matrix, the sound pressure vector and the overall external load vector are combined into a linear algebraic equation system as a two-dimensional sound field model, and the sound field, that is, the sound pressure at each node, is obtained by solving the two-dimensional sound field model;

[0011] Then, the sound pressure value at each point in the space is obtained through the selected basis function and interpolation method.

[0012] As one of the improvements of the above technical solution, the marine environment parameters include: seabed sound velocity, seabed density, seabed topography of the target area, sound source intensity and sound source position.

[0013] As one of the improvements of the above technical solution, according to the acquisition of the marine environment parameters of the target area, the target area is discretized into a finite number of calculation units, and each calculation unit is connected by nodes; according to the unit shape and the number of nodes, the basis function of each calculation unit is constructed; the specific process includes:

[0014] Using three-node triangle units as calculation units, the target area is divided according to the ocean environment parameters of the target area, and the target area is discretized into a finite number of three-node triangle units. The basis function of each calculation unit is N n (x,z),n=1,2,3; where n represents the local number of the node, (x,z) is the established rectangular coordinate system; the z axis passes through the sound source and is vertically downward; the x axis is parallel to the sea surface;

[0015] Assume that the three nodes are represented by A n (x n ,y n ), n = 1, 2, 3, and define the determinant of the Jacobi matrix |J| = (x2-x1)(y3-y1)-(x3-x1)(y2-y1), where x n is the horizontal coordinate of the nth node; y n is the ordinate of the nth node; the expressions of the three basis functions corresponding to the three nodes of the calculation unit are obtained:

[0016]

[0017]

[0018]

[0019] As one of the improvements of the above technical solution, the weak form of the Helmholtz equation is derived using the standard Galerkin variational method; the specific process includes:

[0020] For the infinite line source problem, a rectangular coordinate system (x, z) is established, with the z axis passing through the sound source and pointing vertically downward, and the x axis parallel to the sea surface; at this time, the Helmholtz equation in the inhomogeneous medium is expressed as

[0021]

[0022] Where ρ represents the seawater density and k(x,z) represents the wave number, which can be expressed as Where c(x,z) represents the speed of sound, f is the frequency of the sound source; S ω represents the intensity of the sound source, z s Indicates the depth of the sound source; is the gradient operator; p is the sound pressure; δ is the Dirac function;

[0023] The weak form of the governing equation is derived using the standard Galerkin variational method. The specific process is as follows:

[0024] Multiply both sides of the above Helmholtz equation by the test function q(x,z) and integrate it over the solution domain Ω of the target area to obtain the following expression:

[0025]

[0026] Applying Gauss's theorem to the above expression yields the weak form of the Helmholtz equation:

[0027]

[0028] Among them, Γ represents the boundary, Ω represents the solution domain of the target area;

[0029] The Lagrange interpolation function corresponding to the three-node triangle unit is used to discretize the sound pressure and test function of the target area to obtain the discrete sound pressure p(x,z) and the discrete test function q(x,z);

[0030]

[0031] q(x,z)=N m (x,z),m=1,2,3

[0032] The discretized Helmholtz equation without considering boundary conditions is derived:

[0033]

[0034] Among them, p n is the sound pressure value of the nth node to be determined; Yes n The coefficient is used as the element k of the unit stiffness matrix mn; n and m are local numbers of nodes, n = 1, 2, 3, m = 1, 2, 3; N n and N m are the basis functions corresponding to nodes n and m respectively; and N n and N m The gradient of

[0035] For the external load vector of the unit, the discretized test function q(x,z) is substituted into the right side of the weak form of the Helmholtz equation to obtain Take it as the mth element of the element external load vector.

[0036] As one of the improvements of the above technical solution, the unit stiffness matrix of each calculation unit is established, and the unit stiffness matrices of all calculation units are integrated into the overall stiffness matrix according to the one-to-one correspondence between each local number and each overall number of the node; the specific process includes:

[0037] Establish the unit stiffness matrix of each calculation unit, traverse all calculation units, and add each element of the unit stiffness matrix to each corresponding element of the overall stiffness matrix according to the correspondence between the local number of each node in the calculation unit and the elements in the overall stiffness matrix, and then integrate the elements of the unit stiffness matrix into the corresponding positions in the overall stiffness matrix.

[0038] As one of the improvements of the above technical solution, the overall stiffness matrix, the sound pressure vector and the overall external load vector are combined into a linear algebraic equation group as a two-dimensional sound field model, and the two-dimensional sound field model is solved to obtain the sound field; and then the sound pressure value of each node in the space is obtained by selecting the basis function and the interpolation method; the specific process includes:

[0039] For an infinite area, the sound field satisfies the radiation condition: when r→∞, the sound pressure p→0;

[0040] At this time, it is necessary to truncate the infinite area, use artificial boundary conditions to simulate the emission conditions, and select the perfect matching layer method to simulate the radiation conditions; the perfect matching layer is an absorption layer that is considered to be set to absorb the outgoing waves at various angles;

[0041] Use a perfectly matched layer to simulate an infinite region, perform coordinate transformation, and map the real coordinate axis to the complex coordinate axis. The specific process is as follows:

[0042] For a medium with limited depth, assume that the horizontal distance of the physical domain to be calculated is R, and set a thickness of δ on the right side of the physical domain. x The sound velocity and density in the perfectly matched layer are the same as those of the adjacent left medium;

[0043]

[0044] in, is the complex coordinate after coordinate transformation; i is the imaginary unit; ξ is the integral variable; σ x (ξ) is an artificially defined attenuation function, Among them, σ xm is a constant that determines the maximum attenuation of the perfectly matched layer;

[0045] For a perfectly matched layer, p n The coefficient is denoted as k mn , k mn are the elements of the element stiffness matrix in the perfectly matched layers:

[0046]

[0047] Using the elements of the unit stiffness matrix for the perfectly matched layer obtained above, an overall stiffness matrix is ​​established; then a linear algebraic equation group is formed by combining the overall external load vector and the sound pressure vector as a two-dimensional sound field model, and the two-dimensional sound field model is solved to obtain the sound field;

[0048] Kp=F

[0049] Among them, K represents the overall stiffness matrix, F represents the overall external load vector, and p represents the acoustic field;

[0050] Then, according to the obtained sound field, the sound pressure value of each node in the space is obtained by selecting the basis function and interpolation method.

[0051] The beneficial effects of the present invention compared with the prior art are:

[0052] 1. The method of the present invention has a wide range of applications. Compared with the ray method, it is applicable to high-frequency and low-frequency problems. Compared with the parabolic equation method, it is applicable to calculating the sound field in an ocean environment with drastic level changes. Compared with the simple normal wave and wave number integral method, it is applicable to calculating the sound field in a level change environment.

[0053] 2. The method of the present invention uses the finite element method to simulate underwater sound propagation, which has high accuracy and strong universality, and can solve the underwater sound propagation problem in a two-dimensional complex environment that cannot be solved by existing models; theoretically, the finite element method can obtain calculation results of ideal accuracy by continuously subdividing the grid;

[0054] 3. The method of the present invention is easy to parallelize. With the continuous improvement of computer performance and the continuous development of parallel technology, the calculation efficiency of the finite element model can be effectively improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 It is a method flow chart of a method for acquiring a sound field of a two-dimensional sound field model based on a finite element method of the present invention;

[0056] Figure 2a This is a schematic diagram of the propagation loss of a double-layer medium calculated using the wave number integration method when the sound source frequency is 100 Hz in a waveguide with limited depth;

[0057] Figure 2b It is a schematic diagram of the propagation loss of a double-layer medium calculated using a finite element model when the sound source frequency is 100 Hz in a waveguide with limited depth;

[0058] Figure 3 This is a schematic diagram of the propagation loss at a receiving depth of 60m when the sound source frequency is 100Hz;

[0059] Figure 4a It is the sound field propagation loss in the penetrable seafloor slope waveguide calculated using the simple normal wave model DGMCM when the sound source frequency is 30 Hz, the sound source depth is 20 m, and the horizontal distance is 50 m;

[0060] Figure 4b It is the sound field propagation loss in the penetrable seafloor slope waveguide calculated using the finite element model when the sound source frequency is 30 Hz, the sound source depth is 20 m, and the horizontal distance is 50 m;

[0061] Figure 4c It is the sound field propagation loss in the penetrable seafloor slope waveguide calculated using the simple normal wave model DGMCM when the sound source frequency is 50 Hz, the sound source depth is 20 m, and the horizontal distance is 50 m;

[0062] Figure 4d It is the sound field propagation loss in the penetrable seabed slope waveguide calculated using the finite element model when the sound source frequency is 50 Hz, the sound source depth is 20 m, and the horizontal distance is 50 m. DETAILED DESCRIPTION

[0063] The present invention will now be further described with reference to the accompanying drawings and examples.

[0064] like Figure 1 As shown, the present invention provides a method for acquiring a sound field of a two-dimensional sound field model based on a finite element method, which is a finite element method for calculating a two-dimensional underwater sound propagation problem. The method comprises:

[0065] According to the ocean environment parameters of the target area, the target area is discretized; the target area is discretized into a finite number of calculation units, each of which is connected by nodes; according to the unit shape and the number of nodes, the basis function of each calculation unit is constructed;

[0066] The ocean environment parameters include: seabed sound velocity, seabed density, seabed topography of the target area, sound source intensity and sound source location. Since the triangular unit is more adaptable to the complex ocean environment boundary, especially the sea surface and the seabed will present a complex boundary situation, the triangular unit is used to divide the physical domain (target area).

[0067] Specifically, a three-node triangle unit is used as the basic calculation unit. According to the ocean environment parameters of the target area, the target area is divided and discretized into a finite number of three-node triangle units. The basis function of each basic calculation unit is N n (x,z),n=1,2,3; where n represents the local number of the node, (x,z) is the established rectangular coordinate system; the z axis passes through the sound source and is vertically downward; the x axis is parallel to the sea surface;

[0068] Assume that the coordinates of the three nodes are represented as A n (x n ,y n ), n = 1, 2, 3, and define the determinant of the Jacobi matrix |J| = (x2-x1)(y3-y1)-(x3-x1)(y2-y1), where x n is the horizontal coordinate of the nth node; y n is the ordinate of the nth node; the expressions of the three basis functions corresponding to the three nodes of the calculation unit are obtained;

[0069]

[0070]

[0071]

[0072] Using the standard Galerkin variational method, the weak form of the Helmholtz equation is derived;

[0073] Specifically, for the infinite line source problem, a rectangular coordinate system (x, z) is established, with the z axis passing through the sound source and pointing vertically downward, and the x axis parallel to the sea surface; at this time, the Helmholtz equation in the inhomogeneous medium is expressed as

[0074]

[0075] Where ρ represents the seawater density and k(x,z) represents the wave number, which can be expressed as Where c(x,z) represents the speed of sound, f is the frequency of the sound source; S ω represents the intensity of the sound source, z s Indicates the depth of the sound source; is the gradient operator; p is the sound pressure; δ is the Dirac function;

[0076] The weak form of the governing equation is derived using the standard Galerkin variational method. The specific process is as follows:

[0077] Multiply both sides of the above Helmholtz equation by the test function q(x,z) and integrate it over the solution domain Ω of the target area to obtain the following expression:

[0078]

[0079] Applying Gauss's theorem to the above expression, we obtain the following expression, which is the weak form of the Helmholtz equation (the governing equation to be solved):

[0080]

[0081] Among them, Γ represents the boundary, Ω represents the solution domain of the target area;

[0082] The Lagrange interpolation function (basis function) corresponding to the three-node triangle unit is used to discretize the sound pressure and test function of the target area to obtain the discrete sound pressure p(x, z) and the discrete test function q(x, z);

[0083]

[0084] q(x,z)=N m (x,z),m=1,2,3

[0085] The discretized Helmholtz equation without considering boundary conditions is derived

[0086]

[0087] Among them, p n is the sound pressure value of the nth node to be determined; Yes n The coefficient is used as the element k of the unit stiffness matrix mn ; n and m are local numbers of nodes, n = 1, 2, 3, m = 1, 2, 3; N n and N m are the basis functions corresponding to nodes n and m respectively; and N n and N m The gradient of

[0088] For the unit external load vector, we substitute the discretized test function q(x,z) into the right side of the weak form of the Helmholtz equation to obtain Take it as the mth element of the element external load vector.

[0089] Establish the unit stiffness matrix of each calculation unit, and integrate the unit stiffness matrices of all calculation units into the overall stiffness matrix according to the one-to-one correspondence between the local numbers of the nodes and the overall numbers;

[0090] Specifically, the unit stiffness matrix of each calculation unit is established, and all calculation units are traversed. According to the corresponding relationship between the local number of each node in the calculation unit and the overall number of each element in the overall stiffness matrix, each element of the unit stiffness matrix is ​​added to each corresponding element of the overall stiffness matrix, and then the elements of the unit stiffness matrix are integrated into the corresponding positions in the overall stiffness matrix.

[0091] For example, let us assume that the element stiffness matrix is ​​k, which for a three-node triangular element consists of nine elements: k 11 ,k 12 ,k 13 ,k 21 ,k 22 ,k 23 ,k 31 ,k 32 ,k 33 If the local numbers 1, 2, and 3 of the nodes in the unit correspond to the global numbers 23, 34, and 35 respectively, then the nine elements of the unit stiffness matrix are added to the nine corresponding elements of the global stiffness matrix. Let the global stiffness matrix be K, then these nine elements are K 23,23 ,K 23,34 ,K 23,35 ,K 34,23 ,K 34,34 ,K 34,35 ,K 35,23 ,K 35,34 ,K 35,35 .

[0092] Establish the unit external load vector of each calculation unit, and integrate the unit external load vectors of all calculation units into the overall external load vector according to the one-to-one correspondence between each local number of the node and each overall number;

[0093] Specifically, the method of integrating the unit external load vector into the overall external load vector is the same as the stiffness matrix.

[0094] Establish the unit external load vector of each calculation unit, traverse all calculation units, and add each element of the unit stiffness external load vector to each corresponding element of the overall external load vector according to the correspondence between the local number of each node in the calculation unit and the overall number of each element in the overall external load vector, and then integrate the elements of the unit stiffness external load vector into the corresponding positions in the overall external load vector.

[0095] For example, we assume that the unit external load vector is f. For a three-node triangular unit, the unit external load vector f contains three elements: f1, f2, and f3. If the local numbers 1, 2, and 3 of the unit nodes correspond to the global numbers 23, 34, and 35, respectively, then the three elements of the unit external load vector are added to the three corresponding elements of the global external load vector. If the global external load vector is F, then the three elements are F 23 ,F 34 ,F 35 .

[0096] The influence of boundary conditions is brought into the above-mentioned overall stiffness matrix and overall external load vector. Based on the derived weak form of the Helmholtz equation, the overall stiffness matrix, the sound pressure vector and the overall external load vector are composed into a set of linear algebraic equations as a two-dimensional sound field model (two-dimensional underwater sound propagation model). The two-dimensional sound field model is solved to obtain the sound field, that is, the sound pressure at each node; and then the sound pressure value at each point in the space is obtained through the selected basis function and interpolation method.

[0097] Specifically, the influence of boundary conditions includes: first-type boundary conditions and second-type boundary conditions; among which, the first-type boundary conditions are Dirichlet boundary conditions;

[0098] For the first type of boundary conditions, the sound pressure is equal to a constant, and the "diagonal elements set to one" method is used;

[0099] For the second type of boundary conditions, the common boundary condition is that the normal derivative of the sound pressure is zero. In this case, the first term on the left side of the Helmholtz equation is Zero, no special processing is required.

[0100] For an infinite area, the sound field satisfies the radiation condition, i.e., the extinction condition: when r→∞, the sound pressure p→0;

[0101] At this time, it is necessary to truncate the infinite area and use artificial boundary conditions to simulate the emission conditions. We choose the perfect matching layer method to simulate the radiation conditions; the perfect matching layer is an absorption layer that is considered to be set to absorb the outgoing waves at various angles;

[0102] Use a perfectly matched layer to simulate an infinite area, perform coordinate transformation, map the real coordinate axis to the complex coordinate axis, make the reflection coefficient at the artificial boundary zero, and absorb the outgoing waves at all angles in the perfectly matched layer;

[0103] In the perfectly matched layer, the real coordinates are transformed into the complex domain. For example, for a medium with a finite depth, we set the horizontal distance of the physical domain to be calculated as R, and set a thickness of δ on the right side of the physical domain. xThe sound velocity and density in the perfectly matched layer are the same as those of the adjacent left medium;

[0104]

[0105] in, is the complex coordinate after coordinate transformation; i is the imaginary unit; ξ is the integral variable; σ x (ξ) is an artificially defined attenuation function, Among them, σ xm is a constant that determines the maximum attenuation of the perfectly matched layer;

[0106] For a perfectly matched layer, p n The coefficient is denoted as k mn , k mn are the elements of the element stiffness matrix in the perfectly matched layers:

[0107]

[0108] Using the elements of the unit stiffness matrix for the perfectly matched layer obtained above, an overall stiffness matrix is ​​established; then a linear algebraic equation group is formed by combining the overall external load vector and the sound pressure vector, which is used as a two-dimensional sound field model, and based on the two-dimensional sound field model, the two-dimensional sound field model is solved to obtain the sound field, that is, the sound pressure of each node;

[0109] When we get the overall stiffness matrix and the overall external load vector, we can solve the unknown node acoustic pressure. It can be written as the following equation

[0110] Kp=F

[0111] Among them, K represents the overall stiffness matrix, F represents the overall external load vector, and p represents the acoustic field, that is, the acoustic pressure vector of each node. For the above-mentioned finite element method, the overall stiffness matrix K is a symmetric sparse matrix. Therefore, we use the commonly used linear algebra software package LAPACK to store the matrix K and solve the linear algebraic equations, improve the memory usage, and increase the calculation speed of the linear equations.

[0112] According to the selected basis function and interpolation method, the sound field at any point in the space is calculated to obtain the sound pressure value of each point in the space. Then, the sound pressure value of each point in the space can be obtained through the basis function and interpolation. According to the analytical solution of the standard problem or the calculation results of the existing standard model, the error analysis is carried out. Specifically, the sound field at any point in the space can be expressed as,

[0113] After the above operations and calculations, we can get the sound pressure p at each node n , the sound pressure at any point in space can be obtained by using the sound pressure at each node and the basis function We can get the sound pressure at any point we need, perform necessary processing on them, and get the required data.

[0114] For standard problems with analytical solutions, we can compare the calculation results of the finite element model with the analytical solution and calculate the error; for standard problems without analytical solutions, we can use the commonly used sound field calculation model to compare the calculation results of the finite element model and calculate the error to verify the correctness of the finite element model. The error is the difference between the result of the finite element model and the analytical solution or reference solution.

[0115] like Figure 2a and 2b As shown in Figure 2, it represents the propagation loss of the double-layer medium calculated using different methods when the sound source frequency is 100Hz in a waveguide with limited depth; Figure 2a represents the standard solution calculated by the wavenumber integration method; Figure 2b represents the finite element solution calculated using the finite element model; Figure 2a and 2b As shown, for a waveguide with limited depth, the results calculated by the semi-analytical wavenumber integration method are in very good agreement with those calculated by the finite element model, and the accuracy of the finite element model is very high.

[0116] like Figure 3 As shown in the figure, the propagation loss at a receiving depth of 60m when the sound source frequency is 100Hz, where the dotted line represents the finite element solution using a perfectly matched layer to process the exit boundary, and the dotted line represents the standard solution calculated by the wave number integration method. It can be seen that the sound fields calculated by the two methods are very consistent, with a maximum error of less than 0.1dB, indicating that the accuracy of the finite element model is very high.

[0117] like Figures 4a-4d As shown, the sound field propagation loss in the penetrable seafloor slope waveguide calculated using different methods when the sound source frequency is 30Hz and 50Hz, the sound source depth is 20m, and the horizontal distance is 50m;

[0118] in, Figure 4a It represents the propagation loss calculated by the simple normal wave model DGMCM when the sound source frequency is 30Hz;

[0119] Figure 4b It represents the propagation loss calculated by the finite element model when the sound source frequency is 30Hz;

[0120] Figure 4c It represents the propagation loss calculated by the simple normal wave model DGMCM when the sound source frequency is 50Hz;

[0121] Figure 4d It represents the propagation loss calculated by the finite element model when the sound source frequency is 50Hz.

[0122] These four figures show that at different sound source frequencies, the results calculated by the finite element model and the reference solution calculated by the simple normal wave model DGMCM are highly consistent and accurate. At the same time, it shows that the finite element model is not only suitable for calculating sound fields at different frequencies, but also for calculating sound fields in ocean environments with horizontal changes.

[0123] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit the present invention. Although the present invention is described in detail with reference to the embodiments, it should be understood by those skilled in the art that any modification or equivalent replacement of the technical solutions of the present invention does not depart from the spirit and scope of the technical solutions of the present invention and should be included in the scope of the claims of the present invention.

Claims

1. A method for acquiring a sound field of a two-dimensional sound field model based on a finite element method, the method comprising: According to the ocean environment parameters of the target area, the target area is discretized; The target area is discretized into a finite number of computational units, each of which is connected by nodes; a basis function of each computational unit is constructed according to the unit shape and the number of nodes; Using the standard Galerkin variational method, the weak form of the Helmholtz equation is derived; Establish the unit stiffness matrix of each calculation unit, and integrate the unit stiffness matrices of all calculation units into the overall stiffness matrix according to the one-to-one correspondence between the local numbers of the nodes and the overall numbers; Establish the unit external load vector of each calculation unit, and integrate the unit external load vectors of all calculation units into the overall external load vector according to the one-to-one correspondence between each local number of the node and each overall number; The overall stiffness matrix, the sound pressure vector and the overall external load vector are combined into a linear algebraic equation system as a two-dimensional sound field model, and the sound field, that is, the sound pressure at each node, is obtained by solving the two-dimensional sound field model; Then, the sound pressure value at each point in the space is obtained through the selected basis function and interpolation method.

2. The sound field acquisition method of the two-dimensional sound field model based on the finite element method according to claim 1, characterized in that: The ocean environment parameters include: seabed sound velocity, seabed density, seabed topography of the target area, sound source intensity and sound source location.

3. The sound field acquisition method of the two-dimensional sound field model based on the finite element method according to claim 1, characterized in that: According to the obtained ocean environment parameters of the target area, the target area is discretized into a finite number of calculation units, each of which is connected by nodes; according to the unit shape and the number of nodes, a basis function of each calculation unit is constructed; The specific process includes: Using three-node triangle units as calculation units, the target area is divided according to the ocean environment parameters of the target area, and the target area is discretized into a finite number of three-node triangle units. The basis function of each calculation unit is N n (x,z),n=1,2,3; where n represents the local number of the node, (x,z) is the established rectangular coordinate system; the z axis passes through the sound source and is vertically downward; the x axis is parallel to the sea surface; Assume that the three nodes are represented by A n (x n ,y n ), n = 1, 2, 3, and define the determinant of the Jacobi matrix |J| = (x2-x1)(y3-y1)-(x3-x1)(y2-y1), where x n is the horizontal coordinate of the nth node; y n is the ordinate of the nth node; the expressions of the three basis functions corresponding to the three nodes of the calculation unit are obtained:

4. The sound field acquisition method of the two-dimensional sound field model based on the finite element method according to claim 1, characterized in that: The weak form of the Helmholtz equation is derived using the standard Galerkin variational method; the specific process includes: For the infinite line source problem, a rectangular coordinate system (x, z) is established, where the z axis passes through the sound source and is vertically downward, and the x axis is parallel to the sea surface; at this time, the Helmholtz equation in the inhomogeneous medium is expressed as Where ρ represents the seawater density and k(x,z) represents the wave number, which can be expressed as Where c(x,z) represents the speed of sound, f is the frequency of the sound source; S ω represents the intensity of the sound source, z s Indicates the depth of the sound source; is the gradient operator; p is the sound pressure; δ is the Dirac function; The weak form of the governing equation is derived using the standard Galerkin variational method. The specific process is as follows: Multiply both sides of the above Helmholtz equation by the test function q(x,z) and integrate it over the solution domain Ω of the target area to obtain the following expression: Applying Gauss's theorem to the above expression yields the weak form of the Helmholtz equation: Among them, Γ represents the boundary, Ω represents the solution domain of the target area; The Lagrange interpolation function corresponding to the three-node triangle unit is used to discretize the sound pressure and test function of the target area to obtain the discrete sound pressure p(x,z) and the discrete test function q(x,z); q(x,z)=N m (x,z),m=1,2,3 The discretized Helmholtz equation without considering boundary conditions is derived: Among them, p n is the sound pressure value of the nth node to be determined; Yes n The coefficient is used as the element k of the unit stiffness matrix mn ; n and m are local numbers of nodes, n = 1, 2, 3, m = 1, 2, 3; N n and N m are the basis functions corresponding to nodes n and m respectively; and N n and N m The gradient of For the external load vector of the unit, the discretized test function q(x,z) is substituted into the right side of the weak form of the Helmholtz equation to obtain Take it as the mth element of the element external load vector.

5. The sound field acquisition method of the two-dimensional sound field model based on the finite element method according to claim 1, characterized in that: The unit stiffness matrix of each calculation unit is established, and the unit stiffness matrices of all calculation units are integrated into the overall stiffness matrix according to the one-to-one correspondence between the local numbers of the nodes and the overall numbers; the specific process includes: Establish the unit stiffness matrix of each calculation unit, traverse all calculation units, and add each element of the unit stiffness matrix to each corresponding element of the overall stiffness matrix according to the correspondence between the local number of each node in the calculation unit and the overall number of each element in the overall stiffness matrix, and then integrate the elements of the unit stiffness matrix into the corresponding positions in the overall stiffness matrix.

6. The sound field acquisition method of the two-dimensional sound field model based on the finite element method according to claim 1, characterized in that: The overall stiffness matrix, the sound pressure vector and the overall external load vector are combined into a linear algebraic equation group as a two-dimensional sound field model, and the two-dimensional sound field model is solved to obtain the sound field, that is, the sound pressure of each node; and then the sound pressure value of each point in the space is obtained by selecting a basis function and an interpolation method; The specific process includes: For an infinite area, the sound field satisfies the radiation condition: when r→∞, the sound pressure p→0; At this time, it is necessary to truncate the infinite area, use artificial boundary conditions to simulate the emission conditions, and select the perfect matching layer method to simulate the radiation conditions; the perfect matching layer is an absorption layer that is considered to be set to absorb the outgoing waves at various angles; Use a perfectly matched layer to simulate an infinite region, perform coordinate transformation, and map the real coordinate axis to the complex coordinate axis. The specific process is as follows: For a medium with limited depth, assume that the horizontal distance of the physical domain to be calculated is R, and set a thickness of δ on the right side of the physical domain. x The sound velocity and density in the perfectly matched layer are the same as those of the adjacent left medium; in, is the complex coordinate after coordinate transformation; i is the imaginary unit; ξ is the integral variable; σ x (ξ) is an artificially defined attenuation function, Among them, σ xm is a constant that determines the maximum attenuation of the perfectly matched layer; For a perfectly matched layer, p n The coefficient is denoted as k mn , k mn are the elements of the element stiffness matrix in the perfectly matched layers: Using the elements of the unit stiffness matrix for the perfectly matched layer obtained above, an overall stiffness matrix is ​​established; then a linear algebraic equation group is formed by combining the overall external load vector and the sound pressure vector as a two-dimensional sound field model, and the two-dimensional sound field model is solved to obtain the sound field; Kp=F Among them, K represents the overall stiffness matrix, F represents the overall external load vector, and p represents the acoustic field; Then, according to the obtained sound field, the sound pressure value of each point in the space is obtained by selecting the basis function and interpolation method.

Images