# American Institute of Mathematical Sciences

November  2018, 23(9): 3595-3622. doi: 10.3934/dcdsb.2017216

## A space-time discontinuous Galerkin spectral element method for the Stefan problem

 Department of Mathematics, Florida State University, Tallahassee, FL, 32306, USA

* Corresponding author: Mark Sussman (sussman@math.fsu.edu)

Received  May 2017 Revised  June 2017 Published  September 2017

A novel space-time discontinuous Galerkin (DG) spectral element method is presented to solve the one dimensional Stefan problem in an Eulerian coordinate system. This method employs the level set procedure to describe the time-evolving interface. To deal with the prior unknown interface, a backward transformation and a forward transformation are introduced in the space-time mesh. By combining an Eulerian description, i.e., a fixed frame of reference, with a Lagrangian description, i.e., a moving frame of reference, the issue of dealing with implicitly defined arbitrary shaped space-time elements is avoided. The backward transformation maps the unknown time-varying interface in the fixed frame of reference to a known stationary interface in the moving frame of reference. In the moving frame of reference, the transformed governing equations, written in the space-time framework, are discretized by a DG spectral element method in each space-time slab. The forward transformation is used to update the level set function and then to project the solution in each phase back from the moving frame of reference to the fixed Eulerian grid. Two options for calculating the interface velocity are presented, and both options exhibit spectral accuracy. Benchmark tests indicate that the method converges with spectral accuracy in both space and time for the temperature distribution and the interface velocity. A Picard iteration algorithm is introduced in order to solve the nonlinear algebraic system of equations and it is found that just a few iterations lead to convergence.

Citation: Chaoxu Pei, Mark Sussman, M. Yousuff Hussaini. A space-time discontinuous Galerkin spectral element method for the Stefan problem. Discrete & Continuous Dynamical Systems - B, 2018, 23 (9) : 3595-3622. doi: 10.3934/dcdsb.2017216
##### References:

show all references

##### References:
The Stefan problem on a Cartesian domain $\Omega \in \mathbb{R}^2$. $\Omega^{s}_t$ and $\Omega^{l}_t$ are time-dependent domains at time $t$, while domain $\Omega$ is fixed. $\Gamma(t)$ is the moving interface at time $t$, and $\boldsymbol{n}=(n_1,n_2)$ is the interface normal pointing in the direction of interface movement. $\phi$ is the level set function introduced in section 2.2
A tessellation of the space-time slab ${\mathscr{E}}^n$ of a space-time domain: ${\mathscr{E}}= \Omega \times [t_0,T] \in \mathbb{R}^2$. The coordinates of a point $\bar{\boldsymbol{x}}$ are denoted by $\bar{\boldsymbol{x}}=(x_1,x_{2})$, where $x=x_1$ is the spatial variable and $x_{2}$ denotes the time variable. The square with bold lines is an example of a space-time slab ${\mathscr{E}}^{n}$, such that ${\mathscr{E}}^n={\mathscr{E}} \cap I_n$, where $I_n=[t_n,t_{n+1}]$. We assume that the solid phase is on the left hand side of the interface ${\mathscr{Q}}^n_{\Gamma}$, while the liquid is on the other side. The space-time element ${\mathscr{K}}^{n,s}_j$ is referred to as a regular cell in the solid phase. However, ${\mathscr{K}}^{n,sl}_{j+1}$ and ${\mathscr{K}}^{n,sl}_{j+2}$ are examples of a cut-cell due to the time evolving interface, ${\mathscr{Q}}^n_{\Gamma}$, defined on a nonconforming computational grid. In this example, the curve ${\mathscr{Q}}^n_{\Gamma}$ obtained by 4th order interpolation cuts ${\mathscr{K}}^{n,sl}_{j+1}$ into two subelements, i.e., a pentagon ${\mathscr{K}}^{n,s,\Gamma}_{j+1}$ and a triangle ${\mathscr{K}}^{n,l,\Gamma}_{j+1}$
Examples of a cut-cell in a space-time domain ${\mathscr{E}} \in \mathbb{R}^2$. The moving front is evolving from left to right in each case. The shape of the generated subelements could be triangle, quadrilateral or pentagon
Examples of a mapping technique applied to two kinds of a cut-cell in a space-time domain ${\mathscr{E}} \in \mathbb{R}^2$. In each case, the curve ${\mathscr{Q}}^n_{\Gamma}$ is mapped to a straight bold line ${\tilde{\mathscr{Q}}}^n_{\Gamma}$, which is either aligned with an inter-element boundary or cutting a regular cell into two subelements
Description of the normal probe method applied in a space-time slab ${\mathscr{E}}^n \in \mathbb{R}^2$ in the fixed frame of reference. The polynomial order, $p=(p^{(x)},p^{(t)})$, in each space-time element is chosen to be $(3,2)$. We assume that the solid phase is on the left hand side of the interface ${\mathscr{Q}}^n_{\Gamma}$ while the liquid is on the other side. The fictitious element in the solid phase at time $t_{n+1}$ is denoted as $e^{n+1,s}_{fic}$, and $e^{n+1,l}_{fic}$ is the fictitious element for the liquid phase at time $t_{n+1}$. The length of both fictitious elements is $h$, which is the same as the length of a regular space-time element, i.e., the length of ${\mathscr{K}}^{n,s}_{j-1}$. The number of Gauss-Lobatto nodes in each fictitious element is $p+1$, where $p=p^{(x)}$. The temperature gradient is evaluated at the node with symbol "$\Box$" in each fictitious element
Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order $p^{(x)}$ in the liquid phase. The polynomial order in time direction, $p^{(t)}$, is chosen to be the same as $p^{(x)}$. The comparison is made between the weak form method (46) and the normal probe method (50). The simulation is computed over the time $t=0$ to $t=0.5$ with the time step $\Delta t=3.84 \times 10^{-2}$
Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order $p^{(x)}$ in the liquid phase. The simulation is computed over the time $t=0$ to $t=0.5$. The computational domain is divided into $E^{(x)}=5$ in spatial direction, with $E^{(x),l}=4$ in the liquid in the beginning and $E^{(x),l}=2$ in the liquid at the end. The interface velocity is computed by the weak form method (46). The time step is chosen to be $\Delta t=3.84 \times 10^{-2}$
Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order $p^{(x)}$ in the liquid phase. The simulation is computed over the time $t=0$ to $t=0.5$. The computational domain is divided into $E^{(x)}=5$ in spatial direction, with $E^{(x),l}=4$ in the liquid in the beginning and $E^{(x),l}=2$ in the liquid at the end. The interface velocity is computed by the normal probe method (50). The time step is chosen to be $\Delta t=3.84 \times 10^{-2}$
Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order $p^{(x)}$ in the ice phase. The comparison is made between the weak form method (46) and the normal probe method (50). The interface starts at $x=0.1$ and moves to $x=0.3$. The polynomial order in time direction, $p^{(t)}$, is chosen to be $10$ so that the temporal errors are negligible. The time step is chosen to be $\Delta t=4.9 \times 10^{-2}$. The Stefan number $St$ is set to be $0.02$
Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order $p^{(t)}$ in the ice phase. The interface velocity is computed by the weak form method (46). The interface starts at $x=0.1$ and moves to $x=0.3$. The comparison is made among different numbers of the space-time slab, i.e., $E^{(t)}=41$, $E^{(t)}=81$ and $E^{(t)}=162$. The corresponding time step is $\Delta t=4.9 \times 10^{-2}$, $\Delta t=2.5 \times 10^{-2}$ and $\Delta t=1.2 \times 10^{-2}$. The polynomial order in spatial direction, $p^{(x)}$, is chosen to be $10$ so that the spatial errors are negligible. The computational domain is divided into $E^{(x)}=5$ in the spatial direction. The Stefan number $St$ is set to be $0.02$
Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order $p^{(t)}$ in the ice phase. The interface velocity is computed by the normal probe method (50). The interface starts at $x=0.1$ and moves to $x=0.3$. The comparison is made among different numbers of the space-time slab, i.e., $E^{(t)}=41$, $E^{(t)}=81$ and $E^{(t)}=162$. The corresponding time step is $\Delta t=4.9 \times 10^{-2}$, $\Delta t=2.5 \times 10^{-2}$ and $\Delta t=1.2 \times 10^{-2}$. The polynomial order in spatial direction, $p^{(x)}$, is chosen to be $10$ so that the spatial errors are negligible. The computational domain is divided into $E^{(x)}=5$ in the spatial direction. The Stefan number $St$ is set to be $0.02$
Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order $p^{(x)}$ in the ice phase. The interface velocity is computed by the weak form method (46). The polynomial order in time direction, $p^{(t)}$, is chosen to be the same as $p^{(x)}$. The interface starts at $x=0.1$ and moves to $x=0.4$. The comparison is made among different values of $St$, i.e., $St=0.001$, $St=0.005$, $St=0.01$ and $St=0.05$. The computational domain is divided into $E^{(x)}=6$ in the spatial direction. The time step is chosen to be $\Delta t=0.9 \times 10^{-3}$
Errors in the temperature (left) and the interface velocity (right) as a function of polynomial order $p^{(x)}$ in the ice phase. The interface velocity is computed by the normal probe method (50). The polynomial order in time direction, $p^{(t)}$, is chosen to be the same as $p^{(x)}$. The interface starts at $x=0.1$ and moves to $x=0.4$. The comparison is made among different values of $St$, i.e., $St=0.001$, $St=0.02$, $St=0.05$ and $St=0.08$. The computational domain is divided into $E^{(x)}=6$ in the spatial direction. The time step is chosen to be $\Delta t=0.9 \times 10^{-3}$
Errors of the temperature ($\|Err_{\theta}\|_{\infty}$), interface position ($\|Err_{\Gamma}\|_{\infty}$) and interface velocity ($\|Err_{V_n}\|_{\infty}$) in the last space-time slab. The interface velocity is computed by applying the weak form method (46). The number of the Picard iterations is listed in the last column with $tol=10^{-13}$
 $p=(p^{(x)},p^{(t)})$ $\|Err_{\theta}\|_{\infty}$ $\|Err_{\Gamma}\|_{\infty}$ $\|Err_{V_n}\|_{\infty}$ iter (2, 1) 1.24E-003 1.75E-003 6.12E-003 11 (3, 1) 5.16E-005 6.40E-005 2.68E-004 11 (4, 1) 1.28E-006 1.47E-006 7.20E-006 9 (5, 1) 2.40E-008 2.63E-008 1.44E-007 8 (6, 1) 3.77E-010 4.02E-010 2.31E-009 5
 $p=(p^{(x)},p^{(t)})$ $\|Err_{\theta}\|_{\infty}$ $\|Err_{\Gamma}\|_{\infty}$ $\|Err_{V_n}\|_{\infty}$ iter (2, 1) 1.24E-003 1.75E-003 6.12E-003 11 (3, 1) 5.16E-005 6.40E-005 2.68E-004 11 (4, 1) 1.28E-006 1.47E-006 7.20E-006 9 (5, 1) 2.40E-008 2.63E-008 1.44E-007 8 (6, 1) 3.77E-010 4.02E-010 2.31E-009 5
Errors of the temperature ($\|Err_{\theta}\|_{\infty}$), interface position ($\|Err_{\Gamma}\|_{\infty}$) and interface velocity ($\|Err_{V_n}\|_{\infty}$) in the last space-time slab. The interface velocity is computed by applying the normal probe method (50). The number of the Picard iterations is listed in the last column with $tol=10^{-13}$
 $p=(p^{(x)},p^{(t)})$ $\|Err_{\theta}\|_{\infty}$ $\|Err_{\Gamma}\|_{\infty}$ $\|Err_{V_n}\|_{\infty}$ iter (2, 1) 3.89E-003 5.29E-003 1.34E-002 12 (3, 1) 1.55E-004 1.85E-004 5.08E-004 12 (4, 1) 3.73E-006 4.16E-006 1.25E-005 10 (5, 1) 6.80E-008 7.32E-008 2.29E-007 6 (6, 1) 9.97E-010 1.05E-009 3.47E-009 5
 $p=(p^{(x)},p^{(t)})$ $\|Err_{\theta}\|_{\infty}$ $\|Err_{\Gamma}\|_{\infty}$ $\|Err_{V_n}\|_{\infty}$ iter (2, 1) 3.89E-003 5.29E-003 1.34E-002 12 (3, 1) 1.55E-004 1.85E-004 5.08E-004 12 (4, 1) 3.73E-006 4.16E-006 1.25E-005 10 (5, 1) 6.80E-008 7.32E-008 2.29E-007 6 (6, 1) 9.97E-010 1.05E-009 3.47E-009 5
Order of convergence of the temperature ($\|Err_{\theta}\|_{\infty}$) for both options in the last space-time slab
 Order of convergence $p=(p^{(x)},p^{(t)})$ $E^{(x)}$ Weak form (46) Normal prob (50) (4, 4) 5 — — 10 3.58 3.79 20 4.04 3.85 (5, 5) 5 — — 10 4.65 4.82 20 5.38 4.89
 Order of convergence $p=(p^{(x)},p^{(t)})$ $E^{(x)}$ Weak form (46) Normal prob (50) (4, 4) 5 — — 10 3.58 3.79 20 4.04 3.85 (5, 5) 5 — — 10 4.65 4.82 20 5.38 4.89
Errors of the temperature ($\|Err_{\theta}\|_{\infty}$), interface position ($\|Err_{\Gamma}\|_{\infty}$) and interface velocity ($\|Err_{V_n}\|_{\infty}$) in the last space-time slab. The interface velocity is computed by applying the weak form method (46). The computational domain is divided into $E^{(x)}=6$ in the spatial direction. The interface starts at $x=0.1$ and moves to $x=0.3$. The time step is chosen to be $\Delta t=4.9 \times 10^{-2}$. The number of the Picard iteration is listed in the last column with $tol=10^{-15}$
 $p=(p^{(x)},p^{(t)})$ $\|Err_{\theta}\|_{\infty}$ $\|Err_{\Gamma}\|_{\infty}$ $\|Err_{V_n}\|_{\infty}$ iter (3, 3) 3.86E-006 8.44E-007 1.87E-007 7 (3, 10) 2.46E-006 7.67E-008 1.67E-008 7 (4, 4) 3.82E-009 9.46E-010 2.11E-010 7 (4, 10) 2.41E-009 1.74E-010 3.98E-011 7 (5, 5) 1.48E-009 4.72E-010 1.05E-010 7 (5, 10) 4.18E-010 1.28E-010 2.84E-011 7 (6, 6) 3.07E-012 9.52E-013 2.12E-013 7 (6, 10) 1.14E-012 3.47E-013 7.79E-014 7 (7, 7) 3.19E-013 9.91E-014 2.19E-014 7 (7, 10) 1.46E-013 4.55E-014 1.09E-014 7
 $p=(p^{(x)},p^{(t)})$ $\|Err_{\theta}\|_{\infty}$ $\|Err_{\Gamma}\|_{\infty}$ $\|Err_{V_n}\|_{\infty}$ iter (3, 3) 3.86E-006 8.44E-007 1.87E-007 7 (3, 10) 2.46E-006 7.67E-008 1.67E-008 7 (4, 4) 3.82E-009 9.46E-010 2.11E-010 7 (4, 10) 2.41E-009 1.74E-010 3.98E-011 7 (5, 5) 1.48E-009 4.72E-010 1.05E-010 7 (5, 10) 4.18E-010 1.28E-010 2.84E-011 7 (6, 6) 3.07E-012 9.52E-013 2.12E-013 7 (6, 10) 1.14E-012 3.47E-013 7.79E-014 7 (7, 7) 3.19E-013 9.91E-014 2.19E-014 7 (7, 10) 1.46E-013 4.55E-014 1.09E-014 7
Errors of the temperature ($\|Err_{\theta}\|_{\infty}$), interface position ($\|Err_{\Gamma}\|_{\infty}$) and interface velocity ($\|Err_{V_n}\|_{\infty}$) in the last space-time slab. The interface velocity is computed by applying the normal probe method (50). The computational domain is divided into $E^{(x)}=6$ in the spatial direction. The interface starts at $x=0.1$ and moves to $x=0.3$. The time step is chosen to be $\Delta t=4.9 \times 10^{-2}$. The number of the Picard iteration is listed in the last column with $tol=10^{-15}$
 $p=(p^{(x)},p^{(t)})$ $\|Err_{\theta}\|_{\infty}$ $\|Err_{\Gamma}\|_{\infty}$ $\|Err_{V_n}\|_{\infty}$ iter (3, 3) 1.36E-005 4.65E-006 2.34E-007 7 (3, 10) 1.29E-005 4.42E-006 1.84E-007 8 (4, 4) 8.27E-008 2.67E-008 1.26E-009 8 (4, 10) 8.50E-008 2.75E-008 1.43E-009 8 (5, 5) 1.44E-008 4.58E-009 1.02E-009 7 (5, 10) 1.42E-008 4.52E-009 1.01E-009 7 (6, 6) 1.43E-012 1.62E-014 4.03E-013 10 (6, 10) 1.20E-012 3.86E-013 4.81E-013 8 (7, 7) 4.08E-012 1.27E-012 2.80E-013 10 (7, 10) 4.20E-012 1.31E-012 2.87E-013 9
 $p=(p^{(x)},p^{(t)})$ $\|Err_{\theta}\|_{\infty}$ $\|Err_{\Gamma}\|_{\infty}$ $\|Err_{V_n}\|_{\infty}$ iter (3, 3) 1.36E-005 4.65E-006 2.34E-007 7 (3, 10) 1.29E-005 4.42E-006 1.84E-007 8 (4, 4) 8.27E-008 2.67E-008 1.26E-009 8 (4, 10) 8.50E-008 2.75E-008 1.43E-009 8 (5, 5) 1.44E-008 4.58E-009 1.02E-009 7 (5, 10) 1.42E-008 4.52E-009 1.01E-009 7 (6, 6) 1.43E-012 1.62E-014 4.03E-013 10 (6, 10) 1.20E-012 3.86E-013 4.81E-013 8 (7, 7) 4.08E-012 1.27E-012 2.80E-013 10 (7, 10) 4.20E-012 1.31E-012 2.87E-013 9
 [1] Zexuan Liu, Zhiyuan Sun, Jerry Zhijian Yang. A numerical study of superconvergence of the discontinuous Galerkin method by patch reconstruction. Electronic Research Archive, 2020, 28 (4) : 1487-1501. doi: 10.3934/era.2020078 [2] Yue Feng, Yujie Liu, Ruishu Wang, Shangyou Zhang. A conforming discontinuous Galerkin finite element method on rectangular partitions. Electronic Research Archive, , () : -. doi: 10.3934/era.2020120 [3] Leilei Wei, Yinnian He. A fully discrete local discontinuous Galerkin method with the generalized numerical flux to solve the tempered fractional reaction-diffusion equation. Discrete & Continuous Dynamical Systems - B, 2020  doi: 10.3934/dcdsb.2020319 [4] Shenglan Xie, Maoan Han, Peng Zhu. A posteriori error estimate of weak Galerkin fem for second order elliptic problem with mixed boundary condition. Discrete & Continuous Dynamical Systems - B, 2020  doi: 10.3934/dcdsb.2020340 [5] Mostafa Mbekhta. Representation and approximation of the polar factor of an operator on a Hilbert space. Discrete & Continuous Dynamical Systems - S, 2020  doi: 10.3934/dcdss.2020463 [6] Pierre-Etienne Druet. A theory of generalised solutions for ideal gas mixtures with Maxwell-Stefan diffusion. Discrete & Continuous Dynamical Systems - S, 2020  doi: 10.3934/dcdss.2020458 [7] Hirokazu Ninomiya. Entire solutions of the Allen–Cahn–Nagumo equation in a multi-dimensional space. Discrete & Continuous Dynamical Systems - A, 2021, 41 (1) : 395-412. doi: 10.3934/dcds.2020364 [8] Min Chen, Olivier Goubet, Shenghao Li. Mathematical analysis of bump to bucket problem. Communications on Pure & Applied Analysis, 2020, 19 (12) : 5567-5580. doi: 10.3934/cpaa.2020251 [9] Qingfang Wang, Hua Yang. Solutions of nonlocal problem with critical exponent. Communications on Pure & Applied Analysis, 2020, 19 (12) : 5591-5608. doi: 10.3934/cpaa.2020253 [10] Stefano Bianchini, Paolo Bonicatto. Forward untangling and applications to the uniqueness problem for the continuity equation. Discrete & Continuous Dynamical Systems - A, 2020  doi: 10.3934/dcds.2020384 [11] Lars Grüne, Matthias A. Müller, Christopher M. Kellett, Steven R. Weller. Strict dissipativity for discrete time discounted optimal control problems. Mathematical Control & Related Fields, 2020  doi: 10.3934/mcrf.2020046 [12] Serena Dipierro, Benedetta Pellacci, Enrico Valdinoci, Gianmaria Verzini. Time-fractional equations with reaction terms: Fundamental solutions and asymptotics. Discrete & Continuous Dynamical Systems - A, 2021, 41 (1) : 257-275. doi: 10.3934/dcds.2020137 [13] Marco Ghimenti, Anna Maria Micheletti. Compactness results for linearly perturbed Yamabe problem on manifolds with boundary. Discrete & Continuous Dynamical Systems - S, 2020  doi: 10.3934/dcdss.2020453 [14] Alberto Bressan, Sondre Tesdal Galtung. A 2-dimensional shape optimization problem for tree branches. Networks & Heterogeneous Media, 2020  doi: 10.3934/nhm.2020031 [15] Guido Cavallaro, Roberto Garra, Carlo Marchioro. Long time localization of modified surface quasi-geostrophic equations. Discrete & Continuous Dynamical Systems - B, 2020  doi: 10.3934/dcdsb.2020336 [16] Cuicui Li, Lin Zhou, Zhidong Teng, Buyu Wen. The threshold dynamics of a discrete-time echinococcosis transmission model. Discrete & Continuous Dynamical Systems - B, 2020  doi: 10.3934/dcdsb.2020339 [17] Awais Younus, Zoubia Dastgeer, Nudrat Ishaq, Abdul Ghaffar, Kottakkaran Sooppy Nisar, Devendra Kumar. On the observability of conformable linear time-invariant control systems. Discrete & Continuous Dynamical Systems - S, 2020  doi: 10.3934/dcdss.2020444 [18] Fioralba Cakoni, Pu-Zhao Kow, Jenn-Nan Wang. The interior transmission eigenvalue problem for elastic waves in media with obstacles. Inverse Problems & Imaging, , () : -. doi: 10.3934/ipi.2020075 [19] Marion Darbas, Jérémy Heleine, Stephanie Lohrengel. Numerical resolution by the quasi-reversibility method of a data completion problem for Maxwell's equations. Inverse Problems & Imaging, 2020, 14 (6) : 1107-1133. doi: 10.3934/ipi.2020056 [20] Hoang The Tuan. On the asymptotic behavior of solutions to time-fractional elliptic equations driven by a multiplicative white noise. Discrete & Continuous Dynamical Systems - B, 2020  doi: 10.3934/dcdsb.2020318

2019 Impact Factor: 1.27