[vinz]

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?

[gschaider]

Originally Posted by vinz

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

[vinz]

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 there a way to access functions of openfoam like fvm:ddt or fvm::div inside swak4Foam?

[gschaider]

Originally Posted by vinz

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

[vinz]

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:

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

[gschaider]

Originally Posted by vinz

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:

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

[vinz]

Ok so that's what I was scared about, I didn't miss anything.

To explain a bit more, here is my current momentumSourceDict:

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.

[vinz]

here is the animation of my moving source term:

[gschaider]

Originally Posted by vinz

Ok so that's what I was scared about, I didn't miss anything.

To explain a bit more, here is my current momentumSourceDict:

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

[vinz]

Ok I see.
I worked a bit on the function and it looks much better now:


I use the following functions now:

    "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]

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:

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?

[gschaider]

Originally Posted by vinz

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:

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)

[vinz]

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]

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]

Ok sorry for this stupid question, it looks like deltaT() is the way to get it.

[gschaider]

Originally Posted by vinz

Ok sorry for this stupid question, it looks like deltaT() is the way to get it.

That's right

[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,

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 (+).

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

[gschaider]

Originally Posted by 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, 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]

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.