Hello,

I have been given this UDF by a colleague. It is used for the UDS transport equation that describes a species transport. I have a problem with the UDF, but before taking the issue further I would like some validation that my objections are accurate. My problem is that in the first of the if-else statements dictates weather the face f is a wall face, i.e. (c1 <0). If this is the case, then there should be no flux *through* f, as species cannot travel through a wall. However, the UDF below should allow species to flux *through* face f as the value of "float flux"can evaluate to something other than zero. Am I correct?

Any comments on this matter are, as always, greatly appreciated.

#include "udf.h"

#include "sg.h"

#define ALPHA 4e-5

DEFINE_UDS_FLUX(conv_fluxes, f, thread, eqn) {

cell_t c0, c1;

Thread *tc0, *tc1;

float A[ND_ND];

float dphidx, dphidy;

float elec, conv, flux;

int phi = 0;

c0 = F_C0(f, thread);

tc0 = THREAD_T0(thread);

c1 = F_C1(f, thread);

tc1 = THREAD_T1(thread);

if( c1 < 0 )

{

dphidx = C_UDSI_G(c0,tc0,phi)[0];

dphidy = C_UDSI_G(c0,tc0,phi)[1];

F_AREA(A, f, thread);

elec = ALPHA*(dphidx*A[0]+dphidy*A[1]);

conv = F_FLUX(f, thread);

if( eqn == 4) flux = conv + elec;

else

flux = 0.0;

}

else

{

dphidx = 0.5*(C_UDSI_G(c0,tc0,phi)[0] + C_UDSI_G(c1,tc1,phi)[0]);

dphidy = 0.5*(C_UDSI_G(c0,tc0,phi)[1] + C_UDSI_G(c1,tc1,phi)[1]);

F_AREA(A, f, thread);

elec = ALPHA*(dphidx*A[0]+dphidy*A[1]);

conv = F_FLUX(f, thread);

if( eqn == 4)

flux = conv + elec;

else

flux = 0.0;

}

return flux; }

Hello,

first as a clarification, a face is not neccessary a wall. a face is usually the connection between two cells, namely cell c0 and its neighbour c1. if there is no c1 (c1<0) we are at a boundary wall. the flux function will state how the scalar gets through the face from cell c0 to c1. if you want flux through wall you will have to set is seperatly in the wall-bc.

when setting up a case with uds you can choose one of two default flux functions supplied by fluent. you can either choose "none" which is equal to only diffusive flux (it does not mean no flux at all !). the diffusion coefficient will be set in the material panel as constant or DEFINE_DIFFUSIVITY udf. you second option is "mass-flow-rate" which will compute convective flux and diffusive flux. you third option is to supply your own flux fuction with udf.

now to your udf. the udf uses the macro F_FLUX to get the convective flux through a face and adds a diffusive part, with ALPHA as diffusion coefficient. the computed flux is only applied to scalar number 4 ( if(eqn==4) ) the flux for all other scalars is set to zero. the if( c1<0 ) statement will only look if cell c0 has a neighbour or not. if there is no neighbour it will use the cell c0 values for velocity and scalar gradient to compute the flux, else the mean values between c0 and c1 are taken. further references can be found in the udf manual chapter 11 and chapter 4.3.21(fully commented udf).

good luck RoM

Thank you for the explanation of the UDF. I do understand how it works, my question is:

If c1<0 resolves true, then is there a UDS (species) flux throug a boundary wall face? I believe the answer is yes.

just think of the bounday conditions for velocity and scalars gradients at walls. velocity->no slip->velocity=0 at wall->no convective flux through that wall. scalar gradient=0 at wall->no diffusion through that wall.so your udf will compute zero flux for walls.

RoM

OK, now we are getting somewhere. Why will the scalar gradient=0 at the wall?

ist just an assumption fluent has to make, since its doesnt know any better . fluent does not know whats behind that wall and the simplest approach ist so think that salar concentration outside = concentration inside which results in a gradient=0 bc. there could be oufcourse a gradient and therefor a flux throught that wall but you have to tell fluent. you can do that by setting a fluxes for species, energy and ofcourse scalars as wall-bc. i wound not touch the overall flux function to force fluent to calculate a flux through a wall. this would probably lead to errors in scalar conservation.

RoM

I think that we may have got a little confused at some point, probably by the fact that the diffusive part of the species/UDS transport equation is based not upon species concentration, but another variable. The speices that the flux equation that the UDF applies to is UDS4. The diffusive part of the equation is based on the gradient of UDS0.

However, if a numerical output from Fluent is put to screen, it shows that the variable "ele" is not zero for all of the wall (c1<0) faces. I believe that this is due to the fact that the neighbour cell to the face "f" has a gradient of UDS0. The result is that the UDF is predicting a UDS4 (i.e. Species) flux *through the wall face*. This would be an unphysical situation.

Even if the flux is based on the gradient of an other scalar, the grad=0 bc on walls should apply to all scalar values (energy, species, uds, ...) if not statet otherwise.

Maybe the value of uds-0 is set on some other point (another udf, bc) ?

RoM

"Even if the flux is based on the gradient of an other scalar, the grad=0 bc on walls should apply to all scalar values"

This is true, however, the flux through the wall (c1<0) face isn't determined by the boundary conditions, it is determined by the gradient of UDS0 across the neighbour cell, c0. This means that the flux through the wall face f can be set to a non-zero value regardless of the boundary conditions. To my mind, this makes the UDF unrealistic. Would that be your opinion?

You are right. Its not ok to take just the cell gradient to calculate the diffusive flux. Instead you will have to evaluate the gradient in the direction from center of cell0 to center of cell1 or center of cell0 and center of the boundary face in case of a boundary wall. Equation 5.2-3 in the udf manual wil do this and i added it to your udf. Now everything should work fine. I had not time to test it so it would be nice if could drop me a line how it worked.

RoM

#include "udf.h"

#include "sg.h"

#define ALPHA 4e-5

DEFINE_UDS_FLUX(conv_fluxes, f, thread, eqn)

{

/* cell geometry data */

real ds,A_by_es,diff_flux,cross_flux;

real es[ND_ND];

real dr0[ND_ND];

real dr1[ND_ND];

real A[ND_ND];

cell_t c0, c1;

Thread *tc0, *tc1;

real mean_grad[ND_ND];

real elec, conv, flux;

int phi = 0;

c0 = F_C0(f, thread);

tc0 = THREAD_T0(thread);

c1 = F_C1(f, thread);

tc1 = THREAD_T1(thread);

if( c1 < 0 ) /* boundary face */

{

/* get geometry data */

BOUNDARY_FACE_GEOMETRY(f,thread,A,ds,es,A_by_es,dr 0);

/* are face values of phi available ? */

if( THREAD_STORAGE(thread,SV_UDS_I(phi)) != NULL )

diff_flux = ALPHA*(F_UDSI(f,thread,phi)-C_UDSI (c0,tc0,phi))/ds*A_by_es;

else

/* assume that face value = cell value */

diff_flux=0.0;

elec = diff_flux;

conv = F_FLUX(f, thread);

if( eqn == 4) flux = conv + elec;

else

flux = 0.0;

}

else /* interior cell face */

{

/* get geometry data */

INTERIOR_FACE_GEOMETRY(f,thread,A,ds,es,A_by_es,dr 0,dr1);

/* calculate mean gradient of phi between c0 and c1 */

NV_VS_VS(mean_grad,=,C_UDSI_G(c0,tc0,phi),*,0.5,+, C_UDSI_G(c1,tc1,phi),*,0.5);

/* primary diffusive flux */

diff_flux = ALPHA*(C_UDSI(c1,tc1,phi)-C_UDSI(c0,tc0,phi))/ds*A_by_es;

/* cross flux */

cross_flux = ALPHA*(NV_DOT(mean_grad,A)-(NV_DOT (mean_grad,es))*A_by_es);

elec=diff_flux+cross_flux;

conv = F_FLUX(f, thread);

if( eqn == 4)

flux = conv + elec;

else

flux = 0.0;

}

return flux;

}

Woah, that looks like some complicated code. It may take me a little time to digest that. Thanks for taking the time to look over it for me, I very much appreciate it. Unfortunatly I don't have a great deal of time myself as the project is ending in two weeks. Until now, I have been assuming that as no species can pass through a solid wall, if(c1<0){ele=0;}. I have no idea how this will work when the solid wall is set to have a boundary condition for UDS4.

When I implement and test your code I will get get back to you about how it works.

Once again, thanks for your time.

I have been looking at your UDF for a while now and I have been consulting the UDF manual to try and work out what all the NV macros do. Stumbling across equation 5.2.3 I noticed that this is identical to your UDF, including all the wording such as diff and cross. The text suggests that this is only needed for an unstructured grid. As I am using a structured grid, can I use equation 5.2.2, which is much easier for me to understand and is very much like my original UDF.

Even with the new UDF in place, although I haven't tried it yet, I still think that, if c1<0, a non zero value may be returned through a solid wall, which was the original problem. Species can't travel through walls.

For strucured grids you you need only the first part of equation 5.2.3 .For boundary faces phi1 is the face value of your scalar and phi0 the cell value. If fluent sets the face values correct (equal to cell value at walls) the function returns zero flux. For inlet our outlet boundaries the function could return some flux depending on the bc.

RoM

http://www.shef.ac.uk/suretec/broke.jpg

I tried entering your code into a very simple situation. There is a gradient of UDS0 betwen the black box at the top of the picture and directly below. I supplied an inlet from the left and it flows to the right. Where there is a small part there was a stream of high concentraiton UDS1. This jetted into the main stream. As the iterations proceed, the jet bends towards the top of the screen, as expected. Then after a number of iteration (150) the small high UDS1 concentration zone appears at the bottom of the screen. The value of the concentrations then diverges as the iteration proceed. Any clues? (sorry for low quality jpeg)

I had some time to test the "hand calculated gradients" and the it seems the udf is not very stable. I think it calculates gradients to large in areas with high concentration or coarse grids. If you know how large the max gradients are, you could add a limiter to the udf so it doesnt produce those overshoots. Maybe grid refinement in cirtical areas will help too. If everything failes you could try an unsteady simlulation to see how this large concentration area moves within time. It should move along the wall from left to right and vanish or leave trough the exit.

Good Luck

RoM

I have been doing lots of playing about using Both your model and my original one, trying to increase the stability and physical nature. I think that you might be right about grid sizes, however I am still haveing problems with the original issue. If you look a the new image I uploaded:

http://www.shef.ac.uk/suretec/broke.jpg

This is the result of both my original and your UDF. There is a field across the flow where the black box is. The filed diverts the speices due to the gradient. Because when c<1 and there is a gradient at this face, the flux through the face is a non zero value. In this situation I have a speces (UDS1) flux going THROUGH a wall. This should not occur, as the wall is inert. There is less species flowing from the inlet than to the outlet.

Maybe the problem is not the the description of the flux but something else. You could take a closer look at uds-0 and i added a little udf which will copy cell values and face values of uds-0 for a given boundary into user defined memories (cell values to memory 0 and face values to memory 1). If there is any difference between those two values we probably found the reason for the flux.

#include "udf.h"

#define WALL_ID 1

DEFINE_ON_DEMAND(copy_uds0) {

Domain *d;

Thread *ft,*tc0;

face_t f;

cell_t c0;

d = Get_Domain(1);

ft=Lookup_Thread(d,WALL_ID);

begin_f_loop(f,ft)

{

tc0=THREAD_T0(ft);

c0=F_C0(f,ft);

C_UDMI(c0,tc0,0)=C_UDSI(c0,tc0,0);

if( THREAD_STORAGE(ft,SV_UDS_I(0)) != NULL )

C_UDMI(c0,tc0,1)=F_UDSI(f,ft,0);

else

C_UDMI(c0,tc0,1)=-1;

}

end_f_loop(f,ft) }

I put your code in (thanks, I had not used a define on demand before, so that is another thing I have learned). I put in an output to screen line in the UDF so that there was something more quantifiable to look at and their aren't that many faces. The results return a non-zero value. However, I knew that before executing the code. Allow me to contextualise the problem: There is a field that runs perpendicular to the flow, which is solved by UDS0. This is a standard scalar with no convection term. UDS1 is a species that is affected by gradient of the filed. The stream of species should be attracted to the wall where the field is applied, but not pass through it. The same amount of species that enters the inlet should leave parallel outlet. However, when the (c1<0) evaluates to anything but a zero value there is a flux of the species (UDS1) through the wall, which is unrealistic. When I try setting ele to return 0 when c1<0 I get very odd results.

I tried to set up a case similar to your one and had no unwanted fluxes through walls. The scalar tends to accumulate near the wall for low fluid velocities but it never diffused through the wall. Maybe your problem is somewhere in the defintion of the forcefield. If you want my case send me a private mail and i will send it to you. Maybe it helps to track your problem.

RoM

My email is my andrewgarrard then the at symbol uk2.net

However, I don't see how the simulation you set up could be that much differnt from mine - mine is very simple. The boundary conditions I used to create the force field was the application of a fixed value for UDS0 on both sides of the flow field.