# alpha1 wrong behaviour with LTSInterFoam

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

 July 2, 2014, 08:36 alpha1 wrong behaviour with LTSInterFoam [solved] #1 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 Hi Foamers, I' m setting up (open FOAM v 2.3.0) a case for wave generation and drag calculation using LTSInterFoam. I conduct some basics verification and validation of my set up with a hull, for wich I had the drag values ; and I was able to get a 5% max difference with those value, with a visually coherent wave pattern. Now I'm trying to apply my set-up wih a new hull - wich had quite the same bow as the mini 6.50 747 magnum. The case brake up as the alpha.water get a very wrong behaviour under the hull : outside the boundary layers, alpha seems quite correct but inside the boundary layer it goes to zero near the hull. I've attached some picture with a "cylinder-hull" illustrating this, wich i plotted the interface, colored by z ; and the value of alpha.water on the "cylinder-hull". blue color is full water and red is full air. 1 is the initial situation. 1.jpg 2, 3 and 4 are ploted from the iteration 1600, we can see that under the "cylinder-hull" there is a water-air mix, when it should be full water. 3.jpg 2.jpg the 4 picture show that the matter only occur on the boundary layer - when i put more boundary layers, only the first 4-5 boundary layer on the wall get those wrong alpha.water values. 4.jpg. I also have upload my test case set up here. vSolution : Code: solvers { "alpha.water.*" { nAlphaCorr 0; nAlphaSubCycles 1; cAlpha 0; icAlpha 0; MULESCorr yes; nLimiterIter 10; alphaApplyPrevCorr yes; solver smoothSolver; smoother symGaussSeidel; tolerance 1e-8; relTol 0; minIter 1; } pcorr { solver PCG; preconditioner { preconditioner GAMG; smoother GaussSeidel; agglomerator faceAreaPair; mergeLevels 1; nCellsInCoarsestLevel 100; cacheAgglomeration true; tolerance 1e-6; relTol 0; }; tolerance 1e-6; relTol 0; }; p_rgh { solver GAMG; smoother DIC; agglomerator faceAreaPair; mergeLevels 1; nCellsInCoarsestLevel 100; cacheAgglomeration true; tolerance 1e-8; relTol 0.0; }; p_rghFinal { \$p_rgh; relTol 0; } "(U|k|omega).*" { solver smoothSolver; smoother symGaussSeidel; nSweeps 1; tolerance 1e-7; relTol 0; minIter 1; }; } PIMPLE { momentumPredictor yes; nOuterCorrectors 1; nCorrectors 2; nNonOrthogonalCorrectors 0; maxCo 5; maxAlphaCo 5; rDeltaTSmoothingCoeff 0.05; rDeltaTDampingCoeff 0.5; nAlphaSpreadIter 0; nAlphaSweepIter 0; maxDeltaT 1; } relaxationFactors { fields { p_rgh 0.3; } equations { "(U|k|omega).*" 0.7; "alpha.water.*" 0.7; } } fvScheme Code: ddtSchemes { default localEuler rDeltaT; } gradSchemes { default cellMDLimited Gauss linear 0.5;//Gauss linear;// limitedGrad cellLimited Gauss linear 1; } divSchemes { div(rhoPhi,U) Gauss linearUpwind grad(U); div(phi,alpha) Gauss limitedLinear 1.0 phi ;//Gauss vanLeer; div(phirb,alpha) Gauss interfaceCompression ; div(phi,k) Gauss linearUpwind limitedGrad; div(phi,omega) Gauss linearUpwind limitedGrad; div((muEff*dev(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear limited 0.777; } interpolationSchemes { default linear; } snGradSchemes { default limited 0.777; } fluxRequired { default no; p_rgh; pcorr; alpha.water; } i have already try to change settings there and there, but as I'm unable to diagnose the matter, i've made only hazardous changes. thank you in advance for your help, Quentin Coispeau ps : in my cylinder-test case, the mesh in not high enough in z+ to get the full wave at the bow, but i had the very same alpha problem with my "magnum-like" hull, for wich it was not the case. Last edited by Quentin; July 30, 2014 at 08:48.

 July 3, 2014, 05:06 #2 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 I found out that : calpha can't be set to 0 : the compression speed used to compressed the interface is calculate as |U_c| = min( calpha * |U|, max |U| ). so if calpha is set to 0, you remove the convective compression flux and alpha equation becomes too diffusive. calpha should be equal to 1 or greater. I also had to had a alpha corrector and set the alpha div scheme to vanLeer again. And reduce maxCo and maxAlphaCo

 July 3, 2014, 10:45 #3 Member   laurentL Join Date: Oct 2011 Location: new caledonia Posts: 72 Rep Power: 6 hi Quentin, did it solve all your pb? i just run the beginning of you case. i call this phenomena 'air under hull'. i got this behevior on a catamaran quiet flat hull. i could solve it in Finemarine with a threshold in pressure. i plan to make the same in OF, not done yet... air under hull, alpha.water=1 threshold pressure Laurent

 July 3, 2014, 11:41 #4 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 Hi Laurent, Try to run the same case with these dictionnaries instead of the previous ones : fvSolution Code: solvers { "alpha.water.*" { nAlphaCorr 1; nAlphaSubCycles 1; cAlpha 1; icAlpha 0; MULESCorr yes; nLimiterIter 10; alphaApplyPrevCorr yes; solver smoothSolver; smoother symGaussSeidel; tolerance 1e-8; relTol 0.1; minIter 1; } pcorr { solver PCG; preconditioner { preconditioner GAMG; smoother GaussSeidel; agglomerator faceAreaPair; mergeLevels 2; nCellsInCoarsestLevel 10; cacheAgglomeration false; nPreSweeps 0; nPostSweeps 2; nBottomSweeps 2; tolerance 1e-5; relTol 0; }; tolerance 1e-5; relTol 0.01; maxIter 100; }; p_rgh { solver GAMG; smoother GaussSeidel; agglomerator faceAreaPair; mergeLevels 1; nCellsInCoarsestLevel 10; cacheAgglomeration true; nPreSweeps 0; nPostSweeps 2; nBottomSweeps 2; tolerance 1e-5; relTol 0.05; }; p_rghFinal { solver PCG; preconditioner { preconditioner GAMG; tolerance 1e-07; relTol 0; nVcycles 2; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; cacheAgglomeration on; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; } tolerance 1e-07; relTol 0; maxIter 20; } "(U|k|omega).*" { solver smoothSolver; smoother symGaussSeidel; nSweeps 1; tolerance 1e-7; relTol 0.1; minIter 1; }; } PIMPLE { momentumPredictor yes; nOuterCorrectors 1; nCorrectors 2; nNonOrthogonalCorrectors 0; maxCo 0.2499; maxAlphaCo 0.2499; rDeltaTSmoothingCoeff 0.05; rDeltaTDampingCoeff 0.5; nAlphaSpreadIter 0; nAlphaSweepIter 1; maxDeltaT 1; } relaxationFactors { /* fields { p_rgh 0.3; } equations { "(U|k|omega).*" 0.7; "alpha.water.*" 0.7;*/ } } cache { grad(U); } fvSchemes Code: ddtSchemes { default localEuler rDeltaT; } gradSchemes { default Gauss linear;//cellMDLimited Gauss linear 0.5;// limitedGrad cellLimited Gauss linear 1; } divSchemes { div(rhoPhi,U) Gauss linearUpwind grad(U); div(phi,alpha) Gauss vanLeer; div(phirb,alpha) Gauss interfaceCompression ; div(phi,k) Gauss upwind limitedGrad; div(phi,omega) Gauss upwind limitedGrad; div((muEff*dev(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear limited 0.777; } interpolationSchemes { default linear; } snGradSchemes { default limited 0.777; } fluxRequired { default no; p_rgh; pcorr; alpha.water; } Using calpha = 1 permit to get a much sharper interface and quite solve the problem (behind the bow I still have the first layer of the boundary layer which got wrong values as in the pictures below after 6000 times step 6.jpg 7.jpg But an other matter "the wiggle wave" appears as related in this thread. But as we need the boundary layer to resolve the boundary layer, we can't get a better aspect ratio or am I wrong ?) If you look in fvSolution, i change maxC0 and maxAlphaCo from 5 to 0.2499, so I need 20 times more timeStep, which is quite longer ; i tried to add more alpha subcycles, but the interface becomes quite crazy after not so much timeStep. you're idea of pressure threshold seems quite interesting to me, and should probably be implemented as a new boundary condition. I may give a try in it soon if I can't get satisfactory results without. But i fear it is a bit "cheating" with numerical behaviour, which may false the results - at least it needs validation. I set icalpha to 0.25 because I've seen this value somewhere (don't remember where) but I don't know what's its use and haven't look through the code for it yet.

 July 3, 2014, 11:53 #5 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 Here an illustration of you're pressure treshold idea : the pressure and the interface on the cylinder with the blue/green jump at p = 250 Pa 8.jpg so it seems possible to implement your pressure condition on alpha quite easely.

 July 3, 2014, 17:46 #6 Member   laurentL Join Date: Oct 2011 Location: new caledonia Posts: 72 Rep Power: 6 hi, so, let's say, if p superior at 250 then put alpha=1 but i'am a very beginner in C++... is anyone have a idea of how to do that??? 1- add a variable to declare in transportProperties P0=250 2-add a line in /home/laurent/OpenFOAM/OpenFOAM-2.3.0/applications/solvers/multiphase/interFoam/alphaEqn.H like if (p > P0) { alpha.water=1 } 3- compile could it be SO simple.... i'm dreaming.... any advice will be very welcome! thank to all Laurent

 July 4, 2014, 03:30 #7 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 I think it would be much better to implement a new boundary condition which would be something like Code: if (p < P0) zeroGradient; else fixedValue 1; end I think I will do it today because I still have a one layer thick wrong alpha on my real hull simulation

 July 13, 2014, 12:08 #8 Member   Dave Join Date: Jul 2010 Posts: 97 Rep Power: 8 Dear Quentin, The behavior you are seeing is called "numerical ventilation". The forcing of the mass fraction alpha to be 1 when the pressure is greater than a certain value like can be done in FINE/Marine is something that should be used with caution. With low speed vessels where real ventilation is unlikely, this can be a way to improve your answer. For high speed craft you should confirm ventilation is numerical instead of real ventilation before applying such a correction. Also be careful about the value of pressure you pick being greater than the stagnation pressure of the air or you will create water in places you don't want it. I should point out that you can correct the friction forces in post treatment by correcting the rho used. This is only appropriate in the case that wall functions were used. Both have their downsides. The best is to avoid the issue by refining the mesh near the where the problem starts and minimizing sources of diffusion (schemes, courant number, etc). Easier said than done since they tend to increase the cost.

