Library Subscription: Guest
International Journal for Uncertainty Quantification

Published 6 issues per year

ISSN Print: 2152-5080

ISSN Online: 2152-5099

The Impact Factor measures the average number of citations received in a particular year by papers published in the journal during the two preceding years. 2017 Journal Citation Reports (Clarivate Analytics, 2018) IF: 1.7 To calculate the five year Impact Factor, citations are counted in 2017 to the previous five years and divided by the source items published in the previous five years. 2017 Journal Citation Reports (Clarivate Analytics, 2018) 5-Year IF: 1.9 The Immediacy Index is the average number of times an article is cited in the year it is published. The journal Immediacy Index indicates how quickly articles in a journal are cited. Immediacy Index: 0.5 The Eigenfactor score, developed by Jevin West and Carl Bergstrom at the University of Washington, is a rating of the total importance of a scientific journal. Journals are rated according to the number of incoming citations, with citations from highly ranked journals weighted to make a larger contribution to the eigenfactor than those from poorly ranked journals. Eigenfactor: 0.0007 The Journal Citation Indicator (JCI) is a single measurement of the field-normalized citation impact of journals in the Web of Science Core Collection across disciplines. The key words here are that the metric is normalized and cross-disciplinary. JCI: 0.5 SJR: 0.584 SNIP: 0.676 CiteScore™:: 3 H-Index: 25

Indexed in

PARAMETER SENSITIVITY OF AN EDDY VISCOSITY MODEL: ANALYSIS, COMPUTATION AND ITS APPLICATION TO QUANTIFYING MODEL RELIABILITY

Faranak Pahlevani

Division of Science and Engineering, Penn State Abington, Abington, Pennsylvania 19001, USA

Lisa Davis

Department of Mathematical Sciences, Montana State University, Bozeman, Montana 59717, USA

Abstract

An eddy viscosity model can be used as a computationally tractable alternative to that of the Navier-Stokes equations. Model errors immediately become a concern when considering such an approach, and quantifying this error is essential to understanding and using model predictions within an engineering design process. In this paper, sensitivity analysis is presented for a subgrid eddy viscosity model with respect to variations of the eddy viscosity parameter. We demonstrate the analysis utilizing the sensitivity equation method. Approximating the sensitivity requires the solution of the eddy viscosity model. Therefore, the eddy viscosity model and sensitivity equation are coupled in our analysis and computations. An implicit-explicit time-stepping method is developed and analyzed for this set of equations. Our numerical assessments present the role of the sensitivity in quantifying the modeling error arising from the choice of various values of the eddy viscosity parameter. The sensitivity computation allows one to identify an interval of reliability for the eddy viscosity parameter. This gives the user a range of parameter values for which the eddy viscosity model can be considered to be a reliable approximation to the Navier-Stokes equations. A two-dimensional cavity problem is used to illustrate the ideas. In addition, for the standard model problem of two-dimensional flow around a cylinder, the sensitivity computations are shown to be very useful in improving the flow functional approximations that may be used within an optimal design algorithm.

KEYWORDS: large eddy simulation, eddy viscosity model, sensitivity analysis, uncertainty quantification


1. INTRODUCTION

Turbulence is the center of attention for many scientists because of its practical applications in different branches of science. In predictions of turbulent flows, it is usual to estimate a suitable average instead of pointwise fluid velocities. Large eddy simulation, or LES, is a technique for simulation of the turbulent flow using a filtering procedure on Navier-Stokes equations, ultimately solving the equation for the average velocities. Averaging the Navier-Stokes equations affects the reliability of the solution, leading to model errors and uncertainty in the model predictions. This paper introduces the sensitivity analysis method as a means for assessing the uncertainty of the applied LES model. Assessing model error is only one aspect of the overall framework of error and uncertainty quantification; see [1] for a recent analysis which accounts for errors of various types within a much more formal framework. A paper published recently in this journal [2] shows that a certain method for introducing randomness into systems modeled using differential equations and their subsequent compuational simulations can significantly influence the power spectrum associated with realizations of the random fields. In particular, a parameter distinguishing the type of noise introduced into the model (white, pink, or brown noise) directly affects statistical properties of the output. Hence parametric uncertainty can be introduced into the process during the computational simulation of the system. Moreover, when design algorithms for large-scale engineering models are merged into the framework of uncertainty quantification, a variety of different stochastic variables is encountered and the uncertainty inherent in those parameters must be propagated through deterministic models [3]. Design sensitivity computation is an important element of these algorithms, and the ability to compute analytic moments and their design sensitivities is an important feature that is currently being explored [4]. The focus of this paper is the construction of an algorithm for efficient and accurate computation of sensitivity information for a particular LES model.

Generally speaking, sensitivity analysis of a physical system is the computation of derivatives of its state variables with respect to parameters upon which the response of the system explicitly and/or implicitly depends. By this definition, there are basically two main approaches for numerically approximating the sensitivities. One is using finite differences and the other is to form an equation for the designated sensitivity and then numerically solve it. The latter is called the sensitivity equation method or SEM. The SEM is classified into two different methods, the continuous sensitivity equation method or CSEM and the automatic differentiation method or ADM. This categorization follows from the issue of obtaining the discrete sensitivity equation either by first differentiating the state equation and then discretizing both the state and the sensitivity equations or by first discretizing the state equation and then differentiating the discretized equation to obtain a discrete sensitivity equation. The CSEM corresponds to the differentiate-then-discretize approach, and the ADM corresponds to the discretize-then-differentiate algorithm. In ADM, the discrete sensitivity system is obtained using automatic differentiation software, and this approach has been well-studied and has seen significant developments in the last 20 years; see [5–10] for just a few of the many examples in the literature. CSEM and ADM are not completely mutually exclusive techniques. The two can be combined, and this is illustrated in [11].

Sensitivities computed by CSEM or ADM are not the same, but in both approaches, as grid sizes go to zero the discrete sensitivities converge to the same entity in most cases. There are examples where the sensitivity function is less smooth than the original state variable. If this is not taken into account, each approach can fail to yield the correct sensitivity approximation; see [12] for a simple one-dimensional (1D) example to illustrate this issue. To study advantages and disadvantages of these two strategies in detail, please see [13]. Both CSEM and ADM have been used successfully in many cases, and the approach one may choose depends on which one is easier to handle in the context of its disadvantages. A recent paper by Li et al. [14] demonstrates the use of AD and adjoint method techniques to augment sampling methods with derivative information for efficient estimation of system response under parametric uncertainty. Although the computation of the gradient of the system response may be relatively expensive when compared to the cost of one system simulation for a given choice of parameters, the computational cost is also independent of the dimension of the parameter space. For model simulations with a large number of system parameters, using this approach holds a clear advantage over the CSEM. The paper also shows that obtaining derivative information using that approach must address issues of ill-conditioning in the polynomial regression, and differentiability with respect to uncertainty quantifiers is assumed for the theoretical results.

When one is utilizing a flow solver code, a finite difference quotient approximation to a sensitivity is simple to implement. However, it may not be an efficient method for computing sensitivities (see for example [15]). It produces large errors in addition to being computationally expensive in the sense that the code used for calculating a nonlinear flow has to be run for two different parameter inputs at the very least (see surveys [13, 15]).

In computing the flow sensitivity via CSEM, once the flow solution is obtained only a linear equation needs to be solved, which can be done using the same program as the one used for approximating the flow. Often the flow sensitivity solution can be done simultaneously with the flow solve or as one extra linear solve once the flow solver has converged. Hence, the sensitivity computation is obtained for a fraction of the computational cost of the original flow solve. Therefore, the use of CSEM is preferable to the use of the finite difference method, and one goal of this paper is to illustrate that extra computational effort expended in the sensitivity solve is far outweighed by the additional information that the sensitivity information provides the model practitioner. A comparison between these two methods in calculating sensitivity has been presented for a specific forebody design problem by Borggaard et al. [15] and for the specific model problem in this paper by the author in [16]. CSEM has been used to compute the sensitivity of flows with respect to different flow-related parameters. Much work has been done on this by Borggaard, Godfrey, and others (see the surveys [17–22]).

This paper focuses on the sensitivity of a subgrid eddy viscosity model or EVM with respect to the variation of the eddy viscosity parameter using CSEM. If we consider the eddy viscosity parameter as the filter length scale in the fluid model (see [16]), then sensitivity analysis can be used as a tool to compute the flow solution uncertainty due to the choice of the filter length scale. One subgrid EVM idea was introduced by Layton in [23]. In [24], the authors present the analysis and numerical computations of two first-order semi-implicit schemes for EVM as well as Navier-Stokes equations. An error analysis of EVM using discontinuous polynomial approximations can also be studied in [25].

The subgrid EVM in a bounded, simply connected two- or three-dimensional domain Ω with polygonal boundary ∂Ω is outlined as follows:

(1.1)

Here w : Ω × (0, T] → ℝd, d = 2 or 3, represents the resulting approximation of the large eddy velocity of a viscous incompressible fluid, p : Ω × (0, T] → ℝ the large eddy pressure, , the complement of filtering operator, is considered as Pʹ = IP, where P : L2(Ω) → L is an L2-orthogonal projection, defined on a chosen subspace of L2(Ω) [23], denoted by L, and f, the external force, is in L2[0,T,L2(Ω)]. In this model ν > 0 is the kinematic viscosity, which is inversely proportional to the Reynolds number Re, and α, the eddy viscosity parameter, corresponds to the filter length scale in LES models. Therefore its values vary between 0 and 1 with α = 0 corresponding to the Navier-Stokes equations.

Because different values of α in (1.1) cause different responses of the flow, it is natural to explore the sensitivity of the flow system and also the uncertainty of some computational fluid dynamics predictions which can be affected by changing the flow solution with respect to the variation of this parameter.

Generally, the EVM takes the form of finding the approximate large-scale velocity w and pressure p satisfying ∇ · w = 0 and

(1.2)

where ∇sw = 1/2 (∇w + ∇tw), and νT is the eddy viscosity coefficient that must be specified to select the model. One of the most commonly used EVMs is the Smagorinsky model that is given by the eddy viscosity choice of

(1.3)

Once the Smagorinsky constant Cs ∈ [0.01, 0.1] and the initial boundary conditions are specified, for a given α Eqs. (1.2) and (1.3) determine a solution (w, p) implicitly as a function of α displaying sensitivity with respect to the variations of this parameter. Sensitivity discussion of (1.2) is given in [26]. Other types of EVM are presented and analyzed in [27–29].

The contribution of this paper is twofold. First the sensitivity equation of EVM introduced by Eq. (1.1) with respect to the eddy viscosity parameter is derived. The complete numerical analysis of this equation using a semiimplicit discretization motivated by the work of authors in [24] is provided. Second, the sensitivity information is used to quantify the uncertainty of the EVM solution due to the variations of α. Furthermore, the paper extends the application of the sensitivity computation in improving the predictions of flow functionals.

The paper is organized as follows. Section 2 contains some basic notation and mathematical preliminaries. In Section 3, we derive the sensitivity equation of EVM with respect to the eddy viscosity parameter α. Then we illustrate the discretizations and numerical analysis of both EVM and its corresponding sensitivity equation in Sections 4 and 5. The last section is devoted to numerically demonstrating the use of sensitivity computation in uncertainty analysis of the flow solutions themselves as well as in the improvement of drag estimation. Two standard benchmark problems provide the illustrative examples.

2. PRELIMINARY NOTATION AND FUNCTION SPACES

In this section we briefly introduce some definitions and notations which will be used often throughout this work. For a comprehensive presentation of each concept we refer the interested reader to [30] and [31]. The L2(Ω) inner product, norm, and the seminorm of vectors and tensors are presented by (·, ·), ‖ · ‖, and | · |, respectively. For the Hilbert space ℋr(Ω), the seminorm is denoted by | · |r.

Following is the notation for function spaces that provide us with the natural environment for the variational theory of partial differential equations,

For any triangulation Th of Ω and any positive integer r, if ℙr denotes the set of polynomials of degree r, the finite dimensional subspaces are

The following definitions of norms and function spaces are provided for simplicity in notation.

Definition 2.1.

For Ω ⊆ ℝm the (a, b) weighted norm of a function u : Ω → ℝ is defined by

Definition 2.2.

For p, q ∈ [1, ∞), let for any t ∈ [0, T], u(t) ∈ Lp(Ω). Then for any positive integer N, define

where Δt = T/N, and the function spaces are defined by

For the treatment of the convective term in Eq. (1.1), we remind the reader that the following trilinear form

is well-defined and continuous on these spaces and it is skew-symmetric in its last two arguments. As shown in [30], if u, ∇u, wL2(Ω), ∇ · u = 0 and u · = 0 on ∂Ω, then

