Next Article in Journal
Zinc Plus Biopolymer Coating Slows Nitrogen Release, Decreases Ammonia Volatilization from Urea and Improves Sunflower Productivity
Next Article in Special Issue
Taguchi Optimization of Roundness and Concentricity of a Plastic Injection Molded Barrel of a Telecentric Lens
Previous Article in Journal
Effect of Microcapsules with Waterborne Coating as Core Material on Properties of Coating for Tilia Europaea and Comparison with Other Microcapsules
Previous Article in Special Issue
A Hybrid Cooling Model Based on the Use of Newly Designed Fluted Conformal Cooling Channels and Fastcool Inserts for Green Molds
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Hierarchical Grid Solver for Simulation of Flows of Complex Fluids

by
Antonio Castelo
1,*,†,
Alexandre M. Afonso
2,† and
Wesley De Souza Bezerra
3,†
1
Departamento de Matemática Aplicada e Estatística, Instituto de Ciências Matemáticas e de Computação, Universidade de São Paulo, Cx.P. 668, São Carlos 13560-970, SP, Brazil
2
Centro de Estudos de Fenómenos de Transporte, Departamento de Engenharia Mecânica, Faculdade de Engenharia da Universidade do Porto, 4200-465 Porto, Portugal
3
Faculdade de Ciências Exatas e Tecnologia, Universidade Federal da Grande Dourados, Cx.P. 364, Dourados 79804-970, MS, Brazil
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Polymers 2021, 13(18), 3168; https://doi.org/10.3390/polym13183168
Submission received: 7 July 2021 / Revised: 7 September 2021 / Accepted: 10 September 2021 / Published: 18 September 2021
(This article belongs to the Special Issue Advanced Polymer Simulation and Processing)

Abstract

:
Tree-based grids bring the advantage of using fast Cartesian discretizations, such as finite differences, and the flexibility and accuracy of local mesh refinement. The main challenge is how to adapt the discretization stencil near the interfaces between grid elements of different sizes, which is usually solved by local high-order geometrical interpolations. 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, 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. Numerical tests and application to viscoelastic fluid flow simulations were performed to illustrate the flexibility and robustness of this new approach.

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,3,4,5,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,8,9,10,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 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 kernel-conformation 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 upper-convected 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-lid-driven 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.
Looking at Figure 1c, a second-order approximation to 2 U c y 2 can be given by (we assume the y axis is in the direction bottom → top):
2 U c y 2 1 δ y ( U t 2 U c + U b ) .
Note that U b is not known and must be obtained by interpolation (the same applies to U r ) using the following formula:
U b = k = 1 V b w k b   U k ,
where V b is the number of neighbor cells, which depends on the imposed accuracy of the method.
The weights w k b = w k b ( x ) are obtained through the moving least squares (MLS) method. In a set of n smooth interpolating functions that are linearly independent Φ i : R d R ( d = 2 , 3 ), we want to obtain the interpolated value u such that U b = U ( x ) = k = 1 n c i Φ i ( x ) = c t Φ .
Given m points x 1 , x 2 , , x m R d with m > n and m values u 1 = u ( x 1 ) , u 2 = u ( x 2 ) , , u m = u ( x m ) , to interpolate u in x using MLS consists in minimizing the error E ( c )
E ( c ) = U u 2 2 = i = 1 m ( U ( x i ) u i ) 2 1 x x i 2 .
Or,
E ( c ) = W P c W u 2 2 ,
where W = W ( x ) = δ i j 1 x x i 2 R m × m , P = Φ i ( x j ) R m × n and u = ( u 1 , u 2 , , u m ) .
The solution is given by
c ( x ) = ( W P ) W u
where ( · ) is the Moore–Penrose pseudo-inverse.
Decomposing ( W P ) into Q R we have that
W P = Q R 0 = Q     Q R 0   ,
where Q R m × m , Q R m × n , Q R m × m n , and R R n × n . This decomposition is then used to finally compute
c ( x ) = R 1 Q t W u   .
Then
U ( x ) = c t Φ = u t W Q R t Φ w = k = 1 m w k u k   ,
that is w = w ( x ) = W Q R 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
· u = 0 ,
and conservation of momentum
u t + u · u = p + 1 R e 2 u + · S + 1 F r 2 g + F ,
T = 2 ( 1 β ) R e D + S ,
where u is the velocity field, t is time, p is the pressure, R e is the Reynolds number, F r is the Froude number, g is the gravity force and F is the surface tension force and source force. The symbol D = 1 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 μ 0 = μ S + μ P 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:
T t + u · T u T · T + T · u = 1 D e M T .
where M T is defined by the viscoelastic model
M T = 2 ( 1 β ) R e D T Oldroyd - B , 2 ( 1 β ) R e D T α   R e   D e 1 β T · T Giesekus , 2 ( 1 β ) R e D 1 + ϵ   R e   D e 1 β tr T T ξ   D e T · D + D · T LPTT ,
where D e 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 ϵ , 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
A t + u · A A u + u T A = 1 D e M ( A ) ,
where M ( A ) is a function that depends on the specific constitutive model. The relation between stress tensor T and A is given by
T = 1 β R e   D e ( A I ) ,
that can rewritten as a relation between the tensor S and A given by
S = 1 β R e   D e A I 2 D e   D .
A problem that challenges many researchers in computational rheology is solving Equation (12) or Equation (14) for high values of the Deborah number, D e = λ / t c , where t c is a characteristic time of the flow. This problem occurs because all numerical methods are unstable for certain critical values of D e . 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]
u T = Ω + B + N A 1 ,
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
A t + u · A Ω A A Ω 2 B A = 1 D e M ( A ) ,
where M ( A ) is defined according to the viscoelastic model,
M ( A ) = I A Oldroyd - B , I A α A I · A I Giesekus , 1 + ϵ   R e   D e 1 β tr S I A 2 ξ   D e B B A PTT .
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
D k ( A ) D t = Ω k ( A ) k ( A ) Ω + 2 B + 1 D e M
where B and M are symmetric tensors constructed by the orthogonalization of the diagonal tensors D B and D M , respectively. These tensors can be constructed as
B = O D B O T = O B ˜ Λ J O T M = O D M O T = O M Λ J O T .
In Equations (21), J is the gradient matrix, a diagonal matrix of the form,
J = diag k λ 1 λ 1 ; k λ 2 λ 2 ; k λ 3 λ 3 .

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 T r x and normal stress T x x . 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,40,41]. Essentially, to obtain the viscoelastic component T r x , it is necessary to solve a cubic equation T r x 3 + 3 A T r x 2 B = 0 , whose its solution is given by
T r x = B + A 3 + B 2 1 / 2 1 / 3 + B A 3 + B 2 1 / 2 1 / 3 ,
where A and B depends on the set of know parameters of flow:
A = η p 2 6 ε λ 2 β ,
B = η p 3 u N ε λ 2 R 2 β r .
In Equations (24) and (25), η p 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 N and r is the radial coordinate. After obtaining T r x , one can calculate the normal stress T x x and also integrate the equation of motion to determine the velocity field. The corresponding expressions are given by:
T x x = 2 λ η p T r x 2 ,
u ( r ) = 2 u N β 1 r R 2 + f ( A , B ) ,
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 D e 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 D e = 6.3 , which corresponds to the reference D e N = 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.
Figure 2, Figure 3 and Figure 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 Figure 2, Figure 3 and Figure 4, which corresponds to the solutions given by Equations (23), (27) and (26) for u, T r x and T x x , respectively.

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,44,45,46,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.
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:
k i j = S i j ,
d k i j d x = 1 2 S i j ,
k i j 1 = S i j 2 ,
For the simulations, the Reynolds number was fixed as R e = 0.01 . The proportion of solvent in Oldroyd-B fluid was also fixed as β = 0.5 . Simulations were performed for two different Deborah numbers, D e = 1.0 and D e = 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 ( 8 t 4 ) x 2 1 x 2 2 .
The other cavity walls are stationary and the non-slip condition is imposed over all of them. A regular mesh of 256 × 256 cells was used. The velocity component u and the normal stress T x x 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 Figure 6, Figure 7 and Figure 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.

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 three-dimensional 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 ν 10 4 m 2 /s as the kinematic viscosity of polymer at room temperature, we end up with a Reynolds number of R e =   U in / ν = 1.0 . In this test, we used the PTT model with β = 0.5 , ε = 0.25 and ξ = 0.0 for several values of D e = [ 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 D e = [ 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 T z z 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.
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.

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

  1. Berger, M.J.; Oliger, J. Adaptive mesh refinement for hyperbolic partial differential equations. J. Comput. Phys. 1984, 53, 484–512. [Google Scholar] [CrossRef]
  2. Popinet, S. Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries. J. Comput. Phys. 2003, 190, 572–600. [Google Scholar] [CrossRef] [Green Version]
  3. Losasso, F.; Gibou, F.; Fedkiw, R. Simulating Water and Smoke with an Octree Data Structure. ACM Trans. Graph. 2004, 23, 457–462. [Google Scholar] [CrossRef]
  4. Olshanskii, M.A.; Terekhov, K.M.; Vassilevski, Y.V. An octree-based solver for the incompressible Navier Stokes equations with enhanced stability and low dissipation. Comput. Fluids 2013, 84, 231–246. [Google Scholar] [CrossRef] [Green Version]
  5. Guittet, A.; Theillard, M.; Gibou, F. A stable projection method for the incompressible Navier-Stokes equations on arbitrary geometries and adaptive Quad/Octrees. J. Comput. Phys. 2015, 292, 215–238. [Google Scholar] [CrossRef]
  6. Batty, C. A cell-centred finite volume method for the Poisson problem on non-graded quadtrees with second order accurate gradients. J. Comput. Phys. 2017, 331, 49–72. [Google Scholar] [CrossRef] [Green Version]
  7. Ding, H.; Shu, C.; Yeo, K.S.; Xu, D. Simulation of incompressible viscous flows past a circular cylinder by hybrid FD scheme and meshless least square-based finite difference method. Comput. Methods Appl. Mech. Eng. 2004, 193, 727–744. [Google Scholar] [CrossRef]
  8. Chew, C.S.; Yeo, K.S.; Shu, C. A generalized finite-difference (GFD) ALE scheme for incompressible flows around moving solid bodies on hybrid meshfree-Cartesian grids. J. Comput. Phys. 2006, 218, 510–548. [Google Scholar] [CrossRef]
  9. Min, C.; Gibou, F.; Ceniceros, H.D. A supra-convergent finite difference scheme for the variable coefficient Poisson equation on non-graded grids. J. Comput. Phys. 2006, 218, 123–140. [Google Scholar] [CrossRef]
  10. Ding, H.; Shu, C.; Cai, Q.D. Applications of stencil-adaptive finite difference method to incompressible viscous flows with curved boundary. Comput. Fluids 2007, 36, 786–793. [Google Scholar] [CrossRef]
  11. Wang, X.Y.; Yeo, K.S.; Chew, C.S.; Khoo, B.C. A SVD-GFD scheme for computing 3D incompressible viscous fluid flows. Comput. Fluids 2008, 37, 733–746. [Google Scholar] [CrossRef]
  12. Sousa, F.; Lages, C.; Ansoni, J.; Castelo, A.; Simao, A. A finite difference method with meshless interpolation for incompressible flows in non-graded tree-based grids. J. Comput. Phys. 2019, 396, 848–866. [Google Scholar] [CrossRef]
  13. Bird, R.B.; Armstrong, R.C.; Hassager, O. Dynamics of Polymeric Liquids, 2nd ed.; Fluid Mechanics; John Wiley and Sons Inc.: New York, NY, USA, 1987; Volume 1, p. 784. ISBN 0-471-80245-X. [Google Scholar]
  14. Afonso, A.M.; Pinho, F.T.; Alves, M.A. The kernel-conformation constitutive laws. J. Non-Newton. Fluid Mech. 2012, 167–168, 30–37. [Google Scholar] [CrossRef]
  15. Oldroyd, J.G. On the formulation of rheological equations of state. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 1950, 200, 523–541. [Google Scholar]
  16. Thien, N.P.; Tanner, R.I. A new constitutive equation derived from network theory. J. Non-Newton. Fluid Mech. 1977, 2, 353–365. [Google Scholar] [CrossRef]
  17. Dou, H.S.; Phan-Thien, N. The flow of an Oldroyd-B fluid past a cylinder in a channel: Adaptive viscosity vorticity (DAVSS-ω) formulation. J. Non-Newton. Fluid Mech. 1999, 87, 47–73. [Google Scholar] [CrossRef]
  18. Alves, M.; Pinho, F.; Oliveira, P. Effect of a high-resolution differencing scheme on finite-volume predictions of viscoelastic flows. J. Non-Newton. Fluid Mech. 2000, 93, 287–314. [Google Scholar] [CrossRef]
  19. Harten, A. High resolution schemes for hyperbolic conservation laws. J. Comput. Phys. 1983, 49, 357–393. [Google Scholar] [CrossRef] [Green Version]
  20. Alves, M.; Oliveira, P.; Pinho, F. A convergent and universally bounded interpolation scheme for the treatment of advection. Int. J. Numer. Methods Fluids 2003, 41, 47–75. [Google Scholar] [CrossRef]
  21. Alves, M.A.; Oliveira, P.J.; Pinho, F.T. Benchmark solutions for the flow of Oldroyd-B and PTT fluids in planar contractions. J. Non-Newton. Fluid Mech. 2003, 110, 45–75. [Google Scholar] [CrossRef]
  22. Chinyoka, T.; Renardy, Y.; Renardy, M.; Khismatullin, D. Two-dimensional study of drop deformation under simple shear for Oldroyd-B liquids. J. Non-Newton. Fluid Mech. 2005, 130, 45–56. [Google Scholar] [CrossRef]
  23. Tomé, M.F.; Paulo, G.S.d.; Pinho, F.; Alves, M. Numerical solution of the PTT constitutive equation for unsteady three-dimensional free surface flows. J. Non-Newton. Fluid Mech. 2010, 165, 247–262. [Google Scholar] [CrossRef]
  24. Mompean, G.; Thais, L.; Tomé, M.F.; Castelo, A. Numerical prediction of three-dimensional time-dependent viscoelastic extrudate swell using differential and algebraic models. Comput. Fluids 2011, 44, 68–78. [Google Scholar] [CrossRef]
  25. Tomé, M.F.; Castelo, A.; Afonso, A.M.; Alves, M.A.; Pinho, F.T.d. Application of the log-conformation tensor to three-dimensional time-dependent free surface flows. J. Non-Newton. Fluid Mech. 2012, 175, 44–54. [Google Scholar] [CrossRef]
  26. Oishi, C.M.; Martins, F.P.; Tomé, M.F.; Alves, M.A. Numerical simulation of drop impact and jet buckling problems using the eXtended Pom–Pom model. J. Non-Newton. Fluid Mech. 2012, 169, 91–103. [Google Scholar] [CrossRef]
  27. Paulo, G.S.d.; Oishi, C.M.; Tomé, M.F.; Alves, M.A.; Pinho, F. Numerical solution of the FENE-CR model in complex flows. J. Non-Newton. Fluid Mech. 2014, 204, 50–61. [Google Scholar] [CrossRef]
  28. Tomé, M.F.; Araujo, M.T.d.; Evans, J.D.; McKee, S. Numerical solution of the Giesekus model for incompressible free surface flows without solvent viscosity. J. Non-Newton. Fluid Mech. 2019, 263, 104–119. [Google Scholar] [CrossRef]
  29. Bezerra, W.D.S.; Castelo, A.; Afonso, A.M. Numerical Study of Electro-Osmotic Fluid Flow and Vortex Formation. Micromachines 2019, 10, 796. [Google Scholar] [CrossRef] [Green Version]
  30. Shojaei, A.; Galvanetto, U.; Rabczuk, T.; Jenabi, A.; Zaccariotto, M. A generalized finite difference method based on the Peridynamic differential operator for the solution of problems in bounded and unbounded domains. Comput. Methods Appl. Mech. Eng. 2019, 343, 100–126. [Google Scholar] [CrossRef]
  31. Varchanis, S.; Syrakos, A.; Dimakopoulos, Y.; Tsamopoulos, J. A new finite element formulation for viscoelastic flows: Circumventing simultaneously the LBB condition and the high-Weissenberg number problem. J. Non-Newton. Fluid Mech. 2019, 267, 78–97. [Google Scholar] [CrossRef]
  32. Bertoco, J.; de Araújo, M.S.; Leiva, R.T.; Sánchez, H.A.; Castelo, A. Numerical Simulation of KBKZ Integral Constitutive Equations in Hierarchical Grids. Appl. Sci. 2021, 11, 4875. [Google Scholar] [CrossRef]
  33. Chávez-Negrete, C.; Santana-Quinteros, D.; Domínguez-Mota, F. A Solution of Richards’ Equation by Generalized Finite Differences for Stationary Flow in a Dam. Mathematics 2021, 9, 1604. [Google Scholar] [CrossRef]
  34. McKee, S.; Tomé, M.F.; Ferreira, V.G.; Cuminato, J.A.; Castelo, A.; Sousa, F.; Mangiavacchi, N. The MAC method. Comput. Fluids 2008, 37, 907–930. [Google Scholar] [CrossRef]
  35. Phan-Thien, N.; Tanner, R. A non-linear network viscoelastic model. J. Rheol. 1978, 22, 259–283. [Google Scholar] [CrossRef]
  36. Giesekus, H. A simple constitutive equation for polymer fluids based on the concept of defor mation-dependent tensional mobility. J. Non-Newton. Fluid Mech. 1982, 11, 69–109. [Google Scholar] [CrossRef]
  37. Fattal, R.; Kupferman, R. Constitutive laws for the matrix-logarithm the conformation tensor. J. Non-Newton. Fluid Mech. 2004, 123, 281–285. [Google Scholar] [CrossRef]
  38. Fattal, R.; Kupferman, R. Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation. J. Non-Newton. Fluid Mech. 2005, 126, 23–37. [Google Scholar] [CrossRef]
  39. Oliveira, P.J.; Pinho, F.T. Analytical solution for fully developed channel and pipe flow of Phan-Thien–Tanner fluids. J. Fluid Mech. 1999, 387, 271–280. [Google Scholar] [CrossRef] [Green Version]
  40. Alves, M.A.; Pinho, F.T.; Oliveira, P.J. Study of steady pipe and channel flows of a single-mode Phan-Thien–Tanner fluid. J. Non-Newton. Fluid Mech. 2001, 101, 55–76. [Google Scholar] [CrossRef]
  41. Cruz, D.; Pinho, F.T.d.; Oliveira, P. Analytical solutions for fully developed laminar flow of some viscoelastic liquids with a Newtonian solvent contribution. J. Non-Newton. Fluid Mech. 2005, 132, 28–35. [Google Scholar] [CrossRef] [Green Version]
  42. Pan, F.; Acrivos, A. Steady flows in rectangular cavities. J. Fluid Mech. 1967, 28, 643–655. [Google Scholar] [CrossRef]
  43. Pan, T.W.; Hao, J.; Glowinski, R. On the simulation of a time-dependent cavity flow of an Oldroyd-B fluid. Int. J. Numer. Methods Fluids 2009, 60, 791–808. [Google Scholar] [CrossRef]
  44. Yapici, K.; Karasozen, B.; Uludag, Y. Finite volume simulation of viscoelastic laminar flow in a lid-driven cavity. J. Non-Newton. Fluid Mech. 2009, 164, 51–65. [Google Scholar] [CrossRef]
  45. Poole, R.; Afonso, A.; Pinho, F.; Oliveira, P.; Alves, M. Scaling of purely-elastic instabilities in viscoelastic lid-driven cavity flow. In Proceedings of the XVIth International Workshop for Numerical Methods in Non-Newtonian Flows, Northhampton, MA, USA, 13–16 June 2010. [Google Scholar]
  46. Habla, F.; Tan, M.W.; Haßlberger, J.; Hinrichsen, O. Numerical simulation of the viscoelastic flow in a three-dimensional lid-driven cavity using the log-conformation reformulation in OpenFOAM®. J. Non-Newton. Fluid Mech. 2014, 212, 47–62. [Google Scholar] [CrossRef]
  47. Junior, I.L.P.; Oishi, C.M.; Afonso, A.M.; Alves, M.A.; Pinho, F.T. Numerical study of the square-root conformation tensor formulation for confined and free-surface viscoelastic fluid flows. Adv. Model. Simul. Eng. Sci. 2016, 3, 2. [Google Scholar] [CrossRef] [Green Version]
  48. Grillet, A.M.; Yang, B.; Khomami, B.; Shaqfeh, E.S. Modeling of viscoelastic lid driven cavity flow using finite element simulations. J. Non-Newton. Fluid Mech. 1999, 88, 99–131. [Google Scholar] [CrossRef]
  49. Comminal, R.; Spangenberg, J.; Hattel, J.H. Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation. J. Non-Newton. Fluid Mech. 2015, 223, 37–61. [Google Scholar] [CrossRef]
  50. Martins, F.; Oishi, C.M.; Afonso, A.M.; Alves, M.A. A numerical study of the kernel-conformation transformation for transient viscoelastic fluid flows. J. Comput. Phys. 2015, 302, 653–673. [Google Scholar] [CrossRef] [Green Version]
  51. Dalal, S.; Tomar, G.; Dutta, P. Numerical study of driven flows of shear thinning viscoelastic fluids in rectangular cavities. J. Non-Newton. Fluid Mech. 2016, 229, 59–78. [Google Scholar] [CrossRef]
Figure 1. (a) Example of hierarchical grid. (b) Tree data structure. (c) Finite difference method.
Figure 1. (a) Example of hierarchical grid. (b) Tree data structure. (c) Finite difference method.
Polymers 13 03168 g001
Figure 2. Velocity field for a PTT flow in a pipe.
Figure 2. Velocity field for a PTT flow in a pipe.
Polymers 13 03168 g002
Figure 3. Shear stress T r x for a PTT flow in a pipe.
Figure 3. Shear stress T r x for a PTT flow in a pipe.
Polymers 13 03168 g003
Figure 4. Normal stress for a PTT flow in a pipe.
Figure 4. Normal stress for a PTT flow in a pipe.
Polymers 13 03168 g004
Figure 5. Illustration of lid-driven cavity. The parabolic velocity profile is imposed on the top. The aspect ratio is defined as Λ H / L .
Figure 5. Illustration of lid-driven cavity. The parabolic velocity profile is imposed on the top. The aspect ratio is defined as Λ H / L .
Polymers 13 03168 g005
Figure 6. Velocity field u obtained along the vertical line x = 0.5 : (a) D e = 1.0 ; (b) D e = 2.0 .
Figure 6. Velocity field u obtained along the vertical line x = 0.5 : (a) D e = 1.0 ; (b) D e = 2.0 .
Polymers 13 03168 g006
Figure 7. Velocity field v obtained along the horizontal line y = 0.75 : (a) D e = 1.0 ; (b) D e = 2.0 .
Figure 7. Velocity field v obtained along the horizontal line y = 0.75 : (a) D e = 1.0 ; (b) D e = 2.0 .
Polymers 13 03168 g007
Figure 8. Normal stress T x x obtained along the vertical line x = 0.5 : (a) D e = 1.0 ; (b) D e = 2.0 .
Figure 8. Normal stress T x x obtained along the vertical line x = 0.5 : (a) D e = 1.0 ; (b) D e = 2.0 .
Polymers 13 03168 g008
Figure 9. Geometry for the complex 3D array of channels.
Figure 9. Geometry for the complex 3D array of channels.
Polymers 13 03168 g009
Figure 10. Streamlines for the complex 3D array of channels. The color scale varies from smallest (blue) to largest (red) velocity magnitude.
Figure 10. Streamlines for the complex 3D array of channels. The color scale varies from smallest (blue) to largest (red) velocity magnitude.
Polymers 13 03168 g010
Figure 11. Probe views.
Figure 11. Probe views.
Polymers 13 03168 g011
Figure 12. Tensor components: (a) T x x ; (b) T x y ; (c) T y y ; (d) T y z ; (e) T z z ; (f) T x z .
Figure 12. Tensor components: (a) T x x ; (b) T x y ; (c) T y y ; (d) T y z ; (e) T z z ; (f) T x z .
Polymers 13 03168 g012
Table 1. Previous and current numerical studies concerned with lid-driven cavity flow of constant viscosity viscoelastic fluids.
Table 1. Previous and current numerical studies concerned with lid-driven cavity flow of constant viscosity viscoelastic fluids.
ReferenceAspect RatiosConstitutive EquationDeRegularizationNotes
Grillet et al. [48]0.5, 1.0, 3.0FENE-CR, L 2 =   25 , 100 , 400 0.24 Leakage at corners A and BFE
Fattal and Kupferman [38]1.0Oldroyd-B, β = 0.5 1.0 , 2.0 , 3.0 , 5.0 u ( x ) = 16 U x 2 ( 1 x ) 2 FD, Log conformation technique
Pan et al. [43]1.0Oldroyd-B, β = 0.5 0.5 , 1.0 u ( x ) = 16 U x 2 ( 1 x ) 2 FE, Log conformation technique
Yapici et al. [44]1.0Oldroyd-B, β = 0.3 1.0 NoFV, First-order upwind
Habla et al. [46]1.0Oldroyd-B, β = 0.5 0 to 2 u ( x , z ) = 128 1 + tanh 8 ( t 1 / 2 ) x 2 ( 1 x ) 2 z 2 ( 1 z ) 2 FV, 3D, Log conformation technique, CUBISTA
Comminal et al. [49]1.0Oldroyd-B, β = 0.5 0.25 to 10 u ( x ) = 16 U x 2 ( 1 x ) 2 FD/FV, Log-conformation, stream function
Martins et al. [50]1.0Oldroyd-B, β = 0.5 0.5 , 1.0 , 2.0 u ( x ) = 16 U x 2 ( 1 x ) 2 FD, Kernel-conformation technique
Dalal et al. [51]1.0Oldroyd-B, β = 0.5 1.0 u ( x ) = 16 U x 2 ( 1 x ) 2 FD, Symmetric square root
Palhares Junior et al. [47]1.0Oldroyd-B, β = 0.5 1.0 , 2.0 u ( x , t ) = 8 1 + tanh ( 8 t 4 ) x 2 1 x 2 2 FD, Symmetric square root
Current work1.0Oldroyd-B, β = 0.5 1.0 , 2.0 u ( x , t ) = 8 1 + tanh ( 8 t 4 ) x 2 1 x 2 2 FD, Kernel-conformation technique
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

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

AMA Style

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

Chicago/Turabian Style

Castelo, Antonio, Alexandre M. Afonso, and Wesley De Souza Bezerra. 2021. "A Hierarchical Grid Solver for Simulation of Flows of Complex Fluids" Polymers 13, no. 18: 3168. https://doi.org/10.3390/polym13183168

APA Style

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

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop