Journal of Computational Physics 302 (2015) 653–673 Contents lists available at ScienceDirect Journal of Computational Physics www.elsevier.com/locate/jcp A numerical study of the Kernel-conformation transformation for transient viscoelastic fluid flows F.P. Martins a, C.M. Oishi a,∗, A.M. Afonso b, M.A. Alves b a Departamento de Matemática e Computação, Faculdade de Ciências e Tecnologia, Universidade Estadual Paulista “Júlio de Mesquita Filho”, 19060-900 Presidente Prudente, SP, Brazil b Departamento de Engenharia Química, CEFT, Faculdade de Engenharia da Universidade do Porto, Porto, Portugal a r t i c l e i n f o a b s t r a c t Article history: Received 25 November 2014 Received in revised form 26 August 2015 Accepted 27 August 2015 Available online 2 September 2015 Keywords: High Weissenberg number problem Kernel conformation tensor Complex flows Viscoelastic fluids Die-swell phenomenon This work presents a numerical application of a generic conformation tensor transforma- tion for simulating highly elastic flows of non-Newtonian fluids typically observed in com- putational rheology. In the Kernel-conformation framework [14], the conformation tensor constitutive law for a viscoelastic fluid is transformed introducing a generic tensor trans- formation function. The numerical stability of the application of the Kernel-conformation for highly elastic flows is ultimately related with the specific kernel function used in the matrix transformation, but also to the existence of singularities introduced either by flow geometry or by the characteristics of the constitutive equation. In this work, we imple- ment this methodology in a free-surface Marker-And-Cell discretization methodology im- plemented in a finite differences method. The main contributions of this work are two fold: on one hand, we demonstrate the accuracy of this Kernel-conformation formulation using a finite differences method and free surfaces; on the other hand, we assess the numerical efficiency of specific kernel functions at high-Weissenberg number flows. The numerical study considers different viscoelastic fluid flow problems, including the Poiseuille flow in a channel, the lid-driven cavity flow and the die-swell free surface flow. The numerical re- sults demonstrate the adequacy of this methodology for high Weissenberg number flows using the Oldroyd-B model. © 2015 Elsevier Inc. All rights reserved. 1. Introduction Viscoelastic fluid flows have been successfully simulated using different numerical approaches over the last decades [1]. However, one obstacle usually found in the simulation of complex viscoelastic flows is the lack of convergence above a given (usually modest) Weissenberg number (Wi), and this limitation became known as the High Weissenberg Number Problem (HWNP). This numerical problem is generally related to the fail (or breakdown) in the numerical calculation of high Weissenberg number flows for differential type viscoelastic models (such as UCM, Oldroyd-B, FENE-P) for different computational methodologies (finite differences method – FDM, finite element method – FEM, finite volume method – FVM, etc.). The HWNP can thus be related to a deficiency of the discretization methods, because the exact form of the constitutive equation can admit regular solutions at the continuous level. As pointed out by Fattal and Kupferman [2,3], the HWNP is a numerical instability caused by the failure to balance the exponential growth of the stresses (due to deformation) with their * Corresponding author. E-mail addresses: fpmartins@fct.unesp.br (F.P. Martins), oishi@fct.unesp.br (C.M. Oishi), aafonso@fe.up.pt (A.M. Afonso), mmalves@fe.up.pt (M.A. Alves). http://dx.doi.org/10.1016/j.jcp.2015.08.038 0021-9991/© 2015 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2015.08.038 http://www.ScienceDirect.com/ http://www.elsevier.com/locate/jcp mailto:fpmartins@fct.unesp.br mailto:oishi@fct.unesp.br mailto:aafonso@fe.up.pt mailto:mmalves@fe.up.pt http://dx.doi.org/10.1016/j.jcp.2015.08.038 http://crossmark.crossref.org/dialog/?doi=10.1016/j.jcp.2015.08.038&domain=pdf 654 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 convection, which can be alleviated by means of a linearization of the conformation tensor fields using a tensor logarithm transformation. One key contribution for solving the HWNP in Computational Rheology was presented by Fattal and Kupferman [2,3]. These authors proposed a reformulation of the constitutive laws which describe the non-Newtonian behavior, using a matrix transformation, known as log-conformation, which improves the stability of numerical schemes. The use of the logarithm variable applied to the positive-definite conformation tensor A (formulation for this tensor will be discussed hereafter) drastically changed the view on the numerical simulation for HWNP. Soon after these seminal works, Fattal et al. [4] reported simulations using the log-conformation representation in the context of FEM for solving the flow of viscoelastic fluids past a confined cylinder at high Wi, a classical benchmark flow problem in computational rheology. The matrix transformation approach showed a superior performance and remedied the numerical instabilities of the standard formulation based on the extra-stress tensor, and the authors reported new benchmark results at significantly higher Wi than previous investigations. Since then, many computational techniques based upon this matrix transformation have been published, invariably showing good stability properties for solving challenging problems in viscoelastic flows (e.g. [5–12]). By taking the advantage that the conformation tensor A is a symmetric positive definite tensor, Balci et al. [13] proposed a simple and interesting square root matrix transformation. Briefly, the idea is that the conformation tensor possesses a unique symmetric square root tensor that satisfies a closed form evolution equation which also guarantees that A remains positive definite. Moreover, the authors report that their formulation is easily implemented and does not require the calculation of eigenvectors and eigenvalues at every time-step as happens in the classical log-conformation formulation, thus reducing the expected increase of computational times. Numerical results for a Stokes-viscoelastic system showed accuracy and stability improvements of the square-root transformation in a simple test case. Recently, in order to introduce a generic transformation for the conformation tensor evolution equation, Afonso et al. [14] described the kernel-conformation transformation, which uses an eigendecomposition to represent the conformation tensor A [15]. An interesting study involving numerical methods for HWNP compared different stabilization techniques but has not analyzed the kernel-conformation transformation [16]. More recently, in the context of finite element method, Kwon [17] applied the tensor-logarithmic formulation to investigate elastic instabilities in several benchmark problems. The goal of the present work is to demonstrate the accuracy of the kernel formulation using a finite differences method in free surface flows and assess the numerical efficiency and stability of specific kernel functions at high-Weissenberg number flows. The remainder of this paper starts with the set of modeling equations including the illustration of some kernel functions. Then, we briefly describe the numerical method and proceed to the presentation of the numerical results of the two-dimensional Poiseuille flow, the lid-driven cavity and die-swell free surface flow of an Oldroyd-B fluid. The paper ends with a brief summary of the results and conclusions. 2. Modeling equations Incompressible and isothermal flows are governed by the equation of motion and the mass conservation equation which can be written in dimensionless form as ∂u ∂t + (u · ∇)u = −∇p + β Re ∇2u + ∇ · τ + 1 Fr2 g, (1) ∇ · u = 0, (2) where t is the time, u is the velocity vector field, p is the pressure and g is the gravity field. In Eq. (1), τ is the non- Newtonian part of the extra-stress tensor which is defined by an appropriate constitutive equation. The dimensionless parameters Re = ρLU μ and Fr = U/ √ gL are the Reynolds and Froude numbers, respectively, where L and U are appropriate length and velocity scales, g is the magnitude of the gravity field g, ρ is the fluid density and μ is the total shear viscosity of the fluid. In this work we simulate problems involving the Oldroyd-B model. The solvent viscosity ratio β is defined as the ratio be- tween solvent and total shear viscosities, β = μS/(μS + μP ) = μS/μ. The traditional form to represent the non-Newtonian extra-stress tensor τ is given by the following hyperbolic-type partial differential equation, τ + Wi ∇ τ = 2 1 − β Re D (3) where D is the rate of deformation tensor, D = 1 2 [ (∇u) + (∇u)T ] , (4) and ∇ τ is the upper-convected time derivative of τ , defined as ∇ τ = Dτ Dt − τ · ∇u − ∇uT · τ (5) The Weissenberg number shown in Eq. (3) is defined as Wi = λU/L where λ is the relaxation time of the fluid. F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 655 An alternative form to describe viscoelastic models uses the conformation tensor, A. This tensor is symmetric and positive definite (SPD, see details in [18]) which is an important mathematical property for the construction of matrix transforma- tions and/or decompositions. In general the evolution equation for A can be written as ∂A ∂t + (u · ∇) A − [ A∇u + ∇uTA ] = 1 Wi F(A)H(A), (6) or as used by Fattal and Kupferman [2,3], ∂A ∂t + (u · ∇) A − �A + A� − 2BA = F (A) Wi H (A) , (7) where F(A) and H(A) are scalar and tensor functions, respectively, that depend on the specific constitutive model selected. Note that the conformation tensor can be diagonalized into: A = O�OT , (8) where O is an orthogonal matrix and � is a diagonal matrix. Fattal and Kupferman [2,3] reformulated the velocity gradient and its transpose into a symmetric traceless pure extension component, the tensor B that commutes with A, and an anti-symmetric pure rotation component, tensor �. This local decomposition rule for the traceless second-order transpose of velocity gradient tensor was obtained by setting: M = ∇uT = � + B + NA−1 M̃ = OT ∇uT O = �̃ + B̃ + Ñ�−1, (9) with the tensors ̃B, �̃ and Ñ representing the symmetric and the two anti-symmetric matrices, respectively, used originally by Fattal and Kupferman [2,3]. It can be represented for a two dimensional case as, OT ∇uT O = [ m̃11 m̃12 m̃21 m̃22 ] ⇒ ⎧⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ B̃ = B̃ ii = [ m̃11 0 0 m̃22 ] �̃ = [ 0 � −� 0 ] Ñ = [ 0 n −n 0 ] (10) with � = λ2m̃12+λ1m̃21 λ2−λ1 and n = m̃12+m̃21 λ−1 2 −λ−1 1 . Note that the term NA−1 vanishes in the back substitution of equation (9) onto equation (6), leading to the final form of equation (7) [2,3]. The relation between extra stress tensor τ and A is given by τ = ξ(A − I), (11) where ξ is a scalar function. For the Oldroyd-B model, equation (11) is written as τ = 1 − β Re Wi (A − I). (12) 3. Kernel conformation tensor transformation A generic Kernel-conformation tensor transformation for differential constitutive models was proposed by Afonso et al. [14]. In such framework, several matrix kernel-transformation families were applied to the conformation tensor evolution equation (6). This methodology is generic, and by itself does not guarantee a well-posed numerical solution at high Wi. This will depend on the selected kernel function (as also on the inverse function characteristics), and in particular on its ability to smooth the sharp conformation tensor field zones to better balance numerically the exponential growth and convection of the conformation tensor. Naturally, it is also beneficial when the selected kernel function preserves the positive definiteness of the conformation tensor, as happens for example with the logarithmic and square root transformations. The kernel-conformation constitutive law is obtained by introducing the following kernel-conformation tensor transforma- tion, k (A) = Ok (�) OT , (13) where k() represents a continuous, invertible and differentiable matrix function. The evolution equation for the kernel-conformation tensor, in its eigendecomposed formulation, can be expressed by [14] Dk (�) Dt = 2̃B�J + F Wi H (�) J (14) which is non-zero for the diagonal components of the tensor: 656 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 Dk (λi) Dt = 2B̃iiλi J ii + F Wi H (λi) J ii, (15) no summation over ii is implicit in Eq. (15). The gradient matrix, J, is a diagonal matrix of the form, J = diag ( ∂k (λ1) ∂λ1 ; ∂k (λ2) ∂λ2 ; ∂k (λ3) ∂λ3 ) (16) The evolution equation for k(A) can be expressed in tensorial notation as [14] Dk(A) Dt = �k(A) − k(A)� + 2B+ 1 Wi H (17) where B and H are symmetric tensors constructed by the orthogonalization of the diagonal tensors DB and DH , respectively. These tensors can be constructed as [14] B = ODBOT = OB̃�JOT H = ODHOT = FOH (�) JOT (18) 3.1. Examples of kernel-conformation transformations As already stated, the kernel-conformation transformation function can use any continuous, invertible and differentiable matrix transformation function. Here, for illustrative purposes, we will analyze some specific functions for k (A). For sim- plicity, the Oldroyd-B differential equation is used in this representation, for which F (A) = 1 and H (A) = I − A [14]. The rootk kernel-conformation transformation. In the rootk kernel-conformation transformation family, the kernel relation is defined as [14]: k (A) = A 1 k = O� 1 k OT (19) The inverse rootk kernel function is given by A = [k (A)]k (20) and the gradient matrix J is J = diag ⎛ ⎝∂λ 1 k 1 ∂λ1 ; ∂λ 1 k 2 ∂λ2 ; ∂λ 1 k 3 ∂λ3 ⎞ ⎠= � 1−k k k (21) The evolution equation for the rootk kernel-conformation, is expressed by Dk (A) Dt = �k (A) − k (A)� + 2 k Bk (A) + 1 kWi ( [k (A)]1−k − k (A) ) (22) For k = 1 (k (A) = A), the standard conformation formulation, given by equation (6) is recovered. For k = 2, equation (22) becomes Dk (A) Dt = �k (A) − k (A)� + Bk (A) + 1 2Wi ( [k (A)]−1 − k (A) ) . (23) This transformation is equivalent to the square-root-conformation proposed by Balci et al. [13] (note that k (A) = b = A 1 2 ): Db Dt = ab + b∇u + 1 2Wi (( bT )−1 − b ) (24) where a is an anti-symmetric tensor and the tensor b is symmetric, thus ( bT )−1 = b−1. The log kernel-conformation transformation. In the log kernel-conformation transformation, the functional kernel relation is based in the logarithm of base a of the conformation tensor [14]: k (A) = loga A = O loga (�)OT (25) The inverse kernel function is given by A = ak(A) (26) and the diagonal gradient matrix J is J = diag ( ∂ loga (λ1) ∂λ ; ∂ loga (λ2) ∂λ ; ∂ loga (λ3) ∂λ ) = �−1 ln(a) (27) 1 2 3 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 657 The evolution equation for the loga kernel-conformation tensor is given by Dk (A) Dt = �k (A) − k (A)� + 2 B ln(a) + 1 ln(a)Wi ( a−k(A) − I ) (28) For the natural logarithm (a = e; k(A) = ln(A)), the log-conformation transformation of Fattal and Kupferman [2,3] is recovered, D� Dt = �� − �� + 2B + 1 Wi (e−� − I) (29) 4. The kernel conformation tensor approach in a MAC-type method The kernel conformation tensor approach was originally implemented using a FVM framework [14]. In this section we describe this methodology for a Marker-And-Cell discretization implemented in the scope of the FDM framework. 4.1. Boundary conditions The boundary conditions used in this work for the application of the kernel conformation tensor in a FDM framework are: prescribed inflows, outflows, rigid walls and moving free surfaces. The velocity, prescribed at an inflow, is given by u = uinflow, (30) while at an outflow the homogeneous Neumann condition is employed. On the solid stationary walls, the no-slip condition is used (u = uwall). On the moving free surfaces, surface tension effects are neglected. In this case, we impose the following conditions [19] n · σ · nT = 0, (31) m · σ · nT = 0, (32) where σ is the total stress tensor, given by σ = −pI + β 2 Re D + τ . (33) In equations (31) and (32), n represents a unit vector normal to the free surface, and m is a unit vector tangent to the free surface. Boundary conditions on the extra-stress tensor τ are implemented as a function of the kernel conformation tensor transformation. 4.2. Staggered grid and cell classification The physical domain is discretized using an N × N staggered Cartesian grid with a set of equally spaced points (xi, y j) = ((i + 1/2)h, ( j + 1/2)h) with i, j = 0, . . . , N − 1 where h is the mesh size. The position of the variables in the staggered computational cell assumes that the velocity components are computed at the cell faces while all other variables are calculated at the cell center (i, j). The cells are classified in different types for applying the MAC-type discretization on the governing equations. Briefly, this classification includes: cells that do not contain fluid (abbreviated as [E]); cells that do not have any face in contact with empty cells [F]; cells that contain fluid and have one or more faces in contact with empty cells [S]; cells that belong to rigid boundaries [B]; cells that simulate fluid entrance into the domain [I]; and cells that define fluid exits of the domain [O]. Descriptions of cell labeling procedures, free surface discretization algorithms, and other implementation details were described in detail in previous works [20,21]. 4.3. Numerical method The implicit MAC-type method was used previously to simulate viscoelastic moving free surface flows [22,23]. In sum- mary, the Navier–Stokes equations (1)–(2) are solved using a projection method (or fractional-step algorithm) [24,25] while the non-Newtonian extra-stress tensor τ given by equation (3) is solved using a second-order Runge–Kutta scheme. In order to obtain a stable scheme, an implicit strategy was applied to the normal stress conditions (31) resulting in new equations for the pressure in free surface boundaries. Finally, the ordinary differential equation that models the advection of the particles in the free surface is solved by a second-order Runge–Kutta method. 658 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 In this work, this implicit MAC-type methodology is adapted to the kernel-conformation procedure and the new numer- ical method with corresponding modifications is presented below. The momentum and continuity equations (1)–(2) are discretized in time using a semi-implicit θ -method as u(n+1) δt − θ β Re ∇2u(n+1) = u(n) δt + (1 − θ) β Re ∇2u(n) − 3 2 ((u · ∇)u)(n) + 1 2 ((u · ∇)u)(n−1) − ∇p(n+1) + ∇ · τ (n+ 1 2 ) + 1 Fr2 g, (34) ∇ · u(n+1) = 0, (35) where the parameter θ is equal to 1 2 for the Crank–Nicolson scheme or 1 for the backward Euler method. Using the ideas of the projection method, an intermediate fluid velocity field ũ(n+1) is calculated from (34) resulting in the following equation: ũ(n+1) δt − θ β Re ∇2ũ(n+1) = u(n) δt + (1 − θ) β Re ∇2u(n) − 3 2 ((u · ∇)u)(n) + 1 2 ((u · ∇)u)(n−1) − ∇p(n) + ∇ · τ (n+ 1 2 ) + 1 Fr2 g(n), (36) where the boundary conditions for ̃u(n+1) are the same as those for u(n+1) and p(n) is an approximation to p(n+1) . The term ∇ · τ (n+ 1 2 ) in equation (36) is approximated by ∇ · τ (n+ 1 2 ) = 1 2 [ ∇ · τ (n) + ∇ · τ (n+1) ] . (37) Assuming that in this stage of the numerical scheme τ (n) is a known variable from the previous time step, we need to calculate the non-Newtonian tensor τ (n+1) for equation (37). Here, we introduce the kernel conformation tensor k(A), and the MAC algorithm is modified by the following steps: 1. Evolution equation for the kernel-conformation tensor. Firstly, we need to integrate Eq. (17) in time; for instance, in this work we applied a second order Runge–Kutta scheme. For this purpose, Eq. (17) is re-written in a compact form as ∂k(A) ∂t = f (u,k(A),�,B,H) , (38) where f (u,k(A),�,B,H) = −[(u · ∇)k(A)] + �k(A) − k(A)� + 2B+ 1 Wi H. (39) The next step involves the calculation of an intermediate kernel conformation tensor, k(A)(n+1) , by explicit forward Euler discretization, k(A)(n+1) − k(A)(n) δt = f ( un,k(A)(n),�(n),B(n),H(n) ) . (40) A detailed finite difference approximation for the terms in f () is presented in Appendix A. 2. Initial factorization. In order to obtain the values in the function f in Eq. (40), we begin the phase of matrix factor- izations calculating the eigenvalues and eigenvectors of A to construct the matrix O and the tensor �, following the ideas of Fattal and Kupferman [2,3]. After this, B is determined by the decomposition of the velocity gradient ∇uT . In summary, we have A(n) = Re Wi 1 − β τ (n) + I = (O�OT )(n), (41) M̃(n) = (OT ∇uT O)(n), (42) �(n) = (O�̃OT )(n), (43) B̃(n) = diag(M̃(n)), (44) B(n) = (OB̃OT )(n), (45) where �(n) is the diagonal matrix with the eigenvalues of A(n) , M̃(n) and �̃(n) are constructed according to Eqs. (9) and (10), respectively, and diag(M̃(n)) is a diagonal matrix obtained directly from the diagonal elements of M̃(n) . F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 659 3. Kernel factorization. By choosing a given kernel function, such as those presented in Section 3, and using the eigenval- ues and eigenvectors obtained in the previous step, we decompose: k(A)(n) = (Ok(�)OT )(n), (46) B (n) = (OB̃�JOT )(n) (47) H (n) = (FOH (�) JOT )(n) (48) where J(n) is given according to the kernel function selected (see Eq. (16)), and functions F and H(�) depend on the specific constitutive model used. 4. Solution of the evolution equation. We use Eqs. (43), (46), (47), (48) and the velocity field u to obtain the function f in time level (n). Once k(A)(n) and f are obtained, the intermediate kernel tensor k(A)(n+1) is calculated from Eq. (40). 5. Obtaining the non-Newtonian tensor. At this step, we decompose k(A)(n+1) = (Ok(�)OT )(n+1), (49) and according to kernel-conformation transformation (see Section 3), we apply the inverse kernel functions to calculate the conformation tensor A(n+1) . Finally, the non-Newtonian extra-stress tensor τ (n+1) is obtained from Eq. (12). Once τ (n+1) is obtained by the introduction of the kernel-conformation formulation previously described, we can solve Eq. (36), following the procedure of the MAC methodology. According to the fractional-step algorithm for solving the incompressible Navier–Stokes equations, in general, the velocity field ũ(n+1) , which was calculated from Eq. (36), does not satisfy the incompressibility constraint, i.e., ∇ · ũ(n+1) �= 0. Thus, in order to obtain the corrected velocity field, we use the Helmholtz–Hodge decomposition (cf. [26]) which states that every smooth vector field can be decomposed as a sum of a gradient of a potential function and a divergence-free vector field, u(n+1) = ũ(n+1) − ∇ψ(n+1). (50) Taking the divergence of equation (50) and imposing mass conservation on u(n+1) , we obtain the following equation ∇2ψ(n+1) = ∇ · ũ(n+1). (51) The boundary conditions for this Poisson equation are ψ(n+1) = 0 on outflows while the homogeneous Neumann condition is used for fixed boundaries and inflows. At free surfaces, depending on the direction of the normal vector, we apply the following equation ψ(n+1) δt − 2β Re [(∂2ψ(n+1) ∂ y2 ) n2 x + (∂2ψ(n+1) ∂x2 ) n2 y − 2 (∂2ψ(n+1) ∂x∂ y ) nxny ] − (1 − θ) β Re ∇2ψ(n+1) = 2β Re [ − (∂ ṽ(n+1) ∂ y ) n2 x − (∂ ũ(n+1) ∂x ) n2 y + (∂ ũ(n+1) ∂ y + ∂ ṽ(n+1) ∂x ) nxny ] + [(τ xx)(n+1)n2 x + 2(τ xy)(n+1)nxny + (τ yy)(n+1)n2 y] − p(n). (52) Equation (52) is derived from an implicit approach applied on the solution of normal stress boundary condition (31). This implicit discretization is imposed at free surfaces for obtaining a stable method for viscoelastic free surface flows. More details about this strategy can be found in [22] for two dimensional cases, and in [27] for three-dimensional cases. Following the steps of the projection method, the corrected velocity is calculated from equation (50) while the corrected pressure field is obtained from: p(n+1) = p(n) + ψ(n+1) δt − (1 − θ) β Re ∇2ψ(n+1). (53) The corrected velocity field u(n+1) and the intermediate kernel tensor k(A)(n+1) are employed to calculate the new time-level kernel tensor k(A)(n+1) from the following equation: k(A)(n+1) − k(A)(n) δt = 1 2 [ f ( u(n),k(A)(n),�(n),B(n),H(n) ) + f ( u(n+1),k(A)(n+1),�(n+1),B(n+1),H(n+1) )] . (54) In this stage of the algorithm, we reproduce steps (1)–(4), calculating, from Eq. (41) to (48), the tensors used in the function f (n+1) of Eq. (54). Thus, from the obtained values of k(A)(n+1) , the new time-level conformation tensor A(n+1) is calculated from the inverse kernel functions (see step (5)). Finally, the new time-level extra stress tensor τ (n+1) is determined from Eq. (12). 660 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 For free surface flows, the last step of the algorithm is the advection of the free surface interface. Each particle is convected by the velocity field, from their position x(n) at t = tn to the position xn+1 at t = tn+1, in agreement with dx dt = u. (55) This equation is solved using the second-order RK21 scheme as described in [22]. 5. Numerical results In this section, we present the numerical results for several kernel conformation transformations using three benchmark viscoelastic fluid flows: Poiseuille flow in a channel; lid-driven cavity test; extrudate swell phenomenon. 5.1. Two-dimensional Poiseuille flow In this section, we assess the accuracy of the kernel conformation transformations for solving the time dependent Poiseuille flow between parallel plates. We consider a planar channel of width L and length 10L and impose a parabolic velocity profile at the channel entrance: u(y) = 4y(1 − y), v = 0 . (56) The analytic solutions for the non-zero components of the non-Newtonian extra-stress tensor τ under fully-developed flow are given by τ xx(y) = 2Wi Re (1 − β) ( ∂u ∂ y )2 , τ xy(y) = 1 Re (1 − β) ( ∂u ∂ y ) . (57) The initial conditions are zero-field values for all flow variables. The following data are considered in all simulations: Re = 0.1 and β = 0.5. The characteristic velocity U is taken as the maximum fully-developed velocity at the center of the channel. In addition, for these computations, we used θ = 1 2 for the construction of the numerical method (see Section 4.3). A comparison between analytical and numerical solutions is given at the position x = 5L (middle of the channel) and the relative error is computed using the l2-norm: E(·) = ∥∥∥(·)Mi − (·)Exact ∥∥∥ 2 / ∥∥∥(·)Exact ∥∥∥ 2 , (58) where (·) denotes the solution field and i = 1, 2 indicates the mesh number. To analyze the accuracy of the transformations we performed a mesh refinement study using two uniform meshes: mesh M1 (h = 0.05, 200 × 20 cells) and mesh M2 (h = 0.025, 400 × 40 cells). The final dimensionless time of these simulations is t = 300 which is sufficiently long to achieve steady-state at the monitoring location. Table 1 displays the calculated errors (using u, τ xx and τ xy components) using δt = 10−4 for Wi = 1 and 2.5. For Wi = 1, the kernel conformation transformations offer similar accuracy as the standard extra-stress formulation. In addition, all methods analyzed are stable for this low Wi. For Wi = 2.5, the standard formulation is unstable, resulting in numerical breakdown of the simulation. In a similar manner, the linear transformation is also unstable in the two meshes used in this test problem. As can be seen in Table 1, the stability and accuracy of the results of the root kernel transformation depend on the value of k used. In particular, looking for refined mesh results, we have unstable results for k = −1 and k = −10. On the other hand, at Wi = 2.5, the results of simulations using k = −2, 2, 10 and the Ln function confirm that some kernel trans- formations can stabilize the numerical computations in HWNP. Particularly, analyzing the errors presented in Table 1, we can conclude that, at least for this benchmark problem, the Ln function is more accurate than the root functions analyzed, k = −2, 2, 10. In order to obtain more information about the kernel-conformation transformation at high Wi, the minimum value of the determinant of the conformation tensor, detmin(A), was calculated to assess the time history of the positive definiteness of the conformation tensor. According to Hulsen [18], for the Oldroyd-B model the following condition should hold: det(A) ≥ 1. Results at Wi = 1 and Wi = 2.5 for meshes M1 and M2 are presented in Table 2, which also indicates the time of the numerical divergence in unstable computations. From Table 2, we can verify that for Wi = 1 nearly all methods satisfy the condition detmin(A) ≈ 1 confirming the numerical stability of the results shown in Table 1 for transient calculations. However, for Wi = 2.5, the standard method, linear function and the root kernel transformation for k = −1 and k = −10 did not preserve the positive definiteness of the conformation tensor resulting in numerical divergence. Finally, it is interesting to note that for the more stable kernel transformations (root function with k = 2, 10 and Ln function) the value of detmin(A) converges to 1, as expected for stabilization strategies in HWNP using the Oldroyd-B model. F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 661 Table 1 Numerical results for the two-dimensional Poiseuille flow. Wi Method E(uM1) E(τ M1 xx ) E(τ M1 xy ) E(uM2) E(τ M2 xx ) E(τ M2 xy ) 1.0 Standard 3.733E−03 7.423E−03 3.715E−03 9.441E−04 1.866E−03 9.324E−04 Linear 3.731E−03 7.423E−03 3.715E−03 9.349E−04 1.866E−03 9.324E−04 root (k = −10) 3.733E−03 7.423E−03 3.716E−03 9.330E−04 1.866E−03 9.322E−04 root (k = −2) 3.732E−03 7.421E−03 3.715E−03 9.376E−04 1.865E−03 9.321E−04 root (k = −1) 3.734E−03 7.417E−03 3.714E−03 9.433E−04 1.865E−03 9.315E−04 root (k = 2) 3.729E−03 7.424E−03 3.716E−03 9.270E−04 1.866E−03 9.325E−04 root (k = 10) 3.732E−03 7.424E−03 3.716E−03 9.385E−04 1.866E−03 9.325E−04 Ln 3.734E−03 7.424E−03 3.716E−03 9.321E−04 1.866E−03 9.324E−04 2.5 Standard – – – – – – Linear – – – – – – root (k = −10) 3.690E−03 7.139E−03 4.011E−03 – – – root (k = −2) 4.966E−03 2.134E−02 3.789E−02 1.014E−03 1.339E−02 2.3269E−02 root (k = −1) 4.411E−02 7.902E−01 4.327E−01 – – – root (k = 2) 6.118E−03 1.066E−01 6.020E−02 2.675E−03 7.238E−02 3.707E−02 root (k = 10) 6.110E−03 1.240E−01 7.739E−02 3.108E−03 9.401E−02 5.252E−02 Ln 3.694E−03 7.020E−03 3.919E−03 9.520E−04 1.941E−03 1.351E−03 Table 2 Minimum value of the determinant of the conformation tensor (detmin(A)) and the time of numerical divergence (denoted as tBD). Wi Method tBD(M1) detmin(A(M1)) tBD(M2) detmin(A(M2)) 1.0 Standard 1.000 1.000 Linear 1.000 1.000 root (k = −10) 0.988 0.983 root (k = −2) 0.903 0.903 root (k = −1) 0.795 0.788 root (k = 2) 1.000 1.000 root (k = 10) 1.000 1.000 Ln 1.000 1.000 2.5 Standard 58.34 – 91.60 – Linear 58.35 – 91.60 – root (k = −10) 0.998 0.53 – root (k = −2) 0.899 0.898 root (k = −1) 0.998 18.98 – root (k = 2) 1.000 1.000 root (k = 10) 1.000 1.000 Ln 1.000 1.000 5.2. Two-dimensional lid-driven cavity In this section, we assess the applicability of the kernel functions in the simulation of the lid-driven cavity flow. This benchmark problem has been widely used to assess the stability and accuracy of numerical methods for Newtonian flows. There are also some works that investigate numerically the lid driven cavity flow of viscoelastic fluids (e.g. [3,9,28–35]) typically at low Re flow conditions. We present the results obtained for this benchmark problem using the Oldroyd-B model, and compare with published data from the literature. Following previous studies, we consider a regularized parabolic profile for the top lid: u(x, t) = 8[1 + tanh(8t − 4)]x2(1 − x)2. This regularized velocity profile is used in order to vanish the velocity gradient on the corners and avoid a singularity of the stress field [3,36]. On the remaining cavity walls, we also impose no-slip boundary conditions but the walls are stationary. In all simulations we used the Oldroyd-B model with a solvent viscosity ratio β = 0.5. We consider creeping flow conditions, and set the Reynolds number to Re = 0.01. The Weissenberg number effect is analyzed, and is here defined as Wi = λU/L, where the characteristic velocity U is the maximum wall velocity (at (x, y) = (1/2, 1)) and the characteristic length scale is the cavity width, which is equal to the cavity height. To study the convergence of the numerical method with mesh refinement, the cavity flow was simulated using two meshes: M1 (h = 1 128 , 128 × 128 cells) and M2 (h = 1 256 , 256 × 256 cells). For all figures in this section, the profiles of u-velocity and the non-Newtonian τxx component are plotted along the vertical line x = 0.5 while the v-velocity component is reported at the horizontal line y = 0.75. The origin of the cartesian coordinate system is placed at the lower left corner of the square cavity, and the moving wall is located at y = 1. We start presenting the results for a low Weissenberg number, Wi = 0.5, in order to compare our results with the predictions of the log-conformation finite element method of [9] and the standard extra stress-tensor finite difference formulation of [22]. Figs. 1, 2 and 3 show that the results of Linear, Ln and square-root kernel transformations are very similar with those obtained in [9]. In addition, we present the profiles of the polymeric normal stress component τxx for 662 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 Fig. 1. Numerical simulation of lid-driven cavity flow at t = 20 using the Oldroyd-B model with Re = 0.01, Wi = 0.5 and β = 0.5: (a) u-velocity at (0.5, y); (b) v-velocity at (x, 0.75). Results obtained with the linear transformation and with the standard formulation (Std). the kernel transformations in Fig. 4. The results presented in Figs. 1–4 show that for all kernel functions the simulations are converging with mesh refinement. The time step used in the simulations was fixed, δt = 0.001, both for meshes M1 and M2, and the results are reported at t = 20, which for Wi = 0.5 already corresponds to steady state flow conditions. In order to assess the numerical stability of the kernel conformation transformations, the lid-driven cavity was also simulated for higher Wi flow conditions. In particular, for Wi ≥ 1 and for the Oldroyd-B model, Fattal and Kupferman [3] and Pan et al. [9] have also provided numerical results for this benchmark problem. Figs. 5 and 6 show a mesh refinement study and comparison with literature data for the velocity components at Wi = 1 and Wi = 2, respectively. We present the results only for Ln and square root kernel functions since numerical instabilities are observed for Wi ≥ 1 in the Linear function and standard formulations. Figs. 5 and 6 show a good agreement among our results and those from the literature, showing the potential of the kernel formulations to simulate this benchmark problem at high Wi. In particular, for Wi = 2, results obtained using the Ln kernel function are in good agreement with those reported by Fattal and Kupferman [3]. Results for τxx are shown in Fig. 7 illustrating again the numerical convergence of the solutions with mesh refinement. Finally, in order to further assess the stability of the kernel functions, we simulate this problem for Wi = 3. We compare in Fig. 8(a) and (b) the computed velocity profiles at t = 80 with those presented in [3], demonstrating good agreement among the solutions. The numerical convergence of the normal stress component τxx is depicted in Fig. 9(a). To provide additional results on this benchmark problem, we performed a quantitative comparison of the time evolution of the kinetic energy [3,9] obtained by the kernel functions and the data of Fattal and Kupferman [3]. The kinetic energy is computed as E = ∫ ‖u‖2dA/ ∫ dA [3], although more correctly this quantity represents the average of the square of the A A F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 663 Fig. 2. Numerical simulation of lid-driven cavity flow at t = 20 using the Oldroyd-B model with Re = 0.01, Wi = 0.5 and β = 0.5: (a) u-velocity at (0.5, y); (b) v-velocity at (x, 0.75). Results obtained with the Ln transformation and with the standard formulation (Std). velocity. In order to confirm the numerical convergence of the kernel transformations, we also present the time evolution of the kinetic energy for the square root function in a more refined mesh, defined as M3 (h = 1 320 , 320 × 320 cells). Fig. 9(b) shows that our results evolve to steady state showing a peak of the kinetic energy at E ≈ 0.018 for t ≈ 0.8, followed by a decrease towards an asymptotic value of E ≈ 0.0093 for Ln kernel function in mesh M2 and E ≈ 0.0090 for square root kernel function in mesh M3. In addition, our results did not present an oscillatory behavior of the kinetic energy as observed in the simulations in [3]. In our computations, oscillations were observed for δt = 0.001 resulting in unstable results; to avoid this problem we decreased the time step to δt = 0.0001 for Wi = 3 in order to maintain the stability of the kernel functions at high Wi. 5.3. Application to moving free surface flows: time-dependent extrudate swell problem The numerical simulation of the time-dependent extrudate swell problem for viscoelastic fluids has received substantial attention over the last decades (e.g. [37–41]). However, for high Weissenberg number flows (Wi > 1), numerical results of this phenomenon are very scarce, as a result of numerical instabilities that eventually lead to divergence of the numerical methods. To our best knowledge, the first application of a matrix transformation in the constitutive equation to enhance numerical stability in the simulation of the extrudate swell problem was published recently by Choi and Hulsen [42]. In that work, the authors applied the log-conformation representation in an extended finite element method using the upper-convected 664 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 Fig. 3. Numerical simulation of lid-driven cavity flow at t = 20 using the Oldroyd-B model with Re = 0.01, Wi = 0.5 and β = 0.5: (a) u-velocity at (0.5, y); (b) v-velocity at (x, 0.75). Results obtained with the root (k = 2) transformation and with the standard formulation (Std). Maxwell (UCM) model. The maximum value of the Weissenberg number achieved in the numerical simulations of [42] was Wi > 1.5. More recently, Tomé et al. [43] employed the log-conformation tensor approach in a finite difference method to simulate 3D extrudate swell of UCM fluids and obtained numerical results up to Wi = 2.5. Despite these advances, the understanding of this phenomenon for high Wi is still a challenging task for the computational rheology community. The purpose of this section is to present quantitative data of the time-dependent extrudate swell problem, for bench- marking purposes. We consider a 2D-channel with width L and length 4L and an outflow boundary located at a distance 6L from the channel exit. A domain size of 10L × 3L was employed. On the channel entrance, channel walls and outflow, the boundary conditions were the same as those employed in Section 5.1. The problem was solved using two uniform meshes: M1 (h = 0.05, 200 × 60 cells) and M2 (h = 0.025, 400 × 120 cells). The constant non-dimensional numbers used in all simulations were: Re = 0.1, 1 Fr = 0 (negligible gravitational effects) and β = 0.5. The characteristic velocity scale U is the maximum velocity at the center of the channel. The kernel functions investigated in this problem were the root function with k = 2 and the Ln function, which showed good performance in the previous test-cases. In order to compare our results, we also have included the simulation for low Wi computed with the standard extra-stress formulation. The numerical simulation of the time-dependent extrudate-swell problem is a typical example where the free surface is formed by moving marked particles. In this case, maintaining the global mass conservation is important to ensure the accuracy of the simulations, and in this section we assess the mass conservation. F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 665 Fig. 4. Numerical simulation of lid-driven cavity flow at t = 20 using Oldroyd-B model with Re = 0.01, Wi = 0.5 and β = 0.5. Non-Newtonian normal stress τxx is plotted along the vertical line, x = 0.5. Comparison between different kernel transformations in meshes M1 and M2. Firstly, note that the injected fluid at time t (V exact) can be calculated as V exact = Q t, (59) where Q is the non-dimensional flow rate per unit depth injected at the channel (Q = ∫ 1 0 u(y)dy). In order to calculate the numerical injected fluid (V num2D ) we use the following equation V num2D = 1 2 ∣∣∣∣∣ N−1∑ k=1 (Xk+1 − Xk)(Yk+1 + Yk) ∣∣∣∣∣ , (60) derived from computing the surface integral of the piecewise linear interpolation of the marker particles over the fluid do- main. In Eq. (60), N is the total number of particles defining the free surface profile at time t , with (Xk, Yk) the coordinates of the k-th particle. The relative error between the numerical and exact injected fluid is calculated as E V (t) = V exact − V num2D V exact × 100%. (61) Tables 3 and 4 display the results obtained on meshes M1 and M2 using time-steps δt = 10−4 and δt = 5 × 10−5, respectively. We present the relative error E V (t) and the minimum value of the determinant of A obtained in simulations of the standard and kernel formulations. For the standard formulation, we present the results only for Wi = 0.5 since for Wi ≥ 1 the numerical scheme becomes unstable at later times. Results are presented at t = 15 which represents an intermediate stage of the dynamics of the problem, as illustrated in Fig. 10. Tables 3 and 4 show that, for a fixed time-step, the error in mass conservation usually decreases with mesh refinement. However, the reduction of the time-step has a small influence in the mass conservation showing that for the time-step δt = 10−4 the time integration error is already small enough and the errors arise mostly from spatial discretization. At Wi = 0.5 the three methodologies show oscillating values of E V , thus the asymptotic behavior of the error with mesh refinement was not yet achieved, at least in mesh M1. For Wi ≥ 1, the computed error E V is similar for all formulations, and decreases significantly with mesh refinement. Regarding the positive definiteness of the conformation tensor, Tables 3 and 4 clearly show that the use of both the Ln function and the root function with k = 2 maintains the conformation tensor symmetric and positive definite while evolving on time, and the expected behavior of detmin(A) → 1 is preserved in all cases. An important flow parameter in the numerical simulation of the extrudate-swell problem is the swelling ratio, Sr , given by Sr = Lmax L , (62) where L is the width of channel and Lmax is the maximum width of the jet. For the higher values of Wi used in our study, numerical results available in the literature are scarce. The swell ratio was calculated for all numerical formulations at t = 45, significantly after the jet reaches the outflow, as shown in Fig. 10. 666 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 Fig. 5. Numerical simulation of the lid-driven cavity flow at t = 40 using the Oldroyd-B model with Re = 0.01, Wi = 1 and β = 0.5: (a) u-velocity at (0.5, y) and (b) v-velocity at (x, 0.75). Table 3 Numerical results of the minimum value of the determinant of A and numerical error (61) for the two-dimensional time-dependent extrudate swell problem using time-step δt = 10−4. Results are captured at t = 15. Wi Method detmin(A(M1)) E V (M1) detmin(A(M2)) E V (M2) Wi = 0.5 std 1.000 8.380E–02 1.000 −9.145E–04 root (k = 2) 1.000 −2.839E–02 1.000 4.047E–02 Ln 1.000 −4.233E–02 1.000 3.534E–02 Wi = 1.25 root (k = 2) 1.000 −9.471E–02 1.000 −1.826E–02 Ln 1.000 −8.973E–02 1.000 −3.169E–02 Wi = 1.5 root (k = 2) 1.000 −1.137E–01 1.000 −1.356E–02 Ln 1.000 −1.153E–01 1.000 −2.468E–02 Wi = 1.75 root (k = 2) 1.000 −1.083E–01 1.000 −9.257E–03 Ln 1.000 −1.178E–01 1.000 −7.073E–03 Results for the minimum value of the determinant of A and the swell ratio using meshes M1 and M2 are presented in Table 5. It is clear that all formulations follow the same trend, with an increase of Sr when Wi increases. For low Weissenberg number flows (e.g. Wi = 0.5) we observe a reasonable agreement of the swell ratio calculated using the standard and kernel formulations. For Wi ≥ 1, the standard formulation is numerically unstable. F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 667 Fig. 6. Numerical simulation of the lid-driven cavity flow at t = 80 using the Oldroyd-B model with Re = 0.01, Wi = 2 and β = 0.5: (a) u-velocity at (0.5, y) and (b) v-velocity at (x, 0.75). Table 4 Numerical results of the minimum value of the determinant of A and numerical error (61) for the two-dimensional time-dependent extrudate swell problem using time-step δt = 5 × 10−5. Results are captured at t = 15. Wi Method detmin(A(M1)) E V (M1) detmin(A(M2)) E V (M2) Wi = 0.5 std 1.000 8.377E–02 1.000 −1.600E–03 root (k = 2) 1.000 −2.842E–02 1.000 4.076E–02 Ln 1.000 −4.276E–02 1.000 3.676E–02 Wi = 1.25 root (k = 2) 1.000 −9.669E–02 1.000 −1.810E–02 Ln 1.000 −8.852E–02 1.000 −3.196E–02 Wi = 1.5 root (k = 2) 1.000 −1.131E–01 1.000 −1.579E–02 Ln 1.000 −1.149E–01 1.000 −2.515E–02 Wi = 1.75 root (k = 2) 1.000 −1.081E–01 1.000 −8.910E–03 Ln 1.000 −1.188E–01 1.000 −8.053E–03 To provide further evidence concerning the stability of the kernel formulations, we have investigated the extrudate-swell problem for higher values of Wi. Particularly, in Table 5, one can observe that the positive definiteness of the conformation tensor of the square root and Ln functions at t = 45 maintains the same behavior observed previously at t = 15 (detmin(A) = 1, cf. Tables 3 and 4). In spite of the stable results obtained in the simulations for Wi > 1.75 using the kernel formulations, 668 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 Fig. 7. Numerical simulation of the lid-driven cavity flow using the Oldroyd-B model with Re = 0.01 and β = 0.5: (a) Wi = 1 and (b) Wi = 2. The non- Newtonian normal stress τxx is plotted along the vertical line x = 0.5. Table 5 Numerical results of the minimum value of the determinant of A and the swelling ratio for the two-dimensional time-dependent extrudate swell problem using time-step δt = 5 × 10−5. Results are captured at non-dimensional time t = 45. Wi Method detmin(A(M1)) Sr(M1) detmin(A(M2)) Sr(M2) Wi = 0.5 std 1.000 1.490 1.000 1.400 root (k = 2) 1.000 1.414 1.000 1.369 Ln 1.000 1.422 1.000 1.355 Wi = 1.25 root (k = 2) 1.000 1.875 1.000 1.890 Ln 1.000 1.796 1.000 1.835 Wi = 1.5 root (k = 2) 1.000 1.970 1.000 2.056 Ln 1.000 1.953 1.000 2.052 Wi = 1.75 root (k = 2) 1.000 2.126 1.000 2.273 Ln 1.000 2.342 1.000 2.394 we have not pursued simulation at higher Wi in order to avoid the need of using a longer computational domain and more refined meshes, to maintain good accuracy. F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 669 Fig. 8. Numerical simulation of the lid-driven cavity flow at t = 80 using the Oldroyd-B model with Re = 0.01, Wi = 3 and β = 0.5: (a) u-velocity at (0.5, y) and (b) v-velocity at (x, 0.75). 6. Conclusions This work presented a numerical implementation of the Kernel-conformation analytical framework in a free-surface Marker-And-Cell discretization, implemented in the scope of a finite differences method. The numerical cases included the two-dimensional Poiseuille flow, the lid-driven cavity and die-swell free surface flow. The numerical results demonstrate the ability of this methodology for simulating high Weissenberg numbers flows using the Oldroyd-B model with the appropriate kernel functions. For all problems studied in this work, the limiting Weissenberg number increases using appropriate kernel-conformation formulations (such as the Ln and the square root, k = 2), in com- parison with the standard formulation which usually diverged at modest values of Wi. Acknowledgements C.M. Oishi acknowledges the financial support of FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo), Grants 2009/15892-9 and 2013/07375-0, and CNPq (Conselho Nacional de Desenvolvimento Cientifico e Tecnologico), Grants 473589/2013-3 and 309514/2013-4. A.M. Afonso would also like to thank the Fundação para a Ciência e a Tecnologia (FCT), Portugal, for financial support through the scholarship SFRH/BPD/75436/2010. M.A. Alves acknowledges funding from the FCT, COMPETE and FEDER, through project PTDC/EME-MFE/114322/2009. 670 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 Fig. 9. Numerical simulation of the lid-driven cavity flow at t = 80 using the Oldroyd-B model with Re = 0.01, Wi = 3 and β = 0.5: (a) τxx at (0.5, y) and (b) time evolution of the kinetic energy. Fig. 10. Illustration of the configuration for the extrudate-swell problem using the Ln function for Wi = 1.5, Re = 0.1 and β = 0.5 at: (a) t = 15 and (b) t = 45. Appendix A. Finite difference approximations of evolution equation for kernel-conformation tensor Considering two-dimensional Cartesian coordinates, the finite difference equations of the evolution equation (17) for kernel-conformation tensor are given by F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 671 ∂kxx i, j ∂t = − [conv(ukxx)i, j + conv(vkxx)i, j ]+ 2 ( k xy i, j� xy i, j +B xx i, j ) + 1 Wi H xx i, j, (63) ∂k yy i, j ∂t = − [conv(ukyy)i, j + conv(vkyy)i, j ]+ 2 ( −k xy i, j� xy i, j +B yy i, j ) + 1 Wi H yy i, j , (64) ∂k xy i, j ∂t = − [conv(ukxy)i, j + conv(vkxy)i, j ]+ � xy i, j ( k yy i, j − k xx i, j ) + 2Bxy i, j + 1 Wi H xy i, j, (65) where the convective terms conv(˙)i, j are approximated using the CUBISTA high-resolution scheme [44]. Note that for each kernel function described in Section 3.1 we need to obtain the components of the symmetric tensors B and H, for instance: • rootk kernel-conformation B xx i, j = 1 k ( Bxx i, jk xx i, j + Bxy i, jk xy i, j ) (66) B yy i, j = 1 k ( Bxy i, jk xy i, j + B yy i, jk yy i, j ) (67) B xy i, j = 1 k ( Bxx i, jk xy i, j + Bxy i, jk yy i, j ) (68) H xx i, j = 1 k ( (k1−k)xx i, j − k xx i, j ) (69) H yy i, j = 1 k ( (k1−k) yy i, j − k yy i, j ) (70) H xy i, j = 1 k ( (k1−k) xy i, j − k xy i, j ) (71) • log kernel-conformation transformation B xx i, j = Bxx i, j ln (a) (72) B yy i, j = B yy i, j ln (a) (73) B xy i, j = Bxy i, j ln (a) (74) H xx i, j = 1 ln (a) [( a−k )xx i, j − 1 ] (75) H yy i, j = 1 ln (a) [( a−k )yy i, j − 1 ] (76) H xy i, j = 1 ln (a) ( a−k )xy i, j (77) Remark 1. In order to construct the components of the tensors � and B involved in Eqs. (63)–(65), ∇uT is approximated in the cell center (i, j): ( ∇uT ) i, j = ⎡ ⎢⎢⎣ ui+ 1 2 , j − ui− 1 2 , j δx ui, j+ 1 2 − ui, j− 1 2 δy vi+ 1 2 , j − vi− 1 2 , j δx vi, j+ 1 2 − vi, j− 1 2 δy ⎤ ⎥⎥⎦ (78) where the terms which are not defined in the center are computed using the following averages: ui, j+ 1 2 = 1 4 ( ui+ 1 2 , j + ui− 1 2 , j + ui+ 1 2 , j+1 + ui− 1 2 , j+1 ) , (79) vi+ 1 2 , j = 1 4 ( vi, j+ 1 2 + vi, j− 1 2 + vi+1, j+ 1 2 + vi+1, j− 1 2 ) . (80) 672 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 Remark 2. One important remark in the context of finite difference approximations is the fact that Eq. (36) involves the computation of the stress divergence ∇ · τ at cell edges (for instance, (i + 1 2 , j) and (i, j + 1 2 )). In the present work, these terms are approximated as: ( ∂τ xx ∂x + ∂τ xy ∂ y ) i+ 1 2 , j � τ xx i+1, j − τ xx i, j δx + τ xy i+ 1 2 , j+ 1 2 − τ xy i+ 1 2 , j− 1 2 δy (81) and ( ∂τ yy ∂ y + ∂τ xy ∂x ) i, j+ 1 2 � τ yy i, j+1 − τ yy i, j δy + τ xy i+ 1 2 , j+ 1 2 − τ xy i− 1 2 , j+ 1 2 δx (82) where τ xy i+ 1 2 , j+ 1 2 = 1 4 ( τ xy i+1, j+1 + τ xy i+1, j + τ xy i, j + τ xy i, j+1 ) , τ xy i+ 1 2 , j− 1 2 = 1 4 ( τ xy i+1, j−1 + τ xy i+1, j + τ xy i, j + τ xy i, j−1 ) , τ xy i− 1 2 , j+ 1 2 = 1 4 ( τ xy i−1, j+1 + τ xy i−1, j + τ xy i, j + τ xy i, j+1 ) . (83) References [1] R.G. Owens, T.N. Phillips, Computational Rheology, Imperial College Press, London, 2002. [2] R. Fattal, R. Kupferman, Constitutive laws for the matrix-logarithm of the conformation tensor, J. Non-Newton. Fluid Mech. 123 (2004) 281–285. [3] R. Fattal, R. Kupferman, Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation, J. Non-Newton. Fluid Mech. 126 (2005) 23–37. [4] M.A. Hulsen, R. Fattal, R. Kupferman, Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms, J. Non-Newton. Fluid Mech. 127 (2005) 27–39. [5] Y. Kwon, Finite element analysis of planar 4:1 contraction flow with the tensor-logarithmic formulation of differential constitutive equations, Korea– Aust. Rheol. J. 16 (2004) 183–191. [6] O.M. Coronado, D. Arora, M. Behr, M. Pasquali, A simple method for simulating general viscoelastic fluid flows with an alternate log-conformation formulation, J. Non-Newton. Fluid Mech. 147 (2007) 189–199. [7] R. Guénette, A. Fortin, A. Kane, J.F. Hétu, An adaptive remeshing strategy for viscoelastic fluid flow simulations, J. Non-Newton. Fluid Mech. 153 (2008) 34–45. [8] A.M. Afonso, P.J. Oliveira, F.T. Pinho, M.A. Alves, The log-conformation tensor approach in the finite-volume method framework, J. Non-Newton. Fluid Mech. 157 (2009) 55–65. [9] T.W. Pan, J. Hao, R. Glowinski, On the simulation of a time-dependent cavity flow of an Oldroyd-B fluid, Int. J. Numer. Methods Fluids 60 (2009) 791–808. [10] A. Jafari, N. Fiétier, M.O. Deville, A new extended matrix logarithm formulation for the simulation of viscoelastic fluids by spectral elements, Comput. Fluids 39 (2010) 1425–1438. [11] H. Damanik, J. Hron, A. Ouazzi, S. Turek, A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems, J. Non-Newton. Fluid Mech. 165 (2010) 1105–1113. [12] A.M. Afonso, P.J. Oliveira, F.T. Pinho, M.A. Alves, Dynamics of high-Deborah-number entry flows: a numerical study, J. Fluid Mech. 677 (2011) 272–304. [13] N. Balci, B. Thomases, M. Renardy, C.R. Doering, Symmetric factorization of the conformation tensor in viscoelastic fluid models, J. Non-Newton. Fluid Mech. 166 (2011) 546–553. [14] A.M. Afonso, F.T. Pinho, M.A. Alves, The kernel-conformation constitutive laws, J. Non-Newton. Fluid Mech. 167–168 (2012) 30–37. [15] T. Vaithianathan, L.R. Collins, Numerical approach to simulating turbulent flow of a viscoelastic polymer solution, J. Comput. Phys. 187 (2003) 1–21. [16] X. Chen, H. Marschall, M. Schafer, D. Bothe, A comparison of stabilization approaches for finite-volume simulation of viscoelastic fluid flow, Int. J. Comput. Fluid Dyn. 27 (2013) 229–250. [17] Y. Kwon, Numerical aspects in modeling high Deborah number flow and elastic instability, J. Comput. Phys. 265 (2014) 128–144. [18] M.A. Hulsen, Some properties and analytical expressions for plane flow of Leonov and Giesekus models, J. Non-Newton. Fluid Mech. 30 (1988) 85–92. [19] G.K. Batchelor, An Introduction of Fluid Dynamics, Cambridge University Press, Cambridge, 1967. [20] S. McKee, M.F. Tomé, V.G. Ferreira, J.A. Cuminato, A. Castelo, F.S. Sousa, N. Mangiavacchi, The MAC method, Comput. Fluids 37 (2008) 907–930. [21] M.F. Tomé, S. McKee, GENSMAC: a computational marker-and-cell method for free surface flows in general domains, J. Comput. Phys. 110 (1994) 171–186. [22] C.M. Oishi, F.P. Martins, M.F. Tomé, J.A. Cuminato, S. McKee, Numerical solution of the eXtended Pom–Pom model for viscoelastic free surface flows, J. Non-Newton. Fluid Mech. 166 (2011) 165–179. [23] C.M. Oishi, F.P. Martins, M.F. Tomé, M.A. Alves, Numerical simulation of drop impact and jet buckling problems using the eXtended Pom–Pom model, J. Non-Newton. Fluid Mech. 169–170 (2012) 91–103. [24] A.J. Chorin, Numerical solution of the Navier–Stokes, Math. Comput. 2 (1968) 745–762. [25] R. Témam, Sur l’approximation de la solution des équations de Navier–Stokes par la méthode de pas fractionnaires (ii), Arch. Ration. Mech. Anal. 33 (1969) 377–385. [26] J. Cantarella, D. DeTurck, H. Gluck, Vector calculus and the topology of domains in 3-space, Am. Math. Mon. 109 (2002) 409–442. [27] C.M. Oishi, M.F. Tomé, J.A. Cuminato, S. McKee, An implicit technique for solving 3D low Reynolds number moving free surface flows, J. Comput. Phys. 227 (2008) 7446–7468. [28] A.M. Grillet, B. Yang, B. Khomami, E.S.G. Shaqfeh, Modeling of viscoelastic lid driven cavity flow using finite element simulations, J. Non-Newton. Fluid Mech. 88 (1999) 99–131. [29] E. Mitsoulis, T. Zisis, Flow of Bingham plastics in a lid-driven square cavity, J. Non-Newton. Fluid Mech. 101 (2001) 173–180. http://refhub.elsevier.com/S0021-9991(15)00565-3/bib5068696C6C6970733A32303032s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib46617474616C3A32303034s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib46617474616C3A32303035s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib46617474616C3A32303035s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib48756C73656E3A32303035s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib48756C73656E3A32303035s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4B776F6E3A32303034s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4B776F6E3A32303034s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib436F726F6E61646F3A32303037s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib436F726F6E61646F3A32303037s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4775656E657474653A32303038s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4775656E657474653A32303038s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib41666F6E736F3A32303039s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib41666F6E736F3A32303039s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib50616E3A32303039s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib50616E3A32303039s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4A61666172693A32303130s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4A61666172693A32303130s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib44616D616E696B3A32303130s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib44616D616E696B3A32303130s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib41666F6E736F3A32303131s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib42616C63693A32303131s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib42616C63693A32303131s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib41666F6E736F3A32303132s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib566169746869616E617468616E3A32303033s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4368656E3A32303133s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4368656E3A32303133s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4B776F6E3A32303134s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib48756C73656E3A31393838s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4261746368656C6F723A31393637s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4D634B65653A32303038s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib546F6D653A31393934s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib546F6D653A31393934s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4F697368693A32303131s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4F697368693A32303131s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4F697368693A32303132s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4F697368693A32303132s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib43686F72696E3A31393638s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib54656D616D3A31393639s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib54656D616D3A31393639s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib43616E746172656C6C613A32303032s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4F697368693A3230303862s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4F697368693A3230303862s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4772696C6C65743A31393939s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4772696C6C65743A31393939s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib4D6974736F756C69733A32303031s1 F.P. Martins et al. / Journal of Computational Physics 302 (2015) 653–673 673 [30] D. Vola, L. Boscardin, J.C. Latché, Laminar unsteady flows of Bingham fluids: a numerical strategy and some benchmark results, J. Comput. Phys. 187 (2003) 441–456. [31] F. Belblidia, H. Matallah, B. Puangkird, M.F. Webster, Alternative subcell discretisations for viscoelastic flow: stress interpolation, J. Non-Newton. Fluid Mech. 146 (2007) 59–78. [32] K. Yapici, B. Karasozen, Y. Uludag, Finite volume simulation of viscoelastic laminar flow in a lid-driven cavity, J. Non-Newton. Fluid Mech. 164 (2009) 51–65. [33] J. Zhang, An augmented Lagrangian approach to Bingham fluid flows in a lid-driven square cavity with piecewise linear equal-order finite elements, Comput. Methods Appl. Mech. Eng. 199 (2010) 3051–3057. [34] G. Bohme, A. Muller, Stress analysis for a cavity flow of a memory fluid, Arch. Appl. Mech. 81 (2011) 1807–1826. [35] R. Comminal, J. Spangenberg, J.H. Hattel, Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log- conformation formulation, J. Non-Newton. Fluid Mech. 223 (2015) 37–61. [36] M. Renardy, Stress integration for the constitutive law of the upper convected Maxwell fluid near the corners in a driven cavity, J. Non-Newton. Fluid Mech. 112 (2003) 77–84. [37] M.J. Crochet, R. Keunings, Finite element analysis of die-swell of a highly elastic fluid, J. Non-Newton. Fluid Mech. 10 (1982) 339–356. [38] E. Brasseur, M.M. Fyrillas, G.C. Georgiou, M.J. Crochet, The time-dependent extrudate-swell problem of an Oldroyd-B fluid with slip along the wall, J. Rheol. 42 (1998) 549–566. [39] M. Tomé, N. Mangiavacchi, J.A. Cuminato, A. Castelo, S. McKee, A finite difference technique for simulating unsteady viscoelastic free surface flows, J. Non-Newton. Fluid Mech. 106 (2002) 61–106. [40] G.S. de Paulo, M.F. Tomé, S. McKee, A marker-and-cell approach to viscoelastic free surface flows using the PTT model, J. Non-Newton. Fluid Mech. 147 (2007) 149–174. [41] G. Russo, T.N. Phillips, Spectral element predictions of die-swell for Oldroyd-B fluids, Comput. Fluids 43 (2011) 107–118. [42] Y.J. Choi, M.A. Hulsen, Simulation of extrudate swell using an extended finite element method, Korea–Aust. Rheol. J. 23 (2011) 147–154. [43] M. Tomé, A. Castelo, A.M. Afonso, M.A. Alves, F.T. Pinho, Application of the log-conformation tensor to three-dimensional time-dependent free surface flows, J. Non-Newton. Fluid Mech. 175–176 (2012) 44–54. [44] M.A. Alves, P.J. Oliveira, F.T. Pinho, A convergent and universally bounded interpolation scheme for the treatment of advection, Int. J. Numer. Methods Fluids 41 (2003) 47–75. http://refhub.elsevier.com/S0021-9991(15)00565-3/bib566F6C613A32303033s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib566F6C613A32303033s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib42656C626C696469613A32303037s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib42656C626C696469613A32303037s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib5961706963693A32303039s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib5961706963693A32303039s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib5A616E673A32303130s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib5A616E673A32303130s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib426F686D653A32303131s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib436F6D6D696E616C3A32303135s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib436F6D6D696E616C3A32303135s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib52656E617264793A32303033s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib52656E617264793A32303033s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib63726F6368657431393832s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib42726173736575723A31393938s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib42726173736575723A31393938s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib546F6D653A32303032s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib546F6D653A32303032s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib5061756C6F3A32303037s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib5061756C6F3A32303037s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib527573736F3A32303131s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib43686F693A32303131s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib546F6D653A32303132s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib546F6D653A32303132s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib416C7665736375626973746132s1 http://refhub.elsevier.com/S0021-9991(15)00565-3/bib416C7665736375626973746132s1 A numerical study of the Kernel-conformation transformation for transient viscoelastic fluid flows 1 Introduction 2 Modeling equations 3 Kernel conformation tensor transformation 3.1 Examples of kernel-conformation transformations 4 The kernel conformation tensor approach in a MAC-type method 4.1 Boundary conditions 4.2 Staggered grid and cell classification 4.3 Numerical method 5 Numerical results 5.1 Two-dimensional Poiseuille flow 5.2 Two-dimensional lid-driven cavity 5.3 Application to moving free surface flows: time-dependent extrudate swell problem 6 Conclusions Acknowledgements Appendix A Finite difference approximations of evolution equation for kernel-conformation tensor References