Hasan Misaii
email: hasan.misaii@utt.fr

Firoozeh Haghighi
email: fhaghighi@ut.ac.ir

Mitra Fouladirad
email: mitra.fouladirad@centrale-marseille.fr

Multi-component system maintenance optimisation and masked data

The paper discusses he maintenance optimization of a multi-component system. The exact cause of failure is unknown and only some information about it is available. Opportunistic perfect preventive maintenance is considered and optimized using cost criteria. Numerical implementations illustrate the applicability of results.

I. INTRODUCTION

Maintenance planning for multi-component systems is an important issue and not easy to address because of economic, structural and/or failure dependency between components [START_REF] Li | A condition-based maintenance policy for multi-component systems with Levy copulas dependence[END_REF], [START_REF] Thomas | A survey of maintenance and replacement models for maintainability and reliability of multi-item systems[END_REF]. Maintenance policy optimization depends on the available information at our disposal, this paper discusses maintenance policy optimization in presence of incomplete information. In the case of multi-component systems, the maintenance can be perfect and maintain all the components which is very costly or imperfect in sens that only failed components are replaced. Imperfect maintenance can be applied when the system or the component is not as good as new after maintenance. Sometimes, the system is only repaired in order to continue to operate where the failure rate remains unchanged in comparison to the failure rate before maintenance, that is called minimal maintenance [START_REF] Nakagawa | Maintenance theory of reliability[END_REF]. Moreover, preventive maintenance actions can be planned in order to avoid failure. The preventive action occurs before the failure an is more cost-efficient since it permits to avoid the failure and its induced cost and consequences. Preventive maintenance planning is an important issue since a frequent preventive maintenance could be costly and by reducing the frequency of these actions, the occurence of failures will be more likely which will induce higher costs and a period of unavailability. That is why preventive maintenance scheduling has attracted lot of attentions lately as well as in high-tech industry as in basic production and manufacturing. The preventive maintenance can be planned periodically called period maintenance, based on the age of the system called age-based maintenance, randomly according to the maintenance teams availability or at failure of some components called opportunistic maintenance or based on the health indicator of the component or the system called conditionbased or predictive maintenance, for more details refer to [START_REF] Nakagawa | Maintenance theory of reliability[END_REF], [START_REF] Ilya | Reliability theory: with applications to preventive maintenance[END_REF], [START_REF] Wang | Condition-based maintenance modelling[END_REF], [START_REF] Mobley | An introduction to predictive maintenance[END_REF], [START_REF] Vu | Dynamic opportunistic maintenance planning for multi-component redundant systems with various types of opportunities[END_REF], [START_REF] Yuan | An improved opportunity maintenance model of complex system[END_REF], [START_REF] Misaii | Opportunistic perfect preventive maintenance policy in presence of masked data[END_REF].

The choice of maintenance planning is strongly related to the available data and information. In absence of regular monitoring of a health indicator or sensor measurements on this latter, it is not possible to plan condition-based or predictive maintenance actions. If only historical data is available, periodic or age-based maintenance is carried out. In the case of multi-component systems where it is not sensible to replace the whole system at each maintenance, at components failures some diagnosis can be made on other components and carry out opportunistic maintenance based on observations at maintenance date. However, the cause of the failure of a multi-component system (the failed component(s)) could be not always identified. In this framework, the available data is called masked data. Dealing with this kind of data in the framework of reliability calculation and maintenance planning is a big challenge, [START_REF] Sen | Analysis of masked failure data under competing risks[END_REF], [START_REF] Usher | An iterative approach for estimating component reliability from masked system life data[END_REF]. Since it is important to take advantage of available information, in the framework of masked data, opportunistic maintenance polices could reduce maintenance costs.

In this paper, an opportunistic perfect preventive maintenance is considered in presence of masked data. The preventive maintenance action replaces the failed component as good as new. When the system fails besides the failed component, some endangered components are replaced by a new one thus the perfect preventive replacements are opportunistic. Inspections are periodic and inter-inspection interval is optimized. The rest of the paper is organized as follows. In section 2, the model is explained. The maintenance model is presented in section 3. A numerical example is conducted in order to illustrate the applicability of the proposed method in section 4. Finally, the conclusion is given in section 5.

II. MODEL DESCRIPTION

A series system with J components is considered. When the system fails we observe failure time, t, but the exact cause of failure might be unknown but it is known that it belongs to a subset of {1, 2, ..., J} called MRS. Let M be the MRS corresponding to the failure time t for the system. The set M contains components that are possible to be the cause of system failure and if M = {1, ..., J} then the system cause of failure is called to be completely masked.

Let T l ; l = 1, 2, ..., J be the lifetimes of the lth component (independent components) and assume that the system fails only due to one of the J components, therefore the system failure time T is defined to be T = min(T 1 , ..., T J ). Let be f l (t) and R l (t) probability density and reliability functions denoted by , respectively. The reliability function of T is given by

R(t) = R(t; θ) = P θ (T > t) = J l=1 [1 -F l (t)]
θ = (θ 1 , ..., θ J ) and θ l is the set of parameters related to the lth component and F l is corresponding distribution function of lth component. Let's consider the cause of failure K as random variable. Then the joint probability density function of (T, K) is given by

f T,K (t, l) = f l (t) j≠l [1 -F j (t)] (1)
and also F (j, t) = P(K = j, T ≤ t), or equivalently R(j, t) = P(K = j, T > t) [START_REF] Crowder | Multivariate survival analysis and competing risks[END_REF]. The probability

p t l (M i ) = P(M = M i T = t, K = l
) is called the masking probability, where M i is an observation of M . Some authors such as Mukhopadhyay [START_REF] Mukhopadhyay | Bayesian Analysis of Masked Series System Lifetime Data[END_REF], Kuo and Yang [START_REF] Kuo | Bayesian reliability modeling for masked system lifetime data[END_REF] and Cai and et al. [START_REF] Cai | Robust Bayesian analysis for parallel system with masked data under inverse Weibull lifetime distribution[END_REF], assumed

P(M = M i T = t, K = l) = P(M = M i K = l) = p l (M i ),
that is, the masking probability is independent of failure time, but is dependent to the causes of failure. Similar to a new approach that was presented to model the dependency of the masking probability on the failure time and its exact cause using the multinomial logistic regression model [START_REF] Misaei | Analysis of masked competing risks data with cause and time dependent masking mechanism[END_REF], we also assume that the masking probability depends on the failure time and its exact cause.

Suppose M be the set of all nonempty subsets of {1, ..., J} that have 2 J -1 members. For l = 1, ..., J, define by

M l = {M ∈ M ∶ l ∈ M } the elements of M that include l, thus p t l (M i ) = P(M = M i K = l, T = t) = 0 ∀M i ∈ M c l = M -M l and Mi∈M p t l (M i ) = Mi∈M l t l (M i ) = 1, l = 1, ..., J denote by P l = {p t l (M i ) ∶ M i ∈ M l }, l = 1, 2, .
. . , J then the set of all masking probabilities is P = (P 1 , ..., P J ).

The reliability function of the multi-component system is given by

R(t) = ∞ t J l=1 f l (u) j≠l R j (u)du.
This function depends on the reliability of components and therefore on lifetime distribution of the components which are unknown.

III. MAINTENANCE MODELING

The hypothesis of the maintenance framework are as follows.

• The system is inspected periodically at times kτ ; k = 1, 2, ..., with cost c ins for the system and M k the corresponding masked sets. • The time interval ((k -1)τ, kτ ] is called the kth period.

• Inspections are carried out at the end of each period • The time required for inspection and maintenance actions is negligible.

• The system failure is not self-announced • The components are maintained independently.

• At the kth inspection time, kτ , a maintenance action is performed if the system has been failed during ((k -1)τ, kτ ] interval. The probability of each cause in M k given possibly masked set and interval censored failure time is given by

p jM k = P(K = j M k , u ∈ ((k -1)τ, kτ ]) = ∫ kτ (k-1)τ P(M k j, u)f T,K (u, j)du ∫ kτ (k-1)τ ∑ l ′ ∈M k P(M k l ′ , u)f T,K (u, l ′ )du
where u is the exact failure time. Note that p jM k = 0 for j ∉ M k . • If the system fails at ((k-1)τ, kτ ] a maintenance action is carried out for each component in M k according to a predetermined value of ρ; 0 < ρ < 1, as follows:

-If T l > (k -1)τ and T l < kτ then perfect corrective maintenance (PCM) action is done for component l with cost c lc (that is, the failed component l is replaced by a new one). -If T l > kτ and p lM k > ρ then opportunistic perfect preventive maintenance (OPPM) action is done for component l with cost c lp < c lc (that is, the degraded component l is replaced by a new one).

• A perfect corrective repair (PCR) is done for component l, l ∈ M k , if the system fails at ((k -1)τ, kτ ] and

T l > (k -1)τ & T l < kτ.
Define P cl (kτ ) as probability of perfect corrective repair for component l; l = 1, 2, ..., J, at kth inspection time.

• An opportunistic perfect preventive repair (OPPR) is done for component l, l ∈ M k , if the system fails at ((k -1)τ, kτ ] and

T l > kτ & p lM k > ρ.
Define P pl (kτ ) as the probability of the opportunistic perfect preventive repair for component l; l = 1, 2, ..., J, at kth inspection time, The time from the component installation to its first replacement or the time between two successive replacement of each component is referred to as a renewal cycle. Let L and L j denote the average long-run maintenance cost per unit of time for the system and component j, respectively. Therefore, based on the renewal reward theorem, [START_REF] Ross | Introduction to probability models[END_REF], the expected long-run maintenance cost rate for component j is

L j (τ, ρ) = lim t→∞ C j (t) t = E(C rj ) E(T rj )
where E(C rj ) and E(T rj ) are total expected cost during a replacement cycle and expected length of the replacement cycle for component j, respectively such that

E(C rj ) = ∞ k=1 ( kc ins J +c jp )P pj (kτ )+( kc ins J +c jc )P cj (kτ )
and

E(T rj ) = ∞ k=1
kτ P pj (kτ ) + P cj (kτ ) .

Finally, the total expected long-run maintenance cost rate for the series system until time t is given by [START_REF] Zamani | Proper efficiency, scalar-ization and transformation in multi-objective optimization: unified approaches[END_REF], [START_REF] Ehrgott | Multicriteria optimization[END_REF] L(τ, ρ)

= J j=1 L j (τ, ρ).
The aim is to propose an inspection interval τ which minimises L(τ, ρ):

τ * = argminL(τ, ρ) = J j=1 L j (τ, ρ)
This joint optimisation of costs permits to take into account all the components reliability in the optimisation procedure.

IV. NUMERICAL EXAMPLE A series system with J = 3 components is considered. The system is inspected periodically at times kτ ; k = 1, 2, ..., with cost c ins for the system and M k ; k = 1, 2, ... are corresponding masked sets, thus the collected data are (kτ, M k ); k = 1, 2, .... When the system fails, an opportunistic perfect preventive replacement is carried out for components that are in the masked set with p l > ρ and a perfect corrective replacement is made for failed component with costs c lp and c lc , respectively. Two illustrative examples have been constructed to clarify previous sections. In both of them, it is assumed that the lifetime of components are independent and follow the Weibull distribution with parameter set (α j , β j ); j = 1, 2, 3, as follows:

f Tj (t j ) = α j β j ( t j β j ) (αj -1) exp(-( t j β j ) αj ).
We set parameters as in Table I TABLE I: The optimal inter-inspection interval value (τ ) as decision parameter with different cost rates and ρ

unit costs c ins 1 c 1p 1.5 c 2p 2.5 c 3p 3.5 c 1c 2 c 2c 4 c 3c 5
All masked sets are considered as completely masked sets which is similar to missing setup, that is, M k = {1, 2, ..., J}; k = 1, 2, ..., thus

p jM k =P (K = j M k , u ∈ ((k -1)τ, kτ ]) = ∫ kτ (k-1)τ f T,K (u, j)du ∫ kτ (k-1)τ ∑ l ′ ∈M k f T,K (u, l ′ )du since ∀j ∈ {1, 2, ..., J} and t ∈ [0, ∞) P (M k = M K = j, t) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 if M = {1, 2, ..., J} 0 if M ≠ {1, 2, ..., J}
where the inter-inspection interval, τ , is considered as decision parameter and should be optimized through maintenance optimization problem. It can be noticed in Tables II and III that the cost rate is very sensitive to the failure rate distribution. Indeed with the same maintenance unit costs, the total cost rate can be very high for an increasing failure rate. When the failure rate is increasing a with the same masking probability threshold the impact of a missing failure can be very high between two inspections.

(k -1)τ < T < kτ l ∈ M i Set: l = 1 (k -1)τ < T l < kτ L(τ, ρ) Algorithm 1: Imperfect Corrective Maintenance Policy Yes T l > kτ & p l > ρ T l > kτ & p l ≤ ρ OPPM PCM c lp c lc L l (τ, ρ) l = l + 1 l ≤ J l > J No
Fig. 1: The optimal value of ρ as decision parameter considering (α 1 , α 2 , α 3 ) = (1.25, 2, 2.5) and (β 1 , β 2 , β 3 ) = (2.5, 2.25, 2.75) for the optimal inspection interval τ = 0.67.

In Figure 1, 2 and 3, the cost rate is depicted in function of component cause probability threshold ρ considering the optimal inspection. In Figure 1, since all the failure rates are non-decreasing the cost rate is convex for decision parameters and we can easily identify an optimal value. The optimal inspection interval is relatively small which requires frequent inspections and the cause probability threshold is close to 50%. Since the inspections are frequent the risk of missing the cause is not very high.

In Figure 2 and3, since all the failure rates are decreasing the minimal cost rate exists but the cost rate variations around the optimal value are not very high. Indeed in this framework the system is getting better even though the failures occurs, by relatively frequent inspections and a low cause probability the maintenance operations are occurs frequently. However, since failures are less frequent the maintenance decision rule is more often postponed.

As the failure rate is decreasing (the failure rates in Figure 3 decrease faster than in Figure 2) the flatness around the optimal value is more obvious. The inspections are less frequent and the cause probability threshold still very low.

V. CONCLUSIONS

In this paper an opportunistic maintenance is proposed for a multi-component system with masked data. The maintenance policy is described and a cost criterion optimisation is proposed. A numerical example illustrate the methoddology.

Fig. 2 :

 2 Fig. 2: The optimal value of ρ as decision parameter considering (α 1 , α 2 , α 3 ) = (0.5, 0.25, 0.7) and (β 1 , β 2 , β 3 ) = (1.5, 1.25, 1.75) for the optimal inspection interval τ = 0.88.

Fig. 3 :

 3 Fig. 3: The optimal value of ρ as decision parameter considering α = (0.5, 0.25, 0.7) and (β 1 , β 2 , β 3 ) = (0.25, 0.5, 1.5) under different failure rates, for the optimal inspection interval τ = 1.11.

TABLE II :

 II The optimal inter-inspection interval value (τ ) as decision parameter with different cost rates and ρ considering a constant failure rate (α j = 1 for j = 1, 2, 3)

	τ *	ρ	Cost Rate
	0.03	0.00	154.38
	0.04	0.00	125.23
	0.05	0.10	102.16
	0.06	0.10	92.90
	0.07	0.10	86.11
	0.09	0.10	69.64
	0.11	0.10	62.42
	0.14	0.10	55.72
	0.17	0.10	48.73
	0.26	0.10	36.90
	0.27	0.10	34.14
	0.31	0.10	32.96
	0.42	0.10	24.21
	0.44	0.10	23.16
	0.50	0.10	21.13
	0.53	0.10	20.10
	0.54	0.10	19.90

TABLE III :

 III The inter-inspection interval value (τ ) as decision parameter with different cost rates and ρ considering an increasing failure rate

	τ	ρ	Cost Rate
	0.01	0.00	368.44
	0.01	0.20	445.45
	0.02	0.20	234.48
	0.03	0.10	161.29
	0.04	0.10	123.22
	0.10	0.20	56.15
	0.30	0.20	24.18
	0.40	0.20	19.50
	0.50	0.20	16.30
	0.70	0.20	12.72
	0.80	0.20	11.42
	0.82	0.20	11.20
	0.85	0.20	10.88
	0.86	0.20	10.77
	0.87	0.20	10.68

ACKNOWLEDGMENT This work received support from the French government under the France 2030 investment plan, as part of the Initiative d'Excellence d'Aix-Marseille Université -A*MIDEX (AMX-19-IET-010)


Yue Sun
email: yue.sun@univ-lille.fr

Emmanuel Roubin
email: emmanuel.roubin@univ-grenoble-alpes.fr

Jean-Baptiste Colliat
email: jean-baptiste.colliat@univ-lille.fr

Jianfu Shao
email: jian-fu.shao@polytech-lille.fr

Damage and fracture in brittle materials with enriched finite element method: Numerical study

. As a Finite Element based model, this model performs two kinds of enhancements: i) strong discontinuities, which allows illustration of cracks and fractures; ii) weak discontinuities, which manages to represent heterogeneities explicitly without any need of mesh adaptation. Many existing EFEM models have shown their ability to simulate a lot of main features of brittle/quasi-brittle materials at the macroscopic scale, such as the asymmetric responses in traction/compression. With a large number of finite elements and explicit heterogeneous structures, the mechanical behaviors that applied at the local scale could be very simple. In this chapter, we are looking forward to describing more of the main characters of such material by adding a closure mechanism to the model.

Introduction

The growth and propagation of micro-cracks is a leading cause of fatigue and damage for brittle/quasi-brittle materials, such as shale, concrete-like materials, and rocks. The computational failure mechanics is nowadays a study of significant value. One of the characteristics of the materials referred to above is that they are materials with complex structures. For example, beginning at the mesoscopic scale, concrete exhibits aggregate pieces and macro-pores, and complexity of its structure increases when increasingly finer scales are considered. Another characteristic of these kinds of materials is that they show complex mechanical behaviors, such as the non-symmetric responses in traction and compression as well as the hysteresis phenomenon in cyclic loadings, see Fig. 1. From a general point of view, we can make an assumption that those two complexities are strongly linked. Fig. 1: Hypothesis that complex mechanical behaviors of brittle/quasi-brittle materials [START_REF] Mazars | The unilateral behaviour of damaged concrete[END_REF][START_REF] Terrien | Emission acoustique et" comportement mecanique post-critique" d'un beton sollicite en traction[END_REF][START_REF] Hans | Fracture mechanics of an elastic softening material like concrete[END_REF] can emerge from the upscaling of elements with simple behaviors and explicit heterogeneous structure [START_REF] Stamati | Phase segmentation of concrete x-ray tomographic images at meso-scale: Validation with neutron tomography[END_REF].

The hypothesis stems from the non-linear complex adaptive system [START_REF] Ahmed | An overview of complex adaptive systems[END_REF]. The essence of the complex system is that each individual constituent conforms to (very) simple rules. However, when the system is studied as a whole, emergent responses can be observed after upscaling which are not exist at the single element. Besides, by understanding the behaviors of each element, it may be beneficial for us to understand the entire system.

Naturally, the macroscopic mechanical behaviors of concrete are known a priori from experimental tests, see Fig. 1. The complex heterogeneous structure of concrete is thus defined as the explicit heterogeneous geometry of the material. At the mesoscopic scale, the heterogeneities are referred to as the aggregates and macro-pores, see Fig. 1. It can be constructed within the framework of the finite element method, with a large number of elements and explicit representation of heterogeneities. Moreover, a newly developed technique, namely X-ray tomography, makes it possible to establish a morphological structure of concrete based on real tomographic images [START_REF] Stamati | Phase segmentation of concrete x-ray tomographic images at meso-scale: Validation with neutron tomography[END_REF]. Therefore, it is interesting to look backward to clarify the simple rules of each element.

Based on the previous assumption, attempts to solve the simple mechanical behaviors of each element have been made by former studies. The used method in this chapter is based on the strong discontinuity approach, namely the Enriched Finite Element Method (EFEM) [START_REF] Ortiz | A finite element method for localized failure analysis[END_REF][START_REF] Juan Carlos Simo | An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids[END_REF]. EFEM is an element-based enhancement model, with additional degree of freedom attached inside the element as internal variables. It is capable of yielding the mesh dependence without adding special artificial numerical parameters. As the energy is dissipated over a two-dimension discontinuity interface, by applying a specific kinematic enhancement to the element, the total dissipated energy becomes independent of the mesh size [START_REF] Juan Carlos Simo | An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids[END_REF]. Besides, as the additional interval variable is an element-based enhancement, the increasing number of cracked elements will not affect the size of the assembled stiffness matrix for the total system. Therefore, the model is capable of simulating a system with a large number of elements that carry strong discontinuities. The ratio of the fractured elements to the total elements can be very high. Another benefit of EFEM is that the strong discontinuities can pass through the element in an arbitrary trajectory. Therefore, it is capable of performing simulation of complex and accurate cracks, such as cracks with several branches, cracks stopped by heterogeneities, and a crack that meets another crack.

Moreover, individual to strong discontinuity, the weak discontinuity [START_REF] Ortiz | A finite element method for localized failure analysis[END_REF] can be alternatively included in the model to represent heterogeneities, even for complex geometry [START_REF] Roubin | Modélisation EF et morphologique de milieu hétérogènes à léchelle mésoscopique: applications aux matériaux à matrices cimentaire[END_REF]. Due to the non-adapted meshing method, the heterogeneities can be taken into consideration in an explicit way [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF].

The study performed in this chapter is closely related to the previous work undertaken by [START_REF] Roubin | Modélisation EF et morphologique de milieu hétérogènes à léchelle mésoscopique: applications aux matériaux à matrices cimentaire[END_REF], in which each element is governed by two mechanisms: localization and traction-separation. It has been demonstrated that even though only "opening" mechanism is applied to elements, the model is capable of performing failure behaviors of concrete in traction as well as in compression. Moreover, non-symmetric behaviors can also be observed. However, to the author's knowledge, the behavior performed by crack closure in EFEM has yet to be taken into consideration. Therefore, based on the prior studies, the closure mechanism of the cracks is also taken into account in this chapter.

Kinematics of discontinuities in solids

In this section, we summarize the basic notations that we used in this chapter. The referenced domain Ω ⊂ R n dim is a solid exhibiting heterogeneities and cracks, see Fig. 2. Its smooth boundary Γ ⊂ R n dim -1 can be divided into two disjointed boundaries: the displacement boundary Γ u and the traction boundary Γ t . Within the framework of Finite Element Method, this specimen can be discretized by means of standard isoparametric elements Ω = ∪ n e e=1 Ω e . Dealing with the explicit heterogeneities, the mesh method is non-adapted. Since the heterogeneity surfaces are defined a priori and independent to any mechanical calculation, there will be a set of elements crossed by the interface of heterogeneity S ε , dividing the elements into two parts Ω + e and Ω - e with a unit vector n pointing from Ω - e to Ω + e . At the inner of the dealing body, another set of elements is divided into two parts by the failure path S u . Naturally, the discontinuity surface S ε is defined by the geometry of material, whereas the path and orientation of the strong discontinuity S u usually depend on particular criteria. Fig. 2: Illustration of the heterogeneities and cracks in the referent body and the corresponding enhanced elements.

Dealing with the weak discontinuities and the strong discontinuities, three cases can be present in an element: i) the element contains only a weak discontinuity, ii) the element contains only a strong discontinuity, iii) both of them are present in an element. In the third case, it is considered that the kinematics failure will take place in the heterogeneity interface. It implies that if the element is "close" enough to the heterogeneity, the crack opening will be localized on the heterogeneity surface.

Furthermore, it is assumed that the bulk part has a purely elastic kinematics relationship, and the kinematic behaviors of the strong/weak discontinuity are independent. Hence, the strain field admits an additive form

ε(x) := ∇ sym u(x) + ε(x) + ε(x), (1)
where

∇ sym := 1 2 [∇(•) + ∇ T (•)
] denotes the symmetric gradient operator. Notation • / • represents the weak / strong discontinuity.

Kinematics description of weak discontinuity

In this section, interest is focused only on the elements which exhibit weak discontinuities. The elements of this kind contain different elastic parameters: Young's Modules and Poisson's ratios, thus a strain discontinuity is emerged from this difference, as it shows in Fig. 3. With the normal vector n pointing from one sub-domain Fig. 3: Representation of the strain discontinuity in a tetrahedral element. [START_REF] Vallade | Modélisation multi-échelles des shales: influence de la microstructure sur les propriétés macroscopiques et le processus de fracturation[END_REF] to another, we can construct an orthonormal basis (n, m,t). The continuity of the displacement field should be respected crossing the discontinuity surface, so we have:

u + (x) = u -(x), ∀x ∈ S ε .
(

) 2
The strain jump of the weak discontinuity [|ε|] := ε + (x)ε -(x) can be calculated by the symmetrical gradient of the corresponding displacement field

∇ sym u| S ε = 1 2 (∇u + ∇ T u)| S ε with ∇u =   ∂ u 1,n ∂ u 1,m ∂ u 1,t ∂ u 2,n ∂ u 2,m ∂ u 2,t ∂ u 3,n ∂ u 3,m ∂ u 3,t   (n,m,t) . (3)
It can be deduced that the jump of the strain field takes place in the direction perpendicular to the discontinuity surface but not in the direction along the interface. Accordingly, ∇(ũ

+ ) • m = ∇(ũ -) • m and ∇(ũ + ) •t = ∇(ũ -) •t.
Take this into Eq. ( 3), it can be seen that the last two colons of the matrix ∇u have the same value for both sub-domains. Hence, [ε] n [ε] m [ε] m are three constants that can entirely define the jump of strain field [START_REF] Roubin | Multiscale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]. We can take the displacement field ũ in a first order form [START_REF] Ibrahimbegovic | Multi-scale modelling of heterogeneous structures with inelastic constitutive behavior[END_REF]:

ũ(x) = Θ n • (x -ξ )([ε] n n + [ε] m m + [ε] t t) with Θ = Θ + ∀x ∈ Ω + e Θ -∀x ∈ Ω - e , (4)
with ξ denotes the position of interface surface S ε and n • (xξ ) can be seen as a signed distance to the discontinuity surface. Θ is a still undefined function of Ω e , the explicit expression will be given latter with the help of variational formulations.

We can obtain the weak shape of the enhanced strain field by taking the symmetrical gradient of ũ [START_REF] Benkemoun | Poroelastic two-phase material modeling: theoretical formulation and embedded finite element method implementation[END_REF]:

ε = ∇ sym (ũ) = Θ ([ε] n n ⊗ n + [ε] m 2 (n ⊗ m) sym + [ε] t 2 (n ⊗ t) sym ) (5)

Kinematics description of strong discontinuity

In this part, the main interest is focused on the strong discontinuity kinematics [START_REF] Juan | A class of mixed assumed strain methods and the method of incompatible modes[END_REF][START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. part 1: Fundamentals[END_REF][START_REF] Nathan | Discontinuous modelling of strain localisation and failure[END_REF]. Let us consider an element Ω e exhibits a strong discontinuity dividing itself into two sub-domains, the discontinuity surface is noted as S u . The discontinuity of the displacement field can be decomposed into the following function [START_REF] Roubin | Multiscale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]

u = ū + (H S u -ϕ e )[|u|], (6)
with H S u is the Heaviside function centered on S u , ū is the regular part of the displacement field to impose the standard boundary conditions [START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. part 2: Numerical simulation[END_REF], and [|u|] is a continuous function representing the displacement jump. For the sake of simplicity, [|u|] is considered as a unity function over the finite element. By taking the symmetric gradient of the displacement field, the corresponding strain field can be obtained as:

ε = ∇ sym u = ∇ sym û regular + (H S u -ϕ e )∇ sym ([|u|]) -([|u|] ⊗ ∇ϕ e ) sym bounded enhancement + δ Γ d ([|u|] ⊗ n) unbounded enhancement . (7)
The strain field admits an additive form of three parts. The regular displacement field respects the standard stress-strain constitutive equation. The bounded enhancement can be simplified by taking the jump of the displacement field [|u|] as a constant function. It results into a null value of (H Γ dϕ e )∇ sym ([|u|]). Then we can see that the gradient of the Heaviside function brings an unbounded enhancement part, which bears a Dirac-delta distribution δ . Centered at the discontinuity interface, it carries infinite value at the interface and null value otherwise. Thus far, Eq. ( 7) can be written as:

ε = εb + εu = -([|u|] ⊗ ∇ϕ e ) sym + δ S ([|u|] ⊗ n) (8)
Regarding the unbounded part in the function, some properties such as the traction continuity condition seem difficult to conform. Nevertheless, some solutions are available within the FE context. In this study, the used solution is called Discrete Strong Discontinuity Approach(DSDA) [START_REF] Oliver | On the discrete constitutive models induced by strong discontinuity kinematics and continuum constitutive equations[END_REF], which leads to an underlying discrete model at the discontinuity surface. And a traction vector T which links with the crack opening [|u|] can also be introduced

σ + (x) • n = σ -(x) • n = T , ∀x ∈ S. (9)
It can be then used to model the failure mechanism at the local scale. And beyond the discontinuity, the bulk Ω \S remains an elastic behavior, which can be described by Hooke's law σ = C : ε.

In the case of exhibiting a weak discontinuity, the two sub-domains are isotropic, the operator C +|-is calculated separately by Young's modulus E +|-and Poisson ration ν +|-.

Finite Element implementation

Incompatible modes

Based on existing EFEM studies, a three-field variational formulation is used here [START_REF] Washizu | Variational methods in elasticity and plasticity[END_REF][START_REF] Roubin | Multiscale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF][START_REF] Vallade | Modélisation multi-échelles des shales: influence de la microstructure sur les propriétés macroscopiques et le processus de fracturation[END_REF]:

HW û( û, ε, σ ; η) = Ω ∇ sym η : σ dΩ - Ω η • ρ bdΩ - Γ η • td∂ Ω = 0, (11a) HW σ (û, ε, σ ; τ) = Ω τ : (∇ sym û -ε) = 0, (11b)
HW ε (û, ε, σ ; γ) = Ω γ : ( σ (ε) -σ )dΩ = 0, (11c)
where (û, ε, σ ) represents the standard displacement field, the strain field, and the stress field. Then three corresponding virtual field are noted as ( η, τ, γ), also mutually independent. Then a classic method of incompatible modes is introduced here, mentioned as the Assumed Strain Method. The central idea is assuming that both actual and virtual fields are enhanced in an additional way, making them divided into three parts, the standard part, the enhanced weak discontinuity, and the enhanced strong discontinuity:

ε = ∇ sym û compatible + ε weak + ε strong incompatible , and γ = ∇ sym η compatible + γ weak + γ strong incompatible . (12)
Considering the orthogonal relationship between the enhanced and the virtual part of the strain and stress field [START_REF] Roubin | Multiscale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF], we can derive the following simplified formulation:

Ω ∇ sym η : σ (∇ sym û + ε + ε)dΩ - Ω η • ρ bdΩ - Γ t η • td∂ Ω = 0, (13a) Ω γ : σ (∇ sym û + ε + ε)dΩ = 0, (13b)
Ω γ : σ (∇ sym û + ε + ε)dΩ = 0. ( 13c
)
To ensure the convergence of the method, the patch test should be respected after the orthogonal condition is imposed. By assuming that the discontinuity interface is flat, which makes the normal vector of the interface n a constant over each element, we will take the Θ as a function that depends only on the volume of sub-domains [START_REF] Roubin | Modélisation EF et morphologique de milieu hétérogènes à léchelle mésoscopique: applications aux matériaux à matrices cimentaire[END_REF]:

Θ + = V - V and Θ -= - V + V . (14)
Subsequently, the strong enhanced strain field consists of two parts: the unbounded part δ Γ d ([|u|] ⊗ n) s and the bounded part γb . The later can be explicitly defined by taking into consideration of the zero mean condition. Following the Enhanced Assumed Strain(EAS) method, the complete form of the strong enhancement of the strain field can be obtained by

γ = δ Γ d ([|η|] ⊗ n) sym unbounded - A V ([|η|] ⊗ n) sym bounded . (15)

Finite Element discretization

It is present in this part the resolution of the three-field variational formulation Eq. ( 13). The discretization of the strain field is different depending on whether it is standard or virtual. The standard strain field is enhanced using KES in a kinematics point of view, whereas the virtual strain field is enhanced using the EAS (refers to as • * in the following part) in a statical point of view. Based on the DSDA, the standard discretization of the strain field contains only the bounded enhanced part, while the virtual discretization has both the bounded and unbounded enhancement. The discretization is written as:

ε = ∇ sym û + ε + ε = B d + G w [|ε|] + G s [|u|], (16a)
γ = ∇ sym η + γ + γ = B d + G w [|γ|] + (G * s,b + G * s,u ) [|η|]. (16b)
Several mentioned notations above are: B(= ∂ N) the standard strain interpolation matrix, G w the standard and virtual field corresponding to the weak discontinuity, G s the bounded part of the standard field corresponding to the strong discontinuity, It is worth noting that the strain fields are all expressed in Voigt notation.

G * s := G * s,b + G * s,
Considering the discretization of the system, we can write the global equilibrium equation in residual form

R := n e A e=1 Re = f ext - n e A e=1 Ω e B T σ dΩ = 0, (17)
where

n e A e=1
represents the standard assembly operator which groups all finite elements n e of the global system, and f ext denotes the external force, which writes as:

f ext = n e A e=1 Ω e N T b dΩ + Γ e t N T t dΓ . ( 18
)
We can obtain the equations by injecting parameters of Eq. ( 16) into Eq. ( 13):

Ω G T w σ (d, [|ε|], [|u|]) dΩ = 0, (19a)
Ω (G * ,T s,b + G * ,T s,u ) σ (d, [|ε|], [|u|])dΩ = 0. ( 19b
)
Then Eq. ( 15) can be changed into form:

γ = (G * s,b + G * s,u ) [|η|] = - A V ([|η|] ⊗ n) sym + δ S ([|η|] ⊗ n) sym , (20)
where (• × n) sym is noted as a operator H * s for the rest part, leads to an equivalent Voigt notation for the traction vector

T = σ • n matrix notation = H * ,T s σ Voigt notation . ( 21
)
By considering one of the proprieties of the Dirac-delta distribution δ S , f (x)δ S dx = f (x = S)dx, the unbounded part of the virtual strain field can be further simplified into the form

Ω G * ,T s,u σ (d, [|ε|], [|u|])dΩ = S T (d, [|ε|], [|u|])d∂ Ω . (22)
As for the interpolation matrix corresponding to weak discontinuity, G w is defined separately in two parts Θ +|-as:

G w = G + w = Θ + H w = V - V H w in Ω + e , G - w = Θ -H w = -V + V H w in Ω - e . (23)
Recording to the equation Eq. ( 8), we can see that the bounded part of the strong enhancement of the standard strain field writes as ([|u|] ⊗ ∇ϕ e ) s , which corresponds to the term G s [|u|] in Eq. 16a. Hence G s can be obtained as an equivalent symmetric operator, (• ⊗ ∇ϕ e ) s , with ϕ e is an explicitly defined arbitrary function to separate nodes at Ω + e from Ω - e . The arbitrary function is defined as follows:

ϕ e (x) = n e ∑ a=1 N a p a with p a = 1 if node number a ∈ Ω + e , 0 if node number a ∈ Ω - e . (24)
In summary, the discretized system gives as [START_REF] Roubin | Multiscale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF] R :=

n e A e=1 f int -f ext = n e A e=1 f e ext - Ω e B T σ (d, [|ε|], [|u|])dΩ = 0, (25a)
Re [|ε|] = Ω e G T w σ (d, [|ε|], [|u|]) dΩ = 0 , (25b)
Re [|u|] = Ω e G * ,T s σ (d, [|ε|], [|u|]) dΩ = Ω e \S G * ,T s,b σ (d, [|ε|], [|u|])dΩ + S T d∂ Ω = 0. ( 25c
)
In the above equations, the first one Eq. (25a) represents the global equilibrium equation of a standard Finite Element system, whereas the latter two equations Eq. (25b) and Eq. (25c) are local equations for each element, and carry the enhancement parts. The enhancements [|ε|] and [|u|] are considered as the internal variables and are both solved at the element level. Hence, the global equilibrium equation Eq. (25a) always has the same size no matter how many heterogeneities exist in the system or how many elements begin to crack.

Admissible discrete model with closure mechanism on the discontinuity surface

As it is mentioned in the previous part, a traction vector T is continuous over the element and links the two sub-domains crossing the discontinuity surface S. Our discrete model is formulated based on the description of the relationship between the traction vector and the crack-opening [|u|].

In existing EFEM studies, many discrete models have been proposed and well documented. For instance, the traction-opening fracture [START_REF] Roubin | Multiscale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF], the sliding-opening shear fracture [START_REF] Hauseux | Propagation dincertitudes paramétriques dans les modèles numériques en mécanique non linéaire: applications à des problèmes dexcavation[END_REF], and hydraulic mechanical coupling problem [START_REF] Vallade | Modélisation multi-échelles des shales: influence de la microstructure sur les propriétés macroscopiques et le processus de fracturation[END_REF]. Based on these previous works, the attempt in this section is to propose a discrete model that considers the closure of cracks.

Localization criterion

At the beginning of loading, the element has an elastic behavior, until the stress field, or the strain field reaches a specific value, i.e., the localization criterion. Our localization criterion is stress-based. The yield stress, note as σ y , is considered as a local parameter of the material. The localization criterion wrote as:

Φ l = σ eq -σ y . (26)
In this equation, a negative value of Φ l represents an elastic behavior of the element while a positive value leads to the localization and appearance of the crack. As we can see, the equivalent stress σ eq is a constant value. Two cases should be considered: i) the studied element is isotropic without weak discontinuity, ii) the weak discontinuity is present in the element.

In the first case, the equivalent stress is determined by the major principal stress, σ eq = σ I . The orientation of the discontinuity interface n is the corresponding eigenvector, n I . The localization criterion writes:

Φ strong l = σ I -σ y . (27)
In the second case, the element carries different phases, it is assumed that the crack opens at the heterogeneity surface between the two sub-domains. The normal vector n is defined a priori by the geometric information of the material and thus is independent of the stress state. The equivalent stress is calculated by the projection of the traction vector T on the normal vector n, the localization criterion writes:

Φ weak l = n • T -σ y ( 28
)
We present here an isotropic one-dimensional example. At the state of the stress field reaches the yield stress, the localization takes place. The displacement field and local constitutive behavior are shown in Fig. 4.

Failure criterion -Traction separation law

Following the localization of the element, the failure behavior is driven by the traction-separation law [START_REF] Roubin | Multiscale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

Φ o = σ eq -( σ y -q o ([u]) ), with [|u|] = [u] • n. (29)
A positive value of the traction-separation criterion Φ o means an in-equilibrium state of the element and the crack needs to go further until it reaches an equilibrium state. A negative value presents an elastic loading or unloading, or closure of cracks, which we will discuss later. In this equation, q o ([u]) is referred to as the hardening function. It is a continuous simple decreasing function in terms of the opening value, with its magnitude value is located at [u] = 0. It will approach to zero as the [u] = 0 increases. The hardening function is defined as:

q o ([u]) = σ y 1 -exp - σ y G op [u] , (30)
where σ y and G op are both local material parameters. The latter is referred to as the fracture energy, which governs the amount of necessary energy to create a fully opened fracture. Physically, a bigger value of the fracture energy states for a tougher and more durable material. Assuming a the same one-dimension example as shown in Fig. 4, the tractionseparation phase is depicted in Fig. 5. As the imposed displacement ∆ u increases, the crack will propagate, and lead to decreasing value of the equivalent stress. Since the traction continuity condition is always satisfied, the stress in the bulk volumes Ω e \S will also decrease as shown in Fig. 5b.

Closure criterion

Following the traction-separation law, a closure law is considered in our discrete model if the normal stress is turned into compressive stress. As for the model exhibits no closure mechanism, see Fig. 6a, the crack opening value will not decrease even under high compressive stress. Unlike most phenomenological models that mainly describe stress-strain curves, our model is built at the mesoscopic scale, the complexity consists in depicting the relationship between the traction vector on the discrete discontinuity surface and the crack opening value. Upon physical considerations, the applied closure mechanism is replied on three hypotheses.

σ y [u] σ eq G op [u] max (a) Phenomenological model without closure mechanism σ y [u] σ eq G op G cl [u] max (b) Phenomenological model with exponential closing law σ y σ y [u]
σ eq (c) Phenomenological model with damaged tensile strength Fig. 6: Hypothesis for local constitutive models with non-linear closing laws and damaged tensile strength at the discontinuity interface.

• Hypothesis 1. The closing criterion is driven by an exponential and continuous function, [u] will approach towards zero if the equivalent stress turns to infinite compression but will never close completely, see Fig. 6b. This hypothesis is proposed upon a physical consideration, that an already created crack of a brittle/quasi-brittle material should not disappear even is applied high compressive stress. Furthermore, this driven function also verifies sort of the spirit of "symmetric", that in the opening phase, the material will become progressively fragile, and in the closing stage, the closing of crack will become more and more difficult.

• Hypothesis 2. The amount of required energy for a complete closing G op equals to the energy dissipated at the opening phase G cl .

The design of this hypothesis is an attempt to make the model remain simple and clean. No additional parameters are needed to depict the closure mechanism.

• Hypothesis 3. As the element reopens after the closing phase, the calculation of the new traction-separation law will base on the residual tensile strength σ y and the residual fracture energy G re , see Fig. 6c. This third hypothesis is proposed to deals with the already damaged elements.

Physically, the damaged element should be weaker and more fragile than the sound elements. Even though the crack can be partially closed during the closing phase, the yield stress of the element can not be restored to the value of the sound element. And a certain amount of fracture energy has already been dissipated in the opening phase. Therefore, the calculation of the reopening criterion is defined based on the residual tensile strength and fracture energy. As a result, damages that have occurred to an element throughout history are irreversible.

Unloading procedure

The first procedure that occurs to an element subsequently to the opening phase is the unloading procedure. This procedure takes place on the bulk volumes outside the discontinuity surface, which means that the crack opening value [u] will not change.

After the opening stage, an amount of tensile energy is stored in the bulk volumes.

The discontinuity surface will always under tensile stress until the elastic energy is totally released, i.e., the strain field is null on the bulk volumes. Therefore, the unloading procedure is a pure elastic phase. We can see from Fig. 7 that it happens when ∆ u = [u].

Closure of cracks

After the unloading process, the closure of the crack occurs when the stress on the discontinuous surface switches to compressive stress. The decreasing of the crack opening value is driven by the closing criterion. As it is shown in Fig. 8, the closing criterion is non-linear and continuous. It is assumed that the necessary needed energy for a complete close G cl equals to the amount of energy that dissipated at the opening stage G op . The closing procedure will become more and more difficult and as a result, the crack can never completely closed. It is worth some particular attention at the point of ∆ u = 0, where the imposed displacement equals to zero, whereas the crack opening value is still positive, and the element is under compressive stress.

Reloading procedure

Let us now apply a reloading displacement to the element. Firstly, the element will release the stored compressive energy in bulk volumes. Then the stress on the discontinuity surface will turn to positive as the reloading continues until the equivalent stress reaches the residual critical tensile strength σ y . And the fracture energy is also taken as the residual value G re ,see Fig. 9, which leads to a weaker and more fragile material than the sound one.

Governing equations

In the previous parts, the kinematics behaviors of the model on the discontinuity surface have been decomposed into several different phases. In the closing process, it is assumed that the necessary energy for a total crack closure equals to the dissipated energy in the opening phase. This amount of energy is calculated as:

E op = [u] max 0 σ eq ([u]) d[u] = G op (1 -exp(- σ y [u] max G op )) = G cl . (31)
As a result, the driven function of the closing phase can be completely defined by the material's parameters G op , σ y , and the maxim opening value [u] max , no additional parameters are required. The closing energy should respect the following equation:

E cl = 0 [u] max σ eq ([u]) d[u] = E op = G cl . (32)
The choice made here is :

Φ c = -σ eq + G cl [u] max ln [u] [u] max q c , (33)
where q c is referred to as the hardening function during the closing procedure. A positive value of Φ c means that the crack should continue closing to reach the equilibrium state, and a negative value of Φ c means that the element is under elastic unloading or reloading. Since the maximum crack opening value [u] max is always greater than or equal to the actual crack opening [u], it is impossible to have Φ o and Φ c both positive at the same time.

In the reopening phase, the critical parameters of the traction-separation law are taken as the residual tensile strength (marked as σ y in Fig. 9c) and the residual fracture energy (marked as G re in Fig. 8c and Fig. 9c). The updated traction-separation law writes as

Φ o = σ eq -σ y 1 -exp - σ y G re [u] , with    σ y = σ y exp - σ y ([u] max -[u] cl ) G op G re = G op exp - σ y ([u] max -[u] cl ) G op . ( 34
)

Numerical resolution of the discrete Finite Element system

In this section, the interest is focused on the numerical resolution of the discrete system, i.e., Eq. ( 25). The former equation presents a standard global system, while the latter two equations are formulated at the local scale, and carry non-linear kinematics behaviors. Hence, the system to be solved is non-linear. The linearisation of equations is present in the following part. Then the resolution of the global and local system is also introduced in this section.

Linearisation of equations

The integration of the global system (Eq (25a)) and the weak enhancement (Eq. ( 25b)) are given as the following equations [START_REF] Roubin | Modélisation EF et morphologique de milieu hétérogènes à léchelle mésoscopique: applications aux matériaux à matrices cimentaire[END_REF]:

f e int = Ω e B T σ (d, [|ε|], [|u|])dΩ = K bb d + K bw [|ε|] + K bs [u], (35)
Re [|ε|] = Ω e G T w σ (d, [|ε|], [|u|]) dΩ = K wb d + K ww [|ε|] + K ws [u], (36)
with

K bb = B T V + C + +V -C -B K bw = V + V - V B T (C + -C -)H w K bs = B T (V + C + +V -C -)G s n p K wb = V + V - V H w T (C + -C -)B K ww = V + V - V H w T (V -C + +V + C -)H w K ws = V + V - V H w (C + -C -)G s n p .
Based on the Newton's method, their linearisation formulation gives as [START_REF] Roubin | Multiscale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]

n e A e=1 K bb ∆ d (k+1) n+1 + K bw ∆ [|ε|] (k+1) n+1 + K bs ∆ [u] (k+1) n+1 = - n e A e=1 f e int (k) n+1 -f e ext , K wb ∆ d (k+1) n+1 + K ww ∆ [|ε|] (k+1) n+1 + K ws ∆ [u] (k+1) n+1 = -Re [|ε|] (k) n+1 .
The linearisation of the strong enhancement should be discussed in two cases depending on the element is in the opening stage or the closing stage. Assuming that the discontinuity surface is flat and the stress over the element is constant, the traction vector T can be calculated based on the average value of weighted volumes of sub-domain V +|-:

T = 1 V H * ,T s (V + σ + +V -σ -). ( 37
)
The linearisation of the strong enhancement of the traction-opening phase writes as [START_REF] Roubin | Modélisation EF et morphologique de milieu hétérogènes à léchelle mésoscopique: applications aux matériaux à matrices cimentaire[END_REF]:

K s * b (k) n+1 ∆ d (k+1) n+1 + K s * w (k) n+1 ∆ [|ε|] (k+1) n+1 + (K s * s + K qo ) (k) n+1 ∆ [u] (k+1) n+1 = -Φ o (k) n+1 (38) with K s * b = ∂ σ eq ∂ T 1 V H * ,T s (V + C + +V -C -)B K s * w = ∂ σ eq ∂ T V + V - V H * ,T s (C + -C -)H w K s * s = ∂ σ eq ∂ T 1 V H * ,T s (V + C + +V -C -)G s n p K qo = σ 2 y G op e -σ y [u]/G op
As for the element in the closing stage, the closure criterion Φ c = -σ eq + q c ([u]) have to be taken into consideration. It is a function in terms of two variables: {T , [u]}, thus the increment of Φ c gives as:

∆ Φ c = - ∂ σ eq ∂ T ∆ T + ∂ q c ∂ [u] ∆ [u], (39)
= -K s * b ∆ d -K s * w ∆ [|ε|] -K s * s ∆ [u] + K qc ∆ [u], with K qc = G cl [u] max [u] . (40)
Hence, the linearisation of the strong enhancement that in closing stage is written as:

-K s * b (k) n+1 ∆ d (k+1) n+1 -K s * w (k) n+1 ∆ [|ε|] (k+1) n+1 -(K s * s -K qc ) (k) n+1 ∆ [u] (k+1) n+1 = -Φ c (k) n+1
(41)

Solving the system

Depending on the status of the element is opening or closing, we can write the linearisation of the three equations into a matrix form:

  K bb K bw K bs K wb K ww K ws K s * b K s * w K s * s + K qo   (k) n+1    ∆ d ∆ [|ε|] ∆ [u]    (k+1) n+1 =    -( f e int -f e ext ) -R e [|ε|] -Φ o    (k) n+1 , (42)
or

  K bb K bw K bs K wb K ww K ws -K s * b -K s * w -K s * s + K qc   (k) n+1    ∆ d ∆ [|ε|] ∆ [u]    (k+1) n+1 =    -( f e int -f e ext ) -R e [|ε|] -Φ c    (k) n+1 . ( 43
)
The matrix is solved at two levels, global scale, and local scale. Since the strong kinematics enhancement Φ o and Φ c are both non-linear, Newton's method is practiced here for the resolution of the internal variables [|ε|] and [|u|]. It is assumed that the displacement field d remains constant during the resolution of the local system. So we have :

R e

[|ε|] = 0, (44a)

Φ o = 0 or Φ c = 0. (44b)

R e

[|ε|] = 0 has been developed in Eq. 36, we can see that it is a linear equation. The non-linear aspect of the local system originates from Φ o = 0 or Φ c = 0. Assuming that the equation ( 44) is now solved, the modified stiffness matrix K sc for the element in the opening procedure and in the closing procedure can be written as:

K sc (k) n+1 = K bb -K bw K bs K ww K ws K s * w K s * s + K qo (k) n+1 -1 K wb K s * b (k) n+1 , (45)
and

K sc (k) n+1 = K bb -K bw K bs K ww K ws -K s * w -K s * s + K qc (k) n+1 -1 K wb -K s * b (k) n+1 . (46)
The assembled matrix K sc is then calculated by bringing the stiffness matrix K sc of each element together:

K sc (k) n+1 = n e A e=1 K sc (k) n+1 , (47)
By using the static condensation, the resolving equations at global level writes as

K sc (k) n+1 ∆ d (k+1) n+1 = - n e A e=1 { f e int -f e ext } (k) n+1 . (48)
From Eq. (45) to Eq. ( 48), we can see that the stiffness matrix K sc has the same size of matrix K bb , which equals to the size of the global system no matter how many heterogeneities exist in the system or how many elements start to fail. The degree of freedom of the solving system is constant. In terms of numerical resolution, this feature brings a benefit that the computational solving memory will not increase with the increasing number of failed elements.

Resolution of the cohesive criterion

Both traction-opening process and the closing process in the strong discontinuity of mode-I are described in non-linear forms. In order to solve the equations at the local scale (Eq. 44), let us assume that the local system has already been solved, with R e

[|ε|] = 0 and R e

[|u|] = 0, which give:

K wb d + K ww [|ε|] + K ws [u] = 0, (49)
and

K s * b d + K s * w [|ε|] + K s * s [u] -σ y exp - σ y G op [u] = 0 for opening procedure (50a) -K s * b d -K s * w [|ε|] -K s * s [u] + G cl [u] max ln [u] [u] max = 0 for closing procedure. (50b)
By replacing all terms of [|ε|] by d and [u], the non-linear traction-opening equation can be solved as [START_REF] Roubin | Multiscale failure of heterogeneous materials: A double kinematics enhancement for embedded finite element method[END_REF]:

T e + M[u] = σ y exp(- σ y G op [u]), (51)
with

T e = K s * b -K s * w K -1 ww K wb d, and (52)
M = K s * s -K s * w K -1 ww K ws . (53)
A solution can only be solved if T e > σ y and M < 0. With the help of the Lamber W function [START_REF] Robert M Corless | On the lambertw function[END_REF] W 0 , an analytical solution can be deduced as:

[u] sol = G op σ y   W 0   σ 2 y exp( σ y T e G op M ) G op M   - σ y T e G op M   . (54)
As for the non-linear closing procedure, the cohesive criterion can be reformed as

-T e -M[u] = G cl [u] max ln [u] [u] max , (55)
and the analytical solution for [u] can be solved:

[u] sol = -G cl W 0   - M[u] 2 max exp( T e [u] max G cl ) G cl   (M[u] max ) . (56)
The resolution of the non-linear traction-separation criterion and the closing criterion is plotted separately in Fig. 10a and Fig. 10b.

[u] Once the opening value [u] sol is determined, the value of weak discontinuity [|ε|] sol can be analytically calculated as [START_REF] Roubin | Modélisation EF et morphologique de milieu hétérogènes à léchelle mésoscopique: applications aux matériaux à matrices cimentaire[END_REF] [|ε|] sol = -K -1 ww (K

σ y T e [u] sol Φ o M (a) Traction-opening procedure [u] -T e [u] sol [u] max Φ c -M (
wb d + K ws [u] sol ) (57)
Numerical application to a cubic specimen with heterogeneous structure

In this section, we will focus on the capabilities of the EFEM model by applying it to a heterogeneous specimen. The proposed numerical example is performed on a cube whose length is equal to 100 millimeters. Two phases are modeled, the heterogeneous structure is present by a set of spherical aggregates that disperse randomly in the cube, see Fig. 11. Two groups of aggregates are modeled to represent a total volume fraction of 20%, in which 60% of the spherical aggregates have a radius between 3 and 5 millimeters, and the rest of them are between 8 and 12 millimeters.

Fig. 11: Represent of the projection of the cube to mesh with aggregate disperse in the matrix and weakly enhanced elements between them.

Aiming at modeling the main features of brittle/quasi-brittle materials, Table . 1 lists the parameters of the two phases in the specimen, including the interfaces between them. Each of them contains two elastic parameters and two failure parameters. The aggregates are assumed to be more rigid than the matrix, and remain elastic. The cracks can only initiate and propagate in the matrix and interfaces. For the sake of simplicity, it is considered that the interface element carries the same failure parameters as the matrix element.

Table 1: Considered material parameters for the cube and its rigid inclusions.

Phase E [GPa] ν [-] σ y [MPa] G op [J/m  ] Matrix 22.0 0.2 4.0 1.0 Interface - - 4.0 1.0 Aggregate 78.0 0.2 - -
The specimen is then applied by mechanical loadings. Starting from the simple uniaxial loading. The macroscopic responses of the specimen are displayed in Fig. 12. In order to illustrate the effect of the closure mechanism, the model with and without crack closures are compared in this figure. From Fig. 12, we can see that firstly in the elastic phase, the two models carry the same value of macroscopic Young's module, which is equal to 27.4 GPa. Second, even though the failure criterion at the mesoscopic scale is proposed only in traction, the failure behavior at the macroscopic scale in compression is observed for both models. A reasonable ratio between the traction and compression resistance can be found between 13.05 and 13.49. Third, the difference between the two models' responses is slight in traction, whereas in compression, the model with closure mechanism shows a little higher resistance. The macroscopic response 0 1 2 3

•10 -4 0 1 2 3 (1) (2)
(3) (4)

Strain [-]

Zoom at the tensile responses

Model with closure Model without closure Fig. 12: The comparison between the model with/without closure law for the cube under monotonic loading.

To analysis the failure process in traction and the role of the closure mechanism, we choose here four loading stages, see Fig. 12. The crack patterns of these four stages are displayed in Fig. 13. The crack pattern in the post-peak stage is shown in Fig. 13b. We can see that the number of failed elements is heavily increased. The diffused cracks become denser and begin to coalesce to macroscopic cracks. A very small number of cracks can be found starting to close.

From Fig. 13c to Fig. 13d, the macroscopic response of the material tends to stabilize, whereas the crack pattern continues to change. The major crack continues to propagate and numerous crack closures occur in the specimen. Depending on the evolutions of major cracks, we can roughly classify the closing elements into two main categories. For the sake of clarity, three major cracks are presented here as examples, see Fig. 13d. The first closing category is present by crack 1 and 2. We can see that they are almost entirely closed from the stage (3) to stage (4). They are located at the upper or lower parts of other major cracks. The crack closures of this kind originate from the stress release which is caused by the rapid development of other major cracks (crack 3 in this case). For the model without closure mechanism, this stress-relieving effect can also be observed, see Fig. 14. Even though the cracks can not close for the model without crack closures, it can be seen that the cracks in the stress release area (crack 1 and 2) stop propagating from the stage (3) to stage (4) during the imposed loading increases.

The closing elements in crack 3 belong to another category of closure patterns. They are also triggered by stress release, but the stress release derives itself from the center of cracks, where the failed elements develop faster than the other elements. This observed phenomenon originates from the imposed model at the local scale, in which the equivalent stress on the discontinuity surface will decrease as the crack opening increases, meaning that the failed element propagates more and more easily. The consequence of this type of closure element is that, compared to the model without closure mechanism, the major crack in the model with closure mechanism will become narrower and clearer. The crack evolution of the two models is plotted and compared in Fig. 15. It can be observed from the figure that their main differences appear from the stage (3) to stage (4), i.e., after crack closures start to occur massively in the specimen for the model with closure mechanism. During this period, the number of elements with relatively higher crack opening values gradually increases for both models, corresponding to the propagation of major cracks. While for the number of elements with low crack opening values, the two models show different trends. There are a number of elements that admit a closing behavior for the model with closing mechanism, which leads to an increasing number of elements with small crack opening values. Whereas for the model without closing mechanism, this number remains almost constant.

Next, a loading-unloading displacement is applied to the specimen to illustrate the EFEM model's ability in terms of the unilateral effect. It is also one of the main characters of brittle/quasi-brittle materials. The macroscopic responses of the model are plotted in Fig. 16a. In Fig. 16a, we can see that for both models, the stiffness of the material can be fully recovered even the specimen is very damaged in traction. The reason for this recovery, though, is a little different between the two models. As for the model with-out crack closures, the stiffness of the material is recovered right after the traction loading switches to unloading. This is because the failure behaviors of the element are applied only on the discontinuity interface, the bulk volumes have always pure elastic kinematics. Hence the unloading displacement will be applied to these elastic bulk volumes since the crack openings can not be decreased. These crack openings are also the resource of the permanent plasticity that can be seen in the figure .  As for the model with closure mechanism, the macroscopic response of the specimen under unilateral loading is related to the closure of cracks. In the first phase of unloading, a number of failed elements admit a closure behavior, which leads to a partial recovery of the material stiffness, and also to an amount of energy dissipated in the closing phase. Then, as the imposed displacement recharges, the failed elements reopen, and a hysteresis loop can be observed. This phenomenon shows a good consistency with the experimental observations (Fig. 16b). And when the imposed displacement unloads to compression, we can see that the specimen has its stiffness fully recovered at the point (3). Some permanent plasticity can also be observed at point (2), but much less than for the model without closing mechanism. This is due to crack closures in the sample. Fig. 17: Variation of the cracks of the material during the unloading process in a unilateral test, the points are marked in Fig. 16a.

Interested in the crack patterns of the model with closure mechanism along with the unilateral loadings and its consistency with the experimental observation (Fig. 16b), we label here three stages as shown in Fig. 16a. Their crack patterns are displayed in Fig. 17. It can be seen that:

• At the point (1), a large amount of crack openings apparent in the specimen, accompanied with the loss of stiffness of the material. • At the point (2), the state of null stress, the cracks are partially closed while leading to an amount of plastic strain, the resource of which is the residual openings.

• At the point (3), the cracks are almost entirely closed, leading to the recovery of the stiffness of the specimen.

Finally, a compressive cyclic loading is applied to the material to analyze the performance of the model. The displacement-controlled trajectory is displayed in Fig. 18. and the macroscopic responses of the two models are given in Fig. 19. We can see from Fig. 19 that both models carry an amount of volumetric dilatation, and this phenomenon is more obvious pour the model without closure mechanism. This is because the volumetric dilatation derives itself from the crack openings. And the crack closures in the model with closure mechanism may restrict this effect.

Then in Fig. 20, several differences between the two models can be notices. First, we can see that the first cycle is purely elastic, and starting from the second cycle, the model with closure mechanism losses less of the macroscopic stiffness than the other model. Displayed in the figure by dashed lines, the model with closure mechanism has 42.3% of the initial stiffness at the end of loading, while the model without crack closures has only 18.5%. This can be explained by the crack closures in the unloading phase, which leads to the stiffness of the material partially recovers.

Second, it can be seen that the model with closure mechanism has slightly higher plasticity at the third cycle, which equals to ε p with closure = 2.2 × 10 -4 , while the value for the model without closure mechanism is ε p without closure = 2.0 × 10 -4 . Unlike in traction loadings that the plasticity is directly related to the residual crack opening, the observed plasticity in compression is the result of emergence from the mesoscopic to macroscopic scale since the closing procedure brings additional dissipated energy. It should be noted that at the local scale, no plastic deformation is formulated in the constitutive model.

Finally, the hysteresis phenomenon can be observed for the model with closure mechanism, which is one of the main characteristics of brittle/quasi-brittle materials. The source of this observed phenomenon comes from several ingredients, such as the capability of recovery stiffness, and to dissipate energy during the closing phase. And both of them are part of the effects of the closure mechanism.

Conclusion

In existing EFEM methods, many models have shown their ability to simulate many of the main characters of brittle or quasi-brittle materials. For example, the asymmetric behaviors in traction and compression loadings, the volumetric dilatation, and the progressive loss of stiffness. In this chapter, we applied a closure mechanism to the model to investigate if more behaviors can be observed. The model is tested to a concrete-like specimen with mesoscopic structural heterogeneities. Intricate crack pattern can be observed during the fracture process, and the propagations of cracks may influence each other. The crack closures can be observed not only in cyclic loadings , but also in monotonic loadings. The model shows also similar behaviors as experimental observations in unilateral loadings, including the hysteresis loops, irreversible deformations, and stiffness recovery. In compressional cyclic loading, The hysteresis phenomenon can be observed. Even the friction between the lips of cracks is not taken into concern, the fatigue behaviors can be seen by adding a closure mechanism to mode-I separation crack at the local scale. However, the hysteresis loops are not very significant. This may related with the lack of consideration of the friction between the cracks. In future works, it would be interesting to take this into concern.

  u the full enhancement of the virtual field corresponding to the strong discontinuity, d the nodal displacement field, and [|ε|] (resp.[|u|]) corresponds to the supplement unknown variables brought by weak (resp. strong) discontinuity.

Fig. 4 :

 4 Fig. 4: Local constitutive model at the continuous part and at the discontinuity interface at the moment of localization.

Fig. 5 :

 5 Fig. 5: Local constitutive model at the continuous part and at the discontinuity interface at the traction separation stage.

Fig. 7 :

 7 Fig. 7: Local constitutive model at the continuous part and at the discontinuity interface as the imposed displacement decrease.

Fig. 8 :

 8 Fig. 8: Local constitutive model at the continuous part and at the discontinuity interface at the crack closing phase.

Fig. 9 :

 9 Fig. 9: Local constitutive model at the continuous part and at the discontinuity interface as the imposed displacement reloads.

Fig. 10 :

 10 Fig. 10: Resolution of the strong discontinuity equation at the local scale.

Fig. 13 :

 13 Fig. 13: Cracks (in the middle) and dissipated energy at crack closures (in the right) at four different loading stages for simple traction, they are overlaid in the left, and the rigid inclusions are marked in gray.

Fig. 13a displays

  Fig.13adisplays the crack patterns at the loading stage[START_REF] Ahmed | An overview of complex adaptive systems[END_REF], which corresponds to the material just before the maximum resistance. A number of micro-cracks can be observed that diffuse over the cube, especially beside the inclusions because of the stress concentration. No crack closure can be observed at this phase.The crack pattern in the post-peak stage is shown in Fig.13b. We can see that the number of failed elements is heavily increased. The diffused cracks become denser and begin to coalesce to macroscopic cracks. A very small number of cracks can be found starting to close.From Fig.13cto Fig.13d, the macroscopic response of the material tends to stabilize, whereas the crack pattern continues to change. The major crack continues to propagate and numerous crack closures occur in the specimen. Depending on the evolutions of major cracks, we can roughly classify the closing elements into two main categories. For the sake of clarity, three major cracks are presented here as examples, see Fig.13d.

Fig. 14 :

 14 Fig. 14: Cracks for the model without closure law at the stage (4), and the increasing values of crack opening from the stage (3) to the stage (4).

Fig. 15 :

 15 Fig. 15: Distributions of crack opening values for the model with and without closure mechanism in traction.

Fig. 16 :

 16 Fig. 16: Illustration of the unilateral behaviors of concrete.

  (a) Point (1), ε M = 3.0 × 10 -4 (b) Point (2), ε M = 0.9 × 10 -4 (c) Point (3), ε M = -0.33 × 10 -4

Fig. 18 :

 18 Fig. 18: Proposed displacement path for the cyclic loading.

Fig. 19 :

 19 Fig. 19: Macroscopic response of monotonic/cyclic tests for model with/without adding closure law, in terms of axial (note Ax) transversal (note Tr) and volumetric strain (note Vol).

Fig. 20 :

 20 Fig. 20: Illustration of the loss of stiffness and the plastic deformation in cyclic loading for model with/without crack closures.

Yue SUN, Emmanuel ROUBIN, Jean-Baptiste COLLIAT, Jianfu SHAO