{ "metadata": { "name": "", "signature": "sha256:8730eef1cbe1d2e6004ec6f6b82fb8803efa972d55efebe0059b308167f9be34" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Riemann solvers for the Euler equations\n", "by Mauricio J. Del Razo S. 2014\n", "\n", "We consider the one dimensional Riemann problem for the Euler equations with the Tammann EOS. We want to solve the one dimensional Euler equations,\n", "\n", "$$\n", "\\begin{align*}\n", " \\left[\\begin{array}{c} \\rho \\\\ \\rho u \\\\ E \\end{array} \\right]_t +\n", " \\left[\\begin{array}{c} \\rho u \\\\ \\rho u^2 + p \\\\ u(E+p) \\end{array} \\right]_x\n", "= 0,\n", "\\end{align*}\n", "$$\n", "where $\\rho$ is density, $u$ velocity, $E$ the internal energy and $p$ the pressure and the subcripts $x,t$ denote partial derivatives with respect $x$ and $t$. The Tamann EOS is given by \n", "$$\n", "p=\\rho e(\\gamma_k -1) - \\gamma_k p_{\\infty k},\n", "$$\n", "\n", "where $e$ the specific internal energy and $k =L,R$ determines which coefficients to use for the EOS. The initial conditions are given by the left and right constant states $\\mathbf{Q}_L=[\\rho_L, \\rho_L u_L, E_L]$ and $\\mathbf{Q}_R=[\\rho_R, \\rho_R u_R, E_R]$. Note that using the equation of state, the state of the system can also be written on the primitive variables $[\\rho, u, p]$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. An exact Riemann solver using the Tammann EOS with a jump in the parameters\n", "In order to solve the Riemann problem, we need to find the zero of $\\phi(p)=0$. The $\\phi(p)$ function is defined below and takes as argument the left and right states in the primitive variables and the left and right parameters of the Tammann equation of state. It returns the value of the function $\\phi$ and the corresponding density and velocity in the middle state." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "\n", "def phi_exact(p, *initial_data):\n", "\n", " rho_l,rho_r,ul,ur,pl,pr,gammal,gammar,pinfl,pinfr = initial_data\n", " \n", " global ustar, rhos_l, rhos_r\n", "\n", " # Bar pressure (SEE IVINGS paper)\n", " pbl = pl + pinfl\n", " pbr = pr + pinfr\n", " # Star bar pressure\n", " pbsl = p + pinfl\n", " pbsr = p + pinfr\n", "\n", " # Useful parameters\n", " gl1 = gammal - 1.0\n", " gr1 = gammar - 1.0\n", " bl = (gammal + 1.0)/(gammal - 1.0)\n", " br = (gammar + 1.0)/(gammar - 1.0)\n", " betal = pbl/bl\n", " betar = pbr/br\n", " al = 2.0/((gammal + 1.0)*rho_l)\n", " ar = 2.0/((gammar + 1.0)*rho_r)\n", " # Calculate velocities\n", " cl = np.sqrt(gammal*(pl + pinfl)/rho_l)\n", " cr = np.sqrt(gammar*(pr + pinfr)/rho_r)\n", " # Calculate the rarefaction phi's \n", " phil_m = ul + 2*cl/gl1*(1 - (pbsl/pbl)**(gl1/(2.0*gammal)))\n", " phir_m = ur - 2*cr/gr1*(1 - (pbsr/pbr)**(gr1/(2.0*gammar)))\n", "\n", " # Use shockwave phi's\n", " phil_p = ul - (p - pl)*np.sqrt(al/(pbsl + betal))\n", " phir_p = ur + (p - pr)*np.sqrt(ar/(pbsr + betar))\n", "\n", " # Assign value to phi function depending on the case\n", " # Rarefaction - Rarefaction (use both intergal curves)\n", " if p <= pl and p <= pr:\n", " phi = phir_m - phil_m\n", " ustar = 0.5*(phir_m + phil_m)\n", " rhos_l = rho_l*(pbsl/pbl)**(1.0/gammal)\n", " rhos_r = rho_r*(pbsr/pbr)**(1.0/gammar)\n", " #print('Rarefaction-Rarefaction')\n", " # Shockwave - Rarefaction (use Hugoniot locus and integral curve)\n", " elif pl <= p and p <= pr:\n", " phi = phir_m - phil_p\n", " ustar = 0.5*(phir_m + phil_p)\n", " rhos_l = rho_l*(pbsl/pbl + 1.0/bl)/(pbsl/(pbl*bl) + 1.0)\n", " rhos_r = rho_r*(pbsr/pbr)**(1.0/gammar)\n", " # print('Shockwave-Rarefaction')\n", " # Rarefaction - Shockwave (use integral curve and hugoniot locus)\n", " elif pr <= p and p <= pl:\n", " phi = phir_p - phil_m\n", " ustar = 0.5*(phir_p + phil_m)\n", " rhos_l = rho_l*(pbsl/pbl)**(1.0/gammal)\n", " rhos_r = rho_r*(pbsr/pbr + 1.0/br)/(pbsr/(pbr*br) + 1.0)\n", " #print('Rarefaction-Shockwave')\n", " # Shockwave - Shockwave (use both hugoniot locus)\n", " elif p >= pr and p >= pl:\n", " phi = phir_p - phil_p\n", " ustar = 0.5*(phir_p + phil_p)\n", " rhos_l = rho_l*(pbsl/pbl + 1.0/bl)/(pbsl/(pbl*bl) + 1.0)\n", " rhos_r = rho_r*(pbsr/pbr + 1.0/br)/(pbsr/(pbr*br) + 1.0)\n", " # print('Shockwave-Shockwave')\n", " return phi" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 66 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we choose some left and right states and parameters and find a zero fo the phi function," ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.optimize import fsolve\n", "\n", "# Define a left state\n", "rho_l = 1.0\n", "ul = 0.0\n", "pl = 3.0\n", "\n", "# Define a right state\n", "rho_r = 0.5\n", "ur = 0.0\n", "pr = 1.0\n", "\n", "# Define left and right parameters for the EOS\n", "gammal = 1.4\n", "gammar = 1.4\n", "pinfl = 0.0\n", "pinfr = 0.0\n", "\n", "# Find the pressure that make phi = 0\n", "initial_data = (rho_l,rho_r,ul,ur,pl,pr,gammal,gammar,pinfl,pinfr)\n", "p0 = (pl + pr)/2.0\n", "pstar = fsolve(phi_exact, p0, args=initial_data)\n", "\n", "# Calculate the left and right wave speeds\n", "betal = (pl + pinfl)*(gammal - 1.0)/(gammal + 1.0)\n", "betar = (pr + pinfr)*(gammar - 1.0)/(gammar + 1.0)\n", "alphal = 2.0/(rho_l*(gammal + 1.0))\n", "alphar = 2.0/(rho_r*(gammar + 1.0))\n", "Sl = ul - np.sqrt((pstar + pinfl + betal)/alphal)/rho_l\n", "Sr = ur + np.sqrt((pstar + pinfr + betar)/alphar)/rho_r" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 89 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We just obtained the pressure and velocity middle states $p_*$ and $u_*$ that are continuous across the contact disconitnuity and the middle states for the density $\\rho_{*L}$ and $\\rho_{*R}$ have also been obtained from the function and saved as global variables. We also know the speed the left and right waves are traveling. We can now create a plotting function for our Riemann solution. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "# NOT YET SET TO PLOT RAREFACTIONS CORRECTLY **\n", "%matplotlib inline\n", "from ipywidgets import StaticInteract, RangeWidget, RadioWidget, DropDownWidget\n", "import matplotlib.pyplot as plt\n", "from pylab import *\n", "\n", "# Create plotting function\n", "def plot_riemann(time):\n", " \n", " # Create figure\n", " fig = plt.figure(figsize=(10, 8))\n", " x = np.linspace(-10,10,100)\n", " \n", " # Create subplot for (x,t) plane\n", " tmax = 5\n", " ax = fig.add_subplot(2,2,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", " # \n", " ax.plot(x,x/Sl, '-r', linewidth=2)\n", " ax.plot(x,x/Sr, '-r', linewidth=2)\n", " # Contact discontinuity\n", " ax.plot(x,x/ustar+0.00000001, '--b', linewidth=2)\n", " # Plot time-line in (x,t) plane\n", " ax.plot(x, 0*x+time, 'k', linewidth=3)\n", " \n", " # Create Riemann solution vector\n", " sol = np.zeros((3,len(x)))\n", " for i in range(len(x)):\n", " if x[i] < Sl*time:\n", " sol[0,i] = rho_l\n", " sol[1,i] = ul\n", " sol[2,i] = pl\n", " elif x[i] < ustar*time:\n", " sol[0,i] = rhos_l\n", " sol[1,i] = ustar\n", " sol[2,i] = pstar\n", " elif x[i] < Sr*time:\n", " sol[0,i] = rhos_r\n", " sol[1,i] = ustar\n", " sol[2,i] = pstar\n", " else:\n", " sol[0,i] = rho_r\n", " sol[1,i] = ur\n", " sol[2,i] = pr\n", "\n", " # Create subplot for density\n", " ax2 = fig.add_subplot(2,2,2)\n", " ax2.set_xlabel('x')\n", " ax2.set_ylabel('density')\n", " ax2.axis([min(x),max(x),-0.5,1.5])\n", " ax2.plot(x,sol[0,:], 'k', linewidth = 3)\n", " ax2.fill_between(x,-20, sol[0,:], facecolor='blue', alpha=0.2)\n", " \n", " # Create subplot for velocity\n", " ax3 = fig.add_subplot(2,2,3)\n", " ax3.set_xlabel('x')\n", " ax3.set_ylabel('velocity')\n", " ax3.axis([min(x),max(x),-0.5,1.5])\n", " ax3.plot(x,sol[1,:], 'k', linewidth = 3)\n", " ax3.fill_between(x,-20, sol[1,:], facecolor='blue', alpha=0.2)\n", " \n", " # Create subplot for pressure\n", " ax4 = fig.add_subplot(2,2,4)\n", " ax4.set_xlabel('x')\n", " ax4.set_ylabel('pressure')\n", " ax4.axis([min(x),max(x),-5,5])\n", " ax4.plot(x,sol[2,:], 'k', linewidth = 3)\n", " ax4.fill_between(x,-20, sol[2,:], facecolor='blue', alpha=0.2)\n", " return fig\n", " \n", " \n", "#middle_states = (rhos_l, rhos_r, ustar, pstar)\n", "#all_states = initial_data + middle_states\n", "\n", "#plot_riemann(2,Sl,Sr,*all_states)\n", "# Create interactive static widget to visualize solution \n", "StaticInteract(plot_riemann, time=RangeWidget(0,5,0.25)) \n", " " ], "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": 117, "text": [ "" ] } ], "prompt_number": 117 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }