An Efficient Anti-Optimization Approach for Uncertainty Analysis in Composite Laminates

Abstract

This work presents an efficient approach to quantify uncertainties in composite laminates using the interval analysis, anti-optimization technique, and the α-cut procedure. The solutions are compared with the traditional and robust Monte Carlo method in 3 cases scenarios: natural frequencies, buckling, and strength safe factor. For natural frequencies and buckling loads, the presented Interval based methodology showed 2.5% to 4.5% larger error values when compared to the Monte Carlo method using the same number of function calls. This implies a larger uncertain area, and hence, a better solution. For the strength test using Tsai-Wu failure theory, the error values are even greater: 22% to 46%. A violation of the failure limit was detected by the proposed Interval based approach, but not detected by Monte Carlo method. The solutions show that the presented methodology yields a safer and more precise analysis when compared to the traditional Monte Carlo approach.

Keywords:
Anti-optimization; laminated composites; interval-based uncertainty analysis; convex hull

1. INTRODUCTION

Composite materials are used in all kinds of high-performance structural applications, such as in aircraft, aerospace and the military industry. Over the past few decades, numerous advances have been achieved in the research field of composite materials. These advances have provided designers with powerful tools to more accurately predict the behavior of materials, as well as techniques for producing highly complex parts. Despite the significant advances, there is still a relevant discrepancy between apparently identical laminate specimens when performing experimental measurements. The uncertainties are often pointed as the main cause of the aforementioned divergences in experimental data.

The present uncertainty in composite materials can be a big issue in the design and manufacture phases of the layers that compose the final structure. These uncertainties can be inserted into the system in a number of ways, such as through temperature sensitivity, hygroscopic phenomena, and the own manufacturing process. Even the standardized tests used to characterize the elastic properties of the composite material have some degree of uncertainty, as explained by¹. The lack of accurate material properties (required in order to predict failure mechanisms, dynamic behavior and deformation) may lead to excessively conservative designs^2-4. For instance, in⁵, different fiber orientations and orthotropic proprieties of composite laminates were tested. The authors showed experimentally and analytically that those proprieties affect the interlaminar shear strength and flexural properties, and hence, the overall mechanical performance, justifying the choice of these values as uncertain. An extensive bibliography exists in quantifying the uncertainties like these, as is discussed in²_.

According to⁶, current approaches to uncertainty problems are divided into two major groups: possibilistic modeling of uncertainty and probabilistic modeling of uncertainty. The probabilistic approach requires knowing the probabilistic distribution of the uncertainties and the correlation between variables, which is not generally known in practice. On the other hand, the possibilistic method only requires the range of uncertainty possibilities, which is promptly available. Rao and Mash hour⁷ used the possibilistic approach, interval arithmetic operations and truncation-based interval analysis to perform failure analysis on composite laminates in a fairly simple way. The interval analysis used in this work also fits the possibilistic approach but uses a more advanced methodology than the arithmetic operations aforementioned. In this paper, the developed algorithm is based on the α-cut concept and convex hull from the quick hull algorithm⁸, as explained in⁹.

For the Interval Approach used in the algorithm framework of the present work, the measure of confidence on the uncertainty level of variables can be associated to an $\alpha$-cut level, Figure 1(b), a single value that varies from 0 to 1 (where 0 is the maximum uncertainty in the interval limit and 1 being complete confidence and thus, no uncertainty at all).

Figure 1
(a) Interval uncertainty propagation. (b) Linear approximation of Gaussian distribution for x by a triangular fuzzy-set, with uncertain intervals and α-cuts.

Many methods for evaluating uncertainty in composites have been developed over the past few decades. Several versions of SFEM (Stochastic Finite Element Method) have been proposed to approach the uncertainty problem, such as^10-14. More recently, anti-optimization has emerged as a viable alternative to solve uncertainties problems in laminates^15-17. The method has shown to be well suited in solving such problems, and even in design optimization, where the anti-optimization method is nested within another optimization, like in^18,19. In this approach, while the anti-optimization searches for the worst-case scenario using the uncertainty, the outer optimization handles the design variables, thus leading to a heavy computational load. Chen et al.²⁰ show that it is possible to approach the uncertainty in laminates using Taylor expansion. Another method that has been used in recent years in this area and all kinds of uncertainty problems is the polynomial chaos expansion^21,22 and other meta-models like the Kriging²³.

The aleatory uncertainty in composite materials is also treated through three distinct approaches on account of the model scale⁴: micro-scale, meso-scale, and macro-scale. Some authors, such as^14,24-26, use multiple instances of the aforementioned approaches, called multi-scale analysis. Using a multi-scale analysis (meso and macro-scale)^14,26, the authors applied the Representative Volume Element (RVE) method coupled with finite element analysis to produce better results when compared with Monte Carlo alone. Since the present paper uses ply level properties, it fits into the meso-scale category⁴.

Several simplified methods can be used to evaluate the propagation of the uncertainties intervals, such as gradient-based methods, arithmetic intervals⁷, but none of them is sufficiently accurate for all levels of uncertainty. For this reason, some more elaborated approaches must be used, such as the anti-optimization. For the anti-optimization approach, an optimization algorithm is used to find the limits of the output variables that represent the system's behavior based on combinations of the uncertain variables. The Monte Carlo approach is also useful for this concern, being a robust solution to the uncertainty propagation problem. As explained in², the Monte Carlo simulation technique in conjunction with the finite element (FE) method is widely used for quantifying uncertainties of laminated composite structures, as it is done in^3,14. The comparison between the Monte Carlo and the anti-optimization approach is explored later in this work, where it is shown that Monte Carlo is not the best choice regarding performance and efficiency.

This paper will be presented in the following order: In section 2, the uncertainty propagation and the anti-optimization method will be discussed; in section 3, the formulation for composite materials will be presented; in section 4, the examples and the numerical solutions will be shown (natural frequencies, buckling load, and strength safe factor), as well as some discussion about the results; finally, section 5 concludes the work, summarizing the obtained results.

2. UNCERTAINTY PROPAGATION AND ANTI-OPTIMIZATION

Assume a system with the input vector ${}_{}{}^{m}X={({}_{}{}^{m}x{}_1,{}_{}{}^{m}x{}_2,\dots ,{}_{}{}^{m}x{}_{ni})}^T$, where $ni$ is the number of input variables, $m$ means "measured", and an output vector as ${}_{}{}^{m}Z={({}_{}{}^{m}z{}_1,{}_{}{}^{m}z{}_2,\dots ,{}_{}{}^{m}z{}_{no})}^T$, where no is the number of output variables. Let's assume a number of samples of each one of the input vectors (called realizations)${}_1{}^{m}X$, ${}_2{}^{m}X$,…,${}_{ns}{}^{m}X$ and the corresponding output vectors ${}_1{}^{m}Z$, ${}_2{}^{m}Z$,…,${}_{ns}{}^{m}Z$ , where $ns$ means the number of samples. Furthermore, let's assume a numerical model of the system with the same measured inputs ${}_{}{}^{m}X={({}_{}{}^{m}x{}_1,{}_{}{}^{m}x{}_2,\dots ,{}_{}{}^{m}x{}_{ni})}^T$, adjustable model's parameters vector $\theta ={({\theta }_1,{\theta }_2,\dots ,{\theta }_{np})}^T$, with $np$ representing the number of parameters and predicted vector outputs ${}^{p}Z={({}^{p}z{}_1,{}^{p}z{}_2,\dots ,{}^{p}z{}_{no})}^T=f\left({}^{m}X,\theta \right)$. Regarding the output vector, it is easy to obtain their upper and lower bounds based on the values from all realizations. These output intervals vectors may be put together in a two-column vector representation as:

$\left[{}_{}{}^m\underset{\_}{z},{}_{}{}^m\overline{z}\right]=\left[\underset{i=\mathrm{1,}\dots ,ns}{\min }\left({}_i{}^{m}z\right)\underset{i=\mathrm{1,}\dots ,ns}{\max }\left({}_i{}^{m}z\right)\right]$ (1)

Equation (1) represents a measure of the model's output uncertainty or dispersion. The same can be stated for the predicted output vector (from the numerical model) and the measured input vector:

$$\left[{}^m\underset{¯}{X},{}^m\overline{X}\left]=\right[\underset{i=\mathrm{1,}\dots ,ns}{\min }\left({}_i{}^{m}X\right)\underset{i=\mathrm{1,}\dots ,ns}{\max }\left({}_i{}^{m}X\right)\right]$$ (2)

