An Anisotropic Damage-Plasticity Constitutive Model of Continuous Fiber-Reinforced Polymers

Accurate structural analyses of continuous fiber-reinforced polymers (FRPs) are imperative for diverse engineering applications, demanding efficient material constitutive models. Nonetheless, the constitutive modeling of FRPs is complicated by the nonlinear behavior resulting from internal damages and the inherent plasticity. Consequently, this study presents an innovative anisotropic constitutive model for FRPs, designed to adeptly capture both the damage evolution and plasticity. All requisite parameters can be easily obtained through fundamental mechanical tests, rendering the model practical and user-friendly. The model utilizes the three-dimensional Puck criteria to determine damages, initiating the evolution process through a combination of continuum damage mechanics and linear stiffness attenuation methods. This evolution is coupled with a one-parameter plastic model. Subsequently, the numerical implementation method, integrated into ANSYS, is detailed. This emphasizes the Cauchy stress and consistent tangent stiffness solution strategy. Finally, the effectiveness of the developed model is demonstrated through comprehensive verification, encompassing existing biaxial tension and open-hole-tension tests conducted on carbon and glass FRP laminates. The simulation results exhibit a remarkable correspondence with the experimental data, validating the reliability and accuracy of the proposed model.


Introduction
Continuous fiber-reinforced polymers (FRPs) are progressively replacing metallic materials owing to their advantages in terms of their light weight, high stiffness, exceptional strength, and superior fatigue performance [1][2][3][4].Realizing the full mechanical potential of FRPs necessitates precise structural analyses, which, in turn, require a robust constitutive modeling approach.However, the complexity of constitutive modeling for FRPs arises from distinctive load-dependent failure modes [5,6] and highly nonlinear deformation [7].Concerning the former, numerous effective failure criteria for FRPs have been proposed [8], and predicting failure in FRPs poses minimal challenges in most engineering applications.The primary challenge in constitutive modeling lies in addressing the latter-the accurate representation of nonlinear deformation characteristics.
One contributing factor to the nonlinear behavior of FRPs is the presence of internal damage [9].This damage initially manifests at the micro-level, occurring at a scale comparable to the fiber diameter.It involves phenomena such as fiber breakage, micro-cracking in the matrix, and localized fiber/matrix de-bonding [10].With the progression of damage, these micro-level effects aggregate into macro-scale manifestations, encompassing matrix splitting, the wedge effect, and the scissoring effect [11].This cumulative damage diminishes the load-carrying capacity of the lamina, ultimately leading to structural failure [12].
Polymers 2024, 16, 334 2 of 18 Collectively, both micro-and macro-damages contribute to the macroscopic nonlinear stress-strain response observed in FRPs.
As for micro-damages, the micromechanical modeling method serves as the most direct approach for describing its effects.While capable of providing a detailed micro-stress distribution, the computational expense associated with this method limits its practical application.Consequently, it is often employed specifically for capturing the elastic properties [9,13] and damage parameters [14,15] of FRPs.In contrast, the continuum damage mechanics (CDM) method offers greater computational efficiency.Numerous CDM models exist, among which the Cachan model stands out as one of the most widely utilized [16][17][18][19].This model, utilizing the Gibbs potential, establishes a concise damage evolution based on energy-release rates [20,21].As for macro-damages, a commonly adopted and convenient approach involves defining a stiffness degradation law.Depending on the specific degradation law applied, stiffness degradation methods can be categorized into linear [22][23][24], bilinear [25,26], and exponential degradation methods [27,28].
Another contributing factor of the nonlinear behavior is attributed to the irreversible deformation of the matrix [29,30].Under matrix-dominated loading conditions or in the case of unidirectional laminates, damage has a minor effect on the nonlinear behavior [29].In such cases, the tangent elastic moduli degrade rapidly at the very beginning (less than 10% for in-plane shear).Interestingly, the in situ modulus at the same stress level shows no apparent change during the unloading process [31][32][33].This nonlinearity cannot be solely attributed to damages, as cracks are scarce and invisible [34,35].These observations highlight the necessity for the constitutive model of FRPs to encompass plasticity.
Taking into account the aforementioned considerations, this research aims to propose an anisotropic constitutive model for FRPs, with the capacity to comprehensively capture micro-damage, macro-damage, and plasticity phenomena.This study elucidates the methodologies for damage determination and evolution and the characterization of plasticity.Additionally, detailed descriptions of the numerical implementation integrated into ANSYS and the model verification process are provided.In summary, the constructed model demonstrates satisfactory precision.All intrinsic parameters can be easily derived through fundamental mechanical tests, enhancing the model's practicality and user-friendliness for composite analysis in engineering applications.To the best of the knowledge of the authors, there have been no reports on constructing the constitutive model of FRPs in this way to date.

