流体中非线性超声波数值仿真有限差分

1、A numerical formulation for nonlinear ultrasonic wavespropagation in uidsa,*baESCET,Universidad Rey Juan Carlos,Tulip a n s/n,28933M o stoles,Madrid,SpainbInstituto de Ac u stica,CSIC,Serrano 144,28006Madrid,SpainReceived 17November 2003;received in revised form 15January 2004;accepted 10February 20

2、04Available online 24February 2004AbstractA nite-dierence algorithm is developed for analysing the nonlinear propagation of pulsed and harmonic ultrasonic waves in uid media.The time domain model allows simulations from linear to strongly nonlinear plane waves including weak shock.Eects of absorptio

3、n are included.All the harmonic components are obtained from only one solving process.The evolution of any original signal can be analysed.The nonlinear solution is obtained by the implicit scheme via a fast linear solver.The numerical model is validated by comparison to analytical data.Numerical ex

4、periments are presented and commented.The eect of the initial pulse shape on the evolution of the pressure waveform is especially analysed.Ó2004Elsevier B.V.All rights reserved.Keywords:Nonlinear propagation;Numerical modelling;High-power ultrasoundIn this paper we propose a new numerical algor

5、ithm for studying the nonlinear propagation of transient and periodical signals in uid media.The work is motivated in many applications where high intensity waves,which can not be described by linear laws,are involved (industrial applications of high-power ultrasound,acoustic imaging,biomedical rese

6、arch (test and ther-apy,etc.Several works dealing with analytical and numerical studies of propagation of nonlinear waves can be found in the literature.In particular perturbation methods have obtained second order solutions without geometrical limitations but only in the case of strongly limited wa

7、ve-amplitude 1,2,and two-dimensional models have reached the simulation of nite amplitude wave propagation 3.In Chapter 11of Ref.4a review of computational methods applied to nonlinear propa-gation of acoustic waves is given,and numerical models in the time and frequency domains are commented.These

8、approximations are based on Eulerian coordi-nates and only second order terms of the Mach number are considered in equations.Some authors have used Eulerian coordinates and the retarded timevariable associated with the propagation in the þx direction:t Àx =c 0,which allows them to reduce b

9、y one the order of the dierential equation for wave motion 4,5.The purpose of the present work is to develop a numerical method for studying the nonlinear propaga-tion of plane waves by using the full nonlinear equation derived in Lagrangian coordinates,without any restric-tion about the value of th

10、e acoustic Mach number.The formulation is written in the time domain.Natural spatial and time coordinates are used.Up to the knowledge of the authors,no works exist solving this equation by using natural coordinates.These imply the need of imposing a nonreecting boundary condition.The analysis of th

11、e waveform evolution for any original signal is possible:periodic excitation,Gaussian,rectan-gular pulses,etc.In addition,all the harmonic compo-nents are obtained by only one solving process,with the consequent save in computation time.The solution of the nonlinear problem is obtained via the devel

12、opment of an implicit nite-dierence scheme and by using a fast*Corresponding author.Tel.:+34-91-664-74-82;fax:+34-91-488-73-38.E-mail address:c.vanhilleescet.urjc.es (C.Vanhille.0041-624X/$-see front matter Ó2004Elsevier B.V.All rights reserved. Ultrasonics 42(2004 11231128linear solver.The equ

13、ations of the problem are given in Section 2.Section 3describes the numerical model.In Section 4some numerical experiments are carried out.The conclusions of this work are analysed in the last section of the paper.2.Governing equationsThe full nonlinear one-dimensional wave equation can be derived b

14、y using Lagrangian coordinates.In the case of an ideal gas 6,the following dierential equa-tion is obtained:q 0o 2u o t 2¼c p 011þo uo xc þ1o 2u o x 2þq 0m b o 3u o t o x 2ð1Þwhere p 0is the ambient pressure,c is the specic heat ratio,q 0is the initial state density,u i

15、s the displacement,m is the kinematic shear viscosity,and b is the so-called viscosity number.t and x are,respectively,the time and one-dimensional spatial coordinates.No approxima-tions have been made about the acoustic Mach number value or about the attenuation parameter.The only limitations on pr

16、essure amplitude in the model are those derived from the isentropic approximation 6.However,since we consider the propagation of the wave within an unbounded domain,even in the case of high acoustic Mach number,the isentropic property of the uid can not be questioned.The acoustic pressure waveform a

17、nd distribution are evaluated by using the following expression 6:p ¼p 01þo u o xc Àp 0ð2ÞWe consider progressive plane waves and a source at x ¼0.The uid is assumed to be initially at rest,i.e.,particle displacement and velocity null at t ¼0.The following auxiliar

18、y conditions are written:x ¼0u ð0;t Þ¼f ðt Þx ¼L c 0o u o xðL ;t Þ¼Ào uo t ðL ;t Þt ¼0u ðx ;0Þ¼0o u ðx ;0Þo t¼08x ¼08<:ð3Þwhere c 0is the small-amplitudes value of the sound speed,L

19、is the length of the study domain and f ðt Þis the excitation of the medium dened as a function of time.The full dierential equation (1is solved by using nat-ural space and time variables (x and t ,and the pro-gressive character of the wave is imposed by proposing the nonreecting boundary

20、condition at x ¼L (Eq.(3.3.Numerical formulationThe nite-dierence formulation is developed by considering the dimensionless independent variables X ¼x =k and T ¼x t in Eqs.(1(3,where k is the wavelength of the signal and x is its angular frequency.The X T space is discretized by means

21、 of a uniform grid dened by the steps h in space and s in time.O ðh 2;s 2;h Þnite dierences treat the operators appearing in Eq.(16.O ðh ;s Þnite dierences are developed for Eqs.(2and (3.The discretization leads to a dierence equation for each space grid point.The nonreecting bou

22、ndary condition at x ¼L implies the construction of a new numerical scheme.By setting,A ¼m b s 2xk h 2;B ¼c 20s k x 2h 2;C ¼1k h ;D ¼Àc 0s xk h ;E ¼c 0sxk hþ1and considering u m ;n as the value of the displacement u at the point m ;n of the grid corresponding

23、to the space point m and the time point n ,and M as the number of space points,we obtain the following system:u m ;n ¼u 0f ðn Þm ¼1ð4a ÞÀAu m À1;n þð2A þ1Þu m ;n ÀAu m þ1;n¼ÀAu m À1;n À2þð2A À1&#

24、222;u m ;n À2ÀAu m þ1;n À2þ2u m ;n À1þB u m þ1;n À1ðÀ2u m ;n À1þu m À1;n À1ÞÂ1½þC ðu m ;n À1Àu m À1;n À1ÞÀc À1m ¼2;.;M À1ð4b ÞDu m À1;n 

25、4;Eu m ;n ¼u m ;n À1m ¼M ð4c Þwhere f ðn Þ¼sin ððn À1Þs Þin the case of harmonic excitation of amplitude u 0at the driven frequency f ,orf ðn Þ¼e Àx 2b ððn À1Þsx Àt 0Þ2cos ððn

26、 À1Þs Þwhen working at the source with a pulse of amplitude u 0,driven frequency f ,width band x b ,and centre time t 0.At each time step,the (M ,M set of algebraic equations has to be solved to obtain the grid point values of u .The scheme is implicit.The following economic and fast

27、method,based on a LU decomposition 7valid for the whole j th period,is proposed.Excepting the rst and last equations,the matrix of the system is 3-diagonal,and the Thomas method can be used 8.Moreover the matrix is symmetric,and the method can be improved.The rst and last equations need a specic tre

28、atment.From now on we put M À2¼n .The method is now described by considering a linear (n ,n set of equations Ax ¼b of the form ða b a Þ.We store the elements of the diagonal of A into the n -vector U þ¼ðb b .b b ÞT and a new n -vector L Àis created:L

29、 À¼ð00.00ÞT .A LU-decomposition is now applied by Gaussian elimi-nation without pivoting:1124 C.Vanhille,C.Campos-Pozuelo /Ultrasonics 42(200411231128L Àðk þ1Þ¼a U þðk ÞU þðk þ1Þ¼Àa L Àðk þ1Þ

30、;þb(k ¼1;.;n À1ð5ÞVectors U þand L Àrepresent,respectively,the diagonal of the upper matrix and the lower diagonal of the lower matrix.They constitute the decomposition of the matrix and are used at each time step (A is always the sameby solving two 2-diagonal syst

31、ems.The lower system is solved by forward substitution and the upper one by back-substitution,and the result is saved into a created and initialised vector z :z ð1Þ¼b ð1Þz ði Þ¼b ði ÞÀL Àði Þz ði À1Þi ¼2;.;nz &

32、#240;n Þ¼z ðn ÞU þðn Þz ði Þ¼z ði ÞÀa z ði þ1ÞU þði Þi ¼n À1;.;1ð6ÞThe content of z is then stored on disk before the fol-lowing time step.In addition to the n -vector U þ,scalars

33、a and b ,and the vector b ,this solver only creates two new n -vectors:L Àand z .The operations number is 7n À6for the rst time step and 5n À4for each fol-lowing time step.A Von-Neumann analysis shows the consistent scheme to be conditionally stable,and an optimised relation between h

34、 and s is dened and used to obtain stability 6.The nonreecting condition at L requires an additional analysis of stability by the von Neumann method.The error e i b mh n n at the point m ;n is introduced into Eq.(4c,where b is the frequency of the error and n is the amplication factor.This study lea

35、ds to the inequality:c 0s xk h ð1h Àcos ðb h ÞÞþ1i 2þc 0s xk h2sin 2ðb h ÞP 1:This relation is veried when c 0s xk h P À1,which is al-ways true.Eq.(4cis thus unconditionally stable and does not generate any additional convergence problem.In addition

36、to the creation of the new nite-dierence scheme,the consideration of the propagation of the waves into an open eld domain generates several problems that must have been imperatively solved to make the simulations feasible.First,simulations of strongly nonlinear waves prop-agating in a quite large sp

37、ace implies very high storage and time costs.These aect the possibility of using the model.Computing costs are reduced by sliding and expanding the grid as the disturbance propagates,by taking into account the null displacements in the unperturbed zone (mechanical perturbations reach the i th wavele

38、ngth after i periods.By acting period by period during the whole nite-dierence process,it is possible to avoid their calculation.To this purpose the j th part of the process is realised during the j th periodby considering (j þ1wavelengths and the last two time steps of the previous (the (j

39、92;1thperiod as initial condition of the current (the j thperiod,completed by zero values for the (j þ1th wavelength.This sequence is realised from the rst part of the process (rst period with initial conditions u m ;1¼0and u m ;2¼u m ;1,up to the NP th part of the process (the last N

40、P th period.Fig.1shows a visual representation of the computed values in the X T space.When working with transient signals,each part of the process is realised over a specic number of periods.On the other hand,if a simulation has been carried out without a sucient number of periods,this can be start

41、ed again to increase the scope of the study without calculating again from the beginning.Large CPU time and storage space are thus saved.In this case the initial conditions are the last two time steps of the previous study,completed by zero values for the (NP þ2th wavelength.Once the displaceme

42、nt is known at the grid points,the acoustic pressure values are evaluated at each time step via a classic O ðh Þnite-dierence formula from Eq.(2.4.Numerical experimentsIn this paragraph a number of numerical experiments are presented.All the results are obtained from con-vergent simulation

43、s with a suitable number of spatial points.The numerical model is rst compared,for nite but moderate amplitude waves,to an analytical perturba-tion model previously developed 1.A harmonic source (f ¼20kHz,u 0¼25l mis considered in air.Fig.2displays the analytical and numerical (obtained by

44、 FFTsound pressure distribution of the rst and secondFig.1.Part of the X T space computed during the j th part of the process.C.Vanhille,C.Campos-Pozuelo /Ultrasonics 42(2004112311281125harmonics.The agreement validates the numerical model,as well as the proposed boundary condition at x¼L.A str

45、ongly nonlinear continuous wave in air is now considered(f¼20kHz,u0¼300l m.Fig.3presents the spatial distribution of the fourrst pressure har-monic components obtained by FFT,as well as a linear evaluation of the pressure amplitude.After the initial increase of distortion,a stronger decrea

46、se of the higher frequencies occurs with the distance from the source, giving a more and more sinusoidal signal.The domi-nance of the nonlinear attenuation becomes clear.Application of the numerical tool to the propagation of transient signals is now shown and physical results commented.The source f

47、unction is fðtÞ¼u0eÀ½x bðtÀt0Þ 2cosðx tÞ.A very short(x B¼5Â106,t0¼10À6and a Gaussian(x B¼106,t0¼3Â10À6pulses are chosen.The evolution of wave-shapes and shock formation is studied and the importance of the ini

48、tial wave-shape is analysed.Auid which acoustic proper-ties are similar to those of tissues of the body(except lung,bone,and fatis now considered:c0¼1500m/s, c¼6:2,q0¼1000kg/m3,and a¼11mÀ14.We con-sider the central excitation frequency f¼1MHz,quite usual in medical appl

49、ications.Fig.4shows the evolution of the wave-shape and frequency decomposition for the two cases.The displacement amplitude at the piston is u0¼1:5l m,i.e.,the order of the initial pressure amplitude is15MPa.We observe that the strong har-monic distortion occurs(like for periodical wavesat the

50、 1126 C.Vanhille,C.Campos-Pozuelo/Ultrasonics42(200411231128rst wavelengths from the source in both cases.For the short pulse,the central frequency decreases with the distance to the source.The order of the pulse amplitude is27%of its initial value at7.2cm from the source and its frequency is about

51、three times less than the original; no much harmonic distortion aects this state of the pulse.Thus,the main eect involved is the nonlinear attenuation associated with the harmonic distortion in the nonlinear medium.When dealing with a wider Gaussian pulse,the strong distortion aects as well at the r

52、st wavelengths from the source.The order of the wave amplitude at7.2cm from the source is32%of its original amplitude.The central frequency decrease is only12.5%, in front of300%of the short pulse.This decrease can be interpreted by the stronger presence of low frequencies in the short pulse,which a

53、re quite less attenuated.Another important dierence is the apparition of a low frequency for this kind of signals.This low frequency increases fast with the distance to the source,and corresponds to the self-demodulation of the initial signal.The harmonic decomposition of the response to the Gaussia

54、n pulse with stronger excitation amplitude is shown in Fig.5.A strong demodulation frequency is observed,even more important than the excitation frequency.A new numerical algorithm is developed to study the strongly nonlinear propagation of ultrasonic waves in C.Vanhille,C.Campos-Pozuelo/Ultrasonics

55、42(20041123112811271128 5.00E+6 C. Vanhille, C. Campos-Pozuelo / Ultrasonics 42 (2004 11231128 4.00E+6 simulate strongly nonlinear plane wave propagation of pulsed excitation. Special attention is dedicated to the eect of the initial form of the pulse excitation on the evolution of the signal. 5.5 m

56、icron Pressure (Pa 3.00E+6 Acknowledgements Part of this work has been supported by the CAM 07N/0115/02 and URJC PPR-2003-46 Projects. References 1 C. Campos-Pozuelo, B. Dubus, J.A. Gallego-Jurez, Finite-elea ment analysis of the nonlinear propagation of high intensity acoustic waves, Journal of the Acoustical Society of America 106 (1999 91101. 2 D. Botteldooren, Numerical model for moderately nonlinear sound propagation in three-dimensional structures, Journal of the Acous