In the analysis, we often use the following inequalities:

  • Young′s inequality: for any real numbers a, b, and ϵ > 0:

    with p = q = 2 or p = 4/3 and q = 4.

  • Generalized Hölder′s inequality: Let uLp(Ω), vLq(Ω), and wLr(Ω) so that 1 ≤ p, q, r ≤ 1, and (1/p) + (1/q) + (1/r) = 1. Then
  • Poincare-Friedrich′s inequality: If Ω is bounded and connected, then there is a constant C such that

It is well known that Gronwall′s inequality plays an important role in the study of differential systems of various kind. Here we state the discrete version of Gronwall′s inequality (see [32]).

Lemma 1. If y(k), a(k) : Z+ → ℝ+, C ≥ 0 satisfying the condition y(0) ≤ C and

then

Overall in this paper, C denotes a generic positive constant and is independent of the mesh size, the viscosity ν, and the eddy viscosity parameter α.

It is known from [33] that if the right-hand side function ƒ in Navier-Stokes equations (1.1) with α = 0, belongs to L2(0, T; (ℋ10) ⃰ ) where (ℋ10) ⃰  denotes the dual space of ℋ10, there exists a solution u such that uL2(0, T; ℋ10) ∩ L(0, T; ℋ10). We assume in addition that utL(0, T; ℋ10) and for any integer r ≥ 1,

(2.1)
(2.2)
(2.3)

By assumptions (2.1)–(2.3), we aim to find the error estimates that are of optimal order with respect to the mesh size h, as h → 0. Error estimates of (r + 1)st-order based on (2.1)–(2.3) are obtained if the selected finite element spaces for the discrete velocity and pressure possess (r + 1)st-order approximability properties, i.e., those consisting of piecewise polynomial functions of degree less than r and r − 1, respectively (see [34, 35]).

3. THE CONTINUOUS SENSITIVITY EQUATION OF EVM

Various responses of the turbulence models to the users-elected model parameters have naturally led the scientists to study the sensitivity of models with respect to the variations of those parameters. For example, the EVM in (1.1) displays a sensitivity to the different values of parameter α. We apply the CSEM as a technique for such study and in this section we form a system of the continuous equations for the sensitivity by implicit differentiation of (1.1), with respect to α.

Definition 3.1.

Let (w, p) be the solution of (1.1). The sensitivity of (w, p) to variations of α is defined to be the derivative of (w, p) with respect to α.

Assuming that the L2-orthogonal projection P is differentiable with respect to parameter α, since this operator is a linear operator, using the chain rule it is easy to show that the operator P and then commute with differentiation with respect to α. Therefore, sensitivity of the solution (w, p) of the system (1.1) can be computed from the following sensitivity equation:

(3.1)

As it can be observed in (3.1), w appears in the sensitivity equation. Therefore, to complete the sensitivity analysis we need to couple (3.1) with (1.1). Sensitivities in LES models have been introduced and discussed in [26]. The algorithm of computing sensitivity from (3.1) is presented in [16].

4. THE DISCRETIZATIONS

This section is devoted to introducing the basis for deriving a finite element approximation of (1.1) and (3.1). We first assume that the three basic aspects of the finite element method, i.e., the existence of a triangulation of Ω, the construction of a finite dimensional subspace consisting of piecewise polynomials, and the existence of a basis of functions having small support, hold. Then we apply the following classical steps:

  • Variational Formulation: Equations (1.1) and (3.1) are reformulated in a weak form after multiplication by a suitable set of test functions, vX and λ ∈ Q, and performing an integration upon the domain. At this stage, the integration by parts is used to reduce the order of differentiation for solutions, w and s.
  • Discretization in Space: Let h ∈ (0; 1], tending to zero, be the mesh spacing for finite element, then V h = {v hX h : (λh, ∇ · v h) = 0, for all λhQ h} is a finite dimensional subspace of X h. In general V h is not a subspace of V, but since (V h, Q h) fulfills the inf-sup or Babuska-Brezzi stability condition, the pressure p h can be eliminated from the system in its discrete form (see [30]). By a similar analysis, its discrete sensitivity q h is removed from the sensitivity equation. In the resulting equations, for each t ∈ [0,T], w h and s h are solved in V h. In addition, the nature of the L2-orthogonal projection P leads us to consider a multiscale spatial discretization for the model as follows. Let h and H denote two different mesh widths (h < H). Let LHL2(Ω)d×d; d = 2 or 3, be a finite element space. A discrete analog PH : L2(Ω)d×dLH of the L2-orthogonal projection P is defined by.
    (4.1)

    Note that PH is determined by the choice of LH. The space LH is considered as the space of large scales of the velocity that are numerically solved by EVM since H represents the coarse mesh size. By the properties of P as an L2-orthogonal projection, it is easy to show that H (∇ w h) = (∇ w h)|LH. Therefore, for simplicity in notations we use (∇ w h) instead of H (∇ w h).

  • Discretization in Time: We start with partitioning the time interval [0,T] into N subintervals [tn, tn+1] of length Δt = T/N, with t0 = 0 and tN = T. Then, at each time level tn, an approximation to w and s, denoted by w hn and s hn, respectively, are obtained. To achieve a full discretization of (1.1) and (3.1), we use a first-order implicitexplicit time-stepping scheme. This scheme was first introduced and analyzed for a convection-diffusion type of equation in [36]. A complete analysis of this scheme along with another semi-implicit scheme for EVM can be studied in [24]. The method uses an implicit structure on stabilizing terms and an explicit one on unstabilizing terms. The discretization leads to linearizing w ⋅ ∇w by wn ⋅ ∇wn+1 in Eq. (1.1) and treating P(∇w) in (1.1) and also P(∇s) in (3.1) explicitly.

Note that the time discretizations of the EVM and sensitivity equations are of a basic Euler approach. The details of each step of the discretization are described below.

The variational formulations of EVM and CSEM in X and Q are, respectively, given as

(4.2)

for all vX, λ ∈ Q and

(4.3)

for all vX, λ ∈ Q.

The discretized variational formulation of (4.2) and (4.3) in space, respectively, read as follows:

(4.4)

and

(4.5)

for all v hX h, λhQ h. Let v hV h, then Eqs. (4.4) and (4.5) can be rewritten equivalently as follows:

(4.6)
(4.7)

The fully discrete form of (4.6) and (4.7) reads as follows: Given w hn, we seek w hn+1 satisfying

(4.8)

Similarly, given s hn, we find s hn+1 such that

(4.9)

The algorithm for numerically calculating the L2-projection terms in (4.8) and (4.9), i.e., P(∇w hn ) and P(∇s hn ), is presented in [16]. Next, we study the stability and the convergence of the approximated large eddy velocity w and its sensitivity s resulting from (4.8) and (4.9), respectively.