Modeling Methodology
In light of the aforementioned characteristics of internal damage, the loading process of an FRP unidirectional lamina was delineated into three stages, as illustrated in Figure 1: 1.
Pure elastoplastic stage (OA segment): in this initial stage, the matrix remains devoid of any damage, and the nonlinearity is exclusively attributed to plasticity.

2.
Micro-damage stage (AB segment): As the loading progresses, micro-damages gradually emerge, introducing a subtle deviation in the elastic stress-strain response.The nonlinear behavior of FRPs contributes to both plasticity and internal micro-damages.

3.
Macro-damage stage (BC segment): Advancing to the final stage, macro-damages emerge, resulting in a noticeable deviation in the elastic stress-strain response.The nonlinear behavior of FRPs contributes to both plasticity and macro-damages.
Accordingly, the constitutive model should comprise three components: (1) Damage determination.This component evaluates the onset of micro-damages and their progression into macro-damages, specifically identifying points A and B.

Damage Determination
The three-dimensional (3D) Puck criteria [5,36], acknowledged as being among the most advanced failure criteria since the worldwide failure exercise (WWFE) [7,8,37,38], were employed to assess damage initiation.Rooted in Mohr's fracture hypotheses, Puck's criteria develop from the concept of intrinsically brittle materials.The criteria encompass both fiber fracture (FF) and inter fiber fracture (IFF).The latter specifically refers to a crack running parallel to the fibers through the layer thickness and is characterized by four distinct forms (Figure 2) [5,36].The stress exposure for FF is [5] Figure 1.Three stages of the loading process of an FRP unidirectional lamina.

Damage Determination
The three-dimensional (3D) Puck criteria [5,36], acknowledged as being among the most advanced failure criteria since the worldwide failure exercise (WWFE) [7,8,37,38], were employed to assess damage initiation.Rooted in Mohr's fracture hypotheses, Puck's criteria develop from the concept of intrinsically brittle materials.The criteria encompass both fiber fracture (FF) and inter fiber fracture (IFF).The latter specifically refers to a crack running parallel to the fibers through the layer thickness and is characterized by four distinct forms (Figure 2) [5,36].

Damage Determination
The three-dimensional (3D) Puck criteria [5,36], acknowledged as being among the most advanced failure criteria since the worldwide failure exercise (WWFE) [7,8,37,38], were employed to assess damage initiation.Rooted in Mohr's fracture hypotheses, Puck's criteria develop from the concept of intrinsically brittle materials.The criteria encompass both fiber fracture (FF) and inter fiber fracture (IFF).The latter specifically refers to a crack running parallel to the fibers through the layer thickness and is characterized by four distinct forms (Figure 2) [5,36].The stress exposure for FF is [5] Figure 2. Forms of IFF (Reproduced with permission from [5]).
The stress exposure for FF is [5] where '1' denotes the fiber direction; X t 11 and X c 11 signify the longitudinal tensile and compression strengths, respectively.The stress exposure for IFF is [5] and where '2' and '3' denote the fiber transverse direction; 'θ' is the inclined angle of the action plane; the superscript 't' denotes tension properties; 'c' signifies compression properties; p t,c ⊥⊥ and p t,c ⊥∥ are material-related inclination parameters; and R ⊥ , R ⊥⊥ , and R ⊥∥ denote the resistance against σ n , τ nt , and τ n1 , respectively.In this paper, R t ⊥ and R ⊥∥ were set as X t,0 22 (representing the strength corresponding to a deviation of 0.5 MPa from the linear zone) and S 0 12 (indicating the strength at which the secant in-plane shear modulus deviates by 7% from the elastic modulus) to ascertain the initiation of micro-damages [33,39], denoted as point A in Figure 1.Additionally, the values X t,u 22 , X c,u 22 , and S u 12 , corresponding to the failure strength of a unidirectional lamina under transverse tension, compression, and in-plane strength, were assigned as R t ⊥ , R c ⊥ , and R ⊥∥ respectively.These serve to assess the onset of macro-damages, identified as point B in Figure 1.

Damage Evolution
The Cachan CDM model was utilized to depict the evolution of micro-damages in segment AB.The damaged strain energy density W D was rewritten in terms of stresses on Puck's action plane: where d 2 and d 6 are micro-damage coefficients reflecting transverse and in-plane shear micro-damages, respectively.A quantity was then defined [21]: where 'sup' indicates that Y is the maximum value of the loading history; Y d and Y d ′ are the partial derivations of W D with respect to d 2 and d 6 , respectively: Polymers 2024, 16, 334 5 of 18 The corresponding evolution laws for micro-damages were then written as [21]    The parameters b, Y 0 , Y c , Y ′ 0 , and Y ′ c are material-related parameters and can be characterized through cyclic tensile tests of (±45 • ) 4s and (±67.5 • ) 4s symmetric laminates.Details of the procedure for determining the parameters can be found in Ref. [21].
Furthermore, a linear degradation method was employed to illustrate the evolution of macro-damages in segment BC.The degradation coefficients η E and η G , which represent the degradation of the transverse and in-plane shear elastic modulus due to macro-damages, respectively, are expressed as The superscript '0' signifies the tangent modulus at the onset of macro-damages.η rE and η rG represent the residual stiffness of the transverse and in-plane shear elastic modulus, respectively [40].Furthermore, f E,IFF,max denotes the maximum allowable stress exposure for IFF.

Plasticity Evolution
As the last component of the constitutive model, an anisotropic elastoplastic model, derived from a one-parameter plastic model proposed by Neto [41], was proposed to capture the plasticity of the matrix.The engineering stress and strain were represented in the vector form: where the bold font indicates that the variable is in the form of a matrix.To enhance clarity, the five fundamental components of plasticity, i.e., the elastoplastic strain decomposition, an elastic law, a yield criterion stating the yield function, a plastic flow rule defining the evolution of the plastic strain, and a hardening law characterizing the evolution of the yield limit, were systematically presented one after another.

(a) Elastoplastic strain decomposition
The elastoplastic strain ε was decomposed into the elastic-strain-coupled damage ε ed and the plastic strain ε p .The corresponding rate form is With the assumption of the strain equivalence, the elastic law considering the material damage is where D e is the Voigt form of the four-order stiffness tensor, and ∼ σ is the effective stress tensor, given by Polymers 2024, 16, 334 where M is the damage matrix defined as (c) Yield criterion The equivalent stress ∼ σ y is defined as [42] where a is a material-related coefficient, reflecting the coupling between the transverse plasticity and in-plane shear plasticity.Then, a quadratic anisotropic yield function Φ was utilized: where P is defined as diag(0,3,3,6a,6a,6a). Φ ≤ 0 denotes an elastic state, and Φ > 0 signifies a plastic state.In addition, the relationship between the equivalent stress and the equivalent plastic strain ε p was expressed as the exponent form: where α and β are material-related coefficients.The transverse and in-plane shear linear region of FRPs are small and the exponent function can be approximated as linear in that region.Therefore, the initial yield stress was set to zero.

(d) Plastic flow rule
An associative plasticity was presumed, and the evolution of the plastic strain is expressed as where γ is Lagrange's plastic multiplier; N is the flow vector and equals After definition of the plastic flow rule, the equivalent plastic strain can then be expressed by ε p under the plastic work equivalence: By substituting Equations ( 18), (20), and ( 21) into ( 22), the following was obtained: .

Numerical Implementation Methodology
The material constitutive equations are implemented in Fortran under the Visual Studio framework and linked to ANSYS 19.0 through UserMat, a user-programmable feature of ANSYS.The role of UserMat is to update the Cauchy stress and consistent tangent stiffness by utilizing the received stress, strain, and state variable values of the element at every material integration point during the solution phase.The following sections detail the Cauchy stress solution, consistent tangent stiffness solution, and total numerical algorithm of the developed constitutive model.

Cauchy Stress Solution
The elastoplastic equations present in Section 2.2 are transformed into the backward Euler discretization form: and ∆γ ≥ 0, Φ σn+1 , where the subscripts 'n' and 'n + 1' denote the calculation step of finite element analysis (FEA).To solve this algebraic system, the return mapping algorithm (RMA) is utilized.At first, it is assumed that the material deformation from step 'n' to step 'n + 1' is elastic, and the strain at step 'n + 1' ought to be ε e,trial n+1 = ε e n + ∆ε (29) and ε where the superscript 'trial' denotes a trial quantity.Then, the yield function is calculated with trial strain.If Φ trial ≤ 0, no plasticity occurs and the true quantity is the same as that of the trial.Otherwise, plasticity occurs and a 'plastic correction' should be executed.The strain at step 'n + 1' ought to be and with Multiplying (31) by D e , we can obtain the explicit form of the effective stress: Then, substituting Equations ( 31), (32), and (34) into (33), a single plastic equation with respect to ∆γ of the 'plastic correction' is obtained [41]: Polymers 2024, 16, 334 This equation can be solved by the Newton-Raphson method.The iterative scheme is where '(k)' denotes the iterative step, and it is initialized by The derivative part of Equation ( 36) is where d σn+1 /d∆γ equals Thus far, the Cauchy stress and equivalent plastic strain (used for iterative initialization of the next FEA step) can be updated iteratively until convergence:

