Aurelien Doitrand
email: aurelien.doitrand@insa-lyon.fr

Dominique Leguillon

Rafael Estevez

Experimental determination of generalized stress intensity factors from full-field measurements

Keywords: Generalized stress intensity factors, Digital Image Correlation, Square hole, Mode mixity

Unidirectional tensile tests of PMMA square hole specimens result in mixed mode crack initiation at the square hole corners, highlighted by inclined cracks with respect to the V-notch bisector. Generalized stress intensity factors (GSIF) of the opening and shear modes at a V-notch are derived experimentally in a direct manner from digital image correlation (DIC) displacement and strain fields by means of a path independent integral. With increasing square hole side, the opening mode GSIF increases whereas the shear mode one decreases so that the mode mixity globally decreases, together with the crack deflection with respect to the V-notch bisector. GSIFs predictions obtained by means of finite element calculations are in good agreement with those determined experimentally.

u I : opening mode displacement vector in the V-notch axis system u I r , u I ϕ : opening mode displacement components in the V-notch axis system u II : shear mode displacement vector in the V-notch axis system u II r , u II ϕ : shear mode displacement components in the V-notch axis system u J t , u J n : cartesian displacement components in the V-notch axis system

Introduction

The displacement and stress fields in the vicinity of a singular point can be described by a characteristic exponent λ and a coefficient representing the intensity of the field: the generalized stress intensity factor (GSIF). The GSIF of a singularity is a relevant parameter for failure prediction. In the well-known particular case of a crack in isotropic homogeneous media, λ=0.5 and the GSIF corresponds to the usual stress intensity factor (SIF). The SIF can be computed using J-integral [START_REF] Rice | A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks[END_REF]. However, this method is specifically dedicated to cracks and cannot be extended to compute the GSIF of other singularities such as, for instance, a corner, a V-notch in a homogeneous material. Therefore, alternative approaches have to be employed to compute the GSIF. The numerical computation of 2D GSIF can be performed using a least-squares fitting procedure of Finite Element (FE) nodal displacements to the asymptotic fields in the vicinity of the singular point [START_REF] Barroso | A least squares procedure for the evaluation of multiple generalized stress intensity factors at 2D multimaterial corners by BEM[END_REF][START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF][START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF]. Lazzarin et al. [START_REF] Lazzarin | Rapid calculations of notch stress intensity factors based on averaged strain energy density from coarse meshes: Theoretical bases and applications[END_REF] extracted the 2D opening and shear mode GSIFs of a V-notch based on the strain energy density computed considering the leading order terms of the William's solution [START_REF] Williams | Stress Singularities Resulting from Various Boundary Conditions in Angular Corners of Plates in Extension[END_REF]. They highlighted the possibility of using relatively coarse meshes since this approach involves the mean value of the strain energy density that can be obtained from the nodal displacements. The most accurate method for GSIF computation is based on a path independent integral [START_REF] Labossiere | Calculation of stress intensities at sharp notches in anisotropic media[END_REF][START_REF] Labossiere | Stress intensities at interface corners in anisotropic bimaterials[END_REF][START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF]. Although it can be employed to compute the SIF in the case of a crack (even in anisotropic media), this method differs significantly from J-integral which is strictly dedicated to SIF calculation for a crack. The GSIF use is particularly employed in Finite

Fracture Mechanics (FFM) framework for crack initiation prediction in the vicinity of a singularity or a stress concentrator based on the matched asymptotic (MA) approach of the coupled criterion (CC) [START_REF] Doitrand | Finite element procedure implementation of the coupled criterion: matched asymptotic and full finite element approaches[END_REF][START_REF] Leguillon | Strength or toughness ? A criterion for crack onset at a notch[END_REF][START_REF] Leguillon | Size effects for crack initiation at blunt notches or cavities[END_REF][START_REF] Weissgraeber | A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress-raisers[END_REF]. This method has recently been extended to 3D [START_REF] Doitrand | Computation of generalized stress intensity factors of 3D singularities[END_REF], which was the missing tool for applying the MA approach of the CC in 3D [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF][START_REF] Leguillon | An attempt to extend the 2D coupled criterion for crack nucleation in brittle materials to the 3D case[END_REF].

Although some methods exist to compute from experimental data the stress intensity factors at a crack tip in 2D or in 3D [START_REF] Limodin | Crack closure and stress intensity factor measurements in nodular graphite cast iron using three-dimensional correlation of laboratory X-ray microtomography images[END_REF][START_REF] París | Photoelastic determination of K I and K II : A numerical study on experimental data[END_REF][START_REF] Réthoré | Noise-robust stress intensity factor determination from kinematic field measurements[END_REF][START_REF] Réthoré | Automatic crack tip detection and stress intensity factors estimation of curved cracks from digital images[END_REF][START_REF] Roux | Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches[END_REF][START_REF] Roux | Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2D and 3D cracks[END_REF], few authors estimated GSIFs from experimental data: Based on a least-squares fitting of the asymptotic displacement fields approach (called the overdeterministic method [START_REF] Ayatollahi | Determination of NSIFs and coefficients of higher order terms for sharp notches using finite element method[END_REF]), Torabi et al. recently determined the GSIF at a sharp [START_REF] Torabi | Experimental determination of the notch stress intensity factor for sharp Vnotched specimens by using the digital image correlation method[END_REF] or blunt [START_REF] Torabi | On the use of digital image correlation method for determining the stress field at blunt V-notch neighborhood[END_REF] V-notch. Dunn et al. [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] carried out experiments on V-notch PMMA specimens with various opening angles and depths. They computed numerically 2D V-notch GSIFs for the critical loads measured experimentally and showed that the GSIF of the corner singularity was the relevant parameter to predict crack initiation. The GSIF was computed by least square fitting of the nodal displacements along the V-notch flanks from FE calculations. The critical GSIF was determined to be constant for a given V-notch angle whatever the notch depth, except may be for large V-notch angles. Labossiere and Dunn [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] performed four point bending tests on bimaterial specimens in order to study crack initiation occurring at the bimaterial interface corner under an opening mode. They obtained the 3D interface corner GSIF by means of a least-squares fitting procedure of the FE displacements fields computed on the whole structure in the vicinity of the singular point to the asymptotic displacement fields. The objective of this work is the direct experimental determination of GSIFs from digital image correlation (DIC) displacement and strain fields. Tensile tests on PMMA square hole specimens are presented in Section 3. Section 4 describes the GSIF calculation method that is employed as a postprocessing of either FE calculations or DIC fields. In Section 5, the coupled criterion for crack initiation prediction is recalled. Finally, GSIFs obtained from DIC or FE calculations are compared in Section 6.

