Next Article in Journal
Antioxidant Production in Dunaliella
Previous Article in Journal
Experimental Investigation and Numerical Modeling of Elastic Modulus Variation with Stress during Hydration and Expansion Process of Static Cracking Agent
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Coupled Mode Sound Propagation in Inhomogeneous Stratified Waveguides

1
Acoustic Science and Technology Laboratory, Harbin Engineering University, Harbin 150001, China
2
Key Laboratory of Marine Information Acquisition and Security, Ministry of Industry and Information Technology, Harbin Engineering University, Harbin 150001, China
3
College of Underwater Acoustic Engineering, Harbin Engineering University, Harbin 150001, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2021, 11(9), 3957; https://doi.org/10.3390/app11093957
Submission received: 22 February 2021 / Revised: 24 April 2021 / Accepted: 25 April 2021 / Published: 27 April 2021
(This article belongs to the Section Acoustics and Vibrations)

Abstract

:
An efficient coupled mode method for modeling sound propagation in horizontally stratified inhomogeneous waveguides, in which the seabed is modeled as a (layered) acoustic medium, is presented. The method is based on Fawcett’s coupled mode method and the multimodal admittance method. The acoustic field is expanded onto the unusual local eigenfunctions composed by normal modes in the corresponding one-layer homogeneous waveguides with constant depth equal to the local total depth of the multilayered waveguide. A set of energy-conserving first-order differential equations governing the modal amplitudes of acoustic fields is derived. The admittance method is employed to solve the differential equations in a numerically stable manna. The coupled mode method considers the backscattering effect of inhomogeneities and full coupling between local modes, and offers improvement from the viewpoint of efficiency and computational cost. The acoustic fields predicted by the method agree well with those computed by the commercial finite element software COMSOL Multiphysics. The method can be extended to further establish fast and accurate 3D sound propagation models in complex shallow water environments.

1. Introduction

Sound propagation in shallow water in the presence of inhomogeneities has been studied for several decades, conventionally modeled by the nonseparable Helmholtz equation in horizontally stratified inhomogeneous waveguides [1,2,3,4,5]. The step-wise coupled mode method is commonly used for sound propagation in waveguides in which a step-wise variation of the properties with range occurs [3,4]. However, for the sound propagation problem in waveguides with continuously varying properties, the stepwise approximation of the properties renders the numerical solution procedure cumbersome [6]. The continuous coupled mode method is a classical approach to sound propagation problems in range-dependent waveguides, originally presented by Pierce [1] and Milder [2]. They expanded the sound field in terms of usual local modes with range-dependent modal amplitudes and derived a system of coupled mode equations containing two coupled coefficients. Rutherford and Hawker indicated that the Pierce–Milder coupled mode equations do not conserve the energy among modes because the usual local modes satisfy different boundary conditions than the acoustic pressure [7]. Fawcett derived a coupled system of differential equations containing two coupling matrices and two interface matrices, and treated this coupled mode system by a finite difference scheme [8]. It is noted that, the coupling matrices and interface matrices, which require a knowledge of derivatives of usual local modes, can only be approximately computed [9]. Fawcett’s approach considers the backscattering effect from range-dependent environments and full coupling between modes, and conserves the energy among modes. However, the complexity involved in computing coupling matrices and interface matrices along with the complex numerical scheme for solving differential equations makes this approach impractical, especially in long range and high frequency applications, where the computational requirements increase substantially. In order to reduce computational complexity, several asymptotic coupled mode methods have been developed. The adiabatic coupled mode method neglects mode coupling among modes and provides accurate predictions when the range dependence of a medium is sufficiently weak [5,10]. The one-way approximation coupled mode method is efficient for sound propagation in waveguides where the backscattering effect can be ignored [5,11]. However, asymptotic coupled mode models can provide very poor solutions when the range dependence of a medium is significant. Thus, it is still very important to study an efficient two-way coupled mode method for sound propagation in complex underwater waveguides.
The multimodal admittance method presented by Pagneux [12,13] is a two-way coupled mode method for wave propagation in waveguides with range-dependent cross-sections. This method basically consists in rewriting the Helmholtz equation into a set of first-order evolution equations and then projecting onto the usual local modes to obtain a set of first-order differential equations governing the modal amplitudes of the wave fields. In the case where these differential equations are to be integrated numerically, the problem of numerical divergence owing to the presence of evanescent modes appears. The admittance matrix, namely, the Dirichlet-to-Neumann operator in modal domain, is introduced to avoid numerical problems. Using the admittance matrix, first-order differential equations are reformulated into the Riccati equation governing the admittance matrix and the first-order evolution equation governing the modal amplitudes of acoustic pressure. Both the Riccati equation and the evolution equation can be computed with classical numerical schemes [12,13,14,15]. This method is efficient and numerically stable. Liu and Li [16] derived a coupled system of first-order differential equations for sound propagation in simple one-layer waveguides with range-dependent environments based on Fawcett’s coupled mode method and the multimodal admittance method. We extend the method here to present an efficient coupled mode scheme for sound propagation and wave scattering problems in complex inhomogeneous stratified waveguides.
This paper focuses on the two-way solutions in horizontally stratified waveguides with both bathymetric and volumetric variations. In particular, waveguides with penetrable scatterers are also inspected, modeled by horizontally stratified waveguides with intersecting interfaces. We note that the paper deals with acoustic waveguides in which the seabed is also modeled as a (layered) acoustic medium. In many applications, the modeling of seabed needs to account for the shear rigidity (seabed elasticity); see, e.g., [17] and the references cited therein. In an attempt to reduce the complexity involved in the computation of usual local modes and their derivatives for complex multilayered environments, the local modes for expansion of acoustic pressure in this paper are defined by the transverse eigenfunctions in a one-layer waveguide with constant depth equal to the local total depth of the multilayered waveguide and constant physical parameters. We note that such local modes and their derivatives can be analytically computed. A set of second-order coupled mode equations governing the modal amplitudes is derived by projecting the local mode onto the Helmholtz equation. The energy conservation is guaranteed by introducing interface matrices. Remarkably, we introduce a set of unknown quantities with respect to the modal amplitudes and their horizontal derivatives. Using this unknown vector, the second-order coupled mode equations are reduced into a set of first-order differential equations with respect to the modal amplitudes and the unknown quantities. This first-order coupled system preserves the energy conservation property, but no interface matrix is contained. The required coupling matrices can be computed with numerical integral methods since their integrands are analytical. Moreover, it is straightforward to implement the admittance method for numerically stable solutions. The use of the unknown vector substantially increases the efficiency of the method. Since a more realistic representation of the ocean bottom is a penetrable infinite half space, the application of our coupled mode method to waveguides with a penetrable infinite half space is studied on basis of the elimination of branch cut and branch line integrals [18]. Numerical results for sound propagation in inhomogeneous multilayered waveguides predicted by commercially available finite element software COMSOL multiphysics are also presented to validate the efficiency and accuracy of this method.
The remainder of the paper is organized as follows: Section 2 presents the basic formulation of two-way and one-way asymptotic coupled mode methods for sound propagation in horizontally stratified inhomogeneous waveguides. The sound field expression generated by an arbitrary incident wave in such an environment is provided. By constructing mapping between incident waves and wave fields, optimal incident waves resulting in sound self-focusing at any position [19] in the multilayered waveguide are studied. Section 3 illustrates the sound fields predicted by our coupled mode method and COMSOL for a complex multilayered waveguide, and the sound self-focusing pattern in this numerical example is shown. Section 4 presents the extension of the present coupled mode formalism to waveguides with the effect of internal waves and a penetrable bottom. Section 5 gives the conclusions.

