CFD Online Discussion Forums - Initialization of U by code

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Initialization of U by code (https://www.cfd-online.com/Forums/openfoam-programming-development/67386-initialization-u-code.html)

Sorry for the late reply, I've been really busy!
I'm not sure why you're getting the behaviour that you observe. Possibly U has gone out of scope or something? The write() function just writes the current U field to the 0/U file.
I'd probably need to see it for myself, feel free to post your files in a tarball and I can take a quick look.

Thanks Laurence for your attention even with all your work. Well, this is the code:

Code:

// Velocity field initialization



// Angular velocity

scalar omega=2000;



// Maximum velocity en parabolic profile

scalar Vmax=2.0;



// Tube radius

scalar tubeR=0.05;



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





// Parabolic velocity field calculation

scalarField zComp=Vmax*(1-sqr(mag(palanca)/tubeR));



// Option with rotation and the parabolic

// Rotation

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

// Parabolic in z direction

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



// Write to /0/U

U.write();

The sequence of includes is that follows:

Code:

#   include "createTime.H"

#   include "createMesh.H"

#   include "createFields.H"



#   include "USetUp.H"

(mine is the last one, added in scalarTransportFoam.C)

As you can see, after rewriting in memory the U vector field I call the .write() method to write the field to the hard disk. The problem is that if you run the code once, the initializated U field doesn't remain in memory, i.e. the advection equation is solved with another field, maybe zero, I don't
know. But, in the other hand, if you run the code some steps, then stop it and run it again with the /0/U recently written things go well.

That's all. Regards.

This worried me slightly, so I took a look, and I think I've figured it out.

The scalar transport equation takes phi as the coefficient in fvm::div, not U. So you have overwritten U, but now you have to recreate the flux.

So just put

Code:

#include "createPhi.H"

below

Code:

#include "USetUp.H"

Hope that's the problem.

Laurence, I've been working in another things last months but I've returned to this problem. Your suggestion was right, FOAM was assembling the flux with the original U not with the one calculated by my code. I had to put the

Code:

#include "createPhi.H"

below

Code:

#include "USetUp.H"

Nevertheless it is necessary to comment the original invocation of createPhi.H file in createFields.H because the

Code:

#ifndef createPhi_H

clause within createPhi.H. As you know, this line impedes that definition can be repeated.

Regards.

Dear Santiago,

I am struggling with getting a parabolic internal field for pipe flow in openfoam and I've found your post very useful.
Could you please give me some hints on the line "U.internalField()=omega*(axis^palanca);" ?
Another silly question, where shall I put the code, in 0/U or somewhere else? Do I need to derive my own solver for it?
I really appreciate your reply. Thank you very much.

Best regards,
Tony

Greetings to all!

FYI for future readers: Tony's question has been in essence answered here: http://www.cfd-online.com/Forums/ope...tml#post468380

@Tony: Regarding your question about Santiago's implementation: from what I can figure out from the first few posts, it seems that he created a variant of the solver scalarTransportFoam and added a new header file that includes the code discussed in this thread.

It takes a while to explain it in more detail, so I suggest that you have a look at the following tutorials, in order to get a bit more perspective into the details:

In fact, you'll see that those two tutorials are similar to each other, since the first one implements the explanation given in the second one.

Best regards,
Bruno

Hi, Bruno

Thank you very much for the information. I will have a look at that and get back to you if I have further questions.

Best regards,
Tony