A Backwards-Tracking Lagrangian-Eulerian Method for Viscoelastic Two-Fluid Flows

: A new Lagrangian–Eulerian method for the simulation of viscoelastic free surface ﬂow is proposed. The approach is developed from a method in which the constitutive equation for viscoelastic stress is solved at Lagrangian nodes, which are convected by the ﬂow, and interpolated to the Eulerian grid with radial basis functions. In the new method, a backwards-tracking methodology is employed, allowing for ﬁxed locations for the Lagrangian nodes to be chosen a priori. The proposed method is also extended to the simulation of viscoelastic free surface ﬂow with the volume of ﬂuid method. No unstructured interpolation or node redistribution is required with the new approach. Furthermore, the total amount of Lagrangian nodes is signiﬁcantly reduced when compared to the original Lagrangian–Eulerian method. Consequently, the method is more computationally efﬁcient and robust. No additional stabilization technique, such as both-sides diffusion or reformulation of the constitutive equation, is necessary. A validation is performed with the analytic solution for transient and steady planar Poiseuille ﬂow, with excellent results. Furthermore, the proposed method agrees well with numerical data from the literature for the viscoelastic die swell ﬂow of an Oldroyd-B model. The capabilities to simulate viscoelastic free surface ﬂow are also demonstrated through the simulation of a jet buckling case.


Introduction
Viscoelastic free surface flows are of significant importance for many industrial processes.This includes polymer extrusion, additive manufacturing, seam sealing, and adhesive joining.In such processes, the viscoelastic properties of the flow, such as the flow history, can have major influence on the quality of the final product.Furthermore, production time and raw materials may be subject to a large cost.Therefore, extensive manual effort and physical testing may be necessary to optimize the process in terms of product quality, material consumption, and production cycle time.
Numerical tools can be helpful in reducing the manual preparation time, as they offer the possibility to replace a significant part of the physical testing with computer simulation.Furthermore, numerical simulations may enable testing early in the design phase.For many complex industrial applications, a large demand for suitable numerical simulation tools therefore exists.Furthermore, for such tools to be useful, high demands are put on accuracy, robustness, and computational efficiency.
A common approach for simulating viscoelastic fluid flow is to discretize the governing equations in the Eulerian frame of reference with the finite volume method (FVM) [1][2][3] or the finite element method (FEM) [4,5].The Eulerian frame of reference is suitable for diffusion-dominated problems, including viscous fluid flow as well as heat and mass transfer.On the other hand, viscoelastic constitutive equations are typically hyperbolic and they involve no physical diffusion term.Thus, in this sense, the Lagrangian frame of reference is a more natural description for the constitutive equation.Therefore, Lagrangian or semi-Lagrangian methods are appealing as alternative approaches to the Eulerian frame.
Rasmussen and Hassager [6] developed a Lagrangian finite element method in order to calculate the flow of an upper convected Maxwell (UCM) fluid.The entire flow history was stored and re-meshing was required, due to deformation.Harlen et al. [7] proposed a split Lagrangian-Eulerian method for viscoelastic Stokes flow.The constitutive equation was solved at the nodes of a co-deforming mesh, while the momentum and continuity equations were solved with an Eulerian finite element method.Re-meshing due to mesh distortion was also used in their method.Halin et al. proposed the Lagrangian particle method (LPM) [8].In their method, the momentum and continuity equations were solved with the Eulerian finite element method.The constitutive equation was solved along the trajectories of massless Lagrangian particles.The stresses in the particles were then fitted to local polynomial expressions in order to enable the evaluation of the corresponding finite element integrals in the momentum.Hence, a minimum of three particles per twodimensional element was required for the polynomial approximations to be feasible.Later on, the adaptive Lagrangian particle method (ALPM) [9] was proposed, utilizing adaptive addition and the deletion of particles.However, a fairly large number of particles was required for the simulations to produce stable transient results.Furthermore, Wapperom et al. [10] proposed the backward-tracking Lagrangian particle method (BLPM), in which backwards integration of the velocity field allowed for the choice of fixed particle locations a priori, increasing the efficiency of their particle tracking approach.
The simulation of viscoelastic free surface flow has been approached with different methods.For example, Crochet and Keunings [11] simulated the viscoelastic die swell effect of an Oldroyd-B fluid with a mixed finite element method as early as 1982.More recent similar examples of finite element methods on non-stationary grids can be found in Balemans et al. [12] and Spanjaards et al. [13].The marker-and-cell (MAC) approach is a method that has been used quite extensively for simulation free surface flows [14][15][16][17][18][19][20], in which the governing equations discretized with finite differences and the free surface is tracked while using marker particles.The method has been successfully used for the simulation, for example, of the die swell and viscoelastic jet buckling with a variety of constitutive models.The simulations are typically performed on staggered, uniform grids, which has certain limitations for complex simulation geometries and in terms of computational efficiency when high resolution is required.Furthermore, the existence of a gas phase surrounding the viscoelastic fluid is not taken into account, which may be required for some cases.The front-tracking method is another similar method [21].
The volume of fluid (VOF) method is a popular method for the simulation of Newtonian as well as viscoelastic free surface flow with finite volume discretization.VOF is a diffuse-interface method, in which the presence of two or more fluid phases is represented by their corresponding volume fractions, advected with geometric or algebraic schemes.A sharp fluid interface may then be reconstructed from the solution.A single set of momentum, continuity equations are then solved for the whole domain, and the fluid properties are locally averaged with the fluid volume fraction.Furthermore, the transport of the volume fraction may be solved with the same spatial grid as the momentum and continuity equation, which makes the method suitable for finite volume simulation of free surface flow.Some examples follow.
Habla et al. [22] developed a VOF-solver for viscoelastic free surface flow that is based on the open source solver OpenFOAM [23].The simulations were stabilized by reformulating the viscoelastic constitutive equations with the log-conformation representation (LCR) and with both-sides diffusion (BSD) in the momentum equation.Comminal et al. [24] simulated the viscoelastic die swell effect for Oldroyd-B and Giesekus fluids with the VOF method, combined with two different schemes for the convection of the fluid volume fraction, a geometric scheme, as well as an algebraic scheme that directly solved the transport equation for the volume fraction.Furthermore, their methods were compared to the VOF solver that is available in RheoTool [25], which is an open source toolbox for viscoelastic flow that is based on OpenFOAM.The constitutive equation was also reformulated with LCR in their study.Bonito et al. [26] reported a related method, using a mathematical formulation that is similar to the volume of fluid method, while solving the governing equations with a piecewise linear finite element method, contrary to finite volume discretization.
In previous work, a Lagrangian-Eulerian framework for viscoelastic flow has been proposed [27,28].In the method, the constitutive equation is solved at Lagrangian nodes, which are convected by the flow.The momentum and continuity equations are discretized with the finite volume method on an Eulerian octree grid and then solved with SIMPLEC iterations [29].Boundary conditions that are imposed by objects in the domain are treated using implicit immersed boundary methods [30,31].The contribution to the momentum equation from the viscoelastic stress is established through the unstructured interpolation of the stress tensors from the Lagrangian nodes to the cell centers of the Eulerian grid while using the radial basis function.
The ability to efficiently simulate industrial applications, such as adhesive extrusion, parts assembly, and additive manufacturing, is a main motivator to develop the new method.For this purpose, it is necessary that the numerical method is accurate, computationally efficient, and robust.The ability to model viscoelastic free surface flow in complex geometry with moving objects is also necessary.The proposed method was shown to produce results in good agreement with available analytic and numerical data from the literature [27].It has also been shown that the Lagrangian formulation of the constitutive equation allows for the efficient parallelization of the viscoelastic stress calculation, which makes the method highly suitable for GPU-acceleration [28].Furthermore, the use of an immersed boundary method in combination with the automatically generated grid enables simulation on arbitrary geometry with minimal pre-processing.
In the previous studies [27,28], the proposed Lagrangian-Eulerian method has been validated for single-phase flows with uniform grids.The natural next steps include enabling the simulation of free surface flow as well as simulation on refined grids.Therefore, in the current work, the proposed method is further developed in these two aspects.Firstly, the framework is extended for simulation of viscoelastic free surface flows with the volume of fluid method.Secondly, a backwards-tracking procedure is introduced to the solution of the constitutive equation, partly inspired by the backward-tracking Lagrangian particle method by Wapperom et al. [10].The result is a more robust method, due to the structured relation between the Lagrangian nodes and the Eulerian grid.No stabilization method, i.e., neither both-sides diffusion nor reformulation of the constitutive equation, is required for the studied flows.Furthermore, the computational cost as compared to the original Lagrangian-Eulerian method is reduced, since no unstructured interpolation or node redistribution is needed and because the total number of Lagrangian nodes decreased.In the context of multiphase flows, the viscoelastic constitutive equation is only solved inside the viscoelastic phase, which further increases the computational performance.
The rest of the paper is structured, as follows.In the next section, the governing equations are presented, followed by a detailed description of the proposed numerical method.In the results section, the new method demonstrated and validated with analytic solutions and numerical data from the literature.Finally, conclusions are drawn and future work is outlined.

Governing Equations
Viscoelastic fluid flow is described by the incompressible momentum and continuity equations where ρ is the density, u velocity, µ the solvent viscosity, S = 1 2 (∇u + (∇u) T ) the strain rate tensor, τ the viscoelastic stress tensor, and f a body force.The viscoelastic stress is governed by a constitutive equation of the general form where λ is the relaxation time, η polymeric viscosity, and F a scalar-valued function.The operator (•) denotes the upper-convected derivative of τ, which expands to where is the Lagrangian, or material, time derivative of τ.The Lagrangian time derivative describes the rate of change of τ in an infinitesimal material element, which moves with the flow.Convected time derivatives, such as (•), appear in viscoelastic constitutive equations and they ensure frame-invariance for the transport of tensorial properties.It is remarked that, while the form of (3) is not applicable for all available differential constitutive equations, it covers the constitutive models discussed in this work.
In general, the viscoelastic stress tensor may be modeled as a sum of the contribution of N modes as where τ k is the stress tensor that corresponds to the kth mode.Each stress tensor τ k is then governed by a constitutive equation on the form of (3) with an individual set of model parameters.

Numerical Method
The momentum Equation (1) and the continuity Equation ( 2) are discretized on a collocated Eulerian grid with the finite volume method.The pressure-velocity coupling is accomplished by solving the momentum and continuity equations with SIMPLEC iterations [32].An octree grid is used, in which the cells can be refined in desired locations by recursively dividing them into smaller cells to a given refinement level.The refinements may be static or adaptively updated, typically around moving objects or fluid interfaces.Figure 1 displays an example grid in two dimensions with one refinement level, showing the cell centers and the grid nodes.The grid, including the adaptive refinements, is automatically generated by the solver.Boundary conditions from objects in the domain are imposed while using the mirroring immersed boundary method [30,31].The velocity field is implicitly mirrored across the immersed boundary surface, such that the imposed boundary condition is satisfied for the converged solution.Surface triangulations represent the geometrical objects.
The viscoelastic constitutive equation is solved in the Lagrangian frame of reference in material elements that are represented as Lagrangian nodes.From (3) and ( 4), it follows that the stress in a Lagrangian node is described by the ordinary differential equation (ODE) Furthermore, the trajectory of the node must be known when solving (7), as it involves the local velocity gradient ∇u .The trajectory of the node is given by the ODE where x is the position of the node.Given an initial stress state in a set of Lagrangian nodes, a straightforward approach to obtain the stress state at a future time is to simultaneously solve (7) and ( 8) forwards in time.Indeed, this method has been successfully employed in previous work [27,28].Lagrangian nodes were then initialized with a given distribution density at the start of the simulation.The distribution was maintained in each time step by adding or deleting nodes, if necessary.After solving the constitutive equation at the Lagrangian nodes, the viscoelastic stress was interpolated to the Eulerian cell centers while using radial basis functions (RBF).
In this work, the method is further developed by introducing the concept of backwardstracking.The key idea is to choose the locations of the Lagrangian nodes a priori, at which the viscoelastic stresses are then stored, and then track them backwards in time in order to obtain their initial location in each time step.The constitutive equation may then be solved forwards in time such that the stresses are obtained at for the next time step at the chosen locations.Different choices for the locations are conceivable.In the current work, the locations of the Lagrangian nodes are chosen to be at the Eulerian grid nodes.
Consider the calculation of the viscoelastic stress for the nth simulation step, which corresponds to the time interval I n = [t n , t n+1 ].At this point, it is assumed that momentum equation has been solved, such that the velocity field is known for t ∈ I n .
The Lagrangian constitutive equation is solved along the trajectories of Lagrangian nodes.For a Lagrangian node, the trajectory x(t) can be expressed in terms of the position at the end of the time step, as where x(t n+1 ) is the predefined final location of the Lagrangian node which.The trajectory x(t) for t ∈ I n is first calculated by solving (8) backwards in time.Subsequently, when the trajectory is known, the viscoelastic stress τ(t n , x(t n )) is interpolated from the viscoelastic stresses from the solution of the previous simulation step.The Lagrangian constitutive Equation ( 7) is then solved forwards in time and the viscoelastic stress is thus obtained at the final location x(t n ), i.e., at the predefined location.Generally, the choice of locations for the Lagrangian nodes is somewhat arbitrary.However, by choosing the locations of the Eulerian grid nodes, the connectivity and structure of the octree grid may be utilized when the stress at time t n is interpolated to the position of a Lagrangian node x(t n ).The stress at the node at time t n is obtained through bilinear or trilinear interpolation, respectively, for two-dimensional (2D) and three-dimensiona (3D).Furthermore, the viscoelastic stress contribution to the discretized momentum equation may be calculated directly from node stresses.Therefore, no unstruc-tured interpolation of viscoelastic stresses is necessary.The calculation of this contribution is described in detail further on in this section.
Figure 2 shows a schematic description of the steps that are involved in the algorithm.In summary, the performed steps are: (a) Calculate the Lagrangian node trajectory by solving (8) backwards in time, starting at the predefined location x(t n+1 ).(b) Interpolate the stress τ(t n , x(t n )) to the Lagrangian node from the known node values at time t n .(c) Solve (7) forwards in time along the trajectory x(t), t ∈ I n .
x(t n+1 ) The ordinary differential Equations ( 7) and ( 8) may be solved with any appropriate choice of ODE solution algorithm.In this work, the fourth order Runge-Kutta method, which is commonly referred to as the RK4 method [33], is used.Furthermore, N sub equally sized substeps of length ∆t sub are used, being defined, such that N sub ∆t sub = ∆t = t n+1 − t n .In the current work, N sub = 3 is used.
When solving (7) and ( 8), local quantities that are stored on the Eulerian grid are required along the trajectory of the Lagrangian nodes.More specifically, the velocity u is needed for the backwards-tracking and the velocity gradient ∇u is needed to solve the constitutive equation.The velocities are stored at the Eulerian cell centers.When the local velocity is required, the cell centers containing the Lagrangian node are identified and the velocity is interpolated to its location while using bilinear or trilinear interpolation, respectively, for 2D and 3D.When the velocity gradient is required, it is calculated from the bilinear or trilinear interpolation formula using the same interpolation basis.
For the interpolation of a field φ to a location within a box spanned by the corners (x i , y j , z k ), i, j, k ∈ {0, 1}, the trilinear interpolation formula in order to calculate the interpolant φ reads where φ ijk is the value of φ at the corner (x i , y j , z k ) and with the coefficients where xl is the lth coordinate of the interpolation position relative to the lower corner of the box and ∆x l is the size of the box that coordinate direction.The gradient of φ follows from ( 10) and ( 11), as where the derivatives of the coefficients read The corresponding formulas for bilinear interpolation in two dimensions are obtained by letting the coefficient c k z = 1 and, thus, ∂c k z /∂z = 0, and skipping the summation over k.When interpolating to a Lagrangian node near fluid grid refinements the smallest local cell size is used in order to form the interpolation box.Furthermore, not all corners of the interpolation box coincide with a cell center.The velocities at such corners are obtained through least squares interpolation, by fitting a first order polynomial from the centers of the cells that intersect the interpolation box.The interpolation box is visualized in Figure 3.
The discrete counterpart of ξ is the fluid volume fraction α ∈ [0, 1], which is the local volume average of ξ in an Eulerian cell.The transport of α is described by the convection equation [34] ∂α ∂t By definition, α is the local volume fraction of the viscoelastic phase in a control volume.Hence, if α = 1, the control volume is filled with the viscoelastic fluid and if α = 0 with the Newtonian fluid.If α ∈ (0, 1), this is interpreted as that the control volume intersects the interface between the fluid phases.The transport Equation ( 17) is discretized with the finite volume method on the Eulerian octree grid.In order to minimize the numerical diffusion of α and, thus, avoid smearing the fluid interface, the convection term in ( 17) is discretized with the compact CICSAM scheme [35].
A single set of the Equations ( 1)-( 3) is solved for the whole computational domain.Fluid properties are locally averaged as where φ v and φ N are the properties of the viscoelastic and the Newtonian phase, respectively.The averaging is applied for ρ, µ, η, and λ, as well as for τ when calculating the contribution to the momentum equation.When α is needed at a Lagrangian node, it is interpolated with the same method as used for velocity.Special care for the viscoelastic stress in the interface region α ∈ (0, 1) is necessary due to the non-sharp interface nature of the VOF method.A brief motivation for this follows.In many examples of viscoelastic free surface flow, the viscosity of the Newtonian fluid may very well be on the order of 10 −5 of that of the viscoelastic phase.Hence, the average velocity gradient in the Newtonian phase may be much larger than in the viscoelastic phase.Consequently, if the velocity gradient of the Newtonian phase is allowed to have large influence on the constitutive equation in the interface region; this can result in unphysically large stresses.This is especially the case if not all velocity and length scales in the Newtonian phase are resolved near the interface.Thus, a means of reducing this effect around the free surface is proposed.
A threshold volume fraction is used, such that Lagrangian nodes are only created inside the viscoelastic phase, which are defined by the condition α ≥ α lim = 0.1.Thus, the constitutive equation is only solved in the viscoelastic part of the computational domain and the stress is assumed to vanish outside the viscoelastic phase.Furthermore, when the velocity gradient is calculated at a Lagrangian node, the corners of the interpolation box at which α < 0.1 • α lim are also considered to be outside the viscoelastic phase.The velocities at such corners are excluded from the gradient calculation and they are replaced by 0th order extrapolation from the corners which are inside the viscoelastic phase.
It is remarked that, apart from numerical stability, the described approach is advantageous in terms of computational efficiency, since the constitutive equation only needs to be solved in part of the domain.
When the stress field has been calculated, the term ∇ • (ατ) is integrated and then added to the right hand side of the discretized momentum equation.At this stage, the product rule is applied, such that which can be seen as separating the pure interfacial contribution of the viscoelastic stress to the fluid momentum from the remainder part [36].The formulation of ( 19) is used, since it was found to enhance numerical stability.The first term of ( 19) is integrated over the cells with Gauss' divergence theorem, as c.v.
where c.v. denotes the cell volume, c.s. the cell surface, f.s. the surface of cell face f , and the vector n denotes the surface normal pointing outwards from the cell.In the second step of (20) the surface integral is divided into a sum of the integrals over the respective surface faces of the cell, for which the normal vectors n f are constant.The integral over each face is approximated with the trapezoidal rule while using the stresses at the Eulerian grid nodes.If a cell has neighbors that have a higher refinement level, e.g. in the case shown in Figure 1, the face to each smaller neighbor cell is individually integrated, such that the stress contribution from each grid node is included correctly.The volume integral of the second term of ( 19) is approximated with the cell average stress and volume fraction gradient, as c.v.
where ∆V is the cell volume and (•) denotes volume average.The cells where α < 0.1 • α lim are assumed to be outside the viscoelastic phase and the stress contribution is instead set to zero.The full algorithm for viscoelastic free surface flow can be summarized, as follows, 1.
Calculate and add the viscoelastic stress contribution to the discrete momentum equation.

2.
Solve velocity and pressure fields u, p from the momentum Equation ( 1) and continuity Equation (2) while using SIMPLEC iterations.

3.
Solve the transport of the fluid volume fraction α from (17) 4.
Solve the viscoelastic stress from the constitutive Equation (3) to the new time while using the backwards-tracking procedure.
When compared to the original version of the Lagrangian-Eulerian framework that was proposed in Ingelsten et al. [27], the new method is expected to be more robust due to the structured nature of the Lagrangian nodes as well as the staggered arrangement between the storage locations of the velocity and the viscoelastic stress.Although the computational efficiency is not assessed in detail in this work, the computational cost is, in fact, reduced.As a reference, the computational performance of the original Lagrangian-Eulerian method was studied for the flow of a four-mode PTT fluid over a confined cylinder [27].Four Lagrangian nodes per fluid cell were initialized for the node set.Out of the total stress calculation time, approximately 50% was spent on solving the ODE systems, 30% on the unstructured interpolation, and 10% on the redistribution of the Lagrangian node set.In the new method, the unstructured interpolation, as well as the redistribution, is completely removed.Furthermore, the number of Lagrangian nodes is reduced to be on the same order as the number of Eulerian cells.Thus, a large reduction in computational time should be expected.
The proposed method is implemented in the software IBOFlow ® [37], an in-house CFD code that was developed at the Fraunhofer-Chalmers Research Institute for Industrial Mathematics in Gothenburg, Sweden.In addition to viscoelastic flow, the solver has previously been employed to simulate conjugated heat transfer [38][39][40], and fluid-structure interaction [41], as well as free surface flow of shear-thinning fluids with applications for seam sealing [42,43], adhesive application [44], and 3D-bioprinting [45].

Results
In this section, the proposed method is evaluated for three different flow cases.First, a basic validation is carried out by comparing the numerical results with the analytic solution for a transient channel flow.The method is then compared to the numerical results from the literature for simulations of the viscoelastic die swell effect.The capability to simulate viscoelastic free surface with adaptive mesh refinement is demonstrated for simulations of viscoelastic jet buckling.
The test cases that are reported in this section are simulated while using a single-mode Oldroyd-B model, which has the constitutive equation corresponding to F(τ) = 1 in (3).Furthermore, a viscosity ratio β is defined as By definition of the total viscosity η t = µ + η, it follows from (23) that µ = βη t and η = (1 − β)η t .Note that the Lagrangian-Eulerian method that is proposed in this work is by no means limited to the Oldroyd-B model.All of the constitutive equations in the form of (3) are supported by the framework.

Planar Poiseuille Flow
A basic validation of the proposed method is performed for a viscoelastic planar Poiseuille flow.The numerical results for the transient as well as the steady flow are compared to the corresponding analytic solution.
The computational domain has height and length h and it is shown in Figure 4.The upper boundary is treated as a wall with the no-slip condition u = 0 and the lower boundary has a symmetry condition.At the left and right boundaries, the Dirichlet conditions p = ∆p and p = 0 are imposed, respectively, while the cyclic conditions are used for the velocity and viscoelastic stress.Thus, the numerical model represents an infinitely long channel of height 2h, which is subjected to a constant pressure drop ∆p.The simulation starts from rest, with velocity and viscoelastic stress equal to zero.A Weissenberg number for the flow is defined as Wi = λU/h, where U is the mean steady flow velocity, and a Reynolds number as Re = 2ρUh/(µ + η).The flow is simulated for β = 1/9, 1/18, 1/27 and Wi = 0.1, 1.Here, U = 0.1 m/s, ρ = 1 kg/m 3 and η = 1 Pas are constant, while Wi and β are varied, respectively, by changing λ and µ between simulations.
In order to ensure sufficient temporal resolution, the flow is simulated using different time step lengths.These simulations are performed for the lowest viscosity ratio β = 1/27, for which the transient variations are the largest.A uniform Eulerian grid with the cell size ∆x = H/20 is used.
In Figure 5, the simulated centerline velocity, i.e., at x = h/2, y = 0, obtained with the time step lengths 10 −4 s, 10 −5 s and 10 −6 s, are shown for Wi = 0.1 and Wi = 1.The velocity that si obtained with the longest step length is clearly different from those obtained with the two shorter step lengths.However, the two shorter lengths produce results that practically overlap on the scale of comparison.Thus, the step length ∆t = 10 −5 s is considered to be sufficiently small and it is used to obtain the remaining results reported in this section.[46].In Figure 6, the obtained centerline velocities for the β = 1/9, 1/18, 1/27 are shown for Wi = 0.1 and Wi = 1 and compared to the analytic solution.The results demonstrate the strong influence of the relationship between the solvent viscosity and polymeric viscosity on the transient flow.The numerical results overlap the analytic solution, showing that the proposed method treats the transient flow dynamics correctly.Furthermore, the results confirm that the spatial resolution ∆x/h = 20 is sufficient for this case.In order to assess grid dependency in more detail, the flows for β = 1/27 and Wi = 0.1, 1 are simulated until a steady flow is reached while using varying spatial resolution.The flow being steady is ensured through the condition where ||•|| L 2 denotes the L 2 -norm over the cells of the Eulerian grid, φ is velocity or viscoelastic stress, and the subscript n denotes the quantity at the nth time step.The error measurements with respect to the analytic steady state solution are calculated as where φ s and φ a denote the simulated and analytic solution, respectively.The computed errors for velocity and viscoelastic normal stress are shown as a function of grid size in Figure 7.The errors decrease with second order slope with grid refinements, which is coherent with the second order accuracy that was observed for the original Lagrangian-Eulerian method [27,28].

