CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   bug in cylindricalCS class or fatal setup error?!! (http://www.cfd-online.com/Forums/openfoam/70252-bug-cylindricalcs-class-fatal-setup-error.html)

florian_krause November 18, 2009 11:18

bug in cylindricalCS class or fatal setup error?!!
 
Hi all,

a short while ago, I could finally transform my velocity volVectorField from cartesian to cylindrical coordinates ( thanks to Hrv !! ). I made up a post utility, which basically looks like the code below.

info for my case: LES pipe flow with approx. centreline velocity (in axial pipe / global x-direction) of 1,4m/s when its more or less fully developed.

int main(int argc, char *argv[])
{
timeSelector::addOptions();
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createMesh.H"

forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
mesh.readUpdate();
// defining cylindrical coordinate system

Info<< " creating cylindrical system (r, teta, z) from cartesian system with r = z and z = x" << endl;

cylindricalCS cyl("cylindricalCS", vector(0, 0, 0), vector(1, 0, 0), vector(0, 0, 1));

// transformation of velocity field U from cartesian -> cylindrical

IOobject Uheader
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ
);

// Check U exists
if (Uheader.headerOk())
{
mesh.readUpdate();

Info<< " Reading U" << endl;
volVectorField U(Uheader, mesh);

volVectorField Ucyl
(
IOobject
(
"Ucyl",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector
(
"Ucyl",
Ucyl.dimensions(),
vector::zero
)
);

Info<< " converting U" << endl;

Ucyl.internalField() = cyl.localVector(U.internalField());
forAll (Ucyl.boundaryField(), patchI)
{
Ucyl.boundaryField()[patchI] = cyl.localVector(U.boundaryField()[patchI]);
}

Ucyl.write();
}
else
{
Info<< " No exisiting U field" << endl;
}


What I get out are my velocity components in cyl. CS, lets say (r, phi, z) with z the axial / previous x direction.

As you can see below, the Uphi (tangential velocity) is totally wrong (max Uphi=175m/s) and like half the speed of sound, so I dont know if I should start crying or laughing..... :confused:

Does anyone know if this is a bug or if I did something wrong in the use of the cylindricalCS class ?! How can I fix this?!

dimensions [0 1 -1 0 0 0 0];
internalField nonuniform List<vector>
240000
(
(0.0717613 147.295 1.18858)
(0.0857317 147.196 1.17645)
(0.129587 129.416 1.19515)
....
....


The other problem is, that, if I am not mistaken, with the class I can only convert vector or vectorFields. But what about symm tensors?!?! I would like to transform my symm Reynolds stress tensor.

Does anyone did this before or know how to do such a tensor transformation?!

as always, any help is appreciated!

Best,
Florian

florian_krause April 30, 2010 08:11

Hello,

I want to bring up again my issue, when converting my velocity field from cartesian to cylindrical coordinate system using the cylindricalCS class (see previous post). I am still using OF-1.6.x.

The geometry is a pipe with the global x-axis being the pipe-axis with (0 0 0) being the midpoint/origin of my inlet patch.

I define my cyl. coordinate system as follows:

cylindricalCS cyl("cylindricalCS", vector(0, 0, 0), vector(1, 0, 0), vector(0, 0, 1));

Any ideas what might be wrong with my utility and my cyl. coordinate defintion or how to correctly define a coresponding cyl. coordinate system?

any hint is much appreciated!

Best,
Florian

stevenvanharen November 10, 2010 12:33

Hi Florian,

Did you fix this?

I have tested your utility and it seems the conversion converts the values as if the cartesian velocity vector is just a position vector.

I will combine this with the cell centres utility and then it should not be to hard to calculate the cylindrical velocities using my favourite analytical mechanics book.

Please let me know if you still need this then I can share this with you?

Regards,

Steven

florian_krause November 11, 2010 04:35

Hi Steven,

I never completely fixed this one, mainly due to time. But I would be still interested if you get it working properly... seeing the results of a corresponding thread search, it is or at least it was of general interest.

thanks,
Florian

stevenvanharen November 11, 2010 04:37

It is on my to do list for today, I'll post it here when it is working.

stevenvanharen November 11, 2010 09:17

1 Attachment(s)
Ok, it is working but I am not happy yet.

Be very careful using it! It only works if the x-axis in the cartesian reference frame is aligned with the z-axis of the cylindrical reference frame.

If somebody can show me how to create a volTensorField (or volVectorField) from volScalarField I can improve this quite a lot quite easily.

For now it writes out the three component u_r, u_theta and u_z in different files as volScalarFields. I know this is shit, but if someone helps me with the problem mentioned above I can improve this.

nimasam December 16, 2010 06:07

i look at above code and made some changes i think it works correctly but no guarantee :D now it prints out all values in one volVectorField.
\*---------------------------------------------------------------------------*/
Application
cartToCyl

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "cylindricalCS.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
timeSelector::addOptions();

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

instantList timeDirs = timeSelector::select0(runTime, args);

#include "createMesh.H"

Info<< " Convert cartesian system (x, y, z) to cylindrical system (r, teta, z)" <<endl;

forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
mesh.readUpdate();

cylindricalCS ccs
("ccs",
point (0, 0, 0), // center point of ccs
vector(0, 0, 1), // axis of ccs
vector(1, 0, 0),// base axis for cylindrical angle
false
); //keyword false: Use radians since cos and sin work with radians

//Read the vectors pointing from the origin of the mesh to the cell centers

volVectorField cc = mesh.C();

//The vector field cc will now be transformed into cylindrical coordinates
volVectorField ccCyl
(
IOobject
(
"ccCyl",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
vector(0,0,0)
);

ccCyl.internalField() = ccs.localVector(cc.internalField());

forAll (ccCyl.boundaryField(), patchI)
{
ccCyl.boundaryField()[patchI] = ccs.localVector(cc.boundaryField()[patchI]);
}
//cc.write();
//ccCyl.write();



volScalarField theta = ccCyl.component(vector::Y);
//now we have the all important theta of all cell centre postitions,
//from which we can construct the following transformation tensor:
/*
[cos(theta) sin(theta) 0
-sin(theta) cos(theta) 0
0 0 1]
*/


// transformation of velocity field U from cartesian -> cylindrical
IOobject Uheader
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ
);

// Check U exists
if (Uheader.headerOk())
{
mesh.readUpdate();

Info<< " Reading U" << endl;
volVectorField U(Uheader, mesh);

Info<< " converting U" << endl;

volVectorField Ucyl
(
IOobject
(
"Ucyl",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector("Ucyl",dimensionSet(0,1,-1,0,0,0,0),vector::zero)
);

Ucyl.replace(vector::X, (U.component(vector::X)*cos(theta))+(U.component(v ector::Y)*sin(theta)));//Ur
Ucyl.replace(vector::Y, (-1*U.component(vector::X)*sin(theta))+(U.component( vector::Y)*cos(theta)));//Uteta
Ucyl.replace(vector::Z, U.component(vector::Z)); //Uz
Ucyl.write();

}
else
{
Info<< " No U field exist" << endl;
}
}
return 0;
}



// ************************************************** *********************** //

stevenvanharen December 20, 2010 09:51

1 Attachment(s)
new version

- now writes velocity field as vector

nimasam December 20, 2010 10:26

hi steven
why dont u use "replace" command its easier ;) and some how more pretty!!!

stevenvanharen December 20, 2010 11:13

hi nimasam,

I already changed the code into its current form before your reply. That is why I did not bother to change it. I did not intentionally avoid the replace command.

About the beauty of programming, maybe when you are in the Netherlands we can have a beer because the discussion about the beauty of programming is never ending I think :).

Cheers,

Steven

CaroVandame April 25, 2012 15:23

Hi to all

I'm relatively new to OpenFoam, and especially to the programming part of it. Could someone explain how to use the code that was written above? where and how to compile it for ex?

Also if I understood right, this code transforms the cartesian x-axis into the cylindrical z-axis. Is there a way to change the cartesian system so that if the mesh axis is not aligned with x, we could change the coordinates? like a permutation or something like that?

Thanks, I really appreciate any help

Caroline

tunkers April 26, 2012 19:13

1 Attachment(s)
Hi, A while ago I modified the program posted above by Steven to do what you are asking. The Ux, Uy, Uz from the cartesian XYZ coordinate system transform into Ur, Utheta, Uz in cylindrical coordinates where the Uz axes from both coordinate systems are aligned.

CaroVandame April 26, 2012 19:22

Thanks Eric! This is gonna be helpful :)

could you remind me how to compile it?
Caroline

eagle December 3, 2012 15:38

Hi all !

Im trying to use this utilitie for a body not aligned with the origin of the mesh. Its because Im simulating a bent pipe and the center of bend is the center of origin and not the rotation axis of the pipe.

The Problem is if Im changing the cylindricalCs like

cylindricalCS cyl(

"cylindricalCS",
-> point(0.169, 0, 0),//origin of cs
vector(0, 0, 1), //rotation axis
vector(0, 1, 0), //start of angle
false); //keyword false: Use radians since cos and sin work with radians

The class/or tool doesn't use the point Im giving it. Im relatively new to foam so Im not
sure if the problem is somewhere around here:

ccCyl.internalField() = cyl.localVector(cc.internalField());
forAll (ccCyl.boundaryField(), patchI)
{
ccCyl.boundaryField()[patchI] = cyl.localVector(cc.boundaryField()[patchI]);
}

I checked Theta and its not in accordance to the center of the cylindrical Coordinatesystem.

I know I could just change the origin of the mesh. But as I'm also simulating a pipe helix and need the cylindrical coordinates somewhere in the helix I do need the answer to this anyway.

I could really use some help here ;)

best regards

-----------------------------------------------------

Edit:

by searching for a solution I considered using:

globalToLocal (const vector &, bool translate) const

instead of:

localVector (const vector &global) const


see here:

http://foam.sourceforge.net/docs/cpp/a00294.html

So Im in the opinion during transformation with localVector the translation of the origin is getting lost.
I could imagine using globalToLocal is solving the problem but I cant use it or dont know how because its
a protected function ?


All times are GMT -4. The time now is 13:38.