Consistent Tangent Stiffness Solution
In order to obtain the constituent tangent stiffness, we transfer Equations ( 26)-( 28) into their differential form: From this, the relationship between the stress and trial strain is derived as Therefore, the consistent tangent stiffness can be deduced according to the chain derivation rule:

Numerical Algorithm
The numerical algorithm for the developed constitutive model is illustrated in Figure 3.The algorithm comprises two key steps: plasticity evolution and damage evolution.The process initiates with the computation of effective stress and elastic strain, which is followed Polymers 2024, 16, 334 9 of 18 by solving the elastoplastic equations.Subsequently, the damage status is determined based on the 3D Puck criterion, and finally, the damage parameters are updated.In this solving strategy, the damage and elastoplastic variables are subsequently refreshed.The maintenance of satisfactory calculation accuracy could be ensured by judiciously controlling the time step.In this solving strategy, the damage and elastoplastic variables are subsequently refreshed.The maintenance of satisfactory calculation accuracy could be ensured by judiciously controlling the time step.
The test data from the first two cases were a part of the benchmark data of the WWFE, a famous academic event known as the 'Failure Olympics' of FRPs.A concise summary of the specimen preparation, test equipment, loading scheme, and final results is available in [7], with comprehensive details provided in the respective original articles.Detailed experimental information of the third case can be found in [49].
The parameters for E-glass/epoxy, AS4-carbon/epoxy, and T300-carbon/epoxy have been compiled in Table 1.Specifically, the elastic modulus, strength, and stress-strain curves of the unidirectional laminae were referenced from [37,50].Coefficients relevant to damage determination were collated from [5], while damage evolution coefficients were

Materials Tested and Work Method
To comprehensively validate the constructed constitutive model across various lay-up configurations and stress states, a thorough comparison was conducted between simulations and existing test results of FRP laminates.The assessment focused on three perspectives: 1.

