Comparison between a linear and cubic fiber angle variation on buckling response of variable angle tow composite panels

Abstract

Advances in manufacturing techniques have allowed more flexibility to the design and new possibilities to apply composites materials in lightweight structures. Novel techniques such as the automated fiber placement allow the fibers to follow curvilinear paths, making possible laminate properties that vary within the laminate plane. These types of laminates are known as variable stiffness laminates or variable angle tow. In this work, the maximum critical buckling load of composite panels with variable stiffness through a spatially varying fiber orientation has been analyzed for two different boundary conditions. This works compares the outcomes of in-plane stress and critical buckling load for linear and cubic fiber angle considering four aspect ratios. Manufacturing constraint has been considered in the analysis of the laminates. The finite element method has been applied to solve the system elliptic partial differential equations that govern the in-plane behavior of these panels. The Ritz method has been used to find the buckling loads for the variable stiffness panels. Results for four different aspects ratios are presented. Improvements in the buckling load of up to 18% for cubic fiber angle variation over linear fiber angle variation were found.

Keywords:
Composite material; variable angle tow composites; critical buckling load; in-plane analysis; automated fiber placement

Graphical abstract

1 INTRODUCTION

Composite materials have been widely used in many industries in the last years. The aerospace industry is the one that most uses this material. It is said that the composites are one of the most important classes of engineering materials, being the second one, only behind the steel in industrial importance and range of applications (Clyne and Hull, 2019). Modern aircraft, such as Boeing 787 and Airbus A350 XWB, have been successfully applied advanced composite materials. The last one, for instance, has predominantly a composite airframe with 53% of its structure weight manufactured of composite materials, distributed over the wings, fuselage and other primary aerostructures (Marsh, 2010).

The growing in composite applications motivated by commercial interests has encouraged significant developments in manufacturing processes and structures design. The automated fiber placement (AFP) machines are a technology that has been developed over the past thirty years to attempt the industrial demands for fast manufacturing processes, as well as great flexibility to build composite laminates with complex geometries (Nik et al., 2014). The development of this new technology provided novel design approaches such as variable stiffness laminates. Unlike traditional composite laminates with unidirectional fibers orientation, variable angle tow (VAT) composite laminates have a continuously varying fiber orientation, yielding a variable stiffness laminate. The AFP machines are able to dispose the prepreg fibers following a specified path set by the designer (Brooks and Martins, 2018). These laminates expand the design space to create more adaptable and efficient structures. It has been demonstrated an improvement in the critical buckling load and strength using VAT composites in comparison with conventional composites (Setoodeh et al., 2009a; Wu et al., 2013; Gurdal et al., 2008; IJsselmuiden et al., 2010; Li et al., 2002). The initial work on variable stiffness composites was made by Gurdal and Olmedo (1993). Using a linear fiber angle variable, they showed that improvements on axial stiffness is possible for a laminate subject to uniform end shortening.

Li et al. (2002) applied the fiber steering concept to increase the performance of the bearing strength of bolted holes, and the fibers in this study were steered such that the fibers followed the principle tensile and compressive stress trajectories around the hole. A bolted joint in composite coupons has been strengthened by up to 169% for peak load or 36% for bearing strength using the modified steered pattern. Gurdal et al. (2008) parameterized fiber paths using linear fiber angle variation to find the best configuration regarding the critical buckling load. Their results showed better performance comparing to straight fiber composites. Tatting and Gurdal (2002) designed and manufactured straight and linear fiber angle variation composites plates with and without a central hole. The VAT composite plates showed an improvement of 64.15% and 53.9% on the critical buckling load over the straight fiber plates, respectively for plates with and without a hole. Marouene et al. (2016) also manufactured and tested VAT composite plates with linear fiber angle variation which presented critical buckling load 45% higher as compared to quasi-isotropic straight fiber composite plates. Wu et al. (2012) performed an optimization to enhance the critical buckling load of a composite plate using Lagrangian polynomials to describe the fiber angle variation. Their results showed an increase of 255% increase on the critical buckling over quasi-isotropic laminates.

Huang and Haftka (2005) have shown that increasing the load carrying capability of a composite plate with a hole in the center and loaded in tension results in a set of nearly concentric circles near the hole where the authors chose to limit the tailored region of the composite with a hole to the region close to the hole. Buckling improvement of a component or structure using a variable stiffness composite has been the topic of several research efforts. Although significant increase in the buckling load can be obtained through tailoring laminate stacking sequence of composite panels with traditional straight fibers, the potential of fibrous composites is not fully exploited Setoodeh et al. (2009b). The improvement is attributed to the redistribution of in-plane loads to relatively stiff regions, and then resist buckling in critical regions Hao et al. (2017b). An approach for the design of panels with cutouts was first introduced in the late eighties by Hyer and Charette (1991) in which curvilinear fibers were suggested to improve structural response instead of straight fiber paths.

Setoodeh et al. (2009b) showed that variable stiffness composite is able to withstand more than twice the compressive load before buckling occurs compared to a quasi-isotropic composite panel. The paper has allowed for a better understanding of the load redistribution mechanism responsible for an increased buckling load. A buckling load increase of 33.9% with respect to a constant stiffness composite panel is reported in Lopes et al. (2008) for a variable stiffness composite panel where overlapping tows are present. Similar results for a clamped square variable stiffness composite panel were reported where an improvement of 66% in the buckling load is found for a simply supported composite panel.

Van den Brink et al. (2012) demonstrate that an improvement of buckling load ranging from 35% to 67% can be achieved by using a variable-stiffness design. Hao et al. (2017a) analyzed panel with multiple cutouts using a variable stiffness based on flow field function. They showed that improvements on buckling load is obtained when compared with linear angle variations and straight fibers.

The aim of this work is to compare the maximum critical buckling load VAT laminated plates using a cubic fiber angle variation considering the minimum radius manufacturing constraint. A numerical model based Finite Element Method (FEM) and the Rayleigh-Ritz method formulation is developed to calculate the critical buckling load. The results are then compared to linear fiber angle variation for composite plates.

2 VAT LAMINATES

Unlike of traditional straight fiber composite whose fibers have the same angle orientation on the entire lamina, in the VAT concept, the fiber orientation varies over the lamina based in the reference axis. The fiber orientation can vary along one or two coordinates, depending on the design criteria. Several methods have been proposed to model the fiber orientation. The linear variation is the simplest case as proposed by Gurdal et al. (2008) according to Eq. 1,

$\theta \left(x^{\prime }\right)=\varphi +\left(T_1-T_0\right)\frac{\left|x^{\prime }\right|}{d}+T_0$ (1)

where $\theta$ is the fiber orientation, $\varphi$ is the rotation angle which defines the coordinate system $x^{\prime }-y^{\prime }$, $T_1$ and $T_0$ are the angles at panel edges and center, respectively, and $d$ is a characteristic distance (Fig. 1). The layer representation for linear fiber angle variation is specified by $\varphi T_0|T_1$. There are three parameters to construct the laminates, $\varphi$, $T_0$ and $T_1$. From these parameters is possible to generate laminates with a wide range of tailoring possibilities.

It is also possible to model more complex fiber angle orientations such as non-linear variations with different polynomial degrees, providing high design flexibility, but it has been shown that trajectories represented by high-orders polynomials can be difficult, or even impossible, to manufacture due to limitations of the manufacturing process, especially regarding to the small radii curvature. To avoid possible manufacturing defects, one can make a cubic fiber angle orientation. A cubic fiber angle orientation can be defined as:

$\theta \left(x^{\prime }\right)=\varphi +{\left[\left(c^{-1}a\right)\right]}^T{\left|x\right|}^{\text{'}}$ (2)

Figure 1:
Path definition of a linear angle variation. Adapted from Gurdal et al. (2008).

$c=\left[\begin{array}{cccc}0& 0& 0& 1\\ d^3/216& d^2/36& d/6& 1\\ d^3/27& d^2/9& d/3& 1\\ d^3/8& d^2/4& d/2& 1\end{array}\right],a={\left\{T_0,T_1,T_2,T_3\right\}}^{T}andx={\left\{x^3,x^2,x^1,x^0\right\}}^T$ (3)

where $c$ is the matrix of coefficients $a$ is the angle orientation vector and $x$ is the variable vector related to the $x^{\prime }-y^{\prime }$ coordinate system rotated by an angle $\varphi$. The layer representation is defined as $\varphi \langle T_0\text{|}{T}_1\text{|}{T}_2\text{|}{T}_3\rangle$, similar to the linear angle variation. For $\varphi =0$, $T_0$, $T_1$, $T_2$, $T_3$ are located at $x=0$, $x=a/6$, $x=a/3$, and $x=a/2$, respectively, as shown in Fig. 2. Once is defined the fiber angle orientation, the reference path of steered tow over the plate can be obtained by the following relation:

$\frac{d{y}^{\text{'}}}{d{x}^{\text{'}}}=\tan \theta \left(x^{\prime }\right)$ (4)

Figure 2:
VAT composite plates with $\varphi =0°$. (a) linear angle orientation and (b) cubic angle orientation.

The first fiber path created is the reference fiber path. The next fiber path is made by shifting the reference fiber path a fixed amount in the y direction, for example, when $\varphi =0°$. The remaining fiber paths for the lamina are made in the same manner, with the only difference being the amount that each path is shifted along the shifting axis. This process is followed until all panel to be fully filled, as shown in Fig. 3, considering a linear angle variation and $\varphi =0°$. The same idea can be extended to cubic angle variations and other values of $\varphi$. A 3D schematic view of a VAT composite plate is show in Fig. 4.

Figure 3:
Shifting of reference path over a single lamina surface for a linear angle variation and $\varphi =0°$.

Figure 4:
3D representation of a VAT composite plate and lamina.

The manufacturing constraint of variable stiffness is related to the minimum steering radius that the AFP machine is capable of performing. Defining the path reference, the steering radius can be obtained by inverting the tow path curvature $\kappa$, as defined in Eq. 5.

$r=\frac{1}{\kappa }$ (5)

3 IN-PLANE AND BUCKLING ANALYSIS

The governing equations of the in-plane response of a laminated plates are expressed as a system of differential equations as follows:

$\begin{array}{l}\frac{\partial N_x}{\partial x}+\frac{\partial N_{xy}}{\partial y}=0\\ \frac{\partial N_y}{\partial y}+\frac{\partial N_{xy}}{\partial x}=0\end{array}$ (6)

where $N_x$, $N_y$ and $N_{xy}$ are the in-plane stress resultants. The constitutive relations based on classical theory for thin laminates are (Jones, 1998):

$\left\{\begin{array}{c}N\\ M\end{array}\right\}=\left[\begin{array}{cc}A& B\\ B& D\end{array}\right]\left\{\begin{array}{c}{\int }^0\\ {\kappa }^0\end{array}\right\}$ (7)

where $N$ and $M$ are vector and represent the stress resultant forces and moments, and ${ϵ}^0$ and ${\kappa }^0$ denotes the mid-plane strains and curvatures, respectively. $A$ is the in-plane matrix, $B$ is the coupling matrix and $D$ is the bending matrix. One can notice that for variable angle tow composites these matrices are functions of the panel spatial coordinates. These matrices are defined through the thickness integration

$\left\{A_{ij},B_{ij},D_{ij}\right\}=\underset{-h/2}{\overset{h/2}{{\displaystyle \int }}}{\overline{Q}}_{ij}\left\{\mathrm{1,}z,z^2\right\}dz$ (8)

where $h$ is the panel thickness and ${\overline{Q}}_{ij}$ are the reduced transformed stiffness terms. They are function of panel coordinates and given by

$\begin{array}{l}{\overline{Q}}_{11}=U_1+U_2\cos \left[2\theta \left(x,y\right)\right]+U_3\cos \left[4\theta \left(x,y\right)\right]\\ {\overline{Q}}_{12}=U_4-U_3\cos \left[4\theta \left(x,y\right)\right]\\ {\overline{Q}}_{22}=U_1-U_2\cos \left[2\theta \left(x,y\right)\right]+U_3\cos \left[4\theta \left(x,y\right)\right]\\ {\overline{Q}}_{66}=U_5-U_3\cos \left[4\theta \left(x,y\right)\right]\\ {\overline{Q}}_{16}=-\left(1/2\right)U_2\sin \left[2\theta \left(x,y\right)\right]-U_3\sin \left[4\theta \left(x,y\right)\right]\\ {\overline{Q}}_{26}=-\left(1/2\right)U_2\sin \left[2\theta \left(x,y\right)\right]+U_3\sin \left[4\theta \left(x,y\right)\right]\end{array}$ (9)

where $U_i$ are the invariant properties of composite materials. Considering laminates to be mid-plane symmetric, the matrix $B_{ij}$ is null. Thus, $N$ and $M$ can be decoupled from Eq. 7:

$\left\{\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\end{array}\right\}=\left[\begin{array}{ccc}{A}_{11}\left(x,y\right)& A_{12}\left(x,y\right)& A_{16}\left(x,y\right)\\ A_{12}\left(x,y\right)& A_{22}\left(x,y\right)& A_{26}\left(x,y\right)\\ A_{16}\left(x,y\right)& A_{26}\left(x,y\right)& A_{66}\left(x,y\right)\end{array}\right]\left\{\begin{array}{c}\begin{array}{c}{ϵ}_x^0\\ {ϵ}_y^0\end{array}\\ {\gamma }_{xy}^0\end{array}\right\}$ (10)

$\left\{\begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\right\}=\left[\begin{array}{ccc}{D}_{11}\left(x,y\right)& D_{12}\left(x,y\right)& D_{16}\left(x,y\right)\\ D_{12}\left(x,y\right)& D_{22}\left(x,y\right)& D_{26}\left(x,y\right)\\ D_{16}\left(x,y\right)& D_{26}\left(x,y\right)& D_{66}\left(x,y\right)\end{array}\right]\left\{\begin{array}{c}\begin{array}{c}{\kappa }_x^0\\ {\kappa }_y^0\end{array}\\ {\kappa }_{xy}^0\end{array}\right\}$ (11)

Considering the linear strain displacement relations (Eq. 12) and Eq. 10, Eq. 1 can be written as shown in Eq. 13

${ϵ}_x^0=\frac{\partial u}{\partial x},{ϵ}_y^0=\frac{\partial u}{\partial y},{\gamma }_{xy}^0=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}$ (12)

$\begin{array}{l}{A}_{11}\left(x,y\right)\frac{{\partial }^{2}u}{\partial x^2}+A_{66}\left(x,y\right)\frac{{\partial }^{2}u}{\partial y^2}+\left[A_{12}\left(x,y\right)+A_{66}\left(x,y\right)\right]\frac{{\partial }^{2}v}{\partial x\partial y}\\ +\frac{\partial A_{11}\left(x,y\right)}{\partial x}+\frac{\partial A_{12}\left(x,y\right)}{\partial x}\frac{\partial v}{\partial y}=0\\ A_{66}\left(x,y\right)\frac{{\partial }^{2}v}{\partial x^2}+A_{22}\left(x,y\right)\frac{{\partial }^{2}v}{\partial y^2}+\left[A_{12}\left(x,y\right)+A_{66}\left(x,y\right)\right]\frac{{\partial }^{2}u}{\partial x\partial y}\\ +\frac{\partial A_{66}\left(x,y\right)}{\partial x}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)=0\end{array}$ (13)

These are elliptic partial differential equations (PDEs), coupled in terms of the axial displacement $u$ and transverse displacement $v$. Furthermore, these equations have variable coefficients since the elements of $A_{ij}\left(x,y\right)$ are functions of both the $x$ and $y$. The solution of these equations is made simultaneously yielding the $u$ and $v$ displacement field. Further details of the derivation above can be seen in Whitney (1987) and Sadd (2009). Once $u$ and $v$ displacement field is calculated, the in-plane strains and stresses resultants can be calculated using Eqs. 12 and 10, respectively. In order to solve the buckling problem, it is possible to apply the Ritz method. This method is an energy-based method and it is applied to find the lowest load that will cause buckling (Jones, 1998; Kollar and Springer, 2003; Kumar, 2018). For rectangular plate simply supported along its four edges, the plate deflection expression that satisfies the geometrical boundaries is as follows:

$w\left(x,y\right)={\displaystyle \sum }_{i=1}^I{\displaystyle \sum }_{j=1}^J{C}_{ij}\sin \left(\frac{i\pi x}{a}\right)\sin \left(\frac{j\pi y}{b}\right)$ (14)

where $C_{ij}$ are arbitrary coefficients and $I$ and $J$ are integers indicating the highest number of terms in the $x$ and $y$ direction. According to principle of stationary potential energy, at equilibrium the potential energy, ${\pi }_p$ , must satisfy the following condition

$\frac{\partial {\pi }_p}{\partial C_{ij}}=\frac{\partial \left(U+V\right)}{\partial C_{ij}}=0$ (15)

where $U$ is the strain energy due to bending, and $V$ the potential energy of the in-plane loads. the strain energy of a thin composite panel is given by (Jones, 1998):

$U=\frac{1}{2}\iint \left[D_{11}{w}_{,xx}^2+2D_{12}{w}_{,xx}{w}_{,yy}+D_{22}{w}_{,yy}^2+4\left(D_{16}{w}_{,xx}^2+D_{26}{w}_{,yy}^2\right)+D_{22}{w}_{,xy}^2\right]dxdy$ (16)

where $w$ is a function that represent the plate deflection. The potential energy is related to the applied load in the pre-buckled state and the mid-plane strains that is the effect from the out-plane deflections of the panel. The $V$ expression is given by (Jones, 1998):

$V=\frac{1}{2}\lambda \iint \left(N_x{w}_{,xx}+N_y{w}_{,yy}+2N_{xy}{w}_{,xy}\right)dxdy$ (17)

where $\lambda$ is an arbitrary multiplier. Equation 14 is then substituted into Eqs. 16 and 17. The differentiations are performed and Eq. 15 reduces to

${\displaystyle \sum }_{i=1}^I{\displaystyle \sum }_{j=1}^J\left[{\pi }^4{I}_1+\frac{4{\pi }^{4}ijkl}{a^2{b}^2}{I}_2-2{\pi }^4\left(I_3+I_4\right)+\frac{\lambda {\pi }^{2}ik}{a^2}{I}_5+\frac{\lambda {\pi }^{2}jl}{b^2}{I}_6+\frac{\lambda {\pi }^2}{ab}\left(jk{I}_7+il{I}_8\right)\right]C_{ij}=0$ (18)

where $\overline{x}=\pi x/a$, $\overline{y}=\pi y/b$, the indexes $k$ and $l$ take the values $k=\mathrm{1,}\mathrm{2,}\cdots ,I$ and $l=\mathrm{1,}\mathrm{2,}\cdots ,J$. The terms $I_1$, $I_2$, $I_3$, $I_4$, $I_5$, $I_6$, $I_7$ and $I_8$ are shown in Appendix. Equation 18 can be rewritten as an eigenvalue problem:

$\left(K-\lambda M\right)C_{ij}=0$ (19)

where $K$ is the stiffness matrix and $M$ is the geometric matrix. The lowest value of $\lambda$, ${\lambda }_{cr}$, from Eq. 19 determines the critical buckling load or displacement, as it was analyzed in this work.

4 METHODOLOGY

In this work, it was analyzed panels subject to an end shortening displacement with two different boundary conditions in the edges. Uniform end shortening is a common boundary condition under testing conditions. In this sort of testing, a uniform end displacement is applied while the equivalent load is measured. A schematic representation can be seen in Fig. 5. The Case I consists of an applied end shortening, $u_0$, at $x=\pm a/2$. The transverse edges have no applied forces and they are not free to deform, $v\left(x,y\right)=0$. In the Case II, in the same manner as in Case I, an end shortening is applied, but other edges are free to deform. In both cases, it is not considered a shear stress applied in any edge.

Figure 5:
Panel boundary conditions for Case I and Case II.

To model the fiber angle variation, it was considered a linear and a cubic angle variation. Regarding the angle $\varphi$, it was considered in this work two different values: 0 and 90 deg. Thus, Eqs. 1 and 3 only vary along one direction, $x$ or $y$, depending on the angle. When $\varphi =0°$ , the angle varies on x-axis, $\theta \left(x\right)$. When $\varphi =90°$ , the angle varies on y-axis, $\theta \left(y\right)$. Besides that, the radii of curvature were limited to values less than 635 mm, to avoid manufacturing defects (Tatting and Gurdal, 2002; Marouene et al., 2016). The plate dimensions $a$ and $b$ were considered 1 m and 1 m, respectively, when aspect ratio $a/b$ equals 1. For different aspect ratios analyzed, the dimension $b$ was fixed while dimension $a$ varies.

To calculate the in-plane response of the VAT panels, Finite Element Method (FEM) was applied. The FEM is a technique to solve PDEs numerically. To solve PDEs with the FEM, basically it is needed a discrete representation of a region, known as a mesh, PDEs that model the physical behavior of the phenomena and boundary conditions that coupling the PDEs with the region (Hughes, 2000).

In this work, the software Mathematica (Wolfram et al., 1999) was used to solve the PDEs (Eq. 13). It has a built-in function called NDSolve, which is a powerful solver of differential equations (Mikhailov, 2008). As one of its algorithms, NDSolve has an algorithm based on finite element method. As options, it is possible to choose the element type (MeshElementType) and the element size (MaxElementMeasure). 8-node second-order rectangular element was chosen (Oñate, 2013). The node degrees of freedom are displacement in $x$ and y, displacements $u$ and $v$, respectively.

A mesh convergence was performed for both linear and cubic angle variations to define the number of elements, as shown in Figs. 6 and 7. It can be seen that in the case where the fiber varies linearly, it is needed less elements than a cubic variation. This behavior is expected knowing that a cubic variation is more complex and presents larger displacement gradients. Case II is not shown here since it presents convergence using few elements. Regarding the aspect ratio, it does not show relevant differences between the values 0.5, 1, 2, and 3. For the linear variation, 200 elements were used, and for cubic variation 400 elements. The laminate ${\left[0°+\langle 0°\text{|}90°\rangle \right]}_{3S}$ was evaluated for a linear variation because it presents the biggest displacement gradient and this way the convergence study is valid for angles between 0 and 90 deg. The same idea can be extended to the laminate ${\left[0°+\langle 0°\text{|}90°\text{|}0°\text{|}90°\rangle \right]}_{3S}$ in the cubic variation case.

Regarding the buckling analysis, in order to get reliable results, it is essential to determine the number of terms necessary in the $w\left(x,y\right)$ function. As it was analyzed the buckling load for different aspect ratios (0.5, 1, 2 and 3), it was required to perform a convergence study for all aspect ratios. The convergence of the eigenvalues for Case I for linear and cubic angle variation is shown in Tables 1 and 2, respectively. As the boundary conditions of the Case I involve $N_x$ and $N_y$ in the calculations, it shows most difficulty to converge. Thus, it was only evaluated the convergence for Case I. Regarding the angle $\varphi$, the angles 0 and 90 deg exhibited similar results, so only the results for $\varphi =0°$ are shown. The same happened to others layups. In light of the results, it was adopted $I=J=10$.

Figure 6:
Case I $\left(a/b\right)=1$ with linear variation of the laminate ${\left[0°+\langle 0°\text{|}90°\rangle \right]}_{3S}$.

Figure 7:
Case I $\left(a/b\right)=1$ with linear variation of the laminate ${\left[0°+\langle 0°\text{|}90°\text{|}0°\text{|}90°\rangle \right]}_{3S}$.

Table 1:
Ritz method converge for linear angle variation.

Case I - ${\left[0°\pm \langle 45°\text{|}75°\rangle \right]}_{3S}$
$a/b=0.5$	$a/b=1$	$a/b=2$	$a/b=3$
$u_{cr}$	$u_{cr}$	$u_{cr}$	$u_{cr}$	$I$	$J$	$I\times J$
0.00547088	0.00600322	0.00851569	0.0115659	3	3	9
0.00544822	0.00595684	0.00847652	0.0115544	5	5	25
0.0054434	0.00594918	0.00846769	0.0115511	7	7	49
0.00544329	0.00594825	0.00846567	0.0115501	10	10	100
0.00544305	0.00594791	0.00846501	0.0115496	12	12	144
0.00544302	0.0059478	0.00846472	0.0115493	15	15	225

Table 2:
Ritz method converge for cubic angle variation.

Case I - ${\left[0°\pm \langle 45°\text{|}60°\text{|}75°\text{|}75°\rangle \right]}_{3S}$
$a/b=0.5$	$a/b=1$	$a/b=2$	$a/b=3$
$u_{cr}$	$u_{cr}$	$u_{cr}$	$u_{cr}$	$I$	$J$	$I\times J$
0.00527521	0.00649227	0.0112305	0.0162436	3	3	9
0.00526544	0.00646149	0.0112295	0.0160661	5	5	25
0.00526105	0.00645667	0.0112132	0.0160133	7	7	49
0.00526036	0.00645509	0.0112109	0.0160122	10	10	100
0.00526036	0.00645501	0.0112108	0.0160120	12	12	144
0.00526033	0.00645494	0.0112108	0.0160120	15	15	225

Having in mind that the panels were loaded through a constant displacement at its edges, it is important to define the force necessary to create the uniform displacement. Thus, the total force induced on the panel can be defined as:

$F_x=\underset{0}{\overset{b}{{\displaystyle \int }}}{N}_x\left(a,y\right)dy$ (20)

To obtain an $N_x$ average value, $F_x$ is divided by the length of the panel edge, as follows:

$N_x^{av}=\left(\frac{1}{b}\right)F_x$ (21)

In order to facilitate the comparison between the results, $N_x^{av}$ was normalized as

${\overline{N}}_x^{av}=\frac{N_x^{av}{a}^2}{E_1{h}^3}$ (22)

Finally, the critical normalized average load is obtained by

$N_{cr}={\lambda }_{cr}{\overline{N}}_x^{av}$ (23)

where ${\lambda }_{cr}$ is the lowest eigenvalue of the Eq. 19. To have a parameter related to stiffness, an equivalent panel stiffness was defined to consider an overall panel stiffness

$E_x^{eq}=\frac{F_{x}a}{hb{u}_0}$ (24)

where $h$ is the panel thickness. Substituting Eq. 20 into Eq. 24 yields

$E_x^{eq}=\left(\frac{a}{hb{u}_0}\right)\underset{0}{\overset{b}{{\displaystyle \int }}}{N}_x\left(a,y\right)dy$ (25)

5 RESULTS AND DISCUSSION

The results here are based on a twelve-layer laminate, ${\left[\pm \theta \right]}_{3S}$. The properties of the material are shown in Tab. 3. In order to compare the results from linear and cubic variation, and were normalized as

${\overline{N}}_i=N_i\frac{a}{u_0{A}_{11}^{0°}},i=x,y$ (26)

5.1 Model validation

The model was validated for linear angle variation through both experimental and numerical studies. The work of Marouene et al. (2016) was selected to validate the model. Marouene et al. (2016) manufactured two types of variable-stiffness panels: one variable-stiffness with overlaps and another with gaps, both considering a linear angle variation, shown in Fig. 8. A schematic representation of gaps and overlaps are shown in Fig. 9 The panel's mechanical properties are shown in Tab. 4. The layup and geometric features are shown in Tab. 5. The fiber-steered panels were manufactured using a VIPER AFP machine. To avoid wrinkles and micro-buckling, the authors considered a minimum radius curvature of 635 mm, based on manufacturer requirements. Thus, the same value considered in this work. In this work, it was not considered the thickness variation over the plate, it was considered the panels with constant thickness. The panels shown in Fig. 8 were loaded at horizontal top edge by a uniaxial compression by applying a uniform displacement and the other edge was fixed. More details can be found in Marouene et al. (2016). The results are shown in Tab. 6.

Table 3:
T300/N5208 Carbon/Epoxy properties.

Material Properties	T300/N5208
Longitudinal modulus	180 GPa
Transverse modulus	10.3 GPa
In-plane shear modulus	7.17
Major Poisson’s ratio	0.28

Table 4:
G40-800/5276-1 unidirectional carbon/epoxy prepreg Marouene et al. (2016).

Moduli Properties	G40-800/5276-1
Longitudinal modulus	142.7 GPa
Transverse modulus	9.1 GPa
In-plane shear modulus	4.82 GPa
Major Poisson’s ratio	0.3

Figure 8:
Specimens tested experimentally. Adapted from Marouene et al. (2016).

Figure 9:
Schematic of gaps and overlaps within a ply.

Table 5:
Panel features.

Layup	${\left[\pm \langle 49°\text{|}41°\rangle /\langle 48°\text{|}61°\rangle /\langle 57°\text{|}73°\rangle /\langle 72°\text{|}77°\rangle \right]}_{3S}$
Thickness (average)	Overlap = 2.84 mm and gap = 2.09 mm
Dimensions	a = 254 mm and b = 406 mm

Table 6:
Comparison between the results from Marouene et al. (2016) and the present work.

	Overlap	Gap
Marouene et al. (2016)	18.87 kN	15.06 kN
Present	20.81 kN	8.69 kN
Difference	8.84%	28.71%

The case of panels with gaps, the present model underestimated the result, showing a difference of 28.71% from the experimental result. In fact, to capture all effects involved on this kind of experiment, it is necessary a non-linear formulation, while the present uses a linear formulation. On the other hand, for the case with overlaps, the model overestimated, presenting a difference of 8.84%. These results shown that the gaps present an important role in the results, while overlaps not so much. As the present model considers a constant thickness variation, it approximates to the panel with overlaps. Thus, considering the difference of 8.84%, the model has a good agreement with the experimental results performed by Marouene et al. (2016).

In order to validate the model numerically, the model was compared through the results from Gurdal et al. (2008). Gurdal et al. (2008) evaluated laminates considering two boundary conditions, as evaluated in the work: one in that the transverse edges are restrained (Case I) and another with the transverse edges free to deform (Case II). It was calculated the values of $N_{cr}$ for a family of curves corresponding to various values of $T_0$ from $0°$ to $90°$ with increments of $10°$. Each curve was generated by varying the value of $T_1$ between $0°$ and $90°$ for a given value of $T_0$. The results for the Case I are shown in Fig. 10. For this case, $\varphi =0°$.

According to results calculated, the maximum value of $N_{cr}$ is 1.43, a difference of 0.38% from the results of Gurdal et al. (2008), for the laminate ${\left[0°\pm \langle 0°\text{|}50°\rangle \right]}_{3S}$. The results for the case II are shown in Fig. 11. The Case II was evaluated for $\varphi =90°$. The maximum value of $N_{cr}$, equals to 3.16, happened for the laminate ${\left[0°\pm \langle 0°\text{|}75°\rangle \right]}_{3S}$. This value is only 0.64% higher than the value of 3.14 from the results of Gurdal et al. (2008). In light of the results, one can consider that the model is reliable.

Figure 10:
Buckling and equivalent stiffness of various laminates with transverse edges restrained.

Figure 11:
Buckling and equivalent stiffness of various laminates with transverse edges free.

5.2 Case I results

The results for case I boundary conditions is shown in this section. For linear and cubic angle variation, $T_i\left(i=\mathrm{1,}\mathrm{2,}\mathrm{3,}4\right)$ was allows to vary from $0°$ to $90°$, with increments of $10°$. It was evaluated $N_{cr}$ and $E_x^{eq}$ for aspect ratios equals to 0.5, 1, 2 and 3. Besides that, it was compared the results for $\varphi =0°$ and $\varphi =90°$. The results of the highest values of $N_{cr}$ for each aspect ratio and the corresponding $E_x^{eq}$ are shown in Figs. 12 and 13. The layup of the highest value of $N_{cr}$ for each aspect ratio is shown in Tab. 7. The mode shapes of the best layups and higher critical buckling loads are shown in Fig. 14. As can be seen, they shown basically the same shape.

Table 7:
Best layups for each aspect ratio evaluated for Case I.

	$\varphi =0°$		$\varphi =90°$
$a/b$	Linear	Cubic	Linear	Cubic
0.5	${\left[\pm \langle 0°|20°\rangle \right]}_{3S}$	${\left[\langle 0°|0°|0°|0°\rangle \right]}_{3S}$	${\left[\langle 90°|90°\rangle \right]}_{3S}$	${\left[\langle 90°|90°|90°|90°\rangle \right]}_{3S}$
1	${\left[\pm \langle 0°|40°\rangle \right]}_{3S}$	${\left[\pm \langle 10°|20°|30°|40°\rangle \right]}_{3S}$	${\left[\pm \langle 20°|70°\rangle \right]}_{3S}$	${\left[\pm \langle 30°|40°|50°|80°\rangle \right]}_{3S}$
2	${\left[\pm \langle 0°|50°\rangle \right]}_{3S}$	${\left[\pm \langle 0°|10°|40°|40°\rangle \right]}_{3S}$	${\left[\pm \langle 30°|80°\rangle \right]}_{3S}$	${\left[\pm \langle 30°|40°|60°|90°\rangle \right]}_{3S}$
3	${\left[\pm \langle 0°|40°\rangle \right]}_{3S}$	${\left[\pm \langle 0°|10°|30°|40°\rangle \right]}_{3S}$	${\left[\pm \langle 30°|80°\rangle \right]}_{3S}$	${\left[\pm \langle 30°|40°|50°|80°\rangle \right]}_{3S}$

Figure 12:
Buckling performance comparison between cubic and linear angle variation for Case I.

Figure 13:
Equivalent stiffness performance comparison between cubic and linear angle variation for Case I.

Figure 14:
Buckling mode shape of the best layups for Case I, $a/b=3$ , and $\varphi =90°$.

For $\varphi =0°$ the values of $N_{cr}$ are practically the same for all aspect ratios for both linear and cubic angle variations. The largest difference occurs for $a/b=3$, it being the $N_{cr}$ for cubic variation 2.4% higher. However, the ratio $E_x^{eq}/E_1$ for cubic variation is 5.63%, indicating that the panel with cubic angle variation is more compliant than the linear variation one. Solely for $a/b=0.5$ the ratio $E_x^{eq}/E_1$ for the cubic variation is greater than the linear one, 2.13%. This can be explained due to the fact that the cubic variation layup for $a/b=0.5$ represents a straight fiber in the $x$ direction, while the layup for the linear varies linearly from $20°$ at its edges to $0°$ at its center.

To understand the difference observed for the $N_{cr}$ values between the linear and cubic variation for $a/b=3$, it is compared the behavior of $N_y$, shown in Fig. 15. As can be seen, $N_y$ for cubic variation is flatter at panel center compared to the $N_y$ for linear variation. This impact directly on matrix $M$ of the Eq. 19, yielding a higher value of ${\lambda }_{cr}$. Besides that, for $a/b=3$, $N_x$ of the cubic variation is 5.94% smaller than that of linear variation.

For $\varphi =90°$, the values of $N_{cr}$ for $a/b=0.5$ is identical, as they should be, considering that the linear and cubic variation have the same layup. For $a/b=1$, the cubic layup presents a value of $N_{cr}$ 1.45% greater, followed by a difference of 13% and 18.33% for $a/b=2$ and $3$, respectively. In the same way for $\varphi =0°$, the greater difference occurs when $a/b=3$. However, for $\varphi =90°$, this difference is significantly more pronounced. In Fig. 16 is compared the $N_x$ for linear and cubic variation. It is important to notice that when $\varphi =90°$, $\theta$ is function of y, i.e., $\theta \left(y\right)$. In a similar way to $\varphi =0°$, for a cubic variation $N_x$ is flatter at panel's center, and $N_y$ is 12.7% smaller than the linear variation.

Regarding the ratio $E_x^{eq}/E_1$, for $a/b=0.5$ the value is the same for both linear and cubic variation, as expected. For $a/b=1$, the cubic variation presented a value that is 8.57% greater than that for a linear variation. For $a/b=2$ and $3$, $E_x^{eq}/E_1$ for cubic variation is 12.25% smaller. As the panels for linear and cubic variation present similar angles at its edges and center, the difference in $E_x^{eq}/E_1$ is not so much expressive. Figure 17 shows the values of $N_{cr}$ and the corresponding ratio $E_x^{eq}/E_1$ for $T_0$ ranging from $0°$ to $90°$ with increments of $10°$ for $a/b=2$ and $\varphi =0°$, in the same manner that was made in the work of Gurdal et al. (2008). For $a/b=3$ the curves are similar and for $a/b=0.5$ and $a/b=1$, the curves are not continuous due to ${\kappa }_{\max }>{\kappa }_{cr}$. For this reason, they were omitted. As for $\varphi =90°$ the fiber angle is function of $y$, $\theta \left(y\right)$ , and the dimension $b$ keeps constant, the curves also are not continuous and were omitted.

Figure 15:
${\overline{N}}_y$ of the best laminates for $a/b=3$ and $\varphi =0°$.

Figure 16:
${\overline{N}}_x$ of the best laminates for $a/b=3$ and $\varphi =0°$.

Figure 17:
Buckling and equivalent stiffness of various laminates with $a/b=2$ and $\varphi =0°$ for Case I.

5.3 Case II results

In the same manner as Case I, $T_i\left(i=\mathrm{1,}\mathrm{2,}\mathrm{3,}4\right)$ was allows to vary from $0°$ to $90°$ , with increments of $10°$. It was evaluated $N_{cr}$ and $E_x^{eq}$ for aspect ratios equals to 0.5, 1, 2 and 3. It was compared the results for $\varphi =0°$ and $\varphi =90°$. Figure 18 shows the values of $N_{cr}$ for $\varphi =0°$ and $\varphi =90°$. It can be seen that for $\varphi =0°$ the values of $N_{cr}$ are basically the same for both linear and cubic variations. These results can be addressed through the panels' layups shown in Tab. 8. The mode shapes of the best layups and higher critical buckling loads are shown in Fig. 19. Unlike Case I, they shown different shapes. For the cubic case, there are 3 sine waves, while the linear case presents only one.

Table 8:
Best layups for each aspect ratio evaluated for Case II.

	$\varphi =0°$		$\varphi =90°$
$a/b$	Linear	Cubic	Linear	Cubic
0.5	${\left[\pm \langle 0°|20°\rangle \right]}_{3S}$	${\left[\langle 0°|0°|0°|0°\rangle \right]}_{3S}$	${\left[\langle 90°|90°\rangle \right]}_{3S}$	${\left[\pm \langle 50°|50°|60°|80°\rangle \right]}_{3S}$
1	${\left[\pm \langle 40°|50°\rangle \right]}_{3S}$	${\left[\pm \langle 50°|50°|50°|40°\rangle \right]}_{3S}$	${\left[\pm \langle 20°|70°\rangle \right]}_{3S}$	${\left[\pm \langle 20°|30°|50°|70°\rangle \right]}_{3S}$
2	${\left[\pm \langle 50°|40°\rangle \right]}_{3S}$	${\left[\pm \langle 50°|40°|40°|50°\rangle \right]}_{3S}$	${\left[\pm \langle 20°|70°\rangle \right]}_{3S}$	${\left[\pm \langle 20°|30°|50°|70°\rangle \right]}_{3S}$
3	${\left[\pm \langle 50°|40°\rangle \right]}_{3S}$	${\left[\pm \langle 40°|40°|40°|50°\rangle \right]}_{3S}$	${\left[\pm \langle 20°|70°\rangle \right]}_{3S}$	${\left[\pm \langle 30°|40°|50°|80°\rangle \right]}_{3S}$

Figure 18:
Buckling performance comparison between cubic and linear angle variation for Case II.

Figure 19:
Buckling mode shape of the best layups for Case II, $a/b=3$ , and $\varphi =90°$.

It is noticeable that the difference between the angles is not greater than $10°$, except for the layup ${\left[0°\pm \langle 0°|20°\rangle \right]}_{3S}$ for $a/b=0.5$. Thus, $N_x$ of this laminate presents a small variation when compared to straight fiber laminate ${\left[\langle 0°|0°|0°|0°\rangle \right]}_{3S}$. The difference between the value of $N_{cr}$ is 2.3% greater for linear variation. For $a/b$ equal to 1, 2 and 3, the layups are similar, changing only the angles at the edges and center of the panel from $40°$ to $50°$. Even for the cubic angle variation, the angles between the center and edge, $T_1$ and $T_2$, are equal. Figure 20 shows ${\overline{N}}_x$ for the laminates ${\left[0°\pm \langle 50°|40°\rangle \right]}_{3S}$ and ${\left[0°\pm \langle 40°|40°|40°|50°\rangle \right]}_{3S}$. It can be noticed that ${\overline{N}}_x$ presents, basically, a constant value.

Regarding $E_x^{eq}/E_1$, the results are shown in Fig. 20. When $a/b=0.5$, the laminate for the cubic variation is a straight fiber laminate in the $x$ direction, thus $E_x^{eq}/E_1$ is equal to 1, as expected. For the linear variation case, $E_x^{eq}/E_1$ is 12% smaller than the cubic one. As for $a/b$ equal to 1, 2 and 3 the laminates are similar, so the ratio $E_x^{eq}/E_1$ presents similar values as well. The difference between the linear and cubic variation is 18.18%, 13.33% and 18.75% for $a/b$ equal to 1, 2 and 3, respectively. Figure 22 shows the values of $N_{cr}$ and the corresponding ratio $E_x^{eq}/E_1$ for the linear variation and $a/b=2$. The same explanation about others aspect ratios for Case I is valid here.

Figure 20:
${\overline{N}}_x$ of the best laminates for $a/b=3$ and $\varphi =0°$.

Figure 21:
Equivalent stiffness performance comparison between cubic and linear angle variation for Case I.

For $\varphi =90°$, there are differences between linear and cubic angle variation. When $a/b=0.5$, for the linear variation, the laminate is in the $x$ direction. Therefore, $E_x^{eq}/E_1$ is equal to 1. The laminate for the cubic variation is ${\left[90°\pm \langle 50°|50°|60°|80°\rangle \right]}_{3S}$ and $E_x^{eq}/E_1$ is equal to 0.34. Regarding $N_{cr}$, for cubic variation $N_{cr}$ is 3.45% higher than linear variation. When $a/b=1$, the angles at edges and center are the same for linear and cubic variation. $N_{cr}$ for cubic variation is 6.83% higher than for linear variation, while $E_x^{eq}/E_1$ is 10% smaller.

For $a/b=2$, $N_{cr}$ is 2.97% higher for cubic variation, as for $E_x^{eq}/E_1$ is 11.11% smaller. The angles are the same at center and at edges, as for $a/b=1$. For $a/b=3$, $N_{cr}$ is 5.36% higher for cubic variation, as for $E_x^{eq}/E_1$ is 10% smaller. Figure 23 shows ${\overline{N}}_x$ of the laminates ${\left[90°\pm \langle 20°|70°\rangle \right]}_{3S}$ and ${\left[90°\pm \langle 30°|40°|50°|80°\rangle \right]}_{3S}$ for $a/b=3$. It can be noticed that for the cubic variation ${\overline{N}}_x$ is flatter at the center of the panel as well, in the same manner that for Case I.

Figure 22:
Buckling and equivalent stiffness of various laminates with $a/b=2$ and $\varphi =90°$ for Case II.

Figure 23:
${\overline{N}}_x$ of the best laminates for $a/b=3$ and $\varphi =90°$.

6 CONCLUSIONS

This work presented a numerical investigation that was conducted on how a cubic angle variation affects the in-plane and buckling response when compared to a linear angle variation of panels presenting various aspect ratios. It was imposed a constraint based on manufacturing considerations to avoid fiber defects such as wrinkles and micro-buckling. It was examined panels simply supported under applied uniform end shortening. Two boundary conditions were examined: the transverse edge fixed (Case I) and free to deform (Case II). Besides that, it was considered the fiber angle function varying along $x$ and $y$ direction, i.e., $\theta \left(x\right)$ and $\theta \left(y\right)$.

The results that have been obtained indicate that a cubic angle variation allows more flexibility to tailor the resultant stress caused by the in-plane loads, although it is more prone to achieve high fiber curvature values that can be an important design constraint. Regarding the critical buckling load, the cubic angle variation for the panels with fixed transverse edges exhibit higher values (Case I) when compared with the same panel presenting a linear angle variation. More significant difference between the cubic and linear variation buckling load happens when $a/b$ is equal to 3 and the fiber angles varies along $y$ direction, $\theta \left(y\right)$. For $\varphi =0°$, the gains in the buckling load using a cubic variation are marginal, for both Case I and II, it being less expressive for Case II.

7 ACKNOWLEDGMENTS

To Fundação de Amparo a Pesquisa do Estado de Minas Gerais (FAPEMIG) for financing part of this study.

References

• Brooks, T. R. and Martins, J. R. (2018). On manufacturing constraints for tow-steered composite design optimization. Composite Structures, 204:548 - 559.
• Clyne, T. W. and Hull, D. (2019). An introduction to composite materials. Cambridge university press, Cambridge.
• Gurdal, Z. and Olmedo, R. (1993). In-plane response of laminates with spatially varying fiber orientations-variable stiffness concept. AIAA journal, 31(4):751-758.
• Gurdal, Z., Tatting, B., and Wu, C. (2008). Variable stiffness composite panels: Effects of stiffness variation on the in-plane and buckling response. Composites Part A: Applied Science and Manufacturing, 39(5):911 - 922.
• Hao, P., Liu, C., Yuan, X., Wang, B., Li, G., Zhu, T., and Niu, F. (2017a). Buckling optimization of variable-stiffness composite panels based on flow field function. Composite Structures, 181:240-255.
• Hao, P., Yuan, X., Liu, H., Wang, B., Liu, C., Yang, D., and Zhan, S. (2017b). Isogeometric buckling analysis of composite variable-stiffness panels. Composite Structures, 165:192-208.
• Huang, J. and Haftka, R. (2005). Optimization of fiber orientations near a hole for increased load-carrying capacity of composite laminates. Structural and Multidisciplinary Optimization, 30(5):335-341.
• Hughes, T. J. (2000). The finite element method: linear static and dynamic finite element analysis. Dover Publications, New York.
• Hyer, M. W. and Charette, R. (1991). Use of curvilinear fiber format in composite structure design. AIAA journal, 29(6):1011-1015.
• IJsselmuiden, S. T., Abdalla, M. M., and Gurdal, Z. (2010). Optimization of variable-stiffness panels for maximum buckling load using lamination parameters. AIAA Journal, 48(1):134-143.
• Jones, R. M. (1998). Mechanics of composite materials. Taylor and Francis Group, New York.
• Kollar, L. P. and Springer, G. S. (2003). Mechanics of composite structures. Cambridge university press, Cambridge.
• Kumar, Y. (2018). The rayleigh-ritz method for linear dynamic, static and buckling behavior of beams, shells and plates: A literature review. Journal of Vibration and Control, 24(7):1205-1227.
• Li, R., Kelly, D., and Crosky, A. (2002). Strength improvement by fibre steering around a pin loaded hole. Composite Structures, 57(1-4):377-383.
• Lopes, C., Gürdal, Z., and Camanho, P. (2008). Variable-stiffness composite panels: Buckling and first-ply failure improvements over straight-fibre laminates. Computers & Structures, 86(9):897-907.
• Marouene, A., Boukhili, R., Chen, J., and Yousefpour, A. (2016). Buckling behavior of variable-stiffness composite laminates manufactured by the tow-drop method. Composite Structures, 139:243 - 253.
• Marsh, G. (2010). Airbus A350 XWB update. Reinforced Plastics, 54(6):20-24.
• Mikhailov, M. (2008). Modelling unsaturated flow by using mathematica. Communications in numerical methods in engineering, 24(11):947-959.
• Nik, M. A., Fayazbakhsh, K., Pasini, D., and Lessard, L. (2014). Optimization of variable stiffness composites with embedded defects induced by automated fiber placement. Composite Structures, 107:160-166.
• Oñate, E. (2013). Structural analysis with the finite element method. Linear statics: volume 2: beams, plates and shells. Springer Science & Business Media, Barcelona.
• Sadd, M. (2009). Elasticity: Theory, Applications, and Numerics. Elsevier Science, Boston.
• Setoodeh, S., Abdalla, M. M., IJsselmuiden, S. T., and Gurdal, Z. (2009a). Design of variable-stiffness composite panels for maximum buckling load. Composite Structures, 87(1):109 - 117.
• Setoodeh, S., Abdalla, M. M., IJsselmuiden, S. T., and Gürdal, Z. (2009b). Design of variable-stiffness composite panels for maximum buckling load. Composite structures, 87(1):109-117.
• Tatting, B. F. and Gurdal, Z. (2002). Design and manufacture of elastically tailored tow placed plates (NASA/CR-2002-211919). Technical report, NASA Langley Research Center, Hampton, Virginia.
• Van den Brink, W., Vankan, W., and Maas, R. (2012). Buckling-optimized variable stiffness laminates for a composite fuselage window section. National Aerospace Laboratory NLR.
• Whitney, J. M. (1987). Structural analysis of laminated anisotropic plates. Technomic Publishing Company, Inc., Lancaster.
• Wolfram, S. et al. (1999). The MATHEMATICA® book. Cambridge University Press, Cambridge.
• Wu, Z., Weaver, P. M., Raju, G., and Chul Kim, B. (2012). Buckling analysis and optimisation of variable angle tow composite plates. Thin-Walled Structures, 60:163 - 172.
• Wu, Z., Weaver, P. M., and Raju, G. (2013). Postbuckling optimisation of variable angle tow composite plates. Composite Structures, 103:34-42.

APPENDIX

$I_1=\underset{0}{\overset{a}{{\displaystyle \int }}}\underset{0}{\overset{b}{{\displaystyle \int }}}\left[\frac{i^2{k}^2}{a^4}{D}_{11}+\left(\frac{i^2{k}^2+j^2{k}^2}{a^2{b}^2}\right)D_{12}+\frac{j^2{l}^2}{b^4}{D}_{66}\right]\sin i\overline{x}\sin j\overline{y}\sin k\overline{x}\sin l\overline{y}dxdy$

$I_2=\underset{0}{\overset{a}{{\displaystyle \int }}}\underset{0}{\overset{b}{{\displaystyle \int }}}{D}_{66}\cos i\overline{x}\cos j\overline{y}\cos k\overline{x}\cos l\overline{y}dxdy$

$I_3=\underset{0}{\overset{a}{{\displaystyle \int }}}\underset{0}{\overset{b}{{\displaystyle \int }}}\left(\frac{i{k}^{2}j}{a^{3}b}{D}_{16}+\frac{ij{l}^2}{a{b}^3}{D}_{26}\right)\cos i\overline{x}\cos j\overline{y}\sin k\overline{x}\sin l\overline{y}dxdy$

$I_4=\underset{0}{\overset{a}{{\displaystyle \int }}}\underset{0}{\overset{b}{{\displaystyle \int }}}\left(\frac{i^{2}kl}{a^{3}b}{D}_{16}+\frac{j^{2}kl}{a{b}^3}{D}_{26}\right)\sin i\overline{x}\sin j\overline{y}\cos k\overline{x}\cos l\overline{y}dxdy$

$I_5=\underset{0}{\overset{a}{{\displaystyle \int }}}\underset{0}{\overset{b}{{\displaystyle \int }}}{N}_x\cos i\overline{x}\cos j\overline{y}\cos k\overline{x}\cos l\overline{y}dxdy$

$I_6=\underset{0}{\overset{a}{{\displaystyle \int }}}\underset{0}{\overset{b}{{\displaystyle \int }}}{N}_y\sin i\overline{x}\cos j\overline{y}\sin k\overline{x}\cos l\overline{y}dxdy$

$I_7=\underset{0}{\overset{a}{{\displaystyle \int }}}\underset{0}{\overset{b}{{\displaystyle \int }}}{N}_{xy}\sin i\overline{x}\cos j\overline{y}\sin k\overline{x}\cos l\overline{y}dxdy$

$I_8=\underset{0}{\overset{a}{{\displaystyle \int }}}\underset{0}{\overset{b}{{\displaystyle \int }}}{N}_{xy}\text{cos }i\overline{x}\text{sin }j\overline{y}\text{cos }k\overline{x}\text{sin }l\overline{y}dxdy$

Edited by

Editor:

Marco L. Bittencourt

Publication Dates

• Publication in this collection
  10 Nov 2021
• Date of issue
  2021

History

• Received
  24 Feb 2021
• Reviewed
  29 June 2021
• Accepted
  29 June 2021
• Published
  05 July 2021