denseParticleFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2013-2022 OpenFOAM Foundation
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  denseParticleFoam
26 
27 Description
28  Transient solver for the coupled transport of particle clouds including the
29  effect of the volume fraction of particles on the continuous phase, with
30  optional mesh motion and mesh topology changes.
31 
32 \*---------------------------------------------------------------------------*/
33 
34 #include "NamedEnum.H"
35 
36 namespace Foam
37 {
38  enum class cloudForceSplit
39  {
40  faceExplicitCellImplicit, // Implicit part of the cloud force added to
41  // the cell momentum equation. Explicit part
42  // to the face momentum equation. This is the
43  // least likely to create staggering patterns
44  // in the velocity field, but it can create
45  // unphysical perturbations in cell
46  // velocities even when particles and flow
47  // have the similar velocities.
48 
49  faceExplicitCellLagged, // Entire cloud force evaluated explicitly
50  // and added to the face momentum equation.
51  // Lagged correction (i.e.,
52  // fvm::Sp(cloudSU.diag(), Uc) -
53  // cloudSU.diag()*Uc) added to the cell
54  // momentum equation. This creates physical
55  // cell velocities when particles and flow
56  // have the same velocity, but can also
57  // result in staggering patterns in packed
58  // beds. Unsuitable for MPPIC.
59 
60  faceImplicit // Implicit and explicit parts of the force
61  // both added to the face momentum equation.
62  // Behaves somewhere between the other two.
63  };
64 
65  template<>
67  {
68  "faceExplicitCellImplicit",
69  "faceExplicitCellLagged",
70  "faceImplicit"
71  };
72 
73  const NamedEnum<cloudForceSplit, 3> cloudForceSplitNames;
74 }
75 
76 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
77 
78 #include "fvCFD.H"
79 #include "viscosityModel.H"
81 #include "pimpleControl.H"
82 #include "pressureReference.H"
83 #include "CorrectPhi.H"
84 #include "fvModels.H"
85 #include "fvConstraints.H"
86 #include "parcelCloudList.H"
87 
88 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
89 
90 int main(int argc, char *argv[])
91 {
92  #include "postProcess.H"
93 
94  #include "setRootCaseLists.H"
95  #include "createTime.H"
96  #include "createMesh.H"
97  #include "createDyMControls.H"
98  #include "createFields.H"
99  #include "createUcfIfPresent.H"
100  #include "initContinuityErrs.H"
101 
102  Info<< "\nStarting time loop\n" << endl;
103 
104  while (pimple.run(runTime))
105  {
106  #include "readDyMControls.H"
107  #include "CourantNo.H"
108  #include "setDeltaT.H"
109 
110  // Update the mesh for topology change, mesh to mesh mapping
111  mesh.update();
112 
113  runTime++;
114 
115  Info<< "Time = " << runTime.userTimeName() << nl << endl;
116 
117  // Store the particle positions
118  clouds.storeGlobalPositions();
119 
120  // Move the mesh
121  mesh.move();
122 
123  if (mesh.changing())
124  {
125  if (correctPhi)
126  {
127  #include "correctPhic.H"
128  }
129 
130  if (checkMeshCourantNo)
131  {
132  #include "meshCourantNo.H"
133  }
134  }
135 
136  continuousPhaseViscosity->correct();
138 
139  clouds.evolve();
140 
141  // Update continuous phase volume fraction field
142  alphac = max(1.0 - clouds.theta(), alphacMin);
143  alphac.correctBoundaryConditions();
144  alphacf = fvc::interpolate(alphac);
145  alphaPhic = alphacf*phic;
146 
147  // Cloud forces
148  fvVectorMatrix cloudSU(clouds.SU(Uc));
149  volVectorField cloudSUu
150  (
151  IOobject
152  (
153  "cloudSUu",
154  runTime.timeName(),
155  mesh
156  ),
157  mesh,
159  zeroGradientFvPatchVectorField::typeName
160  );
161  volScalarField cloudSUp
162  (
163  IOobject
164  (
165  "cloudSUp",
166  runTime.timeName(),
167  mesh
168  ),
169  mesh,
171  zeroGradientFvPatchVectorField::typeName
172  );
173 
174  const cloudForceSplit cloudSUSplit =
175  pimple.dict().found("cloudForceSplit")
176  ? cloudForceSplitNames.read(pimple.dict().lookup("cloudForceSplit"))
177  : cloudForceSplit::faceExplicitCellImplicit;
178 
179  switch (cloudSUSplit)
180  {
181  case cloudForceSplit::faceExplicitCellImplicit:
182  cloudSUu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
183  cloudSUu.correctBoundaryConditions();
184  cloudSUp.primitiveFieldRef() = Zero;
185  cloudSUp.correctBoundaryConditions();
186  //cloudSU.diag() = cloudSU.diag();
187  cloudSU.source() = Zero;
188  break;
189 
190  case cloudForceSplit::faceExplicitCellLagged:
191  cloudSUu.primitiveFieldRef() =
192  (cloudSU.diag()*Uc() - cloudSU.source())/mesh.V();
193  cloudSUu.correctBoundaryConditions();
194  cloudSUp.primitiveFieldRef() = Zero;
195  cloudSUp.correctBoundaryConditions();
196  //cloudSU.diag() = cloudSU.diag();
197  cloudSU.source() = cloudSU.diag()*Uc();
198  break;
199 
200  case cloudForceSplit::faceImplicit:
201  cloudSUu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
202  cloudSUu.correctBoundaryConditions();
203  cloudSUp.primitiveFieldRef() = cloudSU.diag()/mesh.V();
204  cloudSUp.correctBoundaryConditions();
205  cloudSU.diag() = Zero;
206  cloudSU.source() = Zero;
207  break;
208  }
209 
210  // --- Pressure-velocity PIMPLE corrector loop
211  while (pimple.loop())
212  {
213  fvModels.correct();
214 
215  #include "UcEqn.H"
216 
217  // --- PISO loop
218  while (pimple.correct())
219  {
220  #include "pEqn.H"
221  }
222 
223  if (pimple.turbCorr())
224  {
225  continuousPhaseTurbulence->correct();
226  }
227  }
228 
229  runTime.write();
230 
231  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
232  << " ClockTime = " << runTime.elapsedClockTime() << " s"
233  << nl << endl;
234  }
235 
236  Info<< "End\n" << endl;
237 
238  return 0;
239 }
240 
241 
242 // ************************************************************************* //
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
pimpleNoLoopControl & pimple
volScalarField muc(IOobject(IOobject::groupName("mu", continuousPhaseName), runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), rhoc *continuousPhaseViscosity->nu())
Info<< "Reading field U\"<< endl;volVectorField Uc(IOobject(IOobject::groupName("U", continuousPhaseName), runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Reading field p\"<< endl;volScalarField p(IOobject("p", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Reading/calculating continuous-phase face flux field phic\"<< endl;surfaceScalarField phic(IOobject(IOobject::groupName("phi", continuousPhaseName), runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), linearInterpolate(Uc) &mesh.Sf());pressureReference pressureReference(p, pimple.dict());mesh.schemes().setFluxRequired(p.name());Info<< "Creating turbulence model\"<< endl;autoPtr< viscosityModel > continuousPhaseViscosity(viscosityModel::New(mesh))
virtual void correct()
Correct the fvModels.
Definition: fvModels.C:353
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
correctPhi
checkMeshCourantNo
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:59
fvMesh & mesh
const char * names[]
Definition: EDC.C:32
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:58
Creates and initialises the continuous phase face velocity field Ufc if required. ...
phic
Definition: correctPhic.H:2
static const zero Zero
Definition: zero.H:97
const dimensionSet dimForce
volScalarField rhoc(IOobject(rhocValue.name(), runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, rhocValue)
static const char nl
Definition: Ostream.H:260
const dimensionSet dimVelocity
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Foam::fvModels & fvModels
Calculates and outputs the mean and maximum Courant Numbers.
Info<< "Creating field alphac\"<< endl;volScalarField alphac(IOobject(IOobject::groupName("alpha", continuousPhaseName), runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), mesh, dimensionedScalar(dimless, 0));Info<< "Constructing clouds"<< endl;parcelCloudList clouds(rhoc, Uc, muc, g);scalar alphacMin(1 - mesh.solution().solverDict(alphac.name()).lookup< scalar >"max"));alphac=max(1.0 - clouds.theta(), alphacMin);alphac.correctBoundaryConditions();surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));surfaceScalarField alphaPhic(IOobject::groupName("alphaPhi", continuousPhaseName), alphacf *phic);autoPtr< phaseIncompressible::momentumTransportModel > continuousPhaseTurbulence(phaseIncompressible::momentumTransportModel::New(alphac, Uc, alphaPhic, phic, continuousPhaseViscosity))
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
const dimensionSet dimVolume
messageStream Info
Execute application functionObjects to post-process existing results.
Namespace for OpenFOAM.