A Second-order Modified Ghost Fluid Method (2nd-MGFM) with Discontinuous Galerkin Method for 1-D compressible Multi-medium Problem with Cylindrical and Spherical Symmetry

Abstract

In this work, we develop the discontinuous Galerkin method to simulate 1-D cylindrical and spherical compressible multi-medium flows with an immiscible interface. To treat the interface with higher-order accuracy, the modified ghost fluid method is extended to a second-order version with source terms, in which linearly distributed ghost fluid states are constructed. A multi-medium generalized Riemann problem with the geometrical source is constructed to predict the states and the spatial derivatives at the interface. The predicted interface states and spatial derivatives are then employed to define the linearly distributed ghost fluid states. Theoretical analysis shows that the proposed second-order modified ghost fluid method (2nd-MGFM) can effectively eliminate the first order major error term occurring to the interface and accumulating with time when there is interface acceleration. Numerical results exhibit the proposed 2nd-MGFM can suppress overheating at the accelerating wall and pressure dislocation at the accelerating interface very well.

Appendices

Appendix

Basic Differential Relations for Multi-medium GRP with the Geometrical Source

Fig. 11

Wave pattern of the multi-medium GRP with the geometrical source (20)

According to the Riemann invariants $\psi$ and $\phi$ corresponding to $\lambda _{-}=u-c$ and $\lambda _{+}=u+c$ waves as Fig. 11, we have

$$\begin{aligned} \begin{aligned} d \psi =d u+\frac{1}{\rho c} d p+\frac{1}{c} T d S,\quad d \phi =d u-\frac{1}{\rho c} d p-\frac{1}{c} T d S. \end{aligned} \end{aligned}$$

(52)

Combining with (4), we obtain

$$\begin{aligned} \begin{aligned} \frac{D\psi }{Dt}=\frac{Du}{Dt}+\frac{1}{\rho c} \frac{Dp}{Dt},\quad \frac{D\phi }{Dt}=\frac{Du}{Dt}-\frac{1}{\rho c} \frac{Dp}{Dt}. \end{aligned} \end{aligned}$$

(53)

Thus, for the state ahead of the left characteristic waves $\lambda _{-}=u-c$ and the state ahead of the right characteristic waves $\lambda _{+}=u+c$, with simple manipulation,

$$\begin{aligned} \begin{aligned}&\left( \frac{D\psi }{Dt}\right) _L=\left( \frac{Du}{Dt}\right) _L+\frac{1}{\rho _L c_L} \left( \frac{Dp}{Dt}\right) _L=-\frac{1}{\rho _{L}} p_{L}^{\prime }-c_{L} u_{L}^{\prime }-\frac{m-1}{x_{c d}} c_{L} u_{L},\\&\left( \frac{D\phi }{Dt}\right) _R=\left( \frac{Du}{Dt}\right) _R-\frac{1}{\rho _R c_R} \left( \frac{Dp}{Dt}\right) _R=-\frac{1}{\rho _{R}} p_{R}^{\prime }+c_{R} u_{R}^{\prime }+\frac{m-1}{x_{cd}} c_{R} u_{R}, \end{aligned} \end{aligned}$$

(54)

which is only related to the initial value $({{{\varvec{U}}}_{L}},{{{\varvec{U}}}^{\prime }_{L}}, {{{\varvec{U}}}_{R}}$, ${{{\varvec{U}}}^{\prime }_{R}})$. According to boundary conditions Theorem 1, there is

$$\begin{aligned} \left( \frac{Du}{Dt} \right) _{*}=\left( \frac{Du}{Dt} \right) _{1*}=\left( \frac{Du}{Dt} \right) _{2*},\quad \left( \frac{Dp}{Dt} \right) _{*}=\left( \frac{Dp}{Dt} \right) _{1*}=\left( \frac{Dp}{Dt} \right) _{2*}. \end{aligned}$$

(55)

For a general EOS $e=e(\rho ,p)$, (5) can be written as

$$\begin{aligned} TdS=\frac{\partial e}{\partial p}dp+(\frac{\partial e}{\partial \rho } -\frac{p}{\rho ^2})d\rho \end{aligned}$$

(56)

In particular, for the EOS (2), (5) can be written as

$$\begin{aligned} TdS = \frac{1}{(\gamma -1)\rho }dp-\frac{c^2}{(\gamma -1)\rho }d{\rho }. \end{aligned}$$

(57)

As described, we can get instantaneous predicted state ${\varvec{U}}({{\rho }_{1*}},{{u}_{1*}},{{p}_{1*}})$, ${\varvec{U}}({{\rho }_{2*}},{{u}_{2*}},{{p}_{2*}})$ via solving the multi-medium Riemann problem (19).

Proof of Theorem 4

We have the following approximation:

$$\begin{aligned} \left( \frac{D\psi }{Dt}\right) _{1*} = \left( \frac{D\psi }{Dt}\right) _{L},\quad \left( \frac{D\phi }{Dt}\right) _{2*} = \left( \frac{D\phi }{Dt}\right) _{R}. \end{aligned}$$

(58)

Thus, (39) and (40) can be obtained by combining (53), (54), (55), (58).

Under approximation (58), there is

$$\begin{aligned} \begin{aligned} (TdS)_{1*}=(TdS)_{L},\quad (TdS)_{2*}=(TdS)_{R},\\ \end{aligned} \end{aligned}$$

(59)

(41) is obtained.

In particular, (41) can be written as (42) for the EOS (2) specifically.

Proof of Theorem 5

Here, assuming that the solution model is a double rarefaction wave structure.

According to [29], the material derivative of $\psi$ corresponding to the state behind of left characteristic wave $\lambda _{-}=u-c$ satisfies

$$\begin{aligned} \left( \frac{D\psi }{Dt}\right) _{1*}&\!=\!\left[ \frac{1+\mu _1^{2}}{1+2 \mu _1^{2}}\left( \frac{c_{1 *}}{c_{L}}\right) ^{1 /\left( 2 \mu _1^{2}\right) }+\frac{\mu _1^{2}}{1+2 \mu _1^{2}}\left( \frac{c_{1 *}}{c_{L}}\right) ^{\left( 1+\mu _1^{2}\right) / \mu _1^{2}}\right] T_{L} S_{L}^{\prime }-c_{L}\left( \frac{c_{1*}}{c_{L}}\right) ^{1 /\left( 2 \mu _1^{2}\right) } \psi _{L}^{\prime }\nonumber \\&\quad +\frac{m-1}{2 x_{cd}} c_{1 *}\left[ Z_{L}\left( c_{1 *}\right) -u_{*}-\left( \frac{c_{1*}}{c_{L}}\right) ^{\frac{1-2 \mu _1^{2}}{2 \mu _1^{2}}} u_{L}\right] , \end{aligned}$$

(60)

where, $\mu _1^{2}=\frac{\gamma _1-1}{\gamma _1+1}$, and

$$\begin{aligned} Z_{L}\left( c_{1^{*}}\right) =\left\{ \begin{array}{ll} c_{L}-c_{1^{*}}+\log \left( \frac{c_{1^{*}}}{c_{L}}\right) \left( u_{L}+\frac{2 c_{L}}{\gamma _1-1}\right) , &{}\gamma _1=3, \\ -2\left[ 3 c_{1*} \log \left( \frac{c_{1 *}}{c_{L}}\right) +\left( 1-\frac{c_{1 *}}{c_{L}}\right) \left( u_{L}+\frac{2 c_{L}}{\gamma _1-1}\right) \right] , &{} \gamma _1=\frac{5}{3}.\\ \frac{\left( \mu _1^{2}-1\right) c_{1 *}}{\mu _1^{2}\left( 4 \mu _1^{2}-1\right) }\left[ 1-\left( \frac{c_{1 *}}{c_{L}}\right) ^{\frac{1-4 \mu _1^{2}}{2 \mu _1^{2}}}\right] +\frac{\left( u_{L}+\frac{2 c_{L}}{\gamma _1-1}\right) }{2 \mu _1^{2}-1}\left[ 1-\left( \frac{c_{1 *}}{c_{L}}\right) ^{\frac{1-2 \mu _1^{2}}{2 \mu _1^{2}}}\right] , &{} \gamma _1\ne 3 \quad and\\ &{}\gamma _1\ne \frac{5}{3}. \end{array}\right. \end{aligned}$$

(61)

Setting $\xi =\frac{c_{1 *}}{c_{L}}-1$, we approximately have

$$\begin{aligned} \left( \frac{D \psi }{D t}\right) _{1 *}=-\frac{1}{\rho _{L}} p_{L}^{\prime }-c_{L} u_{L}^{\prime }-\frac{m-1}{x_{c d}} c_{L} u_{L}+Q_{1} \xi +O\left( \xi ^2\right) , \end{aligned}$$

(62)

via Taylor series expansion. Here,

$$\begin{aligned} Q_{1}=\frac{1}{2 \mu _1^{2}}\left[ -\frac{\gamma _1}{\gamma _1+1} \frac{\rho _{L} p_{L}^{\prime }+p_{L} \rho _{L}^{\prime }}{\rho _{L}^{2}}-u_{L}^{\prime } c_{L}+\frac{m-1}{2 x_{cd}}\left( -2 \mu _1^{2} u_{L} c_{L}+\frac{4 c_{L}^{2}}{\gamma _1+1}\right) \right] . \end{aligned}$$

(63)

Similarly, setting $\eta = \frac{c_{2 *}}{c_{R}}-1$, the material derivative of $\phi$ corresponding to the state of behind of the right characteristic wave $\lambda _{+}=u+c$ satisfies

$$\begin{aligned} \left( \frac{D\phi }{Dt}\right) _{2*} =-\frac{1}{\rho _{R}} p_{R}^{\prime }+c_{R} u_{R}^{\prime }+\frac{m-1}{x_{cd}} c_{R} u_{R}+Q_{2} \eta +O\left( \eta ^2\right) , \end{aligned}$$

(64)

where,

$$\begin{aligned} Q_{2}=\frac{1}{2\mu _2^{2}}\left[ -\frac{\gamma _2}{\gamma _2+1} \frac{\rho _{R} p_{R}^{\prime }+p_{R} \rho _{R}^{\prime }}{\rho _{R}^{2}}+u_{R}^{\prime } c_{R}+\frac{m-1}{2 x_{cd}}\left( 2 \mu _2^{2} u_{R} c_{R}+\frac{4 c_{R}^{2}}{\gamma _2+1}\right) \right] , \end{aligned}$$

(65)

and $\mu _2^{2}=\frac{\gamma _2-1}{\gamma _2+1}$. As a result of (53), (54), (55), (62), (64), we obtain

$$\begin{aligned} \begin{aligned}&\left( \frac{Du}{Dt}\right) _{1*}+\frac{1}{\rho _{1*} c_{1*}}\left( \frac{Dp}{Dt}\right) _{1*}=\left( \frac{D\psi }{Dt}\right) _{1*}=\left( \frac{D\psi }{Dt}\right) _{L}+Q_{1} \xi +O\left( \xi ^{2}\right) ,\\&\left( \frac{Du}{Dt}\right) _{2*}-\frac{1}{\rho _{2*} c_{2*}}\left( \frac{Dp}{Dt}\right) _{2*}=\left( \frac{D\phi }{Dt}\right) _{2*}=\left( \frac{D\phi }{Dt}\right) _{R}+Q_{2} \eta +O\left( \eta ^{2}\right) . \end{aligned} \end{aligned}$$

(66)

Solving above systems, we obtain the acceleration (material derivative), ${u}_{a}=\left( \frac{Du}{Dt} \right) _{*}=\left( \frac{Du}{Dt} \right) _{1*}=\left( \frac{Du}{Dt} \right) _{2*}$ and ${p}_{a}=\left( \frac{Dp}{Dt} \right) _{*}=\left( \frac{Dp}{Dt} \right) _{1*}=\left( \frac{Dp}{Dt} \right) _{2*}$, of the velocity and pressure at the interface $x=x_{cd}(t_0)$. Using (4), we can get the spatial derivatives of velocity and pressure at the interface

$$\begin{aligned} \begin{aligned} \begin{array}{ll} p^{\prime }_{1*}=-{{\rho }_{1*}}{{u}_{a}},\\ \left( x_{cd}^{m-1}u\right) ^{\prime }_{1*}=-\frac{{{p}_{a}}}{{{\rho }_{1*}}c_{1*}^{2}} x_{cd}^{m-1},\\ \end{array}\quad \begin{array}{ll} p^{\prime }_{2*}=-{{\rho }_{2*}}{{u}_{a}},\\ \left( x_{cd}^{m-1}u\right) ^{\prime }_{2*}=-\frac{p_a}{\rho _{2*}c_{2*}^{2}} x_{cd}^{m-1}.\\ \end{array} \end{aligned} \end{aligned}$$

(67)

To obtain the spatial derivative of density at the interface for the EOS (2), we follow the method developed in [29] and have

$$\begin{aligned} \begin{aligned} (TdS)_{1*}=\left( \frac{c_{1*}}{c_L}\right) ^{(1+\mu _1^2)/\mu _1^2}(TdS)_{L},\quad (TdS)_{2*}=\left( \frac{c_{2*}}{c_R}\right) ^{(1+\mu _2^2)/\mu _2^2}(TdS)_{R}, \end{aligned} \end{aligned}$$

(68)

which can be written as

$$\begin{aligned} \begin{aligned} (TdS)_{1*}=\left( 1+\frac{1+\mu _1^2}{\mu _1^2}\xi +O\left( \xi ^2\right) \right) (TdS)_{L},\\ (TdS)_{2*}=\left( 1+\frac{1+\mu _2^2}{\mu _2^2}\eta +O\left( \eta ^2\right) \right) (TdS)_{R}. \end{aligned} \end{aligned}$$

(69)

Based on (57), we thus have

$$\begin{aligned} \begin{aligned} \rho ^{\prime }_{1*}&=\frac{1}{{c_{1*}^{2}}}\left[ {p^{\prime }_{1*}}-\left( 1+\frac{1+\mu _1^2}{\mu _1^2}\xi +O\left( \xi ^2\right) \right) \frac{{\rho }_{1*}}{{\rho }_{L} } ({p^{\prime }_L}-{c^{2}_L}{\rho ^{\prime }_L} ) \right] , \\ \rho ^{\prime }_{2*}&=\frac{1}{{c_{2*}^{2}}}\left[ {p^{\prime }_{2*}}-\left( 1+\frac{1+\mu _2^2}{\mu _2^2}\eta +O\left( \eta ^2\right) \right) \frac{{\rho }_{2*}}{{\rho }_{R} } ({p^{\prime }_R}-{c^{2}_R}{\rho ^{\prime }_R}) \right] . \end{aligned} \end{aligned}$$

(70)

Finally, (43) can be obtained with (67) and (70).