Hysteresis behaviour modelling of woven composite using a collaborative elastoplastic damage model with fractional derivatives

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 extensive use of composite materials in the industrial applications requires a careful examination of their mechanical properties. The material behaviour must be validated experimentally for different types of loadings (static, dynamic, fatigue) depending on the material application. Then, a behaviour model should be developed to simulate structures' response as a supportive tool to the experimental method. As the composites are anisotropic heterogeneous materials, complex models are required to adequately describe their response. Some works referred by [1][2][3][4][5][6] is concerned with the composites behaviour modelling in the mesoscale. Another option is to use the micromechanics models [7,8]. Based on the in-ply damage scenarios at the microscale, a recent work [8] considers the coupling between the delamination and micro cracking. These models take into account only the damage propagation, the material hardening and the in-elastic strains' appearances. Furthermore, the polymer matrix composites have the viscoelastic behaviour due to the nature of matrix and this behaviour is evident in the material's strain-rate sensibility. The dependence of a material's strain rate sensitivity on the viscoelastic effects was first taken into account in the works [9][10][11] for the unidirectional composites and in the works [12][13][14] for the woven composites.
However, the coupling between the cyclic loading damages and the viscoelastic effects has not been considered in the classical damage modelling [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. It must be highlighted that the damage propagation is strongly linked with the material energy dissipation and, thus, with the hysteresis loops appearance. Thus, the precise modelling of hysteresis loops is needed. Therefore, in this paper we propose a collaborative model based on a fractional derivative approach.
A fractional derivative theory is a promising technique to describe the history-dependent phenomena like viscoelasticity. A significant number of fractional viscoelastic models have been developed and successfully applied to model a mechanical response of natural structures and modern heterogenic materials such as elastomers and polymers [15][16][17]. The classical rheological models are rewritten in terms of fractional derivative like the Kelvin-Voigt [18,19], Maxwell [20][21][22] and Zener [23,24] fractional models that allows to obtain creep and relaxation functions adapted for the polymers' response. Rabotnov has developed a hereditary theory where he had introduced a generalized fractional rheological model and a new class of hereditary functions in his works [25][26][27]. The Rabotnov's theory is widely used to describe the behaviour of the polymers, metals and concrete. Bagley and Torvik established the fractional law in the frequency domain to describe the behaviour of polymers and elastomers [28]. The link between molecular structure of polymers and fractional derivative law is given in the works [29][30][31] for the different types of polymers. These relations were obtained by using the Rouse molecular theory [32]. The hysteresis loops' appearance is another viscoelastic effect which can be observed in the viscoelastic materials response and then, it can be treated using the fractional derivative approach. Caputo has described the hysteresis loop using the fractional derivatives in the frequency domain to determine a fatigue limit for a few metals [33]. Mateos [34] proposed a fractional derivative model with a few material parameters to describe the hysteresis composite behaviour in the temporal domain. This model also includes strainrate dependence and can be applied for dynamic loading [35]. The elastic and the irreversible strains, the damage and the strain rate effects are taken into account into [34] and [35] but not the material hardening. Moreover, these phenomena that are cited are inseparable making it difficult to calculate their individual contribution in the material's energy dissipation. Hence this paper proposes a collaborative model to describe the complete composite behaviour during cyclic loading including the representation of the hysteresis phenomena during the loading-unloading paths. The proposed model is composed of the elastoplastic damage behaviour law [3] with possible strain-rate sensitivity [9] and including a newly proposed fractional derivative approach. This model allows us to distinguish the dissipation due to the material damage, plasticity and viscoelasticity.
This paper is organized as follows. In Section 2, the essential mathematical elements of the fractional calculus are introduced and the collaborative model for woven composite materials is proposed in the Section 3. The parameter identification procedure is developed in the Section 4 and the validation of the collaborative model is demonstrated in the Section 5. The Section 6 concludes this paper and notes the various engineering application of this model.

Introduction in fractional calculus
The advent of advanced materials requires new mathematical methods to describe their behaviours. One of the ways is to increase the number of parameters in the existing constitutive laws. This leads to sophisticated models with a significant number of parameters which sometimes require complex identification procedures. On the other hand, new methods may be introduced as the fractional derivative theory.
Fractional operators are an extension of the order of the integrals and derivatives to the fractional order. By their definition, the fractional operators take into account the past loading history and are suitable to describe the viscoelastic phenomena. The solution of fractional differential equations is expressed by the fractional functions [36,37] whereas the solution of ordinary differential equations is an exponential basis. Therefore, the fractional behaviour laws are widely used to describe behaviour of heterogeneous materials and fractal structures such as polymers, alloys and geological strata [27,38,39].
There are a few definitions of the fractional integrals and derivatives which are equivalent for most of the functions but they do not always lead to the same results. In this work the classical Riemann and Liouville definition (6) [40,41] is used. To deduce this definition, the formula of nth repeated integral of the continuous function f ðxÞ (Cauchy's formula) is used: To generalize the Eq. (1), the natural number n is replaced by a real positive number a. It is known that the extension of the factorial function to real and complex numbers is the gamma function CðzÞ defined by: The gamma function CðzÞ satisfies the following recurrence relation: The fractional Riemann-Liouville integral of order a (4) is obtained from the generalized form of the Eq. (1) by using the property of the gamma function (3) ðI a a f Þ RL ðxÞ¼ def 1 CðaÞ where ðI a a f Þ RL ðxÞ is the fractional integral of order a of the function f ðxÞ.
The fractional derivative is an inverse operator to the fractional integral. Let us consider a which is a positive real number and n which is an integer that satisfies the inequality n À 1 < a < n. The fractional derivative of order a is defined as the fractional integral of order n À a derived n times. For example, if a belongs to the open interval ð0; 1Þ, the fractional derivative of order a of the function f ðxÞ is the fractional integral of order 1 À a derived one time: The alternative form of the previous expression is: where ðD a a f Þ RL ðxÞ is the fractional derivative of order a in the sense of Riemann-Liouville.

Numerical evaluation of fractional derivatives
In order to implement the fractional derivatives in a numerical code, different algorithms exist. The L1-approximation [42] of the Riemann-Liouville fractional derivative (6) ðD a 0 f ðxÞÞ L1 has a following form: where Dx is a step interval, N is a number of interval points and C is the gamma function defined by (2). The Riemann-Liouville fractional integral (4) can be rewritten in the alternative forms [43]: The formulas (8) and (9) can be easily implemented in the numerical code if an analytical expression of the function f ðxÞ is known.
According to the expression (6), the fractional derivative is the first derivative of the M1 (8) or M2 (9) integrals. The M1 and M2 fractional derivate approximations are calculated by a central difference scheme: ðD a 0 f ðxÞÞ M2 ¼ In order to estimate the numerical accuracy of the proposed approximations, the fractional derivatives of the linear function f ðxÞ ¼ x are computed by the methods (7), (10) and (11) and compared with the analytical results (12) [36].
The order of fractional derivatives is set a ¼ 0:3. The computational interval x 2 ð0; 1Þ is discretized by N ¼ 20 points. The results are illustrated in the Fig. 1. Accuracy of the proposed methods is estimated by using the relative error d relative i (13) and the total error d (14) functions.
where D a f ðx i Þ analyt is the fractional derivative computed analytically by the Eq. (12), D a f ðx i Þ approx corresponds to the L1, M1 or M2 approximations defined by the Eqs. (7), (10), (11) respectively. The integrals (8) and (9) are computed using the Gaussian quadrature. The computed interval ð0; 1Þ is divided in ten subintervals and five gauss-points are used in each sub-interval.
The results are presented in the Fig. 2 and in the Table 1. The biggest error is observed near the zero point. The total relative errors of M1 and M2 algorithms are less than 1%. The M1-method is the most accurate and is used in further calculations.

Theoretical model for composite ply
The material which is used in this work is a carbon woven fibre epoxy matrix composite. The traction cyclic tests have been carried out on an Instron 3369 machine with a strain rate of 5 mm/min in two different orthotropic directions: fibre and shear directions. The experimental result of the tensile cyclic test for the woven carbon/epoxy composite is illustrated in the Figs. 5 and 6. Subscripts 1 and 2 represent the warp and the weft directions, respectively. The composite is considered to be perfectly balanced and thus the longitudinal and transverse behaviours are considered to be equivalent. Typically, the experimental tests show an elastic brittle response in longitudinal (transverse) direction (Fig. 5). However, shear behaviour has a non-linear character. The irreversible strains, the in-ply damage (which is observed from a regular decrease of shear modulus) and the hysteresis loops are presented in the stress-strain curve (Fig. 6).
The collaborative model is developed to represent the elastoplastic damaged material behaviour as well as the hysteresis loops for woven composites under cyclic loading. The model is based on the collaboration of two sub-models. The first one deals with a composite behaviour during loading path. The elastic and the in-elastic strains are computed as well as the in-ply damage. The      second sub-model involves the fractional derivative approach to represent the hysteresis loops during unloading path. The mesoscale model is developed under a state of plane stress. The constitutive equations are deduced within the framework of irreversible thermodynamics processes using the local state method.

Thermodynamic aspects
Within the framework of irreversible thermodynamics, a Helmholtz potential qw is chosen as a function of internal variables: where e e , d i , p are internal variables associated with elastic strain, damage in the orthotropic directions and cumulated plasticity respectively.

Damage modelling
Under the given loading, the damage appears in the elementary ply. It leads to the degradation of the mechanical properties and then to the material failure. The continuum damage mechanics theory is used to describe the degradation of elementary ply [44]. Using the effective stresses notation [45], the degraded stiffness matrix C is introduced as follows: When the material is not degraded: The constitutive equations are deduced using the second principal of thermodynamics. The elastic strain energy of the damaged material W d e (equivalent to the Helmholtz potential qw) has a following form: where r is the stress vector.
The thermodynamic forces associated with internal variables d ij are defined as following: where Y 11 ; Y 22 and Y 12 are thermodynamic forces associated with damage internal variables d 11 ; d 22 and d 12 .
The associated thermodynamic forces characterize the damage propagation. The state of damage can only grow [1,3,46] and therefore, the threshold of undamaged zone Y ij is defined as a maximal thermodynamic force for all previous time ðsÞ up to the current time ðtÞ [12]: The coupling between the internal damage variables and the associated thermodynamic forces is determined by the experimental data fitting and the functions given below are frequently used for this: Linear law Heaviside function H (This function is used to describe elastic brittle behaviour) where Y 0 ij is initial damage thresholds and Y R ij is failure-damage thresholds in the orthotropic directions. The constants in these laws are the material parameters.
The Eqs. (21) and (22) are not always able to adequately represent the damage propagation especially in the composite material with thermoplastic matrix. Thus, the polynomial law is proposed [47] to increase the performance of the elastoplastic damage model:

Plasticity modelling and damage-plasticity coupling
The experimental data shows the irreversible strains appearing mainly in shear [6,9,12,13]. Thus, plastic flow is considered to be blocked in fibre directions:  12 are the components of the plastic strain. The damage and plasticity coupling is made using the effective stress notationr 12 (25).
The isotropic strain hardening is assumed. The elastic domain is defined by the yield function f : where R 0 is a yield stress and the function RðpÞ is a material characteristic function of the cumulative plastic strain p. Generally, the hardening function RðpÞ is approximated by a power law: where b and k are material parameters identified from the experimental data. The material parameter identification is detailed in the next sections. The resulting stress-strain shear curve is obtained by the elastoplastic damage model [3] for the carbon/epoxy woven composite (Fig. 7). The simulation results are in a good agreement with the experimental data. However, the hysteresis loops cannot be reproduced by this model.

Fractional derivative model to represent hysteresis loops
The woven composites have a hysteresis behaviour for certain loading types such as a cyclic shear test. The reason of this is that the matrix is subjected a bigger part of the shear load and as a result, the viscoelastic response of the polymer matrix can be observed as well as the gradual development of microscopic damage in the matrix or along matrix-fibre interface. This leads to the appearance of hysteresis loops which is associated with material energy dissipation. Since we are dealing with a viscoelastic phenomenon, the loading history has to be taken into account. That is why the fractional derivatives are introduced in the behaviour law to provide the material history dependence.
To deduce a behaviour law during an unloading path let us consider one hysteresis loop. According to the elastoplastic damage model, the plastic strain stays constant and the elastic strain is a linear function of time. However, a non-linearity in strain is observed in the experimental curve during the unloading/reloading path (Fig. 8).
Commonly, an additional viscoelastic strain vector e ve is introduced to take into account the observed non-linearity. In this case, the total strain e t 12 is a sum of three strains: elastic e e 12 , plastic e p 12 and viscoelastic e ve 12 : Generally, this approach leads to the complex models with a large number of parameters.
An alternative approach is proposed in this work. In order to model the viscoelastic effect, the fractional derivatives are introduced in the governing law. The total strain within the hysteresis loop e loop 12 is defined as following:

ðtÞ ð 29Þ
where e e 12 is the elastic strain determined by the elastoplastic damage model, D a is the Riemann-Liouville fractional derivative (6) and A, B and a are the fractional model parameters.
As the plastic flow stays constant, the stress is expressed by the elastic law: The fractional derivative of the elastic strain D a e e 12 ðtÞ provides a strain non-linearity. Hence, the hysteresis loop appears (Fig. 9).
The fractional order a adjusts the loop size. The growth of a leads to the increasing of loop area and consequently the dissipation of viscoelastic energy increases.
The factor G 1 12 governs the slope of the loop (Fig. 10). The damage level is adjusted in respect to the elastoplastic damage model. The term G 0 12 ð1 À d 12 ÞA is a connecting stress that links two submodels: the fractional and the elastoplastic damage models (Fig. 11). Thus the position of the loop is adjusted and the stress becomes a continuous function.

Coupling of two sub-models
In this part, the coupling between two sub-models is considered.  Let us consider the typical response in shear of the carbon/ epoxy woven composite (Fig. 12). The collaboration between two sub-models is ensured by the yield function f (26). There are two possible cases:

Methodology of fractional derivatives implementation
Let us consider the ''unloading-reloading" intervals on the stress-strain curve (Fig. 12) which are indicated by the red circles. One of the intervals is represented for the ''time-strain" axes in the     In other words, to reproduce the hysteresis loop on the ''unloading-reloading" time interval ½t M ; t N , the fractional derivatives are computed on the ''loading-unloading-re loading" interval ½t L ; t N to avoid the mathematical problems.
The fractional operators are approximated by using the M1method expressed by the formulas (8) To use the formulas (32) and (33), the analytical expression of the e e 12 ðtÞ function has to be reconstructed by using the piecewise function.

Parameters identification
In the following sections the identification procedure for the model parameters is described using the experimental data for the carbon/epoxy woven composite (Figs. 5 and 6). The in-ply damage, the isotropic strain hardening and the hysteresis loops modelling are taken into consideration.

Materials characterization in fibres directions
The ½0  Table 2.

Damage identification in fibre directions
Due to the brittle elastic behaviour, the material failure appears instantly. The Heaviside step function (22) is used to describe the damage propagation. The initial damage threshold coincides with the failure threshold: Y 0 11ð22Þ ¼ Y R 11ð22Þ . Damage parameters are provided in the Table 3. The damage propagation is illustrated in the Fig. 14.

Materials characterization in shear
The material parameters are identified from the shear experimental curve for thermoset composite (Fig. 6). The shear undam-aged modulus G 0 12 and the yield stress R 0 are identified from the linear regression of the elastic part in the experimental stressstrain curve. The failure stress r R 12 corresponds to the maximum reached stress during the test. The elastic parameters in shear are summarized in the Table 4.

Damage and plastic strain evaluation in shear
Typical cyclic shear tests show a reduction in the material stiffness. To model the material degradation, the damage variable d 12 is linked with the associated thermodynamic force ffiffiffiffiffiffiffi ffi Y 12 p . The damage level is measured by the shear modulus reduction: where G i 12 is a degraded shear modulus. Using the experimental data we state that damage grows linearly (Fig. 15) in the carbon/epoxy composite in accordance with the Eq. (21). In this case, the coefficient is the speed of damage propagation. The damage parameters are presented in the following Table 5. The cumulative plastic strain p is obtained using the identified damage parameters. The strain hardening function RðpÞ is defined by the Eq. (35).
The experimental data are fitted using the power law (26). The approximation is illustrated in the Fig. 16. The identified parameters are presented in the Table 6.

Identification of fractional derivative model parameters
The fractional model (29) parameters A; B and a are identified by solving a non-linear optimization problem with constrains. The main constrain is imposed on the fractional derivative order a: 0 < a < 1. The objective function is defined as a relative error: Table 3 Longitudinal (or transverse) damage parameters for the carbon/epoxy composite.

Parameter
Carbon/epoxy ffiffiffiffiffiffiffi ffi Y 0   is a shear strain calculated by the fractional model (29) and N is a number of time-points inside the considered interval. For the carbon/epoxy composite, the fractional parameters are represented by linear functions of damage as following: The results are illustrated in the Figs. 17-19. The coefficients of fittings (38), (39) and (40) are collected in the Table 7.

Results
The carbon/epoxy composite behaviour is represented numerically in fibre directions by the collaborative model using the identified parameters in Fig. 20 and the results of simulation are more than satisfying since the exact replica of the elastic brittle behaviour is obtained.
The collaborative model is applied to represent the shear strainstress curve. The results are in good agreement with the experimental data (Fig. 21). Small errors are observed at the first and the last points of the hysteresis loop. The error at the first point is related to the mathematical particularity of fractional calculus [43]. The error at the last point is related to the uncertainty of the optimal solution. Furthermore, the derivation of the piecewise function provides an inaccuracy at the inflection point (the transition from the unloading to the loading path).
The identification of parameters is also made for the thermoplastic woven composite: carbon/PA66 in shear direction. The damage propagation has a non-linear character and is expressed by the third order polynomial law (23). The power law (27) has been chosen to fit the hardening function. The symmetrical hysteresis loops are transformed into the ''bean" shaped loops due to the fibre reorientation (Fig. 22) and this increases the error in damage identification. To improve the simulation results, the reorientation of fibre angle has to be taken into account. Due to the nonlinear damage law, the fractional parameters A and B also has a non-linear evolution and are approximated by the third order polynomial law. The fractional derivative order a has a linear dependence on damage but its values rises because of the loop size increments.
The collaborative model is applied to simulate thermoplastic composite response. The difference between the experimental and simulation results is observed towards the end of the resulting stress-strain curve (Fig. 22) when the level of damage is high. Thus the inexact damage identification provides a supplementary error in addition to the previously cited fractional calculus problems. The asymmetric hysteresis loops cannot be treated by the collaborative model. Despite the unsymmetrical loops form, their area measures correctly and the amount of dissipated energy can be determined precisely.
In conclusion we state that the collaborative model provides more than satisfactory results for the composite materials with different matrixes.
An additional simulation has been made to predict the appearance of hysteresis loops during cyclic loading. A six cycles' experimental shear test is replaced by the fictitious eight cycles loading for the thermoplastic composite (Fig. 23). The simulation is made with the previously identified parameters. The results are shown    The composites with polymer matrix are sensitive to the strain rate. The ability of the collaborative model to represent hysteresis loops in dynamic loading is shown with the example of thermoset composite for two different strain rates: 5 mm/min and 50 mm/ min. The material response for strain rate of 50 mm/min is predicted by the collaborative model from the quasi-static (5 mm/ min) test identification (Fig. 25). The collaborative model is a promising approach to predict the hysteresis loops' appearance and to characterize the damage propagation at different strain rates.

Conclusion
The collaborative model is developed to simulate a complete dynamic behaviour of woven composites. This model consists of two sub-models: the first one is the elastoplastic damage model [3] and the second is the fractional derivative model. The elastoplastic damage model is able to represent the composite material response except the hysteresis loops. To fill this gap, the fractional model is introduced in order to describe the history-depended phenomena. The fractional constitutive law is expressed by a non-standard Kelvin-Voigt model using Riemann-Liouville fractional derivatives. Just a few parameters are required to represent the hysteresis loop in this model. The fractional model is applied during the unloading path which significantly reduces the computation time. Moreover, the fractional model is an optional tool to the elastoplastic damage model and can be used only if the hysteresis modelling is required.
The collaborative model is validated for thermoset and thermoplastic carbon woven composites. Both materials have an elastic brittle behaviour in fibres directions. The simplified identification procedure is used. The material response in shear is characterized by the damage propagation, the irreversible strains appearance and the hysteresis mechanisms. The shear behaviour was simulated by the collaborative model. The results of simulation are in good agreement with the experimental data. Even if the collaborative model is not able to represent unsymmetrical hysteresis loops for thermoplastic composites, their area is computed correctly and the dissipated energy can be estimated precisely. The fibre reorientation during the test has to be taken into account to increase the model performance.
The collaborative model is also able to predict the appearance of additional hysteresis loops. It is an important result which can lead to the simplification of experimental characterization campaign. The model is able to represent the hysteresis loops at different strain rates from the quasi-static experimental data. It is a promising approach to signify damage propagation in crash analysis.
The collaborative model is a simple and efficient approach to make a link between the mechanical and the thermodynamic simulations. The amount of dissipated energy during the cyclic loading can be easily calculated. Therefore, the collaborative model can be used as a numerical support for the self-heating experimental method [48]. A method based on the self-heating tests provides fast fatigue limit identification for the composite materials compared to conventional methods (S-N curves). Using the proposed behaviour model, the thermo-mechanical finite element simulation can be developed in order to compute accurately the dissipated energy.