A Bayesian estimation method for variational phase-field fracture problems

by   Amirreza Khodadadian, et al.

In this work, we propose a parameter estimation framework for fracture propagation problems. The fracture problem is described by a phase-field method. Parameter estimation is realized with a Bayesian framework. Here, the focus is on uncertainties arising in the solid material parameters and the critical energy release rate. A reference value (obtained on a sufficiently small mesh) as the replacement of measurement will be chosen, and their posterior distribution is obtained. Due to time- and mesh dependency of the problem, the computational costs can be high. Using Bayesian inversion, we solve the problem on a relatively coarse mesh and fit the parameters. The obtained load-displacement curve that is usually the target function is matched with the reference values. Finally, our algorithmic approach is substantiated with several numerical examples.



There are no comments yet.


page 5

page 18

page 24

page 28


Bayesian Inversion for Anisotropic Hydraulic Phase-Field Fracture

In this work, a Bayesian inversion framework for hydraulic phase-field t...

Calibrate, Emulate, Sample

Many parameter estimation problems arising in applications are best cast...

Parameter Estimation in Computational Biology by Approximate Bayesian Computation coupled with Sensitivity Analysis

We address the problem of parameter estimation in models of systems biol...

A Bayesian Inference Framework for Procedural Material Parameter Estimation

Procedural material models have been graining traction in many applicati...

Gaussbock: Fast parallel-iterative cosmological parameter estimation with Bayesian nonparametrics

We present and apply Gaussbock, a new embarrassingly parallel iterative ...

A "parallel universe" scheme for crack nucleation in the phase field approach to fracture

Crack nucleation is crucial in many industrial applications. The phase f...

Bayesian Estimations for Diagonalizable Bilinear SPDEs

The main goal of this paper is to study the parameter estimation problem...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

This work is devoted to the quantification of uncertainties in fracture failure problems. To formulate fracture phenomena, the phase-field formulation for the quasi-brittle fracture solid is used. The variational phase-field formulation is a thermodynamically consistent framework to compute the fracture failure process. This formulation emanates from the regularized version of the sharp crack surface function FraMar98 . The regularized fracture phenomena are described with an additional auxiliary smooth indicator function, which is denoted as crack phase-field. Along, with a mechanical field (denoted by ), a minimization problem for the multi-field problem can be formulated. The main feature of variational formulation arises from the crack phase-field is to approximate the discontinuities in the across the lower-dimensional crack topology with the BourFraMar00 .

The resulting, regularized formulation leads to a diffusive transition zone between two phases in the solid, which corresponds to the fractured phase (i.e., ) and intact phase (i.e., ), respectively. The transition zone which is affected by a regularized parameter known as length-scale (denoted by ). Moroever, the regularized parameter is related to the element size . Specifically, hold such that through discretization error estimates AmTo92 ; Braides1998 . Therefore, sufficiently small regularized length-scale is computationally demanding. To date, the focus in such cases was on local mesh adaptivity and parallel computing in order to reduce the computational cost significantly; see for instance noii2019phase ; MaWi19 ; HeWheWi15 ; HeiWi18_pamm ; BuOrSue13 ; BuOrSue10 ; ArFoMiPe15 ; Wi16_dwr_pff ; MaWaWiWo19 . Another recent approach is a global-local technique in which parts of the domain are solved with a simplified approach GeNoiiAllLo18 ; NoAlWiWr19 that also aims to reduce the computational cost.

The main goal in this work is to quantify uncertain parameters of the phase-field formulation. The underlying framework of parameter estimation using Bayesian inference is described in the following. Bayesian inference is a probabilistic method used to estimate the unknown parameters according to the prior knowledge. The observations (experimental or synthetic measurements) can be used to update the prior data and provide the posterior estimation. The distribution provides useful information about the possible range of parameters and their variations and mean. Markov chain Monte Carlo (MCMC)

