Interpolation schemes for valve closure modelling

In transient flow condition the pressure variation will depend, among other factors, on the valve closure law, whose insertion into the water hammer’s software code is an impractical task since there are so many closing laws as valve types, forcing to modify the program whenever a specific valve has to be considered. It is more practical to calculate the valve closure (τ) at a certain time intervals and to transfer this information to the software’s data entry file. As the adopted temporal discretization for the valve closure’s law will not necessarily coincide with the simulation time step (Δt), the τ calculation at each Δt must be done by interpolation of two or more points belonging to the valve closure curve. The results will depend mainly of the valve type, the interpolation scheme used, and the interpolation order (IO) applied. In this article different interpolation methods are applied on two curve types: the first of linear type, characterized by having straight segments with abrupt slope changes; the second with softer forms. The results are compared with the exact solution. It is concluded that Newton-Gregory is the best interpolation method because it generates a lower computational cost and negligible interpolation errors regardless of the valve closure curve’s shape.


INTRODUCTION
The hydraulic transient analysis plays a very important role in the pressurized pipe systems design, and it is essential to the good operation of the pipe networks [12,20,31].The valve closure generates an unsteady flow that affects the pressure head.Of all unsteady flow situations, those caused by the valve's movement are the most common [17,30], and the gate's closure characteristic is one of the most important factors which is rarely investigated [13], more so when it is known that in a pipeline the pressure changes depend primarily on the water velocity, valve closure time, and valve closure curve's shape [11,21,27,29,32].In the valve both the pressure head and the fluid velocity change with the valve's opening change [5].This means that each valve adjustment sets up pressure pulse waves that traverse the system at the wave speed [33].In other words, any valve movement causes pressure waves to propagate through the system.The magnitude of the pressure waves depends on the type of valve and structure [29], the fluid velocity change (Δv) in time [22,31], the hydraulic system properties, and the elastic properties and restraints of the pipe system [17,30], all of which have an influence on the system's transient response [1].Karney and Ruus [11] recognize that during the valve closure the pressure head along the pipe rises and reaches a maximum, and that this maximum can occur during or at the end of the closure operation.The maximum pressure head magnitude and the instant when it occurs largely depend on the valve's openingtime relationship [30], being this a predominant factor in the pressure wave development that characterizes to the water hammer phenomenon [24].To compute the transient state conditions due to the valve closure, the operating rule or relative effective area of the valve opening (τ) in time may be specified either in a tabular form or by an algebraic expression [4].However, it is far more reassuring to be able to calculate the effects for a specific situation [17].Although the charts (tabulated data) or graphs of the τ versus t curves may be supplied by the manufacturers, in many cases the analytical expression or function for the valve closure is not available, and transient simulation must be solved considering a group of tabulated data.Taking into account this last fact, the transient analysis is forced to deal with interpolation techniques for the valve closure calculation in each time step (Δt) because the numerical interpolation is required every time an input data table is used.The conceptual basis of any numerical interpolation scheme is to fit some type of curve or function to a group of tabulated data.In waterhammer analysis, polynomial interpolation is very useful when only is available a data table which represents the real valve closure's curve.There are a large number of numerical techniques aimed at dealing with polynomial interpolation which are based on finding a polynomial that can engage to a certain number of known data points, as for example: Lagrange, Aitken, Newton-Gregory and Neville schemes.The application of these interpolation methods for the valve closure modelling is not abundant because generally is easier try to match the unknown curve with some other different curve previously known.In several cases it is easier to assume that the valve closure function has an approximately linear behaviour, leaving aside any complication associated with the non-linearity of the true curve.This conception of the problem can lead to results totally erroneous expressed as an over-estimation of the results or, in the worse case, it can lead to a sub-estimated solution, especially when a numerical interpolation technique is applied to a valve closure curve with sharp corners and abrupt changes in slope.Since the data points number minus one is referred as the interpolation order (IO), in several cases an IO's increment does not necessarily means an accuracy increment.Unless there is solid evidence that the interpolating function is close to the true function f (x), it is necessary to be cautious about higher IO [23].Some authors [23] have established that a high IO's scheme is more accurate when it is applied on a smooth function, obtaining a poor result when it is applied on a function with sharp corners and abrupt changes in slope.In the next paragraphs, polynomial interpolation schemes related with Lagrange, Aitken, Newton-Gregory and Neville will be shown and they will be applied in order to assess their main features and performance when the valve closure arrangement is represented in different forms: as a series of straight segments with abrupt changes in slope and as smoother continuous curves.A comparison between the results of each polynomial interpolation scheme and the exact result obtained through the analytical curves is shown.

