View Issue Details
ID | Project | Category | View Status | Date Submitted | Last Update |
---|---|---|---|---|---|
0002647 | OpenFOAM | Bug | public | 2017-08-03 10:43 | 2018-07-10 11:25 |
Reporter | Arne | Assigned To | henry | ||
Priority | normal | Severity | minor | Reproducibility | always |
Status | closed | Resolution | no change required | ||
Platform | GNU/Linux | OS | Ubuntu | OS Version | 14.04 |
Product Version | dev | ||||
Summary | 0002647: frictional stress in kinetic theory modeling not taken into account | ||||
Description | Based on my understanding of the thesis of Berend van Wachem the frictional stress has to be taken into account for computing granular temperature equation. In twoPhaseEulerFoam only the granular part is taken into account. To be more specific the file in question is: $FOAM_SOLVERS/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel.C In line 388: the solid shear viscosity is computed based the table in van Wachems thesis (p.47 table 3.2). After that the viscosity is used to compute the solid stress tensor. However, in the thesis is stated (p.52 equation 3.27 and 3.28) that if the solid volume fraction is higher than a threshold value the frictional values have to be added. To my understanding this doesn't happen in the current implementation of the granular temperature equation. The same issue occurs for the solid phase pressure later on in the same function. The pressure is computed due to p.45 eq. 3.22. However, the frictional pressure (p. 52 eq. 3.28) is not taken into account. The frictional terms seem to be ignored intentionally which makes doubt whether I missed something or this part of the implementation is based on more recent publication. Nevertheless, I attached a kineticTheoryModel.C file which includes an in my opinion correct implementation. I ran the fluidized bed test case (RAS and kinetic Theory activated, equilibrium deactivated) with it and encountered no problems. (van Wachem thesis : https://repository.tudelft.nl/islandora/object/uuid:919e2efa-5db2-40e6-9082-83b1416709a6/?collection=research ) | ||||
Tags | multiphase | ||||
|
kineticTheoryModel.C (15,024 bytes)
/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 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 "kineticTheoryModel.H" #include "mathematicalConstants.H" #include "twoPhaseSystem.H" #include "fvOptions.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::RASModels::kineticTheoryModel::kineticTheoryModel ( const volScalarField& alpha, const volScalarField& rho, const volVectorField& U, const surfaceScalarField& alphaRhoPhi, const surfaceScalarField& phi, const transportModel& phase, const word& propertiesName, const word& type ) : eddyViscosity < RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>> > ( type, alpha, rho, U, alphaRhoPhi, phi, phase, propertiesName ), phase_(phase), viscosityModel_ ( kineticTheoryModels::viscosityModel::New ( coeffDict_ ) ), conductivityModel_ ( kineticTheoryModels::conductivityModel::New ( coeffDict_ ) ), radialModel_ ( kineticTheoryModels::radialModel::New ( coeffDict_ ) ), granularPressureModel_ ( kineticTheoryModels::granularPressureModel::New ( coeffDict_ ) ), frictionalStressModel_ ( kineticTheoryModels::frictionalStressModel::New ( coeffDict_ ) ), equilibrium_(coeffDict_.lookup("equilibrium")), e_("e", dimless, coeffDict_), alphaMax_("alphaMax", dimless, coeffDict_), alphaMinFriction_ ( "alphaMinFriction", dimless, coeffDict_ ), residualAlpha_ ( "residualAlpha", dimless, coeffDict_ ), maxNut_ ( "maxNut", dimensionSet(0,2,-1,0,0), coeffDict_.lookupOrDefault<scalar>("maxNut",1000) ), Theta_ ( IOobject ( IOobject::groupName("Theta", phase.name()), U.time().timeName(), U.mesh(), IOobject::MUST_READ, IOobject::AUTO_WRITE ), U.mesh() ), lambda_ ( IOobject ( IOobject::groupName("lambda", phase.name()), U.time().timeName(), U.mesh(), IOobject::NO_READ, IOobject::NO_WRITE ), U.mesh(), dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0) ), gs0_ ( IOobject ( IOobject::groupName("gs0", phase.name()), U.time().timeName(), U.mesh(), IOobject::NO_READ, IOobject::NO_WRITE ), U.mesh(), dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 0.0) ), kappa_ ( IOobject ( IOobject::groupName("kappa", phase.name()), U.time().timeName(), U.mesh(), IOobject::NO_READ, IOobject::NO_WRITE ), U.mesh(), dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0) ), nuFric_ ( IOobject ( IOobject::groupName("nuFric", phase.name()), U.time().timeName(), U.mesh(), IOobject::NO_READ, IOobject::AUTO_WRITE ), U.mesh(), dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0) ) { if (type == typeName) { printCoeffs(type); } } // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // Foam::RASModels::kineticTheoryModel::~kineticTheoryModel() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // bool Foam::RASModels::kineticTheoryModel::read() { if ( eddyViscosity < RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>> >::read() ) { coeffDict().lookup("equilibrium") >> equilibrium_; e_.readIfPresent(coeffDict()); alphaMax_.readIfPresent(coeffDict()); alphaMinFriction_.readIfPresent(coeffDict()); viscosityModel_->read(); conductivityModel_->read(); radialModel_->read(); granularPressureModel_->read(); frictionalStressModel_->read(); return true; } else { return false; } } Foam::tmp<Foam::volScalarField> Foam::RASModels::kineticTheoryModel::k() const { NotImplemented; return nut_; } Foam::tmp<Foam::volScalarField> Foam::RASModels::kineticTheoryModel::epsilon() const { NotImplemented; return nut_; } Foam::tmp<Foam::volSymmTensorField> Foam::RASModels::kineticTheoryModel::R() const { return tmp<volSymmTensorField> ( new volSymmTensorField ( IOobject ( IOobject::groupName("R", U_.group()), runTime_.timeName(), mesh_, IOobject::NO_READ, IOobject::NO_WRITE ), - (nut_)*dev(twoSymm(fvc::grad(U_))) - (lambda_*fvc::div(phi_))*symmTensor::I ) ); } Foam::tmp<Foam::volScalarField> Foam::RASModels::kineticTheoryModel::pPrime() const { const volScalarField& rho = phase_.rho(); tmp<volScalarField> tpPrime ( Theta_ *granularPressureModel_->granularPressureCoeffPrime ( alpha_, radialModel_->g0(alpha_, alphaMinFriction_, alphaMax_), radialModel_->g0prime(alpha_, alphaMinFriction_, alphaMax_), rho, e_ ) + frictionalStressModel_->frictionalPressurePrime ( phase_, alphaMinFriction_, alphaMax_ ) ); volScalarField::Boundary& bpPrime = tpPrime.ref().boundaryFieldRef(); forAll(bpPrime, patchi) { if (!bpPrime[patchi].coupled()) { bpPrime[patchi] == 0; } } return tpPrime; } Foam::tmp<Foam::surfaceScalarField> Foam::RASModels::kineticTheoryModel::pPrimef() const { return fvc::interpolate(pPrime()); } Foam::tmp<Foam::volSymmTensorField> Foam::RASModels::kineticTheoryModel::devRhoReff() const { return tmp<volSymmTensorField> ( new volSymmTensorField ( IOobject ( IOobject::groupName("devRhoReff", U_.group()), runTime_.timeName(), mesh_, IOobject::NO_READ, IOobject::NO_WRITE ), - (rho_*nut_) *dev(twoSymm(fvc::grad(U_))) - ((rho_*lambda_)*fvc::div(phi_))*symmTensor::I ) ); } Foam::tmp<Foam::fvVectorMatrix> Foam::RASModels::kineticTheoryModel::divDevRhoReff ( volVectorField& U ) const { return ( - fvm::laplacian(rho_*nut_, U) - fvc::div ( (rho_*nut_)*dev2(T(fvc::grad(U))) + ((rho_*lambda_)*fvc::div(phi_)) *dimensioned<symmTensor>("I", dimless, symmTensor::I) ) ); } void Foam::RASModels::kineticTheoryModel::correct() { // Local references volScalarField alpha(max(alpha_, scalar(0))); const volScalarField& rho = phase_.rho(); const surfaceScalarField& alphaRhoPhi = alphaRhoPhi_; const volVectorField& U = U_; const volVectorField& Uc_ = refCast<const twoPhaseSystem>(phase_.fluid()).otherPhase(phase_).U(); const scalar sqrtPi = sqrt(constant::mathematical::pi); dimensionedScalar ThetaSmall("ThetaSmall", Theta_.dimensions(), 1.0e-6); dimensionedScalar ThetaSmallSqrt(sqrt(ThetaSmall)); tmp<volScalarField> tda(phase_.d()); const volScalarField& da = tda(); tmp<volTensorField> tgradU(fvc::grad(U_)); const volTensorField& gradU(tgradU()); volSymmTensorField D(symm(gradU)); // Calculating the radial distribution function gs0_ = radialModel_->g0(alpha, alphaMinFriction_, alphaMax_); if (!equilibrium_) { // Particle viscosity (Table 3.2, p.47) nut_ = viscosityModel_->nu(alpha, Theta_, gs0_, rho, da, e_); // Frictional pressure volScalarField pf(frictionalStressModel_->frictionalPressure(phase_,alphaMinFriction_,alphaMax_)); // Frictional viscosity nuFric_ = frictionalStressModel_->nu(phase_,alphaMinFriction_,alphaMax_,pf/rho,D); // Limit viscosity and add frictional viscosity nut_.min(maxNut_); nuFric_ = min(nuFric_, maxNut_ - nut_); nut_ += nuFric_; volScalarField ThetaSqrt("sqrtTheta", sqrt(Theta_)); // Bulk viscosity p. 45 (Lun et al. 1984). lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi; // Stress tensor, Definitions, Table 3.1, p. 43 volSymmTensorField tau ( rho*(2.0*nut_*D + (lambda_ - (2.0/3.0)*nut_)*tr(D)*I) ); // Dissipation (Eq. 3.24, p.50) volScalarField gammaCoeff ( "gammaCoeff", 12.0*(1.0 - sqr(e_)) *max(sqr(alpha), residualAlpha_) *rho*gs0_*(1.0/da)*ThetaSqrt/sqrtPi ); // Drag volScalarField beta ( refCast<const twoPhaseSystem>(phase_.fluid()).drag(phase_).K() ); // Eq. 3.25, p. 50 Js = J1 - J2 volScalarField J1("J1", 3.0*beta); volScalarField J2 ( "J2", 0.25*sqr(beta)*da*magSqr(U - Uc_) /( max(alpha, residualAlpha_)*rho *sqrtPi*(ThetaSqrt + ThetaSmallSqrt) ) ); // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45) volScalarField PsCoeff ( granularPressureModel_->granularPressureCoeff ( alpha, gs0_, rho, e_ ) ); volScalarField Pf ( frictionalStressModel_->frictionalPressure ( phase_, alphaMinFriction_, alphaMax_ ) ); // 'thermal' conductivity (Table 3.3, p. 49) kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_); fv::options& fvOptions(fv::options::New(mesh_)); // Construct the granular temperature equation (Eq. 3.20, p. 44) // NB. note that there are two typos in Eq. 3.20: // Ps should be without grad // the laplacian has the wrong sign fvScalarMatrix ThetaEqn ( 1.5* ( fvm::ddt(alpha, rho, Theta_) + fvm::div(alphaRhoPhi, Theta_) - fvc::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaRhoPhi), Theta_) ) - fvm::laplacian(kappa_, Theta_, "laplacian(kappa,Theta)") == - fvm::SuSp((PsCoeff*I) && gradU, Theta_) + ((tau - Pf*I) && gradU) + fvm::Sp(-gammaCoeff, Theta_) + fvm::Sp(-J1, Theta_) + fvm::Sp(J2/(Theta_ + ThetaSmall), Theta_) + fvOptions(alpha, rho, Theta_) ); ThetaEqn.relax(); fvOptions.constrain(ThetaEqn); ThetaEqn.solve(); fvOptions.correct(Theta_); } else { // Equilibrium => dissipation == production // Eq. 4.14, p.82 volScalarField K1("K1", 2.0*(1.0 + e_)*rho*gs0_); volScalarField K3 ( "K3", 0.5*da*rho* ( (sqrtPi/(3.0*(3.0 - e_))) *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha*gs0_) +1.6*alpha*gs0_*(1.0 + e_)/sqrtPi ) ); volScalarField K2 ( "K2", 4.0*da*rho*(1.0 + e_)*alpha*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0 ); volScalarField K4("K4", 12.0*(1.0 - sqr(e_))*rho*gs0_/(da*sqrtPi)); volScalarField trD ( "trD", alpha/(alpha + residualAlpha_) *fvc::div(phi_) ); volScalarField tr2D("tr2D", sqr(trD)); volScalarField trD2("trD2", tr(D & D)); volScalarField t1("t1", K1*alpha + rho); volScalarField l1("l1", -t1*trD); volScalarField l2("l2", sqr(t1)*tr2D); volScalarField l3 ( "l3", 4.0 *K4 *alpha *(2.0*K3*trD2 + K2*tr2D) ); Theta_ = sqr ( (l1 + sqrt(l2 + l3)) /(2.0*max(alpha, residualAlpha_)*K4) ); kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_); } Theta_.max(0); Theta_.min(100); { // particle viscosity (Table 3.2, p.47) nut_ = viscosityModel_->nu(alpha, Theta_, gs0_, rho, da, e_); volScalarField ThetaSqrt("sqrtTheta", sqrt(Theta_)); // Bulk viscosity p. 45 (Lun et al. 1984). lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi; // Frictional pressure volScalarField pf ( frictionalStressModel_->frictionalPressure ( phase_, alphaMinFriction_, alphaMax_ ) ); nuFric_ = frictionalStressModel_->nu ( phase_, alphaMinFriction_, alphaMax_, pf/rho, D ); // Limit viscosity and add frictional viscosity nut_.min(maxNut_); nuFric_ = min(nuFric_, maxNut_ - nut_); nut_ += nuFric_; } if (debug) { Info<< typeName << ':' << nl << " max(Theta) = " << max(Theta_).value() << nl << " max(nut) = " << max(nut_).value() << endl; } } // ************************************************************************* // |
|
Hi This issue was discussed a couple of years ago, see https://bugs.openfoam.org/view.php?id=1831. |
|
Hi tniemi, thanks for your help. In this case the issue is resolved. Although I understand literature different. However, this has been discussed already. |
|
See https://bugs.openfoam.org/view.php?id=1831 |
Date Modified | Username | Field | Change |
---|---|---|---|
2017-08-03 10:43 | Arne | New Issue | |
2017-08-03 10:43 | Arne | File Added: kineticTheoryModel.C | |
2017-08-03 10:43 | Arne | Tag Attached: multiphase | |
2017-08-03 10:43 | Arne | Tag Attached: kineticTheory | |
2017-08-03 11:34 | tniemi | Note Added: 0008502 | |
2017-08-03 12:48 | Arne | Note Added: 0008503 | |
2017-08-04 11:14 | henry | Assigned To | => henry |
2017-08-04 11:14 | henry | Status | new => closed |
2017-08-04 11:14 | henry | Resolution | open => no change required |
2017-08-04 11:14 | henry | Note Added: 0008509 | |
2017-08-04 11:14 | henry | Relationship added | related to 0001831 |
2018-07-10 11:25 | administrator | Tag Detached: kineticTheory |