A Discontinuous Galerkin–Finite Element Method for the Nonlinear Unsteady Burning Rate Responses of Solid Propellants

Abstract

The unsteady combustion of solid propellants under oscillating environments is the key to understanding the combustion instability inside solid rocket motors. The discontinuous Galerkin–finite element method (DG-FEM) is introduced to provide an efficient yet flexible numerical platform to investigate the combustion dynamics of solid propellants. The algorithm is developed for the classical unsteady model, the Zel’dovich–Novozhilov model. It is then validated based on a special analytical solution. The DG-FEM algorithm is then compared with the classical spectral method based on Laguerre polynomials. It is shown that the DG-FEM works more efficiently than the traditional spectral method, providing a more accurate solution with a lower computational cost.

1. Introduction

The combustion of solid propellants plays a central role in the performance of solid rocket motors (SRMs). In the combustion chamber, the chemical energy of propellants is converted into the mechanical energy of the work fluid. Since the invention of modern SRMs in the early 20th century, its vulnerability to combustion instability has been attracting the attention of scientists and engineers in related fields. This self-sustained instability is induced by the interaction of the combustion dynamics and the flow inside the chamber and may lead to large oscillations of both the flow field and the rocket structures [1].

Unlike liquid propellant rocket engines and gas turbines, the combustion process is confined to a thin region near the surface of the solid propellant. To fully understand the nature of this instability, it is crucial to investigate the transient combustion characteristics of the solid propellants. Due to the complex structures and chemical compositions of the oxidizers, binders and other additives, there are a series of complicated physical and chemical processes in the combustion of solid propellants [1,2]. Fortunately, research shows that some simplified models can capture the dynamic characteristics of the complex system well [3,4,5]. A widely adopted simplification is to assume that the propellants are homogeneous and the combustion is one-dimensional. This is acceptable as the characteristic temporal and spatial scales of the combustion in the solid rocket motor are determined by the acoustic modes of the chamber. In this large spatial and time scale, the heterogeneous structure and the non-uniform flame distribution can be omitted [1]. Additionally, the physical and chemical processes in the gas phase are assumed to be much faster than the heat conduction inside the condensed phase. Based on these assumptions, there are two categories of combustion models, where the condensed phases are treated in the same manner, while the gas phase is treated differently. In the category known as quasi-steady homogeneous one-dimensional (QSHOD) models, the gas phase is treated by flame models [6]. The other category is known as the Zel’dovich–Novozhilov(ZN) model, where the gas phase is treated as a black box via a phenomenological approach [2].

The two categories of models were compared in detail in a previous work [7,8]. It was shown that the two models behave differently in nonlinear regions despite their linear equivalency. Mathematically, the two models have the same form, consisting of exactly the same partial differential equation (PDE) governing the unsteady heat conduction in the solid phase, equipped with a system of boundary conditions derived from the gas processes. There are no analytical solutions to these equations, and the numerical solutions are not straightforward. The analytical treatments include linearization [9], and the nonlinear effects are included later by taking second-order terms into consideration [10]. However, it has been shown that the proper inclusion of the nonlinear effects is necessary for the correct modeling of the combustion instability of SRMs, especially the important triggering phenomenon [11]. In the modeling of thermoacoustic oscillations in SRMs, the unsteady combustion is modeled by the numerical solution [3,4,5].

The numerical treatment of the problem began with the finite difference method (FDM) [12]. The two major drawbacks of FDM are the complexity for higher-order accuracy and its implicit semi-discrete form. In FDM, the order of accuracy can only be improved by introducing higher-order finite difference forms, which leads to a complicated formulation. Additionally, the temperature field in the condensed phase is strongly coupled with the boundary conditions as the boundary is regression in an unsteady manner. This leads to further complexity when the time evolution of the algorithm is implicit. The problem is solved by the spectral method (SM) [13], where the Laguerre polynomials are used to approximate the exact solution. In the spectral method, the order of accuracy can be improved by simply including higher-order polynomials, and the algorithm has an explicit semi-discrete form automatically. However, there are also two major drawbacks of the spectral method. The first is that the retrieving of the temperature field is cumbersome as the calculation of the Laguerre polynomials requires additional effort, especially when the orders are relatively high. The second problem is the inability to extend it to higher dimensions with an irregular geometry. Thus, the present numerical algorithms cannot meet the requirements of the combustion instability modeling of SRMs.

