Zhan Yu

Yue Sun

Minh-Ngoc Vu

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

Modelling of mixed cracks in rock-like brittle materials under compressive stresses by a double phase-field method

Keywords: Damage, cracking, crack coalescence, phase-field method, rock-like materials, Compressive stresses Highlights

A new phase-field method is proposed in this paper for modelling complex cracking process in rock-like brittle materials under compression-dominating stresses. For this purpose, two crack fields are introduced in order to respectively describe tensile and shear cracks, related to two independent dissipation processes. The evolution of these fields is described by using the variational principle and thermodynamics framework. The proposed method is implemented in a finite element code and applied to investigating cracking process in a rocklike sample containing two initial flaws and subjected to uniaxial and bi-axial compression.

Both tensile wing and shear cracks as well as crack coalescence observed in laboratory tests are successfully reproduced by the proposed method. Differently with most previous studies, the effect of confining stress on shear cracking growth is well described.

Introduction

In many engineering materials, macroscopic fractures are generally generated by the coalescence of diffuse micro-cracks. The description of transition from diffuse micro-cracks to localized macroscopic fractures is the key issue of modelling failure process in materials and structures. Further, under multi-axial loading conditions, multiple fractures can initiate, propagate and interact. There is a real challenge for numerical modelling. Different methods have been developed during the last decades to deal with this complex subject. For example, enriched finite element methods have been first developed to account for displacement discontinuities along crack surfaces at the elementary level with the help of enriched shape functions [START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations, part 1: fundamentales[END_REF]. Extended finite element methods (XFEM) have been developed to deal with crack displacement discontinuities by using nodal enrichment techniques [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF]. These methods have successfully been applied to different types of materials and structures, in particular hydraulic fracturing in porous rocks [START_REF] Zeng | Numerical study of hydraulic fracture propagation accounting for rock anisotropy[END_REF][START_REF] Zeng | Study of hydraulic fracturing in an anisotropic poroelastic medium via a hybrid edfm-xfem approach[END_REF]. In these two types of methods, the transition from diffuse damage to localized fractures is still a pending issue. The thick level set has been recently proposed to remediate the issue of damage-fracturing transition [START_REF] Moes | A level set based model for damage growth: The thick level set approach[END_REF][START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF].

The description of multiple cracks in three-dimensional conditions is also a delicate task.

On the other hand, based on homogenization techniques, some micro-mechanical models have also been developed for capturing the initiation and growth of micro-cracks and their transition to localized cracks in a representative elementary volume [START_REF] Zhu | Analytical and numerical analysis of frictional damage in quasi brittle materials[END_REF][START_REF] Zhao | A micromechanics-based plastic damage model for quasi brittle materials under a large range of compressive stress[END_REF][START_REF] Zhao | Analysis of localized cracking in quasi-brittle materials with a micromechanics based friction damage approach[END_REF]. The efficiency of such models at the structure scale still needs to be demonstrated.

More recently, based on the framework of variational principle for fracture mechanics [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and the optimal approximations methods of functionals with jumps [START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via t-convergence[END_REF][START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], the so-called phase-field method has been rapidly developed [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. The sharp topology of crack is approximated by a regularized one which is a non-local function of damage variable and its gradient (phase-field). The phase-field is determined by a specific boundary values problem.

This approach is particularly suitable to deal with the natural transition from diffuse damage to localized cracks, interacting multiple cracks and three dimensional problems. Different extensions have been proposed for various purposes, dynamic brittle fracture [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF], fast computing with a new hybrid formulation [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF], multi-physics problems [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-elasticplastic solids[END_REF], finite deformation [START_REF] Borden | A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects[END_REF], and coupling between damage and plasticity [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF][START_REF] Fang | Phase field fracture in elasto-plastic solids: variational formulation for multi-surface plasticity and effects of plastic yield surfaces and hardening[END_REF], etc.. In most previous studies, the accent was put on tensile cracks which are driven by the elastic strain energy related to tensile eigenvalues. As in a wide range of engineering structures, shear cracks are also an important mechanism, different extensions have been proposed in some recent works. For instance in [START_REF] Zhang | A modification of the phase-field model for mixed mode crack propagation in rock-like materials[END_REF][START_REF] Bryant | A mixed-mode phase field fracture model in anisotropic rocks with consistent kinematics[END_REF], the authors have proposed a specific decomposition of the ratio between driving strain energy and critical fracture energy into a tensile part (mode I) and a shear part (mode II). The phase-field evolution is thus driven by both the tensile and shear strain energy. And each part is controlled by a critical fracture energy (mode I and mode II). The efficiency of this method in modelling mixed crack propagation has been demonstrated. More recently, in [START_REF] Nguyen | Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials[END_REF], the concept of multiple phasefield has been introduced to describe the damage accumulation in different cleavage planes of polycrystalline materials. Similarly, this type of concept has been applied to anisotropic materials in order to capture different cracking patterns in [START_REF] Bleyer | Phase-field modeling of anisotropic brittle fracture including several damage mechanisms[END_REF]. However, those previous studies have mainly focused on tensile and shear loads.

Most rock-like materials are subjected to compression-dominating stresses. Crack initiation and propagation in these materials subjected to compressive loads exhibit complex and mixed patterns [START_REF] Evans | The brittle-ductile transition in rocks: Recent experimental and theoretical progress, The Brittle-Ductile Transition in Rocks[END_REF][START_REF] Wong | The brittle-ductile transition in porous rock: A review[END_REF], such as compression-induced tensile wing and shear cracks. Further, the mechanical behavior of rock-like materials is strongly dependent on compressive mean stress. For instance, the shear crack growth is directly affected by the compressive normal stress. For this purpose, a phase-field model for compressive-shear fracture was proposed by introducing a hybrid approach in [START_REF] Zhou | Phase field modeling of brittle compressive-shear fractures in rock-like materials: A new driving force and a hybrid formulation[END_REF]. Moreover, based on the assumption of frictional contact [START_REF] Fei | A phase-field method for modeling cracks with frictional contact[END_REF] and friction energy [START_REF] Palmer | The growth of slip surfaces in the progressive failure of over-consolidated clay[END_REF], some authors proposed a frictional fracture phase-field model [START_REF] Fei | A phase-field model of frictional shear fracture in geologic materials[END_REF] and a mixed-FE scheme [START_REF] Wang | A stabilized mixed-fe scheme for frictional contact and shear failure analyses in deformable fractured media[END_REF] to deal with compression-induced shear cracks in rock-like materials. On the other hand, the phase-field model for mixed-mode fracture were introduced in [START_REF] Spetz | A modified phase-field fracture model for simulation of mixed mode brittle fractures and compressive cracks in porous rock[END_REF][START_REF] Wang | A phase-field model for mixed-mode fracture based on a unified tensile fracture criterion[END_REF]. Further, some authors focused on the compressive-ductile behavior of mixed-mode phase-field model by considering micro-structure [START_REF] Ulloa | A micromechanicsbased variational phase-field model for fracture in geomaterials with brittle-tensile and compressiveductile behavior[END_REF][START_REF] You | Brittle-ductile failure transition in geomaterials modeled by a modified phase-field method with a varying damage-driving energy coefficient[END_REF].

The so-called multi phase-field were firstly proposed and applied on polycrystal [START_REF] Oshima | Development of multi-phase-field crack model for crack propagation in polycrystal[END_REF].

The conception of "multi" is introduced to describe the anisotropy of damage on cleavage planes. The same conception is also used in the most work of multi phase-field [START_REF] Nguyen | Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials[END_REF][START_REF] Na | Computational thermomechanics of crystalline rock, part i: A combined multi-phasefield/crystal plasticity approach for single crystal simulations[END_REF][START_REF] Dean | A multi phase-field fracture model for long fiber reinforced composites based on the puck theory of failure[END_REF].

Conversely, it is worth to note that the conception of the double phase-field is different from that of the multi phase-field, even if they have a similar call. The principal feature of double phase-field method is the use of two damage variables (or more than two if necessary). This conception is firstly appeared in the study of [START_REF] Bleyer | Phase-field modeling of anisotropic brittle fracture including several damage mechanisms[END_REF], the two damage variables is introduced to identify the different effect of longitudinal and transverse failure for anisotropic materials.

For a different purpose, [START_REF] Yu | Numerical study of thermo-hydro-mechanical responses of in situ heating test with phase-field model[END_REF] proposed a double phase-field with the damage variables defined as tensile damage and shear damage, which are driven by the energy decomposed according to the strain state. Similarly, [START_REF] Fei | Double-phase-field formulation for mixed-mode fracture in rocks[END_REF] proposed a double phase-field method to consider the tensile and compressive-shear crack. The onset/effect of tensile crack and compressive-shear crack are according to crack conditions of "open", "stick" and "slip". The two types of crack are identified based on a complex calculation due to the stress in the direction perpendicular and parallel to the interface of the crack.

In the present study, a new phase-field method is proposed for modelling mixed cracks in rock-like brittle materials subjected to compression-dominating stresses. Two phase-field variables are introduced to describe the tensile and compressive-shear cracks respectively.

Mixed cracks can be represented by the combination of these two phase-field variables by using a spectral decomposition of stress. Different from the work of [START_REF] Yu | Numerical study of thermo-hydro-mechanical responses of in situ heating test with phase-field model[END_REF], the tensile and compressive crack are distinguished by stat of stress instead of strain in order to adopt the rock-like material. Furthermore, with a same purpose of [START_REF] Fei | Double-phase-field formulation for mixed-mode fracture in rocks[END_REF], by using a more simply and also reasonable way, a modified Morh-Coulomb criterion is introduced to identify the compressive-shear crack, which is more suitable to consider the properties of internal friction and cohesion for rock-like material. Another way to deal with the issue of mixed crack is provided. The efficiency of the proposed method is assessed through comparisons between numerical results and experimental evidence. It is shown that different types of crack propagation scenarios can be well described by the proposed method.

A double phase-field method

Regularized crack density

We consider here a solid body occupying the volume Ω, subjected to the body force f b in Ω, the surface force t N on its external boundary ∂Ω f and the prescribed displacement u on its external boundary ∂Ω u . Due to the applied loads, there is initiation and propagation of cracks. In the case of rock like materials under compression-dominating stresses, complex cracking processes have been observed [START_REF] Wong | Micromechanics of faulting in westerly granite[END_REF][START_REF] Wong | Analysis of crack coalescence in rock-like materials containing three flaws-part i: experimental approach[END_REF][START_REF] Paterson | Experimental rock deformation-the brittle field[END_REF]. In general, three types of cracks, namely tensile, shear and mixed, can initiate and propagate depending on loading path and material microstructure. Under direct tension or in triaxial extension conditions, tensile cracks represent the main process. Under triaxial compression conditions, in particular when confining stress is high, frictional shear cracks are the dominant mechanism. In general loading conditions, mixed cracks can be generated by combined tensile and shear strains or stresses. In order to conveniently describe such complex cracking processes, two independent crack fields are here introduced to respectively represent the tensile and shear cracks. Mixed cracks are then described by the combination of these two fields.

In the framework of phase-field method, the real crack surface area is approximated by the regularized one [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. By considering the two crack fields, the total regularized crack surface area Γ l d is expressed as follows:

Γ l d = Γ t l d (d t ) + Γ s l d (d s ) = ∫ Ω {γ t (d t , ∇d t ) + γ s (d s , ∇d s )}dV (1)
where Γ α l d (α = t, s) are the contributions of tensile and shear cracks to the total crack surface area. Two scalar-valued functions γ α (d α , ∇d α ) denote the tensile and shear crack density (surface area per unit volume). There are different forms available to define the crack surface density γ α . A common form was introduced in [START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via t-convergence[END_REF] and it is adopted here:

γ α (d α , ∇d α ) = 1 2l d (d α ) 2 + l d 2 ∇d α .∇d α ; α = t, s (2)
In these functions, a scalar-valued variable d α (x) (α = t, s) is introduced for each type of crack and it is conventionally called the phase-field (or crack) variable. It is worth noticing that the crack surface density γ α is a function of d α and its gradient ∇d α . l d is a length scale parameter controlling the width of regularized smeared crack surfaces. This nonlocal property allows the phase-field method to avoid the pathologic mesh dependency of numerical solution. Further, the value of Further, according to the Griffith theory of linear fracture mechanics, the energy needed for the creation of a crack surface area is characterized by the material toughness or fracture energy. With the regularized smooth representation of cracks adopted here, one can express the energy density per unit volume needed for the creation of tensile and shear cracks as follows:

d α (x) (α = t,
w c (d t , d s , ∇d t , ∇d s ) = g t c γ t (d t , ∇d t ) + g s c γ s (d s , ∇d s ) (3)
g t c and g s c denotes the material toughness for the tensile and shear crack respectively.

Variational framework

According to the variational approach to fracture mechanics of elastic materials proposed in [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], under the isothermal conditions, the total energy functional of a solid domain Ω is composed of the elastic strain energy and that requested to create the cracks, with the regularized crack fields adopted above, one gets:

E(ε, d t , d s ) = ∫ Ω w e (ε, d t , d s )dV + ∫ Ω w c dV ( 4
)
w e is the elastic strain energy density of cracked material. ϵ(u) denotes the second order linear strain tensor with u(x) being the displacement vector such as ϵ = ∇ s u.

Elastic strain energy

In the present study, two crack fields are introduced and affect differently the elastic strain energy. In order to better describe the coupling between the cracking processes and elastic properties of cracked material, the elastic strain energy of undamaged solid is decomposed into a tensile part w 0 e+ + and a compressive part w 0 e-:

w 0 e (ϵ) = w 0 e+ (ϵ) + w 0 e-(ϵ) = 1 2 σ + : ϵ + 1 2 σ -: ϵ (5)
Two symmetric tensors σ ± denote the tensile and compressive parts of the Cauchy stress tensor σ, and they are calculated by using the following fourth order projection operators

P ± σ :    σ + = P + σ : σ σ -= P - σ : σ (6)
These operators are further constructed from the spectral decomposition of the stress tensor proposed in [START_REF] Lubarda | Damage model for brittle elastic solids with unequal tensile and compressive strengths[END_REF][START_REF] Murakami | Continuum damage mechanics: a continuum mechanics approach to the analysis of damage and fracture[END_REF]:

   P + σ = ∑ 3 β=1 ⟨σ β ⟩ + n β ⊗ n β ⊗ n β ⊗ n β P - σ = I -P + σ ( 7
)
I denotes the fourth-order identity tensor. n β (β = 1, 2, 3) define the three orthogonal principal directions of the stress tensor while σ β are the three principal stresses. The bracket ⟨.⟩ + is defined as:

   ⟨a⟩ + = a, a ≥ 0 ⟨a⟩ + = 0, a < 0 (8)
It is worth noticing that the projection operators are widely used for the decomposition of stress or strain tensor into tensile and compressive parts, including some previous phasefield models [START_REF] Miehe | Comparison of two algorithms for the computation of fourth-order isotropic tensor functions[END_REF][START_REF] Miehe | Algorithms for computation of stresses and elasticity moduli in terms of seth-hill's family of generalized strain tensors[END_REF][START_REF] Zhang | A new phase field fracture model for brittle materials that accounts for elastic anisotropy[END_REF]. In general, even for elastically isotropic materials, the projection operators for stress P σ and for strain P ϵ are different. As in the previous work [START_REF] Zhang | A new phase field fracture model for brittle materials that accounts for elastic anisotropy[END_REF], the decomposition of stress tensor is adopted here. Compared with the strain tensor decomposition as that used in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], the stress decomposition avoids the implication of elastic parameters.

It is particularly convenient for future extension to initial anisotropic materials.

The elastic properties of cracked material are affected by the growth of cracks. This is generally described by introducing different forms of degradation functions. In the present study, two degradation functions are defined respectively for the tensile and shear cracks.

The elastic strain energy of cracked material is expressed as follows:

w e (ϵ, d t , d s ) = h t (d t )w 0 e+ + h s (d s )w 0 e- (9)
The common form widely used in previous studies [START_REF] Miehe | Comparison of two algorithms for the computation of fourth-order isotropic tensor functions[END_REF][START_REF] Miehe | Algorithms for computation of stresses and elasticity moduli in terms of seth-hill's family of generalized strain tensors[END_REF][START_REF] Zhang | A new phase field fracture model for brittle materials that accounts for elastic anisotropy[END_REF]] is adopted here for the two degradation functions:

h α (d α ) = (1 -d α ) 2 (1 -k) + k; α = t, s (10)
k is a small positive value that is used to ensure the positive definiteness of elastic stiffness tensor after the material is fully broken (d α = 1).

Governing equations of crack fields

Based on the expressions of w e and w c defined above, the total energy functional E can be detailed in the following form:

E = ∫ Ω w e (ε, d t , d s )dV + ∫ Ω w c (d t , d s , ∇d t , ∇d s )dV = ∫ Ω { h t (d t )w 0 e+ + h s (d s )w 0 e- } dV + ∫ Ω { g t c γ t (d t , ∇d t ) + g s c γ s (d s , ∇d s ) } dV (11)
We assume that the growths of tensile and shear cracks are two independent dissipation processes. Each of them verifies the positiveness condition independently. Further, inspired by the previous studies [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] Wu | Chapiter onephase field modeling of fracture[END_REF], each crack field can be determined by solving the minimization problem of the energy functional E. By calculating the first order variation of E and using the unilateral stationary condition for the total energy functional with respect to each crack variable, namely δE = 0 for δd α > 0 and δE > 0 for δd α = 0 (α = t, s), the following governing equations for the evolution of two crack fields (α = t, s) can be derived:

         -∂we ∂d α -g α c δ d α γ α = 0 , ḋα > 0 , in Ω -∂we ∂d α -g α c δ d α γ α ≤ 0 , ḋα = 0 , in Ω ∂γ ∂∇d α .n = 0 , on δΩ (12)
where

δ d α γ α = ∂γ α ∂d α -div ( ∂γ α ∂∇d α
) denotes the variational derivative of crack field d α [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] Wu | Chapiter onephase field modeling of fracture[END_REF].

By substituting w e and γ α by their specific forms presented above, one obtains the governing equations for two crack fields as follows:

   -h ′ t (d t )w 0 e+ -g t c { 1 l d d t -l d div(∇d t ) } = 0 , ḋt > 0 -h ′ s (d s )w 0 e--g s c { 1 l d d s -l d div(∇d s ) } = 0 , ḋs > 0 (13)
where

h ′ t (d t ) = -2(1 -k)(1 -d t ) and h ′ s (d s ) = -2(1 -k)(1 -d s )
are the derivatives of two degradation functions with respect to two damage variables.

From the threshold function (13)(b), it is observed that the evolution of the shear crack d s is driven by the compressive part of the elastic strain energy associated with the intact state w 0 e-. This is physically not fully justified. In practice, for most rock-like materials under compressive stresses, the shear crack is physically controlled by a generalized shear stress. The later can be well represented by the shear force defined by the widely used Mohr-Coulomb criterion. Based on this physical background and inspired by some previous studies [START_REF] Li | A damage model for hard rock under stress-induced failure mode[END_REF][START_REF] Zhou | Phase field modeling of brittle compressive-shear fractures in rock-like materials: A new driving force and a hybrid formulation[END_REF], an alternative driving energy w s is defined as follows for the description of shear crack growth:

w s (σ) = 1 2G ⟨ σ 1 -σ 3 2cosφ + σ 1 + σ 3 2 tanφ -c ⟩ 2 + (14)
In this relation, σ 1 and σ 3 are the major and minor principal stress respectively. G is the initial shear modulus of intact material. φ denotes the internal friction angle and c the cohesion. Moreover, it is assumed that the shear crack growth occurs only in the compressive domain of principal stresses. Therefore, the driving energy of shear crack growth corresponds to the compressive part of w s , which is defined as:

w s -(σ -) = 1 2G ⟨ ⟨σ 1 ⟩ --⟨σ 3 ⟩ - 2cosφ + ⟨σ 1 ⟩ -+ ⟨σ 3 ⟩ - 2 tanφ -c ⟩ 2 + ( 15
)
with the bracket ⟨.⟩ -defining:

   ⟨a⟩ -= 0, a ≥ 0 ⟨a⟩ -= a, a < 0 (16)
By substituting the physically-based driving energy for the compressive shear crack growth for [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], the new criteria for two crack fields is expressed as:

   -h ′ t (d t )w 0 e+ -g t c { 1 l d d t -l d div(∇d t ) } = 0 , ḋt > 0 -h ′ s (d s )w s --g s c { 1 l d d s -l d div(∇d s ) } = 0 , ḋs > 0 (17)
On the other hand, in order to describe the irreversible process of crack growth in general loading paths with unloading and reloading cycles, the concept of loading history functional has been introduced in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field fe implementations[END_REF]. Based on that concept, two energy history functionals are here defined as follows for the grow description of tensile and shear cracks:

   H t = max[w 0 e+ (σ + , τ )] τ ∈[0,t] H s -= max[w s -(σ -, τ )] τ ∈[0,t] (18)
The time variable τ designates the loading history from the initial state to the current time step t. Using these energy history functionals, the evolution of the tensile and shear crack is determined by the following criteria:

   -h ′ t (d t )H t -g t c { 1 l d d t -l d div(∇d t ) } = 0 -h ′ s (d s )H s --g s c { 1 l d d s -l d div(∇d s ) } = 0 (19)

Numerical implementation with finite element method

The problem to be solved here consists in the determination of displacement field u(x), the tensile crack field d t (x) and the shear crack field d s (x), at each loading step and with the boundary conditions defined above. For this purpose, the total potential energy is expressed as:

Π = E(u, d t , d s ) - ∫ Ω f b • udΩ - ∫ ∂Ω f t N • udA (20)
The stationarity condition of the potential energy δΠ = 0 leads to three coupled systems of equations to be solved in terms of the displacement, the tensile crack and the shear crack fields respectively. However, due to the growth of cracks, it is needed to solve a strongly non-linear problem. For this purpose, the whole loading path is divided into a number of increments or time steps. Starting from the initial known values of all physical quantities, at the end of loading step i, the values of u(x, t i ), d t (x, t i ) and d s (x, t i ) are assumed to be known.

During the loading step t t+1 = t i + ∆t t+1 , the solutions to be determined are the increments of three principal unknowns, namely ∆u(x, ∆t t+1 ), ∆d t (x, ∆t t+1 ) and ∆d s (x, ∆t t+1 ).

The finite element method is adopted in the present work. Based on the geometrical discretization and with the help of the shape functions and their derivatives for the type of element adopted, one obtains three discrete systems of equations to be solved.

For the mechanical problem, the Cauchy stress tensor σ verifies the following equilibrium equations:

   ∇ • σ + f b = 0 in Ω σ • n = t N on ∂Ω t (21)
The constitutive relation of cracked material is defined by:

σ = C(d t , d s ) : ϵ ; C(d t , d s ) = h t (d t )P + σ : C 0 + h s (d s )P - σ : C 0 ( 22
)
where C 0 is the elastic stiffness matrix of intact material. Together with the kinematic relation ϵ = ∇ s u, the discrete system of equations of the mechanical problem is given by:

         K u △U = △F K u = ∫ Ω B T u CB u dΩ △F = ∫ Ω N T u △f b dΩ + ∫ ∂Ω f N T u △t N dA (23)
N u and B u are respectively the matrix of shape functions and related derivatives of displacement components. K u is the global stiffness matrix. △U and △F denote the column matrices of increment nodal displacements and forces during the current loading step (for the sake of simplicity, the step index i + 1 is here omitted).

On the other hand, with the crack growth criteria [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF] in hand, the tensile crack phasefield is determined by the following local equations:

               h ′ t (d t )H t + g t c { 1 l d d t -l d div(∇d t ) } = 0 in Ω ḋt ≥ 0 in Ω d t = 1 on Γ t ∇d t • n = 0 on ∂Ω ( 24
)
while the shear crack phase-field is described by the following ones:

               h ′ s (d s )H s -+ g s c { 1 l d d s -l d div(∇d s ) } = 0 in Ω ḋs ≥ 0 in Ω d s = 1 on Γ s ∇d s • n = 0 on ∂Ω (25)
The same geometrical mesh is used for the phase-field problems as for the mechanical one. Further, as for the displacement field, the two crack fields and their gradients inside each element are also approximated in terms of by the nodal values by using appropriate shape functions and related derivatives. With such approximations, the following discrete systems of equations are obtained for each crack phase-field. Namely, for the tensile cracks field, one gets:

         K d t d t = F d t K d t = ∫ Ω {(g t c /l d + 2H t )N T d N d + g t c l d B T d B d }dΩ F d t = ∫ Ω 2H t N T d dΩ ( 26
)
and for the shear crack field:

         K d s d s = F d s K d s = ∫ Ω {(g s c /l d + 2H s -)N T d N d + g s c l d B T d B d }dΩ F d s = ∫ Ω 2H s -N T d dΩ (27)
N d and B d are the matrices of shape functions and related derivatives for the approximation of two crack phase-fields. K d α (α = t, s) is the global stiffness matrix related to the tensile or shear crack field. F d α denotes the column matrix of nodal forces while d α that of nodal damage variables at the current loading step (for the sake of simplicity, the step index i + 1 is here omitted). It is worth noticing that with the help of two energy history functionals defined in [START_REF] Borden | A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects[END_REF], the values of damage variables at the current loading step are directly determined as functions of the accumulated driving energy for each cracking mechanisms.

It is obvious that three fields are strongly coupled. The calculation of displacement field is influenced by the crack phase-fields due to the fact that the elastic stiffness matrix of cracked material C(x) is a function of two crack fields d t (x) and d s (x). Inversely, the growth of both tensile and shear cracks is controlled by the mechanical energy which is related to the variation of displacement field.

In this way, the so-called Alternate Minimization (AM) solver proposed in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF] is used in this work because of its good robustness performance. At the iteration k + 1 of the loading step n + 1, two solution stages are summarized as:

Stage 1:

K u (C)△U k+1 n+1 = △F n+1 , with C = C([d t , d s ] k n+1 ) (28)
-Solving the displacement △U k+1 n+1 by using the damage values fixed [d t , d s ] k n+1 obtained from the last iteration k.

Stage 2:

   K d t (H t )[d t ] k+1 n+1 = F d t (H t ), with H t = H t (U k+1 n+1 ) K d s (H s -)[d s ] k+1 n+1 = F d s (H s -), with H s -= H s -(U k+1 n+1 ) (29)
-Solving the damage values [d t , d s ] k+1 n+1 by using the updated displacement U k+1 n+1 obtained from the same iteration.

The two steps are repeated until the condition ∥ max([

d t ] k+1 n+1 -[d t ] k n+1 , [d s ] k+1 n+1 -[d s ] k n+1 ) ∥≤ R
is satisfied, with R = 1 × 10 -5 being used in this work.

Experimental verification

In this section, the efficiency of the proposed double phase-field method, in particular its capability to describe complex cracking patterns under compression-dominating loads is verified, through the comparison between numerical results and experimental observations.

Uniaxial compression of a plate with two inclined flaws

In this example, a representative laboratory test is investigated. Uniaxial compression is performed on a rectangular thin plate containing two initial inclined flaws, as illustrated in Figure 2. Different configurations of the initial flaws can be considered. More details about the testing procedures can be found in [START_REF] Bobet | Fracture coalescence in rock materials: Experimental observations and numerical predictions[END_REF][START_REF] Bobet | Numerical modeling of fracture coalescence in a model rock material[END_REF]. This type of test is particularly relevant for the characterization of complex failure process in brittle materials, including tensile, shear and mixed cracks and their bridging. The proposed double phase-field method is now applied to capturing such cracking patterns.

The geometry of the plate and the prescribed boundary conditions are presented in Figure 2. The length of two initial flaws is 2a=12.7mm and their width is 0.1mm. The configuration of two flaws is defined in terms of the combination "flaw angle-spacing-continuity". For example, the case of "45-a-2a" corresponds to the configuration of two flaws inclined with an angle of α = 45 • , spaced by s = a and c = 2a. In this study, three configurations are considered, namely "45-0-2a", "45-a-2a" and "45-2a-2a" are selected for simulation. The values of mechanical parameters involved on the proposed model are selected mainly from previous studies [START_REF] Bobet | Fracture coalescence in rock materials: Experimental observations and numerical predictions[END_REF][START_REF] Bobet | Numerical modeling of fracture coalescence in a model rock material[END_REF]. The set of values is given in Table 1. Therefor, l d is assumed as a pure numerical parameter of the regularized model of brittle fracture which depends on element size h: h < l d /2 [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field fe implementations[END_REF]. In this work, we define l d = 2.5h for both of the tensile and shear damage field. h r is the size of the refined element which are assigned to the region expected to have the propagation of cracks for the sake of computational efficiency.

Following the study of [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF], g t c /l d is considered as a real material parameter, which depends on stiffness and strength:

σ c = 3 16 √ 3 2 √ Eg t c l d .
σ c (MPa) 3.2 φ ( • ) 30
c (MPa) 1
g t c (kN/mm) 4.2 × 10 -6 g s c (kN/mm) 21 × 10 -6 l d (mm) 0.25 k (-) 1 × 10 -9
Mesh:

h r (mm) 0.1
Element number Around 300,000

Node number Around 300,000

Displacement increment:

∆u (mm) 1 × 10 -3
In Figure 3, one shows the distributions of tensile and shear damage (d t and d s ) for the case of "45-0-2a" at six subsequent values of prescribed axial displacement. The red color shows the tensile crack and the blue color shows the shear one. In order to investigate the cracking patterns more clearly, only the crack with the damage value higher than 0.5 is shown. Further, for the purpose of comparison, the main crack patterns observed in the laboratory test and reported in [START_REF] Bobet | Fracture coalescence in rock materials: Experimental observations and numerical predictions[END_REF][START_REF] Bobet | Numerical modeling of fracture coalescence in a model rock material[END_REF] In order to compare with the stress value measured from experiment and identify the stat of displacement loading in the simulation, the stress-displacement curve is given in Figure 6. Furthermore, a quantitative comparison between numerical and experimental results is provided for all of the three tests. The step (b) of Figure 3, 4 and 5 shows the moment of wing crack's onset, and the step (e) shows the one of coalescence crack's onset. Good agreement is achieved from the quantitative comparison for both of wing crack and coalescence crack.

Moreover, it is clear to see that the onset of wing crack does not affect the linear behavior of stress-displacement curves. On the other hand, the onset of coalescence crack between two initial flaws directly leads to the reduction of stress. This reasonable feature verifies the proposed double phase-field model as well. The comparisons between numerical and experimental results for the wing crack and the coalescence crack initial stress

Bi-axial compression of a plate with two inclined flaws

In most previous studies devoted to rock-like materials, only uniaxial compression was generally considered. In practice, these materials are subjected to multi-axial stresses. In this study, the cracking process is investigated in a plate with two initial flaws subjected to bi-axial compression, as shown in Figure 2. Again based on the laboratory tests reported in [START_REF] Bobet | Fracture coalescence in rock materials: Experimental observations and numerical predictions[END_REF][START_REF] Bobet | Numerical modeling of fracture coalescence in a model rock material[END_REF], the configuration "45-a-a" is here considered as an example.

In Figure 7 The stress-displacement curve for the simulation of bi-axial compression is given in Figure 8.An obvious difference between the numerical and experimental results if the friction angle φ = 30 • is used as in the uniaxial compression tests. In this way, the friction angle is assumed to be increased according to the lateral stress. The simulations by using φ = 32 • and 34 • are given. The quantitative stress result is improved while the crack model is not changed (compare with figure 7).

Discussion

The fundamental idea of this double phase-field model is proposed in order to identify the mixed crack problem from the in-situ experiment in the context of radioactive waste disposal. Since the complex damage zone due to the excavation plays an important role on the later thermo-hydromechanical coupling behavior, the concept of our double phase-field model is proposed to distinguish the tensile and compressive-shear damage in [START_REF] Yu | Numerical study of thermo-hydro-mechanical responses of in situ heating test with phase-field model[END_REF].

Recently, Fei and Choo [START_REF] Fei | A phase-field model of frictional shear fracture in geologic materials[END_REF] did a excellent work by using a similar concept of double phasefield model. Following their old work [START_REF] Fei | A phase-field method for modeling cracks with frictional contact[END_REF], they identify the crack condition as: open, stick and slip according to the direction of crack interface. The tensile and compressive-shear crack are properly distinguished with the help of the sign of normal stress and slip criterion of crack. Comparing with the work of Fei and Choo, this work focus on the rock-like material, so that the Mohr-Coulomb criterion is modified to identify the driving force of compressive-shear crack. In this way, the further assumptions of friction crack based on Mohr-Coulomb criterion can be considered in the future study of mixed crack. On the other hand, the tensile and compressive-shear crack can exist at the same time according to the degradation of the stiffness matrix with the help of spectral decomposition. Therefore, both of two double phase-field models work well for the identification of mixed cracks by using several damage variables, but under different physical assumptions, especially for the friction crack.

Conclusion

In this paper, a new phase-field method has been developed for modelling complex cracking modes in rock-like brittle materials subjected to compression-dominating loads. Two interacting crack fields have been introduced for capturing the tensile and shear crack respectively. The emphasis was put on the description of shear crack growth and propagation under compressive stresses through a specific Mohr-Coulomb type frictional sliding criterion.

The combination of two crack fields allows modelling mixed cracking patterns. The proposed method has been implemented in a computer code with the finite element method. It has been applied to describe complex cracking processes in brittle materials under both uniaxial and bi-axial compressive stresses. It was found that the new double phase-field method was able to capture well the main cracking features observed in laboratory tests. Furthermore, the loading stress state corresponding to the initial wing cracks and coalescence cracks is well reproduced. Compared with most previous studies, the new phase-field method is able to well describe the effect of confining stress on the growth of shear cracks commonly observed in rock-like materials. The proposed method can be easily extended to materials with an initial anisotropy and to include plastic deformation.

  s) also indicate the current damage state of material, more precisely d α (x) = 0 being the sound state while d α (x) = 1 the fully damaged one as shown in Figure 1.

Figure 1 :

 1 Figure 1: (a) A sharp crack Γ inside solid body Ω; (b) auxiliary damage variable in sharp crack topology along the 1D line A-A'; (c) and (d) tensile and shear crack phase-field in regularized crack topology along the 1D line A-A'; (e) and (f) regularized representation of tensile and shear cracks with equivalent crack surfaces Γ(d t ) and Γ(d s )

  g s c is calibrated according to the calculated toughness value of tensile crack: g s c =5g t c .

Figure 2 :

 2 Figure 2: Geometry of plate, configuration of initial flaws and boundary conditions for uniaxial (P c = 0) and bi-axial compression tests (P c = 7.5M P a)

  are illustrated in the same Figure. It is seen that from loading step (a) to (b), two external wing cracks and two internal wing cracks are generated from the ends of the initial flaws and they propagate towards the axial compression direction. These wing cracks are all induced by the tensile damage. At the loading step (c) and (d), the wing cracks continue propagating in the axial direction while one shear damage band is generated from the two internal ends of the initial flaws. With the increasing of compression up to step (e) and (f), two external shear cracks are created from the two external ends of the initial flaws. And the internal shear damage band induces a sharp mixed crack which is located at the same position. All these features obtained from the numerical results match very well the experimental evidence.

Figure 3 :

 3 Figure 3: Distribution of tensile and shear damage for the configuration 45-0-2a at three subsequent values of axial displacement (a) u = 0 mm; (b) u = 0.620 mm; (c) u = 0.633 mm; (d) u = 0.635 mm; (e) u = 0.30 mm; (f) u = 0.666 mm

Figure 4 :

 4 Figure 4: Distribution of tensile and shear damage variable for the case of 45-a-2a at three subsequent values of axial displacement (a) u = 0 mm; (b) u = 0.352 mm; (c) u = 0.574 mm; (d) u = 0.620 mm; (e) u = 0.623 mm; (f) u = 0.653 mm

Figure 5 :

 5 Figure 5: Distributions of tensile and shear damage for the configuration of 45-2a-2a at three subsequent steps of axial displacement (a) u = 0 mm; (b) u = 0.370 mm; (c) u = 0.561 mm; (d) u = 0.679 mm; (e) u = 0.685 mm; (f) u = 0.705 mm

Figure 6 :

 6 Figure 6: Left: Curves of stress-displacement for numerical test of '45-0-2a', '45-a-2a' and '45-2a-2a'; Right:

  , the distributions of tensile and shear damage are presented for six different values of axial displacement in a bi-axial compression test with a lateral stress of 7.5 MPa. It is clear that contrary to the uniaxial compression tests considered above, due to the confining effect of lateral stress, the tensile damage is significantly attenuated. The shear cracking is the principal mechanism of damage. The shear crack constitutes the bridge linking two initial flaws. Two second shear cracks are also observed around the external ends of the initial flaws. The cracking patterns provided by the proposed phase-field method agree well with the experimental observations. It seems that the new model is able to describe the effect of confining stress on the kinetics of shear cracks commonly observed in rock-like materials.

Figure 7 :

 7 Figure 7: Distribution of tensile and shear damage for the configuration 45-a-a at three values of prescribed axial displacement at (a) u = 0 mm; (b) u = 0.793 mm; (c) u = 0.801 mm; (d) u = 0.803 mm; (e) u = 0.809 mm; (f) u = 0.850 mm, in a bi-axial compression with a lateral stress of P c =7.5 MPa

Figure 8 :

 8 Figure 8: The curves of stress-displacement simulated by using different friction angle φ = 30 • , 32 • and 34 • comparing with the experimental coalescence crack initial stress, in a bi-axial compression with a lateral stress of P c =7.5 MPa

Table 1 :

 1 Parameter values used in simulations of two inclined flaws test

	Elastic parameters:
	E (GPa)	5.96
	ν (-)	0.24
	Crack fields parameters:

Acknowledgement:

This work has been partially supported by the French National Agency for radioactive waste management (ANDRA).


Keywords: Cell wall, Protein digestion, Polysaccharides, Mechanism, Structure

Pea cell walls have been shown to encapsulate nutrients inside cells, thereby limiting their hydrolysis by digestive enzymes. However, it is unknown how the cell wall performs this barrier function. In particular, this could be due to the presence of specific polysaccharides or, most probably, the organisation of the components within the cell wall. This study aimed to investigate how cell walls prevented protein hydrolysis. To address this objective, isolated cells were obtained using different treatments thought to affect cell walls differently (incubations in water, salt or EDTA solutions) and digested in vitro using a three-phase (oral, gastric and intestinal) model. Purified polysaccharides (cellulose, rhamnogalacturonan I and xyloglucan) and solutions obtained from the incubation of pea fibres and flour were also used in our digestion experiments. We found that protein bioaccessibility (here defined as the amount of protein released from the pea cells) and hydrolysis was lower after the gastric phase for the isolated cells prepared with the salt solution compared to the other treatments. Regardless of the treatment, between 47% and 93% of proteins were released from the cells (bioaccessible) and hydrolysed, respectively. Therefore, after prolonged incubation time proteases seem to be able to penetrate the cells during digestion. In terms of purified polysaccharides, rhamnogalacturonan I had the greatest effect on protein hydrolysis. Incubation solutions made from the pea fibres reduced proteolysis to a greater extent than the pea flour. The present study showed that pea cell walls delay protein digestion mainly via its structural organisation within the cell wall, with purified polysaccharides having a more limited effect.

Pea cell wall polysaccharides and their structural integrity influence protein bioaccessibility and hydrolysis

Maxence Noel a , Frédérique Mayeur-Nickel a , Sandra Wiart-Letort a , Myriam M.-L. Grundy a a PEGASE, INRAE, Institut Agro, 35590, Saint Gilles, France myriam.grundy@inrae.fr, Corresponding author: INRAE, UMR PEGASE, 35590 Saint Gilles, France, Tel: (0033) 2 23 48 50 56; Email: myriam.grundy@inrae.fr (M. Grundy)

Introduction

Pulses such as pea are a potential source of nutrients, such as protein, for both humans and animals. However, previous research has shown that pea cell walls limited the digestibility of nutrients, particularly protein and starch (Grundy, Labarre, Mayeur-Nickel, van Milgen, & Renaudeau, 2023;[START_REF] Junejo | Pea cell wall integrity controls the starch and protein digestion properties in the INFOGEST in vitro simulation[END_REF]. Indeed, the cell wall can influence macronutrient hydrolysis and absorption via a range of mechanisms including encapsulation (Grundy, Edwards, et al., 2016). However, it is still unknown how this effect is achieved and which component(s) of the cell wall is (are) responsible. One hypothesis is that the presence of specific polysaccharides within the cell wall and the cross-links between them and other molecules (proteins and phenolic compounds), and/or their organisation within the cell wall which controls porosity. In order to predict the impact that pulses such as peas can have on health, including digestive functions, it is essential to understand how cell walls modulate nutrient bioaccessibility (i.e., release of a nutrient from the food matrix) and hydrolysis.

Isolated pea cells are an interesting model that can be used to understand the effect of cell walls on nutrient digestion, which is why their study has recently become popular (Bhattarai, Dhital, Mense, Gidley, & Shi, 2018;Edwards et al., 2020;Grundy, Carriere, et al., 2016;[START_REF] Junejo | Pea cell wall integrity controls the starch and protein digestion properties in the INFOGEST in vitro simulation[END_REF][START_REF] Liu | The mechanism of extraction of peanut protein and oil bodies by enzymatic hydrolysis of the cell wall[END_REF][START_REF] Pallares Pallares | Process-induced cell wall permeability modulates the in vitro starch digestion kinetics of common bean cotyledon cells[END_REF][START_REF] Zahir | Effect of soybean processing on cell wall porosity and protein digestibility[END_REF]. Different methods have been used to obtain these isolated cells and even though it is assumed that the cell wall is intact, the impact of the treatment on the polysaccharides within the cell wall remain unclear. Chelating agents, such as cyclohexanediamine tetraacetic acid (CDTA) and ethylenediaminetetraacetic acid (EDTA), can be used to remove the calcium crosslinking the pectic polysaccharides of adjacent cells (located in the middle lamella) together and make cell separation possible, particularly for hard plant tissue [START_REF] Jarvis | The proportion of calcium-bound pectin in plant cell walls[END_REF]. For legumes, cell separation can be achieved with thermal treatment, but this can also result in cell wall swelling (Edwards et al., 2020;[START_REF] Waldron | Plant cell walls and food quality[END_REF]. Finally, soaking in salt solution, often of bicarbonate, combined with thermal treatment, can be used to accelerate the softening of the plant tissue and thereby facilitate cell separation [START_REF] Schoeninger | Processing parameter optimization for obtaining dry beans with reduced cooking time[END_REF].

Cell walls of pea has been previously studied in term of composition, physical properties and functionalities (Dopico, Labrador, & Nicolas, 1986;Grundy et al., 2023;Guevara-Zambrano et al., 2023;[START_REF] Jebalia | Cell walls of extruded pea snacks: Morphological and mechanical characterisation and finite element modelling[END_REF][START_REF] Junejo | Pea cell wall integrity controls the starch and protein digestion properties in the INFOGEST in vitro simulation[END_REF][START_REF] Noguchi | Determination of chemical structure of pea pectin by using pectinolytic enzymes[END_REF][START_REF] Talbott | Molecular size and separability features of pea cell wall polysaccharides: implications for models of primary wall structure[END_REF]. Depending on the variety, pea dietary fibre content ranges from 13 to 31%, with most of it being classified as insoluble fibre (Brummer, Kaviani, & Tosh, 2015;Grundy et al., 2023;[START_REF] Kan | Comparative study on the chemical composition, anthocyanins, tocopherols and carotenoids of selected legumes[END_REF]. Pea cell walls are made predominantly of pectin (mainly rhamnogalacturonan I) and hemicellulose (mainly xyloglucan) (Brillouet & Carré, 1983;Dopico et al., 1986;[START_REF] Noguchi | Determination of chemical structure of pea pectin by using pectinolytic enzymes[END_REF].

One of the difficulties faced when studying cell walls is to obtain precise, reliable values on their composition and structure. Given that cell walls have to be disrupted to extract the polysaccharides they contain, information about their organisation, spatial structure, and certain cross-links is lost. Other, indirect means ought to be employed to compensate for this, for instance particle size measurements and microscopy observations can provide information about the structure and integrity of the cell wall [START_REF] Holland | Plant cell walls: Impact on nutrient bioaccessibility and digestibility[END_REF].

However, these approaches still do not permit to identify the cross-links between the polysaccharides within the cell wall. Using both intact cell walls, individual (extracted) and solubilised polysaccharides can overcome this missing information. Therefore, this study aimed to shed light on the cell wall composition and structure responsible for the observed decrease in protein digestibility. To do so, cells from pea seeds were isolated using different treatments. In parallel, individual polysaccharides and solutions obtained following incubation of pea flour or fibres were added to pea protein to investigate their potential interaction and impact on proteolysis.

Materials and Methods

Materials

Pea (Pisum sativum L.) seeds, protein, and fibres were provided by Roquette (Vic-sur-Aisne, France). As indicated by the supplier, pea proteins were isolated through wet extraction of yellow peas, separated by flocculation and spray-dried. Pea fibres were obtained from the wet extraction of yellow peas, and were composed mainly of internal cell walls from the cotyledons, (the hull having been removed) and residual starch (see Table 1 for composition).

Those fibres were separated by centrifugation and decantation, and then dried. Pea cells were isolated using four different methods based on previous work (Bhattarai, Dhital, Wu, Chen, & Gidley, 2017;Grundy, Wilde, Butterworth, Gray, & Ellis, 2015;Guevara-Zambrano et al., 2023;[START_REF] Junejo | Pea cell wall integrity controls the starch and protein digestion properties in the INFOGEST in vitro simulation[END_REF]: i) incubation of pea seeds in water for 16 h, cooking at 60°C for 1 h and separation with a high shear mixer (IKA® T25 ultra-turrax, treatment 1, Mixer water 60°C), ii) incubation in water for 16 h at room temperature, cooking at 60°C for 1 h and separation gently crushing the peas with a mortar and pestle (treatment 2, Water 60°C), iii) incubation in NaHCO3 (1.5%, w/v) and Na2CO3 (0.5%, w/v) for 16 h at room temperature (pH 9.5), cooking at 60°C for 1 h and separation by gently crushing the peas with a mortar and pestle (treatment 3, Salt 60°C), and iv) incubation in EDTA (50 mM) for 16 h at room temperature (pH 7) without a cooking step, and separation gently crushing the peas with a mortar and pestle (treatment 4, EDTA). Each cell preparation was then loaded on a stack of sieves of 1 000, 150 and 71 µm aperture as well as a sieve base to collect the liquid. After washing the preparation with water to separate the cells from other particles and cell fragments, the cells were collected on the 71 µm sieve. For treatment 4, a dialysis step using a dialysing membrane (Float-ALyzer G2 10 mL, 3.5-5 kDa, Merck) was added to remove the EDTA as the latter can affect enzymes activity. Pea flour was obtained by grinding pea seeds using a hammer mill (40 HP, capacity 3T/H; Stolz, Pontivy, France) fitted with a 2.5 mm grid.

Pea flour and fibres were incubated in deionised water at 4°C on a rotator for periods of 6, 24 and 72 h to obtain solubilised pea fibres. The end point of 72 h was chosen to represent the maximum time that the pea material (flour or extracted fibres) could reside in the gastrointestinal tract before being fermented by colonic microorganisms. This complemented the experiments performed on pure polysaccharides, cellulose, rhamnogalacturonan and xyloglucan, extracted from different plant sources (see below for references).

Pepsin from porcine gastric mucosa (#P6887, 2145 U/mg of solid), bovine bile extract (#B3883), pancreatin from porcine pancreas (#P7545, 3.4 U/mg of solid based on trypsin activity) and cellulose (#310697) were purchased from Merck (Saint Quentin Fallavier, France). Rhamnogalacturonan I (from soybean, #P-RHAGN) and xyloglucan (from tamarind, #P-XYGLN) were purchased from Megazyme (Libios, Vindry-sur-Turdiner, France). All other chemicals, solvents and reagents were purchased from Merck (Saint Quentin Fallavier, France) or Thermo Fisher Scientific (Illkirch-Graffenstaden, France).

Physical and chemical characterisation of materials

Each pea cell preparation and the polysaccharides (cellulose, rhamnogalacturonan and xyloglucan) were analysed for moisture (oven-dried at 102°C), protein (Dumas method with a nitrogen conversion factor of 5.4, LECO, FP828 Carbon/Nitrogen/Protein Determinator, Villepinte, France, AOAC 968.06), and total dietary fibre (method AOAC 991.43) content.

The pea flour and fibres were analysed as described in Grundy et al. (2023) and presented in Table 1. The average particle size of the cells obtained via the different treatments was measured using Malvern laser diffraction particle sizer 3000 equipped with a dispersant unit (Hydro 2000G) filled with water (Malvern Instruments Ltd.). Triplicate measurements were made for each analysis.

Table 1 Composition of the pea flour and fibres (% on a dry weight basis).

In vitro digestion

The cells isolated with treatments 2, 3 and 4 were digested using the three phase Infogest standardised protocol (Brodkorb et al., 2019;Grundy, Tang, van Milgen, & Renaudeau, 2022).

The cells isolated with treatment 1 were not digested in vitro as their cell walls were ruptured (protein 100% bioaccessible and highly hydrolysable).

Pea protein was mixed with the different polysaccharides (either purified or solubilised following incubation) to identify their effect on protein hydrolysis. the extent of protein hydrolysis was assessed after the intestinal phase.

The different materials (cells, pea protein plus polysaccharides, and pea protein plus incubation solutions) were added to simulated salivary fluid so the quantity of protein and dietary fibre (or purified polysaccharides) was equal to 50 mg and 100 mg (dry weight basis), respectively. The composition of the simulated fluids used for each digestion phase can be found elsewhere (Brodkorb et al., 2019). After about 6 h of incubation (2 min of oral phase pH 7, 2 h of gastric phase at pH 3, and 4 h of intestinal phase at pH 7), the enzymatic reactions were stopped by increasing the pH to 9 with 5 N NaOH at the end of the intestinal phase.

Each digestion was performed in triplicate, without (blank sample, to determine protein bioaccessibility) and with (digested sample, determination of protein hydrolysis) enzymes.

Protein and proteolysis products assessment

For the isolated cells, the extent of proteolysis was estimated from the crude protein measurement of the recovered in vitro digesta by the Dumas method (see section 2.2). Briefly, following centrifugation at 2,500 g at 4 °C for 10 min, the pellet containing the undigested material was washed, filtered using a cell-strainer (Falcon®, 40 µm aperture) and dried at 80°C overnight as previously described (Grundy et al., 2022).

For the other samples (digestions with the purified polysaccharides and the incubation solutions), the OPA method was used to quantify the amino groups. Thus, after centrifugation, the supernatant was collected, proteins precipitated in methanol and the small peptides hydrolysed with sulfuric acid as described in [START_REF] Sousa | In vitro digestibility of dietary proteins and in vitro DIAAS analytical workflow based on the INFOGEST static protocol and its validation with in vivo data[END_REF].

Microstructural characterisation

Samples were mounted on microscopy slides, and then visualised with an Apotome microscope and Zen software (Apotome, Zeiss, France). Some of the cells were also dyed with calcofluor white and fast green FCF to identify cell wall and protein components, respectively. Images were captured using 10x and 20x objective lenses.

Statistical analysis

The data were analysed using R studio version 4.1.2. For all tests, the significance level was set at P < 0.05 (2 tailed) and all data were expressed as means of triplicates. The differences in protein bioaccessibility and digestibility between the pea materials were assessed by one-way analysis of variance (ANOVA) followed by Tukey's post-hoc test.

Results and discussion

This study was conducted to better understand how pea cell walls influenced protein digestion, either via its structural organisation (physical barrier) and/or via the properties of its constitutive polysaccharides (e.g., "binding" with proteases, or increase in the viscosity of the gastrointestinal content). To achieve this, we employed an approach similar to our previous work [START_REF] Grundy | The impact of oat structure and β-glucan on in vitro lipid digestion[END_REF] where dietary fibres with different levels of complexity were used:

i) isolated cells obtained with different treatments (it is assumed that the cell walls are intact but each treatment is likely to have a different impact of the cell wall integrity and structure),

ii) pea polysaccharides solubilised following the incubation of either flour or fibres extracted from pea, and iii) purified polysaccharides chemically similar to the ones found in pea cell walls but extracted from other plant sources (commercially available). A set of analyses were then carried out to characterise the pea cells (nutritional composition, particle size, and microscopy) before and after simulated digestion, and the extent of proteolysis estimated from the material recovered after digestion (either the cells or products of proteolysis).

Characterisation of the studied materials

Differences in the content of protein and dietary fibre were observed between the four cells preparations (Table 2). Some protein was lost during the isolation of the cells, notably for EDTA (14.2%) compared to Salt 60°C (18.6%, P < 0.05). However, the dietary fibre content was higher for the EDTA compared to the other treatments (TDF value of 56% for EDTA, between 6.1 and 23.7% for the other cell samples, P < 0.05), in particular the insoluble fraction seemed to have been more impacted by the preparation. It is likely that this fraction corresponded to resistant starches that were not gelatinised due to the absence of cooking. Even though a thermal treatment is included in the TDF analysis (100°c for 4 h in presence of pancreatic α-amylase and amyloglucosidase), it is possible that either the enzymes could not penetrate the cell wall or the starch was not fully gelatinised and/or hydrolysed (Dhital, Bhattarai, Gorham, & Gidley, 2016). Therefore, the analysis did not discriminate between cell wall components and resistant starches.

As expected, the incubations in water did not sufficiently weaken the middle lamella, which resulted in cell wall rupture (Fig. 1.1 and 1.2), losses in dietary fibre (6.08 and 13.83% for Mixer water 60°C and Water 60°C, respectively) and protein (8.6% compared to 18.6 and 14.6% for Salt 60°C and EDTA, respectively). Therefore, differences in the integrity of the cells were visible in the microscopy images (Figure 1). The cells isolated using water and the mortar appeared less damaged than those obtained after mixing. However, a limited number of cells seemed to have intact cell walls with starch (gelatinised) "holding" the cell shape together, some cells also appeared emptied of their content (Fig. 1.2A and 1.2B). Because of the absence of heat treatment, the starch for the EDTA cells were not gelatinised (starch granules appear intact). The cell wall of the cells isolated with Salt 60°C and EDTA seemed to be undamaged as demonstrated by the blue "layer" (cell wall stained with calcofluor White) delimiting the cells (Fig. 1.3a, 3b, 4a and 4b). In addition, for those treatments the protein was only observed within the cells (stained in green) unlike the treatments with water where some of the nutrients (starch and proteins) were "solubilised" with empty cells and cell wall fragments present in the preparation (Fig. 1.1 and 1.2). These observations were confirmed with the particle size distributions (Fig. 2) that showed a mean cell size of 135 µm for the Salt 60°C and EDTA treatments. This value is slightly smaller than those reported in other studies, that is around 150 µm (Edwards et al., 2020;[START_REF] Junejo | Pea cell wall integrity controls the starch and protein digestion properties in the INFOGEST in vitro simulation[END_REF], which could be due to differences in the variety of pea used. A second peak at around 20 µm can be identified, especially for the two water treatments, which is likely to correspond to starch granules, protein aggregates, and other cell fragments solubilised during the analysis in the dispersant unit (as distinguishable in Fig1.1B and 1.2B).

Fig. 2.

Particle size distribution of the isolated cells obtained with the four treatments.

Impact of cell isolation treatment on protein digestion

Figure 3 reveals that more than half of the protein remained inside the cells during the incubation (on average ~52% of protein remained encapsulated), although the amount for the cells obtained with treatment 2 (water 60°C) was actually lower given that some protein was lost during the preparation process (only 8.6 % of protein remained in the sample, see Table 2). This is visible in Fig. 4.1 where the cells appear more damaged for the water treatment than for Salt 60°C and EDTA. The protein hydrolysis after the gastric phase was lower (P < 0.05)

for the cells prepared with the treatment 3 (Salt 60°C) compared to the other treatments.

However, the digestions at the end of the intestinal phase were similar for the three types of cells (~93%). The extent of protein hydrolysis in pea cells reported in the literature is lower than that of the present work (Bhattarai et al., 2017;[START_REF] Junejo | Pea cell wall integrity controls the starch and protein digestion properties in the INFOGEST in vitro simulation[END_REF]. This could be due to the incubation time that was longer (6 h instead of 2 or 4 h in other studies) and to the variety and growing conditions of the peas that may also have been different. Indeed, it is known that the variety and the conditions in which a plant grows can influence its composition, including dietary fibre, and thereby the nutritional value and digestibility of the food products derived from it (Givens, Davies, & Laverick, 2004;Grundy et al., 2020;Hall, Hillen, & Robinson, 2017). Particle size analysis of the recovered cells at different time points of digestion shows an overall increase in the mean volume of particles with a 135 µm diameter (main peak) demonstrating the disappearance of most of the solubilised material (loss of the smaller, soluble and hydrolysable particles such as protein aggregates) (Fig. 5). This is particularly true for the salt treatment, whereas little changes in cell size were observed for the EDTA treatment. Even though the size of the cells was unchanged, their content was hydrolysed as illustrated in Fig 4

( 3A and 3B) where "ghost" cells can be seen. EDTA, which hindered their hydrolysis. More importantly, as mentioned above, having not been hydrothermally treated, the starch was not gelatinised, the compact/semi-crystalline structure thus preventing amylase access to the glycosidic bonds (Fig. 4.2B, 2C and 3C).

Regardless of the treatment, a significant amount of starch remained inside the cells after the gastric phase. The protein however was highly hydrolysed indicating that proteases had access to their substrate. Therefore, while the extra cellular proteins present in the simulated digestive fluids were readily "available" for protease hydrolysis, it is also likely that the enzymes penetrated the cell wall (and hydrolyse the proteins inside the cells) as shown with the diffusion of Fluorescein isothiocyanate (FITC)-labelled dextrans of 10 and 40 kDa in some of the cells (Figure 2S). This size range includes the molecular weight of pepsin (35 kDa) and trypsin (24 kDa). In another study, it was reported that amylase (50 kDa) was unable to penetrate the cell wall of pea (Bhattarai et al., 2017). Here, some amylase must have entered the cell given the degradation of starch in the Salt 60°C samples. The longer incubation time used in this study for the intestinal phase may have resulted in cell wall swelling and thereby facilitating the diffusion of larger molecules through the cell wall. Another recent study also found that chickpea proteins were hydrolysed despite being encapsulated into intact cells at the time of ingestion [START_REF] Perez-Moral | In vitro protein bioaccessibility and human serum amino acid responses to white bread enriched with intact plant cells[END_REF]. Given that pectin is the main polysaccharide thought to influence cell wall porosity, it is not surprising that the EDTA treatment increased the ability of the enzymes to diffuse inside the cell, and consequently protein hydrolysis [START_REF] Mccann | Plant cell wall biology: polysaccharides in architectural and developmental contexts[END_REF]. The preparation of isolated cells would have therefore disrupted the cell wall integrity and, throughout digestion, solubilised their components (pectin and hemicellulose while the cellulosic backbone remained intact).

Effect of pea cell wall and its components on protein digestion

In order to investigate the role of isolated cell wall components on protein hydrolysis, polysaccharides of different chemical natures as well as solubilised pea dietary fibres were utilised in the digestibility experiments. Therefore, the extent of protein hydrolysis was monitored in the presence of three main classes of polysaccharides making up the cell wall composition (i.e., cellulose, hemicellulose and pectin) and dietary fibres obtained after the incubation of pea fibres and flour. These experiments demonstrated that rhamnogalacturonan I (pectin, RhG) and xyloglucan (hemicellulose, XG) reduced pea protein hydrolysis significantly compared to the pea cell walls (Pea). The intestinal digesta containing RhG and XG looked like a "gel" that was generated from these viscous polysaccharides (Fig. 3S image B).

Polysaccharides extracted from pea are not commercially available. Therefore, even though similar in their structure it is likely that the polysaccharides used in this study, especially RhG and XG, had different physico-chemical properties to the pea pectin and hemicellulose solubilised following incubation. Also, pea cell walls may contain polysaccharides of different sizes and chemical structure (branching and sugar composition) (Brillouet & Carré, 1983;[START_REF] Talbott | Molecular size and separability features of pea cell wall polysaccharides: implications for models of primary wall structure[END_REF]. For this reason, pea fibres (consisting of cell wall fragments) and flour were left to incubate in order to solubilise the pea polysaccharides. This enabled us to study the effect of these polysaccharides "outside" the cell wall on protein hydrolysis. Thus, solutions obtained from the incubation of fibres also reduced protein hydrolysis, but to a lesser extent than RhG and XG. Incubation of flour had even less of an impact on protein digestion.

Incubation time did not modulate any of these effects therefore the "active" compounds (probably polysaccharides but maybe also phenolic compounds) must have been released at the early stage of the incubation. Rupture of pea cells appeared to be necessary to release polysaccharides and other cell wall components (e.g., phenolic compounds) that could interact with either the proteases or the proteins to influence protein hydrolysis.

The functionality of extracted, purified compounds is often different to when the same compounds are still present within the food matrix. Indeed, the extraction process is likely to change their chemical composition, molecular size, their (micro)environment and as a consequence, their biological activity. This is well known in the case of β-glucan, which can lose its properties once extracted [START_REF] Wang | Oat beta-glucan: physico-chemical characteristics in relation to its blood-glucose and cholesterol-lowering properties[END_REF]. In this study, we also revealed that cell wall polysaccharides did not have the same impact on protein hydrolysis as part of the cell wall, solubilised in the aqueous environment or extracted. Therefore, the hindering of protein hydrolysis by encapsulation mechanism is time dependent, with the cell wall being degraded as the food travels through the gastrointestinal tract. In particular, swelling and solubilisation of cell wall components (pectin) can increase porosity, thereby allowing the penetration of proteases. Longer incubation time was selected in this work based on in vivo data (Grundy et al., 2023). The Infogest protocol has been developed primarily for dairy sources of proteins, however plant-based foods can be difficult to digest due to the presence of tissue structure containing intact cells, the presence of antinutritional factors and the structure of the proteins (Grundy et al., 2022;[START_REF] Sousa | Protein digestion of different protein sources using the INFOGEST static digestion model[END_REF].

Conclusion

The current study demonstrated that different treatments applied to peas to obtain isolated cells has consequences on the cell wall structure and porosity, which in turn affects the kinetics of protein hydrolysis. Therefore, the way in which a plant tissue is transformed to generate food products is crucial for its behaviour in the gastrointestinal tract and its subsequent health outcome. Specifically, modulating the porosity of cell walls can be a strategy to improve the postprandial absorption of certain nutrients. However, as far as pulses are concerned, it may be of interest to reduce encapsulation of proteins and improve their digestibility. We demonstrated that even though the cell walls appear to maintain their apparent structure, most of the cell contents were hydrolysed after intestinal digestion. Questions still remain regarding the fate of those cell walls in the colonic compartment and the subsequent impact on the composition and functionality of the gut microbiota.

Fig. 1 .

 1 Fig. 1. Microscopy images of the cells at baseline for different preparations: treatment 1 (Mixer

Fig. 3 .

 3 Fig. 3. Amount of protein released from the cells (bioaccessible) and hydrolysed following

Fig. 4 .

 4 Fig. 4. Microscopy images of pea cells recovered following in vitro digestion for the blank (1,

Fig. 5 .

 5 Fig. 5. Particle size distribution of the collected samples after gastric and intestinal incubation

Fig. 6 .

 6 Fig. 6. Protein digestibility expressed as glutamic acid equivalents (in mM) corresponding to

Table 2

 2 Protein and dietary fibre composition of the cells obtained from the four treatments

		Mixer water 60°C Water 60°C	Salt 60°C	EDTA
	DM (%)	24.75	21.07	26.05	24.03
	Protein (%) a	5.09	8.57	18.55	14.20
	TDF (%) a	6.73	13.83	23.66	55.84
	IDF (%)	5.43	13.22	21.97	52.58
	SDF( %)	1.30	1.24	3.39	4.89

a Data expressed on a dry weight basis.

Acknowledgements

This work was funded by INRAE.