CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

implementation of AUSMpw in OpenFOAM

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   October 29, 2021, 06:09
Default implementation of AUSMpw in OpenFOAM
  #1
New Member
 
sanjeev adhikari
Join Date: May 2018
Posts: 13
Rep Power: 8
sanjeev_adhikari is on a distinguished road
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;

}
sanjeev_adhikari is offline   Reply With Quote

Reply

Tags
ausmpw scheme, compressible, openfoam, riemann scheme

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Getting Started with OpenFOAM wyldckat OpenFOAM 25 August 14, 2022 13:55
Map of the OpenFOAM Forum - Understanding where to post your questions! wyldckat OpenFOAM 10 September 2, 2021 05:29
OpenFOAM course for beginners Jibran OpenFOAM Announcements from Other Sources 2 November 4, 2019 08:51
OpenFOAM Training Jan-Jul 2017, Virtual, London, Houston, Berlin CFDFoundation OpenFOAM Announcements from Other Sources 0 January 4, 2017 06:15
Suggestion for a new sub-forum at OpenFOAM's Forum wyldckat Site Help, Feedback & Discussions 20 October 28, 2014 09:04


All times are GMT -4. The time now is 08:39.