src/dbns/numericFlux/numericFlux.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "numericFlux.H"
  27
  28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  29
  30 // Construct from components
  31 template<class Flux, class Limiter>
  32 Foam::numericFlux<Flux, Limiter>::numericFlux
  33 (
  34     const volScalarField& p,
  35     const volVectorField& U,
  36     const volScalarField& T,
  37     basicThermo& thermo
  38 )
  39 :
  40     numericFluxBase<Flux>(),
  41     mesh_(p.mesh()),
  42     p_(p),
  43     U_(U),
  44     T_(T),
  45     thermo_(thermo),
  46     rhoFlux_
  47     (
  48         IOobject
  49         (
  50             "phi",
  51             mesh_.time().timeName(),
  52             mesh_,
  53             IOobject::NO_READ,
  54             IOobject::NO_WRITE
  55         ),
  56         (linearInterpolate(thermo_.rho()*U_) & mesh_.Sf())
  57     ),
  58     rhoUFlux_
  59     (
  60         IOobject
  61         (
  62             "rhoUFlux",
  63             mesh_.time().timeName(),
  64             mesh_,
  65             IOobject::NO_READ,
  66             IOobject::NO_WRITE
  67         ),
  68         rhoFlux_*linearInterpolate(U_)
  69     ),
  70     rhoEFlux_
  71     (
  72         IOobject
  73         (
  74             "rhoEFlux",
  75             mesh_.time().timeName(),
  76             mesh_,
  77             IOobject::NO_READ,
  78             IOobject::NO_WRITE
  79         ),
  80         rhoFlux_*linearInterpolate(thermo.Cv()*T_ + 0.5*magSqr(U_))
  81     ),
  82     gradP_(fvc::grad(p_)),
  83     gradU_(fvc::grad(U_)),
  84     gradT_(fvc::grad(T_))
  85 {}
  86
  87
  88 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  89
  90 template<class Flux, class Limiter>
  91 void Foam::numericFlux<Flux, Limiter>::computeFlux()
  92 {
  93     // Get face-to-cell addressing: face area point from owner to neighbour
  94     const unallocLabelList& owner = mesh_.owner();
  95     const unallocLabelList& neighbour = mesh_.neighbour();
  96
  97     // Get the face area vector
  98     const surfaceVectorField& Sf = mesh_.Sf();
  99     const surfaceScalarField& magSf = mesh_.magSf();
 100
 101     const volVectorField& cellCentre = mesh_.C();
 102     const surfaceVectorField& faceCentre = mesh_.Cf();
 103
 104     // Thermodynamics
 105     const volScalarField Cv = thermo_.Cv();
 106     const volScalarField R  = thermo_.Cp() - Cv;
 107
 108     gradP_ = fvc::grad(p_);
 109     gradP_.correctBoundaryConditions();
 110
 111     gradU_ = fvc::grad(U_);
 112     gradU_.correctBoundaryConditions();
 113
 114     gradT_ = fvc::grad(T_);
 115     gradT_.correctBoundaryConditions();
 116
 117     MDLimiter<scalar, Limiter> scalarPLimiter
 118     (
 119         this->p_,
 120         this->gradP_
 121     );
 122
 123     MDLimiter<vector, Limiter> vectorULimiter
 124     (
 125         this->U_,
 126         this->gradU_
 127     );
 128
 129     MDLimiter<scalar, Limiter> scalarTLimiter
 130     (
 131         this->T_,
 132         this->gradT_
 133     );
 134
 135     // Get limiters
 136     const volScalarField& pLimiter = scalarPLimiter.phiLimiter();
 137     const volVectorField& ULimiter = vectorULimiter.phiLimiter();
 138     const volScalarField& TLimiter = scalarTLimiter.phiLimiter();
 139
 140     // Calculate fluxes at internal faces
 141     forAll (owner, faceI)
 142     {
 143         const label own = owner[faceI];
 144         const label nei = neighbour[faceI];
 145
 146         const vector deltaRLeft = faceCentre[faceI] - cellCentre[own];
 147         const vector deltaRRight = faceCentre[faceI] - cellCentre[nei];
 148
 149         // calculate fluxes with reconstructed primitive variables at faces
 150         Flux::evaluateFlux
 151         (
 152             rhoFlux_[faceI],
 153             rhoUFlux_[faceI],
 154             rhoEFlux_[faceI],
 155             p_[own] + pLimiter[own]*(deltaRLeft & gradP_[own]),
 156             p_[nei] + pLimiter[nei]*(deltaRRight & gradP_[nei]),
 157             U_[own] + cmptMultiply(ULimiter[own], (deltaRLeft & gradU_[own])),
 158             U_[nei] + cmptMultiply(ULimiter[nei], (deltaRRight & gradU_[nei])),
 159             T_[own] + TLimiter[own]*(deltaRLeft & gradT_[own]),
 160             T_[nei] + TLimiter[nei]*(deltaRRight & gradT_[nei]),
 161             R[own],
 162             R[nei],
 163             Cv[own],
 164             Cv[nei],
 165             Sf[faceI],
 166             magSf[faceI]
 167         );
 168     }
 169
 170     // Update boundary field and values
 171     forAll (rhoFlux_.boundaryField(), patchi)
 172     {
 173         const fvPatch& curPatch = p_.boundaryField()[patchi].patch();
 174
 175         // Fluxes
 176         fvsPatchScalarField& pRhoFlux  = rhoFlux_.boundaryField()[patchi];
 177         fvsPatchVectorField& pRhoUFlux = rhoUFlux_.boundaryField()[patchi];
 178         fvsPatchScalarField& pRhoEFlux = rhoEFlux_.boundaryField()[patchi];
 179
 180         // Patch fields
 181         const fvPatchScalarField& pp = p_.boundaryField()[patchi];
 182         const vectorField& pU = U_.boundaryField()[patchi];
 183         const scalarField& pT = T_.boundaryField()[patchi];
 184
 185         const scalarField& pCv = Cv.boundaryField()[patchi];
 186         const scalarField& pR = R.boundaryField()[patchi];
 187
 188         // Gradients
 189         const fvPatchVectorField& pGradP = gradP_.boundaryField()[patchi];
 190         const fvPatchTensorField& pGradU = gradU_.boundaryField()[patchi];
 191         const fvPatchVectorField& pGradT = gradT_.boundaryField()[patchi];
 192
 193         // Limiters
 194         const fvPatchScalarField& pPatchLim = pLimiter.boundaryField()[patchi];
 195         const fvPatchVectorField& UPatchLim = ULimiter.boundaryField()[patchi];
 196         const fvPatchScalarField& TPatchLim = TLimiter.boundaryField()[patchi];
 197
 198         // Face areas
 199         const fvsPatchVectorField& pSf = Sf.boundaryField()[patchi];
 200         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
 201
 202         if (pp.coupled())
 203         {
 204             // Coupled patch
 205             const scalarField ppLeft  =
 206                 p_.boundaryField()[patchi].patchInternalField();
 207
 208             const scalarField ppRight =
 209                 p_.boundaryField()[patchi].patchNeighbourField();
 210
 211             const vectorField pULeft  =
 212                 U_.boundaryField()[patchi].patchInternalField();
 213
 214             const vectorField pURight =
 215                 U_.boundaryField()[patchi].patchNeighbourField();
 216
 217             const scalarField pTLeft  =
 218                 T_.boundaryField()[patchi].patchInternalField();
 219
 220             const scalarField pTRight =
 221                 T_.boundaryField()[patchi].patchNeighbourField();
 222
 223             // Gradients
 224             const vectorField pgradPLeft = pGradP.patchInternalField();
 225             const vectorField pgradPRight = pGradP.patchNeighbourField();
 226
 227             const tensorField pgradULeft = pGradU.patchInternalField();
 228             const tensorField pgradURight = pGradU.patchNeighbourField();
 229
 230             const vectorField pgradTLeft = pGradT.patchInternalField();
 231             const vectorField pgradTRight = pGradT.patchNeighbourField();
 232
 233             // Geometry: call the raw cell-to-face vector by calling
 234             // the base patch (cell-to-face) delta coefficient
 235             // Work out the right delta from the cell-to-cell delta
 236             // across the coupled patch and left delta
 237             vectorField pDeltaRLeft = curPatch.fvPatch::delta();
 238             vectorField pDdeltaRRight = pDeltaRLeft - curPatch.delta();
 239
 240             // Limiters
 241
 242             const scalarField ppLimiterLeft = pPatchLim.patchInternalField();
 243             const scalarField ppLimiterRight = pPatchLim.patchNeighbourField();
 244
 245             const vectorField pULimiterLeft = UPatchLim.patchInternalField();
 246             const vectorField pULimiterRight = UPatchLim.patchNeighbourField();
 247
 248             const scalarField pTLimiterLeft = TPatchLim.patchInternalField();
 249             const scalarField pTLimiterRight = TPatchLim.patchNeighbourField();
 250
 251             forAll (pp, facei)
 252             {
 253                 Flux::evaluateFlux
 254                 (
 255                     pRhoFlux[facei],
 256                     pRhoUFlux[facei],
 257                     pRhoEFlux[facei],
 258
 259                     ppLeft[facei]
 260                   + ppLimiterLeft[facei]*
 261                     (pDeltaRLeft[facei] & pgradPLeft[facei]),
 262
 263                     ppRight[facei]
 264                   + ppLimiterRight[facei]*
 265                     (pDdeltaRRight[facei] & pgradPRight[facei]),
 266
 267                     pULeft[facei]
 268                   + cmptMultiply
 269                     (
 270                         pULimiterLeft[facei],
 271                         pDeltaRLeft[facei] & pgradULeft[facei]
 272                     ),
 273
 274                     pURight[facei]
 275                   + cmptMultiply
 276                     (
 277                         pULimiterRight[facei],
 278                         pDdeltaRRight[facei] & pgradURight[facei]
 279                     ),
 280
 281                     pTLeft[facei]
 282                   + pTLimiterLeft[facei]*
 283                     (pDeltaRLeft[facei] & pgradTLeft[facei]),
 284
 285                     pTRight[facei]
 286                   + pTLimiterRight[facei]*
 287                     (pDdeltaRRight[facei] & pgradTRight[facei]),
 288
 289                     pR[facei],
 290                     pR[facei],
 291                     pCv[facei],
 292                     pCv[facei],
 293                     pSf[facei],
 294                     pMagSf[facei]
 295                 );
 296             }
 297         }
 298         else
 299         {
 300             forAll (pp, facei)
 301             {
 302                 // Calculate fluxes
 303                 Flux::evaluateFlux
 304                 (
 305                     pRhoFlux[facei],
 306                     pRhoUFlux[facei],
 307                     pRhoEFlux[facei],
 308                     pp[facei],
 309                     pp[facei],
 310                     pU[facei],
 311                     pU[facei],
 312                     pT[facei],
 313                     pT[facei],
 314                     pR[facei],
 315                     pR[facei],
 316                     pCv[facei],
 317                     pCv[facei],
 318                     pSf[facei],
 319                     pMagSf[facei]
 320                 );
 321             }
 322         }
 323     }
 324 }
 325
 326
 327 // ************************************************************************* //