|
[Sponsors] |
UDF for mass balances: Integral over faces doesn't equal mass balance? |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
August 11, 2022, 13:22 |
UDF for mass balances: Integral over faces doesn't equal mass balance?
|
#1 |
New Member
Join Date: Dec 2019
Posts: 29
Rep Power: 6 |
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? |
|
September 5, 2022, 10:35 |
|
#2 |
New Member
Join Date: Dec 2019
Posts: 29
Rep Power: 6 |
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, 07:22 |
|
#3 |
New Member
Join Date: Dec 2019
Posts: 29
Rep Power: 6 |
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) } |
|
|
|
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 18:32 |
[snappyHexMesh] snappyHexMesh stuck when snap is turned on | yukuns | OpenFOAM Meshing & Mesh Conversion | 3 | February 2, 2021 13:05 |
[snappyHexMesh] Help with Snappy: no layers growing | GianF | OpenFOAM Meshing & Mesh Conversion | 2 | September 23, 2020 08:26 |
decomposePar problem: Cell 0contains face labels out of range | vaina74 | OpenFOAM Pre-Processing | 37 | July 20, 2020 05:38 |
DecomposePar unequal number of shared faces | maka | OpenFOAM Pre-Processing | 6 | August 12, 2010 09:01 |