Variational Multiscale Modeling with Discontinuous Subscales: Analysis and Application to Scalar Transport

TL;DR

This paper introduces a variational multiscale method with discontinuous subscales for scalar transport, demonstrating stability, optimal convergence, and robustness across different discretizations and flow regimes.

Contribution

It develops a novel multiscale approach using discontinuous Galerkin subscales, providing theoretical analysis and practical validation for scalar transport problems.

Findings

Method is stable and accurate in advective regimes.

Achieves optimal convergence rates in the SUPG norm.

Robust with respect to Peclet number when using rich subscale spaces.

Abstract

We examine a variational multiscale method in which the unresolved fine-scales are approximated element-wise using a discontinuous Galerkin method. We establish stability and convergence results for the methodology as applied to the scalar transport problem, and we prove that the method exhibits optimal convergence rates in the SUPG norm and is robust with respect to the Peclet number if the discontinuous subscale approximation space is sufficiently rich. We apply the method to isogeometric NURBS discretizations of steady and unsteady transport problems, and the corresponding numerical results demonstrate that the method is stable and accurate in the advective limit even when low-order discontinuous subscale approximations are employed.

Figures

Equations

$$\frac{\partial u}{\partial t}+\mathbf{a}\cdot \nabla u-\nabla \cdot \left(\kappa \nabla u\right)$$

$$u$$

$$\kappa \nabla u\cdot \mathbf{n}$$

$$u$$

$$\mathcal{V}:=\left\{v\in H^1\left(\Omega \right):v{|}_{{\Gamma }_D}=0\right\}$$

$$\mathcal{S}:=\left\{u\in H^1\left(\Omega \right):u{|}_{{\Gamma }_D}=g\right\}$$

$${\mathcal{S}}_T:=C([0,T];\mathcal{S})$$

$$\underset{0\le t\le T}{\max }\Vert u(t){\Vert }_{H^1(\Omega )}<\infty .$$

$${\int }_{\Omega }\frac{\partial u}{\partial t}v-{\int }_{\Omega }u\mathbf{a}\cdot \nabla v+{\int }_{\Omega }\kappa \nabla u\cdot \nabla v+{\int }_{{\Gamma }_N}{a}_{n}uv={\int }_{\Omega }fv+{\int }_{{\Gamma }_N}hv$$

$${\mathcal{L}}_t$$

$$\mathcal{V}=\bar{\mathcal{V}}\oplus {\mathcal{V}}^{\prime }$$

$$\mathcal{S}=\bar{\mathcal{S}}+{\mathcal{S}}^{\prime }$$

$${\int }_{\Omega }\frac{\partial \left(\bar{u}+u^{\prime }\right)}{\partial t}\bar{v}-{\int }_{\Omega }\left(\bar{u}+u^{\prime }\right)\mathbf{a}\cdot \nabla \bar{v}+{\int }_{\Omega }\kappa \nabla \left(\bar{u}+u^{\prime }\right)\cdot \nabla \bar{v}+{\int }_{{\Gamma }_N}{a}_n\left(\bar{u}+u^{\prime }\right)\bar{v}={\int }_{\Omega }f\bar{v}+{\int }_{{\Gamma }_N}h\bar{v}$$

$${\int }_{\Omega }\frac{\partial \left(\bar{u}+u^{\prime }\right)}{\partial t}{v}^{\prime }-{\int }_{\Omega }\left(\bar{u}+u^{\prime }\right)\mathbf{a}\cdot \nabla v^{\prime }+{\int }_{\Omega }\kappa \nabla \left(\bar{u}+u^{\prime }\right)\cdot \nabla v^{\prime }+{\int }_{{\Gamma }_N}{a}_n\left(\bar{u}+u^{\prime }\right)v^{\prime }={\int }_{\Omega }f{v}^{\prime }+{\int }_{{\Gamma }_N}h{v}^{\prime }$$

$$\text{Res}\left(\bar{u}\right):=f-{\mathcal{L}}_t\bar{u}$$

$${\mathcal{V}}_{bubble}^{\prime }={\mathcal{S}}_{bubble}^{\prime }:=\left\{v\in \mathcal{V}:v\in H_0^1({\Omega }_K)\text{ for each element }{\Omega }_K\subset \Omega \right\}.$$

$${\int }_{{\Omega }_K}\frac{\partial \left(\bar{u}+u^{\prime }\right)}{\partial t}{v}^{\prime }-{\int }_{{\Omega }_K}\left(\bar{u}+u^{\prime }\right)\mathbf{a}\cdot \nabla v^{\prime }+{\int }_{{\Omega }_K}\kappa \nabla \left(\bar{u}+u^{\prime }\right)\cdot \nabla v^{\prime }={\int }_{{\Omega }_K}f{v}^{\prime }$$

$${\mathcal{L}}_t{u}_K^{\prime }$$

$$u_K^{\prime }$$

$${\mathcal{V}}_{disc}^{\prime }={\mathcal{S}}_{disc}^{\prime }:=\left\{v\in L^2(\Omega ):v{|}_{{\Omega }_K}\in {\mathcal{V}}_{disc,K}^{\prime }\text{ for every }{\Omega }_K\in \mathcal{M}\right\}.$$

$${\int }_{{\Omega }_K}\frac{\partial u^{\prime }}{\partial t}{v}^{\prime }-{\int }_{{\Omega }_K}{u}^{\prime }\mathbf{a}\cdot \nabla v^{\prime }+{\int }_{{\Omega }_K}\kappa \nabla u^{\prime }\cdot \nabla v^{\prime }$$

$$+{\int }_{{\Gamma }_K}{u}^{\prime }{v}^{\prime }{a}_n^{+}-{\int }_{{\Gamma }_K}\kappa \nabla u^{\prime }\cdot \mathbf{n}{v}^{\prime }-{\int }_{{\Gamma }_K}{u}^{\prime }\kappa \nabla v^{\prime }\cdot \mathbf{n}+{\int }_{{\Gamma }_K}\frac{C_{pen}\kappa }{h_K}{u}^{\prime }{v}^{\prime }$$

$${\int }_{\Omega }\frac{\partial \bar{u}}{\partial t}\bar{v}-{\int }_{\Omega }\left(\bar{u}+u^{\prime }\right)\mathbf{a}\cdot \nabla \bar{v}+{\int }_{\Omega }\kappa \nabla \bar{u}\cdot \nabla \bar{v}+{\int }_{{\Gamma }_N}{a}_n\bar{u}\bar{v}={\int }_{\Omega }f\bar{v}+{\int }_{{\Gamma }_N}h\bar{v}$$

$$\text{ for all }\bar{v}\in \bar{\mathcal{V}}$$

$${\int }_{{\Omega }_K}\frac{\partial u^{\prime }}{\partial t}{v}^{\prime }-{\int }_{{\Omega }_K}{u}^{\prime }\mathbf{a}\cdot \nabla v^{\prime }+{\int }_{{\Omega }_K}\left(\kappa +{\kappa }_{art}\right)\nabla u^{\prime }\cdot \nabla v^{\prime }$$

$$+{\int }_{{\Gamma }_K}{u}^{\prime }{v}^{\prime }{a}_n^{+}-{\int }_{{\Gamma }_K}\kappa \nabla u^{\prime }\cdot \mathbf{n}{v}^{\prime }-{\int }_{{\Gamma }_K}{u}^{\prime }\kappa \nabla v^{\prime }\cdot \mathbf{n}+{\int }_{{\Gamma }_K}\frac{C_{pen}\kappa }{h_K}{u}^{\prime }{v}^{\prime }$$

$$\text{ for all }{v}^{\prime }\in {\mathcal{V}}_{disc,K}^{\prime }\text{ and }{\Omega }_K\in \mathcal{M}$$

$${\int }_{{\Omega }_K}{u}^{\prime }\mathbf{a}\cdot \nabla v^{\prime }+{\int }_{{\Omega }_K}\left(\kappa +{\kappa }_{art}\right)\nabla u^{\prime }\cdot \nabla v^{\prime }\cdot \mathbf{n}$$

$$+{\int }_{{\Gamma }_K}{u}^{\prime }{v}^{\prime }{a}_n^{+}-{\int }_{{\Gamma }_K}\kappa \nabla u^{\prime }\cdot \mathbf{n}{v}^{\prime }-{\int }_{{\Gamma }_K}{u}^{\prime }\kappa \nabla v^{\prime }\cdot \mathbf{n}+{\int }_{{\Gamma }_K}\frac{C_{pen}\kappa }{h_K}{u}^{\prime }{v}^{\prime }$$

$$\text{ for all }{v}^{\prime }\in {\mathcal{V}}_{disc,K}^{\prime }\text{ and }{\Omega }_K\in \mathcal{M}$$

$$\mathbf{B}\left({\text{U}}^h,{\text{V}}^h\right):=\bar{B}(\bar{u},\bar{v})+\bar{C}(u^{\prime },\bar{v})+C^{\prime }(\bar{u},v^{\prime })+B^{\prime }(u^{\prime },v^{\prime })$$

$$\bar{B}(\bar{u},\bar{v})$$

$$\bar{C}(u^{\prime },\bar{v})$$

$$C^{\prime }(\bar{u},v^{\prime })$$

$$B^{\prime }(u^{\prime },v^{\prime })$$

$$\mathbf{L}\left({\text{W}}^h\right):=\bar{L}(\bar{v})+L^{\prime }(v^{\prime })$$

$$\bar{L}(\bar{v})$$

Full text

Variational Multiscale Modeling with Discontinuous Subscales:

Analysis and Application to Scalar Transport

Christopher Coley and John A. Evans*∗*

*Ann and H.J. Smead Department of Aerospace Engineering Sciences, University of Colorado Boulder

∗* Corresponding author. E-mail address: john.a.evans@colorado.edu

Abstract

We examine a variational multiscale method in which the unresolved fine-scales are approximated element-wise using a discontinuous Galerkin method. We establish stability and convergence results for the methodology as applied to the scalar transport problem, and we prove that the method exhibits optimal convergence rates in the SUPG norm and is robust with respect to the Péclet number if the discontinuous subscale approximation space is sufficiently rich. We apply the method to isogeometric NURBS discretizations of steady and unsteady transport problems, and the corresponding numerical results demonstrate that the method is stable and accurate in the advective limit even when low-order discontinuous subscale approximations are employed.

1 Introduction

Multiscale phenomena are ubiquitous in science and engineering applications. For instance, turbulent fluid flows are characterized by a continuum of spatial and temporal scales, and even laminar fluid flows exhibit multiscale behavior in the form of boundary and shear layers. Due to the widespread presence and impact of multiscale phenomena, there is a great demand for numerical methods which account for multiscale effects on the numerical solution.

The variational multiscale method was originally introduced as a theoretical framework for incorporating missing fine-scale effects into numerical problems governing coarse-scale behavior [3, 34, 36, 39, 40]. Construction of the variational multiscale method is simple: decompose the solution to a partial differential equation into a sum of coarse-scale and fine-scale components, determine the fine-scale component analytically in terms of the coarse-scale component, and solve for the coarse-scale component numerically. The above scale decomposition is uniquely specified by identifying a projector from the space of all scales onto the coarse-scale subspace. As a consequence, the coarse-scale component is guaranteed to best-fit the solution in a variational sense.

The primary challenge in the variational multiscale method is determining the fine-scale component in terms of the coarse-scale component. Namely, the problem governing the behavior of the fine-scales is infinite-dimensional and nearly as difficult to solve as the original problem of interest. Consequently, in practice, the fine-scale problem must be approximated. The simplest approximations are based on algebraic models which express the fine-scale component in terms of an algebraic expression of the coarse-scale residual. In fact, classical stabilization approaches such as the streamline upwind Petrov Galerkin (SUPG) method [17], the Galerkin Least Squares (GLS) method [37, 38], and the Douglas-Wang method [26] may be viewed as algebraic models. While these models have proven to be quite successful in capturing the mean effect of the unresolved fine-scales on coarse-scale behavior, they typically are unable to account for higher-order moments of fine-scale components [19]. For example, it has been shown that algebraic models yield inaccurate representations of the subgrid stress tensor for turbulent incompressible flows [49].

A more accurate approximation of the fine-scale problem may be obtained using differential models wherein a model differential equation is solved for the fine-scale component. Perhaps the simplest differential model is the method of residual-free bubbles [10, 11, 16, 15, 45]. In this approach, the fine-scale solution space is approximated by the space of bubbles over each element. This yields an element-wise problem which can then be solved for the fine-scale solution field over the element. It should be noted that these element-wise problems are still infinite-dimensional, so they must in turn be solved using a numerical method. Moreover, the fine-scale solution over each element exhibits multiscale features including the presence of boundary layers over element boundaries, so any proposed method must be able to account for such features. In prior work, researchers have turned to Galerkin’s method to solve these problems with either (i) polynomial bubble functions and subgrid viscosity [13, 32] or (ii) piecewise polynomial functions over Shiskin submeshes [12, 14]. The first approach leads to a simple implementation, though the exact value of the required subgrid viscosity is problem-dependent and polynomial bubble functions are unable to accurately represent layers. The second approach leads to a more complex implementation, but it is better able to accurately capture boundary layer phenomena.

In this paper, we examine an alternative approach, which we refer to henceforth as the method of discontinuous subscales, in which the fine-scales in the method of residual-free bubbles are approximated element-wise using a discontinuous Galerkin method [2]. As the discontinuous subscales are not required to meet the residual-free bubble boundary conditions in a strong sense, they are able to represent layers without resorting to a complicated submesh. Herein, we employ the symmetric interior penalty method to approximate diffusive fluxes and the upwind method to approximate advective fluxes on element boundaries in the element-wise fine-scale problems [1, 47]. This yields a stable methodology for transport and incompressible fluid flow even in the advective limit, though it should be mentioned that alternative discontinuous Galerkin methods may be employed [21, 22, 44]. To examine the effectiveness of the method of discontinuous subscales, we conduct a theoretical stability and convergence analysis of the method as applied to the scalar transport problem. Though the scalar transport problem is linear, it exhibits multiscale behavior in the form of boundary and internal layers and thus serves as a simplified vehicle to study the more complicated Navier-Stokes equations. Through our analysis, we find that the method of discontinuous subscales exhibits optimal convergence rates and is robust with respect to the Péclet number if the discontinuous subscale approximation space is sufficiently rich. We further explicitly identify the required richness in the context of affine finite elements. We finish this paper by applying the method to isogeometric discretizations [35] of steady and unsteady transport problems. Our motivation for examining isogeometric discretizations is that they exhibit improved stability and accuracy properties as compared with classical finite elements when applied to difficult transport [5], laminar and turbulent fluid flow [3, 29, 30, 48], and fluid-structure interaction [4] applications. Surprisingly, we find that the method of discontinuous subscales is stable and accurate in the advective limit even when lowest-order discontinuous subscale approximations are employed as a companion fine-scale model for an isogeometric discretization.

It should be mentioned that the method of discontinuous subscales is related to many other multiscale methods in the literature, and in particular, it is closely aligned with the discontinuous residual-free bubble method of Sangalli [46], the multiscale discontinuous Galerkin method of Hughes et al. [8, 18, 41], and the parameter-free variational multiscale method of Cottrell [24]. In the discontinuous residual-free bubble method, one also employs a discontinuous Galerkin method to approximate fine-scale solution behavior, though alternative boundary conditions are employed at the element boundaries. Moreover, no theoretical stability or convergence results exist for the method. In the multiscale discontinuous Galerkin method, variational multiscale analysis and an interscale transfer operator are employed to essentially reduce the computational complexity of a discontinuous Galerkin method to that of a continuous Galerkin method. While the multiscale discontinuous Galerkin method has a firm theoretical foundation, it is fairly difficult to implement, and it is not possible to employ different polynomial degrees for the coarse- and fine-scale components in the method. The parameter-free variational multiscale method is perhaps most closely aligned with the method of discontinuous subscales. In this approach, the discontinuous Galerkin method is employed to approximate the so-called fine-scale Green’s function as opposed to the residual-free bubble directly. No theoretical stability or convergence results exist for the parameter-free variational multiscale method, though we believe that the analysis carried out here can be easily extended to this approach.

An outline of the remainder of the paper is as follows. In Section 2, we provide an overview of the scalar transport problem. In Section 3, we introduce the variational multiscale approach, the method of residual-free bubbles, and the method of discontinuous subscales. In Section 4, we present a theoretical analysis of the method of discontinuous subscales, and in Section 5, we apply the method to a selection of steady and unsteady scalar transport problems. Finally, in Section 6, we provide concluding remarks.

2 The Governing Equations of Scalar Transport

We limit our discussion in this paper to the scalar transport problem, also referred to as the advection-diffusion or drift-diffusion problem in the literature. Though we consider both steady and unsteady transport throughout, we only state the strong and weak forms of the unsteady problem here as the corresponding forms for the steady problem are implied. With the above in mind, the strong form of the unsteady scalar transport problem reads as follows: Find $u:\overline{\Omega}\times\left[0,T\right)\rightarrow\mathbb{R}$ such that:

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$ is the spatial domain of the problem, $d$ is the number of spatial dimensions, $\Gamma$ is the boundary of the domain, $\Gamma_{D}$ is the Dirichlet part of the boundary, $\Gamma_{N}$ is the Neumann part of the boundary, $\bm{n}$ is the unit outward boundary normal, $u$ is the scalar being transported, $\bm{a}:\Omega\rightarrow\mathbb{R}^{d}$ is the advective velocity, $\kappa:\Omega\rightarrow\mathbb{R}_{+}$ is the diffusivity, $f:\Omega\rightarrow\mathbb{R}$ is the applied forcing, $g:\Gamma_{D}\rightarrow\mathbb{R}$ is the prescribed Dirichlet data, $h:\Gamma_{N}\rightarrow\mathbb{R}$ is the prescribed Neumann data, and $u_{0}:\Omega\rightarrow\mathbb{R}$ is the prescribed initial data. In order for the scalar transport problem to be well-posed, we require that $\Gamma=\overline{\Gamma_{D}\cap\Gamma_{N}}$ and that $a_{n}=\bm{a}\cdot\bm{n}\geq 0$ on $\Gamma_{N}$. Note that to simplify our presentation, we have assumed that the advective velocity, diffusivity, forcing, and prescribed data are all independent of time. We further assume throughout that these quantities are all smooth and that the advective velocity is divergence-free.

To establish a weak form of the unsteady scalar transport problem, we must first define a test space and set of trial functions. With this in mind, let:

[TABLE]

and:

[TABLE]

denote the time-instaneous test space and set of trial functions, respectively, where $H^{1}\left(\Omega\right)$ is the Sobolev space of square-integrable functions with square-integrable derivatives. We further define the space-time set of trial functions as:

[TABLE]

which is the set of continuous functions $u:[0,t]\rightarrow\mathcal{S}$ with:

[TABLE]

The weak form of the unsteady advection-diffusion problem then reads as follows: Find $u\in\mathcal{S}_{T}$ such that $u(\cdot,0)=u_{0}(\cdot)$ and:

[TABLE]

for all $v\in\mathcal{V}$ and almost every $t\in(0,T)$. It is well-known that the strong and weak forms of the advection-diffusion problem admit a unique solution that depends smoothly on the prescribed data [31].

To simplify our later discussion, we define the following operator:

[TABLE]

Note that $\mathcal{L}_{t}$ is precisely the differential operator associated with the unsteady scalar transport problem.

3 Variational Multiscale Analysis, Residual-Free Bubbles, and the Method of Discontinuous Subscales

Now that we have defined the governing equations for scalar transport, we are ready to present our numerical method. Before doing so, however, we first review the variational multiscale method followed by the residual-free bubble method.

3.1 The Variational Multiscale Method

In the variational multiscale method, the solution to a partial differential equation is split into a finite-dimensional, coarse-scale component and an infinite-dimensional, fine-scale component through the use of variational projection [40]. As such, the coarse-scale component is a priori guaranteed to best-fit the exact solution in a variational sense.

To proceed forward, we must first define a finite-dimensional, coarse-scale test space $\bar{\mathcal{V}}\subset\mathcal{V}$ and a corresponding continuous, linear projection operator $\mathcal{P}:\mathcal{V}\rightarrow\bar{\mathcal{V}}$. The projection operator naturally splits the test space into coarse-scale and fine-scale components as exhibited by the decomposition:

[TABLE]

where $\mathcal{V}^{\prime}=\ker(\mathcal{P})$ is the infinite-dimensional, fine-scale test space. Consequently, each test function $v\in\mathcal{V}$ is uniquely represented as the sum of a coarse-scale test function $\bar{v}=\mathcal{P}v\in\bar{\mathcal{V}}$ and a fine-scale test function $v^{\prime}=v-\bar{v}\in\mathcal{V}$. Associated with the coarse-scale test space is a corresponding set of trial functions of the form $\bar{\mathcal{S}}=\bar{g}+\bar{\mathcal{V}}$ where $\bar{g}\in\mathcal{S}$ satisfies the non-homogeneous Dirichlet boundary condition $\bar{g}|_{\Gamma_{D}}=g$. This inspires a similar split of the set of trial functions into coarse-scale and fine-scale components:

[TABLE]

where $\mathcal{S}^{\prime}$ is the infinite-dimensional set of fine-scale trial functions. Therefore, the solution to the scalar transport problem $u\in\mathcal{S}$ is uniquely represented as the sum of a coarse-scale component $\bar{u}=\bar{g}+\mathcal{P}u\in\bar{\mathcal{S}}$ and a fine-scale component $u^{\prime}=u-\bar{u}\in\mathcal{S}^{\prime}$. Since the coarse-scale trial functions satisfy the required non-homogeneous boundary condition, the fine-scale trial functions satisfy homogeneous boundary conditions, and thus we have $\mathcal{S}^{\prime}=\mathcal{V}^{\prime}$.

Heretofore, we have discussed how to split the solution to the scalar transport problem into coarse-scale and fine-scale components, but we have not discussed how one may obtain said components via a numerical method. To do so, we simply use the decomposition $\mathcal{V}=\bar{\mathcal{V}}\oplus\mathcal{V}^{\prime}$ and bilinearity to perform a scale splitting of the scalar transport problem. The corresponding variational problem takes the form: Find $\bar{u}\in C\left([0,T];\bar{\mathcal{S}}\right)$ and $u^{\prime}\in C\left([0,T];\mathcal{S}^{\prime}\right)$ such that:

[TABLE]

for all $\bar{v}\in\bar{\mathcal{V}}$ and almost every $t\in(0,T)$ and:

[TABLE]

for all $v^{\prime}\in\mathcal{V}^{\prime}$ and almost every $t\in(0,T)$. Eq. (7) is referred to as the coarse-scale problem and Eq. (8) is referred to as the fine-scale problem. One can solve the fine-scale problem for $u^{\prime}\in\mathcal{S}^{\prime}$ in terms of the so-called coarse-scale residual:

[TABLE]

and insert the resulting solution back into the coarse-scale problem in order to arrive at a final finite-dimensional system for the coarse-scale solution $\bar{u}\in\bar{\mathcal{S}}$.

The primary issue associated with the variational multiscale method is that $\mathcal{V}^{\prime}$ is an infinite-dimensional space and thus solving the fine-scale problem is an intractable task. Fortunately, for most problems of interest, it is sufficient to approximate the effect of the fine-scales on the coarse-scale solution in order to produce stable and accurate numerical solutions [24]. With this in mind, we now turn our attention to differential approaches for modeling the fine-scale problem with specific reference to the method of residual-free bubbles.

3.2 Differential Fine-Scale Modeling with Residual-Free Bubbles

Assume that we have a decomposition of the spatial domain into a mesh of elements $\mathcal{M}=\left\{\Omega_{K}\right\}_{K=1}^{n_{el}}$ satisfying $\overline{\cup_{K}\Omega_{K}}=\overline{\Omega}$. We associate the decomposition with the finite element mesh in the finite element setting, and we associate the decomposition with the Bézier mesh in the context of isogeometric analysis [9]. In the residual-free bubble approach [10, 11, 16], we simply replace the fine-scale test space and set of trial functions in the variational multiscale method by the space of $H^{1}$-conforming bubbles over each element:

[TABLE]

Since the functions in $\mathcal{V}^{\prime}_{bubble}$ and $\mathcal{S}^{\prime}_{bubble}$ are zero-valued on the boundary of each element, the fine-scale problem applies element-by-element. Thus, we have:

[TABLE]

for every $v^{\prime}\in H^{1}_{0}(\Omega_{K})$ and each element $\Omega_{K}\subset\Omega$. By integration-by-parts, the above implies the following element-wise problem:

[TABLE]

where $u^{\prime}_{K}=u^{\prime}|_{\Omega_{K}}$ and $\Gamma_{K}$ is the boundary of the element $\Omega_{K}$. Residual-free bubbles are defined as the functions $u^{\prime}_{K}$ which satisfy the above problem strongly.

Note that the governing equations for the residual-free bubbles are still infinite-dimensional. While these equations are easier to solve than the fine-scale problem associated with the variational multiscale method, they are still nearly as difficult to solve as the original scalar transport problem. Thus, in practice, one must turn to a numerical method to approximate the residual-free bubbles. Additionally, the residual-free bubbles exhibit multiscale behavior such as boundary layers, so any method must be able to account for such behavior. As mentioned in the introduction, a number of approaches have been proposed in prior work based on the continuous Galerkin method. While these approaches have been shown to be closely related to the SUPG method, they exhibit complex features such as problem-dependent subgrid diffusivity [13] or advection-aligned Shishkin subgrid meshes [12] in order to overcome the advective instabilities associated with the continuous Galerkin method and account for the multiscale features in the residual-free bubbles. We instead turn to another means of approximating the residual-free bubbles in the next subsection.

3.3 Approximation of the Residual-Free Bubbles with Discontinuous Subscales

Recall that the residual-free bubbles exhibit the same sorts of multiscale phenomena as the solution to the scalar transport problem, including boundary layers along element boundaries. This inspires us to weaken the enforcement of element-wise homogeneous Dirichlet boundary conditions via a discontinuous Galerkin methodology. We refer to this approach as the method of discontinuous subscales as the corresponding residual-free bubble approximations are discontinuous between adjacent elements.

To set the stage, we need to first define a finite-dimensional test space and set of trial functions for the discontinuous subscales. As with the method of residual-free bubbles, these two sets collapse to the same space, and we simply require that the test and trial functions are $H^{1}$-conforming over each element. We denote the discontinuous subscale approximation space over element $\Omega_{K}$ as $\mathcal{V}^{\prime}_{disc,K}=\mathcal{S}^{^{\prime}}_{disc,K}\subset H^{1}\left(\Omega_{K}\right)$ and the approximation space over the entire domain as:

[TABLE]

There are many potential candidates for the discontinuous subscale approximation space, and we consider polynomial basis functions of a particular degree later in this paper.

It remains now to discretize the governing residual-free bubble equation using our discontinuous subscale approximation space. We turn to a discontinuous Galerkin method in which the symmetric interior penalty method is employed to approximate the diffusive fluxes on element boundaries and the upwind method is employed to approximate the advective fluxes on element boundaries [47]. This approach is known to be stable even in the advective limit, and it results in the following element-wise problem:

[TABLE]

for every $v^{\prime}\in\mathcal{V}^{\prime}_{disc,K}$ where $a^{+}_{n}=\left(\bm{a}\cdot\bm{n}\right)^{+}=\max\left\{\bm{a}\cdot\bm{n},0\right\}$. Above, $h_{K}$ indicates a measure of the size of element $\Omega_{K}$ (e.g., $h_{K}$ could be the diameter of the smallest ball encompassing the element), and $C_{pen}>0$ is a stabilization parameter that must be specifically chosen such that the resulting methodology is stable. It should be noted that the exact value of the stabilization parameter should depend on both the coarse-scale and discontinuous subscale approximation spaces, and we discuss how to choose the stabilization parameter in the next section.

At this juncture, we have arrived at a fully specified method. However, we should highlight a few modifications that yield a more robust and streamlined methodology. First of all, we have found that the optional inclusion of an artificial diffusivity $\kappa_{art}$ in the subscale governing equation generally improves stability as well as the conditioning of the resulting linear system. Guidelines for how to choose the artificial diffusivity are provided in the next section. Second of all, we have found that neglecting the influence of the subscale solution field in the unsteady and diffusive terms appearing in the coarse-scale governing equation also improves stability. It should be mentioned that these terms also vanish if the projection operator is specially chosen as observed in previous works [40]. Finally, one may choose to ignore the time-history of the subscale solution field by removing the time-derivative of the subscales appearing in the subscale governing equation. This approach, which we refer to as the quasi-static model, leads to a simpler implementation than the corrresponding dynamic model. However, both the dynamic and quasi-static models exhibit nearly the same computational cost. It should be noted that the dynamic and quasi-static labels are inspired by the Dynamic Subscales (DSS) and Quasi-Static Subscales (QSS) methods of Codina et al. [23], though the DSS and QSS methods do not attempt to solve for the residual-free bubbles.

Collecting all of our governing equations together, we obtain the following:

Coarse-Scale Governing Equation:

[TABLE]

Dynamic Discontinuous Subscale Model:

[TABLE]

Quasi-Static Discontinuous Subscale Model:

[TABLE]

The governing semi-discrete equations may be discretized in time using any particular time-integrator of interest, yielding a linear system to be solved at each time-step. Note that in both the dynamic and quasi-static discontinuous subscale models, the subscale solution fields on each element are decoupled. Consequently, an element-wise static condensation procedure can be employed to remove the subscale degrees-of-freedom from the discrete system at each time-step, yielding a reduced linear system for the coarse-scale solution field [51]. Moreover, as the subscale solution fields on each element are decoupled, this reduced linear system has exactly the same sparsity pattern as that associated with either Galerkin’s method or the SUPG method.

4 Analysis of the Method of Discontinuous Subscales

Now that we have presented our methodology for solving the scalar transport problem, we establish stability and convergence results. To simplify exposition, we deal solely with the steady advection-diffusion problem with homogeneous Dirichlet boundary conditions. With this in mind, let us define the group vector space $\bm{\mathcal{V}}^{h}=\overline{\mathcal{V}}\times\mathcal{V}^{\prime}_{disc}$, the group variables $\textbf{U}^{h}=\left(\bar{u},u^{\prime}\right)\in\bm{\mathcal{V}}^{h}$ and $\textbf{V}^{h}=\left(\bar{v},v^{\prime}\right)\in\bm{\mathcal{V}}^{h}$, the group bilinear form:

[TABLE]

where:

[TABLE]

and the group linear form:

[TABLE]

where:

[TABLE]

With the above notation established, our discrete problem is as follows: Find $\textbf{U}^{h}\in\bm{\mathcal{V}}^{h}$ such that for all $\textbf{V}^{h}\in\bm{\mathcal{V}}^{h}$:

[TABLE]

