Momentum source
 

Hello everyone!

I created this thread in the common forum aswell, didn't realise there was a special one for UDFs.

I am trying to create an increase of velocity in a flow over a wing in order to simulate a gust wind. The idea is to apply a momentum source in a region upstream of the wing and let the perturbation flow downstream through convection. I found this

https://www.sharcnet.ca/Software/Flu...ug/node309.htm

When I use the UDF below I get strange changes in the velocity field where smaller cells gets a much higher velocity than larger. I suppose it sort of makes sense the source term divides by the volume and I fear I have misunderstood something.

If anyone have experience with this kind of implementation, or has another suggestion of achieving a similar result, it would be of great help. Below you find the udf and a picture of the flowfield.

#include "udf.h"

real A[ND_ND]; /* Declared as global variables*/
real area;

DEFINE_ON_DEMAND(face_area)
{
Domain *d = Get_Domain(1);
Thread *ft;
face_t f;

thread_loop_f(ft,d) /* Finding the area of the cell face with a define on demand macro. Found this in another thread here on CFD-online */
{
begin_f_loop(f,ft)
{
real NV_VEC(A);
F_AREA(A,f,ft);
area = NV_MAG(A);
}
end_f_loop(f,ft)
}

}

DEFINE_SOURCE(source_still_in_sharp_smooth, c,t,dS,eqn)
{

real x[ND_ND];
C_CENTROID(x,c,t)
real x_0 = x[0];
real pi = 3.14159265359;

real Wg = 10; /* amplitude gust */

real start_g = -5; /* start gust source region */
real length_g = 1; /* length gust source region */
real end_g = (start_g + length_g); /* end gust source region */

real real_time = CURRENT_TIME;
real start_time_gust = 0.01; /* start time gust */
real end_time_gust = 0.51; /* end time gust */
real gust_time = (real_time - start_time_gust);
real gust;


if(0 <= gust_time && gust_time <= end_time_gust)
{
if(start_g <= x_0 && x_0 <= end_g)
{
gust = (C_U(c,t) + 5) * C_R(c,t) * area / C_VOLUME(c,t); /* flux through cell, 5 [m/s] is the perturbation */
}
else
{
gust = 0.0;
}
}
else
{
gust = 0.0;
}


real source = gust; /* whetever your source is */
dS[eqn] = C_R(c,t) * area / C_VOLUME(c,t); /* derivative source term */

return source;

}

your problem is the area variable... the initial loop that finds the area, loops over all faces and keeps rewriting over the variable. In the end it is saving the area of the last face into the area variable. Then you area using that same number for all cells. You have to use the area of each cell individually.

Thank you very much for your answer!

I suspected it was something off with using the macros in that way. Do you perhaps have a recommendation of how I could include finding the area in the define_source macro?

I tried to run it with just a static source, divided by the volume for a certain domain. For example:

source = 100/C_VOLUME(c,t);

but the same error with the higher velocity at smaller cells appear. Maybe I am interpreting the momentum source terms unit wrong? As I understood it it is N/m3, meaning I should divide by the volume for each cell?

Regards!

No, you should not divide by the cell volume anymore. Fluent will assign the cell volume-weighted source value when going through the cell loop. What Fluent asks of you as input is indeed a volumetric source term, but that is based on the zone volume, not the cell volume. I.e., if you want to prescribe a constant source of 100N on a volume of 2m^3, you should enter a value of 50. Fluent will then solve the source term on a cell basis (i.e., a smaller cell will have a smaller net source (in N)). In your case, when dividing by the cell-volume, you're basically saying all cells have an identical net source (larger volumetric source for smaller cells), regardless of their size. Thus resulting in a larger velocity in smaller cells

So, I guess, basically remove the area and volume factors?

Thank you for your replies!

I see that helps and explains a lot.

So let's say I want a source term that is dependent on the velocity of each cell. What I'm trying to do is force the momentum equation to converge to the u_target velocity in the region of the source.

S = (u_gust_total * sin(2* pi * t) - C_U(c,t));

say that u_gust_total is 85 m/s and the free stream is 80 m/s, so the C_U values in the source region are around 80 as well but slightly different varying. The expression above would then be sufficient and I wouldn't need a volume or time-step term?

Thank you very much,
regards!

You can do something like that, but make sure your units are correct. Your source needs to be of unit N/m^3, so do something like:

S = C_R(c,t)*(u_gust_total - C_U(c,t))/tau,

Here, tau is a time-scale, giving you some flexibility on how fast it will converge towards your desired velocity.

Thank you so much for your help, I will try it and come back with the result!

Best regards,
Johannes

I'trying to do something similar in my project. I'm trying to set the pressure drop in a tube with a source term UDF. But if I want the pressure drop to be 3 Pa, for example, should I set the source equal to 3 in all cells? I know the source is in units of Pa/m. So I guess my quetion is: if I want my pressure drop to be 3, what value should I set all the cells of the cell-zone.


My geometry is just a tube, with 0-shear walls, and inlet and an outlet, a separate cell-zone in the middle(1 cell thick and 0.02 m thick), and a porous-jump with the the dP at zero(just to capture particles).

What I tried was dividing by the thickness of the cells. This works when the orientation of the geometry is parallel to the x axis. But whenever I try different orientations or geometries (like a cilindrical surface) the values go crazy.

please help

I am still not an expert on UDF:s nor Fluid dynamics, but couldn't you set pressure boundary conditions at inlet and outlet with a difference of 3 pa?

Maybe I have misunderstood what it is you actually want though :)

kkkkkkkk
well, yes, I think I could do that. But I'm trying to simulate a filter, and this source term is supposed to stand in for the resistance to fluid flow that the filter has. So I'm testing my code and the difference in pressure is my tell if the code is working. While I can set the difference in pressure to whatever I want, that won't mean that the code is working.

I'd say your input value is your 3 Pa divided by the thickness of the complete cell zone on which you apply this source, since Pa = N/m^2 already. So dividing by the thickness makes it N/m^3, ensuring you've got a drop of 3 Pa over your cell zone. Dividing by the thickness of the cells would create a drop of 3Pa over EACH cell. This is fine if the zone is one layer thick, but otherwise not.

As for different orientations of your geometry, the used thickness is the one perpendicular to the normal of the areas). Also, I'd guess you need to properly decompose it in x,y and z components if it's not parallel to either axis.

Actually, I got something that works pretty well by just dividing the target pressure by the cell thickness... I have a different problem now: let's say I inject 3 kg/s in an area of 1m2; my pressure is a function of the mass area density (so like dP=C1+(mass density)/C2). The problem is if I want to predict the pressure value after 1 second after first impact, I'd like to plug in 3kg/m2 into this formula and get the dP. That makes sense right? And this works pretty well if the mass distribution is even (when I force it to be completely even the formula works). if there is any variation though, the simulated value is a bit lower than the expected one. And this makes sense I guess: If there is a low density area, the flow will tend to go through there lowering the pressure in the process.

Anyways, my question is how can I predict the global pressure drop. I thought about trying to use a different mass density average (like a harmonic average, maybe?) because I need and average that is lower.

Any help and ideas are appreciated.

Dear Johannes Hall
Can you please update about the final udf used for your above referred problem?
Any help in this regard will be much appreciated.
Thankyou and Regards.

Y momentum source
 

Dear all,

Can anyone please help why is the Y-velocity in source region not coming equal to 15 m/s. I think i made a mistake somewhere in the Y-momentum source UDF?
Thanks in advance
# include "udf.h"

DEFINE_SOURCE(V-gust,cell,thread,dS,eqn)
{
real source=0.0;
real XGS=-7.5; /* X-center of source*/
real YGS=0.0; /* Y-center of source*/
real ZGS=0.0; /* Z-center of source*/
real X= 1.0; /* X, Y, Z dimensions are such so as to make the source-Volume=1m^3*/
real Y= 1.0;
real Z= 1.0;
real t = CURRENT_TIME;

real centroid[3];
C_CENTROID(centroid,cell,thread);
real x_loc=centroid[0];
real y_loc=centroid[1];
real z_loc=centroid[2];
if (t>1.0 && t<1.5)
{
if (fabs(x_loc-XGS)<X && fabs(y_loc-YGS)<Y && fabs(z_loc-ZGS)<Z)
{
real rho=C_R(cell,thread);
source=rho*(15-C_V(cell,thread))/CURRENT_TIMESTEP; /*15(m/s) is perturbation in Y direction*/
dS[eqn]=0.0;
}
}
else
{
source=0.0;
dS[eqn]=0.0;
}
return source;
}