2. Method

2.1. Two-Dimensional Two-Way Coupled Mode Equations

Let us start with sound propagation in the multilayered waveguide sketched in Figure 1 with interfaces h j ( x ) , j = 1 , 2 , , J 1 , separating the layers. The upper water column is confined between a flat pressure-release surface and underlying sediment layers of any shape. The Neumann boundary condition is imposed at the bottom of the lower layer. The whole region is divided into three subregions. Subregion 0 x x R corresponds to the scattering region where both volumetric and bathymetric variations exist. In subregions x < 0 (not shown) and x > x R , mass density and sound speed are assumed to be only depth-dependent, and boundaries and interfaces are horizontal. This assumption allows for consistently formulating the source and radiation conditions by means of normal-mode expansions.
The Helmholtz equation governing the acoustic pressure is
x 1 ρ ( x , z ) p ( x , z ) x + z 1 ρ ( x , z ) p ( x , z ) z + ω 2 ρ ( x , z ) c 2 ( x , z ) p ( x , z ) = 0 ,
supplemented by boundary conditions
p ( x , h 0 ) = 0 , p n x , h J ( x ) = 0 ,
and continuous conditions
p ( x , h j ( x ) ) = p ( x , h j + ( x ) ) , 1 ρ j ( x , h j ( x ) ) p ( x , h j ( x ) ) n = 1 ρ j + 1 ( x , h j + ( x ) ) p ( x , h j + ( x ) ) n ,
where j = 1 , 2 , , J 1 , c ( x , z ) is sound speed, and ρ ( x , z ) is density. c ( x , z ) and ρ ( x , z ) could exhibit sharp discontinuities at the interfaces.
According to the multimodal method, acoustic pressure can be expressed as a sum of transverse eigenfunctions ψ n ( z ; x ) :
p ( x , z ) = n = 0 N 1 P n ( x ) ψ n ( z ; x ) ,
where N is the truncation number, P n ( x ) are modal amplitudes, and ψ n ( z ; x ) are defined to be normal modes of a one-layer homogeneous waveguide with constant depth equal to the local depth of the multilayered waveguide, satisfying the following equations:
2 ψ n ( z ; x ) z 2 + γ n 2 ψ n ( z ; x ) = 0 , ψ n ( h 0 ; x ) = 0 , ψ n ( h J ( x ) ; x ) z = 0 .
Eigenfunctions ψ n ( z ; x ) in the series representation have analytical expressions, i.e.,
ψ n ( z ; x ) = 2 h J ( x ) sin ( n + 0.5 ) π z h J ( x ) , n = 0 , 1 , 2 , , N .
In contrast, the usual local modes that traditional coupled mode method used in the series representation for acoustic pressure, namely, transverse modes at ( x , z ) of the complex multilayered waveguide, can only be numerically computed. The use of eigenfunctions ψ n ( z ; x ) simplifies computations in this sense.
Then, a projection of the Helmholtz equation onto ψ m ( z ; x ) gives
0 h J ( x ) ψ m * x 1 ρ p x + z 1 ρ p z + ω 2 ρ c 2 p d z = 0 .
By substituting Equation (4) into Equation (7) and applying Fawcett’s approach [8], one can obtain second-order coupled-mode equations
n A m n P n + n ( B m n + 2 C m n + D m n ) P n + n ( E m n + F m n + G m n + K m n L m n ) P n = 0 ,
where A m n , K m n , and L m n are correction coefficients resulting from the use of unusual eigenfunctions ψ n ( z ; x ) ; B m n , C m n , F m n , and G m n are mode coupling coefficients; and D m n and E m n are interface coefficients given by, respectively,
A m n = 0 h J ( x ) 1 ρ 1 ψ m * ψ n d z + j = 1 J 1 h j ( x ) h j + 1 ( x ) 1 ρ j + 1 1 ρ 1 ψ m * ψ n d z ,
B m n = 0 h J ( x ) x 1 ρ 1 ψ m * ψ n d z + j = 1 J 1 h j ( x ) h j + 1 ( x ) x 1 ρ j + 1 1 ρ 1 ψ m * ψ n d z ,
C m n = 0 h J ( x ) 1 ρ 1 ψ m * ψ n x d z + j = 1 J 1 h j ( x ) h j + 1 ( x ) 1 ρ j + 1 1 ρ 1 ψ m * ψ n x d z ,
D m n = 1 ρ J ( h J ) h J ( x ) ψ m * ( h J ) ψ n ( h J ) j = 1 J 1 1 ρ j + 1 ( h j ) 1 ρ j ( h j ) h j ψ m * ( h j ) ψ n ( h j ) ,
E m n = 1 ρ J ( h J ) h J ( x ) ψ m * ( h J ) ψ n ( h J ) x j = 1 J 1 1 ρ j + 1 ( h j ) 1 ρ j ( h j ) h j ψ m * ( h j ) ψ n ( h j ) x ,
F m n = 0 h J ( x ) 1 ρ 1 ψ m * 2 ψ n x 2 d z + j = 1 J 1 h j ( x ) h j + 1 ( x ) 1 ρ j + 1 1 ρ 1 ψ m * 2 ψ n x 2 d z ,
G m n = 0 h J ( x ) x 1 ρ ψ m * ψ n x d z + j = 1 J 1 h j ( x ) h j + 1 ( x ) x 1 ρ j + 1 1 ρ 1 ψ m * ψ n x d z ,
K m n = 0 h J ( x ) 1 ρ 1 ω 2 c 1 2 ψ m * ψ n d z + j = 1 J 1 h j ( x ) h j + 1 ( x ) 1 ρ j + 1 ω 2 c j + 1 2 1 ρ 1 ω 2 c 1 2 ψ m * ψ n d z ,
L m n = 0 h J ( x ) 1 ρ 1 ψ m * z ψ n z d z + j = 1 J 1 h j ( x ) h j + 1 ( x ) 1 ρ j + 1 1 ρ 1 ψ m * z ψ n z d z .
Since ψ n ( z ; x ) satisfy different boundary and interface conditions than acoustic pressure when the bottom boundary and interfaces are range dependent, interface coefficients D m n and E m n are introduced in the derivation to guarantee energy conservation among modes by properly applying the boundary and interface conditions. Details of the derivation of D m n and E m n are presented in Appendix A. The interface coefficients can be analytically computed. Other coefficients can be numerically computed by numerical integration method. In this paper, we use the Clenshaw–Curtis method [20] to calculate the correction and coupling coefficients.
Although the fully two-way coupled system Equation (22) conserves the energy among modes, as any classical coupled mode method, the local modal series representation Equation (4) has a slow rate of convergence with the modal amplitudes following the order O ( n 2 ) , where n is the mode number. The poor convergence is attributed to the differences in the interface and boundary conditions satisfied by ψ n ( z ; x ) and acoustic pressure. There are several ways that allow for accelerating the convergence of the local mode series [6,21,22,23,24]. One of them [22,23] uses the idea of boundary modes. It is based on an enriched modal series representation in terms of supplementary boundary modes, which are orthogonal to the local eigenfunctions but satisfy different boundary and interface conditions. Boundary modes have the beneficial effect of acting as evanescent modes and do not change the form of the coupled mode equations. However, the numerical implementation of the boundary modes for multilayered waveguides is complicated and generally requires an amount of computing power. We therefore only consider the nonimproved version of coupled mode equations for sound propagation analysis in multilayered environments.
In order to solve the second-order coupled mode equations, we define a set of unknown quantities S m , m = 0 , 1 , , N 1 , by
S m = n = 0 N 1 A m n P n + n = 0 N 1 C m n P n .
The x-derivative of S m gives
S m = n = 0 N 1 A m n P n + ( A m n + C m n ) P n + C m n P n ,
where
A m n = B m n + C m n + C n m + D m n
C m n = E m n + F m n + G m n + H m n ,
H m n = 0 h J ( x ) 1 ρ 1 ψ m * x ψ n x d z + j = 1 J 1 h j ( x ) h j 1 ( x ) 1 ρ j + 1 1 ρ 1 ψ m * x ψ n x d z .
By substituting Equations (18) and (19) into Equation (8), and rewriting the resulting equation into matrix form, we have
S = C T A 1 S + ( L K ) P ,
where P = ( P n ) , S = ( S m ) , and the elements in matrices A , C , L , and K are A m n , C m n , L m n and K m n , respectively. In Equation (20), we use the identical relation H C T A 1 C 0 . To prove this relation, we expand x ψ n ( z ; x ) by the eigenfunctions ψ n ( z ; x ) as
ψ n ( z ; x ) x = n = 0 N 1 ψ n ( z ; x ) Q n n ( x ) .
Substituting Equation (21) into the expressions of C and H yields C = AQ and H = Q T A Q . Thus, one can obtain H = C T A 1 C .
From Equations (18) and (20), we obtain the first-order coupled mode equations
d d x P S = A 1 C A 1 L K C T A 1 P S .
We note that, to calculate acoustic fields, only four of the nine matrices in Equation (8) remain to be computed. The first-order coupled system Equation (22) preserves the energy conservation property but does not contain any interface matrix. For details of energy conservation of Equation (22), see Appendix B. It is the proper application of S that gives the coupled mode equations a simple form and substantially reduces computing power, especially in long-range application.
The first-order differential equations account for source and radiation conditions. Directly integrating Equation (22) starting from the radiation condition may encounter the problem of numerical divergence owing to the exponential growth of evanescent modes. In order to avoid numerical problems, we introduce admittance matrix Y and propagator matrix M by S ( x ) = Y ( x ) P ( x ) and P ( x ) = M ( x ) P ( 0 ) , respectively. Y satisfies a Riccati equation and M a first-order differential equation coupled to Y:
Y = Y A 1 Y + Y A 1 C + C T A 1 Y + L K ,
M = A 1 C M + A 1 Y M .
Y can be computed by integrating from right to left with the initial value Y ( x R ) = A ( x R ) L ( x R ) K ( x R ) . Eigenvalues of matrix A 1 ( x R ) Y ( x R ) are equal to i k x , where k x are horizontal wavenumbers in region x > x R (for more details, see Appendix C). M ( x ) can be solved from left to right with the initial value M ( 0 ) = I , where I is the identity matrix. In this paper, numerical integration of Equations (23) and (24) is performed using the Magnus method [25,26]. Further, the reflection and transmission matrices that characterize the scattering properties of the scattering region can be computed. Reflection matrix R is defined by P r ( 0 ) = R P i ( 0 ) , and transmission matrix T is defined by P t ( x R ) = T P i ( 0 ) , where P i , P r , and P t denote the modal amplitudes of the incident, reflected and transmitted waves, respectively. We obtain
R = Y 0 + Y ( 0 ) 1 Y 0 Y ( 0 ) ,
and
T = M ( x R ) ( I + R ) ,
where Y ( 0 ) is the admittance matrix at x = 0 , and Y 0 = A ( 0 ) L ( 0 ) K ( 0 ) is the local admittance matrix at x = 0 .
Eventually, one can write the acoustic pressure in the multilayered waveguide as
p ( x , z ) = Ψ T ( z ; x ) M ( x ) ( I + R ) P i ( 0 ) , 0 x x R , Ψ T ( z ; x R ) e A 1 ( x R ) Y ( x R ) ( x x R ) T P i ( 0 ) , x > x R ,
where Ψ is a vector with elements being ψ n . To summarize, once the geometry and medium parameters of the waveguide are identified, admittance matrix Y and propagator matrix M are available for all x. Then, reflection matrix R and transmission matrix T can be deduced. Finally, the sound field generated by any incident wave can be obtained. Although the derivation of the coupled mode system in the paper is on basis of multilayered waveguides in which the properties are continuously varying with range in each layer, the method can be naturally extended to solve sound propagation problems for waveguides in which sharp discontinuities of the properties with range exist; see, e.g., [27,28,29].