$$\left[\underset{\_}{\theta },\overline{\theta }\right]=[\underset{i=\mathrm{1,}\dots ,ns}{\min }\left({}_i\theta \right)\underset{i=\mathrm{1,}\dots ,ns}{\max }\left({}_i\theta \right)$$ (3)

$$\begin{array}{l}\left[{}^p\underset{\_}{z},{}^p\overline{z}\right]=f\left(\left[{}^m\underset{¯}{X},{}^m\overline{X}\right],\left[\underset{\_}{\theta },\overline{\theta }\right]\right)=\\ \left[\underset{i=\mathrm{1,}\dots ,ns}{\min }\left({}_i{}^{p}Z\left(\left[\underset{\_}{\theta },\overline{\theta }\right]\right)\right)\underset{i=\mathrm{1,}\dots ,ns}{\max }\left({}_i{}^{p}Z(\left[\underset{\_}{\theta },\overline{\theta }\right]\right)\right]\end{array}$$ (4)

A very common question is: What are the output intervals $\left[{}_{}{}^p\underset{\_}{z},{}_{}{}^p\overline{z}\right]$ for a given set of measured input intervals $\left[{}_{}{}^m\underset{¯}{X},{}_{}{}^m\overline{X}\right]$ and parameters intervals $\left[\underset{\_}{\theta },\overline{\theta }\right]$? This may be thought of as finding the variability of the response of the system (upper and lower bounds of the output behavior). This also can be understood as an uncertainty propagation in terms of intervals since the uncertainty of the inputs and parameters will be propagated to the outputs. This is as a double optimization (minimization and maximization, also known as anti-optimization) in each of the hypervolume quadrants (2^no) to find the best input ${}_{}{}^{m}X{}^{*}$. and parameter values${\theta }^{*}$ that maximizes (or minimizes) the output norm ${}_{}{}^{p}Z{}_2$ and that colies with the bounded intervals. The final output intervals will be the minimum and maximum of each variable for all hyper-quadrants. Mathematically, this can be stated as:

$\text{Find}{X}^{*}\text{and}{\theta }^{*},\text{at each hyperquadrant},\text{that minimizes and maximizes}{Z}^{*}{}_2$

$$\begin{array}{l}\text{Subject to :}{X}^{*}\in \left[{}^m\underset{¯}{X},{}^m\overline{X}\right],{\theta }^{*}\in \left[\underset{\_}{\theta },\overline{\theta }\right]\text{and}\left[\underset{\_}{z},\overline{z}\right]=\\ \left[\underset{overallquadrants}{min\left(z^{*}\right)}\underset{overallquadrants}{max\left(z^{*}\right)}\right]\end{array}$$ (5)

In this case, where a double optimization happens, the number of design variables is (ni+np). A block diagram can be sketched in this sense, as represented in Figure 1(a).

The Interval upper and lower values form a cloud in the hyperspace and, once obtained, a convenient way to visualize is to construct a convex hull, which is the smallest convex set that contains and encompasses the set of points. This convex hull can also be constructed for the random points obtained by the simulations if the Monte Carlo method is used.

3. COMPOSITE MATERIAL

Laminated composite structures can be challenging to manufacture accurately according to their exact design specifications because of the inherent complexity present in such materials. As a result, undesirable and unavoidable uncertainties arise in composite manufacturing, such as intra-laminate voids, incomplete curing of the resin, excess resin between plies, porosity, excess matrix voids, variations in ply thickness and fiber parameters, all those making the overall behavior prediction of such materials difficult^2,3. The lack of precise tools to evaluate the propagation of these uncertainties may lead to overly conservative or unsafe designs⁴, with a more expensive and less effective final product.

In this study, mass density, ply angle orientation and Elastic modulus values will be considered with uncertain values since those are the most important source of the variability described in the literature².

3.1. Equations for dynamic vibrations of laminated plates

The differential equations for a laminated plate can be derived from analyzing Figure 2, where a force diagram is depicted (comma represents the corresponding partial derivatives).

Figure 2
Force diagram at the middle surface of a laminated plate.

By the resulting forces in x and y directions, we find

$N_{x,x}+N_{xy,y}=\mathrm{0,}$ (6)

$N_{y,y}+N_{xy,x}=0$ (7)

Considering the inertia force (density $\rho$ and height $h$) and the forces in the z-direction, one finds the equilibrium equation

$Q_{x,x}+Q_{y,y}+N_x{w}_{,xx}+2N_{xy}{w}_{,xy}+N_y{w}_{,yy}-\rho h{w}_{,tt}=0$ (8)

If one neglects the third-order terms, the moment equations result in Equation (9) and Equation (10).

$M_{x,x}+M_{xy,y}=Q_x,$ (9)

$M_{y,y}+M_{xy,x}=Q_y,$ (10)

and substituting Equations (9) and (10) into Equation (8), we find

$\begin{array}{l}{M}_{x,xx}+2M_{xy,yy}+M_{y,yy}+N_x{w}_{,xx}+\\ 2N_{xy}{w}_{,xy}+N_y{w}_{,yy}-\rho h{w}_{,tt}=0\end{array}$ (11)

Equations (11), (6) and (7) forms a set of three differential equations of motion. For a laminated plate, the constitutive relations for stress/strains at the local coordinate system for laminate k at height z is (assuming First-Order Shear Deformation Theory, FSDT):

${\left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}}_k={\left[\begin{array}{ccc}{\overline{Q}}_{11}& {\overline{Q}}_{12}& {\overline{Q}}_{16}\\ {\overline{Q}}_{12}& {\overline{Q}}_{22}& {\overline{Q}}_{26}\\ {\overline{Q}}_{16}& {\overline{Q}}_{26}& {\overline{Q}}_{66}\end{array}\right]}_k\left(\left\{\begin{array}{c}{\epsilon }_x^0\\ {\epsilon }_y^0\\ {\gamma }_{xy}^0\end{array}\right\}+z\left\{\begin{array}{c}{k}_x\\ k_y\\ k_{xy}\end{array}\right\}\right)$,(12)

and the internal forces evaluated as:

$\left\{\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\end{array}\right\}={{\displaystyle \int }}_{-t/2}^{+t/2}{\left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}}_{k}dz$and$\left\{\begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\right\}={{\displaystyle \int }}_{-t/2}^{+t/2}{\left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}}_{k}zdz$,(13)

and applying Equation (12), one recovers:

$\left\{\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\end{array}\right\}=A_{ij}\left\{\begin{array}{c}{\epsilon }_x^0\\ {\epsilon }_y^0\\ {\gamma }_{xy}^0\end{array}\right\}+B_{ij}\left\{\begin{array}{c}{k}_x\\ k_y\\ k_{xy}\end{array}\right\}$ and$\left\{\begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\right\}=B_{ij}\left\{\begin{array}{c}{\epsilon }_x^0\\ {\epsilon }_y^0\\ {\gamma }_{xy}^0\end{array}\right\}+D_{ij}\left\{\begin{array}{c}{k}_x\\ k_y\\ k_{xy}\end{array}\right\}$ (14)

The matrices $A$, B and $D$ are evaluated as usual, along each ply, as indicated by Equations (15), (16) and (17).

$A_{ij}={{\displaystyle \sum }}_{k=1}^n\kappa {\left({\overline{Q}}_{ij}\right)}_k\left(z_k-z_{k-1}\right),$ (15)

$B_{ij}=\frac{1}{2}{{\displaystyle \sum }}_{k=1}^n{\left({\overline{Q}}_{ij}\right)}_k\left(z_k{}^2-z_{k-1}{}^2\right),$ (16)

$D_{ij}=\frac{1}{3}{{\displaystyle \sum }}_{k=1}^n{\left({\overline{Q}}_{ij}\right)}_k\left(z_k{}^3-z_{k-1}{}^3\right),$ (17)

where$\kappa =5/6$. Based on the Classical Laminated Theory, taking the strain displacement relationship, one finds Equation (18):

$\left\{\begin{array}{c}{\epsilon }_x^0\\ {\epsilon }_y^0\\ {\gamma }_{xy}^0\end{array}\right\}=+B_{ij}\left\{\begin{array}{c}{u}_{,x}\\ v_{,y}\\ u_{,y}+v_{,x}\end{array}\right\}$and$\left\{\begin{array}{c}{k}_x\\ k_y\\ k_{xy}\end{array}\right\}=\left\{\begin{array}{c}-w_{,xx}\\ -w_{,yy}\\ -2w_{,xy}\end{array}\right\}$,(18)

and finally, the set of the three differential equations, taking into consideration Equation (17), results in Equation (19):

