A Non-Intrusive Stochastic Isogeometric Analysis of Functionally Graded Plates with Material Uncertainty

: A non-intrusive approach coupled with non-uniform rational B-splines based isogeometric ﬁnite element method is proposed here. The developed methodology was employed to study the stochastic static bending and free vibration characteristics of functionally graded material plates with inhered material randomness. A ﬁrst order shear deformation theory with an artiﬁcial shear correction factor was used for spatial discretization. The output randomness is represented by polynomial chaos expansion. The robustness and accuracy of the framework were demonstrated by comparing the results with Monte Carlo simulations. A systematic parametric study was carried out to bring out the sensitivity of the input randomness on the stochastic output response using Sobol’ indices. Functionally graded plates made up of Aluminium (Al) and Zirconium Oxide (ZrO 2 ) were considered in all the numerical examples.


Introduction
Composite materials are widely used in a large variety of structures due to their high strength to weight ratio and excellent thermo-mechanical and corrosion resistance properties [1][2][3]. On the other hand, composite materials can be engineered to achieve a varying spatial composition profile and, hence, arrive at the required properties at specific locations. Such inhomogeneous materials, consisting of two or more different materials, are known as functionally graded materials (FGM). Functionally graded materials can be fabricated by varying the elastic modulus from one phase to another in one or more directions. Some applications of FGM include rocket heat shields, electrical insulation in metal/ceramic joints, dental implants, and bone replacements [4]. The FGM fabrication process results in inherent randomness, since the gradients of the constituent materials are to be varied continuously in the space along with the physical material properties. Therefore, the analysis of FGM using deterministic approaches will not accurately predict actual behaviour. The inherent uncertainties existing in the system are required to be considered while estimating the system response. Thus, structures involving functionally graded materials are required to be modelled and analysed using stochastic analysis [5].
One straight forward approach to quantify the effect of randomness on the output response is to use a direct Monte Carlo simulation (MCS) approach. However, given that the MCS requires a sufficiently large number of samples, typically 10,000, it becomes computationally intractable. To alleviate this issue, the stochastic finite element method was found to be the most common and promising tool to study the stochastic response. Different variants of the stochastic finite element method include the probabilistic finite element method [6], the weighted integral method in stochastic finite element analysis [7], the spectral stochastic finite element method [8], and the stochastic finite element method based on perturbation technique [9], to name a few. A comprehensive discussion of stochastic finite element methods is given in [10].
A stochastic perturbation technique was adopted to study the free vibration characteristics of FGM plates in [11,12]. The spectral stochastic isogeometric analysis (SSIGA) using the first order deformation theory was implemented by Keyen Li et al. to obtain the static response of the FGM plates [13]. Chiba et al. [11], addressed the randomness in the thermal conductivity and coefficient of linear thermal expansion of the FGM plates. The analysis was performed by considering the thermal properties to be distinct for each layer. Stochastic finite element analysis using artificial neural networks (ANN) coupled with higher-order-zigzag theory (HOZT) was presented by R.R. Kumar et al. [14]. In this study, it was shown that the uncertainty in the input parameters has significant effects on the higher order modes of the buckling analysis [14]. The size-dependent effects on thermal buckling and post-buckling behaviours of the imperfect micro FGM plates with porosity is presented using isogeometric analysis in [15]. The seventh order shear deformation plate theory was employed to capture the size-dependent phenomenon. Shaker [12] studied the stochastic free vibrations of FGM plates using the third order shear deformation theory, with uncertainties in the constituent material properties and the volume fraction index. However, a subtle drawback of this approach is that it is limited to small input variance, typically less than 20% [16]. Kumaraian et al. [17] employed a cell-based smoothed discrete shear gap method to study the stochastic free vibration characteristics of functionally graded material plates. The input random field was represented by using the Karhunen-Loéve (KL) expansion, and a polynomial chaos expansion was used to represent the stochastic output response. The critical buckling behaviour of the the plate due to the varying thickness and presence of cracks was addressed using the phase field theory coupled with third order shear deformation theory [18]. Sepahvand et al. [19] studied the stochastic free vibration characteristics of orthotropic plates with uncertainty in elastic moduli using the generalized polynomial chaos expansion (gPCE) to represent both the uncertainty in input as well as in the output. The sensitivity analysis for low-frequency vibration and low-velocity impact of the FGM plates was presented by Karsh et al. [20]. The equilibrium and stability equations of FGM plates for working in the high temperature conditions were formulated by Tran et al. [21]. The optimisation using hybrid evolutionary firefly algorithm of the FGM plates subjected to thermo-mechanical loading was demonstrated by Lieu and Qui X [22].
In this work, the influence of the randomness in the material properties, particularly, Young's modulus and density on the uncertainty in the static bending and free vibration characteristics of FGM plates were studied using a non-intrusive approach. The salient feature of this approach is that it does not require any modification to the existing code. However, to capture the uncertainties, the forward solver has to be used as a black-box, which has to be computationally efficient. Hence, in this paper, we used the non-uniform rational basis spline (NURBS)-based isogeometric analysis framework for the spatial discretization, and the plate kinematics is based on the first order shear deformation theory. The results from the proposed framework are validated against the Monte Carlo Simulations with a sufficiently large number of samples. The sensitivity of the material randomness on the output characteristics are studied using Sobol' indices.
The rest of the paper is organized as follows: an overview of the functionally graded material plates, and the first order shear deformation plate theory is presented in Section 2. Section 3 presents a brief discussion of the isogeometric analysis employed in this study. The non-intrusive spectral approach is presented in Section 4. The influence of the material randomness on the static bending and free vibration of functionally graded plates are presented in Section 5. The results from the present framework are compared with the brute-force Monte Carlo simulations. The sensitivity of the random variables on the stochastic output response was systematically studied. This is followed by conclusions in Section 6.

Effective Modulus and Poisson's Ratio
In this study, a rectangular composite plate made of ceramic and metal is considered for simulations. As shown in Figure 1a, the dimensions along the in-plane directions (x, y) and along the thickness (z) direction are indicated by a, b, and h, respectively. The material on the top surface (z = h/2) of the plate is ceramic rich and is graded to metal towards the bottom surface (z = −h/2) of the plate, according to a power law. The effective properties of the plate made of functionally graded materials can be estimated by employing either the rule of mixtures or the Mori-Tanaka homogenization scheme.
Let V i (i = c, m) denotes the volume fraction of the constituent phase, where the subscripts c and m refer to ceramic and metallic phases, respectively. The total volume fraction of ceramic and metal phases is related as V c + V m = 1. Furthermore, the volume fraction of the ceramic (V c ) is given by Equation (1): where n is known as the volume fraction exponent (n ≥ 0), or the gradient index. The composition of the ceramic and the metal varies linearly when n = 1, whereas, n = 0 indicates a fully ceramic plate. Any other value of n yields a composite material with a smooth transition from ceramic to metal, as shown in Figure 1b. In this study, the Mori-Tanaka homogenization scheme is adopted to compute the effective Young's modulus and Poisson's ratio.
According to the Mori-Tanaka homogenization scheme, the effective Young's modulus (E eff ) and Poisson's ratio (ν eff ) are computed using the effective bulk modulus (K eff ) and the effective shear modulus (G eff ), as [23]: where the effective bulk and shear modulus can be estimated using the relations mentioned below: .

(5)
K m and K c are the bulk modulus and G m and G c are the shear modulus; of the metal and ceramic materials, respectively, and

First Order Shear Deformation Theory
The displacements along the x, y and z directions at a point (x, y, z) in the plate (see Figure 1a), indicated by u, v, w, can be expressed in terms of the mid-plane displacements u 0 , v 0 , w 0 and independent rotations θ x , θ y in the yz and xz planes, respectively, as: u(x, y, z, t) = u 0 (x, y, t) + zθ x (x, y, t), v(x, y, z, t) = v 0 (x, y, t) + zθ y (x, y, t), w(x, y, z, t) = w 0 (x, y, t). The strain field in terms of mid-plane deformation are given by: The mid-plane (ε p ), bending (ε b ), shear (ε s ) strains in Equation (10)  where the subscript "comma" represents the partial derivative with respect to the succeeding spatial coordinate. The resultant membrane stress (N) and the resultant bending stress (M) can be expressed in terms of the membrane strains, ε p and bending strains ε b through the following constitutive relations: where the matrices A = A ij , B = B ij and D b = D ij ; (i, j = 1, 2, 6) indicate the extensional, bending-extensional coupling, and bending stiffness coefficients, which are defined as: Similarly, the transverse shear force Q = {Q xz , Q yz } is related to the transverse shear strains (ε s ) as: is the transverse shear stiffness coefficient, υ i , υ j are the transverse shear coefficients for non-uniform shear strain distribution through the plate thickness. The stiffness coefficients Q ij are defined as: where the modulus of elasticity (E(z)) and the Poisson's ratio (ν) are estimated using Equation (3). The total strain energy (U) of the functionally graded plate can be estimated as: where δ = {u, v, w, θ x , θ y } indicate the vector of degrees of freedom. Furthermore, the strain energy function in Equation (22) can be rewritten as [24]: where K is the linear stiffness matrix. On the other hand, the kinetic energy of the plate is given by: is the mass density which varies along the thickness direction. Moreover, considering the thermal field results in the in-plane stress, N th . The external work due to the in-plane stresses developed in the plate in the presence of a thermal load is given by: The governing equations of motion can be obtained by considering the Lagrange equations of motion mentioned below The governing equations obtained using the minimization of total potential energy in Equation (26), are solved using the Galerkin finite element method. Therefore, the simplified finite element equations can be written as: Free vibration: where are the stiffness matrix, the mass matrix and the force vector, respectively, δ = u o , v o , w o , θ x , θ y is the vector of degrees of freedom and B i , i = 1, 2, 3 are the strain-displacement matrix given by: where n denotes the number of nodes.

Overview of Isogeometric Analysis
In this study, non-uniform rational basis spline (NURBS) basis functions are used in the finite element approximation. A brief introduction to NURBS is presented here, further details on their use in finite element method (FEM) are discussed in [25,26]. The key ingredients in the construction of NURBS basis functions are: the knot vector (a non decreasing sequence of parameter values, ξ i ≤ ξ i+1 , i = 0, 1, · · · , m − 1), the control points, P i , the degree of the curve p, and the weight associated to a control point, w. The i th B-spline basis function of degree p, denoted by N i,p is defined as: A p th degree NURBS curve is defined as follows: where P i are the control points and w i are the associated weights. Figure 2 shows the third order non-uniform rational B-splines for a knot vector, Ξ = {0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5}. NURBS basis functions has the following properties: (i) non-negativity, (ii) partition of unity, ∑ i N i,p = 1; (iii) interpolatory at the end points. As the same function is also used to represent the geometry, the exact representation of the geometry is preserved. It should be noted that the continuity of the NURBS functions can be tailored to the needs of the problem. The B-spline surfaces are defined by the tensor product of basis functions in two parametric dimensions ξ and η with two knot vectors, one in each dimension as: where P i,j is the bidirectional control net and N i,p and M j,q are the B-spline basis functions defined on the knot vectors over an m × n net of control points P i,j . The NURBS surface is then defined by: where w(ξ, η) is the weighting function. The displacement field within the control mesh is approximated by where u oJ , v oJ , w oJ , θ xJ , θ yJ are the nodal variables and C(ξ, η) are the basis functions given by Equation (38). The transverse shear deformations are included in the formulation of the Mindlin theory for thick plates. In the Mindlin theory, the transverse normal to the mid surface of the plate before deformation remain straight but not necessarily normal to the mid surface after deformation. Based on this condition, the continuity requirement on the assumed displacement fields can be relaxed. However, for very thin plates, the following relations must be satisfied: i.e., the shear strain ε s must vanish in the domain as the thickness approaches zero. The lower order NURBS basis functions suffer from shear locking when applied to thin plates. Kikuchi and Ishii [27] introduced an artificial shear correction factor to suppress shear locking in a 4-noded quadrilateral element. In this paper, we employ the same technique to suppress the shear locking syndrome in lower order NURBS basis functions [28]. The modified shear correction factor is given by: where n and β are positive integers, l e is the diameter of the element equal to maximum or diagonal length of the element, h is thickness of the element and υ is the shear correction factor.

Non-Intrusive Spectral Projection
Let Y = f (ζ) represents the response from a computational model, where ζ = {ζ 1 , ζ 2 , · · · , ζ d } denote the random input parameters of a probabilistic model, in which d indicates the number of random parameters. The randomness in the computational model could be due to randomness in the material, geometry, and/or the applied load. In this study, only the randomness in material properties is considered. Therefore, a finite variance of the output Y in the Hilbert space can be represented as: where y i are the polynomial coefficients and Ψ i denote the multivariate basis. An orthonormal multivariate polynomial basis is selected for the polynomial chaos expansion. The orthogonal polynomials can be obtained using the monomials {1, ζ, ζ 2 , · · · }, and the Gram-Schmidt orthogonality condition.
Norbert Wiener first proposed [29] the polynomial chaos expansion (PCE) for Gaussian random input variables in terms of Hermite polynomials, which can be used to expand the random process with finite variance, known as the second order random process. These polynomials can approximate a given functional in the L 2 space and the convergence also happens in L 2 space. Therefore, polynomial chaos expansion theory can be applied to a wide range of real-time applications. The list of random distribution functions and their corresponding optimal orthogonal polynomials are listed in Table 1. The orthogonal polynomials are univariate, and the tensor product in the multi-dimension can provide multivariate orthogonal polynomials. Let us consider the second order random process denoted by the output variable Y, which is a function of random event θ. The polynomial chaos expansion using the Legendre polynomials can be written as: where P n (ζ) denotes the Legendre polynomials of order n and the independent random variable ζ follows the uniform distribution and its dispersion around the mean is specified using the coefficient of variation (CoV). The polynomial chaos expansion is represented using the integrals, which are replaced by the summation in the discrete form. For practical applications, the infinite series in Equation (44) is truncated as mentioned below where p indicates the highest polynomial degree and d denotes the dimensionality of the random input vector. The Legendre polynomials in general can be expressed as where M = n 2 or n−1 2 whichever is an integer. A simplified form of Equation (44) considering N PCE terms is given by The Legendre polynomial basis satisfies the following orthogonality property.
where δ ij represents the Kronecker delta and , indicate the inner product in the Hilbert space determined by the support of the random variable, i.e., a, b = Ω ab dΩ The coefficients in Equation (44) are obtained using the Galerkin projection. This implies Equation (44) is multiplied by every polynomial as shown below: The residual error representing the output response, which is estimated using the PCE can be ensured to be orthogonal to the function space spanned by the PCE basis, through Equation (48). Equation (48) can be evaluated using the Gauss quadrature given by: where W is the weighting functions and N gp is the number of Gauss points. The mean (Y) and variance (V) of the output response mentioned in Equation (46), is estimated using the following orthogonality conditions: Furthermore, sensitivity analysis was performed to identify the input parameters that have greatest influence on the output response of the system. In this study, the global sensitivity analysis was performed, where the global sensitivity index focuses on the variance of the model and determines the impact of input variability on the output variance. The advantage of using PCE is that the global sensitivity functions can be directly obtained from the coefficients of PCE. The sensitivity index of the input parameters can be calculated using the Sobol indices.

Numerical Examples
In this section, stochastic static bending and free vibrations of functionally graded material plate are studied using the NURBS-based isogeometric analysis. An FGM plate made up of Aluminium (Al) and Zirconium Oxide (ZrO 2 ) was considered in the simulations. The plate is a homogeneous ceramic, when the gradient index is n = 0. Four different values for the gradient index (n = 0, 2, 5, 8) were considered.
The mean values of the Young's modulus and density of the constituent phase are given in Table 2. Poisson's ratio, ν is assumed to be deterministic and taken as 0.3 for both metallic and ceramic phases, whilst, the Young's modulus and the density of the plate material are considered to be random with a uniform distribution of the coefficient of variation (CoV) as 20% for both metallic and ceramic phases. Since a uniform distribution is considered, the orthogonal polynomial of the Askey-Wiener scheme becomes the Legendre polynomial. Unless specified, quadratic NURBS were employed in this study, considering only simply supported boundary conditions mentioned below:

Validation
Before proceeding with the stochastic analysis, results from the deterministic analysis were validated using the results from the literature. The plate was modelled with 4, 8, 16, and 24 control points per side for static bending and 8, 16, and 24 control points for free vibration analysis, respectively. Tables 3 and 4 summarize the convergence of the normalised central deflection and the first normalised fundamental frequency with increasing number of control points and for different gradient indexes. In case of the static bending, the plate is subjected to a uniformly distributed load, p = 100 N/m 2 on the top surface of the plate. Results from the present analysis were compared with the edge-based smoothed discrete shear gap method (ES-DSG3) [30], a continuum mechanics based four-node mixed interpolation of tensorial components(MITC4) [30], higher order shear deformation plate theory (HSDT) [31], and kp−Ritz [32,33] methods, where a very good agreement is observed. The analysis was carried out using quadratic NURBS with 16 control points.

Static Bending
The accuracy of the polynomial chaos expansion applied in a non-intrusive manner is illustrated by considering two cases:

•
Case A: E c is random, • Case B: E c and E m is random.
The Young's modulus is considered as the random parameter for static bending, and the normalised central deflection of the plate is obtained for gradient indices n = 0, 5. The probability density function of the normalised central deflection is shown in Figure 3, by considering the PCE of orders one and two and the Monte Carlo simulations. From Figure 3, the results are found to be in good agreement even with the crude Monte Carlo simulations. To perform the MCS, a total of 10,000 iterations are performed. However, for the PCE analysis, 4 and 16 samplings were considered for random variables one and two, respectively, which are significantly less when compared to MCS. The mean and standard deviation of the normalised central deflection are tabulated in Table 5. It is observed that the results obtained from PCE order one are comparable with MCS. This clearly shows the benefits of the non-intrusive PCE analysis over MCS, with only a few evaluations required in case of the latter. Due to the increase in metallic constituent in the plate, the mean central deflection values of the plate found to increase with the increase in the gradient index. The sensitivity analysis can investigate the influence of the metallic fraction on the plate. For this purpose, gradient indices n = 2, 5, 8 were considered. The estimated Sobol' indices are plotted in Figure 4, where the Young's modulus of the metal was found influence significantly, when compared to the ceramic irrespective of the gradient index.

Free Vibration Analysis
The influence of randomness in material properties on the first two fundamental frequencies was studied. The randomness in the Young's modulus and density, for two different gradient indices, n = 0, 5 were studied, considering the following cases:

•
Case A: E c is random.

•
Case B: E c and ρ c is random.

•
Case C: E c and E m is random.

•
Case D: E c , E m and ρ c , ρ m is random. Figures 5-8 show the probability density function for the first two normalised fundamental frequencies for two different gradient indices, n = 0 and 5 for Cases A-D. Results from the present method with PCE orders 1 and 2 are compared with the results from Monte Carlo Simulations runs of 10,000, where a close agreement is achieved. The mean frequency is observed to decrease with increasing gradient index, which is attributed to the increased metallic fraction.     When the number of random variables is four, the total number of evaluations performed using the PCE with Order 2 is 64, which is small when compared to 10,000 Monte Carlo simulations. The MCS took ≈45,062 s, whilst the PCE of order 2 took only around 189 s, which is much lower as compared to the former. In spite of this, it is observed from Figures 5-8, that the proposed framework yields more accurate results. The mean and standard deviation of the first natural frequencies obtained from MCS and the PCE of order equal 1, are tabulated in Table 6. where the results from the PCE agrees well with the results from MCS. For each of the cases considered, the Sobol' index is shown in Figure 9. For a homogeneous plate, i.e., n = 0, both the effective Young's modulus and the density of the material have equal influence When the number of random variables is four, the total number of evaluations performed using the PCE with Order 2 is 64, which is small when compared to 10,000 Monte Carlo simulations. The MCS took ≈45,062 s, whilst the PCE of order 2 took only around 189 s, which is much lower as compared to the former. In spite of this, it is observed from Figures 5-8, that the proposed framework yields more accurate results. The mean and standard deviation of the first natural frequencies obtained from MCS and the PCE of order equal 1, are tabulated in Table 6. where the results from the PCE agrees well with the results from MCS. For each of the cases considered, the Sobol' index is shown in Figure 9. For a homogeneous plate, i.e., n = 0, both the effective Young's modulus and the density of the material have equal influence on the natural frequency of the plate. on the natural frequency of the plate. However, when n > 0, the effective Young's modulus has a significant impact. When both metallic and ceramic Young's modulus are considered to be random, the randomness in Young's modulus of the metallic phase is noticed to have more influence on the normalised fundamental frequency. When all the material properties, Young's modulus, and density of both the constituents are taken as random, density of the metallic phase has a significant influence on the output variable. These observations can be attributed to the increase in the metallic phase of the plate with an increasing gradient index. The influence of CoV of the metallic phase on the mean and standard deviation of the fundamental frequency is shown in Figure 10, where the influence of the gradient index is also depicted. It is inferred that the mean value of the fundamental frequency is found to be constant, and the standard deviation increases with increasing CoV.

Conclusions
In this study, the influence of randomness in the material properties of constituent materials on the static bending and free vibrations of FGM plates is studied using a non-intrusive approach. The stochastic response and its moments are obtained using the non-intrusive technique, with the help of PCE up to two orders. Results from the present analysis are compared with with 10,000 MCS realizations. A global sensitivity analysis is performed to determine the influence of the input Figure 9. Sensitivity study of the random parameters using the first order Sobol' index for the first fundamental frequency of the functionally graded material plate.
However, when n > 0, the effective Young's modulus has a significant impact. When both metallic and ceramic Young's modulus are considered to be random, the randomness in Young's modulus of the metallic phase is noticed to have more influence on the normalised fundamental frequency. When all the material properties, Young's modulus, and density of both the constituents are taken as random, density of the metallic phase has a significant influence on the output variable. These observations can be attributed to the increase in the metallic phase of the plate with an increasing gradient index. The influence of CoV of the metallic phase on the mean and standard deviation of the fundamental frequency is shown in Figure 10, where the influence of the gradient index is also depicted. It is inferred that the mean value of the fundamental frequency is found to be constant, and the standard deviation increases with increasing CoV.

Conclusions
In this study, the influence of randomness in the material properties of constituent materials on the static bending and free vibrations of FGM plates is studied using a non-intrusive approach. The stochastic response and its moments are obtained using the non-intrusive technique, with the help of PCE up to two orders. Results from the present analysis are compared with with 10,000 MCS realizations. A global sensitivity analysis is performed to determine the influence of the input randomness on the output variance, by computing the Sobol' indices, which are determined directly from the PCE's coefficients without any additional sampling. From the numerical analysis, it is inferred that the results estimated using the developed methodology, based on the polynomial chaos expansion of the output variables are found to be in close agreement with the results from Monte Carlo simulations. Therefore, the developed framework can thus be used for studying the influence of other parameters without having to carry out the computationally intensive Monte Carlo simulations. The salient features of the proposed framework are: • It can handle a wider range of CoV, unlike the perturbation based methods.

•
Does not require to modify the existing code like the intrusive based approaches.