FreeFOAM The Cross-Platform CFD Toolkit
rhoCentralDyMFoam.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-2009 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  - rhoCentralDyMFoam [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>
57 
58 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
59 
60 int main(int argc, char *argv[])
61 {
62  #include <OpenFOAM/setRootCase.H>
63 
64  #include <OpenFOAM/createTime.H>
65  #include <OpenFOAM/createMesh.H>
66  #include "createFields.H"
69 
70  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
71 
72  #include "../readFluxScheme.H"
73 
74  dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
75 
76  Info<< "\nStarting time loop\n" << endl;
77 
78  autoPtr<Foam::motionSolver> motionPtr = motionSolver::New(mesh);
79 
80  while (runTime.run())
81  {
82  // --- upwind interpolation of primitive fields on faces
83 
84  surfaceScalarField rho_pos =
85  fvc::interpolate(rho, pos, "reconstruct(rho)");
86  surfaceScalarField rho_neg =
87  fvc::interpolate(rho, neg, "reconstruct(rho)");
88 
89  surfaceVectorField rhoU_pos =
90  fvc::interpolate(rhoU, pos, "reconstruct(U)");
91  surfaceVectorField rhoU_neg =
92  fvc::interpolate(rhoU, neg, "reconstruct(U)");
93 
94  volScalarField rPsi = 1.0/psi;
95  surfaceScalarField rPsi_pos =
96  fvc::interpolate(rPsi, pos, "reconstruct(T)");
97  surfaceScalarField rPsi_neg =
98  fvc::interpolate(rPsi, neg, "reconstruct(T)");
99 
100  surfaceScalarField e_pos =
101  fvc::interpolate(e, pos, "reconstruct(T)");
102  surfaceScalarField e_neg =
103  fvc::interpolate(e, neg, "reconstruct(T)");
104 
105  surfaceVectorField U_pos = rhoU_pos/rho_pos;
106  surfaceVectorField U_neg = rhoU_neg/rho_neg;
107 
108  surfaceScalarField p_pos = rho_pos*rPsi_pos;
109  surfaceScalarField p_neg = rho_neg*rPsi_neg;
110 
111  surfaceScalarField phiv_pos = U_pos & mesh.Sf();
112  surfaceScalarField phiv_neg = U_neg & mesh.Sf();
113 
114  volScalarField c = sqrt(thermo.Cp()/thermo.Cv()*rPsi);
115  surfaceScalarField cSf_pos =
116  fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf();
117  surfaceScalarField cSf_neg =
118  fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf();
119 
120  surfaceScalarField ap =
121  max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero);
122  surfaceScalarField am =
123  min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero);
124 
125  surfaceScalarField a_pos = ap/(ap - am);
126 
127  surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
128 
129  surfaceScalarField aSf = am*a_pos;
130 
131  if (fluxScheme == "Tadmor")
132  {
133  aSf = -0.5*amaxSf;
134  a_pos = 0.5;
135  }
136 
137  surfaceScalarField a_neg = (1.0 - a_pos);
138 
139  phiv_pos *= a_pos;
140  phiv_neg *= a_neg;
141 
142  surfaceScalarField aphiv_pos = phiv_pos - aSf;
143  surfaceScalarField aphiv_neg = phiv_neg + aSf;
144 
145  // Reuse amaxSf for the maximum positive and negative fluxes
146  // estimated by the central scheme
147  amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
148 
149  #include "../compressibleCourantNo.H"
151  #include <finiteVolume/setDeltaT.H>
152 
153  runTime++;
154 
155  Info<< "Time = " << runTime.timeName() << nl << endl;
156 
157  mesh.movePoints(motionPtr->newPoints());
158  phiv_pos = U_pos & mesh.Sf();
159  phiv_neg = U_neg & mesh.Sf();
160  fvc::makeRelative(phiv_pos, U);
161  fvc::makeRelative(phiv_neg, U);
162  phiv_neg -= mesh.phi();
163  phiv_pos *= a_pos;
164  phiv_neg *= a_neg;
165  aphiv_pos = phiv_pos - aSf;
166  aphiv_neg = phiv_neg + aSf;
167 
168  phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
169 
170  surfaceVectorField phiUp =
171  (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
172  + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf();
173 
174  surfaceScalarField phiEp =
175  aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
176  + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
177  + aSf*p_pos - aSf*p_neg;
178 
179  volScalarField muEff = turbulence->muEff();
180  volTensorField tauMC("tauMC", muEff*dev2(fvc::grad(U)().T()));
181 
182  // --- Solve density
184 
185  // --- Solve momentum
186  solve(fvm::ddt(rhoU) + fvc::div(phiUp));
187 
188  U.dimensionedInternalField() =
189  rhoU.dimensionedInternalField()
190  /rho.dimensionedInternalField();
191  U.correctBoundaryConditions();
192  rhoU.boundaryField() = rho.boundaryField()*U.boundaryField();
193 
194  if (!inviscid)
195  {
196  solve
197  (
198  fvm::ddt(rho, U) - fvc::ddt(rho, U)
199  - fvm::laplacian(muEff, U)
200  - fvc::div(tauMC)
201  );
202  rhoU = rho*U;
203  }
204 
205  // --- Solve energy
206  surfaceScalarField sigmaDotU =
207  (
208  (
210  + (mesh.Sf() & fvc::interpolate(tauMC))
211  )
212  & (a_pos*U_pos + a_neg*U_neg)
213  );
214 
215  solve
216  (
217  fvm::ddt(rhoE)
218  + fvc::div(phiEp)
219  - fvc::div(sigmaDotU)
220  );
221 
222  e = rhoE/rho - 0.5*magSqr(U);
224  thermo.correct();
225  rhoE.boundaryField() =
226  rho.boundaryField()*
227  (
228  e.boundaryField() + 0.5*magSqr(U.boundaryField())
229  );
230 
231  if (!inviscid)
232  {
233  volScalarField k("k", thermo.Cp()*muEff/Pr);
234  solve
235  (
236  fvm::ddt(rho, e) - fvc::ddt(rho, e)
237  - fvm::laplacian(turbulence->alphaEff(), e)
238  + fvc::laplacian(turbulence->alpha(), e)
239  - fvc::laplacian(k, T)
240  );
241  thermo.correct();
242  rhoE = rho*(e + 0.5*magSqr(U));
243  }
244 
245  p.dimensionedInternalField() =
246  rho.dimensionedInternalField()
248  p.correctBoundaryConditions();
249  rho.boundaryField() = psi.boundaryField()*p.boundaryField();
250 
251  turbulence->correct();
252 
253  runTime.write();
254 
255  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
256  << " ClockTime = " << runTime.elapsedClockTime() << " s"
257  << nl << endl;
258  }
259 
260  Info<< "End\n" << endl;
261 
262  return 0;
263 }
264 
265 // ************************ vim: set sw=4 sts=4 et: ************************ //