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-2021 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  fvModel,
57  solidificationMeltingSource,
59  );
60  }
61 }
62 
65 
66 
67 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
68 
69 void Foam::fv::solidificationMeltingSource::readCoeffs()
70 {
71  Tsol_ = coeffs().lookup<scalar>("Tsol");
72  Tliq_ = coeffs().lookupOrDefault<scalar>("Tliq", Tsol_);
73  alpha1e_ = coeffs().lookupOrDefault<scalar>("alpha1e", 0.0);
74  L_ = coeffs().lookup<scalar>("L");
75 
76  relax_ = coeffs().lookupOrDefault("relax", 0.9);
77 
78  mode_ = thermoModeTypeNames_.read(coeffs().lookup("thermoMode"));
79 
80  rhoRef_ = coeffs().lookup<scalar>("rhoRef");
81  TName_ = coeffs().lookupOrDefault<word>("T", "T");
82  CpName_ = coeffs().lookupOrDefault<word>("Cp", "Cp");
83  UName_ = coeffs().lookupOrDefault<word>("U", "U");
84  phiName_ = coeffs().lookupOrDefault<word>("phi", "phi");
85 
86  Cu_ = coeffs().lookupOrDefault<scalar>("Cu", 100000);
87  q_ = coeffs().lookupOrDefault("q", 0.001);
88 
89  beta_ = coeffs().lookup<scalar>("beta");
90 }
91 
92 
94 Foam::fv::solidificationMeltingSource::Cp() const
95 {
96  switch (mode_)
97  {
98  case thermoMode::thermo:
99  {
100  const basicThermo& thermo =
101  mesh().lookupObject<basicThermo>(basicThermo::dictName);
102 
103  return thermo.Cp();
104  break;
105  }
106  case thermoMode::lookup:
107  {
108  if (CpName_ == "CpRef")
109  {
110  scalar CpRef = coeffs().lookup<scalar>("CpRef");
111 
112  return volScalarField::New
113  (
114  name() + ":Cp",
115  mesh(),
117  (
119  CpRef
120  ),
121  extrapolatedCalculatedFvPatchScalarField::typeName
122  );
123  }
124  else
125  {
126  return mesh().lookupObject<volScalarField>(CpName_);
127  }
128 
129  break;
130  }
131  default:
132  {
134  << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
135  << abort(FatalError);
136  }
137  }
138 
139  return tmp<volScalarField>(nullptr);
140 }
141 
142 
143 Foam::vector Foam::fv::solidificationMeltingSource::g() const
144 {
145  if (mesh().foundObject<uniformDimensionedVectorField>("g"))
146  {
147  const uniformDimensionedVectorField& value =
148  mesh().lookupObject<uniformDimensionedVectorField>("g");
149  return value.value();
150  }
151  else
152  {
153  return coeffs().lookup("g");
154  }
155 }
156 
157 
158 void Foam::fv::solidificationMeltingSource::update
159 (
160  const volScalarField& Cp
161 ) const
162 {
163  if (curTimeIndex_ == mesh().time().timeIndex())
164  {
165  return;
166  }
167 
168  if (debug)
169  {
170  Info<< type() << ": " << name()
171  << " - updating phase indicator" << endl;
172  }
173 
174  // update old time alpha1 field
175  alpha1_.oldTime();
176 
177  const volScalarField& T = mesh().lookupObject<volScalarField>(TName_);
178 
179  const labelList& cells = set_.cells();
180 
181  forAll(cells, i)
182  {
183  const label celli = cells[i];
184 
185  const scalar Tc = T[celli];
186  const scalar Cpc = Cp[celli];
187  const scalar alpha1New =
188  alpha1_[celli]
189  + relax_*Cpc
190  *(
191  Tc
192  - max
193  (
194  Tsol_,
195  Tsol_
196  + (Tliq_ - Tsol_)*(alpha1_[celli] - alpha1e_)/(1 - alpha1e_)
197  )
198  )/L_;
199 
200  alpha1_[celli] = max(0, min(alpha1New, 1));
201  deltaT_[i] =
202  Tc
203  - max
204  (
205  Tsol_,
206  Tsol_
207  + (Tliq_ - Tsol_)*(alpha1_[celli] - alpha1e_)/(1 - alpha1e_)
208  );
209  }
210 
211  alpha1_.correctBoundaryConditions();
212 
213  curTimeIndex_ = mesh().time().timeIndex();
214 }
215 
216 
217 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
218 
220 (
221  const word& name,
222  const word& modelType,
223  const dictionary& dict,
224  const fvMesh& mesh
225 )
226 :
227  fvModel(name, modelType, dict, mesh),
228  set_(coeffs(), mesh),
229  Tsol_(NaN),
230  Tliq_(NaN),
231  alpha1e_(NaN),
232  L_(NaN),
233  relax_(NaN),
234  mode_(thermoMode::thermo),
235  rhoRef_(NaN),
236  TName_(word::null),
237  CpName_(word::null),
238  UName_(word::null),
239  phiName_(word::null),
240  Cu_(NaN),
241  q_(NaN),
242  beta_(NaN),
243  alpha1_
244  (
245  IOobject
246  (
247  this->name() + ":alpha1",
248  mesh.time().timeName(),
249  mesh,
252  ),
253  mesh,
255  zeroGradientFvPatchScalarField::typeName
256  ),
257  curTimeIndex_(-1),
258  deltaT_(set_.cells().size(), 0)
259 {
260  readCoeffs();
261 }
262 
263 
264 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
265 
267 {
268  switch (mode_)
269  {
270  case thermoMode::thermo:
271  {
272  const basicThermo& thermo =
273  mesh().lookupObject<basicThermo>(basicThermo::dictName);
274 
275  return wordList({UName_, thermo.he().name()});
276  }
277  case thermoMode::lookup:
278  {
279  return wordList({UName_, TName_});
280  }
281  }
282 
283  return wordList::null();
284 }
285 
286 
288 (
289  fvMatrix<scalar>& eqn,
290  const word& fieldName
291 ) const
292 {
293  apply(geometricOneField(), eqn);
294 }
295 
296 
298 (
299  const volScalarField& rho,
300  fvMatrix<scalar>& eqn,
301  const word& fieldName
302 ) const
303 {
304  apply(rho, eqn);
305 }
306 
307 
309 (
310  fvMatrix<vector>& eqn,
311  const word& fieldName
312 ) const
313 {
314  if (debug)
315  {
316  Info<< type() << ": applying source to " << eqn.psi().name() << endl;
317  }
318 
319  const volScalarField Cp(this->Cp());
320 
321  update(Cp);
322 
323  vector g = this->g();
324 
325  scalarField& Sp = eqn.diag();
326  vectorField& Su = eqn.source();
327  const scalarField& V = mesh().V();
328 
329  const labelList& cells = set_.cells();
330 
331  forAll(cells, i)
332  {
333  const label celli = cells[i];
334 
335  const scalar Vc = V[celli];
336  const scalar alpha1c = alpha1_[celli];
337 
338  const scalar S = -Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_);
339  const vector Sb = rhoRef_*g*beta_*deltaT_[i];
340 
341  Sp[celli] += Vc*S;
342  Su[celli] += Vc*Sb;
343  }
344 }
345 
346 
348 (
349  const volScalarField& rho,
350  fvMatrix<vector>& eqn,
351  const word& fieldName
352 ) const
353 {
354  // Momentum source uses a Boussinesq approximation - redirect
355  addSup(eqn, fieldName);
356 }
357 
358 
360 {
361  set_.updateMesh(mpm);
362 }
363 
364 
366 {
367  if (fvModel::read(dict))
368  {
369  set_.read(coeffs());
370  readCoeffs();
371  return true;
372  }
373  else
374  {
375  return false;
376  }
377 
378  return false;
379 }
380 
381 
382 // ************************************************************************* //
defineTypeNameAndDebug(fixedTemperatureConstraint, 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:303
fluidReactionThermo & thermo
Definition: createFields.H:28
Base-class for fluid and solid thermodynamic properties.
Definition: basicThermo.H:77
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
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:323
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Definition: fvMatrix.H:282
dimensionedSymmTensor sqr(const dimensionedVector &dv)
virtual bool read(const dictionary &dict)
Read source dictionary.
Definition: fvModel.C:158
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Finite volume model abstract base class.
Definition: fvModel.H:55
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:636
const dimensionSet dimless
solidificationMeltingSource(const word &name, const word &modelType, const dictionary &dict, const fvMesh &mesh)
Construct from explicit source name and mesh.
static const List< word > & null()
Return a null List.
Definition: ListI.H:118
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
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
Macros for easy insertion into run-time selection tables.
virtual volScalarField & he()=0
Enthalpy/Internal energy [J/kg].
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
const cellShapeList & cells
A class for handling words, derived from string.
Definition: word.H:59
labelList fv(nPoints)
tmp< fvMatrix< Type > > S(const Pair< tmp< volScalarField::Internal >> &, const GeometricField< Type, fvPatchField, volMesh > &)
const Type & value() const
Return const reference to value.
static const word null
An empty word.
Definition: word.H:77
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
errorManip< error > abort(error &err)
Definition: errorManip.H:131
virtual bool read(const dictionary &dict)
Read source dictionary.
static const word dictName
Name of the thermophysical properties dictionary.
Definition: basicThermo.H:129
Field< Type > & source()
Definition: fvMatrix.H:292
const dimensionSet dimEnergy
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 dimMass
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const tmp< volScalarField::Internal > & Sp
Definition: alphaSuSp.H:7
addToRunTimeSelectionTable(fvConstraint, fixedTemperatureConstraint, dictionary)
List< word > wordList
A List of words.
Definition: fileName.H:54
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
label timeIndex
Definition: getTimeIndex.H:4
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
virtual void updateMesh(const mapPolyMesh &)
Update for mesh changes.
A special matrix type and solver, designed for finite volume solutions of scalar equations.
messageStream Info
virtual void addSup(fvMatrix< scalar > &eqn, const word &fieldName) const
Add explicit contribution to enthalpy equation.
const tmp< volScalarField::Internal > & Su
Definition: alphaSuSp.H:6
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
const dimensionSet dimTemperature
Namespace for OpenFOAM.
virtual wordList addSupFields() const
Return the list of fields for which the fvModel adds source term.