| Morice |
June 12, 2020 03:58 |
Divergence & Floating Point Error on using UDF
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;
}
|
