User:Tohline/PGE/RotatingFrame

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 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:

<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>

Proposed Solution Strategy

Preamble:

Specify the three "polar" boundary locations, <math>a, b,</math> and <math>c</math>; specify the direction 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:

<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>

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 positional 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:

<math> \nabla^2 \Phi^* = 4\pi \rho^* . </math>

Step #2:

Use the continuity equation and the curl of the Euler equation to numerically derive the structure 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:

<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>

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:

<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>

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,

<math> \vec{A} \equiv ({\vec{v}}\cdot \nabla){\vec{v}} +2 {\vec{\omega}} \times {\vec{v}} . </math>

The components are:

<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>
<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>
<math> ~~~~~\hat{k}:~~~~~A_z = 0 . </math>

The curl of <math>\vec{A}</math> (needed in Step #2, above) produces a vector with the following three Cartesian components:

<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>
<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>
<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>

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:

<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>
<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>
<math> ~~~~~\hat{k}:~~~~~ \frac{\partial}{\partial x} \biggl[ C_{z1}(x,y) \biggr] = \frac{\partial}{\partial y} \biggl[C_{z2}(x,y) \biggr] , </math>

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,

<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>

The divergence of <math>\vec{A}</math> (providing the right-hand-side of the Poisson-like equation in Step #3, above) generates:

<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>
<math> = \frac{\partial}{\partial x} \biggl[ C_{z2}(x,y) \biggr] + \frac{\partial}{\partial y} \biggl[C_{z1}(x,y) \biggr] . </math>

Riemann Flow-Field

In Riemann S-Type ellipsoids, the adopted planar flow-field as viewed from the rotating reference frame is,

<math> \vec{v} = \hat{i} \biggl( \frac{\lambda a y}{b} \biggr) + \hat{j} \biggl( - \frac{\lambda b x}{a} \biggr) . </math>

Hence,

<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>
<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>

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,

<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>

However, the right-hand-side of our Poisson-like equation is not zero; rather, it is,

<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>

Summary

What can we learn from the Riemann S-Type ellipsoids? Well, let's assume that the structure 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 planar, 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:

Step #1 (simplified):

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:

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

Step #2 (simplified):

Use the continuity equation and the curl of the Euler equation to numerically derive the structure 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:

<math> \frac{\partial}{\partial x}\biggl[ \rho v_x \biggr] + \frac{\partial}{\partial y}\biggl[ \rho v_y \biggr] = 0 . </math>

<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>
<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>
<math> ~~~~~\hat{k}:~~~~~ \frac{\partial}{\partial x} \biggl[ C_{z1}(x,y) \biggr] = \frac{\partial}{\partial y} \biggl[C_{z2}(x,y) \biggr] , </math>

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.

Step #3 (simplified):

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:

<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>

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 algebraic condition:

<math> H + \Phi -\frac{1}{2}\omega^2 (x^2 + y^2) + f(x,y) = C_\mathrm{Bernoulli} , </math>

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)],

<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>

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, Korycansky & Papaloizou (1996, ApJS, 105, 181; hereafter KP96) developed a method to find nontrivial, nonaxisymmetric steady-state flows in a two-dimensional 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 compressible fluid may very well be instructive. What follows is a summary of their approach.

Governing Steady-State Equations

As in our above 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),

<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>

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 surface 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 vector identity,

<math> (\vec{v}\cdot \nabla)\vec{v} = \frac{1}{2}\nabla(v^2) - \vec{v}\times(\nabla\times\vec{v}) , </math>

the above steady-state Euler equation can also be written as,

<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>

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.

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. Among other things, this means that the fluid vorticity as viewed from the rotating frame of reference, <math>\vec{\zeta} \equiv \nabla\times\vec{v}</math>, only has a vertical component. That is, we can write,

<math> \nabla\times\vec{v} = \hat{k}\zeta . </math>

[Note that (unfortunately) KP96 use <math>\omega</math> instead of <math>\zeta</math> to represent the rotating-frame vorticity.] And the steady-state Euler equation becomes,

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

Now, taking the curl of this Euler equation we find the following governing relation:

<math> 0 = \nabla\times [(2\omega + \zeta)\hat{k}\times\vec{v}] . </math>

Using another vector identity, namely,

<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>

and remembering that we are assuming <math>v_z = 0</math>, we see that the curl of this Euler equation has the following Cartesian components, each of which must independently be zero:

<math> ~~~~~\hat{i}:~~~~~ - \biggl[ (2\omega + \zeta)\frac{\partial }{\partial z}v_x + v_x \frac{\partial}{\partial z}(2\omega + \zeta) \biggr] = 0  ; </math>
<math> ~~~~~\hat{j}:~~~~~ - \biggl[ (2\omega + \zeta)\frac{\partial }{\partial z}v_y + v_y \frac{\partial}{\partial z}(2\omega + \zeta) \biggr] = 0 ; </math>
<math> ~~~~~\hat{k}:~~~~~ (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>

© 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