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

Initialization of U by code

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

Like Tree1Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   August 13, 2009, 11:20
Default Initialization of U by code
  #1
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Hi all!. I'm using scalarTransportFoam to solve de transport equation. It's running, but i have a problem. I need to set the initial field of U. I know that there are several utilities to do similar things like funkySetFields, but in this case, i think that i need to write some code.
The mesh is simply a tube, and i want to set up a rigid rotation as a velocity field. I wrote a code using the vector operations provided by OpenFOAM but it crashes with i try to compile it.
Then i tryed putting a portion of H. Jasak code in a .H file, as follows:

USetUp_Jasak.H

// Do cells
const volVectorField& centres = mesh.C();
point origin(1, 1, 0.05);
vector axis(0, 0, -1);
U = axis ˆ (centres.internalField() - origin);
U.write();


This code have no sense to my purposes, but it's only for trying to compile.

Then i added this line in scalarTransportFoam.C file:

int main(int argc, char *argv[])
{

# include "setRootCase.H"

# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "USetUp_Jasak.H" <<<<<<<<<<----------------------

...

It implies that, next U was created and initialized by reading U in /0 directory, the code overwrite the values with the new ones.

When i try to compile, i obtain this:

santiago@elduque:~/OpenFOAM/santiago-1.5/applications/myScalarTransportFoam$ wmake
SOURCE=myScalarTransportFoam.C ; g++ -m32 -Dlinux -DDP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -O3 -DNoRepository -ftemplate-depth-40 -I/home/santiago/OpenFOAM/OpenFOAM-1.5/src/finiteVolume/lnInclude -IlnInclude -I. -I/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude -I/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OSspecific/Unix/lnInclude -fPIC -pthread -c $SOURCE -o Make/linuxGccDPOpt/myScalarTransportFoam.o
In file included from myScalarTransportFoam.C:47:
USetUp_Jasak.H: In function ‘int main(int, char**)’:
USetUp_Jasak.H:5: error: no match for ‘operator=’ in ‘U = Foam:perator^(const Foam::VectorSpace<Form, Cmpt, nCmpt>&, const Foam::tmp<Foam::Field<Type> >&) [with Form = Foam::Vector<double>, Cmpt = double, int nCmpt = 3, Type = Foam::Vector<double>](((const Foam::tmp<Foam::Field<Foam::Vector<double> > >&)((const Foam::tmp<Foam::Field<Foam::Vector<double> > >*)(& Foam:perator-(const Foam::UList<T>&, const Foam::VectorSpace<Form, Cmpt, nCmpt>&) [with Type = Foam::Vector<double>, Form = Foam::Vector<double>, Cmpt = double, int nCmpt = 3](((const Foam::VectorSpace<Foam::Vector<double>, double, 3>&)((const Foam::VectorSpace<Foam::Vector<double>, double, 3>*)(& origin.Foam::Vector<double>::<anonymous>))))))))’
/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/GeometricField.C:1032: note: candidates are: void Foam::GeometricField<Type, PatchField, GeoMesh>:perator=(const Foam::GeometricField<Type, PatchField, GeoMesh>&) [with Type = Foam::Vector<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/GeometricField.C:1057: note: void Foam::GeometricField<Type, PatchField, GeoMesh>:perator=(const Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh> >&) [with Type = Foam::Vector<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/GeometricField.C:1093: note: void Foam::GeometricField<Type, PatchField, GeoMesh>:perator=(const Foam::dimensioned<Type>&) [with Type = Foam::Vector<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
make: *** [Make/linuxGccDPOpt/myScalarTransportFoam.o] Error 1

Any clues?, thanks in advance. Santiago.
santiagomarquezd is offline   Reply With Quote

Old   August 18, 2009, 05:15
Default
  #2
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
In the line:

U = axis ˆ (centres.internalField() - origin);

You have U as a volVectorField, axis as a vector, the internalField of centres and a single point, origin.

Have you tried:

U.internalField() = axis ˆ (centres.internalField() - origin);

Always think about what the fields mean.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   August 18, 2009, 11:20
Default
  #3
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Luarence, your insight was very useful, thank you. I've tryed that with other code but this code was buggy, then i didn't realized that i have to access the U data by this method.
The code of this post was a test code, i wrote the final version and it compiled, but in runtime i have some problems. This is the code:

// Velocity field initialization

// Angular velocity
scalar omega=10;

// Cell centroid coordinates
const volVectorField& centres = mesh.C();

// Coordinate system origin
point origin(0, 0, 0);

// Axis of rotation
vector axis(0, 0, 1);

// Force arm calculation
// Memory allocation
volVectorField palanca
(
IOobject
(
"pal",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh
);


palanca.internalField() = centres - ((centres.internalField()-origin) & axis)*axis;

// Velocity field calculation
U.internalField()=omega*mag(palanca.internalField( ))*axis ^ palanca.internalField();


in runtime i obtain this:

/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Exec : myScalarTransportFoam
Date : Aug 18 2009
Time : 11:52:36
Host : elduque
PID : 6414
Case : /home/santiago/OpenFOAM/santiago-1.5/run/tuboMyTransport
nProcs : 1

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading field T

Reading field U

Reading transportProperties

Reading diffusivity D

Reading/calculating face flux field phi



NO_READ specified for read-constructor of object pal of class IOobject#0 Foam::error:rintStack(Foam::Ostream&) in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Foam::regIOobject::readStream() in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#3 Foam::regIOobject::readStream(Foam::word const&) in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#4 Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>::GeometricField(Foam::IOobject const&, Foam::fvMesh const&) in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#5 main in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#6 __libc_start_main in "/lib/i686/cmov/libc.so.6"
#7 Foam::regIOobject::readIfModified() in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"


From function regIOobject::readStream(const word&)
in file db/regIOobject/regIOobjectRead.C at line 51.

FOAM aborting

Aborted

Really the only reason because i define the volVectorField palanca is to use it as a temporal variable. Because that i set NO_READ and NO_WRITE. It seems that this settings are incorrect. How I can define a temporal vector field to store the intermediate values (with the smallest useful structure)?. If volVectorField is the unique structure that i can use, how i can set it properly?.

Thanks.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   August 18, 2009, 11:38
Default
  #4
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
You could create a temporary object:

tmp<volVectorField> variableName;

In your case your temporary variable is just a field of vectors, so you could just have a tmp<vectorField>. The vectorField is like a volVectorField but without the boundary conditions.
So something like:

tmp<vectorField> palanca = centres - ((centres.internalField()-origin) & axis)*axis;
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   August 18, 2009, 18:36
Default
  #5
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Thanks Laurence!, slowly I'm understanding the structure of FOAM classes. tmp<vectorField> did the trick.

Now i could impose a rigid rotation velocity field with this lines:

tmp<vectorField> palanca;

palanca = centres.internalField() - ((centres.internalField()-origin) & axis)*axis;

// Velocity field calculation
U.internalField()=omega*axis ^ palanca;
U.write();

It compiled and ran properly, i verified the field using Paraview. Now i want to scale the field by the palanca's magnitude (the norm of arm of force vector), then i tryed the following:

U.internalField()=mag(palanca)*omega*axis ^ palanca;

It's good at compilation but don't run, this is the output:

/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Exec : myScalarTransportFoam
Date : Aug 18 2009
Time : 19:38:28
Host : elduque
PID : 9150
Case : /home/santiago/OpenFOAM/santiago-1.5/run/tuboMyTransport
nProcs : 1

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading field T

Reading field U

Reading transportProperties

Reading diffusivity D

Reading/calculating face flux field phi



temporary deallocated#0 Foam::error::printStack(Foam::Ostream&) in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#3 Foam::tmp<Foam::Field<Foam::crossProduct<Foam::Vec tor<double>, Foam::Vector<double> >::type> > Foam::operator^<Foam::Vector<double>, Foam::Vector<double> >(Foam::tmp<Foam::Field<Foam::Vector<double> > > const&, Foam::tmp<Foam::Field<Foam::Vector<double> > > const&) in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#4 main in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#5 __libc_start_main in "/lib/i686/cmov/libc.so.6"
#6 Foam::regIOobject::readIfModified() in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"


From function const T& tmp<T>::operator()() const
in file /home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/tmpI.H at line 190.

FOAM aborting

Aborted

I think something is wrong with the magnitude operator. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   August 18, 2009, 20:19
Default
  #6
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
Maybe forget about the use of tmp, so just define palanca as a vectorField:

Code:
vectorField palanca = centres.internalField() - ((centres.internalField()-origin) & axis)*axis;

U.internalField()=mag(palanca)*omega*axis ^ palanca;
Someone can correct me, but I don't think you can take the mag of a tmp. You could try:

Code:
U.internalField()=mag(palanca())*omega*axis ^ palanca;
But I think for your purposes the first option would make life easier. Unless you're worried about memory usage.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   August 19, 2009, 11:45
Default
  #7
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Laurence, thank you again, your tips were very helpful, I'm imposing the velocity field as I wanted. Now, I'm not completely sure about freeing the temporal memory used, i tryed to invoke the vectorField destructor without success. What is the way in FOAM?. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   August 19, 2009, 14:32
Default
  #8
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
I did a quick search on the forums about tmp, you should be able to find out more about using it. tmp variables can be cleared by using .clear(). You can avoid dereferencing them by putting brackets after the variable when you use it.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   August 21, 2009, 20:25
Default
  #9
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Thanks Laurance for all your support. I'm dealing now with the appropriate algebra to set the vector field, i want to compare the effects of high rotation velocities in advecition-diffusion solvers, specifically between Fluent and OpenFOAM. This is a work in colaboration with a team of chemical engineers with the aim of model polimerization reactors. We aren't completely sure about the influence of false diffusion in RTD determination (using the tracer method). Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   August 22, 2009, 04:59
Default
  #10
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
Sounds interesting, I do population balance modelling and turbulent, reactive, multiphase flow modelling. Numerical diffusion is something I've had to be careful about too.
Do you use population balances when modelling the reaction? Or are you just trying to look at the flow?
We've done a little bit of basic RTD analysis using Lagrangian particle tracking for a stirred tank in both experiments and simulations. The comparison wasn't great, and I think the main problems arose due to:

1) boundary conditions
2) inadequate turbulence modelling

I would "imagine" that numerical diffusion would have less effect on the mean residence time but a greater effect on the shape/width of the distribution?

Keep everyone posted on your progress.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   August 22, 2009, 10:32
Default
  #11
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Laurence, actually we aren't do chemical simulations really. We use a polimer-equivalent materials, most of them with a very high dynamic viscosity (like 1 Pa-seg). All of this lead us to use laminar-regime solvers. Once the velocity field is obtained we run a tracer to obtain the RTD.
In industry, because of the large time of simulation needed, all the people uses very long time steps, usually the RTD obtained looks like a Poisson distribution a lot.
My co-advisor suggested that behavior is due mainly to numerical diffusion in stead of real diffusion, in fact we use small diffusion factors in the solution of adv-diff equation. I'm suspecting now that this sounds more accurate in FEM solvers, because of the streamline diffusion methods (SUPG).
I'm finishing the coding now, when i have the results i'll post them here, so as to retain the topic.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   August 22, 2009, 21:19
Default
  #12
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Hi, Laurence & all the community, I've implemented the initialization for the U field. My new question is how change only one component of a field. Digging in the marvelous forums i found the following method, it compiled without warnings,

// Parabolic distribution in z component
scalarField zComp=Vmax*(1-sqr(mag(palanca)/tubeR));

// Velocity field calculation
U.internalField()=omega*(axis^palanca)+vector(0,0,1)*zComp;

basically it consists in using a mask, but... Is there a more "pure" way to do this, for example by a class' method?

Thanks in advance.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   August 24, 2009, 04:37
Default
  #13
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
You may use Doxygen to find this out. If you search for scalarField in it, you'll find that it's a type definition: Field<type>. The following link gives the class reference for Field:

http://foam.sourceforge.net/doc/Doxy..._1_1Field.html

In the public member functions you'll find component. So to access one component of zComp would be (0=x axis, 1 = yaxis etc.):

Code:
zComp.component(2);
Doxygen is very useful!
dawnrain likes this.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   August 24, 2009, 08:55
Default
  #14
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Yes, it's true, Doxygen is a powerful tool. I forgot to tell it, but I was doing research in Doxygen too. I've found the component method, but probably i used it in a wrong way, so your message push me to try again and now it works.
Actually I think that my previous message was bad written, i wanted to overwrite the z component of the U field (this was really the main topic). Anyway, the final code is the following:

// Velocity field calculation
//Solid rotation
U.internalField()=omega*(axis^palanca);

// Parabolic distribution in z component
U.internalField().component(2)=Vmax*(1-sqr(mag(palanca)/tubeR));

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   August 31, 2009, 12:47
Default
  #15
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Hi folks, as a continuation of my last post I'm copying all the code:

// Velocity field initialization

// Angular velocity
scalar omega=100;

// Maximum velocity en parabolic profile
scalar Vmax=2.0;

// Tube radius
scalar tubeR=0.1;

// Cell centroid coordinates
const volVectorField& centres = mesh.C();

// Coordinate system origin
point origin(0, 0, 0);

// Axis of rotation
vector axis(0, 0, 1);

// Force arm calculation
vectorField palanca;

palanca = centres.internalField() - ((centres.internalField()-origin) & axis)*axis;

// Velocity field calculation
//U.internalField()=omega*(axis^palanca);

// Parabolic distribution in z component
U.internalField().component(2)=Vmax*(1-sqr(mag(palanca)/tubeR));

U.write();

It's OK in compilation & running time but the results are bad. I obtain a zero field in the z component of U in stead of the parabolic field that i wanted. palanca has the distance of each cell center to the axis of rotation, then the equation for component 2 may give me a laminar parabolic profile in all the tube.
I don't know if the problem is in the mag or sqr methods or something else. All answers are welcome. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   August 31, 2009, 13:10
Default
  #16
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Some other important data:

1. I initialized U field to (0, 0, 0).
2. If we try to the following line:

U.internalField().component(2)=mag(palanca);

U field continue having (0, 0, 0)...

I guess that the component method only give us a copy of the data not a reference to overwrite the original.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   August 31, 2009, 17:59
Default
  #17
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
Double check your axis vector. It should be the axis that your flow rotates around. Let's say that 'y' is up, 'x' is right and 'z' is out of the page. You want flow in the 'z' direction, so your axis vector should be your x or y axis. i.e. (0 1 0) or (1 0 0).
Does that make sense?

Edit: Think I misunderstood. Doesn't look like you can overwrite component. Try 'replace' (see Doxygen for GeometriField).
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   August 31, 2009, 19:01
Default
  #18
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Maybe my no so pretty drawing can help us to understand. I want to impose a rotational field with a non-zero tangential velocity.

// Velocity field calculation
//U.internalField()=omega*(axis^palanca);

Due I impose (0,0,0) as the initial condition when I set the rotational velocity the axial velocity remains zero. Then I want to use:

// Parabolic distribution in z component
U.internalField().component(2)=Vmax*(1-sqr(mag(palanca)/tubeR));

to set the axial velocity. Component(2) implies setting the z component of the field, i.e. the axial one.
Actually the problem is that I've checked out with Paraview that this last line don't change the field. I guess that the component method only give us a copy of the component field not a reference, then the assignation is useless, but, I'm only guessing.
When I use this:

// Parabolic distribution in z component
scalarField zComp=Vmax*(1-sqr(mag(palanca)/tubeR));

// Velocity field calculation
U.internalField()=omega*(axis^palanca)+vector(0,0, 1)*zComp;things goes very well, trying to use component() was only to be more "pure" but it seems not to work. Will be another method or correct form of using component() to code this initialization a bit more elegant (and also to learn something more about OpenFOAM coding)?
Regards.
Attached Images
File Type: jpg geometry.jpg (34.7 KB, 19 views)
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   September 1, 2009, 05:01
Default
  #19
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
Yes, I see, try;

scalarField zComp=Vmax*(1-sqr(mag(palanca)/tubeR))
U.internalField().replace(Vector::Z,zComp);

I got replace from the Doxygen entry on GeometricField.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   September 24, 2009, 16:49
Default
  #20
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
Laurence thanks for your answer, the code is running properly, now a detail the line:

Code:
 U.internalField().replace(Vector::Z,zComp);
must be

Code:
 U.internalField().replace(vector::Z,zComp);
or

Code:
 U.internalField().replace(Vector<scalar>::Z,zComp);
in order to run. The first line is calling directly the class Vector without giving all the information to the template, in this case the type of Vector that you want.
Finished this matter, and have another question, this one stopped my from your last post. When I run the case from a ./0/U with all the internalField set in 0 I can't obtain transport within the tube. But once I've ran the case and delete all the timesteps but the first one (0) and then run again things go well. I've checked that OpenFOAM writes the calculated (with our routine) U field in ./0/U field. I guess that FOAM uses the read data in stead of calculated or something like this, then things only run the second time.
Actually I'm not completely sure about what is the function of .write method (use of write is like in my first post, next of the initialization). Firstly I thought that it was to overwrite the values effectively, but now I think that this is the step when U field is being written in ./0 directory. But it remains in memory next to the .write calling?

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Reply

Tags
set up fields, vector manipulation

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
From 2D compressible code to 3D code David Liu Main CFD Forum 22 June 26, 2012 17:59
code initialization - exclude pressure edna CD-adapco 0 January 29, 2008 11:56
CFD code structure (F90) ma Main CFD Forum 4 January 10, 2005 21:47
Design Integration with CFD? John C. Chien Main CFD Forum 19 May 17, 2001 15:56
State of the art in CFD technology Juan Carlos GARCIA SALAS Main CFD Forum 39 November 1, 1999 15:34


All times are GMT -4. The time now is 09:34.