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

AUSM on OpenFOAM

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

Like Tree3Likes
  • 3 Post By Sab

Reply
 
LinkBack Thread Tools Display Modes
Old   October 23, 2012, 08:27
Default AUSM on OpenFOAM
  #1
Sab
New Member
 
M.Sabouri
Join Date: Oct 2012
Posts: 2
Rep Power: 0
Sab is on a distinguished road
Dear All,
I tried to implement the AUSM flux splitting method in the OF. I replaced the flux formulation of rhoCentralFoam with those of basic AUSM method(JCP 107, 23-39, 1993). Running the code for forward step case, after some time steps (about 20), some problems occur in some of boundary cells near the outlet that results in negative temperature. I've compared the computed fluxes with those of a fully explicit version of rhoCentralFoam and there is a good agreement. could any one help me to find out the source of this problem?
Thanks.

Here is the source code.



#include "fvCFD.H"
#include "basicPsiThermo.H"
#include "zeroGradientFvPatchFields.H"
#include "fixedRhoFvPatchScalarField.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

int main(int argc, char *argv[])
{
#include "setRootCase.H"

#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "readThermophysicalProperties.H"
#include "readTimeControls.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#include "readFluxScheme.H"

dimensionedScalar v_zero("v_zero",dimVolume/dimTime, 0.0);

Info<< "\nStarting time loop\n" << endl;

while (runTime.run())
{

volScalarField c = sqrt(thermo.Cp()/(thermo.Cv()*psi));
volScalarField rhoa = rho*c;
volVectorField rhoaU = rhoU*c;
volScalarField rhoah = (rhoE+rho/psi)*c;

volVectorField Mach=U/c;
surfaceScalarField Mach_L=fvc::interpolate(Mach, pos, "reconstruct(M)") & (mesh.Sf()/mag(mesh.Sf())) ;
surfaceScalarField Mach_R=fvc::interpolate(Mach, neg, "reconstruct(M)") & (mesh.Sf()/mag(mesh.Sf())) ;

surfaceScalarField Mach_plus_L=
((sign(1.0-mag(Mach_L))+mag(sign(1.0-mag(Mach_L))))/2)
*(0.25*(Mach_L+1.0)*(Mach_L+1.0))
+((sign(-1.0+mag(Mach_L))+mag(sign(-1.0+mag(Mach_L))))/2)
*(0.5*(Mach_L+mag(Mach_L)));

surfaceScalarField Mach_minus_R=
((sign(1.0-mag(Mach_R))+mag(sign(1.0-mag(Mach_R))))/2)
*(-0.25*(Mach_R-1.0)*(Mach_R-1.0))
+((sign(-1.0+mag(Mach_R))+mag(sign(-1.0+mag(Mach_R))))/2)
*(0.5*(Mach_R-mag(Mach_R)));

surfaceScalarField Mach_1_2=Mach_plus_L+Mach_minus_R;

surfaceScalarField p_L=fvc::interpolate(p, pos, "reconstruct(p)");
surfaceScalarField p_R=fvc::interpolate(p, neg, "reconstruct(p)");
surfaceScalarField p_plus_L=
((sign(1.0-mag(Mach_L))+mag(sign(1.0-mag(Mach_L))))/2)
// *(0.25*p_L*(Mach_L+1.0)*(Mach_L+1.0)*(2.0-Mach_L))
*0.5*p_L*(Mach_L+1.0)
+((sign(-1.0+mag(Mach_L))+mag(sign(-1.0+mag(Mach_L))))/2)
*(0.5*p_L*(100*Mach_L+100*mag(Mach_L))/(100*mag(Mach_L)+VSMALL));
surfaceScalarField p_minus_R=
((sign(1.0-mag(Mach_R))+mag(sign(1.0-mag(Mach_R))))/2)
// *(0.25*p_R*(Mach_R-1.0)*(Mach_R-1.0)*(2.0+Mach_R))
*0.5*p_R*(Mach_R-1.0)
+((sign(-1.0+mag(Mach_R))+mag(sign(-1.0+mag(Mach_R))))/2)
*(0.5*p_R*(-100*Mach_R+100*mag(Mach_R))/(100*mag(Mach_R)+VSMALL));

surfaceScalarField p_1_2=p_plus_L+p_minus_R;

surfaceScalarField Direction=sign(Mach_1_2);

surfaceScalarField rhoa_LR=fvc::interpolate(rhoa, Direction, "reconstruct(rho)");
surfaceVectorField rhoaU_LR=fvc::interpolate(rhoaU, Direction, "reconstruct(U)");
surfaceScalarField rhoah_LR=fvc::interpolate(rhoah, Direction, "reconstruct(T)");

surfaceScalarField c_L=fvc::interpolate(c, pos, "reconstruct(T)");
surfaceScalarField c_R=fvc::interpolate(c, neg, "reconstruct(T)");

surfaceScalarField amaxSf("amaxSf",max(max( mag(Mach_L+1.0)*c_L,mag(Mach_R+1.0)*c_R),max(mag(M ach_L-1.0)*c_L,mag(Mach_R-1.0)*c_R))*mag(mesh.Sf()));

#include "compressibleCourantNo.H"
#include "readTimeControls.H"
#include "setDeltaT.H"

runTime++;

Info<< "Time = " << runTime.timeName() << nl << endl;


surfaceScalarField phi("phi", Mach_1_2*rhoa_LR*mag(mesh.Sf()));

surfaceVectorField phiUp =
Mach_1_2*rhoaU_LR*mag(mesh.Sf())+p_1_2*mesh.Sf();

surfaceScalarField phiEp =
(Mach_1_2*rhoah_LR)*mag(mesh.Sf());
// Info << phiEp;
volTensorField tau("tau", mu*dev(fvc::grad(U)()+fvc::grad(U)().T()));

// --- Solve density
solve(fvm::ddt(rho) + fvc::div(phi));


// Info << Mach_1_2;

// --- Solve momentum


volScalarField rhoBydt(rho/runTime.deltaT());
solve(fvm::ddt(rhoU) + fvc::div(phiUp)- fvc::div(tau));
U.dimensionedInternalField() =
rhoU.dimensionedInternalField()
/rho.dimensionedInternalField();
U.correctBoundaryConditions();
rhoU.boundaryField() = rho.boundaryField()*U.boundaryField();

// --- Solve energy
surfaceScalarField sigmaDotU =
(
(
(mesh.Sf() & fvc::interpolate(tau))
)
& (fvc::interpolate(U))
);
volScalarField k("k", thermo.Cp()*mu/Pr);
solve
(
fvm::ddt(rhoE)
+ fvc::div(phiEp)
- fvc::div(sigmaDotU)
- fvc::laplacian(k, T)
);

e = rhoE/rho - 0.5*magSqr(U);
e.correctBoundaryConditions();

thermo.correct();
rhoE.boundaryField() =
rho.boundaryField()*
(
e.boundaryField() + 0.5*magSqr(U.boundaryField())
);

p.dimensionedInternalField() =
rho.dimensionedInternalField()
/psi.dimensionedInternalField();
p.correctBoundaryConditions();
rho.boundaryField() = psi.boundaryField()*p.boundaryField();
// Info << p;
runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return 0;
}
Sab is offline   Reply With Quote

Old   November 9, 2012, 06:29
Default
  #2
New Member
 
Mohanamuraly
Join Date: May 2009
Posts: 16
Rep Power: 8
mohanamuraly is on a distinguished road
Note that the interpolate operator is exactly not the same as the gradient reconstruction that is popularly used in density based codes. It is some sort of linear interpolation using neighboring cell values on to the face.

So I would not recommend using rhoCentralFoam for this purpose ... Maybe I am wrong do correct me ... And its not straight forward to implement Riemann solvers in rhoCentralFoam. Have a look at AeroFoam it looks like a better solver for such applications.
mohanamuraly is offline   Reply With Quote

Old   November 9, 2012, 15:10
Default
  #3
New Member
 
M. Sabouri
Join Date: Nov 2011
Posts: 22
Rep Power: 5
Moslem is on a distinguished road
Thanks.
I applied the AUSM+up scheme and the problem did not repeat.
In fact I was using the upwind scheme to interpolate the cell properties to the faces. So I think that the interpolation scheme was not the source of problem.
Moslem is offline   Reply With Quote

Old   November 9, 2012, 15:20
Default
  #4
New Member
 
M. Sabouri
Join Date: Nov 2011
Posts: 22
Rep Power: 5
Moslem is on a distinguished road
Although I am not using a Riemann solver, you may find the comments by luca_g intersting:


RhoCentralFoam detail
Moslem is offline   Reply With Quote

Old   November 12, 2012, 08:45
Default
  #5
Member
 
Francesco Capuano
Join Date: May 2010
Posts: 78
Rep Power: 7
francesco_capuano is on a distinguished road
Quote:
Originally Posted by Moslem View Post
Thanks.
I applied the AUSM+up scheme and the problem did not repeat.
In fact I was using the upwind scheme to interpolate the cell properties to the faces. So I think that the interpolation scheme was not the source of problem.
Hi, did you eventually get a good result with your AUSM+up scheme into OpenFOAM?

F.
francesco_capuano is offline   Reply With Quote

Old   November 12, 2012, 09:11
Default
  #6
New Member
 
M. Sabouri
Join Date: Nov 2011
Posts: 22
Rep Power: 5
Moslem is on a distinguished road
Hi,
I fact, Results are not so good. But It can capture the main features of the flow ( forward facing step test in M=3 )

Here are some results for t=4.0.
Attached Images
File Type: jpg P.jpg (45.0 KB, 70 views)
File Type: jpg U_mag.jpg (44.7 KB, 64 views)
Moslem is offline   Reply With Quote

Old   November 12, 2012, 09:22
Default
  #7
Member
 
Francesco Capuano
Join Date: May 2010
Posts: 78
Rep Power: 7
francesco_capuano is on a distinguished road
Thanks for your reply. Your results look qualitatively good. Did you compare them with ones obtained with Kurganov-Tadmor scheme and with the results by Woodward&Colella? What about performances, is your code faster than rhoCentralFoam?
francesco_capuano is offline   Reply With Quote

Old   November 12, 2012, 09:40
Default
  #8
New Member
 
Mohanamuraly
Join Date: May 2009
Posts: 16
Rep Power: 8
mohanamuraly is on a distinguished road
Using "upwind scheme to interpolate " means you are simply using the left/right cell value which makes your scheme first order. You must use vanLeer, vanAlbada and the like to have a limited linear reconstruction using gradients to make your FVM discretization second order.

This would require you to specify the gradScheme in your fvSchemes in addition to interpolationScheme.
mohanamuraly is offline   Reply With Quote

Old   November 12, 2012, 10:41
Default
  #9
New Member
 
M. Sabouri
Join Date: Nov 2011
Posts: 22
Rep Power: 5
Moslem is on a distinguished road
Quote:
Originally Posted by mohanamuraly View Post
Using "upwind scheme to interpolate " means you are simply using the left/right cell value which makes your scheme first order. You must use vanLeer, vanAlbada and the like to have a limited linear reconstruction using gradients to make your FVM discretization second order.

This would require you to specify the gradScheme in your fvSchemes in addition to interpolationScheme.
hi mohanamuraly,
Thanks,
I've used such schemes in obtaining the above results.
I'm somewhat confused. while I'm using expressions like :
Code:
    surfaceScalarField  rhoa_LR=fvc::interpolate(rhoa, Direction, "reconstruct(rho)");
in the code and :
Code:
interpolationSchemes
{
    default         linear;
    reconstruct(rho) vanLeer;
    reconstruct(U)  vanLeerV;
    reconstruct(T)  vanLeer;
}
in the fvSchemes,
does the gradScheme affect the cell to face interpolations?
Moslem is offline   Reply With Quote

Old   November 12, 2012, 10:55
Default
  #10
New Member
 
M. Sabouri
Join Date: Nov 2011
Posts: 22
Rep Power: 5
Moslem is on a distinguished road
Quote:
Originally Posted by francesco_capuano View Post
Thanks for your reply. Your results look qualitatively good. Did you compare them with ones obtained with Kurganov-Tadmor scheme and with the results by Woodward&Colella? What about performances, is your code faster than rhoCentralFoam?
In the present step, rhoCentralFoam that uses the Kurganov-Tadmor scheme has better results. For example, the pressure oscillations that are visible in my results, did not appear in the results obtained by rhoCentralFoam . I should work to remove possible bugs in my code. The rhoCentralFoam uses an implicit predictor corrector approach, while my code is fully explicit. For this test case, my code was about 1.6 times faster.
Moslem is offline   Reply With Quote

Old   November 12, 2012, 23:45
Default
  #11
New Member
 
Mohanamuraly
Join Date: May 2009
Posts: 16
Rep Power: 8
mohanamuraly is on a distinguished road
Yes. The "gradScheme" will affect the interpolation. Have a look at the file

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C : lines (38-152)

......
38 template<class Type, class Limiter, template<class> class LimitFunc>
39 tmp<surfaceScalarField> LimitedScheme<Type, Limiter, LimitFunc>::limiter
40 (
41 const GeometricField<Type, fvPatchField, volMesh>& phi
42 ) const
43 {
.....
67
68 GeometricField<typename Limiter::gradPhiType, fvPatchField, volMesh>
69 gradc(fvc::grad(lPhi));
70
......

for limited interpolation OpenFOAM uses the gradient for reconstruction. I was wrong in my previous post about this ... So the gradScheme that you specify will affect your limited interpolation. If "rho" is the variable you interpolate then in gradScheme specify

gradScheme {
grad(rho) leastSquares;
}

to use leastSquares gradient reconstruction. Just figured this out reading the code last night. Hope this helps ...
mohanamuraly is offline   Reply With Quote

Reply

Tags
ausm, flux spliting, rhocentralfoam

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Issues with OpenFoam sanjibdsharma OpenFOAM 0 August 14, 2009 08:41
Problem installing OpenFOAM 1.5 installation on RHEL 4. vwsj84 OpenFOAM Installation 4 April 23, 2009 04:48
2009 OpenFOAM Summer School in Zagreb, Croatia hjasak OpenFOAM Announcements from Other Sources 0 March 27, 2009 13:08
64bitrhel5 OF installation instructions mirko OpenFOAM Installation 2 August 12, 2008 18:07
OpenFOAM Training and Workshop Hrvoje Jasak Main CFD Forum 0 October 7, 2005 07:14


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