|
[Sponsors] |
Divergence & Floating Point Error on using UDF |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
June 12, 2020, 03:58 |
Divergence & Floating Point Error on using UDF
|
#1 |
Member
Join Date: Oct 2019
Posts: 35
Rep Power: 6 |
Hey guys,
The UDF is an energy and water vapour source UDF and I keep getting a floating-point error and divergence on temperature and water vapour when I try to run the simulation. I have done the math in C++ and Matlab and the math is correct. Maybe there is an issue with the thermodynamics that am not seeing. Honestly, I have no idea what the issue could be. Thanks so much, much appreciated. Code:
#include "udf.h" #include "math.h" #define ra 271 /*aerodynamic resistance*/ #define lad .3033 /*leaf area demsity*/ #define kc .95 /*extinction coefficient of the canopy*/ #define hgt 1.5 /*height of the canopy*/ #define Cp 1006.43 /*specific heat of dry air*/ #define rho 1.225 /*density of dry air*/ #define lai .24271 /*leaf area index*/ #define c3 18.6 #define c4 197.5 #define c5 .31 #define c6 1.2e-6 #define lamda 2.450e6 /*latent heat of vaporisation*/ #define gamma 67.381125 /*psychrometric constant*/ #define ks 0.7 /*coefficient of shading screen*/ /*VARIABLES THAT NEED TO BE MODIFIED BEFORE EVERY SIMULATION*/ #define rgo 43.28 /*global radiation above the canopy*/ DEFINE_INIT(initial_values, domain) { cell_t c; Thread* t; thread_loop_c(t, domain) { begin_c_loop(c, t) { C_UDMI(c, t, 0) = 288.15; /*Ta*/ C_UDMI(c, t, 1) = 1699.81; /*Psat*/ C_UDMI(c, t, 2) = 289.15; /*Temperature of leaf in cropzone1*/ C_UDMI(c, t, 3) = 60.497; /*rs1*/ C_UDMI(c, t, 4) = 331.495; /*rt1*/ C_UDMI(c, t, 5) = 6.19434; /*ET1, Qlat1*/ C_UDMI(c, t, 6) = 292.226; /*Tlr1; Real leaf temp cropzone1*/ C_UDMI(c, t, 7) = 7.66834e-7; /*cropzone1 water vapour mass fractions source*/ C_UDMI(c, t, 8) = 289.15; /*Temperature of leaf in cropzone2*/ C_UDMI(c, t, 9) = 60.4947; /*rs2*/ C_UDMI(c, t, 10) = 331.495; /*rt2*/ C_UDMI(c, t, 11) = 7.66834e-7; /*cropzone2 water vapour mass fractions source*/ C_UDMI(c, t, 12) = 292.226; /*Tlrl; Real leaf temp crop zone2*/ C_UDMI(c, t, 13) = 6.19434; /*ET2, Qlat2*/ } end_c_loop(c, t) } } DEFINE_ADJUST(ta_calc, d) { Thread* t; cell_t c; thread_loop_c(t, d) { if(THREAD_ID(t)==7) begin_c_loop(c,t) { double Ta,Psat; Ta += C_T(c, t) * C_VOLUME(c,t); C_UDMI(c, t, 0) = Ta; Psat = 610.78 * exp((17.25 * (Ta - 273.15)) / (237.8 + (Ta - 273.15))); C_UDMI(c, t, 1) = Psat; } end_c_loop(c,t) } } DEFINE_SOURCE(cropzone1_h20_source, c, t, dS, eqn) { double a, rs1, rt1, hcv1, Tl1, Tl1r, Plsat1, vpd1, ET1, Sw1, St1; real con, source; real x[ND_ND]; C_CENTROID(x, c, t); con = rgo * ks * exp(-kc * lad * (((hgt - x[1]) / hgt))); Tl1 += C_T(c, t) * C_VOLUME(c, t); C_UDMI(c, t, 2) = Tl1; Plsat1 = 610.78 * exp((17.25 * (Tl1 - 273.15)) / (237.8 + (Tl1 - 273.15))); vpd1 = Plsat1 - C_UDMI(c, t, 1); a = con / (2 * lai); rs1 = c3 * ((a + c4) / (a + c5)) * (1 + (c6 * vpd1 * vpd1)); C_UDMI(c, t, 3) = rs1; rt1 = ra + rs1; C_UDMI(c, t, 4) = rt1; ET1 = ((rho * Cp) / gamma) * (vpd1 / rt1); C_UDMI(c, t, 5) = ET1; hcv1 = (rho * Cp) / ra; Tl1r = (con - ET1) / ((2 * hcv1) + C_UDMI(c, t, 0)); C_UDMI(c, t, 6) = Tl1r; source = (rho * Cp * lad * vpd1) / (gamma * lamda * rt1); C_UDMI(c, t, 7) = source; dS[eqn] = (rho * Cp * lad) / (gamma * lamda * rt1); return source; } DEFINE_SOURCE(cropzone1_energy_source, c, t, dS, eqn) { real source; source = - C_UDMI(c, t, 7) * lad; dS[eqn] = -(rho * Cp * lad) / (gamma * C_UDMI(c, t, 4)); return source; } DEFINE_SOURCE(cropzone2_h20_source, c, t, dS, eqn) { double b, rs2, rt2, hcv2, Tl2, Tl2r, Plsat2, vpd2, ET2, Sw2, St2; real x[ND_ND]; real source, con; C_CENTROID(x, c, t); con = rgo * ks * exp(-kc * lad * (((hgt - x[1]) / hgt))); Tl2 += C_T(c, t) * C_VOLUME(c, t); C_UDMI(c, t, 8) = Tl2; Plsat2 = 610.78 * exp((17.25 * (Tl2 - 273.15)) / (237.8 + (Tl2 - 273.15))); vpd2 = Plsat2 - C_UDMI(c, t, 1); b = con / (2 * lai); rs2 = c3 * ((b + c4) / (b + c5)) * (1 + (c6 * vpd2 * vpd2)); C_UDMI(c, t, 9) = rs2; rt2 = ra + rs2; C_UDMI(c, t, 10) = rt2; ET2 = ((rho * Cp) / gamma) * (vpd2 / rt2); C_UDMI(c, t, 11) = ET2; hcv2 = (rho * Cp) / ra; Tl2r = (con - ET2) / ((2 * hcv2) + C_UDMI(c, t, 0)); C_UDMI(c, t, 12) = Tl2r; source = (rho * Cp * lad * vpd2) / (gamma * lamda * rt2); C_UDMI(c, t, 13) = source; dS[eqn] = (rho * Cp * lad) / (gamma * lamda * rt2); return source; } DEFINE_SOURCE(cropzone2_energy_source, c, t, dS, eqn) { real source; source = - C_UDMI(c, t, 13) * lad; dS[eqn] = -(rho * Cp * lad) / (gamma * C_UDMI(c, t, 10)); return source; } |
|
June 15, 2020, 03:53 |
Divergence
|
#2 |
Senior Member
|
Do you get divergence at the very first iteration? In that case, you need to ensure that all the denominators are non-zero and all the sources are comparatively small. Since you are using field values in the UDF, a good initial guess is also important. So, you may try running the solution for a few iterations before hooking the source UDFs.
__________________
Regards, Vinerm PM to be used if and only if you do not want something to be shared publicly. PM is considered to be of the least priority. |
|
June 15, 2020, 23:48 |
Divergence
|
#3 | |
Member
Join Date: Oct 2019
Posts: 35
Rep Power: 6 |
Quote:
Thank you for the response, much appreciated. The divergence is on the very 1st iteration. I realized that a good initial guess is important so the initial values I used are those values that I obtained from running the equations in Matlab making sure that the times are the same. I was thinking of running the solution for a few iterations before hooking the source UDFs but I was not able to find a way to do that. How is that done? do I stop the calculation, hook up the UDFs then continue the calculation? |
||
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Divergence Problem Conjugated Heat Transfer | Christoph.Ha | FLUENT | 27 | December 5, 2017 10:46 |
fluent divergence for no reason | sufjanst | FLUENT | 2 | March 23, 2016 16:08 |
[blockMesh] error EOF in blockMesh | Ahmed Khattab | OpenFOAM Meshing & Mesh Conversion | 7 | May 17, 2012 00:37 |
Primitive error at node- Floating point exception / continuity divergence issues | phales15 | FLUENT | 0 | April 30, 2011 19:05 |
divergence in AMG and floating point | grzes | FLUENT | 0 | February 17, 2006 13:22 |