Nonlinear Medical Ultrasound Tomography: 3D Modeling of Sound Wave Propagation in Human Tissues

: The article aimed to show the fundamental possibility of constructing a computational digital twin of the acoustic tomograph within the framework of a unified physics–mathematical model based on the Navier–Stokes equations. The authors suggested that the size of the modeling area is quite small, sound waves are waves of “small” disturbance, and given that a person consists of more than 60% water, human organs can be modeled using a liquid model, taking into account their density. During numerical experiments, we obtained the pressure registered in the receivers that are located on the side walls of the tomograph. The differences in pressure values are shown depending on the configuration of inclusions in the mannequin imitating internal organs. The results show that the developed technology can be used to probe the human body in medical acoustic tomographs and determine the acoustic parameters of the human body to detect neoplasms.


Introducion
Currently, medical tomography is actively used to diagnose a pathologically altered area of the human organ at an early stage of the development of tumors and diseases.In order to increase its efficiency, new mathematical models, numerical solution methods and computer technologies in ultrasound tomography are being developed, which are based on the propagation of acoustic waves in a (liquid) medium.Therefore, the study of the generation, propagation and registration of acoustic waves in a liquid medium [1,2] is an important problem.
At the moment, models based on systems of equations directly related to conservation laws, as well as models based on a second-order hyperbolic equation, are usually used to simulate ultrasound tomography processes.Models based on the second-order wave equation are usually easier to study, and therefore there are more ways to solve the forward problem efficiently, but there is no guarantee that the numerical solution is close to the physical one.The problem formulation based on conservation laws, which we consider in this article, requires more computational resources.The advantage is the close connection with the physics of wave propagation since the equations can be derived directly from conservation laws.Inverse problems for hyperbolic system were studied in [31].
Recently, there have been a lot of works on nonlinear acoustic tomography [32][33][34][35][36][37][38][39].The development of an adequate mathematical model of medical acoustic tomography describing the propagation of waves in the human body, taking into account most physical effects, is an urgent problem.
In this paper, we use three-dimensional numerical modeling based on the solution of the Navier-Stokes equations to simulate the problem of acoustic wave propagation in human soft tissues [40].The proposed approach allows to describe the process of forming, propagation and registration of the acoustic waves, taking into account non-stationary processes of turbulent mixing [41].
We were inspired by the work in which several aspects of ultrasound in medicine, like acoustic nonlinearity, finite amplitude distortion, shock wave formation, harmonic components, and non-linearly induced absorption were considered [42], as well as the saturation and the effect of these effects on ultrasound.
Let us consider some recent results, related to nonlinear acoustic tomography.The method of measuring the acoustic properties of biological tissue elements at the microscopic level was developed on the basis of acoustic microscopy methods with mechanical scanning in the frequency range from 100 to 200 MHz in [43].
In [44], the authors investigated nonlinear effects, such as progressive distortion of the waveform, the generation of frequency harmonics and acoustic shocks, excessive energy release and acoustic saturation, during the propagation of ultrasound in liquids with relatively low acoustic attenuation.Similar effects are observed in soft tissues, although they are limited by absorption and scattering.
Sufficient conditions were given to recover both the sound speed of the medium being probed and the source [11] in thermoacoustic tomography and photoacoustic tomography, two coupled physics imaging modalitiesthat attempt to combine the high resolution of ultrasound and the high contrast capabilities of electromagnetic waves.
B. Kaltenbacher [45] investigated the issues of well-posedness and the effect of attenuation both for classical models of nonlinear acoustics and for some higher-order equations.
In [46], the authors proposed nonlinear acoustic tomography to obtain both the temperature and velocity fields from the time-of-flight measurements of sonic rays.
In [47], a complex magnetohydrodynamic model was investigated, which follows from the widely accepted dynamic theory of geomagnetism.
The theoretic basis of nonlinear acoustics was presented in [48].
The problem of the propagation of disturbances in a liquid medium from a harmonic oscillation source was considered [49].Estimates of the required spatial and temporal resolutions were obtained to ensure acceptable accuracy of the solution.
In [50], the nonlinear acoustics was used to locate and visualize the defects of materials.To test the developed numerical algorithms, it is very important to have phantom materials with acoustic parameters close to the parameters of the human body.In [51], the acoustic properties of seven phantom materials imitating tissues at different concentrations of their compounds and five phantom materials at room temperature were described and compared to determine the most suitable phantom material for contrast-enhanced ultrasound.
In [52], the 3D ultrasound computer tomography was used to increase the resolution in breast cancer imaging.
The inverse problem for a nonlinear wave equation with a damped and nonlinear term in medical acoustic sounding was investigated in [53].
It was established that 3D low-frequency ultrasound tomography is a safe, lowcost, high-resolution, whole-body medical imaging modality and gives high-resolution quantitatively accurate relevant images, which are potentially clinically useful [54,55].
The article discusses one of the most common numerical methods for modeling hydrodynamics: a finite-volume discretization method together with an iterative algorithm SIMPLE [40,56].The final velocity of acoustic wave propagation is taken into account by using the equation of state for a liquid with a given compressibility coefficient.The numerical method is implemented on the LOGOS software ver.5.3.22.837[56].The results of the numerical simulation of the propagation of sound waves in a digital acoustic tomograph, in the center of which is a digital mannequin of a person with internal inclusions (liver, adipose tissue), are presented in order to identify liver pathologies according to acoustic tomography.The mannequin is submerged in water.It is shown how the speed of sound affects the passage of acoustic waves inside the tomograph.

