Perturbation based stochastic isogeometric analysis for bending of functionally graded plates with the randomness of elastic modulus

author https://doi.org/10.1590/1679-78256066 Abstract We propose, in this paper, stochastic isogeometric analysis (SIGA) is a type of non-statistic approach in which combines the perturbation technique with the standard isogeometric analysis, in particular for static behavior of functionally graded plates with the uncertain elastic modulus. We assume that the spatial random variation of elastic modulus can be modeled as a two dimensional Gaussian random field in the plane of the plate. The random field is discretized to set of random variables using the integration point method. The system equations of SIGA are created using the NURBS functions for approximation displacement fields in conjunction with the first-order and second-order perturbation expansions of random fields, stiffness matrix, displacement fields. Besides the non-statistic approach, Monte Carlo simulation is presented for validation. The accuracy and appropriateness of the non-statistic approach are demonstrated via comparisons of the present results with those given by the stochastic finite element method in the literature and by the Monte Carlo analysis as well. The numerical examples are employed to investigate the effect of the randomness of elastic modulus and system parameters on the first and second statistical moments The first- and second-order perturbation analysis schemes are presented in conjunction with the standard IGA method. Also, Monte Calo simulation is presented by combining the numerical generation of random field using the spectral representation method and standard isogeometric analysis. The validation of the proposed stochastic analysis schemes is determined by comparing the mean, variance, and response variability of the deflection at the center of plates predicted by the proposed method with those from MCS and previous studies. The numerical examples clearly demonstrate that the proposed scheme shows good agreement with those from MCS and the previous research results in the literature, and furthermore give more accurate results. The effect of second-order analysis shown is significant when the degree of uncertainty is high in the random field under consideration. In particular, not only the deterministic analysis results but also the stochastic responses are observed to be improved due to the adoption of the IGA. From the considering various parameters we observe the following: The effects of on the of contrast, the strongly . The effect of correlation ratio on the response COV obvious, and the response COV grows when the correlation distance ratio raises. random material properties provides the advanced analysis of plates. The results of the analysis of structures can use for advanced design.


