Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding post-processing routines to compute surface forces #659

Merged
merged 2 commits into from
Jul 3, 2024
Merged

Adding post-processing routines to compute surface forces #659

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3ca3f96
Adding new post-processing routines to compute surface forces
bernardopacini Jul 2, 2024
3206637
Replacing previous version of forcePerS
bernardopacini Jul 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
172 changes: 83 additions & 89 deletions src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,11 @@

#include "fvCFD.H"
#include "argList.H"
#include "autoPtr.H"
#include "Time.H"
#include "timeSelector.H"
#include "TimePaths.H"
#include "ListOps.H"
#include "fvMesh.H"
#include "OFstream.H"
#include "simpleControl.H"
Expand All @@ -24,40 +28,31 @@ using namespace Foam;

int main(int argc, char* argv[])
{
Info << "Computing forcePerS...." << endl;
Info << "Computing forces...." << endl;

argList::addOption(
"patchNames",
"'(inlet)'",
"'(wing)'",
"List of patch names to compute");

argList::addOption(
"forceDir",
"'(0 0 1)'",
"Force direction");

argList::addOption(
"time",
"1000",
"Tme instance to compute");
"Time instance to compute, if not provided runs all times");

#include "setRootCase.H"
#include "createTime.H"

scalar time;
if (args.optionFound("time"))
if (!args.optionFound("time"))
{
time = readScalar(args.optionLookup("time")());
Info << "Time not set, running all times." << endl;
}
else
{
Info << "time not set! Exit." << endl;
return 1;
}
runTime.setTime(time, 0);

#include "createMesh.H"
#include "createFields.H"
// Create the processor databases
autoPtr<TimePaths> timePaths;
timePaths = autoPtr<TimePaths>::New(args.rootPath(), args.caseName());

const instantList timeDirs(timeSelector::select(timePaths->times(), args));

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

List<scalar> forceDir1;
if (args.optionFound("forceDir"))
{
forceDir1 = scalarList(args.optionLookup("forceDir")());
}
else
forAll(timeDirs, iTime)
{
Info << "forceDir not set! Exit." << endl;
return 1;
}
vector forceDir(vector::zero);
forceDir.x() = forceDir1[0];
forceDir.y() = forceDir1[1];
forceDir.z() = forceDir1[2];

volScalarField forcePerS(
IOobject(
"forcePerS",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE),
mesh,
dimensionedScalar("forcePerS", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
fixedValueFvPatchScalarField::typeName);

// this code is pulled from:
// src/functionObjects/forcces/forces.C
// modified slightly
vector forces(vector::zero);

const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField();
const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField();

volSymmTensorField devRhoReff(
IOobject(
IOobject::groupName("devRhoReff", U.group()),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE),
(-rho * nuEff) * dev(twoSymm(fvc::grad(U))));
runTime.setTime(timeDirs[iTime].value(), 0);

const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField();
#include "createMesh.H"
#include "createFields.H"

forAll(patchNames, cI)
{
// get the patch id label
label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]);
if (patchI < 0)
{
Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl;
return 1;
}
// create a shorter handle for the boundary patch
const fvPatch& patch = mesh.boundary()[patchI];
// normal force
vectorField fN(
Sfb[patchI] * p.boundaryField()[patchI]);
// tangential force
vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
// sum them up
forAll(patch, faceI)
// Initialize surface field for face-centered forces
volVectorField forcePerS(
IOobject(
"forcePerS",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE),
mesh,
dimensionedVector("forcePerS", dimensionSet(1, -1, -2, 0, 0, 0, 0), vector::zero),
fixedValueFvPatchScalarField::typeName);

// this code is pulled from:
// src/functionObjects/forcces/forces.C
// modified slightly
vector forces(vector::zero);

const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField();
const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField();

volSymmTensorField devRhoReff(
IOobject(
IOobject::groupName("devRhoReff", U.group()),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE),
(-rho * nuEff) * dev(twoSymm(fvc::grad(U))));

const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField();

vector totalForce(vector::zero);
forAll(patchNames, cI)
{
forces.x() = fN[faceI].x() + fT[faceI].x();
forces.y() = fN[faceI].y() + fT[faceI].y();
forces.z() = fN[faceI].z() + fT[faceI].z();
scalar force = forces & forceDir;
forcePerS.boundaryFieldRef()[patchI][faceI] = force / magSfb[patchI][faceI];
// get the patch id label
label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]);
if (patchI < 0)
{
Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl;
return 1;
}
// create a shorter handle for the boundary patch
const fvPatch& patch = mesh.boundary()[patchI];
// normal force
vectorField fN(Sfb[patchI] * p.boundaryField()[patchI]);
// tangential force
vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
// sum them up
forAll(patch, faceI)
{
// compute forces
forces.x() = fN[faceI].x() + fT[faceI].x();
forces.y() = fN[faceI].y() + fT[faceI].y();
forces.z() = fN[faceI].z() + fT[faceI].z();

// project force direction
forcePerS.boundaryFieldRef()[patchI][faceI] = forces / magSfb[patchI][faceI];

totalForce.x() += forces.x();
totalForce.y() += forces.y();
totalForce.z() += forces.z();
}
}
}
forcePerS.write();

Info << "Force: " << forces << endl;
forcePerS.write();

Info << "Computing forcePerS.... Completed!" << endl;
Info << "Total force: " << totalForce << endl;

Info << "Computing forces.... Completed!" << endl;
}
return 0;
}

Expand Down
Loading
Loading