Benchmark solutions of the stabilized computations of flows of fluids governed by the Rolie-Poly constitutive model

Recent studies demonstrate that flow induced non-uniformities of concentration can trigger shear banding in the flow of certain viscoelastic fluids. These studies show that the driving mechanisms for such shear banding are related to the coupling of the polymer stresses to an inhomogeneous concentration profile. The Rolie-Poly (RP) viscoelastic constitutive model has been used in such studies since it has been comprehensively subjected to extensive experimental validation with regards to shear banding and has the demonstrated ability to accurately express the rheology of polymer solutions for a wide range of strain rates. The primary aim of this work is to develop an efficient computational methodology that could be used to accurately simulate the flow of complex fluids governed by the Rolie-Poly constitutive equation. The development of such a computational platform is crucially important for the purposes of our follow up studies on the computational analysis of shear banding phenomena by coupling polymer stress with inhomogeneous concentration profile. Our numerical algorithms will be based on the finite volume method (FVM) and will be implemented on the open source software package OpenFOAM®. In this paper, we will present both validation results as well as new benchmark results from our FVM based OpenFOAM® numerical solver for flow of fluids governed by the Rolie-Poly constitutive model. We use two well-known benchmark problems, the lid-driven cavity flow and the 4:1 planar contraction flow problems. In order to stabilize the numerical algorithm at high Weissenberg numbers, we employ either of two stabilization techniques; the Discrete Elastic Viscous Stress Splitting (DEVSS) technique as well as the Log-Conformation Reformulation (LCR) methodology. Validation of our results is done by comparing our (stabilized) numerical results, against data from existing literature. The numerical results obtained for the contraction flow using the LCR stabilization approach are in good agreement with the existing literature for a wider range of Weissenberg numbers. The DEVSS method shows a good agreement only for lower Weissenberg numbers. For the lid-driven cavity flow, good agreement with the existing literature is observed for low Weissenberg numbers using either of the two stabilization techniques.


