On the numerical solution of a time-dependent shape optimization problem for the heat equation

TL;DR

This paper presents a gradient-based numerical method to identify a time-varying inclusion within a domain governed by the heat equation, using boundary measurements to minimize Neumann data mismatch.

Contribution

It introduces a novel approach for solving a time-dependent shape optimization problem involving heat equations, with validated numerical results.

Findings

Successful detection of time-dependent inclusions from boundary data

Validation of the gradient-based optimization method through numerical experiments

Effective minimization of Neumann data mismatch for shape identification

Abstract

This article is concerned with the solution of a time-dependent shape identification problem. Specifically we consider the heat equation in a domain, which contains a time-dependent inclusion of zero temperature. The objective is to detect this inclusion from the given temperature and heat flux at the exterior boundary of the domain. To this end, for a given temperature at the exterior boundary, the mismatch of the Neumann data is minimized. This time-dependent shape optimization problem is then solved by a gradient-based optimization method. Numerical results are presented which validate the present approach.

Equations

$$Q_T=\bigcup _{0<t<T}(\{t\}\times {\Omega }_t).$$

$${\Sigma }_T=\bigcup _{0<t<T}(\{t\}\times {\Gamma }_t)$$

$$\mathbf{\kappa }:[0,T]\times {\mathbb{R}}^d\to {\mathbb{R}}^d,(t,\mathbf{x})\mapsto \mathbf{\kappa }(t,\mathbf{x})$$

$$\Vert \mathbf{\kappa }(t,\mathbf{x}){\Vert }_{C^2([0,T]\times {\mathbb{R}}^d;{\mathbb{R}}^d)},\Vert \mathbf{\kappa }(t,\mathbf{x}{)}^{-1}{\Vert }_{C^2([0,T]\times {\mathbb{R}}^d;{\mathbb{R}}^d)}\le C_{\mathbf{\kappa }}$$

$$0<\underline{\sigma }\le \min \{\sigma (D\mathbf{\kappa })\}\le \max \{\sigma (D\mathbf{\kappa })\}\le \overline{\sigma }<\infty ,$$

$${\partial }_{t}u$$

$$u$$

$$u=f,$$

$$u(0,\cdot )$$

$${\partial }_{t}v$$

$$v$$

$$v(0,\cdot )$$

$$H^{r,s}(Q_{0}):=L^{2}\big{(}(0,T);H^{r}(\Omega_{0})\big{)}\cap H^{s}\big{(}(0,T);L^{2}(\Omega_{0})\big{)},$$

$$H^{r,s}(\Sigma_{0}):=L^{2}\big{(}(0,T);H^{r}(\Sigma_{0})\big{)}\cap H^{s}\big{(}(0,T);L^{2}(\Sigma_{0})\big{)}$$

$${\hat{H}}^{r,s}(Q_0)$$

$${\tilde{H}}^{r,s}(Q_0)$$

$${\hat{H}}^{r,s}({\Sigma }_0)$$

$$\displaystyle\hat{H}_{0}^{r,s}(Q_{0}):=\big{\{}u\in\hat{H}^{r,s}(Q_{0}):u|_{\Sigma_{0}}=0\big{\}},$$

$$\displaystyle\tilde{H}_{0}^{r,s}(Q_{0}):=\big{\{}u\in\tilde{H}^{r,s}(Q_{0}):u|_{\Sigma_{0}}=0\big{\}}.$$

$$\hat{H}^{-r,-s}(Q_{0})=\big{[}\tilde{H}_{0}^{r,s}(Q_{0})\big{]}^{\prime}\quad\text{for }r-\frac{1}{2}\notin\mathbb{Z}.$$

$$V(Q_{0}):=\big{\{}v=U|_{Q_{0}}:U\in C_{0}^{2}\big{(}(-\infty,T)\times\Omega_{0}\big{)}\big{\}}$$

$$H^{r,s}(Q_T):=\{v\in L^2(Q_T):v\circ \mathbf{\kappa }\in H^{r,s}(Q_0)\}$$

$$\Vert v\circ \mathbf{\kappa }{\Vert }_{H^s({\Omega }_0)}\sim \Vert v{\Vert }_{H^s({\Omega }_t)}.$$

$$({\partial }_t-\Delta )u$$

$$u$$

$$S_0:{\hat{H}}^{\frac{1}{2}+s,\left(\frac{1}{2}+s\right)/2}({\Sigma }_0)\to {\hat{H}}^{1,\frac{1}{2}}(Q_0)$$

$$S(v,\phi ):={\int }_0^T{\int }_{{\Omega }_t}\{\nabla v\cdot \nabla \phi +{\partial }_{t}v\phi \}\mathrm{d}\mathbf{x}\mathrm{d}t=0\text{for all}\phi \in {\tilde{H}}_0^{1,\frac{1}{2}}(Q_T).$$

$$({\partial }_t-\Delta )v$$

$$v$$

$$v(0,\cdot )$$

$$S(v,u)={\int }_0^T{\int }_{{\Omega }_t}hud\mathbf{x}dt,$$

