Difference between revisions of "User:Tohline/PGE/Hybrid Scheme"

From VistrailsWiki
Jump to navigation Jump to search
(→‎Preface: Begin to abbreviate preface and create new page with full version of preface)
 
(54 intermediate revisions by the same user not shown)
Line 6: Line 6:
''March 1, 2014'' by Joel E. Tohline
''March 1, 2014'' by Joel E. Tohline


<[[User:Tohline/PGE/Hybrid_Scheme_Preface|Click here for full preface]]> &hellip; Throughout my research career, I have sought new or modified techniques and algorithms that would allow my group to perform more accurate numerical simulations of astrophysical fluid flows. At the beginning &#8212; following the advice of my dissertation co-advisors, Peter Bodenheimer and David Black &#8212; I adopted a cylindrical, rather than Cartesian, computational grid.  When a cylindrical coordinate system is used, one of the components of the equation of motion can be written in a form that fairly naturally conserves angular momentum, and this was quite a desirable feature, given that our investigations were focusing on analyzing the stability of rotating configurations. 
<[[User:Tohline/PGE/Hybrid_Scheme_Preface|full preface]]> &hellip; This seemed too good to be true.  The discovered code modification would allow us to conserve angular momentum very accurately and, at the same time, allow us to use a rotating grid and thereby minimize numerical diffusion &hellip; <[[User:Tohline/PGE/Hybrid_Scheme_Preface|read more]]>


During the mid-90s, when Kimberly Barker New was conducting her dissertation research, we started using a ''rotating'' cylindrical coordinate mesh. On a grid that was spinning at a suitably chosen frequency, advection of fluid ''through'' the grid could be minimized and this, in turn, reduced the undesirable effects of numerical diffusion.  We adopted a fairly standard algorithmic approach very similar to the one used by Norman &amp; Wilson (1980) and acknowledged that we were making a tradeoff:  While the shift to a rotating cylindrical coordinate mesh reduced the effects of numerical diffusion, the shift introduced a rather ugly Coriolis "source" term into two components of the equation of motionThis made it more difficult to ensure conservation of angular momentum.
In his [http://etd.lsu.edu/docs/available/etd-07052010-224427/ dissertation research], Jay Call (see [http://adsabs.harvard.edu/abs/2010CQGra..27q5002C Call, Tohline, &amp; Lehner 2010]) derived a complete description of this hybrid advection scheme in a fully relativistic and generalized coordinate framework.  He showed that one can write the system of fluid equations in a manner that facilitates advection of inertial-frame quantities across a rotating gridIn addition &#8212; and more importantly &#8212; he showed how to write the system of fluid equations to allow advection of inertial-frame angular momentum (generally associated with a cylindrical coordinate mesh) across a rotating ''Cartesian'' gridSince that time, primarily in the context of [http://etd.lsu.edu/docs/available/etd-04022014-210318/ Zach Byerly's doctoral dissertation research], my group has demonstrated that this hybrid advection scheme works extremely well (see [http://adsabs.harvard.edu/abs/2014ApJS..212...23B Byerly, Adelstein-Lelbach, Tohline, &amp; Marcello (2014)])In the discussion that follows here, we derive the Newtonian version of Jay's hybrid scheme.   
 
During my continuing efforts to develop an improved computational fluid dynamics algorithm, one day I noticed that, while imposing some fairly reasonable constraints, the Coriolis term could be removed from the "source" term and folded into the divergence term on the left-hand side of the angular momentum conservation equation.  This manipulation of terms seemed to be saying that the undesirable Coriolis term would disappear while employing a rotating coordinate mesh if the variable that was advected through the grid was the ''inertial-frame'' angular momentum density, rather than the ''rotating-frame'' angular momentum densityThis seemed too good to be true.  The discovered code modification would allow us to conserve angular momentum very accurately and, at the same time, allow us to use a rotating grid and thereby minimize numerical diffusion.  My notes on this topic have been preserved, as they were included in my earliest version of this web-based H_Book; the relevant page can be accessed [http://www.phys.lsu.edu/astro/H_Book.current/Context/PGE/angmom.conserve.html here, which is an html file whose linux time stamp is August 27, 2000].  The "symbol" fonts utilized throughout this old html page seem now to be readable only through Microsoft's Internet Explorer web browserHence, for posterity sake, I have retyped this "year 2000" set of notes into an [[User:Tohline/PGE/AStarScheme#Hybrid_Advection_Scheme_.28Background.29|accompanying page of this wiki]].




Line 405: Line 403:




As foreshadowed above &#8212; see the [[User:Tohline/PGE/Hybrid_Scheme#Recognizing_Statements_of_Conservatio|subsection titled, ''Recognizing Statements of Conservation'']] &#8212; the angular momentum of a Lagrangian fluid element will be conserved if the "source" term, <math>~S = 0</math>.  This situation will arise if, at the fluid element's location, the azimuthal pressure variation, <math>~\partial P/\partial\phi</math>, and the azimuthal variation in the gravitational potential, <math>~\partial \Phi/\partial\phi</math>, are both zero, or if the two balance one another (''i.e.,''<math>~\partial P/\partial\phi=-\rho\partial\Phi/\partial\phi</math>).
As foreshadowed above &#8212; see the [[#Recognizing_Statements_of_Conservation|subsection titled, ''Recognizing Statements of Conservation'']] &#8212; the angular momentum of a Lagrangian fluid element will be conserved if the "source" term, <math>~S = 0</math>.  This situation will arise if, at the fluid element's location, the azimuthal pressure variation, <math>~\partial P/\partial\phi</math>, and the azimuthal variation in the gravitational potential, <math>~\partial \Phi/\partial\phi</math>, are both zero, or if the two balance one another (''i.e.,''<math>~\partial P/\partial\phi=-\rho\partial\Phi/\partial\phi</math>).


Based on the above discussion, we can equally well view the flow from a frame of reference that is rotating with a constant angular velocity, <math>~\Omega_0</math>, and write,  
Based on the above discussion, we can equally well view the flow from a frame of reference that is rotating with a constant angular velocity, <math>~\Omega_0</math>, and write,  
Line 735: Line 733:
In the latter case, <math>~v'_\phi</math> becomes <math>~u_\phi</math>, as it should, when the choice is made to measure the angular momentum density in the "grid" frame, that is, when the choice is made to set <math>~\omega_0 = \Omega_0</math>.
In the latter case, <math>~v'_\phi</math> becomes <math>~u_\phi</math>, as it should, when the choice is made to measure the angular momentum density in the "grid" frame, that is, when the choice is made to set <math>~\omega_0 = \Omega_0</math>.


=TRASH (do not read below this line)=
=Zach's Dissertation Paper=
==Traditional Eulerian Representation (Review)==
== Section 2.3==
Here we review the traditional Eulerian representation of the Euler Equation, as has been [[User:Tohline/PGE/Euler#Eulerian_Representation|discussed in detail earlier]].
Building on our introductory discussion of the [[User:Tohline/PGE/RotatingFrame#Euler_Equation_.28rotating_frame.29|Euler equation]] (see also Appendix 1.D, &sect;3 of [[User:Tohline/Appendix/References|BT87]]), we begin with the,


===in terms of velocity:===
<div align="center">
<font color="#770000">'''Lagrangian Representation'''</font><br />
of the Euler Equation <br />
<font color="#770000">'''as viewed from a Rotating Reference Frame'''</font>


The so-called "Eulerian form" of the Euler equation can be straightforwardly derived from the [[User:Tohline/PGE/Euler#Lagrangian_Representation|standard Lagrangian representation]] to obtain,
<math>\frac{d\mathbf{u'}}{dt} = - \frac{1}{\rho} \nabla P - \nabla \Phi - 2{\vec{\Omega}}_f \times {\mathbf{u'}} - {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x})</math> ,
</div>
where we choose to define the frame rotation by the vector,
<div align="center">
<math>~\vec{\Omega}_f \equiv \hat\mathbf{k} \Omega_0 \, ,</math>
</div>
and,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
~\mathbf{u'}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
~\mathbf{u} - \mathbf{u}_\mathrm{frame} = \mathbf{u} - {\vec{\Omega}}_f \times \vec{x}\, ,
</math>
  </td>
</tr>
</table>
</div>
is the velocity field as viewed from the frame of reference that is rotating at constant angular frequency, <math>~\Omega_0</math>.  (Note that we can retrieve the inertial-frame Euler equation and inertial-frame variables by setting <math>~\Omega_0 = 0</math> at any point in the subsequent derivations.)  Because the velocity field introduced by frame rotation is divergence free, that is, because,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
~\nabla\cdot\mathbf{u}_\mathrm{frame}  =
\nabla\cdot [{\vec{\Omega}}_f \times \vec{x}] = \nabla\cdot [ {\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{e}}_R R + \mathbf{\hat{k}}z) ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
~\nabla\cdot(\hat\mathbf{e}_\varphi R\Omega_0) = \frac{1}{R} \frac{\partial}{\partial\varphi} \biggl(R\Omega_0 \biggr) = 0 \, ,
</math>
  </td>
</tr>
</table>
</div>
the relevant (rotating-frame) continuity equation is identical in form to its inertial-frame counterpart, specifically,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\frac{d\rho}{dt} + \rho\nabla\cdot\mathbf{u'}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
~0 \, .
</math>
  </td>
</tr>
</table>
</div>
We can transform the above rotating-frame Euler equation to a momentum conservation equation by multiplying the equation through by <math>~\rho</math> and using the continuity equation to combine and simplify the left-hand-side, obtaining,


<div align="center">
<div align="center">
<span id="ConservingMomentum:Eulerian"><font color="#770000">'''Eulerian Representation'''</font></span><br />
<math>\frac{d(\rho\mathbf{u'})}{dt} + (\rho\mathbf{u' })\nabla\cdot \mathbf{u'} = - \nabla P - \rho \nabla \Phi - 2\rho {\vec{\Omega}}_f \times {\mathbf{u'}} - \rho {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x})</math> ,
of the Euler Equation,
 
{{User:Tohline/Math/EQ_Euler02}}
</div>
</div>


===in terms of momentum density:===
In order to convert this general-purpose vector equation into the specific set of scalar component equations that embody the desired elements of our hybrid scheme, we need to:
 
* Step 1:  Choose the unit-vector basis set associated with the momentum components that we want to track &#8212; <math>(\mathbf{\hat{i}}, \mathbf{\hat{j}}, \mathbf\hat{k})</math> if tracking linear momentum, or, <math>(\mathbf\hat{e}_R, \mathbf\hat{e}_\varphi, \mathbf\hat{k})</math> if tracking radial and angular momentum &#8212; and break the vector equation into these specified components;
Also, we can multiply this expression through by {{User:Tohline/Math/VAR_Density01}} and combine it with the continuity equation to derive what is commonly referred to as the,
* Step 2:  In all relevant equations, replace the scalar components of the "rotating-frame" momentum density with the scalar components of the "intertial-frame" momentum density, drawing each component relation from the vector transformation, <math>\rho\mathbf{u'} = \rho\mathbf{u} - \rho {\vec{\Omega}}_f \times \vec{x}</math>;
* Step 3:  Write all of the spatial operators in terms of spatial derivatives that are associated with the unit-vector basis set of the desired computational mesh.


===Focus on Tracking Linear Momentum ===
====Step 1====
If the focus is on tracking linear momentum components, then we need to break the vector momentum equation into its <math>(\mathbf\hat{i}, \mathbf\hat{j}, \mathbf\hat{k})</math> components.  This is done by, in turn, "dotting" each unit vector into the vector equation.  It is straightforward once we appreciate that the orientation of these Cartesian unit vectors does not vary in space and that, within the context of the rotating frame on which they are defined, these unit vectors do not vary in time.  Hence, the first term in the vector equation &#8212; the ''material'' time derivative &#8212; can be written as,
<div align="center">
<div align="center">
<span id="ConservingMomentum:Conservative"><font color="#770000">'''Conservative Form'''</font></span><br />
<math>
of the Euler Equation,
\frac{d(\rho\mathbf{u'})}{dt} =  
 
\frac{d}{dt} [ \mathbf{\hat{i}} (\rho u'_x) + \mathbf{\hat{j}} (\rho u'_y) + \mathbf{\hat{k}} (\rho u'_z) ] =
{{User:Tohline/Math/EQ_Euler03}}
\mathbf{\hat{i}} \frac{d(\rho u'_x)}{dt}  + \mathbf{\hat{j}} \frac{d(\rho u'_y)}{dt}  + \mathbf{\hat{k}} \frac{d(\rho u'_z)}{dt}  \, ,
</math>
</div>
</div>
and the process of "dotting" each unit vector into the equation leads to the following set of scalar momentum-component equations:
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="center">
<math>\mathbf{\hat{i}}:</math>
  </td>
  <td align="right">
<math>
\frac{d(\rho u'_x)}{dt} + (\rho u'_x)\nabla\cdot \mathbf{u'}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + 2\rho \mathbf{\hat{i}}\cdot[\mathbf{\hat{i}}\Omega_0 u'_y - \mathbf{\hat{j}}\Omega_0 u'_x] + \rho \mathbf{\hat{i}}\cdot[ \mathbf{\hat{i}}\Omega_0^2 x + \mathbf{\hat{j}}\Omega_0^2 y] \,
</math>
  </td>
</tr>


The second term on the left-hand-side of this last expression represents the divergence of the "dyadic product" of the vector momentum density ({{User:Tohline/Math/VAR_Density01}}{{User:Tohline/Math/VAR_VelocityVector01}}) and the velocity vector {{User:Tohline/Math/VAR_VelocityVector01}} and is sometimes written as, <math>\nabla\cdot [(\rho \vec{v}) \otimes \vec{v}]</math>.
<tr>
  <td align="center">
&nbsp;
  </td>
  <td align="right">
<math>
\Rightarrow ~~~~~ \frac{\partial (\rho u'_x)}{\partial t} + \nabla\cdot [(\rho u'_x)\mathbf{u'}]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + 2\rho \Omega_0 u'_y + \rho \Omega_0^2 x \, ;
</math>
  </td>
</tr>


=Component Forms=
<tr>
  <td align="center">
<math>\mathbf{\hat{j}}:</math>
  </td>
  <td align="right">
<math>
\frac{d(\rho u'_y)}{dt} + (\rho u'_y)\nabla\cdot \mathbf{u'}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{j}}\cdot\nabla P - \rho \mathbf{\hat{j}}\cdot\nabla \Phi + 2\rho \mathbf{\hat{j}}\cdot[\mathbf{\hat{i}}\Omega_0 u'_y - \mathbf{\hat{j}}\Omega_0 u'_x] + \rho \mathbf{\hat{j}}\cdot[ \mathbf{\hat{i}}\Omega_0^2 x + \mathbf{\hat{j}}\Omega_0^2 y] \,
</math>
  </td>
</tr>


Let's split the vector Euler equation into its three scalar components; various examples are identified in Table 1.
<tr>
  <td align="center">
&nbsp;
  </td>
  <td align="right">
<math>
\Rightarrow ~~~~~ \frac{\partial (\rho u'_y)}{\partial t} + \nabla\cdot [(\rho u'_y) \mathbf{u'} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{j}}\cdot\nabla P - \rho \mathbf{\hat{j}}\cdot\nabla \Phi - 2\rho \Omega_0 u'_x + \rho \Omega_0^2 y \, ;
</math>
  </td>
</tr>


<div align="center">
<table border="1" cellpadding="5">
<tr>
<tr>
   <th align="center" rowspan="2">
   <td align="center">
Example #
<math>\mathbf{\hat{k}}:</math>
   </th>
  </td>
   <th align="center" colspan="2">
  <td align="right">
Grid
<math>
   </th>
\frac{d(\rho u'_z)}{dt} + (\rho u'_z)\nabla\cdot \mathbf{u'}
   <th align="center" colspan="2">
</math>
Momentum Vector
   </td>
   </th>
   <td align="center">
<math>~=~</math>
   </td>
   <td align="left">
<math>
- \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi + 2\rho \mathbf{\hat{k}}\cdot[\mathbf{\hat{i}}\Omega_0 u'_y - \mathbf{\hat{j}}\Omega_0 u'_x] + \rho \mathbf{\hat{k}}\cdot[ \mathbf{\hat{i}}\Omega_0^2 x + \mathbf{\hat{j}}\Omega_0^2 y] \,
</math>
   </td>
</tr>
</tr>


<tr>
<tr>
   <th align="center">
   <td align="center">
Basis
&nbsp;
   </th>
   </td>
   <th align="center">
   <td align="right">
Rotating?
<math>
   </th>
\Rightarrow ~~~~~ \frac{\partial (\rho u'_z)}{\partial t} + \nabla\cdot[(\rho u'_z) \mathbf{u'} ]
   <th align="center">
</math>
Basis
   </td>
   </th>
   <td align="center">
   <th align="center">
<math>~=~</math>
Frame
   </td>
   </th>
   <td align="left">
<math>
- \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi \, .
</math>
   </td>
</tr>
</tr>


</table>
</div>
In deriving these expressions, we also (a) have recognized from the start that, when expressed in Cartesian coordinates,
<div align="center">
<table border="0" cellpadding="3">
<tr>
<tr>
  <td align="right">
<math>
~{\vec{\Omega}}_f \times \vec{x}
</math>
  </td>
   <td align="center">
   <td align="center">
1
<math>~=~</math>
  </td>
  <td align="left">
<math>
{\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{i}}x + \mathbf{\hat{j}}y + \mathbf{\hat{k}}z)
=  - \mathbf{\hat{i}}\Omega_0 y + \mathbf{\hat{j}}\Omega_0 x \, ,
</math>
  </td>
</tr>
 
<tr>
  <td align="right">
<math>
{\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x})
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Cartesian
<math>~=~</math>
  </td>
  <td align="left">
<math>
\hat{\mathbf{k}} \Omega_0 \times ( - \mathbf{\hat{i}}\Omega_0 y + \mathbf{\hat{j}}\Omega_0 x )
=  - \mathbf{\hat{i}}\Omega_0^2 x - \mathbf{\hat{j}}\Omega_0^2 y \, ,
</math>
  </td>
</tr>
 
<tr>
  <td align="right">
<math>
{\vec{\Omega}}_f \times {\mathbf{u'}}
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
No
<math>~=~</math>
  </td>
  <td align="left">
<math>
{\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{i}}u'_x + \mathbf{\hat{j}}u'_y + \mathbf{\hat{k}}u'_z)
=  - \mathbf{\hat{i}}\Omega_0 u'_y + \mathbf{\hat{j}}\Omega_0 u'_x\, ,
</math>
  </td>
</tr>
 
</table>
</div>
and (b) have used the familiar operator mapping,
<div align="center">
<math>\frac{d}{dt} \rightarrow \frac{\partial}{\partial t} + \mathbf{u'}\cdot \nabla  \, ,</math>
</div>
to shift from the total (Lagrangian) time derivative to the partial (Eulerian) time derivative, which is usually the more desirable representation for computational simulations.
 
====Step 2====
Next, throughout this set of scalar equations, we replace each component of <math>~\rho\mathbf{u'}</math> with the corresponding component of <math>~(\rho\mathbf{u} - \rho {\vec{\Omega}}_f \times \vec{x})</math>, that is, we perform the following mappings:
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
~\rho u'_x
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Cartesian
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
~\rho (u_x +\Omega_0 y) \, ,
</math>
  </td>
</tr>
 
<tr>
  <td align="right">
<math>
~\rho u'_y
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Inertial
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
~\rho (u_y - \Omega_0 x) \, ,
</math>
   </td>
   </td>
</tr>
</tr>


<tr>
<tr>
  <td align="right">
<math>
~\rho u'_z
</math>
  </td>
   <td align="center">
   <td align="center">
2
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
~\rho u_z \, .
</math>
  </td>
</tr>
 
</table>
</div>
As a result, the first of the three "hybrid" momentum-component equations becomes,
<div align="center">
<table border="0" cellpadding="3">
 
<tr>
  <td align="right">
<math>
\frac{\partial [ \rho (u_x +\Omega_0 y) ] }{\partial t} + \nabla\cdot \{ [\rho (u_x +\Omega_0 y) ]\mathbf{u'}\}
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Cylindrical
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + 2\rho \Omega_0  (u_y - \Omega_0 x) + \rho \Omega_0^2 x
</math>
  </td>
</tr>
</table>
</div>
 
<div align="center">
<table border="0" cellpadding="3">
 
<tr>
  <td align="right">
<math>
\Rightarrow ~~~~~ \frac{\partial (\rho u_x)  }{\partial t} + \nabla\cdot [(\rho u_x) \mathbf{u'} ]
+ \Omega_0 y\biggl[ \frac{\partial \rho }{\partial t} + \nabla\cdot ( \rho  \mathbf{u'} ) \biggr]
+ ( \rho  \mathbf{u'} )\cdot\nabla (\Omega_0 y)
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Yes <math>~(\Omega_0)</math>
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + 2\rho \Omega_0 u_y - \rho \Omega_0^2 x \, .
</math>
  </td>
</tr>
</table>
</div>
Referencing the continuity equation, the middle bracketed term on the left-hand side can be set to zero; and the last term on the left-hand side,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
( \rho  \mathbf{u'} )\cdot\nabla (\Omega_0 y)
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Cylindrical
<math>~=~</math>
  </td>
  <td align="left">
<math>
\rho \Omega_0 u'_y = \rho\Omega_0(u_y - \Omega_0 x) \, ,
</math>
  </td>
</tr>
</table>
</div>
can be combined with terms on the right-hand side &#8212; cutting the Coriolis term in half and canceling the centrifugal acceleration term &#8212; to give,
 
<div align="center">
<table border="0" cellpadding="3">
 
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u_x)  }{\partial t} + \nabla\cdot [(\rho u_x) \mathbf{u'} ]
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Rotating <math>~(\Omega_0)</math>
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + \rho \Omega_0 u_y  \, .
</math>
   </td>
   </td>
</tr>
</tr>
</table>
</div>
Following a similar sequence of steps, the other two "hybrid" momentum conservation relations become,
<div align="center">
<table border="0" cellpadding="3">


<tr>
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u_y)  }{\partial t} + \nabla\cdot [(\rho u_y) \mathbf{u'} ]
</math>
  </td>
   <td align="center">
   <td align="center">
3
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{j}}\cdot\nabla P - \rho \mathbf{\hat{j}}\cdot\nabla \Phi - \rho \Omega_0 u_x  \, ,
</math>
  </td>
</tr>
 
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u_z)  }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ]
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Cylindrical
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi  \, .
</math>
  </td>
</tr>
</table>
</div>
 
====Step 3====
'''<font color="red">Cartesian Grid:</font>''' Assuming the numerical simulation will be conducted on a Cartesian coordinate mesh, the divergence (advection) term on the left-hand-side should be evaluated by breaking the transport velocity into its three Cartesian components,
<div align="center">
<math>~\mathbf{u'} = (u'_x, u'_y, u'_z) \, ,</math>
</div>
and, on the right-hand-side, the projection of the spatial operators should be written in the familiar form,
<div align="center">
<math>\mathbf{\hat{i}}\cdot\nabla \rightarrow \frac{\partial}{\partial x} \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{j}}\cdot\nabla \rightarrow \frac{\partial}{\partial y} \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{k}}\cdot\nabla \rightarrow \frac{\partial}{\partial z} \, .</math>
</div>
In summary, then, the relevant set of momentum conservation equations is,
 
<div align="center">
<table border="1" cellpadding="5" width="80%">
<tr>
  <td align="center" bgcolor="lightgreen">
'''Cartesian Components of the Inertial-Frame Momentum'''<br>
<font size="-1">advected across a</font><br>
'''Rotating, Cartesian Coordinate Mesh'''
  </td>
</tr>
<tr><td align="center">
<table border="0" cellpadding="3">
 
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u_x)  }{\partial t} + \nabla\cdot [(\rho u_x) \mathbf{u'} ]
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Yes <math>~(\Omega_0)</math>
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \frac{\partial}{\partial x}  P - \rho \frac{\partial}{\partial x} \Phi + \rho \Omega_0 u_y
</math>
  </td>
</tr>
 
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u_y)  }{\partial t} + \nabla\cdot [(\rho u_y) \mathbf{u'} ]
</math>  
   </td>
   </td>
   <td align="center">
   <td align="center">
Cylindrical
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \frac{\partial}{\partial y}  P - \rho \frac{\partial}{\partial y} \Phi - \rho \Omega_0 u_x
</math>
  </td>
</tr>
 
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u_z)  }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ]
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Rotating <math>~(\omega_0)</math>
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z}\Phi 
</math>
  </td>
</tr>
 
</table>
</td></tr>
<tr>
  <td align="left">
This is the set of equations that has served as the foundation of the ''Cartesian'' simulations reported in Byerly, Adelstein-Lelbach, Tohline, &amp; Marcello (2014).
   </td>
   </td>
</tr>
</tr>
Line 852: Line 1,274:
</div>
</div>


In the following expressions, we will use <math>~\vec{v}</math> to denote the fluid velocity when it is associated with the rate of fluid transport across the coordinate grid, and we will use <math>~\vec{u}</math> to denote the fluid velocity when it is associated with the momentum density that is being advected.  In all cases, it should be understood that <math>~\vec{v} = \vec{u}</math>, as both vectors refer to the same fluid velocity.  In addition, we will use a "prime" notation to indicate when a velocity is being viewed from a rotating frame of reference; specifically, we will consider rotation about the <math>~z</math>-axis of the coordinate system, that is,
 
'''<font color="red">Cylindrical Grid:</font>''' If, instead, the numerical simulation is to be conducted on a cylindrical coordinate mesh, the spatial operators on both sides of the component momentum equations should be broken down into their cylindrical-coordinate components.  In concert with this, the divergence (advection) term on the left-hand-side should be evaluated by breaking the transport velocity into its three cylindrical components,
<div align="center">
<math>~\mathbf{u'} = (u'_R, u'_\varphi, u'_z) \, .</math>
</div>
Furthermore, recognizing that, when written in cylindrical coordinates, the gradient operator is,
<div align="center">
<math>
\nabla = \mathbf{\hat{e}}_R \frac{\partial}{\partial R} +
\mathbf{\hat{e}}_\varphi \frac{1}{R} \frac{\partial}{\partial \varphi} +\mathbf{\hat{e}}_z \frac{\partial}{\partial z} \, ,
</math>
</div>
and that the unit vectors in cylindrical coordinates can be related to their Cartesian counterparts via the mappings,
<div align="center">
<math>\mathbf{\hat{e}}_R = \mathbf{\hat{i}} \cos\varphi + \mathbf{\hat{j}} \sin\varphi \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{e}}_\varphi = \mathbf{\hat{j}} \cos\varphi - \mathbf{\hat{i}} \sin\varphi \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{e}}_z = \mathbf{\hat{k}}  \, ,</math><br>
</div>
the relevant projections of the gradient operator on the right-hand-sides of the governing equations should take the form,
<div align="center">
<div align="center">
<table border="0" cellpadding="3">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\mathbf{\hat{i}}\cdot\nabla
</math>
  </td>
  <td align="center">
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
\biggl[ \cos\varphi \frac{\partial}{\partial R} -
\frac{\sin\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] \, ,
</math>
  </td>
</tr>
<tr>
  <td align="right">
<math>
\mathbf{\hat{j}}\cdot\nabla
</math>
  </td>
  <td align="center">
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
\biggl[ \sin\varphi \frac{\partial}{\partial R} +
\frac{\cos\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] \, ,
</math>
  </td>
</tr>
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~v'_\phi</math>
<math>
\mathbf{\hat{k}}\cdot\nabla
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~=~</math>
<math>~\rightarrow~</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~v_\phi - R\Omega_0 \, ,</math>
<math>
\frac{\partial}{\partial z} \, .
</math>
   </td>
   </td>
</tr>
</tr>
</table>
</table>
</div>
</div>
and,
 
In summary, then, the relevant set of momentum conservation equations is,
 
<div align="center">
<div align="center">
<table border="1" cellpadding="5">
<tr>
  <td align="center" bgcolor="lightgreen">
'''Cartesian Components of the Inertial-Frame Momentum'''<br>
<font size="-1">advected across a</font><br>
'''Rotating, Cylindrical Coordinate Mesh'''
  </td>
</tr>
<tr><td align="center">
<table border="0" cellpadding="3">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u_x)  }{\partial t} + \nabla\cdot [(\rho u_x) \mathbf{u'} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \biggl[ \cos\varphi \frac{\partial}{\partial R} -
\frac{\sin\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] P -
\rho \biggl[ \cos\varphi \frac{\partial}{\partial R} -
\frac{\sin\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] \Phi + \rho \Omega_0 u_y
</math>
  </td>
</tr>
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~u'_\phi</math>
<math>
\frac{\partial (\rho u_y)  }{\partial t} + \nabla\cdot [(\rho u_y) \mathbf{u'} ]
</math>  
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 879: Line 1,391:
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~u_\phi - R\omega_0 \, ,</math>
<math>
- \biggl[ \sin\varphi \frac{\partial}{\partial R} +
\frac{\cos\varphi}{R} \frac{\partial}{\partial \varphi}\biggr]  P - \rho \biggl[ \sin\varphi \frac{\partial}{\partial R} +
\frac{\cos\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] \Phi - \rho \Omega_0 u_x 
</math>
   </td>
   </td>
</tr>
</tr>
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u_z)  }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z}\Phi 
</math>
  </td>
</tr>
</table>
</td></tr>
</table>
</table>
</div>
</div>
but we will not insist that the two rotation frequencies, <math>~\Omega_0</math> and <math>~\omega_0</math>, have the same value.  Hence, in general, <math>~(\vec{u})' \ne (\vec{v})'</math>.  It is worth emphasizing that, because we will only be considering frame rotation about the <math>z</math>-axis, the cylindrical <math>R</math> and <math>z</math> components of the velocity are interchangeable, that is: <math>~u'_R = v'_R = u_R = v_R</math>; and <math>~u'_z = v'_z = u_z = v_z</math>.


===Example #1===
===Focus on Tracking Angular Momentum ===
This is certainly the most familiar component set.
====Step 1====
If the focus is on tracking angular momentum, then we need to break the vector momentum equation into its <math>(\mathbf\hat{e}_R, \mathbf\hat{e}_\varphi, \mathbf\hat{k})</math> components.  As before, this is done by "dotting" each unit vector into the vector equation. This is less straightforward than in the Cartesian case because the orientation of both the <math>\mathbf\hat{e}_R</math> and <math>\mathbf\hat{e}_\varphi</math> unit vectors vary in space.  As a result, the first term in the vector equation &#8212; the ''material'' time derivative &#8212; generates a couple of extra terms, viz.,
<div align="center">
<div align="center">
<table border="0" cellpadding="3">
<table border="0" cellpadding="3">
Line 893: Line 1,428:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>\boldsymbol{\hat{e}}_x: ~~~\frac{\partial (\rho v_x)}{\partial t} + \nabla\cdot[(\rho v_x) \vec{v}~]</math>
<math>
\frac{d(\rho\mathbf{u'})}{dt}  
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 900: Line 1,437:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial x} - \rho \frac{\partial \Phi}{\partial x} \, ,
\frac{d}{dt} [ \mathbf{\hat{e}}_R (\rho u'_R) + \mathbf{\hat{e}}_\varphi (\rho u'_\varphi) + \mathbf{\hat{k}} (\rho u'_z) ]
</math>
</math>
   </td>
   </td>
Line 907: Line 1,444:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>\boldsymbol{\hat{e}}_y: ~~~\frac{\partial (\rho v_y)}{\partial t} + \nabla\cdot[(\rho v_y) \vec{v}~]</math>
&nbsp;
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 914: Line 1,451:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial y} - \rho \frac{\partial \Phi}{\partial y} \, ,
\mathbf{\hat{e}}_R \frac{d(\rho u'_R)}{dt} + \mathbf{\hat{e}}_\varphi \frac{d(\rho u'_\varphi)}{dt}  + \mathbf{\hat{k}} \frac{d(\rho u'_z)}{dt} 
+ (\rho u'_R) \frac{d}{dt}\mathbf{\hat{e}}_R  + (\rho u'_\varphi) \frac{d}{dt}\mathbf{\hat{e}}_\varphi  \, ,  
</math>
</math>
   </td>
   </td>
Line 921: Line 1,459:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>\boldsymbol{\hat{e}}_z: ~~~\frac{\partial (\rho v_z)}{\partial t} + \nabla\cdot[(\rho v_z) \vec{v}~]</math>
&nbsp;
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 928: Line 1,466:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial z} - \rho \frac{\partial \Phi}{\partial z} \, ,
\mathbf{\hat{e}}_R \frac{d(\rho u'_R)}{dt}  + \mathbf{\hat{e}}_\varphi \frac{d(\rho u'_\varphi)}{dt}  + \mathbf{\hat{k}} \frac{d(\rho u'_z)}{dt} 
+ \mathbf{\hat{e}}_\varphi(\rho u'_R) \frac{u'_\varphi}{R}  - \mathbf{\hat{e}}_R(\rho u'_\varphi) \frac{u'_\varphi}{R}  \, .
</math>
</math>
   </td>
   </td>
Line 934: Line 1,473:
</table>
</table>
</div>
</div>
where, for any one of the three scalar PDEs, advection of the relevant component of the momentum density, <math>~\psi_i</math>, is handled via the operation,
We also recognize that, when expressed in cylindrical coordinates,  


<div align="center">
<div align="center">
Line 941: Line 1,480:
   <td align="right">
   <td align="right">
<math>
<math>
\nabla\cdot[\psi_{i} \vec{v} ]
~{\vec{\Omega}}_f \times \vec{x}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
{\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{e}}_R R + \mathbf{\hat{k}}z)
=  \mathbf{\hat{e}}_\varphi \Omega_0 R \, ,
</math>
  </td>
</tr>
 
<tr>
  <td align="right">
<math>
{\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x})
</math>
</math>
   </td>
   </td>
Line 949: Line 1,505:
   <td align="left">
   <td align="left">
<math>
<math>
\frac{\partial (\psi_i v_x)}{\partial x} + \frac{\partial (\psi_i v_y)}{\partial y} + \frac{\partial (\psi_i v_z)}{\partial z} \, .
\hat{\mathbf{k}} \Omega_0 \times ( \mathbf{\hat{e}}_\varphi \Omega_0 R )
=  - \mathbf{\hat{e}}_R \Omega_0^2 R \, ,
</math>
</math>
   </td>
   </td>
</tr>
</tr>
<tr>
  <td align="right">
<math>
{\vec{\Omega}}_f \times {\mathbf{u'}}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
{\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{e}}_R u'_R + \mathbf{\hat{e}}_\varphi u'_\varphi + \mathbf{\hat{k}}u'_z)
=  \mathbf{\hat{e}}_\varphi \Omega_0 u'_R - \mathbf{\hat{e}}_R \Omega_0 u'_\varphi  \, .
</math>
  </td>
</tr>
</table>
</table>
</div>
</div>


===Example #2===
Hence, the process of "dotting" each unit vector into the equation leads to the following set of scalar momentum-component equations:
This component set has been spelled out in, for example, equations (5) - (7) of [http://adsabs.harvard.edu/abs/1978ApJ...224..497N Norman &amp; Wilson (1978)] and equations (11), (12), &amp; (3) of  [http://adsabs.harvard.edu/abs/1997ApJ...490..311N New &amp; Tohline (1997)].
<div align="center">
<div align="center">
<table border="0" cellpadding="3">
<table border="0" cellpadding="3">
<tr>
  <td align="center">
<math>\mathbf{\hat{e}}_R:</math>
  </td>
  <td align="right">
<math>
\frac{d(\rho u'_R)}{dt} + (\rho u'_R)\nabla\cdot \mathbf{u'} - (\rho u'_\varphi) \frac{u'_\varphi}{R} 
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{e}}_R \cdot\nabla P - \rho \mathbf{\hat{e}}_R \cdot\nabla \Phi
- 2\rho \mathbf{\hat{e}}_R \cdot [ \mathbf{\hat{e}}_\varphi \Omega_0 u'_R - \mathbf{\hat{e}}_R \Omega_0 u'_\varphi  ]
+ \rho \mathbf{\hat{e}}_R\cdot (\mathbf{\hat{e}}_R \Omega_0^2 R) \,
</math>
  </td>
</tr>


<tr>
<tr>
  <td align="center">
&nbsp;
  </td>
   <td align="right">
   <td align="right">
<math>\boldsymbol{\hat{e}}_R: ~~~~~~~\frac{\partial (\rho v_R)}{\partial t} + \nabla\cdot[(\rho v_R) \vec{v}~]</math>
<math>
\Rightarrow ~~~~~ \frac{\partial (\rho u'_R)}{\partial t} + \nabla\cdot [(\rho u'_R)\mathbf{u'}]
</math>  
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 970: Line 1,569:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial R} - \rho \frac{\partial \Phi}{\partial R} + \frac{(\rho R v_\phi)^2}{\rho R^3} + \rho\Omega_0^2 R + \frac{2\Omega_0 (\rho R v_\phi)}{R} \, ,
- \mathbf{\hat{e}}_R \cdot\nabla P - \rho \mathbf{\hat{e}}_R \cdot\nabla \Phi  
+ \frac{\rho}{R} \biggl[ (u'_\varphi)^2 + 2R\Omega_0 u'_\varphi + (\Omega_0 R)^2 \biggr] \,  
</math>
</math>
   </td>
   </td>
Line 976: Line 1,576:


<tr>
<tr>
  <td align="center">
&nbsp;
  </td>
   <td align="right">
   <td align="right">
&nbsp;
&nbsp;
Line 984: Line 1,587:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial R} - \rho \frac{\partial \Phi}{\partial R} + \frac{\rho}{R} (v_\phi + R\Omega_0)^2 \, ,
- \mathbf{\hat{e}}_R \cdot\nabla P - \rho \mathbf{\hat{e}}_R \cdot\nabla \Phi  
+ \frac{\rho}{R} (u'_\varphi + R\Omega_0)^2 \, ;
</math>
</math>
   </td>
   </td>
Line 990: Line 1,594:


<tr>
<tr>
  <td align="center">
<math>\mathbf{\hat{e}}_\varphi:</math>
  </td>
   <td align="right">
   <td align="right">
<math>\boldsymbol{\hat{e}}_\phi: ~~~\frac{\partial (\rho R v_\phi)}{\partial t} + \nabla\cdot[(\rho R v_\phi) \vec{v}~]</math>
<math>
\frac{d(\rho u'_\varphi)}{dt} + (\rho u'_\varphi)\nabla\cdot \mathbf{u'} + (\rho u'_R) \frac{u'_\varphi}{R}
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 998: Line 1,607:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial \phi} - \rho \frac{\partial \Phi}{\partial \phi} - 2\rho (\Omega_0 R )v_R \, ,
- \mathbf{\hat{e}}_\varphi \cdot\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot\nabla \Phi  
- 2\rho \mathbf{\hat{e}}_\varphi \cdot [ \mathbf{\hat{e}}_\varphi \Omega_0 u'_R - \mathbf{\hat{e}}_R \Omega_0 u'_\varphi  ]
+ \rho \mathbf{\hat{e}}_\varphi\cdot (\mathbf{\hat{e}}_R \Omega_0^2 R) \,
</math>
</math>
   </td>
   </td>
Line 1,004: Line 1,615:


<tr>
<tr>
  <td align="center">
(mult. thru by R)
  </td>
   <td align="right">
   <td align="right">
<math>\boldsymbol{\hat{e}}_z: ~~~~~~~~\frac{\partial (\rho v_z)}{\partial t} + \nabla\cdot[(\rho v_z) \vec{v}~]</math>
<math>
\Rightarrow ~~~~~ \frac{d(\rho R u'_\varphi)}{dt} + (\rho R u'_\varphi)\nabla\cdot \mathbf{u'}
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 1,012: Line 1,628:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial z} - \rho \frac{\partial \Phi}{\partial z} \, ,
- \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi  
- 2\rho R \Omega_0 u'_R    \,  
</math>
</math>
   </td>
   </td>
</tr>
</tr>
<tr>
  <td align="center">
&nbsp;
  </td>
  <td align="right">
<math>
\Rightarrow ~~~~~ \frac{\partial (\rho R u'_\varphi)}{\partial t} + \nabla\cdot [(\rho R u'_\varphi) \mathbf{u'} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi
- 2\rho R \Omega_0 u'_R    \, ;
</math>
  </td>
</tr>
<tr>
  <td align="center">
<math>\mathbf{\hat{k}}:</math>
  </td>
  <td align="right">
<math>
\frac{d(\rho u'_z)}{dt} + (\rho u'_z)\nabla\cdot \mathbf{u'}
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{k}} \cdot\nabla P - \rho \mathbf{\hat{k}} \cdot\nabla \Phi
- 2\rho \mathbf{\hat{k}} \cdot [ \mathbf{\hat{e}}_\varphi \Omega_0 u'_R - \mathbf{\hat{e}}_R \Omega_0 u'_\varphi  ]
+ \rho \mathbf{\hat{k}}\cdot (\mathbf{\hat{e}}_R \Omega_0^2 R) \,
</math>
  </td>
</tr>
<tr>
  <td align="center">
&nbsp;
  </td>
  <td align="right">
<math>
\Rightarrow ~~~~~ \frac{\partial (\rho u'_z)}{\partial t} + \nabla\cdot[(\rho u'_z) \mathbf{u'} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi \, .
</math>
  </td>
</tr>
</table>
</table>
</div>
</div>
where, as noted above,
 
----
'''<font color="lightblue">ASIDE:</font>''' If we pause our discussion here and map this set of component equations onto a (rotating) cylindrical coordinate mesh &#8212; that is, if on the right-hand-sides we implement the straightforward operator projections,
<div align="center">
<math>\mathbf{\hat{e}}_R \cdot\nabla \rightarrow \frac{\partial}{\partial R} \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{e}}_\varphi \cdot\nabla \rightarrow \frac{1}{R} \frac{\partial}{\partial \varphi} \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{k}}\cdot\nabla \rightarrow \frac{\partial}{\partial z} \, ,</math>
</div>
we obtain a formulation that is familiar to the astrophysics community.  For example, as the following table of equations illustrates, it is the component set that has been spelled out in equations (5) - (7) of [http://adsabs.harvard.edu/abs/1978ApJ...224..497N Norman &amp; Wilson (1978)] and equations (11), (12), &amp; (3) of  [http://adsabs.harvard.edu/abs/1997ApJ...490..311N New &amp; Tohline (1997)].
 


<div align="center">
<div align="center">
<table border="1" cellpadding="5" width="80%">
<tr>
  <td align="center" bgcolor="lightgreen">
'''Cylindrical Components of the Rotating-Frame Momentum'''<br>
<font size="-1">advected across a</font><br>
'''Rotating, Cylindrical Coordinate Mesh'''
  </td>
</tr>
<tr><td align="center">
<table border="0" cellpadding="3">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u'_R)}{\partial t} + \nabla\cdot [(\rho u'_R)\mathbf{u'}]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \frac{\partial}{\partial R} P - \rho \frac{\partial}{\partial R}\Phi
+ \frac{\rho}{R} \biggl[ (u'_\varphi)^2 + 2R\Omega_0 u'_\varphi + (\Omega_0 R)^2 \biggr] \,
</math>
  </td>
</tr>
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>
<math>
\nabla\cdot[\psi_{i} \vec{v} ]
\frac{\partial (\rho R u'_\varphi)}{\partial t} + \nabla\cdot [(\rho R u'_\varphi) \mathbf{u'} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \frac{\partial}{\partial\varphi} P - \rho \frac{\partial}{\partial\varphi}\Phi
- 2\rho R \Omega_0 u'_R 
</math>
</math>
  </td>
</tr>
<tr>
  <td align="right">
<math>
\frac{\partial (\rho u'_z)  }{\partial t} + \nabla\cdot [(\rho u'_z) \mathbf{u'} ]
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 1,033: Line 1,764:
   <td align="left">
   <td align="left">
<math>
<math>
\frac{\partial (\psi_i v_R)}{\partial R} + \frac{1}{R} \frac{\partial (\psi_i v_\phi)}{\partial\phi} + \frac{\partial (\psi_i v_z)}{\partial z} \, .
- \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z}\Phi 
</math>
</math>
   </td>
   </td>
</tr>
</tr>
</table>
</td></tr>
</table>
</div>
<div align="center" id="NormanWilson78">
<table border="1" cellpadding="5" width="80%">
<tr><td align="center">
Paragraph extracted without modification from p. 499 of [http://adsabs.harvard.edu/abs/1978ApJ...224..497N M. L. Norman &amp; J. R. Wilson (1978)]<p></p>
"''The Fragmentation of Isothermal Rings and Star Formation''"<p></p>
ApJ, vol. 224, pp. 497-511 &copy; [http://aas.org/ American Astronomical Society]
</td></tr>
<tr>
  <td align="center">
<!-- [[File:NormanWilson78D.png|700px|center|Norman &amp; Wilson (1978)]]  -->
[[Image:AAAwaiting01.png|600px|center|Norman &amp; Wilson (1978)]]
  </td>
</tr>
<tr>
  <td align="left">
Note:  For complete correspondence with equations derived herein, set <math>(\gamma-1)E \rightarrow P</math> in all three component equations.
  </td>
</tr>
</table>
</div>
<div align="center" id="NewTohline97">
<table border="1" cellpadding="5" width="80%">
<tr><td align="center">
Equations extracted<sup>&dagger;</sup> from [http://adsabs.harvard.edu/abs/1997ApJ...490..311N K. C. B. New &amp; J. E. Tohline (1997)]<p></p>
"''The Relative Stability Against Merger of Close, Compact Binaries''"<p></p>
ApJ, vol. 490, pp. 311-327 &copy; [http://aas.org/ American Astronomical Society]
</td></tr>
<tr>
  <td align="left">
[[File:NewTohline97C.png|700px|center|New &amp; Tohline (1997)]]
  </td>
</tr>
<tr><td align="left">
<sup>&dagger;</sup>Mathematical expressions displayed here, as a single digital image, with presentation order &amp; layout modified from the original publication.<p></p>
Note:  When comparing this set of equations to the set presented by Norman &amp; Wilson (1978), the definitions of the variables, <math>~\mathcal{S}</math> and <math>~\mathcal{T}</math>, must be swapped.
</td></tr>
</table>
</table>
</div>
</div>


===Example #3===


[<font color="red">Comment by J. E. Tohline (April 7, 2014)</font>]  This is the set of equations that my research group has been using to simulate a wide variety of astrophysical fluid flows over the past twenty years. This is no longer our method of choice, however.  A numerical algorithm based on the hybrid scheme, as summarized below, is far preferable to an algorithm that is based on this more familiar, traditional set of equations for several reasons:
* In the hybrid scheme, the Coriolis term disappears from the source term, so it is much easier to design and implement a computational algorithm that conserves angular momentum conservation.
* Although the hybrid scheme advects ''inertial-frame'' angular momentum, it retains all of the advantages associated with using a ''rotating'' frame of reference; for example, numerical diffusion is less severe and, in general, the Courant-limited time-step is larger than would be the case if you were forced to transport fluid in the inertial frame of reference.
* The hybrid scheme facilitates transport (and conservation) of angular momentum across a (rotating) ''Cartesian'' mesh.  This facilitates the use of adaptive-mesh refinement (AMR) techniques and simplifies load-balancing on distributed memory, high-performance computers.
----
====Step 2====
Next, throughout this set of scalar equations, we replace each component of <math>~\rho\mathbf{u'}</math> with the corresponding component of <math>~(\rho\mathbf{u} - \rho {\vec{\Omega}}_f \times \vec{x})</math><math>= (\rho\mathbf{u} -\mathbf{\hat{e}}_\varphi \rho R\Omega_0)</math>, that is, we perform the following mappings:
<div align="center">
<div align="center">
<table border="0" cellpadding="3">
<table border="0" cellpadding="3">
Line 1,047: Line 1,830:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~\boldsymbol{\hat{e}}_R:</math>
<math>
~\rho u'_R
</math>
   </td>
   </td>
  <td align="center">
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
~\rho u_R \, ,
</math>
  </td>
<tr>
  <td align="right">
<math>
~\rho u'_\varphi
</math>
  </td>
  <td align="center">
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
~\rho (u_\varphi - R\Omega_0 ) \, ,
</math>
  </td>
</tr>
<tr>
  <td align="right">
<math>
~\rho u'_z
</math>
  </td>
  <td align="center">
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
~\rho u_z \, .
</math>
  </td>
</tr>
</table>
</div>
As a result, the first and third of the three "hybrid" momentum-component equations become, respectively,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="center">
<math>\mathbf{\hat{e}}_R:</math>
  </td>
  <td align="right">
<math>
\frac{\partial (\rho u_R)}{\partial t} + \nabla\cdot [(\rho u_R)\mathbf{u'}]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{e}}_R \cdot\nabla P - \rho \mathbf{\hat{e}}_R \cdot\nabla \Phi
+ \frac{\rho u^2_\varphi}{R}  \, ;
</math>
  </td>
</tr>
<tr>
  <td align="center">
<math>\mathbf{\hat{k}}:</math>
  </td>
  <td align="right">
<math>
\frac{\partial (\rho u_z)  }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi  \, .
</math>
  </td>
</tr>
</table>
</div>
The second of the three "hybrid" momentum component equations &#8212; the one governing conservation of angular momentum &#8212; becomes,
<div align="center">
<table border="0" cellpadding="3">
<tr>
   <td align="right">
   <td align="right">
<math>~\frac{\partial (\rho u'_R)}{\partial t} + \nabla\cdot[\rho u'_R (\vec{v})'~]</math>
<math>
\frac{\partial [\rho R (u_\varphi - R\Omega_0 ) ]}{\partial t} + \nabla\cdot \{ [\rho R (u_\varphi - R\Omega_0 ) ] \mathbf{u'} \}
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 1,057: Line 1,936:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial R} - \rho \frac{\partial \Phi}{\partial R} + \frac{\rho}{R} (v'_\phi + R\Omega_0)^2
- \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi
- 2\rho R \Omega_0 u_R    \,
</math>
</math>
   </td>
   </td>
</tr>
</tr>
</table>
</div>
<div align="center">
<table border="0" cellpadding="3">


<tr>
<tr>
   <td align="right">
   <td align="right">
&nbsp;
<math>
\Rightarrow ~~~~~ \frac{\partial (\rho R u_\varphi)  }{\partial t} + \nabla\cdot [(\rho R u_\varphi) \mathbf{u'} ]
- R^2 \Omega_0 \biggl[ \frac{\partial \rho }{\partial t} + \nabla\cdot ( \rho  \mathbf{u'} ) \biggr]
- ( \rho  \mathbf{u'} )\cdot\nabla (R^2 \Omega_0)
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi
- 2\rho R \Omega_0 u_R    \, .
</math>
   </td>
   </td>
</tr>
</table>
</div>
Referencing the continuity equation, as before, the middle bracketed term on the left-hand side can be set to zero; and the last term on the left-hand side,
<div align="center">
<table border="0" cellpadding="3">
<tr>
   <td align="right">
   <td align="right">
&nbsp;
<math>
( \rho  \mathbf{u'} )\cdot\nabla (R^2\Omega_0)
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 1,074: Line 1,982:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial R} - \rho \frac{\partial \Phi}{\partial R} +
2\rho R \Omega_0 u_R \, ,
\frac{\rho (v'_\phi)^2}{R} + 2\rho \Omega_0 v'_\phi + \rho \Omega_0^2 R \, ,
</math>
</math>
   </td>
   </td>
</tr>
</tr>
</table>
</div>
matches and, hence, exactly cancels the Coriolis term on the right-hand side to give,


<div align="center">
<table border="0" cellpadding="3">
<tr>
<tr>
  <td align="center">
<math>\mathbf{\hat{e}}_\varphi:</math>
  </td>
   <td align="right">
   <td align="right">
<math>~\boldsymbol{\hat{e}}_\phi:</math>
<math>
\frac{\partial (\rho R u_\varphi)  }{\partial t} + \nabla\cdot [(\rho R u_\varphi) \mathbf{u'} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi  \, .
</math>
  </td>
</tr>
</table>
</div>
 
====Step 3====
 
'''<font color="red">Cylindrical Grid:</font>''' If the numerical simulation is to be conducted on a cylindrical coordinate mesh, the spatial operators on both sides of the component momentum equations should be broken down into their cylindrical-coordinate components.  That is, as in the above "<font color="lightblue">ASIDE</font>," the appropriate operator projections on the right-hand-sides of the equations are,
<div align="center">
<math>\mathbf{\hat{e}}_R \cdot\nabla \rightarrow \frac{\partial}{\partial R} \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{e}}_\varphi \cdot\nabla \rightarrow \frac{1}{R} \frac{\partial}{\partial \varphi} \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{k}}\cdot\nabla \rightarrow \frac{\partial}{\partial z} \, .</math>
</div>
In concert with this, the divergence (advection) term on the left-hand-side should be evaluated by breaking the transport velocity into its three cylindrical components,
<div align="center">
<math>~\mathbf{u'} = (u'_R, u'_\varphi, u'_z) \, .</math>
</div>
The relevant set of momentum conservation equations is, therefore,
 
<div align="center">
<table border="1" cellpadding="5">
<tr>
  <td align="center" bgcolor="lightgreen">
'''Cylindrical Components of the Inertial-Frame Momentum'''<br>
<font size="-1">advected across a</font><br>
'''Rotating, Cylindrical Coordinate Mesh'''
   </td>
   </td>
</tr>
<tr><td align="center">
<table border="0" cellpadding="3">
<tr>
   <td align="right">
   <td align="right">
<math>~\frac{\partial \{\rho R [u'_\phi + R(\Omega_0 - \omega_0)]\} }{\partial t} +  
<math>
\nabla\cdot[ \{ \rho R [u'_\phi + R(\Omega_0 - \omega_0)] \} (\vec{v})'~]</math>
\frac{\partial (\rho u_R)}{\partial t} + \nabla\cdot [(\rho u_R)\mathbf{u'}]
</math>  
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 1,093: Line 2,050:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial \phi} - \rho \frac{\partial \Phi}{\partial \phi} - 2\rho R\omega_0 v'_R \, ,
- \frac{\partial}{\partial R} P - \rho \frac{\partial}{\partial R} \Phi  + \frac{\rho u^2_\varphi}{R}
</math>
</math>
   </td>
   </td>
Line 1,100: Line 2,057:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~\boldsymbol{\hat{e}}_z:</math>
<math>
\frac{\partial (\rho R u_\varphi)  }{\partial t} + \nabla\cdot [(\rho R u_\varphi) \mathbf{u'} ]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \frac{\partial}{\partial \varphi} P - \rho \frac{\partial}{\partial \varphi} \Phi 
</math>
   </td>
   </td>
</tr>
<tr>
   <td align="right">
   <td align="right">
<math>~\frac{\partial (\rho u'_z)}{\partial t} + \nabla\cdot[\rho u'_z (\vec{v})'~]</math>
<math>
\frac{\partial (\rho u_z) }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ]
</math>  
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 1,110: Line 2,082:
   <td align="left">
   <td align="left">
<math>
<math>
-~\frac{\partial P}{\partial z} - \rho \frac{\partial \Phi}{\partial z} \, ,
- \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z}\Phi
</math>
</math>
   </td>
   </td>
</tr>
</tr>
</table>
</table>
</td></tr>
</table>
</div>
'''<font color="red">Cartesian Grid:</font>''' If, instead, the numerical simulation is to be conducted on a Cartesian coordinate mesh, the spatial operators on both sides of the component momentum equations should be broken down into their Cartesian-coordinate components.  In concert with this, the divergence (advection) term on the left-hand-side should be evaluated by breaking the transport velocity into its three Cartesian components,
<div align="center">
<math>~\mathbf{u'} = (u'_x, u'_y, u'_z) \, .</math>
</div>
Furthermore, recognizing that, when written in Cartesian coordinates, the gradient operator is,
<div align="center">
<math>
\nabla = \mathbf{\hat{i}} \frac{\partial}{\partial x} +
\mathbf{\hat{j}} \frac{\partial}{\partial y} +\mathbf{\hat{k}} \frac{\partial}{\partial z} \, ,
</math>
</div>
</div>
where, as noted above,
and that the unit vectors in Cartesian coordinates can be related to their cylindrical counterparts via the mappings,
<div align="center">
<math>\mathbf{\hat{i}} = \mathbf{\hat{e}}_R \cos\varphi - \mathbf{\hat{e}}_\varphi \sin\varphi \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{j}} = \mathbf{\hat{e}}_R \sin\varphi + \mathbf{\hat{e}}_\varphi \cos\varphi \, ,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\mathbf{\hat{k}} = \mathbf{\hat{k}}  \, ,</math><br>
</div>
the relevant projections of the gradient operator on the right-hand-sides of the governing equations should take the form,
<div align="center">
<div align="center">
<table border="0" cellpadding="3">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>
\mathbf{\hat{e}}_R \cdot\nabla
</math>
  </td>
  <td align="center">
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
\biggl[ \cos\varphi \frac{\partial}{\partial x} +\sin\varphi \frac{\partial}{\partial y} \biggr]
= \biggl[ \frac{x}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial x} +\frac{y}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial y} \biggr]  \, ,
</math>
  </td>
</tr>
<tr>
  <td align="right">
<math>
R\mathbf{\hat{e}}_\varphi \cdot\nabla
</math>
  </td>
  <td align="center">
<math>~\rightarrow~</math>
  </td>
  <td align="left">
<math>
\biggl[ - R\sin\varphi \frac{\partial}{\partial x} +R\cos\varphi \frac{\partial}{\partial y} \biggr]
= \biggl[ - y \frac{\partial}{\partial x} +x \frac{\partial}{\partial y} \biggr]  \, ,
</math>
  </td>
</tr>
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~u'_\phi</math>
<math>
\mathbf{\hat{k}}\cdot\nabla
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~=~</math>
<math>~\rightarrow~</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~u_\phi - R\omega_0 \, ,</math>
<math>
\frac{\partial}{\partial z} \, .
</math>
   </td>
   </td>
</tr>
</tr>
Line 1,133: Line 2,166:
</div>
</div>


and, for any one of the three scalar PDEs, advection of the relevant component of the momentum density, <math>~\psi_i</math>, is handled via the operation,
In summary, then, the relevant set of momentum conservation equations is,
 
<div align="center">
<div align="center">
<table border="1" cellpadding="5">
<tr>
  <td align="center" bgcolor="lightgreen">
'''Cylindrical Components of the Inertial-Frame Momentum'''<br>
<font size="-1">advected across a</font><br>
'''Rotating, Cartesian Coordinate Mesh'''
  </td>
</tr>
<tr><td align="center">
<table border="0" cellpadding="3">
<table border="0" cellpadding="3">
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>
<math>
\nabla\cdot[\psi_{i} (\vec{v})' ]  
\frac{\partial (\rho u_R)}{\partial t} + \nabla\cdot [(\rho u_R)\mathbf{u'}]
</math>
  </td>
  <td align="center">
<math>~=~</math>
  </td>
  <td align="left">
<math>
- \biggl[ \frac{x}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial x} +\frac{y}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial y} \biggr] P
- \rho \biggl[ \frac{x}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial x} +\frac{y}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial y} \biggr] \Phi
+ \frac{\rho u^2_\varphi}{R} 
</math>
</math>
  </td>
</tr>
<tr>
  <td align="right">
<math>
\frac{\partial (\rho R u_\varphi)  }{\partial t} + \nabla\cdot [(\rho R u_\varphi) \mathbf{u'} ]
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 1,147: Line 2,209:
   <td align="left">
   <td align="left">
<math>
<math>
\frac{\partial (\psi_i v'_R)}{\partial R} + \frac{1}{R} \frac{\partial (\psi_i v'_\phi)}{\partial\phi} + \frac{\partial (\psi_i v'_z)}{\partial z}  
\biggl[ y \frac{\partial}{\partial x} - x \frac{\partial}{\partial y} \biggr]  P
+ \rho \biggl[ y \frac{\partial}{\partial x} - x \frac{\partial}{\partial y} \biggr] \Phi 
</math>
</math>
   </td>
   </td>
Line 1,154: Line 2,217:
<tr>
<tr>
   <td align="right">
   <td align="right">
&nbsp;
<math>
\frac{\partial (\rho u_z)  }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ]
</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
Line 1,161: Line 2,226:
   <td align="left">
   <td align="left">
<math>
<math>
\frac{\partial (\psi_i v_R)}{\partial R} + \frac{1}{R} \frac{\partial [\psi_i (v_\phi - R\Omega_0)]}{\partial\phi} + \frac{\partial (\psi_i v_z)}{\partial z} \, .
- \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z} \Phi 
</math>
</math>
  </td>
</tr>
</table>
</td></tr>
<tr>
  <td align="left">
This is the set of equations that has served as the foundation of the so-called ''Hybrid'' simulations reported in [http://adsabs.harvard.edu/abs/2014ApJS..212...23B Byerly, Adelstein-Lelbach, Tohline, &amp; Marcello (2014)].
   </td>
   </td>
</tr>
</tr>
Line 1,169: Line 2,242:


=Related Discussions=
=Related Discussions=
* Euler equation viewed from a [http://www.vistrails.org/index.php?title=User:Tohline/PGE/RotatingFrame rotating frame of reference] or [[Main Page]].
* [http://etd.lsu.edu/docs/available/etd-07052010-224427/ Jay Call's dissertation at LSU]
* An [http://www.vistrails.org/index.php/User:Tohline/PGE/ConservingMomentum earlier draft of this "Euler equation" presentation].
* Euler equation viewed from a [[User:Tohline/PGE/RotatingFrame|rotating frame of reference]], or [[Main Page]].
 
* An [[User:Tohline/PGE/ConservingMomentum|earlier draft of this "Euler equation" presentation]].
* Early draft of Hybrid Scheme discussion:  [[User:Tohline/Appendix/Ramblings/Hybrid_Scheme_old|Appendix/Ramblings/Hybrid Scheme Old]]


{{LSU_HBook_footer}}
{{LSU_HBook_footer}}

Latest revision as of 21:58, 11 March 2019

Hybrid Advection Scheme

Preface

March 1, 2014 by Joel E. Tohline

<full preface> … This seemed too good to be true. The discovered code modification would allow us to conserve angular momentum very accurately and, at the same time, allow us to use a rotating grid and thereby minimize numerical diffusion … <read more>

In his dissertation research, Jay Call (see Call, Tohline, & Lehner 2010) derived a complete description of this hybrid advection scheme in a fully relativistic and generalized coordinate framework. He showed that one can write the system of fluid equations in a manner that facilitates advection of inertial-frame quantities across a rotating grid. In addition — and more importantly — he showed how to write the system of fluid equations to allow advection of inertial-frame angular momentum (generally associated with a cylindrical coordinate mesh) across a rotating Cartesian grid. Since that time, primarily in the context of Zach Byerly's doctoral dissertation research, my group has demonstrated that this hybrid advection scheme works extremely well (see Byerly, Adelstein-Lelbach, Tohline, & Marcello (2014)). In the discussion that follows here, we derive the Newtonian version of Jay's hybrid scheme.


Whitworth's (1981) Isothermal Free-Energy Surface
|   Tiled Menu   |   Tables of Content   |  Banner Video   |  Tohline Home Page   |

Setting the Stage

Recognizing Statements of Conservation

When dealing with the compressible fluid equations, we will often encounter hyperbolic PDEs of the following form:

<math> \frac{d\psi}{dt} + \psi \nabla\cdot \vec{v} </math>

<math>~=~</math>

<math> S \, , </math>

where we are using <math>~\vec{v}</math> to represent the velocity field of the fluid as viewed from an inertial frame of reference, and the total (as opposed to partial) time derivative indicates the time-rate of change of <math>~\psi</math> is being measured in a so-called Lagrangian fashion, that is, at the location of some fluid element and moving along with that fluid element.

When we encounter a situation in which the "source" term, <math>~S</math>, on the right-hand side is zero, we will be able to identify the scalar variable, <math>~\psi</math>, as the volume density of some conserved quantity. For example, the continuity equation — which is a mathematical statement of mass conservation — has the form,

LSU Key.png

<math>\frac{d\rho}{dt} + \rho \nabla \cdot \vec{v} = 0</math>

      or, equivalently,      

<math> \frac{d\ln\rho}{dt} </math>

<math>~=~</math>

<math> ~- \nabla\cdot \vec{v} \, , </math>

where, <math>~\rho</math> is the mass per unit volume or, simply, the mass density of the fluid element. Clearly, when the mass of a Lagrangian fluid element is conserved, the fluid element's mass density changes only in accordance with the divergence of the local velocity field.

Similarly, if we are following the evolution of a fluid that expands and contracts adiabatically, we should expect to encounter an equation of the form,

<math> \frac{ds}{dt} + s\nabla\cdot \vec{v} </math>

<math>~=~</math>

<math> 0 \, , </math>

      or, equivalently,      

<math> \frac{d\ln s}{dt} </math>

<math>~=~</math>

<math> ~- \nabla\cdot \vec{v} \, , </math>

where, <math>~s</math> is the entropy density of a Lagrangian fluid element. Or, if an axisymmetric distribution of fluid is moving in an axisymmetric potential, we should expect the azimuthal component of the fluid's angular momentum to be conserved, in which case we should expect to encounter a dynamical equation of the form,

<math> \frac{d(\rho \varpi v_\phi)}{dt} + (\rho \varpi v_\phi) \nabla\cdot \vec{v} </math>

<math>~=~</math>

<math> 0 \, , </math>

where, <math>~\varpi</math> is the Lagrangian fluid element's (cylindrical radial) distance measured from the symmetry axis of the underlying potential and <math>~v_\phi = \varpi\dot\phi</math> is the azimuthal component of the inertial velocity field, <math>~\vec{v}</math>, at the location of the fluid element.

Alternative Reference Frames

Now, we might want to examine the time-dependent behavior of a fluid system while viewing the flow from a reference frame that is more or less moving along with the fluid. This new frame of reference need not be an inertial frame; for example, when studying a rotating fluid, we may want to view the system's evolution from a rotating frame of reference. This will be accomplished mathematically by adjusting the dynamical equations so that the velocity that appears in the divergence term accounts for the new "frame" velocity field; specifically, we want to replace <math>~\vec{v}</math> with,

<math> \vec{u} </math>

<math>~=~</math>

<math> \vec{v} - \vec{v}_\mathrm{frame} \, . </math>

(Here, we will consider only time-independent functional expressions for the frame velocity, <math>~\vec{v}_\mathrm{frame}</math>.) Of course, switching to the rotating frame must be done in such a way that the resulting, new PDE describes exactly the same physical behavior of the system as was described by the original equation; that is, the new equation must be derivable from the original one.

If <math>~\vec{v}_\mathrm{frame}</math> is a divergence-free velocity field, then the transformation is trivial. For example, if we choose a frame of reference that is rotating uniformly with angular velocity, <math>~\Omega_0</math>, then,

<math> \vec{v}_\mathrm{frame} </math>

<math>~=~</math>

<math> \boldsymbol{\hat{e}}_\phi (\varpi \Omega_0) \, , </math>

and, utilizing cylindrical coordinates,

<math> \nabla\cdot\vec{v}_\mathrm{frame} </math>

<math>~=~</math>

<math> \frac{\partial(0)}{\partial \varpi} + \frac{1}{\varpi}\frac{\partial(\varpi \Omega_0)}{\partial \phi} + \frac{\partial(0)}{\partial z} = 0 \, . </math>

Hence,

<math> \frac{d\psi}{dt} + \psi \nabla\cdot \vec{u} </math>

<math>~=~</math>

<math> \frac{d\psi}{dt} + \psi \nabla\cdot [\vec{v} - \vec{v}_\mathrm{frame}] </math>

<math>~=~</math>

<math> \frac{d\psi}{dt} + \psi \nabla\cdot \vec{v} \, , </math>

so the new generic hyperbolic PDE becomes,

<math> \frac{d\psi}{dt} + \psi \nabla\cdot \vec{u} </math>

<math>~=~</math>

<math> S \, , </math>

and we can be confident that this new PDE represents the physics of the system just as well as the original PDE.

Eulerian Representation

We can shift any of the PDEs from a Lagrangian to an Eulerian representation — and thereby use them to follow the time-rate of change of physical variables at a point in space that is fixed with respect to the chosen frame of reference — by using the following transformation to replace each total time derivative with a partial time derivative:

<math> \frac{d\psi}{dt} </math>

<math>~~~\rightarrow~~~</math>

<math> \frac{\partial \psi}{\partial t} + \vec{u} \cdot \nabla\psi \, . </math>

Hence, the "new" generic hyperbolic PDE derived above can be rewritten as,

<math> \frac{\partial\psi}{\partial t} + \vec{u} \cdot \nabla\psi + \psi \nabla\cdot \vec{u} </math>

<math>~=~</math>

<math> S \, , </math>

or, more succinctly,

<math> \frac{\partial\psi}{\partial t} + \nabla\cdot (\psi \vec{u} ) </math>

<math>~=~</math>

<math> S \, . </math>

This equation also is broadly recognized as a conservation statement because, when <math>~S = 0</math>, the variable <math>~\psi</math> will represent the volume density of a conserved quantity.

We should emphasize that the inertial-frame version of this Eulerian conservation equation can be retrieved straightforwardly by setting <math>~\Omega_0 = 0</math>, which is equivalent to setting <math>~\vec{u} = \vec{v}</math>. It is,

<math> \frac{\partial\psi}{\partial t} + \nabla\cdot (\psi \vec{v}) </math>

<math>~=~</math>

<math> S \, . </math>

The physics of the flow that is being described by this PDE is identical to the physics that is described by the preceding PDE. But an important distinction must be made regarding how the two equations are interpreted. The "inertial frame" version of the equation (explicitly containing <math>~\vec{v}</math>) provides the time-rate of change of <math>~\psi</math> at a fixed point in inertial space, while the "new" version (explicitly containing <math>~\vec{u}</math>) provides the time-rate of change of <math>~\psi</math> at a fixed point in our "new" rotating coordinate frame.

Angular Momentum Conservation

When the three vector components of the Euler equation (of motion) are projected onto a nonrotating cylindrical coordinate grid, the azimuthal component of the Euler equation may be written as,

<math> \frac{d(\rho \varpi v_\phi)}{dt} + (\rho \varpi v_\phi) \nabla\cdot \vec{v} </math>

<math>~=~</math>

<math> -\frac{\partial P}{\partial\phi} - \rho \frac{\partial\Phi}{\partial\phi} \, . </math>

For this equation, the source term is identified as,

<math> ~S </math>

<math>~=~</math>

<math> -\frac{\partial P}{\partial\phi} - \rho \frac{\partial\Phi}{\partial\phi} \, , </math>

and <math>~\psi = (\rho\varpi v_\phi)</math> is the inertial-frame angular momentum density, as measured with respect to the <math>~z</math>-coordinate axis. This corresponds the scalar equation and representation referred to as "Case B (<math>~\eta=3</math>)" in CTL (2010).

From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010)
Case B <math>~(\eta = 3)</math>
with the following replacements:    <math>~(\rho h)_\mathrm{CTL} \rightarrow \rho</math>   ;    <math>~(R)_\mathrm{CTL} \rightarrow \varpi</math>  ;    <math>~(R u^\phi)_\mathrm{CTL} \rightarrow \varpi\dot\phi = v_\phi</math>

<math>~\psi_{(3)}</math>

<math>~S_{(3)}</math>

<math>~\rho \varpi v_\phi</math>

<math>~ - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi}</math>


As foreshadowed above — see the subsection titled, Recognizing Statements of Conservation — the angular momentum of a Lagrangian fluid element will be conserved if the "source" term, <math>~S = 0</math>. This situation will arise if, at the fluid element's location, the azimuthal pressure variation, <math>~\partial P/\partial\phi</math>, and the azimuthal variation in the gravitational potential, <math>~\partial \Phi/\partial\phi</math>, are both zero, or if the two balance one another (i.e.,<math>~\partial P/\partial\phi=-\rho\partial\Phi/\partial\phi</math>).

Based on the above discussion, we can equally well view the flow from a frame of reference that is rotating with a constant angular velocity, <math>~\Omega_0</math>, and write,

<math> \frac{d(\rho \varpi v_\phi)}{dt} + (\rho \varpi v_\phi) \nabla\cdot \vec{u} </math>

<math>~=~</math>

<math> -\frac{\partial P}{\partial\phi} - \rho \frac{\partial\Phi}{\partial\phi} \, , </math>

where, as before,

<math> \vec{u} </math>

<math>~\equiv~</math>

<math> \vec{v} - \boldsymbol{\hat{e}}_\phi \varpi\Omega_0 \, . </math>

Also, following the earlier discussion, if one wants to follow the time-variation of the fluid's inertial-frame angular momentum at a fixed location in inertial space, then the appropriate Eulerian representation of this azimuthal component of the equation of motion is,

<math> \frac{\partial (\rho \varpi v_\phi)}{\partial t} + \nabla\cdot [(\rho \varpi v_\phi) \vec{v}] </math>

<math>~=~</math>

<math> S \, . </math>

If, however, one wants to follow the time-variation of the fluid's inertial-frame angular momentum at a fixed location on a rotating coordinate grid, then the appropriate Eulerian representation of this azimuthal component of the equation of motion is obtained by replacing the "transport" velocity, <math>~\vec{v}</math> with <math>~\vec{u}</math>; specifically,

<math> \frac{\partial (\rho \varpi v_\phi)}{\partial t} + \nabla\cdot [(\rho \varpi v_\phi) \vec{u}] </math>

<math>~=~</math>

<math> S \, . </math>

An Element of the Hybrid Scheme

This last equation displays one subtle, but valuable, element of the hybrid scheme developed by Call, Tohline, & Lehner (2010). The velocity component, <math>~v_\phi</math>, that appears in the formulation of the relevant conserved quantity — the inertial-frame angular momentum density — is drawn from the velocity vector, <math>~\vec{v}</math>, which is different from the transport velocity vector, <math>~\vec{u}</math>, that defines the Eulerian frame from which the dynamical evolution of the system is being viewed. This equation is usually written, instead, in a form such that the angular momentum density is expressed in terms of the azimuthal component of the transport velocity; see, for example, equation (7) in Norman & Wilson (1978) and equation (12) in New & Tohline (1997). In this more familiar formulation, the momentum density and the transport velocity both directly refer to the same frame of reference. But, as a consequence, the source term is more complicated.

The more familiar formulation — including its modified source term — can be derived from our "hybrid" formulation by recognizing that,

<math> ~v_\phi </math>

<math>~=~</math>

<math> ~u_\phi + \varpi\Omega_0 \, . </math>

So we can write,

<math> \frac{\partial [\rho \varpi (u_\phi + \varpi\Omega_0 ) ]}{\partial t} + \nabla\cdot \{[\rho \varpi ( u_\phi + \varpi\Omega_0)] \vec{u} \} </math>

<math>~=~</math>

<math> ~S_{\phi i} \, , </math>

where, as shorthand, we have used,

<math> ~S_{\phi i} </math>

<math>~\equiv~</math>

<math> - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi} \, . </math>

This implies,

<math> \frac{\partial (\rho \varpi u_\phi )}{\partial t} + \nabla\cdot [ (\rho \varpi u_\phi) \vec{u} ] </math>

<math>~=~</math>

<math> S_{\phi i} - \frac{\partial [\rho \varpi (\varpi\Omega_0 ) ]}{\partial t} - \nabla\cdot \{[\rho \varpi (\varpi\Omega_0)] \vec{u} \} </math>

 

<math>~=~</math>

<math> S_{\phi i} - \varpi^2\Omega_0 \biggl\{ \frac{\partial \rho}{\partial t} + \nabla\cdot (\rho \vec{u} ) \biggr\} - \rho \vec{u}\cdot \nabla(\varpi^2 \Omega_0) </math>

 

<math>~=~</math>

<math> S_{\phi i} - 2\rho \varpi u_\varpi \Omega_0 \, . </math>

As we see, all terms involving the velocity now explicitly refer to <math>~\vec{u}</math> and, hence, to the velocity as measured in the rotating reference frame. But the source now includes a Coriolis term. This corresponds the scalar equation and representation referred to as "Case B (<math>~\eta=3'</math>)" in CTL (2010).

From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010)
Case B <math>~(\eta = 3')</math>
as before:    <math>~(\rho h)_\mathrm{CTL} \rightarrow \rho</math>   ;    <math>~(R)_\mathrm{CTL} \rightarrow \varpi</math>  ;    <math>~(R u^\phi)_\mathrm{CTL} \rightarrow \varpi\dot\phi = v_\phi</math>
additional replacements:    <math>~(\bar\omega u^{t'})_\mathrm{CTL} \rightarrow \Omega_0</math>  ;    <math>~u^R \rightarrow v_\varpi = u_\varpi</math>

<math>~\psi_{(3')}</math>

<math>~S_{(3')}</math>

<math>~\rho \varpi (v_\phi - \varpi\Omega_0) = \rho \varpi u_\phi </math>

<math>~ - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi} - 2\rho\varpi u_\varpi \Omega_0</math>

Even Broader Generalization

As Call, Tohline, & Lehner (2010) point out, we are free to measure — and follow the evolution of — the angular momentum density with respect to any of a variety of different rotating frames of reference. Specifically, we are not constrained to choose between the inertial (nonrotating) frame — in which the measured angular momentum density is <math>~(\rho\varpi v_\phi)</math> — and the "grid" frame — in which the measured angular momentum density is <math>~\rho\varpi u_\phi = \rho\varpi(v_\phi - \varpi\Omega_0)</math>. Quite generally, we can choose to measure the angular momentum with respect to a separate "primed" frame that is rotating with angular velocity <math>~\omega_0</math> and in which the measured azimuthal component of the fluid velocity is,

<math> ~v'_\phi ~= ~v_\phi - \varpi \omega_0 \, . </math>

With this definition in hand, we also recognize that,

<math> ~u_\phi = v_\phi - \varpi\Omega_0 = v'_\phi + \varpi (\omega_0 - \Omega_0) \, . </math>

These two substitutions allow us to rewrite the angular momentum evolution equation in the forms that Call, Tohline, & Lehner (2010) label as Case C (<math>~\eta=3</math>) and Case C (<math>~\eta=3'</math>).

From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010)
Case C <math>~(\eta = 3)</math>
with the following replacements:    <math>~(\rho h)_\mathrm{CTL} \rightarrow \rho</math>   ;    <math>~(R)_\mathrm{CTL} \rightarrow \varpi</math>  ;    <math>~(R u^\phi)_\mathrm{CTL} \rightarrow \varpi\dot\phi = v_\phi</math>

<math>~\psi_{(3)}</math>

<math>~S_{(3)}</math>

<math>~\rho \varpi (v'_\phi +\varpi\omega_0) = \rho \varpi v_\phi</math>

<math>~ - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi}</math>

 
 

From Tables 6.1 & 6.2 of Call, Tohline, & Lehner (2010)
Case C <math>~(\eta = 3')</math>
as before:    <math>~(\rho h)_\mathrm{CTL} \rightarrow \rho</math>   ;    <math>~(R)_\mathrm{CTL} \rightarrow \varpi</math>  ;    <math>~(R u^\phi)_\mathrm{CTL} \rightarrow \varpi\dot\phi = v_\phi</math>
additional replacements:    <math>~(\bar\omega u^{t'})_\mathrm{CTL} \rightarrow \Omega_0</math>  ;    <math>~u^R \rightarrow v_\varpi = u_\varpi</math>

<math>~\psi_{(3')}</math>

<math>~S_{(3')}</math>

<math>~\rho \varpi[ v'_\phi + \varpi (\omega_0 - \Omega_0)] = \rho \varpi u_\phi </math>

<math>~ - \frac{\partial P}{\partial\phi} - \rho \frac{\partial \Phi}{\partial\phi} - 2\rho\varpi u_\varpi \Omega_0</math>

In the latter case, <math>~v'_\phi</math> becomes <math>~u_\phi</math>, as it should, when the choice is made to measure the angular momentum density in the "grid" frame, that is, when the choice is made to set <math>~\omega_0 = \Omega_0</math>.

Zach's Dissertation Paper

Section 2.3

Building on our introductory discussion of the Euler equation (see also Appendix 1.D, §3 of BT87), we begin with the,

Lagrangian Representation
of the Euler Equation
as viewed from a Rotating Reference Frame

<math>\frac{d\mathbf{u'}}{dt} = - \frac{1}{\rho} \nabla P - \nabla \Phi - 2{\vec{\Omega}}_f \times {\mathbf{u'}} - {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x})</math> ,

where we choose to define the frame rotation by the vector,

<math>~\vec{\Omega}_f \equiv \hat\mathbf{k} \Omega_0 \, ,</math>

and,

<math> ~\mathbf{u'} </math>

<math>~=~</math>

<math> ~\mathbf{u} - \mathbf{u}_\mathrm{frame} = \mathbf{u} - {\vec{\Omega}}_f \times \vec{x}\, , </math>

is the velocity field as viewed from the frame of reference that is rotating at constant angular frequency, <math>~\Omega_0</math>. (Note that we can retrieve the inertial-frame Euler equation and inertial-frame variables by setting <math>~\Omega_0 = 0</math> at any point in the subsequent derivations.) Because the velocity field introduced by frame rotation is divergence free, that is, because,

<math> ~\nabla\cdot\mathbf{u}_\mathrm{frame} = \nabla\cdot [{\vec{\Omega}}_f \times \vec{x}] = \nabla\cdot [ {\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{e}}_R R + \mathbf{\hat{k}}z) ] </math>

<math>~=~</math>

<math> ~\nabla\cdot(\hat\mathbf{e}_\varphi R\Omega_0) = \frac{1}{R} \frac{\partial}{\partial\varphi} \biggl(R\Omega_0 \biggr) = 0 \, , </math>

the relevant (rotating-frame) continuity equation is identical in form to its inertial-frame counterpart, specifically,

<math> \frac{d\rho}{dt} + \rho\nabla\cdot\mathbf{u'} </math>

<math>~=~</math>

<math> ~0 \, . </math>

We can transform the above rotating-frame Euler equation to a momentum conservation equation by multiplying the equation through by <math>~\rho</math> and using the continuity equation to combine and simplify the left-hand-side, obtaining,

<math>\frac{d(\rho\mathbf{u'})}{dt} + (\rho\mathbf{u' })\nabla\cdot \mathbf{u'} = - \nabla P - \rho \nabla \Phi - 2\rho {\vec{\Omega}}_f \times {\mathbf{u'}} - \rho {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x})</math> ,

In order to convert this general-purpose vector equation into the specific set of scalar component equations that embody the desired elements of our hybrid scheme, we need to:

  • Step 1: Choose the unit-vector basis set associated with the momentum components that we want to track — <math>(\mathbf{\hat{i}}, \mathbf{\hat{j}}, \mathbf\hat{k})</math> if tracking linear momentum, or, <math>(\mathbf\hat{e}_R, \mathbf\hat{e}_\varphi, \mathbf\hat{k})</math> if tracking radial and angular momentum — and break the vector equation into these specified components;
  • Step 2: In all relevant equations, replace the scalar components of the "rotating-frame" momentum density with the scalar components of the "intertial-frame" momentum density, drawing each component relation from the vector transformation, <math>\rho\mathbf{u'} = \rho\mathbf{u} - \rho {\vec{\Omega}}_f \times \vec{x}</math>;
  • Step 3: Write all of the spatial operators in terms of spatial derivatives that are associated with the unit-vector basis set of the desired computational mesh.

Focus on Tracking Linear Momentum

Step 1

If the focus is on tracking linear momentum components, then we need to break the vector momentum equation into its <math>(\mathbf\hat{i}, \mathbf\hat{j}, \mathbf\hat{k})</math> components. This is done by, in turn, "dotting" each unit vector into the vector equation. It is straightforward once we appreciate that the orientation of these Cartesian unit vectors does not vary in space and that, within the context of the rotating frame on which they are defined, these unit vectors do not vary in time. Hence, the first term in the vector equation — the material time derivative — can be written as,

<math> \frac{d(\rho\mathbf{u'})}{dt} = \frac{d}{dt} [ \mathbf{\hat{i}} (\rho u'_x) + \mathbf{\hat{j}} (\rho u'_y) + \mathbf{\hat{k}} (\rho u'_z) ] = \mathbf{\hat{i}} \frac{d(\rho u'_x)}{dt} + \mathbf{\hat{j}} \frac{d(\rho u'_y)}{dt} + \mathbf{\hat{k}} \frac{d(\rho u'_z)}{dt} \, , </math>

and the process of "dotting" each unit vector into the equation leads to the following set of scalar momentum-component equations:

<math>\mathbf{\hat{i}}:</math>

<math> \frac{d(\rho u'_x)}{dt} + (\rho u'_x)\nabla\cdot \mathbf{u'} </math>

<math>~=~</math>

<math> - \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + 2\rho \mathbf{\hat{i}}\cdot[\mathbf{\hat{i}}\Omega_0 u'_y - \mathbf{\hat{j}}\Omega_0 u'_x] + \rho \mathbf{\hat{i}}\cdot[ \mathbf{\hat{i}}\Omega_0^2 x + \mathbf{\hat{j}}\Omega_0^2 y] \, </math>

 

<math> \Rightarrow ~~~~~ \frac{\partial (\rho u'_x)}{\partial t} + \nabla\cdot [(\rho u'_x)\mathbf{u'}] </math>

<math>~=~</math>

<math> - \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + 2\rho \Omega_0 u'_y + \rho \Omega_0^2 x \, ; </math>

<math>\mathbf{\hat{j}}:</math>

<math> \frac{d(\rho u'_y)}{dt} + (\rho u'_y)\nabla\cdot \mathbf{u'} </math>

<math>~=~</math>

<math> - \mathbf{\hat{j}}\cdot\nabla P - \rho \mathbf{\hat{j}}\cdot\nabla \Phi + 2\rho \mathbf{\hat{j}}\cdot[\mathbf{\hat{i}}\Omega_0 u'_y - \mathbf{\hat{j}}\Omega_0 u'_x] + \rho \mathbf{\hat{j}}\cdot[ \mathbf{\hat{i}}\Omega_0^2 x + \mathbf{\hat{j}}\Omega_0^2 y] \, </math>

 

<math> \Rightarrow ~~~~~ \frac{\partial (\rho u'_y)}{\partial t} + \nabla\cdot [(\rho u'_y) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \mathbf{\hat{j}}\cdot\nabla P - \rho \mathbf{\hat{j}}\cdot\nabla \Phi - 2\rho \Omega_0 u'_x + \rho \Omega_0^2 y \, ; </math>

<math>\mathbf{\hat{k}}:</math>

<math> \frac{d(\rho u'_z)}{dt} + (\rho u'_z)\nabla\cdot \mathbf{u'} </math>

<math>~=~</math>

<math> - \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi + 2\rho \mathbf{\hat{k}}\cdot[\mathbf{\hat{i}}\Omega_0 u'_y - \mathbf{\hat{j}}\Omega_0 u'_x] + \rho \mathbf{\hat{k}}\cdot[ \mathbf{\hat{i}}\Omega_0^2 x + \mathbf{\hat{j}}\Omega_0^2 y] \, </math>

 

<math> \Rightarrow ~~~~~ \frac{\partial (\rho u'_z)}{\partial t} + \nabla\cdot[(\rho u'_z) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi \, . </math>

In deriving these expressions, we also (a) have recognized from the start that, when expressed in Cartesian coordinates,

<math> ~{\vec{\Omega}}_f \times \vec{x} </math>

<math>~=~</math>

<math> {\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{i}}x + \mathbf{\hat{j}}y + \mathbf{\hat{k}}z) = - \mathbf{\hat{i}}\Omega_0 y + \mathbf{\hat{j}}\Omega_0 x \, , </math>

<math> {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x}) </math>

<math>~=~</math>

<math> \hat{\mathbf{k}} \Omega_0 \times ( - \mathbf{\hat{i}}\Omega_0 y + \mathbf{\hat{j}}\Omega_0 x ) = - \mathbf{\hat{i}}\Omega_0^2 x - \mathbf{\hat{j}}\Omega_0^2 y \, , </math>

<math> {\vec{\Omega}}_f \times {\mathbf{u'}} </math>

<math>~=~</math>

<math> {\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{i}}u'_x + \mathbf{\hat{j}}u'_y + \mathbf{\hat{k}}u'_z) = - \mathbf{\hat{i}}\Omega_0 u'_y + \mathbf{\hat{j}}\Omega_0 u'_x\, , </math>

and (b) have used the familiar operator mapping,

<math>\frac{d}{dt} \rightarrow \frac{\partial}{\partial t} + \mathbf{u'}\cdot \nabla \, ,</math>

to shift from the total (Lagrangian) time derivative to the partial (Eulerian) time derivative, which is usually the more desirable representation for computational simulations.

Step 2

Next, throughout this set of scalar equations, we replace each component of <math>~\rho\mathbf{u'}</math> with the corresponding component of <math>~(\rho\mathbf{u} - \rho {\vec{\Omega}}_f \times \vec{x})</math>, that is, we perform the following mappings:

<math> ~\rho u'_x </math>

<math>~\rightarrow~</math>

<math> ~\rho (u_x +\Omega_0 y) \, , </math>

<math> ~\rho u'_y </math>

<math>~\rightarrow~</math>

<math> ~\rho (u_y - \Omega_0 x) \, , </math>

<math> ~\rho u'_z </math>

<math>~\rightarrow~</math>

<math> ~\rho u_z \, . </math>

As a result, the first of the three "hybrid" momentum-component equations becomes,

<math> \frac{\partial [ \rho (u_x +\Omega_0 y) ] }{\partial t} + \nabla\cdot \{ [\rho (u_x +\Omega_0 y) ]\mathbf{u'}\} </math>

<math>~=~</math>

<math> - \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + 2\rho \Omega_0 (u_y - \Omega_0 x) + \rho \Omega_0^2 x </math>

<math> \Rightarrow ~~~~~ \frac{\partial (\rho u_x) }{\partial t} + \nabla\cdot [(\rho u_x) \mathbf{u'} ] + \Omega_0 y\biggl[ \frac{\partial \rho }{\partial t} + \nabla\cdot ( \rho \mathbf{u'} ) \biggr] + ( \rho \mathbf{u'} )\cdot\nabla (\Omega_0 y) </math>

<math>~=~</math>

<math> - \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + 2\rho \Omega_0 u_y - \rho \Omega_0^2 x \, . </math>

Referencing the continuity equation, the middle bracketed term on the left-hand side can be set to zero; and the last term on the left-hand side,

<math> ( \rho \mathbf{u'} )\cdot\nabla (\Omega_0 y) </math>

<math>~=~</math>

<math> \rho \Omega_0 u'_y = \rho\Omega_0(u_y - \Omega_0 x) \, , </math>

can be combined with terms on the right-hand side — cutting the Coriolis term in half and canceling the centrifugal acceleration term — to give,

<math> \frac{\partial (\rho u_x) }{\partial t} + \nabla\cdot [(\rho u_x) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \mathbf{\hat{i}}\cdot\nabla P - \rho \mathbf{\hat{i}}\cdot\nabla \Phi + \rho \Omega_0 u_y \, . </math>

Following a similar sequence of steps, the other two "hybrid" momentum conservation relations become,

<math> \frac{\partial (\rho u_y) }{\partial t} + \nabla\cdot [(\rho u_y) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \mathbf{\hat{j}}\cdot\nabla P - \rho \mathbf{\hat{j}}\cdot\nabla \Phi - \rho \Omega_0 u_x \, , </math>

<math> \frac{\partial (\rho u_z) }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi \, . </math>

Step 3

Cartesian Grid: Assuming the numerical simulation will be conducted on a Cartesian coordinate mesh, the divergence (advection) term on the left-hand-side should be evaluated by breaking the transport velocity into its three Cartesian components,

<math>~\mathbf{u'} = (u'_x, u'_y, u'_z) \, ,</math>

and, on the right-hand-side, the projection of the spatial operators should be written in the familiar form,

<math>\mathbf{\hat{i}}\cdot\nabla \rightarrow \frac{\partial}{\partial x} \, ,</math>      <math>\mathbf{\hat{j}}\cdot\nabla \rightarrow \frac{\partial}{\partial y} \, ,</math>      <math>\mathbf{\hat{k}}\cdot\nabla \rightarrow \frac{\partial}{\partial z} \, .</math>

In summary, then, the relevant set of momentum conservation equations is,

Cartesian Components of the Inertial-Frame Momentum
advected across a
Rotating, Cartesian Coordinate Mesh

<math> \frac{\partial (\rho u_x) }{\partial t} + \nabla\cdot [(\rho u_x) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial x} P - \rho \frac{\partial}{\partial x} \Phi + \rho \Omega_0 u_y </math>

<math> \frac{\partial (\rho u_y) }{\partial t} + \nabla\cdot [(\rho u_y) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial y} P - \rho \frac{\partial}{\partial y} \Phi - \rho \Omega_0 u_x </math>

<math> \frac{\partial (\rho u_z) }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z}\Phi </math>

This is the set of equations that has served as the foundation of the Cartesian simulations reported in Byerly, Adelstein-Lelbach, Tohline, & Marcello (2014).


Cylindrical Grid: If, instead, the numerical simulation is to be conducted on a cylindrical coordinate mesh, the spatial operators on both sides of the component momentum equations should be broken down into their cylindrical-coordinate components. In concert with this, the divergence (advection) term on the left-hand-side should be evaluated by breaking the transport velocity into its three cylindrical components,

<math>~\mathbf{u'} = (u'_R, u'_\varphi, u'_z) \, .</math>

Furthermore, recognizing that, when written in cylindrical coordinates, the gradient operator is,

<math> \nabla = \mathbf{\hat{e}}_R \frac{\partial}{\partial R} + \mathbf{\hat{e}}_\varphi \frac{1}{R} \frac{\partial}{\partial \varphi} +\mathbf{\hat{e}}_z \frac{\partial}{\partial z} \, , </math>

and that the unit vectors in cylindrical coordinates can be related to their Cartesian counterparts via the mappings,

<math>\mathbf{\hat{e}}_R = \mathbf{\hat{i}} \cos\varphi + \mathbf{\hat{j}} \sin\varphi \, ,</math>      <math>\mathbf{\hat{e}}_\varphi = \mathbf{\hat{j}} \cos\varphi - \mathbf{\hat{i}} \sin\varphi \, ,</math>      <math>\mathbf{\hat{e}}_z = \mathbf{\hat{k}} \, ,</math>

the relevant projections of the gradient operator on the right-hand-sides of the governing equations should take the form,

<math> \mathbf{\hat{i}}\cdot\nabla </math>

<math>~\rightarrow~</math>

<math> \biggl[ \cos\varphi \frac{\partial}{\partial R} - \frac{\sin\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] \, , </math>

<math> \mathbf{\hat{j}}\cdot\nabla </math>

<math>~\rightarrow~</math>

<math> \biggl[ \sin\varphi \frac{\partial}{\partial R} + \frac{\cos\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] \, , </math>

<math> \mathbf{\hat{k}}\cdot\nabla </math>

<math>~\rightarrow~</math>

<math> \frac{\partial}{\partial z} \, . </math>

In summary, then, the relevant set of momentum conservation equations is,

Cartesian Components of the Inertial-Frame Momentum
advected across a
Rotating, Cylindrical Coordinate Mesh

<math> \frac{\partial (\rho u_x) }{\partial t} + \nabla\cdot [(\rho u_x) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \biggl[ \cos\varphi \frac{\partial}{\partial R} - \frac{\sin\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] P - \rho \biggl[ \cos\varphi \frac{\partial}{\partial R} - \frac{\sin\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] \Phi + \rho \Omega_0 u_y </math>

<math> \frac{\partial (\rho u_y) }{\partial t} + \nabla\cdot [(\rho u_y) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \biggl[ \sin\varphi \frac{\partial}{\partial R} + \frac{\cos\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] P - \rho \biggl[ \sin\varphi \frac{\partial}{\partial R} + \frac{\cos\varphi}{R} \frac{\partial}{\partial \varphi}\biggr] \Phi - \rho \Omega_0 u_x </math>

<math> \frac{\partial (\rho u_z) }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z}\Phi </math>

Focus on Tracking Angular Momentum

Step 1

If the focus is on tracking angular momentum, then we need to break the vector momentum equation into its <math>(\mathbf\hat{e}_R, \mathbf\hat{e}_\varphi, \mathbf\hat{k})</math> components. As before, this is done by "dotting" each unit vector into the vector equation. This is less straightforward than in the Cartesian case because the orientation of both the <math>\mathbf\hat{e}_R</math> and <math>\mathbf\hat{e}_\varphi</math> unit vectors vary in space. As a result, the first term in the vector equation — the material time derivative — generates a couple of extra terms, viz.,

<math> \frac{d(\rho\mathbf{u'})}{dt} </math>

<math>~=~</math>

<math> \frac{d}{dt} [ \mathbf{\hat{e}}_R (\rho u'_R) + \mathbf{\hat{e}}_\varphi (\rho u'_\varphi) + \mathbf{\hat{k}} (\rho u'_z) ] </math>

 

<math>~=~</math>

<math> \mathbf{\hat{e}}_R \frac{d(\rho u'_R)}{dt} + \mathbf{\hat{e}}_\varphi \frac{d(\rho u'_\varphi)}{dt} + \mathbf{\hat{k}} \frac{d(\rho u'_z)}{dt} + (\rho u'_R) \frac{d}{dt}\mathbf{\hat{e}}_R + (\rho u'_\varphi) \frac{d}{dt}\mathbf{\hat{e}}_\varphi \, , </math>

 

<math>~=~</math>

<math> \mathbf{\hat{e}}_R \frac{d(\rho u'_R)}{dt} + \mathbf{\hat{e}}_\varphi \frac{d(\rho u'_\varphi)}{dt} + \mathbf{\hat{k}} \frac{d(\rho u'_z)}{dt} + \mathbf{\hat{e}}_\varphi(\rho u'_R) \frac{u'_\varphi}{R} - \mathbf{\hat{e}}_R(\rho u'_\varphi) \frac{u'_\varphi}{R} \, . </math>

We also recognize that, when expressed in cylindrical coordinates,

<math> ~{\vec{\Omega}}_f \times \vec{x} </math>

<math>~=~</math>

<math> {\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{e}}_R R + \mathbf{\hat{k}}z) = \mathbf{\hat{e}}_\varphi \Omega_0 R \, , </math>

<math> {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x}) </math>

<math>~=~</math>

<math> \hat{\mathbf{k}} \Omega_0 \times ( \mathbf{\hat{e}}_\varphi \Omega_0 R ) = - \mathbf{\hat{e}}_R \Omega_0^2 R \, , </math>

<math> {\vec{\Omega}}_f \times {\mathbf{u'}} </math>

<math>~=~</math>

<math> {\hat\mathbf{k}} \Omega_0\times (\mathbf{\hat{e}}_R u'_R + \mathbf{\hat{e}}_\varphi u'_\varphi + \mathbf{\hat{k}}u'_z) = \mathbf{\hat{e}}_\varphi \Omega_0 u'_R - \mathbf{\hat{e}}_R \Omega_0 u'_\varphi \, . </math>

Hence, the process of "dotting" each unit vector into the equation leads to the following set of scalar momentum-component equations:

<math>\mathbf{\hat{e}}_R:</math>

<math> \frac{d(\rho u'_R)}{dt} + (\rho u'_R)\nabla\cdot \mathbf{u'} - (\rho u'_\varphi) \frac{u'_\varphi}{R} </math>

<math>~=~</math>

<math> - \mathbf{\hat{e}}_R \cdot\nabla P - \rho \mathbf{\hat{e}}_R \cdot\nabla \Phi - 2\rho \mathbf{\hat{e}}_R \cdot [ \mathbf{\hat{e}}_\varphi \Omega_0 u'_R - \mathbf{\hat{e}}_R \Omega_0 u'_\varphi ] + \rho \mathbf{\hat{e}}_R\cdot (\mathbf{\hat{e}}_R \Omega_0^2 R) \, </math>

 

<math> \Rightarrow ~~~~~ \frac{\partial (\rho u'_R)}{\partial t} + \nabla\cdot [(\rho u'_R)\mathbf{u'}] </math>

<math>~=~</math>

<math> - \mathbf{\hat{e}}_R \cdot\nabla P - \rho \mathbf{\hat{e}}_R \cdot\nabla \Phi + \frac{\rho}{R} \biggl[ (u'_\varphi)^2 + 2R\Omega_0 u'_\varphi + (\Omega_0 R)^2 \biggr] \, </math>

 

 

<math>~=~</math>

<math> - \mathbf{\hat{e}}_R \cdot\nabla P - \rho \mathbf{\hat{e}}_R \cdot\nabla \Phi + \frac{\rho}{R} (u'_\varphi + R\Omega_0)^2 \, ; </math>

<math>\mathbf{\hat{e}}_\varphi:</math>

<math> \frac{d(\rho u'_\varphi)}{dt} + (\rho u'_\varphi)\nabla\cdot \mathbf{u'} + (\rho u'_R) \frac{u'_\varphi}{R} </math>

<math>~=~</math>

<math> - \mathbf{\hat{e}}_\varphi \cdot\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot\nabla \Phi - 2\rho \mathbf{\hat{e}}_\varphi \cdot [ \mathbf{\hat{e}}_\varphi \Omega_0 u'_R - \mathbf{\hat{e}}_R \Omega_0 u'_\varphi ] + \rho \mathbf{\hat{e}}_\varphi\cdot (\mathbf{\hat{e}}_R \Omega_0^2 R) \, </math>

(mult. thru by R)

<math> \Rightarrow ~~~~~ \frac{d(\rho R u'_\varphi)}{dt} + (\rho R u'_\varphi)\nabla\cdot \mathbf{u'} </math>

<math>~=~</math>

<math> - \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi - 2\rho R \Omega_0 u'_R \, </math>

 

<math> \Rightarrow ~~~~~ \frac{\partial (\rho R u'_\varphi)}{\partial t} + \nabla\cdot [(\rho R u'_\varphi) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi - 2\rho R \Omega_0 u'_R \, ; </math>

<math>\mathbf{\hat{k}}:</math>

<math> \frac{d(\rho u'_z)}{dt} + (\rho u'_z)\nabla\cdot \mathbf{u'} </math>

<math>~=~</math>

<math> - \mathbf{\hat{k}} \cdot\nabla P - \rho \mathbf{\hat{k}} \cdot\nabla \Phi - 2\rho \mathbf{\hat{k}} \cdot [ \mathbf{\hat{e}}_\varphi \Omega_0 u'_R - \mathbf{\hat{e}}_R \Omega_0 u'_\varphi ] + \rho \mathbf{\hat{k}}\cdot (\mathbf{\hat{e}}_R \Omega_0^2 R) \, </math>

 

<math> \Rightarrow ~~~~~ \frac{\partial (\rho u'_z)}{\partial t} + \nabla\cdot[(\rho u'_z) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi \, . </math>


ASIDE: If we pause our discussion here and map this set of component equations onto a (rotating) cylindrical coordinate mesh — that is, if on the right-hand-sides we implement the straightforward operator projections,

<math>\mathbf{\hat{e}}_R \cdot\nabla \rightarrow \frac{\partial}{\partial R} \, ,</math>      <math>\mathbf{\hat{e}}_\varphi \cdot\nabla \rightarrow \frac{1}{R} \frac{\partial}{\partial \varphi} \, ,</math>      <math>\mathbf{\hat{k}}\cdot\nabla \rightarrow \frac{\partial}{\partial z} \, ,</math>

we obtain a formulation that is familiar to the astrophysics community. For example, as the following table of equations illustrates, it is the component set that has been spelled out in equations (5) - (7) of Norman & Wilson (1978) and equations (11), (12), & (3) of New & Tohline (1997).


Cylindrical Components of the Rotating-Frame Momentum
advected across a
Rotating, Cylindrical Coordinate Mesh

<math> \frac{\partial (\rho u'_R)}{\partial t} + \nabla\cdot [(\rho u'_R)\mathbf{u'}] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial R} P - \rho \frac{\partial}{\partial R}\Phi + \frac{\rho}{R} \biggl[ (u'_\varphi)^2 + 2R\Omega_0 u'_\varphi + (\Omega_0 R)^2 \biggr] \, </math>

<math> \frac{\partial (\rho R u'_\varphi)}{\partial t} + \nabla\cdot [(\rho R u'_\varphi) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial\varphi} P - \rho \frac{\partial}{\partial\varphi}\Phi - 2\rho R \Omega_0 u'_R </math>

<math> \frac{\partial (\rho u'_z) }{\partial t} + \nabla\cdot [(\rho u'_z) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z}\Phi </math>

Paragraph extracted without modification from p. 499 of M. L. Norman & J. R. Wilson (1978)

"The Fragmentation of Isothermal Rings and Star Formation"

ApJ, vol. 224, pp. 497-511 © American Astronomical Society

Norman & Wilson (1978)

Note: For complete correspondence with equations derived herein, set <math>(\gamma-1)E \rightarrow P</math> in all three component equations.

Equations extracted from K. C. B. New & J. E. Tohline (1997)

"The Relative Stability Against Merger of Close, Compact Binaries"

ApJ, vol. 490, pp. 311-327 © American Astronomical Society

New & Tohline (1997)
Mathematical expressions displayed here, as a single digital image, with presentation order & layout modified from the original publication.

Note: When comparing this set of equations to the set presented by Norman & Wilson (1978), the definitions of the variables, <math>~\mathcal{S}</math> and <math>~\mathcal{T}</math>, must be swapped.


[Comment by J. E. Tohline (April 7, 2014)] This is the set of equations that my research group has been using to simulate a wide variety of astrophysical fluid flows over the past twenty years. This is no longer our method of choice, however. A numerical algorithm based on the hybrid scheme, as summarized below, is far preferable to an algorithm that is based on this more familiar, traditional set of equations for several reasons:

  • In the hybrid scheme, the Coriolis term disappears from the source term, so it is much easier to design and implement a computational algorithm that conserves angular momentum conservation.
  • Although the hybrid scheme advects inertial-frame angular momentum, it retains all of the advantages associated with using a rotating frame of reference; for example, numerical diffusion is less severe and, in general, the Courant-limited time-step is larger than would be the case if you were forced to transport fluid in the inertial frame of reference.
  • The hybrid scheme facilitates transport (and conservation) of angular momentum across a (rotating) Cartesian mesh. This facilitates the use of adaptive-mesh refinement (AMR) techniques and simplifies load-balancing on distributed memory, high-performance computers.

Step 2

Next, throughout this set of scalar equations, we replace each component of <math>~\rho\mathbf{u'}</math> with the corresponding component of <math>~(\rho\mathbf{u} - \rho {\vec{\Omega}}_f \times \vec{x})</math><math>= (\rho\mathbf{u} -\mathbf{\hat{e}}_\varphi \rho R\Omega_0)</math>, that is, we perform the following mappings:

<math> ~\rho u'_R </math>

<math>~\rightarrow~</math>

<math> ~\rho u_R \, , </math>

<math> ~\rho u'_\varphi </math>

<math>~\rightarrow~</math>

<math> ~\rho (u_\varphi - R\Omega_0 ) \, , </math>

<math> ~\rho u'_z </math>

<math>~\rightarrow~</math>

<math> ~\rho u_z \, . </math>

As a result, the first and third of the three "hybrid" momentum-component equations become, respectively,

<math>\mathbf{\hat{e}}_R:</math>

<math> \frac{\partial (\rho u_R)}{\partial t} + \nabla\cdot [(\rho u_R)\mathbf{u'}] </math>

<math>~=~</math>

<math> - \mathbf{\hat{e}}_R \cdot\nabla P - \rho \mathbf{\hat{e}}_R \cdot\nabla \Phi + \frac{\rho u^2_\varphi}{R} \, ; </math>

<math>\mathbf{\hat{k}}:</math>

<math> \frac{\partial (\rho u_z) }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \mathbf{\hat{k}}\cdot\nabla P - \rho \mathbf{\hat{k}}\cdot\nabla \Phi \, . </math>

The second of the three "hybrid" momentum component equations — the one governing conservation of angular momentum — becomes,

<math> \frac{\partial [\rho R (u_\varphi - R\Omega_0 ) ]}{\partial t} + \nabla\cdot \{ [\rho R (u_\varphi - R\Omega_0 ) ] \mathbf{u'} \} </math>

<math>~=~</math>

<math> - \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi - 2\rho R \Omega_0 u_R \, </math>

<math> \Rightarrow ~~~~~ \frac{\partial (\rho R u_\varphi) }{\partial t} + \nabla\cdot [(\rho R u_\varphi) \mathbf{u'} ] - R^2 \Omega_0 \biggl[ \frac{\partial \rho }{\partial t} + \nabla\cdot ( \rho \mathbf{u'} ) \biggr] - ( \rho \mathbf{u'} )\cdot\nabla (R^2 \Omega_0) </math>

<math>~=~</math>

<math> - \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi - 2\rho R \Omega_0 u_R \, . </math>

Referencing the continuity equation, as before, the middle bracketed term on the left-hand side can be set to zero; and the last term on the left-hand side,

<math> ( \rho \mathbf{u'} )\cdot\nabla (R^2\Omega_0) </math>

<math>~=~</math>

<math> 2\rho R \Omega_0 u_R \, , </math>

matches and, hence, exactly cancels the Coriolis term on the right-hand side to give,

<math>\mathbf{\hat{e}}_\varphi:</math>

<math> \frac{\partial (\rho R u_\varphi) }{\partial t} + \nabla\cdot [(\rho R u_\varphi) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \mathbf{\hat{e}}_\varphi \cdot R\nabla P - \rho \mathbf{\hat{e}}_\varphi \cdot R \nabla \Phi \, . </math>

Step 3

Cylindrical Grid: If the numerical simulation is to be conducted on a cylindrical coordinate mesh, the spatial operators on both sides of the component momentum equations should be broken down into their cylindrical-coordinate components. That is, as in the above "ASIDE," the appropriate operator projections on the right-hand-sides of the equations are,

<math>\mathbf{\hat{e}}_R \cdot\nabla \rightarrow \frac{\partial}{\partial R} \, ,</math>      <math>\mathbf{\hat{e}}_\varphi \cdot\nabla \rightarrow \frac{1}{R} \frac{\partial}{\partial \varphi} \, ,</math>      <math>\mathbf{\hat{k}}\cdot\nabla \rightarrow \frac{\partial}{\partial z} \, .</math>

In concert with this, the divergence (advection) term on the left-hand-side should be evaluated by breaking the transport velocity into its three cylindrical components,

<math>~\mathbf{u'} = (u'_R, u'_\varphi, u'_z) \, .</math>

The relevant set of momentum conservation equations is, therefore,

Cylindrical Components of the Inertial-Frame Momentum
advected across a
Rotating, Cylindrical Coordinate Mesh

<math> \frac{\partial (\rho u_R)}{\partial t} + \nabla\cdot [(\rho u_R)\mathbf{u'}] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial R} P - \rho \frac{\partial}{\partial R} \Phi + \frac{\rho u^2_\varphi}{R} </math>

<math> \frac{\partial (\rho R u_\varphi) }{\partial t} + \nabla\cdot [(\rho R u_\varphi) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial \varphi} P - \rho \frac{\partial}{\partial \varphi} \Phi </math>

<math> \frac{\partial (\rho u_z) }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z}\Phi </math>


Cartesian Grid: If, instead, the numerical simulation is to be conducted on a Cartesian coordinate mesh, the spatial operators on both sides of the component momentum equations should be broken down into their Cartesian-coordinate components. In concert with this, the divergence (advection) term on the left-hand-side should be evaluated by breaking the transport velocity into its three Cartesian components,

<math>~\mathbf{u'} = (u'_x, u'_y, u'_z) \, .</math>

Furthermore, recognizing that, when written in Cartesian coordinates, the gradient operator is,

<math> \nabla = \mathbf{\hat{i}} \frac{\partial}{\partial x} + \mathbf{\hat{j}} \frac{\partial}{\partial y} +\mathbf{\hat{k}} \frac{\partial}{\partial z} \, , </math>

and that the unit vectors in Cartesian coordinates can be related to their cylindrical counterparts via the mappings,

<math>\mathbf{\hat{i}} = \mathbf{\hat{e}}_R \cos\varphi - \mathbf{\hat{e}}_\varphi \sin\varphi \, ,</math>      <math>\mathbf{\hat{j}} = \mathbf{\hat{e}}_R \sin\varphi + \mathbf{\hat{e}}_\varphi \cos\varphi \, ,</math>      <math>\mathbf{\hat{k}} = \mathbf{\hat{k}} \, ,</math>

the relevant projections of the gradient operator on the right-hand-sides of the governing equations should take the form,

<math> \mathbf{\hat{e}}_R \cdot\nabla </math>

<math>~\rightarrow~</math>

<math> \biggl[ \cos\varphi \frac{\partial}{\partial x} +\sin\varphi \frac{\partial}{\partial y} \biggr] = \biggl[ \frac{x}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial x} +\frac{y}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial y} \biggr] \, , </math>

<math> R\mathbf{\hat{e}}_\varphi \cdot\nabla </math>

<math>~\rightarrow~</math>

<math> \biggl[ - R\sin\varphi \frac{\partial}{\partial x} +R\cos\varphi \frac{\partial}{\partial y} \biggr] = \biggl[ - y \frac{\partial}{\partial x} +x \frac{\partial}{\partial y} \biggr] \, , </math>

<math> \mathbf{\hat{k}}\cdot\nabla </math>

<math>~\rightarrow~</math>

<math> \frac{\partial}{\partial z} \, . </math>

In summary, then, the relevant set of momentum conservation equations is,

Cylindrical Components of the Inertial-Frame Momentum
advected across a
Rotating, Cartesian Coordinate Mesh

<math> \frac{\partial (\rho u_R)}{\partial t} + \nabla\cdot [(\rho u_R)\mathbf{u'}] </math>

<math>~=~</math>

<math> - \biggl[ \frac{x}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial x} +\frac{y}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial y} \biggr] P - \rho \biggl[ \frac{x}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial x} +\frac{y}{(x^2+y^2)^{1/2}} \frac{\partial}{\partial y} \biggr] \Phi + \frac{\rho u^2_\varphi}{R} </math>

<math> \frac{\partial (\rho R u_\varphi) }{\partial t} + \nabla\cdot [(\rho R u_\varphi) \mathbf{u'} ] </math>

<math>~=~</math>

<math> \biggl[ y \frac{\partial}{\partial x} - x \frac{\partial}{\partial y} \biggr] P + \rho \biggl[ y \frac{\partial}{\partial x} - x \frac{\partial}{\partial y} \biggr] \Phi </math>

<math> \frac{\partial (\rho u_z) }{\partial t} + \nabla\cdot [(\rho u_z) \mathbf{u'} ] </math>

<math>~=~</math>

<math> - \frac{\partial}{\partial z} P - \rho \frac{\partial}{\partial z} \Phi </math>

This is the set of equations that has served as the foundation of the so-called Hybrid simulations reported in Byerly, Adelstein-Lelbach, Tohline, & Marcello (2014).

Related Discussions

Whitworth's (1981) Isothermal Free-Energy Surface

© 2014 - 2021 by Joel E. Tohline
|   H_Book Home   |   YouTube   |
Appendices: | Equations | Variables | References | Ramblings | Images | myphys.lsu | ADS |
Recommended citation:   Tohline, Joel E. (2021), The Structure, Stability, & Dynamics of Self-Gravitating Fluids, a (MediaWiki-based) Vistrails.org publication, https://www.vistrails.org/index.php/User:Tohline/citation