2.2. Two-Dimensional One-Way Coupled Mode Equations

For waveguides where the range dependence of medium is weak, one-way approximation coupled-mode formulation provides a speed-up way to solve acoustic pressure by neglecting backscattering energy [30]. This means that elements in reflection matrix R are relatively small in this case. From Equation (25), it follows directly that the admittance matrix at any position approximates the local admittance matrix, and the iterative computation of admittance matrix is no longer necessary. Thus, the admittance matrix takes the following form:
Y ˜ ( x ) = A ( x ) ( L ( x ) K ( x ) ) .
The one-way approximation coupled-mode equation reads
P = ( A 1 C + A 1 Y ˜ ) P ,
with initial condition P i ( 0 ) = P + ( 0 ) and P + ( 0 ) representing the modal amplitudes of forward propagating wave.
Similarly, we introduce propagator matrix M ˜ ( x ) by P ( x ) = M ˜ ( x ) P ( 0 ) . M ˜ ( x ) satisfies the equation M ˜ = ( A 1 C + A 1 Y ˜ ) M ˜ with the initial value M ˜ ( 0 ) = I . Then, acoustic pressure under one-way approximation can be written as
p ( x , z ) = Ψ T ( z ; x ) M ˜ ( x ) P + ( 0 ) .

2.3. Optimal Incident Waves for Sound Focusing

