V Kosin

A Fau

C Jailin

F Hild

T Wick

Parameter identification of a phase-field fracture model using integrated digital image correlation

Keywords: Phase-field fracture modeling, digital image correlation (DIC), Conditioning, Sensitivity analysis

Phase-field fracture (PFF) modeling is a popular approach to model and simulate fracture processes in solids. Accurate material parameters and boundary conditions are of utmost importance to ensure a good prediction quality of numerical simulations. In this work, an Integrated Digital Image Correlation (IDIC) algorithm is proposed to calibrate boundary conditions, Poisson's ratio, fracture energy and internal length, all at once, by using the phase-field model itself and images of a deforming sample. The presented approach is applied to virtual experiments mimicking a single edge notched shear test and implemented in the opensource deal.II-based software pfm-cracks and the digital image correlation library Correli 3.2. The reliability of the results is investigated for different levels of acquisition noise, thereby demonstrating high robustness and accuracy for a wide range of noise levels. The conditioning of the problem is analyzed via sensitivity fields for all parameters and the eigendecomposition of the Hessian matrix used in the IDIC algorithm.

Introduction

Fracture assessment is a current topic in engineering and applied mathematics. To model fracture, various approaches are known such as boundary finite elements [START_REF] Chang | A boundary element method for two dimensional linear elastic fracture analysis[END_REF], displacement discontinuity methods [START_REF] Crouch | Solution of plane elasticity problems by the displacement discontinuity method. i. infinite body solution[END_REF], cohesive-zone models [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF], generalized/extended finite elements [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Fries | The extended/generalized finite element method: An overview of the method and its applications[END_REF][START_REF] Babuska | Stable Generalized Finite Element Method (SGFEM)[END_REF], peridynamics [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF] or (regularized) variational phase-field fracture models [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF][START_REF] Kuhn | A continuum phase field model for fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF]. In this work, the latter ones are of interest, namely, phase-field approaches, which share close relationships with damage models [START_REF] Lemaitre | Mechanics of solid materials[END_REF][START_REF] Lemaitre | A Course on Damage Mechanics[END_REF][START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF][START_REF] Desmorat | Nonlocal models with damage-dependent interactions motivated by internal time[END_REF][START_REF] De Borst | Gradient damage vs phase-field approaches for fracture: Similarities and differences[END_REF]. Based on energy formulations of brittle fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and their numerical implementation [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], this approach has been further developed and applied in numerous situations [START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Wu | Phase field modelling of fracture[END_REF][START_REF] Bourdin | Past and present of variational fracture[END_REF][START_REF] Diehl | A comparative review of peridynamics and phase-field models for engineering fracture mechanics[END_REF][START_REF] Wick | Multiphysics Phase-Field Fracture: Modeling, Adaptive Discretizations, and Solvers[END_REF]. It is intriguing because the crack path, crack direction, and crack tip velocity are obtained within the model in an implicit fashion without explicit need to compute stress intensity factors and crack front positions. Moreover, it is relatively easy to apply (except the computational cost) to two-and three-dimensional spatial settings. Further, crack initiation is included (see, for instance, L-shaped panel tests [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF], where further theoretical and computational analyses about its functionality are carried out [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF][START_REF] Kumar | Revisiting nucleation in the phase-field approach to brittle fracture[END_REF]).

To date, in most studies, the primary focus is on so-called forward models in which the phase-field framework with given model, numerical, and material parameters is employed to compute the material state, which is theoretically, computationally, or experimentally analyzed. Here, one important goal is to use numerical simulations to predict the behavior of components within developments of new materials or work tools. Much less work has been carried in inverse modeling and optimization in which unknown parameters or boundary conditions are to be computed by comparing the numerical simulations to experimental data sets. Some recent work on optimal control with phase-field fracture shows how the phase-field model is driven into some desired states by adapting boundary conditions [START_REF] Neitzel | An optimal control problem governed by a regularized phase-field fracture propagation model[END_REF][START_REF] Neitzel | An optimal control problem governed by a regularized phase-field fracture propagation model. part II: The regularization limit[END_REF][START_REF] Mohammadi | A priori error estimates for a linearized fracture control problem[END_REF][START_REF] Khimin | Space-time formulation, discretization, and computational performance studies for phase-field fracture optimal control problems[END_REF][START_REF] Khimin | Space-time mixed system formulation of phase-field fracture optimal control problems[END_REF]. Topology optimization was designed [START_REF] Desai | Topology optimization of structures undergoing brittle fracture[END_REF], and stochastic phase-field modeling [START_REF] Gerasimov | Stochastic phase-field modeling of brittle fracture: Computing multiple crack patterns and their probabilities[END_REF] was discussed. Parameter identification based on lower-scale simulations (e.g., via molecular dynamics), which are then employed within phase-field fracture models, was proposed in Refs. [START_REF] Hansen-Dörr | Combined molecular dynamics and phase-field modelling of crack propagation in defective graphene[END_REF][START_REF] Patil | A comparative molecular dynamicsphase-field modeling approach to brittle fracture[END_REF]. The closest studies in their objectives to the present one, by addressing parameter identification, are found in Refs. [START_REF] Khodadadian | A bayesian estimation method for variational phase-field fracture problems[END_REF][START_REF] Wu | Parameter identification for phase-field modeling of fracture: a Bayesian approach with sampling-free update[END_REF][START_REF] Noii | Bayesian inversion for unified ductile phase-field fracture[END_REF][START_REF] Noii | Bayesian inversion for anisotropic hydraulic phase-field fracture[END_REF]. However, therein Bayesian inversion [START_REF] Smith | Uncertainty Quantification: Theory, Implementation, and Applications[END_REF] techniques were employed. The main novelty of the current work is to use digital image correlation to estimate unknown parameters, and to the best of the authors' knowledge the integrated combination of images and phase-field simulations has not yet been proposed in the literature.

Phase-field models have been probed against experimentally measured kinematic fields on a qualitative basis essentially. First, simulated strain fields were compared to experimental strain fields measured via digital image correlation [START_REF] Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF][START_REF] Na | Effects of spatial heterogeneity and material anisotropy on the fracture pattern and macroscopic effective toughness of mancos shale in brazilian tests[END_REF][START_REF] Donnini | Uniaxial tensile behavior of ultra-high performance fiberreinforced concrete (uhpfrc): Experiments and modeling[END_REF][START_REF] Hun | Desiccation cracking of heterogeneous clayey soil: Experiments, modeling and simulations[END_REF][START_REF] Xu | Phase-field modeling of crack growth and interaction in rock[END_REF]. Similarly, displacement fields were used instead [START_REF] Wu | Phase-field simulation of interactive mixed-mode fracture tests on cement mortar with full-field displacement boundary conditions[END_REF]. Second, crack patterns were also compared when imaged via X-ray computed tomography [START_REF] Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF][START_REF] Tsitova | Identification of microscale fracture models for mortar with in-situ tests[END_REF] when 3D simulations were driven by measured boundary conditions. In all the aforementioned references, the internal length was not calibrated, even though it is believed to be a material parameter (see [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF] with more references cited therein). On the other hand, from the mathematical viewpoint the length scale was initially introduced as a regularization parameter [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] as in nonlocal damage models. The internal length was identified by using force/strain data for bone [START_REF] Hug | Predicting fracture in the proximal humerus using phase field models[END_REF]. The full strain validations were performed from a purely statistical point of view (and not experimental field against simulated field).

In the context of nonlocal damage, several routes have been considered to calibrate the internal (or also called characteristic) length. It was first estimated experimentally by an energetic approach as the ratio of the fracture energy to the energy dissipated per unit volume of the material [START_REF] Lemaitre | Damage measurements[END_REF][START_REF] Bažant | Measurement of characteristic length of nonlocal continuum[END_REF]. An identification based on several load deflection curves from three point bending tests was proposed [START_REF] Bellégo | Calibration of nonlocal damage model from size effect tests[END_REF]. However, it was shown that the identification of the internal length based on such global curves did not lead to a unique solution [START_REF] Iacono | Estimation of model parameters in nonlocal damage theories by inverse analysis techniques[END_REF]. It was suggested to calibrate a non-local damage model from discrete models [START_REF] Xenos | Calibration of nonlocal models for tensile fracture in quasi-brittle heterogeneous materials[END_REF]. An identification of the damage parameters based on digital image correlation (DIC) data was also considered [START_REF] Geers | Mixed numerical-experimental identification of non-local characteristics of randomfibre-reinforced composites[END_REF]. The internal length was identified without coupling DIC and the identification procedure. The DIC measurements included a restricted area of the sample with a relatively small number of pixels.

In the following, it is proposed to investigate the feasibility of calibrating the internal length of a phasefield model by using displacement field measurements. This parameter is the most delicate [START_REF] Jailin | Virtual hybrid test control of sinuous crack[END_REF]. Other material parameters will also be sought (i.e., Poisson's ratio and fracture energy). To make the analysis even more challenging, the boundary conditions will also be considered unknown. Various methods exist to calibrate material parameters from full-field measurements [START_REF]Full-Field Measurements and Identification in Solid Mechanics[END_REF]. In this study, the displacement fields are determined via digital image correlation [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF][START_REF] Sutton | Computer vision-based, noncontacting deformation measurements in mechanics: A generational transformation[END_REF]. Series of images of two synthetic tests, generated with multiple acquisition noise levels, were constructed by using simulated displacement fields. Instead of following a two-step procedure, namely, first measuring displacement fields via DIC and then post-processing them to extract the unknown parameters, integrated digital image correlation (IDIC) was selected [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF][START_REF] Mathieu | Estimation of elastoplastic parameters via weighted FEMU and integrated-DIC[END_REF].

It performs registrations using mechanically admissible displacement fields (i.e., the two steps are merged into a unique minimization). This approach enables, for instance, meshes to be tailored as fine as needed for numerical reasons [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF][START_REF] Lindner | On the evaluation of stress triaxiality fields in a notched titanium alloy sample via integrated DIC[END_REF][START_REF] Hild | Toward 4d mechanical correlation[END_REF] instead of controlled by measurement uncertainties.

The outline of this paper is as follows. In Section 2, the overall notations of this paper and integrated digital image correlation are introduced. Section 3 is devoted to the governing phase-field model used herein. Section 4 introduces the integrated computational framework applied to applications for parameter identification.

Integrated digital image correlation

Let Ω ⊂ R 2 be the region of interest and I := (0, T ) be the time interval of interest, where T > 0 is the end time value. Moreover, let Ω × I be the so-called space-time cylinder.

Digital image correlation and space-time framework

Digital image correlation has proven to be a powerful tool to measure displacement fields [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF][START_REF] Hild | Digital image correlation[END_REF][START_REF] Sutton | Computer vision-based, noncontacting deformation measurements in mechanics: A generational transformation[END_REF].

Assuming brightness conservation, one can register an image I 0 : Ω → R with a series of images I : Ω×I → R.

Then, there is a displacement field u : Ω × I → R 2 such that

I 0 (x) = I(x + u(x, t), t) =: I u (x, t), ∀(x, t) ∈ Ω × I. (1)
The underlying problem to solve is constructed by minimizing the squared differences of the reference image I 0 and the backdeformed series of images I u over the space-time cylinder Ω × I.

Problem 1 (DIC). Find a displacement field u ∈ V := L 2 (I; L 2 (Ω)) such that

u = arg min v∈V J 1 (v), (2)
with

J 1 (v) := ∥I v -I 0 ∥ V and ∥ • ∥ 2 V := I Ω • 2 dx dt. ( 3
)
Remark 1. In V there are no boundary conditions prescribed for u. We notice that this is usual in DIC since the boundary is considered in the corresponding numerical solution algorithm, and consequently the final solution is unique. A detailed discussion about the well-posedness of the static problem and the use of

L 2 (Ω)
in space is given in Ref. [START_REF] Fedele | Global 2d digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation[END_REF].

Remark 2. Regular DIC algorithms are not formulated in a space-time framework, such that the displacement field is computed for each image I t of the series of images I separately [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF]. But the space-time approach has proven to be more robust with respect to noise and low sensitivities, not only for DIC but also for integrated digital image correlation (IDIC, see next sub-section), which is expected to be important for the calibration of the internal length. A more detailed discussion about the advantages of space-time formulations for DIC can be found in Ref. [START_REF] Besnard | Analysis of image series through digital image correlation[END_REF] and for IDIC in Refs. [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted FEMU and integrated-DIC[END_REF][START_REF] Neggers | Time-resolved integrated digital image correlation[END_REF].

Integrated identification

IDIC is a method to identify parameters of some mechanical models directly from images [START_REF] Hild | Digital image correlation: From measurement to identification of elastic properties -a review[END_REF][START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF]. Based on Problem 1, one is again minimizing the squared differences between a reference image and a backdeformed series of images, but the unknown u is made explicitly dependent on the sought parameters. Let P be a vector space of possible parameterizations and M : P → V be a model that computes a displacement field for the parameterization p ∈ P, then one can substitute the displacement by the model response,

u = M(p). (4)
Therefore the new problem reads as follows:

Problem 2 (IDIC). Let I 0 be a reference image, I t a series of images and M a model. Find a parameterization p ∈ P such that p = arg min q∈P J 2 (q) (5)

with

J 2 (q) := ∥I M(q) -I 0 ∥ 2 V . (6)
Assuming M(p) to be Fréchet differentiable with derivative ∂M(p)(q), then a solution to Problem 2 needs to fulfill the following first-order optimality condition

∂J 2 (p)(q) := I Ω (∂M(p)(q) • ∇I M(p) )(I M(p) -I 0 ) dx dt = 0 ∀q ∈ P. (7)
Solving Equation ( 7) is a nonlinear problem, and Newton's method is proposed as the solution strategy.

Assuming that M(p) is twice Fréchet differentiable with the second-order derivative ∂ 2 M(p)(δp, q), the corresponding Hessian reads [63]

∂ 2 J 2 (p)(δp, q) = I Ω (∂ 2 M(p)(δp, q) • ∇I M(p) )(I M(p) -I 0 ) + (∂M(p)(δp) • ∇ 2 I M(p) • ∂M(p)(q))(I M(p) -I 0 ) + (∂M(p)(δp) • ∇I M(p) )(∂M(p)(q) • ∇I M(p) ) dx dt. ( 8
)
For the computation of both, ∂M(p)(q) and ∂ 2 M(p)(δp, q) a finite-difference scheme is introduced.

Since P is a vector space, there exists a basis {Φ p i } such that q = j a j Φ p j and δp = j b j Φ p j for any (q, δp) ∈ P 2 . Last, the derivatives are approximated as

∂M(p)(q) ≈ i M(p + sΦ p i ) -M(p) s a i =: M ′ s (p)(q), (9)
and

∂ 2 M(p)(δp, q) ≈ i,j M(p + sΦ p i + sΦ p j ) -M(p + sΦ p i ) -M(p + sΦ p j ) + M(p) s 2 a i b j =: M ′′ s (p)(δp, q), (10)
where s > 0 is a given fixed perturbation parameter. Assuming that P is n-dimensional, with n ∈ N being a finite number, M ′ needs n evaluations of M and M ′′ needs n 2 evaluations of M. In the case of the phase-field model, which is described in the next section, the computation of M ′′ is extremely expensive even for a low-dimensional parameter space. Since the term ∂ 2 M includes the residuals

ρ(x, t) := I M(p) (x, t) -I 0 (x), (11)
which tend to zero near the solution, it is neglected. Together with the finite-difference approximation, the new (approximated) Jacobian reads

A s (p)(q) := I Ω (M ′ s (p)(q) • ∇I M(p) )(I M(p) -I 0 ) dx dt, (12)
and the new (approximated) Hessian becomes

∂A s (p)(δp, q) = I Ω (∂M(p)(δp) • ∇ 2 I M(p) • ∂M(p)(q))(I M(p) -I 0 ) + (∂M(p)(δp) • ∇I M(p) )(∂M(p)(q) • ∇I M(p) ) dx dt. (13)
The final simplified Newton algorithm for Problem 2 including line-search iterations, is given by Algorithm 1.

Algorithm 1 Newton's method for IDIC 1: Choose an initial guess p (0) ∈ P, a basis {Φ P k } and s > 0 2: for i = 0, 1, 2, ... do ▷ Newton iteration

3:

Evaluate M(p i )(Φ P k ) for all k to be able to compute M ′ s (p)(q) for all q ∈ P 4:

Solve ∂A s (p (i) )(δp, q) = -A s (p (i) )(q) ∀q ∈ P (14) 5:
for l = 0, 1, ..., l max do ▷ Line-search iteration

6: λ = 0.5 l 7: p (i+1) ← p (i) + λδp 8: if J (p (i+1) ) < J (p (i) ) then 9:
Stop line-search iteration

if k ∥δa k /a (0) k ∥ < TOL with p (0) = (a (0) 1 , a (0)
2 , ...), δp = (δa 1 , δa 2 , ...), then 13:

Stop Newton iteration 14:

end if 15: end for Remark 3. Using again the arguments of vanishing residuals, one could also neglect the term including the Hessian matrix ∇ 2 I M(p) in Equations ( 8) and ( 13). This simplification then results in a Gauss-Newton scheme. The vanishing residuals are achieved by good initial guesses. One can get satisfactory initial guesses by multi-level approaches using for example coarse graining or blurring of the images [START_REF] Hild | Digital image correlation[END_REF]. A more detailed discussion can be found in Refs. [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF][START_REF] Fedele | Global 2d digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation[END_REF].

Remark 4. Instead of using a finite difference scheme for the computation of ∂M(p)(q) and ∂ 2 M(p)(δp, q) one can also use the direct differentiation of the residual equation of the model [START_REF] Grešovnik | Ageneral purpose computational shell for solving inverse and optimisation problems[END_REF].

Phase-field modeling of brittle fracture

In this section, a phase-field approach is introduced to model fracture. One starts from the energy level and obtains, by differentiation, the Euler-Lagrange equations in terms of a coupled variational inequality system.

Modeling

A phase-field model based on the variational formulation of brittle fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and more precisely the related regularization form [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] is briefly introduced to specify the model M used in Section 4.

Assuming small displacements and linear elasticity, the infinitesimal strain tensor is written as

ε := ε(u) = 1 2 (∇u + ∇u T ), (15)
and the elastic energy density function

Ψ 0 (ε) = µ ε : ε + λ 2 tr(ε • ε), (16)
where µ and λ are the Lamé parameters. A smoothed scalar-valued phase-field function

φ : Ω → [0, 1]
is introduced to model damage inside a linear elastic material. The value 0 indicates a fully damaged material, and 1 an undamaged material. Moreover, it is assumed that damage is an irreversible process (i.e., no healing is possible) [START_REF] Lemaitre | A Course on Damage Mechanics[END_REF], which is modeled by ∂ t φ ≤ 0 and known as the crack irreversibility constraint (i.e., an inequality constraint in time). To circumvent unrealistic damage patterns, an additive decomposition of the energy density is introduced [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] with a comparison to other models made in [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF].

Let ε(u) = i ϵ i n i ⊗ n i be the spectral decomposition of the strain tensor and ε ± (u) :

= i ⟨ϵ i ⟩ ± n i ⊗ n i with ⟨•⟩ ± := 1 2 (• ± | • |).
Then, the energy density is decomposed as

Ψ 0 (ε) = Ψ + (ε) + Ψ -(ε), (17)
where

Ψ + (ε) := µ ε + : ε + + λ 2 ⟨tr(ε)⟩ 2 + and Ψ -(ε) := µ ε -: ε -+ λ 2 ⟨tr(ε)⟩ 2 - ( 18
)
are the extensional and contractional energy densities, respectively. Then, the regularized problem is formulated as follows: Find u and φ for almost all times t such that (u, φ) = arg min

(v,ψ) E ℓ ((v, ψ)) and ∂ t φ ≤ 0, (19)
with

E ℓ ((u, φ)) := 1 2 Ω (1 -κ)φ 2 + κ Ψ + (ε(u)) + Ψ -(ε(u)) + G c 1 2ℓ (1 -φ) 2 + ℓ 2 |∇φ| 2 dx, ( 20
)
where G c is the fracture energy, ℓ the internal length, and κ a positive regularization parameter to counteract numerical issues (zero entries on the diagonal where φ = 0 in the discrete system matrix resulting into a singular, non-solvable, linear equation system) when φ tends to 0. Since one works here in a quasistatic regime without time derivatives in the governing equations, the time-continuous crack irreversibility constraint is written in incremental form using a difference-quotient approximation

∂ t φ ≈ φ(t n ) -φ(t n-1 ) t n -t n-1 (21)
with t n , t n-1 ∈ I. Consequently, the phase-field variable is from the following closed convex set

Y in := {φ n ∈ H 1 (Ω)| 0 ≤ φ n ≤ φ n-1 ≤ 1 a.e. in Ω}.
The solutions u(t n ) and φ(t n ) at the time t n are from now on denoted by u n and φ n . For the displacements, one uses the function space

X := {u ∈ H 1 (Ω)| u = u D on Γ ⊂ ∂Ω}.
Remark 5. The Dirichlet boundary conditions u D are part of the unknowns in Section 4 and are going to be computed via IDIC.

The resulting optimization problem is then summarized as follows:

Problem 3. Given φ 0 and for the incremental steps t n , with n = 1, ..., N , find u n ∈ X and φ n ∈ Y in , such that

(u n , φ n ) = arg min (v,ψ)∈X×Yin E ℓ ((v, ψ)) (22)
where E ℓ is defined in Equation [START_REF] Fedele | Global 2d digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation[END_REF].

100
The corresponding optimality condition is given by [START_REF] Wick | Multiphysics Phase-Field Fracture: Modeling, Adaptive Discretizations, and Solvers[END_REF][ Proposition 19]

∂E ℓ ((u n , φ n ))((v, ψ -φ n )) ≥ 0 ∀(v, ψ) ∈ X × Y in , (23)
where

∂E ℓ ((u, φ))((v, ψ)) = Ω (1 -κ)φ 2 + κ ∂Ψ + (ε(u))(ε(v)) + ∂Ψ -(ε(u))(ε(v)) + (1 -κ)φψ Ψ + (ε(u)) + G c - 1 ℓ (1 -φ)ψ + ℓ∇φ • ∇ψ dx (24)
and

∂Ψ + (ε(u))(ε(v)) = µ ε + (u) : ε(v) + λ 2 ⟨tr(ε(u))⟩ + tr(ε(v)), ∂Ψ -(ε(u))(ε(v)) = µ ε -(u) : ε(v) + λ 2 ⟨tr(ε(u))⟩ -tr(ε(v)). (25)
Unfortunately, the energy functional E ℓ ((u, φ)), as defined in Equation ( 20), is not convex simultaneously in each variable such that the numerical solution of the fully monolithic problem becomes difficult. As proposed in Ref. [START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF], the solution variable in the term (1

-κ)φ 2 + κ ∂Ψ + (ε(u))(ε(v)) is approximated by a linear extrapolation in time φ ≈ φ(φ n-1 , φ n-2 ) := φ n-2 + t -t n-2 t n-1 -t n-2 (φ n-1 -φ n-2 ), ( 26
)
where φ n-1 and φ n-2 denote the solutions at the previous times t n-1 and t n-2 , respectively. In order to still have a monolithically-coupled variational formulation, φ is not approximated in the energy functional [START_REF] Fedele | Global 2d digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation[END_REF], but in the optimality condition [START_REF] Geers | Mixed numerical-experimental identification of non-local characteristics of randomfibre-reinforced composites[END_REF], such that ∂E ℓ is approximated by the semilinear form

B ℓ ((u, φ))((v, ψ)) := Ω (1 -κ) φ 2 + κ ∂Ψ + (ε(u))(ε(v)) + ∂Ψ -(ε(u))(ε(v)) + (1 -κ)φψ Ψ + (ε(u)) + G c - 1 ℓ (1 -φ)ψ + ℓ∇φ • ∇ψ dx ( 27
)
and the new optimality condition reads

∂B ℓ ((u n , φ n ))((v, ψ -φ n )) ≥ 0 ∀(v, ψ) ∈ X × Y in . (28)
Equation ( 28) is solved by Newton's method. A detailed discussion about the solution algorithm including the derivation and computation of the Jacobian can be found in Ref. [START_REF] Jodlbauer | Matrix-free multigrid solvers for phase-field fracture problems[END_REF][Section 3].

Discretization

As a next step the solution of the phase-field model, described in the previous section, should be inserted into Algorithm 1. Therefore, a finite-element discretization, defined on a quadrilateral discretization T h of Ω with elements E of size h and piecewise bilinear functions

Q 1 := span{1, x 1 , x 2 , x 1 x 2 } is introduced X h := {v h ∈ C(Ω) ∩ X : v h | E ∈ Q 1 for all E ∈ T h } , Y in,h := {ψ h ∈ C(Ω) ∩ Y in : ψ h | E ∈ Q 1 for all E ∈ T h } . (29)
The discrete optimization problem reads as follows:

Problem 4. Given φ 0 h and for the incremental steps t n , with n = 1, ..., N , find

u n h ∈ X and φ n ∈ Y in , such that (u n h , φ n h ) = arg min (v h ,ψ h )∈X h ×Y in,h B ℓ ((v h , ψ h )), ( 30
)
where B ℓ is defined in Equation (28).

Numerical solution and software

The previous system requires linear and nonlinear solvers. There exist three types of nonlinearities, (i) the variational inequality constraint due to the crack irreversibility condition, (ii) the decomposition of the strain energy, and (iii) the nonlinear coupling between the displacement equation and the phasefield variational inequality. In Ref. [START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF], a combined Newton algorithm was suggested in which the crack irreversibility is addressed in terms of a primal-dual active-set method and the nonlinear system with a classical line-search assisted Newton method. Both steps are combined into one loop and implemented in open-source pfm-cracks [START_REF] Heister | pfm-cracks: A parallel-adaptive framework for phase-field fracture propagation[END_REF] based on the finite-element library deal.II [START_REF] Arndt | The deal.ii finite element library: Design, features, and insights[END_REF][START_REF] Arndt | The deal.ii library, version 9.4[END_REF]. Recently, improvements of the combined active set and Newton solver have been suggested including supporting numerical tests [START_REF] Kolditz | A modified combined active-set newton method for solving phase-field fracture into the monolithic limit[END_REF].

Within the nonlinear combined Newton iterations, the linear systems are solved with GMRES (generalized minimal residual method) [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF], which is preconditioned with a block-diagonal system, where the blocks are approximated with algebraic multigrid from Trilinos [START_REF] Heroux | An overview of the trilinos project[END_REF]. Parallelization is based on MPI, where the adaptive meshes are organized by p4est [START_REF] Burstedde | P4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF].

Application

The previous developments are taken to design an integrated digital image correlation scheme based on the phase-field fracture model. Algorithm 1 has been implemented in the Correli 3.2 framework, which is an extension of the 3.0 version [50], namely, a Matlab based library with C++ kernels for digital image correlation. The phase-field codes to determine fracture propagation are described and referenced in Section 3.3.

Small-sample experiment

The first synthetic experiment to analyze the proposed IDIC approach is a slightly adapted version of the single edge notched shear test presented in Ref. [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF]. The sample is of square shape and each edge is 1 mm long, see Figure 1(a). To generate synthetic images, the phase-field model is used to compute displacement fields for time steps t k = k • 10 -4 second with k = 0, ..., 75. The displacement fields are then applied to the reference image shown in Figure 1(b),

I t (x, t k ) = I 0 (x -u(x, t k )). (31)
The definition (i.e., number of pixels) of the deformed images is identical to that of the reference picture.

Therefore Equation ( 31) is solved for all pixels of each synthetically deformed image. The field of view of the images is slightly larger than the sample to have the fully deformed domain inside of the images.

In the reference image, the sample is identified by the speckle pattern and the added field of view by the black background. The speckle pattern was generated using a Halton sequence placing ellipses of normally distributed intensity with radii uniformly ranging from 5 to 25 px on the sample surface. The pixel size is about 1 µm. The Dirichlet boundary conditions, chosen to be feasible using multi-axial testing machines [START_REF] Carpiuc-Prisacari | Comparison between experimental and numerical results of mixed-mode crack propagation in concrete: Influence of boundary conditions choice[END_REF], are u| Γ0 = 0 [mm] and u| Γ1 = (-2 τ, 1.5

x 1 τ ) T [mm], ( 32
)
where τ is the dimensionless time. Due to the end time of 7.5 • 10 -3 second and a pixel size of about 1 µm, the maximum displacement at the boundary Γ 1 is going to be about 15 px The material and numerical parameters of the synthetic test are gathered in Table 1. The element size h of the structured mesh is of 135 the order of 12 µm (or 11 px) and the internal length ℓ is assumed to be equal to 2 times the element size. The resulting displacement field together with the phase-field solution at frame 75 are shown in Figure 2.

The crack starts propagating from the beginning of the slit at frame 43 toward the bottom/left corner. In practice, one typically performs a first measurement using standard DIC techniques to estimate the boundary conditions, and take them as known for IDIC analyses [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted FEMU and integrated-DIC[END_REF], such that only material parameters need to be optimized. However, numerical models are very sensitive to boundary conditions [START_REF] Repin | Accuracy of mathematical models : dimension reduction, homogenization, and simplification[END_REF] with specific examples related to the current work in [START_REF] Carpiuc-Prisacari | Comparison between experimental and numerical results of mixed-mode crack propagation in concrete: Influence of boundary conditions choice[END_REF][START_REF] Heister | Parallel solution, adaptivity, computational convergence, and open-source code of 2d and 3d pressurized phase-field fracture problems[END_REF], such that slightly incorrect boundary conditions (due to measurement uncertainties in DIC analyses) may influence the quality of IDIC results. Therefore, the boundary conditions as well as the material parameters are optimized simultaneously, such that the space of parameterization is chosen to be P = R 5 with two parameters controlling the boundary condition on Γ 1 and the other three parameters are the Poisson's ratio ν, the fracture energy G c and the internal length ℓ. More precisely, a parameterization p = (a 1 , ..., a 5 ) T ∈ P is defined as

u(x, t)| Γ1×I = (a 1 t, a 2 x 1 t) T [mm], ν = a 3 , G c = a 4 and ℓ = a 5 h, (33)
where x = (x 1 , x 2 ) ∈ Γ 1 and t ∈ I. Using this parameterization, one should be aware that the Poisson's ratio should belong to the interval (-1, 0.5), and the fracture energy should be greater than zero. Furthermore, the internal length ℓ has to be at least as large as the size of one element, e.g. h such that ℓ > h, which is a well-known numerical requirement in the case of low-order finite elements [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][Section 2.1.2] or [82][Section 5.5]. With respect to the selected resolution, the internal length has to be larger than 11 px. Consequently, one has to add inequality constraints -1 < a 3 < 0.5, a 4 > 0 and a 5 ≥ h,

into the minimization problem. In the presented study, a simple active-set scheme was added inside Algorithm 1. The initial parameterization

p 0 = (a (0) 1 , ..., a (0)
5 ) = (-1, 0.5, 0.4, 2, 3) T (35

)
was selected because the corresponding solution did not include a crack, see Figure 3, and is therefore suitable to test the method for non-optimal initialization. This initial guess is also used to define the basis of P

Φ P i = a (0) i e i , (36)
where e i is one in its i-th component and zero else. This definition is especially useful for the finite-difference scheme. By choosing a perturbation parameter s = 0.01, one automatically perturbs the parameters by 1% wrt. its initial values. As one can assume that the initial parameterization is at least of the same order as the real solution, the finite-difference approximation is a reasonably good approximation of the actual derivative of the model.

Sensitivity and conditioning study

The sensitivity fields

S i = M ′ s (p)(Φ P i )
, where M ′ s is defined in Equation ( 9), visually represent at which rate the displacement field changes for a given perturbation of the parameters p at a certain time instant and location. matrices

H inst ij (t) = 1 |Ω| Ω (S i • ∇I M(p) )(S j • ∇I M(p) ) dx, (37)
and the full space-time Hessian matrix

H full ij = 1 |Ω × I| I Ω (S i • ∇I M(p) )(S j • ∇I M(p) ) dx dt. ( 38
)
For both cases, the Hessian matrix of the exact solution is analyzed such that not only the first but also the second term of the true IDIC Hessian (8) were omitted. In the following, the condition numbers for the Hessian matrices are analyzed in more detail. While the condition number indicates good or bad convergence 165 rates [START_REF] Nocedal | Numerical optimization[END_REF], one can also interpret it as an amplification factor of errors existing due to, e.g., acquisition noise [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF][START_REF] Mathieu | Estimation of elastoplastic parameters via weighted FEMU and integrated-DIC[END_REF]. The condition number of the instantaneous Hessian matrix decreases significantly with time (Figure 5) and reaches its lowest level after the damaged zone already propagated and at the end of the test.

Interestingly, the condition number also increased during damage propagation. The full space-time Hessian has a lower condition number than any instantaneous Hessian.

R i = |a i -a real i | |a real i | , (39)
one observes that the first reduction of the residual is related to the boundary conditions. Moreover, the parameters a 1 and a 2 are, on average, the most accurate parameters. Even though the residual barely reduced after iteration 7, the parameters and the displacement field still changed quite a lot. After the fifth iteration, an exponential convergence rate was reached and after the twelfth iteration, the sum of ∥δa i /a (0)

i ∥ was less than 10 -3 such that the iteration was stopped.

Influence of acquisition noise

So far, the virtual images were assumed to be almost optimal i.e., they were constructed such that the gray level conservation held, except for interpolation errors. Real images include much larger errors due to acquisition noise that can be modeled by Gaussian white noise η such that the gray level conservation becomes

I 0 (x) = I M(p) (x, t) + η(x, t). (40)
The influence of noise is analyzed by adding Gaussian white noise with different standard deviation levels to all images, while the identification scheme remains the same. In Figure 9, the relative error R i for each parameter is reported for different levels of noise. Despite the fact that one would typically only expect standard deviations of about 5% of the dynamic range (at the most) i.e., 12.75 GL (Gray Levels), 205 the performance of the identification procedure with standard deviations up to 20% of the dynamic range (51 GL) was investigated. Having also converged for the 20% level indicates that the presented algorithm is robust and remains very accurate. Noticeably, the relative error is less than 1% of the actual solution for 5% noise. The following boundary condition

u| Γ0 = 0 [mm] and u| Γ1 = (-50 τ, 1/4 x 1 τ ) T [mm] (41)
results in a similar crack pattern as in the first test, see Figure 11, but this time with the end time of The parameterization p = (a 1 , ..., a 5 ) T ∈ P is defined as

u| Γ1×I = (10a 1 t, a 2 /20 x 1 t) T , ν = a 3 , G c = a 4 and ℓ = a 5 h, (42)
and the initial parameterization was chosen to be

p 0 = (a (0) 1 , ..., a (0)
5 ) = (-1, 1, 0.4, 0.5, 3) T , (43)
such that, similar to the first test, no crack has propagated.

Influence of noise 225

Gaussian white noise was added up to 20% of the dynamic range. One can see in Figure 13 that the error increased by a factor of about 100 between 1% (0.5 GL) to 3% (1.5 GL) noise. Nevertheless, today's cameras are able to acquire images with noise levels of about 1% of the dynamic range, such that also the internal length is expected to be measurable in real-world settings. Furthermore, the error remained relatively stable even for high noise levels. Thus, the low eigenvalues lead to robustness with respect to Gaussian noise.

Conclusion and Outlook

In this study, a phase-field fracture model was included into an integrated-DIC framework, enabling for the identification of material parameters and boundary conditions. In the first parts of this work, algorithms and conceptional developments were designed in order to utilize numerically a modern phase-field fracture approach with primal-dual active set techniques and parallel computing. Phase-field fracture acted as inner approach of IDIC. Two synthetic experiments were designed based on realistic sample geometries, boundary conditions, image contrast and acquisition noise. The samples were deformed with the phase-field model with realistic material parameters until a crack initiated and propagated.

The results showed the capability to simultaneously identify the boundary conditions and material parameters for the whole acquisition of a single experiment. For both experiments, the parameter sensitivity fields, their relative coupling and amplitudes were investigated at different crack propagation steps. Remarkably, high sensitivities were obtained for Poisson's ratio and the internal length through the proposed space-time formulation for the first test case (i.e., the influence of acquisition noise on the calibration results remained limited over a very large range of standard deviations). The identification of those two parameters is frequently highlighted as challenging in the literature. For the second test case, which was more challenging in terms of noise to signal ratio, good calibration results were obtained for a more limited range of acquisition noise. Such observations show that the sensitivity analyses as discussed herein are very useful to probe IDIC frameworks prior to actual experiments.

Since the present analysis effectively demonstrated the potential for coupling the phase-field fracture model with IDIC using synthetic test cases, a future study is to apply it to real experiments (e.g., the socalled CARPIUC benchmark [START_REF] Carpiuc | Carpiuc benchmark overview: crack advance, reorientation, propagation and initiation under complex loadings[END_REF]), from which the synthetic experiments were inspired. Using projectionbased digital volume correlation (P-DVC), it would also be possible to extend the analysis to 4D measurements [START_REF] Jailin | Fast four-dimensional tensile test monitored via x-ray computed tomography: Elastoplastic identification from radiographs[END_REF][START_REF] Kosin | A projection-based approach to extend digital volume correlation for 4D spacetime measurements[END_REF]. Additionally, even though the presented method converged fast, it may be interesting to further investigate the computation of the sensitivities for each parameter. This step requires a full computation of the phase-field model for each parameter at each Newton iteration, which is in the end a costly operation.

Figure 1 :

 1 Figure 1: Small-sample experiment. Geometry of the sample and applied boundary conditions (a), speckle pattern applied to the sample surface (b). The pixel size is about 1 µm.

Figure 2 :

 2 Figure 2: Solution of the phase-field model of the small-sample test at frame 75 for the reference parameters. Phase-field variable φ (a), displacements u 1 (b) and u 2 (c).

Figure 3 :

 3 Figure 3: Solution of the phase-field model of the small-sample test at frame 75 for the initial parameterization p (0) . Phase-field variable φ (a), displacements u 1 (b) and u 2 (c).

Figure 4 :

 4 Figure 4: Variation of the phase-field variable for 1% perturbation of the parameters (left column), the first component of the displacement field (middle column) and the second component of the displacement field (right column) for the small-sample experiment. All fields are scaled by the perturbation factor s = 0.01.

Figure 5 :

 5 Figure 5: Condition number of the space-time Hessian matrix H full compared to that of the instantaneous Hessian matrix H inst for each image of the small-sample experiment.

Figure 6 :

 6 Figure 6: Solution of the phase-field variable for frames 20 (a), 50 (d) and 75 (g) for the small-sample experiment. Rows one to three correspond to instantaneous Hessian matrices and row four with the full space-time Hessian matrix. The eigendecomposition of the Hessian matrices gives the eigenparameters shown in (b,e,h,j) and the eigenvalues shown in (c,f,i,k).

4. 1 . 2 .

 12 ResultsInstead of using the full sample Ω = [0, 10] 2 for the IDIC problem, a slightly smaller region excluding the slit was used. At the slit, high interpolation errors were observed, which influenced the accuracy of the converged solution. The mask shown in Figure7corresponds to the black area. The residuals of the initial and converged parameterizations are displayed in the region of interest at frames 20, 43 and 75, corresponding to the (essentially) elastic phase, the onset of damage propagation and the end of the test, respectively. The residuals do not vanish completely, but their RMS level is reduced from 38.7 to 1.5, and they have a periodic pattern, which is traced back to gray level (spline) interpolation errors when creating the image series I M(p) . Interestingly, the maximum value of the residuals decreases in time, as indicated by the residuals in Figure7(d-f). This is the result of the increasing sensitivity of the parameters, which was observed in Section 4.1.1.

Figure 7 :

 7 Figure 7: Initial (a-c) and converged (d-f) gray level residuals at frames 20 (a,d), 43 (b,e) and 75 (c,f) for the small-sample experiment.

Figure 8 :

 8 Figure 8: Small-sample experiment. Convergence of the residual J 2 and difference of the approximate and true displacement compared to an exponential convergence rate (dashed line) (a). Convergence of the relative error R defined by Equation (39) for the five identified parameters (b).

Figure 9 :Figure 10 :

 910 Figure 9: Convergence results for the small-sample experiment in terms of parameter errors for different levels of standard deviation of Gaussian white noise.

Figure 11 :

 11 Figure 11: Solution of the phase-field model of the large-sample test at frame 75 for the reference parameters. Phase-field variable φ (a), displacements u 1 (b) and u 2 (c).

220 4 . 2 . 1 .

 421 Conditioning studySimilar to the first experiment, the full space-time Hessian matrix has a lower condition number than any instantaneous Hessian matrix as shown in Figure12. Yet, the eigenvalues are much lower, while the eigenmodes are essentially the same.

Figure 12 :

 12 Figure 12: Condition number of the space-time Hessian matrix H full compared to that of the instantaneous Hessian matrix H inst for each image of the large-sample experiment (a), the eigenparameters of the space-time Hessian (b) and its eigenvalues (c).

Figure 13 :

 13 Figure 13: Convergence results for the large-sample experiment in terms of parameter errors for different levels of standard deviation of Gaussian white noise.

Table 1 :

 1 Reference material and numerical parameters to generate the displacement fields for the deformed configurations for the small-sample experiment.

	Parameter	Definition	Value
	h	Diagonal cell diameter	12 µm (or ≈ 11 px)
	G c	Fracture energy	2.7 kJ/m 2
	E	Young's modulus	250 GPa
	ν	Poisson's ratio	0.2
	ℓ	Internal length	2 h ≈ 24 µm (or ≈ 22 px)
	κ	Regularization parameter	10 -10

Table 2 :

 2 Reference material and numerical parameters to generate the displacement fields for the deformed configurations for the large-sample experiment.

	Parameter	Definition	Value
	h	Diagonal cell diameter	2.4 mm (or ≈ 23 px)
	G c	Fracture energy	0.3 kJ/m 2
	E	Young's modulus	15 GPa
	ν	Poisson's ratio	0.2
	ℓ	Internal length	2 h ≈ 4.7 mm (or ≈ 46 px)
	κ	Regularization parameter	10 -10

Preprint submitted to Computer Methods in Applied Mechanics and EngineeringDecember 29, 2023

The second experiment is used to analyze Algorithm 1 in a more realistic setting. Therefore, the sample is enlarged to have the same size as in the CARPIUC benchmark[START_REF] Carpiuc | Carpiuc benchmark overview: crack advance, reorientation, propagation and initiation under complex loadings[END_REF], see Figure10(a), but the speckle pattern and the overall shape are kept the same as in the previous section, Figure10(b). The material parameters are also adapted to be closer to a quasi-brittle material, see Table2.

τ = 7.5 • 10 -3 and a pixel size of about 104 µm the maximum displacement at the boundary Γ 1 is about 3.6 px. Therefore, the maximum displacement in terms of pixels is reduced by a factor of about 5 compared to the small-sample test. This level allows for a comparison of the two tests and the convergence of the algorithm.

Acknowledgments

The financial support of the French-German University through the French-German Doctoral college "Sophisticated Numerical and Testing Approaches" (CDFA-DFDK 19-04) is acknowledged. Discussions are acknowledged within the framework of the International Research Training Group on Computational Mechanics Techniques in High Dimensions GRK 2657 funded by the German Research Foundation (DFG) under Grant Number 433082294.