FreeFOAM The Cross-Platform CFD Toolkit
rhoCentralFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 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 Usage
32  - rhoCentralFoam [OPTION]
33 
34  @param -case <dir> \n
35  Specify the case directory
36 
37  @param -parallel \n
38  Run the case in parallel
39 
40  @param -help \n
41  Display short usage message
42 
43  @param -doc \n
44  Display Doxygen documentation page
45 
46  @param -srcDoc \n
47  Display source code
48 
49 \*---------------------------------------------------------------------------*/
50 
51 #include <finiteVolume/fvCFD.H>
55 
56 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
57 
58 int main(int argc, char *argv[])
59 {
60  #include <OpenFOAM/setRootCase.H>
61 
62  #include <OpenFOAM/createTime.H>
63  #include <OpenFOAM/createMesh.H>
64  #include "createFields.H"
67 
68  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
69 
70  #include "readFluxScheme.H"
71 
72  dimensionedScalar v_zero("v_zero",dimVolume/dimTime, 0.0);
73 
74  Info<< "\nStarting time loop\n" << endl;
75 
76  while (runTime.run())
77  {
78  // --- upwind interpolation of primitive fields on faces
79 
80  surfaceScalarField rho_pos =
81  fvc::interpolate(rho, pos, "reconstruct(rho)");
82  surfaceScalarField rho_neg =
83  fvc::interpolate(rho, neg, "reconstruct(rho)");
84 
85  surfaceVectorField rhoU_pos =
86  fvc::interpolate(rhoU, pos, "reconstruct(U)");
87  surfaceVectorField rhoU_neg =
88  fvc::interpolate(rhoU, neg, "reconstruct(U)");
89 
90  volScalarField rPsi = 1.0/psi;
91  surfaceScalarField rPsi_pos =
92  fvc::interpolate(rPsi, pos, "reconstruct(T)");
93  surfaceScalarField rPsi_neg =
94  fvc::interpolate(rPsi, neg, "reconstruct(T)");
95 
96  surfaceScalarField e_pos =
97  fvc::interpolate(e, pos, "reconstruct(T)");
98  surfaceScalarField e_neg =
99  fvc::interpolate(e, neg, "reconstruct(T)");
100 
101  surfaceVectorField U_pos = rhoU_pos/rho_pos;
102  surfaceVectorField U_neg = rhoU_neg/rho_neg;
103 
104  surfaceScalarField p_pos = rho_pos*rPsi_pos;
105  surfaceScalarField p_neg = rho_neg*rPsi_neg;
106 
107  surfaceScalarField phiv_pos = U_pos & mesh.Sf();
108  surfaceScalarField phiv_neg = U_neg & mesh.Sf();
109 
110  volScalarField c = sqrt(thermo.Cp()/thermo.Cv()*rPsi);
111  surfaceScalarField cSf_pos = fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf();
112  surfaceScalarField cSf_neg = fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf();
113 
114  surfaceScalarField ap = max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero);
115  surfaceScalarField am = min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero);
116 
117  surfaceScalarField a_pos = ap/(ap - am);
118 
119  surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
120 
121  surfaceScalarField aSf = am*a_pos;
122 
123  if (fluxScheme == "Tadmor")
124  {
125  aSf = -0.5*amaxSf;
126  a_pos = 0.5;
127  }
128 
129  surfaceScalarField a_neg = (1.0 - a_pos);
130 
131  phiv_pos *= a_pos;
132  phiv_neg *= a_neg;
133 
134  surfaceScalarField aphiv_pos = phiv_pos - aSf;
135  surfaceScalarField aphiv_neg = phiv_neg + aSf;
136 
137  // Reuse amaxSf for the maximum positive and negative fluxes
138  // estimated by the central scheme
139  amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
140 
141  #include "compressibleCourantNo.H"
143  #include <finiteVolume/setDeltaT.H>
144 
145  runTime++;
146 
147  Info<< "Time = " << runTime.timeName() << nl << endl;
148 
149  surfaceScalarField phi("phi", aphiv_pos*rho_pos + aphiv_neg*rho_neg);
150 
151  surfaceVectorField phiUp =
152  (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
153  + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf();
154 
155  surfaceScalarField phiEp =
156  aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
157  + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
158  + aSf*p_pos - aSf*p_neg;
159 
160  volTensorField tauMC("tauMC", mu*dev2(fvc::grad(U)().T()));
161 
162  // --- Solve density
164 
165  // --- Solve momentum
166  solve(fvm::ddt(rhoU) + fvc::div(phiUp));
167 
168  U.dimensionedInternalField() =
169  rhoU.dimensionedInternalField()
170  /rho.dimensionedInternalField();
171  U.correctBoundaryConditions();
172  rhoU.boundaryField() = rho.boundaryField()*U.boundaryField();
173 
174  volScalarField rhoBydt(rho/runTime.deltaT());
175 
176  if (!inviscid)
177  {
178  solve
179  (
180  fvm::ddt(rho, U) - fvc::ddt(rho, U)
181  - fvm::laplacian(mu, U)
182  - fvc::div(tauMC)
183  );
184  rhoU = rho*U;
185  }
186 
187  // --- Solve energy
188  surfaceScalarField sigmaDotU =
189  (
190  (
192  + (mesh.Sf() & fvc::interpolate(tauMC))
193  )
194  & (a_pos*U_pos + a_neg*U_neg)
195  );
196 
197  solve
198  (
199  fvm::ddt(rhoE)
200  + fvc::div(phiEp)
201  - fvc::div(sigmaDotU)
202  );
203 
204  e = rhoE/rho - 0.5*magSqr(U);
206  thermo.correct();
207  rhoE.boundaryField() =
208  rho.boundaryField()*
209  (
210  e.boundaryField() + 0.5*magSqr(U.boundaryField())
211  );
212 
213  if (!inviscid)
214  {
215  volScalarField k("k", thermo.Cp()*mu/Pr);
216  solve
217  (
218  fvm::ddt(rho, e) - fvc::ddt(rho, e)
221  - fvc::laplacian(k, T)
222  );
223  thermo.correct();
224  rhoE = rho*(e + 0.5*magSqr(U));
225  }
226 
227  p.dimensionedInternalField() =
228  rho.dimensionedInternalField()
230  p.correctBoundaryConditions();
231  rho.boundaryField() = psi.boundaryField()*p.boundaryField();
232 
233  runTime.write();
234 
235  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
236  << " ClockTime = " << runTime.elapsedClockTime() << " s"
237  << nl << endl;
238  }
239 
240  Info<< "End\n" << endl;
241 
242  return 0;
243 }
244 
245 // ************************ vim: set sw=4 sts=4 et: ************************ //