00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "fvDOM.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028
00029 #include <OpenFOAM/mathematicalConstants.H>
00030 #include <radiation/radiationConstants.H>
00031
00032
00033
00034 namespace Foam
00035 {
00036 namespace radiation
00037 {
00038 defineTypeNameAndDebug(fvDOM, 0);
00039
00040 addToRunTimeSelectionTable
00041 (
00042 radiationModel,
00043 fvDOM,
00044 dictionary
00045 );
00046 }
00047 }
00048
00049
00050
00051
00052 Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
00053 :
00054 radiationModel(typeName, T),
00055 G_
00056 (
00057 IOobject
00058 (
00059 "G",
00060 mesh_.time().timeName(),
00061 mesh_,
00062 IOobject::NO_READ,
00063 IOobject::AUTO_WRITE
00064 ),
00065 mesh_,
00066 dimensionedScalar("G", dimMass/pow3(dimTime), 0.0)
00067 ),
00068 Qr_
00069 (
00070 IOobject
00071 (
00072 "Qr",
00073 mesh_.time().timeName(),
00074 mesh_,
00075 IOobject::NO_READ,
00076 IOobject::AUTO_WRITE
00077 ),
00078 mesh_,
00079 dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
00080 ),
00081 Qem_
00082 (
00083 IOobject
00084 (
00085 "Qem",
00086 mesh_.time().timeName(),
00087 mesh_,
00088 IOobject::NO_READ,
00089 IOobject::NO_WRITE
00090 ),
00091 mesh_,
00092 dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
00093 ),
00094 Qin_
00095 (
00096 IOobject
00097 (
00098 "Qin",
00099 mesh_.time().timeName(),
00100 mesh_,
00101 IOobject::NO_READ,
00102 IOobject::NO_WRITE
00103 ),
00104 mesh_,
00105 dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
00106 ),
00107 a_
00108 (
00109 IOobject
00110 (
00111 "a",
00112 mesh_.time().timeName(),
00113 mesh_,
00114 IOobject::NO_READ,
00115 IOobject::NO_WRITE
00116 ),
00117 mesh_,
00118 dimensionedScalar("a", dimless/dimLength, 0.0)
00119 ),
00120 e_
00121 (
00122 IOobject
00123 (
00124 "e",
00125 mesh_.time().timeName(),
00126 mesh_,
00127 IOobject::NO_READ,
00128 IOobject::NO_WRITE
00129 ),
00130 mesh_,
00131 dimensionedScalar("a", dimless/dimLength, 0.0)
00132 ),
00133 E_
00134 (
00135 IOobject
00136 (
00137 "E",
00138 mesh_.time().timeName(),
00139 mesh_,
00140 IOobject::NO_READ,
00141 IOobject::NO_WRITE
00142 ),
00143 mesh_,
00144 dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
00145 ),
00146 nTheta_(readLabel(coeffs_.lookup("nTheta"))),
00147 nPhi_(readLabel(coeffs_.lookup("nPhi"))),
00148 nRay_(0),
00149 nLambda_(absorptionEmission_->nBands()),
00150 aLambda_(nLambda_),
00151 blackBody_(nLambda_, T),
00152 IRay_(0),
00153 convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
00154 maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50))
00155 {
00156 if (mesh_.nSolutionD() == 3)
00157 {
00158 nRay_ = 4*nPhi_*nTheta_;
00159 IRay_.setSize(nRay_);
00160 scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
00161 scalar deltaTheta = mathematicalConstant::pi/nTheta_;
00162 label i = 0;
00163 for (label n = 1; n <= nTheta_; n++)
00164 {
00165 for (label m = 1; m <= 4*nPhi_; m++)
00166 {
00167 scalar thetai = (2.0*n - 1.0)*deltaTheta/2.0;
00168 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
00169 IRay_.set
00170 (
00171 i,
00172 new radiativeIntensityRay
00173 (
00174 *this,
00175 mesh_,
00176 phii,
00177 thetai,
00178 deltaPhi,
00179 deltaTheta,
00180 nLambda_,
00181 absorptionEmission_,
00182 blackBody_,
00183 i
00184 )
00185 );
00186 i++;
00187 }
00188 }
00189 }
00190 else
00191 {
00192 if (mesh_.nSolutionD() == 2)
00193 {
00194 scalar thetai = mathematicalConstant::piByTwo;
00195 scalar deltaTheta = mathematicalConstant::pi;
00196 nRay_ = 4*nPhi_;
00197 IRay_.setSize(nRay_);
00198 scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
00199 label i = 0;
00200 for (label m = 1; m <= 4*nPhi_; m++)
00201 {
00202 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
00203 IRay_.set
00204 (
00205 i,
00206 new radiativeIntensityRay
00207 (
00208 *this,
00209 mesh_,
00210 phii,
00211 thetai,
00212 deltaPhi,
00213 deltaTheta,
00214 nLambda_,
00215 absorptionEmission_,
00216 blackBody_,
00217 i
00218 )
00219 );
00220 i++;
00221 }
00222 }
00223 else
00224 {
00225 scalar thetai = mathematicalConstant::piByTwo;
00226 scalar deltaTheta = mathematicalConstant::pi;
00227 nRay_ = 2;
00228 IRay_.setSize(nRay_);
00229 scalar deltaPhi = mathematicalConstant::pi;
00230 label i = 0;
00231 for (label m = 1; m <= 2; m++)
00232 {
00233 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
00234 IRay_.set
00235 (
00236 i,
00237 new radiativeIntensityRay
00238 (
00239 *this,
00240 mesh_,
00241 phii,
00242 thetai,
00243 deltaPhi,
00244 deltaTheta,
00245 nLambda_,
00246 absorptionEmission_,
00247 blackBody_,
00248 i
00249 )
00250 );
00251 i++;
00252 }
00253
00254 }
00255 }
00256
00257
00258
00259 forAll(aLambda_, lambdaI)
00260 {
00261 aLambda_.set
00262 (
00263 lambdaI,
00264 new volScalarField
00265 (
00266 IOobject
00267 (
00268 "aLambda_" + Foam::name(lambdaI) ,
00269 mesh_.time().timeName(),
00270 mesh_,
00271 IOobject::NO_READ,
00272 IOobject::NO_WRITE
00273 ),
00274 a_
00275 )
00276 );
00277 }
00278
00279 Info<< "fvDOM : Allocated " << IRay_.size()
00280 << " rays with average orientation:" << nl;
00281 forAll (IRay_, i)
00282 {
00283 Info<< '\t' << IRay_[i].I().name()
00284 << '\t' << IRay_[i].dAve() << nl;
00285 }
00286 Info<< endl;
00287 }
00288
00289
00290
00291
00292 Foam::radiation::fvDOM::~fvDOM()
00293 {}
00294
00295
00296
00297
00298 bool Foam::radiation::fvDOM::read()
00299 {
00300 if (radiationModel::read())
00301 {
00302
00303
00304 coeffs_.readIfPresent("convergence", convergence_);
00305 coeffs_.readIfPresent("maxIter", maxIter_);
00306
00307 return true;
00308 }
00309 else
00310 {
00311 return false;
00312 }
00313 }
00314
00315
00316 void Foam::radiation::fvDOM::calculate()
00317 {
00318 absorptionEmission_->correct(a_, aLambda_);
00319
00320 updateBlackBodyEmission();
00321
00322 scalar maxResidual = 0.0;
00323 label radIter = 0;
00324 do
00325 {
00326 radIter++;
00327 forAll(IRay_, rayI)
00328 {
00329 maxResidual = 0.0;
00330 scalar maxBandResidual = IRay_[rayI].correct();
00331 maxResidual = max(maxBandResidual, maxResidual);
00332 }
00333
00334 Info << "Radiation solver iter: " << radIter << endl;
00335
00336 } while(maxResidual > convergence_ && radIter < maxIter_);
00337
00338 updateG();
00339 }
00340
00341
00342 Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
00343 {
00344 return tmp<volScalarField>
00345 (
00346 new volScalarField
00347 (
00348 IOobject
00349 (
00350 "Rp",
00351 mesh_.time().timeName(),
00352 mesh_,
00353 IOobject::NO_READ,
00354 IOobject::NO_WRITE,
00355 false
00356 ),
00357 4.0*a_*radiation::sigmaSB
00358 )
00359 );
00360 }
00361
00362
00363 Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
00364 Foam::radiation::fvDOM::Ru() const
00365 {
00366
00367 const DimensionedField<scalar, volMesh>& G =
00368 G_.dimensionedInternalField();
00369 const DimensionedField<scalar, volMesh> E =
00370 absorptionEmission_->ECont()().dimensionedInternalField();
00371 const DimensionedField<scalar, volMesh> a =
00372 a_.dimensionedInternalField();
00373
00374 return a*G - E;
00375 }
00376
00377
00378 void Foam::radiation::fvDOM::updateBlackBodyEmission()
00379 {
00380 for (label j=0; j < nLambda_; j++)
00381 {
00382 blackBody_.correct(j, absorptionEmission_->bands(j));
00383 }
00384 }
00385
00386
00387 void Foam::radiation::fvDOM::updateG()
00388 {
00389 G_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
00390 Qr_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
00391 Qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
00392 Qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
00393
00394 forAll(IRay_, rayI)
00395 {
00396 IRay_[rayI].addIntensity();
00397 G_ += IRay_[rayI].I()*IRay_[rayI].omega();
00398 Qr_.boundaryField() += IRay_[rayI].Qr().boundaryField();
00399 Qem_.boundaryField() += IRay_[rayI].Qem().boundaryField();
00400 Qin_.boundaryField() += IRay_[rayI].Qin().boundaryField();
00401 }
00402 }
00403
00404
00405 void Foam::radiation::fvDOM::setRayIdLambdaId
00406 (
00407 const word& name,
00408 label& rayId,
00409 label& lambdaId
00410 ) const
00411 {
00412
00413 size_type i1 = name.find_first_of("_");
00414 size_type i2 = name.find_last_of("_");
00415
00416 rayId = readLabel(IStringStream(name.substr(i1+1, i2-1))());
00417 lambdaId = readLabel(IStringStream(name.substr(i2+1, name.size()-1))());
00418 }
00419
00420
00421