[sanjeev_adhikari]

Hello everyone, I was trying to implement AUSMPW in OpenFOAM. Code is posted below. I can run this code however, i got lumps of low temperature at the inlet while using second method while such lumps occurs along the wall while using first method. I have a confusion over use of plLeft and plRight as it violates the range as mentioned in paper while minxy = 1 which is occurs when pLeft and pRight are equal.


Thanks in Advance !!!


#include "AUSMPWplusFlux.H"


void Foam::AUSMPWplusFlux::evaluateFlux
(
scalar& rhoFlux,
vector& rhoUFlux,
scalar& rhoEFlux,

const scalar& pLeft,
const scalar& pRight,

const vector& rhoULeft,
const vector& rhoURight,

const scalar& rhoLeft,
const scalar& rhoRight,

const scalar& aLeft,
const scalar& aRight,

const scalar& rhoELeft,
const scalar& rhoERight,



const vector& Sf,
const scalar& magSf

) const
{
const scalar alpha = 3.0/16.0;
const scalar beta = 1.0/8.0;


// normal vector
const vector normalVector = Sf/magSf;

const scalar rhoHLeft = rhoELeft + pLeft;



const scalar rhoHRight = rhoERight + pRight;


const vector ULeft=rhoULeft/rhoLeft;
const vector URight=rhoURight/rhoRight;

const scalar qLeft = ULeft & normalVector;
const scalar qRight = URight & normalVector;

const scalar aTilde = 0.5*(aLeft+aRight);

const scalar MaRelLeft = qLeft /aTilde;
const scalar MaRelRight = qRight/aTilde;

const scalar magMaRelLeft = mag(MaRelLeft);
const scalar magMaRelRight = mag(MaRelRight);


// *** Mach number spliting functions ***

const scalar Ma1PlusLeft = 0.5*(MaRelLeft +magMaRelLeft ); //0.5*(M+|M|)
const scalar Ma1MinusRight = 0.5*(MaRelRight-magMaRelRight); //0.5*(M+|M|)

const scalar Ma2PlusLeft = 0.25*sqr(MaRelLeft +1.0);
const scalar Ma2PlusRight = 0.25*sqr(MaRelRight+1.0);
const scalar Ma2MinusLeft = -0.25*sqr(MaRelLeft -1.0);
const scalar Ma2MinusRight = -0.25*sqr(MaRelRight-1.0);

const scalar Ma4BetaPlusLeft = ((magMaRelLeft >= 1.0) ? Ma1PlusLeft : (Ma2PlusLeft *(1.0-16.0*beta*Ma2MinusLeft)));
const scalar Ma4BetaMinusRight = ((magMaRelRight >= 1.0) ? Ma1MinusRight : (Ma2MinusRight*(1.0+16.0*beta*Ma2PlusRight)));

/* const scalar Ma4BetaPlusLeft = ((magMaRelLeft > 1.0) ? Ma1PlusLeft : (Ma2PlusLeft *(1.0-16.0*beta*Ma2MinusLeft)));
const scalar Ma4BetaMinusRight = ((magMaRelRight > 1.0) ? Ma1MinusRight : (Ma2MinusRight*(1.0+16.0*beta*Ma2PlusRight)));
*/
// *** Pressure splitting functions ***

const scalar P5alphaPlusLeft = ((magMaRelLeft >= 1.0) ?
(Ma1PlusLeft/MaRelLeft) : (Ma2PlusLeft *(( 2.0-MaRelLeft) -16.0*alpha*MaRelLeft *Ma2MinusLeft )));
const scalar P5alphaMinusRight = ((magMaRelRight >= 1.0) ?
(Ma1MinusRight/MaRelRight) : (Ma2MinusRight*((-2.0-MaRelRight)+16.0*alpha*MaRelRight*Ma2PlusRight)));

/* const scalar P5alphaPlusLeft = ((magMaRelLeft > 1.0) ?
(Ma1PlusLeft/MaRelLeft) : (Ma2PlusLeft *(( 2.0-MaRelLeft) -16.0*alpha*MaRelLeft *Ma2MinusLeft )));
const scalar P5alphaMinusRight = ((magMaRelRight > 1.0) ?
(Ma1MinusRight/MaRelRight) : (Ma2MinusRight*((-2.0-MaRelRight)+16.0*alpha*MaRelRight*Ma2PlusRight)));

*/


// *** inteface Mach number, pressure, and velocity ***


const scalar pTilde = P5alphaPlusLeft*pLeft + P5alphaMinusRight*pRight;

const scalar minXY = min((pLeft/pRight), (pRight/pLeft));

const scalar plLeft = ((minXY >= 0.0 && minXY < 0.75) ? 0.0 : 4.0*minXY-3.0);

/* scalar plLeft = 0.0;
scalar plRight = 0.0;*/


/* if ((minXY >= 0.0 && minXY < 0.75))
{
plLeft = 0.0;
}
else if ((minXY >= 0.75 && minXY < 1))
{
plLeft = 4.0*minXY-3.0;
}*/
/* else
{
Info << "\t AUSMPW+ flux solver: Error !!! plLeft " << endl;
}*/


const scalar minYX = min((pRight/pLeft),(pLeft/pRight));
const scalar plRight = ((minYX >= 0.0 && minYX < 0.75) ? 0.0 : 4.0*minYX-3.0);

// Info << " Min value min(x,y)" << minXY << endl;
// Info << " Min value min(y,x)" << minYX << endl;

/* if ((minYX >= 0.0 && minYX < 0.75))
{
plRight = 0.0;
}
else if ((minYX >= 0.75 && minYX < 1))
{
plRight = 4.0*minXY-3.0;
}*/
/* else
{
Info << "\t AUSMPW+ flux solver: Error !!! plRight " << endl;
}*/



const scalar vLeft = pow((mag(ULeft)/aTilde), 0.25);
const scalar vRight = pow((mag(URight)/aTilde), 0.25);






scalar fLeft = ((magMaRelLeft > 1.0) ? 0.0 : (pLeft/pTilde -1.0)*plLeft*mag(Ma2PlusLeft)* min(1.0,vLeft));
scalar fRight = ((magMaRelRight > 1.0) ? 0.0 : (pRight/pTilde -1.0)*plRight*mag(Ma2MinusRight)* min(1.0,vRight));

/* const scalar fLeft = ((magMaRelLeft > 1.0) ? 0.0 : (pLeft/pTilde -1.0)*plLeft*mag(Ma2PlusLeft)* min(1.0,vLeft));
const scalar fRight = ((magMaRelRight >1.0) ? 0.0 : (pRight/pTilde -1.0)*plRight*mag(Ma2MinusRight)* min(1.0,vRight));
*/

const scalar weightFn = 1.0 - pow(minXY, 3.0);



/* //first method
const scalar rhoFluxGLeft = (1.0- weightFn)*rhoLeft + weightFn*rhoRight;
const vector rhoUFluxGLeft = (1.0- weightFn)*rhoULeft + weightFn*rhoURight;
const scalar rhoHFluxGLeft = rhoFluxGLeft*rhoHLeft;


const scalar rhoFluxGRight = (1.0- weightFn)*rhoRight + weightFn*rhoLeft;
const vector rhoUFluxGRight = (1.0- weightFn)*rhoURight + weightFn*rhoULeft;
const scalar rhoHFluxGRight = rhoFluxGRight*rhoHRight;*/

// *** interface flux calculation ***

const scalar MaTilde = (Ma4BetaPlusLeft + Ma4BetaMinusRight);

/* //first method
const scalar MaTildeLeft = (1.0+fLeft)*Ma4BetaPlusLeft*aTilde;
const scalar MaTildeRight = (1.0+fRight)*Ma4BetaMinusRight*aTilde;*/

//second method
const scalar weightFnRight = weightFn*(1.0+fRight);
const scalar weightFnLeft = weightFn*(1.0+fLeft);

const scalar MaTildeAvg = fLeft*Ma4BetaPlusLeft+fRight*Ma4BetaMinusRight;


const scalar MaTildePlusLeft = ((MaTilde >= 0.0) ?
(MaTilde- Ma4BetaMinusRight*weightFnRight+ MaTildeAvg) : (Ma4BetaPlusLeft*weightFnLeft));
const scalar MaTildePlusRight = ((MaTilde >= 0.0) ?
(MaTilde - Ma4BetaPlusLeft*weightFnLeft+MaTildeAvg) : (Ma4BetaMinusRight*weightFnRight));

const scalar UTildaPlusLeft = MaTildePlusLeft*aTilde;
const scalar UTildePlusRight = MaTildePlusRight*aTilde;

// *** interface flux calculation ***

/* //first method
rhoFlux = ((MaTilde > 0.0) ? (MaTildeLeft*rhoLeft+MaTildeRight*rhoFluxGLeft)*ma gSf : (MaTildeLeft*rhoFluxGRight+MaTildeRight*rhoRight)* magSf);
rhoUFlux = ((MaTilde > 0.0) ? (MaTildeLeft*rhoULeft+MaTildeRight*rhoUFluxGLeft + pTilde*normalVector)*magSf : (MaTildeLeft*rhoUFluxGRight+MaTildeRight*rhoURight + pTilde*normalVector)*magSf);
rhoEFlux = ((MaTilde > 0.0) ? (MaTildeLeft*rhoHLeft+MaTildeRight*rhoHFluxGLeft)* magSf : (MaTildeLeft*rhoHFluxGRight+MaTildeRight*rhoHRight )*magSf);
*/

//second method

rhoFlux = (UTildaPlusLeft*rhoLeft + UTildePlusRight*rhoRight)*magSf;
rhoUFlux = (UTildaPlusLeft*rhoULeft + UTildePlusRight*rhoURight + pTilde*normalVector)*magSf;
rhoEFlux = (UTildaPlusLeft*rhoHLeft + UTildePlusRight*rhoHRight)*magSf;

//Info << "\t rhoFlux = " << rhoFlux << endl;
//Info << "\t rhoUFlux = " << rhoUFlux << endl;
//Info << "\t rhoEFlux = " << rhoEFlux << endl;

//Info << "\nAUSMPW flux update completed...\n" << endl;

}