[Baum]

Hello all,


I tried finding the sum of one of my species (KOH) flowing into each cell, but got strange results. So to double check my code, I summed up all flows entering and leaving all faces of my cells, and the numbers don't add up at all.

Here's the relevant code snippet:

begin_c_loop(c,t)
{ 

[...]

c_face_loop(c, t, n)        // loops over all faces of the cell
         {
            f = C_FACE(c,t,n);
            tf = C_FACE_THREAD(c,t,n);
            if(!BOUNDARY_FACE_THREAD_P(tf))        // if the face is pointing to a wall etc., it sometimes crashes if we try to find the C_YI of c1. 
            {
                if(F_FLUX(f,tf)<0) // if the flow enters the cell (< 0 = entering, > 0 = leaving). 
                {
                    SumOfFluxes += F_FLUX(f,tf)*F_YI(f,t,index_koh); // in kg/s.
                }
                else     // if the flow leaves the cell
                {
                    SumOfFluxesOut += F_FLUX(f,tf)*F_YI(f,t,index_koh); // outflow in kg/s
                }
                SumOfTotalFluxes += F_FLUX(f,tf);
            }
         }


 C_UDMI(c,t,10) = SumOfFluxes + SumOfFluxesOut;
 C_UDMI(c,t,11) = SumOfTotalFluxes;


 }
  end_c_loop(c,t)

There are no reactions or other sources in the field. Comparing C_UDMI 10 (which in my opinion should be equal to the total mass "lost" in the cell) with the mass imbalance, it shows vastly different values (order of 1e-4 for the UDMI, but < 1e-12 for mass imbalance).

I noticed that F_YI crashes Fluent if I try to access it with the "tf" face thread pointer, but F_FLUX crashes if I give it the "t" thread pointer. Is there a mistake that I'm missing?

[Baum]

Here's a quick update on this: F_FLUX seems to give the right flux over each face, but not always the right sign (for some reason the face normal direction is not always pointing the same way (in/out of the cell)). Therefore, using F_FLUX without an additional direction vector can give the wrong results. Additionally, F_YI does indeed need the "tf" face thread pointer, but as I said, this crashes Fluent. This is due to the fact that F_YI, similar to F_V and others, is only defined on boundary faces, not on internal faces.

I don't know how to fix my problem yet, but maybe someone in the future can use this information.

[Baum]

Second update: At least with the overall mass balance, this code below seems to work (-ish... it still shows slightly different values to the mass imbalance, but it's very close) . It takes into account the fact that the direction of the face normal vector and the F_FLUX sign correlate.

Since F_YI is not defined on the faces, I'd have to add a way to calculate this from here, but for now this can be closed.

thread_loop_c(t,d)            #tdr
    {
        begin_c_loop(c,t)
        {
            SumOfFluxes = 0;
            SumOfFluxesIn = 0;
            SumOfFluxesOut = 0; 
            
            c_face_loop(c, t, n)     
            {
                f = C_FACE(c,t,n);
                tf = C_FACE_THREAD(c,t,n);
                F_CENTROID(centroid_f,f,tf);
                C_CENTROID(centroid_c,c,t);
                F_AREA(A,f,tf);
                NV_VV(direction_vector, =, centroid_f, -, centroid_c);  // Vector FROM cell center TO face center. Used for dot product.
                dotproduct = -NV_DOT(direction_vector,A); // now has the same sign as the dot product of direction vector and F_FLUX.
                    
                if(F_FLUX(f,tf)*dotproduct < 0) // if flow is OUT OF cell, also removes boundary faces with FLUX = 0
                {
                    SumOfFluxes += F_FLUX(f,tf)*dotproduct/ABS(dotproduct);
                    SumOfFluxesOut += F_FLUX(f,tf)*dotproduct/ABS(dotproduct); 
                } 
                else if(F_FLUX(f,tf)*dotproduct > 0) 
                {
                    SumOfFluxes += F_FLUX(f,tf)*dotproduct/ABS(dotproduct); 
                    SumOfFluxesIn += F_FLUX(f,tf)*dotproduct/ABS(dotproduct);
                }
            }
            C_UDMI(c,t,0) = SumOfFluxes;        // for analysis
            C_UDMI(c,t,1) = SumOfFluxesIn;        // for analysis
            C_UDMI(c,t,2) = SumOfFluxesOut;        // for analysis
        }    
        end_c_loop(c,t)

     }