As a byproduct of our coupled mode algorithms, the optimal incident wave that leads to sound focusing at a certain position in a multilayered waveguide is investigated.
From Equation (27), the acoustic pressure at any position ( x 0 , z 0 ) can be rewritten as
p ( x 0 , z 0 ) = Ψ T ( x 0 , z 0 ) M ( x 0 ) ( I + R ) , P i ( 0 ) , 0 x 0 x R , Ψ T ( x 0 , z 0 ) exp A 1 ( x R ) Y ( x R ) ( x 0 x R ) T , P i ( 0 ) , x 0 > x R ,
where < a 1 , a 2 > = a 1 a 2 denotes the scalar product of vectors in Hilbert space, and ‘ ’ denotes the conjugate transpose. The energy flux of the incident wave is assumed to be unity. Then, acoustic pressure modulus at ( x 0 , z 0 ) can reach its maximum value, max p i | p ( x 0 , z 0 ) | , when the modal amplitudes of the incident wave take the following form:
P i ( 0 ; x 0 , z 0 ) = 1 Λ Ψ ( x 0 , z 0 ) T M ( x 0 ) ( I + R ) , 0 x 0 x R , 1 Λ Ψ ( x 0 , z 0 ) T exp A 1 ( x R ) Y ( x R ) ( x 0 x R ) T , x 0 > x R ,
where Λ is fixed by the normalization of the sound energy flux of incident wave E p i = Im ( 0 h 2 ( 0 ) ρ 1 p i * x p i d y ) / 2 ω .

3. Numerical Results

In this section, we validate the present coupled mode method by comparing the results for a numerical example with those obtained using the acoustics module of the commercially available finite element software COMSOL Multiphysics and show the sound focusing pattern by injecting optimal incident waves.
Let us consider sound propagation in a two-layer waveguide with a penetrable scatterer being located in the upper water column. We model such an inhomogeneous waveguide as a virtual four-layer waveguide with intersecting interfaces. In scattering region ( 0 < x < 400 m), sound speed and mass density are c 1 = ( 1500 0.1 z + 0.1 x ) m/s, ρ 1 = ( 1000 + 0.1 z + 0.1 x ) kg/m 3 , c 2 = ( 1700 0.2 z + 0.1 x ) m/s, and ρ 2 = ( 1500 + 0.2 z + 0.1 x ) kg/m 3 . In the transmission region ( x 400 m), parameters are c 1 = ( 1540 0.1 z ) m/s, ρ 1 = ( 1040 + 0.1 z ) kg/m 3 , c 2 = ( 1740 0.2 z ) m/s, and ρ 2 = ( 1540 + 0.2 z ) kg/m 3 . The scatterer has 20m radius, ( x s = 350 , z s = 25 ) m center, and is fluid-filled with c s = 1800 m/s and ρ s = 1600 kg/m 3 . The shapes of the interface and bottom are formulated by Equations (33) and (34), respectively,
h 1 ( x ) = 40 + 10 1 cos ( π x / 400 ) , 0 < x < 400 , 60 . x 400 .
h 2 ( x ) = 100 10 1 cos ( π x / 400 ) , 0 < x < 400 , 80 , x 400 .
The waveguide is excited by a distributed source p i ( z ) = sin ( 1.5 π z / 100 ) from the left at frequency f = 100 Hz . There are 10 propagating modes in the left lead and 12 propagating modes in the right lead. From radiation condition Y ( x = 400 ) , Y , M , R , and T are available. Then, acoustic pressure can be computed using Equation (27) with the source condition P i = ( 0 , 1 / 0.02 , 0 , , 0 ) T . Figure 2a,b, respectively, show the acoustic-pressure moduli calculated by the present two-way coupled-mode method (CMM) and COMSOL. The results were computed on a quad-core 3.6GHz processor running the 64-bit Windows operating system. The truncation number used for the local series representation is N = 20 . The average computation time per output grid [31] is 1.4 × 10 4 s. Clearly, the sound field wildly oscillates in the scattering region, which is expected for the case where the backscattering energy is large and the mode coupling effect is significant. Figure 2c illustrates the acoustic pressure distributions in the horizontal direction at a depth z = 20 m predicted by the two-way CMM, one-way CMM and COMSOL. The spatial discretization and truncation numbers used in the two- and one-way CMM are identical. We see that the calculations by the two-way CMM and COMSOL are in excellent agreement. The slight differences are attributed to the truncation of the infinite modal series Equation (4), and the interpolation difference between COMSOL and our code. It can be minified, of course, by increasing truncation number N and improving mesh quality. However, the results predicted by the one-way CMM are remarkably different, meaning that backscattering energy cannot be neglected in this environment. Since the modal method only necessitates the trace of the the modal amplitudes of the first N local basis functions, the present coupled mode algorithm is much more efficient than COMSOL is, which is a direct Helmholtz solver, for sound propagation over long-range distances.
Next, the acoustic wave self-focusing pattern in a multilayered waveguide is presented. The geometry and medium parameters of the waveguide are kept the same as those in Figure 2; the frequency is f = 100 Hz. Figure 3a shows wave focusing at ( 450 , 20 ) m in the transmission region by injecting the incident wave of Figure 3b, which is computed from Equation (32). Figure 3c plots the corresponding distribution of acoustic pressure along the z direction at x = 450 m. We see that energy in the vicinity of the focus occupies the vast majority of energy along the z direction. This method allows one to enhance the energy density of any position in a multilayered waveguide.

4. Application to Waveguides with Penetrable Infinite Half Space and Internal Wave

In this section, we apply the coupled mode formalism to sound propagation problems in a multilayered waveguide, considering the effect of a penetrable infinite half space and internal waves, which is a common scenario in real ocean environments.
In previous sections, we assumed that the bottom is acoustically rigid. However, a more realistic representation of the ocean bottom is an infinite half space that allows for energy penetrating into the half space. The rigid bottom assumption neglects one or more components of sound fields depending on the choice of branch cut, for instance, leaky modes plus an integral along the Pekeris cut for Pekeris branch cut, or an integral of the continuous spectrum along the EJP cut for EJP branch cut. Difficulties in the application of bounded coupled mode models to a Pekeris waveguide consist of the discrete approximation of continuous spectrum. An available algorithm eliminates the branch cut and branch line integrals by introducing a small complex sound speed gradient in the half space [18]. The branch cut contribution is represented by discrete modes. This method has the advantage that leaky modes do not increase in amplitude with increasing depth. Since branch line modes make a minor contribution to the sound field in the water column, it is possible to ignore the branch line modes in numerical implementation [32]. Therefore, it provides a straightforward way to incorporate with the couple mode algorithms. We extend this approach here to coupled mode sound propagation in multilayered waveguides.
We consider the waveguide shown in Figure 4.
A two-layer system of sound speed in the water column is used. The upper mixed layer has sound speed c m l and water density ρ 1 . The lower water layer has sound speed c 1 and the same density ρ 1 . The shape of the internal solitary wave is formulated by
h i n ( x ) = h + Δ h sech 2 x x w a v e L ,
where Δ h is the amplitude of the internal wave, x w a v e denotes the wave center, and L is the width parameter [33]. The ocean bottom is characterized by a penetrable infinite half space with sound speed c 2 and density ρ 2 . To apply the coupled mode formalism, the penetrable infinite half space is approximately reproduced by a bounded layer. The corresponding sound speed c ˜ 2 has a small gradient,
c ˜ 2 = c 2 1 + i α / α c ,
where the attenuation α is in Gaussian form, α = α 0 e ( z h 2 ( x ) ) 2 / ( 2 λ 2 2 ) , α 0 is in the unit of dB / λ 2 , λ 2 = c 2 / f , and α c = 40 π log 10 e . The thickness of the lower bounded layer, h 2 h 1 ( x ) , is meant to be larger than 3 λ 2 , such that attenuation increases gently in the beginning, then grows rapidly with increasing depth, and finally smoothly reaches its maximum at z = h 2 . We use such a form of α in an attempt to reduce the spurious reflection that can return to the water column and provide accurate solutions in the water column using an acceptable truncation depth.
Figure 5 shows an example of computation of the sound field excited by a point source in the waveguide shown in Figure 4. The physical parameters are ρ 1 = 1000 kg/m 3 , c m l = 1530 m/s, c 1 = 1490 0.1 ( z h i n ( x ) ) m/s, c 2 = 1700 m/s, and ρ 2 = 1500 kg/m 3 . The internal wave shape is
h i n ( x ) = 50 + 50 sech 2 x 0.5 x R 1200 ,
where x R = 16 km. The water–sediment interface is sloping,
h 1 ( x ) = 200 0.005 x , 0 < x < x R , 120 , x x R .
In numerical implementation, a false rigid boundary is placed at h 2 = 200 + 5 λ 2 beneath the water column, and the sound speed in the bottom layer is replaced by c ˜ 2 (Equation (36)) with α 0 = 10 dB / λ 2 . The point source is located in ( x s , z s ) = ( 0 , 30 ) m at a frequency of 50 Hz. The source condition, which is represented by the modal amplitudes of the Green’s function in the waveguide of Figure 4, is [20]
P i = Ψ ( z s ; x s ) 2 ( A 1 ( x s ) Y x s ) 1 exp ( A 1 ( x s ) Y x s | x x s | ) ,
where Ψ is a vector with elements being functions ψ n , and Y x s = A ( x s ) L ( x s ) K ( x s ) . The acoustic pressure can be computed using Equation (27) with the source condition and radiation conditions. In addition, the eigenvalues of matrix i A 1 ( x ) L ( x ) K ( x ) , which fall in the interval ω / c 2 , max ( ω / c 1 ( z ) ) , correspond to the horizontal wavenumbers of the trapped modes. At the chosen frequency, there are six trapped modes at x = 0 and four trapped modes in the transmission region. Figure 5a presents the computed sound field using the present two-way coupled mode method. The truncation number is N = 30 , which has been proved to be enough for numerical convergence. Figure 5b compares the transmission loss predicted by the two-way CMM, one-way CMM, and COMSOL. The numerical model used for COMSOL is a waveguide with a PML beneath the water column. The receiver depth is z r = 60 m. The transmission loss is defined by
TL ( x , z r ) = 10 log 10 | p ( x , z r ) | 2 | p ( x = 1 , z s ) | 2 .
We see that the agreement between the two-way CMM and COMSOL is excellent. Differences are attributed to the deviations in spatial discretization. The results of the one-way CMM are accurate for the range internal from 0 m to about 6 km, implying that the backscattered energy is small in this region. As sound propagates through the internal wave, waves interact with the sloping interface and internal waves, and differences in the results between the one- and two-way CMM increase gradually. Since the computation model used in COMSOL mimics the real waveguide with a penetrable infinite half space, one can conclude that acoustic pressure in the water column of a waveguide with a penetrable infinite half space and internal waves predicted by our two-way coupled mode method is accurate.

