{'fraction_non_alphanumeric': 0.03255230849040291, 'fraction_numerical': 0.0146117931869272, 'mean_word_length': 5.111756935270805, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast and accurate elastic analysis of laminated composite plates via isogeometric collocation and an equilibrium-based stress recovery approach

A novel approach which combines isogeometric collocation and an equilibrium-based stress recovery technique is applied to analyze laminated composite plates. Isogeometric collocation is an appealing strong form alternative to standard Galerkin approaches, able to achieve high order convergence rates coupled with a significantly reduced computational cost. Laminated composite plates are herein conveniently modeled considering only one element through the thickness with homogenized material properties. This guarantees accurate results in terms of displacements and in-plane stress components. To recover an accurate out-of-plane stress state, equilibrium is imposed in strong form as a post-processing correction step, which requires the shape functions to be highly continuous. This continuity demand is fully granted by isogeometric analysis properties, and excellent results are obtained using a minimal number of collocation points per direction, particularly for increasing values of length-to-thickness plate ratio and number of layers.


Introduction
Composite materials consist of two or more materials which combined present enhanced properties that could not be acquired employing any of the constituents alone (see, e.g., [28,41,65,74] and references therein). The interest for composite structures in the engineering community has constantly grown in recent years due to their appealing mechanical properties such as increased stiffness and strength, reduced weight, improved corrosion and wear resistance, just to recall some of them. The majority of man-made composite materials consists of reinforced fibers embedded in a base material, called matrix (see, e.g., [34,65]). The matrix material keeps the fibers together, acts as a load-transfer medium between fibers, process which takes place through shear stresses, and protects those elements from being exposed to the environment, while the resistence properties of composites are given by the fibers which are stiffer and stronger than the soft matrix. In this paper we focus on laminated composite materials, which are formed by a collection of building blocks or plies, stacked to achieve the desired stiffness and thickness.
For this kind of structures also simple loading conditions, such as traction or bending, cause a complex 3D stress state because of the difference in the material properties between the layers, which may lead to delamination, con-sequently requiring an accurate stress evaluation through the thickness (see, e.g., [65,70]). As an alternative to two-dimensional theories, often insufficiently accurate to depict delamination and interlaminar damage, and to layerwise theories, which typically show an high computational cost, a novel method combining an isogeometric analysis (IGA) Galerkin approach with a stress-recovery technique has been recently proposed in [23]. Introduced in 2005 by Hughes et al. [36], IGA aims at integrating design and analysis employing shape functions typically belonging to Computer Aided Design field (such as, B-Splines and NURBS). Using the same shape functions to approximate both geometry and field variables leads to a cost-saving simplification of expensive mesh generation and refinement processes required by standard finite element analyisis. One of the most important features of IGA is the high-regularity of its basis functions leading to superior approximation properties. IGA proved to be successful in a wide variety of problems ranging from solids and structures (see, e.g., [5,11,13,22,24,37,39,49,58]) to fluids (see, e.g., [1,12,30,50]), fluid-structure interaction (see, e.g., [10,35]), opening also the door to geometrically flexible discretizations of higher-order partial differential equations in primal form as in [3,29,43,51]. However, a well-known important issue of IGA is related to the development of efficient integration rules when higher-order approximations are employed (see, e.g., [7,26,38,67]). In attempt to address this problem taking full advantage of the special possibilities offered by IGA, isogeometric collocation (IGA-C) schemes have been proposed in [4]. The aim was to optimize the computational cost still relying on IGA geometrical flexibility and accuracy.
Collocation main idea, in contrast to Galerkin-type formulations, consists in the discretization of the governing partial differential equations in strong form, evaluated at suitable points. Since integration is not required, isogeometric collocation results in a very fast method providing superior performance in terms of accuracy-to-computational effort ratio with respect to Galerkin formulations, in particular when higher-order approximation degrees are adopted (see [68]). Isogeometric collocation has been particularly successful in the context of structural elements, where isogeometric collocation has proven to be particularly stable in the context of mixed methods. In particular, Bernoulli-Euler beam and Kirchhoff plate elements have been proposed [64], while mixed formulations both for Timoshenko initially-straight planar [17] and non-prismatic [9] beams as well as for curved spatial rods [8] have been introduced and studied, and then effectively extended to the geometrically nonlinear case [47,53,54,55,75,77]. Isogeometric collocation has been moreover successfully applied to the solution of Reissner-Mindlin plate problems in [44], and new formulations for shear-deformable beams [45,47], as well as shells [46,56] have been solved also via IGA collocation. Since its introduction, many promising significant works on isogeometric collocation methods have been published also in other fields, including phase-field modeling [31], contact [21,48,76] and poromechanics [59]. Moreover, combinations with different spline spaces, like hierarchical splines, generalized B-Splines, and T-Splines, have been successfully tested in [14,52,68], while alternative effective selection strategies for collocation points have been proposed in [2,31,57]. IGA-Galerkin methods have already been used to solve composite laminate problems, especially relying on high-order theories for enhanced plate and shell theories [27,42,66,69,72]. Recently an isogeometric collocation numerical formulation has been proposed [61] to study Reissner-Mindlin composite plates. Other Galerkin methods [32,33,66] compute instead a full 3D stress state using isogeometric analysis, applying a layerwise technique. This method can also be applied to isogeometric collocation adopting a multipatch approach, which models each layer as a patch (see Figure 1(a)), enforcing normal stress continuity at the inter-patches boundaries [6]. Clearly a layerwise method exploits a number of degrees of freedom directly proportional to the number of layers, inevitably leading to high computational costs. In this paper we apply a single patch 3D isogeometric collocation method to analyze the behavior of composite plates. We adopt a homogenized single-element approach (see Figure 1(b)), which conveniently uses one element through the thickness, coupled with a post-processing technique in order to recover a proper out-of-plane stress state. This method is significantly less expensive compared to a layerwise approach since employs a considerably lower number of degrees of freedom.  The post-processing approach, first proposed in [23], takes inspiration from recovery techniques which can be found in [16,19,25,42,63,73] and is based on the direct integration of the equilibrium equations to compute the out-of-plane stress components from the in-plane ones directly derived from a coarse displacement solution.
The structure of the paper is organized as follows. In Section 2 the fundamental concepts of multivariate B-Splines and NURBS are presented, followed by an introduction to isogeometric collocation and a description of our IGA-C scheme for orthotropic elasticity. In Section 3 we define our isogeometric collocation strategy to study laminated plates, which combines a homogenized single-element approach with an equilibrium-based stress recovery technique.
In Section 4 we present our reference test case and provide results for the single-element approach. Several numerical benchmarks are displayed, which show a significant improvement between non-treated and post-processed outof-plane stress components. Finally we provide some mesh sensitivity tests considering an increasing length-to-thickness ratio and numbers of layers to show the effectiveness of the method. We draw our conclusions in Section 5.

Isogeometric Collocation: Basics and application to orthotropic elasticity
In this section we introduce the notions of multivariate B-Splines and NURBS, provide some details regarding isogeometric collocation and describe our collocation scheme in the context of linear orthotropic elasticity.

Multivariate B-Splines and NURBS
In the following, we introduce the basic definitions and notations about multivariate B-Splines and NURBS. For further details, readers may refer to [15,36,62], and references therein. Multivariate B-Splines are generated through the tensor product of univariate B-Splines. We denote with d p the dimension of the parametric space and therefore d p univariate knot vectors have to be introduced as where p d represents the polynomial degree in the parametric direction d, and m d is the associated number of basis functions. Given the univariate basis functions N d i d ,p d associated to each parametric direction ξ d , the multivariate basis functions B i,p (ξ) are obtained as: where i = {i 1 , ..., i dp } plays the role of a multi-index which describes the con- where the coefficients P i ∈ R ds of the linear combination are the so-called control points (d s is the dimension of the physical space) and the summation is extended to all combinations of the multi-index i. NURBS geometries in R ds are instead obtained from a projective transformation of their B-Spline counterparts in R ds+1 . Defining w i as the collection of weights according to the multi-index i, multivariate NURBS basis functions are obtained as and NURBS multidimensional geometries are built as (see, e.g., [40]), but this represents just the simplest possible option (see, e.g., [18,20] for alternative choices). Along each parametric direction d, Greville abscissae consist of a set of m d points, obtained from the knot vector

An Introduction to Isogeometric collocation
p d being the degree of approximation. Since the approximation is performed through direct collocation of the differential equations, no integrals need to be computed and consequently, evaluation and assembly operations lead to a significantly reduced computational cost.

Numerical formulation for orthotropic elasticity
Once a strategy to select collocation points and compute IGA shape functions is set, a proper description of the equations in strong form for the problem under examination is required, as mentioned in Section 2.2. We therefore recall the classical elasticity problem in strong form considering a small strain regime and detail equilibrium equations using Einstein's notation (7). The where σ ij and u i represent respectively the Cauchy stress and displacement components, while b i and t i the volume and traction forces, n j the outward normal, and u i the prescribed displacements. The elasticity problem is finally completed by the kinematic relations in small strain as well as by the constitutive equations where C ijkm is the fourth order elasicity tensor.
As we described in Section 1, the basic building block of a laminate is a lamina, i.e., a flat arrangement of unidirectional fibers, considering the simplest case, embedded in a matrix. In order to increase the composite resistance properties cross-ply laminates can be employed (i.e., all the plies used to form the composite stacking sequence are piled alternating different fiber layers orientations) in which all unidirectional layers are individually orthotropic.
Since the proposed collocation approach uses one element through the thickness to model the composite plate as a homogenized single building block, we focus in this section on the collocation formulation for a plate formed by only one orthotropic elastic lamina. Considering three mutually orthogonal planes of material symmetry for each ply, the number of elastic coefficients of the fourth order elasticity tensor C ijkm is reduced to 9 in Voigt notation, that can be expressed in terms of engineering constants as The displacement field is then approximate as a linear combination of NURBS multivariate shape functions and control points as follows Having defined τ as the matrix of collocation points, we insert the approximations (11) into kinematics equations (8) and we combine the obtained expressions with the constitutive relations (9). Finally we substitute into where K ij (τ ) cofficients can be expressed as and substituting in (7b) we obtain: withK ij (τ ) components having the following form As we can see from equations (13), Neumann boundary conditions are directly imposed as strong equations at the collocation points belonging to the boundary surface (see, [6,21]), with the usual physical meaning of prescribed boundary traction.
3. An IGA collocation approach to model 3D composite plates In this section we describe our IGA 3D collocation strategy to model composite plates. The proposed method, known as single element approach, relies on a homogenization technique combined with a post-processing approach based on the imposition of equilibrium equations in strong form.

Single-element approach
The single-element approach considers the plate discretized by a single element through the thickness, which strongly reduces the number of degrees of freedom with respect to layerwise methods. The material matrix is therefore homogenized to account for the presence of the layers as Figure 1(b) clearly describes.
Remark 1. Considering a single-element homogenized approach is effective only for through-the-thickness symmetric layer distributions, as for nonsymmetric ply stacking sequences the plate middle plane is not balanced. In the case of non-symmetric layer distributions this technique is still applicable when the stacking sequence can be split into two symmetric piles, using one element per homogenized stack with a C 0 interface.
This method provides accurate results only in terms of displacements and in-plane stress components and, in order to recover a proper out-of-plane stress state, following [23], we propose to couple it with a post-processing technique. To characterize the variation of the material properties from layer to layer, we homogenize the constitutive behavior to create an equivalent single-layer laminate, referring to [71], where explicit expressions for the effective elastic constants of the equivalent laminate are given as ij represents the ij-th component of the fourth order elasticity tensor in Voigt notation (10) for the k-th layer and t k = t k h stands for the volume fraction of the k-th lamina, h being the total thickness and t k the k-th thickness.

Post-processing step: Reconstruction from Equilibrium
As interlaminar delamination and other fracture processes rely mostly on out-of-plane components, a proper through-the-thickness stress description is required. In order to recover a more accurate stress state, we perform a postprocessing step based on the equilibrium equations, following [23], relying on the higher regularity granted by IGA shape functions. This procedure, which takes its roots in [19,25,63,73], has already been proven to be successful for IGA-Galerkin. Inside the plate the stresses should satisfy the equilibrium equation (7a) that can be expanded as Assuming the in-plane stress components to well approximate the laminate behaviour, as it will be shown in Section 4, we can integrate equation 15a and 15b along the thickness, recovering the out-of-plane shear stresses as where ζ represents the coordinate along the thickness direction.
Finally we can insert equations (16a) and (16b) into (15c), recovering the σ 33 component as Following [23], the integral constants are chosen to fulfil the boundary conditions at the top or bottom surfaces.
Recalling that where the homogenized elasticity tensor C is constant, it is clear the necessity of a highly regular displacement solution in order to recover a proper stress state. Such a condition can be easily achieved using isogeometric collocation, due to the possibility to benefit from the high regularity of B-Splines or NURBS. We also remark that the proposed method strongly relies on the possibility to obtain an accurate description (with a relatively coarse mesh) of the in-plane stress state.

Numerical tests
In this section, to assess whether the proposed method can effectively reproduce composite plates behaviour, we consider a classical benchmark problem [60] and we address different aspects such as the effectiveness of the proposed post-processing step, the method sensitivity to parameters of interest (i.e., number of layers and length-to-thickness ratio), and its convergence.

Reference solution: the Pagano layered plate
A square laminated composite plate of total thickness t made of N orthotropic layers is considered. This structure is simply supported and a normal sinusoidal traction is applied on the upper surface, while the lower  Table 1 for 0°-oriented plies. The Neumann boundary conditions on the plate surfaces x 3 = ± t 2 are where σ 0 = 1 MPa.
The simple support edge conditions are taken as All results are then expressed in terms of the following normalized stress

Post-processed out-of-plane stresses
In this section, we comment the results obtained using the proposed IGAcollocation approach, as compared with Pagano's analytical solution [60].
To give an idea of the improvement granted by the post-processing of out-ofplane stress components, in Figures 3 and 4 we compare the reference solution with non-treated and post-processed results for the cases with 3 and 11 layers, considering a length to thickness ratio S = 20. All numerical simulations are carried out using an in-plane degree of approximation p = q = 6 and 10 collocation points for each in-plane parametric direction, while we use an approximation degree r = 4 and one element through the thickness (i.e., r + 1 collocation points). The sampling point where we show results is located at    To show the effect of post-processing at different locations of the plate, in Figures 5-7 the out-of-plane stress state profile is recovered sampling the laminae every quarter of length in both in-plane directions, for the case of a length-to-thickness ratio equal to 20 and 11 layers.  Figure 5: Through-the-thicknessσ 13 profiles for several in plane sampling points. L represents the total length of the plate, that for this case is L = 220 mm (being L = S t with t = 11 mm and S = 20), while the number of layers is 11 ( post-processed solution, × analytical solution [60]).
x 1 = 0  Figure 6: Through-the-thicknessσ 23 profiles for several in plane sampling points. L represents the total length of the plate, that for this case is L = 220 mm (being L = S t with t = 11 mm and S = 20), while the number of layers is 11 ( post-processed solution, × analytical solution [60]).

Convergence behaviour
In order to validate the proposed approach in a wider variety of cases, computations with a different ratio between the thickness of the plate and its length are performed respectively for 3, 11, and 33 layers, considering an increasing number of knot spans. Figures 8 and 9 assess the convergence behaviour of the method, adopting the following error definition Note that relation (22) is used only to estimate the error inside the domain to avoid indeterminate forms. Different combinations of degree of approximations have been also considered. A poorer out-of-plane stress approximation is obtained using a degree equal to 4 in every direction, and, in addition, with this choice locking phenomena may occur for increasing values of length-tothickness ratio. Therefore, we conclude that using a degree of approximation equal to 6 in-plane and equal to 4 through the thickness seems to be a reasonable choice to correctly reproduce the 3D stress state. Using instead uniform approximation degrees p = q = r = 6 does not seem to significantly improve the results (see Figures 8,9,and Table 2). The post-processing method provides better results for increasing values of length-to-thickness ratio and number of layers and therefore proves to be particularly convenient for very large and thin plates. This is clear since a laminae with a large number of thin layers resembles a plate with average properties. What really stands out looking at the displayed mesh sensitivity results, is the fact that collocation perfectly captures the plates behaviour not only using one element through the thickness but also employing only one knot span in the plane of the plate.
A single element of degrees p = q = 6 and r = 4, comprising a total of 7x7x5 collocation points, is able to provide for this example maximum percentage errors of 5% or lower (and of 1% or lower in the cases of 11 and 33 layers) for S = 30 or larger.  Quantitative results are presented in Table 2 for various plate cases, considering a number of layers equal to 11 and 10 collocation points for each in-plane parametric direction. Different number of layers (i.e., 3 and 33) are instead investigated in Appendix A. Increasing length-to-thickness ratios, namely 20, 30, 40, and 50 are considered and the maximum relative error results is reported for a reference point located at x 1 = x 2 = 0.25L. Also different degrees of approximation are investigated. Given these results, we conclude that using an out-of-plane degree of approximation equal to 4 leads to a sufficiently accurate stress state. Furthermore the out-of-plane stress profile reconstruction shows a remarkable improvement for increasing values of number of layers and slenderness parameter S.

Conclusions
In this paper we present a new approach to simulate laminated plates characterized by a symmetric distribution of plies. This technique combines a 3D collocation isogeometric analysis with a post-processing step procedure based on equilibrium equations. Since we adopt a single-element appoach, to take into account variation through the plate thickness of the material properties, we average the constitutive behaviour of each layer considering an homogeneized response. Following this simple approach, we showed that acceptable results can be obtained only in terms of displacements and in-plain stresses. Therefore, we propose to perform a post-processing step which requires the shape functions to be highly continuous. This continuity demand is fully granted by typical IGA shape functions. After the post-processing correction is applied, good results are recovered also in terms of out-of-plane