Die Swell
The die swell effect is a phenomenon that occurs for the free surface flows of viscoelastic fluids flowing in a pipe or channel and exiting through circular or slit dies.Barnes et al. [47] described the origin of the die swell effect by viewing the viscoelastic fluid as a bundle of elastic threads.In the channel flow, the threads are stretched by the normal stress component τ xx , where x is the flow direction.When the fluid exits the channel, the threads are allowed to relax and shorten in length.Consequently, the diameter of the emerging fluid increases.
Flows that exhibit the die swell effect are commonly used to benchmark and validate numerical methods [11,13,14,[16][17][18][19][20]22,24].In addition to the free surface, the flow has extensional characteristics around the channel exit as well as a stress singularity at the channel exit corner.These flow features make the flow suitable for testing robustness and accuracy of new numerical methods.Therefore, in this work, a planar die swell flow is simulated.The case parameters are selected to allow for a meaningful comparison with numerical data from the literature.
A schematic description of the computational domain can be found in Figure 8.A channel of height 2h and length 10h is considered.A symmetry boundary condition is used at y = 0, such that effectively half of the domain is simulated.The expansion zone is 5h high and 12h long, giving the domain a total length of 22h.The length of the channel and size of the expansion zone are chosen to be sufficiently large for boundary effects not to affect the flow near the channel exit.
The exterior boundaries of the expansion zones are treated as outlets with the pressure condition p = 0, and an immersed boundary is used to impose the no-slip condition for the channel walls.Gravity is neglected.The viscoelastic fluid in the channel is an Oldroyd-B fluid with the viscosity ratio β = 1/9.The surrounding fluid is a Newtonian fluid with much lower viscosity and density than the viscoelastic fluid.In this work, these parameters are set to µ air = 10 −6 µ and ρ air = 10 −3 ρ, respectively.Comminal et al. used a similar definition [24], but with ρ air = 10 −2 ρ.

Outlet
The simulation of the flow starts from rest, with the channel filled with viscoelastic fluid.Transient flow is simulated for a sufficiently long physical time for the flow in the channel and the expansion region to fully develop, as well as for the free surface flow to exit through the outlet with a uniform velocity profile.
At the inlet, fully developed flows are imposed for the velocity and viscoelastic stress.For fully developed channel flow of an Oldroyd-B fluid, it can be shown that where U is the mean velocity.A derivation of the expressions can be found in Appendix A.
The channel half height h and velocity U are taken as characteristic length and velocity scales for the flow.Thus, a Weissenberg number for the flow may be defined as and a Reynolds number as A dimensionless time is also defined as t * = tU/h.The swell ratio S r = h max /h is the main quantity of interest, in terms of validation by comparison.The swell ratio is a function of the so-called recoverable shear S R , which is defined as [11] where N 1 = τ xx − τ yy is the first normal stress difference.In the simulations that are reported in this section, the quantities U, η, µ, and ρ are constant, while S R is varied by changing λ between simulations.The Reynolds number is Re = 0.5 for all simulations.
Tanner [48,49] developed a theoretical solution for the amplitude of the die swell as a function of the recoverable shear.For an Oldroyd-B channel emerging from a slit the solution reads The theory does not take the stress singularity at the channel exit corner into account and is, therefore, only expected to yield a good approximation of the swell ratio when the influence of elasticity is small, i.e., for small Wi.Nevertheless, it is often included for comparison in numerical studies of the die swell effect and it is included in this study.
The influence of the spatial discretization on the numerical results is assessed through a grid dependence study and the flow is simulated for S R = 2.5 with three different grids, denoted as M1, M2, and M3.A grid is defined by a base cell size ∆x base and a set of refinements, as described in Section 3. The grids M1, M2, and M3, respectively, have the base cell sizes ∆x base = h/5, h/10, h/20.The largest cells are located in the expansion zone, far from the channel exit.One level of refinement is used in the channel and around the channel opening and two levels of refinement near the exit corner.As an example, grid M1 is shown in Figure 9.The area around the channel exit corner has been zoomed for clarity.Table 1 summarizes the grids used.A constant time step length is used, such that U∆t/∆x min = 0.01, where ∆x min is the smallest cell size of the grid.1.
First, the velocity as well as the viscoelastic normal and shear stress are compared to the analytic steady solution inside the channel.In Figure 10, the corresponding profiles obtained using grid M1 across the channel at x/h = 5, halfway from the inlet to the channel exit.Already at the lowest grid resolution, the simulated profiles overlap the analytic solution.The results that are obtained with M2 and M3 are visually the same and the figures are, therefore, omitted.The results validate that the constitutive equations are solved correctly and that the immersed boundary method used to impose the boundary conditions at the channel wall works as intended for the solution algorithm for the viscoelastic stresses.In Figure 11, the position of the free surface that is visualized by the contour α = 0.5 is shown for S R = 2.5 for the grids that are defined in Table 1.Some small surface oscillations are observed in the results from the two finer grids M2 and M3 in the region x/h < 2. As reference, Comminal et al. [24] reported similar behavior for their corresponding simulations.In their simulations, small self-sustained surface oscillations appeared at a certain level of recoverable shear, in their case S R > 1.5 for the Oldroyd-B fluid, and for sufficiently high resolution.Similar observations were also made for simulations with the Giesekus model in their study.The oscillations were damped out further downstream and they did not influence the calculated swell ratio h max /h.Similar characteristics are observed in the current work, as in Figure 11.Comminal et al. attributed the oscillations to numerical difficulties around the free surface due to the nature of the non-sharp interface in the VOF method.
In the light of this discussion, the results that are obtained on grid M2 and M3 are in relatively good agreement.This is particularly true far downstream of the channel exit, where the free surface produced by M2 and M3 are very close, while M1 produces a different result.The results indicate that M2 is of sufficient spatial resolution to predict the swell ratio S r .Therefore, the results presented in the remainder of this section have been obtained while using grid M2.In addition, the simulation with grid M2 has been repeated with half the time step size, such that U∆t/∆x min = 0.005, which produced equivalent results.The simulations are repeated with grid M2 for S R ∈ {1.0, 1.5, 2.0, 2.5}.In Figure 12, the contour α = 0.5 is shown for S R = 1 at different times.Figure 13 shows the corresponding results for S R = 2.5.Similar behavior is observed in both cases.The emerging viscoelastic fluid increases in diameter directly upon exiting the channel.The magnitude of the increase is strongly influenced by the level of elasticity, as qualitatively predicted by Tanners theory.At a certain distance from the channel exit, approximately for x/h > 3, the emerging viscoelastic fluid reaches a terminal diameter.
The swell ratio is calculated for the simulations as S r = h max /h, where h max is taken as the position of the free surface at x/10.The calculated swell ratios that were obtained with the proposed method are compared to available numerical data from the literature.The data used are the original FEM simulations reported by Crochet and Keunings [11], the marker-and-cell simulations by Tomé et al. [14], the pseudo-VOF simulation by Habla et al. [22], and the three different VOF methods used by Comminal et al. [24].Table 2 briefly summarizes the data.For a full description of the cases and numerical methods, the reader is referred to the respective studies.However, a brief discussion concerning the differences between the results is given below.Table 2. Summary of numerical data as compared for the die swell flow.