Governing Equations and the Numerical Method
Unsteady three-dimensional viscous fluid flows are described by a system of Navier-Stokes equations, which contain the continuity equation, the momentum conservation equation, and the energy conservation equation [40,57].In a conservative form, in Cartesian coordinates, the system has the form: Here, t > 0 denotes the time, is the velocity, ρ is the density, and τ ij is the viscous stress tensor.
The final velocity of the acoustic wave propagation is taken into account by using the equation of state for a liquid with a given compressibility coefficient.The compression and expansion of the gas passes through an adiabatic process, where k is the polytrope index and is defined as: The system (1) and ( 2) is valid for both laminar and turbulent flows [40].Their direct use for modeling flows is reduced to calculating the viscous stress tensor.
In addition, for the numerical solution, the system (1) and ( 2) must be supplemented by boundary conditions that depend, generally speaking, on the considered problem.In the case of incompressible flows, if the boundary of the calculation area is formed by solid walls, it is necessary to put all velocity components at the boundary equal to the corresponding components of the velocity of the solid surface, i.e., it is impossible for the liquid to slip along the "liquid-solid wall" boundary, or to move along the normal to it.On solid walls, the pressure gradient is zero: The speed value is zero: At the "liquid-liquid" boundary, the velocity and shear stresses must be continuous.
In the presence of an air boundary, zero static pressure is fixed on it, and the gradients of velocity and volume fractions are zero: The non-reflective boundary conditions necessary for the free departure of waves from the computational domain are described in detail in [58].Discretization of the system of Equations ( 1) and ( 2) is carried out by the finite volume method on an arbitrary unstructured grid, and for its numerical solution, a completely implicit method [59,60] based on the SIMPLE algorithm [40].
The modeling of flows with a free surface implies certain modifications of the SIMPLE algorithm.The description of the basic formulas of the modified SIMPLE algorithm, boundary conditions and implementation in the LOGOS software is described in detail [59,61].
We briefly present the discretization of the equations in the volume that is needed for further presentation.Let us consider the finite-volume scheme of the discretization of the equations used.The basic equation for solving the system (1) and ( 2) is the scalar quantity transfer equation: The first term in ( 3) is a nonstationary term, the second is a convective term, and then a diffusion term.The equation may also contain sources and sinks represented by the last member of Q.The tensor τ ij contains spatial derivatives of the desired value φ.We assume for simplicity that τ j = µ ∂φ ∂x j .This significantly reduces calculations, but in no way reduces the generality of the methods described below.
We consider an arbitrary unstructured grid, the view of which is shown in Figure 1: , where i = 0, 1, 2 is the number of the vector component.The vector from the center of cell "P" to the center of cell "M" along the k int face is denoted as d i,k m = r 1,M − r i,P , and the vector from the center of "P" to the center of the face is d k S = r k S − r P , where r i is the radius vector.
Time discretization of Equation ( 3) according to the Adams-Beschfort scheme (the second-order scheme) is: To discretize Equation (3) in space, we integrate it in terms of the volume of the cell "P" and proceed to area integration for the convective and diffusion terms: The discrete analog of the diffusion term is written in the following form: Here, n k is the normal of the face k.In the right part, under the sign of the sum in the product, there is a derivative in the direction ∂φ ∂n k , which, in the case of an orthogonal calculation, can be defined as follows: For approximation on a finite-volume grid, the convective term is written as: Here, F k is the volumetric flow through the face k.The value on the face of φ k is determined by the applied sampling scheme of the convective term.There are a large number of sampling schemes that are applicable on arbitrary unstructured grids [62][63][64][65].Among them, several schemes can be distinguished that have the "highest rating" of applicability in solving applied problems: counterflow scheme (Upwind Differences, UDs), counterflow scheme with linear interpolation (Linear Upwind Differences, LUDs), QUICK scheme, central difference scheme (Central Differences, CDs), schemes of the NVD family (Normalized Variable Diagram) and hybrid circuits (the above circuits mixed with a counterflow circuit to increase monotony).
Using this discretization, Equation ( 3) is replaced by a system of linear algebraic equations written for each calculation cell: The resulting system of linear algebraic equations is solved using the AMG multigrid solver [56].
In an acoustic tomograph for the emission of acoustic waves, there are sources and a receiver that record the pressure.Sources and receivers are cylinders made of a material with acoustic parameters equal to water.The receiver records the entire pressure field inside the cylinder, and then the data of the inverse problem are averaged by volume.
To simulate an acoustic tomograph, the LOGOS software implements a technique that allows you to set the sources of acoustic waves, as well as process data in the receiver.The user can specify an arbitrary number of sources/receivers, specifying their location, size, the normal along which they are directed, and their status (that is, whether this "sensor" is a source or receiver).The source of acoustic waves is given as follows (Riker pulse): where ν 0 is the frequency.The acoustic wave source I(t) is included in the source S i , which is added to the momentum Equation ( 2), and looks like this: Here, ρ is the cell density, V is the cell volume, A is the amplitude, and − → n is the normal in the direction of the signal.The pressure in each receiver is calculated by determining the cells that fall under the specified dimensions and location, and calculating the average pressure in the area of each receiver (volume integration).

