Viscoelastic fluid flow simulation using the contravariant deformation formulation

The new formulation for conformation-tensor based viscoelastic fluid models, written in terms of the contravariant deformation tensor as introduced by Hütter et al. (2018), has been approached numerically. It has been implemented in a finite element method framework with the DEVSS-G/SUPG method for stabilisation. A stress-implicit as well as a stress-explicit formulation is presented that allows to integrate the non-linear flow equations, with or without a solvent contribution, forward in time with second-order accuracy. The time-dependent stability in shear flow is maintained using the contravariant deformation, which is tested by solving a perturbed planar Couette flow of an upper-convected Maxwell (UCM) fluid. The stability and accuracy of our implementation has furthermore been tested by solving the flow around a cylinder with confining walls. The new formulation turns out to be more stable as compared to the formulation in terms of the conformation tensor. It enables to simulate the flow of a Giesekus fluid with non-linear parameter α = 0.01 and viscosity ratio ν = ηs/η0 = 0.59 way beyond the Weissenberg number at which the High Weissenberg Number Problem manifests itself using the standard formulation. The stability appears to be similar to the log-conformation representation. The computed stress profiles up to a Weissenberg number of Wi = 0.6 for an Oldroyd-B ∗Corresponding author Email addresses: m.a.carrozza@tue.nl (Mick A. Carrozza), m.a.hulsen@tue.nl (Martien A. Hulsen), m.huetter@tue.nl (Markus Hütter), p.d.anderson@tue.nl (Patrick D. Anderson) Preprint submitted to Journal of Non-Newtonian Fluid Mechanics June 29, 2019 fluid with viscosity ratio ν = 0.59, compare well with benchmark results from other studies, including other finite element, finite volume and spectral element methods.


Introduction
Several numerical stabilisation techniques have already been developed for mixed finite element methods using conformation-tensor based differential models for simulating the flow of viscoelastic fluids [1]. For stabilising the velocitystress interpolation, the DEVSS-G technique [2,3] with smooth velocity gradient 5 interpolation in the constitutive equation [4,5] for additional time-dependent stability in shear flow, and variants thereof are often employed. In addition, stabilisation of the convection in the constitutive equation is obtained with SUPG [6]. The decoupling of the momentum balance equation from the constitutive equation using the implicit stress formulation [7] allows for second-order time- 10 integration of these equations avoiding a fully implicit coupled scheme, even for small or zero solvent viscosities. Still, with all these techniques available, until 2004 [8,9], all methods for numerical simulation of the flow of viscoelastic fluids broke down when the Weissenberg number exceeds a relatively small critical value. This value of O(1) depends on the problem, the constitutive 15 model, the numerical method and the mesh. The High Weissenberg Number Problem (HWNP) is conjectured to be a numerical instability as a result of the failure of the numerical schemes -that make use of polynomial basis functions for approximating the stress -to balance exponential growth of the stress due to deformation with convection [9]. 20 In Fattal et al. [9], the HWNP is restrained by rewriting the equation for the conformation in terms of the logarithm of the conformation, leading to crucial stabilisation of the numerical schemes. This way, the exponential profiles become linear in the logarithm of the conformation. In addition, the positivedefiniteness of the conformation is preserved. Hulsen et al. [10] were the first 25 to implement the log-conformation representation in a finite element method (FEM) context using the DEVSS/DG formulation. However, the required spectral decomposition in this formulation makes numerical implementation difficult, especially when it comes to linearisation required for Newton-Raphson iteration [11,12] to solve the governing equations in a fully coupled manner. Another 30 approach that alleviates the HWNP can be found in [13,14], which use the evolution of the unique symmetric positive-definite square root of the conformation. Finally, in [15] a Cholesky decomposition of the conformation tensor is evolved which keeps positivity of the conformation. All the techniques mentioned in this paragraph enable to simulate and/or obtain improved convergence 35 beyond Weissenberg numbers at which the standard conformation formulation fails. Note, that similar decompositions of the conformation tensor are used in [16,17] for turbulent flow analysis purposes but not for flow computations.
Hütter et al. [18] recently developed a general approach for including thermal fluctuations in conformation-tensor based (non-)linear viscoelastic models 40 in accordance with thermodynamic principles. They generalised these models and reformulated them in terms of a multiplicative decomposition of the conformation tensor which they call the contravariant deformation tensor, and can be interpreted as the elastic part of the deformation gradient tensor. One of the numerical advantages of the reformulation by Hütter et al. [18] is that it 45 is also relevant for deterministic models, i.e. without fluctuations, and we anticipate that the HWNP is naturally avoided. The reason for the avoidance is the resemblance of the evolution equation for the contravariant deformation to the formulation of the constitutive model in the Brownian configurations fields method [19], which is not affected by the instability, regarding deformation and 50 convection [20]. In contrast to approaches [13][14][15], the contravariant deformation tensor (CDT) formulation contains deformation rotation, whereas the other approaches eliminate this rotation of the conformation-related variable to obtain a symmetric or triangular structure. The inherent rotation is a drawback of the CDT formulation, making it difficult to obtain a steady state solution directly 55 as, for example, needed in a linear stability analysis. However, the approaches [13][14][15] contain additional complicated terms in the evolution equation to comply with the restriction to symmetry or triangularity, respectively, making them computationally more expensive and less amenable to the inclusion of thermal fluctuations, in constrast to our CDT formulation.

60
The goal of this paper is to show that the stability and accuracy of the formulation of differential constitutive models for viscoelastic fluids without fluctuations in terms of the contravariant deformation [18] is similar to the logconformation formulation. We aim at providing stable numerical methods to solve time-dependent viscoelastic flows, which at a later stage can be extended 65 with fluctuations according to [18].
The paper is organised as follows. The governing equations including the new formulation of viscoelastic models are described in Section 2. The numerical method is explained in Section 3, where two different second-order timeintegration schemes for viscoelastic flows are presented: a stress-implicit as well 70 as a stress-explicit formulation. Section 4 presents results for homogeneous flows using the new formulation. Results for the planar Couette problem with an initial perturbation are shown in Section 5. In Section 6 the results for the flow around a cylinder with walls at low and high Weissenberg numbers will be shown, comparing the new formulation with the (log-)conformation represen-75 tation and with results from the literature. Conclusions are drawn in Section 7.

General aspects
The flow of incompressible viscoelastic fluids for which inertia is neglected is governed by the momentum and mass balance equations, with pressure p and velocity vector u, η s is the solvent viscosity, D = (∇u + (∇u) T )/2 is the rate-of-deformation tensor with velocity gradient tensor (∇u) ij = ∂u j /∂x i , and τ is the viscoelastic extra stress tensor. Furthermore, a constitutive equation for τ is required, relating it to the deformation history of the fluid. Following Hütter et al. [18], in this paper conformation-tensor based viscoelastic models will be used that are formulated in terms of the contravariant with() the material derivative, L = (∇u) T is the transposed velocity gradient tensor, and g(b) is a tensor function of b whose form depends on the viscoelastic model. Tensor b is related to the conformation tensor c via Note that b is generally an unsymmetric tensor. Moreover, it can be seen that b is not unique since an arbitrary rotation of b according to b = b · R, with R a rotation tensor, does not change the conformation. Hence, no steady state for b is to be expected in steady-state flows. As will be shown later, the tensor b is unsteady in shear flow, even if c is steady. A general form of the evolution with f (c) an isotropic, model specific tensor function of c. Substituting Eqs.
(3)-(4) into Eq. (5) by applying the product rule for differentiation, it can be shown that It has been argued by Hütter et al. [18] that for numerical reasons it is 80 preferable to write viscoelastic models extended with fluctuations in terms of b rather than in c. One argument is that, as a consequence of Eq. (4), c is always (semi-)positive definite. Moreover, the formulation in terms of b is expected to be numerically more stable than the formulation in terms of c, even without fluctuations. That is, the evolution equation of b shows close similarities with 85 the stochastic evolution equation in the method of Brownian configuration fields [19], which is immune to the numerical instability that causes the high Weissenberg number problem (HWNP) [20]. The stability comes from the intrinsic balance between deformation and convection in the evolution equation for b.

90
In this paper, two viscoelastic fluid models will be considered: the upper- the considered models, g(b) and τ are given by [18]:

95
• UCM/Oldroyd-B model: with λ the relaxation time. The extra stress expression for this model is with G the shear modulus.
• Giesekus model: where the dimensionless parameter α determines the magnitude of the anisotropic drag. The extra stress expression for this model is the same as for the UCM/Oldroyd-B model: The Giesekus model (Eq. (9)) reduces to the UCM/Oldroyd-B model (Eq. (7)) for α = 0.

Polar decomposition of the contravariant deformation
The contravariant deformation can be written as follows (see Eq. (4)): where V = √ c is a symmetric, positive definite tensor and R is a rotation tensor.
It should be noted that which is a skew-symmetric tensor. If the material derivative of V vanishes (V = 0), for example in the case of stationary homogeneous flow, Eq. (12) reduces toṘ Setting the right-hand side of Eq. (13) to zero and then left-and right-multiplying where ω is the angular frequency. Note that ω/γ → 1/2 for Wi → 0 (γ = 0), which is the maximum value of ω/γ, and ω/γ → 1/(2Wi) for Wi → ∞.

Numerical disretisation methods
The governing equations of Sec. 2 are discretised in space using the finite element method and the method of lines for the time discretisation applying  [21,22], as will be shown later.
In order to stabilise the convection term in the constitutive equation (Eq. (3)), 125 the SUPG method [6] is applied. This method uses a modified test function d = d + κu · ∇d, with test function d for the standard Galerkin method and κ the SUPG-parameter, which puts more weight to the upwind direction of the streamline.

Weak formulation
where in Eq. (22) the material derivative is split into a local time-derivative term and a convective term, since an Eulerian grid is employed. Furthermore, (., .) indicates the L 2 -inner product on the domain Ω, which is divided into elements of type quadratic triangle. The SUPG-

Time integration
For the time discretisation of the equations, the considered time interval is 145 divided into a number of time steps. The time step ∆t is taken constant, the current time being equal to t n = n∆t, where n is the number of steps already performed. In the following, the evaluation of field variables at a discrete time will be denoted by a subscript; τ n = τ (t = t n ) for example for the viscoelastic extra stress. Two formulations for the time discretisation of the equations will 150 be used: the explicit and implicit stress formulations for viscoelastic fluid flow.

Implicit stress formulation
For flows with a small or zero solvent contribution, the implicit stress formulation [7] is often used to decouple the momentum balance, mass balance and projection equation from the constitutive equation. In this formulation, an expression for τ n+1 is found which still contains unknown velocity terms u n+1 so that the system in Eqs. (19)- (21) is not singular. This results in the set of momentum balance, mass balance and projection equations which is a Stokes-like system in the unknowns u n+1 , p n+1 and G n+1 at the new time t n+1 . An L 2 -projection on the discretised space defined in Eq. (18) is used to obtain c n from b n (see Eq. (4)): with P ∈ B h a test function. After obtaining the solution of these equations, the constitutive equation (Eq. (22)) is solved using a BDF2 time-discretisation scheme with contravariant deformation b prediction: with predictionb n+1 = 2b n − b n−1 . This two-step time-stepping scheme is started by using a first-order semi-implicit scheme in the first time step:

Rotation reinitialisation of b
Very refined meshes to resolve the small length scales due to differences in the rotation of b between different parts of the flow, see Sec. 2.3, can be omitted if R in Eq. (11) is reset to the unit tensor I by computing V , which leaves the conformation unchanged, see Eq. (4). To keep the second-order accuracy of the time-discretisation schemes presented in this section, b at the discrete times t n and t n−1 needs to be reset by the same rotation tensor for reinitialisation.
Therefore, V n is computed in the nodes for resetting the rotation of b n , whereas the rotation of the nodal values of b n−1 is reset according to: where Eq. (11) has been used.
The reinitialisation is performed whenever a measure for the skew-symmetry of b in practice exceeds a critical value anywhere in the entire domain. For this we use the quantity: which is the squared magnitude of the skew-symmetric part of the contravariant  the results of these computations, we use the following notation for the first and second normal-stress difference N 1 and N 2 , and definition of the Weissenberg number Wi:

Convergence in time
For simple shear flow, the exact solutions for the shear stress and first normal-stress difference of the UCM model are given by The relative error in τ xy and N 1 for simple shear at Wi = 1 obtained from the b-form is plotted in Fig. 2(a)       the exact solution, a quantification of the approximation error is calculated For uniaxial extension, the exact solution for the first normal-stress difference of the UCM model reads The relative error in N 1 of the approximations obtained with the b-form for 205 Wi = 0.1 and different time steps ∆t is plotted in Fig. 3(a) as a function of time. It can be seen that after start-up of the flow, the error keeps decreasing with time, because a steady-state b will eventually be obtained for extensional flow, see Eq. (13). In Fig. 3(b) where the quantification of the approximation error is plotted, it can be seen that the slope of −2 corresponding to the accuracy

Time-dependent stability in planar Couette flow
A planar Couette flow of a UCM fluid is considered, as illustrated in Fig. 4, which is known to be a stable flow [29].  re-enters via Γ p,l and that periodic copies of the box in horizontal direction can be imagined. Therefore, periodic boundary conditions for the velocity and contravariant deformation are used: b| Γ p,l = b| Γp,r . extended in the x-direction, are: where t = σ·n is the traction with n the outwardly-directed unit-normal vector.
The flow is assumed to be symmetric and therefore only half of the domain Ω is considered, whereas symmetry boundary conditions for the velocity and traction on the centreline Γ cen are used: where t x is the x-component of the traction. Furthermore, no-slip boundary conditions on the channel walls Γ w and cylinder wall Γ cyl are assumed: Although b is allowed to be an arbitrary rotation tensor at time t = 0 for an initially stress-free fluid (see Eq. (4) and Eqs. (8) and (10) Table 1) and the time step is ∆t = 0.01λ.  Table 1) and a time step of ∆t = 0.01λ are used for which convergence is obtained (see Sec. 6.2.1). Fig. (b) gives a zoom-in on the data.
viscoelastic extra stress tensor in the domain is plotted as a function of time at Wi = 0.5 for different values of a. It can be seen that there is a relatively large deviation in τ xx,max for a = 0.975 from τ xx,max for a = 0 (1% deviation max.) in comparison with τ xx,max for a = 0.97 (0.3% deviation max.). Compared to no reinitialisation at all, the improvement is even more dramatic. This is shown in Fig. 9, where the viscoelastic extra stress components are plotted for the case  Table 1 gives the numerical parameters of the meshes that are used, which are the same meshes as in [21].
In Fig. 10 It can be seen that τ xx converges for meshes M3 and M4 with both the b-( Fig. 10(a)) and log c-representation ( Fig. 10(b)). This holds for the other stress components and field variables as well.  (see Table 1) and a time step of ∆t = 0.01λ are used for which the solution converges (see Sec. 6.2.1).
For checking the time convergence of the computations with the b-representation, the following quantification of the error in viscoelastic extra stress on the cen- Table 1: Characteristics of the meshes generated by Gmsh [30], used for the flow around a cylinder confined between walls. The meshes, which are the same as in [21], are more refined near the cylinder wall and element sizes h are halved at each refinement. N elem and N elem,cyl mean the total number of elements and the number of elements on the cylinder wall, treline and cylinder wall with respect to a reference approximation is used, where τ i and τ ref,i are the nodal values of a stress component for the approxi-290 mation and the reference solution, respectively, and N is the number of nodal points on the centreline and cylinder wall. The convergence is quadratic as is shown in Fig. 11(a), where the steady-state τ −τ ref 2 is plotted as a function of ∆t, with τ ref the approximation for ∆t = 0.0005λ. This is a result of the secondorder time discretisation of the constitutive equation. Fig. 11(b) shows the same 295 time-convergence plot as in Fig. 11(a), but here for the explicit stress formulation. It can be seen that the rate of convergence is the same but the error is somewhat smaller than with the implicit stress formulation, shown in Fig. 11(a).
A quantification of the difference in stress of the b-formulation with respect to the log c-formulation is plotted in Fig. 12 as a function of ∆t. As can be seen,   becomes of the same order of magnitude as the space-discretisation error. The plateau is due to the numerical difference between the band log c-formulation.

305
From the results in this section it can be concluded that the solutions for the explicit and for the implicit stress formulation using the bas well as the    Wi the stress wake downstream of the cylinder extends along the centreline with respect to that obtained at the low Wi = 0.5, the mesh is refined near the entire centreline apart from the refinement near the cylinder wall. The characteristics of the meshes, which are comparable to the meshes in Table 1 and that are also used in [21], are given in Table 2. The front and rear stagnation points are at s = 0 and s = πR, respectively.
Convergence on the cylinder and in the wake are clearly visible. The same Table 2: Characteristics of the meshes generated by Gmsh [30], used for the flow around a cylinder confined between walls at moderate and high Weissenberg number. The meshes, which are the same as in [21], are refined near the cylinder wall and the entire centreline.  Table 2) and the time step is ∆t = 0.01R/U . stress profile obtained with mesh M3r is compared to results from different 325 studies in Fig. 14(a) and good graphical agreement is found. This becomes even  Fig. (b) gives a zoom-in on the profile in the wake. Mesh M3r is used (see Table 2) and the time step is ∆t = 0.01R/U . more apparent in Fig. 14(b), where a zoom-in on the profile of Fig. 14(a) in the wake is shown.
For Wi = 0.6 a similar convergence of τ xx on the cylinder as for Wi = 0.5 is observed in Fig. 15. In the wake, convergence is visible as well, yet τ xx is  Table 2) and the time step is ∆t = 0.01R/U .
profile computed with mesh M3r for Wi = 0.6 with results from the literature is made in Fig. 16(a) and in Fig. 16(b), which gives a zoom-in on the profile of Fig. 16(a) in the wake. It can be seen that the results overlap quite well, 335 even though [35] predicts a lower maximum on the cylinder. The differences are conjectured to be due to mesh dependency of the solutions.
A lack of convergence in the wake emerges when inspecting τ xx on the cylinder wall and centreline for Wi = 0.7, which is revealed in Fig. 17(a), despite convergence on the cylinder. Nevertheless, a comparison is made between our 340 computation of the stress profile with the most refined mesh M3r and in some cases mesh-dependent results from the literature for the most refined mesh, presented in Fig. 17(b) . Good agreement between the calculated values is obtained on the cylinder while a relatively wide variety of values is noticed in the wake. Knechtges et al. [ (b) Figure 16: Comparison of τxx as a function of the path length along the cylinder and centreline with results from the literature for Wi = 0.6. The front and rear stagnation points are at s = 0 and s = πR, respectively. Fig. (b) gives a zoom-in on the profile in the wake. Mesh M3r is used (see Table 2) and the time step is ∆t = 0.01R/U . (a) Component τxx for different meshes (see Table 2). Knechtges et al. [11] Claus et al. [33] Coronado et al. [34] Afonso et al. [35] Damanik et al. [36] Fan et al. [37] Owens et al. [38] (b) Comparison of τxx with results from the literature.  Table 1 and that are also used in [21], are given in Table 2.
The solution will be called stable for a certain Wi if it does not diverge, The stability results are shown in Table 3 and are the same for the explicit and the implicit stress formulations. It seems that there is no limit for stability

Conclusions
The new formulation for conformation-tensor based constitutive models for ample Newton-Raphson or other iteration methods, is straightforward and far less cumbersome than for the log c-representation [11,12], the latter involving a spectral decomposition. Nevertheless, since b remains unsteady in general flows, it is unclear how to search for a steady-state solution with Newton-Raphson iteration. Apart from that, application of the b-formulation to linear stabil-400 ity analysis, which generally requires a linearisation of the equations around a steady-state solution, might not be straightforward; this requires further study.
There are a few remarks to make concerning the application of the b-representation in numerical simulations. Firstly, it entails additional time-discretisation errors in a steady-state solution compared to the log c-formulation, 405 since b keeps rotating in shear. Reinitialisation of the rotation of b in a field is needed to prevent the occurrence of large gradients in b due to the rotation, retaining accuracy of the solution at negligible extra computational cost. Note that reinitialisation has no effect on the time-discretisation error, however, and for this reason it is not at all required for homogeneous flows. Any second-order tensor can be written as the sum of a symmetric part and a skew-symmetric part. Applying this to the contravariant deformation b gives: with b S the symmetric part of b and b A the skew-symmetric part. Substituting this, the quantity in Eq. (35) can be rewritten as: , (B.2) making use of tr(b S · b T A ) = tr(b A · b S ) = −tr(b A · b S ) = 0 since b T A = −b A . It can now be seen that 0 ≤ b A ≤ 1 for b = 0 as tr(b 2 S ) ≥ 0 and tr(b A · b T A ) ≥ 0.