An Energy-Based Concept for Yielding of Multidirectional FRP Composite Structures Using a Mesoscale Lamina Damage Model

Composite structures are made of multidirectional (MD) fiber-reinforced polymer (FRP) composite laminates, which fail due to multiple damages in matrix, interface, and fiber constituents at different scales. The yield point of a unidirectional FRP composite is assumed as the lamina strength limit representing the damage initiation phenomena, while yielding of MD composites in structural applications are not quantified due to the complexity of the sequence of damage evolutions in different laminas dependent on their angle and specification. This paper proposes a new method to identify the yield point of MD composite structures based on the evolution of the damage dissipation energy (DDE). Such a characteristic evolution curve is computed using a validated finite element model with a mesoscale damage-based constitutive model that accounts for different matrix and fiber failure modes in angle lamina. The yield point of composite structures is identified to correspond to a 5% increase in the initial slope of the DDE evolution curve. The yield points of three antisymmetric MD FRP composite structures under flexural loading conditions are established based on Hashin unidirectional (UD) criteria and the energy-based criterion. It is shown that the new energy concept provides a significantly larger safe limit of yield for MD composite structures compared to UD criteria, in which the accumulation of energy dissipated due to all damage modes is less than 5% of the fracture energy required for the structural rupture.


Introduction
The ever-increasing use of advanced polymer composites such as fiber-reinforced polymer (FRP) composites, nanocomposites, etc., as high strength-to-weight ratio and high-stiffness structural materials in advanced industrial applications [1][2][3], presents a unique design challenge with highly anisotropic material responses of the composites [4][5][6][7][8][9]. The FRP composites are essentially brittle with no plastic deformation, thus adequately described using a bilinear stress-strain curves with elastic and the early stage of the loading [10,43,45,51]. In some cases, the MD FRP composite structure was able to sustain up to 10-fold higher loading above the level corresponding to the onset of matrix and interface damages [10,43,45]. Accordingly, considering the UD lamina level yield criteria to estimate the failure of an MD composite structure normally results in the prediction of structural failure at 5%-10% maximum load capacity of the structure [2,10,43,44,52], which prevents the optimum design of light composite structures. Therefore, material and loading-related parameters should be developed to provide a consistent identification of the yield point for the MD FRP composite structures [2,5,10,[43][44][45].
In this study, an energy concept is developed based on the energy dissipated during the inelastic deformation process of lamina using a mesoscale damage model, to estimate the yield of MD FRP composite structures. The critical level of the accumulated damage dissipation energy (DDE) in FRP composite laminate is proposed as the parameter that indicates the yield of the material. The continuum damage mechanics that account for damage initiation and propagation in the material point are used to quantify the material softening process and compute the rate of DDE growth, to establish the critical DDE level. The characteristic evolution of the DDE is established through FE simulations of actual tests using a validated FE model. The yield point is inferred from the DDE curve, when a sudden increasing rate is observed. The approach is illustrated for different types of antisymmetric MD FRP composite structures with the objective of determining the yield strength. The experiments were implemented on carbon and glass fiber-reinforced polymer (CFRP and GFRP) composite structures, such that only lamina damage could occur with negligible interface delamination; therefore, interlaminar damage was not considered [10,43]. A new approach with FE model-based configurations called single-and multi-layer models [10] was used to simulate the FRP composites manufactured using different methods, to model the mesoscale inter-and intralaminar constructions of the composites. The simulation results were validated with the structural response of the composites in the experiments. The method is recommended for determining the yield limit of any type of MD composite structure under different types of load.

Damage Model of FRP Composite Lamina
The response of the FRP composite laminates to applied load such that yielding of the material is achieved is predicted in this study, based on the damage mechanics approach. The damage model of the FRP composite lamina is described below.
The uniaxial behavior of UD FRP composite lamina in orthogonal axes (1-2 axis, Figure 1a) for elastic-damage behavior under tension and compression is shown in Figure 1b. The four bilinear elastic softening curves represent the equivalent stress-displacement behavior of composite lamina in different failure modes of matrix cracking and crushing, and fiber breakage and buckling (following the load arrows of UD lamina, Figure 1b). In an angle lamina under global loading (x-y axis, Figure 1a), the global deformations are mapped into local deformation and used to compute the effective stress parameters. The elastic behavior of the lamina is computed following the classical theory of lamina [53,54]. The stress-displacement relationship of each damage mode is defined in the sections below.

Damage Initiation
The initiation of damage in the lamina for the different failure modes is estimated using Hashin's quadratic stress-based failure model [24]. The model is expressed as a quadratic function of the ratio of the effective stress to strength terms to calculate the values of damage variables for the respective failure mode.
Matrix cracking and crushing: Fiber fracture and buckling/kinking: In the equations above, represents the effective stresses in the lamina, and X T , Y T , X C , Y C , S L , and S T are the strength properties. In Equations (1)-(4), , and , are the internal damage variables in the fiber and matrix phases of the lamina, under tension or compression loadings. Since no plastic deformation is observed in the FRP composite [2,10,55], the permanent deformation of the lamina is considered in the damage evolution processes.

Damage Initiation
The initiation of damage in the lamina for the different failure modes is estimated using Hashin's quadratic stress-based failure model [24]. The model is expressed as a quadratic function of the ratio of the effective stress to strength terms to calculate the values of damage variables for the respective failure mode.
Matrix cracking and crushing: Fiber fracture and buckling/kinking: In the equations above, [σ] represents the effective stresses in the lamina, and X T , Y T , X C , Y C , S L , and S T are the strength properties. In Equations (1) no plastic deformation is observed in the FRP composite [2,10,55], the permanent deformation of the lamina is considered in the damage evolution processes.
Post-damage initiation: Once the onset of damage is predicted in one of the modes, the properties of the material reduce in the other directions/modes and result in early damage. Such effects could be captured by updating the elastic stress tensor (effective stresses in Equations (1)-(4) through internal damage variables. σ ij = σ ij Prior to any damage initiation D σ ij I f any o f the f our damages has initiated (5) where σ is the stress computed using classical lamina theory, [σ] is the effective stress in Equations (1)-(4), and the damage operator, D, is used to consider the effect of early damage initiations. The hypothesis of strain equivalence is used to derive the damage operator as follows [54,56]: where d f , d m , and d s are the fiber, matrix, and shear internal damage variables corresponding to the lamina damage modes in Equations (1)-(4).

Damage Propagation
The evolution of damage to failure at a local material point is obtained through the softening process using energy-based criterion [10,43,45]. In this process, the damage dissipation energy, G DDE (Figure 2a), is employed to determine the constitutive model of the material in each failure mode, which is expressed as the stress-displacement relation. The fracture energy, G C , is the energy that, if dissipated fully, causes the material to fail (G XT C , G XC C , G YT C , and G YC C are the fracture energies in different failure modes; Figure 1b). The value of dissipated energy due to damage is obtained using with the corresponding damage evolution variable, d p , defined as where k 0 eq is the equivalent elastic stiffness, δ 0 eq is the equivalent displacement at the onset of damage in the respective mode d p = 0 , and δ f eq is the equivalent displacement at the separation of the material point d p = 1 . In each failure mode, the critical value of equivalent dissipation energy, G C , is considered as the fracture energy of the lamina. The evolutions of the damage initiation variable (d i in Equations (1)-(4)) and damage propagation variable d p are shown in Figure 2b.
The relation between the equivalent stress-displacement for each failure mode, after onset of damage (dotted lines in Figure 1b,c), is expressed by the equations below [10,[43][44][45].
Matrix tension (σ 22 ≥ 0): Polymers 2020, 12, 157 6 of 18 Matrix compression (σ 22 < 0): Fiber tension (σ 11 ≥ 0): Fiber compression (σ 11 < 0): In these equations, L c is the element characteristic length with magnitude depending on the geometry and the element formulation. For the first-order element, L c is considered as the length of a line across the element. The terms G XT C , G XC C , G YT C , and G YC C are the fiber and matrix fracture energy parameters of the lamina under tension and compression loadings. In Equations (9)-(12), σ o ij , τ o ij , and ε o ij indicate the effective stresses at the onset of damage [10,43].
Polymers 2020, 11, x FOR PEER REVIEW 6 of 18 Fiber compression (σ < 0) : In these equations, is the element characteristic length with magnitude depending on the geometry and the element formulation. For the first-order element, is considered as the length of a line across the element. The terms , , , and are the fiber and matrix fracture energy parameters of the lamina under tension and compression loadings. In Equations (9)-(12), σ , τ , and ε indicate the effective stresses at the onset of damage [10,43].

Damage Dissipation Energy
The energy stored in the FRP composite laminate through elastic-damage deformation, commonly called the internal energy, can be employed to describe the progressive damage process of the composite structure [10,43,45,57]. The internal energy, , can be written for non-viscous composites as = : .
where σ is the stress derived from the constitutive equation of a lamina. The strain rate term is decomposed as where , , are the time rates of elastic, plastic, and creep strains, respectively. Since FRP lamina behaves in the form of elastic-brittle material, thus = = 0, and Equation (13) can be simplified to where is the elastic strain energy. The elastic strain is not recoverable when damage initiates in a

Damage Dissipation Energy
The energy stored in the FRP composite laminate through elastic-damage deformation, commonly called the internal energy, can be employed to describe the progressive damage process of the composite structure [10,43,45,57]. The internal energy, E U , can be written for non-viscous composites as where σ c is the stress derived from the constitutive equation of a lamina. The strain rate term is decomposed as ε cr = 0, and Equation (13) can be simplified to where E S is the elastic strain energy. The elastic strain is not recoverable when damage initiates in a material point. Hence, σ c can be expressed in the following form: where σ u is un-damaged stress, and d is the continuum damage parameter which varies from "zero" for the undamaged state to "one" for the fully damaged state of the material point in the composite lamina. Therefore, substituting σ c into Equation (15) gives elastic strain energy as It is assumed that the damage parameter remains fixed at time t until unloading. Thus, the recoverable strain energy and the dissipated energy during damage can be expressed as follows: Considering the undamaged elastic energy function f u , interchanging the integrals in Equations (18) and (19) yields Since, at time t, d = d t and, at time zero, f u = 0, the first term of the last expression in Equation (21) is zero. By defining the damage strain energy function f c as (1 − d t ) f u , Equations (20) and (21) can be written as follows: The parameter f c can be written for a linear elastic energy function as Substituting Equation (24) into Equations (22) and (23) gets Polymers 2020, 12, 157 The present study employed these equations through FE simulations to establish the characteristic evolution of the DDE and to illustrate the proposed concept for determining the yield limit of MD FRP composite structures [43,45]. The yield point was identified corresponding to a 5% increase in the initial slope of the total DDE evolution curve of the composite structure under specific loading conditions.

Materials and Experimental Procedures
Three different types of antisymmetric MD FRP composite laminate panels were manufactured and machined into beam specimen geometries for flexural tests. These tests were performed in accordance with the ASTM standard [58]. The antisymmetric MD composite lay-up, sample geometry, and load and boundary conditions were selected such that, under the flexural deformation, the various lamina failure modes under tension and compression were activated, while the interlaminar damage and delamination phenomenon were minimized. Subsequently, microscopic fractographic analysis was used to examine the lamina and interface damage events, which indicate dominant lamina damage with negligible interface delamination [10,43,45]. The first group of specimens was fabricated from a glass fiber-reinforced polymer (GFRP) composite panel with thermoplastic resin and eight layers of UD E-glass fiber mats. The GFRP composite was prepared using vacuum-assisted infusion molding (VAIM) process, which resulted in the formation of laminate with no interface between the laminas [10]. The next two groups of the composite samples were made of carbon fiber-reinforced polymer (CFRP) composite laminate. A panel of the CFRP composite was prepared by pre-impregnation of the UD CFRP lamina (M40J fibers and NCHM 6376 resin, Structil France) and cured in an autoclave, resulting in a composite laminate with interface laminas [10]. Details of the manufacturing process of the FRP panels were provided elsewhere [2,10,43]. Microscopic images of the longitudinal cross-section of the FRP composite specimens are shown in Figure 3a. A schematic view of the composite beam in the test set-up and boundary conditions is provided in Figure 3b, in which the lay-ups and dimensions of the beams, and the load configurations are described in Table 1. Several GFRP composite specimens and the first batch of the CFRP composite samples were tested under three-point bending (3PB) conditions, while the second batch of CFRP composite samples was tested under four-point bending (4PB), as mentioned in Table 1. A continuous flexural load was applied to the specimens until a significant amount of degradation in load-deflection response, which represented the occurrence of multiple damages, was observed. The results of the tests in terms of monotonic reaction force versus the deflection of the composite beams were recorded and utilized in the validation procedure of the FE models and simulation processes.
Polymers 2020, 11, x FOR PEER REVIEW 8 of 18 using vacuum-assisted infusion molding (VAIM) process, which resulted in the formation of laminate with no interface between the laminas [10]. The next two groups of the composite samples were made of carbon fiber-reinforced polymer (CFRP) composite laminate. A panel of the CFRP composite was prepared by pre-impregnation of the UD CFRP lamina (M40J fibers and NCHM 6376 resin, Structil France) and cured in an autoclave, resulting in a composite laminate with interface laminas [10]. Details of the manufacturing process of the FRP panels were provided elsewhere [2,10,43]. Microscopic images of the longitudinal cross-section of the FRP composite specimens are shown in Figure 3a. A schematic view of the composite beam in the test set-up and boundary conditions is provided in Figure 3b, in which the lay-ups and dimensions of the beams, and the load configurations are described in Table 1. Several GFRP composite specimens and the first batch of the CFRP composite samples were tested under three-point bending (3PB) conditions, while the second batch of CFRP composite samples was tested under four-point bending (4PB), as mentioned in Table  1. A continuous flexural load was applied to the specimens until a significant amount of degradation in load-deflection response, which represented the occurrence of multiple damages, was observed.
The results of the tests in terms of monotonic reaction force versus the deflection of the composite beams were recorded and utilized in the validation procedure of the FE models and simulation processes.

Finite Element Simulation
Since the manufacturing process of the FRP composite laminates dictates the resulting interface condition between the laminas, the FE models of these cases are different. The GFRP composite laminate, produced by the VAIM method with no apparent interface, was simulated with the single-layer model. The CFRP composite laminate, fabricated via lay-ups of prepreg laminas and autoclave curing with distinct interfaces, was modeled using the multi-layer model. The single-layer and multi-layer FE model constructions are illustrated in Figure 4, while details of these mesoscale FE models were discussed elsewhere [10,43,45].
Polymers 2020, 11, x FOR PEER REVIEW 9 of 18 laminate, produced by the VAIM method with no apparent interface, was simulated with the singlelayer model. The CFRP composite laminate, fabricated via lay-ups of prepreg laminas and autoclave curing with distinct interfaces, was modeled using the multi-layer model. The single-layer and multilayer FE model constructions are illustrated in Figure 4, while details of these mesoscale FE models were discussed elsewhere [10,43,45]. The three-dimensional (3D) FE model geometry of the CFRP composite laminate beam specimen is shown in Figure 5. The damage model (Section 2) was applied in each lamina using the standard model definition step in ABAQUS software [57]. Each lamina was discretized into a layer of eightnode continuum shell elements (SC8R) with reduced integration points for efficient computation [37]. The element mesh was refined with smaller-size elements, defined for the central region of the specimen when the maximum deflection and, thus, damage evolution was anticipated. The loading and support rollers were modeled as rigid bodies and discretized using rigid, four-node continuum elements (R3D4) [57]. The three-dimensional (3D) FE model geometry of the CFRP composite laminate beam specimen is shown in Figure 5. The damage model (Section 2) was applied in each lamina using the standard model definition step in ABAQUS software [57]. Each lamina was discretized into a layer of eight-node continuum shell elements (SC8R) with reduced integration points for efficient computation [37]. The element mesh was refined with smaller-size elements, defined for the central region of the specimen when the maximum deflection and, thus, damage evolution was anticipated. The loading and support rollers were modeled as rigid bodies and discretized using rigid, four-node continuum elements (R3D4) [57].
Frictionless contact was assumed between all contacting bodies. In the multi-layer model (Figure 4b), the interfaces between adjacent laminas were modeled using a surface-to-surface tie with finite displacement interaction of the sharing node pair. This allows relative displacement between adjacent surfaces of the laminas [10,45]. A two-step mesh convergence process was performed to eliminate the effect of element size on the FE-calculated results for the elastic and damage calculations [43,45]. A finer element mesh size than that adequately identified in the elastic analysis is required for element size-independent damage calculations. The resulting element size, at mesh convergence, had an edge length of 0.2 mm. The boundary conditions of the model are illustrated in Figure 3b, while the loading was identical to that used during the test.
The elastic and strength properties of GFRP and CFRP composite laminates used in the FE simulations were obtained through standard tests (ASTM-D4762, [59]) [10,43,45], while the values of fracture energies were extracted from the properties of similar materials in the literature [60][61][62], as listed in Table 2. These properties were utilized rigorously in the FE simulation exercises of the various composite specimen geometries and loading cases, demonstrating comparable load-displacement results with measured responses [2,10,[43][44][45]52].
The three-dimensional (3D) FE model geometry of the CFRP composite laminate beam specimen is shown in Figure 5. The damage model (Section 2) was applied in each lamina using the standard model definition step in ABAQUS software [57]. Each lamina was discretized into a layer of eightnode continuum shell elements (SC8R) with reduced integration points for efficient computation [37]. The element mesh was refined with smaller-size elements, defined for the central region of the specimen when the maximum deflection and, thus, damage evolution was anticipated. The loading and support rollers were modeled as rigid bodies and discretized using rigid, four-node continuum elements (R3D4) [57]. Frictionless contact was assumed between all contacting bodies. In the multi-layer model ( Figure  4b), the interfaces between adjacent laminas were modeled using a surface-to-surface tie with finite displacement interaction of the sharing node pair. This allows relative displacement between adjacent surfaces of the laminas [10,45]. A two-step mesh convergence process was performed to eliminate the effect of element size on the FE-calculated results for the elastic and damage calculations [43,45]. A finer element mesh size than that adequately identified in the elastic analysis is required for element size-independent damage calculations. The resulting element size, at mesh convergence, had an edge

Results and Discussion
The FE-calculated and measured flexural responses of the MD FRP composite beam specimens, expressed in terms of the load-deflection curves, were compared. A comparable response was indicative of the validity of the FE simulation procedures. Subsequently, the calculated deformation and damage responses of the composite laminates were interpreted for the respective failure mechanisms. The characteristic evolution of the DDE could then be established and inferred for the onset of yield of the composite structure.

Structural Response and Damage Evolution of GFRP Composite Beam under Three-Point Bending
The FE-calculated load-deflection response of the MD GFRP composite beam specimen (Case 1, Table 1) was compared with the measured curve, as shown in Figure 6a. The reasonably good prediction of the flexural response rendered the FE simulation valid. The GFRP composite beam showed an initial linear response up to the deflection of about 12 mm, suggesting structural elastic behavior; however, the flexural stiffness of the beam (Figure 6a) indicated a 1.5% reduction of the stiffness compared to initial condition. This was followed by a slight deviation with lower flexural stiffness to the maximum load, likely attributed to damage and possible softening of the structure. The corresponding evolution of the calculated strain energy (SE) and DDE with increasing beam deflection is shown in Figure 6b. A sudden load drop was observed at the maximum load, due to the composite rupture by fiber buckling in the first lamina under compression, as shown in Figure 6c. composite structure, yield began when the rate of the DDE increased to 0.914 N/(mm • s) at a 5% increase in the initial slope DDE evolution curve. Thus, this load (stress) and deflection level of (350 N, 13.5 mm) at which the rate of the DDE abruptly increased was identified as the yield point of the GFRP composite laminates. Deformation beyond the yield limit to failure was dominated by the softening of the structure. The total DDE of the GFRP composite corresponding to the maximum load at failure was 234 N/mm, which was 3.4% of the total SE of the structure at 6960 N/mm (Figure 6b). The flexural stiffness of the GFRP composite was reduced to 3% and 15.3% from the initial value, at the yield point and maximum load level, respectively.

Structural Response and Damage Evolution of CFRP Composite Beam under Three-Point Bending Load
The FE-calculated load-deflection response of the MD CFRP composite beam specimen (Case 2, Table 1) was compared with the measured curve, as shown in Figure 7a. A reasonably good prediction of the flexural response was claimed; thus, the validity of the FE simulation was ensured. Linear elastic response was observed up to a displacement of about 12.5 mm. A noticeable difference between the FE-predicted and measured deformation at larger deflection, with apparent minute load drops predicted along the load-displacement curve. Such a load drop was artificially induced by the relative slip between the CFRP composite laminate beam specimen and the support rollers under the assumed friction-free condition [45]. As described for the GFRP composite (Case 1) above, the observed reduction in the flexural modulus as reflected in the stiffness curve (Figure 7a) was due to the accumulated damage by the multiple modes of the failure of the composite constituents.
The corresponding evolution of the calculated DDE with increasing beam deflection, as shown in Figure 7b, exhibited identical characteristics to that of the GFRP composite (Case 1). Based on the proposed method of specifying the yield point of the composite structure, yield began when the rate of the DDE abruptly increased to 2.1 N/(mm • s) (i.e., 5% increase in the initial slope). Thus, the yield point was identified to correspond to the load and deflection levels of (150 N, 9 mm), when a sudden increase in the rate of the DDE was observed. It is worth mentioning that catastrophic fracture did not occur for this MD CFRP composite beam specimen at the maximum prescribed central displacement of 28 mm. However, extensive softening of the material by the various damage mechanisms was expected to have occurred. The FE model predicted various matrix cracking and crushing phenomena in different laminas with respect to the level of beam deflection, as shown in Results show that a negligible amount of energy was dissipated prior to the beam deflection of 13.5 mm. This suggests that limited damage took place on the structure up to this critical load level. However, the rate of energy dissipation abruptly increased for larger composite beam deflection. Although an MD GFRP composite was considered under complex bending load, only the single damage mode for fiber failure occurred in the first lamina (0 • ) under compression, while other laminas remained elastic until the maximum load was experienced. The fiber damage initiation was predicted at 13.5 mm deflection. Based on the proposed method of specifying the yield point of the composite structure, yield began when the rate of the DDE increased to 0.914 N/(mm·s) at a 5% increase in the initial slope DDE evolution curve. Thus, this load (stress) and deflection level of (350 N, 13.5 mm) at which the rate of the DDE abruptly increased was identified as the yield point of the GFRP composite laminates. Deformation beyond the yield limit to failure was dominated by the softening of the structure. The total DDE of the GFRP composite corresponding to the maximum load at failure was 234 N/mm, which was 3.4% of the total SE of the structure at 6960 N/mm (Figure 6b). The flexural stiffness of the GFRP composite was reduced to 3% and 15.3% from the initial value, at the yield point and maximum load level, respectively.

Structural Response and Damage Evolution of CFRP Composite Beam under Three-Point Bending Load
The FE-calculated load-deflection response of the MD CFRP composite beam specimen (Case 2, Table 1) was compared with the measured curve, as shown in Figure 7a. A reasonably good prediction of the flexural response was claimed; thus, the validity of the FE simulation was ensured. Linear elastic response was observed up to a displacement of about 12.5 mm. A noticeable difference between the FE-predicted and measured deformation at larger deflection, with apparent minute load drops predicted along the load-displacement curve. Such a load drop was artificially induced by the relative slip between the CFRP composite laminate beam specimen and the support rollers under the assumed friction-free condition [45]. As described for the GFRP composite (Case 1) above, the observed reduction in the flexural modulus as reflected in the stiffness curve (Figure 7a) was due to the accumulated damage by the multiple modes of the failure of the composite constituents.
Polymers 2020, 11, x FOR PEER REVIEW 12 of 18 Figure 7c. The time to the onset of matrix cracking was predicted firstly in the bottom most lamina (−45°) of the composite beam when loaded up to 4.7 mm deflection, which represented the elastic deformation limit of the composite structure. The accumulated DDE at this displacement was 680 N/mm, which was 14% of the total strain energy of the composite structure (4825 N/mm, Figure 7b). The initial flexural stiffness of the CFRP composite beam of 17.6 N/mm was reduced to 17.44 and 10.66 N/mm (0.91% and 39.4% reduction) at the yield point and maximum load level, respectively.

Structural Response and Damage Evolution of CFRP Composite Beam under Four-Point Bending
The FE-calculated load-deflection response of the MD GFRP composite beam specimen (Case 3, Table 1) under 4PB was compared with the measured curve, as shown in Figure 8a. Again, a reasonably good comparison was demonstrated, and a valid FE simulation procedure of the test was claimed. A similar trend of the flexural deformation of the MD FRP composite beam specimens was observed in the three different specimens considered in this study. A larger applied total load of 600 N was recorded over a much shorter prescribed displacement of 8 mm, when compared with the 3PB (Case 2, Table 1) and different anti-symmetric lay-ups of the specimen. Linear elastic flexural response was measured up to the central deflection of about 4.6 mm, while the flexural stiffness curve ( Figure  8a) was reduced by 7.1% from the initial condition. The corresponding evolution of SE and DDE in this antisymmetric CFRP composite beam specimen, with the applied displacement of the loading rollers, under the 4PB load, is shown in Figure 8b. Based on the proposed method of identifying the yield point of the composite, yield commenced when the rate of the DDE sharply rises to 11.1 N/(mm • s). The corresponding load and deflection level at yield were 313 N and 3 mm, respectively.
It is worth noting that the load-deflection curve remained fairly linear beyond the yield point up to about 4.7 mm. Within this small deflection range, the rapid evolution of matrix damage in all of the laminas under tensile deformation (Figure 8c) only contributed to a fraction (about 10%) of the total DDE; thus, the effect on the softening of the material remained insignificant. A similar observation applied to the previous different FRP composite specimens and loading conditions, as described above. The accumulated DDE at the end of the prescribed displacement was 598.4 N/mm, which was 23% of the total SE of the composite beam (2567 N/mm, Figure 8b). The flexural stiffness of the CFRP The corresponding evolution of the calculated DDE with increasing beam deflection, as shown in Figure 7b, exhibited identical characteristics to that of the GFRP composite (Case 1). Based on the proposed method of specifying the yield point of the composite structure, yield began when the rate of the DDE abruptly increased to 2.1 N/(mm·s) (i.e., 5% increase in the initial slope). Thus, the yield point was identified to correspond to the load and deflection levels of (150 N, 9 mm), when a sudden increase in the rate of the DDE was observed. It is worth mentioning that catastrophic fracture did not occur for this MD CFRP composite beam specimen at the maximum prescribed central displacement of 28 mm. However, extensive softening of the material by the various damage mechanisms was expected to have occurred. The FE model predicted various matrix cracking and crushing phenomena in different laminas with respect to the level of beam deflection, as shown in Figure 7c. The time to the onset of matrix cracking was predicted firstly in the bottom most lamina (−45 • ) of the composite beam when loaded up to 4.7 mm deflection, which represented the elastic deformation limit of the composite structure. The accumulated DDE at this displacement was 680 N/mm, which was 14% of the total strain energy of the composite structure (4825 N/mm, Figure 7b). The initial flexural stiffness of the CFRP composite beam of 17.6 N/mm was reduced to 17.44 and 10.66 N/mm (0.91% and 39.4% reduction) at the yield point and maximum load level, respectively.

Structural Response and Damage Evolution of CFRP Composite Beam under Four-Point Bending
The FE-calculated load-deflection response of the MD GFRP composite beam specimen (Case 3, Table 1) under 4PB was compared with the measured curve, as shown in Figure 8a. Again, a reasonably good comparison was demonstrated, and a valid FE simulation procedure of the test was claimed. A similar trend of the flexural deformation of the MD FRP composite beam specimens was observed in the three different specimens considered in this study. A larger applied total load of 600 N was recorded over a much shorter prescribed displacement of 8 mm, when compared with the 3PB (Case 2, Table 1) and different anti-symmetric lay-ups of the specimen. Linear elastic flexural response was measured up to the central deflection of about 4.6 mm, while the flexural stiffness curve (Figure 8a) was reduced by 7.1% from the initial condition. The corresponding evolution of SE and DDE in this antisymmetric CFRP composite beam specimen, with the applied displacement of the loading rollers, under the 4PB load, is shown in Figure 8b. Based on the proposed method of identifying the yield point of the composite, yield commenced when the rate of the DDE sharply rises to 11.1 N/(mm·s). The corresponding load and deflection level at yield were 313 N and 3 mm, respectively. It is worth noting that the load-deflection curve remained fairly linear beyond the yield point up to about 4.7 mm. Within this small deflection range, the rapid evolution of matrix damage in all of the laminas under tensile deformation (Figure 8c) only contributed to a fraction (about 10%) of the total DDE; thus, the effect on the softening of the material remained insignificant. A similar observation applied to the previous different FRP composite specimens and loading conditions, as described above. The accumulated DDE at the end of the prescribed displacement was 598.4 N/mm, which was 23% of the total SE of the composite beam (2567 N/mm, Figure 8b). The flexural stiffness of the CFRP composite was reduced by 1.1% and 31% of the initial value at the yield point and maximum load level (8-mm deflection), respectively.

Comparison of the Estimated Yield Limits Based on UD Hashin Criteria and Energy-Based Criterion
Three cases of FRP composite structure were tested under flexural loading condition, while the FE model was used to predict the material behavior and specify the structural yielding based on UD criteria, as well as the energy-based model proposed in this study. A summary of the results in terms of the maximum load and deflection capacities of the structures, yielding according to the UD criteria and energy-based criteria, and the percentage of the yield value to the maximum capacity (MC) of the structure is listed in Table 3.

Comparison of the Estimated Yield Limits Based on UD Hashin Criteria and Energy-Based Criterion
Three cases of FRP composite structure were tested under flexural loading condition, while the FE model was used to predict the material behavior and specify the structural yielding based on UD criteria, as well as the energy-based model proposed in this study. A summary of the results in terms of the maximum load and deflection capacities of the structures, yielding according to the UD criteria and energy-based criteria, and the percentage of the yield value to the maximum capacity (MC) of the structure is listed in Table 3. Table 3. Results of the yield values (UD and energy-based criteria) of the MD composite structure under flexural loading condition.

Comparison of the Estimated Yield Limits Based on UD Hashin Criteria and Energy-Based Criterion
Three cases of FRP composite structure were tested under flexural loading condition, while the FE model was used to predict the material behavior and specify the structural yielding based on UD criteria, as well as the energy-based model proposed in this study. A summary of the results in terms of the maximum load and deflection capacities of the structures, yielding according to the UD criteria and energy-based criteria, and the percentage of the yield value to the maximum capacity (MC) of the structure is listed in Table 3. Results indicated that, in the condition of the single failure mode (Case 1), the UD criteria (i.e., Hashin model) and energy-based criterion would suggest a similar range for yielding of the MD composite structure. However, MD composite structures mostly failed due to multiple damage phenomena (Cases 2 and 3), in which case, considering UD criteria would result in the assumption of Results indicated that, in the condition of the single failure mode (Case 1), the UD criteria (i.e., Hashin model) and energy-based criterion would suggest a similar range for yielding of the MD composite structure. However, MD composite structures mostly failed due to multiple damage phenomena (Cases 2 and 3), in which case, considering UD criteria would result in the assumption of structural yielding at only 10%-20% maximum capacity of the structure (displacement or load). While using the energy-based criterion, the yielding limit could safely (controlled condition in which the accumulation of all damage modes contributed to less than 5% of energy dissipation) increase to 30%-50% of the maximum capacity of the structure. The knowledge of the safe limit of structural yielding could be used for the optimum design of composite structures that are typically implemented under complex loading conditions.

Conclusions
A new criterion was proposed to identify the onset of yield in MD FRP composite laminate structures subjected to any type of loading condition. The new energy concept provides a significantly larger safe limit of yield for MD composite structures that normally fail due to multiple damage phenomena, in which the accumulation of energy dissipated due to all damage modes is less than 5% of the fracture energy required for the structural rupture. The criterion is based on the damage energy dissipated through the multiple softening processes of composite laminate materials under different mesoscale failure modes. The characteristic evolution of the DDE was calculated using a validated FE model with damage-based constitutive equations. The criterion was examined for antisymmetric MD GFRP and CFRP composite laminates under three-and four-point bending test configurations. The following can be concluded:

1.
The yield point of the FRP composite laminate structures could be identified by a 5% increase in the initial slope of the DDE evolution curve with respect to the applied load parameters.

2.
At the yield point, the extent of damage by the various modes depended on material, lay-ups, load, and test configurations.

3.
The yield points of the MD GFRP and CFRP composite laminates (cases 1, 2, and 3) were identified to occur upon flexural loading when the rate of the DDE reached 0.914, 2.1, and 11.1 N/(mm·s), respectively. The corresponding deflections were 13.5, 9, and 3 mm, respectively. 4.
The initial flexural stiffness of the MD GFRP and CFRP composite structures (cases 1, 2, and 3) were measured at 28, 17.6, and 108.26 N/mm, reduced to 27.2, 17.44, and 107.1 N/mm at the yield point, indicating 3%, 0.91%, and 1.1% reductions in the stiffness of the beams, respectively. Therefore, an average 2% reduction in flexural stiffness could be suggested as a mean for the determination of the yield point in MD FRP composite structures under three-and four-point bending loads.

5.
In general, the UD criteria resulted in the assumption of structural yielding at 10%-20% maximum capacity of the structure (displacement or load), whereas, using the energy-based criterion, the yield limit could be safely increased to 30%-50% of the maximum capacity of the structure. Interaction between the nodes-in-contact as separation and sliding motions (mm) N i Interpolation function of the interface segments ρ n Interface segment curvature D ij Stiffness matrix of interface with linear coupled elastic behavior (MPa) F i Force of interface node i th (N) u i Motion of interface node i th (mm)