Pengfei Li

Yi Wu

Julien Yvonnet

Sili Liu

Shuitao Gu

Phase field modeling of dynamic fracture in elastoplastic composites with interfacial debonding

Keywords: Phase field method, Dynamic fracture, Interfacial debonding, Elastoplastic composites, Image-based models

In this work, we extend a phase field formulation for dynamic ductile fracture to consider interfacial debonding in elastoplastic composites. The interfacial weak zone is created through a regularization of the sharp interfaces, and the singular strain part along the interfaces is approximated by using a Taylor expansion. Then, a strain density depending on the displacement jump related to matrix/inclusions decohesion is added to the total energy to take into accout interfacial debonding. The coupling problems (displacement, plasticity and damage problems) are derived within the variational framework and a staggered iterative algorithmic procedure is described to solve the coupling problems. Numerical examples demonstrate that this method can handle the initiation, propagation, and interaction between bulk dynamic fracture and interface cracks, as well as the anisotropic behavior in the complex microstructure of elastic-plastic composite materials. It also indicates that this model is convergent in terms of mesh refinement.

Introduction

The computational modeling of dynamic fracture in elastoplastic composites is of essential importance for a wide range of applications including automotive, electronics, aerospace, and biotechnology engineering, e.g. evaluations of fracture resistance of 3D printed short fiber reinforced composites used in aerospace (see e.g. [START_REF] Yavas | Fracture behavior of 3d printed carbon fiber-reinforced polymer composites[END_REF][START_REF] Islam | Dynamic fracture behaviour of additively manufactured polymers and composites under ballistic impact[END_REF]), where components frequently suffer from dynamic loadings like wind, impact, and blast loads. In order to optimize dynamic fracture processing and evaluate dynamic fracture resistance, thereby maximizing the crack resistance of components, it is crucial to develop an efficient and robust numerical method for the modeling of the dynamic fracture processing.

In the past decades, many models and numerical methods have been proposed to simulate dynamic crack propagation, which can be simply classified into two categories: discrete and continuous approaches. In the former, the displacement field is allowed to be discontinuous on the fracture surface, while in the latter, the displacements are continuous everywhere, but the stresses are gradually decreased to describe the degradation process using some softening material models. The well-known theories for the discrete approaches are the discrete elements [START_REF] Iturrioz | Dynamic fracture analysis of concrete or rock plates by means of the discrete element method[END_REF][START_REF] Orozco | Discrete-element model for dynamic fracture of a single particle[END_REF], peridynamics [START_REF] Ha | Studies of dynamic crack propagation and crack branching with peridynamics[END_REF][START_REF] Kilic | An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory[END_REF][START_REF] Hu | Peridynamic model for dynamic fracture in unidirectional fiber-reinforced composites[END_REF], cohesive models [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Ožbolt | Dynamic fracture of concrete compact tension specimen: Experimental and numerical study[END_REF][START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF], and the extended finite element method (XFEM) [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF][START_REF] Combescure | X-fem a good candidate for energy conservation in simulation of brittle dynamic crack propagation[END_REF][START_REF] Nguyen-Vinh | Extended finite element method for dynamic fracture of piezo-electric materials[END_REF]. Continuum damage mechanics is probably the most widely used theory categorized in the continuous fracture approach [START_REF] Chen | Continuum damage mechanics studies on the dynamic fracture of concrete[END_REF][START_REF] Lee | Continuum damage modeling for dynamic fracture toughness of metal matrix composites[END_REF].

From a computational perspective, most of the above mentioned approaches suffer in situations with complex crack patterns including branching, giving substantial difficulties in numerical implementations, and often restrict the above formulations to simple crack pattern. This difficulty can be overcome by the variational approach to fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF][START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF], also called phase field method in the literature, which regularizes sharp crack discontinuities within a pure continuum formulation, and offers several advantages, such as the possibility to initiate fracture from undamaged zones, to handle arbitrary fracture patterns without specific treatment and use of classical finite elements, a variational framework allowing to involve many models or mechanisms, and a mesh-independence due to a gradient description of damage. Detailed information for the phase field method can be found e.g. in the review papers [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF][START_REF] Wu | Phase-field modeling of fracture[END_REF].

Phase field method has been recently extended to ductile fracture behaviors by coupling gradient damage mechanics with models of elasto-plasticity. These models can be simply classified into two families: weak coupling and strong coupling models. A well-known weak coupling model was proposed in Ambati et al. [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF], where a characteristic degradation function that couples damage to plasticity in a multiplicative format is used. Another typical weak coupling model, one can refer to [24], where the coupling is based on a hardening degradation function related to the equivalent plastic strain which works on the initial fracture toughness. The above mentioned models, in spite of losing the variational consistance, can well simulate the ductile fracture in elastoplastic solids (see [25]). The strong coupling models e.g. in the papers [26][27][28][29], through introducing a damage based plastic energy into the total energy, allow the plastic energy to induce the crack initiation and propagation. The strong coupling models follow the variational framework and allow to include many models or mechanisms. Extensions of the strong coupling model to gradient plasticity at finite strains were considered in [30][31][32][33] based on a rigorous variational principle. A comparative study for different phase field models of fracture coupled with plasticity was outlined in [34].

With respect to dynamic fracture, phase field models have been considered in [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Larsen | Existence of solutions to a regularized model of dynamic fracture[END_REF][START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF][START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF][START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF][START_REF] Nguyen | Modeling dynamic fracture of solids with a phase-field regularized cohesive zone model[END_REF][START_REF] Liu | Virtual modelling integrated phase field method for dynamic fracture analysis[END_REF], just to name a few. These models are based on the assumption that the kinetic energy is not affected by the phase field and the fracture energy is independent of the crack velocity. The dynamic crack branching phenomenon experimentally observed can be captured and mesh independent fracture energies were obtained e.g. in [START_REF] Geelen | A phase-field formulation for dynamic cohesive fracture[END_REF]. Currently, most of these works focused on the numerical simulation of the dynamic crack networks whereas Bleyer et al. [START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF][START_REF] Bleyer | Microbranching instability in phase-field modelling of dynamic brittle fracture[END_REF] presented some explanations of dynamic fracture behavior.

More recently, phase field method has been extended to the coupling of ductile and dynamic to simulate the dynamic crack initiation and propagation in elastoplastic solids. In [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part ii. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF] modeling of dynamic fracture in the logarithmic Lagrangian strain space has been presented with emphasis on analysing the transition of brittle-to-ductile failure mode. In [START_REF] Molnár | An open-source abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation[END_REF], Molnár et al. proposed a robust and versatile open-source Abaqus subroutine for practical engineers as well as research scientists in a widely available commercial finite element code, which can be employed to study both brittle and ductile static and dynamic fracture. To this end, Wu and Shen [START_REF] Wu | An elastoplastic phase-field model for dynamic fracture of nickel-base super-alloys[END_REF] proposed a new phase field model for dynamic fracture in elastoplastic solids, which is achieved through the use of new plastic degradation function that allows for a transition between two dissipative stages.

In most of the above mentioned approaches, interfacial contributions were not included in the constitutive formulation, despite the fact that interfaces in composites may introduce additional weak zones and can strongly affect the damage mechanisms. In this respect, Paggi et al. [START_REF] Paggi | Revisiting the problem of a crack impinging on an interface: a modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model[END_REF] presented a novel combined approach that employs the phase field model to handle bulk brittle fracture and simulates interface debonding using classical cohesive elements. Msekh et al. [START_REF] Msekh | Predictions of j integral and tensile strength of clay/epoxy nanocomposites material using phase field model[END_REF][START_REF] Msekh | Fracture properties prediction of clay/epoxy nanocomposites with interphase zones using a phase field model[END_REF] applied the phase field approach to simulate progressive failure in polymer-matrix composites. Hansen-Dörr et al. [START_REF] Hansen-Dörr | Phase-field modelling of interface failure in brittle materials[END_REF] proposed a phase field approach to simulate interface failure between two dissimilar materials. In [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF][START_REF] Nguyen | Phase field modeling of interfacial damage in heterogeneous media with stiff and soft interphases[END_REF], a simple framework that extends phase field to consider interfacial debonding has been proposed, where a regularized description of the sharp interface and displacement jumps around the interface has been used, and different damage behavior to the bulk and to the interface in the energy is associated with an appropriate weighting of indicator functions.

In this work, we propose an extension of the phase field model as presented in [START_REF] Molnár | An open-source abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation[END_REF] for the simualtion of dynamic fracture in elastoplastic composites. The involved extensions comprise: (i) considering interfacial debonding between the inclusions and the matrix; (ii) modeling interactions between interfacial debonding and bulk dynamic ductile fracture; (iii) handling elastoplastic composite microstructures obtained from X-ray computed tomography (XRCT) technique. The interfacial weak zone is created through a regularization of the sharp interfaces, and the singular strain part along the interfaces is approximated by using a Taylor expansion. Then, a strain density depending on the displacement jump related to matrix/inclusions decohesion is added to the total energy to take into accout interfacial debonding. To maintain the regularized property of the approximation, smooth indicator functions are employed to weight the different terms in the energy with respect to the vicinity of interfaces. Then, the coupling problems, involving displacement, plasticity and damage problems, are derived and a staggered iterative algorithmic procedure is described to solve the coupling problems.

The paper is structured as follows. In Section 2, we introduce a regularization of bulk crack, interfacial damage, displacement jump and plasticity. Section 3 provides the phase field modeling of dynamic fracturing with interfacial damage in elastoplastic composites. Section 4 provides the FEM discretization and the flowchart for the overall algorithm. Section 5 provides the numerical examples. Finally, conclusions are drawn in Section 6.

Regularization of discontinuities and plasticity

Consider an open domain Ω ⊂ R D with D = 2, 3, describing a heterogeneous body which contains internal interfaces between two elasto-plastic phases. The external boundary of Ω is denoted by ∂Ω ∈ R D-1 , composed of Dirichlet and Neumann parts, denoted by ∂Ω u and ∂Ω F , ∂Ω = ∂Ω u ∪ ∂Ω F , ∂Ω u ∩ ∂Ω F = ∅, where displacements and tractions are prescribed, respectively. During the loading, crack may propagate in the matrix and can pass through the interfaces, where the crack surface and the interfaces are collectively denoted by Γ and Γ I (see Fig. 1a), respectively. In this work, we employ the framework presented in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via tconvergence[END_REF] for a regularized representation of discontinuities and extension to interfaces as in [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]. In this regularized framework, the cracks are approximately represented by a scalar phase field d (x, t) (see Fig. 1b) and the interfaces by a fixed scalar function β (x) (see Fig. 1d), the plastic response is regularized by using a gradient plasticity formulation (see Fig. 1c), and the displacement jump around the interface is regularized by using a Taylor expansion (see Fig. 1e).

Regularization of bulk cracks

For a known fixed crack surface Γ (see Fig. 1a), the scalar crack phase field d (x, t) can be determined via solving the following boundary value problem (see [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] for more details):

           d (x, t) -ℓ 2 d ∆d (x, t) = 0 in Ω, d (x, t) = 1 on Γ, ∇d (x, t) • n = 0 on ∂Ω, (1)
where ∆ (•) and ∇ (•) are the Laplacian and gradient operator respectively, ℓ d is a length scale parameter to govern the width of crack regularization zone and gives for ℓ d → 0 the exact sharp crack as shown in Fig. 1b, n is the outward normal on ∂Ω. It can be shown that (1) is the Euler-Lagrange equation for the following variational problem:

d = Arg inf d∈S d Γ d (d) , Γ d (d) = Ω γ d (d) dΩ, S d = {d | d (x) = 1, ∀x ∈Γ} , (2)
where Γ d (d) represents the total length of the crack in 2D and the total crack surface area in 3D, and γ d (d) is the crack surface density function defined by:

γ d (d) = d 2 2ℓ d + ℓ d 2 ∇d • ∇d, (3)
where the second term in γ d (d) penalizes high values of ∇d(x) to prevent well-known related non-convergence issues with respect to the mesh discretization, which are found in local damage model.

Regularization of plasticity

With a similar manner, the plastic response can be regularized (see Fig. 1c) through introducing a gradient part into the initial plastic strain energy density function [30] as

ψ p 0 = σ y p + 1 2 Hp 2 + 1 2 σ 0 ℓ 2 p ∇p • ∇p, (4)
where σ y and H > 0 are the yield stress and the hardening modulus parameters, σ 0 is a parameter to keep unit consistency for the gradient part and σ 0 = 1 MPa used in the following, p denotes an equivalent plastic strain described by the evolution equation

ṗ = 2 3 ∥ εp ∥ , (5)
where ( .) is time derivative and ε p is the plastic strain. In (4), ℓ p is a plastic length scale related to a strain-gradient hardening effect and takes into account size effects to prevent the nonphysical mesh sensitivity of the localized plastic deformation in softening materials.

Regularization of interfaces

Similarily, the scalar interface function β (x) can also be determined through solving the following boundary value problem

           β (x) -ℓ 2 β ∆β (x) = 0 in Ω, β (x) = 1 on Γ I , ∇β (x) • n = 0 on ∂Ω, (6)
where ℓ β is a length scale parameter to govern the width of the regularization zone of the interface and gives for ℓ β → 0 the exact sharp interfaces as shown in Fig. 1d. [START_REF] Kilic | An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory[END_REF] corresponds to the Euler-Lagrange equation for the following variational problem

β = Arg inf β∈S β Γ β (β) , Γ β (β) = Ω γ β (β) dΩ, S β = β | β (x) = 1, ∀x ∈Γ I , (7)
where γ β (β) is defined by

γ β (β) = β 2 2ℓ β + ℓ β 2 ∇β • ∇β. (8)

Regularization of displacement jumps within interfaces

In the present work, a regularized approximation for the displacement jump [[u]] is introduced to consider the interface debonding as presented in [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]. In order to derive a description of the displacement jump [[u]], we first define the indicator function ϕ (x) as

ϕ (x) = β (x) , (9)
then the displacement jump [[u]] can be regularized through using a Taylor expansion of the displacement field around a point x located on the interface (see Fig. 1e)

[[u]] ≃ w (x) = u x+ h 2 n I -u x- h 2 n I = h∇ (u (x)) n I , (10)
where n I is the normal vector to Γ I at the point x and is defined by:

n I = ∇ϕ (x) ∥∇ϕ (x)∥ , (11)
w (x) denotes the regularized displacement jump which is not only defined on the interface but over all the domain, and h is a small scalar parameter.

The traction acting on the interface oriented by n I (see Fig. 1e) is defined by

t (w) = ∂ψ I (w) ∂w , (12)
where ψ I is the interface strain density depending on the regularized displacement jump w. Note that as discussed in [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF], it is not required to introduce history variables into this formulation for the interfaces, and the diffuse damage field is used to describe the irreversibility of the interfacial debonding. For 2D problems, the traction t (w) has the following form

t (w) = t n , t t T , (13)
where t n and t t denote normal and tangential parts of the traction vector t, respectively. In this work, the tangential component in the interface model is neglected (t t = 0), and the normal component has the following form [START_REF] Van Den Bosch | An improved description of the exponential xu and needleman cohesive zone law for mixed-mode decohesion[END_REF]:

t n = g I c w n (δ n ) 2 exp - w n δ n , (14)
where δ n is related to the interface fracture toughness g I c and the interface fracture strength t u by

δ n = g I c t u exp (1) , ( 15
)
w n is the normal displacement jump defined by w n = w•n I . Two examples of relationship ( 14) are shown in Fig. 2.

It should be noted that β (x) and n I do not change during the simulation due to the invariable of interfaces. For a sharp interface, the singular part of the strain along the interface can be defined by then using the above regularized formulation, ( 16) can be approximated as [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]:

ε I (x) = n I (x)⊗ s [[u(x)]] δ(x) ∀x ∈ Γ I , (16)
ε I (x) ≃ n I (x)⊗ s w(x)γ β (x) ∀x ∈ Ω. (17)
It is worth noting that the above three length scale parameters (ℓ d , ℓ p and ℓ β ) do not represent physically the exact crack or plastic response width, but the parameters which are used to regularize the discontinuities and plasticity. In the present work, the length scale parameter ℓ d and ℓ β are simply chosen according to ℓ d = ℓ β = 2h e , where h e is the minimal element size, h is chosen according to h = h e to obtain a minimal error for displacement jump (see [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]), ℓ p ≥ ℓ d is required such that the regularized crack zone lies inside of the plastic zone as outlined in [30], here ℓ p = 3h e is used. Interestingly, it has been shown that ℓ d can be treated as a material parameter associated with the Young's modulus, the tensile strength, and the critical energy release rate of the material in [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Zhang | Numerical evaluation of the phase-field model for brittle fracture with emphasis on the length scale[END_REF]. Possible improvements may imply the use of recent length-free formulations as e.g. proposed in [START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF], where a length-scale insensitive phase field method for brittle fracture has been proposed.

Phase field modeling of elastoplastic dynamic fracture interacting with interfacial damage

In this section, we propose a phase field model for modeling interactions between interfacial debonding and bulk dynamic cracking in elastoplastic composites. The formulation is developed using the variational framework for fracture as introduced in [26].

First, the total strain is defined as:

ε = ε e + ε p + ε I , (18)
where ε I has been defined in [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF] and ε e denotes the elastic strain. The plastic incompressibility is assumed, i.e. Tr(ε p ) = 0, Tr(.) being the trace operator. In the following, we define the deviatoric parts of the stress and of the elastic strain tensors as s = dev(σ) and e e = dev(ε e ), respectively, with dev(.) = (.) -1 3 Tr(.)1, 1 being the first-order identity operator. 7

Energy functional

Let us consider a two-phase medium which is elasto-plastic and suffers from dynamic loading. Small strains are assumed. The action-integral over a time interval I = [t 1 , t 2 ] in the absence of body forces can be defined as

A = I Ω ψ e + ψ p + ψ d + ψ I dΩ - Ω ψ k dΩ - ∂Ω F F • udS dt, (19)
where ψ e denotes elastic strain energy density, ψ p denotes the plastic strain energy density, ψ d is the fracture energy density, ψ I is a strain density depending on the approximated displacement jump across the interfaces and ψ k is the kinetic energy density, F is prescribed traction over the boundary ∂Ω F .

Here, the elastic strain energy density is chosen as

ψ e = g (d) ψ e+ 0 (ε e ) + ψ e- 0 (ε e ) , (20)
with [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]

ψ e+ 0 (ε e ) = 1 2 κ ⟨Tr (ε e )⟩ 2 + + µe e : e e , (21)
ψ e- 0 (ε e ) = 1 2 κ ⟨Tr (ε e )⟩ 2 -, (22)
where, κ and µ denote bulk modulus and shear modulus, respectively,

⟨x⟩ ± = 1 2 (x ± |x|), g(d) = (1 -d) 2 + k,
k being a small numerical parameter to prevent loss of definite posedness of the elastic tensor for fully broken case.

The plastic strain energy density is here chosen as

ψ p = g (d) ψ p 0 , (23)
where ψ p 0 denotes the initial plastic strain energy density as given in [START_REF] Orozco | Discrete-element model for dynamic fracture of a single particle[END_REF]. Note that the presence of g (d) in ( 23) introduces a strong plasticity-damage coupling into the present model, where the plastic strain energy can also induce the crack initiation and propagation.

The fracture energy density is chosen as

ψ d = (1 -β) 2 ψ d 0 , (24)
where ψ d 0 denotes the initial fracture energy density. For ψ d 0 , one choice is

ψ d 0 = g c γ d (d) , (25)
where γ d is defined in (3) and g c denotes the fracture toughness, another choice proposed by Fremond [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF] can provide a threshold

ψ d 0 = ψ c 2d + ℓ 2 d ∇d • ∇d . ( 26
)
Note that compared to (25), in (26) d enters the formulation by a linear term. In (26), ψ c is a specific fracture energy density, which can be further related to a critical fracture stress σ c by

ψ c = 1 2E σ 2 c , ( 27
)
where E is the Young's modulus. In the present work, the initial fracture energy density in ( 26) is used in the following. Following our previous work [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF], we give the energy density of interfacial jump component as

ψ I = ψ I (w) γ β (β) , (28)
where γ β and w are defined in ( 8) and [START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF], respectively. ψ I is the interface strain density given in [START_REF] Combescure | X-fem a good candidate for energy conservation in simulation of brittle dynamic crack propagation[END_REF].

Here, the kinetic energy density has this form

ψ k = 1 2 ρ u • u, ( 29
)
where ρ is the material density and u denotes the velocity.

Variational formulation

In this section, we use the variational framework for fracture as introduced in [26,[START_REF] Rodriguez | A variational approach to the phase field modeling of brittle and ductile fracture[END_REF][START_REF] Li | Anisotropic elastoplastic phase field fracture modeling of 3d printed materials[END_REF]] to derive the evolution criteria of the proposed model. The variational framework involves: irreversibility condition, stability condition, energy balance and alternate minimization.

Irreversibility condition

The irreversibility condition for fracture and plastic variables is expressed as

ḋ ≥ 0, 0 ≤ d ≤ 1, ṗ ≥ 0, (30)
the fracture irreversibility condition is prescribed using a history function [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] (see Section 3.2.4), while the plastic irreversibility condition is applied numerically by simply considering the equivalent plastic strain value corresponding to the previous load step as the minimum admissible level of equivalent plastic strain for a given position in the body.

First-order stability condition

We first define the directional derivative as:

D v f (u) = d dh f (u + hv) h=0 . (31)
The first order stability condition (see [START_REF] Mielke | Evolution of rate-independent systems[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF]) is expressed by

D δu A + D δp A + D δd A ≥ 0, (32)
which yields by using ( 19)

I Ω σ : ε e (δu) + ∂ψ I ∂w •w (δu) + ρü • δu dΩ - ∂Ω F F • δudS dt + I Ω - 3 2 σ : nδp + D δp ψ p dΩ dt+ I Ω ∂ψ e ∂d δd + ∂ψ p ∂d δd + D δd ψ d dΩ dt ≥ 0 ( 33
)
where n is a unit tensor in the direction of the plastic flow and the Cauchy stress is given by

σ = ∂ψ e ∂ε e = g (d) σ + eff + σ - eff , (34)
where σ + eff and σ - eff are the effective tensile and compressive stresses defined as

σ + eff = ∂ψ e+ 0 ∂ε e = κ ⟨Tr (ε e )⟩ + 1 + 2µe e , σ - eff = ∂ψ e- 0 ∂ε e = κ ⟨Tr (ε e )⟩ -1. (35)
From (33), the following results stem out:

• For δp = δd = 0, the time interval I = [t 1 , t 2 ] is arbitrary the integrand must sum to zero so that finally

Ω σ : ε e (δu) dΩ + Ω γ β t•w (δu) dΩ + Ω ρü • δudΩ - ∂Ω F F • δudS = 0 ( 36
)
which is the weak form of the equilibrium equation.

• For δd = 0 and δu = 0, we first express the term:

D δp ψ p = g (d) σ y δp + Hpδp + σ 0 ℓ 2 p ∇p • ∇δp , (37)
considering that the time interval I = [t 1 , t 2 ] is arbitrary, (33) yields:

Ω 3 2 ∥s∥ -g (d) (σ y + Hp) δp -g (d) σ 0 ℓ 2 p ∇p • ∇δp dΩ ≤ 0, (38)
which is the weak form of the plasticity yield criterion. The local form of the plastic yield criterion can be expressed as

F p = 3 2 ∥s∥ -g (d) σ y + Hp -σ 0 ℓ 2 p ∇ • ∇p ≤ 0 in Ω. (39)
• Similarily for δp = 0 and δu = 0, (33) yields:

Ω (1 -d) ψ e+ 0 + σ y p + 1 2 Hp 2 + 1 2 σ 0 ℓ 2 p ∇p • ∇p -(1 -β) 2 ψ c δd -(1 -β) 2 ψ c dδd + ψ c ℓ 2 d ∇d • ∇δd dΩ ≤ 0, (40)
which is the weak form of the damage yield criterion. The local form of the damage yield criterion can be expressed as

F d = (1 -d) ψ e+ 0 + σ y p + 1 2 Hp 2 + 1 2 σ 0 ℓ 2 p ∇p • ∇p -(1 -β) 2 ψ c -(1 -β) 2 ψ c d -ψ c ℓ 2 d ∇ • ∇d ≤ 0 in Ω. (41)

Energy balance

Following a procedure analogous to the treatment of the stability condition, the energy balance condition leads to

I - Ω σ : ε e ( u) + ∂ψ I ∂w •w ( u) + ρü • u dΩ + ∂Ω F F • udS dt + I Ω 3 2 σ : n ṗ -D ṗψ p dΩ dt - I Ω ∂ψ e ∂d ḋ + ∂ψ p ∂d ḋ + D ḋψ d dΩ dt = 0 (42)
Considering that the time interval I = [t 1 , t 2 ] is arbitrary, the following cases are analyzed:

• For u = 0 and ḋ = 0, the plasticity consistency condition is obtained:

F p ṗ = 0 (43)
• For u = 0 and ṗ = 0, the damage consistency condition can be written in the form:

F d ḋ = 0 ( 44
)
where F p and F d are given in ( 39) and ( 41), respectively.

Alternate minimization

In this section, a staggered alternate minimization algorithm for u, p and d is applied. With the total energy [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF] at hand, the alternate minimization follows.

• Minimizing with respect to the displacement field:

D δu A = I Ω σ : ε e (δu) dΩ + Ω γ β t•w (δu) dΩ + Ω ρü • δudΩ - ∂Ω F F • δudS dt = 0 (45)
which fits arbitrary I = [t 1 , t 2 ] and thus leads to

R u = Ω σ : ε e (δu) dΩ + Ω γ β t•w (δu) dΩ + Ω ρü • δudΩ - ∂Ω F F • δudS = 0 (46)
which is the weak form of the mechanical problem to be solved for u, given p and d.

• Minimizing with respect to the equivalent plastic strain, we can then obtain:

R p = Ω 3 2 ∥s∥ -g (d) (σ y + Hp) δp -g (d) σ 0 ℓ 2 p ∇p • ∇δp dΩ = 0, (47)
which is the weak form of the plastic yield criterion to be solved for p with ṗ ≥ 0, given u and d. The increment of the plastic strain tensor can be obtained from the incremental equivalent plastic strain as

εp = 3 2 ṗ s ∥s∥ (48)
• Minimizing with respect to the damage field, we can then obtain:

R d = Ω (1 -d) ψ e+ 0 + σ y p + 1 2 Hp 2 + 1 2 σ 0 ℓ 2 p ∇p • ∇p -(1 -β) 2 ψ c δd -(1 -β) 2 ψ c dδd + ψ c ℓ 2 d ∇d • ∇δd dΩ = 0. ( 49
)
To prescribe damage irreversibility, a history function inspired by [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] is given as

H = max s∈[0,t] ψ e+ 0 (ε e , s) + σ y p(s) + 1 2 Hp(s) 2 + 1 2 σ 0 ℓ 2 p ∇p(s) • ∇p(s) -(1 -β) 2 ψ c + , (50)
then we can rewrite [START_REF] Hansen-Dörr | Phase-field modelling of interface failure in brittle materials[END_REF] as

R d = Ω (1 -d) Hδd -(1 -β) 2 ψ c dδd + ψ c ℓ 2 d ∇d • ∇δd dΩ = 0, (51)
which corresponds to the global damage problem to be solved to find the field d(x), given u and p.

Discretization and numerical implementation

In this section, we detail the finite element discretizations for displacement, plastic and damage problems, and finally provide the overall algorithm. It should be noticed that, in this paper, due to the regularization of plasticity, the plastic problem is handled through solving a global problem for equivalent plastic strain p as shown in [START_REF] Rodriguez | A variational approach to the phase field modeling of brittle and ductile fracture[END_REF], unlike our previous work [START_REF] Li | An extension of the phase field method to model interactions between interfacial damage and brittle fracture in elastoplastic composites[END_REF], where the plastic problem is treated as a local one (at Gauss integration points) and solved by using the return-mapping algorithm.

Displacement problem

In this work, for the sake of clarity, only 2D FEM discretization is detailed. The vector form of second-order tensors are introduced as

[ε] = ε 11 , ε 22 , √ 2ε 12 T , [σ] = σ 11 , σ 22 , √ 2σ 12
T , as well as the FEM approximations

u = N u u e , [ε] = B u u e , w = hNBu e , [ε e + ε p ] = B w u e , (52)
where u e denotes nodal displacement components in one element, N u is a matrix of displacement shape function, B u is a matrix of displacement shape function derivatives, and

N=   n 1 n 2 0 0 0 0 n 1 n 2   , B=   ∂ ∂x 1 ∂ ∂x 2 0 0 0 0 ∂ ∂x 1 ∂ ∂x 2   T N u , (53)
where n 1 and n 2 are the x-and y-components of the normal vector n I in (11), B w is a modified spatial strain-displacement matrix defined by

B w = B u -hγ β QB (54)
in which Q is a matrix expressed by

Q =       n 1 0 0 n 2 1 √ 2 n 2 1 √ 2 n 1       N =       n 2 1 n 1 n 2 0 0 0 0 n 1 n 2 n 2 2 1 √ 2 n 1 n 2 1 √ 2 n 2 2 1 √ 2 n 2 1 1 √ 2 n 1 n 2       . ( 55
)
After rearranging, the discrete version of ( 46) can be expressed globally as

Mü + K u u = F u (56)
with

M = Ω ρN T u N u dΩ, K u = Ω B T w CB w dΩ (57)
F u = Ω B T w Cε p,(k) dΩ + ∂Ω F N T u FdS - Ω γ β (hNB) T t w (k) dΩ (58)
where

C =    g (d) (κII + 2µII d ) Tr (ε e ) ≥ 0 κII + g (d) 2µII d Tr (ε e ) < 0 (59)
II =       1 1 0 1 1 0 0 0 0       , II d =       1 2 -1 2 0 -1 2 1 2 0 0 0 1       . (60)
Above, (.) (k) denotes the value of k-th subiteration for each time step.

Then, an unconditionally stable implicit Newmark scheme is employed to solve [START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF] according to un = un-1 + ∆t 2 ün-1 + ün [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF] with

u n = u n-1 + ∆t un-1 + ∆t 2 2 ün-1 + ün (62) ün = 2 ∆t 2 (u n -ûn ) (63) ûn = u n-1 + ∆t un-1 + ∆t 2 2 ün-1 , (64)
where (.) n denotes the value of n-th time step and ∆t = t n -t n-1 is assumed to be constant.

Introducing (63) into [START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF], the linear problem to be solved for displacement at subiteration k of one time step t n can be obtained as

Kn u u n = Fn u (65) with Kn u = 2 ∆t 2 M + K n u (66) Fn u = F n u + M 2 ∆t 2 u n-1 + ∆t un-1 + ∆t 2 2 ün-1 (67)

Plastic problem

We first rewrite [START_REF] Msekh | Fracture properties prediction of clay/epoxy nanocomposites with interphase zones using a phase field model[END_REF] with small pseudo-time step assumption as

∆ε p = 3 2 ∆p s ∥s∥ (68)
where ∆p = p n -p n-1 is the incremental equivalent plastic strain. Then we can rewrite [START_REF] Msekh | Predictions of j integral and tensile strength of clay/epoxy nanocomposites material using phase field model[END_REF] as

R p = Ω g (d) 3 2 s tr eff -3µp + 3µp n-1 -σ y -Hp δp -σ 0 ℓ 2 p ∇p • ∇δp dΩ = 0, (69)
where s tr eff denotes a non-plastic trial stress. Here, for the sake of simplicity, the superscript n for p n is omitted. Notice that Eq.( 69) is based on an explicit Euler algorithm and thus require a small time step to ensure accuracy and convergence.

The equivalent plastic strain and equivalent plastic strain gradient are approximated in one element by

p = N p p e , ∇p = B p p e , (70)
where N p and B p are matrices of equivalent plastic strain shape function and of equivalent plastic strain shape function derivatives, respectively, and p e denotes nodal equivalent plastic strain in one element. The discretization of plastic problem (69) results into the following discrete system of equations

K p p = F p (71)
with

K p = Ω g (d) (3µ + H) N T p N p + σ 0 ℓ 2 p B T p B p dΩ (72)
F p = Ω g (d) N T p 3 2 s tr eff + 3µp n-1 -σ y dΩ. ( 73
)
It is noted that the present work solves the gradient plasticity as a global problem by finiteelement method, naturally avoiding unphysical oscillations at the elastic-plastic boundary (reported in [31]) inducing by directly using the return-mapping algorithm in finite-element implementations.

Damage problem

The damage and damage gradient are approximated in one element by

d = N d d e , ∇d = B d d e , (74)
where N d and B d denote matrices of damage shape function and of damage shape function derivatives, respectively, and d e denotes nodal damage in one element. The discretization of damage problem (51) results into the following discrete system of equations

K d d = F d (75)
with

K d = Ω H+ (1 -β) 2 ψ c N T d N d + (1 -β) 2 ψ c ℓ 2 d B T d B d dΩ (76)
F d = Ω N T d H dΩ, ( 77
)
where H is given in (50).

Overall algorithm

In the present work, a staggered iterative scheme is employed following [START_REF] Rodriguez | A variational approach to the phase field modeling of brittle and ductile fracture[END_REF]. It consists of solving successively (65), ( 71) and (75). First, (65) is solved to obtain u n,(k+1) with given u n,(k) , p n,(k) , d n,(k) , n and k denoting the time step and subiteration, respectively. In a second step, (71) is solved to obtain p n,(k+1) with given u n,(k+1) , p n,(k) , d n,(k) . In a third step, (75) is solved to obtain d n,(k+1) with given u n,(k+1) , p n,(k+1) , d n,(k) . At each time step, the subiteration is performed until convergence is achieved by using three independent tolerances t u , t p and t d . The detailed algorithm is illustrated in Algorithm 1.

Algorithm 1: Overall algorithm

Initialize u 0 , p 0 and d 0 with assumption of elasticity and undamaged state; for n = 1, . . . , m (Loop over all time steps) do k = 0, erru = errp = errd = 1 u n,(0) = u n-1 , p n,(0) = p n-1 , d n,(0) = d n-1 while erru > t u or errp > t p or errd > t d and k ≤ k max do Given u n,(k) , p n,(k) and d n,(k) , compute u n,(k+1) from (65) Given u n,(k+1) , p n,(k) and d n,(k) , compute p n,(k+1) from (71) Given u n,(k+1) , p n,(k+1) and d n,(k) , compute d n,(k+1) from ( 75)

erru = u n,(k+1) -u n,(k) , errp = p n,(k+1) -p n,(k) , errd = d n,(k+1) -d n,(k) k = k + 1 end u n = u n,(k+1) , p n = p n,(k+1) , d n = d n,

Numerical examples

In this section, all numerical computations are performed within the finite element framework. For the sake of clarity, we recall that the length scale parameter ℓ d and ℓ β are simply chosen according to ℓ d = ℓ β = 2h e , where h e is the minimal element size, h is chosen according to h = h e to obtain a minimal error for displacement jump, ℓ p ≥ ℓ d is required such that the regularized crack zone lies inside of the plastic zone, here ℓ p = 3h e is used.

Dynamic crack branching

In this example, we consider a notched block loaded dynamically in tension. As shown in Fig. 3, a tensile load of σ = 1.0 MPa is applied at the top and bottom edges of the specimen and held constant throughout the course of the simulation. This problem has been widely adopted to study dynamic crack branching in previous works, including [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF] and [START_REF] Molnár | An open-source abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation[END_REF], to name a few.

The material parameters used in the simulations are taken from [START_REF] Molnár | An open-source abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation[END_REF]: Young's modulus 32000 MPa, Poisson's ratio of 0.2, the mass density of ρ = 2450 kg/m 3 , the yield stress σ y = 2 MPa, hardening modulus H = 0 MPa.

2D model with plain strain assumption

In order to verify our implementation, we compare our fracture pattern to ones obtained in [START_REF] Molnár | An open-source abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation[END_REF]. For this purpose, we do not set interface (β = 0) in the specimen as shown in Fig. 3. The spatial discretization comprises 36196 4-node quadrilateral elements, with refinement in the right region (the minimal element size h e = 0.25 mm) where the crack is expected to propagate (see Fig. 3). The time step is set to ∆t = 0.1 µs. Plain strain conditions are assumed.

The final fracture patterns obtained from different σ c as well as reference result in [START_REF] Molnár | An open-source abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation[END_REF] are shown in Fig. 4. We illustrate that the fracture patterns obtained by this proposed model have a good agreement with the reference result especially for σ c = 10 MPa. Thus, σ c = 10 MPa is used in the following unless otherwise stated. The effect of the critical fracture stress σ c on fracture energy and plastic energy is shown in Fig. 5. As can be observed, with the decrease of the critical fracture stress σ c , both the fracture initiation time and the final failure time decrease, the maximum fracture energy and plastic energy tend to decrease, besides the width of the fracture zone shows a decreasing trend. Further, we analyze the influence of the time increment ∆t in the numerical simulation. For this purpose, another two simulations are performed with ∆t = 0.05 µs and ∆t = 0.2 µs, respectively. The corresponding fracture patterns and fracture/plastic energy curves are presented in Fig. 6 and Fig. 7, respectively. The results confirm the stability of the easy-to-implement staggered algorithm as soon as sufficiently small time steps are used. We notice that ∆t = 0.1 µs is small enough for the numerical simulations, thus ∆t = 0.1 µs is used in the following unless otherwise stated.

3D model

To eliminate the influence of plain strain assumption used in 2D model and better investigate the performance of the proposed model, a 3D model with the same geometry as shown in Fig. 3 whereas thickness 0.25 mm in the z direction is investigated. The spatial discretization of the 3D model comprises 237734 tetrahedral elements, with refinement in the right region (h e = 0.25 mm) as shown in Fig. 8d.

Then we continue to investigate the dynamic crack branching by employing the above 3D model and the same parameters used in the 2D model, notice the interface is still not considered. Fig. 8a shows the final fracture pattern, as can be seen, the fracture pattern exhibits an obvious difference with the 2D fracture patterns shown in Fig. 4. This phenomenon can also be found in the reference [START_REF] Molnár | An open-source abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation[END_REF]. Next we investigate the effect of the yield stress σ y on the final crack patterns and fracture/plastic energy. The corresponding results are shown in Fig. 8 and Fig. 9, respectively.

It can be seen that, as the yield stress is increased, the final crack pattern trends to become the result obtained by quasi-brittle dynamic fracture model. In order to better consider plasticity, we continue to use σ y = 2 MPa in the following unless otherwise stated. Further, we investigate the influence of the interface by setting an interface in the 3D model as shown in Fig. 10a. For this purpose, two sets of parameters for the cohesive model (defined in ( 14)), as illustrated in Fig. 2, are considered to illustrate the effect of the interfacial effect on crack patterns and fracture/plastic energy. The mesh model shown in Fig. 8d is used in the numerical simulations. The corresponding results are shown in Fig. 10 and Fig. 11, respectively. We can observe that, compared to no interface, the presence of the interface has an obvious influence on the crack patterns and fracture/plastic energy, and when considering interfacial damage and changing the cohesive parameters, both the fracture patterns and fracture/plastic energy are very similar. Then, we test the convergence of the model with respect to the mesh refinement. For this purpose, the previous simulation used the parameters in model I is conducted with another two refined meshes with 362976 elements and 498246 elements, respectively. Fig. 12 shows the final fracture patterns for the three mesh models. Fig. 13 shows the fracture/plastic energy for the three mesh models. The results show the convergence of the method with mesh refinement.

3D rectangular plate with cylindrical fibers

In this example, we investigate the proposed model through a complex structure which has several cylindrical fibers, as shown in Fig. 14. The dimensions of the 3D structure in x-y plane are 80 × 40 mm 2 , and in order to ensure the used element number within an appropriate range, the thickness of the 3D structure is chosen as 0.25 mm. As can be observed from Fig. 14, there are five cylindrical fibers named by "F1-F5 ", whose diameters are 20 mm (F1), 10 mm (F2-F4) and 8 mm (F5), respectively. The coordinates of the central axis for the five cylindrical fibers in the basis (x, y) are [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF], [START_REF] Wu | An elastoplastic phase-field model for dynamic fracture of nickel-base super-alloys[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF], [START_REF] Mielke | Evolution of rate-independent systems[END_REF]30), [START_REF] Mielke | Evolution of rate-independent systems[END_REF][START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF] and (70, 20), respectively. Similarly, a tensile load of σ = 2.0 MPa is applied at the top and bottom edges of the specimen and held constant throughout the course of the simulation. The spatial discretization of the 3D model comprises 356938 tetrahedral elements with the minimal element size h e = 0.25 mm, as shown in Fig. 14. The material properties are shown in Table 1, and the cohesive parameters for model I as shown in Fig. 2 are used. The numerical simulations are conducted with ∆t = 0.2 µs in the following. Fig. 15 shows the fracture process at different times. We can observe that cracks nucleate from the interface and then kink into the matrix when reaching a certain point. Subsequently, these interface cracks and matrix cracks are interconnected and then lead to a complex fracture pattern. The evolution of fracture energy and plastic energy is shown in Fig. 16. Fracture energy Plastic energy Further, we use this complex structure to further investigate the influence of interface and cohesive parameters. With the similar manner as shown in Section 5.1.2, another two numerical simulations with cohesive model II and no-interface are implemented, respectively. Fig. 17 shows the fracture patterns for the three simulations. We notice that there is not crack initiation until t = 200 µs when interfaces are not considered. Fig. 18 shows the comparison of fracture/plastic energy obtained with model I and model II. The results demonstrate again that, interface has an obvious influence on the crack patterns and fracture/plastic energy, while cohesive parameters seem do not have obvious effect on the fracture patterns and fracture/plastic energy when interface is considered. Finally, we further investigated the convergence of the model with respect to the mesh refinement. For this purpose, the previous simulation is conducted with another two refined meshes with 247364 elements and 485672 elements, respectively. Fig. 19 shows the final fracture patterns for the three mesh models. Fig. 20 shows the fracture/plastic energy for the three mesh models. The results illustrate again the convergence of the method with mesh refinement.

3D complex microstructure obtained from microtomography

In this example, we investigate the capabilities of the proposed method to simulate microcrack propagation in 3D complex microstructure obtained from microtomography such as obtained by experimental imaging techniques, like X-ray computed tomography (XRCT) technique. The microstructure is extracted from a 3D printed fiber reinforced elastoplastic composite sample and has the dimensions 0.178 × 0.178 × 0.178 mm 3 , as shown in Fig. 21a. The spatial discretization of the 3D microstructure comprises 709105 tetrahedral elements (see Fig. 21b) with the minimal element size h e = 0.002 mm, where pore phases are not meshed to save the computational costs. The material parameters are shown in Table 1, and the cohesive parameters for model I as shown in Fig. 2 are used. The boundary conditions are as follows (see Fig. 21c): on the lower end (x = 0), the (x, y, z)-components are fixed, on the upper end, the (y, z)-components are free, while the x-velocities are prescribed by

v (t) =    v 1 t 0 t if t ≤ t 0 v 1 otherwise (78)
for t 0 = 0.02 µs and v 1 = 100 mm/s. The numerical simulations are conducted with ∆t = 0.01 µs and named as "x-loading"in the following. Fig. 22 shows the fracture process at different times, where for the sake of clear visualization, only the crack phase field with values higher than 0.9 is plotted. As can be seen, cracks are initiated in the form of interface debonding and then migrate into the matrix in the form of matrix cracks. Subsequently, these interface cracks and matrix cracks are connected and then lead to the final failure of the microstructure. The similar damage phenomena can also be found in 5.2. The evolution of fracture energy and plastic energy is shown in Fig. 23. Fracture energy Plastic energy Finally, we investigate the capabilities of the proposed method to qualitatively capture the anisotropic behavior of 3D printed microstructure. For this purpose, two additional simulations are performed: (i) in the first one, called "y-loading", the boundary conditions are as follows: on the plane y = 0, the (x, y, z)-components are fixed, on the plane y = max(y), the (x, z)-components are free, while the y-velocities are prescribed by (78), and (ii) in the second one, called "z-loading", the (x, y, z)-components are fixed on the plane z = 0, on the plane z = max(z), the (x, y)-components are free, while the z-velocities are prescribed by (78). Fig. 24 shows the final fracture patterns for the three simulations. Fig. 25 shows the comparison of fracture energy and plastic energy. The results show an obvious difference and demonstrate the capabilities of the proposed method to qualitatively capture the anisotropic behavior of 3D printed microstructure. x-loading y-loading z-loading

Conclusion

In this work, we have extended a dynamic phase field model based on the work [START_REF] Molnár | An open-source abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation[END_REF], to take into account the following features: (i) interfacial debonding in the elastoplastic composites; (ii) interactions between interfacial debonding and bulk dynamic ductile fracture; (iii) the possibility to handle elastoplastic composite microstructures obtained from X-ray computed tomography (XRCT) technique. In order to consider interfacial damage, the total energy has been modified by adding a strain density depending on the displacement jump, where interfaces are regularized by a length scale parameter and singular strain part along the interfaces is approximated by using a Taylor expansion. In that manner, different damage mechanisms can be associated with interfaces as compared to the matrix cracking mechanisms. Different equations (displacement, plasticity and damage) can be obtained by using the variational principle and then be solved in a staggered iterative procedure. As a result, this extension allows simulating initiation, propagation and interactions between both bulk dynamic fracture and interfacial debonding in elastoplastic composites without specific treatment and use of classical finite elements. Several numerical examples, involving complex microstructures e.g. obtained from 3D printed short fiber reinforced composites, have shown the capability of the method to handle complex micro cracks interactions for arbitrary complexity of the microstructures, and the convergence with respect to the mesh refinement. To our best knowledge, the presented simulations involving ductile dynamic fracturing, considering interfacial debonding and realistic microstructure directly obtaining from an experimental XRCT have been presented here for the first time. Then, the present model seems to be very promising for predicting initiation, propagation of complex microcracking in a dynamic-mechanical context in elastoplastic composites, such as 3D printed fiber reinforced composites used in aerospace. One future perspective of this work would be extended to adaptive mesh refinement to deal with the high computaional cost causing by the use of a large number of elements in undamaged zones.

project "Multiscale Modeling and Experimental Investigation of damage in composite components obtained by additive manufacturing process" (MMELED), ANR-16-CE08-0044-04.

Figure 1 :

 1 Figure 1: Regularization of discontinuities and plasticity: (a) a body containing sharp cracks and interfaces; (b) diffused approximation of cracks; (c) regularization of plastic response; (d) diffuse approximation of interfaces; (e) regularization of displacement jump.

Figure 2 :

 2 Figure 2: Two cohesive models for analyzing the influences of the interfacial parameters.

Figure 3 :

 3 Figure 3: Dynamic crack branching: geometry of the edge cracked block (left) and random finite element mesh with a refined zone (he = 0.25 mm) in the expected propagation zone (right).

Figure 4 :Figure 21 : 21 Figure 5 :

 421215 Figure 4: 2D dynamic crack branching: final fracture patterns for (a) σc = 20 MPa, (b) σc = 15 MPa, (c) σc = 10 MPa and (d) reference [44]. [24] T. T. Nguyen, Modeling of complex microcracking in cement based materials by combining numerical simulations based on a phase-field method and experimental 3d imaging, Ph.D. thesis, Paris Est (2015). [25] T. X. Le, Experimental study on the mechanical properties and the microstructure of methane hydrate-bearing sandy sediments, Ph.D. thesis, Paris Est (2019). [26] H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments, Journal of the Mechanics and Physics of Solids 57 (8) (2009) 1209-1229. [27] M. Ambati, T. Gerasimov, L. De Lorenzis, Phase-field modeling of ductile fracture, Computational Mechanics 55 (5) (2015) 1017-1040. [28] C. Miehe, F. Aldakheel, A. Raina, Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory, International Journal of Plasticity 84 (2016) 1-32. [29] R. Alessi, J.-J. Marigo, S. Vidoli, Gradient damage models coupled with plasticity: variational formulation and main properties, Mechanics of Materials 80 (2015) 351-367. [30] P. Rodriguez, J. Ulloa, C. Samaniego, E. Samaniego, A variational approach to the phase field modeling of brittle and ductile fracture, International Journal of Mechanical Sciences 144 (2018) 502-517. [31] C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Computer Methods in Applied Mechanics and Engineering 199 (45-48) (2010) 2765-2778. [32] F. P. Duda, A. Ciarbonetti, P. J. Sánchez, A. E. Huespe, A phase-field/gradient damage model for brittle fracture in elastic-plastic solids, International Journal of Plasticity 65 (2015) 269-296. [33] F. Cuvelier, C. Japhet, G. Scarella, An efficient way to perform the assembly of finite element matrices in matlab and octave, arXiv preprint arXiv:1305.3122. [34] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.

Figure 6 :Figure 21 : 21 Figure 7 :

 621217 Figure 6: 2D dynamic crack branching: final fracture patterns for (a) ∆t = 0.2 µs, (b) ∆t = 0.1 µs and (c) ∆t = 0.05 µs.

Figure 8 :Figure 21 : 21 Figure 9 :

 821219 Figure 8: 3D dynamic crack branching: final fracture patterns for (a) σy = 2 MPa, (b) σy = 3 MPa and (c) σy = 4 MPa; (d) 3D mesh model with a refined zone (he = 0.25 mm).

Figure 10 :Figure 21 : 21 Figure 11 :

 10212111 Figure 10: 3D model with interface: (a) 3D model, (b) fracture pattern obtained by model I, (c) fracture pattern obtained by model II and (d) fracture pattern without considering interface.

Figure 12 :

 12 Figure 12: 3D model with interface. The final fracture patterns for (a) mesh 1 (237734 elements), (b) mesh 2 (362976 elements) and (c) mesh 3 (498246 elements).

Figure 21 : 21 Figure 13 :

 212113 Figure 21: A larger image-based model. Comparison of load-displacement curves for different loading direction.

Figure 14 :

 14 Figure 14: 3D rectangular plate with cylindrical fibers: geometry (left) and 3D mesh model (right).

Figure 15 :

 15 Figure 15: 3D rectangular plate with cylindrical fibers: crack patterns at (a) t = 8 µs, (b) t = 14 µs, (c) t = 18 µs and (d) t = 24 µs.

Figure 22 :

 22 Figure 22: A larger image-based model. Comparison of load-displacement curves for different loading direction.

Figure 16 :

 16 Figure 16: 3D rectangular plate with cylindrical fibers: evolution of fracture energy and plastic energy.

Figure 17 :Figure 21 : 21 Figure 18 :

 17212118 Figure 17: 3D rectangular plate with cylindrical fibers: fracture patterns obtained from (a) model I, (b) model II and (c) no-interface.

Figure 19 :Figure 21 : 21 Figure 20 :

 19212120 Figure 19: 3D rectangular plate with cylindrical fibers. The final fracture patterns for (a) mesh 1 (247364 elements), (b) mesh 2 (356938 elements) and (c) mesh 3 (485672 elements).

Figure 21 :

 21 Figure 21: 3D complex microstructure obtained from microtomography: (a) whole microstructure model (size: 0.178 × 0.178 × 0.178 mm 3 ), (b) mesh model and (c) boundary conditions.

