Euler's equations for a tube with varying cross-section

Question

$\def\pd{\partial}$ $\def\l{\left}\def\r{\right}$ $\def\mdot{{\dot{m}}}$ $\def\eps{\varepsilon}$ Consider a tube with longitudinal coordinate $x$ from $0$ to $l$ and varying cross-section $A(x)$. Derek S. Bale (a former student of Prof. Randall J. LeVeque) gives in his PhD (http://faculty.washington.edu/rjl/students/dbale/thesis.ps.gz) the Euler equations for density $\rho$, impulse density $(\rho u)$ and total energy density $e$ in such a tube in the form $$ \begin{array}{ccccl} \pd_t \rho &+& \frac1{A}\pd_x(A(\rho u)) &=& 0\\ \pd_t (\rho u) &+& \frac1{A}\pd_x\l(A\l(\frac{(\rho u)^2}\rho+p\r)\r) &=& \frac{A'}{A}p\\ \pd_t e &+& \frac1{A}\pd_x\l(A(e+p)\frac{(\rho u)}{\rho}\r)&=& 0 \end{array} $$ as an example for a spatially varying flux function.

Thereby, the pressure is determined by the gamma-law $p =(\gamma-1)\l(e-\frac12 \frac{(\rho u)^2}{\rho}\r)$. Note, that I have slightly reformulated the system such that beside the pressure $p$ only conserved quantities are used.

But, is it not better to change the conserved variables to the mass per length $m':=\rho A$, the mass flow $\mdot:=A\rho u$ and the energy per length $E':=Ae$?

As an intermediate variable one could introduce a pressure force $P:=A p=(\gamma-1)\l(E'-\frac12\frac{\mdot^2}{m'}\r)$.

Euler's equations reformulated with these new conserved quantities are: $$ \begin{array}{ccccl} \pd_t m' &+& \pd_x\mdot &=& 0\\ \pd_t \mdot &+& \pd_x\l(\frac{\mdot^2}{m'}+P\r) &=& \frac{A'}{A}P\\ \pd_t E' &+& \pd_x\l((E'+P)\frac{\mdot}{m'}\r)&=& 0 \end{array} $$ In this way one obtains a system with the structure of the original Euler equations (including a source term) where the flux function does not depend on the spatial variable. Therefore, the numerics should work fine. I could use Roe's method on this problem.

Do I miss here something?

---

EDIT:

Bale uses the first form to calculate supersonic waves through a tube with a narrowing halfway. The speed-up in the narrowing causes a supersonic wave and the corresponding shock wave.

The variation of $A(x)$ is smooth: $$ A(x) = \begin{cases} 1& \text{ for }x<1\\ 1+\frac\eps2\bigl(\cos(\pi(x-1))-1\bigr)&\text{for }1\leq x \leq 3\\ 1&\text{for }x>3 \end{cases} $$

I am aware of the fact that smooth transformations of the conserved quantities, such as $(m',\mdot,E')=(A\rho,A(\rho u),A e)$ have the potential to change the structure of discontinusous solutions.

But, is this really the case here?

Answer 1

$\def\rmC{{\mathrm C}}$ $\def\vr{{\vec r}}$ $\def\vA{{\vec A}}$ $\def\ve{{\vec e}}$ $\def\l{\left}\def\r{\right}$ In the following I try to show that Bale's form is correct. If you complain, leave a comment.

The equations we can trust are the integral equations.

Mass Balance:

$$ \l[\int_0^x \rho(x',t') A(x') dx'\r]_{t'=0}^t + \int_{t'=0}^t \l[(\rho u)(x',t')A(x')\r]_{x'=0}^x d t' =0 $$ with the density $\rho(x,t)$ and the momentum flow $(\rho u)(x,t)$ as averages over $A(x)$.

Momentum balance:

$$ \l[\int_0^x (\rho u)(x',t') A(x')dx'\r]_{t'=0}^t + \ve_x\cdot\int_{t'=0}^t\l( \int_{\partial V}(\rho \vec{v})\vec{v}\cdot d\vec{A} + \int_{\partial V} p d\vec{A}\r)dt =0 $$ with the volume $V=\bigcup_{x'\in[0,x]}A(x')$.

There is no flow out of complementary part $A_\rmC:=\pd V\setminus \bigl(\{0\}\times A(0)\cup \{x\}\times A(x)\bigr)$ of the tube cross-sections $\{0\}\times A(0)$ and $\{x\}\times A(x)$. Furthermore, we approximate $\int_{A(x)} (\rho u)(\vr,t) \vec{v}(\vr,t)\cdot d\vA \approx (\rho u)(x,t) u(x,t) A(x)$ and remark that for integration in x-direction we have $\ve_x\cdot d\vA = d A_x(x')$.

That gives: $$ \begin{array}{l} \l[\int_0^x (\rho u)(x',t') A(x')dx'\r]_{t'=0}^t+ \int_{t'=0}^t \l[(\rho u)(x',t') u(x',t') A(x') + pA(x')\r]_{x'=0}^{x} dt =\\ =\quad-\int_{t'=0}^t \int_{x'=0}^x p(x',t') d A(x') dt' \end{array} $$

Thereby $\int_{x'=0}^x p(x',t') d A(x')$ can be interpreted as Riemann-Stieltjes integral which can be re-written as $\int_{x'=0}^x p(x',t') A'(x') dx'$ in the case that $A(x)$ depends smoothly on $x$.

The Riemann-Stieltjes integral suggests to use $$ p_k\frac{A_{k+\frac12}-A_{k-\frac12}}{A_k} $$ as a 1st-order update for the source term to capture discontinuities of $A(x)$. The 1st-order update of the flux function in the FVM-method is already designed such that discontinuities of $A(x)$ are captured.

Energy balance:

With a similar approximation for the averages over the cross-section as for the momentum balance and with the zero flux out of $A_\rmC$ we get: $$ \l[\int_0^x e(x',t')A(x')dx'\r]_{t'=0}^t + \int_{t'=0}^t \biggl[(e(x',t') + p(x',t'))u(x')A(x')\biggr]_{x'=0}^x dt = 0 $$

Interpretation

LeVeque describes in chapter 6.16 "Capacity-Form Differencing" of his FVM-book that conservation laws of the form $$ \frac{d}{dt}\int_{x_1}^{x_2} \kappa(x) q(x,t) dx = f(q(x_1,t))-f(q(x_2,t)) $$ with a capacity function $\kappa(x)$ can be treated with the help of capacity-form differencing.

The more general integral form is $$ \l[\int_{x_1}^{x_2} \kappa(x) q(x,t) dx\r]_{t_1}^{t_2} = \int_{t_1}^{t_2}(f(q(x_1,t))-f(q(x_2,t))) dt. $$ In the above conservation laws $A(x')$ can be seen as capacity function.

We see that the equations in Bale's form are correct even for shock waves with capacity-form differencing.

The other form I have given in my question is also correct. For that case the conserved quantities stand right under the spatial integral on the left hand side. Everything else is just combination of terms without differentiation.

Now I believe that Bale used his form of the equations to have a simple example for a system with spatially varying flux function.