Mesh Sensitivity Analysis of Axisymmetric Models for Smooth–Turbulent Transient Flows

Abstract

The current paper focuses on the assessment of radial mesh influence on the description of the transient event obtained by an axisymmetric model. The objective is to reduce computational effort while accurately calculating hydraulic transients in smooth–turbulent pressurized pipes. The analyzed pipe system has a reservoir–pipe–valve configuration with an inner diameter of 0.02 m and a total length of 14.96 m, with the initial discharge being equal to 120 × 10⁻³ L/s (Re = 7638). An extensive study is carried out with 80 geometric sequence meshes by varying the total number of cylinders, the geometric common ratio, and the pipe axial discretization. The benefit of increasing the geometric common ratio is highlighted. A detailed comparison between two meshes is presented, in which the best mesh (i.e., the one with the lowest computational effort) has a three-fold higher value of the geometric common ratio. The two meshes show small differences for the instantaneous valve closure, limited to a time interval immediately after the arrival of the pressure surge and only during the first pressure wave. The dynamic characterization of the transient phenomenon demonstrates the in-depth consistency between the model results and the hydraulic transients’ phenomenon in terms of the piezometric head, the wall shear stress, and the mean velocity time-history, in comparison to the results obtained with the shear stress, lateral velocity, and axial velocity profiles.

1. Introduction

Transient flows occur when pressure and flow rate changes are created by manoeuvres in valves or in turbomachines. The numerical modelling of these flows is a complex engineering task that should ideally consider secondary phenomena that are not described by the Joukowsky formulation, such as, fluid–structure interaction, cavitation, frequency-dependent energy dissipation, and pipe wall viscoelasticity [1,2].

A hydraulic transient is not a one-dimensional (1D) phenomenon and has complex radial velocity profile dynamics. It has a different behaviour between the pipe core and the wall region, with the latter having higher velocity gradients. Hence, the unsteady energy dissipation is not in phase with the mean velocity. Therefore, there is no exact and completely satisfying solution for the calculation of unsteady friction in fast transient flows in 1D models [3,4]. The 1D instantaneous acceleration-based formulations lack generality and accuracy [5,6,7,8,9]. The 1D convolution integral-based models require a significant computational effort, and the simplified versions lack both accuracy and generality [10,11,12]. Conversely, the three-dimensional computational fluid dynamics (3D CFD) models are highly accurate in portraying the 3D nature of the transient flow. However, these models are extremely time-consuming for engineering practice [13]. Moreover, the inconsistencies between the outcomes of the various turbulence models raise questions about the rationale for their application [11].

The axisymmetric models, hereafter referred as quasi-2D (Q2D), have proven to be consistent with the phenomenon dynamics, leading to highly accurate results [14]. Their main advantage is to describe solely the most essential features of the transient. The uniform cross-sectional pressure and the flow radial symmetry assumptions ensure the “simplicity” of these mathematical models. These simplifications eliminate the radial and circumferential momentum equations and neglect the circumferential velocity component in the continuity equation [15,16]. Additionally, turbulent models can be adapted for axisymmetric flows, extending beyond the frozen-viscosity approximation, as is the case with 1D convolutional models [17,18]. The drawback of their use is also the numerical effort, resulting in longer computational times than the simplest 1D models. Numerically more efficient versions have been proposed in the past [19,20,21].

Still, Q2D models can be considerably improved with the development and investigation of new optimized meshes. The mesh definition significantly affects the computational effort, but also, more importantly, affects the accuracy of the results. Indeed, the accuracy of a numerical simulation is dependent on the mesh size and geometry [22]. A second important aspect for grid definition refers to the adaption of the mesh to the flow physics to ensure that the velocity variation is of the same order of magnitude along the grid. In practice, in pressurized transient analysis, a non-uniform grid must be generated due to the strong gradients near the wall; finer meshes are needed near the wall to better describe high velocity gradients, whereas coarser meshes are sufficient in the pipe core [23].

This paper aims at the assessment of radial mesh influence on the computation of unsteady energy dissipation in smooth–turbulent flows using a Q2D model. An extensive numerical analysis of the effect of the radial mesh on the computation of transient flows in pressurized pipes is carried out. A compromise between accuracy and computational time is established. A detailed analysis and comparison of the obtained Q2D results for the two best meshes for a theoretical instantaneous valve closure are developed.

The current paper innovatively presents the following: (i) the identification of the necessary improvements in the Q2D model radial meshes to accurately simulate steady and unsteady flows with less computation effort; (ii) guidelines to the definition of radial meshes based on transient simulation for 80 different meshes; (iii) in-depth analysis of the results for two different approaches to radial mesh construction; and (iv) the comprehensive characterization and analysis of transient flow simulation for the best identified mesh.

The numerical simulations consider a pressurized pipe system with a reservoir–pipe–valve configuration. The pipe has an inner diameter, D, of 20 mm, a wall thickness, e, of 0.001 m, a total length, L, of 15.2 m, and a wave velocity, c, of 1250 m/s (experimental and theoretical values of wave velocity yielded similar results). The numerical study considered an initial discharge of 120 × 10⁻³ L/s, a kinematic viscosity of 1.01 × 10⁻⁶ m²/s, corresponding to a Reynolds number of 7638, and the downstream valve’s theoretical instantaneous closure. It corresponds to the experimental facility assembled at the Laboratory of Hydraulics, Water Resources and Environment at Instituto Superior Técnico, which has been used for several experimental and numerical studies with different 1D, Q2D, and 3D models [13,24,25].

2. Numerical Models

Traditional axisymmetric equations, assumptions, and numerical schemes are used. The governing equations constitute a system of two partial differential equations with three unknown parameters (p, $u$, $\upsilon$) with the following form [14,17]:

