Application of the Finite Element Method in the Analysis of Composite Materials: A Review.

The use of composite materials in several sectors, such as aeronautics and automotive, has been gaining distinction in recent years. However, due to their high costs, as well as unique characteristics, consequences of their heterogeneity, they present challenging gaps to be studied. As a result, the finite element method has been used as a way to analyze composite materials subjected to the most distinctive situations. Therefore, this work aims to approach the modeling of composite materials, focusing on material properties, failure criteria, types of elements and main application sectors. From the modeling point of view, different levels of modeling-micro, meso and macro, are presented. Regarding properties, different mechanical characteristics, theories and constitutive relationships involved to model these materials are presented. The text also discusses the types of elements most commonly used to simulate composites, which are solids, peel, plate and cohesive, as well as the various failure criteria developed and used for the simulation of these materials. In addition, the present article lists the main industrial sectors in which composite material simulation is used, and their gains from it, including aeronautics, aerospace, automotive, naval, energy, civil, sports, manufacturing and even electronics.

1. Introduction

Technological advancements have led to an increase in the demand of special materials with unique properties that cannot be found in metal alloys, ceramics or polymers blends [1,2].

To supply these needs, composite materials were developed. They are made of two or more distinctive and immiscible materials with different mechanical, physical and/or chemical properties [1,3,4,5].

Composites are considered heterogeneous and multiphase engineered materials, in which the matrix is responsible for binding the reinforcement together and transferring the loads between the fibers, while the reinforcement adds rigidity and obstructs crack propagation in the structure [1,6,7,8,9].

They can be classified according to the matrix (metallic, polymeric and ceramic) or the type of reinforcement used (fibers or particles) [1,7,8,10]. The ones with a polymeric matrix and continuous fibers have great relevance and significance, due to their excellent mechanical properties, good thermal stability and low density [11].

Fiber reinforcements can be either unidirectional (UD) (Figure 1a) or bidirectional (Figure 1b). Unidirectional fibers have maximum stiffness and strength along the fiber direction and minimal properties in the transverse direction, exhibiting anisotropy. Bidirectional reinforcements have maximum stiffness and strength in the fiber direction [12,13]. Unidirectional fibers can be aligned on a thin plate; pre-impregnated with resin; and used to define the stacking order and layer design for composite laminates (Figure 1c). The bidirectional reinforcements are woven fabrics, and there are several types of weaving [12,13].

Figure 1

Scheme of the most common uses of fiber reinforced composite structures: (a) unidirectional fiber orientation ply, (b) bidirectional fiber orientation ply (woven-ply) and (c) multiorientation laminate, quasi-isotropic laying-up sequence [0°/45°/90°/−45°] 6S. Reproduced with permission [13].

Anisotropic materials show different mechanical properties in each direction; i.e., they are not symmetrical with respect to all their planes or axes. Orthotropic materials are a subset of anisotropic materials that show a symmetry between two planes, in general, the plane parallel to the fibers has significantly superior properties compared to the orthotropic perpendicular plane. Wood is a good example of an orthotropic material; its properties perpendicular to the fiber axis (radial and tangential) are worse than its parallel ones. In this case, the properties of radial and tangential directions are not equal, but similar, being equally inferior to the longitudinal direction.

Composites have a set of performance characteristics that their constituents cannot achieve by themselves individually [5,7]. Due to these combinations, it is possible to obtain lightweight design with high strength and stiffness; some other key characteristics are high-temperature, corrosion and impact resistance. Together, said things make composites more interesting, useful and attractive alternatives [14,15,16,17].

Because of these characteristics, composites are widely applied in automotive, aeronautical, petrochemical, naval, electro-electronic, civil construction, energy, biomedical, sports and manufacturing industries, among others [1,5,14,15,18,19,20]. Their applications can be seen in several industrial sectors; however, they are expensive and difficult to characterize due to their heterogeneity and laminate configurations, which affect their final properties.

Owing to this difficulty and their toughness, to optimize and improve structural projects, as well as to understand better the behaviors of these materials, some researchers resort to computational simulations. Through the use of the finite element method (FEM), it is possible to understand the damages caused in the matrix, the fiber and their interface when the composite is exposed to severe conditions, such as static and dynamic loading, different temperatures and pressures, among others [21,22].

Analytical models are not always able to sufficiently address all failure phenomena that contribute to composite performance [23]. Different failure mechanisms play important roles during the service-lives of composite materials; for example, fracturing of the reinforcement is a partial detachment of the interface, which results in the nucleation and growth of voids, and their coalescence in the matrix.

According to Lasri et al. [24], damage mechanisms in composite materials generally include four types of failure modes: transverse matrix fracture, fiber–matrix interface detachment, fiber rupture and layer delamination. In general, transverse fracture of the matrix is the first damage process to occur, since the matrix has lower failure stress compared to all laminate constituents. Fiber–matrix interface detachment can accompany transverse fracture of the matrix and facilitate its progression [24,25]. Transverse failure can happen without breaking any longitudinal fiber. Such failures are parallel to the fibers and lead to a decrease in stiffness. Thus, damage criteria are required to indicate the onset of failure and damage orientation.

Note that when applying FEM to composites, it is important to consider some particularities of these materials, such as the constitutive law; modelling and failure criteria associated with the composites; and the type of elements used to model the objects.

2. Modelling

For the simulation of composites, three primary approaches are usually applied, those being: (a) a micromechanics-based approach, (b) an equivalent homogeneous material (EHM) based approach and (c) a combination of the two previous approaches. Notice that each method has advantages and disadvantages [23,26,27].

According to Dandekar and Shin [23] the micromechanical based approach describes the material behavior locally, and thereby, it is possible to study local defects, such as fiber–matrix detachment and complex deformation mechanisms, especially in fiber reinforced composites. However, the time required to solve a simulation is very high, because the mesh used is very fine compared to the EHM model.

The EHM approach reduces simulation time, but it is not able to predict local effects; e.g., damage at the fiber–matrix interface [23,28,29]. Dandekar and Shin [23] said that it is possible to take advantage of the two models, and their ability to predict shear force and sub-surface damage.

Venu Gopala Rao, Mahajan and Bhatnagar [30] applied this methodology in their work. They studied the machining of the UD-GFRP composite’s chip formation mechanism. For that, the portion of the work material adjacent to the cutting tool was modelled—fiber and matrix separately—whereas portions away from the cutting tool were modelled as EHM (Figure 2).

Figure 2

Schematic view of the fiber, matrix and equivalent homogeneous material (EHM) domains used in a finite element model for the case of 45° fiber orientation. Reproduced with permission [30].

For them, the most important thing was to study the cutting zone and that is why they modelled this region with the specification of a composite. As for the tool and the region around the cutting zone, single results were not so important to them, those being analyzed with the EHM model.

Jones [31] noted that these approaches define how rigidity and strength properties are chosen for the project materials. A composite material’s behavior can be separated into the micro and the macroscale, defined as follows:

Microscale—study the composite material’s behavior, for which interactions from constituent materials are examined in detail as defined by heterogeneous material behavior.

Macroscale—study the composite material behavior considered to be homogeneous, and the effects of all constituent materials are detected only by the composite material’s mean apparent properties.

Tenek and Argyris [22] went further in their conjectures; they cited that these questions address two fundamental problems: how to define the sheet properties using microscale procedures, and how to apply these properties on a macroscale for a global analysis.

At the microscale there are many difficulties due to the uncertainties that may require stochastic or statistical models. The objective of most approaches is to define the composite modules from all constituent materials or the strengths of the composite in terms of its phases. Therefore, some basic approaches include using the materials’ mechanisms and elasticity theories based on the repetition of a unit cell or some other representative volume, assuming that there is a perfect bond between fibers and matrix, which may not be true most of the time. Frequently, micro-mechanical theories are validated with experimental work [22,31,32].

In Figure 3 is shown a structural evolution scheme for the hierarchy of a unit cell starting from a ply to a global composite structure.

Figure 3

From the microscale to the macroscale. Reproduced with permission [33].

In the simulation study of composites, it is possible to evaluate properties from nanoscale to macroscale, or in other words, to apply the multiscale technique, which consists of simulating the behavior of a composite through multiple time and/or length scales [34,35,36]. Some applications of these techniques are mainly focused on microstructural and mechanical property simulations of many classes of materials, including nanocomposites [35,37,38,39,40].

Usually, both micro and macroscale approaches are extensively applied, but with the increase in nanomaterials as polymer composite reinforcements, nanoscale has been gaining ground. Coupled with the need for a better understanding of damage progression, mesoscale is thus an optimal option (Figure 4).

Figure 4

Schematic view of a four-scale woven fiber composite with polymer matrix: In computational modeling of this structure, each integration point at any scale is a realization of a structure at a finer scale. Due to the delicacy of materials at fine-scales, RVEs at lower scales may embody more uncertainty than those at higher scales. To quantify the uncertainty in a macroscopic quantity of interest, the relevant uncertainty sources at the lower scales should be identified for uncertainty propagation. Reproduced with permission [41].

Multiscale techniques improve classical method solutions and analysis; i.e., solving local problems and taking this information from the smallest scale to the highest level, in a hierarchical process, through homogenization technique [42]. Aboudi, Arnold and Bednarcyk [34] explain that the homogenization technique provides the properties or responses of a "structure" (higher level) when given the properties or responses of the structure "constituents" (smaller scale). In summary, information from lower scale levels is used by higher scale levels to investigate the effects of constituents and microstructure on the mechanical properties, and those of the structural performances of the parts on macroscale properties [43,44].

Shaik and Salvi [45] commented that in studies of large structures, which have millions of components; materials and nonlinear structures; and many joints, it is impossible to model all of these details due to computational limitations. In such a case, a multiscale modelling approach is used when only a particular area of interest needs to be modeled with precision and in detail, saving a lot of time and computing power. However, the authors emphasized that these methodologies require experience and careful judgment by the designer.

