Nonlinear shear wave propagation to assess biomechanical properties in soft
tissueJuan MELCHOR(1,2,3) , Guillermo RUS(1,2,3) , Antonio CALLEJAS(1,2) , Inas FARIS(1,2) , Javier NARANJO-PÉREZ(4) , Miguel
RIVEIRO(1)

Abstract
Ultrasonic-based characterization of soft tissue is an emerging technology with great potential as a clinical diagnostic tool. There is evidence that abnormalities in the structural architecture of soft tissues are linked to abroad range of pathologies including tumors, liver fibrosis, preterm birth...The present work derives different
theories of non-linear formulation of wave propagation using generalized and particular coordinates on a material characterized by different strain energy functions. Afterwards, these formulations are adapted to simulate
the propagation of shear waves on nonlinear elastic materials to explain the properties of soft tissue when ismeasured under these conditions. The resultant system of equations could be solved under a finite differences
framework and probabilistic inverse problem to explore how improves the understanding of soft tissues. The
response of the system varying the linear/non-linear parameters of the strain energy is obtained via Finite Differences Time Domain simulation showing significant variations among the different models. Nonlinear wave
propagation approach implies a new paradigm to model the nonlinear biomechanical behavior of soft tissue
improving the plausibility versus previous elastic or viscoelastic predictions.
Keywords: Sound, Insulation, Transmission

1

INTRODUCTION

The nonlinear elastic field has introduced a great interest in the study of hard and soft tissues mechanics in thelast decades [17, 16]. The research in shear and longitudinal waves propagation in tissue starts at the beginningof the second part of the twentieth century. The ultrasonic imaging in soft tissues begins when the use of
dynamic tests to ultrasonically estimate the compressibility and mobility of breast tumors was analyzed [25].
Ultrasonic shear elastography have been applied by in order to measure the linear and nonlinear elastic responseof hyperelastic tissues. The Static Elastography, originally proposed by Ophir et al. [18], has been successfully
used to measure nonlinear properties of vascular tissues or malingnant and beningn tumors of breast tissues[7, 9]. Nevertheless, the results of this method are not satisfactory for organs or tissues which are deep ordifficult to compress as the strain profile may be uncertain. The Transient Elastography introduced by Stefan
Catheline [4] has been used, for example, to diagnose cirrhosis [5] and for assessment of hepatic fibrosis [21].
The Supersonic Shear Imaging (SSI) technique, presented in 2004 by Bercoff and Tanter [2], has also beenused to obtain the third and fourth order constants, A and D, respectively, of pig brain tissue or human breast
tissue. However, these studies were limited to ex vivo experiments [6, 20, 9]. Nowadays, the consideration of
torsional waves (they propagate radially and in depth) is a potential key due to the limitations of shear andcompressional waves [12]. Firstly, this typology of wave can propagate by quasi-compressible material and,it is more sensitive to quantify the consistency changes, for example, in the case of tumors [23]. Secondly,
the variations of biomechanical parameters are more sensitive in the regime of low energy where this wave isgenerated. Finally, torsional movement does not generate secondary interfering P-waves at the boundary of thetransducer where pure shear waves are difficult to create [12, 11]. A torsional ultrasonic transducer has beenused to measure nonlinear parameters of ligament tissue and the shear modulus of cervical tissue in pregnant
women [3].
The aim of this work is to elaborate a set of nonlinear formulations of a torsional and shear wave propagating
on a hyperelastic tissue defined from the Hamilton?s strain energy function and others [26, 13]. While thequadratically nonlinear propagation of the torsional wave on a hyperelastic material has been studied [19], the647consideration of new nonlinear terms in the equations is an open issue to consider experimentally en the future.

2

MATERIAL AND METHODS

2.1 Fundamentals
The Green-Lagrange strain tensor that governs the elasticity law is defined in terms of displacements in index
notation as,1
εi j = (ui, j + u j,i + uk,i uk, j )(1)2
where the third term is related with geometrical nonlinearity. The physical nonlinearity is here focused on
hyperelastic materials, in which the strain energy function W is defined per unit reference (undeformed) volume
and acts as potential of the stress. The strain energy function is defined following the expression by Landau
and Lifshitz [10],
A
C
λ
(2)
W = I12 + µI2 + I3 + BI1 I2 + I13
2
3
3
where µ and λ are the Lame constants and A, B and C are the TOEC. I1 , I2 and I3 are the invariants of theGreen-Lagrange strain tensor defined by Eringen et al. in 1974 [15],I1 = trε = εii , I2 = trε 2 = εi j ε ji , I3 = trε 3 = εi j ε jl εli(3)

The expansion to fourth order of the energy density is necessary when the nonlinear effects in shear waves
are considered. This is because for incompressible media, nonlinearity at third order is missing in the particle
displacement [8]. This expansion is characterized by four new terms, yielding,
W=

λ 2
A
C
I + µI2 + I3 + BI1 I2 + I13 + DI22 + EI1 I2 + FI12 I2 + GI14
2 1
3
3

(4)

where D, E, F and G are the Fourth Order Elastic Constants (FOEC).
When the case of quasi-incompressible soft tissue media is analyzed, where λ >> µ, the strain energy function
is simplified to,A(5)
W = µI2 + I3 + DI22
3
where A and D are the third and fourth order elastic constants, respectively [8].
The motion equation in the contravariant basis in the Lagrangian configuration is obtained in general coordinates
as,
i
h

∂ uk
∇i δ ji + ∇ j ui S jk = ρ 2(6)∂tThis is the compact form of the momentum equation expressed in the Lagrangian configuration. To expand theequation, several previous steps are required. The tensor including these derivatives multiplied by S jk gives asa result another second kind tensor, ∇ j ui · S jk = C ik . Hence, Equation 6 can be reformulated as the divergenceof a tensor B ik which results,∂ B ik
∂ uk
im k
mk i
+
B
Γ
=
ρ
+
B
Γ
(7)
mi
mi
∂ξi
∂t 2
Note that for this particular case, the general coordinates are ξ 1 = r, ξ 2 = θ and ξ 3 = z, the Christoffel symbols
the configuration depends only on the velocity of the angular coordinate θ .
From other point of view, following where the series expansion concept put forth by Landau [10], only volumetric part is detailed in terms of the nonlinear acoustic parameter β as was explored recently [13].
There are four combinations of nonlinear acoustic parameters β that may explain a different scenario of experimental calculations as in this case where the exploration of quasi-fluids nonlinearity is considered. The concept648

of the non uniqueness for the classical acoustic parameter β have been recently studied by Bender [1] but without a physical explanation. The aforementioned combinations could be expanded as exploring the whole set ofcombinations by quadratic terms as follows,
σi j = −3Kvδi j + 2µDi j −3η v v̇δi j + 2η Ḋi j + +9Kβ vp v2 δi j + 9Kβ d p Dkp D pk δi j + 4µβ ds Dik Dk j + 4µβ cs vDi j| {z } | {z } | {z } | {z } |
{z
} |
{z
}
pressure|pressure

shear

{z

σiLj (Linear)}|pressure

shear

{z

σiVj (Viscous)

(8)

shear

} |

{z

}

σiNL
j (Nonlinear)where v and Di j are the volumetric and deviatoric parts of the stress respectively. Four nonlinear parametersof first order have been defined as β vp , β d p as the volumetric, the deviatoric and the compound denoting thepressure components and β ds , β cs as the deviatoric, and compound denoting the shear components. However inthe case of cuasi-compressible, when compressibility is much higher than shear moduli, it means K >> µ whenthe volumetric nonlinear terms could be neglected,
σi j = −3Kvδi j + 2µDi j −3η v v̇δi j + 2η Ḋi j + 9Kβ d p Dkp D pk δi j + 4µβ ds Dik Dk j
| {z } | {z } | {z } | {z } |
{z
} |
{z
}
pressure|pressure

shear

{z

σiLj (Linear)}|pressure

shear

{z

σiVj (Viscous)

} |

(9)

shear

{z

σiNL
j (Nonlinear)}With the aim to find a relationship between this nonlinear expansion of beta parameters and TOEC, we assuming
now that the strains are separated in volumetric and deviatoric part to the second order, so ε yields,
εik εk j = Dik Dk j ; εi j = Di j ; εk j = Dk j ; εkp = Dkp ; ε pk = D pk
By making use of Cauchy stress described in Equation 9, in nonlinear regime, an equivalence is deducted interms of Third Order Elastic Constants TOEC,σiNLj= Dik Dk j (A + 4µ) + Dkp D pk Bδi j(10)Since the relationship is established in the nonlinear constitutive equation, nonlinear acoustic parameters of firstorder are explicitly deducted as follows in terms of Third Order Elastic Constants,
β dp =

A + 4µB; β ds =
9K
4µ

(11)

2.2 Experimental setup
A schematic view of the experimental setup is shown in Figure 1. The transmitter and receiver are located in
the same transducer as was explained in the previous sections. The measurements was taken as in transmission,
with guarantees perfect alignment along the propagation direction of sound. The response signals were sampledwith a high resolution A/D converter with amplification Fonestar without preamplification stage and saved forpost-processing. The transmitter was a torsional transducer manufactured in our END lab at University ofGranada.

649

to the resulting average of 490 captures of the signal, providing an eﬀective
reduction of noise for the detected response signal, increasing the signal-tonoise ratio around 25 dB. Figure 3.2 depicts the experimental setup used torecord the ultrasonic signals.
Wave generator
Synchronized signal- Acquisition card
A/D converter-

6A
Computer

Excitation signal

ReceiverSpecimen

AmplifierH- H

Response signal

Transmitter-

3.3. Semi-analytical approach

Figure 1. Experimental configuration for a nonlinear torsional transducer
The measurements procedure consists of four steps. Initially, the response

signal was measured in water at two diﬀerent positions corresponding to

We used an arbitrary generator (Agilent 33500) to drive the transmitter at different frequencies. The ultrasoniclocations ahead of and at the back of the specimen (steps 1 and 2 ).
excitation frequencies were arbitrarily chosen. Different frequency combinations were used: f1 = [600, 800, 1000]
Then,
the signal
response
were measured
the specimen
situToalong
Hz. The amplitude of the
input
atsignals
the frequency
was with
increasing,[5, 10]in V.
avoid unwanted interferencethescanning
area
at
both
frequencies(steps3
and4).
effects in the propagation medium, even the risk of burning, we used short trains of pulses for excitation and a 8
ms pulse duration and 3 cycles.
investigatedgeneratedwaves
as a function
ofofdistance. The transmitterIn orderWeto reconstruct
the thelinearand nonlinear
material’s
properties
receiver distance was a total distance of 2.43 mm. At the distance x from the transmitter, the pressure of thethe specimen under inspection, a semi-analytical
approach is proposed by(0)(1)
fundamental, u1 , and that of the generated harmonics, u1 , were determined, for the different frequencies andconsidering the following assumptions: (i) the attenuation in the water layers
excitation levels.
is negligible, (ii) the thickness a3 of water layer 2 is small, and thus thenonlinearity can be assumed not to accumulate over this distance (i.e.

(3)
w ⇡

0) [10], (iii) only second harmonics are considered, and thus frequency-mixing
12

Figure 2. Connective tissue and silicone mold during measurements and oscilloscope view
The material used in this experiment was the silicone mold as a first step. Secondly, a connective tissue was
explored with the transducer prototyped as is shown in Figure 2, with the characteristic amplitudes as in thelongitudinal waves.

3

RESULTS

By adopting the acoustic nonlinear constitutive equation presented above Equations 9 in terms of deviatoric a
volumetric parts, is possible to establish the three dimensional nonlinear equation of motion up to first-ordernonlinearity in terms of four parameters β . This formulation, implies that beta parameter is not unique and canbe defined as the separation between pressure and shear waves and it is solved under perturbation scheme as,
Ba c3sx1 A21 ω1650= 3(9β d p K + 2β ds µ)(12)This results can be compared with the Zarembo ones [26].


1 A
µ
+
+B
βZτ =
ρ0 ρ 2(13)

So, converting β d p and β ds in the corresponding TOEC is posible to obtain a similar expression as,βT =(B + A /2 + 2µ)4(3K + 4µ)(14)

In order to validate experimentally the theoretical background, is necessary to design an experimental procedure, for this purpose these results should be convert into a pressure magnitudes. Then, by making use of theconversion to amplitudes of fundamental and first harmonics into transversal stresses, it is posible to calculate
the nonlinear shear parameters in terms of transversal stresses.
Then, torsional beta, β T could be defined by in terms of pressures and disaggregated in two parts, one due tothe liquid or water part and the other referred to collagen or fiber undulation part,
(1)

βT =

ρc3s T12(B + A /2 + 2µ)=−(0)4(3K + 4µ)π f1 x(T )2

(15)

12

where the relation 2π f = ω1 has been introduced and the relationship between pressure and displacements.
3.1 Experimental results
To obtain the experimental results in the case on torsional nonlinearity coefficient, is necessary to use thefollowing equation calculated from PZT-5 Piezoelectric material where elastic and electric field are coupled,
T12 = CE44 S12 + e24 E2(16)Neglecting elastic field to found a conversion factor between voltage and Pascal, 122.9 Pa was the multiplication
factor for the case of 10V and silicon mold. The conversion factor is formulated from the analytical simplified
model where piezoelectric and design material effects were considered interacting with two layers of tissue. For
this reason, the stress on the piezoelectric ceramic is approximately the stress on the tissue by a correction
factor αc ,piezotissueT12 = αc T12(17)
Parallel to this, to obtain αc transmission coefficient should be calculated, by making use of the energy of thetransmitter and receiver.
4z1 z2energy,tT==(18)(z1 + z2 )2energy,r
where energy,t is the energy of the transmitter and energy, r is the energy of the receiver. Knowing the different
layers of the receiver, the stress and displacement ratios could be deducted as,
Tt = Ti

2zi
2zt
or ut = uizi + zt
zi + zt(19)

The final experimental results are detailed in Table 1 and in Figures 3 where nonlinear torsional acoustic parameter β T was extracted for silicon mold, ligament tissue and liver tissue through the speed of sound of shear
waves cs , density ρ, shear modulus µ and transmission coefficient Z.

651

SampleEnergy
(V)
10
10
10

Silicon moldConnective tissue
Liver tissue

Frequency
(Hz)
800
800
800

Z

cs
(m/s)
13.2
90
4

ρ
(kg/m3 )
1100
1000
1000

βT

µ
(KPa)
0.16
8.1
0.016

13200
90000
4000

-400±3-8000±0.0-88±20

Table 1. Nonlinear torsional results for input frequency and amplitude f1 = 800 Hz and A = 5, 10 V, respectively,
for silicone mold, connective and liver tissue.
Signal [kPa]

Signal [kPa]

0.02
0.01
0
-0.01
-0.02
-0.03
-5

0

5

10

0.15

6

0.1

4

Signal [kPa]

0.03

0.05
0
-0.05
-0.1
0

5

Time [µ s]

-8
-5

15

4

0

-2

10 -4
0.5

1

1.5

2

2.5

5

3

3.5

4

4.5

10

15

Time [µ s]

10 2
10 0
10 -2
10 -4
0

0

0

10 4

Pressure [Pa]

10

10

10 4

10 2
10

-4

Time [µ s]Pressure [Pa]Pressure [Pa]

10

0
-2

-6

-0.15
-5

15

×10 -3

2

0.5

1

1.5

5

2

2.5

3

3.5

4

4.5

5

10 2
10 0
10 -2
10 -4
0

0.5

1

1.5

Frequency, f, [kHz]Frequency, f, [kHz]

2

2.5

3

3.5

4

4.5

5

Frequency, f, [kHz]Figure 3. Nonlinear torsional signal and harmonics from silicon, connective and liver tissue for 10V.
3.2 Numerical results
The difficulty of finding the close form solution makes it necessary to the numerical method Finite DifferenceTime Domain [22]. The boundary conditions are expressed in terms of velocities and stresses and geometricaland mechanical symmetries have been taken into account. The displacements are permitted without any constraint and no stresses are created. In these boundaries the components of the second Piola-Kirchhoff stresstensor are nil. The velocity is u̇θ = 0. The propagation of the signal cannot be simulated indefinitely along thespace. An issue arises when analyzing open space regions where the simulation domain must be limited. In
this work, to ensure that the boundary does not reflect the wave, an Absorbing Boundary Conditions ABC has
been considered in those boundaries where the stress is equal to zero and the wave can be reflected.
It can be noticed that the nonlinear curve does not have a purely sinusoidal behavior as the linear curve does.
Conversely, the ridges and valleys are slightly sloping, which is normal because of its nonlinear nature. The
Fast Fourier Transform (FFT) algorithm has been used to obtain the solution in the frequency domain. The main
remark about the frequency domain nonlinear response is that only odd harmonics are relevant. For instance,
for an excitation frequency of 500 Hz, the first (500 Hz), the third (1500 Hz) and the fifth (2500 Hz) harmonics
are clearly observable [14].
50nonlinearlinear

