The effect of uncertain geometries on advection–diffusion of scalar quantities

The two dimensional advection–diffusion equation in a stochastically varying geometry is considered. The varying domain is transformed into a fixed one and the numerical solution is computed using a high-order finite difference formulation on summation-by-parts form with weakly imposed boundary conditions. Statistics of the solution are computed non-intrusively using quadrature rules given by the probability density function of the random variable. As a quality control, we prove that the continuous problem is strongly well-posed, that the semi-discrete problem is strongly stable and verify the accuracy of the scheme. The technique is applied to a heat transfer problem in incompressible flow. Statistical properties such as confidence intervals and variance of the solution in terms of two functionals are computed and discussed. We show that there is a decreasing sensitivity to geometric uncertainty as we gradually lower the frequency and amplitude of the randomness. The results are less sensitive to variations in the correlation length of the geometry.


Introduction
When solving partial differential equations, uncertain geometry of the computational domain may arise for many reasons. Examples include irregular materials, inaccurate Computer-Aided Design (CAD) software, imprecise manufacturing machines and non-perfect mesh generators. We study the effects of this uncertainty and impose the boundary condition at stochastically varying positions in space. Related techniques are boundary perturbation [26], Lagrangian approach [1] and isoparametric mapping [5]. Other techniques dealing with geometric uncertainty include polynomial chaos with remeshing of geometry [8,9] as well as chaos collocation methods with fictious domains [3,17].
We transform the stochastically varying domain into a fixed one. This procedure has previously been used in Xiu et al. [25] for elliptic problems. Numerical techniques can be employed if the analytical transformation of the geometry is unavailable [4]. In this article it is extended to the analysis of the time-dependent advection-diffusion equation. The continuous problem is analyzed using the energy method, and strong well-posedness is proved [14,15].
We discretize using high-order finite difference methods on summation-by-parts form with weakly imposed boundary conditions, and prove strong stability [23,24]. The statistics of the solution such as the mean, variance and confidence intervals are computed non-intrusively using quadrature rules for the given stochastic distributions [10,12]. As an application, we analyze the heat transfer at rough surfaces in incompressible flow [2,18,22].
The paper will proceed as follows: in Sect. 2 we define the continuous problem in two space dimensions, transform it to the unit square using curvilinear coordinates and derive energy estimates that lead to well-posedness. We formulate a finite difference scheme for the continuous problem and prove stability in Sect. 3. In Sect. 4, we consider a heat transfer problem in incompressible flow. Finally, in Sect. 5 we draw conclusions.

The continuous problem
Consider the advection-diffusion problem on the stochastically varying domain (θ ) In (2.1),ū andv are the known mean velocities in the x-and y-directions satisfying the divergence relationū x +v y = 0 stemming from an incompressible Navier-Stokes solution. Furthermore, = (x, y, t) is a positive diffusion coefficient, u = u(x, y, t, θ ) represents the solution to the problem and θ = (θ 1 , θ 2 , . . . ) is a vector of random variables describing the geometry of the domain. F, g and f are data to the problem. The goal of this study is to investigate the effects of placing the boundary condition Hu = g at the stochastically varying boundary ∂ (θ ).

The transformation
We transform the stochastically varying domain into the unit square by the transformation, where 0 ≤ ξ, η ≤ 1. The Jacobian matrix of the transformation is given by, By applying the chain rule to (2.1) and multiplying by J = x ξ y η − x η y ξ > 0, we obtain The final formulation of the transformed problem is 1]. A more complete derivation of the transformed problem is included in "Appendix A". In (2.4), we have used the notation ∇ = ( ∂ ∂ x , ∂ ∂ y ) T . The transformed fixed domain including normal vectors are given in Fig. 1. Note that the wave speedsã andb depend on the stochastic variables θ .

The energy method
We multiply the transformed problem (2.3) with u, integrate over the domain (ignoring the forcing function F), and apply the Green-Gauss theorem. This yields and n is the outward pointing normal vector from ∂ , see Fig. 1.
The right-hand side (RHS) of (2.5) can be expanded as where for example the fluxes at the boundaries ξ = 1 and η = 1 arẽ respectively. Further, we note thatf andg can also be written in terms of u ξ and u η asf The formulation (2.6) in matrix form can be written The matrices in (2.7) are symmetric, and hence they can be diagonalized as forã,b = 0. By imposing the boundary conditions the (RHS) of (2.8) is bounded by data and hence gives an energy estimate.

Weak imposition of boundary conditions
As a preparation for the numerical approximation, we now impose the boundary conditions weakly using penalty terms. This gives dξ. (2.11) To illustrate the procedure, we assumeã ξ =1 ξ =0 > 0 andb η=1 η=0 > 0 and impose the boundary conditions using the operators in (2.10), which yields dξ.
(2.12) The indefinite terms in (2.12) are canceled by letting dξ, (2.14) which lead directly to an energy estimate. For generalã andb, the choices bounds the (RHS) of (2.11), in a similar way. The special cases withã,b = 0 are treated in a similar way, see "Appendix B".
We can now prove

Proposition 2.1 The problem (2.3) with the boundary conditions (2.9) and the penalty coefficients in (2.15) is strongly well-posed.
Proof Consider the specific case in (2.14). For other values ofã andb, the same general procedure is used. Time integration (from 0 to T ) of (2.14) results in In (2.16), the boundary terms with zero data all give a non-positive contribution, and hence the solution is bounded by data. The bound leads directly to uniqueness, and existence is guaranteed by the fact that we use the correct (i.e., minimal) number of boundary conditions.

The semi-discrete formulation
In this section we consider the numerical approximation of (2.3) formulated by using Summation-By-Parts (SBP) operators with Simultaneous Approximation Terms (SAT), the so called SBP-SAT technique [24]. First, we rewrite our variable coefficient continuous problem (2.3) using the splitting technique described in [13], to obtain, In (3.1), we note that the lower order terms vanish, sinceã ξ +b η = 0. The corresponding semi-discrete version of (3.1) including penalty terms for the boundary conditions isJ In (3.2) and (3.3), P −1 ξ,η Q ξ,η are the finite difference operators, P ξ,η are diagonal positive definite matrices, and Q ξ,η are almost skew-symmetric matrices satisfying The indices i = 0, 1, . . . , N and j = 0, 1, . . . , M correspond to the grid points in ξand η-direction.
To ease the notation we denote which corresponds to the continuous counterparts in (2.10). Finally, The penalty matrices Σ E , Σ W , Σ N and Σ S will be chosen such that the numerical scheme (3.2) becomes stable. For more details on the SBP-SAT techniques, see [24].

Stability
To prove stability (we only consider the west boundary, as the treatment of the other boundaries is similar), we multiply (3.2) with U T (P ξ ⊗ P η ) from the left, add the transpose of the outcome and define the discrete norm The remaining derivations leading to the discrete energy estimate is included in "Appendix C". We can now prove

Proposition 3.1 The numerical approximation (3.2) using the penalty coefficients
is strongly stable.
Proof For ease of presentation we prove the special case whenã,b > 0. By integrating (3.6) in time, considering also the remaining boundaries and using the penalty parameters in (3.7) we find As in the continuous energy estimate (2.14), the RHS of (3.8) consists of boundary data and negative semi-definite dissipative boundary terms which result in a strongly stable numerical approximation.
Remark 3.1 Note the similarity between the discrete energy estimate (3.8) and its continuous counterpart (2.16).

Remark 3.2
The possibility of applying the SBP-SAT technique to the coupled PDEs resulting from the use of polynomial chaos in combination with a stochastic Galerkin projection is shown in Pettersson et al. [19,20].

Rate of convergence for the deterministic case
We useū = sin(x) cos(y),v = − cos(x) sin(y), = 0.01 in order to satisfy the incompressibility condition. The rate of convergence is verified by computing the order of accuracy p defined as In (4.1), u h is the numerical solution, using the grid spacing h, and the manufactured solution is u a = sin(2π(x − t)) + sin(2π(y − t)).
The order of accuracy computed for different number of grid points and SBP-operators, is shown in Table 1. As time-integrator, the classical 4th-order Runge-Kutta method with 5000 grid points was used. The results shown in Table 1 confirm that the scheme is accurate for the 2nd-, 3rd-, 4th-and 5th-order SBP-SAT schemes [23].

Heat transfer at rough surfaces
Equipped with a provably stable scheme, we will now investigate the stochastic properties of a heat distribution problem in incompressible flow. The problem in two dimensions is of the form where we specify the following boundary conditions In (4.2), T is the temperature, (ū,v) the given velocity field, the viscosity. The boundary conditions in (4.3) are a well-posed subset of the general ones derived previously.
To simulate a boundary layer the quantitiesū,v and are chosen as where T ∞ = √ and ∂ T ∂n = n · ∇T . The velocity field is generated on the unit square (Fig. 2), then injected on the corresponding grid points on the varying domain. The simplified velocity field in (4.4) satisfies the divergence relationū x +v y = 0 and has a boundary layer.

