The calculated results weren't output as nonuniform List<Scalar>
Hi guys,
I made a Solver for three in-compressible, isothermal immiscible fluids with phase-change based on interPhaseChangeFoam. I deleted the "interfaceProperties.H" in my main program"CrackingPhaseChangeFoam.C", because I don't need to capture the interface among different phases. : These are head files: #include "fvCFD.H" #include "CMULES.H" #include "subCycle.H" //#include "interfaceProperties.H" #include "myphaseChangeThreePhaseMixture.H" #include "turbulentTransportModel.H" #include "pimpleControl.H" #include "fvOptions.H" The solver could run normally: PIMPLE: iteration 1 smoothSolver: Solving for alpha.a, Initial residual = 0.0125493, Final residual = 5.36539e-09, No Iterations 2 Phase-1 volume fraction = 0.654367 Min(alpha.a) = 0.653536 Max(alpha.a) = 0.655 smoothSolver: Solving for alpha.g, Initial residual = 0.193196, Final residual = 5.17285e-09, No Iterations 2 Phase-2 volume fraction = 0.195087 Min(alpha.g) = 0.15 Max(alpha.g) = 0.196981 MULES: Correcting alpha.a MULES: Correcting alpha.a MULES: Correcting alpha.g MULES: Correcting alpha.g MULES: Solving for alpha.c MULES: Solving for alpha.c Phase-1 volume fraction = 0.654367 Min(alpha.a) = 0.653536 Max(alpha.a) = 0.655 Phase-2 volume fraction = 0.195104 Min(alpha.g) = 0.15 Max(alpha.g) = 0.196131 Phase-3 volume fraction = 0.195122 Min(alpha.c) = 0.195 Max(alpha.c) = 0.195283 GAMG: Solving for p_rgh, Initial residual = 0.00316811, Final residual = 0.000128455, No Iterations 1 GAMG: Solving for p_rgh, Initial residual = 0.00122, Final residual = 9.3244e-05, No Iterations 1 GAMGPCG: Solving for p_rgh, Initial residual = 0.000179942, Final residual = 4.28389e-08, No Iterations 3 ExecutionTime = 4.94 s ClockTime = 11 s Courant Number mean: 0.0583806 max: 0.175259 deltaT = 0.005 Time = 3 But the results of fractions of three phases weren't output as nonuniform List<scalar> but in such format: for example: /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; location "0.6"; object alpha.a; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 0 0 0 0 0 0]; internalField uniform 0.655; boundaryField { inlet { type fixedValue; value uniform 0.655; } outlet { type zeroGradient; } walls { type symmetry; } frontAndBack { type empty; } } // ************************************************** *********************** // The calculated fraction of alpha.a is the same as the alpha.a in the "0" file, please point out the reason of this problem, and what should I do to solve it! Thank you!!! As for the IOdictionary of alphas, I just modified library files "two mixtures" into "three mixtures", which located in $FOAM_SRC/transportmodels, the codes are: #include "threePhaseMixture.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::threePhaseMixture::threePhaseMixture ( const fvMesh& mesh, const dictionary& dict ) : phase1Name_(dict.get<wordList>("phases")[0]),/ phase2Name_(dict.get<wordList>("phases")[1]), phase3Name_(dict.get<wordList>("phases")[2]), alpha1_ ( IOobject ( IOobject::groupName("alpha", phase1Name_), mesh.time().timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ), alpha2_ ( IOobject ( IOobject::groupName("alpha", phase2Name_), mesh.time().timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ), alpha3_ ( IOobject ( IOobject::groupName("alpha", phase3Name_), mesh.time().timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ) {} In alphaEqu.H, the calculating steps of alpha1 are: if (MULES1Corr) { fvScalarMatrix alpha1Eqn ( fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1) + fv::gaussConvectionScheme<scalar> ( mesh, phi, upwind<scalar>(mesh, phi) ).fvmDiv(phi, alpha1) - fvm::Sp(divU, alpha1) == fvm::Sp(-tem,alpha1) + em1 ); alpha1Eqn.solve(); Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(" << alpha1.name() << ") = " << min(alpha1).value() << " Max(" << alpha1.name() << ") = " << max(alpha1).value() << endl; talphaPhi1 = alpha1Eqn.flux(); } if (MULES1Corr) { talphaPhi1Corr.ref() -= talphaPhi1(); volScalarField alpha100("alpha100", alpha10); alpha10 = alpha1; MULES::correct ( geometricOneField(), alpha1, talphaPhi1(), talphaPhi1Corr.ref(), Sp, ( divU*(alpha10 - alpha100) + tem*alpha10 )(), oneField(), zeroField() ); // Under-relax the correction for all but the 1st corrector if (aCorr == 0) { talphaPhi1.ref() += talphaPhi1Corr(); } else { alpha1 = 0.5*alpha1 + 0.5*alpha10; talphaPhi1.ref() += 0.5*talphaPhi1Corr(); } } else { MULES::explicitSolve ( geometricOneField(), alpha1, phi, talphaPhi1Corr.ref(), Sp, (divU*alpha1 + em1)(), oneField(), zeroField() ); talphaPhi1 = talphaPhi1Corr; } } Sp, -tem, em1 are sources. |
Sounds like your code is simply not writing the fields out, but is simply copying the files from the 0 folder. Check the coding for the phases - how did you define them? As an IO-Object with AUTO_WRITE?
You need to share some of your coding if you want us to help you further. |
Please check, thanks
Quote:
As for the IOdictionary of the alpha, I just modified the library files "two mixtures" in to "three mixtures" in $FOAM_SRC/transportmodels, which are following: #include "threePhaseMixture.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::threePhaseMixture::threePhaseMixture ( const fvMesh& mesh, const dictionary& dict ) : phase1Name_(dict.get<wordList>("phases")[0]),/ phase2Name_(dict.get<wordList>("phases")[1]), phase3Name_(dict.get<wordList>("phases")[2]), alpha1_ ( IOobject ( IOobject::groupName("alpha", phase1Name_), mesh.time().timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ), alpha2_ ( IOobject ( IOobject::groupName("alpha", phase2Name_), mesh.time().timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ), alpha3_ ( IOobject ( IOobject::groupName("alpha", phase3Name_), mesh.time().timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ) {} As for the alphaEqu.H: Sp, em1 and -tem are the sources. the calculating steps of the first phase is: if (MULES1Corr) { fvScalarMatrix alpha1Eqn ( fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1) + fv::gaussConvectionScheme<scalar> ( mesh, phi, upwind<scalar>(mesh, phi) ).fvmDiv(phi, alpha1) - fvm::Sp(divU, alpha1) == fvm::Sp(-tem,alpha1) + em1 ); alpha1Eqn.solve(); Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(" << alpha1.name() << ") = " << min(alpha1).value() << " Max(" << alpha1.name() << ") = " << max(alpha1).value() << endl; talphaPhi1 = alpha1Eqn.flux(); } if (MULES1Corr) { talphaPhi1Corr.ref() -= talphaPhi1(); volScalarField alpha100("alpha100", alpha10); alpha10 = alpha1; MULES::correct ( geometricOneField(), alpha1, talphaPhi1(), talphaPhi1Corr.ref(), Sp, ( divU*(alpha10 - alpha100) + tem*alpha10 )(), oneField(), zeroField() ); // Under-relax the correction for all but the 1st corrector if (aCorr == 0) { talphaPhi1.ref() += talphaPhi1Corr(); } else { alpha1 = 0.5*alpha1 + 0.5*alpha10; talphaPhi1.ref() += 0.5*talphaPhi1Corr(); } } else { MULES::explicitSolve ( geometricOneField(), alpha1, phi, talphaPhi1Corr.ref(), Sp, (divU*alpha1 + em1)(), oneField(), zeroField() ); talphaPhi1 = talphaPhi1Corr; } } |
Hi I have already add the key information, please check, thank you very much for your reply.
|
That is very strange ... it is almost as if the version of alpha1 etc that you are solving is different from the one that you have generated in the IOobject. Can you quickly list the coding where you define alpha1, alpha2 & alpha3 (in createFields?). Either that, or something has gone wrong with the MULES implementation - you can try debug that by adding in some write statements into the alphaEqn.H to write out the value of alpha1 say in a few cells, and see how that changes through the alphaEqn.H code. Good luck.
|
Quote:
In createFields.H, I defined the alpha1,alpha2 &alpha3 through: Info<< "Creating myphaseChangeThreePhaseMixture\n" << endl; autoPtr<myphaseChangeThreePhaseMixture> mixture = myphaseChangeThreePhaseMixture::New(U, phi); volScalarField alpha1(mixture->alpha1()); volScalarField alpha2(mixture->alpha2()); volScalarField alpha3(mixture->alpha3()); "myphaseChangeThreePhaseMixture" is a class in which I defined the sources em1,Sp,tem and which inherits the "three mixture" class. |
Quote:
#ifndef threePhaseMixture_H #define threePhaseMixture_H #include "volFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { /*---------------------------------------------------------------------------*\ Class twoPhaseMixture Declaration \*---------------------------------------------------------------------------*/ class threePhaseMixture { protected: // Protected data word phase1Name_; word phase2Name_; word phase3Name_; volScalarField alpha1_; volScalarField alpha2_; volScalarField alpha3_; public: // Constructors //- Construct from components threePhaseMixture ( const fvMesh& mesh, const dictionary& dict ); //- Destructor ~threePhaseMixture() {} // Member Functions const word& phase1Name() const { return phase1Name_; } const word& phase2Name() const { return phase2Name_; } const word& phase3Name() const { return phase3Name_; } //- Return the phase-fraction of phase 1 const volScalarField& alpha1() const { return alpha1_; } //- Return the phase-fraction of phase 1 volScalarField& alpha1() { return alpha1_; } //- Return the phase-fraction of phase 2 const volScalarField& alpha2() const { return alpha2_; } //- Return the phase-fraction of phase 2 volScalarField& alpha2() { return alpha2_; } const volScalarField& alpha3() const { return alpha3_; } }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif // ************************************************** *********************** // |
Quote:
Look back at the createFields for interPhaseChangeFoam and you will see that it defines alpha1 etc as refs to alpha1_ not as copies: Code:
volScalarField& alpha1(mixture.alpha1()); |
Quote:
|
Quote:
But again I faced the same problems of my "sources" in the solver, which is only copying the boundary field files from the beginning to the ending of the calculation. I think I am not wrong with the definitions of sources m1,m2,m3 in "createFields.H", which are: Info<< "Creating myphaseChangeThreePhaseMixture\n" << endl; autoPtr<myphaseChangeThreePhaseMixture> mixture = myphaseChangeThreePhaseMixture::New(U, phi); const volScalarField& m1 = mixture->m1(); const volScalarField& m2 = mixture->m2(); const volScalarField& m3 = mixture->m3(); const volScalarField tem(m1+m2+m3); The codes of sources are following: Foam::volScalarField Foam::myphaseChangeThreePhaseMixtures::cracking::m 1() const { const volScalarField limitedAlpha1 ( min(max(alpha1_, scalar(0)), scalar(1)) ); const volScalarField limitedAlpha2 ( min(max(alpha2_, scalar(0)), scalar(1)) ); const volScalarField limitedAlpha3 ( min(max(alpha3_, scalar(0)), scalar(1)) ); return ( -pow(limitedAlpha1*rho1_o_,2)/((limitedAlpha1*rho1_o_+limitedAlpha2*rho2_o_+limi tedAlpha3*rho3_o_)*rho1_o_) *(K1_+K2_) ); } Foam::volScalarField Foam::myphaseChangeThreePhaseMixtures::cracking::m 2() const { const volScalarField limitedAlpha1 ( min(max(alpha1_, scalar(0)), scalar(1)) ); const volScalarField limitedAlpha2 ( min(max(alpha2_, scalar(0)), scalar(1)) ); const volScalarField limitedAlpha3 ( min(max(alpha3_, scalar(0)), scalar(1)) ); return ( pow(limitedAlpha1*rho1_o_,2)/((limitedAlpha1*rho1_o_+limitedAlpha2*rho2_o_+limi tedAlpha3*rho3_o_)*rho2_o_) *(K1_) ); } Foam::volScalarField Foam::myphaseChangeThreePhaseMixtures::cracking::m 3() const { const volScalarField limitedAlpha1 ( min(max(alpha1_, scalar(0)), scalar(1)) ); const volScalarField limitedAlpha2 ( min(max(alpha2_, scalar(0)), scalar(1)) ); const volScalarField limitedAlpha3 ( min(max(alpha3_, scalar(0)), scalar(1)) ); return ( pow(limitedAlpha1*rho1_o_,2)/((limitedAlpha1*rho1_o_+limitedAlpha2*rho2_o_+limi tedAlpha3*rho3_o_)*rho3_o_) *(K2_) ); } K1_ and K2_ are constants. Please help me with it, thank you!!! |
All times are GMT -4. The time now is 18:52. |