A Parallel Adaptive Viscoelastic Flow Solver with

0 downloads 0 Views 5MB Size Report
A parallel adaptive mesh refinement algorithm has been incorporated into the ... limits the number of approaches for hexahedral mesh refinement and makes it ...
A Parallel Adaptive Viscoelastic Flow Solver with Template Based Dynamic Mesh Refinement Evren ONER † and Mehmet SAHIN ‡ Faculty of Aeronautics and Astronautics, Astronautical Engineering Department, Istanbul Technical University, Maslak, Istanbul, 34469, TURKEY Abstract A parallel adaptive mesh refinement algorithm has been incorporated into the side-centered finite volume method [Sahin, A stable unstructured finite volume method for parallel large-scale viscoelastic fluid flow calculations. J. Non-Newtonian Fluid Mech., 166 (2011) 779–791] in order to obtain highly accurate numerical results for viscoelastic fluid flow problems. The present recursive mesh refinement algorithm is based on a conformal refinement of unstructured quadrilateral/hexahedral elements with templates based on 1:3 refinement of edges. In order to transfer cell-centered data between source and target meshes, a second-order conservative interpolation (remapping) technique similar to the work of Menon and Schmidt [Supermesh construction for conservative interpolation on unstructured meshes: An extension to cell–centered finite-volume variables. Comput. Methods Appl. Mech. Eng., 200 (2011), 2797–2804] is employed and the approach has been extended for side-centered data. The proposed framework has been applied to the classical benchmark problem of an Oldroyd-B fluid past a confined circular cylinder in a rectangular channel and a sphere falling in a circular tube. The calculations confirm that high accuracy can be achieved with the present adaptive mesh refinement.

Keywords: Finite Volume Method; Unstructured Meshes; Adaptive Refinement; Conservative Interpolation; Large-Scale Simulations; Viscoelastic Fluids.

1

INTRODUCTION

Accurate and stable numerical methods are particularly important in viscoelastic flow simulations. Although much progress has been made in recent years toward the accurate and efficient simulation of viscoelastic fluids, many difficulties still persist. One major difficulty in the numerical simulation of viscoelastic fluid flow problems is that it requires to resolve physical features such as geometric singularities, viscoelastic wakes and boundary layers. These features are generally associated with steep gradients in the flow variables, embedded in or adjacent to regions where these flow variables change more smoothly. An adaptive mesh technique takes advantage of these low gradient regions in order to avoid a fine meshing of the entire computational domain. Otherwise, discretized viscoelastic flow problems yield very large linear systems and solving such systems has also proven to be challenging due to computer processing speed and memory limitations. Therefore, considerable effort has been given to the implementation of adaptive procedures and a posteriori error estimates. Jin and Tanner [29] presented an element residual method for obtaining a posteriori error estimate of a finite element method and an adaptive h−refinement procedure. Mutlu et al. [37] used an error estimate based on velocity gradients and an h− refinement based on a Delaunay procedure for an Oldroyd-B and a Phan-Thien-Tanner model. Chauviere and Owens [9] used an error indicator in conjunction with a p−adaptive spectral element method to solve the benchmark problem of an Oldroyd-B fluid past a sphere in a tube. Roquet and Saramito [42] combined a high order mixed finite element approximation, an anisotropic auto-adaptive mesh procedure, and the augmented Lagrangian method to solve the Bingham fluid around a cylinder falling at constant speed between two parallel plates. Cai and Westphal [7] presented an adaptive mixed least-squares finite element method for viscoelastic fluids of Oldroyd-B type. Gu´enette et al. [25] combined a finite element method based † E-mail: ‡ E-mail:

[email protected] [email protected]

1

on log-conformation with an anisotropic adaptive remeshing technique for viscoelastic fluid flows. The authors reported mesh convergence up to W e = 0.7 for the classical benchmark problem of an Oldroyd-B fluid past a cylinder in a confined channel. In this paper we present a framework for parallel adaptive mesh refinement in order to obtain accurate numerical simulations for viscoelastic fluid flows. The present mesh refinement algorithm is based on a conformal refinement of unstructured quadrilateral/hexahedral elements with templates based on 1:3 refinement of edges. The reason to use all-quadrilateral/hexahedral elements is that quadrilateral and hexahedral elements have been proved to be useful for finite element and finite volume methods, and for some applications they are preferred to triangles or tetrahedra [46]. In addition, the template based refinement on quadrilateral/hexahedral elements can be performed with a minimum communication cost which is particularly suitable for parallel large-scale simulations of viscoelastic fluid flow problems. Furthermore, the template based refinement on quadrilateral/hexahedral elements also maintains high quality elements in the refined mesh. In template based refinement, various sets of templates that define possible element divisions are used to avoid hanging nodes and maintain the connectivity between the elements in the coarse and the refined levels. While it is possible to use point insertion (Delaunay-type algorithms) or advancing-front methods to refine triangular and tetrahedral meshes [35], it is not possible to refine quadrilateral and hexahedral meshes by these methods without creating triangular or tetrahedral elements in the process. The inability to use point insertion and advancing-front type algorithms limits the number of approaches for hexahedral mesh refinement and makes it quite difficult to implement. One widely used technique to solve the conformal hexahedral mesh refinement problem is the template-based refinement method proposed by Schneiders [43] using two-refinement (bisection of edges) or three-refinement (trisection of edges) templates. Although the two-refinement process (see, for example, [14, 52]) leads to fewer elements with a better transition region, it is rather complex and introduces conformability issues beyond the unrefined region. On the other hand, the three-refinement process is completely local and it is ideal for parallel implementation. Schneiders’ scheme in its original form has been used by several researchers [51, 33] but its inability to refine concave regions efficiently lead to the introduction of new set of templates by other researchers. This limitation is due to missing or undefined templates which can not be constructed using pure hexahedral elements with sufficient quality [41]. Garimella [23] proposed new set of templates for an edge based marking scheme along with a coarsening procedure for quadrilateral elements. Zhang and Zhao [51] defined new set of hexahedral templates to overcome the concavity problem of Schneiders’ templates. Their approach uses two templates from the work of Schneiders and introduces four new templates. Although their results show that the propagation of refinement regions into unrefined regions is avoided, there is no detailed description of the templates they used. The works of Ito et al. [27] and Sun et al. [44] are similar to that of Zhang and Zao [51] but their templates are described in detail. The multi level refinement implementation in this study uses the templates defined in the works of Schneiders [43] and Ito et al. [27] for two- and three-dimensions, respectively. The use of an adaptive refinement technique requires an interpolation (remapping) algorithm in order to transfer the solution fields from source mesh to target mesh. For a proper coupling between the flow solver and the refinement algorithm, the remapping technique has to be conservative and sufficiently accurate due to stability restrictions. Simple interpolation methods can be accurate, but usually do not guarantee conservation, causing undesirable additional perturbations to the numerical solution beyond what truncation error would indicate. Various approaches have been proposed for conservative remapping in the literature [13, 24, 22, 3, 20, 38]. Grandy [24] geometrically computed the intersection volumes between donor (source) and target meshes and used them to transfer the solution fields, which are defined on a donor mesh, to an unrelated target mesh. However, the remapping technique was first-order accurate. Dukowicz and Padial [13] used a linear function within each cell of the original mesh and presented a second-order conservative global remapping between three-dimensional meshes, each composed of arbitrary hexahedral cells. Margolin and Shashkov [36] applied a conservative interpolation to an Arbitrary Lagrangian-Eulerian (ALE) algorithm, which is based on partitioning cells of the updated mesh into components of cells from the old mesh and swept regions from neighbouring cells which is known as ”swept region remap”. Garimella et al. [22] extended the conservative interpolation (remapping) on three-dimensional generalized unstructured polyhedral meshes for use in continuous rezone ALE simulations. Farrell et al. [20] presented a local Galerkin projection approach within a finite element framework in order to conservatively interpolate between discrete fields. The approach uses an 2

