Direct problem for the hyperbolic first-order system of acoustic waves propagation Текст научной статьи по специальности «Математика»

DIRECT PROBLEM FOR THE HYPERBOLIC FIRST-ORDER SYSTEM OF ACOUSTIC WAVES PROPAGATION

S. I. Kabanikhin1'2'3, D. V. Klyuchinskiy1'2, I. M. Kulikov1, N.S. Novikov1'2, M.A. Shishlenin1'2'3

1 Institute of Computational Mathematics and Mathematical Geophysics SB RAS, 630090, Novosibirsk, Russia 2 Novosibirsk State Univerity, 630090, Novosibirsk, Russia 3 Sobolev Institute of Mathematics, 630090, Novosibirsk, Russia,

UDC 517.968

DOI: 10.24411/9999-016A-2019-10030

We investigate the mathematical modeling of the 2D acoustic waves propagation in homogenious and heterogenious areas. The hyperbolic first order system of partial differential equations is considered and solved by the Godunov method of the first order of approximation and MUSCL-Hancock method of the second order of approximation. Keywords: acoustics, tomography, first-order hyperbolic system, direct problem, inverse problem, Godunov method, MUSCL-Hancock method

Introduction

In this paper we consider the direct problem for the hyperbolic system of the two dimensional acoustic wave propagation. This problem is related to ultrasound tomography [4,5], aimed at the detection of inclusions in the soft human tissue, which is connected with the development of methods and instruments for early breast cancer detection. An important issue is the construction of effective algorithms for solving the inverse problem, and, therefore, direct problem due to their tightly connected efficiency.

We use the hyperbolic first-order system to describe the propagation of the acoustic waves, written in the form of conservation laws of continuum mechanics. This allows to propose more realistic model from the physical point of view. We consider two numerical algorithm for solving direct problems: the classical Godunov scheme of the first-order approximation [6,7,10] and Monotone Upstream-centred Scheme for Conservation Laws of the second-order approximation investigated by Leer [8] and modified by Hancock [9]. MUSCL-Hancock scheme is based on Godunov approach and the main difference from last one is that it uses piece-wise linear approximation of values inside grid cells instead of piece-wise constant approximation by Godunov method. The numerical methods, based on the Godunov approach, allows to find the balance between mathematical modeling of physical process and the effective numerical realization.

The work has been supported by RSCF under grant 19-11-00154 "Creation of conceptionally new mathematical models of acoustic tomography in medicine. Numerical methods' high-performance computing and software systems".

ISBN 978-5-901548-42-4

1 Direct problem

1.1 Main equations

The propagation of the acoustic waves in the 2D medium is described by the following system of integral conservation laws

pudxdy + pdydt = 0, (1)

pvdxdy + pdxdt = 0, (2)

/pdxdi + p^ + ,dxd,) = 0. (3)

Here u = u(x,y, t),v = v(x,y, t) are components of velocity vector, with respect to x and y respectively, p = p(x, y, t) is the acoustic pressure. The parameters of the system describes the properties of the medium: p(x, y) denotes the density of the medium, and c(x, y) is the speed of the wave propagation.

Solution u(x,y, t), v(x,y, t) and p(x,y, t) of integral equations (1)-(3), in general, can be both smooth and discontinuous types. In our problem we will consider only smooth acoustic waves.

1.2 Model of acoustic wave propagation from the source

Let us consider the following direct problem for acoustic wave propagation from the source

^ + 1 TT =0, (x,y) € o,0 <t <T (4)

at p ox

dv 1 dp . . „