3.
Open-hole-tension (OHT) strength prediction.In this case, OHT tests of T300-carbon/ epoxy laminates conducted by Chang [49,50] were chosen as a reference.These tests encompass three different lay-ups ([0/(±45 , each with four different geometric configurations.The test data from the first two cases were a part of the benchmark data of the WWFE, a famous academic event known as the 'Failure Olympics' of FRPs.A concise summary of the specimen preparation, test equipment, loading scheme, and final results is available in [7], with comprehensive details provided in the respective original articles.Detailed experimental information of the third case can be found in [49]. The parameters for E-glass/epoxy, AS4-carbon/epoxy, and T300-carbon/epoxy have been compiled in Table 1.Specifically, the elastic modulus, strength, and stress-strain curves of the unidirectional laminae were referenced from [37,50].Coefficients relevant to damage determination were collated from [5], while damage evolution coefficients were extracted from [21,40].The plastic evolution coefficient a was adopted as recommended in [51].The current emphasis was directed towards the determination of the values associated with the plastic evolution coefficients, denoted as α and β.It is widely recognized that unidirectional laminae exhibit significant plasticity in in-plane shear [29][30][31][32][33][34][35].Drawing upon this recognition, the determination of α and β can be facilitated through successive adjustments, ensuring a close alignment between the simulated in-plane stress-strain curve and the corresponding experimental observations.A quarter FEA model was employed for both biaxial tension and OHT conditions, as illustrated in Figures 4a and 4b, respectively.Symmetry boundaries were constrained in their normal direction.SOLID185, a 3D, eight-node, layered solid element in ANSYS, was utilized for meshing.Layer properties, encompassing the layer orientation angle and thickness, were established through section commands (SECTYPE and SECDATA), with three integration points per layer specified.The previously constructed constitutive model was used as the material model.For the biaxial load condition, the dimensions of the model were set at 100 mm in width and length.The thicknesses of the 90° and 30° E-glass/epoxy layer were 0.172 mm and 0.414 mm [37], respectively, while all layers of the AS4-carbon/epoxy measured 0.1375 mm each [37].The elements were uniformly meshed, and a mesh sensitivity analysis was conducted, revealing minimal mesh dependency in the present numerical model.No stress concentrations were observed, and the stress distribution within each layer was uniform.A mesh density of 50 × 50 × 1 was deemed sufficient and was consequently employed.
For OHT, all layers of the T300-carbon/epoxy possessed the same thickness, and the overall thickness of the laminate was 2.62 mm [50].The laminate geometry (i.e., hole diameter, width, and length, respectively) encompassed four categories: 3.18 mm, 19.05 mm, and 177.8 mm; 6.35 mm, 25.4 mm, and 203.2 mm; 3.18 mm, 12.7 mm, and 177.8 mm; and 6.35 mm, 25.4 mm, and 203.2 mm [49,50].For convenience, these categories are successively labeled A, B, C, and D in the subsequent discussions.Due to the stress concentration near the central hole in the OHT specimen, the meshes within an 8 mm radius of the sample center were refined.Line divisions along the radius and circumference of the quarter circle were set at 20.The longitudinal division of the remaining part of the specimen was established at 30, with a space ratio of 0.6.A mesh sensitivity analysis demonstrated that this mesh density was sufficient.The simulation results were found to be consistent with a model employing twice the mesh density.

Biaxial Tension Failure Envelope
Simulations of biaxial tests were conducted on (90°/+30°/−30°)s E-glass/LY556 laminates.The features of this test are as follows [38]: (1) the biaxial strength under tensioncompression loading is lower than the uniaxial strength values, and (2) there is a tendency for an enhancement in the biaxial tensile strength.As is illustrated in Figure 5, the simulation results effectively capture both features, and the projected strength aligns well with the test results.For the biaxial load condition, the dimensions of the model were set at 100 mm in width and length.The thicknesses of the 90 • and 30 • E-glass/epoxy layer were 0.172 mm and 0.414 mm [37], respectively, while all layers of the AS4-carbon/epoxy measured 0.1375 mm each [37].The elements were uniformly meshed, and a mesh sensitivity analysis was conducted, revealing minimal mesh dependency in the present numerical model.No stress concentrations were observed, and the stress distribution within each layer was uniform.A mesh density of 50 × 50 × 1 was deemed sufficient and was consequently employed.
For OHT, all layers of the T300-carbon/epoxy possessed the same thickness, and the overall thickness of the laminate was 2.62 mm [50].The laminate geometry (i.e., hole diameter, width, and length, respectively) encompassed four categories: 3.18 mm, 19.05 mm, and 177.8 mm; 6.35 mm, 25.4 mm, and 203.2 mm; 3.18 mm, 12.7 mm, and 177.8 mm; and 6.35 mm, 25.4 mm, and 203.2 mm [49,50].For convenience, these categories are successively labeled A, B, C, and D in the subsequent discussions.Due to the stress concentration near the central hole in the OHT specimen, the meshes within an 8 mm radius of the sample center were refined.Line divisions along the radius and circumference of the quarter circle were set at 20.The longitudinal division of the remaining part of the specimen was established at 30, with a space ratio of 0.6.A mesh sensitivity analysis demonstrated that this mesh density was sufficient.The simulation results were found to be consistent with a model employing twice the mesh density.

Biaxial Tension Failure Envelope
Simulations of biaxial tests were conducted on (90 • /+30 • /−30 • ) s E-glass/LY556 laminates.The features of this test are as follows [38]: (1) the biaxial strength under tension-compression loading is lower than the uniaxial strength values, and (2) there is a tendency for an enhancement in the biaxial tensile strength.As is illustrated in Figure 5, the simulation results effectively capture both features, and the projected strength aligns well with the test results.The biaxial test comparison of AS4/3501-6 (90/+45°/−45°/0°)s quasi-isotropic laminates is depicted in Figure 6.The features of this test are as follows [38]: (1) the biaxial strengths are slightly higher than the uniaxial strength values, and (2) once compression stress in the x-direction is applied, the strength in the y-direction is reduced below the uniaxial strength.As is shown in Figure 6, the simulation results effectively represent both features, and the projected strength aligns well with the test results, with the exception of slight non-conservation in the second quadrant.

Biaxial Tension Deformation
Simulations of the biaxial tension under a stress ratio of  : = 1:20 were conducted on the AS4/3501-6 (90/+45°/−45°/0°)s laminate (Figure 7).This biaxial load nearly equaled the y-direction tension.As a result,  exhibited negative results due to the Poisson ratio effect.Figure 7 illustrates that both the simulation curve and experimental test exhibited some nonlinearity.This is attributed to the transverse and shear damages in the 0° and ±45° layers.Simulations also indicated that the 0° layers entered the macro-damage stage at a  of 256 MPa, and the stress-strain response exhibited a slight deviation from the elastic response (the black dashed line in Figure 7).At 352 MPa, the ±45° layers entered the macro-damage stage, and the deviation was more pronounced.The stiffness attenuation of the 0° layers was notable around 540 MPa (with  ,  ,  , and  of the 0° layers reaching 0.7, 0.8, 0.94, and 0.98 at 544 MPa, respectively).The 90° and ±45° layers would carry a greater load, resulting in a gradual change in the curve slope around this stress level.The biaxial test comparison of AS4/3501-6 (90/+45 • /−45 • /0 • ) s quasi-isotropic laminates is depicted in Figure 6.The features of this test are as follows [38]: (1) the biaxial strengths are slightly higher than the uniaxial strength values, and (2) once compression stress in the x-direction is applied, the strength in the y-direction is reduced below the uniaxial strength.As is shown in Figure 6, the simulation results effectively represent both features, and the projected strength aligns well with the test results, with the exception of slight non-conservation in the second quadrant.The biaxial test comparison of AS4/3501-6 (90/+45°/−45°/0°)s quasi-isotropic laminates is depicted in Figure 6.The features of this test are as follows [38]: (1) the biaxial strengths are slightly higher than the uniaxial strength values, and (2) once compression stress in the x-direction is applied, the strength in the y-direction is reduced below the uniaxial strength.As is shown in Figure 6, the simulation results effectively represent both features, and the projected strength aligns well with the test results, with the exception of slight non-conservation in the second quadrant.

Biaxial Tension Deformation
Simulations of the biaxial tension under a stress ratio of  : = 1:20 were conducted on the AS4/3501-6 (90/+45°/−45°/0°)s laminate (Figure 7).This biaxial load nearly equaled the y-direction tension.As a result,  exhibited negative results due to the Poisson ratio effect.Figure 7 illustrates that both the simulation curve and experimental test exhibited some nonlinearity.This is attributed to the transverse and shear damages in the 0° and ±45° layers.Simulations also indicated that the 0° layers entered the macro-damage stage at a  of 256 MPa, and the stress-strain response exhibited a slight deviation from the elastic response (the black dashed line in Figure 7).At 352 MPa, the ±45° layers entered the macro-damage stage, and the deviation was more pronounced.The stiffness attenuation of the 0° layers was notable around 540 MPa (with  ,  ,  , and  of the 0° layers reaching 0.7, 0.8, 0.94, and 0.98 at 544 MPa, respectively).The 90° and ±45° layers would carry a greater load, resulting in a gradual change in the curve slope around this stress level.

Biaxial Tension Deformation
Simulations of the biaxial tension under a stress ratio of σ x :σ y = 1:20 were conducted on the AS4/3501-6 (90/+45 • /−45 • /0 • ) s laminate (Figure 7).This biaxial load nearly equaled the y-direction tension.As a result, ε x exhibited negative results due to the Poisson ratio effect.Figure 7 illustrates that both the simulation curve and experimental test exhibited some nonlinearity.This is attributed to the transverse and shear damages in the 0 • and ±45 • layers.Simulations also indicated that the 0 • layers entered the macro-damage stage at a σ y of 256 MPa, and the stress-strain response exhibited a slight deviation from the elastic response (the black dashed line in Figure 7).At 352 MPa, the ±45 • layers entered the macro-damage stage, and the deviation was more pronounced.The stiffness attenuation of the 0 • layers was notable around 540 MPa (with d 2 , d 6 , η E , and η G of the 0 • layers reaching 0.7, 0.8, 0.94, and 0.98 at 544 MPa, respectively).The 90 • and ±45 • layers would carry a greater load, resulting in a gradual change in the curve slope around this stress level.The features of this test, as documented in [38], are as follows: (1) the measured stress-strain curves exhibited relatively linear behavior up to final failure, and (2) final failure took place due to fiber tension fracture.In general, the simulations demonstrated good agreement with the test results and their features.The projected strength was determined to be 736 MPa, exhibiting only a marginal deviation of approximately 2.5% from the test strength of 718 MPa.The final failure mode was characterized by fiber tension fracture, and the stress-strain response in the 0-500 MPa range closely matched.The primary discrepancy observed was that the strain beyond 500 MPa was slightly higher than that of the test results.The rationale behind this discrepancy lies in the slower damage evolution of the constructed constitutive model compared to the actual process, resulting in the significant stress redistribution of the simulation occurring slightly later.This discrepancy can be attributed to the inaccuracy of parameters related to damage evolution, which were selected based on values suggested for the general carbon/epoxy material.A more accurate determination of these parameters specific to AS4/3501-6 could enhance the simulation accuracy.
The comparison of biaxial tension under a stress ratio of  : = 1:2 for these laminates is illustrated in Figure 8.The simulations revealed that the 0° and ±45° layers entered the macro-damage stage at  values of 245 MPa and 352 MPa, respectively.Due to the increased tension loads along their fiber direction compared to the previous case, the nonlinearity was not pronounced.However, a noticeable stiffness attenuation was observed in the ±45° layers around 660 MPa (where  ,  ,  , and  of the +45° layers reached 0.55, 0.78, 0.84, and 0.94, respectively).The 90° and 0° layers would bear more load, resulting in a gradual change in the curve slope around this stress level.The features of this test, as documented in [38], are as follows: (1) the measured stressstrain curves exhibited relatively linear behavior up to final failure, and (2) final failure took place due to fiber tension fracture.In general, the simulations demonstrated good agreement with the test results and their features.The projected strength was determined to be 736 MPa, exhibiting only a marginal deviation of approximately 2.5% from the test strength of 718 MPa.The final failure mode was characterized by fiber tension fracture, and the stress-strain response in the 0-500 MPa range closely matched.The primary discrepancy observed was that the strain beyond 500 MPa was slightly higher than that of the test results.The rationale behind this discrepancy lies in the slower damage evolution of the constructed constitutive model compared to the actual process, resulting in the significant stress redistribution of the simulation occurring slightly later.This discrepancy can be attributed to the inaccuracy of parameters related to damage evolution, which were selected based on values suggested for the general carbon/epoxy material.A more accurate determination of these parameters specific to AS4/3501-6 could enhance the simulation accuracy.
The comparison of biaxial tension under a stress ratio of σ x :σ y = 1:2 for these laminates is illustrated in Figure 8.The simulations revealed that the 0 • and ±45 • layers entered the macro-damage stage at σ y values of 245 MPa and 352 MPa, respectively.Due to the increased tension loads along their fiber direction compared to the previous case, the nonlinearity was not pronounced.However, a noticeable stiffness attenuation was observed in the ±45 • layers around 660 MPa (where d 2 , d 6 , η E , and η G of the +45 • layers reached 0.55, 0.78, 0.84, and 0.94, respectively).The 90 • and 0 • layers would bear more load, resulting in a gradual change in the curve slope around this stress level.
The features of this test are as follows [38]: (1) both the longitudinal and transverse strain are positive (tensile); (2) the curves exhibited a slight change in the slope at around 400 MPa of transverse stress; and (3) the transverse curves exhibited a change in the stiffness near the final points.Overall, the simulations aligned well with the test results and their features.The projected strength was determined to be 830 MPa, deviating by only approximately 2% from the test strength of 847 MPa, and the final fracture strain was almost the same.The most notable discrepancy pertained to the strength at which the slope changed.This discrepancy could be attributed to inaccuracies in the parameters related to damage evolution, similar to the previous case.This resulted in the damage evolution of the constructed constitutive model being slower than the actual process.The features of this test are as follows [38]: (1) both the longitudinal and transverse strain are positive (tensile); (2) the curves exhibited a slight change in the slope at around 400 MPa of transverse stress; and (3) the transverse curves exhibited a change in the stiffness near the final points.Overall, the simulations aligned well with the test results and their features.The projected strength was determined to be 830 MPa, deviating by only approximately 2% from the test strength of 847 MPa, and the final fracture strain was almost the same.The most notable discrepancy pertained to the strength at which the slope changed.This discrepancy could be attributed to inaccuracies in the parameters related to damage evolution, similar to the previous case.This resulted in the damage evolution of the constructed constitutive model being slower than the actual process.

Open-Hole-Tension Strength
The OHT tensile strength was represented by the boundary-surface tensile load  .The comparison between the simulations and test results are listed in Table 2, alongside the predictions of Chang [50], Tan [52], and Chen [42].The average absolute prediction error was 5.9%, which is lower than Chang's 18.2%, Tan's 18.6%, and Chen's 8.6%.The precision of the constitutive model has been improved.

Open-Hole-Tension Strength
The OHT tensile strength was represented by the boundary-surface tensile load σ m .The comparison between the simulations and test results are listed in Table 2, alongside the predictions of Chang [50], Tan [52], and Chen [42].The average absolute prediction error was 5.9%, which is lower than Chang's 18.2%, Tan's 18.6%, and Chen's 8.6%.The precision of the constitutive model has been improved.indicates that the material is absolutely damaged, while '0' denotes a non-damaged state.Figure 9a illustrates the D e (1,1) attenuation of the first layer, which is primarily subjected to longitudinal tension.Damage first occurs near the hole at 180 MPa, and as the load increases, the damaged zone diffuses towards the edge of the specimen.Figure 9b,c show the D e (4,4) attenuation of the second and seventh layers that are mainly subjected to transverse tension and in-plane shear.Damage initially occurs near the hole at approximately 120 MPa, and with an increasing load, the damaged zone diffuses towards the edge until laminate failure occurs.The attenuation of these two laminae is almost the same.Figure 9d shows the D e (2,2) attenuation of the 10th lamina that is primarily subjected to transverse tension.tudinal tension.Damage first occurs near the hole at 180 MPa, and as the load increases, the damaged zone diffuses towards the edge of the specimen.Figure 9b,c show the De(4,4) attenuation of the second and seventh layers that are mainly subjected to transverse tension and in-plane shear.Damage initially occurs near the hole at approximately 120 MPa, and with an increasing load, the damaged zone diffuses towards the edge until laminate failure occurs.The attenuation of these two laminae is almost the same.Figure 9d shows the De(2,2) attenuation of the 10th lamina that is primarily subjected to transverse tension.

Conclusions
This study presents a novel anisotropic damage-plasticity model for FRPs.All inherent parameters can be readily obtained through basic mechanical tests, making the model practical and user-friendly.Initially, a constitutive modeling method was detailed.The 3D Puck criteria were employed to assess the onset of micro-and macro-damages, initiating the evolution process through a combination of CDM and linear stiffness attenuation methods.The evolution was coupled with a one-parameter plastic model, and five basic components of the plastic model were outlined.Subsequently, the material constitutive equations were implemented in Fortran under the Visual Studio framework and linked to ANSYS 19.0 through UserMat.In a numerical implementation, the Cauchy stress was solved using the RMA, and the consistent tangent stiffness was deduced based on the solved stress and the chain derivation rule.Finally, the developed constitutive model was verified using existing biaxial tension and OHT tests of glass and carbon FRP laminates.The simulations demonstrated a high degree of correspondence with the test results.
Further research will prioritize the precise determination of damage evolution parameters.This study opted for parameters related to damage evolution based on the values suggested by the authors in [21,40].However, as different FRPs possess specific damage evolution values, inaccuracies in these parameters could lead to deviations, as described in

( 2 )
Damage evolution.This component assesses the damage development in the AB and BC segments.(3) Plasticity evolution.This component assesses the evolution of plastic deformation in the entire loading process.Of course, each component should maintain an anisotropic form.Elaborated descriptions of each component are presented below.

Figure 1 .
Figure 1.Three stages of the loading process of an FRP unidirectional lamina.Accordingly, the constitutive model should comprise three components: (1) Damage determination.This component evaluates the onset of micro-damages and their progression into macro-damages, specifically identifying points A and B. (2) Damage evolution.This component assesses the damage development in the AB and BC segments.(3) Plasticity evolution.This component assesses the evolution of plastic deformation in the entire loading process.Of course, each component should maintain an anisotropic form.Elaborated descriptions of each component are presented below.

Figure 1 .
Figure 1.Three stages of the loading process of an FRP unidirectional lamina.Accordingly, the constitutive model should comprise three components: (1) Damage determination.This component evaluates the onset of micro-damages and their progression into macro-damages, specifically identifying points A and B. (2) Damage evolution.This component assesses the damage development in the AB and BC segments.(3) Plasticity evolution.This component assesses the evolution of plastic deformation in the entire loading process.Of course, each component should maintain an anisotropic form.Elaborated descriptions of each component are presented below.

Figure 3 .
Figure 3. Numerical algorithm of the constitutive model.

Figure 3 .
Figure 3. Numerical algorithm of the constitutive model.

Figure 4 .
Figure 4. Numerical models of biaxial load and OHT.

Figure 4 .
Figure 4. Numerical models of biaxial load and OHT.

Figure 9
Figure 9 depicts the stiffness attenuation of the [0/(±45 • ) 3 /90 • 3 ] s laminate labeled B near the open hole.The load direction is in the horizontal direction, and factor '1'indicates that the material is absolutely damaged, while '0' denotes a non-damaged state.Figure9aillustrates the D e (1,1) attenuation of the first layer, which is primarily subjected to longitudinal tension.Damage first occurs near the hole at 180 MPa, and as the load increases, the damaged zone diffuses towards the edge of the specimen.Figure9b,cshow the D e(4,4) attenuation of the second and seventh layers that are mainly subjected to transverse tension and in-plane shear.Damage initially occurs near the hole at approximately 120 MPa, and with an increasing load, the damaged zone diffuses towards the edge until laminate failure occurs.The attenuation of these two laminae is almost the same.Figure9dshows the D e (2,2) attenuation of the 10th lamina that is primarily subjected to transverse tension.

Table 1 .
Parameters of three laminae.
Plastic coefficients b, Y 0 , Y c , Y ′ 0 , Y ′ ′Release rate of damaged strain energy