INTRODUCTION
Functionally graded material is a type of composite material which material components vary continuously along the thickness. There are also similar functionally graded materials in nature such as bones and bamboo trees, which have natural functional grading. Typical functionally graded materials are made from ceramic and metal to improve heatproof properties and achieve high toughness (Koizumi 1997;Shiota and Miyamoto 1997;Miyamoto et al. 1999). Many research was conducted to investigate the behavior of beam, plate and shell made of functionally graded materials in the literature (Ha and Khue 2019;Hafed and Zenkour 2019;Phuc 2019;Vuong and Duc 2019) To compute structural behavior, structural systems are commonly idealized as a mathematical model in which system parameters such as the geometry of structures, material properties, and loading are related to the response of structures such as displacements, strains, and stresses. Researches of computational structural mechanics are divided into two major areas: deterministic analysis and nondeterministic analysis. In general, computational analyses in practice do not deal with the problems involving uncertainties in structural parameters. That is, they ignore the random heterogeneity of materials (e.g., soil, concrete, composites) and loading (e.g., vehicle, wind, wave) and favor using deterministic models with average or extreme values so that it leads to a crude representation of physical behavior with an unignorable error. The advanced structural models, in fact, need to encompass more detailed information on the involved system parameters. Accordingly, the analysis of structures needs to take into account uncertain parameters of structures and to estimate their effect on the responses.
Probabilistic computation on structures with uncertain parameters falls basically into two main categories: methods using a statistical approach and methods using a nonstatistical approach. The statistical approach involves the numerical generation of random parameters and statistical estimation of response due to the random parameters. This process is called Monte Carlo simulation (MCS). Among the nonstatistical approaches, most of the reports in the literature are associated with the stochastic finite element method (SFEM). There are a variety of SFEM schemes depending on how the uncertain parameters are dealt with. For example, Liu et al. (1986) suggested an SFEM called the "probabilistic finite element method" using the interpolation method for computing random fields. Stochastic finite element analysis (Takada 1990;Deodatis 1991;Wall and Deodatis 1994;Noh 2004;Hien 2020) has used the weighted integral of random fields on the finite element area as random variables; Spectral SFEM used Karhunen-Loève expansion series Ghanem and Spanos (1991) or homogeneous chaos expansion method Galal et al. (2008) for the representation of random fields.
Several studies have adopted the SFEM for structural analysis considering random system properties. The SFEM using a weighted integral approach Noh (2011) investigated the response variability of the bending of plate structures. Singh et al. (2008) studied the bending of a laminated composite plate considering uncertain system properties. Zhu and Wu (1991) developed a scheme for both distinct and repeated eigenvalue analysis using SFEM based on a local average technique. Graham and Deodatis (2001) used SFEM to investigate the variability of response of displacements and eigenvalues of structures with multi-uncertain parameters. Shang and Yun (2013) used ABAQUS combined with Karhunen-Loève expansion to simulate the stochastic response of structures under material uncertainties. Jun et al. (2014) analyzed the response of composite beam under stochastic excitations. Kaminski (2001) investigated the response variability of strain and stress tensors in the beam with random material and geometrical parameters. Cavdar et al. Cavdar et al. (2008) applied the perturbation based stochastic finite element method for predicting the response of three-dimension composite frames under earthquake forces. Real et al. (2017) proposed the new approach using Craig-Bampton method to model the uncertain dynamic systems.
Some researchers investigated the static, dynamics and buckling problems of FGM plates with random variables. Talha and Singh (2014) applied a perturbation technique to the finite element method for buckling of FGM plates considering uncertain material properties in thermal environments. Chakraborty and Rahman (2008) studied stochastic multiscale models for fracture analysis of functionally graded materials.
Recently, besides the SFEM, some researchers suggested a mesh-free method in investigating structures with uncertainties. For example, Rahman and Rao (2001) extended the element-free Galerkin method to develop a stochastic meshless method for the linear elasticity problem. Arun et al. (2010) combined Karhunen-Loève expansion with an element-free Galerkin method for the elastoplastic problem. Kim and Inoue (2004) developed the spectral stochastic element-free Galerkin method adopting Karhunen-Loève expansion and Polynomial Chaos series for the linear elastic problems involving random material properties. Hosseini and Shahabian (2014) investigated the stochastic elastic wave propagation in a thick hollow cylinder using a stochastic hybrid mesh-free method.
Latin American Journal of Solids and Structures, 2020, 17(7), e306 3/19 In recent years, isogeometric analysis (IGA) using nonuniform rational B-splines (NURBS) was proposed by Hughes et al. (2005). There have been many studies on the IGA for structural analyses (Bhardwaj et al. 2015b;Tran et al. 2015;Bui et al. 2016;Nguyen et al. 2018;Nguyen et al. 2019) and optimization problems (Ha et al. 2010;Taheri and Hassani 2014;Lieu and Lee 2019;Wang et al. 2019). Most of the previous studies have their focus on deterministic structural systems. For stochastic problems, limited studies on SIGA are available in the literature, for example, Bhardwaj et al. (2015a) used Monte Carlo simulation fatigue crack growth of FGMs using extended isogeometric analysis. Li et al. (2018a) proposed the spectral SIGA for linear elasticity problems. Nguyen et al. (2017) generated random fields by the spectral representation method to investigate the response variability of the buckling load of composite structures. Shahrokhabadi and Vahedifard (2018) combined isogeometric and generation random fields of water in the soil to model seepage in unsaturated soils. Zhang and Shibutani (2019) used polynomial chaos expansions to construct SIGA for uncertainty in shape. Ding et al. (2019) considered higher-order Taylor series of functions of random variables to propose the Isogeometric generalized nth order perturbation-based stochastic method. Eckert et al. (2020) developed a polynomial chaos method for an arbitrary random field in conjugate the standard isogeometric analysis for computational stochastic mechanics.
Almost all previous studies on relation stochastic isogeometric analysis involved random fields, discretization random fields used the Karhunen-Loève expansion or polynomial chaos method. The governing equations used these approach is quite complicated, so it needs to propose other approaches simpler than the approach mentioned above.
In the previous study, Hien and Noh (2017) dealt with the response variability of eigenvalues for the vibration problem of functionally graded plates with uncertain material properties. In this work, we focus on developing the SIGA for the static bending problem of functionally graded plates with uncertain elastic modulus. Based on the integration point method and perturbation technique to formulate the governing equations of SIGA, a strategic solution for the stochastic problem of static bending plate is presented. To demonstrate the appropriateness of the proposed scheme, numerical example analyses are performed using the proposed SIGA, and the results are compared with those given in the previous studies and with the results given by Monte Carlo simulation as well.

