Plot g(r) of planet with two distinct density regions

1. Oct 31, 2013

oddjobmj

1. The problem statement, all variables and given/known data

Consider a simple model for the interior of the Earth: there is a spherical iron core with constant mass density ρ0 and radius a; outside the core is "rock" with constant density ρ1. Use these values for the densities: ρ0= 8.70×103 kg/m3 and ρ1= 4.10×103 kg/m3. The radius of the Earth is R = 6.40×106 m.

Hand in a plot g(r) from 0 to 2R.
The value at r=R must be 9.81 m/s2.

2. Relevant equations

g(r)=$\frac{G*M(r)}{r^2}$

a=4.291*106 m

Total Mass=$\frac{4πρ0}{3}$a^3+$\frac{4πρ}{3}$(R3-a3)

3. The attempt at a solution

Since the mass is a function of r the function will be affected by the density at the current radius. Is there a way to write a single function that will cover this or do I have to do it with a piecewise function where I use ρ0 between 0 < r < a, then ρ between a < r < R, then constant past R?

2. Oct 31, 2013

cepheid

Staff Emeritus
rho(r) is by definition piecewise in this problem.

You just need to find the function M(r) in terms of rho(r). M(r) is the mass enclosed by radius r.

Hopefully when you say "constant past R", you mean constant and zero past R.

3. Oct 31, 2013

oddjobmj

Thank you, cepheid.

When I said constant past r I meant M(r) is constant past r. Yes, the density would be zero.

4. Oct 31, 2013

oddjobmj

Actually, the function for M(r) would be zero if the density were zero. That would lead to 0/r^2 which would basically mean my g(r) wouldn't change at all beyond that point. I don't believe that should be the case. So, at that point I need to simply calculate the total mass of the Earth and plug that in instead of the function for M(r).

Edit:
Can you tell what I did wrong here?

http://imgur.com/0Bp1aaM

Edit 2:

I think the second function should have (r3-a3) instead of (r-a)^3. All the lines are connected now and it looks as if it exhibits the characteristics one would expect from the system.

Last edited: Oct 31, 2013
5. Oct 31, 2013

cepheid

Staff Emeritus
No it wouldn't. Remember that M(r) is the total mass enclosed by the radius r. This is non-zero for r > R. It's just constant, as you pointed out. Hint: for r > a, both densities come into play in calculating M enclosed, not just one.

6. Oct 31, 2013

oddjobmj

I'm sorry, I was unclear. You are correct. I didn't mean to imply that M(r) should actually equal zero. I just meant that if I happened to plug in zero for the density it would spit out zero. Instead I just set the mass at r>R to the total mass of the planet.

7. Oct 31, 2013

cepheid

Staff Emeritus
Cool. Just out of curiosity, what does your new plot look like?

8. Oct 31, 2013

oddjobmj

9. Oct 31, 2013

cepheid

Staff Emeritus
Cool. I got some very similar looking output (attached png).

Code (Text):

#
# g_earth.py -- simple model of earth with const. density iron core for
#               r < a, and const  density rock mantle for a < r < R.
#

import numpy as np
import matplotlib.pyplot as plt
import math

R = 6.40e6      # metres
rho0 = 8.70e3   # kg/m^3
rho1 = 4.10e3   # kg/m^3
a = 4.291e6     # metres
G = 6.67384e-11 # m^3 kg^-1 s^-1

r = np.arange(1.0, 1.0e7, 1.0e3)

M = np.zeros(len(r))
g = np.zeros(len(r))

for i in xrange(len(r)):
if (r[i] < a):
M[i] = rho0*(4.0/3.0)*math.pi*(r[i])**3
g[i] = (G*M[i])/(r[i])**2
elif (r[i] < R):
M[i] = rho1*(4.0/3.0)*math.pi*((r[i])**3 - a**3)\
+ rho0*(4.0/3.0)*math.pi*(a**3)
g[i] = (G*M[i])/(r[i])**2
else:
M[i] = rho1*(4.0/3.0)*math.pi*(R**3 - a**3)\
+ rho0*(4.0/3.0)*math.pi*(a**3)
g[i] = (G*M[i])/(r[i])**2

plt.figure(1, figsize=(12,8), dpi=100)
plt.plot(r,g)
plt.title("Gravitational Field Strength as a Function of Radial Distance for Simple Earth Model")
plt.ylabel("g(r) [N/kg]")
plt.xlabel("r [m]")
plt.savefig("g_earth.png")
plt.show()

Attached Files:

• g_earth.png
File size:
12.7 KB
Views:
59
10. Oct 31, 2013

oddjobmj

Ah, I didn't even think about breaking out IDLE... Also, I'll sleep better tonight knowing you came up with similar results. Thank you!

11. Oct 31, 2013

cepheid

Staff Emeritus
Since g ~ M/r^2, we can predict what g(r) should look like in the three regimes.

For r < a, M(r) ~ r^3, therefore g ~ r^3/r^2. In other words, g ~ r so the linear relationship is totally expected. Here I use ~ to mean "scales as" or "is proportional to."

For r > R, M(r) = const., therefore g ~ 1/r^2, which looks about right on the plot

The middle phase is tricker.