DAInputThermalCoupling.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAInputThermalCoupling.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAInputThermalCoupling, 0);
16 addToRunTimeSelectionTable(DAInput, DAInputThermalCoupling, dictionary);
17 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
18 
19 DAInputThermalCoupling::DAInputThermalCoupling(
20  const word inputName,
21  const word inputType,
22  fvMesh& mesh,
23  const DAOption& daOption,
24  const DAModel& daModel,
25  const DAIndex& daIndex)
26  : DAInput(
27  inputName,
28  inputType,
29  mesh,
30  daOption,
31  daModel,
32  daIndex)
33 {
34 
35  daOption_.getAllOptions().subDict("inputInfo").subDict(inputName_).readEntry("patches", patches_);
36  // NOTE: always sort the patch because the order of the patch element matters in CHT coupling
37  sort(patches_);
38 
39  // check discipline
40  discipline_ = daOption_.getAllOptions().getWord("discipline");
41 
42  // check coupling mode and validate
43  distanceMode_ = daOption_.getAllOptions().getWord("wallDistanceMethod");
44  if (distanceMode_ != "daCustom" && distanceMode_ != "default")
45  {
46  FatalErrorIn(" ") << "wallDistanceMethod: "
47  << distanceMode_ << " not supported!"
48  << " Options are: default and daCustom."
49  << abort(FatalError);
50  }
51 
52  size_ = 0;
53  forAll(patches_, idxI)
54  {
55  word patchName = patches_[idxI];
56  label patchI = mesh_.boundaryMesh().findPatchID(patchName);
57  forAll(mesh_.boundaryMesh()[patchI], faceI)
58  {
59  size_++;
60  }
61  }
62  // we have two sets of variable to transfer
63  size_ *= 2;
64 }
65 
66 void DAInputThermalCoupling::run(const scalarList& input)
67 {
68  /*
69  Description:
70  Assign the input to OF fields
71  */
72 
73  volScalarField& T =
74  const_cast<volScalarField&>(mesh_.thisDb().lookupObject<volScalarField>("T"));
75 
76  // ********* first loop, set the refValue
77  label counterI = 0;
78  forAll(patches_, idxI)
79  {
80  // get the patch id label
81  word patchName = patches_[idxI];
82  label patchI = mesh_.boundaryMesh().findPatchID(patchName);
83 
84  mixedFvPatchField<scalar>& mixedPatch =
85  refCast<mixedFvPatchField<scalar>>(T.boundaryFieldRef()[patchI]);
86 
87  forAll(mixedPatch.refValue(), faceI)
88  {
89  mixedPatch.refValue()[faceI] = input[counterI];
90  mixedPatch.refGrad()[faceI] = 0;
91  counterI++;
92  }
93  }
94 
95  // ********* second loop, set the valueFraction:
96  // neighKDeltaCoeffs / ( neighKDeltaCoeffs + myKDeltaCoeffs)
97  scalar deltaCoeffs = 0;
98 
99  if (discipline_ == "aero")
100  {
101 
102  // for incompressible flow Q = Cp * alphaEff * dT/dz, so kappa = Cp * alphaEff
103  DATurbulenceModel& daTurb = const_cast<DATurbulenceModel&>(daModel_.getDATurbulenceModel());
104  word turbModelType = daTurb.getTurbModelType();
105 
106  if (turbModelType == "incompressible")
107  {
108  volScalarField alphaEff = daTurb.alphaEff();
109 
110  IOdictionary transportProperties(
111  IOobject(
112  "transportProperties",
113  mesh_.time().constant(),
114  mesh_,
115  IOobject::MUST_READ,
116  IOobject::NO_WRITE,
117  false));
118  scalar Cp = readScalar(transportProperties.lookup("Cp"));
119 
120  forAll(patches_, idxI)
121  {
122  // get the patch id label
123  word patchName = patches_[idxI];
124  label patchI = mesh_.boundaryMesh().findPatchID(patchName);
125 
126  mixedFvPatchField<scalar>& mixedPatch =
127  refCast<mixedFvPatchField<scalar>>(T.boundaryFieldRef()[patchI]);
128 
129  forAll(mesh_.boundaryMesh()[patchI], faceI)
130  {
131  if (distanceMode_ == "default")
132  {
133  // deltaCoeffs = 1 / d
134  deltaCoeffs = T.boundaryField()[patchI].patch().deltaCoeffs()[faceI];
135  }
136  else if (distanceMode_ == "daCustom")
137  {
138  label nearWallCellIndex = mesh_.boundaryMesh()[patchI].faceCells()[faceI];
139  vector c1 = mesh_.Cf().boundaryField()[patchI][faceI];
140  vector c2 = mesh_.C()[nearWallCellIndex];
141  scalar d = mag(c1 - c2);
142  deltaCoeffs = 1 / d;
143  }
144  scalar alphaEffBf = alphaEff.boundaryField()[patchI][faceI];
145  scalar myKDeltaCoeffs = Cp * alphaEffBf * deltaCoeffs;
146  // NOTE: we continue to use the counterI from the first loop
147  scalar neighKDeltaCoeffs = input[counterI];
148  mixedPatch.valueFraction()[faceI] = neighKDeltaCoeffs / (myKDeltaCoeffs + neighKDeltaCoeffs);
149  counterI++;
150  }
151  }
152  }
153 
154  if (turbModelType == "compressible")
155  {
156  // for compressible flow Q = alphaEff * dHE/dz, so if enthalpy is used, kappa = Cp * alphaEff
157  // if the internalEnergy is used, kappa = (Cp - R) * alphaEff
158  volScalarField alphaEff = daTurb.alphaEff();
159  // compressible flow, H = alphaEff * dHE/dz
160  const fluidThermo& thermo = mesh_.thisDb().lookupObject<fluidThermo>("thermophysicalProperties");
161  const volScalarField& he = thermo.he();
162 
163  const IOdictionary& thermoDict = mesh_.thisDb().lookupObject<IOdictionary>("thermophysicalProperties");
164  dictionary mixSubDict = thermoDict.subDict("mixture");
165  dictionary specieSubDict = mixSubDict.subDict("specie");
166  scalar molWeight = specieSubDict.getScalar("molWeight");
167  dictionary thermodynamicsSubDict = mixSubDict.subDict("thermodynamics");
168  scalar Cp = thermodynamicsSubDict.getScalar("Cp");
169 
170  // 8314.4700665 gas constant in OpenFOAM
171  // src/OpenFOAM/global/constants/thermodynamic/thermodynamicConstants.H
172  scalar RR = Foam::constant::thermodynamic::RR;
173 
174  // R = RR/molWeight
175  // Foam::specie::R() function in src/thermophysicalModels/specie/specie/specieI.H
176  scalar R = RR / molWeight;
177 
178  scalar tmpVal = 0;
179  // e = (Cp - R) * T, so Q = alphaEff * (Cp-R) * dT/dz
180  if (he.name() == "e")
181  {
182  tmpVal = Cp - R;
183  }
184  // h = Cp * T, so Q = alphaEff * Cp * dT/dz
185  else
186  {
187  tmpVal = Cp;
188  }
189 
190  forAll(patches_, idxI)
191  {
192  // get the patch id label
193  word patchName = patches_[idxI];
194  label patchI = mesh_.boundaryMesh().findPatchID(patchName);
195 
196  mixedFvPatchField<scalar>& mixedPatch =
197  refCast<mixedFvPatchField<scalar>>(T.boundaryFieldRef()[patchI]);
198 
199  forAll(mesh_.boundaryMesh()[patchI], faceI)
200  {
201  if (distanceMode_ == "default")
202  {
203  // deltaCoeffs = 1 / d
204  deltaCoeffs = T.boundaryField()[patchI].patch().deltaCoeffs()[faceI];
205  }
206  else if (distanceMode_ == "daCustom")
207  {
208  label nearWallCellIndex = mesh_.boundaryMesh()[patchI].faceCells()[faceI];
209  vector c1 = mesh_.Cf().boundaryField()[patchI][faceI];
210  vector c2 = mesh_.C()[nearWallCellIndex];
211  scalar d = mag(c1 - c2);
212  deltaCoeffs = 1 / d;
213  }
214  scalar alphaEffBf = alphaEff.boundaryField()[patchI][faceI];
215  scalar myKDeltaCoeffs = tmpVal * alphaEffBf * deltaCoeffs;
216  // NOTE: we continue to use the counterI from the first loop
217  scalar neighKDeltaCoeffs = input[counterI];
218  mixedPatch.valueFraction()[faceI] = neighKDeltaCoeffs / (myKDeltaCoeffs + neighKDeltaCoeffs);
219  counterI++;
220  }
221  }
222  }
223  }
224  else if (discipline_ == "thermal")
225  {
226  // for solid solvers Q = k * dT/dz, so kappa = k
227  IOdictionary solidProperties(
228  IOobject(
229  "solidProperties",
230  mesh_.time().constant(),
231  mesh_,
232  IOobject::MUST_READ,
233  IOobject::NO_WRITE,
234  false));
235  scalar k = readScalar(solidProperties.lookup("k"));
236 
237  forAll(patches_, idxI)
238  {
239  // get the patch id label
240  word patchName = patches_[idxI];
241  label patchI = mesh_.boundaryMesh().findPatchID(patchName);
242 
243  forAll(mesh_.boundaryMesh()[patchI], faceI)
244  {
245  if (distanceMode_ == "default")
246  {
247  // deltaCoeffs = 1 / d
248  deltaCoeffs = T.boundaryField()[patchI].patch().deltaCoeffs()[faceI];
249  }
250  else if (distanceMode_ == "daCustom")
251  {
252  label nearWallCellIndex = mesh_.boundaryMesh()[patchI].faceCells()[faceI];
253  vector c1 = mesh_.Cf().boundaryField()[patchI][faceI];
254  vector c2 = mesh_.C()[nearWallCellIndex];
255  scalar d = mag(c1 - c2);
256  deltaCoeffs = 1 / d;
257  }
258  mixedFvPatchField<scalar>& mixedPatch = refCast<mixedFvPatchField<scalar>>(T.boundaryFieldRef()[patchI]);
259  scalar myKDeltaCoeffs = k * deltaCoeffs;
260  // NOTE: we continue to use the counterI from the first loop
261  scalar neighKDeltaCoeffs = input[counterI];
262  mixedPatch.valueFraction()[faceI] = neighKDeltaCoeffs / (myKDeltaCoeffs + neighKDeltaCoeffs);
263  counterI++;
264  }
265  }
266  }
267  else
268  {
269  FatalErrorIn("DAInputThermalCoupling::run") << " discipline not valid! "
270  << abort(FatalError);
271  }
272 
273  T.correctBoundaryConditions();
274 }
275 
276 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
277 
278 } // End namespace Foam
279 
280 // ************************************************************************* //
Foam::DAInputThermalCoupling::size_
label size_
the total face number for all the patches_ owned by this processor
Definition: DAInputThermalCoupling.H:53
Foam::DAInputThermalCoupling::run
virtual void run(const scalarList &input)
Definition: DAInputThermalCoupling.C:66
Foam::DAInput::inputName_
const word inputName_
name of the input
Definition: DAInput.H:42
Foam::DAOption
Definition: DAOption.H:29
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
Foam::DAInput::daModel_
const DAModel & daModel_
DAIndex object.
Definition: DAInput.H:54
he
volScalarField & he
Definition: EEqnRhoPimple.H:5
Foam::DAInputThermalCoupling::distanceMode_
word distanceMode_
if calculating wallHeatFlux by OpenFOAMs snGrad() or DAFOAM's custom (daCustom) formulation
Definition: DAInputThermalCoupling.H:59
thermo
fluidThermo & thermo
Definition: createRefsRhoPimple.H:6
Foam::DAOption::getAllOptions
const dictionary & getAllOptions() const
return a reference of allOptions_ dictionary
Definition: DAOption.H:56
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAIndex
Definition: DAIndex.H:32
Foam::DAInput::daOption_
const DAOption & daOption_
DAOption object.
Definition: DAInput.H:51
Foam::DAModel
Definition: DAModel.H:57
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
k
dimensionedScalar & k
Definition: createRefsHeatTransfer.H:6
Foam::DAInputThermalCoupling::patches_
wordList patches_
list of patch names for the thermal var
Definition: DAInputThermalCoupling.H:50
Foam
Definition: checkGeometry.C:32
T
volScalarField & T
Definition: createRefsHeatTransfer.H:5
Foam::DAInput
Definition: DAInput.H:30
alphaEff
volScalarField alphaEff("alphaEff", turbulencePtr_->nu()/Pr+alphat)
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DATurbulenceModel::getTurbModelType
word getTurbModelType() const
Definition: DATurbulenceModel.H:260
Foam::DATurbulenceModel
Definition: DATurbulenceModel.H:52
DAInputThermalCoupling.H
Foam::DATurbulenceModel::alphaEff
tmp< volScalarField > alphaEff()
return effective thermal diffusivity
Definition: DATurbulenceModel.C:233
Foam::DAInputThermalCoupling::DAInputThermalCoupling
DAInputThermalCoupling(const word inputName, const word inputType, fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex)
Definition: DAInputThermalCoupling.C:19
Foam::DAModel::getDATurbulenceModel
const DATurbulenceModel & getDATurbulenceModel() const
get a reference to DATurbulenceModel
Definition: DAModel.C:176
Foam::DAInput::mesh_
fvMesh & mesh_
fvMesh
Definition: DAInput.H:48
Foam::DAInputThermalCoupling::discipline_
word discipline_
whether this is flow or solid
Definition: DAInputThermalCoupling.H:56

Generated by   doxygen 1.8.17