Experiments

The material under investigation is a commercial extruded PMMA. Several rectangular specimens were manufactured using laser cutting to obtain square holes which sides are parallel to the specimen sides. The specimens were tested under uniaxial tension on a 20kN Zwick machine. The choice of such a configuration is motivated by the fact that it leads to a mixed mode crack initiation, the mode mixity depending on the square hole size c. It is assumed here that c is small compared to the specimen width w. Experimental tests were performed on a total of 10 configurations with square hole sizes varying from 1.2 to 15.2 mm, with 5.76+0.05 mm thick and 150+0.2 mm long specimens. The other dimensions are reported in Table 1. A speckle pattern, obtained using a white and black paint air spray, was applied on the specimens prior to the tests. The prescribed displacement rate is 0.5 mm/min, which is small enough to ensure an overall quasi-static loading. The test is optically recorded at 1 Hz frequency with 1 s exposure time using a 29M Pixels Allied Vision Prosilica GX camera. Specimens without hole were also tested to measure the material Young's modulus (E =2.7+0.1 GPa) and Poisson's ratio (ν =0.39+0.005) by averaging the longitudinal and transverse strain fields obtained by DIC (Fig. 1).

All DIC computations presented in this work have been performed using UFreckles software [START_REF] Réthoré | UFreckles (Version v 2.0)[END_REF]. They are based on a global DIC method [START_REF] Besnard | Finite-element displacement fields analysis from digital images: application to Portevin-Le Chatelier bands[END_REF][START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF], which requires the definition of a mesh on the reference image. The displacement field minimizing the gray level difference norm between the reference image and the deformed image corrected by this displacement field is computed at the mesh nodes. A regular mesh with 4-node elements is adopted and the mesh size is chosen so that the integration contour radius is larger than five elements.

For each specimen under investigation, cracking occurs at two square hole corners being located at each side of the square. It can be noted that the two cracks on both sides of the square hole initiate quasi-simultaneously. However, it has been observed for some specimen that one crack initiates slightly earlier than the second one which leads to the specimen immediately breaking into two parts. Figs. opening mode loading the crack would initiate along the square hole V-notch bisector, which is not observed. Moreover, the smaller the hole size, the larger the angle between the V-notch bisector and the crack, which reveals that the smaller the hole size, the smaller the amount of opening mode relatively to

(b) (c) (a) (d) (e) (f) (g) (h) (i) (j) 10mm
Loading direction shear mode. Note that, to avoid any ambiguities, mixed mode crack initiation means that the loading in the vicinity of the corner is a mixture of opening and shear. But, the material being homogeneous, the crack selects the direction where it is in only opening mode.

GSIF computation

Singularity exponent, primal and dual fields

First, let us recall the generic form of Williams' expansion in the vicinity of a V-notch [START_REF] Williams | Stress Singularities Resulting from Various Boundary Conditions in Angular Corners of Plates in Extension[END_REF] U

(r, ϕ) = C + kr λ u(ϕ) + ... ( 1
)
where U is the solution to an elastic problem, C holds for the rigid translation, r and ϕ are the polar coordinates emanating from the root of the V-notch. The exponent λ and its associated mode u(ϕ)

depend only on the local geometry and possibly local elastic properties while the GSIF k is a function of the global geometry of the specimen and the intensity of the applied load as well as the way it is applied. In this approach, it is assumed that the local geometry and material properties around the singular point depend only on the polar angle ϕ.

Solutions of the form r λ u(ϕ) are sought under the assumption that there is locally no body forces and that the neighbouring faces undergo vanishing conditions (expressed either in terms of displacement or traction) but disregarding remote boundary conditions. Non homogeneous conditions can be added afterward. the elasticity problem. The analytical and numerical resolution of the eigenvalue problem, of which λ and u are solution, have been widely developed in [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF] and more recently in [START_REF] Doitrand | Finite element procedure implementation of the coupled criterion: matched asymptotic and full finite element approaches[END_REF] dedicated to the implementation of the CC in a commercial code. It is not explained again herein. The primal opening (I) and shear (II) modes displacement and stress fields in polar coordinates in the V-notch frame (Fig. 4) are given by:

u I (r, ϕ) = r λ 1        u I r (ϕ) u I ϕ (ϕ)        and u II (r, ϕ) = r λ 2        u II r (ϕ) u II ϕ (ϕ)        (2)
with (cf., e.g., [START_REF] Seweryn | Elastic stress singularities and corresponding generalized stress intensity factors for angular corners under various boundary conitions[END_REF][START_REF] Vasilopoulos | On the determination of higher order terms of singular elastic stress fields near corners[END_REF])

u I r (ϕ) = [cos((1 + λ 1 )ϕ * ) + λ L +3µ L -λ 1 (λ L +µ L ) (λ L +µ L )(1-λ 1 ) sin(ω(1+λ 1 )/2) sin(ω)(1-λ 1 )/2) cos((1 -λ 1 )ϕ * )]/(2µ L λ 1 γ I-0 ϕϕ ) u I ϕ (ϕ) = [-sin((1 + λ 1 )ϕ * ) -λ L +3µ L +λ 1 (λ L +µ L ) (λ L +µ L )(1-λ 1 ) sin(ω(1+λ 1 )/2) sin(ω)(1-λ 1 )/2) sin((1 -λ 1 )ϕ * )]/(2µ L λ 1 γ I-0 ϕϕ ) u II r (ϕ) = [sin((1 + λ 2 )ϕ * ) + λ L +3µ L -λ 2 (λ L +µ L ) (λ L +µ L )(1+λ 2 ) sin(ω(1+λ 2 )/2) sin(ω)(1-λ 2 )/2) sin((1 -λ 2 )ϕ * )]/(2µ L λ 2 γ II-0 rϕ ) u II ϕ (ϕ) = [cos((1 + λ 2 )ϕ * ) + λ L +3µ L +λ 2 (λ L +µ L ) (λ L +µ L )(1+λ 2 ) sin(ω(1+λ 2 )/2) sin(ω)(1-λ 2 )/2) cos((1 -λ 2 )ϕ * )]/(2µ L λ 2 γ II-0 rϕ ) (3)
with

γ I-0 ϕϕ = (1+λ 1 )sin(ω(1+λ 1 )/2) (1-λ 1 )sin(ω(1-λ 1 )/2) -1 γ II-0 rϕ = 1 -(1-λ 2 )sin(ω(1+λ 2 )/2) (1+λ 2 )sin(ω(1-λ 2 )/2) (4)
where λ L , µ L are the material Lamé coefficients, ϕ * = ϕ -π/2, ω = 3π 2 . The singularity exponents corresponding to primal opening and shear modes, which can be obtained by solving a transcendent equation, are λ 1 =0.545 and λ 2 =0.908. The opening (I) and shear (II) primal stress fields in the vicinity of the V-notch can be expressed in polar coordinates as:

σ I (r, ϕ) =                σ rr σ ϕϕ σ rϕ                = r λ 1 -1                σ I rr (ϕ) σ I ϕϕ (ϕ) σ I rϕ (ϕ)                and σ II (r, ϕ) = r λ 2 -1                σ II rr (ϕ) σ II ϕϕ (ϕ) σ II rϕ (ϕ)                (5)
with

σ I rr (ϕ) = [cos((1 + λ 1 )ϕ * ) + 3-λ 1 1-λ 1 sin(ω(1+λ 1 )/2) sin(ω)(1-λ 1 )/2) cos((1 -λ 1 )ϕ * )]/γ I-0 ϕϕ σ I ϕϕ (ϕ) = [-cos((1 + λ 1 )ϕ * ) + 1+λ 1 1-λ 1 sin(ω(1+λ 1 )/2) sin(ω(1-λ 1 )/2) cos((1 -λ 1 )ϕ * )]/γ I-0 ϕϕ σ I rϕ (ϕ) = [-sin((1 + λ 1 )ϕ * ) + sin(ω(1+λ 1 )/2) sin(ω(1-λ 1 )/2) sin((1 -λ 1 )ϕ * )]/γ I-0 ϕϕ σ II rr (ϕ) = [sin((1 + λ 2 )ϕ * ) + 3-λ 2 1+λ 2 sin(ω(1+λ 2 )/2) sin(ω(1-λ 2 )/2) sin((1 -λ 2 )ϕ * )]/γ II-0 ϕϕ σ II ϕϕ (ϕ) = [-sin((1 + λ 2 )ϕ * ) + sin(ω(1+λ 2 )/2) sin(ω(1-λ 2 )/2) sin((1 -λ 1 )ϕ * )]/γ II-0 ϕϕ σ II rϕ (ϕ) = [cos((1 + λ 2 )ϕ * ) -1-λ 2 1+λ 2 sin(ω(1+λ 2 )/2) sin(ω(1-λ 2 )/2) cos((1 -λ 2 )ϕ * )]/γ II-0 rϕ (6)
The primal displacements and stresses corresponding to opening (J = I) or shear (J = II) modes can finally be expressed using Cartesian coordinates in the V-notch frame (O , t, n) (cf. corresponds to the V-notch corner) as:

       u J t u J n        = P        u J r (ϕ) u J ϕ (ϕ)        (7)
and

    σ J tt σ J tn σ J tn σ J nn     = P     σ J rr σ J rϕ σ J rϕ σ J ϕϕ     P T (8)
where

P =     cos ϕ -sin ϕ sin ϕ cos ϕ     (9)
P is the matrix that allows transforming polar to Cartesian coordinates in the V-notch frame.

Moreover, it is emphasized that the dual displacement and stress fields can be obtained by replacing λ 1 by -λ 1 and λ 2 by -λ 2 in the expressions given for the primal displacement and stress fields.

The contour integral calculation (cf. Section 4.2) requires that the DIC fields together with the primal and dual fields are expressed in the same frame. In practice, we choose to transfer the DIC fields to the V-notch frame. Denoting γ the angle between the DIC frame (O, x, y) and the V-notch frame (O , t, n), we can express the DIC displacement and stress fields in the V-notch frame as:

       u DIC t u DIC n        = Q        u DIC x u DIC y        (10)