hoang2013complexity is a common computational approach for extracting information of the inverse problem (posterior distribution). Metropolis-Hastings (MH) algorithm hastings1970monte is the most popular MCMC method to generate a Markov chain employing a proposal distribution for new steps. In practice, a reliable estimation of influential parameters is not possible or needs significant efforts. In khodadadian2019bayesian ; mirsian2019new , the authors used the Metropolis-Hastings algorithm to estimate the unknown parameters in field-effect sensors. It enables authors to estimate PT-density of the target molecules which can not be experimentally estimated. The interested reader refers to minson2013bayesian ; cardiff2009bayesian for more application of Bayesian estimation in applied science. In the same line, other optimization approach can be used to determine intrinsic material properties of the specimen from experimental load-displacement curve, see e.g., noii2019characterization .

The current work follows two goals. First, using Bayesian inversion to determine the unknown and effective phase-field fracture parameters, i.e., Lamé 1st parameter and shear modulus as well as Griffith’s critical elastic energy release rate. To this end, a reference value (obtained on a sufficiently refined mesh) as the replacement of measurement will be chosen, and their posterior distribution is obtained. The second goal is related to the temporal and spatial mesh dependencies of the problem. The computational costs can be high, specifically when an appropriate estimation is required. Specifically, when we aim to have a variational phase-field modeling within multi-physics framework, see e.g. miehe2016phase ; heider2017phase ; heider2018modeling . Using Bayesian inversion; we strive to solve the problem with a coarser mesh and fit the parameters. The obtained load-displacement curve (as an important characteristic) is matched with the reference value. In spite of using coarser meshes and therefore significantly lower computational costs (in terms of CPU timings), the accuracy of the solution is reliable (crack initiation and material fracture time estimated precisely).

The paper is organized as follows: In Section 2, we describe the variational isotropic phase-field formulation for the brittle fracture that is a thermodynamically consistent framework to compute the fracture failure process. In section 3, the Bayesian inference is explained. We describe how the MH algorithm will be used to estimate the unknown parameters in phase-field fracture. Also we point out the critical points in the load-displacement curve, which must be estimated precisely with the Bayesian approach. In Section 4, the Bayesian framework is adopted to estimate unknown parameters in the phase-field fracture approach. In Section 5

, three specific numerical examples with different parameters and geometry will be given. We will use two proposal distributions (uniform and normal distribution) to sample the candidates and estimate the unknown parameters with different mesh sizes. Finally, in

6 we will draw paper conclusions and explain our future planes for employing Bayesian inversion in heterogeneous materials.

2 Variational isotropic phase-field brittle fracture

2.1 The primary fields for the variational phase-field formulation

We consider a smooth, open and bounded domain , . In this computational domain, a lower dimensional fracture can be indicated by . In the following, Dirichlet boundaries conditions indicated as , and Neumann boundaries conditions are given on , where is the outer boundary of and is the crack boundary. The geometric setup including notations is illustrated in Figure 1a. The surface fracture is estimated in . A region without any fracture (i.e., an intact region) is indicated by such that and .

Figure 1: (a) Geometric setup: the intact region indicated by and is the crack phase-field surface. The entire domain is denoted by . The crack phase-field is approximated in the domain . The fracture boundary is and the outer boundary of the domain is . is represented by means of such that the transition area is with thickness . (b) Regularized crack phase-field profile for a different length scale. A smaller value for the length scale lets the crack phase-field profile converge to a delta distribution.

The variational phase-field formulation is a thermodynamically consistent framework to compute the fracture process. Due to the presence of the crack surface, we formulate the fracture problem as a two-field problem including the displacement field and the crack phase-field . The crack phase-field function interpolates between , which indicates undamaged material, and , which indicates a fully broken material phase.

For stating the variational formulations, the spaces


are used. Herein, denotes a closed, non-empty and convex set which is a subset of the linear function space (see e.g., KiStam00 ).

2.2 Variational formulation for the mechanical contribution

In the following, a variational setting for quasi-brittle fracture in bulk materials with small deformations is formulated. To formulate the bulk free energy stored in the material, we define the first and second invariants as


with the second-order infinitesimal small strain tensor defined as


The isotropic scalar valued free-energy function reads


with the elastic Lamé constants denoted by and . A stress-free condition for the bulk energy-density function requires . Hence, the bulk free-energy functional including the stored internal energy and the imposed external energy is



is the imposed traction traction vector on


Following FraMar98 , we define the total energetic functional which includes the stored bulk-energy functional and fracture dissipation as


where is the so called the Griffith’s critical elastic-energy release rate Also, refers to the -dimensional Hausdorff measure (see e.g. BourFraMar00 ). Following BourFraMar00 , is regularized (i.e. approximated) by the crack phase-field (see e.g. BourFraMar00 ). Doing so, a second-order variational phase-field formulation is employed; see Section2.3. Additional to that, a second-order stress degradation state function (intacted-fractured transition formulation) is used as a monotonically decreasing function which is lower semi-continuous order; see Section 2.5.

2.3 Crack phase-field formulation in a regularized setting

Let us represent a regularized (i.e., approximated) crack surface for the sharp-crack topology (which is a Kronecker delta function) thorough the exponential function , which satisfies at as a Dirichlet boundary condition and as . This is explicitly shown in Figure 1b for different length scales. A regularized crack surface energy functional is


based on the crack surface density function per unit volume of the solid. The equation is so-called AT2 model because of the quadratic term in PDE.

We set sharp crack surfaces as Dirichlet boundary conditions in . Hence, the crack phase-field is obtained from the minimization of the regularized crack density function as


Figure 2 gives the numerical solution that arises from the minimization Eq. 8 and demonstrates the effect of different regularized length scales on the numerical solution. Clearly, a smaller length scale leads to a narrower transition zone (see Figure 2c). That is also in agreement with the crack phase-field profile shown in Figure 1b.

Figure 2: Effect of different length scales on the crack phase-field resolution as calculated by the minimization problem in Eq. 8 such that .

2.4 Strain-energy decomposition for the bulk free-energy

Fracture mechanics is the process which results in the compression free state. As a result, a fracture process behaves differently in the tension phase and in compression phase, see e.g. aldakheel2014towards . In the following, an additive split for the strain energy density function to distinguish the tension and compression phases is used. Instead of dealing with a full linearized strain tensor , the additive decomposition

of the strain tensor based on its eigenvalues is used

MieWelHof10b ; HeWheWi15 . Herein, refers to the a Macaulay brackets for . Furthermore, and refer to the tension and compression parts of the strain, respectively. The are the principal strains (i.e., the eigenvalues of the ) and the

are the principal strain directions (i.e., the eigenvectors of the

). To determine the tension and compression parts of total strain , a positive-negative fourth-order projection tensor is


such that the fourth-order projection tensor projects the total linearized strain onto its tension-compression counterparts, i.e., . Hence an additive formulation of the strain-energy density function consisting of the tension and the compression parts reads


Here, the scalar valued principal invariants in the positive and negative modes are


2.5 Energy functional for the isotropic crack topology

Due to the physical response of the fracture process, it is assumed that the degradation of the bulk material due to the crack propagation depends only on the tensile and isochoric counterpart of the stored bulk energy density function. Thus, there is no degradation of the bulk material in compression mode, see MieWelHof10b . Hence, the degradation function denoted as acts only on the positive part of bulk energy given in Eq. 10, i.e.,


This function results in degradation of the solid during the evolving crack phase-field parameter . Due to the transition between the intact region and the fractured phase, the degradation function has flowing properties, i.e.,


Following MieWelHof10b , the small residual scalar is introduced to prevent numerical instabilities. It is imposed on the degradation function, which now reads


The stored bulk density function is denoted as . Together with the fracture density function , it gives the the total density function



Formulation 2.1 (Energy functional for isotropic crack topology).

We assume that and are given as well as initial conditions and . For the loading increments , find and such that the functional

is minimized.

Herein, to make sure that phase-field quantity  lies in the interval , we define to map negative values of to positive values. In Formulation 2.1, the stationary points of the energy functional are determined by the first-order necessary conditions, namely the Euler-Lagrange equations, which can be found by differentiation with respect to and .

Formulation 2.2 (Euler-Lagrange equations).

Let , be given as well as the initial conditions and . For the loading increments , find and such that


Herein, and are the first directional derivatives of the energy functional given in Formulation 2.1 with respect to the two fields, i.e., and , respectively. Furthermore, is the deformation test function and is the phase-field test function.

Furthermore, the second-order constitutive stress tensor with respect to Eq. 16 reads




2.6 Crack driving forces for brittle failure

Following aldakheel2018 ; teichtmeister2017 , we determine the crack driving state function to couple between two PDEs. Hence, crack driving state function acts as a right hand side for te phase-field equation. To formulate the crack driving state function, we consider the crack irreversibility condition, which is the inequality constraint imposed on our variational formulation. The first variation of the total pseudo-energy density with respect to the crack phase-field given in (15) reads


Herein, the functional derivative of with respect to is


Maximization the inequality given in Eq. 20 with respect to the time history reads


We multiply Eq. 22 by . Then Eq. 22 can be restated as


Here, denotes a positive crack driving force that is used as a history field from initial time up to the current time. Note that the crack driving state function is affected by the length-scale parameter and hence depends on the regularization parameter.

Input:    loading data on ,
            solution from time step .
[0.1cm] Initialization of alternate minimization scheme ():    set FLAG:=true while FLAG do
        1. set
        2. given , solve for , set ,
        3. given , solve for , set ,
        4. define alternate minimization residual for the obtained pair
 5. then
         end if
end while
Algorithm 1 Alternate minimization scheme for Formulation (2.3) at a fixed loading step .
Formulation 2.3 (Final Euler-Lagrange equations).

Let us assume that , are given as well as the initial condition and . For the loading increments , find and such that


The multi-field problem given in Formulation (2.3) depending on and implies alternately fixing and , which is a so called alternate minimization scheme, and then solving the corresponding equations until convergence. The alternate minimization scheme applied to the Formulation (2.3) is summarized in Algorithm 1.

3 Uncertainty quantification

In this section, we explain the source of uncertainty in phase-field models. Then, we introduce a computationally effective numerical technique to estimate the unknown parameters.

In the phase-field model, the uncertainty arises from the Lamé parameters including the shear modulus and Lamé’s first parameter as well as Griffith’s critical elastic energy release rate (material stiffness parameter) , which are assumed to be random fields. We represent the uncertainty by a spatially-varying Gamma random field. The field

can be characterized by its expectation and covariance using a Karhunen-Loéve expansion (KLE). Considering the probability density function

, the expansion is


which leads to the decomposition


Here the first term is the mean value,

are the orthogonal eigenfunctions,

are the corresponding eigenvalues of the eigenvalue problem ghanem2003stochastic


and the

are mutually uncorrelated random variables satisfying


where indicates the expectation of the random variables.

The infinite series can be truncated to a finite series expansion (i.e., an -term truncation) by ghanem2003stochastic


For the Gaussian random field, we employ an exponential covariance kernel as


where is the correlation length as well as

is the standard deviation.

For a random field, we describe the parameters using an KL-expansion. Considering the Gaussian field , a lognormal random field can be generated by the transformation . For instance, for the parameter , the truncated KL-expansion can be written as


3.1 Bayesian inference

We consider Formulation (2.3) as the forward model , where . The forward model explains the response of the model to different influential parameters (here , , and ). We can write the statistical model in the form smith2013uncertainty


where indicates a vector of observations (e.g., measurements). The error term

arises from uncertainties such as measurement error due experimental situations. More precisely, it is due to the modeling and the measurements and is assumed to have a Gaussian distribution of the form

with known covariance matrix . The error is independent and identically distributed and independent from the realizations. Here, for sake of simplicity, we assume (for a positive constant ).

For a realization of the random field corresponding to a realization of the observations , the posterior distribution is given by


where is the prior density (prior knowledge) and is the space of parameters (the denominator is a normalization constant) stuart2010inverse . The likelihood function can be defined as smith2013uncertainty


As an essential characteristic of the phase-field model, the load-displacement curve (i.e., the global measurement) in addition to the crack pattern (i.e., the local measurement) are appropriate quantities to show the crack propagation as a function of time. Figure 3 indicates the load-displacement curve during the failure process. Three major points are the following.

\⃝raisebox{-0.9pt}{1}First stable position. This point corresponds to the stationary limit such that we are completely in elastic region ().

\⃝raisebox{-0.9pt}{2}First peak point. Prior to this point crack nucleation has occurred and now we have crack initiation. Hence, this peak point corresponds to the critical load quantity such that the new crack surface appears (i.e., there exist some elements which have some support with ).

\⃝raisebox{-0.9pt}{3}Failure point. At this point failure of the structure has occurred and so increasing the load applied to the material will not change the crack surface anymore.

The interval between point 1 and point 2 in Figure 3 typically refers to the primary path where we are almost in the elastic region. The secondary path (sometimes referred to as the softening path) starts with crack initiation occurring at point 2. The whole process recapitulates the load-deflection curve in the failure process.

Figure 3: The schematic of load-deflection response for the failure process including primary path (prior to the crack initiation, i.e., between point 1 and 2) and secondary path (during crack propagation, i.e., between point 2 and 3).

The main aim of solving the inverse problem followed here is to determine the random field to satisfy (33). We strive to find a posterior distribution of suitable values of the parameters , , and in order to match the simulated values (arising from (25)) with the observations. The distribution provides all useful statistical information about the parameter.

Remark 3.1.

Note that the principal parameters , , and are linked. Specifically, and hold such that and through discretization error estimates. A simplified analysis for a decoupled, linearized, problem is provided in MaWi19 , and we set with and with . Hence, the equations to be solved through Formulation 2.3 are strongly related to the element size , i.e., the degradation function for displacement equation the and crack driving state in the phase-field equation. Hence, to resolve the crack phase-field which also tends to the sharp crack surface, a sufficiently small is chosen to obtain the reference solution (due to the lankness of the experiments). We also notice that from a theoretical side, the inequality should hold, which can be justified through convergence results, see e.g., AmTo92 ; Braides1998 .

The crack pattern is a time-dependent process (more precisely in a quasi-static regime, the cracking process is load-dependent), i.e., after initiation it is propagated through time. In order to approximate the parameters precisely, we estimate the likelihood during all time steps. Therefore, the posterior distribution maximizes the likelihood function for all time steps, and therefore we have an exacter curve for all crack nucleation and propagation times.

MCMC is a suitable technique to calculate the posterior distribution. When the parameters are not strongly correlated, the MH algorithm smith1993bayesian is an efficient computational technique among MCMC methods. We propose a new candidate (a value of ) according to a proposal distribution and calculate its acceptance/rejection probability. The ratio indicates how likely the new proposal is with regard to the current sample. In other words, by using the likelihood function (35), the ratio determines whether the proposed value is accepted or rejected with respect to the observation (here the solution of Formulation 2.3 with a very fine mesh). As mentioned, fast convergence means that the parameters are uncorrelated. A summary of the MH algorithm is given below.

Initialization: set prior data and number of samples .
for  do
        1. Propose a new candidate based on the proposal distribution .
        2. Compute the acceptance/rejection probability .
        3. Generate a random number .
        4. then
          accept the proposed candidate and set
          reject the proposed candidate and set
         end if
end for
Algorithm 2 The Metropolis-Hastings algorithm.

4 Bayesian inversion for phase-field fracture

In this key section, we combine the phase-field algorithm from Section 2 with the Bayesian framework presented in Section 3.

First, we define two sampling strategies as follows:

  • One-dimensional Bayesian inversion. We first use -samples (according to the proposal distribution) and extract the posterior distribution of the first unknown (e.g., ) where other parameters are according to the mean value. Then obtained information is used to estimate the posterior distribution of next unknowns (e.g., ); where here still is according to the mean. Finally, the expected value of obtained parameters (here and ) will be employed in the Bayesian inference (see Algorithm 2) to estimate .

  • Multi-dimensional Bayesian inversion. A three-dimensional candidate (, as -th candidate) is proposed (based on the proposal distribution) and the algorithm computes its acceptance/rejection probability.

To make the procedures more clarified we explain the multi-dimensional approach in Algorithm 3. Clearly, for the one-dimensional setting; for each parameter (e.g., ), it can be reproduced separately. We will study both techniques in the first example and the more efficient method will be used for other simulations.

for  do
        1. Propose the -th candidate according to the proposal distribution (unform or normal).
        2.  set FLAG=true
       while FLAG do
                (i) solve the Formulation (2.3)by Algorithm 1 considering and
                      the proposed candidate .
                (ii) approximate
                (iii) estimate the crack pattern at the loading stage by
   (iv) if   or   then
                    set FLAG=false
                   end if
       end while
       3. Calculate the likelihood function (35) for (during all -steps, until ) with respect to where
         indicates the reference value at the -th loading step.
        4. Compute the acceptance/rejection probability .
        5. Use Algorithm 2 to determine (i.e., is accepted/rejected).
end for
Algorithm 3 The multi-dimensional Bayesian inversion for phase-field fracture.

Here, is the sufficiently large value that is set by the user. Also, is a sufficiently small value to guarantee that the crack phase-field model reached to the material failure time. Note, in part (iv) for the while-loop step, the criteria in the secondary path (i.e., during crack propagation state) guarantees that reaction force under imposed Dirichlet boundary surface is almost zero. Hence, no more force exists to produce a fractured state. We now term this as a complete failure point. But, in some cases, e.g., shear test as reported in MieWelHof10b , by increasing the monotonic displacement load, is not reached to zero. For this type of problem, if holds then phase-field step (i.e., while-loop step) in Algorithm 3 will terminate.

5 Numerical examples

In this section we consider numerical test problems to determine the unknown parameters using given Bayesian inference. Specifically, we investigate:

  • Example 1. the single edge notch tension (SENT) test;

  • Example 2. double edge notch tension (DENT) test;

  • Example 3. tension test with two voids.

The observations can be computed by very fine meshes (here the reference values) as an appropriate replacement of the measurements (see Remark 3.1). Regarding the observational noise, is assumed. The main aim here is to estimate the effective parameters (, , and ) in order to match the load-displacement curve with the reference value. Furthermore, we assume that there is no correlation between , , and . To characterize the random fields, we can use the KL-expansion with and the correlation length as well.

The phase-field parameters set by , and regularized length scale (respecting the condition ). The stopping criterion for the iterative Newton method scheme, i.e. the relative residual norm that is


is chosen to . Here, indicates a discretized setting of weak forms described in Formulation (2.3). Regarding alternate minimization scheme we set for all numerical examples. Finally is chosen to guarantee that the we solve the model only until the material failure time.

5.1 Example 1. The single edge notch tension (SENT) test

This example considers the single edge notch tension. The specimen is fixed at the bottom. We have traction-free conditions on both sides. A non-homogeneous Dirichlet condition is applied at the top. The domain includes a predefined single notch (as an initial crack state imposed on the domain) from the left edge to the body center, as shown in Figure 4a. We set hence , hence the predefined notch is in the plane and is restricted to . This numerical example is computed by imposing a monotonic displacement at the top surface of the specimen in a vertical direction. The finite element discretization corresponding to is indicated in Figure 4b.

Figure 4: Schematic of SENT (Example 1) (left) and its corresponding mesh with (right).

For the shear modulus, we assume the variation range . Regarding the Lamé modulus , the parameter varies between and . Finally, we consider the interval between and for .

We solved the PDE model (Formulation 2.3) with , , and MieWelHof10b and the displacement during the time (as the reference solution) with was obtained. The main goal is to obtain the suitable values of , and such that the simulations match the reference value.

For this example, we use a uniformly distributed prior distribution and the symmetric proposal distribution




indicates the characteristic function of the interval


First, we describe the effect of each parameter on the displacement. As the Lamé constants (i.e., and ) become larger, the material response becomes stiffer; crack initiation takes longer to occur. Additionally, a larger crack release energy rate (as an indicator for the material resistance against the crack driving force) delays crack nucleation and hence crack dislocation. All these facts are illustrated in Figure 8.

Figure 8: The force/displacement curve for different values of (top left), (top right) and (bottom) in the SENT example (Example 1).
Figure 15: The prior (green line), the mean value (red line), and the normalized probability density function (pdf) of posterior distribution of (top left), (top right) and (bottom) for the SENT example. Here we compare the distributions obtained by the three-dimensional Bayesian inversion (first row) and the one-dimensional Bayesian inversion (second row).
one-dimensional three-dimensional
mean () ratio (%) mean () ratio (%)
84.1 24 85 24.5
120.2 30.3 119.1 30
29 28.5
Table 1: The mean values of the posterior distributions obtained by one-dimensional and three-dimensional Bayesian inversion in addition to their acceptance ratios for the SENT (Example 1). The units are in .
Figure 16: The load-displacement curve for the one-dimensional (black) and three-dimensional (red) posterior distributions in addition to the ones for the prior distribution (green) and the reference value (blue) for the SENT example (Example 1) with .

The posterior distributions obtained by both one- and three-dimensional Bayesian inversions are shown in Figure 15. The mean values of the distributions are and , and the acceptance rates are 24% and 24.5%, respectively. As we already expected (see Figure 8), the 1st parameter Lamé is not as influential as the shear modulus. This fact is validated by the posterior distribution, where wider distributions (with higher acceptance rates, approximately 30%) are obtained. Regarding the material stiffness parameter , the acceptance rates are near 29%. The values are summarized in Table 1.

To verify the parameters obtained by the Bayesian approach, we solved the forward model using the mean values of the posterior distributions. Figure 16 shows the load-displacement diagram according to prior and posterior distributions. As expected, during the nucleation and propagation process, using Bayesian inversion results in better agreement (compared to the prior). Furthermore, a better estimation is achieved by simultaneous multi-dimensional Bayesian inversion. From now onward, this approach will be used for Bayesian inference.

Figure 17: The autocorrelation function for one- and multidimensional Bayesian inference in the SENT example (Example 1).

5.1.1 The convergence of MCMC

A customary method to assess the convergence of the MCMC is the calculation of its autocorrelation. The lag- autocorrelation function (ACF) is defined as

where is the -th element of the Markov chain and indicates the mean value. For the Markov chains, is positive and strictly decreasing. Also, a rapid decay in the ACF indicates the samples are uncorrelated and mixing well. Figure 17 shows the convergence of the MCMC where the Lamé constants and the crirical elastic energy release rate are estimated with the one-dimensional MH algorithm. Also, we estimated the convergence observed in the multi-dimensional approach. For all parameters, MCMC converges fast indicating that they are not correlated. As expected, the multi-dimensional approach converges slower than the one-dimensional one.

Figure 24: The effect of the mesh size on the crack propagation in the SENT example (Example 1). The mesh sizes are (from the left) , , , , and (the reference). The effective parameters are chosen according to the prior values.
Figure 25: The load-displacement curve of in the SENT example (Example 1) for different mesh sizes, where the parameters are chosen according to the prior.

As noted above, the phase-field solution depends on and . A detailed computational analysis was for instance performed in HeWheWi15 ; HeiWi18_pamm . In general, for smaller (and also smaller ) the crack path is better resolved, but leads to a much higher computational cost.

Figure 24 illustrates the crack pattern using different mesh sizes varying between and . For these mesh sizes, we also show the load-displacement diagram in Figure 25.

Figure 26: The load-displacement curves for the reference, prior, and posterior values (left panel) for in the SENT example (Example 1). The marginal posterior distributions of the effective parameters are shown in the right panel.
Figure 27: The load-displacement curves for the reference, prior, and posterior values (left panel) for in the SENT example (Example 1). The marginal posterior distributions of the effective parameters are shown in the right panel.
Figure 28: The displacement curves for the reference, prior, and posterior values (left panel) for in the SENT example (Example 1). The marginal posterior distributions of the effective parameters are shown in the right panel.

Here we strive to solve the problem using a coarse mesh and employ MCMC to find parameters that make the solution more precise compared with the reference value. Figure 26 shows the obtained displacement with both prior and posterior distribution for . The efficiency of the Bayesian estimation is pointed out here since the peak point and the failure point are estimated precisely. The marginal posterior distributions are shown in the right panel as well. The estimation can also be performed for finer meshes: Figure 27 and Figure 28 illustrate the load-displacement curves for and , respectively. In both cases, in addition to the precise estimation of the crack-initiation point and the material-failure point, the curve is closer to the reference value. Again, the posterior distributions are shown on the right panels. Finally, the mean values of the posterior distributions in addition to their acceptance rates are indicated in Table 2.

rate (%) rate (%) rate (%)
84.0 24 134.5 15 0.00220 20
92.5 26 136.1 16 0.00270 25
88.0 25 119.8 28 0.00268 28
Table 2: The mean of the marginal posterior distributions of , , and in the SENT example (Example 1) for , , and . All units are in .