intermediate mesh, known as the ”supermesh”, consisting of the intersections of the target and donor elements. Then the numerical integrations are performed on the supermesh, which is similar to that of the ”intersection based remap” [24, 13]. Menon and Schmidt [38] used the same supermesh approach for conservatively transferring, or remapping, cell-centered variable fields from one mesh to another with second-order accuracy. The method is applicable to any polyhedral source or target mesh. In the current study, the approach of Menon and Schmidt [38] is used to transfer the cell centered discrete solution fields. However, the approach has been extended for the face/edge centered data in order to interpolate the face/edge centered discrete solution fields. Conservative interpolations have been also considered for transferring loads along the interface in coupled fluid-solid simulations involving non-matching meshes [19, 28] in the literature. To the authors’ best knowledge, the conservative interpolations have not been used to transfer the solution fields for viscoelastic fluid flow problems. The present refinement and remapping algorithms are incorporated into our ViscoSolve algorithm [47, 48] which is based on the side-centered finite volume method where the velocity vector components are defined at the mid-point of each cell face, while the pressure term and the extra stress tensor are defined at element centroids. The present arrangement of the primitive variables leads to a stable numerical scheme and it does not require any ad-hoc modifications in order to enhance the pressure-velocity-stress coupling. The time stepping algorithm used decouples the calculation of the polymeric stress by solution of a hyperbolic constitutive equation from the evolution of the velocity and pressure fields by solution of a generalized Stokes problem as in the work of Caola et al. [8]. The numerical algorithm also uses the log-conformation representation proposed in [21] in order improve the limiting Weissenberg numbers due to the classical high Weissenberg number problem [31]. The resulting large algebraic linear systems are solved using the FGMRES(m) Krylov iterative method [45] with the restricted additive Schwarz preconditioner for the extra stress tensor and the geometric non-nested multilevel preconditioner for the Stokes system [47]. The present multilevel preconditioner for the Stokes system is essential for parallel scalable viscoelastic flow computations because the solution of Stokes flow is a serious bottleneck in performing parallel large-scale viscoelastic numerical simulations. This is because, as it is well known, one-level methods lead to non-scalable solvers since they cause an increase in the number of iterations as the number of sub-domains is increased. However, it is not possible to apply the standard multigrid methods with classical smoothing techniques (e.g., Jacobi, Gauss-Seidel) for the coupled iterative solution of the momentum and continuity equations because of the zero-block in the saddle point problem. In order to avoid the zero-block in the saddle point problem, we use an upper triangular right preconditioner which results in a scaled discrete Laplacian instead of a zero block in the original system. The implementation of the preconditioned Krylov subspace algorithm, matrix-matrix multiplication and the multilevel preconditioner were carried out using the PETSc software package [4] developed at the Argonne National Laboratories for improving the efficiency of the parallel code. The METIS library [30] is used to partition the computational domain into sub-domains for a balanced domain decomposition. The numerical method has been successfully used to simulate the purely viscoelastic wake instabilities behind a confined circular cylinder for the first time [48]. The current paper is organized as follows: Section 2 provides some details on the unstructured sidecentered finite volume method, parallel adaptive mesh refinement and conservative interpolation (remapping). In Section 3, the interpolation algorithm is initially tested on several cell- and side-centered data. Then the proposed method is applied to an Oldroyd-B fluid past a confined circular cylinder in a rectangular channel and a sphere falling in a circular tube using several different adaptation criterions. Concluding remarks are provided in Section 4.

2

PROBLEM DEFINITION and NUMERICAL METHOD

In this section, the implementation details of the unstructured finite volume algorithm, template based recursive adaptive mesh refinement algorithm and conservative interpolation (remapping) algorithm for cell- and side-centered discrete solution fields are described in detail.

3

2.1

Finite Volume Algorithm

The governing equations for three-dimensional unsteady flow of an incompressible and isothermal OldroydB fluid can be written in dimensionless form as follows: the continuity equation −∇ · u = 0 the momentum equations

[ Re

] ∂u + (u · ∇)u + ∇p = β∇2 u + ∇ · T ∂t

and the constitutive equation for the Oldroyd-B model [ ] ∂T ⊤ We + (u · ∇)T − (∇u) · T − T · ∇u = (1 − β)(∇u + ∇u⊤ ) − T ∂t

(1)

(2)

(3)

In these equations u represents the velocity vector, p is the pressure and T is the extra stress tensor. The dimensionless parameters are the Reynolds number Re, the Weissenberg number W e and the viscosity ratio β. The governing equations above are discretized using the side-centered finite volume method on conformal unstructured quadrilateral/hexahedral meshes as described in [47] in detail.

2.2

Template Based Mesh Refinement Algorithm

The most commonly used technique to solve the conformal hexahedral mesh refinement problem is the template-based hierarchical refinement method proposed by Schneiders [43]. In template-based refinement, various sets of templates that define possible element divisions are used to avoid hanging nodes and maintain the connectivity between the elements in the refined level and the original coarse level. In order to perform multi-level refinement and coarsening, the following hierarchical data structure is stored in an array with the following data: • NODES: Stores the node indices of the element. • PARENT : Stores the parent element of the element. • CHILDREN: Stores the children of the element. • LEVEL: Stores the refinement level of the element. • REGULAR: Stores if the element is suitable for further refinement. Set to 1 if the element can be refined, 0 otherwise. • ACTIVE: Stores if the element is needed for the final level mesh. • MARK: Stores if the element is marked for refinement. The storage of the data above requires an integer array of length 18 for each quadrilateral element in two-dimensions and an integer array of length 50 for each hexahedral element in three-dimensions. Template based refinement procedure starts with setting target refinement levels (LEVEL) by an adaptation criterion. The target element level of an element can be above or below its current level. These target refinement levels are then adjusted so that there is at most one level difference between neighboring elements. The one-level difference rule ensures that the mesh size is smoothly graded. After the level adjustment, the mesh adaptation algorithm performs coarsening of elements before refinement. Coarsening is performed on elements whose current refinement level is higher than the target refinement level. The coarsening algorithm is based on the method described in [23]. In this approach, coarsening is performed using the knowledge of the hierarchical data structure of the adapted mesh, i.e., if an element is to be deleted then its children are also deleted and the parent element is restored. Therefore, the coarsening process can not proceed beyond the initial coarse level. Before actual deletion of elements, the target levels of elements are adjusted again to ensure that there is not more than one level of difference between two adjacent elements and the target levels of children are consistent. After coarsening, element nodes are marked by propagating the target LEVEL of the element to all its nodes and an appropriate 4

[a]

