# Buckling Analysis of Variable-Angle Tow Composite Plates through Variable Kinematics Hierarchical Models

Gaetano Giunta, Domenico Andrea Iannotta, Levent Kirkayak, Marco Montemurro

► To cite this version:

Gaetano Giunta, Domenico Andrea Iannotta, Levent Kirkayak, Marco Montemurro. Buckling Analysis of Variable-Angle Tow Composite Plates through Variable Kinematics Hierarchical Models. Journal of Composites Science, 2024, 8 (8), pp.320. 10.3390/jcs8080320. hal-04719624

HAL Id: hal-04719624

https://hal.science/hal-04719624v1

Submitted on 3 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Buckling Analysis of Variable-Angle Tow Composite Plates through Variable Kinematics Hierarchical Models

Gaetano Giunta, Domenico Andrea Iannotta, Levent Kirkayak and Marco Montemurro

Special Issue

Characterization and Modelling of Composites, Volume III

Edited by
Dr. Stelios K. Georgantzinos

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

Article

# Buckling Analysis of Variable-Angle Tow Composite Plates through Variable Kinematics Hierarchical Models

Gaetano Giunta <fnref n="1" />,* , Domenico Andrea Iannotta <fnref n="1" />,<fnref n="2" />,<fnref n="3" /> , Levent Kirkayak <fnref n="1" /> and Marco Montemurro <fnref n="3" />

**Abstract:** Variable-Angle Tow (VAT) laminates can improve straight fiber composites' mechanical properties thanks to the application of curvilinear fibers. This characteristic allows one to achieve ambitious objectives for design and performance purposes. Nevertheless, the wider design space and the higher number of parameters result in a more complex structural problem. Among the various approaches that have been used for VAT study, Carrera's Unified Formulation (CUF) allows one to obtain multiple theories within the same framework, guaranteeing a good compromise between the results' accuracy and the computational cost. In this article, the linear buckling behavior of VAT laminates is analyzed through the extension of CUF 2D plate models within Reissner's Mixed Variational Theorem (RMVT). The results show that RMVT can better approximate the prebuckling nonuniform stress field of the plate when compared to standard approaches, thus improving the prediction of the linear buckling loads of VAT composites.

**Keywords:** buckling analysis; finite element method; variable-angle tow plates; Carrera's unified formulation; Reissner's mixed variational theorem

## 1. Introduction

Composite materials have attracted considerable interest in various fields of application over the last decades. Thanks to their high stiffness-to-weight ratio, these materials gained great importance in those contexts where light structures with good mechanical properties are needed: this is the case in aerospace, automotive, and construction industries. Despite this, a prevailing belief is that the optimal use of laminate structures can still be improved in order to fully unlock the potential of fiber properties. For instance, when utilizing composites with straight fibers and a constant thickness, it becomes challenging to fully profit from the directional properties of this kind of material. This limitation becomes especially problematic for complex geometries featuring geometrical discontinuities such as cut-outs. VAT plates are distinguished by the variation in fiber angles within the structure plane, which significantly enhances the design possibilities for a specific structure. Originally, VATs were produced using automated tape placement and automated fiber- placement techniques [1]. Modern production techniques, such as Additive Manufacturing (AM) technologies, enable the overcoming of limitations related to automated processes, such as the existence of defects like overlaps. AM entails the incremental deposition of materials layer by layer to fabricate a three-dimensional object. In the context of VAT com- posites, AM processes are employed to deposit and cure layers of composite materials with varying fiber orientations, as explained by Zhuo et al. [2]. A key challenge in employing VATs lies in the complexity of the analysis. In such scenarios, it becomes essential to account

<figure><img src="image_5.png" /></figure>

### check for updates

**Citation:** Giunta, G.; Iannotta, D.A.; Kirkayak, L.; Montemurro, M. Buckling Analysis of Variable-Angle Tow Composite Plates through Variable Kinematics Hierarchical Models. *J. Compos. Sci.* **2024**, *8*, 320. https://doi.org/10.3390/jcs8080320

Academic Editor: Stelios K. Georgantzinos

Received: 14 May 2024

Revised: 6 August 2024

Accepted: 9 August 2024

Published: 13 August 2024

<figure><img src="image_6.png" /></figure>

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Luxembourg Institute of Science and Technology, 5 Avenue des Hauts-Fourneaux, L-4362 Esch-sur-Alzette, Luxembourg; domenico-andrea.iannotta@list.lu (D.A.I.); levent.kirkayak@list.lu (L.K.)</footnote>

<footnote n="2">² Doctoral School in Science and Engineering (DSSE), University of Luxembourg, 2 Avenue de l'Université, L-4365 Esch-sur-Alzette, Luxembourg</footnote>

<footnote n="3">³ Arts et Métiers Institute of Technology, Université de Bordeaux, CNRS, INRA, Bordeaux INP, HESAM Université, I2M UMR 5295, F-33405 Talence, France; marco.montemurro@ensam.eu</footnote>

<footnote>* Correspondence: gaetano.giunta@list.lu</footnote>
for a higher number of variables, and during the optimization process, it is possible to end
up with unfeasible fiber patterns.

The following is a concise summary of typical techniques employed in the analysis of the mechanical behavior of VATs, with particular emphasis on buckling investigations. One of the early methods made the assumption of a constant fiber angle within each element in a Finite Element (FE) framework. This approach approximates the continuous variation in fiber direction discretely in a stepwise manner. Hyer and Lee in [3] used this technique to maximize the buckling load of a composite plate through a sensitivity analysis. The fiber angle corresponding to the maximum buckling load was identified for each element. Through a gradient-search technique, the influence of orientation changes on the buckling load was assessed. The buckling load was then optimized through recursive modifications of the elements' orientations, also considering the interactions between elements. However, a key limitation of these approaches arises from their reliance on an element-wise representation of the fiber paths. With this approach, the optimized solutions might exhibit fragmented fiber paths, making them impractical for manufacturing and necessitating a computationally expensive postprocessing phase for adjustment [4]. A discrepancy between the structural responses computed based on the optimized fiber path and those computed after the optimized solution underwent postprocessing (recovery) has often been observed.

Gürdal et al. [5] utilized the classical lamination theory to analyze the mechanical response of VAT plates. They used the Rayleigh–Ritz method to conduct buckling analyses. The same method was applied by Oliveri and Milazzo [6] for analyzing the postbuckling behavior of VAT-stiffened panels through the first-order shear deformation theory. The structures were modeled through the domain discretization into plate-like subdomains by imposing boundary conditions for each component. Hao et al. [7] used the Mindlin plate theory within an isogeometric analysis to perform a linear buckling study of variable-stiffness panels, ensuring the continuity of the fiber angle on the whole structure plane. Hao et al. [8] introduced a flow field function for the representation of fiber paths through a reduced number of variables. The buckling optimization of variable-stiffness composite panels with multiple cut-outs was performed through a bi-level optimization framework. In the first level, straight fiber paths with optimized orientation angles were first established to achieve a quasi-optimal design for maximizing buckling load and, in the second level, the flow field function was utilized to represent the fiber path. Raju et al. [9] used Airy’s stress function to accurately predict the prebuckling behavior of VAT plates. They used the differential quadrature method and the classical laminated plate theory to analyze the buckling response of the VATs modeled. Sciascia et al. [10] studied VAT shells considering buckling, free vibrations, and prestressed vibrations. The eigenvalue problem was solved through the Ritz formulation, while the shell kinematics was described with the first-order shear deformation theory. Chen et al. [11] presented an analytical model to predict the global, mixed, and local buckling response of VAT plates with delaminations. The Rayleigh–Ritz approach was used to solve the prebuckling and buckling problems. Wu et al. [12] introduced a two-level optimization framework for the buckling of VAT plates. VATs were defined through lamination parameters represented by B-spline entities.

The Multiscale Two-Level (MS2L) approach leverages a dual-scale analysis of compos-
ite structures. It considers both micro and macro levels, enabling a detailed examination
of material properties at the microscopic level while accounting for the overall structural
behavior. This method allows one to split the optimization problem into two levels. In
the first step, the composite is considered as an equivalent homogeneous anisotropic plate,
aiming to determine the optimal distribution of geometric and mechanical design variables
governing VAT structure behavior at this scale. Among various macroscopic anisotropy
representations, the most efficient one utilizes polar parameters (see Montemurro [13,14]).
In the second step, the focus shifts to the mesoscopic scale, and the objective is to identify at
least one stack that aligns with the optimized arrangement of the polar parameters obtained
from the initial optimization phase. The MS2L technique was applied by Montemurro
and Catapano [15] with the aim of maximizing the buckling load of VAT plates. In the initial step, B-spline surfaces are used to represent the distribution of the polar parameters across the structure, while in the subsequent step, manufacturing constraints are taken into consideration. The initial framework for the gradient-based optimization of VAT structures was introduced by Montemurro and Catapano in [16,17]. Its application to eigenvalue buckling problems was conducted by Fiordilino et al. [18], who introduced the analytical expression of the gradient of objective and constraint functions in buckling problems.

Carrera's Unified Formulation (CUF) has shown accurate results in the analysis of VATs. Considering plate geometries, CUF allows one to choose a priori the approximation functions along the thickness (Carrera [19,20]). Early works concerning CUF presented the development of a Navier closed-form solution, which has been applied by Carrera et al. [21] to perform the static and failure analysis of isotropic plates under various loading conditions. Giunta et al. [22] used CUF in order to study the linear buckling of thin-walled beams. Hui et al. [23] conducted a multiscale nonlinear analysis of composite beam structures using a series of one-dimensional (1D) CUF models. The asymptotic numerical method was employed to study the impact of microscopic imperfections on the macroscale response. Specifically, the influence of fiber defects was investigated by introducing sinusoidal geometries. By considering these imperfections, this work aimed to understand how they affect the overall behavior of the composite beam structures at a larger scale. A data-driven computational mechanics approach has also been employed with CUF models. In this context, Hui et al. [24] used 1D CUF models for a static analysis of beam structures and found that the accuracy is dependent on the number of layers in the database. This demonstrates the potential for CUF models to be integrated with modern computing techniques. Viglietti et al. [25] and Fallahi et al. [26] employed a 1D CUF model to conduct free-vibration and buckling analyses of VAT structures, respectively. Fallahi et al. [27] further investigated buckling optimization in VATs using a genetic algorithm, applied to a one-dimensional CUF model. Sánchez-Majano et al. [28] studied the effects of manufacturing defects, such as varying fiber volume and fiber misalignments, on the buckling response of VAT plates. The Monte Carlo method was applied to perform stochastic buckling analyses. Pagani and Sánchez-Majano [29] incorporated Monte Carlo simulations and layer-wise CUF models in order to explore the impact of mesoscale fiber misalignments as defects in VAT laminates.

The extension of CUF allows for the use of different expansions for each component of the displacement vector. This approach was demonstrated by Demasi et al. [30] for the study of VAT plates using triangular elements. Vescovini and Dozio [31] used Ritz's method within CUF for vibrational and buckling analyses of VATs. CUF offers an additional advantage through its compatibility with various variational formulations. Reissner's Mixed Variational Theorem, an alternative to the traditional Principle of Virtual Displacements, considers both displacements and out-of-plane transverse stresses as unknowns, demonstrating the flexibility of CUF in its applications. RMVT was widely used within CUF for the study of straight fiber composite structures. Carrera and Demasi [32] showed how to apply RMVT to CUF from a theoretical point of view, whereas in the second part of the article [33], static analyses were performed on straight fiber plates.

The first work regarding the application of RMVT within CUF for the analysis of VAT plates was published by Giunta et al. [34], showing an improvement in the accuracy of results when the mixed formulation is used for the dynamic analysis of this kind of structure. Subsequently, the same approach was used by Iannotta et al. [35] to enhance the prediction of VAT composites' static response. The goal of this paper is to extend the previous works to the eigenvalue buckling analysis when curvilinear fibers are present. The structure of the paper is as follows: Section 2 introduces the CUF approach and demonstrates its application to the two variational statements mentioned earlier. Section 3 considers three benchmark problems and performs linear buckling analyses with various loading conditions to discuss the differences between PVD and RMVT models and the reference solutions obtained in Abaqus using three-dimensional (3D) elements. Finally, Section 4 draws the conclusive remarks.
## 2. Carrera's Unified Formulation

A plate is a flat and thin structural element having in-plane dimensions significantly greater than its thickness. A global Cartesian reference system is considered to describe displacements, strains, and stresses. The x-axis and y-axis are aligned with the in-plane sides of the plate, whose dimensions are indicated as a and b, respectively. The z-axis is oriented perpendicular to the plane of the plate, in the same direction as the thickness h, which is negligible in comparison with a and b. The problem geometry and global reference system are illustrated in Figure 1.

<figure><img src="image_7.png" /><figcaption>Figure 1. Plate geometry and reference system.</figcaption></figure>

The displacement field is expressed as:

$$ \mathbf{u} = \begin{cases} u_x \\ u_y \\ u_z \end{cases}. \qquad (1) $$

The strain vector can be represented through Voigt's notation and then divided into its in-plane and out-of-plane components:

$$ \epsilon_p = \begin{cases} \epsilon_{xx} \\ \epsilon_{yy} \\ \epsilon_{xy} \end{cases}, \quad \epsilon_n = \begin{cases} \gamma_{xz} \\ \gamma_{yz} \\ \gamma_{zz} \end{cases}. \qquad (2) $$

The assumption of small displacements permits the utilization of a linear relationship:

$$ \epsilon_p = D_p \mathbf{u}, \\ \epsilon_n = (D_{n\Omega} + D_{nz}) \mathbf{u}, \qquad (3) $$

where $D_p, D_{n\Omega}$, and $D_{nz}$ are the following differential operators:

$$ D_p = \begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0 \\ 0 & \frac{\partial}{\partial y} & 0 \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \end{bmatrix}, \quad D_{n\Omega} = \begin{bmatrix} 0 & 0 & \frac{\partial}{\partial x} \\ 0 & 0 & \frac{\partial}{\partial y} \\ 0 & 0 & 0 \end{bmatrix}, \quad D_{nz} = \begin{bmatrix} \frac{\partial}{\partial z} & 0 & 0 \\ 0 & \frac{\partial}{\partial z} & 0 \\ 0 & 0 & \frac{\partial}{\partial z} \end{bmatrix}. \qquad (4) $$

The nonlinear part of the strains can be introduced in a Green-Lagrange sense (see [36]):

$$ \delta\epsilon_{ij}^{nl} = \frac{1}{2} \left( \frac{\partial u_x}{\partial i} \frac{\partial u_x}{\partial j} + \frac{\partial u_y}{\partial i} \frac{\partial u_y}{\partial j} + \frac{\partial u_z}{\partial i} \frac{\partial u_z}{\partial j} \right) \text{ with } i, j = x, y, z. \qquad (5) $$

The complete Cartesian components of strain can be obtained by summing the linear and nonlinear terms. The in-plane and out-of-plane stress components can be written as

$$ \sigma_p = \begin{cases} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{cases}, \quad \sigma_n = \begin{cases} \sigma_{xz} \\ \sigma_{yz} \\ \sigma_{zz} \end{cases}. \qquad (6) $$
Hooke's law reads

$$ 
\begin{aligned}
\sigma_p &= \tilde{C}_{pp}\epsilon_p + \tilde{C}_{pn}\epsilon_n, \\
\sigma_n &= \tilde{C}_{np}\epsilon_p + \tilde{C}_{nn}\epsilon_n,
\end{aligned}
\qquad (7) $$

where the terms $\tilde{C}_{pp}$, $\tilde{C}_{pn}$, $\tilde{C}_{np}$, and $\tilde{C}_{nn}$ are the components of the material stiffness matrix.

## 2.1. Variable-Angle Tow Composite Plates

VAT laminates are characterized by a point-wise variation in the material stiffness matrix components along the in-plane directions. Because of the laminate stacking sequence, the material stiffness matrix is also subjected to a layer-wise variation along the thickness of the plate. The equation enabling the rotation of the material stiffness matrix C by a designated angle $\theta$ around the z-axis is expressed as follows:

$$ \tilde{C} = TCT^T. \qquad (8) $$

Here, C denotes the material stiffness matrix in the material reference system, whereas $\tilde{C}$ represents the matrix after a rotation. The rotation matrix T is a function of the angle $\theta$. For conciseness, the specific components of C and T are omitted in this context (refer to Reddy [36] for detailed information). A linear variation law can be formulated as

$$ \theta(\alpha) = \Phi + T_0 + \frac{T_1 - T_0}{d}|\alpha|. \qquad (9) $$

The angle $\Phi$ represents the initial direction of variation for $\theta$, and $\alpha$ is a spatial variable given as

$$ \alpha = x' \cos(\Phi) + y' \sin(\Phi). \qquad (10) $$

where $x'$ and $y'$ denote the axes of the angle reference system. $T_0$ indicates the initial fiber angle when $\alpha = 0$, and $T_1$ represents the fiber angle when $\alpha = d$, as shown in Figure 2.

<figure><img src="image_8.png" /><figcaption>Figure 2. Example of in-plane fiber path.</figcaption></figure>

Within this paper (and without a loss of generality), the fiber angle is measured with respect to the x'-axis. The direction of variation of $\theta$ can be along $x'$, $y'$, or a combination of them, depending on the specific case. The notation $\Phi < T_0, T_1 >$ is used to describe the in-plane path of the fibers, which is based on the parameters introduced earlier. More details about the variation law of the local fiber orientation can be found in Gürdal et al. [5].

## 2.2. Variational Formulation

The governing equations are here derived through the consideration of PVD and RMVT variational statements. The key distinction lies in PVD focusing solely on displacements as unknowns, whereas RMVT incorporates both displacements and transverse stresses, denoted as $\sigma_n$, as primary unknowns. To effectively address the buckling problem, a prebuckling analysis is essential for determining the nonuniform stress distribution resulting from an external load. For the PVD case, the virtual internal work can be expressed as

$$ \delta L_i^{\text{PVD}} = \int_{\Omega} \int_{h} \left( \delta \epsilon_{pG}^{T} \sigma_{pH} + \delta \epsilon_{nG}^{T} \sigma_{nH} \right) dz \, d\Omega. \qquad (11) $$

The subscript 'G' refers to the components derived from the geometrical relations outlined in Equation (3), whereas the subscript 'H' corresponds to components obtained through
Hooke's law as specified in Equation (7). The term $\Omega$ represents the mid surface of the plate in the in-plane directions, and the subscript 'T' stands for the transpose of a vector/matrix. The virtual internal work for the RMVT case is

$$ \delta L_i^{\text{RMVT}} = \int_{\Omega} \int_{h} \left[ \delta \epsilon_{pG}^T \sigma_{pH} + \delta \epsilon_{nG}^T \sigma_{nM} + \delta \sigma_{nM}^T (\epsilon_{nG} - \epsilon_{nH}) \right] dz \, d\Omega, \quad (12) $$

where the 'M' subscript refers to the transverse stress components considered as primary unknowns in the mixed formulation. In the RMVT formulation, Hooke's law is expressed as follows:

$$ \begin{aligned} \sigma_{pH} &= \hat{C}_{pp} \epsilon_{pG} + \hat{C}_{pn} \sigma_{nM}, \\ \epsilon_{nH} &= \hat{C}_{np} \epsilon_{pG} + \hat{C}_{nn} \sigma_{nM}, \end{aligned} \quad (13) $$

where $\hat{C}_{pp}$, $\hat{C}_{pn}$, $\hat{C}_{np}$, and $\hat{C}_{nn}$ are obtained by the following relations (see Carrera and Demasi [32]):

$$ \begin{aligned} \hat{C}_{pp} &= \tilde{C}_{pp} - \tilde{C}_{pn} \tilde{C}_{nn}^{-1} \tilde{C}_{np}, \\ \hat{C}_{pn} &= \tilde{C}_{pn} \tilde{C}_{nn}^{-1}, \\ \hat{C}_{np} &= -\tilde{C}_{nn}^{-1} \tilde{C}_{np}, \\ \hat{C}_{nn} &= \tilde{C}_{nn}^{-1}. \end{aligned} \quad (14) $$

### 2.2.1. Prebuckling Problem

In the context of a prebuckling analysis, the following balance equation applies:

$$ \delta L_i = \delta L_e, \quad (15) $$

where the virtual work associated with external loads reads

$$ \delta L_e = \int_{\Omega_p} \delta \mathbf{u}^T \mathbf{p} \, d\Omega, \quad (16) $$

where $\mathbf{p}$ is a surface load applied on the surface $\Omega_p$ at z coordinate $z_p$.

### 2.2.2. Linear Buckling Problem

The stress field obtained through the prebuckling analysis can be integrated for the computation of the prestresses virtual work $\delta L_{\sigma^0}$, which reads

$$ \delta L_{\sigma^0} = \int_{\Omega} \int_{h} \left( \delta \epsilon_{xx}^{nl} \sigma_{xx}^{0} + \delta \epsilon_{yy}^{nl} \sigma_{yy}^{0} + \delta \epsilon_{xy}^{nl} \sigma_{xy}^{0} + \delta \epsilon_{xz}^{nl} \sigma_{xz}^{0} + \delta \epsilon_{yz}^{nl} \sigma_{yz}^{0} + \delta \epsilon_{zz}^{nl} \sigma_{zz}^{0} \right) dz \, d\Omega, \quad (17) $$

where $\sigma_{ij}^0$ is the prestresses energetically conjugate to the nonlinear strains $\epsilon_{ij}^{nl}$, with $i, j = x, y, z$. The linear buckling problem can be written for both PVD and RMVT cases as

$$ \delta L_i = \delta L_{\sigma^0}. \quad (18) $$

### 2.3. Kinematic Assumption and Finite Element Approximation

CUF allows one to introduce an axiomatic approximation in order to mathematically represent the primary unknowns along the direction parallel to the thickness (see Carrera [20]). Considering $f = f(x, y, z)$ as a generic unknown component, the following expansion can be introduced:

$$ f(x, y, z) = F_{\tau}(z)g_{\tau}(x, y), \quad \tau = 0, 1, \dots, N. \quad (19) $$

In this context, $f$ can represent solely a displacement component within a formulation derived by the PVD or a displacement or an out-of-plane stress component in the case of an RMVT formulation. $F_{\tau}$ serves as an approximation function along the thickness, while $g_{\tau}$ represents an unknown two-dimensional function that approximates the in-plane variation. The used Einstein's notation assumes that a twice-repeated index implies a sum
over that index range. This notation simplifies mathematical expressions and calculations. N represents the approximation order, which can be chosen a priori. Similarly, $F_τ$ can also be imposed beforehand. The flexibility of the CUF enables the formulation of multiple theories within the same framework, making it a versatile tool for VAT analysis. Equivalent single layer (ESL) and layer-wise (LW) theories can be obtained through the adequate choice of $F_τ$ and a coherent implementation of the layers' stiffness matrices assembly.

### 2.3.1. Equivalent Single-Layer Theories

In the ESL framework, 1D polynomials of the type $z^τ$ are used as $F_τ$ in order to obtain Taylor's expansion for the representation of primary variables:

$$F_{\tau}(z) = z^{\tau}, \quad \tau = 0, 1, \dots, N. \qquad (20)$$

In this case, the number of layers does not affect the total number of unknowns. The components of the stiffness matrix are computed by summing the integrals of the thickness functions for each layer, each multiplied by its corresponding stiffness coefficient. Even if ESL models are characterized by a diminished computational complexity/cost and show a good prediction of thin laminate behavior, they lack accuracy when thick plates are considered. Due to their reliance on $C^\infty$ approximation functions, ESL approaches are unable to accurately capture the zigzag displacement effect. However, it might be feasible to incorporate this feature by introducing Murakami's function, as detailed in Carrera [37].

### 2.3.2. Layer-Wise Theories

Lagrange or Legendre polynomials can be used to approximate the unknown fields independently for each layer, resulting, in this way, in a LW model for which it is possible to write the following approximation along the z-axis:

$$f^k(x,y,z) = F_b(z)g_b^k(x,y) + F_r(z)g_r^k(x,y) + F_t(z)g_t^k(x,y), \quad r = 2, \dots, N. \qquad (21)$$

The superscript 'k' denotes a generic layer within the structure, with k ranging from one to $N_l$, where $N_l$ represents the total number of layers. Subscripts 't' and 'b' correspond to the top and bottom faces of the generic layer, respectively. In the case of Legendre polynomials, the approximating functions along the z-axis are

$$F_t(z(\xi_k)) = \frac{P_0 + P_1}{2}, \quad F_b(z(\xi_k)) = \frac{P_0 - P_1}{2}, \quad F_r(z(\xi_k)) = P_r - P_{r-2}, \quad r = 2, \dots, N, \qquad (22)$$

where $P_i = P_i(\xi_k)$ represents the ith-order Legendre polynomial defined within the domain of the kth layer with the dimensionless through-the-thickness local coordinate $\xi_k$ bounded within -1 and +1. LW models can accurately predict the zigzag through-the-thickness behavior of the displacement field. However, they come with a higher computational cost, as the number of unknowns depends on the number of layers.

### 2.3.3. Finite Element Formulation

When an FE approximation is applied, it is necessary to incorporate the shape functions into the formulation. In the case of a 2D model, Equation (19) becomes

$$f(x,y,z) = F_{\tau}(z)N_{i}(x,y)g_{\tau i}, \quad \tau = 0, 1, \dots, N, \quad i = 1, \dots, N_{n}, \qquad (23)$$

where $N_i$ denotes the shape functions of the element and $N_n$ corresponds to the number of nodes utilized for the domain discretization. Classical Lagrange shape functions are used. They are not explicitly detailed here for the sake of brevity, and interested readers can refer to [38] for further information.
## 2.4. Acronym System

A comprehensive acronym system is implemented to identify the derived theories. The first letter designates the employed approximation level: 'E' corresponds to ESL models, while 'L' corresponds to LW models. The second letter denotes the variational statement: 'D' or 'M' signifies PVD or RMVT, respectively. The final number indicates the expansion order utilized along the plate thickness. If the first number is present, it signifies the number of virtual layers employed for the LW model to represent each physical layer. If the acronym begins without a numerical value, it is implicitly understood that only one virtual layer has been used for each physical layer. This system is shown in Figure 3.

<figure><img src="image_9.png" /><figcaption>Figure 3. Acronym system.</figcaption></figure>

For instance, the displacement field in EDN models corresponds to the following vectorial form:

$$u = F_0u_0 + F_1u_1 + \dots + F_Nu_N = F_{\tau}u_{\tau}, \quad \tau = 0, 1, \dots, N, \qquad (24)$$

where $F_{\tau} = z^{\tau}$ and $u_{\tau} = u_{\tau}(x, y)$. Classical theories arise as specific cases of the ED1 solution: classical lamination theory (CLT) and first-order shear deformation theory (FSDT) are denoted as CLT and FSDT, respectively. FSDT is derived by penalizing the $u_{z1}$ term in a first-order through-the-thickness approximation, whereas CLT requires the penalization of the terms coming from the work of the transverse shear stresses. The material stiffness matrix is reduced to account for a plane stress state in both CLT and FSDT and avoid thickness locking.

For LDN solutions, only displacements are considered as unknowns:

$$u^k = F_0u_0^k + F_1u_1^k + \dots + F_Nu_N^k = F_{\tau}u_{\tau}^k, \quad \tau = 0, 1, \dots, N, \quad k = 1, 2, \dots, N_l. \qquad (25)$$

For LMN solutions, transverse stresses are also included among the unknowns:

$$\sigma_n^k = F_b\sigma_b^k + F_r\sigma_r^k + F_t\sigma_t^k, \quad r = 2, \dots, N, \quad k = 1, 2, \dots, N_l. \qquad (26)$$

In both cases, $N$ refers to the approximation order employed in each layer. Noticeably, ESL theories can be viewed as specific instances of LW theories. In ESL, integration along the thickness is carried out to represent composite properties through an equivalent single layer. In contrast, in LW theories, integration is computed layer by layer. This allows for the individual representation of the kinematics of each layer in LW models. Unless otherwise specified, LDN solutions are derived using Lagrange polynomials with equally spaced nodes, while LMN solutions utilize Legendre polynomials.
## 2.5. Stiffness Matrices Expression

Considering PVD, the displacement field groups the primary unknowns. In a PVD context, displacements from Equation (23) can be written as follows:

$$ \mathbf{u} = F_{\tau} N_{i} \begin{cases} q_{x\tau i} \\ q_{y\tau i} \\ q_{z\tau i} \end{cases} = F_{\tau} N_{i} \mathbf{q}_{\tau i}. \qquad (27) $$

The PVD governing equations can be written in a compact vectorial form:

$$ \delta \mathbf{q}_{\tau i}^{T} \mathbf{K}^{\tau sij} \mathbf{q}_{sj} = \delta L_{\sigma 0}, \qquad (28) $$

where $\mathbf{K}^{\tau sij}$ is a 3 × 3 Fundamental Nucleus (FN). The loops on the indices $\tau$, $s$, $i$, and $j$ allow one to build the stiffness matrix of the whole plate element expanding over the kinematic approximation and the Finite Element approximation. In the RMVT case, transverse stresses also constitute an unknown field:

$$ \sigma_n = F_{\tau} N_i \begin{cases} g_{x\tau i} \\ g_{y\tau i} \\ g_{z\tau i} \end{cases} = F_{\tau} N_i g_{\tau i}. \qquad (29) $$

The RMVT governing equations can be written in the following compact form:

$$ \begin{gathered} \delta \mathbf{q}_{\tau i}^{T} \mathbf{K}_{uu}^{\tau sij} \mathbf{q}_{sj} + \delta \mathbf{q}_{\tau i}^{T} \mathbf{K}_{u\sigma}^{\tau sij} \mathbf{g}_{sj} = \delta L_{\sigma 0}, \\ \delta \mathbf{g}_{\tau i}^{T} \mathbf{K}_{u\sigma}^{\tau sij} \mathbf{q}_{sj} + \delta \mathbf{g}_{\tau i}^{T} \mathbf{K}_{\sigma\sigma}^{\tau sij} \mathbf{g}_{sj} = 0. \end{gathered} \qquad (30) $$

In this case, four fundamental nuclei are obtained. Gauss quadrature was used to compute the in-plane integrals. Due to the varying stiffness coefficients of the materials, the number of Gauss points must be adjusted for a given analytical formula describing the local fiber orientation to ensure converged integration. A reduced integration is used to correct the shear-locking phenomenon. A 4 × 4 grid of Gauss points is used for the fully integrated terms, while a 2 × 2 grid is used for the reduced ones. The prestresses' virtual work reads

$$ \delta L_{\sigma 0} = \delta \mathbf{q}_{\tau i}^{T} \mathbf{K}_{\sigma 0}^{\tau sij} \mathbf{q}_{sj}, \qquad (31) $$

where $\mathbf{K}_{\sigma 0}^{\tau sij}$ is the fundamental nucleus of the geometric stiffness matrix. More details on the substitutions needed to obtain Equations (28), (30), and (31) are reported in Appendix A, whereas Appendix B shows the components of the stiffness matrices' fundamental nuclei.

## 3. Numerical Results

This section presents some numerical investigations to assess the proposed Finite Element formulation. Three benchmark cases are considered: a monolayer plate, a multilayer plate, and a multilayer plate with a central circular cut-out. For all the cases, a square plate ($a = b$) is considered, whereas different boundary and load conditions are applied in order to obtain a variety of results as wide as possible. Two materials are used whose properties are represented in Table 1 where 'L' and 'T' stand for longitudinal and transverse direction, respectively.

<table><caption>**Table 1.** Material properties.</caption><thead><tr><th>Case</th><th><i>E</i><sub>L</sub> [GPa]</th><th><i>E</i><sub>T</sub> [GPa]</th><th><i>G</i><sub>LT</sub> [GPa]</th><th><i>G</i><sub>TT</sub> [GPa]</th><th><i>ν</i><sub>LT</sub>, <i>ν</i><sub>TT</sub></th></tr></thead><tbody><tr><td>1</td><td>50</td><td>10</td><td>5</td><td>5</td><td>0.25</td></tr><tr><td>2, 3</td><td>181</td><td>10.27</td><td>7.17</td><td>3.78</td><td>0.28</td></tr></tbody></table>
Abaqus 3D models are developed as reference solutions. For these models, the mesh is constituted by quadratic solid elements with reduced integration and three degrees of freedom per node (C3D20R). A refined in-plane mesh is needed to obtain accurate results in Abaqus due to the element-wise constant orientation of the fibers. Nine-node square elements (QUAD9) are used for CUF solutions. The results are presented in terms of an equivalent critical force $F_i^{cr}$ for the ith buckling mode defined as

$$F_i^{cr} = \lambda_i ah, \qquad (32)$$

where $\lambda_i$ is the ith buckling eigenvalue representative of an applied lateral surface load (whose units are [Pa]).

The following equation is used to compute the percentage errors:

$$\delta_i = \frac{|F_{i,\text{CUF}}^{cr} - F_{i,\text{Abaqus3D}}^{cr}|}{|F_{i,\text{Abaqus3D}}^{cr}|} \cdot 100. \qquad (33)$$

### 3.1. Monolayer Plate

The first benchmark problem is a monolayer plate characterized by the following dimensions: $a = b = 1$ m, $h = 0.01$ m. The fiber angle is represented as a function of $y'$. In this case, the local frame of the fiber path is aligned with the global reference frame of the structure. Hence, the characteristic length in Equation (9) is set as $d = b$. Figure 4 represents the law of the fibers' angle, which can be written as $90 < 0, 90 >$. This law is taken from Viglietti et al. [25], where it is applied on a rectangular plate for vibration analyses.

<figure><img src="image_10.png" /><figcaption>Figure 4. In-plane fiber variation path, case 1.</figcaption></figure>

The plate is clamped in correspondence to the side where $y/b = 0$, while only the displacement $u_y$ is free at $y = b$. A constant pressure $P_{yy}^0 = 1$ Pa is applied at $y = b$, as shown in Figure 5.

<figure><img src="image_11.png" /><figcaption>Figure 5. Boundary and loading conditions, case 1.</figcaption></figure>
The Abaqus reference solution contains 80 elements along each side and 16 elements along the thickness. Table 2 shows the results of the preliminary convergence analysis, considering the ED2 model. The first four buckling loads are shown for each mesh, together with the number of degrees of freedom (DOF). It is possible to observe that by increasing the number of in-plane elements, the results get progressively closer to the reference solution. The 8 × 8 mesh is considered the best compromise between the computational cost and results' accuracy, and for this reason, it is used in the following analyses in order to study the behavior of higher-order CUF theories.

<table><caption>**Table 2.** Comparison of DOF and critical loads $F_i^{cr}$ [N] between Abaqus 3D and ED2 mesh, case 1.</caption><thead><tr><th></th><th>DOF</th><th>$F_1^{cr} \cdot 10^{-3}$</th><th>$F_2^{cr} \cdot 10^{-3}$</th><th>$F_3^{cr} \cdot 10^{-3}$</th><th>$F_4^{cr} \cdot 10^{-3}$</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>1,310,499</td><td>28.926</td><td>57.366</td><td>85.585</td><td>122.498</td></tr><tr><td>2 × 2</td><td>225</td><td>29.901</td><td>57.475</td><td>136.926</td><td>170.607</td></tr><tr><td>4 × 4</td><td>729</td><td>28.126</td><td>55.545</td><td>88.768</td><td>125.994</td></tr><tr><td>6 × 6</td><td>1521</td><td>28.333</td><td>56.023</td><td>85.550</td><td>122.138</td></tr><tr><td>8 × 8</td><td>2601</td><td>28.473</td><td>56.323</td><td>85.196</td><td>121.744</td></tr><tr><td>10 × 10</td><td>3969</td><td>28.565</td><td>56.524</td><td>85.160</td><td>121.739</td></tr></tbody></table>

Table 3 shows the DOF for different theories and expansion orders. FSDT and CLT show the smallest number of DOF. Because of the way in which they are modeled through CUF, FSDT and CLT theories have the same number of DOF. It is possible to observe that mixed CUF models, which can be considered the computationally most expensive ones, are characterized by a number of DOF that is two magnitude orders smaller than the Abaqus reference solution.

<table><caption>**Table 3.** Number of degrees of freedom, case 1.</caption><thead><tr><th>Model</th><th>DOF</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>1,310,499</td></tr><tr><td>3LM4</td><td>22,542</td></tr><tr><td>3LM2</td><td>12,138</td></tr><tr><td>3LD4</td><td>11,271</td></tr><tr><td>3LD2</td><td>6069</td></tr><tr><td>ED4</td><td>4335</td></tr><tr><td>ED2</td><td>2601</td></tr><tr><td>FSDT</td><td>1734</td></tr><tr><td>CLT</td><td>1734</td></tr></tbody></table>

In order to retrieve the plate nonuniform stress field when a compression load is applied, a prebuckling analysis is performed. Figure 6 shows the in-plane contour plots of the in-plane displacement $u_y^0$ and the in-plane normal stresses $\sigma_{xx}^0$ and $\sigma_{yy}^0$. The comparison between the 3LM4 and Abaqus 3D results can be observed, showing a good agreement between the two approaches. In Figure 6a, the displacement along $y$ is negative because the plate is compressed. The maximum value is zero in correspondence to the clamped side, while the minimum value equal to $-7.37 \times 10^{-12}$ m is reached at the top-right corner. As shown in Figure 6b, the maximum and minimum values of $\sigma_{yy}$ are $-0.585$ Pa and $-4.2$ Pa, respectively. It is possible to observe that the normal stress along $y$ is almost constant because of the loading and boundary conditions that are considered in this case. Figure 6b shows that the maximum and minimum values of $\sigma_{xx}$ are $0.083$ Pa and $-0.331$ Pa, respectively. Since the plate is compressed along $y$, the normal stress along $x$ is one order of magnitude smaller than $\sigma_{yy}$. It is possible to observe that the minimum value in this case is reached at the top side of the plate.
<figure><img src="image_12.png" /><figcaption>(a) Displacement along y-axis $u_y^0$ [m].</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) In-plane normal stress $\sigma_{yy}^0$ [Pa].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) In-plane normal stress $\sigma_{xx}^0$ [Pa].</figcaption></figure>

<caption>**Figure 6.** Contour plots, comparison between the 3LM4 model (**left**) and Abaqus 3D (**right**), case 1.</caption>

Table 4 shows the first four critical loads for different theories. The CLT is able to predict the buckling loads of the plate, but it shows the lowest level of accuracy among the considered theories. The 3LD2 and 3LD4 models show a similar approximation since they both show an error of 0.2% on the first buckling load. The best approximation of the first critical load is given by the 3LM2 and 3LM4 models, which show an error of 0.1%. As far as the higher critical loads are concerned, the LM models are closer to the reference solution.
<table><caption>**Table 4.** Critical loads $F_i^{cr}$ [N], case 1.</caption><thead><tr><th></th><th>F<sub>1</sub><sup>cr</sup> · 10<sup>-3</sup></th><th>F<sub>2</sub><sup>cr</sup> · 10<sup>-3</sup></th><th>F<sub>3</sub><sup>cr</sup> · 10<sup>-3</sup></th><th>F<sub>4</sub><sup>cr</sup> · 10<sup>-3</sup></th></tr></thead><tbody><tr><td>Abaqus</td><td>28.926</td><td>57.366</td><td>85.585</td><td>122.498</td></tr><tr><td>3D</td><td>28.962</td><td>57.319</td><td>85.796</td><td>122.756</td></tr><tr><td>3LM4</td><td>28.958</td><td>57.309</td><td>85.788</td><td>122.743</td></tr><tr><td>3LM2</td><td>28.979</td><td>57.478</td><td>85.891</td><td>122.928</td></tr><tr><td>3LD4</td><td>28.979</td><td>57.479</td><td>85.892</td><td>122.929</td></tr><tr><td>3LD2</td><td>28.468</td><td>56.296</td><td>85.154</td><td>121.683</td></tr><tr><td>ED4</td><td>28.473</td><td>56.323</td><td>85.196</td><td>121.744</td></tr><tr><td>ED2</td><td>28.973</td><td>57.492</td><td>85.918</td><td>122.976</td></tr><tr><td>FSDT</td><td>29.008</td><td>57.681</td><td>86.179</td><td>123.361</td></tr><tr><td>CLT</td><td>28.962</td><td>57.319</td><td>85.796</td><td>122.756</td></tr></tbody></table>

## 3.2. Multilayer Plate

The second benchmark problem that is considered is a multilayer plate, and it was taken from Hao et al. [7]. It has also been studied by Fallahi et al. [26]. The plate is square and composed of sixteen layers with equal thickness. The in-plane dimensions are $a = b = 254$ mm, and each ply has a thickness of 0.15 mm. The fiber angle is expressed as a function of the $x'$-axis. Axes $x'$ and $y'$ of the local reference system of the fiber path are parallel to axes $x$ and $y$ of the global reference system of the plate, but they are translated such that their origin is located at the center of the plate ($a/2, b/2$). In this case, the characteristic length in Equation (9) is $a/2$, hence $d = a/2$. A symmetric and balanced stack is considered with the following fiber path parameters: $[0 < 60, 15 > |0 < -60, -15 > |0 < -60, -15 > |0 < 60, 15 >]_4$. The stacking sequence is represented in Figure 7 for the representative layers $0 < 60, 15 >$ and $0 < 60, -15 >$ in the symmetric pattern.

<figure><img src="image_15.png" /><figcaption>**Figure 7.** In-plane fiber path and stacking sequence, case 2.</figcaption></figure>

As for the previous case, the Abaqus reference solution contains 80 elements along each side and 16 elements along the thickness. For the CUF results, an 8 × 8 mesh is considered. At each side of the plate, the displacement along the z-axis is equal to zero ($u_z = 0$). At the lower-left corner (0,0), the displacement along $y$ is also zero ($u_y = 0$), while at the lower-right corner ($a,0$), all the displacements are constrained ($u_x = u_y = 0$). Two loading conditions are considered: a pure compression case along the x-axis ($P_{xx}^0 = 1$ Pa) and a combined compression-shear case ($P_{xx}^0 = P_{xy}^0 = 1$ Pa). The boundary and loading conditions are shown in Figure 8.

Table 5 shows the first four buckling loads for the pure compression and compression-shear cases. Considering pure compression, the model using CLT is able to approximate the results given by Abaqus 3D. Nevertheless, this theory shows a percentage error of 1.6% for the first critical load, which grows up to 3.1% for the fourth critical load. The error on $F_i^{cr}$ can be reduced to 0.04% considering the ED4 model. Through LW models, this error can be lower than 0.001%. In this case, the best approximation of $F_i^{cr}$ is given by the LM4 model.
<figure><img src="image_16.png" /><figcaption>**Figure 8.** Boundary and loading conditions, case 2.</figcaption></figure>

<table><caption>**Table 5.** Critical loads $F_i^{cr}$ [N], pure compression ($P_{xx}^0 = 1$ Pa), and compression-shear ($P_{xx}^0 = P_{xy}^0 = 1$ Pa) loads, case 2.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Pure compression P<sup>0</sup><sub>xx</sub> = 1 Pa</th></tr><tr><th>F<sub>1</sub><sup>cr</sup> · 10<sup>-3</sup></th><th>F<sub>2</sub><sup>cr</sup> · 10<sup>-3</sup></th><th>F<sub>3</sub><sup>cr</sup> · 10<sup>-3</sup></th><th>F<sub>4</sub><sup>cr</sup> · 10<sup>-3</sup></th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>13.63</td><td>21.57</td><td>35.42</td><td>54.46</td></tr><tr><td>Ref. [7]</td><td>13.63</td><td>21.64</td><td>35.41</td><td>54.56</td></tr><tr><td>Ref. [26]</td><td>13.67</td><td>21.68</td><td>35.69</td><td>54.60</td></tr><tr><td>LM4</td><td>13.63</td><td>21.57</td><td>35.62</td><td>54.47</td></tr><tr><td>LM2</td><td>13.63</td><td>21.57</td><td>35.62</td><td>54.47</td></tr><tr><td>LD4</td><td>13.63</td><td>21.57</td><td>35.62</td><td>54.47</td></tr><tr><td>LD2</td><td>13.63</td><td>21.57</td><td>35.62</td><td>54.47</td></tr><tr><td>ED4</td><td>13.64</td><td>21.71</td><td>35.80</td><td>54.60</td></tr><tr><td>ED2</td><td>13.67</td><td>21.74</td><td>35.85</td><td>54.81</td></tr><tr><td>FSDT</td><td>13.68</td><td>21.75</td><td>35.86</td><td>54.83</td></tr><tr><td>CLT</td><td>13.85</td><td>22.04</td><td>36.42</td><td>56.17</td></tr></tbody></table>

<table><thead><tr><th rowspan="2"></th><th colspan="4">Compression-shear P<sup>0</sup><sub>xx</sub> = P<sup>0</sup><sub>xy</sub> = 1 Pa</th></tr><tr><th>F<sub>1</sub><sup>cr</sup> · 10<sup>-3</sup></th><th>F<sub>2</sub><sup>cr</sup> · 10<sup>-3</sup></th><th>F<sub>3</sub><sup>cr</sup> · 10<sup>-3</sup></th><th>F<sub>4</sub><sup>cr</sup> · 10<sup>-3</sup></th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>12.04</td><td>18.46</td><td>30.86</td><td>41.87</td></tr><tr><td>Ref. [7]</td><td>12.04</td><td>18.49</td><td>30.83</td><td>41.84</td></tr><tr><td>LM4</td><td>12.05</td><td>18.51</td><td>31.03</td><td>42.39</td></tr><tr><td>LM2</td><td>12.05</td><td>18.51</td><td>31.03</td><td>42.39</td></tr><tr><td>LD4</td><td>12.05</td><td>18.51</td><td>31.03</td><td>42.39</td></tr><tr><td>LD2</td><td>12.05</td><td>18.51</td><td>31.03</td><td>42.39</td></tr><tr><td>ED4</td><td>12.07</td><td>18.55</td><td>31.14</td><td>42.56</td></tr><tr><td>ED2</td><td>12.08</td><td>18.58</td><td>31.20</td><td>42.75</td></tr><tr><td>FSDT</td><td>12.10</td><td>18.59</td><td>31.22</td><td>42.73</td></tr><tr><td>CLT</td><td>12.26</td><td>18.84</td><td>31.75</td><td>43.70</td></tr></tbody></table>

Considering the compression-shear loading condition, CLT shows a higher maximum error of 4.4% in correspondence to $F_4^{cr}$. This error is reduced to 1.2% with LD4 and LM4 theories. The minimum error on the first buckling load is 0.090% and is obtained by an LM4 model. It is possible to notice that, because of the combined load case, lower-order theories show higher errors in comparison with the pure compression load. Furthermore, the addition of the shear load generates a reduction in the eigenvalues obtained by solving the buckling problem. This can be justified considering that, in the compression-shear case, there are two loads applied to the structure, which makes it easier to reach the buckling state. The eigenvalues' reduction is more evident for higher critical loads: considering the Abaqus 3D solution, the first critical load is 11.7% smaller in the compression-shear case, while the fourth critical load has a reduction of 23.1%.
### 3.3. Multilayer Plate with a Central Cut-Out

The third case of analysis is represented by a square multilayer plate with a central circular cut-out. The plate is composed of eight layers, each of them 0.15 mm thick. The in-plane dimensions are $a = b = 254$ mm. The material properties are the same as in the second case. The center of the cut-out is placed at $(a/2, b/2)$, and its radius is $r = 40$ mm. The fiber angle distribution is the same as in the previous case, and the following stack is considered: $[0 < 60, 15 > |0 < -60, -15 > |0 < -60, -15 > |0 < 60, 15 >]_2$. In this case, the Abaqus reference solution is made of 100'352 elements: 6272 elements are defined along the plane of the plate and 16 elements are defined along the thickness. For the CUF results, 72 plate elements were used. The boundary and loading conditions are the same as those presented for case 2. Table 6 shows the first four buckling loads for the pure compression and compression-shear loads. Considering pure compression, it is possible to observe that LM theories give the best approximation of $F_1^{cr}$ and $F_4^{cr}$, while $F_2^{cr}$ and $F_3^{cr}$ are better approximated by LD ones. LM2 shows a minimum error of 0.2% in correspondence to $F_1^{cr}$, while LM4 shows a minimum error of 0.3% in correspondence to $F_4^{cr}$. CLT and FSDT models can predict the buckling loads for this case, even though they present higher errors in comparison with other theories.

<table><caption>**Table 6.** Critical loads $F_i^{cr}$ [N], pure compression ($P_{xx}^0 = 1$ Pa), and compression-shear ($P_{xx}^0 = P_{xy}^0 = 1$ Pa) loads, case 3.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Pure compression P<sup>0</sup><sub>xx</sub> = 1 Pa</th></tr><tr><th>F<sub>1</sub><sup>cr</sup> &middot; 10<sup>-2</sup></th><th>F<sub>2</sub><sup>cr</sup> &middot; 10<sup>-2</sup></th><th>F<sub>3</sub><sup>cr</sup> &middot; 10<sup>-2</sup></th><th>F<sub>4</sub><sup>cr</sup> &middot; 10<sup>-2</sup></th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>13.322</td><td>22.522</td><td>31.974</td><td>36.974</td></tr><tr><td>LM4</td><td>13.294</td><td>22.386</td><td>31.352</td><td>37.071</td></tr><tr><td>LM2</td><td>13.298</td><td>22.390</td><td>31.360</td><td>37.077</td></tr><tr><td>LD4</td><td>13.378</td><td>22.534</td><td>31.661</td><td>37.280</td></tr><tr><td>LD2</td><td>13.383</td><td>22.539</td><td>31.672</td><td>37.288</td></tr><tr><td>ED4</td><td>13.369</td><td>22.554</td><td>31.763</td><td>37.301</td></tr><tr><td>ED2</td><td>13.420</td><td>22.633</td><td>31.916</td><td>37.432</td></tr><tr><td>FSDT</td><td>13.464</td><td>22.658</td><td>31.920</td><td>37.486</td></tr><tr><td>CLT</td><td>13.539</td><td>22.862</td><td>32.353</td><td>37.790</td></tr></tbody><thead><tr><th rowspan="2"></th><th colspan="4">Compression-shear P<sup>0</sup><sub>xx</sub> = P<sup>0</sup><sub>xy</sub> = 1 Pa</th></tr><tr><th>F<sub>1</sub><sup>cr</sup> &middot; 10<sup>-2</sup></th><th>F<sub>2</sub><sup>cr</sup> &middot; 10<sup>-2</sup></th><th>F<sub>3</sub><sup>cr</sup> &middot; 10<sup>-2</sup></th><th>F<sub>4</sub><sup>cr</sup> &middot; 10<sup>-2</sup></th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>11.152</td><td>19.396</td><td>23.956</td><td>32.302</td></tr><tr><td>LM4</td><td>11.146</td><td>19.322</td><td>23.491</td><td>32.337</td></tr><tr><td>LM2</td><td>11.149</td><td>19.326</td><td>23.496</td><td>32.343</td></tr><tr><td>LD4</td><td>11.204</td><td>19.436</td><td>23.638</td><td>32.519</td></tr><tr><td>LD2</td><td>11.207</td><td>19.440</td><td>23.646</td><td>32.526</td></tr><tr><td>ED4</td><td>11.228</td><td>19.490</td><td>23.730</td><td>32.611</td></tr><tr><td>ED2</td><td>11.253</td><td>19.536</td><td>23.808</td><td>32.709</td></tr><tr><td>FSDT</td><td>11.264</td><td>19.534</td><td>23.808</td><td>32.697</td></tr><tr><td>CLT</td><td>11.326</td><td>19.713</td><td>24.097</td><td>33.020</td></tr></tbody></table>

Also in the compression-shear case, results are closer to the reference solution when LD and LM theories are applied. Even if CLT and FSDT models show a worse approximation of reference solutions, they can correctly match the buckling modes obtained through Abaqus 3D. The buckling modes for the compression-shear case are presented in Figure 9, where a comparison between the Abaqus 3D and LM4 model is shown. In both the presented solutions and the Abaqus model, the eigenvectors are normalized such that for each eigenvalue, the maximal displacement value of the corresponding mode is equal to one and the scale used in all the presented modal shapes is the same. In comparison with the pure compression case, these mode shapes are more complex because of the combined loading condition. Therefore, higher-order models are needed in order to predict the buckling in a more accurate way.
<figure><img src="image_17.png" /><figcaption>**Figure 9.** Buckling modes, compression-shear load $P_{xx}^0 = P_{xy}^0 = 1$ Pa, and comparison between LM2 model (**left**) and Abaqus 3D (**right**), case 3.</figcaption></figure>

**4. Conclusions**

In this study, linear buckling analyses of VAT plates were performed by developing two-dimensional CUF Finite Elements within Reissner's Mixed Variational Theorem. Thanks to its adaptability to a wide range of structural problems and the possibility to choose the approximation theory a priori, CUF shows a good prediction of buckling loads with an acceptable computational cost. Even though only linear laws were considered in the numerical simulations for the variation in the fiber angle, the implemented numerical framework allows one to address a generic variation law. Three numerical examples of plates containing curvilinear fibers were analyzed, and the results were compared to Abaqus 3D and literature reference solutions. The comparison among the various results allows us to distinguish the developed models according to their accuracy and numerical complexity. Even though CLT and FSDT show a reduced amount of degrees of freedom, allowing one to correctly predict
buckling loads for the considered benchmark structures 2 and 3, they show a higher loss of accuracy when plates that are not symmetric and balanced are considered, as observed in the first benchmark case. Increasing the through-the-thickness polynomial order, general ESL displacement theories are obtained, which allow one to improve the results' accuracy. A further improvement in the prediction of buckling loads is obtained through the development of LW displacement models. ESL and LW theories can correctly predict Abaqus 3D results with a limited number of DOF, but being derived from a PVD variational statement, they show a lack of accuracy for the out-of-plane stresses, which strongly influence the evaluation of the geometric stiffness matrix needed for the linear buckling analysis. Hence, the derivation of LW models to RMVT allows one to overcome this problem, enhancing the solution accuracy, especially in the case of the first critical load. Even if LM models allow one to improve the prediction of the through-the-thickness behavior, they are characterized by the maximum computational cost among the two-dimensional theories developed in this study, and yet their application is still advantageous mainly when plates characterized by complex stacking sequences are considered, especially when compared to the cumbersome three-dimensional simulations. In summary, the utilization of RMVT in the context of CUF has demonstrated the potential to enhance the precision and efficiency of modeling VAT plates for buckling analyses. Nevertheless, the scope of this approach extends beyond plate analysis, as there are promising prospects for employing, as a future perspective, this framework in the examination of more challenging VAT structures such as shells.

**Author Contributions:** Methodology, G.G., D.A.I., L.K. and M.M.; software, D.A.I.; validation, D.A.I.; writing—original draft, D.A.I.; writing—review and editing, G.G., L.K. and M.M.; supervision, G.G. and M.M.; funding acquisition, G.G. and M.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Luxembourg National Research Fund (FNR), grant reference INTER/ANR/21/16215936 GLAMOUR-VSC. For the purpose of open access, and in fulfillment of the obligations arising from the grant agreement, the authors have applied a Creative Commons Attribution 4.0 International (CC BY 4.0) license to any author-accepted manuscript version arising from this submission. M. Montemurro is grateful to the French National Research Agency for supporting this work through the research project GLAMOUR-VSC (Global-LocAl two-level Multi-scale optimisation strategy accOUnting for pRocess-induced singularities to design Variable Stiffness Composites) ANR-21-CE10-0014.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data will be made available upon request.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## Appendix A. Governing Equations

Through the substitution of Equations (3), (7), (27) and (11) into Equation (18), the PVD governing equations are obtained:

$$ \begin{aligned} & \int_{\Omega} \delta q_{\tau i}^{T} (D_{p}^{T}(N_{i} \mathbf{I}) \tilde{Z}_{pp}^{\tau s} D_{p}(N_{j} \mathbf{I}) + D_{p}^{T}(N_{i} \mathbf{I}) \tilde{Z}_{pn}^{\tau s} D_{n \Omega}(N_{j} \mathbf{I}) + D_{p}^{T}(N_{i} \mathbf{I}) \tilde{Z}_{pn}^{\tau s, z}(N_{j} \mathbf{I}) + \\ & + D_{n \Omega}^{T}(N_{i} \mathbf{I}) \tilde{Z}_{np}^{\tau s} D_{p}(N_{j} \mathbf{I}) + D_{n \Omega}^{T}(N_{i} \mathbf{I}) \tilde{Z}_{nn}^{\tau s} D_{n \Omega}(N_{j} \mathbf{I}) + D_{n \Omega}^{T}(N_{i} \mathbf{I}) \tilde{Z}_{nn}^{\tau s, z}(N_{j} \mathbf{I}) + \\ & + (N_{i} \mathbf{I}) \tilde{Z}_{np}^{\tau, z s} D_{p}(N_{j} \mathbf{I}) + (N_{i} \mathbf{I}) \tilde{Z}_{nn}^{\tau, z s} D_{n \Omega}(N_{j} \mathbf{I}) + (N_{i} \mathbf{I}) \tilde{Z}_{nn}^{\tau, z s, z}(N_{j} \mathbf{I})) q_{sj} d\Omega = \delta L_{\nu 0}, \end{aligned} \quad (A1) $$

where

$$ (\tilde{Z}_{wr}^{\tau s}, \tilde{Z}_{wr}^{\tau, z s}, \tilde{Z}_{wr}^{\tau s, z}, \tilde{Z}_{wr}^{\tau, z s, z}) = (\tilde{C}_{wr} E_{\tau s}, \tilde{C}_{wr} E_{\tau, z s}, \tilde{C}_{wr} E_{\tau s, z}, \tilde{C}_{wr} E_{\tau, z s, z}) : w, r = p, n, \quad (A2) $$

$$ (E_{\tau s}, E_{\tau, z s}, E_{\tau s, z}, E_{\tau, z s, z}) = \int_{h} (F_{\tau} F_{s}, F_{\tau, z} F_{s}, F_{\tau} F_{s, z}, F_{\tau, z} F_{s, z}) dz. \quad (A3) $$
Index 'z', when preceded by a comma, refers to the derivative along the z-axis.

In a similar way, through the substitution of Equations (3), (13), (27), (29), and (12) into Equation (18), the governing equations of the RMVT can be written as follows:

$$ \int_{\Omega} \delta \boldsymbol{q}_{\tau i}^{\mathrm{T}} [\boldsymbol{D}_{p}^{\mathrm{T}}(N_{i} \boldsymbol{I}) \hat{\boldsymbol{Z}}_{pp}^{\mathrm{ts}} \boldsymbol{D}_{p}(N_{j} \boldsymbol{I})] \boldsymbol{q}_{sj} + \delta \boldsymbol{q}_{\tau i}^{\mathrm{T}} [\boldsymbol{D}_{p}^{\mathrm{T}}(N_{i} \boldsymbol{I}) \hat{\boldsymbol{Z}}_{pn}^{\mathrm{ts}}(N_{j} \boldsymbol{I}) + \boldsymbol{D}_{n \Omega}^{\mathrm{T}}(N_{i} \boldsymbol{I})(E_{\tau s} \boldsymbol{I})(N_{j} \boldsymbol{I}) + \\ + (N_{i} \boldsymbol{I})(E_{\tau z s} \boldsymbol{I})(N_{j} \boldsymbol{I})] \boldsymbol{g}_{sj} + \delta \boldsymbol{g}_{\tau i}^{\mathrm{T}} [(N_{i} \boldsymbol{I})(E_{\tau s} \boldsymbol{I}) \boldsymbol{D}_{n \Omega}(N_{j} \boldsymbol{I}) + (N_{i} \boldsymbol{I})(E_{\tau s z} \boldsymbol{I})(N_{j} \boldsymbol{I}) + \\ -(N_{i} \boldsymbol{I}) \hat{\boldsymbol{Z}}_{np}^{\mathrm{ts}} \boldsymbol{D}_{p}(N_{j} \boldsymbol{I})] \boldsymbol{q}_{sj} - \delta \boldsymbol{g}_{\tau i}^{\mathrm{T}}(N_{i} \boldsymbol{I}) \hat{\boldsymbol{Z}}_{nn}^{\mathrm{ts}}(N_{j} \boldsymbol{I}) \boldsymbol{g}_{sj} d\Omega = \delta L_{\sigma 0}, \quad (A4) $$

where

$$ (\hat{\boldsymbol{Z}}_{wr}^{\mathrm{ts}}, \hat{\boldsymbol{Z}}_{wr}^{\tau z s}, \hat{\boldsymbol{Z}}_{wr}^{\tau s z}, \hat{\boldsymbol{Z}}_{wr}^{\tau z s z}) = (\hat{\boldsymbol{C}}_{wr} E_{\tau s}, \hat{\boldsymbol{C}}_{wr} E_{\tau z s}, \hat{\boldsymbol{C}}_{wr} E_{\tau s z}, \hat{\boldsymbol{C}}_{wr} E_{\tau z s z}) : w, r = p, n. \quad (A5) $$

Through the substitution of Equations (5) and (27) into Equation (17), it is possible to obtain the prestresses' virtual work:

$$ \delta L_{\sigma 0} = \int_{\Omega} \int_{h} \boldsymbol{q}_{\tau i}^{\mathrm{T}} [F_{\tau} F_{s} N_{i,x} N_{j,y} \sigma_{xx}^{0} + F_{\tau} F_{s} N_{i,y} N_{j,y} \sigma_{yy}^{0} + F_{\tau} F_{s} (N_{i,x} N_{j,y} + N_{i,y} N_{j,x}) \sigma_{xy}^{0} + \\ + (F_{\tau} F_{s,z} N_{i,x} N_{j} + F_{\tau,z} F_{s} N_{i} N_{j,x}) \sigma_{xz}^{0} + (F_{\tau} F_{s,z} N_{i,y} N_{j} + F_{\tau,z} F_{s} N_{i} N_{j,y}) \sigma_{yz}^{0} + \\ + F_{\tau,z} F_{s,z} N_{i} N_{j} \sigma_{zz}^{0}] \boldsymbol{I} \boldsymbol{q}_{sj} dz d\Omega. \quad (A6) $$

## Appendix B. Expression of the Fundamental Nuclei

This appendix reports the FN of the structure stiffness matrix for the PVD and RMVT variational statements. The components of the FN for the PVD case can be written as follows in the case of orthotropic material:

$$ \begin{aligned} K_{xx}^{\mathrm{tsij}} &= \int_{\Omega} (\tilde{Z}_{pp11}^{\mathrm{ts}} N_{j,x} N_{i,x} + \tilde{Z}_{pp16}^{\mathrm{ts}} N_{j,y} N_{i,x} + \tilde{Z}_{pp16}^{\mathrm{ts}} N_{j,x} N_{i,y} + \tilde{Z}_{pp66}^{\mathrm{ts}} N_{j,y} N_{i,y} + \tilde{Z}_{nn44}^{\tau,z s,z} N_{j} N_{i}) d\Omega, \\ K_{xy}^{\mathrm{tsij}} &= \int_{\Omega} (\tilde{Z}_{pp12}^{\mathrm{ts}} N_{j,y} N_{i,x} + \tilde{Z}_{pp16}^{\mathrm{ts}} N_{j,x} N_{i,x} + \tilde{Z}_{pp26}^{\mathrm{ts}} N_{j,y} N_{i,y} + \tilde{Z}_{pp66}^{\mathrm{ts}} N_{j,x} N_{i,y} + \tilde{Z}_{nn45}^{\tau,z s,z} N_{j} N_{i}) d\Omega, \\ K_{xz}^{\mathrm{tsij}} &= \int_{\Omega} (\tilde{Z}_{pn13}^{\mathrm{ts,z}} N_{j} N_{i,x} + \tilde{Z}_{pn36}^{\mathrm{ts,z}} N_{j} N_{i,y} + \tilde{Z}_{nn44}^{\tau,z s} N_{j,x} N_{i} + \tilde{Z}_{nn45}^{\tau,z s} N_{j,y} N_{i}) d\Omega, \\ K_{yx}^{\mathrm{tsij}} &= \int_{\Omega} (\tilde{Z}_{pp12}^{\mathrm{ts}} N_{j,x} N_{i,y} + \tilde{Z}_{pp26}^{\mathrm{ts}} N_{j,y} N_{i,y} + \tilde{Z}_{pp16}^{\mathrm{ts}} N_{j,x} N_{i,x} + \tilde{Z}_{pp66}^{\mathrm{ts}} N_{j,y} N_{i,x} + \tilde{Z}_{nn45}^{\tau,z s,z} N_{j} N_{i}) d\Omega, \\ K_{yy}^{\mathrm{tsij}} &= \int_{\Omega} (\tilde{Z}_{pp22}^{\mathrm{ts}} N_{j,y} N_{i,y} + \tilde{Z}_{pp26}^{\mathrm{ts}} N_{j,x} N_{i,y} + \tilde{Z}_{pp26}^{\mathrm{ts}} N_{j,y} N_{i,x} + \tilde{Z}_{pp66}^{\mathrm{ts}} N_{j,x} N_{i,x} + \tilde{Z}_{nn55}^{\tau,z s,z} N_{j} N_{i}) d\Omega, \quad (A7) \\ K_{yz}^{\mathrm{tsij}} &= \int_{\Omega} (\tilde{Z}_{pn23}^{\mathrm{ts,z}} N_{j} N_{i,y} + \tilde{Z}_{pn36}^{\mathrm{ts,z}} N_{j} N_{i,x} + \tilde{Z}_{nn45}^{\tau,z s} N_{j,x} N_{i} + \tilde{Z}_{nn55}^{\tau,z s} N_{j,y} N_{i}) d\Omega, \\ K_{zx}^{\mathrm{tsij}} &= \int_{\Omega} (\tilde{Z}_{nn44}^{\mathrm{ts,z}} N_{j} N_{i,x} + \tilde{Z}_{nn45}^{\mathrm{ts,z}} N_{j} N_{i,y} + \tilde{Z}_{np13}^{\tau,z s} N_{j,x} N_{i} + \tilde{Z}_{np36}^{\tau,z s} N_{j,y} N_{i}) d\Omega, \\ K_{zy}^{\mathrm{tsij}} &= \int_{\Omega} (\tilde{Z}_{nn45}^{\mathrm{ts,z}} N_{j} N_{i,x} + \tilde{Z}_{nn55}^{\mathrm{ts,z}} N_{j} N_{i,y} + \tilde{Z}_{np23}^{\tau,z s} N_{j,y} N_{i} + \tilde{Z}_{pp36}^{\tau,z s} N_{j,x} N_{i}) d\Omega, \\ K_{zz}^{\mathrm{tsij}} &= \int_{\Omega} (\tilde{Z}_{nn44}^{\mathrm{ts}} N_{j,x} N_{i,x} + \tilde{Z}_{nn45}^{\mathrm{ts}} N_{j,y} N_{i,x} + \tilde{Z}_{nn45}^{\mathrm{ts}} N_{j,x} N_{i,y} + \tilde{Z}_{nn55}^{\mathrm{ts}} N_{j,y} N_{i,y} + \tilde{Z}_{nn33}^{\tau,z s,z} N_{j} N_{i}) d\Omega. \end{aligned} $$

The subscripts 'x' and 'y' when preceded by a comma refer to the derivative versus the corresponding in-plane direction. The components of the FN for the RMVT case can be written as follows:
$$
\begin{align*}
K_{uuxx}^{\tau sij} &= \int_{\Omega} (\hat{Z}_{pp11}^{\tau s} N_{j,x} N_{i,x} + \hat{Z}_{pp31}^{\tau s} N_{j,x} N_{i,y} + \hat{Z}_{pp13}^{\tau s} N_{j,y} N_{i,x} + \hat{Z}_{pp33}^{\tau s} N_{j,y} N_{i,y}) d\Omega, \\
K_{uuxy}^{\tau sij} &= \int_{\Omega} (\hat{Z}_{pp12}^{\tau s} N_{j,y} N_{i,x} + \hat{Z}_{pp32}^{\tau s} N_{j,y} N_{i,y} + \hat{Z}_{pp13}^{\tau s} N_{j,x} N_{i,x} + \hat{Z}_{pp33}^{\tau s} N_{j,x} N_{i,y}) d\Omega, \\
K_{uuyx}^{\tau sij} &= \int_{\Omega} (\hat{Z}_{pp21}^{\tau s} N_{j,x} N_{i,y} + \hat{Z}_{pp31}^{\tau s} N_{j,x} N_{i,x} + \hat{Z}_{pp23}^{\tau s} N_{j,y} N_{i,y} + \hat{Z}_{pp33}^{\tau s} N_{j,y} N_{i,x}) d\Omega, \\
K_{uuyy}^{\tau sij} &= \int_{\Omega} (\hat{Z}_{pp22}^{\tau s} N_{j,y} N_{i,y} + \hat{Z}_{pp32}^{\tau s} N_{j,y} N_{i,x} + \hat{Z}_{pp23}^{\tau s} N_{j,x} N_{i,y} + \hat{Z}_{pp33}^{\tau s} N_{j,x} N_{i,x}) d\Omega, \\
K_{uuxz}^{\tau sij} &= 0, \quad K_{uuyz}^{\tau sij} = 0, \quad K_{uuzx}^{\tau sij} = 0, \quad K_{uuzy}^{\tau sij} = 0, \quad K_{uuzz}^{\tau sij} = 0, \\
K_{u\sigma xx}^{\tau sij} &= \int_{\Omega} (E_{\tau,zs} N_j N_i) d\Omega, \quad K_{u\sigma xz}^{\tau sij} = \int_{\Omega} (\hat{Z}_{pn13}^{\tau s} N_j N_{i,x} + \hat{Z}_{pn33}^{\tau s} N_j N_{i,y}) d\Omega, \\
K_{u\sigma yy}^{\tau sij} &= \int_{\Omega} (E_{\tau,zs} N_j N_i) d\Omega, \quad K_{u\sigma yz}^{\tau sij} = \int_{\Omega} (\hat{Z}_{pn23}^{\tau s} N_j N_{i,y} + \hat{Z}_{pn33}^{\tau s} N_j N_{i,x}) d\Omega, \\
K_{u\sigma zx}^{\tau sij} &= \int_{\Omega} (E_{\tau,s} N_j N_{i,x}) d\Omega, \quad K_{u\sigma zy}^{\tau sij} = \int_{\Omega} (E_{\tau,s} N_j N_{i,y}) d\Omega, \quad K_{u\sigma zz}^{\tau sij} = \int_{\Omega} (E_{\tau,zs} N_j N_i) d\Omega, \\
K_{u\sigma xy}^{\tau sij} &= 0, \quad K_{u\sigma yx}^{\tau sij} = 0, \\
K_{\sigma uuxx}^{\tau sij} &= \int_{\Omega} (E_{\tau,s,z} N_j N_i) d\Omega, \quad K_{\sigma uuzx}^{\tau sij} = \int_{\Omega} (E_{\tau,s} N_{j,x} N_i) d\Omega, \quad K_{\sigma uuyy}^{\tau sij} = \int_{\Omega} (E_{\tau,s,z} N_j N_i) d\Omega, \\
K_{\sigma uuyz}^{\tau sij} &= \int_{\Omega} (E_{\tau,s} N_{j,y} N_i) d\Omega, \quad K_{\sigma uuzx}^{\tau sij} = -\int_{\Omega} (\hat{Z}_{np31}^{\tau s} N_{j,x} N_i - \hat{Z}_{np33}^{\tau s} N_{j,y} N_i) d\Omega, \\
K_{\sigma uzy}^{\tau sij} &= -\int_{\Omega} (\hat{Z}_{np32}^{\tau s} N_{j,y} N_i - \hat{Z}_{np33}^{\tau s} N_{j,x} N_i) d\Omega, \quad K_{\sigma uzz}^{\tau sij} = \int_{\Omega} (E_{\tau,s,z} N_j N_i) d\Omega, \\
K_{\sigma uxy}^{\tau sij} &= 0, \quad K_{\sigma uyx}^{\tau sij} = 0, \\
K_{\sigma\sigma uxx}^{\tau sij} &= -\int_{\Omega} (\hat{Z}_{nn11}^{\tau s} N_j N_i) d\Omega, \quad K_{\sigma\sigma xy}^{\tau sij} = -\int_{\Omega} (\hat{Z}_{nn12}^{\tau s} N_j N_i) d\Omega, \\
K_{\sigma\sigma yx}^{\tau sij} &= -\int_{\Omega} (\hat{Z}_{nn21}^{\tau s} N_j N_i) d\Omega, \quad K_{\sigma\sigma xx}^{\tau sij} = -\int_{\Omega} (\hat{Z}_{nn22}^{\tau s} N_j N_i) d\Omega, \\
K_{\sigma\sigma zx}^{\tau sij} &= 0, \quad K_{\sigma\sigma yz}^{\tau sij} = 0, \quad K_{\sigma\sigma zx}^{\tau sij} = 0, \quad K_{\sigma\sigma zy}^{\tau sij} = 0, \quad K_{\sigma\sigma zz}^{\tau sij} = 0.
\end{align*}
\tag{A8}
$$

References

1. Dirk, H.J.L.; Ward, C.; Potter, K.D. The engineering aspects of automated prepreg layup: History, present and future. *Compos. Part B Eng.* **2012**, *43*, 997–1009.

2. Zhuo, P.; Li, S.; Ashcroft, I.A.; Jones, A.I. Material extrusion additive manufacturing of continuous fiber reinforced polymer matrix composites: A review and outlook. *Compos. Part B Eng.* **2021**, *224*, 109143. [CrossRef]

3. Hyer, M.W.; Lee, H.H. *Innovative Design of Composite Structures: The Use of Curvilinear Fiber Format to Improve Buckling Resistance of Composite Plates with Central Circular Holes*; Technical Report; NASA Langley Research Center: Hampton, VA, USA, 1990.

4. Brampton, C.J.; Wu, K.C.; Kim, H.A. New optimization method for steered fiber composites using the level set method. *Struct. Multidiscip. Optim.* **2015**, *52*, 493–505. [CrossRef]

5. Gürdal, Z.; Tatting, B.F.; Wu, C.K. Variable stiffness composite panels: Effects of stiffness variation on the in-plane and buckling response. *Compos. Part A Appl. Sci. Manuf.* **2008**, *39*, 911–922. [CrossRef]

6. Oliveri, V.; Milazzo, A. A Rayleigh-Ritz approach for postbuckling analysis of variable angle tow composite stiffened panels. *Comput. Struct.* **2018**, *196*, 263–276. [CrossRef]

7. Hao, P.; Yuan, X.; Liu, H.; Wang, B.; Liu, C.; Yang, D.; Zhan, S. Isogeometric buckling analysis of composite variable-stiffness panels. *Compos. Struct.* **2017**, *165*, 192–208. [CrossRef]

8. Hao, P.; Liu, C.; Yuan, X.; Wang, B.; Li, G.; Zhu, T.; Niu, F. Buckling optimization of variable-stiffness composite panels based on flow field function. *Compos. Struct.* **2017**, *181*, 240–255. [CrossRef]

9. Raju, G.; Wu, Z.; Kim, B.C.; Weaver, P.M. Prebuckling and buckling analysis of variable angle tow plates with general boundary conditions. *Compos. Struct.* **2012**, *94*, 2961–2970. [CrossRef]

10. Sciascia, G.; Oliveri, V.; Weaver, P.M. Eigenfrequencies of prestressed variable stiffness composite shells. *Compos. Struct.* **2021**, *270*, 114019. [CrossRef]

11. Chen, X.; Wu, Z.; Nie, G.; Weaver, P. Buckling analysis of variable angle tow composite plates with a through-the-width or an embedded rectangular delamination. *Int. J. Solids Struct.* **2018**, *138*, 166–180. [CrossRef]

12. Wu, Z.; Raju, G.; Weaver, P.M. Framework for the buckling optimization of variable-angle tow composite plates. *AIAA J.* **2015**, *53*, 3788–3804. [CrossRef]
13. Montemurro, M. An extension of the polar method to the first-order shear deformation theory of laminates. *Compos. Struct.* **2015**, *127*, 328–339. [CrossRef]

14. Montemurro, M. The polar analysis of the third-order shear deformation theory of laminates. *Compos. Struct.* **2015**, *131*, 775–789. [CrossRef]

15. Montemurro, M.; Catapano, A. On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates. *Compos. Struct.* **2017**, *161*, 145–159. [CrossRef]

16. Catapano, A.; Montemurro, M.; Balcou, J.-A.; Panettieri, E. Rapid prototyping of variable angle-tow composites. *Aerotec. Missili Spaz.* **2019**, *98*, 257–271. [CrossRef]

17. Montemurro, M.; Catapano, A. A general b-spline surfaces theoretical framework for optimisation of variable angle-tow laminates. *Compos. Struct.* **2019**, *209*, 561–578. [CrossRef]

18. Fiordilino, G.A.; Izzi, M.I.; Montemurro, M. A general isogeometric polar approach for the optimisation of variable stiffness composites: Application to eigenvalue buckling problems. *Mech. Mater.* **2021**, *153*, 103574. [CrossRef]

19. Carrera, E. Theories and finite elements for multilayered, anisotropic, composite plates and shells. *Arch. Comput. Methods Eng.* **2002**, *9*, 87–140. [CrossRef]

20. Carrera, E. Theories and finite elements for multilayered plates and shells: A unified compact formulation with numerical assessment and benchmarking. *Arch. Comput. Methods Eng.* **2003**, *10*, 215–296. [CrossRef]

21. Carrera, E.; Giunta, G.; Brischetto, S. Hierarchical closed form solutions for plates bent by localized transverse loadings. *J. Zhejiang-Univ.-Sci. A* **2007**, *8*, 1026–1037. [CrossRef]

22. Gaetano, G.; Salim, B.; Fabio, B.; Erasmo, C. Hierarchical theories for a linearised stability analysis of thin-walled beams with open and closed cross-section. *Adv. Aircr. Spacecr. Sci.* **2014**, *1*, 253.

23. Hui, Y.; Xu, R.; Giunta, G.; De Pietro, G.; Hu, H.; Belouettar, S.; Carrera, E. Multiscale CUF-FE2 nonlinear analysis of composite beam structures. *Comput. Struct.* **2019**, *221*, 28–43. [CrossRef]

24. Hui, Y.; Bai, X.; Yang, Y.; Yang, J.; Huang, Q.; Liu, X.; Huang, W.; Giunta, G.; Hu, H. A data-driven CUF-based beam model based on the tree-search algorithm. *Compos. Struct.* **2022**, *300*, 116123. [CrossRef]

25. Viglietti, A.; Zappino, E.; Carrera, E. Analysis of variable angle tow composites structures using variable kinematic models. *Compos. Part B Eng.* **2019**, *171*, 272–283. [CrossRef]

26. Fallahi, N.; Viglietti, A.; Carrera, E.; Pagani, A.; Zappino, E. Effect of fiber orientation path on the buckling, free vibration, and static analyses of variable angle tow panels. *Facta Univ. Ser. Mech. Eng.* **2020**, *18*, 165–188. [CrossRef]

27. Fallahi, N.; Carrera, E.; Pagani, A. Application of ga optimization in analysis of variable stiffness composites. *Int. J. Mater. Metall. Eng.* **2021**, *15*, 65–70.

28. Sanchez-Majano, A.R.; Pagani, A.; Petrolo, M.; Zhang, C. Buckling sensitivity of tow-steered plates subjected to multiscale defects by high-order finite elements and polynomial chaos expansion. *Materials* **2021**, *14*, 2706. [CrossRef] [PubMed]

29. Pagani, A.; Sanchez-Majano, A.R. Influence of fiber misalignments on buckling performance of variable stiffness composites using layerwise models and random fields. *Mech. Adv. Mater. Struct.* **2022**, *29*, 384–399. [CrossRef]

30. Demasi, L.; Ashenafi, Y.; Cavallaro, R.; Santarpia, E. Generalized unified formulation shell element for functionally graded variable-stiffness composite laminates and aeroelastic applications. *Compos. Struct.* **2015**, *131*, 501–515. [CrossRef]

31. Vescovini, R.; Dozio, L. A variable-kinematic model for variable stiffness plates: Vibration and buckling analysis. *Compos. Struct.* **2016**, *142*, 15–26. [CrossRef]

32. Carrera, E.; Demasi, L. Classical and advanced multilayered plate elements based upon pvd and rmvt. part 1: Derivation of finite element matrices. *Int. J. Numer. Methods Eng.* **2002**, *55*, 191–231. [CrossRef]

33. Carrera, E.; Demasi, L. Classical and advanced multilayered plate elements based upon pvd and rmvt. part 2: Numerical implementations. *Int. J. Numer. Methods Eng.* **2002**, *55*, 253–291. [CrossRef]

34. Giunta, G.; Iannotta, D.A.; Montemurro, M. A FEM Free Vibration Analysis of Variable Stiffness Composite Plates through Hierarchical Modeling. *Materials* **2023**, *16*, 4643. [CrossRef] [PubMed]

35. Iannotta, D.A.; Giunta, G.; Montemurro, M. A mechanical analysis of variable angle-tow composite plates through variable kinematics models based on Carrera’s unified formulation. *Compos. Struct.* **2024**, *327*, 117717. [CrossRef]

36. Reddy, J.N. *Mechanics of Laminated Composite Plates and Shells: Theory and Analysis*; CRC Press: Boca Raton, FL, USA, 2003.

37. Carrera, E. On the use of the murakami's zig-zag function in the modeling of layered plates and shells. *Comput. Struct.* **2004**, *82*, 541–554. [CrossRef]

38. Bathe, K.-J. *Finite Element Procedures*; Klaus-Jurgen Bathe: Watertown, MA, USA, 2007.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.