Statistical results
In the calculations below, the 4th-order Runge-Kutta method is used together with 3rd-order SBP-operators on a grid with 50 and 100 grid points in the x-and y-direction and 9000 grid points in time, in order to minimize the time discretization error.
We start by enforcing the following stochastic variation on the south boundary of the geometry, (see Fig. 3) where θ 1 ∼ N (−1, 1) and θ 2 ∼ U (2, 10) are stochastic variables controlling the amplitude and frequency of the periodic variation respectively. In order to study the influence of different correlation lengths, we use y S long (x, θ 1 ) = 0.05θ 1 sin(2π x) y S short (x, θ 1 ) = 0.05θ 1 sin(2π x) + 0.005 sin(16π x), where we let y S short and y S long represent short and long correlation lengths respectively.
As typical measures of the results, we compute statistics of the integral of the solution and squared solution over the domain, that is u(x, y, t, θ 1 , θ 2 ) dx dy, and u 2 (x, y, t, θ 1 , θ 2 ) dx dy.
To compute the integrals in the stochastic analysis, we have used 20 grid points in both the θ 1 -and θ 2 -direction. For high-dimensional problems adaptive sparse grid techniques or multilevel Monte Carlo methods can be used to improve the efficiency of calculations like these, see for example [6,11]. However, in this particular case, with only two stochastic dimensions, straightforward quadrature is efficient enough. Figures 4 and 5 show the variance with respect to θ 2 of the integral of the solution and squared solution respectively, as a function of time for different realizations of θ 1 with y S as south boundary. Figures 4 and 5 both illustrate the fact that the variance increase with increasing amplitude, as could be expected. Figures 6 and 7 depict the variance with respect to θ 1 of the integral of the solution and squared solution respectively, as a function of time for fixed values of θ 2 using y S as south boundary. As can be seen, an increased frequency leads to an increased variance. Hence high-frequency random variation in the geometry affects the solution more than low-frequency random variation. Figures 8 and 9 illustrate the effects of correlation length on the variance. The variance as a function of time is shown for two different correlation lengths (one short y S short and one long y S long ). The figures show no significant difference between the two cases, and we conclude that the correlation length has a minor impact on the variance of the solution.

Conclusions and future work
We have studied how the solution to the advection-diffusion equation is affected by imposing boundary data on a stochastically varying geometry. The problem was transformed to the unit square resulting in a formulation with stochastically varying wave speeds. Strong well-posedness and strong stability were proven.  . 5 The variance of the non-linear functional u 2 (x, y, t, θ) dx dy with respect to θ 2 (frequency) for different realizations of θ 1 (amplitude) As an application, the two-dimensional heat transfer problem in incompressible flow with a given velocity field was studied. One of the boundaries was assumed to be stochastically varying. The geometry of the boundary was prescribed to have a periodic behaviour with stochastic variations in both amplitude and frequency. The variances were computed for different fixed realizations of θ 1 (when varying θ 2 ) and θ 2 (when varying θ 1 ) controlling the amplitude and frequency respectively. A tentative conclusion is that an increased frequency of the randomness in the geometry  Short correlation length Long correlation length Fig. 9 The variance of the non-linear functional u 2 (x, y, t, θ) dx dy with respect to θ 1 (amplitude) for different correlation lengths leads to an increased variance in the solution. Also, as expected, the variance of the solution grows as the amplitude of the randomness in the geometry increases. Finally, computational results suggests that the correlation length of the geometry has no significant impact on the variance of the solution.
In the next paper we will extend the analysis using polynomial chaos combined with stochastic Galerkin projection to incompletely parabolic systems, including calculations using the Navier-Stokes equations.
Due to symmetry, the matrices in (B.1) can be diagonalized as dξ.

(B.2)
To bound the RHS of (B.2) we impose the following boundary conditions As in the continuous case, see (2.15), we cancel the indefinite terms, by the choice The use of (C. and hence, the RHS of (3.6) is bounded by data with the initial assumption (as in the continuous case) thatÃ > 0. Note the resemblance between (C.4) and the related continuous estimate in (2.14) considering only the west boundary. When considering all the boundaries, the following choices give us a discrete energy estimate.