**1. Introduction**

Many researchers are constantly working on improving numerical solution techniques for partial differential equations that govern the flow of Newtonian and non-Newtonian fluids. One of the major problems faced is the part that generates the geometry of the problem to be simulated.

Cartesian hierarchical grids, or tree-based grids are the most common choices for discretizing the spatial domain. This choice allows the implementation of the finite difference method, while avoiding working with more complicated stencils, which occurs for example in curved meshes. Thus, it becomes easier to process the flow properties when a refinement of the mesh is desired in a determined region of the domain, since in a Cartesian grid, the flows are calculated in facets parallel to the Cartesian axes, favoring the implementation of the numerical method [1]. In the literature, quadtree and octree are 2D and 3D meshes, respectively, generated to perform these problems. One can say that a hierarchical grid is a generalization of quadtree and octree. In this sense, the choice of hierarchical grids is convenient to address the problem of flows in complex geometries [2–6].

Regarding the mesh refinement, one of the difficulties is to calculate the flow property value on these interfaces. High-order interpolations are commonly used. Several improvements of the interpolation techniques have been studied [7–11], in order to optimize the number of cells used in calculations, since this influences the computational time and storage over simulations.

In this way, the HiG-Fow system makes interpolations using the method of moving least squares, adapting the stencil according to the interface between the fine and coarse

**Citation:** Castelo, A.; Afonso, A.M.; De Souza Bezerra, W. A Hierarchical Grid Solver for Simulation of Flows of Complex Fluids. *Polymers* **2021**, *13*, 3168. https://doi.org/10.3390/ polym13183168

Academic Editors: Roland G. Winkler and Carlos A. García-González

Received: 7 July 2021 Accepted: 10 September 2021 Published: 18 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

grids. Sousa et al. developed this methodology and compared it with non-graded methods by using the new system to simulate Newtonian flows [12].

Our interest is to use the HiG-Flow for the simulation of non-Newtonian flows. In this way, a code module for simulations of non-Newtonian flows was implemented, taking into account considerations shown in Section 3.

Depending on the temperature or mixture in liquid solvents, polymeric materials behave similar to viscoelastic fluids [13]. In this work we show that a new computer system is able to perform numerical simulations of viscoelastic fluid flows in two and three dimensions in channels with complex geometries. In one of the most common applications, polymers are used to construct electronic devices. Thus, the study of viscoelastic fluids is important due to applications in science and technology and the use of numerical simulators can be useful for support in important decisions in the engineering design of any device. In general, the behavior of viscoelastic fluid can be described using an appropriate constitutive model. So, in addition to using HiG-Flow in Newtonian flows, the system has implemented a module to solve the constitutive equations through the kernelconformation technique. [14]. Different constitutive models are implemented, among them the Oldroyd-B [15] and Phan-Thien–Tanner (PTT) [16], which we used as reference in this work. In the last 20 years several works involving the solution of these constitutive models have been published. In 1999, Dou and Phan-Thien [17] used the finite volume method to solve the flow in a channel of an Oldroyd-B fluid past a circular cylinder. Then, Alves et al. [18] showed the effect of a high-resolution scheme MINMOD [19] on an upperconvected Maxwell fluid solution, improving accuracy and increasing the convergence rate of the finite volume method and then they proposed a new high resolution scheme [20]. Later the article was published [21] with benchmark solutions for the flow of Oldroyd-B and PTT fluids in planar contractions. In the year 2005 Chinyoka et al. [22] studied the deformation of a circular drop of an Oldroyd-B fluid by applying the volume-of-fluid method for two-dimensional interfaces. Later, Tomé et al. [23] applied the finite difference method to simulate free surface flow of PTT fluid in three dimensional geometry. Then, Mompean et al. [24] investigated fluid flows using the Upper-Convected Maxwell (UCM) constitutive equation and an explicit algebraic model to develop an approximation that could be applied to the extrudate-swell problem. In 2012, Tom é et al. [25] applied the log-conformation technique to study three-dimensional viscoelastic flows for jet buckling analysis and later Oishi et al. [26] and Paulo et al. [27] continued studies in this same way.

In 2019, Tomé et al. [28] presented a solution method for the Giesekus model flow and proposed a new analytical solution for this problem. In 2019, Bezerra et al. [29] used HiG-Flow to perform the solution of electro-osmotic flow of a viscoelastic fluid, where they proposed an approximation for the vortices simulation in a nozzle. Shojaei et al. [30] investigated a generalized finite difference method using the weighted moving least squares procedure, in the same way of our proposed numerical solution. Corresponding with one of the proposals of this work, [31] used stabilization techniques in 2D and 3D viscoelastic fluid flows. In 2020, Guan et al. proposed a improved finite difference method and they checked its convergence. Recently, [32] presented a implementation and computational verification of KBKZ integral constitutive equations in hierarchical grids. More recently, ref. [33] performed a generalized finite differences method for flows in a dam.

The finite difference method was used in the discretization of equations. The HiG-Flow system was also implemented taking into account advances in the MAC-Marker and Cell method [34], allowing the implementation of several solution methods for the different terms of the equation of motion as well as the constitutive model solution. Convective terms in equations can be solved by high-accuracy methods. Moreover, we can say that the main novelty for the simulation of viscoelastic fluids is the kernel-conformation technique. The technique is already known, however, the differential is the manner it was implemented in, which allows the user to choose a numerical stabilizer easily—one just needs to enter the desired mathematical function, the derivative of this function and its inverse function. More details can be found in the Governing Equations section. Here, numerical stabilizers

were used for Oldroyd-B flow solution in a 2D cavity and for a PTT fluid in a complex 3D geometry.

In Section 2, we show the finite differences method of the approximation used. Then, the governing equations and the constitutive models are presented in Section 3, as well as the explanation of the kernel-conformation technique. In Section 4, we present the validation tests for a PTT fluid flow in a pipe and to an Oldroyd-B fluid flow in a 2D-liddriven cavity. Finally, we performed simulations of a PTT flow in a complex 3D geometry and the results are shown in Section 5.

#### **2. Finite Difference Approximation in Tree-Based Grids**

In the HiG-Flow code, the equations are solved using finite difference approach in hierarchical meshes. Figure 1a shows a representative type of mesh and the dependency structure (tree data structure) is presented in Figure 1b. In this approach, cells can be partitioned into different geometrical shapes. Such generalization leads to the difficult task of finding an accurate approximation to the different differential and integral operators.

**Figure 1.** (**a**) Example of hierarchical grid. (**b**) Tree data structure. (**c**) Finite difference method.

Looking at Figure 1c, a second-order approximation to *<sup>∂</sup>* <sup>2</sup>*U<sup>c</sup> ∂y* <sup>2</sup> can be given by (we assume the *y* axis is in the direction bottom → top):

$$\frac{\partial^2 \mathcal{U}\_c}{\partial y^2} \approx \frac{1}{\delta y} (\mathcal{U}\_l - 2\mathcal{U}\_c + \mathcal{U}\_b). \tag{1}$$

Note that *U<sup>b</sup>* is not known and must be obtained by interpolation (the same applies to *Ur*) using the following formula:

$$\mathcal{U}l\_b = \sum\_{k=1}^{V\_b} w\_k^b \mathcal{U}\_{k\prime} \tag{2}$$

where V*<sup>b</sup>* is the number of neighbor cells, which depends on the imposed accuracy of the method.

The weights *w b <sup>k</sup>* = *w b k* (**x**) are obtained through the moving least squares (MLS) method. In a set of *n* smooth interpolating functions that are linearly independent **Φ***i* : <sup>R</sup>*<sup>d</sup>* <sup>→</sup> <sup>R</sup> (*<sup>d</sup>* <sup>=</sup> 2, 3), we want to obtain the interpolated value *<sup>u</sup>* such that *<sup>U</sup><sup>b</sup>* <sup>=</sup> *<sup>U</sup>*(**x**) = ∑ *n k*=1 *ci***Φ***i*(**x**) = **c** *<sup>t</sup>***Φ**.

