Analytical and numerical modeling of front propagation and interaction of fronts in nonlinear thermoviscous fluids including dissipation
Anders R. Rasmussen, Mads P. Sørensen, Yuri B. Gaididei, Peter L. Christiansen
AAnalytical and numerical modeling of front propagation and interaction of fronts innonlinear thermoviscous fluids including dissipation
Anders R. Rasmussen ∗ and Mads P. Sørensen Department of Mathematics, Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark
Yuri B. Gaididei
Bogolyubov Institute for Theoretical Physics, 03680 Kiev, Ukraine
Peter L. Christiansen
Department of Informatics and Mathematical Modelling and Department of Physics,Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark (Dated: November 7, 2018)A wave equation, that governs finite amplitude acoustic disturbances in a thermoviscous Newto-nian fluid, and includes nonlinear terms up to second order, is proposed. In contrast to the modelknown as the Kuznetsov equation, the proposed nonlinear wave equation preserves the Hamilto-nian structure of the fundamental fluid dynamical equations in the non-dissipative limit. An exacttraveling front solution is obtained from a generalized traveling wave assumption. This solution is,in an overall sense, equivalent to the Taylor shock solution of the Burgers equation. However, incontrast to the Burgers equation, the model equation considered here is capable to describe wavespropagating in opposite directions. Owing to the Hamiltonian structure of the proposed modelequation, the front solution is in agreement with the classical Rankine-Hugoniot relations. The ex-act front solution propagates at supersonic speed with respect to the fluid ahead of it, and subsonicspeed with respect to the fluid behind it, similarly to the fluid dynamical shock. Linear stabilityanalysis reveals that the front is stable when the acoustic pressure belongs to a critical interval,and is otherwise unstable. These results are verified numerically. Studies of head-on colliding frontsdemonstrate that the front propagation speed changes upon collision.
PACS numbers: 43.25.Cb, 43.25.Jh, 43.25.TsKeywords: thermoviscous fluids, traveling fronts, Rankine-Hugoniot relations, shocks
I. INTRODUCTION
The “classical” equation of nonlinear acoustics [1], theso-called Kuznetsov equation [2], governs finite amplitudeacoustic disturbances in a Newtonian, homogeneous, vis-cous, and heat conducting fluid. The model equation andits paraxial approximation, the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation [2, 3], are occasionally en-countered within studies related to nonlinear wave prop-agation. See e.g. the recent works by Jordan [4] whopresented the derivation and analysis of an exact trav-eling wave solution to the one-dimensional Kuznetsovequation, and by Jing and Cleveland [5] who described athree-dimensional numerical code that solves a general-ization of the KZK equation, and the references cited inthe introductory sections of those papers. Other recentworks based on the Kuznetsov equation include: anal-ysis of energy effects accompanying a strong sound dis-turbance [6], studies of generation of higher harmonicsand dissipation based on a 3D finite element formulation[7], and studies of nonlinear wave motion in cylindricalcoordinates [8]. The derivations of the Kuznetsov equa-tion [2, 9, 10] and related model equations [11, 12, 13] ∗ Electronic address: anders r [email protected] are based on the complete system of the equations offluid dynamics. It has been demonstrated that this sys-tem of equations is of Hamiltonian structure in the non-dissipative limit [14]. However, in the non-dissipativelimit, the Kuznetsov equation does not retain the Hamil-tonian structure.In this paper we propose a nonlinear wave equation,which, in the non-dissipative limit, preserves the Hamil-tonian structure of the fundamental equations. Further-more, we present the derivation and analysis of an ex-act traveling front solution, which applies equally well tothe proposed nonlinear wave equation and the Kuznetsovequation. The derivation of the exact solution is basedon a generalized traveling wave assumption, which leadsto a wider class of exact solutions compared to the onereported by Jordan [4, 15]. Furthermore, the introduc-tion of the generalized assumption is necessary in order tointerpret the results of numerical simulations of head-oncolliding fronts presented in this paper. In order to re-late our results to the classical literature, we demonstratethat the exact front solution retains a number of prop-erties of the fluid dynamical shock. The paper is struc-tured as follows: The proposed equation and its Hamilto-nian structure are discussed in Section II. Section III con-tains the derivation of our exact traveling front solutionand analysis of its stability properties. In Section IV wedemonstrate that the front is related the classical shock.Section V presents numerical investigations of the front, a r X i v : . [ phy s i c s . f l u - dyn ] J un while Section VI contains our conclusions. II. NONLINEAR WAVE EQUATIONS
Equations governing finite amplitude acoustic distur-bances in a Newtonian, homogeneous, viscous and heatconducting fluid may be derived from the following fourequations of fluid dynamics: the equation of motion ρ (cid:18) ∂ u ∂t + ( u · ∇ ) u (cid:19) = − ∇ p + η ∆ u + (cid:16) η ζ (cid:17) ∇ ( ∇ · u ) , (1) the equation of continuity ∂ρ∂t + ∇ · ( ρ u ) = 0 , (2) the heat transfer equation ρT (cid:18) ∂s∂t + ( u · ∇ ) s (cid:19) = η (cid:18) ∂u i ∂x j + ∂u j ∂x i − δ ij ∂u k ∂x k (cid:19) + ζ ( ∇ · u ) + κ ∆ T, (3)and the equation of state p = p ( ρ, s ) . (4)Here x = ( x, y, z ) are the spatial (Cartesian) coordinatesand t denotes time. u = ( u, q, w ) is the fluid particlevelocity, ρ is the density of the medium, p , s , and T are the thermodynamic variables pressure, entropy andtemperature, respectively. η and ζ are the coefficients ofshear and bulk viscosity, and κ is the heat conductivitycoefficient. ∆ is the Laplace operator.To obtain a nonlinear wave equation all dependentvariables except one are eliminated from the system (1)–(4), resulting in a nonlinear wave equation for that singlevariable [2, 9, 10, 11, 12, 13]. The deviations of ρ , p , s , and T from their equilibrium values, ρ , p , s , and T are assumed to be small, as well as the fluid parti-cle velocity, | u | . The heat conductivity coefficient κ andthe viscosities η and ζ are also treated as small quanti-ties. In order to obtain a second order approximation, allequations are written retaining terms up to second orderin the small quantities. It is assumed that the flow isrotation free, ∇ × u = 0, thus u ≡ − ∇ ψ, (5)where ψ is the velocity potential. Furthermore, it hasbecome customary to use the following approximationfor the equation of state [16] p − p = c ( ρ − ρ ) + c ρ B/A ρ − ρ ) + (cid:18) ∂p∂s (cid:19) ρ,s = s ( s − s ) , (6) TABLE I: Values of c , B/A , and b for three different sub-stances. The values for b are rough estimates obtained fromEq. (8) neglecting the influence of bulk viscosity and thermallosses.Substance c (m s − ) B/A b (m s − )Water 1483 (20 (cid:137) ) a . (cid:137) ) b . × − (20 (cid:137) ) c Air 343 (20 (cid:137) ) a . (cid:137) ) b,d × − (27 (cid:137) ) c Soft tissue 1540 e . (cid:137) ) b,f N/A a Ref. 19 b Ref. 17 c Values for ρ and η are obtained from Ref. 19. d Diatomic gas e Ref. 20 f Human breast fat where
B/A is the fluid nonlinearity parameter [17] and c ≡ ( ∂p/∂ρ ) s,ρ = ρ is the small-signal sound speed.Then, from Eqs. (1)–(6) we obtain the following non-linear wave equation ∂ ψ∂t − c ∆ ψ = ∂ψ∂t ∆ ψ + ∂∂t (cid:32) b ∆ ψ + ( ∇ ψ ) + B/A − c (cid:18) ∂ψ∂t (cid:19) (cid:33) , (7)where b is the diffusivity of sound [18] b ≡ ρ (cid:26) η + ζ + κ (cid:18) C V − C p (cid:19)(cid:27) , (8)and C V and C p denote the heat capacities at constantvolume and pressure, respectively. Typical values of thephysical parameters c , B/A , and b are given in Table I.In the first order approximation, Eq. (7) reduces to ∂ ψ∂t = c ∆ ψ. (9)Introducing Eq. (9) in the first term on the right handside of Eq. (7), the Kuznetsov equation [2] ∂ ψ∂t − c ∆ ψ = ∂∂t (cid:32) b ∆ ψ + ( ∇ ψ ) + B/A c (cid:18) ∂ψ∂t (cid:19) (cid:33) , (10)is obtained.In absence of dissipation, i.e. η = ζ = 0, Eqs. (1) and(2) possess Hamiltonian structure [14]. This property is,however, not retained in Eq. (10) with b = 0, i.e. the non-dissipative limit of the Kuznetsov equation is not Hamil-tonian. In contrast, Eq. (7) does retain the Hamiltonianstructure in the non-dissipative limit. Accordingly, theequation may be derived from the Lagrangian density L = ( ψ t ) − c ( ∇ ψ ) − B/A − c ( ψ t ) − ψ t ( ∇ ψ ) , (11)using the Euler-Lagrange equation . From the Legendretransformation [22] we obtain the corresponding Hamil-tonian density as H = c ( ∇ ψ ) ψ t ) − B/A − c ( ψ t ) , (12)which may be integrated to yield the total Hamiltonian H = (cid:90) + ∞−∞ (cid:90) + ∞−∞ (cid:90) + ∞−∞ H dx dy dz. (13)Taking the time derivative of H in Eq. (13) with H re-placed by Eq. (12), and using Eq. (7), one can obtain asimple expression for dH/dt . Doing this in one spatialdimension we obtain after some calculations dHdt = (cid:104) c ψ t ψ x + ( ψ t ) ψ x (cid:105) + ∞−∞ − b (cid:90) + ∞−∞ ( ψ xt ) dx. (14)In Eq. (14), which is sometimes called the energy balanceequation, we observe that the first terms on the righthand side correspond to energy in- and output at the twoboundaries, and that the last term accounts for energydissipation inside the system.In the remaining portion of this paper we shall limitthe analysis to one-dimensional plane fields, in which casethe proposed model equation (7) reduces to ψ tt − c ψ xx = ψ t ψ xx + ∂∂t (cid:18) bψ xx + ( ψ x ) + B/A − c ( ψ t ) (cid:19) , (15)where subscripts indicate partial differentiation.Finally, for later reference we give the second-orderexpressions for the acoustic density, ρ − ρ , and acousticpressure, p − p , in terms of the velocity potential, ψ .From the equations of motion (1) and state (6), subjectto the basic assumptions of the derivation of the twomodel equations (7) and (10), we obtain ρ − ρ = ρ c (cid:32) ψ t − ( ψ x ) − B/A − c ( ψ t ) − bψ xx (cid:33) , (16) Letting η = ζ = 0, u = − ∇ ψ , and p = ρ γ /γ in Eqs. (1–2), and(4), respectively, one can derive the Lagrangian density L PEE = c γ „ γ − c „ ψ t − ( ∇ ψ ) «« γγ − , corresponding to the potential Euler equation (PEE) given inRef. 21. Expanding L PEE to third order and letting γ = B/A +1we obtain Eq. (11). and p − p = ρ (cid:32) ψ t − ( ψ x ) c ( ψ t ) (cid:33) − (cid:18) η + ζ (cid:19) ψ xx , (17)respectively. It should be noted that Eqs. (16) and (17)are derived from the fundamental equations, thus, theexpressions are not specific to any of the two model equa-tions (7) and (10). III. EXACT TRAVELING FRONT SOLUTION
Recently, a standard traveling wave approach wasapplied to the one-dimensional approximation of theKuznetsov equation (10) to reveal an exact traveling wavesolution [4, 15]. In this section we extend the standardapproach by introducing a generalized traveling wave as-sumption and analyze the stability properties of the so-lution.
A. Generalized traveling wave analysis
We introduce the following generalized traveling waveassumption ψ ( x, t ) = Ψ( x − vt ) − λx + σt ≡ Ψ( ξ ) − λx + σt, (18)where λ and σ are arbitrary constants, v denotes thewave propagation velocity, and ξ ≡ x − vt is a wavevariable. The inclusion of − λx + σt in Eq. (18) leadsto a wider class of exact solutions, compared to the oneobtained from the assumption ψ = Ψ( x − vt ), which isthe standard one. Furthermore, the introduction of thegeneralized assumption is necessary in order to interpretthe results of numerical simulations of head-on collidingfronts presented in Section V B. Inserting Eq. (18) intothe nonlinear wave equation (15) we obtain the ordinarydifferential equation (cid:0) v − c (cid:1) Ψ (cid:48)(cid:48) = ( − v Ψ (cid:48) + σ ) Ψ (cid:48)(cid:48) − v ddξ (cid:26) b Ψ (cid:48)(cid:48) + (Ψ (cid:48) − λ ) + B/A − c ( − v Ψ (cid:48) + σ ) (cid:27) , (19)where prime denotes ordinary differentiation with respectto ξ . Integrating once and introducing Φ ≡ − Ψ (cid:48) , Eq. (19)reduces to C = vb Φ (cid:48) − (cid:18)
32 +
B/A − c v (cid:19) v Φ + (cid:26)(cid:18) − B/A − c σ (cid:19) v − λv − c − σ (cid:27) Φ , (20)where C is a constant of integration. Requiring that thesolution satisfy Φ (cid:48) → ξ → ±∞ , and eitherΦ → (cid:26) θ, ξ → + ∞ , ξ → −∞ or Φ → (cid:26) , ξ → + ∞ θ, ξ → −∞ , (21)where θ is an arbitrary constant, lead us to C = 0 and B/A − c θv − (cid:18) − B/A − c σ (cid:19) v + (cid:18) θ + 2 λ (cid:19) v + c + σ = 0 . (22)In order to obtain our traveling wave solution, we sep-arate the variables in Eq. (20) subject to C = 0, then,using Eq. (22), we find the solution to be the travelingfront Φ = θ (cid:26) − tanh (cid:18) ξ − x ) l (cid:19)(cid:27) , (23) l ≡ b (cid:18) B/A − c v + 32 (cid:19) θ , (24)where x is an integration constant, | l | is the front thick-ness, and 0 < Φ < θ . Finally, using Φ = − Ψ (cid:48) and in-serting Eq. (23) into Eq. (18) we obtain (apart from anarbitrary constant of integration) ψ ( x, t ) = − θ (cid:26) ξ − l (cid:18) cosh 2( ξ − x ) l (cid:19)(cid:27) − λx + σt, (25)which is the exact solution for the velocity potential.Traveling tanh solutions, such as the front solu-tion (23), are often called Taylor shocks. The existenceof an exact solution of this type to the classical Burgersequation is a well known result [23]. However, the Burg-ers equation is restricted to wave propagation either tothe left or to the right. The model equation consideredin this paper does not suffer from this limitation, as shallbe illustrated in Section V B.Regarding the exact solution derived above, the phys-ical properties of the flow associated with the travelingfront are obtained from the partial derivatives of Eq. (25),which are given by − ψ x = Φ + λ and ψ t = v Φ + σ. (26)According to Eq. (5) the fluid particle velocity is ob-tained as u = − ψ x , and the first order approximation ofEq. (17) yields the acoustic pressure as p − p ≈ ρ ψ t .The boundary conditions of the front are obtained fromEqs. (23) and (26) as − ψ x → (cid:26) θ + λ, ξ → ∓∞ λ, ξ → ±∞ , (27a) ψ t → (cid:26) vθ + σ, ξ → ∓∞ σ, ξ → ±∞ , (27b) where upper (lower) signs apply for l > l < v , θ , λ , and σ , that was introducedin the derivation of the exact solution, determine the fourboundary conditions of the front. From these boundaryconditions we find that θ and vθ correspond to the heightsof the jump across the front measured in − ψ x and ψ t ,respectively, see Fig. 1. At this point it is appropriate toemphasize that, in order for the exact solution to exist,the four parameters v , θ , λ , and σ must satisfy the cubicequation (22). Furthermore, the allowable values of thewave propagation velocity correspond to the real roots ofthis equation. A noticeably property of Eq. (22), whichwill prove useful later on, is that the equation is invariantunder the transformation v → v, θ → − θ, λ → θ + λ, σ → vθ + σ. (28)Also the boundary conditions (27) are invariant, since,according to Eq. (24), the above transformation leads to l → − l . σνθ+σ ξ ψ t l λθ+λ ξ − ψ x l FIG. 1: The exact solution of Eq. (15) represents a travelingfront. In order for the solution to exist, the wave propagationvelocity, v , the front height, θ , and the two constants, λ and σ , must satisfy Eq. (22). The plot shows a front with l > In order to investigate the relationship between thefront height, θ , and the front propagation velocity, v , wesolve Eq. (22) with respect to θ to obtain θ = (cid:18) − B/A − c σ (cid:19) v − λv − c − σv (cid:18)
32 +
B/A − c v (cid:19) . (29)For B/A < θ ( v ) has singularities at v = v s ≡ ± (cid:18) c − B/A (cid:19) / , (30)and for B/A > at ( v, θ ) =( v max , θ max ), where v max is obtained as v max = c (cid:32) B/A + (cid:112) B/A ) + 12( B/A − B/A − (cid:33) / , (31) The critical point ( v, θ ) = ( v max , θ max ) was identified by Jordan[4] as the solution bifurcation point. when λ = 0 and σ = 0. These two characteristic proper-ties of the curve are illustrated in Fig. 2. B/A ≤ θ / c v/c B/A ≤
10 10.8
B/A ≤
10 1 5 10 150123 θ / c v/c B/A ≥
11 1.11.55
B/A ≥
11 1.11.55
B/A ≥
11 1.11.55
B/A ≥ FIG. 2: The plots show the relationship between the frontheight, θ , and the front propagation velocity, v , given byEq. (29) with λ = 0, σ = 0, and B/A = { , . , , . , . , } (see labels on the plots). The dashed lines indicate the singu-larity at v = v s , which is defined in Eq. (30), and crosses in-dicate the maximum ( θ, v ) = ( θ max , v max ) defined in Eq. (31). Finally, it should be emphasized that the general-ized traveling wave analysis conducted above also appliesto the one-dimensional approximation of the Kuznetsovequation (10). In this case Eq. (22) is replaced by B/A c θv − (cid:18) − B/Ac σ (cid:19) v + ( θ + 2 λ ) v + c = 0 , (32)and Eq. (24) by l = 4 b (cid:18) B/A c v + 1 (cid:19) θ . (33)Apart from these changes, a generalized traveling waveanalysis of the Kuznetsov equation is basically identicalto that of Eq. (15).The Hamiltonian structure, however, is unique tothe proposed nonlinear wave equation (7) and its one-dimensional approximation Eq. (15). In order to estab-lish a relationship between the exact solution, derived inthis section, and the Hamiltonian structure of the gov-erning equation, we insert Eq. (25) into Eq. (14) and theone-dimensional approximations of Eqs. (12) and (13).Then, after some calculations, we find that Eq. (14) re-duces to the cubic equation (22). Hence, the exact trav-eling front solution of the proposed Hamiltonian modelequation (15) satisfies the energy balance equation (14). B. Linear stability analysis
In order to gain insight into the stability propertiesof the traveling front solution we initially consider the Eliminating λ and σ from Eq. (32) makes the equation equivalentto the previously reported result [4] constant solution − ψ x = K and ψ t = L, (34)which satisfies the nonlinear wave equation (15). The twoconstants K and L are arbitrary. In order to investigatethe linear stability properties of this solution, we addsmall perturbation terms to the constant values as − ψ x = K − εχ x and ψ t = L + εχ t , (35)where χ = χ ( x, t ) and ε (cid:28)
1. Then, inserting Eqs. (35)into Eq. (15) and keeping terms up to first order in ε weobtain the following linear perturbation equation (cid:18) − B/A − c L (cid:19) χ tt − (cid:0) c + L (cid:1) χ xx = bχ xxt − Kχ xt . (36)Inserting the single Fourier mode χ ( x, t ) = D e i( kx − ωt ) , (37)where D is the amplitude, k is the wave number, and ω is the angular frequency, into Eq. (36), we obtain thefollowing dispersion relation (cid:18) B/A − c L − (cid:19) ω + (cid:0) Kk − ibk (cid:1) ω + (cid:0) c + L (cid:1) k = 0 . (38)The constant solution (34) is asymptotically stable onlyif all solutions of Eq. (36) approach zero as t → ∞ .This is the case when the imaginary part of both rootsin Eq. (38), ω and ω , are negative. It can be shownthat for B/A > − c < L < c B/A − . (39)When B/A < − c < L < ∞ . (40)Hence, the stability properties of the constant solution(34) are determined exclusively by L , i.e. the constantvalue of ψ t . Recall that the acoustic pressure is pro-portional to ψ t , thus, the level of the acoustic pressuredetermines the stability properties of the solution.In order for the front solution to be stable, it is a neces-sary condition that both left and right asymptotic valuesof ψ t , given by Eq. (27b), belong to the interval (39)when B/A >
1, and the interval (40) when
B/A <
IV. FRONT-SHOCK RELATIONSHIP
Within fluid dynamics, a shock denotes a sharp changeof the physical quantities. A shock propagates at super-sonic speed with respect to the fluid ahead of it, whileit remains subsonic with respect to the fluid behind it.The physical quantities of the flow on each side of theshock are connected by the Rankine-Hugoniot relations,which are conservation equations for mass, momentumand energy. In the following we shall demonstrate thatthe front solution of the proposed Hamiltonian modelequation (15) retains these properties.
A. The Rankine-Hugoniot relations
Using square brackets to denote the change in value ofany quantity across a shock, e.g.[ ρ ] = ρ a − ρ b , (41)where b denotes the value behind the shock and a de-notes the value ahead of it, the Rankine-Hugoniot rela-tions may be written as [24]mass : [ ρ ( u − v )] = 0 , (42)momentum : (cid:104) p + ρ ( u − v ) (cid:105) = 0 , (43)energy : (cid:104) h + ( u − v ) / (cid:105) = 0 , (44)where v is the shock propagation velocity and h is theenthalpy.We now replace u , ρ , p , and h with expressions in termsof ψ x and ψ t , and write all equations retaining terms upto second order. Upon setting u = − ψ x and substitutingEqs. (16) and (17) into Eqs. (42) and (43) we thus obtain (cid:34)(cid:32) ( ψ x ) B/A − c ( ψ t ) (cid:33) v − ψ t ( ψ x + v ) − c ψ x (cid:35) = 0 , (45)and (cid:34) B/A − c ( ψ t ) v − ψ t v − ψ t ψ x v − ( ψ t ) − c ψ x v − c ψ t (cid:35) = 0 , (46)respectively. The dissipative terms involving κ , ζ , and η do not appear in Eqs. (45) and (46), since ψ xx → ψ x and ψ t acrossthe front are obtained from the boundary conditions (27).Assuming that l > v >
0, and using the notationintroduced in Eq. (41) we may write[ ψ x ] = θ, [ ψ t ] = − vθ. (47a) Furthermore, changes in products of ψ x and ψ t are (cid:104) ( ψ x ) (cid:105) = − θ − θλ, (47b) (cid:104) ( ψ t ) (cid:105) = − v θ − vθσ, (47c) (cid:104) ψ x ψ t (cid:105) = vθ + vθλ + θσ. (47d)Inserting Eqs. (47) into Eqs. (45) and (46) both conser-vation equations reduce to the cubic equation (22). Thisstriking result leads to the conclusion, that Eq. (22) im-plies conservation of mass and momentum. At this pointit should be noted that the generalized traveling waveanalysis of the Kuznetsov equation (10) leads to the cu-bic equation (32), which is not in agreement with theconservation equations for mass and momentum.In order to handle the enthalpy in the condition for en-ergy conservation (44) we shall make use of the followingfundamental thermodynamic relationship [24] ∇ h = ∇ pρ . (48)Using the equation of motion (1), subject to the basicassumptions of the derivation of the model equations inSection II, we obtain from Eq. (44)[ ψ t + vψ x ] = 0 . (49)Alternatively, Eq. (49) follows directly from the general-ized traveling wave assumption (18). Hence, the travelingwave assumption implies energy conservation in the flow. B. Sub-/supersonic speeds of propagation
In order to determine whether the traveling front so-lution, derived in Section III A, propagates at sub- or su-personic speed with respect to the fluid ahead of it andthe fluid behind it, we need to introduce the speed ofsound in these regions of the fluid. Without loss of gen-erality, we may consider only fronts propagating in thepositive direction, v >
0, since Eq. (15) is invariant underthe transformation x → − x . Furthermore, we shall limitthe analysis to stable fronts, i.e. ψ t must belong to theinterval (39) when B/A >
1, and the interval (40) when
B/A <
1. Then, letting θ → v yields the small signal propagation velocity, whichis equivalent to the speed of sound, c . Introducing λ = K and σ = L we obtain v = c ( K, L ) ≡ K + (cid:115) K + ( L + c ) (cid:18) − B/A − c L (cid:19) − B/A − c L , (50)where K and L denote the constant levels of − ψ x and ψ t ,respectively, at which the speed of sound (50) is evalu-ated. Inserting the boundary conditions of the front intoEq. (50), i.e. substituting Eq. (27a) for K and Eq. (27b)for L , we obtain the speed of sound ahead of, c a , andbehind, c b , the front c ab = c ( λ, σ ) , (51) c ba = c ( θ + λ, v θ + σ ) , (52)where upper (lower) subscripts apply for l > l < v ,to c a and c b we make the following observations. Insert-ing Eq. (29) into Eq. (24) yields l = 4 bv (cid:18) − B/A − c σ (cid:19) v − λv − c − σ . (53)The denominator in Eq. (53) becomes zero when v = c ( λ, σ ), where c ( λ, σ ) is given by Eq. (50). Then, giventhat v >
0, we obtain from Eq. (53) that v > c ( λ, σ ) ⇔ l > v < c ( λ, σ ) ⇔ l < . (54)Finally, from Eqs. (51) and (54) it follows that v > c a and v < c b . (55)Hence, in all cases, the propagation velocity of the exacttraveling front solution is supersonic with respect to thefluid ahead of the front, and subsonic with respect to thefluid behind it. V. NUMERICAL RESULTS
All numerical calculations rely on a commercially avail-able software package , which is based on the finite ele-ment method. For convenience we introduce the follow-ing non-dimensional variables, denoted by tilde˜ ψ (˜ x, ˜ t ) = 1 b ψ ( x, t ) , ˜ x = c b x, ˜ t = c b t. (56)Under this transformation we may write Eq. (15) as ψ tt − ψ xx = ψ t ψ xx + ∂∂t (cid:18) ψ xx + ( ψ x ) + B/A −
12 ( ψ t ) (cid:19) , (57)where the tildes have been omitted. From a comparisonof Eqs. (15) and (57), we find that the results of theprevious sections subject to b = 1 and c = 1 apply to COMSOL version 3.2a, (2005)
Eq. (57). Non-dimensional versions of the parameters, v , θ , λ , and σ , also indicated by tilde, become˜ λ = λc , ˜ σ = σc , ˜ θ = θc , ˜ v = vc . (58)In the following analysis we consider only the non-dimensional formulation of the problem. For notationalsimplicity we shall omit the tildes. A. Investigation of the front stability criterion
In order to investigate, numerically, the stability prop-erties of the front, we chose as initial condition for thenumerical solution, the exact solution given by Eqs. (25)and (26), and choose v , θ , λ , σ , and B/A such thatEq. (22) is satisfied. For the sake of clarity, we shalllimit the numerical investigations to the specific case of λ = 0, σ = 0, v >
0, and l >
0, which, according toEq. (27), corresponds to fronts propagating to the rightinto an unperturbed fluid.A first numerical simulation is presented in Fig. 3. Ev-idently, the numerical algorithm successfully integratesthe initial condition forward in time. This finding indi-cates that, for the specific choice of parameters, v = 1 . θ = 0 . B/A = 5, the front exists and is stable. Asecond numerical simulation is presented in Fig. 4. Thisinitial condition is given a larger velocity, v = 1 .
7, anda larger height, θ = 0 . B/A = 5 remains un-changed compared to the first example. The parametersare chosen such that Eq. (22) remains satisfied. Appar-ently, the numerical algorithm fails when integrating thesolution forward in time, which indicates that the frontis unstable for the specific choice of parameters. Giventhat σ = 0 and l >
0, the left and right asymptotes ofthe front are given by ψ t = vθ and ψ t = 0, respectively,according to Eq. (27b). Clearly, the right value belongsto the interval (39), thus, it does not causes instabilityof the front. However, if the left value, vθ , lies outsidethe interval (39), it causes instability of the front. Insert-ing the value of B/A from the two examples above intoEq. (39), we find that in the first and second example, vθ lies inside and outside the interval (39), respectively.Hence, the behavior observed in Figs. 3 and 4 agrees withthe stability criterion introduced in Section III B.A large number of numerical simulations have beenperformed in order to systematically investigate the sta-bility properties of the front. Within each simulationthe parameters in the initial condition are, again, chosensuch that Eq. (22) is satisfied. The result of this inves-tigation is presented in Fig. 5. Still, σ = 0 and l > ψ t = vθ . For B/A >
B/A, v )-plane is obtained when vθ equals the upperbound of the interval (39). Using Eq. (29) we obtain vθ = 1 B/A − ⇒ v = (cid:112) ( B/A − B/A + 1)
B/A − . (59)For B/A <
1, the stability threshold curve is given byEq. (30), since vθ lies within (outside) the interval (40)when v < v s ( v > v s ), according to Eq. (29). The twostability threshold curves are included in Fig. 5. The fineagreement between the numerical results and the stabil-ity threshold curves indicates that the stability criterion,introduced in Section III B, is both necessary and suffi-cient in order for the front solution to be stable. −20 −10 0 10 2000.050.1 − ψ x x −20 −10 0 10 2000.050.10.15 ψ t x FIG. 3: The initial condition at t = 0 (bold line) is obtainedfrom Eqs. (25) and (26) subject to v = 1 . θ = 0 . λ = 0, σ = 0, and B/A = 5, which satisfy Eq. (22). The numericalsolutions are shown over the time interval 0 ≤ t ≤ −30 −20 −10 0 10 00.050.10.150.20.25 x ψ t −30 −20 −10 0 100.25930.25940.25950.25960.2597 FIG. 4: The numerical algorithm fails to integrate this solu-tion forward in time (insert shows magnification). See captionof Fig. 3. v = 1 . θ = 0 . λ = 0, σ = 0, and B/A = 5.The numerical solution is shown at time t = 3 . × − . B/A v FIG. 5: Each point in the (
B/A, v )-plane represents a nu-merical simulation, with the initial condition obtained fromEqs. (25) and (26) subject to λ = 0, σ = 0, and θ given byEq. (29). Crosses and circles indicate stable and unstable so-lutions, respectively (compare with Figs. 3 and 4). Solid linesrepresent the stability threshold curves given by Eqs. (30) and(59). B. Head-on colliding fronts
The numerical simulation presented in Fig. 6 shows theresult of a head-on collision between two fronts. From thesimulation we observe that two new fronts emerge uponthe collision. The contour plot reveals that these frontstravel at a higher speed, compared to the speed of thefronts before the collision. For other choices of initialcondition, we found the outcome of the head-on collisionto be fronts traveling at lower speed, compared to thatof the fronts before the collision.In order to analyze solutions of Eq. (15) that com-prise two fronts, we assume that these fronts belong tothe class of exact front solutions derived in Section III Aabove. Investigations of the fronts that emerge upon ahead-on collision have made it clear that this assumptionis true, only when the generalized traveling wave assump-tion is considered, in contrast to the standard travelingwave assumption. Then, for each of the two fronts in thesolution we introduce four new parameters, v , θ , λ , and σ , which must satisfy Eq. (22) as B/A − θ i v i − (1 − ( B/A − σ i ) v i + (cid:18) θ i + 2 λ i (cid:19) v i + σ i + 1 = 0 , i = 1 , , (60)where subscript 1 and 2 denote parameters associatedwith waves positioned to the left and right, respectively.Furthermore, we require that solutions comprising twofronts are continuous and satisfy the following set of ar-bitrary boundary conditions − ψ x → (cid:26) P, x → −∞
Q, x → + ∞ , ψ t → (cid:26) R, x → −∞
S, x → + ∞ . (61)Assuming that l > l <
0, we find, using Eq. (27),that the these requirements lead to the following condi-tions λ = λ = P − θ = Q − θ , (62a) σ = σ = R − v θ = S − v θ , (62b) θ + θ = P − Q, v θ + v θ = R − S. (62c)Then, we substitute the boundary values found in Fig. 6for P , Q , R , and S in Eqs. (62), and substitute thevalue of B/A into Eq. (60). Finally, solving the systemof equations (60) and (62), we obtain the results listedin Table II. The solution in the first row of the tablecorresponds to the initial fronts found in Fig. 6. Thesolution in the second row corresponds to two unstablefronts, according the stability criterion discussed above.The two fronts that emerge upon the head-on collisionare defined by the values found in the third row of thetable. Hence, the fronts after the collision travel at thevelocities − v = v = 1 .
76, which is in agreement withthe velocities determined from the slope of the contourlines in Fig. 6.
VI. CONCLUSIONS
A nonlinear wave equation that governs finite ampli-tude acoustic disturbances in a thermoviscous Newtonianfluid, and includes nonlinear terms up to second order,has been presented. The single dependent variable isthe velocity potential. It has been demonstrated that,in the non-dissipative limit, the equation preserves theHamiltonian structure of the fundamental fluid dynam-ical equations, hence, the model equation is associatedwith corresponding Lagrangian and Hamiltonian densi-ties. Furthermore, we found that the Kuznetsov equationis an approximation of the proposed nonlinear wave equa-tion. However, in the non-dissipative limit the Kuznetsovequation is not Hamiltonian. Exact traveling front solu-tions, for the partial derivatives with respect to space andtime of the dependent variable, has been obtained usinga generalized traveling wave assumption. This general-ized assumption leads to a wider class of exact solutionscompared the one obtained from a standard travelingwave assumption, since the generalized assumption in-cludes two arbitrary constants, which are added to thepartial derivatives. As a result of the generalized trav-eling wave analysis we found that, in order for the frontto exist, its boundary values, its propagation velocity,and the physical parameters of the problem must satisfya given cubic equation in the front propagation veloc-ity. The derivation of the exact solution applies equally t x − ψ x t x ψ t t x FIG. 6: The initial condition (bold lines in the two topmostplots), corresponds to two fronts that make a head-on collisionat t = 42. The initial fronts are defined by v = − v = 1 . θ = − θ = 8 . × − , λ = λ = 0, σ = σ = 0, and B/A = 5, where subscript 1 and 2 relate to fronts positionedthe left and right, respectively. For each of the two fronts theparameters satisfy Eq. (22). Lowermost: contour lines givenby − ψ x = Z , where Z takes 4 equidistantly spaced valuesacross each front. well to the proposed Hamiltonian model equation andthe Kuznetsov equation. Results for both equations havebeen given.It has been demonstrated that the overall stabilityproperties of the front are determined by the stability ofthe two asymptotic tails of the front. A linear stabilityanalysis of these steady parts of the solution revealed thatthe front is stable when the partial derivative with respectto time, which is proportional to the acoustic pressure,belongs to a critical interval, and is otherwise unstable.This stability criterion has been verified numerically, by0 TABLE II: Solution of Eqs. (60) and (62) subject to P = − Q = 8 . × − , R = S = 9 . × − , and B/A = 5 (comparewith Fig. 6).Solution v θ v θ λ = λ σ = σ .
19 8.07 × − − .
19 -8.07 × − − .
25 8.07 × − .
25 -8.07 × − × − − .
76 8.07 × − .
76 -8.07 × − × − using the exact front solution as initial condition in anumber of numerical simulations.It has been demonstrated that, in all cases, the frontpropagates at supersonic speed with respect to the fluidahead of it, while it remains subsonic with respect tothe fluid behind it. The same properties have been re-ported for the classical fluid dynamical shock. Further-more, it has been demonstrated that the cubic equation,mentioned above, is equivalent to the well establishedRankine-Hugoniot relations, which connect the physicalquantities on each side of a shock. However, this resultwas accomplished only when considering the cubic equa-tion obtained from the analysis of the proposed Hamil-tonian wave equation. The generalized traveling waveanalysis based on the Kuznetsov equation is not in agree-ment with the Rankine-Hugoniot relations. Estimates ofthe front thickness may be obtained using the values forthe diffusivity of sound listed in Table. I. In water andair front thicknesses are found to be of the order 10 − and 10 − meters, respectively. However, caution shouldbe taken with these estimates, as the small length scalesviolates the continuum assumption of the governing equa-tions.Numerical simulations of two head-on colliding fronts have demonstrated that two new fronts emerge upon thecollision, and that these fronts, in the general case, travelat speeds, which are different from the speeds of thefronts before the collision. It has been demonstrated thatthe velocities of the fronts after the collision may be cal-culated, based on information about the fronts before thecollision. However, in order to accomplish this calcula-tion, it has proven necessary to introduce the generalizedtraveling wave assumption in the derivation of the frontsolution.In future studies, it would be rewarding to further in-vestigate a variety of interacting fronts, other than thehead-on collision reported in this paper. Also a searchfor other types of wave solutions, might learn us moreabout the properties of the proposed Hamiltonian modelequation and the Kuznetsov equation. Acknowledgments
One of the authors (YuBG) would like to thank theMIDIT Center and Civilingeniør Frederik Leth Chris-tiansens Almennyttige Fond for financial support. [1] S. I. Aanonsen, T. Barkve, J. N. Tjøtta, and S. Tjøtta,“Distortion and harmonic generation in the nearfield ofa finite amplitude sound beam”, J. Acoust. Soc. Am. ,749–768 (1984).[2] V. P. Kuznetsov, “Equations of nonlinear acoustics”, Sov.Phys. Acoust. , 467–470 (1971).[3] E. A. Zabolotskaya and R. V. Khokhlov, “Quasi-planewaves in the nonlinear acoustics of confined beams”, Sov.Phys. Acoust. , 35–40 (1969).[4] P. M. Jordan, “An analytical study of Kuznetsov’s equa-tion: diffusive solitions, shock formation, and solutionbifurcation”, Phys. Lett. A , 77–84 (2004).[5] Y. Jing and R. O. Cleveland, “Modeling the propagationof nonlinear three-dimensional acoustic beams in inho-mogeneous media”, J. Acoust. Soc. Am. , 1352–1364(2007).[6] J. W´ojcik, “Conservation of energy and absorption inacoustic fields for linear and nonlinear propagation”, J.Acoust. Soc. Am. , 2654–2663 (1998).[7] J. Hoffelner, H. Landes, M. Kaltenbacher, and R. Lerch,“Finite element simulation of nonlinear wave propaga-tion in thermoviscous fluids including dissipation”, IEEETrans. Ultrason., Ferroelect., Freq. Contr. , 779–786 (2001).[8] A. Shermenev, “Separation of variables for the nonlin-ear wave equation in cylindrical coordinates”, Physica D , 205–215 (2005).[9] B. O. Enflo and C. M. Hedberg, Theory of Nonlin-ear Acoustics in Fluids (Kluwer Academic, Dordrecht)(2002).[10] S. Makarov and M. Ochmann, “Nonlinear and thermo-viscous phenomena in acoustics, part II”, Acustica ,197–222 (1996).[11] M. V. Aver’yanov, V. A. Khokhlova, O. A. Sapozhnikov,P. Blanc-Benon, and R. O. Cleveland, “Parabolic equa-tion for nonlinear acoustic wave propagation in inho-mogeneous moving media”, Acoust. Phys. , 623–632(2006).[12] K. Naugolnykh and L. Ostrovsky, Nonlinear Wave Pro-cesses in Acoustics (Cambridge University Press, Cam-bridge) (1998).[13] L. H. S¨oderholm, “A higer order acoustic equation forthe slightly viscous case”, Acustica , 29–33 (2001).[14] V. E. Zakharov and E. A. Kuznetsov, “Hamiltonian for-malism for nonlinear waves”, Physics-Uspekhi , 1087–1116 (1997). [15] P. M. Jordan, “Bifurcations of diffusive soliton solu-tions to Kuznetsov’s equation”, J. Acoust. Soc. Am. ,2283–2283 (2003).[16] S. Makarov and M. Ochmann, “Nonlinear and thermo-viscous phenomena in acoustics, part I”, Acustica ,579–606 (1996).[17] R. T. Beyer, “The parameter B/A ”, in
Nonlinear Acous-tics , edited by M. F. Hamilton and D. T. Blackstock,chapter 2, 25–39 (Academic Press, San Diego) (1998).[18] M. F. Hamilton and C. L. Morfey, “Model equations”,in
Nonlinear Acoustics , edited by M. F. Hamilton andD. T. Blackstock, chapter 3, 41–64 (Academic Press, SanDiego) (1998).[19] D. R. Lide, ed.,
CRC Handbook of Chemistry andPhysics, Internet Version 2007, (87th Edition), (Taylor and Francis, Boca Ra- ton) (2007).[20] R. Gent,
Applied Physics and Technology of DiagnosticUltrasound (Milner Publishing, Prospect) (1997).[21] I. Christov, C. I. Christov, and P. M. Jordan, “Model-ing weakly nonlinear acoustic wave propagation”, Q. J.Mechanics. Appl. Math. , 473–495 (2007).[22] H. Goldstein, C. Poole, and J. Safko, Classical Mechan-ics , 3 edition (Addison Wesley) (2001).[23] M. F. Hamilton, D. T. Blackstock, and A. D. Pierce,“Progressive waves in lossless and lossy fluids”, in
Nonlin-ear Acoustics , edited by M. F. Hamilton and D. T. Black-stock, chapter 3, 65–150 (Academic Press, San Diego)(1998).[24] L. D. Landau and E. M. Lifshitz,