Numerical Experiments
In this section, we consider a numerical simulation of the propagation of acoustic waves in the human body in order to study the possibility of determining liver pathologies with acoustic tomography.
Numerical calculations were performed for two cases: in the first case, two cubic inclusions are contained inside the mannequin symmetrically relative to the center-on the right, the inclusion parameters correspond to liver tissue, and on the left, adipose tissue; in the second case, two cubic inclusions are also located inside the mannequin, but the liver is surrounded by adipose tissue 0.5 cm thick.An acoustic tomograph is a cylinder filled with water.The parameters of the mannequin's body correspond to the muscle tissue (see Table 1).Liver pathologies include many diseases, including neoplasms.We considered fatty liver disease, in which there is an excessive accumulation of fat in the liver.To simplify the situation, we assumed that the liver is covered with an additional layer of adipose tissue.In Figure 2, we provide the scheme of the numerical experiment-the mannequin is located inside the tomograph (the cylinder): We use a block grid with a cell size of 0.004 m, consisting of 15 million cells, during the experiment.Figure 3 shows the resulting grid (in vertical and horizontal sections).At the initial moment of time, there are no disturbances, and the liquid in the tomograph is at rest.The digital mannequin parameters are a height of 1.95 m, shoulder width of 0.56 m, and thickness of 0.34 m.Parameters of the tomograph (cylinder) are a height of 2 m and diameter of 0.7 m.Table 1 shows the acoustic parameters of the human body.
The mean parameters of the human body, presented in Table 1, were taken in accordance with [43,51,66,67].
At the first stage, the calculation was carried out with two cubic inclusions, 0.05 × 0.05 × 0.05 m, located symmetrically relative to the center inside the mannequinliver (right) and adipose tissue (left).Eight transducers (sources and receivers) are located on the side of the tomograph at the level of the internal organs (Figure 4).The location of the source of the sounding wave in the calculations is marked with the number 3; all the others are receivers and record the pressure as a function of time.The source (receiver) is a cylinder with a diameter of 0.01 m and a height of 0.05 m made of a material with acoustic parameters equal to those of water.In the receiver, the pressure function captures the entire pressure field inside the cylinder and then averages by volume.During the computations, we used the Ricker wavelet (4) with the frequency ν 0 = 20,000 Hz.In Figure 5, we demonstrate the pressure field at t = 0.0004 s in a section located at a height in the center of the internal organs.At the second stage, a calculation was carried out in which a cube with liver parameters is surrounded by a boundary layer with a thickness of 0.5 cm.Thus, there are two cubic inclusions in the mannequin, located symmetrically relative to the center inside the mannequin: adipose tissue (left) with a size of 0.05 × 0.05 × 0.05 m, and liver (right) with a size of 0.06 × 0.06 × 0.06 m. Figure 6 shows the pressure field at time 0.0004 s in a cross section.In Figure 7, one can see the pressure, registered in the receivers for both experiments.Significant differences in pressure in the two calculations are observed in the 7th receiver.When the wave passes through the liver surrounded by adipose tissue, the pressure amplitude is increased compared to calculation 1, and the passage time of the acoustic waves is decreased due to the fact that the speed of sound in adipose tissue is higher than in water.Minor differences are also observed in receiver 6 (Figure 8).In order to test the efficiency of the technology, we carried out another numerical simulation, similar to calculation 2. The speed of sound in the adipose tissue surrounding the liver we set to 1000 m/s. Figure 10 shows graphs of pressure for three calculations in receivers 6 and 7, where the greatest difference is observed.When a wave passes through the liver surrounded by adipose tissue at the speed of sound of 1000 m/s (experiment 3), the passage time of the acoustic waves to the receivers is increased compared to calculations 1 and 2. The pressure amplitude in receiver 6 is increased compared to experiments 1 and 2, and in receiver 7, it is lower compared to experiment 2.

Discussion
The article presents the results of the numerical simulation of acoustic wave propagation in the human body with nonlinear acoustic tomography based on the Navier-Stokes equations.This mathematical model takes into account most of the nonlinear effects of acoustic wave propagation, including shock waves.The computational technology presented in the article shows the fundamental possibility of using mathematical modeling methods based on the numerical solution of the complete hydrodynamic system of the Navier-Stokes equations to study the properties of acoustic waves during tomography of the human body.The introduction of such methods in the future will make it possible to understand almost all physical aspects of the functioning of such an ultrasound medical tomograph.
The article aimed to show the fundamental possibility of constructing a computational digital twin of the acoustic tomograph within the framework of a unified physicsmathematical model based on the Navier-Stokes equations.The authors proceeded from the fact that the size of the modeling area is quite small, sound waves are waves of "small" disturbance, and given that a human is more than 60% water, human organs can be modeled by a liquid model with their density.This is a kind of approximation in which, as in any model, there are assumptions, but it has the right to life and, in general, will allow us to determine the areas of change in the parameters of the sound wave to identify the required parameters.During numerical experiments, we obtained the pressure registered in the receivers that were located on the side walls of the tomograph.The differences in pressure values are shown depending on the configuration of inclusions in the mannequin imitating internal organs.The results show that the developed technology can be used to probe the human body in medical acoustic tomographs and determine the acoustic parameters of the human body to detect neoplasms.Numerical solution methods are currently being developed and implemented to solve inverse problems [68][69][70][71][72], as well as methods based on deep learning [73][74][75].

Figure 1 .
Figure 1.The calculated grid k is the set of faces of the cell "P", which is composed of a set of inner faces of k int and a set of outer faces of k s .The neighboring cell along the inner face of k int is denoted by M k int .The vector is the area of the face k-S i,k , where i = 0, 1, 2 is the number of the vector component.The vector from the center of cell "P" to the center of cell "M" along the k int face is denoted as d i,k m = r 1,M − r i,P , and the vector from the center of "P" to the center of the face is d k S = r k S − r P , where r i is the radius vector.

Figure 2 .
Figure 2. The digital mannequin of the human body inside the acoustic tomograph.

Figure 3 .
Figure 3. Constructed block grid: (a) calculated grid in the horizontal section, (b) calculated grid in the vertical section.

Figure 4 .
Figure 4.The location of sources and receivers.

Figure 7 .
Figure 7.Comparison of two calculations (calculation 1-there is no layer of adipose tissue around the liver; calculation 2-the liver is surrounded by a layer of adipose tissue 0.01 m thick).

Figure 8 .
Figure 8.Comparison of pressure in receiver 6.

Figure 9
Figure 9 represents the quantitative difference of the registered pressure for the two experiments.

Figure 9 .
Figure 9.The registered pressure for the two experiments.

Figure 10 .
Figure 10.Comparison of two calculations (experiment 1-there is no layer of adipose tissue around the liver; experiment 2-the liver is surrounded by a layer of adipose tissue at a sound speed of 1460 m/s; experiment 3-the liver is surrounded by a layer of adipose tissue at a sound speed of 1000 m/s.

Table 1 .
Acoustic parameters of the human body.