L Augustins
email: louis.augustins@mpsa.com

R Billardon
email: rene.billardon@safranmbd.com

F Hild
email: hild@lmt.ens-cachan.fr

F Hild Lmt

Constitutive model for flake graphite cast iron automotive brake disks From macroscopic multiscale models to a 1D rheological description

Keywords: Flake Graphite Cast Iron, Multiscale model, Constitutive model

published or not. The documents may come

Constitutive model for flake graphite cast iron automotive brake disks: From macroscopic multiscale models to a 1D rheological description L Augustins, René R Billardon, François

Introduction

The automotive industry is increasingly facing major challenges related to higher standards decreasing the pollution caused by the emission of greenhouse gases, and the resulting need for lighter structures. Optimal design implies having efficient tools to accurately predict the life of structures under any load and especially fatigue loadings. One possible strategy is based on the fundamental assumption that it is possible to decouple the cyclic mechanical response of the material from fatigue damage [START_REF] Constantinescu | A computational approach to thermomechanical fatigue[END_REF]. It is usually referred to as a locally coupled approach [START_REF] Benallal | Continuum damage mechanics and local approach to fracture: Numerical procedures[END_REF]. The latter requires a number of assumptions to be checked, which are mainly related to the stabilization of the local mechanical variable histories under cyclic loadings.

This work focuses on the thermomechanical fatigue behavior of automotive brake disks. One of the critical points of the design process, as part of the previously described approach, is the correct description of the anisothermal thermomechanical behavior of the material, namely, flake graphite cast iron (FGCI) in most applications. The thermomechanical properties of FGCI depend on the quality and the distribution of flakes and on the matrix composition. Although graphite flakes, with a length ranging from 0.03 to 1 mm, have very little influence on the compressive strength of the matrix, they considerably reduce its tensile strength. Owing to the very low strength of graphite, the flakes behave like cracks [START_REF] Clough | The deformation and rupture of grey cast iron[END_REF]. This phenomenon induces a notch sensitivity, causing FGCI to become more brittle than spheroidal graphite cast irons [START_REF] Meyersberg | Sur la signification du module d'élasticité de la fonte[END_REF]. When macrocracks initiate, they propagate almost instantaneously. In FGCIs, flakes can be seen as microcracks leading to tensile failure (without large elongations or necking) following the graphite flake distribution. Flake debonding under tensile loadings is one of the causes of the nonlinear stress-strain response for relatively low stress values [START_REF] Gilbert | An evaluation of the stress-strain properties of flake graphite cast iron in tension and compression[END_REF], which induces a strong asymmetry of the tension/compression behavior.

Several multiscale models of FGCIs, which are based on finite element simulations, have been devised [START_REF] Russell | Finite element simulation of the microstructure of gray cast iron[END_REF][START_REF] Downing | Stress/strain simulation model for grey cast iron[END_REF][START_REF] Haenny | The role of the matrix graphite interaction in the tensile behavior of grey cast iron[END_REF] to study the influence of the microstructure, in particular the size of flakes, on the mechanical response. The authors were able to simulate cyclic tests with good agreements. The major drawback of these models, which postulate an initially damaged condition, is that they do not allow the simulation of a progressive debonding of graphite flakes; possibly, the modeling of this phenomenon is very complex and cannot be easily parameterized.

Further, many models, either specific or not to FGCIs [START_REF] Mohr | Welche Umstaende bedingen die Elastizitaetsgrenze und den Bruch eines Materials[END_REF][START_REF] Drucker | Soil mechanics and plastic analysis or limit design[END_REF][START_REF] Raghava | The macroscopic yield behavior of polymers[END_REF][START_REF] Hjelm | Yield surface for grey cast iron under biaxial stress[END_REF][START_REF] Josefson | A non linear kinematic hardening model for elastoplastic deformations in grey cast iron[END_REF][START_REF] Altenbach | On elastoplastic deformation of grey cast iron[END_REF][START_REF] Brooks | The development of a material model for cast iron that can be used for brake system analysis[END_REF], account for a dissymmetry between the behavior in tension and in compression through the modification of the yield surface, by introducing stress invariants, such as the hydrostatic pressure, in the formulation. The most commonly used model is that proposed by Hejlm [START_REF] Hjelm | Yield surface for grey cast iron under biaxial stress[END_REF], which is implemented in the commercial finite element code Abaqus. These models make it possible to represent the overall dissymmetry but fail to account for either the loss of stiffness related to flake debonding, or the stiffness recovery. As an illustration of this lack of precision, the simulation of a symmetric strain controlled cyclic test was achieved using a two-kinematic hardening elastoplastic model with a Drucker-Prager yield surface. Results of this simulation are shown in Figure 1, where the limits of such models are exhibited.

In this paper, after a brief presentation of the experimental database, other multiscale models are proposed leading to a precise description of FGCI behavior. An energetic analysis of these models allows their thermodynamic accuracy to be validated and leads to the formulation of a one dimensional rheological model, which is able to describe precisely the material response under cyclic loading histories.

Experimental database

In order to study the complex behavior of brake discs gray cast iron, a comprehensive experimental campaign was conducted. The material used to make automotive brakes is a pearlitic FGCI whose chemical composition is presented in Table 1. The flake distribution is uniform, with sizes ranging between 60 and 250 µm. It is therefore of type IA4-5 [START_REF]AFNOR: Fonderie -fonte à graphite lamellaire[END_REF]. Figure 2 shows a micrograph of the studied grade. The experimental database includes monotonic tensile and compressive tests (Figure 3), cyclic tests with dwell phases in tension and compression (Figure 4). All the tests were performed from 20 • C up to 600 • C on cylindrical samples, which are machined from the ring of a brake disk, in order to be representative of actual components. For the monotonic tests shown in Figure 3, Young's moduli in tension and compression are the same, but a nonlinear stress-strain relationship appears in tension for low stress values (i.e. less than 100 MPa). Furthermore, stress levels achieved in compression are much higher than those achieved in tension, for symmetric strain levels. Cyclic tests consist of a succession of cycles performed at a given strain level ε 1 . This level varies during the test and its initial value is ε 1 = 0.1%. A cycle is composed of: -a compressive phase from ε = 0 to ε = -ε 1 of 2.5-s duration; -a dwell phase in compression when ε = -ε 1 of a temperature related duration; -a compressive/tensile phase when ε = -ε 1 to ε = ε 1 of 5-s duration; -a dwell phase in tension when ε = ε 1 of a temperature related duration; -an unloading phase back to zero strain of 2.5-s duration.

After six cycles at this strain level ε 1 (leading to plastic shakedown), the next strain level is considered. During the test, the following values are successively reached for ε 1 : 0.1%, 0.2%, 0.3%, 0.4% and 0.5%. A strong dissymmetry is observed between tensile and compressive phases (Figure 4), not only in terms of maximum stresses but also in terms of observed elastic modulus. Another particularity is the form of the hysteresis loop, which is characterized by the presence of an inflection point.

Some observations, which are summarized in Figure 5, are derived from the study of the experimental results: -a strong dissymmetry exists in terms of maximum stress level between tension and compression; -during loading phases, the current Young's modulus is identical to the initial modulus; -Young's moduli during the unloading phases are lower than the initial modulus; -during unloading phases, there is an inflection point corresponding to the gradual closure of cracks induced by matrix/graphite debonding during tensile phases, which gives rise to an increased macroscopic stiffness. The detailed understanding of strain and damage mechanisms in FGCI is the key step in the development of a constitutive model. As mentioned in the introduction, multiscale approaches offer a number of answers but fail to simulate accurately the behavior of the material. In this section, unit cells representative of the material are introduced to quantify accurately the phenomena of debonding, opening and closure of the graphite flakes. All the computations are conducted with a general finite element solver (i.e. Abaqus) that uses an implicit integration scheme.

Constitutive model of the matrix

With the experimental database, a constitutive model for the matrix is defined. Cyclic tests show that there is no cyclic softening nor hardening. This result suggests that hardening is purely kinematic (i.e. no isotropic component). Monotonic compressive tests exhibit a highly nonlinear portion in the stress-strain relationship, which is followed by a linear response up to 2% strain (Figure 6). The choice is made to use two kinematic hardening variables [START_REF] Chaboche | Constitutive equations for cyclic plasticity and cyclic viscoplasticity[END_REF], one nonlinear X 1 [START_REF] Armstrong | A mathematical representation of the muultiaxial Bauschinger effect[END_REF] and another one X 2 , which is linear [START_REF] Prager | A new method of analyzing stress and strains in work hardening plastic solids[END_REF]. Since this section only deals with simulations of the behavior at room temperature, an elastoplastic law (i.e. no viscoplasticity) is considered. The representation of the behavior at higher temperatures is addressed in Section 4. Table 2 presents the parameters identified at 20 • C. The simulation shown in Figure 6 describes very well the experimental data and validates the choice of a two-kinematic hardening model. The first model is inspired by that proposed by Russel [START_REF] Russell | Finite element simulation of the microstructure of gray cast iron[END_REF] but displays significant differences, notably in terms of the description of flake debonding. An axisymmetric cell having a flat central crack is defined thereby allowing the three-dimensional effects of stress concentrations at the crack tip to be accurately represented. This crack is an idealization of the graphite flakes perpendicular to the loading direction. The crack surfaces are in contact accounting for closure during unloading phases. The crack propagation is controlled by the applied stress or strain through the use of the *Debond option in Abaqus, which allows the release of initially linked nodes. Thus, the progressive debonding of the graphite flakes, depending upon the applied load, is only parameterized by the crack size.

The boundary conditions applied to the cell are (Figure 7):

the upper horizontal face remains horizontal during loading: y HF = y p1 ; -the lateral vertical face remains vertical during loading x LF = x p2 . The cell mesh is fine enough to simulate accurately crack opening, namely, elements near the crack tip are 50 µm in size for a maximum cell radius of 7 mm. An attempt is then made to identify the crack growth law as a function of the applied stress, in order to obtain a simulated stress-strain curve as close as possible to the stress-strain curve recorded during the monotonic tensile test. The crack growth law is determined "manually" in this exploratory phase, namely, a 20 point tabular evolution of the crack length a as a function of the applied stress is determined step by step from the monotonic tension test. The principle is to fit the experimental tensile test with the appropriate crack size evolution. It is then easy to simulate the cyclic test, incrementing the crack size to its corresponding value, using the previously defined crack opening law, when switching to the new strain level amplitude. The simulation of the cyclic test begins with a strain range of 0.2% and a pre-existing crack. The response of the model is shown in Figure 8. It is found that the correlation between the model and the experimental data is very good not only during loading phases but also during unloading (i.e. elastic) phases, with a reduced stiffness very close to the measured one. Because the constitutive law is only elastoplastic, the dwell phases with stress relaxation are not simulated. The inflection point that corresponds to crack closure is represented accurately. However, the model predicts that this point, which is experimentally stationary, shifts slightly in the area of negative strains for each increment of the strain amplitude. One hypothesis is to consider that this point is experimentally stationary owing to a wide variety both in size and distribution of the flakes, allowing a gradual closure. Another hypothesis is to consider that the single crack model induces significant stress concentrations at the crack tip, most likely higher than those present in the actual material. The tensile stress concentration generates compressive residual stresses and a large area of positive plastic strains, which then "delay" the complete closure of the crack.

Damage identification The numerical model allows the concept of damage to be addressed. The crack is related to a damage variable D through a decrease of the effective section. The effective stress σ introduced by Kachanov [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF] is defined as

σ = σ (1 -D) (1)
The damage variable then reads

D = 1 - Ẽ E ( 2
)
where Ẽ is the current Young's modulus, calculated from an elastic simulation.

Knowing the crack size a as a function of the applied stress σ, it is possible to directly link D to σ. The change of D with the applied stress is thus obtained as shown in Figure 9. This result indicates that the flake debonding only occurs after a threshold value and then follows a logarithmic evolution. The yield stress in compression is about 180 MPa, which means that damage in tension occurs during the elastic regime. Along the same lines as the elastoplastic bar model of Downing [START_REF] Downing | Stress/strain simulation model for grey cast iron[END_REF], and given the knowledge of the damage growth law as a function of the applied load, a rheological model made of five bars (i.e. two-dimensional plane stress elements) with different sections is considered to represent the cyclic test shown in Section 2, to the exclusion of any other test. Figure 10 schematically shows the way the fivebar system is working. Bar number 1 with the largest section accounts for the uncracked ligament whose vertical displacements are forbidden on its lower side. Its section S 1 is directly linked to the maximum damage level corresponding to the strain level of 0.5 %, denoted by D 0.5 % . The section of bar 1 is defined as

S 1 = S T (1 -D 0.5% )
where S T is the total section of the k bars, S T = k S k . The other four bars correspond to the cracked area. They are in contact with a rigid support, and each bar corresponds to a discretization of damage, which is linked to a strain level. The contact can be open or fully bonded. For a strain level of 0.5 %, the contact is open for all bars except bar number 1. For a strain level of 0.4 %, the contact is open for bar number 5, 4 and 3 and closed elsewhere. Section S 2 of bar 2 reads

S 2 = S T (1 -D 0.4 % ) -S 1
For a strain level of 0.3%, the contact is open for bars number 5 and 4 and closed elsewhere. Section S 3 of bar 3 becomes

S 3 = S T (1 -D 0.3 % ) -S 1 -S 2
For a strain level of 0.2%, the contact is open for bar number 5 and closed elsewhere so that section S 4 reads

S 4 = S T (1 -D 0.2 % ) -S 1 -S 2 -S 3
Last, the section of bar number 5 reads

S 5 = S T -S 1 -S 2 -S 3 -S 4
Displacements of the upper faces of all the bars are kinematically connected to a single master node that can be driven in a load or displacement controlled mode. During the simulation of the cyclic test at 20 • C, the contact condition of the cracked bars 2, 3, 4 and 5 changes from bonded to open, depending on the simulated strain level. The displacement of the master node, as in the test, is "strain" controlled.

The simulations of the rheological model are shown in Figure 11. The predictions are very close to both the experimental data and the results obtained with the axysimmetric crack model. There is still a slight shift of the inflection point toward negative strains at each strain level increment, but the amplitude of the shift is much smaller. This shift can be explained by the absence of effects related to the high stress concentrations in the vicinity of the crack tips, a pattern that is more representative of the actual microstructure characterized by a high density of flakes isotropically distributed.

Energy balance analysis

In order to compare the two previous models from an energetic point of view, an energy balance analysis is performed on both models.

Thermodynamic framework

The thermodynamic framework defined below is that used for the undamaged matrix representation. For an elastoplastic model with one nonlinear and one linear kinematic hardening, state variables are ε, εp and α associated with the thermodynamic forces σ, Λε p and X.

Von Mises' yield surface is written as

f (σ, X) = (σ -X) eq -σ 0 (3)
where X = X 1 + X 2 , and σ 0 the initial yield stress. The free energy is chosen as the state potential [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF]

ρω = 1 2 (ε -εp) : C : (ε -εp) + ωp(α 1 , α 2 ) ( 4
)
where ρ is the mass density, ω the specific free energy, C Hooke's tensor and

ωp(α i ) = 1 3 i C i α i : α i , (i = 1, 2
) the energy stored related to the hardenings.

The state laws become

σ = ρ ∂ω ∂ε ; Λε p = -ρ ∂ω ∂εp = σ ; X i = -ρ ∂ω ∂α i , i = 1, 2 (5)
The dissipation potential F , is different from the yield criterion f (i.e. non-associated plasticity [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF])

F (σ, X) = (σ -X) eq -σ 0 + 3γ 4C X 1 : X 1 ( 6
)
where C is the hardening modulus, and γ the inverse of a characteristic strain associated with the nonlinearity of the hardening. The growth laws read

εp = λ ∂F ∂σ = λ 3 2 (s -X) (σ -X) eq (7) α1 = λ ∂F ∂X 1 = -εp + λ 3γ 2C 1 X 1 ; α2 = λ ∂F ∂X 2 = -εp (8)
The cumulated plastic strain rate ṗ is deduced from energy equivalence (σ -X) eq ṗ = (σ -X) : εp

and corresponds to the plasticity multiplier ṗ = λ = 2 3 εp : εp [START_REF] Chaboche | Constitutive equations for cyclic plasticity and cyclic viscoplasticity[END_REF] The intrinsic dissipation is then expressed as

D i = σ : ε -ρ ∂ω ∂A i : Ȧi
where A i are all the internal variables except the plastic strain tensor so that

D i = (σ 0 + 3γ 2C X 1 : X 1 ) ṗ (11)
Under isothermal conditions, the supplied mechanical work is decomposed into

• the free energy: ρω = ωe + ωs, itself made of:

the recoverable elastic energy: ωe = 1 2 εe : C : εe and the stored energy:

ωs(α) = 1 3 (C 1 α 1 : α 1 + C 2 α 2 : α 2 ) • and the dissipated energy D 1 , which is obtained by integrating σ 0 + 3γ 2C 1 X 1 : X 1 ṗ
over the load history.

Energetic calculations

The analysis of energy balance is carried out on a tensile loading. When a cyclic test is simulated on the axisymmetric crack model, Abaqus straightforwardly delivers a series of data:

the so-called plastic work ψp through ELPD output variable [START_REF]Simulia: Abaqus/Strandard output variable identifiers (section 4.2.1). Abaqus 6.10 Analysis User's Manual[END_REF], which gives access to the total plastic work, integrated over the whole model, -the total strain energy ω eT through ELSE output variable [START_REF]Simulia: Abaqus/Strandard output variable identifiers (section 4.2.1). Abaqus 6.10 Analysis User's Manual[END_REF], which is also integrated over the whole model.

Consequently, other data have to be determined:

the stored energy ψs, including the part stored by hardening ω h and the energy related to the residual stresses ωs. -The dissipated energy D 1 .

The ω h is obtained from an Abaqus UVARM user subroutine allowing to integrate over the entire model the energy stored by the hardenings

