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

#1 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 281
Rep Power: 10 
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: 3,920
Rep Power: 41 
Quote:


January 6, 2012, 05:03 

#3 
Senior Member
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 281
Rep Power: 10 
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: 281
Rep Power: 10 
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: 3,920
Rep Power: 41 
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: 281
Rep Power: 10 
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: 3,920
Rep Power: 41 
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: 281
Rep Power: 10 
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: 281
Rep Power: 10 
here is the animation of my moving source term:


January 12, 2012, 19:51 

#10  
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,920
Rep Power: 41 
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: 281
Rep Power: 10 
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: 281
Rep Power: 10 
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: 3,920
Rep Power: 41 
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: 281
Rep Power: 10 
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: 281
Rep Power: 10 
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: 281
Rep Power: 10 
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: 3,920
Rep Power: 41 

Thread Tools  
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 05:39 
How to write source term into scalar Fiel  JimKnopf  OpenFOAM Programming & Development  0  March 23, 2011 06:59 
Compiling gmshFoam with OpenFOAM1.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 