key: cord-0058702-824ciypn
authors: Azevedo, Ramoni Z. S.; Santos, Isaac P.
title: Multiscale Finite Element Formulation for the 3D Diffusion-Convection Equation
date: 2020-08-20
journal: Computational Science and Its Applications - ICCSA 2020
DOI: 10.1007/978-3-030-58808-3_33
sha: 095eff0f4aba562f24370b93e0c7f6f93af1760e
doc_id: 58702
cord_uid: 824ciypn

In this work we present a numerical study of the Dynamic Diffusion (DD) method for solving diffusion-convection problems with dominant convection on three dimensional domains. The DD method is a two-scale nonlinear model for convection-dominated transport problems, obtained by incorporating to the multiscale formulation a nonlinear dissipative operator acting isotropically in both scales of the discretization. The standard finite element spaces are enriched with bubble functions in order to add stability properties to the numerical model. The DD method for two dimensional problems results in good solutions compared to other known methods. We analyze the impact of this methodology on three dimensional domains comparing its numerical results with those obtained using the Consistent Approximate Upwind (CAU) method.

The diffusive-convective transport equation plays an important role in numerous physical phenomena. It models several problems in engineering, such as the process of pollution on rivers, seas and the atmosphere; oil reservoir flow; heat transfer, and so forth [1] . However, the accurate and reliable numerical simulation of this problem is still a challenging task [2] [3] [4] .

It is well known that the standard Galerkin Finite Element method (FEM) is not appropriate for solving convection-dominated problems [2] , since it experiences spurious oscillations in the numerical solution, that grow as the diffusion coefficient tends to zero. The problem is caused by its loss of stability with respect to the standard H 1 -norm [5] . Since the early 1980s, various stabilization techniques have been developed to overcome the weaknesses of the FEM to solve such problems.

One of the most known stabilized formulation to solve convection-dominated transport problem is the Streamline Upwind Petrov-Galerkin (SUPG) method (also known as the Streamline Diffusion Finite Element method (SDFEM)), proposed in [6] and analyzed in [7] . This methodology consists of adding residualbased terms to the FEM, preserving the consistency of the variational formulation. Other methodologies were developed to circumvent the difficulty of the FEM to solve problems of this nature, such as, the Galerkin-Least Square method (GLS) [8] (a variant of the SUPG), edge-based stabilization (the Continuous Interior Penalty (CIP) method [9] ), and local projections (the Local Projection Stabilization (LPS) [4, 10] ). All these numerical formulations are capable to stabilize most of the unphysical oscillations from the FEM, caused by dominating convection. However, localized oscillations may remain in the neighborhood of sharp layers for nonregular solutions [11] . As a remedy, a number of so-called Spurious Oscillations at Layers Diminishing (SOLD) schemes have been proposed; see the excellent reviews [12, 13] . SOLD schemes add additional, in general nonlinear, discontinuity (or shock-capturing) terms to these formulations in order to reduce the localized oscillations.

Residual-based stabilized methods have been reinterpreted within the context of the Variational Multiscale (VMS) framework by Hughes and collaborators in the 1990s [14] [15] [16] [17] . The starting point of the VMS framework is the splitting of the finite element spaces into two parts: resolved (or coarse) and fine (or subgrid) subspaces. Consequently, the variational problem is decomposed into two subproblems, one associated to the resolved scales and the other associated to the fine ones. The fine scales are approximated or modeled, and then, incorporeted into the coarse scale problem, yielding in an enriched problem for the resolved scales. Many classical stabilization formulation for convection-dominated problems, such as SUPG and GLS can be reinterpreted as residual-based fine-scale models [16, 18] . Also, other known multiscale methods, such as the residual-free bubbles method [19, 20] and subgrid viscosity methods [21, 22] , can be placed within this framework.

Over the last years, some nonlinear two-scale variational methods have been developed, such as the Nonlinear Subgrid Stabilization (NSGS) [23, 24] and Dynamic Diffusion (DD) methods [11, 25] to solve convection-dominated problems; the Nonlinear Multiscale Viscosity (NMV) method [26, 27] to solve the system of compressible Euler equations. In [28] the NSGS method was applied to solve incompressible Navier-Stokes equations. In this class of methods, bubble functions are used to enrich the standard finite element spaces and residualbased nonlinear artificial diffusion terms are added to the resulting formulation, in order to improve the stability properties of the numerical model. In particular, the DD method is obtained by incorporating to the VMS formulation, a nonlinear dissipative operator acting isotropically in both scales of the discretization. This method results in good solutions for two dimensional problems, compared to other known stabilized methods, as can be seen in [11, 25] .

The aim of this work is to present a numerical study of the DD method for solving steady-state diffusion-convection problems with dominant convection on three dimensional domains. We compare its computational results with the SOLD method, Consistent Approximate Upwind Petrov-Galerkin (CAU) [29] , in the solution of two convection dominated problems with internal and external layers.

The remainder of this work is organized as follows. Section 2 briefly addresses the governing equation, the numerical formulations of the Galerkin finite element, CAU and DD methods. The numerical experiments are conducted in Sect. 3 and Sect. 4 concludes this paper.

Let Ω ⊂ R 3 be a bounded domain with a Lipschitz boundary Γ with an outward unit normal n. The linear and steady-state diffusion-convection equation, with a Dirichlet boundary condition, is modeled by

where u represents the quantity being transported, κ > 0 is the diffusivity coefficient, β ∈ [L ∞ (Ω)] 3 is the velocity field such that ∇ · β = 0, f ∈ L 2 (Ω) is a given source term and g ∈ L 2 (Γ ) is a given function.

Consider the set of all admissible functions S and the space of test functions V defined as

The continuous variational problem associated to the boundary value problem (1) and (2) can be stated as: find u ∈ S such that

with

We consider a regular partition T h = {Ω 1 , Ω 2 , . . . , Ω nel } of the domain Ω into tetrahedral elements. The discrete spaces associated to S and V are defined as

in which P 1 (Ω e ) is the set of first order polynomials on local coordinates. The Galerkin formulation is obtained by restricting the problem (5) to the discrete spaces S h ⊂ S and V h ⊂ V [2, 30] . In other words, the variational statement of the approximate problem consists of finding u h ∈ S h , such that,

where

The Galerkin finite element approximation for the variable u h is given at the linear tetrahedral element level as follows,

where N j is the local interpolation function associated with the nodal point j of the element Ω e . These functions expressed in global variables x = (x, y, z) are of the form

where |Ω e | is the volume of the element, calculated by

whereas the coefficients of the shape functions are given by

The transformation from global variables x = (x, y, z) to local variables ξ = (ξ, η, ζ), represented in Fig. 1 , results in the following interpolation functions [31] , 

and the inverse transformation matrix is determined by 

The integrals present in the formulations can be easily solved analytically, using the formula for integrating products of powers of linear interpolation functions over tetrahedral elements [31] ,

The FEM method results in stable and accurate approximate solutions when the transport equation is diffusive dominant [6, 12] . However, when the problem is convective dominant, the numerical solution may presents spurious oscillations, mainly when boundary layers (narrow regions with steep gradients) appear in the solution, requiring the use of stabilized methods [6, 12, 23, 32] . In the following subsections we describe two different stabilization methods used to control the oscillations caused by the dominant convection.

The CAU is a stabilized method proposed in [29] for solving convectiondominated problems. It is a SOLD (discontinuity-capturing) model that consists of splitting the stabilization formulation into two parts, where the first is a linear stabilization term, the standard SUPG operator [6] , and the second one is a nonlinear operator aiming to control the solution gradient, preventing localized oscillations around boundary layers. In this work we follow the CAU formulation described in [33] . The CAU formulation consists of finding u h ∈ S h such that,

is given by (11) , and

The stabilization parameters τ SUP G and τ CAU are defined as follow,

where · is the Euclidian norm in R 3 , P e is the local Péclet number and h is the characteristic length of mesh:

on Ω e ,P e is the local Péclet number andh the caracteristic length of mesh modifieds in termos of an auxiliary vector fieldβ,

The velocity fieldβ, given bȳ

determines a new upwind direction, β −β, calculated in each iteration in terms of the approximate solution of the previous iteration. It is important to note that the β −β is in the direction of ∇u h , allowing the CAU method to possess greater stability than the SUPG method in regions near to internal and/or external boundary layers [33] . The method is solved by using an iterative procedure defined as: given u n h , we find u n+1 h satisfying,

where the initial solution u 0 h is obtained by the standard Galerkin finite element method (10) . The nonlinear convergence is checked under a given prescribed relative error with tolerance tol CAU and a maximum number of iterations (itmax CAU ).

As discussed in the introduction, the DD method is a nonlinear multiscale method that consists of adding a nonlinear artificial diffusion operator acting on both scales of the discretization [11, 25] . To present the DD method, we define the enriched spaces S E and V E , spaces of admissible functions and test functions, as

where S B is the space spanned by bubble functions. By denoting N b ∈ H 1 0 (Ω e ) as the bubble function defined in each element Ω e , we define S B | Ωe = span(N b ) and S B = ⊕| Ωe S B | Ωe for all Ω e in T h . Here, we use the following function [34] ,

where N j is the local interpolation function associated with the nodal point j of the element Ω e .

The discrete variational formulation of the DD method can be statement as:

with , v) and F G (v) are given by (11) and (12) , respectivelly. The nonlinear operator A DD (u h ; u, v) in (21) is defined as

where τ DD (u h ) = μ(h)C b (u h ) represents the amount of artificial diffusion, with

where R(u h ) = β · ∇u h − f on Ω e , and the caracteristic length of mesh,

In the last expression, |Ω e | is the volume of the element Ω e , Γ + = {x ∈ Ω e ; β · n > 0} is the outflow part of Γ , and n is the unit outward normal vector to the boundary Γ . The method is solved by using an iterative procedure defined as: given u n , we find u n+1 satisfying

where the initial solution is u 0 h = 0. To improve convergence, the following weighting rule to determine (23) is used with w = 0.5:

Using the multiscale decomposition process, the Eq. (21) results in two equations (subproblems) in each element Ω e ∈ T h ,

where we have used the known results

. The Eqs. (24) and (25) result in the following local system of algebraic equations, associated with each element Ω e ,

where the local matrices and vectors in Eq. (26) are defined as follow,

and U e h and U e b are the vectors that store the variables u h | Ωe and u b , respectively. Performing a static condensation to eliminate the unknowns U e b at each element, the system (26) can be written in terms of the macro solution U e h , as follows,

Remark 1. Since κ > 0, the local matrix A e bb , associated to the term

is always nonsingular.

The global system is obtained by assembling all local systems (27) calculated on each element Ω e ,

In each iteration of (28)-(29), the linear system resulting is solved by the Generalized Minimal Residual (GMRES) method [35] . The nonlinear convergence is checked under a given prescribed relative error,

and a maximum number of iterations (itmax DD ).

In this section we present two numerical experiments in order to evaluate the quality of the numerical solution obtained with the DD method in three dimensional domains. The computational code was developed in C language in a computer with the following characteristics: Intel Corel i5-2450M, CPU 2.50 GHz × 4 and 8 GB of RAM. The linear systems are solved with the GMRES method considering 150 Krylov vectors in the base, the maximum number of 100 cycles and a tolerance of 10 −7 . For the nonlinearity of the CAU and DD methods, we use the following tolerance parameters: tol CAU = tol DD = 10 −2 ; and the following maximum number of iterations: itermax CAU = itermaxDD = 50.

This problem simulates a three dimensional convection-dominated diffusionconvection problem, with internal and external layers, in a unit cube Ω =]0, 1[ 3 with the following coefficients: κ = 10 −9 , β = (0, 1, 0) T and f = 1. The Dirichlet boundary conditions are defined as u = 0 on all faces except in the plane y = 0, where a jump discontinuity is introduced in the inflow boundary,

The domain Ω was partitioned in 32 divisions in each direction, so the mesh has 35937 nodes and 196608 elements. Due to the difficulty in representing results in 3D domains, the solutions are exhibited at the xy plane located at z = 0.5. Figure 2 shows the solutions obtained by the CAU and DD methods. The CAU solution ( Fig. 2(a) ) appears more diffusive in the external boundary layers (narrow regions close to Γ where ∇u is large) at x = 0, x = 1 and y = 1, whereas the DD method ( Fig. 2(b) ) presents some nonphysical under/overshoots in the external layer, at y = 1, mainly for x > 0.5. In the internal layer, both methods present some diffusion in the region next to the outflow boundary. For a better understanding of the behavior of the CAU and DD methods, we made three cuts (cross sections) in the solutions described in Fig. 2 , at x = 0.25, x = 0.75 and y = 0.9. Figure 3(a) shows the cross section at x = 0.25 (lower graph) and at x = 0.75 (upper graph). We can observe the diffusive behavior of the CAU method in the external boundary layer (at y = 1). In this region the DD solution is better, despite the slight under/overshoots. Figure 3(b) shows the cross section at y = 0.9. The CAU method is more diffusive in the three layers at x = 0, x = 0.5 and x = 1. 

In this experiment we consider a convection-dominated diffusion-convection problem defined on Ω =]0, 1[ 3 with the following coefficients: κ = 10 −6 , β(1, 1, 1) T and f = 1. The Dirichlet boundary conditions are defined as u = 0 on faces x = 0, y = 0 and z = 0 and u = 1 on faces x = 1, y = 1 and z = 1. The mesh used in this experiment is the same one as the Problem 1.

The solutions are exposed at the xy plane located at z = 0.5. Figure 4 shows the solutions obtained by the CAU and DD methods. The solution obtained with the CAU method appears to be more diffusive at the external boundary layer at x = 1 and at y = 1 than that of Fig. 4(b) . However, some non physical under/overshoots localized at the external layer remain in the DD solution. These numerical solutions are better understood, analyzing the cuts (cross section) of them shown in the Fig. 5, Fig. 5(a) shows the cross section at x = 0.7 and Fig. 5(b) shows the cross section at y = 0.9, and Fig. 6 shows the cross section along the diagonal of the xy plane. All of them show the more diffusive behavior of the CAU method in the boundary layer regions and the small oscillations presented by the DD method solution. 

This work deals with the evaluation of the DD method in the solution of diffusionconvection problems with dominant convection on three dimensional domains. We studied two problems with internal and external layers, comparing the results obtained with those from the CAU method. In both experiments, the DD method obtained solutions less diffusive than the CAU, but with small under/overshoot in the external layers.

As in two dimensional problems, the DD method proved to be a promising methodology for solving convection-dominated problems on three dimensional domains. However, more numerical experiments need to be carried out in order to strengthen our conclusions. As future work, we plan to apply this methodology for solving 3D unsteady transport problems, including the Navier-Stokes equations, and also to investigate the convergence rates for the errors in norms L 2 (Ω) and H 1 (Ω).

On the numerical solution of the one-dimensional convectiondiffusion equation

Finite Element Method for Flow Problems

Stabilized finite element methods with shock-capturing for nonlinear convection-diffusion-reaction models

A local projection stabilization finite element method with nonlinear crosswind diffusion for convection-diffusionreaction equations

Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems

Streamline upwind Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations

Finite element methods for linear hyperbolic problems

A new finite element formulation for computational fluid dynamics: VIII. the Galerkin/Least-Squares method for advective-diffusive equations

Edge stabilization for Galerkin approximations of convection-diffusion-reaction problems

Blending low-order stabilised finite element methods: a positivity-preserving local projection method for the convection-diffusion equation

A parameter-free dynamic diffusion method for advection-diffusionreaction problems

On spurious oscillations at layers diminishing (SOLD) methods for convection-diffusion equations: part I -a review

On spurious oscillations at layers diminishing (SOLD) methods for convection-diffusion equations: part II -analysis for P1 and Q1 finite elements

Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods

A space-time formulation for multiscale phenomena

The variational multiscale method -a paradigm for computational mechanics

Multiscale and stabilized methods

Choosing bubbles for advection-diffusion problems

On the stability of residual-free bubbles for convection-diffusion problems and their approximation by a two-level finite element method

Stabilization of Galerkin approximations of transport equation by subgrid modeling

Subgrid stabilization of Galerkin approximations of linear monotone operators

A nonlinear subgrid method for advection-diffusion problems

Numerical analysis of the nonlinear subgrid scale method

Dynamic diffusion formulations for advection dominated transport problems

A nonlinear multiscale viscosity method to solve compressible flow problems

Comparative studies to solve compressible problems by multiscale finite element methods

A multiscale finite element formulation for the incompressible Navier-Stokes equations

A consistent approximate upwind Petrov-Galerkin method for convection-dominated problems

Introductionà l'analyse numérique deséquations aux dérivées partielles

The Finite Element Method: Linear Static and Dymanic Finite Element Analysis

Finite element analysis of convection dominated reaction-diffusion problems

A new stabilized finite element formulation for scalar convection-diffusion problems: the streamline and approximate upwind/Petrov-Galerkin method

Theory and Practice of Finite Elements

GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems