implementation of AUSMpw in OpenFOAM

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

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

 Tags ausmpw scheme, compressible, openfoam, riemann scheme

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Linear Mode

 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post wyldckat OpenFOAM 25 August 14, 2022 13:55 wyldckat OpenFOAM 10 September 2, 2021 05:29 Jibran OpenFOAM Announcements from Other Sources 2 November 4, 2019 08:51 CFDFoundation OpenFOAM Announcements from Other Sources 0 January 4, 2017 06:15 wyldckat Site Help, Feedback & Discussions 20 October 28, 2014 09:04

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

 Contact Us - CFD Online - Privacy Statement - Top