commit 7bdfdbc395b9910d65a3907f247c988490e800bb
parent aa08810adb2a925ea0ec2a28fadd480ef3a16d5c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 28 Mar 2023 17:08:11 +0200
Reverse the convention of flux density calculations
Until now, the calculated flux density was the flux leaving the solid.
This is the inverse of an imposed flux which is usually defined as the
net flux entering the solid. This commit corrects this contradiction
by reversing the sign of the flux density returned by the
sdis_solve_[probe_]_boundary_flux functions. The net imposed flux and
the calculated flux densities are therefore now based on the same
convention: they are positive if energy enters the solid.
Diffstat:
5 files changed, 25 insertions(+), 18 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -1342,12 +1342,18 @@ sdis_solve_boundary
const struct sdis_solve_boundary_args* args,
struct sdis_estimator** estimator);
+/* Calculate the flux density in W/m² _entering_ the solid through the given
+ * boundary position, i.e. the flux density is positive or negative if the
+ * solid gains or loses energy, respectively. */
SDIS_API res_T
sdis_solve_probe_boundary_flux
(struct sdis_scene* scn,
const struct sdis_solve_probe_boundary_flux_args* args,
struct sdis_estimator** estimator);
+/* Calculate the average flux density in W/m² _entering_ the solid through the
+ * given boundary surfaces, i.e. the flux density is positive or negative if
+ * the solid gains or loses energy, respectively. */
SDIS_API res_T
sdis_solve_boundary_flux
(struct sdis_scene* scn,
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -817,13 +817,13 @@ XD(solve_boundary_flux)
continue;
} else if(res_simul == RES_OK) { /* Update accumulators */
const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
- /* Convective flux from solid to fluid */
- const double w_conv = hc * (result.Tboundary - result.Tfluid);
- /* Radiative flux from solid to ambient */
+ /* Convective flux from fluid to solid */
+ const double w_conv = hc * (result.Tfluid - result.Tboundary);
+ /* Radiative flux from ambient to solid */
const double w_rad = (result.Tradiative < 0) ?
- 0 : hr * (result.Tboundary - result.Tradiative);
- /* Imposed flux that goes _out_ of the solid */
- const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? -imposed_flux : 0;
+ 0 : hr * (result.Tradiative - result.Tboundary);
+ /* Imposed flux that goes _into_ the solid */
+ const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
const double w_total = w_conv + w_rad + w_imp;
/* Temperature */
acc_temp->sum += result.Tboundary;
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -611,13 +611,14 @@ XD(solve_probe_boundary_flux)
continue;
} else if(res_simul == RES_OK) { /* Update accumulators */
const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
- /* Convective flux from solid to fluid */
- const double w_conv = hc * (result.Tboundary - result.Tfluid);
- /* Radiative flux from solid to ambient */
+ /* Convective flux from fluid to solid */
+ const double w_conv = hc * (result.Tfluid - result.Tboundary);
+ /* Radiative flux from ambient to solid */
const double w_rad = (result.Tradiative < 0) ?
- 0 : hr * (result.Tboundary - result.Tradiative);
- /* Imposed flux that goes _out_ of the solid */
- const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? -imposed_flux : 0;
+ 0 : hr * (result.Tradiative - result.Tboundary);
+ /* Imposed flux that goes _into_ the solid */
+ const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
+ /* Total flux */
const double w_total = w_conv + w_rad + w_imp;
/* Temperature */
acc_temp->sum += result.Tboundary;
diff --git a/src/test_sdis_flux_with_h.c b/src/test_sdis_flux_with_h.c
@@ -259,7 +259,7 @@ main(int argc, char** argv)
printf("Imposed flux (left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE);
OK(sdis_estimator_get_total_flux(estimator, &mc));
printf("Total flux (left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE);
- CHK(eq_eps(PHI2, mc.E, 3*mc.SE));
+ CHK(eq_eps(-PHI2, mc.E, 3*mc.SE));
OK(sdis_estimator_ref_put(estimator));
probe_flux_args.nrealisations = N;
@@ -272,7 +272,7 @@ main(int argc, char** argv)
printf("Imposed flux (probe on left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE);
OK(sdis_estimator_get_total_flux(estimator, &mc));
printf("Total flux (probe on left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE);
- CHK(eq_eps(PHI2, mc.E, 3*mc.SE));
+ CHK(eq_eps(-PHI2, mc.E, 3*mc.SE));
OK(sdis_estimator_ref_put(estimator));
probe_args.nrealisations = N;
diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c
@@ -32,8 +32,8 @@
* with Hrad = 4 * BOLTZMANN_CONSTANT * Tref^3 * epsilon
* T(-X) = Tb
*
- * CF = H * (T - Tf)
- * RF = Hrad * (T - Trad)
+ * CF = H * (Tf - T)
+ * RF = Hrad * (Trad - T)
* TF = CF + RF
*
* with Tf the temperature of the surrounding fluid, lambda the conductivity of
@@ -373,8 +373,8 @@ main(int argc, char** argv)
OK(sdis_interface_ref_put(interf_H));
analyticT = (H*Tf + Hrad*Trad + LAMBDA * Tb) / (H + Hrad + LAMBDA);
- analyticCF = H * (analyticT - Tf);
- analyticRF = Hrad * (analyticT - Trad);
+ analyticCF = H * (Tf - analyticT);
+ analyticRF = Hrad * (Trad - analyticT);
analyticTF = analyticCF + analyticRF;
#define SOLVE sdis_solve_probe_boundary_flux
