Bayesian inversion for unified ductile phase-field fracture

The prediction of crack initiation and propagation in ductile failure processes are challenging tasks for the design and fabrication of metallic materials and structures on a large scale. Numerical aspects of ductile failure dictate a sub-optimal calibration of plasticity- and fracture-related parameters for a large number of material properties. These parameters enter the system of partial differential equations as a forward model. In this work, we develop a step-wise Bayesian inversion framework for ductile fracture to provide accurate knowledge regarding the effective mechanical parameters. To this end, synthetic and experimental observations are used to estimate the posterior density of the unknowns. To model the ductile failure behavior of solid materials, we rely on the phase-field approach to fracture, for which we present a unified formulation that allows recovering different models on a variational basis. In the variational framework, incremental minimization principles for a class of gradient-type dissipative materials are used to derive the governing equations. The overall formulation is revisited and extended to the case of anisotropic ductile fracture. Three different models are subsequently recovered by certain choices of parameters and constitutive functions, which are later assessed through Bayesian inversion techniques. To estimate the posterior density function of ductile material parameters, three common Markov chain Monte Carlo (MCMC) techniques are employed: (i) the Metropolis-Hastings algorithm, (ii) delayed-rejection adaptive Metropolis, and (iii) ensemble Kalman filter combined with MCMC. To examine the computational efficiency of the MCMC methods, we employ the R-convergence tool.


Introduction
Fracture in the form of evolving crack surfaces in ductile solid materials exhibits dominant plastic deformation. In comparison to brittle materials, the crack evolves at a slow rate and is accompanied by a huge plastic distortion. The prediction of such failure mechanisms due to crack initiation and growth coupled with elastic-plastic deformations is an intriguingly challenging task and plays an extremely important role in various engineering applications.
Recently, in the setting of continuum mechanics, a new perspective was proposed for embedding microscopic mechanisms into the macromechanical continuum formulation, based on a multi-field incremental variational framework for gradient-extended standard dissipative solids [1,2]. Typical examples are theories of gradient-enhanced damage [3,4,5,6], phase-field models [7,8,9], and strain gradient plasticity [10,11,12]. Such models incorporate non-local effects based on length scales, which reflect properties of the material micro-structure size with respect to the macro-structure size. In this context, the term size effects is used to describe the influence of the macro-structure size on the mechanical response during inelastic deformations. Thus, micro-structure interaction effects are introduced through the so-called local length-scale, which describes the gradient information of the quantity of interest within neighboring material points (e.g., the damage or ductility zones). From a mathematical point of view, local length-scales regularize both the plastic response as well as the crack discontinuities. Hence, it resolves the loss of ellipticity of the governing equations and avoids pathological mesh-dependence in postcritical ranges, as well documented in [13,14,15]. Within the variational framework for gradient-extended dissipative phenomena, the modeling challenge is two-fold.
• First, the derivation of well-posed theoretical formulations for describing the forward model. Hereby, variational phase-field modeling is considered, which is a regularized approach to fracture with a strong capability to simulate complicated failure processes. This includes crack initiation (also in the absence of a crack tip singularity) [16,17,18], propagation, coalescence, and branching, without additional ad-hoc criteria [8,19]. A summary of multiphysics phase-field fracture models is outlined in [20].
• The second challenge is to elucidate the backward model to estimate the model parameters and other univariate quantities of interest. A Bayesian estimation model (as an inverse model ) is here used for the ductile fracture problem to provide accurate knowledge regarding the effective mechanical parameters.
The majority of the ductile phase-field models found in the literature are based on local plasticity. In this setting, a strong localization of plastic strains may occur during the post-critical regime, while the damage gradient, as well as the displacement field, suffer jumps [41,42]. These occurrences are particularly relevant in the case of perfect plasticity due to the absence of a plastic regularization mechanism. Thus, from a numerical perspective, the use of local plasticity in phase-field models does not ensure mesh-objective simulations in the post-critical regime and may lead to non-realistic localized responses, such as ductile fracture with damage growth in non-plasticized regions [43]. To address these problems, phase-field models coupled to gradient-extended plasticity have been proposed in the literature, which incorporate a plastic internal length scale, in the spirit of [44]. The resulting formulation allows for a physically meaningful description of the coupled plasticity-damage evolution and mesh-objective finite element simulations. Models of this class were considered in [45,40], where a variationally consistent energetic formulation was adopted to derive the coupled system of partial differential equations (PDEs) that governs the gradient-extended elastic-plastic damage response. This approach is consistent with the models proposed in [46] in a finite-strain setting and the extensions to micromorphic regularization [14,47,43], where the governing equations were derived from rate-type variational principles, namely, the principle of virtual power.
In this study, we present a unified formulation for ductile phase-field fracture based on variational principles, rooted in incremental energy minimization, for gradient-extended dissipative solids [2,48]. The coupling of plasticity to the crack phase-field is achieved by a constitutive work density function, which is characterized by a degraded stored elastic energy and the accumulated dissipated energy due to plasticity and damage. Three different models are subsequently recovered by certain choices of parameters and constitutive functions. Specifically, two phase-field models coupled to local plasticity are derived, followed by a model that considers gradient extended plasticity. The overall formulation is revisited and extended to the case of anisotropic ductile fracture. Thereby, at a specific material point, the stress state relates to the given direction (resembling solids enhanced with stiff fibers), which entails a deformation-direction-dependent solid material. Hence, similar to [49], a stiffness parameter is introduced to enforce the crack phase-field evolution according to the preferred fiber orientation.

Bayesian inversion as a backward model
Providing reliable mechanical parameters is essential in computational mechanics to construct models with accurate predictive ability. Many of these a priori unknown parameters cannot be estimated directly through experimental procedures, and a significant effort is often needed to obtain reliable values. Furthermore, material parameters fluctuate randomly in space, giving rise to spatial uncertainty through the geometry. Consequently, the development of a sound statistical framework emerges as an interesting approach to reliably estimate mechanical properties.
Bayesian inversion is a probabilistic technique used to identify unknown parameters. Hereby, a forward model (e.g., a system of PDEs) obtains a set of given data according to the prior density (the initial/prior information of the parameters) and gives a response related to the given unknowns. The output of the inverse problem is the posterior density, which is related to a reference observation (e.g., from experimental or synthetic measurements). This distribution provides very useful information concerning the parameter range, its standard deviation, and expectation.
Markov Chain Monte Carlo (MCMC) methods are frequently employed to extract the posterior distribution of a parameter of interest [50]. We propose several candidates according to the prior density (e.g., uniform/Gaussian) and determine whether the proposed one is rejected or accepted. The Metropolis-Hastings algorithm is one of the most common MCMC techniques due to its efficiency and easiness. In [51], a Bayesian frame-work according to this algorithm was developed to identify the material parameters in brittle fracture. The method suffers from slow convergence, since most of the candidates are rejected. To enhance its efficiency, several techniques have been recently developed, such as delayed rejection adaptive Metropolis (DRAM) [52], differential evolution adaptive Metropolis (DREAM) [53], and ensemble Kalman filter with MCMC (EnKF-MCMC) [54]. The reader is referred to [51,55,56,57,58] for the application of Bayesian inversion in applied sciences.
A Bayesian approach (as a backward model) to estimate model parameters for brittle fracture in elastic solids has recently been proposed in [51]. Specifically, the Metropolis-Hastings algorithm is devised therein to approximate model parameters based on synthetic measurements which are obtained through a sufficiently refined discretization space as the replacement of experimental observations. Thereafter, a Bayesian inversion framework towards hydraulic phase-field fracture was designed for transversely isotropic and layered orthotropic poroelastic materials [59]. Specifically, the DRAM algorithm was extended for parameter identification.
In this study, we develop a Bayesian inversion framework to identify the effective mechanical parameters in ductile fracture. To this end, we first introduce a methodology to estimate the unknowns in different stages, i.e., elastic, plastic, and fracturing responses for isotropic and anisotropic materials. Three specific MCMC techniques are used to estimate the material parameters using synthetic and experimental measured data. Afterwards, a fair comparison is drawn between two improved models to determine the better convergence rate. As previously mentioned, having accurate information regarding the material parameters will enhance the accuracy of the PDE-based model. For instance, in the present case of ductile fracture, the goal is to predict the dissipative response in different stages, anticipating crack initiation and its propagation during time. We will thus employ the inferred parameters in the model equations and compare the response with the initial knowledge, highlighting the role of the probabilistic approach in improving the model's performance.

Physical interpretation of the ductile parameters
In ductile fracture, crack propagation is affected by several material properties. Figure  1 shows the range of the different parameters for a wide variety of materials. The effective parameters required in the models are introduced below.
• The bulk modulus K indicates how much the solid will compress as a result of an applied external pressure, and denotes the relation between a change in pressure and the resulting decrease/increase in fractional volume compression. See, e.g., [60,61,62,63].
• The shear modulus µ is a positive constant, smaller than K, which indicates the response of the solid to shear stress (the ratio of shear stress to shear strain). Large shearing stresses give rise to flow and permanent deformation or fracture. See, e.g., [63,64,61].
The elastic properties of the solid can be alternatively described in terms of the Young's modulus and the Poisson's ratio, or any other pair of Lamé's parameters. In our previous work [51], it was reported that due to the boundness of the Poisson's • The yield stress σ Y represents the maximum stress that can be applied without exceeding a specified value of permanent strain. For a solid, it denotes the stress related to the yield point (the starting point of plasticity) where the material starts to deform in the plastic regime [64].
• The critical value α crit stems from a physical assumption that fracture evolution is promoted once a threshold value for the accumulated plastic strain has been reached [67].
• The specific fracture energy ψ c characterizes the dissipated energy during a complete damage process in a homogeneous volume element [71]. This property can be interpreted as the amount of strain energy density (strain on a unit volume of material) that a given material can absorb before it fractures [72]. A study for different materials can be found in [72].
In Section 3.4, the role of these quantities in different stages of the deformation process will be clarified.
The paper is structured as follows. In Section 2, a unified version for ductile phase-field fracture models is presented in a variational setting, making use of incremental energy minimization. In Section 3, we first introduce the MCMC techniques and explain the three specific Bayesian estimation methods. Afterwards, a parameter identification setting for ductile fracture is presented based on MCMC. Thereafter a review of different possibilities (specific observations) for the implementation of the Bayesian inversion framework is introduced. In Section 4, we employ the presented techniques to precisely estimate the effective parameters in different stages of the deformation process. This information allows us to enhance the accuracy of the models, as evidenced by very good agreements between simulations and experiments, highlighting the noticeable efficiency of the MCMC techniques. Furthermore, a fair comparison between the proposed models is outlined. Finally, the conclusions are drawn in Section 5.
2. Phase-field modeling of ductile fracture in anisotropic elasticplastic materials In this section, we summarize the material models considered for the phase-field approach to ductile fracture. Three models found in the literature are revisited and extended to anisotropic fracture, considering the case of transversely isotropic materials. To this end, a unified formulation is first provided in Sections 2.1-2.3. The three examined models are then recovered in Sections 2.4.1-2.4.3. These models will be analyzed in subsequent sections using Bayesian inversion techniques, aiming for parameter identification in anisotropic elastic-plastic fracturing materials.

Basic continuum mechanics
Let B ⊂ R δ be an arbitrary solid domain, δ = {2, 3} with a smooth boundary ∂B ( Figure 2). We assume Dirichlet boundary conditions on ∂ D B and Neumann boundary conditions on ∂ N B := Γ N ∪C, where Γ N denotes the outer domain boundary and C ∈ R δ−1 is the crack boundary, as illustrated in Figure 2b.
The response of the fracturing solid at material points x ∈ B and time t ∈ T = [0, T ] is described by the displacement field u(x, t) and the crack phase-field d(x, t) as u : Intact and fully fractured states of the material are characterized by d(x, t) = 0 and d(x, t) = 1, respectively. In order to derive the variational formulation, the following space is first defined. For an arbitrary A ⊂ R δ , we set We also denote the vector valued space Concerning the crack phase-field, we set where d n is the damage value in a previous time instant. Note that W d dn is a non-empty, closed and convex subset of W d , and introduces the evolutionary character of the phasefield, incorporating an irreversibility condition in incremental form.
Focusing on the isochoric setting of von Mises plasticity theory, we define the plastic strain tensor ε p (x, t) and the hardening variable α(x, t) as and α : where R δ×δ dev := {e ∈ R δ×δ : e T = e, tr [e] = 0} is the set of symmetric second-order tensors with vanishing trace. The plastic strain tensor is considered as a local internal variable, while the hardening variable is a possibly non-local internal variable. In particular, α may be introduced to incorporate phenomenological hardening responses and/or non-local effects, for which the evolution equatioṅ is considered. As such, α can be viewed as the equivalent plastic strain, which starts to evolve from the initial condition α(x, 0) = 0 . Concerning function spaces, we assume sufficiently regularized plastic responses, i.e., endowed with hardening and/or non-local effects, for which we assume ε p ∈ Q := L 2 (B; R δ×δ dev ). Moreover, in view of (6), it follows that α is irreversible. Assuming in this section the setting of gradient-extended plasticity, we define the function spaces where W α = L 2 (B) for local plasticity, while W α = H 1 (B) for gradient plasticity. The hardening law (6) is thus enforced in incremental form by setting α ∈ W α αn, ε p −ε p n . The gradient of the displacement field defines the symmetric strain tensor of the geometrically linear theory as In view of the small strain hypothesis and the isochoric nature of the plastic strains, the strain tensor is additively decomposed into an elastic part ε e and a plastic part ε p as For simplicity, the anisotropic material is assumed to be strengthened by a single family of fibers, whose direction is described by a unit vector field a ( Figure 2). Consequently, the direction-dependent response is characterized by the second-order structural tensor The introduction of additional preferred directions can be easily incorporated in future work, following, e.g., [73]. The solid B is loaded by prescribed deformations and external tractions on the boundary, defined by time-dependent Dirichlet conditions and Neumann conditions u = u on ∂ D B and σ · n = τ on ∂ N B, where n is the outward unit normal vector on the surface ∂B. The stress tensor σ is the thermodynamic dual to ε andτ is the prescribed traction vector. Finally, the stress equilibrium is defined as the quasi-static form of the balance of linear momentum where dynamic effects are neglected and f is a given body force.

Energy quantities and variational principles
Let C denote the set of constitutive state variables. In the most general setting considered in this study, one has A pseudo-energy density per unit volume is then defined as W := W (C), which is additively decomposed into an elastic contribution W elas , a plastic contribution W plas , and a (regularized) fracture contribution W f rac : We note that W is a state function that contains both energetic and dissipative contributions. With this function at hand, a pseudo potential energy functional can be written as where E ext denotes the work of external loads: In variationally consistent models, the governing equations of the fracturing elastoplastic solid can be derived from knowledge of the energy functional (15) by invoking ratetype variational principles [1,47] in agreement with the principle of virtual power [74,75]. In such cases, a global rate potential of the form is defined, where Φ vis denotes the dissipative power density due to viscous resistance forces. In line with previous works [7], the function is considered, where η f and η p are material parameters that characterize the viscous response of the fracture and plasticity evolutions, respectively. Then, minimization of (17) with respect tou, the plasticity variables (ε p ,α) subject to the hardening law (6), and the crack phase-fieldḋ subject to the irreversibility conditionḋ ≥ 0 provide the governing equations for the elasticity problem, the plasticity problem, and the fracture problem, respectively. Such a variational structure results in a convenient numerical implementation based on incremental energy minimization, for which an algorithmic representation of the energy functional (15) is defined as where ∆t := t − t n denotes the time step. The coupled evolution problem then follows as the incremental minimization principle Remark 2.1. From equation (18), it is clear that the rate-independent case is recovered by letting η f → 0 and η p → 0. In this case, the coupled evolution problem can be equivalently derived in variational form using the energetic formulation for rate-independent systems [76,77], based on notions of energy balance and stability. This path is followed, for instance, in references [41,24,71,78,45,40]. Moreover, the fact that (14) is a state function implies that the incremental rate-independent problem exactly recovers the continuous counterpart.
For the variational formulation setting, it suffices to define the constitutive energy density functions W elas , W plas , and W f rac to establish the multi-field evolution problem in terms of (20). As we shall recall in the sequel, such a variational structure is not always present in phase-field models for ductile fracture, resulting in greater flexibility at the cost of a convenient mathematical structure.
2.2.1. Elastic contribution. The elastic energy density W elas in (14) is expressed in terms of the effective strain energy density ψ e . For transversely isotropic materials, ψ e is defined in terms of the elastic strain tensor ε e and the structural tensor M . In our formulation, in order to preclude fracture in compression, a decomposition of the effective strain energy density into damageable and undamageable parts is employed. Thus, we perform additive decomposition of the strain tensor into volume-changing (volumetric) and volume-preserving (deviatoric) counterparts: where the volumetric strain is denoted as ε e,vol (u) := 1 3 (ε e (u) : I)I and the deviatoric strain is denoted as ε e,dev (u) := P : ε e . The fourth-order projection tensor P := I − 1 3 I ⊗ I is introduced to map the full strain tensor onto its deviatoric component. Therein, I ijkl := 1 2 δ ik δ jl + δ il δ jk is the fourth-order symmetric identity tensor.
The anisotropic strain energy function. To complete the formulation, the anisotropic strain energy function reads where the stiffness parameter χ a characterizes the anisotropic deformation response with preferred direction a. The pseudo-invariant I 4 is defined as Extending the anisotropic energy into damageable and undamageable parts (see also [59]), admits the following splits where with the Macaulay bracket x ± := (x ± |x|)/2.
The total elastic strain energy function. The elastic contribution to the pseudoenergy density (14) finally reads where g e (d, α) is the elastic degradation function.
Following the Coleman-Noll procedure, the stress tensor is obtained from the potential W elas in (31) as where σ iso and σ aniso are the effective stress tensors, given by 2.2.2. Fracture contribution. The phase-field contribution W f rac is expressed in terms of the crack surface energy density γ l and the fracture length-scale parameter l f that governs the regularization. In particular, the sharp-crack surface topology C is regularized by a functional C l , as outlined in [14] and [48]. This geometrical perspective is in agreement with the framework of [79], which was conceived as a Γ-convergence regularization of the variational approach to Griffith fracture [80]. For the case of isotropic materials, the regularized functional reads In this work, following [49,59,81], anisotropic effects are introduced by means of the structural tensor M . In particular, we assume that γ l admits the additive decomposition In line with standard phase-field models, a general surface density function for the isotropic part γ iso l is defined as where ω(d) is a monotonic and continuous local fracture energy function such that ω(0) = 0 and ω(1) = 1. A variety of suitable choices for ω(d) are available in the literature [82,83,84]. Here, the widely adopted linear and quadratic formulations are considered, which yield, respectively, models with and without an elastic stage. Specifically, we define On the other hand, the anisotropic part γ aniso l reads γ aniso Finally, the fracture contribution to the pseudo-energy density (14) reads where g f is a parameter that allows to recover different models found in the literature, as will become apparent in the sequel.

Plastic contribution.
The plastic contribution W plas is expressed in terms of an effective plastic energy density ψ p , whose form will depend on the adopted phenomenological model. In line with previous works [43,45,40] let us consider a function in the context of gradient-extended von Mises plasticity: with the initial yield stress σ Y , the isotropic hardening modulus H ≥ 0 and the plastic length-scale l p . The plastic contribution to the pseudo-energy density (14) then reads where g p (d) is the plastic degradation function. The models presented in Sections 2.4.1 and 2.4.2 are restricted to local plasticity, for which l p = 0, while the model presented in Section 2.4.3 will include non-local effects, with l p > 0. For a variational treatment, it is convenient to invoke the energetic-dissipative decomposition of the plastic energy (41). Thus, the plastic energy density can be further decomposed as

Stationarity conditions and governing equations
Let us now derive the variationally consistent equations for the multi-field coupled problem. To this end, we seek to find the stationarity conditions for the minimization problem (20). The models presented in Sections 2.4.1-2.4.3 shall take the developments below as canonical forms, and will then deviate from the variationally consistent expressions in favor of greater flexibility.
2.3.1. Elasticity. The minimization with respect to the displacement field in the variational principle (20) yields which corresponds to the weak form of the mechanical balance equations (11), and W u 0 denotes the function space for the virtual displacement fields, i.e., with homogeneous kinematic boundary conditions.

2.3.2.
Fracture. The directional derivative of (15) with respect to the crack phasefield can be written as where the indicator function I + : R → R ∪ {+∞} has been introduced to impose the irreversibility condition embedded in d ∈ W d dn . Let us now define the fracture yield function The strong form of (44) can then be written as Recalling that the strong form yields, for the rate-dependent case, the evolution equation On the other hand, for the rate-independent case, we obtain the KKT conditions The main challenge in solving this evolution problem lies on imposing the irreversibility condition d ≥ d n , which allows to replace the set-valued expressions (44) or (46) by equalities. Several alternatives are available in the literature to tackle this problem, including simple penalization methods [85], augmented Lagrangian penalization [86], the primal-dual active set method [87], interior point methods [88], and the complementary system with Lagrange multipliers [89]. In this work, we employ the maximum crackdriving state function method based on the history field, as outlined in [48,90] and related works.

Plasticity.
The three models considered in this study employ von Mises plasticity in the rate-independent case, such that η p = 0 in (17). Moreover, as will become clear, the plasticity problem does not follow from the incremental minimization problem (20) in all three models. In particular, a variational formulation for the plasticity problem that is consistent with the governing equations is only possible if the elastic degradation function introduced in (31) does not depend on the hardening variable α, that is, g e := g e (d), such that W elas := W elas (ε, ε p , d; M ). Let us now summarize the variational formulation of such a model in the general gradient-extended case. A free energy density function for ductile phase-field fracture can be defined as Applying the Coleman-Noll procedure to the free energy density function (50) yields the following thermodynamic conjugate variables: In agreement with the classical setting of elasto-plasticity, the yield function is defined as With the yield function at hand, the strong form of the evolution problem follows from the principle of maximum plastic dissipation where Φ p is the plastic dissipation potential. The Euler equations of the maximization principle (53) follow as the flow rule and hardening laẇ together with the KKT conditions conditions Equations (54) and (55) constitute the so-called dual form of the elasto-plastic problem in strong form. To arrive at a primal formulation [91,92], the dissipation potential is evaluated from (53) as The Legendre transformation of Φ p then reads which yields, as a necessary condition, the primal representation of the plasticity evolution problem in the form of a Biot-type equation: From standard arguments of convex analysis [93,91], this expression implies the associative flow relations (54) as well as the loading/unloading conditions (55).
To derive the above governing equations from the incremental minimization problem (20), we make use of equations (42) and (56), such that the functional derivative of (19) with respect to {ε p , α} can be written as In view of equations (50) and (51), the strong form of (59) can be written as Equation (60) represents an incremental version of the compact evolution equation (58) for the plasticity model, recovered in a variationally consistent manner from the incremental minimization principle (20). Recall that in the present formulation, this was achieved by assuming an elastic degradation function that does not depend on the plastic variables. Model property

Specific models revisited
In this section, three benchmark phase-field models for ductile fracture, hereinafter labeled M 1 , M 2 , and M 3 , are revisited within the framework elaborated in the previous sections. The material parameters and constitutive functions that allow to recover each model from the general formulation are presented in Table 1.

2.4.1.
Local plasticity with G c based fracture criteria: Model 1 (M 1 ). The first model considered in this study takes the work from [23] as a point of departure. Therein, an extension of the model proposed in [22] was considered by further coupling the fracture process to plasticity through dependence of the elastic degradation function on the hardening variable α. The model was subsequently extended to finite strains and presented with experimental verification in [67]. As discussed in Section 2.3.3, dependence of the elastic degradation function on α results in lack of variational consistency for the plasticity evolution problem, in favor of greater flexibility. In this case, consider the Global Primary Fields : and the Constitutive State Variables : With the constitutive choices shown in Table 1 for M 1 , the following forms of the governing equations presented in Section 2.3 are obtained.
The strong form of the crack phase-field evolution (46) takes the form Herein, α crit is a threshold material parameter introduced in [23] to calibrate the softening response. To enforce the crack irreversibility condition, and, therefore, to cast this inequality constrained boundary value problem (BVP) as an equality constrained BVP, the history field is introduced. Here, ζ ≥ 0 is a scaling parameter that introduces further flexibility in the formulation, allowing to tune the post-critical range (cf. [14]). Equation (63) is then restated as With the last expression, and in view of (43), the global primary fields are found as the solution of the following coupled problem: find u ∈ W u u and d ∈ W d such that Thus, Remark 2.2. The introduction of the history field in the displacement Euler-Lagrange equation (63), finally yielding (M 1 ) 1 , results in a loss of variational consistency with respect to the energy functional (19) due to the filtering of the maximum history value of D and the scaling factor ζ for ζ = 1. The upside of this choice is a convenient numerical strategy for solving the original inequality-constrained PDE, and greater flexibility in the model.
The role of the parameter ζ, i.e., tuning the post-critical range by scaling the driving force, is already achieved in the present model by means of α crit . Consequently, ζ = 1 is assumed hereafter for M 1 .
Remark 2.4. At this point, it is worth noting that the crack driving force in (M 1 ) 2 , i.e., H, is scaled by the hardening variable α, such that the crack driving force vanishes for α → 0. As a consequence, fracture cannot occur outside the ductility zone, and a response corresponding to elastic damage followed by plastic damage is not possible in this model due to the strong coupling between damage and plasticity. For a detailed discussion of different possible elastic-plastic-damage evolution response, see [24].
Concerning the plasticity evolution problem, a variational derivation in the sense of (59) is not possible in the present model due to the dependence of the elastic degradation function g p on α. In this case, the local evolution of the plasticity variables {ε p , α} according to equations (54) and (55) (alternatively, (58) or (60) in the incremental form) is postulated in a non-variational context. For the present model, the yield function (52) takes the form 2.4.2. Local plasticity with ψ c based fracture criteria: Model 2 (M 2 ). The second model is based on the geometrically conceived approach to the phase-field modeling of ductile fracture, conceptually based on the local plasticity theory described in [48] and considered in subsequent works [35,34]. The original model is constructed within a variationally consistent framework, in agreement with the incremental energy minimization principle (20). In this case, consider the Global Primary Fields : and the Constitutive State Variables : With the constitutive choices shown in Table 1 for M 2 , the following forms of the governing equations described in Section 2.3 are obtained.
Letting l d := √ 2 l f (cf. [71]), and after simple manipulations, the strong form of the crack phase-field evolution (46) can be written as (cf. [48,14]): where the role of ψ c as a specific critical fracture energy density is clearly reflected. To enforce the crack irreversibility condition, and, therefore, to cast this inequality constrained BVP as an equality constrained BVP, the history field is introduced. Letting η d := η f /(2ψ c ), (69) is restated as With the last expression, and in view of (43), the global primary fields are found as the solution of the following coupled problem: find u ∈ W u u and d ∈ W d , such that Thus, Note that, in light of Remark 2.2, the introduction of the history field and the scaling parameter ζ in (70) results in a loss of variational consistency with respect to the energy functional (19) for the fracture problem.
As opposed to M 1 , the local plasticity evolution problem in the present model is variationally consistent (see Section 2.3.3). The local evolution of the plasticity variables {ε p , α} according to the evolution equation (60), which represents an incremental, primal version of equations (54) and (55), is then a necessary condition of the minimization principle (20). For the present model, the yield function (52) takes the form 2.4.3. Non-local plasticity with w 0 based fracture criteria: Model 3 (M 3 ). The third model considered in this study is inspired by the variational phase-field models coupled to gradient plasticity proposed in [48,14]. The modeling framework adopted therein and in subsequent studies [47,46,31] is consistent with the rate-type variational framework of [1]. In the small-strain rate-independent case, similar models were proposed in [27,45,40], where a variationally consistent energetic formulation was adopted to derive the governing equations. Consider, in this case, the Global Primary Fields : U := {u, d, α}, representing a combination of a first-order gradient plasticity model and a first-order gradient damage model. With the constitutive choices shown in Table 1 for M 3 , the following forms of the governing equations described in Section 2.3 are obtained.
With a slight change of parameters, the fracture problem in the present model admits the same formulation of M 2 in Section 2.4.2. In this case, according to Table 1, the strong form of the crack phase-field evolution (46) can be written as: where w 0 is a critical fracture energy density. One can show the identity of w 0 = 2ψ c holds for brittle fracture, but in the present gradient plasticity model, w 0 = 2ψ c due to the non-local term in ψ p . Indeed, the main difference of the present model with respect to M 1 and M 2 is that the plastic free energy ψ p , defined in (40), is considered here with l p > 0, and thus introduces non-local effects in the fracture driving force. To enforce the crack irreversibility condition, we define the history field Thus, (75) is restated as As before, in light of Remark 2.2, the introduction of the history field and the scaling parameter ζ in (70) results in a loss of variational consistency with respect to the energy functional (19) for the fracture problem.
As in M 2 , the plasticity evolution problem for the present model is variationally consistent (see Section 2.3.3). Moreover, the problem now includes non-local effects modulated by the plastic length-scale l p > 0, where the yield function (52) reads To derive the global PDE governing the evolution of the non-local field α, we take the weak form (59) as a point of departure, such that, for |ε p − ε p n | > 0: wheren := (ε p − ε p n )/|ε p − ε p n | is the direction of the plastic flow. Note that in (79), we have considered virtual fields such that the direction of the plastic flow is fixed, while the virtual equivalent plastic strain δα ∈ W α 0, δε p is allowed to vary. To solve (79), we must enforce the constraint embedded in α ∈ W α αn, ε p −ε p n (Equation (7)), such that Recalling that (79) is a weak representation of (60) for |ε p −ε p n | > 0, and that (60) implies the incremental version of the plastic flow rule (54), a possible way to proceed is to replace the local field ε p by setting, in agreement with (54), where, from standard arguments of von Mises plasticity, F p,trial := F p (ε, ε p n , d; M ). With this expression at hand, (79) may be left as a function of the non-local field α, subject to the irreversibility condition α ≥ α n .
Finally, the global primary fields are found as the solution of the following coupled problem: find u ∈ W u u , α ∈ W α , and d ∈ W d , such that Thus, Note that in (M 3 ) 2 , the indication function has been introduced to impose irreversibility of α, resulting in a multivalued expression. To show the consistency of this equation with the governing equations derived in Section 2.3.3, we first recall that the incremental flow rule has been enforced by means of (82). Then, we note that the strong form of (M 3 ) 2 yields β ∈ ∂ d I + (α − α n ), with β given in (78). In view of (47), it is easy to see that this expression represents the incremental version of the KKT conditions (55). To handle the inequality constraint and eliminate the multivalued term ∂ d I + (α − α n ), a constrained optimization technique is required. For instance, an interior-point method has been recently proposed in [88].
In the sequel, the phase-field models for ductile fracture formulated in M 1 , M 2 , and M 3 will be taken as inputs for the Bayesian inversion framework described in a detail in Section 4.

Parameter estimation based on Bayesian inference
In this section, we review different parameter estimation techniques based on MCMC to identify the mechanical parameters involved in ductile fracture. First, some basic statistical principals are briefly recalled.
In Bayesian estimation, a parametric forward model (e.g., a PDE-based model or a coupled variational inequality system) is used to update the available data (considered as random variables) based on the available information (denoting the prior knowledge). The posterior information is then provided as output [94,50,95].
Bayes' formula prescribes the probability of an event according to related prior information and is given by where P (A|B) denotes the conditional probability of event A happening when B has happened (likewise for P (A|B)), and P (·) is the probability of observations A and B.
Using a probability density function π, we can rewrite (83) as Here, π 0 (χ) is the prior distribution which indicates the available information regarding the parameter χ. For the ductile fracture case, the set of parameters χ is indicated in (107). Moreover, π(χ|m) denotes the posterior density, i.e., the probability density of the parameter χ considering the measurement m. The probability of the parameter χ with respect to the observation/measurement is described by the likelihood function π(m|χ). The denominator π(m) is a constant normalization factor, such that π(χ|m) ∝ π(m|χ)π 0 (χ).
In Bayesian inversion, the solution of the inverse problem is the posterior density giving the distribution of the unknown parameter values based on the sampled observations. MCMC is a popular method to calculate this distribution, where a Markov chain is constructed whose stationary distribution is the sought posterior distribution in Bayes' theorem.
In order to identify the unknown parameters, we introduce the following statistical model: where M is an n-dimensional vector that indicates the measurement, f denotes the PDEbased model, and χ = {χ 1 , χ 2 . . . , χ k } is a k-dimensional vector denoting the model parameters.
The model output f (x, χ) is the response quantity of interest, collected in an n-dimensional vector, where n = n T n C , with n C denoting the number of components of the response variable and n T denoting the number of time steps. In the present work, the force-displacement curve is taken as the response variable f (χ), such that n C = 1. For the measurement error ε, we employ a Gaussian independent and identically distributed error ε ∼ N (0, σ 2 I), where σ 2 is a fidelity parameter.
The inverse problem in the Bayesian framework can thus be stated as follows: given a measurement m, find the posterior density π(χ|m).
To this end, one makes use of Bayes' theorem of inverse problems [50], which can be stated as follows: Proposition 1 (Bayes' theorem for parameter estimation). We consider random parameter variables χ and a specific prior distribution π 0 (χ), and we consider m to be a realization of the random observation variable (denoting the measurement or virtual observation). The posterior distribution considering the measurement m follows as When using the above relation, one implicitly assumes that observed data is used to construct the posterior density. We should note that in our problem of interest, in case that an experimental measurement is not available, a virtual observation obs resulting from a fine spatial discretization is alternatively employed. Obviously, the observation is more valuable when a real experiment exists.
If we employ the statistical model (86) with the assumption that errors are Gaussian independent and identically distributed and ε i ∼ N (0, σ 2 ), where σ 2 is fixed, then the likelihood function is where is the sum of square errors.
Considering the given likelihood function (89), the posterior distribution has the following form: where SS ξ is the sum of squares defined by the integration variable; see (90). From a numerical point of view, we can approximate the integral as where the quadrature points and weights are denoted, respectively, by ξ i and w i .
In statistics, MCMC methods comprise a class of algorithms for sampling from a probability distribution. By constructing a Markov chain with the desired distribution as its equilibrium distribution, one can obtain a sample of the desired distribution by observing the chain after a number of steps. The more steps there are, the more closely the distribution of the sample matches the actual desired distribution.
In Bayesian statistics, the recent development of MCMC methods has been a key step in making it possible to compute large hierarchical models that require integration over hundreds or even thousands of unknown parameters. In rare event sampling, they are also used for generating samples that gradually populate the rare failure region.
Below, different popular MCMC methods are reviewed. These methods will be used to identify the parameters in ductile fracture in Section 4. A detailed comparison between the performance of the methods will be given to clarify their efficiency.

Metropolis and Metropolis-Hasting Algorithms
The Metropolis-Hastings (MH) algorithm is one of the most common techniques among the MCMC methods due to its simplicity for implementation and also its ability to handle different scientific/engineering problems (specifically when the parameters are not strongly correlated) [50]. In order to estimate the posterior distribution, in each iteration, a new candidate parameter value is proposed based on the current sample value according to a proposal distribution. Then, the acceptance ratio is calculated to decide whether the candidate value is accepted or rejected. The acceptance ratio points out how probable the new candidate value is with respect to the current sample.
The method was first introduced by Metropolis [96] based on a random walk. The algorithm starts from the initial guess (the prior value) χ 0 . Afterwards, according to the chosen proposal distribution a new candidate χ j−1 is proposed, which possibly depends on the previous candidates. Having the new candidate χ , the acceptance rate is calculated as As the next step, a random variable R ∼ Uniform (0, 1) is produced. If R < λ the candidate is accepted; otherwise, we reject the new proposal and keep the previous candidate in the chain. We follow this procedure for a sufficiently high number of replications. As seen, the algorithm is simple and efficient, specifically when a suitable proposal density is chosen and a large sampling is used. However, an inappropriate proposal results in a significant decrease in performance. If the proposal is very large, many of the candidates will be rejected; therefore, a good convergence to the target density (posterior distribution) will not be achieved. In contrast, if the proposal is too narrow, although many of the candidates are accepted, the chain movement is very slow, and many of the targets will not be captured.
In the Metropolis algorithm, a symmetric proposal density φ(χ |χ j−1 ) = φ(χ j−1 |χ ) is assumed. According to this condition, a movement towards the proposed candidate from the current point is equal to a backward movement (from the current candidate to the proposed point). The use of a non-symmetric proposal distribution was proposed as an efficient improvement by Hastings [97]. Considering N number of samples, the algorithm is summarized in Algorithm 1.
We can draw the following conclusions: • A proposal χ that results in π(m|χ ) > π(m|χ j−1 ) entails a small sum of squared error and thus leads to candidate acceptance.
• A proposal χ that leads to π(m|χ ) < π(m|χ j−1 ) entails a higher sum of squared error and the proposal may be rejected.
Regarding the proposal functions and how they affect the posterior distribution, if the variance is too large, a large percentage of the candidates will be rejected, since they will have smaller likelihoods, and hence the chain will stagnate for long periods. The acceptance ratio will be high if the variance is small, but the algorithm will be slow to explore the parameter space.
There are different measures to determine if the Markov chain is efficiently sampling from the posteriori density. A good criterion is the acceptance rate (the percentage of accepted candidates). The ratio can be used to tune the proposal density, i.e., reduce its variance. Another efficiency test is the autocorrelation function. The lag-τ autocorrelation function ACF : N → [−1, 1] is estimated by Here, χ j denotes the j-th element of the Markov chain andχ is the mean value. Note that ACF (τ ) is positive and monotonically decreasing. The interested readers can refer to [51], where the authors studied the effect of ACF on different parameters in phase-field modeling of brittle fracture. A more advanced convergence analysis such asR-statistics can be implemented when multiple MCMCs with different initial values are used. In Section 4, we will use such a diagnostic tool to compare the performance of the Bayesian techniques.

Delayed Rejection Adaptive Metropolis (DRAM)
At this point, it is worth discussing some improvements in the MH algorithm based on the proposal distribution. The main disadvantage of the model is that the covariance of the proposal should be tuned manually. To improve the efficiency, an alternative to using a fixed proposal distribution in each iteration is to update the distribution according to the available samples (adaptive Metropolis). This approach is useful since the posterior distribution is not sensitive to the proposal distribution.
To adapt the proposal function according to the obtained information, Haario [98] proposed a technique where the current point is chosen as the proposal center and the covariance function is updated using the estimated data. To this end, one can use the following proposal estimation where the parameter is chosen very small (close to zero) and S p = 2.38 2 j (as the scaling parameter). The covariance function is calculated by whereχ j = 1 j+1 j i=0 χ j [50]. The proposal adaptation can be done after a specific number of steps (e.g., 1000) instead of all steps. The efficiency of the algorithm can be further enhanced by adding a delayed rejection step. Mira [99] proposed that instead of a rejected candidate, a second stage is used to propose it from another proposal density. As the first step we propose the new candidate using the Cholesky decomposition of the covariance function (96): where U ∼ Uniform (0, I j ) and I j is the j-dimensional identity matrix. The alternative proposal χ * * is chosen using the proposal function where V j is the covariance matrix estimated by the adaptive algorithm [99]. The essential parameter is γ 2 , which will be chosen less than one so that the next stage has a narrower proposal function (normally, γ 2 = 1/5 is chosen). We use the following acceptance ratio: Then, similar to the MH algorithm, we follow the Markov chain to accept/reject the candidate. A summary of the process is given in Algorithm 2.
while j < N 1. Propose a new candidate χ * = χ j−1 + R j Z j where R j is the Cholesky decomposition of V j and Z j ∼ Uniform (0, I j ) where I j denotes the identity matrix.

MCMC with ensemble-Kalman filter
As previously mentioned, a good detection of the proposal density will enhance the Markov chain movement to the target density. Here, we introduce another proposal distribution detection technique using an ensemble-Kalman filter to obtain where ∆χ denotes the jump of Kalman-inspired proposal. In order to update the proposal, we can separate (99) into The first term indicates the so-called Kalman gain, i.e., where C χM is the covariance matrix between the inferred parameters and the PDE-based model, C M M is the covariance matrix of the PDE response, and C M denotes the measurement noise covariance matrix [100]. In (100), y j−1 is the residual of candidates with respect to the model. In other words, consideringm an observation/measurement, y j−1 =m − f (χ j−1 ) and s j−1 ∼ N (0, R), related to the density of measurement.
Considering the ductile fracture process, crack propagation is modeled until full fracture has occurred. To highlight this procedure, we have added criterion (iv) to Algorithm 3. For the two other algorithms (MH and DRAM techniques), the same condition can be considered.

