|
[Sponsors] |
October 23, 2012, 08:27 |
AUSM on OpenFOAM
|
#1 |
New Member
M.Sabouri
Join Date: Oct 2012
Posts: 2
Rep Power: 0 |
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; } |
|
Tags |
ausm, flux spliting, rhocentralfoam |
|
|
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 12: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 |