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

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

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

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

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

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.

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;
}



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

new version

- now writes velocity field as vector

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

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

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

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.

Thanks Eric! This is gonna be helpful :)

could you remind me how to compile it?
Caroline

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 ?

I'm sorry for reopening an old post, but has someone managed to compile the cartToCyl utility for OpenFOAM-2.x.x (I'm interested in 2.2 and later versions)?

Thanks

V.

Ok, my fault. It actually compiles without a problem on 2.2.x, needs only minor modifications in 2.3.x and 2.4.x.

V.

Hello Krastev,

I have difficulty compiling it in OpenFOAM 4x. I am new to OpenFOAM. Please, it would be helpful if you state the modifications those are needed for the later versions as you have mentioned.

Thanks,
Padmini Shivaraj

Hello Krastev,

I have difficulty compiling it in OpenFOAM 4x. I am new to OpenFOAM. Please, it would be helpful if you state the modifications those are needed for the later versions as you have mentioned.

Thanks,
Padmini Shivaraj

cartToCyl-Uz for OpenFOAM4.0
 

Hello Foamers,

it is very impressive to see the beautiful coding with tensor operation here!

While compiling this utility in OpenFOAM4.0, I got the error that
......
cartToCyl.C:98:24: error: no match for ‘operator=’ (operand types are ‘const Internal {aka const Foam::DimensionedField<Foam::Vector<double>, Foam::volMesh>}’ and ‘Foam::tmp<Foam::Field<Foam::Vector<double> > >’)
ccCyl.internalField() = cyl.localVector(cc.internalField());^

In file included from /opt/openfoam4/src/OpenFOAM/lnInclude/DimensionedField.H:366:0
......

I checked the related files but couldn't figure it out as new programmer to OpenFOAM. I would appreciate greatly if anyone can help to update this utility to OpenFOAM4.0.

Thank you very much in advance!

sitajeje

cartToCyl-Uz for OpenFOAM4.0
 

Dear Foamers,

I compared the GeometricField.H of OF4.1 and OF2.1, and found that the type of internalField() used in this utility is changed:
OpenFOAM4.1
typedef DimensionedField<Type, GeoMesh> Internal;
....
inline const Internal& internalField() const;

OpenFOAM2.1
typedef Field<Type> InternalField;
....
inline const InternalField& internalField() const;

This made the error that I posted in the last thread. I wonder how should I amend this code so that it can be used for OF4.0? I appreciate greatly for any suggestions and hints!

Thank you very much in advance!

Best regards,
sitajeje