$$\begin{array}{l}{A}_{11}{u}_{,xx}+2A_{16}{u}_{,xy}+A_{66}{u}_{,yy}+A_{16}{v}_{,xx}+\left(A_{11}+A_{66}\right)v_{,xy}+\\ A_{26}{v}_{,yy}-B_{11}{w}_{,xxx}-3B_{16}{w}_{,xxy}-\left(B_{12}+2B_{66}\right)w_{,xyy}-B_{26}{w}_{,yyy}=0\\ A_{16}{u}_{,xx}+\left(A_{12}+A_{66}\right)u_{,xy}+A_{26}{u}_{,yy}+A_{66}{v}_{,xx}+\\ 2A_{26}{v}_{,xy}+A_{22}{v}_{,yy}-B_{22}{w}_{,yyy}=0\\ D_{11}{w}_{,xxxx}+4D_{16}{w}_{,xxxy}+2\left(D_{12}+2D_{66}\right)w_{,xxyy}+4D_{26}{w}_{,xyyy}+\\ D_{22}{w}_{,yyyy}-B_{11}{u}_{,xxx}-3B_{16}{u}_{,xxy}-\left(B_{12}+2B_{66}\right)u_{,xyy}\\ -B_{26}{u}_{,yyy}-B_{16}{v}_{,xxx}-\left(B_{12}+2B_{66}\right)v_{,xxy}-3B_{26}{v}_{,xyy}-\\ B_{22}{v}_{,yyy}+N_x{w}_{,xy}+N_y{w}_{,yy}+\rho h{w}_{,tt}=0\end{array}$$ (19)

In the case of symmetrically angle-ply oriented composite plates, $A_{16}=A_{26}=D_{16}=D_{26}=B_{ij}=0$ and the previous equations simplify. As shown in²⁷, to solve such equation, one can use the separability of space and time with Fourier series in the form

$w\left(x,y,t\right)=[{\alpha }_{mn}\cos \left(\omega t\right)+{\beta }_{mn}\sin \left(\omega t\right)]X_m\left(x\right)Y_n\left(y\right)$,(20)

and considering the corresponding boundary conditions, results in a system of equations that will give the natural frequencies and mode shapes only if a nontrivial solution exists. In this case, for a simply supported plate with edges a and b, the following equation for natural frequencies results:

${\omega }_{m,n}^2=\frac{{\pi }^4}{\rho h}[D_{11}{(\frac{m}{a})}^4+2\left(D_{12}+2D_{66}\right){(\frac{m}{a})}^2{(\frac{n}{b})}^2+D_{22}\left(\frac{n}{b}{)}^4\right]$ ,(21)

where $m$ and $n$ are positive integers that should be tried to get each of the mode frequencies.

3.2. Equations for buckling of laminated plates

Assuming that the only applied load are those from in-plane (no shear forces), we can split the effects in the total potential energy ($\Pi$) into two parts, as shown in Equation (22).The first part is the strain energy due to bending ($U_b$), and the second one, the work done by external forces ($U_p)$, represented by Equation (22)

$\text{Π}=U_b-U_p$ (22)

One can evaluate those parts ($U_b$ and $U_p$) for an assumed field solution that accounts for boundary conditions of the form of sine series. The critical buckling load will be evaluated for a clamped plate under uniaxial load in the second example presented in this work. In this case, it is assumed a solution in the form

$w\left(x,y\right)=w_{m,n}\sin \left(m\pi x/a\right)\sin \left(m\pi y/b\right),$ (23)

with $w_{m,n}$ being the displacement coefficients and $m$, $n$, positive integers. Substituting this and assuming displacement field into the total potential energy, forcing a stationary condition ($\partial \Pi /\partial w_{mn}$=0) and solving for $\lambda$=$N_{crit}/N_x$ (critical to the reference applied load ratio), the buckling loads are computed. In fact, the first buckling load is of interest, and one should search for the combinations of $m$ and $n$that gives the lower $\lambda$ (first buckling load).In case of simply supported composite plates, the solution is in the form

${\lambda }_{m,n}=-(D_{11}{m}^4+2\left(D_{12}+D_{66}\right)m^2{n}^2(\frac{a}{b}){}^2+2D_{22}{n}^4{(\frac{a}{b})}^4/\left(m^2{N}_x+n^2{(\frac{a}{b})}^2{N}_y\right)$ (24)

In case of a clamped plate with uniaxial load, N_x, the resulting expression is given by Equations (25) and (26)

${\lambda }_{m,n}=-(G_1{G}_5+G_2{G}_4)-\sqrt{\begin{array}{l}(G_1{G}_5+G_2{G}_4)²-4\\ \left(G_2{G}_5-G_3^2\right)G_1{G}_4\end{array}}/\left[2\left(G_2{G}_5-G_3^2\right)\right]$ $form=1|\text{n}\in \text{N}$ (25)

${\lambda }_{m,n}=-(H_1{H}_5+H_2{H}_4-\frac{\sqrt{\begin{array}{l}(H_1{H}_5+H_2{H}_4)²-4\\ \left(H_2{H}_5-H_3^2\right)H_1{H}_4\end{array}}}{\left[2\left(H_2{H}_5-H_3^2\right)\right]}$ $form\ne 1|\left\{\text{m},\text{n}\right\}\in \text{N}$ (26)

where

$G_1=\frac{{\pi }^4\left(3D_{11}{b}^4+2D_{12}{a}^2{b}^2+3D_{22}{a}^4+4D_{66}{a}^2{b}^2\right)}{4a^3{b}^3}$ (27)

$G_2=-\frac{{\pi }^{2}b\left(3n^4-12n^2\right)}{16a\left(n^4-4n^2\right)}$ (28)

$G_3=\frac{8b\left(-n^3+{\left(-1\right)}^{n+2}{n}^3\right)}{3a\left(n^4-4n^2\right)}$ (29)

$G_4=\frac{k{\pi }^2\left(A_{44}{a}^2+A_{55}{n}^2{b}^2\right)}{4ab}$ (30)

$G_5=-\frac{{\pi }^{2}b\left(n^6-4n^4\right)}{4a\left(n^4-4n^2\right)}$ (31)

$H_1=\frac{{\pi }^4\left(\begin{array}{c}{D}_{11}{b}^4\left(18m^2+3+3m^4\right)+D_{12}{a}^2{b}^2\left(8+8m^2\right)+\\ 16D_{22}{a}^4+D_{66}{a}^2{b}^2\left(16+16m^2\right)\end{array}\right)}{32a^3{b}^3}$ (32)

$H_2=-\frac{{\pi }^{2}b\left(\begin{array}{l}3m^6+3+3n^4-3m^4-3m^2+\\ 3m^2{n}^4-12m^2{n}^2-6n^2-4m^4{n}^2\end{array}\right)}{32a\left(1-2n^2+m^4-2m^2{n}^2-2m^2+n^4\right)}$ (33)

$H_3=-\frac{8b\left(-m{n}^3+{\left(-1\right)}^{m+n+1}m{n}^3\right)}{3a\left(1-2n^2+m^4-2m^2{n}^2-2m^2+n^4\right)}$ (34)

$H_4=\frac{k{\pi }^2\left(A_{44}{a}^2+A_{55}{n}^2{b}^2\right)}{4ab}$ (35)

$H_5=-\frac{{\pi }^{2}b\left(n^2+n^6-2m^2{n}^4+m^4{n}^2-2n^4-2m^2{n}^2\right)}{4a\left(1-2n^2+m^4-2m^2{n}^2-2m^2+n^4\right)}$ (36)

3.3. Equations for strength evaluation in composite materials

For strength evaluation in composite materials, this paper uses the Tsai-Wu failure criteria²⁸. According to²⁹, Tsai-Wu failure theory is based on the total strain energy failure theory of Beltrami applied to lamina in a plane stress state. The following equation is used to evaluate the failure criteria in the lamina. If the equation is violated, the lamina is considered to be failed.

$H_1{\sigma }_1+H_2{\sigma }_2+H_6{\tau }_{12}+H_{11}{\sigma }_1^2+H_{22}{\sigma }_2^2+H_{66}{\tau }_{12}^2+2H_{12}{\sigma }_1{\sigma }_2-1<0$ (37)

or, in terms of a safety factor, $\alpha$:

$\begin{array}{l}(H_1{\sigma }_1+H_2{\sigma }_2+H_6{\tau }_{12})\alpha +(H_{11}{\sigma }_1^2+H_{22}{\sigma }_2^2+\\ H_{66}{\tau }_{12}^2+2H_{12}{\sigma }_1{\sigma }_2){\alpha }^2-1=0\end{array}$,(38)

where $\lambda =\alpha -1$ lower than 0 represents failure.

The components $H_1,H_2,H_6,H_{11},H_{22}$ and $H_{66}$ can be obtained based on the strength limit parameters of the material as follows

$H_1=\frac{1}{{\left({\sigma }_1^T\right)}_{ult}}-\frac{1}{{\left({\sigma }_1^C\right)}_{ult}}$ (39)

$H_{11}=\frac{1}{{\left({\sigma }_1^T\right)}_{ult}{\left({\sigma }_1^C\right)}_{ult}}$ (40)

$H_2=\frac{1}{{\left({\sigma }_2^T\right)}_{ult}}-\frac{1}{{\left({\sigma }_2^C\right)}_{ult}}$ (41)

$H_{22}=\frac{1}{{\left({\sigma }_2^T\right)}_{ult}{\left({\sigma }_2^C\right)}_{ult}}$ (42)

$H_6=0$ (43)

$H_{66}=\frac{1}{{\left({\tau }_{12}\right)}_{ult}^2}$ (44)

$H_{12}=\frac{1}{2}\sqrt{\frac{1}{{\left({\sigma }_1^T\right)}_{ult}{\left({\sigma }_1^C\right)}_{ult}{\left({\sigma }_2^T\right)}_{ult}{\left({\sigma }_2^C\right)}_{ult}}}$ (45)

Note that $H_{12}$ is found using some empirical suggestions, where the one presented here is found by Mises-Hencky criterion²⁹.

4. NUMERICAL EXAMPLES

4.1. Example 1 - Natural frequency of a simply supported plate

The following example shows a composite with the total length and breadth of the laminated plate with orthotropic cross-ply [90°/0°/0°/90°] configuration, dimensions a=0.6 m and b=0.6 m and laminate thickness of t=1.25×10^-4 m. The material properties of the laminate are E₁₁=135 GPa, E₂₂=8.8 GPa, G₁₂=4.8 GPa, ν₁₂=0.33 and ρ=1380kg/m³ (carbon-epoxy properties²⁹). The boundary conditions are simply supported on each side (SSSS). In this example it is assumed an uncertainty interval of ±10%, for the worst-case scenario (corresponding to α-cut level 0) in E₁₁, ρ, G₁₂ and an uncertainty ±5^o for θ(very common range of uncertainty found in²). This uncertainty setup is assumed to be independent for each layer, so the number of uncertainties is 16. The proposed algorithm is tested using the interval approach proposed in this paper and compared with simple MC (Monte Carlo simulations). For fair comparisons, the number of MC simulations equals the number of function calls used in the Interval-based proposed algorithm (that runs first). The uncertainty in the dynamic behavior is sought in this example, so the first four natural frequencies are evaluated. The analytical solutions for the deterministic situation for this particular case are given by Equation (21). The following Figure 3 is obtained for the convex hulls for the first 4 natural frequencies (Hz).

Figure 3
Output uncertain intervals and convex hulls for the interval-based procedure and MC Simulation for mode frequencies, where $f_i$ are the mode frequencies values in Hz.

Table 1 shows the relative differences between the Interval-based Method (using anti-optimization) proposed here and the traditional Monte Carlo Method based on the Error (E) values defined as the ratio between twice the Interval Radius of the uncertain output variables (first 4 natural frequencies) for the zero α-cut level and the Interval Center. The Interval Radius for an uncertain variable $z$ is defined as $IR=\left(\overline{z}-\underset{\_}{z}\right)/2$, while the Interval Center is defined as $IC=\left(\overline{z}+\underset{\_}{z}\right)/2$, so $E=2IR/IC$.

Table 1
Error-values for interval output variables obtained by the Interval-based Method and the MC method.

Mode Frequency (Hz)	Interval-based Method			MC Method
	Lower bound $\underset{\_}{f}$	Upper bound $\overline{f}$	E (%)	Lower bound $\underset{\_}{f}$	Upper bound $\overline{f}$	E (%)
f 1	${2.963510}^2$	${3.550310}^2$	$18.22$	${3.001210}^2$	${3.468610}^2$	$14.51$
f 2	${5.702510}^2$	${6.751810}^2$	17.07	${5.808310}^2$	${6.577510}^2$	$12.51$
f 3	${9.992910}^2$	${1.214010}^3$	$19.51$	${1.008910}^3$	${1.195410}^3$	$16.95$
f 4	${1.097010}^3$	${1.302110}^3$	$17.26$	${1.111810}^3$	${1.274510}^3$	$13.69$

Analyzing Table 1, one can notice the spread of the results obtained with the Interval-based Method when compared to the MC method, since E values are greater (3.6% on average). This also can be noticed in the graphs of Figure 3, where the convex hull and interval limits for the Interval-based Method encompass the corresponding values for the MC method. The numerical evaluations took about 260 seconds in an i5-3.7 GHz computer with 16 GB RAM.

4.2. Example 2 - Buckling load of a clamped plate

The second example is intended to investigate the uncertainty in the first four buckling loads in the case of a clamped plate with a one-dimensional compressive load N_x. In this case, the mass density is not assumed uncertain since it does not collaborate with the buckling load. The material properties are the same as the ones used in example 1. It is assumed the same uncertainty of ±5^o for θ at each layer and uncertainty of ±10%, (for α-cut level 0) in E₁₁, G₁₁ and ν₁₂. The dimensions for the plate are a=0.6 and b=0.6, and the laminate thickness is t=2.4×10^-2m with orthotropic cross-ply [0°/90°/0/90°]. Figure 4 is obtained for the convex hulls for the first 4 buckling loads, which are dimensionless.

Figure 4
Output uncertain intervals and convex hulls for the Interval-based procedure and MC Simulation for buckling loads (dimensionless buckling load ratio, $\lambda =N_{crit}/N_x$).

For comparisons purposes, Table 2 shows the differences between the Interval and the simple Monte Carlo results in terms of the Error values.

Table 2
Error-values for interval output variables obtained by the Interval-based Method and the MC method.

Buckling load ratio	Interval-based Method			MC Method
	Lower bound $\underset{\_}{\lambda }$	Upper bound $\overline{\lambda }$	E (%)	Lower bound $\underset{\_}{\lambda }$	Upper bound $\overline{\lambda }$	E (%)
${\lambda }_1$	${1.673010}^7$	${2.028710}^7$	$9.24$	${1.693910}^7$	${1.994910}^7$	$6.28$
${\lambda }_2$	${1.707710}^7$	${2.072910}^7$	$9.34$	${1.729010}^7$	${2.040310}^7$	$6.49$
${\lambda }_3$	${1.783210}^7$	${2.168910}^7$	$9.52$	${1.809810}^7$	${2.135510}^7$	$6.49$
${\lambda }_4$	${1.814310}^7$	${2.209910}^7$	$9.68$	${1.832810}^7$	${2.177510}^7$	$7.15$

Again, Analyzing Table 2, one can notice that the spread of the results obtained with the Interval-based Method is better when compared to the MC method since E values are greater (2.8% on average). The numerical evaluations took about 90 seconds in an i9-3.6 GHz computer with 32 GB RAM.

4.3. Example 3 - Strength safety factor using Tsai-Wu failure criterion

The third example presented in this paper is the evaluation of the Tsai Wu failure criterion. Each laminate thickness is t=4.75×10^-3m with orthotropic cross-ply [0°/90°/0/90°]. Like in example 2, mass density is not assumed uncertain since it does not contribute to the Tsai-Wu failure criterion. The material properties used are: E₁₁=181 GPa, E₂₂=10.3 GPa, G₁₂=7.17 GPa, and ν₁₂=0.28. It is assumed the same uncertainty used in example 2,±5^o for θ at each layer and uncertainty of ±10%, (for α-cut level 0) in E₁₁, G₁₁ and ν₁₂. The load configuration used is $N_x=10000$ N/m, $N_y=-10000$ N/m, $M_x=10000$ Nm/m, $M_y=-10000$Nm/m. The limit strength parameters of the material used (carbon-epoxy²⁹) are ${\left({\sigma }_1^T\right)}_{ult}=1.5$×${10}^9\text{Pa},{\left({\sigma }_1^C\right)}_{ult}=1.5$×${10}^9\text{Pa},{\left({\sigma }_2^T\right)}_{ult}=4.0$×${10}^7\text{Pa},{\left({\sigma }_2^C\right)}_{ult}=2.46$×${10}^8\text{Pa},$ and ${\left({\tau }_{12}\right)}_{ult}=6.8$×${10}^7\text{Pa}$.

Convex hulls presented in Figure 5 is obtained for the maximum and minimum safety factor, $\lambda$, found using the Tsai Wu failure criterion. Table 3 shows the differences in the Interval based method and Monte Carlo method, as well as the Error values. The time consumed in the simulation was 80 seconds in an i5-3.7 GHz computer with 16 GB RAM.

Figure 5
Output uncertain intervals and convex hull for the interval-based procedure and MC simulation for Tsai-Wu criterion, maximum and minimum$\lambda$.

Table 3
Error-values for interval output variables obtained by the Interval-based Method and the MC method.

Buckling load ratio	Interval-based Method			MC Method
	Lower bound $\underset{\_}{\lambda }$	Upper bound $\overline{\lambda }$	E (%)	Lower bound $\underset{\_}{\lambda }$	Upper bound $\overline{\lambda }$	E (%)
${\lambda }_{\text{min}}^{\text{Tsai}-\text{Wu}}$	$-{6.298310}^{-2}$	${3.900710}^{-1}$	$170.95$	${2.424510}^{-2}$	${3.548710}^{-1}$	$124.75$
${\lambda }_{\text{max}}^{\text{Tsai}-\text{Wu}}$	${2.073810}^0$	${8.148310}^0$	$91.81$	${3.012210}^0$	${7.601910}^0$	$69.37$

Note that the convex hull obtained by the method presented in this paper shows a violation of the failure criteria $\lambda =0$, while Monte Carlo (10740 samples) finds no configuration that violates this limit.

Similar to examples 1 and 2, example 3 shows that the error found in the Interval based method was 46% and 22% higher when compared to the Monte Carlo method. This indicates that the Interval-based method found wider limits for the uncertain response, or in other words, a safer and more precise approach to uncertainty evaluation.

5. CONCLUSION

It was observed that the convex hulls as the output intervals are not well-defined when using traditional Monte Carlo simulation. This is attributed to the fact that the worst scenario is a precise combination of values of uncertain input variables that simple random simulations cannot obtain. Despite the considerable amount of simulations, the extreme (worst-case scenario) is only correctly defined by the Interval based methodology proposed in this paper. In the examples analyzed, mode frequencies and buckling loads, the Interval based methodology presented 2,5% to 4,5% larger error values when compared to the Monte Carlo method, which implies a larger uncertain area. For the Tsai-Wu failure theory, the error values are even greater (22% to 46%), where the Monte Carlo method did not detect a violation of the failure limit, which was found by the Interval based approach presented in this paper.

In future studies, one can apply this analysis to more complex manufacturing processes such as ply-off laminates studied by³⁰, which are quite susceptible to imperfections during their manufacture, according to the authors. These imperfections are mainly due to the ply orientation uncertainties that are already included in the present methodology.

6. ACKNOWLEDGMENTS

The authors thank CNPq and CAPES for the partial financial support for this research.

7. References

• ¹
  He T, Liu L, Makeev A. Uncertainty analysis in composite material properties characterization using DIC and finite element model updating. Comp. Struct. 2018;184:337-51. http://dx.doi.org/10.1016/j.compstruct.2017.10.009
  » http://dx.doi.org/10.1016/j.compstruct.2017.10.009
• ²
  Dey S, Mukhopadhyay T, Adhikari S. Stochastic free vibration analysis of angle-ply composite plates – A RS-HDMR approach. Compos Struct. 2015;122:526-36. http://dx.doi.org/10.1016/j.compstruct.2014.09.057
  » http://dx.doi.org/10.1016/j.compstruct.2014.09.057
• ³
  Dey S, Mukhopadhyay T, Adhikari S. Uncertainty quantification in laminated composites: a meta-model based approach. Boca Raton: CRC Press; 2018.
• ⁴
  Sriramula S, Chryssanthopoulos MK. Quantification of uncertainty modelling in stochastic analysis of FRP composites. Composites Part A: Applied Science and Manufacturing. 2009;40(11):1673-84. https://doi.org/10.1016/j.compositesa.2009.08.020
  » https://doi.org/10.1016/j.compositesa.2009.08.020
• ⁵
  Kakur N, Krishnapillai S, Ramachandran V. Effect of fiber orientation on carbon/epoxy and glass/epoxy composites subjected to shear and bending. Diffus Defect Data Solid State Data Pt B Solid State Phenom. 2017;267:103-8. http://dx.doi.org/10.4028/www.scientific.net/SSP.267.103
  » http://dx.doi.org/10.4028/www.scientific.net/SSP.267.103
• ⁶
  Wehrle EJ. Design optimization of lightweight space-frame structures considering crashworthiness and parameter uncertainty [dissertation]. Munich: Technical University of Munich; 2015.
• ⁷
  Rao SS, Alazwari MA. Failure modeling and analysis of composite laminates: interval-based approaches. J Reinf Plast Compos. 2020;39(21-22):817-36. http://dx.doi.org/10.1177/0731684420932645
  » http://dx.doi.org/10.1177/0731684420932645
• ⁸
  Barber CB, Dobkin DP, Huhdanpaa H. The quickhull algorithm for convex hulls. ACM Trans Math Softw. 1996;22(4):469-83. http://dx.doi.org/10.1145/235815.235821
  » http://dx.doi.org/10.1145/235815.235821
• ⁹
  Faes M, Moens D. Identification and quantification of spatial interval uncertainty in numerical models. Comp. Struct. 2017;192:16-33. http://dx.doi.org/10.1016/ j.compstruc.2017.07.006
  » http://dx.doi.org/10.1016/ j.compstruc.2017.07.006
• ¹⁰
  Venini P, Mariani C. Free vibrations of uncertain composite plates via stochastic Rayleigh-Ritz approach. Computers & Structures. 1997;64(1-4):407-23. https://doi.org/10.1016/S0045-7949(96)00161-7
  » https://doi.org/10.1016/S0045-7949(96)00161-7
• ¹¹
  Liu WK, Mani A, Belytschko T. Finite element methods in probabilistic mechanics. Probab Eng Mech. 1987;2:201-13. http://dx.doi.org/10.1016/0266-8920(87)90010-5
  » http://dx.doi.org/10.1016/0266-8920(87)90010-5
• ¹²
  Cherng RH, Wen YK. Stochastic finite element analysis of uncertain nonlinear plane trusses under random excitation [dissertation]. Urbana-Champaign, IL: Civil Engineering Studies, University of Illinois; 1992.
• ¹³
  Wall FJ, Bucher CG. Sensitivity of expected exceedance rate of SDOF system response to statistical uncertainties of loading andsystem parameters. Probab Enann Mech. 1987;2:138-46. http://dx.doi.org/10.1016/0266-8920(87)90004-X
  » http://dx.doi.org/10.1016/0266-8920(87)90004-X
• ¹⁴
  Ma D, Amico S, Giglio M, Manes A. Effect of fibre bundle uncertainty on the tensile and shear behaviour of plain-woven composites. Compos Struct. 2020;259. http://dx.doi.org/10.1016/j.compstruct.2020.113440
  » http://dx.doi.org/10.1016/j.compstruct.2020.113440
• ¹⁵
  Kim T, Hwang I, Sin H. Optimal design of composite laminate with uncertainty in loading and material properties considered. ASME AMD. 2005;256:557-64.
• ¹⁶
  Lee J, Haftka RT, Griffin OH. Detecting delaminations in a composite beam using anti-optimization. Struct Optim. 1994;8:93-100. http://dx.doi.org/10.1007/BF01743304
  » http://dx.doi.org/10.1007/BF01743304
• ¹⁷
  Sarp A, Francoise L, Georges D, Vincent C. Optimization of Laminated Composites Under Buckling Uncertainties via Anti-Optimization. In: 9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization ; Atlanta. Atlanta: AIAA/ISSMO; 2002. https://doi.org/10.2514/6.2002-5417
  » https://doi.org/10.2514/6.2002-5417
• ¹⁸
  Lombardi M, Haftka RT. Anti-optimization technique for structural design under load uncertainties. Comput Methods Appl Mech Eng. 1998;157(1-2):19-31. http://dx.doi.org/10.1016/S0045-7825(97)00148-5.
  » https://doi.org/. http://dx.doi.org/10.1016/S0045-7825(97)00148-5
• ¹⁹
  Alfredo F. Buckling optimization and antioptimization of composite plates: uncertain loading combinations. Int J Numer Methods Eng. 2002;53(3):719-32. http://dx.doi.org/10.1002/nme.309
  » http://dx.doi.org/10.1002/nme.309
• ²⁰
  Chen X, Wang X, Qiu Z, Wang L, Li X, Shi Q. A novel reliability-based two-level optimization method for composite laminated structures. Compos Struct. 2018;192:336-46. http://dx.doi.org/10.1016/j.compstruct.2018.03.016
  » http://dx.doi.org/10.1016/j.compstruct.2018.03.016
• ²¹
  Chen X, Wang X, Qiu Z. A novel uncertainty analysis method for composite structures with mixed uncertainties including random and Interval variables. Compos Struct. 2018;184:400-10. http://dx.doi.org/10.1016/j.compstruct.2017.09.068
  » http://dx.doi.org/10.1016/j.compstruct.2017.09.068
• ²²
  Peng X, Li D, Wu H, Liu Z, Li J, Jiang S, et al. Uncertainty analysis of composite laminated plate with data-driven polynomial chaos expansion method under insufficient input data of uncertain parameters. Compos Struct. 2019;209:625-33. http://dx.doi.org/10.1016/j.compstruct.2018.11.015
  » http://dx.doi.org/10.1016/j.compstruct.2018.11.015
• ²³
  Chen G, Wang T, Lu C, Yang Y, Li L, Yin Z, et al. Uncertainty representation of natural frequency for laminated composite cylindrical shells considering probabilistic and interval variables. Appl Sci. 2021;11(4):1883. http://dx.doi.org/10.3390/app11041883
  » http://dx.doi.org/10.3390/app11041883
• ²⁴
  Jin JW, Jeon BW, Choi CW, Kang KW. Multi-scale probabilistic analysis for the mechanical properties of plain weave carbon/epoxy composites using the homogenization technique. Appl Sci. 2020;10:6542. http://dx.doi.org/10.3390/app10186542
  » http://dx.doi.org/10.3390/app10186542
• ²⁵
  Shi D, Teng X, Jing X, Lyu S, Yang X. A multi-scale stochastic model for damage analysis and performance dispersion study of a 2.5D fiber-reinforced ceramic matrix composites. Compos Struct. 2020;248:112549. http://dx.doi.org/10.1016/j.compstruct.2020.112549
  » http://dx.doi.org/10.1016/j.compstruct.2020.112549
• ²⁶
  Huang T, Gao J, Sun Q, Zeng D, Su X, Liu WK, et al. Stochastic nonlinear analysis of unidirectional fiber composites using image-based microstructural uncertainty quantification. Compos Struct. 2021;260:113470. http://dx.doi.org/10.1016/j.compstruct.2020.113470
  » http://dx.doi.org/10.1016/j.compstruct.2020.113470
• ²⁷
  Luo Y, Hong M, Liu Y. Analytical solutions to the fundamental frequency of arbitrary laminated plates under various boundary conditions. J Marine Sc Appl. 2015;14:46-52. http://dx.doi.org/10.1007/s11804-015-1294-x
  » http://dx.doi.org/10.1007/s11804-015-1294-x
• ²⁸
  Tsai SW, Wu EM. A General theory of strength for anisotropic materials. J Compos Mater. 1971;5(1):58-80. http://dx.doi.org/10.1177/002199837100500106
  » http://dx.doi.org/10.1177/002199837100500106
• ²⁹
  Kaw AK. Mechanics of composite materials. 2nd ed. Boca Raton: Taylor & Francis; 2006.
• ³⁰
  Méndez DC, Almeida Y, Cunha F Jr, Gomes SG. Optimum design of composite structures with ply drop-offs using response surface methodology. Eng Comp. 2021;[ahead-of-print]. https://doi.org/10.1108/EC-07-2020-0354
  » https://doi.org/10.1108/EC-07-2020-0354

Publication Dates

• Publication in this collection
  22 Oct 2021
• Date of issue
  2021

History

• Received
  02 July 2021
• Reviewed
  10 Sept 2021
• Accepted
  19 Sept 2021