Introduction
The recent studies by Cromer et al [1] demonstrate that flow induced non-uniformities of concentration can trigger shear banding in shear flow of fluids governed by the Rolie-Poly constitutive model. In this work we develop and test a computational methodology that could be used to accurately simulate the flow of complex fluids governed by the Rolie-Poly constitutive equation.
There are a number of benchmark flow problems that are often used as test cases when developing new numerical methods in computational rheology, primarily for validating the new numerical methodologies. Typical benchmark flows used include the lid-driven cavity flow, flows through axisymmetric and planar contractions and flows around cylinders or spheres. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
One advantage of the lid-driven cavity flow is that the geometry is simple and thus easily yields to numerical methodology. Notwithstanding its simple geometry, complex flow structures such as eddies and vortices can be visualized. Basically, a fluid is contained in a rectangular geometry whose top wall is translated horizontally at some specified velocity. Fluid motion is subsequently generated, similarly to wall driven shear flows. Discontinuity in the boundary conditions at the two top corners where the side walls meet the lid are inevitable, leading to the so called corner singularities [2]. The experimental and numerical results for the planar contraction flow, in particular the 4:1 contraction flow problem, have been extensively documented in the literature for both Newtonian and viscoelastic fluids. The abundance of such flow results in the literature gives an obvious advantage to using the planar contraction flows as benchmarks for newly developed numerical methodologies on the simulation of, say, the flow of viscoelastic fluids.
In order to eliminate the effects of corner singularities, different techniques have been employed. One such technique involves incorporating a controlled amount of leakage at the upstream and downstream corners of the cavity as given in Grillet et al [3]. Another method involves regularizing the lid velocity. This is done by imposing a parabolic velocity distribution along the moving lid so as to have a vanishing velocity and velocity gradient at the two upper corners as given in the works of Poole et al [4], Habla et al [5] and Comminal et al [6] among others.
For contraction flow benchmark problems, the focus is on predicting the formation, development and dynamics of the corner and lip vortices. The vortex lengths are then used in the comparison of numerical and/or experimental results. However, the numerical simulation of such flows of viscoelastic fluids present a huge challenge, in particular, the High Weissenberg Number Problem (HWNP). The HWNP is the breakdown of numerical schemes beyond some critical value of the Deborah (or Weissenberg) number [7]. The Deborah number, determines the elastic character of the flow and is given as the ratio of the polymer relaxation time to the time scale of the flow. It has been observed that when the Deborah number exceeds some critical value, numerical methods for viscoelastic fluid flow computations break down. This critical value is dependent on the viscoelastic constitutive model, the flow type, the numerical method and mesh used.
The HWNP is a result of steep exponential profiles which arise as due to the viscoelastic stresses experiencing a combination of deformation and convection. These exponential profiles are caused by the inappropriateness of polynomial based approximations that are used to represent the viscoelastic stress tensor. In particular, such polynomial based approximations are exponential in regions of high deformation rates, or near stagnation points [8].
The Discrete Elastic Viscous Stress Splitting (DEVSS) and LCR techniques are some of the methods, among many others, that have been developed to mitigate against the HWNP. The DEVSS method involves reexpressing the constitutive equation to include an explicit viscous stress and introducing the velocity gradient as an additional variable [9]. A variant of such method is the Elastic Viscous Stress Splitting (EVSS) technique [10]. Other methods aim to stabilise the advection terms in the constitutive equations. Such methods include, the upwinding techniques for convective equations, for example the Discontinuous Galerkin (DG) method [11] and the Streamline Upwind Petrov Galerkin (SUPG) method [12].
The LCR approach for numerical stabilization was recently proposed by Fattal and Kupferman [13]. This technique reformulates the constitutive law in terms of a matrix logarithm and improves the representation of large stress gradients by linearizing the stress profile. The LCR technique introduces a better polynomial interpolation via logarithmic variables and preserves the positive-definiteness of the conformation tensor. It is worth noting that the loss of the positive-definiteness of the stress tensor is a precursor to the HWNP.
The LCR method has been implemented in various studies such as in the finite volume method (FVM) framework for creeping flows of viscoelastic fluids in steady and unsteady flows around a confined cylinder [14], the numerical simulation of viscoelastic flow in three-dimensional lid-driven cavity flow in OpenFoam® [15] and in the robust simulations of viscoelastic flows at high Weissenberg numbers [6]. The LCR representation has also been used with the finite element method for simulating lid-driven cavity flow of Oldroyd-B fluids at high Weissenberg numbers [16].
A combination of the LCR and pure stream-function formulations was used in [17] for the two-dimensional flow of an Oldroyd-B fluid inside a lid-driven cavity simulated over a wide range of Weissenberg numbers while [18] employed a combination of LCR representation and an operator splitting lie scheme to simulate timedependent cavity flow of an Oldroyd-B fluid using a regularized velocity to remove the discontinuities at the two upper corners.
In this work, both the lid-driven cavity flow and 4:1 planar contraction flow will be used to benchmark the Rolie-Poly viscoelastic flow solver that we have developed. The solver is implemented on the open-source, OpenFoam®platform. Some of the advantages of using OpenFoam®are that, as an open source platform, it has no limiting aspects regarding licensing fees and that it has the ability to flexibly deal effectively with complex geometries. The OpenFoam®software is based on the finite volume method and has the ability of parallelization. OpenFoam®uses the C++ object oriented programming, making it convenient for users to incorporate their own models since existing solvers can be modified. For example, Favero et al [19] developed and incorporated a viscoelastic fluid flow solver into OpenFoam®which was later extended by Florian et al to handle two-phase flows of viscoelastic fluids. Silva and Lage [20] extended the two-phase flow solver to include a multi-phase flow formulation.
The purpose of the present work is to develop a viscoelastic flow solver for fluids governed by the Rolie-Poly constitutive model. The solver implements both the LCR Reformulation and DEVSS techniques. We assess the efficiency and accuracy of the solver by comparing the simulation results, for values of the material parameters corresponding to the widely studied Oldroyd-B fluid model, with those Oldroyd-B results in the existing literature. We then present benchmark simulation results and predictions obtained for the full Rolie-Poly constitutive model using the lid-driven cavity and 4:1 planar contraction benchmark flows.
The momentum equations read, We will use the Rolie-Poly constitutive equation of Likhtman and Graham [21] to model the dynamics of the viscoelastic stress t . At present, the Rolie-Poly model remains the most advanced differential constitutive formulation of the Doi-Edwards tube models for linear entangled polymer melts and it includes the processes of reptation, convective and reptation-driven constraint release, chain stretch and contour length fluctuation. The Rolie-Poly constitutive equation may be written as, where λ is the reptation relaxation time, λ R is the Rouse relaxation time, β * is the convective-constraint release (CCR) parameter and often varies from 0 to 1. I is the identity tensor, †denotes the matrix transpose and tr denotes the trace. The parameter δ is obtained from experimental data and takes the value of −1/2. We note that when l  ¥ R the Rolie-Poly model reduces to the Oldroyd-B model.

DEVSS approach
The DEVSS technique is used to improve the numerical stability by introducing an additional diffusion term on each side of the momentum equations [22]. The momentum equations (2) are then rewritten as, where k is a positive number that is related to the parameters of the constitutive model and, according to Jovani [22], a good choice is k=η p where η p is the polymer viscosity.

Log-conformation reformulation (LCR) approach
To improve the numerical stability, the log transformation known as the Log-Conformation Reformulation (LCR) approach is used. The approach consists in a change of variable for the polymeric extra stress which is related to the conformation tensor (c) by the equation where λ is the relaxation time and η p is the polymer viscosity. The viscoelastic constitutive equation can be written in terms of the conformation tensor as where ( ) c f R is the relaxation function which is a polynomial of the conformation tensor. The process outlined in this section is as suggested by Fattal and Kupferman [13]. Since the conformation tensor c is a positive-definite matrix, it can be diagonalized according to where Λ is a diagonal matrix consisting of the eigenvalues of c and where the matrix R is orthogonal and contains the corresponding eigenvectors. Instead of solving the equation (3), for the stress t , this equation is reformulated in terms of the natural logarithm of the conformation tensor c, T Y then becomes our new variable and the velocity gradient U is decomposed as 9 1 where Ω and N , which both account for rotations are anti-symmetric tensors. B is a diagonal tensor which accounts for pure extensions. Substituting equation (8) and velocity gradient decomposition 9 into equation (6) yields the log-conformation equation, where B and Ω are pure extension (symmetric, traceless) and pure rotation (anti-symmetric) matrices respectively that are obtained from the projection of the velocity gradient into the base of the stress tensor. The eigen decomposition of the conformation tensor in a two-dimensional flow is given as where λ 1 and λ 2 are the eigenvalues and R the orthogonal matrix containing the eigenvectors. The change of the base of velocity gradient is given as The pure extension and rotation matrices are obtained as, where ζ=(m 12 +m 21 )/(λ 2 −λ 1 ). For the first time step, Ω is set to be zero and c I leads to λ 2 =λ 1 would result in an undefined division by zero. The conformation tensor c is recovered from the matrix exponential of Y, once equation (10) has been solved. To verify the positive-definiteness of the conformation tensor ( ) c det must be greater than zero.

Solution algorithm and numerical method
We employ the finite volume method implemented on the open source software package OpenFoam®. The algorithm used to calculate the pressure field is the semi-implicit method for the pressure linked equation consistent (SIMPLEC) which is a modification of the SIMPLE algorithm since the SIMPLEC algorithm has been proved to have better convergence properties than the SIMPLE algorithm [23,24]. The Convergent and  Universally Bounded Interpolation Scheme for the Treatment of Advection (Cubista) scheme, a high resolution scheme with improved iterative convergence properties, is used to discretize the convective terms and temporal derivatives appearing in the momentum and transport equations. The Cubista scheme was devised by incorporating total-variation diminishing constraints, appropriate for unsteady problems, into an implicit timemarching method used for steady flow problems [25].
The computational steps followed by the solver are: 1. Initialize the variables; 2. For the LCR approach, solve for Y via equation (10), or alternatively;    The solver for the stress constitutive equations as well as for the momentum equations is the bi-conjugate gradient stabilized (BiCGstab) solver with a Cholesky preconditioner. A preconditioned conjugate gradient (PCG) solver is used to solve the pressure equation in conjunction with a simplified diagonal-based incomplete LU (DILU) preconditioner.

Geometry and boundary conditions
For the lid-driven cavity flow, see figure 1, a square cavity is used with the upper lid moving with a regularized velocity. The flow is assumed to be two dimensional (2D), being limited to the xy-plane. No-slip boundary and  For validation and becnmarking purposes of our Rolie-Poly solver, we take advantage of the fact that the Rolie-Poly model reduces to the Oldroyd-B model when l  ¥ R as well as the fact that benchmark Oldroyd-B results widely exist in the literature. We therefore validate and benchmark our solver by recovering the Oldroyd-B results from our Rolie-Poly solver and comparing these with existing literature on Oldroyd-B lid-driven cavity, say [6]. The maximum lid velocity is set to u=1. A retardation ratio β of 1/2 is used and the Reynolds number is set to negligible value of Re=ρuL/η 0 =5×10 −4 which is considered as creeping flow. The velocity components and the components of the log-conformation tensor are obtained. As in the existing literature, the three components of the tensor, the velocity and log-conformation profiles along vertical line x=1/2 and along horizontal line y=3/4 including the history of the specific Kinetic energy E k were considered.
To validate the DEVSS method, we only consider the velocity profile, say, along lines x=1/2 and y=3/4 and history of specific kinetic energy since there are no log-conformation tensors.
Our obtained results were all in satisfactory agreement with those in the existing literature, say [6]. A sample summary of our graphical results are shown in figures 2-10.
Given the satisfactory agreement between our results and the benchmarking results from the literature for the lid-driven cavity geometry, we therefore proceed to present the lid-driven cavity results for the full Rolie-Poly constitutive model as predicted and simulated by our solver.

Lid-driven cavity results for the full Rolie-Poly model
In this section we present the results for the numerical simulation for the lid-driven cavity flow for the full Rolie-Poly constitutive model with l < < ¥ 0 R . As in the previous sections, we will also make comparisons between the DEVSS and the LCR stabilization methods. The velocity components and the components of the logconformation tensor are computed and plotted. The components of the tensor, the velocity and logconformation profiles along vertical line x=1/2 and along horizontal line y=3/4 are presented in figure 11. Figure 11 shows no apparent difference between the LCR method and the DEVSS results for Deborah numbers 1 and 2. There is a slight difference in the velocity profiles from Deborah numbers 3. Due to the better convergence of the LCR method to existing literature as demonstrated in previous section, we expect the LCR results to be the more accurate ones. The lid-driven cavity flow results were tested and confirmed for mesh convergence.

Vortex formation and growth
The flow of viscoelastic fluids in lid driven cavity flow geometries, such as the one under investigation, usually leads to the formation of corner vortices. We conclude this section by presenting the streamlines diagrams for the flow of Rolie-Poly fluids in the lid driven cavity geometry for the purposes of investigating the formation (and growth) of corner vortices.      For comparison between the LCR and DEVSS, the geometry similar to the one used for validation and has an upstream thickness of 2 H=0.025 6 m and downstream thickness of 2 h=0.006 4 m with an upstream length of 80 h and a downstream length of 50 h.
To differentiate between the geometries, we will refer to the former geometry as planar 1 and the later as planar 2.
At the inlet, a uniform velocity profile is enforced. No-slip boundary conditions are enforced at the walls. Due to the long downstream channel, we assume that the velocity profile is fully developed and hence that the inflow and outflow conditions don't affect the flow in the contraction region.
For pressure, a zero-gradient condition is imposed at the inlet and normal to the wall. At the outlet, pressure is assigned a fixed value of zero. The stress is assigned a fixed value at the inlet and a zero-gradient condition at the outlet. As in similar 2D computation studies of viscoelastic fluid dynamics, say [26][27][28] the stresses will be linearly extrapolated at the walls.

Mesh convergence
Three different hexahedral meshes were used to analyse the mesh convergence. The characteristics of the three meshes are listed in table 1 for Planar 1 and table 2 for Planar 2. All the meshes used are such that there is a higher refinement near the walls and in the vicinity of the contraction since these regions are known to possess the largest velocity and stress gradients.
Our simulated results for the velocity profile and stress profile τ xx for the full Rolie-Poly constitutive model obtained via the LCR approach and using the three meshes are compared in figure 14. The greatest differences are observed in the vicinity the contraction region. The differences are more significant for the stress profile tau xx when comparing meshes 1 and 3.
To quantitatively compare the predictions of the velocity and stress τ xx for the different meshes for geometry Planar 2, the percentage normalized error was calculated in the following way  where N is number of discretization points, X j is the value of the variable to be considered at a mesh point of the line of symmetry and X j ref are the values obtained in Mesh 3. The values of the percentage normalized error (PNE) for Deborah number De=3 calculated using equation (15) are shown in table 3.
The PNE values for Mesh 2 in relation to Mesh 3 for both the velocity and stress profiles τ xx are lower than 5% implying that greater accuracy is obtained by using Mesh 3. Due to high computational costs related to using Mesh 3, Mesh 2 will be used with an error of less than 5%.

Numerical validation
We assess the accuracy of the Rolie-Poly viscoelastic solver by benchmarking our results against the known results for the Oldroyd-B model. Similarly to the lid-driven cavity validation, we recover the Oldroyd-B constitutive model directly from our Rolie-Poly model, essentially achieved by switching off the λ R terms The numerical results for the Oldroyd-B model from the existing literature on 4:1 planar contraction flow, [29] are then used to benchmark and validate the results from our solver.
The Rolie-Poly constitutive model has a retardation constant of β=1/9 and the simulations are done at low Reynolds number flow (Re=0.01). The flow is simulated for a wide range of Deborah numbers 1De12. The evolution of the non-dimensional corner vortex size X R and the non-dimensional lip vortex size X L as functions of Deborah number, De, are also compared with results from the existing literature. Tables 4 and 5 display results from existing literature for the Oldroyd-B model [29] compared with our simulated results of the recovered Oldroyd-B model from the constitutive Rolie-Poly model, for the different mesh sizes. Table 4 displays the non-dimensional corner vortex and table 5 displays the results for the non-dimensional lip-vortex.
The results are displayed for the range of Deborah numbers from 1 to 12 for the corner vortex size and for the lip vortex size, Deborah numbers 1 to 6. The lip vortex size is a more sensitive parameter than the corner vortex size.
Generally, for the corner vortex size, good agreement is observed between our, LCR based, results and those of Pimenta and Alves [29] for Deborah numbers 1 to 5 and for the Meshes 1 and 2 and for Deborah numbers 1 up to 4 for the third mesh. The relative error for the values for Deborah numbers 1 to 4 for LCR approach is lower than 0.5 percent.
For the DEVSS method good agreement is observed for Deborah numbers 1 and 2 only with a relative percentage error of less than 3.5 percent for the first two meshes. The lip-corner vortex is a sensitive parameter and there is agreement for LCR for Deborah numbers 1 to 4 whereas for DEVSS, agreement is only observed for Deborah numbers 1 and 2.
The results indicate that our solver can be used with relative accuracy for Deborah numbers of up to 5, especially using the LCR stabilization method.

Comparison between DEVSS and LCR
Comparison of Rolie-Poly LCR and DEVSS stabilization approaches was done using Mesh 2 for Planar 2 due to the high computational cost of using Mesh 3. The axial velocity and axial normal stress as a function of the axial position at the line of symmetry for the Rolie-Poly model for different Deborah numbers were obtained. No significant difference was observed between the LCR and DEVSS methods in the axial velocities and as such only the velocity profiles for Deborah numbers 4 and 5 were presented, see figure 15. Figure 15 shows that the velocity profiles are indistinguishable, because they do not display substantial variations as with the stress profiles. The figure also shows that the flow rapidly accelerates from the fully developed flow upstream into the smaller contraction channel leading to large extensional stresses near the contraction plane as was observed by [30] in their experimental work. This causes the maximum stress τ xx to occur just upstream of the contraction plane as is seen in figure 16. The stresses then recede from this maximum.
Unlike the velocity profiles, the polymer stress profiles τ xx and τ yy , on the other hand, do display significant differences for the different Deborah numbers as illustrated in figure 16. We also notice, from figure 16, that increasing the Deborah number leads to a decrease in the maximum value of stress value τ xx as is reported in [30].
The velocity profile and stress profiles τ xx and τ yy are plotted at different positions as shown in figure 17. No significant difference is observed in the velocities, differences are however observed in the stress profiles.  It is expected that there will be slight differences between the LCR and DEVSS method since the velocity profiles have insubstantial variations as opposed to the stress profiles. Figures 18 and 19 which present the stress profiles at different x and y positions, therefore display a greater disparity between the LCR and DEVSS methods with increasing Deborah number. This could be due to the LCR method's ability to mitigate the HWNP better than the DEVSS.

Vortex formation and growth
In contraction flow geometries, such as the 4:1 contraction flow geometry, there is fundamental observable difference between the flow of Newtonian fluids on the one hand and that of viscoelastic fluids on the other; namely the formation of corner vortices in the viscoelastic flow case. Such vortices are completely absent in corresponding Newtonian flows. In this section, we therefore present the streamlines for the flow of Rolie-Poly fluids in the 4:1 contraction geometry for the purposes of investigating the formation (and growth) of corner vortices. Figures 20(a)-(g), show the formation and growth of the corner vortices for various values of the Deborah number and also for both the LCR and DEVSS stabilization methodologies. We observe a reduction in the size of the corner vortices with an increase in the Deborah number. Similar behaviour was observed by both   Alfonso et al [31] and Raphaël et al [17] in their 4:1 contraction simulations of the Oldroyd-B model. We also observe that the sizes of the corner vortices obtained via the LCR formulation are smaller than those calculated with the DEVSS method.

Discussion and conclusion
We have developed a numerical solver for the computation of viscoelastic fluid flows. The solver is based on the finite volume method and implemented in the open source OpenFOAM®software platform. Additionally the solver is developed specifically for the simulation of flows of viscoelastic fluids that are governed by the Rolie-Poly constitutive model. The Rolie-Poly model reduces to the Oldroyd-B model for certain values of the material constants. We take advantage of this relationship to validate our solver based on the existing benchmark solutions for Oldroyd-B fluids. We in particular focus attention on the two benchmark problems; lid-driven cavity flow as well as the 4:1 contraction flow. Such benchmark comparison of our data with the existing literature, [6] and [29], shows generally good agreement.
We report a significant improvement in dealing with the High Weissenberg Number Problem (HWNP) when using the LCR approach. The DEVSS stabilization method is also reliable, and also agrees with the results from the existing literature, but at lower values of the Deborah numbers than those achievable with the LCR approach. Based on our results, we can therefore confirm than the LCR stabilization technique is more reliable than widely used DEVSS stabilization method. We therefore recommend the LCR technique for the numerical stabilization of viscoelastic flow computations going forward.
In conclusion, this work has presented a new solver for the simulation of flows of viscoelastic fluids governed by the Rolie-Poly constitutive model. The results show that the solver is consistent with expectations as can be observed from the favourable comparison with results from the existing literature. Our results also demonstrate that the LCR method better mitigates against the HWNP as compared to the widely employed DEVSS method. Further studies however need to be conducted to investigate the disparities in the results for Deborah numbers greater than 7 especially for the LCR representation.