A shape optimization method for moving interface problems governed by the heat equation☆
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 has been computed, the Stefan condition is used to determine the velocity at , 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 in the liquid phase The goal is to determine the
Adjoint formulation
Choosing an -dimensional vector space and a basis for the discretization of and reduces (2.4) to a nonlinear optimization problem in variables. The partial derivatives with respect to the basis coefficients are obtained by replacing in (2.6) by each basis function. Thus the computation of the complete gradient for a given involves solving (2.5) 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
Discretization of thermal potentials
The integral operator in (4.5) is a generalized Abel integral operator because it can be written in the form with the kernel Since it follows that as long as . However, the limit of as approaches zero is because results in a non-trivial contribution of the second term. Thus is only smooth in the interior of the triangle
Numerical results
We give more details about the discretization of the optimization problem (2.4). We seek an approximation of the unknown function in a finite dimensional vector space of functions that satisfy the initial condition . If is a basis for this vector space, then and the perturbation are linear combinations of the basis functions Thus the problem (2.4) is replaced by the optimization problem in
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)
- et al.
The numerical solution of one-phase classical Stefan problem
J. Comput. Appl. Math.
(1997) - et al.
An integral equation method for spherical Stefan problems
Appl. Math. Comput.
(2012) - et al.
Crystal growth and dendritic solidification
J. Comput. Phys.
(1992) - et al.
A simple level set method for solving Stefan problems
J. Comput. Phys.
(1997) Nyström Discretization of parabolic boundary integral equations
Appl. Numer. Math.
(2009)Free and Moving Boundary Problems
(1984)The Classical Stefan Problem, Basic Concepts, Modelling and Analysis with Quasi-Analytical Solutions and Methods
(2017)The Stefan Problem
(1971)One-dimensional parabolic free boundary problems
SIAM Rev.
(1977)- et al.
A computationally efficient solution technique for moving-boundary problems in finite media
IMA J. Appl. Math.
(1999)
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
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.