## 1 Introduction

Laminated composite materials are widely used in the aerospace and automotive as well as construction industries. Their lightness, strength and stiffness make them a very appealing solution in multiple situations, and those are also the reasons behind their continuous spread to many other consumer goods during the last decades.

Composite laminates are made of stacked material layers, also denoted as plies, that present different mechanical properties. These laminae are frequently made of stiff reinforcement fibers embedded in a soft matrix material that glues them together, as it is the case of glass or carbon fiber materials. The constituent plies are usually stacked in such a way that their main fiber direction follows different orientations. Classical stack-up sequences are the and layer setups, in which plies of the same material are stacked using 2 or 4 different in-plane fiber orientation angles, respectively (see, e.g., reddy_practical_1995).

Due to this heterogeneous stack of layers, laminated structures present a quite complex mechanical response when loaded. In order to correctly assess the possible failure modes it is quite important to characterize precisely the strain and stress distributions in the different material layers and at their interfaces. We refer the interested reader to the classical references gibson_principles_2016; jones_mechanics_2018; vinson_behavior_2004; reddy_mechanics_2004 for a deeper insight on the many mechanical aspects of composite materials.

Beyond academic cases for which closed form solutions exist (e.g., pagano_exact_1970; varadan_bending_1991), the complexity of the laminates elastic behavior calls for the use of numerical methods. And, without doubt, the finite element method (FEM) is nowadays, and has been during the last decades, the most popular tool for the analysis of their mechanical response. A review of the different finite element methods available for their study is out of the scope of this paper, but we refer the interested reader to the numerous reviews on the topic, e.g. carrera_theories_2002; carrera_theories_2003; reddy_theories_1994; carrera_survey_2008; zhang_recent_2009; s_qatu_review_2012; kreja_literature_2011; liew_overview_2019, and the many references therein.

From that vast literature it is easy to realize that a high percentage of the published finite element models correspond to 2D plate and shell models zhang_recent_2009; kreja_literature_2011; liew_overview_2019. The two-dimensional FEM models are primarily founded on two main composite theories: the equivalent single-layer kreja_literature_2011; kulkarni_review_2018 and layerwise reddy_mechanics_2004; carrera_multilayered_1999; ferreira_analysis_2005 approaches. As discussed in the review liew_overview_2019

, layerwise based FEM models, that present good compromise between accuracy and computational cost, combine an in-plane 2D FEM discretization with additional degrees of freedom for describing the displacements and/or transverse stresses of all the material layers. In these approaches the number of required degrees of freedom scales linearly with the number of plies.

Similarly, 3D elasticity models are able to accurately reproduce the laminate’s strain and stress profiles by representing each ply with one or more elements along its thickness. As in the case of layerwise approaches, the number of degrees of freedom scales linearly with the number of layers, but, additionally, the computational cost of the matrix assembly is very high compared to 2D models liew_overview_2019. In this paper we propose a new methodology that significantly reduces this assembly cost, becoming analogous to the complexity of a 2D problem and independent of the total number of material layers.

The introduction of the Isogeometric Analysis (IGA) concept by Hughes et al. in hughes_isogeometric_2005; cottrell_isogeometric_2009 constituted a great success. Based on the high continuity of spline basis functions and their superior approximation properties compared to classical finite elements methods evans_n-widths_2009; beirao_da_veiga_estimates_2011, IGA has proven to be a powerful tool for a broad spectrum of applications: from solid mechanics cottrell_isogeometric_2006; elguedj_$bar_2008; kiendl_isogeometric_2009 to fluid dynamics bazilevs_variational_2007; hsu_dynamic_2015 and other fields gomez_isogeometric_2008; buffa_isogeometric_2010. This success is also extensive to the analysis of laminated composites, where different authors proposed reduced dimension methods for plates, shells and beams bazilevs_3d_2010; kapoor_interlaminar_2013; nguyen-xuan_isogeometric_2013; thai_isogeometric_2015; remmers_isogeometric_2015; pavan_bending_2017; farzam_new_2018; ghafari_isogeometric_2019.

On the other hand, despite their high accuracy, 3D solid IGA models in which at least one element along the thickness is used for every material ply guo_layerwise_2014; guo_layerwise_2015, present a reduced interest due to their very high computational cost, as in the case of similar FEM approaches. Even if in practice IGA models are able to achieve the same level of accuracy as FEM models, but using much coarser in-plane discretizations, the assembly cost per degree of freedom is often a major bottleneck. To reduce this cost is a quite active research field hughes_efficient_2010; auricchio_simple_2012; schillinger_reduced_2014; antolin_efficient_2015; barton_optimal_2016; barton_gaussgalerkin_2017; calabro_fast_2017; mantzaflaris_low_2017; hiemstra_fast_2019.

In order to overcome this burden, in dufour_cost-effective_2018 we proposed the use of 3D isogeometric models with a single element through the full thickness. In a post-processing stage a very high-fidelity stress profile is recovered from the coarse simulation result. This approach has been recently extended in patton_fast_2019 to the case of collocation methods. In dufour_cost-effective_2018, despite the use of single element through the thickness, the stiffness matrix assembly cost is still high, due to the use of layerwise quadrature schemes along the lamina thickness (the number of quadrature points scales linearly with the number of layers).

In this work, inspired by the sum-factorization technique melenk_fully_2001; ainsworth_bernsteinbezier_2011; antolin_efficient_2015, we propose a new method for the assembly of 3D laminated composites for both finite element and isogeometric methods. Our approach reduces drastically their assembly cost and is very effective for laminated structures that present many plies with the same configuration repeated along the laminate stack-up sequence (e.g., for and layer setups). By decomposing the computation of the required 3D integrals into its in-plane and out-of-plane contributions, the complexity of our assembly procedure scales as the assembly cost of a 2D problem, multiplied by the number of different ply configurations (for instance, 2 configurations for and 4 for ), being independent of the total number of laminae.

The computed stiffness matrices with this new technique are identical, up to machine precision, to the ones obtained using a standard assembly procedure.

The assembly times obtained with the proposed method are even faster than the ones achieved with the standard method, using a simplified 3D model with a single element through the whole lamina’s thickness, and combined with a material homogenization approach (e.g., sun_three-dimensional_1988). In addition, the homogenization approach would lead to less accurate results.

The rest of this paper is structured as follows. Initially, Galerkin methods for the analysis of laminate composite structures are introduced in Section 2, using finite element and isogeometric discretizations. The proposed fast assembly method is presented in Section 3, as well as a discussion of its theoretical computational complexity. In Section 4 numerical experiments that consider different isogeometric discretizations and laminate configurations illustrate the performance of the proposed method confronted with classical assembly techniques. Conclusions are drawn in Section 5. Finally, in A an alternative implementation that overcomes the use of Voigt’s notation is detailed.

## 2 Galerkin methods for 3D laminated composites

In this section we begin by introducing the linear elasticity problem for 3D composite laminates in the context of finite element and isogeometric discretizations. We initially describe the continuous problem and its notation in Section 2.1, followed by the family of used discretizations in Section 2.2, that we particularize, in Section 2.3, to the case of multi-layered materials. Finally, the stiffness matrix’s computational complexity is discussed in Section 2.4.

### 2.1 Continuous elasticity problem

Let us start by setting the continuous three-dimensional elasticity problem in strong form, namely

(1) |

where is displacement of the elastic body at every point, i.e., the problem unknown, is the applied loading and

the small strain tensor, computed as

. For small strains and linear materials the stress tensor is computed as , where is the fourth order elasticity tensor. Without constituting any limitation, and for the sake of simplicity, only homogeneous Dirichlet boundary conditions are considered.The domain occupied by the elastic body is the image of the parametric domain mapped with , which we assume to be a bi-Lipschitz homeomorphism. Thus, a parametric point is mapped into its physical image as .

In a classical way, considering the functional space , the variational form associated to the strong problem (1) can be formulated as: find such that:

(2) |

where the bilinear for is:

(3) |

### 2.2 Discrete elasticity problem

In order to discretize the strong problem (2), we follow a standard Galerkin approach and choose a finite dimensional space , such that and:

(4a) | |||

(4b) |

where are the space basis functions, is the space dimension and is a generic partition of the parametric domain.

###### Assumption 1.

In this work, the basis functions are chosen such that

(5) |

with , and . I.e., the basis functions are created as the combination of in-plane functions and out-of-plane functions . Therefore, the space dimension is .

Finite element spaces composed of Lagrangian hexahedron and wedge elements carry basis functions that fall into the category defined in (5). However, this is not the case of tetrahedral finite element meshes. We refer the interested reader to the classical references hughes_finite_1987; zienkiewicz_finite_2013; brenner_mathematical_2008; ciarlet_finite_2002 for further details on the definition of finite element spaces.

Regarding isogeometric discretizations, the basis functions definition (5) allows to consider different in-plane discretizations combined with non-rational B-spline basis functions along the third direction. Thus, standard B-spline basis functions hughes_isogeometric_2005; cottrell_isogeometric_2009, including NURBS, can be used in-plane, but also non-tensor product schemes such as HR-splines giannelli_thb-splines:_2012; vuong_hierarchical_2011, LR-splines dokken_polynomial_2013; bressan_properties_2013 or T-splines beirao_da_veiga_analysis-suitable_2013; bazilevs_isogeometric_2010.

Thus, based on the above defined finite space , we discretize the trial and test functions of the problem (2) by means of

where are their control point coefficients and . Plugging them into (3) the bilinear form becomes

(6) |

that can be expressed as:

(7) |

where is the problem’s stiffness matrix, and is built in the same way as . Using Voigt’s notation, the matrices can be computed as (see, e.g., hughes_finite_1987; zienkiewicz_finite_2013)

(8) |

where the strain-displacement matrices are calculated as

(9) |

and , for , is the

-th component of the gradient vector

. Accordingly, the material matrix , Voigt’s representation of the , is computed as:(10) |

By pulling-back the gradient to the parametric domain through , where

(11) |

it is possible to compute the integral (8) in the parametric domain as:

(12) |

where

(13) |

and is the Voigt’s representation of the pulled-back tensor .

### 2.3 Geometric structure of laminated composites

Let us now consider the formation of stiffness matrices for 3D laminated composite structures, that present a multi-layered structure along the shell thickness, as represented in Figure 1.

Here we assume that, as shown in that figure, the third parametric coordinate corresponds to the thickness direction. Then, the parametric domain can be decomposed as into its in-plane and out-of-plane sub-domains.

Along the thickness direction, layers are considered, being the coordinates of the layer interfaces , such that , and , for . Thus, the parametric domain of a single layer is defined as

(14) |

with and . By an abuse of notation we construct the image of every layer in the physical domain as .

While three parametric coordinates were used for the parametric domain , we now split them into the in-plane coordinates for , and the out-of-plane coordinate for (cf. Figure 1).

###### Assumption 2.

We assume that the map gradient does not depend on the third parametric direction, but only on the in-plane ones, i.e. .

This is a reasonable assumption in the case of 3D shell-like structures, whose parametrization is frequently created as the extrusion of a 2D freeform manifold along a constant direction:

(15) |

where is the 2D manifold and the extrusion direction. It is easy to realize in (15) that the gradient does not depend on .

We consider, for every layer along the thickness, a different material. Thus, a different elasticity tangent tensor is associated to every layer:

(16) |

Then, for the case of laminated composites, the stiffness matrix (12) can be computed like:

(17) |

where is the Voigt’s representation of every tensor .

### 2.4 Computational cost of classical stiffness matrix assembly

The standard way of computing the integrals above is to calculate their contribution element by element, considering only the functions that have support in every element, and then assembling all of them together. For that purpose, a quadrature based on the multi-layered structure must be aapplied.

Let us now, for the sake of exposition’s clarity and without constituting any limitation, assume that the structure’s mesh presents only one element through the thickness, as sketched in Figure 1. This is the case of “solid shell” discretizations. Thus, if the basis functions used in the discretization have degree along the three parametric directions, we consider an in-plane quadrature rule with points per element, while, in order to integrate precisely along the thickness, at least quadrature points must be used for every material layer. Hence, the total number of quadrature points per in-plane element is . In addition, there are non-vanishing basis functions in each element, therefore, the assembly’s computational complexity of every in-plane element is (see antolin_efficient_2015 for a further discussion).

As it can be seen, the number of quadrature points, and therefore the computational cost, scales linearly with the number of layers . When a high number of layers is considered, the assembly of the stiffness matrix becomes very expensive in terms of floating point operations.

###### Remark 1.

Even if for exposition purposes we limited our discussion to the case of a single element through the thickness, the previous result can be extended in a straightforward manner to the case in which the discretization has more than one element along the third parametric direction, according to the Assumption 1. In this case, as long as

quadrature points are used along the whole laminate’s thickness, independently of the discretization along the out-of-plane direction, the above complexity’s estimation above is still valid.

###### Remark 2.

In the case of wedge elements both the number of in-plane basis functions per element and quadrature points scale as . While they present out-of-plane functions and quadrature points. Therefore, the previous discussion regarding the method’s complexity is still valid. The same applies to the complexity results of the proposed fast method in Section 3.1.

In the next section we propose an alternative strategy for assembling the stiffness matrix. This new methodology presents a lower computational cost and guarantees its exactness up to machine precision.

## 3 Fast assembly of stiffness matrices for 3D laminate composite structures

Let us first introduced the last assumption in which this work relies upon.

###### Assumption 3.

We assume that for every single layer, the material only presents in-plane variations, and not out of plane. I.e., the tensor can be expressed as:

(18) |

In addition to purely homogeneous material plies, the assumption above is fulfilled by numerous 3D laminates, as it is the case of some fiber reinforced plies (e.g., carbon fiber) or honeycomb panels, that can present an heterogenous (and anisotropic) in-plane behavior, but their mechanical properties can be assumed homogeneous along the ply’s thickness.

We now explicitly introduce the three components of the gradient vectors present in (11):

(19) |

Based on Assumptions 1 and 2, the basis function gradients can be split in their in-plane and out-of-plane components as:

(20) |

where

Then, can be expressed as

(21) |

with

(22a) | ||||

(22b) |

Thus, the strain-displacement matrix in (13) can be also split in its in-plane and out-of-plane contributions:

(23) |

where the in-plane matrices and are built by substituting the gradient in (13) with and , respectively. Hence, considering Assumptions 2 and 3 and plugging the splitting (23) of into the computation of the stiffness matrix (17), we obtain:

(24) |

Finally, gathering the in-plane and the out-of-plane terms, the previous expression can be reformulated as

(25) |

where

(26) |

and:

(27a) | ||||

(27b) | ||||

(27c) | ||||

(27d) |

Thus, the 3D integrals involved in the computation of the stiffness matrix (17) are decomposed in (25) as combinations of 2D (26) and 1D (27) integrals. It is worth noting that the stiffness matrices computed by applying (25) will be identical, up to machine precision, to the ones obtained with (17), no approximation is introduced.

The operators are computed as 2D integrals for every material layer, however, from layer to layer, only the material term changes. Thus, in the cases in which the same material is used with different orientations, the matrices are needed to be computed only once for each of them. E.g., in the classical cross-ply design in which the same fibered material is stacked with orientations, the operators will be computed only for two different layers, one corresponding to the orientation, and another one for . In the also common case of orientations, the matrices will be computed for four different cases only.

###### Remark 3.

As it is detailed in, e.g., (reddy_practical_1995, Section 3.3-4), the material matrix of an orthotropic material (satisfying the Assumption 3) with a given in-plane orientation angle can be split as:

(28) |

where the matrices only depend on the material properties, but not on the orientation . The same decomposition applies to .

Therefore, in the frequent case in which the same material is used for the different layers, but changing its orientation from layer to layer, a maximum of five different matrices , one for every operator , must be computed.

On the other hand, in A we propose an alternative way for computing the matrices without the use of Voigt’s notation.

### 3.1 Computational cost of fast stiffness matrix assembly

Studying Equations (26) and (27) it is simple to realize that the computational cost of is negligible compared to the cost of . In fact, can be even pre-computed analytically, without the use of a numerical quadrature. Thus, the complexity of the proposed assembly technique is bounded by the computational cost of the matrices .

The computation of is carried out by assembling its contribution for every in-plane element. As before, considering degree along all the parametric directions, non-zero in-plane basis functions and in-plane quadrature points must be used in the integration of every single in-plane element. Therefore, the computational complexity of every matrix , considering all the involved layers, is , where is the number of different operators considered in the laminated structure. E.g., for the stack design , and for . As stated in the Remark 3, for the common case in which the same anisotropic material is used, but stacked with different orientations, it holds .

Thus, the computational complexity of the proposed method is much smaller than the one of the standard integration procedure: against , with (cf. Section 2.4). Additionally, for a small value and a high number of layers, it holds , what makes the difference between both complexities even higher.

###### Remark 4.

As far as the computational cost of the stiffness matrix assembly is bounded by the calculation of the 2D operators , the assembly complexity is independent of the number of elements or functions used along the thickness. Therefore, it is possible to use any degree, number of elements or continuity for the discretization along thickness at the same computational cost for the assembly.

This can be useful, for instance, in the case of solid shells discretizations, where the use of more than one element through the thickness would alleviate possible locking phenomena.

Moreover, our approach makes the use of 3D models viable: as pointed out in liew_overview_2019, the assembly burden is one the main obstacles that prevents their use in the analysis of laminated composites with a large number of laminae. However, it is worth reminding that other operations, such as the solution of the linear system of equations, are still dominated by the number of degrees of freedom and function’s continuity (see, e.g., the discussion in collier_cost_2012; collier_cost_2013 for isogeometric discretizations).

###### Remark 5.

In the computation of , in Equation (26), it is possible to apply the sum-factorization technique proposed in antolin_efficient_2015, that would allow to reduce the computation complexity to . This improvement has not been detailed in this work for the sake of brevity. By directly applying the sum-factorization method as described in antolin_efficient_2015, without splitting the assembly in its in-plane/out-of-plane contributions, a theoretical computational complexity is expected.

In addition, considering the weighted-quadrature techniques proposed in calabro_fast_2017; hiemstra_fast_2019 it would be potentially possible to improve, even further, the dependency on of the computational complexity of in the case of highly continuous spline basis functions. Another possibility for reducing that cost would be to apply the low-rank approximations techniques, as proposed in mantzaflaris_low_2017 also in the context of isogeometric analysis.

## 4 Numerical experiment

In this section we aim at illustrating the performance of the proposed assembly method (compared to the standard one) by computing the stiffness matrix associated of the classical Pagano plate problem pagano_exact_1970 in the context of tensor-product isogeometric discretizations. This problem consists in a flat rectangular plate composed of a group of stacked layers of the same linear elastic orthotropic material distributed with different orientations. The material mechanical properties are: , , , and , and the material layers are stacked in two different configurations, and .

The same discretization degree is used along all the parametric directions, and different values of are considered. Regarding the mesh discretization, different number of elements along the in-plane directions are considered, while one element through the thickness is used. In all the cases, non-rational continuous B-splines are used.

All the results shown below were obtained with implementations of the standard assembly method and the fast one proposed in this work based on the isogeometric analysis library GeoPDEs version 3.0 vazquez_new_2016.

In Figures 2 and 3 the assembly times of the standard method and the fast one proposed in this work are compared for the configuration. Figure 2 shows the assembly times as a function of the number of elements along the in-plane direction, for different numbers of materials and degrees; and, on the other hand, the assembly times as a function of the number of layers are shown in Figure 3. Analogously, the dependence of the assembly times respect to the number of layers is shown in Figure 4 for the configuration and different discretizations.

As it can be seen in Figures 2 to 4, the fast assembly method (solid lines) outperforms the standard one (dashed lines) for all the degrees, meshes and number of layers considered. For high degree and large number of layers, the proposed method reduces the assembly time of the standard procedure by more than two orders of magnitude. It is also worth mentioning that the matrices computed using the standard and fast methods are equal, term by term, up to machine precision.

As shown in Figures 3 and 4, and explained in Section 3.1, the assembly time of the standard method scales linearly with the number of layers , whereas for the fast method is constant and independent of the number of layers. The latter being true as far as the number of different material configurations is . E.g., as it can be appreciated in Figure 3, for the laminate design (), the layer case () is faster than the layers case (): the assembly time initially scales with the number of layers. Nevertheless, this is no longer true for , as the condition holds.

In addition, it is worth mentioning that in Figure 3 the assembly time for is less than twice the time for . This effect is caused by the overhead of evaluating basis functions, allocating data structures, etc. For a single layer this time overhead is not negligible respect to the computing time for evaluating the matrices . This is even clearer for and , where the cost of is negligible respect to other costs.

On the other hand, as it can be observed in Figures 3 and 4, for small number of elements ( and meshes, black and red solid lines, respectively) the assembly time is not exactly constant: it is initially flat, but presents a slight increment as the number of layers grows. This is due to the fact that the constant cost of the computation of the matrices (for a small number of elements) dominates for a small number of layers, whereas, as the number of layers increases, the cost of combining those matrices for layers (Equation (25)), with , starts to be dominant.

Finally, we would like to remark that, as highlighted in Section 1, the proposed method, considering multiple laminae, is even faster than a simplified 3D model that uses a single layer with homogeneized material properties, assembled with the standard procedure. The assembly times of the simplified model are equivalent to the ones obtained for the one layer cases shown in the numerical experiments (i.e., dashed black lines in Figure 2 and first point of all dashed lines in Figure 3).

## 5 Conclusions

In this paper we present a new method for assembling stiffness matrices of 3D laminate composite structures, in a linear elasticity framework, using isogeometric and finite element discretizations. This method relies on three main assumptions: the basis functions can be multiplicatively split into its in-plane and out-of-plane components; the gradient of the geometric parametrization map is constant along the thickness direction (as it is the case of geometries built as extrusion of 2D freeforms); and the mechanical properties of each material layer can only vary in-plane, and not out-of-plane (this is the case of many materials, e.g., carbon fiber plies, that present a 2D in-plane reinforcement, being homogeneous along the ply thickness). Based on these assumptions, we develop an assembly method in which the 3D stiffness matrix is computed as a combination of 2D and 1D integrals and that results in a matrix identical, up to machine precision, to the one obtained using the standard assembly procedure. The theoretical assembly cost of the proposed method per in-plane element is , where is the discretization degree and the number of different material configuratioins used in the cross-ply laminated composite. is such that , and in general it holds , where is the number of layers in the structure. Thus, the cost of the proposed method is virtually independent of . On the other hand, the standard integration method, in which quadrature points along the thickness for every material layer are considered, presents a computational cost per in-plane element , that scales linearly with .

The performance of the proposed method is illustrated by means of the classical Pagano plate problem in the context of isogeometric discretizations, in which we confront the obtained assembly times against the ones obtained with the standard method. The results show how the proposed method outperforms the standard one for any considered discretization: for high values of and a difference higher than two orders of magnitude, in terms of assembly times between both methods, was observed. A further improvement of the proposed method’s cost, in the context of isogeometric discretizations, may be achieved by applying weighted-quadrature techniques calabro_fast_2017; hiemstra_fast_2019 for the computation of the involved 2D integrals.

## Acknowledgements

The author gratefully acknowledges the support of the European Research Council, through the ERC AdG n. 694515 - CHANGE grant, and also thanks Annalisa Buffa for her helpful insights and suggestions.

## Appendix A Fast assembly of stiffness matrices for laminated composites without using Voigt’s notation

Planas et al. introduced in planas_b_2012 an alternative methodology for the assembly of stiffness matrices in the context of solid mechanics, denoted as “ free”, that overcomes the use of Voigt’s notation. Based on that approach, we detail in this Appendix the construction of the operators defined in (26) avoiding the use of Voigt’s notation, i.e, without building the material matrix (10) and the strain-displacement operators (23).

Using the bracket operator defined in planas_b_2012, the stiffness matrix (8) can be computed as

(29) |

where the contraction is such that the vector is contracted with the second component of , and with the fourth one, i.e.:

(30) |

Thus, being the Cartesian orthonormal basis, the contraction for small strain linear isotropic materials can be computed as planas_b_2012:

(31) |

where and are the Lamé coefficients and is the identity tensor. For orthotropic materials the contraction reads:

(32) |

with , where and are the main orthonormal in-plane material directions, and , , , with , are the nine material coefficients (see, e.g., (schroder_simple_2002, Section 3.3) for further details).

Let us now write the map gradient as a function of the covariant basis :

(33) |

where the Assumption 2 was considered. In the same way, its inverse can be expressed as

(34) |

where is the contravariant basis, such that and if , for .

Pulling-back the computation of to the parametric domain, as in (12), we obtain:

(35) |

that can be rewritten as

(36) |

where is the pull-back of with for the second and fourth components.

On the other hand, can be split according to its Cartesian components as (see Assumption 1):