5. Conclusions

An efficient coupled mode method is presented for sound propagation through multilayered inhomogeneous waveguides based on Fawcett’s coupled mode method and the multimodal admittance method. The inhomogeneous environments include range-dependent physical parameters and range-dependent interfaces (boundaries). The method uses the unusual local transverse eigenfuctions in the series representation of the acoustic pressure, which are defined to be the transverse modes in a one-layer homogeneous waveguide with constant depth equal to the local total depth of the multilayered waveguide and constant physical parameters. The unusual local modes have analytical solutions. Projecting the Helmholtz equation onto the unusual local modes yields a set of second-order differential equations with respect to the modal amplitudes of acoustic pressure. Energy conservation among modes is guaranteed by the introduction of two interface matrices in the case of a lossless acoustic waveguide. By properly applying a set of unknown quantities S in terms of the modal amplitudes of acoustic pressure and their horizontal derivatives, the second-order coupled mode equations are reduced into a set of first-order differential equations. This first-order coupled system preserves the energy conservation property but requires no knowledge of interface matrices. The admittance method was employed to solve the first-order differential equations for avoiding numerical problems. The reflection and transmission matrices that characterize the scattering properties are also given. Our coupled mode method considers the full coupling between modes and the backscattering effect and substantially reduces computing power compared with fully numerical methods, especially in long-range and high-frequency applications. Numerical results show that our two-way coupled mode method provides us with accurate solutions for multilayered inhomogeneous waveguides with acoustically rigid bottom or a penetrable infinite half space. Our coupled system is appropriate for two-way solutions in long-range/high-frequency propagation and scattering problems and can be naturally extended to treat sound propagation through real ocean environments.

Author Contributions

Conceptualization, J.L. and Q.L.; methodology, J.L.; software, J.L.; validation, J.L.; formal analysis, J.L.; investigation, J.L.; resources, Q.L.; data curation, J.L.; writing—original draft preparation, J.L.; writing—review and editing, Q.L.; visualization, J.L. and Q.L.; supervision, Q.L.; project administration, Q.L.; funding acquisition, Q.L. Both authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Interface Matrices D and E

Here, we present the derivation of D m n and E m n . Integrating by part for the projection of the term / z ( ρ 1 p / z ) in Equation (1), and using the following boundary and interface conditions,
p ( h J ) y = h J p ( h J ) x , 1 ρ j p ( h j ) y h j p ( h j ) x = 1 ρ j + 1 p ( h j + ) y h j p ( h j + ) x ,
where j = 1 , 2 , , J 1 , one can obtain
0 h J ψ m * z 1 ρ p z d z = j = 0 J 1 ψ m * 1 ρ j + 1 p z h j + h j + 1 j = 0 J 1 h j + h j + 1 ψ m * z 1 ρ j + 1 p z d z = 1 ρ J ( h J ) h J ( x ) ψ m * ( h J ) p ( h J ) x j = 1 J 1 1 ρ j + 1 ( h j ) 1 ρ j ( h j ) h j ( x ) ψ m * ( h j ) p ( h j ) x 0 h J ( x ) 1 ρ 1 ψ m * z p z d z j = 1 J 1 h j ( x ) h j + 1 ( x ) 1 ρ j + 1 1 ρ 1 ψ m * z p z d z .
Substituting the mode expansion Equation (4) into Equation (A2), we have
0 h J ψ m * z 1 ρ p z d z = n D m n P n + n ( E m n L m n ) P n .

Appendix B. Energy Conservation

Let us start with the energy conservation property of the boundary value problem of Helmholtz equation Equation (1). Supposing that p and p * are two solutions, we have
x 1 ρ p x + z 1 ρ p z + ω 2 ρ c p = 0 ,
x 1 ρ p * x + z 1 ρ p * z + ω 2 ρ c p * = 0 ,
By applying the operator 0 h J ( x ) (Equation (A4) × p * Equation (A5) × p ) d z , one can obtain
0 h J ( x ) x p * 1 ρ p x p 1 ρ p * x + z p * 1 ρ p z p 1 ρ p * z d z = 0 .
Equation (A6) can be reduced to
x 0 h J ( x ) p * 1 ρ p x p 1 ρ p * x d z + j = 0 J 1 p * 1 ρ j + 1 p n p 1 ρ j + 1 p * n h j + ( x ) h j + 1 ( x ) = 0 ,
where the second term is equal to zero accounting for the interface and boundary conditions, and the first term corresponds to the conservation of energy flux E ( x ) ,
x 0 h J ( x ) p * 1 ρ p x p 1 ρ p * x d z E ( x ) = d d x 1 2 ω Im 0 h J ( x ) 1 ρ p * p x d z = 0 .
Then, translating Equation (A8) on our modal amplitudes, we have
x P S S P d d x ( P S ) = 0 .
The modal amplitudes satisfy the coupled mode equations Equation (22) with ( A 1 C ) T + C T A 1 = 0, A 1 = ( A 1 ) T , and L K = ( L K ) T , meaning that our numerical scheme preserves the energy conservation property. It follows directly that without the introduction of interface matrices D and E , the coupled mode equations will not conserve the energy among modes.

Appendix C. Radiation Condition

The initial value of the admittance matrix (Equation (23)) corresponds to the radiation condition, which is represented by the admittance matrix in the range-independent region x > x R . The transverse normal modes ϕ i ( z ) in a range-independent multilayered waveguide are solutions of the eigenproblem
d d z 1 ρ ( z ) d ϕ i ( z ) d z + ω 2 ρ ( z ) c 2 ( z ) ϕ i ( z ) = 1 ρ ( z ) k x i 2 ϕ i ( z ) ,
with boundary conditions ϕ i ( 0 ) = 0 , ϕ i ( h J ) = 0 , and continuous conditions ϕ i ( h j ) = ϕ i ( h j + ) and ϕ i ( h j ) / ρ j ( h j ) = ϕ i ( h j + ) / ρ j + 1 ( h j + ) , j = 1 , 2 , , J 1 , where k x i are the axial wavenumbers. To solve ϕ i , we first employ a mode expansion,
ϕ i ( z ) = n = 0 N 1 w i , n ψ n ( z ) ,
where w i , n are the modal amplitudes, and ψ n ( z ) are transverse normal modes in a homogeneous waveguide with a constant depth equal to h J , i.e., ψ n ( z ) = 2 / h J sin ( n + 1 / 2 ) π z / h J . Applying the operator
0 h J ( x ) ψ m * ( · ) d z ,
one can obtain
A 1 ( K L ) W = W K x 2 ,
where the elements in W are w i , n , and K x is a diagonal matrix with diagonal elements being k x i . Therefore, the eigenvalue problem governed by Equation (A10) is transformed into a matrix eigenvalue problem of Equation (A12). The axial wavenumbers k x i and the modal amplitudes w i , n can be obtained by eigen-decomposing the matrix A 1 ( K L ) .
Alternatively, for environments of which the properties are invariant in the horizontal, the coupling among modes is lost. The second-order coupled mode equation is then reduced to
P = A 1 ( L K ) P .
Since S = Y P and P = A 1 S , the admittance matrix in the transmission region takes the following form
Y = A ( L K ) .

References

  1. Pierce, A.D. Extension of the Method of Normal Modes to Sound Propagation in an Almost-Stratified Medium. J. Acoust. Soc. Am. 1965, 37, 19–27. [Google Scholar] [CrossRef]
  2. Milder, D.M. Ray and Wave Invariants for SOFAR Channel Propagation. J. Acoust. Soc. Am. 1969, 46, 1259–1263. [Google Scholar] [CrossRef]
  3. Evans, R.B. A coupled mode solution for acoustic propagation in a waveguide with stepwise depth variations of a penetrable bottom. J. Acoust. Soc. Am. 1983, 74, 188–195. [Google Scholar] [CrossRef]
  4. Luo, W.; Schmidt, H. Three-dimensional propagation and scattering around a conical seamount. J. Acoust. Soc. Am. 2009, 125, 52–65. [Google Scholar] [CrossRef] [Green Version]
  5. Jensen, F.B.; Kuperman, W.A.; Porter, M.B.; Schmidt, H. Computational Ocean Acoustics, 2nd ed.; Springer: New York, NY, USA, 2011. [Google Scholar]
  6. Belibassakis, K.A.; Athanassoulis, G.A.; Papathanasiou, T.K.; Filopoulos, S.P.; Markolefas, S. Acoustic wave propagation in inhomogeneous, layered waveguides based on modal expansions and hp-FEM. Wave Motion 2014, 51, 1021–1043. [Google Scholar] [CrossRef]
  7. Rutherford, S.R.; Hawker, K.E. Consistent coupled mode theory of sound propagation for a class of nonseparable problems. J. Acoust. Soc. Am. 1981, 70, 554–564. [Google Scholar] [CrossRef]
  8. Fawcett, J.A. A derivation of the differential equations of coupled-mode propagation. J. Acoust. Soc. Am. 1992, 92, 290–295. [Google Scholar] [CrossRef]
  9. Abawi, A.T. An energy-conserving one-way coupled mode propagation model. J. Acoust. Soc. Am. 2002, 111, 160–167. [Google Scholar] [CrossRef] [PubMed]
  10. Pierce, A.D. Augmented adiabatic mode theory for upslope propagation from a point source in variable-depth shallow water overlying a fluid bottom. J. Acoust. Soc. Am. 1983, 74, 1837–1847. [Google Scholar] [CrossRef]
  11. Chiu, C.S.; Miller, J.H.; Lynch, J.F. Forward coupled-mode propagation modeling for coastal acoustic tomography. J. Acoust. Soc. Am. 1996, 99, 793–802. [Google Scholar] [CrossRef] [Green Version]
  12. Pagneux, V. Multimodal admittance method in waveguides and singularity behavior at high frequencies. J. Comput. Appl. Math. 2010, 234, 1834–1841. [Google Scholar] [CrossRef] [Green Version]
  13. Pagneux, V.; Amir, N.; Kergomard, J. A study of wave propagation in varying cross-section waveguides by modal decomposition. Part I. Theory and validation. J. Acoust. Soc. Am. 1996, 100, 2034–2048. [Google Scholar] [CrossRef]
  14. Félix, S.; Pagneux, V. Sound propagation in rigid bends: A multimodal approach. J. Acoust. Soc. Am. 2001, 110, 1329–1337. [Google Scholar] [CrossRef]
  15. Félix, S.; Pagneux, V. Multimodal analysis of acoustic propagation in three-dimensional bends. Wave Motion 2002, 36, 157–168. [Google Scholar] [CrossRef]
  16. Liu, J.; Li, Q. A coupled mode method for sound propagation in range-dependent waveguides. Acta Phys. Sin. 2021, in press. [Google Scholar] [CrossRef]
  17. Tsouvalas, A. Underwater Noise Emission Due to Offshore Pile Installation: A Review. Energies 2020, 13, 3037. [Google Scholar] [CrossRef]
  18. Westwood, E.K.; Koch, R.A. Elimination of branch cuts from the normal-mode solution using gradient half spaces. J. Acoust. Soc. Am. 1999, 106, 2513–2523. [Google Scholar] [CrossRef]
  19. Guo, W.; Yang, D.S. Sound focusing in inhomogeneous waveguides. Acta Phys. Sin. 2020, 69, 074301. [Google Scholar] [CrossRef]
  20. Li, Q.; Liu, J.; Guo, W. Sound propagation in inhomogeneous waveguides with sound-speed profiles using the multimodal admittance method. Chin. Phys. B 2020, 29, 14303. [Google Scholar] [CrossRef]
  21. Athanassoulis, G.A.; Belibassakis, K.A. Rapidly-Convergent Local-Mode Representations for Wave Propagation and Scattering in Curved-Boundary Waveguides; Springer: Berlin/Heidelberg, Germany, 2003; pp. 451–456. [Google Scholar]
  22. Maurel, A.; Mercier, J.F.; Pagneux, V. Improved multimodal admittance method in varying cross section waveguides. Proc. R. Soc. A 2014, 470, 20130448. [Google Scholar] [CrossRef]
  23. Mercier, J.F.; Maurel, A. Improved multimodal method for the acoustic propagation in waveguides with a wall impedance and a uniform flow. Proc. R. Soc. A 2016, 472, 20160094. [Google Scholar] [CrossRef]
  24. Bi, W.; Pagneux, V.; Lafarge, D.; Aurégan, Y. An improved multimodal method for sound propagation in nonuniform lined ducts. J. Acoust. Soc. Am. 2007, 122, 280–290. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Lu, Y.Y.; McLaughlin, J. The Riccati method for the Helmholtz equation. J. Acoust. Soc. Am. 1996, 100, 1432–1446. [Google Scholar] [CrossRef]
  26. Iserles, A.; Marthinsen, A.; Nørsett, S.P. On the Implementation of the Method of Magnus Series for Linear Differential Equations. BIT Numer. Math. 1999, 39, 281–304. [Google Scholar] [CrossRef]
  27. Maurel, A.; Mercier, J.F.; Félix, S. Wave propagation through penetrable scatterers in a waveguide and through a penetrable grating. J. Acoust. Soc. Am. 2014, 135, 165–174. [Google Scholar] [CrossRef]
  28. Pagneux, V.; Maurel, A. Lamb wave propagation in inhomogeneous elastic waveguides. Proc. R. Soc. A 2002, 458, 1913–1930. [Google Scholar] [CrossRef]
  29. Tsouvalas, A.; Metrikine, A.V. Noise reduction by the application of an air-bubble curtain in offshore pile driving. J. Sound Vib. 2016, 371, 150–170. [Google Scholar] [CrossRef]
  30. Doc, J.B.; Félix, S.; Lihoreau, B. Coarse-grid computation of the one-way propagation of coupled modes in a varying cross-section waveguide. J. Acoust. Soc. Am. 2013, 133, 2528–2532. [Google Scholar] [CrossRef] [PubMed]
  31. Peng, Y.; Tsouvalas, A.; Stampoultzoglou, T.; Metrikine, A. A fast computational model for near- and far-field noise prediction due to offshore pile driving. J. Acoust. Soc. Am. 2021, 149, 1772–1790. [Google Scholar] [CrossRef]
  32. Knobles, D.P.; Stotts, S.A.; Koch, R.A. Low frequency coupled mode sound propagation over a continental shelf. J. Acoust. Soc. Am. 2003, 113, 781–787. [Google Scholar] [CrossRef] [PubMed]
  33. Shmelev, A.A.; Lynch, J.F.; Lin, Y.T.; Schmidt, H. Three-dimensional coupled mode analysis of internal-wave acoustic ducts. J. Acoust. Soc. Am. 2014, 135, 2497–2512. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. (Color online) Geometry of 2D multilayered waveguide.
