Difference between revisions of "User:Tohline/AxisymmetricConfigurations/PoissonEq"

From VistrailsWiki
Jump to navigation Jump to search
Line 56: Line 56:


Also, following the lead of [http://adsabs.harvard.edu/abs/1975ApJ...199..619B Black & Bodenheimer (1975)], a hybrid scheme was developed to solve the Poisson equation.  As defined immediately above, the ''integral'' expression was evaluated at grid locations along the outermost radial and vertical boundaries by integrating over the mass distribution everywhere inside of these grid boundaries (see, for example, the schematically illustrated "pink" mass distribution in the annotated figure, below).  Separately, the ''differential'' expression was used to evaluate the potential at all ''interior'' grid locations; this differential expression was supplemented by the implementation of ''Neumann'' boundary conditions (reflection symmetry) along the equatorial plane, and by using the values of the potential just determined along the outermost grid boundaries to provide ''Direchlet'' boundary conditions along those gird boundaries.
Also, following the lead of [http://adsabs.harvard.edu/abs/1975ApJ...199..619B Black & Bodenheimer (1975)], a hybrid scheme was developed to solve the Poisson equation.  As defined immediately above, the ''integral'' expression was evaluated at grid locations along the outermost radial and vertical boundaries by integrating over the mass distribution everywhere inside of these grid boundaries (see, for example, the schematically illustrated "pink" mass distribution in the annotated figure, below).  Separately, the ''differential'' expression was used to evaluate the potential at all ''interior'' grid locations; this differential expression was supplemented by the implementation of ''Neumann'' boundary conditions (reflection symmetry) along the equatorial plane, and by using the values of the potential just determined along the outermost grid boundaries to provide ''Direchlet'' boundary conditions along those gird boundaries.
An extension from 2D to 3D was accomplished by using separate, finite Fourier series expansion to represent the azimuthal dependence of <math>~\rho</math> and <math>~\Phi</math>.  This allowed us to decouple the ''differential'' expression for the Poisson equation into a finite number of independent Fourier modes, and to straightforwardly solve a finite set of 2D Helmholtz equations.


==Constructing Two-Dimensional, Axisymmetric Structures==
==Constructing Two-Dimensional, Axisymmetric Structures==

Revision as of 05:06, 10 April 2018

Solving the Poisson Equation Numerically

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


The set of Principal Governing Equations that serve as the foundation of our study of the structure, stability, and dynamical evolution of self-gravitating fluids contains an equation of motion (the Euler equation) that includes an acceleration due to local gradients in the (Newtonian) gravitational potential, <math>~\Phi</math>. As has been pointed out in an accompanying chapter that discusses the origin of the Poisson equation, the mathematical definition of this acceleration is fundamentally drawn from Isaac Newton's inverse-square law of gravitation, but takes into account that our fluid systems are not point-mass sources but, rather, are represented by a continuous distribution of mass via the function, <math>~\rho(\vec{x},t)</math>. As indicated, in our study, <math>~\rho</math> may depend on time as well as space. The acceleration felt at any point in space is obtained by integrating over the accelerations exerted by each differential mass element. As has been explicitly demonstrated in, respectively, Step 1 and Step 3 of the same accompanying chapter, at any point in time the spatial variation of the gravitational potential, <math>~\Phi(\vec{x})</math>, may be determined from <math>~\rho(\vec{x})</math> via either an integral or a differential equation as follows:


Poisson Equation
Integral Expression Differential Expression

<math>~ \Phi(\vec{x})</math>

<math>~=</math>

<math>~ -G \int \frac{\rho(\vec{x}^{~'})}{|\vec{x}^{~'} - \vec{x}|} d^3x^' \, .</math>

LSU Key.png

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


While it is possible in some restricted situations to determine analytic expressions for the matched pair of functions, <math>~\Phi </math> and <math>~\rho</math>, that satisfy the Poisson equation, modeling the vast majority of interesting astrophysical problems requires the develop of a numerical scheme to solve the Poisson equation. In what follows, our aim is twofold: (a) To recount — in a reasonable amount of detail — the steps that we have taken over the past, approximately forty years to develop more and more accurate and efficient ways to solve the Poisson equation in full three-dimensional generality; and (b) To summarize alternative techniques that have been successfully employed by other research groups over the years.

Our Approach

Influenced by Black & Bodenheimer (1975)

D. C. Black & P. Bodenheimer (1975, ApJ, 199, 619 - 632) published a detailed description of the numerical techniques that they implemented in order to model, in two dimensions, the dynamical collapse of axisymmetric interstellar gas clouds. Among the techniques was their method of solving, in cylindrical coordinates, the two-dimensional (2D) Poisson equation. As is acknowledged in the last paragraph of their paper, Black & Bodenheimer were introduced "to the pitfalls of two-dimensional hydrodynamics" by Drs. J. LeBlanc and J. Wilson who, at the time, were both staff scientists at the Lawrence Livermore Laboratory. In 1976, having completed my formal graduate-level course work in the astronomy program at UC, Santa Cruz, Peter Bodenheimer (on the faculty of Lick Observatory and UC, Santa Cruz) asked me if I would be interested in working with him and David Black (a planetary scientist at NASA, Ames Research Center) on the development of a fully three-dimensional hydrodynamics code to study, not only the collapse, but also the fragmentation of self-gravitating gas clouds. I jumped at the opportunity. As a result, an effort to model the process of spontaneous cloud fragmentation became the focus of my doctoral dissertation research.

Borrowing from the Black & Bodenheimer (1975) work, I decided to solve the set of 3D principal governing equations on a cylindrical coordinate <math>~(R, \theta, Z)</math> mesh:

  • whose outermost radial and vertical boundaries were placed entirely outside of the mass distribution (by way of illustration, see the green-dashed lines in the annotated figure, below);
  • that exhibited reflection symmetry through the <math>~(Z = 0)</math> equatorial plane ;
  • that allowed for non-uniform (logarithmically stretched) grid spacing in both the radial and vertical directions; and,
  • (extending the Black & Bodenheimer work from 2D to 3D) with strictly uniform spacing in the azimuthal-coordinate direction.

Also, following the lead of Black & Bodenheimer (1975), a hybrid scheme was developed to solve the Poisson equation. As defined immediately above, the integral expression was evaluated at grid locations along the outermost radial and vertical boundaries by integrating over the mass distribution everywhere inside of these grid boundaries (see, for example, the schematically illustrated "pink" mass distribution in the annotated figure, below). Separately, the differential expression was used to evaluate the potential at all interior grid locations; this differential expression was supplemented by the implementation of Neumann boundary conditions (reflection symmetry) along the equatorial plane, and by using the values of the potential just determined along the outermost grid boundaries to provide Direchlet boundary conditions along those gird boundaries.

An extension from 2D to 3D was accomplished by using separate, finite Fourier series expansion to represent the azimuthal dependence of <math>~\rho</math> and <math>~\Phi</math>. This allowed us to decouple the differential expression for the Poisson equation into a finite number of independent Fourier modes, and to straightforwardly solve a finite set of 2D Helmholtz equations.

Constructing Two-Dimensional, Axisymmetric Structures

As has been explained in an accompanying discussion, our objective is to solve an algebraic expression for hydrostatic balance,

<math>~H + \Phi + \Psi = C_0</math> ,

in conjunction with the Poisson equation in a form that is appropriate for two-dimensional, axisymmetric systems, namely,

<math>~ \frac{1}{\varpi} \frac{\partial }{\partial\varpi} \biggl[ \varpi \frac{\partial \Phi}{\partial\varpi} \biggr] + \frac{\partial^2 \Phi}{\partial z^2} = 4\pi G \rho . </math>

Steps to Follow

Annotation of
Figure 1 from I. Hachisu (1986)

HSCF Meridional Plane
One quadrant of a meridional-plane cross-section (pink) through: (a) spheroidal structure; (b) toroidal structure. A green, dashed rectangular grid boundary is also illustrated.
  1. Choose a particular barotropic equation of state.   More specifically, functionally define the density-enthalpy relationship, <math>~\rho(H)</math>, and identify what value, <math>~H_\mathrm{surface}</math>, the enthalpy will have at the surface of your configuration. For example, if a polytropic equation of state is adopted, <math>~H_\mathrm{surface} = 0</math> is a physically reasonable prescription.
  2. Choosing from, for example, a list of astrophysically relevant simple rotation profiles, specify the corresponding functional form of the centrifugal potential, <math>~\Psi(\varpi)</math>, that will define the radial distribution of specific angular momentum in your equilibrium configuration. If the choice is uniform rotation, then <math>~\Psi = - \varpi^2 \omega_0^2/2 \, ,</math> where <math>~\omega_0</math> is a constant to be determined.
  3. On your chosen computational lattice — for example, on a cylindrical-coordinate mesh — identify two boundary points, A and B, that will lie on the surface of your equilibrium configuration. These two points should remain fixed in space during the HSCF iteration cycle and ultimately will confine the volume and define the geometry of the derived equilibrium object. Note that, by definition, the enthalpy at these two points is, <math>~H_A = H_B = H_\mathrm{surface}</math>.
  4. Throughout the volume of your computational lattice, guess a trial distribution of the mass density, <math>~\rho(\varpi,z)</math>, such that no material falls outside a volume defined by the two boundary points, A and B, that were identified in Step #3. Usually an initially uniform density distribution will suffice to start the SCF iteration.
  5. Via some accurate numerical algorithm, solve the Poisson equation to determine the gravitational potential, <math>~\Phi(\varpi,z)</math>, throughout the computational lattice corresponding to the trial mass-density distribution that was specified in Step #4 (or in Step #9).
  6. From the gravitational potential determined in Step #5, identify the values of <math>~\Phi_A</math> and <math>~\Phi_B</math> at the two boundary points that were selected in Step #3.
  7. From the "known" values of the enthalpy (Step #3) and the gravitational potential (Step #6) at the two selected surface boundary points A and B, determine the values of the constants, <math>~C_0</math> and <math>~\omega_0</math>, that appear in the algebraic equation that defines hydrostatic equilibrium.
  8. From the most recently determined values of the gravitational potential, <math>~\Phi(\varpi,z)</math> (Step #5), and the values of the two constants, <math>~C_0</math> and <math>~\omega_0</math> just determined (Step #7), determine the enthalpy distribution throughout the computational lattice.
  9. From <math>~H(\varpi,z)</math> and the selected barotropic equation of state (Step #1), calculate an "improved guess" of the density distribution, <math>~\rho(\varpi,z)</math>, throughout the computational lattice.
  10. Has the model converged to a satisfactory equilibrium solution? (Usually a satisfactory solution has been achieved when the derived model parameters — for example, the values of <math>~C_0</math> and <math>~\omega_0</math> — change very little between successive iterations and the viral error is sufficiently small.)
    • If the answer is, "NO":   Repeat steps 5 through 10.
    • If the answer is, "YES":   Stop iteration.

Related Discussions

Reviews

Solution Methods

Early Eriguchi Applications

Other Example Applications

Henyey Technique for Nonrotating Stars


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

Retrieved from "https://www.vistrails.org//index.php?title=User:Tohline/AxisymmetricConfigurations/PoissonEq&oldid=15100"