^ + -piy =0, (x,^ € n, (5)

d ( du dv \

+ pc2{ Vx + Yy ) = v)m> {x> y) e ^ (6)

Pl(x,y)edQ =0. (7)

u,v,plt=0 =0 (8)

Here n = (x, y) € [0, L] x [0, L]. Function 9q describes the location of the source, which has the form of Ricker wavelet: ^

1 - 2Ut- I e"""1^"^ (9)

For solving the direct problem (4)—(8), in the next section we consider numerical method, proposed by S.K. Godunov in [7]. The main element of this method is the solving of Riemann problem, which consists of an initial value problem composed of a conservation equation together with piece-wise constant data having a single discontinuity. This problem can be solve explicitly or by using some approximation [10].

1.3 Godunov finite-difference scheme

Let us introduce the mesh in R2, provided by the lines x = x« and y = yj, and denote by hx, hy parameters of the mesh. We will also denote by pair of indices (i + 1/2, j + 1/2) the grid cell Ii+1/2,j+1/2, corresponded to points (xj, yj), (xj, yj+i), (xi+i, yj) ,(xi+\, yj+i). We suppose, that values of functions u(x, y, t),v(x,y, t),p(x,y, t), which describes the state of the medium, are constant in each cell, and denoted this values as ui+1/2^+1/2, vi+1/2^+1/2, pi+1/2,j+1/2. Now let us consider the ratio (1) for the cell Ii+1/2,j+1/2:

iii \ ito +T iVi+1 p[ul+1/2,3+1/2 -ui+1 /2d+! ¡2) hxhy + I I [p (xi+i ,y, t) - p (xi ,y, t)]dydt = 0, (10)

The upper index indicates the values of u, v,p for the moment tQ + t (t is the time step), and the lower index corresponds to the values of functions at the moment t q. We approximate the functions p(xi,y, t), p(xi+\,y, t),

corresponded to values at the boundary of the cell, by some constant values, which we denote as Pij+i/2 and Pi+i,j+i/2 respectively. Thus, we can obtain:

p (vi+i/2,j+i/2 - Ui+i/2j+i/2) hxhy + Thy (Pi+ij+i/2 - Pi,j+i/2) = 0. (11) We use the conservation laws (2), (3) to obtain two more ratios:

p («i+i/2'j+i/2 - vi+i/w/2) hxhy + rhx (Pi+i/2,j+i - P+i/u) = 0, (12) (pl+l/2,3+l/2 - Pi+i/2,j+i/2) hxhy +tpc2hy (Ui+i,j+i/2 - Ui,j+i/2) +

+Tp(?hx (Vi+i/2,j+i - Vi+i/2,j) = T0ni+1/2J+1/2Ik. (13)

where Ik - approximation of the source on time step to. We use the solution of the problem of "decay of discontinuity" to obtain values U, P at the boundary x = Xi

ui-i/2,j+i/2 + ui+i/2,j+i/2 1 Pi+i/2,j+i/2 - Pi-i/2,j + i/2 , ,

U = u^+i/2 =-2--^-2--(14)

p_n _ Pi-i/2,j+i/2 + Pi+i/2,j+i/2 ui+i/2,j+i/2 - ui-i/2,j+i/2 ( .

P = Pi,3 + i/2 = -2--pC-2-

The same ratios for V and P can be obtained at the boundary y = yj

V = V = Vi+i/2,j-i/2 + vi+i/2,j + i/2 1 Pi+i/2,j + i/2 - Pi+i/2,j-i/2 (16)

i+i/2'j 2 pc 2 ()

ft+i/2,j-i/2 + ft+i/2,j + i/2 ^+i/2,j+i/2 - Vj+i/2,j-i/2 , ,

P = Pi+i/2,j =-2--pc-2--()

Therefore, according to (11), (12), (13), we can write the final ratios for computing the state of the medium for "upper" layer = o + :

u*+i/2,3 + i/2 = u. + i/2,. + i/2 - -h- (Pi+ij + i/2 - Pi,o + i/2) , (18)

^+i/2,^ + i/2 = V+i/2,3 + i/2 - ph~y (Pi+i/U+i - Pi+i/2,j) , (19)

r+i/2,3 + i/2 = pr+i/2iJ + i/2 - -pc2 (U+^ + i/2 - Uitj + i/2)

lx

T „2

- TiJc2 (Vi+i/2,3 + i - Vi+i/2,j) - ^ni+i/2,i + i/2 ^ - (20)

hy

As was mentioned previously, we solve the problem in the domain Q = (x, y) G [0, L] x [0, L]. We would like to have somehow modeled non-reflecting conditions near the boundaries x = 0, y = 0, x = L,y = 0. To model this situation we do not set exact boundary conditions at "big" values P0j+i/2,U0^+i/2, PN,j+i/2,U^,j+i/2, Pi+i/2,0, Ui+i/2,0, Pi+i/2,N, Ui+i/2,N. Instead of this we just rewrite "small" values in the outside cells of the area in the next way:

ui/2,j+i/2 = u3/2,j+i/2, uN-i/2,j + i/2 = uN-3/2,j+i/2 vi/2,j+i/2 = V3/2,j+i/2, VN-i/2,j+i/2 = VN-3/2,j + i/2 Pi/2,j+i/2 = P3/2,j+i/2, PN-i/2,j + i/2 = PN-3/2,j+i/2

for boundaries x = 0 and x = L, and the equivalent formulas for boundaries y = 0 and y = L.

These conditions guaranties us an attenuation of waves near the domain boundaries and physically means, that waves just leave the considered area.

Here values U, V, P are defined by (14)—(17). The stability of the scheme can be achieved by using the CFL condition [7,10]:

T i^hx + hy ) < 1

The constant C0 here is the upper-bound estimate for the speed of wave propagation in the computational domain.

1.4 MUSCL-Hancock finite-difference scheme

In this section for solving (4)—(8) we use the MUSCL-Hancock finite-difference scheme of the second order. This scheme has a classical "predictor-corrector" type. The idea of the method is to use piece-wise linear values of variables in cells instead piece-wise constant values as in classical Godunov method.

For simplicity, let us represent our system in conservative form in 1D case without right-hand side

dU dF(U) , ,

m +~ir = °, (21)

Here is U - is a vector of three components U = (u, , ), F(U) - nonlinear function of U. The MUSCL-Hancock approach achieves a second order extension of the Godunov first order upwind method if the intercell flux Fi is computed according to the following steps:

Step 1. Predictor step. Data reconstruction.

In the data reconstruction step, data cell average values U™+1/2 are locally replaced by piece-wise linear functions in each cell Ii+1/2 = [xj,xj+1], namely

(x — xj+1 /2) U+i/2(x) = U+1/2 + ^ a+1/2) At

Ai +1/2 is a suitably chosen slope vector of Ui+1/2(x) in cell Ii+1/2. The extreme points x = 0 and x = Ax, in local coordinates, correspond to the intercell boundaries xi and xi+1 in global coordinates. The values Ui+1/2(x) at the extreme points are

Ui+1/2 = Ui+1/2 - +1/2 Ui+1/2 = Ui+1/2 + 1Ai+1/2

where super indexes L and R denotes left and right extreme point, respectively. These values, also, are usually called boundary extrapolated values. Note, that U and Ai+1/2 are vectors of three components of the Euler equations and thus there are six extrapolated values. Step 2. Corrector step. Evolution.

For each cell Ii+1/2, the boundary extrapolated values UL+1/2 and UR+1/2 are corrected (evolved) by a time ^, according to

uL =U L + 1 A

i+1/2 = U*+1/2 + 2 Ax

TTr - TTR 1 At

i+1/2 = Ui+1/2 + 2Ax

F(UL+1/2) -F(UR+1/2) F(UL+1/2) -F(UR+1/2)

Note that this evolution step is entirely contained in each cell Ii+1/2, as the intercell fluxes are evaluated at the boundary extrapolated values of each cell. At each intercell boundary i there are two fluxes, namely F(URR1/2) and F(UL+1/2), which are in general distinct. This does not really affect the conservative character of the overall method, as this step is only an intermediate step; the intercell flux Fi to be used is yet to be evaluated. Step 3. The Riemann Problem.

To compute the intercell flux Fi one now solves the conventional Riemann problem with data

_R _L

UL = U i-1/2, UR = U i+1/2

to obtain the similarity solution Ui(x/t). The intercell flux Fi is now computed in exactly the same way as in the Godunov first order upwind method from previous section, namely

Fi = F (Ui(0))

Here Uj(0) denotes the value of Ui(x/t) at x/t = 0, i.e. the value of Ui(x/t) along the t-axis. There the ten possible wave patterns in the solution of the Riemann problem that must be taken into account when computing the Godunov flux.

A possible choice for the slope vector Ai+1/2 is

Ai+1/2 = 2(1 + u)Ai + 2(1 - u)Ai+1

where

A,

Tin Tin

Ui+i/2 - Ui-i/2,

A

■i+i

Tin Tin

Ui+3/2 - Ui+i/2

and w G [-1,1].

Also, in the cases of high values of gradient (discontinuous solutions) a technology of slope limiter can be used. There are a lot of different limiter in literature, we present only two of them.

Ai+i/2 = max(0, min(^Aj, Aj+i), min(Ai; fiAi+i)), if Am > 0

Ai+i/2 = min(0, max(^Aj, Aj+i), max(Aj, fiAi+i)), if Aj+i < 0

Here fl =1 corresponds to MINMOD limiter, and fl = 2 to SUPERB EE limiter. In our problem there are no discontinuous solutions, that's why the choose of special limiter has a very little influence on the final solution. Step 4. Conservation laws.

Computed "big" values of flux Fi are used in finite-difference conservation laws to obtain "small" values of vector U = (u, v,p) inside cells at time step t + t0

U = U - T- (Fi+i - Fi) h

2 Numerical results

In this section we present results of solving direct problem on the grid 400 x 400 points. Physical size of area [0 : 0.3] x [0 : 0.3]. CFL number was taken equal to 0.5. We consider heterogenious area, which simulates some object inside a water. This object has two insertions with higher density and velocity than object has. Fig. 1 shows the distribution of density and velocity in this area.

Figure 1: Density and velocity distribution in the area. Fig. 2 shows the modeling results of pressure for different time points.

Conclusion

We have developed a direct solver for the hyperbolic first-order system of acoustics equations via Godunov method and Muscl-Hancock methods (with two types of limiters). Numerical experiments show wave propagation from the source without oscillations and quite similar results for these methods due to smooth type of wave. Each finite-difference scheme has been developed via OpenMP multi-threading technology. The results of speedup is good enough, but its description is beyond the scope of this article.

We have considered the mathematical model of acoustic tomograph in the form of the system of hyperbolic equations of the first order because in this form written conservation laws of continuum mechanics. It allows us to control the preservation of the basic invariants during the solution of direct and inverse problems. This is important for solving unstable problems, as the conservation laws of the main invariants are the only criterion

t = 0.195 t = 0.285

Figure 2: Acoustic wave. Pressure at different time points.

of the well-posedness of the solution. Using the problem formulation in the hyperbolic system we control the numerical solution and we guarantee that numerical solution is close to the physical solution of the process.

The method of the solution of the direct problem for the acoustic tomograph which is based on the the two dimensional hyperbolic system of acoustic equations is developed in [15]. The approach for direct problem solution, proposed in this article, will be used to solve coefficient inverse problem for acoustics wave propagation.

References

[1] Kang H., Tanuma K., Inverse problems for scalar conservation laws, Inverse Problems, 21, P. 1047-1059 (2005)

[2] Berres S, Burger R., Coronel A., Sepulveda M., Numerical identification of parameters for a strongly degenerate convection-diffusion problem modelling centrifugation of flocculated suspensions, Applied Numerical Mathematics, 52, P. 311-337 (2005)

[3] Holden H., Priuli F.S., Risebro N.H., On an inverse problem for scalar conservation laws, Inverse Problems, 30, No. 035015, 35pp (2014)

[4] Burov V.A., Vecherin S.N., Morozov C.M., Rumyantseva O.D., Modeling of the exact solution of the inverse scattering problem by functional methods, Acoustic Journal, 56, 4, p. 516-536 (2010) (in Russian).

[5] Burov V.A., Voloshinov V.B., Dmitriev K.V., Polikarpova N.V., Acoustic waves in metamaterials, crystals, and anomalously refracting structures, Advances in Physical Sciences, 54, 1165-1170 (2011)

[6] Godunov S.K., Differential method for numerical computation of noncontinuous solutions of hydrodynamics equations, Matematicheskiy Sbornik, 47, 3, p. 271-306 (1959) (in Russian)

[7] Godunov S.K., Zabrodin A.V., Ivanov M.Y., Kraikov A.N., Prokopov G.P.: Numerical solution for multidimensional problems of gas mechanics, Moscow, Nauka, 1976.

[8] B. van Leer. MUSCL, A New Approach to Numerical Gas Dynamics. In Computing in Plasma Physics and Astrophysics, Max-Planck-Institute Plasma Physik, Garchung, Germany, April 1976.

[9] B. van Leer. On the Relation Between the Upwind-Differencing Schemes of Godunov, Enguist-Osher and Roe. SIAM J. Sci. Stat. Comput., 5(1):1-20, 1985.

10] Godunov S.K., Kulikov I.M., Computation of Discontinuous Solutions of Fluid Dynamics Equations with Entropy Nondecrease Guarantee, Computational Mathematics and Mathematical Physics 54, 6, p. 1008-1021 (2014)

11] S.I. Kabanikhin, D.B. Nurseitov, M.A. Shishlenin, B.B. Sholpanbaev, Inverse problems for the ground penetrating radar, Journal of Inverse and Ill-Posed Problems, 21, No. 6. P. 885-892 (2013)

12] S.I. Kabanikhin, M.A. Shishlenin, Quasi-solution in inverse coefficient problems, Journal Inverse and Ill-Posed Problems, 16, No. 7. P. 705-713 (2008)

13] S.I. Kabanikhin, M.A. Shishlenin, Numerical Methods for Solving Inverse Hyperbolic Problems, Computational Methods for Applied Inverse Problems (Ed. by Y. Wang, A. Yagola, and C. Yang). Berlin-Boston-Beijing: De Gruyter and Higher Education Press, P. 369-393 (2012)

14] S.I. Kabanikhin, M.A. Shishlenin, Acoustic sounding methods of linearization and conversion of the wave field , Siberian Electronic Mathematical Reports, 7, P. C.199-C.206 (2010) Mathematical Notes. 28 (4).P. 723—727 (1981)

15] Kulikov I.M., Novikov N.S., Shishlenin M.A. Mathematical modelling of ultrasound wave propagation in two dimensional medium: direct and inverse problem. Siberian Electronic Mathematical Reports, 12, Pp. C.219-C.228 (2015)

Kabanikhin Sergey — Dr. Sc., senior fellow of the Institute of Computational Mathematics and Mathematical Geophysics SB RAS;

e-mail: ksi52@mail.ru;

Klyuchinskiy Dmitriy — student of the Novosibirsk State University;

e-mail: dmitriy_klyuchinskiy@mail.ru; Kulikov Igor — Dr. Sc., senior fellow of the Institute of Computational Mathematics and Mathematical Geophysics SB RAS;

e-mail: kulikov@ssd.sscc.ru; Novikov Nikita — PhD, senior fellow of the Institute of Computational Mathematics and Mathematical Geophysics SB RASy;

e-mail: novikov-1989@yandex.ru; Shishlenin Maxim — Dr. Sc., Head of the laboratory of the Institute of Computational Mathematics and Mathematical Geophysics SB RAS; Novosibirsk State University;

e-mail: mshishlenin@ngs.ru. Received — April 30, 2019