M A Laribi

S Tamboura

Joseph Fitoussi

A Tcharkhtchi

H Ben Dali

R Tié Bi

Fast fatigue life prediction of short fiber reinforced composites using a new hybrid damage approach: Application to SMC

Keywords:

come

Fast fatigue life prediction of short fiber reinforced composites using a new hybrid damage approach:

Application to SMC

Introduction

For an efficient design of short fiber reinforced composites SFRC automotive components, several modes of failure have to be taken into account by the manufacturers. Particularly, fatigue loading may be considered as one of the most critical failure mode which strongly affects the component lifetime. In fact, poor knowledge and control of this mode of failure lead to specifications provided by the manufacturers which often need to be clarified. Moreover, the existent experimental databases, which should allow at least a pre-design, remain relatively poor especially for new material formulations. On the other hand, the development times become increasingly short and must also integrate the material formulation optimization. Several iterations are often necessary in order to identify the correct couple formulationprocess suitable for each new application. Thus, it seems essential to build a rapid methodology of pre-dimensioning of composite structures allowing, within a short time, to identify the appropriate material formulation for each new automotive SFRC component. This methodology must therefore be pragmatic for industrial purpose. In the case of components submitted to fatigue loading, this is a critical issue because of the too long and too costly experimental campaigns needed to established SN curves.

Moreover, because of their heterogeneous architecture, anisotropic behavior and damage development, short fiber reinforcement composites SFRC structure design is very complex compared to homogeneous materials such as metals. Indeed, the mechanical properties of these materials were found to be greatly affected by their composition [START_REF] Degrieck | Fatigue damage modelling of fibre-reinforced composite materials[END_REF]: matrix type (thermoplastic or thermoset), volume fraction of reinforcing phase and their geometry and spatial distribution [START_REF] Simoes | Influence of concrete strength and steel fibre geometry on the fibre/matrix interface[END_REF][START_REF] Rolland | In situ X-ray tomography investigation on damage mechanisms in short glass fibre reinforced thermoplastics: effects of fibre orientation and relative humidity[END_REF][START_REF] Mortazavian | Effects of fiber orientation and anisotropy on tensile strength and elastic modulus of short fiber reinforced polymer composites[END_REF][START_REF] Roundi | Experimental and numerical investigation of the effects of stacking sequence and stress ratio on fatigue damage of glass/epoxy composites[END_REF]. The diversity of composite material architectures and their complexity has for each application its own specificity. This wide variety of microstructures induces a significant dispersion of the mechanical properties, in particular under fatigue loading where the prediction of thresholds and kinetics of damage, their effect on the residual mechanical behavior and on the lifetime of the materials can present a very important challenge. These scattered results on composite S-N curves have been demonstrated by Barnard and al [6]. The authors admit the change of the damage mechanism as a result of microstructure variation and damage during the different stages of the fatigue process. They have also shown that there is a direct relationship between the strength statistical T variation and the scatter observed in the composite S-N curves. One can expect that in the case of SFRC materials, this effect should be largely exacerbated because of the high variability of the microstructure generated during their molding process [START_REF] Shirinbayan | High strain rate visco-damageable behavior of Advanced Sheet Molding Compound (A-SMC) under tension[END_REF]. In fact, even for a single formulation, numerous different microstructures are present in a real SFRC automotive component. Consequently, a predictive modeling must imperatively integrate the microstructure as an input data.

Therefore, considering the complexity of the local mechanisms occurring during loading, the chosen approach must be based on a solid knowledge of these different phenomena. In the case of short fiber reinforced composites such as Sheet Molding Compound, it has been widely demonstrated that local damage must be a central focus of the design methodology [START_REF] Degrieck | Fatigue damage modelling of fibre-reinforced composite materials[END_REF][START_REF] Arif | Multiscale fatigue damage characterization in short glass fiber reinforced polyamide-66[END_REF]. Unlike thermoplastic based SFRC, non-linear behavior due to matrix visco-plasticity does not have to be taken into account at ambient temperature and low frequency.

However, in the case of fatigue behavior, environmental and loading conditions effects (temperature [START_REF] Bellenger | Thermal and mechanical fatigue of a PA66/ glass fibers composite material[END_REF], stress ratio [START_REF] Khan | Effect of stress ratio or mean stress on fatigue delamination growth in composites: critical review[END_REF][START_REF] Andersons | Empirical model for stress ratio effect on fatigue delamination growth rate in composite laminates[END_REF], cycling frequency [START_REF] Kotik | Frequency effect in short-beam shear fatigue of a glass fiber reinforced polyester composite[END_REF][START_REF] Eftekhari | On the strengthening effect of increasing cycling frequency on fatigue behavior of some polymers and their composites[END_REF]) are considered to be of first importance. Thus, the chosen methodology should be able to be adapted to predict the lifetime of these materials under these various complex conditions. On the other hand, fatigue damage and failure behavior of SFRC materials appears to be a very complicated and various process. Mainly, the existing predictive fatigue models for SFRC are generally classified into three categories [START_REF] Degrieck | Fatigue damage modelling of fibre-reinforced composite materials[END_REF][START_REF] Passipoularidis | Fatigue evaluation algorithms: review[END_REF]; empirical, based on S-N curves; phenomenological, based on residual strength or on residual stiffness evolutions and finally micromechanical, based on microscopic damage mechanisms description.

Great efforts have been expended to develop new methods allowing reducing substantially cost and time dedicated to characterization of fatigue behavior [START_REF] Palumbo | A new rapid thermographic method to assess the fatigue limit in GFRP composites[END_REF][START_REF] Rajaneesh | Long-term life prediction of woven CFRP laminates under three point flexural fatigue[END_REF][START_REF] Jegou | Fast prediction of the Wöhler curve from heat build-up measurements on Short Fiber Reinforced Plastic[END_REF][START_REF] Marco | Prediction of fatigue properties of natural rubber based on the descriptions of the cracks population and of the dissipated energy[END_REF][START_REF] Abello | Fast prediction of the fatigue behavior of short fiber reinforced thermoplastics from heat build-up measurements[END_REF][START_REF] Masquelier | Thermal measurements on elastomeric materials: from the characterization of the dissipation gradients to the prediction of the fatigue properties[END_REF][START_REF] Marco | Fast prediction of the Wohler curve from thermal measurements for a wide range of NR and SBR[END_REF]. L. Jegou and al. proposed a new approach [START_REF] Jegou | Fast prediction of the Wöhler curve from heat build-up measurements on Short Fiber Reinforced Plastic[END_REF] based on the temperature measurements. An energetic fatigue criteria, identified from measurement results, is obtained from only one single microstructure in order to predict the fatigue curve from the heat build-up curve. The duration of the whole procedure is less than two days. Yann Marco and al [START_REF] Marco | Prediction of fatigue properties of natural rubber based on the descriptions of the cracks population and of the dissipated energy[END_REF] proposed in another recent paper a similar approach by relating the fatigue lifetime to the crack population evaluated from interrupted fatigue tests and the energy dissipation for a Natural Rubber compound filled with carbon black. This approach [START_REF] Jegou | Fast prediction of the Wöhler curve from heat build-up measurements on Short Fiber Reinforced Plastic[END_REF][START_REF] Marco | Prediction of fatigue properties of natural rubber based on the descriptions of the cracks population and of the dissipated energy[END_REF][START_REF] Abello | Fast prediction of the fatigue behavior of short fiber reinforced thermoplastics from heat build-up measurements[END_REF][START_REF] Masquelier | Thermal measurements on elastomeric materials: from the characterization of the dissipation gradients to the prediction of the fatigue properties[END_REF][START_REF] Marco | Fast prediction of the Wohler curve from thermal measurements for a wide range of NR and SBR[END_REF] proposes an efficient Wohler curves prediction based only on thermal measurements and an energetic criterion and has been applied on a wide variety of thermoplastic polymers and rubbers [START_REF] Marco | Fast prediction of the Wohler curve from thermal measurements for a wide range of NR and SBR[END_REF]. However, the authors do not specify if this approach should be effective in the case of thermoset polymer or their composites presenting a lower thermal effect especially with low operating frequency and amplitude [START_REF] Shirinbayan | Coupled effect of loading frequency and amplitude on the fatigue behavior of advanced sheet molding compound (A-SMC)[END_REF].

Jain and al [START_REF] Jain | The Master SN curve approach -A hybrid multi-scale fatigue simulation of short fiber reinforced composites[END_REF][START_REF] Jain | A feasibility study of the Master SN curve approach for short fiber reinforced composites[END_REF] have proposed an approach based on a micromechanical model relating the macroscopic damage state in quasi-static loading to that reached as a certain point on a givenS-N curve. The authors propose a methodology allowing prediction of the S-N curves for other composites microstructures from the knowledge of the so called "master curve" which corresponds to the S-N curve established for one specific microstructure taken as a reference. However, the proposed methodology is based on the experimental statement that for similar number of cycles to failure, the difference in the loss of stiffness curves for different microstructures is not significant and can be statistically treated as being the same. Experimental evidence shows that, in the case of SMC materials, the loss of stiffness kinetic and its critical value just before failure is strongly dependent of the microstructure.

In this paper, a hybrid modeling approach is proposed. In order to take into account microstructure and damage effect, this model is based on a Mori and Tanaka predictive approach integrating a damage criterion at the microstructure scale. This model is identified in the case of a simple tensile loading for one single microstructure and is able to predict the monotonic mechanical response for each other microstructure. Then, using the output of this micromechanical model, a phenomenological formulation allows a rapid evaluation of the SN curve for different other microstructures. In order to drive the model, limited experimental data from laboratory tests are needed.