Given *<sup>m</sup>* points **<sup>x</sup>**1, **<sup>x</sup>**2, . . . , **<sup>x</sup>***<sup>m</sup>* <sup>∈</sup> <sup>R</sup>*<sup>d</sup>* with *<sup>m</sup>* <sup>&</sup>gt; *<sup>n</sup>* and *<sup>m</sup>* values *<sup>u</sup>*<sup>1</sup> <sup>=</sup> *<sup>u</sup>*(**x**1), *<sup>u</sup>*<sup>2</sup> <sup>=</sup> *<sup>u</sup>*(**x**2), . . ., *u<sup>m</sup>* = *u*(**x***m*), to interpolate *u* in **x** using MLS consists in minimizing the error *E*(**c**)

$$E(\mathbf{c}) = \|\mathcal{U} - \boldsymbol{u}\|\_2^2 = \sum\_{i=1}^m (\mathcal{U}(\mathbf{x}\_i) - \boldsymbol{u}\_i)^2 \frac{1}{\|\mathbf{x} - \mathbf{x}\_i\|\_2}. \tag{3}$$

Or,

$$E(\mathbf{c}) = \|\boldsymbol{W}\mathbf{P}\mathbf{c} - \boldsymbol{W}\mathbf{u}\|\_{\mathcal{D}'}^2 \tag{4}$$

where *<sup>W</sup>* <sup>=</sup> *<sup>W</sup>*(**x**) = <sup>n</sup> *<sup>δ</sup>ij*<sup>q</sup> <sup>1</sup> k**x**−**x***i*k<sup>2</sup> o <sup>∈</sup> <sup>R</sup>*m*×*m*, *<sup>P</sup>* <sup>=</sup> **Φ***i*(**x***j*) <sup>∈</sup> <sup>R</sup>*m*×*<sup>n</sup>* and **<sup>u</sup>** <sup>=</sup> (*u*1, *u*2, . . . , *um*).

The solution is given by

$$\mathbf{c}(\mathbf{x}) = (\mathcal{W}P)^{\dagger}\mathcal{W}\mathbf{u} \tag{5}$$

where (·) † is the Moore–Penrose pseudo-inverse. Decomposing (*WP*) into *QR* we have that

$$WP = Q\begin{bmatrix} \ R \\ \mathbf{0} \end{bmatrix} = \begin{bmatrix} \mathbf{Q}\_{\parallel} & \mathbf{Q}\_{\perp} \end{bmatrix} \begin{bmatrix} \ R \\ \mathbf{0} \end{bmatrix} \tag{6}$$

where *<sup>Q</sup>* <sup>∈</sup> <sup>R</sup>*m*×*m*, *<sup>Q</sup>*<sup>k</sup> <sup>∈</sup> <sup>R</sup>*m*×*<sup>n</sup>* , *<sup>Q</sup>*<sup>⊥</sup> <sup>∈</sup> <sup>R</sup>*m*×*m*−*<sup>n</sup>* , and *<sup>R</sup>* <sup>∈</sup> <sup>R</sup>*n*×*<sup>n</sup>* . This decomposition is then used to finally compute

$$\mathbf{c}(\mathbf{x}) = \mathbb{R}^{-1} \mathbb{Q}\_{\parallel}^{t} \mathcal{W} \mathbf{u} \; . \tag{7}$$

Then

$$\mathcal{U}I(\mathbf{x}) = \mathbf{c}^t \boldsymbol{\Phi} = \mathbf{u}^t \underbrace{\mathcal{W}Q\_{\parallel}\mathcal{R}^{-t}\boldsymbol{\Phi}}\_{\mathbf{w}} = \sum\_{k=1}^{m} w\_k u\_k \,. \tag{8}$$

that is **<sup>w</sup>** = **<sup>w</sup>**(**x**) = *WQ*k*<sup>R</sup>* <sup>−</sup>*t***Φ**.

The procedure to calculate **w**(**x**) must be performed for each approximation *U*(**x**). This is performed only once since we are using a static mesh.

#### **3. Governing Equations**

The flow is assumed to be isothermal, laminar and the fluids incompressible. The governing equations are those expressing conservation of mass

$$\nabla \cdot \mathbf{u} = \mathbf{0},\tag{9}$$

and conservation of momentum

$$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u} + \nabla \cdot \mathbf{S} + \frac{1}{Fr^2} \mathbf{g} + \mathbf{F} \tag{10}$$

$$\mathbf{T} = \frac{2(1-\beta)}{Re}\mathbf{D} + \mathbf{S} \,\tag{11}$$

