
[Sponsors] 
January 5, 2012, 06:12 
Moving source term

#1 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
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? 

January 5, 2012, 07:22 

#2  
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,184
Rep Power: 46 
Quote:


January 6, 2012, 05:03 

#3 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
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! 

January 6, 2012, 07:13 

#4 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
Is there a way to access functions of openfoam like fvm:ddt or fvm::div inside swak4Foam?


January 7, 2012, 19:01 

#5 
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,184
Rep Power: 46 
The functions from the fvmnamespace 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


January 12, 2012, 03:03 

#6 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
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" ; } ); 

January 12, 2012, 05:37 

#7  
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,184
Rep Power: 46 
Quote:
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 manipulateFieldfunctionObject (which has a similar condition) might be your friend 

January 12, 2012, 05:53 

#8 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
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>=Z0.01) ? U_t : 0);" "U_x=U_t2*sin(3.1415/2theta);" "U_y=U_t2*cos(3.1415/2theta);" "vectU=vector(U_x,U_y,0);" ); expression "rho*vectU"; dimensions [1 2 2 0 0 0 0]; } 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>=Z0.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. 

January 12, 2012, 06:59 

#9 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
here is the animation of my moving source term:


January 12, 2012, 19:51 

#10  
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,184
Rep Power: 46 
Quote:
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 Gaussfunction. I'll call it like that) Steps would be 1. Calculate "Gauss" 2. Volumeintegrate 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 nosource to fullsource (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 errorthing 

January 13, 2012, 05:25 

#11 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
Ok I see.
I worked a bit on the function and it looks much better now: I use the following functions now: Code:
"deltaZ=((Zpos().z)>0 ? Zpos().z : pos().zZ);" "U_t2=((time()>0.5 && deltaZ<0.01) ? (1deltaZ/0.01)*U_t : 0);" 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! 

February 3, 2012, 11:15 

#12 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
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() ); 

February 3, 2012, 13:53 

#13  
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,184
Rep Power: 46 
Quote:
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) 

February 3, 2012, 14:12 

#14 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
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. 

February 15, 2012, 09:14 

#15 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
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. 

February 15, 2012, 09:28 

#16 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 283
Rep Power: 14 
Ok sorry for this stupid question, it looks like deltaT() is the way to get it.


February 15, 2012, 09:48 

#17 
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,184
Rep Power: 46 

October 1, 2019, 08:20 
Movement of boiling plane with CONDITION

#18  
Senior Member

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:
Quote:
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 

October 2, 2019, 17:01 

#19  
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,184
Rep Power: 46 
Quote:
If the solver you're using supports fvOptions then you can add a sourceterm with the scalarSwakExplicitSourcefvOption
__________________
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 

October 3, 2019, 06:00 

#20 
Senior Member

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. 

Thread Tools  Search this Thread 
Display Modes  


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 06:39 
How to write source term into scalar Fiel  JimKnopf  OpenFOAM Programming & Development  0  March 23, 2011 06:59 
[Gmsh] Compiling gmshFoam with OpenFOAM1.5  BlGene  OpenFOAM Meshing & Mesh Conversion  10  August 6, 2009 05:26 
UDF Source Term Units?  Brian  FLUENT  1  October 24, 2005 10:15 