The text has not been modified since this page was first written in 1996.
Copyright © 1996-2007 by
L. M. Stockman
- All Rights Reserved
A solar model is the solution to a set of equations describing the physical processes occurring within the Sun. These equations are often divided into two subsets, the structure equations and the chemical evolution equations. In the very simplest cases, these equations are coupled nonlinear first-order ordinary differential equations, but more realistic models usually involve coupled nonlinear first-order partial differential equations. In either case, the equations cannot be solved analytically and thus must be solved numerically. Therefore, a solar model is usually comprised of a set of tables describing the conditions (chemical composition, density, luminosity, mass, pressure, temperature, etc.) at different depths in the star.
The physical variables involved in these equations are luminosity, mass, pressure, radius, and temperature. Other quantities such as density, energy generation rate, and opacity, are derived from these five variables. The equations of stellar structure are usually written with mass as the independent variable, not radius. One reason for this is that most stellar models assume that the star has constant mass over a given time but not constant radius. Indeed, we expect the radius of a star to vary over its lifetime.
The equations of stellar structure in my model are
dL ( dT gamma - 1 T dP )
-- = epsilon - Cp ( -- - --------- - -- )
dM ( dt gamma P dt )
dP - G M
-- = -------
4
dM 4 pi R
dR 1
-- = -----------
2
dM 4 pi R rho
dT - 3 kappa L
-- = ---------------- (radiative energy transfer)
2 4 3
dM 64 pi a c R T
dT gamma - 1 T dP
-- = --------- - -- (convective energy transfer)
dM gamma P dM
where L is the luminosity, M is the mass, P is the pressure, R is the radius, T is the temperature, t is the time, epsilon is the energy generation rate (from nuclear fusion), kappa is the Rosseland mean opacity, and rho is the density. Other physical constants have their usual meanings.
The equations of chemical evolution depend solely upon which nuclear processes are running in the interior of the star. In the Sun, the most important fusion process is the proton-proton chain (pp chain) which fuses four hydrogen atoms to one helium atom. However, another reaction network called the carbon-nitrogen-oxygen cycle (CNO cycle) is also active in the Sun, although it is much less important than the pp chain. In my model I keep track of five elements: hydrogen, two species of helium, carbon, and nitrogen. After making certain simplifying assumptions, the equations of chemical evolution in my model are
d X1 A1 ( 2 2
---- = -- ( -3 R11 X1 + 2 R33 X3 - R34 X2 X4 - 2 R112 X1 X12
dt N (
) d ( d X1 )
- 2 R114 X1 X14 ) + -- ( sigma ---- )
) dM ( dM )
d X3 A3 ( 2 2 ) d ( d X3 )
---- = -- ( R11 X1 - 2 R33 X3 - R34 X3 X4 ) + -- ( sigma ---- )
dt N ( ) dM ( dM )
d X4 A4 ( 2 ) d ( d X4 )
---- = -- ( R33 X3 - R34 X3 X4 + R114 X1 X14 ) + -- ( sigma ---- )
dt N ( ) dM ( dM )
d X12 A12 ( ) d ( d X12 )
----- = --- ( -R112 X1 X12 + R114 X1 X14 ) + -- ( sigma ----- )
dt N ( ) dM ( dM )
d X14 A14 ( ) d ( d X14 )
----- = --- ( R112 X1 X12 - R114 X1 X14 ) + -- ( sigma ----- )
dt N ( ) dM ( dM )
where RIJ is the reaction rate between species I and species J, X1 is the hydrogen abundance, X3 is the helium 3 abundance, X4 is the helium 4 abundance, X12 is the carbon 12 abundance, X14 is the nitrogen 14 abundance, and sigma is a collection of constants which includes a diffusive term.
The total abundance (the sum of all the abundances) must equal one so each abundance is a number between zero and one.
The energy generation rate can be described by something as simple as a power law to more complicated expressions involving nuclear reation rates. My more advanced solar model uses the equation
-5 2 -5 2
epsilon = 1.069 × 10 R11 X1 + 2.060 × 10 R33 X3
-5 -5
+ 3.043 × 10 R34 X3 X4 + 1.763 × 10 R112 X1 X12
-5
+ 2.246 × 10 R114 X1 X14
Like the energy generation rate, opacity can be described by a power law but more complicated models use opacity tables. I use the OPAL opacity tables and then perform a four-dimensional interpolation procedure across the logarithms of the hydrogen abundance, metals abundance, density, and temperature.
For the purposes of my model, the ideal gas equation of state
mu P
rho = -- -
Rg T
is used in the calculation of density. The variable mu is the mean molecular weight and is given by
84
mu = -----------------------------------------------
168 X1 + 84 X3 + 63 X4 + 49 X12 + 48 X14 + 42 Z
where Z is the metals abundance (not including carbon 12 or nitrogen 14).
Because no theory of convection exists, the problem of convective versus radiative transport of energy is treated very simply. The Schwarzschild condition is used to determine whether the radiative or convective transport derivative is used for each calculation. The amount of overshooting past the edge of any convective zone is handled simply: a percentage of the mass of the convective zone. The more overshooting there is, the longer the initial convective core lasts.
A very simple diffusion function is used. The diffusion constant is taken to be a very large value in convective regions and zero in radiative regions. The diffusion function is
2 2 4
16 pi rho R nu
sigma = -----------------
2
Msun
where nu is the diffusion constant.
Effects from magnetism, rotation, and non-sphericity are not considered.
The central boundary conditions for the stellar structure equations are taken to be L = 0, P = central pressure, R = 0, and T = central temperature. The central pressure and temperature are unknown and their values must be estimated. The surface boundary conditions for the stellar structure equations are taken to be L = surface luminosity, P = 0, R = surface radius, and T = 0. Again, the surface luminosity and radius are unknown and their values must be estimated.
For the chemical evolution equations, the diffusion is taken to be zero at the center and surface.
The initial conditions that need to be specified are the initial chemical abundances, and estimates for the central pressure and temperature, and surface luminosity and radius. A homogeneous chemical abundance distribution is the most usual configuration. My model begins with a homogeneous distribution of X1 = 0.71, X3 = 0.00001, X4 = 0.27, X12 = 0.2Z, and X14 = 0.05Z where Z = 1.0 - X1 - X3 - X4.
For a given set of initial conditions, the equations of stellar structure are solved to give a model at a given time. Once the model has been calculated, the chemical abundances are evolved one time step. Then a new model can be calculated for the evolved chemistry. This new model is then used to improve the chemical abundances for that time step and a final model for the given time step is calculated. This predictor-corrector procedure continues until the desired age of the star is reached.
It is clear that for the given boundary conditions, the equations of stellar structure are singular at both the center and the surface. Therefore, it is necessary to use special techniques at both boundaries to safely integrate away from them. My model uses truncated series approximations to the solutions at the center and surface.
The two-point boundary value problem involving the equations of stellar structure is solved using a shooting method. My ordinary differential equation integrator is the classical fourth-order Runge-Kutta algorithm. The initial model is calculated using adaptive stepsize control so that the mass grid has a fine mesh near the center and the surface. (The adaptive scheme is related to the scale heights of the four physical variables.) Since the integrator cannot be used at the center or the surface, special routines described above are used to take the first step out from the center and in from the surface. The shooting method involves integrating out from the center and in from the surface to some fitting point. If the correct values for central pressure and temperature, and surface luminosity and radius have been used, then the values for luminosity, pressure, radius, and temperature will agree at the fitting point for both integrations. However, if the values do not agree, a multidimensional Newton-Raphson method is used to find the corrections to the unknown boundary conditions that will hopefully improve the fit. This process of "shooting," comparing the fit, and correcting the boundary conditions continues until the values agree at the fitting point.
Technically, the luminosity equation is not an ordinary differential equation but a partial differential equation when entropy terms are included. The time derivatives are simply written as backward differences.
The evolution equations are nonlinear partial differential equations. The nuclear source terms are linearized according to the scheme
k+1 k+1 k+1 k k k+1 k k X X ~ X X + X X - X X i j i j i j i j
where k denotes the time step and i and j denotes the species. When diffusion is included, tridiagonal systems result for each species. The evolved value of helium 3 is solved for first, followed by carbon 12. Nitrogen 14 is calculated next in such a way as to guarantee that the total abundance of carbon 12 and nitrogen 14 is constant. Finally, the hydrogen abundance is computed. Helium 4 is obtained by subtracting the known quantities from one. Once approximations for the new abundances are known, the equations are solved again and again until the improvement is negligible. Only then is the next model calculated.
The text has not been modified since this page was first written in 1996.
Copyright © 1996-2007 by
L. M. Stockman
- All Rights Reserved