ConicalShockRelations.js

//******************************************************************************
// Stephen Krauss
//
// NAME: ConicalShockRelations.js
//
// PURPOSE: Provides functions to calculate properties of the shock wave
//          resulting from supersonic compressible flow around a right circular
//          cylinder. Uses common compressible flow relations and the
//          Taylor-Maccoll differential equations for flow around a cylinder.
//******************************************************************************

//--------------------------------------------------------------------------
// Name:        getFlowDeflectionAngle
// Description: Uses the delta-beta-m equation to calculate flow deflection
//              angle caused by flow through an oblique shock wave.
// Arguments:   - gamma, ratio of specific heats for the flow
//              - m1, flow Mach number
//              - beta, oblique shock wave angle in radians
// Returns:     - delta, flow deflection angle in radians
//--------------------------------------------------------------------------
function getFlowDeflectionAngle(gamma, m1, beta) 
{
    
    var delta;
    
    // Calculate the flow deflection angle from the delta-beta-m equation
    delta =  Math.atan((m1*m1 * Math.sin(2.0*beta) - 2.0/Math.tan(beta))
                       /(2.0 + m1*m1 * (gamma + Math.cos(2.0*beta))));
   
    return delta;
    
} // getDelta

//--------------------------------------------------------------------------
// Name:        getShockAngle
// Description: Numerically solves the Taylor-Maccoll differential equations
//              to calculate the shock angle generated by supersonic flow
//              over a cone.
// Arguments:   - gamma, ratio of specific heats for the flow
//              - m1, flow Mach number
//              - sigma, cone angle in radians
// Returns:     - beta, shock angle in radians
//--------------------------------------------------------------------------
function getShockAngle(gamma, m1, sigma) 
{

    // Angle increment in radians. Controls speed and accuracy of 
    // the calculation.
    var ANGLE_INC = 0.0005;

    // Error check the inputs
    if(gamma <= 1.0) 
    {
        alert("ERROR: Gamma must be greater than 1!");
        return;
    } // if
    
    if(m1 <= 1.0) 
    {
        alert("ERROR: m1 must be greater than 1!");
        return;
    } // if
    
    if(sigma > Math.PI/2.0) 
    {
        alert("ERROR: Cone angle must be less than 90 degrees!");
        return;
    } // if
    
    if(sigma < 0.0) {
        alert("ERROR: Cone angle must be greater than 0 degrees!");
        return;
    } // if
        
    // Guess a wave angle to start searching from,
    // use Mach angle since (Mach angle)<(wave angle)<(90deg)
    beta = Math.asin(1.0/m1);
    
    // Get the cone angle resulting from this guess
    sigma0 = getConeAngle(gamma, m1, beta);
    
    // Increase wave angle guess by the angle increment and calculate 
    // the resulting cone angle until the calculated cone angle matches 
    // the input cone angle.
    while(sigma0 < sigma) {

        beta = beta + ANGLE_INC;
        sigma0 = getConeAngle(gamma, m1, beta);
        
    } // while
        
    return beta
        
} // getShockAngle

//--------------------------------------------------------------------------
// Name:        getdVr
// Description: Calculates the dVr component of the Taylor-Maccoll
//              differential equations.
// Arguments:   - gamma, ratio of specific heats for the flow
//              - beta, shock angle in radians
//              - Vr, non-dimensional Vr flow velocity component
//              - Vtheta, non-dimensional Vtheta flow velocity component
// Returns:     - dVr, evaluated Taylor-Maccoll equation
//--------------------------------------------------------------------------
function getdVr(gamma, beta, Vr, Vtheta) 
{

    dVr = Vtheta;
    
    return dVr;
    
} // getdVr

//--------------------------------------------------------------------------
// Name:        getdVtheta
// Description: Calculates the dVtheta component of the Taylor-Maccoll
//              differential equations.
// Arguments:   - gamma, ratio of specific heats for the flow
//              - beta, shock angle in radians
//              - Vr, non-dimensional Vr flow velocity component
//              - Vtheta, non-dimensional Vtheta flow velocity component
// Returns:     - dVtheta, evaluated Taylor-Maccoll equation
//--------------------------------------------------------------------------
function getdVtheta(gamma, beta, Vr, Vtheta) 
{
     
    dVtheta = (Vtheta*Vtheta*Vr-(gamma-1.0)/2.0*(1.0-Vr*Vr-Vtheta*Vtheta)
               *(2.0*Vr+Vtheta*(1.0/Math.tan(beta)))) 
               /((gamma-1.0)/2.0*(1.0-Vr*Vr-Vtheta*Vtheta)-Vtheta*Vtheta);

              
    return dVtheta
    
} // getdVtheta

//--------------------------------------------------------------------------
// Name:        getConeAngle
// Description: Numerically solves the Taylor-Maccoll differential equations
//              to calculate the cone angle that would cause a supersonic
//              flow to create a given shock wave angle.
// Arguments:   - gamma, ratio of specific heats for the flow
//              - m1, flow Mach number
//              - beta, shock angle in radians
// Returns:     - sigma, cone angle in radians
//--------------------------------------------------------------------------
function getConeAngle(gamma, m1, beta) 
{
    
    // Angle increment in radians. Controls speed and accuracy of 
    // the calculation.
    var ANGLE_INC = -0.0005;
    
    // Error check the inputs
    if(gamma <= 1.0) 
    {
        alert("ERROR: Gamma must be greater than 1!");
        return;
    } // if
    
    if(m1 <= 1.0) 
    {
        alert("ERROR: m1 must be greater than 1!");
        return;
    } // if
    
    if(beta > Math.PI/2.0) 
    {
        alert("ERROR: Shock angle must be less than 90 degrees!");
        return;
    } // if
    
    if(beta < 0.0) {
        alert("ERROR: Shock angle must be greater than 0 degrees!");
        return;
    } // if

    var theta0 = beta;

    // Calculate the flow deflection angle right at the shock       
    var delta = getFlowDeflectionAngle(gamma, m1, beta);
        
    // Calculate the Mach number after the shock using Mach number components
    // normal to the shock and the normal shock relations
    var m1n = m1*Math.sin(beta);
    var m2n = Math.sqrt((1.0 + 0.5 * (gamma - 1.0) * m1n * m1n) / (gamma * m1n * m1n - 0.5 * (gamma - 1.0)));
    var m2 = m2n/Math.sin(beta-delta);
    
    // Calculate non-dimensional components of the velocity 
    // in theta and r direction. These are initial values for the ODE
    var V0 = Math.pow((1.0+2.0/((gamma-1.0)*m2*m2)),(-0.5));
    var Vr0 = V0 * Math.cos(theta0-delta);
    var Vtheta0 = -V0 * Math.sin(theta0-delta);
    
    // Using a Runge-Kutta 4th order ODE estimation, step down the shock angle
    // by the angle increment until the the Vtheta component is zero,
    // which means we've reached the cone surface.
    while(Vtheta0 < 0.0)
    {
                
        // Calculate the Runge-Kutta 4th order constants        
        var k1 = getdVr(gamma, theta0, Vr0, Vtheta0);
        var j1 = getdVtheta(gamma, theta0, Vr0, Vtheta0);
        
        var k2 = getdVr(gamma, theta0+ANGLE_INC/2.0, Vr0+(ANGLE_INC/2.0)*k1, Vtheta0+(ANGLE_INC/2.0)*j1);
        var j2 = getdVtheta(gamma, theta0+ANGLE_INC/2.0, Vr0+(ANGLE_INC/2.0)*k1, Vtheta0+(ANGLE_INC/2.0)*j1);
        
        var k3 = getdVr(gamma, theta0+ANGLE_INC/2.0, Vr0+(ANGLE_INC/2.0)*k2, Vtheta0+(ANGLE_INC/2.0)*j2);
        var j3 = getdVtheta(gamma, theta0+ANGLE_INC/2.0, Vr0+(ANGLE_INC/2.0)*k2, Vtheta0+(ANGLE_INC/2.0)*j2);
        
        var k4 = getdVr(gamma, theta0+ANGLE_INC, Vr0+ANGLE_INC*k3, Vtheta0+ANGLE_INC*j3);
        var j4 = getdVtheta(gamma, theta0+ANGLE_INC, Vr0+ANGLE_INC*k3, Vtheta0+ANGLE_INC*j3);
        
        // Use the Runge-Kutta constants to estimate Vr and Vtheta
        Vr0 = Vr0 + (ANGLE_INC/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4);
        Vtheta0 = Vtheta0 + (ANGLE_INC/6.0) * (j1 + 2.0*j2 + 2.0*j3 + j4);
        
        theta0 = theta0 + ANGLE_INC;
        
    } // while

    // The angle at which Vtheta equals zero is the cone angle
    var sigma = theta0;
    
    return sigma
    
} // getConeAngle