{ "metadata": { "name": "", "signature": "sha256:4b88382b5a09a2b6942dfcdeccdf53bbf1eba6cc8151c227802212830ea76e16" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Acoustics equations and their Riemann solver in 1D and 2D\n", "by Mauricio J. Del Razo S. 2014\n", "\n", "A brief intuitive derivation of the acoustics equations is presented. We derive the normal Riemann solver for one and two dimensional acoustics and the transverse Riemann solver for two dimensional acoustics in a cartesian grid. These solvers are further extended for general quadrilateral mapped grids. Along some of these examples, we will show how to use PyClaw to solve the acoustics equations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Linear Acoustics equations\n", "The linear acoustic equations can be obtained by linearizing the conservation of mass and momentum for an element of fluid. The conservation of mass and momentum are given by\n", "\n", "$$\n", "\\frac{\\partial \\rho}{\\partial t} + \\nabla \\cdot \\left(\\rho \\bar{u}\\right) = 0 \\\\\n", "\\frac{d \\rho \\bar{u}}{dt} = -\\nabla P\n", "$$\n", "\n", "where $d/dt$ denotes the material derivative: $\\frac{d}{dt} = \\frac{\\partial}{\\partial t} + \\bar{u} \\cdot \\nabla $. In one dimension, these can be easily rewritten as,\n", "\n", "$$\n", "\\rho_t + (\\rho u)_x = 0, \\\\\n", "(\\rho u)_t + (\\rho u^2 + P(\\rho))_x =0,\n", "$$\n", "\n", "where we assumed the pressure is a function of the density. Linearizing the equation around $\\rho_0$ and $u_0 =0$, we obtain the simpler system for the perturbations,\n", "\n", "$$\n", "\\rho_t + (\\rho u)_x = 0, \\\\\n", "(\\rho u)_t + P(\\rho_0)\\rho_x =0.\n", "$$\n", "\n", "where $\\rho$ and $u$ are now the perturbations around $\\rho_0$ and $u_0 = 0$. As perturbations on the pressure and density satisfy, \n", "\n", "$$\n", "p\\approx P'(\\rho_0)\\rho, \\\\\n", "\\rho u \\approx \\rho_0 u,\n", "$$\n", "\n", "we can rewrite the system of linear acoustic equations as\n", "\n", "$$\n", " \\left[ \\begin{array}{c}\n", "p \\\\\n", "u \n", "\\end{array} \\right]_t\n", "+ \\underbrace{\\left[ \\begin{array}{cc}\n", "0 & K_0 \\\\\n", "1/\\rho_0 & 0 \\\\\n", "\\end{array} \\right]}_{\\mathbf{A}}\n", "\\left[ \\begin{array}{c}\n", "p \\\\\n", "u \\end{array} \\right]_x = 0,\n", "$$\n", "\n", "with $K_0=\\rho_0 P'(\\rho_0)$ the bulk modulus of compressibility that will tell us how compressible is a material. This can be easily extended to two dimensions,\n", "\n", "$$\n", " \\left[ \\begin{array}{c}\n", "p \\\\\n", "u \\\\\n", "v\n", "\\end{array} \\right]_t\n", "+ \\underbrace{\\left[ \\begin{array}{ccc}\n", "0 & K_0 & 0 \\\\\n", "1/\\rho_0 & 0 & 0 \\\\\n", "0 & 0 & 0 \\\\\n", "\\end{array} \\right]}_{\\mathbf{A}}\n", "\\left[ \\begin{array}{c}\n", "p \\\\\n", "u \\\\\n", "v \\end{array} \\right]_x\n", "+\\underbrace{\\left[ \\begin{array}{ccc}\n", "0 & 0 & K_0 \\\\\n", "0 & 0 & 0 \\\\\n", "1/\\rho_0 & 0 & 0 \\\\\n", "\\end{array} \\right]}_{\\mathbf{B}}\n", "\\left[ \\begin{array}{c}\n", "p \\\\\n", "u \\\\\n", "v \\end{array} \\right]_y = 0\n", "$$\n", "\n", "and three dimensions\n", "\n", "\n", "$$\n", " \\left[ \\begin{array}{c}\n", "p \\\\\n", "u \\\\\n", "v \\\\\n", "w\n", "\\end{array} \\right]_t\n", "+ \\underbrace{\\left[ \\begin{array}{ccc}\n", "0 & K_0 & 0 & 0 \\\\\n", "1/\\rho_0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0\n", "\\end{array} \\right]}_{\\mathbf{A}}\n", "\\left[ \\begin{array}{c}\n", "p \\\\\n", "u \\\\\n", "v \\\\\n", "w\\end{array} \\right]_x\n", "+\\underbrace{\\left[ \\begin{array}{ccc}\n", "0 & 0 & K_0 & 0 \\\\\n", "0 & 0 & 0 & 0 \\\\\n", "1/\\rho_0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0\\\\\n", "\\end{array} \\right]}_{\\mathbf{B}}\n", "\\left[ \\begin{array}{c}\n", "p \\\\\n", "u \\\\\n", "v \\\\\n", "w\\end{array} \\right]_y \n", "+\\underbrace{\\left[ \\begin{array}{ccc}\n", "0 & 0 & 0 & K_0 \\\\\n", "0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 \\\\\n", "1/\\rho_0 & 0 & 0 & 0\\\\\n", "\\end{array} \\right]}_{\\mathbf{C}}\n", "\\left[ \\begin{array}{c}\n", "p \\\\\n", "u \\\\\n", "v \\\\\n", "w\\end{array} \\right]_z= 0\n", "$$\n", "where the velocity in $x,y,z$ is $u,v,w$ and the matrices $A$, $B$ and $C$ have those forms depending on the dimensionality." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. The Riemann problem in 1D (extended to discontinuous density and bulk modulus)\n", "In one dimension we need to solve,\n", "\n", "$$\n", "\\bar{q}_t + \\mathbf{A} \\bar{q}_x = 0, \\\\\n", "$$\n", "with\n", "$$\n", "\\bar{q} = \\left[ \\begin{array}{c}\n", "p \\\\\n", "u\n", "\\end{array}\\right]\n", "\\hspace{20mm}\n", "\\mathbf{A} = \n", "\\left[ \\begin{array}{cc}\n", "0 & K_0 \\\\\n", "1/\\rho_0 & 0\n", "\\end{array} \\right],\n", "$$\n", "and initial condition\n", "$$\n", "\\bar{q}(x,0) = \\begin{cases}\n", "\\bar{q_L} & \\text{if } x \\le 0, \\\\\n", "\\bar{q_R} & \\text{if } x > 0.\n", "\\end{cases}\n", "$$\n", "The solution of the Riemann problem can be obtained by transforming the system into two uncoupled advection equations, see [LeVeque 2002](http://depts.washington.edu/clawpack/book.html). When transformed back to the original coordinates the solution will have the following structure, \n", "\n", "\n", "\n", "Two propagating acoustic waves (one to the left and one to the right) with speeds $c_L$ and $c_R$, the two initial conditions $q_L$ and $q_R$ in the far field and one more extra state in between the waves $q_M$. If we can find $c_L$, $c_R$ and $q_M$, we have solved the problem.\n", "\n", "The following section show how to obtain this solution to the Riemann problem, for details on the methodology, see [LeVeque 2002](http://depts.washington.edu/clawpack/book.html). We will begin by obtaining the eigenvector and eigenvalues of $\\mathbf{A}$ symbolically (although this might be very simple to do by hand, on more complicated systems of equations might become handy)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Initialize sympy for symbolic computations\n", "import sympy as sy\n", "# Mathjax printing with simpy\n", "sy.init_printing(use_latex='mathjax')\n", "\n", "# Declare symbolic variables\n", "rho_0, K_0 = sy.symbols('rho_0 K_0', real=True)\n", "R = sy.Matrix(2, 2, lambda i,j: 0)\n", "Rval = sy.Matrix(2, 1, lambda i,j: 0)\n", "\n", "# Define symbolic matrix\n", "An = sy.Matrix([[0, K_0],\n", " [1/rho_0, 0]])\n", "\n", "# Calculate eigenvalues and eigenvectors of matrix An\n", "# eigvec[i][0] => ith eigenvalue\n", "# eigvec[i][2] => ith eigenvector\n", "eig_vec = An.eigenvects()\n", "\n", "\n", "# Simplify eigenvalues and eigenvectors and save in symbolic arrays\n", "for i in range(2):\n", " Rval[i] = sy.simplify(eig_vec[i][0])\n", " for j in range(2):\n", " R[j,i] = sy.simplify(eig_vec[i][2][0][j])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "# Show eigenvalue and corresponding eigenvectors in a widget\n", "\n", "from IPython.html.widgets import interact, interactive\n", "from IPython.display import display, HTML, Latex\n", "\n", "# Function to print the eigenvalue and eigenvector\n", "def print_eigen(Eigenvalue):\n", " display(Latex('Eigenvalue =' + '$' + sy.latex(Rval[Eigenvalue]) + '\\hspace{20mm} $' + \n", " 'Eigenvector =$' + sy.latex(R[:,Eigenvalue]) + '$'))\n", " \n", "# Interactive widget to display eigenvalue and eigenvector\n", "interact(print_eigen, Eigenvalue= {'one':0,'two':1});" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "Eigenvalue =$- \\sqrt{\\frac{K_{0}}{\\rho_{0}}}\\hspace{20mm} $Eigenvector =$\\left[\\begin{smallmatrix}- \\frac{K_{0}}{\\sqrt{\\frac{K_{0}}{\\rho_{0}}}}\\\\1\\end{smallmatrix}\\right]$" ], "metadata": {}, "output_type": "display_data", "text": [ "" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The eigenvalues and eigenvectors obtained are the same than those in [LeVeque 2002](http://depts.washington.edu/clawpack/book.html), except for constant factors in the eigenvectors. The speeds of the left and right going acoustic waves are given by the eigenvalues: $s_L = -c_0$, $s_R = c_0$ with: $c_0 = \\sqrt{K_0/\\rho_0}$.\n", "However, we would like to have different materials in the left and right side, which means different density ($\\rho_L$ and $\\rho_R$) and bulk modulus ($K_L$ and $K_R$), therefore\n", "$$\n", "s_L = -c_L \\hspace{20mm} s_R = c_R \\hspace{20mm} \\mathrm{with:} \\hspace{5mm} c_M = \\sqrt{\\frac{K_{M}}{\\rho_{M}}},\n", "$$\n", "with $M=L,R$. Using the fact that $K_{M}=c_{M}^2 \\rho_{M}$ and recalling the specific acoustic impedance is given by $Z_{M}=\\rho_{M} c_{M} = K_M/\\sqrt{K_M/\\rho_M}$, we can write the matrix of column eigenvectors $\\mathbf{R}=[\\bar{r_L}, \\bar{r_R}]$ as,\n", "$$\n", "\\mathbf{R} = \n", "\\left[ \\begin{array}{ccccc}\n", "-Z_{L} & Z_{R} \\\\\n", " 1 & 1 \\\\\n", "\\end{array} \\right],\n", "$$\n", "\n", "Let $\\delta \\bar{Q} = [\\Delta_p, \\Delta_u]^T$ be the jump of $\\bar{q}$ across the discontinuity. In order to solve the Riemann problem, we need to expand the jump accross the discontinuity $\\Delta \\bar{q}$ as a linear combination of the eigenvectors,\n", "$$\n", "\\alpha_L \\bar{r_L} + \\alpha_R \\bar{r_R} = \\Delta \\bar{q},\n", "$$\n", "or more compactly\n", "$$\n", "\\mathbf{R} \\bar{\\alpha} = \\Delta \\bar{q}.\n", "$$\n", "We will need to obtain the value of $\\bar{\\alpha}$, solving again this equation symbolically for illustrative purposes," ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Declare symbolic variables\n", "ZL, ZR, dp, du = sy.symbols('Z_L Z_R Delta_p Delta_u', real=True)\n", "R = sy.Matrix(2, 2, lambda i,j: 0)\n", "\n", "# Define symbolic matrix and vector\n", "R = sy.Matrix([[-ZL, ZR],\n", " [1, 1]])\n", "dq = sy.Matrix([[dp],[du]])\n", "\n", "# Solve system\n", "alpha = R.inv()*dq\n", "\n", "#Simplify and print\n", "al = sy.simplify(alpha[0])\n", "ar = sy.simplify(alpha[1])\n", "display(Latex('$\\\\alpha_L =$' + '$' + sy.latex(al) + '$' ))\n", "display(Latex('$\\\\alpha_R =$' + '$' + sy.latex(ar) + '$' ))\n" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$\\alpha_L =$$\\frac{- \\Delta_{p} + \\Delta_{u} Z_{R}}{Z_{L} + Z_{R}}$" ], "metadata": {}, "output_type": "display_data", "text": [ "" ] }, { "latex": [ "$\\alpha_R =$$\\frac{\\Delta_{p} + \\Delta_{u} Z_{L}}{Z_{L} + Z_{R}}$" ], "metadata": {}, "output_type": "display_data", "text": [ "" ] } ], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We already know the wave speeds $c_L$ and $c_R$ and the initial states $q_L$ and $q_R$. We only need to know the middle state $q_M$, which will be given by\n", "$$\n", "\\bar{q_M} = \\bar{q_L} + \\alpha_L \\bar{r_L} = \\bar{q_R} - \\alpha_R \\bar{r_R}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Implementing a sample solution into python\n", "Now we can plot and construct the solution to the Riemann problem. We will give an outline of the algorithm followed by the code:\n", "
    \n", "
  1. Define piecewise constant initial condition with the left and right states and corresponding bulk modulii,\n", "$$\n", "q_L = [p_L, u_L] \\ \\ \\text{with} \\ \\ K_L, \\\\\n", "q_R = [p_R, u_R] \\ \\ \\text{with} \\ \\ K_R. \n", "$$\n", "
    2. \n", "
  2. Calculate the wave speeds given by the eigenvalues and the corresponding eigenvectors of $\\mathbf{A}$.\n", "
  3. Compute $\\bar{\\alpha}$ that are the weights of the jumps in the direction of each of the eigenvectors
    4. \n", "
  4. Compute the middle state $q_{m}$ using the following formulas for the middle states,\n", "\n", "$$\n", "q_{m} = q_{L} + \\alpha_L \\bar{r_L} \\ \\ \\ \\mathrm{or} \\\\\n", "q_{m} = q_{R} - \\alpha_R \\bar{r_R}\n", "$$\n", "
    5. \n", "Now, we will implement this into a python script and plot the solution." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# SOLUTION IMPLEMENTATION\n", "#import numpy\n", "import numpy as np\n", "\n", "# Set left and right state and deltaq (q = [sig11,sig22,sig12,u,v])\n", "qL = np.array([-1.0, 0.0])\n", "qR = np.array([1.0, 0.0])\n", "dq = np.array([qR[0]-qL[0], qR[1]-qL[1]])\n", "\n", "# Set bulk and density (left and right)\n", "bulkL = 1; rhoL = 1 \n", "bulkR = 4; rhoR = 2\n", "\n", "# Define speeds and impedance (left and right)\n", "cL = np.sqrt(bulkL/rhoL); ZL = rhoL*cL\n", "cR = np.sqrt(bulkR/rhoR); ZR = rhoR*cR\n", "\n", "# Define the 2 eigenvectors (from columns of Matrix R)\n", "r1 = np.array([-ZL, 1])\n", "r2 = np.array([ZR, 1])\n", "\n", "# Compute the 2 alphas\n", "det = ZL + ZR\n", "alL = (-dq[0] + dq[1]*ZR)/det\n", "alR = (dq[0] + dq[1]*ZL)/det\n", "\n", "# Compute middle state qm\n", "qm = qL + alL*r1 \n", "## Should be equivalent to\n", "#qms = qR - alR*r2 #it is!\n", " \n", "# Compute waves characteristics for plotting\n", "x = np.linspace(-5,5,50)\n", "Wc = np.zeros((2,len(x)))\n", "Wc[0][:] = -x/cL\n", "Wc[1][:] = x/cR" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "#SOLUTION PLOTTING\n", "# Required for plotting\n", "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pylab import *\n", "from IPython.html.widgets import interact\n", "from ipywidgets import StaticInteract, RangeWidget, RadioWidget, DropDownWidget\n", "\n", "# Plot Riemann solution\n", "def plot_Riemann(time):\n", " \n", " fig = plt.figure(figsize=(15, 5))\n", " \n", " # Create subplot for (x,t) plane\n", " tmax = 5\n", " ax = fig.add_subplot(1,3,1)\n", " ax.set_xlabel('x')\n", " ax.set_ylabel('time')\n", " ax.axis([min(x),max(x),0,tmax])\n", " \n", " # Plot characteristic lines\n", " # Acoustic waves in red\n", " ax.plot(x,Wc[0][:], '-r', linewidth=2)\n", " ax.plot(x,Wc[1][:], '-r', linewidth=2)\n", " # Plot time-line in (x,t) plane\n", " ax.plot(x, 0*x+time, 'k', linewidth=3)\n", "\n", " # Create pressure subplot for Riemann solution\n", " ax2 = fig.add_subplot(1,3,2)\n", " ax2.set_xlabel('x')\n", " ax2.set_ylabel('Pressure')\n", " ax2.axis([min(x),max(x),-2,2])\n", " \n", " # Create Riemann solution vector and plot\n", " sol = np.zeros(len(x))\n", " for i in range(len(x)):\n", " if x[i] < -cL*time:\n", " sol[i] = qL[0]\n", " elif x[i] < cR*time:\n", " sol[i] = qm[0]\n", " else:\n", " sol[i] = qR[0]\n", " ax2.plot(x,sol, 'k', linewidth = 3)\n", " ax2.fill_between(x,-20, sol, facecolor='blue', alpha=0.2)\n", " \n", " # Create velocity subplot for Riemann solution\n", " ax3 = fig.add_subplot(1,3,3)\n", " ax3.set_xlabel('x')\n", " ax3.set_ylabel('Velocity')\n", " ax3.axis([min(x),max(x),-2,2])\n", " \n", " # Create Riemann solution vector and plot\n", " sol = np.zeros(len(x))\n", " for i in range(len(x)):\n", " if x[i] < -cL*time:\n", " sol[i] = qL[1]\n", " elif x[i] < cR*time:\n", " sol[i] = qm[1]\n", " else:\n", " sol[i] = qR[1]\n", " ax3.plot(x,sol, 'k', linewidth = 3)\n", " ax3.fill_between(x,-20, sol, facecolor='blue', alpha=0.2)\n", " \n", " return fig\n", " \n", "# Create interactive widget to visualize solution \n", "#interact(plot_Riemann, time=(0,5,0.1), qvar={'sigma11':0, 'sigma22':1, 'sigma12':2, \n", " #'normal velocity u':3, 'transverse velocity v':4});\n", "StaticInteract(plot_Riemann, time=RangeWidget(0,5,0.25)) " ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", " \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", " time: \n", "
    \n", " " ], "metadata": {}, "output_type": "pyout", "prompt_number": 6, "text": [ "" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The solution shows the $x-t$ to show the charactersitics and a horizontal line to show the current time. The other two plots show the solution for the pressure $p$ and velocity $u$ as a function of $x$, and the time slider lets you adjust\n", "the time to observe the time evolution of the solution.\n", "\n", "One way to understand how the middle state $q_M$ is chosen is to observe the phase plane. For instance, in the plot below we show the phase plane where the pressure $p$ corresponds to the $x$ axis and the velocity $u$ corresponds to the $y$ axis. In this plot, the initial left and right states $q_L$ and $q_M$ are points in the plane as well as the new middle state $q_M$. We show the phase plane for some specific parameters. The eigenvector directions are shown in the bottom of the interactive plot. The location of the middle state in the phase plane is obtain by the intersection of the parallelogram created by the eigenvectors and the initial states, which of the two intersections will depend on the location of the initial states. The solution will always move from $q_L$ to $q_M$ in the direction of its corresponding eigenvector $\\bar{r_1}$ and the to $q_R$ in the direction of $\\bar{r_2}$. This is well illustarted in the following interactive app. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "from IPython.display import HTML\n", "HTML('')" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "" ], "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. The normal Riemann problem in 2D\n", "When solving a 2D problem, we can use dimensional splitting or unsplit methods. In any of those cases, we will need to solve a Riemann solver in the normal direction. In the $x$ direction, this will involve solving $\\bar{q}_t + \\mathbf{A}\\bar{q}_x=0$ and in the $y$ direction it will solve $\\bar{q}_t + \\mathbf{B}\\bar{q}_y=0$; we will only solve it in the $x$ direction, since it's completely analogous to solve it for the $y$ direction. The procedure to solve this problem will be the same than for the one-dimensional Riemann problem already explained in detail, so we will present the solution very briefly. We will start with the $\\mathbf{A}$ matrix,\n", "\n", "$$\n", "\\mathbf{A} =\\left[ \\begin{array}{ccc}\n", "0 & K_0 & 0 \\\\\n", "1/\\rho_0 & 0 & 0 \\\\\n", "0 & 0 & 0 \\\\\n", "\\end{array} \\right].\n", "$$\n", "\n", "The matrix of column eigenvectors $\\mathbf{R_x}=[\\bar{r_L}, \\bar{r_0}, \\bar{r_R}]$ is\n", "\n", "$$\n", "\\mathbf{R_x} = \n", "\\left[ \\begin{array}{ccc}\n", "-Z_{0} & 0 & Z_{0} \\\\\n", " 1 & 0 & 1 \\\\\n", " 0 & 1 & 0\n", "\\end{array} \\right],\n", "$$\n", "\n", "with $Z_0=\\rho_0 c_0$ and $c_0=\\sqrt{K_0/\\rho_0}$.\n", "Solving the system $\\mathbf{R}\\bar{\\alpha}=\\delta \\bar{q} =[\\Delta_p, \\Delta_u, \\Delta_v]^T$, we obtain\n", "\n", "$$\n", "\\alpha_L = \\frac{-\\Delta_p + Z_0 \\Delta_u}{2Z_0},\\\\\n", "\\alpha_0 = \\Delta_v,\\\\\n", "\\alpha_R = \\frac{\\Delta_p + Z_0 \\Delta_u}{2Z_0},\n", "$$\n", "\n", "and the speed $s_L = -c_0$, $s_0=0$ and $s_R = c_0$ correspond to the eigenvalues of the matrix $\\mathbf{A}$. Note in this case we assumed the same material on both sides. With this information, we can reconstruct the solution to the Riemann problem as before." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Clawpack implementation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Clawpack has all the framework required to solve hyperbolic equations by solving Riemann problems on each grid cell. If we can provide the Riemann solver to Clawpack, it can do all the hard work for us. In order to implement the solution into Clawpack, we need to rewrite it in its wave formulation, where we update the value at cell $\\bar{Q}_{i,j}^n$, with $i,j$ the indexes in space and $n$ in time following,\n", "\n", "$$\n", "\\bar{Q}_{i,j}^{n+1} = \\bar{Q}_{i,j}^{n} - \\frac{\\Delta t}{\\Delta x}\n", "\\left[A^+\\Delta Q_{i-1/2,j} + A^-\\Delta Q_{i+1/2,j}\\right] + \\text{High order corrections & limiters},\n", "$$ \n", "\n", "when sweeping in the $i$ direction (normal direction) of the computational domain. The same process is analogous in the $j$ direction. Also note this equation requires output of two Riemann problems: the one in edge $i-1/2$ and the one in $i+1/2$. The output of the Riemann solver routine for Clawpack should be given by the waves speed and wave fluctuactions. The wave speed we already know from the eigenvalues in the previous section. In order to calculate the wave fluctuations $A^{\\pm}\\Delta Q_{i-1/2,j}$ from the left edge, we first need the waves given by\n", "\n", "$$\n", "W_i = \\alpha_i \\bar{r_i}\n", "$$\n", "\n", "with $i=L,0,R$ and $\\bar{r_i}$ the corresponding eigenvectors calculated previously. Solving the Riemann problem in the right edge is completely analogous. The positive and negative wave fluctuations $A^\\pm \\Delta Q_{i- 1/2,j}$ will be given in terms of the speeds(eigenvalues) and the waves,\n", "\n", "$$\n", "A^+ \\Delta Q_{i-1/2,j} = s_R W_R \\\\\n", "A^- \\Delta Q_{i-1/2,j} = s_L W_L,\n", "$$\n", "Note the wave $W_0$ doesn't really make any difference in the fluctuations since it has speed zero." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. The transverse Riemann problem in 2D" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to obtain second order accuracy using unsplit methods, we will need to solve a transverse Riemann solver in addition to the normal one. This will require to split the normal wave fluctuactions $A^\\pm \\Delta Q_{i- 1/2,j}$ at edge $i-1/2$ into transverse wave fluctuactions $B^\\pm A^+\\Delta Q_{i- 1/2,j}$ and $B^\\pm A^-\\Delta Q_{i- 1/2,j}$; the former descomposition is shown in the following figure and explained in detail on [LeVeque 2002](http://depts.washington.edu/clawpack/book.html).\n", "\n", "\n", "We will only focuse on $B^\\pm A^+\\Delta Q_{i- 1/2,j}$, since the other case is analogous. We begin by decomposing the normal wave fluctuations into a linear combination of transverse waves, which will follow the same equations and jacobian as before, so we obtain\n", "\n", "$$\n", "A^+\\Delta Q_{i- 1/2,j} = \\mathbf{R_y} \\bar{\\beta} = \\beta_D \\bar{r_D} + \\beta_M \\bar{r_M} + \\beta_U \\bar{r_U},\n", "$$\n", "\n", "where the subindexes $B,M$ and $U$ denote bottom, middle and upper direction. Where $\\bar{r_i}$ are the column eigenvectors of $\\mathbf{B}$, if sweeping in the $x$ direction or of $\\mathbf{A}$ if sweeping on the $y$ direction. In this case, we will only consider the case when sweeping in the $x$ direction, since the other is analogous, so the matrix of eigenvectors of $\\mathbf{B}$ is now $\\mathbf{R_y}=[\\bar{r_B},\\bar{r_M},\\bar{r_U}]$, i.e.\n", "\n", "$$\n", "\\mathbf{R_y} = \n", "\\left[ \\begin{array}{ccc}\n", "-Z_{0} & 0 & Z_{0} \\\\\n", " 0 & 1 & 0 \\\\\n", " 1 & 0 & 1\n", "\\end{array} \\right].\n", "$$\n", "\n", "We solve the system in the same way than the system for $\\bar{\\alpha}$. When we add the right input parameters for the eigenvectors, we obtain $\\bar{\\beta}$ is given by,\n", "\n", "$$\n", "\\beta_B = \\frac{-A^+_1 + z_0 A^+_3}{2 Z_0}, \\\\\n", "\\beta_M = A^+_2, \\\\\n", "\\beta_U = \\frac{A^+_1 + z_0 A^+_3}{2 Z_0}, \\\\\n", "$$\n", "\n", "where the fluctuation vector $A^+\\Delta Q_{i- 1/2,j}=[A^+_1, A^+_2,A^+_3]$, where the speeds \n", "\n", "The speeds (eigenvalues) are \n", "$$\n", "s_B = -c_{0} \\ \\ s_M = 0 \\ \\ s_U = c_{0},\n", "$$\n", "\n", "with $c_0=\\sqrt{K_0/\\rho_0}$. The transverse waves are given by $W_i = \\beta_i \\bar{r_i}$ with $i=B,M,U$, so the fluctuations that will be required to be implemented into Clawpack are given by,\n", "\n", "$$\n", "B^+A^+ \\Delta Q_{i-1/2,j} = s_U W_U \\\\\n", "B^-A^+ \\Delta Q_{i-1/2,j} = s_B W_B.\n", "$$\n", "\n", "Note once again the middle wave has speed zero, so it doesn't affect these fluctuations. \n", "\n", "Now we will show how the normal and transverse flucutuations propagate from the grid edges to the grid cells. We first\n", "obtain the solution for the normal and transverse fluctuation and then we plot them. Note this plot doesn't really show the difference in the vlaues for the pressure and velocities, but it does illustrate the wave speeds and how the normal and transvere fluctuations propagate." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# NORMAL AND TRANSVERSE RIEMANN SOLVER SOLUTION (cartesian grid)\n", "#import numpy\n", "import numpy as np\n", "\n", "# Set left and right state and deltaq (q = [sig11,sig22,sig12,u,v])\n", "qL = np.array([1.0, 0.0, 0.0])\n", "qR = np.array([0.0, 0.0, 0.0])\n", "dq = np.array([qR[0]-qL[0], qR[1]-qL[1],qR[2]-qL[2]])\n", "\n", "# Set normal for normal solver\n", "nx = 1.0\n", "ny = 0.0\n", "# Set normal for transverse solver\n", "nxt = 0.0\n", "nyt = 1.0\n", "\n", "# Set grid size (dx2=dx/2) and time step\n", "dx2 = 0.01\n", "dy2 = 0.01\n", "dt = 0.01\n", "\n", "# Set bulk and density (left and right)\n", "bulkL = 1; rhoL = 1 \n", "bulkR = 8; rhoR = 2\n", "\n", "# Define speeds and impedance (left and right)\n", "cL = np.sqrt(bulkL/rhoL); ZL = rhoL*cL\n", "cR = np.sqrt(bulkR/rhoR); ZR = rhoR*cR\n", "\n", "## NORMAL SOLVER\n", "\n", "# Define the 2 eigenvectors (from columns of Matrix R)\n", "rL = np.array([-ZL, nx, ny])\n", "rR = np.array([ZR, nx, ny])\n", "# Define eigenvalues\n", "sL = -cL\n", "sR = cR\n", "\n", "# Compute the 2 alphas\n", "det = ZL + ZR\n", "alL = (-dq[0] + (nx*dq[1] + ny*dq[2])*ZR)/det\n", "alR = (dq[0] + (nx*dq[1] + ny*dq[2])*ZL)/det\n", "\n", "# Compute wave fluctuations\n", "amdq = alL*rL*sL\n", "apdq = alR*rR*sR\n", "\n", "## TRANSVERSE SOLVER (assuming same material on top and bottom cells)\n", "# Define the 2 eigenvectors (from columns of Matrix R)\n", "rB = np.array([-ZR, nxt, nyt])\n", "rU = np.array([ZR, nxt, nyt])\n", "# Define eigenvalues\n", "sB = -cR\n", "sU = cR\n", "\n", "det = 2.0*ZR\n", "beB = (-apdq[0] + (nxt*apdq[1] + nyt*apdq[2])*ZR)/det\n", "beU = (apdq[0] + (nxt*apdq[1] + nyt*apdq[2])*ZR)/det\n", "\n", "# Compute transverse wave fluctuations\n", "bmdq = beB*rB*sB\n", "bpdq = beU*rU*sU" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 24 }, { "cell_type": "code", "collapsed": false, "input": [ "# CREATE INTERACTIVE PLOT\n", "# Required for plotting\n", "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib.path import Path\n", "import matplotlib.patches as patches\n", "from matplotlib.patches import Ellipse, Polygon\n", "from pylab import *\n", "from IPython.html.widgets import interact\n", "from ipywidgets import StaticInteract, RangeWidget, RadioWidget, DropDownWidget\n", "\n", "# Plot diagram for transverse solver\n", " \n", "def plot_trans_diag(dt,transverse_solver_up,transverse_solver_down):\n", " x = np.linspace(-2*dx2,2*dx2,50);\n", " y1 = 0.0*x + dy2;\n", " y2 = 0.0*x - dy2;\n", "\n", " y = np.linspace(-2*dy2,2*dy2,50);\n", " x1 = 0.0*y + dx2;\n", " x2 = 0.0*y - dx2;\n", " \n", " fig = plt.figure()\n", " ax = fig.add_subplot(111,aspect='equal') \n", " \n", " ax.plot(x,y1,'-k',linewidth=2)\n", " ax.plot(x,y2,'-k',linewidth=2)\n", " ax.plot(x1,y,'-k',linewidth=2)\n", " ax.plot(x2,y,'-k',linewidth=2)\n", " ax.axis([min(x),max(x),min(y),max(y)])\n", " ax.axis('off')\n", " \n", " xplus = sR*dt -dx2\n", " xminus = sL*dt -dx2\n", " yplus = sU*dt\n", " yminus = sB*dt\n", " \n", " # Main patches (A+DQ and A-DQ)\n", " verts = [-dx2,-dy2], [xplus,-dy2], [xplus,dy2], [-dx2,dy2] \n", " poly = patches.Polygon(verts, facecolor='blue', alpha=0.5, linewidth=3) \n", " ax.add_patch(poly)\n", " \n", " verts = [-dx2,-dy2], [xminus,-dy2], [xminus,dy2], [-dx2,dy2] \n", " poly = patches.Polygon(verts, facecolor='blue', alpha=0.5, linewidth=3) \n", " ax.add_patch(poly)\n", " \n", " #Transverse patches\n", " if (transverse_solver_up=='On'):\n", " verts = [-dx2,-dy2], [xplus,-dy2+yplus], [xplus,dy2+yplus], [-dx2,dy2] \n", " poly = patches.Polygon(verts, facecolor='yellow', alpha=0.5, linewidth=3) \n", " ax.add_patch(poly) \n", " \n", " if (transverse_solver_down=='On'):\n", " verts = [-dx2,-dy2], [xplus,-dy2+yminus], [xplus,dy2+yminus], [-dx2,dy2] \n", " poly = patches.Polygon(verts, facecolor='red', alpha=0.5, linewidth=3) \n", " ax.add_patch(poly) \n", " \n", " return fig\n", " \n", "StaticInteract(plot_trans_diag, dt=RangeWidget(0.0,0.008,0.0004), transverse_solver_up=RadioWidget(['On','Off']),\n", " transverse_solver_down=RadioWidget(['On','Off']))" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", " \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", " dt: \n", "
    \n", "transverse_solver_down: On: Off: \n", "
    \n", "transverse_solver_up: On: Off: \n", "
    \n", " " ], "metadata": {}, "output_type": "pyout", "prompt_number": 25, "text": [ "" ] } ], "prompt_number": 25 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. The normal and transverse Riemann problem for quadrilateral mapped grids (and varying density and bulk modulus)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to solve the normal and transverse acoustics Riemann problems for general quadrilateral mapped grids, we will need to solve them in the normal and transverse direction to a grid cell. We can do this by the same process as before, but using a matrix that represents the problem in the normal or transverse direction. The Riemann problem in the normal direction of a grid edge is $\\bar{q}_t + \\mathbf{A_n}\\bar{q}_n=0$, where $\\bar{q}_n$ is the derivative in the normal to the edge direction, $\\mathbf{A_n}=n_x \\mathbf{A} + n_y\\mathbf{B}$, and $\\hat{n}=(n_x,n_y)$ is the normal to the edge where the --normal or transverse-- Riemann problem is being solved. The jacobian matrix is,\n", "\n", "$$\n", "\\mathbf{A_n} =\n", "\\left[ \\begin{array}{ccccc}\n", "0 & n_x K_0 & n_y K_0 \\\\\n", "n_x/\\rho_0 & 0 & 0 \\\\\n", "n_y/\\rho_0 & 0 & 0 \\\\\n", "\\end{array} \\right],\n", "$$\n", "\n", "Using different materials on the left and right side, the matrix of eigenvectors $\\mathbf{R}=[\\bar{r_L},\\bar{r_0},\\bar{r_R}]$ is now\n", "\n", "$$\n", "\\mathbf{R} = \n", "\\left[ \\begin{array}{ccc}\n", "-Z_{L} & 0 & Z_{R} \\\\\n", " n_x & -n_y & n_x \\\\\n", " n_y & n_x & n_y\n", "\\end{array} \\right],\n", "$$\n", "\n", "with $Z_i=\\rho_i c_i$ and $c_i = \\sqrt{K_i/\\rho_i}$, where we now have two different values for the impedance, $i=L,R$,\n", "on the left and right side indicating different materials. The respective eigenvalues are again given by $s_L=-c_L$, $s_0=0$ and $s_R=c_R$.In order to solve the Riemann problem, we need to solve again $\\mathbf{R}\\bar{\\alpha} = \\delta \\bar{q}=[\\Delta_p, \\Delta_u, \\Delta_v]^T$,\n", "\n", "$$\n", "\\alpha_L = \\frac{-\\Delta_p + Z_R(n_x\\Delta_u + n_y\\Delta_v)}{Z_L+Z_R} \\\\\n", "\\alpha_0 = n_x\\Delta_v - n_y \\Delta_u \\\\\n", "\\alpha_R = \\frac{\\Delta_p + Z_L(n_x\\Delta_u + n_y\\Delta_v)}{Z_L+Z_R}\n", "$$\n", "\n", "\n", "The clawpack implementation will when employing the mapped grid we will need to scale the wave speeds $s$ in the wave fluctuations by the factor $\\gamma_x = dy_{phy}/dy_{com}$ if sweeping in the $x$ direction or $\\gamma_y = dx_{phy}/dx_{com}$ when sweeping in the $y$ direction, where $dx_{phy}$ denotes mesh grid spacing in the physical domain and $dx_{com}$ in the computational one. The fluctuation for the normal solver in the mapped grid should then be,\n", "\n", "$$\n", "A^+ \\Delta Q_{i-1/2,j}^{map} = \\gamma s_R \\alpha_R \\bar{r_R}, \\\\\n", "A^- \\Delta Q_{i-1/2,j}^{map} = \\gamma s_L \\alpha_L \\bar{r_L},\n", "$$\n", "\n", "where $\\gamma$ should be $\\gamma_x$ or $\\gamma_y$ depending on the sweeping direction.\n", "Analogously, the up and bottom eigenvectors for the transverse solver will be in principle the same, but with their corresponding normals $\\bar{r_B} = [Z_B,n{Bx},n_{By}]^T$, $\\bar{r_U} = [Z_B,n{Ux},n_{Uy}]^T$, and eigenvalues $s_U=c_U$ and $s_B=-c_B$. Solving the usual system results in,\n", "\n", "$$\n", "\\beta_B = \\frac{-A^+_1 + Z_U(n_{Mx} A^+_2 + n_{My} A^+_3)}{Z_U+Z_B} \\\\\n", "\\beta_U = \\frac{A^+_1 + Z_B(n_{Ux} A^+_2 + n_{Uy} A^+_3)}{Z_U+Z_B},\n", "$$\n", "\n", "where $A^+\\Delta Q_{i- 1/2,j}=[A^+_1, A^+_2,A^+_3]$, $\\hat{n_B}=(n_{Mx},n_{My})$ is the normal to the lower edge of\n", "the middle cell and $\\hat{n_U}=(n_{Ux},n_{Uy})$ is the normal to lower edge of the upper cell.\n", "\n", "\n", "and for the transverse solver,\n", "\n", "$$\n", "B^+A^+ \\Delta Q_{i-1/2,j}^{map} = \\gamma_U s_{U} \\beta_U \\bar{r_U} \\\\\n", "B^-A^+ \\Delta Q_{i-1/2,j}^{map} = -\\gamma_M s_{B} \\beta_B \\bar{r_B},\n", "$$\n", "\n", "where $\\gamma_U$ correspond to the scaling ratio of the bottom edge of the upper cell and $\\gamma_M$ to the scaling ratio of the bottom edge of the middle cell. The same process can be repeated for the left normal fluctuation $A^-\\Delta Q_{i-1/2,j}$. Note the direction of the normal solver is the direction in which we are sweeping, which could be $x$ or $y$ direction of our computational domain.\n" ] } ], "metadata": {} } ] }