# 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  November 2018 Early access  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] Yanzhao Cao, Li Yin. Spectral Galerkin method for stochastic wave equations driven by space-time white noise. Communications on Pure & Applied Analysis, 2007, 6 (3) : 607-617. doi: 10.3934/cpaa.2007.6.607 [2] Yinhua Xia, Yan Xu, Chi-Wang Shu. Efficient time discretization for local discontinuous Galerkin methods. Discrete & Continuous Dynamical Systems - B, 2007, 8 (3) : 677-693. doi: 10.3934/dcdsb.2007.8.677 [3] Yuming Zhang. On continuity equations in space-time domains. Discrete & Continuous Dynamical Systems, 2018, 38 (10) : 4837-4873. doi: 10.3934/dcds.2018212 [4] Marita Holtmannspötter, Arnd Rösch, Boris Vexler. A priori error estimates for the space-time finite element discretization of an optimal control problem governed by a coupled linear PDE-ODE system. Mathematical Control & Related Fields, 2021, 11 (3) : 601-624. doi: 10.3934/mcrf.2021014 [5] Zeyu Xia, Xiaofeng Yang. A second order accuracy in time, Fourier pseudo-spectral numerical scheme for "Good" Boussinesq equation. Discrete & Continuous Dynamical Systems - B, 2020, 25 (9) : 3749-3763. doi: 10.3934/dcdsb.2020089 [6] Vincent Astier, Thomas Unger. Galois extensions, positive involutions and an application to unitary space-time coding. Advances in Mathematics of Communications, 2019, 13 (3) : 513-516. doi: 10.3934/amc.2019032 [7] Dong-Ho Tsai, Chia-Hsing Nien. On space-time periodic solutions of the one-dimensional heat equation. Discrete & Continuous Dynamical Systems, 2020, 40 (6) : 3997-4017. doi: 10.3934/dcds.2020037 [8] Susanne Pumplün, Thomas Unger. Space-time block codes from nonassociative division algebras. Advances in Mathematics of Communications, 2011, 5 (3) : 449-471. doi: 10.3934/amc.2011.5.449 [9] Gerard A. Maugin, Martine Rousseau. Prolegomena to studies on dynamic materials and their space-time homogenization. Discrete & Continuous Dynamical Systems - S, 2013, 6 (6) : 1599-1608. doi: 10.3934/dcdss.2013.6.1599 [10] Dmitry Turaev, Sergey Zelik. Analytical proof of space-time chaos in Ginzburg-Landau equations. Discrete & Continuous Dynamical Systems, 2010, 28 (4) : 1713-1751. doi: 10.3934/dcds.2010.28.1713 [11] Frédérique Oggier, B. A. Sethuraman. Quotients of orders in cyclic algebras and space-time codes. Advances in Mathematics of Communications, 2013, 7 (4) : 441-461. doi: 10.3934/amc.2013.7.441 [12] Grégory Berhuy. Algebraic space-time codes based on division algebras with a unitary involution. Advances in Mathematics of Communications, 2014, 8 (2) : 167-189. doi: 10.3934/amc.2014.8.167 [13] David Grant, Mahesh K. Varanasi. Duality theory for space-time codes over finite fields. Advances in Mathematics of Communications, 2008, 2 (1) : 35-54. doi: 10.3934/amc.2008.2.35 [14] Montgomery Taylor. The diffusion phenomenon for damped wave equations with space-time dependent coefficients. Discrete & Continuous Dynamical Systems, 2018, 38 (11) : 5921-5941. doi: 10.3934/dcds.2018257 [15] Na An, Chaobao Huang, Xijun Yu. Error analysis of discontinuous Galerkin method for the time fractional KdV equation with weak singularity solution. Discrete & Continuous Dynamical Systems - B, 2020, 25 (1) : 321-334. doi: 10.3934/dcdsb.2019185 [16] Konstantinos Chrysafinos, Efthimios N. Karatzas. Symmetric error estimates for discontinuous Galerkin approximations for an optimal control problem associated to semilinear parabolic PDE's. Discrete & Continuous Dynamical Systems - B, 2012, 17 (5) : 1473-1506. doi: 10.3934/dcdsb.2012.17.1473 [17] Qingjie Hu, Zhihao Ge, Yinnian He. Discontinuous Galerkin method for the Helmholtz transmission problem in two-level homogeneous media. Discrete & Continuous Dynamical Systems - B, 2020, 25 (8) : 2923-2948. doi: 10.3934/dcdsb.2020046 [18] Xiaomao Deng, Xiao-Chuan Cai, Jun Zou. A parallel space-time domain decomposition method for unsteady source inversion problems. Inverse Problems & Imaging, 2015, 9 (4) : 1069-1091. doi: 10.3934/ipi.2015.9.1069 [19] David Grant, Mahesh K. Varanasi. The equivalence of space-time codes and codes defined over finite fields and Galois rings. Advances in Mathematics of Communications, 2008, 2 (2) : 131-145. doi: 10.3934/amc.2008.2.131 [20] Xiaopeng Zhao. Space-time decay estimates of solutions to liquid crystal system in $\mathbb{R}^3$. Communications on Pure & Applied Analysis, 2019, 18 (1) : 1-13. doi: 10.3934/cpaa.2019001

2020 Impact Factor: 1.327