Another algorithm, the discontinuous Galerkin–finite element method (DG-FEM), may provide a solution. The DG-FEM was first proposed for the neutron transportation problem in the 1970s [14]. The method is a combination of the finite volume method (FVM) and the finite element method (FEM), where the solution is approximated by polynomials inside each element as in FEM, while the inter-element interaction is treated by the flux function as in FVM. By choosing a higher-order polynomial inside elements, the high-order accuracy can be reached as in the traditional FEM. On the other hand, the discontinuous nature of DG-FEM breaks the band mass matrices in the traditional FEM. The mass matrices in DG-FEM are block diagonal, which can be inverted straightforwardly, yielding an explicit semi-discrete form. The DG-FEM can be used in various problems in different fields. There are general computational fluid dynamics platforms based on DG-FEM [15]. Additionally, its high fidelity makes it especially appropriate for the calculation of thermoacoustic oscillations [16]. For the unsteady combustion problem in the present paper, the algorithm should be in the explicit form and can achieve high-order accuracy easily.

In present paper, the DG-FEM is adopted to solve these problems for the first time. To derive a DG-FEM algorithm for the ZN model, the governing equation is converted to a system of first-order equations by introducing the heat flux as a new variable. The performance of the FDM, the SM, and the newly constructed algorithm is then compared. The first benchmark case is an exact solution of a modified ZN model, of which the governing equation is the same. The efficiencies of the three algorithms are compared in this case. The SM and the DG-FEM are further compared in the calculation of the nonlinear responses of the ZN model. The applicability and efficiency of the DG-FEM are investigated.

After the present introductory part, the mathematical formulations of the combustion model and the algorithm are discussed in Section 2. The performance of the existing algorithms is compared systematically in Section 3. Finally, the concluding remarks are given in Section 4.

2. The Formulation of the Algorithm

As discussed in a previous work [8], the unsteady behavior of the combustion models can be unified by a phenomenological approach. Thus, the phenomenological ZN model is used in the present paper as the benchmark.

2.1. The Physical Problem

The ZN model is a phenomenological approach without detailed modeling of the flame, following the notion of Zel’dovich [2,10,13]. The governing equation depicts the unsteady heat conduction inside the condensed phase:

$$\frac{\partial \theta }{\partial t}-R\frac{\partial \theta }{\partial x}-\frac{{\partial }^2\theta }{\partial x^2}=0,x\in [0,\infty ),t\in [0,\infty ),$$

(1)

where non-dimensional variables are used. The propellant is a semi-infinity set in $[0,\infty )$, and $\theta (x,t)$ is its non-dimensional temperature distribution. R is the burning rate and the regression rate of the burning surface, determined by the surface and gas phase conditions.

The gas phase process are described by the burning rate and the heat flux boundary condition:

$$\left\{\begin{array}{c}{\displaystyle R=P^{n}exp\left[k_{ZN}\left({\theta }_s-\frac{\phi }{R}\right)\right],R=exp\frac{k_{ZN}({\theta }_s-1)}{r_{ZN}}}\hfill \\ {\displaystyle \frac{\partial \theta (0,t)}{\partial x}=-\frac{R}{k_{ZN}}(nlnP+k_{ZN}{\theta }_s-lnR)=-\phi ,}\hfill \\ {\displaystyle \theta (+\infty ,t)=0}\hfill \end{array}\right.$$

(2)

where P is the ambient pressure, ${\theta }_s$ is the temperature of the burning surface, and the heat flux at the burning surface is $\phi$. There are three parameters $k_{ZN},r_{ZN},n$ in the model, which have been discussed in detail in [2,10]. n is the pressure component of the propellant. $k_{ZN},r_{ZN}$ are the sensitivity parameters of the burning rate and the burning surface temperature in the ambient environment, respectively.

The PDE (1) and the boundary condition (2) define a coupled problem of the unsteady temperature field and its regressing boundary. It should be remarked that the infinity boundary condition $\theta (+\infty ,t)=0$ can only be exactly fulfilled by the spectral method, as discussed in [7,8]. In the FDM or the DG-FEM in the present paper, an approximate Dirichlet boundary condition $\theta (x_R,t)=0$ is used, which has little influence on the key region as the temperature will decay exponentially away from the burning surface in the manner of $e^{-x}$.

2.2. The DG-FEM Algorithm

Based on the algorithm in [14], the DG-FEM algorithm can be constructed as follows. To begin with, the governing equation can be rewritten as a first-order system by introducing a new quantity q,

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial \theta }{\partial t}-\frac{\partial }{\partial x}(R\theta +q)=0,}\hfill \\ {\displaystyle q=\frac{\partial \theta }{\partial x},}\hfill \end{array}\right.$$

(3)

where q is the non-dimensional heat flux.

Suppose that the calculation domain $\mathsf{\Omega }=[x_L,x_R]$, where $x_L=0$, is approximated by a collection of elements ${\mathsf{D}}^k=[x_l^k,x_r^k]$, with $x_l^1=x_L,x_r^K=x_R$. The approximate solution is polynomials of order N,

$$\left[\begin{array}{c}{\displaystyle \theta (x,t)}\\ {\displaystyle q(x,t)}\end{array}\right]\simeq \left[\begin{array}{c}{\displaystyle {\theta }_h(x,t)}\\ {\displaystyle q_h(x,t)}\end{array}\right]=\underset{k=1}{\overset{K}{⨁}}\left[\begin{array}{c}{\displaystyle {\theta }_h^k(x,t)}\\ {\displaystyle q_h^k(x,t)}\end{array}\right]=\underset{k=1}{\overset{K}{⨁}}\sum _{i=1}^{N_p}\left[\begin{array}{c}{\displaystyle {\theta }_h^k(x_i,t)}\\ {\displaystyle q_h^k(x_i,t)}\end{array}\right]{\ell }_i^k\left(x\right),$$

(4)

where $N_p=N+1$ is the number of terms of a polynomial with order N. The subscript “h” represents the approximation, which is the direct sum of K local representations ${\theta }_h^k(x,t),q_h^k(x,t)$. For each ${\theta }_h^k(x,t),q_h^k(x,t)$, it is further expanded by $N_p$ Lagrangian polynomials ${\ell }_i^k\left(x\right)$, where $x_i$ are the interpolation points. This representation of the approximation solution is thus named the nodal form. It can be rewritten as the modal form

$$\left[\begin{array}{c}{\displaystyle {\theta }_h^k(x,t)}\\ {\displaystyle q_h^k(x,t)}\end{array}\right]=\sum _{n=1}^{N_p}\left[\begin{array}{c}{\displaystyle {\widehat{\theta }}_n^k\left(t\right)}\\ {\displaystyle {\widehat{q}}_n^k\left(t\right)}\end{array}\right]{\psi }_n\left(x\right),$$

where the spatial and temporal dependences are separated by the modal functions ${\psi }_n\left(x\right)$ and their coefficients.

We introduce a globally defined space ${\mathsf{V}}_h$ of test functions, ${\varphi }_h$, as ${\mathsf{V}}_h={⨁}_{k=1}^K{\mathsf{V}}_h^k$, where the locally defined spaces are ${\mathsf{V}}_h^k=\mathrm{span}{\left\{{\psi }_n\left({\mathsf{D}}^k\right)\right\}}_{n=1}^{N_p}$. A locally defined member ${\varphi }_h^k$ can be expressed as

$${\varphi }_h^k\left(x\right)=\sum _{n=1}^{N_p}{\widehat{\varphi }}_n^k{\psi }_n\left(x\right).$$

The approximation of the governing equation is assumed to fulfill

$$\left\{\begin{array}{c}{\displaystyle {\int }_{{\mathsf{D}}^k}\left[\frac{\partial {\theta }_h}{\partial t}-\frac{\partial }{\partial x}(R\left(t\right){\theta }_h+q_h)\right]{\psi }_n\left(x\right)\mathrm{d}x=0,1\le n\le N_p,}\hfill \\ {\displaystyle {\int }_{{\mathsf{D}}^k}\left[q_h-\frac{\partial {\theta }_h}{\partial x}\right]{\psi }_n\left(x\right)\mathrm{d}x=0,1\le n\le N_p,}\hfill \end{array}\right.$$

or equivalently

$$\left\{\begin{array}{c}{\displaystyle {\int }_{x_l^k}^{x_r^k}\left[\frac{\partial {\theta }_h}{\partial t}-\frac{\partial }{\partial x}(R\left(t\right){\theta }_h+q_h)\right]{\psi }_n\left(x\right)\mathrm{d}x=0,1\le n\le N_p,}\hfill \\ {\displaystyle {\int }_{x_l^k}^{x_r^k}\left[q_h-\frac{\partial {\theta }_h}{\partial x}\right]{\psi }_n\left(x\right)\mathrm{d}x=0,1\le n\le N_p.}\hfill \end{array}\right.$$

Spatial integration by parts leads to the semi-discrete form

$$\left\{\begin{array}{c}{\displaystyle {\int }_{{\mathsf{D}}^k}\left[\frac{\partial {\theta }_h^k}{\partial t}{\psi }_n+(R\left(t\right){\theta }_h^k+q_h^k)\frac{\mathrm{d}{\psi }_n}{\mathrm{d}x}\right]\mathrm{d}x={\int }_{\partial {\mathsf{D}}^k}\widehat{\mathit{n}}·{(R\left(t\right){\theta }_h^k+q_h^k)}^{*}{\psi }_n\mathrm{d}x,1\le n\le N_p,}\hfill \\ {\displaystyle {\int }_{{\mathsf{D}}^k}\left[q_h{\psi }_n+{\theta }_h\frac{\mathrm{d}{\psi }_n}{\mathrm{d}x}\right]\mathrm{d}x={\int }_{\partial {\mathsf{D}}^k}\widehat{\mathit{n}}·{\left({\theta }_h^k\right)}^{*}{\psi }_n\mathrm{d}x,1\le n\le N_p,}\hfill \end{array}\right.$$

(5)

or equivalently

$$\left\{\begin{array}{c}{\displaystyle {\int }_{x_l^k}^{x_r^k}\left[\frac{\partial {\theta }_h^k}{\partial t}{\psi }_n+(R\left(t\right){\theta }_h^k+q_h^k)\frac{\mathrm{d}{\psi }_n}{\mathrm{d}x}\right]\mathrm{d}x={(R\left(t\right){\theta }_h^k+q_h^k)}^{*}{\psi }_n{|}_{x_l^k}^{x_r^k},1\le n\le N_p,}\hfill \\ {\displaystyle {\int }_{x_l^k}^{x_r^k}\left[q_h{\psi }_n+{\theta }_h\frac{\mathrm{d}{\psi }_n}{\mathrm{d}x}\right]\mathrm{d}x={\left({\theta }_h^k\right)}^{*}{\psi }_n{|}_{x_l^k}^{x_r^k},1\le n\le N_p,}\hfill \end{array}\right.$$

which is the weak form of the original equation. The terms with superscript stars denote the flux on the boundary of two elements. It is a central point of the DG-FEM, which is the “FVM-like” part of the algorithm. In the present problem, the central flux is adopted for the temperature ${\theta }_h$ and the upwind flux is selected for $(R\left(t\right){\theta }_h-q_h)$,

$$\begin{array}{c}{\displaystyle {\left({\theta }_h\right)}^{*}=\frac{{\left({\theta }_h\right)}^{+}+{\left({\theta }_h\right)}^{-}}{2},}\hfill \\ {\displaystyle {(R\left(t\right){\theta }_h+q_h)}^{*}}\hfill \\ {\displaystyle =\frac{{(R{\theta }_h+q_h)}^{+}+{(R{\theta }_h+q_h)}^{-}}{2}+\left|R\right|\left({\widehat{\mathit{n}}}^{+}{(R{\theta }_h+q_h)}^{+}+{\widehat{\mathit{n}}}^{-}{(R{\theta }_h+q_h)}^{-}\right)}\hfill \\ {\displaystyle =\frac{{(R{\theta }_h+q_h)}^{+}+{(R{\theta }_h+q_h)}^{-}}{2}+\left|R\right|\left({(R{\theta }_h+q_h)}^{+}-{(R{\theta }_h+q_h)}^{-}\right),}\hfill \end{array}$$

(6)

where the superscripts “+” and “−” denote the values on the right and the left side of the boundary, respectively.

Next, by mapping the elements ${\mathsf{D}}^k$ to a standard reference interval $\mathsf{I}=[-1,1]$, the modal functions and nodal functions can be treated in the same manner in different elements. The coefficients of the derivatives constitute a matrix, the $(i,j)$-term of which is

$${\mathcal{M}}_{ij}^k={\int }_{x_l^k}^{x_r^k}{\ell }_i^k\left(x\right){\ell }_j^k\left(x\right)\mathrm{d}x=\frac{h^k}{2}{\int }_{-1}^1{\ell }_i\left(r\right){\ell }_j\left(r\right)\mathrm{d}r=\frac{h^k}{2}{({\ell }_i,{\ell }_j)}_{\mathsf{I}}=\frac{h^k}{2}{\mathcal{M}}_{ij}.$$

(7)

It can be seen that the matrices for different elements have the same form, except for a length scale $h^k/2$. The semi-discrete form can then be converted to an explicit one by multiplying the inverse mass matrix on both sides of the equations. Another class of coefficients comprises the spatial derivatives of modal functions:

$${\mathcal{S}}_{ij}^k={\int }_{x_l^k}^{x_r^k}{\ell }_i^k\left(x\right)\frac{\mathrm{d}{\ell }_j^k\left(x\right)}{\mathrm{d}x}\mathrm{d}x={\int }_{-1}^1{\ell }_i\left(x\right)\frac{\mathrm{d}{\ell }_j\left(x\right)}{\mathrm{d}r}\mathrm{d}r={\mathcal{S}}_{ij},$$

(8)

which constitute the stiffness matrices. Substituting the mass and stiffness matrices into the weak form of Equation (5) leads to

$$\left\{\begin{array}{c}{\displaystyle \frac{\mathrm{d}{\theta }_h^k}{\mathrm{d}t}=-{\left({\mathcal{M}}^k\right)}^{-1}{\left({\mathcal{S}}^k\right)}^{T}R\left(t\right)({\theta }_h^k+q_h^k)+{\left({\mathcal{M}}^k\right)}^{-1}{(R\left(t\right){\theta }_h^k+q_h^k)}^{*}{\ell }^k\left(x\right){|}_{x_l^k}^{x_r^k},1\le n\le N_p,}\hfill \\ {\displaystyle q_h=-{\left({\mathcal{M}}^k\right)}^{-1}{\left({\mathcal{S}}^k\right)}^T{\theta }_h+{\left({\mathcal{M}}^k\right)}^{-1}{\left({\theta }_h^k\right)}^{*}{\ell }^k\left(x\right){|}_{x_l^k}^{x_r^k},1\le n\le N_p.}\hfill \end{array}\right.$$

(9)

The semi-discrete weak form (9), together with the flux functions (6) defined, is the semi-discrete form of the original Equation (1), where the spatial derivatives are discretized and the temporal derivatives are kept.

Finally, a fourth-order Runge–Kutta method is used for the temporal discretization following [14]. The overall algorithm will thus provide spatial accuracy of the N-th order and temporal accuracy of the fourth order. The determination of the size of the time step follows

$$\Delta t=C_t\underset{i,k}{min}{(\Delta x_i^k)}^2,$$

(10)

where $\Delta x_i^k$ is the spatial mesh size of the i-th nodes in the k-th element, and $C_t$ is the Courant–Friedrichs–Leray-like number. In the present work, $C_t$ is set to 0.5 to guarantee the robustness of the algorithm, according to [14]. It should be remarked that for higher efficiency, the time step size should be determined adaptively based on the physical features in the time evolution [17,18]. In the application of the ZN model in the combustion instability study, the temporal characteristics are governed by the acoustic modes of the chamber, and a fixed time step is adequate. However, if the algorithm is further used for more complicated models of solid propellants where the microscopic structures are resolved, an adaptive time step algorithm will be needed to maximize the calculation efficiency.

3. Results and Discussion

The performance of the DG-FEM algorithm will be compared with that of the two traditional algorithms in two typical problems. The FDM treatment, using the backward time central space (BTCS) algorithm, is proposed in [12], which has been used in the combustion instability research in solid rocket motors [3]. The spectral method treatment, using the Laguerre polynomial expansion, is proposed by [13], which has also been used in the combustion instability research in solid rocket motors [4,5].

3.1. An Exact Solution

The ZN model does not yield analytical solutions, but we can define an analytical solution directly with a modified boundary condition as

$$\begin{array}{c}\theta (x,t)=\exp (-x+\sin (t\left)\right)\hfill \\ {\displaystyle \frac{\partial \theta (0,t)}{\partial x}=-\exp \left(\sin t\right),R\left(t\right)=1-\cos t}\hfill \end{array}$$

(11)

To evaluate the errors of each algorithm, the $L^2$-norm of the deviations to the exact surface temperature ${\theta }_s\left(t\right)=\theta (0,t)$ for the calculated time interval is used.

$$ϵ={\left({\int }_0^T{\parallel {\theta }_h(0,t)-\theta (0,t)\parallel }^2\mathrm{d}t\right)}^{1/2}$$

The performance of the two traditional treatments, the FDM with the BTCS algorithm and the SM with Laguerre polynomial expansion, and the newly proposed DG-FEM is compared. The calculation is carried out on an Intel i7-12700h platform.

3.1.1. The FDM

The algorithm used here is based on [3,12]. A spatial coordinate transformation $\eta =e^{-kx}$ is used to improve the resolution near the burning surface with a high temperature gradient [3]. The parameter k controls the density of the meshes near the burning surface, but a larger k will also lead to the poorer robustness of the algorithm. In the present work, $k=0.5$ is adopted, and the transformation also maps the interval $(0,+\infty )$ to $[0,1]$.

$$\begin{array}{c}{\displaystyle \frac{\partial \theta }{\partial t}+(Rk\eta -k^2\eta )\frac{\partial \theta }{\partial \eta }-{\left(k\eta \right)}^2\frac{{\partial }^2\theta }{\partial {\eta }^2}=0,x\in [0,1],t\in [0,\infty ),}\hfill \\ {\displaystyle \frac{\partial \theta (1,t)}{\partial \eta }=\frac{1}{k}exp(sint),}\hfill \\ {\displaystyle \theta (0,t)=0.}\hfill \end{array}$$

(12)

The BTCS scheme for the inner grids is

$$\frac{{\theta }_i^{(n+1)}-{\theta }_i^{\left(n\right)}}{\Delta t}+C_1\frac{{\theta }_{i+1}^{(n+1)}-{\theta }_{i-1}^{\left(n\right)}}{2\Delta \eta }-C_2\frac{{\theta }_{i+1}^{(n+1)}-2{\theta }_i^{(n+1)}+{\theta }_{i-1}^{(n+1)}}{{(\Delta \eta )}^2}=0,i=2,\dots ,N_{\eta }-1,$$

where $\Delta t,\Delta \eta$ are the step sizes for temporal and spatial variables; the subscript and superscript denote the number of nodes in the spatial and temporal directions, and $C_1=R\left(n\right)k{\eta }_i-k^2{\eta }_i,C_2={\left(k{\eta }_i\right)}^2$.

The treatment for the left boundary is trivial; it is the Dirichlet zero boundary condition and is already included in the equations for the inner grids. The right boundary is a Neumann boundary condition. By introducing and extrapolating nodes outside the boundary ${\theta }_{N_{\eta +1}}$, the last iteration equation for the right boundary can be solved as

$$\frac{2C_2}{0.5C_1-C_2}{\theta }_{N_{\eta }-1}^{(n+1)}-\frac{1+2C_2}{0.5C_1-C_2}{\theta }_{N_{\eta }}^{(n+1)}=2\Delta \eta \frac{\partial \theta (1,t)}{\partial \eta }-\frac{1}{0.5C_1-C_2}{\theta }_{N_{\eta }}^{\left(n\right)}$$

Thus, the original PDE is converted to a system of $(N_{\eta }-1)$-finite difference equations. In the system of equations, the unknowns, which are the temperatures in the $(n+1)$-time step, cannot be directly represented by the already known temperatures in the $\left(n\right)$-time step. This is the implicit nature of the BTCS scheme. As a system of linear equations with a sparse matrix (a tridiagonal matrix), a numerical solution by iteration performs better. The successive over-relaxation (SOR) iteration method can then be used to solve the system of equations. By requiring the temperature field variation between two successive iterations to be controlled by a predefined value, convergence in one time step is reached.

Five different mesh numbers $N_{\eta }=10,50,100,150,200$ are selected for calculation, with corresponding time steps $\Delta t=C(\Delta \eta ),C=0.2$, following [3]. It should be remarked that the calculation of the problem (11) by the BTCS FDM is much faster than the real unsteady combustion of solid propellants. This is because the burning rate $R\left(t\right)$ is given explicitly in each time step by Equation (11), reducing the problem to a linear one, while a loop for burning rate iteration is needed in the ZN model (1) or other unsteady combustion models.

The error for the burning surface temperature is used for the evaluation of the algorithm. For the non-dimensional time interval $[0,1000]$, the results for calculation are summarized in Figure 1.

It can be seen from the figure that the calculation errors are periodic in phase with the given oscillations. The errors decrease with the increasing number of meshes, as expected. However, the calculation costs improve rapidly as the number of time steps increases, as shown in

Table 1. The calculation errors and costs of the FDM.

Number of Meshes ${\mathit{N}}_{\mathit{\eta }}$	${\mathit{L}}^2$-Error $\mathit{ϵ}$	Normalized Calculation Time
10	2.582	1
50	0.4218	54.42
100	0.2006	344.52
150	0.1325	1076.08
200	0.0989	2454.0

, where the calculation time has been normalized based on the time for the first case, $2.4158$ s.

It should be remarked that the temporal and spatial steps are actually correlated by the relation $\Delta t=C(\Delta \eta )$. Once the two important parameters k and C are fixed, the time step size is fixed for each $N_{\eta }$. There may exist further potential for higher efficiency with varied $k,C$, or, equivalently, more appropriate temporal and spatial step size combinations. However, this will lead to the lower robustness of the algorithm. It has been found that divergences are more frequently encountered at a slightly larger time step size.

3.1.2. The SM

The algorithm used here is based on [7,8], where the temperature is expanded by the Laguerre polynomials $L_i\left(x\right)$ as

$$\theta (t,x)\sim e^{-x}\sum _{i=0}^{N_{\mathrm{sp}}}{u}_i\left(t\right)L_i\left(x\right),L_i\left(x\right)=\frac{e^x}{i!}\frac{d^i}{d{x}^i}{e}^{-x}{x}^i,$$

where $N_{\mathrm{sp}}$ is the order of the spectral expansion. By substituting the above expansion into the governing Equation (1) and using the properties of the Laguerre polynomials, the PDE can be semi-discretized as

$$\left\{\begin{array}{c}{\displaystyle {\dot{u}}_0=\phi -R{\theta }_s,}\hfill \\ {\displaystyle {\dot{u}}_1=\phi -{\theta }_s+R(u_0-{\theta }_s),}\hfill \\ {\displaystyle {\dot{u}}_i=\phi -i{\theta }_s+\sum _{j=0}^{i-2}(i-j-1)u_j}\hfill \\ {\displaystyle +R\left(\sum _{j=0}^{i-1}{u}_j-{\theta }_s\right),i=2,3,\dots ,N,}\hfill \end{array}\right.$$

(13)

The calculation task is the same as for the FDM. There is a special characteristic for the SM where there is no spatial mesh, and there is no general dependence of the time step size on the number of polynomials. As a result, the time step size $\Delta t$ is set as a calculation parameter along with $N_{\mathrm{sp}}$. The calculation errors behave similarly as in the FDM, and the error plot is omitted here. The performance of the SM is shown in

Table 2. The calculation errors and costs of the SM.

Order of Polynomials ${\mathit{N}}_{\mathbf{sp}}$	Time Step Size $\mathit{\Delta }\mathit{t}$	${\mathit{L}}^2$-Error $\mathit{ϵ}$	Normalized Calculation Time
5	${10}^{-2}$	$7.31\times {10}^{-3}$	$0.3132$
10	${10}^{-2}$	$3.55\times {10}^{-3}$	$0.3253$
15	${10}^{-2}$	$3.5\times {10}^{-3}$	$0.3425$
15	${10}^{-3}$	$3.49\times {10}^{-4}$	$4.363$
20	${10}^{-3}$	$3.48\times {10}^{-4}$	$4.723$

, where the calculation time has been normalized based on the time for the first case in FDM, $2.4158$ s. It can be seen from the table that the SM yields much more accurate results than the FDM, with a considerably shorter calculation time. Additionally, it also can be seen from the results that the performance of the SM is highly influenced by the time step size, especially for $N_{\mathrm{sp}}>10$. From $N_{\mathrm{sp}}=5$ to $N_{\mathrm{sp}}=10$, the calculation error can be reduced by 50% with only a slightly longer calculation time. However, for a larger $N_{\mathrm{sp}}$, the error becomes almost invariant with increasing $N_{\mathrm{sp}}$. which is in accordance with the results demonstrated in the eigenvalue calculation in [7].

3.1.3. The DG-FEM

The DG-FEM for solid propellants proposed in present work has two calculation parameters, the number of elements K and the order of polynomials $N_{\mathrm{DG}}$. The time step is determined by the minimum sub-mesh in each element according to (10). A similar exponentially growing mesh as for the FDM case is generated with a dense distribution near the boundary. The calculation results are collected in

Table 3. The calculation errors and costs of the DG-FEM.

Order of Polynomials ${\mathit{N}}_{\mathbf{DG}}$	Number of Elements K	${\mathit{L}}^2$-Error $\mathit{ϵ}$	Normalized Calculation Time
2	10	$0.9419$	$0.062$
4	10	$0.0204$	$0.4793$
6	10	$1.22\times {10}^{-3}$	$2.016$
8	10	$1.46\times {10}^{-4}$	$6.002$
10	10	$2.79\times {10}^{-5}$	$14.14$
10	5	$2.81\times {10}^{-5}$	$12.86$
10	20	$2.78\times {10}^{-5}$	$16.24$

, where the calculation time has been normalized based on the time for the first case in FDM, $2.4158$ s. It can be seen that the DG-FEM performs slightly better than the SM. Another important feature is that the calculation errors are mainly determined by $N_{\mathrm{DG}}$. In the $N_{\mathrm{DG}}=10$ case, the three values $K=5,10,20$ yield similar accuracies with similar time.

It is noticeable from the calculation that the influence of the number of elements K is much smaller than the order of polynomials, while, in the FDM, the performance is quite different for different meshes. The reason for the difference is that the “effective” mesh size is actually the inner mesh inside each element determined by the interpolation polynomials. For example, the time step size is determined by (10), where the minimum mesh is much smaller than the size of an element. Thus, for the same $N_{\mathrm{DG}}$, the sizes of the minimum mesh in different Ks may be close. Additionally, the solution of the governing equation decays exponentially in the spatial direction, and the dynamics are determined by the grids near the boundary. This also contributes to the weak dependence on the number of elements K for the algorithm.

3.1.4. The Overall Comparison

Summarizing the above calculations leads to a direct comparison of the three methods in terms of their performance on the simplified problem with an exact solution, as shown in Figure 2. The results of the three algorithms generally fit well on lines in the double-log diagram, which implies the linear dependence of the order of the calculation error on the order of the calculation time. The slope (in absolute value) indicates the efficiency of the algorithm. A steeper line implies a faster reduction in the error with the increasing calculation cost.

It can be seen from the figure that the FDM is far less efficient than the other two algorithm. Confined by the low order of accuracy, the error of the FDM can only be reduced by decreasing the size of the meshes, which is computationally expensive. The slope is smaller (in the absolute value) than the other two lines. The efficiency of the DG-FEM and the SM is generally of the same level, of which the SM is slightly better. Indeed, the major advantage of the SM is its ability to gain a high order of accuracy. However, as a more flexible algorithm, the performance of the DG-FEM is remarkable, reaching similar accuracy as the SM with a similar computation time.

3.2. The Capture of Complex Dynamics

It has been found that the ZN model yields chaotic behavior in some specific parameter ranges [2,8]. The chaotic motion, determined by some strange attractor, is known to be an extremely subtle structure. In the previous subsection, it is shown that the FDM is far less efficient than the SM and the DG-FEM, and the latter two are compared in terms of the complicated nonlinear behavior of the ZN model.

The PDE (1) equipped with boundary conditions (2) is investigated with the parameters $k_{ZN}=2.09,r_{ZN}=0.33$, at which the system will experience an intrinsic chaotic combustion oscillation.

3.2.1. The SM

When using the SM with $N_{\mathrm{sp}}=15,\Delta t=1\times {10}^{-3}$, the calculation time is 89.52 s. We can see in Figure 3 that the burning rate experiences an irregular oscillation in the truncated $[500,1000]$ time interval. The irregularity can be further confirmed by its fast Fourier transform (FFT), where the frequency spectrum is distributed continuously, except for the eigenfrequency and its higher-order resonances.

A further investigation is carried out via the phase space reconstruction in Figure 4. The phase space reconstruction is carried out by extracting vectors from the single time series $R\left(t\right)$. A common reconstruction method is through the time delay [19], where the three time series are $R\left(t\right),R(t-\Delta \tau ),R(t-2\Delta \tau ),\Delta \tau =0.1$.

A typical strange attractor can be seen in the reconstructed phase portrait in Figure 4a. The Poincaré map is a cross-section view of the attractor, where the intersected points become a line on the Poincaré section. This also indicates the chaotic behavior of the solution.

To investigate the resolution limit of the SM algorithm, the major influence parameter $\Delta t$ is increased gradually. When $\Delta t$ grows to $1\times {10}^{-2}$, the calculation time reduces to 7.37 s, and the burning rate time series enters a different state, as shown in Figure 5. The burning rate oscillates in a much more regular manner, and the frequency distribution becomes separated.

The differences in the phase space are clearer, as shown in Figure 6. The trajectory becomes less dense than in Figure 4, which indicates that the oscillation becomes periodic. This is also validated by the Poincaré section. Although there are a number of intersected points, they are distributed separately in the section.

It should be remarked that the above case is a critical one, where the solution is on the edge of chaotic and periodic motions. Thus, $\Delta t=1\times {10}^{-2}$ can be seen as the critical value for the SM to accurately capture the complicated chaotic dynamics of the ZN model. Indeed, it can be verified that a smaller $\Delta t$ yields chaotic motion and a larger $\Delta t$ may lead to the divergence of the algorithm. At the critical time step size, the calculation time is 7.37 s.

3.2.2. The DG-FEM

It is not surprising that the DG-FEM with sufficient accuracy also yields a correct chaotic burning rate response as the SM does. The calculation results, as well as the phase reconstruction treatments, are collected in Figure 7. The results are almost the same as the results for the SM with $\Delta t=1\times {10}^{-3}$. Recall that the major calculation parameter of the DG-FEM is the order of the polynomials, and an $N_{\mathrm{DG}}=8,K=10$ setting yields the results in 18.4 s.

By gradually reducing $N_{\mathrm{DG}}$, we find that the solution becomes inaccurate when $N_{\mathrm{DG}}=5$, as shown in Figure 8. A typical feature of this solution is that there exists a series of small-amplitude oscillations near the unity between large-amplitude oscillations, although the overall behavior still seems chaotic. The small-amplitude part will grow with smaller $N_{\mathrm{DG}}$, and the oscillation will disappear at $N_{\mathrm{DG}}=2$.

The calculation time for $N_{\mathrm{DG}}=5$ case is 3.2 s, which is remarkable compared with the SM with $\Delta t=1\times {10}^{-2}$. The DG-FEM uses a much smaller calculation time to reach slightly better results than the SM.

4. Conclusions

A high-accuracy DG-FEM algorithm is constructed for the unsteady burning rate response model for solid propellants. It is then shown that the DG-FEM algorithm performs better than the traditionally used FDM and SM algorithms.

Based on a modified problem of the ZN model with an exact solution, the performance of the three algorithms is compared thoroughly. The mesh size, the time step size, and the order of the interpolation polynomials of the FDM, the SM, and the DG-FEM are found to be the key parameters for the three algorithms, respectively. They all yield higher accuracy with a longer calculation time. From the slope in the double-log diagram of the accuracy versus calculation time, it is shown that the SM and DG-FEM are far better than the FDM, and the DG-FEM is the best.

A further comparison is carried out for the intrinsic instability of the ZN model in the chaotic regime. It is shown that at the “inaccurate limit”, where the algorithm is no longer able to capture the subtle chaotic feature, the DG-FEM requires far less calculation time than the SM.

The newly constructed DG-FEM algorithm is more suitable for unsteady burning rate response calculations than the two existing algorithms. Additionally, the DG-FEM also has an important advantage in that it can be easily extended to higher-dimensional structures with complex geometry.