Skip to content

Commit 8f6be26

Browse files
ChrisPsenicaChris Psenica
ChrisPsenica
and
Chris Psenica
authored
Added two flags to the heat flux Obj Func: "byUnitArea" & "formulation" (#771)
* Adding flags * Removing .VScode * Adding in DAFunction test for the total and daCustom HFX obj func. * removing VScode * Correcting issue in scalar for HFX_custom test function * Changed the byUnitArea scheme to take a bool as the input rather than string. I expanded the daCustom formulation to all regimes (incompressible, compressible, and solid domains). Lastly, added compressible and incompressible unit tests for HFX using both new flags. * Delete .vscode directory * Fixing formatting --------- Co-authored-by: Chris Psenica <[email protected]>
1 parent b8b1b97 commit 8f6be26

File tree

3 files changed

+168
-25
lines changed

3 files changed

+168
-25
lines changed

‎src/adjoint/DAFunction/DAFunctionWallHeatFlux.C

100755100644
+140-23
Original file line numberDiff line numberDiff line change
@@ -23,11 +23,11 @@ DAFunctionWallHeatFlux::DAFunctionWallHeatFlux(
2323
const DAIndex& daIndex,
2424
const word functionName)
2525
: DAFunction(
26-
mesh,
27-
daOption,
28-
daModel,
29-
daIndex,
30-
functionName),
26+
mesh,
27+
daOption,
28+
daModel,
29+
daIndex,
30+
functionName),
3131
wallHeatFlux_(
3232
IOobject(
3333
"wallHeatFlux",
@@ -40,10 +40,15 @@ DAFunctionWallHeatFlux::DAFunctionWallHeatFlux(
4040
"calculated")
4141
{
4242

43+
// check and assign values for scheme and formulation
44+
formMode_ = functionDict_.lookupOrDefault<word>("formulation", "default");
45+
calcMode_ = functionDict_.lookupOrDefault<bool>("byUnitArea", true);
46+
4347
if (mesh_.thisDb().foundObject<DATurbulenceModel>("DATurbulenceModel"))
4448
{
4549
DATurbulenceModel& daTurbModel =
4650
const_cast<DATurbulenceModel&>(daModel_.getDATurbulenceModel());
51+
4752
if (daTurbModel.getTurbModelType() == "incompressible")
4853
{
4954
// initialize the Prandtl number from transportProperties
@@ -89,7 +94,7 @@ DAFunctionWallHeatFlux::DAFunctionWallHeatFlux(
8994
}
9095
}
9196

92-
/// calculate the value of objective function
97+
//---------- Calculate Objective Function Value ----------
9398
scalar DAFunctionWallHeatFlux::calcFunction()
9499
{
95100
/*
@@ -101,85 +106,197 @@ scalar DAFunctionWallHeatFlux::calcFunction()
101106
functionValue: the sum of objective, reduced across all processors and scaled by "scale"
102107
*/
103108

104-
// always calculate the area of all the heat flux patches
105-
areaSum_ = 0.0;
106-
forAll(faceSources_, idxI)
109+
// only calculate the area of all the heat flux patches if scheme is byUnitArea
110+
if (calcMode_)
107111
{
108-
const label& functionFaceI = faceSources_[idxI];
109-
label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
110-
const label patchI = daIndex_.bFacePatchI[bFaceI];
111-
const label faceI = daIndex_.bFaceFaceI[bFaceI];
112-
areaSum_ += mesh_.magSf().boundaryField()[patchI][faceI];
112+
areaSum_ = 0.0;
113+
forAll(faceSources_, idxI)
114+
{
115+
const label& functionFaceI = faceSources_[idxI];
116+
label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
117+
const label patchI = daIndex_.bFacePatchI[bFaceI];
118+
const label faceI = daIndex_.bFaceFaceI[bFaceI];
119+
areaSum_ += mesh_.magSf().boundaryField()[patchI][faceI];
120+
}
121+
reduce(areaSum_, sumOp<scalar>());
113122
}
114-
reduce(areaSum_, sumOp<scalar>());
115123

116124
// initialize objFunValue
117125
scalar functionValue = 0.0;
118126

119127
volScalarField::Boundary& wallHeatFluxBf = wallHeatFlux_.boundaryFieldRef();
120128

129+
// calculate HFX for fluid domain
121130
if (mesh_.thisDb().foundObject<DATurbulenceModel>("DATurbulenceModel"))
122131
{
123132
DATurbulenceModel& daTurbModel =
124133
const_cast<DATurbulenceModel&>(daModel_.getDATurbulenceModel());
134+
135+
// calculate HFX for incompressible flow (no he for incompressible -> HFX = Cp * alphaEff * dT/dz)
125136
if (daTurbModel.getTurbModelType() == "incompressible")
126137
{
127-
// incompressible flow does not have he, so we do H = Cp * alphaEff * dT/dz
128138
const objectRegistry& db = mesh_.thisDb();
129139
const volScalarField& T = db.lookupObject<volScalarField>("T");
130140
volScalarField alphaEff = daTurbModel.alphaEff();
131141
const volScalarField::Boundary& TBf = T.boundaryField();
132142
const volScalarField::Boundary& alphaEffBf = alphaEff.boundaryField();
143+
133144
forAll(wallHeatFluxBf, patchI)
134145
{
135146
if (!wallHeatFluxBf[patchI].coupled())
136147
{
137-
wallHeatFluxBf[patchI] = Cp_ * alphaEffBf[patchI] * TBf[patchI].snGrad();
148+
// use OpenFOAM's snGrad()
149+
if (formMode_ == "default")
150+
{
151+
wallHeatFluxBf[patchI] = Cp_ * alphaEffBf[patchI] * TBf[patchI].snGrad();
152+
}
153+
// use DAFOAM's custom formulation
154+
else if (formMode_ == "daCustom")
155+
{
156+
forAll(wallHeatFluxBf[patchI], faceI)
157+
{
158+
scalar T2 = TBf[patchI][faceI];
159+
label nearWallCellIndex = mesh_.boundaryMesh()[patchI].faceCells()[faceI];
160+
scalar T1 = T[nearWallCellIndex];
161+
vector c1 = mesh_.Cf().boundaryField()[patchI][faceI];
162+
vector c2 = mesh_.C()[nearWallCellIndex];
163+
scalar d = mag(c1 - c2);
164+
scalar dTdz = (T2 - T1) / d;
165+
wallHeatFluxBf[patchI][faceI] = Cp_ * alphaEffBf[patchI][faceI] * dTdz;
166+
}
167+
}
168+
// error message incase of invalid entry
169+
else
170+
{
171+
FatalErrorIn(" ") << "formulation: "
172+
<< formMode_ << " not supported!"
173+
<< " Options are: default and daCustom."
174+
<< abort(FatalError);
175+
}
138176
}
139177
}
140178
}
179+
// calculate HFX for compressible flow (HFX = alphaEff * dHe/dz)
141180
else if (daTurbModel.getTurbModelType() == "compressible")
142181
{
143-
// compressible flow, H = alphaEff * dHE/dz
144182
fluidThermo& thermo = const_cast<fluidThermo&>(
145183
mesh_.thisDb().lookupObject<fluidThermo>("thermophysicalProperties"));
146184
volScalarField& he = thermo.he();
147185
const volScalarField::Boundary& heBf = he.boundaryField();
148186
volScalarField alphaEff = daTurbModel.alphaEff();
149187
const volScalarField::Boundary& alphaEffBf = alphaEff.boundaryField();
188+
150189
forAll(wallHeatFluxBf, patchI)
151190
{
152191
if (!wallHeatFluxBf[patchI].coupled())
153192
{
154-
wallHeatFluxBf[patchI] = alphaEffBf[patchI] * heBf[patchI].snGrad();
193+
// use OpenFOAM's snGrad()
194+
if (formMode_ == "default")
195+
{
196+
wallHeatFluxBf[patchI] = alphaEffBf[patchI] * heBf[patchI].snGrad();
197+
}
198+
// use DAFOAM's custom formulation
199+
else if (formMode_ == "daCustom")
200+
{
201+
forAll(wallHeatFluxBf[patchI], faceI)
202+
{
203+
scalar He2 = heBf[patchI][faceI];
204+
label nearWallCellIndex = mesh_.boundaryMesh()[patchI].faceCells()[faceI];
205+
scalar He1 = he[nearWallCellIndex];
206+
vector c1 = mesh_.Cf().boundaryField()[patchI][faceI];
207+
vector c2 = mesh_.C()[nearWallCellIndex];
208+
scalar d = mag(c1 - c2);
209+
scalar dHedz = (He2 - He1) / d;
210+
wallHeatFluxBf[patchI][faceI] = alphaEffBf[patchI][faceI] * dHedz;
211+
}
212+
}
213+
// error message incase of invalid entry
214+
else
215+
{
216+
FatalErrorIn(" ") << "formulation: "
217+
<< formMode_ << " not supported!"
218+
<< " Options are: default and daCustom."
219+
<< abort(FatalError);
220+
}
155221
}
156222
}
157223
}
158224
}
225+
226+
// calculate HFX for solid domain (HFX = k * dT/dz, where k = DT / rho / Cp)
159227
else
160228
{
161-
// solid. H = k * dT/dz, where k = DT / rho / Cp
162229
const objectRegistry& db = mesh_.thisDb();
163230
const volScalarField& T = db.lookupObject<volScalarField>("T");
164231
const volScalarField::Boundary& TBf = T.boundaryField();
232+
165233
forAll(wallHeatFluxBf, patchI)
166234
{
167235
if (!wallHeatFluxBf[patchI].coupled())
168236
{
169-
wallHeatFluxBf[patchI] = k_ * TBf[patchI].snGrad();
237+
238+
// use OpenFOAM's snGrad()
239+
if (formMode_ == "default")
240+
{
241+
wallHeatFluxBf[patchI] = k_ * TBf[patchI].snGrad();
242+
}
243+
// use DAFOAM's custom formulation
244+
else if (formMode_ == "daCustom")
245+
{
246+
forAll(wallHeatFluxBf[patchI], faceI)
247+
{
248+
scalar T2 = TBf[patchI][faceI];
249+
label nearWallCellIndex = mesh_.boundaryMesh()[patchI].faceCells()[faceI];
250+
scalar T1 = T[nearWallCellIndex];
251+
vector c1 = mesh_.Cf().boundaryField()[patchI][faceI];
252+
vector c2 = mesh_.C()[nearWallCellIndex];
253+
scalar d = mag(c1 - c2);
254+
scalar dTdz = (T2 - T1) / d;
255+
wallHeatFluxBf[patchI][faceI] = k_ * dTdz;
256+
}
257+
}
258+
// error message incase of invalid entry
259+
else
260+
{
261+
FatalErrorIn(" ") << "formulation: "
262+
<< formMode_ << " not supported!"
263+
<< " Options are: default and daCustom."
264+
<< abort(FatalError);
265+
}
170266
}
171267
}
172268
}
269+
173270
// calculate area weighted heat flux
271+
scalar val = 0;
174272
forAll(faceSources_, idxI)
175273
{
176274
const label& functionFaceI = faceSources_[idxI];
177275
label bFaceI = functionFaceI - daIndex_.nLocalInternalFaces;
178276
const label patchI = daIndex_.bFacePatchI[bFaceI];
179277
const label faceI = daIndex_.bFaceFaceI[bFaceI];
180-
181278
scalar area = mesh_.magSf().boundaryField()[patchI][faceI];
182-
scalar val = scale_ * wallHeatFluxBf[patchI][faceI] * area / areaSum_;
279+
280+
// represent wallHeatFlux by unit area
281+
if (calcMode_)
282+
{
283+
val = scale_ * wallHeatFluxBf[patchI][faceI] * area / areaSum_;
284+
}
285+
// represent wallHeatFlux as total heat transfer through surface
286+
else if (!calcMode_)
287+
{
288+
val = scale_ * wallHeatFluxBf[patchI][faceI] * area;
289+
}
290+
// error message incase of invalid entry
291+
else
292+
{
293+
FatalErrorIn(" ") << "byUnitArea: "
294+
<< calcMode_ << " not supported!"
295+
<< " Options are: True (default) and False."
296+
<< abort(FatalError);
297+
}
298+
299+
// update obj. func. val
183300
functionValue += val;
184301
}
185302

‎src/adjoint/DAFunction/DAFunctionWallHeatFlux.H

100755100644
+4-1
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,10 @@ protected:
4242
scalar areaSum_ = -9999.0;
4343

4444
/// if calculating flux per unit area or total, which mode to use
45-
word calcMode_;
45+
bool calcMode_;
46+
47+
/// if calculating wallHeatFlux by OpenFOAMs snGrad() or DAFOAM's custom (daCustom) formulation
48+
word formMode_;
4649

4750
public:
4851
TypeName("wallHeatFlux");

‎tests/runUnitTests_DAFunction.py

+24-1
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,14 @@
6565
"patches": ["walls"],
6666
"scale": 0.001,
6767
},
68+
"HFX_custom": {
69+
"type": "wallHeatFlux",
70+
"byUnitArea": False,
71+
"formulation": "daCustom",
72+
"source": "patchToFace",
73+
"patches": ["walls"],
74+
"scale": 1.0,
75+
},
6876
"PMean": {
6977
"type": "patchMean",
7078
"source": "patchToFace",
@@ -175,6 +183,7 @@
175183
"CMZ": 7.1768354637131795,
176184
"TP1": 2.5561992647012914,
177185
"HFX": 1.5730364723151125,
186+
"HFX_custom": 5502.596904327625,
178187
"PMean": 77.80996323506456,
179188
"UMean": 15.469850053816028,
180189
"skewness": 1.3140396235456926,
@@ -234,6 +243,14 @@
234243
"patches": ["walls"],
235244
"scale": 1.0,
236245
},
246+
"HFX_custom": {
247+
"type": "wallHeatFlux",
248+
"byUnitArea": False,
249+
"formulation": "daCustom",
250+
"source": "patchToFace",
251+
"patches": ["walls"],
252+
"scale": 1.0,
253+
},
237254
"TTR": {
238255
"type": "totalTemperatureRatio",
239256
"source": "patchToFace",
@@ -268,7 +285,13 @@
268285
if gcomm.rank == 0:
269286
print(funcs)
270287

271-
funcs_ref = {"HFX": 15.494946599784122, "TTR": 1.000000415241806, "MFR": 128.22449523986853, "TPR": 0.9934250594975235}
288+
funcs_ref = {
289+
"HFX": 15.494946599784122,
290+
"HFX_custom": 54.15546957805019,
291+
"TTR": 1.000000415241806,
292+
"MFR": 128.22449523986853,
293+
"TPR": 0.9934250594975235,
294+
}
272295

273296
fail = 0
274297
failedFunc = []

0 commit comments

Comments
 (0)