where **u** is the velocity field, *t* is time, *p* is the pressure, *Re* is the Reynolds number, *Fr* is the Froude number, **g** is the gravity force and **F** is the surface tension force and source force. The symbol **D** = <sup>1</sup> 2 ∇**u** + (∇**u**) *T* is the rate of deformation tensor, **T** is the elastic stress. The amount of Newtonian solvent is controlled by the dimensionless solvent viscosity coefficient, *β* = *µS µ*0 , where *µ*<sup>0</sup> = *µ<sup>S</sup>* + *µ<sup>P</sup>* denotes the total shear viscosity. Several polymeric constitutive equations are implemented in the current version of the solver: the upper-convected Maxwell model, the Oldroyd-B model, the linear form of the Phan-Thien/Tanner (LPTT) model [35] and the Giesekus model [36]. For an isothermal flow, these five rheological equations of state can be written in a compact form as:

$$\frac{\partial \mathbf{T}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{T} - \left[\left(\nabla \mathbf{u}\right)^{T} \cdot \mathbf{T} + \mathbf{T} \cdot \nabla \mathbf{u}\right] = \frac{1}{De} \mathbf{M}(\mathbf{T}).\tag{12}$$

where **M**(**T**) is defined by the viscoelastic model

$$\mathbf{M}(\mathbf{T}) = \begin{cases} \frac{2(1-\beta)}{Re}\mathbf{D} - \mathbf{T} & \text{Oldroyd-B} \\ \frac{2(1-\beta)}{Re}\mathbf{D} - \mathbf{T} - \frac{a \operatorname{Re} De}{1-\beta}\mathbf{T} \cdot \mathbf{T} & \text{Gieseks}, \\ \frac{2(1-\beta)}{Re}\mathbf{D} - \left(1 + \frac{\varepsilon \operatorname{Re} De}{1-\beta}\operatorname{tr}(\mathbf{T})\right)\mathbf{T} - \xi^{\tau}De(\mathbf{T}\cdot\mathbf{D} + \mathbf{D}\cdot\mathbf{T}) & \text{LPTt}, \end{cases} \tag{13}$$

where *De* is the Deborah number. The stress coefficient function of the LPTT model depends on the trace of **T**, tr(**T**) and introduces the dimensionless parameter *e*, which is closely related to the steady-state elongational viscosity in extensional flows. The slip parameter, *ξ*, takes into account the non-affine motion between the polymer molecules and the continuum. The polymer strands embedded in the medium may slip with respect to the deformation of the macroscopic medium, thus each strand may transmit only a fraction of its tension to the surrounding continuum. When *ξ* = 0 there is no slip and the motion becomes affine. Parameter *ξ* is responsible for a non-zero second normal-stress difference in shear, leading to secondary flows in ducts having non-circular cross-sections, which is superimposed on the streamwise flow. In the non-linear term of the Giesekus model, *α* represents a dimensionless "mobility factor".

An alternative form to describe viscoelastic models is by using the conformation tensor, **A**. This tensor is Symmetric and Positive Definite (SPD), which is an important mathematical property for the construction of matrix transformations and/or decompositions. In general, the equation for **A** can be written as

$$\frac{\partial \mathbf{A}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{A} - \left[\mathbf{A}\nabla \mathbf{u} + \nabla \mathbf{u}^{\mathsf{T}} \mathbf{A}\right] = \frac{1}{D\epsilon} \mathcal{M}(\mathbf{A}),\tag{14}$$

where M(**A**) is a function that depends on the specific constitutive model. The relation between stress tensor **T** and **A** is given by

$$\mathbf{T} = \frac{1-\beta}{\mathrm{Re}\,De}(\mathbf{A} - \mathbf{I}) \,. \tag{15}$$

that can rewritten as a relation between the tensor **S** and **A** given by

$$\mathbf{S} = \frac{1-\beta}{\mathrm{Re}\,D\epsilon} (\mathbf{A} - \mathbf{I} - 2\mathrm{De}\,\mathbf{D}).\tag{16}$$

A problem that challenges many researchers in computational rheology is solving Equation (12) or Equation (14) for high values of the Deborah number, *De* = *λ*/*tc*, where *t<sup>c</sup>* is a characteristic time of the flow. This problem occurs because all numerical methods are unstable for certain critical values of *De*. In order to overcome such failure, Fattal and Kupferman [37] proposed a reformulation of the differential constitutive equations into a equation for the matrix logarithm of the conformation tensor. Extending the ideas proposed by [37,38], ref. [14] presented a generic kernel-conformation tensor transformation that allows us to apply different kernel functions to the matrix transformation.

The reformulation of the tensor conformation was possible by the decomposition of the velocity gradient proposal by [37,38]

$$\nabla \mathbf{u}^{\mathsf{T}} = \boldsymbol{\Omega} + \mathbf{B} + \mathbf{N} \mathbf{A}^{-1} \tag{17}$$

where **Ω** and **N** are anti-symmetric tensors, **B** is symmetric and commutes with **A**. Thus, the constitutive equation based on the conformation tensor can be rewritten using the decomposition (17) as

$$\frac{\partial \mathbf{A}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{A} - (\Omega \mathbf{A} - \mathbf{A}\Omega) - 2\mathbf{B}\mathbf{A} = \frac{1}{De} \mathcal{M}(\mathbf{A}) \tag{18}$$

where M(**A**) is defined according to the viscoelastic model,

$$\mathcal{M}(\mathbf{A}) = \begin{cases} \mathbf{I} - \mathbf{A} & \text{Oldroyd-B}, \\ \mathbf{I} - \mathbf{A} - \mathfrak{a}(\mathbf{A} - \mathbf{I}) \cdot (\mathbf{A} - \mathbf{I}) & \text{Gieseksus}, \\ \left(1 + \frac{\varepsilon \operatorname{Re} D \varepsilon}{1 - \beta} \text{tr}(\mathbf{S})\right)(\mathbf{I} - \mathbf{A}) - 2\xi \operatorname{De}(\mathbf{B} - \mathbf{B}\mathbf{A}) & \text{PTT}. \end{cases} \tag{19}$$

Fattal and Kupferman showed that the matrix logarithm of the conformation tensor is a linear transformation of **A** and derived a constitutive equation from the Equation (18) in the function of the matrix logarithm. Afonso et al. proposed a generic *kernel-conformation* tensor transformation for a large class of differential constitutive models, in which the evolution equation for k(**A**), can be expressed in its tensorial formulations as

$$\frac{\text{Dk}(\mathbf{A})}{\text{Dt}} = \mathbf{\Omega} \mathbf{k}(\mathbf{A}) - \mathbb{k}(\mathbf{A})\mathbf{\Omega} + 2\mathbb{B} + \frac{1}{De}\mathbb{M} \tag{20}$$

where B and M are symmetric tensors constructed by the orthogonalization of the diagonal tensors **D**<sup>B</sup> and **D**M, respectively. These tensors can be constructed as

$$\begin{aligned} \mathbb{B} &= \mathbf{O} \mathbf{D}\_{\mathbb{B}} \mathbf{O}^{T} = & \mathbf{O} \widetilde{\mathbf{B}} \mathbf{A} \mathbf{J} \mathbf{O}^{T} \\ \mathbb{M} &= \mathbf{O} \mathbf{D}\_{\mathbb{M}} \mathbf{O}^{T} = & \mathbf{O} \mathcal{M}(\Lambda) \mathbf{J} \mathbf{O}^{T} . \end{aligned} \tag{21}$$

In Equations (21), **J** is the gradient matrix, a diagonal matrix of the form,

$$\mathbf{J} = \text{diag}\left(\frac{\partial \mathbb{k}(\lambda\_1)}{\partial \lambda\_1}; \frac{\partial \mathbb{k}(\lambda\_2)}{\partial \lambda\_2}; \frac{\partial \mathbb{k}(\lambda\_3)}{\partial \lambda\_3}\right). \tag{22}$$

#### **4. Verification Tests**

In this section, we address two test problems in terms of checking the HiG-Flow code for simulations of viscoelastic flows. One of the problems is the flow of a Phan-Thien– Tanner model fluid in a circular cross-section channel. The other test problem concerns the constitutive model of Oldroyd-B. The geometry used for this test was a driven cavity in two dimensions.

#### *4.1. Phan-Thien–Tanner Model Fluid Flow in a Pipe*

We consider a flow into a circular cylinder of radius *R*, where there exists only the axial velocity component *u*, which depends on the radial coordinate *r*. In addition, we consider that the fluid obeys the PTT fluid model [16] and the flow occurs in the *x* direction, the same as the cylinder axis. Here, we consider the known solutions available of the literature to the flow properties, namely velocity *u*, shear stress *Trx* and normal stress *Txx*. More detailed treatment about the analytical solution to this problem in a steady state, as well as the results verified here, can be found in [39–41]. Essentially, to obtain the viscoelastic component *Trx*, it is necessary to solve a cubic equation *T* 3 *rx* + 3*ATrx* − 2*B* = 0, whose its solution is given by

$$T\_{r\ge} = \left[\mathcal{B} + \left(A^3 + \mathcal{B}^2\right)^{1/2}\right]^{1/3} + \left[\mathcal{B} - \left(A^3 + \mathcal{B}^2\right)^{1/2}\right]^{1/3} \tag{23}$$

where *A* and *B* depends on the set of know parameters of flow:

$$A = \frac{\eta\_p^2}{6\varepsilon\lambda^2 \beta'}\tag{24}$$

$$B = -\frac{\eta\_p^3 \mu\_N}{\varepsilon \lambda^2 R^2 \beta} r.\tag{25}$$

In Equations (24) and (25), *η<sup>p</sup>* is the polymer viscosity, *ε* is the PTT parameter, *λ* is the relaxation time and *R* is the cylinder radius. The amount of solvent contribution is given by *β* = *ηs*/*η*0, the reference velocity is *u<sup>N</sup>* and *r* is the radial coordinate. After obtaining *Trx*, one can calculate the normal stress *Txx* and also integrate the equation of motion to determine the velocity field. The corresponding expressions are given by:

$$T\_{\rm xx} = \frac{2\lambda}{\eta\_p} T\_{\rm rx}^2 \tag{26}$$

$$u(r) = \frac{2u\_N}{\beta} \left[ 1 - \left(\frac{r}{R}\right)^2 \right] + f(A, B),\tag{27}$$

where *f*(*A*, *B*) is a function that depends on the parameters *A* and *B* given in (24) and (25), respectively. These simulation parameters can be adjusted when the polymer viscosity is fixed, then by varying *β* it is possible to control the amount of the solvent contribution. In addition, just as *β*, *ε* and *De* are input viscoelastic parameters, *λ* is adjusted by the Deborah number. To verify that the results are in agreement with [41], we set *ε* = 0.25 and *De* = 6.3, which corresponds to the reference *De<sup>N</sup>* = 1.0. Non-slip boundary conditions were used for the velocity in the cylinder wall. At the channel inlet, we imposed a parabolic velocity profile and at the outlet, the homogeneous Neummann boundary condition, that is, spatial variations in velocity are not allowed at the outflow. For pressure, a zero gradient was imposed on the wall and at the channel inlet while the outflow was fixed at a constant value. The initial conditions for the bulk domain is zero velocity.

Figures 2–4 show the velocity field, shear stress and normal stress, respectively, as a function of the amount of solvent, which is controlled by *β* parameter. When *β* ≈ 1, the polymer concentration is approximately zero and the fluid has Newtonian behavior. On the other hand, if *β* ≈ 0, the behavior of the PTT fluid resembles that of the Oldroyd-B model. The curves represented by down triangles corresponds to *β* = 0.9, up triangles to *β* = 0.5, circles to *β* = 0.2 and squares to *β* = 0.01. All these results are obtained by HiG-Flow simulation. They are in perfect agreement with the analytical curves represented by solid lines in Figures 2–4, which corresponds to the solutions given by Equations (23), (27) and (26) for *u*, *Trx* and *Txx*, respectively.

**Figure 2.** Velocity field for a PTT flow in a pipe.

**Figure 3.** Shear stress *Trx* for a PTT flow in a pipe.

**Figure 4.** Normal stress for a PTT flow in a pipe.

#### *4.2. 2D-Driven Cavity with Oldroyd-B Flow*

Flows in rectangular cavities have been studied since 1967 when the article [42] was published. In the 21st century, several studies of this type for viscoelastic fluids have been published [38,43–47]. The problem studied here has no analytical solution; however, there are results obtained by the cited authors that can be used for comparison. The data used for comparison in this study were provided by Palhares Junior et al. [47].

Figure 5 shows the schematic of lid-driven cavity. A parabolic profile velocity is imposed on the top. The aspect ratio is defined as Λ ≡ *H*/*L*. Some concerned works are listed in Table 1.

**Figure 5.** Illustration of lid-driven cavity. The parabolic velocity profile is imposed on the top. The aspect ratio is defined as Λ ≡ *H*/*L*.

The use of stabilization methods within the HiG-Flow system can be considered simple from the coder point of view because the code has been implemented in such a way that one can make a choice directly in the main simulation file. In this sense, the kernel conformation tensor is used to perform this operation, as previously described in the Section 3. Generically, the user simply writes the kernel k, the kernel derivative *d*k/*d***x** for the Jacobian transformation calculation and the kernel inverse k −1 correspondents. For the square root stabilizer used in these simulations, one can write Equations (28)–(30), respectively, as follows:

$$\mathbf{k}\_{\mathrm{ij}} = \sqrt{\mathbf{S}\_{\mathrm{ij}}} \tag{28}$$

$$\frac{d\mathbb{k}\_{ij}}{d\boldsymbol{\omega}} = \frac{1}{\mathbb{2}\sqrt{\mathbb{S}\_{ij}}},\tag{29}$$

$$\mathbf{k}\_{\mathrm{ij}}^{-1} = \mathbf{S}\_{\mathrm{ij}}^{2} \tag{30}$$

**Table 1.** Previous and current numerical studies concerned with lid-driven cavity flow of constant viscosity viscoelastic fluids.


For the simulations, the Reynolds number was fixed as *Re* = 0.01. The proportion of solvent in Oldroyd-B fluid was also fixed as *β* = 0.5. Simulations were performed for two different Deborah numbers, *De* = 1.0 and *De* = 2.0. On the top lid-driven section of the cavity, we imposed a parabolic velocity profile given by

$$u(x,t) = 8[1 + \tanh(8t - 4)]x^2 \left(1 - x^2\right)^2. \tag{31}$$

The other cavity walls are stationary and the non-slip condition is imposed over all of them. A regular mesh of 256x256 cells was used. The velocity component *u* and the normal stress *Txx* were plotted along the vertical line *x* = 0.5 while the velocity component *v* was obtained on the horizontal line *y* = 0.75. The (*x*, *y*) coordinates are scaled by the cavity side size *L* = 1 unit of length. The results are shown in Figures 6–8. The curves are represented by squares and circles corresponding to the HiG-Flow and Palhares Junior et al. results, respectively. The graphs indicate that the results are in good agreement.

**Figure 6.** Velocity field *u* obtained along the vertical line *x* = 0.5: (**a**) *De* = 1.0; (**b**) *De* = 2.0.

**Figure 7.** Velocity field *v* obtained along the horizontal line *y* = 0.75: (**a**) *De* = 1.0; (**b**) *De* = 2.0.

**Figure 8.** Normal stress *Txx* obtained along the vertical line *x* = 0.5: (**a**) *De* = 1.0; (**b**) *De* = 2.0.

#### **5. Simulation in Complex 3D Array of Channels**

In this section, we present the results obtained with the HiG-Flow system for flows in complex domains. We simulated an incompressible viscoelastic fluid flow in a complex array of microchannels, introducing some level of geometric complexity in the threedimensional flow domain.

The geometry, as well as boundary conditions, can be seen in Figure 9. The total width, length and height are set to be *W* = 0.8 mm, *L* = 2.4 mm and *H* = 0.4 mm, respectively. The inlet is a channel of 0.1 mm × 0.1 mm, where polymer at temperature is injected with a constant velocity of *U*in = 0.1 mm/s. Scaling this geometry by ` = 0.1 mm, and using *<sup>ν</sup>* <sup>≈</sup> <sup>10</sup>−<sup>4</sup> <sup>m</sup>2/s as the kinematic viscosity of polymer at room temperature, we end up with a Reynolds number of *Re* = ` *U*in/*ν* = 1.0. In this test, we used the PTT model with *β* = 0.5, *ε* = 0.25 and *ξ* = 0.0 for several values of *De* = [0 − 500] as viscoelastic parameters.

Streamlines for the flow of a Newtonian fluid can be observed in Figure 10. The result is qualitative, but demonstrates the robustness and applicability of this newly developed methodology. Several simulations using viscoelastic fluids for *De* = [0 − 500] were performed on the 3D complex geometry. We analyzed the complex fluid flow by observing the profiles of the polymeric stresses along the probe line near the 3D channel exits, as shown in Figure 11. The probe is aligned on the *y* direction at half channel height (along the *z* direction), orthogonal to the main flow direction near the channel exits.

The increasing values of elasticity, reflected on the value of Deborah number represented in Figure 12, affects the six components of the non-dimensional extra stress tensor along the probe line, with higher impact for the normal components, as the *Tzz* profiles. Nevertheless, given that no geometrical singularity is presented along the probe line, the maximum value for all extra stress components is not significant and slightly affected by the increase in elasticity.

**Figure 9.** Geometry for the complex 3D array of channels.

**Figure 10.** Streamlines for the complex 3D array of channels. The color scale varies from smallest (blue) to largest (red) velocity magnitude.

We used a computer with a 3.1 GHz Intel Core i7 Quad-Core Processor and 16 GB 2133 MHz LPDDR3 memory. The HiG-Flow software was used with four cores for all the calculation, and each simulation took 14 h of processing.

**Figure 11.** Probe views.

**Figure 12.** Tensor components: (**a**) *Txx*; (**b**) *Txy*; (**c**) *Tyy*; (**d**) *Tyz*; (**e**) *Tzz*; (**f**) *Txz*.

#### **6. Conclusions**

Tree-based grids bring the advantage of using fast Cartesian discretizations, such as finite differences, and the flexibility and accuracy of local mesh refinement. Most methods usually avoid this by limiting the mesh configuration (usually to graded quadtree/octree grids), reducing the number of cases to be treated locally. In this work, we employ a moving least squares meshless interpolation technique, allowing for more complex mesh configurations, while still keeping the overall order of accuracy. This technique was implemented in the HiG-Flow code to simulate Newtonian, generalized Newtonian and viscoelastic fluids flows. The code verification and testing was performed using numerical stabilizers for the Oldroyd-B flow solution in a 2D cavity and for a PTT fluid in a complex 3D geometry.

**Author Contributions:** Funding acquisition, A.C.; Methodology, A.C., A.M.A. and W.D.S.B.; Project administration, A.C.; Software, A.C., A.M.A. and W.D.S.B.; Supervision, A.C.; Validation, A.C., A.M.A. and W.D.S.B.; Writing—Original draft, A.C., A.M.A. and W.D.S.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Fundação para a Ciência e a Tecnologia I.P. (FCT), CEFT (Centro de Estudos de Fenómenos de Transporte) and by FEDER funds through COMPETE2020. Grants: PTDC/EMS-ENE/3362/2014, POCI-01-0145-FEDER-016665, UIDB/00013/2020 and UIDP/00013/ 2020. All authors are grateful for the financial support from the Brazilian Petroleum Agency (ANP)/Petrobras, grant 0050.0075367.12.9, and from the São Paulo Research Foundation (FAPESP), grants 2013/07375-0, 2017/21105-6 and 2020/02990-1. A. Castelo is also grateful for the financial support from the National Council for Scientific and Technological Development (CNPq), grant 307483/2017-7. Research was carried out using the computational resources of the Center for Mathematical Sciences Applied to Industry (CeMEAI), funded by FAPESP grant 2013/07375-0.

**Acknowledgments:** Castelo, A. acknowledges the support of the ICMC-Instituto de Ciências Matemáticas e de Computação, Departamento de Matemática e Estatística, USP, São Carlos, SP. Afonso, A.M. acknowledges the support of the Faculdade de Engenharia da Universidade do Porto and Bezerra. W. S. acknowledges the support of the UFGD-Universidade Federal da Grande Dourados, Instituto de Ciências Exatas e Tecnologia.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

