
[Sponsors] 
November 18, 2009, 11:18 
bug in cylindricalCS class or fatal setup error?!!

#1 
Senior Member
Florian Krause
Join Date: Mar 2009
Location: Munich
Posts: 103
Rep Power: 8 
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 xdirection) 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..... 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 

April 30, 2010, 08:11 

#2 
Senior Member
Florian Krause
Join Date: Mar 2009
Location: Munich
Posts: 103
Rep Power: 8 
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 OF1.6.x. The geometry is a pipe with the global xaxis being the pipeaxis 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 

November 10, 2010, 12:33 

#3 
Senior Member
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 7 
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 

November 11, 2010, 04:35 

#4 
Senior Member
Florian Krause
Join Date: Mar 2009
Location: Munich
Posts: 103
Rep Power: 8 
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 

November 11, 2010, 04:37 

#5 
Senior Member
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 7 
It is on my to do list for today, I'll post it here when it is working.


November 11, 2010, 09:17 

#6 
Senior Member
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 7 
Ok, it is working but I am not happy yet.
Be very careful using it! It only works if the xaxis in the cartesian reference frame is aligned with the zaxis 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. 

December 16, 2010, 06:07 

#7 
Senior Member

i look at above code and made some changes i think it works correctly but no guarantee 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; } // ************************************************** *********************** // 

December 20, 2010, 09:51 

#8 
Senior Member
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 7 
new version
 now writes velocity field as vector 

December 20, 2010, 11:13 

#10 
Senior Member
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 7 
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 

April 25, 2012, 15:23 

#11 
New Member
Caroline Vandame
Join Date: Aug 2010
Posts: 26
Rep Power: 7 
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 xaxis into the cylindrical zaxis. 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 

April 26, 2012, 19:13 

#12 
New Member
Eric
Join Date: Aug 2010
Posts: 14
Rep Power: 7 
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.


April 26, 2012, 19:22 

#13 
New Member
Caroline Vandame
Join Date: Aug 2010
Posts: 26
Rep Power: 7 
Thanks Eric! This is gonna be helpful
could you remind me how to compile it? Caroline 

December 3, 2012, 15:38 

#14 
New Member
KH
Join Date: Oct 2012
Posts: 1
Rep Power: 0 
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 ? Last edited by eagle; December 3, 2012 at 17:56. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
How to get the max value of the whole field  waynezw0618  OpenFOAM Running, Solving & CFD  4  June 17, 2008 05:07 
Which is the fatal error when foamInstallationTest  eddy  OpenFOAM Installation  2  March 30, 2008 21:13 
Installation problem with GCC  Norma McKee (Mckee)  OpenFOAM Installation  10  March 4, 2007 08:09 
a question of open ".cas" and ".dat" files  fanzhong Meng  FLUENT  4  May 15, 2006 11:40 
Fatal error error writing to tmp No space left on device  maka  OpenFOAM Installation  2  April 3, 2006 08:48 