July 15, 2014, 03:07
#9
New Member

Quentin Coispeau
Join Date: Nov 2013
Posts: 17
Rep Power: 4
Hi Dave,

thank you very much for your answer. knowing the name of the issue will help me find documentation about it.

I agree with you saying that the pressure threshold trick is more a bandage than a solution.

Quote:
 I should point out that you can correct the friction forces in post treatment by correcting the rho used. This is only appropriate in the case that wall functions were used. Both have their downsides.
I understand what you mean but i don't understand in which way i could correct rho and to correct what error ? hyper-ventilation ?

EDIT : at the bow, the prismatic boundar layers are oriented ~ 30-40° of the flow, which don't help with the diffusivity. Refining the mesh seems to limit the problem but can't avoid this angle difference between flow and layers... by the way is it correct that the central differencing scheme, i.e, Gauss linear in OF, is the less diffusive scheme (in fact with no numerical diffusion) ?

Last edited by Quentin; July 15, 2014 at 04:44.

 July 15, 2014, 04:41 #10 Member   laurentL Join Date: Oct 2011 Location: new caledonia Posts: 72 Rep Power: 6 hi Quentin, did you manager to change the OF's code to do the trick, with the pressure? i also did a lots of tryings on better refinement on flat multihull but this trick was working perfectly well, once you know the pressure threshold (mostly depending of the speed ...) could you please help me to make the trick on OF? thanks LL

 July 15, 2014, 05:03 #11 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 Hi Laurent, yes I manage to do the trick, using the inletOutlet boundary condition as a start. I'll give you my code but i'm not very sure of what I've done and it's definitively not very clean, but it seems to work ; so be careful. It didn't wok in my case because I have both hyper ventilation of air under the hull and water at the bow ; so i'm working again with Dave's advices. You have to recompile the inletOutlet chaging all "inletOutlet" for "alphaHull", and modifying the content of 1) alphaHullfvPatchField.C Code: /*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . \*---------------------------------------------------------------------------*/ #include "alphaHullFvPatchField.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template Foam::alphaHullFvPatchField::alphaHullFvPatchField ( const fvPatch& p, const DimensionedField& iF ) : mixedFvPatchField(p, iF), phiName_("p"), pthreshold_(800.0) { this->refValue() = pTraits::zero; this->refGrad() = pTraits::zero; this->valueFraction() = 0.0; } template Foam::alphaHullFvPatchField::alphaHullFvPatchField ( const alphaHullFvPatchField& ptf, const fvPatch& p, const DimensionedField& iF, const fvPatchFieldMapper& mapper ) : mixedFvPatchField(ptf, p, iF, mapper), phiName_(ptf.phiName_) {} template Foam::alphaHullFvPatchField::alphaHullFvPatchField ( const fvPatch& p, const DimensionedField& iF, const dictionary& dict ) : mixedFvPatchField(p, iF), phiName_(dict.lookupOrDefault("phi", "p")), pthreshold_(dict.lookupOrDefault("pthreshold", 800.0)) { this->refValue() = Field("inletValue", dict, p.size()); fvPatchField::operator=(this->patchInternalField()); this->refGrad() = pTraits::zero; this->valueFraction() = 1.0; } template Foam::alphaHullFvPatchField::alphaHullFvPatchField ( const alphaHullFvPatchField& ptf ) : mixedFvPatchField(ptf), phiName_(ptf.phiName_) {} template Foam::alphaHullFvPatchField::alphaHullFvPatchField ( const alphaHullFvPatchField& ptf, const DimensionedField& iF ) : mixedFvPatchField(ptf, iF), phiName_(ptf.phiName_) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template void Foam::alphaHullFvPatchField::updateCoeffs() { if (this->updated()) { return; } const Field& phip = this->patch().template lookupPatchField ( "p" ); this->valueFraction() = pos(phip - pthreshold_); mixedFvPatchField::updateCoeffs(); } template void Foam::alphaHullFvPatchField::write(Ostream& os) const { fvPatchField::write(os); if (phiName_ != "phi") { os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl; } this->refValue().writeEntry("inletValue", os); this->writeEntry("value", os); } // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template void Foam::alphaHullFvPatchField::operator= ( const fvPatchField& ptf ) { fvPatchField::operator= ( this->valueFraction()*this->refValue() + (1 - this->valueFraction())*ptf ); } // ************************************************************************* // and 2)alphaHullFvPatchField.H Code: /*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . Class Foam::alphaHullFvPatchField Group grpOutletBoundaryConditions Description This boundary condition provides a generic outflow condition, with specified inflow for the case of return flow. \heading Patch usage \table Property | Description | Required | Default value phi | flux field name | no | phi alphaValue | inlet value for reverse flow | yes | \endtable Example of the boundary condition specification: \verbatim myPatch { type alphaHull; phi phi; alphaValue uniform 0; value uniform 0; } \endverbatim The mode of operation is determined by the sign of the flux across the patch faces. Note Sign conventions: - positive flux (out of domain): apply zero-gradient condition - negative flux (into of domain): apply the user-specified fixed value SeeAlso Foam::mixedFvPatchField Foam::zeroGradientFvPatchField SourceFiles alphaHullFvPatchField.C \*---------------------------------------------------------------------------*/ #ifndef alphaHullFvPatchField_H #define alphaHullFvPatchField_H #include "mixedFvPatchField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { /*---------------------------------------------------------------------------*\ Class alphaHullFvPatchField Declaration \*---------------------------------------------------------------------------*/ template class alphaHullFvPatchField : public mixedFvPatchField { protected: // Protected data //- Name of flux field word phiName_; //- Value of the pressure threshold scalar pthreshold_; public: //- Runtime type information TypeName("alphaHull"); // Constructors //- Construct from patch and internal field alphaHullFvPatchField ( const fvPatch&, const DimensionedField& ); //- Construct from patch, internal field and dictionary alphaHullFvPatchField ( const fvPatch&, const DimensionedField&, const dictionary& ); //- Construct by mappin"phi"g given alphaHullFvPatchField onto a new patch alphaHullFvPatchField ( const alphaHullFvPatchField&, const fvPatch&, const DimensionedField&, const fvPatchFieldMapper& ); //- Construct as copy alphaHullFvPatchField ( const alphaHullFvPatchField& ); //- Construct and return a clone virtual tmp > clone() const { return tmp > ( new alphaHullFvPatchField(*this) ); } //- Construct as copy setting internal field reference alphaHullFvPatchField ( const alphaHullFvPatchField&, const DimensionedField& ); //- Construct and return a clone setting internal field reference virtual tmp > clone ( const DimensionedField& iF ) const { return tmp > ( new alphaHullFvPatchField(*this, iF) ); } // Member functions //- Update the coefficients associated with the patch field virtual void updateCoeffs(); //- Write virtual void write(Ostream&) const; // Member operators virtual void operator=(const fvPatchField& pvf); }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #ifdef NoRepository # include "alphaHullFvPatchField.C" #endif // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif // ************************************************************************* // then recompile and add your library to controlDict dictionnary and the call of the BC would be : Code: hull { type alphaHull; inletValue uniform 1.0; } the important fonction is updatecoeff() in alphaHullFvPatchField.C. I hope it will help you but again my code is not clean at all, It's at test phase.

 July 16, 2014, 03:53 Adding source term to avoid numerical ventilation #12 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 Hi foamers, Dave, With the knowledge of the name of the problem, i found his very intesting publication about resistance calculation for AC33 hulls. In paragraph 4.4 they are speaking about numerical ventilation, explaining how you can correct your skin resistance calculation with rho, but also telling that you could solve the problem adding a source term in the transport equations of ad \alpha_{air} wich would be something like with being an arbitrary distance from the wall and an arbitrary number of iterations: in the transport equation of the air phase : if && , then , else in the transport equation of the water phase : if && , then , else As the interFoam solver only solve one transport equation for alpha, i.e : with ; I think I could rewrite the source terme as : if && , then , else if && , then , else . but i'm no quite sure of this as if i write the transport equation with a source term for : this give me if i transform into the transport equation for : so a gradient of U appears. What do you think about it ? EDIT : it's not a gradient, it's a divergence ... Last edited by Quentin; July 17, 2014 at 04:39.

 July 18, 2014, 12:48 #13 Member   Dave Join Date: Jul 2010 Posts: 97 Rep Power: 8 Dear Quentin, Yes, this paper was what I had in mind though again it is only good if you have wall functions. Regarding source terms, I would suggest looking at the cavitatingFoam and/or interPhaseChangeFoam to see how they have handled sources in the alpha equation.

 July 30, 2014, 05:12 Solved #14 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 Dear Dave, Laurent & foamers, I managed to imlement this source term in LTSinterFoam solver. 1) the source I finally implemented is : if && if && else 2) I implemented implicitely he negative part and explicitely the positive one, so the disctretised equation has the form (with Euler explicit time discretisation) : if && if && usual form else. I choose usually n = 200. 3) so Here is the new alphaEqn.H ; please note that i'm not completly sure of what i did in MULS::LTScorrect and MULES::LTSexplicitsolve : Code: { word alphaScheme("div(phi,alpha)"); word alpharScheme("div(phirb,alpha)"); // Standard face-flux compression coefficient surfaceScalarField phic(interface.cAlpha()*mag(phi/mesh.magSf())); // Add the optional isotropic compression contribution if (icAlpha > 0) { phic *= (1.0 - icAlpha); phic += (interface.cAlpha()*icAlpha)*fvc::interpolate(mag(U)); } // Do not compress interface at non-coupled boundary faces // (inlets, outlets etc.) forAll(phic.boundaryField(), patchi) { fvsPatchScalarField& phicp = phic.boundaryField()[patchi]; if (!phicp.coupled()) { phicp == 0; } } //** Creating Source Terms volScalarField::DimensionedInternalField Su ( IOobject ( "Su", runTime.timeName(), mesh ), mesh, dimensionedScalar("Su", dimensionSet(0,0, -1,0,0,0,0), 0.0) ); volScalarField::DimensionedField Sp ( IOobject ( "Sp", runTime.timeName(), mesh ), mesh, dimensionedScalar("Sp", dimensionSet(0,0, -1,0,0,0,0), 0.0) ); //** Evaluating Source Terms forAll (alpha1, celli) { if(y[celli] <= d && alpha1[celli] <= 0.5 && alpha1[celli] >= 0.0) { Sp[celli] -= 1.0/n * rsubDeltaT.field()[celli] ; } else if(y[celli] <= d && alpha1[celli] >= 0.5 && alpha1[celli] <= 1.0) { Su[celli] += 1.0/n * rsubDeltaT.field()[celli] ; Sp[celli] -= 1.0/n * rsubDeltaT.field()[celli] ; } } tmp tphiAlpha; if (MULESCorr) { //** Modified Alpha Eqn fvScalarMatrix alpha1Eqn ( #ifdef LTSSOLVE fv::localEulerDdtScheme(mesh, rDeltaT.name()).fvmDdt(alpha1) #else fv::EulerDdtScheme(mesh).fvmDdt(alpha1) #endif + fv::gaussConvectionScheme ( mesh, phi, upwind(mesh, phi) ).fvmDiv(phi, alpha1) ); solve(alpha1Eqn == Su + fvm::Sp(Sp,alpha1)); Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(alpha1) = " << min(alpha1).value() << " Max(alpha1) = " << max(alpha1).value() << endl; tmp tphiAlphaUD(alpha1Eqn.flux()); tphiAlpha = tmp ( new surfaceScalarField(tphiAlphaUD()) ); if (alphaApplyPrevCorr && tphiAlphaCorr0.valid()) { //** evaluation of alpha_old for LTScorrect volScalarField alpha10(alpha1); Info<< "Applying the previous iteration compression flux" << endl; #ifdef LTSSOLVE //** modified LTScorrect - /!\ unsure MULES::LTScorrect ( geometricOneField(), alpha1, tphiAlpha(), tphiAlphaCorr0(), Sp, (-Sp*alpha10)(), 1, 0 ); #else MULES::correct(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0); #endif tphiAlpha() += tphiAlphaCorr0(); } // Cache the upwind-flux tphiAlphaCorr0 = tphiAlphaUD; alpha2 = 1.0 - alpha1; interface.correct(); } for (int aCorr=0; aCorr tphiAlphaUn ( fvc::flux ( phi, alpha1, alphaScheme ) + fvc::flux ( -fvc::flux(-phir, alpha2, alpharScheme), alpha1, alpharScheme ) ); if (MULESCorr) { tmp tphiAlphaCorr(tphiAlphaUn() - tphiAlpha()); volScalarField alpha10(alpha1); #ifdef LTSSOLVE //** modified LTScorrect - /!\ unsure MULES::LTScorrect( geometricOneField(), alpha1, tphiAlphaUn(), tphiAlphaCorr(), Sp, (-Sp*alpha10)(), 1, 0); #else MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0); #endif // Under-relax the correction for all but the 1st corrector if (aCorr == 0) { tphiAlpha() += tphiAlphaCorr(); } else { alpha1 = 0.5*alpha1 + 0.5*alpha10; tphiAlpha() += 0.5*tphiAlphaCorr(); } } else { tphiAlpha = tphiAlphaUn; #ifdef LTSSOLVE //** modified explicitLTSSolve - /!\ unsure MULES::explicitLTSSolve(geometricOneField(), alpha1, phi, tphiAlpha(), Sp, Su, 1, 0); #else MULES::explicitSolve(alpha1, phi, tphiAlpha(), 1, 0); #endif } alpha2 = 1.0 - alpha1; interface.correct(); } rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2; if (alphaApplyPrevCorr && MULESCorr) { tphiAlphaCorr0 = tphiAlpha() - tphiAlphaCorr0; } Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(alpha1) = " << min(alpha1).value() << " Max(alpha1) = " << max(alpha1).value() << endl; } I also added in createField.H HTML Code:  //** Wall distance Info<< "Calculating wall distance field" <("d", 0.0) ); scalar n ( alphaControls.lookupOrDefault("n", 1.0) ); so the Source terms n and d can be set in fvSolution in alpha.water solver's parameters. I havn't tested it with nalphasubcycles > 1, but otherwise, it seems to work for me. Last edited by Quentin; July 30, 2014 at 07:57.

December 10, 2014, 07:54
#15
Member

james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 37
Rep Power: 3
Quote:
 Originally Posted by Quentin Dear Dave, Laurent & foamers, I managed to imlement this source term in LTSinterFoam solver. 1) the source I finally implemented is : if && if && else 2) I implemented implicitely he negative part and explicitely the positive one, so the disctretised equation has the form (with Euler explicit time discretisation) : if && if && usual form else. I choose usually n = 200. 3) so Here is the new alphaEqn.H ; please note that i'm not completly sure of what i did in MULS::LTScorrect and MULES::LTSexplicitsolve : Code: { word alphaScheme("div(phi,alpha)"); word alpharScheme("div(phirb,alpha)"); // Standard face-flux compression coefficient surfaceScalarField phic(interface.cAlpha()*mag(phi/mesh.magSf())); // Add the optional isotropic compression contribution if (icAlpha > 0) { phic *= (1.0 - icAlpha); phic += (interface.cAlpha()*icAlpha)*fvc::interpolate(mag(U)); } // Do not compress interface at non-coupled boundary faces // (inlets, outlets etc.) forAll(phic.boundaryField(), patchi) { fvsPatchScalarField& phicp = phic.boundaryField()[patchi]; if (!phicp.coupled()) { phicp == 0; } } //** Creating Source Terms volScalarField::DimensionedInternalField Su ( IOobject ( "Su", runTime.timeName(), mesh ), mesh, dimensionedScalar("Su", dimensionSet(0,0, -1,0,0,0,0), 0.0) ); volScalarField::DimensionedField Sp ( IOobject ( "Sp", runTime.timeName(), mesh ), mesh, dimensionedScalar("Sp", dimensionSet(0,0, -1,0,0,0,0), 0.0) ); //** Evaluating Source Terms forAll (alpha1, celli) { if(y[celli] <= d && alpha1[celli] <= 0.5 && alpha1[celli] >= 0.0) { Sp[celli] -= 1.0/n * rsubDeltaT.field()[celli] ; } else if(y[celli] <= d && alpha1[celli] >= 0.5 && alpha1[celli] <= 1.0) { Su[celli] += 1.0/n * rsubDeltaT.field()[celli] ; Sp[celli] -= 1.0/n * rsubDeltaT.field()[celli] ; } } tmp tphiAlpha; if (MULESCorr) { //** Modified Alpha Eqn fvScalarMatrix alpha1Eqn ( #ifdef LTSSOLVE fv::localEulerDdtScheme(mesh, rDeltaT.name()).fvmDdt(alpha1) #else fv::EulerDdtScheme(mesh).fvmDdt(alpha1) #endif + fv::gaussConvectionScheme ( mesh, phi, upwind(mesh, phi) ).fvmDiv(phi, alpha1) ); solve(alpha1Eqn == Su + fvm::Sp(Sp,alpha1)); Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(alpha1) = " << min(alpha1).value() << " Max(alpha1) = " << max(alpha1).value() << endl; tmp tphiAlphaUD(alpha1Eqn.flux()); tphiAlpha = tmp ( new surfaceScalarField(tphiAlphaUD()) ); if (alphaApplyPrevCorr && tphiAlphaCorr0.valid()) { //** evaluation of alpha_old for LTScorrect volScalarField alpha10(alpha1); Info<< "Applying the previous iteration compression flux" << endl; #ifdef LTSSOLVE //** modified LTScorrect - /!\ unsure MULES::LTScorrect ( geometricOneField(), alpha1, tphiAlpha(), tphiAlphaCorr0(), Sp, (-Sp*alpha10)(), 1, 0 ); #else MULES::correct(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0); #endif tphiAlpha() += tphiAlphaCorr0(); } // Cache the upwind-flux tphiAlphaCorr0 = tphiAlphaUD; alpha2 = 1.0 - alpha1; interface.correct(); } for (int aCorr=0; aCorr tphiAlphaUn ( fvc::flux ( phi, alpha1, alphaScheme ) + fvc::flux ( -fvc::flux(-phir, alpha2, alpharScheme), alpha1, alpharScheme ) ); if (MULESCorr) { tmp tphiAlphaCorr(tphiAlphaUn() - tphiAlpha()); volScalarField alpha10(alpha1); #ifdef LTSSOLVE //** modified LTScorrect - /!\ unsure MULES::LTScorrect( geometricOneField(), alpha1, tphiAlphaUn(), tphiAlphaCorr(), Sp, (-Sp*alpha10)(), 1, 0); #else MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0); #endif // Under-relax the correction for all but the 1st corrector if (aCorr == 0) { tphiAlpha() += tphiAlphaCorr(); } else { alpha1 = 0.5*alpha1 + 0.5*alpha10; tphiAlpha() += 0.5*tphiAlphaCorr(); } } else { tphiAlpha = tphiAlphaUn; #ifdef LTSSOLVE //** modified explicitLTSSolve - /!\ unsure MULES::explicitLTSSolve(geometricOneField(), alpha1, phi, tphiAlpha(), Sp, Su, 1, 0); #else MULES::explicitSolve(alpha1, phi, tphiAlpha(), 1, 0); #endif } alpha2 = 1.0 - alpha1; interface.correct(); } rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2; if (alphaApplyPrevCorr && MULESCorr) { tphiAlphaCorr0 = tphiAlpha() - tphiAlphaCorr0; } Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(alpha1) = " << min(alpha1).value() << " Max(alpha1) = " << max(alpha1).value() << endl; } I also added in createField.H HTML Code:  //** Wall distance Info<< "Calculating wall distance field" <("d", 0.0) ); scalar n ( alphaControls.lookupOrDefault("n", 1.0) ); so the Source terms n and d can be set in fvSolution in alpha.water solver's parameters. I havn't tested it with nalphasubcycles > 1, but otherwise, it seems to work for me.

Quentin,

Thanks for the info. With little knowledge of OpenFOAM and the discretization of these equations, Implementing source terms is quite the challenge. I have a few questions.

Using these conditions:

"
1) the source I finally implemented is :
if &&
if &&
else
"

I see this code:

Code:
 //** Evaluating Source Terms
forAll (alpha1, celli)
{
if(y[celli] <= d && alpha1[celli] <= 0.5 && alpha1[celli] >= 0.0)
{

Sp[celli] -=  1.0/n * rsubDeltaT.field()[celli] ;
}
else if(y[celli] <= d && alpha1[celli] >= 0.5 && alpha1[celli] <= 1.0)
{
Su[celli] += 1.0/n * rsubDeltaT.field()[celli] ;

Sp[celli] -=  1.0/n * rsubDeltaT.field()[celli] ;
}
}
What is rsubDeltaT.field()? can you explain the input and output of this field, how its defined or what values and units it takes?

Could you explain for both conditions above (alpha <.5 and alpha > .5) why Sp is defined in both conditions and Su only defined in the latter?

Thank you so much for your help, James

 January 10, 2015, 10:49 #16 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 Hi James, I havn't work on this for a few month so I don't remember exactly but : Su is the part of the source term not linked to [MATH]\alpha[MATH\] and Sp is the one linked to. It allows you to implement the source term explicitely or implicitely (Su => explicit ; Sp => implicit). there are some conditions regarding the stability of the solution whether the implicit part is positiv or negativ ( I don't remember exactly). I would add that using this source term seems very dangerous to me because it brakes the physic and it was much more correct to activate the compression term already implemented in OpenFOAM through the parameter calpha. in my case I finally used for the alpha solver : calpha = 1; nalphacorr = 0 for the first iterations (before the pression equation was stabilized (~ 50-100iterations) ) then 1 nalphasubcycle = 1 and i wasn't using a PIMPLE loop (noutercorrector = 1) hopfully it will help you, Quentin

 January 10, 2015, 10:56 #17 New Member   Quentin Coispeau Join Date: Nov 2013 Posts: 17 Rep Power: 4 rsubdeltaT is the field containing the inverse of local time step (you have rdeltaT which is the inverse of main local time step and r subdeltaT is rDeltaT times bnalphasubcycle). Regarding Su and Sp, in the "> 0.5" condition, the first part is implemented explicitely (Su) and the second implicitely (Sp). For the "<0.5" one, there is only an implicit part Sp and Su = 0.

 Tags alpha1, ltsinterfoam

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

 Similar Threads Thread Thread Starter Forum Replies Last Post Wokl OpenFOAM Running, Solving & CFD 4 January 9, 2014 03:52 JFM OpenFOAM Running, Solving & CFD 0 November 17, 2013 01:15 Kopy OpenFOAM Running, Solving & CFD 3 March 15, 2013 12:27 sharonyue OpenFOAM Running, Solving & CFD 13 January 2, 2013 23:40 RedAdmiral Main CFD Forum 0 November 18, 2012 10:51

All times are GMT -4. The time now is 06:33.