THE NURBS FUNCTIONS
In this section, the B-spline basis functions and nonuniform rational B-spline (NURBS) functions are briefly reviewed (Piegl and Tiller 1997;Cottrell et al. 2009 The B-spline basis functions are defined recursively starting with piecewise constants (p = 0,1, 2, …) as follows: The B-spline surfaces create by the tensor product of basis functions in two parametric dimensions of η and ξ with two knot vectors The non-uniform rational B-spline curve is constructed by combining the rational basis functions and coefficients at control points: In general, the NURBS surface of order p in the ξ direction and order q in the η direction can be expressed as:

ISOGEOMETRIC FORMULATIONS FOR FGM PLATES BASED ON REFINED PLATE THEORY
Consider a rectangular plate made of functionally graded material as shown in Figure 1. In this work, the refined plate theory (RPT) (Shimpi and Patel 2006) which enables taking into account the shear deformation effect is used. Thus, the displacements field at (x, y, z) in the FGM plates can be expressed as follows: where 0 0 , , , b s u v w w are the displacement components on the mid-plane.
By differentiating Eq. (7), the linear strains can be obtained as: x h z b a y Perturbation based stochastic isogeometric analysis for bending of functionally graded plates with the randomness of elastic modulus Hien Duy Ta et al.
Latin American Journal of Solids and Structures, 2020, 17 (7), e306 5/19 where: We consider an FGM plate made from a mixture of ceramic and metal with the subscripts m and c, respectively. The elastic modulus is assumed to be a random field in the plane of the plate: where ( ) , r x y is the homogeneous random field with zero mean and ( ) 0 E z is assumed as: where p is the power index of the volume fraction.
The auto-correlation function of the random fields ( ) The linear constitutive relations of FGM plates: ( ) Latin American Journal of Solids and Structures, 2020, 17 (7) The displacement field  U of the plate is approximated by the NURBS functions and the vector of nodal P U at the control point P: where: The strains can be computed as: The weak form for static for bending of FGM plates can be written as to form: where the stiffness matrices , b s A D are given as follows: Perturbation based stochastic isogeometric analysis for bending of functionally graded plates with the randomness of elastic modulus where , ,  = KU F (23) where K, U, and F denote stiffness matrix, generalized displacement, and force vector, respectively.
The global stiffness matrix K including the random field of elastic modulus is given by: and the load vector is given by:

STOCHASTIC ISOGEOMETRIC ANALYSIS FOR STATIC BENDING OF PLATES
Similar to numerical methods applied in the stochastic analyses such as SFEM and stochastic meshless method, the SIGA also necessitates the discretization of random fields into a vector of random variables. Various methods for the discretization of random fields have been developed in the literature, for example, the nodal point method (Liu, Belytschko et al. 1986;Michael Kleiber and Hien 1992), the integration point method (Brenner and Bucher 1995;Matthies et al. 1997), the local averaging method (Vanmarcke and Grigoriu 1983), the Karhunen-Loève expansion (Ghanem and Spanos 1991;Ngah and Young 2007;Chen and Guedes Soares 2008;Druesne et al. 2016;Li et al. 2018b) and the weighted integral approach (Takada 1990;Deodatis 1991;Noh 2005). In this work, the random field of elastic modulus is discretized following the integration point method, i.e., the random field is approximated at Gauss points of each element.
Computation of the stiffness matrix in Eq.(24) using Gaussian integration quadrature as form: where W Gauss is weights of the Gaussian quadrature. In Eq. (27), the number of Gauss points equals the number of random variables. The set of N random variables by discretization random field can be written in the form: We assume that all components of the random field r are small. Taking the Taylor series expansion at r = 0, i.e., at the mean of r, gives: In which: Perturbation based stochastic isogeometric analysis for bending of functionally graded plates with the randomness of elastic modulus For the first-order perturbation equation: For the second-order perturbation equation: sign the mean and covariance of the displacement U, respectively.
Using the expectation operator on Eq. (30), the solutions of the first-order perturbation technique for the displacement are: where ij R is the covariance matrix of a random variable r and N is the number of random variables.
The second-order perturbation solutions are: The generalized displacement vector, U, represents nodal parameters at control points; it is not the actual displacements at nodes. Using the relation between U and  U in Eq. (16), the mean vector and covariance matrix can be obtained as: where Θ is a transformation matrix which is calculated via the relation between U and  U in Eq.(16).
The response variability can be represented by means of the coefficient of variation (COV), which is defined as a ratio of the standard deviation to the mean of displacement as follows: The numerical solution procedure for a nonstatistical approach is shown on the flow chart in Figure 2. The main steps of the process solving can briefly be described as follows: 1) The first step is to input data and discretize the structural domain.
2) The second step is to calculate the deterministic part of the stiffness matrix and load vector.
3) The third step is to solve the zeroth perturbation equation in Eq. (32) to find out displacements, strains, stresses.

