PDE Project

Garrett Shuldes
April 17, 2021

Stars are commonplace throughout the universe. Whether you look towards
the plane of our galaxy, towards the galaxies that surround us, or towards the
bright glowing orb that rises and sets every day, there is no question that there
is an abundance of stars. Thanks to our very own Sun, we know quite a bit
about how these stars live and die, and what processes happen within the outer
bounds of the star to make it so hot, and to make it live as long as it does. But
what would happen if stars worked in a different fashion?

1 Introduction
Stars such as our Sun go through cyclical phases of activity levels. These cycles
are created through the continuous fluctuation of the magnetic poles within the
core. Every 11 years, about, the magnetic alignment of the Sun reverses, and so
every 11 years there is a change in the activity levels of the Sun. There are far
more sunspots that appear on the surface due to a more active magnetic field,
there are more coronal mass ejections, and there are far more coronal loops
within this active period. No matter what, however, the process of nuclear
fusion occurs within the immensely dense, high pressured core of the Sun, and
this process both heats the Sun, and creates all of the energy that is put out,
whether it be in the form of radio, micro, infrared, visible, or ultraviolet light.
The question that I aim to focus on throughout this project is, then, what
would happen if the processes were combined? What would happen if the cycles
of solar activity were caused by a fluctuating energy source in the core? What
would happen to the core if it were to be fluctuating at a consistent rate? Would
the different layers of the Sun stay at a temperature consistent with what we
believe them to be? Do all of the layers tend toward the same temperature? Is
there an equilibrium temperature?
In order to answer any of these questions, a model was needed that would
be able to look at the center source, the temperatures of each layer, the abil-
ity to create a distinct barrier between the layers, and the ability to track the
trends that a specific point of the modeled star takes as time progresses. Using a
Jupyter Notebook and a partial differential heat equation with non-uniform ma-
terials and layers, a model of star very similar to our own Sun was created, with

1
a matching activity cycle period, matching layer temperatures, and matching
layer radii.

2 The Model
There are three parts to this stellar model. The first part contains the formation
and structure of the model itself. This includes the defined function, the setup
for that function, and for the general visualization of what the modeled star
looks like, the second part contains the actual calculations and data recording
that is required to answer the proposed questions, and the third part contains
the analysis and representation of the results of the computations.

2.1 Formation and Structure

Stars are large. Large enough to make the numbers being worked with unap-
pealing, and so for the purposes of this model, the physical size of the Star was
reduced by a factor of 100,000, and the physical representation of the model
looks as such,

There are 3 distinct layers shown, with an additional 2 bunched together

at the very edges. This plot does nothing if not to show what the system as a
whole looks like, without the additional data that will be presented.

2
2.1.1 The Function
Following this visual of what the system looks like, then, is the function that will
be run over a specified time period, and will provide the data to be analyzed and
interpreted. The function is of the 2-dimensional heat equation in a non-uniform
domain. One possible representation of this in code is as follows,

def fheat(s,t, dx,nx,dy,ny):

k = layer_boundaries
ka = zeros((ny+2,nx+2))
ka[1:-1,1:-1] = k

T = s.reshape((ny,nx))
Ta = zeros((ny+2,nx+2))
Ta[1:-1,1:-1] = T

gradientk = (ka[2:,1:-1]-ka[:-2,1:-1])/2*dx + (ka[1:-1,2:]-ka[1:-1,:-2])/2*dy

gradientT = (Ta[2:,1:-1]-Ta[:-2,1:-1])/2*dx + (Ta[1:-1,2:]-Ta[1:-1,:-2])/2*dy
laplacianT = k*((Ta[2:,1:-1]-2*Ta[1:-1,1:-1]+Ta[:-2,1:-1])/dy**2
+ (Ta[1:-1,2:]-2*Ta[1:-1,1:-1]+Ta[1:-1,:-2])/dx**2)
source = S*(2+cos(t/11))

dTdt = (gradientk * gradientT) + laplacianT + source

return dTdt.reshape(ny*nx)

where there are 6 arguments passed into the function including; the state of the
function, the time that the function will be run over, and then the scope of the
system in both the x and y directions, as well as the difference between each
point to be considered. Because the domain is non-uniform, the boundaries of
each layer depend on position, the chain rule dictates that the gradient of both
the state of the system (the temperature) and the boundaries must be added to
the heat equation. This requires that both the state and the boundaries have
an augmented matrix, so that the finite difference approximations for each can
be calculated. The returned value, then, is the calculated gradients, laplacian,
and added source term.

2.1.2 System Setup and Structure

In order for this function to work, however, there is background information and
context that the function requires. To begin with, the heat transfer boundaries
were setup, and look like such,

nx = 41; ny = 41; halfnx = int(nx/2); halfny = int(ny/2)

x = linspace(-10,10,nx); y = linspace(-10,10,ny)
xx,yy = meshgrid(x,y)
dx = x[1]-x[0]; dy = y[1]-y[0]

3
#Each layer of a Star is made up of different materials/densities
rr = sqrt(xx**2+yy**2)
core = tanh(5*(rr-1.5))
radiative = tanh(5*(rr-4.5))
convective = tanh(5*(rr-6.5))
photosphere = tanh(5*(rr-6.505))
chromosphere = tanh(5*(rr-6.605))

layers = -core - 0.5*radiative - 0.1*convective - 0.00025*photosphere

+ 0.0005*chromosphere

#Loop to bring tanh functions to a minimum of 0

minimum_layer = zeros(nx)
for i in range(nx):
minimum_layer[i] = min(layers[i])
minimum = min(minimum_layer)
if sign(minimum) == -1.0:
minimum = -minimum

layer_boundaries = -core - 0.5*radiative - 0.1*convective - 0.00025*photosphere

+ 0.0005*chromosphere + minimum

where the hyperbolic tangent functions are used to determine the difference in
heating rates between each of the layers, as near the core, heat is transferred
largely through radiation, and near the outer layers of the Sun, heat is trans-
ferred largely through convection. The value of layer boundaries is larger at the
core because radiation typically tends to transfer heat faster than either con-
vection or conduction. The image below gives a 2-dimensional representation of
the boundaries.

4
 Following the setup of the heat transfer boundaries comes the setup of the
initial conditions, which appear as such,

arr = zeros((ny,nx))
corearr = zeros((ny,nx))

#Create indices for temperatures

chromosphere_mask = xx**2 + yy**2 < chromosphere_radius**2
photosphere_mask = xx**2 + yy**2 < photosphere_radius**2
convective_mask = xx**2 + yy**2 < convective_radius**2
radiative_mask = xx**2 + yy**2 < radiative_radius**2
core_mask = xx**2 + yy**2 < core_radius**2

#Initial temperatures of each layer

arr[chromosphere_mask] = 6300
arr[photosphere_mask] = 5800
arr[convective_mask] = 2e6
arr[radiative_mask] = 7.5e6
arr[core_mask] = 1.5e7

#Initial Conditions

5
corearr[core_mask] = 1.5e7
S = corearr*exp(-xx**4 - yy**4)
ts = linspace(0,100,101)
ts_2 = linspace(0,1100,101)
ts_3 = linspace(0,2200,101)
T0 = arr

where arr is simply the array that contains the initial temperature of the star.
The masks that are created represent a circle of indices within the radius of each
of the specified layers, and the values in arr are then set to the temperature
of the layer, from outermost to innermost. The temperatures must be set from
chromosphere first to core last because this allows each of the layers to be set
to their specific temperature, without overlap from an outer layer.
For the initial conditions used, corearr is defined so that the source, S
comes specifically from the core, and nowhere else, the three time series allow
for multiple lengths of time to be analyzed at the same time, and T 0 is used
for clarity as the initial temperature state of the system. Now that we have all
of the required background systems set up, and the structure of our model has
been created, the function can be ran, and the results visualized.

2.2 Results and Visualization

Using the odeint function, imported from scipy, the above function, f heat, is
able to be ran through for each of the given time series, and the evolution of
the system is able to be visualized. Passing the function through odeint and
assigning the results to a variable appears as such,

ans = odeint(fheat,T0.reshape(ny*nx),ts,(dx,nx,dy,ny))
ans2 = odeint(fheat,T0.reshape(ny*nx),ts_2,(dx,nx,dy,ny))
ans3 = odeint(fheat,T0.reshape(ny*nx),ts_3,(dx,nx,dy,ny))

where the f heat is called as the function to be run through, T 0 is called as the
variable to represent the initial state of the system, ts is called as the time series
that is used to determine the length of time that the system is integrated for,
and the variables in parentheses are extra arguments passed into the function.
The results of odeint are best represented through a 2-dimensional image
plot, that contains a colorbar to help interpret the values of the system. The
following figures shows the state of the system at the final time steps, along with
two points along the surface of the system. The first image shows the final state
of the first time series, the second image shows the final state of the second time
series, and the third image shows the final state of the third time series.

6
 For this first plot, the time that the system was ran for ranged from 0
to 100 intervals. Here, we can see that there is a distinct difference between
the temperatures of the system, but that there are distinct boundaries, as was
described by layer boundaries, though there is a distinct point where the system
ends and the space around it begins. The system has just reached it’s minimum
in activity, and so is the coolest that it will ever be. Here, we can see that there
seems to be something a bit weird going on, as the top right and bottom left
seem to be a bit warmer than the top left and bottom right, even though the
system was set up completely symmetric. The two scatter points are located
at points that are on the very edge, and have two very different values for
temperature. The temperature of the system at the point in red is 465, 031,
whereas the temperature of the system at the point in orange is 986, 402, more
than double that of the point in red.

7
 This second plot ran for 1100 intervals, and we can see that the final state
of the system is far different than the final state of the previous system. The
system at this state, is just beginning to come back from a minimum, and so is
colder than it’s average, but not as cold as it could be. The final temperatures,
even in this colder state, however, are higher even than the starting conditions,
and the disparity between the two identified points is far more pronounced. The
heat seems to be transferring from the core to the upper right and lower left far
more effectively than it is to the upper left and lower right, which again, is quite
puzzling as the system was set up completely symmetric. The temperature of
the system at the point in red is 800, 188 and the temperature of the system at
the point in orange in 3, 352, 748, now more than four times as hot. What, then,
happens if the time interval is increased even further? Does the heat continue
to spread evenly to the upper right and lower left? Does the upper left and
lower right eventually begin to be heated as well?

8
 As it turns out, no, the upper left and lower right do not begin to be heated
equally, though the heat system is trending towards equilibrium as the the upper
right and lower left are rapidly rising in temperature. The system shown was
ran for 2200 intervals, and at this point, is approaching a minimum in activity.
Again, however, though the maximum temperature is lower than the previous
image, the overall temperature of the system is increasing steadily. The disparity
between the two points identified in the system are now drastic, as the point in
red has a temperature of 1, 003, 187 and the point in orange has a temperature
of 4, 680, 577.
The system is obviously trending upward in terms of temperature, but is it
actually approaching an equilibrium? Because we have set up the system as an
array of each individual point, we can look at a specific point throughout the
evolution of the system. The following three images show the evolution of the
point indicated by the orange dot, along with a logarithmic curve fit of the form
a ∗ xb .

9
 The above image shows the evolution of the orange point for the first time
series of 100 intervals, and goes through almost two full activity cycles. We
can see from the image that the temperature of the point is increasing, and
that the logarithmic fit works very, which indicates that there is no equilibrium
temperature. The parameters for the logarithmic function are:

• a = 3.740132e+05 +/- 2.955532e+04

• b = 1.868421e-01 +/- 2.000152e-02

however, this is a very limited sample to fit a curve to. What happens if we go
further?

10
 Given a few more activity cycles, the fit appears to fit far better. This
plot uses the system that was given 1100 intervals to run through, and there
is a corresponding increase in the number of cycles. Again, however, as the
logarithmic fit seems to fit so wonderfully, there seems to be no equilibrium
point. The parameters for this logarithmic function are:

• a = 7.352039e+04 +/- 4.548315e+03

• b = 5.488182e-01 +/- 9.446914e-03

where, again, we can expand the scope of the data set that can be fit.

11
 This final plot, then, shows the evolution of the system given 2200 inter-
vals. Again, there is a corresponding increase in activity cycles, and, again, the
logarithmic fit seems to be a very close match. The combination of all three
plots showing a logarithmic fit so closely matching each of the different systems
indicates very strongly that there is no equilibrium point for the temperature
of the system. We know that the core will have no equilibrium either because
we know that the outer surfaces are always increasing in temperature, but the
source is not located in the outer surfaces, thus requiring that everywhere in
the system be heated continuously forever. The parameters for this logarithmic
fit are:

• a = 8.627996e+04 +/- 3.888326e+03

• b = 5.242409e-01 +/- 6.233676e-03

where we can see that the uncertainty in the values has decreased from one
system to the next, showing that fit is highly accurate, as the more data you
provide it, the closer it gets to being completely correct (though, given the
oscillatory nature of the system and the non-oscillatory nature of the fit, it will
never be completely correct).
Those three plots showed the evolution for the point represented by the
orange dot, but what about the point represented by the red dot? The behavior

12
of this point turns out to be far different. Rather than steadily increase over
time to become far hotter than the initial temperature, as the orange point did,
the red point experiences a sharp decrease in temperature immediately.

Using an exponential decay fit, scipy’s curvef it was unable to produce values
for the parameters, so, once again, the logarithmic function was used for the
fit. The plot shows this stark contrast in behavior between the two points,
where there is a sharp decrease immediately as the system evolves, but unlike
the previous cases, it is unclear if the system is increasing in temperature, or if
it is remaining constant. The parameters for the fit are:

• a = 7.007377e+05 +/- 9.262445e+04

• b = 1.557041e-06 +/- 3.521911e-02.

Does this point of the system remain constant? The b parameter would cer-
tainly suggest so. Being so near to 0, the fit remains relatively constant. What
happens, however, if we increase the time interval to 1100?

13
 Here there is behavior very similar to what was shown of the previous point.
Thought there is an immediate decrease in temperature, but, just as before,
there is a steady increase of the temperature. Again, as the logarithmic fit seems
to be so closely matched, it appears that this point as well has no equilibrium
point, which would be expected since the system as a whole has no equilibrium.
The parameters for the fit are:

• a = 2.722024e+05 +/- 5.515687e+04

• b = 2.159175e-01 +/- 3.191535e-02.

Unlike before, we can now see from both the plot and the b parameter that the
fit is no longer linear and constant. b is still small, but it is not so small that it
keeps the fit relatively constant. What happens, then, if we increase the time
interval to 2200?

14
 With this final plot, there is a clear and steady increase in the temperature,
though, the temperature of the point has not yet reached its initial state. Given
that there is no equilibrium point, it can be assumed that the system will even-
tually surpass that temperature, but unfortunately, it requires a time interval
far longer than my computer can handle. The parameters for the fit are:

• a = 2.091795e+05 +/- 4.098213e+04

• b = 2.593508e-01 +/- 2.769410e-02.

There are a couple of interesting aspects of the parameters to point out.

The uncertainties are high, but that is not because the fits are inaccurate (at
least, towards the end) but because the initial behavior of the system does not
match any other fit. The behavior of the system towards the end most closely
resembles a logarithmic fit, and so rather than look at the uncertainties for
what they are, they must be looked at comparatively. In each of the systems,
as the time interval is increased, and more activity cycles are run through, the
uncertainties decrease, suggesting that the fit is accurate and adequate. Again,
these three plots suggest that the system as a whole has no equilibrium, but
even though it was set up symmetrically, there are differences in how the system
acts based upon where the system is being observed at.

15
3 Conclusion
When the system was initially set up, there were multiple conditions which were
considered, and a single behavior that was expected to occur. The first condi-
tion was that the system would be made symmetrically, and that the system
would have a symmetrical heat source, meaning that there would be no disparity
between where the heat would be transferred. The second condition was that
the system would be set up as a crude model of the Sun. The layers and their
radii correspond to the Sun’s layers and their radii, and the temperature of each
layer is the same as the expected temperature of the Sun’s layers. The third
condition was that the Sun’s core pulsated in energy output, and that therefore
caused the activity cycles, rather than the reversal of the Sun’s magnetic poles.
The behavior that was expected of the system was that the system would be-
have symmetrically. Whether the surface of the system grew as time progressed,
or whether it oscillated with the core, but remained at a constant range of tem-
peratures, the system was expected to act symmetrically. As the results showed,
however, the system did not act symmetrically, but rather, it had a mirrored
effect. The upper right and lower left portions of the system grew symmetrically
with each other, but asymmetrically with the upper left and lower right portions
of the system, which correspondingly grew symmetrically.
An aspect of the system that was unclear until the function was ran was
whether or not the system would find an equilibrium temperature, or if it would
increase continuously. As was found out, the system was insulated very well,
therefore meaning that because heat was unable to escape the system, it was
contained, and the continually rose, with no equilibrium temperature available.
The feasibility of the oscillatory nature of a core within a star was proven
to be ineffective by this model, as the star would be increasing in temperature
forever. There are, however, a couple of simplifying features that need to be
considered with this conclusion. In reality, stars radiate heat into interstellar
space, through this radiation, they lose energy that was once held within the
star. This model, however, allows for no radiation into space, instead creating
a fully insulated star. Also in reality, stars transfer heat through both radia-
tion and convection, whereas this model simply uses a change in material heat
conductivity to represent that. Furthermore, in reality, a star will have rota-
tion and magnetic field, providing a plethora of different ways in which heat is
transferred and released.
If this model were to be recreated or improved, the key areas to focus on
would be expanding the physical size of the model, and expanding the time
intervals with which the model was ran. My computer, unfortunately, is unable
to handle the amount of information that would be required to run this model to
its full extent, and therefore the results suffered. I believe that the reason that
the model grew mirrored was because the boundaries of the outer layers were
pushed so close together that the convective layer sometimes reached all the
way to the surface, instead of having both the photosphere and chromosphere
buffering it. Also, a magnetic field, the rotation of material, or both could be
introduced to make the system more sophisticated.

16