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

Abstract

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.

Graphical abstract

Keywords

uncertain composite plate / stochastic assume mode method / Karhunen−Loève theorem / polynomial chaos approach / plate buckling / irregularity in buckling mode shapes

1 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–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–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–27].

Previous studies have considered the effect of random but uniform fluctuations of mechanical properties on the buckling behavior of composite plates [19,20,28,29]. Recent developments in neural network-based approaches, manifold methods, and NURBS-based methods enhance the performance of classic approaches in plate stability and dynamics [30–36]. Different practical stochastic frameworks and sensitivity analyses can be found in literature and applied in stochastic polymer composites simulations [37–46].

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 semi-analytical, non-intrusive approach, which is computationally efficient with a reduced size of the mathematical model for composite plates under compressive loading.

2 Materials and methods

2.1 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.

2.2 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).

(1)C_{{\text{E}}_{\text{11}}}={\sigma }_{{\text{E}}_{\text{11}}}^2\exp (-|(x-x^{\prime })/l_x|-|(y-y^{\prime })/l_y|),

(2)C_{{\text{G}}_{\text{12}}}={\sigma }_{{\text{G}}_{\text{12}}}^2\exp (-|(x-x^{\prime })/l_x|-|(y-y^{\prime })/l_y|),

(3)C_{{\text{G}}_{\text{23}}}={\sigma }_{{\text{G}}_{\text{23}}}^2\exp (-|(x-x^{\prime })/l_x|-|(y-y^{\prime })/l_y|),

(4)\frac{{\nu }_{12}{E}_2}{1-{\nu }_{12}{\nu }_{21}}=\frac{{\nu }_{21}{E}_1}{1-{\nu }_{12}{\nu }_{21}},G_{12}=G_{13}.

where C is kernel intended for auto-correlation function; l_x and l_y are correlation length in x and y direction, respectively; \sigma 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].

(5)x\in \left[-\frac{a}{2},\frac{a}{2}\right],y\in \left[-\frac{a}{2},\frac{a}{2}\right],

(6)E_{11}(x,y,\mathrm{\Omega })={\overline{E}}_{11}+{\sigma }_{{\text{E}}_{\text{11}}}\sum _{i=1}^n\sqrt{{\lambda }_{x_i}{\lambda }_{y_i}}{\xi }_i{\varphi }_i(x){\varphi }_i(y).

Equations (7)–(9) were used to obtain eigenvalues (\lambda) and eigenvectors ({\varphi }_i(x)) 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].

(7)\{\begin{array}{c}c-{\omega }_i\tan ({\omega }_{i}d)=0,i:\text{odd},\\ {\omega }_i+c\tan ({\omega }_{i}d)=0,i:\text{even},\end{array}\begin{array}{c}c=1/l_x,1/l_y,\\ d=a,\end{array}

(8){\lambda }_{\mathit{\text{z}}}=\frac{2c}{{\omega }_{\mathit{\text{z}}}^2+c^2},{\lambda }_{\mathit{\text{z}}}={\lambda }_{x_i},{\lambda }_{y_i},

(9){\varphi }_i(v)=\{\begin{array}{c}{\displaystyle \frac{\cos ({\omega }_{i}v)}{\sqrt{a+\sin (2{\omega }_{i}a)/2{\omega }_i}}},\text{for}i:\text{odd},\\ {\displaystyle \frac{\sin ({\omega }_{i}v)}{\sqrt{a-\sin (2{\omega }_{i}a)/2{\omega }_i}}},\text{for}i:\text{even},\end{array}v=x,y.

The number of terms in the K–L expansion (n) was determined using the criteria of Eq. (10), where \tau 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.

(10)\frac{{\displaystyle \sum _{i=1}^n{\lambda }_i}}{{\displaystyle \sum _{i=1}^{\mathrm{\infty }}{\lambda }_i}}⩾\tau .

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.

(11)\begin{array}{c}u(x,y,\Omega )=u_0+\mathit{\text{z}}{\varphi }_x-c_1{\mathit{\text{z}}}^3({\varphi }_x+w_{0,x}),\\ v(x,y,\Omega )=v_0+\mathit{\text{z}}{\varphi }_y-c_1{\mathit{\text{z}}}^3({\varphi }_y+w_{0,y}),\\ w(x,y,\Omega )=w_0,\end{array}

where u_0, v_0, and w_0 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.

(12)\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{xy}\end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_{xx}^{(0)}\\ {\epsilon }_{yy}^{(0)}\\ {\epsilon }_{xy}^{(0)}\end{array}\right\}+\mathit{\text{z}}\left\{\begin{array}{c}{\epsilon }_{xx}^{(1)}\\ {\epsilon }_{yy}^{(1)}\\ {\epsilon }_{xy}^{(1)}\end{array}\right\}+{\mathit{\text{z}}}^3\left\{\begin{array}{c}{\epsilon }_{xx}^{(3)}\\ {\epsilon }_{yy}^{(3)}\\ {\epsilon }_{xy}^{(3)}\end{array}\right\},

(13)\left\{\begin{array}{c}{\gamma }_{y\mathit{\text{z}}}\\ {\gamma }_{x\mathit{\text{z}}}\end{array}\right\}=\left\{\begin{array}{c}{\gamma }_{y\mathit{\text{z}}}^{(0)}\\ {\gamma }_{x\mathit{\text{z}}}^{(0)}\end{array}\right\}+{\mathit{\text{z}}}^2\left\{\begin{array}{c}{\gamma }_{y\mathit{\text{z}}}^{(2)}\\ {\gamma }_{x\mathit{\text{z}}}^{(2)}\end{array}\right\},

(14)\left\{\begin{array}{ccc}{\epsilon }_{xx}^{(0)}& {\epsilon }_{yy}^{(0)}& {\epsilon }_{xy}^{(0)}\end{array}\right\}=\left\{\begin{array}{ccc}{u}_{0,x}& v_{0,y}& u_{0,y}+v_{0,x}\end{array}\right\},

(15)\left\{\begin{array}{ccc}{\epsilon }_{xx}^{(1)}& {\epsilon }_{yy}^{(1)}& {\epsilon }_{xy}^{(1)}\end{array}\right\}=\left\{\begin{array}{ccc}{\varphi }_{x,x}& {\varphi }_{y,y}& {\varphi }_{x,y}+{\varphi }_{y,x}\end{array}\right\},

(16)\begin{array}{rl}& \left\{\begin{array}{ccc}{\epsilon }_{xx}^{(3)}& {\epsilon }_{yy}^{(3)}& {\epsilon }_{xy}^{(3)}\end{array}\right\}\\ =& -c_1\left\{\begin{array}{ccc}{\varphi }_{x,x}+w_{0,xx}& {\varphi }_{y,y}+w_{0,yy}& {\varphi }_{x,y}+{\varphi }_{y,x}+2w_{0,xy}\end{array}\right\},\end{array}

(17)\left\{\begin{array}{cc}{\gamma }_{x\mathit{\text{z}}}^{(0)}& {\gamma }_{y\mathit{\text{z}}}^{(0)}\end{array}\right\}=\left\{\begin{array}{cc}{\omega }_{0,x}+{\phi }_x& {\omega }_{0,y}+{\phi }_y\end{array}\right\},

(18)\left\{\begin{array}{cc}{\gamma }_{y\mathit{\text{z}}}^{(2)}& {\gamma }_{x\mathit{\text{z}}}^{(2)}\end{array}\right\}=-c_2\left\{\begin{array}{cc}{w}_{0,x}+{\varphi }_x& w_{0,y}+{\varphi }_y\end{array}\right\},

where c_1 and c_2 are constants that are equal to zero for FSDT and c_1=4/(3h^2) and c_2=3c_1 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. N.layer is for the number of layers. The shear correction factor (K_{\text{s}}) for FSDT and HSDT is 5/6 and 1, respectively. Potential energy caused by work of external compressive force (U_{\text{f}}) was extracted according to Eq. (20), where N_{\text{c}x} and N_{\text{c}y} compressive loads are as shown in Fig.1 and d_{ij} 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.

(19)\begin{array}{rl}{U}_{\text{p}}=\sum _{i=1}^{N.layer}\frac{1}{2}\underset{V}{\int }& ({\sigma }_{xx}{\epsilon }_{xx}+{\sigma }_{yy}{\epsilon }_{yy}+{\sigma }_{xy}{\epsilon }_{xy}\\ & +K_{\text{s}}{\sigma }_{y\mathit{\text{z}}}{\epsilon }_{y\mathit{\text{z}}}+K_{\text{s}}{\sigma }_{x\mathit{\text{z}}}{\epsilon }_{x\mathit{\text{z}}}{)}_i\text{d}V,\end{array}

(20)U_{\text{f}}=\frac{1}{2}\underset{A}{\int }(N_{\text{c}x}{d}_{xx}+N_{\text{c}y}{d}_{yy})\text{d}x\text{d}y,

(21)d_{ij}=u_{,i}{u}_{,j}+v_{,i}{v}_{,j}+w_{,i}{w}_{,j},

(22)\begin{array}{rl}& Q_{11}=\frac{E_{11}(x,y,\mathrm{\Omega })}{1-{\nu }_{12}{\nu }_{21}},Q_{12}=\frac{{\nu }_{12}{E}_{22}(x,y,\mathrm{\Omega })}{1-{\nu }_{12}{\nu }_{21}}=\frac{{\nu }_{21}{E}_{11}(x,y,\mathrm{\Omega })}{1-{\nu }_{12}{\nu }_{21}},\\ & Q_{22}=\frac{E_{22}(x,y,\mathrm{\Omega })}{1-{\nu }_{12}{\nu }_{21}},Q_{66}=G_{12}(x,y,\mathrm{\Omega }),\\ & Q_{44}=G_{23}(x,y,\mathrm{\Omega }),Q_{55}=G_{13}(x,y,\mathrm{\Omega }),\end{array}

(23){\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}}_m^k={\left[\begin{array}{ccc}{\overline{Q}}_{11}& {\overline{Q}}_{12}& {\overline{Q}}_{16}\\ {\overline{Q}}_{21}& {\overline{Q}}_{22}& {\overline{Q}}_{26}\\ {\overline{Q}}_{16}& {\overline{Q}}_{26}& {\overline{Q}}_{66}\end{array}\right]}^k\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{xy}\end{array}\right\},

