applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     rhoCentralFoam
  26
  27 Description
  28     Density-based compressible flow solver based on central-upwind schemes of
  29     Kurganov and Tadmor
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "fvCFD.H"
  34 #include "basicPsiThermo.H"
  35 #include "turbulenceModel.H"
  36 #include "zeroGradientFvPatchFields.H"
  37 #include "fixedRhoFvPatchScalarField.H"
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 int main(int argc, char *argv[])
  42 {
  43     #include "setRootCase.H"
  44
  45     #include "createTime.H"
  46     #include "createMesh.H"
  47     #include "createFields.H"
  48     #include "readThermophysicalProperties.H"
  49     #include "readTimeControls.H"
  50
  51     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  52
  53     #include "readFluxScheme.H"
  54
  55     dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
  56
  57     Info<< "\nStarting time loop\n" << endl;
  58
  59     while (runTime.run())
  60     {
  61         // --- upwind interpolation of primitive fields on faces
  62
  63         surfaceScalarField rho_pos
  64         (
  65             fvc::interpolate(rho, pos, "reconstruct(rho)")
  66         );
  67         surfaceScalarField rho_neg
  68         (
  69             fvc::interpolate(rho, neg, "reconstruct(rho)")
  70         );
  71
  72         surfaceVectorField rhoU_pos
  73         (
  74             fvc::interpolate(rhoU, pos, "reconstruct(U)")
  75         );
  76         surfaceVectorField rhoU_neg
  77         (
  78             fvc::interpolate(rhoU, neg, "reconstruct(U)")
  79         );
  80
  81         volScalarField rPsi(1.0/psi);
  82         surfaceScalarField rPsi_pos
  83         (
  84             fvc::interpolate(rPsi, pos, "reconstruct(T)")
  85         );
  86         surfaceScalarField rPsi_neg
  87         (
  88             fvc::interpolate(rPsi, neg, "reconstruct(T)")
  89         );
  90
  91         surfaceScalarField e_pos
  92         (
  93             fvc::interpolate(e, pos, "reconstruct(T)")
  94         );
  95         surfaceScalarField e_neg
  96         (
  97             fvc::interpolate(e, neg, "reconstruct(T)")
  98         );
  99
 100         surfaceVectorField U_pos(rhoU_pos/rho_pos);
 101         surfaceVectorField U_neg(rhoU_neg/rho_neg);
 102
 103         surfaceScalarField p_pos(rho_pos*rPsi_pos);
 104         surfaceScalarField p_neg(rho_neg*rPsi_neg);
 105
 106         surfaceScalarField phiv_pos(U_pos & mesh.Sf());
 107         surfaceScalarField phiv_neg(U_neg & mesh.Sf());
 108
 109         volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
 110         surfaceScalarField cSf_pos
 111         (
 112             fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf()
 113         );
 114         surfaceScalarField cSf_neg
 115         (
 116             fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf()
 117         );
 118
 119         surfaceScalarField ap
 120         (
 121             max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
 122         );
 123         surfaceScalarField am
 124         (
 125             min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
 126         );
 127
 128         surfaceScalarField a_pos(ap/(ap - am));
 129
 130         surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
 131
 132         surfaceScalarField aSf(am*a_pos);
 133
 134         if (fluxScheme == "Tadmor")
 135         {
 136             aSf = -0.5*amaxSf;
 137             a_pos = 0.5;
 138         }
 139
 140         surfaceScalarField a_neg(1.0 - a_pos);
 141
 142         phiv_pos *= a_pos;
 143         phiv_neg *= a_neg;
 144
 145         surfaceScalarField aphiv_pos(phiv_pos - aSf);
 146         surfaceScalarField aphiv_neg(phiv_neg + aSf);
 147
 148         // Reuse amaxSf for the maximum positive and negative fluxes
 149         // estimated by the central scheme
 150         amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
 151
 152         #include "compressibleCourantNo.H"
 153         #include "readTimeControls.H"
 154         #include "setDeltaT.H"
 155
 156         runTime++;
 157
 158         Info<< "Time = " << runTime.timeName() << nl << endl;
 159
 160         surfaceScalarField phi("phi", aphiv_pos*rho_pos + aphiv_neg*rho_neg);
 161
 162         surfaceVectorField phiUp
 163         (
 164             (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
 165           + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
 166         );
 167
 168         surfaceScalarField phiEp
 169         (
 170             aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
 171           + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
 172           + aSf*p_pos - aSf*p_neg
 173         );
 174
 175         volScalarField muEff(turbulence->muEff());
 176         volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
 177
 178         // --- Solve density
 179         solve(fvm::ddt(rho) + fvc::div(phi));
 180
 181         // --- Solve momentum
 182         solve(fvm::ddt(rhoU) + fvc::div(phiUp));
 183
 184         U.dimensionedInternalField() =
 185             rhoU.dimensionedInternalField()
 186            /rho.dimensionedInternalField();
 187         U.correctBoundaryConditions();
 188         rhoU.boundaryField() = rho.boundaryField()*U.boundaryField();
 189
 190         volScalarField rhoBydt(rho/runTime.deltaT());
 191
 192         if (!inviscid)
 193         {
 194             solve
 195             (
 196                 fvm::ddt(rho, U) - fvc::ddt(rho, U)
 197               - fvm::laplacian(muEff, U)
 198               - fvc::div(tauMC)
 199             );
 200             rhoU = rho*U;
 201         }
 202
 203         // --- Solve energy
 204         surfaceScalarField sigmaDotU
 205         (
 206             (
 207                 fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
 208               + (mesh.Sf() & fvc::interpolate(tauMC))
 209             )
 210             & (a_pos*U_pos + a_neg*U_neg)
 211         );
 212
 213         solve
 214         (
 215             fvm::ddt(rhoE)
 216           + fvc::div(phiEp)
 217           - fvc::div(sigmaDotU)
 218         );
 219
 220         e = rhoE/rho - 0.5*magSqr(U);
 221         e.correctBoundaryConditions();
 222         thermo.correct();
 223         rhoE.boundaryField() =
 224             rho.boundaryField()*
 225             (
 226                 e.boundaryField() + 0.5*magSqr(U.boundaryField())
 227             );
 228
 229         if (!inviscid)
 230         {
 231             volScalarField k("k", thermo.Cp()*muEff/Pr);
 232             solve
 233             (
 234                 fvm::ddt(rho, e) - fvc::ddt(rho, e)
 235               - fvm::laplacian(turbulence->alphaEff(), e)
 236               + fvc::laplacian(turbulence->alpha(), e)
 237               - fvc::laplacian(k, T)
 238             );
 239             thermo.correct();
 240             rhoE = rho*(e + 0.5*magSqr(U));
 241         }
 242
 243         p.dimensionedInternalField() =
 244             rho.dimensionedInternalField()
 245            /psi.dimensionedInternalField();
 246         p.correctBoundaryConditions();
 247         rho.boundaryField() = psi.boundaryField()*p.boundaryField();
 248
 249         turbulence->correct();
 250
 251         runTime.write();
 252
 253         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 254             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 255             << nl << endl;
 256     }
 257
 258     Info<< "End\n" << endl;
 259
 260     return 0;
 261 }
 262
 263 // ************************************************************************* //