CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > ANSYS

UDF for mass balances: Integral over faces doesn't equal mass balance?

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 11, 2022, 14:22
Post UDF for mass balances: Integral over faces doesn't equal mass balance?
  #1
New Member
 
Join Date: Dec 2019
Posts: 29
Rep Power: 5
Baum is on a distinguished road
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:

Code:
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 is offline   Reply With Quote

Old   September 5, 2022, 11:35
Default
  #2
New Member
 
Join Date: Dec 2019
Posts: 29
Rep Power: 5
Baum is on a distinguished road
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 is offline   Reply With Quote

Old   September 7, 2022, 08:22
Default
  #3
New Member
 
Join Date: Dec 2019
Posts: 29
Rep Power: 5
Baum is on a distinguished road
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.

Code:
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)

     }
Baum is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[Other] Can't Shake Erros: patch type 'patch' not constraint type 'empty' BrendaEM OpenFOAM Meshing & Mesh Conversion 12 April 3, 2022 19:32
[snappyHexMesh] snappyHexMesh stuck when snap is turned on yukuns OpenFOAM Meshing & Mesh Conversion 3 February 2, 2021 14:05
[snappyHexMesh] Help with Snappy: no layers growing GianF OpenFOAM Meshing & Mesh Conversion 2 September 23, 2020 09:26
decomposePar problem: Cell 0contains face labels out of range vaina74 OpenFOAM Pre-Processing 37 July 20, 2020 06:38
DecomposePar unequal number of shared faces maka OpenFOAM Pre-Processing 6 August 12, 2010 10:01


All times are GMT -4. The time now is 17:12.