(24){\left\{\begin{array}{c}{\sigma }_{y\mathit{\text{z}}}\\ {\sigma }_{x\mathit{\text{z}}}\end{array}\right\}}_m^k={\left[\begin{array}{cc}{\overline{Q}}_{44}& {\overline{Q}}_{45}\\ {\overline{Q}}_{45}& {\overline{Q}}_{55}\end{array}\right]}^k\left\{\begin{array}{c}{{\gamma }^0}_{y\mathit{\text{z}}}\\ {{\gamma }^0}_{x\mathit{\text{z}}}\end{array}\right\}.

2.3 Motion equations

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 q(\theta ) 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).

(25)\begin{array}{c}\pi =(U_{\text{P}}+U_{\text{f}}),\\ \frac{\mathrm{\partial }\pi }{\mathrm{\partial }{q}_i}(q(\mathrm{\Omega }))=0,i=1,...,n\end{array}

(26)\begin{array}{c}{u}_0(\mathrm{\Omega },x,y)=\sum _{m=1}^M\sum _{n=1}^N{U}_{mn}(\mathrm{\Omega })\sin \left(\frac{n\text{π}}{l}x\right)\cos \left(\frac{m\text{π}}{l}y\right),\\ v_0(\mathrm{\Omega },x,y)=\sum _{m=1}^M\sum _{n=1}^N{V}_{mn}(\mathrm{\Omega })\cos \left(\frac{n\text{π}}{l}x\right)\sin \left(\frac{m\text{π}}{l}y\right),\\ w_0(\mathrm{\Omega },x,y)=\sum _{m=1}^M\sum _{n=1}^N{W}_{mn}(\mathrm{\Omega })\sin \left(\frac{n\text{π}}{l}x\right)\sin \left(\frac{m\text{π}}{l}y\right),\\ {\varphi }_x(\mathrm{\Omega },x,y)=\sum _{m=1}^M\sum _{n=1}^N{X}_{mn}(\mathrm{\Omega })\cos \left(\frac{n\text{π}}{l}x\right)\sin \left(\frac{m\text{π}}{l}y\right),\\ {\varphi }_y(\mathrm{\Omega },x,y)=\sum _{m=1}^M\sum _{n=1}^N{Y}_{mn}(\mathrm{\Omega })\sin \left(\frac{n\text{π}}{l}x\right)\cos \left(\frac{m\text{π}}{l}y\right),\end{array}

(27)\begin{array}{c}{\mathrm{\Gamma }}_m(x)=\left(\cosh {\displaystyle \frac{{\lambda }_{m}x}{a}}-\cos {\displaystyle \frac{{\lambda }_{m}x}{a}}\right)-{\gamma }_i\left(\sinh {\displaystyle \frac{{\lambda }_{m}x}{a}}-\sin {\displaystyle \frac{{\lambda }_{m}x}{a}}\right),\\ {\mathrm{\Pi }}_m(x)={\displaystyle \frac{\mathrm{\Gamma }(x)}{{\lambda }_m/a}},\\ u_0(\mathrm{\Omega },x,y)={\displaystyle \sum _{m=1}^M{\displaystyle \sum _{n=1}^N{U}_{mn}(\mathrm{\Omega }){\mathrm{\Gamma }}_m(x){\mathrm{\Pi }}_n(x)},}\\ v_0(\mathrm{\Omega },x,y)={\displaystyle \sum _{m=1}^M{\displaystyle \sum _{n=1}^N{V}_{mn}(\mathrm{\Omega }){\mathrm{\Gamma }}_m(x){\mathrm{\Pi }}_n(x)},}\\ w_0(\mathrm{\Omega },x,y)={\displaystyle \sum _{m=1}^M{\displaystyle \sum _{n=1}^N{W}_{mn}(\mathrm{\Omega }){\mathrm{\Gamma }}_m(x){\mathrm{\Pi }}_n(x)},}\\ {\phi }_x(\mathrm{\Omega },x,y)={\displaystyle \sum _{m=1}^M{\displaystyle \sum _{n=1}^N{X}_{mn}(\mathrm{\Omega }){\mathrm{\Gamma }}_m(x){\mathrm{\Pi }}_n(x)},}\end{array}

where U_{mn}, V_{mn}, W_{mn}, X_{mn}, and Y_{mn} are generalized coordinates and Ω is the probabilistic vector of standard random variables.

Using Eq. (25) and by placing stochastic properties (Eq. (22)) (after rotation, for rotated layers) and the assumed displacement field (Eq. (11)) in the energy Eqs. (19) and (20), the motion Equations can be extracted as Eq. (28).

(28)\left\{{\overline{K}}_{\text{e}}+\sum _{i=1}^n{K}_{e_i}{\xi }_i+({N_{\text{c}}}_x+N_{\text{c}y})\left({\overline{K}}_g+\sum _{i=1}^n{K}_{g_i}{\xi }_i\right)\right\}\left\{X\right\}(\mathrm{\Omega })=0,

where {\overline{K}}_{\text{g}} and \overline{K} are the average of the geometric and elastic stiffness matrices, respectively. \xi 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 non-sampling-based approach (Pc approach). Here, the non-intrusive Pc approach was used to obtain the probability space of critical buckling loads.

2.4 Non-intrusive polynomial chaos approach

PC expansion is used for defining random fields and random variables. The second-order variable X:\mathrm{\Theta }\to \mathbb{R} was considered, which can be described as possibly infinite variables of {\xi }_1(\theta ), {\xi }_2(\theta ), {\xi }_3(\theta ),… The independent variables were gathered into \mathrm{\Omega }(\theta )=({\xi }_1(\theta ),{\xi }_2(\theta ),{\xi }_3(\theta ),\dots ) and the variables, including stochastic eigenvalues and eigenvectors, can be obtained in the form of X=X(\mathrm{\Omega }), {N_{\text{c}}}_{\text{r}}=N(\mathrm{\Omega }). Equation (29) was utilized for defining stochastic eigenvectors related to Eq. (28), where the set of orthonormal polynomials is defined as \left\{{\psi }_0,{\psi }_1,{\psi }_2,\dots \right\}, which is a principle for L_2(\mathrm{\Theta },P). Also, deterministic coefficients x_0,x_1,x_2,\dots introduce X vectors. It can be obtained as Eqs. (29) and (30), if {\psi }_0(\xi )=1 based on orthonormal characteristic {⟨{\psi }_m,{\psi }_n⟩}_{\mathit{\text{θ}}}={\delta }_{mn}. Consequently, the variance of X is defined as Eq. (31). A similar definition is used for the definition of stochastic eigenvalues.

(29)X(\xi )=\sum _{n=0}^{+\mathrm{\infty }}{x}_n{\psi }_n(\xi ),

(30){\mu }_X=E\left\{X\right\}=\sum _{n=0}^{+\mathrm{\infty }}{x}_n{⟨{\psi }_0,{\psi }_n⟩}_{\text{θ}}=x_0,

(31){\sigma }_X^2=E\left\{X^2\right\}-E{\left\{X\right\}}^2=E\left\{{\left(\sum _{n=1}^{+\mathrm{\infty }}{x}_n{\psi }_n(\xi )\right)}^2\right\}=\sum _{n=1}^{+\mathrm{\infty }}{x}_n^2.

Since all numerical calculations require parameterization of the random variable X with a limited number of {\nu }_{\mathrm{r}\nu } of independent random variables, X can be defined as Eq. (32) if \xi =({\xi }_1,{\xi }_2,\dots ,{\xi }_{{\nu }_{\mathrm{r}\nu }}). {\nu }_{\mathrm{r}\nu } and p_{\mathrm{P}\mathrm{C}} are the result of random modeling and a compromise between accuracy and efficiency, respectively.

(32)X^{{\nu }_{\text{PC}}}(\xi )=\sum _{n=0}^{{\nu }_{\text{PC}}}{x}_n{\psi }_n(\xi ),

where p_{\text{PC}},{\nu }_{\mathrm{r}\nu }\to +\mathrm{\infty }. 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 X:D\times \mathrm{\Theta }\to \mathbb{R} and a set of random variables was placed in interval of X(u,0)\in L_2(\mathrm{\Theta },P). 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 non-intrusive method with linear regression was used so that random variables were considered as a finite number of {\nu }_{\text{gq}} and accordingly, the coefficients were obtained using x={\left({\mathrm{\Psi }}^{\text{T}}\mathrm{\Psi }\right)}^{-1}{\mathrm{\Psi }}^{\text{T}}X. Then, the final equation was extracted as Eq. (35), which can be applied to find the basis (x_i) of stochastic space of eigenvalue (critical load) and corresponding eigenvectors (buckling mode shapes) in v_{\text{pc}} design points instead of a full run based on the MC simulation.

(33)X(u,\xi )=\sum _{n=0}^{+\mathrm{\infty }}{x}_n(u){\psi }_n(\xi ),

(34)X^{{\nu }_{\text{pc}}}(u,\xi )=\sum _{n=0}^{{\nu }_{\text{pc}}}{x}_n(u){\psi }_n(\xi ),

(35)\underset{\mathrm{\Psi }}{\underbrace{\left[\begin{array}{ccc}{\psi }_0({\xi }_1)& \cdots & {\psi }_{{\nu }_{\text{pc}}}({\xi }_1)\\ ⋮& ⋮& ⋮\\ {\psi }_0({\xi }_{{\nu }_{\text{gq}}})& \cdots & {\psi }_{{\nu }_{\text{pc}}}({\xi }_{{\nu }_{\text{gq}}})\end{array}\right]}}\underset{x}{\underbrace{\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{{\nu }_{\text{pc}}}\end{array}\right]}}\approx \underset{X}{\underbrace{\left[\begin{array}{c}X(u,{\xi }_1)\\ ⋮\\ X(u,{\xi }_{{\nu }_{\text{gq}}})\end{array}\right]}}.

2.5 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.

2.6 Deterministic properties

The square composite with simple supported stacking sequence [\theta ,-\theta ] and lamina properties, according to Tab.1 was analyzed under uniaxial loading. The normalized buckling critical load was defined according to {\overline{\lambda }}_{\text{cr}}={\lambda }_{\text{cr}}{a}^2/E_2{h}^3. The comparison of our results with previous studies showed a high level of agreement, verifying our equations and codes for plates with deterministic properties (Tab.2).

2.7 Stochastic properties

The square composite plate with simply supported boundary condition, stacking sequence [0,90], and lamina properties (Tab.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.

3 Results and discussion

3.1 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 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. Fig.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.

3.2 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. Fig.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.

3.3 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 Tab.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 Tab.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.

4 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.

References

[1] Al-Jumaili S K, Eaton M J, Holford K M, Pearson M R, Crivelli D, Pullin R. Characterisation of fatigue damage in composites using an acoustic emission parameter correction technique. Composites. Part B, Engineering, 2018, 151: 237–244

[2] Maleki S, Rafiee R, Hasannia A, Habibagahi M R. Investigating the influence of delamination on the stiffness of composite pipes under compressive transverse loading using cohesive zone method. Frontiers of Structural and Civil Engineering, 2019, 13(6): 1316–1323

[3] Rajak D K, Pagar D D, Menezes P L, Linul E. Fiber-reinforced polymer composites: Manufacturing, properties, and applications. Polymers, 2019, 11(10): 1667

[4] Sakata S, Okuda K, Ikeda K. Stochastic analysis of laminated composite plate considering stochastic homogenization problem. Frontiers of Structural and Civil Engineering, 2015, 9(2): 141–153

[5] Nikbakht S, Kamarian S, Shakeri M. A review on optimization of composite structures Part II: Functionally graded materials. Composite Structures, 2019, 214: 83–102

[6] SchenkC ASchuëllerG I. Uncertainty assessment of large finite element systems. Lecture Notes in Applied and Computational Mechanics. Vol. 24 Series. New York: Springer Berlin Heidelberg, 2005

[7] FishmanG. Monte Carlo: Concepts, Algorithms, and Applications. Springer Series in Operations Research and Financial Engineering. New York: Springer Science & Business Media, 2013

[8] NayfehA H. Perturbation Methods. New York: John Wiley & Sons, 2008

[9] Chow P L. Perturbation methods in stochastic wave propagation. SIAM Review, 1975, 17(1): 57–81

[10] HosderSWaltersRPerezR. A non-intrusive polynomial chaos method for uncertainty propagation in CFD simulations. In: 44th AIAA aerospace sciences meeting and exhibit. Nevada: American Institute of Aeronautics and Astronautics, 2006

[11] GelfandA EDeyD KChangH. Model Determination Using Predictive Distributions with Implementation via Sampling-Based Methods. Technical Report 462. Department of Statistics, Stanford University. 1992

[12] Gelfand A E. Model Determination Using Sampling-Based Methods. Markov Chain Monte Carlo in Practice, Chapter 9. London: Chapman & Hall, 1996, 145–161

[13] Sepahvand K, Marburg S, Hardtke H J. Uncertainty quantification in stochastic systems using polynomial chaos expansion. International Journal of Applied Mechanics, 2010, 2(2): 305–353

[14] Bisagni C. Numerical analysis and experimental correlation of composite shell buckling and post-buckling. Composites. Part B, Engineering, 2000, 31(8): 655–667

[15] Telford R, Peeters D, Rouhi M, Weaver P M. Experimental and numerical study of bending-induced buckling of stiffened composite plate assemblies. Composites. Part B, Engineering, 2022, 233: 109642

[16] Rozylo P, Teter A, Debski H, Wysmulski P, Falkowicz K. Experimental and numerical study of the buckling of composite profiles with open cross section under axial compression. Applied Composite Materials, 2017, 24(5): 1251–1264

[17] Zhang Y, Tao W, Chen Y, Lei Z, Bai R, Lei Z. Experiment and numerical simulation for the compressive buckling behavior of double-sided laser-welded Al–Li alloy aircraft fuselage panel. Materials (Basel), 2020, 13(16): 3599

[18] Ly H B, Desceliers C, Minh Le L, Le T T, Thai Pham B, Nguyen-Ngoc L, Doan V T, Le M. Quantification of uncertainties on the critical buckling load of columns under axial compression with uncertain random materials. Materials (Basel), 2019, 12(11): 1828

[19] Sharma N, Nishad M, Maiti D K, Sunny M R, Singh B N. Uncertainty quantification in buckling strength of variable stiffness laminated composite plate under thermal loading. Composite Structures, 2021, 275: 114486

[20] Kharghani N, Soares C. Effect of uncertainty in the geometry and material properties on the post-buckling behavior of a composite laminate. Maritime Technology and Engineering, 2016, 3: 497–503

[21] Nguyen H X, Duy Hien T, Lee J, Nguyen-Xuan H. Stochastic buckling behaviour of laminated composite structures with uncertain material properties. Aerospace Science and Technology, 2017, 66: 274–283

[22] Hu L, Feng P, Meng Y, Yang J. Buckling behavior analysis of prestressed CFRP-reinforced steel columns via FEM and ANN. Engineering Structures, 2021, 245: 112853

[23] Dey S, Mukhopadhyay T, Spickenheuer A, Gohs U, Adhikari S. Uncertainty quantification in natural frequency of composite plates—An Artificial neural network based approach. Advanced Composites Letters, 2016, 25(2): 43–48

[24] Sasikumar P, Venketeswaran A, Suresh R, Gupta S. A data driven polynomial chaos based approach for stochastic analysis of CFRP laminated composite plates. Composite Structures, 2015, 125: 212–227

[25] Dey S, Mukhopadhyay T, Sahu S, Li G, Rabitz H, Adhikari S. Thermal uncertainty quantification in frequency responses of laminated composite plates. Composites. Part B, Engineering, 2015, 80: 186–197

[26] Chandrashekhar M, Ganguli R. Damage assessment of composite plate structures with material and measurement uncertainty. Mechanical Systems and Signal Processing, 2016, 75: 75–93

[27] Swain P R, Dash P, Singh B N. Stochastic nonlinear bending analysis of piezoelectric laminated composite plates with uncertainty in material properties. Mechanics Based Design of Structures and Machines, 2021, 49(2): 194–216

[28] Singh B N, Iyengar N, Yadav D. Effects of random material properties on buckling of composite plates. Journal of Engineering Mechanics, 2001, 127(9): 873–879

[29] Kalfountzos C D, Bikakis G S, Theotokoglou E E. Deterministic and probabilistic buckling response of fiber–metal laminate panels under uniaxial compression. Aircraft Engineering and Aerospace Technology, 2022, 94(5): 745–759

[30] Zhuang X, Guo H, Alajlan N, Zhu H, Rabczuk T. Deep autoencoder based energy method for the bending, vibration, and buckling analysis of Kirchhoff plates with transfer learning. European Journal of Mechanics. A, Solids, 2021, 87: 104225

[31] Guo H, Zhuang X, Rabczuk T. A deep collocation method for the bending analysis of Kirchhoff plate. Computers, Materials & Continua, 2019, 59(2): 433–456

[32] Guo H, Zheng H, Zhuang X. Numerical manifold method for vibration analysis of Kirchhoff’s plates of arbitrary geometry. Applied Mathematical Modelling, 2019, 66: 695–727

[33] Guo H, Zheng H. The linear analysis of thin shell problems using the numerical manifold method. Thin-walled Structures, 2018, 124: 366–383

[34] Vu-Bac N, Duong T X, Lahmer T, Zhuang X, Sauer R A, Park H S, Rabczuk T. A NURBS-based inverse analysis for reconstruction of nonlinear deformations of thin shell structures. Computer Methods in Applied Mechanics and Engineering, 2018, 331: 427–455

[35] Vu-Bac N, Duong T, Lahmer T, Areias P, Sauer R A, Park H S, Rabczuk T. A NURBS-based inverse analysis of thermal expansion induced morphing of thin shells. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 480–510

[36] Vu-Bac N, Rabczuk T, Park H, Fu X, Zhuang X. A NURBS-based inverse analysis of swelling induced morphing of thin stimuli-responsive polymer gels. Computer Methods in Applied Mechanics and Engineering, 2022, 397: 115049

[37] Vu-Bac N, Lahmer T, Keitel H, Zhao J, Zhuang X, Rabczuk T. Stochastic predictions of bulk properties of amorphous polyethylene based on molecular dynamics simulations. Mechanics of Materials, 2014, 68: 70–84

[38] Vu-Bac N, Lahmer T, Zhang Y, Zhuang X, Rabczuk T. Stochastic predictions of interfacial characteristic of polymeric nanocomposites (PNCs). Composites. Part B, Engineering, 2014, 59: 80–95

[39] Vu-Bac N, Silani M, Lahmer T, Zhuang X, Rabczuk T. A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites. Computational Materials Science, 2015, 96: 520–535

[40] AdhikariS. Free vibration analysis of angle-ply composite plates with uncertain properties. In: 17th AIAA Non-Deterministic Approaches Conference. Florida: American Institute of Aeronautics and Astronautics, 2015

[41] Vu-Bac N, Rafiee R, Zhuang X, Lahmer T, Rabczuk T. Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters. Composites. Part B, Engineering, 2015, 68: 446–464

[42] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

[43] Vu-Bac N, Zhuang X, Rabczuk T. Uncertainty quantification for mechanical properties of polyethylene based on fully atomistic model. Materials (Basel), 2019, 12(21): 3613

[44] Liu B, Vu-Bac N, Zhuang X, Rabczuk T. Stochastic multiscale modeling of heat conductivity of Polymeric clay nanocomposites. Mechanics of Materials, 2020, 142: 103280

[45] Liu B, Vu-Bac N, Zhuang X, Fu X, Rabczuk T. Stochastic full-range multiscale modeling of thermal conductivity of Polymeric carbon nanotubes composites: A machine learning approach. Composite Structures, 2022, 289: 115393

[46] Liu B, Vu-Bac N, Rabczuk T. A stochastic multiscale method for the prediction of the thermal conductivity of Polymer nanocomposites through hybrid machine learning algorithms. Composite Structures, 2021, 273: 114269

[47] Fakoor M, Parviz H. Uncertainty propagation in dynamics of composite plates: A semi-analytical non-sampling-based approach. Frontiers of Structural and Civil Engineering, 2020, 14(6): 1359–1371

[48] Fakoor M, Parviz H, Abbasi A. Uncertainty propagation analysis in free vibration of uncertain composite plate using stochastic finite element method. Amirkabir Journal of Mechanical Engineering, 2019, 52(12): 3503–3520

[49] Sriramula S, Chryssanthopoulos M K. An experimental characterisation of spatial variability in GFRP composite panels. Structural Safety, 2013, 42: 1–11

[50] Ghanem R G, Spanos P D. Stochastic finite element method: Response statistics. In: Stochastic Finite Elements: A Spectral Approach. New York: Springer, 1991, 101–119

[51] Fish J, Wu W. A nonintrusive stochastic multiscale solver. International Journal for Numerical Methods in Engineering, 2011, 88(9): 862–879

[52] ReddyJ N. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. 2nd ed. New York: CRC Press, 2003

[53] Kim S E, Thai H T, Lee J. A two variable refined plate theory for laminated composite plates. Composite Structures, 2009, 89(2): 197–205

[54] Tran L V, Thai C H, Le H T, Gan B S, Lee J, Nguyen-Xuan H. Isogeometric analysis of laminated composite plates based on a four-variable refined plate theory. Engineering Analysis with Boundary Elements, 2014, 47: 68–81

RIGHTS & PERMISSIONS

The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn