Energy Analysis of a Quasi-Two-Dimensional Friction Model for Simulation of Transient Flows in Viscoelastic Pipes

Abstract

Quasi-two-dimensional (quasi-2D) friction models have been widely investigated in transient pipe flows. In the case of viscoelastic pipes, however, the effect of different values of the Reynolds number (Re) on pressure fluctuations (which can lead to water hammer) have not been considered in detail. This study establishes a quasi-2D friction model employing an integral total energy method and investigates the work due to frictional and viscoelastic terms at different Re values. The results show that viscoelastic work (W_P) and frictional work (D_f) increase with an increase in Re. However, when the initial Re values are high, the D_f values are much larger than the W_P values. In addition, for Re < 3 × 10⁵, the 1D model underestimated the viscoelastic terms. There was no significant difference between the two models for Re > 3 × 10⁵. In the case of different initial Re values, the two models produced almost the same values for W_P. This study provides a theoretical basis for investigating transient flow from the perspective of energy analysis.

1. Introduction

In recent years, viscoelastic pipes have been increasingly used in urban water-supply systems [1]. The transient-flow pressure fluctuation generated by such pipes is different from the pressure fluctuation of traditional elastic pipes because of their different viscoelastic properties; when pumps suddenly stop or valves close quickly, a water-hammer accident result. The peak values of pressure fluctuation are relatively small; however, pressure-fluctuation attenuation and phase delay are large. Therefore, a correct understanding of the role of viscoelastic characteristics in the transient-flow pressure-fluctuation process is necessary.

Accurate numerical calculation of transient flow in viscoelastic pipes must consider the effect of not only unsteady friction (UF) of the pipe wall but also wall viscoelasticity. Some researchers have studied the pressure fluctuation of transient flows in viscoelastic pipes. Rieutord and Blanchard [2] pointed out that the acceleration of pressure attenuation is due to a time offset between the pressure and the retarded strain of the pipe wall. Covas et al. [3] proposed that, owing to retarded deformation, the maximum pressure of the transient flow of viscoelastic pipes is smaller than that of elastic pipes, and the attenuation of the pressure wave is faster. In addition, Covas et al. [4] found that the viscoelastic effect is more evident in a transient flow when the retarded time of the viscoelastic pipes is less than the propagation time of the pressure wave along the pipes.

Numerous studies on the influence of UF on the transient flow of elastic pipes have been conducted. The results show that a one-dimensional (1D) quasi-steady friction model can accurately simulate the maximum value of pressure fluctuation in elastic pipelines, but in most cases, it cannot accurately describe the peak-pressure damping [1]. One-dimensional UF model [5] and two-dimensional (2D) friction models [6] can accurately simulate not only the maximum peak of pressure fluctuations but also the attenuation of pressure fluctuation. Pezzinga [6] compared the error values of 1D- and 2D-model-simulated maximum and minimum pressure heads with experimental results and confirmed the accuracy of the 2D-model-simulated shear stress. In a study of different friction models of the transient flow through a simulated viscoelastic pipeline, Firkowski et al. [7] examined the effects of UF and experimentally obtained creep functions. Urbanowicz et al. [8] investigated quasi-steady and unsteady frictions in viscoelastic pipes and presented a simplified effective numerical-convolution integral to describe retarded strain. Sun et al. [9] concluded that a 1D quasi-steady friction model could accurately simulate the pressure-fluctuation attenuation of viscoelastic pipes at a certain water temperature.

In previous studies, the viscoelastic parameter was mistakenly thought to be energy dissipation, but an increasing number of studies have shown that there is energy transfer between the fluid and the pipe wall [10]. Karney [11] divided the energy involved in transient flow into four parts: the kinetic energy of the fluid, the internal energy associated with fluid compressibility and pipeline-elasticity effects, the energy dissipated by friction, and the work done at the ends of the pipeline. Through integral processing, the basic momentum and continuity equations of transient flow can be used to analyse the energy of the entire pipeline. Duan et al. [10] used Fourier methods and energy analysis to show that energy transfer was related to the ratio of the pressure wave period to the retarded time (T/τ). Duan et al. [12] used local transient analysis (LTA) and integral total energy (ITE) methods to study the influence of friction on transient flow in a 1D UF model and 2D k-ε turbulence model. The results showed that the friction dissipation calculated by the 1D model was smaller than that calculated by the 2D model after the first pressure-fluctuation period. Meniconi et al. [13] introduced nondimensional parameters to analyse the energy dissipation of the transient flow in viscoelastic pipes, finding numerically that the energy dissipation was related to the Reynolds number (Re) and valve opening. Riasi et al. [14] employed a 2D k-ε turbulence model to infer nondimensional parameters related to Re and the viscoelasticity of the pipe. The influence of these parameters on the pressure fluctuation, wall-shear stress, velocity profiles, turbulence production, and dissipation was investigated.

