
[Sponsors] 
July 2, 2014, 08:36 
alpha1 wrong behaviour with LTSInterFoam [solved]

#1 
New Member
Quentin Coispeau
Join Date: Nov 2013
Posts: 17
Rep Power: 12 
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 setup 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 "cylinderhull" illustrating this, wich i plotted the interface, colored by z ; and the value of alpha.water on the "cylinderhull". 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 "cylinderhull" there is a waterair 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 45 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 1e8; relTol 0; minIter 1; } pcorr { solver PCG; preconditioner { preconditioner GAMG; smoother GaussSeidel; agglomerator faceAreaPair; mergeLevels 1; nCellsInCoarsestLevel 100; cacheAgglomeration true; tolerance 1e6; relTol 0; }; tolerance 1e6; relTol 0; }; p_rgh { solver GAMG; smoother DIC; agglomerator faceAreaPair; mergeLevels 1; nCellsInCoarsestLevel 100; cacheAgglomeration true; tolerance 1e8; relTol 0.0; }; p_rghFinal { $p_rgh; relTol 0; } "(Ukomega).*" { solver smoothSolver; smoother symGaussSeidel; nSweeps 1; tolerance 1e7; 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 { "(Ukomega).*" 0.7; "alpha.water.*" 0.7; } } 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; } thank you in advance for your help, Quentin Coispeau ps : in my cylindertest 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 "magnumlike" 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: 12 
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: 73
Rep Power: 14 
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... http://www.cfdonline.com/Forums/ope...pressure.html Laurent 

July 3, 2014, 11:41 

#4 
New Member
Quentin Coispeau
Join Date: Nov 2013
Posts: 17
Rep Power: 12 
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 1e8; 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 1e5; relTol 0; }; tolerance 1e5; 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 1e5; relTol 0.05; }; p_rghFinal { solver PCG; preconditioner { preconditioner GAMG; tolerance 1e07; relTol 0; nVcycles 2; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; cacheAgglomeration on; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; } tolerance 1e07; relTol 0; maxIter 20; } "(Ukomega).*" { solver smoothSolver; smoother symGaussSeidel; nSweeps 1; tolerance 1e7; 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 { "(Ukomega).*" 0.7; "alpha.water.*" 0.7;*/ } } cache { grad(U); } 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; } 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: 12 
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: 73
Rep Power: 14 
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 2add a line in /home/laurent/OpenFOAM/OpenFOAM2.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: 12 
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 

July 13, 2014, 12:08 

