solidificationMeltingSource.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) 2014-2018 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 \*---------------------------------------------------------------------------*/
25 
27 #include "fvMatrices.H"
28 #include "basicThermo.H"
33 #include "geometricOneField.H"
34 
35 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  template<>
40  const char* NamedEnum
41  <
43  2
44  >::names[] =
45  {
46  "thermo",
47  "lookup"
48  };
49 
50  namespace fv
51  {
52  defineTypeNameAndDebug(solidificationMeltingSource, 0);
53 
55  (
56  option,
57  solidificationMeltingSource,
59  );
60  }
61 }
62 
65 
66 
67 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
68 
70 Foam::fv::solidificationMeltingSource::Cp() const
71 {
72  switch (mode_)
73  {
74  case mdThermo:
75  {
76  const basicThermo& thermo =
77  mesh_.lookupObject<basicThermo>(basicThermo::dictName);
78 
79  return thermo.Cp();
80  break;
81  }
82  case mdLookup:
83  {
84  if (CpName_ == "CpRef")
85  {
86  scalar CpRef = readScalar(coeffs_.lookup("CpRef"));
87 
88  return volScalarField::New
89  (
90  name_ + ":Cp",
91  mesh_,
93  (
95  CpRef
96  ),
97  extrapolatedCalculatedFvPatchScalarField::typeName
98  );
99  }
100  else
101  {
102  return mesh_.lookupObject<volScalarField>(CpName_);
103  }
104 
105  break;
106  }
107  default:
108  {
110  << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
111  << abort(FatalError);
112  }
113  }
114 
115  return tmp<volScalarField>(nullptr);
116 }
117 
118 
119 Foam::vector Foam::fv::solidificationMeltingSource::g() const
120 {
121  if (mesh_.foundObject<uniformDimensionedVectorField>("g"))
122  {
123  const uniformDimensionedVectorField& value =
124  mesh_.lookupObject<uniformDimensionedVectorField>("g");
125  return value.value();
126  }
127  else
128  {
129  return coeffs_.lookup("g");
130  }
131 }
132 
133 
134 void Foam::fv::solidificationMeltingSource::update(const volScalarField& Cp)
135 {
136  if (curTimeIndex_ == mesh_.time().timeIndex())
137  {
138  return;
139  }
140 
141  if (debug)
142  {
143  Info<< type() << ": " << name_ << " - updating phase indicator" << endl;
144  }
145 
146  // update old time alpha1 field
147  alpha1_.oldTime();
148 
149  const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
150 
151  forAll(cells_, i)
152  {
153  const label celli = cells_[i];
154 
155  const scalar Tc = T[celli];
156  const scalar Cpc = Cp[celli];
157  const scalar alpha1New =
158  alpha1_[celli]
159  + relax_*Cpc
160  *(
161  Tc
162  - max
163  (
164  Tsol_,
165  Tsol_
166  + (Tliq_ - Tsol_)*(alpha1_[celli] - alpha1e_)/(1 - alpha1e_)
167  )
168  )/L_;
169 
170  alpha1_[celli] = max(0, min(alpha1New, 1));
171  deltaT_[i] =
172  Tc
173  - max
174  (
175  Tsol_,
176  Tsol_
177  + (Tliq_ - Tsol_)*(alpha1_[celli] - alpha1e_)/(1 - alpha1e_)
178  );
179  }
180 
181  alpha1_.correctBoundaryConditions();
182 
183  curTimeIndex_ = mesh_.time().timeIndex();
184 }
185 
186 
187 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
188 
190 (
191  const word& sourceName,
192  const word& modelType,
193  const dictionary& dict,
194  const fvMesh& mesh
195 )
196 :
197  cellSetOption(sourceName, modelType, dict, mesh),
198  Tsol_(readScalar(coeffs_.lookup("Tsol"))),
199  Tliq_(coeffs_.lookupOrDefault("Tliq", Tsol_)),
200  alpha1e_(coeffs_.lookupOrDefault("alpha1e", 0.0)),
201  L_(readScalar(coeffs_.lookup("L"))),
202  relax_(coeffs_.lookupOrDefault("relax", 0.9)),
203  mode_(thermoModeTypeNames_.read(coeffs_.lookup("thermoMode"))),
204  rhoRef_(readScalar(coeffs_.lookup("rhoRef"))),
205  TName_(coeffs_.lookupOrDefault<word>("T", "T")),
206  CpName_(coeffs_.lookupOrDefault<word>("Cp", "Cp")),
207  UName_(coeffs_.lookupOrDefault<word>("U", "U")),
208  phiName_(coeffs_.lookupOrDefault<word>("phi", "phi")),
209  Cu_(coeffs_.lookupOrDefault<scalar>("Cu", 100000)),
210  q_(coeffs_.lookupOrDefault("q", 0.001)),
211  beta_(readScalar(coeffs_.lookup("beta"))),
212  alpha1_
213  (
214  IOobject
215  (
216  name_ + ":alpha1",
217  mesh.time().timeName(),
218  mesh,
221  ),
222  mesh,
224  zeroGradientFvPatchScalarField::typeName
225  ),
226  curTimeIndex_(-1),
227  deltaT_(cells_.size(), 0)
228 {
229  fieldNames_.setSize(2);
230  fieldNames_[0] = UName_;
231 
232  switch (mode_)
233  {
234  case mdThermo:
235  {
236  const basicThermo& thermo =
237  mesh_.lookupObject<basicThermo>(basicThermo::dictName);
238 
239  fieldNames_[1] = thermo.he().name();
240  break;
241  }
242  case mdLookup:
243  {
244  fieldNames_[1] = TName_;
245  break;
246  }
247  default:
248  {
250  << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
251  << abort(FatalError);
252  }
253  }
254 
255  applied_.setSize(fieldNames_.size(), false);
256 }
257 
258 
259 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
260 
262 (
263  fvMatrix<scalar>& eqn,
264  const label fieldi
265 )
266 {
267  apply(geometricOneField(), eqn);
268 }
269 
270 
272 (
273  const volScalarField& rho,
274  fvMatrix<scalar>& eqn,
275  const label fieldi
276 )
277 {
278  apply(rho, eqn);
279 }
280 
281 
283 (
284  fvMatrix<vector>& eqn,
285  const label fieldi
286 )
287 {
288  if (debug)
289  {
290  Info<< type() << ": applying source to " << eqn.psi().name() << endl;
291  }
292 
293  const volScalarField Cp(this->Cp());
294 
295  update(Cp);
296 
297  vector g = this->g();
298 
299  scalarField& Sp = eqn.diag();
300  vectorField& Su = eqn.source();
301  const scalarField& V = mesh_.V();
302 
303  forAll(cells_, i)
304  {
305  const label celli = cells_[i];
306 
307  const scalar Vc = V[celli];
308  const scalar alpha1c = alpha1_[celli];
309 
310  const scalar S = -Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_);
311  const vector Sb = rhoRef_*g*beta_*deltaT_[i];
312 
313  Sp[celli] += Vc*S;
314  Su[celli] += Vc*Sb;
315  }
316 }
317 
318 
320 (
321  const volScalarField& rho,
322  fvMatrix<vector>& eqn,
323  const label fieldi
324 )
325 {
326  // Momentum source uses a Boussinesq approximation - redirect
327  addSup(eqn, fieldi);
328 }
329 
330 
331 // ************************************************************************* //
defineTypeNameAndDebug(option, 0)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const word & name() const
Return name.
Definition: IOobject.H:295
Abstract base-class for fluid and solid thermodynamic properties.
Definition: basicThermo.H:52
zeroField Su
Definition: alphaSuSp.H:1
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
static const NamedEnum< thermoMode, 2 > thermoModeTypeNames_
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Definition: fvMatrix.H:282
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
addToRunTimeSelectionTable(option, fixedTemperatureConstraint, dictionary)
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
rhoReactionThermo & thermo
Definition: createFields.H:28
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
Macros for easy insertion into run-time selection tables.
solidificationMeltingSource(const word &sourceName, const word &modelType, const dictionary &dict, const fvMesh &mesh)
Construct from explicit source name and mesh.
virtual volScalarField & he()=0
Enthalpy/Internal energy [J/kg].
const word dictName() const
Return the local dictionary name (final part of scoped name)
Definition: dictionary.H:123
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Mesh &, const dimensionSet &, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return a temporary field constructed from name, mesh, dimensionSet.
dynamicFvMesh & mesh
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
A class for handling words, derived from string.
Definition: word.H:59
labelList fv(nPoints)
const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0)
Definition: dimensionSets.H:52
const Type & value() const
Return const reference to value.
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
Definition: doubleScalar.H:68
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Field< Type > & source()
Definition: fvMatrix.H:292
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
virtual tmp< volScalarField > Cp() const =0
Heat capacity at constant pressure [J/kg/K].
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const dimensionSet dimEnergy
virtual void addSup(fvMatrix< scalar > &eqn, const label fieldi)
Add explicit contribution to enthalpy equation.
dimensionedScalar pow3(const dimensionedScalar &ds)
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
scalar Cp(const scalar p, const scalar T) const
Definition: EtoHthermo.H:2
A special matrix type and solver, designed for finite volume solutions of scalar equations.
messageStream Info
Cell-set options abstract base class. Provides a base set of controls, e.g.:
Definition: cellSetOption.H:69
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
A class for managing temporary objects.
Definition: PtrList.H:53
const dimensionedVector & g
scalarField & diag()
Definition: lduMatrix.C:186
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.
zeroField Sp
Definition: alphaSuSp.H:2