From e9f35a9d19bd8b4959bfc5c482adb1c37fbe667e Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Fri, 12 Jun 2015 15:35:43 +0100 Subject: [PATCH] Lagrangian radiation models and fvDOM: Apply updates Applying patches provided by Timo Niemi Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1636 --- .../Templates/ThermoParcel/ThermoParcel.C | 16 +++++++++---- .../radiationModels/fvDOM/fvDOM/fvDOM.C | 12 ++++++---- .../radiativeIntensityRay.C | 23 ++++++++++++++----- 3 files changed, 36 insertions(+), 15 deletions(-) diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C index 74ece70250..e5b469d4e3 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -210,6 +210,9 @@ void Foam::ThermoParcel::calc // Sum Ni*Cpi*Wi of emission species scalar NCpW = 0.0; + // Store T for consistent radiation source + const scalar T0 = this->T_; + // Calculate new particle temperature this->T_ = this->calcHeatTransfer @@ -255,7 +258,7 @@ void Foam::ThermoParcel::calc if (td.cloud().radiation()) { const scalar ap = this->areaP(); - const scalar T4 = pow4(this->T_); + const scalar T4 = pow4(T0); td.cloud().radAreaP()[cellI] += dt*np0*ap; td.cloud().radT4()[cellI] += dt*np0*T4; td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4; @@ -304,7 +307,7 @@ Foam::scalar Foam::ThermoParcel::calcHeatTransfer htc = max(htc, ROOTVSMALL); const scalar As = this->areaS(d); - scalar ap = Tc_ + Sh/As/htc; + scalar ap = Tc_ + Sh/(As*htc); scalar bp = 6.0*(Sh/As + htc*(Tc_ - T_)); if (td.cloud().radiation()) { @@ -313,8 +316,11 @@ Foam::scalar Foam::ThermoParcel::calcHeatTransfer const scalar sigma = physicoChemical::sigma.value(); const scalar epsilon = td.cloud().constProps().epsilon0(); - ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T_)/htc); - bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T_))); + // Assume constant source + scalar s = epsilon*(Gc/4.0 - sigma*pow4(T_)); + + ap += s/htc; + bp += 6.0*s; } bp /= rho*d*Cp_*(ap - T_) + ROOTVSMALL; diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C index dab8a9a184..9c6781cef8 100644 --- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C +++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C @@ -482,7 +482,8 @@ Foam::tmp Foam::radiation::fvDOM::Rp() const IOobject::NO_WRITE, false ), - 4.0*a_*physicoChemical::sigma //absorptionEmission_->a() + // Only include continuous phase emission + 4*absorptionEmission_->aCont()*physicoChemical::sigma ) ); } @@ -494,12 +495,15 @@ Foam::radiation::fvDOM::Ru() const const DimensionedField& G = G_.dimensionedInternalField(); + const DimensionedField E = absorptionEmission_->ECont()().dimensionedInternalField(); - const DimensionedField a = - a_.dimensionedInternalField(); - return a*G - E; + // Only include continuous phase absorption + const DimensionedField a = + absorptionEmission_->aCont()().dimensionedInternalField(); + + return a*G - E; } diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C index 64900e4bed..64e42b85bb 100644 --- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C +++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C @@ -30,7 +30,6 @@ License using namespace Foam::constant; - const Foam::word Foam::radiation::radiativeIntensityRay::intensityPrefix("ILambda"); @@ -225,9 +224,15 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct() + fvm::Sp(k*omega_, ILambda_[lambdaI]) == 1.0/constant::mathematical::pi*omega_ - * ( - k*blackBody_.bLambda(lambdaI) - + absorptionEmission_.ECont(lambdaI)/4 + *( + // Remove aDisp from k + (k - absorptionEmission_.aDisp(lambdaI)) + *blackBody_.bLambda(lambdaI) + + + absorptionEmission_.ECont(lambdaI) + + // Add EDisp term from parcels + + absorptionEmission_.EDisp(lambdaI) ) ); } @@ -240,8 +245,14 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct() == 1.0/constant::mathematical::pi*omega_ * ( - k*blackBody_.bLambda(lambdaI) - + absorptionEmission_.ECont(lambdaI)/4 + // Remove aDisp from k + (k - absorptionEmission_.aDisp(lambdaI)) + *blackBody_.bLambda(lambdaI) + + + absorptionEmission_.ECont(lambdaI) + + // Add EDisp term from parcels + + absorptionEmission_.EDisp(lambdaI) ) ); }