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

Moving source term

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

Like Tree6Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 5, 2012, 06:12
Default Moving source term
  #1
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 18
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?
missios likes this.
vinz is offline   Reply With Quote

Old   January 5, 2012, 07:22
Default
  #2
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,225
Rep Power: 51
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 http://www.cfd-online.com/Forums/ope...tml#post329700 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: 283
Rep Power: 18
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: 283
Rep Power: 18
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: 4,225
Rep Power: 51
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: 283
Rep Power: 18
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?
Kummi likes this.
vinz is offline   Reply With Quote

Old   January 12, 2012, 05:37
Default
  #7
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,225
Rep Power: 51
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: 283
Rep Power: 18
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.
Kummi likes this.

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: 283
Rep Power: 18
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: 4,225
Rep Power: 51
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: 283
Rep Power: 18
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: 283
Rep Power: 18
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: 4,225
Rep Power: 51
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: 283
Rep Power: 18
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: 283
Rep Power: 18
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: 283
Rep Power: 18
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: 4,225
Rep Power: 51
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

Old   October 1, 2019, 08:20
Question Movement of boiling plane with CONDITION
  #18
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 346
Rep Power: 11
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Hello Foamers,
I know this is old thread, however I feel that I can get my answer here based on the above discussion. Kindly share the ideas please.

This work is a part of main problem. The geometry is a simple rectangular closed block full of liquid (moisutre). The heat (temperature) is supplied from one wall end and at the saturation temperature of 100deg, evaporation occurs (@T=100deg, alpha (liquid) = 0].

As vinz mentioned above,
Quote:
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.
This statement highly relates my problem. In my work, the boiling plane is defined with the position (where CONDITION: @T=100, alpha =0) followed by the calculation of mass and heat fluxes, where the mass flux is calculated by the plane movement with respect to time [ATTACHMENT 1 gives clear overview]. In addition, depending upon boiling plane distance (xb) - the calculations are made for high (-) and low temperature sides (+).
Quote:
ENERGY EQU: rho*Cp*(dT/dt) = del/delx[K*dT/dx ] + rho*dQ/dT + r*cp*dT/dx (W/m3)
Unsteady conduction = Diffusion term + SOURCE TERM 1 + SOURCE TERM 2
Mass flux is quantified as rate of evaporation (r). The calculated (r) is multiplied with Cp and dT/dx to give source term.

Source term explained in main ENERGY equ.... r*cp*dt/dx ==> here (r) evaporation rate possess the dxb/dt (rate of progression of boiling plane) --> so this kind of problem can be moving source term problem? Correct me if Im wrong please.

Can the function objects help in achieving it (or) something like swak4FOAM ??
Kindly guide me Foamers. Thank you
Attached Images
File Type: jpg ATTACHMENT_1_Boiling_ .jpg (105.7 KB, 21 views)
Kummi is offline   Reply With Quote

Old   October 2, 2019, 17:01
Default
  #19
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,225
Rep Power: 51
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by Kummi View Post
Hello Foamers,
I know this is old thread, however I feel that I can get my answer here based on the above discussion. Kindly share the ideas please.

This work is a part of main problem. The geometry is a simple rectangular closed block full of liquid (moisutre). The heat (temperature) is supplied from one wall end and at the saturation temperature of 100deg, evaporation occurs (@T=100deg, alpha (liquid) = 0].

As vinz mentioned above, This statement highly relates my problem. In my work, the boiling plane is defined with the position (where CONDITION: @T=100, alpha =0) followed by the calculation of mass and heat fluxes, where the mass flux is calculated by the plane movement with respect to time [ATTACHMENT 1 gives clear overview]. In addition, depending upon boiling plane distance (xb) - the calculations are made for high (-) and low temperature sides (+).Mass flux is quantified as rate of evaporation (r). The calculated (r) is multiplied with Cp and dT/dx to give source term.

Source term explained in main ENERGY equ.... r*cp*dt/dx ==> here (r) evaporation rate possess the dxb/dt (rate of progression of boiling plane) --> so this kind of problem can be moving source term problem? Correct me if Im wrong please.

Can the function objects help in achieving it (or) something like swak4FOAM ??
Kindly guide me Foamers. Thank you

If the solver you're using supports fvOptions then you can add a source-term with the scalarSwakExplicitSource-fvOption
Kummi likes this.
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request
gschaider is offline   Reply With Quote

Old   October 3, 2019, 06:00
Default
  #20
Senior Member
 
Kumaresh
Join Date: Oct 2016
Posts: 346
Rep Power: 11
Kummi is on a distinguished road
Send a message via Yahoo to Kummi
Thank you for your response Bernhard Gschaider! fvOptions is not supported in my work.

Actually My problem is coal pyrolysis using region Model in OpenFOAM. The code is constructed in a region for pyrolysis of DRY COAL~ no issues with it.. As a next step, I want to construct code for WET COAL PYROLYSIS. In this process of wet coal pyrolysis, moisture evaporates at 100deg before pyrolysis occurs (at 400deg). Inside pyrolysis code - I should introduce this Moisture evaporation code.
.
So my discussion here is not about pyrolysis. Discussion is about "Evaporation of moisture".

In my work, the Moisture evaporation is based on some simple technique.

(1) by defining the position of the plane
(2) based on it, mass and heat fluxes are calculated

My question is that how the code can be constructed defining position of plane and make this plane move with respect to time.
Kummi is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem of SOURCE term gradient in UDS wind Fluent UDF and Scheme Programming 6 December 1, 2022 15:21
transient source term strakakl OpenFOAM 38 November 19, 2013 02:18
How to write source term into scalar Fiel JimKnopf OpenFOAM Programming & Development 0 March 23, 2011 06:59
[Gmsh] Compiling gmshFoam with OpenFOAM-1.5 BlGene OpenFOAM Meshing & Mesh Conversion 10 August 6, 2009 05:26
UDF Source Term Units? Brian FLUENT 1 October 24, 2005 10:15


All times are GMT -4. The time now is 01:30.