1 Introduction
Fracture propagation using variational approaches and phasefield 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 phasefield fracture, more recently some work started on parameter identification employing Bayesian inversion [18, 36, 28, 29], stochastic phasefield 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 phasefield 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 phasefield 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 nonconvexity of the governing energy functional of the forward phasefield 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], quasimonolithic 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 controltostate operator. Therein, Newtontype methods require the evaluation of the adjoint, tangent, and adjoint Hessian equations. The latter requires the evaluation of secondorder 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 timestepping 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 nonpropagating fracture is subject to an optimal boundary control.
The outline of this paper is as follows: In Section 2, the phasefield 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 Phasefield 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 displacementfield and the phasefield , and for the control , where
Moreover we consider a bounded time interval and introduce the spaces
On respectively we define the scalar products
with induced norms and , and furthermore the restricted inner products
with induced seminorms and . We also notice that we later work with , defined like , and with a semilinear form in which the first argument is nonlinear and the second argument is linear.
2.2 Energy functional of quasistatic 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
(1) 
where denotes a force that is applied in orthogonal direction to ,
is the elasticity tensor and
the symmetric gradient. Then, we havewhere are the Lamé parameters and
is the identity matrix. The socalled 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 asIn order to ensure differentiability up to second order, an alternative is to work with a fourthorder 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 rateindependent damage model [19].
Finally the forward problem consists in finding that solves the following optimization problem for given intial data and given control :
(2) 
with penalty parameter and convexification parameter .
Remark 2.1 (Initial condition ).
Note that we are concerned with quasistatic 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 quasistatic case that for sufficiently large values of the controltostate 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 timestepping 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
(3)  
for every test function .
2.4 Galerkin time discretization
Using a time grid we first partition the interval into leftopen subintervals ,
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 timediscretized function space in order to distinguish from the continuous space . To this end, we have
Remark 2.3.
Since we work with , i.e., constant functions in time, we have
To work with the discontinuities in , we introduce the notation
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 energyrelated terms of (3) with a semilinear form ,
(4)  
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
(5a)  
(5b)  
(5c)  
(5d) 
The time integral in (5c) has been approximated by the rightsided 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
(6) 
Now, since we are employing a dG(0) scheme, our test functions satisfy
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
(7)  
2.5 Timestepping scheme
We begin the solution process by solving the last line of (7):
(8)  
or equivalently . Then we proceed and solve for and every the following equation:
(9)  
Remark 2.4 (Spatial discretization).
We notice that in this work we are only interested in the temporal discretization. A full spacetime 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 followup paper.
3 Optimization with phasefield fracture
We formulate the following separable NLP. For given we seek a solution of
(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 nonnegative and weakly semicontinuous.
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
(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).

State equation: given , find such that for all (3) holds:
(12) 
Adjoint equation: given and , find such that for all
(13) 
Tangent equation: given , and a direction , find such that for all
(14)
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:
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
(16)  
Note that we have scaled the initial conditions with two different parameters and . For the phasefield 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 .
4.2 Adjoint
In the adjoint for dG(0) we seek such that
The first interesting part is the calculation of the derivative of . We formulate it directly in the weak form
(17)  
Herein the partial derivative of reads
(18)  
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
(19)  
At this point we have to decide how to approximate the time derivative . While for is easily approximated as
this procedure will not work for the first mesh point . The forward difference
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):
(20)  
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 phasefield 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:
(21)  
4.3 Adjoint timestepping scheme
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 :
(22)  
Now we collect what is left, multiply by and use the property ():
(23)  
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
Now we solve for the equation