1. Introduction
The water hammer in pipes is a well-known phenomenon that has been studied by numerous authors for over a century (Michaud Reference Michaud1878; Allievi Reference Allievi1903). Water hammer occurs when pipe flow is abruptly altered due to the opening or closing of valves. In such cases, the boundary condition at the valve couples velocity and hydraulic head through a nonlinear relation based on Bernoulli’s principle. Owing to this complexity, to the best of our knowledge, no analytical solutions are currently available. Depending on whether the valve closing or opening time is shorter or longer than the critical time, the closure is classified as sudden or slow. While for sudden closure the maximum pressure surge reaches the well-known Joukowsky maximum value (Joukowsky Reference Joukowsky1898), for more realistic slow manoeuvres the maximum surge is smaller and needs to be determined numerically.
In his pioneering work, Allievi developed a method for determining the exact time-dependent behaviour of the hydraulic head at the valve for both sudden and slow closures. This technique remains widely used in the design of piping systems and is taught in standard hydraulics courses. Another well-established approach is the method of characteristics for solving parabolic systems (Mei Reference Mei1997). However, these techniques rely on the numerical solution of equations that couple hydraulic head and velocity, and therefore do not yield practical, explicit expressions.
Due to the inherent complexity of the governing equations, some authors have focused on solving the water hammer problem by prescribing the velocity at the valve – a reasonable assumption only in the case of sudden closures (Di Nucci et al. Reference Di Nucci, Urbanowicz, Michele, Celli, Pasquali and Di Risio2025). In such cases, a linearised, space-dependent solution can be readily obtained via modal expansions. Further extensions include the analysis of the effects such as frictional damping (Urbanowicz et al. Reference Urbanowicz, Bergant, Stosiak, Karpenko and Bogdevicius2023b , Reference Urbanowicz, Jing, Bergant, Stosiak and Lubeckic ), accelerating pipe flow (Urbanowicz et al. Reference Urbanowicz, Bergant, Stosiak, Deptula and Karpenko2023a ), pipe viscoelasticity (Bergant et al. Reference Bergant, Tijsseling, Vitkowsky, Covas, Simpson and Lambert2008a ; Mei et al. Reference Mei, Li, Michele, Sammarco and McBeth2021), leakage (Bergant et al. Reference Bergant, Tijsseling, Vitkowsky, Covas, Simpson and Lambert2008b ), column separation (Bergant, Simpson & Tijsseling Reference Bergant, Simpson and Tijsseling2006), turbulence (Di Nucci, Michele & Di Risio Reference Di Nucci, Michele and Di Risio2024), polymer pipes (Yao, Kember & Hansen Reference Yao, Kember and Hansen2016) and two-dimensional geometries (Ghidaui et al. Reference Ghidaui, Zhao, McInnis and Axworthy2005).
In this paper, we address the water hammer problem under slow valve manoeuvres by applying a perturbation scheme, where the small parameter is the ratio of the critical time to the valve closing or opening time. This enables us to decompose the nonlinear problem into a sequence of linearised boundary-value problems, which we solve using the Laplace transform. We extend the expansion to second order and obtain very good agreement between the exact numerical solution derived from Allievi’s coupled equations and the analytical results developed herein. Furthermore, we derive explicit space- and time-dependent expressions for the hydraulic head at each order. These expressions enable the estimation of maximum and minimum hydraulic head values as a function of Allievi’s number and could likely be incorporated into design standards for engineering of pipe systems.
2. Governing equations
Referring to figure 1, we consider a straight cylindrical pipe of length
$L$
and constant cross-sectional area
$\varOmega$
, connected to a reservoir of constant fluid depth. Let
$t$
denote time and
$x$
be the spatial coordinate oriented along the axis of the pipe, such that
$x=0$
corresponds to the pipe inlet, while
$x=L$
marks the location of a nozzle. The nozzle terminates with a cross-sectional area
$\omega (t)$
, where
$\omega (t) \ll \varOmega$
, and its length is assumed negligible relative to
$L$
. A valve located at the end of the nozzle regulates
$\omega (t)$
, thus controlling the cross-sectional average velocity
$U(x,t)$
within the pipe, which is defined as positive in the rightward direction. The horizontal reference plane
$z=0$
is set at the valve elevation, while the free surface elevation within the reservoir is fixed at
$z=H$
. We assume that the pipe diameter is small compared with its length
$L$
, that damping effects due to fluid viscosity, turbulence and pipe-wall viscoelasticity are negligible, and that
$ U \ll c$
, where
$ c$
is the constant wave speed in the pipe. The latter assumption is equivalent to assuming small Mach numbers,
$ {Ma} \ll 1$
(White Reference White1991). Under these conditions, the continuity and momentum equations can be expressed as
\begin{equation} \frac {\partial U}{\partial t} + g\frac {\partial h}{\partial x}= 0,\quad \frac {\partial U}{\partial x}c^2 + g\frac {\partial h}{\partial t}= 0, \end{equation}
where the hydraulic head
$h(x,t)$
is defined as
$h(x,t) = z(x) + p(x,t)/(\rho g)$
, with
$g$
denoting the acceleration due to gravity,
$z(x)$
the elevation of the pipe axis,
$\rho$
the fluid density and
$p(x,t)$
the pressure.
Sketch representing the reservoir–pipe system.
By combining (2.1), we obtain the classical wave equations governing
$U$
and
$h$
\begin{equation} c^2 \frac {\partial ^2 U}{\partial x^2} - \frac {\partial ^2 U}{\partial t^2} = 0,\quad c^2 \frac {\partial ^2 h}{\partial x^2} - \frac {\partial ^2 h}{\partial t^2} = 0. \end{equation}
At the inlet (
$x=0$
), a constant hydraulic head
$h = H$
is imposed. Therefore
\begin{align} h(0,t) = H \quad \rightarrow \quad \left.\frac {\partial h(x,t)}{\partial t}\right |_{x=0} = 0 \quad \rightarrow \quad \left.\frac {\partial U(x,t)}{\partial x}\right |_{x=0} = 0. \end{align}
In this work, we address both valve opening and closure problems, and thus the initial conditions differ based on the scenario considered. For the valve closure problem, we assume a steady-state flow condition prior to closure, leading to the initial conditions
\begin{equation} U(x,0) = U_0 = \frac {\omega _0}{\varOmega } \sqrt {2\textit{gH}}, \quad h(x,0) = H,\quad \left.\frac {\partial U(x,t)}{\partial t}\right |_{t=0} = \left.\frac {\partial h(x,t)}{\partial t}\right |_{t=0} = 0, \end{equation}
where
$U_0$
is constant and
$\omega _0$
is the value of the cross-sectional area when the valve is fully open. For the valve opening problem, we instead assume
\begin{align} U(x,0) = 0, \quad h(x,0) = H,\quad \left.\frac {\partial U(x,t)}{\partial t}\right |_{t=0} = \left.\frac {\partial h(x,t)}{\partial t}\right |_{t=0} = 0, \end{align}
i.e. a zero initial velocity and a hydraulic head equal to the free surface elevation in the reservoir. We now turn to the boundary condition at
$x=L$
. It is worth noting that most classical models of water hammer phenomena assume valve closures or openings, where the outlet velocity is directly prescribed. This modelling approach permits a complete analytical solution, as presented in Di Nucci et al. (Reference Di Nucci, Urbanowicz, Michele, Celli, Pasquali and Di Risio2025), and enables the recovery of the well-known Joukowsky value as the maximum attainable pressure in the pipe (Joukowsky Reference Joukowsky1898). However, this assumption is often unrealistic for two main reasons: (i) in practice, valve operations occur over a finite time interval; (ii) the outlet velocity is nonlinearly related to the local hydraulic head
$h(L,t)$
according to Bernoulli’s principle. To this end, we assume that the cross-sectional area of the nozzle,
$\omega (t)$
, evolves over time according to
\begin{align} \omega (t) &= \omega _0\eta (t), \quad t \in (0,T], \end{align}
\begin{align} \omega (t) &= \omega _0\eta (T), \quad t \gt T, \end{align}
where
$\eta$
represents the time sequence of the manoeuvre and
$T$
is the duration of valve closure or opening manoeuvre, i.e. the opening ratio, and is defined by
\begin{equation} \eta (t)= \alpha +\beta f(t), \quad t \in (0,T], \end{equation}
where
$\alpha =1$
,
$\beta =-1$
for the valve closure problem and
$\alpha =0$
,
$\beta =1$
for the valve opening problem, whereas
$f\in [0;1]$
is a positive function of time.
By applying Bernoulli’s principle, the boundary condition at the outlet
$x=L$
is given by
\begin{equation} U(L,t)\varOmega = \omega (t)\sqrt {2gh(L,t)},\,\,\rightarrow \,\,\frac {U(L,t)}{U_0} = \eta (t)\sqrt {\frac {h(L,t)}{H}}. \end{equation}
Equation (2.9) is nonlinear and couples the variables
$U$
and
$h$
. The solution to (2.2), subject to the initial and boundary conditions (2.3)–(2.9), can be obtained through the method of characteristics (Mei Reference Mei1997). The pioneering numerical Allievi method enables the determination of velocity and pressure at the valve
$x=L$
(Allievi Reference Allievi1903). However, these approaches typically require the implementation of a numerical scheme to resolve
$U$
and
$h$
at each time step, thus losing physical insights into the parametric dependence. Furthermore, to the best of the authors’ knowledge, a fully analytical solution to the aforementioned problem has not yet been reported in the literature.
To find the analytical solution, we first introduce the following non-dimensional quantities:
\begin{equation} t^{\prime} = \frac {2t}{\theta }=\frac {ct}{L},\quad x^{\prime} = \frac {x}{L},\quad U^{\prime} = \frac {U}{U_0},\quad h^{\prime} = \frac {h}{H}, \end{equation}
where
$\theta = 2L/c$
denotes the time required for a pressure wave to travel along the pipe and return, commonly referred to as the critical time (‘durata di fase’ in Allievi Reference Allievi1903). Assuming that
$f(t)$
varies slowly with time, we introduce the small parameter
\begin{equation} \epsilon = \frac {{St}}{2}= \frac {\theta }{T} \ll 1, \end{equation}
where
${St}$
denotes the Strouhal number (White Reference White1991). The assumption (2.11) implies a slow valve manoeuvre, a situation of practical interest in many engineering applications. Substituting (2.10) into (2.1) we get continuity and momentum equations in non-dimensional form
\begin{equation} 2{Al}\frac {\partial U^{\prime}}{\partial t^{\prime}} + \frac {\partial h^{\prime}}{\partial x^{\prime}}= 0,\quad 2{Al}\frac {\partial U^{\prime}}{\partial x^{\prime}} + \frac {\partial h^{\prime}}{\partial t^{\prime}}= 0, \end{equation}
where
$ {Al} = c U_0 / (2\textit{gH})$
is the Allievi number (Allievi Reference Allievi1903). To provide a clearer physical interpretation of the Allievi number, we rewrite
$ {Al}$
as
\begin{equation} {Al} = \frac {\rho c U_0^2}{2 \rho \textit{gH} U_0} = \frac {F_e}{2F_s}, \end{equation}
where
$ F_e = \rho U_0^2 c$
and
$ F_s = \rho \textit{gH} U_0$
denote the energy flux densities of the pressure waves (Landau & Lifshitz Reference Landau and Lifshitz1989, Chapter 8) and of the base flow, respectively. Large values of
$ {Al}$
therefore are expected to give rise to substantial hydraulic head oscillations, a result that will be confirmed mathematically later in this section.
Equations (2.12) allow us to write the following boundary-value problem for
$h^{\prime}$
:
\begin{align} \frac {\partial ^2 h^{\prime}}{\partial {x^{\prime}}^2} - \frac {\partial ^2 h^{\prime}}{\partial {t^{\prime}}^2} &= 0, \qquad x^{\prime}\in (0,1), t^{\prime}\gt 0, \end{align}
\begin{align} h^{\prime} &= 1, \qquad x^{\prime}\in (0,1), t^{\prime}=0, \end{align}
\begin{align} h^{\prime}&= 1, \qquad x^{\prime}=0, t^{\prime}\gt 0, \end{align}
\begin{align} \frac {\partial \eta }{\partial t^{\prime}}\sqrt {h^{\prime}}+\frac {\eta }{2\sqrt {h^{\prime}}}\frac {\partial h^{\prime}}{\partial t^{\prime}}+\frac {1}{2{Al}}\frac {\partial h^{\prime}}{\partial x^{\prime}} &=0, \qquad x^{\prime}=1, t^{\prime}\gt 0, \end{align}
where the boundary condition (2.17) is obtained by combining the first of (2.12) and the time derivative of (2.9). Since
$f$
is a slowly varying function of time we can write
\begin{equation} \eta = \alpha +\beta f^{\prime}(t^{\prime})\epsilon , \quad t \in (0,2/\epsilon ], \end{equation}
where
$f=f^{\prime}\epsilon$
. We are now in a condition to adopt a weakly nonlinear approach based on a perturbation expansion in
$\epsilon$
. Specifically, we expand
$h^{\prime}$
as
\begin{equation} h^{\prime} = h^{\prime}_{0} + \epsilon h^{\prime}_{1} + \epsilon ^2 h^{\prime}_{2} + O(\epsilon ^3). \end{equation}
Substituting the expansion (2.19) into the governing equation and boundary conditions (2.14)–(2.17), and Taylor expanding (2.17) about
$\epsilon =0$
, we obtain a sequence of linearised boundary-value problems up to order
$O (\epsilon ^2 )$
.
We now focus on the analytical solution of the hydraulic head
$h^{\prime}$
at each order.
2.1. Leading-order problem
The boundary-value problem at the leading-order is given by
\begin{align} \frac {\partial ^2 h_0^{\prime}}{\partial {x^{\prime}}^2} - \frac {\partial ^2 h_0^{\prime}}{\partial {t^{\prime}}^2} &= 0, \qquad x^{\prime}\in (0,1), t^{\prime}\gt 0, \end{align}
\begin{align} h_0^{\prime} &= 1, \qquad x^{\prime}\in (0,1), t^{\prime}=0, \end{align}
\begin{align} h_0^{\prime}&= 1, \qquad x^{\prime}=0, t^{\prime}\gt 0, \end{align}
\begin{align} \frac {\alpha }{2\sqrt {h_0^{\prime}}}\frac {\partial h_0^{\prime}}{\partial t^{\prime}}+\frac {1}{2{Al}}\frac {\partial h_0^{\prime}}{\partial x^{\prime}} &=0, \qquad x^{\prime}=1, t^{\prime}\gt 0. \end{align}
Since there is no forcing at this order, we obtain the straightforward solution
$h^{\prime}_0 = 1$
, which corresponds to the steady-state hydraulic head in the absence of wave-like perturbations.
2.2. First-order problem
At the first order we obtain
\begin{align} \frac {\partial ^2 h_1^{\prime}}{\partial {x^{\prime}}^2} - \frac {\partial ^2 h_1^{\prime}}{\partial {t^{\prime}}^2} &= 0, \qquad\qquad\qquad x^{\prime}\in (0,1), t^{\prime}\gt 0, \end{align}
\begin{align} h_1^{\prime} &= 0, \qquad\qquad\qquad x^{\prime}\in (0,1), t^{\prime}=0, \end{align}
\begin{align} h_1^{\prime}&= 0, \qquad\qquad\qquad x^{\prime}=0, t^{\prime}\gt 0, \end{align}
\begin{align} \alpha \frac {\partial h_1^{\prime}}{\partial t^{\prime}}+\frac {1}{{Al}}\frac {\partial h_1^{\prime}}{\partial x^{\prime}} &=-2\beta \frac {\partial f^{\prime}}{\partial t^{\prime}}, \quad\qquad x^{\prime}=1, t^{\prime}\gt 0. \end{align}
The boundary-value problem (2.24)–(2.27) is forced by the known time-dependent term
$\partial f^{\prime}/\partial t^{\prime}$
, which is the operator imposed speed at which the cross-sectional area changes at the valve. We expect the first-order solution to be significantly different between the valve opening and closing cases because the term
$\alpha$
is null or equal to unity, respectively. The time-dependent solution of (2.24)–(2.27) can be found by first applying the Laplace transform and its inverse (Mei Reference Mei1997)
\begin{equation} \overline {h_1^{\prime}}(x^{\prime},s)=\int _{0}^\infty \mathrm{e}^{-st^{\prime}}h_1^{\prime}\,\mathrm{d}t^{\prime},\quad h_1^{\prime}(x^{\prime},t^{\prime})=\frac {1}{2\pi \mathrm{i}}\int _{a-\mathrm{i}\infty }^{a+\mathrm{i}\infty }\mathrm{e}^{st^{\prime}}\overline {h_1^{\prime}}\,\mathrm{d}s, \end{equation}
where
$a$
is a positive constant,
$s$
is a complex variable and the path of integration is a vertical line parallel and to the right of the imaginary axis in the complex
$s$
plane. Applying (2.28) to the boundary-value problem (2.24)–(2.27), we obtain
\begin{align} \frac {\partial ^2 \overline {h_1^{\prime}}}{\partial {x^{\prime}}^2} - s^2\overline {h_1^{\prime}} &= 0, \qquad\qquad\quad x^{\prime}\in (0,1), \end{align}
\begin{align} \overline {h_1^{\prime}}&= 0, \qquad\qquad\quad x^{\prime}=0, \end{align}
\begin{align} \alpha s\overline {h_1^{\prime}}+\frac {1}{{Al}}\frac {\partial \overline {h_1^{\prime}}}{\partial x^{\prime}} &=-2\beta s\overline {f^{\prime}}, \qquad x^{\prime}=1, \end{align}
where we applied Laplace transform properties to time derivatives and resorted to the initial condition
$h_1^{\prime}(x,0)=0, f^{\prime}(0)=0$
. Solution of (2.29)–(2.31) yields
\begin{equation} h_1^{\prime}=-\beta \frac {{Al}}{\pi \mathrm{i}}\int _{a-\mathrm{i}\infty }^{a+\mathrm{i}\infty }\frac {\overline {f^{\prime}}\mathrm{e}^{st^{\prime}}\sinh (sx^{\prime})}{\cosh (s)+\alpha {Al}\sinh (s)}\,\mathrm{d}s. \end{equation}
The complex integral (2.32) can be solved once
$\overline {f^{\prime}}$
is determined. We consider a linear variation in the valve complete closure/opening manoeuvre and obtain
\begin{equation} f=\frac {t}{T},\,\,\rightarrow \,\,f^{\prime}=\frac {t^{\prime}}{2},\,\,\rightarrow \,\,\overline {f^{\prime}}=\frac {1}{2}\int _{0}^\infty \mathrm{e}^{-st^{\prime}}t^{\prime}\,\mathrm{d}t^{\prime}=\frac {1}{2s^2}. \end{equation}
Substituting the latter into (2.32) yields
\begin{equation} h_1^{\prime}=-\beta \frac {{Al}}{2\pi \mathrm{i}}\int _{a-\mathrm{i}\infty }^{a+\mathrm{i}\infty }\frac {\mathrm{e}^{st^{\prime}}\sinh (sx^{\prime})}{s^2\left [\cosh (s)+\alpha {Al}\sinh (s)\right ]}\,\mathrm{d}s. \end{equation}
Let us first assume
${Al} \neq 1$
. We find that there is a pole at the origin,
$s = 0$
, and a series of poles at
\begin{equation} s_n =\text{Re}\{s_n\}+\mathrm{i}\,\text{Im}\{s_n\}= \frac {1}{2} \log \left (\left | \frac {1 - \alpha {Al}}{1 + \alpha {Al}} \right |\right ) + \mathrm{i}\, \pi \times \left \{ \begin{aligned} & \frac {1}{2} + n, \quad \text{if }\, \alpha {Al} \lt 1 ,\\ & n, \quad\qquad \text{if }\, \alpha {Al} \gt 1, \end{aligned} \right. \quad n \in \mathbb{Z}, \end{equation}
which are the complex roots of
$\cosh (s) + \alpha {Al} \sinh (s) = 0$
. In the case of valve closure (
$\alpha = 1$
), we have
$\text{Re}\{s_n\} \lt 0$
, indicating that each pole corresponds to an oscillatory and damped mode; the
$s_n$
lie to the left of the imaginary axis. For the valve opening problem (
$\alpha = 0$
), we obtain undamped oscillatory modes since
$\text{Re}\{s_n\} = 0$
. Moreover, we note that
$\text{Im}\{s_n\}$
depends on the product
$\alpha {Al}$
, which has important implications for the behaviour of the solution, as will be shown later. To solve (2.34), we introduce a large semicircle to the left of the integration path and apply both Jordan’s lemma and Cauchy’s residue theorem to obtain
\begin{equation} h_1^{\prime} = -\beta {Al} \left [ x^{\prime} + \sum _{n=-\infty }^{\infty } \frac {\mathrm{e}^{s_n t^{\prime}} \sinh (s_n x^{\prime})}{s_n^2 \left [ \sinh (s_n) + \alpha {Al} \cosh (s_n) \right ]} \right ]. \end{equation}
We note that this leading-order solution comprises two components: a series of oscillatory modes with sinusoidal shape, and a time-independent term representing the contribution from
$s = 0$
, obtained by expanding the integrand of (2.34) as a Laurent series. In addition, from (2.35) it follows that, as
$ {Al} \rightarrow 0$
, the damping of the oscillatory modes decreases. However, in this limit the hydraulic head oscillations are also reduced, since the solution (2.36) scales with
$ {Al}$
.
In the case of valve closure (
$\alpha = 1$
) with
${Al} = 1$
, a different approach is required because the poles
$s_n$
shift towards
$\text{Re}\{s\} \to -\infty$
(see (2.35)). Therefore, we rewrite (2.34) as
\begin{equation} h_1^{\prime} = \frac {1}{4\pi \mathrm{i}} \left ( \int _{a - \mathrm{i} \infty }^{a + \mathrm{i} \infty } \frac {\mathrm{e}^{s(t^{\prime} + x^{\prime} - 1)}}{s^2} \, \mathrm{d}s - \int _{a - \mathrm{i} \infty }^{a + \mathrm{i} \infty } \frac {\mathrm{e}^{s(t^{\prime} - x^{\prime} - 1)}}{s^2} \, \mathrm{d}s \right ). \end{equation}
Since there is a single pole at
$s = 0$
, we apply Jordan’s lemma and obtain the following solution:
\begin{equation} h_1^{\prime} = \frac {1}{2} \left [ (t^{\prime} + x^{\prime} - 1)\, \text{H}(t^{\prime} + x^{\prime} - 1) - (t^{\prime} - x^{\prime} - 1)\, \text{H}(t^{\prime} - x^{\prime} - 1) \right ], \end{equation}
where
$\text{H}$
denotes the Heaviside step function. Interestingly, when Allievi’s number is exactly equal to unity, oscillatory modes are absent, and the hydraulic head increases linearly with time and then remains constant for
$t^{\prime} \gt 2$
.
The first-order solution, given by (2.36) and (2.38), is now complete. To better understand its behaviour, we examine the values of
$h_1^{\prime}$
at the valve location,
$x^{\prime} = 1$
, and at integer times
$t^{\prime} = 1, 2, 3, \ldots$
. This analysis also enables us to derive an explicit expression for (2.36) via series expansions of the hyperbolic functions (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007).
For the valve closure problem, we obtain:
\begin{equation} h_1^{\prime}(1,t^{\prime}) = \left \{ \begin{aligned} &{Al} \left [ 1 - {Al} ({Al} - 1)^{\frac {t^{\prime} - 1}{2}} ({Al} + 1)^{-\frac {t^{\prime} + 1}{2}} \right ], \qquad t^{\prime} = 1,3,5,\ldots. \\ &{Al} \left [ 1 - \left ( \frac {{Al} - 1}{{Al} + 1} \right )^{\frac {t^{\prime}}{2}} \right ], \qquad\qquad\qquad\qquad t^{\prime} = 2,4,6,\ldots. \end{aligned} \right. \end{equation}
For
${Al} \lt 1$
, expression (2.39) indicates that the maximum hydraulic head occurs at
$t^{\prime} = 2$
(i.e.
$t=\theta$
), when the first pressure wave reflects back to the valve. In this case, we find
$h_1^{\prime}(1,2) = 2{Al} / ({Al} + 1)$
. Conversely, for
${Al} \gt 1$
, the expression shows that the maximum value increases with each subsequent wave reflection at
$x^{\prime} = 1$
, asymptotically approaching the limiting value
$h_1^{\prime} = {Al}$
. When
${Al} = 1$
, the hydraulic head reaches its maximum value
$h_1^{\prime} = {Al}$
at
$t^{\prime} = 2$
and remains constant thereafter, consistent with the behaviour described above.
A similar procedure, applied to the valve opening problem, yields
\begin{equation} h_1^{\prime}(1,t^{\prime}) = \left \{ \begin{aligned} &{Al}, \qquad\qquad\qquad\quad\qquad t^{\prime} = 1,3,5,\ldots. \\ & -{Al} \left [ 1 - (-1)^{\frac {t^{\prime}}{2}} \right ], \qquad t^{\prime} = 2,4,6,\ldots.\end{aligned} \right. \end{equation}
Expression (2.40) shows that the hydraulic head reaches its minimum value,
$h_1^{\prime} = -2{Al}$
, at
$t^{\prime} = 2,6,\ldots$
, while at
$t^{\prime} = 4,8,\ldots$
it returns to
$h_1^{\prime} = 0$
, which corresponds to a local maximum.
2.3. Second-order problem
The second-order boundary-value problem reads
\begin{align} \frac {\partial ^2 h_2^{\prime}}{\partial {x^{\prime}}^2} - \frac {\partial ^2 h_2^{\prime}}{\partial {t^{\prime}}^2} &= 0, \!\!\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad x^{\prime}\in (0,1), t^{\prime}\gt 0, \end{align}
\begin{align} h_2^{\prime} &= 0, \!\!\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad x^{\prime}\in (0,1), t^{\prime}=0, \end{align}
\begin{align} h_2^{\prime}&= 0, \!\!\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad x^{\prime}=0, t^{\prime}\gt 0, \end{align}
\begin{align} \alpha \frac {\partial h_2^{\prime}}{\partial t^{\prime}}+\frac {1}{{Al}}\frac {\partial h_2^{\prime}}{\partial x^{\prime}} &=-\beta \left (h_1^{\prime}\frac {\partial f^{\prime}}{\partial t^{\prime}}+f^{\prime}\frac {\partial h_1^{\prime}}{\partial t^{\prime}}\right )+\frac {\alpha }{2}h_1^{\prime}\frac {\partial h_1^{\prime}}{\partial t^{\prime}}, \qquad x^{\prime}=1, t^{\prime}\gt 0. \end{align}
Given the time dependence of the forcing term in (2.44), we proceed as in § 2.2 and apply the Laplace transform to
$ h_2^{\prime}$
. The solution for the second-order hydraulic head in integral form is expressed as
\begin{equation} h_2^{\prime} = \frac {{Al}}{2 \pi \mathrm{i}} \int _{a - \mathrm{i}\infty }^{a + \mathrm{i}\infty } \frac {\overline {F^{\prime}} \, \mathrm{e}^{s t^{\prime}} \sinh (s x^{\prime})}{s \left [ \cosh (s) + \alpha {Al} \sinh (s) \right ]} \, \mathrm{d}s, \end{equation}
where the transformed forcing term is given by
\begin{equation} \overline {F^{\prime}} = \int _0^\infty \mathrm{e}^{-s t^{\prime}} \left [ -\beta \left ( h_1^{\prime} \frac {\partial f^{\prime}}{\partial t^{\prime}} + f^{\prime} \frac {\partial h_1^{\prime}}{\partial t^{\prime}} \right ) + \alpha \frac {h_1^{\prime}}{2} \frac {\partial h_1^{\prime}}{\partial t^{\prime}} \right ] \, \mathrm{d}t^{\prime}, \end{equation}
and is known from the first-order solution obtained previously. For the moment, we assume
${Al} \neq 1$
. The integrand in (2.45) can be rewritten as the product of two functions
\begin{equation} h_2^{\prime} = \frac {1}{2 \pi \mathrm{i}} \int _{a - \mathrm{i}\infty }^{a + \mathrm{i}\infty } \overline {F^{\prime}} \, \overline {\mathcal{F}^{\prime}} \, \mathrm{d}s, \end{equation}
where
\begin{equation} \overline {\mathcal{F}^{\prime}} = \frac {{Al} \, \mathrm{e}^{s t^{\prime}} \sinh (s x^{\prime})}{s \left [ \cosh (s) + \alpha {Al} \sinh (s) \right ]} \quad \rightarrow \quad \mathcal{F}^{\prime} = {Al} \sum _{n = -\infty }^{\infty } \frac {\mathrm{e}^{s_n t^{\prime}} \sinh (s_n x^{\prime})}{s_n \left [ \sinh (s_n) + \alpha {Al} \cosh (s_n) \right ]}. \end{equation}
Applying the convolution theorem (Mei Reference Mei1997) leads to
\begin{align} h_2^{\prime} &= \int _0^{t^{\prime}} \mathcal{F}^{\prime}(t^{\prime} - \tau ^{\prime}) F^{\prime}(\tau ^{\prime}) \, \mathrm{d}\tau ^{\prime} \nonumber \\ &= {Al}\! \sum _{n = -\infty }^{\infty }\! \int _0^{t^{\prime}}\!\! \frac {\mathrm{e}^{s_n (t^{\prime} - \tau ^{\prime})} \sinh (s_n x^{\prime})}{s_n \left [ \sinh (s_n) + \alpha {Al} \cosh (s_n) \right ]} \left [ -\beta \left ( h_1^{\prime} \frac {\partial f^{\prime}}{\partial \tau ^{\prime}} + f^{\prime} \frac {\partial h_1^{\prime}}{\partial \tau ^{\prime}} \right ) + \alpha \frac {h_1^{\prime}}{2} \frac {\partial h_1^{\prime}}{\partial \tau ^{\prime}} \right ] \mathrm{d}\tau ^{\prime}. \end{align}
Expression (2.49) represents the generalised second-order solution of the hydraulic head. An explicit solution of this integral can be derived by considering the linear variation of
$f$
with time (2.33). By substituting (2.36) in (2.49) we obtain
\begin{align} h_2^{\prime}&=\frac {{Al}^2}{2}\sum _{n=-\infty }^{n=\infty }A_ns_n\left \{x^{\prime}\frac {e^{s_n t^{\prime}}-1}{s_n}+\sum _{j=-\infty }^{j=\infty }A_j\left \{\left (\frac {e^{s_jt^{\prime}}-e^{s_nt^{\prime}}}{s_j-s_n}\delta _{nj}+e^{s_nt^{\prime}}t^{\prime}\gamma _{nj}\right )\right.\right.\nonumber \\ &\quad +s_j\left [\frac {e^{s_nt^{\prime}}-e^{s_jt^{\prime}}\left (s_jt^{\prime}-s_nt^{\prime}-1\right )}{(s_j-s_n)^2}\delta _{nj}+\frac {e^{s_nt^{\prime}}{t^{\prime}}^2}{2}\gamma _{nj}\right ]\nonumber \\ &\quad\left.\left.+\frac {\alpha {Al}s_j}{2}\left [x^{\prime}\left (\frac {e^{s_jt^{\prime}}-e^{s_nt^{\prime}}}{s_j-s_n}\delta _{nj}+e^{s_nt^{\prime}}t^{\prime}\gamma _{nj}\right )+\sum _{k=-\infty }^{k=\infty }A_k\frac {e^{(s_j+s_k)t^{\prime}}-e^{s_nt^{\prime}}}{s_k+s_j-s_n}\right ]\right \}\right \}, \end{align}
where
\begin{equation} A_n=\frac {\sinh (s_nx^{\prime})}{s_n^2\left (\sinh (s_n)+\alpha {Al}\cosh (s_n)\right )},\quad \left \{ \begin{aligned} & \delta _{nj}=1,\,\gamma _{nj}=0,\qquad \text{if}\,\,n\neq j.\\ & \delta _{nj}=0,\,\gamma _{nj}=1,\qquad \text{if}\,\,n=j. \end{aligned} \right. \end{equation}
If the Allievi number is equal to unity and the valve is closing
$(\alpha =1, \beta =-1)$
, the convolution integral (2.49) modifies into
\begin{equation} h_2^{\prime} = \frac {1}{2} \int _0^{t^{\prime}} \left [ \mathrm{H}(t^{\prime} - \tau ^{\prime}) - \mathrm{H}(t^{\prime} - \tau ^{\prime} - 2) \right ] \left [ \left ( h_1^{\prime} \frac {\partial f^{\prime}}{\partial \tau ^{\prime}} + f^{\prime} \frac {\partial h_1^{\prime}}{\partial \tau ^{\prime}} \right ) + \frac {h_1^{\prime}}{2} \frac {\partial h_1^{\prime}}{\partial \tau ^{\prime}} \right ] \, \mathrm{d}\tau ^{\prime}, \end{equation}
where we used a similar approach for the solution of (2.37) because the only pole is
$ s=0$
. Substituting (2.38) into (2.52), we obtain the second-order hydraulic head for
${Al} = 1$
and linear valve closure
\begin{equation} h_2^{\prime} = \frac {1}{32} \left [ 5 {t^{\prime}}^2 + (t^{\prime} - 4)(5 t^{\prime} - 8) \mathrm{H}(t^{\prime} - 4) - 2 (t^{\prime} - 2)(5 t^{\prime} - 4) \mathrm{H}(t^{\prime} - 2) \right ]. \end{equation}
Given the second-order solutions (2.50) and (2.53), we proceed as in the previous section by evaluating the hydraulic head at
$ x^{\prime} = 1$
. Let us first consider the case of a closing valve. Since
$ h_2^{\prime}$
represents the second-order correction, we evaluate (2.50) and (2.53) at
$ t^{\prime} = 2$
and
$ t^{\prime} \to \infty$
. Using series expansion formulas for hyperbolic functions (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007), and after some algebraic manipulation, we obtain
\begin{equation} h_2^{\prime}(1,2) = \frac {{Al}^2(2 + 3{Al})}{(1 + {Al})^3}, \quad h_2^{\prime}(1, \infty ) = \frac {{Al}^2}{2}. \end{equation}
For the opening valve case, we can proceed similarly to obtain
\begin{equation} h_2^{\prime}(1,2) = 2 {Al}^2, \end{equation}
showing that the second-order contribution invariably leads to an increase in the hydraulic head since expressions (2.54) and (2.55) are always positive.
2.4. Maximum and minimum values of hydraulic head at the valve
Let us now add the leading-order solution
$h_0=H$
, the first-order solutions (2.39)–(2.40) and the second-order solutions (2.54)–(2.55) and evaluate corresponding maxima and minima of
$h(L,t)$
for the closure and opening valve cases, respectively. We obtain in dimensional variables
\begin{align} \frac {h_{{max}}}{H}& =\frac {h(L,\theta )}{H}=1+\frac {\theta {Al}}{T(1+{Al})}\left [2+\frac {\theta {Al}(2+3{Al})}{T(1+{Al})^2}\right ], \quad \text{valve closure},\,\,{Al}\in [0,\mathcal{A}], \end{align}
\begin{align} \frac {h_{{max}}}{H} & =\frac {h(L,\infty )}{H}=1+\frac {\theta {Al}}{T}\left [1+\frac {\theta {Al}}{2T}\right ], \qquad\qquad\qquad\quad \text{valve closure},\,\,{Al}\gt \mathcal{A}, \end{align}
\begin{align} \frac {h_{{min}}}{H} & =\frac {h(L,\theta )}{H}=1-\frac {\theta {Al}}{T}\left [2-\frac {2\theta {Al}}{T}\right ], \qquad\qquad\qquad\quad \text{valve opening}, \end{align}
where
$\mathcal{A}$
is the Allievi number satisfying the condition
$ h(L,\theta ) = h(L,\infty )$
approximately given by
\begin{equation} {Al} \sim \mathcal{A} = \frac {8T + 7\theta }{8T + 5\theta }, \end{equation}
which is slightly greater than unity and increases monotonically with
$\theta /T$
. In other words, if the valve is closing, maximum of hydraulic head occurs at
$t\sim \theta$
if
${Al}\in [0,\mathcal{A}]$
and at
$t\rightarrow \infty$
if
${Al}\gt \mathcal{A}$
.
Note that (2.56)–(2.58) are multiplied by
$ \theta /T = \epsilon = {St}/2$
; hence, the larger the Strouhal number, the greater the hydraulic head fluctuations, provided that
$ {St} \ll 1$
. A similar result was recently reported by Kurihara, Kiyama & Tagawa (Reference Kurihara, Kiyama and Tagawa2025), who conducted an experimental campaign to investigate the pressure fluctuations of liquids in a column subjected to short-time acceleration as a function of
$ {St}$
.
Expressions (2.56)–(2.58) are of practical relevance and can be used to estimate the maximum and minimum hydraulic heads and pressures at the valve in pipes subject to slow valve manoeuvres. These results are novel and exhibit notable differences compared with existing expressions found in the literature. For instance, (2.56) and (2.57) are much more accurate than the widely known Allievi–Michaud approximation for a closing valve (Michaud Reference Michaud1878; Allievi Reference Allievi1903; Marchi & Rubatta Reference Marchi and Rubatta2004)
\begin{equation} \frac {h_{{max}}}{H} = 1 + \frac {2\theta {Al}}{T}, \end{equation}
which is extensively used in civil engineering design standards to obtain a first estimate of the maximum hydraulic head in pipes in presence of water hammer effects caused by slow valve closure.
3. Results and discussion
In this section, we first investigate the behaviour of the non-dimensional hydraulic head
$ h/H = 1 + \epsilon h^{\prime}_1 + \epsilon ^2 h^{\prime}_2$
evaluated at
$ x = L$
in the case of valve closure and a linear time-dependent valve manoeuvre. Figure 2(a) shows the evolution of
$ h(L,t)/H$
versus non-dimensional time
$ 2t/\theta$
for
$ \epsilon = \theta /T=0.1$
and three Allievi numbers
$ {Al} = [0.5, 1, 1.5]$
. In the same figure, we compare the analytical solution with numerical results obtained using the scheme developed by Allievi (Allievi Reference Allievi1903), which gives the exact value of
$ h/H$
at the valve. The agreement between the curves is excellent, with almost complete overlap, demonstrating the reliability of our analytical model and the correctness of (2.56)–(2.57). The same figure shows that, when
$ {Al} \lt 1$
, the hydraulic head reaches its maximum at
$ t = \theta$
, i.e. when the first shock wave returns to the valve. Using the theoretical approximation (2.56), we estimate the corresponding value to be
$ h_{{max}}/H \approx 1.069$
, which agrees well with Allievi’s approach. The curve then oscillates about a value that can be approximated by taking the limit of (2.36) and (2.50) as
$ t \to \infty$
. We obtain
\begin{equation} \lim _{t^{\prime} \to \infty } \frac {h}{H} = 1 + \epsilon {Al} + \frac {\epsilon ^2 {Al}^2}{2}, \end{equation}
which for
$ {Al} = 0.5$
gives
$ h/H \approx 1.051$
. When
$ {Al} = 1$
,
$ h/H$
increases until
$ t = \theta$
, then slightly decreases until
$ t = 2\theta$
and remains constant afterwards. This behaviour is also predicted by the explicit solutions (2.38) and (2.53). In this case, the maximum value is
$ h_{{max}}/H = 1.106$
, and the asymptotic value is
$ h/H = 1.105$
, as obtained from (3.1). For a greater Allievi number,
$ {Al} = 1.5$
, the curve increases monotonically and reaches a maximum at large times. Using the theoretical approximation (2.57), we find
$ h_{{max}}/H = 1.161$
, which is again consistent with figure 2(a).
In the valve opening case shown in figure 2(b), the results for
$ {Al} = [0.5, 1, 1.5]$
agree with the analytical prediction only for
$ t \lt 2\theta$
. At later times, the solution diverges owing to the undamped oscillatory modes arising for
$ \alpha = 0$
in the first-order mode solution (2.36). This indicates that, for valve opening at sufficiently large
$ {Al}$
, alternative methods are required to capture the long-time dynamics. Nevertheless, the minima occurring at
$ t = \theta$
are very well captured, allowing the minimum hydraulic head to be estimated using (2.58). We obtain
$ h_{{min}}/H = [0.905, 0.820, 0.745]$
for
$ {Al} = [0.5, 1, 1.5]$
, values that closely match those from Allievi’s numerical method.
Non-dimensional hydraulic head at the valve versus non-dimensional time for
$ \epsilon = \theta /T= 0.1$
and various Allievi numbers
$ {Al} = [0.5, 1, 1.5]$
: (a) valve closure and (b) valve opening cases.
We now analyse the maximum and minimum hydraulic head at the valve by comparing the analytical results (2.56)–(2.58) with those obtained via Allievi’s procedure. Figure 3(a) displays the maximum hydraulic head for the valve closure case using (2.56)–(2.57), for
$ \epsilon = \theta /T=[0.05, 0.1, 0.15, 0.2]$
and
$ {Al} \in [0, 4]$
. The agreement between analytical and numerical curves is excellent, further confirming the validity of the second-order theory. In the valve opening case, figure 3(b) shows that, within the considered range of
$ {Al}$
, the agreement between the curves (2.58) and Allievi’s method is good only for small
$ \epsilon \leqslant 0.1$
. For
$ 0.2\geqslant \epsilon \geqslant 0.1$
, the theoretical model yields accurate results only when
$ {Al} \lt 2$
. We point out that this result is consistent with the perturbation expansion (2.19), which remains accurate provided that
$ \epsilon \ll O(1)$
(Mei Reference Mei1997).
(a) Maximum values of non-dimensional hydraulic head versus Allievi number
${Al}$
for the closing valve case and (b) minimum values of hydraulic head versus
${Al}$
when the valve is opening. The small parameters are
$\epsilon =\theta /T=[0.05,0.1,0.15,0.2]$
.
4. Conclusions
We have presented a mathematical model for water hammer induced by slow valve manoeuvres. The nonlinear governing equations were solved using a perturbation scheme combined with the Laplace transform, leading to explicit space- and time-dependent solutions for the hydraulic head in the pipe. The leading-order problem corresponds to the steady-state solution, while the first-order boundary-value problem yields wave-like modes that are either damped or undamped, depending on whether valve closure or opening is considered.
The analysis was extended to second order and compared with numerical solutions obtained from Allievi’s model, which provides the exact solution at the valve. For closure times larger than the critical time, the agreement between the analytical and numerical results is excellent, confirming the reliability of the theoretical approach. In contrast, during valve opening with large
$ {Al}$
, the analytical solution remains accurate only at short time scales owing to the undamped character of the first-order solution.
At this order of approximation, we also derived practical expressions for estimating the maximum and minimum values of the hydraulic head at the valve. These expressions are of significant engineering interest and may be adopted in the design of pipeline systems as they provide improved accuracy over other widely used approximated formulas, such as the classical Allievi–Michaud approximation.
The current model assumes no fluid damping or viscoelastic behaviour in the pipes, both of which are known to attenuate oscillations in the hydraulic head, particularly over longer time scales. While these limitations are significant, they are not critical for design purposes, as the first peak surge typically governs the design criteria. Furthermore, system properties such as pipe cross-section and wave speed are treated as constant, although they may vary spatially in practice. The impact of these variations will be examined in future work.
Declaration of interests
The authors report no conflict of interest.
Open access
