Comparison of viscoelastic models with a different number of parameters for transient simulations

The numerical models for the simulation of transients are an important tool for the design and the diagnosis of pressurized pipe systems (Covas et al. 2004; Duan et al. 2010b, 2012; Meniconi et al. 2011; Gong et al. 2013; Massari et al. 2014; Duan 2017; Laucelli et al. 2016; Sun et al. 2016). Such models require the definition of a functional relationship between stresses and strains in the pipe material, i.e., its rheological model. The spreading of pipes of polymeric materials, such as high density polyethylene (HDPE) and polyvinyl chloride (PVC), and their use in pressurized pipe systems instead of metallic pipes, increased the complexity of the rheological models, from elastic to viscoelastic. The increase in the complexity of the models corresponded to an increase in the number of the parameters they depend on and to the need of reliable and efficient algorithms to calibrate them.

For the specific purpose of transient modeling in polymeric pipes, the evaluation on specimens of the rheological parameters, i.e., by creep and relaxation tests in laboratory, is not successful (Ghilardi & Paoletti 1986; Covas et al. 2005). For this reason, the most common procedure used so far is the calibration of the parameters of the implemented rheological models using transient tests.

In the literature, several case studies and calibration procedures have been presented. As an example, Pezzinga & Scandura (1995) used a trial and error procedure to calibrate a rheological model consisting of one Kelvin–Voigt (KV) element with two parameters, by means of transient tests on a HDPE pipe. In the following papers, Pezzinga (2000, 2002) and Pezzinga et al. (2014) calibrated the same rheological model by means of a micro-genetic algorithm. The calibration procedure is based on the maximization of a function evaluated as the inverse of the sum of the squared residuals (RSS) between simulated and experimental data. The numerical model is obtained by the integration of the governing equations with the method of the characteristics (MOC).

Covas et al. (2005) used the experimental data collected from a PE pipe-rig to calibrate a rheological model composed by a different number of KV elements, fixing some of the parameters and varying the others. The calibration is based on the minimization of the RSS by means of a Levenberg–Marquardt algorithm. A sensitivity analysis on the RSS suggests that the optimal results are obtained with five KV elements and five calibration parameters. The need of a proper choice of the time step in the MOC model is discussed.

Soares et al. (2011) used a genetic algorithm to calibrate different rheological models on data coming from tests carried out in an experimental hydraulic circuit composed of PVC pipes. Numerical MOC models implementing different rheological models, with up to three KV elements, are calibrated minimizing the RSS by means of both a genetic and a Levenberg–Marquardt algorithm. The algorithms are used in series since the genetic algorithm is considered as a good tool for the detection of local minima but not accurate enough for the global minimum determination. The considered optimal value is obtained by one KV element and fixing some parameters.

Although in the available literature different choices are made with regard to the complexity of the model, the number of the involved parameter is not considered as a parameter itself. In other words, the complexity level of the used viscoelastic models is not investigated as a variable of the problem and determined by means of an objective and rational criterion.

In this paper, several aspects of the calibration of the viscoelastic parameters are investigated. In the first part, the implementation of the viscoelastic elements in a frequency domain model is described. The choice of the optimization function and algorithm is discussed, giving some general remarks regarding the calibration procedure. A method for the estimation of the appropriate number of the viscoelastic model parameters is then presented. The proposed approach is applied to synthetic data, obtained corrupting with a white noise the numerical results of an implemented viscoelastic model, and then to experimental data, obtained by tests carried out on a HDPE pipe. The obtained results and the original contributions of the paper are summarized in the conclusions.

Many authors (e.g., Ghilardi & Paoletti 1986, 1987; Pezzinga & Scandura 1995; Covas et al. 2005; Lee et al. 2009; Duan et al. 2012; Keramat et al. 2012; Pezzinga 2014; Evangelista et al. 2015) have presented the implementation of S rheological components in MOC numerical models, to relate hoop stresses and strains in the pipe. Other authors have presented the implementation of the same rheological components in frequency domain integration models (e.g., Duan et al. 2010a, 2012; Vítkovský et al. 2011; Lee et al. 2013; Capponi et al. 2017a; Ferrante & Capponi 2017b). With respect to the time domain, the frequency domain integration is faster and easy when the variation in time at one section of the system is sought. The frequency domain integration requires a linearization of the governing equations but limited to the second order of the steady-friction resistance term, when the boundary conditions are linear (Capponi et al. 2017c). Since the viscoelastic term in the governing equations can be represented as a convolution in time (Weinerowska-Bords 2015) it does not require a linearization and, due to the well-known properties of Fourier and Laplace transforms of such integrals, it is easier to be handled in the frequency domain than in the time domain. Similarly, the unsteady-friction term is a convolution integral in the time domain as well (Weinerowska-Bords 2015) and can be considered and easily implemented in the frequency domain models. With reference to the considered conditions and to the results presented in the following, the effects of viscoelasticity prevail on those introduced by the friction term (Duan et al. 2010a, 2012). Hence, for the sake of simplicity and considering that the main purpose of this paper is the comparison of the viscoelastic models, a frequency domain numerical model without an unsteady-friction term is used.

The function depends on the chosen rheological model of the pipe material (Ferrante & Capponi 2017b). In the fifth column of Table 1 the formulae are given for the S models, for different values of m. In these formulae, assuming that the thickness of the pipe wall ⁠, and ⁠, it is ⁠, where is the pipe constraint coefficient and is the water density. In the same table the model parameters are also specified. Since all the used models require a spring in series with the rest of the element arrangement, the Young modulus of this spring, ⁠, should be added directly to the list or indirectly, i.e., considering the wave speed a as a further parameter.

In Equation (6), is introduced as a parameter, without any need of further space discretization and compatibility conditions between space and time grids. If a reduced number of measurement sections is used in the simulation, the large amount of calculations needed by the time domain integration, at each time step and for each node of the grid, can be strongly reduced by the frequency domain integration. Since for the calibration procedures the dependent variables are evaluated at a few sections, frequency domain models are much more efficient than their time domain correspondent. The reliability and the accuracy of frequency domain solutions rely on other parameters (Lee & Vítkovský 2010; Lee 2013).

With respect to the approximations introduced by the steady-friction term linearization, Capponi et al. (2017c) have shown that, with reference to the test conditions of this investigation, their effects are negligible when compared to those due to viscoelasticity. Furthermore, the computational burden is reduced of about two orders of magnitude, with respect to the time domain (e.g., Capponi et al. 2017a, 2017b; Ferrante & Capponi 2017b). Hence, considering that the increase in the simulation efficiency justifies the introduced approximations, the frequency domain model is used in the following.

The values of the parameters of the numerical model can be evaluated by means of a calibration procedure, i.e., by the comparison of the results of the numerical model with the experimental data. For a given model, the values of the parameters are varied to minimize the distance between the model and the actual phenomenon that produced the experimental data. Three different aspects can be considered: (i) the definition of a measure for the distance between the chosen model and the actual phenomenon, (ii) the algorithm to minimize such a distance, and (iii) the definition of a measure to compare the distance from the actual phenomenon of models with a different number of parameters.

Several functions can be used to measure the distance between the actual phenomenon and the model, or at least the distance between the acquired pressure head values during a transient, and the corresponding simulated values, H, produced by a model (Savic et al. 2009). The calibration procedures aim to estimate the parameter values that minimize or maximize the value of this function, defined as the optimization function.

Various algorithms can be used to minimize the optimization function in the parameter space. The optimization by trial and error is not a reliable and reproducible process while the optimization by brute force, i.e., by the direct scrutiny of the function on a grid of the parameter space, is numerically cumbersome and can be used only for efficient models with a reduced number of parameters. A genetic algorithm and two other nonlinear minimization techniques, i.e., the Nelder–Mead or the Levenberg–Marquardt algorithms were used for the presented calibration procedure. All the mentioned algorithms are available as libraries of interpreters and compilers. A freeware collection of minimization functions is provided in Octave (Eaton et al. 2015).

The minimization of the optimization function leads the chosen model as close as possible to the actual phenomenon but it does not provide enough information about the relative distance of different models from the actual phenomenon, at least when the comparison is among calibrated models with a different number of parameters. In this case, the measure of the relative distance cannot be solely based on the residual analysis and on RSS but needs also to consider the number of the parameters involved.

The model assessed to be the best among the others in a set is the one with the minimum value of ⁠, i.e., ⁠. For all the other fitted models, the larger the value of ⁠, the less plausible it is that they can be considered as the best model. As a rule of thumb (Burnham & Anderson 2003), all the models with have still substantial support to be estimated as the best, models with have considerably less support while models with have no essential support.

It is worth noting that Equation (12) depends on the value of ⁠. If the parameter estimation procedure is carried out perfectly, the obtained value of for a chosen model is the minimum. Hence, the efficiency and reliability of the calibration algorithm in reaching the minimum of the objective function plays a crucial role in the comparison of the models. This aspect is of great importance for the fitting of the viscoelastic models in the case of diagnosis by means of transients, because of the low gradients of the optimization functions in the parameter space region close to the minimum. This flat region implies that the sensitivities of the optimization function with respect to the parameters are low, and hence, that the minimization algorithm has to be efficient since large variations in the parameters correspond to small variations in the optimization function values.

To verify the reliability of the calibration procedure and of the proposed criterion of model comparison, a synthetic data set was generated by means of a numerical model based on Equation (6). Two RPV systems (Figure 1), with pipes of a total length, ⁠, of 102.58 m, and diameter, D of 96.8 mm, were considered, made of two different pipe materials, similar to HDPE. For the first material, A, the rheology was considered completely defined by an S3 model with = 353.73 m/s, = 4.500 10⁹ Pa, = 4.500 10¹⁰ Pa, = 4.500 10¹¹ Pa, and = 9.000 10⁹ Pa s. As a result, the KV elements were characterized by = 0.020 s, = 0.200 s, and = 2.000 s. The second material, B, was still defined by an S3 model and had the same value of a as the first one, but = 4.500 10¹⁰ and = 9.00 10⁸ Pa s, = 9.00 10⁹ Pa s, and = 9.000 10¹⁰ Pa s. As a result, the KV elements were characterized by the same values of of the previous model, i.e., = 0.020 s, = 0.200 s, and = 2.000 s.

The steady-friction term was evaluated during the simulation with the Blasius formula for turbulent flow in smooth pipes. The pressure head variations in time, or pressure signals, generated in the two systems with the two different pipe materials, were then corrupted with a white noise. The same generated noise was used for both signals, with a zero mean and a standard deviation of 10⁻¹ m, to simulate the acquisition by an actual pressure transducer. The sample variance of the generated random noise was 1.005 10⁻² m².

Assuming that the generated synthetic signal, affected by the superimposed random noise, represents the truth and that the model based on the implemented S3 contains all and only the actual governing equations of the phenomenon, a blind calibration procedure was applied as in cases where the actual model is unknown. Since, in reality, also the number of the parameters is considered as unknown, models with a different number of KV elements and parameters were considered.

The calibration procedure for the five models (S1, S2, S3, S4, and S5) and the two pipe materials (A and B), required the use of a nonlinear algorithm to minimize ⁠. The first 20 seconds of the pressure signal were considered for the evaluation of ⁠. The upper and lower bounds of 10¹² and 10⁶ were used for the viscoelastic parameters. A genetic algorithm, with populations of 200 individuals, a number of populations equal to the number of the parameters and up to 100 generations, produced values of always greater than 3.0, far from the known solution. As a consequence, the heuristic Nelder–Mead algorithm was used instead, with a minimum tolerance of 10⁻¹⁶ on the function and on the parameter variations and a maximum number of 1,000 function evaluations. Since references to the use of the Levenberg–Marquardt algorithm can be found in the literature (Covas et al. 2005; Soares et al. 2011), a check was made to verify that the Nelder–Mead solution was the same as the Levenberg–Marquardt algorithm, with the same tolerances. At the end of the calibration procedure, the optimal set of parameters and the minimized value of were available for each combination of model and material.

In Figure 2 the synthetic signals are compared to those simulated implementing the viscoelastic models S1, S2, S3, S4, and S5 and considering the system made of the two materials, A and B (Figure 2(a) and 2(b), respectively). Only small differences can be appreciated between the synthetic signals and those simulated by S1 and S2, while other simulated signals overlap for both materials. Differences among results of the models S1, S2, and S3 are much more evident in Figure 3, where the absolute values of the functions of Equation (5) are plotted for the pipe materials A and B (Figure 3(a) and 3(b), respectively). The results of the models S4 and S5 were very close to those of S3 and are not plotted.

The wave speed was included as a parameter in the calibration procedure. The calibrated values of a were between 353.2 and 354.1 m/s, and hence with differences of less than ±0.5 m/s with respect to the value, used to generate the synthetic data.

The minima reached for decrease with the number of KV elements, m, and tend towards the sample variance of the added random noise.

Figure 4 presents the calibration results in the parameter space. Due to their significance, J and are chosen instead of E and ⁠. Dashed lines and solid circles denote the true solutions.

With reference to the values, we can conclude that, for both materials, values of S1 and S2 fall in the range of the variation of the true values. For the model S3, the values corresponding to the obtained minimum are very close to the true ones in the case of material A (Figure 4(a)) while for material B both and are clearly far from the true ones (Figure 4(b)). For S4 and S5, the missing solutions overlap with others or are outside the limits of logarithmic axis.

No correlation can be clearly indicated between the calibrated values of and the maneuver duration (about 0.15 s), the pipe characteristic time (0.58 s), and the duration of the synthetic signals used for the calibration (20 s).

Figure 4(a) shows an interesting correlation between estimated and J values, which leads the calibration to end along the line defined by the true solutions in the parameter space.

The variation of with the number m of the KV elements of the used models is shown in Figure 5. While is monotonically decreasing with m, the minima of are obtained for ⁠, corresponding to the value of the true model for both materials, A and B. The models S1 and S2 are characterized by a and hence can reasonably be excluded if compared to S3. The obtained values of ⁠, close to 4 and 7 for models S4 and S5, respectively, suggest that these models are significantly more distant from the true model than S3.

The crosses in Figure 5 denote the theoretical values obtained by assuming = 1.005 10⁻² m in Equation (12), i.e., considering that the differences between calibrated and actual data are completely explained by the random noise. A calibration that does not reach the theoretical minimum produces a higher value of and affects the measure of the distance from the true model. In this sense, a reliable calibration procedure is crucial. Since hollow circles and crosses coincide in Figure 5, the used algorithms, with the given tolerances and settings, can be considered sufficiently accurate for the aims of the comparison.

To reduce the number of the parameters, some authors (Covas et al. 2005; Soares et al. 2011) prefer to fix the values of and calibrate the others. To assess the use of this strategy, five reduced S models, with fixed values of ⁠, were also considered in the set of model candidates, with m ranging from 1 to 5. The fixed values of were the first m values of 0.01, 0.10, 1.00, 10.00, and 100.0 s. In this case, the obtained value of for the reduced S2 model, with fixed values of = 0.01 s and = 0.1 s, was less than that obtained by the best calibration of the S3 model, estimating all the parameters. This result suggests that, following this strategy, the number of parameters could be further decreased.

The calibration procedure based on the Nelder–Mead algorithm and the model comparison methodology verified on the synthetic data were applied to the experimental data to determine the optimal number of the KV elements needed to model transients in an actual HDPE pipe.

The tests were carried out at the Water Engineering Laboratory of the University of Perugia, Italy, on the same RPV system used to generate the synthetic data, with the same and D. The pipe was a HDPE DN110 PN10, according to UNI EN 12201 and UNI EN ISO 15494. An electromagnetic flowmeter was used to measure the discharge in the steady-state initial conditions, with an accuracy of 0.2% of the measured value. A piezo resistive pressure transducer, PT⁠, with a full scale of 6 bar and an accuracy of 0.15% of the full scale was used to measure the pressure upstream of the maneuver valve.

At the downstream end of the pipe, a hand-operated ball valve, DV, discharging into the air and a remotely controlled butterfly valve, MV, immediately upstream of DV (Figure 1) were used to have an initial condition of = 3.64 l/s and to produce the maneuver of Equation (13), with the same parameters used for the generation of the synthetic data. The obtained variation in time of the pressure head at PT⁠, ⁠, is shown in Figure 6.

Six models differing for the number m of the KV elements were implemented and calibrated to the experimental data. The Nelder–Mead algorithm was used to evaluate the parameters of the models minimizing the value of ⁠, with the same tolerances and limits used for the synthetic data. The comparison of the simulated and experimental data in Figure 6 and in Figure 7 confirms the results obtained for the calibrations by synthetic data. While the differences in the pressure signal are small, differences in the function are evident, at least for models S1, S2, and S3. Results of the S4, S5, and S6 models, not shown in the figures, are very close to the results of the S3 model.

The values of the parameters obtained for the six models by the calibration procedure are shown in Figure 8.

As expected, the values of obtained at the end of the calibration for each model decrease monotonically with the number of used KV elements.

The values of shown in Figure 9 allow assessment of whether the decrease is justified by the increase in the model complexity. Based on ⁠, S3 is the model closest to the actual phenomenon that generated the experimental data, among the set of models containing S1, S2, S3, S4, S5, and S6. No other model in the set shows a and hence there is no substantial support to estimate another model in the set as the best one.

Assuming that the value of of the model S3 corresponds to that of the noise of the experimental data, the values of can be estimated for S4, S5, and S6. The comparison of these values (crosses) with the data obtained by the calibration confirms that the increase of the number of the parameters, for ⁠, cannot be justified by the decrease in ⁠. On the other hand, the decrease in from S1 and S2 to S3 justifies the increase in the number of the parameters.

The results shown allow some remarks on the calibration of viscoelastic models for transient test in polymeric pipes. Several issues need to be considered to achieve reliable results.

The first issue is related to the choice of the viscoelastic model. Although other models could be considered (Ferrante & Capponi 2017b), the most common choice is a series of one or more KV elements with a spring. A set of viscoelastic models with up to six KV elements is considered in this paper.

The parameters of the models can be estimated by means of transient tests, since creep and relaxation tests on specimens do not provide reliable results. This leads to the second issue, that is to define a measure of the estimation reliability of a set of model parameters. This measure is usually expressed as a function of the differences between simulated and experimental data. The estimate of the parameters requires the minimization of this function. The sum of the squares of the residuals is one of the most common choices and is used in this paper. Unfortunately, this function presents wide flat regions in the parameter space close to the optimal solutions and this leads to the third issue, that is the choice of the algorithm to be used to estimate the parameters minimizing the distance between simulated and measured data. The most used algorithms in the literature are also applied in this paper to synthetic and experimental data. The genetic algorithm with the shown settings is confirmed to be not reliable and efficient in reaching the global optimal solution. In contrast, the Nelder–Mead algorithm is able to attain the optimal solution in a reasonable computational time, even for low tolerances. This algorithm is simpler and more robust than the Levenberg–Marquardt algorithm, at least for the considered cases.

All the models in a set, with a different number of viscoelastic elements, provide optimal solutions and this leads to the fourth issue, that is, the comparison of models with a different number of parameters. The need of the trade-off between the reduction of the squared sum of the residuals and the increase in the number of the parameters can be handled using different approaches. The Akaike information criterion is used in this paper to introduce into the calibration procedure the number of the parameters as a further parameter, although other techniques can also be used for an efficient investigation of the behavioral model space (Beaumont et al. 2002; Sadegh & Vrugt 2014). The use of this criterion is tested on synthetic data and is applied to estimate the best model for experimental data. The obtained results show that the criterion can be successful and that a viscoelastic model with three KV elements, S3, is the best estimated model in the set of models with up to six KV elements.

The calibrations presented in this paper for the experimental data did not show any correlation of with the duration of the maneuver, the duration of the acquisition, and the characteristic time of the system. For the synthetic data, these durations are taken into account in other model components and are not correlated to the parameters in the used viscoelastic model equations. The difficulties in the calibration procedure suggest that the sensitivity of the sum of the squared residuals to is low. Variations of one order of magnitude of can cause small variations in the sum of the squared residuals. Probably due to this low sensitivity, the strategy of a reduction of the number of parameters by fixing values of can be successful, as for the shown application of a reduced S2 model. The reduction of the number of the parameters by fixing the values with a different order of magnitude (e.g., 0.01, 0.10, 1.00, 10.00 s) is based on a subjective choice of the fixed values and unfortunately it does not allow a general conclusion to be drawn here.

With reference to the applications, different model uses, such as for the preliminary design of a simple pipe or for the diagnosis of a complex system, require different levels of accuracy. In this paper, the obtained minima are comparable to the accuracy of a pressure transducer, as required by the diagnosis of pressurized pipe systems by transient tests. Furthermore, the initial 20 s of the experimental pressure signal are used to evaluate ⁠. This duration should also be changed according to the aims of the calibration and of the model use, since the design of a simple system could require the simulation of the first characteristic time, while the long-term analysis could require longer durations, with respect to that used here.

The implementation of the rheology of a polymeric material pipe in the transient governing equations requires the estimation of the viscoelastic parameters. The most common way to define such parameters is by means of a calibration procedure on experimental data. The main contribution of this paper is to address the problem of overfitting of viscoelastic models for transient analysis.

The implementation of the viscoelastic component in a frequency domain model reduces the computational time and speeds up the calibration. Furthermore, the presence of the viscoelastic parameters in the governing equations is limited to a complex function, ⁠. The analysis of this function for the considered calibrations shows that it could be used to enhance the differences from one model to another.

The standard procedure used, based on the minimization of the sum of the squares of the residuals by means of the combined use of a genetic and the Nelder–Mead algorithm, allows it to be shown that the Akaike criterion can be used to define the optimal number of parameters.

Since all the used algorithms and methodologies are easily available, this criterion, or other available criteria, can be easily used to avoid the overfitting. Additional investigations are needed to confirm some of the results shown and to asses the proper choice regarding different or more sophisticated viscoelastic models, optimization functions, calibration algorithms, and criteria for the comparison of different models.