As mentioned, Re is an important parameter in transient flow. Hence, an understanding of how it influences viscoelastic and frictional effects on the energy dissipation of transient flow is critical for the design and operation of a pipeline system. In this study, transient flow with different initial Re values was investigated for both 1D and 2D models. First, the energy equation of the 2D model was derived using the ITE method. The contributions of the viscoelastic and friction terms to the work done in both friction models were analysed and discussed, and the change rule of the work done under different Re values were clarified. The trends of work done and energy dissipation at different Re values were examined from the perspective of energy analysis, and the laws governing fluctuation of friction and viscoelastic terms at different stages of transient flow were studied. Our study established a quasi-2D model based on energy analysis to analyse the friction and viscoelastic effects of transient flow at different Reynolds number.

2. Materials and Methods

2.1. Governing Equations

The continuity and momentum equations of the quasi-2D transient flow model are expressed as [4,5]:

$$\frac{g}{a^2}\frac{\partial H}{\partial t}+\frac{\partial u}{\partial x}+\frac{1}{r}\frac{\partial \left(rv\right)}{\partial r}+2\frac{\partial {\epsilon }_r}{\partial t}=0,$$

(1)

$$\frac{\partial u}{\partial t}+g\frac{\partial H}{\partial x}-\frac{1}{\rho r}\frac{\partial \left(r\tau \right)}{\partial r}=0,.$$

(2)

where $H$ is the pressure head, $Q$ is the discharge,$a$ is the wave speed, r is the radial distance from the pipe centre, $x$ is the axial coordinate along the pipe, $t$ is the time, u is the longitudinal velocity, v is the radial velocity, $\tau$ is the shear stress, ${\epsilon }_r$ is the retarded strain, $g$ is the gravitational acceleration, and $\rho$ is the density.

By integrating Equations (1) and (2) over the cross-sectional area of the pipeline, the continuity and momentum equations for the 1D transient flow in viscoelastic pipes can be expressed as [4]:

$$\frac{\partial H}{\partial t}+\frac{a^2}{gA}\frac{\partial Q}{\partial x}+\frac{2a^2}{g}\frac{\partial {\epsilon }_r}{\partial t}=0,$$

(3)

$$\frac{\partial Q}{\partial t}+gA\frac{\partial Q}{\partial x}+\frac{\pi D}{\rho }{\tau }_w=0,$$

(4)

where $A$ is the cross-sectional area of the pipeline, and ${\tau }_w$ is the pipe-wall shear stress evaluated using the Darcy–Weisbach formula.

2.2. Kelvin-Voigt Model

The behaviour of the viscoelastic pipes was simulated using the Kelvin-Voight (K-V) model [15,16]. Each of the K-V components consists of a spring and sticky pot, representing the instantaneous strain and retarded strain components, respectively. The total creep compliance is given by:

$$J\left(t\right)=J_0+{\displaystyle \sum }_{k=1}^N{J}_k\left(1-e^{-\frac{t}{{\tau }_k}}\right),$$

(5)

where J₀ = 1/E₀, E₀ is the elastic modulus, $J_k=1/E_k$ is the creep compliance of the k-th element, ${\tau }_k$ is the retarded time of the k-th element, and $E_k$ is the elasticity modulus of the k-th element.

For viscoelastic pipes, the total strain ε is given by the sum of the instantaneous strain (ε_e) and retarded component (ε_r):

$$\epsilon ={\epsilon }_e+{\epsilon }_r.$$

(6)

Using the K-V model with N elements, the retarded strain is the sum of the single-element deformations [3]:

$$\partial {\epsilon }_r\left(t\right)={\displaystyle \sum }_{k=1\cdots N}\partial {\epsilon }_{rk}\left(t\right).$$

(7)

The retarded strain is [4]:

$${\epsilon }_r\left(x,t\right)={\displaystyle \sum }_{k=1\cdots N}{\epsilon }_{rk}\left(x,t\right)={\displaystyle \sum }_{k=1\cdots N}\frac{\gamma \alpha D}{2e}{{\displaystyle \int }}_0^t\left[H\left(x,t-t^{\prime }\right)-H_0\left(x\right)\right]\frac{J_k}{{\tau }_k}{e}^{\frac{-t^{\prime }}{{\tau }_k}}d{t}^{\prime },$$

(8)

where t’ is the time variation.

The retarded strain rate is algebraically expressed as [17,18]:

$$\frac{\partial {\epsilon }_r\left(t\right)}{\partial t}={\displaystyle \sum }_{k=1\cdots N}\frac{\partial {\epsilon }_{rk}\left(t\right)}{\partial t}={\displaystyle \sum }_{k=1}^N\left(AH-V_E\right),$$

(9)

with

$$V_E=B{H}_0+\left(A-B\right)H\left(x,t-\Delta t\right)+C{\epsilon }_{rk}\left(x,t-\Delta t\right),$$

(10)

$$A=\frac{\alpha D}{2e}\gamma \frac{J_k}{\Delta t}\left(1-e^{-\frac{\Delta t}{{\tau }_k}}\right),B=\frac{\alpha D}{2e}\gamma \frac{J_k}{\Delta t}{e}^{-\frac{\Delta t}{{\tau }_k}},C=\frac{e^{-\frac{\Delta t}{{\tau }_k}}}{T_k},$$

(11)

where $\gamma$ is the bulk weight, $\alpha$ is the constraint coefficient, $D$ is the pipe diameter, $e$ is the wall thickness, $H$ is the pressure head, and $\Delta t=L/\left(a\times N_x\right)$. Here, $N_x$ is the number of pipe meshes.

2.3. Numerical Scheme of 1D and Quasi-2D Models

The solutions for transient flows in viscoelastic pipes were calculated by quasi-2D and corresponding 1D models using the method of characteristics (MOC). The characteristic forms of the quasi-2D and corresponding 1D transient-flow equations are:

$$\frac{dH}{dt}\pm \frac{a}{g}\frac{du}{dt}+\frac{a^2}{gr}\frac{\partial \left(rv\right)}{\partial r}\mp \frac{a}{g}\frac{1}{\rho r}\frac{\partial \left(r\tau \right)}{\partial r}+\frac{2a^2}{g}{\displaystyle \sum }_{k=1}^N\left(AH-V_E\right)=0,$$

(12)

$$\frac{dH}{dt}\pm \frac{a}{g}\frac{dQ}{dt}\pm \frac{a^2}{g}\frac{4{\tau }_w}{\rho D}+\frac{2a^2}{g}{\displaystyle \sum }_{k=1}^N\left(AH-V_E\right)=0.$$

(13)

Because of the steep radial velocity gradient, a gradually encrypted radial grid was applied from the pipe centre to the pipe wall. This is shown in Figure 1, where j is the subscript representing the radial grid number, r_j is the radial distance between the outer surface of cylinder j and the pipe centre, and r_cj is the radial distance between the centre of cylinder j and the pipe centre.

Integrating Equation (11) on the characteristic lines between times nΔt and (n + 1)Δt (Figure 1) results in the following discretised expressions:

$$\left(1+F\right)H_i^{n+1}-\theta C_{q,j}{q}_{i,j-1}^{n+1}+\theta C_{q,j}{q}_{i,j}^{n+1}-\omega C_{u1,j}{u}_{i,j-1}^{n+1}+\left(\frac{a}{g}+\omega C_{u2,j}\right)u_{i,j}^{n+1}-\omega C_{u3,j}{u}_{i,j+1}^{n+1}=K_{pi,j},$$

(14)

$$\left(1+F\right)H_i^{n+1}-\theta C_{q,j}{q}_{i,j-1}^{n+1}+\theta C_{q,j}{q}_{i,j}^{n+1}+\omega C_{u1,j}{u}_{i,j-1}^{n+1}-\left(\frac{a}{g}+\omega C_{u2,j}\right)u_{i,j}^{n+1}+\omega C_{u3,j}{u}_{i,j+1}^{n+1}=K_{ni,j},$$

(15)

with

$$F=\frac{2a^2\Delta t}{g}{\displaystyle \sum }_{k=1}^N\frac{\alpha D}{2e}\gamma \frac{J_k}{\Delta t}{e}^{-\Delta t/{\tau }_k},$$

(16)

$$K_{pi,j}=H_{i-1}^n+\left(1-\theta \right)C_{qj}\left(q_{i-1,j-1}^n-q_{i-1,j}^n\right)+\left(1-\epsilon \right)C_{u1,j}{u}_{i-1,j-1}^n$$

$$+\left[a/g-\left(1-\epsilon \right)C_{u2,j}\right]u_{i-1,j}^n+\left(1-\epsilon \right)C_{u3,j}{u}_{i-1,j+1}^n+\frac{2a^2\Delta t}{g}{\displaystyle \sum }_{k=1}^N{V}_{E },$$

(17)

$$K_{ni,j}=H_{i+1}^n+\left(1-\theta \right)C_{qj}\left(q_{i+1,j-1}^n-q_{i+1,j}^n\right)-\left(1-\epsilon \right)C_{u1,j}{u}_{i+1,j-1}^n$$

$$-\left[\frac{a}{g}-\left(1-\epsilon \right)C_{u2,j}\right]u_{i+1,j}^n-\left(1-\epsilon \right)C_{u3,j}{u}_{i+1,j+1}^n+\frac{2a^2\Delta t}{g}{\displaystyle \sum }_{k=1}^N{V}_E ,$$

(18)

$$C_{qj}=\frac{a\Delta t^2}{g{r}_{cj}\left(r_j-r_{j-1}\right)}, C_{u1,j}=\frac{a\Delta t{v}_{Tj-1}{r}_{j-1}}{g{r}_{cj}\left(r_{cj}-r_{cj-1}\right)\left(r_j-r_{j-1}\right)} ,$$

$$C_{u3,j}=\frac{a\Delta t{v}_{Tj}{r}_j}{g{r}_{cj}\left(r_{cj+1}-r_{cj}\right)\left(r_j-r_{j-1}\right)}, C_{u2,j}=C_{u1,j}+C_{u3,j},$$

(19)

where ε and θ are the weighting coefficients, q(=rv) is the radial flux, and the source terms K_pi,j, and K_ni,j are known values, the elements of which depend on H, u, and q at the previous time level.

When Equation (14) is subtracted from Equation (15), we obtain:

$$\omega C_{u1,j}{u}_{i,j-1}^{n+1}-\left(\frac{a}{g}+\omega C_{u2,j}\right)u_{i,j}^{n+1}+\omega C_{u3,j}{u}_{i,j+1}^{n+1}=0.5\left(K_{ni,j}-K_{pi,j}\right).$$

(20)

Instead, if they are added, we obtain:

$$\left(1+F\right)H_i^{n+1}-\theta C_{q,j}{q}_{i,j-1}^{n+1}+\theta C_{q,j}{q}_{i,j}^{n+1}=0.5\left(K_{ni,j}+K_{pi,j}\right).$$

(21)

These equations enable us to calculate the axial velocity, nodal pressure head, and radial flux.

Integrating Equation (12) on the characteristic lines between times nΔt and (n + 1)Δt, the discretised forms of the 1D characteristic equations are expressed as:

$$\left(1+F\right)H_i^{n+1}+B{Q}_i^{n+1}+\frac{4\epsilon a\Delta t}{\rho gD}{\tau }_{w,i}^{n+1}=C_p,$$

(22)

$$\left(1+F\right)H_i^{n+1}-B{Q}_i^{n+1}-\frac{4\epsilon a\Delta t}{\rho gD}{\tau }_{w,i}^{n+1}=C_{m },$$

(23)

with

$${\tau }_{w,i}^{n+1}=\rho v_T\frac{u_{i,N_{R+1}}^{n+1}}{\Delta r_{c{N}_{r+1}}}=\rho v_T\frac{u_{i,N_{R+1}}^{n+1}}{\left(D/2-r_{c{N}_{r+1}}\right)}$$

(24)

$$C_p=H_{i-1}^n+\frac{a}{gA}{Q}_{i-1}^n-\frac{4a\Delta t}{\rho gD}{\tau }_{w,i-1}^n+\frac{2a^2\Delta t}{g}{\displaystyle \sum }_{k=1}^N{V}_E,$$

(25)

$$C_m=H_{i+1}^n-\frac{a}{gA}{Q}_{i+1}^n+\frac{4a\Delta t}{\rho gD}{\tau }_{w,i+1}^n+\frac{2a^2\Delta t}{g}{\displaystyle \sum }_{k=1}^N{V}_{E }.$$

(26)

A five-region turbulent flow model can be introduced to calculate the turbulent eddy viscosity ($v_T$). Combining Equations (21) and (22), the pressure head and discharge may be obtained as:

$$H_i^{n+1}=\frac{C_p+C_m}{2\left(1+F\right)},$$

(27)

$$Q_i^{n+1}=\frac{C_p-C_m}{2B}-\frac{4\epsilon a\Delta t}{\rho gDB}{\tau }_{wi}^{n+1}.$$

(28)

2.4. ITE Method Based on Quasi-2D Friction Model

The energy equation corresponding to the quasi-2D model is obtained by integrating the basic differential equation. However, in addition to axial flow, radial flow is considered in the quasi-2D model; therefore, radial integration should be considered when the equation associated with the flow rate is integrated [10].

$${{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^{D/2}\left[Eq. \left(3\right)\times {\rho }^{2}gH\right]drdx;\begin{array}{cc}& \end{array}{{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^{D/2}\left[Eq. \left(4\right)\times u\right]drdx$$

(29)

The energy equation of the quasi-2D model is as follows (a detailed derivation is provided in Appendix A):

$$\frac{dU}{dt}+\frac{dT}{dt}+D_f+W_E+W_P=0,$$

(30)

where U is the total internal energy of the system, T is the total kinetic energy of the system, D_f is the total rate of frictional dissipation, W_E is the total rate of work from the ends of the pipe, and W_P is the total rate of work from the pipe wall, expressed as follows:

$$U=\frac{\rho g^{2}A}{2a^2}{{\displaystyle \int }}_0^L{H}^2\left(x,t\right)dx,$$

(31)

$$T\left(t\right)=\frac{\pi \rho D}{4}{{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^{D/2}{u}^{2}drdx,$$

(32)

$$D_f\left(t\right)=-\frac{\pi D}{2}{{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^{D/2}\frac{u}{r}\frac{\partial \left(r\tau \right)}{\partial r}drdx,$$

(33)

$$W_E=\frac{\pi \rho Dg}{2}{{\displaystyle \int }}_0^{D/2}\left[u\left(L,t\right)H\left(L,t\right)-u\left(0,t\right)H\left(0,t\right)\right]dr,$$

(34)

$$W_P=2\rho Ag{{\displaystyle \int }}_0^{L}H\frac{\partial \epsilon }{\partial t}dx.$$

(35)

3. Results: Experimental Setup and Validation

The experimental data in this study were obtained from the Water Engineering Laboratory at the University of Perugia, Italy [19]. High-density polyethylene (HDPE) pipe was used in the experiment. A pressure tank was used to maintain the pressure upstream of the pipes, and a fast-closing pneumatic valve was used to generate excitation of the transient flow downstream of the pipes. The inner diameter (D) of the pipe was 93.3 mm, and the wall thickness (e) was 8.1 mm. The length of the experimental pipes (L) used in this study was 128.6 m.

Five groups of experimental data, each with a different Re, initial flow rate, and initial pressure, were analysed in this study. The specific experimental parameters are listed in

Table 1. Experimental settings [19].

Case	Q (L/s)	Re	H0 (m)	Tc (s)
1	1	13,380	21.63	0.0875
2	2.04	28,320	21.13	0.0752
3	2.95	40,950	20.74	0.1188
4	4.03	55,940	20.34	0.1575
5	5.02	69,680	19.82	0.1533

. The limiting coefficient of the pipe was 0.97, and the relative roughness of the pipe was 8.57 × 10⁻⁴. The K-V parameters were calibrated according to literature [19], and their values were as follows: 1D model (E₀ = 1556 N/mm², E₁ = 7820 N/mm², ${\tau }_1$ = 582.6 ms, E₂ = 18,370 N/mm², ${\tau }_2$ = 59.76 ms, E₃ = 6842 N/mm², and ${\tau }_3$ = 21,570 ms); 2D model (E₀ = 1563 N/mm², E₁ = 9596 N/mm², ${\tau }_1$ = 562.3 ms, E₂ = 19,490 N/mm², ${\tau }_2$ = 52.18 ms, E₃ = 5834 N/mm², and ${\tau }_3$ = 19,680 ms).

The numerical results of the 1D and quasi-2D models at Re = 55,940 are compared to the experimental results in Figure 2a. The figure shows that the simulation results of the quasi-2D model are very close to the experimental results in terms of peak value, valley value, and phase, especially in the second half-circle of transient flow (30–50 L/a); the same is true for the 1D model. Figure 2b shows the simulation results for the pressure head at different Re values, in which the peak pressure increases with Re. The pressure head decreases with the increase of Reynolds number, whereas the phase is unchanged at different Reynolds numbers.

Figure 3 shows the variation trend of the total, retarded, and instantaneous strain in the two cases. In the critical region of transient flow, the retarded strain accounted for a much smaller proportion of the total strain than the instantaneous strain. Although the retarded strain was relatively small, it can be seen from Figure 3 that, with an increase in L/a, the retarded strain decreased further and the instantaneous strain also showed a decreasing trend. At the later stage of transient flow (30 L/a), the total and instantaneous strain overlapped substantially.

From Figure 3, it can also be seen that the instantaneous strain reached its maximum value at the initial moment of the transient flow, whereas the retarded strain reached its maximum value in the third phase of transient flow (3 L/a); thus, the total strain reached its maximum value at 3 L/a. This suggests that the pressure fluctuation in viscoelastic pipes exhibits a delay compared with that in elastic pipes.

Figure 4 depicts the changes in the instantaneous and retarded strains at five Re values. Both tended to increase with Re. At the same time, there was a little effect on the peak value of the retarded strain at different Re values. The maximum peak values of the instantaneous and retarded strains changed almost linearly with Re. At the same Re value, the maximum instantaneous strain was approximately five times the maximum retarded strain.

Figure 5 illustrates the time difference ($\Delta t$) between the extreme value (peak no.) of retarded strain and instantaneous strain (the time of retarded strain/the time of instantaneous strain). According to the line graph, the time difference is the biggest at the first few peak numbers of the Re values, whereas it is closer at the end peak number. This means that different Reynolds numbers have greater influence on the first few cycles of strain, that is, viscoelasticity is obvious in the early phase of transient flow. During the whole period of transient flow, the time difference between retarded strain and instantaneous strain not only shows an upward trend with the peak number monotonic increasing but increases as Re values increase. It would be expected that the retarded strain effect would become smaller with time because the velocity decreases with time; and the difference of Re values also affects viscoelasticity during the transient flow, though not greatly.

To further illustrate the viscoelastic effect, the ratio of retarded strain to total strain (${\epsilon }_r/\epsilon$) is shown in Figure 6. As can be seen in the line graph, the ratio at the valley decreases from 0.12 in the first valley to 0.05 in the last valley, whereas there is a minor fluctuation in this peak numbers (from 0.175 in the first peak to 0.125 in the last peak, but the smallest value is 0.11 in the five peak). The ratios are consistent at different Re values, which means that the viscoelastic term has no significant effect on transient flow at different Reynolds numbers.

4. Discussion: Energy Analysis at Different Reynolds Numbers

To understand the influence of different Re values on the energy variations of viscoelastic pipes, distinctions between the viscoelasticity and friction terms were analysed.

4.1. D_f in 1D and 2D Models

The energy changes in the viscoelastic and friction terms of the 1D and 2D friction models under different initial Re values were further analysed based on the ITE method (i.e., Equation (29)). Combined with Equation (33), the simulation results of the work done by the friction term D_f are shown in Figure 7. The results show that larger D_f values correspond to greater values of the cumulative sum of the wall-shear stress of the entire pipeline. This indicates that the work done by the friction term becomes larger as Re values increase.

In Figure 7a,b, the D_f values calculated using the 2D model are larger than those calculated using the 1D model. It can be seen that the 1D model underestimated the instantaneous wall shear stress under flow conditions with initial Re values in the range of 1.0 × 10⁵–3.0 × 10⁵.

As shown in Figure 7, with an increase in the initial Re, D_f differed only slightly between the 1D and 2D models. This means that the instantaneous wall-shear stress calculated by the 1D model was close to that calculated by the 2D model, with an initial Re value of 4.0 × 10⁵–7.0 × 10⁵.

For all Re values, D_f was maximal near the initial time and gradually decayed with time in both the 2D and 1D models.

Furthermore, the D_f values were positive for both models in the critical region of turbulence. Thus, the D_f values show the dissipation of the friction energy in the transient flow. This is consistent with the description of transient flow energy in pipelines [10]. The D_f values also increased with an increasing initial Re.

4.2. W_P in 1D and 2D Models

The simulation results for W_P (the work done by the viscoelastic term of the pipe wall per unit time) are shown in Figure 8. The W_P values gradually increased with the initial Re values.

The W_P values from the 1D model were similar at all Re values to those from the 2D model because the governing equations of both models were consistent in calculating the retarded strain. The work due to the viscoelastic term of the pipe wall in the two friction models was essentially the same under different initial Re values. In particular, the maximum value of W_P calculated by both friction models appeared at approximately 2 L/a for all Re values, and the W_P values gradually decayed over time. In the critical region of turbulence, W_P had both positive and negative values. This indicates that the interaction between the fluid and the pipe wall during transient flow involves both energy transfer and energy dissipation, which agrees with the description of the energy variation of transient flows in pipelines in the literature [10].

By comparing Figure 7 and Figure 8, it can be seen that the work done by the viscoelastic and friction terms was similar at low Re. With an increase in the initial Re, D_f became significantly greater than W_P. This shows that the work done by the friction term had a significant influence on the energy dissipation when the initial Re values were relatively high [13].

To further illustrate the dependence of the viscoelastic effect at different flow conditions, the W_P values of viscoelastic term were compared with the 1D model at different Reynolds number (Figure 9). From these figures, the difference of W_P values between 1D model considering both friction and viscoelastic effects, and only considering the viscoelastic effect, increase with Re values, but compared with the increase of D_f values at different Reynolds number, there was a slight increase in W_P values in orders of magnitude. This is because the viscoelastic effect does not dominate under different initial Reynolds numbers [10,20].

4.3. Energy integral in 2D Models

Figure 10 shows the cumulative energies of the viscoelastic and friction terms, along with their sum (the total dissipational energy change) for the entire pipeline under different initial Re values.

The cumulative work done by the friction term increased with time, most steeply between 0 and 5 L/a; after 20 L/a, the energy dissipation gradually stabilised. This indicates that energy dissipation was the dominant mode of energy change.

The positive energy change of the viscoelastic term was the largest between 0 and 2.5 L/a, and the negative energy change was largest between 2.5 and 5 L/a, as shown in Figure 10. The energy change in the viscoelastic term had a sinusoidal fluctuation and an overall upward trend over time. This indicates that the energy change between the viscoelastic pipes and the fluid in the pipe not only led to energy conversion but also energy dissipation in the initial stage of the transient flow in the pipes. The energy changes became predominantly dissipative as time elapsed.

It can also be seen from Figure 10 that the energy changes of both the viscoelastic and friction terms increased with Re. At different initial Re values, the energy variation trend of the friction term was essentially the same. However, for the energy change of the viscoelastic term, the fluctuation range was larger and the duration longer when the initial Re values were small.

In particular, the proportion of energy dissipation generated by the friction term in the total energy change increased continuously with an increase in the initial Re values. Accordingly, the proportion of energy change generated by the viscoelastic term in the total energy change decreased. More specifically, the proportions of the viscoelastic terms in the total energy variation were 22.1% (Re = 13,880), 18.2% (Re = 28,320), 15.1% (Re = 40,950), 11.9% (Re = 55,940), and 10.1% (Re = 69,680) when the transient flow time was approximately 60 L/a.

5. Conclusions

In this study, the ITE method was used to derive the energy equation of a quasi-2D model of viscoelastic pipes, and the work changes of the friction term and viscoelastic term in 1D and 2D models under different initial Re values were compared. The results are as follows:

(1) When the initial Re was < 3.0 × 10⁵, the 1D (but not the 2D) model underestimated the work due to the friction term D_f. For Re > 3.0 × 10⁵, this error was reduced.

(2) The work done by the viscoelastic term of the pipe wall remained approximately constant under different initial Re values for both the 1D and 2D models.

(3) With an increase in Re values, both the viscoelastic and frictional work in the 1D and 2D models increased with time. However, D_f > W_P for a large initial Re.

(4) The energy dissipation of the friction term increased significantly initially then slowed down, and finally approached a constant value.

(5) The energy change of the viscoelastic term exhibited sinusoidal fluctuation during the early stage of the transient flow. At a smaller initial Re, the fluctuation in energy lasted for a long time, presented an overall upward trend, and finally approached a constant value.

(6) With an increase in initial Re values, the proportion of the energy dissipation generated by the friction term in the total energy change increased continuously, whereas the proportion of the energy change generated by the viscoelastic term in the total energy change decreased correspondingly.