Editing ThreeDimensionalConfigurations/Challenges (section)

==Compressible Structures==

Here we draw heavily on the published work of [http://adsabs.harvard.edu/abs/1996ApJS..105..181K Korycansky &amp; Papaloizou] (1996, ApJS, 105, 181; hereafter KP96) that we have reviewed in a  [[Apps/Korycansky_Papaloizou_1996#Korycansky_and_Papaloizou_.281996.29|separate chapter]], and on the doctoral dissertation of [https://digitalcommons.lsu.edu/gradschool_disstheses/6650/ Saied W. Andalib (1998)].

===Part I===
Returning to the above-mentioned,

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

<math>\frac{\partial \bold{u}}{\partial t} + (\bold{u}\cdot \nabla) \bold{u} = - \frac{1}{\rho} \nabla P - \nabla \Phi_\mathrm{grav} 
- {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x}) - 2{\vec{\Omega}}_f \times \bold{u} \, ,</math> </div>

we next note &#8212; as we have done in our [[PGE/Euler#in_terms_of_the_vorticity:|broader discussion of the Euler equation]] &#8212; that,
<div align="center">
<math> 
(\bold{u} \cdot\nabla)\bold{u} = \frac{1}{2}\nabla(\bold{u} \cdot \bold{u}) - \bold{u} \times(\nabla\times\bold{u})
= \frac{1}{2}\nabla(u^2) + \boldsymbol\zeta \times \bold{u} ,
</math>
</div>
where, as above, <math>\boldsymbol\zeta \equiv \nabla\times \bold{u}</math> is the [https://en.wikipedia.org/wiki/Vorticity vorticity].  Making this substitution, we obtain what is essentially equation (7) of [http://adsabs.harvard.edu/abs/1996ApJS..105..181K KP96], that is, the,

<div align="center">
Euler Equation<br />
written <font color="#770000">'''in terms of the Vorticity'''</font> and<br />
<font color="#770000">'''as viewed from a Rotating Reference Frame'''</font>

<math>\frac{\partial \bold{u}}{\partial t} + (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}}= - \frac{1}{\rho} \nabla P - \nabla \biggl[\Phi + \frac{1}{2}u^2 - \frac{1}{2}|{\vec{\Omega}}_f \times \vec{x}|^2 \biggr]</math> .
</div>
Hence, in steady-state, the Euler equation becomes,
<div align="center">
<math>
\nabla F_B  + \vec{A} = 0 ,
</math>
</div>
where, the scalar "Bernoulli" function,
<div align="center">
<math>
F_B  \equiv  \frac{1}{2}u^2 + H + \Phi - \frac{1}{2}|\Omega\hat{k} \times \vec{x}|^2 ;
</math>
</div>
and,
<div align="center">
<math>
\vec{A} \equiv ({\boldsymbol\zeta}+2{\vec\Omega}_f) \times {\bold{u}} .
</math>
</div>


For later use &hellip;
<table border="1" align="center" cellpadding="10" width="80%"><tr><td align="left">
<ol>
<li><font color="red">Curl of steady-state Euler equation:</font>  
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\nabla F_B + \bold{A}</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\nabla \times \biggl[ \nabla F_B + \bold{A} \biggr]</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\nabla \times  \bold{A} \, .</math>
  </td>
</tr>
</table>
This last step is justified because the [https://en.wikipedia.org/wiki/Vector_calculus_identities#Curl_of_gradient_is_zero curl of any gradient is zero].
</li>
<li>
<font color="red">KP96 only deal with two-dimensional motion in the equatorial plane and, hence, there is no vertical motion:</font><br />
Hence, <math>~\bold{u}</math> lies in the equatorial plane; both <math>~\vec\zeta</math> and <math>~\vec\Omega_f</math> only have z-components; and, <math>~\bold{A}</math> is perpendicular to both <math>~\vec\Omega_f</math> and <math>~\bold{u}</math>.  Also, given that <math>~\bold{A}</math> necessarily lies in the equatorial plane, its curl will only have a z-component, that is,
<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>~\nabla \times  \bold{A} = 0</math>
  </td>
  <td align="center">
&nbsp; &nbsp; &nbsp; <math>~\Leftrightarrow</math>&nbsp; &nbsp; &nbsp; 
  </td>
  <td align="left">
<math>~[\nabla \times  \bold{A}]_z = 0 \, .</math>
  </td>
</tr>
</table>
</li>
<li>
<font color="red">"Dot" <math>~\bold{u}</math> into the steady-state Euler equation:</font>
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\nabla F_B + \bold{A}</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\bold{u} \cdot \biggl[ \nabla F_B + \bold{A} \biggr]</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\bold{u} \cdot \nabla F_B  \, .</math>
  </td>
</tr>
</table>
This last step is justified because <math>~\bold{A}</math> is necessarily always perpendicular to <math>~\bold{u}</math>.
</li>
</ol>
</td></tr></table>


----


We will leave discussion of the Euler equation, for the moment, and instead look at the continuity equation.  As viewed from the rotating frame of reference,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{\partial (\rho \bold{u})}{\partial t} + \nabla\cdot (\rho\bold{u})</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~0 \, .</math>
  </td>
</tr>
</table>

If we are able to write the momentum density (in the rotating frame) in terms of a stream-function, <math>~\Psi</math>, such that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\rho\bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\nabla \times (\boldsymbol{\hat{k}} \Psi) 
=
\boldsymbol{\hat\imath} \biggl[ \frac{\partial \Psi}{\partial y} \biggr]
- \boldsymbol{\hat\jmath} \biggl[  \frac{\partial \Psi}{\partial x}\biggr] \, ,</math>
  </td>
</tr>
</table>
then satisfying the steady-state continuity equation is guaranteed because the [https://en.wikipedia.org/wiki/Vector_calculus_identities#Divergence_of_curl_is_zero divergence of a curl is always zero].  Note, as well, that when written in terms of this stream-function, the z-component of the vorticity will be,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\zeta_z </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial }{\partial x}\biggl[- \frac{1}{\rho} \frac{\partial \Psi}{\partial x}  \biggr] - \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y}  \biggr] \, .
</math>
  </td>
</tr>
</table>

Note that the steady-state continuity equation may be rewritten in the form,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\nabla\cdot \bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- \bold{u} \cdot \nabla [ \ln \rho] \, .
</math>
  </td>
</tr>
</table>



----



It can also be shown that the condition, <math>~[\nabla \times \bold{A}]_z = 0</math> can be rewritten as,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\nabla\cdot \bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- \bold{u} \cdot \nabla [ \ln(2\Omega_f + \zeta_z] \, .
</math>
  </td>
</tr>
</table>

By combining these last two expressions, we appreciate that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
\bold{u} \cdot \nabla \ln \biggl[ \frac{(2\Omega_f + \zeta_z)}{\rho} \biggr] 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~0 \, .</math>
  </td>
</tr>
</table>
This means that, in the steady-state flow whose spatial structure we are seeking, the velocity vector, <math>~\bold{u}</math> (and also the momentum density vector, <math>~\rho \bold{u}</math>) must everywhere be tangent to contours of constant ''vortensity'', <math>~[(2\Omega_f + \zeta_z)/\rho]</math>.

We need a function <math>~g(\Psi) </math> such that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~g(\Psi) </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{\rho} \biggl\{
\zeta_z + 2\Omega_f
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{\rho} \biggl\{
2\Omega_f
-
\frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x}  \biggr] 
- \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y}  \biggr] 
\biggr\} \, .
</math>
  </td>
</tr>
</table>

Let's try, <math>~\Psi = \rho^2</math>, and
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\rho</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\rho_c \biggl\{1 - \biggl[ \frac{y^2}{b^2} + \frac{x^2}{a^2}\biggr]\biggr\} </math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~\frac{\partial^2 \rho}{\partial x^2} = -\frac{\partial}{\partial x}\biggl\{ \frac{2\rho_c x}{a^2}\biggr\} = - \frac{2\rho_c}{a^2}</math>
  </td>
  <td align="center">
&nbsp; &nbsp; and, &nbsp; &nbsp;
  </td>
  <td align="left">
<math>~\frac{\partial^2 \rho}{\partial y^2} = - \frac{\partial}{\partial y}\biggl\{ \frac{2\rho_c y}{b^2} \biggr\} = - \frac{2\rho_c}{b^2} \, .</math>
  </td>
</tr>
</table>

Then,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~g(\Psi) </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{\rho} \biggl\{
2\Omega_f
-
\frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \rho^2}{\partial x}  \biggr] 
- \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \rho^2}{\partial y}  \biggr] 
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{\rho} \biggl\{
2\Omega_f
-
2 \biggl[ \frac{\partial^2 \rho}{\partial x^2}  \biggr] 
- 2 \biggl[\frac{\partial^2 \rho}{\partial y^2}  \biggr] 
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{\rho} \biggl\{
2\Omega_f
+
4\rho_c \biggl[ \frac{1}{a^2} + \frac{1}{b^2}  \biggr] 
\biggr\} \, .
</math>
  </td>
</tr>
</table>

Hence,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~g(\Psi) </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl\{
2\Omega_f
+
4\rho_c \biggl[ \frac{1}{a^2} + \frac{1}{b^2}  \biggr] 
\biggr\} \Psi^{-1 / 2} \, .
</math>
  </td>
</tr>
</table>

Next, given that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{dF_B}{d\Psi}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- g(\Psi) \, ,</math>
  </td>
</tr>
</table>
we conclude that, to within an additive constant,

<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>~ F_B(\Psi)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- \int g(\Psi) d\Psi
= - g_0 \int \Psi^{-1 / 2} d\Psi
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- 2g_0 \Psi^{1 / 2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- 2g_0 \rho \, , 
</math>
  </td>
</tr>

</table>
where,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~g_0</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\biggl\{ 2\Omega_f + 4\rho_c \biggl[ \frac{1}{a^2} + \frac{1}{b^2}  \biggr] \biggr\} \, .
</math>
  </td>
</tr>
</table>

Here's what to do:

Given <math>~g(\Psi)</math>, write out the functional forms of <math>~\bold{A}</math> and <math>~F_B(\Psi)</math>.  Then see if <math>~\nabla F_B = - \bold{A}</math>.

===Part II===

Consider a steady-state configuration that is the compressible analog of a Riemann S-type ellipsoid; even better, give the configuration a "peanut" shape rather than the shape of an ellipsoid.  As viewed from a frame that is spinning with the configuration's overall angular velocity, <math>~\vec\Omega_f = \boldsymbol{\hat{k}} \Omega_f</math>, generally we expect the configuration's internal (and surface) flow to be represented by a set of nested stream-lines and at every <math>~(x, y)</math> location the fluid's velocity (and its momentum-density vector) will be tangent to the stream-line that runs through that point.  It is customary to represent the stream-function by a scalar quantity, <math>~\Psi(x, y)</math>, appreciating that each stream-line will be defined by a curve, <math>~\Psi = \mathrm{constant}</math>; also, the local spatial gradient of <math>~\Psi(x,y)</math> will be perpendicular to the local stream-line, hence it will be perpendicular to the local velocity vector as well.  If we specifically introduce the stream-function via the relation,

<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>~\rho\bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\nabla \times (\boldsymbol{\hat{k}} \Psi) 
=
\boldsymbol{\hat\imath} \biggl[ \frac{\partial \Psi}{\partial y} \biggr]
- \boldsymbol{\hat\jmath} \biggl[  \frac{\partial \Psi}{\partial x}\biggr] \, ,</math>
  </td>
</tr>
</table>

it will display all of the just-described attributes and we are also guaranteed that the steady-state continuity equation will be satisfied everywhere, because the divergence of a curl is always zero.

We also have demonstrated that the vector, <math>~\bold{A}</math>, has the right properties if,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
\bold{u} \cdot \nabla \ln \biggl[ \frac{(2\Omega_f + \zeta_z)}{\rho} \biggr] 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~0 \, .</math>
  </td>
</tr>
</table>

This means that, at every location in the plane of the fluid flow, the gradient of the ''vortensity'' also must be perpendicular to the velocity vector.  This constraint can be immediately satisfied if we simply demand that the vortensity be a function of the stream-function, <math>~\Psi</math>, that is, we need,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{(2\Omega_f + \zeta_z)}{\rho}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~g(\Psi) \, .</math>
  </td>
</tr>
</table>
In other words, the scalar vortensity is constant along each stream-line.  And, once we have determined the mathematical expression for this function, we will know that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\bold{A}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~[ \boldsymbol{\hat{k}} g(\Psi) ] \times \rho\bold{u} \, ;</math>
  </td>
</tr>
</table>
Furthermore, we should be able to determine the mathematical expression for <math>~F_B(x,y)</math> because,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{dF_B}{d\Psi}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- g(\Psi) \, .</math>
  </td>
</tr>
</table>
As a check, we should find that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\nabla F_B + \bold{A}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~0 \, .</math>
  </td>
</tr>
</table>

===Part III===

Here we closely follow Chapter 4 of [https://digitalcommons.lsu.edu/gradschool_disstheses/6650/ Saied W. Andalib (1998)].  

From &sect;4.1 (p. 80):  "<font color="orange">Euler's equation, the equation of continuity, the Poisson equation and the equation of state &hellip; govern the dynamics and evolution of these equilibrium configurations.</font>"

====Equation of Continuity====
In steady state, <math>~\partial (\rho\bold{u})/\partial t = 0</math>. Hence the rotating-frame-based continuity equation becomes,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\nabla \cdot (\rho\bold{u})</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~0 \, .</math>
  </td>
</tr>
</table>
If we insist that the momentum-density vector be expressible in terms of the curl of a vector &#8212; for example,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\rho\bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\nabla \times \boldsymbol{\mathfrak{J}} \, ,</math>
  </td>
</tr>
<tr><td align="center" colspan="3">[https://digitalcommons.lsu.edu/gradschool_disstheses/6650/ Saied W. Andalib (1998)], &sect;4.1, p. 80, Eq. (4.1)
</td></tr>
</table>

then satisfying this steady-state continuity equation is guaranteed because the [https://en.wikipedia.org/wiki/Vector_calculus_identities#Divergence_of_curl_is_zero divergence of a curl is always zero].   "<font color="orange">The task</font> of satisfying the steady-state equation of continuity <font color="orange">then shifts to identifying an appropriate expression for the vector potential, <math>~\boldsymbol{\mathfrak{J}} \, .</math></font>"  In the most general case, in terms of this vector potential the three Cartesian components of the momentum-density vector are,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\rho u_x</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial \mathfrak{J}_z}{\partial y}
-
\frac{\partial \mathfrak{J}_y}{\partial z} \, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\rho u_y</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial \mathfrak{J}_x}{\partial z}
-
\frac{\partial \mathfrak{J}_z}{\partial x} \, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\rho u_z</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial \mathfrak{J}_y}{\partial x}
-
\frac{\partial \mathfrak{J}_x}{\partial y} \, .
</math>
  </td>
</tr>
</table>
Here, <font color="red">we will follow Andalib's lead and only look for fluid flows with no vertical motions</font>.  That is to say, we will set <math>~\rho u_z = 0</math>, in which case this last expression establishes the constraint,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{\partial \mathfrak{J}_y}{\partial x}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial \mathfrak{J}_x}{\partial y} \, .
</math>
  </td>
</tr>
</table>
"<font color="orange">A general solution to this equation can be found only if there exists a scalar function <math>~\Gamma(x, y, z)</math> such that &hellip;</font>"

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\mathfrak{J}_y = \frac{\partial \Gamma}{\partial y}</math>
  </td>
  <td align="center">
&nbsp; &nbsp; &nbsp; and, &nbsp; &nbsp; &nbsp; 
  </td>
  <td align="left">
<math>~\mathfrak{J}_x = \frac{\partial \Gamma}{\partial x} \, ;</math>
  </td>
</tr>
</table>
note that this adopted functional behavior works because the constraint becomes,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{\partial^2 \Gamma}{\partial x \partial y}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{\partial^2 \Gamma}{\partial y \partial x} \, .</math>
  </td>
</tr>
</table>

Hence, the expressions for the x- and y-components of the momentum-density vector may be rewritten, respectively, as,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\rho u_x</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial \mathfrak{J}_z}{\partial y}
-
\frac{\partial}{\partial z} \biggl[ \frac{\partial \Gamma}{\partial y} \biggr]
=
+ \frac{\partial}{\partial y}\biggl[ \mathfrak{J}_z - \frac{\partial\Gamma}{\partial z} \biggr]
\, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\rho u_y</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial }{\partial z}\biggl[ \frac{\partial \Gamma}{\partial x} \biggr]
-
\frac{\partial \mathfrak{J}_z}{\partial x} 
=
- \frac{\partial}{\partial x}\biggl[ \mathfrak{J}_z - \frac{\partial\Gamma}{\partial z} \biggr]
\, .
</math>
  </td>
</tr>
</table>

If we again <font color="red">follow Andalib's lead and only look for models in which the x-y-plane flow is independent of the vertical coordinate, z</font>, then, <font color="orange"><math>~\mathfrak{J}_z</math> and <math>~\partial\Gamma/\partial z</math> must be functions of x and y only.  Therefore, <math>~\mathfrak{J}_z</math> is independent of z and <math>~\Gamma</math> is at most linear in z.</font>  Now, rather than focusing on the determination of <math>~\Gamma(x,y)</math>, we can just as well define the scalar function, 
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\Psi(x,y)</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\mathfrak{J}_z - \frac{\partial\Gamma}{\partial z} \, ,
</math>
  </td>
</tr>
</table>
in which case "<font color="orange">&hellip; the components of the momentum density may be written as:</font>"
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\rho \bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\boldsymbol{\hat\imath} \frac{\partial \Psi}{\partial y} - \boldsymbol{\hat\jmath} \frac{\partial \Psi}{\partial x} \, .
</math>
  </td>
</tr>
</table>
It is straightforward to demonstrate that this expression for the momentum-density vector does satisfy the steady-state continuity equation.  "<font color="orange">The function <math>~\Psi(x, y)</math> will serve a similar role as the velocity potential for incompressible fluids.</font>"

====Related Useful Expressions====

Given that, by our design, the fluid motion will be confined to the x-y-plane, the fluid vorticity will have only a z-component; that is, 
<div align="center">
<math>~\vec\zeta = \nabla \times \bold{u} = \boldsymbol{\hat{k}} \zeta_z</math>,  
</div>
where,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\zeta_z </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \, .
</math>
  </td>
</tr>
</table>

And when it is written in terms of <math>~\Psi(x,y)</math>, this z-component of the vorticity will be obtained from the expression,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\zeta_z </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\partial }{\partial x}\biggl[- \frac{1}{\rho} \frac{\partial \Psi}{\partial x}  \biggr] - \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y}  \biggr] \, .
</math>
  </td>
</tr>
</table>
<span id="LHS">This is useful to know because</span>, in the Euler equation (see immediately below) we will encounter a term that involves the cross product of the vector, <math>~(\vec\zeta + 2\vec\Omega) </math>, with the rotating-frame-based velocity vector.  Appreciating as well that the vector, <math>~\vec\Omega = \boldsymbol{\hat{k}} \Omega_f</math>, only has a nonzero z-component, we recognize that this term may be written as,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~(\vec\zeta + 2\vec\Omega) \times \bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\boldsymbol{\hat{k}} \biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr] \times \biggl[ \boldsymbol{\hat\imath} \frac{\partial \Psi}{\partial y} - \boldsymbol{\hat\jmath} \frac{\partial \Psi}{\partial x} \biggr]</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr] \biggl[ \boldsymbol{\hat\jmath} \frac{\partial \Psi}{\partial y} + \boldsymbol{\hat\imath} \frac{\partial \Psi}{\partial x} \biggr]</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr] \nabla\Psi \, .</math>
  </td>
</tr>
<tr><td align="center" colspan="3">[https://digitalcommons.lsu.edu/gradschool_disstheses/6650/ Saied W. Andalib (1998)], &sect;4.2, p. 83, Eq. (4.13)</td></tr>
</table>

====Euler Equation====
We begin with the,

<div align="center">
Euler Equation<br />
written <font color="#770000">'''in terms of the Vorticity'''</font> and<br />
<font color="#770000">'''as viewed from a Rotating Reference Frame'''</font>

<math>\frac{\partial \bold{u}}{\partial t} + (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}}= - \frac{1}{\rho} \nabla P - \nabla \biggl[\Phi + \frac{1}{2}u^2 - \frac{1}{2}|{\vec{\Omega}}_f \times \vec{x}|^2 \biggr]</math> .
</div>

Next, we rewrite this expression to incorporate the following three realizations:
<ul>
<li>For a barotropic fluid, the term involving the pressure gradient can be replaced with a term involving the enthalpy via the relation, <math>~\nabla H = \nabla P/\rho</math>.</li>
<li>The expression for the centrifugal potential can be rewritten as, <math>~\tfrac{1}{2}|\vec\Omega_f \times \vec{x}|^2 = \tfrac{1}{2}\Omega_f^2 (x^2 + y^2)</math>.</li>
<li>In steady state, <math>~\partial \bold{u}/\partial t = 0</math>.</li>
</ul>
This means that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~ (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>- \nabla \biggl[H + \Phi_\mathrm{grav} + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] \, .</math>
  </td>
</tr>
</table>

If the term on the left-hand-side of this equation can be expressed in terms of the gradient of a scalar function, then it can be readily grouped with all the other terms on the right-hand-side, which already are in the gradient form.

<table border="1" align="center" width="80%" cellpadding="10"><tr><td align="left">
<div align="center">
'''Striving for Gradient Form'''
</div>

As we have already demonstrated [[#LHS|above]], the term on the left-hand-side of the Euler equation can be rewritten as,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~(\vec\zeta + 2\vec\Omega) \times \bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr] \nabla\Psi \, .</math>
  </td>
</tr>
</table>

If the term inside the square brackets on the right-hand-side were a constant &#8212; that is, independent of position &#8212; then it could immediately be moved inside the gradient operator and we will have accomplished our objective.  But, while <math>~\Omega_f</math> ''is'' constant "<font color="orange">&hellip; generally the vorticity and density are both functions of <math>~x</math> and <math>~y</math>.</font>"  As Andalib has explained, "<font color="orange">The expression &hellip; can be cast in the form of a gradient only if</font>
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr]</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~g(\Psi) \, ,</math>
  </td>
</tr>
</table>
<font color="orange">where <math>~g</math> is an arbitrary function.</font>"  Specifically in this case, the term on the left-hand-side of the Euler equation may be written as,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~(\vec\zeta + 2\vec\Omega) \times \bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\nabla F_B(\Psi) </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl[ \frac{dF_B(\Psi)}{d\Psi} \biggr] \nabla \Psi \, .</math>
  </td>
</tr>
</table>
That is, we accomplish our objective by recognizing that the sought-after function, <math>~F_B(\Psi)</math>, is obtained from <math>~g(\Psi)</math> via the relation,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{dF_B(\Psi)}{d\Psi} </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~g(\Psi) \, .</math>
  </td>
</tr>
</table>

----

For example, try &hellip;
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~F_B(\Psi)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~C_0 \Psi + \frac{1}{2} C_1 \Psi^2 </math>
  </td>
</tr>
<tr><td align="center" colspan="3">[https://digitalcommons.lsu.edu/gradschool_disstheses/6650/ Saied W. Andalib (1998)], &sect;4.3, p. 85, Eq. (4.22)</td></tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ g(\Psi)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~C_0 + C_1 \Psi \, .</math>
  </td>
</tr>
</table>

This means that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr]</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~~C_0 + C_1 \Psi</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ 
2\Omega_f - C_0 \rho
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~~C_1 \rho \Psi - \zeta_z </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~~C_1 \rho \Psi + \frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x}  \biggr] + \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y}  \biggr] \, .</math>
  </td>
</tr>
<tr><td align="center" colspan="3">[https://digitalcommons.lsu.edu/gradschool_disstheses/6650/ Saied W. Andalib (1998)], &sect;4.3, p. 85, Eq. (4.24)</td></tr>
</table>

</td></tr></table>


<span id="AndalibBernoulli">Having transformed the left-hand-side term</span> into the gradient of the scalar function, <math>~F_B(\Psi)</math>, the Euler equation can now be written as,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\nabla \biggl[H + \Phi_\mathrm{grav} + F_B(\Psi) + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ 0</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~
H + \Phi_\mathrm{grav} + F_B(\Psi) + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2)  
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>C_B \, , </math>
  </td>
</tr>
</table>
where we will refer to <math>~C_B</math> as the Bernoulli constant.

====Strategy====

<font color="red">'''STEP 0:'''</font>&nbsp; Choose the pair of model-sequence parameters, <math>~(C_0, C_1)</math>, that are associated with the function, <math>~F_B(\Psi)</math>.  Hold these fixed during iterations.

<font color="red">'''STEP 1:'''</font>&nbsp; Guess a density distribution, <math>~\rho(x,y)</math>.  For example, pick the equatorial-plane (uniform) density distribution of a Riemann S-type ellipsoid with an equatorial-plane axis-ratio, <math>~b/a</math> and meridional-plane axis-ratio, <math>~c/a</math>; use the same <math>~b/a</math> ratio to define two points on the configuration's surface throughout the iteration cycle.

<font color="red">'''STEP 2:'''</font>&nbsp; Given <math>~\rho(x,y)</math>, solve the Poisson equation to obtain, <math>~\Phi_\mathrm{grav}(x,y)</math>.  In the first iteration, this should exactly match the <math>~A_1, A_2, A_3</math> values associated with the chosen Riemann S-type ellipsoid.

<font color="red">'''STEP 3:'''</font>&nbsp; Guess a value of <math>~\Omega_f</math> &#8212; perhaps the spin-frequency associated with your "initial guess" Riemann ellipsoid &#8212; then solve the following two-dimensional, elliptic PDE to obtain <math>~\Psi(x,y)</math> &hellip;
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~ 
2\Omega_f - C_0 \rho
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~~C_1 \rho \Psi + \frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x}  \biggr] + \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y}  \biggr] \, .</math>
  </td>
</tr>
</table>
<table border="1" cellpadding="10" width="70%" align="center"><tr><td align="left">
<div align="center">'''Boundary Condition'''</div>

Moving along various rays from the center of the configuration, outward, the surface is determined by the location along each ray where <math>~H(x,y)</math> goes to zero for the first time.  We set <math>~\Psi = 0</math> at these various surface locations.  At each of these locations, the velocity vector must be tangent to the surface.  This requirement also, then, sets the value of <math>~\partial \Psi/\partial y</math> and <math>~\partial \Psi/\partial x</math> at each location.
</td></tr></table>

<font color="red">'''STEP 4:'''</font>&nbsp; Determine (rotating-frame) velocity from knowledge of <math>~\Psi(x,y)</math> and <math>~\rho(x,y)</math>.
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{\rho}\biggl\{ \boldsymbol{\hat\imath} \biggl[ \frac{\partial \Psi}{\partial y} \biggr]
- \boldsymbol{\hat\jmath} \biggl[  \frac{\partial \Psi}{\partial x}\biggr] \biggr\} </math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~ u^2 = \bold{u} \cdot \bold{u}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{\rho^2}\biggl\{ \biggl[ \frac{\partial \Psi}{\partial y} \biggr]^2
+ \biggl[  \frac{\partial \Psi}{\partial x}\biggr]^2 \biggr\} \, .</math>
  </td>
</tr>
</table>

<font color="red">'''STEP 5:'''</font>&nbsp; Using the "scalar Euler equation,"

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
H + \Phi_\mathrm{grav} + C_0 \Psi + \frac{1}{2} C_1 \Psi^2 + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2)  
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>C_B \, , </math>
  </td>
</tr>
<tr><td align="center" colspan="3">[https://digitalcommons.lsu.edu/gradschool_disstheses/6650/ Saied W. Andalib (1998)], &sect;4.3, p. 85, Eq. (4.23)</td></tr>
</table>
<ul>
<li>Set <math>~H = 0</math> at two different points on the surface of the configuration &#8212; usually at <math>~(x,y) = (a,0)</math> and <math>~(x,y) = (0,b)</math> &#8212; to determine values of the two constants, <math>~\Omega_f^2</math> and <math>~C_B</math>.</li>
<li>At all points inside the configuration, determine <math>~H(x,y)</math>.</li>
</ul>
<font color="red">'''STEP 6:'''</font>&nbsp; Use the barotropic equation of state to determine the "new" mass-density distribution from the knowledge of the enthalpy, <math>~H(x,y)</math>.