From f2ba82e3b567488be556b582577c8b70392b70f5 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Wed, 6 Dec 2017 11:06:00 +0000 Subject: [PATCH] fvDOM: Extended to any axis aligned 1-D and 2-D geometry --- ...iffusiveRadiationMixedFvPatchScalarField.C | 2 +- .../radiationModels/fvDOM/fvDOM/fvDOM.C | 38 +++++++----------- .../radiationModels/fvDOM/fvDOM/fvDOM.H | 11 +++-- .../radiativeIntensityRay.C | 40 ++++++++++++++++++- 4 files changed, 61 insertions(+), 30 deletions(-) diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C index ec4f0ba7d..7d0b2aeb6 100644 --- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C +++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C @@ -192,7 +192,7 @@ updateCoeffs() scalarField& qem = ray.qem().boundaryFieldRef()[patchi]; scalarField& qin = ray.qin().boundaryFieldRef()[patchi]; - const vector& myRayId = dom.IRay(rayId).d(); + const vector& myRayId = dom.IRay(rayId).dAve(); // Use updated Ir while iterating over rays // avoids to used lagged qin diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C index 94d365081..f0cb9d0e3 100644 --- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C +++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C @@ -87,14 +87,6 @@ void Foam::radiation::fvDOM::initialise() // 2D else if (mesh_.nSolutionD() == 2) { - // Currently 2D solution is limited to the x-y plane - if (mesh_.solutionD()[vector::Z] != -1) - { - FatalErrorInFunction - << "Currently 2D solution is limited to the x-y plane" - << exit(FatalError); - } - scalar thetai = piByTwo; scalar deltaTheta = pi; nRay_ = 4*nPhi_; @@ -127,14 +119,6 @@ void Foam::radiation::fvDOM::initialise() // 1D else { - // Currently 1D solution is limited to the x-direction - if (mesh_.solutionD()[vector::X] != 1) - { - FatalErrorInFunction - << "Currently 1D solution is limited to the x-direction" - << exit(FatalError); - } - scalar thetai = piByTwo; scalar deltaTheta = pi; nRay_ = 2; @@ -187,20 +171,28 @@ void Foam::radiation::fvDOM::initialise() ); } - Info<< "fvDOM : Allocated " << IRay_.size() - << " rays with average orientation:" << nl; + // Calculate the maximum solid angle forAll(IRay_, rayId) { if (omegaMax_ < IRay_[rayId].omega()) { omegaMax_ = IRay_[rayId].omega(); } - Info<< '\t' << IRay_[rayId].I().name() << " : " << "omega : " - << '\t' << IRay_[rayId].omega() << nl; } - Info<< endl; + + Info<< typeName << ": Created " << IRay_.size() << " rays with average " + << "directions (dAve) and solid angles (omega)" << endl; + Info<< incrIndent; + forAll(IRay_, rayId) + { + Info<< indent + << "Ray " << IRay_[rayId].I().name() << ": " + << "dAve = " << IRay_[rayId].dAve() << ", " + << "omega = " << IRay_[rayId].omega() << endl; + } + Info<< decrIndent << endl; } @@ -242,8 +234,8 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T) "qem", mesh_.time().timeName(), mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE ), mesh_, dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0) diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H index cbcbf4d64..ba4fffe8b 100644 --- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H +++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H @@ -49,11 +49,14 @@ Description solverFreq 1; // Number of flow iterations per radiation iteration \endverbatim - The total number of solid angles is 4*nPhi*nTheta. + In 1-D the ray directions are bound to one of the X, Y or Z directions. The + total number of solid angles is 2. nPhi and nTheta are ignored. - In 1D the direction of the rays is X (nPhi and nTheta are ignored) - In 2D the direction of the rays is on X-Y plane (only nPhi is considered) - In 3D (nPhi and nTheta are considered) + In 2-D the ray directions are within one of the X-Y, X-Z or Y-Z planes. The + total number of solid angles is 4*nPhi. nTheta is ignored. + + In 3D the rays span all directions. The total number of solid angles is + 4*nPhi*nTheta. SourceFiles fvDOM.C diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C index 60ac75c60..d87cb7b1e 100644 --- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C +++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C @@ -137,6 +137,42 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay 0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta) ); + // Transform directions so that they fall inside the bounds of reduced + // dimension cases + if (mesh_.nSolutionD() == 2) + { + vector meshDir(vector::zero); + for (direction cmpt=0; cmpt IDefaultPtr; forAll(ILambda_, lambdaI) @@ -210,12 +246,12 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct() scalar maxResidual = -GREAT; + const surfaceScalarField Ji(dAve_ & mesh_.Sf()); + forAll(ILambda_, lambdaI) { const volScalarField& k = dom_.aLambda(lambdaI); - const surfaceScalarField Ji(dAve_ & mesh_.Sf()); - fvScalarMatrix IiEq ( fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")