Figure 1. (Color online) Geometry of 2D multilayered waveguide.
Applsci 11 03957 g001
Figure 2. (Color online) Sound field (acoustic pressure modulus) in a two-layer waveguide with a penetrable scatterer. (a) Sound field computed by the present two-way coupled mode method with N = 20 . (b) Sound field computed by COMSOL. (c) Acoustic pressure modulus distribution along x direction at z = 20 . The waveguide is forced by a distributed source at 100 Hz.
Figure 2. (Color online) Sound field (acoustic pressure modulus) in a two-layer waveguide with a penetrable scatterer. (a) Sound field computed by the present two-way coupled mode method with N = 20 . (b) Sound field computed by COMSOL. (c) Acoustic pressure modulus distribution along x direction at z = 20 . The waveguide is forced by a distributed source at 100 Hz.
Applsci 11 03957 g002
Figure 3. (Color online) Acoustic waves self-focusing in a two-layer range-dependent waveguide with a penetrable scatterer. The focus is at ( 450 , 20 ) m. (a) Sound field. (b) Acoustic pressure modulus of the optimal incident wave. (c) Acoustic pressure modulus along z-direction at x = 450 m.
Figure 3. (Color online) Acoustic waves self-focusing in a two-layer range-dependent waveguide with a penetrable scatterer. The focus is at ( 450 , 20 ) m. (a) Sound field. (b) Acoustic pressure modulus of the optimal incident wave. (c) Acoustic pressure modulus along z-direction at x = 450 m.
Applsci 11 03957 g003
Figure 4. Approximation model for a waveguide with an internal solitary wave and a penetrable infinite homogeneous half space.
Figure 4. Approximation model for a waveguide with an internal solitary wave and a penetrable infinite homogeneous half space.
Applsci 11 03957 g004
Figure 5. (Color online) Sound field (acoustic pressure modulus) generated by a point source in a waveguide with a penetrable infinite half space and an internal wave. (a) Sound field computed by the present two-way coupled mode method with N = 30 . (b) Transmission loss versus range predicted using the present two-way coupled mode method, one-way coupled mode method, and COMSOL. The point source is located at ( 0 , 30 ) m. The frequency is 50 Hz.
Figure 5. (Color online) Sound field (acoustic pressure modulus) generated by a point source in a waveguide with a penetrable infinite half space and an internal wave. (a) Sound field computed by the present two-way coupled mode method with N = 30 . (b) Transmission loss versus range predicted using the present two-way coupled mode method, one-way coupled mode method, and COMSOL. The point source is located at ( 0 , 30 ) m. The frequency is 50 Hz.
Applsci 11 03957 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, J.; Li, Q. Coupled Mode Sound Propagation in Inhomogeneous Stratified Waveguides. Appl. Sci. 2021, 11, 3957. https://doi.org/10.3390/app11093957

AMA Style

Liu J, Li Q. Coupled Mode Sound Propagation in Inhomogeneous Stratified Waveguides. Applied Sciences. 2021; 11(9):3957. https://doi.org/10.3390/app11093957

Chicago/Turabian Style

Liu, Juan, and Qi Li. 2021. "Coupled Mode Sound Propagation in Inhomogeneous Stratified Waveguides" Applied Sciences 11, no. 9: 3957. https://doi.org/10.3390/app11093957

APA Style

Liu, J., & Li, Q. (2021). Coupled Mode Sound Propagation in Inhomogeneous Stratified Waveguides. Applied Sciences, 11(9), 3957. https://doi.org/10.3390/app11093957

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop