Hegedus Aerodynamics

Note: the equations below require MathML to be displayed.

Artificial Dissipation

Since central differencing is used for Aero Troll CFD, artificial dissipation is added to control the odd-even uncoupling of the grid and to control nonlinear effects such as shocks (1).  The nonlinear scalar model is adapted from reference (2) and from information provided in reference (3).

$$\frac{\left(1+\varphi \right)}{J}\frac{\Delta Q^{n+1}}{\Delta t}+{\left(\frac{\partial A}{\partial \xi }+\frac{\partial B}{\partial \eta }+\frac{\partial C}{\partial \zeta }-\frac{\partial D_{\xi }}{\partial Q}-\frac{\partial D_{\eta }}{\partial Q}-\frac{\partial D_{\zeta }}{\partial Q}\right)}^n\Delta Q^{n+1}=\frac{\varphi }{J}\frac{\Delta Q^n}{\Delta t}-{\left(\frac{\partial \widehat{F}}{\partial \xi }+\frac{\partial \widehat{G}}{\partial \eta }+\frac{\partial \widehat{H}}{\partial \zeta }-D_{\xi }-D_{\eta }-D_{\zeta }\right)}^n$$ $$D_{\xi }={\nabla }_{\xi }\left(\mathrm{d2nd}{\Delta }_{\xi }{Q}_D-\mathrm{d4th}{\Delta }_{\xi }{\nabla }_{\xi }{\Delta }_{\xi }{Q}_D\right)$$ $$D_{\eta }={\nabla }_{\eta }\left(\mathrm{d2nd}{\Delta }_{\eta }{Q}_D-\mathrm{d4th}{\Delta }_{\eta }{\nabla }_{\eta }{\Delta }_{\eta }{Q}_D\right)$$ $$D_{\zeta }={\nabla }_{\zeta }\left(\mathrm{d2nd}{\Delta }_{\zeta }{Q}_D-\mathrm{d4th}{\Delta }_{\zeta }{\nabla }_{\zeta }{\Delta }_{\zeta }{Q}_D\right)$$ $${\mathrm{d2nd}}_{i+\frac{1}{2}}={\lambda }_{i+\frac{1}{2}\text{,}j\text{,}k}{ϵ}_{i+\frac{1}{2}\text{,}j\text{,}k}^{\left(2\right)}$$ $${\mathrm{d4th}}_{i+\frac{1}{2}}={\lambda }_{i+\frac{1}{2}\text{,}j\text{,}k}{ϵ}_{i+\frac{1}{2}\text{,}j\text{,}k}^{\left(4\right)}$$ The numerical implementation of the spectral radius, $\lambda$, in Aero Troll CFD deviates from what seems to be the standard implementation and instead uses the default implementation as stated in the OVERFLOW manual (3), i.e. the spectral radius in the direction of differentiation. $${\lambda }_{i+\frac{1}{2}\text{,}j\text{,}k}=\frac{1}{2}\left({\left({\lambda }_{\xi }\right)}_{i+1\text{,}j\text{,}k}+{\left({\lambda }_{\xi }\right)}_{i\text{,}j\text{,}k}\right)$$ $$\left({\lambda }_{\xi }\right)=2^{\left(\text{*}\right)}{J}^{-1}\left(\left|\overrightarrow{\xi }·\overrightarrow{V}\right|+c\left|\overrightarrow{\xi }\right|\right)$$ Note(*) a "2" has been added in front of the spectral radius equation above.  The value exists because of historical reasons related to testing Aero Troll CFD.  The original idea was to keep the dissipation constants close to OVERFLOW to try to parallel some published OVERFLOW results.  However, it turns out that some cases go unstable using the published dissipation constants, therefore the dissipation constants were internally doubled by inserting the "2."  The decision was made not to remove the "2" and not redefine the dissipation constants for the public release of Aero Troll CFD so old test cases would not break.

The formulation for ${ϵ}^{\left(\text{2}\right)}$ for Aero Troll CFD varies from what is provided in literature.  Literature provides the following formulation: $${ϵ}_{i+\frac{1}{2}\text{,}j\text{,}k}^{\left(\text{2}\right)}={\kappa }^{\left(\text{2}\right)}\max \left\{{\nu }_{i+1}\text{,}\nu \right\}$$ Instead, Aero Troll CFD uses the following formulation for ${ϵ}^{\left(\text{2}\right)}$: $${ϵ}_{i+\frac{1}{2}\text{,}j\text{,}k}^{\left(\text{2}\right)}=\frac{{\kappa }^{\left(2\right)}}{2}\left({\nu }_{i+1}+{\nu }_i\right)$$ The reason the different formulation is used is that the formulation provided in literature has problems converging for some supersonic cases due to discrete switching between the ${\nu }_{i+1}$ and ${\nu }_i$ values.  On such case for a NACA 0012 airfoil at Mach 1.4 can be seen here. $${ϵ}_{i+\frac{1}{2}\text{,}j\text{,}k}^{\left(4\right)}=\text{max}\left\{0\text{,}\left({\kappa }^{\left(4\right)}-{ϵ}_{i+\frac{1}{2}\text{,}j\text{,}k}^{\left(2\right)}\right)\right\}$$ $${\kappa }^{\left(2\right)}=2.0\text{,}{\kappa }^{\left(4\right)}=0.04$$ The values for the dissipation constants ${\kappa }^{\left(2\right)}$ and ${\kappa }^{\left(4\right)}$ shown above come from the OVERFLOW manual (3). $${\nu }_i=\frac{\left|p_{i+1\text{,}j\text{,}k}-2p_{i\text{,}j\text{,}k}+p_{i-1\text{,}j\text{,}k}\right|}{p_{i+1\text{,}j\text{,}k}+2p_{i\text{,}j\text{,}k}+p_{i-1\text{,}j\text{,}k}}$$ $$Q_D=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho w\\ \rho h_0\end{array}\right]=Q+\left[\begin{array}{c}0.0\\ 0.0\\ 0.0\\ 0.0\\ p\end{array}\right]$$

---

Discretization

$$D_{\xi }=\left({\mathrm{d2nd}}_{i+\frac{1}{2}}\left(Q_{D\text{,}i+1}-Q_{D\text{,}i}\right)-{\mathrm{d4th}}_{i+\frac{1}{2}}\left(Q_{D\text{,}i+2}-3Q_{D\text{,}i+1}+3Q_{D\text{,}i}-Q_{D\text{,}i-1}\right)\right)-\left({\mathrm{d2nd}}_{i-\frac{1}{2}}\left(Q_{D\text{,}i}-Q_{D\text{,}i-1}\right)-{\mathrm{d4th}}_{i-\frac{1}{2}}\left(Q_{D\text{,}i+1}-3Q_{D\text{,}i}+3Q_{D\text{,}i-1}-Q_{D\mathrm{,}i-2}\right)\right)$$

---

Axisymmetric Artificial Dissipation

The terms $\mathrm{d2nd}$, $\mathrm{d4th}$, $\rho$, and $\rho h_0$ are constant in the η direction, therefore:

$$D_{\eta }={\mathrm{d2nd}}_{\eta }\frac{{\partial }^2\left(\rho \overrightarrow{V}\right)}{\partial {\eta }^2}-{\mathrm{d4th}}_{\eta }\frac{{\partial }^4\left(\rho \overrightarrow{V}\right)}{\partial {\eta }^4}$$ Where:

$$\frac{\partial \left(\rho \overrightarrow{V}\right)}{\partial \eta }=\Delta {\overrightarrow{\rho V}}_{\eta }={\overrightarrow{\rho V}}_{\mathrm{norm}}\times \Delta \overrightarrow{\Theta }$$ and:

$$\frac{{\partial }^2\left(\rho \overrightarrow{V}\right)}{\partial {\eta }^2}={\Delta }^2\overrightarrow{\rho V_{\eta }}=-\left(\Delta \overrightarrow{\Theta }·\Delta \overrightarrow{\Theta }\right){\overrightarrow{\rho V}}_{\mathrm{norm}}$$ and:

$$\frac{{\partial }^3\left(\rho \overrightarrow{V}\right)}{\partial {\eta }^3}={\Delta }^3\overrightarrow{\rho V_{\eta }}=-\left(\Delta \overrightarrow{\Theta }·\Delta \overrightarrow{\Theta }\right)\left({\overrightarrow{\rho V}}_{\mathrm{norm}}\times \Delta \overrightarrow{\Theta }\right)$$ and:

$$\frac{{\partial }^4\left(\rho \overrightarrow{V}\right)}{\partial {\eta }^4}={\Delta }^4\overrightarrow{\rho V_{\eta }}={\left(\Delta \overrightarrow{\Theta }·\Delta \overrightarrow{\Theta }\right)}^2{\overrightarrow{\rho V}}_{\mathrm{norm}}$$

---

Artificial Dissipation Jacobian

$$\frac{\partial Q_D}{\partial Q}=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ \frac{1}{2}\left(\gamma -1\right)\overrightarrow{V}·\overrightarrow{V}& -\left(\gamma -1\right)u& -\left(\gamma -1\right)v& -\left(\gamma -1\right)w& \gamma \end{array}\right]$$ The max eigenvalue is $\gamma$.

---

Diagonal Artificial Dissipation Jacobian

$$\frac{\partial Q_D}{\partial Q_k}=T_k\left(I+\left(\gamma -1\right)T_k^{-1}\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ \frac{1}{2}\overrightarrow{V}·\overrightarrow{V}& -u& -v& -w& 1\end{array}\right]T_k\right)T_k^{-1}$$ $$\frac{\partial Q_D}{\partial Q_k}=T_k\left(I+\frac{\rho c}{\sqrt{2}}{T}_k^{-1}\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 1\end{array}\right]\right)T_k^{-1}$$ $$\frac{\partial Q_D}{\partial Q_k}=T_k\left(I+\left[\begin{array}{ccccc}0& 0& 0& -{\overrightarrow{k}}_x\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}& -{\overrightarrow{k}}_x\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}\\ 0& 0& 0& -{\overrightarrow{k}}_y\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}& -{\overrightarrow{k}}_y\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}\\ 0& 0& 0& -{\overrightarrow{k}}_z\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}& -{\overrightarrow{k}}_z\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}\\ 0& 0& 0& \frac{1}{2}\left(\gamma -1\right)& \frac{1}{2}\left(\gamma -1\right)\\ 0& 0& 0& \frac{1}{2}\left(\gamma -1\right)& \frac{1}{2}\left(\gamma -1\right)\end{array}\right]\right)T_k^{-1}$$ $$\frac{\partial Q_D}{\partial Q_k}=T_k\left[\begin{array}{ccccc}1& 0& 0& -{\overrightarrow{k}}_x\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}& -{\overrightarrow{k}}_x\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}\\ 0& 1& 0& -{\overrightarrow{k}}_y\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}& -{\overrightarrow{k}}_y\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}\\ 0& 0& 1& -{\overrightarrow{k}}_z\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}& -{\overrightarrow{k}}_z\left(\gamma -1\right)\frac{\rho }{c\sqrt{2}}\\ 0& 0& 0& \frac{1}{2}\left(\gamma +1\right)& \frac{1}{2}\left(\gamma -1\right)\\ 0& 0& 0& \frac{1}{2}\left(\gamma -1\right)& \frac{1}{2}\left(\gamma +1\right)\end{array}\right]T_k^{-1}$$ It can be seen that the last two equations are uncoupled from the first three, therefore the eigenvalues for the last two equations are:

$$\left[\frac{1}{2}\left(\gamma +1\right)-\lambda \right]\left[\frac{1}{2}\left(\gamma +1\right)-\lambda \right]-\left[\frac{1}{2}\left(\gamma -1\right)\right]\left[\frac{1}{2}\left(\gamma -1\right)\right]=0$$ $${\lambda }_4=1\text{,}{\lambda }_5=\gamma$$ The spectral radius, i.e. the maximum eigenvalue, for the last two equations is γ.

The eigenvalues for the first three equations are:

$${\lambda }_1={\lambda }_2={\lambda }_3=1$$ Therefore the spectral radius for the first three equations is 1.

$$\frac{\partial Q_D}{\partial Q_k}\approx T_k\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& \gamma & 0\\ 0& 0& 0& 0& \gamma \end{array}\right]T_k^{-1}$$ The above diagonalized artificial dissipation Jacobian was compared to results for one were the diagonal had all γs.  In general, the two forms of the diagonalized Jacobian behave similar with the one with all γs converging slightly faster on some results.  Therefore the diagonalized artificial dissipation Jacobian in Aero Troll contains all γs on the diagonal.

References

1) Pulliam, T. H., "Artificial Dissipation Models for the Euler Equations," AIAA Journal, Vol. 24, No. 12, Dec. 1986.
2) Slooff, J. W, and Schmidt, W., "Computational Aerodynamics Based on the Euler Equations," AGARD-AG-325, AGARD, Sept. 1984.
3) Nichols, R. H., and Buning P. G., "User's Manual for OVERFLOW 2.1, Version 2.1t Aug 2008," Aug 2008.