CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Moving source term

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

Reply
 
LinkBack Thread Tools Display Modes
Old   January 5, 2012, 06:12
Default Moving source term
  #1
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Dear foamers,

I would like to try to modify a solver in order to be able to generate a moving source term inside a domain.

In details, I would like to create a transient solver, add vorticity at some point in the domain, and move this point with respect to time on a prescribed trajectory in my domain.

Does anyone have an idea on how to proceed to realize something like that?
vinz is offline   Reply With Quote

Old   January 5, 2012, 07:22
Default
  #2
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Dear foamers,

I would like to try to modify a solver in order to be able to generate a moving source term inside a domain.

In details, I would like to create a transient solver, add vorticity at some point in the domain, and move this point with respect to time on a prescribed trajectory in my domain.

Does anyone have an idea on how to proceed to realize something like that?
One way to accomplish this with a minimum of C++-programming would be swak4Foam. See transient source term for an example on how you have to modify the solver in two places and then use any function for the source term
gschaider is offline   Reply With Quote

Old   January 6, 2012, 05:03
Default
  #3
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Thank you for your reply Bernhard!

As always, you developped the perfect tool

I downloaded swak4Foam, installed it. I also modified pimpleFoam in order to be able to add a momentumSource term. Now I have to find the way to create an expression which will create me some swirl in my domain!
vinz is offline   Reply With Quote

Old   January 6, 2012, 07:13
Default
  #4
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Is there a way to access functions of openfoam like fvm:ddt or fvm::div inside swak4Foam?
vinz is offline   Reply With Quote

Old   January 7, 2012, 19:01
Default
  #5
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Is there a way to access functions of openfoam like fvm:ddt or fvm::div inside swak4Foam?
The functions from the fvm-namespace are used for building a coefficient matrix. So: no. But of course you can use the fvc::div (just say "div" in an expression). ddt is currently not implemented
gschaider is offline   Reply With Quote

Old   January 12, 2012, 03:03
Default
  #6
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Ok I see, thank you for your reply.

I have another question. I saw in one of your presentation about swak4Foam that it is possible to work with cellSet. You give an example which looks like this:
Code:
// Name o f s e t t o o p e r a t e on
name i n B a l l ;
// One o f c l e a r /new/ i n v e r t / add / d e l e t e | s u b s e t / l i s t
a c t i o n new ;

topoSetSources
(
expressionToCell
{
e x p r e s s i o n " alpha1 > 0.5" ;
}
);
Is it possible to recreate this cellSet every time step or every "n" time steps? Let say I want to define a plane in my domain on which I want to define some parameter, and I want to make this plane move with respect to time. Is there a way to do that? How would you proceed?
vinz is offline   Reply With Quote

Old   January 12, 2012, 05:37
Default
  #7
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Ok I see, thank you for your reply.

I have another question. I saw in one of your presentation about swak4Foam that it is possible to work with cellSet. You give an example which looks like this:
Code:
// Name o f s e t t o o p e r a t e on
name i n B a l l ;
// One o f c l e a r /new/ i n v e r t / add / d e l e t e | s u b s e t / l i s t
a c t i o n new ;

topoSetSources
(
expressionToCell
{
e x p r e s s i o n " alpha1 > 0.5" ;
}
);
Is it possible to recreate this cellSet every time step or every "n" time steps? Let say I want to define a plane in my domain on which I want to define some parameter, and I want to make this plane move with respect to time. Is there a way to do that? How would you proceed?
I've played with the idea, but havn't implemented it yet. It shouldn't be too hard to implement. Reason why I havn't done so is a) I never needed it until now b) if another program (functionObject or solver) that isn't aware that a cellSet might change its size uses this then this feature would have the potential to create segmentation faults.

Where would you use that cellSet? If it is the region where you want your source term defined then the (required) "condition" parameter is exactly what you need. If you want to modify an existing field in memory then the manipulateField-functionObject (which has a similar condition) might be your friend
gschaider is offline   Reply With Quote

Old   January 12, 2012, 05:53
Default
  #8
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Ok so that's what I was scared about, I didn't miss anything.

To explain a bit more, here is my current momentumSourceDict:
Code:
UEqn {
variables (
    //trajectory of my center point
    "alpha=0;"
    "Z=0.1*(time()-0.5)+0.1;"
    "Y=0;"
    "X=0;"

    //position with respect to a center point
    "xp=pos().x - X;"
    "yp=pos().y - Y;"
    "r=sqrt(xp*xp+yp*yp);"
    "theta=((xp>0) ? acos(-yp/r) : -acos(-yp/r));"  

    //tangential speed to add as a momentum source term
    "R=0.02;"
    "gama=10000*3.1415*R;" 
    "U_t=((r<R) ? gama*r/(2*3.1415*R*R) : gama/(2*3.1415*r));" 
    "U_t2=((time()>0.51 && pos().z<=Z+0.01 && pos().z>=Z-0.01) ? U_t : 0);"
    "U_x=-U_t2*sin(3.1415/2-theta);"
    "U_y=-U_t2*cos(3.1415/2-theta);"
    "vectU=vector(U_x,U_y,0);"
);

expression   "rho*vectU";

dimensions [1 -2 -2 0 0 0 0];
}
Don't really look at the values, they are only for testing.
With this dictionnary, I managed to add a source term moving in my field. The problem is that I face a constraint which is the size of my mesh. You see I use the condition "pos().z<=Z+0.01 && pos().z>=Z-0.01" to impose my source term. The "+-0.01" is comming from the size of my cells which are 0.01m. If I use a smaller value, the source term is null everywhere in the field. But this way I capture two cells.

The problem is that with this value, my term jumps from one cell to the next one, and I would like to get something smoother.
I will try to make a small movie to show you if needed.

Last edited by vinz; January 12, 2012 at 06:09.
vinz is offline   Reply With Quote

Old   January 12, 2012, 06:59
Default
  #9
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
here is the animation of my moving source term:

vinz is offline   Reply With Quote

Old   January 12, 2012, 19:51
Default
  #10
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Ok so that's what I was scared about, I didn't miss anything.

To explain a bit more, here is my current momentumSourceDict:
Code:
UEqn {
variables (
    //trajectory of my center point
    "alpha=0;"
    "Z=0.1*(time()-0.5)+0.1;"
    "Y=0;"
    "X=0;"

    //position with respect to a center point
    "xp=pos().x - X;"
    "yp=pos().y - Y;"
    "r=sqrt(xp*xp+yp*yp);"
    "theta=((xp>0) ? acos(-yp/r) : -acos(-yp/r));"  

    //tangential speed to add as a momentum source term
    "R=0.02;"
    "gama=10000*3.1415*R;" 
    "U_t=((r<R) ? gama*r/(2*3.1415*R*R) : gama/(2*3.1415*r));" 
    "U_t2=((time()>0.51 && pos().z<=Z+0.01 && pos().z>=Z-0.01) ? U_t : 0);"
    "U_x=-U_t2*sin(3.1415/2-theta);"
    "U_y=-U_t2*cos(3.1415/2-theta);"
    "vectU=vector(U_x,U_y,0);"
);

expression   "rho*vectU";

dimensions [1 -2 -2 0 0 0 0];
}
Don't really look at the values, they are only for testing.
With this dictionnary, I managed to add a source term moving in my field. The problem is that I face a constraint which is the size of my mesh. You see I use the condition "pos().z<=Z+0.01 && pos().z>=Z-0.01" to impose my source term. The "+-0.01" is comming from the size of my cells which are 0.01m. If I use a smaller value, the source term is null everywhere in the field. But this way I capture two cells.

The problem is that with this value, my term jumps from one cell to the next one, and I would like to get something smoother.
I will try to make a small movie to show you if needed.
I understand the problem, but I fail to see how a cellSet would help here.

I think one solution could be to have instead of the step you're using a smoother function that covers more cells but drops of quite quickly to 0 outside of the region of interest (something that looks like a Gauss-function. I'll call it like that)

Steps would be
1. Calculate "Gauss"
2. Volume-integrate it and compare to the desired integral (due to the discretization it won't be exactly the same)
3. Scale whole Gauss so that it fits the integral
4. Use as source term

The scaling takes care of the "now you see it, now you don't"-effect you described above

The continuous nature of the function makes sure that cells do not immediately switch from no-source to full-source (in time) which may lead to weird behaviour at the "boundary" of the source

Finding a "good" function is a bit of a trial and error-thing
gschaider is offline   Reply With Quote

Old   January 13, 2012, 05:25
Default
  #11
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Ok I see.
I worked a bit on the function and it looks much better now:


I use the following functions now:
Code:
    "deltaZ=((Z-pos().z)>0 ? Z-pos().z : pos().z-Z);"
    "U_t2=((time()>0.5 && deltaZ<0.01) ? (1-deltaZ/0.01)*U_t : 0);"
0.01 being in my case the size in Z direction of my mesh.
As a first step to study feasability it is enough for me. I can work on a Gaussian function if needed a bit later.

Is there a way to access cell size to replace this ugly 0.01?

Thanks for your help!
vinz is offline   Reply With Quote

Old   February 3, 2012, 11:15
Default
  #12
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Now that I know how to add a source term in my U equation and how to make it move, I have to know how to impose the correct value.

I have a value of the circulation Gama of the vortex I want to create. From that, I can also compute Vtheta and eventually Vy and Vz if my vortex is in a plane orthogonal to X.

My question is: how can I compute the value of the momentumSource I add in the following way to my equation using the above mentionned variables Gama, Vtheta, Uy, Uz:

Code:
tmp<fvVectorMatrix> UEqn
(
    fvm::ddt(rho, U)
  + fvm::div(phi, U)
  + turbulence->divDevRhoReff(U)
    ==
    momentumSource()

);
or is there a way to directly add Uy and Uz induced by my vortex to my flow field?
vinz is offline   Reply With Quote

Old   February 3, 2012, 13:53
Default
  #13
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Now that I know how to add a source term in my U equation and how to make it move, I have to know how to impose the correct value.

I have a value of the circulation Gama of the vortex I want to create. From that, I can also compute Vtheta and eventually Vy and Vz if my vortex is in a plane orthogonal to X.

My question is: how can I compute the value of the momentumSource I add in the following way to my equation using the above mentionned variables Gama, Vtheta, Uy, Uz:

Code:
tmp<fvVectorMatrix> UEqn
(
    fvm::ddt(rho, U)
  + fvm::div(phi, U)
  + turbulence->divDevRhoReff(U)
    ==
    momentumSource()

);
or is there a way to directly add Uy and Uz induced by my vortex to my flow field?
There is an object forceEquation in swak that allows you to fix the values in a certain region (see the interFoamWithFixed-example). Of course doing that has to be done with care as it may lead to situations were conservation is violated.

The problem I see for your application is that you want one component of U to be "free". That is currently possible (I'm not totally sure if it is even feasible without a lot of dirty workarounds)
gschaider is offline   Reply With Quote

Old   February 3, 2012, 14:12
Default
  #14
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Yes the basic idea is just that I want to add some velocity to my velocity components in one plane, here Vy and Vz. The other component is usually 0 unless I impose a inlet value in my domain.

I will have a look at the forceEquation.

The momentumSource as I have added it in the equation above is supposed to be a massic force, right? I am currently adding some kind of vector which qualitatively makes the flow spin as want it to spin. But I have not found the physical equation allowing me to define mathematically the value of this source term based on the circulation for instance.
vinz is offline   Reply With Quote

Old   February 15, 2012, 09:14
Default
  #15
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Dear Bernhard,

I have one more question. I noticed there is parser for the function "time()" Is there a way to access the current time step value into swak4Foam?

If it is not the case, what is/are the file(s) in the library I should try to modify to create a parser for a function which could be named "timeStep()" for example?

thanks in advance for your help.
vinz is offline   Reply With Quote

Old   February 15, 2012, 09:28
Default
  #16
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Ok sorry for this stupid question, it looks like deltaT() is the way to get it.
vinz is offline   Reply With Quote

Old   February 15, 2012, 09:48
Default
  #17
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Ok sorry for this stupid question, it looks like deltaT() is the way to get it.
That's right
gschaider is offline   Reply With Quote

Reply

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
transient source term strakakl OpenFOAM 38 November 19, 2013 02:18
Problem of SOURCE term gradient in UDS wind Fluent UDF and Scheme Programming 5 June 21, 2013 05:39
How to write source term into scalar Fiel JimKnopf OpenFOAM Programming & Development 0 March 23, 2011 06:59
Compiling gmshFoam with OpenFOAM-1.5 BlGene Open Source Meshers: Gmsh, Netgen, CGNS, ... 10 August 6, 2009 04:26
UDF Source Term Units? Brian FLUENT 1 October 24, 2005 09:15


All times are GMT -4. The time now is 16:44.