$${\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{c^2}{g}}{\displaystyle \frac{\partial u}{\partial x}}=-{\displaystyle \frac{c^2}{gr}}{\displaystyle \frac{\partial \left(r\upsilon \right)}{\partial r}}$$

(1)

$${\displaystyle \frac{1}{g}}{\displaystyle \frac{\partial u}{\partial t}}+{\displaystyle \frac{\partial H}{\partial x}}={\displaystyle \frac{1}{r\rho }}{\displaystyle \frac{\partial r\tau }{\partial r}}$$

(2)

$$\tau =\rho \nu {\displaystyle \frac{\partial u}{\partial r}}-\rho \overline{u^{\prime }}\overline{{\upsilon }^{\prime }}$$

(3)

where x is the distance along the pipe; r is the distance from the axis in the radial direction; t is time; H is the piezometric head; $u$(x, r, t) is the local axial velocity; $\upsilon$(x, r, t) is the local radial velocity; g is the gravitational acceleration; c is the wave velocity; $\rho$ is the liquid density; $\tau$ is the shear stress; $\nu$ is the fluid kinematic viscosity; and $\rho \overline{u^{\prime }}\overline{{\upsilon }^{\prime }}$ represents the Reynolds stress term. The Reynolds stress calculation used the five-layer algebraic and time-constant turbulence model [14,26].

The main code was developed in VBA and some parts in C++, the parts corresponding to a higher number of computations. The implemented numerical scheme was the one proposed by Zhao and Ghidaoui [19], which allows to obtain the solution in N_C operations instead of the N_C³ necessary for solving the original Vardy and Hwang [14] model, with N_C being the total number of cylinders considered in the radial discretization. The numerical scheme, the mesh setup definition, the steady- and unsteady-state flow calculations, the boundary conditions, the truncation error, and the implemented five-region turbulence model are described in detail in Ferreira and Covas [22].

The radial mesh (or grid) generation used a geometric sequence (GS) defined by fixing the first cylinder thickness near the pipe wall, $\mathsf{\Delta }{r}_{j=N_C-1},$ and the geometric sequence common ratio, C_R. The thickness of cylinder j is calculated as follows:

$$\mathsf{\Delta }{r}_j=\mathsf{\Delta }{r}_{j=N_C-1}.{(1+C_R)}^{N_C-j}$$

(4)

A particular GS mesh is referred to as “GS _CR N_C”. For example, a GS_5%40 mesh, depicted in Figure 1, corresponds to a radial grid comprising a total of 40 cylinders (N_C = 40), with a 5% increase in cylinder thickness from the wall to the centre of the pipe (C_R = 5%).

3. Mesh Analysis

The most efficient mesh in the context of the current research stands for the one that ensures the best compromise between accuracy and computational time. For this purpose, a mesh sensitivity analysis is conducted, including two components as follows: (i) the steady-state analysis, which should assure a minimum initial mean velocity error; and (ii) the unsteady-state analysis, focusing on the unsteady wall shear stress computation.

3.1. Steady-State Mesh Assessment

The implemented Q2D model uses the method of characteristics to solve the transient flow propagation along the pipe axis and a finite difference numerical scheme for the shear stress and volume transfer calculation between adjacent cylinders in each pipe cross section. The latter numerical scheme has a truncation error proportional to ${\mathsf{\Delta }r}^2.{\partial }^{3}ru /{\partial }^{3}r$, which is proportional to ${\mathsf{\Delta }r}^2$, $O\left({\mathsf{\Delta }x}^2\right)$, and to the third derivative of the axial velocity.

The performance of the used finite difference scheme is analyzed based on the capacity to correctly represent the shear stress variation (i.e., high-gradient areas corresponding to higher mesh resolutions). In the case of a uniform shear stress variation, as in laminar flow conditions, an equal area mesh has been demonstrated to be a good approximation [22]. However, in turbulent flows, the difference between the total viscosity (i.e., the summation of the molecular or kinematic viscosity and the turbulent viscosity) near the wall and in the pipe core implies a steeper velocity gradient near the pipe wall. The high velocity gradient in this region requires a higher concentration of grid points near the wall, which can be achieved by using the geometric sequence (GS) mesh.

Since the refinement in the wall region is crucial, the adopted GS mesh is assessed for different numbers of cylinders, N_C, and of cylinder increase ratio values, C_R. The variation in the relative error of the mean velocity obtained by the GS mesh as a function of the N_C values is shown in Figure 2, corresponding each line to different C_R values. Note that the reference velocity, U_ref, corresponds to the mean velocity to which the velocity converges (when increasing the N_C) for each C_R value.

Two features should be highlighted. The first is that the mesh error decreases exponentially with N_C and the decay variation increases with C_R. For low N_C meshes (e.g., N_C = 10), the accuracy is almost C_R-independent. For intermediate N_C meshes, the accuracy increases with the C_R value. For N_C = 60, the C_R = 9% mesh attains a high accuracy with an error of 0.3%, while the C_R = 1% mesh leads to a much higher error of 6.65%. To achieve the same accuracy with the lowest C_R value (1%), a mesh with N_C = 120 should be considered (i.e., twice the radial grid points and computation time). The second feature is that, after reaching a certain accuracy value (N_C = 100), higher C_R meshes (C_R = 9% and C_R = 7%) present a significant decrease in the accuracy increase with N_C. However, low C_R value meshes are far from this pattern for the analyzed N_C range (Figure 2b).

