Hydrodynamics Systems with a General EOS

This notebook explores the eigensystem of the Euler equations augmented with an additional thermodynamic variable to describe a general equation of state in the reconstruction of interface states


In [1]:
from sympy import init_session
init_session()


IPython console for SymPy 1.0 (Python 3.5.2-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.0/

In [2]:
from sympy.abc import rho, tau, alpha
rho, tau, c, h, p = symbols("rho tau c h p", real=True, positive=True)
re = symbols(r"(\rho{}e)", real=True, positive=True)
ge = symbols(r"\gamma_e", real=True, positive=True)
alpha, u = symbols("alpha u", real=True)

The two routines below simplify the analysis of the eigensystem by simultaneously finding the set of orthonormal left and right eigenvectors of a matrix from the primitive form of the hydro equations


In [3]:
plus = u + c
zero = u
minus = u - c

class Eigenvector(object):
    def __init__(self, name, ev, r, l=None):
        self.name = name
        if name == "minus":
            self.d = 0
        elif name == "zero":
            self.d = 1
        elif name == "plus":
            self.d = 2
        else:
            self.d = None

        self.ev = ev
        self.l = l
        self.r = r

    def __lt__(self, other):
        return self.d < other.d

    def __str__(self):
        return "{} wave, r = {}, l = {}".format(self.eigenvalue, self.r, self.l)

def eigensystem(A, suba=None, subb=None):
    # get the left and right eigenvectors that diagonalize the system.
    # it is best to use sympy diagonalize() for this purpose than getting
    # the left and right eigenvectors independently.

    e = []

    R, D = A.diagonalize()

    # the columns of R are the right eigenvectors and the diagonal
    # element of D is the corresponding eigenvalues

    for n in range(A.shape[0]):
        r = R.col(n)
        ev = D[n,n]
                    
        #print("here", r, ev)
        if suba is not None and subb is not None:
            ev = ev.subs(suba, subb)
            
        # which eigenvalue are we?
        if simplify(ev - minus) == 0:
            name = "minus"
        elif simplify(ev - plus) == 0:
            name = "plus"
        elif simplify(ev - zero) == 0:
            name = "zero"
        else:
            return None
                                                                                                            
        # normalize the right eigenvector
        v = r[0]
        if v != 0:
            r = r/v

        if suba is not None and subb is not None:
            r = simplify(r.subs(suba, subb))
            
        e.append(Eigenvector(name=name, ev=ev, r=r))

    # now sort the system from smallest (u-c) to largest (u+c)
    e.sort()

    # now let's construct the R with this sorting
    for n in range(A.shape[0]):
        R[:,n] = e[n].r

    # the left eigenvector matrix, L, is just the inverse
    L = R**-1

    for n in range(A.shape[0]):
        e[n].l = L.row(n)

    return e

Euler Equations with $(\rho e)$

The Euler equations in primitive variable form, $q = (\rho, u, p, (\rho e))^\intercal$ appear as: \begin{align} \frac{\partial \rho}{\partial t} &= -u \frac{\partial \rho}{\partial x} - \rho \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \frac{1}{\rho} \frac{\partial p}{\partial x} \\ \frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \rho c^2 \frac{\partial u}{\partial x} \\ \frac{\partial (\rho e)}{\partial t} &= -u \frac{\partial (\rho e)}{\partial x} - \rho h \frac{\partial u}{\partial x} \end{align}

In vector form, we have: $$q_t + A(q) q_x = 0$$

with the matrix $A(q)$:

$$A(q) = \left ( \begin{array}{cccc} u & \rho & 0 & 0\\ 0 & u & 1/\rho & 0\\ 0 & \rho c^2 & u & 0\\ 0 & \rho h & 0 & u\end{array} \right ) $$

The sound speed is related to the adiabatic index, $\Gamma_1$, as $c^2 = \Gamma_1 p /\rho$.

We can represent this matrix symbolically in SymPy and explore its eigensystem.


In [4]:
q = Matrix([rho, u, p, re]).transpose()
A = Matrix([[u, rho, 0, 0], [0, u, rho**-1, 0], [0, c**2 * rho, u, 0], [0, rho*h, 0, u]])
A


Out[4]:
$$\left[\begin{matrix}u & \rho & 0 & 0\\0 & u & \frac{1}{\rho} & 0\\0 & c^{2} \rho & u & 0\\0 & h \rho & 0 & u\end{matrix}\right]$$

The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.


In [5]:
A.eigenvals()


Out[5]:
$$\left \{ u : 2, \quad - c + u : 1, \quad c + u : 1\right \}$$

We see that there are 2 eigenvalues $u$ -- the addition of $(\rho e)$ to the system adds this degeneracy.

Eigenvectors

The right eigenvectors are defined for a given eigenvalue, $\lambda$, as: $$A r = \lambda r$$ and the left eigenvectors satisfy: $$l A = \lambda l$$

Note that the left and right eigenvectors are orthogonal to those corresponding to a different eigenvalue, and usually normalized so: $$l^i \cdot r^j = \delta_{ij}$$


In [6]:
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A)

In [7]:
# printing them out for inspection
from sympy.printing.mathml import mathml
for e in eigen:
    print(e.name)
    display(e.r, e.l)


minus
$$\left[\begin{matrix}1\\- \frac{c}{\rho}\\c^{2}\\h\end{matrix}\right]$$
$$\left[\begin{matrix}0 & - \frac{\rho}{2 c} & \frac{1}{2 c^{2}} & 0\end{matrix}\right]$$
zero
$$\left[\begin{matrix}1\\0\\0\\0\end{matrix}\right]$$
$$\left[\begin{matrix}1 & 0 & - \frac{1}{c^{2}} & 0\end{matrix}\right]$$
zero
$$\left[\begin{matrix}0\\0\\0\\1\end{matrix}\right]$$
$$\left[\begin{matrix}0 & 0 & - \frac{h}{c^{2}} & 1\end{matrix}\right]$$
plus
$$\left[\begin{matrix}1\\\frac{c}{\rho}\\c^{2}\\h\end{matrix}\right]$$
$$\left[\begin{matrix}0 & \frac{\rho}{2 c} & \frac{1}{2 c^{2}} & 0\end{matrix}\right]$$

$\beta$'s and final update

The final interface state is writen by projecting the jump in primitive variables, $\Delta q$, into characteristic variables (as $l \cdot \Delta q$), and then adding up all the jumps that reach the interface.

The convention is to write $\beta^\nu = l^\nu \cdot \Delta q$, where the superscript identifies which eigenvalue (and corresponding eigenvectors) we are considering. Note, that often a reference state is used, and the jump, $\Delta q$, will be the difference with respect to this reference state. For PPM, the $\Delta q$ will take the form of the integral under the parabola over the range that each wave can reach.

The final interface state is then: $$ q_\mathrm{int} = q - \sum_\nu \beta^\nu r^\nu $$

The tracing projects the primitive variables into characteristic variables by defining $$ \Delta q = \left ( \begin{array}{c} \Delta \rho \\ \Delta u \\ \Delta p \\ \Delta (\rho e) \end{array} \right ) $$ and then $\beta^\nu = l^\nu \cdot \Delta q$


In [8]:
from sympy.abc import delta

dr = symbols(r"\Delta\rho")
du = symbols(r"\Delta{}u")
dp = symbols(r"\Delta{}p")
dre = symbols(r"\Delta(\rho{}e)")

rhoi = symbols(r"\rho_\mathrm{int}")
ui = symbols(r"u_\mathrm{int}")
pri = symbols(r"p_\mathrm{int}")
rei = symbols(r"(\rho{}e)_\mathrm{int}")

# this is the jump
dq = Matrix([[dr, du, dp, dre]]).transpose()

# this is the interface state vector
qint = Matrix([[rhoi, ui, pri, rei]]).transpose()

Now compute the $\beta$s


In [10]:
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
    betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]

for n, e in enumerate(eigen):
    print(e.name)
    beta = e.l.dot(dq)
    display(Eq(betas[n], simplify(beta)))


minus
$$\beta^- = \frac{1}{2 c^{2}} \left(\Delta{}p - \Delta{}u c \rho\right)$$
zero
$$\beta^0_0 = \Delta\rho - \frac{\Delta{}p}{c^{2}}$$
zero
$$\beta^0_1 = \Delta(\rho{}e) - \frac{\Delta{}p h}{c^{2}}$$
plus
$$\beta^+ = \frac{1}{2 c^{2}} \left(\Delta{}p + \Delta{}u c \rho\right)$$

and now the final interface states


In [11]:
for n in range(len(eigen)):
    rhs = q[n]
    for m in range(len(eigen)):
        rhs -= betas[m]*eigen[m].r[n]
    display(Eq(qint[n],rhs))


$$\rho_\mathrm{int} = - \beta^+ - \beta^- - \beta^0_0 + \rho$$
$$u_\mathrm{int} = - \frac{\beta^+ c}{\rho} + \frac{\beta^- c}{\rho} + u$$
$$p_\mathrm{int} = - \beta^+ c^{2} - \beta^- c^{2} + p$$
$$(\rho{}e)_\mathrm{int} = (\rho{}e) - \beta^+ h - \beta^- h - \beta^0_1$$

Euler Equations with $(\gamma_e)$

We can define $\gamma_e = p/(\rho e) + 1$ and differentiate it to get: $$ \frac{\partial \gamma_e}{\partial t} = -u \frac{\partial \gamma_e}{\partial x} + (\gamma_e - 1)(\gamma_e - \Gamma_1) \frac{\partial u}{\partial x} $$ The original CG paper used $\tau = 1/\rho$ in place of density. With this, the continuity equation becomes: $$ \frac{\partial \tau}{\partial t} = -u\frac{\partial \tau}{\partial x} + \tau \frac{\partial u}{\partial x} $$

The Euler equations with this set of primitive variables, $q = (\tau, u, p, \gamma_e)^\intercal$ appear as: \begin{align} \frac{\partial \tau}{\partial t} &= -u\frac{\partial \tau}{\partial x} + \tau \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \frac{1}{\rho} \frac{\partial p}{\partial x} \\ \frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \rho c^2 \frac{\partial u}{\partial x} \\ \frac{\partial \gamma_e}{\partial t} &= -u \frac{\partial \gamma_e}{\partial x} + (\gamma_e - 1)(\gamma_e - \Gamma_1) \frac{\partial u}{\partial x} \end{align}

For convenience, we define $$ \alpha = (\gamma_e - 1)(\gamma_e - \Gamma_1) $$ and then in vector form, we have: $$q_t + A(q) q_x = 0$$

with the matrix $A(q)$:

$$A(q) = \left ( \begin{array}{cccc} u & -\tau & 0 & 0\\ 0 & u & \tau & 0\\ 0 & c^2/\tau & u & 0\\ 0 & -\alpha & 0 & u\end{array} \right ) $$

We can represent this matrix symbolically in SymPy and explore its eigensystem.


In [12]:
q = Matrix([[tau, u, p, ge]]).transpose()

A = Matrix([[u, -tau, 0, 0], [0, u, tau, 0], [0, c**2/tau, u, 0], [0, -alpha, 0, u]])
A


Out[12]:
$$\left[\begin{matrix}u & - \tau & 0 & 0\\0 & u & \tau & 0\\0 & \frac{c^{2}}{\tau} & u & 0\\0 & - \alpha & 0 & u\end{matrix}\right]$$

The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.


In [13]:
A.eigenvals()


Out[13]:
$$\left \{ u : 2, \quad - c + u : 1, \quad c + u : 1\right \}$$

We see that there are 2 eigenvalues $u$ -- the addition of $\gamma_e$ to the system adds this degeneracy.

Eigenvectors


In [14]:
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A)

In [15]:
# printing them out for inspection
from sympy.printing.mathml import mathml
for e in eigen:
    print(e.name)
    display(e.r, e.l)


minus
$$\left[\begin{matrix}1\\\frac{c}{\tau}\\- \frac{c^{2}}{\tau^{2}}\\\frac{\alpha}{\tau}\end{matrix}\right]$$
$$\left[\begin{matrix}0 & \frac{\tau}{2 c} & - \frac{\tau^{2}}{2 c^{2}} & 0\end{matrix}\right]$$
zero
$$\left[\begin{matrix}1\\0\\0\\0\end{matrix}\right]$$
$$\left[\begin{matrix}1 & 0 & \frac{\tau^{2}}{c^{2}} & 0\end{matrix}\right]$$
zero
$$\left[\begin{matrix}0\\0\\0\\1\end{matrix}\right]$$
$$\left[\begin{matrix}0 & 0 & \frac{\alpha \tau}{c^{2}} & 1\end{matrix}\right]$$
plus
$$\left[\begin{matrix}1\\- \frac{c}{\tau}\\- \frac{c^{2}}{\tau^{2}}\\\frac{\alpha}{\tau}\end{matrix}\right]$$
$$\left[\begin{matrix}0 & - \frac{\tau}{2 c} & - \frac{\tau^{2}}{2 c^{2}} & 0\end{matrix}\right]$$

$\beta$'s and final update

The tracing projects the primitive variables into characteristic variables by defining $$ \Delta q = \left ( \begin{array}{c} \Delta \tau \\ \Delta u \\ \Delta p \\ \Delta \gamma_e \end{array} \right ) $$ and then $\beta^\nu = l^\nu \cdot \Delta q$


In [16]:
dtau = symbols(r"\Delta\tau")
dge = symbols(r"\Delta\gamma_e")

dq = Matrix([[dtau, du, dp, dge]]).transpose()

taui = symbols(r"\tau_\mathrm{int}")
gei = symbols(r"(\gamma_e)_\mathrm{int}")

qint = Matrix([[taui, ui, pri, gei]]).transpose()

In [18]:
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
    betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]

for n, e in enumerate(eigen):
    print(e.name)
    beta = e.l.dot(dq)
    display(Eq(betas[n], simplify(beta)))


minus
$$\beta^- = \frac{\tau}{2 c^{2}} \left(- \Delta{}p \tau + \Delta{}u c\right)$$
zero
$$\beta^0_0 = \Delta\tau + \frac{\Delta{}p \tau^{2}}{c^{2}}$$
zero
$$\beta^0_1 = \Delta\gamma_e + \frac{\Delta{}p \alpha}{c^{2}} \tau$$
plus
$$\beta^+ = - \frac{\tau}{2 c^{2}} \left(\Delta{}p \tau + \Delta{}u c\right)$$

and now the final interface states


In [19]:
for n in range(len(eigen)):
    rhs = q[n]
    for m in range(len(eigen)):
        rhs -= betas[m]*eigen[m].r[n]
    display(Eq(qint[n],rhs))


$$\tau_\mathrm{int} = - \beta^+ - \beta^- - \beta^0_0 + \tau$$
$$u_\mathrm{int} = \frac{\beta^+ c}{\tau} - \frac{\beta^- c}{\tau} + u$$
$$p_\mathrm{int} = \frac{\beta^+ c^{2}}{\tau^{2}} + \frac{\beta^- c^{2}}{\tau^{2}} + p$$
$$(\gamma_e)_\mathrm{int} = - \frac{\beta^+ \alpha}{\tau} - \frac{\beta^- \alpha}{\tau} - \beta^0_1 + \gamma_e$$

Gray FLD Radiation Euler Equations with $(\rho e)$

Following Zhang et al. (2011), the equations of gray radiation hydrodynamics with primitive variables $q = (\rho, u, p, (\rho e)_g, E_r)^\intercal$ are: \begin{align} \frac{\partial \rho}{\partial t} &= -u\frac{\partial \rho}{\partial x} - \rho \frac{\partial u}{\partial x} \ \frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \frac{1}{\rho} \frac{\partial p}{\partial x}

- \frac{\lambda_f}{\rho} \frac{\partial E_r}{\partial x}\\

\frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \rho c_g^2 \frac{\partial u}{\partial x} \ \frac{\partial (\rho e)_g}{\partial t} &= -u \frac{\partial (\rho e)_g}{\partial x} - \rho h \frac{\partial u}{\partial x} \ \frac{\partial E_r}{\partial t} &= -\frac{3-f}{2} E_r \frac{\partial u}{\partial x} - \left ( \frac{3-f}{2} - \lambda_f\right ) u \frac{\partial E_r}{\partial x} \end{align} where $(\rho e)_g$ is the gas internal energy density, $h_g$ is the gas specific enthalpy, $c_g$ is the gas sound speed (obeying $c_g^2 = \Gamma_1 p /\rho$), $E_r$ is the radiation energy density, $f$ is the Eddington factor, and $\lambda_f$ is the flux limiter.

Following Zhang et al., we make the approximation that $$ \frac{3-f}{2} = \lambda_f + 1 $$ and then in vector form, we have: $$q_t + A(q) q_x = 0$$

with the matrix $A(q)$:

$$A(q) = \left ( \begin{array}{ccccc} u & \rho & 0 & 0 & 0 \\ 0 & u & 1/\rho & 0 & \lambda_f/\rho\\ 0 & \rho c_g^2 & u & 0 & 0 \\ 0 & \rho h_g & 0 & u & 0 \\ 0 & (3-f)E_r/2 & 0 & 0 & u\end{array} \right ) $$

We can represent this matrix symbolically in SymPy and explore its eigensystem.


In [20]:
cg = symbols(r"c_g", real=True, positive=True)
hg = symbols(r"h_g", real=True, positive=True)
Er = symbols(r"E_r", real=True, positive=True)
lf = symbols(r"\lambda_f", real=True)
f = symbols("f", real=True, positive=True)

q = Matrix([[rho, u, p, re, Er]]).transpose()
A = Matrix([[u, rho, 0, 0, 0], 
            [0, u, rho**-1, 0, lf/rho], 
            [0, rho*cg**2, u, 0, 0], 
            [0, rho*hg, 0, u, 0],
            [0, (lf+1)*Er, 0, 0, u]])
A


Out[20]:
$$\left[\begin{matrix}u & \rho & 0 & 0 & 0\\0 & u & \frac{1}{\rho} & 0 & \frac{\lambda_f}{\rho}\\0 & c_{g}^{2} \rho & u & 0 & 0\\0 & h_{g} \rho & 0 & u & 0\\0 & E_{r} \left(\lambda_f + 1\right) & 0 & 0 & u\end{matrix}\right]$$

The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.


In [21]:
A.eigenvals()


Out[21]:
$$\left \{ u : 3, \quad \frac{1}{\sqrt{\rho}} \left(\sqrt{\rho} u - \sqrt{E_{r} \lambda_f^{2} + E_{r} \lambda_f + c_{g}^{2} \rho}\right) : 1, \quad \frac{1}{\sqrt{\rho}} \left(\sqrt{\rho} u + \sqrt{E_{r} \lambda_f^{2} + E_{r} \lambda_f + c_{g}^{2} \rho}\right) : 1\right \}$$

We see that there are 3 eigenvalues $u$. We identify the total sound speed (radiation + gas) as: $$c^2 = c_g^2 + (\lambda_f + 1)\frac{\lambda_f E_r}{\rho}$$ We can simplify these by substituting in that relationship


In [22]:
cc = c**2 - (lf +1)*lf*Er/rho

evs = A.eigenvals()
for e in evs.keys():
    display(powsimp(simplify(e.subs(cg**2, cc))))


$$u$$
$$c + u$$
$$- c + u$$

We see that these are the same form of the eigenvalues for the pure hydrodynamics system

Eigenvectors


In [23]:
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A, suba=cg, subb=sqrt(cc))

In [24]:
# printing them out for inspection
for e in eigen:
    print(e.name)
    display(e.r)
    display(simplify(e.l))


minus
$$\left[\begin{matrix}1\\- \frac{c}{\rho}\\\frac{1}{\rho} \left(- E_{r} \lambda_f^{2} - E_{r} \lambda_f + c^{2} \rho\right)\\h_{g}\\\frac{E_{r}}{\rho} \left(\lambda_f + 1\right)\end{matrix}\right]$$
$$\left[\begin{matrix}0 & - \frac{\rho}{2 c} & \frac{1}{2 c^{2}} & 0 & \frac{\lambda_f}{2 c^{2}}\end{matrix}\right]$$
zero
$$\left[\begin{matrix}1\\0\\0\\0\\0\end{matrix}\right]$$
$$\left[\begin{matrix}1 & 0 & - \frac{1}{c^{2}} & 0 & - \frac{\lambda_f}{c^{2}}\end{matrix}\right]$$
zero
$$\left[\begin{matrix}0\\0\\0\\1\\0\end{matrix}\right]$$
$$\left[\begin{matrix}0 & 0 & - \frac{h_{g}}{c^{2}} & 1 & - \frac{\lambda_f h_{g}}{c^{2}}\end{matrix}\right]$$
zero
$$\left[\begin{matrix}0\\0\\- \lambda_f\\0\\1\end{matrix}\right]$$
$$\left[\begin{matrix}0 & 0 & - \frac{E_{r} \left(\lambda_f + 1\right)}{c^{2} \rho} & 0 & \frac{1}{c^{2} \rho} \left(- E_{r} \lambda_f^{2} - E_{r} \lambda_f + c^{2} \rho\right)\end{matrix}\right]$$
plus
$$\left[\begin{matrix}1\\\frac{c}{\rho}\\\frac{1}{\rho} \left(- E_{r} \lambda_f^{2} - E_{r} \lambda_f + c^{2} \rho\right)\\h_{g}\\\frac{E_{r}}{\rho} \left(\lambda_f + 1\right)\end{matrix}\right]$$
$$\left[\begin{matrix}0 & \frac{\rho}{2 c} & \frac{1}{2 c^{2}} & 0 & \frac{\lambda_f}{2 c^{2}}\end{matrix}\right]$$

$\beta$'s and final update

The tracing projects the primitive variables into characteristic variables by defining $\beta^\nu = l^\nu \cdot \Delta q$


In [25]:
dEr = symbols(r"\Delta{}E_r")
Eri = symbols(r"{E_r}_\mathrm{int}")
dq = Matrix([[dr, du, dp, dre, dEr]]).transpose()

qint = Matrix([[rhoi, ui, pri, rei, Eri]])

In [27]:
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
    betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]

for n, e in enumerate(eigen):
    print(e.name)
    beta = e.l.dot(dq)
    display(Eq(betas[n], simplify(beta)))


minus
$$\beta^- = \frac{1}{2 c^{2}} \left(\Delta{}E_r \lambda_f + \Delta{}p - \Delta{}u c \rho\right)$$
zero
$$\beta^0_0 = \frac{1}{c^{2}} \left(\Delta\rho c^{2} - \Delta{}E_r \lambda_f - \Delta{}p\right)$$
zero
$$\beta^0_1 = \frac{1}{c^{2}} \left(\Delta(\rho{}e) c^{2} - \Delta{}E_r \lambda_f h_{g} - \Delta{}p h_{g}\right)$$
zero
$$\beta^0_2 = \frac{1}{c^{2} \rho} \left(- E_{r} \Delta{}E_r \lambda_f^{2} - E_{r} \Delta{}E_r \lambda_f - E_{r} \Delta{}p \lambda_f - E_{r} \Delta{}p + \Delta{}E_r c^{2} \rho\right)$$
plus
$$\beta^+ = \frac{1}{2 c^{2}} \left(\Delta{}E_r \lambda_f + \Delta{}p + \Delta{}u c \rho\right)$$

and now the final interface states


In [28]:
for n in range(len(eigen)):
    rhs = q[n]
    for m in range(len(eigen)):
        rhs -= betas[m]*eigen[m].r[n]
    display(Eq(qint[n],rhs))


$$\rho_\mathrm{int} = - \beta^+ - \beta^- - \beta^0_0 + \rho$$
$$u_\mathrm{int} = - \frac{\beta^+ c}{\rho} + \frac{\beta^- c}{\rho} + u$$
$$p_\mathrm{int} = - \frac{\beta^+}{\rho} \left(- E_{r} \lambda_f^{2} - E_{r} \lambda_f + c^{2} \rho\right) - \frac{\beta^-}{\rho} \left(- E_{r} \lambda_f^{2} - E_{r} \lambda_f + c^{2} \rho\right) + \beta^0_2 \lambda_f + p$$
$$(\rho{}e)_\mathrm{int} = (\rho{}e) - \beta^+ h_{g} - \beta^- h_{g} - \beta^0_1$$
$${E_r}_\mathrm{int} = - \frac{E_{r} \beta^+}{\rho} \left(\lambda_f + 1\right) - \frac{E_{r} \beta^-}{\rho} \left(\lambda_f + 1\right) + E_{r} - \beta^0_2$$

In [ ]:

Gray FLD Radiation Euler Equations with $(\gamma_e)$

We now look at the same system with a different auxillary thermodynamic variable (as we did with pure hydro), using $q = (\tau, u, p, {\gamma_e}_g, E_r)^\intercal$: \begin{align} \frac{\partial \tau}{\partial t} &= -u\frac{\partial \tau}{\partial x} + \tau \frac{\partial u}{\partial x} \ \frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \tau \frac{\partial p}{\partial x}

- \tau \lambda_f \frac{\partial E_r}{\partial x}\\

\frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \frac{c_g^2}{\tau} \frac{\partial u}{\partial x} \ \frac{\partial {\gamma_e}_g}{\partial t} &= -u \frac{\partial {\gamma_e}_g}{\partial x} + \alpha \frac{\partial u}{\partial x} \ \frac{\partial E_r}{\partial t} &= -\frac{3-f}{2} E_r \frac{\partial u}{\partial x} - \left ( \frac{3-f}{2} - \lambda_f\right ) u \frac{\partial E_r}{\partial x} \end{align} here, ${\gamma_e}_g$ is defined solely in terms of the gas pressure and energy, the remaining variables have the same meaning as above. We again make the approximation that $$ \frac{3-f}{2} = \lambda_f + 1 $$ and then in vector form, we have: $$q_t + A(q) q_x = 0$$

with the matrix $A(q)$:

$$A(q) = \left ( \begin{array}{ccccc} u & -\tau & 0 & 0 & 0 \\ 0 & u & \tau & 0 & \tau \lambda_f\\ 0 & c_g^2/\tau & u & 0 & 0 \\ 0 & -\alpha & 0 & u & 0 \\ 0 & (3-f)E_r/2 & 0 & 0 & u\end{array} \right ) $$

We can represent this matrix symbolically in SymPy and explore its eigensystem.


In [29]:
geg = symbols(r"{\gamma_e}_g", real=True, positive=True)
Er = symbols(r"E_r", real=True, positive=True)
lf = symbols(r"\lambda_f", real=True, positive=True)
cg = symbols(r"c_g", real=True, positive=True)

q = Matrix([[tau, u, p, geg, Er]]).transpose()
A = Matrix([[u, -tau, 0, 0, 0], 
            [0, u, tau, 0, tau*lf], 
            [0, cg**2/tau, u, 0, 0], 
            [0, -alpha, 0, u, 0],
            [0, (lf+1)*Er, 0, 0, u]])
A


Out[29]:
$$\left[\begin{matrix}u & - \tau & 0 & 0 & 0\\0 & u & \tau & 0 & \lambda_f \tau\\0 & \frac{c_{g}^{2}}{\tau} & u & 0 & 0\\0 & - \alpha & 0 & u & 0\\0 & E_{r} \left(\lambda_f + 1\right) & 0 & 0 & u\end{matrix}\right]$$

The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.


In [30]:
A.eigenvals()


Out[30]:
$$\left \{ u : 3, \quad u - \sqrt{E_{r} \lambda_f^{2} \tau + E_{r} \lambda_f \tau + c_{g}^{2}} : 1, \quad u + \sqrt{E_{r} \lambda_f^{2} \tau + E_{r} \lambda_f \tau + c_{g}^{2}} : 1\right \}$$

We see that there are 3 eigenvalues $u$. We identify the total sound speed (radiation + gas) as: $$c^2 = c_g^2 + (\lambda_f + 1)\lambda_f E_r\tau$$ We can simplify these by substituting in that relationship


In [31]:
cc = c**2 - (lf +1)*lf*Er*tau

evs = A.eigenvals()
for e in evs.keys():
    display(powsimp(simplify(e.subs(cg**2, cc))))


$$c + u$$
$$u$$
$$- c + u$$

We see that these are the same form of the eigenvalues for the pure hydrodynamics system

Eigenvectors


In [32]:
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A, suba=cg, subb=sqrt(cc))

In [33]:
# printing them out for inspection
for e in eigen:
    print(e.name)
    display(e.r)
    display(simplify(e.l))


minus
$$\left[\begin{matrix}1\\\frac{c}{\tau}\\\frac{1}{\tau^{2}} \left(E_{r} \lambda_f^{2} \tau + E_{r} \lambda_f \tau - c^{2}\right)\\\frac{\alpha}{\tau}\\- \frac{E_{r}}{\tau} \left(\lambda_f + 1\right)\end{matrix}\right]$$
$$\left[\begin{matrix}0 & \frac{\tau}{2 c} & - \frac{\tau^{2}}{2 c^{2}} & 0 & - \frac{\lambda_f \tau^{2}}{2 c^{2}}\end{matrix}\right]$$
zero
$$\left[\begin{matrix}1\\0\\0\\0\\0\end{matrix}\right]$$
$$\left[\begin{matrix}1 & 0 & \frac{\tau^{2}}{c^{2}} & 0 & \frac{\lambda_f \tau^{2}}{c^{2}}\end{matrix}\right]$$
zero
$$\left[\begin{matrix}0\\0\\0\\1\\0\end{matrix}\right]$$
$$\left[\begin{matrix}0 & 0 & \frac{\alpha \tau}{c^{2}} & 1 & \frac{\lambda_f \alpha}{c^{2}} \tau\end{matrix}\right]$$
zero
$$\left[\begin{matrix}0\\0\\- \lambda_f\\0\\1\end{matrix}\right]$$
$$\left[\begin{matrix}0 & 0 & - \frac{E_{r} \tau}{c^{2}} \left(\lambda_f + 1\right) & 0 & \frac{1}{c^{2}} \left(- E_{r} \lambda_f^{2} \tau - E_{r} \lambda_f \tau + c^{2}\right)\end{matrix}\right]$$
plus
$$\left[\begin{matrix}1\\- \frac{c}{\tau}\\\frac{1}{\tau^{2}} \left(E_{r} \lambda_f^{2} \tau + E_{r} \lambda_f \tau - c^{2}\right)\\\frac{\alpha}{\tau}\\- \frac{E_{r}}{\tau} \left(\lambda_f + 1\right)\end{matrix}\right]$$
$$\left[\begin{matrix}0 & - \frac{\tau}{2 c} & - \frac{\tau^{2}}{2 c^{2}} & 0 & - \frac{\lambda_f \tau^{2}}{2 c^{2}}\end{matrix}\right]$$

$\beta$'s and final update

The tracing projects the primitive variables into characteristic variables by defining $\beta^\nu = l^\nu \cdot \Delta q$


In [34]:
dgeg = symbols(r"\Delta{\gamma_e}_g")
dEr = symbols(r"\Delta{}E_r")
dq = Matrix([[dtau, du, dp, dgeg, dEr]]).transpose()

Eri = symbols(r"{E_r}_\mathrm{int}")
qint = Matrix([[taui, ui, pri, gei, Eri]])

In [36]:
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
    betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]

for n, e in enumerate(eigen):
    print(e.name)
    beta = e.l.dot(dq)
    display(Eq(betas[n], simplify(beta)))


minus
$$\beta^- = \frac{\tau}{2 c^{2}} \left(- \Delta{}E_r \lambda_f \tau - \Delta{}p \tau + \Delta{}u c\right)$$
zero
$$\beta^0_0 = \frac{1}{c^{2}} \left(\Delta\tau c^{2} + \Delta{}E_r \lambda_f \tau^{2} + \Delta{}p \tau^{2}\right)$$
zero
$$\beta^0_1 = \frac{1}{c^{2}} \left(\Delta{\gamma_e}_g c^{2} + \Delta{}E_r \lambda_f \alpha \tau + \Delta{}p \alpha \tau\right)$$
zero
$$\beta^0_2 = \frac{1}{c^{2}} \left(- E_{r} \Delta{}E_r \lambda_f^{2} \tau - E_{r} \Delta{}E_r \lambda_f \tau - E_{r} \Delta{}p \lambda_f \tau - E_{r} \Delta{}p \tau + \Delta{}E_r c^{2}\right)$$
plus
$$\beta^+ = - \frac{\tau}{2 c^{2}} \left(\Delta{}E_r \lambda_f \tau + \Delta{}p \tau + \Delta{}u c\right)$$

and now the final interface states


In [37]:
for n in range(len(eigen)):
    rhs = q[n]
    for m in range(len(eigen)):
        rhs -= betas[m]*eigen[m].r[n]
    display(Eq(qint[n],rhs))


$$\tau_\mathrm{int} = - \beta^+ - \beta^- - \beta^0_0 + \tau$$
$$u_\mathrm{int} = \frac{\beta^+ c}{\tau} - \frac{\beta^- c}{\tau} + u$$
$$p_\mathrm{int} = - \frac{\beta^+}{\tau^{2}} \left(E_{r} \lambda_f^{2} \tau + E_{r} \lambda_f \tau - c^{2}\right) - \frac{\beta^-}{\tau^{2}} \left(E_{r} \lambda_f^{2} \tau + E_{r} \lambda_f \tau - c^{2}\right) + \beta^0_2 \lambda_f + p$$
$$(\gamma_e)_\mathrm{int} = - \frac{\beta^+ \alpha}{\tau} - \frac{\beta^- \alpha}{\tau} - \beta^0_1 + {\gamma_e}_g$$
$${E_r}_\mathrm{int} = \frac{E_{r} \beta^+}{\tau} \left(\lambda_f + 1\right) + \frac{E_{r} \beta^-}{\tau} \left(\lambda_f + 1\right) + E_{r} - \beta^0_2$$

In [ ]: