A shape optimization method for moving interface problems governed by the heat equation

doi:10.1016/j.cam.2020.113266

Journal of Computational and Applied Mathematics

Volume 390, July 2021, 113266

https://doi.org/10.1016/j.cam.2020.113266 Get rights and content

Abstract

The one dimensional Stefan problem is reformulated as a shape optimization problem for the position of the phase transition as a function of time. The functional to be minimized is the mismatch of the Dirichlet to Neumann map at the moving interface. We show that the minimizer is the only stationary point of the shape functional. A gradient based optimization method is derived using shape calculus. The state and adjoint equations of the heat equation are solved with integral equation techniques which avoid a discretization in the domain. A Nyström quadrature method is analyzed and numerical results are presented.

Introduction

Melting and solidification problems with either planar or spherical geometry arise in many applications of material science. Often, such processes are accurately modeled by a one dimensional Stefan problem. This is a well-studied topic and we refer to, e.g., [1], [2], [3] for an overview of the engineering and mathematical literature. Here, the evolution of the phase change is governed by the Stefan condition which states that the heat flux at the interface is proportional to the interface velocity.

The Stefan problem is an important example of a free surface problem where in addition to the temperature distribution the position of the phase transition is a-priori unknown. The literature on analytical and numerical methods for this type of problem is too extensive to review here in detail. To put this paper in context we only mention two different classes of numerical techniques and a few representative papers.

The first class is the front-tracking method. Assuming the solution for some time $t > 0$ has been computed, the Stefan condition is used to determine the velocity at $t$ , and hence the future position of the interface. This results in an ODE-type solver for the interface. While explicit solvers are fairly straightforward to implement, they are generally stable only for low order time stepping. Implicit solvers have better stability but are more difficult to implement because a nonlinear equation must be solved in each time step. Several variations of front tracking methods have appeared, the main difference is how the heat equation in the bulk is solved. This can be either finite differences [4], [5], [6], the heat balance integral method [7] or boundary integral equations [8], [9]. In this work we also employ integral equation methods to solve the parabolic PDE. This approach eliminates the need to discretize the domain and avoids dealing with remeshing the domain in each time step.

The second class of techniques for free surface problems is the level set method. Here the interface is implicitly defined by the zero-level of an auxiliary function. In the case of moving geometries, this function satisfies the Hamilton–Jacobi equation, a hyperbolic PDE, which is coupled with the parabolic equation for the temperature. This idea has generated considerable interest [10], [11], however, we do not consider it in this work as the Hamilton–Jacobi equation is not tractable with integral equation methods.

This paper approaches the Stefan problem from the viewpoint of shape optimization. Here the goal is to determine a domain that minimizes an objective function which is often constrained by a PDE [12], [13], [14]. There has been considerable interest to use shape optimization for stationary free boundary problems of elliptic problems, and in the convergence of their discretization. The papers [15], [16], [17] are representative examples. The literature on shape optimization for parabolic equations is much sparser. However, there are recent papers that deal with the reconstruction of the shape of an interior inclusion from time dependent, overdetermined Cauchy data on an exterior boundary [18], [19], [20].

Unlike the inverse shape reconstruction, the Stefan problem is not ill-posed, but is nevertheless challenging because of the time dependence of the geometry. In [21] an error functional was defined that assigns to an arbitrary evolution of the interface the mismatch between the Neumann data and the Stefan condition. The solution of the Stefan problem is the minimizer of this functional. In this article, we present a numerical method to compute the minimum via gradient-based optimization techniques. Hence we will derive the shape calculus to determine the derivative of the solution to the heat equation with respect to changing the time-dependent interface. Further, we will derive an adjoint equation to compute the shape gradient efficiently. With the help of the adjoint, we will show that the only stationary point of the error functional is the minimum solution. The practical consequence is that a steepest-descent type algorithm cannot terminate in a local minimum that is not the desired solution.

To solve the state and adjoint equation efficiently and with high accuracy we will derive an integral equation reformulation for boundary values of the heat equation. Here special attention must be given to additional terms that arise because of time dependent geometries. We present a Nyström quadrature method that handles singularities in the integral operators and the solution. Finally, we present numerical examples to confirm the practicality of our approach.

Section snippets

Mathematical formulation

We consider the classical Stefan problem, which describes the melting of a block of ice which initially occupies a half-space. Melting is induced by applying a heat source at the end. After non-dimensionalization of the space and time variables, the problem can be stated in terms of the heat equation for temperature $U (x, t)$ in the liquid phase $\partial_{t} U (x, t) - \partial_{x}^{2} U (x, t) = 0, x \in (0, R (t)), t \in (0, T),$ $\partial_{x} U (0, t) = - f (t), t \in (0, T),$ $U (R (t), t) = 0, t \in (0, T),$ $R^{'} (t) + \partial_{x} U (R (t), t) = 0, t \in (0, T),$ $R (0) = 0 .$ The goal is to determine the

Adjoint formulation

Choosing an $n$ -dimensional vector space and a basis for the discretization of $r$ and $\hat{r}$ reduces (2.4) to a nonlinear optimization problem in $n$ variables. The partial derivatives with respect to the basis coefficients are obtained by replacing $\hat{r}$ in (2.6) by each basis function. Thus the computation of the complete gradient for a given $r (t)$ involves solving (2.5) $n$ times, which has to be performed in each step of a gradient-based optimization method.

Fortunately this large computational overhead

Integral formulations

Since the state and adjoint equations have to be solved repeatedly in a minimization procedure it is important to use efficient solution methods. We will consider an integral approach, which is equivalent to the boundary integral method for parabolic problems in higher space dimensions. We will see that in the case of one space dimension the heat equation can be reduced to an Abel-type integral equation in time.

The free space Green’s function for the one dimensional heat equation is $G_{f} (x, y, t, τ) =$

Discretization of thermal potentials

The integral operator $V$ in (4.5) is a generalized Abel integral operator because it can be written in the form $V q (t) = \int_{0}^{t} \frac{1}{\sqrt{t - τ}} k (t, τ) q (τ) d τ$ with the kernel $k (t, τ) = \frac{1}{\sqrt{4 π}} [exp (- \frac{{(r (t) - r (τ))}^{2}}{t - τ}) + exp (- \frac{{(r (t) + r (τ))}^{2}}{t - τ})] .$ Since $r (t) - r (τ) = (t - τ) r_{1} (t, τ),$ $r (t) + r (τ) = 2 r (t) - (t - τ) r_{1} (t, τ),$ it follows that $k (t, t) = 1 ∕ \sqrt{4 π}$ as long as $t > 0$ . However, the limit of $k (t, 0)$ as $t$ approaches zero is $1 ∕ \sqrt{π}$ because $r (t) \sim t$ results in a non-trivial contribution of the second term. Thus $k (t, τ)$ is only smooth in the interior of the triangle $σ_{T} ≔ {(t$

Numerical results

We give more details about the discretization of the optimization problem (2.4). We seek an approximation $r_{N} (t)$ of the unknown function $r (t)$ in a finite dimensional vector space of functions that satisfy the initial condition $r (0) = 0$ . If $φ_{i} (t)$ is a basis for this vector space, then $r_{N} (t)$ and the perturbation $\hat{r} (t)$ are linear combinations of the basis functions $r_{N} (t) = \sum_{j = 0}^{N} r_{j} φ_{j} (t), {\hat{r}}_{N} (t) = \sum_{j = 0}^{N} a_{i} φ_{j} (t) .$ Thus the problem (2.4) is replaced by the optimization problem in $R^{N + 1}$ $R_{N} \equiv min_{[r_{0}, \dots, r_{N}] \in R^{N + 1}} J (r_{N}),$

Conclusion

We have formulated the classical Stefan problem as a shape optimization problem and solved the state and adjoint equation using an integral equation formulation. We were able to show that the shape functional is only stationary at the optimal solution, a property that is rarely seen in shape optimization. Our implementation indicates that the convergence of the interface position is consistent with the convergence of ansatz function space and quadrature order of the Nyström discretization of

References (25)

KutluayS. et al.
The numerical solution of one-phase classical Stefan problem
J. Comput. Appl. Math.
(1997)
CaseE. et al.
An integral equation method for spherical Stefan problems
Appl. Math. Comput.
(2012)
SethianJ.A. et al.
Crystal growth and dendritic solidification
J. Comput. Phys.
(1992)
ChenS. et al.
A simple level set method for solving Stefan problems
J. Comput. Phys.
(1997)
TauschJ.
Nyström Discretization of parabolic boundary integral equations
Appl. Numer. Math.
(2009)
CrankJ.
Free and Moving Boundary Problems
(1984)
GuptaS.C.
The Classical Stefan Problem, Basic Concepts, Modelling and Analysis with Quasi-Analytical Solutions and Methods
(2017)
RubinsteinL.I.
The Stefan Problem
(1971)
MeyerG.H.
One-dimensional parabolic free boundary problems
SIAM Rev.
(1977)
LiuF. et al.
A computationally efficient solution technique for moving-boundary problems in finite media
IMA J. Appl. Math.
(1999)

MitchellS.L. et al.

Application of standard and refined heat balance integral methods to one-dimensional Stefan problems

SIAM Rev.

(2010)

ShaH. et al.

Computation of the solidification of pure metals in plate geometry using Green’s function method

Int. J. Heat Mass Transfer

(1998)

Cited by (1)

On The Reformulation Of The Classical Stefan Problem As A Shape Optimization Problem
2022, SIAM Journal on Control and Optimization

^☆: This material is based upon work supported by the National Science Foundation, United States under grant DMS-1720431.

View full text

A shape optimization method for moving interface problems governed by the heat equation☆

Abstract

Introduction

Section snippets

Mathematical formulation

Adjoint formulation

Integral formulations

Discretization of thermal potentials

Numerical results

Conclusion

J. Comput. Appl. Math.

Appl. Math. Comput.

J. Comput. Phys.

J. Comput. Phys.

Appl. Numer. Math.

Free and Moving Boundary Problems

The Classical Stefan Problem, Basic Concepts, Modelling and Analysis with Quasi-Analytical Solutions and Methods

The Stefan Problem

One-dimensional parabolic free boundary problems

SIAM Rev.

A computationally efficient solution technique for moving-boundary problems in finite media

IMA J. Appl. Math.

Application of standard and refined heat balance integral methods to one-dimensional Stefan problems

SIAM Rev.

Computation of the solidification of pure metals in plate geometry using Green’s function method

Int. J. Heat Mass Transfer