stardis-solver

Solve coupled heat transfers
git clone git://git.meso-star.fr/stardis-solver.git
Log | Files | Refs | README | LICENSE

commit 04ce087a7761b2dab196f0c95f2a5b9047523ee8
parent a7f6b5a129924c0e6e598402cbb8a1d83ebd04e4
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  1 Apr 2022 12:07:32 +0200

Add net flux in combination with other limit conditions

Diffstat:
Msrc/sdis_heat_path_boundary_Xd.h | 15---------------
Msrc/sdis_heat_path_boundary_Xd_c.h | 44++++++++++++++++++++++++++++++++++++++++++++
Msrc/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h | 18+++++++++++++++++-
Msrc/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h | 41+++++++++++++++++++++++++++++++++++++++++
Msrc/sdis_heat_path_boundary_c.h | 31++++++++++++++++++++++++++++++-
Msrc/sdis_interface.c | 9+--------
Msrc/sdis_medium_c.h | 2+-
7 files changed, 134 insertions(+), 26 deletions(-)

diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h @@ -39,7 +39,6 @@ XD(boundary_path) struct sdis_interface* interf = NULL; struct sdis_medium* mdm_front = NULL; struct sdis_medium* mdm_back = NULL; - struct sdis_medium* mdm = NULL; double tmp; res_T res = RES_OK; ASSERT(scn && ctx && rwalk && rng && T); @@ -70,20 +69,6 @@ XD(boundary_path) goto exit; } - /* Check if the boundary flux is known. Note that currently, only solid media - * can have a flux as limit condition */ - mdm = interface_get_medium(interf, frag.side); - if(sdis_medium_get_type(mdm) == SDIS_SOLID) { - const double phi = interface_side_get_flux(interf, &frag); - if(phi != SDIS_FLUX_NONE) { - res = XD(solid_boundary_with_flux_path) - (scn, ctx, &frag, phi, rwalk, rng, T); - if(res != RES_OK) goto error; - - goto exit; - } - } - mdm_front = interface_get_medium(interf, SDIS_FRONT); mdm_back = interface_get_medium(interf, SDIS_BACK); diff --git a/src/sdis_heat_path_boundary_Xd_c.h b/src/sdis_heat_path_boundary_Xd_c.h @@ -872,4 +872,48 @@ error: goto exit; } +res_T +XD(handle_net_flux) + (const struct sdis_scene* scn, + const struct handle_net_flux_args* args, + struct XD(temperature)* T) +{ + double flux_term; + double phi; + res_T res = RES_OK; + CHK(scn && T); + CHK(args->interf && args->frag); + CHK(args->h_cond >= 0); + CHK(args->h_conv >= 0); + CHK(args->h_radi >= 0); + CHK(args->h_cond + args->h_conv + args->h_radi > 0); + + phi = interface_side_get_flux(args->interf, args->frag); + if(phi == SDIS_FLUX_NONE) goto exit; /* No flux. Do nothig */ + + if(args->picard_order > 1 && phi != 0) { + log_err(scn->dev, + "%s: invalid flux '%g' W/m^2. Could not manage a flux != 0 when the " + "picard order is not equal to 1; Picard order is currently set to %lu.\n", + FUNC_NAME, phi, (unsigned long)args->picard_order); + res = RES_BAD_ARG; + goto error; + } + + flux_term = 1.0 / (args->h_cond + args->h_conv + args->h_radi); + T->value += phi * flux_term; + + /* Register the net flux term */ + if(args->green_path) { + res = green_path_add_flux_term + (args->green_path, args->interf, args->frag, flux_term); + if(res != RES_OK) goto error; + } + +exit: + return res; +error: + goto exit; +} + #include "sdis_Xd_end.h" diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h @@ -119,6 +119,9 @@ XD(solid_fluid_boundary_picard1_path) struct ssp_rng* rng, struct XD(temperature)* T) { + /* Input argument used to handle the net flux */ + struct handle_net_flux_args handle_net_flux_args = HANDLE_NET_FLUX_ARGS_NULL; + /* Input/output arguments of the function used to sample a reinjection */ struct XD(sample_reinjection_step_args) samp_reinject_step_args = XD(SAMPLE_REINJECTION_STEP_ARGS_NULL); @@ -150,6 +153,7 @@ XD(solid_fluid_boundary_picard1_path) double lambda; /* Solid conductivity */ double delta_boundary; /* Orthogonal reinjection dst at the boundary */ double delta; /* Orthogonal fitted reinjection dst at the boundary */ + double delta_m; /* delta in meter */ double r; struct sdis_heat_vertex hvtx = SDIS_HEAT_VERTEX_NULL; @@ -209,10 +213,11 @@ XD(solid_fluid_boundary_picard1_path) /* Define the orthogonal dst from the reinjection pos to the interface */ delta = reinject_step.distance / sqrt(DIM); + delta_m = delta * scn->fp_to_meter; /* Compute the convective, conductive and the upper bound radiative coef */ h_conv = interface_get_convection_coef(interf, frag); - h_cond = lambda / (delta * scn->fp_to_meter); + h_cond = lambda / delta_m; h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * ctx->That3 * epsilon; /* Compute a global upper bound coefficient */ @@ -222,6 +227,17 @@ XD(solid_fluid_boundary_picard1_path) p_conv = h_conv / h_hat; p_cond = h_cond / h_hat; + /* Handle the net flux if any */ + handle_net_flux_args.interf = interf; + handle_net_flux_args.frag = frag; + handle_net_flux_args.green_path = ctx->green_path; + handle_net_flux_args.picard_order = get_picard_order(ctx); + handle_net_flux_args.h_cond = h_cond; + handle_net_flux_args.h_conv = h_conv; + handle_net_flux_args.h_radi = h_radi_hat; + res = XD(handle_net_flux)(scn, &handle_net_flux_args, T); + if(res != RES_OK) goto error; + /* Fetch the last registered heat path vertex */ if(ctx->heat_path) hvtx = *heat_path_get_last_vertex(ctx->heat_path); diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h @@ -26,6 +26,41 @@ #include "sdis_Xd_begin.h" /******************************************************************************* + * Non generic helper functions + ******************************************************************************/ +#ifndef SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARD_N_H +#define SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARD_N_H + +static INLINE res_T +check_net_flux + (struct sdis_scene* scn, + const struct sdis_interface* interf, + const struct sdis_interface_fragment* frag, + const size_t picard_order) +{ + double phi; + res_T res = RES_OK; + ASSERT(scn && interf && frag && picard_order > 1); + + phi = interface_side_get_flux(interf, frag); + if(phi != SDIS_FLUX_NONE && phi != 0) { + log_err(scn->dev, + "%s: invalid flux '%g' W/m^2. Could not manage a flux != 0 when the " + "picard order is not equal to 1; Picard order is currently set to %lu.\n", + FUNC_NAME, phi, (unsigned long)picard_order); + res = RES_BAD_ARG; + goto error; + } + +exit: + return res; +error: + goto exit; +} + +#endif /* SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARD_N_H */ + +/******************************************************************************* * Generic helper functions ******************************************************************************/ static INLINE res_T @@ -139,6 +174,7 @@ XD(solid_fluid_boundary_picardN_path) ASSERT(scn && rwalk && rng && T && ctx); ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag)); + /* Fetch the Min/max temperature */ Tmin = ctx->Tmin; Tmin2 = ctx->Tmin2; @@ -159,6 +195,11 @@ XD(solid_fluid_boundary_picardN_path) ASSERT(fluid->type == SDIS_FLUID); } + /* Check that no net flux is set for this interface since the provided + * picardN algorithm does not handle it */ + res = check_net_flux(scn, interf, frag, get_picard_order(ctx)); + if(res != RES_OK) goto error; + /* Setup a fragment for the fluid side */ frag_fluid = *frag; frag_fluid.side = fluid_side; diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h @@ -68,7 +68,7 @@ SAMPLE_REINJECTION_STEP_ARGS_NULL_3d = SAMPLE_REINJECTION_STEP_ARGS_NULL___3d; #define REINJECTION_STEP_NULL___2d {S2D_HIT_NULL__, {0,0}, 0} #define REINJECTION_STEP_NULL___3d {S3D_HIT_NULL__, {0,0,0}, 0} -static const struct reinjection_step_2d +static const struct reinjection_step_2d REINJECTION_STEP_NULL_2d = REINJECTION_STEP_NULL___2d; static const struct reinjection_step_3d REINJECTION_STEP_NULL_3d = REINJECTION_STEP_NULL___3d; @@ -140,6 +140,35 @@ solid_reinjection_3d struct solid_reinjection_args_3d* args); /******************************************************************************* + * Handle net flux + ******************************************************************************/ +struct handle_net_flux_args { + struct sdis_interface* interf; + const struct sdis_interface_fragment* frag; + struct green_path_handle* green_path; + + size_t picard_order; + double h_cond; /* Convective coefficient, i.e. lambda/delta */ + double h_conv; /* Condutive coefficient */ + double h_radi; /* Radiative coefficient */ +}; +#define HANDLE_NET_FLUX_ARGS_NULL__ {NULL,NULL,NULL,0,0,0,0} +static const struct handle_net_flux_args HANDLE_NET_FLUX_ARGS_NULL = + HANDLE_NET_FLUX_ARGS_NULL__; + +extern LOCAL_SYM res_T +handle_net_flux_2d + (const struct sdis_scene* scn, + const struct handle_net_flux_args* args, + struct temperature_2d* T); + +extern LOCAL_SYM res_T +handle_net_flux_3d + (const struct sdis_scene* scn, + const struct handle_net_flux_args* args, + struct temperature_3d* T); + +/******************************************************************************* * Boundary sub-paths ******************************************************************************/ extern LOCAL_SYM res_T diff --git a/src/sdis_interface.c b/src/sdis_interface.c @@ -73,6 +73,7 @@ check_interface_shader FOR_EACH(i, 0, 2) { switch(type[i]) { + case SDIS_FLUID: /* No constraint */ break; case SDIS_SOLID: if(shaders[i]->emissivity || shaders[i]->specular_fraction @@ -84,14 +85,6 @@ check_interface_shader caller_name); } break; - case SDIS_FLUID: - if(shaders[i]->flux) { - log_warn(dev, - "%s: the interface side toward a fluid can't have a flux property. " - "The shader's pointer function for this attribute should be NULL.\n", - caller_name); - } - break; default: FATAL("Unreachable code.\n"); break; } } diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h @@ -176,7 +176,7 @@ fluid_get_properties /******************************************************************************* * Solid local functions ******************************************************************************/ -DEFINE_MDM_CHK_PROP_FUNC(solid, calorific_capacity,0, INF, 0, 1) +DEFINE_MDM_CHK_PROP_FUNC(solid, calorific_capacity, 0, INF, 0, 1) DEFINE_MDM_CHK_PROP_FUNC(solid, thermal_conductivity, 0, INF, 0, 1) DEFINE_MDM_CHK_PROP_FUNC(solid, volumic_mass, 0, INF, 0, 1) DEFINE_MDM_CHK_PROP_FUNC(solid, delta, 0, INF, 0, 1)