Figure 22 :

 22 Figure 22: 3D complex microstructure obtained from microtomography: crack patterns at (a) t = 0.25 µs, (b) t = 0.30 µs and (c) t = 0.45 µs.For the sake of clear visualization, only the crack phase field with values higher than 0.9 is plotted.

Figure 22 :

 22 Figure 22: A larger image-based model. Comparison of load-displacement curves for different loading direction.

22

 22

Figure 23 :

 23 Figure 23: 3D complex microstructure obtained from microtomography: evolution of fracture energy and plastic energy.

Figure 24 :

 24 Figure 24: 3D complex microstructure obtained from microtomography: final fracture patterns for (a) x-loading, (b) y-loading and (c) z-loading. For the sake of clear visualization, only the crack phase field with values higher than 0.9 is plotted.

Figure 21 :

 21 Figure 21: A larger image-based model. Comparison of load-displacement curves for different loading direction.

21

 21

Figure 25 :

 25 Figure 25: 3D complex microstructure obtained from microtomography: comparison of fracture energy (left) and plastic energy (right) for three different loading directions.

Table 1 :

 1 Material parameters used in the numerical simulations

	Name	Symbol Matrix Fiber	Unit
	Young's modulus	E	32	60	GPa
	Poisson's ratio	ν	0.2	0.25	[-]
	Mass density	ρ	2450	4750 kg/m 3
	Yield stress	σ y	4	8	MPa
	Hardening modulus	H	2	4	MPa
	Critical fracture stress	σ c	10	30	MPa

Acknowledgements

Pengfei Li gratefully acknowledges the financial support from Jiangsu Open University Fourteenth Five-Year Plan Project with project number 2022XK008 and from Jiangsu Education Department with project number 23KJB130002. Yi Wu acknowledges support from the China Postdoctoral Science Foundation (Grant No. 2023M730217). This work has also benefited from the

Conflict of interest statement

On behalf of all authors, the corresponding author states that there is no conflict of interest.


Joseph Doumit

Marie-Josepha Youssef

Abdel Charbel

Joumana Nour

Catherine Farah

Douillard

Uncoordinated Spectrum Access for Massive Connectivity in IoT Applications

Boosted by the Internet of Things (IoT), future cellular communication networks will be able to support new services with high social impact. Corresponding applications are expected to span almost every domain, especially the most critical ones, such as the medical field or smart factories. However, their adoption is conditioned by strict compliance with stringent quality of service (QoS) requirements for the underlying network. In particular, the achieved throughput and latency are of utmost importance for efficient monitoring with multiple sensors in a variety of application areas. In this paper, we propose an intelligent self-configurable system structure, achieved through resource allocation techniques using machine learning, and capable of flexibly fulfilling such requirements. In addition, special care was taken to minimize the transmit power of the IoT devices while efficiently supporting different application-dependent traffic types, an essential prerequisite for any wide adoption of such systems with massive network access.

INTRODUCTION

Internet of Things (IoT) devices are expected to account for the largest part of mobile traffic in the fifth generation (5G) and beyond 5G (B5G) networks [START_REF] Al-Sarawi | Internet of Things Market Analysis Forecasts, 2020-2030[END_REF]. Following a trend already started with 5G, the sixth-generation (6G) is expected to go beyond people-to-people communication to connect devices, sensors, vehicles, computing resources, and robotic agents [START_REF] Giordani | Toward 6G Networks: Use Cases and Technologies[END_REF]. Therefore, 6G will require technologies enabling the support of a massive deployment of IoT devices with different requirements in terms of data rate, latency, and reliability. To do so, various techniques can be smartly combined to pave the way for scalable and self-configurable massive IoT systems.

To increase system capacity for massive systems, the concept of non-orthogonal multiple access (NOMA) has been introduced to support more devices than the number of available orthogonal resources in time and frequency. NOMA allows multiple devices to share the same time and frequency resources via power domain or code domain multiplexing [START_REF] Dai | A Survey of Non-Orthogonal Multiple Access for 5G[END_REF]. Regardless of the adopted NOMA technique, it increases the receiver complexity, required to recover the non-orthogonal signals. In this paper, we focus on power-domain NOMA in the uplink and will be denoted by NOMA in the following for brevity.

In addition to NOMA, relaying can be useful in massive communication settings. Indeed, it can reduce network congestion at the base station (BS) level, extend cellular coverage and increase achieved rates thanks to shorter point-to-point transmissions. Mobile relays that can be deployed quickly and efficiently, such as drones, can also restore communications in damaged network infrastructures [START_REF] Youssef | Full-Duplex and Backhaul-Constrained UAV-Enabled Networks Using NOMA[END_REF]. To achieve high data rates, reduce latency and increase spectral efficiency, fullduplex relays, simultaneously transmitting to the BS and receiving from IoT devices on the same frequency resource, Joseph Doumit, Charbel Abdel Nour and Catherine Douillard are with IMT Atlantique, Lab-STICC, UMR CNRS 6285, F-29238 Brest, France. Marie-Josepha Youssef is with Capital Fund Management, Paris, France. Joumana Farah is with INSA Rennes, France. can be implemented. However, full-duplex relays suffer from self-interference, which must be mitigated [START_REF] Youssef | Full-Duplex and Backhaul-Constrained UAV-Enabled Networks Using NOMA[END_REF].

On a different note, the development of scalable networks having variable numbers of devices and traffic types is essential for unleashing the potential of IoT. Uncoordinated or grantfree spectrum access strategies can be suited to this context. Indeed, since such strategies do not require any coordination at the level of the nodes in the network, signaling overhead and communication latency can be significantly reduced, compared to conventional access schemes [START_REF] Zucchetto | Uncoordinated Access Schemes for the IoT: Approaches, Regulations, and Performance[END_REF]. Moreover, to reduce planning and maintenance costs, such networks can be equipped with self-organizing capabilities, enabling efficient resource sharing through power and spectrum allocation in a distributed manner. To that end, techniques such as reinforcement learning (RL) [START_REF] Sutton | Reinforcement Learning, Second Edition -An Introduction[END_REF] can be particularly appealing.

Hence, by using NOMA in a grant-free communication setting with relays, systems supporting scalable and selfconfigurable massive access can be efficiently devised. However, such complex systems need to be thoroughly designed in order to meet the quality of service (QoS) requirements for IoT applications. The design of such a system represents the main focus of this work.

Prerequisites and System Components

1) NOMA: For an efficient uplink NOMA, the multiplexed devices are allocated different power levels according to their channel conditions and are directly superimposed after classic channel coding and modulation [START_REF] Choi | NOMA-Based Random Access With Multichannel ALOHA[END_REF]. Figure 1 illustrates the principle of NOMA for the uplink. At the receiver, multiuser detection algorithms such as successive interference cancellation (SIC) are used to cancel multi-user interference and separate the different signals. SIC exploits the signal-tointerference-plus-noise ratio (SINR) difference between users to recover the transmitted signals. To ensure successful SIC, the signals of devices transmitting on the same channel must be received with different power levels. Therefore, if two or more devices transmitting on the same channel choose the same power level, a collision occurs, leading to a failed signal decoding [START_REF] Doumit | Resource Allocation in Full-Duplex Uncoordinated Communication Systems with NOMA[END_REF]. In the case of a collision, devices with signals having weaker received power levels than the colliding signals cannot be recovered either. These penalizing repercussions can impact upper-level signals if their target SINR is not met. Therefore, signal collisions should be minimized.

2) RL: RL is a framework where one or multiple learning agents interact with an environment in order to achieve a certain goal, usually represented by the maximization of a numerical reward. Knowing the environment status, the learning agent chooses its action based on a given strategy. Depending on the chosen action, it receives a reward and moves to a new state. This cycle is then repeated. That said, actions that achieve the highest rewards are not known beforehand. This motivates a two-phase approach consisting of exploration and exploitation phases [START_REF] Sutton | Reinforcement Learning, Second Edition -An Introduction[END_REF]. The exploration phase intends to progressively adapt and define a rewards mapping based on the actions undertaken. This mapping depends on the average of previous rewards obtained during a sliding window and the strategy of the agent. Then, the exploitation phase is the phase where the agent continuously chooses the action having the highest reward for future plays. Many studies resort to the multi-armed bandits (MAB), a RL framework applied in a single-state environment [START_REF] Youssef | Stochastic Multi-Player Multi-Armed Bandits with Multiple Plays for Uncoordinated Spectrum Access[END_REF], to solve the problem of channel allocation in ad-hoc networks and treat the exploration and exploitation phases separately to learn the channels [START_REF] Wang | Decentralized Learning for Channel Allocation in IoT Networks Over Unlicensed Bandwidth as a Contextual Multi-Player Multi-Armed Bandit Game[END_REF] and to manage the strategy for edge users [START_REF] Wang | Decentralized Task Offloading in Edge Computing: A Multi-User Multi-Armed Bandit Approach[END_REF].

3) Relaying and RL: One of the main drawbacks of the MAB framework resides in the necessity to perform an exploration phase for each variation of one of the essential system parameters, e.g., the number of active devices at the relay level or the transmission conditions including the channel model. To circumvent the need for an exploration phase when the number of active devices varies, we propose to introduce relays equipped with queuing capabilities to help prevent or at least reduce data loss in case of system congestion, hence efficiently managing sporadic traffic. This can be of particular interest since it enables a higher level of flexibility and scalability for the MAB framework in managing the underlying network congestion. However, since it requires memory buffers, queuing tends to increase the transmission link latency, which can considerably hinder the performance of some IoT applications. To this effect, we propose to introduce specific latency-related design constraints for the scheduling algorithm deployed at the relays. To evaluate the efficiency of the proposed algorithms, the age of information (AoI) is investigated at the receivers. AoI is a metric used to measure the information freshness when it reaches the destination [START_REF] Talak | Improving Age of Information in Wireless Networks With Perfect Channel State Information[END_REF]. It is defined as the time elapsed since the latest useful piece of information, that reached its intended destination, has been generated at its source. AoI is particularly relevant in applications in which information is time-varying and where the timeliness of information is crucial, such as vehicular monitoring systems or industrial sensor networks [START_REF] Ni | Vehicular beacon broadcast scheduling based on age of information (AoI)[END_REF]. The target framework accounts for the AoI constraints, depending on the type of considered traffic.

Contributions

Previous research has not considered the utilization of fullduplex relays in an uncoordinated uplink NOMA system with transmit power limitations and AoI restrictions for devices with heterogeneous traffic types. Thus, in this work, we investigate such a system. The focus is on enabling IoT devices to autonomously organize their communications with available relays in an uncoordinated spectrum access environment. The system includes a mixture of traffic types, including continuous, sporadic, and periodic transmissions. Furthermore, each IoT device operates under a limited power budget. We propose a multi-step algorithm for uplink massive IoT communication systems capable of supporting different traffic types while leveraging new intelligent self-configurable network architectures. The main contribution of this article concerns the proposal of efficient methods to:

• enable a flexible and scalable communication network that handles devices having multiple traffic types with latency and power constraints while reducing the data loss and the congestion at the level of the relays in a distributed communication system, • allocate the available spectrum and power resources to the devices in an uncoordinated manner while maximizing achieved rates and minimizing transmit power, • reduce the AoI for the IoT devices having latency constraints by using the MAB framework.

SYSTEM MODEL

In the following, we consider an uplink communication system with a BS, full-duplex relays uniformly deployed over the cell and a certain number of IoT devices, as shown in Fig. devices with sporadic traffic are assigned a probability of access to the spectrum. Moreover, it is assumed that the totality of a sporadic-type message can be transmitted during one time slot. The relevance of these three types of traffic can be illustrated with an application example in an e-health context: the remote monitoring of a medical act. One can imagine that the intervention is filmed by cameras that send a continuous stream, whereas specific sensors send periodic information about the patient (e.g., blood pressure) and vital sign changes are remotely managed by sensors sending sporadic control signals.

We also consider that, due to inconvenient radio conditions between the BS and the IoT devices with respect to the available transmit power, the devices cannot transmit their signals directly to the BS. Hence relays are introduced to forward the messages to the BS. A two-hop communication network is, therefore, set up with two types of links: links from IoT devices to relays (relay uplinks) and links from relays to BS (backhaul links). We assume that all transmission links undergo large-scale fading, due to path loss and shadowing, and the corresponding channel gain is assumed to be perfectly estimated for all links.

As mentioned earlier, spectrum access is fully uncoordinated: IoT devices must therefore arrange communications without exchanging information with each other and without any coordination at the level of the BS or the relays. To do so, each IoT device is allowed to communicate with one relay over one frequency subband.

For relay uplinks, multiple access is provided through NOMA. At the relay, a limited number of received power levels is allowed and SIC decoding is performed in the decreasing order of power levels. The transmitting devices are assumed to have an exact knowledge of their channel gains with the relays as well as the received power levels at the relays in order to control their transmit power level. The number of received power levels at the relays is the result of a trade-off. On one hand, IoT devices are power-limited. On the other hand, there is an interest in having as many power levels as possible in order to satisfy the requirements in terms of spectral efficiency and massive connectivity. Let L be the number of received power levels available per channel at the relay. We also consider that continuous traffic has priority over periodic and sporadic traffic. For this purpose, the L -1 highest power levels are dedicated to devices delivering continuous traffic and the L-th weakest level is only available for non-continuous (sporadic/periodic) traffic. This NOMA model ensures that the minimum rate required for continuoustraffic devices is guaranteed even if multiple collisions occur on the lowest power level. The maximum number of handled collisions depends on the average number of active devices generating non-continuous traffic.

As already mentioned, each relay operates in full-duplex mode. It receives the NOMA superimposed signals from the IoT devices, decodes these signals using SIC, and stores them in the relay memory buffer. The contents of the buffer are then transmitted to the BS on the same frequency band over a number of time slots that depends on the QoS requirements of the IoT devices, the number of correctly decoded signals at the relay, the backhaul channel conditions and the power budget of the relay.

For backhaul links, each relay is assigned a frequency subband, and a conventional orthogonal multiple access is implemented. Thereby, backhaul links do not suffer from interrelay interference. Moreover, since the relays operate in a fullduplex mode, the same subbands are used in the uplink and backhaul links.

With the proposed system model, the overall instantaneous AoI of an IoT device is computed as the summation of the AoIs for the relay uplinks and backhaul links that carry the data sent by the device. For the relay uplink, the AoI depends on the number of collisions the message undergoes before being successfully decoded at the relay. It can be reduced by optimizing the subband and power level allocations and by minimizing the number of signal collisions. For the backhaul link, AoI is directly related to the congestion level of the buffer. It is measured as the time that a message has to be stored in the buffer before being sent to the BS.

The main features of this system model are summarized in Table I. A major problem to be solved in uncoordinated access systems is the allocation of communication resources to the different devices of the system, especially frequency subbands and power levels, taking the specificity of their traffic type into account. The QoS requirements are different for devices generating continuous traffic and those transmitting data in periodic or sporadic mode. The proposed resource allocation approach is formulated using the MAB framework, and solved with one of the most common algorithms called the upper confidence bound (UCB) algorithm [START_REF] Adjif | Joint Channel Selection and Power Control for NOMA: A Multi-Armed Bandit Approach[END_REF]. In the MAB framework, agents are called players and executing an action consists in playing an arm. In the UCB algorithm, one has first to define the set of players, the set of arms and the rewards. In our system, the players are the IoT devices. Each device has to transmit on one frequency subband/power level combination. Hence, the set of arms is the set of available combinations of frequency subbands/power levels or, equivalently, relays/power levels. The nature of the rewards depends on the type of devices. Since devices generating continuous traffic are mainly dedicated to the more rate-hungry applications, their reward is the rate assigned to the device. For non-continuous traffic, we focus on applications requiring relatively low throughput but strict latency. In this case, the reward also includes the AoI. In addition, other constraints have to be taken into account, such as the power budget of the devices and relays. Thus, the algorithms are proposed in uncoordinated communication systems to achieve system stability and minimize the AoI. The allocation of the resources for continuous traffic-type devices guarantees system stability with no continuous overflow, whereas the allocation for non-continuous traffic-type devices is designed to meet the AoI constraints.

CONTINUOUS TRAFFIC: SYSTEM STABILITY

The first objective is to set up a stable system, where the relay buffers do not suffer from persistent overflows. To do so, the rate requirements of devices generating continuous traffic should not exceed the capacity of the backhaul link.

The proposed solution is based on the UCB approach, where the IoT devices with continuous traffic form the set of players. In practice, since each IoT device has a specific power budget and because of possible bad radio conditions, some power levels may not be available for all devices. Therefore, each device forms its own set of available arms from the L -1 highest power levels according to its power budget and the radio channel conditions to the relays. In our system, the reward received by a device when it plays an arm is the rate achieved by the device, given by Shannon's formula as developed in [START_REF] Doumit | Resource Allocation in Full-Duplex Uncoordinated Communication Systems with NOMA[END_REF].

In case two or more devices play the same arm, a collision occurs, resulting in a zero reward for devices with nondecodable signals. The main challenge is then the choice of the decision function that makes a device choose a frequency subband/power level pair. Every IoT device with continuous traffic aims to achieve a transmission rate at least equal to the requested one while spending the least amount of transmit power. Thus, two parameters should be considered in the design of the decision function: the rate achieved by the device and its transmit power which is the chosen received power level inversely proportional to the square of the channel gain. The minimum value of the transmit power is reached when the device selects the radio channel with the highest gain, i.e., the best radio channel. However, if each device chooses its best arm with respect to its channel gain, system convergence cannot always be guaranteed. Hence, IoT devices should choose their subband/power level pair to strike a compromise between the transmit power and the achieved rate. In practice, the trade-off is obtained via a weight β, whose value depends Start Coordinated Setup? Ranking of subband/power levels in terms of the channel conditions Ranking of IoT devices in terms of the queue status at the relays Application of Gale-Shapley algorithm [START_REF] Gale | College Admissions and the Stability of Marriage[END_REF] for matching resources to devices Create set of actions for each device Select arm with the highest decision value using (1)

Receive reward computed as a function of the ratio Achieved Rate/AoI Update UCB parameters Timeslot = T?

End

Yes No No Yes Fig. 3: Flowchart of the applied algorithms in coordinated and uncoordinated settings.

on system requirements and which can be adjusted by the system operator:

decision = (1 -β) × Reward -β × P ower. (1)
In addition, in order to prevent the overflow of buffers at the relay level due to a mismatch between the aggregated continuous traffic rate and the capacity of the backhaul link and, therefore, to guarantee the stability of the system, a specific procedure has been implemented. When buffering occurs and the memory contents reach a specified threshold, a certain number of IoT devices with continuous traffic type are invited to transmit to a less congested relay to avoid any overflow. However, since the system is uncoordinated, no information is shared between the devices and the involved relay. Therefore, the devices need to learn alone that they need to choose another relay. To do so, when the relay has reached its filling threshold, a zero reward (instead of the achieved rate) is assigned to all the connected devices for a limited period, called timeout period. Its duration is a function of the system parameters and the number of connected devices. During the timeout period, the IoT devices connected to the relay are invited to run a new exploration phase and to discover new suitable arms. By playing new arms associated with less congested relays, the devices reach a new stable state of the system able to satisfy the required constraints. Once reached, a new exploitation phase can begin.

NON-CONTINUOUS TRAFFIC: AOI MINIMIZATION

The proposed solution is also based on the UCB approach, where the IoT devices with non-continuous traffic types form the set of players.

The set of arms is different from that of devices with continuous traffic because only the L-th and weakest power level of each subband is available for IoT devices with non-continuous traffic. However, some subbands may not be available for some devices since periodic and sporadic transmissions need to efficiently share the available resources. In fact, collision management, channel quality and available device power budget may introduce restrictions related to subband access.

As a major shift from the case of continuous traffic devices, the reward for non-continuous traffic devices also accounts for the AoI. In practice, we have implemented a reward function proportional to the achieved rate and inversely proportional to the AoI. Hence, if a relay is congested, the AoI of connected devices tends to increase over successive time slots. Then, the corresponding reward decreases during these time slots and the device is encouraged to play new arms associated with less congested relays.

Significant improvement can still be achieved by devising a reward dependent upon a predicted AoI for the next time slot instead of the measured current value. In fact, with a deterministic access periodicity for continuous and periodic traffic devices, this prediction aims at estimating the access probability of devices with sporadic traffic. When coupled with the filling status of the relay, estimating the number of noncontinuous devices connected in the next time slot allows the relay to predict whether it can avoid congestion or whether it should encourage the connected devices with sporadic traffic to connect to other relays. This incentive is achieved by controlling the value of the reward. The algorithms developed throughout the article are summarized in Fig. 3.