In the next section, the different steps of the proposed methodology are presented. Then, the micromechanical damage model is described until the identification of an equation of state applicable under monotonic tensile loading. The model is applied to a SMC composite material which is presented together with the experimental methods and results in the third section. The last section is dedicated to write an equation of state under fatigue loading to be solved in order to establish SN curves with minimum experimental data. The methodology is validated for three kinds of microstructure: a Randomly Oriented SMC (RO) and a Highly Oriented SMC under two fatigue loading directions (HO-0°and HO-90°).

Presentation of the new hybrid methodology

In this study, we are looking for a pragmatic methodology of rapid design of SFRC materials submitted to fatigue loading. This approach is suitable for thermoset based composite materials such as SMC composites for which non-linear behavior is essentially attributed to damage phenomenon. Therefore, this is a stiffness based model based on a micromechanical approach in which a damage criterion in introduced at the local scale. This model is presented in the next section. Further development of our approach will be proposed in the case of thermoplastic based composites for which matrix non-linearity should be taken into account.

The key idea of this approach is to introduce the output of an accurate micromechanical damage model into a phenomenological approach. Such a hybrid model combines the advantages of both micromechanical and phenomenological approaches: In one hand, the proposed methodology is pragmatic and well adapted to industrial design. On the other hand, our model keeps the reliability of the micromechanical approach through the description of local damage evolution. As it has been shown in previous studies [START_REF] Tamboura | Evaluation de la tenue en fatigue du composite SMC R42[END_REF][START_REF] Jendli | Micromechanical analysis of strain rate effect on damage evolution in sheet molding compound composites[END_REF], the non-linear behavior of standard SMCs submitted to quasi-static, dynamic or fatigue loading is mainly governed by the same phenomenology: fibre-matrix interface failure development. Under tensile or fatigue loading, accumulation and propagation of local damage lead to a progressive loss of stiffness which can be considered as a good indicator of damage state. Therefore, independently of the loading type, any local damage state can be directly related to the corresponding degraded macroscopic properties. In other words, the relation between local damage and macroscopic loss of stiffness should be understood as an intrinsic relationship. Consequently, it is obvious to consider that fatigue damage and tensile damage can be described by the same type of relationship.

Previous experimental studies [START_REF] Tamboura | Evaluation de la tenue en fatigue du composite SMC R42[END_REF][START_REF] Jendli | Micromechanical analysis of strain rate effect on damage evolution in sheet molding compound composites[END_REF][START_REF] Jendli | Multi-scales modelling of dynamic behaviour for discontinuous fibre SMC composites[END_REF][START_REF] Fitoussi | Determination of a tridimentional failure criterion at the fiber-matrix interface of an organic matrix/discontinuous reinforcement composite[END_REF] showed that final failure appears when the micro-cracks density reaches a certain critical value. This critical value can be represented as a local damage rate, d dc where d is related to the current amount of fibre-matrix interface damage and d c is a critical value. In SMC composite materials, one can assume that the pseudo-delamination phenomenon which always appears just before failure [START_REF] Shirinbayan | High strain rate visco-damageable behavior of Advanced Sheet Molding Compound (A-SMC) under tension[END_REF][START_REF] Shirinbayan | Mechanical characterization of a Low Density Sheet Molding Compound (LD-SMC): Multi-scale damage analysis and strain rate effect[END_REF], takes place when the local damage rate is equal to 1. Therefore, one can consider the local failure criterion as:

= d d 1 c (1)
In the proposed fatigue life prediction model, this criterion is supposed to be an intrinsic characteristic of the material independently of the applied loading type: quasi-static or fatigue. Physically, it means that failure always occurs for the same critical amount of local damage independently of the loading history when the non-damage volume of the material becomes too small to withstand the applied load. Consequently, when the critical local damage rate is reached, coalescence of micro-cracks rapidly leads to final failure. Therefore, one can consider that the failure criterion described in equation ( 1) can be identified using quasi-static response and then be used for the prediction of fatigue degradation and failure. It is then obvious to consider that fatigue life prediction model for elastic damageable composites should be based on a quasi-static loading micromechanical approach taking into account microstructure and damage phenomenon at the local scales. The same idea is used in recent publications and applied to thermoset matrix composites [START_REF] Jain | The Master SN curve approach -A hybrid multi-scale fatigue simulation of short fiber reinforced composites[END_REF][START_REF] Jain | A feasibility study of the Master SN curve approach for short fiber reinforced composites[END_REF].

Continuous efforts have been dedicated in the last two decades to the experimental analysis and the multi-scale modeling of SMC composites submitted to quasi-static, high strain rate and fatigue behavior [START_REF] Jain | The Master SN curve approach -A hybrid multi-scale fatigue simulation of short fiber reinforced composites[END_REF][START_REF] Tamboura | Evaluation de la tenue en fatigue du composite SMC R42[END_REF][START_REF] Krause | Physically based micromechanical approach to model damage initiation and evolution of fiber reinforced polymers under fatigue loading conditions[END_REF][START_REF] Oskay | Numerical Modelling of Failure in Advanced Composite Materials[END_REF]. Thus, the present approach is based on the use of an existing micromechanical model predicting quasi-static and dynamic behavior of SMC composite [START_REF] Jendli | Multi-scales modelling of dynamic behaviour for discontinuous fibre SMC composites[END_REF]. In this model, a fibre-matrix interface failure is introduced in the Mori and Tanaka approach [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF]. Micromechanical parameters describing fibre-matrix interface failure development and the critical local damage value (d c ) can be identified by inverse engineering using tensile and loading-unloading tensile test results (see section 5).

Therefore, the local damage rate evolution can be determined and related to the macroscopic loss of stiffness in the case of a simple tensile loading. The relation between the local damage rate and the macroscopic loss of stiffness derived from the micromechanical damage modeling of this simple tensile test can be considered as an equation of state which can be then adapted and used as an intrinsic relationship for fatigue life prediction. The introduction of this equation of state into a stiffness based fatigue phenomenological formulation leads to the expression of the evolution of the local damage rate during fatigue loading. Finally, the failure criterion expressed in equation ( 1) allows predicting S-N curves of SMC materials for different microstructures. The whole methodology is represented schematically in Fig. 1.

One can summarize the rapid evaluation of the Wöhler curve by the four following stages: 1) Performing tensile and loading-unloading tests in order to determine the damage threshold σ s and the loss of stiffness kinetics under monotonic tensile loading represented by the parameter a (see section 5). 2) Performing tension-tension fatigue tests without the necessity to reach the failure in order to identify the kinetics of reduction in stiffness in fatigue loading represented by the parameter B (see section 6).

3) Identification of the micromechanical model for monotonic loading in order to determine the tensile equation of state, 4) Introduction of the identified loss of stiffness kinetics into the equation of state and resolution in order to plot the Wohler curves (cf. equation ( 14)).

This procedure can be used for each new microstructure in order to obtain a rapid evaluation of the corresponding Whöler curve.

Material and experimental methods

Our methodology is applied to an industrial SMC material. Hereafter are described material and experimental methods.

Material

The material in study is polyester based Sheet Molding Compound composite (SMC). Glass fibers and CaCO3 particles weight contents are of 28% and 37% respectively are introduced in the formulation.

Glass fibers are presented in the form of bundles, with constant length L = 25 mm. Each bundle contains approximately 250 glass fibers of 15 μm diameter. Bundles are originally randomly oriented in the plane of the non-reticulated SMC flanks to be processed by thermocompression.

SMC plates are obtained by thermo-compression under a pressure of 60 Bars and an average temperature of 165 °C of the mold. Reticulation time is less than two minutes. Two molding configurations have been used. In the first one, non-reticulated SMC flanks have been placed in the center of the rectangular mold (120 × 250 mm 2 ) and recover more than 80% of its surface. Consequently, limited flowing during compression molding leads to a quit random orientation of the fibers (RO material). On the other hand, the same amount of SMC flanks has been placed in one edge of the mold in such a way that only 50% of the mold was recovered. This leads to a highly oriented distribution of fibers in the flow direction (HO material).

Mechanical characterization methods

Overall mechanical characterizations of the SMC (quasi-static and fatigue) behavior have been performed on a MTS 830 hydraulic machine.

Specimen geometry

Tensile and fatigue specimens were cut from RO and HO plates. Fig. 2 shows the used specimen geometry.

In the case of HO microstructure, two specimen orientations were chosen. HO-0°and HO-90°denote specimens cut in the flow direction and perpendicularly to it respectively.

Test configurations

Three test configurations have been used:

-Tensile tests until failure (in the two main directions for the oriented material). Thus, average final failure stresses are provided for the three test configurations: RO, HO-0°and HO-90°. -Loading-unloading tensile tests with progressive increase of the maximum load. The main aim of these tests was to plot the relative loss of stiffness evolution ( E E0

) versus the applied stress. -Tension-tension stress controlled fatigue tests until failure at several applied maximum stress. The chosen stress-ratio and frequency were R = 0.1 and f = 10 Hz. All fatigue tests were preceded by a tensile loading-unloading-elastic reloading cycle allowing a better identification of the non-damaged material Young's modulus, E 0 and E 1 , corresponding to the residual stiffness after the first cycle. The relative loss of stiffness evolution during fatigue is determined and Whöler curves are plotted for each microstructure configuration. Note that the systems cannot reach the requested cycling immediately at the first applied cycle. For an applied frequency of 10 Hz, it takes generally N s cycles to reach the goal (around 20 cycles in our case). So, it is considered that the required experimental conditions are properly applied only from the N s cycle. We verify for all performed fatigue tests that no evolution of the relative loss of stiffness is observed until cycle N s . Thus, one can assume that = The micromechanical multi-scale modeling proposed here is derived from the Jendli & al. work [START_REF] Tamboura | Evaluation de la tenue en fatigue du composite SMC R42[END_REF]. This predictive modeling involves the introduction of local damage into standard homogenization technique. The latter allows predicting the overall behavior of heterogeneous materials from their constituent's properties. The damage model considered in this study is based on the theory proposed by Mori and Tanaka [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] and adapted by Benveniste [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory[END_REF]. Mori and Tanaka's theory gives access to the whole elasticity tensor of the composite as a function of microstructure parameters such as the fiber content, distribution of orientation and length of fibers. Composite stiffness tensor is given by the expression:

= + + - -- C C I f Q I f S I Q [ ( ( ) ) ] Comp m 1 1 (2)
where C m and C i are the stiffness tensor of the matrix and the re- inforcement respectively, f represents the reinforcement volume frac- tion and Q is a pseudo-localization tensor defined for each family of reinforcement "i" as:

= - -- - - Q ((C C )(S I) C ) (C C ) i m i i i 1 i m (3)
This expression integrates the Eshelby tensor S [START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF]. Moreover, for a given macroscopic applied stress ∑, average stress and strain fields in each phase of the material (matrix or fiber) can be determined. Several criteria of local damage can then be introduced and put in competition at the local scale. Because fiber-matrix interface debonding is the predominant damage mechanism occurring in SMC composites, local stresses on the interface must be calculated. For a given macroscopic applied stress ∑, Mori and Tanaka theory gives ac- cess to the average stress tensor for each family of fibers oriented in the "i" direction as:

∑ = + - + - - σ C ( I (S I)Q )(I f(S I)Q )C i m i i m 1 (4)
Using the continuity condition of the normal stress at the interface, the average stress in the fiber given by equation ( 4) is used to calculate the normal and shear interfacial stresses, σ n and τ respectively, at each point located at the fiber-matrix interface:

= → → → = → = → - σ T•n, with T σ • n and τ T (σ ) n i i n 2 (5)
where → n is the normal vector at the considered point of the interface see Fig. 3 (a).

Fitoussi & al [START_REF] Fitoussi | A statistical micromechanical model of anisotropic damage for S.M.C. composites[END_REF] propose the introduction of a statistical local interface failure criterion in order to predict threshold and kinetic of damage at the local scale. Jendli & al [START_REF] Jendli | Multi-scales modelling of dynamic behaviour for discontinuous fibre SMC composites[END_REF] generalized this approach for dynamic loading. At each calculation step, a local failure probability is calculated at each point located on the fiber-matrix interface and defined by two angles, θ and ϕ. The first one is related to the in-plane orientation of the considered fiber and ϕ defines the position of the interfacial point around the fiber equator. The probability function is written according to Weibull law [START_REF] Fitoussi | Determination of a tridimentional failure criterion at the fiber-matrix interface of an organic matrix/discontinuous reinforcement composite[END_REF]:

⎜ ⎟ ⎜ ⎟ ⎜ ⎟ = - ⎛ ⎝ ⎜ -⎛ ⎝ ⎛ ⎝ ⎞ ⎠ + ⎛ ⎝ ⎞ ⎠ ⎞ ⎠ ⎞ ⎠ ⎟ P (θ, φ) 1 exp σ σ τ τ r 0 2 0 2 m (6)
where τ and σ are the shear and normal stresses derived from equation (5). σ 0 and τ 0 are interfacial strengths under normal and shear stress at the interface respectively. "m" is a statistical parameter directly function of the local distribution of microstructure.

Local damage is introduced in the composite material through a density of zero stiffness penny shape heterogeneities describing the presence of the micro cracks coming from the interfacial damage. A damaged state should be described by micro-cracks and fibers. Two kinds active fibers participating in the composite reinforcement should be distinguished: the non-debonded fibers and the part of the debonded fibers which remain active (see Fig. 3 (b)).

Therefore, at each calculation step, n, local interface failure probabilities are calculated and gives access to the amounts of remaining active fiber, f n act , including the undamaged fibers volume fraction, f n ND , the part of the debonded fiber which continue participating in the composite reinforcement, and the micro-cracks volume fraction to be introduced. Therefore, the local damage state is described at each step, n, by the following equations:

= - - f ( 1 P )* f n ND r n n 1 ND ∑ = + = - f f k P. f n act n ND i 1 n r i 1 ND = + - - f f h . P . f n mc n 1 mc r n n 1 ND ( 7
)
where k is a reduction coefficient applied to partially debonded fibers. This parameter has been evaluated by finite element calculations performed on a representative cell containing a partially debonded fiber and is equal to 0,5. h is the ratio between the volume of the introduced penny shape (representing interfacial micro-crack) and the fiber evaluated by geometric considerations. Therefore, the homogenization of the composite material containing remaining active fibre, matrix and micro-cracks can be achieved according to a two successive stages as described in Jendli & al. paper [START_REF] Jendli | Multi-scales modelling of dynamic behaviour for discontinuous fibre SMC composites[END_REF].

Experimental results and micromechanical model identification

Loss of stiffness during tensile loading

Tensile loss of stiffness is measured by the loading-unloading procedure described in section 3.2.2. The micromechanical model described above allows predicting the evolution of the properties of the equivalent homogeneous material. As a result, one can easily plot the evolution of the progressive reduction of stiffness during a tensile loading. The key issue of the micromechanical model lies in the identification of the local fiber-matrix interface failure criterion parameters. The latter is identified by reverse engineering on the basis of the reduction in stiffness measured in the testing microstructure. The identification procedure is described in Ref. [START_REF] Jendli | Multi-scales modelling of dynamic behaviour for discontinuous fibre SMC composites[END_REF]. Comparison between model and experimental data is shown in Fig. 4.

For the sake of simplification, the average evolution of the relative loss of stiffness should be considered. In a first approach, a linear function can be chosen to represent the average evolution of the relative loss of stiffness. This linear function is defined in the range of stress values corresponding to the evaluated acceptable fatigue applied stresses in real structures.

Thus, relative loss of stiffness during a tensile loading can be expressed by:

⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = + - E E 1 a(σ σ ). 0 S ( 8
)
where σ S is the average damage threshold stress, a is a kinetic para- meter and σ is the applied stress. Therefore, when considering another microstructure.

Determination of an accurate local damage indicator

The hybrid approach proposed in this paper relies on the definition of a local damage indicator and its critical value corresponding to a critical damage state leading to failure. Thus, the definition of the local damage indicator, d, and its associated critical value, d c , should be related to the local parameters describing the local damage state; namely, the micro-crack density f n mc , and the remaining active fibers, f n act which are the output of the micromechanical modeling.

The chosen indicator is:

= * d f f mc act (9)
Final failure stress is needed to identify, through a reverse engineering method, the critical value of the local damage rate, d c . Indeed, the micromechanical model allows plotting the evolution of the local damage amount, d, until the measured failure stress for the considered microstructure. The corresponding value to failure is equal to d c . For the sake of confidentiality normalized average failure stress values is defined by the ratio between the failure stress and the average failure stress of RO samples, Thus, the normalized average failure stress of RO, HO-90°and HO-0°samples are 1( ± 0,074), 0,83( ± 0,076) and 1,69( ± 0,15). It is obvious to note that for the three studied microstructure (RO, HO-0°and HO-90°) identified values of d c are of the same order of magnitude. This indicates that the critical local configuration leading to the final failure is independent of the microstructure.

Equation of state under tensile loading

The micromechanical model allows plotting the evolution of the local damage rate, d dc as a function of the stiffness reduction, E E0 (see Fig. 5). As we mentioned above, this relation is considered as an equation of state relating microscopic damage state to macroscopic properties degradation. Therefore, this relationship is supposed to be accurate even for other loading scheme such as fatigue.

It is easy to describe the equation of state as a polynomial function of the second degree:

⎜ ⎟ ⎜ ⎟ = ⎡ ⎣ ⎢ ⎛ ⎝ ⎞ ⎠ ⎤ ⎦ ⎥ + ⎡ ⎣ ⎢ ⎛ ⎝ ⎞ ⎠ ⎤ ⎦ ⎥ + d d α E E β E E γ c i 0 2 i 0 i (10)
α, β and γ are micromechanical parameters, "i" index indicates the considered microstructure, (RO, HO-0°and H0-90°).

Hybrid model and fatigue life prediction

In order to emphasize the degradation kinetic, it is interesting to plot the evolution of the relative loss of stiffness E E1 , where E is the residual stiffness for the current cycle and E 1 is the residual stiffness after the first cycle. Comparison between relative loss of stiffness E E1 , plotted for different values of applied fatigue stress highlights an average damage kinetic which can be described by the following power function:

⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = E E AN 1 N B (11)
One can notice that, in the case of a SMC composite material, the variation of the slope around the average kinetic does not depend on the applied stress level. This experimental evidence has been also confirmed in the case of other microstructure configurations such as HO-0°and HO-90°. The scatter observed on the damage kinetic and on the fatigue life time is essentially due to the dispersion of microstructure which has been demonstrated to be very variable in the case of SMC composites. Like in tension loading, this topic is not addressed and will be the subject of a further paper.

In Fig. 6, all the curves are plotted from the N s cycle which should be considered equivalent to the first cycle as mentioned above.

Besides, relative loss of stiffness during fatigue loading can be decomposed into two terms:

⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = ⎛ ⎝ ⎞ ⎠ ⎛ ⎝ ⎞ ⎠ E E E E * E E 0 1 1 0
The first term, is related to the stiffness reduction consecutive to the first cycle and is determined as described in section 3.2.2. Therefore, relative loss of stiffness during fatigue can be written as:

⎜ ⎟ ⎛ ⎝ ⎞ ⎠ = + - ⎢ ⎣ ⎢ ⎥ ⎦ ⎥ E E (1 a(σ σ ))* N N 0 N S s B ( 12
)
Arising from the introduction of equation ( 12) into equation ( 10), an analytical expression of the evolution of the local damage rate The resolution of equation ( 13) leads to the expression of the number of cycles to failure as a function of the applied stress:

= ⎡ ⎣ ⎤ ⎦ = + - -- - - N N* with G r s G (1 a (σ σ )) 1 B i i β β 4α (γ 1) 2α i i i S i i 2 i i i (14)
This expression allows plotting SN curve for any microstructure configurations. Damage is indirectly taken into account though the micromechanical parameters, α i , β i and γ i . The initial damage is conditioned by two parameters; a i and σ i S and fatigue kinetics of damage depends on the parameter B i . The five first parameters are directly derived from the micromechanical model as a function of the considered microstructure. B i can be obtained by measuring the slope of the loss of stiffness under fatigue loading. Practically, one can perform a fatigue test at a maximal stress equal to half of the ultimate stress. While the experimental data expected from this test is the B i parameter (corresponding to the loss of stiffness slope), it is not needed to complete the test until final failure. Note that this test should be repeated several times for scattering materials such as SMC composites.

The proposed procedure was applied for three different microstructure configurations: RO, HO-0°and HO-90°. A micromechanical tensile modeling was performed for each configuration in order to identify α , β , γ a , σ . i i i, i i S B i was identified from a fatigue test as described above. For the sake of validation of our methodology, a fatigue campaign was completed for the three considered microstructures. Finally, comparison of simulated Wohler curves to experimental results for RO, HO-0°and HO-90°samples is illustrated in Fig. 7. One can note a very good agreement which demonstrates the efficiency of the proposed methodology.

Highlights and actual limitations of the approach

Because of the character of novelty of the proposed approach, it is necessary to discuss about its potential and actual limitations. The main assumption is that there is a link between monotonic and fatigue mechanical behavior. Some authors also recently highlighted a relation between fatigue limit and quasi-static damage threshold [START_REF] Lomov | Fatigue limit: is there a link to quasi-static damage? Conference paper[END_REF]. Indeed, experimental evidence shows that the same damage and failure mechanisms are observed for the two types of loading. Therefore, the relationship between the crack density and the loss in stiffness can be determined in the case of a simple monotonic tensile test and be used in the case of tensile fatigue loading.

It has been demonstrated that this approach is efficiency in the case of constant tensile loading. On the other hand, real structures are submitted to variable amplitude loading. However, under tensile loading, it has been shown that the macroscopic fatigue damage rate (characterized by B) can be considered as independent of the amplitude in the case of standard SMC composites. Therefore, it can be assumed that each microscopic damage state can be represented by an overall local damage value (here d dc ) independently of the tensile history which has led to this specific damage state. Thus, the increment of local damage ( Δd dc ) is calculated at the local scale as a function of the actual damaged microstructure for the actual macroscopic stress (σ in equation ( 13)) independently of the tensile loading history. A further paper will illustrate the efficiency of this approach even for variable amplitude fatigue loading. Several variable amplitude cases will be illustrated for validation of the approach in the case of tensile fatigue (to be submitted to Composites Part B).

Many other fatigue parameters influences should be evaluated in order to improve the efficiency of our approach. For example, the effect of compression and shear stresses and more generally the effect of complex loading and that of the mean stress should be studied. Actually, the model only takes into account the effect of the maximum applied stress. In a first approximation, this assumption can be considered valuable since the damage threshold under compression is 50% higher than tensile damage threshold. In the case of shear stress, the maximum absolute value must be considered. For each of these cases, it should be also verified that local damage scenario are similar in monotonic and fatigue loading. Moreover, even if the micromechanical approach allows the identification of a , σ , α , β , γ i i S i i i , the macroscopic damage rate, B i , must be evaluated experimentally together with its sensitivity to amplitude for any fatigue loading scheme.

One should also consider the limitation of the approach in the case of thermoplastic matrix. Indeed, the model should be improved in order to take into account the effect of the local plastic deformation accumulated under fatigue loading.

Conclusion and perspectives

In this paper, we propose a new hybrid approach allowing fast Wöhler curves prediction for SMC composite materials as a function of the considered microstructure. As fatigue tests require long and costly experimental campaigns, it is of first importance to speed the fatigue life evaluation which should be determined several times in the design phases which involve several compounding and processing steps.

On one hand, this approach can be qualified as a "stiffness based fatigue life prediction model". On the other hand, the model is based on a multi-scale description of damage. A micromechanical model is used and integrates, at the local scale, the main damage phenomenon observed in the case of SMC submitted to monotonic or fatigue loading: fiber-matrix interface failure. This pragmatic approach supposes that there is an equation of state which remains valid for both monotonic and fatigue loading. In other words, we consider that for a given residual stiffness state, the global quantity of local damage cumulated in the material is equivalent for tensile or fatigue loading. In light of this hypothesis, the equation of state can be established using a micromechanical modeling whose parameters can be identified in the case of a simple tensile test and be used in the case of fatigue loading. Moreover, final failure is assumed to always occur when a critical amount of local damage is reached. The latter is also identified by the micromechanical model by calculating the amount of local damage reached at the measured tensile failure. It has been shown that, in the case of tensile loading, a direct relation between the local damage rate and the number of applied cycles can be derived. The resolution of this equation for the identified critical amount of local damage allows plotting the SN curve for the considered microstructure.

The predictions of our hybrid approach are validated with experimental results in the case of constant tensile fatigue loading. Indeed, the SN curves are in very good agreement for three different microstructure configurations: RO, HO = 0°and HO = 90°.

The originality of this approach lies in using local damage state as variable in the fatigue behavior description. The use of a micromechanical approach which integrates microstructure as an input data is a major asset of this fast prediction method. Moreover, the model was implemented in industrial software and used as a damage approach for composite system's reliability in the case of the new Peugeot 3008 tailgate reinforced by SMC composite material [START_REF] Tiebi | Damage approach for composite system's reliability: application on SMC component in automotive industry[END_REF]. Application of this approach under a thermo mechanical fatigue loading required by Peugeot has been also validated [START_REF] Tiebi | Damage approach for composite system's reliability: application on SMC component in automotive industry[END_REF]. However, it should be noticed that the presented approach is relevant only for elastic damageable composite materials such as SMC. Ongoing work concerns the development of this approach for thermoplastic matrix composites for which matrix plasticity should be coupled to damage at the local scale. An evaluation of the approach should be performed in the case of complex loading (at least two dimensional loading for industrial applications) and variable fatigue loading. In the same way, in order to generalize the approach for other materials, one should also consider other local damage mechanisms such as fiber breakage or matrix micro-cracking which can be placed in competition with fiber-matrix interface failure.

Fig. 1 .

 1 Fig. 1. Representation of the proposed methodology.

Fig. 2 .

 2 Fig. 2. SMC specimen geometry used for tensile and fatigue tests.

Fig. 3 .

 3 Fig. 3. (a) Projection of the stress tensor (normal and shear interface stresses, σ and τ respectively). (b) Micromechanical damage description. (From Jendli & al. [27]).

Fig. 4 .

 4 Fig. 4. Evolution of the relative loss of stiffness during tensile loading.

E

  E1 defines the fatigue loss of stiffness rate described by the parameter B whereas E E 1 0

  the evolution of the local damage state d dc , for each considered number of cycle, N, applied at the considered maximum stress, σ. As defined above, macroscopic failure occurs when , where N r is the number of cycle to failure.

Fig. 5 .Fig. 6 . 7 .

 567 Fig. 5. Evolution of the local damage rate in function of the loss in stiffness.


Rodrigo Méndez Rojano
email: rodrigo.mendez-rojano@umontpellier.fr

Simon Mendez

Franck Nicoud

Rodrigo Méndez

Introducing the pro

Keywords: Thrombosis, Medical devices, Contact System, Hemodynamics

Introduction

The use of blood-coated devices is a common medical practice to treat several cardiovascular conditions. However, the poor thrombotic performance of some types of medical devices remains an issue in clinical practice [START_REF] Chan | Catheter thrombosis and percutaneous coronary intervention:fundamental perspectives on blood, artificial surfaces and antithrombotic drugs[END_REF]; [START_REF] Mehra | The vexing problem of thrombosis in long-term mechanical circulatory support[END_REF]; [START_REF] Hastings | Sources of circuit thrombosis in pediatric extracorporeal membrane oxygenation[END_REF]). For instance, thrombus formation may lead to device malfunction, thrombo-embolism or increased risk of a hemorrhagic stroke due to the anti-coagulant therapy. Thrombosis is a complex phenomenon that involves several processes such as coagulation reactions, platelet activation/aggregation, blood flow dynamics, among others. Furthermore, these mechanisms interact with each other at different time and spatial scales as explained in [START_REF] Furie | Mechanisms of thrombus formation[END_REF]. When it comes to device-related thrombosis, the artificial wall of the device interacts with the thrombotic processes as pointed out by [START_REF] Gorbet | Biomaterial-associated thrombosis: roles of coagulation factors, complement, platelets and leukocytes[END_REF]. Wall protein adsorption (Vroman effect) drives contact activation that initiates the coagulation reactions. At the same time, platelets adhesion and activation take place [START_REF] Jaffer | Medical device-induced thrombosis: what causes it and how can we prevent it?[END_REF]), constituting the platelet-fibrin clot that will form the stable thrombus. In recent studies, computational fluid dynamics has been used to predict the flow in medical devices, whith the long-term objective of improving their thrombotic performance [START_REF] Alemu | Design optimization of a mechanical heart valve for reducing valve thrombogenicity -a case study with ats valve[END_REF]; [START_REF] Piatti | Hemodynamic and thrombogenic analysis of a trileaflet polymeric valve using a fluid-structure interaction approach[END_REF]). However, most of the thrombus formation models are dedicated to platelet activity or in-vivo situations in which no medical device is present. For instance, numerical studies that consider platelet activity are notably focused on the characterization of mechanical activation of platelets [START_REF] Hellums | whitaker lecture: Biorheology in thrombosis research[END_REF]; [START_REF] Kroll | Platelets and shear stress[END_REF]; [START_REF] Jetsy | Platelet activation in a circulating flow loop: combined effects of shear stress and exposure time[END_REF]; [START_REF] Nobili | Platelet activation due to hemodynamic shear stresses: Damage accumulation model and comparison to in vitro measurement[END_REF]; [START_REF] Shadden | Potential fluid mechanic pathways of platelet activation[END_REF]). In contrast, in-vivo models are based on a series of reaction equations modeling the coagulation cascade, from the exposure of tissue factor to the generation of thrombin [START_REF] Sorensen | Computational simulation of platelet deposition and activation: I. model development and properties[END_REF]; [START_REF] Kuharsky | Surface-mediated control of blood coagulation: the role of binding site densities and platelet deposition[END_REF]; [START_REF] Fogelson | Immersed-boundary-type models of intravascular platelet aggregation[END_REF]; [START_REF] Anand | A model for the formation, growth, and lysis of clots in quiescent plasma. a comparison between the effects of antithrombin iii deficiency and protein c deficiency[END_REF]; [START_REF] Seo | A coupled chemo-fluidic computational model for thrombogenesis in infarcted left ventricles[END_REF]; [START_REF] Yazdani | A general shear-dependent model for thrombus formation[END_REF]). Generally, in vivo models, prescribe an injury zone determined by the user producing tissue factor and thus, initiating the coagulation reactions [START_REF] Biasetti | An integrated fluid-chemical model toward modeling the formation of intra-luminal thrombus in abdominal aortic aneurysms[END_REF]; [START_REF] Seo | A coupled chemo-fluidic computational model for thrombogenesis in infarcted left ventricles[END_REF]). This type of approach cannot be used for device-related thrombosis due to the lack of a priori knowledge of the thrombotic zone. For this reason, the models dedicated to evaluate the device thrombotic performance have focused on platelet activity. For instance, [START_REF] Goodman | Computational model of device-induced thrombosis and thromboembolism[END_REF] developed a model to predict device-induced thrombosis based on platelet activation (agonist induced platelet activation as in [START_REF] Sorensen | Computational simulation of platelet deposition and activation: I. model development and properties[END_REF] and mechanical platelet activation by shear stress). Furthermore, the model accounts for platelet adhesion to artificial surfaces and thrombus growth by substantially increasing the fluid viscosity. [START_REF] Taylor | Development of a computational model for macroscopic predictions of device-induced thrombosis[END_REF] developed a macroscopic model that includes platelet activation and thrombus growth in regions related to high and low wall shear stress (WSS), respectively. Their model was derived from the work of [START_REF] Fogelson | Continuum models of platelet aggregation: Formulation and mechanical properties[END_REF] with some modifications to account for device-related thrombosis. While platelet activity has been extensively considered, other major players like the complement and contact activation systems must be included in the modeling, as pointed out by [START_REF] Gorbet | Biomaterial-associated thrombosis: roles of coagulation factors, complement, platelets and leukocytes[END_REF]. Moreover, considering contact activation could provide significant information to platelet models since one mechanism of their activation is due to thrombin as explained in [START_REF] Fogelson | Fluid mechanics of blood clot formation[END_REF], the latter being partially produced by reactions initiated by the contact system [START_REF] Larsson | A factor xiia inhibitory antibody provides thromboprotection in extracorporeal circulation without increasing bleeding risk[END_REF]). To the best of our knowledge, the only model considering the interaction of an idealized flow configuration and the contact activation system is the one of [START_REF] Gregory | An analysis of the contact phase of blood coagulation: Effects of shear rate and surface are intertwined[END_REF]. In their work, a steady-state solution of a kinetic model is analyzed. The coagulation reactions in this model include contact activation of factor XII and the common pathway of the coagulation cascade (until factor XI a ). The model results showed significant sensitivity to the flow dynamics and surface activity. A more recent kinetic model that considers the contact system in a platelet plasma kinetic network was developed by [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF]. Their model predicts the coagulation time of quiescent human blood samples without the action of tissue factor (TF). The results show that coagulation is due to the combined action of factors XII a and XI a . In parallel, recent clinical studies have pointed at the contact activation phase as a promising therapeutic target providing an anticoagulant strategy without increasing the risk of bleeding [START_REF] Jaffer | Medical device-induced thrombosis: what causes it and how can we prevent it?[END_REF]; [START_REF] Larsson | A factor xiia inhibitory antibody provides thromboprotection in extracorporeal circulation without increasing bleeding risk[END_REF]). A model that incorporates the contact activation system and the coagulation reactions until thrombin formation coupled with a realistic representation of flow dynamics is thus needed. To this respect, the objective of this study is to evaluate the potential of the contact system to initiate thrombin production. The role of the anionic phospholipids (platelet activity), fibrin formation and fluid-structure interaction of a growing thrombus were thus neglected, since these phenomena, although prominent in the general thrombus formation mechanism, are not specificaly involved in the surface production of thrombin. A surface-mediated thrombin production model is presented in a computational hemodynamics framework. The model accounts for interaction of the blood flow with the contact activation system and the rest of the blood coagulation cascade; it is detailed in Section 2, together with the numerical method. Furthermore, the model was applied to a device-characteristic flow configuration presented in the experimental work of [START_REF] Taylor | In vitro quantification of time dependant thrombus size using magnetic resonance image and computational simulations of thrombus surface shear stresses[END_REF]. To this end, a specific boundary condition was developed to mimic the contact activation at the device surface. The results discussed in section 3 illustrate the potential of the approach.

Materials and Methods

Governing equations

Fluid dynamics equations

Flow dynamics is described by the incompressible Navier Stokes equations that read

∇ • u = 0, (1)
ρ( ∂u ∂t + u • ∇u) = -∇p + µ∇ 2 u, ( 2
)
where u is the velocity field, p is the pressure field, ρ stands for the density and µ is the dynamic viscosity (both assumed constant-valued). The blood was considered as a Newtonian fluid (constant viscosity) since the objective of the present study was not focused on the effects due to Non-Newtonian behavior of blood.

Convection-Diffusion-Reaction equations

To account for the coagulation reactions, a set of evolution equations (one for each coagulation factor) is solved. The equation that captures the spatial and temporal evolution of the biochemical species is widely known as the convection-diffusion-reaction equation

∂C i ∂t + ∇ • (C i u) = ∇ • (D i ∇C i ) + R i , (3)
where C i (M = mol L -1 = 10 3 mol m -3 ) and D i (m 2 s -1 ) are the molar concentration and diffusion coefficients for each species i, R i (M s -1 ) is the reaction source term (computed according to the reactions rates presented in Table 1).

Coagulation model

A wide variety of models has been proposed to predict thrombin formation inside the blood stream [START_REF] Anand | A model incorporating some of the mechanical and biochemical factors underlying clot formation and dissolution in flowing blood[END_REF]; [START_REF] Leiderman | Grow with the flow: A spatial-temporal model of platelet deposition and blood coagulation under flow[END_REF]; [START_REF] Biasetti | An integrated fluid-chemical model toward modeling the formation of intra-luminal thrombus in abdominal aortic aneurysms[END_REF]). In those models, coagulation reactions start at the extrinsic pathway by the exposure of TF.

In the present work, the presence of an artificial surface in contact with blood is accounted for by including the contact activation system [START_REF] Basmadjian | Coagulation on biomaterials in flowing blood: Some theoretical considerations[END_REF]).

For this purpose the kinetic model of Chatterjee et al. ( 2010) is adopted. This set of reactions predicts thrombin production in the absence of exogenous TF and starts the coagulation reactions by the contact activation of factor XII. Figure 1 shows a schematic of the coagulation model used for the simulations including the contact activation system and the rest of the coagulation cascade reactions.

In Table 1, the 26 reactions used are listed with their respective reaction rates. The types of reaction kinetics are of first order, second order and Michaelis-Menten. Reaction rates for the different reactions are listed in three different columns according to their type. The source terms R i appearing in Eq. 3 are deduced from the coagulation reactions. It is important to mention that the 53 reaction rates values do not consider the platelet activity that is also reported by [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF].

Diffusion coefficients were estimated with the correlation of [START_REF] Young | Estimation of diffusion coefficients of proteins[END_REF]. The molecular weight of each chemical species is needed to estimate the diffusion coefficient along with the temperature and the dynamic viscosity. For the present study a temperature of T = 25 • C and a dynamic viscosity of µ = 3.71×10 -3 Pa s were used. The correlation used considers a partial specific volume of 0.73 cm 3 g -1 which is assumed to be common for all the species. The expression proposed by Young et al writes:

D = 8.34 × 10 -8 (T /µM 1/3 ) (4)
The diffusion values estimated with Eq. 4 are extremely small (of the order of 10 -11 m 2 s -1 ). If these values were to be used, numerical instabilities and large diffusion characteristic times would render the convergence to a steady state solution extremely difficult. This problem has been emphasized in studies using numerical models for coagulation. This difficulty is either overcome by refining the mesh which is hardly feasible except for 1D simulations, or by increasing artificially the diffusion coefficients as in [START_REF] Biasetti | An integrated fluid-chemical model toward modeling the formation of intra-luminal thrombus in abdominal aortic aneurysms[END_REF]. Here, we follow the latter approach and multiply all diffusion coefficients by a factor of 1000.

The main objective of this work is to evaluate the potential of the contact activation system to initiate thrombin formation. Consequently, the model did not include inhibitors, TF and fibrin-related reactions that are included in [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF]. Finally the model assumes an infinite number of binding sites (infinite platelet supply), as in [START_REF] Biasetti | An integrated fluid-chemical model toward modeling the formation of intra-luminal thrombus in abdominal aortic aneurysms[END_REF]. The reader is referred to [START_REF] Kuharsky | Surface-mediated control of blood coagulation: the role of binding site densities and platelet deposition[END_REF]; [START_REF] Anand | A model incorporating some of the mechanical and biochemical factors underlying clot formation and dissolution in flowing blood[END_REF] for a discussion on the role of a finite number of binding sites. It is important to note that the scope of this article is not to model the complete and complex dynamics of thrombus formation. A full model must include the role of anionic phospholipids, platelet adhesion and aggregation as in [START_REF] Leiderman | Grow with the flow: A spatial-temporal model of platelet deposition and blood coagulation under flow[END_REF]. Other important players for a robust model should be included such as fibrin production [START_REF] Neeves | Thrombin flux and wall shear rate regulate fibrin fiber deposition state during polymerization under flow[END_REF]), coagulation inhibitors [START_REF] Zhu | Mathematical modeling of blood coagulation cascade: kinetics of intrinsic and extrinsic pathways in normal and deficient conditions[END_REF]; [START_REF] Rice | Differential contributions of intrinsic and extrinsic pathways to thrombin generation in adult, maternal and cord plasma samples[END_REF]) or thrombus growth (solid-fluid interaction) as model by [START_REF] Ngoepe | Computational modelling of clot development in patient-specific cerebral aneurysm cases[END_REF]. Such a comprehensive model may provide a good prediction tool to contribute to the optimization of medical devices. However, the potential of the contact activation system to initiate coagulation reactions is better assessed by focusing on this mechanism alone, without including all the other phenomena relevant to thrombus formation. This is the methodology followed in this study.

A boundary condition for the contact activation

To model the contact activation of factor XII, a wall boundary condition is derived from the species mass balance at the wall where protein adsorption takes place followed by contact activation. The balance is obtained by integrating Eq. 3 over a rectangular control volume Ω (which includes a volume of fluid and a volume of wall, see Fig. 2) of surface Σ. Using the divergence theorem to reformulate the conservative terms and integrating over Ω leads to:

XII XI P K XIIa V II V IIa XIa Kal IX IXa V IIIa V III IXa\V IIIa X Xa Va V Xa\Va Prothrombin Thrombin Contact Activation System
# Reaction Forward M -1 s -1 Reverse s -1 Forward s -1 Xa + V II → Xa + V IIa k 1 = 1.3 × 10 7 IIa + V II → IIa + V IIa k 2 = 2.3 × 10 4 II + Xa → IIa + Xa k 3 = 7.5 × 10 3 IIa + V III → IIa + V IIIa k 4 = 2.0 × 10 7 V IIIa + IXa ↔ IXa = V IIIa k 5 = 1.0 × 10 7 k 6 = 5.0 × 10 -3 IXa = V IIIa + X ↔ IXa = V IIIa = X → IXa = V IIIa + Xa k 7 = 1.0 × 10 8 k 8 = 1.0 × 10 -3 k 9 = 8.2 V IIIa ↔ V IIIa 1 • L + V IIIa 2 k 10 = 6.0 × 10 -3 k 11 = 2.2 × 10 4 IXa = V IIIa = X → V IIIa 1 • L + V IIIa 2 + X + IXa k 12 = 1.0 × 10 -3 IXa = V IIIa → V IIIa 1 • L + V IIIa 2 + IXa k 13 = 1.0 × 10 -3 IIa + V → IIa + V a k 14 = 2.0 × 10 7 Xa + V a ↔ Xa = V a k 15 = 4.0 × 10 8 k 16 = 0.2 Xa = V a + II ↔ Xa = V a = II → Xa = V a + mIIa k 17 = 1.0 × 10 8 k 18 = 103 k 19 = 63.5 Xa = V a + mIIa → Xa = V a + IIa k 20 = 1.5 × 10 7 14* XII → XIIa ks = 5.0 × 10 -7 ms -1 XIIa + XII ↔ XIIa = XII → XIIa + XIIa k 22 = 1 × 10 8 k 23 = 750 k 24 = 3.3 × 10 -2 XIIa + P K ↔ XIIa = P K → XIIa + K k 25 = 1 × 10 8 k 26 = 3.6 × 10 3 k 27 = 40 XII + K ↔ XII = K → XIIa + K k 28 = 1 × 10 8 k 29 = 45.3 k 30 = 5.7 P K + K → K + K k 31 = 2.7 × 10 4 XI + IIa ↔ XI = IIa → XIa + IIa k 32 = 1 × 10 8 k 33 = 5 k 34 = 1.3 × 10 -4 XIIa + XI ↔ XIIa = XI → XIIa + XIa k 35 = 1.0 × 10 8 k 36 = 200 k 37 = 5.7 × 10 -4 XIa + XI → XIa + XIa k 38 = 3.19 × 10 6 XIa + IX ↔ XIa = IX → XIa + IXa k 39 = 1.0 × 10 8 k 40 = 41.0 k 41 = 7.7 IXa + X ↔ IXa = X → IXa + Xa k 42 = 1.0 × 10 8 k 43 = 0.64 k 44 = 7.0 × 10 -4 Xa + V III ↔ Xa = V III → Xa + V IIIa k 45 = 1.0 × 10 8 k 46 = 2.1 k 47 = 0.023 V IIa + IX ↔ V IIa = IX → V IIa + IXa k 48 = 1.0 × 10 8 k 49 = 0.9 k 50 = 3.6 × 10 -5 V IIa + X ↔ V IIa = X → V IIa + Xa k 51 = 1.0 × 10 8 k 52 = 210 k 53 = 1.6 × 10 -6
Table 1 Coagulation reactions. Kinetic parameters and coagulation reactions used in the simulations. Reaction 14 is taken into account only at the surface by means of the wall boundary condition. Parameter k s units (m s -1 ) correspond to a surface reaction; its value was obtained by adapting the volume reaction rate from [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF].

Ω ∂C i ∂t dΩ + Ω -R i dΩ = Σ (D i ∇C i ) • ndΣ - Σ C i u • ndΣ (5)
The reaction source term R i can be decomposed into the volume (blood-phase) reaction source term ωi and the surface source term ṡi multiplied by a Dirac function δ (which is no zero at the device surface). The mass balance becomes:

Ω ∂C i ∂t dΩ + Ω -( ωi + δ ṡi )dΩ = Σ (D i ∇C i ) • ndΣ - Σ C i u • ndΣ (6)
The surface Σ can be separated in lateral, inferior and superior surfaces respectively Σ = 2Σ l + Σ in + Σ s . If the height h tends to zero (h → 0) the volume integrals over Ω and the surface integrals over the lateral surfaces Σ l are equal to zero; the inferior surface Σ in collapse with the device wall and the surface reaction source term can be transformed in a surface integral over the device wall thanks to the Dirac function. In addition the convection fluxes at the device surface are equal to zero (assuming the wall is at rest and impermeable), thus Eq.5 can be rewritten as:

-

Σs (D i ∇C i ) • n s dΣ s = Σin (D i ∇C i ) • n in dΣ in + Σin ṡi dΣ in (7)
The diffusion flux inside the solid domain is zero since no species is present in the device itself; thus one can write:

-

Σs (D i ∇C i ) • n s dΣ s = Σin ṡi dΣ in (8)
A boundary condition can be obtained by applying Eq. 8 to factor XII a . The surface source term is obtained assuming an infinitely fast adsorption and an infinite number of binding sites for factor XII, thus the reaction that takes places is the conversion of the near-wall factor XII (indexed by w ) to XII a . Finally the following wall boundary condition is obtained:

D XIIa ∂C XIIa ∂n = -k s C XIIw (9)
where k s is a material property (surface reaction rate) and the direction n is the normal direction pointing to the blood-phase of the domain. In this study the value published in [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF] for the activation of factor XII in a perfectly mixed environment (volumetric reaction rate) was converted to a surface reaction rate k s . From a straightforward dimensional analysis, one obtains

k s ∼ k v V batch
S batch where V batch and S batch are the characteristic volume and the active surface of the experimental set up respectively. The surface reaction rate k s = 5.0 × 10 -7 ms -1 is retrieved by using V batch ∼ 10 µL and S batch ∼ 0.1 cm 2 following the experimental values reported in [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF].A range of values was also considered to investigate the effect of this physical quantity on the results (see section 3).

Numerical implementation and verification

Equations 1 and 2 were solved with a fourth-order finitevolume method on an unstructured grid, implemented in the in-house YALES2BIO 1 solver. YALES2BIO is based on a massivel parallel finite volume solver for incompressible flows [START_REF] Moureau | Design of a massively parallel CFD code for complex geometries[END_REF]) and has been extensively validated for several biomedical applications at micro and macro scale (Chnafa et al. (2014b,a) 2017)). It uses a spatial fourth-order scheme with an explicit four-step Runge-Kutta scheme in time. The pressure term is handle with the Chorin's projectioncorrection method. Equation 3 was also solved using a fourth-order scheme in space and time for each species. Several test cases were developed to verify the implementation of Eq. 3 in YALES2BIO. Test cases with reactions of the Michaeles-Menten type [START_REF] Michaelis | The original michaelis constant: translation of the 1913 michaelismenten paper[END_REF]), as well as first and second order reactions were compared against analytic solutions. Furthermore, the results of the kinetic model introduced by [START_REF] Hockin | A model for the stoichiometric regulation of blood coagulation[END_REF] were reproduced. In this test case the reactions start with an initial condition of exogenous TF. This allows to validate the different types of reactions that are shared by the different coagulation kinetic schemes; for instance, reactions 1 to 13 in Table 1 are shared with the scheme of [START_REF] Hockin | A model for the stoichiometric regulation of blood coagulation[END_REF].

Moreover, a simple 1D diffusive test case involving a volume and a surface reaction (contact activation boundary condition) without convection was specifically designed during the course of this study. In this academical configuration illustrated in Fig. 3, a species A can be either activated at the bottom wall and become A a (surface reaction rate k s ) or transformed into species B in the flow domain (volume reaction rate k v ), whereas the upper wall acts like an infinite reservoir at concentration C A (h) = α. In the calculations the domain is 2D, with periodic conditions in the directions tangential to the walls. The analytic steady state solution for the concentration of species A reads:

C A (y) = α e β(h-y) (k s (e 2βy -1) + βD A (1 + e 2βy )) k s (e 2βh -1) + βD A (1 + e 2βh ) ( 10
)
where k s and k v are the surface and volume reaction rates, respectively, β = kv D A and h is the channel height. Table 2 shows the three different operating points that were used to verify the results.

Cases k s m s -1 k v s -1 Case 1 1.0 × 10 -4 1.0 × 10 4 Case 2 1.0 × 10 -4 1.0 × 10 3 Case 3 1.0 × 10 -2 1.0 × 10 2
Table 2 Operating points. The different set of parameters (k s and k v ) used for the verification of the code. The corresponding steady state solutions are shown in Fig. 3.

Figure 3 shows the steady state solution of species A for both the numerical and analytic cases. An excellent agreement is obtained for all of the operating points. , Case 3 ).

Backward facing step configuration

Numerical simulations are performed in a backward facing step (BFS) configuration presented in Fig. 4, advancing the set of reactions presented in Table 1. This type of configuration involves a flow separation that several medical devices could feature and which has been related to thrombus formation [START_REF] Hastings | Sources of circuit thrombosis in pediatric extracorporeal membrane oxygenation[END_REF]; [START_REF] Taylor | In vitro quantification of time dependant thrombus size using magnetic resonance image and computational simulations of thrombus surface shear stresses[END_REF]). The step height h s = 2.5 mm and vessel height H = 10 mm match the experimental dimensions presented by [START_REF] Taylor | In vitro quantification of time dependant thrombus size using magnetic resonance image and computational simulations of thrombus surface shear stresses[END_REF].

The inlet and outlet conditions were placed at a distance of 10H to avoid boundary effects on the region of interest located downstream of the step: the expected recirculation zone for the flow regime considered (Reynolds number Re = 430 , based on the inlet channel height and the upstream bulk velocity) is Xr = 8.9h s according the results of [START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF]. The kinematic viscosity was set to a constant value ν = 3.5 × 10 -6 m 2 s -1 .

Boundary and initial conditions. The inlet velocity was set to U x = 0.2 m s -1 , consistently with the experimental conditions described in [START_REF] Taylor | In vitro quantification of time dependant thrombus size using magnetic resonance image and computational simulations of thrombus surface shear stresses[END_REF]. Table 3 shows the initial and inlet boundary conditions factors that were taken from [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF] and which correspond to physiological values. The boundary and initial conditions for the rest of the factors were set to zero. One of the most important features of the simulation is the reactive surface boundary condition, Eq. 9.

The boundary condition for the contact activation of factor XII was applied at the entire wall of the computational domain (see Fig. 4).

Factor Concentration (M) Spatial and temporal resolution. Several grids were used to ensure that the simulations yield results independent from the grid used in the calculation. Table 4 shows the representative cell size and the number of elements used for the computations. The simulation time step is given by the Courant-Friedrichs-Levy stability condition which is define using the time step ∆t, the local velocity field u and the mesh size ∆x as CF L = u∆t ∆x ≤ 1. The value used in the current simulations is CF L = 0.9 which corresponds to a time step of the order 10 -4 s. Table 4 Mesh properties. ∆ h is the representative cell size computed as in [START_REF] Celik | Procedure for estimation and reporting of un-certainty due to discretization in cfd applications[END_REF] and the total number of grid elements (triangles) are displayed.

V II 1.0 × 10 -8 V IIa 1.0 × 10 -10 X 1.6 × 10 -7 IX 9.0 × 10 -8 II 1.4 × 10 -6 V III 7.0 × 10 -10 V 2.0 × 10 -8 XII 3.4 × 10 -7 P K 4.5 × 10 -7 XI 3.1 × 10 -8
3 Results: Backward facing step

Flow Dynamics

The velocity field reached a steady state solution before 1.0 s as shown in Fig. 5 for different spatial point probes located inside and outside the recirculation zone. ) Figure 6 shows the steady state velocity field as well as the streamlines that allow the visualization of recirculation region after the step. The reattachment length X r is in good agreement with the numerical and experimental results published in [START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF]; Velocity Magnitude (m s -1 ) Fig. 6 Velocity magnitude with white stream lines. The reattachment length X r is shown. [START_REF] Armaly | Experimental and theoretical investigation of backward-facing step flow[END_REF] which lead to X r = 8.9h. The comparison was performed for all the grids used in this study. The percentage error of the coarser mesh was e rel = 3.1% and for the medium and fine mesh of about e rel = 2.9% which is a fair approximation considering the main objective of the current study.

P 1 U x P 1 U y P 2 U x P 2 U y P 3 U x P 3 U y

Coagulation cascade

In order to quantify the outcome of the coagulation cascade, the value of total thrombin is considered (II a + 1.2mII a is a quantitative measure of thrombin activity in the system, according to [START_REF] Hockin | A model for the stoichiometric regulation of blood coagulation[END_REF]).

First, a mesh convergence study was performed using the mean value of the total thrombin in the computational domain. The methodology published in [START_REF] Celik | Procedure for estimation and reporting of un-certainty due to discretization in cfd applications[END_REF] was use. The convergence grid index is used to quantify the distance of a computed value (simulation result) to the asymptotic value estimated by an extrapolation method using the different levels of discretization. For the finest mesh the convergence grid index obtained was CGI f inest = 0.41% with a relative extrapolated error of e ext = 0.33%, the quantity used for the analysis was the total thrombin concentration value of a point probe inside the recirculation region. These values ensure that the total thrombin concentration field would not substantially change with even a finer mesh. The surface rate used for the mesh convergence study was k s = 5.0 × 10 -4 m s -1 , in order to reduce the computational cost of the study.

Figure 7 shows the time evolution of total thrombin at point probes P 1 , P 2 and P 3 (see Fig. 4 for the location of the probes). The time response is characterized by a lag time of about 100 s. Then, a propagation phase is observed until a small peak which evolves into a more stable plateau close to the inlet value of II = 1.4 µM. time (s) IIa + mIIa (M)

P 1 P 2 P 3
Fig. 7 Time evolution of II a + 1.2mII a for each point probe (see Fig. 4). P 1 is near the wall region, P 2 inside the recirculation zone and P 3 outside the recirculation region.

Figure 8 shows the concentration field evolution of total thrombin at different time values. It can be observed that the region of largest total thrombin concentration is located inside the recirculation zone downstream of the step. Large concentrations of total thrombin extend further downstream at the lower wall.

The local concentration of coagulation factors inside the recirculation region is shown in Fig. 9. The plot is extracted along a profile located at the half height of the step (h s ), y = 1.25 mm. It can be observed that the largest concentration of total thrombin is located at the step wall.

Thrombin concentration levels at the recirculation region are large enough to form a stable thrombus (Biasetti et al. ( 2012)). Furthermore, the zone with the largest concentration of thrombin aligns well with the experimental zone in which thrombus is formed [START_REF] Taylor | In vitro quantification of time dependant thrombus size using magnetic resonance image and computational simulations of thrombus surface shear stresses[END_REF]). In macroscopic applications, flow metrics such as the wall shear stress are commonly used to evaluate the thrombogenic performance of medical devices. Such flow quantities are practical to use due to their low computational cost compared to a reacting flow study. It is thus interesting to investigate how the outcome from the contact system and the wall shear stress correlate. Figure 10 shows a scatter plot of the wall shear stress against the near wall concentrations of factors XII a and XI a . The quantities were normalized by the maximum wall shear stress and maximum concentrations of species XII a and XI a . Factor XII a is clearly not correlated to WSS. Even though a clear correlation is not observed neither for XIa, this factor takes its largest values in areas where WSS is small. At the end, these results demonstrate that accounting for the contact system generates a physical behaviour which cannot be properly represented by WSS-based flow analysis. 3.3 Parametric studies

Influence of surface kinetics k s

As already mentioned, there is a lack of reliable data about the surface reaction rate (k s contact activation) and this quantity may depend on the material type.

The nominal value (k s = 5.0 × 10 -7 m s -1 ) was obtained from the conversion of the volume reaction rate published in [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF], as previously explained, using a dimensional analysis that involves the characteristic volume and the active surface from the experimental set up reported in [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF].

To investigate how the outcome of the simulations depends on this parameter, a parametric study was performed. Six values were considered in the range k s = [5.0 × 10 -5 , 5.0 × 10 -10 ] m s -1 spanning six orders of magnitude. Figure 11 shows the evolution at point probe P 1 of factors XII a , XI a and II a + 1.2mII a over time. In addition, inserts in the plots display the variation of specific times of interest with respect to k s .

For II a + 1.2mII a and XII a , the inflection point of the time evolution graphs is extracted by fitting a hyperbolic tangent function. For factor XI a , the insert shows the time at which the maximum concentration is reached. The thrombin generation starts earlier as the parameter k s is increased. However, while five orders of magnitude are explored for the value of k s , the impact on the thrombin generation time is moderate and the dependence appears to be logarithmic. In addition, the total thrombin concentration fields (steady state solution) were compared, showing no significant difference. It can be observed from the evolution of factors XII a and XI a that the time delay present in the dynamics of these factors follow the same trend as the one observed for II a + 1.2mII a . In other words, the initial delay in the activation of factor XII is directly responsible for the delay in thrombin generation. While the influence of k s on the time needed for generation of factor XII a is very large, this time has to be compared with the total time needed for the rest of the coagulation cascade, which is 10 times larger in the present case. It is important to note that when k s was set to zero the thrombin concentration levels remained close zero (not shown).

Parametric study: Influence of factor XII inlet boundary condition.

A parametric study on the influence of of the concentration of factor XII at the inlet of the domain was also performed. This allows to evaluate the variability of the results to different samples with deficient factor XII. Table 5 shows the four different concentration values that were used for this purpose.

Case XII (µM) 1 0.306 2 0.272 3 0.17 4 0.102

Table 5 Factor XII inlet concentrations used in the parametric study. These values correspond to 90%, 70%, 50% and 30% of the nominal concentration reported in Table 3.

Figure 12 shows the time evolution of factors XII a , XI a and II a + 1.2mII a at the point probe P 1 . For each case, factor XII a reaches a plateau equal to the inlet concentration of factor XII, as expected. The results show longer lag times for thrombin propagation as the concentration of factor XII is decreased. Interestingly the maximum production of factor XI a is increased as factor XII concentration is reduced. This can be explained by the fact that a larger proportion of the available factor XI is activated by factor II a (reaction 19) in this case, instead of being generated by factor XII a (reaction 20). This behavior is confirmed by the time evolution of the source term of reactions 19 and 20 (not shown) and is also coherent with the reverse rate of reaction 19 (k 33 in Table 1) being very small compared to the reverse rate of reaction 20 (k 36 in Table 1). The results suggest that a small amount of activated factor XII (contact activation) is sufficient to initiate the thrombin propagation phase.

Discussion and Conclusion

Previous studies have modeled device-related thrombosis with an approach focussing on platelet activity [START_REF] Goodman | Computational model of device-induced thrombosis and thromboembolism[END_REF]; [START_REF] Taylor | Development of a computational model for macroscopic predictions of device-induced thrombosis[END_REF]). However, the contact activation system is a key component of the thrombotic network [START_REF] Larsson | A factor xiia inhibitory antibody provides thromboprotection in extracorporeal circulation without increasing bleeding risk[END_REF][START_REF] Jaffer | Medical device-induced thrombosis: what causes it and how can we prevent it?[END_REF]) which up to date has not been included in a complete model of thrombus formation. The current model aims to serve as a first step towards the integration of the contact activation in a complete macroscopic ), XII 50% = 0.17 µM ( ), XII 30% = 0.102 µM (

).

thrombus formation model specific to medical devices.

The current model was applied to a backward facing step configuration typical of medical devices. One of the main features of the model is that it allows the initiation of the coagulation reactions without an a priori knowledge of the thrombotic zone. The results show that a significant amount of thrombin is generated in the recirculating region behind the backward facing step. This could be explained by comparing the characteristic residence and chemical times. Large residence time at the recirculating regions enhances the thrombin production reactions. Qualitatively, the model leads to higher total thrombin concentrations in regions where the thrombus formation was observed experimentally.

The results also show a poor correlation between wall shear stress and near wall concentrations of factor XI a , and no correlation for factor XIIa that is present all along the device surface. This suggests that the WSS should be complemented with other flow metrics that quantify diffusive and convective transports to characterize the thrombogenic performance of medical devices.

Molecular diffusivity of the species has been modified for numerical reasons (from D i ∼ 10 -11 m 2 s -1 to D i ∼ 10 -8 m 2 s -1 ), because the proper resolution of the species boundary layer would make the calculations prohibitive. However, it is most probably that if this modification has a consequence, it is only in the sense of underestimating the thrombin generation by the contact system. First, in the main stream, convection clearly dominates, the Peclet number is of the order 10 7 and 10 4 for the molecular diffusivity of D i ∼ 10 -11 m 2 s -1 to D i ∼ 10 -8 m 2 s -1 respectively. In contrast, in the recirculation zone the results are controlled by the competition between reactions and diffusion. Large D i values actually underestimate the local concentration of chemical species since they are diffused outside the recirculation region faster. The higher the diffusion coefficients, the lower the concentration of thrombin in the recirculation region. This is illustrated by two additional simulations with diffusion coefficient values of the order D i ∼ 10 -7 m 2 s -1 and D i ∼ 10 -6 m 2 s -1 , which were performed to compare the concentration values inside the recirculation region. Figure 13 shows the underestimation of thrombin as the molecular diffusivity is augmented. IIa + 1.2mIIa M Fig. 13 Maximal total thrombin concentration inside the recirculation region for three orders of magnitude of the diffusion coefficients (D i ∼ 10 -8 , 10 -7 , 10 -6 m 2 s -1 ).

An interesting feature of the present model is that it includes a material property k s (typical to each medical device) which is a novel feature for thrombus predictive models representing the kinetics of the contact sur-face reaction. A parametric study was performed on k s showing a moderate influence of this parameter on the total thrombin generation, although it strongly modifies the generation time of factor XII a . This is because the characteristic time of this later mechanism (reaction 14 in Table 1) is very small compared to the whole time required to generate thrombin. In the end, the current study focuses on the contact system as one potential player in thrombin production. This approach could complement more complex models that consider platelet activity, fibrin formation and fluid structure interaction of the growing thrombus, providing a proper initiation of the coagulation reactions in medical devices.

Fig. 1

 1 Fig. 1 Schematic of the coagulation reaction network. The kinetic network is inspired from the platelet-plasma model of Chatterjee et al. (2010) without considering TF, inhibitors and fibrinogen reactions. The dashed lines correspond to the feedback loop reactions related to thrombin.

Fig. 2

 2 Fig. 2 Mass balance around the device wall. Ω Control of volume, Σ Surface envelope of the control volume, Sub index in, s and l stand for the inferior, superior and lateral surfaces respectively, h height of the control volume.

  ; Mendez et al. (2014); Sigüenza et al. (2017); Chnafa et al. (2016); Lanotte et al. (2016); Zmijanovic et al. (

Fig. 3

 3 Fig. 3 1D reaction-diffusion test case. C A = α = 1.0 µM, initial volume concentration of all the species C i = 0 M, diffusion coefficient of species D A = D B = 1.0 × 10 -2 m 2 s -1 , channel height h = 0.01 m. Numerical results (Case 1 , Case 2 , Case 3 ) and analytic Eq. 10 data (Case 1 , Case 2, Case 3 ).

Fig. 4

 4 Fig. 4 2D Backward facing step computational domain. P i are the locations of the probes used in Figs. 5, 7 and Fig. 11 (P 1 = (0.5, 0.5), P 2 = (4.0, 1.5), P 3 = (4.0, 4.0) mm); the step and channel heights match the experimental values of Tayloret al. (2014). The thick black line denotes the solid boundary where boundary condition Eq. 9 is applied and activation of factor XII is allowed.

Fig. 5

 5 Fig. 5 Temporal evolution of x and y velocity components in different point probes.

Fig. 8 Fig. 9

 89 Fig. 8 Total thrombin concentration field at t = 100, 125, 150 and 175 s. The black dashed line represents over which total thrombin is plotted in Fig. 9. The red line (t = 100) represents approximately the experimental wall region of thrombus formation of Taylor et al. (2014).

Fig. 10

 10 Fig. 10 Scatter plot of wall shear stress and wall concentrations of factors XII a and XI a .

Fig. 11

 11 Fig. 11 Time evolution at point probe P 1 of total thrombin (upper), activated factor XII (medium) and activated factor XI(bottom). k s = 5 × 10 -5 ( ), k s = 5 × 10 -6 ( ), k s = 5 × 10 -7 ( ), k s = 5 × 10 -8 ( ), k s = 5 × 10 -9 ( ), k s = 5 × 10 -10 ( ).The insets of II a + 1.2mII a and XII a factor show the Inflection Point Time (IPT) for each k s value. For factor XI a , the inset displays the maximum concentration time (MCT) as a function of k s .

Fig. 12

 12 Fig.12Time evolution of point probe P 1 for total thrombin (upper), activated factor XII (medium) and activated factor XI (bottom) at point probe P 1 . The inlet and initial concentrations are XII 90% = 0.306 µM ( ), XII 70% = 0.272 µM (), XII 50% = 0.17 µM ( ), XII 30% = 0.102 µM ().

Table 3

 3 Initial and inlet factor concentrations. Physiological values used in the simulations taken from Chatterjee et al. (2010).

http://www.math.univ-montp2.fr/∼yales2bio/

Acknowledgements We are grateful to the CONACyT, Mexico scholarship and the LabEx Numev (convention ANR-10-LABX-20) for their financial support. Drs. Vincent Moureau and Ghislain Lartigue from CORIA laboratory, and the SUC-CESS scientific group are acknowledged for providing the YALES2 solver which constitutes the basis of the YALES2BIO tool. Professor Jean-François Schved and MD. Muriel Giansily from the University Hospital of Montpellier are thanked for helpful discussions about the biological aspects of thrombosis. This work was performed using HPC resources from GENCI-CINES (Grants 2016-c2015037194 and 2017-A0020307194) and with the support of the High Performance Computing Platform HPC@LR.

Conflict of Interest:

The authors declare that they have no conflict of interest.