Throughout, we assume that $\overline{\mathcal{V}}$ is a finite-dimensional space of continuous piecewise polynomial or tensor-product polynomial functions of degree $p$ defined over a given mesh $\mathcal{M}$ of simplices (triangles and tetrahedra) or parallelotopes (quadrilaterals and hexahedra). We also assume that $\mathcal{V}^{\prime}_{disc}$ is a finite dimensional-space of discontinuous piecewise polynomial or tensor-product polynomials of degree $p_{f}$ defined over the same mesh $\mathcal{M}$. While our analysis only strictly covers the setting of affinely-mapped finite elements, it easily extends to the more general settings of curvilinear finite elements and isogeometric analysis.

Throughout this section, we make use of the classical Lebesgue spaces $L^{q}(D)$ endowed with the norm $\|\cdot\|_{L^{q}(D)}$ where $1\leq q\leq\infty$ and $D\subset\mathbb{R}^{d}$ is a generic open domain for integer $d\geq 1$. We will also utilize the Sobolev spaces $W^{k,q}(D)$ for $k$ a non-negative integer and $1\leq q\leq\infty$, endowed with the norm:

[TABLE]

and semi-norm:

[TABLE]

In the setting when $q=2$, the Sobolev spaces $W^{k,q}(D)$ become $H^{k}(D)$ and we use the simplified notation $\|\cdot\|_{k,D}\equiv\|\cdot\|_{k,2,D}$ and $|\cdot|_{k,D}\equiv|\cdot|_{k,2,D}$. Sobolev norms are defined over boundaries of open domains in an analogous manner.

To proceed forward, we must make certain assumptions regarding the form of the penalty constant $C_{pen}$ and the artificial diffusivity $\kappa_{art}$. In particular, we assume the following:

Assumption 1: The penalty constant $C_{pen}$ satisfies $C_{pen}\geq 8C_{trace}$ where $C_{trace}>0$ is a sufficiently large positive constant such that:

[TABLE]

*for every $\bar{v}\in\overline{\mathcal{V}}$, $v^{\prime}\in\mathcal{V}^{\prime}_{disc}$, and element $\Omega_{K}\in\mathcal{M}$.

Assumption 2: The artificial diffusivity $\kappa_{art}$ takes the form:

[TABLE]

where $C_{art}\geq 0$ is an arbitrary non-negative constant and $C_{inv}>0$ is a sufficiently large positive constant such that:

[TABLE]

*for every $\bar{v}\in\overline{\mathcal{V}}$, $v^{\prime}\in\mathcal{V}^{\prime}_{disc}$, and element $\Omega_{K}\in\mathcal{M}$.

One can readily find the trace and inverse constants associated with Assumptions 1 and 2 by solving element-wise eigenproblems. Alternatively, explicit bounds for these constants are available both in the setting of finite elements and isogeometric analysis [5, 28, 33, 50].

The following result shows that the group bilinear form is coercive, and hence the group variable solution is unique.

Theorem 4.1

The coercivity result

[TABLE]

holds for all $\textup{{V}}^{h}\in\bm{\mathcal{V}}^{h}$ provided Assumption 1 is satisfied.

• Proof

For the coarse-scale bilinear form, we evaluate:

[TABLE]

Similarly, for the fine-scale bilinear form, we evaluate:

[TABLE]

For the bilinear forms coupling the coarse-scales and fine-scales, we evaluate:

[TABLE]

To continue, we recognize that, for every element $\Omega_{K}\in\mathcal{M}$:

[TABLE]

by Young’s inequality, and:

[TABLE]

by the trace inequality. By a similar argument, we have:

[TABLE]

for every element $\Omega_{K}\in\mathcal{M}$. Finally, by the triangle inequality, we have:

[TABLE]

for every element $\Omega_{K}\in\mathcal{M}$. Collecting our results contained in (23)-(28), we obtain:

[TABLE]

The desired expression follows by adding the above inequalities and invoking Assumption 1.

The coercivity result provided in Theorem 4.1 ensures that the group variable solution is unique, but it does not ensure that the corresponding methodology is stable. In fact, note that in the limit $\kappa\rightarrow 0$, the coercivity result suggests that the methodology loses control of the coarse-scale solution entirely. However, we are able to show that the methodology satisfies a different notion of stability, namely inf-sup stability, in this limit. To proceed forward, let us define the norm:

[TABLE]

for every $\textbf{V}^{h}\in\bm{\mathcal{V}}^{h}$ where $\mathcal{P}^{\prime}_{K}$ is the $L^{2}$-projector onto the space of discontinuous subscales $\mathcal{V}^{\prime}_{disc,K}$ on element $\Omega_{K}\in\mathcal{M}$ and $\tau_{K}$ is the element-wise constant:

[TABLE]

where $C_{\tau}>0$ is an arbitrary positive constant. With the above norm defined, we have the following inf-sup stability result.

Theorem 4.2

The inf-sup stability result

[TABLE]

holds provided Assumptions 1 and 2 are satisfied where $\beta$ is a positive constant independent of the problem parameters $\bm{a}$ and $\kappa$ and the mesh size $h=\max_{K}h_{K}$.

• Proof

Let $\textbf{U}^{h}=(\bar{u},u^{\prime})$ be an arbitrary member of $\bm{\mathcal{V}}^{h}$. It suffices to show that there exists some $\textbf{V}^{h}=(\bar{v},v^{\prime})\in\bm{\mathcal{V}}^{h}$ such that $\textbf{B}(\textbf{U}^{h},\textbf{V}^{h})\geq\beta\|\textbf{U}^{h}\|_{S}\|\textbf{V}^{h}\|_{S}$. In this direction, let $z^{\prime}\in\mathcal{V}^{\prime}_{disc}$ be defined such that $z^{\prime}|_{\Omega_{K}}=\tau_{K}\mathcal{P}_{K}^{\prime}\left(\bm{a}\cdot\nabla\left(\bar{u}+u^{\prime}\right)\right)$. Note that we have the following inequalities bounding the element-wise $L^{2}$-norm of $z^{\prime}$:

[TABLE]

the element-wise $L^{2}$-norm of the gradient of $z^{\prime}$:

[TABLE]

again the element-wise $L^{2}$-norm of the gradient of $z^{\prime}$:

[TABLE]

the element-wise $L^{2}$-norm of the trace of $z^{\prime}$:

[TABLE]

again the element-wise $L^{2}$-norm of the trace of $z^{\prime}$:

[TABLE]

and the element-wise $L^{2}$-norm of the normal derivative trace of $z^{\prime}$:

[TABLE]

We need one more inequality for $z^{\prime}$. Notably, observe:

[TABLE]

Let us define $\textbf{Z}^{h}=(0,z^{\prime})$. Note immediately that as a consequence of the above inequalities, the following inequality holds:

[TABLE]

where:

[TABLE]

We now seek a lower bound for the quantity $\textbf{B}(\textbf{U}^{h},\textbf{Z}^{h})$. We observe that:

[TABLE]

We deal with each of the expressions appearing above one-by-one. First, we note that, by definition:

[TABLE]

For the next term, we have:

[TABLE]

where $\gamma_{1}>0$ is an arbitrarily chosen positive number. In analogous fashion, we have:

[TABLE]

[TABLE]

[TABLE]

where again $\gamma_{2},\gamma_{3},\gamma_{4}>0$ are arbitrarily chosen positive numbers. The remaining terms require slightly more care. For the first remaining term, we have:

[TABLE]

where $\gamma_{5}>0$ is another arbitrarily chosen positive number. For the second remaining term, we have:

[TABLE]

where $\gamma_{6}>0$ is yet another arbitrarily chosen positive number. For the last remaining term, we integrate by parts, resulting in:

[TABLE]

The two terms on the right-hand-side are then easily bound as before, yielding the inequality:

[TABLE]

where $\gamma_{7}>0,\gamma_{8}>0$ are two final arbitrarily chosen positive numbers. Collecting all of the above inequalities, we obtain the composite inequality:

[TABLE]

wherein:

[TABLE]

[TABLE]

We now assume that $\gamma_{1}$ through $\gamma_{8}$ are chosen sufficiently small to guarantee that $C_{1}>0$. Note that we can choose such constants independent of the problem parameters and the mesh size. This choice in turn defines the constants $C_{2}$ through $C_{6}$. We are now in a position to define a suitable group test function $\textbf{V}^{h}=(\bar{v},v^{\prime})\in\bm{\mathcal{V}}^{h}$. Namely, we select $\bar{v}=\bar{u}$ and $v^{\prime}=u^{\prime}+C_{lin}z^{\prime}$ where:

[TABLE]

Then, by Theorem 4.1 and (52), it follows that:

[TABLE]

with:

[TABLE]

and:

[TABLE]

with $C_{v}=1+C_{lin}C_{z}$. Thus the desired condition holds with $\beta=C_{bound}/C_{v}$ independent of the problem parameters $\bm{a}$ and $\kappa$ and the mesh size $h$.

We now introduce one more assumption. This assumption guarantees that our methodology is at least as stable as the SUPG method for steady scalar transport, as is shown in Corollary 4.3.

Assumption 3: The following inequality holds for each element $\Omega_{K}\in\mathcal{M}$:

[TABLE]

where $C_{SUPG}>0$ is a positive constant independent of the problem parameters $\bm{a}$ and $\kappa$ and the mesh size $h$.

Corollary 4.3

Provided Assumptions 1, 2, and 3 hold, the following inf-sup stability result is satisfied:

[TABLE]

where:

[TABLE]

for every $\textbf{{V}}^{h}\in\bm{\mathcal{V}}^{h}$ and $\beta_{SUPG}$ is a positive constant independent of the problem parameters $\bm{a}$ and $\kappa$ and the mesh size $h$.

There remains the question of whether or not we can expect Assumption 3 to be satisfied for a given finite element discretization. The following lemma demonstrates that Assumption 3 is indeed satisfied if the discontinuous subscale solution space is sufficiently rich and the artificial diffusivity is chosen in an intelligent manner.

Lemma 4.4

Assumption 3 is satisfied provided that one of the following two conditions is satisfied:

C1:* It holds that $\bm{a}\cdot\nabla\left(\bar{v}+v^{\prime}\right)\in\mathcal{V}^{\prime}_{disc}$ for every $(\bar{v},v^{\prime})\in\bm{\mathcal{V}}^{h}$.*

C2:* It holds that $\bm{a}\cdot\nabla\bar{v}\in\mathcal{V}^{\prime}_{disc}$ for every $\bar{v}\in\overline{\mathcal{V}}$ and Assumption 2 is satisfied with $C_{art}>0$.*

• Proof

We write:

[TABLE]

where $\mathcal{I}$ is the identity operator. If Condition C1 holds, it follows that:

[TABLE]

and consequently the lemma is satisfied with $C_{SUPG}=1$. If Condition C2 holds, we instead have:

[TABLE]

As Assumption 2 is satisfied with $C_{art}>0$, it follows that:

[TABLE]

and consequently the lemma is satisfied with $C_{SUPG}=\max\left\{1,\frac{C_{\tau}}{C_{art}}\right\}$.

Suppose that the imposed velocity field is a polynomial function of degree $q$ over each element $\Omega_{K}\in\mathcal{M}$. Then, for the case of a finite element mesh of simplices, we see that $\bm{a}\cdot\nabla\bar{v}$ is a polynomial function of degree $p+q-1$ over each element and $\bm{a}\cdot\nabla v^{\prime}$ is a polynomial function of degree $p_{f}+q-1$ over each element. Consequently, if $p_{f}\geq p+q-1$ and $q\leq 1$, then Condition C1 of Lemma 4.4 is satisfied, and if only $p_{f}\geq p+q-1$ is satisfied but $C_{art}$ is chosen to be a positive number, then Condition C2 of Lemma 4.4 is satisfied. For a continuous piecewise linear finite element discretization with discontinuous piecewise linear subscales, we see that Condition C1 of Lemma 4.4 is satisfied for a piecewise linear velocity field. Alternately, for a continuous piecewise quadratic finite element discretization with discontinuous piecewise linear subscales, we see that Condition C1 of Lemma 4.4 is satisfied for a piecewise constant velocity field. By enriching the subscale space to discontinuous piecewise quadratic subscales, we find Condition C1 of Lemma 4.4 is satisfied again for a piecewise linear velocity field.

Lemma 4.4 provides a general guideline for how to choose the polynomial degree of the subscale space. Nonetheless, we have observed our methodology often returns accurate and stable results even when the conditions of the lemma are not satisfied. In particular, we have observed our methodology is stable if we employ smooth splines for our coarse-scale solution and discontinuous piecewise bi-linear finite elements for our subscale solution. We anticipate that Lemma 4.4 holds in this case, though such an analysis is beyond the scope of the current work.

Our final theorem demonstrates that our method exhibits optimal convergence rates with respect to the SUPG norm provided Assumptions 1, 2, and 3 hold.

Theorem 4.5

Suppose that the exact solution satisfies the smoothness condition $u\in H^{p+1}(\Omega)$ and Assumptions 1, 2, and 3 hold. Then, the error $\textup{{E}}=(u-\bar{u},-u^{\prime})$ satisfies the a priori estimate:

[TABLE]

where $C_{apriori}$ is a positive constant independent of the problem parameters $\bm{a}$ and $\kappa$ and the mesh size $h$.

• Proof

Let $\tilde{u}\in\overline{\mathcal{V}}$ be an interpolation function which we will define later, and let us split the error into two components, a method error defined as $\textbf{E}^{h}=(\tilde{u}-\bar{u},-u^{\prime})=(e^{h},e^{\prime})$ and an interpolation error defined as $\bm{\eta}=(u-\tilde{u},0)=(\eta,0)$. Our first objective is to bound the method error by the interpolation error. By Corollary 4.3, we have that:

[TABLE]

Furthermore, it is readily shown that our method is consistent, that is, $\textup{{B}}(\textup{{E}},\textup{{V}}^{h})=0$, so we further have that:

[TABLE]

We now require a bound on the term $\textup{{B}}(\bm{\eta},\textup{{V}}^{h})$. Direct substitution results in:

[TABLE]

We can bound the first term on the right hand side of (59) as:

[TABLE]

The second term on the right hand side of (59) requires more care. We first integrate by parts to obtain:

[TABLE]

To proceed, we require bounds for each of the four terms on the right hand side of (61). The first and third terms are easily bounded as before, yielding:

[TABLE]

Similarly, we can bound the second and fourth terms like:

[TABLE]

Collecting (59)-(65) and applying the Cauchy-Schwarz inequality, we obtain the composite inequality:

[TABLE]

Combining the above expression with (58), we finally obtain the method error bound:

[TABLE]

Now suppose we have chosen the interpolant $\tilde{u}$ to be a “best” interpolant such that the following local interpolation estimates hold for every element $\Omega_{K}\in\mathcal{M}$:

[TABLE]

where $C_{shape}>0$ is a positive constant independent of the mesh size but possibly dependent on the mesh regularity and coarse-scale polynomial degree $p$ [20]. It then follows that:

[TABLE]

Combining (67)-(70) with the method error bound (66) yields:

[TABLE]

where:

[TABLE]

The interpolation error may be bounded in a similar manner, resulting in:

[TABLE]

where $C_{interpolation}>0$ is a positive constant independent of the mesh size but possibly dependent on the mesh regularity and coarse-scale polynomial degree $p$. The desired result follows by combining (71) and (72) with the decomposition $\textbf{E}=\textbf{E}^{h}+\bm{\eta}$.

Note remarkably that the error estimate characterized by Theorem 5.5 is completely independent of the subscale polynomial degree. Consequently, the subscales act purely to stabilize the coarse-scale solution and have virtually no impact on solution accuracy. This is one of the main messages of the current work: it is possible to stabilize a high-order numerical method with a low-order, but stable, numerical method.

It should also be remarked that while the constants appearing in all of the above estimates are independent of the problem parameters $\bm{a}$ and $\kappa$ and the mesh size $h$, they are possibly dependent on mesh regularity, the coarse-scale polynomial degree $p$, and the subscale polynomial degree $p_{f}$. A more delicate analysis is required to obtain dependencies with respect to these variables. This would be useful, for instance, in the context of a boundary layer application wherein a highly skewed mesh should be employed near solid boundaries.

5 Numerical Results

We finish this paper by applying our methodology to a collection of two-dimensional steady and unsteady transport problems. All of the following results were obtained by applying the method of discontinuous subscales to isogeometric Non-Uniform Rational B-splines (NURBS) discretizations on uniform grids [35]. These discretizations employ $C^{p-1}$-continuous piecewise tensor-product polynomials or rational functions of degree $p$. For $p=1$, isogeometric NURBS discretizations correspond to standard bilinear finite elements, while for $p>1$, isogeometric NURBS discretizations exhibit enhanced continuity as compared to standard tensor-product finite elements. For more information on NURBS-based isogeometric analysis, refer to [25]. Throughout, we used the following values for $C_{pen}$ and $\kappa_{art}$:

[TABLE]

and:

[TABLE]

where $p_{f}$ is the subscale polynomial degree and $h_{K}$ is the mesh size for a given element $\Omega_{K}\in\mathcal{M}$. These are inspired by theoretical estimates for the trace and inverse constants associated with a discontinuous subscale discretization [33, 50]. Note that these values are simpler than the ones proposed in the previous section, but we have found that they still yield stable and accurate results.

5.1 Steady manufactured solution

We begin by considering a simple manufactured solution $u\left(x,y\right)=\sin\left(\pi x\right)\sin\left(\pi y\right)$ to the steady scalar transport problem. This solution is realized on the domain $\Omega=\left[0,1\right]\times\left[0,1\right]$ by setting the velocity to $\bm{a}=\left(\sqrt{2}/2,\sqrt{2}/2\right)$, setting the forcing to:

[TABLE]

and applying homogeneous Dirichlet boundary conditions along the entire domain boundary. To ensure the problem is advection-dominated, the diffusion is set to $\kappa=10^{-6}$. The problem is then characterized by the Péclet number $Pe=\frac{|\bm{a}|L}{\kappa}=10^{6}$ where $L$ is the length of the domain. Numerical solutions are calculated on grids of $2\times 2$, $4\times 4$, $8\times 8$, $16\times 16$, $32\times 32$, $64\times 64$, and $128\times 128$ elements for coarse-scale polynomial degrees $p=1,2$ and subscale polynomial degrees $p_{f}=1,2$. The $L^{2}$-norms of the error in the coarse-scale solution and the subscale solution are presented in Fig. 1. We observe that both the coarse-scale and subscale solutions converge optimally under mesh refinement.

5.2 Steady advection skew to the mesh

We next consider the classical two-dimensional advection skew to mesh problem. This problem is graphically depicted in Fig. 2. The Péclet number for the problem is $Pe=10^{6}$, so the problem is advection-dominated. Throughout this subsection, we select $\theta=45^{\circ}$.

First, we explore the effect of mesh refinement. Solutions on uniform NURBS meshes with $32\times 32$, $64\times 64$, and $128\times 128$ elements for both $p=1$ and $p=2$ are computed and compared. For each of the cases, the subscale polynomial degree is chosen to be $p_{f}=1$. A representative solution for $p=2$, $p_{f}=1$, and a $64\times 64$ mesh is shown in Fig. 3. The effect of mesh refinement is illustrated by comparison of the obtained solutions along the slice $y=0.7$ in Fig. 4. The results are as expected. Refinement of the mesh allows for a more accurate representation of the boundary and internal layers. The results are polluted by the presence of oscillations near the boundary and internal layers, but for $p=1$, these oscillations are limited to a region of one or two elements away from the layers. For $p=2$, the oscillations do infiltrate further into the domain, but it is our experience that these oscillations are not due to method instability but rather the fact that we are trying to fit the sharp gradients present in the boundary and internal layers. As isogeometric NURBS discretizations exhibit a high level of continuity, any oscillations due to overfitting such gradients are expected to extend into the domain in analogy with Gibbs’ phenomena. Nonetheless, the oscillations decay in magnitude away from the layers and are still limited to a region of a fixed number of elements away from the layers. As we will later see, these oscillations can be greatly suppressed through the use of weak enforcement of Dirichlet boundary conditions.

We now explore the effect of degree elevation. First, we examine the effect of coarse-scale degree elevation. Solutions for $p=1,2,3$ with $p_{f}=1$ on a uniform NURBS mesh with $64\times 64$ elements are computed and compared. The effect of coarse-scale degree elevation is illustrated by comparison of the obtained solutions along the slice $y=0.7$ in Fig. 5. Note that each of the coarse-scale solutions are able to capture the internal and boundary layers, and the coarse-scale solutions associated with higher polynomial degrees exhibit sharper representations of the layers. However, the coarse-scale solutions associated with higher polynomial degrees also exhibit oscillations in the regions near the internal and boundary layers due to overfitting and Gibbs’ phenomena. Next, we explore the effect of subscale degree elevation. Solutions for $p=2$ with $p_{f}=1,2,3$ are computed on a uniform NURBS mesh with $64\times 64$ elements and compared. The effect of fine-scale degree elevation is illustrated by comparison of the obtained solutions along the slice $y=0.7$ in Fig. 5. From the figure, it is seen that the choice of subscale polynomial degree does not greatly affect the coarse-scale solution. Consequently, it is recommended that one choose the lowest possible subscale polynomial degree in order to stabilize the higher-order coarse-scale solution at minimal computational cost.

Heretofore, we have enforced the Dirichlet boundary conditions on the coarse-scale solution field in a strong manner. However, it has been noted in recent work that greatly enhanced results may be obtained by instead imposing these in a weak manner using a combination of upwinding and Nitsche’s method [6, 7]. The resulting coarse-scale governing equation takes the form:

Coarse-Scale Governing Equation with Weak Boundary Condition Treatment:

[TABLE]

In the above equation, $C_{Nitsche}$ is a penalty constant that must be chosen sufficiently large for method stability, $a^{+}_{n}=\left(\bm{a}\cdot\bm{n}\right)^{+}=\max\left\{\bm{a}\cdot\bm{n},0\right\}$, and $a^{-}_{n}=\left(\bm{a}\cdot\bm{n}\right)^{-}=\min\left\{\bm{a}\cdot\bm{n},0\right\}$. Herein, we choose $C_{Nitsche}=4p$. It should be mentioned that the discontinuous subscale model remains untouched if one elects to weakly enforce Dirichlet boundary conditions on the coarse-scale solution.

To assess the effect of weak boundary condition enforcement, we have computed solutions using both using strongly-enforced and weakly-enforced Dirichlet boundary conditions with $p=2$ and $p_{f}=1$ on a uniform NURBS mesh with $64\times 64$ elements. The effect of boundary condition enforcement is illustrated by comparison of the obtained solutions along the slice $y=0.7$ in Fig. 6. Note from the figure that there no longer remain any oscillations near the boundary layer if one employs a weak boundary condition enforcement. Moreover, while there remain oscillations near the internal layer, these oscillations are small in magnitude and limited to a region of a small number of elements away from the layer. Consequently, it is strongly advised that one weakly enforces Dirichlet boundary conditions along portions of the boundary where a layer is expected.

5.3 Steady advection in a quarter-annulus

We next consider a problem posed on a non-square geometry, namely the advection of a Gaussian curve in a quarter annulus. This problem is graphically depicted in Fig. 7. Here, $r_{i}$ is the inner radius of the annulus and $r_{o}$ is the outer radius of the annulus and they are chosen to be $r_{i}=1$ and $r_{o}=2$. The flow field for this problem is chosen to be $\bm{a}=\left(y,-x\right)$ and the diffusivity is chosen as $\kappa=1\times 10^{6}$ such that the problem is advection-dominated. Dirichlet boundary conditions are strongly enforced on the lower boundary of the annulus with:

[TABLE]

where $\sigma=0.05$ and $\mu=-1.5$. Traction-free Neumann boundary conditions are prescribed on the other three boundaries. A typical solution corresponding to $p=2$, $p_{f}=1$, and a NURBS mesh of $64\times 64$ elements is displayed in Fig. 8. From the figure, it is clear that the method is able to very accurately capture the solution free of oscillations on the curved geometry.

5.4 Unsteady advection of a Gaussian hill

We now consider an unsteady scalar transport problem, namely the advection of a Gaussian hill in a rotating flow field. For this problem, the domain is $\Omega=[0,1]\times[0,1]$ and the flow field is given as $\bm{a}=\left(y-0.5,x-0.5\right)$. The initial scalar field is given as:

[TABLE]

where $\sigma=0.065$, and Dirichlet boundary conditions are strongly enforced along the entire boundary of the domain. Numerical solutions are calculated on a $64\times 64$ NURBS mesh with coarse-scale polynomial degree $p=2$ and subscale polynomial degree $p_{f}=1$. Both the dynamic and quasi-static discontinuous subscale models were employed to stabilize the solution. Time-integration is carried out with the generalized-alpha method using a time step of $\Delta t=0.005$ and $\rho_{\infty}=1$ such that $\alpha_{m}=0.5$, $\alpha_{f}=0.5$, and $\gamma=0.5$. The small time step is chosen such that the only source of error is spatial error.

Establishing the initial condition for this problem takes some care, and the approach differs between the dynamic and quasi-static models. For the quasi-static model, we do not require a time-history of the discontinuous subscale solution. Thus, only an initial condition for the coarse-scale solution field is required, and this is established through an $L^{2}$-projection of the exact initial condition. For the dynamic model, a time-history of the discontinuous subscale solution is required. Consequently, after an initial condition for the coarse-scale solution field is obtained, an initial condition for the discontinuous subscale solution is established through an $L^{2}$-projection of the difference between the exact initial condition and the coarse-scale initial condition.

5.4.1 Results for the Quasi-static Discontinuous Subscale Model

We first display results obtained using the quasi-static discontinuous subscale model. In Fig. 9, the concentration along the centerline of the hill (found at $y=0.25$) is plotted for the initial solution and the solution after 1, 5, and 10 full revolutions ($t=6.28,31.40,62.80$, respectively). It is clear that the method is able to preserve the hill with very high accuracy even after several rotations. After one rotation, the peak is reduced to $0.997$, representing a numerical dispersion of the peak of 0.3%. After ten rotations, the peak is reduced to $0.969$, representing a numerical dispersion of the peak of 3.1%. There is also some slight skewing of the solution in the downwind direction.

5.4.2 Results for the Dynamic Discontinuous Subscale Model

We next display results obtained using the dynamic discontinuous subscale model. In Fig. 9, the concentration along the centerline of the hill (found at $y=0.25$) is plotted for the initial solution and the solution after 1, 5, and 10 full revolutions ($t=6.28,31.40,62.80$, respectively). Like the quasi-static model, it is clear that the method is superb in preserving the hill. After one rotation, the peak is reduced to $0.996$, representing a numerical dispersion of the peak of 0.4%. After ten rotations, the peak is reduced to $0.965$, representing a numerical dispersion of the peak of 3.5%. Again, like the quasi-static model, there is some slight skewing of the solution in the downwind direction. Curiously, the magnitude of the subscale solution is three orders less for the dynamic model than the quasi-static model, but the two models produce very similar results both quantitatively and qualitatively.

5.5 Unsteady advancing front

We finish by returning back to the advection skew to the mesh problem but in an unsteady context. Namely, we initialize the scalar field to be zero everywhere except at the Dirichlet boundary where the initial condition meets the specified boundary condition. We then let the solution evolve according to the problem specifications in Fig. 2. Since we expect a boundary layer to form for this problem, we elect to enforce the Dirichlet conditions in a weak fashion. For all our numerical results, a uniform NURBS mesh with $64\times 64$ square elements was employed with with a coarse-scale polynomial degree of $p=2$ and a subscale polynomial degree of $p_{f}=1$. Time-integration is carried out using the generalized-alpha method using a time step of $\Delta t=0.01$ and $\rho_{\infty}=1$ such that $\alpha_{m}=0.5$, $\alpha_{f}=0.5$, and $\gamma=0.5$. We compute and compare solutions using both the dynamic and quasi-static discontinuous subscale models.

5.5.1 Results for the Quasi-static Discontinuous Subscale Model

Solutions for the advancing front problem solved with the quasi-static discontinuous subscale model at representative time steps of $t=0,0.51,2.01$ are shown in Figs. 10 and 11. The coarse-scale solution converges to the steady-state solution as expected. The coarse-scale solution is quite accurate in the eyeball norm, though at each time step, minor oscillations are present near the internal layer similar to what was observed in the steady skew to the mesh problem. Additionally, the magnitude of the oscillations in the fine-scale solution at the final time-step correspond to roughly the same magnitude of the oscillations observed in the steady problem. This is expected, since we expect the steady state subscale solution to coincide with the subscale solution in the steady problem.

5.5.2 Results for the Dynamic Discontinuous Subscale Model

Solutions for the advancing front problem solved with the dynamic discontinuous subscale model at representative time steps of $t=0,0.51,2.01$ are shown in Figs. 10 and 11. Again, the coarse-scale solution is quite accurate and converges to the steady-state solution as expected, and minor oscillations are present near the internal layer similar to what was observed in the steady skew to the mesh problem. Surprisingly, the magnitude of the subscale solution is much smaller than the subscale solution computed using the quasi-static discontinuous subscale model at intermediate times, yet the obtained numerical results for the coarse-scale solution field are qualitatively and quantitatively similar. The coarse-scale solution field for the dynamic model does exhibit fewer oscillations than the coarse-scale field for the quasi-static model, particularly away from the boundary layer (e.g., near the spatial location $(x,y)=(1,1)$ at time $t=2.01$).

6 Conclusions

In this paper, we examined a variational multiscale method, which we refer to as the method of discontinuous subscales, that is based on approximating the residual-free bubbles using a discontinuous Galerkin formulation. We established stability and convergence results for the methodology and demonstrated its applicability to scalar transport problems through a collection of numerical examples. We demonstrated that the method is accurate, displaying optimal convergence rates with respect to mesh refinement, and stable in the advective-limit. We also demonstrated that, somewhat surprisingly, lowest-order discontinuous subscale approximations are sufficient to stabilize high-order coarse-scale approximations in the context of isogeometric analysis.

There are two main avenues of research that we plan to pursue in future work. The first avenue is the design of discontinuous subscale approaches which yield coarse-scale solutions satisfying positivity and monotonicity constraints. This can be achieved through the use of discontinuity capturing operators [42, 43] or by directly embedding the constraints within the variational multiscale framework [27]. The second avenue of research we plan to pursue is the design of discontinuous subscale approaches for LES-type turbulence modeling. While lowest-order discontinuous subscale approximations were found to stabilize coarse-scale approximations in the context of scalar transport, we expect that such approximations will not be sufficient to yield accurate representations of the subgrid stress tensor for turbulent incompressible flows [49]. We anticipate that improved LES turbulence models may be attained through adaptive refinement of the discontinuous subscale approximation space.

Acknowledgements

This material is based upon work supported by the Air Force Office of Scientific Research under Grant No. FA9550-14-1-0113. The authors would also like to acknowledge early conversations with Thomas J.R. Hughes and J. Austin Cottrell which motivated the subject of this paper.