Skip to content

Commit 7441216

Browse files
Adding post-processing routines to compute surface forces (#659)
* Adding new post-processing routines to compute surface forces * Replacing previous version of forcePerS
1 parent cd2b115 commit 7441216

File tree

2 files changed

+177
-174
lines changed

2 files changed

+177
-174
lines changed

src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C

Lines changed: 83 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,11 @@
1111

1212
#include "fvCFD.H"
1313
#include "argList.H"
14+
#include "autoPtr.H"
1415
#include "Time.H"
16+
#include "timeSelector.H"
17+
#include "TimePaths.H"
18+
#include "ListOps.H"
1519
#include "fvMesh.H"
1620
#include "OFstream.H"
1721
#include "simpleControl.H"
@@ -24,40 +28,31 @@ using namespace Foam;
2428

2529
int main(int argc, char* argv[])
2630
{
27-
Info << "Computing forcePerS...." << endl;
31+
Info << "Computing forces...." << endl;
2832

2933
argList::addOption(
3034
"patchNames",
31-
"'(inlet)'",
35+
"'(wing)'",
3236
"List of patch names to compute");
3337

34-
argList::addOption(
35-
"forceDir",
36-
"'(0 0 1)'",
37-
"Force direction");
38-
3938
argList::addOption(
4039
"time",
4140
"1000",
42-
"Tme instance to compute");
41+
"Time instance to compute, if not provided runs all times");
4342

4443
#include "setRootCase.H"
4544
#include "createTime.H"
4645

47-
scalar time;
48-
if (args.optionFound("time"))
46+
if (!args.optionFound("time"))
4947
{
50-
time = readScalar(args.optionLookup("time")());
48+
Info << "Time not set, running all times." << endl;
5149
}
52-
else
53-
{
54-
Info << "time not set! Exit." << endl;
55-
return 1;
56-
}
57-
runTime.setTime(time, 0);
5850

59-
#include "createMesh.H"
60-
#include "createFields.H"
51+
// Create the processor databases
52+
autoPtr<TimePaths> timePaths;
53+
timePaths = autoPtr<TimePaths>::New(args.rootPath(), args.caseName());
54+
55+
const instantList timeDirs(timeSelector::select(timePaths->times(), args));
6156

6257
List<wordRe> patchNames;
6358
if (args.optionFound("patchNames"))
@@ -70,83 +65,82 @@ int main(int argc, char* argv[])
7065
return 1;
7166
}
7267

73-
List<scalar> forceDir1;
74-
if (args.optionFound("forceDir"))
75-
{
76-
forceDir1 = scalarList(args.optionLookup("forceDir")());
77-
}
78-
else
68+
forAll(timeDirs, iTime)
7969
{
80-
Info << "forceDir not set! Exit." << endl;
81-
return 1;
82-
}
83-
vector forceDir(vector::zero);
84-
forceDir.x() = forceDir1[0];
85-
forceDir.y() = forceDir1[1];
86-
forceDir.z() = forceDir1[2];
87-
88-
volScalarField forcePerS(
89-
IOobject(
90-
"forcePerS",
91-
runTime.timeName(),
92-
mesh,
93-
IOobject::NO_READ,
94-
IOobject::NO_WRITE),
95-
mesh,
96-
dimensionedScalar("forcePerS", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
97-
fixedValueFvPatchScalarField::typeName);
98-
99-
// this code is pulled from:
100-
// src/functionObjects/forcces/forces.C
101-
// modified slightly
102-
vector forces(vector::zero);
103-
104-
const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField();
105-
const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField();
106-
107-
volSymmTensorField devRhoReff(
108-
IOobject(
109-
IOobject::groupName("devRhoReff", U.group()),
110-
mesh.time().timeName(),
111-
mesh,
112-
IOobject::NO_READ,
113-
IOobject::NO_WRITE),
114-
(-rho * nuEff) * dev(twoSymm(fvc::grad(U))));
70+
runTime.setTime(timeDirs[iTime].value(), 0);
11571

116-
const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField();
72+
#include "createMesh.H"
73+
#include "createFields.H"
11774

118-
forAll(patchNames, cI)
119-
{
120-
// get the patch id label
121-
label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]);
122-
if (patchI < 0)
123-
{
124-
Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl;
125-
return 1;
126-
}
127-
// create a shorter handle for the boundary patch
128-
const fvPatch& patch = mesh.boundary()[patchI];
129-
// normal force
130-
vectorField fN(
131-
Sfb[patchI] * p.boundaryField()[patchI]);
132-
// tangential force
133-
vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
134-
// sum them up
135-
forAll(patch, faceI)
75+
// Initialize surface field for face-centered forces
76+
volVectorField forcePerS(
77+
IOobject(
78+
"forcePerS",
79+
runTime.timeName(),
80+
mesh,
81+
IOobject::NO_READ,
82+
IOobject::NO_WRITE),
83+
mesh,
84+
dimensionedVector("forcePerS", dimensionSet(1, -1, -2, 0, 0, 0, 0), vector::zero),
85+
fixedValueFvPatchScalarField::typeName);
86+
87+
// this code is pulled from:
88+
// src/functionObjects/forcces/forces.C
89+
// modified slightly
90+
vector forces(vector::zero);
91+
92+
const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField();
93+
const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField();
94+
95+
volSymmTensorField devRhoReff(
96+
IOobject(
97+
IOobject::groupName("devRhoReff", U.group()),
98+
mesh.time().timeName(),
99+
mesh,
100+
IOobject::NO_READ,
101+
IOobject::NO_WRITE),
102+
(-rho * nuEff) * dev(twoSymm(fvc::grad(U))));
103+
104+
const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField();
105+
106+
vector totalForce(vector::zero);
107+
forAll(patchNames, cI)
136108
{
137-
forces.x() = fN[faceI].x() + fT[faceI].x();
138-
forces.y() = fN[faceI].y() + fT[faceI].y();
139-
forces.z() = fN[faceI].z() + fT[faceI].z();
140-
scalar force = forces & forceDir;
141-
forcePerS.boundaryFieldRef()[patchI][faceI] = force / magSfb[patchI][faceI];
109+
// get the patch id label
110+
label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]);
111+
if (patchI < 0)
112+
{
113+
Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl;
114+
return 1;
115+
}
116+
// create a shorter handle for the boundary patch
117+
const fvPatch& patch = mesh.boundary()[patchI];
118+
// normal force
119+
vectorField fN(Sfb[patchI] * p.boundaryField()[patchI]);
120+
// tangential force
121+
vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
122+
// sum them up
123+
forAll(patch, faceI)
124+
{
125+
// compute forces
126+
forces.x() = fN[faceI].x() + fT[faceI].x();
127+
forces.y() = fN[faceI].y() + fT[faceI].y();
128+
forces.z() = fN[faceI].z() + fT[faceI].z();
129+
130+
// project force direction
131+
forcePerS.boundaryFieldRef()[patchI][faceI] = forces / magSfb[patchI][faceI];
132+
133+
totalForce.x() += forces.x();
134+
totalForce.y() += forces.y();
135+
totalForce.z() += forces.z();
136+
}
142137
}
143-
}
144-
forcePerS.write();
145-
146-
Info << "Force: " << forces << endl;
138+
forcePerS.write();
147139

148-
Info << "Computing forcePerS.... Completed!" << endl;
140+
Info << "Total force: " << totalForce << endl;
149141

142+
Info << "Computing forces.... Completed!" << endl;
143+
}
150144
return 0;
151145
}
152146

0 commit comments

Comments
 (0)