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