7.3. Compressible Euler Equations

7.3.1. Introduction

The compressible Euler equations are equations for perfect fluid. Perfect fluids have no heat conduction (T^{i0} = T^{0i} = 0) and no viscosity (T^{ij} = p\one), so in the comoving frame the stress energy tensor is:

T^{\alpha\beta} = \diag(\rho c^2, p, p, p) = \left(\rho+{p\over c^2}\right)u^\alpha u^\beta + p g^{\alpha\beta}

(we use g^{\mu\nu} = \diag(-1, 1, 1, 1)). Relativistic Euler equations are given by the conservation of the stress energy tensor and the particle number conservation:

\partial_\nu T^{\mu\nu} = 0 \partial_\mu(nu^\mu) = 0

By doing the nonrelativistic limit (see Perfect Fluids for a detailed derivation), we get the following Euler equations:

{\partial\rho\over\partial t} + \nabla\cdot(\rho{\bf u}) = 0 {\partial(\rho{\bf u})\over\partial t} + \nabla\cdot(\rho{\bf u}{\bf u}^T) + \nabla p - {\bf f} = 0 {\partial E\over\partial t} + \nabla\cdot({\bf u}(E+p)) = 0

where

E = \rho e + \half \rho u^2

is the total energy per unit volume, composed of the kinetic energy per unit volume (\half \rho u^2) and the internal energy per unit volume (\rho e), where e is the internal energy per unit mass (e={U\over nM}). The energy E doesn’t contain the rest mass energy, but all other energies are hidden in the internal energy.

We use the ideal gas equations, so:

e = T c_v p = {n\over V} \bar RT = {n M\over V} {\bar R\over M}T = \rho RT = \rho R {e\over c_v} = {R\over c_v} (E-\half \rho u^2)

where n is the number of moles of gas, M is the molar mass of the gas (i.e. a mass of one mole of the gas, e.g. for dry air we get M=28.956\rm\,g/mol, as it is mainly composed of 20% of oxygen with atomic mass 16 and 78% of nitrogen with atomic mass 14, both form diatomic molecules, so the molecular mass is twice the atomic mass giving the total of 0.2\cdot2\cdot16+0.78\cdot2\cdot14 = 28.24, the rest is given by the other components and one also has to average over all isotopes), \bar R=N_A k_B\doteq8.3145\rm\,{J\over K\, mol} is the ideal gas constant (N_A is the Avogadro constant and k_B is the Boltzmann constant), R={\bar R\over M} is the specific ideal gas constant (e.g. for dry air we get R = {8.3145\over 28.956}\rm\,{J\over g\,K}\doteq 287.14\rm\,{J\over kg\,K}), \rho={n M\over V}={p\over RT} is the density of the gas (e.g. for dry air at the pressure 10^5\rm\,Pa and temperature 22\rm\,^\circ C we get \rho = {10^5\over 287.14 \cdot (22+273.15)}\rm\,{kg\over m^3} =1.18\,{kg\over m^3}), c_v is the specific heat capacity at constant volume (i.e. the amount of energy needed to heat one kg by one Kelvin at constant volume, e.g. for dry air the experimental value is about c_v = 717.5\rm\,{J\over kg\,K}), V is the volume and T is the temperature of the gas. Of those, V, n, M, R, \bar R are constants, \rho, e, E and T are functions of (t, x, y, z).

Here are the SI units of the various terms in the Euler equations:

[u] = \rm m\,s^{-1} [\rho] = \rm kg\,m^{-3} \rm N = kg \, m \, s^{-2} \rm J = N\, m = kg \, m^2 \, s^{-2} [p] = \rm N\, m^{-2} = kg\, m^{-1}\,s^{-2} [\half \rho u^2] = [\rho][u]^2 = \rm kg\,m^{-3}\,m^2\,s^{-2} = kg\, m^{-1}\,s^{-2} [E] = \rm J\,m^{-3} = kg\, m^{-1}\,s^{-2} [R] = \rm J\,kg^{-1}\,K^{-1} = m^2\,s^{-2}\,K^{-1} [c_v] = \rm J\,kg^{-1}\,K^{-1} = m^2\,s^{-2}\,K^{-1} [e] = {[E]\over[\rho]} = \rm {kg\,m^{-1}\,s^{-2}\over kg\,m^{-3}} = m^2\,s^{-2}

In order to calculate the specific heat ratio \kappa, we use R=c_p-c_v:

\kappa = {c_p\over c_v} = {c_v+R\over c_v} = 1 + {R\over c_v}

and the speed of sound is:

c = \sqrt{\kappa {p \over \rho}}

7.3.2. Dimensionless Euler Equations

We choose 3 constants l_r, u_r and \rho_r - characteristic length of the domain, velocity and density. Now we multiply the Euler equations with proper combinations of these constants as follows:

\left[{\partial\rho\over\partial t} + \nabla\cdot(\rho{\bf u})\right] {l_r\over\rho_r u_r} = 0 \left[{\partial(\rho{\bf u})\over\partial t} + \nabla\cdot(\rho{\bf u}{\bf u}^T) + \nabla p - {\bf f}\right] {l_r\over\rho_r u_r^2} = 0 \left[{\partial E\over\partial t} + \nabla\cdot({\bf u}(E+p))\right] {l_r\over\rho_r u_r^3} = 0

This is equal to:

{\partial\tilde\rho\over\partial \tilde t} + \tilde\nabla\cdot(\tilde\rho \tilde{\bf u}) = 0 {\partial(\tilde \rho\tilde{\bf u})\over\partial \tilde t} + \tilde\nabla\cdot(\tilde\rho\tilde{\bf u}\tilde{\bf u}^T) + \tilde\nabla\tilde p - \tilde{\bf f} = 0 {\partial\tilde E\over\partial\tilde t} + \tilde\nabla\cdot( \tilde{\bf u}(\tilde E+\tilde p)) = 0

where:

t_r = {l_r\over u_r} \tilde t = {t\over t_r} \tilde \rho = {\rho\over\rho_r} \tilde{\bf u} = {{\bf u}\over u_r} \tilde\nabla = l_r\nabla \tilde E = {E\over \rho_r u^2_r} \tilde p = {p\over \rho_r u^2_r} \tilde {\bf f} = {\bf f}{l_r\over\rho_r u^2_r}

In particular, if {\bf f}=(0, 0, -\rho g), then

\tilde{\bf f}=(0, 0, -\tilde\rho \tilde g) \tilde g = g{l_r\over u^2_r} = g{t_r^2\over l_r}

So the dimensionless Euler equations look exactly the same as the original ones, we just need to rescale all the quantities using the relations above.

7.3.3. Conservative Form of the Euler Equations

We can write the Euler equations as:

{\partial{\bf w}\over \partial t} + {\partial{\bf f}_x\over \partial x} + {\partial{\bf f}_y\over \partial y} + {\partial{\bf f}_z\over \partial z} + {\bf g}= 0

where:

{\bf w} = \left( \begin{array}{c} \varrho\\ \rho u_1\\ \rho u_2\\ \rho u_3\\ E \end{array} \right) = \left( \begin{array}{c} w_0 \\ w_1 \\ w_2 \\ w_3 \\ w_4 \\ \end{array} \right) {\bf f}_x = \left( \begin{array}{c} \rho u_1\\ \rho u_1^2 + p\\ \rho u_1 u_2\\ \rho u_1 u_3\\ u_1(E+p) \end{array} \right) = \left( \begin{array}{c} w_1\\ \frac{w_1^2}{w_0} + p\\ \frac{w_1w_2}{w_0}\\ \frac{w_1w_3}{w_0}\\ \frac{w_1}{w_0}(w_4+p) \end{array} \right) {\bf f}_y = \left( \begin{array}{c} \rho u_2\\ \rho u_2 u_1\\ \rho u_2^2 + p\\ \rho u_2 u_3\\ u_2(E+p) \end{array} \right) = \left( \begin{array}{c} w_2\\ \frac{w_2w_1}{w_0}\\ \frac{w_2^2}{w_0} + p\\ \frac{w_2w_3}{w_0}\\ \frac{w_2}{w_0}(w_4+p) \end{array} \right) {\bf f}_z = \left( \begin{array}{c} \rho u_3\\ \rho u_3 u_1\\ \rho u_3 u_2\\ \rho u_3^2 + p\\ u_3(E+p) \end{array} \right) = \left( \begin{array}{c} w_3\\ \frac{w_3w_1}{w_0}\\ \frac{w_3w_2}{w_0}\\ \frac{w_3^2}{w_0} + p\\ \frac{w_3}{w_0}(w_4+p) \end{array} \right) {\bf g} = \left( \begin{array}{c} 0\\ -f_x\\ -f_y\\ -f_z\\ 0\\ \end{array} \right) p = {R\over c_v} \left(E-\half \rho\left(u_1^2 + u_2^2 + u_3^3\right)\right) = {R\over c_v} \left(w_4-{w_1^2+w_2^2+w_3^2\over2w_0}\right)

We solve for the unknowns w_0, w_1, w_2, w_3 and w_4 as functions of (t, x, y, z), the rest (R, c_v, f_x, f_y, f_z) are either constants or depend on the unknowns. In order to convert from the physical quantities \rho, u_1, u_2, u_3, E and p to w_0, …, w_4, we use:

w_0 = \rho w_1 = \rho u_1 w_2 = \rho u_2 w_3 = \rho u_3 w_4 = E = p {c_v \over R} + \half \rho \left(u_1^2 + u_2^2 + u_3^2\right)

the opposite conversion is:

\rho = w_0 u_1 = {w_1\over w_0} u_2 = {w_2\over w_0} u_3 = {w_3\over w_0} E = w_4 p = {R\over c_v} \left(w_4-{w_1^2+w_2^2+w_3^2\over2w_0}\right)

Sometimes people also use u, v and w instead of u_1, u_2 and u_3.

Note: \rho {\bf u}\equiv{\bf j}, where {\bf j} is the fluid density current (it’s a 3-vector) and also w^\mu \equiv j^\mu (here w^\mu is the same as w_\mu, e.g. we are a bit sloppy about the notation), where j^\mu is the density 4-current (e.g. the first 4 components of {\bf w} are exactly the components of the 4-current j^\mu):

j^\mu =\rho v^\mu = \rho\gamma (c, {\bf u}) = \gamma \left( \begin{array}{c} c\rho\\ \rho u_1\\ \rho u_2\\ \rho u_3\\ \end{array} \right)

where as usual \mu = 0, 1, 2, 3 is the relativistic index, c is the speed of light, and in the nonrelativistic limit (c\to\infty) we get \gamma\to1 and the remaining c in j^0 will cancel with c in \partial_0 = {1\over c}{\partial\over\partial t}, so it will not be present in the final equations (that involve terms like \partial_\mu j^\mu). We can also just set c=1 as usual in relativistic physics.

7.3.4. Weak Formulation

The Euler equations:

{\partial{\bf w}\over \partial t} + {\partial{\bf f}_x\over \partial x} + {\partial{\bf f}_y\over \partial y} + {\partial{\bf f}_z\over \partial z} + {\bf g}= 0

are nonlinear. The time-derivative is approximated using the implicit Euler method

{{\bf w}^{n+1}-{\bf w}^n\over \tau} + {\partial{\bf f}_x({\bf w}^{n+1})\over \partial x} + {\partial{\bf f}_y({\bf w}^{n+1})\over \partial y} + {\partial{\bf f}_z({\bf w}^{n+1})\over \partial z} + {\bf g}= 0

The vector-valued test functions for the above system of equations have the form:

\left( \begin{array}{c} \varphi^0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ \varphi^1 \\ 0 \\ 0 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ 0 \\ \varphi^2 \\ 0 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \varphi^3 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ \varphi^4 \\ \end{array} \right)

After multiplying the equation system with the test functions and integrating over the domain \Omega, we obtain (here the index i=0, 1, 2, 3, 4 is numbering the 5 equations, so we are not summing over it):

\int_{\Omega} {w_i^{n+1}-w_i^n\over\tau}\varphi^i +{\partial\left({\bf f}_x({\bf w}^{n+1})\right)_i\over \partial x}\varphi^i +{\partial\left({\bf f}_y({\bf w}^{n+1})\right)_i\over \partial y}\varphi^i +{\partial\left({\bf f}_z({\bf w}^{n+1})\right)_i\over \partial z}\varphi^i + g_i \varphi^i \ \d^3 x =0

Now we integrate by parts:

\int_{\Omega} {w_i^{n+1}-w_i^n\over\tau}\varphi^i - \left({\bf f}_x({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial x} - \left({\bf f}_y({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial y} - \left({\bf f}_z({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial z} + g_i \varphi^i \ \d^3 x + +\int_{\partial\Omega} \left({\bf f}_x({\bf w}^{n+1})\right)_i \varphi^i\, n_x + \left({\bf f}_y({\bf w}^{n+1})\right)_i \varphi^i\, n_y + \left({\bf f}_z({\bf w}^{n+1})\right)_i \varphi^i\, n_z \ \d^2 x =0

where {\bf n} = (n_x, n_y, n_z) is the outward surface normal to \partial\Omega. Rearranging:

\int_{\Omega} {w_i^{n+1}\over\tau}\varphi^i - \left({\bf f}_x({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial x} - \left({\bf f}_y({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial y} - \left({\bf f}_z({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial z} \ \d^3 x + +\int_{\partial\Omega} \left({\bf f}_x({\bf w}^{n+1})\right)_i \varphi^i\, n_x + \left({\bf f}_y({\bf w}^{n+1})\right)_i \varphi^i\, n_y + \left({\bf f}_z({\bf w}^{n+1})\right)_i \varphi^i\, n_z \ \d^2 x = \int_{\Omega} {w_i^n\over\tau}\varphi^i - g_i\varphi^i \ \d^3 x

We can then linearize this for example by taking the flux jacobians {\bf A}_x({\bf w}^{n+1}) on the previous time level {\bf A}_x({\bf w}^n).

The finite element formulation is obtained from here by replacing in the standard way the unknown solution w^{n+1} by a piecewise-polynomial unknown function

w_h^{n+1} = \sum_{k=1}^N y_k \psi_k,

where \psi_k are the basis functions of the piecewise-polynomial finite element space. This turns the above weak formulation into a finite number of nonlinear algebraic equations of the form F(Y) = 0 that will be solved using the Newton’s method.

Explicit Method

We also derive the weak formulation for the explicit method. Euler equations:

{\partial{\bf w}\over \partial t} + {\partial{\bf f}_x\over \partial x} + {\partial{\bf f}_y\over \partial y} + {\partial{\bf f}_z\over \partial z} + {\bf g}= 0

The time-derivative is approximated using the explicit Euler method

{{\bf w}^{n+1}-{\bf w}^n\over \tau} + {\partial{\bf f}_x({\bf w}^{n})\over \partial x} + {\partial{\bf f}_y({\bf w}^{n})\over \partial y} + {\partial{\bf f}_z({\bf w}^{n})\over \partial z} + {\bf g}= 0

The vector-valued test functions for the above system of equations have the form:

\left( \begin{array}{c} \varphi^0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ \varphi^1 \\ 0 \\ 0 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ 0 \\ \varphi^2 \\ 0 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \varphi^3 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ \varphi^4 \\ \end{array} \right)

After multiplying the equation system with the test functions and integrating over the domain \Omega, we obtain (here the index i=0, 1, 2, 3, 4 is numbering the 5 equations, so we are not summing over it):

\int_{\Omega} {w_i^{n+1}-w_i^n\over\tau}\varphi^i +{\partial\left({\bf f}_x({\bf w}^{n})\right)_i\over \partial x}\varphi^i +{\partial\left({\bf f}_y({\bf w}^{n})\right)_i\over \partial y}\varphi^i +{\partial\left({\bf f}_z({\bf w}^{n})\right)_i\over \partial z}\varphi^i + g_i \varphi^i \ \d^3 x =0

Now we integrate by parts:

\int_{\Omega} {w_i^{n+1}-w_i^n\over\tau}\varphi^i - \left({\bf f}_x({\bf w}^{n})\right)_i {\partial \varphi^i\over\partial x} - \left({\bf f}_y({\bf w}^{n})\right)_i {\partial \varphi^i\over\partial y} - \left({\bf f}_z({\bf w}^{n})\right)_i {\partial \varphi^i\over\partial z} + g_i \varphi^i \ \d^3 x + +\int_{\partial\Omega} \left({\bf f}_x({\bf w}^{n})\right)_i \varphi^i\, n_x + \left({\bf f}_y({\bf w}^{n})\right)_i \varphi^i\, n_y + \left({\bf f}_z({\bf w}^{n})\right)_i \varphi^i\, n_z \ \d^2 x =0

where {\bf n} = (n_x, n_y, n_z) is the outward surface normal to \partial\Omega. Rearranging:

\int_{\Omega} {w_i^{n+1}\over\tau}\varphi^i \ \d^3 x = \int_{\Omega}{w_i^n\over\tau}\varphi^i + \left({\bf f}_x({\bf w}^{n})\right)_i {\partial \varphi^i\over\partial x} + \left({\bf f}_y({\bf w}^{n})\right)_i {\partial \varphi^i\over\partial y} + \left({\bf f}_z({\bf w}^{n})\right)_i {\partial \varphi^i\over\partial z} - g_i\varphi^i \ \d^3 x + -\int_{\partial\Omega} \left({\bf f}_x({\bf w}^{n})\right)_i \varphi^i\, n_x + \left({\bf f}_y({\bf w}^{n})\right)_i \varphi^i\, n_y + \left({\bf f}_z({\bf w}^{n})\right)_i \varphi^i\, n_z \ \d^2 x

7.3.5. Flux Jacobians

Now we write the spatial derivatives using the so called flux Jacobians {\bf A}_x, {\bf A}_y and {\bf A}_z:

{\partial{\bf f}_x\over \partial x} = {\partial{\bf f}_x\over \partial {\bf w}} {\partial{\bf w}\over \partial x} \equiv {\bf A}_x {\partial{\bf w}\over \partial x} {\bf A}_x={\bf A}_x({\bf w})\equiv{\partial{\bf f}_x\over \partial {\bf w}}

Similarly for y and z, so we get:

{\partial{\bf w}\over \partial t} + {\bf A}_x {\partial{\bf w}\over \partial x} + {\bf A}_y {\partial{\bf w}\over \partial y} + {\bf A}_z {\partial{\bf w}\over \partial z} + {\bf g}= 0

One nice thing about these particular {\bf f}_x, {\bf f}_y and {\bf f}_z functions is that they are homogeneous of degree 1:

{\bf f}_x(\lambda{\bf w}) =\lambda\,{\bf f}_x({\bf w})

so the Euler equation/formula for the homogeneous function is:

{\bf w}\cdot {\partial {\bf f}_x({\bf w})\over\partial {\bf w}} ={\bf f}_x({\bf w}) {\bf w}\cdot {\bf A}_x ={\bf f}_x({\bf w})

So both the {\bf f}_x and it’s derivative can be nicely factored out using the flux Jacobian:

{\bf f}_x = {\bf A}_x\, {\bf w} {\partial{\bf f}_x\over \partial x} = {\bf A}_x {\partial{\bf w}\over \partial x}

by differentiating the first equation and substracting the second, we get:

{\partial {\bf A}_x\over\partial x}\, {\bf w} = 0

similarly for y and z. To calculate the Jacobians, we’ll need:

{\partial p\over \partial {\bf w}}= {R\over c_v} \left( \begin{array}{ccccc} {w_1^2+w_2^2+w_3^2\over 2w_0^2} & -{w_1\over w_0} & -{w_2\over w_0} & -{w_3\over w_0} & 1\\ \end{array} \right)

then we can calculate the Jacobians (and we substitute for p):

{\bf A}_x({\bf w}) = {\partial{\bf f}_x\over \partial {\bf w}}= \left( \begin{array}{ccccc} 0 & 1 & 0 & 0 & 0\\ -{w_1^2\over w_0^2} +{R\over c_v}{w_1^2+w_2^2+w_3^2\over 2 w_0^2} & {2w_1\over w_0}-{R\over c_v}{w_1\over w_0} & -{R\over c_v}{w_2\over w_0} & -{R\over c_v}{w_3\over w_0} & {R\over c_v}\\ -{w_1w_2\over w_0^2} & {w_2\over w_0} & {w_1\over w_0} & 0 & 0\\ -{w_1w_3\over w_0^2} & {w_3\over w_0} & 0 & {w_1\over w_0} & 0 \\ -{w_1w_4\over w_0^2}-{w_1\over w_0^2}{R\over c_v} \left(w_4-{w_1^2+w_2^2+w_3^2\over 2 w_0}\right) +{w_1\over w_0}{R\over c_v}{w_1^2+w_2^2+w_3^2\over 2 w_0^2}& {w_4\over w_0}+{1\over w_0}{R\over c_v} \left(w_4-{w_1^2+w_2^2+w_3^2\over 2 w_0}\right) -{R\over c_v}{w_1^2\over w_0^2} & -{R\over c_v}{w_1w_2\over w_0^2} & -{R\over c_v}{w_1w_3\over w_0^2} & {w_1\over w_0}+{R\over c_v}{w_1\over w_0} \\ \end{array} \right) {\bf A}_y({\bf w}) = {\partial{\bf f}_y\over \partial {\bf w}}= \left( \begin{array}{ccccc} 0 & 0 & 1 & 0 & 0\\ -{w_2w_1\over w_0^2} & {w_2\over w_0} & {w_1\over w_0} & 0 & 0\\ -{w_2^2\over w_0^2} +{R\over c_v}{w_1^2+w_2^2+w_3^2\over 2 w_0^2} & -{R\over c_v}{w_1\over w_0} & {2w_2\over w_0}-{R\over c_v}{w_2\over w_0} & -{R\over c_v}{w_3\over w_0} & {R\over c_v}\\ -{w_2w_3\over w_0^2} & 0 & {w_3\over w_0} & {w_2\over w_0} & 0 \\ -{w_2w_4\over w_0^2}-{w_2\over w_0^2}{R\over c_v} \left(w_4-{w_1^2+w_2^2+w_3^2\over 2 w_0}\right) +{w_2\over w_0}{R\over c_v}{w_1^2+w_2^2+w_3^2\over 2 w_0^2}& -{R\over c_v}{w_2w_1\over w_0^2} & {w_4\over w_0}+{1\over w_0}{R\over c_v} \left(w_4-{w_1^2+w_2^2+w_3^2\over 2 w_0}\right) -{R\over c_v}{w_2^2\over w_0^2} & -{R\over c_v}{w_2w_3\over w_0^2} & {w_2\over w_0}+{R\over c_v}{w_2\over w_0} \\ \end{array} \right) {\bf A}_z({\bf w}) = {\partial{\bf f}_z\over \partial {\bf w}}= \left( \begin{array}{ccccc} 0 & 0 & 0 & 1 & 0\\ -{w_3w_1\over w_0^2} & {w_3\over w_0} & 0 & {w_1\over w_0} & 0 \\ -{w_3w_2\over w_0^2} & 0 & {w_3\over w_0} & {w_2\over w_0} & 0 \\ -{w_3^2\over w_0^2} +{R\over c_v}{w_1^2+w_2^2+w_3^2\over 2 w_0^2} & -{R\over c_v}{w_1\over w_0} & -{R\over c_v}{w_2\over w_0} & {2w_3\over w_0} -{R\over c_v}{w_3\over w_0} & {R\over c_v}\\ -{w_3w_4\over w_0^2}-{w_3\over w_0^2}{R\over c_v} \left(w_4-{w_1^2+w_2^2+w_3^2\over 2 w_0}\right) +{w_3\over w_0}{R\over c_v}{w_1^2+w_2^2+w_3^2\over 2 w_0^2}& -{R\over c_v}{w_3w_1\over w_0^2} & -{R\over c_v}{w_3w_2\over w_0^2} & {w_4\over w_0}+{1\over w_0}{R\over c_v} \left(w_4-{w_1^2+w_2^2+w_3^2\over 2 w_0}\right) -{R\over c_v}{w_3^2\over w_0^2} & {w_3\over w_0}+{R\over c_v}{w_3\over w_0} \\ \end{array} \right)

7.3.6. 2D Version of the Equations

{\partial{\bf w}\over \partial t} + {\partial{\bf f}_x\over \partial x} + {\partial{\bf f}_y\over \partial y} + {\bf g}= 0

where:

{\bf w} = \left( \begin{array}{c} \varrho\\ \rho u_1\\ \rho u_2\\ E \end{array} \right) = \left( \begin{array}{c} w_0 \\ w_1 \\ w_2 \\ w_3 \\ \end{array} \right) {\bf f}_x = \left( \begin{array}{c} \rho u_1\\ \rho u_1^2 + p\\ \rho u_1 u_2\\ u_1(E+p) \end{array} \right) = \left( \begin{array}{c} w_1\\ \frac{w_1^2}{w_0} + p\\ \frac{w_1w_2}{w_0}\\ \frac{w_1}{w_0}(w_3+p) \end{array} \right) {\bf f}_y = \left( \begin{array}{c} \rho u_2\\ \rho u_2 u_1\\ \rho u_2^2 + p\\ u_2(E+p) \end{array} \right) = \left( \begin{array}{c} w_2\\ \frac{w_2w_1}{w_0}\\ \frac{w_2^2}{w_0} + p\\ \frac{w_2}{w_0}(w_3+p) \end{array} \right) {\bf g} = \left( \begin{array}{c} 0\\ -f_x\\ -f_y\\ 0\\ \end{array} \right) p = {R\over c_v} \left(E-\half \rho\left(u_1^2 + u_2^2\right)\right) = {R\over c_v} \left(w_3-{w_1^2+w_2^2\over2w_0}\right)

Discretizing the time derivative:

{{\bf w}^{n+1}-{\bf w}^n\over \tau} + {\partial{\bf f}_x({\bf w}^{n+1})\over \partial x} + {\partial{\bf f}_y({\bf w}^{n+1})\over \partial y} + {\bf g}= 0

The vector-valued test functions for the above system of equations have the form:

\left( \begin{array}{c} \varphi^0 \\ 0 \\ 0 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ \varphi^1 \\ 0 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ 0 \\ \varphi^2 \\ 0 \\ \end{array} \right),\ \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \varphi^3 \\ \end{array} \right)

After multiplying the equation system with the test functions and integrating over the domain \Omega, we obtain:

\int_{\Omega} {w_i^{n+1}-w_i^n\over\tau}\varphi^i +{\partial\left({\bf f}_x({\bf w}^{n+1})\right)_i\over \partial x}\varphi^i +{\partial\left({\bf f}_y({\bf w}^{n+1})\right)_i\over \partial y}\varphi^i + g_i \varphi^i \ \d^2 x =0

Now we integrate by parts:

\int_{\Omega} {w_i^{n+1}-w_i^n\over\tau}\varphi^i - \left({\bf f}_x({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial x} - \left({\bf f}_y({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial y} + g_i \varphi^i \ \d^2 x + +\int_{\partial\Omega} \left({\bf f}_x({\bf w}^{n+1})\right)_i \varphi^i\, n_x + \left({\bf f}_y({\bf w}^{n+1})\right)_i \varphi^i\, n_y \ \d x =0

where {\bf n} = (n_x, n_y) is the outward surface normal to \partial\Omega. Rearranging:

\int_{\Omega} {w_i^{n+1}\over\tau}\varphi^i - \left({\bf f}_x({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial x} - \left({\bf f}_y({\bf w}^{n+1})\right)_i {\partial \varphi^i\over\partial y} \ \d^2 x + +\int_{\partial\Omega} \left({\bf f}_x({\bf w}^{n+1})\right)_i \varphi^i\, n_x + \left({\bf f}_y({\bf w}^{n+1})\right)_i \varphi^i\, n_y \ \d x = \int_{\Omega} {w_i^n\over\tau}\varphi^i - g_i\varphi^i \ \d^2 x

The 2D flux Jacobians are:

{\bf A}_x({\bf w}) = {\partial{\bf f}_x\over \partial {\bf w}}= \left( \begin{array}{cccc} 0 & 1 & 0 & 0\\ -{w_1^2\over w_0^2} +{R\over c_v}{w_1^2+w_2^2\over 2 w_0^2} & {2w_1\over w_0}-{R\over c_v}{w_1\over w_0} & -{R\over c_v}{w_2\over w_0} & {R\over c_v}\\ -{w_1w_2\over w_0^2} & {w_2\over w_0} & {w_1\over w_0} & 0 \\ -{w_1w_3\over w_0^2}-{w_1\over w_0^2}{R\over c_v} \left(w_3-{w_1^2+w_2^2\over 2 w_0}\right) +{w_1\over w_0}{R\over c_v}{w_1^2+w_2^2\over 2 w_0^2} & {w_3\over w_0}+{1\over w_0}{R\over c_v} \left(w_3-{w_1^2+w_2^2\over 2 w_0}\right) -{R\over c_v}{w_1^2\over w_0^2} & -{R\over c_v}{w_1w_2\over w_0^2} & {w_1\over w_0}+{R\over c_v}{w_1\over w_0} \\ \end{array} \right) {\bf A}_y({\bf w}) = {\partial{\bf f}_y\over \partial {\bf w}}= \left( \begin{array}{cccc} 0 & 0 & 1 & 0\\ -{w_2w_1\over w_0^2} & {w_2\over w_0} & {w_1\over w_0} & 0 \\ -{w_2^2\over w_0^2} +{R\over c_v}{w_1^2+w_2^2\over 2 w_0^2} & -{R\over c_v}{w_1\over w_0} & {2w_2\over w_0} -{R\over c_v}{w_2\over w_0} & {R\over c_v}\\ -{w_2w_3\over w_0^2}-{w_2\over w_0^2}{R\over c_v} \left(w_3-{w_1^2+w_2^2\over 2 w_0}\right) +{w_2\over w_0}{R\over c_v}{w_1^2+w_2^2\over 2 w_0^2}& -{R\over c_v}{w_2w_1\over w_0^2} & {w_3\over w_0}+{1\over w_0}{R\over c_v} \left(w_3-{w_1^2+w_2^2\over 2 w_0}\right) -{R\over c_v}{w_2^2\over w_0^2} & {w_2\over w_0}+{R\over c_v}{w_2\over w_0} \\ \end{array} \right)

7.3.7. Sea Breeze Modeling

In our 2D model we make the following assumptions:

f_x = 0 f_y = -\rho g = -w_0 g

and the boundary condition is as follows:

T'(x, t) = \left(A\over2\right) \sin \left(\pi (t-t_0)\over 24\right) \left(1+\tanh\left(S(x)\over L\right)\right) T(x) = T_0 + T'(x, t)

The weak formulation in 2D is (here i = 0, 1, 2, 3):

\int_{\Omega} {w_i^{n+1}\over\tau}\varphi^i - \left({\bf A}_x({\bf w}^n)\right)_{ij} w_j^{n+1} {\partial \varphi^i\over\partial x} - \left({\bf A}_z({\bf w}^n)\right)_{ij} w_j^{n+1} {\partial \varphi^i\over\partial z} \ \d^2 x + +\int_{\partial\Omega} \left({\bf A}_x({\bf w}^n)\right)_{ij}w_j^{n+1} \varphi^i\, n_x + \left({\bf A}_z({\bf w}^n)\right)_{ij}w_j^{n+1} \varphi^i\, n_z \ \d x = \int_{\Omega} {w_i^n\over\tau}\varphi^i - g_i \varphi^i \ \d^2 x

In order to specify the input forms for Hermes, we’ll write the weak formulation as:

B_{00}(w_0, \varphi^0) + B_{01}(w_1, \varphi^0) + B_{02}(w_2, \varphi^0)+ B_{03}(w_3, \varphi^0) = l_0(\varphi^0) B_{10}(w_0, \varphi^1) + B_{11}(w_1, \varphi^1) + B_{12}(w_2, \varphi^1)+ B_{13}(w_3, \varphi^1) = l_1(\varphi^1) B_{20}(w_0, \varphi^2) + B_{21}(w_1, \varphi^2) + B_{22}(w_2, \varphi^2)+ B_{23}(w_3, \varphi^2) = l_2(\varphi^2) B_{30}(w_0, \varphi^3) + B_{31}(w_1, \varphi^3) + B_{32}(w_2, \varphi^3)+ B_{33}(w_3, \varphi^3) = l_3(\varphi^3)

where the forms are (we write w_i instead of w_i^{n+1}):

l_0(\varphi^0) = \int_\Omega {w_0^n\varphi^0\over\tau} \,\d^2 x l_1(\varphi^1) = \int_\Omega {w_1^n\varphi^1\over\tau} \,\d^2 x l_2(\varphi^2) = \int_\Omega {w_2^n\varphi^2\over\tau} + \rho g \varphi^2 \,\d^2 x l_3(\varphi^3) = \int_\Omega {w_3^n\varphi^3\over\tau} \,\d^2 x B_{ij}(w_j, \varphi^i) = \int_{\Omega} {w_i\over\tau}\varphi^i \delta_{ij} - \left({\bf A}_x({\bf w}^n)\right)_{ij} w_j {\partial \varphi^i\over\partial x} - \left({\bf A}_z({\bf w}^n)\right)_{ij} w_j {\partial \varphi^i\over\partial z} \ \d^2 x

In the last expression we do not sum over i nor j. In particular:

B_{00}(w_0, \varphi^0) = \int_{\Omega} {w_0\over\tau}\varphi^0 - \left({\bf A}_x({\bf w}^n)\right)_{00} w_0 {\partial \varphi^0\over\partial x} - \left({\bf A}_z({\bf w}^n)\right)_{00} w_0 {\partial \varphi^0\over\partial z} \ \d^2 x = \int_{\Omega} {w_0\over\tau}\varphi^0 \ \d^2 x B_{01}(w_1, \varphi^0) = \int_{\Omega} - \left({\bf A}_x({\bf w}^n)\right)_{01} w_1 {\partial \varphi^0\over\partial x} - \left({\bf A}_z({\bf w}^n)\right)_{01} w_1 {\partial \varphi^0\over\partial z} \ \d^2 x = \int_{\Omega} - \left({\bf A}_x({\bf w}^n)\right)_{01} w_1 {\partial \varphi^0\over\partial x} \ \d^2 x B_{02}(w_2, \varphi^0) = \int_{\Omega} - \left({\bf A}_x({\bf w}^n)\right)_{02} w_2 {\partial \varphi^0\over\partial x} - \left({\bf A}_z({\bf w}^n)\right)_{02} w_2 {\partial \varphi^0\over\partial z} \ \d^2 x = \int_{\Omega} - \left({\bf A}_z({\bf w}^n)\right)_{02} w_2 {\partial \varphi^0\over\partial z} \ \d^2 x B_{03}(w_3, \varphi^0) = \int_{\Omega} - \left({\bf A}_x({\bf w}^n)\right)_{03} w_3 {\partial \varphi^0\over\partial x} - \left({\bf A}_z({\bf w}^n)\right)_{03} w_3 {\partial \varphi^0\over\partial z} \ \d^2 x =0 B_{10}(w_0, \varphi^1) = \int_{\Omega} - \left({\bf A}_x({\bf w}^n)\right)_{10} w_0 {\partial \varphi^1\over\partial x} - \left({\bf A}_z({\bf w}^n)\right)_{10} w_0 {\partial \varphi^1\over\partial z} \ \d^2 x B_{11}(w_1, \varphi^1) = \int_{\Omega} {w_1\over\tau}\varphi^1 - \left({\bf A}_x({\bf w}^n)\right)_{11} w_1 {\partial \varphi^1\over\partial x} - \left({\bf A}_z({\bf w}^n)\right)_{11} w_1 {\partial \varphi^1\over\partial z} \ \d^2 x \cdots

7.3.8. Boundary Conditions

We rewrite the boundary integral by rotating coordinates, so that the flow is only in the x direction (thus we only have {\bf f}_x):

\int_{\partial\Omega} \left({\bf f}_x({\bf w})\right)_i \varphi^i\, n_x + \left({\bf f}_y({\bf w})\right)_i \varphi^i\, n_y + \left({\bf f}_z({\bf w})\right)_i \varphi^i\, n_z \ \d^2 x = = \int_{\partial\Omega} T^{-1} {\bf f}_x(T {\bf w}) \varphi^i \ \d^2 x

Now we need to approximate {\bf f}_x(T {\bf w}) somehow. We do that by solving the following 1D problem (Riemann problem):

{\partial {\bf w}\over\partial t} + {\partial\over\partial x} {\bf f}({\bf w}) = 0

or:

(7.3.8.1){\partial {\bf w}\over\partial t} + {\bf A}({\bf w}) {\partial{\bf w}\over\partial x} = 0

{\bf w}(x, t) = \left( \begin{array}{c} w_0 \\ w_1 \\ w_2 \\ w_3 \\ w_4 \\ \end{array} \right)

And we approximate {\bf f}_x({\bf w})={\bf f}({\bf w}(0, t)). The initial condition is:

{\bf w}(x, 0) = \begin{cases}{\bf w}_L&x<0\cr {\bf w}_R & x > 0\cr \end{cases} = {\bf w}_L(1-H(x)) + {\bf w}_R H(x)

Now we write:

{\bf w}(x, t) = \sum_i \xi^i(x, t) {\bf r}_i {\bf w}_L = \sum_i \alpha_i {\bf r}_i {\bf w}_R = \sum_i \beta_i {\bf r}_i \xi^i(x, 0) = \begin{cases} \alpha_i & x < 0\cr \beta_i & x > 0\cr \end{cases}

and substitute into (7.3.8.1):

\sum_i\left( {\partial \xi^i\over\partial t} + {\bf A}({\bf w}) {\partial \xi^i\over\partial x} \right) {\bf r}_i = 0 \sum_i\left( {\partial \xi^i\over\partial t} + \lambda_i({\bf w}) {\partial \xi^i\over\partial x} \right) {\bf r}_i = 0

so we get:

{\partial \xi^i\over\partial t} + \lambda_i({\bf w}) {\partial \xi^i\over\partial x} = 0

This is a nonlinear problem, that cannot be solved exactly. First, let {\bf A} doesn’t depend on {\bf w}. Then also \lambda_i are constants:

{\partial \xi^i\over\partial t} + \lambda_i {\partial \xi^i\over\partial x} = 0

and the solution is constant along the characteristic x(t) = \lambda_i t + c for t>0 and we get:

\xi_i(x, t) = \xi^i(x-\lambda_i t, 0) = \begin{cases} \alpha_i & x-\lambda_i t < 0\cr \beta_i & x-\lambda_i t > 0\cr \end{cases} =\alpha_i (1-H(x-\lambda_i t)) + \beta_i H(x-\lambda_i t)

and

{\bf w}(x, t) = \sum_i \xi^i(x, t) {\bf r}_i = \sum_i \left( \alpha_i (1-H(x-\lambda_i t)) + \beta_i H(x-\lambda_i t) \right){\bf r}_i {\bf w}(0, t) = \sum_i \left( \alpha_i (1-H(-\lambda_i t)) + \beta_i H(-\lambda_i t) \right){\bf r}_i = = \sum_i \left( \alpha_i H(\lambda_i t) + \beta_i H(-\lambda_i t) \right){\bf r}_i= = \sum_i \left( \alpha_i H(\lambda_i) + \beta_i H(-\lambda_i) \right){\bf r}_i= = \sum_{i=k+1}^n \alpha_i {\bf r}_i + \sum_{i=1}^k \beta_i {\bf r}_i

so:

{\bf f}({\bf w}(0, t)) = {\bf A}{\bf w}(0, t) = \sum_{i=k+1}^n {\bf A}\alpha_i {\bf r}_i + \sum_{i=1}^k {\bf A}\beta_i {\bf r}_i = \sum_{i=k+1}^n \lambda_i\alpha_i {\bf r}_i + \sum_{i=1}^k \lambda_i\beta_i {\bf r}_i= = {\bf A}^+\sum_{i=1}^n \alpha_i {\bf r}_i + {\bf A}^-\sum_{i=1}^n \beta_i {\bf r}_i= = {\bf A}^+{\bf w}_L + {\bf A}^-{\bf w}_R

In the nonlinear case we cannot solve it exactly, but we can approximate the solution by:

{\bf f}({\bf w}(0, t)) = {\bf f}^+({\bf w}_L) + {\bf f}^-({\bf w}_R) = = {\bf f}({\bf w}_R) - \int_{{\bf w}_L}^{{\bf w}_R} {\bf A}^+({\bf w}) \d {\bf w} = = {\bf f}({\bf w}_L) + \int_{{\bf w}_L}^{{\bf w}_R} {\bf A}^-({\bf w}) \d {\bf w} \approx

(7.3.8.2)\approx {\bf f}({\bf w}_L) + {\bf A}^-({\bf w}_R) {\bf w}_R - {\bf A}^-({\bf w}_L) {\bf w}_L

Let’s say the domain is for x<0 and we are applying the BC condition from x>0. Then {\bf w}_L is taken from the solution and {\bf w}_R is prescribed, for example at the bottom it could be:

{\bf w}_R = \left( \begin{array}{c} \rho \\ \rho u_1 \\ 0 \\ 0 \\ E \\ \end{array} \right)

Now we need to calculate {\bf A}^-. First we write:

{\bf A}_x = {\bf R}{\bf D}_x{\bf R}^{-1} {\bf A}_x^- = {\bf R}{\bf D}_x^-{\bf R}^{-1} {\bf D}_x({\bf w}) = {w_1\over w_0}\one + \diag(-c, 0, 0, 0, c) = \left( \begin{array}{ccccc} u_1-c & 0 & 0 & 0 & 0 \\ 0 & u_1 & 0 & 0 & 0 \\ 0 & 0 & u_1 & 0 & 0 \\ 0 & 0 & 0 & u_1 & 0 \\ 0 & 0 & 0 & 0 & u_1 + c \\ \end{array} \right) {\bf D}_x({\bf w})^- =\begin{cases} \diag({w_1\over w_0}-c, {w_1\over w_0}, {w_1\over w_0}, {w_1\over w_0}, 0) & w_1 < 0\cr \diag({w_1\over w_0}-c, 0, 0, 0, 0) & w_1 > 0\cr \end{cases}

Explicit forms of the matrices:

{\bf R} = \left(\begin{smallmatrix}1 & 1 & 1 & 1 & 1\\u - c & u & u & u & c + u\\v & v & v - c & v & v\\w & w & w & w - c & w\\- c u - \frac{c^{2}}{1 - \kappa} + \frac{1}{2} u^{2} + \frac{1}{2} v^{2} + \frac{1}{2} w^{2} & \frac{1}{2} u^{2} + \frac{1}{2} v^{2} + \frac{1}{2} w^{2} & - c v + \frac{1}{2} u^{2} + \frac{1}{2} v^{2} + \frac{1}{2} w^{2} & - c w + \frac{1}{2} u^{2} + \frac{1}{2} v^{2} + \frac{1}{2} w^{2} & c u - \frac{c^{2}}{1 - \kappa} + \frac{1}{2} u^{2} + \frac{1}{2} v^{2} + \frac{1}{2} w^{2}\end{smallmatrix}\right) {\bf R}^{-1} = {1\over c^2} \left(\begin{smallmatrix}\frac{1}{2} c u - \frac{1}{4} u^{2} - \frac{1}{4} v^{2} - \frac{1}{4} w^{2} + \frac{1}{4} \kappa u^{2} + \frac{1}{4} \kappa v^{2} + \frac{1}{4} \kappa w^{2} & \frac{1}{2} u - \frac{1}{2} c - \frac{1}{2} \kappa u & \frac{1}{2} v - \frac{1}{2} \kappa v & \frac{1}{2} w - \frac{1}{2} \kappa w & - \frac{1}{2} + \frac{1}{2} \kappa\\- c v - c w + c^{2} + \frac{1}{2} u^{2} + \frac{1}{2} v^{2} + \frac{1}{2} w^{2} - \frac{1}{2} \kappa u^{2} - \frac{1}{2} \kappa v^{2} - \frac{1}{2} \kappa w^{2} & - u + \kappa u & c - v + \kappa v & c - w + \kappa w & 1 - \kappa\\c v & 0 & - c & 0 & 0\\c w & 0 & 0 & - c & 0\\- \frac{1}{2} c u - \frac{1}{4} u^{2} - \frac{1}{4} v^{2} - \frac{1}{4} w^{2} + \frac{1}{4} \kappa u^{2} + \frac{1}{4} \kappa v^{2} + \frac{1}{4} \kappa w^{2} & \frac{1}{2} c + \frac{1}{2} u - \frac{1}{2} \kappa u & \frac{1}{2} v - \frac{1}{2} \kappa v & \frac{1}{2} w - \frac{1}{2} \kappa w & - \frac{1}{2} + \frac{1}{2} \kappa\end{smallmatrix}\right) {\bf D}_x = \left(\begin{smallmatrix}u - c & 0 & 0 & 0 & 0\\0 & u & 0 & 0 & 0\\0 & 0 & u & 0 & 0\\0 & 0 & 0 & u & 0\\0 & 0 & 0 & 0 & c + u\end{smallmatrix}\right) {\bf A}_x = \left(\begin{smallmatrix}0 & 1 & 0 & 0 & 0\\- \frac{3}{2} u^{2} - \frac{1}{2} v^{2} - \frac{1}{2} w^{2} + \frac{1}{2} \kappa u^{2} + \frac{1}{2} \kappa v^{2} + \frac{1}{2} \kappa w^{2} & 3 u - \kappa u & v - \kappa v & w - \kappa w & -1 + \kappa\\- u v & v & u & 0 & 0\\- u w & w & 0 & u & 0\\\frac{- 2 u v^{2} - 2 u w^{2} + 2 u c^{2} - u \kappa^{2} v^{2} - u \kappa^{2} w^{2} + 3 \kappa u v^{2} + 3 \kappa u w^{2} - 2 u^{3} - \kappa^{2} u^{3} + 3 \kappa u^{3}}{2 - 2 \kappa} & \frac{v^{2} + w^{2} - 2 c^{2} + 3 u^{2} - \kappa v^{2} - \kappa w^{2} - 5 \kappa u^{2} + 2 \kappa^{2} u^{2}}{2 - 2 \kappa} & u v - \kappa u v & u w - \kappa u w & \kappa u\end{smallmatrix}\right)

For u_1<0:

{\bf D}_x^- = \left(\begin{smallmatrix}u - c & 0 & 0 & 0 & 0\\0 & u & 0 & 0 & 0\\0 & 0 & u & 0 & 0\\0 & 0 & 0 & u & 0\\0 & 0 & 0 & 0 & 0\end{smallmatrix}\right) {\bf A}_x^- = \left(\begin{smallmatrix}\frac{2 c v^{2} + 2 c w^{2} + 2 u v^{2} + 2 u w^{2} + 4 u c^{2} + 6 c u^{2} - 2 c \kappa u^{2} - 2 c \kappa v^{2} - 2 c \kappa w^{2} - 2 \kappa u v^{2} - 2 \kappa u w^{2} + 2 u^{3} - 2 \kappa u^{3}}{8 c^{2}} & \frac{- 2 c u + c \kappa u + c^{2} - u^{2} + \kappa u^{2}}{2 c^{2}} & \frac{- c v - u v + c \kappa v + \kappa u v}{2 c^{2}} & \frac{- c w - u w + c \kappa w + \kappa u w}{2 c^{2}} & \frac{c + u - c \kappa - \kappa u}{2 c^{2}}\\\frac{- 2 c^{2} u^{2} - 2 c^{2} v^{2} - 2 c^{2} w^{2} + 2 u^{2} v^{2} + 2 u^{2} w^{2} - 2 \kappa u^{2} v^{2} - 2 \kappa u^{2} w^{2} + 2 \kappa c^{2} u^{2} + 2 \kappa c^{2} v^{2} + 2 \kappa c^{2} w^{2} + 4 c u v^{2} + 4 c u w^{2} - 4 c \kappa u v^{2} - 4 c \kappa u w^{2} + 4 u c^{3} + 8 c u^{3} - 4 c \kappa u^{3} + 2 u^{4} - 2 \kappa u^{4}}{8 c^{2}} & \frac{- 3 c u^{2} + 3 u c^{2} - \kappa u c^{2} + 2 c \kappa u^{2} - c^{3} - u^{3} + \kappa u^{3}}{2 c^{2}} & \frac{- 2 c u v + 2 c \kappa u v + v c^{2} - v u^{2} + \kappa v u^{2} - \kappa v c^{2}}{2 c^{2}} & \frac{- 2 c u w + 2 c \kappa u w + w c^{2} - w u^{2} + \kappa w u^{2} - \kappa w c^{2}}{2 c^{2}} & \frac{2 c u - 2 c \kappa u + u^{2} - c^{2} + \kappa c^{2} - \kappa u^{2}}{2 c^{2}}\\\frac{- 4 u v c^{2} + 2 c v w^{2} + 2 u v w^{2} + 6 c v u^{2} - 2 c \kappa v u^{2} - 2 c \kappa v w^{2} - 2 \kappa u v w^{2} + 2 c v^{3} + 2 u v^{3} + 2 v u^{3} - 2 c \kappa v^{3} - 2 \kappa u v^{3} - 2 \kappa v u^{3}}{8 c^{2}} & \frac{- 2 c u v + c \kappa u v + v c^{2} - v u^{2} + \kappa v u^{2}}{2 c^{2}} & \frac{- c v^{2} - u v^{2} + 2 u c^{2} + c \kappa v^{2} + \kappa u v^{2}}{2 c^{2}} & \frac{- c v w - u v w + c \kappa v w + \kappa u v w}{2 c^{2}} & \frac{c v + u v - c \kappa v - \kappa u v}{2 c^{2}}\\\frac{- 4 u w c^{2} + 2 c w v^{2} + 2 u w v^{2} + 6 c w u^{2} - 2 c \kappa w u^{2} - 2 c \kappa w v^{2} - 2 \kappa u w v^{2} + 2 c w^{3} + 2 u w^{3} + 2 w u^{3} - 2 c \kappa w^{3} - 2 \kappa u w^{3} - 2 \kappa w u^{3}}{8 c^{2}} & \frac{- 2 c u w + c \kappa u w + w c^{2} - w u^{2} + \kappa w u^{2}}{2 c^{2}} & \frac{- c v w - u v w + c \kappa v w + \kappa u v w}{2 c^{2}} & \frac{- c w^{2} - u w^{2} + 2 u c^{2} + c \kappa w^{2} + \kappa u w^{2}}{2 c^{2}} & \frac{c w + u w - c \kappa w - \kappa u w}{2 c^{2}}\\\frac{- 2 c^{3} u^{2} - 2 c^{3} v^{2} - 2 c^{3} w^{2} + 2 u^{3} v^{2} + 2 u^{3} w^{2} - 6 u c^{2} v^{2} - 6 u c^{2} w^{2} - 4 \kappa u^{3} v^{2} - 4 \kappa u^{3} w^{2} - 2 \kappa c^{3} u^{2} + 2 c v^{2} w^{2} + 2 \kappa c^{3} v^{2} + 2 \kappa c^{3} w^{2} + 2 u v^{2} w^{2} + 2 \kappa^{2} u^{3} v^{2} + 2 \kappa^{2} u^{3} w^{2} + 6 c u^{2} v^{2} + 6 c u^{2} w^{2} - 10 c \kappa u^{2} v^{2} - 10 c \kappa u^{2} w^{2} - 4 c \kappa v^{2} w^{2} - 4 \kappa u v^{2} w^{2} - 2 u c^{2} \kappa^{2} v^{2} - 2 u c^{2} \kappa^{2} w^{2} + 2 c \kappa^{2} v^{2} w^{2} + 2 u \kappa^{2} v^{2} w^{2} + 4 c \kappa^{2} u^{2} v^{2} + 4 c \kappa^{2} u^{2} w^{2} + 8 \kappa u c^{2} v^{2} + 8 \kappa u c^{2} w^{2} - 2 c^{2} u^{3} - 2 c^{2} \kappa^{2} u^{3} + 4 \kappa c^{2} u^{3} + c v^{4} + c w^{4} + u v^{4} + u w^{4} + 4 u c^{4} + 5 c u^{4} + c \kappa^{2} v^{4} + c \kappa^{2} w^{4} + u \kappa^{2} v^{4} + u \kappa^{2} w^{4} - 8 c \kappa u^{4} - 2 c \kappa v^{4} - 2 c \kappa w^{4} - 2 \kappa u v^{4} - 2 \kappa u w^{4} + 3 c \kappa^{2} u^{4} + u^{5} + \kappa^{2} u^{5} - 2 \kappa u^{5}}{8 c^{2} - 8 \kappa c^{2}} & \frac{c^{2} v^{2} + c^{2} w^{2} - u^{2} v^{2} - u^{2} w^{2} + 3 c^{2} u^{2} - \kappa c^{2} v^{2} - \kappa c^{2} w^{2} - \kappa^{2} u^{2} v^{2} - \kappa^{2} u^{2} w^{2} - 5 \kappa c^{2} u^{2} - 2 c u v^{2} - 2 c u w^{2} + 2 \kappa u^{2} v^{2} + 2 \kappa u^{2} w^{2} + 2 c^{2} \kappa^{2} u^{2} - c u \kappa^{2} v^{2} - c u \kappa^{2} w^{2} + 3 c \kappa u v^{2} + 3 c \kappa u w^{2} - 4 c u^{3} + 2 u c^{3} - 3 c \kappa^{2} u^{3} + 7 c \kappa u^{3} - u^{4} - 2 c^{4} - \kappa^{2} u^{4} + 2 \kappa u^{4}}{4 c^{2} - 4 \kappa c^{2}} & \frac{- 6 c v u^{2} - 2 c v w^{2} - 2 u v w^{2} + 8 u v c^{2} - 4 \kappa u v c^{2} + 2 c \kappa v w^{2} + 2 \kappa u v w^{2} + 6 c \kappa v u^{2} - 2 c v^{3} - 2 u v^{3} - 2 v u^{3} + 4 v c^{3} + 2 c \kappa v^{3} + 2 \kappa u v^{3} + 2 \kappa v u^{3}}{8 c^{2}} & \frac{- 6 c w u^{2} - 2 c w v^{2} - 2 u w v^{2} + 8 u w c^{2} - 4 \kappa u w c^{2} + 2 c \kappa w v^{2} + 2 \kappa u w v^{2} + 6 c \kappa w u^{2} - 2 c w^{3} - 2 u w^{3} - 2 w u^{3} + 4 w c^{3} + 2 c \kappa w^{3} + 2 \kappa u w^{3} + 2 \kappa w u^{3}}{8 c^{2}} & \frac{2 c v^{2} + 2 c w^{2} + 2 u v^{2} + 2 u w^{2} + 6 c u^{2} - 6 c \kappa u^{2} - 2 c \kappa v^{2} - 2 c \kappa w^{2} - 2 \kappa u v^{2} - 2 \kappa u w^{2} + 4 \kappa u c^{2} - 4 c^{3} + 2 u^{3} - 2 \kappa u^{3}}{8 c^{2}}\end{smallmatrix}\right)

For u_1>0:

{\bf D}_x^- = \left(\begin{smallmatrix}u - c & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{smallmatrix}\right) {\bf A}_x^- = \left(\begin{smallmatrix}\frac{- 4 u c^{2} - 2 u v^{2} - 2 u w^{2} + 2 c v^{2} + 2 c w^{2} + 6 c u^{2} - 2 c \kappa u^{2} - 2 c \kappa v^{2} - 2 c \kappa w^{2} + 2 \kappa u v^{2} + 2 \kappa u w^{2} - 2 u^{3} + 2 \kappa u^{3}}{8 c^{2}} & \frac{- 2 c u + c \kappa u + c^{2} + u^{2} - \kappa u^{2}}{2 c^{2}} & \frac{u v - c v + c \kappa v - \kappa u v}{2 c^{2}} & \frac{u w - c w + c \kappa w - \kappa u w}{2 c^{2}} & \frac{c - u + \kappa u - c \kappa}{2 c^{2}}\\\frac{- 10 c^{2} u^{2} - 2 c^{2} v^{2} - 2 c^{2} w^{2} - 2 u^{2} v^{2} - 2 u^{2} w^{2} + 2 \kappa c^{2} u^{2} + 2 \kappa c^{2} v^{2} + 2 \kappa c^{2} w^{2} + 2 \kappa u^{2} v^{2} + 2 \kappa u^{2} w^{2} + 4 c u v^{2} + 4 c u w^{2} - 4 c \kappa u v^{2} - 4 c \kappa u w^{2} + 4 u c^{3} + 8 c u^{3} - 4 c \kappa u^{3} - 2 u^{4} + 2 \kappa u^{4}}{8 c^{2}} & \frac{- 3 c u^{2} + 3 u c^{2} - \kappa u c^{2} + 2 c \kappa u^{2} + u^{3} - c^{3} - \kappa u^{3}}{2 c^{2}} & \frac{- 2 c u v + 2 c \kappa u v + v c^{2} + v u^{2} - \kappa v c^{2} - \kappa v u^{2}}{2 c^{2}} & \frac{- 2 c u w + 2 c \kappa u w + w c^{2} + w u^{2} - \kappa w c^{2} - \kappa w u^{2}}{2 c^{2}} & \frac{2 c u - 2 c \kappa u - c^{2} - u^{2} + \kappa c^{2} + \kappa u^{2}}{2 c^{2}}\\\frac{- 4 u v c^{2} - 2 u v w^{2} + 2 c v w^{2} + 6 c v u^{2} - 2 c \kappa v u^{2} - 2 c \kappa v w^{2} + 2 \kappa u v w^{2} - 2 u v^{3} - 2 v u^{3} + 2 c v^{3} - 2 c \kappa v^{3} + 2 \kappa u v^{3} + 2 \kappa v u^{3}}{8 c^{2}} & \frac{- 2 c u v + c \kappa u v + v c^{2} + v u^{2} - \kappa v u^{2}}{2 c^{2}} & \frac{u v^{2} - c v^{2} + c \kappa v^{2} - \kappa u v^{2}}{2 c^{2}} & \frac{u v w - c v w + c \kappa v w - \kappa u v w}{2 c^{2}} & \frac{c v - u v + \kappa u v - c \kappa v}{2 c^{2}}\\\frac{- 4 u w c^{2} - 2 u w v^{2} + 2 c w v^{2} + 6 c w u^{2} - 2 c \kappa w u^{2} - 2 c \kappa w v^{2} + 2 \kappa u w v^{2} - 2 u w^{3} - 2 w u^{3} + 2 c w^{3} - 2 c \kappa w^{3} + 2 \kappa u w^{3} + 2 \kappa w u^{3}}{8 c^{2}} & \frac{- 2 c u w + c \kappa u w + w c^{2} + w u^{2} - \kappa w u^{2}}{2 c^{2}} & \frac{u v w - c v w + c \kappa v w - \kappa u v w}{2 c^{2}} & \frac{u w^{2} - c w^{2} + c \kappa w^{2} - \kappa u w^{2}}{2 c^{2}} & \frac{c w - u w + \kappa u w - c \kappa w}{2 c^{2}}\\\frac{- 2 c^{3} u^{2} - 2 c^{3} v^{2} - 2 c^{3} w^{2} - 2 u^{3} v^{2} - 2 u^{3} w^{2} - 2 \kappa c^{3} u^{2} - 2 u c^{2} v^{2} - 2 u c^{2} w^{2} - 2 u v^{2} w^{2} - 2 \kappa^{2} u^{3} v^{2} - 2 \kappa^{2} u^{3} w^{2} + 2 c v^{2} w^{2} + 2 \kappa c^{3} v^{2} + 2 \kappa c^{3} w^{2} + 4 \kappa u^{3} v^{2} + 4 \kappa u^{3} w^{2} + 6 c u^{2} v^{2} + 6 c u^{2} w^{2} - 10 c \kappa u^{2} v^{2} - 10 c \kappa u^{2} w^{2} - 4 c \kappa v^{2} w^{2} - 2 u c^{2} \kappa^{2} v^{2} - 2 u c^{2} \kappa^{2} w^{2} - 2 u \kappa^{2} v^{2} w^{2} + 2 c \kappa^{2} v^{2} w^{2} + 4 c \kappa^{2} u^{2} v^{2} + 4 c \kappa^{2} u^{2} w^{2} + 4 \kappa u c^{2} v^{2} + 4 \kappa u c^{2} w^{2} + 4 \kappa u v^{2} w^{2} - 6 c^{2} u^{3} - 2 c^{2} \kappa^{2} u^{3} + 8 \kappa c^{2} u^{3} + c v^{4} + c w^{4} - u v^{4} - u w^{4} + 4 u c^{4} + 5 c u^{4} + c \kappa^{2} v^{4} + c \kappa^{2} w^{4} - u \kappa^{2} v^{4} - u \kappa^{2} w^{4} - 8 c \kappa u^{4} - 2 c \kappa v^{4} - 2 c \kappa w^{4} + 2 \kappa u v^{4} + 2 \kappa u w^{4} + 3 c \kappa^{2} u^{4} - u^{5} - \kappa^{2} u^{5} + 2 \kappa u^{5}}{8 c^{2} - 8 \kappa c^{2}} & \frac{c^{2} v^{2} + c^{2} w^{2} + u^{2} v^{2} + u^{2} w^{2} + 3 c^{2} u^{2} + \kappa^{2} u^{2} v^{2} + \kappa^{2} u^{2} w^{2} - \kappa c^{2} v^{2} - \kappa c^{2} w^{2} - 5 \kappa c^{2} u^{2} - 2 c u v^{2} - 2 c u w^{2} - 2 \kappa u^{2} v^{2} - 2 \kappa u^{2} w^{2} + 2 c^{2} \kappa^{2} u^{2} - c u \kappa^{2} v^{2} - c u \kappa^{2} w^{2} + 3 c \kappa u v^{2} + 3 c \kappa u w^{2} - 4 c u^{3} + 2 u c^{3} - 3 c \kappa^{2} u^{3} + 7 c \kappa u^{3} + u^{4} - 2 c^{4} + \kappa^{2} u^{4} - 2 \kappa u^{4}}{4 c^{2} - 4 \kappa c^{2}} & \frac{- 6 c v u^{2} - 2 c v w^{2} + 2 u v w^{2} - 4 \kappa u v c^{2} - 2 \kappa u v w^{2} + 2 c \kappa v w^{2} + 6 c \kappa v u^{2} - 2 c v^{3} + 2 u v^{3} + 2 v u^{3} + 4 v c^{3} - 2 \kappa u v^{3} - 2 \kappa v u^{3} + 2 c \kappa v^{3}}{8 c^{2}} & \frac{- 6 c w u^{2} - 2 c w v^{2} + 2 u w v^{2} - 4 \kappa u w c^{2} - 2 \kappa u w v^{2} + 2 c \kappa w v^{2} + 6 c \kappa w u^{2} - 2 c w^{3} + 2 u w^{3} + 2 w u^{3} + 4 w c^{3} - 2 \kappa u w^{3} - 2 \kappa w u^{3} + 2 c \kappa w^{3}}{8 c^{2}} & \frac{- 2 u v^{2} - 2 u w^{2} + 2 c v^{2} + 2 c w^{2} + 6 c u^{2} - 6 c \kappa u^{2} - 2 c \kappa v^{2} - 2 c \kappa w^{2} + 2 \kappa u v^{2} + 2 \kappa u w^{2} + 4 \kappa u c^{2} - 4 c^{3} - 2 u^{3} + 2 \kappa u^{3}}{8 c^{2}}\end{smallmatrix}\right)

Boundary Conditions for the Sea Breeze Model

In the boundary (line) integral we prescribe w_3^{n+1} using a Dirichlet condition and calculate it at each iteration using:

w_3^{n+1} = E = \rho T c_v + \half \rho u^2 = w_0 T c_v + {w_1^2+w_2^2\over 2w_0}

where T(t) is a known function of time (it changes with the day and night) and also prescribe w_1^{n+1}=0 on the left and right end of the domain and w_2^{n+1}=0 at the top and bottom.

All the surface integrals turn out to be zero. On the top and bottom edges we have {\bf n} = (n_x, n_y) = (0, \pm 1) respectively and we prescribe w_2=0, so we get (remember we do not sum over i):

\int_{\partial\Omega} \left({\bf A}_x({\bf w}^n)\right)_{ij}w_j \varphi^i\, n_x + \left({\bf A}_y({\bf w}^n)\right)_{ij}w_j \varphi^i\, n_y \ \d x = = \int_{\partial\Omega} \left({\bf f}_x({\bf w}^n)\right)_i \varphi^i\, n_x + \left({\bf f}_y({\bf w}^n)\right)_i \varphi^i\, n_y \ \d x = = \pm\int_{\partial\Omega} \left({\bf f}_y({\bf w}^n)\right)_i \varphi^i \ \d x

where:

{\bf f}_y = \left( \begin{array}{c} w_2\\ \frac{w_2w_1}{w_0}\\ \frac{w_2^2}{w_0} + p\\ \frac{w_2}{w_0}(w_3+p) \end{array} \right) = \left( \begin{array}{c} 0\\ 0\\ p\\ 0 \end{array} \right)

So all the components i\neq 3 of the surface integral are zero, and for i=3 the test function \varphi^3 is not there, because we prescribe the Dirichlet BC w^3=0, so the surface integral vanishes for all i.

Similarly on the left and right edges we have {\bf n} = (n_x, n_y) = (\pm1, 0) respectively and we prescribe w_1=0, so we get (remember we do not sum over i):

\int_{\partial\Omega} \left({\bf A}_x({\bf w}^n)\right)_{ij}w_j \varphi^i\, n_x + \left({\bf A}_y({\bf w}^n)\right)_{ij}w_j \varphi^i\, n_y \ \d x = = \int_{\partial\Omega} \left({\bf f}_x({\bf w}^n)\right)_i \varphi^i\, n_x + \left({\bf f}_y({\bf w}^n)\right)_i \varphi^i\, n_y \ \d x = = \pm\int_{\partial\Omega} \left({\bf f}_x({\bf w}^n)\right)_i \varphi^i \ \d x

where:

{\bf f}_x = \left( \begin{array}{c} w_1\\ \frac{w_1^2}{w_0} + p\\ \frac{w_1w_2}{w_0}\\ \frac{w_1}{w_0}(w_3+p) \end{array} \right) = \left( \begin{array}{c} 0\\ p\\ 0\\ 0 \end{array} \right)

So all the components i\neq 1 of the surface integral are zero, and for i=1 the test function \varphi^1 is not there, because we prescribe the Dirichlet BC w^1=0, so the surface integral vanishes for all i.

7.3.9. Newton Method

The residual is:

F_{i,m}(Y^{n+1}) = \int_\Omega {w_{i,m}(y_m^{n+1}) - w_{i,m}(y^n)\over\tau} \varphi_{i, m} -f_{x,m}(w(y^n)){\partial \varphi_{i, m}\over\partial x} -f_{y,m}(w(y^n)){\partial \varphi_{i, m}\over\partial y} +\delta_{3, m} g \varphi_{i, m} \,\d x\, \d y + - \int_{\partial\Omega} f_{x,m}(w(y^n))\varphi_{i, m}\nu_x +f_{y,m}(w(y^n))\varphi_{i, m}\nu_y \, \d S = 0

where m = 0, 1, 2, 4 numbers the equations, i = 1, 2, ..., M numbers the finite element basis functions, N = 4M, Y = (y_0^1, y_1^1, y_2^1, y_3^1, y_0^2, y_1^2, ...). The Jacobian is:

J(Y^n) = {\partial F_{i, m}\over\partial y_{r, s}} (Y^n) = \int_\Omega {\varphi_{r, s}\over\tau} \varphi_{i, m} -A_{x, m, s}(w(y^n)) \varphi_{r, s} {\partial\varphi_{i, m}\over\partial x} -A_{y, m, s}(w(y^n)) \varphi_{r, s} {\partial\varphi_{i, m}\over\partial y} \,\d x\,\d y + \int_{\partial\Omega} A_{x, m, s}(w(y^n)) \varphi_{r, s} \varphi_{i, m}\nu_x + A_{y, m, s}(w(y^n)) \varphi_{r, s} \varphi_{i, m}\nu_y \,\d S

And the Newton method then is:

J(Y^n) \delta Y^{n+1} = -F(Y^n)

7.3.10. Older notes

Author: Pavel Solin

Governing Equations and Boundary Conditions

(7.3.10.1)\frac{\partial}{\partial t} \left( \begin{array}{c} \varrho\\ U\\ W\\ \theta \end{array} \right) + \frac{\partial}{\partial x} \left( \begin{array}{c} U\\ \frac{U^2}{\varrho} + R\theta\\ \frac{UW}{\varrho}\\ \frac{\theta U}{\varrho} \end{array} \right) + \frac{\partial}{\partial z} \left( \begin{array}{c} W\\ \frac{UW}{\varrho}\\ \frac{W^2}{\varrho} + R\theta\\ \frac{\theta W}{\varrho} \end{array} \right) + \left( \begin{array}{c} 0\\ 0\\ \varrho g\\ \frac{R\theta}{c_v}\mbox{div}{\bf v} \end{array} \right) = \left( \begin{array}{c} 0\\ 0\\ 0\\ 0 \end{array} \right),

where \varrho is the air density, {\bf v} = (u,w) is the velocity, U = \varrho u, W = \varrho w, T is the temperature, \theta = \varrho T, and g is the gravitational acceleration constant. We use the perfect gas state equation p = \varrho R T = R \theta for the pressure.

Boundary conditions are prescribed as follows:

  • edge a: \partial \varrho / \partial \nu = 0, \partial U / \partial \nu = 0, W = 0, \theta = \mbox{tanh}(x)*\mbox{sin}(\pi t /86400)

  • edges b, c: \partial \varrho / \partial \nu = 0, U = 0, \partial W / \partial \nu = 0, \partial \theta/ \partial \nu = 0

  • edge d: \partial \varrho / \partial \nu = 0, \partial U / \partial \nu = 0, W = 0, \partial \theta/ \partial \nu = 0

Initial conditions have the form

\begin{eqnarray*} p(z) &=& p_0 - 11476\frac{z}{1000} + 529.54 \left(\frac{z}{1000} \right)^2 - 9.38 \left(\frac{z}{1000} \right)^3,\\ T(z) &=& T_0 - 8.3194 \frac{z}{1000} + 0.2932 \left(\frac{z}{1000} \right)^2 - 0.0109 \left(\frac{z}{1000} \right)^3,\\ \varrho(z) &=& \frac{p(z)}{R T(z)},\\ \theta(z) &=& \varrho(z)T(z),\\ U(z) &=& 0, \\ W(z) &=& 0. \end{eqnarray*}

Discretization and the Newton’s Method

We will use the implicit Euler method in time, i.e.,

\frac{\partial \varrho}{\partial t} \approx \frac{\varrho^{n+1} - \varrho^n}{\tau}

etc. Let’s discuss one equation of (7.3.10.1) at a time:

Continuity equation: The weak formulation of

\frac{\varrho^{n+1} - \varrho^n}{\tau} + \frac{\partial U^{n+1}}{\partial x} + \frac{\partial W^{n+1}}{\partial z} = 0

reads

(7.3.10.2)F_i^{\varrho}(Y^{n+1}) = \int_{\Omega} \frac{\varrho^{n+1}}{\tau} \varphi^{\varrho}_i - \int_{\Omega} \frac{\varrho^{n}}{\tau} \varphi^{\varrho}_i + \int_{\Omega} \frac{\partial U^{n+1}}{\partial x} \varphi^{\varrho}_i + \int_{\Omega} \frac{\partial W^{n+1}}{\partial z} \varphi^{\varrho}_i = 0

The global coefficient vector Y^{n+1} consists of four parts Y^{\varrho}, Y^{U}, Y^{W} and Y^{\theta} corresponding to the fields \varrho, U, W and \theta, respectively. The same holds for the vector function F which consists of four parts F^{\varrho}, F^{U}, F^{W} and F^{\theta}. Thus the global Jacobi matrix will have a four-by-four block structure. We denote

(7.3.10.3)\varrho^{n+1} = \sum_{k=1}^{N^{\varrho}} y^{\varrho}_k \varphi^{\varrho}_k, \ \ \ U^{n+1} = \sum_{k=1}^{N^{U}} y^{U}_k \varphi^{U}_k, \ \ \ W^{n+1} = \sum_{k=1}^{N^{W}} y^{W}_k \varphi^{W}_k, \ \ \ \theta^{n+1} = \sum_{k=1}^{N^{\theta}} y^{\theta}_k \varphi^{\theta}_k.

It follows from (7.3.10.2) and (7.3.10.3) that

\frac{\partial F^{\varrho}_i}{\partial y^{\varrho}_j} = \int_{\Omega} \frac{\varphi^{\varrho}_j}{\tau} \varphi^{\varrho}_i, \ \ \ \frac{\partial F^{\varrho}_i}{\partial y^{U}_j} = \int_{\Omega} \frac{\partial \varphi^{U}_j}{\partial x} \varphi^{\varrho}_i, \ \ \ \frac{\partial F^{\varrho}_i}{\partial y^{W}_j} = \int_{\Omega} \frac{\partial \varphi^{W}_j}{\partial z} \varphi^{\varrho}_i, \ \ \ \frac{\partial F^{\varrho}_i}{\partial y^{W}_j} = 0.

First momentum equation: The second equation of (7.3.10.1) has the form

\frac{\partial U}{\partial t} + \frac{2U}{\varrho}\frac{\partial U}{\partial x} - \frac{U^2}{\varrho^2} \frac{\partial \varrho}{\partial x} + R\frac{\partial \theta}{\partial x} + \frac{W}{\varrho}\frac{\partial U}{\partial z} + \frac{U}{\varrho}\frac{\partial W}{\partial z} - \frac{UW}{\varrho^2}\frac{\partial \varrho}{\partial z} = 0.

After applying the implicit Euler method, we obtain

\frac{\partial U^{n+1}}{\tau} - \frac{\partial U^{n}}{\tau} + \frac{2U^{n+1}}{\varrho^{n+1}}\frac{\partial U^{n+1}}{\partial x} - \frac{(U^{n+1})^2}{(\varrho^{n+1})^2} \frac{\partial \varrho^{n+1}}{\partial x} + R\frac{\partial \theta^{n+1}}{\partial x}

+ \frac{W^{n+1}}{\varrho^{n+1}}\frac{\partial U^{n+1}}{\partial z} + \frac{U^{n+1}}{\varrho^{n+1}}\frac{\partial W^{n+1}}{\partial z} - \frac{U^{n+1}W^{n+1}}{(\varrho^{n+1})^2}\frac{\partial \varrho^{n+1}}{\partial z} = 0.

Thus we obtain

\frac{\partial F^{U}_i}{\partial y^{\varrho}_j} = - \int_{\Omega}\frac{2U}{\varrho^2}\frac{\partial U}{\partial x} \varphi^{\varrho}_j \varphi^{U}_i - \int_{\Omega} U^2 \left[(-2)\frac{1}{\varrho^3}\frac{\partial \varrho}{\partial x} \varphi^{\varrho}_j + \frac{1}{\varrho^2}\frac{\partial \varphi^{\varrho}_j}{\partial x}\right] \varphi^U_i

+ \int_{\Omega} \frac{W}{\varrho^2}\frac{\partial U}{\partial z}(-1)\varphi^{\varrho}_j \varphi^U_i + \int_{\Omega} \frac{U}{\varrho^2}\frac{\partial W}{\partial z}(-1)\varphi^{\varrho}_j \varphi^U_i - \int_{\Omega} UW \left[(-2)\frac{1}{\varrho^3}\frac{\partial \varrho}{\partial z} \varphi^{\varrho}_j + \frac{1}{\varrho^2}\frac{\partial \varphi^{\varrho}_j}{\partial z} \right] \varphi^{U}_i.

Analogously,

\frac{\partial F^{U}_i}{\partial y^{U}_j} = \int_{\Omega}\frac{\varphi^U_j}{\tau}\varphi^U_i + \int_{\Omega}\frac{2}{\varrho} \left[ \frac{\partial U}{\partial x}\varphi^U_j + U \frac{\partial \varphi^U_j}{\partial x} \right] \varphi^U_i - \int_{\Omega} \frac{2U}{\varrho^2}\frac{\partial \varrho}{\partial x} \varphi^U_j \varphi^U_i

+ \int_{\Omega} \frac{W}{\varrho}\frac{\partial \varphi^U_j}{\partial z} \varphi^U_i + \int_{\Omega} \frac{1}{\varrho}\frac{\partial W}{\partial z} \varphi^U_j \varphi^U_i - \int_{\Omega} \frac{W}{\varrho^2}\frac{\partial \varrho}{\partial z} \varphi^U_j \varphi^U_i,

\frac{\partial F^{U}_i}{\partial y^{W}_j} = \int_{\Omega} \frac{1}{\varrho}\frac{\partial U}{\partial z} \varphi^W_j \varphi^U_i + \int_{\Omega} \frac{U}{\varrho}\frac{\partial \varphi^W_j}{\partial z} \varphi^U_i - \int_{\Omega} \frac{U}{\varrho^2}\frac{\partial \varrho}{\partial z} \varphi^W_j \varphi^U_i,

\frac{\partial F^{U}_i}{\partial y^{\theta}_j} = \int_{\Omega} R \frac{\partial \varphi^{\theta}_j}{\partial x} \varphi^U_i.

Second momentum equation: The third equation of (7.3.10.1) reads

\frac{\partial W}{\partial t} + \frac{W}{\varrho}\frac{\partial U}{\partial x} + \frac{U}{\varrho}\frac{\partial W}{\partial x} - \frac{UW}{\varrho^2}\frac{\partial \varrho}{\partial x} + \frac{2W}{\varrho}\frac{\partial W}{\partial z} - \frac{W^2}{\varrho^2} \frac{\partial \varrho}{\partial x} + R\frac{\partial \theta}{\partial z} + \varrho g= 0.

After applying the implicit Euler method, we obtain

\frac{\partial W^{n+1}}{\tau} - \frac{\partial W^{n}}{\tau} + \frac{W^{n+1}}{\varrho^{n+1}}\frac{\partial U^{n+1}}{\partial x} + \frac{U^{n+1}}{\varrho^{n+1}}\frac{\partial W^{n+1}}{\partial x} - \frac{U^{n+1}W^{n+1}}{(\varrho^{n+1})^2}\frac{\partial \varrho^{n+1}}{\partial x}

+ \frac{2W^{n+1}}{\varrho^{n+1}}\frac{\partial W^{n+1}}{\partial z} - \frac{(W^{n+1})^2}{(\varrho^{n+1})^2} \frac{\partial \varrho^{n+1}}{\partial x} + R\frac{\partial \theta^{n+1}}{\partial z} + \varrho^{n+1} g= 0.

Thus we obtain

\frac{\partial F^{W}_i}{\partial y^{\varrho}_j} = + \int_{\Omega} \frac{W}{\varrho^2}\frac{\partial U}{\partial x}(-1)\varphi^{\varrho}_j \varphi^W_i + \int_{\Omega} \frac{U}{\varrho^2}\frac{\partial W}{\partial x}(-1)\varphi^{\varrho}_j \varphi^W_i - \int_{\Omega}\frac{2W}{\varrho^2}\frac{\partial W}{\partial x} \varphi^{\varrho}_j \varphi^{W}_i

- \int_{\Omega} UW \left[(-2)\frac{1}{\varrho^3}\frac{\partial \varrho}{\partial x} \varphi^{\varrho}_j + \frac{1}{\varrho^2}\frac{\partial \varphi^{\varrho}_j}{\partial x} \right] \varphi^{W}_i - \int_{\Omega} W^2 \left[(-2)\frac{1}{\varrho^3}\frac{\partial \varrho}{\partial z} \varphi^{\varrho}_j + \frac{1}{\varrho^2}\frac{\partial \varphi^{\varrho}_j}{\partial z}\right] \varphi^W_i + \int_{\Omega}g \varphi^{\varrho}_j \varphi^{W}_i.

Analogously,

\frac{\partial F^{W}_i}{\partial y^{U}_j} = \int_{\Omega} \frac{W}{\varrho}\frac{\partial \varphi^U_j}{\partial x} \varphi^W_i + \int_{\Omega} \frac{1}{\varrho}\frac{\partial W}{\partial x} \varphi^U_j \varphi^W_i - \int_{\Omega} \frac{W}{\varrho^2}\frac{\partial \varrho}{\partial x} \varphi^U_j \varphi^W_i,

\frac{\partial F^{W}_i}{\partial y^{W}_j} = \int_{\Omega}\frac{\varphi^W_j}{\tau}\varphi^W_i + \int_{\Omega} \frac{1}{\varrho}\frac{\partial U}{\partial x} \varphi^W_j \varphi^W_i + \int_{\Omega} \frac{U}{\varrho}\frac{\partial \varphi^W_j}{\partial x} \varphi^W_i - \int_{\Omega} \frac{U}{\varrho^2}\frac{\partial \varrho}{\partial x} \varphi^W_j \varphi^W_i

+ \int_{\Omega}\frac{2}{\varrho} \left[ \frac{\partial W}{\partial z}\varphi^W_j + W \frac{\partial \varphi^W_j}{\partial z} \right] \varphi^W_i - \int_{\Omega} \frac{2W}{\varrho^2}\frac{\partial \varrho}{\partial z} \varphi^W_j \varphi^W_i,

\frac{\partial F^{W}_i}{\partial y^{\theta}_j} = \int_{\Omega} R \frac{\partial \varphi^{\theta}_j}{\partial z} \varphi^W_i.

Internal energy equation: The last equation of (7.3.10.1) has the form

\frac{\partial \theta}{\partial t} + \mbox{div}(\theta {\bf v}) + \frac{R \theta}{c_v} \mbox{div}{\bf v} = 0

where \theta = \varrho T. This can be written equivalently as

\frac{\partial \theta}{\partial t} + \nabla \theta \cdot {\bf v} + \gamma \theta \mbox{div} {\bf v} = 0.

Written in terms of single derivatives, this is

\frac{\partial \theta}{\partial t} + \frac{\partial \theta}{\partial x} \frac{U}{\varrho} + \frac{\partial \theta}{\partial z} \frac{W}{\varrho} + \gamma \theta \frac{\partial}{\partial x}\left(\frac{U}{\varrho} \right) + \gamma \theta \frac{\partial}{\partial z}\left(\frac{W}{\varrho} \right) = 0,

i.e.,

\frac{\partial \theta}{\partial t} + \frac{\partial \theta}{\partial x} \frac{U}{\varrho} + \frac{\partial \theta}{\partial z} \frac{W}{\varrho} + \gamma \frac{\theta}{\varrho} \frac{\partial U}{\partial x} - \gamma \frac{\theta U}{\varrho^2}\frac{\partial \varrho}{\partial x} + \gamma \frac{\theta}{\varrho} \frac{\partial W}{\partial z} - \gamma \frac{\theta W}{\varrho^2}\frac{\partial \varrho}{\partial z} = 0.

Weak formulation:

F^{\theta}_i(Y) = \int_{\Omega} \frac{\theta^{n+1}}{\tau} \varphi^{\theta}_i - \int_{\Omega} \frac{\theta^{n}}{\tau} \varphi^{\theta}_i + \int_{\Omega} \frac{\partial \theta^{n+1}}{\partial x} \frac{U^{n+1}}{\varrho^{n+1}}\varphi^{\theta}_i + \int_{\Omega} \frac{\partial \theta^{n+1}}{\partial z} \frac{W^{n+1}}{\varrho^{n+1}} \varphi^{\theta}_i

+ \int_{\Omega} \gamma \frac{\theta^{n+1}}{\varrho^{n+1}} \frac{\partial U^{n+1}}{\partial x}\varphi^{\theta}_i - \int_{\Omega} \gamma \frac{\theta^{n+1} U^{n+1}}{(\varrho^{n+1})^2}\frac{\partial \varrho^{n+1}}{\partial x}\varphi^{\theta}_i + \int_{\Omega} \gamma \frac{\theta^{n+1}}{\varrho^{n+1}} \frac{\partial W^{n+1}}{\partial z}\varphi^{\theta}_i -\int_{\Omega} \gamma \frac{\theta^{n+1} W^{n+1}}{(\varrho^{n+1})^2}\frac{\partial \varrho^{n+1}}{\partial z} \varphi^{\theta}_i= 0.

For the derivatives of the weak form we obtain:

\frac{\partial F^{\theta}_i}{\partial y^{\varrho}_j} = - \int_{\Omega} \frac{\partial \theta}{\partial x} \frac{U}{\varrho^2}\varphi^{\varrho}_j\varphi^{\theta}_i - \int_{\Omega} \frac{\partial \theta}{\partial z} \frac{W}{\varrho^2}\varphi^{\varrho}_j\varphi^{\theta}_i - \int_{\Omega} \gamma \frac{\theta}{\varrho^2} \frac{\partial U}{\partial x}\varphi^{\varrho}_j\varphi^{\theta}_i - \int_{\Omega} \gamma \frac{\theta}{\varrho^2} \frac{\partial W}{\partial z}\varphi^{\varrho}_j\varphi^{\theta}_i

+ \int_{\Omega} 2\gamma \frac{\theta U}{\varrho^3}\frac{\partial \varrho}{\partial x}\varphi^{\varrho}_j\varphi^{\theta}_i - \int_{\Omega} \gamma \frac{\theta U}{\varrho^2}\frac{\varphi^{\varrho}_j}{\partial x}\varphi^{\theta}_i + \int_{\Omega} 2\gamma \frac{\theta W}{\varrho^3}\frac{\partial \varrho}{\partial z}\varphi^{\varrho}_j\varphi^{\theta}_i - \int_{\Omega} \gamma \frac{\theta W}{\varrho^2}\frac{\varphi^{\varrho}_j}{\partial z}\varphi^{\theta}_i.

\frac{\partial F^{\theta}_i}{\partial y^{U}_j} = \int_{\Omega} \frac{\partial \theta}{\partial x} \frac{1}{\varrho} \varphi^{U}_j\varphi^{\theta}_i + \int_{\Omega} \gamma \frac{\theta}{\varrho}\frac{\varphi^{U}_j}{\partial x}\varphi^{\theta}_i - \int_{\Omega} \gamma \frac{\theta}{\varrho^2}\frac{\partial \varrho}{\partial x}\varphi^{U}_j\varphi^{\theta}_i.

\frac{\partial F^{\theta}_i}{\partial y^{W}_j} = \int_{\Omega} \frac{\partial \theta}{\partial z} \frac{1}{\varrho} \varphi^{W}_j\varphi^{\theta}_i + \int_{\Omega} \gamma \frac{\theta}{\varrho}\frac{\varphi^{W}_j}{\partial z}\varphi^{\theta}_i - \int_{\Omega} \gamma \frac{\theta}{\varrho^2}\frac{\partial \varrho}{\partial z}\varphi^{W}_j\varphi^{\theta}_i.

\frac{\partial F^{\theta}_i}{\partial y^{\theta}_j} = \int_{\Omega} \frac{1}{\tau} \varphi^{\theta}_j\varphi^{\theta}_i + \int_{\Omega} \frac{U}{\varrho}\frac{\varphi^{\theta}_j}{\partial x}\varphi^{\theta}_i + \int_{\Omega} \frac{W}{\varrho}\frac{\varphi^{\theta}_j}{\partial z}\varphi^{\theta}_i

+ \int_{\Omega} \frac{\gamma}{\varrho} \frac{\partial U}{\partial x} \varphi^{\theta}_j\varphi^{\theta}_i + \int_{\Omega} \frac{\gamma}{\varrho} \frac{\partial W}{\partial z} \varphi^{\theta}_j\varphi^{\theta}_i - \int_{\Omega} \frac{\gamma U}{\varrho^2} \frac{\partial \varrho}{\partial x} \varphi^{\theta}_j\varphi^{\theta}_i - \int_{\Omega} \frac{\gamma W}{\varrho^2} \frac{\partial \varrho}{\partial z} \varphi^{\theta}_j\varphi^{\theta}_i.