Flux Calculations in Cylindrical Coordinates

[hwhitehead]

I'm looking at the calculations performed to include rotating frame effects in cylindrical polar coordinates, specifically the function void HydroSourceTerms::RotatingSystemSourceTerms.

When calculating the momentum updates in the X1 ($r$) and X2 ($\phi$) directions, the two appear to be given different treatments.

Notably, when calculating the momentum change in X2, the fluxes appear to be marginalized over the cell edges and center. Specifically, the term $\rho v_r$ is represented here as an average of mom1 = den * prim(IVX,k,j,i) (the cell center radial flux) and flux_c = 0.5*(flux[X1DIR](IDN,k,j,i)+flux[X1DIR](IDN,k,j,i+1)) (the cell face radial flux).

Conversely, when calculating the momentum change in X1, the face fluxes are not considered and $\rho v_\phi$ is written here as den * prim(IVY,k,j,i).

Would anyone be able to explain why face flux average is performed in the $\phi$ direction but not the $r$ direction? I assume it has to do with conservation in some form, but I cannot see why.

[c-white]

A partial answer: The form of the rotational terms is somewhat consistent with that of the nonrotational curvilinear source terms in that an effort is made to exactly conserve angular momentum $R \rho v^\phi$ rather than azimuthal momentum $\rho v^\phi$, but not to do anything special for radial momentum $\rho v^R$.

The radial momentum equation in differential form is

$$\partial_t (\rho v^R) + \frac{1}{R} \partial_R (R (\rho v^R v^R + p)) + \frac{1}{R} \partial_\phi (\rho v^R v^\phi) + \partial_z (\rho v^R v^z) = \frac{1}{R} (\rho v^\phi v^\phi + p) + \Omega^2 R \rho + 2 \Omega \rho v^\phi.$$

In integral form, this is exactly

$$\partial_t \langle \rho v^R \rangle_{R\phi z}$$

$$+ \frac{1}{V_{R\phi z}} [\langle \rho v^R v^R + p \rangle_{\phi z} A_{\phi z}]_{R^-}^{R^+}$$

$$+ \frac{1}{V_{R\phi z}} [\langle \rho v^R v^\phi \rangle_{Rz} A_{Rz}]_{\phi^-}^{\phi^+}$$

$$+ \frac{1}{V_{R\phi z}} [\langle \rho v^R v^z \rangle_{R\phi} A_{R\phi}]_{z^-}^{z^+}$$

$$= \langle \frac{1}{R} (\rho v^\phi v^\phi + p) \rangle_{R\phi z} + \Omega^2 \langle R \rho \rangle_{R\phi z} + 2 \Omega \langle \rho v^\phi \rangle_{R\phi z}.$$

The code is evaluating the first source term as

$$\langle \frac{1}{R} (\rho v^\phi v^\phi + p) \rangle_{R\phi z} = \frac{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} (\rho v^\phi v^\phi + p)\ dR\ d\phi\ dz}{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} R\ dR\ d\phi\ dz} \approx \frac{(R_+ - R_-)}{(R_+^2 - R_-^2) / 2} [\rho v^\phi v^\phi + p]_0,$$

centrifugal term as

$$\Omega^2 \langle R \rho \rangle_{R\phi z} = \Omega^2 \frac{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} R^2 \rho\ dR\ d\phi\ dz}{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} R\ dR\ d\phi\ dz} \approx \Omega^2 \frac{(R_+ - R_-)}{(R_+^2 - R_-^2) / 2} [R^2 \rho]_0,$$

and the Coriolis term as

$$2 \Omega \langle \rho v^\phi \rangle_{R\phi z} = 2 \Omega \frac{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} R \rho v^\phi\ dR\ d\phi\ dz}{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} R\ dR\ d\phi\ dz} \approx 2 \Omega \frac{(R_+ - R_-)}{(R_+^2 - R_-^2) / 2} [R \rho v^\phi]_0.$$

All three approximations are consistent. While the appearance of $\rho v^\phi$ in the Coriolis term suggests one could write this in terms of $\phi$-fluxes of mass, there is little to be gained by doing this.

The azimuthal momentum equation in differential form is

$$\partial_t (\rho v^\phi) + \frac{1}{R} \partial_R (R \rho v^R v^phi) + \frac{1}{R} \partial_\phi (\rho v^\phi v^\phi + p) + \partial_z (\rho v^\phi v^z) = -\frac{1}{R} \rho v^R v^\phi - 2 \Omega \rho v^R.$$

In integral form, this is exactly

$$\partial_t \langle \rho v^\phi \rangle_{R\phi z}$$

$$+ \frac{1}{V_{R\phi z}} [\langle \rho v^R v^\phi \rangle_{\phi z} A_{\phi z}]_{R^-}^{R^+}$$

$$+ \frac{1}{V_{R\phi z}} [\langle \rho v^\phi v^\phi + p\rangle_{Rz} A_{Rz}]_{\phi^-}^{\phi^+}$$

$$+ \frac{1}{V_{R\phi z}} [\langle \rho v^\phi v^z \rangle_{R\phi} A_{R\phi}]_{z^-}^{z^+}$$

$$= -\langle \frac{1}{R} (\rho v^R v^\phi) \rangle_{R\phi z} - 2 \Omega \langle \rho v^R \rangle_{R\phi z}.$$

The code is evaluating the first source term as

$$-\langle \frac{1}{R} (\rho v^R v^\phi) \rangle_{R\phi z} = -\frac{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} \rho v^R v^\phi\ dR\ d\phi\ dz}{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} R\ dR\ d\phi\ dz} \approx -\frac{(R_+ - R_-)}{(R_+^2 - R_-^2) / 2}$$

$$\times (\frac{R_-}{R_- + R_+} [\langle \rho v^R v^\phi \rangle_{\phi z}]_{R^-}$$

$$+ \frac{R_+}{R_- + R_+} [\langle \rho v^R v^\phi \rangle_{\phi z}]_{R^+}),$$

which indirectly keeps the evolution of $R \rho v^\phi$ exact at the expense of $\rho v^\phi$. The derivation for this can be found in Chapter 2 here: https://dataspace.princeton.edu/handle/88435/dsp017s75df84b

Analogously, the code is evaluating the Coriolis term as

$$-2 \Omega \langle \rho v^R \rangle_{R\phi z} = -2 \Omega \frac{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} R \rho v^R\ dR\ d\phi\ dz}{\int_{z^-}^{z^+} \int_{\phi^-}^{\phi^+} \int_{R^-}^{R^+} R\ dR\ d\phi\ dz} \approx -2 \Omega \frac{(R_+ - R_-)}{(R_+^2 - R_-^2) / 2} R_0$$

$$\times (\frac{1}{4} [\langle \rho v^R \rangle_{\phi z}]_{R^-}$$

$$+ \frac{1}{2} [\rho v^R]_0$$

$$+ \frac{1}{4} [\langle \rho v^R \rangle_{\phi z}]_{R^+}).$$

Why this exact form, which looks like a three-point numerical quadrature? Why not use only the two flux terms and omit the middle term? I'm not sure, and I don't know any source that has exactly derived these equations.

[hwhitehead]

Thanks for laying this all out so clearly, this is incredibly helpful. A preference for conserving $R \rho v^\phi$ makes sense, though I confess I am somewhat out of my depth with the discretised derivatives. I'll try to do some more reading around the numerical theory.