Bayesian inversion for ductile phase-field fracture
A proper knowledge about the mechanical parameters that influence the behavior of fracturing solids is crucial to observe the model response and precisely predict crack initiation and propagation during different stages of the deformation process. Bayesian inversion techniques are convenient tools to monitor the crack behavior using observations (e.g., the measured data) and solving the inverse problem considering the forward model (here M 1 , M 2 , and M 3 ). Below, we review different possibilities for the implementation of Bayesian inversion in the context of elastic-plastic fracturing solids governed by phasefield models.
• Based on the load-displacement curve: this approach allows us to observe the crack behavior in all time steps up to complete failure. At time-step n, the load-displacement curve can be computed as where n is the outward unit normal on the surface, defined in (11). The main advantage of working with this curve is its easiness, since it involves a one-dimensional parameter. However, it is sensitive to the mesh size and the length scale; therefore, a sufficiently small (and thus more computationally expensive) mesh size is needed.

Algorithm 3 Bayesian inversion with ensemble-Kalman filter.
Initialization (j = 0): initiate the samples according to the prior density χ 0 While FLAG do (ii) update the Kalman gain 3. Set j = j + 1.
• Based on the point-wise primary fields: this approach monitors the crack behavior, the displacement, and the equivalent plastic strain in the entire geometry.
Here, the Bayesian setting strives to find the inferred parameters χ which minimize where L 2 -norm can be used, andū,d, andᾱ are the experimental data throughout geometry with respective measurement errors ε 1 , ε 2 , and ε 3 . This method is informative and provides precise information since the displacement and phase-filed in the entire geometry are considered. However, it is difficult and perhaps even impos-sible to obtain the measured data from actual experiments in a point-wise manner. Furthermore, a small mesh size must be chosen in numerical simulations to guarantee accurate estimations in the whole geometry, further rendering the method computationally prohibitive.
• Based on the point-wise phase-field propagation: this approach is less complex than the previous method. Reliable experimental values of the crack path can be obtained using X-ray or µ-CT scan in two-or three-dimensional problems; see, e.g., [101]. However, the computational issue regarding mesh sensitivity persists.
• Based on snapshots of proper orthogonal decomposition (POD): this approach employs a reduced order method (ROM) to reduce the computational complexity. Here, the snapshots of the solution (using measurements of the crack phasefield) are used to construct the POD basis [102]. If an efficient ROM is used, the computational complexity can be reduced significantly. Similarly, a Global-Local approach [103] can be employed to reduce the computational complexity of the forward model.
• Based on the effective stress-strain response: this approach explains the relation between ε and σ. It provides useful information regarding different material properties such as bulk modulus, hardening, and yield strength. Therefore, considering the availability of measurements, it entails an instructive procedure. Nevertheless, the computational costs, i.e., the effect of the mesh size and length scale, must be taken into account.
In this work, we choose the load-displacement curve as the observation, and sufficiently small mesh sizes are used to model the crack propagation. Nevertheless, it is worth noting that POD-ROM approaches have significant simulation advantages (noticeable computational cost reduction), while methods based on the stress-strain response are very informative. These procedures will be addressed in future works. Now, we proceed to establish a Bayesian inversion (BI) setting to identify the different parameters in ductile fracture. Let us assume that the response of ductile phase-field fracture is either elastic, followed by elastic-plastic, followed by elastic-damage (hereafter E-P-D); or elastic, followed by elastic-plastic, followed elastic-plastic-damage (hereafter E-P-DP). Next, we aim to determine the candidate χ ∈ (µ, K, H, σ Y , ψ c , G c , w 0 , l p ) as follows: (1) To findμ andK, we set H 0 → ∞, l 0 p → 0 (in case of M 3 ), χ 0 a → 0 (in the anisotropic case), and G 0 c → ∞ in M 1 , ψ 0 c → ∞ in M 2 , and w 0 0 → ∞ in M 3 , thus reflecting an elastic response. We then have (2) To findσ Y , we set H 0 → 0, l 0 p → 0 (in case of M 3 ), χ 0 a → 0 (in the anisotropic case), and G 0 c → ∞ in M 1 , ψ 0 c → ∞ in M 2 , and w 0 0 → ∞ in M 3 , thus reflecting an ideal plastic response. We then havẽ Step-wise Bayesian inversion method to determine the posterior density of the material unknowns for ductile phase-field fracture models.
(3) To findl p (in case of M 3 ) andH, we set χ 0 a → 0 (in the anisotropic case), G 0 c → ∞ in M 1 , ψ 0 c → ∞ in M 2 , and w 0 0 → ∞ in M 3 , thus reflecting an elastic-plastic response prior to fracture. We then have (4) To findχ a (in the anisotropic case),G c in M 1 ,ψ c in M 2 , andw 0 in M 3 , which reflect a ductile anisotropic fracture response, we have (5) Finally, we obtain the following parameter estimation: (µ, K, σ Y , H, l p , G c , ψ c , w 0 , χ a ) = BI μ,K,σ Y ,H,l p ,G c ,ψ c ,w 0 ,χ a . Figure 3 shows the overall procedure, indicating all stages of the deformation process. Note that, from the implementation point of view, we set the limit ∞ as 10 8 × E, where E refers to Young's modulus, while the lower limit is set to 0. From the statistical point of view, we employ the MCMC techniques (Algorithms 1-3) to identify the parameters in Step (1); then, we follow Step (2) to estimate σ Y and pursue the parameter identification procedure until we determine the whole set of material parameters in Step (4).
In a similar manner, if the response of ductile phase-field fracture is elastic, followed by elastic-damage, followed by elastic-plastic-damage, i.e., E-D-PD (see [24] for a detailed discussion on different possible evolutions), we first determine the elastic moduli, followed by the the anisotropic fracture properties by assuming that the response is brittle, and finally, we determine the plastic proprieties in the final dissipative stage.

Numerical examples
This section demonstrates the performance of the proposed Bayesian inversion approaches for parameter estimation within the ductile phase-field fracture models presented earlier. We investigate four numerical examples. To validate the numerical method, the last two examples are concerned with experimental observations, in which the posterior responses are compared with experimental load-displacement curves. The material parameters listed in Table 1 are considered, which are initialized based on [67,23]. In the MCMC the observational noise σ 2 = 10 −3 is used.
Space discretization. In the numerical simulations, the global primary variables are discretized using finite element basis functions, with bilinear quadrilateral Q 1 elements for the two-dimensional problems and trilinear hexahedral H 1 elements for the threedimensional problems.
Solution of the nonlinear problems. A staggered scheme is used for solving the variational equations resulting from the ductile phase-field fracture models (Section 2.4). For model 1 (M 1 ), and model 2 (M 2 ), we alternately solve for d/u by fixing u/d until convergence is reached. Accordingly, for model 3 (M 3 ), we alternately solve for u by fixing (α, d), and then solve for α by fixing (u, d). Next, we obtain the plastic strain tensor ε p though the incremental plastic evolution equation (81), and lastly, we find d by fixing (u, α), repeating the procedure until convergence is reached.
At this point, it is necessary to remark on the convergence criteria for the staggered scheme. Let n and k represent the loading time step and iteration counter of the alternate minimization scheme, respectively. At the fixed loading time step n, we obtain a converged state if the following holds: M • (u n,k , α n,k , d n,k ) ≤ Tol stag. with • ∈ {1, 2, 3} and Tol stag. ≈ 10 −3 . (108) Additionally, an iterative Newton solver is used in which the nonlinear equation systems are solved. The stopping criterion of the single scale and local Newton methods is Tol N-R = 10 −10 . Specifically, the relative residual norm is given by Residual : Here, F refers to the residual of the equilibrium equation of the nonlinear single scale and local BVPs. The interested reader can refer to [104,105,106,107,108,109,110] for the developed linear/nonlinear solvers for phasefield fracture.

Example 1: Asymmetrically I-shaped specimen under tensile loading
To gain a first insight into the performance of the Bayesian inversion approach, the following numerical example is concerned with the asymmetrically notched I-shaped specimen under tension. The configuration is shown in Figure 4a. The geometrical dimensions are set as H 1 = 110 mm, H 2 = 25 mm, r 1 = 3.625 mm, w 1 = 22 mm, and w 2 = 14.8 mm, with half-circular notches of radius r 2 = 2.5 mm. The two notches are placed at a vertical distance from the center of 10 mm.
A monotonic displacement increment ∆ū y = 2 × 10 −3 mm is applied in a vertical direction at the top boundary of the specimen. The minimum finite element size in the domain is 0.3 mm, for which the heuristic requirement h < l/2 inside the localization zone is fulfilled. Consequently, the I-shaped domain partition contains 21598 elements. The material and numerical parameters are given in Table 1 and Table 2, respectively.
For all three models, the shear modulus µ, the bulk modulus K, the hardening modulus H, and the yield stress σ Y are common. First, we study the effect of the common parameters on the load-displacement curve. Figure 5 shows the diagrams where M 1 is used to obtain the solutions. To monitor different critical values α crit and energy release   rates G c , as well as fracture energies ψ c , M 1 and M 2 are respectively employed. Moreover, we used M 3 to observe how the curve is affected by specific values of w 0 and the parameter ζ, as shown in Figure 5.
Next, we proceed to identify the effective parameters in the ductile fracture process using the Metropolis-Hastings algorithm introduced in Section 3.1. The Bayesian framework for ductile fracture is presented in Section 3.4. We employ a uniform distribution to estimate the parameters more accurately, as listed in Table 2 (the prior densities     initial values) and use N = 10 000 number of candidates. Regarding the reference values, a synthetic measurement is used, using a total number of degrees of freedom N dof = 28 380; the rest of the initial values are summarized in Table 2. The posterior density of the parameters using the three models is shown in Figure 6. The mean values of the posterior distributions are used to verify the parameter estimation. The inferred information is listed in Table 3. To verify the accuracy of the data, we employ the parameters in all three models and compute the load-displacement curve until the fracture point. Figure  7 shows the curves resulting from the different models and the reference observation. An excellent agreement indicates that the Bayesian inversion framework identified the parameters correctly, showing a consistent behavior for all three models in all stages.
The accuracy of the Bayesian inversion for all models enables us to provide equivalence for the model parameters. As previously mentioned, all models have four common parameters, but each model is also characterized by its own features. Here, we strive to find an equivalent value for different fracture energies. This allows us to use M 1 and M 3 and derive similar quantities in M 2 , and vice versa. To that end, we select the diagram estimated by M 2 as the reference observation, where all parameters are chosen according to the estimated values (see Table 3). However, ψ c varies between 25 and 65. We again use the MH algorithm to identify the equivalence of ψ c in M 1 (i.e, G c and α crit ) and M 3 (i.e., l p and w 0 ). The estimated quantities are summarized in Table 4. Figure 9 presents   Table 3 and Table 6 are used for Example 1 (left) and Example 2 (right).
the results obtained by the inferred values, where, again, Bayesian inversion provides a very good agreement.
The resulting equivalent plastic strain (α) and crack phase-field (d) at complete failure are shown in Figure 8. The solutions are based on the posterior density of the material parameters for different models, which are given in Table 3. Accordingly, the fracture path initiates within the maximum equivalent plastic regions, which appear near the notches. Next, the crack propagates in the plastic localization band, in which two cracks merge at the specimen center. It can be observed that even though the load-displacement curves shown are practically identical in all models, the corresponding phase-field profiles, and thus, hardening profiles, are not; see Figure 8. This can be explained, first of all, by the solution non-uniqueness of the phase-field fracture problem, and, secondly, by the fact that the different phase-field models in fact provide only the approximation of the fracture problem. Thus, the necessity of comparing the results with an experimental observation is crucial. Hereby, based on the experimental test provided in [67] (second experiment), a sharp crack transition between two notches is expected. Thus M 2 and M 3 seem to yield a more accurate fracture pattern.

Example 2: I-shaped tensile specimen for anisotropic ductile fracture
The main objective of this example is the adoption of Bayesian inversion for an anisotropic ductile phase-field fracture process. The BVP depicted in Figure 4b  Herein, we assume that the material constituents are not distributed uniformly through the continuum domain, and thus, the material is divided into several phases. Hence, heterogeneity in strength from one area of the domain to another one is expected. Note, however, that by means of the Bayesian inversion framework, we aim to determine the effective mechanical parameters. Here, we consider the parameters as a random field (with given mean and variation). Figure 10 illustrates the fluctuation of different material parameters (on the element-wise basis) with spatial correlation where a 10% variation is included. For instance, for the parameter K, the expectation is assumed as K=75 000 MPa, with a variation between 71 700 MPa and 78 800 MPa. This fluctuation will be used to provide the reference observation. Specifically, we will replicate 500 simulations (with a specific mesh size) to estimate the reference observation considering the mentioned variation. The distribution of the parameters on the geometry is shown in Figure 10.
The numerical example is performed by applying a monotonic displacement increment ∆ū y = 2 × 10 −3 mm in the vertical direction at the top boundary of the specimen ( Figure  4). To remove the rigid body motion, the bottom edge is fixed in the x−y directions. The minimum finite element size is 0.45 mm. The two-dimensional I-shaped domain partition contains 17038 elements.
In this example, to determine the effective mechanical parameters, the DRAM algorithm is used. Due to the anisotropic structure of the solid, in addition to the identified parameters in Example 1, the stiffness parameter χ a must be estimated. Here, we select φ = 45 o for the parameter identification and propose N = 10 000 candidates. A synthetic reference value using N dof = 26 870 is employed as the reference observation considering the already mentioned parameter variation. The prior (uniform) densities of the parameters are listed in Table 5. For each parameter, the inferred values by different models are relatively similar, which highlights the robustness of the Bayesian setting. The posterior distributions are depicted in Figure 11. The mean values of the posterior distributions are given in Table 6 using M 1 , M 2 , and M 3 . Again, we solve the model equations employing   Parameter the identified parameters to verify the effectiveness of the Bayesian framework, as shown in Figure 7 (right). All diagrams show that implementing the DRAM algorithm gives rise to a reasonable agreement between the models and also the reference observation.
For further investigation of the sensitivity of the inferred parameters obtained through the DRAM algorithm, two additional preferential fiber directions, namely φ = 30 o and φ = 60 o , are employed. Again, Figure 12 illustrates the robustness of the Bayesian setting, showing for all models and different orientations a consistent behavior in the in elastic, plastic, and fracture stages.  Table 8. An important observation is that both the equivalent plastic strain and the crack phase-field evolve in the direction of the preferred fiber orientation. Note that here, to avoid a strong deviation of the surrounding points due to the random normal distribution of the material properties, one could use an additional random distribution length-scale to achieve mesh objectivity, that is, the so-called heterogeneity length-scale; see [111]. In this study, we have not used a heterogeneity length-scale since we have assumed that the material distribution provides the synthetic observations.     Model Both DRAM and ensemble-Kalman filter (EKF) are efficient MCMC techniques and have shown their computational performance reasonably. However, a fair comparison can determine which method will be more advantageous in ductile fracture.

Convergence performance of the MCMC methods.
In the already mentioned examples, the MCMC techniques have been used to identify the mechanical parameters. The main advantage of the MH algorithm is its ease of implementation. However, it suffers from slow convergence, and the starting value may affect the convergence status. The DRAM and EKF variants, as more advanced techniques, show more a productive performance. In this part, we strive to study their efficiency in the context of ductile fracture.
Convergence diagnostics is essential in MCMC methods since it determines the accuracy of the parameter, and with how many iterations the chain converges to the target distribution. Here, we use multiple chains with different initialized values, expecting that a significantly large number of Markov chains gives rise to the same results. In other words, the candidate distribution from chains should be similar using multiple chains initial starting values. R-convergence diagnostics [112,113] is an efficient tool to monitor the convergence of the MCMC by comparing the between and within chain estimates for model parameters and other univariate quantities of interest. Assuming m parallel chains, we determine the variance between the chain means B/N and calculate the average of the within chain variances W . The target variance is given by where N is the length of the chain. Then, we calculate the potential scale reduction factor, or PSRF (also calledR-statistics) bŷ If the MCMC method converges appropriately, the chains are not affected by the starting point, andR reduces to 1. In other words, we can conclude that all chains are close to the target distribution [112]. In order to verify the method efficiency, a threshold can be defined, e.g., values less than 1.5 or 1.2 indicating a good convergence performance.
In the I-shaped example, in order to draw a comparison between DRAM and EKF algorithms, we study their convergence to conclude which model shows a faster convergence taking all inferred parameters into account. We use N = 1 000 and five parallel MCMCs (m = 5) with a uniform distribution indicated in Table 7. Using M 2 , Figure  15 shows that by employing EKF in all parameters, fewer candidates are necessary to converge to the posterior density. Indeed, for all parameters excluding K, after 1 500 samples,R-statistics converges to 1, showing a high level of accuracy. The performance of the DRAM is acceptable, since most of the variables after 2 500 samples are below the threshold, although again, the bulk modulus shows more variation (probably due to the large chosen prior density).
According to the above-mentioned discussion, we choose the EKF technique for the rest of the examples. We use five parallel MCMCs with 2 000 samples. For the observation, we use experimental data taken from [67]. Figure 16 shows the posterior density of estimated parameters, and the mean values are summarized in Table 8. We then employ   Table 8 and Table 10 are used for Example 3 (left) and Example 4 (right). The results are also compared with the simulation results of [67].
the identified quantities in all three models. Figure 18 shows a comparison between the load-displacement curve obtained by the models and the experiments. Interestingly, by employing the Bayesian framework, all crack propagation stages (until fracture reaches the boundary) are modeled accurately.
Next, we investigate the ductile failure response employing the posterior density of the material parameters given in Table 8. The evolution of the crack phase-field d is provided in Figure 18 at three deformation stages up to complete failure. Additionally, the equivalent plastic strain α at final failure is shown. It can be grasped that, regardless of the formulations, fracture initiation appears at the center of the specimen, and then evolves towards the two edges of the specimen until complete failure. The first important observation is that the simulation results are in well-agreement with the experimental failure pattern shown in Figure 14c. Another important observation is that the crack phase-field in M 2 and M 3 is more diffuse than M 1 . The main reason for this is that the crack driving force of M 1 is scaled by α; thus, the phase-field diffusivity is strongly coupled to the ductile response; see Remark 2.4.

Example 4: Sandia fracture challenge
The last example aims at estimating the posterior density of the unknown material properties for a specimen frequently used in the literature, namely, a Sandia fracture challenge [114]. The 2014 fracture challenge problem launched by the Sandia National Lab [114] has provided an ideal platform to assess the computational capability and limitations of each participating team [115]. Specifically, this challenge aims to evaluate the computational ability to predict crack initiation and propagation of ductile fracture with respect to the experimental observation. A recent comparative literature overview was conducted in [116]. As reported in [67], standard phase-field formulations without estimating accurate ductile material properties quantitatively overestimate the post-yielding load-displacement response, which can be improved by performing an accurate calibration of the plasticity and phase-field parameters. Hence, we aim at reproducing the experimentally observed Sandia fracture challenge through the proposed Bayesian inversion framework.
The experiments are based on the material Al-5052 H34, which experimentally induces a complex failure mode; see [65,67]. The configuration is shown in Figure 19a, while the experimental observations are shown in Figure 19c. The geometrical configuration includes two pins. The top pin is displaced vertically, while the lower pin is fixed in all directions. The two pins are considered to be rigid (here taken 10 times stiffer than the rest of the domain). The geometrical dimensions are set as H 1 = 80 mm, H 2 = 35 mm, H 3 = 15 mm, H 4 = 6.5 mm, w 1 = 22 mm, w 2 = 36 mm, and w 3 = 12.5 mm. The pins have an identical radius of r 1 = 6 mm, while the horizontal notch is rounded with a radius of r 5 = 3.25 mm. The specimen includes three voids with centers and radii c 2 = (x 2 , y 2 ) = (27, 32) mm and r 2 = 3 mm, c 3 = (x 3 , y 3 ) = (24, 45) mm and r 3 = 1.75 mm, and c 4 = (x 4 , y 4 ) = (22, 38) mm and r 4 = 1.75 mm, respectively. The specimen domain has a 2 mm thickness, as shown in Figure 19b.
The numerical example is performed by applying a monotonic displacement increment ∆ū y = 0.02 mm in the vertical direction at the top pin for 400 time steps. The minimum finite element size is 0.32 mm. Consequently, the Sandia specimen domain partition contains 17980 hexahedron linear elements.
Due to its efficiency, the EKF technique is again considered for this example to identify the parameters. We use four parallel MCMCs with 2 500 samples. The experimental data (reference observation) is taken from [67] and the prior densities are indicated in Table 9. The posterior distributions are shown in Figure 20, while the mean values are summarized  Parameter H µ K σ Y G c α crit ψ c w 0 l p ζ  Table 10. Finally, once again, we solve M 1 , M 2 , and M 3 using the inferred values, and compare the results with those obtained in [67]. As shown in Figure 17, by employing the proposed Bayesian inversion framework, the parameters are estimated accurately, showing a very good agreement between the simulated data in all models and the experiments. This implementation also enhanced the computational capabilities of M 1 , showing an improved model accuracy compared to the results obtained in [67]. For a better insight into the fracture process in all models, the evolution of the crack phase-field d is provided in Figure 21 at three deformation stages up to complete failure. Additionally, the equivalent plastic strain α at final failure is shown. Note that the solutions are based on the posterior density of the material parameters, which are given in Table 10.
It can be grasped that, for all three models, the fracture path first initiates at the void located in the middle of the specimen, and afterwards, evolves towards the right edge of the domain. In addition, a secondary crack initiates from the central void, but this time from its left side, and then propagates towards the left edge of the specimen until complete failure. Note further that the simulation results are in well-agreement with the experimental failure pattern shown in Figure 19c.
Based on our numerical result, it is worth noting that the stage of the secondary crack was no longer predicted by model 2. To estimate the posterior density of the material parameters through Bayesian inversion, several candidates are required. Thus, a stable forward method is crucial. Otherwise, deviating material properties will result in unstable solutions. In this context, we highlight that model 3 provides the most stable solutions, at the cost of an additional PDE that must be solved to obtain the plastic response (as opposed to local plasticity in model 1 and model 2).

Conclusion
In this work, we have proposed a robust and efficient step-wise Bayesian inversion method for ductile fracture problems using phase-field models. In particular, a Bayesian inversion framework (as a probabilistic technique) based on MCMC is developed to identify unknown ductile fracture parameters. Three common MCMC methods, namely the MH algorithm, DRAM algorithm, and EKF-MCMC have been used to estimate the effective parameters in ductile fracture. The posterior density results from the inverse problem are evaluated with synthetic measurements (for the first two examples) as well as experimental observations (for the last two examples). To approximate ductile failure, a phase-field fracture formulation is used for a ductile material exhibiting J 2 -plasticity in a quasi-static kinematically linear regime. To do so, we have presented a unified formulation for phase-field modeling of ductile fracture, which is resolved through an incremental energy minimization approach. The overall formulation is revisited and extended to the case of anisotropic ductile fracture. Three different specific models are subsequently recovered by certain choices of parameters and constitutive functions.
In the first numerical example, the equivalence of the parameters in the different models is provided. In the second example (anisotropic ductile fracture), we investigate the evolution of the failure response through the posterior density function obtained by the proposed step-wise Bayesian inversion method. We have shown that the equivalent plastic strain α as well as the crack phase-field d evolve in the direction of the preferred fiber orientation. The last two examples are concerned with the experimental observations to estimate the posterior density of the material unknowns. We observed that, although the MH algorithm can be implemented easily, it is sensitive to the initial guess, leading to slow convergence, and may depend on the prior density. As more advanced techniques, we compared the convergence of the DRAM and EKF-MCMC methods by employing a reliable convergence diagnostic tool, namelyR-convergence. Using the Kalman filter improves considerably the convergence of different MCMCs, i.e., fewer iterations are needed to obtain a high level of accuracy. We conclude that this method is more efficient compared to the DRAM algorithm.
Through our findings, the coupling scheme between the step-wise Bayesian inversion framework and ductile fracture simulations results in an accurate and a reliable information related to the model parameters. As a consequence, an excellent agreement was obtained between the results of all examined models and the experimental observations.