RESOURCE ALLOCATION USING COORDINATED MATCHING

In order to assess the efficiency of the proposed uncoordinated access approach, a coordinated allocation technique based on matching theory is introduced to benchmark the performance of the proposed algorithm. While lacking the self-configuration capability of the proposed approach, this reference technique is feasible in practice at the cost of comparatively high signaling overhead.

Matching theory is a mathematical framework in economics [START_REF] Gale | College Admissions and the Stability of Marriage[END_REF] that describes the formation of mutually beneficial relationships. In particular, it is commonly used to solve assignment problems between two disjoint sets of players. The matching problem consists in assigning resources from one set to the players of the other set. To do so, each player uses a preference relation based on a computed metric to rank the players of the other set from the viewpoint of a target utility function. By maximizing its utility, each player maximizes its objective function.

In the context of this paper, the resource allocation problem consists in assigning the network resources, i.e., frequency subband/power level pairs, to the IoT devices given a set of constraints from both the devices and the resources. It can therefore be modeled as a one-to-one matching game, meaning that each device can connect to one resource and each resource can be matched to one IoT device. We consider that the system is fully coordinated by the BS and that the BS knows all devices and relay power budgets, relay buffers status and device traffic types. With this information and the centralized ability to devise a suitable resource allocation, this coordinated matching-based approach is expected to provide the best results in terms of QoS requirements. Being power limited, the IoT devices rank the subband/power level pairs according to their channel gain with the relay. In fact, the subband/power level pair requiring the least amount of transmit power is preferred by the device. In contrast, resource ranking mainly depends on the status of the corresponding relay buffer and the conditions between the relay and the BS. The fuller the relay buffer, the less likely it is to accept devices. The same applies to the channel gain between relay and BS: the better the radio conditions, the more likely the relay will accept devices.

NUMERICAL RESULTS

We consider an uplink communication system with one BS, one relay and 60 IoT devices divided into 36 devices generating continuous traffic, 12 devices having a sporadic traffic type with a spectrum access probability equal to 0.3 and 12 devices with a periodic traffic type. The periodic traffic devices are, in turn, divided into 4 groups of 3 devices and a time offset is applied for network access of the devices in the different groups to avoid mutual repeated collisions.

The target data rate is 0.75 Mb/s for the continuous traffic devices and 0.3 Mb/s for the ones with non-continuous traffic. The number of available subbands is 9. Five power levels are available at each subband, hence the total number of available resources is 45. This allows us to address the problem of congestion on a small scale, first. The objective is to prove the superiority of the proposal in supporting more linked devices than a conventional system while using the same finite number of resources. Then, by adding network nodes, the proposal can be scaled up to enable massive connectivity. Signals undergo a distance-dependent path-loss with a decay factor of 3.76, and a zero-mean log-normal shadowing with an 8 dB variance. The noise power spectral density is equal to 4 × 10 -21 W/Hz. The power budget of the IoT devices is equal to 0.2 W and the one of the relays is 5 W. The congestion threshold of the relays is set to 50 percent of their whole buffer capacity (300 packets) and the duration of the timeout period is taken equal to 50 time slots.

Benchmarks

The performance of the proposed algorithm with and without AoI prediction is evaluated for 1000 time slots using Matlab and compared with the fully coordinated approach based on matching theory and with two more basic uncoordinated solutions:

• random access NOMA (RA-NOMA) method, where devices choose the frequency subband/power level pair randomly, • an allocation technique (No-AoI) where the devices solely allocate their arms using UCB without taking the AoI effect into consideration.

Evaluation

The efficiency of the arm assignment process is evaluated in Fig. (4a), where the instantaneous successful transmission percentage, i.e., the proportion of devices having their signals correctly decoded by the relay, is plotted as a function of time. We can observe that the basic RA-NOMA approach cannot achieve a success rate above 40 percent with an average of less than 30 percent. In contrast, since the matching technique is fully coordinated by the BS, convergence towards the full satisfaction state is reached in the very first time slots. Similarly, the proposed UCB-based technique also reaches convergence with a 100 percent satisfaction rate for all IoT devices after an exploration phase of 80 time slots.

The average transmit power per device as a function of the requested rate is evaluated in Fig. (4b). Simulations showed that a balanced scheme with β = 0.5 is the best choice for minimizing power consumption. Moreover, for the same average transmit power, the proposed UCB-based method with β = 0.5 is the only solution within 0.15 Mb/s of the fully coordinated benchmark solution for the same average transmit power.

Figures (4c) and (4d) show the average AoI value when varying the total number of sporadic players or the access probability. In Fig. (4d), the number of sporadic devices is 9. The proposed approach significantly outperforms basic random access approaches.

REMAINING CHALLENGES

The next step prior to any wide adoption resides in a thorough investigation of system robustness to imperfect or partial channel state information (CSI). Following this step, it would be possible to evaluate the benefits of the reduced signaling overhead. Additional impactful fine-tuning includes the trade-off between exploration and exploitation phases and higher-level auto-configuration capabilities. The latter involves devising adaptive methods for unsupervised traffic rerouting to avoid relay failures, or even, intelligent automated allocations to improve system spectral and/or power efficiencies. Intelligent reflective surfaces and on-demand mobile relays can represent appealing complementary means to these same effects and may be considered for future developments. On a different note, some physical layer techniques, such as waveform design, are considered essential enablers for uncoordinated spectrum access and should be thoroughly investigated. Indeed, robustness to impairments including carrier frequency offsets, timing offsets and efficient multi-user management are key for unleashing the potentials of such decentralized intelligent allocation methods.

CONCLUSION

In this article, the uncoordinated decentralized joint channel and power allocation problem was addressed and formulated using the MAB framework. A solution based on the UCB algorithm was then proposed. In addition to increased spectral efficiency, the proposed solutions minimize AoI for IoT devices in a two-hop mixed-traffic NOMA uplink communication system with relays. Thanks to their self-configuration capabilities, they are able to address relay-level congestion, provide the scalability for mixed traffic support and manage AoI for the coexistence of multiple services and applications. This was achieved through the UCB algorithm augmented by both the introduction of buffering and an AoI-controlling scheduler at the level of the relays. AoI prediction for sporadic traffic was shown to further improve performance for delay-sensitive applications. The resulting self-configurable and scalable proposal was shown to outperform competing uncoordinated techniques while approaching the performance of a centralized solution with perfect knowledge of all system parameters. Future work should focus on key enablers for such a system including imperfect or partial CSI and physical layer support.

Fig. 1 :

 1 Fig. 1: Illustration example of uplink NOMA principle with four power levels per subband. The devices' signals are superimposed in the power domain on the frequency resources at each time slot. Each device is associated with a power level on one subband and the height of each box represents the value of the received power.

2 .Fig. 2 :

 22 Fig. 2: Uplink communication system model with one BS and two relays, with a certain number of devices having different traffic nature.

Fig. 4 :

 4 Fig. 4: Algorithm evaluation in terms of success rate, transmit power and AoI in function of time slots. (4a) Instantaneous percentage of successfully connected devices. (4b) Average transmit power for a varying requested rate. (4c) Average AoI for a varying number of sporadic traffic devices. (4d) Average AoI for a varying access probability of sporadic traffic devices.

TABLE I :

 I System features

	Feature	Type or nature
	Type of system	Uplink with uncoordinated spec-
		trum access and NOMA
	Traffic	Mixed: continuous, periodic and
		sporadic
	Links	Relay uplink: IoT device → relay
		Backhaul link: relay → BS
	Relay mode	Full duplex
	Channel state infor-	Perfectly estimated at the receiver
	mation
	Optimization prob-	Subband/power level pair allocation
	lem
	Objective function	Traffic type-dependent device satis-
		faction with best possible spectral
		efficiency
	THE RESOURCE ALLOCATION PROBLEM

  Joseph Doumit received the B.E. and M.S. degrees in computer and communications engineering from Lebanese University, Lebanon, in 2020. He is currently pursuing a Ph.D. degree in information and communication engineering with IMT Atlantique, France. His current research interests include resource allocation, non-orthogonal multiple access, reconfigurable intelligent surfaces, and reinforcement learning.Marie-Josepha Youssef received the B.E. and M.S. degrees in computer and communications engineering from Lebanese University in 2016, and the Ph.D. degree in information and communication engineering from the IMT Atlantique, France, in 2020. She is currently a quantitative researcher at Capital Fund Management. Her current research interests include resource allocation, non-orthogonal multiple access, matching theory, reinforcement learning, unmanned-aerial vehicles, and grant-free communications.Charbel AbdelNour obtained his PhD from Telecom Bretagne in 2008 and his Habilitation thesis from the University of South Brittany in 2020. Since Dec. 2022, Charbel holds a Professor position at the Mathematical and Electrical Engineering Department at IMT Atlantique. His interests concern coded modulations, resource allocation for NOMA, waveform design, MIMO and iterative receivers. Lately, he presented several contributions to the H2020 METIS, FANTASTIC5G and EPIC projects and to the 3GPP consortium related to coding solutions for 5G. He also served as co-chair of the technical program committee of ISTC (International Symposium on Topics in Coding) in 2023. Joumana Farah received a B.E. in Electrical Engineering from the Lebanese University, in 1998, the M.E. and Ph.D. degrees in 1999 and 2002 from Grenoble INP, France. In 2010, she obtained the Habilitation to Direct Research from Université Paris VI. Until 2022, she was a Professor at the Lebanese University. Since January 2023, she has joined the Institut National des Sciences Appliquées (INSA) at Rennes, France. She has nine registered patents and has coauthored over a hundred journal and conference papers. Her current research interests include resource allocation techniques, channel coding and interference management techniques. Catherine Douillard is a full Professor at IMT Atlantique/Lab-STICC laboratory, Brest, France. Her main research interests are error correcting codes, coded modulations and diversity techniques for multi-carrier, multiantenna and multiple access transmission systems. Between 2007 and 2012, she participated in Digital Video Broadcasting Technical Modules for the definition of DVB-T2, DVB-NGH and DVB-RCS2 standards. She served several times as chair of the technical program committee and general chair of ISTC (International Symposium on Topics in Coding). In 2009, she received the SEE/IEEE Glavieux Award for her contribution to standards and related industrial impact.

ACKNOWLEDGMENT

This work was carried out in the context of Beyond5G, a project funded by the French government as part of the economic recovery plan and the investments for the future program.