Space-Time Mixed System Formulation of Phase-Field Fracture Optimal Control Problems

In this work, space-time formulations and Galerkin discretizations for phase-field fracture optimal control problems are considered. The fracture irreversibility constraint is formulated on the time-continuous level and is regularized by means of penalization. The optimization scheme is formulated in terms of the reduced approach and then solved with a Newton method. To this end, the state, adjoint, tangent, and adjoint Hessian equations are derived. The key focus is on the design of appropriate function spaces and the rigorous justification of all Fréchet derivatives that require fourth-order regularizations. Therein, a second-order time derivative on the phase-field variable appears, which is reformulated as a mixed first-order-in-time system. These derivations are carefully established for all four equations. Finally, the corresponding time-stepping schemes are derived by employing a dG(r\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r$$\end{document}) discretization in time.


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 [6,12] and overview articles and monographs include [7,8,11,31,32] 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,28,29,33], topology optimization with shape derivatives obtained with the adjoint method [10], stochastic phase-field modeling [14], and optimal control [15][16][17][25][26][27].Due to the irreversibility constraint on the fracture growth, optimization problems subject to such an evolution become mathematical programs with complementarity constraints (MPCC) [4,21,22] so that standard constraint qualifications like [30,34] cannot hold.In [17,[25][26][27] the complementarity constraint was replaced with a suitable penalty term.
Specifically, based on [26,27], a computational space-time framework for phasefield fracture optimal control problems was designed in [17], including detailed tests with six numerical examples.In particular, we employ Galerkin formulations in time and discuss in detail how the crack irreversibility constraint is formulated using a penalization [23,26] and an additional viscous regularization [19,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's method requires the evaluation of the state, adjoint, tangent, and adjoint Hessian equations.The latter requires the evaluation of second-order derivatives; see, e.g., [5], for parabolic optimization problems.Based on these settings, concrete time-stepping schemes are derived.As usual, the state and tangent equations run forward in time, whereas the adjoint and adjoint Hessian equations run backward in time.
However, in [17] the directional derivative of the irreversibility constraint was numerically approximated to obtain an efficient numerical scheme.The rigorous mathematical treatment is the key objective in the current work.As it will turn out, the justification of differentiability requires higher-order regularity in corresponding space-time function spaces.To this end, we start from an energy functional and propose suitable function spaces for the state and control variables.In addition, higher-order derivatives of the penalty functional are established.Using these results, the weak formulation of the state equation is derived.However, due to the crack irreversibility regularization, a second-order time derivative on the phase-field variable arises, which was not the case in [17] due to the previously mentioned numerical approximation.Dealing with this (higher-order) time derivative, we follow a classical approach (see, e.g., [3]), and formulate a mixed first-order-in-time system, which yields an additional equation.For that mixed system, the space-time formulation and the resulting timestepping scheme are derived.Finally, the weak formulation serves as constraint in the overall optimal control problem and, using the Lagrangian, the three additional equations are derived in a rigorous space-time fashion along with the resulting timestepping schemes.
The outline of this paper is as follows: In Sect.2, we introduce the phase-field fracture forward model and derive the derivatives of the fracture irreversibility penalization.In Sect.3, a Galerkin space-time discretization for the mixed system in its weak form is provided.Next, in Sect.4, we state the optimization problem, including the reduced space approach.In the key Sect.5, the Lagrangian and three auxiliary equations are derived in full detail.Our work is summarized in Sect.6.

Phase-Field Fracture Forward Model
To formulate the state equation, we first introduce some basic notation and then proceed with the energy functional.In addition, we establish higher-order derivatives for the penalty functional representing the crack irreversibility and we finally obtain a regularized energy functional.

Notation
We consider a bounded domain ⊂ R for the displacement field u and V ϕ := H 1 ( ) for the phase-field ϕ, and Q := L 2 ( N ; R 2 ) for the control force q, where Moreover we consider a compact time interval I = [0, T ] and introduce the spaces We denote natural scalar products and norms with the space as index, such as ( , L 2 scalar products and norms with the domain, such as

Energy Functional of Quasi-Static Variational Fracture Modeling
In the next step we introduce a functional E γ η ε : W × X → R from which we derive our forward problem.Here E γ η ε (q, u) is defined as the sum of the regularized total energy of a crack plus penalty and convexification terms for the time dependent irreversibility constraint φ ≤ 0. The regularized total energy of a crack is given by where q denotes a force that is applied in orthogonal direction to N ⊂ ∂ , C ∈ R 2×2×2×2 is the elasticity tensor and e(u) the symmetric gradient.Then, we have under certain assumptions where μ > 0, λ > −2/3μ are the Lamé parameters and I ∈ R 2×2 is the identity matrix.The so-called degradation function g(ϕ) := (1−κ)ϕ 2 +κ > 0 helps to extend the displacements to the entire domain .Due to the bulk regularization parameter κ > 0, the first term in E ε (q, u) is bounded away from zero even for ϕ = 0.The term is a regularized form of the Hausdorff measure [2].Herein G c > 0 denotes the critical energy release rate: for larger G c more energy is required to form a new fracture.The phase-field regularization parameter ε represents the width of the transition zone 0 ≤ ϕ ≤ 1, i.e., of the zone between fully broken material (ϕ = 0) and fully intact material (ϕ = 1).So far the problem consists in finding a function u ∈ X that minimizes the regularized total energy (1) subject to the irreversibility constraint φ(t, x) ≤ 0 a.e. in I × .

Crack Irreversibility Penalization and its Differentiability
In the sequel, the constraint φ(t, x) ≤ 0 is being replaced by a penalty term, which is defined as To ensure Fréchet differentiability in X ϕ up to third order, which will be necessary for the adjoint Hessian equation, we work with max{ φ, 0} 4 , although max{ φ, 0} 3+ρ for some ρ > 0 might suffice.Observe that f 4 L 4 (I × ) exists for every f ∈ X φ : from f (t) ∈ V ϕ it follows that f (t) ∈ L q ( ) for any q ∈ [1, ∞) due to the continuous embedding V ϕ → L q ( ), and from f ∈ X φ we similarly have f ∈ C(I , V ϕ ).In summary, any f ∈ X φ is also in L 4 (I , L 4 ( )).We notice that this part differs from the other penalizations used in [26] and [17].In addition, we notice other theoretical and numerical work on penalizations for phase-field fracture in [13,24].Later we will minimize (1) by solving the corresponding necessary optimality conditions of first order.In order to compute the Fréchet derivatives of R, we first prove an estimate for the L 2 norm of a product and then we use the directional derivative of the maximum operator h( f ) := max{ f , 0} in direction g: Differentiating R(ϕ) for given directions ψ i leads to L 2 norms of products of time derivatives ψi .In order to prove Fréchet differentiability of R we have to bound these L 2 norms in a proper way.The following proposition gives the required estimate with respect to the X ϕ norm.
Proof For m = 1, ψ I × ≤ ψ X ϕ holds by definition of the norms: For m > 1, set q := 2m and employ Hölder's inequality ( . The last norm is finite: ψi V ϕ ∈ H 1 (I ) holds by definition since ψi ∈ H 1 (I , V ϕ ), and the continuous embeddings . Now we apply Hölder's inequality (again with q = 2m) to the outer norm and bound the resulting L q (I ) norms with the L ∞ (I ) norm where Finally, we apply again the embedding which gives the desired result with

Proposition 2.2
The penalty functional R is three times Fréchet differentiable in X ϕ .For directions ψ 1 , ψ 2 , ψ 3 ∈ X ϕ the derivatives are given by Proof We have to show that the -linear functionals R ( ) (ϕ) are bounded with respect to the X ϕ norm and that the resulting remainder terms of the Taylor expansion satisfy r (ψ) ∈ o( ψ X ϕ ), where After deriving the given expressions for the functionals and showing their boundedness, we will actually prove the stronger estimates r (ψ) ∈ O( ψ +1 X ϕ ) using integral representations of r (ψ).

Convexification
One final modification of E ε is needed: we add the convexification term η 2 φ 2 I × for some η > 0. Indeed, in [27], the term η(ϕ i − ϕ i−1 , ψ) in the time steps t i−1 , t i was considered for η > 0, where ϕ i − ϕ i−1 is the finite difference approximation of φ and ψ is the test function.This term corresponds to a potential viscous regularization of a rate-independent damage model [19].We note that we extend this here to the continuous time case.

Regularized Energy Functional
Finally, the forward problem consists in finding u ∈ X that solves the following optimization problem for given initial data ϕ 0 , φ0 ∈ V ϕ and given control q ∈ W : with penalty parameter γ > 0 and convexification parameter η > 0. We notice that having fourth powers in R(ϕ) will yield third-order terms in time in the first-order necessary conditions on the PDE level.Conceptionally, the procedure is then similar to the elastic wave equation, see e.g., [3], and we borrow ideas for our formulations that will follow in the remainder of this paper.We notice that it holds 0 ≤ ϕ ≤ 1 (without imposing those bounds explicitly), by arguments from [27, Section 2].

Remark 2.1 (Convexification)
We notice that strict positivity η > 0 improves the numerical solution process of (4).For the quasi-static case of a related formulation with a fourth-order regularization and time-discrete convexification, it was proved in [27] that for sufficiently large values of η the control-to-state mapping associated with ( 3) is single valued due to strict convexity of the energy corresponding to the equation.However, the convexification term η 2 φ 2 I × also penalizes crack growth.To ensure the dominance of the physically motivated term γ R(ϕ), a careful weighting of γ and η is required with γ η.As it is often the case in numerics, the weighting is chosen heuristically or by experience since no rigorous numerical analysis exists.It is indeed a drawback of the (simple) penalization approach that several regularization parameters interact with material parameters and discretization parameters.

Space-Time Weak Formulation and Discretization
In this section, a weak formulation is first derived.Due to the second-order time derivative in the phase-field variable, we formulate an equivalent mixed first-orderin-time system, which differs from our prior work [17].These ideas are based on [3].Afterwards, our emphasis is on the Galerkin in time finite element discretization for this mixed system.Finally, we briefly account for the spatial discretization by employing finite elements, too.

Weak Formulation
Before we continue with the spatial discretization and the concrete time-stepping scheme, we state the weak form of (3).To this end, we replace (3) by its first-order optimality conditions ∂ u E γ η ε (q, u) = 0, see e.g., [26], yielding a coupled nonlinear PDE system: given ϕ 0 , φ0 ∈ V ϕ and q ∈ W , find u ∈ X such that every test function In order to derive the state equation from (4), we first combine the two equations into a single equation, To simplify the notation, let us collect the energy-related terms of (5) as follows.
Definition 3.1 The spatial semi-linear form a : For every subinterval J ⊆ I , a corresponding semi-linear form a J : W × X × X → R is defined as time integral over J , Next, we have to perform integration by parts in order to move the time derivatives in (5) from the test functions to the phase-field ϕ: and similarly Thus (5) becomes Additional focus has to be put on the second-order time derivative of ϕ.Following the approach of [3, Sect.3.1], we introduce an additional variable ϕ ∈ X φ that represents the first-order time derivative of the phase-field, ϕ (t, x) := φ(t, x).
(8) By this we replace every time derivative of ϕ in (7) with the new variable ϕ and attach the weak form of (8) as an additional equation.Consequently, ( 7) is replaced by a first-order mixed system in weak form: Given ϕ 0 , φ0 ∈ V ϕ and q ∈ W , find u ∈ X and ϕ ∈ X φ such that the following equations hold for all = ( u , ϕ ) ∈ X and ψ ∈ X ϕ : Remark 3.2 In contrast to [3, Eq. (3.16)], the starting point for our investigations is given by a weak PDE and not a strong formulation.Consequently, we define the mixed system and already include the initial values as in [3, Eq. (3.17)].The initial conditions for ϕ are specified in the first equation of ( 9) and the initial conditions for ϕ = φ in the second equation.

Galerkin Time Discretization
Using a time grid 0 By using the discontinuous Galerkin method, here dG(r ), we seek for a solution u in the space X r k of piecewise polynomials of degree r .The subindex k denotes the time-discretized function space in order to distinguish it from the continuous space X .To this end, we have To work with the discontinuities in X r k , we introduce the standard notation We keep following the approach of [3] with special focus on [3, Eq. (3.27)- (3.30)].Notice that in (9) we work with two terms involving time derivatives, and one of them even involves a maximum.Moreover, all time derivatives are only applied to ϕ and not to the full variable u as in [3].Since ϕ is in Y := X φ , we define the corresponding discretized space as Consequently, the semi-discretized state equation consists in finding u ∈ X r k and ϕ ∈ Y r k for a given control q such that for every ∈ X r k and ψ ∈ Y r k it holds Here the first line (10a) represents the dG(r ) discretization of the term while (10c) and (10d) result from the dG(r ) discretization of In (10b) and (10f) we have the respective initial conditions for ϕ and ϕ , and in (10e) the physical terms (energy).The last two lines (10g) and (10h) contain the temporal boundary terms due to the integration by parts that we did before.

Remark 3.3
In (10) the sums of jump terms start with m = 0, whereas in [3] they start with m = 1.The reason is that our later derivations follow [5] where the sums also start with m = 0.Both formulations are in fact equivalent, which follows from [20,Remark 3.2].
By shifting the index of the jump terms by one, all summations range from m = 1 to m = M and can therefore be combined into a single sum for each equation.This yields our final form of the full state equation.

Proposition 3.1 The state equation in dG(r ) form reads: Find u ∈ X r
k and ϕ ∈ Y r k for a given control q such that for every ∈ X r k and The discontinuous in time trial and test functions ( , ψ) ∈ X r k × Y r k allow for sequential decoupling from which we obtain the time-stepping scheme.

Proposition 3.2
The time-stepping for the state equation (10) starts with solving the initial conditions at m = 0: Then the equations for m = 1, . . ., M − 1 need to be solved in a forward recursion: Finally, we have to solve the terminal conditions at m = M:

Spatial Discretization
For the spatial discretization, we employ again a Galerkin finite element scheme by introducing H 1 conforming discrete spaces.We consider two-dimensional shape-regular meshes with quadrilateral elements K forming the mesh T h = {K }; see [9].The spatial discretization parameter is the diameter h K of the element K .On the mesh T h we construct a finite element space V h := V uh × V ϕh as usual: Herein Q s (K ) consists of shape functions that are obtained as bilinear transformations of functions defined on the master element K = (0, 1) 2 , where Qs ( K ) is the space of tensor product polynomials up to degree s in dimension d, defined as Specifically, for s = 1 and d = 2 we have Similar to the state variables, the control variables need to be discretized in time and space.The polynomial degrees of the discrete spaces may differ, and for the control we denote them by r and s .To this end, we employ a dG(r ) discretization in time and a cG(s ) discretization in space.The semi-discrete control space is denoted by W k and the fully discrete space is denoted by W kh .
We notice for the next two sections that the following derivations are independent of the specific spatial discretization and for this reason the subindex h is omitted.

Optimization with Phase-Field Fracture
We formulate the following nonlinear optimal control problem with a standard tracking type cost functional.For given ϕ 0 , φ0 ∈ V ϕ we seek a solution (q, u, ϕ ) 12), (13), and ( 14), (15) where ϕ d ∈ V ϕ is some desired time-independent phase-field and q d ∈ W is a suitable nominal control that can be used for numerical stabilization.The second norm represents a common Tikhonov regularization with the Tikhonov parameter α > 0.
Remark 4. 1 We notice that the existence of a global solution in W × X has been proved for a related problem in [26,Thm. 4.3].That problem differs from (15) in that it is posed in the primal form with a tracking type functional for the displacements u and a penalization in discrete time (also of order four).

Reduced Optimization Problem and Solution Algorithm
We solve ( 15) by a reduced space approach.To this end, we assume the existence of a solution operator S : W → X × Y , q → (u, ϕ ), via equation (9).With this solution operator the cost functional J (q, u) reduces to j : W → R, j(q) := J (q, S(q)).

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 ( 20)- (22).Specific emphasis is on the regularization terms for the crack irreversibility and the convexification.The Lagrangian is defined based on the dG(r ) formulation, and this Galerkin discretization yields naturally the adjoint, tangent and adjoint Hessian.Furthermore, this approach exhibits the property that optimization and discretization interchange.

Lagrangian
The Lagrangian L (10) is defined directly in the dG(r ) form as L r k (q, ū, z) := J (q, u) (24)

Adjoint Equation
In this section we derive a time-stepping scheme for the adjoint equation ( 20) from the discrete time Lagrangian (24).In the adjoint for dG(r ) we seek z = (z, When calculating the partial derivative of L with respect to ū, some care has to be taken with the temporal boundary values from prior integration by parts, i.e., with the last two lines in (24).We first formulate the adjoint equation in the continuous time weak form and later in the dG(r ) setting: Here a I u denotes the integral of a u over I , and the terms (25c), (25e), and (25f) arise from the integration by parts (i.b.p.) applied to the partial derivatives where BT denotes the corresponding boundary terms in (25c), (25e), and (25f).Note that the integration by parts is necessary in these cases because after computing the partial derivative the time derivative is applied to the direction ϕ or ψ and not to the adjoint z.Now we see that the last term in (25b) cancels with the second term in (25c), and so do the terms in (25e), (25f) with those in (25i), (25j).Consequently we obtain the adjoint equation in dG(r ) form as Herein a I u is the time integral of the partial derivative Next we exploit the separability property of the cost functional J to expand the time integrals (ψ, z ϕ ) I × and a I ū (q, u)( ¯ , z) into m subintegrals.After that we distinguish between the terms for m = 0, . . ., M − 1 and m = M and combine them.Finally, we obtain the desired dG(r ) form of the adjoint equation wherein J u,m denotes the partial derivative of J m with respect to u.

Proposition 5.1 The adjoint equation in dG(r ) form (where we write t M for T ) reads
This leads immediately to the resulting time-stepping scheme.

Proposition 5.2
The adjoint time-stepping for (27) starts with solving the terminal conditions at m = M for z(t M ) and z ϕ (t M ): Then the equations for m = M − 1, . . ., 1 need to be solved in a backward recursion: Finally, we have to solve the initial conditions at m = 0:

Tangent Equation
The second auxiliary equation is the tangent equation.In this equation we seek δ ū = (δu, δϕ ) ∈ X r k × Y r k such that Here we will apply the same procedure as for the state equation.Recall that L(q, ū, z) contains the integral a I (q, u)(z) with z entering linearly.Hence the partial derivative required for L ū z (q, ū, z)(δ ū, ¯ ) is simply a I u (q, u)(δu, ), and the partial derivative required for L q z (q, ū, z)(δq, ¯ ) is derived from (6) as the integral a I q (q, u)(δq, ): Furthermore, J (q, u) does not depend on z or z ϕ , hence J q z and J ū z vanish.In contrast to the adjoint we formulate the tangent equation directly including the dG(r ) jump terms: (δϕ , ψ) I m × − (δϕ(0), ψ(0)) As for the state equation we shift the summation index for the jump terms by one and combine all sums into a single one.

Adjoint Hessian Equation
The third and last auxiliary equation is the adjoint Hessian equation.In this equation we seek δ z = (δz, δz ϕ ) ∈ X r k × Y r k such that for all ¯ ∈ X r k × Y r k the following equation holds: First we see that L q ū (q, ū, z)(δq, ¯ ) = 0 since q and ū are decoupled.The derivative of a I in L z ū (q, ū, z)(δ z, ¯ ) is given by a I u (q, u)( , δz) due to the linearity of z in a I .However, a I uu arises as a genuine second-order derivative in L ū ū (q, ū, z)(δ ū, ¯ ) where Now we obtain the adjoint Hessian equation The following result guarantees that the least regular term (34b) is well-defined.

Proposition 5.5 The partial derivative ∂
Proof Similar to (2), an elementary calculation for ϕ ∈ The product rule together with the estimate |max{c, 0}| ≤ |c| for c ∈ R then yields The factors φ , δ φ are in L 2 (I , V ϕ ); all other factors in the right-hand side norms are in X φ → L ∞ (I , V ϕ ) while V ϕ = H 1 ( ) → L q ( ) for every q ≥ 2. Arguing as in Proposition 2.1, this gives the finiteness of the first product norm (and similarly of the two other ones): Remark 5.1 This proposition together with Proposition 2.2 shows why we chose a fourth order penalization for the fracture irreversibility constraint: this provides just sufficient regularity for the derivative ∂ t (max{ϕ , 0}ψδϕ ) to be in L 2 (I , L 2 ( )).A third order penalization would instead lead to nonexistent derivatives.
Finally, we apply dG(r ) to the terms involving δ żϕ and δ ż ϕ .We do not need to apply dG(r ) to (max{ϕ , 0} 2 ψδϕ , żϕ ) I × since it is part of L ū ū (q, ū, z)(δ ū, ¯ ): this acts as a right hand side as it does not involve the solution variable δ z.
with weights w m ≥ 0 can easily be added (or be used instead).In the four equations to be solved and in the associated time-stepping schemes this will produce objective terms similar to the current ones for m = 1, . . ., M plus additional terms at t 0 = 0.The Tikhonov regularization might also be formulated in discrete time, with contributions q(t m ) − q d (t m ) 2 .However, then we would have to replace W with C(I , Q) or a subspace thereof so that the point evaluations of q and q d are well-defined.

Conclusions
In this paper, we established rigorous space-time formulations for the state, adjoint, tangent, and adjoint Hessian equations of phase-field fracture optimal control problems.Due to the crack irreversibility constraint a higher-order penalization was adopted that allowed us to prove Fréchet differentiability.On the other hand, this higher-order term yielded a second-order time derivative of the phase-field variable.To this end, a mixed first-order-in-time system was formulated.Based on this system the Lagrangian and auxiliary equations were established.By using a discontinuous Galerkin discretization in time, the time integrals decouple and we derived the resulting time-stepping schemes.Possible extensions of the current work include the implementation as this mixed phase-field system has not been addressed in the numerical realization so far.Along with the implementation, comparisons of gradient methods, quasi-Newton approaches (BFGS), and the proposed full Newton method could be undertaken.However, the main computational challenges are the linear and nonlinear solvers of the forward problem as they become severely expensive due to numerous forward and backward runs within the optimization loop.This holds true for two-dimensional configurations, as shown in our prior work [17], as well as three-dimensional extensions.Conceptionally, numerically and implementation-wise three-dimensional configurations can be addressed by our model, while some mathematical statements must be carefully revisited, and the major challenge is the numerical cost.Another interesting extension benefits from the space-time formulation as the temporal part could be solved with higher-order in time Galerkin finite elements as for instance dG(1) methods in which linear elements per time interval are employed.A final modification can be the replacement of the penalty method by an augmented Lagrangian or primal-dual active set method, wherein the mathematical theory will be more involved.Numerically, we have implemented in other works (see references in the monograph [31,Chapter 5]) such alternative penalizations of the crack irreversibility constraint.The broader impact of the current work is reflected by the fact that in phase-field fracture numerous papers on forward modeling have been published to date, but due to mathematical and numerical challenges only a few have been published on optimization (which are cited in our introduction).Consequently, the present work provides one possibility of formulating phase-field fracture optimization problems.

2 .
The boundary is partitioned as ∂ = D .∪ N where both D and N have nonzero Hausdorff measure.Next we define Hilbert spaces