5. THE STABILITY AND ERROR ANALYSIS

This section illustrates the stability and error estimate of the approximated eddy viscosity and its sensitivity in (4.8) and (4.9). The stability and error bounds, and the rate of convergence for EVM in (4.8) is stated in Proposition 5.1, Theorem 3, and Corollary 5. The proofs of these results are similar to the ones for the sensitivity equation in (4.9) and can be also found in [24]. Therefore the authors present the proofs for the sensitivity equation that are more technical.

Our error analysis requires us to assume that for a constant C, the finite dimensional space LH satisfies an inverse inequality of the form found in [23] and given by

(5.1)

The inverse inequality (5.1) is satisfied if the projection space LH is chosen so that LH = ∇XH. As an example, consider selecting the function spaces X h and LH to be continuous piecewise linear polynomials and L2-piecewise constant functions, respectively.

In addition, we list an equality and some inequalities that are used throughout the given proofs in the following lemma.

Lemma 2. Let u h, v hX h, then the following equality and inequalities hold:

(5.2)
(5.3)
(5.4)
(5.5)

Proof. The equality (5.2) is obtained using the fact that P and are complement projections, i.e., P + = I. The inequality (5.3) is an immediate conclusion of (5.2). Inequalities (5.4) and (5.5) are basic inequalities in algebra.

Proposition 5.1.

For any positive integer N, the solution w hN of (4.8) is bounded and

Proof. See [24].

The proofs of Proposition 5.2 and Theorem 4 are largely technical; therefore we give in advance a brief sketch of those proofs step by step.

  • Step 1: In variational formulation (4.9), move the convective and the projection terms to the right side.
  • Step 2: Choose the test function w h = s hn+1 and modify the convective term as following:
    (5.6)
  • Step 3: Apply the Cauchy-Schwarz inequality then Young′s inequality on each term of the right-hand side.
  • Step 4: Sum the inequality over the time and apply Gronwall′s Lemma 1.

Proposition 5.2.

Given any strong solution w of (1.1) and its approximated solution w h of (4.8) on the interval [0,T] such that ∇(ww h) ∈ l4 (L2), for any positive integer N, the sensitivity solution s hN of (4.9) is stable. Specifically,

(5.7)

Proof. First take the convective and the projection terms in (4.9) to the other side of the equation. Then set v h = s hn+1, replace b(s hn+1, w hn+1, s hn+1) according to (5.6), and apply the inequality (5.4) to the first term on the left side to get,

(5.8)

To find the stability bound, each term on the right side of (5.8) should be bounded. The standard method is to apply the Cauchy-Schwartz inequality followed by the Young′s inequality. This implies

Sum the above inequality from 0 to N, then multiply by 2Δt to get

(5.9)

Using the inequality (5.3) on the right-hand side of (5.9) and by Gronwall′s inequality, one finds

where C ⃰ = exp{1/2(6/ν)3C(|ww h|4l4(L2) + T|w|4L(L4))}. The bound in (5.7) is obtained simply by replacing the stability bound for w h from Proposition 5.1 in the last inequality.

Theorem 3. Let (w, p) be a strong solution of (1.1) on the interval [0,T]. Suppose

Let en = wnw hn denote the global error. Assume that e0 = 0. Then, the approximated solution w hn is convergent and there is a constant C ⃰  independent of α, h, and H such that

(5.10)

where

Proof. See [24].

Following is the algorithm of the proof of Theorem 4 on sensitivity error estimation.

  • Step 1: Subtract the variational formulation (4.9) from (4.3) to obtain the error equation (5.14). Notice that the local truncation error is in the following form:
    (5.11)
  • Step 2: Take (, ) ∈ (V h, Q h) to be any interpolation and split the error e sn = snn − (s hnn) = ηn − φ hn. Set v h = φ hn+1 and λ h = n+1 and replace P by IP′ to get (5.15).
  • Step 3: Use (5.2), (5.4) and (5.5). Rewrite the trilinear terms as follows:
    (5.12)
  • Step 4: On the right side of (5.16), use inequalities (2.1)–(2.3) along with Cauchy-Schwarz inequality, then apply Young′s inequality.
  • Step 5: Add up over time steps and apply Gronwall′s Lemma 1 to (5.18) to obtain an error estimate for φ h.
  • Step 6: Compute the error estimate ss h by triangle inequality. The terms containing ´ are bounded by the interpolation error estimate of the finite element spaces.

Theorem 4. Suppose that the assumption of 3 holds and in addition

as the mesh size h → 0. Let (s, q) be a sensitivity strong solution such that

Let e sn = sns hn denote the global error. Assume that e s0 = 0. Then, the approximated sensitivity s hn is convergent and there is a constant C ⃰  independent of α, h, and H such that

(5.13)

where

Proof. Subtracting (4.9) from (4.3) yields the following error equation:

(5.14)

for all v hV h, λhQ h. For any interpolation , the error is decomposed into e sn = snn − (s hnn) = ηn − φ hn. Take v h = φhn+1, and λh = n+1 in (5.14). Replace P by IP′ in the following terms: P(wt,n), P(wttt,n), P(wtt,n), and P(∇v h). The resulting equation is

(5.15)

Applying (5.2) and (5.4) to the left side and (5.5) to the first term on the right side of (5.15), then rewriting the trilinear terms according to (5.12), gives

(5.16)

Technically, we want to bound all the terms on the right-hand side of (5.16). This is done first by the Cauchy-Schwartz inequality along with inequalities (2.1) to (2.3), then by applying Young′s inequality (see the bounds in the Appendix). Inserting all the bounds into (5.16) yields

(5.17)

In the next step, multiplying by Δt and summing from 0 to N gives

(5.18)

Define

and let

The result is obtained by applying Gronwall′s inequality and the fact that ⫼e2 ≤ ⫼η⫼2 + ⫼φh2.

Corollary 5. The optimal rate of convergence for the left-hand side of (5.10) is obtained with the choice of 0 ≤ α = ha ≤ Δt, with a > 0, and hH = hb < 1 for b < 1. The error estimate becomes O(hr + Δt).

Proof. See [24].

Corollary 6. Under the assumption of Corollary 5, the rate of convergence for the left-hand side of (5.13) is O(hr + ha−bΔt + Δt), where ba.

Proof. To find the rate of convergence for the error estimate in (5.13), we compute the rate of convergence for ℬ1(ν, α, H); ℬ2(ν, α, H), and ℬ3(ν, e). By Corollary 5, given 0 ≤ α = ha ≤ Δt, with a > 0, and hH = hb < 1 for b < 1, ℬ3(ν, e) is O(hr + Δt).

Let ℬ1(ν, α, H), and ℬ2(ν, α, H) be of order hk1 and hk2Δtk3tk4, respectively, with k1, k2, k3, k4 ≥ 0. The coefficients of ℬ2(ν, α, H) imply that k3 = k4 = 0 and balancing α and H in α2H−2 = h2(a−b) enforces a condition on selecting a and b such that ba. Thus ℬ2(ν, α, H) is of order ha−b.

In ℬ1(ν, α, H), αH−2h2 and α2H−2 are the key coefficients. Considering αH−2h2 = ha−2b+2 and α2H−2 = h2(a−b), they are equal if a = 2. For a < 2, α2H−2 = h2(a−b) has the smallest power and for a > 2, αH−2h2 = ha−2b+2. In all the cases, ℬ1(ν, α, H) yields an order of hk1 with k ≥ 0. Therefore, the right-hand side of (5.13) converges with a least rate of O(hr + ha−bΔt + Δt).

Remark 7. An important conclusion of Theorem 4 is that the convergence of approximated sensitivity depends on the convergence of approximated eddy viscosity velocity wh. Therefore, an accurate sensitivity requires the sensitivity computations to be carried out using a finer mesh. This is a computational challenge when obtaining a good velocity approximation becomes difficult for small values of viscosity. In [26], the authors show this fact using a different approach by analyzing the boundary layers. Achieving the optimal rate of convergence for the sensitivity as stated in Corollary 6 may not be computationally feasible in cases of very small viscosities.

6. NUMERICAL RESULTS

In this section, our numerical studies illustrate the convergence results for the fully discrete sensitivity equation in (4.9) and two applications of the sensitivity computation for the subgrid EVM given by (1.1), one in assessing the uncertainty of the approximated solution due to the variation of the eddy viscosity parameter and the second one in improving the flow functionals. In these studies, we have used the two-dimensional cavity problem and the twodimensional flow around a cylinder as our test problems.

Section 6.1 presents numerical studies in support of the error estimates and the convergence rate stated in Theorem 4 and Corollary 6. The experiment contained in Section 6.2 provides an example for numerically specifying an interval of reliability which indicates a range of values of the eddy viscosity parameter that yield the most reliable approximation for a given problem. The computations in Section 6.3 present a more accurate estimation of the drag values in the classical example of flow around a cylinder by using the flow variable sensitivities as the first-order correction terms in the Taylor expansion of flow variables.

Numerical results contained in this paper are obtained by choosing the finite element spaces X h, Q h, and LH to be spaces of piecewise polynomials of degree 2, piecewise linear polynomials, and piecewise constant functions, respectively. This set of finite element spaces ensures the stability of the computed L2-orthogonal projection P from the system of (4.1) at each time step. Furthermore, the inequality (5.1) holds for P. All calculations presented in this paper were carried out using the software package FreeFEM++; see [37] for details and examples.

6.1 Convergence Results

Convergence results for the numerical scheme in (4.9) are given in this section. These results are shown using the two-dimensional cavity problem for the spatial convergence and the two-dimensional flow in a channel with cylinder for the convergence in time. In [24], the authors provide the analysis and numerical computations of the scheme in (4.8) for the approximated velocity w h.

6.1.1 Convergence in Space

The numerical experiment is performed on the two-dimensional cavity problem [38]; see also [16] for a comparison of the sensitivity computations for the cavity problem via two different strategies.We numerically solve Eq. (1.1) with the domain defined by Ω = [0, 1]×[0, 1]. The upper boundary moves with the velocity w(t, x, y) = [16x2(1−x)2, 0]T. The initial condition is chosen to be w(0, x, y) = (3y2 − 2y, 0)T in Ω. It is clear that since the initial and boundary conditions for w do not depend on α, we have zero initial and boundary condition for sensitivity s.

The sensitivity sh is computed at each time step ti using the scheme in (4.9) through computing w h(ti+1) from (4.8). We use the numerical approximation obtained with a mesh size h ⃰  = 1/49 and a coarse mesh size H ⃰  = 1/7 as the sensitivity true solution of (4.9). For the error computations, a sequence of grids is generated with the structure (h, H)k = (1/k, 1/√k), for k = 3, 4, 5, 6. These computations are obtained using a fixed eddy viscosity parameter of α = 0,05 and the uniform time step Δt = 0.001 for 1000 steps. Given a grid size, denoted by h, the error in the calculation of {s hn } Nn=1 is defined as follows:

Selecting different values for ν = 1, 0.5, 0.25, 0.125 in (4.8) and (4.9), the error computations for ‖e sL(L2) are given in Fig. 1. According to data shown in this figure, ‖e sL(L2) decreases with the mesh refinement for all ν, but its value is larger for smaller values of viscosity.

FIG. 1: Sensitivity spatial error and its log-log plot.

The convergence rate of the error is determined by calculating the exponent a in the expression

where h1 = 1/25 and h2 = 1/36. Mathematically, this is equivalent to measuring the slope of the line connecting the leftmost two points in the log-log error plots shown in Fig. 1. The rate coefficient a for different viscosity values are listed in Table 1. Note that in this experiment α = 0.05 ≈ h0.835975 and h = h0.5. Thus, with our choice of spaces according to Corollary 6, we expect a rate of convergence between 0.335975 and 2 for all ν.

TABLE 1: Rate of convergence for different values of ν

ν Rate of convergence a
1 1.253
0.5 1.196
0.25 1.280
0.125 1.549

6.1.2 Convergence in Time

For this numerical experiment, we consider the standard test problem of two-dimensional flow in a channel around a cylinder. The geometry is indicated in Fig. 2. The channel height and width are, respectively, 0.41 and 2.2 m. The cylinder, denoted by D, has radius 0.05 m, and it has been placed in the geometry so that its center is given by the coordinates (0.2, 0.2).

FIG. 2: Geometry of two-dimensional flow around cylinder.

We seek to generate a numerical approximation to the solution of (1.1) where for 0 ≤ t ≤ 1, the inflow conditions are given below and represent a parabolic inflow that is periodic in time:

(6.1)

A free condition is used for the outflow boundary condition, and the remaining boundary and initial conditions are given by

For the triangulation of the domain in Fig. 2, note that each mesh is nonuniform in the sense that the mesh is finer around the cylinder D. A given mesh is constructed using a maximum size, denoted by hmax, for the sides of the channel and a minimum size, denoted by hmin, for the boundary of D. We identify these maximum and minimum sizes of the mesh using the ordered pair (hmax, hmin). For the numerical computation of the projection operator, the applied coarse mesh has the same structure and is always chosen as (√hmax, √hmin). Figure 3 displays an example of a mesh of size (1/36, 1/49).

FIG. 3: Mesh in a channel of size (1/36, 1/49).

In this experiment, the spatial grid size is set to be (1/49, 1/64). The computations are done using uniform time steps Δt = 0.1, 0.02, 0.01, 0.002 and the errors are compared to the fixed time step Δt = 0.001 for the final time T = 1. The convergence in time is tested for different values of viscosity ν = 1, 0.1, 0.01, 0.001 with α = 0.025. The sensitivity error in l2(L2) for the scheme in (4.9) along with its corresponding log-log plot are displayed in Fig. 4. This figure shows that the sensitivity components for all viscosity values converge as the time steps are chosen smaller and a convergence rate of about 0.8 is obtained for ν < 1.

FIG. 4: Sensitivity error in time and its log-log plot.

6.2 The Interval of Reliability

The aim of this numerical study is to show that the flow sensitivity obtained from (4.9) can be used to quantify the reliability of a flow solution computed using (4.8) for different values of α. The idea is simply based on the following difference quotient for the sensitivity:

(6.2)

Considering w as an implicit function of α, w(0) indicates the true solution of Navier-Stokes equations while w(α) for α > 0 denotes the corresponding EVM approximation. According to (6.2), the EVM solution is an accurate approximation to the Navier-Stokes solution when ‖w(α) − w(0)‖ is small, and the accuracy of the EVM approximation can be estimated by measuring α‖s‖. As noted in Section 1, the sensitivity calculation can be coupled with that of the original EVM simulation. Hence, one can compute a model simulation with a given set of parameters and, for a nominal extra cost, the sensitivity computation can be obtained and a quantitative measure of reliability can be computed. This is illustrated for a test problem in the following paragraphs.

The test problem used in this numerical study is the two-dimensional cavity problem defined on a unit square with the same boundary and initial conditions as described in Section 6.1.1. All computations are carried out with a fixed mesh size h = 1/36 and coarse mesh size H = 1/6, and the uniform step Δt = 0.001 for 1000 steps. Table 2 presents α‖sL∞(0,T; L2) values for variations of the parameter α and different flow viscosities corresponding to Reynolds numbers of 1000, 5000, 10,000, and 50,000. This table demonstrates that the large eddy velocity w is highly sensitive with respect to the small values of α, i.e., α ≤ 0.025, for all the tested Reynolds numbers and there is a large increase in the sensitivity values as α takes on values closer to 0 (see Fig. 5).

TABLE 2: Values of α‖sL∞(0,T; L2) at the final time

α Re = 1000 Re = 5000 Re = 10, 000 Re = 50, 000
0.75 0.085492 0.0855515 0.0855589 0.0855649
0.5 0.0815713 0.0816575 0.0816682 0.0816769
0.25 0.0773889 0.0775531 0.0775736 0.0775901
0.075 0.0774764 0.0780792 0.0781552 0.0782162
0.05 0.0866153 0.08830347 0.0882171 0.0883635
0.025 3.03762 3.63664 3.72093 3.79002
0.0075 4355 8156.78 8867.12 9487.9
0.005 31871.8 76437.2 85985.1 94611.4

FIG. 5: Sensitivity at t = 0.01 and 0.1, respectively.

Note that by selecting large values for the eddy viscosity parameter α, e.g., larger or equal to 0.25, all the velocity scales that are less than or equal to the parameter α are removed from the flow structure. Hence, an approximated EVM flow solution corresponding to large values of α is not considered to be a reliable approximation to a solution of the Navier-Stokes model because too much of the small-scale structure has been lost. This situation is especially tenuous for the case of high Reynolds numbers where the velocity contains a large number of small scales. Therefore, the practitioner must seek to find the optimal balance between choosing a value of α that is “small enough” to provide a reliable approximation to the Navier-Stokes flow while choosing a value of α that is “large enough” so that the computation of the large eddy velocity w is relatively insensitive to small changes in α. That is, we seek to identify a range of values for the interval of reliability so that both α and α‖s‖ are small. To determine the interval of reliability, one can use (6.2) along with the Taylor expansion,

(6.3)

One first obtains a number of flow and sensitivity approximations for a variety of values of α as seen in Table 2. For the calculations of the current example, a reliable solution is obtained for α values when α‖s‖ from (6.2) is minimized and O2) in (6.3) is of order less than 0.01.

Table 2 suggests that an approximate interval of reliability for this model problem is [0.05,0.075]. That is, for α in the interval [0.05,0.075], the approximated flow response using the EVM (4.8) is a reliable approximation to the Navier-Stokes flow for the given Reynolds numbers in this experiment.

6.3 Improving Flow Functionals

Let J(w(0)) = J(u) be a flow functional. Assume that w(0) is extremely computationally expensive to obtain directly, and one would like to compute a less expensive approximation. Then, expanding J(u) around a nonzero α implies that

(6.4)

Notice that for linear functionals, J′ = J; therefore, the approximation (6.4) is rewritten as

(6.5)

This idea was proposed by Anitescu and Layton for LES models and was tested on the Smagorinsky model in [26].

The lift and drag functional for Navier-Stokes equations is given by

(6.6)

where denotes the normal vector on the cylinder boundary D directing into the channel, ∇su presents the deformation tensor and is 1/2(∇u + ∇uT ), the unit vector â in the positive direction of x axis or negative direction of y axis yields the drag or lift flow functional.

The computations in this experiment are carried out using the two-dimensional flow around a cylinder given in Section 6.1.2. The reference value of drag for this test problem is obtained by performing the direct numerical solution (DNS) method to a uniform fine mesh of size (1/100, 1/121) for 0 ≤ t ≤ 4.

Table 3 indicates the reference values of maximum drag, and the error in its estimation using the approximated large eddy velocity and pressure (w, p) and (w − αs, p − αq) in (6.6) for different values of ν. In this experiment, the approximated flow variables and their sensitivities are obtained from (4.8) and (4.9) with α = 0.00125 and a mesh size of (1/49, 1/64). This table presents more accurate drag values by (w − αs, p − αq) for all ν, especially for ν ≤ 0.01.

TABLE 3: Maximum drag values and the errors

Max. drag Error using Error using
ν(Re) J(u, p) |J(u, p) − J(w, p)| |J(u, p) − J(w − αs, p − αq)|
1 (1) 63.7703 0.4037 0.3702
0.1 (10) 41.1958 0.3628 0.3555
0.01 (100) 36.0677 0.389 0.0152
0.001 (1000) 35.29035 0.28095 0.02095
0.0001 (10,000) 35.1186 0.4354 0.0154

Examining the norm of the sensitivity quantities for a range of values of ν in Table 4, one observes that the norm of the sensitivity is negligible for ν ≥ 0.1. Therefore, the flow is insensitive to small changes in ν over that range of values. This is also reflected in Table 3, as we see that the computed drag values using (w − αs, p − αq) show only a small improvement in comparison to the ones computed using (w, p) when ν ≥ 0.1. However, for ν ≤ 0.01, the errors incurred by using (w − αs, p − αq) improve (or decrease) by a full order of magnitude. According to Table 4, for small values of ν, i.e., ν ≤ 0.01, the flow becomes more sensitive, and using sensitivity information improves the estimated values of the drag functional significantly.

TABLE 4: Sensitivity for different values of ν

ν(Re) α‖sL∞(0,T; L2)
1 (1) 7.19057 × 10−6
0.1 (10) 2.88244 × 10−4
0.01 (100) 0.00483735
0.001 (1000) 0.0155576
0.0001 (10,000) 0.0201101

Note that this type of information can also be obtained graphically, and the reduction of error that one sees by using sensitivity information error can fluctuate depending on the boundary and initial conditions as well as the time at which one examines the errors. For this particular model problem, we present a few graphical representations of error that exhibit such fluctuations. A comparison between the error estimates of drag values using (w, p) and (w − αs, p − αq) at time T = 1, 2, 3, and 4 for different values of ν is shown in Figs. 6 and 7. For the particular boundary and initial conditions in this example, note that the maximum value of drag occurs at time T = 2, and the drag values at T = 4 are very small since the inflow vanishes at this specific time. Hence, the amount of improvement in flow functional estimation that one obtains using the sensitivity information can fluctuate over the range of parameter values as well as over the time interval of interest. In particular, reductions of error are the most significant where the sensitivities are larger than the norm. Here, that corresponds to ν ≤ 0.01. This is of course a crude estimation in the sense that we have only examined a small number of parameter values between 10−4 and 1. However, even the limited information provided by such an experiment provide beneficial insight as to when the EVM is a reliable model to use within the framework of an optimal design process.

FIG. 6: Drag error estimates at T = 1 and 2.

FIG. 7: Drag error estimates at T = 3 and 4.

Note that the sensitivity information used here can be extracted in a very inexpensive manner once the numerical method for (4.8) is constructed. The incorporation of the sensitivity equation in (4.9) into a numerical algorithm that computes both the eddy velocities and their sensitivities is very straightforward as the bulk of the work for the computation is in the implementation of Eq. (4.8). Once that is accomplished, computations for Eq. (4.9) are easily added since all of the data structures and projection calculations are virtually the same or very similar.

7. CONCLUDING REMARKS

In this paper, we have developed and analyzed a numerical method for simulating an eddy viscosity model and its corresponding sensitivity with respect to the eddy viscosity parameter. The method is a first-order implicit-explicit time-stepping method. The fully discrete approximations are constructed for velocity and its sensitivity utilizing the finite element method in space. We provided the optimal rate of convergence for these estimations using piecewise polynomial functions of degree r. Once the numerical algorithm for solving the EVM is implemented, the sensitivity calculations can be appended to it for a nominal amount of effort. Convergence studies for both time and space were included for standard test problems. Subsequent numerical experiments illustrate that the sensitivity calculations are quite useful in quantifying the uncertainty of model errors incurred by using the EVM as an approximation to the full-blown Navier-Stokes model. In addition, the sensitivity information is shown to provide increased accuracy in flow functionals for a relatively inexpensive cost of the sensitivity calculation. This increase in accuracy may also lead to a more efficient optimal design process.

REFERENCES

1. Liang, B. and Mahadevan, S., Error and uncertainty quantification and sensitivity analysis in mechanics computational models, Int. J. Uncertainty Quantif., 2:147–161, 2012.

2. Stoyanov, M., Gunzburger, M., and Burkardt, J., Pink noise and their effect on solutions of differential equations, Int. J. Uncertainty Quantif., 1:257–278, 2011.

3. Wojtkiewicz, S., Eldred, M., Field, R., Urbina, A., and Red-horse, J., Uncertainty quantification in large computational engineering models, In AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA Paper no. 2001-1455, 2001.

4. Eldred, M., Design under uncertainty employing stochastic expansion methods, Int. J. Uncertainty Quantif., 2:119–146, 2011.

5. Abdel-Khalik, H., Hovland, P., Lyons, A., Stover, T., and Utke, J., A low rank approach to automatic differentiation, Advances in automatic differentiation, Lect. Notes Comput. Sci. Eng., 64:55–65, 2008.

6. Averick, B., Moré, J., Bischof, C., Carle, A., and Griewank, A., Computing large sparse Jacobian matrices using automatic differentiation, SIAM J. Sci. Comput., 15:285–294, 1994.

7. Bischof, C. and Griewank, A., Computational differentiation and multidisciplinary design, Inverse problems and optimal design in industry (Philadelphia, PA, 1993), Eur. Consort. Math. Indust., 10:187–211, 1994.

8. Griewank, A., A mathematical view of automatic differentiation, Acta Num., 12:321–398, 2003.

9. Griewank, A. and Walther, A., Evaluating Derivatives, Principles and Techniques of Algorithmic Differentiation, SIAM, Philadelphia, 2008.

10. Hovland, P., Mohammadi, B., and Bischof, C., Automatic differentiation and Navier-Stokes computations, Computational methods for optimal design and control, Progr. Syst. Control Theory, 24:265–284, 1998.

11. Borggaard, J. and Verma, A., On efficient solutions to the continuous sensitivity equation using automatic differentiation, SIAM J. Sci. Comput., 22:39–62, 2000.

12. Davis, L. and Singler, J., Computational issues in sensitivity analysis for 1D interface problems, J. Computat. Math., 29:108–130, 2011.

13. Gunzburger, M., Perspectives in Flow Control and Optimization, SIAM, Philadelphia, 2003.

14. Li, Y., Anitescu, M., Roderick, O., and Hickernell, F., Orthogonal bases for polynomial regression with derivative information in uncertainty quantification, Int. J. Uncertainty Quantif., 1:97–112, 2011.

15. Borggaard, J., Burns, J., Cliff, E., and Gunzburger, M., Sensitivity calculations for a 2D, inviscid, supersonic forebody problem, In AMS-IMS-SIAM, Philadelphia, pp. 14–24, 1993.

16. Pahlevani, F., Sensitivity computations of eddy viscosity models with an application in drag computations, Int. J. Num. Methods Fluids, 52:381–392, 2006.

17. Borggaard, J. and Burns, J., A sensitivity equation approach to shape optimization in fluid flows, In IMA Conf. on Fluid Flow Control, Springer-Verlag, pp. 49–78, 1995.

18. Borggaard, J., Pelletier, D., and Turgeon, E., A continuous sensitivity equation method for flows with temperature dependent properties, In AIAA/USAF/NASA/ISSMO Symp. on Multidisciplinary Analysis and Design, AIAA Paper no. 2000–0881, 2000.

19. Borggaard, J., Pelletier, D., and Turgeon, E., Sensitivity and uncertainty analysis for variable property flows, In AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper no. 2001–0140, 2001.

20. Godfrey, A. and Cliff, E., Direct calculation of aerodynamic force derivatives: A sensitivity equation approach, In AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper no. 1998–0393, 1998.

21. Godfrey, A. and Cliff, E., Sensitivity equations for turbulent flows, In AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper no. 2001–1060, 2001.

22. Godfrey, A., Cliff, E., and Eppard,W., Using sensitivity equations for chemically reacting flows, In AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, AIAA Paper no. 1998–4805, 1998.

23. Layton,W., A connection between subgrid scale eddy viscosity and mixed methods, Appl. Math. Comput., 133:147–157, 2002.

24. Davis, L. and Pahlevani, F., Semi-implicit schemes for transient for navier-stokes equations and eddy viscosity models, J. Num. Methods Partial Differ. Eq., 25:212–231, 2009.

25. Kaya, S. and Riviere, B., Analysis of a discontinuous Galerkin and eddy viscosity method for Navier-Stokes equations, SIAM J. Num. Anal., 43:1572–1595, 2005.

26. Anitescu, M. and Layton, W., Sensitivities in large eddy simulation and improved estimates of turbulent flow functionals, SIAM J. Sci. Comput., 29:1650–1667, 2007.

27. John, V., Large Eddy Simulation of Turbulent Incompressible Flows. Analytical and Numerical Results for a Class of LES Models, Springer, Berlin, 2004.

28. Iliescu, T. and Layton, W., Approximating the large eddied in fluid motion III: The Boussinesq model for turbulent fluctuations, An. Stintf. Univ., 44:245–261, 1998.

29. Layton, W. and Lewandowski, R., Analysis of an eddy viscosity model for large eddy simulation of turbulent flows, J. Math. Fluid Mech., 2:374–399, 2002.

30. Girault, V. and Raviart, P., Finite Element Approximation of the Navier-Stokes Equations, Springer, Berlin, 1979.

31. Brenner, S. and Scott, R., The Mathematical Theory of Finite Element Methods, Springer, Berlin, 2002.

32. Niamsup, P. and Phat, V., Asymptotic stability of nonlinear control systems described by difference equations with multiple delays, J. Differ. Eq., 2000(11):1–17, 2000.

33. Temam, R., Navier-Stokes Equations and Nonlinear Functional Analysis, SIAM, Philadelphia, 1995.

34. Girault, V. and Scott, R., A quasi-local interpolation operator preserving the discrete divergence, Calcolo, 40:1–19, 2003.

35. Heywood, J. and Rannacher, R., Finite element approximation of the nonstationary Navier-Stokes problem III: Smoothing property and higher order estimates for spatial discretization, SIAM J. Num. Anal., 25:489–512, 1988.

36. Anitescu, M., Pahlevani, F., and Layton,W., Implicit for local effects and explicit for nonlocal effects is unconditionally stable, Electr. Trans. Num. Anal., 18:174–187, 2004.

37. Hecht, F., Pironneau, O., and Ohtsuka, K., Software FreeFem++, www.freefem.org, 2003.

38. Weinan, E. and Liu, J.-G., Vorticity boundary condition and related issues for finite difference schemes, J. Comput. Phys., 124:368–382, 1996.

39. Bardenhagen, S., Brydon, A., Williams, T., and Collet, C., Coupling grain scale and bulk mechanical response for PBXs using numerical simulations of real microstructures, Shock Compression of Condensed Matter, Am. Inst. Phys. Conf. Ser., 845:479–482, 2006.

40. Dolbow, J., Khaleel, M., and Mitchell, J., Multiscale Mathematics Initiative: A Roadmap, Prepared for the U.S. Department of Energy under contract DE-AC05-76RL01830, 2004.

41. Foo, J., Multi-Element Probabilistic Collocation in High Dimensions: Applications to Systems Biology and Physical Systems, Brown University, Division of Applied Mathematics, 2008.

42. Ghanem, R. and Spanos, P., Stochastic Finite Elements: A Spectral Approach, Springer, New York, 1991.

43. Klimke, A., Uncertainty Modeling Using Fuzzy Arithmetic and Sparse Grids, Universität Stuttgart, 2006.

44. Xiang, M. and Zabaras, N., An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations, J. Comput. Phys., 228(8):3084–3113, 2009.

APPENDIX A.

In the process of proving Theorem 4, assign the terms on the right-hand side of (5.16) by T1, T2, …, T16. Note that the local truncation error in (5.16) consists of three terms as shown in (5.11) and here T14, T15, and T16 are their corresponding terms. To bound these terms, we first use Cauchy-Schwarz inequality on all terms. In addition, assumptions (2.1)–(2.3) were used on T1, T2, T5, T8, T9, T10, T11 and T12, and the inequality (5.1) was applied in bounding T11, T15, and T16. The final bounds are obtained by applying the Young′s inequality.

Note that for the bound on T6, we assume that ‖eL∞(L2) |e|L∞(L2) can be bounded by a small number for a small h. The coefficients ai and bi for i = 1, 2, …, 16 are selected as follows.

TABLE A.5: Coefficients ai and bi

i ai bi
1 6/ν ν/12
2 6/ν ν/24
3 153/32ν3 ν/10
4 6/ν ν/24
5 6/ν ν/24
6 ν/24
7 6/ν ν/24
8 6/ν ν/24
9 5/2ν ν/10
10 α2(5/2ν) ν/10
11 α2(6/ν) ν/24
12 5/2ν ν/10
13 5/2ν ν/10
14 6/ν ν/24
15 α2(6/ν) ν/24
16 α2(6/ν) ν/24
Begell Digital Portal Begell Digital Library eBooks Journals References & Proceedings Research Collections Prices and Subscription Policies Begell House Contact Us Language English 中文 Русский Português German French Spain