src/dbns/dbnsFlux/betaFlux/betaFlux.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 "betaFlux.H"
  27
  28 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  29
  30 void Foam::betaFlux::evaluateFlux
  31 (
  32     scalar& rhoFlux,
  33     vector& rhoUFlux,
  34     scalar& rhoEFlux,
  35     const scalar& pLeft,
  36     const scalar& pRight,
  37     const vector& ULeft,
  38     const vector& URight,
  39     const scalar& TLeft,
  40     const scalar& TRight,
  41     const scalar& RLeft,
  42     const scalar& RRight,
  43     const scalar& CvLeft,
  44     const scalar& CvRight,
  45     const vector& Sf,
  46     const scalar& magSf
  47 ) const
  48 {
  49     // Step 1: decode rho left and right:
  50     scalar rhoLeft = pLeft/(RLeft*TLeft);
  51     scalar rhoRight = pRight/(RRight*TRight);
  52
  53     // Decode left and right total energy:
  54     scalar eLeft = CvLeft*TLeft+0.5*magSqr(ULeft);
  55     scalar eRight = CvRight*TRight+0.5*magSqr(URight);
  56
  57     // Adiabatic exponent is constant for ideal gas but if Cp=Cp(T)
  58     // it must be computed for each cell and evaluated at each face
  59     // through reconstruction
  60     const scalar kappaLeft = (CvLeft+RLeft)/CvLeft;
  61     const scalar kappaRight = (CvRight+RRight)/CvRight;
  62
  63     // normal vector
  64     vector normalVector = Sf/magSf;
  65
  66     // Compute left and right contravariant velocities:
  67     const scalar contrVLeft  = (ULeft & normalVector);
  68     const scalar contrVRight = (URight & normalVector);
  69
  70     // Compute left and right total enthalpies:
  71     const scalar hLeft = eLeft + pLeft/rhoLeft;
  72     const scalar hRight = eRight + pRight/rhoRight;
  73
  74     // Compute left and right velocity square
  75     const scalar qLeftSquare  = magSqr(ULeft);
  76     const scalar qRightSquare = magSqr(URight);
  77
  78     // compute left and right speed of sound
  79     const scalar cLeft =
  80         sqrt(max((kappaLeft - 1)*(hLeft - 0.5*qLeftSquare), SMALL));
  81     const scalar cRight =
  82         sqrt(max((kappaRight - 1)*(hRight - 0.5*qRightSquare), SMALL));
  83
  84     const scalar magULeft  = mag(ULeft);
  85     const scalar magURight = mag(URight);
  86
  87     const scalar MLeft = mag(magULeft/cLeft);
  88     const scalar MRight = mag(magURight/cRight);
  89
  90     // Compute beta parameter - this should be done in multidimensional way
  91     // similarly to multidimensional limiters
  92     const scalar Mmax = max(MLeft, MRight);
  93     const scalar Mmin = min(MLeft, MRight);
  94
  95     scalar beta = 0;
  96
  97     if (Mmin > 0)
  98     {
  99         scalar alpha = 0.1;
 100         scalar r = sqrt(sqr(Mmax) - sqr(Mmin))/Mmin;
 101
 102         if (r >= alpha)
 103         {
 104             beta = min(r,1);
 105         }
 106     }
 107
 108     // Step 2: compute Roe averaged quantities for face:
 109     const scalar rhoTilde = sqrt(max(rhoLeft*rhoRight, SMALL));
 110
 111     // Some temporary variables:
 112     const scalar rhoLeftSqrt = sqrt(max(rhoLeft, SMALL));
 113     const scalar rhoRightSqrt = sqrt(max(rhoRight, SMALL));
 114
 115     const scalar wLeft = rhoLeftSqrt/(rhoLeftSqrt + rhoRightSqrt);
 116     const scalar wRight = 1 - wLeft;
 117
 118     const vector UTilde = ULeft*wLeft + URight*wRight;
 119     const scalar hTilde = hLeft*wLeft + hRight*wRight;
 120     const scalar qTildeSquare = magSqr(UTilde);
 121     const scalar kappaTilde = kappaLeft*wLeft + kappaRight*wRight;
 122
 123     // Speed of sound
 124     const scalar cTilde =
 125         sqrt(max((kappaTilde - 1)*(hTilde - 0.5*qTildeSquare), SMALL));
 126
 127     // Roe averaged contravariant velocity
 128     const scalar contrVTilde = (UTilde & normalVector);
 129
 130     // Step 3: compute primitive differences:
 131     const scalar deltaP = pRight - pLeft;
 132     const scalar deltaRho = rhoRight - rhoLeft;
 133     const vector deltaU = URight - ULeft;
 134     const scalar deltaContrV = (deltaU & normalVector);
 135
 136     // Step 4: compute wave strengths:
 137
 138     // Roe and Pike - formulation
 139     const scalar r1 =
 140         (deltaP - rhoTilde*cTilde*deltaContrV)/(2.0*sqr(cTilde));
 141     const scalar r2 = deltaRho - deltaP/sqr(cTilde);
 142     const scalar r3 =
 143         (deltaP + rhoTilde*cTilde*deltaContrV)/(2.0*sqr(cTilde));
 144
 145     // Step 5: compute l vectors
 146
 147     // rho row:
 148     const scalar l1rho = 1;
 149     const scalar l2rho = 1;
 150     const scalar l3rho = 0;
 151     const scalar l4rho = 1;
 152
 153     // first U column
 154     const vector l1U = UTilde - cTilde*normalVector;
 155
 156     // second U column
 157     const vector l2U = UTilde;
 158
 159     // third U column
 160     const vector l3U = deltaU - deltaContrV*normalVector;
 161
 162     // fourth U column
 163     const vector l4U = UTilde + cTilde*normalVector;
 164
 165     // E row
 166     const scalar l1e = hTilde - cTilde*contrVTilde;
 167     const scalar l2e = 0.5*qTildeSquare;
 168     const scalar l3e = (UTilde & deltaU) - contrVTilde*deltaContrV;
 169     const scalar l4e = hTilde + cTilde*contrVTilde;
 170
 171     // Step 6: compute eigenvalues
 172
 173     scalar lambda1 = mag(contrVTilde - cTilde);
 174     scalar lambda2 = mag(contrVTilde);
 175     scalar lambda3 = mag(contrVTilde + cTilde);
 176
 177     scalar lambdaMax = max(max(lambda1, lambda2), lambda3);
 178
 179     scalar lambda1Max = beta*lambda1 + (1 - beta)*lambdaMax;
 180     scalar lambda2Max = beta*lambda2 + (1 - beta)*lambdaMax;
 181     scalar lambda3Max = beta*lambda3 + (1 - beta)*lambdaMax;
 182
 183     // Step 7: check for Harten entropy correction
 184
 185
 186     // Components of deltaF1
 187     const scalar diffF11  = lambda1Max*r1*l1rho;
 188     const vector diffF124 = lambda1Max*r1*l1U;
 189     const scalar diffF15  = lambda1Max*r1*l1e;
 190
 191     // Components of deltaF2
 192     const scalar diffF21  = lambda2Max*(r2*l2rho + rhoTilde*l3rho);
 193     const vector diffF224 = lambda2Max*(r2*l2U + rhoTilde*l3U);
 194     const scalar diffF25  = lambda2Max*(r2*l2e + rhoTilde*l3e);
 195
 196     // Components of deltaF3
 197     const scalar diffF31  = lambda3Max*r3*l4rho;
 198     const vector diffF324 = lambda3Max*r3*l4U;
 199     const scalar diffF35  = lambda3Max*r3*l4e;
 200
 201     // Step 8: compute left and right fluxes
 202
 203     // Left flux 5-vector
 204     const scalar fluxLeft11 = rhoLeft*contrVLeft;
 205     const vector fluxLeft124 = ULeft*fluxLeft11 + normalVector*pLeft;
 206     const scalar fluxLeft15 = hLeft*fluxLeft11;
 207
 208     // Right flux 5-vector
 209     const scalar fluxRight11 = rhoRight*contrVRight;
 210     const vector fluxRight124 = URight*fluxRight11 + normalVector*pRight;
 211     const scalar fluxRight15 = hRight*fluxRight11;
 212
 213     // Step 9: compute face flux 5-vector
 214     const scalar flux1 =
 215         0.5*(fluxLeft11 + fluxRight11 - (diffF11 + diffF21 + diffF31));
 216     const vector flux24 =
 217         0.5*(fluxLeft124 + fluxRight124 - (diffF124 + diffF224 + diffF324));
 218     const scalar flux5 =
 219         0.5*(fluxLeft15 + fluxRight15 - (diffF15 + diffF25 + diffF35));
 220
 221     // Compute private data
 222     rhoFlux  = flux1*magSf;
 223     rhoUFlux = flux24*magSf;
 224     rhoEFlux = flux5*magSf;
 225 }
 226
 227 // ************************************************************************* //