5.2 Example 2. Double edge notch tension (DENT) test

This numerical example is a fracture process that occurs through the coalescence and merging of two cracks in the domain. We consider the tension test with a double notch located on the left and right edge. The specimen is fixed on the bottom. We have traction-free conditions on both sides. A non-homogeneous Dirichlet condition is applied to the top-edge. The domain has a predefined two-notch located in the left and right edge in the body as shown in Figure 29a. We set and hence . For the double-edge-notches, let and with the predefined crack length of (Figure 29a). This numerical example is computed by imposing a monotonic displacement at the top surface of the specimen in a vertical direction. The finite element discretization uses is indicated in Figure 29b.

Figure 29: Schematic diagram for the DENT example (Example 2) (left) and its corresponding mesh with (right).
Figure 33: The force/displacement curve for different values of (top left), (top right) and (bottom) for the DENT example (Example 2).
Figure 43: The prior (green curve), the mean value (red line), and posterior (histogram) distribution of (the first line), (the second line) and (the third line) for the DENT example (Example 2). For the posterior distribution, we used samples (the first column), samples (the second column), and samples (the third column).

According to the truncated KL-expansion, for the Lamé modulus , Eq. 32 gives the mean value of and the standard deviation of . Therefore, the parameter varies between and . For the shear modulus, the expectation of and the standard deviation of leads to the variation range . Similarly, by using a KL-expansion for , we obtained the variation range between and . Figure 33 illustrates the effect of their different values on the curve, where similar to SENT the influence of the shear modulus and material stiffness parameter is larger than the Lamé modulus.

We assumed the symmetric proposal distribution, namely the normal distribution


Here we plan to study the effect of the number of samples on the posterior distribution. Figure 43 shows the distributions of , , and in addition to their prior distribution using , , and . The calculations are done with , and is used as the reference. As shown, with a larger number of samples, the distribution is close to a normal distribution. Table 3 points out the mean values in addition to the acceptance rate. For both Lamé constants, the suitable chosen prior distribution gives rise to a higher acceptance rate (more than for all candidates). However, for , the candidates tend to lower values (to converge more to the reference value), and therefore more candidates are rejected (around acceptance rate). Table 3 summarizes the mean values of the posterior distributions of all influential parameters.

Figure 51: The effect of the mesh size on the crack propagation in the DENT example (Example 2). The mesh sizes are (from the left) , , , , , and (the reference).
=3 000 =15 000 =50 000
mean () rate (%) mean () rate (%) mean () rate (%)
8.25 31.1 8.22 32.0 8.21 32.8
11.91 31.0 11.93 31.8 11.65 36.2
15.0 15.4 15.3
Table 3: The mean value and the acceptance rate of the posterior distributions of , , and with , , and in the DENT example (Example 2). All units are in .

As the next step, we use different mesh sizes for the Bayesian inversion. Figure 51 shows the crack pattern using different mesh sizes changing from to . Finer meshes lead to a smoother and more reliable pattern. Figure 54 depicts the load-displacement diagram using the prior values. With coarse meshes, the curve is significantly different from the reference including crack initiation. Using Bayesian inversion (see Figure 43) enables us to predict the crack propagation and initiation more precisely. As the figure shows, even for the coarsest mesh (compare to ) the peak and fracture points are estimated precisely. For finer meshes (e.g., ) the diagram is adjusted tangibly compared to the reference value. Finally, a summary of the mean values (of posterior distributions) and their respective acceptance rate is given in Table 4.

rate (%) rate (%) rate (%)
8.91 29.0 12.2 37.3 11.2
8.11 36.8 12.1 38.2 14.2
8.25 34.2 12.0 39.1 24.2
Table 4: The mean of the posterior distributions of , , and for different mesh sizes in the DENT example. The units are in

The significant advantage of the developed Bayesian inversion is a significant computational cost reduction. As shown, for SENT and DENT, by using Bayesian inference for coarser meshes, the estimated load-displacement curve is very close to the reference values. We should note that the needed CPU time for is approximately 4.5 hours; however the solution with is obtained in less than 10 minutes. This fact pronounces the computational efficiency provided by Bayesian inversion, i.e., obtaining a relatively preicse solution in spite of using much coarser meshes.