Joel2: Created page with "FORCETOC =Compressible Analogs of Riemann S-Type Ellipsoids= Here we attempt to develop a self-consistent-field type, iterative technique that will permit the construction of steady-state structures that are compressible analogs of Riemann S-Type (incompressible) ellipsoids. We will build upon the recent work of [http://adsabs.harvard.edu/abs/2006ApJ...639..549O Ou (2006)]. ==Standard Stea..."

2023-12-11T22:40:25Z

Created page with "__FORCETOC__   =Compressible Analogs of Riemann S-Type Ellipsoids= Here we attempt to develop a self-consistent-field type, iterative technique that will permit the construction of steady-state structures that are compressible analogs of Riemann S-Type (incompressible) ellipsoids. We will build upon the recent work of [http://adsabs.harvard.edu/abs/2006ApJ...639..549O Ou (2006)]. ==Standard Stea..."

New page

__FORCETOC__ 


=Compressible Analogs of Riemann S-Type Ellipsoids=
Here we attempt to develop a self-consistent-field type, iterative technique that will permit the construction of steady-state structures that are compressible analogs of Riemann S-Type (incompressible) ellipsoids. We will build upon the recent work of [http://adsabs.harvard.edu/abs/2006ApJ...639..549O Ou (2006)].

==Standard Steady-State Governing Relations==
As viewed from a rotating frame of reference and written in Eulerian form, the steady-state version of the three-dimensional principal governing equations are:
<div align="center">
<math>
\nabla\cdot(\rho \vec{v}) = 0
</math>

<math>
(\vec{v}\cdot \nabla)\vec{v} = -\nabla \biggl[H + \Phi -\frac{1}{2}\omega^2 R^2 \biggr] -2\vec{\omega}\times\vec{v}
</math>

<math>
\nabla^2 \Phi = 4\pi G \rho
</math>
</div>

==Proposed Solution Strategy==
===Preamble:===
Specify the three "polar" boundary locations, <math>a, b,</math> and <math>c</math>; specify the <i>direction</i> but not the magnitude of the rotating frame's angular velocity, for example, <math>(\vec{\omega}/\omega) = \hat{k}</math>; pin the central density to the value <math>\rho_c = 1</math>. Define the following dimensionless density, velocity vector, angular velocity, enthalpy, gravitational potential, and position vector:
<div align="center">
<math>
\rho^* \equiv \frac{\rho}{\rho_c} ; ~~~~~{\vec{v}}^* \equiv \frac{\vec{v}}{[a^2G\rho_c]^{1/2}} ; ~~~~~\omega^* \equiv \frac{\omega}{[G\rho_c]^{1/2}} ;
</math>

<math>
H^* \equiv \frac{H}{[a^2G\rho_c]} ; ~~~~~\Phi^* \equiv \frac{\Phi}{[a^2G\rho_c]} ; ~~~~~{\vec{x}}^* \equiv \frac{\vec{x}}{a} .
</math>
</div>
From here, on, spatial operators are assumed to be in terms of the dimensionless coordinates.

===Step #1:===
Guess a 3D density distribution — such as a uniform-density ellipsoid, or one of the converged models from Ou (2006) — that conforms to a selected set of <i>positional</i> boundary conditions, that is, where the density goes to zero along the three principal axes at <math>x=a</math>, <math>y = b</math>, and <math>z = c</math>. Solve the Poisson equation in order to derive values for <math>\Phi</math> everywhere inside and on the surface of the 3D configuration:
<div align="center">
<math>
\nabla^2 \Phi^* = 4\pi \rho^* .
</math>
</div>

===Step #2:===
Use the continuity equation and the curl of the Euler equation to numerically derive the <i>structure</i> but not the overall magnitude of the velocity flow-field throughout the 3D configuration. Take advantage of the fact that the direction, <math>(\vec{\omega}/\omega)</math>, has been specified; and assume that the 3D density distribution is known. Here are the relevant equations:
<div align="center">
<math>
\nabla\cdot(\rho^* {\vec{v}}^*) = 0 ;
</math>

<math>
\nabla\times \biggl[({\vec{v}}^*\cdot \nabla){\vec{v}}^* +2 {\vec{\omega}}^* \times {\vec{v}}^* \biggr] = 0 .
</math>
</div>

The first of these is a scalar equation; the second is a vector equation and it will presumably provide two useful scalar equations (perhaps constraining the two components of <math>{\vec{v}}^*</math> that are perpendicular to <math>\hat{k}</math> ?). Since the left-hand-side of the second equation is obviously nonlinear in the velocity, we may have to linearize this set of equations and look for small "corrections" <math>\delta\vec{v}</math> to an initial "guess" for the velocity field, such as the flow field in Riemann S-type ellipsoids, which is also the flow-field adopted by Ou (2006).

===Step #3:===
Take the divergence of the Euler equation and use it to solve for <math>H</math> throughout the configuration, given the structure of the flow-field obtained in Step #2. Boundary conditions at the three "poles" of the configuration may suffice to uniquely determine <math>\omega</math>, the overall normalization factor for the flow-field <math>\vec\zeta</math> — hopefully this is analogous to solving for the vorticity parameter <math>\lambda</math> in Ou (2006) — and the Bernoulli constant (or something equivalent). The relevant "Poisson"-like equation is:

<div align="center">
<math>
\nabla^2 \biggl[H^* + \Phi^* -\frac{1}{2}(\omega^*)^2 \biggl(\frac{R}{a}\biggr)^2 \biggr] = - \nabla\cdot [({\vec{v}}^*\cdot \nabla){\vec{v}}^* + 2 {\vec{\omega}}^*\times {\vec{v}}^* ] .
</math>
</div>

==Example of Riemann S-Type Ellipsoids==
===Preamble===

First, set <math>{\vec{\omega}} = \hat{k}\omega</math> and <math>v_z = 0</math>, and write out the Cartesian components of the vector,
<div align="center">
<math>
\vec{A} \equiv ({\vec{v}}\cdot \nabla){\vec{v}} +2 {\vec{\omega}} \times {\vec{v}} .
</math>
</div>
The components are:
<div align="left">
<math>
~~~~~\hat{i}:~~~~~A_x = v_x \frac{\partial v_x}{\partial x} + v_y \frac{\partial v_x}{\partial y} -2\omega v_y ;
</math><br />
<math>
~~~~~\hat{j}:~~~~~A_y = v_x \frac{\partial v_y}{\partial x} + v_y \frac{\partial v_y}{\partial y} +2\omega v_x ;
</math><br />
<math>
~~~~~\hat{k}:~~~~~A_z = 0 .
</math>
</div>
The curl of <math>\vec{A}</math> (needed in Step #2, above) produces a vector with the following three Cartesian components:
<div align="left">
<math>
~~~~~\hat{i}:~~~~~[\nabla\times\vec{A}]_x = \frac{\partial}{\partial y} \biggl[0 \biggr] - \frac{\partial}{\partial z} \biggl[ v_x \frac{\partial v_y}{\partial x} + v_y \frac{\partial v_y}{\partial y} +2\omega v_x \biggr] ;
</math><br />
<math>
~~~~~\hat{j}:~~~~~[\nabla\times\vec{A}]_y = \frac{\partial}{\partial z} \biggl[ v_x \frac{\partial v_x}{\partial x} + v_y \frac{\partial v_x}{\partial y} -2\omega v_y \biggr] - \frac{\partial}{\partial x} \biggl[0 \biggr] ;
</math><br />
<math>
~~~~~\hat{k}:~~~~~[\nabla\times\vec{A}]_z = \frac{\partial}{\partial x} \biggl[ v_x \frac{\partial v_y}{\partial x} + v_y \frac{\partial v_y}{\partial y} +2\omega v_x \biggr] - \frac{\partial}{\partial y} \biggl[ v_x \frac{\partial v_x}{\partial x} + v_y \frac{\partial v_x}{\partial y} -2\omega v_y \biggr] .
</math>
</div>
Hence, demanding (as in Step #2) that <math>\nabla\times\vec{A} = 0</math> means that each of these components independently must be zero. This, in turn, implies:
<div align="left">
<math>
~~~~~\hat{i}:~~~~~ v_x \frac{\partial v_y}{\partial x} + v_y \frac{\partial v_y}{\partial y} +2\omega v_x = C_{z1}(x,y);
</math><br />
<math>
~~~~~\hat{j}:~~~~~ v_x \frac{\partial v_x}{\partial x} + v_y \frac{\partial v_x}{\partial y} -2\omega v_y = C_{z2}(x,y) ;
</math><br />
<math>
~~~~~\hat{k}:~~~~~ \frac{\partial}{\partial x} \biggl[ C_{z1}(x,y) \biggr] = \frac{\partial}{\partial y} \biggl[C_{z2}(x,y) \biggr] ,
</math>
</div>
where the integration "constants" <math>C_{z1}</math> and <math>C_{z2}</math> may be functions of <math>x</math> and/or <math>y</math> but they must be independent of <math>z</math>.

Generically, the continuity equation demands,
<div align="center">
<math>
0 = \frac{\partial}{\partial x}\biggl[ \rho v_x \biggr] + \frac{\partial}{\partial y}\biggl[ \rho v_y \biggr] + \frac{\partial}{\partial z}\biggl[ \rho v_z \biggr] .
</math>
</div>

The divergence of <math>\vec{A}</math> (providing the right-hand-side of the Poisson-like equation in Step #3, above) generates:
<div align="center">
<math>
\nabla\cdot\vec{A} = \frac{\partial}{\partial x} \biggl[ v_x \frac{\partial v_x}{\partial x} + v_y \frac{\partial v_x}{\partial y} -2\omega v_y \biggr] + \frac{\partial}{\partial y} \biggl[ v_x \frac{\partial v_y}{\partial x} + v_y \frac{\partial v_y}{\partial y} +2\omega v_x \biggr] + \frac{\partial}{\partial z} \biggl[ 0 \biggr]
</math><br />
<math>
= \frac{\partial}{\partial x} \biggl[ C_{z2}(x,y) \biggr] + \frac{\partial}{\partial y} \biggl[C_{z1}(x,y) \biggr] .
</math>
</div>

===Riemann Flow-Field===
In Riemann S-Type ellipsoids, the adopted planar flow-field as viewed from the rotating reference frame is,
<div align="center">
<math>
\vec{v} = \hat{i} \biggl( \frac{\lambda a y}{b} \biggr) + \hat{j} \biggl( - \frac{\lambda b x}{a} \biggr) .
</math>
</div>
Hence,
<div align="center">
<math>
C_{z1}(x,y) = \biggl( \frac{\lambda a y}{b} \biggr) \frac{\partial}{\partial x}\biggl( - \frac{\lambda b x}{a} \biggr) + \biggl( - \frac{\lambda b x}{a} \biggr) \frac{\partial}{\partial y}\biggl( - \frac{\lambda b x}{a} \biggr) +2\omega \biggl( \frac{\lambda a y}{b} \biggr) = \biggl[2\omega\biggl( \frac{\lambda a }{b} \biggr)- \lambda^2 \biggr]y,
</math><br />
<math>
C_{z2}(x,y) = \biggl( \frac{\lambda a y}{b} \biggr) \frac{\partial}{\partial x}\biggl( \frac{\lambda a y}{b} \biggr) + \biggl( - \frac{\lambda b x}{a} \biggr) \frac{\partial}{\partial y}\biggl( \frac{\lambda a y}{b} \biggr) -2\omega \biggl( - \frac{\lambda b x}{a} \biggr) = \biggl[2\omega\biggl( \frac{\lambda b }{a} \biggr) - \lambda^2 \biggr]x .
</math>
</div>
Because <math>C_{z1}</math> is independent of <math>x</math> and <math>C_{z2}</math> is independent of <math>y</math>, we see that <math>[\nabla\times\vec{A}]_z = 0</math>, trivially. With this specified velocity flow-field and the appreciation that Riemann S-type ellipsoids also have uniform density, the continuity equation is also trivially satisfied; specifically,
<div align="center">
<math>
\frac{\partial}{\partial x}\biggl[ \rho v_x \biggr] + \frac{\partial}{\partial y}\biggl[ \rho v_y \biggr] + \frac{\partial}{\partial z}\biggl[ \rho v_z \biggr] = \rho \biggl[ \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \biggr] = 0 .
</math>
</div>
However, the right-hand-side of our Poisson-like equation is not zero; rather, it is,
<div align="center">
<math>
\nabla\cdot\vec{A} = \biggl[2\omega\biggl( \frac{\lambda b }{a} \biggr) - \lambda^2 \biggr] + \biggl[2\omega\biggl( \frac{\lambda a }{b} \biggr)- \lambda^2 \biggr] = 2 \biggl[\omega\lambda \biggl( \frac{b}{a} + \frac{a}{b} \biggr) - \lambda^2 \biggr].
</math>
</div>

===Summary===
What can we learn from the Riemann S-Type ellipsoids? Well, let's assume that the <i>structure</i> of our equilibrium flow-field will be more complicated than simply "flow along elliptical paths", but let's continue to assume that a solution can be found in which the flow remains <i>planar</i>, that is, <math>v_z=0</math> everywhere. Also, let's continue to align the rotation axis of the frame with the <math>z</math>-axis of the configuration. The three steps in our proposed solution strategy become:

====<font color="darkblue">Step #1 (simplified):</font>====
Guess a 3D density distribution where the density goes to zero along the three principal axes at <math>x=a</math>, <math>y = b</math>, and <math>z = c</math>. Solve the Poisson equation in order to derive values for <math>\Phi</math> everywhere inside and on the surface of the 3D configuration:
<div align="center">
<math>
\nabla^2 \Phi = 4\pi G \rho .
</math>
</div>

====<font color="darkblue">Step #2 (simplified):</font>====
Use the continuity equation and the curl of the Euler equation to numerically derive the <i>structure</i> but not the overall magnitude of the velocity flow-field throughout the 3D configuration, assuming <math>\vec{\omega}=\hat{k}\omega</math> and the 3D density distribution is known. Here are the relevant equations:
<div align="center">
<math>
\frac{\partial}{\partial x}\biggl[ \rho v_x \biggr] + \frac{\partial}{\partial y}\biggl[ \rho v_y \biggr] = 0 .
</math>
</div>
<div align="left">
<math>
~~~~~\hat{i}:~~~~~ v_x \frac{\partial v_y}{\partial x} + v_y \frac{\partial v_y}{\partial y} +2\omega v_x = C_{z1}(x,y);
</math><br />
<math>
~~~~~\hat{j}:~~~~~ v_x \frac{\partial v_x}{\partial x} + v_y \frac{\partial v_x}{\partial y} -2\omega v_y = C_{z2}(x,y) ;
</math><br />
<math>
~~~~~\hat{k}:~~~~~ \frac{\partial}{\partial x} \biggl[ C_{z1}(x,y) \biggr] = \frac{\partial}{\partial y} \biggl[C_{z2}(x,y) \biggr] ,
</math>
</div>

Now, for the adopted flow-field of the Riemann ellipsoids, <math>C_{z1}(x,y)</math> is only a function of <math>y</math> and <math>C_{z2}(x,y)</math> is only a function of <math>x</math> — that is, <math>C_{z1}(x,y) \rightarrow C_{z1}(y)</math> and <math>C_{z2}(x,y) \rightarrow C_{z2}(x)</math>. Hence, the third (<math>\hat{k}</math>) condition is automatically satisfied. I don't know whether or not we will find that our more general velocity fields exhibit the same nice character.

====<font color="darkblue">Step #3 (simplified):</font>====
Take the divergence of the Euler equation and use it to solve for <math>H</math> throughout the configuration, given the structure of the flow-field obtained in Step #2. The relevant "Poisson"-like equation is:

<div align="center">
<math>
\nabla^2 \biggl[H + \Phi -\frac{1}{2}\omega^2 (x^2 + y^2) \biggr] = - \frac{\partial}{\partial x} \biggl[ C_{z2}(x,y) \biggr] - \frac{\partial}{\partial y} \biggl[C_{z1}(x,y) \biggr] .
</math>
</div>
In the Riemann ellipsoids, the RHS of this expression is a constant. More importantly, in the Riemann case, it is possible to bring the constants from Step #2 inside the spatial operator on the LHS and establish the following <i>algebraic</i> condition:
<div align="center">
<math>
H + \Phi -\frac{1}{2}\omega^2 (x^2 + y^2) + f(x,y) = C_\mathrm{Bernoulli} ,
</math>
</div>
where the function <math>f(x,y)</math> contains only quadratic terms in <math>x</math> and <math>y</math>. Specifically [see Eq. (6) of Ou (2006)],
<div align="center">
<math>
f(x,y) = x^2 \biggl[\omega\lambda \frac{b}{a}-\frac{\lambda^2}{2} \biggr] + y^2 \biggl[ \omega\lambda \frac{a}{b}-\frac{\lambda^2}{2} \biggr] .
</math>
</div>
This makes sense because in the Riemann case the equipotential surfaces were perfect ellipsoids, just like the shape of the centrifugal potential term. So the remaining "source" terms could adopt the same "shape". But in our more general case, surfaces of constant <math>\Phi</math> will have more complicated shapes, so it will be necessary for the function <math>f(x,y)</math> (or equivalent) to incorporate higher order terms in <math>x</math> and <math>y</math>. Is it possible that equilibrium configurations can be found by including additional terms just in <math>x^4</math> and <math>y^4</math>?

==A Related Two-dimensional Treatment==
===Preamble===
In a related work, [http://adsabs.harvard.edu/abs/1996ApJS..105..181K Korycansky & Papaloizou] (1996, ApJS, 105, 181; hereafter KP96) developed a method to find nontrivial, nonaxisymmetric steady-state flows in a <i>two-dimensional</i> setting. Specifically, they constructed infinitesimally thin steady-state disk structures in the presence of a time-independent, nonaxisymmetric perturbing potential. While their problem was only two-dimensional and they did not seek a self-consistent solution of the gravitational Poisson equation, the approach they took to solving the 2D Euler equation in tandem with the continuity equation for a <i>compressible</i> fluid may very well be instructive. What follows is a summary of their approach.

===Governing Steady-State Equations===

As in our above [http://www.vistrails.org/index.php/User:Tohline/PGE/RotatingFrame#Example_of_Riemann_S-Type_Ellipsoids preamble to discussion of Riemann S-Type Ellipsoids], KP96 set <math>{\vec{\omega}} = \hat{k}\omega</math>. Hence, their steady-state Euler equation and steady-state continuity equation become (see their Eq. 1 or their Eq. 7),
<div align="center">
<math>
(\vec{v}\cdot \nabla)\vec{v} + 2\omega\hat{k}\times\vec{v} + \nabla \biggl[H + \Phi -\frac{1}{2}\omega^2 R^2 \biggr] = 0 ,
</math>

<math>
\nabla\cdot(\rho \vec{v}) = \vec{v}\cdot\nabla\rho + \rho\nabla\cdot\vec{v} = 0 .
</math>
</div>
Note that the KP96 notation is slightly different from ours:
* <math>\Sigma</math> is used in place of <math>\rho</math> to denote a two-dimensional <i>surface</i> density;
* <math>\Omega</math> is used instead of <math>\omega</math> to denote the angular frequency of the rotating reference frame;
* <math>W</math> is used instead of <math>H</math> to denote the fluid enthalpy; and
* <math>\Phi_g</math> represents the combined, time-independent gravitational and centrifugal potential, that is, <math>\Phi_g = (\Phi - \omega^2 R^2/2)</math>.

Using the <font color="darkgreen">vector identity</font>,
<div align="center">
<math>
(\vec{v}\cdot \nabla)\vec{v} = \frac{1}{2}\nabla(v^2) - \vec{v}\times(\nabla\times\vec{v}) ,
</math>
</div>
the above steady-state Euler equation can also be written as,
<div align="center">
<math>
2\omega\hat{k}\times\vec{v} - \vec{v}\times(\nabla\times\vec{v}) + \nabla \biggl[\frac{1}{2}v^2 + H + \Phi -\frac{1}{2}\omega^2 R^2 \biggr] = 0 .
</math>
</div>
Up to this point, no assumptions have been made regarding the behavior of the vector flow-field; we have only chosen to align the <math>\vec{\omega}</math> with the coordinate unit vector, <math>\hat{k}</math>. In particular, these derived forms for the steady-state Euler and continuity equations can serve to describe a fully 3D problem.

Before proceeding further we should emphasize that, in the context of the Euler equation written in this form (i.e., the form preferred by KP96), the vector <math>\vec{A}</math> defined in the preamble, above, should be written,
<div align="center">
<math>
\vec{A} = 2\omega\hat{k}\times\vec{v} +(\nabla\times\vec{v})\times\vec{v} + \nabla \biggl[\frac{1}{2}v^2 \biggr] .
</math>
</div>

===No Vertical Motions===
Now we restrict the flow by setting <math>v_z = 0</math>, that is, from here on we will assume that all the motion is planar. Also, following the lead of KP96, we define the vorticity of the fluid,
<div align="center">
<math>
\vec{\zeta} \equiv \nabla\times\vec{v} = \hat{i}\zeta_x + \hat{j}\zeta_y + \hat{k}\zeta_z .
</math>
</div>
[Note that (unfortunately) KP96 use <math>\omega</math> instead of <math>\zeta</math> to represent the rotating-frame vorticity.] In terms of the components of the vorticity vector, the steady-state Euler equation therefore becomes,
<div align="center">
<math>
(2\omega + \zeta_z)\hat{k}\times\vec{v} + (\hat{i}\zeta_x + \hat{j}\zeta_y)\times\vec{v} + \nabla \biggl[\frac{1}{2}v^2 + H + \Phi -\frac{1}{2}\omega^2 R^2 \biggr] = 0 .
</math>
</div>
Continuing to follow the lead KP96, we next <font color="red">take the curl of this Euler equation</font>. Because the curl of a gradient is always zero, this leads us to the same condition discussed above — but this time written in terms of the components of the vorticity — namely,
<div align="center">
<math>
\nabla\times\vec{A} = 0 = \nabla\times [(2\omega + \zeta_z)\hat{k}\times\vec{v} + (\hat{i}\zeta_x + \hat{j}\zeta_y)\times\vec{v}] .
</math>
</div>
Using another <font color="darkgreen">vector identity</font>, namely,
<div align="center">
<math>
\nabla\times(\vec{C} \times \vec{B}) = (\vec{B}\cdot\nabla)\vec{C} - (\vec{C}\cdot\nabla)\vec{B} + \vec{C}(\nabla\cdot\vec{B}) - \vec{B}(\nabla\cdot\vec{C}),
</math>
</div>
and remembering that we are assuming <math>v_z = 0</math>, we see in this case that the vector condition <math>\nabla\times\vec{A}=0</math> leads to the following three independent scalar constraints:
<div align="left">
<math>
~~~~~\hat{i}:~~~~~ [\nabla\times\vec{A}]_x = - \frac{\partial }{\partial z}\biggl[ (2\omega + \zeta)v_x \biggr] + \frac{\partial}{\partial y} \biggl[ \zeta_x v_y - \zeta_y v_x \biggr] = 0 ;
</math><br />
<math>
~~~~~\hat{j}:~~~~~ [\nabla\times\vec{A}]_y = - \frac{\partial }{\partial z} \biggl[ (2\omega + \zeta)v_y \biggr] - \frac{\partial}{\partial x} \biggl[ \zeta_x v_y - \zeta_y v_x \biggr] = 0 ;
</math><br />
<math>
~~~~~\hat{k}:~~~~~ [\nabla\times\vec{A}]_z = (2\omega + \zeta)\nabla\cdot\vec{v} + \biggl[ v_x \frac{\partial}{\partial x} + v_y \frac{\partial}{\partial y} \biggr](2\omega + \zeta) = 0 .
</math>
</div>
With the understanding that, by definition,
<div align="center">
<math>
\zeta_x \equiv - \frac{\partial v_y}{\partial z} , ~~~~~
\zeta_y \equiv + \frac{\partial v_x}{\partial z} , ~~~~~ \mathrm{and} ~~~~~
\zeta_z \equiv \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y} ,
</math>
</div>
it can be shown that these three constraints are identical to the ones presented in the preamble, above.

===Solution Strategy===

<font color="darkblue"><b>Constraint #1:</b></font>
For their two-dimensional disk problem, KP96 focused on the constraint provided by the z-component of the curl of the Euler equation, which can be rewritten as (see above derivation, or Eq. 2 of KP96),
<div align="center">
<math>
\nabla\cdot\vec{v} =-\vec{v} \cdot \biggl[ \frac{\nabla(2\omega + \zeta_z)}{(2\omega + \zeta_z)} \biggr]
= -\vec{v} \cdot \nabla[\ln(2\omega + \zeta_z)].
</math>
</div>
<font color="darkblue"><b>Constraint #2:</b></font>
But from the continuity equation they also know that,
<div align="center">
<math>
\nabla\cdot\vec{v} = -\vec{v}\cdot\biggl[\frac{\nabla\rho}{\rho} \biggr] = -\vec{v} \cdot \nabla[\ln\rho] .
</math>
</div>
Hence,
<div align="center">
<math>
\vec{v} \cdot \nabla[\ln(2\omega + \zeta_z)] = \vec{v} \cdot \nabla[\ln\rho] ,
</math>
</div>
that is,
<div align="center">
<math>
\vec{v} \cdot \nabla\ln\biggl[ \frac{(2\omega + \zeta_z)}{\rho} \biggr] = 0 .
</math>
</div>
This is essentially KP96's Eq. (3).

<font color="darkblue"><b>Introduce stream function:</b></font>
The constraint implied by the continuity equation also suggests that it might be useful to define a stream function in terms of the momentum density — instead of in terms of just the velocity, which is the natural treatment in the context of incompressible fluid flows. KP96 do this. They define the stream function, <math>\Psi</math>, such that (see their Eq. 4),
<div align="center">
<math>
\rho\vec{v} = \nabla\times(\hat{k}\Psi) .
</math>
</div>
in which case,
<div align="center">
<math>
v_x = \frac{1}{\rho} \frac{\partial \Psi}{\partial y} ~~~~~\mathrm{and}~~~~~
v_y = - \frac{1}{\rho} \frac{\partial \Psi}{\partial x} .
</math>
</div>
This implies as well that the z-component of the fluid vorticity can be expressed in terms
of the stream function as follows (see Eq. 5 of KP96):
<div align="center">
<math>
\zeta_z = - \nabla\cdot \biggl( \frac{\nabla\Psi}{\rho} \biggr) = - \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>
</div>

According to KP96, this expression, taken in combination with the conclusion drawn above from the second constraint — that is, Eq. (3) taken in combination with Eq. (4) from KP96 — "tell us that the 'vortensity' <math>(\zeta_z + 2\omega)/\rho</math> is constant along streamlines which are lines of constant <math>\Psi</math>." The vortensity is therefore a function of <math>\Psi</math> alone, so we can write,
<div align="center">
<math>
\frac{\zeta_z + 2\omega}{\rho} = g(\Psi) .
</math>
</div>

<font color="darkblue"><b>Constraint #3:</b></font>

Taking the scalar product of <math>\vec{v}</math> and the following form of the steady-state Euler equation,
<div align="center">
<math>
2\omega\hat{k}\times\vec{v} - \vec{v}\times(\nabla\times\vec{v}) + \nabla \biggl[\frac{1}{2}v^2 + H + \Phi -\frac{1}{2}\omega^2 R^2 \biggr] = 0 ,
</math>
</div>
we obtain the constraint,
<div align="center">
<math>
\vec{v}\cdot\nabla \biggl[\frac{1}{2}v^2 + H + \Phi -\frac{1}{2}\omega^2 R^2 \biggr] = 0 .
</math>
</div>
When tied with our earlier discussion, this means that the Bernoulli function also must be constant along streamlines. Hence, we can write,
<div align="center">
<math>
\frac{1}{2}v^2 + H + \Phi -\frac{1}{2}\omega^2 R^2 = F(\Psi) .
</math>
</div>
KP96 then go on to demonstrate that the relationship between the functions <math>g(\Psi)</math> and <math>F(\Psi)</math> is,
<div align="center">
<math>
\frac{dF}{d\Psi} = -g(\Psi) ,
</math>
</div>
which allows the determination of <math>F</math> up to a constant of integration.

===Summary===
In summary, KP96 constrain their flow as follows:
# They use the z-component of the curl of the Euler equation;
# They use the compressible version of the continuity equation;
# Instead of taking the divergence of the Euler equation to obtain a Poisson-like equation, they obtain an algebraic constraint on the Bernoulli function (as in our traditional SCF technique) by simply "dotting" <math>\vec{v}</math> into the Euler equation.

{{ SGFfooter }}

Apps/RiemannEllipsoidsCompressible - Revision history