research/OpenFOAM-5.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

lagrangian::ThermoSurfaceFilm: Updated adhesion->splash transition Weber number

Browse source 
according to

        Bai et al, `Modelling of gasoline spray impingement', Atom. Sprays,
        vol 12, pp 1-27, 2002

Resolves bug-report https://bugs.openfoam.org/view.php?id=2478
This commit is contained in:
Henry Weller 2017-03-24 14:56:47 +00:00
parent 86bf44de10
commit 78e40da8ac
2 changed files with 36 additions and 36 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

68
src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
View file

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -111,13 +111,13 @@ Foam::vector Foam::ThermoSurfaceFilm<CloudType>::splashDirection
const vector& nf const vector& nf
) const ) const
{ {
// azimuthal angle [rad] // Azimuthal angle [rad]
const scalar phiSi = twoPi*rndGen_.sample01<scalar>(); const scalar phiSi = twoPi*rndGen_.sample01<scalar>();
// ejection angle [rad] // Ejection angle [rad]
const scalar thetaSi = pi/180.0*(rndGen_.sample01<scalar>()*(50 - 5) + 5); const scalar thetaSi = pi/180.0*(rndGen_.sample01<scalar>()*(50 - 5) + 5);
// direction vector of new parcel // Direction vector of new parcel
const scalar alpha = sin(thetaSi); const scalar alpha = sin(thetaSi);
const scalar dcorr = cos(thetaSi); const scalar dcorr = cos(thetaSi);
const vector normal = alpha*(tanVec1*cos(phiSi) + tanVec2*sin(phiSi)); const vector normal = alpha*(tanVec1*cos(phiSi) + tanVec2*sin(phiSi));
@ -226,7 +226,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::drySplashInteraction
const vector& Up = this->owner().U().boundaryField()[pp.index()][facei]; const vector& Up = this->owner().U().boundaryField()[pp.index()][facei];
const vector& nf = pp.faceNormals()[facei]; const vector& nf = pp.faceNormals()[facei];
// local pressure // Local pressure
const scalar pc = thermo_.thermo().p()[p.cell()]; const scalar pc = thermo_.thermo().p()[p.cell()];
// Retrieve parcel properties // Retrieve parcel properties
@ -247,13 +247,13 @@ void Foam::ThermoSurfaceFilm<CloudType>::drySplashInteraction
// Critical Weber number // Critical Weber number
const scalar Wec = Adry_*pow(La, -0.183); const scalar Wec = Adry_*pow(La, -0.183);
if (We < Wec) // adhesion - assume absorb if (We < Wec) // Adhesion - assume absorb
{ {
absorbInteraction(filmModel, p, pp, facei, m, keepParticle); absorbInteraction(filmModel, p, pp, facei, m, keepParticle);
} }
else // splash else // Splash
{ {
// ratio of incident mass to splashing mass // Ratio of incident mass to splashing mass
const scalar mRatio = 0.2 + 0.6*rndGen_.sample01<scalar>(); const scalar mRatio = 0.2 + 0.6*rndGen_.sample01<scalar>();
splashInteraction splashInteraction
(filmModel, p, pp, facei, mRatio, We, Wec, sigma, keepParticle); (filmModel, p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
@ -282,7 +282,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
const vector& Up = this->owner().U().boundaryField()[pp.index()][facei]; const vector& Up = this->owner().U().boundaryField()[pp.index()][facei];
const vector& nf = pp.faceNormals()[facei]; const vector& nf = pp.faceNormals()[facei];
// local pressure // Local pressure
const scalar pc = thermo_.thermo().p()[p.cell()]; const scalar pc = thermo_.thermo().p()[p.cell()];
// Retrieve parcel properties // Retrieve parcel properties
@ -305,31 +305,31 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
// Critical Weber number // Critical Weber number
const scalar Wec = Awet_*pow(La, -0.183); const scalar Wec = Awet_*pow(La, -0.183);
if (We < 1) // adhesion - assume absorb if (We < 2) // Adhesion - assume absorb
{ {
absorbInteraction(filmModel, p, pp, facei, m, keepParticle); absorbInteraction(filmModel, p, pp, facei, m, keepParticle);
} }
else if ((We >= 1) && (We < 20)) // bounce else if ((We >= 2) && (We < 20)) // Bounce
{ {
// incident angle of impingement // Incident angle of impingement
const scalar theta = pi/2 - acos(U/mag(U) & nf); const scalar theta = pi/2 - acos(U/mag(U) & nf);
// restitution coefficient // Restitution coefficient
const scalar epsilon = 0.993 - theta*(1.76 - theta*(1.56 - theta*0.49)); const scalar epsilon = 0.993 - theta*(1.76 - theta*(1.56 - theta*0.49));
// update parcel velocity // Update parcel velocity
U = -epsilon*(Un) + 5/7*(Ut); U = -epsilon*(Un) + 5/7*(Ut);
keepParticle = true; keepParticle = true;
return; return;
} }
else if ((We >= 20) && (We < Wec)) // spread - assume absorb else if ((We >= 20) && (We < Wec)) // Spread - assume absorb
{ {
absorbInteraction(filmModel, p, pp, facei, m, keepParticle); absorbInteraction(filmModel, p, pp, facei, m, keepParticle);
} }
else // splash else // Splash
{ {
// ratio of incident mass to splashing mass // Ratio of incident mass to splashing mass
// splash mass can be > incident mass due to film entrainment // splash mass can be > incident mass due to film entrainment
const scalar mRatio = 0.2 + 0.9*rndGen_.sample01<scalar>(); const scalar mRatio = 0.2 + 0.9*rndGen_.sample01<scalar>();
splashInteraction splashInteraction
@ -371,24 +371,24 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
const vector& posC = mesh.C()[p.cell()]; const vector& posC = mesh.C()[p.cell()];
const vector& posCf = mesh.Cf().boundaryField()[pp.index()][facei]; const vector& posCf = mesh.Cf().boundaryField()[pp.index()][facei];
// total mass of (all) splashed parcels // Total mass of (all) splashed parcels
const scalar mSplash = m*mRatio; const scalar mSplash = m*mRatio;
// number of splashed particles per incoming particle // Number of splashed particles per incoming particle
const scalar Ns = 5.0*(We/Wec - 1.0); const scalar Ns = 5.0*(We/Wec - 1.0);
// average diameter of splashed particles // Average diameter of splashed particles
const scalar dBarSplash = 1/cbrt(6.0)*cbrt(mRatio/Ns)*d + ROOTVSMALL; const scalar dBarSplash = 1/cbrt(6.0)*cbrt(mRatio/Ns)*d + ROOTVSMALL;
// cumulative diameter splash distribution // Cumulative diameter splash distribution
const scalar dMax = 0.9*cbrt(mRatio)*d; const scalar dMax = 0.9*cbrt(mRatio)*d;
const scalar dMin = 0.1*dMax; const scalar dMin = 0.1*dMax;
const scalar K = exp(-dMin/dBarSplash) - exp(-dMax/dBarSplash); const scalar K = exp(-dMin/dBarSplash) - exp(-dMax/dBarSplash);
// surface energy of secondary parcels [J] // Surface energy of secondary parcels [J]
scalar ESigmaSec = 0; scalar ESigmaSec = 0;
// sample splash distribution to determine secondary parcel diameters // Sample splash distribution to determine secondary parcel diameters
scalarList dNew(parcelsPerSplash_); scalarList dNew(parcelsPerSplash_);
scalarList npNew(parcelsPerSplash_); scalarList npNew(parcelsPerSplash_);
forAll(dNew, i) forAll(dNew, i)
@ -399,26 +399,26 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
ESigmaSec += npNew[i]*sigma*p.areaS(dNew[i]); ESigmaSec += npNew[i]*sigma*p.areaS(dNew[i]);
} }
// incident kinetic energy [J] // Incident kinetic energy [J]
const scalar EKIn = 0.5*m*magSqr(Urel); const scalar EKIn = 0.5*m*magSqr(Urel);
// incident surface energy [J] // Incident surface energy [J]
const scalar ESigmaIn = np*sigma*p.areaS(d); const scalar ESigmaIn = np*sigma*p.areaS(d);
// dissipative energy // Dissipative energy
const scalar Ed = max(0.8*EKIn, np*Wec/12*pi*sigma*sqr(d)); const scalar Ed = max(0.8*EKIn, np*Wec/12*pi*sigma*sqr(d));
// total energy [J] // Total energy [J]
const scalar EKs = EKIn + ESigmaIn - ESigmaSec - Ed; const scalar EKs = EKIn + ESigmaIn - ESigmaSec - Ed;
// switch to absorb if insufficient energy for splash // Switch to absorb if insufficient energy for splash
if (EKs <= 0) if (EKs <= 0)
{ {
absorbInteraction(filmModel, p, pp, facei, m, keepParticle); absorbInteraction(filmModel, p, pp, facei, m, keepParticle);
return; return;
} }
// helper variables to calculate magUns0 // Helper variables to calculate magUns0
const scalar logD = log(d); const scalar logD = log(d);
const scalar coeff2 = log(dNew[0]) - logD + ROOTVSMALL; const scalar coeff2 = log(dNew[0]) - logD + ROOTVSMALL;
scalar coeff1 = 0.0; scalar coeff1 = 0.0;
@ -427,7 +427,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
coeff1 += sqr(log(dNew[i]) - logD); coeff1 += sqr(log(dNew[i]) - logD);
} }
// magnitude of the normal velocity of the first splashed parcel // Magnitude of the normal velocity of the first splashed parcel
const scalar magUns0 = const scalar magUns0 =
sqrt(2.0*parcelsPerSplash_*EKs/mSplash/(1.0 + coeff1/sqr(coeff2))); sqrt(2.0*parcelsPerSplash_*EKs/mSplash/(1.0 + coeff1/sqr(coeff2)));
@ -448,7 +448,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
pPtr->typeId() = splashParcelType_; pPtr->typeId() = splashParcelType_;
} }
// perturb new parcels towards the owner cell centre // Perturb new parcels towards the owner cell centre
pPtr->position() += 0.5*rndGen_.sample01<scalar>()*(posC - posCf); pPtr->position() += 0.5*rndGen_.sample01<scalar>()*(posC - posCf);
pPtr->nParticle() = npNew[i]; pPtr->nParticle() = npNew[i];
@ -466,7 +466,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
nParcelsSplashed_++; nParcelsSplashed_++;
} }
// transfer remaining part of parcel to film 0 - splashMass can be -ve // Transfer remaining part of parcel to film 0 - splashMass can be -ve
// if entraining from the film // if entraining from the film
const scalar mDash = m - mSplash; const scalar mDash = m - mSplash;
absorbInteraction(filmModel, p, pp, facei, mDash, keepParticle); absorbInteraction(filmModel, p, pp, facei, mDash, keepParticle);
@ -613,11 +613,11 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
} }
} }
// transfer parcel/parcel interactions complete // Transfer parcel/parcel interactions complete
return true; return true;
} }
// parcel not interacting with film // Parcel not interacting with film
return false; return false;
} }

4
src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
View file

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -36,7 +36,7 @@ Description
Bai and Gosman, `Mathematical modelling of wall films formed by Bai and Gosman, `Mathematical modelling of wall films formed by
impinging sprays', SAE 960626, 1996 impinging sprays', SAE 960626, 1996
Bai et al, `Modelling off gasoline spray impingement', Atom. Sprays, Bai et al, `Modelling of gasoline spray impingement', Atom. Sprays,
vol 12, pp 1-27, 2002 vol 12, pp 1-27, 2002