The spectral element method as an efficient tool for transient simulations of hydraulic systems

TL;DR

This paper demonstrates that the spectral element method provides highly accurate and computationally efficient transient simulations of hydraulic systems, enabling real-time control and optimization in engineering applications.

Contribution

It introduces the use of the spectral element method for hydraulic system simulations, highlighting its accuracy, efficiency, and suitability for real-time control applications.

Findings

SEM yields highly accurate approximations at modest computational costs.

Boundary conditions are imposed weakly using a stable upwind discretization.

The method effectively simulates water hammer effects and pressure wave evolution.

Abstract

This paper presents transient numerical simulations of hydraulic systems in engineering applications using the spectral element method (SEM). Along with a detailed description of the underlying numerical method, it is shown that the SEM yields highly accurate numerical approximations at modest computational costs, which is in particular useful for optimization-based control applications. In order to enable fast explicit time stepping methods, the boundary conditions are imposed weakly using a numerically stable upwind discretization. The benefits of the SEM in the area of hydraulic system simulations are demonstrated in various examples including several simulations of strong water hammer effects. Due to its exceptional convergence characteristics, the SEM is particularly well suited to be used in real-time capable control applications. As an example, it is shown that the time evolution of pressure waves in a large scale pumped-storage power plant can be well approximated using a low-dimensional system representation utilizing a minimum number of dynamical states.

Figures

Tables

Table 1. Table 1: Numerical flux components q ∗ superscript 𝑞 q^{*} and h ∗ superscript ℎ h^{*} corresponding to different Dirichlet boundary conditions.

boundary condition	numerical flux component $q^{\ast }$	numerical flux component $h^{\ast }$
$h(0,t)=h_0\left(t\right)$	$q_1^{\ast }=q_1-\frac{gA\left(0\right)}{c\left(0\right)}\left(h_1-h_0\left(t\right)\right)$	$h_1^{\ast }=h_0\left(t\right)$
$q(0,t)=q_0\left(t\right)$	$q_1^{\ast }=q_0\left(t\right)$	$h_1^{\ast }=h_1-\frac{c\left(0\right)}{gA\left(0\right)}\left(q_1-q_0\left(t\right)\right)$
$h(L,t)=h_L\left(t\right)$	$q_J^{\ast }=q_J-\frac{gA\left(L\right)}{c\left(L\right)}\left(h_L\left(t\right)-h_J\right)$	$h_J^{\ast }=h_L\left(t\right)$
$q(L,t)=q_L\left(t\right)$	$q_J^{\ast }=q_L\left(t\right)$	$h_J^{\ast }=h_J-\frac{c\left(L\right)}{gA\left(L\right)}\left(q_L\left(t\right)-q_J\right)$

Equations

$$h=\frac{p}{\rho g}+x_2,$$

$$\frac{\partial }{\partial t}\left[\begin{array}{c}h\\ q\end{array}\right]+\left[\begin{array}{cc}0& c^2/(gA)\\ gA& 0\end{array}\right]\frac{\partial }{\partial z}\left[\begin{array}{c}h\\ q\end{array}\right]=-\left[\begin{array}{c}0\\ Rq|q|\end{array}\right],$$

$$R=\frac{f_{\lambda }}{2DA}.$$

$$\frac{\partial}{\partial t}\begin{bmatrix}\epsilon h\\ \mu q\end{bmatrix}+\frac{\partial}{\partial z}\bigg{(}\begin{bmatrix}0&1/\mu\\ 1/\epsilon&0\end{bmatrix}\begin{bmatrix}\epsilon h\\ \mu q\end{bmatrix}\bigg{)}=-\begin{bmatrix}0\\ rq|q|\end{bmatrix},$$

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}(\mathbf{U})}{\partial z}=-\mathbf{f}(\mathbf{U}),$$

$$\mathbf{f}(\mathbf{U})=\left[\begin{array}{c}0\\ r{\mu }^{-1}{U}_2|{\mu }^{-1}{U}_2|\end{array}\right]$$

$$\mathbf{U}=\left[\begin{array}{c}{U}_1\\ U_2\end{array}\right]=\left[\begin{array}{c}ϵh\\ \mu q\end{array}\right]$$

$$\mathbf{F}(\mathbf{U})=B\mathbf{U}=\left[\begin{array}{cc}0& 1/\mu \\ 1/ϵ& 0\end{array}\right]\left[\begin{array}{c}{U}_1\\ U_2\end{array}\right].$$

$$\int_{0}^{L}\Big{[}\frac{\partial\bm{U}}{\partial t}+\frac{\partial\bm{F}(\bm{U})}{\partial z}+\bm{f}(\bm{U})\Big{]}v\,dz=\bm{0}\;\;\;\textrm{for all}\;\;\;v\in\mathcal{V}.$$

$$\int_{0}^{L}\frac{\partial\bm{U}}{\partial t}v\,dz+\bm{F}(\bm{U})v\Big{|}_{0}^{L}-\int_{0}^{L}\bm{F}(\bm{U})\frac{\partial v}{\partial z}\,dz+\int_{0}^{L}\bm{f}(\bm{U})v\,dz=\bm{0}\;\;\;\textrm{for all}\;\;\;v\in\mathcal{V}.$$

$$\bm{F}^{*}(\bm{U}^{-},\bm{U}^{+})=\frac{1}{2}\big{[}\bm{F}(\bm{U}^{-})+\bm{F}(\bm{U}^{+})\big{]}-\frac{1}{2}|\lambda_{\mathrm{max}}(B)|\big{(}\bm{U}^{+}-\bm{U}^{-}\big{)},$$

$${\mathbf{U}}^{-}(z)=\underset{\delta \to 0^{+}}{\lim }\mathbf{U}(z-\delta ),{\mathbf{U}}^{+}(z)=\underset{\delta \to 0^{+}}{\lim }\mathbf{U}(z+\delta )$$

$${\mathbf{F}}^{\ast }({\mathbf{U}}^{-},{\mathbf{U}}^{+})=\frac{1}{2}\left[\begin{array}{c}(U_2^{-}+U_2^{+})/\mu \\ (U_1^{-}+U_1^{+})/ϵ\end{array}\right]-\frac{c}{2}\left[\begin{array}{c}{U}_1^{+}-U_1^{-}\\ U_2^{+}-U_2^{-}\end{array}\right]=\frac{1}{2}\left[\begin{array}{c}{q}^{-}+q^{+}\\ h^{-}+h^{+}\end{array}\right]-\frac{c}{2}\left[\begin{array}{c}ϵh^{+}-ϵh^{-}\\ \mu q^{+}-\mu q^{-}\end{array}\right]$$

$${\int }_0^L\frac{\partial }{\partial t}\left[\begin{array}{c}ϵh\\ \mu q\end{array}\right]vdz={\int }_0^L\left[\begin{array}{c}q\\ h\end{array}\right]\frac{\partial v}{\partial z}dz-{\int }_0^L\left[\begin{array}{c}0\\ rq|q|\end{array}\right]vdz+\left[\begin{array}{c}{F}_1^{\ast }(0,t)\\ F_2^{\ast }(0,t)\end{array}\right]v(0)-\left[\begin{array}{c}{F}_1^{\ast }(L,t)\\ F_2^{\ast }(L,t)\end{array}\right]v(L)$$

$$\Omega =\bigcup _{m=1}^M{Z}_m$$

$${\mathcal{V}}_h=\{v\in C^0(\Omega ):v_{|Z_m}\in {\mathbb{P}}_{N_m},m=1,\dots ,M\}.$$

$${\ell }_i^{(m)}(\xi )=\prod _{k=0,k\ne i}^{N_m}\frac{(\xi -{\xi }_k)}{({\xi }_i-{\xi }_k)},\xi \in \hat{\Omega },i=0,\dots ,N_m$$

$${\ell }_i^{(m)}({\xi }_j)={\delta }_{i,j},j=0,\dots ,N_m.$$

$${\Gamma }_m(\xi )=\frac{{\bar{z}}_m-{\bar{z}}_{m-1}}{2}\xi +\frac{{\bar{z}}_m+{\bar{z}}_{m-1}}{2},\xi \in \hat{\Omega }$$

$${\phi }_i^{(m)}(z)={\ell }_i^{(m)}({\Gamma }_m^{-1}(z)),i=0,\dots ,N_m$$

$$J=\sum _{m=1}^M{N}_m+1$$

$$\mathcal{Z}=\{z_1,\dots ,z_J\},$$

$${\varphi }_i(z_j)={\delta }_{i,j},z_j\in \mathcal{Z},i=1,\dots ,J.$$

$$h(z,t)\approx h_h(z,t)=\sum _{j=1}^J{h}_j(t){\varphi }_j(z)\text{and}q(z,t)\approx q_h(z,t)=\sum _{j=1}^J{q}_j(t){\varphi }_j(z).$$

$$Q_h(z,t)=\sum _{j=1}^J{Q}_j(t){\varphi }_j(z),Q_j(t)=q_j(t)|q_j(t)|,j=1,\dots ,J$$

$$(M_{\chi }{)}_{i,j}={\int }_0^L\chi (z){\varphi }_i{\varphi }_{j}dz,$$

$$S_{i,j}={\int }_0^L\frac{\partial {\varphi }_i}{\partial z}{\varphi }_{j}dz,$$

$$\mathbf{h}={\left[\begin{array}{c}{h}_1,\dots ,h_J\end{array}\right]}^{\top },\mathbf{q}={\left[\begin{array}{c}{q}_1,\dots ,q_J\end{array}\right]}^{\top },\mathbf{Q}={\left[\begin{array}{c}{q}_1|q_1|,\dots ,q_J|q_J|\end{array}\right]}^{\top }$$

$${\mathbf{e}}_1={\left[\begin{array}{c}1,0,\dots ,0\end{array}\right]}^{\top },{\mathbf{e}}_J={\left[\begin{array}{c}0,\dots ,0,1\end{array}\right]}^{\top }$$

$$M_{ϵ}\frac{d\mathbf{h}}{dt}$$

Full text

1 Introduction

Fast variations of the volume flow in a pipeline system give rise to pressure waves propagating through the pressure lines and causing problems like noise, vibrations, or in the worst case, even pipe collapse. These transient hydraulic effects constitute a major risk in various areas like cooling-water systems, pumped-storage power plants or even in the feed lines of liquid rocket engines [1, 2, 3, 4].

Wave propagation in pressure lines is well described by a system of hyperbolic partial differential equations [5, 1, 6, 7, 2, 8]. Various numerical methods for the simulation of hyperbolic systems are known in literature, ranging from finite difference to finite volume and spectral methods [9, 10, 11]. Moreover, discontinuous Galerkin methods [12] represent an attractive alternative as they try to combine the advantages of the finite volume method with those of the finite element method.

However, the predominant method for the simulation of pressure waves in pipeline systems is the method of characteristics (MOC) [1, 6, 13, 7, 2]. It is based on the concept of the same name used to analyze and solve systems of hyperbolic partial differential equations. The MOC works best when all parameters characterizing the pressure lines are piecewise constant. Under this condition, the MOC represents an extremely robust method which is able to propagate shock waves with arbitrary steep gradients. The main disadvantage results from the fact that the time step size, the spatial grid size and the local wave propagation speed are required to satisfy a (strict) Courant-Friedrichs-Lewy (CFL) condition. To meet these requirements it is common practice to adjust the lengths or the local wave speeds of individual pipeline segments which may cause significant modeling errors. Moreover, in case of non-constant coefficient functions a large number of spatial grid points is needed in order to compute satisfactory numerical approximations.

Besides the MOC, finite difference, finite volume, discontinuous Galerkin and spectral methods have been used to simulate transient hydraulic systems [14, 15, 16, 17, 18, 19]. In this context, it should be noted that numerous variations of spectral methods can be found in literature [11, 20]. Ideally, spectral methods yield highly accurate solutions requiring only a minimum number of degrees of freedom, which makes them attractive from a control applications point of view. However, as spectral methods are formulated on a single interval or on a multidimensional cuboid, they turn out to be too inflexible in many applications.

The spectral element method (SEM) takes this yet a step further and applies a spectral method to every element of a properly partitioned computational domain. In this manner, it is possible to combine the adaptivity of an ordinary finite element method with the accuracy of a spectral method [21]. In fact, the SEM could be interpreted as a special variation of a high-order finite element method. However, a unique feature of the SEM is the quadrature rule, which relies on the same nodes as the definition of the element shape functions. As a consequence, the mass matrix is a diagonal matrix, which enables the application of fast explicit time stepping methods. Furthermore, thanks to the quadrature rule and the nodal interpolation property the assembling process of nonlinear expressions can often be simplified significantly.

The SEM is widely used in elastic wave propagation and fluid dynamical problems [22, 23, 24]. In this paper, we apply the SEM to a system of hyperbolic partial differential equations frequently used to describe hydraulic systems in engineering applications. As for all hyperbolic equations, special care needs to be taken with regard to the implementation of the boundary conditions. Here, we choose to impose the boundary conditions weakly using a numerically stable upwind discretization. As outlined above, the resulting mass matrix is a diagonal matrix. Moreover, we show that the nonlinear pipe friction term can be easily taken into account. Hence, the application of explicit time stepping strategies is very simple and yields extremely short computation times as needed in real-time control applications.

A minor drawback of the SEM might be seen in the fact that, in contrast to the MOC or discontinuous Galerkin methods, the SEM cannot propagate shock waves with arbitrary steep gradients. Nonetheless, we will demonstrate that the SEM yields highly accurate solutions to rather strong water hammer problems using a modest number of spatial variables. Additionally, it is shown that the underlying spatial semi-discretization can be integrated in time utilizing high-order explicit time-integration methods in combination with comparatively mild time step size restrictions.

Our aim is to demonstrate the benefits of the SEM in the area of transient hydraulic simulations, cf. also [25]. Moreover, we want to show that, due to its exceptional convergence properties, the SEM is particularly well suited to be used in optimal and model predictive control applications utilizing direct transcription (full discretization) methods [26, 27].

The paper is organized as follows: In Section 2, we introduce the underlying hyperbolic system which describes wave propagation effects on a single pressure line. Section 3 covers the weak formulation, the spatial semi-discretization and the implementation of Dirichlet boundary conditions. In Section 4, we provide a careful numerical convergence analysis using several examples of increasing complexity. In this context, we also consider the above-mentioned strong water hammer problems. Finally, Section 5 is devoted to the simulation of wave effects in a pumped-storage power plant. In particular, it is demonstrated that the SEM can be used to obtain low-dimensional semi-discretizations which are perfectly suited to be used in real-time capable control applications.

2 Mathematical Modelling

Wave propagation in pressure lines is well described in the literature [1, 6, 7, 2, 5]. Using mass and momentum balance equations it is possible to derive a system of hyperbolic partial differential equations. A basic assumption is that the liquid flow is one-dimensional, i.e., the characteristic quantities are cross-section averaged [13]. In most applications, the average velocity of the fluid is small in comparison to the characteristic wave speeds $c$, which is why certain convective terms can be neglected. Moreover, in many cases it is advantageous to consider the piezometric head $h$ instead of the pressure $p$. Both quantities are related to each other in the form

[TABLE]

where $x_{2}=x_{2}(z)$ denotes the height of the pressure line as function of the arclength $z$. Furthermore, all simulations in this paper are based on the numerical values $g=9.81\,\mathrm{m}/\mathrm{s}$ and $\rho=1000\,\mathrm{kg}/\mathrm{m}^{3}$ for the gravitational acceleration and the density of water, respectively.

The time-evolution of the piezometric head $h$ and volume flow $q$ on a single pressure line is well described by the hyperbolic system [2, 1]

[TABLE]

where pipe friction is taken into account by the source term on the right hand side using a position-dependent function

[TABLE]

Here, $f_{\lambda}$ denotes the Darcy-Weisbach friction factor and $D=2(A/\pi)^{1/2}$ is the diameter of the pressure line with $A$ the cross-sectional area. Strictly speaking, this description of pipe friction is valid only for stationary pipe flow [6]. However, it is frequently used as a good approximation in transient simulations as well [1, 2]. A description of pipe frictional losses for non-stationary simulations is given in [5, 6]. In many situations of practical interest, the effect of pipe frictional losses is very small. In particular, this will apply if fast transient simulations on a relatively short time-scale are considered. However, in stationary calculations, e.g., calculation of optimal stationary operating points of pumped-storage power plants, even small contributions of pipe frictional losses might be important and thus need to be taken into account.

3 Numerical Method

We now show how to apply the spectral element method to the hyperbolic system (2). To this end, we first restrict ourselves to the case of a single pressure line described by the parameter functions $A,c$ and $f_{\lambda}$ given on the computational domain $\Omega=[0,L]$.

3.1 Weak formulation

Introducing $\epsilon=gA/c^{2}$, $\mu=1/(gA)$ and $r=\mu R$, system (2) can be written in conservative form (with source)

[TABLE]

i.e.,

[TABLE]

where

[TABLE]

denotes the friction term as a function of the new variables

[TABLE]

and the flux function is given by

[TABLE]

Let $\mathcal{V}$ denote a suitable space of test functions on the interval $[0,L]$. The weak formulation of (4) reads

[TABLE]

Using integration by parts, (6) becomes

[TABLE]

Next, we replace the flux at the boundaries by an appropriate numerical approximation. In order to obtain a numerically stable and convergent scheme, we choose a so-called upwind flux. In particular, we choose the Lax-Friedrichs flux [20] given as

[TABLE]

where

[TABLE]

and $|\lambda_{\mathrm{max}}(B)|$ denotes the maximum (in absolute value) eigenvalue of the matrix $B$ defined in (5). Assuming that $\epsilon^{-}=\epsilon^{+}$ and $\mu^{-}=\mu^{+}$ holds true at the boundary points $z=0$ and $z=L$, we find

[TABLE]

and hence, (7) yields

[TABLE]

for all $v\in\mathcal{V}$. Here, $F_{1}^{*}(z,t)$ and $F_{2}^{*}(z,t)$ denote the components of the numerical flux $\bm{F}^{*}$ evaluated at the position $z$ and time $t$.

3.2 Spatial discretization

In the following, we consider a partition of the computational domain into $M$ elements

[TABLE]

with $Z_{m}=[\bar{z}_{m-1},\bar{z}_{m}]$ and $0=\bar{z}_{0}<\bar{z}_{1}<\dots<\bar{z}_{M-1}<\bar{z}_{M}=L$. Moreover, we define the finite-dimensional space

[TABLE]

A precise description of a basis of $\mathcal{V}_{h}$ in the sense of the SEM can be found in the literature, see, e.g., [20]. As in the case of the classical FEM the definition of the basis functions is given by means of the element shape functions [20, 28], which in turn are defined in terms of the Legendre-Gauss-Lobatto (LGL for short) nodes $\xi_{0},\dots,\xi_{N_{m}}$ on the reference domain $\hat{\Omega}=[-1,1]$. The LGL-nodes are the zeros of $(1-\xi)^{2}P_{N_{m}}^{\prime}(\xi),$ where $P_{N_{m}}^{\prime}$ denotes the derivative of the Legendre polynomial $P_{N_{m}}$ of degree $N_{m}$, cf. [11]. Using the LGL nodes the element shape functions can be written as

[TABLE]

with

[TABLE]

As an example, the element shape functions for $N_{m}=4$ are shown in Fig. 1.

By means of the map

[TABLE]

we may define the element basis functions

[TABLE]

on each element $Z_{m}=(\bar{z}_{m-1},\bar{z}_{m})$ with $m=1,\dots,M$. The element basis functions are now used to construct the global basis functions $\phi_{1},\dots,\phi_{J}$. To this end, all element basis functions are extended by zero outside of their domains. Moreover, in order to ensure inter-element continuity at the nodes $\bar{z}_{1},\dots,\bar{z}_{M-1}$, adjacent element basis functions are suitably combined into one single basis function. Fig. 2 depicts an example with three elements and the polynomial degrees $N_{1}=2$, $N_{2}=4$ and $N_{3}=3$. In general, we construct

[TABLE]

global basis functions $\phi_{j}$, $j=1,\dots,J$ representing a basis for the ansatz space $\mathcal{V}_{h}$.

Using (12) the LGL nodes of the element shape functions can be mapped to the domains of the corresponding elements. Subsequently, we collect all nodes in

[TABLE]

where each inter-element node appears only once. An important feature of the SEM is the nodal interpolation property:

[TABLE]

In the next step, we replace $h$ and $q$ in the weak formulation (9) by the approximations

[TABLE]

Moreover, we temporarily introduce

[TABLE]

to approximate the nonlinear expression $q(z,t)|q(z,t)|$ in the second component of (9). In the sense of a Galerkin method, the space of test functions $\mathcal{V}$ is replaced by the space of ansatz functions $\mathcal{V}_{h}$. With the help of the mass matrix

[TABLE]

the stiffness matrix

[TABLE]

as well as the vectors

[TABLE]

and

[TABLE]

we obtain

[TABLE]

and $\bm{q}|\bm{q}|$ denotes a pointwise multiplication of two column vectors. In (17), we have introduced the notations

[TABLE]

which are motivated by the observation that $F_{1}^{*}(z,t)$ and $F_{2}^{*}(z,t)$ are approximations of $q(z,t)$ and $h(z,t)$, respectively, cf. (5). The values of $q_{1}^{*},h_{1}^{*},q_{J}^{*}$ and $h_{J}^{*}$ are used to realize the boundary conditions and will be specified in the subsequent sections.

Like in the case of classical finite element methods, the assembling process of the mass and stiffness matrices is realized by means of the corresponding element mass and element stiffness matrices [29, 28, 30]. In case of the SEM all integrals are approximated by the quadrature rule

[TABLE]

Hence, the assembling process is based on the same nodes $\xi_{0},\dots,\xi_{N_{m}}$ as the definition of the element shape functions. In fact, it can be easily verified that, as a result of the nodal interpolation property (14) and the quadrature rule (18), any mass matrix $M_{\chi}$ (irrespective of the coefficient function $\chi$) is a diagonal matrix [11]. This property allows for the application of fast explicit time stepping methods to the system of ordinary differential equations (ODEs) given in (17). The simplicity of the nonlinear pipe frictional loss term $M_{r}\bm{q}|\bm{q}|$ is even more remarkable. Indeed, the mass matrix $M_{r}$ has to be computed only once. This is in strong contrast to classical finite element methods which in the worst case scenario would require to assemble a new mass matrix in each time step. Finally, we would like to mention that the expression in (16) coincides with the nonlinear expression $q_{h}(z,t)|q_{h}(z,t)|$ for $z=z_{j}$, $j=1,\dots,J$. Therefore, $Q_{h}(z,t)$ and $q_{h}(z,t)|q_{h}(z,t)|$ are equivalent with respect to the quadrature rule (18) and hence (16) does not introduce an additional disretization error.

3.3 Dirichlet boundary conditions

One of the most frequently used boundary conditions are Dirichlet boundary conditions

[TABLE]

for some given (potentially time dependent) functions $h_{0}(t),q_{0}(t),h_{L}(t)$ and $q_{L}(t)$. In general, boundary conditions are realized using the values $q_{1}^{*},h_{1}^{*}$ and $q_{J}^{*},h_{J}^{*}$ in (17) corresponding to the components of the numerical flux (8)

[TABLE]

evaluated at the boundary points $z=0$ and $z=L$, respectively.

As an example let us consider the Dirichlet boundary condition $h(L,t)=h_{L}(t)$. A simple strategy reported in literature of Discontinuous Galerkin methods (cf. [12, 31]) is given by

[TABLE]

i.e.,

[TABLE]

where $q^{-}$ and $h^{-}$ have to be replaced by the last nodal values $q_{J}$ and $h_{J}$, respectively. In the same way, all remaining boundary conditions of (19) can be realized, see Tab. 1.

3.4 General remarks concerning the numerical simlations

All simulations of this article are realized using the Matlab programming language. In most examples, reference solutions are provided by the MOC. In this context, it is important to note that due to the simple structure of the underlying algorithm, the MOC can be implemented very efficiently using Matlab-specific vectorization techniques. However, as the MOC is a low-order numerical method, very small time step sizes $(\triangle t)_{\mathrm{MOC}}$ are required. In its most basic form (without interpolations) the MOC requires the spatial step size $(\triangle z)_{\mathrm{MOC}}$ to satisfy the (strict) CFL condition $(\triangle z)_{\mathrm{MOC}}/(\triangle t)_{\mathrm{MOC}}=c$. The situation becomes even worse in pressure lines consisting of $K$ segments characterized by different wave speeds $c_{k}$, $k=1,\dots,K$ where the MOC needs to satisfy

[TABLE]

simultaneously for all segments $k=1,\dots,K$. As a result, only certain combinations of spatial mesh sizes and time step sizes are possible.

A crucial step in the implementation of the SEM is the computation of the element mass and element stiffness matrices, which require the nodes $\xi_{0},\dots,\xi_{N_{m}}$ and the weights $w_{0},\dots,w_{N_{m}}$ corresponding to the quadrature rule given in (18). Moreover, the first derivatives of the element shape functions $\ell_{0}^{(m)},\dots,\ell_{N_{m}}^{(m)}$ need to be evaluated at $\xi_{0},\dots,\xi_{N_{m}}$ in order to assemble the stiffness matrix. Computing these numerical values with high precision is a nontrivial task, especially if the polynomial degree $N_{m}$ is large. We therefore make use of the Matlab routines given in [32], which are numerically stable and efficient at the same time.

Using (14) and (15) we immediately find $h_{h}(z_{j},t)=h_{j}$ and $q_{h}(z_{j},t)=q_{j}$ for all $j=1,\dots,J$. However, in order to compare solutions of the SEM with those of the MOC we need to evaluate $h_{h}$ and $q_{h}$ at the (densely sampled) grid points associated to the MOC. This interpolation process is realized using barycentric interpolation [33], which is a fast and stable variant of Lagrange polynomial interpolation. It is important to note that a direct evaluation of $h_{h}$ and $q_{h}$ based on formula (10) is numerically unstable and would ruin all numerical convergence rates reported below.

Finally, we would like to mention that in the graphical representations we always show the spatial distributions of the pressure $p$ instead of the piezometric head $h$. However, since the numerical implementation of the SEM and the MOC are based on (2) which is formulated in terms of the piezometric head $h$ and the volume flow $q$, we show numerical errors and relative differences in terms of the piezometric head $h$. The numerical differences in terms of $q$ would look very similar and yield the same convergence rates.

4 Numerical convergence analysis

In the following, we consider several numerical examples of increasing complexity. The simulations in this section are limited to a single horizontal pressure line. Without loss of generality, we set $x_{2}(z)=0$ for all $z\in\Omega$, $\Omega=[0,L]$. The length of the pressure line $L$ and its diameter $D$ are similar to those used in [19].

4.1 A mathematical benchmark problem

In the first example, we consider constant coefficient functions $D,c>0$ and $f_{\lambda}=0$. In particular, we are interested in computing the time evolution corresponding to the initial values

[TABLE]

for $z\in\mathbb{R}$ and $t>0$. Using the characteristic decomposition

[TABLE]

it follows immediately that

[TABLE]

solve (2) for all $z,t\in\mathbb{R}$. Here, $w_{1,s}$ and $w_{2,s}$ denote the initial values given by $w_{1,s}(z)=[q_{s}(z)+(gA/c)h_{s}(z)]/2$ and $w_{2,s}(z)=[q_{s}(z)-(gA/c)h_{s}(z)]/2$, respectively.

At the boundaries of $\Omega$ we implement the following pair of transparent boundary conditions:

[TABLE]

Obviously, these boundary conditions are of different nature compared to those listed in Tab. 1. However, we may define the functions $q_{0}(t)=-(gA/c)h_{1}$ and $q_{L}(t)=(gA/c)h_{J}$ with the first and last nodal value $h_{1}$ and $h_{J}$, respectively. Subsequent application of the second and fourth line of Tab. 1 yields

[TABLE]

We note that these equations coincide exactly with the Silver-Müller boundary conditions (reduced to one spatial dimension) presented in [31]. Below we will see that (24) can be used in combination with explicit time stepping methods and reasonable time step sizes $\triangle t$. However, we found that a simple modification of (24) results in a stable discretization for even larger time step sizes. To this end, we require (23) to be satisfied for the numerical flux values $q^{*}$ and $h^{*}$, i.e. we define $q_{0}(t)=-(gA/c)h_{1}^{*}$ and $q_{L}(t)=(gA/c)h_{J}^{*}$. Application of the second and fourth line in Tab. 1 yields

[TABLE]

For the numerical experiment we choose $D=0.01$m, $c=1200\,\mathrm{m}/\mathrm{s}$ and $L=12$ m. Furthermore, we consider the initial value distributions

[TABLE]

for $z\in\Omega$ with $\alpha=100$ m and $\beta=1$.

In the first numerical simulation, the computational domain is partitioned into only $M=10$ elements of equal size and a polynomial degree of $N_{m}=5$ is used for all elements $m=1,\dots,M$. In order to propagate the numerical solution in time we employ the classical Runge-Kutta method of order $4$ [34]. An illustration of the time evolution of $p$ and $q$ for $t\in[0,T]$ with $T=8$ ms is given in Fig. 3. Using the numerical flux components given in (25) we may utilize a time step size of $\triangle t=0.2$ ms which corresponds to only $40$ time steps for the whole simulation. In contrast, a time step size of $\triangle t=0.05$ ms is needed to ensure numerical stability if the numerical flux values are chosen according to (24).

As a next step we want to compute numerical convergence rates with respect to the size of the elements. To this end, we divide the computational domain into $M$ elements of the same size $\triangle z$. Moreover, we employ the same polynomial degree $N_{m}=N$ for all $m=1,\dots,M$. During a series of simulations we keep the polynomial degree constant whereas the computational domain is divided into decreasing element sizes $\triangle z$. In order to keep the time discretization error as low as possible we use the Matlab function ode45 ($\texttt{RelTol}=2.5\times 10^{-14}$). In this specific example, an exact reference solution is available and can be computed using (22). The relative errors are evaluated at time $T=5$ ms, i.e., exactly when half of the wave packet has already left the computational domain. In this manner, it is guaranteed that the effect of the boundary conditions is fully included in the numerical convergence rates depicted in the left column of Fig. 4. It can clearly be seen that the convergence rates are of order $N+1$ for odd degree and $N$ for even degree. In this context, it is worth noting that such odd-even phenomena have already been reported in literature of discontinuous Galerkin methods (DGM) [35, 36] and that the way we implement the boundary conditions is based on the same principles as in the DGM.

The effect reported above might be interesting from a theoretical point of view, however, we will now demonstrate that for realistic element sizes $\triangle z$ it makes hardly any difference whether one uses odd or even polynomial degrees. For this purpose, we partition the computational domain into $M=10$ elements of the same size. In each simulation, the polynomial degree $N_{m}=N$ is successively increased for all elements $m=1,\dots,M$ and the relative errors are gathered in the right column of Fig. 4. The resulting graph perfectly agrees with the anticipated spectral, i.e., exponential convergence rate. Furthermore, there is a significant improvement of the accuracy irrespective of whether one increases $N$ starting from an even or an odd polynomial degree.

4.2 Water hammer simulation 1

We consider the same pressure line as in the previous example. However, the boundary condition at the left boundary point $z=0$ is replaced with a nonhomogenous Dirichlet boundary condition $h(0,t)=h_{0}$, where $h_{0}=p_{0}/(\rho g)$ for some given pressure $p_{0}$. At the same time the boundary condition corresponding to the right boundary point $z=L$ is replaced with the boundary condition

[TABLE]

which describes the turbulent volume flow through a valve. Here, $p_{v}$ denotes the pressure at the valve outlet, $A_{v}$ is the cross-sectional area, $\alpha_{v}$ is referred to as the contraction coefficient, and $u(t)$ describes the control input related to the opening of the valve. Assuming that $p(L,t)>p_{v}$ holds true for all times $t$, we may write

[TABLE]

with $C_{v}=\alpha_{v}\sqrt{2g}A_{v}$ and $h_{v}=p_{v}/(\rho g)$.

In a water hammer simulation, the initial values of $h$ and $q$ are given by the stationary solution of (2) and the control input $u$ is decreased from a positive value $u_{0}$ at time $t=0$ to $u_{T_{c}}=0$ at time $t=T_{c}$ in a strictly monotonous manner. In order to prevent a distortion of the convergence rates reported below, the control input needs to be a sufficiently smooth function. Therefore, we employ

[TABLE]

where $\sigma$ denotes the eighth order sharpened raised cosine filter [11]

[TABLE]

with

[TABLE]

We note that $u$ is seven times continuously differentiable in $\mathbb{R}$.

The first boundary condition at the left boundary point $z=0$ can be implemented easily using the first line of Tab. 1. Similar to the previous example we consider two methods to realize the second boundary condition in (26). To this end, we consider the components of the numerical flux

[TABLE]

corresponding to the last line of Tab. 1. In case of the first method, (26) is evaluated using the last nodal value of $h$:

[TABLE]

In case of the second approach, we require (26) to be satisfied for the numerical flux values $q_{J}^{*}$ and $h_{J}^{*}$ itself. Therefore, $q_{L}(t)=q_{J}^{*}$ is given implicitly via

[TABLE]

In a numerical implementation, $h_{J}^{*}$ is replaced according to the second expression in (28) and hence we obtain a quadratic equation in $q_{J}^{*}$. However, as we will see below, the numerical effort caused by solving an additional quadratic equation is negligibly small compared to the improvements regarding the stability of the resulting numerical scheme.

The following numerical experiments are based on the parameter values $p_{0}=120$ bar, $p_{v}=100$ bar, $\alpha_{v}=0.7$ and $A_{v}=A_{L}/5$, where $A_{L}$ denotes the cross-sectional area at $z=L$. The control input is given according to (27) using a valve closure time of $T_{c}=5$ ms, see Fig. 5. All other parameters $D=0.01$m, $c=1200\,\mathrm{m}/\mathrm{s}$, $L=12$ m and $f_{\lambda}=0$ remain unchanged.

In the first numerical experiment, we partition the computational domain into $M=10$ elements of equal size and a polynomial degree of $N_{m}=5$ is employed on each element $m=1\dots,10$. For the implementation of the valve boundary condition, we employ the numerical flux values according to the implicit formulation (30). The time integration is realized using the classical Runge Kutta method in combination with a time step size of $\triangle t=0.2$ ms which corresponds to only $200$ time steps for the whole time evolution shown in Fig. 6. In contrast, if the valve boundary condition is realized using (29), a time step size of $\triangle t=0.025$ ms is needed to ensure numerical stability. Thus, a minor modification of the discretization scheme improves the computational efficiency by a factor of almost ten.

Next, we take a look at the water hammer signature, i.e., the time evolution of the pressure at $z=L$. In particular, we compare the water hammer signatures obtained by means of the SEM with those obtained using the MOC. In order to satisfy the CFL condition of the MOC (for all considered element sizes), a comparatively small time step size of $\triangle t=0.025$ ms is used for the time integration. Note that this restriction only applies to the MOC. In fact it would be possible to use a much larger time step size for the time integration of the SEM. As before, the time integration of the SEM is based on the classical Runge-Kutta method.

In Fig. 7, we show the results of the SEM using a different number of equally sized elements $M$ with fixed polynomial degrees $N_{m}=3$ for all $m=1,\dots,M$. While for small $M$ we observe significant overshoots similar to the Gibbs phenomenon, these overshoots decrease rapidly with increasing $M$. In fact, for $M=40$ we perfectly recover the water hammer signature computed by the MOC.

The results of a complementary situation are depicted in Fig. 8. Here, we consider a fixed number of elements $M=10$ but different polynomial degrees $N_{m}=N$ with $m=1,\dots,M$. It can clearly be seen that, as the polynomial degree is increased, the solutions corresponding to the SEM converge quickly towards the solution obtained by the MOC.

Considering very similar simulations in [19], the results presented above are quite remarkable. Using only $M=10$ elements in combination with a polynomial degree of $N=5$, we are able to perfectly recover the water hammer signature computed by the MOC. This is in strong contrast to the Chebyshev super spectral viscosity method [19] which, by means of an additional viscosity operator, is able to counteract the Gibbs-like phenomenon. However, at the same time, the viscosity operator leads to a significant smoothing of the original water hammer signature. In this context, it is important to note that the system sizes of both discretizations are comparable. Moreover, the water hammer effect considered here is very similar in terms of the ratio of the valve reflection time $(2L/c)$ to the valve closure time $T_{c}$.

Finally, we compute numerical convergence rates corresponding to the SEM as a function of the element size $\triangle z$ and the polynomial degree $N_{m}$. The time evolution is again based on the Matlab function ode45 ($\texttt{RelTol}=2.5\times 10^{-14}$) which allows to keep the time discretization error as low as possible.

Due to the nonlinear valve boundary condition (26), an exact reference solution would be extremely difficult to obtain. We therefore compute a numerical reference solution using the MOC. The relative differences will be evaluated at time $T=20$ ms. Since the MOC is only a low order method, an extremely small time step size of $(\triangle t)_{\mathrm{MOC}}=0.00025$ ms is needed in order to obtain sufficiently accurate solutions, and as a result of the CFL condition (21), the number of spatial grid points $M_{\mathrm{MOC}}=40\,000$ is quite large. Thus, computing a reference solution at time $T=20$ ms is a rather time-consuming process.

The left column of Fig. 9 shows the relative differences as a function of the element size $\triangle z$ for fixed polynomial degrees $N_{m}=N$ for all $m=1,\dots,M$. As we can see, we obtain the same numerical convergence rates as in the mathematical benchmark problem (cf. Fig. 4) exhibiting the same odd-even phenomenon. However, we will again show that in realistic scenarios, where the number of elements is much smaller, this odd-even phenomenon is not an issue. To this end, let us consider a subdivision of the computational domain into $M=10$ equally sized elements. While the number of elements is kept constant we compute the relative differences for increasing polynomial degrees $N_{m}=N$ with $m=1,\dots,M$. The results are gathered in the right column of Fig. 9. Like in the mathematical benchmark problem we observe a nearly perfect spectral convergence rate. Moreover, no staircase effect is visible and, therefore, the improvement of accuracy is not dependent on whether we increase $N$ starting from an odd or an even polynomial degree.

4.3 Water hammer simulation 2

So far, we neglected the effect of pipe friction losses. Furthermore, the diameter of the pressure line considered in the previous sections was assumed to be constant. Here, we will focus our attention on a more general situation, namely a pressure line which is characterized by the coefficient functions $D$, $c$ and $f_{\lambda}$ as depicted in Fig. 10. All the other parameters are chosen like in the previous water hammer problem and also the boundary conditions and the time evolution of the control input $u$ remain the same.

The pressure line consists of three segments (separated by dashed orange lines) as needed to resolve the abrupt changes in the coefficient functions. In the first numerical simulation, the segments are subdivided into $M_{1}=10$, $M_{2}=1$ and $M_{3}=10$ elements, see Fig. 10. Moreover, a polynomial degree of $N=3$ is assigned to all elements. Please note that the element sizes vary slightly from one segment of the pressure line to another. In order to propagate the numerical solution in time we utilize the classical Runge Kutta method in combination with a time step size of $\triangle t=0.2$ ms. The corresponding time evolution of $p$ and $q$ is shown in Fig. 11 and exhibits much more complicated patterns compared to the first water hammer problem (cf. Fig. 6).

Taking a closer look at the derivation presented in Sec. 3, it is quite obvious that it makes hardly any difference whether the SEM is used in combination with constant or piecewise smooth coefficient functions as long as abrupt changes in $D$, $c$, or $f_{\lambda}$ are taken into account by a properly adapted spatial mesh. This is in strong contrast to the MOC which in case of a variable diameter $D$ performs significantly worse. Also a non-zero pipe friction coefficient $f_{\lambda}$ deteriorates the accuracy of the MOC. However, since pipe frictional losses are comparatively small, this effect becomes only visible at very fine spatio-temporal discretizations.

For the reasons indicated above, it is basically impossible to obtain a sufficiently accurate solution by means of the MOC which will then serve as a reference solution to investigate the numerical convergence rates of the SEM. Instead, we will show that for decreasing step sizes $(\triangle t)_{\mathrm{MOC}}$ the solutions of the MOC converge towards a reference solution provided by the SEM. Here, it should be noted that the spatial mesh sizes $(\triangle z)_{\mathrm{MOC}}^{(k)}$ are completely determined by the wave speeds on the different segments $k=1,2,3$ via the CFL condition (21).

In order to compute a reference solution by means of the SEM, the first, second and third segment is partitioned into $40$, $4$ and $40$ elements, respectively. Moreover, a polynomial degree of $N_{m}=8$ is used for all elements $m=1,\dots,M$ and the time evolution is based on the Matlab function ode45 ($\texttt{RelTol}=2.5\times 10^{-14}$). It is worth noting that a further refinement of the spatio-temporal discretization would not change the relative differences shown in Fig. 12.

It can clearly be seen that, for practically relevant time step sizes $(\triangle t)_{\mathrm{MOC}}$, the numerical solutions corresponding to the MOC converge at a rate of approximately $2$ towards the reference solution provided by the SEM. However, for small time step sizes the convergence rate decreases to a value of approximately $1$. This behavior agrees well with the numerical approximations underlying the MOC. The derivation of the MOC is based on an integration of the characteristic equations along the positive and negative characteristics $dz/dt=c$ and $dz/dt=-c$, respectively. To this end, variable cross-sectional areas are replaced by average values corresponding to two consecutive grid points, which causes the second order convergence rate shown in Fig. 12. In order to integrate the nonlinear pipe frictional loss term, the expression $q|q|$ is approximated by its value at the beginning of each time step. While this simple approximation is necessary in order to make the MOC an explicit method, it is also responsible for the first order convergence rate seen in Fig. 12. However, since pipe frictional losses are very small the effect is only visible at a low level of relative differences. In fact, if pipe frictional losses are neglected ($f_{\lambda}=0$) we obtain a uniform quadratic convergence rate, cf. Fig. 12.

By means of the above simulations we have, at least implicitly, demonstrated that the SEM works as expected also in the case of variable coefficient functions $D$, $c$ and $f_{\lambda}$. Furthermore, we have investigated the numerical convergence behavior of the MOC. To the best of our knowledge, and despite the fact that the MOC is the most widely used method in the context of transient hydraulic simulations, the convergence analysis of the MOC presented above cannot be found in literature in this form. We would like to emphasize that in the numerical example, the wave speeds and the lengths of the segments were chosen in such a way that the CFL conditions in (21) can be easily satisfied. In a real-world example the lengths of the segments often have to be modified with respect to the original data, and hence, the numerical solutions of the MOC converge towards the solution of a (slightly) modified problem.

5 Simulation of hydraulic transients in a pumped-storage power plant

Having demonstrated the accuracy and effectiveness of the SEM for the numerical simulation of wave propagation effects on a single pressure line, we will now turn our attention to a coupled system of several pressure lines. In particular, we want to simulate the transient behavior of pressure waves in a large-scale pumped-storage power plant. A detailed mathematical model of a variable speed pumped-storage power plant including a comprehensive description of the underlying hydraulic and electrical systems is presented in [37].

The main reason for considering the pumped-storage power plant example is to study the applicability and the advantage of the SEM for a real-world control problem. For the application of modern optimization-based control techniques like model predictive control, appropriate numerical approximations which are computationally very efficient and at the same time capable of accurately describing the time evolution of the pressure waves are a key element.

5.1 Model description

The spatial positions of the pressure lines, the bifurcations and the turbines are illustrated in Fig. 13 using a $(x_{1},x_{2},x_{3})$-coordinate frame. Here, $x_{2}$ is oriented in opposite direction to the gravitational acceleration $g$. Each of the pressure lines will be parameterized by its arclength $z_{i}\in[0,L_{i}]$, with $L_{i}$ denoting the length of the $i$th pressure line, $i=1,\dots,6$. Moreover, we assume that the 1st pressure line is connected to a reservoir at $z_{1}(0)$ and the 6th pressure line is connected to a reservoir at $z_{6}(L_{6})$.

Fig. 14 shows the coefficient functions $\alpha$, $D$, $c$ and $f_{\lambda}$ corresponding to the 1st, 2nd, 4th and 6th pressure line. The coefficient functions of the 3rd and 5th pressure line coincide with those of the 2nd and the 4th pressure line, respectively. We note that the inclination $\alpha_{i}$ determines the height profile

[TABLE]

of the $i$th pressure line as a function of the local arclength $z_{i}$. However, these height profiles are only needed for the conversion of pressure values to piezometric heads (and vice versa), see (1). The time evolution of the piezometric head $h_{i}$ and the volume flow $q_{i}$ on the $i$th pressure line is governed by (2), i.e.

[TABLE]

for $i=1,\dots,6$. In order to complete the description of the mathematical model, we need to specify twelve boundary conditions.

Given the relatively short simulation periods considered below, it is reasonable to assume that the filling levels of the reservoirs are virtually constant for all $t\in[0,T]$. Therefore, the piezometric heads at the reservoirs are fixed by Dirichlet boundary conditions

[TABLE]

where

[TABLE]

with $p_{\mathrm{res,t}}$ and $p_{\mathrm{res,b}}$ being the pressure at the top and bottom reservoir, respectively.

Furthermore, since the lengths of the bifurcations are short in comparison to the lengths of the pressure lines, it is fair to assume that the bifurcations can be well described by the relations

[TABLE]

The remaining four boundary conditions are

[TABLE]

where the volume flows through the turbines $q_{\mathrm{I}}$ and $q_{\mathrm{II}}$ are given implicitly by the nonlinear algebraic equations [37]

[TABLE]

Here, $\omega_{\mathrm{I}}$ and $\omega_{\mathrm{II}}$ denote angular velocities corresponding to the two turbines and

[TABLE]

are dimensionless variables which, together with the guide vane positions $\chi_{\mathrm{I}}$ and $\chi_{\mathrm{II}}$, are used to evaluate the characteristic map $W_{H}$ shown in Fig. 15. We note that the above equations describe two Francis turbines which are assumed to be of the same type and, therefore, identical parameters $h_{\mathrm{ref}}$, $q_{\mathrm{ref}}$ and $\omega_{\mathrm{ref}}$ are used both in (35a) and (35b). Moreover, in order to keep the model as simple as possible, we assume that the angular velocities of the turbines are held constant at $\omega_{\mathrm{I}}=\omega_{\mathrm{II}}=\omega_{\mathrm{ref}}$ by means of a subordinate controller.

5.2 Semi-discretization

According to (17), the semi-discretizations of the systems given in (31) read

[TABLE]

Here, it is important to note that $h_{2}(L_{2},t),h_{4}(0,t)$ in (35a) and $h_{3}(L_{3},t),h_{5}(0,t)$ in (35b) have been replaced by $h_{2,J_{2}}^{*}$, $h_{4,1}^{*}$ and $h_{3,J_{3}}^{*}$, $h_{5,1}^{*}$, respectively. Another possibility would be to use the nodal values $h_{2,J_{2}}$, $h_{4,1}$ and $h_{3,J_{3}}$, $h_{5,1}$ directly. However, as in the case of the nonlinear valve boundary condition considered in Sec. 4.2, the latter choice requires a much smaller time step size $\triangle t$ in order to ensure numerical stability in combination with explicit time stepping methods.

5.3 Numerical simulation

We now consider a numerical simulation of the pumped-storage power plant introduced above. In particular, we want to demonstrate that the SEM is able to reproduce a solution provided by the MOC using a minimum number of spatial grid points. In this context, it is important to note that the relative errors reported in Sec. 4 are disproportionally small when compared to the expected errors of the underlying model. Our long-term objective is to operate pumped-storage power plants using model predictive control (MPC). In such applications, it is crucial to have a low-dimensional system representation available, which describes the dynamic behavior sufficiently well while remaining numerically efficient enough to be evaluated in real time.

The lengths of the pressure lines shown in Fig. 14 represent realistic values but have been modified slightly in order to be able to satisfy the CFL conditions (21) required by the MOC. The largest time step size which satisfies all CFL conditions simultaneously is $(\triangle t)_{\mathrm{MOC}}=2\times 10^{-3}$ s. Since the spatial mesh sizes of the MOC are determined by (21) it follows that $J_{\mathrm{MOC}}=2\times 1036$ spatial grid points are needed for the discretization of the physical quantities $h$ and $q$. Note that an even finer spatio-temporal discretization would not change the numerical differences reported below.

In case of the SEM, we employ a single element per segment, see Fig. 14. Thus, in the given example, segments and elements coincide. While our implementation allows for arbitrary subdivisions of the segments, it turned out that single elements of comparatively high order are more efficient when it comes to reaching medium level accuracies using a minimum number of spatial grid points. The polynomial degrees of the element basis functions corresponding to the different pressure lines and segments are $\mathcal{N}_{1}=(3,1,3,1,2)$, $\mathcal{N}_{2}=\mathcal{N}_{3}=1$, $\mathcal{N}_{4}=\mathcal{N}_{5}=2$ and $\mathcal{N}_{6}=(1,1,8)$. Thus, $J_{\mathrm{SEM}}=2\times 32$ degrees of freedom are used to represent the approximations of $h$ and $q$, cf. (13). Since (36) constitutes a semi-explicit system of differential algebraic equations (DAE) of index 1, there are basically two established approaches to realize the time integration. The first approach consists of solving the algebraic equations (36c) locally for the algebraic variables $q_{\mathrm{I}}$ and $q_{\mathrm{II}}$. Hence, the remaining system of equations (36a) and (36b) can be treated using explicit time stepping methods like the classical Runge Kutta method or ode45 as implemented in the Matlab scripting language. Alternatively, the whole system is treated as a single DAE system which can be integrated in time using implicit methods like ode15s which is also an integral part of the Matlab programming environment.

The numerical simulation is based on the following set of physical parameters: $p_{\mathrm{res,t}}=5\,\mathrm{bar}$, $p_{\mathrm{res,b}}=2\,\mathrm{bar}$, $h_{\mathrm{ref}}=313.4\,\mathrm{m}$, $q_{\mathrm{ref}}=44.7\,\mathrm{m}^{3}/\mathrm{s}$, and $\omega_{\mathrm{ref}}=42.1\,\mathrm{rad}/\mathrm{s}$. Moreover, we assume that the initial state of the system is given by a stationary solution of (36) corresponding to the initial guide vane positions $\chi_{\mathrm{I}}(t=0)$ and $\chi_{\mathrm{II}}(t=0)$. The guide vane positions $\chi_{\mathrm{I}}$ and $\chi_{\mathrm{II}}$ will be varied as shown in the last row of Fig. 16 for the simulation period $[0,T]$ with $T=20$ s.

Due to the restrictions imposed by the CFL conditions, the MOC requires $10\,000$ time steps in order to compute the complete time evolution. The elapsed calculation time amounts to approximately $4.4$ seconds. In contrast, if the time evolution is computed by means of the SEM and ode45 ($\texttt{RelTol}=2.5\times 10^{-6}$) the required calculation time amounts to approximately $1.1$ seconds. Using ode15s ($\texttt{RelTol}=2.5\times 10^{-6}$), the computing time can be even further reduced to roughly $0.75$ seconds.

The first row of Fig. 16 illustrates the time evolution of the pressure at the upstream side of the first and second turbine, respectively. Due to the short lengths of the 2nd, 3rd, 4th and 5th pressure line the differences of both graphs are comparatively small. The second row shows the time evolution of the volume flows through the turbines. It can be clearly seen that even within time intervals where the guide vane positions are held constant the volume flows in the turbines are far from being constant as well. In fact, it would take much longer time intervals ($\approx 20$ s) until the system becomes (approximately) stationary again. Additionally, the third and fourth row of Fig. 16 demonstrate that the corresponding (absolute) differences between the solutions of the SEM and the MOC are well below $5\times 10^{-3}$ for all times $t\in[0,T]$.

So far, we have only considered the differences of the numerical solutions corresponding to the SEM and the MOC at or in the two Francis turbines. However, in many applications the numerical solutions are required to be accurate along the whole pipeline system. For this reason, the relative differences between the numerical solutions corresponding to the SEM and the MOC are given in Fig. 17. More precisely, the relative differences of the piezometric heads and the volume flows are shown as a function of time and separately for all six pressure lines. Quite obviously, the relative differences of the piezometric heads $h_{1},\dots,h_{6}$ and the volume flows $q_{1},\dots,q_{6}$ are well below $1\times 10^{-2}$ for all times $t\in[0,T]$ and hence the given level of accuracy is significantly smaller than the expected overall modeling error. We note that, in the present case, the time evolution corresponding to the SEM was calculated using the Matlab function ode45. However, the second approach based on ode15s yields practically the same results.

Of course, one could object that faster variations of the guide vane positions would cause stronger pressure fluctuations which could not be resolved by the low-dimensional discretization of the SEM outlined above. However, it turns out that the variations of $\chi_{\mathrm{I}}$ and $\chi_{\mathrm{II}}$ already represent a possible worst case scenario. To make this clear, we consider the time evolution of the pressure on the 6th pressure line as illustrated in Fig. 18. Due to its length the 6th pressure line shows the most pronounced variations of the pressure. Moreover, due to its spatial position, the pressure on the 6th pressure line is most likely to approach zero, which, with respect to cavitation, is a critical scenario. In fact, we find that the time evolution of $p$ and $q$ is comparatively smooth. However, Fig. 19 unveils that at $t=10$ s the pressure becomes very close to zero (locally). In a real system, the pressure should never fall below $1$ bar (lower dashed line) and due to model inaccuracies even larger lower bounds are realistic. At the same time, the pressure along the pressure lines is required to remain below a reasonable upper bound, which is defined in terms of the pressure distribution of the stationary state (orange dashed line). However, Fig. 19 shows that even a generously dimensioned upper bound (upper dashed line) is dramatically exceeded at $t=16$ s.

In optimal or model predictive control applications, these lower and upper bounds for the pressure distributions must be systematically taken into account. In this context, the most rigorous approach is based on direct transcription methods [26, 27]. However, in order to keep the numerical effort as low as possible, direct transcription methods are dependent on a low dimensional system representation. Comparing the system size of the MOC ($J_{\mathrm{MOC}}=2\times 1036$) with the system size of the SEM ($J_{\mathrm{SEM}}=2\times 32$), it is obvious that the SEM will give a significant advantage in this field of applications.

6 Conclusions and outlook

In summary, this paper shows that the spectral element method (SEM) represents an efficient and flexible method for transient simulations of hydraulic systems, in particular also for the use in optimization-based control applications. By means of the SEM it is possible to compute highly accurate numerical solutions of the water hammer problem using modest computational resources. As a high order method, the SEM requires a minimum number of dynamical states in order to yield excellent approximations of sufficiently smooth pressure and volume flow fluctuations. In fact, we demonstrated that the hydraulic system of a large scale pumped-storage power plant can be well approximated using a low-dimensional system representation which consists of only a few dozen dynamical states. The corresponding semi-discretization can be integrated in time using standard time stepping strategies and is therefore well suited for optimal and model predictive control applications. In particular, model predictive control of pumped-storage power plants, while respecting lower and upper pressure bounds, is within reach and will be presented in a subsequent publication.

Acknowledgements

The authors would like to thank M. Egretzberger and M. Meusburger from Andritz Hydro GmbH for the fruitful cooperation. This work was supported by the Austrian Research Promotion Agency (FFG), Grant No.: 849267. Furthermore, the financial support by the Austrian Federal Ministry of Science, Research and Economy and the National Foundation for Research, Technology and Development is gratefully acknowledged.