Strong Stability Preserving Runge-Kutta Methods Applied to Water Hammer Problem

ABSTRACT

The characteristic method of lines is the most used numerical method applied to the water hammer problem. It transforms a system of partial differential equations involving the independent variables time and space in two ordinary differential equations along the characteristics curves and then solve it numerically. This approach, although showing great stability and quick execution time, creates ∆x-∆t dependency to properly model the phenomenon. In this article we test a different approach, using the method of lines in the usual form, without characteristics curves and then applying strong stability preserving Runge-Kutta Methods aiming to get stability with greater ∆t.

Keywords:
method of characteristic; method of lines; strong stability preserving methods; Water Hammer

1 INTRODUCTION

Pipe networks may have abrupt changes in flow rate due to valves being closed or the action of hydraulic pumps, causing also abrupt pressure variations. This phenomenon is commonly called in literature by ”Water Hammer”. The Water Hammer phenomenon is highly studied in engineering. According ²⁰, high variations of hydraulic head in a duct may damage it in the short term. Ducts should de designed with correct diameter, wall size thickness or allow the use of water hammer control devices. We need then to have a good prediction of the pressures the duct will suffer. The equations that model the phenomenon is also used to pre-locate water loss in pipe networks due to local damages in ducts as seen in ⁸.

In this work, we apply strong stability preserving Runge-Kutta methods (SSPRK) ^{5), (7} in to a system of partial differential equations modeling the water hammer problem. A simple model of a reservoir connected to a duct and this to a valve that closes abruptly was used (Fig. 1). At the left side of the duct, hydraulic head is kept constant while at the right side, flow rate abruptly goes to 0. Hydraulic head and flow rate is assumed to be known at time $t=0$, as well as expressions for the boundary conditions. It is then desired to predict the flow rate and the hydraulic head at each point of the duct over time.

Figure 1:
Hydraulic Sistem.

Based on ^{2), (4), (8), (19}, the System of Differential Equation 1.1 represent the most common equations used to model the problem. Reference ¹⁴ makes a physical deduction for these equations. Wichowski ²⁰ work with a general form of these equations, changing the term |Q|Q to |Q|Q ^m−1 , but we use the simple form.

$\left\{\begin{array}{cc}{H}_t+\frac{a^2}{gA}{Q}_x& =0\\ Q_t+gA{H}_x+\frac{f}{2DA}Q\left|Q\right|& =0\end{array}\right.$ (1.1)

In the System 1.1, $H(x,t)$ is the hydraulic head and $Q(x,t)$ the flow rate in a given position x of the duct at time t. Constant f represents the friction factor, A is the cross-sectional area of the duct, g the gravitational acceleration, D the diameter of the duct and a a constant called pressure velocity. The initial and boundary conditions are described in Section 2.

Currently, the most used numerical method apllied to the problem is the characteristic method of lines, usualy called simply method of characteristics (MOC) ^{1), (4), (19), (20}. There are several different kind of MOC methods, each one with its particularities and own characteristics curves. In this work we implement a simple form of the method just to comparison. This is detailed in Section 3. The MOC method establish a dependency between temporal ∆t and spatial ∆x stepsizes discretization to properly model the phenomenon. To achieve more freedom in ∆x and ∆t choice, a simple alternative to MOC, proposed in this work, is to use the method of lines (MOL) in the usual form (without the characteristics curves) to transform the partial differential equations in two variables in a system of differential equations in just one variable, and then use a known numerical method to solve it. This is described in Section 4. Trying to get good numerical solutions with greater ∆t stepsizes than in MOC scheme, we choose to adapt Strong Stability Preserving Runge-Kutta Methods. These methods are based on ^{5), (6} and it is described in Section 5. Section 6 present SSPRK schemes based on ^{5), (6), (7), (11), (12} and compares MOL and MOC methods. Section 7 concludes the work.

2 INITIAL AND BOUNDARIES CONDITIONS

We made tests for the implemented numerical methods with two examples in literature. For the first example, the initial and boundary conditions were based on an experiment described in Wichowiski’s work ²⁰ and also explored in ¹⁴. The length of the duct was set 41 meters (L = 41m), the diameter of the cross section 42 millimeters (D = 0.042m). The pressure velocity was taken as a = 1260m/s. The friction factor was considered f = 0.025 and the gravity acceleration g = 9.81m/s ². Initial flow rate was $Q(x,0)=0.000453014m^3/s$, over the entire length of the duct as well as the initial hydraulic head was $H(x,0)=50m$. The hydraulic head located on the left side of the duct remains constant at 50m, that is, $H(0,t)=50m$ for all t. At the right side of the duct, a valve is placed, with fast closing time T _c = 0.034s, beginning at time t = 0.16s. This closing time was considered to influence linearly the flow rate on the right side of the duct trough the following equation

$Q(L,t)=\left\{\begin{array}{ccc}Q0& \text{if}& t\le 0.16\\ Q0+Q0\frac{\left(0.16-t\right)}{T_c}& \text{if}& 0.16<t\le 0.16+T_c\\ 0& \text{if}& t>0.16+T_c\end{array}\right..$ (2.1)

For the second example, based in ¹³, we consider a duct with lenght L = 12000 ft and cross section diameter D = 2 ft. Pressure velocity was a = 3000 ft/s. Initial head was taken varying linearly from H0 = 600 ft in the beginning of the duct to Hf = 530 ft at the end, close to the valve. The friction fator was f = 0.02. The flow rate at the valve was described with the formula $Q(L,t)=\left(C_d{A}_v\right)H(L,t)$, where C _d was considered the coefficient of discharge and A _v the area of the valve at time t. This area changes linearly from time t = 0 up to completely closed at time t = 4. The initial flow rate was Q(x, 0) = 20.022713 ft ³ /s. The converted gravity acceleration was g = 32.18504 ft/s ².

3 METHOD OF CHARACTERISTICS - MOC

Fox ⁴ list some advantages of the characteristic method to solve the problem when compared to other methods. The characteristic method, directly applied to System 1.1 works as follows:

Multiply the first equation in System 1.1 by a number λ and add both equations. Take λ = gA/a, to get Equation 3.1.

$\frac{gA}{a}(H_t+a{H}_x)+(Q_t+a{Q}_x)+\frac{f}{2DA}Q\left|Q\right|=0.$ (3.1)

Consider a straight line x(t), with $\frac{dx}{dt}=a$ represented in Figure 2. Equation 3.1 may be seen as an equation involving total derivatives

$\frac{gA}{a}\frac{dH}{dt}+\frac{dQ}{dt}+\frac{f}{2DA}Q\left|Q\right|=0.$ (3.2)

Figure 2:
Characteristics Lines.

Equation 3.2 is valid on the straight line x(t). Replace then the total derivatives of the hydraulic head and flow rate by simple differences to get the first equation of the Linear System 3.3. In this system, H _i,j and Q _i,j are, respectively, numerical approximations for the hydraulic head and flow rate values at a point P _i,j of the discretized domain. The j index represents time variations while the i index spatial variations along the duct, as shown in Figure 2. The second equation of System 3.3 is obtained analogously, with $\lambda =-gA/a$ and $\frac{dx}{dt}=-a$. Solving System 3.3, allows us to determine the hydraulic head and flow rate of a point P _i,j+1 , indicated by $H_{i,j+1}^P$ and $Q_{i,j+1}^P$, as long as the hydraulic head and flow rate at the points P _i−1,j and P _i+1,j are known.

$\left\{\begin{array}{l}{H}_{i,j+1}^P-H_{i-1,j}+\frac{a}{gA}(Q_{i,j+1}^P-Q_{i-1,j})+\frac{af}{2gD{A}^2}\left|Q_{i-1,j}\right|Q_{i-1,j}\Delta t=0\\ H_{i,j+1}^P-H_{i+1,j}-\frac{a}{gA}(Q_{i,j+1}^P-Q_{i+1,j})-\frac{af}{2gD{A}^2}\left|Q_{i+1,j}\right|Q_{i+1,j}\Delta t=0\end{array}\right.$ (3.3)

To consider the boundary conditions, note that if we know the hydraulic head in P _i,j+1 , P _i+1,j , and the flow rate in P _i+1,j , we can find out the flow rate in P _i,j+1 using only the second equation of System 3.3. If we know the flow rate in P _i,j+1 , P _i−1,j , and the hydraulic head in P _i−1,j , we can find out the head in P _i,j+1 using only the first equation of System 3.3.

The Continuous line in Figure 3 represents the graphic of the hydraulic head obtained by MOC related to the first example and trough a discretized scheme with 20 subdivisions of the duct i.e, an amount of points n = 21. We locate the sample point at position $\overline{x}$ = 30.75m, close to the valve. With this scheme we have ∆x = 2.05m and consequently, to maintain the relationship $\frac{∆x}{∆t}=a$, ∆t = 0, 001626984. Changing this relationship changes the problem and obviously also the graphic. In Figure 4 we used n = 81. With correspoding ∆x and ∆t stepsizes parameters. Related to the second example, in Figure 5 we choose ∆x = 600 ft and located the sample point at position $\overline{x}$ = 9000 ft.

Figure 3:
First Example. Numerical Result. Hydraulic Head at $\overline{x}$ = 30.75m. ∆x = 2.05m and ∆t = 0.001626984

Figure 4:
First Example. Numerical Result. Hydraulic Head at $\overline{x}$ = 30.75m. ∆x = 0.5125m and ∆t = 0.000406746

Figure 5:
Second Example. Numerical Result. Hydraulic Head at $\overline{x}$ = 9000 f t. ∆x = 600 f t and ∆t = 0.2

4 METHOD OF LINES - MOL

The method of lines (MOL), transforms a system of differential equations in two variables, to a system of differential equations in just one variable (t). This occurs by discretizing variable x.

Consider the differential equations in System 1.1. Fix $\overline{x}$ and use the first two terms of taylor’s series expansion on x for H. We get Equation 4.1.

$H_x(\overline{x},t)\approx \frac{H(\overline{x}+\Delta x,t)-H(\overline{x},t)}{\Delta x}$ (4.1)

Placing Equation 4.1 in the the second equation of System 1.1, we get Equation 4.2

$Q_t(\overline{x},t)+gA\frac{H(\overline{x}+\Delta x,t)-H(\overline{x},t)}{\Delta x}+\frac{f}{2DA}Q(\overline{x},t)\left|Q\right(\overline{x},t\left)\right|=0$ (4.2)

Analougsly, using the first two terms of the taylor series expansion for Q, with −∆x it follows

$Q_x(\overline{x},t)\approx \frac{Q(\overline{x},t)-Q(\overline{x}-\Delta x,t)}{\Delta x}$ (4.3)

Replace then in the first equation of System 1.1 and we get Equation 4.4.

$H_t(\overline{x},t)=-\frac{a^2}{gA}\frac{Q(\overline{x},t)-Q(\overline{x}-\Delta x,t)}{\Delta x}$ (4.4)

Based on boundary conditions, we have initial values of hydraulic head in $x=0,H_0=H(0,t)$, and flow rate in $x=L,Q_L\left(t\right)=Q(L,t)$. Making x vary in the discretized points of a partition $\left\{x_0=0,x_1,x_2,\dots ,x_n=L\right\}$, we have the System 4.5. We must discover the following functions: $Q(0,t),Q(\Delta x,t),Q(2\Delta x,t),\dots ,Q\left(\right(n-1)\Delta x,t)$ and $H(\Delta x,t),H(2\Delta x,t),\dots ,H(L,t)$. That is, a system of differential equations with 2n equations and 2n functions to be discovered.

$\left\{\begin{array}{ccc}{Q}_t(0,t)& =& -\frac{f}{2DA}Q(0,t)\left|Q\right(0,t\left)\right|-gA\frac{H(\Delta x,t)-H_0}{\Delta x}\\ Q_t(\Delta x,t)& =& -\frac{f}{2DA}Q(\Delta x,t)\left|Q\right(\Delta x,t\left)\right|-gA\frac{H(2\Delta x,t)-H(\Delta x,t)}{\Delta x}\\ ⋮& ⋮& ⋮\\ Q_t\left(\right(n-1)\Delta x,t)& =& -\frac{f}{2DA}Q\left(\right(n-1)\Delta x,t)\left|Q\right((n-1)\Delta x,t\left)\right|-gA\frac{H(L,t)-H\left(\right(n-1)\Delta x,t)}{\Delta x}\\ H_t(\Delta x,t)& =& -\frac{a^2}{gA}\frac{Q(\Delta x,t)-Q(0,t)}{\Delta x}\\ H_t(2\Delta x,t)& =& -\frac{a^2}{gA}\frac{Q(2\Delta x,t)-Q(\Delta x,t)}{\Delta x}\\ ⋮& ⋮& ⋮\\ H_t(L,t)& =& -\frac{a^2}{gA}\frac{Q_L\left(t\right)-Q\left(\right(n-1)\Delta x,t)}{\Delta x}\end{array}\right.$ (4.5)

Take $y\left(t\right)=\left[\begin{array}{c}Q(0,t)\\ Q(\Delta x,t)\\ ⋮\\ Q\left(\right(n-1)\Delta x,t)\\ H(\Delta x,t)\\ H(2\Delta x,t)\\ ⋮\\ H(L,t)\end{array}\right]$, and we get then an initial value problem $y\text{'}=F(t,y),y\left(0\right)=y_0$ where

$y_0=\left[\begin{array}{c}Q(0,0)\\ Q(\Delta x,0)\\ ⋮\\ Q\left(\right(n-1)\Delta x,0)\\ H(\Delta x,0)\\ H(2\Delta x,0)\\ ⋮\\ H(L,0)\end{array}\right].$

We can then choose a numerical method to solve the equation.

5 SSPRK SCHEMES

Numerical methods applied to differential equations involving non smooth solutions often shows problems as spurious oscillations. One way to attack these problems is through strong stability preserving methods (SSP) ^{6), (7), (12), (16), (18}. The concept of SSP methods was introduced for autonomous systems by ³, but in ¹⁷ and ⁹ the concept is generalized for non-autonomous systems. These ones aims to preserve the stability of the forward Euler method but showing high accuracy order. Suppose we have a differential equation $y\text{'}=F(t,y)$, and there are $\overline{\Delta t}$ such that for all $t\in \mathrm{ℝ},y\in {\mathrm{ℝ}}^n$ and ∆t in the range $0\le \Delta t\le \overline{\Delta t}$ the forward Euler Method satisfies:

$\left|\right|y+\Delta tF(t,y)\left|\right|\le \left|\right|y\left|\right|$

then if $y_{t_n}$ and $y_{t_{n+1}}$ are two consecutive steps of a numerical method, we say it is a SSP method if there is a constant C such that for $\Delta t\le C\overline{\Delta t}$ we have

$\Vert y_{t_{n+1}}\Vert \le \Vert y_{t_n}\Vert$

In case of hypebolic autonomous equations, the constant C is called CFL coefficient or SSP coefficient.

Special kind of SSP methods are the strong stability preserving Runge-kutta methods (SSPRK). These methods are formulated for autonomous equations. In this work, MOL framework is most of time autonomous, being non-autonomous just during valve closing time. We present here a formulation for time dependent equations. Adapted from ⁵, our SSPRK scheme will be writen in the form

$y^i=y_{t_n}+\Delta t\sum _{j=1}^m{a}_{ij}F(t_n,y^j)(1\le i\le m)$ (5.1)

$y_{t_{n+1}}=y_{t_n}+\Delta t\sum _{j=1}^m{b}_{j}F(t_n,y^j).$ (5.2)

where y ⁱ represent auxiliary values such that, together with $y_{t_n}$, the next step $y_{t_{n+1}}$ are calculated. Note that time t _n in expressions $F(t_n,y^j)$ may be generalized, but this work do not make use of this generalization. The equations may also be written in the form

$y^0=y_{t_n}$ (5.3)

$y^i=\sum _{j=0}^{i-1}\left({\alpha }_{ij}{y}^j+\Delta t{\beta }_{ij}F(t_n,y^j)\right)(1\le i\le m)$ (5.4)

$y_{t_{n+1}}=y^m.$ (5.5)

To consistency purpose, we need $\sum _{j=0}^{i-1}{\alpha }_{ij}=1$. Furthermore, we will take ${\alpha }_{ij}\ge 0\forall i,j$. Equations 5.3, 5.4 and 5.5 are called the Shu-Usher representation.

According Equation 5.4, we have

$$\begin{gathered} \Vert y^i\Vert =||\sum _{j=0}^{i-1}\left({\alpha }_{ij}{y}^j+\Delta t{\beta }_{ij}F(t_n,y^j)\right)||=||\sum _{j=0}^{i-1}{\alpha }_{ij}\left(y^j+\Delta t\frac{{\beta }_{ij}}{{\alpha }_{ij}}F(t_n,y^j)\right)|| \\ \le \sum _{j=0}^{i-1}||{\alpha }_{ij}\left(y^j+\Delta t\frac{{\beta }_{ij}}{{\alpha }_{ij}}F(t_n,y^j)\right)||\le \sum _{j=0}^{i-1}{\alpha }_{ij}||\left(y^j+\Delta t\frac{{\beta }_{ij}}{{\alpha }_{ij}}F(t_n,y^j)\right)||. \end{gathered}$$

Considering forward Euler stability we have $||\left(y^j+\Delta t\frac{{\beta }_{ij}}{{\alpha }_{ij}}F(t_n,y^j)\right)||\le ||y_{t_n}||.$ if $\frac{{\beta }_{ij}}{{\alpha }_{ij}}\Delta t\le \overline{\Delta t}$, so $||y^i||\le ||y_{t_n}||$ for all i, in particular $||y^m(=y_{t_{n+1}})||\le ||y_{t_n}||$. If time step ∆t satisfies $\Delta t\le C\overline{\Delta t}$ where $C=min\left\{\frac{{\alpha }_{ij}}{{\beta }_{ij}}\right\}$, we have then a SSP method.

The SSPRK schemes will be described trough tables as presented in ¹². For example, our adapted SSPRK(3, 3) method ^{5), (12} is given by

$$\begin{gathered} y^0=y_{t_n} \\ {y}^1=y^0+\Delta tF(t_n,y^0) \\ {y}^2=\frac{3}{4}{y}^0+\frac{1}{4}{y}^1+\frac{1}{4}\Delta tF(t_n,y^1) \\ {y}_{t_{n+1}}=\frac{1}{3}{y}^0+\frac{2}{3}{y}^2+\frac{2}{3}\Delta tF(t_n,y^2). \end{gathered}$$

The first parameter in the description of the method refers to number of stages used, and the second to the order of the method. This order was deduced for autonomous systems. We will make reference to this order, but note that during valve closing time system is not autonomous and the order of the SSPRK schemes may be lower than specified by the parameters. Our purpose is to get stability with greater ∆t stepsize. The SSPRK(3, 3) scheme may be represented as Table 1.

Table 1:
SSPRK(3, 3).

stages m = 3
1	0	0	αij
$\frac{3}{4}$	$\frac{1}{4}$	0
$\frac{1}{3}$	0	$\frac{2}{3}$
1	0	0	β ij
0	$\frac{1}{4}$	0
0	0	$\frac{2}{3}$	⁝
C = 1

6 RESULTS

For comparison purpose, considering the first example, we made several tables with MOC scheme and SSPRK’s schemes with different values of hydraulic head for different positions and time. For each time, we get the values of hydraulic head truncating the discrete values of time to the nearest value specified in tables. Tables 2 and 3 shows that MOC scheme has litlle changes related with variations in ∆x stepsize. The graphics for $\overline{x}$ = 30.75 are presented in Figures 3 and 4 with continuous lines. It was also collected some values of MOC scheme implemented by ²⁰ and for an experiment also presented in ²⁰ for comparison.

Table 2:
MOC - First Example n = 21, ∆t = 0.001626984.

	t = 0.4002	t = 0.6003	t = 1.2007	t = 1.4008	t = 1.9995
H(2.05,t)	46.368557844	53.453654047	47.228802928	50.075671819	49.936437764
H(10.25,t)	31.885577986	67.198397789	41.047034021	54.465531060	49.682221321
H(20.50,t)	18.116331293	83.967954550	33.546768788	61.514505102	49.814750124
H(30.75,t)	12.332717431	85.401035193	26.242311040	68.614126035	49.574290268

Table 3:
MOC - First Example n = 81, ∆t = 0.000406746.

	t = 0.4002	t = 0.6003	t = 1.2007	t = 1.4008	t = 1.9995
H(2.05,t)	46.365612181	53.455680166	47.228095566	50.075712447	49.936399875
H(10.25,t)	31.870830369	67.208632451	41.044319136	54.466820459	49.682032796
H(20.50,t)	18.093539172	83.990549701	33.541124933	61.518467815	49.814686930
H(30.75,t)	12.308769424	85.423426506	26.229078974	68.621308310	49.579782151

We made tests with different SSPRK schemes using method of lines described in section 5. To show the results, although several tests made, we choose only two methods of third order, SSPRK(3, 3) and SSPRK(8, 3) respectively, and two methods of fourth order, SSPRK(5, 4) and SSPRK(10, 4). Some numerical results are presented in Tables 7, 8, 9, 10, 11 and 12. Figures 3 and 4 shows graphics comparing MOC scheme implemented in this work with SSPRK(8,3) for n = 21 and n = 81.

Table 4:
SSPRK(8, 3).

stages m = 8
1.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	α ij
0.00000	1.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
0.00000	0.00000	1.00000	0.00000	0.000000	0.00000	0.00000	0.00000
0.00000	0.00000	0.00000	1.00000	0.00000	0.00000	0.00000	0.00000
0.42137	0.00595	0.00000	0.00000	0.57268	0.00000	0.00000	0.00000
0.00000	0.00425	0.00000	0.00000	0.00000	0.99575	0.00000	0.00000
0.00000	0.00000	0.10438	0.24327	0.00000	0.00000	0.65235	0.00000
0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	1.00000
0.19580	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	β ij
0.00000	0.19580	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
0.00000	0.00000	0.19580	0.00000	0.00000	0.00000	0.00000	0.00000
0.00000	0.00000	0.00000	0.195804	0.00000	0.00000	0.00000	0.00000
0.00000	0.00000	0.00000	0.00000	0.11213	0.00000	0.00000	0.00000
0.00000	0.00000	0.00000	0.00000	0.00000	0.19497	0.00000	0.00000
0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.12773	0.00000
0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.19580
C = 5:10714756443533

Table 5:
SSPRK(5, 4).

stages m = 5
1.00000	0.00000	0.00000	0.00000	0.00000	α ij
0.44437	0.55563	0.00000	0.00000	0.00000
0.62010	0.00000	0.37990	0.00000	0.000000
0.17808	0.00000	0.00000	0.82192	0.00000
0.00000	0.00000	0.51723	0.09606	0.38671
0.39175	0.00000	0.00000	0.00000	0.00000	β ij
0.00000	0.36841	0.00000	0.00000	0.00000
0.00000	0.00000	0.25189	0.00000	0.00000
0.00000	0.00000	0.00000	0.54497	0.00000
0.00000	0.00000	0.00000	0.06369	0.22601
C = 1:50818004918983

Table 6:
SSPRK(10, 4).

stages m = 10
1	0	0	0	0	0	0	0	0	0	α ij
0	1	0	0	0	0	0	0	0	0
0	0	1	0	0	0	0	0	0	0
0	0	0	1	0	0	0	0	0	0
$\frac{3}{5}$	0	0	0	$\frac{2}{5}$	0	0	0	0	0
0	0	0	0	0	1	0	0	0	0
0	0	0	0	0	0	1	0	0	0
0	0	0	0	0	0	0	1	0	0
0	0	0	0	0	0	0	0	1	0
$\frac{1}{25}$	0	0	0	$\frac{9}{25}$	0	0	0	0	$\frac{3}{5}$
$\frac{1}{6}$	0	0	0	0	0	0	0	0	0	β ij
0	$\frac{1}{6}$	0	0	0	0	0	0	0	0
0	0	$\frac{1}{6}$	0	0	0	0	0	0	0
0	0	0	$\frac{1}{6}$	0	0	0	0	0	0
0	0	0	0	$\frac{1}{15}$	0	0	0	0	0
0	0	0	0	0	$\frac{1}{6}$	0	0	0	0
0	0	0	0	0	0	$\frac{1}{6}$	0	0	0
0	0	0	0	0	0	0	$\frac{1}{6}$	0	0
0	0	0	0	0	0	0	0	$\frac{1}{6}$	0
0	0	0	0	$\frac{3}{50}$	0	0	0	0	$\frac{1}{10}$
C = 6

Table 7:
SSPRK(8,3) - First Example. n = 21, ∆t = 0.001626984.

	t = 0.4002	t = 0.6003	t = 1.2007	t = 1.4008	t = 1.9995
H(2.05,t)	46.491927125	52.910010319	47.828096312	51.961796722	49.211801783
H(10.25,t)	34.803247547	63.352693665	39.854139270	59.269416159	44.603100490
H(20.50,t)	25.918028306	72.753482081	31.491379443	67.136901218	36.012250860
H(30.75,t)	17.585926684	80.381591854	24.442660837	75.038739004	29.189659818

Table 8:
SSPRK(8,3) - First Example. n = 81, ∆t = 0.000406746.

	t = 0.4002	t = 0.6003	t = 1.2007	t = 1.4008	t = 1.9995
H(2.05,t)	46.399909872	53.286353811	47.231083889	52.665247115	49.914436402
H(10.25,t)	31.579783388	66.912987762	36.202602196	61.399626254	47.351201254
H(20.50,t)	19.995471282	81.770435903	28.348938355	67.897351090	41.470434276
H(30.75,t)	12.166379527	85.276889212	21.201382198	75.118232740	35.302671433

Table 9:
SSPRK(5,4) - First Example. n = 21, ∆t = 0.001626984.

	t = 0.4002	t = 0.6003	t = 1.2007	t = 1.4008	t = 1.9995
H(2.05,t)	46.622290500	52.787097816	47.859869461	52.001547901	49.400676570
H(10.25,t)	34.730774346	63.430241751	39.763815865	59.265371250	44.648718727
H(20.50,t)	25.902774758	72.900989346	31.627179567	67.085087086	35.949792153
H(30.75,t)	17.579794809	80.412269778	24.307952188	75.283937377	29.189659818

Table 10:
SSPRK(5,4) - First Example. n = 81, ∆t = 0.000406746.

	t = 0.4002	t = 0.6003	t = 1.2007	t = 1.4008	t = 1.9995
H(2.05,t)	46.392902495	53.248890600	47.286185082	52.662221091	49.958119629
H(10.25,t)	31.554043905	66.933341395	36.251208880	61.423049272	47.363253108
H(20.50,t)	19.996401835	81.768488200	28.400733481	67.866603775	41.527717899
H(30.75,t)	12.176423261	85.247735385	21.205533682	75.042390649	35.337746007

Table 11:
SSPRK(10,4) - First Example. n = 21, ∆t = 0.001626984.

	t = 0.4002	t = 0.6003	t = 1.2007	t = 1.4008	t = 1.9995
H(2.05,t)	46.630537958	52.874441670	47.782888375	51.946800931	49.325031239
H(10.25,t)	34.794634255	63.429623872	39.757019818	59.358695457	44.655587575
H(20.50,t)	26.040282068	72.960496815	31.666308626	67.200823483	36.204486777
H(30.75,t)	17.683457455	80.461548650	24.689195290	75.169506541	29.157117808

Table 12:
SSPRK(10,4) - First Example. n = 81, ∆t = 0.000406746

	t = 0.4002	t = 0.6003	t = 1.2007	t = 1.4008	t = 1.9995
H(2.05,t)	46.423072842	53.241602318	47.278751691	52.673792916	49.935388076
H(10.25,t)	31.564857169	66.932428451	36.247620786	61.417326849	47.399552813
H(20.50,t)	19.986709858	81.750450806	28.358903054	67.921668891	41.568215757
H(30.75,t)	12.173207506	85.290344991	21.236086361	75.100513047	35.308647206

Description of coeficients of each SSPRK schemes are presented in Tables 1, 4, 5 and 6.

We also test the SSPRK schemes with another example found in the literature. This one is described in ¹³. There, the authors considered a flow rate at the valve that depends also on the hidraulic head $Q(L,t)=\left(C_d{A}_v\right)H(L,t)$, then the boundary condition on Q at the end of the valve had to be calculated dinamicaly in both schemes. For MOC case, to calculate Q _n,j and H _n,j we have to so solve the system formed by the first equation of System 3.3 with i = n and $Q_{n,j}=\left(C_d{A}_v\right)H_{n,j}$. In SSPRK schemes, before calculating H _n,j , we just use formula $Q_{n,j}=\left(C_d{A}_v\right)H_{n,j}$. Figure 5 compares SSPRK(5, 4) with n = 21 with MOC. It was also collected some values of the numerical solution in ¹³ for comparison.

MOC scheme is more stable then SSPRK ones. For some values of ∆t stepsize, SSPRK numerical solutions presents spurious errors and asymptotic behavior. We test then ∆t stepsizes with four digit precision after comma to find, for each SSPRK schemes, the largest ∆t that makes numerical solution to not present asymptotic behavior. We use these values to compare SSPRK schemes and MOC schemes. For the MOC scheme, we used ∆t satisfying the relation $\frac{\Delta x}{\Delta t}=a$. In graphic of Figure 6 we show the relation between these ∆t stepsizes and the number of points n of the spacial discretization. In Graphic of Figure 7, we show the execution time. We use Scilab and a Lenovo ideapad notebook, Intel@ Core i5 − 8265U 1.60GHz ×8 (quad core, two threads each core). The operational system used was Linux ubuntu 19.10

Figure 6:
Comparison between the number of points n and corresponding ∆t for each scheme.

Figure 7:
Comparison between the number of points n and Execution Time, using ∆t of Figure 6.

7 CONCLUSION

The watter hammer problem is largely studied in engineering and mathematics. There are several numerical methods applied in it’s solution. In the past years, works involving applications on the problem usually uses system 1.1 and the methods of characteristics (MOC) to solve it. The MOC method shows great stability, considering different values for ∆x, furthermore, it also shows good execution time. One of the disavantages of the MOC method is that it creates a $\Delta x-\Delta t$ dependency that results on lack of choice in the discrete scheme, i.e, if we determine ∆x step, we also determine ∆t step. The goal of this work was to use strong stability preserving runge-kutta schemes (SSPRK) with a formulation of the usual method of lines on System 1.1. On autonomous systems, The SSPRK schemes acts to increase the order of the method maintaining the forward Euler stability. System 4.5 is not autonomous, so the order of the SSPRK schemes may be lower than specified, but still theoretically maintaining SSP stability for a ∆t choice satisfying $\Delta t\le C\overline{\Delta t}$, where $\overline{\Delta t}$ is an upper limit for stability to the forward euler method, and C is the SSP coefficient. Higher values of C, theoreticaly allows greater ∆t, decreasing then the execution time.

Based on Figures 3, 4 and 5, we conclude that all schemes shows similar hydraulic head amplitudes, i.e, they predict properly pressures the duct will suffer, but SSPRK schemes tend to show considerable difference in phase trough time. This difference seems to decrease for refined spatials discretization. In terms of amplitude prediction, SSPRK schemes also behave weel comparing to the experiment and MOC implementations presented in ²⁰ for the first example, and with MOC implementation presented in ¹³ for the second example.

Graphics of Figure 6, shows a gain in ∆t stepsize with the most part of the tested SSPRK schemes. This gain depends on the order of the SSPRK scheme and on the SSP coefficient C. Although allowing considerably increase of ∆t stepsize, MOC scheme still show less execution time as seen in graphic of Figure 7. To achieve greater ∆t, and consequently less execution time, we may increase the order or increase the constant C but this often increases also the stages necessary for the SSPRK methods, acting in the opposite direction. This is notably noted comparing SSPRK(5, 4) and SSPRK(10, 4). Although SSPRK(10, 4) admit greater ∆t, the difference in execution time is not so big.

To improve the execution time, we still may try other numerical methods apllied on System 4.5 or other SSPRK scheme. We also may consider multistep SSPRK methods ⁵, different space discretizations schemes ¹⁵ or even downwinding schemes ¹⁰.

REFERENCES

• ¹
  M.B. Abbott. “An Introduction to The Method of Characteristics”. Thames and Hudson London (1966).
• ²
  M.H. Chaudry. “Applied Hydraulic Transients”. Litton Educacional Publishing, New York (1979).
• ³
  Chi-Wang Shu & S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77(2) (1988), 439-471.
• ⁴
  J.A. Fox. “Hydraulic Analysis of Unsteady Flow in Pipe Networks”. The Macmillian Press LTD (1977).
• ⁵
  S. Gottlieb, D. Ketcheson & C.W. Shu. “Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations”. World Scientific (2011).
• ⁶
  S. Gottlieb , C. Shu & E. Tadmor. Strong Stability Preserving High Order Time Discretization Methods. SIAM Review, 43(1) (2001), 89-112.
• ⁷
  S. Gottlieb & C.W. Shu. Total Variation Diminishing Runge-Kutta Schemes. ICASE Report, 201591 (1996).
• ⁸
  C.C. Gumier. “Aplicação de Modelo Matemático de Simulação-Otimização na Gestão de Perda de água em Sistemas de Abastecimento”. Ph.D. thesis, Faculdade de Engenharia Civil, Arquitetura e Urbanismo, Universidade Estadual de Campinas, Campinas, SP (2005).
• ⁹
  Z. Horváth, H. Podhaisky & R. Weiner. Strong stability preserving explicit peer methods. Journal of Computational and Applied Mathematics, 296 (2016), 776-788.
• ¹⁰
  L. Isherwood, Z. Grant & S. Gottlieb . “Downwinding for Preserving Strong Stability in Explicit Integrating Factor Runge-Kutta Methods” (2018).
• ¹¹
  D. Ketcheson , S. Gottlieb & C. Macdonald. Strong Stability Preserving Two-step Runge-Kutta Methods. SIAM Journal on Numerical Analysis, 49 (2011).
• ¹²
  S.J. Ruuth. Global Optimization of Explicit Strong-Stability-Preserving Runge-Kutta Methods. Mathematics of Computation, 75(253) (2005), 183-207.
• ¹³
  M.D. Saikia & A.K. Sarma. Simulation of water hammer flows with unsteady friction factor. ARPN Journal of Engineering and Applied Sciences, 1(4) (2006), 776-788.
• ¹⁴
  D.F.G. Santiago, D.H. Paula, F.J. da Silva, D.B. Santos, E. da Conceição Silva, A.F. Antunis, N. da Silva & D.R. Trindade. Modelagem Matemática da Carga Hidráulica e Vazão em Condutos. ForScience, 7(1) (2019).
• ¹⁵
  C. Shu . “Differential Quadrature”. Springer-Verlag London Limited (1962).
• ¹⁶
  C. Shu . Total-Variation-Diminishing Time Discretizations. SIAM J. Sci. Statist. Comput, 9(6) (1988), 1073-1084.
• ¹⁷
  M.N. Spijker. Stepsize Conditions for General Monotonicity in Numerical Initial Value Problems. SIAM Journal on Numerical Analysis , 45(3) (2007), 1226-1245.
• ¹⁸
  R.J. Spiteri & S.J. Ruuth . A new Class of Optimal High-Order Strong-Stability-Preserving Time-Stepping Schemes,. SIAM J. Numer. Anal, 40(2) (2002), 469-491.
• ¹⁹
  V.L. Streeter. Unsteady Flow Calculations by Numerical Methods. Journal of Basic Engineering, (1972), 457-465.
• ²⁰
  R. Wichowski. Hydraulic Transients Analysis in Pipe Networks by the Method of Characteristics (MOC). Archives of Hydro-Engineering and Environmental Mechanics, 53(3) (2006), 267-291.

Publication Dates

• Publication in this collection
  04 Apr 2022
• Date of issue
  Jan-Mar 2022

History

• Received
  06 May 2020
• Accepted
  23 June 2021