A Virtual Testing Approach for Laminated Composites Based on Micromechanics

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Introduction
The last quarter-century has witnessed considerable research efforts in the mechanics of composites in order to understand and predict the behavior of these materials. The ultimate goal is the design of the materials/structures/manufacturing processes. The scientific literature on composites is immense: numerous books and a dozen of international reviews have dealt with the understanding, modeling, and identification of the thermomechanical behavior of composites, with particular emphasis on long-fiber laminated composites, which are the most important composites for the aeronautical and space industries and now automotive industries. Even in the case of laminated composites, the prediction of the evolution of damage up to, and including, final fracture remains a major challenge which is at the heart of today's virtual structural testing revolution. Virtual structural testing consists, whenever possible, in replacing the numerous experimental tests used today by virtual tests. In the early 2000s, there were two competing approaches: the micromechanics of stratified composites, see [1][2][3][4], in which one counts cracks. Another very active area was that of the structural analysis itself, but it is unfortunately limited to elasticity, especially delamination analysis [5][6][7][8] and optimization [9]. Over the past few years, numerous computational models have been introduced outside of these two families to simulate and interpret experimental tests but very few could P. Ladevèze ( ) • D. Néron • H. Bainier LMT (ENS Cachan, CNRS, Université Paris-Saclay), 61 avenue du Président Wilson, F-94235 Cachan Cedex, France e-mail: ladeveze@lmt.ens-cachan.fr; neron@lmt.ens-cachan.fr; bainier@lmt.ens-cachan.fr be seen today as predictive ones; among them, one can cite [10,11] and the LMT 1 -Cachan's damage-fracture mesomodel. This last answer to Virtual Structural Testing is in-between the "fracture criterion" vision, which is rather coarse, and the "micro" vision, which is very satisfactory in terms of physics, but beyond the capabilities of today's computation means.
More, LMT had been working to create a micro-meso bridge [12][13][14][15][16] in order to develop a complete synergy between micromechanics and macromechanics. The first application was the Computational Damage Micromodel introduced in [17][18][19], which is compatible with all the preexisting knowledge on the microand mesolevels and, contrary to its "micro" rival models, capable of simulating complex structures. A major difference compared to existing models is that all cracks are modeled in extenso regarding both their initiation and their propagation thanks to what one calls "Finite Fracture Mechanics" [3,20]. Therefore, the reproduction of experiments does not require the introduction of predefined sets of cracks based on knowledge derived from tests. However such a model leads to considerable computation costs and then cannot be used now in structural engineering calculations [21]. So, in the LMT-approach, this micromodel constitutes a foundation that can be seen as a Reference Virtual Material of a whole family of materials and not just of a single material. Using the micro-meso bridge, the LMT-mesomodel appears as an homogenized version of the micromodel and is sound in terms of physics [22]. This version has been implemented or is pending implementation into industrial codes for Virtual Testing. It constitutes a significant improvement over the standard version developed over some 20 years and available in many industrial structural analysis programs [23][24][25][26]. A last step has been done recently in [27] with the mesomodeling of the interaction between microcracking in plies and delamination which is a quite important observed phenomenon [28][29][30][31][32][33][34][35][36][37].
The paper is a revisit of our multiscale modeling approach. We will emphasize on the practical aspects related to the different stones: • a unified micromodel also involving fatigue, high-velocity, and oxidation, i.e., a Virtual Reference Material; • a general bridge between micro-and meso-mechanics; • the damage mesomodel used for structure computations as a homogenized model.
Some complements regarding identification and also kinking and crack initiation are given.
An another main question in Virtual Structural Testing is to compute engineering composite structures with such a mesomodel. This is not so easy because in the case of composites, damage and fracture involve many local instabilities; there are many cracks in the different plies and interfaces in particular splits, which are initiated and then propagated. For quasi-static loadings, to master the used localization limiters is not a simple task. The last part of paper deals with this main question; we will give our recent results to maximize the robustness of the computation strategy; in particular, we will describe a new technique to handle splits. Finally, the current capabilities and limits of this multiscale approach are pointed out as well as computational challenges that accompany Virtual Structural Testing.

The Main Damage Mechanisms
Damage mechanics of laminated composites is today well understood thanks to the numerous experimental and theoretical works relative to the so-called Micromechanics of Laminates (see the review papers [3,4,38] and in particular the books [1,2]). Figure 23.1 shows, apart fiber breaking, the different scenarios on the microscale (note that scenarios 3 and 4 are generally missing in micromechanics).
In most practical cases, the chain of scenarios follows Fig. 23.2. Scenarios 3 and 4 start, leading to a rather diffuse damage inside the plies and interfaces. Through a percolation phenomenon, transverse microcracks appear and then scenario 1 is active. The competition between transverse microcracking and local delamination ends with the saturation of scenario 1 and is relayed by the catastrophic development of scenario 2. Finally, the final fracture appears with fiber breaking and delamination.

Basic Aspects
The working scale in micromechanics is between the dimension of the structure and the diameter of a fiber ( Fig. 23.3). The state of the structure is in fact described as an assembly of cracked interfaces and cracked layers made with a "fiber-matrix" material assumed to be homogeneous or quasi-homogeneous. The difficulty comes from the great number of cracks. One starts with the initial state, where residual stresses occur and that can be calculated from the process simulation. A more pragmatic and standard approach simply consists in introducing a uniform negative variation of temperature and then the corresponding residual stresses are computed in elasticity.
where b 2 is a material constant which balances the influence of the transverse energy and the shear energy. For small damage rates and quasi-static loadings, we get The damage evolution is then defined by a linear function. Generally, d does not exceed a value close to 0.5 associated with an instability and there is a saturation phenomenon: microcracking starts, thanks to a percolation phenomenon. However for thermoplastic materials, this instability phenomenon can appear later, thanks to fiber rotations. The model remains valid for a rather large temperature range [40] and, at room temperature, a typical shear damage material function is given in Fig. 23.4.

Coupling Between Damage and (Visco)Plasticity
The microcracks lead to sliding with friction, and thus to inelastic strains, the matrix can also involve (visco)plasticity. The effective stress and inelastic strain are defined by: where " ij for i; j 2 f1; 2; 3g denotes the usual inelastic strain. It has been shown that very simple high fidelity material models can be obtained using such effective quantities. A very simple plasticity model describing inelastic strains is defined by the following elastic domain: Hardening is assumed to be isotropic, which means that the threshold R is a function of the cumulated strain p. p ! R.p/ is a material function, p being defined by: where a is a material coupling constant. The yield conditions are An example of such a hardening curve is given for the IM6-914 material in Fig. 23.5.

Identification
The main test is OE45 n ; 45 n 2S . To get the coupling coefficients a, b 2 , and b 3 , one performs, for example, the test OE67; 5 n ; 67; 5 n 2S .

Modeling of Delamination and Microcracking
One introduces minimum cracked surfaces, the characteristic length being the thickness of the elementary ply (see Fig. 23.6). All the finite energy release rates are associated with these surfaces. and

Propagation Criterion of Existing Transverse Microcracks
The propagation is driven by: The energy release rate related to microcracking could be computed simply by using the tunneling value. More generally, the computation of the different energy release rates could be greatly simplified by using analytical expressions. When an elementary surface is cracked, unilateral contact conditions with friction occur. The critical values are stochastic fields which are replaced, after discretization, by independent stochastic variables for which a modified normal law is introduced. Results are quasi-independent of the chosen probability law [17,18].

Delamination Cracking
The material interface, i.e., a small matrix layer between two layers of different orientations, is supposed to be independent on the fiber direction angle. This assumption has been experimentally confirmed taking into account intra-ply damages [41]. The interface is made of matrix described as an elastic cohesive interface with the following fracture criterion involving finite energy release rates: The material constants to identify herein are the classical critical energy release rates and the "transition" thickness N h which is between 2 and 5 times the thickness of the elementary ply. For delamination cracking, only one interface is investigated. A quick identification could be considering first that microcracking propagation could be interpreted as delamination for a [0, 0] interface. Experimental results [41] leads to: More, one has also:

Fiber Breaking
It is assumed to be brittle for traction behavior. A minimum volume to fracture is introduced as a cube of height h. So, we introduce where k is a material constant. Denoting the mean value over the cube of height h by h h h i i i, the criteria for traction and compression are Function`, which depends on the parameters˛C and m C , is detailed in Appendix. It describes the kinking damage in compression. Parameters˛C and m C can be identified, thanks to experimental and theoretical results given in [42,43]. Approximated as a brittle mechanism, the second relation of (23.16) could be written: A least, one has to identify the two thresholds Y t d f c , Y c d f c and the coupling coefficient k.

Structure Computation
Structure computations have been done in [21,44,45] using the LATIN multiscale computational approach and parallel computing. Results are very close to experiment and Fig. 23.7 shows such a comparison for a composite plate OE0 2 =90 2 s Fig. 23.7 Comparison experiment and computation for a holed composite plate OE0 2 ; 90 2 in tension [45] in tension [45]. In particular, one reproduces perfectly the splits in the 0 ı -ply. Unfortunately, the computation cost is today much too prohibitive for its use for structure computation.

Toward a Unified Model
Extension to fatigue has been done in [46,47] for fatigue impacts the fiber-matrix material. As the diffuse damage, although it exists, appears to be small, we have preferred to introduce in [46] a new internal damage variable Á f such that for cyclic loading, the critical energy release rates: where Á f follows classical fatigue laws in term of damage forces over the "minimal" surfaces introduced in the last paragraph. That is the only change in the model. An illustration is given in Fig. 23.8. Several items are less mature. Viscoelasticity has been added considering effective quantities [45] and in relation with impact problems in [48,49]. Numerous experimental results can be founded in [50]. Oxidation has been studied in [46,51,52]. Electrical behavior of laminated composites is studied in [53].

Extension
A first attempt to extend the concept of Reference Virtual Material to woven composite materials has been done in [54]. The RVM is then characterized by a precise description of the woven ply architecture (warp and weft yarns) and by the introduction of discrete microcracks and local delaminations. More, in woven   [46] composites, the additional scale of the yarn is present. This work extends the multiscale computational approach previously developed for laminates made of unidirectional plies to woven composites.

The Method
A rather complete bridge has been built in [12,13,15,16] and Fig. 23.9 describes the used two-scale computational scheme. At the left, the real structure submitted to a given loading is defined at the microscale. It is made with the fiber-matrix material and there are cracks, delamination, and transverse ones. The problem to solve is clearly a two-scale one, and then the solution consists in two parts: the large wavelength part for which the characteristic length is the structure dimension and the small wavelength part which has a characteristic length equal to the ply thickness.
A classical scheme to solve this two-scale problem is to separate the calculation of the two parts. In a first step, one determines the large wavelength part solving the so-called homogenized problem where the real structure is replaced by the homogenized structure. Its solution defines the mesoquantities. In a second step, the microquantities that make the small wavelength part of the solution are determined The equivalence should hold for any value of the residual which can be written in term of mesoquantities. The fundamental micro-meso link which defines the socalled homogenized structure holds exactly for the two basic problems. It could be written: where h h i i D 1 mes. / R dS. is any cross-section orthogonal to N 3 and compatible to the periodicity associated with the layer or interface containing . Practically for large crack densities, could be replaced by any large cross-section with respect to the plate thickness. is the projector on the plane orthogonal to N 3 .
A very large amount of calculations have been performed considering the practical ranges of the parameters and in particular of the microdamage variables. We have proved that this cell represents all engineering situations and that it can be approximated by two basic problems: the ply basic problem and the interface basic problem. The homogenized operators are quasi-intrinsic. They do not depend practically on the parameters of the upper and lower parts.

The Ply Basic Problem
The mesomodel of the fiber-matrix material should take into account the matrix transverse microcracking through the additional mesodamage variables N d 22 , N d 12 , and N d 23 . It follows that the diagonal terms of the energy written in (23.1) are modified considering the new moduli: To compute the additional mesodamage variables in terms of the micro ones, , C , and , one considers the microcell defined in Fig. 23.11. The central ply is a 90 ı -ply and the upper and lower parts are 0 ı -layers. The microdamage dimensionless parameters are Illustrations of this approach are given in [12,13]. In practice, one can neglect the influence of C and . Then, we get the intrinsic material functions:  The microcracking damage functions for a T700 M21 carbon-epoxy composite at room temperature [55] that are plotted in Fig. 23.12 for a T700 M21 carbon-epoxy composite. One has to use these functions and their derivatives. So here, we propose to simplify their use considering approximations over the practical range [0,0.8]: where A, A 12 , A 22 , and L are coefficients computed from the original functions.
Let us go further with the evolution law at the mesoscale. First, for a given , one associates the following damage force: and with an added coefficient . The role of this coefficient is to take into account the "not quite" periodical nature of the transverse cracking network. This is a shifting coefficient comparable to that used in [12,15] and very close to that introduced by Nairn and Hu [38]. Its identification could be performed considering the comparison between random and periodic behaviors done in [17][18][19], Another technique, proposed in [55], is based on a multiple cracking tests on cross-ply laminates. The identified value for classical laminates is 1:4 (23.28)

Remark
The external plies, and also the parts of plies adjacent to completely delaminated interfaces, behave regarding microcracking, as half a ply [15,38].

The Interface Basic Problem
At the mesoscale, the interface concept is modified: microcracking should be taken into account. However, it stays a surfacic material model and we have proved in [16,27] that it can be determined, thanks to the microcell defined in Fig. 23.13. The central interface is between two cracked layers: a Â /2-ply and a Â /2-ply. The upper and lower parts are 0 ı -layers. The interface is described at the microscale as a classical cohesive interface. The microdamage dimensionless parameters are • transverse microcracking rates: D h=D, 0 D h 0 =D 0 ; • local delamination rates: The problem to be solved involves periodic conditions. It is a 3D-problem which can be solved, thanks to two 2D-problems [16]. We have proved in [27] that its energy is equal to: where h I is the interface thickness, E and G the Young and shear moduli of the matrix.
One has also: where . / is a material function computed with an interface with Â D 90 ı . One can approximate it by: where I is a coefficient. For classical carbon-epoxy composites I is close to 0.5. Let us note that even if d 33 is equal to zero, interface shear stiffness can decrease strongly. The effective damage force, which is responsible for the increase of the interface's damage, is where m D and˛D are related to delamination initiation. They characterize the cohesive zone dimension which is larger than the thickness of the ply.

Delamination Criterion for Inplane Loading
The microdelamination initiation corresponds to the microcracking saturation. We have proved in [27] that the microdelamination (scenario 2) is instable for classical carbon-epoxy composites. So, in [27], we have introduced a brittle criteria at the mesoscale computed from the micro-meso bridge. However, the dimensionless microcracking density at saturation can be seen as a material constant. Numerous studies showed that this value is approximately intrinsic regardless of the cross-ply sequence chosen [13,29]. In practice, this parameter enables one to express the limitation of the cracking phenomenon due to the microdelamination initiation at matrix crack tips. Consequently, we propose the following criterion: where and 0 are the dimensionless microcracking densities of the adjacent plies to the interface. S is the saturation value which is around 0.8.

The Damage Mesomodel
An approach aiming at being able to answer Virtual Testing while maintaining a thorough and predictive description of the physics of composite degradation is the damage mesomodel developed at the LMT-Cachan since the 1980s [23][24][25]56]. Recent studies allowed to move from the standard damage model to a more complete description of the degradation mechanisms at each scale, based on micromechanical considerations. So now, the damage mesomodel proposed for Virtual Testing is the homogenized damage model of the Virtual Reference Material detailed described in Sect. 23.2. At the mesoscale, characterized by the thickness of the ply, the laminate structure is described as a stacking sequence of homogeneous layers through the thickness and of interlaminar interfaces (see Fig. 23.14). The main damage mechanisms are described in the 3D-continuum mechanics framework as: fiber breaking, matrix microcracking, and debonding of adjacent layers. The single-layer model includes both damage and inelasticity. The interlaminar interface is defined as a two-dimensional mechanical model which ensures traction and displacement transfer from one ply to the next. Its mechanical behavior depends on the angle between the fibers of two adjacent layers. A priori, 0 ı =0 ı interfaces are not introduced. The damage mechanisms are taken into account by means of internal damage variables. The mesomodel is then defined by adding another property: a uniform damage state is prescribed throughout the thickness of the elementary ply. This point plays a major role when trying to simulate a macrocrack with a damage model. As a complement, damage models with delay effect are introduced.

Fiber Breaking
Let us consider the traction and compression damage forces introduced in Sect. 23.2.2.4 but now the stresses are the meso-ones: dz. The damage evolution law is defined thanks to a damage model with delay effect. One has where q D sup.r T d f ;`.r C ;˛C; m C / d f / and with Parameters˛C, m C , Y t d f c , and Y c d f c have already been introduced in Sect. 23.2.2.4. Let us note that the compression behavior could be approximated by a brittle mechanism. The localization limiter parameters a c and c will be discussed in the next paragraph.

The Interface
The interface model which is a non-standard cohesive interface is completely defined in Sect. 23.3.2.2 as well as its interaction with the microcracking densities of the adjacent plies. This is an orthotropic modeling, the axis being the bisectors of the angle related to the fiber directions of the adjacent plies.

Localization Limiters and Numerical Parameters
An important issue is the objective prediction of final fracture which is rather well-understood nowadays. It is well known that classical damage models are nonconsistent. A visible lack is the abnormal sensitivity to the mesh of the finite element solution. Among the several remedies which have been proposed, we have followed an approach specific to laminate composites introduced in [25,56]. It is prescribed that the damage state is piecewise constant in the thickness of the laminate. Moreover, as a complement, we use damage models with delay effect combined with a dynamic analysis (further developments can be found in [26,57] and [58]): with q D ! d, ! being a function of the local strain.
There are two material constants c and a c . The variation of the force ! does not lead to instantaneous variations of the damage variable d. There is a certain delay, defined by the characteristic time c . Moreover, a maximum damage rate, equal to 1= c , does exist. A first identification consists of taking half the Rayleigh wave speed combined with the critical value of the energy release rate. Let us also point out here that a clear distinction can be made between this damage model with delay effects and viscoelastic or viscoplastic models: the characteristic time introduced in the damage model with delay effects is several orders of magnitude less than in the viscous case. This characteristic time is, in fact, related to the fracture process. For quasi-static conditions, the time discretization is not able to follow the c -scale; then, such a damage model with delay effect should be seen as a mathematical regularization.
Numerical tests done in [59] lead to: where T is the discretization time. An illustration done in [60] is given in Fig. 23.15 where a classical OE0=90 4 S tension test is simulated. It shows that the simulation reproduces the damage physics correctly. Until 1), transverse microcracking development is observed. Diffuse damage remains weak and is not shown in the damage charts. From (1) to (2), delamination develops very quickly and, in the end, the specimen fails by fiber failure. For this test case, a finite element calculation carried out with a classical cohesive interface would not reproduce the interface damage physics correctly. Indeed, in this type of model, the delamination is activated only by out-of-plane stresses which are really small and would not be sufficient to activate the damage mechanism. Experimental damage mechanisms, stress/strain curve, and damage prediction in a cross-ply tensile test OE0=90 4 S with the proposed interface model [27] Another approach, implemented in SAMCEF (the finite element code used here), is the well-known gradient regularization introduced in [61] and developed for laminates in [11]. The scalar "regularized force" Q Y is defined from the initial one Y by the equation:

23.43)
where l c is a characteristic length of the order of the ply thickness. Practically, the parameters l c and c are identified to fit crack initiation instabilities and they define a certain volume to fracture. Illustrations are given in Figs. 23.16 and 23.17 showing a classical OE0 2 =90 2 S tension test for a typical carbon-epoxy composite [59].

Split Detection and Propagation
The damage mesomodel could be too much conservative in some cases such as composite structures involving extensive splitting. Here, the mesomodel has been found to present a shortcoming: due to numerical constants, the calculated splits, i.e., the completely damaged zones due to matrix microcracking, are too thick. To model and simulate splits is a big problem. A first approach is the high fidelity micromechanic model described in Sect. 23.2 where every single discontinuity in the plies and in the interfaces is described and therefore, one has to solve numerically a very large scale problem involving thousands of cracks. Even with a high performance computational tools, such a micromechanic model leads to prohibitive computational efforts and, thus, is far from meeting the virtual structural testing requirements. However, fundamentally, a split appears to be as a particular microcrack relatively isolated and mainly due to shear [44]. From purely microscopic and mesoscopic approaches, intermediate approaches have recently been proposed in the literature. Classical cohesive interfaces are a priori introduced for transverse microcracking and combined with delamination interfaces. Two ways can be distinguished. The first one uses a priori experimental information about the cracking pattern (e.g., the position of the splits) and then is not a predictive approach [62]. The second one considers a set of equally space potential cohesive interfaces in the different plies. To avoid prohibitive computations, the distance between two cohesive interfaces cannot be too small and then is not necessarily compatible with the micromechanics [37,63]. Anyway these approaches do not give the "continuous" part of the microcracking density.
Our understanding is quite different. The difficulty to get "thin" splits with the mesomodel is only due to numerical problems associated with localization phenomena involving very small time and space discretization. The mesomodel, where microcracking is described by microcracking density, is correct and a split

Limits
A first limitation of the proposed mesomodel is that material fracture is described by means of only two types of macrocracks: • delamination cracks within the interfaces; • cracks orthogonal to the laminate's mid-plane, each cracked layer being completely cracked through its thickness.
The layers-in our sense-are assumed not to be too thick. Another limitation is that very severe dynamic loadings cannot be studied as the dynamic wavelength must be larger than the thickness of the plies. Another point, which is not completely satisfying, is the today modeling of fiber fracture initiation in the cases where one has instabilities. Today, one identifies from tests the localization parameters but they are not true material parameters; another way is to use one of the different average techniques developed in [10,72,73]. A better model resulting from an analysis at the fiber scale is at the present time in progress.

Conclusion
The multiscale material model described in this paper for laminate composites is quite general and should be extended to other composite materials. The version of the mesomodel presented here is pending implementation into industrial codes for Virtual Testing. It constitutes a significant improvement over the standard version developed over some 20 years and available in many industrial structural analysis programs, especially regarding the interaction between microcracks and delamination and split modeling and computation.
Taking into account all Lack-of-Knowledge (LoK) sources, essentially here the defects, constitutes one of the scientific challenges involved by Virtual Testing [74]. Not only does one have to simulate reality, but one must also position the result of the model quantitatively with respect to that reality which, itself, presents some variability, both at the micro-, meso-, and macroscales. Moreover, computational costs being prohibitive for designers, another challenge-probably the main onehas to be tackled using the damage mesomodel presented here: the building of virtual charts, i.e., reduced models including the description of uncertainties [75].

Appendix: The Basic Damage Law
Let us introduce z D .1 d/x with d D`.x/. Function`is plotted in Fig. 23.20 and is defined as follows: • for x 2 OE0;˛0, • for x 2 OE˛0; 1, z is linear and d D`.x/ D 1 1 x 1 x 1 ˛0 z.˛0/ So, function`over OE0; 1 is maximum for x D˛. The first derivative of`is continue over the complete interval OE0; 1. z can be interpreted as a stress and x as a strain. Let us note that the slope of the secant is: 1 d. Another main property is related to the maximum of`over OE0; 1. This point can also be interpreted as an instability point or an initiation point. Let us consider that the strain x is the square root of the damage force; one has at the point: