CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Wall shear stress and cell area code in SimpleFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   July 8, 2017, 05:19
Post Wall shear stress and cell area code in SimpleFoam
  #1
New Member
 
Francisco Almeida
Join Date: Jul 2017
Posts: 1
Rep Power: 0
Fran6 is on a distinguished road
Hi all,

this is my first post. I have been reading a lot about wall shear stress (WSS), I need to calculate it for my thesis, however, I want to obtain the values for wss and area to calculate the average WSS on an artery.

I found different codes looking on previous posts like this,

Improve simpleFoam convergence

write cell volumes

I modified the code (from previous posts) to add the cell area and their locations in the mesh (x y z) but when I try to add the wall shear stress I am unable to. Here is the code I am using, it doesn't contain the WSS yet, I compiled a new simpleFoam solver to run this .H file

Quote:
{
Info << "Writing Area of the cells on the patch!" << endl;

word wallPatchName = "walls";
Info << "Searching for patch " << wallPatchName << endl;

// Check, if the patch name exists:
label patchID = mesh.boundaryMesh().findPatchID(wallPatchName);
if (patchID < 0)
{
Info << "No patch found with name "" << wallPatchName << ""!"
<< endl;
}

// Pointer to the patch we want to have shear stress for:
const polyPatch& cPatch = mesh.boundaryMesh()[patchID];
label nFaces = cPatch.size();

pointField faceAreas(nFaces, vector(0,0,0));

label runFaces = 0;

// collecting the interesting data
for (int i=0; i<nFaces; i++)
{
faceAreas[runFaces] = mesh.Sf().boundaryField()[patchID][i];
runFaces++;
}

// for parallel
reduce(nFaces, sumOp<label>());
Info << "Number of global faces: " << nFaces << endl;

// make a cute filename
fileName casePath = cwd();//casePath.name();
fileName caseName = casePath.name();
fileName myFile = caseName+"_wallShearStress.csv";
Info << "Data is written to file " << myFile << endl;

// create a new file
OFstream *myStream = NULL;
if (Pstream::master())
{
myStream = new OFstream(myFile);
}

// slaves send data to master, master writes data to file
if (Pstream:arRun())
{
if (Pstream::myProcNo() != Pstream::masterNo())
{
OPstream toMaster(Pstream::blocking, Pstream::masterNo());
toMaster << faceAreas;
}
else // master part
{
// first write own data
for (int i=0; i<faceAreas.size(); i++)
{
*myStream << i << ";"
<< mag(faceAreas[i]) << endl;

}

label run = faceAreas.size();

// then receive slave data and write
pointField bufferFaceAreas;

for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
{
IPstream fromSlave(Pstream::blocking, slave);
fromSlave >> bufferFaceAreas;

for (int i=0; i<bufferFaceAreas.size(); i++)
{
*myStream << run << "; "
<< mag(bufferFaceAreas[i]) << endl;
run++;
}
}
}
}
else // single processor only
{
for (int i=0; i<nFaces; i++)
{
*myStream << i << "; "
// coordinates on face centers:
<< mesh.Cf().boundaryField()[patchID][i].component(0) << "; "
<< mesh.Cf().boundaryField()[patchID][i].component(1) << "; "
<< mesh.Cf().boundaryField()[patchID][i].component(2) << "; "
<< mag(mesh.Sf().boundaryField()[patchID][i]) << endl;
}
}

}
when I try to use the following code I can't compile it, something is wrong with the wall shear stress

Quote:
#include "RASModel.H"
{
Info << "Writing wall shear stress for matlab!" << endl;

word wallPatchName = "bottom";
Info << "Searching for patch " << wallPatchName << endl;

// Check, if the patch name exists:
label patchID = mesh.boundaryMesh().findPatchID(wallPatchName);
if (patchID < 0)
{
Info << "No patch found with name "" << wallPatchName << ""!"
<< endl;
}
else
{
// calculating tau

# include "createFields2.H"
mesh.readUpdate();
volSymmTensorField Reff(RASModel->devReff());
volVectorField tau
(
IOobject
(
"tau",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector
(
"tau",
Reff.dimensions(),
vector::zero
)
);

forAll(tau.boundaryField(), patchi)
{
tau.boundaryField()[patchi] =
(
-mesh.Sf().boundaryField()[patchi]
/mesh.magSf().boundaryField()[patchi]
) & Reff.boundaryField()[patchi];
}



// Pointer to the patch we want to have shear stress for:
const polyPatch& cPatch = mesh.boundaryMesh()[patchID];
label nFaces = cPatch.size();

pointField faceCenterCoords(nFaces, vector(0,0,0));
pointField faceAreas(nFaces, vector(0,0,0));
pointField wallShearStress(nFaces, vector(0,0,0));

label runFaces = 0;

// collecting the interesting data
for (int i=0; i<nFaces; i++)
{
faceCenterCoords[runFaces] = mesh.Cf().boundaryField()[patchID][i];
faceAreas[runFaces] = mesh.Sf().boundaryField()[patchID][i];
wallShearStress[runFaces] = tau.boundaryField()[patchID][i];
runFaces++;
}

// for parallel
reduce(nFaces, sumOp<label>());
Info << "Number of global faces: " << nFaces << endl;


// make a cute filename
fileName casePath = cwd();//casePath.name();
fileName caseName = casePath.name();
fileName myFile = "tau.dat";
Info << "Data is written to file " << myFile << endl;

// create a new file
OFstream *myStream = NULL;
if (Pstream::master())
{
myStream = new OFstream(myFile);
}

// slaves send data to master, master writes data to file
if (Pstream:arRun())
{
if (Pstream::myProcNo() != Pstream::masterNo())
{
OPstream toMaster(Pstream::blocking, Pstream::masterNo());
toMaster << faceCenterCoords;
toMaster << faceAreas;
toMaster << wallShearStress;
}
else // master part
{
// first write own data
for (int i=0; i<faceCenterCoords.size(); i++)
{
*myStream << i << ";"
<< faceCenterCoords[i].component(0) << "; "
<< faceCenterCoords[i].component(1) << "; "
<< faceCenterCoords[i].component(2) << "; "

<< wallShearStress[i].component(0) << "; "
<< wallShearStress[i].component(1) << "; "
<< wallShearStress[i].component(2) << "; "

<< mag(wallShearStress[i]) << endl;
}

label run = faceCenterCoords.size();

// then receive slave data and write
pointField bufferFaceCenterCoords;
pointField bufferFaceAreas;
pointField bufferWallShearStress;


for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
{
IPstream fromSlave(Pstream::blocking, slave);
fromSlave >> bufferFaceCenterCoords;
fromSlave >> bufferFaceAreas;
fromSlave >> bufferWallShearStress;

for (int i=0; i<bufferFaceCenterCoords.size(); i++)
{
*myStream << run << ";"
<< bufferFaceCenterCoords[i].component(0) << "; "
<< bufferFaceCenterCoords[i].component(1) << "; "
<< bufferFaceCenterCoords[i].component(2) << "; "

<< bufferWallShearStress[i].component(0) << "; "
<< bufferWallShearStress[i].component(1) << "; "
<< bufferWallShearStress[i].component(2) << "; "

<< mag(bufferWallShearStress[i]) << endl;
run++;
}
}
}
}
else // single processor only
{
for (int i=0; i<nFaces; i++)
{
*myStream << i << "; "
// coordinates on face centers:
<< mesh.Cf().boundaryField()[patchID][i].component(0) << "; "
<< mesh.Cf().boundaryField()[patchID][i].component(1) << "; "
<< mesh.Cf().boundaryField()[patchID][i].component(2) << "; "

// wall shear stress as a vector:
<< tau.boundaryField()[patchID][i].component(0) << "; "
<< tau.boundaryField()[patchID][i].component(1) << "; "
<< tau.boundaryField()[patchID][i].component(2) << "; "

// magnitude of wall shear stress:
<< mag(tau.boundaryField()[patchID][i]) << endl;
}
}

}
}
I get this error

Quote:
In file included from sempleFoam.C:61:0:
write.H:19:34: fatal error: createFields2.H: No such file or directory
compilation terminated.
/opt/openfoam4/wmake/rules/General/transform:8: recipe for target 'Make/linux64GccDPInt32Opt/sempleFoam.o' failed
changing that createfields2.H to createfields.H i get the following,

Quote:
In file included from sempleFoam.C:61:0:
write.H: In function ‘int main(int, char**)’:
write.H:21:41: error: missing template arguments before ‘->’ token
volSymmTensorField Reff(RASModel->devReff());
^
write.H:43:41: error: passing ‘const Foam::fvPatchField<Foam::Vector<double> >’ as ‘this’ argument discards qualifiers [-fpermissive]
tau.boundaryField()[patchi] =
^
In file included from /opt/openfoam4/src/finiteVolume/lnInclude/fvPatchField.H:592:0,
from /opt/openfoam4/src/finiteVolume/lnInclude/volFields.H:40,
from /opt/openfoam4/src/finiteVolume/lnInclude/surfaceInterpolationScheme.C:27,
from /opt/openfoam4/src/finiteVolume/lnInclude/surfaceInterpolationScheme.H:304,
from /opt/openfoam4/src/finiteVolume/lnInclude/surfaceInterpolate.H:41,
from /opt/openfoam4/src/finiteVolume/lnInclude/fvc.H:39,
from /opt/openfoam4/src/finiteVolume/lnInclude/fvCFD.H:8,
from sempleFoam.C:33:
/opt/openfoam4/src/finiteVolume/lnInclude/fvPatchField.C:395:6: note: in call to ‘void Foam::fvPatchField<Type>:perator=(const Foam::UList<T>&) [with Type = Foam::Vector<double>]’
void Foam::fvPatchField<Type>:perator=
^
/opt/openfoam4/wmake/rules/General/transform:8: recipe for target 'Make/linux64GccDPInt32Opt/sempleFoam.o' failed
I just want to obtain the following, cell area, wall shear stress, location of these two and a way to check them in paraView. At the moment I can generate a .csv file with the area and location (thanks to the previous posts)

Is there a way to fix this or to use the existing wallShearStress function in this code?, I am carrying out a Laminar flow simulation on an artery, so turbulence is not important in this case, thanks for your help,

Francisco Almeida

Last edited by Fran6; July 8, 2017 at 05:26. Reason: title change
Fran6 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Wall Shear Stress Sampling Probes Mahe88 OpenFOAM Post-Processing 0 December 14, 2016 11:09
Wall Shear Stress prediction varies with fix y+ Jeeloong FLUENT 0 November 24, 2016 18:17
Exporting wall shear stress in relative frame of reference gotang FLUENT 0 September 15, 2014 10:35
Natural convection in a closed domain STILL NEEDING help! Yr0gErG FLUENT 3 June 12, 2013 02:12
Re: About wall shear stress Mike FLUENT 9 November 17, 2003 15:41


All times are GMT -4. The time now is 12:27.