Work Method Re
Current work VOF, Lagrangian-Eulerian 0.5 Crochet & Keunings [11] Mixed FEM 0 Tomé et al. [14] GENSMAC 0.5 Habla et al. [22] pseudo-VOF 0.5 Comminal et al. [24] (CCU) VOF, Geometric scheme 0 Comminal et al. [24] (HRIC) VOF, Algebraic scheme 0 Comminal et al. [24] (RheoTool) VOF, Algebraic scheme (MULES) 0.01 The calculated swell ratios as a function of the recoverable shear for the flows simulated with the proposed method are compared to the results from the literature in Figure 14.It is remarked that Comminal et al. employed a different definition of the Weissenberg number and swell ratio.However, all of the data that are presented in Figure 14 have been adopted to the definition used in the current work.
Indeed, a certain spread among the data is observed.The differences are not unexpected, though, as there are slight differences between how the different data have been obtained.This is in terms of the numerical method used as well as the simulation setup.Furthermore, Comminal et al. even obtained different results while using three slightly different numerical methods for exactly the same flow.
Crochet and Keunings and Tomé et al. did not consider a fluid surrounding the viscoelastic phase, as the per construction of their numerical methods.Comminal et al. used a similar definition of the surrounding fluid as in the current work, but with larger density.Habla et al. stated that the surrounding fluid is treated as air with µ air → 0 and ρ air → 0, but the magnitudes were not reported.Furthermore, the value of viscosity ratio of the viscoelastic fluid β varied between their simulations.This is in contrast to a constant β = 1/9, which has been used in the other simulations included in this work.
Given the above discussion, a certain variance in the data is to be expected.However, the different results do follow the same general trend, including the swell ratios predicted with the proposed Lagrangian-Eulerian method.
Another aspect is that of grid dependence.The simulation performed by Crochet and Keunings had a grid resolution which, by today's standards, was very coarse.They used six triangular finite elements across the half-width of the channel.Therefore, the resulting swell ratios should be treated with caution when comparing to more recent results.Tomé et al. used a uniform grid of ∆x/h = 1/10 which is on the order of the coarse grid M1.
Habla et al. reported the number of control volumes to be 4165, but did not specify the local cell sizes in detail.The number of control volumes is again comparable to grid M1 in the current work.Comminal et al. did conduct a grid dependence study with three different grids.Their finest cells were closer to the grids that were used in this work as compared to the other authors.Therefore, could be expected for the results reported in this work to be closest to the results of Comminal et al. given the similarities in grid resolution as well as the use of the VOF method.Their computed swell ratios are generally larger than the other literature results, with the exception of those from Crochet and Keunings.The same is also true for the swell ratios that were computed with the proposed method, particularly for increasing S R .
In conclusion, while there is a spread between the data that are produced in different numerical studies, they follow the same general trend.The results that were obtained with the method proposed in this work follow the same trend and with the swell ratios comparable to the earlier published works.

Jet Buckling
The simulation of jet buckling is a common flow case for testing the capabilities of numerical methods for free surface flow.The phenomenon has been numerically investigated for different viscoelastic constitutive models by several authors [14][15][16][17][18]20,26].
The phenomenon occurs for a fluid jet that flows onto a rigid plate under certain conditions.Cruickshank [50] proposed a condition that is based on experimental and theoretical observations, stating that a planar Newtonian jet will buckle if H/D > 3π and Re < 0.56, where H is the distance from the inlet to the plate, D the inlet width, and Re = ρUD/µ, where U is the inlet velocity.A modified yet approximate condition, based on numerical investigation, was later proposed by Tomé It is remarked that, in simulations, numerical round-off errors are responsible for triggering the buckling, since the flow setup itself is actually symmetric.
Figure 15 shows a schematic of the fluid jet buckling domain that is used for simulation with the Lagrangian-Eulerian method proposed in this work.The domain is 50 mm wide and 100 mm high.An inlet of width D = 5 mm is located at the center of the upper boundary.At the inlet, a uniform velocity U = 0.5 m/s is imposed.The upper boundaries to the left and right, respectively, are outlets with the Dirichlet pressure condition p = 0.The remaining domain boundaries are treated as walls by imposing the no-slip condition.Gravity is acting in the negative y-direction with g = 9.81 m/s 2 .
The viscoelastic fluid is a single-mode Oldroyd-B fluid with the viscosity ratio β = 0.1, density ρ = 1000 kg/m 3 , and the total viscosity η t = 10 Pa • s.The surrounding Newtonian viscosity is air with viscosity µ air = 1.8205 × 10 −5 Pa • s and density ρ air = 1.204 kg/m 3 .
The Reynolds number for the flow is defined as Re = ρUD/(µ + η), which, for the chosen parameters, yields Re = 0.25.A Weissenberg number is defined as Wi = λU/D.The flow is simulated for the Oldroyd-B fluid with λ = 0.1 s, corresponding Wi = 10.As reference, a simulation is also performed with a Newtonian fluid, for which β = Wi = 0.In the viscoelastic jet buckling simulations that were performed in this work, the octree grid is adaptively refined with two levels around the viscoelastic fluid phase.As the simulation progresses, the grid is therefore updated and the cells are refined and coarsened, where needed.A grid dependence study is performed while using three adaptive grids of increasing resolution.The grids have the respective base cell sizes ∆x base = D/2, D/4, D/8, which correspond to the finest cell sizes ∆x min = D/8, D/16, D/32.The simulations are performed with a constant time step ∆t satisfying ∆tU/∆x min = 0.08.
The grid dependence is evaluated with respect to two simulation properties.Firstly, the free surface α = 0.5 is compared for the three adaptive grids.The first normal stress difference N 1 = τ yy − τ xx along the jet at x = 5D is also compared.Because the buckling itself is triggered by numerical round-off errors, the exact conditions at which the buckling starts in the simulations may vary between grids.Therefore, the grid study comparison is made at a simulation time before buckling is initiated.
In Figure 16, the free surface is compared for the three grids at the simulation times t = 130 ms and t = 140 ms.The results obtained with the three grids are in relatively good agreement, particularly those that are obtained with the two finer grids.In Figure 17, the first normal stress differences along the jet are compared at the same simulation times.These results indicate even more strongly that the two finer grids produce stresses that are very close, while a slight discrepancy for the coarsest grid is observed.The remainder of the results reported in this section have been obtained while using ∆x min = D/16.In Figures 18 and 19, a series of snapshots from the simulation of the Newtonian and the viscoelastic jets are shown, respectively.For the Newtonian case, only a very slight tendency to buckle is observed as a slight asymmetry in the fluid on the rigid surface.However, for the viscoelastic fluid, the buckling phenomenon is very apparent.When the viscoelastic jet impacts the rigid surface, the material first builds upwards until, inevitably, the buckling is initiated.The results demonstrate the strong influence of elasticity to the fluid jet buckling phenomenon.Furthermore, it confirms the capability of the proposed method to capture the viscoelastic effects and predict the phenomenon.

Conclusions
In this work, a new Lagrangian-Eulerian method for the simulation of viscoelastic free surface flow has been proposed.The method was developed from a previously proposed method, in which the fluid momentum and continuity equations were solved on a stationary Eulerian octree grid, while the viscoelastic constitutive equation was solved along the trajectories of Lagrangian nodes which were convected by the flow.The main improvements in this work was the introduction of a backwards-tracking procedure for solving the viscoelastic constitutive equation, as well an extension to the simulation of viscoelastic free surface flow with the volume of fluid method.
The backwards-tracking procedure allowed for the storage of the viscoelastic stresses in a structured arrangement.Consequently, the unstructured interpolation of stresses as well as the addition and deletion of Lagrangian nodes, which was required in the previous method, was eliminated in favor of a significant reduction of the computational cost and increased robustness.Furthermore, no additional stabilization method, such as both sides diffusion or reformulation of the constitutive equation, was employed for any of the flows studied in this work.
The proposed method was validated with analytic solutions for a planar Poiseuille flow of Oldroyd-B fluids.The transient, as well steady flow, solutions were in excellent agreement with analytic solutions and the steady solution converged to the analytic solution with second order accuracy.Furthermore, the proposed method was evaluated for two types of viscoelastic free surface flow.The viscoelastic die swell effect was simulated with the Oldroyd-B model and the predicted swell ratios were in good agreement with the numerical data from the literature.Planar jet buckling of a highly viscoelastic fluid was also simulated, demonstrating the capability of the method for an additional flow case as well as the use of adaptive grid refinements.
The computational performance was not subject to detailed study in this work.However, as per the construction of the new method, the computational cost for the viscoelastic stress calculation was significantly reduced when compared to the previous Lagrangian-Eulerian method.Furthermore, it has previously been shown that the Lagrangian formulation of the constitutive equation is highly suitable for parallel calculations, including for GPU-acceleration.This also applies for the new method.
In conclusion, the new method is an important step in the development of efficient and useful tools for simulating industrial processes involving viscoelastic fluids.In future research, the proposed Lagrangian-Eulerian method will be used for the numerical investigation of viscoelastic free surface flows with moving objects, including for viscoelastic adhesive application and parts assembly, as well as for additive manufacturing processes.where C 1 is an integration constant.Due to the symmetry condition ∂u ∂y | y=0 = 0, which by insertion of y = 0 into (A7) leads to C 1 = 0. Inserting C 1 = 0 to (A7) and integrating in y again yields where C 2 is an integration constant.Using the no-slip condition u y=h = 0 gives

Figure 2 .
Schematic description of the steps that are involved in the backwards-tracking algorithm.

Figure 3 .
Basis for interpolating properties stored at the fluid grid to Lagrangian nodes in areas with uniform grid spacing (a) and near refinements (b).The volume of fluid (VOF) method is used to model the presence of a viscoelastic fluid phase and a Newtonian phase.A color function ξ, is then defined, such that ξ = 1 In the viscoelastic phase 0 In the Newtonian phase .

1 Figure 5 .
Figure 5. Centerline velocity in channel flow simulated using different time steps.Xue et al. presented an analytic solution of the transient solution for the flow considered [46].In Figure 6, the obtained centerline velocities for the β = 1/9, 1/18, 1/27 are shown for Wi = 0.1 and Wi = 1 and compared to the analytic solution.The results demonstrate the strong influence of the relationship between the solvent viscosity and polymeric viscosity on the transient flow.The numerical results overlap the analytic solution, showing

Figure 7 .
Figure 7. Computed errors with respect to the analytic solution for velocity (a) and viscoelastic normal stress (b).

Figure 10 .
Figure 10.Simulated profiles for S R = 2.5 across the channel at x = 5h, obtained with grid M1, compared to the analytic solution.

Figure 11 .
Figure 11.Free surface position in the expansion zone of the die swell geometry, simulated for grids M1, M2, and M3 for S R = 2.5.

28 Figure 13 .
Figure 13.Snapshots of die swell simulation with S R = 2.5, interface between the viscoelastic phase (green) and Newtonian phase (white) visualized by α = 0.5.

Figure 14 .
Figure 14.Simulated swell ratio compared to data from the literature.

Table 1 .
Summary of grids that are used to study the grid dependence of the die swell simulations.Grid M1 used in the grid dependence study, as defined in Table 3U/(2h 2 ).The velocity can then be expressed in terms of the mean velocity asu(y) = 3U 2h 2 (h 2 − y 2 ) =which gives the final expressions for τ xy and τ xx .