1.6 Inclined Plane Jupyter Notebook

Instructions

• To use this page, read the text between each of the python code windows, then press RUN to execute the code in the box.
• It is necessary to run each of the boxes in order.
• Don't click restart in the cells unless you come back to the top of the page and start over.
• If you modify the code, this modification will not remain after you leave or refresh this page.
• You must run simulations/cells that require user input to completion. Make sure you see the Message informing you a cell is complete before continuing to the next cell. Here's why you must do this:
  • These cells may not be re-run until the simulation is completed, otherwise no output will be generated.
  • All following cells will not run until the user input cell is completed.
• For Dropdown Menus: You do not need to rerun the cells with the menu to change your selection. However, you must rerun all of the following cells to implement the change.

Be patient, some times it can take 1-2 minutes for the juypter kernel to start.

An interactive example of how to calculate the velocity and strain rate for flow down an inclined plane. For this example, we'll consider the case of magma flow down an inclined plane. The viscosity of magma primarily depends on temperature, and can range from \(10^2\) to \(10^6\) Pa s. For this example we will consider values of \(10^2\) to \(10^4\) Pa s. Hill slopes on the side of volcanoes vary depending on the type of volcano and location relative to the (summit versus flanks). Slopes on shield volcanos are typically 5 degrees, while the slopes on stratovolcanos are 5 to 30 degrees. Very close to the vent at the peak of a stratovolcano, the slopes can locally reach up to 70-80 degrees.

Key Questions: Consider these as you work your way through this page.

• What conditions (angle and viscosity) produce the highest flow velocities?
• How do velocity and strain rate change within the flow with depth?
• Where, in the flow profile, do the highest magnitude strain rates occur?
• Where, in the flow profile, do the highest magnitude flow velocities occur?

First, we import the necessary packages and set physical constants. Note: When you first run this box upon loading or refreshing the page it will say "Waiting for kernel...". This is supposed to happen and may take a minute or so.

import numpy as np
import matplotlib.pyplot as plt

g = 9.81 #m/s

print('Done. Go to next box.')

Now, let's set some parameters and define functions to calculate the desired quantities, this way it's easier if we want to use the formula multiple times. This function takes 5 input parameters: alpha, rho, eta, y, and h. Alpha is the angle of the inclined plane, rho is the density of the flowing material, eta is the viscosity of the material, and h is the thickness of the flow. y is a range of depths within the flow, starting at 0, and ending at h with an increment defined by step. The density used here is the density of basalt. Then, the function calculates and returns the velocity and strain rate versus depth. The calculations are done using the equations derived in 1.5 Viscous Deformation.

h = 50 #m
rho = 3000 # kg / m^3
step = 1
y = np.arange(0,h+step,step) #m

def inclined_plane_velocity(alpha,rho,eta,h,y): #alpha in degrees, rho in kg/m^3, eta in Pa*s, and y and h in m
    v = (rho * g * np.sin(np.radians(alpha)) / (2*eta)) * (h**2 - y**2) # m/s velocity
    return v

def inclined_plane_strain_rate(alpha,rho,eta,y): #rho in kg/m^3, alpha in degrees,  eta in Pa*s, and y and h in m
    return -(rho * g * np.sin(np.radians(alpha)) * y / (2*eta)) # (1/s)

print('The parameters and functions are defined, move on to the next cell')

Throughout this entire exercise, the flow thickness and density will remain constant, and we will consider different angles and viscosities. We need to define the angle and viscosity for this cell, let's vary the angle and hold viscosity constant (eta). This cell uses the functions above to calculate and plot the velocity and strain rate versus depth for each angle in alpha_range. Run the simulation and consider the key questions. You may enter any angle you like between 0 and 80.

eta = 1e2 # Pa*s
handle = display(None, display_id=True)

x = np.zeros(len(y)) 

bbox = dict(boxstyle ="round", fc = '1') 
fig, (ax,ax1) = plt.subplots(1,2,figsize=(14,5))
line, = ax.plot(x,y,color='red',label = 'Velocity')
line1, = ax1.plot(x,y,color='blue',label = 'Strain Rate')
ax.set_ylabel('Depth (m)')
ax.set_xlabel('Velocity (km/s)')
ax.set_xlim(0,350)
ax.set_ylim(0,h+step)
ax.invert_yaxis()
ax.invert_xaxis()
ax.set_title('Inclined Plane Flow Velocity versus Depth: Angle Variation')
ax.legend(loc=3)
ax1.set_ylabel('Depth (m)')
ax1.set_xlabel('Strain Rate (1/s)')
ax1.set_title('Inclined Plane Flow Strain Rate versus Depth: Angle Variation')
ax1.set_ylim(0,h+step)
ax1.invert_yaxis()
ax1.set_xlim(-7000,0)

ax1.legend(loc=6)
print('Enter a angle for the inclined plane between 0 and 80 degrees: ')
alpha_current = float(input())
while((alpha_current<=80) & (alpha_current>=0)):
    print('Now Plotting:', alpha_current, 'degrees')
    v_current = inclined_plane_velocity(alpha_current,rho,eta,h,y)/1000 # km / s
    sr_current = inclined_plane_strain_rate(alpha_current,rho,eta,y) # 1/s
    line.set_xdata(v_current)
    line1.set_xdata(sr_current)
    
    annotation = ax.annotate('Angle = ' + str(alpha_current) + '$ ^{\circ} $',(300,35),bbox=bbox,color = 'red',fontsize = 12)
    annotation1 = ax1.annotate('Angle = ' + str(alpha_current) + '$ ^{\circ} $',(-6800,20),bbox=bbox,color = 'blue',fontsize = 12)
    handle.update(fig)
    print('Enter a number greater than 80 to stop plotting or Enter an angle (0-80 deg) to continue: ')
    alpha_current = float(input())
    annotation.remove()
    annotation1.remove()
plt.close()

print('The simulation is complete. Consider the key questions related to these figures, then move on.')

Now let's hold the angle of the plane constant and vary viscosity. The first lines in the cell take the angle to be 15 degrees, and create a range of viscosities (eta_range). The rest of the cell is set up exactly like the previous one, the only difference is that now the for-loop is iterating over the viscosities and not angles. Run the simulation and observe the results. Press Enter or Return to change the viscosity in the plot below after you run the cell.

alpha = 15 # degrees
eta_range = np.arange(1e2,1.05e4,(1.05e4-1e2)/5) #Pa*s

handle = display(None, display_id=True)

x = np.zeros(len(y)) 

bbox = dict(boxstyle ="round", fc = '1') 
fig, (ax,ax1) = plt.subplots(1,2,figsize=(15,5))
line, = ax.plot(x,y,color='red',label = 'Velocity')
line1, = ax1.plot(x,y,color='blue',label = 'Strain Rate')
ax.set_ylabel('Depth (m)')
ax.set_xlabel('Velocity (km/s)')
ax.set_xlim(0,120)
ax.set_ylim(0,h+step)
ax.invert_yaxis()
ax.invert_xaxis()
ax.set_title('Inclined Plane Flow Velocity versus Depth: Viscosity Variation')
ax.legend(loc=3)
ax1.set_ylabel('Depth (m)')
ax1.set_xlabel('Strain Rate (1/s)')
ax1.set_title('Inclined Plane Flow Strain Rate versus Depth: Viscosity Variation')
ax1.set_ylim(0,h+step)
ax1.invert_yaxis()
ax1.set_xlim(-2800,0)

ax1.legend(loc=6)
for eta_current in eta_range:
    print('Now Plotting:', eta_current, 'Pa*s')
    v_current = inclined_plane_velocity(alpha,rho,eta_current,h,y)/1000 # km / s
    sr_current = inclined_plane_strain_rate(alpha,rho,eta_current,y) # 1/s
    line.set_xdata(v_current)
    line1.set_xdata(sr_current)
    
    annotation = ax.annotate('Viscosity = ' + '{:.2E}'.format(eta_current) + ' Pa*s',(110,35),bbox=bbox,color = 'red',fontsize = 12)
    annotation1 = ax1.annotate('Viscosity = ' + '{:.2E}'.format(eta_current) + ' Pa*s',(-2750,20),bbox=bbox,color = 'blue',fontsize = 12)
    handle.update(fig)
    if eta_current == eta_range[0]:
    	print('Click enter or return to continue: ')
    wait = input()
    annotation.remove()
    annotation1.remove()
plt.close()

print('The simulation is complete. Consider the key questions related to these figures, then move on.')

Now let's create a physical model of the flow by rotating the previous velocity graph such that the depth axis is at the angle of the inclined plane. For this example, we will look at three different angles (alpha_range_2) and will hold viscosity constant (eta_2). We will use a rotation matrix to rotate the velocity profile. This involves some basic linear algebra, so you will not be required to come up with something like this on your own for this course, however, pay attention to the output. This cell defines some functions to accomplish this, calls them, and graphs the result. This cell creates normal graphs, not simulations like the above cells.

eta_2 = 1e6
alpha_range_2 = np.array([10,40,70]) #degrees
def rotation(alpha,x ,y):
    alpha_rad = np.radians(alpha)
    r = np.matrix([[np.cos(alpha_rad), -np.sin(alpha_rad)], [np.sin(alpha_rad), np.cos(alpha_rad)]])
    Mesh = np.matrix([x,y])
    Rot = r * Mesh
    x_rot = np.array(Rot[0])[0]
    y_rot = np.array(Rot[1])[0]
    return x_rot, y_rot

def physical_model_ip(alpha,rho,eta,h):
    step = 1
    y = np.arange(0,h+step,step)
    Velocities = (rho * g * np.sin(np.radians(alpha)) / (2*eta)) * (h**2 - y**2) # m/s velocity
    vel_rotated, y_rotated = rotation(alpha,Velocities,y)
    q = np.arange(-1,2*np.max(Velocities)) #X-Coordinates for horizontal Lines
    top = np.zeros(len(q)) #Top line
    bottom = np.ones(len(q))*h #Bottom Line
    q_rot_1, top_rot = rotation(alpha,q,top)
    q_rot_2, bot_rot = rotation(alpha,q,bottom)
    
    # Calculate the aspect ratio of the rotated box so you can rest this to 
    # to match an aspect ratio of 1 for the unrotated box.
    asp2 = (np.max(q_rot_1)-np.min(q_rot_2))/(np.max(bot_rot)-np.min(top_rot))
    return vel_rotated, y_rotated, q_rot_1, q_rot_2, top_rot, bot_rot, asp2
    
vel_rotated, y_rotated, q_rot_1, q_rot_2, top_rot, bot_rot, asp2 = physical_model_ip(alpha_range_2[0],rho,eta_2,h)
    
fig, (ax,ax1,ax2) = plt.subplots(1,3,figsize = (15,5))


ax.plot(vel_rotated,y_rotated,color='blue',label = 'velocity')
ax.plot(q_rot_1,top_rot,color='black',label = 'Top of Flow')
ax.plot(q_rot_2,bot_rot,color='gray',label = 'Bottom of Flow')
ax.invert_xaxis()
ax.invert_yaxis()
ax.axis('off')
ax.legend(loc=9, bbox_to_anchor=(1, 1))
ax.set_aspect(asp2)
ax.set_title('Shape of the flow profile for alpha = ' + str(alpha_range_2[0])+'$ ^{\circ} $')
ax.annotate('$u_{max}$',(np.max(vel_rotated),y_rotated[np.where(vel_rotated == np.max(vel_rotated))[0]]),xytext = (np.max(vel_rotated)+3,y_rotated[np.where(vel_rotated == np.max(vel_rotated))[0]] + 15),color = 'red',size = 15,arrowprops = {'arrowstyle':'->','color':'red'})

vel_rotated, y_rotated, q_rot_1, q_rot_2, top_rot, bot_rot, asp2 = physical_model_ip(alpha_range_2[1],rho,eta_2,h)
ax1.plot(vel_rotated,y_rotated,color='blue',label = 'velocity')
ax1.plot(q_rot_1,top_rot,color='black',label = 'Top of Flow')
ax1.plot(q_rot_2,bot_rot,color='gray',label = 'Bottom of Flow')
ax1.invert_xaxis()
ax1.invert_yaxis()
ax1.axis('off')
ax1.legend(loc=9, bbox_to_anchor=(1, 1))
ax1.set_aspect(asp2)
ax1.set_title('Shape of the flow profile for alpha = ' + str(alpha_range_2[1])+'$ ^{\circ} $')
ax1.annotate('$u_{max}$',(np.max(vel_rotated),y_rotated[np.where(vel_rotated == np.max(vel_rotated))[0]]),xytext = (np.max(vel_rotated)+5,y_rotated[np.where(vel_rotated == np.max(vel_rotated))[0]] + 15),color = 'red',size = 15,arrowprops = {'arrowstyle':'->','color':'red'})

vel_rotated, y_rotated, q_rot_1, q_rot_2, top_rot, bot_rot, asp2 = physical_model_ip(alpha_range_2[2],rho,eta_2,h)
ax2.plot(vel_rotated,y_rotated,color='blue',label = 'velocity')
ax2.plot(q_rot_1,top_rot,color='black',label = 'Top of Flow')
ax2.plot(q_rot_2,bot_rot,color='gray',label = 'Bottom of Flow')
ax2.invert_xaxis()
ax2.invert_yaxis()
ax2.axis('off')
ax2.legend(loc=9, bbox_to_anchor=(1, 1))
ax2.set_aspect(asp2)
ax2.set_title('Shape of the flow profile for alpha = ' + str(alpha_range_2[2])+'$ ^{\circ} $')
ax2.annotate('$u_{max}$',(np.max(vel_rotated),y_rotated[np.where(vel_rotated == np.max(vel_rotated))[0]]),xytext = (np.max(vel_rotated)+5,y_rotated[np.where(vel_rotated == np.max(vel_rotated))[0]] + 15),color = 'red',size = 15,arrowprops = {'arrowstyle':'->','color':'red'})



plt.tight_layout()
plt.show()
    
print('The figures are plotted, consider the key questions.')

Using the figures above, how would you answer the key questions:

• What conditions (angle and viscosity) produce the highest flow velocities?
• How do velocity and strain rate change within the flow with depth?
• Where, in the flow profile, do the highest magnitude strain rates occur?
• Where, in the flow profile, do the highest magnitude flow velocities occur?