src/dbns/dbnsFlux/hllcFlux/hllcFlux.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 "hllcFlux.H"
  27
  28 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  29
  30 void Foam::hllcFlux::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 left and right:
  50     // normal vector
  51     const vector normalVector = Sf/magSf;
  52
  53     // Ratio of specific heat capacities
  54     const scalar kappaLeft = (RLeft + CvLeft)/CvLeft;
  55     const scalar kappaRight = (RRight + CvRight)/CvRight;
  56
  57     // Compute conservative variables assuming perfect gas law
  58
  59     // Density
  60     const scalar rhoLeft = pLeft/(RLeft*TLeft);
  61     const scalar rhoRight = pRight/(RRight*TRight);
  62
  63     // DensityVelocity
  64     const vector rhoULeft = rhoLeft*ULeft;
  65     const vector rhoURight = rhoRight*URight;
  66
  67     // DensityTotalEnergy
  68     const scalar rhoELeft = rhoLeft*(CvLeft*TLeft+0.5*magSqr(ULeft));
  69     const scalar rhoERight = rhoRight*(CvRight*TRight+0.5*magSqr(URight));
  70
  71     // Compute left and right total enthalpies:
  72     const scalar HLeft = (rhoELeft + pLeft)/rhoLeft;
  73     const scalar HRight = (rhoERight + pRight)/rhoRight;
  74
  75     // Compute qLeft and qRight (q_{l,r} = U_{l,r} \bullet n)
  76     const scalar qLeft = (ULeft & normalVector);
  77     const scalar qRight = (URight & normalVector);
  78
  79     // Speed of sound, for left and right side, assuming perfect gas
  80     const scalar aLeft =
  81         Foam::sqrt(max(0.0,kappaLeft * pLeft/rhoLeft));
  82
  83     const scalar aRight =
  84         Foam::sqrt(max(0.0,kappaRight * pRight/rhoRight));
  85
  86
  87     // Step 2:
  88     // needs rho_{l,r}, U_{l,r}, H_{l,r}, kappa_{l,r}, Gamma_{l,r}, q_{l,r}
  89
  90     // Compute Roe weights
  91     const scalar rhoLeftSqrt = Foam::sqrt(max(0.0,rhoLeft));
  92     const scalar rhoRightSqrt = Foam::sqrt(max(0.0,rhoRight));
  93
  94     const scalar wLeft = rhoLeftSqrt
  95         /stabilise((rhoLeftSqrt + rhoRightSqrt),VSMALL);
  96
  97     const scalar wRight = 1 - wLeft;
  98
  99     // Roe averaged velocity
 100     const vector UTilde = wLeft*ULeft + wRight*URight;
 101
 102     // Roe averaged contravariant velocity
 103     const scalar contrUTilde = (UTilde & normalVector);
 104
 105     // Roe averaged total enthalpy
 106     const scalar HTilde = wLeft*HLeft + wRight*HRight;
 107
 108     // Roe averaged kappa
 109     // TODO: needs to be verified!
 110     const scalar kappaTilde = wLeft*kappaLeft + wRight*kappaRight;
 111
 112     // Speed of sound with Roe reconstruction values
 113     // TODO: not sure if the correct (flow speed) and kappa is used here
 114     const scalar aTilde =
 115         Foam::sqrt(max(0 ,(kappaTilde - 1)*(HTilde - 0.5*sqr(contrUTilde))));
 116
 117     // Step 3: compute signal speeds for face:
 118     const scalar SLeft  = min(qLeft-aLeft, contrUTilde-aTilde);
 119     const scalar SRight = max(contrUTilde+aTilde, qRight+aRight);
 120
 121     const scalar SStar = (rhoRight*qRight*(SRight-qRight)
 122       - rhoLeft*qLeft*(SLeft - qLeft) + pLeft - pRight )/
 123         stabilise((rhoRight*(SRight-qRight)-rhoLeft*(SLeft-qLeft)),VSMALL);
 124
 125     // Compute pressure in star region from the right side
 126     const scalar pStarRight =
 127         rhoRight*(qRight - SRight)*(qRight - SStar) + pRight;
 128
 129     // Should be equal to the left side
 130     const scalar pStarLeft  =
 131         rhoLeft*(qLeft -  SLeft)*(qLeft - SStar) + pLeft;
 132
 133     // Give a warning if this is not the case
 134     if (mag(pStarRight - pStarLeft) > 1e-6)
 135     {
 136         Info << "mag(pStarRight-pStarLeft) > VSMALL " << endl;
 137     }
 138
 139     // Use pStarRight for pStar, as in theory, pStarRight == pStarLeft
 140     const scalar pStar = pStarRight;
 141
 142     // Step 4: upwinding - compute states:
 143     scalar convectionSpeed = 0.0;
 144     scalar rhoState = 0.0;
 145     vector rhoUState = vector::zero;
 146     scalar rhoEState = 0.0;
 147     scalar pState = 0.0;
 148
 149     if (pos(SLeft))
 150     {
 151         // compute F_l
 152         convectionSpeed = qLeft;
 153         rhoState  = rhoLeft;
 154         rhoUState = rhoULeft;
 155         rhoEState = rhoELeft;
 156         pState = pLeft;
 157     }
 158     else if (pos(SStar))
 159     {
 160         scalar omegaLeft = scalar(1.0)/stabilise((SLeft - SStar), VSMALL);
 161
 162         // Compute left star region
 163         convectionSpeed = SStar;
 164         rhoState  = omegaLeft*(SLeft - qLeft)*rhoLeft;
 165         rhoUState = omegaLeft*((SLeft - qLeft)*rhoULeft
 166         + (pStar - pLeft)*normalVector);
 167         rhoEState = omegaLeft*((SLeft - qLeft)*rhoELeft
 168         - pLeft*qLeft + pStar*SStar);
 169         pState = pStar;
 170     }
 171     else if (pos(SRight))
 172     {
 173         scalar omegaRight = scalar(1.0)/stabilise((SRight - SStar), VSMALL);
 174
 175         // compute right star region
 176         convectionSpeed = SStar;
 177         rhoState  = omegaRight*(SRight - qRight)*rhoRight;
 178         rhoUState = omegaRight*((SRight - qRight)*rhoURight
 179         + (pStar - pRight)*normalVector);
 180         rhoEState = omegaRight*((SRight - qRight)*rhoERight
 181         - pRight*qRight + pStar*SStar);
 182         pState = pStar;
 183     }
 184     else if (neg(SRight))
 185     {
 186         // compute F_r
 187         convectionSpeed = qRight;
 188         rhoState  = rhoRight;
 189         rhoUState = rhoURight;
 190         rhoEState = rhoERight;
 191         pState = pRight;
 192     }
 193     else
 194     {
 195         Info << "Error in HLLC Riemann solver" << endl;
 196     }
 197
 198     rhoFlux  = (convectionSpeed*rhoState)*magSf;
 199     rhoUFlux = (convectionSpeed*rhoUState+pState*normalVector)*magSf;
 200     rhoEFlux = (convectionSpeed*(rhoEState+pState))*magSf;
 201 }
 202
 203
 204 // ************************************************************************* //