AP02_02.vp Introduction When solving a wide range of problems related to nonlin- ear acoustics, we may describe the nonlinear sound waves in fluids by using the Burgers model equation. This equation is named after Johannes Martinus Burgers, who published an equation of this form in his paper concerning turbulent phenomena modelling in 1948 (see [1]). However, this type of equation used in continuum mechanics was first presented in a paper published in a meteorology journal by H. Bateman in 1915 (see [2]). Thanks to the fact that this journal was not read by experts in continuum mechanics, the equation has become known as the Burgers equation. The possibilities of using this equation in nonlinear acoustics were probably discovered by J. Cole in 1949. Subsequently E. Hopf (in 1950) and J. Cole (in 1951) published the transformation inde- pendently in their papers (see [3, 4]). The transformation enables us to find the general analytic solution of the Burgers equation, as a result of which this equation plays a major role in nonlinear acoustics. The Burgers equation can be used for modelling both travelling and standing nonlinear plane waves [5–8]. The equation is the simplest model equation that can describe the second order nonlinear effects connected with the propaga- tion of high-amplitude (finite-amplitude waves) plane waves and, in addition, the dissipative effects in real fluids [9–12]. There are several approximate solutions of the Burgers equa- tion (see [13]). These solutions are always fixed to areas before and after the shock formation. For an area where the shock wave is forming no approximate solution has yet been found. It is therefore necessary to solve the Burgers equation numerically in this area (see [14,15]). Numerical solutions themselves have difficulties with stability and accuracy. A nu- merical solution of the Burgers equation is discussed in this paper. Burgers equation A complete system of equations which describes the prop- agation of waves in thermo-viscous fluids consists of the Navier-Stokes equation, the equation of continuity, the heat- -transfer equation and the thermodynamic state equations. Unfortunately, no general solution of this system has yet been found, and if we were to consider all these equations in their full form, many problems would arise even with a numerical solution (the large number of operations required, instability, etc.). For these reasons, it is advisable to approximate this system by an appropriate model equation in a given approxi- mation, e.g., the quasi-linear homogeneous Burgers equation (see [13,16]). Homogeneous Burgers equation � � � � � � � � � � v x c v v b c v � � � 0 2 0 0 3 2 22 0, (1) here v is acoustic velocity, � is retarded time, � � �t x c0, where t is time, x is the distance from the sound source, c0 is the velocity of sound propagation in the linear approxima- tion, b is the dissipation coefficient, �0 is density of medium, � is nonlinearity coefficient. For the calculations that follow it is advisable to rewrite the equation in the following dimensionless (normed) form: � � � � � � V s V V y G V y � � � 1 0 0 2 2 (2) where y � �� ; s v c x� � �m 0 2 ; V v v � m ; G v c b 0 0 02 � � � � m , G0 represents the Goldberg number, vm is the velocity ampli- tude of the sound source and � is the annular frequency of the sound source. The Burgers equation (1) is a parabolic equation which takes a hyperbolic form in the case of an ideal fluid. The ho- mogeneous Burgers equation (1) describes nonlinear plane travelling waves. Cole-Hopf transformation V G U U y � 2 0 � � . (3) Using transformation (3) we can transform the Burgers equation (2) to the linear diffusion equation [17]: � � � � U s G U y � 1 0 2 2 (4) © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 71 Acta Polytechnica Vol. 42 No. 2/2002 Solution of the Burgers Equation in the Time Domain M. Bednařík, P. Koníček, M. Červenka This paper deals with a theoretical description of the propagation of a finite amplitude acoustic waves. The theory based on the homogeneous Burgers equation of the second order of accuracy is presented here. This equation takes into account both nonlinear effects and dissipation. The method for solving this equation, using the well-known Cole-Hopf transformation , is presented. Two methods for numerical solution of these equations in the time domain are presented. The first is based on the simple Simpson method, which is suitable for smaller Goldberg numbers. The second uses the more advanced saddle point method, and is appropriate for large Goldberg numbers. Keywords: Burgers equation, Cole-Hopf transformation, Saddle point method, nonlinearity. Solution of the Burgers equation The solution of the diffusion equation (4) is given by the Poisson integral, thus the solution of the Burgers equation can be written as the fraction of two integrals with infinite limits: � � � � V y z s G y z s G V y y z G z � � � � � � � � � � � �� � exp , exp 0 2 0 04 2 0 d d � � � �0 2 0 04 2 0 y z s G V y y z z� � � � � � � �� � , d d (5) where V (0, y) is the given boundary condition. In most cases the integrals in equation (5) cannot be solved analytically. Numerical integration It is necessary to integrate Eq. (5) only in intervals � � � �q q i i 1 2, where the value of an integrand is less than a cho- sen parameter � . In general, it is necessary to determine these intervals numerically. When we can express the approximate solution of the Burgers equation in the following form � � � � � � � � V y z s G y z s G V y y z z q q i i i � � � � � � � � � � � exp ,0 2 0 04 2 0 1 2 d d � � � � � � � � 1 0 2 0 0 1 4 2 0 1 2 n z q q i n G y z s G V y y z i i � � � � � � � � � � � exp , d d .(6) In the special cases G0 1� or roughly s < 10, the integra- tion limits analytically can be determined, using the following speculation: The integrand contains the exponential function which for values of z drawing apart from the value of y falls to zero very quickly. Let us choose a small parameter 1 0� � � , � � exp � �� � � � � � � y z G s 2 0 4 . (7) Then we are able to calculate the limits from the equation as q y s G 12 0 2, ln� � � �. (8) The integration itself can be performed by means of the Simpson method, which is based on approximating the inte- grated function in the segments by parts of a parabola � � � � � � � �f z z q q n f z f z f z q q d 1 2 2 1 0 1 2 3 4 2 � � � � � �� [ � � � � � � �� � �� �2 42 1f z f z f zn n n ], (9) where n � 2k, k N� is the number of equidistant intervals and � �z q i q q ni � � �1 2 1 . However, there is a problem with floating-point arithmetic, especially for large values of the number G0. When calculating the integral, we add terms according to (9), and after a short time (especially after inte- grating one of the maxima) the partial is result much larger than the next added terms. Thus these terms are not taken into account, and the result is not accurate. If the accuracy of 18 valid digits is used, we can calculate for values of G0 up to approximately 600. If it were necessary to use this method for larger values of G0, we would have to use some suitable soft- ware that allows calculations of high precision, e.g., the free library GNU MP. Saddle point method If we study the integrand properties, the calculation can be significantly simplified and the problems caused by pro- cessing numbers in high precision can be avoided [17]. When using the saddle point method, the accuracy of the values of the integrals obtained in equation (5) improves with increas- ing Goldberg number G0. For easier orientation we express equation (5) in the fol- lowing form � � � � � � V f z s y G F z s y z G F z s y � � � � �� � � � �� � ; , exp ; , exp ; , 0 0 2 2 d ���� � dz (10) where � �f z s y y z s ; , � � , � � � � � �F z s y y z s V y y z ; , ,� � � 2 0 2 0 d , (11) here s and y represent the parameters. Figs. 1–3 show the courses of the integrands for different Goldberg numbers (G0 � 10, 50, 500) at the various distances (s � 1, 10, 50) for a harmonic boundary condition. It is obvious from the shape of these three-dimensional graphs that the courses of the integrands have significant extrema (max- ima). Consequently we can restrict the integration limits to their close neighbourhood. Further, we can see that when the Goldberg number increases, the values of the integrand extrema also increase and diminish to zero from these ex- 72 Acta Polytechnica Vol. 42 No. 2/2002 Fig. 1: Space evolution of integrand, G0 � 10, s � 1, 10, 50 trema more rapidly (the course of the integrand is narrower). These pictures also show that, on the contrary, in the case of increasing dimensionless distance s near to the extreme, the integrand falls more slowly and the course of integrand spreads (the extrema are not so significant). Consequently we need to take into account two maxima for certain (ap- pointed) values of dimensionless time to do the integration (the maxima superimpose in the direction of the z axis). Let us now analyse the relation (5) for G0 1� and constant s and y. In this case the largest contributions of both integrands come from the extrema of the function F ( z; s, y). From the condition � � � � � � � � � � F z s y z z y z s V y y z ; , ,� � � � � � � � � � 2 0 2 0 0d (12) we obtain � � y z s V z � � �0 0, . (13) Let us mark the solutions of equation (13) z0 and select only those which represent maxima of the integrands � � � � � � � � � � � � � � �F z s y z s V z z 0 2 00 1 0 0 ; , , . (14) Now the contribution of the neighbourhood of this ex- treme can be expressed if we expand the function F ( z; s, y) into the Taylor series at the point z � z0 and we consider only the first three terms (the quadratic order). For the function F ( z; s, y) we shall omit writing parameters s and y; � ��F z and � ���F z will have the meaning of the first and second derivation of this function with respect to variable z. � � � � � �� � � �� � � � � �� � F z F z F z z z F z z z F z F z z z � , � � � � � �� � � � � �� � 0 0 0 0 0 2 0 0 0 2 1 2 1 2 (15) because at the stationary point there is always � �� �F z0 0. Now we substitute this expansion into the integral � �exp �� � �� �� � G F z z0 2 d (16) which we rearrange by means of the well-known Laplace- -Gauss integral e da x 2� �� � � � 2 0z a a , (17) to the form � � � � 4 1 20 0 0 0 G F z G F z �� � � � �� exp . (18) Now let us suppose that there exists only one maximum, z0, that fulfils condition (13), and let us substitute this repre- sentation of integral (16) into the numerator of relation (5). We can see that this term is multiplied by the function f(z) here. The value of this function does not vary significantly in the neighbourhood of the extreme. Then we can treat this function as a constant with value f(z0) (this assumption is satisfied because these extrema are narrow in comparison with the behaviour of the function). On the basis of the foregoing reasoning we can write the numerator of term (5) as � � � � y z s G F z z y z s G F z G � � � � �� � � � �� � �� � exp � � exp 0 0 0 0 0 2 4 1 d � � 2 0F z � � �� (19) and the denominator can be written as � � � � � �exp � exp�� � �� � �� � � � �� �� � G F z z G F z G F z0 0 0 0 0 2 4 1 2 d . (20) Then we have V y z s �� � 0 , (21) where z0 is given by expressions (13) and (14). © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 73 Acta Polytechnica Vol. 42 No. 2/2002 Fig. 2: Space evolution of integrand, G0 � 50, s � 1, 10, 50 Fig. 3: Space evolution of integrand, G0 � 500, s � 1, 10, 50 However, relation (21) can give wrong solutions in certain cases. The explanation is that there are more stationary points satisfying equation (13) and condition (14), usually two. Now let us suppose two maxima z1 and z2, when z z1 2� . Their contributions are determined by relations (19) and (20), and thus: � � � � � � V y z s G F z G F z G F z i� � �� � � � �� �� � � i i i i 4 1 2 4 1 0 0 1 2 0 exp ex � �p �� � �� � � G F z i 0 1 2 i 2 . (22) This expression can be extended to an arbitrary number of maxima zi. When � � � �F z F z1 2� , the contribution of one extreme becomes very high due to large G0 – for example in the case G0 � 20000 and � � � �F z F z1 2 1� � this rate is equal to e 10000 1: . For this reason, we can use the approximate solution in some circumstances. � � � � � � � � V y z s F z F z V y z s F z F z � , � , . � � � � � � 1 1 2 2 2 1 when and when (23) Examples of the solution of the Burgers equation Figure 4 a) shows waveforms with the harmonic boundary condition when G0 � 1500. The waveforms are calculated us- ing the saddle point method at dimensionless distances which were chosen in order to show the development of a nonlinear wave. Figure 4 b) shows waveforms for the identical boundary condition at the same distances as in Figure 4 a), but for G0 � 50. Direct integration was used in this case. We can see that the dissipation effect is stronger in this case, and there- fore the width of the shock front is more significant. Figure 5 compares the time domain solution (saddle point method) with the frequency domain solution of the Burgers equation (see [18]). Conclusion This paper analyses a method for solving the Burgers equation in the time domain by means of the Cole-Hopf transformation. However, this method poses problems con- nected with the numerical solution of the Poisson integral for the case G0 � 1. The first way to solve these problems is based on the application of an appropriate software instrument (e.g., the 74 Acta Polytechnica Vol. 42 No. 2/2002 a) b) Fig. 4: Waveforms at distances s � 0, 1, 4, 10, 20 for a) G0 � 1500 b) G0 � 50 Fig. 5: Comparison between the time domain solution and the frequency domain solution of the Burgers equation, G0 � 1000, s � 1.1 and s � 3, M represents number of harmonics free library GNU MP), which enables calculation with accu- racy to an arbitrary number of valid digits. The next way involves applying the saddle-point method, which becomes more accurate the larger the Goldberg numbers are. The solution in the time domain enables us to find the proper parameters of the numerical frequency correction that is often used for solving the Burgers equation in the frequency domain. When we compare the results of the solutions of the Burgers equation obtained in the time domain and in the frequency domain, it is evident that they are almost identical in the case when the number of used harmonics in the frequency domain exceeds 250, but the calculation is notice- ably longer. Acknowledgment The research was supported by CTU research program 16:CEZ:J04/98:212300016 and by GACR grant 202/01/1372. References [1] Burgers, J. M.: A mathematical model illustrating the theory of turbulence. Advances in Applied Mechanics, 1948, Vol.1, p. 171–199. [2] Bateman, H.: Some recent researches on the motions of fluids. Monthly Weather Review, 1915, Vol. 43, p. 163–170. [3] Hopf, E.: The partial differential equation u uu ut x xx� � . Comm. Pure and Appl. Math., 1950, Vol. 3, p. 201–230. [4] Cole, J. D.: On a quasi-linear parabolic equation occurring in aerodynamics. Q. Appl. Math., 1951, Vol. 9, p. 225–236. [5] Koníček, P., Bednařík, M., Červenka, M.: Finite-ampli- tude acoustic waves in a liquid-filled tube. Hydroacoustics (Editors: E. Kozaczka, G. Grelowska), Gdynia: Naval Academy, 2000, p. 85–88. [6] Bednařík, M.: Description of plane finite-amplitude traveling waves in air-filled tubes. Acoustica, 1996, Vol. 82, p. 196. [7] Bednařík, M.: Description of high-intensity sound fields in fluids. Proceedings of the 31st Conference on Acoustics, Prague, 1994, p. 134–138. [8] Ochmann, M.: Nonlinear resonant oscilllations in closed tubes – An application of the averaging method. J. Acoust. Soc. Am., 1985, Vol. 77, p. 61–66. [9] Bednařík, M.: Description of nonlinear waves in gas-filled tubes by the Burgers and the KZK equations with a fractional derivative. Proceedings of the 16th ICA, ICA, Seattle, 1998, p. 2307–2308. [10] Blackstock, D. T.: Connection between the Fay and Fubuini Solution for Plane Sound Waves of Finite Amplitude. J. Acoust. Soc. Am., 1966, Vol. 39, p. 1019–1026. [11] Blackstock, D. T.: Thermoviscous Attenuation of Plane, Peri- odic, Finite-Amplitude Sound Waves. J. Acoust. Soc. Am., 1964, Vol. 36, p. 534–542. [12] Webster, D. A., Blackstock, D. T.: Finite-amplitude satura- tion of plane sound waves in air. J. Acoust. Soc. Am., 1977, Vol. 62, p. 518–523. [13] Hamilton, M. F., Blackstock, D. T.: Nonlinear Acoustic. USA: Academic Press, 1998, p. 117–139. [14] Bednařík, M.: Nonlinear propagation of waves in hard- -walled ducts. Proceedings of the CTU Workshop, Prague: CTU, 1995, p. 91–92. [15] Mitome, H.: An exact solution for finite-amplitude plane sound waves in a dissipative fluid. J. Acoust. Soc. Am., 1989, Vol. 86, No. 6, p. 2334–2338. [16] Makarov, S., Ochmann, M.: Nonlinear and thermoviscous phenomena in acoustics 2. Acustica, 1997, Vol. 83, p. 197–222. [17] Whitham, G. B.: Linear and Nonlinear Waves. New York: John Willey, 1974, p. 152–157. [18] Trivett, D. H., Van Buren, A. L.: Propagation of plane, cy- lindrical, and spherical finite amplitude waves. J. Acoust. Soc. Am., 1981, Vol. 69, No. 4, p. 943–949. Dr. Ing. Michal Bednařík e-mail: bednarik@fel.cvut.cz Dr. Mgr. Petr Koníček e-mail: konicek@fel.cvut.cz Ing. Milan Červenka e-mail: cervenm3@fel.cvut.cz Department of Physics Czech Technical University in Prague Faculty of Electrical Engineering Technická 2 166 27 Praha 6, Czech Republic © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 75 Acta Polytechnica Vol. 42 No. 2/2002 Table of Contents Concept Design and Reliability 3 G. Cooper, G. Thompson The Impact of Odors 13 M. V. Jokl Optimization of the Odor Microclimate 23 M. V. Jokl A Low-complexity Wavelet Based Algorithm for Inter-frame Image Prediction 30 A. K. Haghi Sampling by Fluidics and Microfluidics 41 V. Tesaø Mixing Suspensions in Slender Tanks 50 F. Rieger, E. Rzyski Formation of Haloforms during Chlorination of Natural Waters 56 A. Grünwald, B. Š•astný, K. Slavíèková, M. Slavíèek Reducing Ultrasonic Signal Noise by Algorithms based on Wavelet Thresholding 60 M. Kreidl, P. Houfek Optimisation and Just-in-Time 66 V. Beran Solution of the Burgers Equation in the Time Domain 71 M. Bednaøík, P. Koníèek, M. Èervenka