# Space-time formulation and time discretization of phase-field fracture optimal control problems

The purpose of this work is the development of space-time discretization schemes for phase-field optimal control problems. First, a time discretization of the forward problem is derived using a discontinuous Galerkin formulation. Here, a challenge is to include regularization terms and the crack irreversibility constraint. The optimal control setting is formulated by means of the Lagrangian approach from which the primal part, adjoint, tangent and adjoint Hessian are derived. Herein the overall Newton algorithm is based on a reduced approach by eliminating the state constraint. From the low-order discontinuous Galerkin discretization, adjoint time-stepping schemes are finally obtained.

• 2 publications
• 3 publications
• 25 publications
03/28/2022

### Computational performance studies for space-time phase-field fracture optimal control problems

The purpose of this work are computational demonstations for a newly dev...
04/20/2020

### Space-time finite element discretization of parabolic optimal control problems with energy regularization

We analyze space-time finite element methods for the numerical solution ...
10/01/2020

### A Space-Time Variational Method for Optimal Control Problems

We consider a space-time variational formulation of a PDE-constrained op...
03/23/2020

### Discontinuous Galerkin method for a distributed optimal control problem governed by a time fractional diffusion equation

This paper is devoted to the numerical analysis of a control constrained...
03/31/2020

### Unstructured Space-Time Finite Element Methods for Optimal Sparse Control of Parabolic Equations

We consider a space-time finite element method on fully unstructured sim...
07/05/2016

### Optimal control for a robotic exploration, pick-up and delivery problem

This paper addresses an optimal control problem for a robot that has to ...
02/28/2021

### A Certified Reduced Basis Method for Linear Parametrized Parabolic Optimal Control Problems in Space-Time Formulation

In this work, we propose to efficiently solve time dependent parametrize...

## 1 Introduction

Fracture propagation using variational approaches and phase-field methods is currently an important topic in applied mathematics and engineering. The approach was established in [13, 6] and overview articles and monographs include [7, 8, 35, 34, 12] with numerous further references cited therein. While the major amount of work concentrates on forward modeling of phase-field fracture, more recently some work started on parameter identification employing Bayesian inversion [18, 36, 28, 29], stochastic phase-field modeling [15], and optimal control [26, 27, 25].

The main objective of this work is to design a computational framework for the last topic mentioned, namely phase-field fracture optimal control problems. In prior work [26, 27] the emphasis was on mathematical analysis and a brief illustration in terms of a numerical simulation for a fixed fracture. However, computational details have not yet been discussed therein, but are necessary in order to apply and investigate the methodology for more practical applications such as propagating fractures. Due to the irreversibility constraint on the fracture growth, optimization problems subject to such an evolution become mathematical programs with complementarity constraints (MPCC) [2, 22, 23] so that standard constraint qualifications like [30, 37] cannot hold. Our computational approach requires stronger regularity and hence we replace the complementarity constraint with a suitable penalty term.

Designing a computational framework for phase-field fracture optimal control is novel and challenging because robust forward and optimization solvers are required. For the forward solver, as intensively discussed in the literature, the linear and nonlinear solutions are demanding because of the non-convexity of the governing energy functional of the forward phase-field fracture model and the relationship of discretization and regularization parameters. For the nonlinear solution various methods were proposed such as alternating minimization (staggered solution) [5, 10], quasi-monolithic solutions [16, 34], and fully monolithic schemes [14, 32, 33, 20, 31]. Nonetheless, monolithic solutions remain difficult and we add an additional viscous regularization term as originally proposed in [19] and used in our governing model from [27]. The optimization problem is formulated in terms of the reduced approach by eliminating the state variable with a control-to-state operator. Therein, Newton-type methods require the evaluation of the adjoint, tangent, and adjoint Hessian equations. The latter requires the evaluation of second-order derivatives; see, e.g., [3] for parabolic optimization problems. The last paper serves as point of departure for our approach in the current work. Specifically, we employ Galerkin formulations in time and discuss in detail how the crack irreversibility constraint is formulated using a penalization [24, 26] and an additional viscous regularization [27, 19]. Based on these settings, concrete time-stepping schemes are derived. As usual, the primal and tangent problem run forward in time whereas the adjoint and adjoint Hessian equations run backward in time. We notice that some preliminary results are published in the book chapter [17] wherein a non-propagating fracture is subject to an optimal boundary control.

The outline of this paper is as follows: In Section 2, the phase-field fracture forward model is introduced. Furthermore, a Galerkin time discretization is provided. Next, in Section 3, the optimization problem is stated, including the reduced space approach. In the key Section 4 the Lagrangian and three auxiliary equations are carefully derived in great detail. Our work is summarized in Section 5.

## 2 Phase-field fracture forward model

To formulate the forward problem, we first introduce some basic notation.

### 2.1 Notation

We consider a bounded domain . The boundary is partitioned as where both and have nonzero Hausdorff measure. Next we define two function spaces, for the displacement-field and the phase-field , and for the control , where

 H1(Ω;\R2) \coloneqq\defsetv∈L2(Ω;\R2)Dαv∈L2(Ω;\R2) ∀α∈N20, \absα≤1, H1D(Ω;\R2) \coloneqq\defsetv∈H1(Ω;\R2)v|ΓD=0.

Moreover we consider a bounded time interval and introduce the spaces

 X\coloneqq\defsetu=(u,φ)u∈L2(I,V),∂tφ∈L2(I,H−1(Ω)),W\coloneqqC(I,Q).

On respectively we define the scalar products

 (u,v) \coloneqq∫Ωu⋅vdx∀u,v∈V, (u,v)I \coloneqq∫I∫Ωu⋅vdxdt=∫I(u(t),v(t))dt∀u,v∈X,

with induced norms and , and furthermore the restricted inner products

 (u(t),v(t))\set∂tφ(t)>0 \coloneqq{(u(t),v(t)),∂tφ(t)>0,0,else, (u,v)\set∂tφ>0,I \coloneqq∫I(u(t),v(t))\set∂tφ(t)>0dt∀u,v∈X,

with induced semi-norms and . We also notice that we later work with , defined like , and with a semi-linear form in which the first argument is nonlinear and the second argument is linear.

### 2.2 Energy functional of quasi-static variational fracture modeling

In the next step we introduce a functional from which we derive our forward problem. Here is defined as the sum of the regularized total energy of a crack plus a penalty term for the time dependent irreversibility constraint . The regularized total energy of a crack is given by

 Eε(q;u,φ)\coloneqq12(g(φ)\Ce(u),e(u))I−(q,u)ΓN,I+GcΓε(φ), (1)

where denotes a force that is applied in orthogonal direction to ,

is the elasticity tensor and

the symmetric gradient. Then, we have

 \Ce(u)=σ(u)=2μe(u)+λtr(e(u))I,

where are the Lamé parameters and

is the identity matrix. The so-called degradation function

helps to extend the displacements to the entire domain . The term is a regularized form of the Hausdorff measure [1]. So far the problem consists in finding a function that minimizes the regularized total energy subject to the irreversibility constraint . In the sequel, the constraint is being replaced by a penalty term, which is defined as

 R(φ)\coloneqq\norm∂tφ2\set∂tφ>0,I.

In order to ensure differentiability up to second order, an alternative is to work with a fourth-order penalization [26]. One final modification of is necessary. We add the convexification term for some . Indeed, in [27], the term in time steps was considered for . This term corresponds to a potential viscous regularization of a rate-independent damage model [19].

Finally the forward problem consists in finding that solves the following optimization problem for given intial data and given control :

 minu Eγε(q;u,φ)\coloneqqEε(q;u,φ)+γ2R(φ)+η2\norm∂tφ2I, (2)

with penalty parameter and convexification parameter .

###### Remark 2.1 (Initial condition u0).

Note that we are concerned with quasi-static brittle fracture without explicit time derivative in the displacement equation. Nonetheless, we introduce for formal reasons . First, we can develop in an analogous fashion time discretization schemes for the overall forward model. Second, it facilitates the extension to problems in which the displacement equation does have a time derivative, such as dynamic fracture [9, 4]. Third, having allows for a monolithic implementation structure, and the system matrix for the initial condition is regular.

###### Remark 2.2 (Convexification).

We notice that strict positivity ensures the required regularity in time, . Moreover, it improves the numerical solution process of (3). In fact, one can show for the quasi-static case that for sufficiently large values of the control-to-state mapping associated with (2) is single valued due to strict convexity of the energy corresponding to the equation. However, the convexification term also penalizes crack growth. To ensure the dominance of the physically motivated term we have to choose .

### 2.3 Weak formulation

Before we continue with the spatial discretization and the concrete time-stepping scheme, we state the weak form of (2). To this end we replace (2) by its first order optimality conditions, see e.g., [26], yielding a coupled nonlinear PDE system: given and , find such that

 (g(φ)\Ce(u),e(Φu\optindex))I−(q,Φu\optindex)ΓN,I =0, (3) Gcε(∇φ,∇Φφ\optindex)I−Gcε(1−φ,Φφ\optindex)I+(1−κ)(φ\Ce(u):e(u),Φφ\optindex)I +γ(∂tφ,Φφ\optindex)\set∂tφ>0,I+η(∂tφ,Φφ\optindex)I =0,

for every test function .

### 2.4 Galerkin time discretization

Using a time grid we first partition the interval into left-open subintervals ,

 I=\set0∪I1∪⋯∪IM.

By using the discontinuous Galerkin method, here dG(0), we seek for a solution in the space of piecewise polynomials of degree . The subindex denotes the time-discretized function space in order to distinguish from the continuous space . To this end, we have

 X0k\coloneqq\defsetv∈Xv|Im∈\polynomials0(Im,V),m=1,…,M and v(0)∈V.
###### Remark 2.3.

Since we work with , i.e., constant functions in time, we have

 ∂tv=v−m−v+m−1=0∀v∈X0k and m=1,…,M.

To work with the discontinuities in , we introduce the notation

 v+m\coloneqqv(tm+),v−m\coloneqqv(tm−)=v(tm),[v]m\coloneqqv+m−v−m.

Now the discretized state equation can be derived from (3) by combining the two equations into a single equation (5). To simplify the notation let us replace the energy-related terms of (3) with a semi-linear form ,

 a(q,u)(Φ) \coloneqqg(φ)⋅(\Ce(u),e(Φu\optindex)) (4) +Gcε(∇φ,∇Φφ\optindex)−Gcε(1−φ,Φφ\optindex) +(1−κ)(φ⋅\Ce(u):e(u),Φφ\optindex)−(q,Φu:y)ΓN.

Here denotes the component of in . Now the fully discretized state equation consists of finding a function for a given control such that for every

 0 =M∑m=1[γ(∂tφ,Φφ\optindex)\set∂tφ>0,Im+η(∂tφ,Φφ\optindex)Im] (5a) +M−1∑m=0[γ([φ]m,Φ+φ\optindexm)\setφ−m+1>φ−m+η([φ]m,Φ+φ\optindexm)] (5b) +M∑m=1a(q(tm),u(tm))(Φ(tm))Δtm (5c) +(u−0−u0,Φ−u\optindex0)+(φ−0−φ0,Φ−φ\optindex0). (5d)

The time integral in (5c) has been approximated by the right-sided box rule, where . Since the functions in might be discontinuous, we have to add jump terms in the typical dG(0) manner, which are contained in (5b). By expanding these jump terms, (5b) (with index shifted by one) becomes

 M∑m=1[γ(φ+m−1−φ−m−1,Φ+φ\optindexm−1)\setφ−m>φ−m−1+η(φ+m−1−φ−m−1,Φ+φ\optindexm−1)]. (6)

Now, since we are employing a dG(0) scheme, our test functions satisfy

 Φ+m−1=Φ−m∀m=1,…,M.

Therefore the two terms containing in (6) cancel and (5a) vanishes entirely by creftype 2.3. Combining the resulting expression with (5c) and (5d), we finally rewrite (5) as

 0 =M∑m=1(γ[(φ−m,Φ−φ\optindexm)\setφ−m>φ−m−1−(φ−m−1,Φ+φ\optindexm−1)\setφ−m>φ−m−1] (7) +η[(φ−m,Φ−φ\optindexm)−(φ−m−1,Φ+φ\optindexm−1)] +a(q(tm),u(tm))(Φ(tm))Δtm) +(u−0−u0,Φ−u\optindex0)+(φ−0−φ0,Φ−φ\optindex0).

### 2.5 Time-stepping scheme

We begin the solution process by solving the last line of (7):

 (u−0,Φ−u\optindex0) =(u0,Φ−u\optindex0), (8) (φ−0,Φ−φ\optindex0) =(φ0,Φ−φ\optindex0),

or equivalently . Then we proceed and solve for and every the following equation:

 0 =γ(φ(tm),Φφ\optindex(tm))\setφ(tm)>φ(tm−1)+η(φ(tm),Φφ\optindex(tm)) (9) −γ(φ(tm−1),Φφ\optindex(tm−1))\setφ(tm)>φ(tm−1)−η(φ(tm−1),Φφ\optindex(tm−1)) +a(q(tm),u(tm))(Φ(tm))Δtm.
###### Remark 2.4 (Spatial discretization).

We notice that in this work we are only interested in the temporal discretization. A full space-time Galerkin finite element discretization with continuous finite elements using quadrilaterals in space (see, e.g., the classical textbook [11] on finite elements) will be discussed in a follow-up paper.

## 3 Optimization with phase-field fracture

We formulate the following separable NLP. For given we seek a solution of

 minq,u J(q,u)\qstq(q,u) solves (???) and (???) for m=1,…,M, (10)

where is some separable functional, . To simplify the notation we assume that for . The existence of a global solution of in has been shown in [26, Theorem 4.3] for functions that are non-negative and weakly semi-continuous.

### 3.1 Reduced optimization problem and solution algorithm

We solve (10) by a reduced space approach. To this end, we assume the existence of a solution operator via equation (3). With this solution operator the cost functional can be reduced to , . As a result we can replace (10) by the unconstrained optimization problem

 minq j(q). (11)

The reduced problem is solved by Newton’s method applied to , and hence we need computable representations of the derivatives and . The established approach in [3] requires the solution of the following four equations for the Lagrangian ; the concrete form is defined in (16).

1. State equation: given , find such that for all (3) holds:

 L′z(q,u,z)(Φ)=0. (12)
2. Adjoint equation: given and , find such that for all

 L′u(q,u,z)(Φ)=0. (13)
3. Tangent equation: given , and a direction , find such that for all

 L′′qz(q,u,z)(δq,Φ)+L′′uz(q,u,z)(δu,Φ)=0. (14)
4. Adjoint Hessian equation: given , , from (13), from (14), and a direction , find such that for all

 L′′qu(q,u,z)(δq,Φ)+L′′uu(q,u,z)(δu,Φ)+L′′zu(q,u,z)(δz,Φ)=0. (15)

Solving these equations in a special order (see for instance [3, 21]) leads to the following representations of the derivatives that we need for Newton’s method:

 j′(q)(δq) =L′q(q,u,z)(δq)∀δq∈W, j′′(q)(δq1,δq2) =L′′qq(q,u,z)(δq1,δq2)+L′′uq(q,u,z)(δu,δq2) +L′′zq(q,u,z)(δz,δq2)∀δq1,δq2∈W.

## 4 Lagrangian and auxiliary equations

In the following main section, we specify the previously given abstract formulations in detail. We first derive the Lagrangian and then the three auxiliary equations (13)–(15). Specific emphasis is on the regularization terms for the crack irreversibility and the convexification.

### 4.1 Lagrangian

We formulate the Lagrangian within the dG(0) setting as

 L(q,u,z) \coloneqqJ(q,u) (16) −γ(∂tφ,zφ\optindex)\set∂tφ>0,I−η(∂tφ,zφ\optindex)I −∫Ia(q(t),u(t))(z(t))dt −η0(u(0)−u0,zu\optindex(0))−η(φ(0)−φ0,zφ\optindex(0)).

Note that we have scaled the initial conditions with two different parameters and . For the phase-field variable we use the convexification parameter of its time derivative to obtain . This is common in the context of a dG(0) setting as it produces desired cancelations with the jump terms resulting from the discontinuities of the test functions. In contrast, the initial condition for has no physical meaning. Therefore we use a separate parameter to obtain . Later we choose .

In the adjoint for dG(0) we seek such that

 L′u(q,u,z)(Φ)=0∀Φ∈X0k.

The first interesting part is the calculation of the derivative of . We formulate it directly in the weak form

 L′u(q,u,z)(Φ) =J′u(q,u)(Φ) (17) −γ(∂tΦφ\optindex,zφ\optindex)\set∂tφ>0,I−η(∂tΦφ\optindex,zφ\optindex)I −∫Ia′u(q(t),u(t))(Φ(t),z(t))dt −η0(Φu\optindex(0),zu\optindex(0))−η(Φφ\optindex(0),zφ\optindex(0)).

Herein the partial derivative of reads

 a′u(q,u)(Φ,z) =((1−κ)φ2+κ)⋅(\Ce(Φu\optindex),e(zu\optindex)) (18) +2φ(1−κ)Φφ\optindex(\Ce(u),e(zu\optindex)) +Gcε(∇Φφ\optindex,∇zφ\optindex)+Gcε(Φφ\optindex,zφ\optindex) +(1−κ)(Φφ\optindex⋅\Ce(u):e(u),zφ\optindex) +2φ(1−κ)(\Ce(Φu\optindex):e(u),zφ\optindex).

Now the main problem is that the time derivatives are applied to the test function  as usual in the adjoint. Therefore we use integration by parts to shift the time derivatives over to . Then the second line in (17) becomes

 γ(Φφ\optindex,∂tzφ\optindex)\set∂tφ>0,I+η(Φφ\optindex,∂tzφ\optindex)I (19) +γ(Φφ\optindex(0),zφ\optindex(0))\set∂tφ(0)>0+η(Φφ\optindex(0),zφ\optindex(0)) −γ(Φφ\optindex(T),zφ\optindex(T))\set∂tφ(T)>0−η(Φφ\optindex(T),zφ\optindex(T)).

At this point we have to decide how to approximate the time derivative . While for is easily approximated as

 ∂tφ(tm)≈φ(tm)−φ(tm−1)tm−tm−1,

this procedure will not work for the first mesh point . The forward difference

 ∂tφ(0)≈φ(t1)−φ(t0)t1−t0

is a good choice because it simplifies the condition to and leads to desired cancelations in (20). Now we will repeat the procedure that we applied to the state equation. We approximate the time derivatives and add the jump terms (with shifted index) as we did in (5), obtaining expressions similar to (6):

 L′u(q,u,z)(Φ) =J′u(q,u)(Φ) (20) +M∑m=1[γ(Φ−φ\optindexm,z−φ\optindexm−z+φ\optindexm−1)\setφ−m>φ−m−1 +η(Φ−φ\optindexm,z−φ\optindexm−z+φ\optindexm−1)] −γ(Φ−φ\optindexM,z−φ\optindexM)\setφ(tM)>φ(tM−1)−η(Φ−φ\optindexM,z−φ\optindexM) +γ(Φ−φ\optindex0,z−φ\optindex0)\setφ(t1)>φ(t0)+η(Φ−φ\optindex0,z−φ\optindex0) +M∑m=1[γ(Φ−φ\optindexm−1,z+φ\optindexm−1−z−φ\optindexm−1)\setφ−m>φ−m−1 +η(Φ−φ\optindexm−1,z+φ\optindexm−1−z−φ\optindexm−1)] −M∑m=1a′u(q(tm),u(tm))(Φ(tm),z(tm))Δtm −η0(Φ−u\optindex0,z−u\optindex0)−η(Φ−φ\optindex0,z−φ\optindex0).

Since , we have and see that the first sum vanishes entirely. We also see that the terms in the fifth and the last line of (20) cancel. Moreover, we assume that in the initial step, and hence the term in the fifth line vanishes as well.

###### Remark 4.1 (Projection of the initial solution).

The assumption is numerically justified since at some initial phase-field solution is prescribed. From to an projection of the initial conditions is employed that conserves the crack irreversibility constraint.

By the above arguments we eliminate the second, third and fifth line of (20) and the second term of the last line, whereas the initial values for are still present:

 L′u(q,u,z)(Φ) =J′u(q,u)(Φ) (21) −γ(Φ−φ\optindexM,z−φ\optindexM)\setφ(tM)>φ(tM−1)−η(Φ−φ\optindexM,z−φ\optindexM) +M∑m=1[γ(Φ−φ\optindexm−1,z+φ\optindexm−1−z−φ\optindexm−1)\setφ−m>φ−m−1 +η(Φ−φ\optindexm−1,z+φ\optindexm−1−z−φ\optindexm−1)] −M∑m=1a′u(q(tm),u(tm))(Φ(tm),z(tm))Δtm −η0(Φ−u\optindex0,z−u\optindex0).

From here on we exploit the separable structure of . We start the solution process by pulling out from (21) every term associated with the last time point :

 a′u (q(tM)u(tM))(Φ(tM),z(tM))ΔtM (22) +γ(Φ−φ\optindexM,z−φ\optindexM)\setφ−m>φ−m−1+η(Φ−φ\optindexM,z−φ\optindexM) =J′u(q(tM),u(tM))(Φ(tM))∀Φ∈X0k.

Now we collect what is left, multiply by and use the property ():

 0 =M∑m=1[γ(Φ−φ\optindexm−1,z−φ\optindexm−1−z−φ\optindexm)\setφ−m>φ−m−1+η(Φ−φ\optindexm−1,z−φ\optindexm−1−z−φ\optindexm)] (23) +M−1∑m=1a′u(q(tm),u(tm))(Φ(tm),z(tm))Δtm −M−1∑m=1J′u(q(tm),u(tm))(Φ(tm)) +η0(Φ−u\optindex0,z−u\optindex0)∀Φ∈X0k.

To formulate the equations that are actually solved in every time step we want to rewrite the entire equation as a single sum. Therefore we shift down the index of the first sum (the jump terms), take out the terms for , and obtain

 0 =M−1∑m=1([γ(Φ−φ\optindexm,z−φ\optindexm−z−φ\optindexm+1)\setφ−m+1>φ−m+η(Φ−φ\optindexm,z−φ\optindexm−z−φ\optindexm+1)] +a′u(q(tm),u(tm))(Φ(tm),z(tm))Δtm −J′u(q(tm),u(tm))(Φ(tm))) +γ(Φ−φ\optindex0,z−φ\optindex0−z−φ\optindex1)\setφ−1>φ−0+η(Φ−φ\optindex0,z−φ\optindex0−z−φ\optindex1) +η0(Φ−u\optindex0,z−u\optindex0).

Now we solve for the equation

 a′u (q(tm),u(tm))(Φ(tm),z(tm))Δtm +γ(Φ−φ\optindexm,z−φ\optindexm−z−φ\optindexm+1)\setφ−m+1>φ