A multiscale approach to liquid flows in pipes I: The single pipe

A multiscale approach to liquid flows in pipes I: the single pipe Andrea Corli ∗ Ingenuin Gasser Arne Roggensack † § Mária Lukáčová-Medvid’ová Ulf Teschke ‡ ¶ July 8, 2012 AMS(MOS) subject classifications: 76Q05, 35L02, 35L05, 65M08 Keywords: pipe flow, Saint-Venant equations, multiscale analysis, water hammer, pressure waves Abstract In the present paper we study the propagation of pressure waves in a barotropic flow through a pipe, with a possibly varying cross-sectional area. The basic model is the SaintVenant system. We derive two multiscale models for the cases of weak and strong damping, respectively, which describe the time evolution of the piezometric head and the velocity. If the damping is weak, then the corresponding first-order hyperbolic system is linear but contains an additional integro-differential equation that takes into account the damping. In the case of strong damping, the system is nonlinear. The full and multiscale models are compared numerically; we also discuss results obtained by a largely used commercial software. The numerical experiments clearly demonstrate the efficiency of the multiscale models and their ability to yield reliable numerical approximations even for coarse grids, that is not the case for the full model. 1 Introduction In this paper we study a mathematical model for the liquid flow in a single pressurized pipe; a forthcoming paper shall be concerned with pressurized pipe networks. The model is based on a hyperbolic system of two balance laws, the Saint-Venant or shallow water equations, in one spatial dimension; it takes into account both the friction of the liquid at the pipe walls, the elastic behavior of the pipes and the gravity effects on the flow. Initial and boundary conditions are prescribed. ∗ Department of Mathematics, University of Ferrara, Ferrara, Italy Fachbereich Mathematik, Universität Hamburg, Hamburg, Germany ‡ Institut für Mathematik, Johannes Gutenberg-Universität Mainz, Mainz, Germany § Fachbereich Mathematik, Universität Hamburg, Hamburg, Germany ¶ Department of Mechanical Engineering, University of Applied Science Hamburg and IMS Ingenieurgesellschaft mbH, Hamburg † 1 The standard approach, when one only considers the flow behaviour of a liquid in a pipe, is to use an incompressible model. For that purpose, the effects due to compressibility can be neglected in typical applications. However, an incompressible model assumes an infinite sound speed and therefore it is inappropriate for studying pressure (sound) waves in liquids. In this paper we consider pressure waves in liquids. The focus is on a multiscale approach, in order to emphasize the different time-regimes related to the fluid motion and to the pressure waves. The latter are well known in hydraulic engineering because they may cause the so-called water hammer effect, leading even to the breakage of the pipes, see [5, 10, 16, 28]. This effect is the consequence either of a sudden closure of a valve downstream or of a rapid variation of the pressure upstream: a pressure wave is then generated, travelling at a speed usually much larger than that of the fluid. In case of a closed pipe, the pressure wave may bounce many times before being damped by the fluid. In our multiscale approach we model such phenomenon by deducing an asymptotic system from the complete one. The starting point of our analysis is the initial-boundary value problem in a stripe 0 ≤ x ≤ 1 for the scaled Saint-Venant system  (Aρ)t + (Aρu)x = 0 (1.1) (Aρu)t + (Aρu2 )x + 12 Apx = −ζAρu|u| − Fr12 Aρzx . We refer to Section 2 for more information on the scaling. Above, ρ is the density and u the velocity; A = A(t, x) and z = z(x) denote the cross-sectional area and the height profile of the axis of the pipe, respectively. With  and ζ we denote two parameters lumping together several physical data; roughly speaking, they account for dynamical and friction effects, respectively. We consider a turbulent flow, and then the parameter ζ does not depend on the state variables. Moreover, Fr is the Froude number. As pressure law we consider p ∼ p0 + C ln ρ, for a positive constant C (see below). In the case that both A and z are constant, the existence and uniqueness of global W 1,∞ solutions to the Cauchy problem for (1.1) was first proved in [20], for a general pressure law and suitably small initial data; moreover, the friction parameter ζ was allowed to depend on the momentum q = ρu, in order to take into account both laminar (i.e., when ζ ∼ 1/|q|) and turbulent flows (when ζ is constant). Subsequently, the global existence of weak solutions for the initial-boundary value problem was proved in [21], in the case of a linear pressure law and boundary conditions ρ(0, t) = ρB (t) and q(L, t) = 0; a suitable bound on the initial data guaranteed that the flow was subsonic. The case of laminar flows has been long studied both for blood flows in arteries, see for instance [13], and for applications to porous media equations [22]. About related results on the initial-value problem for the gas flow in a pipe we quote [19]; for simpler models based on a p-system, we refer to [15] for pipes with elbows and [3, 2], [6, 7, 8] for problems with junctions. In this paper, we introduce a fast-time and a slow-time scale for (1.1), in order to separate the propagation of the sound waves, occurring in the former, from the fluid dynamics, occurring in the latter. Moreover, two different asymptotics are considered, according to the relevance of the damping effects. In the case of weak damping (ζ << −1 ) we obtain two wave equations coupled only through the initial and boundary conditions; when the damping is sufficiently strong (ζ ∼ −1 ) we obtain instead a semilinear system. Both cases are then studied; remark that we are always in the turbulent regime prescribed by (1.1). All of that is the content of Section 3; the previous Section 2 introduces the full models in the different settings. No rigorous 2 result about the global existence or the decay of the solutions is given, since this is not the aim of the paper. Indeed, both problems are rather stiff and many issues are still left open; for instance, see Section 3.2 for a short discussion of a damped wave equation arising in the asymptotics. As a consequence, also the asymptotic expansions are merely formal. Several numerical simulations are provided; they are contained in Section 4. All of them show a very good agreement between the full and the asymptotic models. Moreover, we compare some results of ours with the simulations obtained by the commercial software package MIKE URBAN [11]. Once more, the agreement is very good. Other numerical experiments can been found, for instance, in [24, 14], in a simplified setting, and in [25]. Let us mention that many commercial softwares are at disposal for this problem, see e.g. AFT-Impulse [1] or MOUSE (MOdel for Urban SEwers) [11]. 2 The model and the scaling The model. We consider the general equations of fluid-dynamics in a pipe in one spatial dimension. We assume that the pipe has circular cross-section with internal diameter D̃ and area Ã; different cross-sections may be treated analogously by introducing the notion of hydraulic diameter. The diameter (and the area) may vary both because of different sizes of the pieces composing the pipe and because of elastic effects due to the pressure of the flow. The balance equations for the mass and momentum read [5] ( (Ãρ̃)t̃ + (Ãρ̃ũ)x̃ = 0, (2.1)
(Ãρ̃ũ)t̃ + (à ρ̃ũ2 + p̃) x̃ = Ãx̃ p̃ − 2λD̃ Ãρ̃ũ|ũ| − Ãg ρ̃z̃x̃ , where t̃, x̃ are real variables and ρ̃, ũ, p̃ are the unknown fluid density, velocity and pressure; moreover, we denote q̃ = ρ̃ũ. The function z̃ = z̃(x̃) is the height profile of the axis of the pipe while the constants λ and g are the Darcy friction factor and the gravity acceleration, respectively. The parameter λ, which depends on the diameter, also accounts for viscosity and it is usually computed by means of the Colebrook equation or using a Moody diagram. The system (2.1) is closed by a constitutive law p̃ = p̃(ρ̃) for the pressure, with p̃0 (ρ̃) > 0. In the case that à is a known function of x̃, possibly constant, and p p̃00 (ρ̃) 6= 0, the system p (2.1) is strictly hyperbolic with two genuinely nonlinear eigenvalues ũ ± p̃0 (ρ̃); the term p̃0 (ρ̃) is called celerity in the applied literature. In the current case, where the fluid is a liquid, we assume the barotropic law [28], [16], ρ̃ (2.2) p̃ = pr + K ln , ρr where K is the elasticity module (or adiabatic compressibility module) of the fluid and ρr , pr are constant p reference quantities. Under (2.2) and assuming that A is constant, the eigenvalues are ũ ± K/ρ̃ for ρ̃ > 0. Remark that the greater is ρ̃, the less is the celerity; this depends on the fact that the pressure law (2.2) is concave (as several stress-strain relations in elasticity) instead of being convex (as for gases). Other constitutive laws can be considered, for instance the empiric (convex) power law [27, (1.19)]  n ρ̃ − Bpr . (2.3) p̃ = pr (B + 1) ρr The pressure laws (2.2) and (2.3) coincides at first order if n = 1 and B = K/pr − 1. 3 In addition to the initial conditions, in the case of a pipe modelled by the interval [0, L] boundary conditions must be taken into account. For simplicity, we always consider the subsonic p case |ũ| < p̃0 (ρ̃); other cases can be treated as well. In such a case, two boundary conditions at the ends of the pipe must be assigned. Usually one assigns either the pressure p̃, to model inflow or outflow at a reservoir, or the discharge Q̃ = Ãũ, to model drawings. In particular, a null drawing models a pipe closed at that end. The solutions under consideration are intended to be smooth; compatibility conditions between the initial and the boundary data are always assumed. The scaling. In addition to ρr and pr introduced above we use reference quantities xr t r , x r , ur = tr and for a generic function f we define f˜(x̃, t̃) = fr · f ( xx̃r , tt̃r ). Reference quantities and their 2 typical values in the case of water are given in Table 1. Remark that since à = √π(D̃) /4, and the same relation holds for Ar and Dr , for the scaled quantities we have D = A. Moreover, we introduce the parameters K 1 , = 2 ρr u2r ζ=λ L , 2Dr 1 gL sin θr . = u2r Fr2 Here, θr is a reference angle between the axis of the pipe and a horizontal line while Fr is the Froude number, the ratio between the inertial and gravitational forces; remark that ur > 0, see Table 1. If the height profile z̃ is constant, then the term containing z̃x in (2.1) is missing and the introduction of Fr is not needed; if z̃ is not constant we assume that zr > 0. We refer to Table 2 for the characteristic values of the dimensionless parameters. We then obtain the scaled model (see (1.1)) ( (Aρ)t + (Aρu)x = 0, (Aρu)t + (Aρu2 )x + 1 2 Apx √ = −ζ Aρu|u| − 1 Aρzx , Fr2 (2.4) and (2.2) becomes p = 1 + ln ρ. (2.5) Analogously for (2.3) we find p = (B + 1)ρn − B. The scaled system (2.4) can be written in terms of p and u by taking into account (2.5): ( pt + ux + upx = − A1 (At + uAx ) , (2.6) ut + uux + 12 ρ1 px = − √ζA u|u| − Fr12 zx . A different writing is obtained by introducing the piezometric head H =p−1+ 2 2 2 z = ln ρ + z, Fr2 Fr2 (2.7) 2  so that ρ = eH− Fr2 z and Ht = pt , Hx = px + Fr 2 zx . This gives the equations in terms of H, u: ( 2 Ht + ux + uHx = − A1 (At + uAx ) +  uzx  , Fr2 (2.8) ζ 1 1 1 1 ut + uux + 2 ρ Hx = − √A u|u| − Fr2 1 − ρ zx . Since we deal with smooth solutions, boundary conditions for (2.8) are straightforwardly deduced by those imposed to (2.4). 4 Inelastic pipes. We simplify the previous systems by assuming that the pipe is perfectly inelastic. This means that A = A(x) is a known function, possibly constant. Then (2.4) reads ( ρt + (ρu)x = − AAx ρu, (ρu)t + (ρu2 )x + 1 p 2 x = − √ζA ρu|u| − 1 ρzx , Fr2 (2.9) to be studied, for instance, under the initial-boundary conditions (recall that we assume a subsonic flow)   p(0, t) = pl (t) , p(1, t) = pr (t) , ρ(x, 0) = ρ∗ (x) , u(x, 0) = u∗ (x) , and or (2.10)  p(0, t) = pl (t) , u(1, t) = 0 . The first set of boundary conditions prescribes inflows (outflows) at the boundary, while the condition on the velocity represents a pipe closed at that end. In turn, by writing Q for u, in order to use the same notation as in many hydraulic textbooks, the system (2.8) reads  2  Ht + Qx + QHx = − AAx Q + Fr 2 Qzx ,   (2.11)  Qt + QQx + 12 1 Hx = − √ζ Q|Q| − 12 1 − 1 zx . ρ  ρ Fr A The eigenvalues of (2.14) are Q ± 1 √  ρ. H(x, 0) = H∗ (x) , Q(x, 0) = Q∗ (x) , Initial-boundary conditions for (2.14) are as above:   H(0, t) = Hl (t) , H(1, t) = Hr (t) , or  H(0, t) = Hl (t) , Q(1, t) = 0 . and (2.12) . Remark that the introduction of the discharge Q̃ = ρu, instead of Q, as is common in the applied hydraulic literature, makes the equations more cumbersome. In the special case of pipes of constant cross section we take A = 1. Then, (2.9) reads ( ρt + (ρu)x = 0, (2.13) (ρu)t + (ρu2 )x + 12 px = −ζρu|u| − Fr12 ρzx , while system (2.11) becomes   Ht + Qx + QHx =  Qt + QQx + 1 1 H 2 ρ x 2 Qzx , Fr2 = −ζQ|Q| − 1 Fr2   1 − ρ1 zx . (2.14) Elastic pipes. The assumption of constant cross section of the pipes is not realistic in some cases; for simplicity we discuss only the case when the area of the pipe at rest is constant. The elasticity of the pipe plays a fundamental role in the study of the waterhammer effect, which is mostly concerned with pipes in hydraulic plants. Such effect, however, has some interest also in urban networks, for instance when a pump providing water to a reservoir is switched off or when a valve in a reservoir is closed. Moreover, the waterhammer effect has been proved to be a useful tool for leak detection in closed conduits [4]. 5 The case of elastic pipes is modelled for instance in [29, eqn. (34), page 17], where the Mariotte law 1 D̃ 1 dà . dp̃ =φ , · , (2.15) for φ = ER s̃ dt̃ à dt̃ is postulated to hold for thin iron pipes. Here, ddt̃ = ∂t̃ + ũ∂x̃ denotes the total derivative, ER the elasticity module of the pipe material and s̃ the thickness of the pipe wall. Moreover, it is assumed that only the derivatives of the area are sufficiently large; therefore, D̃ is assumed here to equal Dr . For the same reason, also the parameter λ is assumed to depend only on Dr . Denoting 1 , (2.16) σ= 1 + φK then, under (2.15), the system (2.1) may be written as ( ρ̃t̃ + (ρ̃ũ)x̃ − (1 − σ)ρ̃ũx̃ = 0 (2.17) (ρ̃ũ)t̃ + (ρ̃ũ2 + σ p̃)x̃ − (1 − σ)ρ̃(ũt + (u2 )x ) = −σ 2λD̃ ρ̃ũ|ũ| − σgρ̃z̃x̃ . √ p Under (2.2), the eigenvalues of (2.17) are ũ± σ · K/ρ̃ see [29, eqn. (37), page 18]: the celerity √ is diminished of the factor σ with respect to the case when à is constant. In this case the scaled system becomes, compare with (2.4), ( ρt + (ρu)x − (1 − σ)ρux = 0, (2.18) (ρu)t + (ρu2 )x + σ2 px − (1 − σ)ρ(ut + 2uux ) = −σζρu|u| − Frσ2 ρzx . In turn, introducing p, the system (2.6) becomes ( pt + σux + upx = 0 , ut + uux + 1 1  2 ρ px At last, system (2.8) now reads   Ht + σQx + QHx = 3  Qt + QQx + 1 1 2 ρ Hx = −ζu|u| − 1 z . Fr2 x 2 Qzx , Fr2 = −ζQ|Q| − 1 Fr2   1 − ρ1 zx . (2.19) (2.20) Asymptotics In this section we analyze the asymptotic behaviour of the system (2.14). We have the x-space scale, and introduce a slow-time scale t and a fast-time scale τ = t/; then we look for H and Q solving (2.14) under the form H = H(x, t, τ ), Q = Q(x, t, τ ). We distinguish the case of weak damping, where ζ is small, from that of strong damping, where ζ ∼ α/. Two different asymptotic expansion are proposed. 3.1 The case of weak damping In this subsection we consider the system (2.14) in the case the damping coefficient ζ is fixed and, in particular, does not depend on . We expect the sound (or pressure) waves on the fast scale and the (incompressible) fluid dynamics on the slow scale. For simplicity, we first give the detail about the asymptotics in the case A ≡ 1; then we provide the equations for the general inelastic case A = A(x) and the elastic case (when the area at rest is constant). 6 The case A ≡ 1. We introduce the following asymptotic expansions into the system (2.14): H(x, t, τ ) = H1 (x, τ ) + 2 H2 (x, t) + O(3 ) , 2 Q(x, t, τ ) = Q0 (x, t, τ ) + Q1 (x, t) + O( ) . (3.1) (3.2) We now make some remarks about these expansions. First, no term H0 is present in (3.1), because in most real-world applications the scaled piezometric head is not expected to be of order 1; indeed, in our scaling the reference value for H is 0. Second, we model the fast dynamics at the higher orders in  on a τ scale and the slow dynamics at lower order in  on a t scale. This is why H1 depends on τ , while H2 and Q1 depend on t. However, we look for Q0 under the form Q0 (x, t, τ ) = Q0 (x, τ ) + S(t) , (3.3) allowing the fluid motion to react to sudden changes in the pressure on a long-time scale at first order. The determination of the functions Q0 and S shall be done through the initial conditions and an averaging argument, see also the Remark 3.1. Indeed, by plugging the expansions (3.1)– (3.2) into (2.14), it is possible, for any for k = 0, 1, . . ., to deduce recursively a hierarchy of hyperbolic systems coupling Hk+1 (x, t, τ ) with Qk (x, t, τ ), where the already determined terms Hk , Hk−1 , . . . appear as well as Qk−1 , Qk−2 , . . .. Our aim is not in the determination of all these profiles, that would be easy from an analytical point of view but uninteresting for numerical simulations. Here, we focus mainly on the fast-scale terms H1 and Q0 ; the role of the slow-scale term S, which solves an approximate equation, is of lumping together some second-order effects and an averaged behaviour of both H1 and Q0 . We assume that the initial data in (2.12) are decomposed into H∗ (x) = H10 (x) + 2 H20 (x) + O(3 ) , 2 Q∗ (x) = Q00 (x) + Q10 (x) + O( ) , (3.4) (3.5) for 0 ≤ x ≤ 1. Concerning the boundary data we focus on the case of a pipe with an inflow at one end and an outflow at the other end; the other case follows as well. We assume that the boundary data in (2.12) satisfy Hl (t) = H1l (τ ) + 2 H2l (t) + O(3 ) , Hr (t) = H1r (τ ) + 2 H2r (t) + O(3 ) . Moreover, we fix the value of the constant initial velocity of the (slow) fluid-dynamic flow, S(0) = S0 . (3.6) 2 Next, we use ρ = eH− Fr2 z = 1 + H1 + O(2 ) and plug (3.1)-(3.2) into (2.14); we obtain  2 3  H1τ + Q0x +  Q0 H1x = Fr 2 Q0 zx + O( ) ,     1 1 (3.7) Q + Q + Q + Q Q + (1 − H ) = H + H 1τ 1 1x 2x 0t 0τ 0 0x      = −ζQ0 |Q0 | − Fr12 H1 (1 − H1 )zx + O(2 ) . This gives a hierarchy of equations. We only extract the order 0 equation from the first equation in (3.7) and order −1 , 0 equations from the second; we obtain   H1τ + Q0x = 0 , (3.8) Q + H1x = 0 ,  0τ Q0t + Q0 Q0x − H1 H1x + H2x = −ζQ0 |Q0 | . 7 Remark that Q1 does not appear here because of the ansatz on Q. Let us first focus on the first two equations in (3.8), which write, by (3.3),  H1τ + Q0x = 0 , (3.9) Q0τ + H1x = 0 . They form a hyperbolic system with propagation speeds ±1. By recalling (3.3) and (3.6), we prescribe the initial and boundary conditions for (3.9) as H1 (x, 0) = H10 (x) , H1 (0, τ ) = H1l (τ ) , Q0 (x, 0) = Q00 (x) − S0 , H1 (1, τ ) = H1r (τ ) . The first two equations in (3.8) lead to a wave equation for H1 :   H1τ τ = H1xx , H1 (0, τ ) = H1l (τ ) , H1 (1, τ ) = H1r (τ ) ,  H1 (x, 0) = H10 (x) , H1τ (x, 0) = −Q00x (x) . (3.10) (3.11) (3.12) The initial value for the fast-time derivative of H1 is given by the first equation (3.9). As for the problem (2.14), standard compatibility conditions must be assigned on the initial and boundary data in (3.12), which are deduced by those imposed on the data of (2.14). In the same way, from the first two equations in (3.8) we are led to a wave equation for Q0 :   Q0τ τ = Q0xx , Q0x (0, τ ) = −H1lτ (τ ) , Q0x (1, τ ) = −H1rτ (τ ) , (3.13)  Q0 (x, 0) = Q00 (x) − S0 , Q0τ (x, 0) = −H10x (x) . The initial value for the fast time derivative of Q0 is given by the second equation in (3.9). Notice that the damping coefficient ζ is missing in both (3.12) and (3.13). Therefore H1 is not damped in this asymptotic regime; on the other hand, S shall depend on ζ, and this gives a long-time damping for Q0 . The existence of solutions to the initial-boundary value problems (3.12) and (3.13), is classical and can be proved by a straightforward integration along the characteristics. In turn, this establishes the existence of solutions to (3.9) with initial data (3.10) and boundary data (3.11). Now, we exploit the third equation in (3.8) to deduce an averaged approximate equation for S. We integrate (3.8)3 with respect to x from 0 to 1; denoting ∆f = f (1) − f (0), for any function depending on x, we obtain Z 1

1 2 2 ∆H1 − ∆Q0 − St = −∆Q0 S − ∆H2 + ζ(Q0 + S)|Q0 + S| dx , (3.14) 2 0 where we still have τ dependent terms. However, notice that the functions H1 and Q0 solving (3.12) and (3.13) are not damped. In the current setting, where only boundary conditions which are strongly localized in time are given, both H1 and Q0 are represented by single waves bouncing back and forth in the pipe. Since the propagation speed of both wave equations is 1 and the domain is the interval [0, 1], we obtain an approximate integro-differential equation for S by averaging (3.14) with respect to the τ variable:  Z Z Z Z 1
1 2 2 (Q0 + S)|Q0 + S| dx dτ (3.15) St = −∆H2 + − ∆H1 − ∆Q0 dτ − −∆Q0 dτ S − ζ − 2 0 R2 R combined with (3.6). Here, we denoted − dτ = 12 0 dτ . To summarize, we solve 8 1. the initial-boundary value problem (3.12) for H1 , 2. the initial-boundary value problem (3.13) for Q0 , 3. the initial value problem for the integro-differential equation (3.15) for S, and obtain Q0 = Q0 + S. Remark 3.1 The solution process outlined above is invariant with respect to the transformation Q0 (x, τ ) → Q0 (x, τ ) − T S(t) → S(t) + T . This can easily be checked using the properties of (3.13) and (3.17). This is important and tells us that it is not important how the constant part of the initial velocity is divided between fast and slow problem. In other words, Q0 is invariant by changes of S0 . Remark 3.2 If instead of the expansion (3.1) we simply considered H(x, t, τ ) = 2 H2 (x, t) + O(3 ) , then, in the leading order, we should obtain the incompressible equations  Q0x = 0 , Q0t + (H2 )x = −ζQ0 |Q0 | . (3.16) However, this asymptotics is too rough because in (3.16) information travels with infinite speed (since the pressure correction H2 satisfies the linear elliptic equation H2xx = 0). Example 3.1 (Incompressible flow) Assume A ≡ 1 and that for some S0 ∈ R we have H10 (x, 0) ≡ 0, H1l (τ ) ≡ H1r (τ ) ≡ 0, Q00 (x) ≡ S0 . Then Q00x (x) ≡ 0 and H1 (x, τ ) ≡ 0, Q0 (x, τ ) ≡ 0, Q0 (x, t, τ ) = S(t), where S solves the following initial-value problem for the ordinary differential equation  St = −ζS|S| − ∆H2 , (3.17) S(0) = S0 . Lemma 3.1 For every S0 and ∆H2 ∈ R the initial-value problem (3.17) has a unique solution S ∈ C 1 (R). Moreover, for t → ∞,  q q √ |∆H2 |  2|  − sgn(∆H2 ) |∆H 6 0, + O(1) · e−2 |∆H2 | t , if ∆H2 = ζ ζ S(t) ∼ 1  , if ∆H2 = 0 .  ζt The proof follows by a straightforward computation and is therefore omitted. Example 3.2 (Sound waves in a fluid at rest) Assume again A ≡ 1 and H2l (t, τ ) ≡ H2r (t, τ ) ≡ 0, Q00 (x) ≡ S0 = 0. Then S(t) ≡ 0. Moreover, H1 solves the problem (3.12) with initial datum H1τ ≡ 0 and does not depend on Q0 . In turn, Q0 ≡ Q0 solves the problem (3.13) with initial datum Q00 (x) ≡ 0. In this case, the fluid velocity is different from zero only where the sound waves are supported and vanish otherwise. 9 The case A = A(x). We consider now the case of a general inelastic pipe, namely A = A(x). If Ax 6= 0, the damping effect of this term appears at order zero and does not allow the presence of the term S. We consider then the simplified asymptotics H(x, t, τ ) = H1 (x, τ ) + O(2 ) , Q(x, t, τ ) = Q0 (x, τ ) + O() , Then, the first two equations in (3.8) become  H1τ + Q0x = − AAx Q0 , Q0τ + H1x = 0 . (3.18) However, in this case the wave equations obtained as in (3.12) and (3.13) do not decouple, and then we must consider directly (3.18). The elastic case. We now consider the case of an elastic pipe, namely the system (2.20). In this case the full asymptotic expansions (3.1)–(3.2). We first deduce the system  H1τ + σQ0x = 0 , (3.19) Q0τ + H1x = 0 , whence the wave equations   H1τ τ = σH1xx , H (0, τ ) = H1l (τ ) , H1 (1, τ ) = H1r (τ ) ,  1 H1 (x, 0) = H10 (x) , H1τ (x, 0) = −σQ00x (x) , and   Q0τ τ = σQ0xx , Q (0, τ ) = − σ1 H1lτ (τ ) , Q0x (1, τ ) = − σ1 H1rτ (τ ) ,  0x Q0 (x, 0) = Q00 (x) − S0 , Q0τ (x, 0) = −H10x (x) . The equation (3.15) for S now reads   Z  Z Z Z 1 1 1 1 ζ 2 2 St = −∆H2 + − ∆H1 − ∆Q0 dτ − −∆Q0 dτ S − − (Q0 + S)|Q0 + S| dx dτ . 2 σ σ σ 0 (3.20) 3.2 The case of strong damping In this subsection we still consider the system (2.14) but in the case of strong damping ζ= α ,  (3.21) for some positive α. Indeed, in some cases, for instance for long pipes, the parameter ζ cannot be discarded, and a different asymptotics must be introduced. We first focus on the case A ≡ 1 and introduce then the asymptotic expansions, analogous to (3.1)–(3.2), H(x, t, τ ) = H1 (x, τ ) + O(2 ) , (3.22) Q(x, t, τ ) = Q0 (x, τ ) + O() , (3.23) 10 into the system (2.14). We do not consider higher-order asymptotic expansions as (3.1)–(3.2) since the equations for the higher-order terms seems difficult to justify; moreover, asymptotic expansions as (3.1)–(3.2) lead to the same system (3.24) below. We plug (3.22)-(3.23) into (2.14), keeping the terms of order 0 in (2.14)1 and those of order −1 in (2.14)2 . We find  H1τ + Q0x = 0 , (3.24) Q0τ + H1x = −α|Q0 |Q0 . This system is again considered under initial and boundary conditions (3.10)-(3.11). When it is written in term of two wave equations, the system (3.24) decouples into a damped wave equation for Q0 and a wave equation with source term for H1 . By using (3.4)–(3.6) to deduce the initial and boundary conditions we obtain   Q0τ τ − Q0xx = −2α|Q0 |Q0τ , 0 (τ ) , Q (1, τ ) = −H 0 (τ ) , (3.25) Q (0, τ ) = −H1l 0x 1r  0x 0 (x) − α|Q (x)|Q (x) , Q0 (x, 0) = Q00 (x) , Q0τ (x, 0) = −H10 00 00 and   H1τ τ − H0xx = 2α|Q0 |Q0x , H (0, τ ) = −H1l (τ ) , H1 (1, τ ) = H1r (τ ) ,  1 H1 (x, 0) = H10 (x) , H1τ (x, 0) = −Q000 (x) . (3.26) Differently from the previous case of weak damping, one first solves (3.25) and then (3.26). From a theoretical point of view, to our best knowledge no result seems at disposal about the global existence and the decay of solutions to (3.25). In the case that one boundary data is zero, the proof in [23] applies without changes, but gives no information on the decay of the solutions. If both boundary data vanish, we quote [12], where however the decay is proved. We refer also to [9] for the linear case, when the damping is a(x)Q0τ for a suitable positive function a, and zero Dirichlet conditions. Moreover, the theory developed in [18], see e.g. §5, Theorem 1.1, for smooth solutions in a general nonlinear case does not apply to our problem, even for small data. In fact, in our case the minimal characterizing number of the boundary data equals 1 (instead of being strictly less than 1), meaning that the boundary conditions are not sufficiently dissipative. The case A = A(x). As above, we plug (3.22)–(3.23) into (2.11) and obtain  A   H1τ + Q0x = − x Q0 , A 1   Q0τ + H1x = −α √ |Q0 |Q0 , A whence we deduce the corresponding equation   1 Ax Ax Q0x − 2α √ |Q0 |Q0τ . Q0 + Q0τ τ − Q0xx = A x A A The elastic case. The elastic case follows by plugging (3.22)–(3.23) into (2.20). We obtain  H1τ + σQ0x = 0 , Q0τ + H1x = −α|Q0 |Q0 , and then Q0τ τ − σQ0xx = −2α|Q0 |Q0τ . 11 4 Numerical experiments In this section, through several numerical simulations we show the behavior of the full models (2.14) and (2.20), as well as that of the multiscale models considered above, both in the case of weak and strong damping. First, we describe a numerical method used for the full system. The nonlinear system (2.14) is nonconservative and will be discretized by the second order Lax-Wendroff method in the following way un+1 := unj − j  
 ∆t2 2 n
 n ∆t n n n A unj unj+1 − unj−1 + A u − 2u + u u j+1 j j−1 + ∆tS(uj ). j 2∆x 2∆x2 Here, unj = (Hjn , Qnj ) denotes average values of the piezometric head and the velocity in the mesh cell j, at time step n, while S(u) is the source term on the right hand side of (2.14). The matrix A(u) is the Jacobian matrix of the system, i.e.,   Q 1 A(u) = 1 . (4.1) Q 2 ρ The resulting scheme is explicit in time and requires a restriction on the time step according to the CFL stability condition ∆t ≤ CFL; ∆x max |λ(A(unj ))| j CFL ∈ (0, 1]. In all our calculations we have taken the CFL number to be 0.9. We set z ≡ 0 and take the following initial conditions H∗ (x) = 0, Q∗ (x) = 0, (4.2) and boundary conditions (see (2.12))  open pipe  Hr (t) = 0, Hl (t) = bto (t), or  Qr (t) = 0, closed pipe (water hammer problem). (4.3) Here above, for to = 0.1 · , the term bto (t) = ( − e 0 (t−to )2 2 t2 o −(t−to ) if 0 ≤ t < 2to , otherwise, models a smooth bump of duration 2to with the maximum at to . In some simulations below we will vary the cross section. To this aim we consider two cases, corresponding to an increasing and a decreasing cross section:   if x ∈ [0, 31 ], if x ∈ [0, 31 ],   2 1 A2 (x) = f2 (x) if x ∈ ( 13 , 32 ), (4.4) A1 (x) = f1 (x) if x ∈ ( 13 , 32 ),     2 2 2 if x ∈ [ 3 , 1], 1 if x ∈ [ 3 , 1], for    1 3 1 , f1 (x) = − cos 3π x − 2 2 3    2 3 1 −x . f2 (x) = − cos 3π 2 2 3 12 We consider pipes whose length is either L = 100m (hydraulic plant) or L = 1000m (urban network). The computational domain [0, 1] has been divided into N = 1000 grid cells, except in a few examples where we try to reduce this number in order to test the quality of the scheme. We set the final time either to Tf =  (the waves travel once through the pipe) or to Tf = 10  for long-time simulations (in this case the waves travel ten times through the pipe). The size of  is determined by the typical velocity, which is different in the cases of hydraulic plants and urban networks. Now, we provide details both for the numerical schemes used in the multiscale models and the numerical experiments that we performed; moreover, we compare the results obtained for the full and multiscale models. In addition, for some examples we were able to compare our results with those obtained by using the commercial software package MIKE URBAN [11]1 . Example 4.1 (Weak damping) In this experiment we focus on the case of weak damping for both the full model (2.14) (or (2.20) in the elastic case) and the multiscale model (3.8) with (3.15) (or (3.20) in the elastic case). The initial and boundary conditions, which are deduced by (4.2)–(4.3), are the following, for τo = 0.1: H10 (x) = H20 (x) = 0, H1l (τ ) = bτ0 (τ ), Q00 (x) = 1, H1r (τ ) = 0, H2l (t) = H2r (t) = 0, S0 = 1. Notice that the system (3.9) is linear hyperbolic; as a consequence, we apply a suitable second order finite-volume method [17]. In particular, the so-called flux-vector splitting is used in order to approximate the flux function F := (Q0 , H1 )T , cf. (3.9). Thus, the numerical flux function is H(uL , uR ) := A+ uL + A− uR , where A+ and A− are, respectively, the positive and negative parts of the Jacobian   0 1 A := . 1 0 Moreover, uL = ((H1 )L , (Q0 )L ) and uR = ((H1 )R , (Q0 )R ) are the left and right vector states. The finite-volume approximation then yields := unj − un+1 j ∆τ (H(ûj , ûj+1 ) − H(ûj−1 , ûj )) , ∆x (4.5) where the index j ∈ {1, . . . , J} denotes the mesh cell and n ∈ {1, . . . , N } the time step. The linear reconstruction is used in order to obtain the second order approximation. Slopes are limited by the so-called minmod limiter [17] to overcome over- and undershoots. Having thus obtained H1 (x, τ ) and Q0 (x, τ ), we can pass to the equation (3.15) with the initial condition (3.6) to obtain the additional term S. The right-hand side of (3.6) is approximated by the trapezoidal rule. In Figure 1 we compare the full model with the asymptotic model at different times t = 14 Tf , 12 Tf , 43 Tf , Tf (identified by different colors), both in the inelastic and elastic case. The long-time behaviour in the inelastic case is shown in Figure 2; there, the negative values for H are caused by the reflections at the boundary due to the Dirichlet type conditions for 1 We used MIKE URBAN 2009 with the module WH (Water Hammer). 13 H. Both Figures show a very good agreement of the full and asymptotic model when the same number of grid points is used. In Figure 3 we show the reflections of Q, which satisfies the asymptotic model, in the xt plane. An analogous simulation for the full model provides results indistinguishable from these ones, as we saw in Figure 2. In Figure 4 we compare the water-hammer case (pipe closed on the right, Qr = 0) with the open-pipe case (given pressure), see (4.3), by using the asymptotic model. The difference lies in the reflections of the pressure H and of the velocity Q at the right hand side (blue and yellow colors indicate positive and negative values, respectively). In Figure 5 we reduce the number of grid points in order to investigate efficiency of the multiscale model. We approximate now the first order hyperbolic system (3.8) by the LaxWendroff scheme in order to suppress the effects of different numerical discretization. The results obtained by the asymptotic and full model on a coarse grid are comparable, although some shift error can be noticed. Clearly, the modelling error is amplified by the larger discretization error, i.e. when a coarse grid is used. In the next experiments, the second order upwind finite volume method has been used for the approximation of the asymptotic model. The case of varying cross-sections (nozzles) is investigated in Figure 6. In that figure, above the increasing cross-section is shown, below the decreasing one, see (4.4). We again compare the full with the asymptotic model for H (figures on the left): the agreement is still good, although small phase shift became visible. Note the reflections caused by the variation of the cross-section. The figures on the right provide a contour plot of Q in the xt space. In all the examples provided up to now we fixed ζ = 0.4; such a value justified the weak damping assumption. As we emphasized, we found a very good agreement between the full and the asymptotic model (although ζ was not that small). We could not expect such good agreements a priori if we use the weak-damping asymptotic model for large values of ζ. As a test, however, the very large value ζ = 62 is considered in Figure 7. Surprisingly, the agreement is still quite good. Example 4.2 (Strong damping) In this experiment we consider the same initial-boundary value problems as above in the case of strong damping, i.e., ζ = 1/. We again compare the behavior of the full (2.6) and multiscale (3.24) models. In order to approximate the full model we use the Lax-Wendroff method described at the beginning of this section. Note that (3.24) is a semilinear hyperbolic system of balance laws. For the hyperbolic part we use the same method as for the case of weak damping, cf. (4.5). For the source term −α|Q0 |Q0 , the cell centered approximation at the n-th time step is used. Due to the explicit approximation in time, the CFL stability condition ∆τ /∆x ≤ CF L, CF L = 0.9 is used to limit the time step. The strong damping regime is characterized by large values of ζ. Note also that ζ = λL/(2Dr ), where λ is determined by the Colebrook formula, cf. Appendix A. Thus we either determine ζ or the absolute roughness of the pipe ar . In Figure 8 we take the value ζ = 62, already considered in Figure 7. In spite that such value is already very high, in Figures 8–10 we further increase the value of ζ to stress the effects of damping. In every case we compare simulations of the full model (dashed lines) with simulations of the asymptotic model (solid lines). In addition, for H we show results obtained by the commercial software package MIKE URBAN [11] 2 (dotted lines). 2 In MIKE URBAN the governing equations are based on the full model (2.14) but the nonlinear terms QHx 14 As in Example 4.1, in every figure we represent both H and Q, in the elastic (σ = 1) and in the inelastic case (σ = 0.84). Solutions are drawn at four different times t = 41 Tf , 21 Tf , 34 Tf , Tf . We see, as expected, the reduction of the wave speed in the inelastic case. Moreover, an increasing roughness damps the waves more and more. The results show very good agreement between the simulations of the full and the asymptotic model. The agreement with the MIKE URBAN package is also very good. Figure 11 focus on the long-time behavior Tf = 10  in the cases ζ = 62 and ζ = 380; the waves have already traveled few times back and forward. Again, the agreement between the asymptotic and the full model is very good. The effects of the roughness of the pipe on the damping can be seen very well. Notice that the values of Q at the left of the peak are negative in the case of the extremely large value ζ = 380. This effect can already be observed earlier in time (see Figure 10) but in the long-time behaviur it becomes more evident. This is due to the fact that the damping term “introduces” an additional negative derivative Qt which finally leads to a negative value of Q. In Figures 12 we reduce the number of grid points to illustrate the efficiency of the multiscale model. For short time-scales the solution obtained for the asymptotic model on a coarse grid agrees very well with the reference solution (full model on the fine grid). Analogously as in Figure 5 using coarse grids and long time-scales the effects of the phase shift errors are more visible. Finally, in Figure 13 we perform simulations in the case of varying cross-sections. The crosssections vary as in the weak damping case. Once more the reflections due to the variation of the cross-section can be observed. 5 Conclusions We presented two different asymptotic models for the description of sound waves in liquid flows in pipes. These multiscale models are motivated by two typical applications, namely, hydraulic plants and urban networks. The numerical simulations show a very good agreement, on one side, between the full models and the asymptotic models and, on the other, with the results provided by the commercial software package MIKE URBAN. Our approach opens the way to simulations of sound phenomena of fluid flow through curved (either S- or U -shaped, for instance) pipes or, more generally, in pipe networks. There, a good understanding of the phenomena on the single pipe level will be crucial. This will be a goal of our further research. Acknowledgements. We thank the DHI-WASY GmbH for providing a license for the software package MIKE URBAN. Andrea Corli thanks for the kind hospitality the Department of Mathematics of the University of Hamburg, where most of this research took place. He also thanks Stefano Alvisi and Alessandro Valiani for several enlightenments on flows through pipes. and QQx are neglected. The numerical method is a finite-difference scheme; it is fourth-order accurate, spacecompact and implicit [26]. The implicit finite-difference formulation is based on a non-staggered grid in time and space. 15 A Tables A table of reference quantities and typical values now follows. Concerning the water velocity in pipes and diameter of pipes, we distinguish two different framework: the case of an urban network (briefly, u.n.) and that of an hydraulic plant (h.p.). These essentially give the extremal ranges for the quantities under consideration. Quantity u ρ p x t D s z g ER φ σ Reference quantity ur ρr pr = K xr = L tr = L/ur Dr zr = L sin θr Typical value 1–10 m s−1 103 kg m−1 2 · 109 Pa 100 m 100–10 s 8 · 10−2 –1 m 4 · 10−3 –0.01 m 0–2 m 9.81 m s−2 2 · 1011 Pa (for iron) 10−10 Pa−1 0.84 Table 1: Scaling table with typical reference values. The quantities refer to those denoted ˜· in the system (2.1). The first number in the column “Typical value” refers to an urban network, the second to an hydraulic plant. The Darcy friction factor λ is computed as follows. Let µ be the dynamic viscosity for water, whose value at 20 ◦ C is about 10−3 N s m−3 . The Reynolds number for pipes is given by Re = ρr uµr Dr . Recall that a flow is said to be laminar when Re < 2300, transient when 2300 < Re < 4000 and turbulent when 4000 < Re. We are therefore largely in a turbulent regime, and this justify the quadratic term in the second equation of (2.1). Let ar the absolute −5 roughness of the pipe; this value ranges from about 10−3 m for rusty steel  or concrete  to 10 ar 2.51 + Re·λ we find m for new steel. Then, solving the Colebrook equation √1λ = −2 log10 3.7·D r the values reported in Table 2. References [1] AFT. AFT Impulse. http://www.aft.com/. [2] M. K. Banda, M. Herty, and A. Klar. Coupling conditions for gas networks governed by the isothermal Euler equations. Netw. Heterog. Media, 1(2):295–314, 2006. [3] M. K. Banda, M. Herty, and A. Klar. Gas flow in pipeline networks. Netw. Heterog. Media, 1(1):41– 56, 2006. [4] B. Brunone and M. Ferrante. Pressure waves as a tool for leak detection in closed conduits. Urban Water J., 2(2):145–155, 2004. [5] M. H. Chaudry. Applied Hydraulic Transients. Van Nostrand Reinhold, New York, second edition, 1987. 16 1 2 1   Re 1 (Fr)2 λ ζ u.n. h.p. 2 · 106 1.4 · 103 7.1 · 10−4 8 · 104 0 0.041–0.019 25–12 2 · 104 1.4 · 102 7.1 · 10−3 107 0–9.8 0.019–0.008 0.95–0.4 Table 2: Comparison of different values in the different regimes of an urban network, u.n., and hydraulic plant, h.p., L = 100. The first value for λ and ζ corresponds to rusty steel or concrete, the second one to new steel [6] R. M. Colombo and M. Garavello. A well posed Riemann problem for the p-system at a junction. Netw. Heterog. Media, 1(3):495–511, 2006. [7] R. M. Colombo and M. Garavello. On the p-system at a junction. In Control methods in PDEdynamical systems, volume 426 of Contemp. Math., pages 193–217. Amer. Math. Soc., Providence, RI, 2007. [8] R. M. Colombo and M. Garavello. On the Cauchy problem for the p-system at a junction. SIAM J. Math. Anal., 39(5):1456–1471, 2008. [9] S. Cox and E. Zuazua. The rate at which energy decays in a damped string. Comm. Partial Differential Equations, 19(1-2):213–243, 1994. [10] J. A. Cunge, F. M. Holly, and A. Verwey. Practical aspects of computational river hydraulics. Prentice Hall, Boston, 1980. [11] DHI. MOUSE. http://www.dhigroup.com/. [12] E. Feireisl. Strong decay for wave equations with nonlinear nonmonotone damping. Nonlinear Anal., 21(1):49–63, 1993. [13] L. Formaggia, D. Lamponi, and A. Quarteroni. One-dimensional models for blood flow in arteries. J. Engrg. Math., 47(3-4):251–276, 2003. [14] P. Garcia-Navarro and A. Priestley. A conservative and shape-preserving semi-Lagrangian method for the solution of the shallow water equations. Internat. J. Numer. Methods Fluids, 18(3):273–294, 1994. [15] H. Holden and N. H. Risebro. Riemann problems with a kink. SIAM J. Math. Anal., 30(3):497–515, 1999. [16] B. E. Larock, R. W. Jeppson, and G. Z. Watters. Hydraulics of pipeline systems. CRC Press, Boca Raton, FL, 1999. [17] R. J. LeVeque. Finite volume methods for hyperbolic problems. Mathematics. Cambridge University Press, Cambridge, 2002. Cambridge Texts in Applied [18] T. T. Li. Global classical solutions for quasilinear hyperbolic systems, volume 32 of RAM: Research in Applied Mathematics. Masson, Paris, 1994. [19] T. P. Liu. Transonic gas flow in a duct of varying area. Arch. Rational Mech. Anal., 80(1):1–18, 1982. [20] M. Luskin. On the existence of global smooth solutions for a model equation for fluid flow in a pipe. J. Math. Anal. Appl., 84(2):614–630, 1981. 17 [21] M. Luskin and B. Temple. The existence of a global weak solution to the nonlinear waterhammer problem. Comm. Pure Appl. Math., 35(5):697–735, 1982. [22] R. Pan and K. Zhao. Initial boundary value problem for compressible Euler equations with damping. Indiana Univ. Math. J., 57(5):2257–2282, 2008. [23] M. A. Rammaha and T. A. Strei. Global existence and nonexistence for nonlinear wave equations with damping and source terms. Trans. Amer. Math. Soc., 354(9):3621–3637, 2002. [24] I. A. Sibetheros and E. R. Holley. Spline interpolations for water hammer analysis. J. Hydraul. Eng., 117:1332–1351, 1991. [25] R. Szymkiewicz and M. Mitosek. Numerical aspects of improvement of the unsteady pipe flow equations. Internat. J. Numer. Methods Fluids, 55(11):1039–1058, 2007. [26] A. Verwey and J. H. Yu. A space-compact high-order implicit scheme for water hammer simulations. Proceedings of XXVth IAHR Congress, Tokyo, 1993. [27] F. M. White. Fluid Mechanics. McGraw-Hill, Boston, fifth edition, 2003. [28] E. B. Wylie and V. L. Streeter. Fluid transients. McGraw-Hill, New York, 1978. [29] W. Zielke. Elektronische Berechnung vor Rohr- und Gerinneströmungen. Erich Schmidt Verlag, Berlin, 1974. 18 −3 8 x 10 1.2 1 6 0.8 0.6 H Q 4 0.4 2 0.2 0 0 0 0.2 0.4 0.6 0.8 −0.2 1 0 0.2 0.4 x 0.6 0.8 1 0.8 1 x (a) piezometric head H, σ = 1 (b) velocity Q, σ = 1 −3 8 x 10 1.2 1 6 0.8 0.6 H Q 4 0.4 2 0.2 0 0 0 0.2 0.4 0.6 0.8 −0.2 1 0 0.2 0.4 x 0.6 x (c) piezometric head H, σ = 0.84 (d) velocity Q, σ = 0.84 Figure 1: Smooth pipe, hydraulic plant: ζ = 0.4, L = 100m, D = 1m, Tf =  = 7.1 · 10−3 . Full model (dashed line) vs. asymptotic model for weak damping (solid line) at four different times. −3 8 x 10 1.2 6 1 4 0.8 2 Q H 0.6 0 0.4 −2 0.2 −4 0 −6 −8 0 0.2 0.4 0.6 0.8 −0.2 1 x 0 0.2 0.4 0.6 0.8 1 x (a) piezometric head H, σ = 1 (b) velocity Q, σ = 1 Figure 2: As in Figure 1 but for long time Tf = 10 = 7.1 · 10−2 . Full model (dashed line) vs. asymptotic model for weak damping (solid line). 19 (a) velocity Q, σ = 1 Figure 3: As in Figure 1 but for long time Tf = 10 = 7.1 · 10−2 . Velocity Q from the asymptotic model for the weak damping. (a) piezometric head H, water hammer (b) velocity Q, water hammer (c) piezometric head H, open pipe (d) velocity Q, open pipe Figure 4: As in Figure 1, σ = 1, but for long time Tf = 10 = 7.1 · 10−2 . Closed pipe (water hammer) vs. open pipe case for the asymptotic model for the weak damping. 20 −3 8 x 10 1.2 1 6 0.8 4 H Q 0.6 0.4 2 0.2 0 0 −2 0 0.2 0.4 0.6 0.8 −0.2 1 0 0.2 0.4 x 0.6 0.8 1 x (a) piezometric head H, Tf =  = 7.1 · 10−3 (b) velocity Q, Tf =  = 7.1 · 10−3 −3 8 x 10 1.2 6 1 4 0.8 2 Q H 0.6 0 0.4 −2 0.2 −4 0 −6 −8 0 0.2 0.4 0.6 0.8 −0.2 1 x 0 0.2 0.4 0.6 0.8 1 x (c) piezometric head H, Tf = 10 = 7.1 · 10−2 (d) velocity Q, Tf = 10 = 7.1 · 10−2 Figure 5: As in Figure 1, σ = 1: full model with N = 1000 grid cells (solid line), full model with N = 100 (dotted line), asymptotic model for weak damping with N = 100 (dashed line). 21 −3 8 x 10 6 H 4 2 0 −2 0 0.2 0.4 0.6 0.8 1 x (a) piezometric head H (b) velocity Q, asymptotic model −3 x 10 8 H 6 4 2 0 0 0.2 0.4 0.6 0.8 1 x (c) piezometric head H (d) velocity Q, asymptotic model Figure 6: As in Figure 1, σ = 1. Above: increasing cross-section A1 , below: decreasing crosssection A2 . Full model (dashed line) vs. asymptotic model for weak damping (solid line). −4 8 x 10 1.2 6 1 4 0.8 2 Q H 0.6 0 0.4 −2 0.2 −4 0 −6 −8 0 0.2 0.4 0.6 0.8 −0.2 1 x 0 0.2 0.4 0.6 0.8 1 x (a) piezometric head H (b) velocity Q Figure 7: Long-time behaviour, urban network, ζ = 62 or ar ≈ 0.0028mm, L = 1000m, D = 0.08 m, σ = 1, Tf = 10 = 7.1 · 10−3 . Full model (dashed line) vs. asymptotic model for weak damping (solid line). 22 −4 8 x 10 1.2 1 6 0.8 0.6 H Q 4 0.4 2 0.2 0 0 0 0.2 0.4 0.6 0.8 −0.2 1 0 0.2 0.4 x 0.6 0.8 1 0.8 1 x (a) piezometric head H, σ = 1 (b) velocity Q, σ = 1 −4 8 x 10 1.2 1 6 0.8 0.6 H Q 4 0.4 2 0.2 0 0 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 x x (c) piezometric head H, σ = 0.84 (d) velocity Q, σ = 0.84 Figure 8: Urban network, ζ = 62 or ar ≈ 0.0028mm, L = 1000m, D = 0.08m, Tf =  = 7.1 · 10−4 . Full model (dashed line), asymptotic model for strong damping (solid line), MIKE URBAN (dotted line, only for H). 23 −4 8 x 10 1.2 1 6 0.8 0.6 H Q 4 0.4 2 0.2 0 0 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.8 1 x x (a) piezometric head H, σ = 1 (b) velocity Q, σ = 1 −4 8 x 10 1.2 1 6 0.8 0.6 H Q 4 0.4 2 0.2 0 0 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 x x (c) piezometric head H, σ = 0.84 (d) velocity Q, σ = 0.84 Figure 9: As in Figure 8 but for a very rough pipe: ζ = 170 or ar ≈ 0.28mm. 24 −4 8 x 10 1.2 1 6 0.8 0.6 H Q 4 0.4 2 0.2 0 0 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.8 1 x x (a) piezometric head H, σ = 1 (b) velocity Q, σ = 1 −4 8 x 10 1.2 1 6 0.8 0.6 H Q 4 0.4 2 0.2 0 0 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 x x (c) piezometric head H, σ = 0.84 (d) velocity Q, σ = 0.84 Figure 10: As in Figure 8 but for an extremely rough pipe: ζ = 380 or ar ≈ 2.8mm. 25 −4 8 x 10 1.2 6 1 4 0.8 2 Q H 0.6 0 0.4 −2 0.2 −4 0 −6 −8 0 0.2 0.4 0.6 0.8 −0.2 1 0 0.2 0.4 x 0.6 0.8 1 0.8 1 x (a) piezometric head H, ζ = 62 (b) velocity Q, ζ = 62 −4 8 x 10 0.8 0.7 6 0.6 4 0.5 2 Q H 0.4 0 0.3 −2 0.2 −4 0.1 −6 −8 0 0 0.2 0.4 0.6 0.8 −0.1 1 x 0 0.2 0.4 0.6 x (c) piezometric head H, ζ = 380 (d) velocity Q, ζ = 380 Figure 11: As in Figure 8, long-time behaviour, σ = 1, Tf = 10 = 7.1·10−3 . Full model (dashed line) vs. asymptotic model for strong damping (solid line). 26 −4 8 x 10 1.2 1 6 0.8 0.6 H Q 4 0.4 2 0.2 0 0 0 0.2 0.4 0.6 0.8 −0.2 1 0 0.2 0.4 x 0.6 0.8 1 x (a) piezometric head H, Tf =  = 7.1 · 10−4 (b) velocity Q, Tf =  = 7.1 · 10−4 −4 8 x 10 0.8 0.7 6 0.6 4 0.5 Q H 0.4 2 0.3 0 0.2 0.1 −2 0 −4 0 0.2 0.4 0.6 0.8 −0.1 1 x 0 0.2 0.4 0.6 0.8 1 x (c) piezometric head H, Tf = 10 = 7.1 · 10−3 (d) velocity Q, Tf = 10 = 7.1 · 10−3 Figure 12: As in Figure 10. Full model with N = 1000 grid cells (solid line), full model with N = 100 (dotted line), asymptotic model for strong damping with N = 100 (dashed line). 27 −4 8 x 10 6 H 4 2 0 0 0.2 0.4 0.6 0.8 1 x (a) piezometric head H (b) velocity Q, asymptotic model −4 10 x 10 8 H 6 4 2 0 0 0.2 0.4 0.6 0.8 1 x (c) piezometric head H (d) velocity Q, asymptotic model Figure 13: As in Figure 10, Tf =  = 7.1 · 10−4 . Above: increasing cross-section A1 , below: decreasing cross-section A2 . Full model (dashed line) vs. asymptotic model for strong damping (solid line). 28