flowRateExact.C

Go to the documentation of this file.

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>

using namespace std;


double flowRate
(
    const double D,
    const double rho0,
    const double p0,
    const double p2Byp0
)
{
    const double p2 = p0*p2Byp0;
    const double kappa = 1.4;

    const double etac = pow(2.0/(kappa+1.0), kappa/(kappa-1.0));

    double p1 = etac*p0;
    if(p2 > p1)
    {
        p1 = p2;
    }

    // for a nozzle
    const double Cd0 = 0.84;
    const double Cd1 = 0.66;
    // Not 100% sure this is correct - eqn. (10) in Lbd 5 uses misleading
    // nomenclature
    const double Cdg =
        Cd0 - Cd1*pow(p1/p0, 2.0) + (2*Cd1-Cd0)*pow(p1/p0, 3.0);

    const double Phi =
        sqrt
        (
            kappa/(kappa-1.0)
            *(pow(p1/p0, 2.0/kappa) - pow(p1/p0, (kappa+1.0)/kappa))
        );

    return M_PI*pow(D, 2.0)*Cdg*Phi*sqrt(2.0*rho0*p0);
}


int
main (int argc, char *argv[])
{
    ifstream fi;
    fi.open ("in.txt");
    if(!fi.is_open())
    {
        cerr << "Could not open in.txt!" << endl;
        return 1;
    }

    double D, rho0, p0, p1Byp0;
    fi >> D >> rho0 >> p0 >> p1Byp0;

    double mdot = flowRate(D, rho0, p0, p1Byp0);

    double p1 = p0*p1Byp0;

    cout << "D = " << D
        << " rho0 = " << rho0
        << " p0 = " << p0
        << " p1/p0 = " << p1Byp0
        << " p1 = " << p1
        << " mdot = " << mdot
        << endl;

    ofstream fo;
    fo.open ("out.txt");

    fo << mdot << endl;

    fi.close();
    fo.close();
}