ω h = 2π n i=1 3 4C 1 X i 1 : X i 1 + 3 4C 2 X i 2 : X i 2 r i S i (12)
where S i is the area of an element, r i the radius of its centroid, and n the total number of elements in the model. The dissipated energy D 1 is obtained by subtracting ω h from the plastic work ψp:

D 1 = ψp -ω h .
To determine the stored energy associated with the presence of residual stresses ωs, an elastic unloading must be performed at each increment [START_REF] Aravas | On the calculations of the strored energy of cold works[END_REF]; the stored energy (i.e. residual part of ω eT ) is that associated with the residual stresses. The recoverable elastic energy ωe is then equal to ω eTωs. The free energy ρω is finally equal to ψs + ωe.

Working along similar lines, it is possible to conduct an energetic analysis of the 5-bar model. A Finite Element calculation is carried out for each level of damage (i.e. for strains ranging from 0.2% to 0.5%). The energy stored in hardenings is then calculated as

ω h = n i=1 3 4C 1 X i 1 : X i 1 + 3 4C 2 X i 2 : X i 2 S i (13)
It is worth noting that there is no energy related to residual stresses in this model. Since the different bars are independent of each other and the stress field in each bar is homogeneous, no residual stresses are generated in the case of a tensile loading followed by an elastic unloading. For the axisymmetric crack model, residual stresses are related to the plastic zone surrounding the crack tip.

The energy balance, per unit area, is shown in Figure 12. An equivalence between the two models is observed, thus allowing a constitutive model based on the 5-bar model to be developed. The 5-bar model introduced in Section 3.3 helped us understand the strain mechanisms of FGCI as well as the flake graphite opening and closure conditions from a macroscopic point of view. Based on this model, a one-dimensional two-branch rheological model displayed in Figure 13 is proposed in this section. The model consists of one branch that describes the matrix, i.e. the undamaged part of the material (the right branch, called the matrix branch) and a second branch corresponds to the areas surrounding the debonded graphite flakes (the left branch, called the crack branch), such that the latter, depending on the loading conditions may or may not transmit a load. When the flakes are open, the crack branch must remain inactive and the entire load is carried by the matrix branch. Conversely, when the flakes are closed, the load is carried by both branches.

The case considered hereafter is that of a non-constrained medium, with a preexisting crack, whose length corresponds to a specific level of damage D. During a tensile loading, applying a weight (1 -D) to several coefficients of the constitutive model is sufficient to bring about the damaged material behavior. During a compressive loading, the behavior is that of the undamaged material since the crack is closed. A weight D is then applied to the same coefficients in the crack branch in order to equate the sum of the stresses in the two branches and the stress of the undamaged material.

Under cyclic loadings, the flakes opening and closure conditions need to be defined. They are schematically shown in Figure 13 by the element with two horizontal bars on the crack branch. The deactivation condition is the same as that described in Section 3.2, namely, the crack closes at the same strain level as that for which it had opened. The opening condition remains to be defined. In the present rheological model, this condition is quite simple. No tensile stresses can be carried by the crack branch. The opening (and therefore deactivation) condition of the crack branch corresponds to the instant when the stress becomes positive. During this period of deactivation, the state variables do not change and the elastic strain is equal to zero.

In order to demonstrate which coefficients should be multiplied by D and (1 -D), the same case is studied, namely, that of an unconstrained medium with a pre-existing crack submitted to a monotonic compressive loading. The overall stress is equal to the sum of the stresses in the crack branch, denoted by σc, and that of the matrix branch, denoted by σm. In the elastic case, the Young's modulus is affected by damage and the overall stress reads

σ = σc + σm = EDε + E(1 -D)ε = Eε (14)
In the viscous-elastic case, only Norton's law K parameter is affected, thereby leading to the overall viscous stress

σv = σvc + σvm = KDε 1/m v + K(1 -D)ε 1/m v = Kε 1/m v (15)
Last, for the elastoplastic case with nonlinear kinematic hardening, only the hardening modulus C is influenced by damage, leading to the overall plastic stress

σ = σc + σm = CD γ 1 -e -γεp + C(1 -D) γ 1 -e -γεp = C γ 1 -e -γεp (16)
The thermodynamic framework chosen to describe in one dimension the undamaged material is defined below. This framework is that of generalized standard materials [START_REF] Germain | Cours de mécanique des milieux continus[END_REF], with the unified elastoviscoplastic model proposed by Chaboche [START_REF] Chaboche | Description phénoménologique de la viscoplasticité cyclique avec endommagement[END_REF], which takes into account viscous effects present at high temperatures. The state potential is the free energy density

ρω = 1 2 (ε -ε th -εvp)E(ε -ε th -εvp) + i 1 2 C i α 2 i ( 17
)
where ε th is the thermal strain and εvp the viscoplastic strain. The state laws are defined as

σ = ρ ∂ω ∂ε = Eεe ; X i = -ρ ∂ω ∂α i = -C i α i (18)
The first kinematic hardening X 1 is assumed to be nonlinear, the second, X 2 , being linear. The function F is here defined by

F = σ - i X i -σ 0 + γX 2 1 2C 1 = f + γX 2 1 2C 1 (19)
In viscoplasticity, the viscoplastic multiplier λ is obtained via Norton's law λ = f K m . The growth laws are determined by applying the generalized normality rule

εvp = λ ∂F ∂σ = f K m ; α1 = λ ∂F ∂X 1 = -εvp(1 + γ) ; α2 = λ ∂F ∂X 2 = -εvp (20)
Finally, the cumulative viscoplastic strain rate is deduced from the energy equivalence (Equation ( 9))

ṗ = λ = f K m = εvp (21)

Numerical integration of the constitutive equations

The 1D two-branch constitutive model introduced in Section 4 was implemented using an implicit integration scheme [START_REF] Ortiz | Accuracy and stability of integration algorithms for elastoplastic constitutive relations[END_REF] in ZeBuLoN [START_REF]Z-Set: Non-linear material and structure analysis suite[END_REF]. The implicit scheme is a generalized middle point θ-method [START_REF] Lemaitre | Engineering Damage Mechanics[END_REF]. It consists of assessing the rate of change of variables v i at time t + ∆t such that

vi (t + θ∆t) = v i (t + ∆t) -v i (t) ∆t = ∆v i ∆t (22)
The quantities vi are estimated from the growth laws with physical quantities expressed at time t + θ∆t. This scheme is a generalization of the explicit scheme (i.e. θ = 0) and the fully implicit Euler method (i.e. θ = 1). In the present case, choosing θ = 1 is sufficient. The solution of the numerical problem in structural mechanics consists of calculating the stress increment ∆σ induced by a total strain increment ∆ε prescribed during ∆t.

Hereafter and unless explicitly stated, variables v i are always expressed at time t + θ∆t, all parameters of the material being measured at temperature T (t + θ∆t) and ∆v i is the increment of v i between times t and t + ∆t. The first order approximation v i ≡ v i (t + θ∆t) = v i (t)+θ∆v i is used, where v i (t) is assumed to be known. To further simplify the notations, in what follows the integration scheme is described on the undamaged material. [START_REF] Prager | A new method of analyzing stress and strains in work hardening plastic solids[END_REF] It is written as

J =     

Stress increment

The stress increment ∆σ is calculated from the elastic strain ∆εe when Newton's method has converged. The variation of the elastic modulus E must also be taken into account. The implicit integration scheme leads to ∆σ = σ (t + θ∆t) ∆t [START_REF] Russell | Finite element simulation of the microstructure of gray cast iron[END_REF] which is developed as

∆σ = (Eεe)∆t = E∆εe + ∆Eεe (31)

Consistent tangent operator

It is useful to define the consistent tangent matrix (here a scalar) related to the integration scheme to optimize the global convergence and avoid too small time increments. This operator is defined by L = ∂∆σ ∂∆ε [START_REF]Z-Set: Non-linear material and structure analysis suite[END_REF] It can be shown that [START_REF] Besson | Implementation of material constitutive equations in finite element codes[END_REF] L = EJ * ee (33)

where J * ee is defined from the inverse of the Jacobian J -1 . It is the block that corresponds to the subscripts related to the elastic strain εe.

Specific implementation in ZeBuLoN

The solution of the two branches of the rheological model is sought in parallel with a unique total strain increment ∆εe, and by assigning factors D or (1 -D), as previously defined, to the model parameters. The principle of the algorithm, schematically shown in Figure 14, is the following (subscripts c correspond to the crack branch and subscripts m to the matrix branch):

-The crack branch is deactivated if the stress is greater than or equal to zero.

This deactivation occurs at a strain level ε t called transition strain. The crack branch is activated when the total strain reaches ε t . -When the branch is deactivated, the internal variables do not change and the elastic strain increment and stress increment, respectively denoted by ∆εec and ∆σc, are equal to zero. The total stress increment ∆σ is then equal to the stress increment in the matrix branch ∆σm. -When the crack branch is activated, the total stress is the sum of the contributions of the two branches: ∆σ=∆σm+∆σc.

Parameter identification

To identify the parameters of the model, two steps are needed. First, the growth of the damage variable is determined from the monotonic tests. Then, for a set evolution, the other parameters are tuned. In both cases, ZeBuLoN software with a Levenberg-Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] is used. The identification method thus involves the optimization of a cost function that evaluates, for a given set of parameters, the correlation between the simulation and the experiment. The error is then expressed in terms of stresses.

Damage growth

Using the monotonic compressive test for each temperature, a first set of parameters related to the undamaged material is tuned. The damage growth as a function of the applied stress is then approximated by a second order polynomial whose coefficients become the parameters of a new identification problem. The Fig. 14: Algorithmic implementation of the rheological model latter involves simulating a monotonic tensile test using the previously identified parameters and acting only on the damage variable. The procedure is shown in Figure 15 at room temperature. The damage growth law is easy to use in nonisothermal conditions as it assumes a linear interpolation between two temperatures. For instance, if the damage growth law is defined for 300 • C and 400 • C, for a given stress, the damage at 350 • C is a linear interpolation of the damage value at 300 • C and 400 • C. The resulting polynomial coefficient values as a function of the temperature are summarized in Table 3. The calibration of a new set of parameters based on cyclic tests is then performed, for which the previously obtained growth law of the damage variable must be set. This new identification is necessary because the previous one was conducted on monotonic tests with a poor description of the viscous behavior.

The dwell phases of the cyclic tests enable for a robust identification of the viscous parameters. Elastic modulus and hardening modulus of the linear hardening X 2 are identified "manually" on the linear parts of the monotonic tension tests and the five other parameters (Yield stress, two non linear kinematic hardening terms, two viscosity terms) are identified all together on the cyclic test. The same procedure is carried out for the different temperatures. Results are presented in Table 4. There is a good description of the behavior under different loading conditions. The loss of stiffness due to tensile loadings and stiffness recovery due to compressive loadings are accurately represented for all the investigated temperatures. The stress relaxation is also very well described, even at high temperatures. For such temperatures (i.e. over 500 • C) the description of the inelastic behavior is less accurate than for lower temperatures. This effect is probably due to the use of a relatively simple viscous (i.e. Norton) potential in the constitutive model.

The shift of the inflexion point in the negative strain area is also exacerbated at high temperatures. This phenomenon highlights the limitations of the proposed rheological model and the associated flake opening/closure conditions. Owing to the presence of graphite flakes, the behavior of FGCI strongly differs under tensile and compressive loadings. This asymmetry makes the development of a constitutive model very challenging. The simplest approaches based on the modification of the yield surface do not account for complex phenomena involved during the debonding and closure of the flake/matrix interfaces during cyclic loadings.

A rheological model composed of two-dimensional elements with varying sections has helped us define the basis of a constitutive model, which is established for the time being in one dimension, under tensile/compressive loadings. This model accurately represents the principal mechanisms of FGCI. A comparison of this model and the cyclic tests at temperatures ranging from 20 • C and 600 • C shows that it is able to describe cyclic loadings at different strain rates, with dwell phases under both tensile and compressive loadings.

It now remains to extend the model proposed to a three-dimensional representation in order to simulate the behavior of brake disks subjected to complex loadings. This issue is addressed in a companion paper [START_REF] Augustins | A constitutive model for flake graphite cast iron automotive brake disks. Induced anisotropic damage model under complex loadings[END_REF].

Fig. 1 :

 1 Fig. 1: Simulation of a cyclic test with an elastoplastic model and Drucker-Prager criterion

Fig. 2 :

 2 Fig. 2: Micrography of flake graphite cast iron

Fig. 3 :

 3 Fig. 3: Monotonic tests performed at different temperatures for the considered FGCI

Fig. 4 :Fig. 5 :

 45 Fig. 4: Cyclic test at room temperature on the considered FGCI

Fig. 6 :

 6 Fig. 6: Simulation of the compression test with the identified elastoplastic model

Fig. 7 :

 7 Fig. 7: Axisymmetric crack model

Fig. 8 :

 8 Fig. 8: Comparison of the axisymmetric crack model and the cyclic test

Fig. 9 :

 9 Fig. 9: Damage growth law

Fig. 10 :

 10 Fig. 10: Five-bar rheological model

Fig. 11 :

 11 Fig. 11: Prediction of the 5-bar model and experimental data of the cyclic test at 20 • C

( a )Fig. 12 :

 a12 Fig. 12: Various energetic contributions in a tensile test when two different models are considered

Fig. 13 :

 13 Fig. 13: Rheological model consisting of two branches: the left part is related to the flakes and the right part defines the matrix

Fig. 15 :

 15 Fig. 15: Damage identification procedure based upon monotonic tests

4. 5

 5 Results and discussion The simulations of the cyclic tests with the tuned parameters are shown in Figure 16 for four different temperatures. Simulations between 20 • C and 400 • C are not presented because the behavior of the considered cast iron is very similar for temperatures up to 300 • C [Unpublished test report, Exova, 2011].

(a) 20 Fig. 16 :

 2016 Fig. 16: Comparison of the model predictions and the cyclic tests at four different temperatures

Table 1 :

 1 Chemical composition (wt.%) of the studied FGCI

	C	Si	Mn	S	P	Ni	Cr	Mo	Cu	Sn	Ti	Fe

3.4 2.45 0.40 0.15 0.15 0.25 0.20 0.10 0.60 0.10 0.03 Bal.

Table 2 :

 2 Matrix parameters -20 • C

	E (GPa) ν	σy	C 1 (MPa) γ 1	C 2 (MPa)
	119	0,28 181 230000	1634 17300

  

	∂Rε e	∂Rε e	∂Rε e	∂Rε e	
	∂∆εe	∂∆α1	∂∆α2	∂∆p
	∂Rα 1 ∂∆εe ∂Rα 2 ∂∆εe ∂R p	∂Rα 1 ∂∆α1 ∂Rα 2 ∂∆α1 ∂R p	∂Rα 1 ∂∆α2 ∂Rα 2 ∂∆α2 ∂R p	∂Rα 1 ∂∆p ∂Rα 2 ∂R p ∂∆p	    	(29)
	∂∆εe	∂∆α1	∂∆α2	∂∆p

Table 3 :

 3 Polynomial coefficients

	Température ( • C)	a 0	a 1	a 2
	20 • C	-1.15E-05 6.23E-03 4.46E-01
	200 • C	-1.18E-05 5.25E-03 1.95E-01
	400 • C	-1.18E-05 5.25E-03 1.95E-01
	500 • C	-2.20E-05 7.55E-03 2.70E-01
	600 • C	-4.52E-05 1.19E-02 3.88E-01
	800 • C	-1.39E-04 2.07E-02 3.63E-01

Table 4 :

 4 Model parameters

	Température ( • C)	E (GPa) σy (MPa)	K	m	C 1 (M P a) γ	C 2 (M P a)
	20 • C	119	158	8.5 110	151300	1010 11000
	200 • C	119	145	9.5 110	127900	1140 11000
	400 • C	110	61	9	306	124400	1210 9000
	500 • C	110	10	6	652	86000	2940 8500
	600 • C	87	1	3.7 1295 10000	5000 7800
	800 • C	62	1	4.7 304	500	3940 100

Acknowledgements This work has been performed within the Executive Doctorate in Industry framework, and has been supported by PSA Peugeot Citroën.

Residual vector

Newton's method numerically minimizes the local residual defined as

The chosen integration variables are {εe, α 1 , α 2 , p}, so that the residual vector components are expressed as

Jacobian matrix

The nonlinear problem to be solved is a 4 × 4 dimension system. The solution using Newton's method requires knowing the Jacobian matrix J of the system, determined here analytically, and defined by components such that


L Augustins
email: louis.augustins@mpsa.com

R Billardon
email: rene.billardon@safranmbd.com

F Hild
email: hild@lmt.ens-cachan.fr

F Hild Lmt

Constitutive model for flake graphite cast iron automotive brake discs Induced anisotropic damage model under complex loadings

Keywords: Flake graphite cast iron, Constitutive model, Brake discs, Anisotropic damage, Damage deactivation

The present paper details an elasto-viscoplastic constitutive model for automotive brake discs made of flake graphite cast iron. In a companion paper [2], the authors proposed a one-dimensional setting appropriate for representing the complex behavior of the material (i.e. asymmetry between tensile and compressive loadings) under anisothermal conditions. The generalization of this 1D model to 3D cases on a volume element and the associated challenges are addressed. A direct transposition is not possible and an alternative solution without unilateral conditions is first proposed. Induced anisotropic damage and associated constitutive laws are then introduced. The transition from the volume element to the real structure and the numerical implementation requires a specific basis change. Brake disc simulations with this constitutive model show that unilateral conditions are needed for the friction bands. A damage deactivation procedure is therefore defined.

Introduction

The present paper is the second part of a study devoted to the simulation of the thermomechanical behavior of automotive brake discs made of flake graphite (grey) cast iron (i.e. FGCI). An overview of the industrial and scientific issuesmainly related to the thermomechanical fatigue design of automotive brake discs and to the correct description of the properties of their constitutive material -is proposed in the first part [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF].

The model developed earlier is based on an elasto-viscoplastic postulate with the formalism proposed by Chaboche [START_REF] Chaboche | Description phénoménologique de la viscoplasticité cyclique avec endommagement[END_REF], and on a precise description of the graphite flakes opening and closure phases, which are responsible for the strong dissymmetry between tensile and compressive response of FGCI [START_REF] Gilbert | An evaluation of the stress-strain properties of flake graphite cast iron in tension and compression[END_REF]. Such a model allows for an accurate simulation of the cyclic behavior of FGCI in isothermal or anisothermal conditions for a one-dimensional tension/compression loading.

Automotive brake discs are subjected to severe thermomechanical loadings [START_REF] Dufrenoy | Etude du comportement thermomécanique des disques de freins vis à vis des risques de défaillance. application au domaine ferroviaire[END_REF][START_REF] Abdi | Effet de la distribution du flux dans un disque de frein sous sollicitations thermomécaniques[END_REF][START_REF] Vallet | Etude de la fissuration d'un disque de frein à partir de l'analyse de son comportement thermomécanique[END_REF] leading to cyclic plasticity in the grooves and potential crack initiation [START_REF] Vallet | Etude de la fissuration d'un disque de frein à partir de l'analyse de son comportement thermomécanique[END_REF][START_REF] Dufrenoy | Damage mechanisms and thermomechanical loading of brake discs[END_REF]. Three-dimensional finite element calculations are necessary to assess the physical quantities that control the fatigue life of such a structure. As a one-dimensional description of the material was proposed [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF], a 3D extension applicable under complex loadings needs to be formulated.

In the present paper, the generalization of the 1D model to 3D cases on a volume element is addressed as a first step. The transition from this volume element to the real structure and the corresponding numerical implementation are then described.

2 From a 1D description to a 3D constitutive model

Extension of the rheological model

A complete 3D braking simulation requires the 1D model proposed by the authors [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF] to be extended to 3D situations, which are valid under complex loading histories. This generalization is not straight forward:

-Superposition of 1D models composed of two 1D branches (i.e. one for the damaged matrix and one for the debonded flakes) in the three principal directions. When the stress tensor has more than one component, reducing the three-dimensional formulation to a one-dimensional expression proves impossible. Since there is no reason for the stress tensor to have only one component in a 3D structure subjected to complex loads, this option must be disregarded.

-Extension of the 1D model [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF] to a constitutive law containing two 3D branches. The objective of the present work is to implement the constitutive model in the finite element solver Abaqus using a user material subroutine (UMAT [START_REF]Simulia: Abaqus/Strandard subroutines[END_REF]). The input data to this type of routine are the total strain tensor increment, and the internal variables at the beginning of the increment. If, for instance, a tensile loading in direction 1 is considered, the global strain increments will have three non zero terms, namely, ∆ 11 , and ∆ 22 = ∆ 33 = -ν∆ 11 . The Poisson's ratio ν is differently affected by damage in the two branches (D and (1 -D)). As a consequence, there is an incompatibility between the transverse strains ( 22 and 33 ) of the two branches. This option entails an over or underestimation of the stresses in one of the branches.

It is therefore worth asking whether managing the opening/closure mechanisms is crucial to a specific application on brake discs, and to what extent accurately simulating the behavior of the disc without this phenomenon is possible. To answer this question, it is necessary to study the stress/strain paths in the critical areas of the disc, and to understand the conditions under which the flakes open and close within the structure. A 3D constitutive model that can take into account flake debonding induced by loading and the associated anisotropic damage description is required. This model will represent flake opening but not its closure, a pattern that corresponds to the matrix branch in the 1D model [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF], and will make it possible to say whether debonded flakes fulfill the conditions under which they will close on a real structure.

2.2 3D Constitutive model without flake closure

Anisotropic damage

Any damage orientation distribution can be theoretically considered, resulting in an anisotropic damage description [START_REF] Leckie | Tensorial nature of damage measuring internal variables[END_REF][START_REF] Lubarda | Damage tensors and crack density distribution[END_REF]. However, a preferred orientation of microcracking, more or less perpendicular to the maximum principal stress, is more often the case [START_REF] Lemaitre | Mécanique des Matériaux Solides[END_REF]. The simplest way to represent damage without referring to a particular basis (i.e. intrinsic writing) is to resort to even order tensors [START_REF] Leckie | Tensorial nature of damage measuring internal variables[END_REF].

The most general form that allows the basis change invariance for any damage distribution to be satisfied requires the use of an eighth order tensor [START_REF] Chaboche | Description phénoménologique de la viscoplasticité cyclique avec endommagement[END_REF]. Many authors prefer simpler models and favor fourth order [START_REF] Chaboche | Le concept de contrainte effective appliquée à l'élasticité et à la viscoplasticité en présence d'un endommagement anisotrope[END_REF][START_REF] Ortiz | A constitutive theory for the inelastic behavior of concrete[END_REF][START_REF] Krajcinovic | Damage Mechanics[END_REF] or second order [START_REF] Murakami | A continuum theory of creep and creep damage[END_REF][START_REF] Kachanov | Continuum model of medium with cracks[END_REF][START_REF] Cordebois | Endommagement anisotrope en élasticité et plasticité[END_REF][START_REF] Chow | An anisotropic theory of continuum damage mechanics for ductile fracture[END_REF][START_REF] Halm | A model of anisotropic damage by mesocrack growth; unilateral effect[END_REF][START_REF] Desmorat | Non local anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF] tensors, by which complex anisotropy problems and simple damage symmetries such as orthotropy are described.

In this study, damage is induced by the loading and corresponds to graphite flakes debonding, which is mostly oriented perpendicular to the loading direction. Figure 1 illustrates this effect by showing an idealized distribution composed of planar flakes. Given that the loading on a brake disc is known, the selected option is to represent damage by a second order tensor via the introduction of the degradation of the elastic behavior [START_REF] Cordebois | Endommagement anisotrope en élasticité et plasticité[END_REF]. The form of the damaged material compliance tensor S is derived from Refs. [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF][START_REF] Forquin | A probabilistic damage model of the dynamic fragmentation process in brittle materials[END_REF], itself derived from Fabrikant's work [START_REF] Fabrikant | Complete solutions to some mixed boundary value problems in elasticity[END_REF]. The authors studied unidirectional crack distributions perpendicular to the eigen directions in an isotropic elastic matrix. If only the distribution perpendicular to the 1-axis is activated, the following results are obtained

E 1 E 0 = (1 -D) , E 2 E 0 = E 3 E 0 = ν 23 ν 0 = 1
(1)

ν 12 ν 0 = ν 13 ν 0 = E 1 E 0 , G 12 G 0 = G 13 G 0 = (1 -D) α (2)
where subscripts 0 correspond to the initial isotropic matrix, E, G and ν to Young's modulus, shear modulus and Poisson's ratio, and D to the damage variable. In Equation (2), a power α is introduced and is equal to 0.5 [START_REF] Petry | Caractérisation et modélisation du comportement et de l'endommagement d'alliages métalliques sur une grande plage de température[END_REF]. The superposition principle then makes it possible to define a compliance tensor for the damaged material involving three scalar damage variables associated with three crack distributions

S = C -1 = 1 E           1 1-D1 -ν -ν 0 0 0 1 1-D2 -ν 0 0 0 1 1-D3 0 0 0 sym. 1+ν √ (1-D1)(1-D2) 0 0 1+ν √ (1-D2)(1-D3) 0 1+ν √ (1-D1)(1-D3)           (3)
where C is Hooke's tensor of the damaged material. In the sequel, the notation C refers to a fourth-order tensor. This relationship is written in the eigen stress basis, which corresponds to the damage basis. As shown in the sequel, this setting proves sufficient for brake discs because damage anisotropy is induced by loading whose eigen directions do not change with time. This result is true in the specific problem discussed in the present work because of several features. In the most general case of a disc subjected to acceleration/breaking cycles, body forces induce shear stresses. The friction between the pad and the disc also results in a shear stress on the friction band. Models were proposed to solve such cases with an orthotropic damage growth rule applied to variable principal directions [START_REF] Skrzypek | Application of the orthotropic damage growth rule to variable principal directions[END_REF]. However, in the case of automotive grey cast iron brake discs, the stress levels generated by the heat flux on the friction band are significantly higher than the aforementioned stresses that are neglected herein.

Effective stress, yield criterion

The effective stress tensor σ is defined by the strain equivalence principle [START_REF] Chaboche | Sur l'utilisation des variables d'état internes pour la description du comportement viscoplastique et de la rupture par endommagement[END_REF][START_REF] Lemaitre | Aspect phénoménologique de la rupture par endommagement[END_REF] and corresponds to the stress that must be applied on an undamaged elastic volume element to obtain the same strain as that of a damaged element. This definition makes explicit the effective stress versus stress tensor σ

σ = C : S : σ ( 4
)
where C is Hooke's tensor for the undamaged material.

The equivalent Von Mises stress, for the undamaged material, is then defined by

σeq = 3 2 s : s (5)
where s is the deviatoric stress tensor

s = K : σ ( 6
)
where K is the projector operator onto the deviatoric space

K = I - 1 3 1 ⊗ 1 ( 7
)
1 and I the second and fourth order identity tensors, respectively. The equivalent stress of the damaged material is defined on the assumption that the first yielding state of the damaged material under a stress σ is equivalent to that of the undamaged material under a stress σ

σeq = 3 2 s : s = 3 2 (K : σ) : (K : σ) (8)
with M = K : C : S and given Equation ( 4), it leads to

σeq = 3 2 (M : σ) : (M : σ) (9)

Thermodynamic framework

In this section, the thermodynamic framework for the undamaged material is defined within the generalized standard material setting [START_REF] Germain | Cours de mécanique des milieux continus[END_REF] with a unified viscoplasticity model [START_REF] Chaboche | Description phénoménologique de la viscoplasticité cyclique avec endommagement[END_REF]. The state potential is defined by the free energy density

ρω = 1 2 ( -th -vp) : C : ( -th -vp) + ωp(α) (10)
where ωp(α) = 1 3 Cα : α is the energy stored in hardenings. The dual dissipation potential is chosen as

φ = K m + 1 f K m+1 (11)
where K and m are Norton's viscous parameters, and f the yield surface.

In order to simplify these equations, only one nonlinear kinematic hardening term is considered in the present section. This hypothesis implies a flow potential different from f (i.e. non-associated plasticity [START_REF] Lemaitre | Mécanique des Matériaux Solides[END_REF])

F = f + 3γ 4C X : X ( 12
)
with f = (σ -X) eq -σ 0 ( 13
)
where C is the hardening modulus, and γ the inverse of a characteristic strain associated with the nonlinearity of the hardening. Applying the normality rule, the growth laws read

˙ vp = ∂φ ∂σ = λ ∂F ∂σ = f K m 3 2 (s -X) (σ -X) eq ( 14
) α = ∂φ ∂X = λ ∂F ∂X = -˙ vp -ṗγα (15)
The cumulated plastic strain rate ṗ is deduced from energy equivalence (s -X) : ˙ vp = (σ -X) eq ṗ ( 16)

and corresponds to the viscoplastic multiplier

ṗ = λ = f K m = 2 3 ˙ vp : ˙ vp (17)

State and growth laws

In this subsection, state and growth laws are introduced for the damaged material on the assumption of a state decoupling between plasticity and damage [START_REF] Marquis | Phénoménologie et thermodynamique -Application au couplage entre la thermoélasticité, la plasticité, le vieillissement et l'endommagement[END_REF]. The thermodynamic potential is written in an analogous way as Equation ( 10), with the difference that C is substituted for

C ρω = 1 2 ( -th -vp) : C : ( -th -vp) + ωp(α i ) (18)
where the viscoplasticity related terms remain unchanged. The state laws then read

σ = ρ ∂ω ∂ (19)
X i = -ρ ∂ω ∂α i (20)
The first back-stress X 1 is nonlinear and the second one X 2 is linear. The thermodynamic force Y associated with D is defined as

Y = -ρ ∂ω ∂D (21)
It should be emphasized that the damage variable is introduced in the model to simulate the stabilized behavior of the material and not to calculate continuous damage right up to material failure. The growth of D must therefore be bounded by a maximum value -potentially reached after a first tensile phase -associated with the volume fraction of debonded graphite flakes. Under these conditions, in order to have more freedom in the definition of the damage growth law, the latter is written in a non standard form detailed in Ref. [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF]. The flow potential F now reads

F = σ - i X i eq -σ 0 + 3γ 4C 1 X 1 : X 1 - 4 9 C 2 1 α 1 : α 1 (22) with σ - i X i eq = 3 2 M : σ - i X i : M : σ - i X i = J (23)
and

X i = -2 3 Cα i .
The growth laws for the damaged material become

˙ vp = ∂φ ∂σ = λ ∂F ∂σ = f K m 3 2 M T : M : σ -M T : i X i σ - i X i eq = f K m n (24)
with

n = 3 2 M T : M : σ -M T : i X i σ - i X i eq α1 = ∂φ ∂X 1 = λ ∂F ∂X 1 = - f K m      3 2 M : σ - i X i σ - i X i eq + γα 1      = f K m (n X -γα 1 ) (25)
and

n X = - 3 2 M : σ - i X i σ - i X i eq α2 = ∂φ ∂X 2 = λ ∂F ∂X 2 = - f K m      3 2 M : σ - i X i σ - i X i eq      = f K m n X (26)
"Auxiliary" variables are introduced via the use of "effective" viscoplastic strain and cumulated viscoplastic strain rate that help, in the anisotropic case, to significantly simplify the numerical implementation of the model. The "effective" viscoplastic strain is defined as

˙ vp = ∂φ ∂ σ = λ ∂F ∂ σ = f K m 3 2 s - i X i σ - i X i eq (27)
The "effective" cumulated viscoplastic strain rate ṗ is deduced from energy equivalence

s - i X i : ˙ vp = σ - i X i eq ṗ ( 28
) such that ṗ = λ = f K m = 2 3 ˙ vp : ˙ vp (29)
Equation [START_REF] Lemaitre | Eng. Damage Mech[END_REF] shows that the viscous stress is a function of the "effective" cumulated viscoplastic strain rate rather than the cumulated viscoplastic strain rate.

Numerical integration

The constitutive model was implemented using an implicit integration scheme [START_REF] Ortiz | Accuracy and stability of integration algorithms for elastoplastic constitutive relations[END_REF] in an Abaqus UMAT routine [START_REF]Simulia: Abaqus/Strandard subroutines[END_REF]. It is a generalized middle point θ-method [START_REF] Lemaitre | Eng. Damage Mech[END_REF] with θ = 1 (i.e. fully implicit Euler method). The solution to the numerical problem in structural mechanics consists of calculating the stress increment ∆σ induced by a total strain increment ∆ prescribed during ∆t.

Newton's method numerically minimizes the local residual defined as

{R} = {∆v i } -∆t{ vi (t + θ∆t)} (30)
The chosen integration variables are { e, α 1 , α 2 , p} so that the residual vector components are expressed as

R e = ∆ e -(∆ -∆pn) (31)
R p = ∆p -∆t J -σ 0 K m ( 32
)
Rα 1 = ∆α 1 -∆p (n X -γα 1 ) (33)
Rα 2 = ∆α 2 -∆pn X (34)
The nonlinear problem to be solved is a 19×19 dimension system. The solution using Newton's method requires knowing the Jacobian matrix J of the system, which is determined analytically and defined by components such that

J ij = ∂R i ∂∆v j (35)
It is written as

J =       ∂R e ∂∆ e ∂R
      (36)
In order to express each term of the matrix, several quantities are defined

N = ∂n ∂σ = 1 J 3 2 M T : M -n ⊗ n (37)
S = ∂n ∂X = - 1 J 3 2 M T + n ⊗ n X (38) N X = ∂n X ∂σ = - 1 J 3 2 M -n X ⊗ n (39) S X = ∂n X ∂X = 1 J 3 2 K -n X ⊗ n X (40)
The terms derived from Equation [START_REF] Lubarda | Damage tensors and crack density distribution[END_REF] read

∂R e ∂∆ e = I + θ∆pN : C (41) ∂R e ∂∆p = n (42) ∂R e ∂∆α 1 = - 2 3 C 1 θ∆pS (43) ∂R e ∂∆α 2 = - 2 3 C 2 θ∆pS (44)
When the viscoplastic flow is activated (J > σ 0 ), the terms derived from Equation [START_REF] Marquis | Phénoménologie et thermodynamique -Application au couplage entre la thermoélasticité, la plasticité, le vieillissement et l'endommagement[END_REF] become

∂R p ∂∆ e = - m K θ∆p J -σ 0 K m-1 n : C ( 45
)
∂R p ∂∆p = 1 ( 46
)
∂R p ∂∆α 1 = m K θ∆p J -σ 0 K m-1 2C 1 3 n X ( 47
)
∂R p ∂∆α 2 = m K θ∆p J -σ 0 K m-1 2C 2 3 n X (48)
and for hardening variables (see Equations. ( 33) and ( 34))

∂Rα 1 ∂∆ e = -θ∆pN X : C (49) ∂Rα 1 ∂∆p = -n X + γα 1 ( 50
)
∂Rα 1 ∂∆α 1 = I + θ∆p 2C 1 3 S X + γI (51) ∂Rα 1 ∂∆α 2 = 2 3 C 2 θ∆pS X (52)
and

∂Rα 2 ∂∆ e = -θ∆pN X : C (53) ∂Rα 2 ∂∆p = -n X ( 54
)
∂Rα 2 ∂∆α 1 = 2 3 C 1 θ∆pS X ( 55
)
∂Rα 2 ∂∆α 2 = I + 2 3 C 2 θ∆pS X (56)
The stress increment ∆σ is calculated from the elastic strain increment ∆ e when Newton's method has converged. The variation of the damaged Hooke's tensor C must also be considered. The implicit integration scheme leads to

∆σ = σ (t + θ∆t) ∆t (57)
which is developed as

∆σ = C : ∆ e + ∆ C : e (58)
It is useful to define the consistent tangent matrix linked to the integration scheme in order to optimize the global convergence and avoid too small time increments. This operator is defined by

L = ∂∆σ ∂∆ (59)
It can be shown that [START_REF] Besson | Implementation of material constitutive equations in finite element codes[END_REF]

L = C : J * ee ( 60
)
where J * ee is defined from the inverse of the Jacobian J -1 . It is the block that corresponds to the subscripts linked to the elastic strain e tensor.

Although the application of non-associated plasticity (f = F ) results in a lack of symmetry of the tangent operator, no particular problems of convergence were observed on the complex 3D FE problems considered herein. Furthermore, convergence with the calculated elastoplastic tangent matrix is far better than that observed with a simple elastic one.

3 From a volume element to a structure

The constitutive model presented in Section 2.2, which is valid for monotonic loading histories (i.e. without flake closure) on a volume element, must now be extended to complex loadings in a structure. The main challenge is to deal with the damage directions, namely, the description of the damaged compliance tensor S defined in Equation (3) implies that the constitutive equations be written in the damage basis, which corresponds to the eigen basis assumed to remain unchanged with time.

In this section, the architecture of a brake disc system and the damage mechanisms on critical areas are briefly described. Then, following a presentation of the heat flux modeling, the implementation of the model is described taking into account the previous challenges. Brake disc simulations, on which the flake closure conditions are studied, are then carried out. Model parameters are those identified in Ref. [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF]. A solution is finally proposed to deal with flake opening and closure in specific areas of the disc.

Architecture of a brake disc system

Brake disc systems appeared in the 1950s to overcome the drawbacks of drum brakes. Such systems are now widely used in the automotive as well as in the rail and aeronautical industries. The system, which is shown in Figure 2, includes the following elements:

-A plain or ventilated disc, generally made of cast iron, connected to the rim by means of the hub. As a rotating part, its role is to absorb and dissipate the mechanical energy transformed into heat during braking. -Two brake pads, made of frictional material (the friction lining) mounted on a metallic support. The pads are brought into contact with the disc during braking.

-A caliper, a motionless part generally connected to the steering knuckle, on which the pads are attached. Yoke shaped, it covers an angular sector of the disc. Its role is to press the pads against the disc to slow it down. The caliper is equipped with one or more hydraulic pistons. Damage inducing fatigue cracks on a brake disc is mainly due to the presence of high thermal gradients that occur during braking. As a consequence, the disc tends to expand in the hot zones, but neighboring cold zones prevent free expansion. There are two main types of thermal gradients:

the gradient between the friction ring, which is directly affected by the heat flux, and the bowl; -the gradient between the friction bands.

The gradient between the friction ring and the bowl (see Figure 3) gives rise to "coning" [START_REF] Lee | Temperature and coning analysis of brake rotor using axisymmetric finite element technique[END_REF], namely, an axial distortion of the brake disc planar geometry, which results in the disc becoming conical in shape, owing to the different thermal expansions. The groove is thus subjected to high tensile stresses during braking leading to local plastic strains. Upon cooling, the friction ring straightens without returning to its original position in conjunction with the occurrence of compressive plasticity in the same zone. The area is thus subjected to cyclic plasticity and can lead, on severe bench tests, to circumferential cracks revealed in Figure 4(a) with dye penetrant. The gradient between the friction bands (see Figure 3) is explained by the following sequence. As breaking generally lasts only a few seconds, the friction bands are heated very quickly. The surface temperature can reach up to 600 • C in 4 or 5 s, while the center of the friction ring only reaches about 200 or 300 • C during the same time. This difference in temperature with depth induces high levels of radial and hoop compressive stresses. Upon cooling, residual tensile stresses following plastic compressive strains appear. On severe bench tests, this phenomenon leads to thermal striping [START_REF] D'cruz | Surface crack initiation in ventilated disc brakes under transient thermal loading[END_REF][START_REF] Bagnoli | Thermal fatigue cracks of fire fighting vehicles gray iron brake discs[END_REF] As a first approximation, the intensity of the heat flux experienced by the disc directly depends on the lining pressure field P , the friction coefficient µ between the disc and the linings, the wheel angular velocity ω and r the radius on which the pressure is applied: φ = µrωP . However, this heat flux is very unstationary on a brake disc due to thermoelastic instabilities (TEI) [START_REF] Cristol | Coupling between thermal localisation and friction mechanisms: heat accumulation effect[END_REF] and fading phenomenon [START_REF] Limpert | Brake design and safety[END_REF] (i.e. drop of the friction coefficient for high temperatures) so that a more global definition of the heat flux derived from an energetic framework is preferred. Mechanical stresses due to the caliper are deemed negligible with respect to the thermal stresses, namely, the pressure of the pads on the disc results from the hydraulic actuator, which rarely exceeds 80 bar, leading to compressive stresses of about a few MPa. As a consequence, the model used for braking simulations is an angular sector (Figure 2(b)) on which an axisymmetric flux is applied.

A test bench with inertia J spinning with an initial velocity ω 0 is considered on which braking with constant deceleration (i.e. the speed decreases linearly during braking) can be carried out several times. Thus, if ∆t is the braking duration, the angular velocity ω at time t reads

ω(t) = ω 0 1 - t ∆t (61)
The kinetic energy is then deduced at time t

E k (t) = 1 2 Jω(t) 2 = 1 2 Jω 2 0 1 - t ∆t 2 (62)
It is considered that all the kinetic energy is converted into heat and that 95 % of the latter is transmitted to the disc, while the remaining 5% is transmitted to the linings. The heat flux transmitted to the disc is then equal to the derivative of the kinetic energy with respect to time, weighted by a factor ρ = 0.95. On the contact area with each pad, the expression of heat flux per unit area as a function of time t then reads

φ(t) = -ρ 1 Sp dE k dt = ρ Jω 0 ω(t) Sp∆t ( 63
)
where Sp is the contact surface between the pad and the disc. For an axisymmetric flux modeling, Sp corresponds to the whole frictional surface

Sp = π(R 2 out -R 2 in ) ( 64
)
where R in et R out are the inner and outer friction surface radii.

Implementation strategy

An implementation strategy of the three-dimensional model presented in Section 2 is proposed. As already mentioned, the main difficulties pertain to the management of the damage directions when solving the problem through a UMAT routine.

Principal directions and damage basis

The prerequisite is that the damage directions do not change with time. As damage occurs in planes normal to the eigen directions, it requires that the eigen directions must also remain unchanged. In order to check that hypothesis, a braking simulation was carried out with the proposed constitutive model without damage.

The eigenvector associated with time t corresponding to the maximum eigen stress is calculated. This vector n corresponds to the normal to the plane that will experience the most part of induced damage. The traction T = σ • n is then calculated for all times. If the 1-component of T is always equal to the maximum eigen stress, it means that the eigen directions do not change with time.

Figure 5 shows for an element in the groove the eigen stresses σ I , σ II and σ III as functions of the 1-component of T . The latter exactly corresponds to σ I for positive values, then is successively equal to σ II and σ III . This means that the eigen directions do not change, but the magnitude of the eigen stresses does change for negative values. The choice made in Section 2.2.1 for the representation of anisotropic damage is therefore validated. From a numerical point of view, a major difficulty arises from the fact that the UMAT input data ∆ are defined in the global basis. The latter is necessarily different from the eigen (or damage) basis of each element. As the numerical solution must take place in the stationary damage basis, a specific procedure was developed in the UMAT routine (Figure 6):

-During the first increment, the eigen basis e i ⊗ e j is defined once and for all.

The transformation matrix P ij = e i • E j is also defined, where E j corresponds to the global basis. -The first increment is treated in the global basis, which is of no consequence as the damage is initialized to zero in all directions. -During the second increment, the strain tensor increment ∆ D and the stress tensor σ D are defined in the damage basis e i ⊗ e j such that ∆ D kl = P T ki ∆ ij P jl and σkl D = P T ki σ ij P jl .

-The problem is solved in the damage basis.

-After convergence of Newton's scheme and updating of the stress and internal variables, an inverse basis change is carried out, from σ D toward σ, in order to enable Abaqus to reach global equilibrium. The internal variables are stored in the damage basis in the UMAT routine. -From the third increment onward, the steps are those described for the second one. As shown in Equation (60), the tangent operator L is calculated as L = C : J * ee , such that it is defined in the damage basis, and denoted by L D . In order to ensure global convergence, this matrix must be defined in the global basis of the finite element model. It is therefore necessary to operate a basis change, from the damage basis toward the structure basis. On a fourth order tensor, this operation is time consuming and Bond matrices [START_REF] Auld | Acoustic Fields and Waves in Solids[END_REF]

Similarly, M and N matrices are defined as With the constitutive model previously described it is possible to carry out braking simulations in which stress/strain paths in critical areas (i.e. the groove and the friction bands) are studied. Two configurations were tested, namely, the first one is a small plain brake disc 266-mm in diameter (denoted D266) mounted on a light vehicle; the second one is a large ventilated disc 330-mm in diameter (denoted D330) mounted on a much heavier vehicle. A series of severe brakings (i.e. deceleration from 80% Vmax to 40% Vmax at 0.5g) was simulated for the two configurations. The constitutive model accounts for anisotropic damage but without flake closure. Stress/strain paths in the radial and hoop directions in the groove are shown in Figure 7(a). The response is quickly characterized by plastic shakedown. The "coning" of the disc during the first braking induces a global deformation of the friction ring and permanent strains in the groove. As a result, even at the end of the cooling phase, there is still a positive residual strain. Under these conditions, flake closure will never occur in the groove and the present constitutive model is sufficient to describe the actual behavior of that area.

M =         a 2
The strong thermal gradient on the friction bands gives rise to compressive stresses and plasticity during braking, and residual tensile stresses after cooling (Figure 7(b)). The latter ones are below the threshold (about 90 MPa) at which damage occurs [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF]. This result means that in the friction bands graphite flakes do not debond, and damage is equal to zero in all directions. Again, the present constitutive model is sufficient to describe the actual behavior of the studied area.

Outside these critical areas, no other part of the disc requires flake closure be taken into account; the flakes are either debonded or closed, and remain in that state. For this disc configuration and the considered braking conditions, the Conclusions are rather different for D330 brake disc, for which the applied thermal flux is more important, owing to the higher inertia of the corresponding vehicle. The stress/strain path on the friction bands shown in Figure 8 with the current constitutive model shows a non-physical drift of the mechanical response toward positive stresses. This phenomenon is explained by the following reasons:

the tensile residual stress after the first cooling is greater than 90 MPa, damage is therefore activated in this direction; -as a consequence, Hooke's tensor is affected by damage; -given that no closure condition or damage deactivation is considered, components of Hooke's tensor remain underestimated, resulting in an overestimation of the minimum compressive stress reached at the end of the braking phase; -the subsequent tensile stress is thus overestimated, damage increases and components of Hooke's tensor decrease; -the process intensifies for each braking and results in a totally non physical response. The solution chosen to counteract this phenomenon consists of a damage deactivation strategy. The procedure described in Section 3.6 is applied on elements for which a mechanical strain value of -0.2% is observed. This choice is made because the friction bands alone are concerned, namely, the only area where this value is reached. Moreover it is assumed that only severe brakings will cause damage in tension, resulting in large negative mechanical strains.

Damage deactivation

The damage deactivation process consists of deactivating damage instantaneously when flakes close. Consequently, it must be treated direction by direction. As shown in Ref. [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF], flake closure occurs approximately at the strain level at which the flakes opened. The strain level is easy to define with the 1D two-branch model, but there is no direct way to define it with the above 3D model. However, for the loading considered in the friction bands, initial flake debonding (i.e. for the first braking) and subsequent flake openings (i.e. for subsequent brakings) are most likely to occur when the stress shifts from negative to positive levels, at a corresponding strain level t . This condition defines activation and the deactivation condition can therefore be defined as the shift from > t to < t .

Such activation and deactivation conditions were tested with unsatisfactory results (Figure 9). Simple damage deactivation does not physically represent the mechanisms involved, namely, the material actually includes two separate phases represented by the two-branch rheological model [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF]. Deactivating damage instantaneously upon unloading thus results in non-smooth paths and non physical behaviors with a drift even more pronounced than in the undamaged case. This bilinear path is a typical inconsistency associated with damage activation/deactivation conditions. More realistic models based on continuous damage deactivation are proposed, in which microcracks close gradually [START_REF] Ganczarski | Damage deactivation of engineering materials ans structures[END_REF]. Yet the specific conditions under which flakes open and close on a 3D structure make the application of such approaches quite complex in the present case. The solution finally adopted to counteract this phenomenon is to enforce damage deactivation always at the same strain level, thereby avoiding drift. The simplest and most pragmatic choice is to deactivate damage as early as the first unloading increment. This condition makes no sense on fully reversed cyclic tests and will result in a very bad description of the overall unloading behavior. Conversely, the error is minimal on the friction band (i.e. the only area where this condition is in fact activated), because subsequent closure occurs very quickly after unloading. The corresponding stress path is shown in Figure 10 where a stabilized behavior is observed, a reasonable maximum stress estimation and a low error during unloading. The present paper describes an attempt to generalize a one-dimensional setting [START_REF] Augustins | Constitutive model for flake graphite cast iron automotive brake discs. From macroscopic multiscale models to a 1D rheological description[END_REF] to describe the thermomechanical behavior of flake graphite cast iron to a three dimensional constitutive law valid for complex loadings on a brake disc. Directly transposing the 1D model to a 3D model was examined but proved to be inconsistent. Consequently, a pragmatic solution was proposed (Section 2.2), which is based on the following points:

-No area other than the friction bands experiences flake opening/closure.

-Damage deactivation is introduced on the friction bands.

-The opening condition is defined by σ ≥ 0 in the considered direction.

-Flake closure is enforced at the very beginning of the unloading phase.

The model allows an accurate simulation of the mechanical behavior in the groove. In the friction bands, the asymptotic global behavior is also well described. However, a small error during the unloading phase, which is inherent to the damage deactivation procedure, still exists. This error gives rise to a small overestimation of the mechanical energy dissipated during a braking cycle.

This model must now be tested on a large number of brake discs with varied designs to be completely validated. Low cycle fatigue criteria, based on quantities such as mechanical dissipated energy will also be evaluated, in order to work out a complete brake disc fatigue design framework.

Fig. 1 :

 1 Fig. 1: Schematic drawing of induced damage

Fig. 2 :

 2 Fig. 2: Brake disc system

Fig. 3 :

 3 Fig. 3: Temperature field at maximum temperature time of the friction bands

  in the form of radial cracks on friction bands (Figure 4(b)).

Fig. 4 :

 4 Fig. 4: Cracked areas on a brake disc subjected to severe bench test

Fig. 5 :

 5 Fig. 5: Eigen stresses versus 1-component of the traction on the critical plane

Fig. 6 :

 6 Fig. 6: Block diagram of the basis change procedure in the UMAT routine

Fig. 7 :

 7 Fig. 7: Stress/strain paths on the D266 disc without damage

Fig. 8 :

 8 Fig. 8: Stress/strain path on the D330 friction band without closure condition

Fig. 9 :

 9 Fig. 9: Stress/strain path in the D330 friction band with damage deactivation at σ = 0

Fig. 10 :

 10 Fig. 10: Stress/strain path in the D330 friction band with damage deactivation upon early unloading

  are preferred.

	Matrix P defined above reads
			a 11 a 12 a 12	
	P =		a 21 a 22 a 23	
			a 31 a 32 a 33

  33 2a 33 a 31 2a 31 a 32 a 21 a 31 a 22 a 32 a 23 a 33 a 22 a 33 + a 23 a 32 a 21 a 33 + a 23 a 31 a 22 a 31 + a 21 a 32 a 31 a 11 a 32 a 12 a 33 a 13 a 12 a 33 + a 13 a 32 a 13 a 31 + a 11 a 33 a 11 a 32 + a 12 a 31 a 11 a 21 a 12 a 22 a 13 a 23 a 12 a 23 + a 13 a 22 a 13 a 21 + a 11 a 23 a 11 a 22 + a 12 a 21 21 a 31 2a 22 a 32 2a 23 a 33 a 22 a 33 + a 23 a 32 a 21 a 33 + a 23 a 31 a 22 a 31 + a 21 a 32 2a 31 a 11 2a 32 a 12 2a 33 a 13 a 12 a 33 + a 13 a 32 a 13 a 31 + a 11 a 33 a 11 a 32 + a 12 a 31 2a 11 a 21 2a 12 a 22 2a 13 a 23 a 12 a 23 + a 13 a 22 a 13 a 21 + a 11 a 23 a 11 a 22 + a 12 a 21

			11 a 2 21 a 2 31	a 2 12 a 2 22 a 2 32	a 2 13 a 2 23 a 2 33	2a 12 a 13 2a 22 a 23 2a 32 a  2a 13 a 11 2a 11 a 12 2a 23 a 21 2a 21 a 22       
									(66)
	and
	N =	       	a 2 11 a 2 21 a 2 31	a 2 12 a 2 22 a 2 32	a 2 13 a 2 23 a 2 33	a 12 a 13 a 22 a 23 a 32 a 33	a 13 a 11 a 23 a 21 a 33 a 31	a 11 a 12 a 21 a 22 a 31 a 32	       
									(67)
	Definitions (66) and (67) imply that N -1 = M T , which helps defining the tangent
	matrix L in the global basis
						L = M L D N T			(68)
	3.5 Stress/strain paths in the structure

2a

Acknowledgements This work was carried out within the Executive Doctorate in Industry framework, and was supported by PSA Peugeot Citroën.