flux.c

Go to the documentation of this file.

 
#include "flux.h"
 
#ifndef CONST_H
    #define gamma   1.4
    #define pmin    1e-15
#endif
 
#define max(a,b) (((a) > (b)) ? (a) : (b))
 
void ns_flux_container(int nfvars, int ndims, double *u, double *nf, double *f)
{   
    double rho = u[0];
    double et = u[nfvars-1];
    double temp = 0.0;
    double contrav = 0.0;
    int i;
 
    for (i=0; i<ndims; i++) {
        contrav += u[i+1]*nf[i];
        temp += u[i+1]*u[i+1];
    }
    contrav /= rho;
 
    // Apply lower bound of pressure value
    double p = (gamma - 1.0)*(et - 0.5*temp/rho);
    if (p < pmin) {
        p = pmin;
        et = p/(gamma-1.0) + 0.5*temp/rho;
        u[nfvars-1] = et;
    }
    
    // Total enthalpy
    double ht = et + p;
 
    // Computes flux array
    f[0] = rho*contrav;
    for (int i=0; i<ndims; i++) {
        f[i+1] = u[i+1] * contrav + nf[i]*p;
    }
    f[nfvars-1] = ht*contrav;
}
 
void rans_flux_container(int nfvars, int ndims, int nturbvars, double *u, double *nf, double *f)
{
    double rho = u[0];
    double contrav = 0.0;
 
    for (int i=0; i<ndims; i++) {
        contrav += u[i+1] * nf[i];
    }
    contrav /= rho;
 
    for (int i=0; i<nturbvars; i++) {
        f[i] = u[nfvars+i]*contrav;
    }
}
 
 
int rans_source_jacobian(int nvars, int ntvars, double betast, \
                         double *uf, double *tmat, double *dsrc)
{
    /* 1-equation RANS model (Spalart-Allmaras) */
    if (ntvars == 1) {
        tmat[0] += dsrc[nvars-1];
    }
 
    /* 2-equations RANS model (kw-SST) */
    else if (ntvars == 2) {
        double k = uf[nvars-2] / uf[0];
        tmat[0] += dsrc[nvars-2];
        tmat[1] += max(betast*k, 0.0);
        tmat[3] += dsrc[nvars-1];
    }
 
    else {
        return -1;
    }
 
    return 0;
}