Numerical Simulation of KBKZ Integral Constitutive Equations in Hierarchical Grids

: In this work, we present the implementation and veriﬁcation of HiGTree-HiGFlow solver (see for numerical simulation of the KBKZ integral constitutive equation. The numerical method proposed herein is a ﬁnite difference technique using tree-based grids. The advantage of using hierarchical grids is that they allow us to achieve great accuracy in local mesh reﬁnements. A moving least squares (MLS) interpolation technique is used to adapt the discretization stencil near the inter-faces between grid elements of different sizes. The momentum and mass conservation equations are solved by an implicit method and the Chorin projection method is used for decoupling the velocity and pressure. The Finger tensor is calculated using the deformation ﬁelds method and a three-node quadrature formula is used to derive an expression for the integral tensor. The results of velocity and stress ﬁelds in channel and contraction-ﬂow problems obtained in our simulations show good agreement with numerical and experimental results found in the literature.


Introduction
Over the years, several software programs have been developed to solve problems involving complex viscoelastic fluid flows. Due to a lack of generality, some challenges can arise when trying to solve problems with specific characteristics. In most works that develop numerical methods for simulating viscoelastic flows, the constitutive equations are approximated by differential equations, such as the Oldroyd-B [1][2][3], Upper-Convected-Maxwell (UCM) [4,5], Phan-Thien-Tanner [6,7], eXtended Pom-Pom [8,9] models, among others. However, advances in computational resources have motivated researchers to consider more sophisticated rheological models that are expressed in integral form instead of differential equations. In this sense, integral models allow a better approximation of the behavior of viscoelastic fluids. However, they require a greater computational effort and this is because, at each moment of the simulation, it is necessary to store and access the history of the entire deformation of the fluid (since it previously started to be deformed). Among the integral models that we found in the literature, the constitutive equation KBKZ-PSM has been considered by many researchers who study numerical methods for this kind of fluid. A detailed discussion of the importance of the KBKZ-PSM integral constitutive model and the development of numerical techniques to approximate integral models can be found in the works of Tanner [10] and Mitsoulis [11]. The vast majority of problems using the KBKZ-PSM integral model involve confined flows, such as channel-flows [12,13] and flows in abrupt contractions [3,14,15]. Flow problems possessing free surface(s) have also been considered by some researchers. More interesting flow problems that involve transient free surface(s) and integral models are the filament stretching [16] and a numerical study of the die swell phenomenon [17][18][19][20][21].
On the other hand, a numerical solution of partial differential equations in general grids has been questioned by many researchers in recent decades. Many schemes try to combine efficiency and simplicity with the flexibility of unstructured mesh networks. A major advantage of using such meshes is the ability to refine the mesh locally, improving accuracy in specific regions without dramatically increasing the number of unknowns. Among all possible ways of discretizing the spatial domain (simplified meshes, curvilinear meshes, among others), hierarchical meshes based on Cartesian trees are a common choice. They allow the development of finite difference methods, without the hassle of mapping and transforming distorted elements or dealing with general and complicated stencils, as in non-Cartesian grids. Since flows are generally computed on facets aligned with the Cartesian axis, numerical schemes are generally simpler to derive. However, these facets are generally shared by different numbers of elements on each side, which is the main challenge in implementing numerical methods. Different techniques to deal with this problem have been developed in the literature, most of them restricted to quadtree meshes (in 2D) or octree (3D) meshes, which are special cases of hierarchical grids represented by data structures quadtree/octree. Despite this restriction, these tree-based data structures are generally good enough and still an adequate choice for adaptive grids and moving borders [22]. Thus, we intend to implement in the present work the transient KBKZ-PSM model through a method of finite differences in hierarchical meshes that employ interpolations using the moving least squares (MLS) method [23]. The developed numerical method is verified by using mesh refinement in channel flow and we show results from the simulation of the 4:1 contraction problem using a KBKZ fluid. Our results are compared to the ones obtained using the OpenFOAM system [24], which uses finite volumes in the discretization of Navier-Stokes equations. We used the OpenFOAM v2006 version to implement the equations with the finite volume method.

Governing Equations
The governing equations for transient, isothermal and incompressible flows are the mass conservation and the equation of motion, which, in dimensionless form, can be written as follows (for details, see Tomé et al. [21]): Using the EVSS transformation [25], the extra-stress tensor τ is written as where S is a non-Newtonian tensor, v is the velocity field, p is the kinematic pressure and t is the time. In these equations, F represents the external forces, ǫ is a stability parameter (as shown in Araújo et al. [26]), Re = ρ 0 UL η 0 is the Reynolds number, η 0 is the zeroshear-rate viscosity, ρ 0 is the fluid density and U and L are the velocity and length scales, respectively. In this work, the rheological model that defines the behavior of fluid flow is the KBKZ-PSM [11] integral constitutive equation, which is shown below: where B t ′ (t) is the Finger tensor and M is the memory function, which adopts the following form: and H is the Papanastasiou-Scriven-Macosko [25] damping function, which is calculated using the following equation: The parameters λ k , a k , k are relaxation times, relaxation moduli and the number of relaxation modes, respectively. The quantities are the first and second invariants of B t ′ (t), respectively. The parameters a k , λ k , α, β are obtained from a curve fitting to the rheological properties of the fluid.
is the average relaxation time and the zero-shear-rate viscosity is written as η 0 = ∑ a k λ k .

Numerical Method
In this section, we present the methodology used in this work, the HiGTree/HiGFlow and the OpenFoam systems.

HiGTree/HiGFlow System
The HiGFlow system is a C language software, developed at ICMC-USP, which brings together a series of methods for the numerical simulation of flow of single-phase and multiphase fluids, using the finite difference technique. This system is being developed in a modular way, allowing new techniques and methods to be easily tested and added to the system. One feature is that the user chooses the dimension and the modules to be used in the program (such as single-phase, Newtonian, generalized Newtonian, viscoelastic) at compile time. In the same way, the user specifies the numerical techniques to be used in the input files: projection method, numerical scheme for the convective term, model of the constitutive equation for viscoelastic flows, in addition to the various parameters for simulation. In this work, all tests were performed in two dimensions (2D), and the following numerical techniques were chosen: an implicit Euler method to compute the velocity, the CUBISTA scheme to discretize the convective terms and an explicit Euler method for the convection of the Finger tensor.
On the other hand, the HiGTree system is responsible for creating the data structure, domains, linear and non-linear system solvers, as well as carrying out the interpolations schemes. Parallelization strategies are also implemented through the PETSc library (Portable, Extensible Toolkit for Scientific Computation), which contains a set of functions implementing the best-known methods for representing matrices, vectors and data in parallel, solution of linear systems with pre-conditioning, solution of linear and non-linear systems, ordinary differential equations, etc.

Hierarchical Grids
Equations (1) and (2) are approximated using finite differences in hierarchical Cartesian meshes. An illustrative representation of the mesh is given in Figure 1a and its structure of dependencies is illustrated in Figure 1b. In this data structure, each cell can be partitioned into distinct geometric shapes. Such generalization imposes difficulties in the numerical approximation in finite differences. For example, considering Figure 2, suppose that we are interested in approximating the second derivative in y direction centered on U c . Using second-order finite differences, we have: We can notice that, for this case, U b does not match known values in the mesh (recalling that the components of the velocity field are computed in the facet centers), but it can be calculated by interpolation using values of neighboring cells as shown in the following equation: The number of neighbors V b is defined according to the imposed precision. The weights w b k = w k (x) are calculated using the moving least squares (MLS) method [23].

Calculation of v(x, t n+1 ) and p(x, t n+1 )
Upon discretizing Equation (2) in time using, for instance, a first-order explicit discretization, the idea of the incremental projection method is to use the newest previous pressure field, which yields an explicitly computed velocity field v * that is not divergencefree, through the solution of the following equation: The corrected velocity field can be computed from the decomposition itself: where ϕ = −δt(p n+1 − p n ), which is obtained by solving the Poisson equation: with n.∇ϕ = 0 on the boundaries ∂Ω. Equation (10) can be easily derived by obtaining the divergence of Equation (9) (for more details, see [23,27]).

Calculation of the Extra-Stress Tensor τ(x, t n+1 )
We follow the methodologies described in Tomé et al. [21] and Araújo et al. [26] to calculate the extra-stress tensor τ(x, t n+1 ). The constitutive Equation (3) can be written as follows: The parameter s c (s c is a time interval) depends on the relaxation parameter λ re f . This methodology is called s-approach and is described in more detail in Hulsen et al. [28].
where N is a fixed number. Then, the integral equation can be written as: where and, therefore, the first integral becomes: which can be solved without any further issues. Regarding the integrals within the summation operator in Equation (12), we use the method of undetermined coefficients (with a second-order quadrate formula) for their calculation (for details, see Tomé et al. [21]). In the following sections, we describe the method used to compute the tensor B t ′ (t n+1 ) (t n+1 ) and how the points t ′ j (t n+1 ) are calculated.
One of the key issues of the deformation fields method is how the integration nodes , because such distribution can affect the accuracy of the results when solving complex flows. In Araújo et al. [26], the authors presented one discretization using a function that allowed them to determine the distribution of the time-integration points, which showed excellent results in some of the specific flow cases studied (such as extensional flows). However, care must be taken if we plan to generalize these results to more complex flows. In this work, we decide to use the more generic methodology presented in Tomé et al. [21]. We consider time-dependent flows so that the integration nodes are calculated using a geometric progression at time t n+1 as follows: One of the difficulties in the numerical simulation of viscoelastic flows using integral constitutive models is how to calculate accurately the strain history. In finite elements, this can be accomplished by a particle-tracking method based on the velocity field (see [12]), but, here, a different approach is taken. We follow the ideas of the deformation fields method [28] in which the Finger tensor is obtained by solving an is calculated using the Euler method and the highorder upwind scheme CUBISTA [29] is used to discretize the convective terms. We point out that the Finger tensor B t ′ (t) (x, t n+1 ) is calculated at the past times t ′ (t). The updated Finger tensor B t ′ (t n+1 ) (x, t n+1 ) is evaluated using a second-order interpolation method that is discussed in detail by Tomé et al. [3].

OpenFOAM System
All numerical experiments carried out in the present work will be compared with the results obtained using the OpenFOAM solver for integral models implemented by Araujo et al. [26]. The meshes were adapted in order to have simulations with similar conditions (and as close as possible) to the HiGFlow meshes. For the simulation of the contraction problem, for instance, the mesh shown in Figure 3 was used, where five regions with different refinements in the x direction can be observed. Notice that the upstream and downstream regions of the contraction geometry have volumes with exactly the same dimensions used in the HiGFlow simulations. On the other hand, a regular mesh was used for the channel-flow case. It is worth noting that the simulations were performed using the PISO method and half of the computational domain, considering the flow symmetry and the lower computational cost.
The coupling between stress and velocity was performed using the Improved Both Sides Diffusion (iBSD) [30] method, which adds a diffusive term on both sides of the momentum equation. For the solution of the linear systems resulting from the discretization of the velocity, the Bi-CGSTAB (BiConjugate Gradient Stabilized) method [31] was used with DILU (Simplified Diagonal-based Incomplete LU preconditioner) preconditioner, and, for the pressure, the conjugated preconditioned gradients (PCG) method was used with DIC (Simplified Diagonal-based Incomplete Cholesky) preconditioner.
In OpenFOAM, it is possible to choose the methods of discretization for some terms of an equation-for instance, diffusive or convective terms. Regarding this work, the numerical schemes used are described in Table 1.
where t c = min{t, s max }, N is the number of integration points and ξ is a parameter that depends on the value of s 1 (for more details, see [26]). All the simulations were performed using s 1 = ∆t and N = 51. The Finger tensors B(t n , t n − s k ) are convected according to Equation (14). We use a Euler explict scheme to obtain the fields B(t n+1 , t n − s k ). These fields are then interpolated, allowing us to calculate the fields B(t n+1 , t n+1 − s k ).

Results
In this section, we present a verification of the methodology described in Section 3.1. Initially, the methodology is applied to the channel-flow problem. Using several meshes (uniform and non-uniform), the HiGFlow system showed good agreement with the solution obtained with the OpenFOAM system. Results of meshes' orders and errors are also shown. Lastly, the numerical simulation of contraction flows is presented. The results are compared with solutions of the OpenFOAM system [26], Freeflow system [3], Mitsoulis [14] and Quinzani [15].

Mesh Independence in Channel-Flow
The numerical method described in Section 3 was applied to simulate the flow of a KBKZ fluid in a 2D planar channel (see Figure 4) of length 10L and height L, where L = 0.006 m. At the channel entrance, a dimensionless parabolic velocity profile given by u(y) = 4y(1 − y) was used. The scaling parameters were the centerline velocity, U = 0.1 ms −1 , and the fluid simulated was FLUID S1, whose parameters are described in Table 2. In this flow, we had Re = 0.34, De = 1, ǫ = 0.1, s c = 0.1s and the number of deformation fields was N = 100.   In order to verify the mesh convergence of the results, the flow was simulated using several meshes (see Tables 3 and 4 and Figure 5).  Figure 5 shows the non-uniform meshes, where we can see the structure of the mesh. In Figure 5a, the mesh MRI (with two levels of refinement) is depicted, and in Figure 5b, we can see the mesh MRV I (with three levels of refinement). The u-profiles are illustrated in Figure 6, where the mesh convergence can be seen. We adopted mesh MV as a reference mesh (black line) and the solutions of the refined meshes (full symbols) and uniform meshes (empty symbols) are shown. In this figure, we also show the OpenFoam system profile using the mesh MV. We saw good agreement between the solutions obtained in both systems.
Our results for the tensor components τ xx and τ yy are presented in Figure 7, where we can also see good agreement between the numerical solutions of the HiGFlow and the OpenFOAM systems. To verify the convergence, we show the errors (L 1 , L 2 , L ∞ ) and the orders in Table 5. The errors are calculated using the following equations:

Numerical Simulation of 4:1 Abrupt Planar Contraction Problem
In this section, we show the simulations for the 4:1 abrupt planar contraction flow. This problem is interesting because, for instance, the flow near the contraction is a complex mixture of shear and elongation, and secondary fluid motions might exist, even in the Newtonian limit (see [15]). For this reason, contraction flows have been extensively studied previously in the literature (see [3,14,15,26]). Figure 8 shows the domain representation, where we adopted a dimensionless parabolic inlet velocity profile u(y) = 3 8 1 4 (2 − y)(2 + y). The scaling parameter L = 0.0064 m is the height of the small channel, and we used N = 50 deformation fields, ǫ = 0.1 and the time interval s c = 0.1 s. In Table 6, we report the scaling parameters of average velocityū used in all simulations and the dimensionless parameters Re = ρLU η 0 , De = λU L , De(γ) (see [15]) and the characteristic shear rateγ.  The mesh M 1 used in these simulations is shown in Figure 9. We used three levels of refinement, with the most refined part near the contraction region. We also used one uniform mesh M 2 for De = 0.94 in order to check the convergence solutions in two meshes. In M 2 , d x = d y = 0.03125 m, and in M 1 , we use small values of dx = dy = 0.03125 m as well as larger values, dx = dy = 0.125 m. In Figure 10, we illustrate the centerline axial profile velocity u x (y) solutions for different values of the Deborah number De = 0.41, 0.94, 1.45 and 2.07 for the HiGFlow (green lines) and OpenFOAM (black lines) systems using mesh M1. The experimental results of Quinzani et al. [15] and numerical results of Mitsoulis [14] and Tomé et al. [3] (for De = 1.45) are also shown for comparison purposes. For the case with De = 0.94, we show two solutions using our methodology in mesh M 1 (non-uniform mesh) and mesh M 2 , where we can see good agreement between the solutions. For this reason, we adopted the M 1 mesh to simulate the other cases with different values of number De. For the profile u x (y), the HiGFlow system showed good agreement with the OpenFOAM solutions in the regions before and after the contraction. Near to contraction (see Figure 10 Figure 10. Centerline axial velocity profiles obtained using HiGFlow (green lines) and OpenFOAM (black lines). Experimental [15] and numerical results [3,14] found in the literature are also shown for comparison purposes.
In Figure 11, we show the numerical solution using HiGFlow (green lines) and Open-FOAM (black lines) systems for the tensor components τ xx and τ yy with the same flow parameter values reported in Table 6 using M 1 and M 2 for the case with De = 0.94. The methodologies presented in this work are different. OpenFOAM uses the finite volume method while HiGFlow approximates the equations using finite differences. Therefore, the solutions obtained will not be equal but should be comparable. Outside the contraction region, the OpenFOAM and HiGFlow solutions are very similar for all values of De. However, in the region close to the contraction, the solutions obtained using OpenFOAM showed a higher peak (x ≈ 0) but this is mostly seen in the cases with the highest values of De. In Figure 12a, we show the comparison between the first normal stress difference values N1 = τ xx − τ yy for two different cases of Deborah number, De = 0.41 and De = 1.45. For better visualization, the other two values of De (see Table 6) are illustrated in Figure 12b, where we can see that the values of HiGFlow (green lines) have good agreement with experimental data (orange triangles) reported by Quinzani [15], while the solution using OpenFOAM was similar to the solution presented by Mitsoulis [14].
Quinzani [15] Mitsoulis [14] Tomé [3] HF Quinzani [15] Mitsoulis [14] HF In Figure 13, we compare the streamlines obtained using OpenFOAM (a) and HiGFlow (b) with a fixed value of De = 2.07. We can see that there is vortex formation in both cases and that the solutions are relatively comparable.

Discussion
The work aims to present the numerical simulation of the KBKZ integral constitutive equations for incompressible and transient complex flows. We used a new solver HiGTree/HiGFlow lately developed by Souza et al. [23]. In this solver, we have implemented the methodology described in Section 3.1 to simulate viscoelastic flows modeled by integral constitutive equations. Initially, the numerical technique was verified by refined mesh in channel flows. Using the FLUID S1 (see Table 2), we performed nine simulations using non-refined and refined meshes. For comparison purposes, a mesh was chosen and the simulation using the OpenFOAM solver [26] was performed. In these simulations, we can see that, although the methodologies used in HiGFlow and OpenFOAM are quite different (the first uses finite differences and the second uses finite volume), we obtained very similar results in both systems. We also verified that the errors decrease with the mesh refinement and that the order of convergence of the velocity was around two, as expected.
A classic problem in the simulation of integral viscoelastic flows is known as 4:1 abrupt contraction and, thus, the literature for this problem is extensive. We chose to check our methodology for the four values of De presented in Mitsoulis [14]. In addition to the comparison with the results of this author, we performed the simulations using the OpenFOAM solver [26] and also compared our results with the experimental ones from Quinzani [15] and with the numerical results from FreeFlow [3]. We know that, in the contraction region, there are singularities and numerical techniques that might exhibit different behaviors. Although the values obtained by us in this work differ somewhat from the values obtained by Mitsoulis [14] or OpenFOAM [26], they were quite comparable to the experimental values of Quinzani [15]. Thus, we verified that the methodology presented here is capable of simulating complex flows in transient fluid regimes governed by integral constitutive models using the rapid technique of finite differences in hierarchical meshes with local refinement.
The computational efficiency of the models has been previously studied [3,21,26]. However, it is worth mentioning that the methods used in the present work used fewer integration points (deformation fields) compared to the points used in the early work of Hulsen et al. [28] and are similar to those used more recently by Hulsen and Anderson [32]. This improvement is due mainly to the distinct methodologies adopted to obtain the integration points, which allow efficient distribution of the elapsed time. Although the results presented here are two-dimensional, our methodology is still able to simulate flows in higher dimensions (the user just has to specify the dimension (2, 3 or higher) in the input data file). Our future work will be to present simulations in three dimensions for classic problems. In three dimensions, for example, integral models are still computationally expensive, as there is a need to store and connect a fixed N number of fields to each cell. Therefore, we will also work on ways to improve or modify the integral calculation-for instance, as was done in the recent work of Hulsen [32].