4)
The fourth step is to find the first-order perturbation solution for determining the mean and covariance matrix of the displacement parameter.

5)
The fifth step is to find the second-order perturbation solution for determining the mean and covariance matrix of the displacement parameter.

6)
The final step is to find the transformation matrix and then employ a loop to calculate the mean and covariance matrix of real displacements in Eq. (38).

MONTE CARLO SIMULATION
In verifying the proposed SIGA, the crude Monte Carlo simulation is adopted and performed. The spectral representation method reviewed by Shinozuka and Deodatis (Shinozuka and Deodatis 1991;Shinozuka and Deodatis 1996) is well suited for the numerical generation of random sample functions of the Gaussian zero-mean homogeneous random fields.
The Gaussian zero-mean homogeneous random fields ( ) , r x y can be represented by the summation expression for the random sample Shinozuka and Deodatis (1996): Here, ε is set to be 0.001 in the numerical generation.
With each random sample, the stiffness matrix is determined to contain the elastic modulus which is expressed as a cosine series as given in Eq.(40). The computation of the stiffness matrix is more complex due to the variable elastic modulus. It can approximately be calculated by the Gaussian quadrature for the integration of the stiffness matrix. Repeated solutions are needed for each sample to obtain the statistical properties of responses of the corresponding structure.