Figure 3 depicts the dimensionless velocity, $u^{+}$(=$U/u_{\tau }$), with the dimensionless distance from the wall, $y^{+}$ (=$u_{\tau }y/\nu$), for two meshes, both with 60 cylinders: the first mesh with a C_R = 1% (Figure 3a) and the second with C_R = 9% (Figure 3b). The second mesh (C_R = 9%) allows to reduce the thickness of the first cylinder by almost two orders of magnitude (compared to C_R = 1%) by increasing the number of mesh points near the wall. This can be observed in Figure 3a,b by comparing the first mesh point near the pipe wall in both meshes. Results for C_R = 9% have a better fitting with the theoretical law-of-the-wall $u^{+}=y^{+}$ than those for C_R = 1%.

The previous examples justify that the error increases with the C_R decrease, due to the fact that the grid interval undergoes a notable expansion near the wall, where higher velocity gradients are observed. The effect of the steady-state velocity error is more significant in the calculation of the first pressure wave, in which the energy dissipation is small and less significant in the subsequent peaks. For the defined flow conditions (Q₀ = 120.6 × 10⁻³ Ls⁻¹), when using the Joukwosky formula ($∆H=c∆U/g$, where ∆H is the change in the pressure head, c is the fluid wave speed, ∆U is the change in the mean velocity, and g is the gravity acceleration), the initial mean velocity errors of 1% and 10% correspond to a maximum head of 0.49 and 4.91 m, respectively. These values correspond to 2 and 20 times, respectively, that of the initial steady-state friction head losses. The use of more sophisticated numerical models, such as those based on Q2D or 3D CFD approaches, is not a worthwhile option, unless a high level of accuracy in the determination of the steady-state velocity profile can be achieved.

3.2. Unsteady-State Mesh Assessment

3.2.1. The Rational

The optimization of the N_C and C_R parameters of the GS mesh is investigated based on transient simulations for 80 different meshes. The total number of cylinders, N_C, varies between 60, to control the steady-state velocity profile error, and 120, to avoid very time-consuming simulation times. The N_C increment is constant and equal to 20 (i.e., grids with 60, 80, 100, and 120 cylinders are used). The geometric rate, C_R, varies between 1% and 9%, with a step of 2%, as previously considered for steady-state evaluations (i.e., 1%, 3%, 5%, 7%, and 9%). For every combination of N_C and C_R, four different axial discretizations of the pipe, ∆x, are analyzed (i.e., 0.01, 0.02, 0.05 to 0.10 m).

Three parameters have been used to describe previous mesh characteristics (

Table 1. Meshes’ parameters considered in the 80 transient flow simulations.

${\mathit{N}}_{\mathit{C}}$	${\mathsf{\Lambda }}_{\mathit{R}}$	${\mathsf{\Lambda }}_{\mathit{X}}$
60	0.50, 1.00, 2.50, 5.00	0.18, 0.30, 0.43, 0.55, 0.63
80		0.20, 0.36, 0.53, 0.63, 0.71
100		0.22, 0.43, 0.60, 0.70, 0.75
120		0.23, 0.49, 0.67, 0.77, 0.80

): N_C, the total number of concentric cylinders; ${\mathsf{\Lambda }}_X$, the dimensionless mesh size in the pipe axial direction (${\mathsf{\Lambda }}_X=∆x/D$); and ${\mathsf{\Lambda }}_R$, the ratio between the number of cylinders in the three outers region (laminar and buffer I and II), N_3R, and the total number of cylinders, N_C (${\mathsf{\Lambda }}_R=N_{3R}/N_C$). The Courant–Friedrich–Lewy condition ($∆\mathrm{x}/\mathrm{c}∆\mathrm{t}=1$) is ensured when implementing the MOC along the pipe direction. Thus, by considering different ∆x values (${\mathsf{\Lambda }}_X$), the influence of ∆t is implicitly analyzed.

Since the energy dissipated decays exponentially with time, only the results for the first wave period are presented, but the results for five wave periods have also been calculated and analyzed, leading to similar conclusions. The reference mesh is GS_5%120 with ∆x = 0.01 m, which corresponds to the highest N_C, ${\mathsf{\Lambda }}_X$ (5.0) and ${\mathsf{\Lambda }}_R$ (0.67) values without numerical instabilities. The results obtained for this mesh have been compared with those from Zielke’s convolution model (i.e., lower ${\mathsf{\Lambda }}_R$ or N_C values do not ensure the same accuracy), and an exceptional accuracy in transient pressures, shear stresses, and mean velocities has been observed.

The influence of the mesh parameters on numerical accuracy is examined according to two main aspects: (i) the mesh capability to correctly calculate the velocity profile after the pressure surge arrival (i.e., properly describing the high-gradient velocity profile near the wall area), referred to as the short-term accuracy; and (ii) the grid capability to correctly describe the phenomenon when the created flow inversion (“vortex”) moves into the pipe core, referred to as the long-term accuracy.

3.2.2. Short-Term Accuracy

Figure 4 shows the comparison of the results obtained for the reference mesh (GS_5%120) with those from two grids with a lower C_R value (C_R = 1%): the first grid (GS_1%60) with half of the grid points, corresponding to half of the computational effort; and the second mesh (GS_1%120) with the same number of grid points as the benchmark. The results of the GS_5%120 mesh show the highest initial peak and a considerable decay with time, with an almost 80% amplitude reduction during the three-wave periods (Figure 4a).

The GS_1%120 mesh results only depict the reference grid behaviour after a short interval after each wall shear stress peak. The peaks tend to approximate to those of the reference mesh almost at the end of the simulation (t/T = 3.125), but, in the first peak, this mesh leads to a variation in the amplitude of almost 15% (Figure 4a).The GS_1%60 grid shows a lower accuracy in the description of the unsteady wall shear stress variation, with low and constant ${\tau }_{wu}$ peaks and a linear decay after each velocity profile inversion that does not describe the real exponential energy dissipation. Also, this grid does not correctly simulate the velocity profile changes, as shown by the inexistence of shear stress peaks after the wave reflection in the upstream reservoir (see t/T = 2.125 and t/T = 2.625 in Figure 4b).

Figure 5 illustrates the C_R influence on the shear stress peaks for a constant N_C value. The benchmark mesh results are compared with those from two meshes with half the N_C value, but with higher or equal C_R values (GS_5%60 and GS_9%60).

The effect of the mesh refinement in the wall area strongly depends on the C_R value. Higher C_R values correspond to more refined meshes in the wall area and courser meshes in the pipe core (for the same N_C). A low resolution near the wall meshes leads to lower and constant ${\tau }_{wu}$ peaks and a linear decay after each velocity profile inversion. This decay does not describe the real exponential energy dissipation, with it being responsible for the obtained rectangular shape of the pressure wave, unlike the rounded shape observed in the benchmark mesh.

The GS_9%60 mesh shows approximately the same performance as the benchmark mesh (Figure 5a). It allows a very precise energy dissipation in each pressure wave arrival and practically eliminates the discrepancy over time (Figure 5b). This excellent energy dissipation calculation is responsible for a pressure history identical to that of the benchmark mesh. The GS_5%60 simulation results are less accurate than those for GS_9%60 without the exponential ${\tau }_{wu}$ peak decay over time, and, for that reason, this grid does not describe the rounder ${\tau }_{wu}$ peak observed in the third wave period (Figure 5c). Nevertheless, GS_5%60 mesh produces an almost perfect pressure time variation.

The previous examples highlight the importance of C_R increasing to reproduce the real energy dissipation in unsteady flows (compare GS_1%60 and GS_9%60 simulations) and the higher effect of C_R than of N_C on the accuracy of the results (compare GS_1%120 and GS_9%60 results).

Figure 6 depicts the effect of pipe axial discretization for three values of ${\mathsf{\Lambda }}_R$ (0.22, 0.60, and 0.75) with a fixed C_R value (100), considering three spatial discretizations (${\mathsf{\Lambda }}_X$ = 5.0, 1.0, and 0.5). Results of the mesh with the lowest ${\mathsf{\Lambda }}_R$ values (${\mathsf{\Lambda }}_R$ = 0.22) are independent of ${\mathsf{\Lambda }}_X$ (no benefit from ∆x reduction). The mesh with an intermediate ${\mathsf{\Lambda }}_R$ value (${\mathsf{\Lambda }}_R$ = 0.60) benefits considerably from the first reduction of ${\mathsf{\Lambda }}_X$ (from 5.0 to 1.0) and presents similar results to the ${\mathsf{\Lambda }}_R$= 0.75 mesh for ${\mathsf{\Lambda }}_X$= 1.0. The highest resolution mesh at the wall (${\mathsf{\Lambda }}_R$ = 0.75) needs the maximum pipe discretization (the highest values for ${\mathsf{\Lambda }}_X$= 0.5) to ensure the accurate calculation of the highest peaks. Therefore, low-resolution meshes (lower ${\mathsf{\Lambda }}_R$ values) do not benefit from further pipe axial discretization. High-resolution meshes (higher ${\mathsf{\Lambda }}_R$ values), without proper pipe axial discretization, present numerical instabilities and no improvement in terms of accuracy.

Table 2. Wall shear stress peak value at the pipe midsection for all simulations (80 meshes).

${\mathbf{N}}_{\mathit{C}}$	${\mathsf{\Lambda }}_{\mathit{R}}$	${\mathbf{C}}_{\mathit{R}}$	${\mathsf{\Lambda }}_{\mathit{X}}$
			5.0	2.5	1.0	0.5
60	0.18	1	5.22 (−95.1%)	5.25 (−95.1%)	5.27 (−95.1%)	5.27 (−95.1%)
	0.30	3	11.00 (−89.7%)	11.20 (−89.5%)	11.40 (−89.3%)	11.40 (−89.3%)
	0.43	5	21.40 (−79.9%)	23.10 (−78.3%)	24.20 (−77.3%)	24.80 (−76.7%)
	0.55	7	32.80 (−69.3%)	39.90 (−62.5%)	45.40 (−57.4%)	51.20 (−52.0%)
	0.63	9		50.50 (−52.6%)	63.00 (−40.9%)	84.40 (−20.5%)
80	0.20	1	7.90 (−92.6%)	8.00 (−92.5%)	8.10 (−92.4%)	8.10 (−92.4%)
	0.36	3	19.60 (−81.6%)	20.80 (−80.4%)	21.70 (−79.6%)	22.10 (−79.3%)
	0.53	5		42.50 (−60.1%)	51.90 (−51.3%)	56.70 (−46.8%)
	0.63	7				94.80 (−11.0%)
	0.71	9
100	0.22	1	11.40 (−89.3%)	11.60 (−89.1%)	11.80 (−88.9%)	11.80 (−88.9%)
	0.43	3	29.30 (−72.5%)	33.90 (−68.2%)	38.10 (−64.3%)	39.80 (−62.6%)
	0.60	5				92.90 (−12.8%)
	0.70	7
	0.75	9
120	0.23	1	15.10 (−85.8%)	15.70 (−85.3%)	16.10 (−84.9%)	16.20 (−84.8%)
	0.49	3		45.40 (−57.4%)	57.01 (−46.4%)	63.70 (−40.2%)
	0.67	5				106.6
	0.77	7
	0.80	9
Numerical instabilities did not allow to obtain a solution
Short-term accuracy was not influenced by a higher pipe discretization
A higher conduit discretization had a limited influence on short term-accuracy
The benchmark (or reference) mesh

presents the calculated ${\tau }_w$ peak in the first (t/T = 0) wave period and compares these values with those from the benchmark mesh (GS_5%120). The values presented between brackets represent the relative error between each simulation and the reference mesh simulation (“green” cell). The results of the simulations with numerical instabilities are not presented in this table (see “grey” cells).

For meshes with lower ${\mathsf{\Lambda }}_R$ and N_C values (i.e., low refinement in the wall area), there is no reason to decrease the ${\mathsf{\Lambda }}_X$ parameter (i.e., increase the pipe discretization); this is shown in Table 2 marked by “orange” shade for N_C = 60 with ${\mathsf{\Lambda }}_R$ = 0.18 and 0.30; for N_C = 80 meshes with ${\mathsf{\Lambda }}_R$ = 0.20; and for N_C = 100 meshes with ${\mathsf{\Lambda }}_R$ = 0.22. A second set of meshes is flagged with a higher ${\mathsf{\Lambda }}_R$ value, for which the ${\mathsf{\Lambda }}_X$ reduction corresponds to a moderate accuracy improvement. These meshes are marked with a “red” shade. Finally, for higher ${\mathsf{\Lambda }}_R$ meshes, the ${\tau }_w$ simulation considerably benefits from a ${\mathsf{\Lambda }}_X$ decrease. All simulations benefit from a ${\mathsf{\Lambda }}_R$ increase for the same ${\mathsf{\Lambda }}_X$, independently of the N_C value.

For the same accuracy, the computational effort can be optimized by combining ${\mathsf{\Lambda }}_{\mathrm{X}}$ and ${\mathsf{\Lambda }}_{\mathrm{R}}$. The GS_7%60 mesh (${\mathsf{\Lambda }}_R$ = 0.55 and N_C = 60) with ${\mathsf{\Lambda }}_X$ = 5.0 ensures approximately six times a higher shear stress first peak at t/T = 0 than a lower resolution mesh (${\mathsf{\Lambda }}_R$ = 0.18) with the same N_C. For ${\mathsf{\Lambda }}_R$ = 0.55, the ${\mathsf{\Lambda }}_X$ decrease significantly improves the accuracy, leading to a ten-fold higher peak for ${\mathsf{\Lambda }}_X$ = 0.5 than ${\mathsf{\Lambda }}_{\mathrm{X}}$ = 5.0. To improve the accuracy of GS_7%60 mesh, a careful optimization of the parameters must be carried out. For instance, for N_C = 80, the ${\mathsf{\Lambda }}_X$ value must be less than 1.0 and ${\mathsf{\Lambda }}_R$ should be the highest value. For N_C = 100 only, with the highest ${\mathsf{\Lambda }}_R$ value (0.75), only the benchmark with N_C = 120 has shown better results.

3.2.3. Long-Term Accuracy

Long-term accuracy is assessed by comparing the value of the wall shear stress immediately before the arrival of the pressure wave reflected from the upstream reservoir for the 80 meshes. The objective is to determine whether the increase in ${\mathsf{\Lambda }}_R$ threatens the long-term accuracy of the simulation due to the inability of the mesh to correctly calculate the shear stress as the flow vortex created by the valve closure moves towards the pipe core due to insufficient grid points.

The accuracy is generally high (compared to the benchmark), with relative errors lower than 1.5% in most cases. The worst long-term accuracy is observed for the lowest N_C and ${\mathsf{\Lambda }}_R$ meshes (N_C = 60 with ${\mathsf{\Lambda }}_R$ = 0.18), since the initial wall shear peak values are underestimated, though with relative errors being relatively low (6.2–6.9%). For N_C = 60 meshes using higher ${\mathsf{\Lambda }}_R$ values, the long-term accuracy is high. Overall, long-term accuracy tends to benefit from an increase in ${\mathsf{\Lambda }}_R$ (especially for low N_C meshes) and a decrease in ${\mathsf{\Lambda }}_X$, the latter with less impact.

3.3. Numerical Stability

For a specific N_C value, the preceding simulations have demonstrated numerical instabilities as a consequence of an increase in ${\mathsf{\Lambda }}_R$. A reduction in ${\mathsf{\Lambda }}_X$ mitigates this effect. As the N_C value increases, the numerical instabilities occur for lower ${\mathsf{\Lambda }}_R$ values (see Table 2 where simulations with numerical instabilities are highlighted in “grey”).

An analysis for the complete simulation set is carried out to identify the correlation between N_C, ${\mathsf{\Lambda }}_R$, and ${\mathsf{\Lambda }}_X$ and the occurrence of numerical instabilities. No correlation is observed between the numerical instability and a particular parameter, rather being a combination of parameters. In fact, numerical instabilities tend to appear for a higher number of cylinders in the outer layers, N_3R (i.e., $N_C\times {\mathsf{\Lambda }}_R$), because these correspond to a higher mesh resolution near the wall and, thus, lead to higher axial velocity gradients immediately after the pressure surge. On the other hand, to maintain the numerical stability for higher ${\tau }_w$ peak values, the time interval must decrease (lower ${\mathsf{\Lambda }}_X$ values) since the method of characteristics used to calculate the transient in the axial direction is required to comply with the Courant–Friedrich–Lewy condition.

Figure 7 plots the complete set of simulations in terms of dimensionless parameters, identifying the results with or without numerical instabilities. A relation between N_3R/N_X (N_X being the number of sections in the axial direction) and ${\mathsf{\Lambda }}_R$ and simulations with numerical instability is observed.

For meshes with a ${\mathsf{\Lambda }}_R$ ratio smaller than 0.5, no numerical instability is observed in the 80 simulations; for ${\mathsf{\Lambda }}_X$ = 0.5 meshes (“yellow” marks), the ${\mathsf{\Lambda }}_R$ can increase until 0.7 without any numerical instability; for ${\mathsf{\Lambda }}_X$ = 1.0 meshes (“green” marks), the ${\mathsf{\Lambda }}_R$ is limited to 0.6; for ${\mathsf{\Lambda }}_X$ = 2.5 meshes (“red” marks), the ${\mathsf{\Lambda }}_R$ is limited to 0.6; and for ${\mathsf{\Lambda }}_X$ = 5.0 (“blue” marks), the ${\mathsf{\Lambda }}_R$ is limited to 0.5. Figure 7 demonstrates that the numerical stability is a direct consequence of the mesh resolution in the wall area and is not solely dependent on ${\mathsf{\Lambda }}_R$. Additionally, the stability dependence on ${\mathsf{\Lambda }}_X$ is also recognized. Vardy and Hwang [8] found an optimal ratio N_X/N_C = 1.5 for the accuracy and stability of the solutions, showing that a N_X/N_C higher than 1.5 reduces the accuracy and a N_X/N_C lower than 1.5 leads to numerical instabilities. This aspect is true for the equal area cylinder meshes used by Vardy and Hwang because there is a constant relation between N_3R and N_C. Zhao and Ghidaoui [19] argued that, for a geometric sequence mesh (GS) with a fixed N_C value, a correlation between ${\mathsf{\Lambda }}_R$ and ${\mathsf{\Lambda }}_X$ must be ensured. The presented stability analysis adds the importance of the N_3R parameter and, consequently, the correlation between C_R and N_C.

The observed limits for ${\mathsf{\Lambda }}_R$ herein refer to this particular pipe system and the instantaneous valve closure. Should an S-shaped valve closure be simulated, the velocity gradient near the wall would be reduced, and the majority of the previous numerical instabilities would not be observed. In that case, higher ${\mathsf{\Lambda }}_X$ values could be used for the same ${\mathsf{\Lambda }}_R$.

3.4. The Most Efficient Mesh

The accuracy of transient simulation enhances significantly with a ${\mathsf{\Lambda }}_R$ increase. Yet, increasing the ${\mathsf{\Lambda }}_R$ is limited by the simulation numerical stability and, thus, must be followed by a ${\mathsf{\Lambda }}_X$ decrease. To ensure an accurate simulation with a considerable reduction in the computer effort, the simultaneous optimization of both parameters, ${\mathsf{\Lambda }}_R$ and ${\mathsf{\Lambda }}_X$, is carried out for the 48 meshes without numerical instabilities (Table 2). The computational effort is measured by the ratio between the total number of grid points (N_C × N_X) and the lowest number of grid points of the 48 meshes (i.e., 9000, corresponding to N_C = 60 and Δx = 0.10 m). Since the Zhao and Ghidaoui [19] model is used, the computational effort increases linearly with the total number of grid points. For example, a mesh with N_C = 120 and Δx = 0.01 m (total number of grid points equal to 180,000) has a computational effort 20 times that of the mesh with the lowest number of grid points.

The accuracy is assessed comparing the piezometric head results for each mesh with those of the highest resolution mesh (GS_5%120 with ${\mathsf{\Lambda }}_R$ = 0.49 and ${\mathsf{\Lambda }}_X$ = 0.5) by means of the mean absolute percentage error (MAPE) in three-wave periods using the following equation:

$$\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{E}={\displaystyle \frac{100}{n}}\sum _{t=0}^n\left|{\displaystyle \frac{X_i-{\widehat{X}}_i}{X_i}}\right|$$

(5)

where $X_i$ is the true value (given by benchmark mesh value) for time t, ${\widehat{X}}_i$ is the estimated value (simulated value) for time i, and n is the number of time steps in three-wave periods. The difference between $X_i$ and ${\widehat{X}}_i$ is divided by the benchmark value $X_i$ again.

3.4.1. Transient Results Calculation for the Piezometric Head in the Dimensionless Form

The dimensionless form of the piezometric head, $(H–H_0)/\left(a{U}_0/g\right)$, is used herein. This approach only accounts for the unsteady-state error (i.e., it cancels the steady-state simulation error associated with $H_0$ and $U_0$). Figure 8a presents computational effort versus the MAPE for the 48 simulations, whereas Figure 8b only shows the solutions obtained with MAPE < 0.20%. The meshes are divided into four groups according to the ${\mathsf{\Lambda }}_{\mathrm{R}}$ parameter: (i) the “blue” group with ${\mathsf{\Lambda }}_R$ ≤ 0.2; (ii) the “orange” group with ${\mathsf{\Lambda }}_R$ between 0.2 and 0.3; (iii) the “grey” group with ${\mathsf{\Lambda }}_R$ between 0.3 and 0.5; and (iv) the “yellow” group with the highest ${\mathsf{\Lambda }}_R$ (≥0.50) meshes.

The lowest ${\mathsf{\Lambda }}_R$ meshes (the blue and orange groups) present the most unsatisfactory results, since they have a similar computational effort as the other groups but with lower accuracy (i.e., higher MAPE). The grey group contains meshes with the best results (Meshes 1, 2, and 3, see Figure 8b), though having a lower ${\mathsf{\Lambda }}_R$ value than the yellow group. This result is explained by the fact that the meshes with higher N_C values (yellow meshes) were removed from the analysis due to numerical instabilities. For instance, the yellow group does not have a valid mesh for N_C = 120 and has one single mesh with C_R = 5% (Mesh 5) for N_C = 100.

Mesh 1 (see Figure 8b) ensures the best relationship between computational effort and accuracy. If a higher accuracy is needed, the best option is to decrease the axial discretization interval to Δx = 0.02 m, as obtained by Mesh 2 (the same radial mesh, but a five-fold higher axial discretization). On the contrary, if a lower resolution is accepted, the best option is to decrease N_C and to increase C_R, as represented by Mesh 3. This analysis highlights that the Δx increase is a good solution, when no numerical instability is observed, and its decrease is only justified to control the occurrence of numerical instabilities. The dashed grey line, containing Meshes 1, 2, and 3, defines the Pareto front (i.e., the set of non-dominated solutions) of this two-objective optimization problem, and all the other meshes correspond to dominated solutions.

3.4.2. Transient Results Calculation for the Piezometric Head in the Dimensional Form

A high-accuracy transient simulation is only justified if a high-accuracy steady-state velocity profile is ensured by the Q2D model, as previously stated (see Section 3.1). Figure 9 shows the computational effort versus MAPE for the 48 simulations using the piezometric head results to account for the steady-state velocity profile error.

The MAPE shows a significant increase (compared with Figure 8), highlighting the importance of correctly calculating the steady-state velocity profile and the relatively low transient error provided by the Q2D models, as previously referred. The yellow group presents the lowest transient error since it has the best steady-state simulation, and, consequently, it presents better results than the grey meshes. The best result is obtained for Mesh 7 (corresponding to GS_9%60 with Δx = 0.05 m mesh). The previous mesh has the highest C_R value, and, consequently, the lowest steady-state error. The dashed yellow line corresponds to the Pareto front in this case.

In conclusion, the Q2D model with the proper selection of ${\mathsf{\Lambda }}_R$ and ${\mathsf{\Lambda }}_X$ correctly simulates the transient flow energy dissipation with the lowest N_C value (i.e., N_C = 60). The steady-state error is commonly ignored, though it considerably influences the transient state results, which can be improved by using higher C_R values.

Mesh 7 (GS_9%60) is the best option when considering the piezometric head (i.e., considering steady-state mean velocity error), and Mesh 1 (GS_3%120) is the best option when the dimensionless piezometric head is considered (i.e., when the steady-state error is neglect).

4. Performance of the Two Best Meshes

4.1. Aim and Meshes Characterization

The results obtained for the two best meshes, GS_9%60 and GS_3%120, are compared for a five-period pressure wave simulation with an instantaneous valve closure. The comparison aims to identify differences in the results associated with the C_R optimization.

The two best meshes, GS_9%60 and GS_3%120, have a ${\mathsf{\Lambda }}_X$ value equal to 2.5 (i.e., ∆x = 0.05 m). The GS_9%60 mesh has a higher ${\mathsf{\Lambda }}_X$value (0.63 compared with the 0.49 value of the GS_3%120 mesh) to compensate the smaller N_C value. The GS_9%60 mesh reduces to half the first cylinder thickness. Comparing the two meshes number of grid points, the GS_9%60 mesh has approximately the same number of grid points in the viscous layer, half the number of points in the buffer I and II layers, and one third of the points in the central regions (logarithmic and core regions).

4.2. The Steady-State Flow Simulation

Figure 10 presents the dimensionless velocity profile for the steady-state flow for the two meshes ($u^{+}$ as a function of $y^{+}$), with both velocity profiles showing the same increase with y+. The turbulent steady-state dimensionless velocity profile is particularly sensitive to the radial mesh interval between the viscous and log-layer due to the sharp increase in velocity observed. The GS_9%60 mesh shows a correct transition between the viscous and log-layer with a courser grid. Both meshes show an excellent performance, demonstrated by (i) the velocity profile matching the law of the wall ($u^{+}={\displaystyle \frac{1}{\kappa }}ln{y}^{+}+C)$ for values of $y^{+}$ > 30; (ii) the velocity profile variation close to the wall (viscous sublayer) being linear with $y^{+}$ ($u^{+}=y^{+}$); and (iii) the gradual velocity profile variation observed between the viscous and log-layer. Yet, GS_9%60 has half the grid points of GS_3%120 for the same accuracy (performance).

4.3. The Unsteady Flow Simulation

The time-history of piezometric head at the pipe mid and end sections in the first five periods of the transient are simulated by the Q2D model and are presented in Figure 11 for both meshes (GS_9%60 and GS_3%120). Results are presented in the dimensionless form. Hardly any differences are observed: the two grids lead to an identical pressure wave damping, associated with the correct simulation of energy dissipation after the velocity shift, and the development of a similar S-shape wave form associated with the energy dissipation immediately after the pressure surge.

The wall shear stress results for both meshes are presented in Figure 12. Once again, differences between the results of the two meshes are minor and restricted to the first pressure peaks and during a short time interval after the arrival of each wave front, after which the results converge very quickly. The GS_9%60 mesh results have slightly higher wall shear stress peaks.

The correlation between the wall shear stress and the mean velocity at the pipe midsection for both meshes (GS_3%120 and GS_9%60) is shown in Figure 13. Results are presented for the first, the second, and the third wave periods. In the first pressure surge (t/T=0.125, see Figure 13a), the mean velocity reduction is abrupt (with a rectangular pressure wave shape, as observed in Figure 11). The two meshes show different peak values, with the GS_9%60 mesh generating a higher ${\tau }_W$value, which reduces the maximum mean velocity variation compared with GS_3%120. Higher ${\tau }_W$ values lead to a smoother pressure wave. The GS_9%60 mesh increases the pressure wave smoothing, and doing so tends to attenuate the difference between the two meshes results. Consequently, the results in the third period are almost identical for both meshes (Figure 13c). As the simulation proceeds, there is a rounding of the pressure wave and a progressive decay of the maximum mean velocity variation between time intervals, as observed by the elliptical shape in Figure 13c. The two meshes present similar performances, as previously verified for the wall shear stress analysis, because the mesh accuracy requirement in the wall region decreases.

The velocity profiles are also evaluated, and minimal differences are observed between the two meshes, with the exception of regions in close proximity to the wall and immediately following the arrival of the initial pressure wave. These observations are coherent with the previous results.

The two meshes yielded comparable results. However, the use of a low-resolution mesh on the wall (low N_C combined with low C_R) will result in an inability to accurately calculate the peak value of wall shear stress, and the transformation from a rectangular shape into an S-type shape pressure wave will not occur. The time-history of the piezometric head will continue to exhibit abrupt variations without the associated dissipation of transient state energy.

5. Transient Dynamic Characterization

The dynamic characterization of the transient phenomenon is further developed using the GS_9%60 mesh. The dimensionless mean velocity time-history is presented in Figure 14. The profiles t₀ to t₈ describe the instants immediately before and after each pressure variation in the first period (0.125 ≤ t/T ≤ 1.125) at the pipe midsection. Instants t₉ to t₁₇ describe the homologous instants, but for the fifth wave period (4.125 ≤ t/T ≤ 5.125).

Figure 15 shows the local acceleration and the wall shear stress time-history, both in absolute values. The critical features include the following: (i) the wall shear stress peak values directly correspond to the local acceleration peaks; and (ii) the switch from an instantaneous to gradual variation, as previously observed concerning the piezometric head time-history. However, the wall shear stress presents a slower decay and tends to a non-zero value, while the local acceleration presents an initial peak, followed by a fast decay and almost immediately tending to zero. This difference justifies the instantaneous acceleration-based model’s main shortcoming.

For the analyzed instants (t₀ to t₁₆), the shear stress profile and the dimensionless velocity profile are presented in Figure 16 and Figure 17, respectively, both plotted with y+ in a logarithmic scale to better distinguish the five turbulent regions.

The shear stress profile allows to analyze the mesh influence on the simulation accuracy because the numerical scheme truncation error at each grid point corresponds to the local shear stress gradient multiplied by the grid dimension. Thus, a high shear stress gradient area should correspond to a higher mesh refinement, and, conversely, a linear profile allows for the use of courser meshes. The shear stress profile reacts to the pressure surge and the no-slip condition with an immediate change in the wall region (for an illustration of this, see t₁ and t₅ in Figure 16). In the equilibrium phase, the peak reduces as the “vortex” sheet moves into the pipe core. The wall shear stress peaks diminish with time due to the pressure surge smoother variation, observed by the rounder wall shear stress profile (the previous S-shape effect observed in Figure 14). Nevertheless, the surge influence is limited to the viscous and buffer I layers and is not considerably changed along the simulation. Inside the buffer I layer, the wall shear stress presents a lower and linear gradient that is accurately described by a coarser mesh.

The velocity profile is the result of the equilibrium process up to that moment. As the transient propagates, successive vortices are created that move toward the pipe core, and their effect is strongly dissipated along time. In Figure 17, the first peak (t₁) effect during the initial quarter-wave period spreads until the buffer I layer (t₂). At the end of the first half-wave period (t₄), the steady profile (t₀) is mainly affected from the wall to the buffer II layer, whereas, at the end the fifth period (t₁₆), the effect reaches the log layer. Nevertheless, during the complete analysis, the high shear stress gradients are limited to the viscous and buffer I layers. This explains why the two meshes present similar results (GS_9%60 has 31 grid points and GS_3%120 has 43 grids points in the two outer layers).

Figure 18 illustrates the axial velocity time-history at different distances, r, from the pipe centerline described as a percentage of the pipe radius, R, from r = 50%R to 92.5%R. Figure 18a presents four grid points that are situated within the core, log, and buffer II layer (the log and buffer II layers’ upper limits correspond to r = 85%R and r = 92.5%R). Figure 18b provides a detailed description of the axial velocity at several points in the viscous layer (99%R to 99.9%R). All points within the pipe core (Figure 18a) follow the mean velocity trace and are not affected by the high-velocity gradient in the pipe wall (with no spikes and quite smooth changes). Figure 18b illustrates a contrasting behaviour amongst the core points, which can be readily distinguished from the wall region points with a fast response to the pressure variations. This change is evident, with the viscous layer points exhibiting a performance analogous to that of the wall shear stress. The closer the points are to the pipe wall, the similar the velocity traces to the wall shear stress are, as observed in Figure 18b.

The preceding evaluations have highlighted the need to increase the mesh refinement in the outer layers (predominantly in the viscous and buffer I layers), given that this is the region where the flow dynamics undergo the most significant and rapid changes. This provides further evidence to support the hypothesis that an increase in C_R values may be beneficial. Although the results show some variation in the core over time, the fluctuations are markedly less pronounced when this occurs. The necessity for controlling the finite difference truncation error indicates that a reduction in the mesh may be a viable approach in instances where higher shear stress variations are observed.

6. Conclusions

The utilization of a geometric sequence grid for the simulation of transient turbulent flow is supported by the high concentration of mesh points distributed between the pipe wall and the log layer. The accuracy of the unsteady simulation is enhanced by higher C_R values, which allows the use of a high-resolution mesh in the wall region, while maintaining a low N_C value and a lower computation time. The value of the conduit axial discretization should be sufficiently high to prevent numerical instabilities.

The comparison between the two best performance meshes—GS_9%60 and GS_3%120—with a total number of 60 and 120 cylinders and C_R of 9% and 3%, respectively, has demonstrated the benefit of the C_R optimization. The two meshes show small differences for the instantaneous valve closure, which are limited to a reduced number of time steps after the pressure surge arrival and solely during the first pressure waves. The GS_9%60 mesh has the advantage of requiring half of the computational time. The transient dynamic characterization reveals the in-depth consistency between the model results and the hydraulic transient phenomenon. The high shear stress gradient is confined to the viscous and buffer I layer (the inner part of the flow is in phase with the mean velocity), thereby justifying the refinement of the radial mesh in this specific area.