The deformation fields method revisited : stable simulation of instationary viscoelastic fluid flow using integral models

The implementation of the deformation fields method for integral models within a finite element context [1, 2] has been updated with various techniques to have a numerical stability that is comparable to state-of-the-art implementations of differential models. In particular, the time-dependent stability in shear flow, decoupled schemes for zero or small solvent viscosities and the log-conformation representation now have counterparts in the numerical implementation of integral models leading to similar numerical stability. The new techniques have been tested in transient shear flow and the flow around a cilinder confined between two plates for the integral version of upper-convected Maxwell model and for integral models having a non-constant damping function.


Introduction
Integral models are frequently used for modelling the constitutive behaviour of polymer liquids. In particular the factorable K-BKZ equation, where the linear rheology is specified by the memory function and the non-linear rheology is given by the damping function, is very popular. For an extensive review of the various damping functions that have been used see [3]. Also, new constitutive models based on molecular considerations are often formulated in terms of memory integrals, see for example [4]. Furthermore, recent developments in modelling complex fluids with a broad power-law relaxation spectrum are also of the integral type [5].
Various techniques have been proposed to simulate flows of viscoelastic fluids described by integral models. For a review of the various finite element method implementations for integral models see [6]. As has been concluded in that review, currently only two methods are capable of performing transient simulations: the Lagrangian Integral Method (LIM) [7,8]

and the Deformation Fields
Method (DFM) [1,2]. LIM is suitable in a Lagrangian formulation, whereas DFM can be used in any formulation, but is used mainly in an Eulerian (fixed domain) or Arbitrary Lagrange Euler (ALE) formulation.
Despite the many good constitutive models available, numerical simulation of viscoelastic fluid flows using integral models is less popular than simulations with differential models. The reason is that implementation of integral models is more involved because of the complete deformation history that needs to be determined. As a result, new stabilization techniques have been specifically developed for differential models and the numerical techniques for integral models are lagging behind in this regard. In particular, the time-dependent stability in shear flow of the velocity-gradient projection technique [9,10], decoupled schemes for zero or small solvent viscosities [11] and the log-conformation representation (LCR) [12] do not have a similar counterpart in the numerical implementation of integral models.
In this paper we propose and implement various techniques within a finite element context to give the DFM implementation of integral models stability properties comparable with state-of-the-art finite element techniques for differential models. To show stability and accuracy, we apply the techniques to a perturbed single shear flow and the standard benchmark geometry of the flow around a cylinder confined between plates.
The paper is structured as follows. In Section 2 we present the balance equations and the constitutive models used. Also the DFM is described and rewritten in terms of the deformation gradient instead of the Finger deformation tensor. In Section 3 the numerical discretization based on the finite element method is described. In Section 4 simulation results are presented. Finally, in Section 5 we end with a discussion and final conclusions.

Balance equations and constitutive models
Consider a flow in a fixed domain Ω (see Fig. 1). We assume that inertia in ⌦ \ in u n Figure 1: Flow domain Ω with boundary Γ where on part of the boundary (Γ in ) the velocity u is directed inwards (u · n < 0).
can be neglected and that the fluid is incompressible. Therefore the balance of momentum and mass become where u is the velocity vector field, p is the pressure field, τ is the extra-stress tensor field, η s is the solvent viscosity of the fluid and D = (L + L T )/2 is the symmetric part of the velocity gradient tensor L = (∇u) T . All fields are dependent on time t and position in space x, i.e. u(t, x), but function arguments will only be added when needed. The extra-stress tensor field τ is given by a memory integral of the K-BKZ type: where B = B(t, t , x) is the Finger deformation tensor of a material particle at the current position x and time t, with respect to the reference configuration is the memory function and h(I 1 ) is the damping function depending on the first invariant I 1 = tr B.
For the purpose of this paper, but without loss of generality, we restrict m(τ ) to be a single exponential where G is the modulus and λ is the relaxation time. For the damping function we will use as introduced in [13] and also used in [5] for the fractional K-BKZ integral model.
Note, that for a = 0, η s = 0 the constitutive model reduces to the integral form of the upper-convected Maxwell (UCM) model.
For verification of the results, we use the UCM model in differential form, given by where c is the conformation tensor. For the UCM model η s = 0. Note, that in equilibrium c = I and τ = GI. For integral models it is quite common to drop isotropic terms from the extra-stress expression (see, for example [14]).
To easily compare the integral and differential formulation of UCM, we did the same for the latter here also.

Initial and boundary conditions
It is assumed that before the start time t = 0 the flow is at rest, which means that B(t, t , x) = I for t ≤ t < 0 and all x ∈ Ω. The stress tensor in this equilibrium state is τ = GI. The boundary conditions on Γ are similar to the (Navier-) Stokes equations, for example specification of the velocity vector on one part of Γ (Dirichlet) and the traction vector on the remaining part of Γ (Neumann). If the velocity vector on the boundary is such that there is an inflow boundary Γ in (see Fig. 1), the fluid deformation history at that boundary needs to be specified. However, if the solvent viscosity is zero (η s = 0) it is not allowed to specify the fluid history completely, similar to the restrictions on imposing stress components at the inflow boundary for differential models [15,16]. In order to avoid complications related to inflow boundaries, we assume periodical boundary conditions for the actual problems considered in this paper. In case of periodical boundary conditions the inflow boundary is coupled to outflow boundary, therefore the flow in the domain Ω provides its own deformation history for the fluid particles entering the domain.

Deformation fields
For the evaluation of the memory integral in Eq. (3), we need the deformation history B(t, t , x), −∞ < t ≤ t for the particle that is currently at the position x. In a Lagrangian description, having the positions of the material particles as a primary unknown, it is relatively easy to obtain the deformation gradient F (t, t , x) = ∂x/∂x , −∞ < t ≤ t and then compute the Finger tensor via B = F · F T . This is the approach taken in [7,8] and works well for problems where a fixed amount of fluid is being studied, such as stretching a filament.
However, for problems in a fixed domain the Eulerian description is preferred.
For that purpose, the method of deformation fields [1,2] has been developed.
Assume that at t = 0 the flow starts and there has been no deformation of the material before that. One way of obtaining F (t, t , x) for t ≤ t is defining a (large) number of discrete reference times t i (both for t < 0 and t > 0) and solve the following equatioṅ forward in time for all discrete reference times t i ≤ t, i.e. a discrete time point 'kicks in' as soon as time t passes the time point. In an Eulerian frame, this equation is written as follows, The solutions for F (t, t i , x) are called deformation fields. There is a problem with this approach: the number of 'active fields', i.e. the fields with t i ≤ t keeps increasing as time increases. This can be overcome by removing the 'oldest field' as a new field becomes active, hence keeping the number of active fields constant. This is basically the first generation of the deformation fields method as introduced in [1], except that in [1] B is computed directly (see below) and not with F as an intermediate field. A drawback of this approach, in contrast to differential models, is that the fields are unsteady for steady flows, which leads to additional time discretization errors. Also the labelling of the fields by the reference time leads to difficult bookkeeping, including the removal and creation of fields at each time step. As a result, also the distribution of fields over the history is not very flexible and basically coupled with the time step ∆t.
The second-generation deformation fields [2] removes much of the drawbacks of the first generation by labeling the deformation gradient tensor not by an absolute (fixed) reference time t , but by the reference time relative to the current time, i.e. by the age τ = t − t : F (t, τ, x). For that, Eq. (7) is rewritten in terms of the new set of independent variables as follows: which is a 4D hyperbolic equation (assuming 3D space) with a 'convection speed' of 1 in the τ -axis direction. The boundary condition (for the τ -axis) is and the initial condition (fluid at rest for t < 0) is In principal the age axis is unlimited, i.e. τ ∈ [0, ∞). However, in practice we are dealing with 'fading memory fluids' and the age axis is taken to be finite i.e.
Note, that in [2] the Finger tensor B(t, τ, x), τ ∈ [0, τ c ] is computed directly with boundary condition (for the τ -axis) and initial condition (fluid at rest for t < 0) We call this the B -approach, which has the advantage that the number of components to solve is less, since B is symmetric. A disadvantage is that the equation to solve contains an additional term as compared to the F -approach given by Eq. (9), possibly increasing the computational effort. More importantly however, it turns out that the B -approach is less stable than the F -approach and that the latter is the preferred approach, as we will show later.

Numerical discretization
The unknown field variables in the F -approach are the velocity u(t, x), the pressure p(t, x) and the deformation gradient F (t, τ, x). The governing equations are Eqs. (1) to (3) and Eq. (9). We will use the finite element method to discretize in real space (x) and in 'age space' (τ ). Time discretization will be based on second-order backward differencing (BDF2). We start with the discretization of the deformation gradient.

Weak formulation and discretization of the deformation gradient
We will apply the finite element method to discretize Eq. (9) in both the τdirection and in real space. For that, we write the weak formulation as follows: for all test functions H(t, τ, x). Here, (·, ·) is the usual inner product for tensors in real space. Note, that we have limited integration in the τ -direction to the finite interval with a cut-off time, i.e. [0, τ c ].
Although quite general four-dimensional elements could be used for discretization, we will not do that. First, we divide the integration domain [0, This division is global and independent of t and x, thus creating '∆τ ' slabs. We will employ 'stretched' grids, as in [2], where ∆τ k = f ∆τ k−1 (k = 0) with constant factor f . The factor f is uniquely determined by N , ∆τ 0 and τ c and can be found by iteration. In this work, we have used ∆τ 0 = 3λ/N , which is close to the optimal value for a single exponential memory function, as determined in [2]. When appropriate, the value of f used in the simulations is included in the caption of the figures.
We apply the discontinuous Galerkin method in τ -space, using piecewise continuous base functions, that can jump across the slab boundaries. In real coordinate space we employ a standard mesh with continuous base functions.
Since we solve a hyperbolic equation, we apply SUPG stabilization. We can now write the discretized weak form on a single slab: where is the jump of F at the slab boundary τ k and τ is the SUPG time scale parameter. Here, the superscripts + and − denote the value at τ k in interval k and k − 1, respectively. Note, that the use of the DG method together with the division into intervals (slabs) makes it possible to solve F in the τ -direction 'slab by slab' like a time step scheme. First the entry condition F − (τ 0 ) = I is imposed, then slab 0 is solved, which gives F − (τ 1 ) as input for slab 1, etc.
Now we introduce the base functions for F (t, τ, x) in τ -direction by writing the test and trial functions as follows and the indices i and j are limited to fields in slab k.
The system Eq. (18) can be solved 'slab for slab' following the sequence k = 0, . . . , N − 1, where the solution values at the end of the slab, i.e. F − (t, τ k+1 , x), are input for the next slab. This is similar to a time-stepping scheme. In fact, the DG-method using order q polynomials leads to a 2q + 1 order accurate implicit Runge-Kutta method [17].
Time discretization is based on backward differencing, where we combine implicit BDF2 for all terms, except −L · F j , which is treated with extrapolated BDF [18] of the same order. In this way, the components of the fields are decoupled and the system matrix for each component is the same. The linear system is non-symmetric and a direct solver from the HSL-library (MA41) [19] is used in this work to obtain the solution for F j .
After computing the deformation fields, memory integrals, like in Eq. (3), are numerically evaluated by applying a standard Gauss rule on each slab in τ -direction and adding the results of all N slabs.
The B -approach is similar to the F -approach, applying the same discretization methods and base functions to Eq. (12) for the tensor field B, instead of Eq. (9) for the tensor field F . Note, that in the F -approach the B-tensor field is computed by F · F T and is therefore guaranteed to be positive semidefinite everywhere, whereas in the B -approach positive semidefiniteness of B might be locally violated due to poor resolution of the mesh. The latter is also true for the stress tensor τ in the B -approach, although to a lesser degree, since it is a weighted sum of many B-fields.

Discretization of the momentum and mass balance
Our aim is to formulate a discretization of the momentum and mass balance that is able to simulate time-dependent flows having zero or low-solvent viscosity contribution. For that, we proceed similarly as in [11], where the differential In order to apply a similar procedure for integral models, we need to obtain a 'differential form' of the integral model Eq. (3). For that, we differentiate the integral model with respect to time. This is performed in Appendix A and the result is where Note, that the second-order tensors τ , A and the fourth-order tensor K can be computed once the deformation history B(τ ) is known.
The discretization of the momentum and mass balance is now performed as follows. At the next discrete time (t n+1 ) the momentum and mass balance are evaluated: An expression for τ n+1 is found by integrating the differential expression Eq. (21) from t n to t n+1 = t n +∆t using a combined explicit-implicit Euler scheme, where τ , A and K are evaluated at t n and u at t n+1 : After evaluation of τ n+1 and substitution into Eq. (25) we get the following 'Stokes' like system for u n+1 and p n+1 : Note, that the error in τ n+1 is O(∆t 2 ). Therefore, we can expect an error of also O(∆t 2 ) in the time discretization of the momentum balance. Since the deformation at time t n , i.e. B n , is known and a continuous field, the system is well-defined and can be discretized in space using standard velocitypressure finite-elements. For that, we use Taylor-Hood elements (P 2 /P 1 on triangles/tetrahedra, Q 2 /Q 1 on quadrilaterals/hexahedra).
One practical problem in the finite element discretization is the evaluation of the stress convection term u n+1 · ∇τ n . We avoid taking spatial derivatives of τ n , which is not in a P 1 (triangles/tetrahedra) or Q 1 (quadrilaterals/hexahedra) space, by first projecting τ n on a P 1 or Q 1 space and using the projected stress field in evaluating the stress gradients.
We note that in this work we use the DEVSS-G method [10,20] for additional stabilization of the momentum balance. This involves adding an additional term −α∇·(∇u n+1 −G T n+1 ) to the left-hand side of Eq. (28) and solving an additional equation The auxiliary variable G is interpolated using continuous basis functions (P 1 on triangles/tetrahedra, Q 1 on quadrilaterals/hexahedra), which means that G is the least-squares projection of the velocity gradient (∇u) T on the P 1 or Q 1 continuous finite element space. For the coefficient α we take the zero-shear-rate polymer viscosity, which is Gλ for the memory function given in Eq. (4).
The linear system is non-symmetric and a direct solver from the HSL-library (MA41) [19] is used in this work to obtain the solution for u n+1 , p n+1 and G n+1 .
The solution for u n+1 is used to obtain the deformation fields at the new time t n+1 , i.e. F j (t n+1 , x), using the discretization discussed in Sec. 3.1 with L = (∇u n+1 ) T . Similar to differential models [9,21], it is also possible to use L = G n+1 instead. It will be shown below that the latter option leads to much better temporal stability in shear flow.

Planar Couette flow of an upper-convected Maxwell fluid
We where rand is 2 × 2 tensor where each of the components in the (x, y) Cartesian system is a random number on the interval [0, 1], set independently for all nodes and all fields. Note, thatγτ is the shear strain, which is equal to Wi over a time span equal to λ. Therefore, the perturbation is assumed to be small enough for the difference from the steady state to be considered a linear perturbation for high values of Wi. Results for the integral UCM at Wi = 10 for a 20 × 20 mesh (with quadrilateral elements) are shown in Fig. 3. The numerical technique applied is described in [11] (without LCR). The numerical parameters  obtained using a differential UCM model (see Eq. (6)) in Fig. 4. We see that our implementation of the deformation fields behaves quite similar to the differential model implementation. We get the typical stable curve ( [9,21]) when L = G is used in the evolution of F . After an initial increase, the perturbation will die out and eventually becomes a straight line in the log plot. The form of the transient for the continuous spectrum, as predicted by Kupferman [22], is also shown. Eventually, there is an even slower decaying eigenvalue, which is probably related to the two discrete Gorodtsov & Leonov (GL) modes [23]. The discrete GL-modes decays as exp(−t/(2λ) for Wi 1 and are located near the walls. Whether we resolve the GL-modes with this rather coarse mesh, is beyond the scope of the current paper.
Finally, we note that using L = (∇u) T , i.e. obtaining the velocity gradient directly from differentiating the finite element velocity field, leads to a numerically unstable flow for the integral UCM model using a deformation fields method, as shown in Fig. 3. This instability, is similar to the well-known behaviour when using the differential UCM with L = (∇u) T ( [9,21]), as shown in Fig. 4.
For the remaining part of this paper, we will use L = G throughout.

Problem description
We consider the planar flow past a cylinder of radius R, centrally positioned between two flat plates (see Fig. 5). The distance between the two plates is 2H, with H = 2R and the length of the flow domain is 30R. In the following we will use an (x, y) coordinate system with the origin at the center of the cylinder. We assume the flow domain to be periodically extended in x-direction with period

Mesh
We use two series of uniformly refined meshes, created using Gmsh [24].
The element sizes are imposed by specifying the size in geometrical points on the boundary of the model in Gmsh. The characteristics of the first series of pbc pbc L = 30R x R y H = 2R 1 Figure 5: Geometry of the flow past a cylinder confined between two walls. Periodical boundary conditions (pbc) for velocity and stress (differential model) or deformation fields are applied between the inflow and outflow boundary.
meshes used (M0-M4) are given in Table 1. In Fig. 6 the starting mesh M0 is shown. The other meshes are refinements (starting from M0) using a factor of two in each coordinate direction, as compared to the previous refinement. The characteristics of the second series of meshes (M0r-M4r) are given in Table 2 and are similar to the corresponding meshes in Table 1, except that the level of refinement near the cylinder is extended to the complete center line. In Fig. 7 the starting mesh M0r is shown.    First, we solve the problem for a differential UCM model at Wi = 0.5. The numerical technique used is described in [11] (without LCR). After starting the flow we simulate up to t = 20λ with a time step ∆t = 0.1λ. Over the interval 19λ ≤ t ≤ 20λ the relative change in the stress in the domain is O(10 −6 ), which we assume is sufficient to have reached the steady state solution. In In this subsection we will be studying the numerical stability of the proposed techniques. We will use the meshes M0-M2 of Table 1 for that, which already show the large difference in numerical stability of the two approaches: the Fapproach and the B -approach.
We study a model with a damping function having a = 0.0001. No corresponding differential model is available for comparison. The viscosity function is given in Fig. 11 and the planar extensional function in Fig. 12. The difference with the UCM model is that the viscosity is now shear thinning and the extensional viscosity is finite for all˙ . It should be noted, that at λγ = 129.8 the shear stress goes through a maximum value (of 24.79G) for this single mode model, which turns up as a slope less than -1 in the viscosity curve. Obviously (generalized) shear rates beyond this value should be avoided, since it would lead to extreme strain-rate localization at the cylinder wall when approaching the steady state flow. The planar extensional stress difference goes monotonically to a constant value (of 10000G), and therefore the slope goes to −1 for large extensional rates. Note, that the vertical asymptote for the UCM model is at λ˙ = 0.5.
We study the stability for higher Wi of the two approaches using B and F . The discretization in τ -direction is: N = 50, ∆τ 0 = 0.06λ, τ c = 15λ, f = 1.05546. The time step used is ∆t = 0.1λ. The results are given in Table 3.
It is clear that the limiting Wi is quite low for the B -approach for this model and mesh-refinement does not seem to improve the stability limit. It should be noted, that the flow seems to become steady for a while and then suddenly there is a catastrophic crash (at t = 12.9λ, 18.    For simplicity, we want to restrict the memory function to a single discrete mode. In order to obtain a stable steady flow for high Wi for this case, the damping parameter a needs to be increased to lower the stress growth in elongation. At the same time, this also results in a lower shear rate for the stress maximum in shear, invoking severe stress localisation for lower Wi. Therefore, we need to introduce a solvent viscosity η s to remove the stress maximum altogether. A broader spectrum (e.g. using multiple discrete modes) is also a possibility to achieve that, but this is beyond the scope of the current paper.
We will use a damping parameter of a = 0.01 and a solvent content of β = η s /η 0 = 0.5. The shear viscosity and planar extensional viscosity functions are plotted in Figures 14 and 15, respectively.
We solve the flow with the F -approach using the meshes in Table 2  requires 4000 steps to reach t = 10λ, which can be quite expensive for refined meshes. The time step limitation due to the strain is related to our choice to take the L · F term explicit in the time discretization.
For not too high values of Wi (we tested up to Wi=10) the flow becomes steady with relatively smooth stress fields, confirming the stability of the Fapproach, resembling the stability of the log-conformation representation in differential models [25]. For high Wi (we tested Wi = 100), the flow becomes close to steady after a few relaxation times. Very small stress fluctuations remain, which are probably related to the small scale spatial stress variations remaining present. We did not study the high Wi case further, because it is way beyond the validity of a model with a single discrete mode.
To test the convergence with mesh refinement, we consider the case of Wi = 5. The flow is computed using a time step of ∆t = 0.05λ using 200 steps (end time is 10λ). After starting the flow, there is a stress build-up within 1-4 relaxation times and then a transition to a steady state. This can be seen in In Fig. 17 we show the distribution of τ xx in the steady state. Two boundary layers, one on the cylinder and one on the wall, and a long stress wake are present. In Fig. 18 we show the value of τ xx on the center line and on the cylinder wall for all the meshes used (M0r-M4r). We clearly see mesh convergence on the cilinder wall. However, convergence in the wake has not been achieved yet even for the very refined mesh Mr4. This confirms that convergence in the wake is difficult for high Weissenberg numbers, similar to results found in [25] for a Giesekus fluid. The stress development around the cylinder has an effect on the velocity near the cylinder. In Fig. 19 the velocity component u x /U is shown on the center line for all the meshes used (M0r-M4r). Also here convergence has not yet been achieved in the wake. Due to the viscoelastic stresses, there  are velocity overshoots both in front and in the wake of the cylinder, where the latter is more significant. Note, that the Newtonian value of u x /U in a channel without the cylinder present would be 1.5, hence a little flattening of the profile is seen due the mild shear thinning.

Discussion and conclusions
In this paper we have 'upgraded' the deformation fields method in [2]  methods for differential viscoelastic fluid models. Flows in standard geometries (time-dependent shear flow between two plates and flow around a cylinder confined in a channel) have been simulated to study the stability and accuracy of the new implementation of the deformation fields method. The following techniques/methods have been introduced to obtain the improved stability: • Spatial discretization for the deformation fields is using SUPG for stabilizing convection in real space (x).
• The convection in real space and in 'fluid history' space (τ ) is taken implicitly to allow for larger time steps.
• A novel scheme, based on a differential form of the integral model, that decouples the momentum and mass balance from the constitutive model, allowing for low or zero-solvent viscosity.
• For obtaining the deformation history the 'projected' velocity gradient G is used instead of the (∇u) T , leading to better transient stability.
• The deformation history is represented (and solved) by using the deformation gradient F instead of the (Finger) deformation tensor B.
The first four techniques are more or less extensions to integral models of well known similar techniques for differential models. The last technique is different and calls for some more discussion. It basically replaces the log-conformation representation [12] (LCR) for differential models and gives similar stability. The reason for the stability is the same as the remarkable stability of stochastic models based on the method of Brownian configuration fields (BCF) versus conformation tensor models (without LCR). In [26] it has been argued that the stability of BCF comes from the formulation of the models in terms contravariantly deforming vectors (Q): ∂Q ∂t = −u · ∇Q + L · Q + . . .
In elements near stagnation points (u = 0) the first term on the right-hand side (convection) balances the second term (deformation) quite naturally and numerical stability is obtained even on coarse grids. Conformation tensor models (without LCR) are formulated in terms of contravariantly deforming tensors: ∂c ∂t = −u · ∇c + L · c + c · L T + . . .
In elements near stagnation points the first term on the right-hand side (convection) balances the second term (deformation), however there is an additional deformation term leading to numerical instability for underresolved grids. Since the deformation gradient F behaves like contravariantly deforming vectors (see Eq. (9)) and the Finger tensor B like contravariantly deforming tensors (see Eq. (12)), we can expect much better stability in the F -approach than in the B -approach. Note, that when applying the LCR [12] both deformation terms are basically 'neutralized', but from a stability point of view, this is not needed.
An additional benefit of the 'vector' based methods (BCF, F -approach) versus the 'tensor' based methods (conformation tensor, B -approach) is the guaranteed positive-(semi)definiteness of the stress tensor. Since positive-definiteness is required in the physical model, instabilities due to negative eigenvalues of the stress tensor can be avoided. Note, that LCR also leads to a positive definite stress tensor.
The next step is extending the range of constitutive models. In particular, we are interested in integral models with long-range memory functions, such as powerlaw functions [27] and powerlaw-like functions as appear in fractional K-BKZ models [5]. This will be part of future work. where we transformed the memory integrals from the domain t ∈ (−∞, t] to τ ∈ [0, ∞). Note, that we have also splitted the material derivative of the stress tensor into a local and convective part: 10) since in this paper we solve all equations in the Eulerian frame.