HYPAD-UQ: A Derivative-Based Uncertainty Quantification Method Using a Hypercomplex Finite Element Method

Balcer, Matthew; Aristizabal, Mauricio; Rincon-Tabares, Juan-Sebastian; Montoya, Arturo; Restrepo, David; Millwater, Harry

doi:10.1115/1.4062459

Abstract

A derivative-based uncertainty quantification (UQ) method called HYPAD-UQ that utilizes sensitivities from a computational model was developed to approximate the statistical moments and Sobol' indices of the model output. Hypercomplex automatic differentiation (HYPAD) was used as a means to obtain accurate high-order partial derivatives from computational models such as finite element analyses. These sensitivities are used to construct a surrogate model of the output using a Taylor series expansion and subsequently used to estimate statistical moments (mean, variance, skewness, and kurtosis) and Sobol' indices using algebraic expansions. The uncertainty in a transient linear heat transfer analysis was quantified with HYPAD-UQ using first-order through seventh-order partial derivatives with respect to seven random variables encompassing material properties, geometry, and boundary conditions. Random sampling of the analytical solution and the regression-based stochastic perturbation finite element method were also conducted to compare accuracy and computational cost. The results indicate that HYPAD-UQ has superior accuracy for the same computational effort compared to the regression-based stochastic perturbation finite element method. Sensitivities calculated with HYPAD can allow higher-order Taylor series expansions to be an effective and practical UQ method.

1 Introduction

Conducting uncertainty quantification (UQ) with high-fidelity models is challenging due to computational demands. A cheap-to-evaluate surrogate model such as polynomial chaos expansion [1], Gaussian process [2], or the Taylor series expansion [3] can be used as a substitute. The Taylor series expansion has shown to be a viable surrogate model for many types of physics [4]. However, it is mainly used when the variance in random variables is low to ensure high accuracy with a first-order or second-order expansion. Expanding the Taylor series to higher orders can improve the accuracy [5], but sensitivities higher than second-order are often difficult to compute with traditional methods such as finite differences or symbolic differentiation. Numerical differentiation methods have been used to construct higher-order Taylor series expansions such as the adjoint method [6] and modified forward finite difference (ModFFD) [7]. However, these methods have disadvantages when calculating higher-order derivatives, specifically, orders past second-order. In adjoint differentiation, partial derivatives need to be derived for each mathematical formulation and for each variable. Higher-order adjoint-based derivatives exist, however, the formulations become intractable as the order increases. Higher-order partial derivatives can be calculated with ModFFD in a straight-forward manner, however, the step size highly influences the accuracy of the higher-order derivatives. Hence, ModFFD can become computationally expensive when a step size convergence study is conducted and accuracy is not guaranteed. For these reasons, this work utilized higher-order partial derivatives calculated with an accurate, efficient, and easy-to-use numerical differentiation method called hypercomplex automatic differentiation (HYPAD) to construct Taylor series expansions of finite element analysis (FEA) nodal outputs in order to conduct UQ, known as HYPAD-UQ.

The partial derivatives that make up the Taylor series expansion can be calculated by a variety of methods. Symbolic derivatives of governing FEA equations have been developed for many different physics including elasticity, plasticity, heat transfer, dynamics, and many others [4,8]. This method is very efficient [9], however, derivative expressions for new mathematical models, which are not necessarily finite element-based, need to be derived for each variable of interest. Alternatively, adjoint-based differentiation has been successful in calculating first-order and second-order derivatives very efficiently [6]. Reference [10] used adjoint-based derivatives to approximate first-order and second-order Taylor series expansions of the expected value, variance, and skewness of the uncollided particle leakage. However, similar to symbolic differentiation, adjoint-based derivative expressions for each variable of interest need to be derived for different mathematical formulations.

Numerical differentiation methods exist which do not require new mathematical formulation. One widely used method is finite difference which is nonintrusive, i.e., source code does not need to be altered in any way, and only requires evaluations of the model at multiple input parameters. One drawback is that the finite difference method is subject to subtraction cancelation error and truncation error [11]. Therefore, a convergence study for the optimal step size is needed which makes this method computationally demanding, especially for higher-order derivatives. The complex-step differentiation method is analogous to the forward finite difference method in that the variable of interest is perturbed along an imaginary direction in order to calculate first-order partial derivatives [12]. This method yields the most accurate numerical derivative of the model. Although source code must be modified to use complex algebra, no new formulation is required for the derivative expressions. Complex-step differentiation was used in Ref. [13] to construct a first-order Taylor series to approximate variance of residual stress in a thick-walled sphere. Arbitrary-order partial derivatives can be calculated by generalizing the complex-step differentiation method using HYPAD.

Hypercomplex automatic differentiation is a highly accurate numerical differentiation method that uses hypercomplex algebra within the formulation of a mathematical model to calculate arbitrary-order partial derivatives. Analogous to the forward finite difference method, derivatives are calculated in HYPAD by perturbing variables of interest along independent imaginary or nonreal directions using a hypercomplex algebra [14–16]. The two most common HYPAD methods were developed with multicomplex numbers [17] and multidual numbers [18] which are generalizations of complex numbers and dual numbers, respectively, to any number of imaginary directions. Unlike the finite difference method, there is no subtractive cancelation error; therefore, step sizes can be very small to yield numerically accurate sensitivities.

Hypercomplex automatic differentiation has been implemented in FEA, known as the hypercomplex finite element method (ZFEM) [19]. A single nth-order derivative and all combinations of lower-order derivatives with respect to the perturbed variables can be computed in one analysis using multicomplex or multidual algebra. Hence, the number of model evaluations needed to compute a complete nth-order Taylor series is equivalent to the number of nth-order derivatives in the Taylor series. The block forward substitution method has been applied to the solution of the linear finite element system of equations to allow the reuse of the real solution for multiple derivative computations [20]. However, recalculation of the imaginary parts of hypercomplex matrices is still required when calculating multiple higher-order derivatives. A recently developed hypercomplex algebra called order truncated imaginary (OTI) algebra was developed for the purpose of calculating all combinations of partial derivatives with respect to multiple variables in a single analysis without any repetitive calculations. This method has shown to yield a dramatic increase in efficiency compared to multidual-based differentiation [16]. Further details on HYPAD are given in Sec. 2.2.

The Taylor series expansion of an output can be computed nonintrusively using linear regression. Reference [21] showed that Taylor series expansions approximated by regression can yield accurate central moments for various single variable functions. They also showed that higher-order expansions were necessary to achieve a sufficient estimate of central moments as the range of variance in the random variable increased. A disadvantage of linear regression, however, is the potential need for a large number of samples to obtain sufficient accuracy of the regression model.

Conducting global sensitivity analysis to identify random variables that significantly contribute to the uncertainty in the output is important in engineering design. This can be accomplished by computing the Sobol' variance-based sensitivity indices, herein referred to as Sobol' indices, of the output which measures the contribution of the random variables, or interactions between random variables, to the variance of the output [22]. Calculating these Sobol' indices with random sampling can become impractical to compute due to high computational demands. Therefore, surrogate models, such as polynomial chaos expansion [23], can be used in place of the high-fidelity model to approximate Sobol' indices. Computing Sobol' indices with the Taylor series are rarely employed in the literature due to concerns about the accuracy and the need for higher-order derivatives to compute interaction effects. However, previous work has shown the effectiveness of the Taylor series in global sensitivity analysis. In Ref. [24], main effect Sobol' indices were computed from a first-order Taylor series expansion constructed with adjoint-based derivatives. They showed that out of 16 random variables, about 5–8 were important for structural reliability of a shell structure. Reference [25] showed that the Taylor series coincides with the analysis of variance decomposition when the model is multilinear.

HYPAD-UQ is not the only method that utilizes derivatives for UQ. Reference [26] computed multiple first-order Taylor series expansions at different expansion points to create a piecewise linear surrogate model. Derivatives evaluated at multiple sampling points are also used in gradient enhanced surrogate modeling [27,28] which improves the accuracy of sampling-based surrogate models and potentially allows for fewer sampling points. While these methods are attractive because only lower-order derivatives are required, selection and distribution of the sample points are an open question. In optimization problems, Ref. [29] proposed a methodology to approximate the computation of first-order sensitivities using the “hyperdifferential sensitivity analysis” method. This method approximates first-order sensitivities of partial differential equation-constrained optimization solution with respect to random inputs. The advantages of this method are that it can be used to compute local and global sensitivities and that it is highly parallelizable; however, the method requires appropriate selection of sample points and thus multiple solutions to the optimization problem. HYPAD-UQ, the method presented here, only requires one FEA model evaluation at the mean of the random variables and is general for arbitrary-order derivatives.

This paper is organized as follows: Sobol' indices and central moments of Taylor series expansions for any random variable distribution type are derived in Sec. 2.1. An overview of HYPAD is given in Sec. 2.2. The implementation of HYPAD in finite elements is given in Sec. 2.3. Section 3 shows an application of HYPAD-UQ on a transient heat transfer analysis in a fin where accuracy and computational cost is compared to the linear regression-based stochastic perturbation method. A discussion on the effect of error and computational cost of HYPAD-UQ is presented in Sec. 4. Finally, a conclusion of this work is provided in Sec. 5.

2 Methodology

HYPAD-UQ creates a Taylor series expansion of an output of a mathematical model in order to approximate probability densities, central moments, and Sobol' indices, see Sec. 2.1. This Taylor series expansion about the mean of the output is constructed from partial derivatives of the output with respect to the random variables using HYPAD, see Sec. 2.2. The implementation of OTI-based HYPAD in finite elements in the context of a thermal analysis is given in Sec. 2.3.

2.1 Taylor Series Expansion of Central Moments and Sobol' Indices.

It is convenient to derive the Taylor series expansions of variance using the high-dimensional model representation (HDMR) [30] so that variance-based sensitivity indices, known as Sobol' indices, can be readily derived. Sobol' indices and the expected value of the first-order through third-order Taylor series expansions are derived below for any distribution of the random variables. Skewness and kurtosis of Taylor series expansions are also shown for any distribution type through second-order. Skewness and kurtosis of third-order expansions have been derived and were used in Sec. 3; however, they are not shown due to the large number of terms.

A function

f (x)

of independent random variables

x = [x_{1}, \dots, x_{r}]

can be decomposed into orthogonal component functions that depend on combinations of variables in x. In this work, lower case symbols in bold are vectors and upper case symbols in bold are matrices. The HDMR of

f (x)

is expressed as

f (x) = f_{0} + \sum_{i = 1}^{r} f_{i} (x_{i}) + \sum_{1 \leq i < j \leq r} f_{i j} (x_{i}, x_{j}) + \dots + f_{12 \dots r} (x_{1}, x_{2}, \dots, x_{r})

(1)

where r is the number of variables in x, f₀ is a constant, and

f_{i j \dots}

are component functions that depend only on the variables denoted in the subscript [30]. The component functions in Eq. can be expressed in terms of conditional expectations

E [\cdot]

⁠, as

f_{0} = E [f (x)] = \int f (x) \cdot p_{x} d x

(2)

f_{i} = E [f (x) | x_{i}] - f_{0} = \int f (x) \cdot p_{\sim i} d x_{\sim i} - f_{0}

(3)

f_{i j} = E [f (x) | x_{i}, x_{j}] - f_{i} - f_{j} - f_{0} = \int f (x) \cdot p_{\sim i, j} d x_{\sim i, j} - f_{i} - f_{j} - f_{0}

(4)

and so on [22]. The subscript $\sim i, j, \dots$ represents the exclusion of the random variables denoted by the indices and $p_{x}$ is the joint probability distribution of all random variables in x.

Assuming

f (x)

is square integrable and all random variables in x are independent, the variance of Eq. gives the contribution of variance that depends on a single random variable and combinations of random variables, known as the analysis of variance decomposition

V = \sum_{i = 1}^{r} V_{i} + \sum_{1 \leq i < j \leq r} V_{i j} + \dots + V_{12 \dots r}

(5)

where

V_{i} = \int f_{i}^{2} \cdot p_{i} d x_{i}

(6)

V_{i j} = \int f_{i j}^{2} \cdot p_{i j} d x_{i} d x_{j}

(7)

and so on [22,31]. p_ij denotes the joint probability distribution of x_i and x_j. The partial variances in Eq. can be normalized by the total variance V to yield

1 = \sum_{i = 1}^{r} S_{i} + \sum_{1 \leq i < j \leq r} S_{i j} + \dots + S_{12 \dots r}

(8)

where the Sobol' indices are expressed as

S_{i} = \frac{V_{i}}{V}

(9)

S_{i j} = \frac{V_{i j}}{V}

(10)

and so on.

The following derives Sobol' indices of first-order through third-order Taylor series expansions for any distribution of random variables. The nth-order Taylor series expansion Y_n around the mean of the random variables

μ

is expressed as

f (x) \approx Y_{n} (x) = f (μ) + \sum_{j_{1} = 1}^{r} D_{j_{1}}^{(1)} Δ_{j_{1}} + \frac{1}{2} \underset{j_{1}, j_{2} = 1}{\sum^{r}} D_{j_{1} j_{2}}^{(2)} Δ_{j_{1}} Δ_{j_{2}} + \dots + \frac{1}{n!} \sum_{j_{1}, \dots, j_{n} = 1}^{r} D_{j_{1} \dots j_{n}}^{(n)} Δ_{j_{1}} \dots Δ_{j_{n}}

(11)

where

Δ_{j_{k}} = (x_{j_{k}} - μ_{j_{k}})

and

D_{j_{1} \dots j_{s}}^{(n)}

denote the nth-order partial derivative of

f (x)

with respect to

x_{j_{1}}, \dots, x_{j_{n}}

evaluated at

μ

⁠, i.e.,

D_{j_{1} \dots j_{s}}^{(n)} = \partial^{n} f (μ) / \partial x_{j_{1}} \dots \partial x_{j_{s}}

⁠. The total number of terms in the nth-order Taylor series expansion of r variables is

d = \frac{(r + n)!}{r! n!}

(12)

and the number of terms in the Taylor series that contain pth-order derivatives is

d_{p} = \frac{(p + r - 1)!}{(r - 1)! p!}

(13)

Using Eqs. – and noting that

E [Δ_{k}^{l}] = μ_{l k}

is the lth central moment of random variable x_k, defined as

μ_{l k} = E [{(x_{k} - μ_{k})}^{l}] = \int {(x_{k} - μ_{k})}^{l} p_{k} d x_{k},

(14)

the HDMR of the first-order and second-order Taylor series expansions is

Y_{1} (x) = \underset{f_{0}}{\underset{︸}{f (μ)}} + \sum_{i = 1}^{r} \underset{f_{i}}{\underset{︸}{D_{i}^{(1)} Δ_{i}}}

(15)

Y_{2} (x) = \underset{f_{0}}{\underset{︸}{f (μ) + \frac{1}{2} \sum_{j_{1} = 1}^{r} D_{j_{1} j_{1}}^{(2)} μ_{2 j_{1}}}} + \sum_{i = 1}^{r} \underset{f_{i}}{\underset{︸}{{D_{i}^{(1)} Δ_{i} + \frac{1}{2} D_{i i}^{(2)} (Δ_{i}^{2} - μ_{2 i})}}} + \sum_{1 \leq i < j \leq r} \underset{f_{i j}}{\underset{︸}{D_{i j}^{(2)} Δ_{i} Δ_{j}}}

(16)

where the indices i and j represent indices of specific variables as in Eq. and index j_k represents a generic index. The expected value of the first-order and second-order Taylor series expansions is

E [Y_{1}] = f (μ)

(17)

E [Y_{2}] = E [Y_{1}] + \sum_{j_{1} = 1}^{r} \frac{1}{2} D_{j_{1} j_{1}}^{(2)} μ_{2 j_{1}}

(18)

Sobol' indices of Taylor series expansions are derived by integrating Y_n using Eqs. and to yield the partial variances. The total variance of Y_n is calculated by the summation of all partial variances, as in Eq. . The resulting Sobol' indices of first-order and second-order Taylor series expansions are

First - order S_{i} [Y_{1}] = \frac{{(D_{i}^{(1)})}^{2} μ_{2 i}}{V [Y_{1}]}

(19)

Second - order S_{i} [Y_{2}] = \frac{{(D_{i}^{(1)})}^{2} μ_{2 i} + \frac{1}{4} {(D_{i i}^{(2)})}^{2} (μ_{4 i} - μ_{2 i}^{2}) + D_{i}^{(1)} D_{i i}^{(2)} μ_{3 i}}{V [Y_{2}]}

(20)

S_{i j} [Y_{2}] = \frac{{(D_{i j}^{(2)})}^{2} μ_{2 i} μ_{2 j}}{V [Y_{2}]}

(21)

where the total variance of first-order and second-order Taylor series expansions is

V [Y_{1}] = \sum_{j_{1} = 1}^{r} {(D_{j_{1}}^{(1)})}^{2} μ_{2 j_{1}}

(22)

V [Y_{2}] = V [Y_{1}] + \sum_{j_{1} = 1}^{r} {\frac{1}{4} {(D_{j_{1} j_{1}}^{(2)})}^{2} (μ_{4 j_{1}} - μ_{2 j_{1}}^{2}) + D_{j_{1}}^{(1)} D_{j_{1} j_{1}}^{(2)} μ_{3 j_{1}}} + \sum_{1 \leq j_{1} < j_{2} \leq r} {(D_{j_{1} j_{2}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}}

(23)

respectively. Note that second-order mixed partial derivatives are needed to approximate any interaction effects with the Taylor series. See Appendix for the expected value and partial variances of a third-order Taylor series expansion.

Approximations for higher-order statistical moments are useful for characterizing the distribution of the output. Skewness and kurtosis are defined as

skewness = \frac{μ_{3} (f (x))}{μ_{2} {(f (x))}^{3 / 2}}

(24)

kurtosis = \frac{μ_{4} (f (x))}{μ_{2} {(f (x))}^{2}}

(25)

where

μ_{l} (f (x))

is the lth central moment of a function which is defined in Eq. by substituting x_k with

f (x)

⁠. Similarly, the lth central moment of an nth-order Taylor series expansion is found by substituting x_k with Y_n, represented as

μ_{l} (f (x)) \approx μ_{l} (Y_{n} (x)) = E [{(Y_{n} (x) - E [Y_{n} (x)])}^{l}]

(26)

The third central moment of a first-order Taylor series expansion is

μ_{3} (Y_{1}) = \underset{j_{1} = 1}{\sum^{r}} {(D_{j_{1}}^{(1)})}^{3} μ_{3 j_{1}}

(27)

and the fourth central moment of a first-order Taylor series expansion is

μ_{4} (Y_{1}) = \underset{j_{1} = 1}{\sum^{r}} {(D_{j_{1}}^{(1)})}^{4} μ_{4 j_{1}} + \underset{j_{1} \neq j_{2}}{\sum^{r}} 3 {(D_{j_{1}}^{(1)})}^{2} {(D_{j_{2}}^{(1)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}}

(28)

See Appendix for third and fourth central moments of the second-order Taylor series expansion.

2.2 Hypercomplex-Based Automatic Differentiation.

The methodology of HYPAD is introduced with multidual (also known as “hyperdual”) algebra in Sec. 2.2.1. A summary of OTI-based HYPAD is given in Sec. 2.2.2, which is an extension of multidual-based differentiation that is more efficient when calculating multiple higher-order derivatives.

2.2.1 Differentiation With Complex and Hypercomplex Algebra.

First-order partial derivatives can be calculated using complex algebra, known as the complex-step method [12]. It is a well-established method and used in many areas such as aerospace [32], computer graphics [33], and astrophysics [34]. The complex-step method is similar to the forward finite difference method except the perturbation is applied to the variable of interest along an imaginary direction instead of the real direction. The first-order partial derivative can be obtained from a first-order Taylor series expansion of a function

f (x)

⁠, with a perturbation applied to the variable of interest x_k in the imaginary direction. The resulting expression for a first-order partial derivative with respect to x_k is

D_{x_{k}}^{(1)} = \frac{\partial f (x)}{\partial x_{k}} = \frac{I m [f (x + i h e_{k})]}{h} + O (h^{2})

(29)

where $i^{2} = - 1$ ⁠, h is step size of the perturbation, and $e_{k}$ is a unit vector with a value of one in the kth position. This unit vector which indicates that the perturbation is applied to the variable of interest x_k. A step size h of the order of $10^{- 10}$ times the variable of interest usually yields an accurate derivative value for the given model [35]. The complex-step method does not add any extra error within the formulation. The accuracy of the derivative is completely dependent on the numerical algorithm in which it is employed.

First-order partial derivatives are computed with the complex-step method by applying a perturbation to x_k in the imaginary direction as $x_{k} + i h$ ⁠. The function is then evaluated with complex algebra. The real part contains the function output and the derivative is extracted from the imaginary part of the output, following Eq. (29).

Higher-order partial derivatives can be derived in a similar manner by applying perturbations to variables of interest along multiple imaginary or nonreal directions. Two subsets of hypercomplex algebra have been developed to achieve this and both methods compute accurate arbitrary-order derivatives in a similar manner. Reference [17] used multicomplex numbers, which is an extension of complex numbers to multiple imaginary directions. A more efficient method uses a generalization of dual numbers, which contains a nonreal part represented by $ϵ^{2} = 0$ ⁠, called hyperdual numbers, herein referred to as multidual numbers to emphasize the parallel construction with multicomplex numbers [18]. Multidual-based differentiation has no truncation error and is step-size independent, unlike the multicomplex method. However, a small enough step size in multicomplex should yield equivalent derivative values. Hypercomplex differentiation can be implemented in source code by altering real-valued algebra with a hypercomplex algebra. External libraries have been developed that overload real-valued algebraic operations with hypercomplex algebra [14,16,18,36,37] and the Julia programing language intrinsically supports multidual algebra [38].

Derivatives are calculated with multidual algebra by perturbing the n number of variables of interest along independent imaginary directions

ϵ_{1}, \dots, ϵ_{n}

⁠. For example, the second-order mixed partial derivative of

f (x)

with respect to x_i and x_j can be calculated by applying perturbations along the ϵ₁ and ϵ₂ directions. The resulting perturbation is represented as

x + ϵ_{1} h_{1} e_{i} + ϵ_{2} h_{2} e_{j} + 0 ϵ_{1} ϵ_{2}

⁠, where

e_{i}

and

e_{j}

are unit vectors that contain one in the ith and jth positions, respectively, in x which indicate the variables being perturbed and 0 is a vector of zeros of size r. The function is then evaluated with these perturbations using multidual algebra and the first-order and second-order partial derivatives are extracted from the imaginary parts of the output as [39]

D_{x_{i}}^{(1)} = \frac{\partial f (x)}{\partial x_{i}} = \frac{I m_{ϵ_{1}} [f (x + ϵ_{1} h_{1} e_{i} + ϵ_{2} h_{2} e_{j} + 0 ϵ_{1} ϵ_{2})]}{h_{1}},

(30)

D_{x_{j}}^{(1)} = \frac{\partial f (x)}{\partial x_{j}} = \frac{I m_{ϵ_{2}} [f (x + ϵ_{1} h_{1} e_{i} + ϵ_{2} h_{2} e_{j} + 0 ϵ_{1} ϵ_{2})]}{h_{2}},

(31)

D_{x_{i} x_{j}}^{(2)} = \frac{\partial^{2} f (x)}{\partial x_{i} \partial x_{j}} = \frac{I m_{ϵ_{1} ϵ_{2}} [f (x + ϵ_{1} h_{1} e_{i} + ϵ_{2} h_{2} e_{j} + 0 ϵ_{1} ϵ_{2})]}{h_{1} h_{2}}

(32)

Arbitrary-order partial derivatives with respect to n variables can be computed in a similar manner by applying perturbations along n independent imaginary directions. The function is evaluated with a hypercomplex algebra and partial derivatives are extracted from the imaginary parts of

f (x)

by

D_{x_{i} \dots x_{k}}^{(n)} = \frac{\partial^{n} f (x)}{\partial x_{i} \dots \partial x_{k}} = \frac{I m_{ϵ_{1} \dots ϵ_{n}} [f (x + ϵ_{1} h_{1} e_{i} + \dots + ϵ_{n} h_{n} e_{k} + \dots + 0 ϵ_{1} \dots ϵ_{n})]}{h_{1} \dots h_{n}}

(33)

where h_i is step sizes and $e_{k}$ is a unit vector with a value of one in the kth position [39]. Note that all combinations of lower-order partial derivatives can be extracted from the same function evaluation, similar to Eqs. (30)–(32). This is accomplished by extracting the imaginary parts of the respective variables of interest and dividing by the corresponding step sizes.

2.2.2 Hypercomplex Automatic Differentiation With Order Truncated Imaginary Algebra.

Multidual numbers are an effective method to compute higher-order derivatives. However, they increase in size by $2^{n}$ each time a new basis (⁠ $ϵ_{1}, \dots, ϵ_{n}$ ⁠) is added, where n is the total number of bases. A multidual number with n bases is capable of computing a single nth-order derivative and all combinations of lower-order partial derivatives. Repeated computations of lower-order derivatives cause multidual-based and multicomplex-based differentiation to become inefficient, especially for a large number of variables. Moreover, to build a complete nth-order Taylor series expansion of r variables using multidual or multicomplex algebra, d_p evaluations are required, see Eq. (13). An alternative algebra called OTI algebra forms the full nth-order Taylor series expansion without repeating derivative computations [16].

Similar to multidual numbers, OTI numbers contain r imaginary bases

ϵ_{1}, \dots, ϵ_{r}

where r is the number of variables of interest. Multidual numbers follow a fixed truncation of

ϵ_{k}^{2} = 0

⁠; however, in OTI numbers, the truncation condition is set to truncate those imaginary directions with order greater than n, where n is the highest order of derivative desired, that is,

ϵ_{k}^{p} = 0

where p > n. The set of OTI numbers with r bases and truncation order n is represented as

O T I_{r}^{n}

⁠. An example of an OTI number with truncation order n = 3 and number of bases r = 2,

a^{*} \in O T I_{2}^{3}

⁠, is

a^{*} = \underset{order 0}{\underset{︸}{a_{R}}} + \underset{order 1}{\underset{︸}{a_{ϵ_{1}} ϵ_{1} + a_{ϵ_{2}} ϵ_{2}}} + \underset{order 2}{\underset{︸}{a_{ϵ_{1}^{2}} ϵ_{1}^{2} + a_{ϵ_{1} ϵ_{2}} ϵ_{1} ϵ_{2} + a_{ϵ_{2}^{2}} ϵ_{2}^{2}}} + \underset{order 3}{\underset{︸}{a_{ϵ_{1}^{3}} ϵ_{1}^{3} + a_{ϵ_{1}^{2} ϵ_{2}} ϵ_{1}^{2} ϵ_{2} + a_{ϵ_{1} ϵ_{2}^{2}} ϵ_{1} ϵ_{2}^{2} + a_{ϵ_{2}^{3}} ϵ_{2}^{3}}}

(34)

where superscript “ $*$ ” indicates hypercomplex parts exist, subscript R denotes the real part of $a^{*}$ ⁠, and $a_{α}$ is the real coefficient along imaginary direction $α$ ⁠. The total number of terms in an OTI number, d, is the same as a Taylor series expansion truncated to order n, see Eq. (12).

Derivative computation with OTI numbers is as follows. Variables of interest are perturbed along independent basis directions. The truncation order n is set to the maximum order of derivative required. For example, for third-order partial derivatives with respect to two variables (r = 2, n = 3), the perturbations are

x^{*} = x_{0} + ϵ_{1}, x^{*} \in O T I_{2}^{3}

(35)

y^{*} = y_{0} + ϵ_{2}, y^{*} \in O T I_{2}^{3}

(36)

The function of interest is evaluated with the perturbed inputs,

f (x^{*}, y^{*})

⁠. The result is an OTI number with a total number of bases r and order n containing the partial derivatives of the function. Following the perturbations as in Eqs. and , the result is

\begin{matrix} f (x^{*}, y^{*}) = f (x_{0}, y_{0}) + \frac{1}{1!} D_{x}^{(1)} ϵ_{1} + \frac{1}{1!} D_{y}^{(1)} ϵ_{2} \\ + \frac{1}{2!} D_{x^{2}}^{(2)} ϵ_{1}^{2} + \frac{1}{1! 1!} D_{x y}^{(2)} ϵ_{1} ϵ_{2} + \frac{1}{2!} D_{y^{2}}^{(2)} ϵ_{2}^{2} + \frac{1}{3!} D_{x^{3}}^{(3)} ϵ_{1}^{3} \\ + \frac{1}{2! 1!} D_{x^{2} y}^{(3)} ϵ_{1}^{2} ϵ_{2} + \frac{1}{1! 2!} D_{x y^{2}}^{(3)} ϵ_{1} ϵ_{2}^{2} + \frac{1}{3!} D_{y^{3}}^{(3)} ϵ_{2}^{3} . \end{matrix}

(37)

Equation contains all derivatives up to third-order with respect to x and y. OTI numbers avoid computation of repeated imaginary directions using algebraic simplification. The coefficients of the form

1 / p!

⁠, referred to as scaling factors, are used in OTI numbers so that algebraic operations can be conducted without tracking the coefficient of each basis [16]. These scaling factors are identical to the coefficients that appear in the Taylor series expansion. Derivatives are retrieved by extracting the coefficient along the corresponding direction and multiplying the result by the corresponding scaling factor as

f (x_{0}, y_{0}) = R e [f (x^{*}, y^{*})]

(38)

\begin{matrix} D_{x}^{(1)} = \frac{\partial f}{\partial x} (x_{0}, y_{0}) = 1! \cdot I m_{ϵ_{1}} [f (x^{*}, y^{*})] \\ ⋮ \end{matrix}

(39)

D_{xyy}^{(3)} = \frac{\partial^{3} f}{\partial x \partial y^{2}} (x_{0}, y_{0}) = 1! 2! \cdot I m_{ϵ_{1} ϵ_{2}^{2}} [f (x^{*}, y^{*})]

(40)

D_{yyy}^{(3)} = \frac{\partial^{3} f}{\partial y^{3}} (x_{0}, y_{0}) = 3! \cdot I m_{ϵ_{2}^{3}} [f (x^{*}, y^{*})]

(41)

The pth-order derivative with respect to x is associated with the $ϵ_{1}^{p}$ imaginary direction as ϵ₁ was the imaginary direction used to perturb x. The coefficient along the $ϵ_{1}^{p}$ direction has a scaling factor $p!$ ⁠. A sparse implementation of OTI algebra is used so that only the nonzero imaginary coefficients are accounted for in the computations [16]. This reduces the computational cost and may lead to a decrease in memory usage.

2.3 Hypercomplex Finite Element Method for Transient Heat Transfer Problems.

Hypercomplex algebra has been implemented in the finite element formulation, known as ZFEM. This method was used to calculate partial derivatives for the numerical example in Sec. 3. The implementation of OTI-based HYPAD within the finite element formulation is described in the context of linear transient heat transfer in Sec. 2.3.1 and the solution to the OTI-based finite element system of equations is described in Sec. 2.3.2.

2.3.1 Hypercomplex Finite Element Method Formulation in Transient Heat Transfer Problems.

The following describes ZFEM for transient heat transfer problems in two dimensions with implicit Euler time integration. The differential equation for the transient heat transfer problem is defined as

C_{p} ρ \frac{\partial T}{\partial t} - \nabla \cdot (k \nabla T) = s in Ω

(42)

T (x, t = 0) = T_{0} in Ω

(43)

T (x, t) = \bar{T} on Γ_{1}

(44)

k \nabla T \cdot \hat{n} = h (T (x, t) - T_{\infty}) on Γ_{2}

(45)

where $T (x, t)$ is the temperature field valid for any point x in the domain Ω, s is the inner heat generation rate per unit volume, t is time, C_p is the specific heat, ρ is the material density, k is the thermal conductivity, T₀ is the initial temperature distribution, $\bar{T}$ is the prescribed temperature on boundary $Γ_{1}$ ⁠, h is the convective heat transfer coefficient occurring on boundary $Γ_{2}, T_{\infty}$ is the convection fluid temperature, and $\hat{n}$ is the outward unit normal of boundary $Γ_{2}$ [40].

Hypercomplex finite element method consists primarily in elevating the conventional operations to hypercomplex algebra in order to enable derivative computation. That is, input parameters are perturbed along imaginary directions as in Eqs. and which may correspond to material, geometrical, and/or boundary condition parameters. The temperature field

T (x, t)

becomes a hypercomplex temperature field and by nodal temperature solutions that are interpolated by the element basis functions as

T^{*} (x, t) = N u_{t}^{*}

(46)

where $u_{t}^{*}$ is the nodal temperature solution vector at time t and N is the matrix of real-valued element shape functions. In this work, the superscript “ $*$ ” of a matrix denotes that the matrix contains OTI numbers. The solution vector becomes hypercomplex because the input parameters are perturbed along imaginary directions that contain derivative information.

Additional degrees-of-freedom are created in order to account for the OTI imaginary directions. The number of integration points remains the same, however, these integrals now contain hypercomplex parts. Figure 1 shows the additional degrees-of-freedom required for an OTI analysis for r = 2 and n = 3 for a quadrilateral element with four integration points.

Fig. 1

View large Download slide

Additional degrees-of-freedom required for a hypercomplex four-noded quadrilateral element using OTI numbers with number of bases r = 2 and truncation order n = 3 (⁠ $O T I_{2}^{3}$ ⁠) in order to calculate first-order through third-order partial derivatives with respect to two variables

The nodal temperature is found by solving a linear system of hypercomplex equations as

K^{*} u_{t_{i + 1}}^{*} = f^{*}

(47)

where $K^{*}$ is the global system matrix and $f^{*}$ is the global source vector where both are formed in an element-by-element approach as in traditional finite element computations. The vector $u_{t_{i + 1}}^{*}$ represents the unknown nodal temperatures for time $t_{i + 1} = t_{i} + Δ t$ ⁠. The solution is known at time t_i and $Δ t$ is the time increment.

The elemental system matrix

K^{e}^{*}

for an implicit Euler integration scheme is formed at the elemental level, denoted by the superscript e, as

K^{e}^{*} = \frac{1}{Δ t} M^{e}^{*} + K_{c}^{e}^{*} + K_{h}^{e}^{*}

(48)

where contributions to the elemental system matrix are denoted by

M^{e}^{*}

for the capacitance,

K_{c}^{e}^{*}

for the conduction, and

K_{h}^{e}^{*}

for the convection. The matrices in Eq. are defined as

M^{e}^{*} = \int_{Ω^{e}} C_{p}^{*} ρ^{*} N^{T} N | J^{*} | d Ω

(49)

K_{c}^{e}^{*} = \int_{Ω^{e}} k^{*} B^{*}^{T} B^{*} | J^{*} | d Ω

(50)

K_{h}^{e}^{*} = \int_{Γ^{e}} h^{*} N^{T} N | J^{*} | d Γ

(51)

where

J^{*}

is the Jacobian matrix which is a function of the hypercomplex nodal coordinates [15],

| J^{*} |

is the determinant of the Jacobian matrix, and

B^{*}

is defined as

B^{*} = J^{* - 1} [\begin{matrix} \partial N_{1} / \partial ξ & \partial N_{2} / \partial ξ & \partial N_{3} / \partial ξ & \partial N_{4} / \partial ξ \\ \partial N_{1} / \partial η & \partial N_{2} / \partial η & \partial N_{3} / \partial η & \partial N_{4} / \partial η \end{matrix}]

(52)

The transpose of an OTI matrix is carried out in the same manner as the transpose of a real-valued matrix. Integrals in Eqs. (49)–(51) are solved using conventional 2 × 2 Gauss integration.

The elemental load vector

f^{e}^{*}

is defined as

f^{e}^{*} = \frac{1}{Δ t} M^{e}^{*} u_{t_{i}}^{e}^{*} + f_{s}^{e}^{*} + f_{h}^{e}^{*}

(53)

where

f_{s}^{e}^{*}

accounts for the heat source effects,

f_{h}^{e}^{*}

adds the convective contributions, and

u_{t_{i}}^{e}^{*}

is the known nodal temperature solution at time t_i. The load vector terms are computed in a similar manner by

f_{h}^{e}^{*} = \int_{Γ^{e}} h^{*} T_{\infty}^{*} N | J^{*} | d Γ

(54)

f_{s}^{e}^{*} = \int_{Ω^{e}} s^{*} N | J^{*} | d Ω

(55)

The global system of equations is found by assembling the elemental matrices into the global system of equations.

2.3.2 Solution to the Order Truncated Imaginary System of Equations.

The OTI global system of equations contains many degrees-of-freedom due to the imaginary parts. Solving the full system of equations becomes computationally demanding in terms of memory storage. Therefore, the block forward substitution method, which is based on the Cauchy–Riemann matrix form of the OTI algebra, was used to solve the OTI system of equations in a computationally efficient manner [20]. The block forward substitution method only needs to be implemented once and is independent of the time discretization used. A comprehensive method to form the OTI matrix form is presented in Ref. [16]. As an example, the OTI system of equations for r = 2 and n = 3 is

[\begin{matrix} K_{R} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ K_{ϵ_{1}} & K_{R} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ K_{ϵ_{2}} & 0 & K_{R} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ K_{ϵ_{1}^{2}} & K_{ϵ_{1}} & 0 & K_{R} & 0 & 0 & 0 & 0 & 0 & 0 \\ K_{ϵ_{1} ϵ_{2}} & K_{ϵ_{2}} & K_{ϵ_{1}} & 0 & K_{R} & 0 & 0 & 0 & 0 & 0 \\ K_{ϵ_{2}^{2}} & 0 & K_{ϵ_{2}} & 0 & 0 & K_{R} & 0 & 0 & 0 & 0 \\ K_{ϵ_{1}^{3}} & K_{ϵ_{1}^{2}} & 0 & K_{ϵ_{1}} & 0 & 0 & K_{R} & 0 & 0 & 0 \\ K_{ϵ_{1}^{2} ϵ_{2}} & K_{ϵ_{1} ϵ_{2}} & K_{ϵ_{1}^{2}} & K_{ϵ_{2}} & K_{ϵ_{1}} & 0 & 0 & K_{R} & 0 & 0 \\ K_{ϵ_{1} ϵ_{2}^{2}} & K_{ϵ_{2}^{2}} & K_{ϵ_{1} ϵ_{2}} & 0 & K_{ϵ_{2}} & K_{ϵ_{1}} & 0 & 0 & K_{R} & 0 \\ K_{ϵ_{2}^{3}} & 0 & K_{ϵ_{2}^{2}} & 0 & 0 & K_{ϵ_{2}} & 0 & 0 & 0 & K_{R} \end{matrix}] {\begin{matrix} u_{R} \\ u_{ϵ_{1}} \\ u_{ϵ_{2}} \\ u_{ϵ_{1}^{2}} \\ u_{ϵ_{1} ϵ_{2}} \\ u_{ϵ_{2}^{2}} \\ u_{ϵ_{1}^{3}} \\ u_{ϵ_{1}^{2} ϵ_{2}} \\ u_{ϵ_{1} ϵ_{2}^{2}} \\ u_{ϵ_{2}^{3}} \end{matrix}} = {\begin{matrix} f_{R} \\ f_{ϵ_{1}} \\ f_{ϵ_{2}} \\ f_{ϵ_{1}^{2}} \\ f_{ϵ_{1} ϵ_{2}} \\ f_{ϵ_{2}^{2}} \\ f_{ϵ_{1}^{3}} \\ f_{ϵ_{1}^{2} ϵ_{2}} \\ f_{ϵ_{1} ϵ_{2}^{2}} \\ f_{ϵ_{2}^{3}} \end{matrix}}

(56)

The new system of equations can be solved by blocks in a forward substitution scheme due to the lower-triangular shape of OTI matrix form. First, the real-only system of equations is solved as

K_{R} u_{R} = f_{R}

(57)

This system of equations corresponds to the traditional finite element system of equations and is equivalent to the conventional finite element problem. Following the solution of the real part of the system of equations, coefficients of order 1 are solved for the first-order derivatives as

K_{R} u_{ϵ_{1}} = f_{ϵ_{1}} - K_{ϵ_{1}} u_{R}

(58)

K_{R} u_{ϵ_{2}} = f_{ϵ_{2}} - K_{ϵ_{2}} u_{R}

(59)

where both systems depend only on the real part of the solution vector

u_{R}

and the first-order terms of

K^{*}

and

f^{*}

⁠. They can be solved simultaneously since both first-order systems are independent of each other and share the same system matrix

K_{R}

⁠. This is accomplished by stacking all vector terms into a single matrix

u^{(1)}

and the corresponding right hand sides into a single matrix as

K_{R} [\begin{matrix} u_{ϵ_{1}} & u_{ϵ_{2}} \end{matrix}] = [\begin{matrix} f_{ϵ_{1}} - K_{ϵ_{1}} u_{R} & f_{ϵ_{2}} - K_{ϵ_{2}} u_{R} \end{matrix}]

(60)

The system for the terms of order 2, for the second-order derivatives, arising from the block forward substitution is

K_{R} u_{ϵ_{1}^{2}} = f_{ϵ_{1}^{2}} - K_{ϵ_{1}} u_{ϵ_{1}} - K_{ϵ_{1}^{2}} u_{R}

(61)

K_{R} u_{ϵ_{1} ϵ_{2}} = f_{ϵ_{1} ϵ_{2}} - K_{ϵ_{1}} u_{ϵ_{2}} - K_{ϵ_{2}} u_{ϵ_{1}} - K_{ϵ_{2}^{2}} u_{R}

(62)

K_{R} u_{ϵ_{2}^{2}} = f_{ϵ_{2}^{2}} - K_{ϵ_{2}} u_{ϵ_{2}} - K_{ϵ_{2}^{2}} u_{R}

(63)

which can also be solved simultaneously after solving first-order terms, as they are independent. The third-order terms are solved as

K_{R} u_{ϵ_{1}^{3}} = f_{ϵ_{1}^{3}} - K_{ϵ_{1}} u_{ϵ_{1}^{2}} - K_{ϵ_{1}^{2}} u_{ϵ_{2}} - K_{ϵ_{1}^{3}} u_{R}

(64)

K_{R} u_{ϵ_{1}^{2} ϵ_{2}} = f_{ϵ_{1}^{2} ϵ_{2}} - K_{ϵ_{1}} u_{ϵ_{1} ϵ_{2}} - K_{ϵ_{2}} u_{ϵ_{1}^{2}} - K_{ϵ_{1}^{2}} u_{ϵ_{2}} - K_{ϵ_{1} ϵ_{2}} u_{ϵ_{1}} - K_{ϵ_{1}^{2} ϵ_{2}} u_{R}

(65)

K_{R} u_{ϵ_{1} ϵ_{2}^{2}} = f_{ϵ_{1} ϵ_{2}^{2}} - K_{ϵ_{1}} u_{ϵ_{2}^{2}} - K_{ϵ_{2}} u_{ϵ_{1} ϵ_{2}} - K_{ϵ_{1} ϵ_{2}} u_{ϵ_{2}} - K_{ϵ_{2}^{2}} u_{ϵ_{1}} - K_{ϵ_{1} ϵ_{2}^{2}} u_{R}

(66)

K_{R} u_{ϵ_{2}^{3}} = f_{ϵ_{2}^{3}} - K_{ϵ_{2}} u_{ϵ_{2}^{2}} - K_{ϵ_{2}^{2}} u_{ϵ_{2}} - K_{ϵ_{2}^{3}} u_{R}

(67)

which can also be solved simultaneously after solving for the second-order terms. In general, the solution to the order p terms is accomplished by

K_{R} [u_{α_{1}^{p}}, \dots, u_{α_{i}^{p}}, \dots] = [RH S_{α_{1}^{p}}, \dots, RH S_{α_{i}^{p}}, \dots]

(68)

where the solution to the coefficient of the $α_{i}^{p}$ imaginary direction is found by forming the right hand side term $RH S_{α_{i}^{p}}$ that depends on the $α_{i}^{p}$ term of $f^{*}$ and all terms from $K^{*}$ and $u^{*}$ such that when multiplied, result in $α_{i}^{p}$ direction. Since all systems depend on $K_{R}$ ⁠, Lower-upper (LU) decomposition or Cholesky factorization can be used to reduce the computational effort in solving the hypercomplex system of equations.

Hypercomplex automatic differentiation can be implemented in nonlinear finite element simulations as well. This topic is not covered in this work; however, the interested reader can see examples of nonlinear HYPAD implementations in Refs. [20], [41], and [42].

3 Numerical Example: Uncertainty Quantification of Transient Heat Transfer in a Fin

HYPAD-UQ was applied to a transient linear thermal finite element analysis of a heated fin. The stochastic perturbation finite element method was also conducted to compare accuracy with HYPAD-UQ for the same computational time. The analytical solution of the fin temperature in Sec. 3.1 and a description of the FEA model are presented in Sec. 3.2. The computational time and the error of partial derivatives calculated with HYPAD are presented in Sec. 3.3. Numerical results are presented in Sec. 3.4 for two cases of random variable parameters.

3.1 Analytical Solution.

A fin is initially at ambient temperature

T_{\infty}

and is heated by a wall of constant temperature T_W, see Fig. 2. The length b of the fin is long enough to assume negligible heat loss through the tip of the fin. In addition, the thickness δ is small enough so that heat only flows along the x axis. The transient normalized temperature distribution θ along the x axis is

θ (X, τ) = \frac{T (X, τ) - T_{\infty}}{T_{W} - T_{\infty}} = θ_{S S} (X) + θ_{τ} (X, τ)

(69)

where

T (X, τ)

is the temperature distribution of the fin that depends on a length parameter

X = x / b

and a time parameter

τ = t k / b^{2} ρ C_{p}

⁠, t is time, k is thermal conductivity, ρ is density, and C_p is specific heat [43]. The normalized steady-state temperature distribution is expressed as

θ_{S S} (X) = \frac{\cosh (ω X)}{\cosh (ω)}

(70)

where

ω^{2} = 2 h_{U} b^{2} / k δ L

⁠, L is plate depth, and h_U is the convection heat transfer coefficient. The normalized transient temperature distribution is expressed as

θ_{τ} (X, τ) = 2 \sum_{j = 1}^{\infty} {(- 1)}^{j + 1} (\frac{θ_{0}}{λ_{j}} - \frac{λ_{j}}{ω^{2} + λ_{j}^{2}}) cos (λ_{j} X) e^{- (ω^{2} + λ_{j}^{2}) τ}

(71)

where $λ_{j} = π (2 j - 1) / 2$ and $θ_{0} = (T_{0} - T_{\infty}) / (T_{W} - T_{\infty})$ ⁠. This analysis assumed the initial temperature $T_{0} = T_{\infty}$ ⁠, therefore, $θ_{0} = 0$ ⁠. In addition, the infinite sum in Eq. (71) was solved numerically using a truncation of j = 100.

Fig. 2

View large Download slide

Diagram of a heated fin

3.2 Finite Element Model.

The fin was modeled in 2D with 10,000 four-noded quadrilateral elements and 11,011 nodes using an in-house Python finite element code. The analysis was conducted using implicit Euler time integration for a total of 450 s in order to reach steady-state. A time-step of 1 s was used which was determined from a convergence study [44]. A real-variable version of the FEA simulation was used to generate samples to train the linear regression models. A hypercomplex version of the finite element code was used to calculate partial derivatives of temperature evaluated at the mean using OTI-based HYPAD, as described in Secs. 2.3.1 and 2.3.2.

The 2D FEA formulation does not account for out-of-plane heat dissipation. Therefore, a volumetric sink term U_V was implemented in Eq. . The relation between U_V and s is defined as

s = - U_{V} h_{U} (T - T_{\infty})

(72)

A dependence on h_U was included in Eq. (72) to ensure derivatives with respect to h_U account for heat transfer in the out-of-plane direction. U_V was approximated by minimizing the difference of T between the FEA simulation and the analytic solution, both modeled at the mean values of random variables, see Table 1. The resulting value was $U_{V} \approx$ 735/114 $\approx$ 6.447 $m^{- 1}$ [44].

Table 1

Random variable parameters for the two case studies

			Case 1		Case 2
Random variable	Symbol	Mean	Distribution	δ_x	Distribution	δ_x
Thermal conductivity	k	7.1 W/m K	Normal	0.10	Log-normal	0.20
Specific heat	C_p	580 J/kg K	Normal	0.03	Log-normal	0.20
Density	ρ	4430 kg/m³	Normal	0.03	Log-normal	0.20
Heat transfer coefficient	h_U	114 W/m² K	Normal	0.10	Log-normal	0.20
Ambient temperature	$T_{\infty}$	283 $K$	Normal	0.001	Triangle	0.01
Temperature of heat source	T_W	389 $K$	Normal	0.05	Uniform	0.20
Length of fin	b	51 × 10⁻³ $m$	Normal	0.01	Uniform	0.20

			Case 1		Case 2
Random variable	Symbol	Mean	Distribution	δ_x	Distribution	δ_x
Thermal conductivity	k	7.1 W/m K	Normal	0.10	Log-normal	0.20
Specific heat	C_p	580 J/kg K	Normal	0.03	Log-normal	0.20
Density	ρ	4430 kg/m³	Normal	0.03	Log-normal	0.20
Heat transfer coefficient	h_U	114 W/m² K	Normal	0.10	Log-normal	0.20
Ambient temperature	$T_{\infty}$	283 $K$	Normal	0.001	Triangle	0.01
Temperature of heat source	T_W	389 $K$	Normal	0.05	Uniform	0.20
Length of fin	b	51 × 10⁻³ $m$	Normal	0.01	Uniform	0.20

3.3 High-Order Order Truncated Imaginary-Based Derivatives.

Hypercomplex automatic differentiation was used to compute all first-order through seventh-order partial derivatives of the nodal temperature at the tip of the fin with respect to the seven random variables evaluated at the mean, see Table 1, at each time-step $Δ t = 1$ s for the entire analysis $1 \leq t \leq 450$ s. Although this implementation of OTI differentiation requires derivatives to be calculated at every node in the FEA model for each time-step, this study only analyzed temperature at the tip of the fin; however, UQ results for each node could have been calculated from theses derivatives. For example, variance estimates could have been computed for each node in the model. At each node, a total of 3431 derivatives were computed in addition to the conventional FEA solution. In particular, a total of seven first-order, 28 second-order, 84 third-order, 210 fourth-order, 462 fifth-order, 924 sixth-order, and 1716 seventh-order partial derivatives were calculated. Of these derivatives, seven first-order, 25 second-order, 65 third-order, 140 fourth-order, 266 fifth-order, 462 sixth-order, and 750 seventh-order derivatives were nonzero. Due to the use of the sparse OTI implementation, only these nonzero imaginary coefficients were accounted for within the algebraic operations in the finite element procedure. In summary, a total of 18,894,876 nonzero first-order through seventh-order derivatives were calculated and stored in the model for each time-step but were overwritten during the next time-step.

The real solution and the HYPAD partial derivatives are in good agreement with the analytical solution, see Fig. 3 for the normalized root mean square error (NRMSE) of HYPAD-based derivatives of the FEA model compared to derivatives of the analytical solution, which were computed numerically with HYPAD. The NRMSE measures the error across all time steps and is defined as

NRMSE = \frac{1}{max (ϕ_{analytic})} \sqrt{\frac{\sum_{i = 1}^{N} {(ϕ_{approx .}^{(i)} - ϕ_{analytic}^{(i)})}^{2}}{M}}

(73)

where $ϕ_{approx .}^{(i)}$ and $ϕ_{analytic}^{(i)}$ are the HYPAD-based derivatives of the FEA model and analytical derivative values, respectively, corresponding to time-step i, M is the total number of time steps, and $max (ϕ_{analytic})$ is the maximum of all analytical values. The conventional FEA solution for T has a NRMSE error below 1 × 10⁻⁴. However, the HYPAD derivatives have a slightly larger error which increases with the order of derivative.

Fig. 3

View large Download slide

Box plot of the NRMSE of T calculated with the real-variable FEA simulation (order = 0) and first-order through seventh-order OTI-based HYPAD partial derivatives of T

Figure 3 shows that the error in derivatives increase as the order of derivative increases. This is due to the numerical discretization error that accrues with each order since higher-order derivatives depend on values of lower-order derivatives. Furthermore, there is an additional source of error in the FEA simulation due to the numerical approximation in the out-of-plane heat dissipation, see Sec. 3.2.

A real-variable analysis and seven independent OTI-based FEA simulations were performed to analyze the computational performance, see Fig. 4. Each simulation calculates the traditional real-valued temperature and all combinations of first-order through nth-order derivatives. The real solution and all combinations of first-order through nth-order partial derivatives were calculated in a single OTI run. Higher-order HYPAD derivatives require more CPU time; however, the number of partial derivatives being calculated also increases. The time to compute the elemental matrices and assemble the global system of equations are denoted by the assembly time. Solution time refers to the CPU time to solve the system of equations.

Fig. 4

View large Download slide

CPU time to compute first-order through nth-order derivatives with OTI-based HYPAD relative to the real solution CPU time

3.4 Uncertainty Quantification of Temperature at Tip of Fin.

Partial derivatives calculated with HYPAD were used in the Taylor series to calculate central moments, probability distributions, and Sobol' indices of temperature T at the tip of a heated fin. The accuracy of HYPAD-UQ was compared to random sampling of the analytical solution following a latin hypercube sampling (LHS) design. Computational time was compared to a Taylor series expansion constructed by regression, similar to the method described in Ref. [21]. The unknown coefficients of the complete polynomial basis of the Taylor series were solved for the regression analysis using ordinary least squares (OLSs). An additional comparison was made to the Taylor series created from partial derivatives of the analytical solution, also calculated with OTI-based HYPAD, to show the effect of the error in the HYPAD derivatives on the UQ metrics. Herein, derivatives calculated from the analytical solution are referred to as the analytical derivatives while the derivatives calculated from the FEA model will be referred to as HYPAD derivatives.

Uncertainty in T was quantified for two cases of parameters of the random variables. The first case, presented in Sec. 3.4.1, assumed all variables follow a Gaussian distribution with low coefficient of variations, $δ_{x} = σ_{x} / μ_{x}$ where σ_x is the standard deviation of random variable x. The coefficient of variations δ_x for the material properties k, C_p, and ρ were assumed to be the experimental uncertainty values listed in Ref. [45] for Ti6-Al4V and δ_x was assumed for the remaining variables, see Table 1. The mean values were assumed to be the values listed in Ref. [43]. The second case, presented in Sec. 3.4.2, assumed non-Gaussian distributions with a larger δ_x for each variable. The mean values remained unchanged, therefore, the same Taylor series expansion was used for both cases. Standard deviations were set to 20% of the mean values except $T_{\infty}$ ⁠. The coefficient of variation of $T_{\infty}$ was chosen so that other random variables significantly contribute to the variance of temperature at the tip of the fin, see Table 1. The thickness δ and depth L of the fin were assumed to be deterministic values of $δ = 4.75$ mm and L = 1 m, respectively, for both cases [43].

3.4.1 Gaussian Random Variables With Low Coefficient of Variations.

HYPAD-UQ was applied to quantify the uncertainty of temperature at the tip of heated fin following the parameters listed in Table 1 under case 1. The HYPAD-UQ approximations for the first four central moments are in close agreement with those estimated from 1 × 10⁶ evaluations of the analytical solution following an LHS design, see Fig. 5. HYPAD-UQ central moments were calculated algebraically using Eq. (26) for the first-order through third-order Taylor series expansions, as described in Sec. 2.1. Moments approximated with fourth-order and higher-order expansions were calculated from 1 × 10⁶ samples of the Taylor series following an LHS design. The computational cost of the seventh-order Taylor series using HYPAD-based derivatives was about 206 times a single traditional FEA run, see Fig. 4. Therefore, moments calculated with a third-degree OLS model trained with m = 206 samples was compared to the HYPAD-UQ moments, shown in Fig. 5. Central moments were also calculated algebraically from the third-degree OLS model using Eq. (26). The 95% confidence intervals of the OLS-based central moments and analytical central moments were calculated from the 2.5 and 97.5 percentiles produced from bootstrap resampling with 100 replications [46].

Fig. 5

View large Download slide

Case 1—central moments of T (K) calculated by LHS of the analytical solution, HYPAD-UQ, and a third-degree OLS model trained with m = 206 samples. Shaded regions on the analytical and OLS central moments represent the 95% confidence interval: (a) expected value, (b) variance, (c) skewness, and (d) kurtosis.

Figure 6 shows the NRMSE of central moments approximated by HYPAD-UQ and OLS models against the CPU time relative to a single real-variable analysis. HYPAD-UQ showed a significant computational advantage over the OLS models. The NRMSE was calculated with Eq. (73) assuming $ϕ_{analytic}$ to be the central moments approximated from 1 × 10⁶ evaluations of the analytical solution. To compare accuracy based on CPU time, the number of training samples of the OLS models were chosen to match the CPU time required to calculate HYPAD derivatives, see Fig. 4. The first-degree through third-degree OLS models were trained with 11, 22, 65, 134, and 206 samples of the real-variable FEA simulation following an LHS design. OLS models were only trained for cases where the number of samples exceeded the number of unknown coefficients.

Fig. 6

View large Download slide

Case 1—NRMSE of central moments approximated by HYPAD-UQ, the Taylor series constructed with analytical derivatives, and OLS models plotted against CPU time relative to a single real-valued FEA run. Error bars on the OLS NRMSE values represent the 95% confidence interval: (a) NRMSE of the expected value, (b) NRMSE of variance, (c) NRMSE of skewness, and (d) NRMSE of kurtosis.

Figure 6 also shows the NRMSE of moments calculated from the Taylor series constructed from partial derivatives of the analytical solution. This shows a comparison of HYPAD-UQ and a numerically accurate Taylor series using the analytical solution. The error in HYPAD derivatives limits the accuracy of the central moment approximations. However, HYPAD-UQ is still in very good agreement with random sampling of the analytic solution, see Fig. 5.

The probability density function (PDF) generated from sampling the HYPAD-based Taylor series expansion is in close agreement with the PDF estimated from 1 × 10⁶ evaluations of the analytical solution and the third-degree OLS model trained with m = 206 samples, see Fig. 7. Probability densities were estimated using 1 × 10⁶ evaluations of each model with the Epanechnikov method for nonparametric kernel regression with a window size of about 0.076 for each model [47].

Fig. 7

View large Download slide

Case 1—PDF of T (K) at t = 450 s calculated from LHS of the analytical solution, HYPAD-UQ, and a third-degree OLS model trained with m = 206 samples

The accuracy in the tails of the distributions approximated with HYPAD-UQ converge to the distribution approximated by the analytical solution, see Fig. 8 for evaluations of the HYPAD-based Taylor series and the analytical solution plotted on a Gaussian probability scale at a steady-state time of t = 450 s. The first-order Taylor series expansion produces a Gaussian distribution due to all random variables being Gaussian. However, higher-order Taylor series expansions converge to the non-Gaussian distribution of the analytical solution of T.

Fig. 8

View large Download slide

Case 1—Cumulative probability of T at t = 450 s calculated from evaluations of the analytical solution and HYPAD-UQ plotted on a Gaussian probability scale

Sobol' indices calculated with HYPAD-UQ are in good agreement with sampling of the analytical solution following an LHS design, see Fig. 9 for the values of Sobol' indices through time. HYPAD-UQ Sobol' indices were calculated algebraically as described in Sec. 2.2. The sum of all interaction effects was calculated by subtracting the sum of all main effects from one, at each time-step. Sobol' indices were calculated from the analytical solution by the sampling method described in Ref. [31], with a total of 1 × 10⁸ samples for each variable at each time-step.

Fig. 9

View large Download slide

Case 1—Cumulative values of main and interaction Sobol' indices of T calculated by LHS of the analytical solution and HYPAD-UQ: (a) LHS of the analytic solution, (b) first-order HYPAD-UQ, (c) second-order HYPAD-UQ, and (d) third-order HYPAD-UQ

Table 2 shows the maximum absolute error ε of Sobol' indices calculated by each method, measured by

ε = max [| ϕ_{approx .}^{(i)} - ϕ_{analytic}^{(i)} |]

(74)

where $ϕ_{approx .}^{(i)}$ is the Sobol' index approximated by HYPAD-UQ and $ϕ_{analytic}^{(i)}$ is the Sobol' index calculated from sampling of the analytic solution at each time point i. The first-order HYPAD-UQ approximations of Sobol' indices are within 3.1% maximum absolute error compared to LHS of the analytic solution across all time points, even though there is no approximation for any interaction effects. Higher-order HYPAD-based Taylor series expansions improve the approximations to be within 1.7% absolute error for second-order and 0.92% absolute error for third-order.

Table 2

Case 1—maximum absolute error ε of Sobol' indices across all time points between LHS of the analytical solution and HYPAD-UQ, the analytic Taylor series expansion, and OLS models trained with m = 206 samples

		Maximum absolute error ε (%)
Method	Order/degree	$S_{k}$	$S_{C_{p}}$	$S_{ρ}$	$S_{h_{U}}$	$S_{T_{\infty}}$	$S_{T_{w}}$	S_b	Interaction
HYPAD-UQ	First-order	2.8	0.15	0.083	1.5	3.1^a	1.1	0.12	3.1^a
	Second	1.4	0.083	0.078	1.7^a	1.5	0.88	0.065	0.43
	Third	0.43	0.11	0.030	0.92^a	0.43	0.42	0.043	0.42
Analytic derivatives	First-order	2.5	0.099	0.074	1.8	2.6	1.6	0.15	3.1^a
	Second	1.2	0.047	0.11	2.2^a	1.0	0.64	0.094	0.41
	Third	0.92	0.071	0.052	1.2^a	1.0	0.88	0.074	0.41
OLS, m = 206	First-degree	1.8	0.78	0.50	2.6	1.6	4.9^a	0.32	3.1
	Second	0.80	0.12	0.073	2.0^a	1.1	1.7	0.15	0.40
	Third	0.80	0.074	0.055	1.9^a	1.0	1.1	0.074	0.40

		Maximum absolute error ε (%)
Method	Order/degree	$S_{k}$	$S_{C_{p}}$	$S_{ρ}$	$S_{h_{U}}$	$S_{T_{\infty}}$	$S_{T_{w}}$	S_b	Interaction
HYPAD-UQ	First-order	2.8	0.15	0.083	1.5	3.1^a	1.1	0.12	3.1^a
	Second	1.4	0.083	0.078	1.7^a	1.5	0.88	0.065	0.43
	Third	0.43	0.11	0.030	0.92^a	0.43	0.42	0.043	0.42
Analytic derivatives	First-order	2.5	0.099	0.074	1.8	2.6	1.6	0.15	3.1^a
	Second	1.2	0.047	0.11	2.2^a	1.0	0.64	0.094	0.41
	Third	0.92	0.071	0.052	1.2^a	1.0	0.88	0.074	0.41
OLS, m = 206	First-degree	1.8	0.78	0.50	2.6	1.6	4.9^a	0.32	3.1
	Second	0.80	0.12	0.073	2.0^a	1.1	1.7	0.15	0.40
	Third	0.80	0.074	0.055	1.9^a	1.0	1.1	0.074	0.40

a

Largest ε for each method.

Both of the Taylor series expansions using HYPAD derivatives and analytic derivatives have similar accuracy. This shows that the error in the HYPAD derivatives do not significantly affect the accuracy of Sobol' indices. The OLS model trained with m = 206 samples also achieves high accuracy even with a first-degree polynomial basis. The lack of interaction effects is one reason for high accuracy of first-order and first-degree approximations. Section 3.4.2 repeats this analysis for random variable parameters that produce a substantial interaction effect.

3.4.2 Non-Gaussian Random Variables With Large Coefficient of Variations.

The analysis in Sec. 3.4.1 was repeated for independent random variables with larger coefficient of variations and non-Gaussian distributions, see Table 1 under case 2. The derivatives from Sec. 3.4.1 were reused in this analysis because the mean values of the random variables did not change.

Central moments approximated with HYPAD-UQ converge to the LHS results of the analytical solution using 1 × 10⁷ samples, see Fig. 10. However, fifth-order and higher-order approximations for skewness and kurtosis diverge from the analytic solution in the transient region between t = 20 s and t = 100 s.

Fig. 10

View large Download slide

Case 2—central moments of T (K) calculated by LHS of the analytical solution, HYPAD-UQ, and a second-degree OLS model trained with m = 206 samples. Shaded regions on the analytical and OLS central moments represent the 95% confidence interval: (a) expected value, (b) variance, (c) skewness, and (d) kurtosis.

HYPAD-UQ central moments were compared to a second-degree OLS model trained with m = 206 samples, which has approximately the same computational cost as the seventh-order HYPAD-based Taylor series expansion. Third-degree OLS results are not shown due to large confidence intervals in the moments which is a consequence of the limited number of training points. The fourth-order HYPAD-UQ expansion yielded the best results for the central moments, see Fig. 10, which took about 22 times a single real-valued analysis. These results achieved greater accuracy than the OLS model which required 206 times a single real-variable analysis.

HYPAD-UQ central moments achieved low NRMSE values; however, higher-order expansions do not guarantee a monotonic decrease in error, as shown by Fig. 11. This behavior does not solely come from the error in the derivatives as shown by the nonmonotonic convergence of the NRMSE of central moments calculated using the Taylor series expansion constructed using analytic derivatives.

Fig. 11

View large Download slide

Case 2—NRMSE of central moments approximated by HYPAD-UQ, the Taylor series constructed with analytical derivatives, and OLS models plotted against CPU time relative to a single real-valued FEA run. Error bars on the OLS NRMSE values represent the 95% confidence interval: (a) NRMSE of the expected value, (b) NRMSE of variance, (c) NRMSE of skewness, and (d)NRMSE of kurtosis.

Probability densities were estimated by the analytical solution, HYPAD-UQ, and a second-degree OLS, trained with m = 206 samples, using 1 × 10⁶ evaluations of each model following an LHS design, see Fig. 12. The Epanechnikov method for nonparametric kernel regression with a window size of about 0.33 was used to approximate the PDF of each model. The PDF approximated with HYPAD-UQ converges to the analytical solution at fourth-order. The PDF estimated from the OLS model is in close agreement with LHS of the analytic solution, however, third-order through seventh-order HYPAD-UQ methods are in closer agreement.

Fig. 12

View large Download slide

Case 2—PDF of T at t = 450 s calculated from LHS of the analytical solution, HYPAD-UQ, and a third-degree OLS model trained with m = 206 samples

Although the HYPAD-UQ approximations of the PDF converge to the analytical solution as the order of expansion increases, higher-order Taylor series expansions deviate from the analytical solution near low and high cumulative probabilities, see Fig. 13. An empirical cumulative distribution was created using 1 × 10⁶ evaluations of the analytical solution and first-order through seventh-order HYPAD-based Taylor series expansions. Figure 13(a) shows the cumulative distribution function (CDF) at a transient time of t = 35 s, which is the time point of maximum error in the seventh-order HYPAD-UQ approximation for kurtosis and Fig. 13(b) shows the CDF at a steady-state time of t = 450 s.

Fig. 13

View large Download slide

Case 2—Cumulative probability of T at t = 35 s and t = 450 s calculated from evaluations of the analytical solution and HYPAD-UQ plotted on a Gaussian probability scale: (a) cumulative probability at t = 35 s and (b) cumulative probability at t = 450 s

Higher-order Taylor series expansions can diverge from the solution, see Fig. 14. The same evaluations of the HYPAD-based Taylor series expansions used to generate the CDF, see Fig. 13, were compared to the analytical solution at t = 35 s, see Fig. 14(a), and t = 450 s, see Fig. 14(b). A majority of the Taylor series evaluations tend to approach the analytical solution as the expansion order increases; however, large errors in higher-order Taylor series predictions can occur. Figure 14(a) shows that at t = 450 s, a few combinations of random variable parameters cause the third-order, fifth-order, and seventh-order HYPAD-based Taylor series to significantly under predict T and the second-order, fourth-order, and sixth-order HYPAD Taylor series to significantly over predict T. The low values predicted by the fifth-order and seventh-order Taylor series expansions at t = 450 s cause the empirically generated CDF to deviate from the analytical solution near low cumulative probabilities, see Fig. 13(b). The even-ordered Taylor series expansions appear to yield a better CDF because these outliers get grouped with data of similar magnitude. A similar conclusion can be drawn for the CDF at t = 35 s.

Fig. 14

View large Download slide

Case 2—Values of T predicted by HYPAD-based Taylor series expansions and the analytical solution at t = 35 s and t = 450 s: (a) values of T predicted by HYPAD-based Taylor series expansions and the analytical solution at t = 35 s and (b) values of T predicted by HYPAD-based Taylor series expansions and the analytical solution at t = 450 s

Sobol' indices approximated with second-order and third-order HYPAD-UQ methods agree well with sampling of the analytic solution, see Fig. 15. Sobol' indices were calculated from the analytic solution the pick-and-freeze method, following Ref. [48]. This method required a total of $(r + 1) m$ samples to calculate all main effect Sobol' indices, where r is the number of random variables and $m =$ 1 × 10⁷ was chosen, following an LHS design.

Fig. 15

View large Download slide

Case 2—Cumulative values of main and interaction Sobol' indices of T calculated by LHS of the analytical solution and HYPAD-UQ: (a) LHS of the analytic solution, (b) first-order HYPAD-UQ, (c) second-order HYPAD-UQ, and (d) third-order HYPAD-UQ

The cumulative interaction effect is large with a maximum of about 27.8% of the total variance. Figure 16 shows the main and interaction effect Sobol' indices at a steady-state time of t = 450 s. At this time point, second-order HYPAD-UQ Sobol' indices are in good agreement with analytical Sobol' indices even though second-order HYPAD-UQ underestimates the total variance, see Fig. 10. Sobol' indices were omitted from Fig. 16 for values less than 0.02. HYPAD-UQ interaction effects were calculated algebraically using Eq. (21) for second-order and Eqs. (A3) and (A4) for third-order expansions. Interaction Sobol' indices were calculated from the analytical solution using the nested-loop sampling method in Ref. [31] with 1 × 10⁴ samples in each loop for a total of 1 × 10⁸ samples.

Fig. 16

View large Download slide

Case 2—Main and two-variable interaction Sobol' indices at t = 450 s calculated by HYPAD-UQ and LHS of the analytical solution

As shown in Fig. 16, main effects approximated by HYPAD-UQ improve with the order of the expansion. The third-order HYPAD-UQ Sobol' indices have a maximum of 13% absolute error across all time points, see Table 3 for the maximum absolute error of Sobol' indices using Eq. (74). Sobol' indices approximated with analytic derivatives are also within 13% maximum absolute error at third-order. Table 3 also shows that the first-degree OLS model trained with m = 206 samples approximates main effect Sobol' indices better than the first-order HYPAD-UQ method. However, second-order and higher-order HYPAD-UQ methods perform better than the first-degree and second-degree OLS models.

Table 3

Case 2—maximum absolute error ε of Sobol' indices across all time points between sampling of the analytical solution and HYPAD-UQ, the analytic Taylor series expansion, and OLS models trained with m = 206 samples

		Maximum absolute error ε (%)
Method	Order/degree	$S_{k}$	$S_{C_{p}}$	$S_{ρ}$	$S_{h_{U}}$	$S_{T_{\infty}}$	$S_{T_{w}}$	S_b	Interaction
HYPAD-UQ	First-order	3.2	1.1	1.7	3.4	46^a	7.6	12	28
	Second	2.9	1.0	1.7	0.86	34^a	9.1	11	13
	Third	1.2	0.86	0.69	2.1	13^a	1.9	2.9	9.5
Analytic derivatives	First-order	3.2	1.1	1.7	3.6	46^a	7.7	12	28
	Second	2.9	1.0	1.7	0.97	34^a	9.1	11	13
	Third	1.3	1.0	0.73	2.2	13^a	2.1	2.8	9.3
OLS, m = 206	First-degree	1.9	2.0	1.2	3.1	26	3	9.7	28^a
OLS, m = 206	Second	3.1	6.5	6.6	5.6	49	50^a	25	23

		Maximum absolute error ε (%)
Method	Order/degree	$S_{k}$	$S_{C_{p}}$	$S_{ρ}$	$S_{h_{U}}$	$S_{T_{\infty}}$	$S_{T_{w}}$	S_b	Interaction
HYPAD-UQ	First-order	3.2	1.1	1.7	3.4	46^a	7.6	12	28
	Second	2.9	1.0	1.7	0.86	34^a	9.1	11	13
	Third	1.2	0.86	0.69	2.1	13^a	1.9	2.9	9.5
Analytic derivatives	First-order	3.2	1.1	1.7	3.6	46^a	7.7	12	28
	Second	2.9	1.0	1.7	0.97	34^a	9.1	11	13
	Third	1.3	1.0	0.73	2.2	13^a	2.1	2.8	9.3
OLS, m = 206	First-degree	1.9	2.0	1.2	3.1	26	3	9.7	28^a
OLS, m = 206	Second	3.1	6.5	6.6	5.6	49	50^a	25	23

a

Largest ε of each method.

4 Discussion

4.1 Effect of Error on Hypercomplex Automatic Differentiation-Uncertainty Quantification.

Hypercomplex automatic differentiation calculates the most accurate numerical derivatives possible for a given model. However, since FEA models contain epistemic error such as mesh and time discretization error, the HYPAD-based derivatives will not be as accurate as the analytical derivatives. Furthermore, the error will accumulate with the order of the derivative since higher-order derivatives depend on lower-order derivatives. As a result, the error in the derivatives will limit the accuracy of the Taylor series expansion and consequently the statistical moments, see Fig. 6, for example. In this case, HYPAD-based derivatives higher than second-order do not increase the accuracy in the expected value and variance approximations. In contrast, the analytical derivatives show a substantial increase in accuracy of these moments for higher-order expansions of the analytical Taylor series.

Increasing the order of expansion of the Taylor series does not guarantee a monotonic decrease in error, even if the derivative values are exact. The polynomial fit produced by the Taylor series may yield large errors for certain functions and/or parameters of random variables. For example, the Taylor series predictions of temperature yield large errors compared to the analytical solution at steady-state, see Fig. 14(b), for a few random variable parameters. The error increases with order of expansion for these points, although, the majority of evaluations approach the analytical solution as the order increases.

4.2 Computational Cost of Hypercomplex Automatic Differentiation-Uncertainty Quantification.

The computational cost of HYPAD-UQ depends on the number of sensitivities to be computed which depends on the number of variables r and the order of Taylor series expansion n. The number of partial derivatives in an nth-order Taylor series expansion is $(r + n)! / r! n!$ ⁠. As a result, calculation of a high-order Taylor series expansions for a large number of variables may become infeasible.

Increasing the order of the expansion can increase the accuracy of the Taylor series. However, the amount of additional computational resources needed to calculate all higher-order partial derivatives required to increase the order of expansion may not outweigh the increase in accuracy. As shown in Fig. 4, the CPU time to calculate all partial derivatives with OTI-based HYPAD increases by a factor between 1.5 and 3 as the order of derivative increases. The HYPAD-UQ approximations of central moments calculated in Sec. 3 show that fifth- and higher-order expansions resulted in a negligible increase in accuracy for the low variance case, see Fig. 6, and caused a decrease in accuracy for the large variance case, see Fig. 11. Note that this is caused by the fit of the Taylor series expansion and not the error in FEA-based HYPAD derivatives, as shown by the comparison to the analytical Taylor series.

Although HYPAD was implemented in a simple FEA simulation in this work, the method can be applied in an identical fashion to more complex simulations as the method is general and can be applied to any numerical model. However, as the complexity of the simulation increases, the additional time to compute derivatives will increase as well.

Only a few variables will be important in many engineering problems. As a consequence, there will be terms in the Taylor series that do not need to be computed. Furthermore, the number of nonimportant terms in the Taylor series will increase as the order of expansion increases. Therefore, calculating only the important derivatives can lead to a substantial increase in computational efficiency of HYPAD-UQ if multicomplex-based or multidual-based differentiation is used. As an example, consider the variance approximation of T at steady-state for the large variance case presented in Sec. 3.4.2. The main effect and two-variable interaction effect Sobol' indices indicate that ρ and C_p are noninfluential, see Fig. 16. In addition, only six out of the 28 two-variable interaction effects that have more than a 2% contribution to the total variance. A sparse second-order HYPAD-UQ approximation of the total variance at steady-state, using only the five first-order and 11 second-order sensitivities corresponding to the Sobol' indices shown in Fig. 16, yielded 62.52% of the total variance approximated from random samples of the analytical solution. In comparison, a full second-order HYPAD-UQ approximation using all seven first-order and 35 second-order derivatives resulted in 62.77% of the total variance. Furthermore, a sparse third-order HYPAD-UQ approximation of the total variance at steady-state, using only the five first-order, 11 second-order, and 21 third-order sensitivities corresponding to the Sobol' indices shown in Fig. 16, yielded 91.57% of the total variance. In comparison, a full third-order HYPAD-UQ approximation using seven first-order, 28 second-order, and 84 third-order derivatives resulted in 90.71% of the total variance.

5 Conclusion

A new UQ method called HYPAD-UQ was developed which uses higher-order derivatives from a finite element formulation to construct Taylor series surrogate models. These surrogates are used to quickly compute the statistical moments, PDF, and Sobol' indices of the model outputs. The basis for the method is the use of hypercomplex OTI algebra to compute accurate derivatives in an efficient manner.

The HYPAD-UQ method was applied to a linear transient thermal finite element analysis to demonstrate the methodology and compare the computational performance and accuracy against the regression-based stochastic perturbation finite element method in which unknown coefficients of a complete nth-degree polynomial basis were solved using the ordinary least squares method. First-order through seventh-order HYPAD-based Taylor series expansions were used to approximate central moments, the PDF, and Sobol' indices of temperature in a heated fin for two case studies of random variable parameters. The first case study assumed Gaussian independent random variables with low variation. For this example, HYPAD-UQ achieved similar accuracy compared to the regression-based stochastic perturbation finite element method in approximately an order of magnitude less CPU time. The second case assumed non-Gaussian independent random variables with a larger variation than the previous case. HYPAD-UQ achieved significantly better accuracy than the regression-based method in approximately an order of magnitude less CPU time. Future work will compare the results of HYPAD-UQ with other state-of-the-art surrogate models such as polynomial chaos expansion and random Gaussian processes.

Sobol' indices calculated with HYPAD-UQ were shown to have high accuracy in both case studies of the numerical example. HYPAD-UQ correctly identified important main and interaction effects for the non-Gaussian large variance case study with only a second-order Taylor series expansion. The specific implementation of HYPAD used in this work took an additional CPU time of five times a traditional FEA run to compute all partial derivatives in the complete second-order Taylor series expansion. HYPAD-UQ approximations of Sobol' indices achieved far greater accuracy than the regression-based stochastic perturbation finite element method using 206 FEA evaluations.

In addition, a comparison was made to the Taylor series constructed with analytical derivatives in order to show the effect of the error in derivatives obtained from the FEA model on the accuracy of HYPAD-UQ. In the low variance case, error in the HYPAD-based derivatives limited the accuracy in central moments. However, in the large variance case, the analytic Taylor series showed only marginal improvement over HYPAD-UQ.

Funding Data

The United States Army Research Office (Grant No. W911NF2010315; Funder ID: 10.13039/100000183).
The Department of Energy/National Nuclear Administration (Award No. DE-NA0003948; Funder ID: 10.13039/100000015).

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Appendix: Central Moments and Sobol' Indices of Higher-Order Taylor Series Expansions

The following shows the expressions the expected value and Sobol' indices of third-order Taylor series expansions as well as third and fourth central moments of second-order expansions, for any distribution of the random variables in x. The expected value of a third-order Taylor series expansion is

E [Y_{3}] = E [Y_{2}] + \sum_{j_{1} = 1}^{r} \frac{1}{6} D_{j_{1} j_{1} j_{1}}^{(3)} μ_{3 j_{1}}

(A1)

The Sobol' indices of a third-order Taylor series expansion are expressed as

S_{i} = \frac{V_{i} [Y_{3}]}{V [Y_{3}]}

(A2)

S_{i j} = \frac{V_{i j} [Y_{3}]}{V [Y_{3}]}

(A3)

S_{ijk} = \frac{V_{ijk} [Y_{3}]}{V [Y_{3}]}

(A4)

where partial variances of a third-order Taylor series expansion are

\begin{matrix} V_{i} [Y_{3}] = V_{i} [Y_{2}] + \frac{1}{3} D_{i}^{(1)} D_{iii}^{(3)} μ_{4 i} + \frac{1}{6} D_{i i}^{(2)} D_{iii}^{(3)} μ_{5 i} \\ - \frac{1}{6} D_{i i}^{(2)} D_{iii}^{(3)} μ_{3 i} μ_{2 i} + \frac{1}{36} {(D_{iii}^{(3)})}^{2} μ_{6 i} - \frac{1}{36} {(D_{iii}^{(3)})}^{2} μ_{3 i}^{2} \\ + \underset{j_{1} \neq i}{\sum^{r}} {D_{i}^{(1)} D_{j_{1} j_{1} i}^{(3)} μ_{2 i} μ_{2 j_{1}} + \frac{1}{2} D_{i i}^{(2)} D_{j_{1} j_{1} i}^{(3)} μ_{3 i} μ_{2 j_{1}} + \frac{1}{6} D_{iii}^{(3)} D_{j_{1} j_{1} i}^{(3)} μ_{4 i} μ_{2 j_{1}} \\ + \frac{1}{4} {(D_{j_{1} j_{1} i}^{(3)})}^{2} μ_{2 i} μ_{2 j_{1}}^{2}} + \frac{1}{2} \sum_{j_{1} \neq i}^{r - 1} \underset{j_{2} \neq i}{\underset{j_{2} > j_{1}}{\sum^{r}}} D_{j_{1} j_{1} i}^{(3)} D_{j_{2} j_{2} i}^{(3)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 i} \end{matrix}

(A5)

\begin{matrix} V_{i j} [Y_{3}] = V_{i j} [Y_{2}] + D_{i j}^{(2)} D_{iij}^{(3)} μ_{3 i} μ_{2 j} + D_{i j}^{(2)} D_{ijj}^{(3)} μ_{2 i} μ_{3 j} + \frac{1}{4} {(D_{iij}^{(3)})}^{2} μ_{4 i} μ_{2 j} \\ - \frac{1}{4} {(D_{iij}^{(3)})}^{2} μ_{2 j} μ_{2 i}^{2} + \frac{1}{2} D_{iij}^{(3)} D_{ijj}^{(3)} μ_{3 i} μ_{3 j} + \frac{1}{4} {(D_{ijj}^{(3)})}^{2} μ_{2 i} μ_{4 j} - \frac{1}{4} {(D_{ijj}^{(3)})}^{2} μ_{2 i} μ_{2 j}^{2} \end{matrix}

(A6)

V_{ijk} [Y_{3}] = {(D_{ijk}^{(3)})}^{2} μ_{2 i} μ_{2 j} μ_{2 k}

(A7)

and the total variance of a third-order Taylor series expansion is the sum of all partial variances of the third-order Taylor series expansion, expressed as

V [Y_{3}] = \sum_{i = 1}^{r} V_{i} [Y_{3}] + \sum_{1 \leq i < j \leq r} V_{i j} [Y_{3}] + \sum_{1 \leq i < j < k \leq r} V_{ijk} [Y_{3}]

(A8)

The third central moment of a second-order Taylor series expansion is

\begin{matrix} μ_{3} (Y_{2}) = μ_{3} (Y_{1}) + \sum_{j_{1} = 1}^{r} {\frac{1}{8} {(D_{j_{1} j_{1}}^{(2)})}^{3} μ_{6 j_{1}} + \frac{1}{4} {(D_{j_{1} j_{1}}^{(2)})}^{3} μ_{2 j_{1}}^{3} - \frac{3}{2} {(D_{j_{1}}^{(1)})}^{2} D_{j_{1} j_{1}}^{(2)} μ_{2 j_{1}}^{2} \\ + \frac{3}{2} {(D_{j_{1}}^{(1)})}^{2} D_{j_{1} j_{1}}^{(2)} μ_{4 j_{1}} + \frac{3}{4} D_{j_{1}}^{(1)} {(D_{j_{1} j_{1}}^{(2)})}^{2} μ_{5 j_{1}} - \frac{3}{8} {(D_{j_{1} j_{1}}^{(2)})}^{3} μ_{2 j_{1}} μ_{4 j_{1}} - \frac{3}{2} D_{j_{1}}^{(1)} {(D_{j_{1} j_{1}}^{(2)})}^{2} μ_{2 j_{1}} μ_{3 j_{1}}} \\ + \underset{j_{1} \neq j_{2}}{\sum^{r}} {\frac{1}{2} {(D_{j_{1} j_{2}}^{(2)})}^{3} μ_{3 j_{1}} μ_{3 j_{2}} + 3 D_{j_{1}}^{(1)} {(D_{j_{1} j_{2}}^{(2)})}^{2} μ_{2 j_{2}} μ_{3 j_{1}} + \frac{3}{2} D_{j_{1} j_{1}}^{(2)} {(D_{j_{1} j_{2}}^{(2)})}^{2} μ_{2 j_{2}} μ_{4 j_{1}} \\ - \frac{3}{2} S_{j_{1} j_{1}}^{(2)} {(D_{j_{1} j_{2}}^{(2)})}^{2} μ_{2 j_{1}}^{2} μ_{2 j_{2}} + 3 D_{j_{1}}^{(1)} D_{j_{1} j_{2}}^{(2)} D_{j_{2}}^{(1)} μ_{2 j_{1}} μ_{2 j_{2}} + 3 D_{j_{1} j_{1}}^{(2)} D_{j_{1} j_{2}}^{(2)} D_{j_{2}}^{(1)} μ_{2 j_{2}} μ_{3 j_{1}} \\ + \frac{3}{4} D_{j_{1} j_{1}}^{(2)} D_{j_{1} j_{2}}^{(2)} D_{j_{2} j_{2}}^{(2)} μ_{3 j_{1}} μ_{3 j_{2}}} + \underset{j_{1} \neq j_{2} \neq j_{3}}{\sum^{r}} D_{j_{1} j_{2}}^{(2)} D_{j_{1} j_{3}}^{(2)} D_{j_{2} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} \end{matrix}

(A9)

The fourth central moment of a second-order Taylor series expansion is

\begin{matrix} μ_{4} (Y_{2}) = μ_{4} (Y_{1}) + \sum_{j_{1} = 1}^{r} {\frac{1}{16} {(D_{j_{1} j_{1}}^{(2)})}^{4} μ_{8 j_{1}} - \frac{3}{16} {(D_{j_{1} j_{1}}^{(2)})}^{4} μ_{2 j_{1}}^{4} + \frac{3}{2} {(D_{j_{1}}^{(1)})}^{2} {(D_{j_{1} j_{1}}^{(2)})}^{2} μ_{6 j_{1}} \\ + \frac{3}{8} {(D_{j_{1} j_{1}}^{(2)})}^{4} μ_{2 j_{1}}^{2} μ_{4 j_{1}} + 2 {(D_{j_{1}}^{(1)})}^{3} D_{j_{1} j_{1}}^{(2)} μ_{5 j_{1}} + \frac{1}{2} D_{j_{1}}^{(1)} {(D_{j_{1} j_{1}}^{(2)})}^{3} μ_{7 j_{1}} \\ + \frac{3}{2} {(D_{j_{1}}^{(1)})}^{2} {(D_{j_{1} j_{1}}^{(2)})}^{2} μ_{2 j_{1}}^{3} - 2 {(D_{j_{1}}^{(1)})}^{3} D_{j_{1} j_{1}}^{(2)} μ_{2 j_{1}} μ_{3 j_{1}} - \frac{3}{2} D_{j_{1}}^{(1)} {(D_{j_{1} j_{1}}^{(2)})}^{3} μ_{2 j_{1}} μ_{5 j_{1}} \\ + \frac{3}{2} D_{j_{1}}^{(1)} {(D_{j_{1} j_{1}}^{(2)})}^{3} μ_{2 j_{1}}^{2} μ_{3 j_{1}} - 3 {(D_{j_{1}}^{(1)})}^{2} {(D_{j_{1} j_{1}}^{(2)})}^{2} μ_{2 j_{1}} μ_{4 j_{1}} - \frac{1}{4} {(D_{j_{1} j_{1}}^{(2)})}^{4} μ_{2 j_{1}} μ_{6 j_{1}}} \\ + \underset{j_{1} \neq j_{2}}{\sum^{r}} {\frac{1}{2} {(D_{j_{1} j_{2}}^{(2)})}^{4} μ_{4 j_{1}} μ_{4 j_{2}} + 4 {(D_{j_{1} j_{2}}^{(2)})}^{3} D_{j_{2}}^{(1)} μ_{3 j_{1}} μ_{4 j_{2}} \\ + 2 {(D_{j_{1} j_{2}}^{(2)})}^{3} D_{j_{2} j_{2}}^{(2)} μ_{3 j_{1}} μ_{5 j_{2}} + \frac{3}{16} {(D_{j_{1} j_{1}}^{(2)})}^{2} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{1}}^{2} μ_{2 j_{2}}^{2} \\ + \frac{3}{2} {(D_{j_{1}}^{(1)})}^{2} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{1}} μ_{4 j_{2}} + 6 {(D_{j_{1} j_{2}}^{(2)})}^{2} {(D_{j_{2}}^{(1)})}^{2} μ_{2 j_{1}} μ_{4 j_{2}} \\ - \frac{3}{2} {(D_{j_{1}}^{(1)})}^{2} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}}^{2} - \frac{3}{8} {(D_{j_{1} j_{1}}^{(2)})}^{2} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{1}}^{2} μ_{4 j_{2}} \end{matrix}

\begin{matrix} + \frac{3}{16} {(D_{j_{1} j_{1}}^{(2)})}^{2} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{4 j_{1}} μ_{4 j_{2}} + \frac{3}{2} {(D_{j_{1} j_{2}}^{(2)})}^{2} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}}^{3} \\ - 3 {(D_{j_{1} j_{2}}^{(2)})}^{2} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{4 j_{2}} + \frac{3}{2} D_{j_{1} j_{1}}^{(2)} {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{1}}^{2} μ_{2 j_{2}}^{2} \\ + 6 D_{j_{1}}^{(1)} {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{2}}^{(1)} μ_{3 j_{1}} μ_{3 j_{2}} + 6 D_{j_{1} j_{1}}^{(2)} D_{j_{1} j_{2}}^{(2)} {(D_{j_{2}}^{(1)})}^{2} μ_{3 j_{1}} μ_{3 j_{2}} \\ + 3 D_{j_{1}}^{(1)} D_{j_{1} j_{2}}^{(2)} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{1}} μ_{5 j_{2}} + 6 D_{j_{1}}^{(1)} {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{2} j_{2}}^{(2)} μ_{3 j_{1}} μ_{4 j_{2}} \\ + \frac{3}{2} D_{j_{1}}^{(1)} D_{j_{1} j_{1}}^{(2)} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{3 j_{1}} μ_{4 j_{2}} + \frac{3}{2} D_{j_{1} j_{1}}^{(2)} D_{j_{1} j_{2}}^{(2)} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{3 j_{1}} μ_{5 j_{2}} \\ + \frac{3}{2} D_{j_{1} j_{1}}^{(2)} {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{2} j_{2}}^{(2)} μ_{4 j_{1}} μ_{4 j_{2}} + 6 {(D_{j_{1}}^{(1)})}^{2} D_{j_{2}}^{(1)} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{1}} μ_{3 j_{2}} \\ - 2 {(D_{j_{1} j_{2}}^{(2)})}^{3} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{2}} μ_{3 j_{1}} μ_{3 j_{2}} - 3 D_{j_{1}}^{(1)} D_{j_{1} j_{1}}^{(2)} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{2}}^{2} μ_{3 j_{1}} \\ + 6 {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{2}}^{(1)} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{1}} μ_{5 j_{2}} - 6 D_{j_{1}}^{(1)} {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{2}}^{2} μ_{3 j_{1}} \\ - 3 D_{j_{1} j_{1}}^{(2)} {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{1}}^{2} μ_{4 j_{2}} + \frac{3}{2} {(D_{j_{1} j_{1}}^{(2)})}^{2} D_{j_{2}}^{(1)} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{1}}^{2} μ_{3 j_{2}} \\ - 3 D_{j_{1} j_{1}}^{(2)} D_{j_{1} j_{2}}^{(2)} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{2}} μ_{3 j_{1}} μ_{3 j_{2}} - 6 {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{2}}^{(1)} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{3 j_{2}} \\ + 3 D_{j_{1}}^{(1)} D_{j_{1} j_{1}}^{(2)} D_{j_{2}}^{(1)} D_{j_{2} j_{2}}^{(2)} μ_{3 j_{1}} μ_{3 j_{2}} + 12 D_{j_{1}}^{(1)} D_{j_{1} j_{2}}^{(2)} D_{j_{2}}^{(1)} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{1}} μ_{4 j_{2}} \\ - 6 D_{j_{1}}^{(1)} D_{j_{1} j_{2}}^{(2)} {(D_{j_{2} j_{2}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{3 j_{2}} + 6 D_{j_{1} j_{1}}^{(2)} D_{j_{1} j_{2}}^{(2)} D_{j_{2}}^{(1)} D_{j_{2} j_{2}}^{(2)} μ_{3 j_{1}} μ_{4 j_{2}} \\ - 12 D_{j_{1}}^{(1)} D_{j_{1} j_{2}}^{(2)} D_{j_{2}}^{(1)} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}}^{2} - 6 D_{j_{1} j_{1}}^{(2)} D_{j_{1} j_{2}}^{(2)} D_{j_{2}}^{(1)} D_{j_{2} j_{2}}^{(2)} μ_{2 j_{2}}^{2} μ_{3 j_{1}} \\ + 12 D_{j_{1}}^{(1)} D_{j_{1} j_{2}}^{(2)} {(D_{j_{2}}^{(1)})}^{2} μ_{2 j_{1}} μ_{3 j_{2}}} \\ + \underset{j_{1} \neq j_{2} \neq j_{3}}{\sum^{r}} {3 {(D_{j_{1} j_{3}}^{(2)})}^{2} {(D_{j_{2} j_{3}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{4 j_{3}} + 3 {(D_{j_{1} j_{2}}^{(2)})}^{2} {(D_{j_{3}}^{(1)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} \\ + \frac{3}{4} {(D_{j_{1} j_{2}}^{(2)})}^{2} {(D_{j_{3} j_{3}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{4 j_{3}} - \frac{3}{4} {(D_{j_{1} j_{2}}^{(2)})}^{2} {(D_{j_{3} j_{3}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}}^{2} \\ + \frac{15}{4} D_{j_{1} j_{1}}^{(2)} D_{j_{1} j_{2}}^{(2)} {(D_{j_{2} j_{3}}^{(2)})}^{2} μ_{3 j_{1}} μ_{2 j_{3}} μ_{3 j_{2}} + \frac{9}{4} D_{j_{1} j_{2}}^{(2)} {(D_{j_{1} j_{3}}^{(2)})}^{2} D_{j_{2} j_{2}}^{(2)} μ_{3 j_{1}} μ_{2 j_{3}} μ_{3 j_{2}} \\ + 6 {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{1} j_{3}}^{(2)} D_{j_{2} j_{3}}^{(2)} μ_{3 j_{1}} μ_{2 j_{3}} μ_{3 j_{2}} + 12 {(D_{j_{1} j_{3}}^{(2)})}^{2} D_{j_{2}}^{(1)} D_{j_{2} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{3 j_{3}} \\ + 3 D_{j_{1} j_{1}}^{(2)} D_{j_{2} j_{2}}^{(2)} {(D_{j_{3}}^{(1)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} - 3 D_{j_{1} j_{1}}^{(2)} {(D_{j_{2}}^{(1)})}^{2} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} \\ - 2 D_{j_{1} j_{1}}^{(2)} {(D_{j_{2} j_{3}}^{(2)})}^{2} D_{j_{3}}^{(1)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{3 j_{3}} + 2 {(D_{j_{1} j_{3}}^{(2)})}^{2} D_{j_{2} j_{2}}^{(2)} D_{j_{3}}^{(1)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{3 j_{3}} \\ - D_{j_{1} j_{1}}^{(2)} {(D_{j_{2} j_{3}}^{(2)})}^{2} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{4 j_{3}} + {(D_{j_{1} j_{3}}^{(2)})}^{2} D_{j_{2} j_{2}}^{(2)} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{4 j_{3}} \\ + 3 {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{3}}^{(1)} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{3 j_{3}} + \frac{3}{4} D_{j_{1} j_{1}}^{(2)} D_{j_{2} j_{2}}^{(2)} {(D_{j_{3} j_{3}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}}^{2} \\ - \frac{3}{4} D_{j_{1} j_{1}}^{(2)} {(D_{j_{2} j_{2}}^{(2)})}^{2} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}}^{2} μ_{2 j_{3}} + D_{j_{1} j_{1}}^{(2)} {(D_{j_{2} j_{3}}^{(2)})}^{2} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}}^{2} \\ - {(D_{j_{1} j_{3}}^{(2)})}^{2} D_{j_{2} j_{2}}^{(2)} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}}^{2} - 6 D_{j_{1} j_{2}}^{(2)} D_{j_{1} j_{3}}^{(2)} D_{j_{2} j_{3}}^{(2)} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}}^{2} \\ + 3 D_{j_{1} j_{1}}^{(2)} D_{j_{1} j_{3}}^{(2)} D_{j_{2} j_{2}}^{(2)} D_{j_{2} j_{3}}^{(2)} μ_{3 j_{1}} μ_{2 j_{3}} μ_{3 j_{2}} - 6 D_{j_{1}}^{(1)} D_{j_{1} j_{2}}^{(2)} D_{j_{2}}^{(1)} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} \\ + 12 D_{j_{1} j_{2}}^{(2)} D_{j_{1} j_{3}}^{(2)} D_{j_{2}}^{(1)} D_{j_{3}}^{(1)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} + 12 D_{j_{1} j_{2}}^{(2)} D_{j_{1} j_{3}}^{(2)} D_{j_{2}}^{(1)} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{3 j_{3}} \\ + 12 D_{j_{1} j_{2}}^{(2)} D_{j_{1} j_{3}}^{(2)} D_{j_{2} j_{3}}^{(2)} D_{j_{3}}^{(1)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{3 j_{3}} + 6 D_{j_{1} j_{2}}^{(2)} D_{j_{1} j_{3}}^{(2)} D_{j_{2} j_{3}}^{(2)} D_{j_{3} j_{3}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{4 j_{3}} \\ + 6 D_{j_{1} j_{1}}^{(2)} D_{j_{2}}^{(1)} D_{j_{2} j_{3}}^{(2)} D_{j_{3}}^{(1)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}}} \\ + \underset{j_{1} \neq j_{2} \neq j_{3} \neq j_{4}}{\sum^{r}} {\frac{3}{4} {(D_{j_{1} j_{2}}^{(2)})}^{2} {(D_{j_{3} j_{4}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} μ_{2 j_{4}} + \frac{9}{8} D_{j_{1} j_{1}}^{(2)} D_{j_{2} j_{2}}^{(2)} {(D_{j_{3} j_{4}}^{(2)})}^{2} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} μ_{2 j_{4}} \\ - \frac{3}{4} D_{j_{1} j_{1}}^{(2)} {(D_{j_{2} j_{3}}^{(2)})}^{2} D_{j_{4} j_{4}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} μ_{2 j_{4}} - \frac{3}{8} {(D_{j_{1} j_{2}}^{(2)})}^{2} D_{j_{3} j_{3}}^{(2)} D_{j_{4} j_{4}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} μ_{2 j_{4}} \\ + \frac{3}{2} D_{j_{1} j_{2}}^{(2)} D_{j_{1} j_{3}}^{(2)} D_{j_{2} j_{4}}^{(2)} D_{j_{3} j_{4}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} μ_{2 j_{4}} + \frac{3}{2} D_{j_{1} j_{2}}^{(2)} D_{j_{1} j_{4}}^{(2)} D_{j_{2} j_{3}}^{(2)} D_{j_{3} j_{4}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} μ_{2 j_{4}} \\ + \frac{1}{2} D_{j_{1} j_{3}}^{(2)} D_{j_{1} j_{4}}^{(2)} D_{j_{2} j_{2}}^{(2)} D_{j_{3} j_{4}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} μ_{2 j_{4}} + D_{j_{1} j_{1}}^{(2)} D_{j_{2} j_{3}}^{(2)} D_{j_{2} j_{4}}^{(2)} D_{j_{3} j_{4}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} μ_{2 j_{4}} \\ - \frac{3}{2} D_{j_{1} j_{2}}^{(2)} D_{j_{1} j_{3}}^{(2)} D_{j_{2} j_{3}}^{(2)} D_{j_{4} j_{4}}^{(2)} μ_{2 j_{1}} μ_{2 j_{2}} μ_{2 j_{3}} μ_{2 j_{4}}} \end{matrix}

(A10)

References

1.

Xiu

,

D.

,

2010

,

Numerical Methods for Stochastic Computations

,

Princeton University Press

,

Princeton, NJ

.

2.

DiazDelaO

,

F. A.

, and

Adhikari

,

S.

,

2011

, “

Gaussian Process Emulators for the Stochastic Finite Element Method

,”

Int. J. Numer. Methods Eng.

,

87

(

6

), pp.

521

–

540

.10.1002/nme.3116

Google Scholar

Crossref

3.

Tukey

,

J. W.

,

1957

, “

The Propagation of Errors, Fluctuations, and Tolerances Basic Generalized Formulas

,” Statistical Techniques Research Group, Princeton, NJ, Technical Report No. 10.

4.

Kaminski

,

M.

,

2013

,

The Stochastic Perturbation Method for Computational Mechanics

,

Wiley

,

Hoboken, NJ

.

Google Scholar

Crossref

5.

Kaminski

,

M.

,

2015

, “

On the Dual Iterative Stochastic Perturbation-Based Finite Element Method in Solid Mechanics With Gaussian Uncertainties

,”

Int. J. Numer. Methods Eng.

,

104

(

11

), pp.

1038

–

1060

.10.1002/nme.4976

Google Scholar

Crossref

6.

Cacuci

,

D. G.

,

2003

,

Sensitivity and Uncertainty Analysis

, Vol.

1

,

CRC Press

,

Boca Raton, FL

.

Google Scholar

Crossref

7.

Paudel

,

A.

,

Gupta

,

S.

,

Thapa

,

M.

,

Mulani

,

S. B.

, and

Walters

,

R. W.

,

2022

, “

Higher-Order Taylor Series Expansion for Uncertainty Quantification With Efficient Local Sensitivity

,”

Aerosp. Sci. Technol.

,

126

, p.

107574

.10.1016/j.ast.2022.107574

Google Scholar

Crossref

8.

Hien

,

T. D.

, and

Kleiber

,

M.

,

1997

, “

Stochastic Finite Element Modelling in Linear Transient Heat Transfer

,”

Comput. Methods Appl. Mech. Eng.

,

144

(

1–2

), pp.

111

–

124

.10.1016/S0045-7825(96)01168-1

9.

Tsay

,

J. J.

, and

Arora

,

J. S.

,

1990

, “

Nonlinear Structural Design Sensitivity Analysis for Path Dependent Problems. Part 1: General Theory

,”

Comput. Methods Appl. Mech. Eng.

,

81

(

2

), pp.

183

–

208

.10.1016/0045-7825(90)90109-Y

Google Scholar

Crossref

10.

Cacuci

,

D. G.

, and

Favorite

,

J. A.

,

2018

, “

Second-Order Sensitivity Analysis of Uncollided Particle Contributions to Radiation Detector Responses

,”

Nucl. Sci. Eng.

,

190

(

2

), pp.

105

–

133

.10.1080/00295639.2018.1426899

Google Scholar

Crossref

11.

Press

,

W. H.

,

Teukolsky

,

S. A.

,

Vetterling

,

W. T.

, and

Flannery

,

B. P.

,

2007

,

Numerical Recipes: The Art of Scientific Computing

, 3rd ed.,

Cambridge University Press

, Cambridge, UK.

12.

Squire

,

W.

, and

Trapp

,

G.

,

1998

, “

Using Complex Variables to Estimate Derivatives of Real Functions

,”

SIAM Rev.

,

40

(

1

), pp.

110

–

112

.10.1137/S003614459631241X

Google Scholar

Crossref

13.

Fielder

,

R.

,

Millwater

,

H.

,

Montoya

,

A.

, and

Golden

,

P.

,

2019

, “

Efficient Estimate of Residual Stress Variance Using Complex Variable Finite Element Methods

,”

Int. J. Pressure Vessels Piping

,

173

, pp.

101

–

113

.10.1016/j.ijpvp.2019.05.004

Google Scholar

Crossref

14.

Aguirre-Mesa

,

A. M.

,

Garcia

,

M. J.

, and

Millwater

,

H.

,

2020

, “

Multiz: A Library for Computation of High-Order Derivatives Using Multicomplex or Multidual Numbers

,”

ACM Trans. Math. Software

,

46

(

3

), pp.

1

–

30

.10.1145/3378538

Google Scholar

Crossref

15.

Aristizabal

,

M.

,

Ramirez-Tamayo

,

D.

,

Garcia

,

M.

,

Aguirre-Mesa

,

A.

,

Montoya

,

A.

, and

Millwater

,

H.

,

2019

, “

Quaternion and Octonion-Based Finite Element Analysis Methods for Computing Multiple First Order Derivatives

,”

J. Comput. Phys.

,

397

, p.

108831

.10.1016/j.jcp.2019.07.030

Google Scholar

Crossref

16.

Cano

,

M. A.

,

2020

, “

Order Truncated Imaginary Algebra for Computation of Multivariable High-Order Derivatives in Finite Element Analysis

,” Ph.D. thesis,

Universidad EAFIT, Medellín, Colombia

.

17.

Lantoine

,

G.

,

Russell

,

R. P.

, and

Dargent

,

T.

,

2012

, “

Using Multicomplex Variables for Automatic Computation of High-Order Derivatives

,”

ACM Trans. Math. Software

,

38

(

3

), pp.

1

–

21

.10.1145/2168773.2168774

Google Scholar

Crossref

18.

Fike

,

J. A.

, and

Alonso

,

J. J.

,

2011

, “

The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations

,”

AIAA

Paper No. 2011-886. 10.2514/6.2011-886

19.

Aguirre-Mesa

,

A. M.

,

Ramirez-Tamayo

,

D.

,

Garcia

,

M. J.

,

Montoya

,

A.

, and

Millwater

,

H.

,

2019

, “

A Stiffness Derivative Local Hypercomplex-Variable Finite Element Method for Computing the Energy Release Rate

,”

Eng. Fract. Mech.

,

218

, p.

106581

.10.1016/j.engfracmech.2019.106581

Google Scholar

Crossref

20.

Aguirre-Mesa

,

A. M.

,

Garcia

,

M. J.

,

Aristizabal

,

M.

,

Wagner

,

D.

,

Ramirez-Tamayo

,

D.

,

Montoya

,

A.

, and

Millwater

,

H.

,

2021

, “

A Block Forward Substitution Method for Solving the Hypercomplex Finite Element System of Equations

,”

Comput. Methods Appl. Mech. Eng.

,

387

, p.

114195

.10.1016/j.cma.2021.114195

Google Scholar

Crossref

21.

Kaminski

,

M.

,

2022

, “

Uncertainty Analysis in Solid Mechanics With Uniform and Triangular Distributions Using Stochastic Perturbation-Based Finite Element Method

,”

Finite Elem. Anal. Des.

,

200

, p.

103648

.10.1016/j.finel.2021.103648

Google Scholar

Crossref

22.

Sobol

,

I.

,

1993

, “

Sensitivity Estimates for Nonlinear Mathematical Models

,”

Math. Modell. Comput. Exp.

,

1

(

4

), pp.

407

–

414

.

23.

Sudret

,

B.

,

2008

, “

Global Sensitivity Analysis Using Polynomial Chaos Expansions

,”

Reliab. Eng. Syst. Saf.

,

93

(

7

), pp.

964

–

979

.10.1016/j.ress.2007.04.002

Google Scholar

Crossref

24.

Carneiro

,

G. D. N.

, and

Antonio

,

C. C.

,

2019

, “

Sobol' Indices as Dimension Reduction Technique in Evolutionary-Based Reliability Assessment

,”

Eng. Comput.

,

37

(

1

), pp.

368

–

398

.10.1108/EC-03-2019-0113

Google Scholar

Crossref

25.

Borgonovo

,

E.

, and

Smith

,

C. L.

,

2011

, “

A Study of Interactions in the Risk Assessment of Complex Engineering Systems: An Application to Space PSA

,”

Oper. Res.

,

59

(

6

), pp.

1461

–

1476

.10.1287/opre.1110.0973

Google Scholar

Crossref

26.

Sakata

,

S.-I.

, and

Torigoe

,

I.

,

2015

, “

A Successive Perturbation-Based Multiscale Stochastic Analysis Method for Composite Materials

,”

Finite Elem. Anal. Des.

,

102–103

, pp.

74

–

84

.10.1016/j.finel.2015.05.001

27.

Bouhlel

,

M. A.

,

Hwang

,

J. T.

,

Bartoli

,

N.

,

Lafage

,

R.

,

Morlier

,

J.

, and

Martins

,

J. R. R. A.

,

2019

, “

A Python Surrogate Modeling Framework With Derivatives

,”

Adv. Eng. Software

,

135

, p.

102662

.10.1016/j.advengsoft.2019.03.005

Google Scholar

Crossref

28.

Guo

,

L.

,

Narayan

,

A.

, and

Zhou

,

T.

,

2018

, “

A Gradient Enhanced l1-Minimization for Sparse Approximation of Polynomial Chaos Expansions

,”

J. Comput. Phys.

,

367

, pp.

49

–

64

.10.1016/j.jcp.2018.04.026

Google Scholar

Crossref

29.

Hart

,

J.

,

Waanders

,

B. V. B.

, and

Herzog

,

R.

,

2019

, “

Hyper-Differential Sensitivity Analysis of Uncertain Parameters in PDE-Constrained Optimization

,”

Int. J. Uncertainty Quantif.

, 10(3), pp.

225

–

248

.10.1615/Int.J.UncertaintyQuantification.2020032480

30.

Rabitz

,

H.

, and

Aliş

,

Ö. F.

,

1999

, “

General Foundations of High-Dimensional Model Representations

,”

J. Math. Chem.

,

25

(

2/3

), pp.

197

–

233

.10.1023/A:1019188517934

Google Scholar

Crossref

31.

Sobol

,

I. M.

,

2001

, “

Global Sensitivity Indices for Nonlinear Mathematical Models and Their Monte Carlo Estimates

,”

Math. Comput. Simul.

,

55

(

1

), pp.

271

–

280

.10.1016/S0378-4754(00)00270-6

32.

Alvarez

,

E. J.

, and

Ning

,

A.

,

2020

, “

High-Fidelity Modeling of Multirotor Aerodynamic Interactions for Aircraft Design

,”

AIAA J.

,

58

(

10

), pp.

4385

–

4400

.10.2514/1.J059178

Google Scholar

Crossref

33.

Luo

,

R.

,

Xu

,

W.

,

Shao

,

T.

,

Xu

,

H.

, and

Yang

,

Y.

,

2019

, “

Accelerated Complex-Step Finite Difference for Expedient Deformable Simulation

,”

ACM Trans. Graphics

,

38

(

6

), pp.

1

–

16

.10.1145/3355089.3356493

Google Scholar

Crossref

34.

Olivares

,

H.

,

Porth

,

O.

,

Davelaar

,

J.

,

Most

,

E. R.

,

Fromm

,

C. M.

,

Mizuno

,

Y.

,

Younsi

,

Z.

, and

Rezzolla

,

L.

,

2019

, “

Constrained Transport and Adaptive Mesh Refinement in the Black Hole Accretion Code

,”

Astron. Astrophys.

,

629

, p.

A61

.10.1051/0004-6361/201935559

Google Scholar

Crossref

35.

Martins

,

J. R. R. A.

,

Sturdza

,

P.

, and

Alonso

,

J. J.

,

2003

, “

The Complex-Step Derivative Approximation

,”

ACM Trans. Math. Software

,

29

(

3

), pp.

245

–

262

.10.1145/838250.838251

Google Scholar

Crossref

36.

Oberbichler

,

T.

, Wuchner, R., and Bletzinger, K.-U., 2021, “

Efficient Computation of Nonlinear Isogeometric Elements Using the Adjoint Method and Algorithmic Differentiation

,”

Comput. Methods Appl. Mech. Eng.

, 381, p.

113817

.10.1016/j.cma.2021.113817

37.

Casado

,

J. M. V.

, and

Hewson

,

R.

,

2020

, “

Algorithm 1008: Multicomplex Number Class for Matlab, With a Focus on the Accurate Calculation of Small Imaginary Terms for Multicomplex Step Sensitivity Calculations

,”

ACM Trans. Math. Software

,

46

(

2

), pp.

1

–

26

.10.1145/3378542

Google Scholar

Crossref

38.

Revels

,

J.

,

Lubin

,

M.

, and

Papamarkou

,

T.

,

2016

, “

Forward-Mode Automatic Differentiation in Julia

,” e-print

arXiv:1607.07892

.https://arxiv.org/abs/1607.07892

39.

Balcer

,

M. R.

,

Millwater

,

H.

, and

Favorite

,

J. A.

,

2021

, “

Multidual Sensitivity Method in Ray-Tracing Transport Simulations

,”

Nucl. Sci. Eng.

,

195

(

9

), pp.

907

–

936

.10.1080/00295639.2021.1883949

Google Scholar

Crossref

40.

Fish

,

J.

, and

Belytschko

,

T.

,

2007

,

A First Course in Finite Elements

,

Wiley

, Chichester, UK.

Google Scholar

Crossref

41.

Ramirez-Tamayo

,

D.

,

Soulami

,

A.

,

Gupta

,

V.

,

Restrepo

,

D.

,

Montoya

,

A.

, and

Millwater

,

H.

,

2021

, “

A Complex-Variable Cohesive Finite Element Subroutine to Enable Efficient Determination of Interfacial Cohesive Material Parameters

,”

Eng. Fract. Mech.

,

247

, p.

107638

.10.1016/j.engfracmech.2021.107638

Google Scholar

Crossref

42.

Montoya

,

A.

,

Fielder

,

R.

,

Gomez-Farias

,

A.

, and

Millwater

,

H.

,

2015

, “

Finite Element Sensitivity for Plasticity Using Complex Variable Methods

,”

J. Eng. Mech.

,

141

(

2

), p.

04014118

.10.1061/(ASCE)EM.1943-7889.0000837

Google Scholar

Crossref

43.

Kraus

,

A. D.

,

Aziz

,

A.

, and

Welty

,

J.

,

2001

,

Transient Heat Transfer in Extended Surfaces

,

Wiley

, Hoboken, NJ, pp.

754

–

818

.

44.

Rincon-Tabares

,

J.-S.

,

Velasquez-Gonzalez

,

J. C.

,

Ramirez-Tamayo

,

D.

,

Montoya

,

A.

,

Millwater

,

H.

, and

Restrepo

,

D.

,

2022

, “

Sensitivity Analysis for Transient Thermal Problems Using the Complex-Variable Finite Element Method

,”

Appl. Sci.

,

12

(

5

), p.

2738

.10.3390/app12052738

Google Scholar

Crossref

45.

Mills

,

K. C.

,

2002

,

Recommended Values of Thermophysical Properties for Selected Commercial Alloys

,

Woodhead Publishing

,

Buckingham, UK

.

Google Scholar

Crossref

46.

Efron

,

B.

, and

Tibshirani

,

R. J.

,

1994

,

An Introduction to the Bootstrap

,

CRC Press

, Boca Raton, FL.

Google Scholar

Crossref

47.

Epanechnikov

,

V. A.

,

1969

, “

Non-Parametric Estimation of a Multivariate Probability Density

,”

Theory Probab. Its Appl.

,

14

(

1

), pp.

153

–

158

.10.1137/1114019

Google Scholar

Crossref

48.

Saltelli

,

A.

,

Ratto

,

M.

,

Andres

,

T.

,

Campolongo

,

F.

,

Cariboni

,

J.

,

Gatelli

,

D.

,

Saisana

,

M.

, and

Tarantola

,

S.

,

2008

,

Global Sensitivity Analysis: The Primer

,

Wiley

,

Hoboken, NJ

.

HYPAD-UQ: A Derivative-Based Uncertainty Quantification Method Using a Hypercomplex Finite Element Method

Abstract

1 Introduction

2 Methodology

2.1 Taylor Series Expansion of Central Moments and Sobol' Indices.

2.2 Hypercomplex-Based Automatic Differentiation.

2.2.1 Differentiation With Complex and Hypercomplex Algebra.

2.2.2 Hypercomplex Automatic Differentiation With Order Truncated Imaginary Algebra.

2.3 Hypercomplex Finite Element Method for Transient Heat Transfer Problems.

2.3.1 Hypercomplex Finite Element Method Formulation in Transient Heat Transfer Problems.

2.3.2 Solution to the Order Truncated Imaginary System of Equations.

3 Numerical Example: Uncertainty Quantification of Transient Heat Transfer in a Fin

3.1 Analytical Solution.

3.2 Finite Element Model.

3.3 High-Order Order Truncated Imaginary-Based Derivatives.

3.4 Uncertainty Quantification of Temperature at Tip of Fin.

3.4.1 Gaussian Random Variables With Low Coefficient of Variations.

3.4.2 Non-Gaussian Random Variables With Large Coefficient of Variations.

4 Discussion

4.1 Effect of Error on Hypercomplex Automatic Differentiation-Uncertainty Quantification.

4.2 Computational Cost of Hypercomplex Automatic Differentiation-Uncertainty Quantification.

5 Conclusion

Funding Data

Data Availability Statement

Appendix: Central Moments and Sobol' Indices of Higher-Order Taylor Series Expansions

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

HYPAD-UQ: A Derivative-Based Uncertainty Quantification Method Using a Hypercomplex Finite Element Method

Abstract

1 Introduction

2 Methodology

2.1 Taylor Series Expansion of Central Moments and Sobol' Indices.

2.2 Hypercomplex-Based Automatic Differentiation.

2.2.1 Differentiation With Complex and Hypercomplex Algebra.

2.2.2 Hypercomplex Automatic Differentiation With Order Truncated Imaginary Algebra.

2.3 Hypercomplex Finite Element Method for Transient Heat Transfer Problems.

2.3.1 Hypercomplex Finite Element Method Formulation in Transient Heat Transfer Problems.

2.3.2 Solution to the Order Truncated Imaginary System of Equations.

3 Numerical Example: Uncertainty Quantification of Transient Heat Transfer in a Fin

3.1 Analytical Solution.

3.2 Finite Element Model.

3.3 High-Order Order Truncated Imaginary-Based Derivatives.

3.4 Uncertainty Quantification of Temperature at Tip of Fin.

3.4.1 Gaussian Random Variables With Low Coefficient of Variations.

3.4.2 Non-Gaussian Random Variables With Large Coefficient of Variations.

4 Discussion

4.1 Effect of Error on Hypercomplex Automatic Differentiation-Uncertainty Quantification.

4.2 Computational Cost of Hypercomplex Automatic Differentiation-Uncertainty Quantification.

5 Conclusion

Funding Data

Data Availability Statement

Appendix: Central Moments and Sobol' Indices of Higher-Order Taylor Series Expansions

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only