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

LinkBack  Thread Tools  Search this Thread  Display Modes 
August 11, 2022, 14:22 
UDF for mass balances: Integral over faces doesn't equal mass balance?

#1 
New Member
Join Date: Dec 2019
Posts: 29
Rep Power: 5 
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 1e4 for the UDMI, but < 1e12 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? 

September 5, 2022, 11:35 

#2 
New Member
Join Date: Dec 2019
Posts: 29
Rep Power: 5 
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. 

September 7, 2022, 08:22 

#3 
New Member
Join Date: Dec 2019
Posts: 29
Rep Power: 5 
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) } 

Thread Tools  Search this Thread 
Display Modes  


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 PreProcessing  37  July 20, 2020 06:38 
DecomposePar unequal number of shared faces  maka  OpenFOAM PreProcessing  6  August 12, 2010 10:01 