objectiveForce.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2007-2019 PCOpt/NTUA
9  Copyright (C) 2013-2019 FOSS GP
10  Copyright (C) 2019-2020 OpenCFD Ltd.
11 -------------------------------------------------------------------------------
12 License
13  This file is part of OpenFOAM.
14 
15  OpenFOAM is free software: you can redistribute it and/or modify it
16  under the terms of the GNU General Public License as published by
17  the Free Software Foundation, either version 3 of the License, or
18  (at your option) any later version.
19 
20  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23  for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include "objectiveForce.H"
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37 
38 namespace objectives
39 {
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
45 (
49 );
50 
51 
52 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
53 
55 (
56  const fvMesh& mesh,
57  const dictionary& dict,
58  const word& adjointSolverName,
59  const word& primalSolverName
60 )
61 :
62  objectiveIncompressible(mesh, dict, adjointSolverName, primalSolverName),
63  forcePatches_
64  (
65  mesh_.boundaryMesh().patchSet
66  (
67  dict.get<wordRes>("patches")
68  ).sortedToc()
69  ),
70  forceDirection_(dict.get<vector>("direction")),
71  Aref_(dict.get<scalar>("Aref")),
72  rhoInf_(dict.get<scalar>("rhoInf")),
73  UInf_(dict.get<scalar>("UInf")),
74  stressXPtr_
75  (
76  Foam::createZeroFieldPtr<vector>
77  (
78  mesh_, "stressX", dimLength/sqr(dimTime)
79  )
80  ),
81  stressYPtr_
82  (
83  Foam::createZeroFieldPtr<vector>
84  (
85  mesh_, "stressY", dimLength/sqr(dimTime)
86  )
87  ),
88  stressZPtr_
89  (
90  Foam::createZeroFieldPtr<vector>
91  (
92  mesh_, "stressZ", dimLength/sqr(dimTime)
93  )
94  )
95 {
96  // Sanity check and print info
97  if (forcePatches_.empty())
98  {
100  << "No valid patch name on which to minimize " << type() << endl
101  << exit(FatalError);
102  }
103  if (debug)
104  {
105  Info<< "Minimizing " << type() << " in patches:" << endl;
106  for (const label patchI : forcePatches_)
107  {
108  Info<< "\t " << mesh_.boundary()[patchI].name() << endl;
109  }
110  }
111 
112  // Allocate boundary field pointers
113  bdJdpPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
114  bdSdbMultPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
115  bdxdbMultPtr_.reset(createZeroBoundaryPtr<vector>(mesh_));
116  bdJdStressPtr_.reset(createZeroBoundaryPtr<tensor>(mesh_));
117 }
118 
119 
120 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
121 
123 {
124  vector pressureForce(Zero);
125  vector viscousForce(Zero);
126  vector cumulativeForce(Zero);
127 
128 
129  const volScalarField& p = vars_.pInst();
132 
133  volSymmTensorField devReff(turbulence->devReff());
134 
135  for (const label patchI : forcePatches_)
136  {
137  pressureForce += gSum
138  (
139  mesh_.Sf().boundaryField()[patchI] * p.boundaryField()[patchI]
140  );
141  // Viscous term calculated using the full tensor derivative
142  viscousForce += gSum
143  (
144  devReff.boundaryField()[patchI]
145  & mesh_.Sf().boundaryField()[patchI]
146  );
147  }
148 
149  cumulativeForce = pressureForce + viscousForce;
150 
151  scalar force = cumulativeForce & forceDirection_;
152 
153  // Intentionally not using denom - derived might implement virtual denom()
154  // function differently
155  scalar Cforce = force/(0.5*UInf_*UInf_*Aref_);
156 
157  DebugInfo
158  << "Force|Coeff " << force << "|" << Cforce << endl;
159 
160  J_ = Cforce;
161 
162  return Cforce;
163 }
164 
165 
167 {
168  for (const label patchI : forcePatches_)
169  {
170  bdJdpPtr_()[patchI] = forceDirection_/denom();
171  }
172 }
173 
174 
176 {
177  // Compute contributions with mean fields, if present
178  const volScalarField& p = vars_.p();
179  const volVectorField& U = vars_.U();
181  turbVars = vars_.RASModelVariables();
182  const singlePhaseTransportModel& lamTransp = vars_.laminarTransport();
183 
184  tmp<volSymmTensorField> tdevReff = turbVars->devReff(lamTransp, U);
185  const volSymmTensorField& devReff = tdevReff();
186 
187  for (const label patchI : forcePatches_)
188  {
189  bdSdbMultPtr_()[patchI] =
190  (
191  (
192  forceDirection_& devReff.boundaryField()[patchI]
193  )
194  + (forceDirection_)*p.boundaryField()[patchI]
195  )
196  /denom();
197  }
198 }
199 
200 
202 {
203  const volScalarField& p = vars_.p();
204  const volVectorField& U = vars_.U();
205 
207  turbVars = vars_.RASModelVariables();
208  const singlePhaseTransportModel& lamTransp = vars_.laminarTransport();
209 
210  //tmp<volSymmTensorField> tdevReff = turbVars->devReff(lamTransp, U);
211  //const volSymmTensorField& devReff = tdevReff();
212 
213  volScalarField nuEff(lamTransp.nu() + turbVars->nutRef());
214  volTensorField gradU(fvc::grad(U));
215  volTensorField::Boundary& gradUbf = gradU.boundaryFieldRef();
216 
217  // Explicitly correct the boundary gradient to get rid of
218  // the tangential component
219  forAll(mesh_.boundary(), patchI)
220  {
221  const fvPatch& patch = mesh_.boundary()[patchI];
222  if (isA<wallFvPatch>(patch))
223  {
224  tmp<vectorField> nf = patch.nf();
225  gradUbf[patchI] = nf*U.boundaryField()[patchI].snGrad();
226  }
227  }
228 
229  volTensorField stress(nuEff*(gradU + T(gradU)));
230 
231  stressXPtr_().replace(0, stress.component(0));
232  stressXPtr_().replace(1, stress.component(1));
233  stressXPtr_().replace(2, stress.component(2));
234 
235  stressYPtr_().replace(0, stress.component(3));
236  stressYPtr_().replace(1, stress.component(4));
237  stressYPtr_().replace(2, stress.component(5));
238 
239  stressZPtr_().replace(0, stress.component(6));
240  stressZPtr_().replace(1, stress.component(7));
241  stressZPtr_().replace(2, stress.component(8));
242 
243  volTensorField gradStressX(fvc::grad(stressXPtr_()));
244  volTensorField gradStressY(fvc::grad(stressYPtr_()));
245  volTensorField gradStressZ(fvc::grad(stressZPtr_()));
246 
247  // the notorious second-order derivative at the wall. Use with caution!
248  volVectorField gradp(fvc::grad(p));
249 
250  for (const label patchI : forcePatches_)
251  {
252  const fvPatch& patch = mesh_.boundary()[patchI];
253  tmp<vectorField> tnf = patch.nf();
254  const vectorField& nf = tnf();
255  bdxdbMultPtr_()[patchI] =
256  (
257  (
258  (
259  -(forceDirection_.x() * gradStressX.boundaryField()[patchI])
260  -(forceDirection_.y() * gradStressY.boundaryField()[patchI])
261  -(forceDirection_.z() * gradStressZ.boundaryField()[patchI])
262  ) & nf
263  )
264  + (forceDirection_ & nf)*gradp.boundaryField()[patchI]
265  )
266  /denom();
267  }
268 }
269 
270 
272 {
273  for (const label patchI : forcePatches_)
274  {
275  const fvPatch& patch = mesh_.boundary()[patchI];
276  tmp<vectorField> tnf = patch.nf();
277  const vectorField& nf = tnf();
278  bdJdStressPtr_()[patchI] = (forceDirection_ * nf)/denom();
279  }
280 }
281 
282 
283 scalar objectiveForce::denom() const
284 {
285  return 0.5*UInf_*UInf_*Aref_;
286 }
287 
288 
290 {
291  return forceDirection_;
292 }
293 
294 
295 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
296 
297 } // End namespace objectives
298 } // End namespace Foam
299 
300 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::Vector::x
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
Foam::GeometricField::component
tmp< GeometricField< cmptType, PatchField, GeoMesh > > component(const direction) const
Return a component of the field.
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::fvc::grad
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
turbulence
Info<< "Reading field U\n"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);volScalarField rho(IOobject("rho", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), thermo.rho());volVectorField rhoU(IOobject("rhoU", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *U);volScalarField rhoE(IOobject("rhoE", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *(e+0.5 *magSqr(U)));surfaceScalarField pos(IOobject("pos", runTime.timeName(), mesh), mesh, dimensionedScalar("pos", dimless, 1.0));surfaceScalarField neg(IOobject("neg", runTime.timeName(), mesh), mesh, dimensionedScalar("neg", dimless, -1.0));surfaceScalarField phi("phi", fvc::flux(rhoU));Info<< "Creating turbulence model\n"<< endl;autoPtr< compressible::turbulenceModel > turbulence(compressible::turbulenceModel::New(rho, U, phi, thermo))
Definition: createFields.H:94
Foam::dimLength
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::objectives::objectiveForce::Aref_
scalar Aref_
Definition: objectiveForce.H:68
Foam::objectiveIncompressible::vars_
const incompressibleVars & vars_
Definition: objectiveIncompressible.H:62
Foam::objectives::addToRunTimeSelectionTable
addToRunTimeSelectionTable(objectiveIncompressible, objectiveForce, dictionary)
Foam::objectives::objectiveForce::update_dxdbMultiplier
void update_dxdbMultiplier()
Update delta(x)/delta b multiplier.
Definition: objectiveForce.C:201
Foam::objectives::objectiveForce
Definition: objectiveForce.H:58
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::gSum
Type gSum(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:594
objectiveForce.H
Foam::singlePhaseTransportModel
A simple single-phase transport model based on viscosityModel.
Definition: singlePhaseTransportModel.H:57
Foam::objectives::objectiveForce::stressZPtr_
autoPtr< volVectorField > stressZPtr_
Definition: objectiveForce.H:74
Foam::incompressibleVars::p
const volScalarField & p() const
Return const reference to pressure.
Definition: incompressibleVars.C:305
Foam::objectives::objectiveForce::update_dSdbMultiplier
void update_dSdbMultiplier()
Update delta(n dS)/delta b multiplier.
Definition: objectiveForce.C:175
Foam::Vector::z
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::incompressibleVars::RASModelVariables
const autoPtr< incompressible::RASModelVariables > & RASModelVariables() const
Return const reference to the turbulence model variables.
Definition: incompressibleVars.C:444
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
Foam::objectives::defineTypeNameAndDebug
defineTypeNameAndDebug(objectiveForce, 0)
Foam::Field< vector >
Foam::objective::J_
scalar J_
Objective function value and weight.
Definition: objective.H:77
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::fvPatch
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:65
Foam::incompressibleVars::laminarTransport
const singlePhaseTransportModel & laminarTransport() const
Return const reference to transport model.
Definition: incompressibleVars.C:418
Foam::T
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Definition: FieldFieldFunctions.C:58
Foam::objectives::objectiveForce::UInf_
scalar UInf_
Definition: objectiveForce.H:70
Foam::incompressibleVars::U
const volVectorField & U() const
Return const reference to velocity.
Definition: incompressibleVars.C:331
Foam::objectiveIncompressible
Abstract base class for objective functions in incompressible flows.
Definition: objectiveIncompressible.H:54
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::objectives::objectiveForce::stressXPtr_
autoPtr< volVectorField > stressXPtr_
Definition: objectiveForce.H:72
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
Foam::objectives::objectiveForce::forcePatches_
labelHashSet forcePatches_
Definition: objectiveForce.H:66
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:85
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::incompressibleVars::turbulence
const autoPtr< incompressible::turbulenceModel > & turbulence() const
Return const reference to the turbulence model.
Definition: incompressibleVars.C:431
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::fvMesh::boundary
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:685
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
Foam::objectives::objectiveForce::forceDirection
const vector & forceDirection() const
Return force direction.
Definition: objectiveForce.C:289
U
U
Definition: pEqn.H:72
Foam::Vector::y
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
Foam::objectives::objectiveForce::objectiveForce
objectiveForce(const fvMesh &mesh, const dictionary &dict, const word &adjointSolverName, const word &primalSolverName)
Construct from components.
Definition: objectiveForce.C:55
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
DebugInfo
#define DebugInfo
Report an information message using Foam::Info.
Definition: messageStream.H:382
Foam::objectives::objectiveForce::update_boundarydJdp
void update_boundarydJdp()
Update values to be added to the adjoint wall velocity.
Definition: objectiveForce.C:166
Foam::singlePhaseTransportModel::nu
virtual tmp< volScalarField > nu() const
Return the laminar viscosity.
Definition: singlePhaseTransportModel.C:72
Foam::foamVersion::patch
const std::string patch
OpenFOAM patch number as a std::string.
Foam::objectives::objectiveForce::J
scalar J()
Return the objective function value.
Definition: objectiveForce.C:122
Foam::Vector< scalar >
Foam::type
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:590
Foam::objectives::objectiveForce::denom
virtual scalar denom() const
Return denominator, without density.
Definition: objectiveForce.C:283
Foam::objective::mesh_
const fvMesh & mesh_
Definition: objective.H:67
Foam::wordRes
A List of wordRe with additional matching capabilities.
Definition: wordRes.H:51
Foam::objectives::objectiveForce::update_dJdStressMultiplier
void update_dJdStressMultiplier()
Update dJ/dStress multiplier.
Definition: objectiveForce.C:271
Foam::objectives::objectiveForce::stressYPtr_
autoPtr< volVectorField > stressYPtr_
Definition: objectiveForce.H:73
Foam::objectiveIncompressible::bdJdpPtr_
autoPtr< boundaryVectorField > bdJdpPtr_
Adjoint (intlet,wall) velocity.
Definition: objectiveIncompressible.H:88
Foam::GeometricField< scalar, fvPatchField, volMesh >
Foam::objectives::objectiveForce::forceDirection_
vector forceDirection_
Definition: objectiveForce.H:67
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:62
Foam::incompressibleVars::pInst
const volScalarField & pInst() const
Return const reference to pressure.
Definition: incompressibleVars.C:382
Foam::fvMesh::Sf
const surfaceVectorField & Sf() const
Return cell face area vectors.
Definition: fvMeshGeometry.C:319