[b] Figure 1: Two-dimensional refinement templates used by Schneiders (1996) [a] and Ito et al. (2009) [b] . template is applied based on the number of marked nodes. It turns out that the templates used to refine an element can be uniquely defined based on its node markings and the refinement process can be done in a purely local element-by-element manner. For example, if all the nodes of the element are marked, uniform refinement template is used. Figure 1 shows all possible node markings and corresponding twodimensional templates proposed by Schneiders [43] and Ito et al. [27]. During the refinement process, care must be taken to avoid duplicate node creation along element edges. For this purpose, each refinement operation stores edge indices of the element, along with the indices of newly created nodes on that edge in EDGEMARK array. Before each element refinement, EDGEMARK array is checked to see if there are nodes already created on the element edges. After refining an element, the hierarchical data structure is updated. Refined element is marked as INACTIVE, as this element is no longer part of the mesh at this refinement level. Newly created elements are set as CHILDREN of the refined element and the original element is marked as PARENT of the children. Current refinement LEVEL of newly created elements is increased by one. Depending on the applied template, new elements can be marked as REGULAR or IRREGULAR. The application of uniform refinement template produces REGULAR elements. All other templates produce IRREGULAR elements. This classification is used to prevent successive refinement of low quality IRREGULAR elements that belong to the transition region, as further refinement of these elements can degrade mesh quality. The present refinement procedure continues to iterate over the mesh elements until all elements have reached their target level. For the two-dimensional parallel refinement algorithm, the templates proposed by Schneiders [43] are implemented due to their overall mesh quality. The main problem in three-dimensional template based refinement is the high number of possible template configurations. There are a total of 22 different node marking configurations after considering symmetry and rotations [43]. However, it is not possible define all-hexahedral refinement templates for all node marking configurations due to geometric restrictions. Schneiders [43] proposed only four templates for node, edge, face and volume refinement as shown in Figure 2-[a]. In order to refine using only these four templates, a preprocessing step is needed to convert arbitrary node marking configurations to one of these four templates. However, by using only four templates, the propagation of refinement region into unrefined regions is unavoidable for concave domains, known as the concavity problem, and this causes refinement of the entire domain or most parts of the domain. Zhang and Zhao [51] defined new set of hexahedral templates to overcome the concavity problem of Schneiders templates [51]. Their approach uses two templates from the work of Schneiders and introduces four new templates. Although their results show that the propagation of refinement regions into unrefined regions is avoided, there is no detailed description of the templates they used. The works of Ito et al. [27] and Sun et al. [44] are similar to that of Zhang and Zao [51] but their templates are described in detail and have face-diagonal symmetry which makes implementation easier. The templates used by Ito et al. [27] are partially able to handle concave regions without expanding into unrefined regions. In their approach, there are a total of five templates and arbitrary node marking configurations have to be converted to one of these templates

5

[a]

[b] Figure 2: Three-dimensional refinement templates used by Schneiders [43] [a] and Ito et al. [27] [b] . as shown in Table 1. The application of refinement algorithm given by Ito et al. [27] leads a lower number of hexahedral elements but the overall quality of the resulting mesh will be slightly lower than that of Schneiders algorithm [43]. Therefore, the three-dimensional refinement algorithm is based on the templates proposed by Ito et al. [27]. Three-dimensional refinement templates proposed by Ito et al. are shown in Figure 2-[b]. The refinement algorithm is very similar to the two-dimensional case. However, the size of NODES and CHILDREN data arrays must be increased in order to store additional element data in three-dimensions. The EDGEMARK array used to avoid duplicate node creation on the element edges needs to be supplemented with a FACEMARK array in order to avoid duplicate node creation on element faces. The algorithm for the three-dimensional template based refinement is provided in Algorithm 1. The three-dimensional comparison of the algorithms based on Schneiders [43] and Ito et al. [27] is shown in Figure 3.

2.3

Conservative Interpolation Algorithm

The use of adaptive refinement requires an interpolation of discrete solution fields from the source mesh to adapted target mesh. For a proper coupling between the flow solver and the refinement module, the remapping procedure needs to be sufficiently accurate and conservative. Various conservative interpolation approaches have been proposed in the literature [13, 24, 22, 3, 20, 38]. In the current study, the conservative approach similar to that of Menon and Schmidt [38] is used to transfer cell-centered data.

[a]

[b]

[c]

Figure 3: Marked elements for refinement (red) [a], refined mesh with Schneiders algoritm [43] [b] and refined mesh with Ito et al. algorithm [27] [c].

6

Algorithm 1 Template based refinement algorithm for hexahedral elements. • Set element refinement levels according to a refinement criterion. • Repeat until there is no change in element levels: – For all elements marked for coarsening, apply from highest to lowest level: ∗ Check and adjust neighboring element levels so that there is at most one level difference. ∗ Find the maximum target level of all siblings of the element and set target level as the maximum level. ∗ If the element is IRREGULAR or marked for coarsening: · set PARENT as ACTIVE. · set CHILDREN of the PARENT as INACTIVE. · set target level of PARENT as the maximum target level of of its children. • Repeat until there is no change in element levels: – For all elements marked for refinement, apply from lowest to highest level: ∗ Check and adjust neighboring element levels so that there is at most one level difference. ∗ Convert arbitrary node mark configurations to one of the predefined template. • Propagate element levels to nodes. • Generate edge numbering information. • For all elements, apply from lowest to highest level: – Select appropriate template: ∗ Count marked nodes of element to decide the template. ∗ Check marked nodes to decide how to rotate the template. · Rotate template. – Apply template: ∗ ∗ ∗ ∗ ∗ ∗

Check FACEMARKS to see if there are previously created nodes on the element faces. Check EDGEMARKS for previously created nodes on the element edges. Create new nodes. Mark newly created nodes on the element faces. Mark newly created nodes on the element edges. Create new elements.

– Update tree structure: ∗ ∗ ∗ ∗ ∗ ∗

Set refined element as INACTIVE. Set newly created elements ACTIVE. Set newly created elements as REGULAR if it is fully refined. If not it is IRREGULAR. Set newly created elements as CHILDREN of the refined element. Set the refined element as PARENT to newly created elements. Update refinement LEVEL of elements.

7

Node marking

Template used

1

NA

2A

2B

T2

3A

T3

4A

T4

2C

3B

3C

4B

4C

4D

4E

5A

5B

6A

T6

4F

5C

6B

6C

7

8

T8

Table 1: Node markings and corresponding refinement templates used by Ito et al. [27].

The work of Menon and Schmidt [38] is based on the supermesh approach proposed by Farrell et al. [20]. In the supermesh approach, an intermediate mesh is constructed in order to facilitate the interpolation procedure with the following features: • Supermesh must contain all nodes and preserve the edges of both source and target meshes. • Intersection area of each cell in supermesh must be either zero or equal the area of supermesh cell. The supermesh acts as an intermediate mesh for solution transfer from the source mesh to adapted target mesh. Figure 4 illustrates the construction of a super mesh from quadrilateral source and target meshes. The supermesh construction in three-dimensions is also shown in Figure 5. In here, a hexahedral cell is divided into 24 tetrahedra. However, the construction of supermesh is not unique and it is proven to be rather challenging in three-dimensions using a constrained Delaunay triangulation algorithm [20]. Therefore, instead of constructing the supermesh, the current interpolation procedure employs local element by element intersection calculations between the source and target elements. For this purpose, the quadrilateral source and target elements are initially split into 4 triangular elements. Then the possible intersections for triangular source and target elements are identified. This identification step is straightforward since the present template based mesh refinement procedure is hierarchical in nature. For the intersection of a pair of triangular elements, each edge of the target triangle is considered for clipping the source triangular element and the remaining source element region is formed by a new set of triangular elements. Once the new triangular elements are constructed for the intersection region, a least square procedure is used for evaluating the gradient values within each quadrilateral source element. These values are stored to evaluate the values at the centroids of the new triangles. For the second order 8

=⇒

=⇒

[a]

[b]

=⇒

=⇒ [e]

[c] [d] Figure 4: Two-dimensional quadrilateral source mesh [a], triangulated source mesh [b], quadrilateral target mesh [c], triangulated target mesh [d] and supermesh construction [e]. accuracy, an interpolation method must be capable of preserving a linearly-varying function q(x). q(xj ) = qi + (xj − xi ) · (∇q)i + O(xj − xi )2

(4)

where xj is the centroid of new triangles, xi is the centroid of the quadrilateral source elements and (∇q)i is the gradient vector, which is assumed to be constant, as shown in Figure 6. The centroid values of the target quadrilateral elements are computed by summation of the new triangle values over their areas Aj qT =

1 ∑∑ Aj q(xj ) AT i j

(5)

where AT is the target quadrilateral element area. The second summation is due to the contribution from the neighbouring quadrilateral source elements. In the current interpolation algorithm, the source element values are directly copied if the element template does not change after the refinement process. Menon and Schmidt [38] demonstrated that the present second-order interpolation will exactly conserve the source and target integral values. Because the integral value QT over the target element can be written as ne ∑ ∑ QT = [qi + (xj − xi ) · (∇q)i ] Aj (6) i=1

j

If we assume that the gradient value is constant within the source quadrilateral element QT =

ne ∑

qi

i=1



Aj +

j

ne ∑ ∑ (∇q)i · (xj − xi )Aj i=1

(7)

j

∑ ∑ Noting that Ai = j Aj and j (xj − xi )Aj = 0 since xi is the quadrilateral source element centroid, the integral source value will be recovered. QT =

ne ∑ i=1

qi



Aj =

j

ne ∑

qi Ai = QS

(8)

i=1

In here, it should be noted that the interpolation is globally conservative only if the source and target mesh boundaries are conformal. This will be inherently satisfied with the current template based mesh refinement algorithm since the new vertices are introduced on the coarse mesh straight edges. The present approach is also locally conservative for the conformal refinement since the summation of CHILDREN 9

=⇒

=⇒

[a]

[b]

=⇒

=⇒ [e]

[c]

[d]

Figure 5: Three-dimensional hexahedral source mesh [a], tetrahedral source mesh [b], hexahedral target mesh [c], tetrahedral target mesh [d] and supermesh construction [e]. values for a refined element has to be equal to the PARENT element value in the source element. The above cell-centered algorithm is used to transfer the discrete solution fields involving the extra stress tensor components in the present viscoelastic simulations in section 3. In the current work, the method of Menon and Schmidt [38] has been extended for the side-centered data. The same idea is also applicable to vertex-based data. For this purpose, the quadrilateral element as shown in Figure 7 is split into four triangles and their second-order cell center values qci are approximated as follows [ ] 2 1 q1 + q2 + q3 + q4 qci = qi + (9) 3 3 4 Then, the problem turns into transferring the cell-centered data from the triangular source mesh to the triangular target mesh, which can be easily handled with the cell-centered remapping algorithm. However, the new cell-centered data has to be interpolated back to the side-centered data. Therefore, the right hand side of the equation (9) is implicitly constructed in a matrix form, which is actually the mass matrix in [47] for the momentum part, while the left hand side is summed over the triangular elements which belong to the same common edge qi . The resulting system of algebraic equations becomes [MT ]{qT } = RHS

(10)

where [MT ] is the mass matrix for the target mesh. Although we have to solve an algebraic equation for the side-centered data, the mass matrix is highly well conditioned and it converges only within several iterations. The details of the conservative remapping algorithm for the side-centered data is provided in Algorithm 2. The present side-centered algorithm is used to transfer the velocity vector components between the source and target meshes. It is important to note that for an overall conservative algorithm, both the interpolation algorithm and the flow solver should be conservative. Furthermore, the mass matrix used in the time integration scheme in the flow solver should be the same with the mass matrix used in the conservative interpolation. The second point is generally not well understood in the literature. In the current algorithm, the refinement and interpolation algorithms are implemented in 10

xi xj

Figure 6: Two-dimensional cell-centered conservative interpolation from quadrilateral source element centroid (xi ) to triangular supermesh centroid (xj ). parallel and computational mesh after each refinement is partitioned into sub-domains for a balanced domain decomposition using the METIS library [30]. q3 qc3

q4

qc4

qc2

q2

qc1 q1

Figure 7: Two-dimensional side-centered data interpolation from quadrilateral element edge midpoints to triangular element centroids.

Algorithm 2 Solution remapping algorithm for side-centered formulation • Decompose source mesh into simplicial elements (triangle/tetrahedra). • Decompose target mesh into simplicial elements (triangle/tetrahedra). • Compute cell-centered quantities qci for each source element. • Apply cell-centered remapping algorithm between the source and target meshes. • Construct mass matrix MT over the target mesh and construct RHS. • Solve the linear system to compute face-centered data on the target mesh.

11

3

COMPUTATIONAL RESULTS

In this section, the template based refinement algorithm and conservative remapping are initially tested for cell- and side-centered data. Then the proposed method is applied to an Oldroyd-B fluid past a confined circular cylinder in a rectangular channel and a sphere falling in a circular tube using several different adaptation criterions. The simulations indicate that the present mesh refinement algorithm allows to reproduce highly accurate numerical results.

3.1

Conservative Interpolation

In order to validate two-dimensional remapping algorithms for cell- and side-centered data, three different analytical functions including linear, exponential and sinusoidal terms are considered in a domain of [−1, 1] × [−1, 1]. • Linear field: q(x, y, z) = 2x + 3y + z • Exponential field: q(x, y, z) = e−( r2 ) where r = 1

√ x2 + y 2 + z 2

• Sinusoidal field: q(x, y, z) = 1 + sin(2πx)sin(2πy) These analytical values are defined on the initial mesh cell-centers and edge-midpoints. Then the values are transferred to the refined target mesh and vice versa. This process is repeated 100 times in our case. The computational initial and refined meshes are shown in Figure 8. The initial mesh consists of 441 nodes and 400 quadrilateral elements meanwhile the refined mesh consists of 1090 nodes and 1049 active quadrilateral elements. The difference in the global integral quantity evaluated on the initial source mesh and the final refined mesh is provided in Table 2. The results of cell-centered data show that the remapping of linear field is almost exact (up to machine precision) while exponential and sinusoidal fields show only very mild increase in the error. However, there is a slight increase in the difference magnitude for the side-centered data because the interpolation of the side-centered data requires to solve an algebraic equation and the present solver tolerance is set to |b − Ax|2 = ϵ where ϵ = 1 × 10−32 .

[a] [b] Figure 8: Two-dimensional initial source [a] and refined target [b] meshes used for 2D conservative remapping. In addition to two-dimensional cases, similar three-dimensional validation cases are also performed in a domain of [−1, 1] × [−1, 1] × [−1, 1]. The used computational meshes are shown in Figure 9. The initial mesh consists of 800 nodes and 6859 hexahedral elements, meanwhile the refined target mesh consists of 16652 nodes and 14966 active hexahedral elements. The remapping process is again repeated 100 times. The global integral errors are provided in Table 3. The magnitude of the errors is very similar to that of 12

Field function q(x, y) = 2x + 3y 1 q(x, y) = e−( r2 ) q(x, y) = 1 + sin(2πx)sin(2πy)

Global integral error for cell-centered data -3.552713678800E-015 3.885780586188E-016 5.995204332975E-015

Global integral error for side-centered data -6.661338147750E-015 1.9428902930940E-015 -1.0547118733938E-014

Table 2: Global integral errors for several different analytical functions after 100 remappings for 2D. two-dimensional cases. For the intersection tests in the present conservative interpolation, we not only considered the CHILDREN of common hexahedral elements at the same LEVEL in the source and target meshes but also considered all the CHILDREN of the common hexahedral element neighbours in the source mesh due to non-planar surfaces in three-dimensions. However, there is no need to consider the common element neighbours in two-dimensions. The comparison of initial value contours and the final interpolated value contours is provided in Figure 10 for the sinusoidal function. In here, a new polyhedral mesh is constructed for the face-centered data in order to present the function values at exactly face mid-points. The present results clearly indicate that the cell- and face-centered data remapping are conservative.

