1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Oct 31, 2013 #1
    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)=[itex]\frac{G*M(r)}{r^2}[/itex]

    a=4.291*106 m

    Total Mass=[itex]\frac{4πρ0}{3}[/itex]a^3+[itex]\frac{4πρ}{3}[/itex](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. jcsd
  3. Oct 31, 2013 #2

    cepheid

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    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.
     
  4. Oct 31, 2013 #3
    Thank you, cepheid.

    When I said constant past r I meant M(r) is constant past r. Yes, the density would be zero.
     
  5. Oct 31, 2013 #4
    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
  6. Oct 31, 2013 #5

    cepheid

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    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.
     
  7. Oct 31, 2013 #6
    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.
     
  8. Oct 31, 2013 #7

    cepheid

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    Cool. Just out of curiosity, what does your new plot look like?
     
  9. Oct 31, 2013 #8
  10. Oct 31, 2013 #9

    cepheid

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    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:

  11. Oct 31, 2013 #10
    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!
     
  12. Oct 31, 2013 #11

    cepheid

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Plot g(r) of planet with two distinct density regions
  1. Planet density (Replies: 2)

  2. Planet density (Replies: 2)

  3. Density of a planet (Replies: 6)

  4. Plot the range r (Replies: 11)

Loading...