EX9_POISEUILLE_FLOW.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Viscous flow\n",
    "\n",
    "### Author(s), contact(s), and dates:\n",
    "- Author: Valentina Espinoza Fernández (University of Copenhagen)\n",
    "- Email:  vf@ign.ku.dk\n",
    "- Date:   4th of March 2024"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 9: Asthenospheric flow underneath a plate\n",
    "\n",
    "During the lectures on viscous flow dynamics, one of the exercises (Example 2) showed that the flow down a conduit is the result of gravity working against viscous resistance. In such a case the dynamic balance equation took the form\n",
    "\n",
    "$$ \\mu \\frac{\\delta^2 v_y}{\\delta x^2} + \\rho g = 0$$\n",
    "\n",
    "and the envelope of the resulting flow is a parabola (i.e., a second-degree polynomial function).\n",
    "\n",
    "The case of flow within a **horizontal** conduit, where a constant horizontal pressure gradient ($\\delta p/\\delta x$ = constant) is present, yields a similar equation:\n",
    "\n",
    "$$ \\mu \\frac{\\delta^2 v_x}{\\delta y^2} - \\frac{\\delta p}{\\delta x}  = 0$$\n",
    "\n",
    "and the resulting flow envelope is also a polynomial of second-degree.\n",
    "\n",
    "With reference to Figure 1 (in the exercise paper), imagine a tectonic plate (in orange) of length $L$ that moves at a velocity $v_p$ on top of a Newtonian asthenosphere (in green) of viscosity $\\mu_a$ and thickness $H_a$. \n",
    "\n",
    "Perform the following exercise:\n",
    "\n",
    "1. Obtain an expression for the velocity of asthenospheric horizontal flow underneath the plate ($v_x(y)$) by solving the flow equation (equation 2). **HELP**: under the specific conditions of this exercise, eq. 2 is a relatively simple differential equation. In order to solve it, notice you have the second derivative equal to a term deemed as constant ($\\frac{\\delta p}{\\delta x}$). Think then what are the velocity conditions at the top ($y = H_a$) and bottom ($y = 0$) boundaries of the asthenosphere (i.e., the boundary conditions).\n",
    "\n",
    "2. Imagine also that the horizontal pressure difference within the asthenosphere ($\\delta p$) across the plate length is equal to the pressure generated by a 1-km thick block of crust (density $\\rho$ is 2700 kg/m3). The gradient of pressure along the x-axis can therefore be expressed as:\n",
    "$$-\\rho_c g h / L$$ \n",
    "where $\\rho_c$ is the density of the crust and $h$ is the thickness of the block of crust.\n",
    "Plot a graph of the flow velocity against asthenospheric depth using parameter values (i.e., L, $\\mu_a$, etc.) of your choice.\n",
    "\n",
    "3. Plot the flow velocity against asthenospheric depth for the *Couette* part of the flow, i.e. the flow that is solely induced by plate motion (and not pressure gradients). **HELP**: solve the second derivative equation for the case where $\\delta p = 0$.\n",
    "\n",
    "4. Find the expression for the *Poiseuille* part of the flow, i.e. the flow that is solely induced by the pressure gradient. **HELP**: Isolate the *Couette* part expression from the general solution obtained in (1).\n",
    "\n",
    "5. Obtain an expression for the generic horizontal shear stress within the asthenosphere. **HELP**: recall that the shear stress for a 2D scenario is given by \n",
    "\n",
    "$$\\tau = 2\\mu \\cdot \\frac{1}{2} \\left( \\frac{\\delta v_x}{\\delta y} + \\frac{\\delta v_y}{\\delta x} \\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 2. Horizontal pressure gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To set up a function that calculate the velocity of the asthenosphere at a given depth (y)\n",
    "def velocity(y, v_p, H_a, rho_c, h, mu, L):\n",
    "    g = 9.81        # Gravity acceleration [m/s^2]\n",
    "    vel_m_s = 0.0   # Write an expression for the velocity in [m/s] !!\n",
    "    \n",
    "    return vel_m_s * (3600 * 24 * 365 * 100)    # To return velocity in [cm/yr]\n",
    "\n",
    "\n",
    "# Given values\n",
    "rho_c = 2700    # Density of the crust [kg/m^3]\n",
    "h = 1e3         # Height of the crust [m]\n",
    "\n",
    "\n",
    "# Fill with your own values !!\n",
    "v_p = 0.0   # Velocity of the plate [m/s]\n",
    "H_a = 0.0   # Thickness of the asthenosphere [m]\n",
    "mu = 0.0    # Viscosity of the asthenosphere [Pa s]\n",
    "L = 0.0     # Length of the plate [m]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To set the basic figure, labels and grid\n",
    "fig, ax = plt.subplots() \n",
    "ax.set(xlabel='Velocity (cm/yr)', ylabel='Height (km)')\n",
    "ax.grid(alpha=0.5)\n",
    "\n",
    "y_array = np.linspace(0, H_a, 1000)   # To create an array of 1000 y-values from 0 to H_a\n",
    "x_array = velocity(y_array, v_p, H_a, rho_c, h, mu, L)\n",
    "\n",
    "\n",
    "ax.plot(x_array, y_array/1e3, label='Total flow')       # To plot velocity [cm/yr] vs height [km]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 3: Plot the *Couette* velocity profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def couette_flow(y, v_p, H_a):\n",
    "    vel_m_s = 0.0   # Write an expression for the velocity in [m/s] \n",
    "    return vel_m_s * (3600 * 24 * 365 * 100)    # To return velocity in units cm/yr\n",
    "\n",
    "\n",
    "# To re-use the figure created above (ax), and plot the Couette flow on it\n",
    "ax.plot(couette_flow(y_array, v_p, H_a), y_array/1e3, label='Couette flow')\n",
    "ax.legend()\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 4: Plot the *Poiseuille* velocity profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def poiseuille_flow(y, rho_c, h, mu_a, L, H_a):\n",
    "    g = 9.81        # Gravity acceleration [m/s^2]\n",
    "    vel_m_s = 0.0   # Write an expression for the velocity in [m/s] !!\n",
    "    return vel_m_s * (3600 * 24 * 365 * 100)    # To return velocity in units cm/yr\n",
    "\n",
    "\n",
    "ax.plot(poiseuille_flow(y_array, rho_c, h, mu, L, H_a), y_array/1e3, label='Poiseuille flow')\n",
    "ax.legend()\n",
    "fig"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py38",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}