[a]

[b]

Figure 9: Three-dimensional initial source [a] and refined target [b] meshes used for 3D conservative remapping.

Field function q(x, y, z) = 2x + 3y + z 1 q(x, y, z) = e−( r2 ) q(x, y, z) = 1 + sin(2πx)sin(2πy)sin(2πz)

Global integral error for cell-centered data 5.551115123125E-017 -3.885780586188E-016 -2.220446049250E-016

Global integral error for face-centered data 2.973905572755E-015 4.607425552194E-015 3.552713678800E-015

Table 3: Global integral errors for several different analytical functions after 100 remappings for 3D.

13

[a]

[b]

Figure 10: Three-dimensional initial value [a] and interpolated value after 100 cycles [b] on y = 0 plane for side-centered data using q(x, y, z) = 1 + sin(2πx)sin(2πy)sin(2πz).

3.2

Oldroyd-B Fluid Past a Confined Circular Cylinder

The present parallel refinement algorithm is applied to the classical two-dimensional benchmark problem of an Oldroyd-B fluid flow past a confined circular cylinder in a channel. The problem has been extensively studied by many researchers [2, 1, 10, 12, 25, 26, 32, 39, 47, 18] in the literature. For this case, we consider a circular cylinder of radius R positioned symmetrically in a channel of width 2H. The blockage ratio of the computational domain R/H is taken to be 0.5 and the domain extends 12R in the upstream and downstream of the cylinder. The dimensionless parameters are the Reynolds number Re = ⟨U ⟩R/η, the Weissenberg number W e = λ⟨U ⟩/R and the viscosity ratio β = ηs /η. The physical parameters are the fluid density ρ, the average velocity at the inlet ⟨U ⟩, the relaxation time λ, the zero-shear-rate viscosity of the fluid η and the solvent viscosity ηs . The viscosity ratio β is set to 0.59, which is the value used in benchmarks for the Oldroyd-B model. The physical boundary conditions are the fully developed velocity profile at the inlet, the natural (traction-free) boundary condition at the outlet and the no-slip boundary conditions on the solid walls. The extra stress boundary conditions at the inlet are introduced through the fluxes using their analytical values. A solution based adaptive mesh refinement technique requires an adaptation criterion to determine elements for refinement in order to improve local mesh resolution. In order to capture the sharp extra stress tensor gradients within the viscoelastic wake along the channel symmetry line behind the cylinder and the viscoelastic boundary layer over the cylinder surface, an initial sensor function which is based on the von Mises stress value is defined as √ f1 (σ) = A [max(|λ1 |, |λ2 |)] (11) where A is quadrilateral element area and λi values are the eigenvalues of the Hessian matrix ] [ ] [ ∂2σ ∂2σ σxx σxy ∂x2 ∂x∂y H= = ∂2σ ∂2σ σxy σyy 2 ∂x∂y

(12)

∂y

and von Mises stress value σ is given as σ

=

1 √ 2 + T 2 + 6T √ (Txx − Tyy )2 + Txx xy yy 2

(13)

In order to evaluate the second-order derivations for the Hessian matrix, a third-order least square approach is constructed at the element centroid to determine the best values of first- and second-order derivatives in a least square sense. min||Ax − b||2 14

(14)

where 

σx  σy  x=  σxx  σyy σxy and A matrix is the N × 5 matrix  x1 − xi y1 − yi  x2 − xi y2 − yi   . . A=  . .   . . xN − xi yN − yi

   ,  

− xi )2 − xi )2 . . . 1 (x − xi )2 N 2 1 2 (x1 1 2 (x2

    b=   

σi − σ1 σi − σ2 . . . σi − σN

− yi )2 − yi )2 . . . 1 (y − yi )2 N 2 1 2 (y1 1 2 (y2

       

(x1 − xi )(y1 − yi ) (x2 − xi )(y2 − yi ) . . . (xN − xi )(yN − yi )

(15)

       

(16)

The number N includes its adjacent elements and their neighbours (first- and second-neighbours). LAPACK driver routines are used for both the solution of the least square equation problem and the computation of eigenvalues in equation (11). In here, the sensor function values f1 (σ) above a certain threshold value, which is 0.5, are marked for mesh refinement. The present sensor function is very similar to the work of Bo and Shashkov [50] in which the maximum absolute eigenvalue of the Hessian matrix is used. Finding a sensor function which is appropriate for the present benchmark problem is rather difficult. For example, the use of a rather coarse mesh requires a relatively small threshold value in order to capture the viscoelastic wake. However, this leads an excessive number of elements within the viscoelastic boundary layer. Therefore, the elements within a refinement box of [0, 8R] × [−0.1R, 0.1R] 1 are flagged for further refinement in the wake region and the maximum refinement LEVEL outside the refinement box is limited to 2. This is because the solution within the viscoelastic boundary layer does converge within the first several refinements and we are particularly interested in the convergence within the viscoelastic wake. In addition, the maximum refinement LEVEL within the refinement box is also limited to 5 because the time stepping algorithm which decouples the calculation of the polymeric stress from the evolution of the velocity and pressure fields requires excessively small time steps with further mesh refinements due to the Courant-Friedrichs-Lewy restriction. Hence it is rather difficult to reach the steady state solution within a reasonable time for further refinements. The present refinement algorithm is initiated only if the RMS value of the extra stress tensor is lower than 1 × 10−6 . Another problem with the present benchmark problem is the curved boundary on the cylinder surface. Although the present conformal refinement algorithm leads to a conservative interpolation, the edges on the initial mesh curved boundaries lead to corners or geometric singularities with further mesh refinement and the sensor function picks the regions around these points causing excessive number of elements with non-smooth solutions. In order to overcome this problem, the new vertices on the boundary has to be projected to the cylinder surface. Although this modification leads to the loss of global conservation property, the effect is only limited to the coarse elements on the cylinder surface due to the local conservation property. The present two-dimensional calculations are performed on the Maslak machine (Intel(R) Xeon(R) CPU E5-2683, 2.00GHz) using 24 nodes and the Karadeniz (Intel Xeon 5550) machine using 128 nodes at the National Center for High Performance Computing. The initial calculation is presented for the flow of an Oldroyd-B fluid past a confined circular cylinder in a channel at a Weissenberg number of 0.7. The time accurate simulations are carried out using the first-order accurate Euler implicit method in time and the global time step is reduced after each refinement process based on minimum element area. We also tried several other sensor functions for a better accuracy and our best result is obtained with the following sensor function √ f2 (σ) = A [max(|λ1 |, |λ2 |)] (17) 1 The bracketed numbers refer to distances along the x− and y−axis while assuming assuming that the cylinder center and its rear stagnation point are located at (0,0) and (R,0), respectively.

15

Refinement Level 0 1 2 3 4 5

Vertex # 10305 33025 175132 202818 407734 1877014

Element # 10025 35466 195145 226527 457841 2113179

Edge # 20330 65577 349238 404574 814380 3752920

Active Element # 10025 32552 174106 201756 406646 1875906

min(Ai ) 3.510E-4 3.743E-5 4.079E-6 5.292E-7 5.864e-8 6.498E-9

max(Txx ) 17.34 22.84 29.06 35.10 39.29 41.23

Table 4: Computational mesh details for each refinment level. min(Ai ) shows the minimum element area and max(Txx ) indicates the maximum value in the wake at W e = 0.7 and β = 0.59. with a threshold value of 5. The computational meshes for the above refinement criteria f2 (σ) are provided in Figure 11 along the computed Txx contours. The further details of the meshes for each level are also provided in Table 4. The final mesh has 1877014 vertices and 1875906 active quadrilateral elements with 15009464 DOF. The significant decrease in the minimum element area is evident from Table 4. The details of the final refined computational mesh are also provided in Figure 12 just behind the cylinder surface. The convergence of Txx with mesh refinement is provided in Figure 13 on the cylinder surface and along the center line in the wake and the numerical result for the final level is compared with the several results available in the literature. Although the mesh convergence trend is obvious at W e = 0.7 as seen in Figure 13-[a], there is some small difference for the maximum value of Txx along the center line in the wake. On the other hand, the initial refinement criteria f1 (σ) leads to a lower value by giving a maximum value of 39.97 in the wake for LEVEL 5. The details of the convergence of the present extreme values of Txx in the wake are provided in Table 4. Nevertheless, the comparison of the present final result indicates very good agreement with the numerical results of Fan et al. [18], Hulsen et al. [26] and Afonso et al. [1]. The numerical results of Fan et al. [18] were obtained using a Galerkin/least-squares hp finite element method. Hulsen et al. [26] used the DEVSS/DG formulation in a FEM context with the log-conformation formulation. Afonso et al. [1] employed a structured collocated FVM based on the log-conformation formulation on highly refined meshes. The comparison in Figure 13-[b] shows a very good agreement. Although the maximum value of Txx is lower on the cylinder surface for mesh M60W R in reference [1], the authors present favorable comparisons on the cylinder surface using higher mesh resolutions (see, for example, figure 15 of [2] for mesh M120). In addition, there is a remarkable agreement with the one-dimensional DG calculation of Hulsen et al. [26] in Figure 6 of [26] (M4-1D) along the center line in the wake. More recently, Gu´enette et al. [25] combined a finite element approach based on a Discrete Elastic Viscous Split Stress (DEVSS) formulation with an anisotropic adaptive remeshing strategy and confirmed the numerical results of Hulsen et al. [26] at W e = 0.7 which is in accord with the present results. The refinement criteria above leads to 13067686 active elements with a time step of 5.9974 × 10−5 for LEVEL 6. Therefore, it is prohibitively expensive to obtain the steady-state solution for this level. However, the refinement and remapping algorithms are sufficiently fast. As an example, the refinement and remapping steps for LEVEL 6 requires 14.79 and 43.59 seconds, respectively, on the Maslak machine at the National Center for High Performance Computing using 24 nodes. The second calculation is presented at a Weissenberg number of 0.8 with the same refinement criterion. The convergence of Txx with mesh refinement is illustrated in Figure 14-[a] on the cylinder surface and along the center line in the wake. The final mesh refinement leads to 2246211 vertices and 2245065 active quadrilateral elements with 17962812 DOF. For this case, no mesh convergence tendency is observed in the stress profiles. In addition to the above simulation, we reduced the sensor function threshold value by half and applied a prescribed refinement within the refinement box along the symmetry axis behind the cylinder in order to increase the mesh resolution close to the symmetry axis. The convergence of Txx profile with mesh refinement is provided in Figure 14-[b]. This case leads to 407213 vertices and 406134 active quadrilateral elements with 3251230 DOF. Although this modification reduces the number of quadrilateral elements significantly, we could not observe any significant change in the stress profile. Therefore, the convergence is very sensitive to the mesh resolution just next to the symmetry line in the wake. Even though further mesh refinement is possible, it becomes prohibitively expensive to obtain the converged steady state solution due to excessively small time step. In the future, we will implement the fully coupled version of the viscoelastic solver in order to avoid the Courant-Friedrichs-Lewy restriction 16

level = 0

level = 1

level = 2

level = 3

level = 4

level = 5 Figure 11: Adaptive mesh refinement (left) and computed Txx contours (right) at W e = 0.7 and β = 0.59. and apply the present template based refinement algorithm to higher Weissenberg numbers. Nevertheless, the adaptive local mesh refinement significantly improves the employed iterative solver solve time since the method uses the interpolated previous coarse mesh solution as an initial guess rather than zero solution. 17

Figure 12: Details of final adapted mesh (level = 5) just behind the cylinder at W e = 0.7 and β = 0.59. 120

120 Level=5 Level=4 Level=3 Level=2 Level=1 Level=0

100

80

+++ + + + + + + + + + + + + + + + + +

100

80

+

+ +

+

60

40

+

+

Txx

Txx

60

+

+

+ +

+

+

40

+ + +++ + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + ++ + + ++ + +++ ++++ ++ +++++ +++ +++++++

+ + +

20

20

0

0

-20 -1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

+ + + + + + + ++ +++ +++++++++++++++++++++++++++++++++++

-20 -1.5

3

Level=5 (Present) Yurun et al. Hulsen et al. Afonso et al.

+

-1

-0.5

0

0.5

x/R

1

1.5

2

+ + + + + +

2.5

3

x/R

140

140

120

120

100

100

80

80

60

Txx

Txx

[a] Figure 13: Convergence of Txx with mesh refinement level [a] and the comparison of Txx with several results in the literature [b] at W e = 0.7 and β = 0.59.

Level=5 Level=4 Level=3 Level=2 Level=1 Level=0

40

20

Level=5 Level=4 Level=3 Level=2 Level=1 Level=0

40

20

0

-20 -1.5

60

[b]

0

-1

-0.5

0

0.5

1

x/R

Figure 14: Convergence of Txx W e = 0.8 and β = 0.59.

1.5

2

2.5

-20 -1.5

3

-1

-0.5

0

0.5

1

x/R

1.5

2

2.5

[a] with adaptive mesh refinement [a] and prescribed mesh refinement [b] at

18

3

[b]

3.3

Oldroyd-B Fluid around a Sphere Falling in a Cylindrical Tube

The problem of a viscoelastic fluid around a sphere falling in a cylindrical tube has been studied extensively by various researchers [34, 6, 40, 17, 9, 15, 16]. For this problem, a rigid sphere of radius RS is falling with a terminal velocity Us along the axis of a cylindrical tube of radius RT . The ratio of the sphere radius to the tube radius is taken to be 0.5 and the downstream and upstream distance from the center of the sphere are 12RS . The dimensionless parameters are the Reynolds number Re = US RS /η, the Weissenberg number W e = λUS /RS and the viscosity ratio β = ηs /η. The physical parameters are the density ρ, the terminal velocity US , the relaxation time λ, the zero-shear-rate viscosity of the fluid η and the solvent viscosity ηs . The viscosity ratio β is chosen to be 0.50, which is the classical value used in the benchmark problem for the Oldroyd-B fluid. The present three-dimensional boundary conditions are the prescribed unit uniform velocity along the tube axis at the inflow and on the circular tube solid wall, no-slip boundary condition on the sphere surface and natural (traction-free) boundary conditions at the outflow. The initial mesh consists of 234620 nodes and 226181 hexahedral element leading to 3643853 DOF. The initial all-hexahedral mesh is created by using the mapping, paving and sweeping algorithms available within the CUBIT mesh generation environment [5]. The three-dimensional mesh refinement calculation is carried out for an Oldroyd-B fluid around a sphere falling in a cylindrical tube at a Weissenberg number of 1.2. Our best result is obtained for the following sensor function in three-dimensions √ f (σ) = V [max(|λ1 |, |λ2 |, |λ3 |)] (18) which is similar to that of the two-dimensional case but it uses three-dimensional von Mises stress and Hessian matrix. In here, V represents the element volume and the threshold value is also equal to 5. In three-dimensions, the number of elements increases rapidly with mesh refinement due to 1:3 refinement of edges. Even though the process can be highly parallelized, the excessive number of elements makes the approach impractical because of computational resource restraints. Therefore, the maximum refinement LEVEL is limited to 5 in the wake region and 1 in the rest of the domain since the solution around the sphere converges rapidly with mesh refinement. In this case only the hexahedral elements exactly on the symmetry axis behind the sphere are marked for further refinement in order to reduce the required number of elements as in the final two-dimensional case at W e = 0.8. The calculations are carried on the Levrek (Xeon E5-2690 2.90GHz) machine at TUBITAK ULAKBIM, High Performance and Grid Computing Center, using 96 nodes and the mesh refinement is proceeded to next LEVEL when the RMS values of the extra stress tensor is lower than 1 × 10−6 . The computational meshes for each refinement level are provided in Figure 15 along with the computed Txx contours for the z ≤ 0 domain. The further details of the refined meshes are provided in Table 5. The final refined mesh leads to 2340392 vertices and 2314365 active hexahedral elements with 37095687 DOF. For the present three-dimensional case, the second refinement leads to the refinement of most parts of the domain behind the sphere and significantly increases the number of hexahedral elements. The details of the final adapted mesh are also shown in Figure 16 close to the symmetry line behind the sphere. The convergence of the computed Txx values on the sphere surface and along the center line in the wake is provided in Figure 17 for each level and the result for the final level is compared with the axisymmetric simulations of Fan et al. [16] and Chauviere and Owens [9]. The results of Fan et al. [16] were obtained using a Galerkin/least-square hp finite element method meanwhile Chauviere and Owens [9] used an adaptive spectral element method based on a posteriori error indicator. The comparison shows a relatively good agreement keeping in mind that Refinement Level 0 1 2 3 4

Vertex # 234620 897878 1921140 2007012 2340392

Element # 226181 903899 1967604 2058899 2412123

Face # 686862 2650402 5718306 5973144 6965044

Active Element # 226181 876395 1898847 1983775 2314365

min(Vi ) 1.893E-6 8.406E-8 9.232E-9 3.370E-10 1.242E-11

max(Txx ) 25.23 35.77 44.29 51.51 54.88

Table 5: Computational mesh details with refinement level. min(Vi ) shows the minimum element volume and max(Txx ) indicates the maximum value in the wake at W e = 1.2 and β = 0.5. 19

level = 0

level = 1

level = 2

level = 3

level = 4 Figure 15: Adaptive mesh refinement (left) and computed Txx contours for z ≤ 0 domain (right) at W e = 1.2 and β = 0.5. the present simulation is three-dimensional. The present three-dimensional refinement and remapping algorithms requires 54.45 and 716.81 seconds, respectively, for LEVEL 4 on the Levrek machine at TUBITAK ULAKBIM using 96 nodes. Although the conservative remapping algorithm requires relatively longer time in three-dimensions, this is mainly due to the large number of possible tetrahedral element volume-volume intersection calculations between source and target meshes.

20

[a]

[b]

Figure 16: Partial view of final adapted mesh (level = 4) around a sphere [a] and enlarged view of mesh close to symmetry axis [b] at W e = 1.2 and β = 0.5. 80

60

60

40

40

Txx

Txx

80

20

20

Level=4 (Present) Chauvière&Owens Yurun

Level=4 Level=3 Level=2 Level=1 Level=0

0

-20 -1.5

-1

-0.5

0

0.5

1

x/R

1.5

2

2.5

0

-20 -1.5

3

[a]

-1

-0.5

0

0.5

1

1.5

2

2.5

x/R

Figure 17: Convergence of Txx with mesh refinement level [a] and the comparison of Txx with several results in the literature [b] at W e = 1.2 and β = 0.5.

4

CONCLUSIONS

A framework for parallel adaptive mesh refinement has been described for the large-scale simulation of viscoelastic fluid flows. The recursive dynamic mesh adaptation algorithm performs conformal allquadrilateral/hexahedral refinement and uses the templates given in the works of Schneiders [38] and Ito et al. [23] for two- and three-dimensions, respectively. A conservative second-order solution remapping technique similar to that of Menon and Schmidt [38] is used to transfer the cell-centered data and the approach has been extended for the side-centered data for the finite volume formulation. To the authors’ best knowledge, this work is the first application of the parallel conformal template based refinement algorithm for viscoelastic fluid flow problems. Although the proposed refinement algorithm produces sufficiently accurate numerical results for the present benchmark problems, the employed time stepping algorithm, which decouples the calculation of the polymeric stress from the evolution of the velocity and pressure fields, requires very small time steps at higher Weissenberg numbers with further mesh refinement while degrading the overall efficiency of the proposed method. In the future, we will implement the fully coupled version of the present solver with the Newton’s method as in the recent work of Damanika et al. [11] in order to apply the present approach to higher Weissenberg numbers.

21

3

[b]

5

ACKNOWLEDGMENTS

The authors gratefully acknowledge the use of the Ferrel machine at the Faculty of Aeronautics and Astronautics at ITU, the computing resources provided by the National Center for High Performance Computing of Turkey (UYBHM) under Grant No. 10752009 and the computing facilities at TUBITAK ULAKBIM, High Performance and Grid Computing Center.

References [1] A. Afonso, P.J. Oliveira, F.T. Pinho, M.A. Alves, The log-conformation tensor approach in the finite-volume method framework, J. Non-Newton. Fluid Mech. 157, (2009), 55–65. [2] M.A. Alves, F.T. Pinho, P.J. Oliveira, The flow of viscoelastic fluids past a cylinder: finite-volume high-resolution methods, J. Non-Newton. Fluid Mech. 97, (2001), 207–232. [3] B. Azarenok, A method for conservative remapping on hexahedral meshes. Mathematical Models and Computer Simulations, 1, (2009), 51–63. [4] S. Balay, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith and H. Zhang, PETSc Users Manual. ANL-95/11, Mathematic and Computer Science Division, Argonne National Laboratory, (2004). http://www-unix.mcs.anl.gov/petsc/petsc-2/ [5] T. D. Blacker, S. Benzley, S. Jankovich, R. Kerr, J. Kraftcheck, R. Kerr, P. Knupp, R. Leland, D. Melander, R. Meyers, S. Mitchell, J. Shepard, T. Tautges and D. White, CUBIT Mesh Generation Enviroment Users Manual: Volume 1, Sandia National Laboratories, Albuquerque, NM, (1999). [6] C. Bodart and M. J. Crochet, The time-dependent flow of a viscoelastic fluid around a sphere. J. Non-Newtonian Fluid Mech. 54, (1994), 303–329. [7] Z. Cai and C. R. Westphal, An adaptive mixed least-squares finite element method for viscoelastic fluids of Oldroyd type, J. Non-Newton. Fluid Mech. 159, (2009), 72–80. [8] A. E. Caola, Y. L. Joo, R. C. Armstrong and R. A. Brown, Highly parallel time integration of viscoelastic flows. J. Non-Newton. Fluid Mech. 100, (2001), 191–216. [9] C. Chauviere and R.G. Owens, How accurate is your solution? Error indicators for viscoelastic flow calculations, J. Non-Newton. Fluid Mech. 95, (2000), 1–33. [10] O.M. Coronado, D. Arora, M. Behr, M. Pasquali, A simple method for simulating general viscoelastic fluid flows with an alternate log-conformation formulation, J. Non-Newton. Fluid Mech. 147, (2007), 189–199. [11] H. Damanika, J. Hronb, A. Ouazzia and S. Turek, A monolithic FEM approach for the logconformation reformulation (LCR) of viscoelastic flow problems. J. Non-Newton. Fluid Mech. 165, (2010), 1105–1113. [12] H.-S. Dou, N. Phan-Thien, Parallelisation of an unstructured finite volume code with PVM: viscoelastic flow around a cylinder, J. Non-Newton. Fluid Mech. 77, (1998), 21–51. [13] J. K. Dukowicz and N. Padial, REMAP3D: A conservative three dimensional remapping code. Tech. Rep., LA-12136-MS, Report of Los Alamos National Laboratory, Los Alamos, NM, USA, (1991) [14] M. S. Ebeida, A. Patney, J. D. Owens and E. Mestreau, Isotropic conforming refinement of quadrilateral and hexahedral meshes using two-refinement templates. Int. J. Numer. Meth. Eng. 88, (2011), 974–985. [15] Y. Fan, Solution behavior of the falling sphere problem in viscoelastic flows. Acta Mechanica Sinica 19, (2003), 394–408.

22

[16] Y. Fan, Limiting behavior of the solutions of a falling sphere in a tube filled with viscoelastic fluids. J. Non-Newton. Fluid Mech. 110, (2003), 77–102. [17] Y. Fan, M. J. Crochet, High-order finite element methods for steady viscoelastic flows. J. NonNewton. Fluid Mech. 57, (1995), 283–311. [18] Y. Fan, R.I. Tanner, N. Phan-Thien, Galerkin/least-square finite element methods for steady viscoelastic flows. J. Non-Newton. Fluid Mech. 84, (1999), 233–256. [19] C. Farhat, M. Lesoinne and P. L. Tallec, Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: momentum and energy conservation, optimal discretization and application to aeroelasticity. Comput. Methods Appl. Mech. Engrg. 157, (1998), 95–114. [20] P. Farrell, M. Piggott, C. Pain, G. Gorman and C. Wilson, Conservative interpolation between unstructured meshes via supermesh construction. Comput. Methods Appl. Mech. Eng. 198, (2009), 33–36. [21] R. Fattal and R. Kupferman, Constitutive laws for the matrix-logarithm of the conformation tensor. J. Non-Newton. Fluid Mech. 123, (2004), 281–285. [22] R. Garimella, M. Kucharik and M. Shashkov , An efficient linearity and bound preserving conservative interpolation (remapping) on polyhedral meshes. Comput. Fluids, 36, (2007), 224–237. [23] R. Garimella, Conformal renement of unstructured quadrilateral meshes. Proceedings of the 18th International Meshing Roundtable, (2009), 31–44. [24] J. Grandy, Conservative remapping and region overlays by intersecting arbitrary polyhedra. J. Comput. Physics, 148, (1999), 433–466. [25] R. Gu´enette, A. Fortin, A. Kane and J.-F. H´etu, An adaptive remeshing strategy for viscoelastic fluid flow simulations. J. Non-Newton. Fluid Mech. 153, (2008), 34–45. [26] M.A. Hulsen, R. Fattal, R. Kupferman, Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms. J. Non-Newton. Fluid Mech. 127, (2005), 27–39. [27] Y. Ito, A. M. Shih and B. K. Soni, Octree-based reasonable quality hexahedral mesh generation using a new set of refinement templates. Int. J. Numer. Meth. Eng. 77, (2009), 1809–1833. [28] R. K. Jaiman, X. Jiao, P. H. Geubelle and E. Loth, Conservative load transfer along curved fluid-solid interface with non-matching meshes. J. Comput. Physics 218, (2006), 372–397. [29] H. Jin and R. I. Tanner, An a posteriori error estimate and adaptive refinement procedure for viscoelastic fluid flows: The modified upper-convected Maxwell model. Comput. Mech. 13, (1994), 400–413. [30] G. Karypis and V. Kumar, A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM J. Sci. Comput. 20, (1998), 359–392. [31] R. Keunings, On the high Weissenberg number problem. J. Non-Newton. Fluid Mech. 20, (1986), 209–226. [32] J.M. Kim, C. Kim, K.H. Ahn, S.J. Lee, An efficient iterative solver and highresolution computations of the Oldroyd-B fluid flow past a confined cylinder. J. Non-Newton. Fluid Mech. 123, (2004), 161–173. [33] D. Kwak and Y. Im, Hexahedral mesh generation for remeshing in three-dimensional metal forming analyses. J. Mater. Process. Tech. 138, (2003), 531–537.

23

[34] W. J. Lunsmann, L. Genieser, R. C. Armstrong, R. A. Brown, Finite element analysis of steady viscoelastic flow around a sphere in a tube: calculations with constant viscosity models. Int. J. Numer. Meth. Eng. 48, (1993), 63–99. [35] D. J. Mavriplis, Unstructured mesh generation and adaptivity. ICASE Report No. 95–26, (1995). [36] L. G. Margolina and Mikhail Shashkov, Second-order sign-preserving conservative interpolation (remapping) on general grids. J. Comput. Physics 184, (2003), 266–298. [37] I. Mutlu, P. Townsend and M. F. Webster, Adaptive solutions for viscoelastic flows, Comm. Numer. Methods Eng. 12, (1996), 643–655. [38] S. Menon and D. P. Schmidt, Supermesh construction for conservative interpolation on unstructured meshes: An extension to cell–centered finite-volume variables. Comput. Methods Appl. Mech. Eng., 200, (2011), 2797–2804. [39] R.G. Owens, C. Chauvire, T.N. Philips, A locally-upwinded spectral technique (LUST) for viscoelastic flows. J. Non-Newton. Fluid Mech. 108, (2002) 49–71. [40] R. G. Owens and T. N. Phillips, Steady viscoelastic flow past a sphere using spectral elements. Int. J. Numer. Meth. Eng. 39, (1996), 1517–1534. [41] M. Parrish, M. Borden, M. Staten and S. Benzley, A selective approach to conformal refinement of unstructured hexahedral finite element meshes. Proceedings of the 16th International Meshing Roundtable, (2008), 251–268. [42] N. Roquet and P. Saramito, An adaptive finite element method for Bingham fluid flows around a cylinder. Comput. Methods Appl. Mech. Eng. 192, (2003), 3317–3341. [43] R. Schneiders, Refining quadrilateral and hexahedral element meshes. 5th International Meshing Roundtable, (1996), 679–688. [44] L. Sun, G. Zhao and X. Ma, Adaptive generation and local refinement methods of three-dimensional hexahedral element mesh. Finite Elem. Anal. Des. 50, (2011), 184–200. [45] Y. Saad, A flexible inner-product preconditioned GMRES algorithm. SIAM J. Sci. Statist. Comput. 14, (1993), 461–469. [46] M. Sahin, A preconditioned semi-staggered dilation-free finite volume method for the incompressible Navier-Stokes equations on all-hexahedral elements. Int. J. Numer. Meth. Fluids. 49, (2005), 959– 974. [47] M. Sahin, A stable unstructured finite volume method for parallel large-scale viscoelastic fluid flow calculations. J. Non-Newton. Fluid Mech. 166, (2011), 779–791. [48] M. Sahin, Parallel large-scale numerical simulation of purely-elastic instabilities behind a confined circular cylinder in a rectangular channel. J. Non-Newton. Fluid Mech. 195, (2013), 46–56. [49] I. E. Sutherland and G. W. Hodgman, Reentrant polygon clipping, Commun. ACM, 17, (1974), 32–42. [50] W. Bo and M. Shashkov, R-adaptive reconnection-based arbitrary Lagrangian Eulerian method − R-ReALE. J. Math. Study 48, (2015), 125–167 [51] H. Zhang and G. Zhao, Adaptive hexahedral mesh generation based on local domain curvature and thickness using a modified grid-based method. Finite Elem. Anal. Des. 43, (2007), 43–691. [52] Y. Zhang, X. Liang and G. Xu, A robust 2-refinement algorithm in octree or rhombic dodecahedral tree based all-hexahedral mesh generation. Comput. Methods Appl. Mech. Eng. 256, (2013), 88-100.

24