流體中非線性超聲波數(shù)值仿真有限差分

1、A numerical formulation for nonlinear ultrasonic wavespropagation in uidsC.Vanhillea,*,C.Campos-PozuelobaESCET,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

2、2004;accepted 10February 2004Available online 24February 2004AbstractKeywords:Nonlinear propagation;Numerical modelling;High-power ultrasound1.IntroductionIn this paper we propose a new numerical algorithm for studying the nonlinear propagation of transient and periodical signals in uid media.The wo

3、rk 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 research (test and ther-apy,etc.Several works dealing with analytical and numerical studies of propagat

4、ion 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 wave-amplitude 1,2,and two-dimensional models have reached the simulation of nite amplitude wave prop

5、agation 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 approximations are based on Eulerian coordi-nates and only second order terms of the Mach number ar

6、e 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 by one the order of the dierential equation for wave motion 4,5.The purpose of the present work is t

7、o 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 the acoustic Mach number.The formulation is written in the time domain.Natural spatial and time coord

8、inates 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 the waveform evolution for any original signal is possible:periodic excitation,Gaussian,rectan-gular

9、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 development of an implicit nite-dierence scheme and by using a fast*Corresponding author.Tel.:+34-91-664

10、-74-82;fax:+34-91-488-73-38. Ultrasonics 42(2004 11231128linear solver.The equations 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.Govern

11、ing equationsThe full nonlinear one-dimensional wave equation can be derived by 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 a

12、mbient pressure,c is the specic heat ratio,q 0is the initial state density,u is 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 M

13、ach number value or about the attenuation parameter.The only limitations on pressure 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 isentro

14、pic property of the uid can not be questioned.The acoustic pressure waveform and 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.

15、e.,particle displacement and velocity null at t ¼0.The following auxiliary 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

16、<:ð3Þwhere c 0is the small-amplitudes value of the sound speed,L 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-gr

17、essive character of the wave is imposed by proposing the nonreecting boundary 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 t

18、he signal and x is its angular frequency.The X T space is discretized by means 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 discretizat

19、ion leads to a dierence equation for each space grid point.The nonreecting boundary 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 a

20、s the value of the displacement u at the point m ;n of the grid corresponding 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

21、ÀAu m þ1;n¼ÀAu m À1;n À2þð2A À1Þ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 À

22、;1ÞÀc À1m ¼2;.;M À1ð4b ÞDu m À1;n þ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 Þ

23、8;e Àx 2b ððn À1Þsx Àt 0Þ2cos ððn À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 g

24、rid point values of u .The scheme is implicit.The following economic and fast 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 symmetr

25、ic,and the method can be improved.The rst and last equations need a specic treatment.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

26、 U þ¼ðb b .b b ÞT and a new n -vector L Àis created:L À¼ð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

27、 ÞU þðk þ1Þ¼Àa L Àðk þ1Þþ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

28、are used at each time step (A is always the sameby solving two 2-diagonal systems.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 

29、40;i ÞÀL Àði Þz ði À1Þi ¼2;.;nz ð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

30、before the fol-lowing time step.In addition to the n -vector U þ,scalars 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 co

31、nsistent scheme to be conditionally stable,and an optimised relation between h 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

32、b is the frequency of the error and n is the amplication factor.This study leads 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 unconditionall

33、y stable and does not generate any additional convergence problem.In addition 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.

34、First,simulations of strongly nonlinear waves prop-agating in a quite large space 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 displac

35、ements in the unperturbed zone (mechanical perturbations reach the i th wavelength 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 periodOn the other h

36、and,if a simulation has been carried out without a sucient number of periods,this can be started 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 th

37、e previous study,completed by zero values for the (NP þ2th wavelength.Once the displacement 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 o

38、f numerical experiments are presented.All the results are obtained from con-vergent simulations with a suitable number of spatial points. secondC.Vanhille,C.Campos-Pozuelo /Ultrasonics 42(2004112311281125harmonics.The agreement validates the numerical model,as well as the proposed boundary condition

39、 at x¼L. 1126 C.Vanhille,C.Campos-Pozuelo/Ultrasonics42(2004112311285.ConclusionsA new numerical algorithm is developed to study the strongly nonlinear propagation of ultrasonic waves in C.Vanhille,C.Campos-Pozuelo/Ultrasonics42(20041123112811271128 5.00E+6 C. Vanhille, C. Campos-Pozuelo / Ultr

40、asonics 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 micron Pressure (Pa 3.00E+6 Acknowledgements Part of this work has bee

41、n 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, Nu

42、merical model for moderately nonlinear sound propagation in three-dimensional structures, Journal of the Acoustical Society of America 100 (1996 13571367. 3 V.W. Sparrow, R. Raspet, A numerical method for general nite amplitude wave propagation in two dimensions and its application to spark pulses, Journal of the Acoustical Society of America 90 (1991 26832691. 4 M.H. Hamilton, D.T. Blackstock (Eds., Nonlinear Acoustics, Academic Press, New York, 1998. 5 D.T. Blackstock, Propagation of plane sound waves of nite amplitude in nondissipative uids, Journal of the A