and

    σ DIC tt σ DIC tn σ DIC tn σ DIC nn     = Q     σ DIC xx σ DIC xy σ DIC xy σ DIC yy     Q T (11)
where

Q =     cos γ sin γ -sin γ cos γ     (12)
Q is the rotation matrix between the DIC frame and the V-notch frame. In case of an isotropic homogeneous material (as PMMA in the present case), λ is independent of the elastic coefficients of the material. The exponent and the discretized values of the primal and dual modes (displacements and stresses) are tabulated to be used in the calculations described in the sequel.

The contour integral for GSIF calculation

SIF computation for a crack can be performed using, for instance, J-integral, which is not applicable in the case of other singularities such as V-notches, which require an alternative approach for GSIF calculation. The goal of this section is to propose a method to extract the GSIF k from U or more likely from a known approximation, a measured DIC field U DIC or a computed FE field U FE . To this aim, let us consider the following contour integral for two elastic solutions U and V fulfilling locally the homogeneous equilibrium equations and boundary conditions

Ψ(U , V ) = Γ [σ(U ).n.V -σ(V ).n.U ]dl, ( 13
)
where Γ is any contours encompassing the singular point and finishing on the neighbouring traction free (or vanishing displacements) faces, σ(U ) is the stress field associated to the displacement field U through Hooke's law, n is the normal to Γ pointing towards the origin (the singular point). It can be proved that the integral ( 13) is path independent.

Remark : Two important properties are recalled here (demonstrations can be found in [START_REF] Doitrand | Finite element procedure implementation of the coupled criterion: matched asymptotic and full finite element approaches[END_REF] and [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF]).

(i) If λ is a solution to the eigenvalue problem mentioned in Section 4.2, then -λ is a solution too. The condition λ > 0 is necessary to have a solution with a finite energy, nevertheless, -λ is a mathematical solution with its own eigenvector u -(ϕ). The functions u(ϕ) and u -(ϕ) are called respectively primal and dual modes. This property plays an essential role in the calculation of the 2D GSIFs.

(ii) If α and β are two singular exponents with their eigenmodes u α (ϕ) and u β (ϕ) then

β = -α =⇒ Ψ(r α u α , r β u β ) = 0 (14)
The above integral vanishes except for β = -α, i.e. if β is the dual exponent to α. Then, noting

U DIC = i k i r α i u α i (ϕ)
, we obtain:

Ψ(U DIC , r -α u -α ) = i k i Ψ(r α i u α i (ϕ), r -α u -α ) = k α Ψ(r α u α (ϕ), r -α u -α ) (15)
Then, from this property and (1), it comes:

k λ = Ψ(U DIC , r -λ u -) Ψ(r λ u, r -λ u -) (16)
Note that the denominator in ( 16) can be calculated once and for all independently of the contours as a consequence of [START_REF] Irwin | Analysis of stresses and strains near the end of a crack traversing a plate[END_REF].

Remark : The proposed approach for GSIF calculation can be applied for any singular point, including V-notches in the present analysis. The particular case of SIF calculation for a crack can thus be assessed employing the proposed approach.

Implementation of GSIF calculation from DIC fields

The computation of the GSIF is based on a post-processing of displacement and stress fields. On the one hand, in the case of FE calculations, it is straightforward to compute both quantities. On the other hand, DIC allows the measurement of displacement and strain full fields. It is thus necessary to assume a constitutive law in order to determine the missing stress field. Herein, we assume a linear elastic isotropic material behavior so as to compute the stresses σ from the strains ε using Hooke's law:

σ = E 1 + ν (ε + ν 1 -2ν tr(ε)I) ( 17
)
where E is the Young modulus, ν the Poisson ratio of the material and I the identity operator. Then, computing the integral [START_REF] Kotousov | Three-dimensional stress states at crack tip induced by shear and anti-plane loading[END_REF] requires the definition of a contour Γ. We choose to discretize a circular contour that ends on the V-notch flanks. Of course, it is ensured that all the integration points of the contour lies within the zone where the DIC fields are computed from full field measurements. Once the contour geometry is defined, the second step consists in extracting the quantities of interest along the contour (displacements and stresses). DIC displacements and stresses along the contour are obtained through interpolating the calculated quantities at the contour integration points and primal and dual modes are obtained from expressions given in Section 4.1.

Finally, the integral Ψ(u DIC , u J-) (where u J-denotes the dual displacement field, J = I for opening mode and J = II for shear mode) is computed on the discrete basis of the chosen contour, which allows determining the corresponding GSIF (Eqn. ( 16)). The proposed approach is used in next sections to compute the square hole V-notch GSIFs from displacement and stress fields obtained by FE calculations and from displacement and strain fields obtained by DIC (Section 6).

FE implementation of GSIF calculation

A 2D FE model of square hole specimens is set-up using Abaqus TM . Both plain stress and plain strain calculations have been performed, its influence on GSIF calculation is discussed in Section 6.

Linear elastic isotropic material behavior with E=2.7GPa Young's modulus and ν =0.39 Poisson's ratio is used. The material strength and fracture toughness, determined by inverse identification on circular hole specimens [START_REF] Doitrand | Experimental characterization of mixed mode failure of square hole specimens[END_REF], are σ c =43MPa and G c =90J/m 2 . As stated in introduction, the studied material is extruded PMMA, which presents lower strength and toughness than cast PMMA. It can be noted that the material under investigation exhibits a nonlinear behavior for stress levels larger than 20 MPa (Fig. 1). Nevertheless, the imposed stress levels do not exceed 15MPa for the specimen containing the smallest hole and are even smaller for larger holes. The stress is expected to be higher than 20 MPa in a localized zone in the vicinity of the V-notch singular point. Indeed, from theoretical considerations based on the coupled criterion, the size of this zone is of the same order of magnitude as EGc σ 2 c =130µm [START_REF] Martin | Understanding the tensile strength of ceramics in the presence of small critical flaws[END_REF]. In practice, we checked from numerical simulations that the size of the zone for which the stresses are larger than 20MPa for the maximum imposed stress levels lies between 100 and 200 µm for holes comprised between 1.2mm and 15.2mm. It is much smaller than the contour integral path radius. Therefore, It can be concluded that we compute the contour integral on a path for which the material behavior remains linear elastic.

A unit displacement is imposed on top and bottom faces of the specimen. Opening and shear mode GSIFs are computed according to the method presented in previous section for ten different circular contours with radius lying between c/10 and c/3. Maximum differences of 1% are obtained using different contours, which demonstrates the path independence of the integral [START_REF] Griffith | The phenomenon of rupture and flow in solids[END_REF] and thus the GSIF path independence [START_REF] Kotousov | Three-dimensional stress states at crack tip induced by shear and anti-plane loading[END_REF]. Fig. 5a shows the opening (k I ) and shear (k II ) mode GSIF variations as a function of the square hole size obtained numerically as well as the mode mixity defined as

ψ = tan -1 ( k II k I L λ II -λ I mat ) [6]
, where the characteristic length L mat is defined as

L mat = EGc σ 2 c
. It can be seen that a mixed mode configuration is obtained whatever the hole size. Moreover, the mode mixity is decreasing with increasing hole size, which is consistent with experimental observations of crack deflection with respect to the V-notch bisector. As explained in Section 3, for some specimens the initiation of the second crack occurred slightly after the first one. The presence of a first crack may have an influence on second crack initiation. As a matter of example, and since we did not quantify the first crack length when the second crack initiates, we study numerically the extreme case of a second crack initiating when the first crack already fully propagated to the specimen side (Fig. 5b-c).

Fig. 5b-c shows the GSIF and mode mixity variations as a function of the square hole side. It is clearly observed that the presence of a first fully propagated crack tends to decrease the mode mixity at the opposite corner, hence resulting in a lower crack deflection. Note that the experiments actually corresponds to an intermediate situation between the simultaneous initiation of both cracks and the initiation of the second crack when the first crack has fully propagated.

The coupled criterion

The criterion

The coupled criterion (CC) aims at filling the gap of crack nucleation prediction in brittle material [START_REF] Leguillon | Strength or toughness ? A criterion for crack onset at a notch[END_REF], which is not covered by classical Linear Elastic Fracture Mechanics [START_REF] Griffith | The phenomenon of rupture and flow in solids[END_REF][START_REF] Irwin | Analysis of stresses and strains near the end of a crack traversing a plate[END_REF]. This approach, which allows the determination of the initiation loading level, crack length and crack angle, is briefly recalled here. A detailed description of this approach for mixed mode crack initiation can be found in [START_REF] Sapora | A Finite Fracture Mechanics approach to V-notched elements subjected to mixed-mode loading[END_REF][START_REF] Doitrand | Experimental characterization of mixed mode failure of square hole specimens[END_REF]. Basically, it lies on the simultaneous fulfilment of a stress and an energy conditions. The former states that prior to crack initiation, the tensile stress acting on the presupposed crack plane σ nn must overcome the material strength σ c σ nn (r) ≥ σ c for any points along the presupposed crack path. [START_REF] Labossiere | Stress intensities at interface corners in anisotropic bimaterials[END_REF] The latter requires that the potential energy change due to crack nucleation (-∆W = W (0) -W (l),

where l is the crack length) is larger than the energy required to create a crack surface G c l (written here in 2D without ambiguity)

-∆W l ≥ G c . (19)
The CC can be applied either through the matched asymptotic (MA) expansions method which is quasi-analytic or through a full finite element (FFE) approach. Both presented in detailed for instance in [START_REF] Doitrand | Finite element procedure implementation of the coupled criterion: matched asymptotic and full finite element approaches[END_REF] and recalled in the sequel are equivalent provided the initiation crack length is small with respect to the other structural characteristic lengths (this is checked afterward). Under small deformations and linear elastic framework, the stress and the potential energy are respectively proportional to the applied load and the square of applied load, which allows both conditions to be combined into a single equation. Solving this equation allows computing the initiation crack extension and loading level.

Under mixed mode loading, the crack angle is a priori unknown, thus this procedure can be repeated for several crack angles. The initiation crack angle is obtained as the one minimizing the loading level for crack initiation, this minimum loading actually being the critical loading level.

The full finite element (FFE) approach

It is the most straightforward way to implement the CC and there is no restriction on its use. A first computation is carried out without a crack in order to determine the tensile stress acting on the presupposed crack paths. The plural is used because the crack direction is not known initially, its determination requires applying the CC for different directions and selecting the one minimizing the load at initiation. Next, for each crack angle, the FE implementation of the CC requires to release double nodes (sometimes called unbuttoning) in order to build the function W (l) and as a consequence to check the inequality [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF]. Therefore, this procedure requires several linear FE calculations to be executed. Nevertheless, it can be noted that these calculations can be carried out simultaneously so that the approach remains efficient numerically.

The matched asymptotic (MA) approach

This approach, described for instance in [START_REF] Doitrand | Finite element procedure implementation of the coupled criterion: matched asymptotic and full finite element approaches[END_REF][START_REF] Leguillon | Strength or toughness ? A criterion for crack onset at a notch[END_REF][START_REF] Leguillon | Size effects for crack initiation at blunt notches or cavities[END_REF], can be carried out under the assumption that the crack length at initiation remains small compared to any dimensions of the structure (here the square hole size and the ligament width). It is based on a two scale modelling under the form of outer (far field) and inner (near field) expansions that must match in an intermediate area. Outer means out of a vicinity of the notch root and inner inside this vicinity. The elastic solution involving a crack with length l is expanded in terms of this crack extension length

U l (x 1 , x 2 ) = U 0 (x 1 , x 2 ) + f 1 (l)U 1 (x 1 , x 2 ) + ... (20)
where U 0 is the solution to the unperturbed problem (i.e. without crack) and f 1 (l) → 0 as l → 0. Note that the terms U 0 , U 1 , etc. of the expansion are defined on a simplified domain where the small crack is not visible.

Similarly, after having stretched the space variables y i = x i /l and considered the limit as l → 0, the inner expansion takes the form

U l (x 1 , x 2 ) = U l (ly 1 , ly 2 ) = F 0 (l)V 0 (y 1 , y 2 ) + F 1 (l)V 1 (y 1 , y 2 ) + ... (21)
where F 1 (l)/F 0 (l) → 0 as l → 0. Parallel to the outer expansion, the terms of the inner expansion are defined on another kind of simplified domain, unbounded and where the remote boundaries have disappear. The matching conditions impose to this expansion to match at infinity with the behavior of the outer expansion when approaching 0 defined at the leading order by the Williams expansion (1) of U 0 , this leads to [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF] U l (ly 1 , ly

2 ) = C + k l λ V 1 (y 1 , y 2 ) + ... (22)
Using a superposition principle, it is shown that V 1 is solution to a well-posed problem, it is independent of the applied load and the global geometry of the specimen.

Thus, two descriptions are available prior to (l = 0) and following crack initiation (l > 0) and the change in potential energy can be calculated

∆W = A k l 2λ + ... ( 23
)
where A is a scaling coefficient that depends on the crack direction (but not its length). It can be calculated once and for all for each crack direction. Eqn. [START_REF] Leguillon | Strength or toughness ? A criterion for crack onset at a notch[END_REF] together with the stress condition derived from (1) and the Hooke's law [START_REF] Labossiere | Calculation of stress intensities at sharp notches in anisotropic media[END_REF] give the crack length at initiation l c and the critical value

of the GSIF k c l c = G c A 1 σ 2 c and k c = G c A 1-λ σ 2λ-1 c (24)
The selected crack direction corresponds to the one maximizing A, i.e. minimizing k c .

Results and discussion

The method presented in Section 4 is now applied to the experimental determination of GSIFs from displacement and strain fields measured by DIC.

GSIF variation during the test

A practical interest of the proposed approach is the ability to compute the GSIF variations during the test. minimum and maximum GSIF obtained on 10 different contours are presented. It can be observed that opening and shear mode GSIFs vary in an almost linear way, the shear mode sign changing depending on the square hole corner. As shown in Fig. 6, it can be seen that the uncertainties on GSIFs calculation are larger for smaller square hole size, and also larger for shear than opening mode. Moreover, from a sufficiently high load level (σ 0 = 2MPa for c = 6.16mm and σ 0 = 5MPa for c = 1.6mm), the mode mixity reaches a constant value until crack initiation. As explained in the previous section, we are not able to quantify the GSIF change in the case where one cracks initiates slightly before the other, which could possibly be undertaken using a rapid camera in order to catch the very first crack initiation.

Critical GSIF

From Figs. 7 and8, it can be seen that at least one of the initiated cracks correspond to the larger opening mode GSIF among the four square corners, which supports crack initiation prediction based on GSIF. Fig. 9 shows the failure opening and shear mode GSIFs as well as the mode mixity as a function of the square hole size obtained experimentally and numerically (employing either plane strain or plane stress conditions) for the loading level predicting using either FFE or MA approaches of the CC (the reader is referred to Section 5, [START_REF] Cornetti | Mode mixity and size effect in V-notched structures[END_REF] and [START_REF] Doitrand | Experimental characterization of mixed mode failure of square hole specimens[END_REF] for details about the CC theory and numerical implementation in the case of square holes).

First, it can be observed that both MA and FFE approaches of the CC give similar opening and shear failure GSIF predictions, which could be expected since the material characteristic length L mat is small with respect to any other structural dimensions [START_REF] Doitrand | Finite element procedure implementation of the coupled criterion: matched asymptotic and full finite element approaches[END_REF]. It can be noted that both opening and shear mode GSIFs as well as mode mixity predicted employing plane stress conditions are slightly smaller than those obtained with plane strain conditions. A slight better agreement with experimental data is obtained employing plain stress conditions, which can be explained by the fact that the GSIF are computed based on displacement and strain fields obtained on the specimen surface. Qualitative trends of increasing opening mode and decreasing shear mode GSIF are well captured by the numerical predictions, a satisfying agreement between numerical and experimental GSIF being obtained.

Note that the agreement between the failure GSIFs determined experimentally and computed numerically was not expected to be perfect since the proposed approach relies on several assumptions: (i) The FE (either FFE or MA) approach assumes a perfectly sharp notch whereas this cannot be physically achieved, there is a small rounding and this increases the apparent critical GSIFs. Following an identical procedure to that developed in [START_REF] Leguillon | Determination of the length of a short crack at a v-notch from a full field measurement[END_REF], it can be shown that the apparent GSIFs k 1 and k 2 vary respectively with ρ 2λ 1 and ρ 2λ 2 where ρ is the root radius (assumed to be small) that blunts the V-notch.

(ii) More difficult to control, flaws due to machining inevitably exist at the root of the notch which tend, unlike the previous case, to decrease the apparent critical GSIFs. The laser cut also induces a thermal treatment at the very surface of the specimen hole, which may affect the material properties in a localized layer around the hole.

(iii) We considered 2D displacements and strain fields obtained by DIC on the surface of the specimens. Therefore, a source of uncertainties comes from the fact that the plane in which the displacement fields are measured is the specimen surface whereas the 2D FE model is more inline with the middle plane of the specimen [START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF]. In practice, we observed the specimen surface where the singularity of the V-notch tip is slightly different than along the notch edge (this one actually corresponds to the 2D singularity).

It is the singularity of a V-notch reaching a free surface. Computing the exponents and the associated modes of such a singularity can be carried out, but would require 3D DIC fields in order to compute experimentally the corresponding GSIFs (using the approach described in [START_REF] Doitrand | Computation of generalized stress intensity factors of 3D singularities[END_REF][START_REF] Leguillon | Computation of 3D-singularities in elasticity[END_REF]). The 2D contour integral allows computing 2D opening and shear mode GSIFs from 2D displacement fields. It has to be noted that there also exist an out of plane mode which magnitude may be significant for a long crack loaded in shear [START_REF] Kotousov | Three-dimensional stress states at crack tip induced by shear and anti-plane loading[END_REF] and depends on the plate thickness [START_REF] Kotousov | Effect of plate thickness on stress state at sharp notches and the strength paradox of thick plates[END_REF], which we do not consider since it has only a weak influence on the initiation crack direction. This mode corresponds to mode III in the case of a crack, which leads to a fragmentation of its front [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF][START_REF] Pham | On the growth of cracks under mixed-mode I + III loading[END_REF] which orientation is determined by mode I and II [START_REF] Buchholz | Fracture analyses and experimental results of crack growth under general mixed mode loading conditions[END_REF]. At interior points along the straight V-notch (i.e. except at the two surface ends), the leading terms of the William's expansion are the 2D plane strain singularities [START_REF] Apel | Edge singularities and structure of the 3-D Williams expansion[END_REF][START_REF] Yosibash | Edge stress intensity functions in polyhedral domains and their extraction by a quasidual function method[END_REF]. It is not the case at the two surface ends where the singularity is not the same as for interior points. Indeed, for the 3D corner, there are two singular modes: a symmetric one that extend the 2D symmetric (opening) mode and a second mode which is a mixture of the two 2D shear modes. 2D and 3D modes have slightly different singular exponents (for instance for opening mode, the singularity exponent is 0.545 in 2D and 0.669 in 3D). Moreover, it can be shown that the trace of the 3D modes on the free surface is similar to the 2D modes. Therefore, the error made in exploiting information extracted from the surface remains acceptable. We also checked based on 3D FE calculations that there is a small difference (less than 3%) between the 2D GSIFs computed based on the specimen surface or middle plane displacement fields. Fig. 10 depicts the initiation crack deflection with respect to the V-notch bisector measured experimentally and obtained numerically using either MA or FFE approach. Similarly to the GSIF prediction, the crack deflection predictions obtained using either one or the other approach are similar and are in good agreement with the crack deflection measured experimentally. The crack deflection of a second crack in presence of a first crack that would have fully propagated to the specimen edge predicted using FFE approach of the CC is also reported. As explained in a previous section, it is obviously an extreme case that does not occur experimentally, which however gives a lower bound of the mode mixity and of the crack deflection with respect to the V-notch bisector.

Conclusion

Mixed mode failure is observed in PMMA specimens containing square holes under uniaxial tensile loading. The larger the hole size, the smaller the mode mixity defined as a normalized ratio between the GSIFs of the shear and opening modes. By means of DIC, we compute the displacement and strain fields around the square hole corners and determine the corresponding stress fields assuming a linear elastic isotropic material behavior which properties are also determined using DIC. By means of a path independent integral, the square hole corner GSIFs are determined experimentally without relying on any FE models. The uncertainty on the GSIF calculation is lower than 20 % in most cases, except for too small holes for which the uncertainty linked to the displacement field measurements using DIC are higher. Opening and shear mode GSIFs variations as a function of the imposed stress are recorded during the test, which allows determining the mode mixity variation on the fly, it reaches a constant value from a certain imposed loading level. Crack initiation occurs at one of the corners exhibiting the largest GSIFs. Numerical predictions of the GSIF based on the same path independent integral and on crack initiation prediction obtained using the CC are in reasonable agreement with the GSIFs determined experimentally. Future work will cover the extension of the proposed approach to case of 3D singularities, which will require DIC fields obtained by Digital Volume Correlation in order to compute the GSIF of a 3D singularity experimentally.

A

  : scaling coefficient c : square hole side C : rigid translation constant E : Young's modulus G c : fracture toughness I : identity operator k : generalized stress intensity factor k I , k II : opening and shear mode generalized stress intensity factor l : crack length L mat : characteristic length n : normal to the integral contour P : polar to cartesian coordinate transformation matrix Q : DIC to V-notch frame transformation matrix r : polar coordinate u : singular mode u + , u -: primal and dual modes

:

  DIC displacement vector in the V-notch frame U , V : displacement vectors U DIC : DIC displacement vector U F E : FE displacement vector w : specimen width W : potential energy γ : angle between DIC and V-notch frames Γ : integral contour λ, α, β : singularity exponent λ I : opening mode singularity exponent λ II : shear mode singularity exponent λ L , µ L : Lamé's coefficients ν : Poisson's ratio ϕ : polar coordinate ψ : mode mixity Ψ : contour integral ω : V-notch angle ρ : notch radius σ c : strength σ 0 : prescribed stress σ I , σ II : opening and shear mode stress in the V-notch axis system σ I rr , σ I ϕϕ , σ I rϕ : opening mode stress components in the V-notch axis system σ J tt , σ J tn , σ J nn : cartesian stress components in the V-notch axis system σ II rr , σ II ϕϕ , σ II rϕ : shear mode stress components in the V-notch axis system σ DIC xx , σ DIC xy , σ DIC yy : cartesian stress components in the global frame σ DIC tt , σ DIC tn , σ DIC nn : DIC stress components in the V-notch frame σ : stress tensor σ yy : stress tensor component ε : strain tensor ε xx , ε yy : strain tensor components θ c : crack deflection DIC : Digital Image Correlation FE : Finite Element FFE : Full Finite Element FFM : Finite Fracture Mechanics GSIF : Generalized Stress Intensity Factors MA : Matched Asymptotic SIF : Stress Intensity Factors

Figure 1 :Figure 2 :

 12 Figure 1: Stress-Strain curves of the material under investigation.

Figure 3 :

 3 Figure 3: Photographs of the tested specimens with square hole size from (a) 15.2mm to (j) 1.23mm after failure.

Fig. 4 Figure 4 :

 44 Figure 4: Photograph of the square hole specimen indicating the local V-notch frames (blue, red, black and green) and the DIC frame (purple).

Figure 5 :

 5 Figure 5: GSIF variations obtained from FE calculations (for a +1mm imposed displacement on the top and bottom edges of the specimen) for the top right corner as a function of the square hole size in the presence of (a) no crack or (b-c) one crack located at the (b) top left or (c) bottom left corner.

  Figs. 7a-c and 8a-c show the opening and shear mode and mode mixity variations at each corner as a function of the prescribed stress σ 0 for c = 6.16 mm and c = 1.6 mm hole size. For each

Figure 7 :

 7 Figure 7: (a) Opening (k I ) and (b) shear (k II ) mode GSIFs and (c) mode mixity ψ at the four square hole corners (c=6.16mm) as a function of the prescribed stress σ 0 . Photograph of the square hole (d) before (highlighting the corners with corresponding colors in (a-c)) and (e) after failure.

  10mm

Figure 8 :

 8 Figure 8: (a) Opening (k I ) and (b) shear (k II ) mode GSIFs and (c) mode mixity ψ of the four square hole corners (c=1.6mm) as a function of the imposed stress σ 0 . Photograph of the square hole (d) before (highlighting the corresponding corners in (a-c)) and (e) after failure.

Figure 9 :

 9 Figure 9: Failure (a) Opening (k I ) and (b) shear (k II ) mode GSIFs and (c) mode mixity ψ of the four square hole corners as a function of the square hole side.

Figure 10 :

 10 Figure 10: Crack deflection with respect to the square hole V-notch bisector as a function of the square hole size measured experimentally (error bars) and obtained numerically in presence of (red square) or without (blue circles) an initial crack.

  Vicentini et al. [43] determined the opening and shear mode GSIFs at a bimaterial closed corner in Brazilian disk specimens by means of FE calculations based on the failure load measured experimentally. It can be noted that the previously cited work proposed indirect GSIF estimates based on experimental data, since these approaches rely on a FE model. Their main limitation is the representativeness of the FE model with respect to the experimental test in terms of specimen geometry or boundary

conditions for instance. Moreover, these works only concerned pure opening mode configurations but such approaches would not allow identifying the GSIFs corresponding to different modes in a mixed mode loading configuration.

Table 1 :

 1 Dimensions of the tested specimens.

	2 and 3 depict the specimens with different square hole sizes respectively before the tests and after
	failure. The crack orientation, measured from microscope observations of the square hole V-notches
	after failure, indicates that failure occurs under mixed mode loading. Indeed, in the case of a pure
	Square hole side c (+0.05mm) 1.23 1.60 2.35 2.70 3.60 5.05 6.16 9.23 12.28 15.2
	Specimen width w (+0.2 mm) 9.80 14.80 19.80 24.80 29.80 29.80 29.80 49.80 49.80 49.80

Integral path independence

The GSIFs are calculated for all square hole configurations tested experimentally using ten different circular contours for a prescribed stress σ 0 =5 MPa (the ratio between the applied force and the initial surface of the specimen cross section). Fig. 6 depicts the relative variation (i.e. the GSIF range to mean value ratio obtained for the ten contours) of the calculated GSIF at each corner as a function of the square hole size. It can be observed that for square hole sizes c ≥ 2 mm, the opening mode GSIF is computed within less than 20% uncertainty in the worst case. The uncertainty on shear mode GSIF calculation is slightly larger than for the opening mode. For smaller holes, the uncertainty drastically increases, which can be explained by the fact that similar speckle patterns were used for all specimens.

The smaller the hole, the smaller the DIC mesh size and hence the higher the uncertainties in the calculation of the displacement fields [START_REF] Besnard | Finite-element displacement fields analysis from digital images: application to Portevin-Le Chatelier bands[END_REF].