Displacement u θ [µm]Displacement u θ [µm]Displacement u θ [µm]

10
5
0
-5
nonlinearlinear

-10

5

0

-5

nonlinearlinear

-5
0

5

10

15
Time [ms]

20

25

0

b)

4

6
Time [ms]

8

10

0

c)

10-6
Power spectrum

Power spectrum

10-6

2

10-8

10-10

nonlinearlinear2

4

6Time [ms]

8

10

10-6
Power spectrum

a)

10-8

nonlinearlinear

10-10

10-8

nonlinearlinear

10-10

-12

10

0

5

10

15

20

0

1

2

3

4

5

6

0

5

10

15

20

Frequency [kHz]Frequency [kHz]Frequency [kHz]a)b)c)Figure 4. Displacement uθ at the receiver (left column) and Fourier transform (right column) for a frequency
of: a) 500 Hz, b) 1500 Hz and c) 2000 Hz. The parameters of the hydrogel are: A = 40 kPa and D = 3000 kPa.

652

4

CONCLUSIONSIn this work, a set of formulations haven been developed from classical theories of nonlinear elasticity andits derivation to cylindrical coordinates of the propagation of the torsional wave on a hyperelastic material.
For this reason, a breakthrough has been achieved in the field of ultrasound.The nonlinear classical acoustics
extension may be also applied to derive a solution that separates the nonlinear terms that are interpreted as
nonlinearity from the matrix and the nonlinearity from the fibers. This is an important conclusion with a mainnovelty in tissue microstructure. The results suggest that these nonlinear terms that from now can be obtainedat real time describe the behavior of the tissue in a new scale. Although this works have contributed in some
way to advancing the understanding of ultrasonic nonlinearity, there are several aspects to be studied in the
future with potential impacts on biomedical research and understanding the mechanics of tissues. However, soft
tissues should be modeled as transverse isotropic, orthotropic or, more accurately, anisotropic materials. Anotherimportant consideration is the viscosity of soft tissue both in linear and nonlinear regime.

ACKNOWLEDGEMENTS
This research was supported by Ministry of Education DPI2017-83859- R, DPI2014-51870-R, and UNGR15-CE3664, Ministry of Health PI16/00339 Carlos III/FEDER and Junta de Andalucia PI-0107-2017, PIN-0030-2017.

REFERENCES[1] F. A. Bender, J.-Y. Kim, L. J. Jacobs, and J. Qu. The generation of second harmonic waves in an isotropic
solid with quadratic nonlinearity under the presence of a stress-free boundary. Wave Motion, 50(2):146–161,
2013.
[2] J. Bercoff, M. Tanter, and M. Fink. Supersonic shear imaging: A new technique for soft tissue elasticity
mapping. IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 51(4):396–409, 2004.
[3] A. Callejas, A. Gomez, J. Melchor, M. Riveiro, P. Massó, J. Torres, M. López-López, and G. Rus. Performance study of a torsional wave sensor and cervical tissue characterization. Sensors, 17(9), 2017.
[4] S. Catheline, F. Wu, and M. Fink. A solution to diffraction biases in sonoelasticity: The acoustic impulse
technique. The Journal of the Acoustical Society of America, 105(5):2941–2950, 1999.
[5] N. Ganne-Carrié, M. Ziol, V. de Ledinghen, C. Douvin, P. Marcellin, L. Castera, D. Dhumeaux, J. Trinchet,and M. Beaugrand. Accuracy of liver stiffness measurement for the diagnosis of cirrhosis in patients withchronic liver diseases. Hepatology, 44(6):1511–1517, 2006.
[6] J. Gennisson, M. Rénier, S. Catheline, C. Barrière, J. Bercoff, M. Tanter, and M. Fink. Acoustoelasticity
in soft solids: Assessment of the nonlinear shear modulus with the acoustic radiation force. Journal of the
Acoustical Society of America, 122(6):3211–3219, 2007.
[7] S. Goenezen, J. Dord, Z. Sink, P. E. Barbone, J. Jiang, T. J. Hall, and A. A. Oberai. Linear and nonlinear
elastic modulus imaging: An application to breast cancer diagnosis. IEEE Transactions on Medical Imaging,
31(8):1628–1637, 2012.
[8] M. F. Hamilton, Y. A. Ilinskii, and E. A. Zabolotskaya. Separation of compressibility and shear deformation in the elastic energy density (L). Journal of the Acoustical Society of America, 116(1):41–44, 2004.
[9] Y. Jiang, G. Li, L. Qian, S. Liang, M. Destrade, and Y. Cao. Measuring the linear and nonlinear elastic
properties of brain tissue with shear waves and inverse analysis. Biomechanics and modeling in mechanobiology, 14(5):1119–1128, 2015.

653

[10] L. Landau and E. M. Lifshitz. Theory of elasticity, volume 7 of Course of Theoretical Physics. Pergamon
Pr, 1 edition, 1959.
[11] J. Melchor, R. Muñoz, and G. Rus. Torsional ultrasound sensor optimization for soft tissue characterization.
Sensors, 17(6):1402, 2017.
[12] J. Melchor and G. Rus. Torsional ultrasonic transducer computational design optimization. Ultrasonics,
54(7):1950–1962, 2014.
[13] R. Muñoz and J. Melchor. Nonlinear classical elasticity model for materials with fluid and matrix phases.
Mathematical Problems in Engineering, 2018, 2018.
[14] J. Naranjo-Pérez, M. Riveiro, A. Callejas, G. Rus, and J. Melchor. Nonlinear torsional wave propagationin cylindrical coordinates to assess biomechanical parameters. Journal of Sound Vibration, 445:103–116,
2019.
[15] A. Norris. From nonlinear acoustics, ed. M.F. Hamilton and D.T. Blackstock (academic press, san diego,
1998), chap. 9, pp. 263-277.
[16] A. Oberai, N. Gokhale, S. Goenezen, P. Barbone, T. Hall, A. Sommer, and J. Jiang. Linear and nonlinear elasticity imaging of soft tissue in vivo: demonstration of feasibility. Physics in Medicine & Biology,
54(5):1191, 2009.
[17] R. W. Ogden. Nonlinear elasticity, anisotropy, material stability and residual stresses in soft tissue. Springer,
2003.
[18] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li. Elastography: a quantitative method for imaging
the elasticity of biological tissues. Ultrasonic imaging, 13(2):111–134, 1991.
[19] J. Rushchitsky and Y. Simchuk. Quadratic nonlinear torsional hyperelastic waves in isotropic cylinders:Primary analysis of evolution. International Applied Mechanics, 44(3):304–312, 2008.
[20] M. Rénier, J. . Gennisson, M. Tanter, S. Catheline, C. Barrière, D. Royer, and M. Fink. Nonlinear shear
elastic moduli in quasi-incompressible soft solids. In Proceedings - IEEE Ultrasonics Symposium, pages
554–557, 2007.
[21] L. Sandrin, B. Fourquet, J. . Hasquenoph, S. Yon, C. Fournier, F. Mal, C. Christidis, M. Ziol, B. Poulet,
F. Kazemi, M. Beaugrand, and R. Palau. Transient elastography: A new noninvasive method for assessmentof hepatic fibrosis. Ultrasound in Medicine and Biology, 29(12):1705–1713, 2003.
[22] J. Schneider. Understanding the finite-difference time-domain method. School of electrical engineering andcomputer science Washington State University., 2010.
[23] F. Sebag, J. Vaillant-Lombard, J. Berbis, V. Griset, J. Henry, P. Petit, and C. Oliver. Shear wave elastography: a new ultrasound imaging mode for the differential diagnosis of benign and malignant thyroid
nodules. The Journal of Clinical Endocrinology & Metabolism, 95(12):5281–5288, 2010.
[24] H. von Gierke, H. L. Oestreicher, E. Franke, H. Parrack, and W. von Wittern. Physics of vibrations inliving tissues. Journal of applied physiology, 4(12):886–900, 1952.
[25] P. Wells and H. Liang. Medical ultrasound: imaging of soft tissue strain and elasticity. Journal of the Royal
Society Interface, 8(64):1521–1549, 2011.
[26] L. Zarembo and Krasilnikov. An Introduction to Nonlinear Acoustics. Moscow, 1966.

654