Development Length of Fluids Modelled by the gPTT Constitutive Differential Equation

: In this work, we present a numerical study on the development length (the length from the channel inlet required for the velocity to reach 99% of its fully-developed value) of a pressure-driven viscoelastic ﬂuid ﬂow (between parallel plates) modelled by the generalised Phan–Thien and Tanner (gPTT) constitutive equation. The governing equations are solved using the ﬁnite-difference method, and, a thorough analysis on the effect of the model parameters α and β is presented. The numerical results showed that in the creeping ﬂow limit ( Re = 0), the development length for the velocity exhibits a non-monotonic behaviour. The development length increases with Wi . For low values of Wi , the highest value of the development length is obtained for α = β = 0.5; for high values of Wi , the highest value of the development length is obtained for α = β = 1.5. This work also considers the inﬂuence of the elasticity number. number seems to be neglected by the ﬂuid, since we obtain the same development length for El = 1 and El = 10. For El = 0.1, the results are quite different, which is due to the low values of Wi compared to the El = 1 and El = 10 cases. For β = 0.5, the rate of destruction of junctions increases and the development length decreases by about half. Again, the El = 1 and El = 10 cases show similar development lengths, which is due to the similar values of Wi and the almost creeping ﬂow conditions.


Introduction
A variety of functional applications are based on the premise that the flow is fully developed. It is assumed that after a certain time the fluid has travelled a certain length (development length-L) along the channel, after which the flow no longer changes in the direction of flow. This is used, for example, in extrusion dies, lab-on-a-ship, etc. [1][2][3].
The development length of Newtonian flows in channels and pipes (see Figure 1) is well understood [4]. Durst et al. [4] developed two correlations between L (the distance the fluid travels to become fully developed) and the Reynolds number, Re = ρUH η , where U is the imposed mean inlet velocity, ρ is the fluid density, H is the width of the channel (for pipes, one and predict well the development length for channel and pipe flows, respectively. For other works on the development length of Newtonian fluids please consult the following references [5][6][7][8][9][10]. For generalised Newtonian fluids (with varying viscosity), several works have been proposed in the literature [11][12][13][14][15][16][17][18][19]. We would like to highlight the works of Fernandes et al. [19] and Poole and Ridley [18], in which they presented two correlations for the development length in channel and pipe flows of power-law fluids (the viscosity is a function of the second invariant of the deformation tensor,γ (for simple flows,γ is simply the shear rate). The viscosity is then given by η = kγ n−1 ). The correlations are given by, for channel and pipe flows, respectively. Here, Re gen = 6ρU 2−n H n k n 4n+2 n , Re MR = 8ρU 2−n D n k n 6n+2 n , and f (n) = −0.355 1+2 exp (0.553−4.273n) . Note the increasing complexity in the correlations when going from a Newtonian to a power-law fluid.
In the case of viscoelastic fluids, the number of papers on this topic is smaller. This is due to the complexity of viscoelastic flows, such as the presence of singularities at the entrance of the channel, overshoots in the velocity profile, and the high Weissenberg number problem.
We would like to highlight the work of Na and Yoo [20] in which they perform numerical simulations to determine the development length of an Oldroyd-B fluid and conclude that the development length (for a fixed Re) increases slightly with the Weissenberg number, Wi = λU H (where λ is the relaxation time of the fluid in seconds), but is more strongly affected by Re. Liang [1] proposed a theoretical work for the development length of viscoelastic fluids entering an extrusion die. They presented an expression for estimating the length of the entrance region, which has applications in the extrusion industry. In the work by Philippou et al. [10] the authors present a study on the flow development of a Bingham plastic fluid in tubes and channels considering the Papanastasiou regularisation and the finite element method. They considered the Navier's slip law at the wall and concluded that as slip increases, the development length initially increases exhibiting a global maximum before vanishing rapidly above the critical point corresponding to sliding flow. More recently, Yapici et al. [21] presented a study on the development length of steady flows of Oldroyd-B and Phan-Thien-Tanner (PTT) fluids through a two-dimensional rectangular channel and concluded that the development length determined for the Oldroyd-B fluid varies exponentially with Wi and linearly for the linear PTT model; they also concluded that higher entry lengths are predicted with increasing Wi (at fixed Re).
To remove the unstable numerical effect of the singularity at the entrance corner, a continuous inlet velocity profile is used in both works of Na et al. [20] and Yapici et al. [21]. This regularised profile can affect the true development length, so in the work of Guilherme [22] the log-conformation formulation [23][24][25] is used, which reduces the rate of increase of the stresses and thus avoids the need to introduce artificial inlet velocity profiles.
It should be noted that the development of a correlation for the prediction of the development length of such complex fluids is still difficult due to the high number of parameters involved and the fact that it is model dependent.
Here we follow the work of Guilherme [22], where a detailed analysis of the development length of the linear PTT model is performed. We extend his work to the exponential and generalised PTT models [22,26,27].
This work is organised as follows. First, we present the differential equations to be solved and their numerical solution. In Section 3, we present the geometry and the meshes. In Section 4, we perform a validation of the numerical method and the meshes, using Newtonian benchmark results. In Section 5, we present and discuss the results for viscoelastic fluids. The paper ends with the main conclusions in Section 6.

Governing Equations
The equations governing the flow of an incompressible fluid, under isothermal conditions, are the continuity, ∇ · u = 0 (5) and the momentum equations, where u is the velocity vector, p is the pressure, σ is the stress tensor (to be defined later), ρ is the mass density, and F represents the external forces. Note that all variables are dimensionless, with: x = x * /H, u = u * /U, t = t * U/H, p = p * /(ρU 2 ), σ = σ * /(ρU 2 ) (the * represents the dimensional variable).
In order to achieve a closed system of equations, a constitutive equation for the extrastress tensor, σ, is required. Recently, Ferrás et al. [27] proposed a new differential model based on the Phan-Thien-Tanner constitutive equation [26] (see also [28]). This new model considers a more general function for the rate of destruction of junctions, the Mittag-Leffler function, where one or two fitting parameters are included, in order to achieve additional fitting flexibility.
The Mittag-Leffler function is defined as, with α, β being real and positive. Γ(·) is the gamma function, given by: when α = β = 1, the Mittag-Leffler function reduces to the exponential function, and when β = 1 the original one-parameter Mittag-Leffler function, E α is obtained. The constitutive equation is given by: where σ kk is the trace of the stress tensor, Wi = λU/H is the Weissenberg number, Re = UH/ν is the Reynolds number (ν = µ 0 /ρ is the kinematic viscosity), D = 1 2 ∇u + (∇u) t is the rate of deformation tensor, σ is the elastic stress, and ζ = µ S µ 0 is the viscosity coefficient, where µ 0 = µ S + µ P is the total shear viscosity (µ S is the solvent/Newtonian viscosity, µ P is the polymer viscosity) and σ represents the Gordon-Schowalter derivative.
The stress function, K(σ kk ), is given by a new formulation that imparts more flexibility and accuracy to the model predictions, as discussed in [27,29,30]. It is given by: where ε represents the extensibility parameter, Γ is the Gamma function, and the normalisation Γ(β) is used to ensure that K(0) = 1, for all choices of β.

Numerical Method
The numerical method used in this work is based on finite differences. It can deal with tree-like mesh grids (see Figure 2b) and allows fast Cartesian discretizations, flexibility and accuracy, and local mesh refinement. In order to fit the discretization stencil near the interfaces between grid elements of different sizes, a robust method based on a moving least squares meshless interpolation technique is used to compute the weights of the finite difference approximation in a given hierarchical grid, allowing complex mesh configurations and preserving the overall accuracy of the resulting method. Figure 2a shows a schematic representation of the mesh refinement. Note that some of the points (variables) of the computational cells (red dots) must be approximated because they are not located in the center of the cell (e.g., the center of computational cell 1 is not the same as the location of the red dot used to compute the derivative of the property being evaluated). To solve this problem, we use a special adaptive least square interpolation (MLS). The method is known as HiG-Flow, and a detailed mathematical explanation can be found in [31,32]. For the numerical solution of the Navier-Stokes equations together with the constitutive equation given by the gPTT model, the momentum Equation (6) is rewritten: From Equation (9), the rheological constitutive equation can be written as: where M(σ) is given by Equation (14), and the parameter ξ (0 ≤ ξ ≤ 1) accounts for the slip between the molecular network and the continuous medium. For ξ = 0 there is no slip and the motion becomes affine.
The parameter ξ leads to a non-zero second normal-stress difference in shear, resulting in secondary flows in ducts having non-circular cross-sections. Since in this work we only consider 2D flows, and, due to the high number of parameters involved in the numerical simulations and the gPTT model itself, we have only considered the case when ξ = 0.
3.1. Calculation of u (n+1) and p (n+1) To calculate the velocity u (n+1) and pressure p (n+1) fields, we use the incremental projection method by Chorin [33], that uncouples the mass conservation and momentum equations, given by Equations (5) and (6), respectively. This method allows one to obtain an intermediate velocity field u n+1 from Equation (11). In the HiG-Flow methodology, this Equation (11) can be approximated using an explicit Euler method, Runge-Kutta RK-2 or RK-4, or, the semi-implicit Euler methods, Cranck-Nicolson, and BDF2. One can also choose a spatial discretization orders of 2 or 4. One can use the the convective central schemes or Upwind (order 1), or, schemes of order 2 like the Cubista [34] and Quick [35].
In this work an Implicit Euler scheme together with a second order spatial approximation and an Upwind Cubista scheme for the convective terms, is used: here, δt is the time step, n represents the known values of velocity, stress, and pressure at instant n, and n + 1 represents the new velocity field values (unknown) to be obtained from the solution of the equation. At the inlet, (see Figure 3) we consider a constant velocity profile, u(y) = 1 (the stress components are set to 0) and at the outlet, we assume fully developed boundary conditions (Neumann boundary conditions) for the velocity and stress (the pressure is imposed). Finally, at the walls (y = 0 and y = 1), we have the empirical no-slip boundary condition (u = 0). Using the projection method, it is well known that the velocity field u n+1 obtained from Equation (15) may not satisfy the mass conservation equation. Therefore, in order to solve this problem, the equation for the potential ψ (n+1) = δt(p n − p (n+1) ) is solved, and the Helmholtz-Hodge decomposition (see [31,32,36] for more details) is used to correct the previous non-conservative velocity field u (n+1) , The new velocity field u (n+1) satisfies the mass conservation equation. Finally, the pressure is updated p (n+1) = p n + ψ (n+1) δt .

Calculation of the Extra-Stress Tensor σ (n+1)
The velocity and pressure fields were obtained in the previous subsection. We now aim to obtain the extra-stress tensor σ (n+1) field. Equation (13) is solved using the Explicit Euler method, and, to calculate M(σ) (see Equation (14)), the Mittag-Leffler function and the term Γ(β)E α,β are obtained numerically from Equation (18) and the approximations presented in the work by R. Gorenflo, J. Loutchko, and Y. Luchko [37], The numerical implementation of the Mittag-Leffler function follows the work by Davide Verotta and Eduardo Mendes [38] (developed in Fortran), which is adapted in this work to C++. The original fortran code is based on a Matlab function developed by Igor Podlubny and Martin Kacena [39] which, in turn, was based on the reference [37].

Geometry
Due to the low Re values considered in this work, and based on the few literature results on the development length of viscoelastic fluids, we considered a geometry where the length of the channel is fixed at 10 times its width ( Figure 3).

Meshes
We performed simulations considering more than 8 levels of mesh refinement. After some numerical experiments, the following meshes were considered:  Numerical simulations were also performed considering a mesh with additional refinement in the centerline of the channel, as shown in Figure 5. A total number of 16,000 cells was used.

Validation of the Numerical Method
The validation of the numerical method is performed in two steps. First, the numericallydetermined fully-developed velocities and stresses are compared with the analytical solution [27]. Then, the development length obtained for the gPTT model with Wi = 0.001 (almost Newtonian fluid) is compared with the benchmark results of Durst et al. [4]. (a) α = 0.5 and β = 0.5 (b) α = 1.5 and β = 1.5  It can be seen that an excellent agreement is obtained between the analytical and numerical solutions for all the considered values of Wi, which underlines the robustness of the numerical method and the meshes.

Comparison with the Analytical Solution
For lower values of α and β, we obtain a higher destruction rate of the junctions in the gPTT model. Note that the typical viscoelastic velocity profile is more flattened for lower values of α and β. In this case, the different values of Wi have a stronger impact on the model's behaviour. This result is similar to those found in the literature comparing linear and exponential functions of the trace of the stress tensor.

Comparison with the Development Length of a Newtonian Fluid
In the limiting case of Wi → 0 we obtain a Newtonian fluid. Therefore, we considered Wi = 0.001 and performed simulations for the development length of a gPTT fluid, using the geometry shown in Figure 4. We considered a Reynolds number in the range [0, 100], where the nonlinear variation of the development length with Re is more pronounced. The other parameters of the model were set as follows: Figure 7a shows a comparison between the development length obtained with the gPTT model, a Newtonian fluid, and that obtained by the Durst et al. [4] correlation for the variation of the development length with Re (see Equation (1)). The three results practically overlap, proving once again the robustness of the numerical method.
As Re increases, the results for the gPTT model in the coarse mesh are slightly higher than those obtained for the Newtonian fluid and the correlation. However, in the nonlinear domain the results are quite accurate. Figure 7b shows the velocity profiles obtained in the fully developed region of the channel (mesh M 3 ) considering the gPTT and Newtonian models for Re = 0.001 and Re = 100. 9 of 19 (a) (b) 1 10 100 Newtonian -Re = 10 2 gPTT -Re = 10 2   Again, there is excellent agreement between the two solutions for the two different values of Re. This shows that the value of Wi = 0.001 is a good approximation for the Newtonian fluid.
Based on these results, the numerical code is now able to predict the development length of the fluid modelled by the gPTT model considering a wider range of Wi numbers. This gives a total number of 132 simulations. The simulations with the finest mesh took about 15 h each.

Development Length of a gPTT Fluid
The first set of 72 simulations allowed conclusions to be drawn about the convergence of the numerical method and the error in calculating the development length using the Richardson extrapolation technique. Based on the results of these simulations, a second set of simulations was performed for higher values of Wi using meshes M 2 (see Figure 4) and M r (see Figure 5). These meshes were chosen based on a trade-off between accuracy and computational time. The development length, determined as the length from the channel inlet required for the velocity to reach 99% of its fully developed value, and denoted here as L 99% , is shown in Table 1. The results are shown only for Wi up to 0.4, since convergence problems for finer meshes are observed for higher values of Wi. The main problem arises from the singularity at the entrance corner of the channel, which generates an error that propagates along the channel (for more details, see [24,25]). Note that the error is higher at the lowest and highest values of α and β, being more pronounced when α and β are low. The maximum error was 3.6% and was obtained, as expected, for Wi = 0.4 and α = β = 0.5. It should be noted that the errors are quite low and therefore these solutions can be used as benchmarks. Figure 8 shows the development lengths for mesh M 2 presented in Table 1. A nonlinear variation of the development length with α, β, and Wi is observed. these parameters on the development length is small. These results can be justified by the low value  259 Therefore, to capture the essence of the development length for higher Wi values, we considered 260 another development length, L 98% (length from channel entry required for the velocity to reach 98% 261 of its fully-developed value), which is less restrictive.

262
The results obtained for L 98% are shown in Figure 9 for mesh M 2 and Table 2 for the three meshes As Wi increases. the highest value of development length is reached with a low rate of the junctions decreases. the information is transmitted more slowly due to the small number of new The development length increases with Wi, with viscoelastic effects delaying the diffusion and convection of information from the walls to the center of the channel. This diffusion and convection is also strongly influenced by the parameters of the Mittag-Leffler function. For low values of Wi, the highest value of the development length is obtained for α = β = 0.5; for high values of Wi, the highest value of the development length is obtained for α = β = 1.5. A molecular continuum explanation of this phenomenon is not an easy task. At high α and β values, the rate of destruction of the junctions is lower than at low α and β values. This means that when Wi values are low and the rate of junction destruction is high, information travels slowly from the wall to the center of the channel (compared to when the rate of junction destruction is low). The opposite was expected. Note that in this case the development lengths are very similar for all tested values of α and β, and therefore the influence of these parameters on the development length is small. These results can be justified by the low value of Wi.
As Wi increases, the highest value of development length is reached with a low rate of destruction of junctions. This result can be justified by the fact that as the rate of destruction of the junctions decreases, the information is transmitted more slowly due to the small number of new contacts between the strands representing the molecules.

Case of L 98%
In addition to the error arising at the entrance corner due to a singularity, we also have the problem of the development length, which takes into account 99% of the fully developed maximum velocity and may not work so well in an intermediate mesh as M 2 , for higher Wi. This leads to an increased difficulty for the numerical method to capture L 99% for the mesh M 2 .
Therefore, to capture the essence of the development length for higher Wi values, we considered another development length, L 98% (length from channel entry required for the velocity to reach 98% of its fully-developed value), which is less restrictive.
The results obtained for L 98% are shown in Figure 9 for mesh M 2 and Table 2 for the three meshes (along with the extrapolated development length value).  The results obtained follow the same trend observed in the case L 99% , with some minor differences. The development length is initially higher for the case α = β = 0.5 than for the case α = β = 1, until Wi = 0.2, where the growth rate of L 98% becomes larger with Wi for α = β = 1. For L 99% , the two development lengths are quite similar.
For Wi = 1, we obtain a development length of 2.679 for α = β = 1.5 and a development length of 1.313 for α = 1.5, β = 0.5 (and 1.368 for α = β = 0.5). Again, these results are consistent with the idea that the higher the rate of destruction of junctions, the smaller the development length (information travels faster across the channel).    Figure 10 shows the reference velocity profiles obtained for the classical case of an exponential PTT model (α = β = 1). The velocity profile evolves from a plug profile in the center of the channel (for profiles near the inlet) to the typical parabolic profile (in the fully developed region). Note the overshoots near the walls that occur when the fluid is still developing. This is due to the different characteristic times of the fluid and the diffusion of information moving from the walls (y/H = 0 and y/H = 1) to the center of the channel (y/H = 0.5). Figure 11 shows the velocity profiles obtained for α = 0.5, β = 0.5 and α = 0.5, β = 1.5. It can be seen that the overshoots are stronger near the inlet (compared to the exponential PTT model) and that a lower maximum velocity is obtained when the fluid is fully developed. The influence of the parameter β is residual.  Figure 12 shows the velocity profiles obtained for α = 1.5, β = 0.5 and α = 1.5, β = 1.5. In Figure   259 12a it can be seen that the overshoots near the inlet resemble the case of the exponential PTT model, 260 and that the maximum velocity increases again. The main difference is that the plug profile is now  Figure 12 shows the velocity profiles obtained for α = 1.5, β = 0.5 and α = 1.5, β = 1.5. In Figure 12a it can be seen that the overshoots near the inlet resemble the case of the exponential PTT model, and that the maximum velocity increases again. The main difference is that the plug profile is now less pronounced. When β increases (Figure 12b), one can observe a dramatic change in the evolution of the velocity profiles. The velocity overshoots are almost suppressed and the maximum velocity increases. This means that the rate of destruction of junctions improves the diffusion of information.

The Influence of the Elasticity Number
In this subsection we study the influence of the elasticity number, El = Wi Re , on the development length of the gPTT model. We consider three different values for El (0.1, 1, 10), α = 0.5, and two different values for β, 0.5 and 1.5. The mesh used for the simulation is M 2 .
The results are shown in Tables 3 and 4 and Figure 13 for L 98% and L 99% . As expected, the results for low values of Re are consistent with those obtained earlier in this work (see previous sections for more details). The results obtained for the different definitions of the development length are qualitatively similar. However, higher values of development length are obtained for L 99% , as expected.
For high values of β (1.5), the influence of the elasticity number seems to be neglected by the fluid, since we obtain the same development length for El = 1 and El = 10. For El = 0.1 the results are quite different, which is due to the low values of Wi compared to the El = 1 and El = 10 cases.
For β = 0.5, the rate of destruction of junctions increases and the development length decreases by about half. Again, the El = 1 and El = 10 cases show similar development lengths, which is due to the similar values of Wi and the almost creeping flow conditions.
The results show that for the tested ranges of El, Wi, and Re, no critical value is found for El. This is due to the fact that Mach's Elastic number is less than 1.

282
In this work, we present a numerical study on the development length of a pressure-driven  The numerical results showed that the in the creeping flow limit (i.e., Re = 0), the development

Conclusions
In this work, we present a numerical study on the development length of a pressuredriven viscoelastic fluid flow (between parallel plates) modelled by the generalised Phan-Thien and Tanner (gPTT) constitutive equation. The governing equations are solved using the finite-difference method, and, a thorough analysis on the effect of the model parameters α and β is presented. We consider two different definition of the development length: The length from the channel inlet required for the velocity to reach 99% (and 98%) of its fully-developed value.
The numerical results showed that the in the creeping flow limit (i.e., Re = 0), the development length for the velocity exhibits a non-monotonic behaviour. The development length increases with Wi, with viscoelastic effects delaying the diffusion and convection of information from the walls to the center of the channel. For low values of Wi, the highest value of the development length is obtained for α = β = 0.5; for high values of Wi, the highest value of the development length is obtained for α = β = 1.5. At high α and β values, the rate of destruction of the junctions is lower than at low α and β values. This means that when Wi values are low and the rate of junction destruction is high, information travels slowly from the wall to the center of the channel (compared to when the rate of junction destruction is low). The opposite was expected. Note that in this case, the development lengths are very similar for all tested values of α and β, and therefore the influence of these parameters on the development length is small. These results can be justified by the low value of Wi.
As Wi increases, the highest value of development length is reached with a low rate of destruction of junctions. This result can be justified by the fact that as the rate of destruction of the junctions decreases, the information is transmitted more slowly due to the small number of new contacts between the strands representing the molecules.
We also studied the influence of the elasticity number, El, on the development length of the gPTT model. As expected, the results for low values of Re are consistent with those obtained for creeping flow. For high values of β (1.5), the influence of the elasticity number seems to be neglected by the fluid, since we obtain the same development length for El = 1 and El = 10. For El = 0.1, the results are quite different, which is due to the low values of Wi compared to the El = 1 and El = 10 cases. For β = 0.5, the rate of destruction of junctions increases and the development length decreases by about half. Again, the El = 1 and El = 10 cases show similar development lengths, which is due to the similar values of Wi and the almost creeping flow conditions.