Uncertainty Quantification in Three Dimensional Natural Convection using Polynomial Chaos Expansion and Deep Neural Networks
UUncertainty Quantification in Three DimensionalNatural Convection using Polynomial Chaos Expansionand Deep Neural Networks
Shantanu Shahane , Narayana R. Aluru, Surya Pratap Vanka Department of Mechanical Science and EngineeringUniversity of Illinois at Urbana-ChampaignUrbana, Illinois 61801
Abstract
This paper analyzes the effects of input uncertainties on the outputs of athree dimensional natural convection problem in a differentially heated cubi-cal enclosure. Two different cases are considered for parameter uncertaintypropagation and global sensitivity analysis. In case A, stochastic variationis introduced in the two non-dimensional parameters (Rayleigh and Prandtlnumbers) with an assumption that the boundary temperature is uniform. Be-ing a two dimensional stochastic problem, the polynomial chaos expansion(PCE) method is used as a surrogate model. Case B deals with non-uniformstochasticity in the boundary temperature. Instead of the traditional Gaus-sian process model with the Karhunen-Lo` e ve expansion, a novel approachis successfully implemented to model uncertainty in the boundary condition.The boundary is divided into multiple domains and the temperature imposedon each domain is assumed to be an independent and identically distributed(i.i.d) random variable. Deep neural networks are trained with the boundary Corresponding Author
Preprint submitted to International Journal of Heat and Mass Transfer January 9, 2019 a r X i v : . [ c s . NA ] J a n emperatures as inputs and Nusselt number, internal temperature or veloc-ities as outputs. The number of domains which is essentially the stochasticdimension is 4, 8, 16 or 32. Rigorous training and testing process showsthat the neural network is able to approximate the outputs to a reasonableaccuracy. For a high stochastic dimension such as 32, it is computationallyexpensive to fit the PCE. This paper demonstrates a novel way of using thedeep neural network as a surrogate modeling method for uncertainty quan-tification with the number of simulations much fewer than that required forfitting the PCE, thus, saving the computational cost. Keywords:
Deep Neural Networks, Polynomial Chaos Expansion, NaturalConvection, Uncertainty Quantification
1. Introduction
Flow due to natural convection has been studied extensively in the lit-erature [1–7] since it is practically useful in cooling or heating systems forapplications like electronics, nuclear reactors, computing servers etc [8–10].De Vahl Davis [1], Shu et al. [6] and Le Qu´er´e [3] simulated natural convec-tion in a two dimensional differentially heated square cavity with gravity ina direction orthogonal to the applied temperature difference. De Vahl Davis[1] used the stream function-vorticity formulation for laminar flow. Accuratebenchmark solution was obtained using mesh refinement and extrapolation.Le Qu´er´e [3] used a pseudo-spectral algorithm combining spatial expressionsof Chebyshev polynomial series with a finite difference time marching scheme.The results for Rayleigh number upto 10 are presented. Shu et al. [6] solvedthe same problem with local radial basis function based differential quadra-2ure (RDF-DQ) method. This method is a mesh-free approach with the RBFsas test functions to estimate the derivatives at any node as a weighted sumof values at the neighboring nodes. They discussed the effects of the RBFshape parameter and its fine tuning to get accurate solutions. Fusegi et al.[2] presented the results for three dimensional differentially heated cubicalenclosure. They used the finite difference discretization with SIMPLE algo-rithm [11] for laminar flow at Rayleigh numbers in the range of 10 to 10 .Rayleigh-B´enard is another class of natural convection problems in whichthe temperature difference is applied parallel to the direction of gravity withthe lower wall heated and upper wall cooled. Hu et al. [4], Li et al. [5] andYigit et al. [7] studied the Rayleigh-B´enard convection numerically in twoand three dimensional cavities with cubical and cylindrical shapes havingvarious aspect ratios.The numerical and experimental results of natural convection are utilizedfor verification and validation of the numerical software packages. Due toinaccuracies in the measurement and control, the experimental results areprone to errors. These errors can be estimated by introducing stochasticvariations in the inputs and propagating these to the outputs. It is possiblethat the stochastic mean of an output parameter is different from its value atthe input mean. Deterministic simulations alone cannot estimate the shift ofmean. Thus in the recent years, there has been a growing interest in analysisof the effects of stochastic variations in the inputs on the outputs. There aremultiple examples in the literature in which the uncertainty propagation tech-niques are combined with the deterministic numerical simulations [12–23]. Itis popular to use the polynomial chaos expansion (PCE) for uncertainty prop-3gation in which the output is approximated as a summation of polynomialbasis which are functions of the stochastic inputs. The two main classes ofmethods to estimate the coefficients of the PCE are stochastic Galerkin pro-jection [12–15] and collocation [17–19]. Stochastic Galerkin method requiresmodification of the underlying deterministic code since it requires solution ofa new set of equations and thus, is called as intrusive method. This becomesa significant additional effort of software development and it is difficult tocouple with the legacy codes. Hence, recently non-intrusive stochastic collo-cation methods have gained popularity. The basic idea is to have multipleevaluations of the deterministic simulation at predefined collocation pointswhich are samples from the underlying probability distribution function ofthe input parameters. The PCE coefficients are then estimated from theoutput values obtained from the deterministic solution at these input sam-ples. The coefficients can be used for post-processing operations like outputstatistics estimation, response surface plotting and sensitivity analysis.Uncertainty quantification (UQ) for various types of natural convectionproblems has been studied in the literature [15–17, 20]. Maitre et al. [15] usedthe zero-Mach-number model to simulate natural convection in a two dimen-sional differentially heated square cavity with uncertainty in the cold walltemperature. The random component of the cold wall temperature is mod-eled using the Gaussian process with an auto-correlation function which isapproximated by the truncated Karhunen-Lo` e ve (KL) expansion. PCE coef-ficients are estimated by the stochastic Galerkin projection. Output statisticsfor various values of the non-Boussinesq parameter (cid:15) are presented. Ganap-athysubramanian and Zabaras [17] presented an adaptive refinement based4pproach for reducing the number of deterministic simulations in high di-mensional stochastic simulations. The adaptive sampling method is appliedto the two dimensional natural convection problem with random boundarycondition which is modeled by the KL expansion method. Venturi et al. [16]studied the stability of the two dimensional Rayleigh-B´enard convection sub-ject to stochastic boundary temperatures. The random boundary conditionis assumed to be a non-uniform Gaussian random processes approximatedby the KL expansion. It is found that the stochastic wall temperatures canextend the stability range of quasi-conduction states beyond the classicalbifurcation point. Fajraoui et al. [20] analyzed the natural convection ofporous media in a two dimensional differentially heated square cavity withuncertainty in the Rayleigh number, permeability anisotropy ratio, disper-sion coefficients and heterogeneity variation. PCE method is used to estimatethe statistics and sensitivity of the output parameters such as temperatureand Nusselt number distributions.The PCE method is extremely useful for low dimensional uncertaintyquantification. But at higher dimensions, it faces the problem known as‘curse of dimensionality’ i.e., for a linear increase in the stochastic dimen-sions, the number of samples grows exponentially. The Smolyak algorithm[24] addresses this problem to some extent by reducing the number of samplesin high dimensions without compromising the interpolation accuracy. Evenwith the use of the Smolyak algorithm, number of samples of the order of10 − are required for five or more input dimensions. For instance, aneight and sixteen dimensional problem needs 3905 and 51073 samples respec-tively, for the accuracy level of five [25]. Practically, it is computationally5xpensive to simulate the deterministic solution thousands of times. Thus, analternate method is required for uncertainty propagation. The Monte Carlomethod is a simple approach which approximates the statistics of the outputby running the deterministic simulations at pseudo random samples of theinputs [26]. Since the error using the Monte Carlo method is O ( / √ n ), thenumber of samples is practically too high which makes using the Monte Carlomethod directly with the deterministic simulation difficult. Thus, it is popu-lar to use a surrogate model which is trained and tested using deterministicsimulations. A good surrogate model can be trained with a small numberof deterministic simulations and its evaluation is cheap. A well tested surro-gate model is further used to estimate the outputs at multiple sample inputs.Since the surrogate model evaluation is cheap, there is practically no limiton the number of input samples for the Monte Carlo method. Note that thePCE is also a surrogate model which is ideal for low stochastic dimensionswith a possibility of direct estimation of the output statistics without the useof the Monte Carlo method.In order to simulate a high dimensional stochastic problem, a neural net-work (NN) is used as a surrogate model. Hornik et al. [27] showed thatmultilayer feed forward networks are universal approximators i.e., with mildassumptions on the underlying function to be approximated, the network canachieve any desired degree of accuracy by choosing suitable number of neu-rons. The NNs can handle the ‘curse of dimensionality’ by multiple nonlinearactivation functions. In recent years, the NNs have been extremely popularin many fields of work as discussed by the review paper by Schmidhuber [28].Here, only the applications related to surrogate modeling for numerical sim-6lations are discussed [29–36]. Sablani et al. [30] used a NN as a surrogatemodel for inverse heat conduction problem of estimation of heat transfer co-efficient from the temperature-time history at different locations. Since theNN is trained for the inverse problem directly using forward deterministicsimulations, the estimation can be done non-iteratively. Cz´el et al. [32] simi-larly used NN for non-iterative estimation of heat capacity and temperaturedependent thermal conductivity using the experimental transient tempera-ture histories. A radial basis function type NN is trained using the numericalsolution of the direct heat conduction problem. Both the above publicationsshow that computational time is saved by non-iterative estimation due tothe NN surrogate model coupled with the forward numerical simulations.Gholami et al. [29] trained a NN for a three dimensional two fluid flow in a90 o curved channel and compared both the numerical simulations and NNpredictions with experimental data. It is reported that the NN model isreasonably accurate and significantly faster compared to the full numericalsimulation. Tripathy and Bilionis [34] trained a NN to solve a steady statetwo dimensional diffusion process with spatially varying uncertainty in thediffusion coefficient. This uncertainty is modeled as a log normal randomfield with mean and covariance functions of the Gaussian random field whichis approximated by the Karhunen-Lo` e ve (KL) expansion. Using the trainedNN as a surrogate with diffusion coefficient as the input, statistics of theoutput parameter are estimated. Zhang et al. [36] have shown the utilityof the physics informed neural networks (PINNs) for uncertainty quantifica-tion in direct and inverse stochastic problems. The basic idea of a PINN isto minimize the residual when the NN is substituted in the model differen-7ial equation together with the standard loss function of the NN. Automaticdifferentiation is used to estimate the residual. The paper claims that min-imizing the residual along with loss function enhances the accuracy of theprediction.This paper presents results of input uncertainty propagation for a threedimensional natural convection problem in a differentially heated cubical en-closure. Two different cases of input uncertainties are considered. Case Aassumes that the boundary conditions are uniform. Thus, uncertainty isintroduced in the two non-dimensional parameters (Rayleigh and Prandtlnumbers). For this case, the polynomial chaos expansion (PCE) method isused as a surrogate model with stochastic collocation to estimate the PCEcoefficients. Since this is a two dimensional stochastic problem, the numberof samples required is small enough and thus, the estimation of the PCEcoefficients is feasible. Case B deals with non-uniform stochastic boundarycondition with deterministic material properties of the fluid. Since the tem-perature difference between the opposite walls drives the natural convectionflow, the cold wall is held at a constant temperature and uncertainty is in-troduced in the hot wall temperature. Although the conventional methodto deal with boundary condition uncertainties is to use the Gaussian processmodel with the Karhunen-Lo` e ve expansion [15–17, 34], in this work, a novelapproach is successfully implemented. The basic idea is to divide the hotwall into multiple domains and impose a temperature boundary conditionon each domain. It is assumed that each domain temperature is an inde-pendent and identically distributed (i.i.d) random variable. Sets of forwarddeterministic simulations are used to train deep neural networks (DNNs) with8he boundary temperatures as inputs and Nusselt number, internal temper-ature or velocities as outputs. The number of domains which is equal to thestochastic dimension ranges from 4 to 32. The DNN is successfully trainedand tested with less number of samples compared to those required for PCEcoefficient estimation for the high dimensional problem. This new approachto deal with stochastic boundary conditions with DNN as a surrogate modelis found to be much better than those methods presented in the literature sofar.
2. Deterministic Problem Description
In this paper, we consider the three dimensional natural convection in adifferentially heated cube. A temperature gradient is applied on two oppositefaces ( X = 0 and X = 1) of a cube with sides of length L . The remainingfour faces are thermally insulated. Gravity is imposed in the Y directionwhich is orthogonal to the direction of the temperature gradient. Because ofthe thermal expansion of the fluid inside the cube due to the temperaturevariation, a buoyant force causes the lighter fluid to move upwards thuscreating currents.The flow field can be described by three dimensional incompressible Navier-Stokes and energy equations. For moderate density variations, natural con-vection can be modeled using the Boussinesq approximation. The system ofequations is written in terms of non-dimensional variables as follows [3]: ∇ · u = 0 (1)9 u ∂t + ( u · ∇ ) u = P rRa . ∇ u − ∇ P − ˆ g P r
Θ (2) ∂ Θ ∂t + ∇ · u Θ = 1 Ra . ∇ Θ (3)where, u is the velocity vector, Θ is the temperature, t is time, ˆ g is theunit vector in the direction of gravity, P is the pressure, P r = ν / α is thePrandtl number and Ra = gβ ∆ T L / να is the Rayleigh number. Characteristicvalues for non-dimensionalization are as follows: velocity u c = ( α / L ) Ra . ,time t c = ( L / α ) Ra − . , pressure P c = ρu c and L is the cavity length. Non-dimensional temperature is defined as Θ = ( T − T m ) / ( T h − T c ) where, T h and T c are hot and cold wall temperatures respectively and T m = ( T h + T c ) / is themean temperature. The governing equations (1-3) are solved using the software OpenCast[37] with finite volume method on a collocated grid. The fractional stepmethod [38] is used to integrate the equations. An intermediate velocityfield ( u* ) is first estimated by solving the modified momentum equation (4)without the pressure gradient. The diffusion term is discretized implicitlyusing second order Crank-Nicolson method whereas an explicit second orderAdams-Bashforth is used for the convection term. u* − u n ∆ t = − Conv ( u n , u n − ) + Dif f ( u* , u n ) + Buoy ( T n ) (4)10he velocity correction equation (5) is obtained by subtracting the modifiedmomentum equation from the original equation. u n +1 = u* − ( ∇ Φ) n +1 ∆ tρ (5)Imposing divergence free condition on the ( n + 1) velocity field gives thePoisson equation for Φ (6). ∇ · (cid:18) ∇ Φ ρ (cid:19) n +1 = ∇ · u* ∆ t (6)The overall algorithm for marching from time step n to n + 1 can be sum-marized as follows:1. Solve for u* using eq. (4)2. Solve the Poisson equation for Φ (6) iteratively to estimate Φ n +1
3. Correct the velocities ( u n +1 ) using eq. (5)4. Solve the energy equation (3)Note that the pressure P can be estimated from Φ if required: P = Φ − µ ∆ tρ ∇ Φ (7)11 . Grid Independence Study and Verification (a) Temperature vs X: Ra = 10 (b) Temperature vs X: Ra = 10 (c) Y vs X Velocity: Ra = 10 (d) Y vs X Velocity: Ra = 10 e) Y Velocity vs X: Ra = 10 (f) Y Velocity vs X: Ra = 10 Figure 1: Temperature and Velocity Contours
The numerical results from the software are verified using the three di-mensional natural convection simulations of Fusegi et al. [2]. For Rayleighnumber up to 10 , the flow remains laminar [2]. Hence Rayleigh numbers of10 and 10 are used for verification. The cube is meshed with three differ-ent grid sizes with 32 , 64 and 128 structured hexahedrons. Steady statesolution is estimated by time marching. For all the scalar fields φ , when thenon-dimensional steady state error computed over the entire domain definedas max ( | φ new − φ old | ) / max ( | φ new | ) is less than 10 − , it is assumed that the steadystate is reached. The temperatures at boundary faces X = 0 and X = 1 areset to 0.95 and 1.05 respectively.The temperature and velocities are plotted along a centerline for boththe Rayleigh numbers (fig. 1). Note that the characteristic velocity used byFusegi et al. [2] is different compared to the discussion in section 2.1. Hence13nly for verification, the velocities are scaled by u c = (cid:112) gβL ( T h − T c ). Ineach figure, estimates from the three grid levels computed by OpenCast [37]are superimposed with the results from Fusegi et al. [2] whenever available.The numerical estimates from OpenCast match well with the published re-sults thus, verifying the code. For the Ra = 10 case (figs. 1a, 1c and 1e),it can be seen that all the three grid results from OpenCast overlap witheach other. For Ra = 10 (figs. 1b, 1d and 1f), the coarsest grid (32 ) plot isslightly off but the remaining two finer grid plots overlap. This shows thatgrid independence is achieved and for all further computations, a grid of size101 is used.
4. Parameter Uncertainty Quantification
Engineering problems have uncertainties in input parameters due to diffi-culty in precise measurement and control. It is necessary to propagate the in-put stochasticity to the output. Let w ( x , ξ ) be a function which maps inputsto an output. Here, the vector x denotes all the deterministic inputs whereas,the vector ξ denotes all the stochastic parameters. It is assumed that thestochastic variable follows a known probability distribution: ξ ∼ f ( ξ ). Theaim is to estimate the stochastic mean of the output w ( x , ξ ) defined as: w f ( x ) = ˆ w ( x , ξ ) f ( ξ ) d ξ (8)Since the function w cannot be expressed in a closed form, above integral hasto be approximated numerically. Importance sampling based Monte Carlo14pproximation of the integral is given by [26]: w f ( x ) ≈ n i = n (cid:88) i =1 w ( x , ξ i ) f ( ξ i ) p ( ξ i ) (9)where, ξ i are n samples drawn from the probability distribution p ( ξ ). It iseffective to set p ( ξ ) as f ( ξ ) to reduce variance i.e., ξ i ∼ f ( ξ ) [26]. Thus,eq. (9) is simplified to: w f ( x ) ≈ n i = n (cid:88) i =1 w ( x , ξ i ) (10)The error in the integral estimated using the Monte Carlo method is O ( / √ n ). Thus, the number of samples ( n in eq. (10)) can be of the order ofthousands or more depending on the integrand. It is practically not possibleto have so many samples as each sample corresponds to one deterministicsimulation which is typically computationally costly. Thus, a surrogate modelwhich is trained and tested using deterministic simulations is further used toestimate the outputs at multiple sample inputs. Since the surrogate modelevaluation is cheap, there is practically no limit on the number of inputsamples.The input parameters affecting the simulation of natural convection areboundary temperature, domain length and material properties like viscos-ity, thermal diffusivity and coefficient of thermal expansion. Two separatecases of input uncertainties are analyzed in this paper. Case A assumesthat the boundary temperature is uniform. Thus, the governing equations(1-3) show that the physics can be parametrized just using Rayleigh and15randtl numbers. The uniform boundary temperature, domain length andthe material properties all are included in the Rayleigh and Prandtl numbers.Hence, from the perspective of parameter uncertainty propagation, this is atwo dimensional problem. Case B considers the possibility of uncertainty inthe non-uniform boundary temperature with deterministic material proper-ties. Case A is analyzed using the polynomial chaos expansion (Section 4.1)whereas, deep neural networks (Section 4.2) are used as surrogate models forthe case B.Such an analysis is practically important as there are stochastic variationsin the boundary conditions due to inaccuracy in measurement and control.The fluid material properties also vary stochastically due to the presenceof impurities. Hence, the effect of these uncertainties on the temperatureand velocity distribution and Nusselt number is studied in this work. Thefollowing sections summarize the surrogate modeling strategies. Literature on uncertainty quantification describes various methods to es-timate the relationship between stochastic inputs and outputs. Most of thesemethods rely on the idea of expanding the outputs as a linear combinationof polynomial basis functions in the stochastic dimension. Orthogonal poly-nomial is a popular choice as basis since orthogonality helps in convergence.Xiu and Karniadakis [39] showed that the Wiener’s polynomial chaos [40]with orthogonal polynomials of Askey family lead to optimal convergence ofthe interpolation error. They identified which orthogonal polynomial basis issuitable depending on the probability distribution function followed by thestochastic variable. For example the Hermite polynomials are orthogonal to16ach other when the standard normal distribution is used as weighting func-tion. Thus, Hermite polynomial basis function is recommended when thestochastic variable follows normal distribution. A polynomial chaos series(eq. (11)) is used to expand a second order random field [39]. The series istruncated to order n for all practical purposes. w ( x , ξ ( θ )) = ∞ (cid:88) i =0 w i ( x ) Ψ i ( ξ ( θ )) ≈ n (cid:88) i =0 w i ( x ) Ψ i ( ξ ( θ )) (11)where, x is the spatial variable, ξ = ( ξ , ξ , . . . ξ r ) is the random variablevector, θ is an elementary event, Ψ i is a multi-dimensional orthogonal poly-nomial of order i , w is the output to be estimated and w i are the coefficientsof the series. In this case, ξ is a two dimensional vector following normaldistribution and Ψ i is a two dimensional Hermite polynomial.In this paper, stochastic collocation method is used to estimate the de-terministic coefficients ( w i ) of the polynomial chaos expansion. Collocationacts as a wrapper over the existing deterministic software since it is a non-intrusive method. Thus, no modification of the deterministic software isrequired. After deterministic simulations at M sample points ( ξ m ), a con-straint w ( x , ξ m ) = w sim ( x , ξ m ) is imposed. The left hand side is estimatedfrom polynomial chaos expansion (eq. (11)) and right hand side from eachdeterministic simulation. These M constraints can be written in the matrixvector form [41]. For accuracy, it is recommended to have more samples thanthe number of basis functions ( M > n +1) and thus the Vandermonde system17eq. (12)) is overdetermined. Ψ ( ξ ) · · · Ψ n ( ξ )... ... Ψ (cid:0) ξ M (cid:1) · · · Ψ n (cid:0) ξ M (cid:1) w ( x )... w n ( x ) = w sim ( x , ξ )... w sim (cid:0) x , ξ M (cid:1) (12)Sampling strategy plays a vital role in the accuracy and stability ofstochastic collocation. Uniformly distributed samples lead to highly oscil-latory interpolation and hence, poor convergence. Thus, for one dimensionalstochastic problems, the roots of the basis orthogonal polynomials is a popu-lar choice of sample points ( ξ m ) [41]. For higher stochastic dimensions, tensorproduct of the single dimensional samples can be used. The sample size growsexponentially with dimensions if tensor product is used. This is a problemas each sample corresponds to a deterministic simulation and thus, the costof computations grows exponentially. An algorithm to reduce the number ofsamples in high dimensions without spoiling the interpolation accuracy waspublished by Smolyak [24]. It is found that for stochastic dimensions lessthan three, the Smolyak algorithm is not effective in reducing the samplesize [25]. Thus, in this work, a tensor product of single dimensional samples(roots of the Hermite polynomials) are used [25]. UQLab, a MATLAB basedtool developed by Marelli and Sudret [42] is used for estimation of polyno-mial chaos coefficients and response surfaces. The Polynomial Chaos-Krigingmodule of UQLab is used as it is found to be more effective than the basic18olynomial chaos method. Practically, a heat exchanger setup with a closed loop feedback systemis used to maintain the cold and hot wall temperatures. Due to the er-rors in measurement and control, there are stochastic variations in the settemperature. For a large wall, there would be multiple heat exchangers incontact with the wall. Since the objective is to maintain a uniform temper-ature, the design and operation of all the heat exchangers would be similar.Thus, it is safe to assume that the temperature achieved by each of themis a random variable following the normal distribution with mean as the ex-pected temperature and error modeled as the standard deviation. It is alsoassumed that the heat exchangers are independently controlled. Thus, theset temperatures are independent and identically distributed (i.i.d) randomvariables.In the literature, there are examples of using Gaussian process modelwith the Karhunen-Lo` e ve expansion to model uncertainty in the boundarycondition [15–17]. In the present work, the wall is subdivided into multipledomains with different values of temperatures imposed as boundary condi-tions. From fig. 9, it can be seen that the variation of Nusselt number isstronger in the direction of gravity ( Y ) compared to the orthogonal direction( Z ). Thus, the wall is divided into strips along gravity ( Y ). Figure 2 showssamples of the temperature boundary condition with 4, 8 and 16 number ofstrips. Each strip temperature is assumed to be an i.i.d. random variablefollowing a normal distribution ( µ = 1 . , σ = 0 . a) Number of Strips: 4 (b) Number of Strips: 8 (c) Number of Strips: 16 Figure 2: Samples of Temperature Boundary Condition
A neural network is a set of interconnected nodes such that the informa-tion flows from inputs to outputs. Each node is known as a neuron. Figure 3ashows a single neuron which has n scalar inputs ( x , x , ..., x n ) and single out-put ( y ). Each neuron performs the following two operations in sequence:1. Linear transformation: a = (cid:80) ni =1 w i x i + b ; where, w i are the weightsand b is a bias term2. Element-wise nonlinear transformation: y = σ ( a ); where, σ is the acti-vation functionA neural network is formed by stacking single neurons in a layer and con-necting multiple layers as shown in fig. 3b. It depicts an input (layer L1),an output (layer L4) and two hidden layers (layers L2 and L3). The arrowsindicate the direction of information flow from input to output layer throughthe hidden layers. A deep neural network (DNN) is essentially a neural net-work with multiple hidden layers. Adding multiple hidden layers increasesthe nonlinearity of the network and thus, the network can approximate morecomplex functions successfully. The number of neurons in the input and out-put layers is specified by the problem definition whereas, number of hidden20ayers and neurons has to be fine tuned. (a) Single Neuron (b) Deep Neural Network Figure 3: Neural Network Schematics
The linear transformation followed by the nonlinear activation functionof each neuron can be written in a single matrix vector equation: y ( j ) = x if j = 1 σ ( W ( j ) y ( j − + b ( j ) ) ∀ j ∈ { , , ..., L } (13)where, x ∈ R l is the input vector, W ( j ) ∈ R l j × l j − is the matrix of theweights, y ( j ) ∈ R l j is the activation produced by the j th layer and b ( j ) ∈ R l j is the bias. L is the total number of layers including input, output and hiddenlayers. Number of neurons in the j th layer is denoted by l j . For instance, infig. 3b, L = 4, l = 4, l = 8 and l = l = 6. Applying eq. (13) sequentially21tarting from the input layer is known as forward propagation. This opera-tion estimates the output vector ( y ( L ) ) from the input vector ( y (1) = x ) ifthe weights and bias are known. Logistic sigmoid, hyperbolic tangent andrectified linear unit (ReLU) are some of the popular activation functions [43].In this work, the ReLU function defined by σ ( y ) = max { , y } is used for allthe hidden layers. For output layer, in order to allow negative values, theidentity function σ ( y ) = y is used.The process of estimation of weights and bias using a given set of inputs( x i , ≤ i ≤ m ) and the corresponding outputs ( z i , ≤ i ≤ m ) is knownas training. For the given set of m training samples, mean squared errorbetween the neural network estimate ( ˆ z i ) and the true value ( z i ) of theoutput is defined as the loss function: L ( W , b ; x i , z i ) = 1 m i = m (cid:88) i =1 || z i − ˆ z i || (14)It is commonly seen that the neural network performs well on the trainingdata but performs poorly on the unseen test data. This phenomenon isknown as overfitting and is controlled with regularization. Goodfellow et al.[43] discuss various regularization methods in detail. Here, the L2 weightregularization with parameter λ is used in which, the loss function (eq. (14))is modified: L ( W , b ; x i , z i ) = 1 m i = m (cid:88) i =1 || z i − ˆ z i || + λ L (cid:88) l =1 || W ( l ) || (15)The weights and bias which minimize the loss function are estimated using22 numerical optimization algorithm. The gradient of the loss function withrespect to the weights and bias is required in the optimization algorithmslike gradient descent. The gradient is estimated by the backpropagationalgorithm [44]. The Adam optimizer is used in this work with the parameters β and β as suggested by Kingma and Ba [45]. Other hyper-parameters likelearning rate, regularization constant, number of hidden layers and numberof hidden units are tuned using a validation set. All the implementationdetails are given in section 5.3.1.
5. Uncertainty Propagation Results
Uncertainty quantification analyzes the effects of small stochasticity inthe input on the output. Since the stochasticity in the input is of the orderof a small percentage of its mean, a similar variation is expected in the out-put. Thus, in this paper, a comparison of the stochastic mean of each outputvariable is done with the value of that variable at the input mean. The tem-perature and velocities are non-dimensionalized as discussed in section 2.1.Figures 4 and 5 plot the Nusselt number at hot wall and temperature, Xand Y velocities along the Z midplane for Rayleigh number of 10 and 10 ,respectively. The stochastic means of the Nusselt number, temperature andvelocities are expected to follow trends similar to the figs. 4 and 5. Thus,in the following sections, contour plots of the difference between the outputstochastic mean and the deterministic values are plotted for comparison. Thedifference gives an estimate of the effect of uncertainty in the input on theoutputs compared to the case with deterministic inputs.23 a) Hot Wall Nusselt Number (b) Temperature at Z = 0 . Z = 0 . Z = 0 . Figure 4: Deterministic Results for Ra = 10 and P r = 7 . a) Hot Wall Nusselt Number (b) Temperature at Z = 0 . Z = 0 . Z = 0 . Figure 5: Deterministic Results for Ra = 10 and P r = 7 . The non-dimensionalized governing equations (1-3) show that the naturalconvection problem is parametrized by two parameters viz. the Rayleighnumber and the Prandtl number. It is assumed that both of them follow anormal distribution with a 2% standard deviation with respect to mean: • Ra ∼ N ( µ = 10 , σ = 0 . µ ) or Ra ∼ N ( µ = 10 , σ = 0 . µ )25 P r ∼ N ( µ = 7 . , σ = 0 . µ )Results from two different Rayleigh numbers (10 and 10 ) for which the flowis known to remain laminar are presented here [2]. High Rayleigh number im-plies higher buoyancy compared to viscous forces and thus, higher velocitiesand Nusselt number are observed. The fluid inside the cube is assumed to bewater and hence, the Prandtl number is set to 7.5. As mentioned in section4.1, a tensor product of roots of Hermite polynomial scaled with mean andstandard deviation is chosen as samples. Hundred uniform Latin hypercube samples are used as test points toverify the convergence of the stochastic collocation method. Two indepen-dent estimates of the same output parameter are obtained using polynomialchaos expansion and deterministic simulation. The root mean square of thedifference between these two estimates normalized by the maximum valueof the parameter is defined as the non-dimensional error estimate. Spatialmean Nusselt number (eq. (16)) over the hot face is used to estimate thecollocation error.
N u mean = ˆ N u ( y, z ) (cid:12)(cid:12)(cid:12) x =1 dydz = ˆ ∂T ( y, z ) ∂x (cid:12)(cid:12)(cid:12) x =1 dydz (16)First column of table 1 is the accuracy level of the sample points used forinterpolation. Accuracy level l integrates polynomials up to degree 2 l − M in eq. (12)). The last two columns listthe non-dimensional RMS error in computation of the spatial mean Nusselt26ccuracy Level Table 1: Stochastic Collocation Error Analysis number for both the Rayleigh numbers. Although the error increases slightlyat the first two levels, eventually the error drops with higher accuracy level.This proves the convergence of the stochastic method. At the highest accu-racy level, the error is of order 10 − or 10 − which shows that the polynomialchaos is reasonably accurate and can be used for further analysis. (a) Ra = 10 (b) Ra = 10 Figure 6: Mean Nusselt Number Estimate from Numerical Simulation and PolynomialChaos Expansion
For visual inspection, spatial mean Nusselt number estimates from thenumerical simulation and the polynomial chaos expansion are plotted to-27ether in fig. 6 for the hundred test points. Ideally, all the points should lieon the Y = X line but due to the stochastic interpolation error, some pointsare off the line. Since most of the points follow the expected trend of the Y = X line, it can be concluded that the polynomial chaos is accurate. Response surface gives a visual representation of the variation of an out-put parameter with input stochastic parameters. For a two dimensionalstochastic problem, the response surface can be plotted as a contour. Fig-ure 7 plots the response surfaces of Nusselt number averaged over the hot wall(eq. (16)) for both the Rayleigh numbers. In each plot, X and Y axes denoteRayleigh and Prandtl numbers, respectively which are the stochastic inputparameters. Since both the input variables are assumed to follow normal dis-tribution, they are plotted in the range ( µ − σ, µ +3 σ ). For example, Rayleighand Prandtl numbers are plotted in the range (10 − × , + 3 × . − × . , . × . a) Ra = 10 (b) Ra = 10 Figure 7: Spatial Mean Nusselt Number Response Surface (a) Diff. of StochasticMean and DeterministicValue (b) Stochastic Std. Dev. (c) Stochastic Std. Dev.by Stochastic Mean
Figure 8: Local Nusselt Number at Hot Wall for Ra = 10 The local Nusselt number on the wall (
N u ( y, z )) varies due to the stochas-ticity in the Rayleigh and Prandtl numbers. Figures 8 and 9 plot its statisticsfor both the Rayleigh numbers. The stochastic mean plots look visually sim-ilar to the plots at mean Rayleigh and Prandtl numbers (figs. 4a and 5a).Thus, figs. 8a and 9a plot contours of difference between the stochastic meanand deterministic value of the Nusselt number at hot wall (figs. 4a and 5a).29 a) Diff. of StochasticMean and DeterministicValue (b) Stochastic Std. Dev. (c) Stochastic Std. Dev.by Stochastic Mean Figure 9: Local Nusselt Number at Hot Wall for Ra = 10 Figures 8b and 9b plot the standard deviation due to the stochastic inputvariation. Higher standard deviation is observed in the region of higher mean.Thus, figs. 8c and 9c are plotted to annihilate the effect of mean. The localratio of standard deviation to mean shows higher values on the left and rightsides near Y = 0 and Y = 1 as gravity is acting in Y direction. It can beseen that the difference between the stochastic mean and the deterministicvalues (fig. 8a) is one order of magnitude higher than the stochastic standarddeviation (fig. 8b) for Rayleigh number of 10 . This implies that the inputstochasticity shifts the deterministic mean of the output more than its stan-dard deviation. On the other hand, for Rayleigh number of 10 , both thedifference (fig. 9a) and stochastic standard deviation (fig. 9b) are of similarorders of magnitude. Thus, the shifting of mean and the standard deviationon the shift are of similar orders of magnitude for higher Rayleigh number.The difference between the stochastic mean and the deterministic value is ofthe order of 0 .
3% of the deterministic value for both the Rayleigh numbers.30 .2.3. Velocity and Temperature (a) Temperature (b) X Velocity (c) Y Velocity
Figure 10: Difference between Stochastic Mean and Deterministic Value at Z = 0 . Ra = 10 (a) Temperature (b) X Velocity (c) Y Velocity Figure 11: Stochastic Standard Deviation at Z = 0 . Ra = 10 For the case of temperature gradient in the X direction and gravity inthe Y direction, temperature and velocity contours at the Z mid-plane arequite informative. All the stochastic mean contour plots look visually similarto those of the deterministic natural convection problem (figs. 4 and 5) andthus, are not plotted. Instead, the difference between the stochastic mean anddeterministic values is plotted here. Figures 10 to 13 plot the difference andstandard deviation of temperature and velocities for the Rayleigh numbersof 10 and 10 , respectively. It is observed that both the difference and31tandard deviations are of similar orders of magnitude. The difference andstandard deviation in the temperature are three orders of magnitude smallerthan the mean. On the other hand, the difference and standard deviation invelocities are two orders of magnitude smaller than the mean. Thus, it can beconcluded that the effect of uncertainty is significant on the velocities thanon the temperature. The uncertainty has higher impact on the temperatureat Ra = 10 simulation than at Ra = 10 . Similar to the Nusselt number, thestandard deviation is higher when the mean value is higher. The differencebetween the stochastic mean and the deterministic value for temperature isof the order of 0 . − .
3% of the deterministic value for low and high Rayleighnumbers respectively. The difference is of the order of 1 −
3% for both thevelocities of low and high Rayleigh numbers respectively. Thus, the effect ofinput stochasticity is higher on the velocities compared to temperature. (a) Temperature (b) X Velocity (c) Y Velocity
Figure 12: Difference between Stochastic Mean and Deterministic Value at Z = 0 . Ra = 10 a) Temperature (b) X Velocity (c) Y Velocity Figure 13: Stochastic Standard Deviation at Z = 0 . Ra = 10 (a) Ra = 10 (b) Ra = 10 Figure 14: Sensitivity: 6 Outputs, 2 Inputs
Sensitivity analysis quantifies the variation in an output due to the vari-ation in a particular input. In this work, the global sensitivity is estimatedusing the Sobol indices based on the Sobol decomposition [46]. Partial Sobolindex measures the contribution of a subset of inputs to the total variance.This includes the variance due to a coupling between the inputs. Total Sobolindex for each input is defined as the sum of all the partial indices involv-ing that input parameter. Here, the total Sobol indices are estimated fromthe polynomial chaos coefficients using the sensitivity analysis tool of the33oftware UQLab [42].Wall Nusselt number and X and Y velocities (u,v) have been chosenas representative outputs of the natural convection problem to study thesensitivity. The mean and maximum are taken over the hot wall and entirecube for the Nusselt number and velocities, respectively. Figure 14 plotsthe sensitivity (total Sobol index) of the mean and maximum of Nusseltnumber and velocities with respect to each stochastic input parameter viz.Rayleigh and Prandtl numbers. It can be seen that the Nusselt number ismore sensitive to Prandtl number for the case of Ra = 10 whereas, it ismore sensitive to Rayleigh number for Ra = 10 . On the other hand, thevelocities are more sensitive to the Prandtl number in both the cases. Butfor the higher Rayleigh number case, their sensitivity towards the Prandtlnumber increases. As described briefly in Section 4.2, the hot wall is divided into strips ina direction orthogonal to gravity (fig. 2). Different boundary temperaturesare prescribed on each strip. It is assumed that each strip temperature is anindependent and identically distributed (i.i.d) random variable and follows anormal distribution with µ = 1 .
05 and 3 σ = 0 .
01. Since the cold wall is heldat a constant temperature of 0.95, the mean temperature difference drivingthe natural convection flow is still 0.1. This implies that a 3 σ error of 10%is specified in the input stochasticity. The number of strips in this study isvaried from 4 to 32 in multiples of 2. The material properties of the fluidare kept constant. In order to estimate the statistics of the outputs, a deep34eural network (DNN) surrogate model is used. The boundary temperatureson each strip are inputs to the DNN. A 99 finite volume mesh is used forthe numerical simulation. Four separate DNNs are trained with the followingoutputs:1. Nusselt number along the hot wall (99 = 9801 outputs)2. Temperature along the Z = 0 . = 9801 outputs)3. X velocity along the Z = 0 . = 9801 outputs)4. Y velocity along the Z = 0 . = 9801 outputs)Latin hypercube samples (LHS) are generated using the python packagepyDOE [47]. The uniformly distributed LHS are transformed into normaldistribution using the inverse cumulative distribution function (ppf) of thestatistical functions module of scipy [48]. Separate sets of LHS are generatedfor training, validation and testing. Cases with 4, 8 and 16 strips are trainedwith 500 samples whereas, for 32 strips, 1000 samples are required. For eachcase, two different sets of 100 samples are used for validation and testing.The number of neurons in the input and output layers is specified by thenumber of inputs and outputs. The learning rate, regularization constant,optimizer and the number of hidden layers and neurons are highly problemspecific and are chosen so that both the training and validation error aresimultaneously minimized. The prediction accuracy is then checked on anunseen testing set. This overall procedure helps in fitting a DNN with lowbias and low variance [43]. The DNNs are implemented in the Python libraryTensorflow [49] with a high level API Keras [50]. Among the various opti-mizers available in Keras, the Adam optimizer [45] is found most suitable inthis work. Settings of Adam optimizer are as follows: learning rate of 10 − ,35 = 0 . β = 0 .
999 and ‘amsgrad’ option switched on. ReLU and identityare the activation functions for all the hidden and output layers respectively.Other hyperparameters specific to each of the four DNNs are as follows:1. Wall Nusselt number: λ = 0 . L h = 5, n h = 3002. Temperature: λ = 0 . L h = 4, n h = 3003. X velocity: λ = 0 . L h = 4, n h = 3004. Y velocity: λ = 0 . L h = 4, n h = 300where, λ is L regularization constant, L h is the number of hidden layersand n h is the number of neurons in each hidden layer. All the DNNs aretrained for 100 iterations on the entire dataset known as epochs. All thehyperparameters given above are tuned using the validation set with an ob-jective to minimize bias and variance. Figure 15 plots the loss versus epochsduring training for each of the four DNNs for the case with 4 strips. Lossesfor DNNs of 8, 12 and 16 strips are similar and hence are not plotted here.Since both the training and validation losses are close to each other, it showsthat the variance is low.Table 2 documents the training and testing errors for all the 16 DNNs: 4DNNs each for 4, 8, 16, and 32 strips. The relative average percent error isdefined as hundred times the L norm of the difference between true values(numerical simulation) and DNN estimates divided by the maximum absolutevalue of the output. The training and testing errors are small and closeenough thus, implying low bias and variance. For visual inspection, estimatesfrom the numerical simulation and the DNN are plotted together in fig. 16for the testing samples for the case of 4 strips. Each output is normalizedby subtracting its mean and dividing by its standard deviation and thus, is36on-dimensional. Ideally, all the points should lie on the Y = X line butdue to the interpolation error, some points are off the line. Since most of thepoints follow the expected trend of the Y = X line, it can be concluded thatthe neural network surrogate is accurate. (a) Hot Wall Nusselt Number (b) Temperature at Z = 0 . Z = 0 . Z = 0 . Figure 15: DNN Training and Validation Loss for 4 Strips
37 Strips 8 Strips 16 Strips 32 StripsTrain Test Train Test Train Test Train TestWall Nusselt No. 0.302 0.390 0.205 0.512 0.602 2.260 1.393 2.373Temperature 0.029 0.030 0.024 0.025 0.022 0.023 0.022 0.024X Velocity 1.114 1.144 0.864 0.903 0.704 0.731 0.644 0.722Y Velocity 0.372 0.387 0.297 0.311 0.248 0.260 0.233 0.230
Table 2: Relative Average Percent Error: DNN Training and Testing (a) Hot Wall Nusselt Number (b) Temperature at Z = 0 . Z = 0 . Z = 0 . Figure 16: Estimate from Numerical Simulation and DNN for 4 Strips .3.2. Nusselt Number (a) 4 Strips (b) 8 Strips(c) 16 Strips (d) 32 Strips Figure 17: Ra = 10 Hot Wall Nusselt Number: Difference between Stochastic Mean andDeterministic Value
Figure 17 plots the difference between the stochastic mean and the de-terministic value of the Nusselt number on the hot wall. Since this wallis subjected to the stochastic boundary condition, the demarcations of thestrips can be seen. For example, figs. 17a and 17b have four and eight stripsrespectively. Similar strips are also observed on the stochastic standard de-viation contours (fig. 18). The maximum deterministic value of the Nusselt39umber over the hot wall is 18.71 (fig. 5a). The maximum values of thedifference are 1.62, 1.91, 1.93 and 2.17 for 4, 8, 16 and 32 strips respectively.Hence, this shift in the stochastic mean from the deterministic value is com-parable to the input uncertainty of 10%. On the other hand, the stochasticstandard deviation is of the order of 10 − . Thus, it is seen that the shift inmean is more pronounced than its deviation similar to case A. (a) 4 Strips (b) 8 Strips(c) 16 Strips (d) 32 Strips Figure 18: Ra = 10 Hot Wall Nusselt Number: Stochastic Standard Deviation .3.3. Velocity and Temperature (a) 4 Strips (b) 8 Strips(c) 16 Strips (d) 32 Strips Figure 19: Ra = 10 Z Midplane Temperature: Difference between Stochastic Mean andDeterministic Value
Figures 19 to 21 plot the difference between the stochastic mean and thedeterministic value of the temperature, X and Y velocities respectively alongthe Z = 0 . −
4% of thedeterministic value. For velocities the difference is around 4 −
5% of thedeterministic value. This is smaller than that of the Nusselt number which is8 − (a) 4 Strips (b) 8 Strips c) 16 Strips (d) 32 Strips Figure 20: Ra = 10 Z Midplane X-Velocity: Difference between Stochastic Mean andDeterministic Value (a) 4 Strips (b) 8 Strips c) 16 Strips (d) 32 Strips Figure 21: Ra = 10 Z Midplane Y-Velocity: Difference between Stochastic Mean andDeterministic Value (a) Temperature (b) X-Velocity (c) Y-Velocity
Figure 22: Stochastic Standard Deviation at Z Mid-plane Ra = 10 (32 Strips)
6. Conclusions
This paper presents the input uncertainty propagation results for a threedimensional natural convection problem in a differentially heated cubicalenclosure with two different cases. Case A assumes that the boundaryconditions are uniform. Thus, uncertainty is introduced in the two non-dimensional parameters (Rayleigh and Prandtl numbers). For this case, the44olynomial chaos expansion (PCE) method is used as a surrogate modelwith stochastic collocation to estimate the PCE coefficients. Case B dealswith non-uniform stochastic boundary condition with deterministic materialproperties of the fluid. Since the temperature difference between the oppositewalls drives the natural convection flow, the cold wall is held at a constanttemperature and uncertainty is introduced in the hot wall temperature. Adeep neural network based surrogate model is used here for estimating theoutput statistics with Monte Carlo method.It is observed that the mean value of an output parameter averaged overthe stochastic variation of the input can be different from the deterministicoutput value. The difference normalized by the deterministic value is definedas the relative shift of mean. Table 3 summarizes the shift of mean relativeto deterministic values of each of the output. Mean Rayleigh numbers of10 and 10 are considered in case A whereas, case B deals only with amean of 10 . The relative values of input standard deviation with respectto mean is 2% and 3.33% for case A and case B respectively. It can be seenthat for case A, the stochastic effect is negligible in the Nusselt number andtemperature whereas, it is of similar order as input in the velocities. Onthe other hand, for case B, the stochastic effect is much higher in Nusseltnumber. For both the cases, the standard deviation is much lower than theshift of mean. In general, it is concluded that the effect of uncertainties in theboundary conditions are much higher than the uncertainties in the Rayleighand Prandtl numbers with uniform boundary conditions.45 nput Nusselt No. Temperature Velocities Case A Ra = 10 σ = 2% µ Ra = 10 σ = 2% µ Ra = 10 σ = 3 . µ
12 % 4 % 5 %
Table 3: Shift of Mean Relative to Deterministic Values
This paper also demonstrates the use of deep neural network for uncer-tainty quantification in natural convection problem for a high dimensionalstochastic problem. The novel approach of dividing the boundary surfaceinto domains and treating each domain value as a stochastic input is demon-strated here for the first time (to the best of our knowledge). This approachtries to mimic the experiments, when multiple feedback control systems areused to impose the boundary conditions. Note that although in this work, it isassumed that the domain values are independent and identically distributedvariables, this assumption can be relaxed during sampling if additional dataof the feedback control system is available. The approach can also be usedfor other class of problems like optimization and inverse problems.
Acknowledgments
This work was funded by the Digital Manufacturing and Design Inno-vation Institute with support in part by the U.S. Department of the Army.Any opinions, findings, and conclusions or recommendations expressed in thismaterial are those of the author(s) and do not necessarily reflect the viewsof the Department of the Army. 46 eferencesReferences [1] G. De Vahl Davis, Natural convection of air in a square cavity: a benchmark numerical solution, International Journal for Numerical Methodsin Fluids 3 (3) (1983) 249–264.[2] T. Fusegi, J. Hyun, K. Kuwahara, B. Farouk, A numerical study ofthree-dimensional natural convection in a differentially heated cubicalenclosure, International Journal of Heat and Mass Transfer 34 (6) (1991)1543–1557.[3] P. Le Qu´er´e, Accurate solutions to the square thermally driven cavityat high Rayleigh number, Computers & Fluids 20 (1) (1991) 29–41.[4] Y. Hu, Y. Li, C. Wu, S. Li, M. Li, Flow pattern and heat transfer inRayleigh-B´enard convection of cold water near its density maximum ina rectangular cavity, International Journal of Heat and Mass Transfer107 (2017) 1065–1075.[5] Y. Li, H. Zhang, L. Zhang, C. Wu, Three-dimensional numerical sim-ulation of double-diffusive Rayleigh–B´enard convection in a cylindri-cal enclosure of aspect ratio 2, International Journal of Heat and MassTransfer 98 (2016) 472–483.[6] C. Shu, H. Ding, K. Yeo, Local radial basis function-based differentialquadrature method and its application to solve two-dimensional incom-pressible Navier–Stokes equations, Computer Methods in Applied Me-chanics and Engineering 192 (7-8) (2003) 941–954.477] S. Yigit, R. Poole, N. Chakraborty, Effects of aspect ratio on lami-nar Rayleigh–B´enard convection of power-law fluids in rectangular en-closures: A numerical investigation, International Journal of Heat andMass Transfer 91 (2015) 1292–1307.[8] T. Icoz, Y. Jaluria, Design of cooling systems for electronic equipmentusing both experimental and numerical inputs, Journal of ElectronicPackaging 126 (4) (2004) 465–471.[9] A. Sharma, K. Velusamy, C. Balaji, S. Venkateshan, Conjugate turbu-lent natural convection with surface radiation in air filled rectangularenclosures, International Journal of Heat and Mass Transfer 50 (3-4)(2007) 625–639.[10] M. Herrlin, C. Belady, Gravity-assisted air mixing in data centers andhow it affects the rack cooling effectiveness, in: Thermal and Thermo-mechanical Phenomena in Electronics Systems, 2006. ITHERM’06. TheTenth Intersociety Conference on, IEEE, 5–pp, 2006.[11] S. Patankar, Numerical heat transfer and fluid flow, CRC Press, 1980.[12] D. Xiu, G. Karniadakis, Modeling uncertainty in steady state diffusionproblems via generalized polynomial chaos, Computer Methods in Ap-plied Mechanics and Engineering 191 (43) (2002) 4927–4948.[13] D. Xiu, G. Karniadakis, Modeling uncertainty in flow simulations viageneralized polynomial chaos, Journal of Computational Physics 187 (1)(2003) 137–167. 4814] O. Maitre, O. Knio, H. Najm, R. Ghanem, A stochastic projectionmethod for fluid flow: I. basic formulation, Journal of ComputationalPhysics 173 (2) (2001) 481–511.[15] O. Maitre, M. Reagan, B. Debusschere, H. Najm, R. Ghanem, O. Knio,Natural convection in a closed cavity under stochastic non-Boussinesqconditions, SIAM Journal on Scientific Computing 26 (2) (2004) 375–394.[16] D. Venturi, M. Choi, G. Karniadakis, Supercritical quasi-conductionstates in stochastic Rayleigh–B´enard convection, International Journalof Heat and Mass Transfer 55 (13-14) (2012) 3732–3743.[17] B. Ganapathysubramanian, N. Zabaras, Sparse grid collocation schemesfor stochastic natural convection problems, Journal of ComputationalPhysics 225 (1) (2007) 652–685.[18] M. Carnevale, F. Montomoli, A. DAmmaro, S. Salvadori, F. Martelli,Uncertainty quantification: A stochastic method for heat transfer pre-diction using LES, Journal of Turbomachinery 135 (5) (2013) 051021.[19] P. Marepalli, J. Murthy, B. Qiu, X. Ruan, Quantifying uncertaintyin multiscale heat conduction calculations, Journal of Heat Transfer136 (11) (2014) 111301.[20] N. Fajraoui, M. Fahs, A. Younes, B. Sudret, Analyzing natural convec-tion in porous enclosure with polynomial chaos expansions: Effect ofthermal dispersion, anisotropic permeability and heterogeneity, Interna-tional Journal of Heat and Mass Transfer 115 (2017) 205–224.4921] K. Fezi, M. Krane, Uncertainty quantification in modeling metal alloysolidification, Journal of Heat Transfer 139 (8) (2017) 082301.[22] S. Hosder, R. Walters, R. Perez, A non-intrusive polynomial chaosmethod for uncertainty propagation in CFD simulations, in: 44th AIAAAerospace Sciences Meeting and Exhibit, 891, 2006.[23] D. Kumar, M. Raisee, C. Lacor, An efficient non-intrusive reduced basismodel for high dimensional stochastic problems in CFD, Computers &Fluids 138 (2016) 67–82.[24] S. Smolyak, Quadrature and interpolation formulas for tensor productsof certain classes of functions, in: Soviet Math. Dokl., vol. 4, 240–243,1963.[25] F. Heiss, V. Winschel, Likelihood approximation by numerical integra-tion on sparse grids, Journal of Econometrics 144 (1) (2008) 62–80.[26] R. Caflisch, Monte carlo and quasi-monte carlo methods, Acta Numerica7 (1998) 1–49.[27] K. Hornik, M. Stinchcombe, H. White, Multilayer feedforward networksare universal approximators, Neural Networks 2 (5) (1989) 359–366.[28] J. Schmidhuber, Deep learning in neural networks: An overview, NeuralNetworks 61 (2015) 85–117.[29] A. Gholami, H. Bonakdari, A. Zaji, A. Akhtari, Simulation of openchannel bend characteristics using computational fluid dynamics and50rtificial neural networks, Engineering Applications of ComputationalFluid Mechanics 9 (1) (2015) 355–369.[30] S. Sablani, A. Kacimov, J. Perret, A. Mujumdar, A. Campo, Non-iterative estimation of heat transfer coefficients using artificial neuralnetwork models, International Journal of Heat and Mass Transfer 48 (3-4) (2005) 665–679.[31] A. Santra, N. Chakraborty, S. Sen, Prediction of heat transfer due topresence of copper–water nanofluid using resilient-propagation neuralnetwork, International Journal of Thermal Sciences 48 (7) (2009) 1311–1318.[32] B. Cz´el, K. Woodbury, G. Gr´of, Simultaneous estimation oftemperature-dependent volumetric heat capacity and thermal conduc-tivity functions via neural networks, International Journal of Heat andMass Transfer 68 (2014) 1–13.[33] Z. Zhang, K. Duraisamy, Machine learning methods for data-driven tur-bulence modeling, in: 22nd AIAA Computational Fluid Dynamics Con-ference, 2460, 2015.[34] R. Tripathy, I. Bilionis, Deep UQ: Learning deep neural network sur-rogate models for high dimensional uncertainty quantification, arXivpreprint arXiv:1802.00850 .[35] Y. Khoo, J. Lu, L. Ying, Solving parametric PDE problems with artifi-cial neural networks, arXiv preprint arXiv:1707.03351 .5136] D. Zhang, L. Lu, L. Guo, G. Karniadakis, Quantifying total uncer-tainty in physics-informed neural networks for solving forward and in-verse stochastic problems, arXiv preprint arXiv:1809.08327 .[37] S. Shahane, N. Aluru, P. Ferreira, S. Kapoor, S. Vanka, Finite VolumeSimulation Framework for Die Casting with Uncertainty Quantification,arXiv preprint arXiv:1810.08572 .[38] F. Harlow, J. Welch, Numerical calculation of time-dependent viscousincompressible flow of fluid with free surface, The physics of fluids 8 (12)(1965) 2182–2189.[39] D. Xiu, G. Karniadakis, The Wiener–Askey polynomial chaos forstochastic differential equations, SIAM Journal on Scientific Comput-ing 24 (2) (2002) 619–644.[40] N. Wiener, The homogeneous chaos, American Journal of Mathematics60 (4).[41] R. Smith, Uncertainty quantification: theory, implementation, and ap-plications, 2013.[42] S. Marelli, B. Sudret, UQLab: A framework for uncertainty quantifica-tion in Matlab, in: Vulnerability, Uncertainty, and Risk: Quantification,Mitigation, and Management, 2554–2563, 2014.[43] I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, , 2016.5244] Y. Chauvin, D. Rumelhart, Backpropagation: theory, architectures, andapplications, Psychology Press, 2013.[45] D. Kingma, J. Ba, Adam: A method for stochastic optimization, arXivpreprint arXiv:1412.6980 .[46] I. Sobol, Sensitivity estimates for nonlinear mathematical models, Math-ematical modelling and computational experiments 1 (4) (1993) 407–414.[47] pyDOE: The experimental design package for python, URL https://pythonhosted.org/pyDOE/ , 2018.[48] Statistical functions (scipy.stats), URL https://docs.scipy.org/doc/scipy/reference/stats.html , 2018.[49] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro,G. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfel-low, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser,M. Kudlur, J. Levenberg, D. Man´e, R. Monga, S. Moore, D. Mur-ray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Tal-war, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Vi´egas, O. Vinyals,P. Warden, M. Wattenberg, M. Wicke, Y. Yu, X. Zheng, TensorFlow:Large-Scale Machine Learning on Heterogeneous Systems, URL , software available from tensorflow.org, 2015.[50] F. Chollet, et al., Keras, https://keras.iohttps://keras.io