Numerical tests
The accuracy of the proposed analysis scheme can be verified if the method predicts well the previous reports given in the literature, in particular for the structures composed of isotropic materials. Even though the formulation is given for the FGM plates, it can also be applied to the analysis of the isotropic plates if we use corresponding parameters.
The first example is a square plate subject to a uniformly distributed load having a magnitude of 1.0 ( Figure 3). The data of this example follows the first example in the previous study by Choi and Noh (1996). The boundary conditions considered are two cases: simple support and clamped support. Because of the symmetry in the shape and the applied load of the structure, a quarter model is adopted (Figure 3). The mean Young's modulus is chosen to have the same value as used in the literature Choi and Noh (1996): E 0 = 10.29 × 10 3 . Poisson's ratio is 0.30. The thickness of the plate is assumed as 1.0. The auto-correlation function is assumed as: where, , x y ξ ξ are components of the separation vector ξ between two points in the domain of the structure. The corresponding correlation distances are denoted by 1 2 , d d . The coefficient of variation of the random field σ is assumed as 0.1.
We provide the stochastic response at the center point of the structure (point P in Figure 3). In the numerical analyses, we use 64 quadratic NURBS elements (8 × 8 mesh).
In Table 1, we compare the mean and standard deviation of center displacements of the plate obtained by three analyses: present approach -the first-order solution of the SIGA scheme, Monte Carlo simulation with 10,000 samples based on the IGA scheme, weighted integral method, and Monte Carlo simulation Choi and Noh (1996). It is observed that all the results of stochastic isogeometric analyses, irrespective of boundary conditions, show larger values than those obtained by the weighted integral method. It is noted that the MCS results are even larger for SIGA than for FEM, which shows the superiority of the IGA scheme over conventional finite element analysis. Table 1 Comparison of SIGA results with Noh Choi and Noh (1996) We consider a normalized displacement at the center of the plate: where the stiffness parameter:  . It is clear that all the responses of the second-order SIGA are larger than those of the first-order SIGA. It is also observed that the mean and the standard deviation of normalized displacement increases, with a decreasing rate, as the power index p increases. However, we can note that the response COV is not affected by the order of the scheme adopted (first or second order), as shown in Figure 8, because the denominator (mean) and the numerator (standard deviation) of the response COV increase simultaneously with almost the same degree. It is clear that the effect of the power index p on response COV in Fig. 8b is very small. A similar trend is observed also for the clamped plate. The power index p affects the mean and the standard deviation of the normalized displacement as shown in Figs. 6 and 7, however, it does not affect the responses COV.

The effect of correlation distance ratio on the static bending of plates
The effect of correlation distance ratio between two directions, 2 1 d d , on the response variability of displacements is considered in this part. We provide results obtained by using the first-order SIGA for the same plate used in Section 5.2.
The dimensions of FGM plates: a = 0.2(m), b = 0.3(m), h = 0.01(m). The power index p is fixed as 1.0. The effect of the correlation distance ratio 2 1 d d on the response COV of the normalized displacement at the center of FGM plates with simple and fixed boundary conditions is represented in Figure 9 with the coefficient 0 1 . σ = . It is obvious that the response COV goes up as the ratio 2 1 d d increases. SIGA is proposed for the static bending analysis of FGM plates with randomness in elastic modulus. In constructing the governing equations for the FGM plate, the effect of higher-order shear deformation is taken into account. Two dimensional random fields of elastic modulus are discretized using the integration point method to formulate the firstand second-order perturbation expansion of displacements and stiffness matrix. This approach is simpler than the Karhunen-Loève expansion approach or polynomial chaos approach in the Spectral stochastic isogeometric analysis.
The first-and second-order perturbation analysis schemes are presented in conjunction with the standard IGA method. Also, Monte Calo simulation is presented by combining the numerical generation of random field using the spectral representation method and standard isogeometric analysis. The validation of the proposed stochastic analysis schemes is determined by comparing the mean, variance, and response variability of the deflection at the center of plates predicted by the proposed method with those from MCS and previous studies. The numerical examples clearly demonstrate that the proposed scheme shows good agreement with those from MCS and the previous research results in the literature, and furthermore give more accurate results. The effect of second-order analysis shown is significant when the degree of uncertainty is high in the random field under consideration. In particular, not only the deterministic analysis results but also the stochastic responses are observed to be improved due to the adoption of the IGA.
From the analyses considering various parameters of the given FGM plate problems, we observe the following: The effects of power index on the response COV of displacement is negligible. In contrast, the influence of power index p on the standard deviation of displacement is prominent depending strongly on the value of p. The effect of correlation distance ratio on the response COV is obvious, and the response COV grows when the correlation distance ratio raises.
In this work, the SIGA which calculates the stochastic response considering the random field of material properties provides the advanced stochastic analysis of plates. The results of the stochastic analysis of structures can use for advanced design.