src/dbns/dbnsFlux/rusanovFlux/rusanovFlux.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 "rusanovFlux.H"
  27
  28 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  29
  30 void Foam::rusanovFlux::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     // Step 2: compute Roe averged quantities for face:
  75     const scalar rhoTilde = sqrt(max(rhoLeft*rhoRight, SMALL));
  76
  77     // Some temporary variables:
  78     const scalar rhoLeftSqrt = sqrt(max(rhoLeft, SMALL));
  79     const scalar rhoRightSqrt = sqrt(max(rhoRight, SMALL));
  80
  81     const scalar wLeft = rhoLeftSqrt/(rhoLeftSqrt + rhoRightSqrt);
  82     const scalar wRight = 1 - wLeft;
  83
  84     const vector UTilde = ULeft*wLeft + URight*wRight;
  85     const scalar hTilde = hLeft*wLeft + hRight*wRight;
  86     const scalar qTildeSquare = magSqr(UTilde);
  87     const scalar kappaTilde = kappaLeft*wLeft + kappaRight*wRight;
  88
  89     // Speed of sound
  90     const scalar cTilde =
  91         sqrt(max((kappaTilde - 1)*(hTilde - 0.5*qTildeSquare), SMALL));
  92
  93     // Roe averaged contravariant velocity
  94     const scalar contrVTilde = (UTilde & normalVector);
  95
  96     // Step 3: compute primitive differences:
  97     const scalar deltaP = pRight - pLeft;
  98     const scalar deltaRho = rhoRight - rhoLeft;
  99     const vector deltaU = URight - ULeft;
 100     const scalar deltaContrV = (deltaU & normalVector);
 101
 102     // Step 4: compute wave strengths:
 103
 104     // Roe and Pike - formulation
 105     const scalar r1 =
 106         (deltaP - rhoTilde*cTilde*deltaContrV)/(2.0*sqr(cTilde));
 107     const scalar r2 = deltaRho - deltaP/sqr(cTilde);
 108     const scalar r3 =
 109         (deltaP + rhoTilde*cTilde*deltaContrV)/(2.0*sqr(cTilde));
 110
 111     // Step 5: compute l vectors
 112
 113     // rho row:
 114     const scalar l1rho = 1;
 115     const scalar l2rho = 1;
 116     const scalar l3rho = 0;
 117     const scalar l4rho = 1;
 118
 119     // first U column
 120     const vector l1U = UTilde - cTilde*normalVector;
 121
 122     // second U column
 123     const vector l2U = UTilde;
 124
 125     // third U column
 126     const vector l3U = deltaU - deltaContrV*normalVector;
 127
 128     // fourth U column
 129     const vector l4U = UTilde + cTilde*normalVector;
 130
 131     // E row
 132     const scalar l1e = hTilde - cTilde*contrVTilde;
 133     const scalar l2e = 0.5*qTildeSquare;
 134     const scalar l3e = (UTilde & deltaU) - contrVTilde*deltaContrV;
 135     const scalar l4e = hTilde + cTilde*contrVTilde;
 136
 137     // Step 6: compute eigenvalues
 138
 139     // derived from algebra by hand, only for Euler equation usefull
 140     scalar lambda1 = mag(contrVTilde - cTilde);
 141     scalar lambda2 = mag(contrVTilde);
 142     scalar lambda3 = mag(contrVTilde + cTilde);
 143
 144     scalar lambdaMax = max(max(lambda1,lambda2),lambda3);
 145
 146     // Step 7: Compute flux differences
 147
 148     // Components of deltaF1
 149     const scalar diffF11 = lambdaMax*r1*l1rho;
 150     const vector diffF124 = lambdaMax*r1*l1U;
 151     const scalar diffF15 = lambdaMax*r1*l1e;
 152
 153     // Components of deltaF2
 154     const scalar diffF21 = lambdaMax*(r2*l2rho + rhoTilde*l3rho);
 155     const vector diffF224 = lambdaMax*(r2*l2U + rhoTilde*l3U);
 156     const scalar diffF25 = lambdaMax*(r2*l2e + rhoTilde*l3e);
 157
 158     // Components of deltaF3
 159     const scalar diffF31 = lambdaMax*r3*l4rho;
 160     const vector diffF324 = lambdaMax*r3*l4U;
 161     const scalar diffF35 = lambdaMax*r3*l4e;
 162
 163     // Step 8: compute left and right fluxes
 164
 165     // Left flux 5-vector
 166     const scalar fluxLeft11 = rhoLeft*contrVLeft;
 167     const vector fluxLeft124 = ULeft*fluxLeft11 + normalVector*pLeft;
 168     const scalar fluxLeft15 = hLeft*fluxLeft11;
 169
 170     // Right flux 5-vector
 171     const scalar fluxRight11 = rhoRight*contrVRight;
 172     const vector fluxRight124 = URight*fluxRight11 + normalVector*pRight;
 173     const scalar fluxRight15 = hRight*fluxRight11;
 174
 175     // Step 10: compute face flux 5-vector
 176     const scalar flux1 =
 177         0.5*(fluxLeft11 + fluxRight11 - (diffF11 + diffF21 + diffF31));
 178     const vector flux24 =
 179         0.5*(fluxLeft124 + fluxRight124 - (diffF124 + diffF224 + diffF324));
 180     const scalar flux5 =
 181         0.5*(fluxLeft15 + fluxRight15 - (diffF15 + diffF25 + diffF35));
 182
 183     // Compute private data
 184     rhoFlux  = flux1*magSf;
 185     rhoUFlux = flux24*magSf;
 186     rhoEFlux = flux5*magSf;
 187 }
 188
 189 // ************************************************************************* //