In this approach, a geometry, with a repeating unit cell (RUC) and all of it constituents, is developed to build a representative volume element (RVE), along with its components and their interactions, to be modeled. The structure is shaped based on RUC/RVE and analyzed on different length scales with the desired confidence level, adding damage and failure aspects. The results undergo qualitative and quantitative evaluations from the material, configurational and architectural perspectives [34,45].

Shaik and Salvi [45] explained that the project can be divided into several scale sizes, depending on the confidence of the associated scale theory and the level of interest; i.e., the scale could be divided into micro (RUC), meso (RVE) and macroscale. At the micro level, the study focuses on the fiber: composition, geometry and orientation within the RUC. At meso level, mechanical characteristics of the material built from many of RUCs are studied, which produce homogeneous properties independent of any final effect and influence from the properties of the components. Finally, macro level includes mesoscale properties when the laws of continuous mechanics are applied. By connecting these scales, the performances of macroscale structures can be related to microscale individual constituents, such as the fiber, the resin and their interface.

Figure 5 presents a schematic of multiscale modeling used to analyze the blade of a wind turbine, showing the use of the microscale, mesoscale and macroscale, and the example of the application of RCU and RVE of a unidirectional laminate.

Figure 5

Schematic of multiscale modelling of engineering composite structures. Reproduced with permission [46].

Kwon, Allen and Talreja [47] created flowcharts (Figure 6 and Figure 7) exemplifying what happens at each of the scales, where the microscale basically approaches the individual fiber and matrix–fiber interfaces, the mesoscale approaches the individual layers and the macroscale considers effect of the complete laminate homogeneously [48].

Figure 6

Hierarchy of multiscale analysis for a unidirectional fiber reinforced composite.

Figure 7

Multiscale analysis hierarchy for a two-way fiber reinforced composite.

Initially, Kwon, Allen and Talreja [47] developed the flowchart for unidirectional composites (UD) (Figure 6). First, in the "Fiber Module," the properties of the fiber and matrix materials, and the geometric properties, are correlated to defining the composite properties. These properties are used for each blade with its fiber orientation relative to the global coordinate system [47].

The "Laminating Module" calculates the properties of the laminated composite, so these properties are used for finite element analysis of the structure, completing the stiffness loop. The inverse order is used to decompose the stress and deformations from the macro level to the micro-level; i.e., stress and deformations in the fiber and matrix materials.

After the calculation of microscale stresses and strains, the damage and/or failure criteria are applied. Because damage and failure are described at the constituent level, damage and failure modes are simplified and based on physics. At the microscale, there are three possible damages and failures: fiber rupture, matrix failure and interface detachment. Different failure or damage criteria can be applied to these three damage modes. At the macro level, there are more complex damage modes, such as delamination. Only the location and orientation of the damage or failure will dictate the difference between the macroscale failure modes. As a result, the damage and failure modes can be understood in unified and simplified concepts [47].

For the woven fibers (MD), the "Fabric Module" is added which relates the properties of the UD fiber to the effective properties of the fabric (Figure 7). The purpose of this module is to calculate the properties of the fabric using information from the UD fibers and woven fabric and the decomposition of the deformations and tensions in the fibers [47].

Mao et al. [43] demonstrates in his work how to apply this methodology. The authors started the modelling from the micromechanical computation (Figure 8a), for which a unit cell with a square arrangement was used to represent the behavior of the UD material of the fiber yarn.

Figure 8

(a) Square arrangement of the microscale unit cell and (b) mesoscale, macroscopic unit cell.

It is assumed that the fibers are arranged in a uniform distribution with the measured volumetric fraction and the same average filament diameter. In addition, the mesoscale unit cell (Figure 8b) is modelled to describe the woven architecture of the fiber yarns and the composite resin pocket. Two types of fiber yarn are modelled: weft yarn (longitudinal direction) and warp yarn (transverse direction). In FE analysis, periodic boundary conditions are used to eliminate boundary effects.

As shown in Figure 9, the strained material properties of the macroscale model are obtained from micromechanical and mesomechanical analyses.

Figure 9

Multiscale modelling strategy for woven composite laminates.

According Tian et al. [49], mechanical behaviors of heterogeneous materials are often described by using RVEs in the FE. The author mentioned two theories for RVE: first, Hill’s theory, for which the RVE must be large enough to contain a large number of fibers in the heterogeneous materials and be a statistical representation of the heterogeneous materials. The effective properties derived from the RVEs represent the real properties of the material on the macroscopic scale, which is commonly known as the micro-meso-macro principle, since scale separation is necessary [50]. Alternatively, in Drugan and Willis’s theory, the RVE must contain the smallest volume of composites for which the mean mechanical responses remain constitutively valid [51].

Tian et al. [49] pointed out that FE with RVEs is not common for modeling composites reinforced by fibers randomly distributed on a microscale, because these micro-architectures are much more difficult to model sometimes due to the high-volume fractions and large fiber aspect ratios to consider [52,53]. Therefore, in order to numerically model composites reinforced by spatially randomly discontinued fibers, it is important to generate RVEs with high fiber volume fractions and large fiber proportions.

The literature mentions two approaches as the most usual ones, namely, the random sequential adsorption (RSA) algorithm and the Monte Carlo (MC) procedure, for generating the artificial RVEs with randomly distributed fibers [49,51,52,53]. However, according to Lu, Yuan and Liu [50], with these approaches it is difficult to generate RVEs with high fiber aspect ratios (FARs) and fiber volume fractions (FVFs). A new approach to try to circumvent this situation is the use of the automatic searching and coupling (ASC) technique, in which it is possible to generate the 3D RVE to analyze the composite with random fibers with a wide range of FARs [50,51,54].

Lu, Yuan and Liu [50] say: “compared with the conventional model, the present model is easier to generate and more time-saving as it eliminates the drawback of free meshing. In addition, the ASC technique can remove the additional stiffness introduced by the embedded element technique, and hence can improve precision and convergence. Moreover, our technique facilitates the direct application of the 3D periodic boundary conditions to the RVE.” See Figure 10.

Figure 10

The sketch for the fiber architectures in the 3D model, (a) overall spatial view, (b) top vie, (c) side view, and (d) a fiber is described by the center point C, and two Euler angles θ and φ. Reproduced with permission [50].

Cohesive Zone Model

Another method that has been widely applied in composite materials is the cohesive zone model (CZM). According to Barbero [27] the CZM is based on the assumption that the stress transfer capacity between the two separating faces of delamination is not lost completely at damage initiation, but rather is a progressive event governed by progressive stiffness reduction of the interface between the two separating faces (Figure 11).

Figure 11

Cohesive zone model to simulate crack propagation.

The models are typically expressed as a function of normal and tangential tractions in terms of separation distances. The forms of the functions and parameters change from model to model [55,56,57,58].

CZM has been used previously to study crack tip plasticity and creep under static and fatigue loading conditions, and in polymer cracks, adhesive joints, interface cracks in biomaterials and crack bridging due to fibers and ductile particles in composites [55,58,59]. Currently, CZMs are increasingly being used to simulate discrete fracture processes in various systems of homogeneous and non-homogeneous materials [55,56,57,60], with great emphasis on understanding the evolution of delamination in laminates [55,59,61,62,63,64,65].

According to Roy [58], in the presence of a large fracture process zone near the crack tip, the basic assumptions of the mechanics of linear elastic fracture are no longer valid. Specifically, in some polymers, such as hardened epoxies, the occurrence of void nucleation and growth ahead of the crack tip results in a damage zone that is not free from traction. Additionally, for a crack in a composite with a fiber-reinforced polymer matrix, the fiber bridge may be present within the damage zone. Therefore, in these cases, a cohesive layer modelling approach would be more accurate at accounting for non-linear processes that occur within the "damage zone" [58]. Furthermore, the interface modelling using CZM has a distinct advantage compared to other global approaches (e.g., shear lag model), in that it is based on a micromechanical method [55].

The cohesive zone model was proposed by Barenblatt in the 1960s based on Griffith’s theory of fracturing, in order to investigate the crack propagation in brittle materials. He assumed that finite molecular cohesion forces exist near the crack faces and described the crack propagation in perfectly brittle materials using his model. Then, Dugdale extended this concept to the perfectly plastic materials by postulating the existence of a process zone at the crack tip [55,57,58,66,67]. Roy [58] mentions that in the coming years, many uses of CZMs were used to better understand the functioning of cracking in laminated composites, with the objective of capturing the Burridge-Andrew mechanism using the material point method.

Using the cohesive modelling, no additional properties are necessary to simulate crack growth. Only the cohesive law is needed to analyze both the initiation and growth of a crack. Typically, cohesive elements in FEM codes follow a predefined traction separation law that simulates the crack initiation and propagation. Another advantage of CZMs is that these models can simulate different types of failure mechanisms, such as fiber-matrix debonding and interlaminar delamination [58,61,62,65,68,69].

Chaboche [57] mentions that the cohesive model considers the presence of a process zone at the tip of the crack, with an appropriate constitutive law, relating the normal tensile stress (T) and the relative displacement (u) between the two sides of the crack. The relation between T and u is characterized by a softening law (decreasing function) and the area under the stress-displacement response corresponds to the fracture energy Gc of the material.

According to Khoramishad et al. [67] the cohesive zone model (Figure 12), combines a strength-based failure criterion to predict the damage initiation and a fracture mechanics-based criterion to determine the damage propagation.

Figure 12

Schematic damage process zone and corresponding bi-linear traction–separation law in an adhesively bonded joint. Reproduced with permission [67].

3. Constitutive Laws of a Composite Material

In general, the constitutive equation of a linear elastic solid is known as the generalized Hooke’s Law, which relates nine Cauchy stress components with nine deformation components, giving a total of eighty-one constants [31,70,71,72].

Equations (1) and (2) present the combination of elasticity constants, where represents the stress component, the strain components, the stiffness matrix and the compliance matrix of the material, those being inversely related as follows: .

According to Bednarcyk, Aboudi and Arnol [73], and Bauchau and Craig [74] the physical phenomena such as heat conduction, diffusion, electric permittivity, magnetic permeability and electric conductivity are governed by the material constitutive laws. These constitutive laws characterize the mechanical behavior of a material and consist of a set of mathematical idealizations of the observed behavior [74].

3.1. Anisotropic Material

Due to the growing importance of composite materials, the linearly elastic behavior of anisotropic materials must be understood. The physical properties of anisotropic materials are directional; i.e., the physical response of the material depends on the direction in which it acts [74].

According to Soriano [75] and Daniel and Ishai [72], a linear anisotropic material has a matrix with independent elastic properties, which, in turn, make the characterization very difficult. In general, the stiffness matrix has 36 independent coefficients (Equation (3)), owed to the symmetry between and , and between and ; hence the reduction from 81 to 36 elastic constants [3,31,70].

However, the symmetry requirement for anisotropic materials is reducing the elastic components to 21, by the relation [3,31,70,74]. Azevedo [71] pointed out that elastic properties required to define an anisotropic material can be represented by stress–strain ratios, the main coefficients being: longitudinal and transverse moduli of elasticity and the Poisson’s coefficient. Equation (4) shows the generalised mathematical representation of Hooke’s law for anisotropic materials, where E is the longitudinal elasticity modulus or Young’s modulus, G is the transverse modulus of elasticity or shear modulus, υ is the Poisson’s coefficient and is the angular/linear deformation.

For Vanalli [76] the failure analysis of structures made out of anisotropic materials is complex. The author emphasized that in these cases it must be assumed that the failure is caused by normal and shear stresses, since failure can occur due to different sets of stress acting on the element.

3.2. Orthotropic Material

Comparing an orthotropic material with a generally isotropic theory, it is observed that the first one presents three symmetry planes orthogonal to each other—x1 x2, x1 x3 e x2 x3, thereby reducing from 21 to 9 the number of independent coefficients [3,31,70,71,72,77]. The reasons for that are:

The angular deformations are independent of normal stress;

Linear deformations are independent of tangential stresses;

Each tangential tension causes only angular deformation in the plane in which it acts.

Equation (5) shows the generalised representation by Hooke’s law for orthotropic materials, where E is the modulus of longitudinal elasticity or Young’s modulus, G is the cross modulus of elasticity or shear modulus and υ is the Poisson coefficient [71].

When considering a composite reinforced with unidirectional fibers, it would be automatically classified as an anisotropic material because there is no symmetry between the planes. If we analyse the thickness, we can see that it is much smaller in size than any other plane dimension. Because of this, many researchers consider unidirectional laminates as orthotropic materials, taking into account the plane stress state, according to the hypotheses shown in Equation (6) [70].

3.3. Transverse Isotropic Material

A transverse isotropic material can be defined as an orthotropic material that presents isotropy in one of the planes of symmetry, which means it has the same properties in all directions in this plane [31,71,72,74].

Comparing a transverse isotropic material with a generally orthotropic material, it is observed that in transverse isotropic materials there is a symmetry between the planes x1 x3 and x1 x2, reducing nine to five independent coefficients. According to Azevedo [71] beyond the considerations for orthotropic materials, it can be added:

The linear deformations in the plane x2 x3 caused by the normal stress σ11 are equal;

The linear deformations ε22 and ε33 caused by the normal stress σ22 are equal to the deformations ε33 and ε 22, respectively, caused by a tension σ22 = σ33;

Each tangential tension only causes angular deformation in the plane in which it acts;

The angular strain γ23 caused by a stress σ 23 is equal to an angular strain γ13 caused by stress σ13 = σ23.

Equation (7) presents the generalised representation by Hooke’s law for transverse isotropic material, where E is the modulus of longitudinal elasticity or Young’s modulus, G is the transverse modulus of elasticity or shear modulus and υ is the Poisson coefficient [71].

4. Failure Criteria

According to Kaw [3] the success when using a composite structure is related to its efficiency and safety. For this, some criteria were adopted to identify possible failures associated with a component. For Ochoa and Reddy [78] and Jones [31] and Kaw [3], a failure criterion aims to provide a comprehension of the effects caused by combined loads (double or triple stress state) in the structure, indicating when there is a local or global failure. For Kaw [3], in general, the theories are related to normal and shear forces of the laminate, defining the stress states in which the failure occurs.

According to Jones [31] and Kaw [3], failure criteria were initially created for isotropic materials, where maximum normal and shear stresses of the material were found when the maximum stress was greater than the last force, indicative of material failure.

Among many failure criteria relevant to isotropic materials are the maximum normal stress (Rankine), maximum shear stress (Tresca), maximum normal strain (Saint–Venant) and maximum strain theory (Von Mises) [3,31,72,79]. Feng [80] and Kaw [3] cite that based on these theories, the failure criteria for anisotropic and orthotropic materials were developed; based on the orientation of the fibers, four parameters of normal resistance and one of shear resistance are considered, making a total of five fundamental resistance parameters for the use of the failure criterion (Figure 13).

Figure 13

Basic strength parameters of unidirectional lamina for in-plane loading, (a) longitudinal tensile, (b) longitudinal compressive, (c) transverse tensile, (d) transverse compressive, and (e) in-plane or interlaminar shear.

The failure theory problem of an orthotropic sheet to a certain extent is identical to the isotropic one, which in this case is the prediction of when a sheet is submitted to a biaxial or triaxial stress state, using resistance data obtained from uniaxial experiments [81].

Currently, there are several failure criteria for composites: Hill, Tsai–Hill, Tsai–Wu, Hashin-Rotem, Hashin, maximum stress, Hoffman, maximum strain, Hou, Puck–Schürmann, Chang-Chang, Linde, LaRC03, LaRC04, Maimí, Hart-Smith, Yeh-Stratton and others [16,82,83,84,85,86,87,88].

According to Mendonça [81], some of these criteria are of common use and well-established in the literature. Normally, these criteria are characterized by ignoring any aspect from any physical process involved in the failure, considering only macroscopic effects observed in the standard specimen. The most widely-used failure criteria are the maximum stress, Tsai–Hill, Tsai–Wu, Hashin and Puck–Schürmann (Table 1).

Table 1

Most common failure criteria for composite materials.

Failure Criteria	Formula	References
Maximum Stress	−X_c<σ_1<X_t; −Y_c<σ_2<Y_t; |τ_12|<S	[32,65,72,83,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109]
Tsai–Hill	σ12X^2+σ22Y^2−σ_1 σ_2X^2+τ122S^2≤1	[26,28,82,85,87,105,106,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124]
Tsai–Wu	F_1σ_1+F_2σ_2+F_6τ_12+F_11σ12+F_22σ22+F_66τ122+2F_12σ_1σ_2+F_16σ_1τ_12+F_26σ_2τ_12≤1	[82,83,85,87,106,111,122,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147]
Hashin	(σ_1X_T)^2+(τ_12S_12)^2=1; (σ_2Y_T)^2+(τ_12S_12)^2=1; σ_1=−X_c; (σ_22S_23)^2+[(Y_C2S_23)^2−1]σ_2Y_C+(τ_12S_12)^2=1	[24,83,85,87,93,104,105,107,122,127,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163]
Puck–Schürmann	∅_MC=(t_LS_L−η_Lt_N)^2+(t_TS_T−η_Tt_N)^2	[84,106,155,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178]

The maximum stress criterion, based on Rankine’s theory, is not an interactive criterion; i.e., it does not consider the combined effects of the various components of the tensor. This criterion provides for the rupture when one of the tensor components arrives at the corresponding tensile stress [179,180].

Hill, in 1950, based on the von Mises criterion, established one of the first failure criteria for anisotropic materials, which is a generalisation of the flow behavior for isotropic materials [181]. Although it is a more general criterion, it has as a drawback: the determination of several parameters to establish the complete equation of the model. In 1965, Tsai proposed a modified Hill criterion where he quantified the traction and compression inequality for orthotropic materials, which was called the Tsai–Hill criterion [181]. For several authors, the Tsai–Hill criterion is one of the best and most widely used failure criteria for laminates because it considers the interactions between the stress components; however, this criterion is not invariant in relation to the coordinate system; therefore, only orthotropic materials should be applied [3,31,81,179,182,183]. Despite being considered one of the best failure criteria, many consider that this criterion has several deficiencies in its theoretical basis [81,179]. Mendonça (2005) cites that there are basically three deficiencies in Tsai–Hill’s theory; namely:

It does not intrinsically consider differences in tensile and compressive strength;

It does not present good results in the state of loading by compression in the three main axes;

It supposes that a hydrostatic state of stresses cannot cause failure—in the case of anisotropic materials, a hydrostatic state of stress causes shear deformation and failure.

Tsai and Wu in 1971 presented a criterion based on the Tsai–Hill criterion, aiming to increase the number of terms in the Tsai–Hill failure criterion equation, to better approximate the experimental data, considering a two-dimensional stress state [81,181,182]. The Tsai–Wu criterion is an interactive criterion, which provides for component rupture due to the combination of tensions acting on the part [179,180,183]. In addition, this criterion, in its three-dimensional form, takes into account the effect of the hydrostatic component of the stresses differently from the previously-described criteria [184]. The interactions between the stress components are independent of the material properties. However, since it is not a failure criterion based on physical phenomena, it can predict the occurrence of the damage, but cannot distinguish between the different failure modes; it can only predict whether or not the failure occurs in the structure [81,180,183,184]. The Tsai–Wu criterion became one of the most used criteria, and to this day several works are developed based on the same. [181].

Hashin, in 1980, proposed a failure criterion divided into subcriteria, for failure in unidirectional fiber reinforced sheets—transversely isotropic—based on the quadratic polynomial of tensions [81,149,150,180,181,182,185]. Differently from the Tsai–Hill and Tsai–Wu criteria, which do not allow an identification of failure modes; the Hashin criterion considers modes of failure of the fiber and matrix, distinguishing between tensile and compression loads, addressing four main modes: traction and compression of the fibers and matrix [81,183,185]. According to Laurin, Carrere and Maire [186], the historical importance of this criterion is that it started a different way of designing failure criteria for composite materials. Hashin [150] first identified the predominant failure modes, and subsequently the variables associated with these modes, and then proposed the interactions between the variables involved in each failure mode. Despite the wide use of failure criteria, they present many difficulties regarding the accuracy of the results, because of undesired failure modes; plastic deformations and geometric nonlinearity of the parts; the effect of the residual tensions of the composite fabrication; and dispersions in the experimental results due to the heterogeneous nature of the materials [179].

Puck followed the failure theory framework of Hashin and proposed an elaborate scheme for implementation of his theory in Puck and Schürmann. As in Hashin, Puck’s theory (Puck, and Puck and Schürmann) recognizes a failure in UD composites to be in fiber failure (FF) and inter-fibre failure (IFF) modes [167,187,188,189]. For fiber failure, there are two modes of compression and traction. In the case of inter-fiber failure, there are Modes A, B and C, which include matrix fracture or fiber-matrix displacement. The inter-fiber failure modes were based on the Coulomb–Mohr’s fracture hypothesis which is appropriate for brittle fracture behavior of composite materials, wherein failure on a plane occurs when certain resistances, related to its cohesion and internal friction, are overcome [93,164,187,190]. Acoording to Wang and Zhao [165], Puck’s criterion has become a mainstream failure criterion for predicting responses of a composite subjected to impact loads.

Failure criteria are used to determine when material failure will occur. This concept can be applied in several areas, such as for mechanical tests, manufacturing processes, corrosive environments, etc.; for mechanical tests (Figure 14) and machining processes (Figure 15), they are the most used ones, where the criterion is combined with the material property degradation rule for the failure analysis.

Figure 14

Comparison of theoretical and experimental results. Reproduced with permission [191].

Figure 15

Progressive failure analysis with Hashin failure criteria for 45° fiber orientation. (a) Fiber–matrix interface detachment, (b) matrix fracture and (c) fiber rupture. Reproduced with permission [24].

5. Types of Elements Applied in Composite Modelling

For the analysis of a component, it is necessary to create a mesh on it [192]. The mesh is composed of elements and nodes. While the elements are subdivisions of the analyzed structure, the nodes are the connections of these subdivisions [21,193]. There are several types of elements; for example: bar, beam, columnar, triangular, quadrilateral, plate, shell, solid, etc. [21,194,195]. However, regarding an efficient and effective analysis for composite materials, there are four types mostly chosen: solids, beam, plate and shell.

Solid elements are the least used one for composites, because they require a model with many layers or a costly and time-consuming full-size structure, becoming consequently, unfeasible [196]. Besides these reasons, if the laminate thickness is very thin, layers constructed with solid elements can result in ill-conditioned equations. These factors lead to the use of other elements with lower computational demand and well-conditioned equations.

The beam element can be defined as having one of its dimensions larger than the others. One of the axes is defined along the longer dimension, and a cross-section taken perpendicular to this axis is assumed to vary smoothly along the beam length [74].

According to Bauchau and Craig [74] civil engineering structures often consist of an assembly or grid of beams with T or I-shaped cross-sections. A large number of machine parts are beam-like structures as well: lever arms, shafts, etc. Finally, several aeronautical structures such as wings and fuselages can also be treated as thin-walled beams.

Bauchau and Craig [74] cited that long and slender aircraft wings can be analyzed, as a first approximation, like beam structures, but a more refined and detailed analysis should treat separately the upper and lower skins of the wings as thin plates or shells supported by ribs and longerons, or stiffeners. Nevertheless, aircraft wings with small aspect ratios cannot be treated as beams because two their dimensions are larger than their thicknesses. But most of the time they can often be represented as plates. The aircraft fuselage is also constructed of thin-walled structures stiffened with ribs and longerons, and the thin-walled portions between the stiffeners can be drawn as thin plates. Last but not least, thin-walled beams can be modelled as plates when considering a localized behavior induced by attachments or supports.

Both plate and shell are considered two-dimensional or surface elements because two their 2D dimensions (length and width) are much larger than their thicknesses, which are given by the number of layers in their laminates [21,22,27,75,77,196,197,198,199,200,201]. Because of this, even with classical theory mathematically differentiating these two elements, the terms plate and shell are often used interchangeably, assuming that a plate element is flat, but when curved, it would become a shell element [21,74,197,198].

According to Barbero [77] whenever the thickness coordinates are eliminated from the general equation, they create a 3D problem in a 2D design. The author even mentioned that modelling laminate composites differ from any conventional materials modelling in three aspects:

The constitutive equations of each layer are orthotropic;

The constitutive equations of the element depend on the kinematic considerations of the plate/shell theory employed and its implementation on the element;

The symmetry of the material is as important as the geometry and symmetry of the loading when trying to use conditions of symmetry in the models.

Regarding structural composites, plate and shell elements are the most used types [74,202]. However, according to Tenek and Argyris [22] the use of a 2D element may neglect the flexural stiffness of a material. This issue does not exist in solid elements once they are a 3D type of structure; i.e., the thickness is included in the general equation.

5.1. Plate Element

A plate is a 2D plane solid element whose thickness (h or t), usually measured in the z-axis direction, is a lot smaller than its length and width, which are located in the xy plane [21,203,204,205]. The plate element can support actions that promote transverse flexion; in addition, it has two bending moments and one torsional moment [21,75,196,201].

The most classic examples of plates are slabs; buildings’ floor slabs; bridge decks; sides of rectangular water tanks and other fluid retaining structures; and tables. They transmit the loads that act on the normal direction of their midplane (z-axis) [75,199,201]. The two main theories that describe a plate element behavior are Kirchhoff theory and Mindlin theory, both being based on kinematic hypotheses [21,75,196,201,206,207,208,209].

5.1.1. Elements of Kirchhoff Theory

Classic Kirchhoff theory, also called classical theory, is analogous to Euler–Bernoulli’s beam theory; it is employed in the study of thin plates [74], whose relationship between the smallest plate gap and the thickness (t) is less than twenty [75,196,201,208]. This theory considers the thickness to be inextensible and neglects the shear strain deformations, assuming that a normal line segment at the mean surface remains rectilinear and perpendicular to the surface after deforming the plate (Figure 16) [75,206,208,210,211]. Altenbach and Eremeyev [203] added that Kirchhoff theory considers that the plate is made of a homogeneous, isotropic, linear elastic material. It assumes the validity of generalized Hooke’s law. According to Vaz [201], Altenbach and Eremeyev [203], Saliba et al. [212] and Schneider et al. [213], the kinematic hypotheses of Kirchhoff’s theory for plates with total isotropy are:

Figure 16

(a) Section of plate of thickness t, under transverse loading per unit area, where is a transverse displacement of a point of the mean surface, and are the rotations normal to the same point according to the x and y axes. Reproduced with permission [75]. (b) Field of displacements according to Kirchhoff’s plate theory. Reproduced with permission [212].

Any point P (x, y) on the average surface of the plate moves only in the z direction—that is, it has only vertical displacement w (x, y);

The normal stress in the z-direction (σz) is negligible;

The longitudinal strain is zero at any point on the plate, i.e., = 0;

A straight and normal line to the average surface before loading and cutting the median plane of the plate at point P (x, y) remains straight and normal to the plane tangent to the average surface at that point after loading, and therefore, the shear deformations e are zero.

Based on the third hypothesis above, it is possible to define expressions that describe the displacement fields of the plates, which are represented in Equation (8) [199,201].

According to Vaz [201] it is possible to obtain deformations for an infinitesimal plane element from z-dimension parallel to the mean plane of the plate (Equation (9)), ε being the deformation vector of a given point of the plate, and is the vector containing the curvatures of Kirchhoff’s theory relative to a point in the middle plane of the plate which is on the same vertical line as the point where ε was calculated.

Vaz [201] also pointed out that for plate element, there are actions from bending moments and shearing forces, and the shear moment can lead to vertical shear stresses, and consequently, distortions and , which are zero, as observed in Equation (10). That is the reason why this theory can only be applied to thin plates.

Soriano [75] and Vaz [201] cited that through the constitutive relations, it is possible to obtain the general equation for the element (Equations (11) and (12)). where:

In Equation (13), the summarized form of the equation is given, where {M} is the vector of moments at a point on the middle surface of the plate, [E] is the plate bending stiffness matrix by Kirchhoff theory and {ε} is the deformation vector at a given point [75,201].

Considering the previous equations, it is possible to obtain the deformation energy for an element (Equation (14)); the total potential energy (Equation (15)), where is the transverse force per unit of positive area in the z-direction; and the general equilibrium equation or the minimum potential energy principle (Equation (16)) [22,75,202,204].

Through Kirchhoff theory, two elements were developed, those being the rectangular and triangular [75,199,201,212,214]. For rectangular elements there are two types; one element is nonconforming and has three degrees of freedom per node: a vertical displacement (w) and two rotations (θx and θy) (Figure 17). The other has four degrees of freedom per node: a displacement vertical (w), two rotations (θx and θy) and a curvature (w, xy) [75,212,215,216,217,218,219].

Figure 17

Rectangular plate element by Kirchhoff’s theory.

Soriano [75] cite that the rectangular element does not have a constant shear state. Because of this, Kirchhoff’s triangular element was created, and it can have six or nine nodal displacements. Saliba [214] and Saliba et al. [212] cite that Kirchhoff’s nonconforming triangular element with nine terms was developed by Cheung et al. in 1968 and that the element has the three degrees of freedom per node: a vertical displacement (w) and two rotations (θx and θy) (Figure 18a). Soriano [75] argues that Morley in 1971 developed the nonconforming but convergent Kirchhoff triangular element with six terms, which has only vertex transverse displacements and normal rotations at the midpoints of the sides (Figure 18b).

Figure 18

Non-conforming triangular elements. Reproduced with permission [75].

Soriano [75] and Saliba [214] report the existence of other elements with discrete constraints and conforming elements, those being three triangular and one rectangular element (Figure 19).

Figure 19

Compliant thin-plate elements. Reproduced with permission [75].

5.1.2. Elements of Mindlin Theory

Mindlin or Reissner–Mindlin theory for plates is equivalent to Timoshenko theory for beams, in which the main hypothesis is that the cross section of the beam remains flat, but not necessarily perpendicular to the tangent of the elastic line when deformed [21,203,207,209,210,213]. Mindlin theory is hierarchically superior to the Kirchhoff because it presents a three-dimensional solution and can be applied for both thin plates and spines [75,206,207,208]. Mindlin theory is a shear-deformable plate theory [203].

According to Soriano [75], Altenbach and Eremeyev [203], Saliba et al. [212] and Vaz [201], the kinematic hypotheses of Mindlin’s theory are:

Any point P (x, y) on the average surface of the plate moves only in the z direction—that is, it has only vertical displacement w (x, y);

The normal stress in the z direction (σz) is negligible;

The vertical longitudinal strain is zero at any point on the plate—i.e., εz = 0;

A straight and normal line to the average surface before loading and cutting the median plane of the plate at point P (x, y), remains straight after loading, and straight but not necessarily normal to this plane, after deformation.

The change in the fourth hypothesis reflects in the displacement field of the plate (Figure 20), even though it does not link the rotation of the vertical line passing through P (x, y) to the derivatives of the vertical displacement w (x, y) [201,212,220].

Figure 20

Field of displacements according to Reissner–Mindlin plate theory Reproduced with permission [212].

According to Vaz [201], based on the third hypothesis, it is possible to define the expressions describes the plate displacement fields (Equation (17)).

Vaz [201] also mentioned that is possible to obtain the element deformation (Equation (18)), in which in the summarized equation, the deformation vector ε is subdivided into two vectors: associated with the bending moments and associated with the shear forces; is the transformation matrix subdivided into Tb and Ts; M is Mindlin curvature vector, divided into kb and ks.

Vaz [201] points out that now that the curvatures associated with the deformations and are only null if presented in Equation (19). The author cites that when it is valid, by the Kirchhoff hypothesis the rotations θ are given by the derivatives of w. As the γ distortions are not necessarily zero in the Mindlin theory, shear stresses and shear stresses will also not be zero.

Soriano [75] and Vaz [201] cite that through constitutive law, it is possible to obtain a general equation for the element (Equations (20) and (21)). where

Equation (22) presents the summarized form of the equation, where {M} is the momentum vector at a point on the mean surface of the plate, [E] is the plate bending stiffness matrix by Kirchhoff theory and {} is the bending deformation vector of thin plates [75].

For an isotropic element given by Equation (23), where K is the shear factor and G is the transverse modulus of elasticity, β is the shear deformation and KA represents a reduced area. In Equation (24), a new summarized equation for Mindlin theory is shown, where represents the general shear stress rotation or deformation K [75].

Considering the previous equations, it is possible to obtain the deformation energy for the element (Equation (25)) and the total potential energy (Equation (26)), where is the transverse force per unit of a positive area in z-direction and is the identity matrix of the elastic coefficients in the general matrix [E] [75].

Soriano [75], Saliba et al. [212] and Vaz [201] mentioned that unlike Kirchhoff’s theory, the rotations θx and θy are independent of the displacement w (x, y). This independence allows us to formulate C0-continuity elements.

Through Reissner–Mindlin theory, several elements have been developed, which can either be curved or not. In Figure 21 some of the isoparametric elements applied to the plate theory are presented—independent approximations for w (x, y), θx (x, y) and θy (x, y), are easily written due to general parametric FEM formulation [75,212,214]. For elements Q4, Q8, Q9, Q16, T3, T6 and T10, all nodes have three degrees of freedom: one vertical displacement and two rotations. However, the element Q9H, called heterosis, has the peripheral nodes with three degrees of freedom and the central node with only two rotations

Figure 21

Finite elements based on the Reissner–Mindlin theory.

Soriano [75] mentioned that for some reason, low order elements are subjected to locking or convergence ratio reduction. In order to identify the susceptibility to locking and the quality relation between the elements, the heuristic beam restraint ratio is generalized. For this, two constraints of shear stress are associated with each point of stiffness matrix integration; one related to and another to . However, support constraints on the spatial distribution are not considered, these being dependent on the geometric distortion of the element; i.e., this distortion influences the ability of the element to represent constant or zero shear strain deformations.

5.1.3. Theory of Kirchhoff versus Theory of Mindlin

Kirchhoff theory is suitable for a thin plate; Reissner–Mindlin plates; and thin and thick plates (also called semi-thick), but its application to thin plates requires special attention in order to make the elements capable of representing real-life cases [75,196,206,207,212,221]. Schneider, Kienzler and Böhm [213] cited that well-established standard theories for (linear geometry) homogeneous isotropic plate bending problems are: Kirchhoff theory, for neglecting the influence of shear deformations only suitable for very thin plates, and Reissner–Mindlin theory, which considers the influence of shear deformations and is used for thick plates.

However, there are more factors to be considered, such as static or dynamic behavior, a plate made from a single material and layers coming from distinct materials (sandwich or laminated). Shear strain consideration is the most important issue in a dynamic behavior and/or in a sandwich plate [75].

5.2. Shell Element

A shell is a two-dimensional planar solid whose thickness (h or t), usually measured along the z-axis, is much smaller than its length and width, both located in the xy plane [21,22,75,198,203]. This element is curved and can withstand bending and membrane effects, consisting of an average surface deformation for the element located on the same surface [21,75,201,222,223]. Examples of shell structures include acoustic shells, stadiums, large-span rooves, cooling towers, piping systems, pressure vessels, aircraft fuselages, rockets, water tanks, arch dams and many more. Even in the field of biomechanics, shell elements are used for the analyses of the skull, crustaceans’ shapes, red blood cells, etc. [75,109,199,200,224,225].

The shell element has probably generated more academic work in FEM technology than any other topic; however, it shows more computational barriers among all continuous structural elements, due to its curved geometry and the larger number of parameters involved [75,199,200,224].

According to Barbero [77], most of the composite structures are modelled using plate and shell elements. According to the author, this happens because, beyond reducing the numbers of nodes and elements, when compared to the solid element, it makes the modelling of thick laminates easy (Figure 22).

Figure 22

(a) 20-node isoparametric solid element; (b) reduction to eight nodes using a shell element; and (c) b-node representation.

Its shape allows certain membrane tensile systems to act parallel to its tangential plane and become primary deformation carriers. In fact, the analysis of many fine elements is based solely on shell membrane theory, neglecting their flexural stiffness [22].

Mathematically, the shell element model is similar to the plate element, since it is common to consider null the transverse normal stress component [21,22,75,198]. The shell geometry can be defined by its average surfaces or just one of its outer surfaces, called the reference surface, along with the thickness of each point. In general, the average surface is used as the reference surface [75].

According to Soriano [75], in the case of a shell element, bending is associated with the resultant loading forces (Mx, My, Mxy, Qx and Qy). In the case of small thicknesses, shell curvature radius expressions are identical to those from plate elements. Tensile components from a membrane effect, are the same as those occurring in the plane stress state, although they are considered by their results per unit length of the reference surface (Figure 23).

Figure 23

(a) Free body diagram for the equilibrium of in-plane forces and (b) free body diagram for the equilibrium of transverse shear forces.

5.2.1. Shell Theories

There are basically three coherent approaches for the shell element analysis. (a) Shell structure is faceted with flat elements, (b) via elements formulated on the basis of curved-shell theory or (c) degenerated three-dimensional elements [221,226,227].

The Theory of Flat Plate

The general theory of thin shell or flat shell was presented by H. Aron in 1874 and by A.E. Love in 1888, but it was only applied to solve engineering problems a century later. Similarly, to plate theory, plane shell theory basically differs from the idea of shear stress deformations, because in this case the analysis of shells with these elements is performed by superposing the membrane stiffness due to the plate elements [75,228]. Depending on the type of problem analyzed, the solutions obtained may depend on the discretization degree.

In spite of the presented difficulties (discontinuity in the momentum of interface elements) these elements are applied in linear and nonlinear shell analyses [226]. Plane shell theory can be divided into sub-theories under the assumptions of Reissner–Mindlin (first order theories), higher order theories and discrete layer theories [75,222,229].

In this theory, shear stress deformations are neglected, assuming that a line segment perpendicular to the average surface of the shell remains straight and normal to this surface after its deformation. Differently from the Reissner–Mindlin hypothesis, which assumes that the segment remains straight, not normal to the middle surface though [75,229].

In higher-order theory, nonlinear polynomial laws are adopted to define the segment bending after deformation; physical models, though, are better represented in thematic models than in other theories. Discrete layers theory is suitable for laminated shells, with a linear displacement field adopted by segments and thickness imposing appropriate contact conditions at the ply interfaces [75].

Three-Dimensional Elements

When a shell analysis is done with three-dimensional finite elements, many numerical difficulties may occur due to the discretization along with thickness, leading to an equation system with a greater number of unknowns compared to the degenerate model using a reference surface [75,226].

All shell models have uncertainties when compared to three-dimensional elasticity theory; nonetheless, they have the advantage of operating over very small magnitudes, allowing easy calculation of stresses [75].

Degenerate Shell Element of the Three-Dimensional Element

In the three-dimensional, degenerated shell element approach, simply known as the degenerated shell element, the element behavior, towards independent displacements and rotations, is degenerated from three-dimensional tensions and deformations [226].

These elements have the advantage of requiring only a C0-continuous function, once equilibrium equations are second-order differential equations. The degenerative concept of finite element formulation was extended by several authors into a linear and nonlinear analysis of anisotropic laminated composite structures [226].

Modelling using the degenerated shell is adopted because it reduces computational time, referring specifically to data provision and analysis, as well for numerical reasons [75,226].

5.2.2. Shell Element Types

The classical shell formulation requires displacements of a fifth-degree polynomial; consequently, a high number of nodal parameters are needed for both thin and shallow shells. The most practical solution is to develop shell elements starting with associations of plane stress elements and plate flexion, or by using the degenerated, three-dimensional curved element on a surface, adopting kinematic and mechanical constraints (Figure 24), resulting in in-plane elements, curved elements (with Reissner–Mindlin hypotheses) and an axisymmetric shell with an asymmetric loading [75,221,222].

Figure 24

(a) Association of elements and (b) degenerate shell.

Flat Elements

According to Soriano [75], it is possible to combine plate elements with flat plane stress elements; thus, it is worked out with the plane stress resultants, considered plane elements. The plane element generates an approximation of the geometry in the curved shell discretization, replacing it by a set of facet elements, aside from displacement field approximations, inherent in finite element analysis.

This type of discretization requires a large number of elements with a refined polyhedral surface approaching the original mean surface. In this case, the triangular elements (Figure 25) better represent a double-curved shell geometry than the quadrilateral elements, which in this case are more interesting for single-curved shells and flat-shell discretization [75].

Figure 25

Triangle element in xy plane. Reproduced with permission [75].

Considering many contributions from various elements, it is possible to determine a global stiffness matrix (Equation (27)), with [λ] representing a three-dimensional rotation matrix [75].

Curved elements (with Reissner–Mindlin Hypotheses)

According to Soriano [75], in 1968, Ahmad, Irons and Zienkiewicz developed a degenerated shell element, starting from a three-dimensional curved element. This degeneration became widely used with the knowledge of reduced/selective integration and evolved with complete integration of the mixed formulation, aiming to obtain robust elements.

Cook [21] cited that for curved elements, one can start from the middle surface, with the geometric definition and displacement field, adopting Reissner–Mindlin hypotheses, instead of three-dimensional explicitly degenerated element on its average surface by imposing a normal condition which maintains it straight, but not necessarily normal to this surface according to Reissner–Mindlin theory (Figure 26). In general, any other reference surface which is not necessarily the average one may be used as an outer surface. However, the average surface is usually adopted in the case of single-layer shells with homogeneous thickness [75].

Figure 26

Curved element. Reproduced with permission [75].

According to Soriano [75] the main advantages of using the degenerate shell are:

Working with the shell hypothesis from the beginning, obtaining, in a simple way, a wide range of elements;

Developing curved elements that only need C0 continuity;

Using only linear displacements and rotations as degrees of freedom, making it possible the use shell elements to discretize beam and plate elements;

Considering the effect of shear strain on a wide variety of thicknesses.

Equation (28) shows the parametric form for a three-dimensional geometry, starting from the mean surface of a curved element, where ζ is the dimensionless coordinate of the z-axis with values at the outer surfaces of ±1; is the thickness at the nodal point i; are the z-axis directional cosines also at point i, components of the vector .

Assigning the displacements according to the local axes x, y and z for u’, v’and w’ respectively, there are deformation vectors in a global reference and deformation components, respectively [75]. Based on these considerations and through generalized Hooke’s law, (assuming ), one can obtain the general equation for a degraded shell (Equation (29)).

From the curved shell element formulation with Mindilin theory for plates, Soriano [75] highlights the following differences:

For plates, initially, shear rotations were separated from the plate and worked on the tension-deformation relationships with the resultant stress, excluding, consequently, integration along with thickness in the rigidity matrix and nodal forces equivalence expressions;

For shells, working with total rotation and stress components results in expressions of rigidity matrix and equivalent nodal forces that require integration along with thickness. Note that plate elements could also be formulated the same way.

Asymmetric Shell with Asymmetric Loading

The use of structural axisymmetric (geometry and support conditions) and asymmetric loading leads to a simpler discrete model than the corresponding three-dimensional one. This simplicity is linked to a specific geometry and a smaller number of variables to be determined [75,230].

In 1963, Grafton and Strome presented the axisymmetric shell discretization, with asymmetric loading, in truncated cones corresponding to rectilinear finite elements according to a meridian shell, which has two nodal points and three displacements per node. Jones and Stone in 1966 modified the work from Grafton and Strome, considering curved elements according to the meridian; both authors considered the thin shell theory though [75].

These elements of revolution (conical and curved elements) have nodal circles and not nodal points, as for plate elements, and in general, there are two nodal circles per element, which has two translations (radial and axial) and one rotation [21]. Figure 27 represents the axisymmetric shell element, where x and z axes are respectively tangent and normal to a meridian of a mean surface at each point r, θ and Z with radial direction and Z to axial [75].

Figure 27

Asymmetric shell element. Reproduced with permission [75].

The equation that defines the displacements’ interpolation is represented in Equation (30), where u is the radial displacement, w is the axial displacement and is the rotation of the nodal point i, according to the circumferential direction [75].

Through all these considerations regarding stress and strain, the local and global referential are obtained by exclusion of , and by the exchange of y with θ, resulting in the axisymmetric shell general equation (Equation (31)) [75].

5.3. Cohesive Elements

Cohesive elements, also called decohesion elements or interface elements, are useful in modelling adhesives, bonded interfaces, delamination and rock fracturing [58,65,231,232,233,234,235]. The constitutive response of these elements depends on the specific application and is based on certain assumptions about the stress and strain states that are appropriate for each application area. The nature of the mechanical constitutive response can be broadly classified based on [232,236,237]:

Continuum-based modeling;

Laterally unconstrained adhesive patche;

Traction-separation based modeling.

In these approaches, cohesive elements are used to represent the behavior of a fracture, while traditional volumetric elements represent deformations of the continuous medium. Cohesive elements are inserted at the interfaces between pairs of adjacent volumetric elements in the finite element mesh (Figure 28) [238].

Figure 28

Interface element analysis of delamination in a double cantilever beam. Reproduced with permission [64].

5.3.1. Continuum-Based Modeling

According to Joshi, Pal and Chakraborty [237], continuum-based modeling is used when the cohesive zone has a finite thickness such as a joining of two surfaces with the help of adhesive material such as glue (Figure 29). The thickness, stiffness and strength of a cohesive zone can be estimated using experimental methods. In the case of continuum modelling, one directs stress in the direction of thickness, and two shear stresses mutually perpendicular and along the plane of the adhesive are present.

Figure 29

Progressive failure process of the adhesive layer in the lap-shear joint, (a) damage initiation at the overlap edges (SDEG-% = 0), (b) propagation towards the joint centre (SDEG-% ≈ 40), (c) joint failure (SDEG-% = 100). Reproduced with Creative Common License [239].

The cohesive elements model the initial loading, the initiation of damage and the propagation of damage leading to eventual failure in the material [232].

The continuum-based can be applied in 2D and 3D problems. In 2D problems the continuum-based constitutive model assumes one direct strain (through-thickness), one transverse shear strain and all stress components to be active at a material point. In 3D problems it assumes one direct strain (through-thickness), two transverse shear strains, and all stress components to be active at a material point (Figure 30) [232,237].

Figure 30

Spatial representation of CH3D8 (eight-node three-dimensional) cohesive element. Reproduced with Creative Common License [240].

5.3.2. Laterally Unconstrained Adhesive Patche

This approach is appropriate for modelling joints with gaskets (Figure 31). The macroscopic properties of the gasket, such as strength and stiffness, are used for the analysis. Only unidirectional stress along the through-thickness direction is considered in the analysis. The nonlinear and hyperelastic behavior of the materials used for gaskets—rubber, foam, etc.—can be captured in the constitutive relations used for the modelling techniques for laterally unconstrained adhesive patches [237]. The constitutive responses of gaskets modelled with cohesive elements can be defined using only macroscopic properties such as stiffness and strength [232].

Figure 31

Typical application involving gaskets, (a) Reproduced with permission [241] and (b) Reproduced with permission [242].

5.3.3. Traction-Separation-Based Modeling

The modelling of bonded interfaces in composite materials often involves situations wherein the intermediate glue material is very thin and for all practical purposes may be considered to be of zero thickness. Therefore when the macroscopic properties of the material, such as the rigidity and strength of the adhesive material, are not important, traction-separation-based modelling can be used (Figure 32) [60,63,232,237,243,244].

Figure 32

Damage of cohesive elements in the wake of delamination front. Reproduced with permission [246].

In cases of the macroscopic material properties are not relevant directly, the analysis should be based on concepts derived from fracture mechanics—such as the amount of energy required to create new surfaces [232,245].

The cohesive elements model (CZM) models the initial loading, the initiation of damage and the propagation of damage leading to eventual failure at the bonded interface. The behavior of the interface prior to initiation of damage is often described as linearly elastic in terms of a penalty stiffness that degrades under tensile and/or shear loading but is unaffected by pure compression [232,247].

According to Abena, Soo and Essa [60], the limitation of this approach is the inability to represent the beginning of the damage and the propagation of the failure under compression and the inability to produce any stress related to a membrane response. In contrast, elements representing the surrounding phases (matrix and fibre) are able to fail under compression and a membrane response, and are consequently deleted during the analysis. Therefore, the cohesive elements could remain in the model even if their surrounding elements fail. When this happens, the cohesive elements lose their aim, since they are not linking matrix and fibre any more, and they also usually experience excessive distortion since their nodes become free to move [60].

During modelling of the joints under traction-separation technique, before damage initiation, linear elastic behavior is assumed [60]. In 2D problems, the traction-separation-based model assumes two components of separation (one normal to the interface and the other parallel to it), and the corresponding stress components are assumed to be active at a material point. In 3D problems there are three components of separation (one normal to the interface and two parallel to it), and the corresponding stress components are assumed to be active at a material point [60,232].

The linear elastic behavior before the initiation of damage can thus be governed by the constitutive relations, as given by Equation (32) [60], where is the normal stress along the local direction 3 (through-thickness), and and are the shear stress components along the local directions 1 and 2, respectively. The is the normal strain along the local direction 3, and and are the shear strain components along the local directions 1 and 2, respectively.

The failure mechanism consists of damage initiation criterion and damage evolution law. The damage initiation can be governed by the criteria of maximum stress and maximum strain [60]. Regarding damage evolution law, different cohesive laws have been proposed in the literature, but normally, assumptions of zero adhesive thickness are made (Figure 33) [234,248,249].

Figure 33

Typical traction separation laws according to: (a) Needleman 1987, (b) Needleman 1990, (c) Hillerborg 1976, (d) Bazant 2002, (e) Scheider and Brocks 2003, (f) Tvergaard and Hutchinson 1992.

Once the damage criteria are met in any one of the modes, then the stiffness starts degrading, causing a gradual failure. This is also called as softening [237,250].

According to Arafah [248], Schwalbe, Scheider and Cornec [251] and Budiman et al. [252], the traction separation law (TSL) can be described by the following parameters: the critical separation ()—that is, the maximum displacement jump across the crack at which the cohesive element becomes completely broken; the cohesive strength (), which is the maximum traction at the crack plane; and the cohesive energy , which is the amount of energy consumed to create new crack surfaces (i.e., separation energy similar to Griffith’s fracture concept). The cohesive energy can be calculated from the area under the traction separation law T (), as in Equation (33).

The initial response is assumed to be linear, and once a damage initiation criterion is met, damage can occur according to a user-defined damage evolution law (Figure 34) [249].

Figure 34

Typical traction–separation response. Reproduced with Creative Common License [249].

In the case of mixed mode loading, a tangential separation mode, usually designated Mode II and Mode III, accompanies the normally considered crack opening (Mode I) [233,251]. In linear elastic fracture mechanics, a phase angle can be define by Equation (34) Figyre, where KI and KII denote the stress intensity factors for crack opening Modes I and II respectively [250].

In the context of the cohesive model, a tangential displacement represents the additional shear mode and is superimposed on the displacement normal to the crack plane (or plane of expected damage in the absence of a pre-existing crack) (Equation (35)) [250].

5.3.4. Cohesive Element Types

The cohesive elements can be defined as 2D or 3D elements. The 2D has four nodes and two integration points with a linear displacement formulation. The 3D can be constructed using eight nodes and four integration points with linear displacement formulation (Figure 35). The local coordinate system of the cohesive element could be defined with respect to the initial configuration or the actual configuration (i.e., moving coordinate system) (Figure 35) [240,248]. For the aim of calculating the stresses and separations of the cohesive element, they are connected to the adjacent continuum elements by sharing the respective common nodes (Figure 35) [248].

Figure 35

Cohesive element nodes and integration point positions, and connecting cohesive and continuum elements.

Modelling using the 2D cohesive elements has two options: the plane strain/stress and shell model. The difference between them is that the shell element is defined in the three-dimensional space. Therefore, any separation may be in-plane or out of the plane, and the in-plane direction must be defined by the user, which can be done by a fifth node, as shown in Figure 36 [250].

Figure 36

Cohesive element (a) plane stress/strain models and (b) shell models.

6. Main Applications of Finite Elements in the Study of Composite Materials

Every industrial sector feels over the years an increasing demand for innovative products that outperform competitors and meet market needs. For many design applications, resistant and lightweight materials are required, which makes laminates the ideal solution. However, new product developments and launches, and new technologies, cannot compromise product quality, reliability and speed [253].

The use of composites is growing rapidly in many industries, with new technologies becoming necessary for design, analysis and optimization. However, mechanical designers who work with steel and metals in general do not have many problems predicting the behaviors of these materials in use, because they are considered isotropic materials. Composite materials do not work the same way, once they are made from distinctive constituent materials [254]. Furthermore, their manufacturing is very expensive, making impractical the construction of real-size or scale-size prototypes.

In this context, numerical simulations evolved, mainly using the finite element method (FEM), and became a valuable tool to reduce project costs, enabling real-size components and equipment, boundary and entry conditions to be easily and quickly modelled. The usage of FEM to help engineering projects has increased over the past decades, due to computational capacity and data-processing performance, reduction project elaboration time and low costs.

Conventional software packages that use FEM in structure analysis, temperature, air and fluid fluxes, were not developed for composite materials, because of the fiber weaving and the microscale simulation approach needed. Acknowledging this deficiency, software providers started to include and add specific tools making composite analysis possible [254], owing to huge market demand for specialized tools to solve the constructive and operational particularities of these materials.

With this in mind, we can cite customized algorithms developed to work with well-known software such as Abaqus and HyperSizer, and add-on modules; for example, Ansys Composite Prep/Post, NISAII/COMPOSITE from Cranes Software Inc, FiberSIM from Siemens, Helius-Composite from Autodesk and GENOA from Altair and Laminate Tools from Anaglyph [255,256,257,258,259,260,261].

As already mentioned, the field of application of composite materials is very wide (Figure 37); consequently, the FEM software for laminate analyses has to meet all these specificities [262,263,264,265,266].

Figure 37

Some fields of FEM’s application to composites: (a) wind power, (b) sports, (c) automotive, (d) construction, (e) naval. Reproduced with Creative Common License [267], (f) aeronautics. Reproduced with permission [268] and (g) space. Reproduced with Creative Common License [269] industries.

6.1. Aronautical

As a way to demonstrate the importance of the use of FEM analysis in composite materials, there is the aeronautical sector, which could increase the percentage of laminates in aircraft—an advance impossible to be made with prototypes because of their high costs and very long manufacturing time.

After the 2000s, the aeronautical industry increased composites’ presence by structural weight by ≈40% [270] (Figure 38). This coincides with the increase in the processing capacity of computers and the evolution of numerical simulation software dedicated to composites.

Figure 38

Use of composite materials in aircraft over the years.

According to Vijaykumar Rayavarapu [271], Hindustan Aeronautics Limited R&D—manager, a major advancement of simulation for this sector is the possibility to predict the best geometries and materials for withstanding bird impacts, which according to Zhou, Sun and Huang [272] is a recurring and problematic phenomenon in aviation worldwide, because a bird impact causes structural, serious damage, particularly threatening the safety of aircraft and passengers; e.g., loss of aircraft and even the loss of lives.

Rayavarapu said, “In the past, the only way to determine whether composite aircraft components could withstand bird strikes was with time-consuming physical tests. Now, engineers use simulation to get the design right the first time. Bird strike simulation saves the companies design time and thousands of dollars per test of composite helicopter/aircraft components.”

According to Bouvet, Rivallant and Barrau [273] and Sun et al. [249], another problem for companies that use composites in their products is the impact caused by the tools falling on the components during their assembly or maintenance (low velocity impact). Composite laminates have brittle behavior and can undergo significant damage in terms of matrix cracks, fiber breakages or delamination. This damage is particularly dangerous because it drastically reduces the structural integrity of the structure, and at the same time can leave very limited visible marks on the surface impacted [274,275].

Through the start-up Carbon, another great advancement was achieved by the use of FEM to reduce cargo weight to increase aircraft fuel efficiency. For them, “Lightweighting is one of the most important trends in the aerospace industry today, as jet manufacturers and their suppliers work to reduce the overall weight of planes and improve their fuel efficiency, include reducing the weight of the cargo” [276].

Huang [277] and Yang [278] still commented on the introduction of nanostructured materials in aircraft, which will help to increase the resistance and decrease their structural weight. These composites taken with the multifunctional ones have been the subjects of several studies in the last decade, as they present significant improvements in mechanical, chemical, electrical and thermal properties; in other words, not only mechanical support functions, but also integrated functions [278,279,280,281,282,283].

In general, aeronautical companies use FEM in structural analyses, vibrations, dynamic buckling and fluid flows in the fuselage, wing, turbine and other parts (Figure 39) [284,285,286]. These analyses help determine the best geometry and material for each part of the aircraft, and understand how it will behave in service; for example, how the pressure and drag force will affect the aircraft during a flight [285,287].

Figure 39

Some applications of FEM in composites in the aeronautical industry: (a) fluid–structure. Reproduced with Creative Common License [288], (b) wing geometry. Reproduced with permission [289] and (c) tube rotor blade improvement. Reproduced with Creative Common License [290].

6.2. Space

The space sector, and the aeronautical sector, have benefited greatly from the advances of numerical simulation, because nowadays they can simulate a structure or component of a vehicle and test what the best geometry or material for it will be. In addition, there is the possibility of simulating various types of environments, such as the vacuum of space, solar radiation and even conditions from other planets (Figure 40) [90,291,292,293]. According to John Thornton, chief engineer of Astrobotic Technology Inc. [294], “Using design and simulation tools is possible designed and refined a lightweight spacecraft able to withstand static acceleration and dynamic random vibration loads of launch while maintaining an acceptable level of safety.” He still complements simulations for helping to reduce costs related to prototypes and physical testing.

Figure 40

Some applications of FEM in composites in the space sector: (a) satellite. Reproduced with Creative Common License [269], and (b) coupled launcher–satellite structure. Reproduced with permission [295].

6.3. Automotive

Another sector that has been applying composites to its products is the automotive one; they use FEM for the bodies of the cars and engines, mainly for impact tests. According to Jung [296] they predict damage under impact loads on automotive components, as impacts can reduce the structural integrity of composite components. Impact damage induces a variety of failure modes: matrix cracking, fiber breakage, and delamination. Among these failure modes, delamination is the most significant because the interlaminar mechanical properties are much weaker than those the rest of the composite, and it is difficult to visually detect this damage from the outside [296,297,298].

One of the applications of FEM in this scenario is to understand the behavior of the material during an impact, as in the work developed by Chiacchiarelli, Cerrutti and Flores-Johnson [299], in which they used layers with the same density but different compressive strengths in the direction of load, in a multilayered, rigid polyurethane foam block subjected to impact, since, by its excellent energy absorption capacity, it can be used as a high-level coating resistant to collisions, which makes it attractive for the automotive sector.

In the area of super-sport and competition cars, such as those for Formula One racing, it is very advantageous to use simulations to find the best geometry of a car, since a few seconds can define the winner (Figure 41) [253,282,287,300,301,302].

Figure 41

Some applications of FEM in composites in the automotive sector: (a) impact. Reproduced with Creative Common License [266], and (b) car racing.

According to Sandeep Sovani, director of Global Automotive Industry-ANSYS [253], “Simulation is key to solving issues upfront in the design phase. Companies that fine-tune auto body and chassis can reduce fuel consumption and build in reliability upfront in the design process with simulation.”

As examples the Sovani cite, “The Tier 1-supplier DENSO embeds CAE into all phases of its product development process, improving quality and reducing time to market along the way. Advanced virtual analysis enables such pacesetters to create category-changing; the KTM Technologies incorporated radical composites into a sports car, which called for new design, analysis and optimization technologies. Created using simulation, the product struck a fine balance between requirements, performance and costs while exceeding customers’ requirements; the ZF-TRW study the friction-induced brake squeal, for grow important as other vehicle noise sources are mitigated because of this the ZF-TRW engineers accurately simulated squeal and automated the simulation process while reducing time and money spent on validation testing. Performing simulation early in the design process helps to avoid costs associated with multiple prototypes, rework and tooling changes; and Valeo used nonlinear best practices to simulate thermoplastic snap- fits, leveraging HPC that shrunk simulation time by 50 percent” [253].

6.4. Naval

In the naval sector, the use of composites has been growing with the aid of FEM, with the main research being for marine propellers, ships and submarines (Figure 42). This research field uses FEM either to determine the mechanical properties or to predict the vibrations of its vessels, besides optimizing its structural weight [69,258,303,304,305,306,307].

Figure 42

Some applications of FEM in composites in the naval sector: (a) propellers. Reproduced with permission [303], (b) submarines. Reproduced with Creative Common License [267] and (c) ships.

6.5. Energy

In the energy field, composites have been gaining space in the wind sector (Figure 43a), mainly due to the fabrication of the blade turbines, and with the help of FEM it is possible to predict how the wind; the heat; and the sun’s radiation, centrifugal force and gravity load will affect the blades—and the ideal geometry to optimize the turbine [46,308,309,310].

Figure 43

Some applications of FEM in composites for the energy sector: (a) wind turbines. Reproduced with permission [310] and (b) pipes for the oil and gas industry.

In addition, Oliveira [311] cited a major innovation of the WEG company, which “wanted to replace the steel cover that contains the rotor coil head with an alternative material in a new line of turbogenerators. The engineering team explored the use of a pre-impregnated composite material in the form of a banding tape instead of a retaining ring. The result was a fully validated component with a lower rotor mass that is also cost- effective to manufacture. New materials, like the composites in this case, can reduce feedstock costs by as much as 77 percent and wound rotor manufacturing costs by 18 to 20 percent.”

Another branch of the energy sector that has turned its attention to laminates is the oil and gas industry (Figure 43b), as they have been looking for deeper oil fields, and therefore, alternatives to the traditional armor wires are needed, due to a more severe environment. The traditional armor wires made from carbon steel are capable of operating at sea depths of 2.4 km but pose limitations, such as tensile fatigue, buckling and corrosion and excessive weight; thus, one option is the flexible hybrid composite pipe [312,313,314,315,316,317]. According to Anderson et al. [314] the combination of ultra-deepwater (>1500 m) applications and large pipe diameter requirements presents severe engineering challenges for rigid pipe technology. To address those needs, General Electric, with the support of a research partnership with Secure Energy for America, has a development program to generate a flexible pipe with an internal diameter of greater than seven inches for ultra-deepwater applications, and for this they use the FEM.

6.6. Civil Construction

In civil construction, FEM has been used over the years for the structural calculations of buildings, bridges loads, etc., both of brickwork and wood (Figure 44). For designers, the advance of FEM in composites has helped to elaborate hybrid constructions, wherein fiber-reinforced polymeric composites are used as reinforcements, for example, of bridges and columns, besides helping analyze the effects of natural events such as earthquakes, tornados and hurricanes on builds [142,318,319,320,321]. In addition, as in the aerospace sector, the new nanostructured composites have been helping to reduce the structural weights of the constructions [277,322].

Figure 44

Some applications of FEM for composites in construction: (a) Strengthening flexural solid timber beams with CFRP. Reproduced with permission [323] and (b) bridges.

6.7. Sports

Composite materials have not only been gaining a place in the large industrial sectors, but today they have been sports sector favorites as well, especially in the area of high performance or high-level sports. This field has been using composite materials to reduce equipment weight; take as examples, bicycles, tennis racquets, skis and even racing prosthetics (Figure 45). This field uses FEM in its projects as a way to determine geometries that have high resistance and lower weight to improve the athletes’ performances [253,324,325,326,327,328,329,330].

Figure 45

Some applications of FEM for composites in the sport sector: (a) a tennis racket, (b) a bicycle and (c) a prosthetic. Reproduced with Creative Common License [331].

For KTM Technologies GmbH company, “Designers must predict how well the finished product will perform in the real world, such as on a race track or road. Predicting failure, delamination, ultimate strength and other development variables is critical before prototype and manufacturing stages” [253].

Scholz et al. [324] wrote: “Composites typically possess a superior strength to weight characteristic compared to monolithic materials and offer excellent biocompatibility. They are, therefore, favorable for both hard- and soft-tissue applications and the design of prostheses. In particular, the development of specifically designed carbon fiber sports prostheses now allows lower-limb amputees to actively participate in competitive sports. Sensory feedback systems, porous composite materials for tissue engineering and functional coatings for metallic implants are further developments anticipated to be introduced in next-generation orthopedic medicine.”

6.8. Manufacturing

Due to so many composite applications, the manufacturing sector has studied ways to minimize the damage generated during manufacturing processes, especially in machining (Figure 46). This concern is due to the damages that the process causes to laminates, such as delamination, interlaminar fissures, fiber–matrix displacement and thermal damage [332,333,334,335].

Figure 46

Example of damage cause by machining, (a) Reproduced with Creative Common License [333] and (b) Reproduced with permission [334].

Another concern is the quality shown in the machined pieces; in this case, these pieces are perforated because the drilling process is the most used one for laminates, since most of the composite joints are bolted joints [332,336,337]. For this, static and dynamic mechanical tests are simulated in order to verify the damages caused by drilling and propagated by the test (Figure 47).

Figure 47

(a) Damage occurred in the FEM (left) and in the mechanical test (right). (b) Digital Image Correlation with overlapped shear strain (above) and FE damage (below) for the OHT coupon at 1,000,000 cycles. Reproduced with permission [338].

6.9. High-Performance Electronics

One of the most recent fields of composite application is high-performance electronics, because modern electronic devices rely on novel composite materials to achieve superior performance [339,340].

According to Sudhir Sharma, director of high-tech industry strategy and marketing of ANSYS, “Today, high-tech companies turn to advanced lightweight, yet strong, materials to create flexible mobile and wearable electronics. However, a range of complex issues must be considered when evaluating new materials—including electrical conduction properties, structural strength, dimensional stability over time and resistance to thermal build-up. The high-tech engineers simulate the electrical performance and thermal performance of electronic systems and devices. Design for manufacturability is also important” [313].

Sharma cited two examples of research using FEM in this area; the first case involved 3M, which published a groundbreaking study on how a novel embedded-capacitance composite material affected the electrical performance of a printed circuit board. The second is a study by the University of Pittsburgh and Carnegie Mellon University, wherein the engineers use FEM to assess the performances of new nanocomposites that have the potential to revolutionize power transformer technology.

Another research area that is growing in electronics regards polymer composites of high thermal conductivity, commonly used in many industries for renewable energy systems and electronic systems, such as solar cells, light-emitting diodes, Li-ion batteries and microelectronic packaging. This research is due to the potential applications as flexible polymers in electronic packaging and encapsulations, and satellite devices, because the local overheating and heat accumulation that occurs due to increasing power density may lead to the degradation and failure of functional systems [341,342,343].

7. Conclusions

In this article, the main contributions in the field of composite simulation were reviewed and presented, such as the specific characteristics, theories and constitutive relationships of the composite materials, and the types of elements and failure criteria used in each case. In addition, the main industrial sectors that have been using composite simulation as an important part of project development and study with these materials were also described.

At first, the definitions and characteristics that these materials can have were briefly commented on, as were the most usual configurations that are used in the various industrial sectors.

After that, the different types of models used in the simulations of composite materials, called micromechanics, homogeneous equivalent and the combination between them, were detailed. Starting with the definitions and concepts involved in micromechanical modeling, which individually approach the properties of each constituent material (fiber and matrix) and the interface between them, we presented a model with which it is possible to study local defects and failure mechanisms between the components; however, it demands high computational performance. Continuing with the definition of homogeneous equivalent modeling—it disregards the difference between constituents and treats the material as homogeneous, decreasing the demand for computational performance and making it impossible to study the effects locally.

After that, in the same section, the different levels of simulation that can be used in composite materials simulations were approached. Starting with simulations at the micro level, which focus on the interaction and interface between a fiber filament and the matrix, we then moved on to simulations at the meso level, approaching units of material volume and the relationships between the layers that make up the material; and finally, we described macro level simulations, which consider these materials as homogeneous.

Then, the constitutive relationships that can govern mechanical properties of these materials were discussed, presenting the differences between them one by one. We started with the constitutive relationship of an anisotropic material, followed that by an orthotropic material and finished with a transversely isotropic material.

In sequence, the article reviewed the various failure mechanisms that these materials can present, those being the failures in the fibers and in the matrix, due to loads that can be of compression, tension and/or shear. We also approached the various failure criteria proposed in research, such as Hill, Tsai–Hill, Tsai–Wu, Hashin-Rotem, Hashin, maximum stress, Hoffman, maximum stress, Hou, Puck–Schürmann, Chang-Chang, Linde, LaRC03, LaRC04, Maimi, Hart-Smith and Yeh-Stratton—each trying their best to get closer to the experimental results. Among all of those, maximum stress, Tsai–Hill, Tsai–Wu, Hashin and Puck–Schürmann were discussed more deeply, as they are the most commonly used.

In the ensuing section, the types of elements and theories that govern these were discussed, as were their differences and applications in the simulation of composites. Among all types of available elements, the elements of plate, shell and cohesive were highlighted, since these are the most commonly used.

The last section presents reports from companies and industrial sectors, such as aeronautics, aerospace, automotive, naval, energy, civil, sports, manufacturing and even electronics, which use finite element simulations with composite materials in very different cases. These cases start from simple mechanical test simulations, and go to more complex studies, such as topology optimization; impact; tribology; vibrations; fluid flow; environmental conditioning; effects of natural events (earthquake, tornado and drilling); prediction of damages resulting from the machining; and even electrical properties, such as electrical conduction and capacitance. That came together to demonstrate the great importance that simulations with finite elements of composite materials has presented in recent years.