#8 
Senior Member
Dave
Join Date: Jul 2010
Posts: 100
Rep Power: 16 
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: 12 
Hi Dave,
thank you very much for your answer. knowing the name of the issue will help me find documentation about it. I had already tried to have less diffusive numerical scheme but it wasn't enough. I'll follow your advice about refining the mesh where the problem appears. I agree with you saying that the pressure threshold trick is more a bandage than a solution. Quote:
EDIT : at the bow, the prismatic boundar layers are oriented ~ 3040° 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: 73
Rep Power: 14 
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: 12 
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 <http://www.gnu.org/licenses/>. \**/ #include "alphaHullFvPatchField.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template<class Type> Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField ( const fvPatch& p, const DimensionedField<Type, volMesh>& iF ) : mixedFvPatchField<Type>(p, iF), phiName_("p"), pthreshold_(800.0) { this>refValue() = pTraits<Type>::zero; this>refGrad() = pTraits<Type>::zero; this>valueFraction() = 0.0; } template<class Type> Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField ( const alphaHullFvPatchField<Type>& ptf, const fvPatch& p, const DimensionedField<Type, volMesh>& iF, const fvPatchFieldMapper& mapper ) : mixedFvPatchField<Type>(ptf, p, iF, mapper), phiName_(ptf.phiName_) {} template<class Type> Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField ( const fvPatch& p, const DimensionedField<Type, volMesh>& iF, const dictionary& dict ) : mixedFvPatchField<Type>(p, iF), phiName_(dict.lookupOrDefault<word>("phi", "p")), pthreshold_(dict.lookupOrDefault<scalar>("pthreshold", 800.0)) { this>refValue() = Field<Type>("inletValue", dict, p.size()); fvPatchField<Type>::operator=(this>patchInternalField()); this>refGrad() = pTraits<Type>::zero; this>valueFraction() = 1.0; } template<class Type> Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField ( const alphaHullFvPatchField<Type>& ptf ) : mixedFvPatchField<Type>(ptf), phiName_(ptf.phiName_) {} template<class Type> Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField ( const alphaHullFvPatchField<Type>& ptf, const DimensionedField<Type, volMesh>& iF ) : mixedFvPatchField<Type>(ptf, iF), phiName_(ptf.phiName_) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class Type> void Foam::alphaHullFvPatchField<Type>::updateCoeffs() { if (this>updated()) { return; } const Field<scalar>& phip = this>patch().template lookupPatchField<volScalarField, scalar> ( "p" ); this>valueFraction() = pos(phip  pthreshold_); mixedFvPatchField<Type>::updateCoeffs(); } template<class Type> void Foam::alphaHullFvPatchField<Type>::write(Ostream& os) const { fvPatchField<Type>::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<class Type> void Foam::alphaHullFvPatchField<Type>::operator= ( const fvPatchField<Type>& ptf ) { fvPatchField<Type>::operator= ( this>valueFraction()*this>refValue() + (1  this>valueFraction())*ptf ); } // ************************************************************************* // 2)alphaHullFvPatchField.H Code:
/**\ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  \\ / A nd  Copyright (C) 20112012 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 <http://www.gnu.org/licenses/>. 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 zerogradient condition  negative flux (into of domain): apply the userspecified 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 Type> class alphaHullFvPatchField : public mixedFvPatchField<Type> { 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<Type, volMesh>& ); // Construct from patch, internal field and dictionary alphaHullFvPatchField ( const fvPatch&, const DimensionedField<Type, volMesh>&, const dictionary& ); // Construct by mappin"phi"g given alphaHullFvPatchField onto a new patch alphaHullFvPatchField ( const alphaHullFvPatchField<Type>&, const fvPatch&, const DimensionedField<Type, volMesh>&, const fvPatchFieldMapper& ); // Construct as copy alphaHullFvPatchField ( const alphaHullFvPatchField<Type>& ); // Construct and return a clone virtual tmp<fvPatchField<Type> > clone() const { return tmp<fvPatchField<Type> > ( new alphaHullFvPatchField<Type>(*this) ); } // Construct as copy setting internal field reference alphaHullFvPatchField ( const alphaHullFvPatchField<Type>&, const DimensionedField<Type, volMesh>& ); // Construct and return a clone setting internal field reference virtual tmp<fvPatchField<Type> > clone ( const DimensionedField<Type, volMesh>& iF ) const { return tmp<fvPatchField<Type> > ( new alphaHullFvPatchField<Type>(*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<Type>& pvf); }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #ifdef NoRepository # include "alphaHullFvPatchField.C" #endif // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif // ************************************************************************* // Code:
hull { type alphaHull; inletValue uniform 1.0; } 

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: 12 
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:
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 
Senior Member
Dave
Join Date: Jul 2010
Posts: 100
Rep Power: 16 
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: 12 
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 faceflux 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 noncoupled 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<surfaceScalarField> tphiAlpha; if (MULESCorr) { //** Modified Alpha Eqn fvScalarMatrix alpha1Eqn ( #ifdef LTSSOLVE fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1) #else fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1) #endif + fv::gaussConvectionScheme<scalar> ( mesh, phi, upwind<scalar>(mesh, phi) ).fvmDiv(phi, alpha1) ); solve(alpha1Eqn == Su + fvm::Sp(Sp,alpha1)); Info<< "Phase1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(alpha1) = " << min(alpha1).value() << " Max(alpha1) = " << max(alpha1).value() << endl; tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux()); tphiAlpha = tmp<surfaceScalarField> ( 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 upwindflux tphiAlphaCorr0 = tphiAlphaUD; alpha2 = 1.0  alpha1; interface.correct(); } for (int aCorr=0; aCorr<nAlphaCorr; aCorr++) { surfaceScalarField phir(phic*interface.nHatf()); tmp<surfaceScalarField> tphiAlphaUn ( fvc::flux ( phi, alpha1, alphaScheme ) + fvc::flux ( fvc::flux(phir, alpha2, alpharScheme), alpha1, alpharScheme ) ); if (MULESCorr) { tmp<surfaceScalarField> 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 // Underrelax 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<< "Phase1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(alpha1) = " << min(alpha1).value() << " Max(alpha1) = " << max(alpha1).value() << endl; } HTML Code:
//** Wall distance Info<< "Calculating wall distance field" <<endl; volScalarField y = wallDist(mesh).y(); Code:
//** Added Source Terms controls scalar d ( alphaControls.lookupOrDefault<scalar>("d", 0.0) ); scalar n ( alphaControls.lookupOrDefault<scalar>("n", 1.0) ); 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, 06:54 

#15 
Member
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 39
Rep Power: 12 
Quote:
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] ; } } 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, 09:49 

#16 
New Member
Quentin Coispeau
Join Date: Nov 2013
Posts: 17
Rep Power: 12 
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 (~ 50100iterations) ) then 1 nalphasubcycle = 1 and i wasn't using a PIMPLE loop (noutercorrector = 1) hopfully it will help you, Quentin 

January 10, 2015, 09:56 

#17 
New Member
Quentin Coispeau
Join Date: Nov 2013
Posts: 17
Rep Power: 12 
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  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
InterFoam channel Flow runnig very slowly  Wokl  OpenFOAM Running, Solving & CFD  4  January 9, 2014 02:52 
LTSInterFoam Crashing  JFM  OpenFOAM Running, Solving & CFD  0  November 17, 2013 00:15 
Unexpected behaviour of timeVaryingMappedFixedValue  Kopy  OpenFOAM Running, Solving & CFD  3  March 15, 2013 11:27 
same geometry,structured and unstructured mesh,different behaviour.  sharonyue  OpenFOAM Running, Solving & CFD  13  January 2, 2013 22:40 
Particle behaviour in cfd  RedAdmiral  Main CFD Forum  0  November 18, 2012 09:51 