INTERPOLATION PROBLEM STATEMENT
Suppose that Table 1 is generated from some experiment, where τ k , k = 0, …N, are known to be the values of a certain function τ(t), evaluated at τ k , k = 0, …, N in an interval containing these data points.
Table 1.Values of t i for the function τ i .
Note that only the functional values are known, but not the function τ(t).The problem is to find the functional value τ i corresponding to a non-tabulated intermediate value t=t i .Such a problem is called an interpolation problem [15].The numbers t 0 , t 1 ,…, t N are called nodes or interpolation data points [18].Thus, given: • The N+1 data points: The objective of the interpolation technique is predict τ i value at t = τ i .The interpolant or interpolating polynomial [18] for Lagrange (with an example of application below), Aitken, Newton-Gregory and Neville will be briefly described in the following paragraphs.

LAGRANGE INTERPOLATION SCHEME
Lagrange's method is a well-known classical technique for interpolation.Suppose the data set consists of N data points: (t 0 , τ 0 ), (t 1 , τ 1 ), …, (t N , τ N ).The interpolation polynomial τ(t) will have degree N−1.It is given by [3,6,7]: Where the functions φ i (t) i = 0,1, …, N are given by: The numerator of φ i (t) contains the entire sequence of factors (t -t 0 ), (t -t 1 )…(t -t N ), with the exception of the single factor (t -t i ).Likewise, the denominator contains the entire sequence of factors (t i -t 1 ), (t i -t 2 ), … (t i -t N ), with exception of the single factor (t i -t i ).
Notice that φ i (t j ) = 1 (numerator = denominator), but φ i (t j ) = 0 [numerator = 0, since it contains the factor (t j -t j ) for any j not equal to i].This means that P(t i ) = τ i , which is the exactly the wanted value.

Example of application
Consider a simple example, where the data for τ are (see Table 2): Table 2. τ vs. t data.
Where t = time (s) and τ = (K L0 /K L ) ½ , with τ = nondimensional loss coefficient which varies between 0 and 1, and K L = valve loss coefficient, where the subscript "0" in the numerator denotes "valve fully open" [30].In this case there are N = 4 data points, so we will create a polynomial of degree N -1 = 3. Then: t 0 = 0,t 1 = 10,t 2 = 20,t 3 = 30, and τ 0 = 1.0,τ 1 = 0.7,τ 2 = 0.3,τ 3 = 0.0.Besides: Multiplying each of these expressions by the corresponding τ i , i = 0, 1, 2, 3, and adding together the terms of like power then gives: A quick plot of the data together with the polynomial shows that it indeed passes through each of the data points (see Figure 1).For an intermediate τ value, for example, τ = 6 (s), the value for τ is equal to 0.843.
• Lagrange advantages: the interpolation polynomial can be written down without the solution of a linear system of equations; the basis polynomials φ i (t) depend only on t 0 , …, t N , and not on the values τ i .• Lagrange disadvantages: each τ(t) evaluation requires calculating N 2 additions and multiplications; adding a new data pair (t, τ) requires a new computation from scratch; and the computation is numerically unstable.
It is pertinent to mention that there are two more efficient and computationally attractive alternative schemes that help to mitigate the Lagrange disadvantages, especially in relation with its numerical instability: the Lagrange Modified Form (LMF) and the Lagrange Barycentric Form (LBF).
Because their level of complexity exceeds the objectives of this article, these methods will not be described here, although it is possible to find more information about them in Higham [9].

AITKEN INTERPOLATION SCHEME
The Aitken interpolation method is well-known and it is considered as the state-of-the-art for the interpolation functions over real numbers.Furthermore, Aitken interpolation method is constructive in a way that permits the addition of a new interpolation data point with low computational cost even in arbitrary space grids [19,25].The interpolation formula (1) needs all ordinates t i to find an estimated solution for the interpolation problem.For large problems this method is not efficient because it cannot make use of the previous interpolating result when we add more interpolating data points to get better result.To achieve a more accurate solution by adding more data points, Aitken uses the previous interpolating result to save multiplications.The method can be generally defined as follows [7,23]: ( ) In which i = 0, …, N and k = 0, …, N.This algorithm generates a sequence I i0, i1, …, iN which converges to a approximated function τ(t) for a given τ.
• Aitken advantages: it is easy to program; it reduces the general point interpolation to a succession of linear interpolations; it does not require the use of auxiliary interpolation tables; it is very useful for interpolating when a new node is added, with low computational cost.• Aitken disadvantage: it gives best results only when it is applied to evaluate the interpolator polynomial in a single point, or in a very small number of points.
Kumar [16] presents interpolation formulas in terms of divided differences through Aitken and Neville schemes which can be used to iterate even Hermite interpolation problems.The numerical examples suggest that the Neville divided difference form is a better choice than Newton and Aitken's divided difference forms even in Chebyshev points.Due to space reasons no further details will be given here.

NEWTON-GREGORY INTERPOLATION SCHEME
A polynomial interpolation is a Newton-Gregory type if the data points t 1 , t 2 , t 3 , … are equispaced.The Newton-Gregory's method is based upon the finite differences method and it can build up what is called a difference table for a data point series.Rather than defining a functions set for a given data points set as Lagrange did [the set of Li(t) functions], Newton-Gregory takes things one node at a time.It is an interpolation method based upon the Taylor series.The unique polynomial satisfying some N + 1 value is of degree N [3,[6][7][8]: The N+1 constants a i being found from the N+1 linear equations obtained by data substitution.Suppose the N+1 values are (t i , τ i ), i = 0, 1, 2,…, N, and suppose that equal intervals between successive values is h, then t r -t s = (r-s) • h.Forming successive differences and taking into account that t-t 0 = kh, with 0 ≤ k ≤ N, then: The expression ( 11) is known as the forward interpolation formula and it is appropriate when the required value (t 0 , τ 0 ) lies near the beginning of the tabulated data; otherwise, when the value (t 0 , τ 0 ) is near the end of the table, it is necessary to apply the backward interpolation formula.
• Newton-Gregory advantages: the addition of more interpolation points is enough to calculate the new coefficients to find the polynomial; it is possible to know the maximum polynomial degree for a specific precision; and the errors are easy to detect using the coefficient table.• Newton-Gregory disadvantages: the polynomial does not adapt to an inverse interpolation, unless it is linear; if there are more than one ordinate values set, new coefficients must be calculated for each one.

NEVILLE INTERPOLATION SCHEME
Let P 1 be the value at t of the unique zero degree polynomial passing through the node (s 1 , τ 1 ); so P 1 = τ 1 .Likewise define P 2 , P 3 , P 4 , … P N .Now let P 12 be the value at t of the unique one degree polynomial passing through both (s 1 , τ 1 ) and (s 2 , τ 2 ) values.Similarly define P 23 , P 34 , …, P (N-1)N for the higher-order polynomials up to P 123…N , which is the value of the unique interpolating polynomial through all N data points (i.e., it is the desired answer).The various P form a tableau with ancestors on the left leading to a single descendant at the extreme right.The Neville's algorithm is a recursive method which is based on the relationship between a daughter P and its two parents [3,23]: This recurrence works because the two parents already agree at data points s k+1, …, s k+m-1 .An improvement on the recurrence ( 12) is to keep track of the small differences between parents and daughters through the definition of corrections factors C and D. At each level m, the C's and D's are the corrections that make the interpolation one order higher.The final answer P 1, …, N is equal to the sum of any τ i plus a set of C's and/or D's that form a path through the family tree to the rightmost daughter.Because the advantages and disadvantages of Neville are similar to those of Aitken, so do not give further details here.

EXAMPLE OF APPLICATION 1: VALVE CLOSURE CURVE WITH STRAIGHT SEGMENTS AND ABRUPT CHANGES IN SLOPE
Each interpolation scheme mentioned above will be applied for the interpolation of the curve that describes the closure of a control valve (see Figure 2).The valve closure decreases linearly from τ = 0.6 to τ = 0.2 in 10 (s).This valve setting is maintained for 15 (s).The valve then reopens linearly to its initial value of τ = 0.6 in 5 (s).The valve closure arrangement has 7 data points.Figures 3 to 8 show the numerical performance of each interpolation method.Table 3 shows the average error generated by each interpolation method according to the applied IO.
In this case the average error E(t) is calculated as follows: Where T max = maximum simulation time = 30 (s), τ(t) is the exact value obtained from τ vs t equations shown in Figure 2, Δt = 0.5 (s), and p(t) is the interpolated value taking into account in the calculation 2, 4 and 7 data points of the valve's closure curve.Table 3 and Figures 3 to 8 show that Newton-Gregory interpolation method is the only one that presents a coincident result with the exact solution (average error = 0%) when IO = 1.The remaining interpolation schemes have different error ranges depending on the used IO, revealing that the Lagrange shows the poorest performance when IO is equal to 1 and 3 (Figures 5 and 6).On the other hand, Aitken, Neville and Newton-Gregory schemes register the same error when IO = 6.Except for Newton-Gregory, all the schemes have the smallest error when IO = 6.

APPLICATION EXAMPLE 2: VALVE CLOSURE CURVE WITH SOFT SHAPE
In this case the following equation will be applied [2]: Figure 2. Valve closure arrangement [10].
Table 3. Errors when the curve of Figure 2 is interpolated using different interpolation schemes with IO = 1, 3 and 6.

Interpolation Method
Average Error E(t) (%) Newton-Gregory 0.0 −4.Where τ(t) is the exact value for the valve closure obtained from equation ( 14); t = time (s); T = valve's closure time = 5 (s), and n is an exponential parameter whose value depends on the valve's type (see Figure 9 and Table 4).In this case it is assumed that valve's discharge coefficient does not vary over time.The τ(t) value from equation ( 14) will be calculated for different valve's types, every each of them with a different closure law composed by 11 data points.Then Lagrange, Aitken, Newton-Gregory and Neville interpolation schemes will be applied using every algebraic curve with IO equal to 1, 2, 3 and 10.The main goal is to determine which combination between interpolation scheme and IO is better suited for each valve closure curve (ball, butterfly, etc.).In order to calculate the errors, the analytical solution given by equation ( 14) will be considered as the exact solution.
For example, the Figures 10, 11, 12 and 13 show the error generated by Lagrange interpolation scheme when it is applied to calculate τ from equation (14) with IOs equal to 1, 2, 3 and 10, respectively, taking into account that the superscript n is valve's type dependant, such as is shown in Table 4. From Figures 10 to 13 it is possible to see that Lagrange scheme's results are mainly dependant of: • Type of valve.
Table 5 shows a results summary, where the average error E(t) is calculated using the equation ( 13), with T max = 5 (s) and Δt = 0.5 (s) in this case.
In Table 5 Lagrange always presents negative values when it is applied to interpolate the ball or butterfly closure valve curves (Figures 10 and 11), with extreme errors ranging between -592% (ball valve) and -2783% (butterfly valve), in both cases at time 4.5 (s).Furthermore, in the case of gate (square) and needle valves, Lagrange always registers positive average errors (Figures 12 and 13).In the first case, the maximum positive error is +75%; in the second case, the maximum positive error is +76%, in both cases at 4.5 (s).An interesting point concerns to the IO.For example, in all the cases where Lagrange was applied, the best result was obtained with IO Table 4. Values for n according to the valve's type [2].

COMPUTATIONAL COST
An important interpolation's issue is referred to the computational cost (or number of algebraic operations) that each scheme must perform for obtaining some result, which directly depends on the number of data points (N) used in the interpolation process.For example, considering only the multiplications and divisions (MDs)'s number, the algebraic operations quantity used by each method (in each time step) is shown in Table 6.Table 6.Number of MDs according to interpolation scheme [14,28].shows algebraic operations number for Lagrange, Aitken, Neville and Newton-Gregory when N = 2, 3, 4 and 11 in each case, where it is verified that Newton-Gregory is the most efficient method no matter the chosen N. Aitken and Neville are more efficient than Lagrange when N≤5.When N = 11, Newton-Gregory calculates 77 MDs at each time step.In contrast, Lagrange and Aitken / Neville must calculate 168 and 198 MDs, respectively, which mean an increase of 118% and 157% in comparison with Newton-Gregory, respectively.

CONCLUSIONS
The valve closure curve interpolation is not a trivial task because the results depend mainly on the interpolation method adopted and the chosen IO, being difficult to choose the best method / IO combination without prior careful analysis to obtain results with a minimum error in the interpolation process.The way the value of τ(t) varies in time differs from one type of valve to another, so it is not always possible to have a generic interpolation method that shows stability and accuracy for all cases.The details become more relevant depending of the closure curve's form, what can condition the relationship between pressure and flow [26] and strongly influence the transient conditions of the pipe, such as was discovered by Joukowsky at the end of the century XIX.For these reasons it is important to calculate (or rather interpolate) the τ(t) value in every computational time step as accurately as possible.By observing the results shown in Table 3 it is concluded that when the valve opening curve has straight segments and abrupt changes in slope as is shown in the Figure 2, the best interpolation method is Newton-Gregory, which it should always be applied with IO = 1.On the other hand, when the valve opening curve has smoother forms (as is shown in Figure 9), the best option remains Newton-Gregory regardless both the IO and the valve's type chosen (ball, butterfly, etc.).In this case Lagrange registers the poorest results, especially when IO= 1 (see Table 5).Besides, Lagrange, Aitken and Neville register a zero average error only when IO = 10.This last point is relevant because according to Figure 20, as N grows, the number of MDs increases in each Δt.Therefore, the adoption of a greater IO necessary to guarantee the Lagrange, Aitken or Neville's numerical accuracy, unfailingly increases these methods' inefficiency level in computational cost terms.Finally, when the valve closure curve is interpolated, the Δt magnitude must fulfil with two requirements at least: (1) it must be a multiple of the valve's closure time (T), so that the calculation allows incorporate exactly the time in which the last τ vs. t curve's node is calculated; (2) its value must allow calculate the τ corresponding to the special nodes.For example, the curve of Figure 2 was discretized with Δt = 0.5 (s), multiple value of T = 30 (s).With this the nodes located at t = 0, 10, 25 and 30 (s) were included into the τ calculation.
Obviously, these conditions are difficult to achieve when the transient flow is modelling in complex pipe networks, where Δt should be calculated not only according to the τ vs. t curve shape, but also according to the pipe network's characteristics (lengths, wave speeds, system discretization, etc.), which may force to alter some initial values in order to fulfil with the Courant condition.

Figure 19 .
Figure 19.Results obtained by Neville according to different valve's type.IO = 10.

2 Figure 20 .
Figure 20.Number of MDs of every interpolation scheme according to adopted N.