Brief description:
This is a description of the equations I used in the star modelling program
and the methods used to solve them. We did not go into this in much detail
in class, so about all I can offer is what we did do. The program takes
the user's input for the following values of the star: mass, surface
luminosity, surface pressure, surface temperature, surface density, radius,
surface fraction of H, surface fraction of He, surface fraction of metals,
core fraction of H, core fraction of He, core fraction of metals, ratio of
specific heats, number of steps to take, whether or not to use CNO cycle.
The options of the core values for H, He, and metals and use of the CNO
cycle are added in order to help increase the accuracy of the program.
However, sometimes better results may be achieved without using them.
Using these values, the program uses sets of equations in order to proceed
inward towards the center of the star, calculating the values of mass M(r),
luminosity L(r), pressure P(r), temperature T(r), density rho(r), energy
generated epsilon(r), and opacity kappa(r) at each shell and using the new
values to calculate the values for the next shell in. The process repeats
until the center has been reached or the mass and/or luminosity become less
than zero before the center is reached. The rest of this file describes
the equations and the way they were given to us in class and then goes
into a more detailed description of how the calculations are made.
Conservation of mass:
---------------- M(r + delta_r)
/|\
| delta_r
\|/
---------------- M(r)
M(r + delta_r) - M(r) = mass inside one shell
M(r + delta_r) - M(r) = rho(r) * V <-- Substitute in for volume V
M(r + delta_r) - M(r) = rho(r) * 4 * pi * r^2 * delta_r
M(r + delta_r) - M(r)
--------------------- = rho(r) * 4 * pi * r^2 <-- left side is definition
delta_r of the derivative
dM(r)
----- = rho(r) * 4 * pi * r^2
dr
Conservation of energy:
---------------- L(r + delta_r)
/|\
| delta_r
\|/
---------------- L(r)
L(r + delta_r) = L(r) + (epsilon(r) * rho(r) * V)
solving as in the mass equation above,
dL(r)
----- = epsilon(r) * rho(r) * 4 * pi * r^2
dr
Hydrostatic equilibrium:
Fa |
\|/
/----------/
/m area=a/|
/----------/ |delta_r
| |f | |
| \|/ | /
|----------|/
/|\
| Fb
Fa + f = Fb <-- looking at forces acting on a slab of the shell
Fb - Fa = f = G * M(r) * m / r^2 <-- forces are gravitational
f = G * M(r) * rho(r) * V / r^2
using P = F / a and F = P * a
f = G * M(r) * rho(r) * a * delta_r / r^2 = a * (Pa - Pb)
Pb - Pa G * M(r) * rho(r)
------- = -----------------
delta_r r^2
dP(r) -G * M(r) * rho(r) must be negative because as r increases,
----- = ------------------ P decrease and equation must match reality
dr r^2
Temperature:
{ -3 kappa(r) * rho(r) L(r)
{ --------- * ----------------- * ---- <-- radiative
dT(r) { 16*pi*a*c T^3 r^2
----- = {
dr { ( 1 ) T(r) dP(r)
{ (1 - -----) * ------ * ----- <-- convective
{ ( gamma) rho(r) dr
Here, a is the cgs equivalent of the Stefan-Boltzman constant sigma
and gamma is the ratio of specific heats. The lower value of the two
must be used, so the program calculates both, checks to see which is
lower, and uses the lower value in the rest of the calculations.
Perfect gas law:
P(r) = N * k * T(r) N = # particles / volume ; k = Boltzmann's constant
X = fraction by mass of H Y = fraction by mass of He,
Z = fraction by mass of metals and X + Y + Z = 1
H He metals
# nuclei X*rho/H Y*rho/4H Z*rho/AH <-- negligible
# electrons X*rho/H 2Y*rho/4H AZ*rho/2HA = Z*rho/2H
k * T(r) * rho(r) ( 3 1 )
P(r) = ----------------- * ( 2X + - * Y + - * Z )
H ( 4 2 )
Other Equations:
Z * (1 + X) rho(r)
kappa(r) = 4.34E+25 * ----------- * -------
3.162 T^(3.5)
( 1E+6 )^(2/3) ( (1E+6)^(1/3))
epsilon(r) = 2.5E+6 * rho(r) * X^2 * ( ---- ) * exp (-33.8*(----) )
[pp chain] ( T ) ( ( T ) )
XZ ( 1E+6 )^(2/3) ( (1E+6)^(1/3))
epsilon(r) = 9.5E+28 * rho(r) * -- * ( ---- ) * exp (-152.3*(----) )
[CNO cycle] 3 ( T ) ( ( T ) )
Method of solution:
The process of solving the equations is started by calculating the delta_r
for the equations. This is simply the radius of the star divided by the
number of steps the user wants to take (delta_r = radius/#steps). The
initial values for epsilon and kappa are also calculated. A loop is
started which counts the number of steps taken. The inner radius of the
shell is calculated by (outer radius - delta_r) and it is used to calculate
the new values for mass, luminosity, pressure, temperatue, density, energy
generated, and opacity. Most of these equations are differential ones and
in order to solve them, we used the definition of the derivative rather
than integrating the expressions. That is, we have an expression for
dQ/dr. By definition, dQ/dr = ( Q(r + delta_r) - Q(r) ) / delta_r and
solving for Q(r) gives, Q(r) = -((delta_r * dQ/dr) - Q(r + delta_r)).
These values are printed to a file and to the screen after which they are
used as the new initial values. The program checks to see if the energy
generated has become greater than zero which indicates the core region has
been reached. When the program reaches the core region, it switches from
using the surface values for H, He, and metals to the core values. This
gives the program better accuracy. It then checks to see if the mass or
luminosity has become less than zero before the center of the star is
reached. If so, the computer stops calculating and informs the user that
the center was not reached and that the data file is not fully complete.
Otherwise, the program loops back to the point where the inner radius is
calculated and continues running. If the center of the star is reached,
the program notifies the user that it was able to complete the calculations
and a full set of data is in the file. The program pauses when the screen
gets full so that the user can see the values being calculated before they
scroll past. The program prints the step number it is at along with the
rest of the values, so the user can tell how far into the interior the
program has reached. Using the default values for the Sun, which give the
best results (that I was able to come up with), the program can get 95% of
the way to the center. At step number 96, both the mass and luminosity
become negative. My professor at the time I wrote this program for class
told me that everyone who has written this program for him has obtained
approximately the same results -- about 95% of the way into the Sun.
I hope that this file has given you enough information to understand how the
program works and how a simple stellar model can be built. This is all of
the information I have on this subject right now, but if I can be of help,
please feel free to contact me at my current Internet address below. I will
do my best to help you, but remember that I am a student and still have a
lot to learn. Enjoy the program!
View standard disclaimer for my personal web pages.