Probabilistic stability of uncertain composite plates and stochastic irregularity in their buckling mode shapes: A semi-analytical non-intrusive approach

In this study, the mechanical properties of the composite plate were considered Gaussian random fields and their effects on the buckling load and corresponding mode shapes were studied by developing a semi-analytical non-intrusive approach. The random fields were decomposed by the Karhunen-Loève method. The strains were defined based on the assumptions of the first-order and higher-order shear-deformation theories. Stochastic equations of motion were extracted using Euler–Lagrange equations. The probabilistic response space was obtained by employing the non-intrusive polynomial chaos method. Finally, the effect of spatially varying stochastic properties on the critical load of the plate and the irregularity of buckling mode shapes and their sequences were studied for the first time. Our findings showed that different shear deformation plate theories could significantly influence the reliability of thicker plates under compressive loading. It is suggested that a linear relationship exists between the mechanical properties’ variation coefficient and critical loads’ variation coefficient. Also, in modeling the plate properties as random fields, a significant stochastic irregularity is obtained in buckling mode shapes, which is crucial in practical applications.


Introduction
Due to their outstanding mechanical and physical properties, composite materials have been increasingly applied in various industry sectors. The high specific strength, high stiffness, excellent fatigue behavior, and corrosion resistance of composite materials have increased their utilization in weight-critical applications, such as aerospace structures [1,2]. Specific fabrication techniques (layering, curing process, etc.) and the uncertainty in service life caused by environmental conditions have led to a considerable statistical dispersion in the properties of composite materials compared to their metallic counterparts [3,4]. The wide scattering observed in the mechanical properties and geometrical parameters (such as fiber angle and ply thickness), as well as the exposure of composite structures to scattered thermal and mechanical loads have posed a fundamental and unavoidable challenge to the analysis of their structural reliability [5]. Recently, attention to the stochastic nature of engineering structures and associated uncertainty assessment has become an essential topic in the research community.
Several methods regarding the uncertainty assessment and quantification of composite materials, including the Finite Element Method (FEM) [6], Monte Carlo (MC) [7], Perturbation method [8,9], and Polynomial Chaos (PC) [10], have been proposed and applied in numerous engineering problems. These methods can generally be divided into sampling-and non-sampling-based approaches. In sampling-based methods, the problem must be run from the beginning for each new input value [11,12].
For example, the MC method is a straightforward sampling-based approach for low-scale models [7]. On the other hand, in non-sampling-based techniques, such as the PC method, the problem is solved for a limited number of inputs, and a polynomial is subsequently fitted to the probabilistic response space. The obtained PC expansion is run for each new value rather than the entire model, significantly reducing the computational cost compared to sampling-based approaches [13]. The advantage of non-sampling methods is noteworthy, especially for large-scale mathematical models where sampling-based strategies do not have the appropriate computational performance.
Several studies have investigated the buckling of composite plates, in which the experimental buckling load is generally different from the predicted value. The primary justification for this discrepancy is related to the imperfections observed in the geometry and properties of the modeled structures [14][15][16][17]. Although geometrical uncertainty is not usually discussed for flat plates, scattered mechanical properties have been frequently observed in experimental tests and identified as an effective critical parameter in buckling [18][19][20]. Nguyen et al. [21] reported non-uniform mechanical properties across a composite plate spatially varying with a completely stochastic nature that affects the buckling loading. Despite these observations, the buckling phenomenon is commonly simulated, assuming uniform and deterministic mechanical properties throughout the plate [22]. Few studies have considered the stochastic nature of composites, focusing on the mechanical behavior (linear and nonlinear), boundary conditions, or the uncertainty modeling type [23][24][25][26][27].
Experimental observations, however, indicate a probabilistic spatial distribution of geometrical and mechanical properties in any composite component. Therefore, this study considers the elastic properties of composite plates based on the mean values with spatial uncertainty and randomly oriented imperfections to investigate their effects on the buckling loads and corresponding mode shapes. Studying the stochastic bulking mode shapes is crucial for practical buckling tests and is addressed for the first time in the present work. For this purpose, all plate properties are modeled as a Gaussian random fields using the exponential kernel, which is commonly used in reliability analysis and validated using experimental results [47,48]. The probabilistic response space was modeled using the PC method to investigate and quantify uncertainty in buckling load and corresponding mode shapes. Applying the assumed stochastic modes (SAM) and Polynomial Chaos method introduces a new semianalytical, non-intrusive approach, which is computationally efficient with a reduced size of the mathematical model for composite plates under compressive loading.

Problem statement
The present study analyzes the buckling behavior of a square composite plate ( Fig. 1) with non-deterministic properties under compression loading. The mechanical properties, modeled as Gaussian random fields (not random variables), consist of orthotropic tensile and shear moduli, as identified in Ref. [49]. The sources of uncertainty in the mechanical properties include spatially variable fiber volume fraction and fiber orientation, along with the randomness from non-ideal cure processes. An exponential autocorrelation function was considered for the random fields. Note that the significant and minor Poisson's ratios were assumed to be constant based on their minimal effect on buckling behavior. The stochastic field discretization was performed by utilizing the Karhunen-Loève (K-L) method theory. In addition, the stochastic equations of motion were derived by employing the Euler-Lagrange method using the SAM and the assumption of First Order or Higher Order Shear Deformation Theories, FSDT and HSDT, respectively. The statistical properties of the critical load were obtained using the elastic and geometric stiffness matrix and through the non-intrusive polynomial chaos (Pc) approach.
Two boundary conditions were considered, including all simply supported or clamped edges. The uncertainty propagation in critical loads was studied, while the plate was subjected to a uniaxial or biaxial compressive load by employing different plate deformation theories. Furthermore, the effects of uncertainty propagation, corresponding to the plate's properties, on the uncertainty of critical loads and buckling mode shapes were assessed.

Constitutive equations
First, stochastic properties were defined. The random field's autocovariance intended for spatially variable and stochastic properties was considered with an exponential kernel, which is commonly used in structural reliability problems [49]. The exponential kernels corresponding to random properties were extracted through Eqs. (1)-(3) [49]. Three autocovariance were required according to the assumptions of Eq. (4).
where C is kernel intended for auto-correlation function; and are correlation length in x and y direction, respectively; is the standard deviation of the desired properties.
For defining autocovariance related to stochastic properties, the K−L expansion was used to discretize the continuous random field of stochastic properties. For example, the random field of the tensile modulus in the fiber direction and on the ranges defined for x and y in Eq. (5) is obtained as Eq. (6) through K−L expansion and based on the product of eigenvalues and eigenfunctions, as well as standard random variables [50].
Equations (7)-(9) were used to obtain eigenvalues ( ) and eigenvectors ( ) in the K-L expansion for the problem with exponential autocovariance related to Eqs.
(1)-(3), depending on the correlation length of the autocovariance function in different directions. This process has been described in detail in Ref. [50]. τ The number of terms in the K-L expansion (n) was determined using the criteria of Eq. (10), where is usually considered to be 90% [51]. Based on Eq. (10), a shorter correlation length increases the number of K-L expansion terms in Eq. (6) to satisfy this criterion.
The displacement field was defined to derive the stochastic motion equations, which were performed through the FSDT and HSDT (Eq. (11)) using the displacement of the middle points of the plate in different directions.
u(x, y, where , , and are the displacements of the plate midpoints in the directions of x, y, z, respectively, as shown in Fig. 1. Ω is a stochastic vector of standard random variables. Through linear definitions of strain and assuming small displacements, linear strains can be extracted according to Eqs. (12)-(18) [52], where ε and γ are in-plane and out-of-plane strains. Superscripted index zero was defined for the strains caused by the displacement of the middle plate, and superscripted indices 1 and 3 are for the strains caused by rotation. { where and are constants that are equal to zero for FSDT and and for HSDT. Assuming that the structure has a linear elastic behavior, the strain energy caused by the elastic stresses is defined using Eq. (19) through the sum of the potential energies of the different layers of the plate.
is for the number of layers. The shear correction factor ( ) for FSDT and HSDT is 5/6 and 1, respectively. Potential energy caused by work of external compressive force ( ) was extracted according to Eq. (20), where and compressive loads are as shown in Fig. 1 and was defined as Eq. (21). The on-axis spatially variable and stochastic stiffness matrix elements were defined based on Eq. (22). Also, the constitutive equation of a composite layer (layer k) can be defined as Eqs. (23) and (24), where it is necessary to extract the stiffness matrix under rotation for each layer depending on the angle of fiber to the principal coordinate direction. , y, Ω),

Motion equations q(θ)
Based on the definitions of strain energy, the principle of minimum potential energy was employed to obtain motion equations of the structure under compressive loading (Eq. (25)), where is the stochastic degree of freedom in the global coordinate. The SAM, commonly used in dynamic problems, was employed to discretize the equations into generalized coordinates using Eq. (25) [47]. As beam functions, suitable admissible modes for simply supported and clamped boundary conditions are respectively given in Eqs. (26) and (27).
, and are generalized coordinates and Ω is the probabilistic vector of standard random variables.
K gK ξ where and are the average of the geometric and elastic stiffness matrices, respectively.
is the probabilistic standard variable corresponding to each stochastic property in the expansion of Eq. (6). Equation (28) is also defined for predefined stochastic properties with arbitrary (predefined) standard deviation and can easily be converted into an eigenvalue problem. In addition, the statistical properties of the eigenvalue (i.e., the critical load of the plate) can be extracted by generating random numbers for random variables based on the sampling-based approach (MC method) or the nonsampling-based approach (Pc approach). Here, the nonintrusive Pc approach was used to obtain the probability space of critical buckling loads.

Non-intrusive polynomial chaos approach
PC expansion is used for defining random fields and random variables. The second-order variable was considered, which can be described as possibly infinite variables of , , , … The independent variables were gathered into and the variables, including stochastic eigenvalues and eigenvectors, can be obtained in the form of , . Equation (29) was utilized for defining stochastic eigenvectors related to Eq. (28), where the set of orthonormal polynomials is defined as , which is a principle for . Also, deterministic coefficients introduce X vectors. It can be obtained as Eqs. (29) and (30), if based on orthonormal characteristic . Consequently, the variance of X is defined as Eq. (31). A similar definition is used for the definition of stochastic eigenvalues.
Since all numerical calculations require parameterization of the random variable X with a limited number of of independent random variables, X can be defined as Eq. (32) if . and are the result of random modeling and a compromise between accuracy and efficiency, respectively.
where . The PC polynomial sequence is selected so that low-order approximations remain constant to increase PC values. The expansion was extended to a second-order equation and a set of random variables was placed in interval of . u is a parameterizing variable that can vary by changing the PC expansion's definite coefficients (Eqs. (33) and (34)). The need to change the model is eliminated using the non-intrusive method and by considering the desired model as a black box. Therefore, the estimated expansion coefficients were achieved using a sparse-grid-based method. Here, a form of the nonintrusive method with linear regression was used so that random variables were considered as a finite number of and accordingly, the coefficients were obtained using . Then, the final equation was extracted as Eq. (35), which can be applied to find the basis ( ) of stochastic space of eigenvalue (critical load) and corresponding eigenvectors (buckling mode shapes) in design points instead of a full run based on the MC simulation.

Verification of equations
To verify the equations and codes, the plate buckling problem was solved in two scenarios, with deterministic and stochastic properties. Then, the obtained results were compared with previous studies.

Deterministic properties
The square composite with simple supported stacking sequence and lamina properties, according to Table 1 was analyzed under uniaxial loading. The normalized buckling critical load was defined according to . The comparison of our results with previous studies showed a high level of agreement, verifying our equations and codes for plates with deterministic properties (Table 2).

Stochastic properties
The square composite plate with simply supported boundary condition, stacking sequence [0,90], and lamina properties ( Table 1) was analyzed so that the modulus of elasticity in the direction of the fiber was considered as a random variable. The obtained results for different coefficients of variation were compared with previous studies (Fig. 2). A 5th order PC was used in non-intrusive simulations. As shown in Fig. 1, the accuracy of the present method was optimal and better than the methods defined based on the perturbation method. In addition, our results showed a high level of agreement with previous studies defined based on the MC method.

Effects of different shear deformation plate theories on the uncertainty propagation of critical loads
In the first step, the uncertainty propagation in the critical buckling load was studied by applying different shear property quantity deformation plate theories in uniaxial and biaxial loading scenarios for two plates with symmetric [0,90,0] and antisymmetric [45,−45] stacking sequences. The properties of the plate, including the modulus in the direction of the fiber and perpendicular to the fiber, as well as the shear modulus, were considered Gaussian random fields. The sides to thickness ratio was considered variable to evaluate the critical load sensitivity of the induction of shear stresses in the cross-section of the plate. Figure 3 shows the PDF of the critical buckling load in different loading scenarios and the ratios of different sides to thickness for the plate with all edges simply supported as boundary conditions. As shown in Fig. 3, the probability distribution of the critical load was Gaussian. As the thickness of the plate increased, the sensitivity of the effect of shear stresses on the critical load increased. Also, a significant difference was observed between the results of the FSDT and HSDT theories, in particular at higher thicknesses. In higher thickness plates, the HSDT theory is preferred over the FSDT for simulating the shear flow in the plate section due to the higher accuracy. It is due to shear locking and inaccurate approximation of shear loads that FSDT makes. Overall, the critical load of antisymmetric plates is more sensitive to shear load accuracy due to asymmetric shear distribution, which is why the HSDT theory was used for analysis in the present study. In addition, the observed increase in the critical buckling load accompanied by a decrease in thickness is due to the process of buckling load normalization.

Effects of boundary conditions on the uncertainty propagation in the plate buckling loads
The symmetrical and antisymmetric plates defined in the previous section were examined under different boundary conditions. The graphs corresponding to the changes in the normalized critical load standard deviation were compared to the increase in the random field variation coefficient under uniaxial and biaxial loads, as shown in Fig. 4. Our results showed that the standard deviation sensitivity of critical load to the uncertainty propagation in the properties with clamped boundary conditions was higher than the standard deviation sensitivity of plates with simply supported boundary conditions. Figure 5 illustrates the results of variation coefficient sensitivity of critical buckling load to variation coefficient of the properties in the Gaussian random field for different thicknesses, different loading scenarios, and different boundary conditions. As seen in Fig. 5, the critical load variation coefficient in uniaxial and biaxial loading conditions, as well as different boundary conditions and   different thicknesses resulted in a similar, yet unique graph. This graph can be used as the main graph or reference graph to calculate the buckling reliability of plates under compressive loading.

Investigating the buckling mode shapes
Due to the composite plates' spatially varying properties and the random variations, it is of special importance to study the buckling mode shapes in practical buckling tests. For example, the arrangement of strain gauges in plate buckling tests varies depending on the expected mode shape. Therefore, the deterministic mode shape and four stochastic samples of mode shapes for the first three buckling modes are shown in Table 3 for symmetric composite plates with different boundary conditions, regardless of their occurrence likelihood and with assuming the Gaussian random properties introduced in the previous sections.
As seen in Table 3, the occurrence sequence of the buckling modes was completely irregular and stochastic, which has commonly been observed in practical tests that are well developed as buckling of imperfect plates. This can be due to the creation of areas with variable and random stiffness in the plate that affect the buckling mode shapes. The PDF of higher-order modes was studied to investigate the phenomenon and the occurrence likelihood of these modes and quantify the uncertainty in the sequence of buckling modes. Therefore, the PDF corresponding to the three first critical buckling loads for the plate with symmetric stacking sequence and a/h = 50 under uniaxial loading was extracted. As shown in Fig. 6, there was a chance of probabilistic event interference for modes 2 and 1, as well as for modes 2 and 3. In addition, the occurrence likelihood of the buckling mode shape corresponding to mode 2 is greater than mode 1, which is lower for modes 1 and 3. Therefore, the sequence of buckling modes is completely random according to spatially variable and stochastic properties (Fig. 6). Moreover, it is necessary to consider the non-uniform distribution of composite properties in the plate coordinates as the statistical scattering is strongly observed in the critical load value and in the corresponding buckling mode shapes. This can be crucial in the buckling tests of composite plates, where it cannot be predicted by ignoring the spatial variable of properties in the plate coordinates. With that said, simulating the properties of composites as a random field and not a random variable is crucial in calculating the reliability of composite plates under buckling failure modes.

Conclusions
This work investigated the buckling of composite plates with symmetrical and antisymmetric stacking sequences, different boundary conditions, and different loading scenarios. The elastic properties of the plate were considered as a Gaussian random field. The motion equation was obtained using FSDT and HSDT theories, as well as energy equations. Then, the buckling eigenvalue problem was extracted by employing the stochastic assume mode method. The equations in random space were solved using the non-intrusive Pc approach.
Our results showed that the high-ordered shear deformation theory offers better accuracy in probabilistic space, especially for thicker plates. This is a key issue in calculating the reliability of plates under buckling loading because the PDF tends to shift to larger values, as commonly observed in plates with anti-symmetrical stacking sequences.
Furthermore, the uncertainty propagation in the critical load was investigated for plates with simply supported and clamped boundary conditions. The results showed that the standard deviation of critical load variation for plates with clamped supported boundary conditions is greater than those of similar plates with simply supported boundary conditions. The stability of the obtained graphs corresponding to the critical load variation coefficient relative to the variation coefficient of properties indicates that these graphs can be used to calculate the reliability of similar plates under buckling loading. Finally, the irregularity of the mode shapes and the sequence of the buckling modes were discussed. The geometrical irregularity and stochastic sequences observed in the modes were justified by applying PDF for higher-order modes. The spatially varying stochastic properties significantly affect the critical buckling load and should be considered in practical buckling tests due to the effect on the plate reliability in buckling failure mode.