$${\int }_0^T⟨{\partial }_t{v}^t(t),u^t(t){⟩}_{L^2({\Omega }_0)}+a(t;v^t(t),u^t(t))\mathrm{d}t={\int }_0^T⟨h^t(t),u^t(t){⟩}_{L^2({\Omega }_0)}\mathrm{d}t,$$

$$a(t;u^t,u^t)\ge \alpha \Vert u^t{\Vert }_{H^1({\Omega }_0)}^2-\lambda \Vert u^t{\Vert }_{L^2({\Omega }_0)}^2$$

$$a(t;u^t,u^t)$$

$$\displaystyle\quad-\int_{\Omega_{0}}\bigg{\|}\Big{(}\underbrace{(\operatorname{D}\bm{\kappa})^{-\intercal}\frac{1}{\det(\operatorname{D}\bm{\kappa})}\nabla\big{(}\det(\operatorname{D}\bm{\kappa})\big{)}}_{:=a_{1}}+\underbrace{\partial_{t}\bm{\kappa}}_{:=a_{2}}\Big{)}u^{t}\bigg{\|}\|(\operatorname{D}\bm{\kappa})^{-\intercal}\nabla u^{t}\|\,\mathrm{d}\mathbf{x}.$$

$$a(t;u^t,u^t)$$

Full text

1. Introduction

Shape optimization appears in a wide range of problems from engineering, especially for designing and constructing industrial components or in non-destructive testing. Many practical problems from engineering amount to partial differential equations for an unknown function, which needs to be computed to obtain the quantity of interest. Shape optimization is then concerned with the minimization of this quantity of interest. While shape optimization in case of elliptic partial differential equations is a well studied topic in literature, see for example [4, 22] and the references therein, not so much is known about shape optimization in case of parabolic partial differential equations.

Theoretical results for parabolic shape optimization problems with time-independent shapes can be found in [21, 22, 26], while practical results are found for example in [1, 2, 12]. This is in contrast to the results for parabolic shape optimization problems with time-dependent shapes. Theoretical results are for example available in [6, 7, 18], but to the best of our knowledge, no results about efficient computations of such time-dependent shape optimization problems exist.

This article is based on the previous article [12] by two of the authors, where a parabolic shape optimization problem is considered for a time-independent shape. The goal therein was to detect a fixed inclusion or void of zero temperature inside a three-dimensional solid or liquid body by measurements of the temperature and the transient heat flux at the accessible outer boundary. Since the underlying shape calculus turned out to be rather standard due to the stationarity of the inclusion, the focus has been on the development of an efficient solver for the underlying heat equation. In contrast, in the present article, we now consider an inclusion, which changes its shape during time. Therefore, the shape calculus becomes the focus, while the numerical experiments are performed in two space dimensions and serve as a proof of concept.

The problem under consideration is reformulated as a shape optimization problem by means of a tracking-type functional for the Neumann data. Therefore, for given temperature at the exterior boundary, the mismatch of the Neumann data is minimized in a least-squares sense. Since we intend to apply a gradient-based optimization algorithm, we compute the shape gradient of this functional by means of the adjoint approach, which is known to reduce the computational effort. Then, we make a parametric ansatz for the inclusion and use a boundary element method to solve the heat equations for the primal state and the adjoint state. Numerical results validate that the present approach is feasible, leading to meaningful reconstructions.

The remainder of the article is organized as follows. In Section 2, we state the problem under consideration. Section 3 is dedicated to the time-dependent shape calculus of our functional. Section 4 shows how we can discretize our problem in the case of a void which is star-shaped for all points of time. In order to solve the heat equation on the current domain, Section 5 explains how to do this by using a boundary element method. Since the method parallels that of [12], this section only discusses the changes for the moving boundaries considered in this article. In order to illustrate the developed techniques, they are applied to the example shown in Section 6. Finally, in Section 7, we give some concluding remarks.

2. Problem formulation

2.1. Model problem

Let $D\subset\mathbb{R}^{d}$ with $d=2,3$ be a simply connected, spatial domain with boundary $\Gamma^{f}=\partial D$. Moreover, we have a time component, and thus the domain $(0,T)\times D$ forms a cylindrical domain, called the space-time cylinder. At every time $t\in[0,T]$, a simply connected subdomain $S_{t}\subset D$ with boundary $\Gamma_{t}=\partial S_{t}$ lies inside $D$ such that it holds $\operatorname{dist}(\Gamma^{f},\Gamma_{t})>0$ for all $t$. The difference domain is called $\Omega_{t}:=D\setminus\overline{S_{t}}$. Taking into account the time again, we thus consider tubes (i.e., non-cylindrical domains), which contain a void and are represented as

[TABLE]

The interior boundary of the tube $Q_{T}$ is called

[TABLE]

and the exterior boundary of the tube is called $\Sigma^{f}=(0,T)\times\Gamma^{f}$.111We assume that the exterior boundary $\Gamma^{f}$ does not depend on time, but this is no necessity for the shape calculus presented in the subsequent chapter. The topological setup is illustrated in Figure 1. It is in analogy to [12], but we consider an interior boundary $\Gamma_{t}$ which moves in time instead of a fixed, interior boundary $\Gamma_{0}$.

For every time step $t$, we assume to have a smooth $C^{2}$-diffeomorphism $\bm{\kappa}$, which maps the initial domain $\Omega_{0}$ onto the time-dependent domain $\Omega_{t}$. In accordance with [18], we write

[TABLE]

to emphasize the dependence of the mapping $\bm{\kappa}$ on the time, where we have $\bm{\kappa}(t,\Omega_{0})=\Omega_{t}$. Here, $\bm{\kappa}\in C^{2}([0,T]\times\mathbb{R}^{d})$ and, as in [11, pg. 826], we assume the uniformity condition

[TABLE]

for some constant $C_{\bm{\kappa}}\in(0,\infty)$. To reduce the technical level of the ensuing discussion, we assume that $\Omega_{0}$ has $C^{2}$-smooth boundaries which implies that the boundaries of $\Omega_{t}$ have the same regularity.

Remark 2.1**.**

Notice that, due to the uniformity condition (2.2), we have as in [11]

[TABLE]

where $\sigma(.)$ denote the singular values. Moreover, as in [11, Remark 1, pg. 827], we assume $\det(\operatorname{D}\bm{\kappa})$ to be positive. The smoothness of the mapping also implies that the time derivative $\partial_{t}\bm{\kappa}$ is uniformly bounded.

We shall consider the following, overdetermined initial boundary value problem for the heat equation, where $f$ and $g$ are defined at the fixed exterior boundary $\Sigma^{f}$

[TABLE]

Here, $\mathbf{n}$ denotes the normal pointing outward of the domain $\Omega_{t}$. In what follows, we assume that $f$ vanishes for $t=0$, which implies the compatibility with the initial condition. We then seek the free boundary $\Sigma_{T}$, such that the overdetermined problem (2.3) allows for a solution $u$. In [1, Theorem 1.1], the uniqueness of such a boundary $\Sigma_{T}$ is proven in the case of a time-independent boundary. In view of the bijective mapping $\bm{\kappa}$ (2.1), this uniqueness result also holds in the time-dependent case.

2.2. Reformulation as a shape optimization problem

The task of finding the unknown boundary $\Sigma_{T}$ is reformulated as a shape optimization problem by introducing the function $v$ as the solution of the initial boundary value problem with Dirichlet boundary conditions for the heat equation

[TABLE]

We set $Q_{0}=(0,T)\times\Omega_{0}$, which has two time-independent boundaries denoted by $\Sigma_{0}:=(0,T)\times\partial\Omega_{0}$. The appropriate function spaces for parabolic problems in time invariat domains are the anisotropic Sobolev spaces, defined by

[TABLE]

see, e.g., [1, 3, 16]. Likewise, the corresponding boundary spaces are

[TABLE]

which are defined for $C^{2}$-boundary when $r\leq 2$. With these definitions at hand, we can moreover define

[TABLE]

As in the elliptic case, we can include also (spatial) zero boundary conditions into the function spaces by setting

[TABLE]

The dual spaces are denoted by $r,s\leq 0$ and we especially have

[TABLE]

Finally, we introduce the test space

[TABLE]

as in [1], which is a dense subspace of $\tilde{H}_{0}^{1,\frac{1}{2}}(Q_{0})$ [1] (for a $C^{\infty}$-boundary, see for example [14, Remark 2.2 on pg. 8]).

We are now in the position to introduce the non-cylindrical analogues of the above spaces by setting

[TABLE]

and likewise for all the other spaces, where the composition with $\bm{\kappa}$ only acts on the spatial component. Due to the chain rule, $v\circ\bm{\kappa}$ and $v$ have the same Sobolev regularity, provided that the mapping $\bm{\kappa}$ is smooth enough, see for example [17, Theorem 3.23] for the elliptic case. We especially have the equivalence of norms for $|s|\leq 2$

[TABLE]

For the cylindrical case it is well known that the solution operator $f\mapsto S_{0}f:=u$ of Dirichlet problem of the heat equation

[TABLE]

with homogeneous initial conditions is an isomorphism between the spaces

[TABLE]

for $s>-\frac{1}{2}$ when $\Omega_{0}$ is smooth and for $|s|<\frac{1}{2}$ when $\Omega_{0}$ is Lipschitz, see [16, Theorem 5.3] and [3, Proposition 4.13].

For the existence, uniqueness and regularity of solutions to (2.4), we have to make sure the analogous result also holds on a non-cylindrical domain. The main techique of the argument is to transport the heat equation to a parabolic problem with variable coefficients in the space-time cylinder $Q_{0}$ and apply the same functional analytic tools of the above references there.

Theorem 2.2**.**

There exists a unique solution $v\in\hat{H}^{1,\frac{1}{2}}(Q_{T})$ satisfying the boundary condition in (2.4) and

[TABLE]

Proof.

The assertion follows if we can show existence and uniqueness of the solution to the following generalization of problem (2.4)

[TABLE]

Its weak formulation reads

[TABLE]

where $S$ is given by (2.7). We set $u^{t}=u\circ\bm{\kappa}$ and similarly for $v^{t}$ and $h^{t}$.

Transforming (2.9) back to $Q_{0}$ by using Lemma A.2 with $\bm{\xi}=\bm{\kappa}$, $Q_{\varsigma}=Q_{T}$ and $Q_{\tau}=Q_{0}$ gives

[TABLE]

where $a$ is defined in Lemma A.2.

To show solvabilty of (2.9) we apply [15, Chapter 3, Theorem 4.1.], which requires boundedness and coercivity of $a$. The boundedness follows easily from Remark 2.1. It remains to show coercivity, that is, there exist some constants $\alpha>0$, $\lambda\in\mathbb{R}$, such that for almost all $t\in(0,T)$

[TABLE]

holds for all $u^{t}\in H^{1}_{0}(\Omega_{0})$. With the help of the Cauchy-Schwarz inequality, we have

[TABLE]

Completing the square gives

[TABLE]

Discarding the positive term and due to Remark 2.1, we have

[TABLE]

and, therefore, by using the parallelogram law

[TABLE]

Now we can apply again Remark 2.1 to $a_{1}$ and $a_{2}$ and the Poincaré-Friedrichs inequality to the first term to arrive at the desired estimate (2.10).

Secondly, following the lines of [3], the analogue of [3, Lemma 2.8] reads: For every $h\in\hat{H}^{-1,-\frac{1}{2}}(Q_{T})$, there exists a unique $v\in\hat{H}^{1,\frac{1}{2}}_{0}(Q_{T})$ satisfying $(\partial_{t}-\Delta)v=h$ in $Q_{T}$. For the proof, we can straightforwardly modify the proof [3, Lemma 2.8], which uses the adjoint operator and interpolation results.

Thirdly, due to the surjectivity of the trace operator, we can then follow the proof of [3, Theorem 2.9] to obtain the statement in the theorem. ∎

For the given state equation (2.4), we introduce the tracking-type functional for the Neumann data at the fixed boundary $\Sigma^{f}$

[TABLE]

This objective functional should be minimized in the space of admissible boundaries $\Sigma_{T}$. It is nonnegative, and it is zero and hence minimal if and only if $v=u$. The objective functional measures the $L^{2}$-error of the data mismatch and thus corresponds to the minimization in the least-squares sense.

3. Computation of the shape derivative

3.1. Shape calculus

In order to minimize the objective functional (2.11), we apply a gradient-based optimization method. To this end, we shall compute the shape derivative of the functional.

The shape calculus for time-dependent problems has been formulated by means of the speed method in [7] and [18]. The speed method allows for deformations which are not only small perturbations of the domain. One intends to find a velocity field $\mathbf{V}$, which generates the optimal tube. The solution $\mathbf{T}(t,\mathbf{\cdot}):\mathbf{x}\mapsto\mathbf{x}_{t}=\mathbf{T}(t,\mathbf{x})$ of the differential equation [27, pg. 6]

[TABLE]

describes the pathline of an individual particle being exposed to the velocity field $\mathbf{V}$. Hence, if we would inject a drop of dye at a certain point and time, and we do a time-lapse photography, we would see the pathline [23]. In other words, when considering $t$ as the trajectory parameter, a fixed point $\mathbf{x}$ gets moved along the trajectory $\mathbf{x}_{t}=\mathbf{T}(t,\mathbf{x})$. The point $\mathbf{x}$ can be thought of as the Lagrangian (or material) coordinate, while $\mathbf{x}_{t}$ is the Eulerian (field) coordinate [22, pg. 49]. The speed method is favorable when considering the Eulerian setting [18].

For the Lagrangian setting, which we consider here, the perturbation of identity is preferable. The shape calculus for the perturbation of identity is briefly stated in [18] as well. For our computations, we shall exploit the bijective mapping $\bm{\kappa}$ from (2.1), which implies the mapping scheme displayed in Figure 2. With the mapping $\bm{\kappa}$ we can associate the velocity field

[TABLE]

which could be used for the speed method. Since the outer boundary $\Sigma^{f}$ of the tube is fixed, this vector field is zero in normal direction there.

In order to apply the traditional shape calculus, we would like to perturb the tube. To this end, we consider a vector field $\mathbf{Z}(t,\mathbf{x})$, which generates the perturbation of identity $\mathbf{I}+s\mathbf{Z}$. It yields a new tube

[TABLE]

Notice that the perturbations under consideration are horizontal, meaning that we consider perturbations of $(t,\bm{\kappa})$ of the type $(0,\mathbf{Z})$, compare [18]. Moreover, $\mathbf{I}+s\mathbf{Z}$ should satisfy a uniformity condition as in (2.2).

3.2. Local shape derivative

As in the time-independent case, we can define non-cylindrical material and local shape derivatives. The material derivative $\dot{v}[\mathbf{Z}]$ is defined as

[TABLE]

while the local shape derivative $\delta v=\delta v[\mathbf{Z}]$ in the direction $\mathbf{Z}$ is given by

[TABLE]

Here, $v_{t,s}$ denotes the state computed on the perturbed domain $Q_{T}^{s}$ and $v_{t}$ the state computed on $Q_{T}$, see [18, pg. 166]. These two non-cylindrical derivatives are connected by the relation

[TABLE]

Theorem 3.1**.**

The local shape derivative of the state $v$ from (2.4) can be computed as the solution of the partial differential equation

[TABLE]

The proof of the local shape derivative is presented in Appendix A, where we reformulate the time-independent proof found in [1] for the time-dependent setting.

3.3. Shape derivative of the objective functional

With the local shape derivative at hand, we are in the position to compute the shape derivative of the objective functional (2.11), which is defined by

[TABLE]

Theorem 3.2**.**

The shape derivative of the objective functional (2.11) in the direction $\mathbf{Z}(t,\mathbf{x})\in C^{2}(\Sigma_{T})$ reads

[TABLE]

where the adjoint state $p$ satisfies also the heat equation, but reversal in time:

[TABLE]

Proof.

Since we are not perturbing the exterior boundary, we have $\mathbf{Z}={\bf 0}$ in a neighborhood of $\Sigma^{f}$. Therefore, we conclude

[TABLE]

In view of the adjoint state equation (3.4), we can reformulate the derivative of $J$ by

[TABLE]

To derive (3.3), we apply Green’s theorem and obtain

[TABLE]

Since the integrands are smooth enough, we can apply the Reynolds transport theorem (see [10, pg. 78] for example) to treat the domain integral. Recall that the velocity $\mathbf{V}$, which transports the initial domain through the space-time tube, is given by (3.1). In combination with the end and initial conditions of $p$ and $\delta v$, respectively, we thus obtain

[TABLE]

that is

[TABLE]

Hence, by inserting the boundary condition for $\delta v$, we finally arrive at the desired result (3.3). ∎

Note that the tracking-type functional for the Dirichlet data has been considered in the setting of the speed method in [18, pg. 36–46]. It leads also to the same local shape derivative and shape gradient as in the time-independent case derived in [12]. This is thus consistent with the formulae stated here in case of the tracking-type functional for the Neumann data.

4. Discretization of the shape optimization problem

For our numerical experiments, we consider a two-dimensional spatial domain with a star-shaped void. As only its boundary is of interest and the shape gradient is also defined as a boundary integral, it suffices to parametrize just the interior boundary. Moreover, we consider only boundary perturbation fields $\mathbf{Z}$, because these are the only relevant perturbation fields as (3.3) shows.

Our choice of parametrization of the interior moving boundary $\Sigma_{T}$ of $Q_{T}$ is

[TABLE]

where the time-dependent parametrization $\bm{\gamma}(t,\cdot):[0,2\pi)\to\Gamma_{t}$ employs polar coordinates

[TABLE]

Here, $w(t,\phi)$ denotes the time- and angle-dependent radius, given by

[TABLE]

with $L_{\ell}(t)$ being appropriate dilations and translations of the Legendre polynomials of degree $\ell$.

Finding the optimal tube now corresponds to determining the unknown coefficients $\alpha_{k,\ell}$ and $\beta_{k,\ell}$ of the parametrization. Hence, we have the following finite dimensional problem:

[TABLE]

Here, $Z_{N}$ is the finite dimensional ansatz space of parametrizations. To compute the discrete shape gradient, we hence have to consider the directions

[TABLE]

for all $\ell=0,\dots,N_{L}$ and $k=0,\dots,N_{K}$, and

[TABLE]

for all $\ell=0,\dots,N_{L}$ and $k=1,\dots,N_{K}-1$.

With the specific parametrization at hand, the discrete shape gradient with respect to the parameters $t$ and $\phi$ reads

[TABLE]

compare (3.3), where we plugged in the choices for the perturbation fields (4.1) and (4.2), respectively, and used the parametrization $\bm{\gamma}$ to compute the normal $\mathbf{n}$.

The integral in the shape gradient (3.3) is computed by using a trapezoidal rule in space and a trapezoidal rule with a singularity correction at the endpoint $t=T$ in time (see the next section for details). The Legendre polynomials are computed by using their three term recurrence formula as described in [20], and are normalized afterwards while the Fourier series is evaluated efficiently by the fast Fourier transform.

The gradient-based method of our choice is the quasi Newton method, updated by the inverse BFGS rule without damping, cf. [9]. A second order line search is applied to find an appropriate step size in the quasi Newton method. For an overview of possible other optimization algorithms in general, see [5, 8].

5. Solving parabolic boundary value problems

We briefly describe the numerical method for solving the state and adjoint equation by using a boundary integral formulation. Since this is the approach that was already taken in [12] for a fixed boundary, we focus in this section on the changes for the time dependent case.

Both, the state and the adjoint equation, are Dirichlet problems of the heat equation with homogeneous initial conditions. In the case of the adjoint equation this becomes apparent after the change of variables $t\mapsto T-t$.

The boundary integral approach has distinct advantages over domain based approaches, because it is not necessary to mesh a time dependent domain or consider the transported problem in a cylindrical domain. Instead, we solve the Green’s integral equation. For a time-dependent boundary, it has the form

[TABLE]

Here, $\mathcal{V}$ and $\mathcal{K}$ are the thermal single and double layer operators defined below, and $\phi$ is a solution to the source-free heat equation with homogeneous initial conditions. Time dependence of the surface appears in the normal trace, which is defined as

[TABLE]

where $\langle{\bf V},{\bf n}\rangle$ is the normal velocity of the surface. The extra term in the definition of $\gamma_{1}^{\pm}$ arises from the Reynolds transport theorem in the derivation of (5.1). Details can be found in [25].

For the discretization of (5.1), it is desirable to have a method that can be easily adapted to time-dependent geometries, hence we use the Nyström discretization method of [24, 25]. To that end, we write the thermal layer potentials in the form

[TABLE]

where

[TABLE]

and $\Gamma_{\tau}\cup\Gamma^{f}=\partial\Omega_{\tau}$, i.e., the union of the free and the fixed boundary.

The kernel in the above time-dependent surface potentials is the Green’s function of the ($d-1$)-dimensional heat equation. Thus, they may be regarded as Poisson-Weierstrass integrals defined on a surface instead of the usual plane. As in the planar case, these integrals are smooth functions in all variables when $0\leq\tau\leq t$. The limiting behavior of these functions as $\tau\to t$ is

[TABLE]

where $H(\cdot)$ is the mean curvature of the surface $\Gamma_{t}\cup\Gamma^{f}$, see [25].

Since the functions $V\phi$ and $K\phi$ are smooth, the integral operators in (5.3) and (5.4) have a $(t-\tau)^{-1/2}$ singularity, which suggests to use the trapezoidal rule with a singularity correction at the endpoint $t=\tau$. It is shown in [24] that the rule

[TABLE]

where $h$ is the time step length, $t_{j}=hj$ and

[TABLE]

has a quadrature error of $\epsilon_{h}=\mathcal{O}(h^{3/2})$. Here, the prime at the summation sign indicates that the $j=0$ term in the sum is multiplied by the factor $1/2$. For the double layer analogous result holds when the $\mu_{n}$-term is mulitplied by the curvature. A fully discrete version is obtained by approximating the surface integrals in (5.5) and (5.6) by a surface quadrature rule, usually a composite rule that integrates polynomials on triangular patches exactly. If the spatial mesh width $h_{s}$ satisfies $\sqrt{h_{s}}\sim h$ and the spatial rule has at least degree of precision two then the quadrature error in (5.8) can be preserved, see [24]. In the time dependent case, these rules are constructed on $\Gamma_{0}\cup\Gamma^{f}$ and then mapped to $\Gamma_{t}\cup\Gamma^{f}$.

For the state equation, the solution is smooth and the Nyström method based on the above quadrature is used to computed the normal trace of the solution. Thus the Neumann data at the quadrature nodes is computed from (5.1) by substituting the given Dirichlet data of (2.4). This gives approximate values of the shape functional (2.11) and the boundary condition in the adjoint state (3.4).

The next task is to compute the Neumann data in the shape gradient (3.3) by solving the adjoint state. As already observed in [12], the adjoint equation (after time transformation $t\mapsto T-t$) has a singularity at $\tau=0$ because the homogeneous initial condition is not compatible with the in general non vanishing Dirichlet condition at $t=0$.

It can be concluded from (5.7) and Green’s integral equation that the Neumann data has a $t^{-1/2}$-singularity at $t=0$. To preserve the $\mathcal{O}(h^{\frac{3}{2}})$ accuracy, the time quadrature rule (5.8) must be modified with singularity corrections on both endpoints. Since the normal velocity of the boundary does not appear in (5.7), the derivation and the weights of this rule are identical to the case of a steady boundary. Since this can be found in [12], it is not repeated here.

6. Numerical experiments

The exterior, fixed boundary is chosen as the mantle of the cylinder with radius 1, where its height corresponds to the time interval $(0,T)=(0,1)$. We choose $N_{t}=90$ time intervals and, for every time step, $N_{\bf x}=80$ spatial points. The void is depicted in Figure 3, where the time corresponds to the $z$-axis. It is discretized by the same number of time intervals and spatial points as the exterior boundary.

We first solve the forward problem to construct the desired Neumann data $g$. To this end, we choose the desired shape found in Figure 3 and choose the Dirichlet data $f(t,\cdot)=t$, which matches with the initial data $u(0,\cdot)=0$ in $\Omega_{0}$. In order to avoid an inverse crime, we use an indirect boundary element approach by solving the thermal single layer equation and then recover the Neumann data by applying the thermal adjoint operator. In addition, we add $1\%$ random noise to the synthetic data.

Now, we can tackle the inverse problem. For the parametrization of the interior boundary, we choose $16$ Fourier coefficients in space ($N_{K}=8)$ and $10$ Legendre polynomials in time ($N_{L}=9$), leading to $160$ design parameters in total. As an initial guess for the free inner boundary, we choose the cylinder of radius $0.3$. We perform 100 iterations in the optimization procedure and use a quasi Newton method updated by the limited memory inverse BFGS rule, where 10 updates are stored, see [19] for example.

[math]20$$40$$60$$80$$100[math]1$$2$$\cdot 10^{-3}IterationValue of functional

[math]20$$40$$60$$80$$100[math]1$$2$$3$$\cdot 10^{-2}Iteration$\ell^{\infty}$-norm of gradient

The histories of the functional (left) and of the shape gradient (right).

In Figure 4 on the left, the evolution of the shape gradient during the course of the minimization algorithm is shown, while on the right the evolution of the functional is displayed. In Figure 5, we can see the $\ell^{2}$-error in the shape coefficients corresponding to the shape error. We clearly observe convergence of the minimization algorithm.

[math]20$$40$$60$$80$$100$$2\cdot 10^{-2}$$4\cdot 10^{-2}$$6\cdot 10^{-2}$$8\cdot 10^{-2}$$0.1$$8\cdot 10^{-2}Iteration$\ell^{2}$-error of the shape

In Figure 6, we present the final reconstruction of the shape, where the wireframe corresponds to the exact shape and the solid shape is its reconstruction. It can be figured out that the final reconstruction of the shape is not very good at the starting time $t=0$ and the stopping time $T=1$. But in between, the shape is very well reconstructed.

7. Conclusion

In this article, we solved a time-dependent shape reconstruction problem by means of shape optimization. We computed the shape derivative of the tracking-type functional for the Neumann data with the help of the perturbation of identity. It turned out that this shape derivative coincides with the one when the void is time-independent. We also demonstrated by numerical experiments that it is indeed possible to reconstruct a time-dependent shape by the proposed approach. By restricting to star-shaped voids, we have been able to compute the error between the desired shape and the reconstructed shape. The convergence of the minimization algorithm has clearly been observed.

Appendix A Local shape derivative

The proof of the local shape derivative follows the lines of [1]. We state here the adjustment to the time-dependent setting.

We first present two general lemmas, which are used later. We consider a mapping $\bm{\xi}$, which maps a domain $\Omega_{\tau}$ to a domain $\Omega_{\varsigma}$ and satisfies a uniformity condition as in (2.2). We will use the lemmas for $\bm{\xi}=\bm{\kappa}$ and $\bm{\xi}=\mathbf{I}+s\mathbf{Z}$. Let us denote $Q_{\tau}=\cup_{\tau}\{\tau\}\times\Omega_{\tau}$ and analogously for $Q_{\varsigma}$ and the lateral area by $\Sigma_{\tau}$ or $\Sigma_{\varsigma}$, respectively.

Lemma A.1**.**

For $v$ smooth enough it holds

[TABLE]

and

[TABLE]

Proof.

By the chain rule, we can compute

[TABLE]

from where (A.1) follows immediately. Moreover, the multivariable chain rule yields

[TABLE]

since only the spatial component is affected by the composition with $\bm{\xi}$. Using (A.1), we get (A.2). ∎

Notice that the identities (A.1) and (A.2) are also stated in [18, pg. 43] in the setting of the speed method.

Lemma A.2**.**

Let $v\in\hat{H}^{1,\frac{1}{2}}(Q_{T})$ and $\varphi\in\tilde{H}_{0}^{1,\frac{1}{2}}(Q_{T})$. Then, the transport of

[TABLE]

from $Q_{\varsigma}$ to $Q_{\tau}$ gives the parabolic problem

[TABLE]

with

[TABLE]

where $v^{\tau,\varsigma}=v\circ\bm{\xi}$ and similarly for $\varphi^{\tau,\varsigma}$ and $h^{\tau,\varsigma}$.

Proof.

With the aid of Lemma A.1, the transport of (A.3) from $Q_{\varsigma}$ onto $Q_{\tau}$ gives

[TABLE]

Using Green’s first identity and the zero boundary condition yields

[TABLE]

Thus, in the strong formulation, we have when dividing by $\det(\operatorname{D}\bm{\xi})$ that

[TABLE]

Rewriting gives

[TABLE]

Testing again with a function $\varphi^{\tau,\varsigma}$ gives the weak formulation

[TABLE]

which can be reformulated by using the divergence theorem with vanishing boundary terms to

[TABLE]

From here, the claim follows immediately. ∎

In order to compute the local shape derivative, we first introduce the material derivative to (2.4) as stated in the following lemma.

Lemma A.3**.**

The material derivative of (2.4), which is defined as the limit

[TABLE]

exists in $\hat{H}^{1,\frac{1}{2}}_{0}(Q_{T})$ and satisfies

[TABLE]

where $S$ is given by (2.7) and

[TABLE]

Proof.

Let $v_{t,s}$ be the solution of the perturbed problem on $Q_{T}^{s}$, satisfying the same boundary conditions as stated in (2.4). As an immediate consequence of [13, Chapter IV, Theorem 9.1], the solution $v_{t,s}$ lies in $\hat{H}^{2,1}(Q_{T}^{s})$ under our smoothness assumptions. Notice that the increased regularity of the solution of the differential equation is needed for the boundary condition of the local shape derivative (3.2).

We have for the perturbed bilinear form

[TABLE]

that $S_{s}(v_{t,s},\varphi)=0$ for all $\varphi\in\tilde{H}_{0}^{1,\frac{1}{2}}(Q_{T}^{s})$. The existence and uniqueness of a solution follows as in Theorem 2.2 by using that the transformation $\bm{\kappa}+s\mathbf{Z}\circ\bm{\kappa}$ satisfies again a uniformity condition as stated in (2.2). With similar computations as in the proof of Lemma A.2, when setting $\bm{\xi}=\mathbf{I}+s\mathbf{Z}$, $\Omega_{\tau}=\Omega_{t}$ and $\Omega_{\varsigma}=\Omega_{t,s}$, the transformation of the integral in (A.8) back onto $\Omega_{t}$ reads

[TABLE]

where we have set $v^{t,s}:=v_{t,s}\circ(\mathbf{I}+s\mathbf{Z})$ and $\varphi^{s}$ analogously. We define this bilinear form on the unperturbed domain as

[TABLE]

where

[TABLE]

Note that the last term in the definition of $S^{s}(w,\varphi)$ is new in comparison with [1].

We conclude the following statement:

[TABLE]

for $v_{t,s}\in\hat{H}^{2,1}(Q_{T}^{s})$ is equivalent to

[TABLE]

for $v^{t,s}\in\hat{H}^{2,1}(Q_{T})$. Integration by parts, where we use the zero boundary values of the test function, and dividing by \det\big{(}\operatorname{D}(\mathbf{I}+s\mathbf{Z})\big{)} verifies that (A.9) is equivalent to the formulation

[TABLE]

Because of $S(v,\varphi)=0$ and $S^{s}(v^{t,s},\varphi)=0$, it holds

[TABLE]

We can therefore consider

[TABLE]

for the computation of the material derivative, where

[TABLE]

Herein, the second line is new in comparison with [1].

We reformulate the expression for $G_{s}(\varphi)$ the same way as in [1] and we arrive at

[TABLE]

where the last line is new in this time-dependent setting in comparison with the proof given in [1]. We now need to show that $G_{s}$ converges to $G$ stated in (A.7).

Clearly, $\varphi\mapsto G_{s}(\varphi)$ is a bounded linear functional on $\tilde{H}_{0}^{1,\frac{1}{2}}(Q_{T})$, i.e. G_{s}\in\Big{(}\tilde{H}_{0}^{1,\frac{1}{2}}(Q_{T})\Big{)}^{\prime}. Therefore, we can interchange the integration and the limes $s\to 0$. Especially, as in [1], we have

[TABLE]

and

[TABLE]

as $s\to 0$. Thus, it remains to compute

[TABLE]

By using the Neumann series, we have

[TABLE]

and therefore

[TABLE]

In order to conclude the convergence $G_{s}\to G$ as $s\to 0$, we need that $v^{t,s}$ converges to $v$ in $H^{1,0}(Q_{T})$. To this end, we transform the equations for $v$ and for $v^{t,s}$ to $Q_{0}$ by using the transformation $\bm{\kappa}$, yielding two differential equations similar to (A.10). Applying [13, Theorem 4.5 on pg. 166] implies the convergence of $v_{t,s}\circ(\mathbf{I}+s\mathbf{Z})\circ\bm{\kappa}$ to $v\circ\bm{\kappa}$ and thus, with the uniformity condition (2.2), also $v^{t,s}$ converges to $v$. Therefore, we have convergence of $G_{s}\to G$ as $s\to 0$ in the dual space of $\tilde{H}^{1,\frac{1}{2}}_{0}(Q_{T})$ as in [1], with $G(\varphi)$ as in (A.7).

Now, we can argue as in [1]: Since the solution operator is an isomorphism from $\hat{H}^{-1,-\frac{1}{2}}(Q_{T})$ to $\hat{H}^{1,\frac{1}{2}}_{0}(Q_{T})$ (see Theorem 2.2), the statement in Lemma A.3 is true. ∎

Having the material derivative for (2.4) at hand, we are finally in the position to prove the local shape derivative posed in Theorem 3.1.

Proof of Theorem 3.1.

Starting from the material derivative, we would like to compute the local shape derivative $\delta v$.

If we consider $v\in\hat{H}^{2,1}(Q_{T})$, we also have $\nabla v\in H^{1,\frac{1}{2}}(Q_{T})$ and $\Delta v\in L^{2}(Q_{T})$, as in [1]. This follows from $\bm{\kappa}$ being a diffeomorphism and from the time-independent case in [14, Proposition 2.3 on pg. 14 with $r=2$, $s=1$, $j=2$ and $k=0$]. Then for $\varphi\in V(Q_{T})$ (see (2.5) for the definition of the space), we have the same identity as in [1, pg. 859], namely

[TABLE]

Applying this identity and the divergence theorem to (A.7) yields

[TABLE]

where the boundary terms vanish due to the compact support of $\varphi$. Note that only the last term of the integrand differs from the computations in [1]. It holds

[TABLE]

and, therefore, we can apply the divergence theorem again to get

[TABLE]

Taking the two time derivatives together yields

[TABLE]

This is the same expression as in [1]. Thus, the local shape derivative satisfies the same partial differential equation as in [1] except for being in a space-time tube $Q_{T}$ instead a space-time cylinder $Q_{0}$. ∎