1 /*============================================================================
2 * Handle the solidification module with CDO schemes
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <float.h>
35
36 /*----------------------------------------------------------------------------
37 * Local headers
38 *----------------------------------------------------------------------------*/
39
40 #include <bft_error.h>
41 #include <bft_mem.h>
42 #include <bft_printf.h>
43
44 #include "cs_cdofb_scaleq.h"
45 #include "cs_navsto_system.h"
46 #include "cs_parall.h"
47 #include "cs_post.h"
48 #include "cs_solid_selection.h"
49
50 /*----------------------------------------------------------------------------
51 * Header for the current file
52 *----------------------------------------------------------------------------*/
53
54 #include "cs_solidification.h"
55
56 /*----------------------------------------------------------------------------*/
57
58 BEGIN_C_DECLS
59
60 /*=============================================================================
61 * Additional doxygen documentation
62 *============================================================================*/
63
64 /*!
65 \file cs_solidification.c
66
67 \brief Structure and functions handling the solidification module
68 (modified Navier-Stokes + thermal module + transport equations)
69
70 */
71
72 /*=============================================================================
73 * Local macro definitions
74 *============================================================================*/
75
76 #define CS_SOLIDIFICATION_DBG 0
77
78 static const char _state_names[CS_SOLIDIFICATION_N_STATES][32] = {
79
80 "Solid",
81 "Mushy",
82 "Liquid",
83 "Eutectic"
84
85 };
86
87 /*============================================================================
88 * Type definitions
89 *============================================================================*/
90
91 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
92
93 /*============================================================================
94 * Static variables
95 *============================================================================*/
96
97 static cs_solidification_t *cs_solidification_structure = NULL;
98 static cs_real_t cs_solidification_forcing_eps = 1e-10;
99 static cs_real_t cs_solidification_eutectic_threshold = 1e-4;
100
101 static const double cs_solidification_diffusion_eps = 1e-16;
102 static const char _err_empty_module[] =
103 " Stop execution.\n"
104 " The structure related to the solidification module is empty.\n"
105 " Please check your settings.\n";
106
107 /*============================================================================
108 * Private static inline function prototypes
109 *============================================================================*/
110
111 /*----------------------------------------------------------------------------*/
112 /*!
113 * \brief Compute the liquidus temperature knowing the bulk concentration
114 * Assumption of the lever rule.
115 *
116 * \param[in] alloy pointer to a binary alloy structure
117 * \param[in] conc value of the bulk concentration
118 *
119 * \return the computed liquidus temperature
120 */
121 /*----------------------------------------------------------------------------*/
122
123 static inline cs_real_t
_get_t_liquidus(const cs_solidification_binary_alloy_t * alloy,const cs_real_t conc)124 _get_t_liquidus(const cs_solidification_binary_alloy_t *alloy,
125 const cs_real_t conc)
126 {
127 return fmax(alloy->t_eut, alloy->t_melt + alloy->ml * conc);
128 }
129
130 /*----------------------------------------------------------------------------*/
131 /*!
132 * \brief Compute the solidus temperature knowing the bulk concentration
133 * Assumption of the lever rule.
134 *
135 * \param[in] alloy pointer to a binary alloy structure
136 * \param[in] conc value of the bulk concentration
137 *
138 * \return the computed solidus temperature
139 */
140 /*----------------------------------------------------------------------------*/
141
142 static inline cs_real_t
_get_t_solidus(const cs_solidification_binary_alloy_t * alloy,const cs_real_t conc)143 _get_t_solidus(const cs_solidification_binary_alloy_t *alloy,
144 const cs_real_t conc)
145 {
146 if (conc < alloy->cs1)
147 return alloy->t_melt + alloy->ml * conc * alloy->inv_kp;
148 else
149 return alloy->t_eut;
150 }
151
152 /*----------------------------------------------------------------------------*/
153 /*!
154 * \brief Compute the value of eta (Cliq = eta * Cbulk) knowing the bulk
155 * concentration and the phase diagram.
156 * Assumption of the lever rule.
157 *
158 * \param[in] alloy pointer to a binary alloy structure
159 * \param[in] conc value of the bulk concentration
160 *
161 * \return the value of eta
162 */
163 /*----------------------------------------------------------------------------*/
164
165 static inline cs_real_t
_get_eta(const cs_solidification_binary_alloy_t * alloy,const cs_real_t conc)166 _get_eta(const cs_solidification_binary_alloy_t *alloy,
167 const cs_real_t conc)
168 {
169 /* Update eta */
170
171 if (conc > alloy->cs1)
172 /* In this case Cl = C_eut = eta * Cbulk--> eta = C_eut/Cbulk */
173 return alloy->c_eut/conc;
174 else
175 return alloy->inv_kp;
176 }
177
178 /*----------------------------------------------------------------------------*/
179 /*!
180 * \brief Determine in which state is a couple (temp, conc)
181 * Assumption of the lever rule.
182 *
183 * \param[in] alloy pointer to a binary alloy structure
184 * \param[in] temp value of the temperature
185 * \param[in] conc value of the bulk concentration
186 *
187 * \return the state among (liquid, solid, mushy or eutectic)
188 */
189 /*----------------------------------------------------------------------------*/
190
191 static inline cs_solidification_state_t
_which_state(const cs_solidification_binary_alloy_t * alloy,const cs_real_t temp,const cs_real_t conc)192 _which_state(const cs_solidification_binary_alloy_t *alloy,
193 const cs_real_t temp,
194 const cs_real_t conc)
195 {
196 const cs_real_t t_liquidus = _get_t_liquidus(alloy, conc);
197
198 if (temp > t_liquidus)
199 return CS_SOLIDIFICATION_STATE_LIQUID;
200
201 else { /* temp < t_liquidus */
202
203 const cs_real_t t_solidus = _get_t_solidus(alloy, conc);
204 if (temp > t_solidus)
205 return CS_SOLIDIFICATION_STATE_MUSHY;
206
207 else { /* temp < t_solidus */
208
209 if (conc < alloy->cs1 || temp < alloy->t_eut_inf)
210 return CS_SOLIDIFICATION_STATE_SOLID;
211 else
212 return CS_SOLIDIFICATION_STATE_EUTECTIC;
213
214 } /* solidus */
215 } /* liquidus */
216
217 }
218
219 /*----------------------------------------------------------------------------*/
220 /*!
221 * \brief Determine in which state is a tuple (temp, conc, gl) from the
222 * evaluation of its enthalpy. The calling code has to be sure that the
223 * tuple is consistent.
224 * Assumption of the lever rule.
225 *
226 * \param[in] alloy pointer to a binary alloy structure
227 * \param[in] latent_heat value of the latent heat coefficient
228 * \param[in] cp value of the heat capacity
229 * \param[in] temp value of the temperature
230 * \param[in] conc value of the bulk concentration
231 * \param[in] gliq value of the liquid fraction
232 *
233 * \return the state among (liquid, solid, mushy or eutectic)
234 */
235 /*----------------------------------------------------------------------------*/
236
237 static inline cs_solidification_state_t
_which_state_by_enthalpy(const cs_solidification_binary_alloy_t * alloy,const cs_real_t latent_heat,const cs_real_t cp,const cs_real_t temp,const cs_real_t conc,const cs_real_t gliq)238 _which_state_by_enthalpy(const cs_solidification_binary_alloy_t *alloy,
239 const cs_real_t latent_heat,
240 const cs_real_t cp,
241 const cs_real_t temp,
242 const cs_real_t conc,
243 const cs_real_t gliq)
244 {
245 const cs_real_t h_liq = cp*_get_t_liquidus(alloy, conc) + latent_heat;
246 const cs_real_t h = cp*temp + gliq*latent_heat;
247
248 if (h > h_liq)
249 return CS_SOLIDIFICATION_STATE_LIQUID;
250
251 else {
252
253 if (conc > alloy->cs1) { /* Part with eutectic */
254
255 const cs_real_t h_sol = cp*alloy->t_eut;
256 const cs_real_t gleut = (conc - alloy->cs1)*alloy->dgldC_eut;
257 const cs_real_t h_eut = cp*alloy->t_eut + gleut*latent_heat;
258
259 if (h > h_eut)
260 return CS_SOLIDIFICATION_STATE_MUSHY;
261 else if (h > h_sol)
262 return CS_SOLIDIFICATION_STATE_EUTECTIC;
263 else
264 return CS_SOLIDIFICATION_STATE_SOLID;
265
266 }
267 else { /* Part without eutectic */
268
269 const cs_real_t h_sol = cp*(alloy->t_melt+alloy->ml*conc*alloy->inv_kp);
270 if (h > h_sol)
271 return CS_SOLIDIFICATION_STATE_MUSHY;
272 else
273 return CS_SOLIDIFICATION_STATE_SOLID;
274
275 } /* Eutectic or not that is the question ? */
276
277 } /* Liquid ? */
278
279 }
280
281 /*----------------------------------------------------------------------------*/
282 /*!
283 * \brief Compute the derivatives of g_l w.r.t. the temperature and the
284 * bulk concentration when the current state is MUSHY
285 * Assumption of the lever rule.
286 *
287 * \param[in] alloy pointer to a binary alloy structure
288 * \param[in] temp value of the temperature
289 * \param[in] conc value of the bulk concentration
290 * \param[out] dgldT value of the derivative of g_l w.r.t. the temperature
291 * \param[out] dgldC value of the derivative of g_l w.r.t. the concentration
292 */
293 /*----------------------------------------------------------------------------*/
294
295 static inline void
_get_dgl_mushy(const cs_solidification_binary_alloy_t * alloy,const cs_real_t temp,const cs_real_t conc,cs_real_t * dgldT,cs_real_t * dgldC)296 _get_dgl_mushy(const cs_solidification_binary_alloy_t *alloy,
297 const cs_real_t temp,
298 const cs_real_t conc,
299 cs_real_t *dgldT,
300 cs_real_t *dgldC)
301 {
302 const double _dTm = temp - alloy->t_melt;
303 const double _kml = alloy->ml * alloy->inv_kpm1;
304
305 *dgldT = _kml * conc/(_dTm*_dTm);
306 *dgldC = -_kml / _dTm;
307 }
308
309 /*============================================================================
310 * Private function prototypes
311 *============================================================================*/
312
313 /*----------------------------------------------------------------------------*/
314 /*!
315 * \brief Check the convergence of the non-linear algorithm
316 *
317 * \param[in] nl_algo_type type of non-linear algorithm
318 * \param[in] pre_iter previous iterate values
319 * \param[in, out] cur_iter current iterate values
320 * \param[in, out] algo pointer to a cs_iter_algo_t structure
321 *
322 * \return the convergence state
323 */
324 /*----------------------------------------------------------------------------*/
325
326 static cs_sles_convergence_state_t
_check_nl_cvg(cs_param_nl_algo_t nl_algo_type,const cs_real_t * pre_iter,cs_real_t * cur_iter,cs_iter_algo_t * algo)327 _check_nl_cvg(cs_param_nl_algo_t nl_algo_type,
328 const cs_real_t *pre_iter,
329 cs_real_t *cur_iter,
330 cs_iter_algo_t *algo)
331 {
332 assert(algo != NULL);
333
334 if (nl_algo_type == CS_PARAM_NL_ALGO_ANDERSON && algo->n_algo_iter > 0) {
335
336 /* TODO */
337
338 } /* Anderson acceleration */
339
340 algo->prev_res = algo->res;
341 algo->res = cs_cdo_blas_square_norm_pcsp_diff(pre_iter, cur_iter);
342 assert(algo->res > -DBL_MIN);
343 algo->res = sqrt(algo->res);
344
345 if (algo->n_algo_iter < 1) /* Store the first residual to detect a
346 divergence */
347 algo->res0 = algo->res;
348
349 /* Update the convergence members */
350
351 cs_iter_algo_update_cvg(algo);
352
353 if (algo->param.verbosity > 0) {
354
355 switch (nl_algo_type) {
356
357 case CS_PARAM_NL_ALGO_ANDERSON:
358 if (algo->n_algo_iter == 1)
359 cs_log_printf(CS_LOG_DEFAULT,
360 "### SOLIDIFICATION %12s.It Algo.Res Tolerance\n",
361 "Anderson");
362 cs_log_printf(CS_LOG_DEFAULT,
363 "### SOLIDIFICATION %12s.It%02d %5.3e %6.4e\n",
364 "Anderson", algo->n_algo_iter, algo->res, algo->tol);
365 break;
366
367 case CS_PARAM_NL_ALGO_PICARD:
368 if (algo->n_algo_iter == 1)
369 cs_log_printf(CS_LOG_DEFAULT,
370 "### SOLIDIFICATION %12s.It Algo.Res Tolerance\n",
371 "Picard");
372 cs_log_printf(CS_LOG_DEFAULT,
373 "### SOLIDIFICATION %12s.It%02d %5.3e %6.4e\n",
374 "Picard", algo->n_algo_iter, algo->res, algo->tol);
375 break;
376
377 default:
378 break;
379
380 }
381
382 } /* verbosity > 0 */
383
384 return algo->cvg;
385 }
386
387 /*----------------------------------------------------------------------------*/
388 /*!
389 * \brief Create the structure dedicated to the management of the
390 * solidification module
391 *
392 * \return a pointer to a new allocated cs_solidification_t structure
393 */
394 /*----------------------------------------------------------------------------*/
395
396 static cs_solidification_t *
_solidification_create(void)397 _solidification_create(void)
398 {
399 cs_solidification_t *solid = NULL;
400
401 BFT_MALLOC(solid, 1, cs_solidification_t);
402
403 /* Default initialization */
404
405 solid->model = CS_SOLIDIFICATION_N_MODELS;
406 solid->options = 0;
407 solid->post_flag = 0;
408 solid->verbosity = 1;
409
410 /* Properties */
411
412 solid->mass_density = NULL;
413 solid->viscosity = NULL;
414
415 /* Quantities related to the liquid fraction */
416
417 solid->g_l = NULL;
418 solid->g_l_field = NULL;
419
420 /* State related to each cell */
421
422 solid->cell_state = NULL;
423
424 /* Monitoring */
425
426 for (int i = 0; i < CS_SOLIDIFICATION_N_STATES; i++)
427 solid->n_g_cells[i] = 0;
428 for (int i = 0; i < CS_SOLIDIFICATION_N_STATES; i++)
429 solid->state_ratio[i] = 0;
430
431 /* Plot writer related to the solidification process */
432
433 solid->plot_state = NULL;
434
435 /* Structure related to the thermal system solved as a sub-module */
436
437 solid->temperature = NULL;
438 solid->thermal_reaction_coef = NULL;
439 solid->thermal_reaction_coef_array = NULL;
440 solid->thermal_source_term_array = NULL;
441
442 /* Structure cast on-the-fly w.r.t. the modelling choice */
443
444 solid->model_context = NULL;
445
446 /* Quantities/structure related to the forcing term treated as a reaction term
447 in the momentum equation */
448
449 solid->forcing_mom = NULL;
450 solid->forcing_mom_array = NULL;
451 solid->forcing_coef = 0;
452 solid->first_cell = -1;
453
454 return solid;
455 }
456
457 /*----------------------------------------------------------------------------*/
458 /*!
459 * \brief Compute the number of cells in a given state for monitoring purpose
460 *
461 * \param[in] connect pointer to a cs_cdo_connect_t structure
462 * \param[in] quant pointer to a cs_cdo_quantities_t structure
463 * \param[in,out] solid pointer to the main cs_solidification_t structure
464 */
465 /*----------------------------------------------------------------------------*/
466
467 static void
_monitor_cell_state(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,cs_solidification_t * solid)468 _monitor_cell_state(const cs_cdo_connect_t *connect,
469 const cs_cdo_quantities_t *quant,
470 cs_solidification_t *solid)
471 {
472 for (int i = 0; i < CS_SOLIDIFICATION_N_STATES; i++) solid->n_g_cells[i] = 0;
473
474 for (cs_lnum_t c = 0; c < quant->n_cells; c++) {
475
476 if (connect->cell_flag[c] & CS_FLAG_SOLID_CELL)
477 solid->n_g_cells[CS_SOLIDIFICATION_STATE_SOLID] += 1;
478 else
479 solid->n_g_cells[solid->cell_state[c]] += 1;
480
481 }
482 }
483
484 /*----------------------------------------------------------------------------*/
485 /*!
486 * \brief Perform the monitoring dedicated to the solidification module
487 *
488 * \param[in] quant pointer to a cs_cdo_quantities_t structure
489 */
490 /*----------------------------------------------------------------------------*/
491
492 static void
_do_monitoring(const cs_cdo_quantities_t * quant)493 _do_monitoring(const cs_cdo_quantities_t *quant)
494 {
495 cs_solidification_t *solid = cs_solidification_structure;
496 assert(solid->temperature != NULL);
497
498 for (int i = 0; i < CS_SOLIDIFICATION_N_STATES; i++)
499 solid->state_ratio[i] = 0;
500
501 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
502
503 const cs_real_t vol_c = quant->cell_vol[c_id];
504
505 switch (solid->cell_state[c_id]) {
506 case CS_SOLIDIFICATION_STATE_SOLID:
507 solid->state_ratio[CS_SOLIDIFICATION_STATE_SOLID] += vol_c;
508 break;
509 case CS_SOLIDIFICATION_STATE_LIQUID:
510 solid->state_ratio[CS_SOLIDIFICATION_STATE_LIQUID] += vol_c;
511 break;
512 case CS_SOLIDIFICATION_STATE_MUSHY:
513 solid->state_ratio[CS_SOLIDIFICATION_STATE_MUSHY] += vol_c;
514 break;
515 case CS_SOLIDIFICATION_STATE_EUTECTIC:
516 solid->state_ratio[CS_SOLIDIFICATION_STATE_EUTECTIC] += vol_c;
517 break;
518
519 default: /* Should not be in this case */
520 break;
521
522 } /* End of switch */
523
524 } /* Loop on cells */
525
526 /* Finalize the monitoring step*/
527
528 cs_parall_sum(CS_SOLIDIFICATION_N_STATES, CS_REAL_TYPE, solid->state_ratio);
529 const double inv_voltot = 100./quant->vol_tot;
530 for (int i = 0; i < CS_SOLIDIFICATION_N_STATES; i++)
531 solid->state_ratio[i] *= inv_voltot;
532
533 cs_log_printf(CS_LOG_DEFAULT,
534 "### Solidification monitoring: liquid/mushy/solid states\n"
535 " * Solid | %6.2f\%% for %9lu cells;\n"
536 " * Mushy | %6.2f\%% for %9lu cells;\n"
537 " * Liquid | %6.2f\%% for %9lu cells;\n",
538 solid->state_ratio[CS_SOLIDIFICATION_STATE_SOLID],
539 solid->n_g_cells[CS_SOLIDIFICATION_STATE_SOLID],
540 solid->state_ratio[CS_SOLIDIFICATION_STATE_MUSHY],
541 solid->n_g_cells[CS_SOLIDIFICATION_STATE_MUSHY],
542 solid->state_ratio[CS_SOLIDIFICATION_STATE_LIQUID],
543 solid->n_g_cells[CS_SOLIDIFICATION_STATE_LIQUID]);
544
545 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY)
546 cs_log_printf(CS_LOG_DEFAULT,
547 " * Eutectic | %6.2f\%% for %9lu cells;\n",
548 solid->state_ratio[CS_SOLIDIFICATION_STATE_EUTECTIC],
549 solid->n_g_cells[CS_SOLIDIFICATION_STATE_EUTECTIC]);
550
551 }
552
553 /*----------------------------------------------------------------------------*/
554 /*!
555 * \brief Add a source term to the solute equation derived from an explicit
556 * use of the advective and diffusive operator
557 * Generic function prototype for a hook during the cellwise building
558 * of the linear system
559 * Fit the cs_equation_user_hook_t prototype. This function may be
560 * called by different OpenMP threads
561 *
562 * \param[in] eqp pointer to a cs_equation_param_t structure
563 * \param[in] eqb pointer to a cs_equation_builder_t structure
564 * \param[in] eqc context to cast for this discretization
565 * \param[in] cm pointer to a cellwise view of the mesh
566 * \param[in, out] mass_hodge pointer to a cs_hodge_t structure (mass matrix)
567 * \param[in, out] diff_hodge pointer to a cs_hodge_t structure (diffusion)
568 * \param[in, out] csys pointer to a cellwise view of the system
569 * \param[in, out] cb pointer to a cellwise builder
570 */
571 /*----------------------------------------------------------------------------*/
572
573 static void
_fb_solute_source_term(const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const void * eq_context,const cs_cell_mesh_t * cm,cs_hodge_t * mass_hodge,cs_hodge_t * diff_hodge,cs_cell_sys_t * csys,cs_cell_builder_t * cb)574 _fb_solute_source_term(const cs_equation_param_t *eqp,
575 const cs_equation_builder_t *eqb,
576 const void *eq_context,
577 const cs_cell_mesh_t *cm,
578 cs_hodge_t *mass_hodge,
579 cs_hodge_t *diff_hodge,
580 cs_cell_sys_t *csys,
581 cs_cell_builder_t *cb)
582 {
583 CS_UNUSED(mass_hodge);
584 CS_UNUSED(eqb);
585
586 if (cb->cell_flag & CS_FLAG_SOLID_CELL)
587 return; /* No solute evolution in permanent solid zone */
588
589 const cs_cdofb_scaleq_t *eqc = (const cs_cdofb_scaleq_t *)eq_context;
590
591 cs_solidification_t *solid = cs_solidification_structure;
592 cs_solidification_binary_alloy_t *alloy
593 = (cs_solidification_binary_alloy_t *)solid->model_context;
594
595 cs_real_t *cl_c = alloy->c_l_cells;
596 cs_real_t *cl_f = alloy->c_l_faces;
597
598 /* Diffusion part of the source term to add */
599
600 cs_hodge_set_property_value_cw(cm, cb->t_pty_eval, cb->cell_flag,
601 diff_hodge);
602
603 /* Define the local stiffness matrix: local matrix owned by the cellwise
604 builder (store in cb->loc) */
605
606 eqc->get_stiffness_matrix(cm, diff_hodge, cb);
607
608 /* Build the cellwise array: c - c_l
609 One should have c_l >= c. Therefore, one takes fmin(...,0) */
610
611 for (short int f = 0; f < cm->n_fc; f++)
612 cb->values[f] = fmin(csys->val_n[f] - cl_f[cm->f_ids[f]], 0);
613 cb->values[cm->n_fc] = fmin(csys->val_n[cm->n_fc] - cl_c[cm->c_id], 0);
614
615 /* Update the RHS with the diffusion contribution */
616
617 cs_sdm_update_matvec(cb->loc, cb->values, csys->rhs);
618
619 /* Define the local advection matrix */
620
621 /* Open hook: Compute the advection flux for the numerical scheme and store
622 the advection fluxes across primal faces */
623
624 eqc->advection_open(eqp, cm, csys, eqc->advection_input, cb);
625
626 eqc->advection_main(eqp, cm, csys, eqc->advection_scheme, cb);
627
628 /* Build the cellwise array: c - c_l
629 One should have c_l >= c. Therefore, one takes fmin(...,0) */
630
631 for (short int f = 0; f < cm->n_fc; f++)
632 cb->values[f] = fmin(csys->val_n[f] - cl_f[cm->f_ids[f]], 0);
633 cb->values[cm->n_fc] = fmin(csys->val_n[cm->n_fc] - cl_c[cm->c_id], 0);
634
635 /* Update the RHS with the convection contribution */
636
637 cs_sdm_update_matvec(cb->loc, cb->values, csys->rhs);
638 }
639
640 /*----------------------------------------------------------------------------*/
641 /*!
642 * \brief Build the list of (local) solid cells and enforce a zero-velocity
643 * for this selection
644 *
645 * \param[in] connect pointer to a cs_cdo_connect_t structure
646 * \param[in] quant pointer to a cs_cdo_quantities_t structure
647 */
648 /*----------------------------------------------------------------------------*/
649
650 static void
_enforce_solid_cells(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant)651 _enforce_solid_cells(const cs_cdo_connect_t *connect,
652 const cs_cdo_quantities_t *quant)
653 {
654 cs_solidification_t *solid = cs_solidification_structure;
655
656 cs_gnum_t n_solid_cells = solid->n_g_cells[CS_SOLIDIFICATION_STATE_SOLID];
657
658 /* List of solid cells */
659
660 cs_solid_selection_t *scells = cs_solid_selection_get();
661
662 if (n_solid_cells > (cs_gnum_t)scells->n_cells)
663 BFT_REALLOC(scells->cell_ids, n_solid_cells, cs_lnum_t);
664
665 scells->n_cells = n_solid_cells;
666
667 if (n_solid_cells > 0) {
668
669 n_solid_cells = 0;
670 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
671 if (solid->cell_state[c_id] == CS_SOLIDIFICATION_STATE_SOLID)
672 scells->cell_ids[n_solid_cells++] = c_id;
673 }
674
675 assert(n_solid_cells == solid->n_g_cells[CS_SOLIDIFICATION_STATE_SOLID]);
676
677 }
678
679 /* Parallel synchronization of the number of solid cells. Update the
680 structure storing the list of solid cells and faces */
681
682 cs_solid_selection_sync(connect);
683
684 /* Enforce a zero velocity inside solid cells (enforcement of the momentum
685 equation) */
686
687 cs_navsto_system_set_solid_cells(scells->n_cells, scells->cell_ids);
688 }
689
690 /*----------------------------------------------------------------------------*/
691 /*!
692 * \brief Compute the enthalpy at each cell centers
693 *
694 * \param[in] quant pointer to a cs_cdo_quantities_t structure
695 * \param[in] t_eval physical time at which evaluation is performed
696 * \param[in] temp array of temperature values at each cell
697 * \param[in] g_l array of the liquid fraction values at each cell
698 * \param[in] t_ref reference temperature
699 * \param[in] latent_heat value of the latent heat coefficient
700 * \param[in] rho property related to the mass density
701 * \param[in] cp property related to the heat capacity
702 * \param[in, out] enthalpy array of enthalpy values at each cell
703 */
704 /*----------------------------------------------------------------------------*/
705
706 static void
_compute_enthalpy(const cs_cdo_quantities_t * quant,cs_real_t t_eval,const cs_real_t temp[],const cs_real_t g_l[],const cs_real_t temp_ref,const cs_real_t latent_heat,const cs_property_t * rho,const cs_property_t * cp,cs_real_t enthalpy[])707 _compute_enthalpy(const cs_cdo_quantities_t *quant,
708 cs_real_t t_eval,
709 const cs_real_t temp[],
710 const cs_real_t g_l[],
711 const cs_real_t temp_ref,
712 const cs_real_t latent_heat,
713 const cs_property_t *rho,
714 const cs_property_t *cp,
715 cs_real_t enthalpy[])
716 {
717 assert(temp != NULL && g_l != NULL && enthalpy != NULL);
718
719 if (quant->n_cells < 1)
720 return;
721
722 cs_real_t rho_c, cp_c;
723
724 bool rho_is_uniform = cs_property_is_uniform(rho);
725 bool cp_is_uniform = cs_property_is_uniform(cp);
726
727 /* Use cell with id 0 to evaluate the properties */
728
729 if (rho_is_uniform)
730 rho_c = cs_property_get_cell_value(0, t_eval, rho);
731
732 if (cp_is_uniform)
733 cp_c = cs_property_get_cell_value(0, t_eval, cp);
734
735 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
736 for (cs_lnum_t c = 0; c < quant->n_cells; c++) {
737
738 /* Retrieve the value of the properties if non uniform */
739
740 if (!rho_is_uniform)
741 rho_c = cs_property_get_cell_value(c, t_eval, rho);
742 if (!cp_is_uniform)
743 cp_c = cs_property_get_cell_value(c, t_eval, cp);
744
745 enthalpy[c] = rho_c *
746 /* part linked to the variation of | part linked to the phase change
747 temperature | */
748 ( cp_c * (temp[c] - temp_ref) + latent_heat * g_l[c] );
749
750 } /* Loop on cells */
751 }
752
753 /*----------------------------------------------------------------------------*
754 * Update functions for the Voller & Prakash modelling
755 *----------------------------------------------------------------------------*/
756
757 /*----------------------------------------------------------------------------*/
758 /*!
759 * \brief Update/initialize the liquid fraction, the cell state and the array
760 * used to compute the forcing term in the momentum equation.
761 * This corresponds to the methodology described in the paper
762 * written by Voller and Prakash (87).
763 *
764 * \param[in] mesh pointer to a cs_mesh_t structure
765 * \param[in] connect pointer to a cs_cdo_connect_t structure
766 * \param[in] quant pointer to a cs_cdo_quantities_t structure
767 * \param[in] ts pointer to a cs_time_step_t structure
768 */
769 /*----------------------------------------------------------------------------*/
770
771 static void
_update_gl_voller_legacy(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)772 _update_gl_voller_legacy(const cs_mesh_t *mesh,
773 const cs_cdo_connect_t *connect,
774 const cs_cdo_quantities_t *quant,
775 const cs_time_step_t *ts)
776 {
777 CS_UNUSED(mesh);
778
779 cs_solidification_t *solid = cs_solidification_structure;
780 cs_solidification_voller_t *v_model
781 = (cs_solidification_voller_t *)solid->model_context;
782
783 assert(solid->temperature != NULL);
784 assert(v_model != NULL);
785
786 cs_real_t *g_l = solid->g_l_field->val;
787 cs_real_t *temp = solid->temperature->val;
788 assert(temp != NULL);
789
790 /* 1./(t_liquidus - t_solidus) = \partial g_l/\partial Temp */
791
792 const cs_real_t dgldT = 1./(v_model->t_liquidus - v_model->t_solidus);
793 const cs_real_t inv_forcing_eps = 1./cs_solidification_forcing_eps;
794
795 assert(cs_property_is_uniform(solid->viscosity));
796 const cs_real_t viscl0 = cs_property_get_cell_value(solid->first_cell,
797 ts->t_cur,
798 solid->viscosity);
799 const cs_real_t forcing_coef = solid->forcing_coef * viscl0;
800
801 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
802
803 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL) {
804
805 g_l[c_id] = 0;
806 solid->cell_state[c_id] = CS_SOLIDIFICATION_STATE_SOLID;
807
808 }
809
810 /* Update the liquid fraction */
811
812 else if (temp[c_id] < v_model->t_solidus) {
813
814 g_l[c_id] = 0;
815 solid->cell_state[c_id] = CS_SOLIDIFICATION_STATE_SOLID;
816
817 /* Update the forcing coefficient treated as a property for a reaction
818 term in the momentum eq. */
819
820 solid->forcing_mom_array[c_id] = forcing_coef*inv_forcing_eps;
821
822 }
823 else if (temp[c_id] > v_model->t_liquidus) {
824
825 g_l[c_id] = 1;
826 solid->cell_state[c_id] = CS_SOLIDIFICATION_STATE_LIQUID;
827
828 /* Update the forcing coefficient treated as a property for a reaction
829 term in the momentum eq. */
830
831 solid->forcing_mom_array[c_id] = 0;
832
833 }
834 else { /* Mushy zone */
835
836 const cs_real_t glc = (temp[c_id] - v_model->t_solidus) * dgldT;
837
838 g_l[c_id] = glc;
839 solid->cell_state[c_id] = CS_SOLIDIFICATION_STATE_MUSHY;
840
841 /* Update the forcing coefficient treated as a property for a reaction
842 term in the momentum eq. */
843
844 const cs_real_t glm1 = 1 - glc;
845 solid->forcing_mom_array[c_id] =
846 forcing_coef * glm1*glm1/(glc*glc*glc + cs_solidification_forcing_eps);
847
848 }
849
850 } /* Loop on cells */
851 }
852
853 /*----------------------------------------------------------------------------*/
854 /*!
855 * \brief Update/initialize the liquid fraction and the cell state. The array
856 * used to compute the forcing term in the momentum equation is not
857 * considered in the following function since there is no velocity
858 * field. This corresponds to the methodology described in the paper
859 * written by Voller and Prakash (87).
860 *
861 * \param[in] mesh pointer to a cs_mesh_t structure
862 * \param[in] connect pointer to a cs_cdo_connect_t structure
863 * \param[in] quant pointer to a cs_cdo_quantities_t structure
864 * \param[in] ts pointer to a cs_time_step_t structure
865 */
866 /*----------------------------------------------------------------------------*/
867
868 static void
_update_gl_voller_legacy_no_velocity(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)869 _update_gl_voller_legacy_no_velocity(const cs_mesh_t *mesh,
870 const cs_cdo_connect_t *connect,
871 const cs_cdo_quantities_t *quant,
872 const cs_time_step_t *ts)
873 {
874 CS_UNUSED(mesh);
875 CS_UNUSED(ts);
876
877 cs_solidification_t *solid = cs_solidification_structure;
878 cs_solidification_voller_t *v_model =
879 (cs_solidification_voller_t *)solid->model_context;
880
881 assert(solid->temperature != NULL);
882 assert(v_model != NULL);
883
884 cs_real_t *g_l = solid->g_l_field->val;
885 cs_real_t *temp = solid->temperature->val;
886 assert(temp != NULL);
887
888 /* 1./(t_liquidus - t_solidus) = \partial g_l/\partial Temp */
889
890 const cs_real_t dgldT = 1./(v_model->t_liquidus - v_model->t_solidus);
891
892 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
893
894 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL) {
895
896 g_l[c_id] = 0;
897 solid->cell_state[c_id] = CS_SOLIDIFICATION_STATE_SOLID;
898
899 }
900
901 /* Update the liquid fraction */
902
903 else if (temp[c_id] < v_model->t_solidus) {
904
905 g_l[c_id] = 0;
906 solid->cell_state[c_id] = CS_SOLIDIFICATION_STATE_SOLID;
907
908 }
909 else if (temp[c_id] > v_model->t_liquidus) {
910
911 g_l[c_id] = 1;
912 solid->cell_state[c_id] = CS_SOLIDIFICATION_STATE_LIQUID;
913
914 }
915 else { /* Mushy zone */
916
917 const cs_real_t glc = (temp[c_id] - v_model->t_solidus) * dgldT;
918
919 g_l[c_id] = glc;
920 solid->cell_state[c_id] = CS_SOLIDIFICATION_STATE_MUSHY;
921
922 }
923
924 } /* Loop on cells */
925 }
926
927 /*----------------------------------------------------------------------------*/
928 /*!
929 * \brief Update/initialize the reaction and source term for the thermal
930 * equation. This corresponds to the methodology described in the paper
931 * written by Voller and Prakash (87)
932 *
933 * \param[in] mesh pointer to a cs_mesh_t structure
934 * \param[in] connect pointer to a cs_cdo_connect_t structure
935 * \param[in] quant pointer to a cs_cdo_quantities_t structure
936 * \param[in] ts pointer to a cs_time_step_t structure
937 */
938 /*----------------------------------------------------------------------------*/
939
940 static void
_update_thm_voller_legacy(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)941 _update_thm_voller_legacy(const cs_mesh_t *mesh,
942 const cs_cdo_connect_t *connect,
943 const cs_cdo_quantities_t *quant,
944 const cs_time_step_t *ts)
945 {
946 CS_UNUSED(mesh);
947 CS_UNUSED(connect);
948
949 cs_solidification_t *solid = cs_solidification_structure;
950 cs_solidification_voller_t *v_model
951 = (cs_solidification_voller_t *)solid->model_context;
952
953 assert(v_model != NULL);
954 assert(solid->temperature != NULL);
955
956 const cs_real_t *temp = solid->temperature->val;
957 assert(temp != NULL);
958
959 /* 1./(t_liquidus - t_solidus) = \partial g_l/\partial Temp */
960
961 const cs_real_t rho0 = solid->mass_density->ref_value;
962 const cs_real_t dgldT = 1./(v_model->t_liquidus - v_model->t_solidus);
963 const cs_real_t dgldT_coef = rho0*solid->latent_heat*dgldT/ts->dt[0];
964
965 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
966
967 if (solid->cell_state[c_id] == CS_SOLIDIFICATION_STATE_MUSHY) {
968
969 solid->thermal_reaction_coef_array[c_id] = dgldT_coef;
970 solid->thermal_source_term_array[c_id] =
971 dgldT_coef*temp[c_id]*quant->cell_vol[c_id];
972
973 }
974 else {
975
976 solid->thermal_reaction_coef_array[c_id] = 0;
977 solid->thermal_source_term_array[c_id] = 0;
978
979 }
980
981 } /* Loop on cells */
982 }
983
984 /*----------------------------------------------------------------------------*/
985 /*!
986 * \brief Update/initialize the reaction and source term for the thermal
987 * equation. One considers the state at the previous time step and that
988 * the kth sub-iteration to determine the solidification path and
989 * compute the related quantities.
990 *
991 * \param[in] mesh pointer to a cs_mesh_t structure
992 * \param[in] connect pointer to a cs_cdo_connect_t structure
993 * \param[in] quant pointer to a cs_cdo_quantities_t structure
994 * \param[in] ts pointer to a cs_time_step_t structure
995 */
996 /*----------------------------------------------------------------------------*/
997
998 static void
_update_thm_voller_path(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)999 _update_thm_voller_path(const cs_mesh_t *mesh,
1000 const cs_cdo_connect_t *connect,
1001 const cs_cdo_quantities_t *quant,
1002 const cs_time_step_t *ts)
1003 {
1004 CS_UNUSED(mesh);
1005 CS_UNUSED(connect);
1006
1007 cs_solidification_t *solid = cs_solidification_structure;
1008 cs_solidification_voller_t *v_model
1009 = (cs_solidification_voller_t *)solid->model_context;
1010
1011 /* Sanity checks */
1012
1013 assert(v_model != NULL);
1014 assert(solid->temperature != NULL);
1015
1016 const cs_real_t *temp = solid->temperature->val;
1017 const cs_real_t *temp_pre = solid->temperature->val_pre;
1018 assert(temp != NULL && temp_pre != NULL);
1019
1020 /* 1./(t_liquidus - t_solidus) = \partial g_l/\partial Temp */
1021
1022 const cs_real_t dgldT = 1./(v_model->t_liquidus - v_model->t_solidus);
1023 const cs_real_t coef =
1024 solid->mass_density->ref_value*solid->latent_heat/ts->dt[0];
1025 const cs_real_t dgldT_coef = coef * dgldT;
1026
1027 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1028
1029 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL) {
1030
1031 /* Keep the solid state during all the computation */
1032
1033 solid->thermal_reaction_coef_array[c_id] = 0;
1034 solid->thermal_source_term_array[c_id] = 0;
1035
1036 }
1037 else if (temp[c_id] < v_model->t_solidus) {
1038
1039 if (temp_pre[c_id] > v_model->t_liquidus) {
1040
1041 /* Liquid --> solid state */
1042
1043 solid->thermal_reaction_coef_array[c_id] = dgldT_coef;
1044 solid->thermal_source_term_array[c_id] =
1045 dgldT_coef*v_model->t_liquidus*quant->cell_vol[c_id];
1046
1047 }
1048 else if (temp_pre[c_id] < v_model->t_solidus) {
1049
1050 /* Solid --> Solid state */
1051
1052 solid->thermal_reaction_coef_array[c_id] = 0;
1053 solid->thermal_source_term_array[c_id] = 0;
1054
1055 }
1056 else { /* Mushy --> solid state */
1057
1058 solid->thermal_reaction_coef_array[c_id] = dgldT_coef;
1059 solid->thermal_source_term_array[c_id] =
1060 dgldT_coef*temp_pre[c_id]*quant->cell_vol[c_id];
1061
1062 /* Strictly speaking this should not be divided by 1/dt but with a
1063 smaller time step (Tsolidus is reached before the end of the time
1064 step) */
1065
1066 }
1067
1068 }
1069 else if (temp[c_id] > v_model->t_liquidus) {
1070
1071 if (temp_pre[c_id] > v_model->t_liquidus) {
1072
1073 /* Liquid --> liquid state */
1074
1075 solid->thermal_reaction_coef_array[c_id] = 0;
1076 solid->thermal_source_term_array[c_id] = 0;
1077
1078 }
1079 else if (temp_pre[c_id] < v_model->t_solidus) {
1080
1081 /* Solid --> liquid state */
1082
1083 solid->thermal_reaction_coef_array[c_id] = dgldT_coef;
1084 solid->thermal_source_term_array[c_id] =
1085 dgldT_coef*v_model->t_solidus*quant->cell_vol[c_id];
1086
1087 }
1088 else { /* Mushy --> liquid state */
1089
1090 solid->thermal_reaction_coef_array[c_id] = dgldT_coef;
1091 solid->thermal_source_term_array[c_id] =
1092 dgldT_coef*temp_pre[c_id]*quant->cell_vol[c_id];
1093
1094 }
1095
1096 }
1097 else {
1098
1099 if (temp_pre[c_id] > v_model->t_liquidus) {
1100
1101 /* Liquid --> mushy state */
1102
1103 solid->thermal_reaction_coef_array[c_id] = dgldT_coef;
1104 solid->thermal_source_term_array[c_id] =
1105 dgldT_coef*v_model->t_liquidus*quant->cell_vol[c_id];
1106
1107 }
1108 else if (temp_pre[c_id] < v_model->t_solidus) {
1109
1110 /* Solid --> mushy state */
1111
1112 solid->thermal_reaction_coef_array[c_id] = dgldT_coef;
1113 solid->thermal_source_term_array[c_id] =
1114 dgldT_coef*v_model->t_solidus*quant->cell_vol[c_id];
1115
1116 }
1117 else { /* Mushy --> mushy state */
1118
1119 solid->thermal_reaction_coef_array[c_id] = dgldT_coef;
1120 solid->thermal_source_term_array[c_id] =
1121 dgldT_coef*temp_pre[c_id]*quant->cell_vol[c_id];
1122
1123 } /* State for the previous temp (n-1) */
1124
1125 } /* State for the current temp (n+1,k+1) */
1126
1127 } /* Loop on cells */
1128
1129 }
1130
1131 /*----------------------------------------------------------------------------*
1132 * Update functions for the binary alloy modelling
1133 *----------------------------------------------------------------------------*/
1134
1135 /*----------------------------------------------------------------------------*/
1136 /*!
1137 * \brief Update the state associated to each cell in the case of a binary
1138 * alloy. No MPI synchronization has to be performed at this stage.
1139 *
1140 * \param[in] connect pointer to a cs_cdo_connect_t structure
1141 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1142 * \param[in] ts pointer to a cs_time_step_t structure
1143 */
1144 /*----------------------------------------------------------------------------*/
1145
1146 static void
_update_binary_alloy_final_state(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)1147 _update_binary_alloy_final_state(const cs_cdo_connect_t *connect,
1148 const cs_cdo_quantities_t *quant,
1149 const cs_time_step_t *ts)
1150 {
1151 CS_UNUSED(ts);
1152
1153 cs_solidification_t *solid = cs_solidification_structure;
1154 cs_solidification_binary_alloy_t *alloy
1155 = (cs_solidification_binary_alloy_t *)solid->model_context;
1156
1157 /* Update the cell state (at this stage, one should have converged between
1158 * the couple (temp, conc) and the liquid fraction */
1159 const cs_real_t *t_bulk = solid->temperature->val;
1160 const cs_real_t *c_bulk = alloy->c_bulk->val;
1161 const cs_real_t *g_l = solid->g_l_field->val;
1162
1163 for (int i = 0; i < CS_SOLIDIFICATION_N_STATES; i++) solid->n_g_cells[i] = 0;
1164
1165 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1166
1167 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL) {
1168
1169 solid->cell_state[c_id] = CS_SOLIDIFICATION_STATE_SOLID;
1170 solid->n_g_cells[CS_SOLIDIFICATION_STATE_SOLID] += 1;
1171
1172 }
1173 else {
1174
1175 cs_solidification_state_t
1176 state = _which_state_by_enthalpy(alloy,
1177 solid->latent_heat,
1178 solid->cp->ref_value,
1179 t_bulk[c_id],
1180 c_bulk[c_id],
1181 g_l[c_id]);
1182
1183 solid->cell_state[c_id] = state;
1184 solid->n_g_cells[state] += 1;
1185
1186 }
1187
1188 } /* Loop on cells */
1189
1190 }
1191
1192 /*----------------------------------------------------------------------------*/
1193 /*!
1194 * \brief Update the Darcy term (acting as a penalization) in the momentum
1195 * equation and enforce solid cells by setting a zero mass flux.
1196 * The parallel reduction on the cell state is performed here (not
1197 * before to avoid calling the enforcement if no solid cell is locally
1198 * detected).
1199 *
1200 * \param[in] mesh pointer to a cs_mesh_t structure
1201 * \param[in] connect pointer to a cs_cdo_connect_t structure
1202 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1203 * \param[in] ts pointer to a cs_time_step_t structure
1204 */
1205 /*----------------------------------------------------------------------------*/
1206
1207 static void
_update_velocity_forcing(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)1208 _update_velocity_forcing(const cs_mesh_t *mesh,
1209 const cs_cdo_connect_t *connect,
1210 const cs_cdo_quantities_t *quant,
1211 const cs_time_step_t *ts)
1212 {
1213 CS_UNUSED(mesh);
1214
1215 cs_solidification_t *solid = cs_solidification_structure;
1216
1217 /* At this stage, the number of solid cells is a local count
1218 * Set the enforcement of the velocity for solid cells */
1219
1220 _enforce_solid_cells(connect, quant);
1221
1222 /* Parallel synchronization of the number of cells in each state
1223 * This should be done done now to avoid going to the cell enforcement whereas
1224 * there is nothing to do locally */
1225
1226 cs_parall_sum(CS_SOLIDIFICATION_N_STATES, CS_GNUM_TYPE, solid->n_g_cells);
1227
1228 assert(cs_property_is_uniform(solid->viscosity));
1229 const cs_real_t viscl0 = cs_property_get_cell_value(0, ts->t_cur,
1230 solid->viscosity);
1231 const cs_real_t forcing_coef = solid->forcing_coef * viscl0;
1232 const cs_real_t *g_l = solid->g_l_field->val;
1233
1234 /* Set the forcing term in the momentum equation */
1235
1236 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1237
1238 if (g_l[c_id] < 1.) { /* Not fully liquid */
1239
1240 const cs_real_t gsc = 1 - g_l[c_id];
1241 const cs_real_t glc3 = g_l[c_id]*g_l[c_id]*g_l[c_id];
1242
1243 solid->forcing_mom_array[c_id] =
1244 forcing_coef * gsc*gsc/(glc3 + cs_solidification_forcing_eps);
1245
1246 }
1247 else
1248 solid->forcing_mom_array[c_id] = 0;
1249
1250 } /* Loop on cells */
1251
1252 }
1253
1254 /*----------------------------------------------------------------------------*/
1255 /*!
1256 * \brief Update the concentration of solute in the liquid phase at the cell
1257 * center. This value is used in the buoyancy term in the momentum
1258 * equation.
1259 *
1260 * \param[in] mesh pointer to a cs_mesh_t structure
1261 * \param[in] connect pointer to a cs_cdo_connect_t structure
1262 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1263 * \param[in] ts pointer to a cs_time_step_t structure
1264 */
1265 /*----------------------------------------------------------------------------*/
1266
1267 static void
_update_clc(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)1268 _update_clc(const cs_mesh_t *mesh,
1269 const cs_cdo_connect_t *connect,
1270 const cs_cdo_quantities_t *quant,
1271 const cs_time_step_t *ts)
1272 {
1273 CS_UNUSED(mesh);
1274 CS_UNUSED(ts);
1275
1276 cs_solidification_t *solid = cs_solidification_structure;
1277 cs_solidification_binary_alloy_t *alloy
1278 = (cs_solidification_binary_alloy_t *)solid->model_context;
1279
1280 const cs_real_t *c_bulk = alloy->c_bulk->val;
1281 const cs_real_t *t_bulk = solid->temperature->val;
1282 const cs_real_t *g_l_pre = solid->g_l_field->val_pre;
1283
1284 cs_real_t *c_l = alloy->c_l_cells;
1285
1286 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1287
1288 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL) {
1289 c_l[c_id] = 0.;
1290 continue;
1291 }
1292
1293 const cs_real_t conc = c_bulk[c_id];
1294 const cs_real_t temp = t_bulk[c_id];
1295
1296 switch (_which_state(alloy, temp, conc)) {
1297
1298 case CS_SOLIDIFICATION_STATE_SOLID:
1299 /* If this is the first time that one reaches the solid state for this
1300 * cell (i.e previously with g_l > 0), then one updates the liquid
1301 * concentration and one keeps that value */
1302 if (g_l_pre[c_id] > 0) {
1303 if (conc < alloy->cs1)
1304 c_l[c_id] = conc * alloy->inv_kp;
1305 else
1306 c_l[c_id] = alloy->c_eut;
1307 }
1308 break;
1309
1310 case CS_SOLIDIFICATION_STATE_MUSHY:
1311 c_l[c_id] = (temp - alloy->t_melt) * alloy->inv_ml;
1312 break;
1313
1314 case CS_SOLIDIFICATION_STATE_LIQUID:
1315 c_l[c_id] = conc;
1316 break;
1317
1318 case CS_SOLIDIFICATION_STATE_EUTECTIC:
1319 c_l[c_id] = alloy->c_eut;
1320 break;
1321
1322 default:
1323 bft_error(__FILE__, __LINE__, 0,
1324 " %s: Invalid state for cell %ld\n", __func__, (long)c_id);
1325 break;
1326
1327 } /* Switch on cell state */
1328
1329 } /* Loop on cells */
1330
1331 }
1332
1333 /*----------------------------------------------------------------------------*/
1334 /*!
1335 * \brief Update the liquid fraction in each cell
1336 * This function reproduces the same process as the one used in the
1337 * legacy FV scheme.
1338 * This corresponds to the case of a binary alloy model with no
1339 * advective source term for the solute transport.
1340 *
1341 * \param[in] mesh pointer to a cs_mesh_t structure
1342 * \param[in] connect pointer to a cs_cdo_connect_t structure
1343 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1344 * \param[in] ts pointer to a cs_time_step_t structure
1345 */
1346 /*----------------------------------------------------------------------------*/
1347
1348 static void
_update_gl_legacy(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)1349 _update_gl_legacy(const cs_mesh_t *mesh,
1350 const cs_cdo_connect_t *connect,
1351 const cs_cdo_quantities_t *quant,
1352 const cs_time_step_t *ts)
1353 {
1354 CS_UNUSED(mesh);
1355 CS_UNUSED(ts);
1356
1357 cs_solidification_t *solid = cs_solidification_structure;
1358 cs_solidification_binary_alloy_t *alloy
1359 = (cs_solidification_binary_alloy_t *)solid->model_context;
1360
1361 const cs_real_t *c_bulk = alloy->c_bulk->val;
1362 const cs_real_t *t_bulk = solid->temperature->val;
1363 const cs_real_t *g_l_pre = solid->g_l_field->val_pre;
1364 cs_real_t *g_l = solid->g_l_field->val;
1365
1366 /* Update g_l values in each cell as well as the cell state and the related
1367 count */
1368
1369 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1370
1371 cs_real_t eta_new, gliq;
1372
1373 const cs_real_t eta_old = alloy->eta_coef_array[c_id];
1374 const cs_real_t conc = c_bulk[c_id];
1375 const cs_real_t temp = t_bulk[c_id];
1376
1377 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL)
1378 continue; /* No update */
1379
1380 /* Knowing in which part of the phase diagram we are and we then update the
1381 * value of the liquid fraction: g_l and eta (the coefficient between the
1382 * concentration in the liquid phase and the bulk concentration */
1383
1384 switch (_which_state(alloy, temp, conc)) {
1385
1386 case CS_SOLIDIFICATION_STATE_SOLID:
1387 gliq = 0.;
1388 if (g_l_pre[c_id] > 0) /* Not in a solid state */
1389 eta_new = _get_eta(alloy, conc);
1390 else
1391 eta_new = eta_old;
1392 break;
1393
1394 case CS_SOLIDIFICATION_STATE_MUSHY:
1395 gliq = alloy->inv_kpm1* (alloy->kp - alloy->ml*conc/(temp-alloy->t_melt));
1396
1397 /* Make sure that the liquid fraction remains inside physical bounds */
1398
1399 gliq = fmin(fmax(0, gliq), 1.);
1400
1401 eta_new = 1/( gliq * (1-alloy->kp) + alloy->kp );
1402 break;
1403
1404 case CS_SOLIDIFICATION_STATE_LIQUID:
1405 gliq = 1;
1406 eta_new = 1;
1407 break;
1408
1409 case CS_SOLIDIFICATION_STATE_EUTECTIC:
1410 gliq = (conc - alloy->cs1)*alloy->dgldC_eut;
1411
1412 /* Make sure that the liquid fraction remains inside physical bounds */
1413
1414 gliq = fmin(fmax(0, gliq), 1.);
1415
1416 eta_new = _get_eta(alloy, conc);
1417 break;
1418
1419 default:
1420 bft_error(__FILE__, __LINE__, 0,
1421 " %s: Invalid state for cell %ld\n", __func__, (long)c_id);
1422 break;
1423
1424 } /* Switch on cell state */
1425
1426 /* Update the liquid fraction and apply if needed a relaxation */
1427
1428 if (alloy->gliq_relax > 0)
1429 g_l[c_id] = (1 - alloy->gliq_relax)*gliq + alloy->gliq_relax*g_l[c_id];
1430 else
1431 g_l[c_id] = gliq;
1432
1433 /* Update eta and apply if needed a relaxation */
1434
1435 if (alloy->eta_relax > 0)
1436 alloy->eta_coef_array[c_id] =
1437 (1-alloy->eta_relax)*eta_new + alloy->eta_relax*eta_old;
1438 else
1439 alloy->eta_coef_array[c_id] = eta_new;
1440
1441 } /* Loop on cells */
1442
1443 }
1444
1445 /*----------------------------------------------------------------------------*/
1446 /*!
1447 * \brief Update the liquid fraction in each cell
1448 * This function reproduces the same process as the one used in the
1449 * legacy FV scheme.
1450 * This corresponds to the case of a binary alloy model with an
1451 * advective source term for the solute transport.
1452 *
1453 * \param[in] mesh pointer to a cs_mesh_t structure
1454 * \param[in] connect pointer to a cs_cdo_connect_t structure
1455 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1456 * \param[in] ts pointer to a cs_time_step_t structure
1457 */
1458 /*----------------------------------------------------------------------------*/
1459
1460 static void
_update_gl_legacy_ast(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)1461 _update_gl_legacy_ast(const cs_mesh_t *mesh,
1462 const cs_cdo_connect_t *connect,
1463 const cs_cdo_quantities_t *quant,
1464 const cs_time_step_t *ts)
1465 {
1466 CS_UNUSED(mesh);
1467 CS_UNUSED(ts);
1468
1469 cs_solidification_t *solid = cs_solidification_structure;
1470 cs_solidification_binary_alloy_t *alloy
1471 = (cs_solidification_binary_alloy_t *)solid->model_context;
1472
1473 const cs_real_t *c_bulk = alloy->c_bulk->val;
1474 const cs_real_t *t_bulk = solid->temperature->val;
1475 const cs_real_t *g_l_pre = solid->g_l_field->val_pre;
1476 cs_real_t *g_l = solid->g_l_field->val;
1477 cs_real_t *c_l = alloy->c_l_cells;
1478
1479 /* Update g_l values in each cell as well as the cell state and the related
1480 count */
1481
1482 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1483
1484 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL)
1485 continue; /* No update */
1486
1487 cs_real_t gliq = 1; /* Initialization as liquid */
1488
1489 const cs_real_t conc = c_bulk[c_id];
1490 const cs_real_t temp = t_bulk[c_id];
1491
1492 /* Knowing in which part of the phase diagram we are and we then update
1493 * the value of the liquid fraction: g_l and the concentration of the
1494 * liquid "solute" c_l */
1495
1496 switch (_which_state(alloy, temp, conc)) {
1497
1498 case CS_SOLIDIFICATION_STATE_SOLID:
1499 gliq = 0.;
1500
1501 /* If this is the first time that one reaches the solid state for this
1502 * cell (i.e previously with g_l > 0), then one updates the liquid
1503 * concentration and one keeps that value */
1504
1505 if (g_l_pre[c_id] > 0) {
1506 if (conc < alloy->cs1)
1507 c_l[c_id] = conc * alloy->inv_kp;
1508 else
1509 c_l[c_id] = alloy->c_eut;
1510 }
1511 break;
1512
1513 case CS_SOLIDIFICATION_STATE_MUSHY:
1514 gliq = alloy->inv_kpm1* (alloy->kp - alloy->ml*conc/(temp-alloy->t_melt));
1515 c_l[c_id] = (temp - alloy->t_melt) * alloy->inv_ml;
1516 break;
1517
1518 case CS_SOLIDIFICATION_STATE_LIQUID:
1519 c_l[c_id] = conc;
1520 break;
1521
1522 case CS_SOLIDIFICATION_STATE_EUTECTIC:
1523 gliq = (conc - alloy->cs1)*alloy->dgldC_eut;
1524 c_l[c_id] = alloy->c_eut;
1525 break;
1526
1527 default:
1528 bft_error(__FILE__, __LINE__, 0,
1529 " %s: Invalid state for cell %ld\n", __func__, (long)c_id);
1530 break;
1531
1532 } /* Switch on cell state */
1533
1534 /* Make sure that the liquid fraction remains inside physical bounds */
1535
1536 gliq = fmin(fmax(0, gliq), 1.);
1537
1538 /* Relaxation if needed for the liquid fraction */
1539
1540 if (alloy->gliq_relax > 0)
1541 g_l[c_id] = (1 - alloy->gliq_relax)*gliq + alloy->gliq_relax*g_l[c_id];
1542 else
1543 g_l[c_id] = gliq;
1544
1545 } /* Loop on cells */
1546
1547 /* Update c_l at face values */
1548
1549 const cs_equation_t *tr_eq = alloy->solute_equation;
1550 const cs_real_t *c_bulk_f = cs_equation_get_face_values(tr_eq, false);
1551 const cs_real_t *t_bulk_f = alloy->temp_faces;
1552
1553 for (cs_lnum_t f_id = 0; f_id < quant->n_faces; f_id++) {
1554
1555 const cs_real_t conc = c_bulk_f[f_id];
1556 const cs_real_t temp = t_bulk_f[f_id];
1557
1558 /* Knowing in which part of the phase diagram we are, we then update
1559 * the value of the concentration of the liquid "solute" */
1560
1561 switch (_which_state(alloy, temp, conc)) {
1562
1563 case CS_SOLIDIFICATION_STATE_SOLID:
1564 if (conc < alloy->cs1)
1565 alloy->c_l_faces[f_id] = conc * alloy->inv_kp;
1566 else
1567 alloy->c_l_faces[f_id] = alloy->c_eut;
1568 break;
1569
1570 case CS_SOLIDIFICATION_STATE_MUSHY:
1571 alloy->c_l_faces[f_id] = (temp - alloy->t_melt) * alloy->inv_ml;
1572 break;
1573
1574 case CS_SOLIDIFICATION_STATE_LIQUID:
1575 alloy->c_l_faces[f_id] = conc;
1576 break;
1577
1578 case CS_SOLIDIFICATION_STATE_EUTECTIC:
1579 alloy->c_l_faces[f_id] = alloy->c_eut;
1580 break;
1581
1582 default:
1583 bft_error(__FILE__, __LINE__, 0,
1584 " %s: Invalid state for face %ld\n", __func__, (long)f_id);
1585 break;
1586
1587 } /* Switch on face state */
1588
1589 } /* Loop on faces */
1590
1591 }
1592
1593 /*----------------------------------------------------------------------------*/
1594 /*!
1595 * \brief Update the source term for the thermal equation.
1596 * This function reproduces the same process as the one used in the
1597 * legacy FV scheme. This corresponds to the binary alloy model.
1598 *
1599 * \param[in] mesh pointer to a cs_mesh_t structure
1600 * \param[in] connect pointer to a cs_cdo_connect_t structure
1601 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1602 * \param[in] ts pointer to a cs_time_step_t structure
1603 */
1604 /*----------------------------------------------------------------------------*/
1605
1606 static void
_update_thm_legacy(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)1607 _update_thm_legacy(const cs_mesh_t *mesh,
1608 const cs_cdo_connect_t *connect,
1609 const cs_cdo_quantities_t *quant,
1610 const cs_time_step_t *ts)
1611 {
1612 CS_UNUSED(mesh);
1613 cs_solidification_t *solid = cs_solidification_structure;
1614 cs_solidification_binary_alloy_t *alloy
1615 = (cs_solidification_binary_alloy_t *)solid->model_context;
1616
1617 const cs_real_t *c_bulk = alloy->c_bulk->val;
1618 const cs_real_t *c_bulk_pre = alloy->c_bulk->val_pre;
1619 const cs_real_t *t_bulk_pre = solid->temperature->val_pre;
1620
1621 const cs_real_t rhoL = solid->mass_density->ref_value * solid->latent_heat;
1622 const cs_real_t rhoLovdt = rhoL/ts->dt[0];
1623
1624 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1625
1626 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL)
1627 continue; /* No update: 0 by default */
1628
1629 const cs_real_t conc = c_bulk[c_id];
1630 const cs_real_t conc_pre = c_bulk_pre[c_id];
1631 const cs_real_t temp_pre = t_bulk_pre[c_id];
1632
1633 /* Knowing in which part of the phase diagram we are, then we update
1634 * the value of the concentration of the liquid "solute" */
1635
1636 switch (_which_state(alloy, temp_pre, conc_pre)) {
1637
1638 case CS_SOLIDIFICATION_STATE_SOLID:
1639 case CS_SOLIDIFICATION_STATE_LIQUID:
1640 solid->thermal_reaction_coef_array[c_id] = 0;
1641 solid->thermal_source_term_array[c_id] = 0;
1642 break;
1643
1644 case CS_SOLIDIFICATION_STATE_MUSHY:
1645 {
1646 cs_real_t dgldC, dgldT;
1647 _get_dgl_mushy(alloy, temp_pre, conc_pre, &dgldT, &dgldC);
1648
1649 solid->thermal_reaction_coef_array[c_id] = dgldT * rhoLovdt;
1650 solid->thermal_source_term_array[c_id] =
1651 quant->cell_vol[c_id] * rhoLovdt * ( dgldT * temp_pre +
1652 dgldC * (conc_pre - conc) );
1653 }
1654 break;
1655
1656 case CS_SOLIDIFICATION_STATE_EUTECTIC:
1657 solid->thermal_reaction_coef_array[c_id] = 0;
1658 solid->thermal_source_term_array[c_id] = quant->cell_vol[c_id] *
1659 rhoLovdt * alloy->dgldC_eut * (conc_pre - conc);
1660 break;
1661
1662 default:
1663 bft_error(__FILE__, __LINE__, 0,
1664 " %s: Invalid state for cell %ld\n", __func__, (long)c_id);
1665 break;
1666
1667 } /* Switch on cell state */
1668
1669 } /* Loop on cells */
1670
1671 }
1672
1673 /*----------------------------------------------------------------------------*/
1674 /*!
1675 * \brief Update the liquid fraction in each cell and related quantities.
1676 * This corresponds to the case of a binary alloy model with no
1677 * advective source term for the solute transport.
1678 *
1679 * \param[in] mesh pointer to a cs_mesh_t structure
1680 * \param[in] connect pointer to a cs_cdo_connect_t structure
1681 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1682 * \param[in] ts pointer to a cs_time_step_t structure
1683 */
1684 /*----------------------------------------------------------------------------*/
1685
1686 static void
_update_gl_taylor(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)1687 _update_gl_taylor(const cs_mesh_t *mesh,
1688 const cs_cdo_connect_t *connect,
1689 const cs_cdo_quantities_t *quant,
1690 const cs_time_step_t *ts)
1691 {
1692 CS_UNUSED(mesh);
1693 CS_UNUSED(ts);
1694 cs_solidification_t *solid = cs_solidification_structure;
1695 cs_solidification_binary_alloy_t *alloy
1696 = (cs_solidification_binary_alloy_t *)solid->model_context;
1697
1698 const double cpovL = solid->cp->ref_value/solid->latent_heat;
1699
1700 const cs_real_t *c_bulk = alloy->c_bulk->val;
1701 const cs_real_t *c_bulk_pre = alloy->c_bulk->val_pre;
1702 const cs_real_t *t_bulk_pre = solid->temperature->val_pre;
1703 cs_real_t *t_bulk = solid->temperature->val;
1704 cs_real_t *g_l = solid->g_l_field->val;
1705
1706 /* Update g_l values in each cell as well as the cell state and the related
1707 count */
1708
1709 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1710
1711 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL)
1712 continue; /* No update */
1713
1714 const cs_real_t conc = c_bulk[c_id]; /* conc_{n+1}^{k+1} */
1715 const cs_real_t temp = t_bulk[c_id]; /* temp_{n+1}^{k+1} */
1716 const cs_real_t conc_pre = c_bulk_pre[c_id];
1717 const cs_real_t temp_pre = t_bulk_pre[c_id];
1718
1719 cs_real_t dgldC, dgldT, gliq, eta_new;
1720 cs_solidification_state_t state, state_pre;
1721
1722 /* gliq, temp and conc iterates may not be related with the gliq(temp, conc)
1723 * function until convergence is reached. So one needs to be careful. */
1724
1725 state = _which_state(alloy, temp, conc);
1726 state_pre = _which_state(alloy, temp_pre, conc_pre);
1727 eta_new = alloy->eta_coef_array[c_id]; /* avoid a warning */
1728 gliq = g_l[c_id]; /* avoid a warning */
1729
1730 /* Knowing in which part of the phase diagram we are and we then update
1731 * the value of the liquid fraction: g_l and the concentration of the
1732 * liquid "solute" */
1733
1734 switch (state) {
1735
1736 case CS_SOLIDIFICATION_STATE_SOLID:
1737
1738 if (state_pre == CS_SOLIDIFICATION_STATE_LIQUID) {
1739
1740 /* Liquid --> Solid transition */
1741
1742 const cs_real_t t_liquidus = _get_t_liquidus(alloy, conc_pre);
1743
1744 _get_dgl_mushy(alloy, t_liquidus, conc_pre, &dgldT, &dgldC);
1745
1746 const cs_real_t t_star =
1747 ( cpovL*temp + dgldT*t_liquidus + dgldC*(conc_pre - conc) ) /
1748 ( cpovL + dgldT );
1749
1750 t_bulk[c_id] = t_star;
1751
1752 gliq = 1 + (dgldT*(t_star - t_liquidus) + dgldC*(conc-conc_pre));
1753
1754 /* Make sure that the liquid fraction remains inside physical bounds */
1755
1756 gliq = fmin(fmax(0, gliq), 1.);
1757
1758 if (t_star > alloy->t_eut_sup) /* Mushy or liquid */
1759 eta_new = 1/( gliq * (1-alloy->kp) + alloy->kp );
1760 else /* Eutectic or solid */
1761 eta_new = _get_eta(alloy, conc);
1762
1763 }
1764 else {
1765 gliq = 0;
1766 eta_new = _get_eta(alloy, conc);
1767 }
1768 break;
1769
1770 case CS_SOLIDIFICATION_STATE_MUSHY:
1771 if (state_pre == CS_SOLIDIFICATION_STATE_LIQUID) {
1772 /* Liquid --> Mushy transition */
1773 const cs_real_t t_liquidus = _get_t_liquidus(alloy, conc_pre);
1774
1775 _get_dgl_mushy(alloy, t_liquidus, conc_pre, &dgldT, &dgldC);
1776
1777 const cs_real_t t_star =
1778 ( cpovL*temp + dgldT*t_liquidus + dgldC*(conc_pre - conc) ) /
1779 ( cpovL + dgldT );
1780
1781 gliq = 1 + (dgldT*(t_star - t_liquidus) + dgldC*(conc-conc_pre));
1782
1783 t_bulk[c_id] = t_star;
1784
1785 }
1786 else
1787 gliq = alloy->inv_kpm1 *
1788 ( alloy->kp - alloy->ml*conc / (temp - alloy->t_melt) );
1789
1790 /* Make sure that the liquid fraction remains inside physical bounds */
1791
1792 gliq = fmin(fmax(0, gliq), 1.);
1793
1794 eta_new = 1/( gliq * (1-alloy->kp) + alloy->kp );
1795 break;
1796
1797 case CS_SOLIDIFICATION_STATE_LIQUID:
1798 gliq = 1;
1799 eta_new = 1;
1800 break;
1801
1802 case CS_SOLIDIFICATION_STATE_EUTECTIC:
1803 if (state_pre == CS_SOLIDIFICATION_STATE_LIQUID) {
1804
1805 /* Liquid --> Eutectic transition */
1806
1807 const cs_real_t t_liquidus = _get_t_liquidus(alloy, conc_pre);
1808
1809 _get_dgl_mushy(alloy, t_liquidus, conc_pre, &dgldT, &dgldC);
1810
1811 const cs_real_t t_star =
1812 ( cpovL*temp + dgldT*t_liquidus + dgldC*(conc_pre - conc) ) /
1813 ( cpovL + dgldT );
1814
1815 t_bulk[c_id] = t_star;
1816
1817 gliq = 1 + (dgldT*(t_star - t_liquidus) + dgldC*(conc-conc_pre));
1818
1819 /* Make sure that the liquid fraction remains inside physical bounds */
1820
1821 gliq = fmin(fmax(0, gliq), 1.);
1822
1823 if (t_star > alloy->t_eut_inf)
1824 eta_new = 1/( gliq * (1-alloy->kp) + alloy->kp );
1825 else
1826 eta_new = _get_eta(alloy, conc);
1827
1828 }
1829 else {
1830
1831 const cs_real_t temp_k = alloy->tk_bulk[c_id]; /* temp_{n+1}^k */
1832
1833 /* g_l[c_id] is the value at the iterate k */
1834
1835 gliq = g_l[c_id] + cpovL * (temp_k - alloy->t_eut);
1836
1837 /* Make sure that the liquid fraction remains inside physical bounds */
1838
1839 gliq = fmin(fmax(0, gliq), 1.);
1840
1841 /* In this case Cl = C_eut = eta * Cbulk--> eta = C_eut/Cbulk */
1842
1843 eta_new = _get_eta(alloy, conc);
1844
1845 }
1846 break;
1847
1848 default:
1849 bft_error(__FILE__, __LINE__, 0, " %s: Invalid state for cell %ld\n",
1850 __func__, (long)c_id);
1851 break;
1852
1853 } /* Switch on cell state */
1854
1855 /* Update the liquid fraction and apply if needed a relaxation */
1856
1857 if (alloy->gliq_relax > 0)
1858 g_l[c_id] = (1 - alloy->gliq_relax)*gliq + alloy->gliq_relax*g_l[c_id];
1859 else
1860 g_l[c_id] = gliq;
1861
1862 /* Update eta and apply if needed a relaxation */
1863
1864 if (alloy->eta_relax > 0) {
1865 cs_real_t eta_old = alloy->eta_coef_array[c_id];
1866 alloy->eta_coef_array[c_id] =
1867 (1-alloy->eta_relax)*eta_new + alloy->eta_relax*eta_old;
1868 }
1869 else
1870 alloy->eta_coef_array[c_id] = eta_new;
1871
1872 } /* Loop on cells */
1873
1874 }
1875
1876 /*----------------------------------------------------------------------------*/
1877 /*!
1878 * \brief Update the source term for the thermal equation.
1879 * This corresponds to the binary alloy model.
1880 *
1881 * \param[in] mesh pointer to a cs_mesh_t structure
1882 * \param[in] connect pointer to a cs_cdo_connect_t structure
1883 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1884 * \param[in] ts pointer to a cs_time_step_t structure
1885 */
1886 /*----------------------------------------------------------------------------*/
1887
1888 static void
_update_thm_taylor(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)1889 _update_thm_taylor(const cs_mesh_t *mesh,
1890 const cs_cdo_connect_t *connect,
1891 const cs_cdo_quantities_t *quant,
1892 const cs_time_step_t *ts)
1893 {
1894 CS_UNUSED(mesh);
1895 cs_real_t dgldC, dgldT;
1896 cs_solidification_state_t state_k;
1897
1898 cs_solidification_t *solid = cs_solidification_structure;
1899 cs_solidification_binary_alloy_t *alloy
1900 = (cs_solidification_binary_alloy_t *)solid->model_context;
1901
1902 const cs_real_t *c_bulk = alloy->c_bulk->val;
1903 const cs_real_t *c_bulk_pre = alloy->c_bulk->val_pre;
1904 const cs_real_t *t_bulk_pre = solid->temperature->val_pre;
1905 const cs_real_t *g_l_pre = solid->g_l_field->val_pre;
1906
1907 const cs_real_t rhoL = solid->mass_density->ref_value * solid->latent_heat;
1908 const cs_real_t rhoLovdt = rhoL/ts->dt[0];
1909 const double cpovL = solid->cp->ref_value/solid->latent_heat;
1910
1911 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1912
1913 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL)
1914 continue; /* No update */
1915
1916 const cs_real_t conc = c_bulk[c_id];
1917 const cs_real_t conc_pre = c_bulk_pre[c_id];
1918 const cs_real_t temp_pre = t_bulk_pre[c_id];
1919 const cs_real_t gliq_pre = g_l_pre[c_id];
1920
1921 const cs_real_t rhocvolLovdt = quant->cell_vol[c_id] * rhoLovdt;
1922
1923 state_k = _which_state(alloy, alloy->tk_bulk[c_id], alloy->ck_bulk[c_id]);
1924
1925 /* Knowing in which part of the phase diagram we are, then we update
1926 * the value of the concentration of the liquid "solute" */
1927
1928 switch (_which_state(alloy, temp_pre, conc_pre)) {
1929
1930 case CS_SOLIDIFICATION_STATE_LIQUID:
1931 /* From the knowledge of the previous iteration, try something
1932 smarter... */
1933
1934 if (state_k == CS_SOLIDIFICATION_STATE_LIQUID) {
1935 solid->thermal_reaction_coef_array[c_id] = 0;
1936 solid->thermal_source_term_array[c_id] = 0;
1937 }
1938 else { /* Liquid --> Mushy transition */
1939 /* Liquid --> Solid transition */
1940 /* Liquid --> Eutectic transition */
1941
1942 const cs_real_t t_liquidus = _get_t_liquidus(alloy, conc_pre);
1943
1944 _get_dgl_mushy(alloy, t_liquidus, conc_pre, &dgldT, &dgldC);
1945
1946 solid->thermal_reaction_coef_array[c_id] = dgldT * rhoLovdt;
1947 solid->thermal_source_term_array[c_id] = rhocvolLovdt *
1948 ( dgldT * t_liquidus + dgldC * (conc_pre - conc) );
1949
1950 }
1951 break;
1952
1953 case CS_SOLIDIFICATION_STATE_MUSHY:
1954 {
1955 _get_dgl_mushy(alloy, temp_pre, conc_pre, &dgldT, &dgldC);
1956
1957 solid->thermal_reaction_coef_array[c_id] = dgldT * rhoLovdt;
1958 solid->thermal_source_term_array[c_id] = rhocvolLovdt *
1959 ( dgldT * temp_pre + dgldC * (conc_pre - conc) );
1960 }
1961 break;
1962
1963 case CS_SOLIDIFICATION_STATE_EUTECTIC:
1964 {
1965 const cs_real_t temp_k = alloy->tk_bulk[c_id]; /* temp_{n+1}^k */
1966
1967 solid->thermal_reaction_coef_array[c_id] = 0;
1968
1969 /* Estimate the variation of liquid fraction so that the physical
1970 bounds are satisfied for the liquid fraction) */
1971
1972 cs_real_t dgl = cpovL * (temp_k - alloy->t_eut);
1973
1974 if (dgl + gliq_pre < 0)
1975 dgl = -gliq_pre;
1976 else if (dgl + gliq_pre > 1)
1977 dgl = 1 - gliq_pre;
1978
1979 solid->thermal_source_term_array[c_id] = rhocvolLovdt * dgl;
1980
1981 }
1982 break;
1983
1984 case CS_SOLIDIFICATION_STATE_SOLID:
1985 solid->thermal_reaction_coef_array[c_id] = 0;
1986 solid->thermal_source_term_array[c_id] = 0;
1987 break;
1988
1989 default:
1990 bft_error(__FILE__, __LINE__, 0,
1991 " %s: Invalid state for cell %ld\n", __func__, (long)c_id);
1992 break;
1993
1994 } /* Switch on cell state */
1995
1996 } /* Loop on cells */
1997
1998 }
1999
2000 /*----------------------------------------------------------------------------*/
2001 /*!
2002 * \brief Update the liquid fraction in each cell and related quantities.
2003 * This corresponds to the case of a binary alloy model with no
2004 * advective source term for the solute transport.
2005 *
2006 * \param[in] mesh pointer to a cs_mesh_t structure
2007 * \param[in] connect pointer to a cs_cdo_connect_t structure
2008 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2009 * \param[in] ts pointer to a cs_time_step_t structure
2010 */
2011 /*----------------------------------------------------------------------------*/
2012
2013 static void
_update_gl_binary_path(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)2014 _update_gl_binary_path(const cs_mesh_t *mesh,
2015 const cs_cdo_connect_t *connect,
2016 const cs_cdo_quantities_t *quant,
2017 const cs_time_step_t *ts)
2018 {
2019 CS_UNUSED(mesh);
2020 CS_UNUSED(ts);
2021 cs_solidification_t *solid = cs_solidification_structure;
2022 cs_solidification_binary_alloy_t *alloy
2023 = (cs_solidification_binary_alloy_t *)solid->model_context;
2024
2025 const double L = solid->latent_heat;
2026 const cs_real_t cp0 = solid->cp->ref_value;
2027 const double cpovL = cp0/L;
2028
2029 const cs_real_t *c_bulk = alloy->c_bulk->val;
2030 const cs_real_t *c_bulk_pre = alloy->c_bulk->val_pre;
2031 cs_real_t *t_bulk = solid->temperature->val;
2032 const cs_real_t *t_bulk_pre = solid->temperature->val_pre;
2033 cs_real_t *g_l = solid->g_l_field->val;
2034 const cs_real_t *g_l_pre = solid->g_l_field->val_pre;
2035
2036 /* Update g_l values in each cell as well as the cell state and the related
2037 count */
2038
2039 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
2040
2041 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL)
2042 continue; /* No update */
2043
2044 const cs_real_t conc = c_bulk[c_id]; /* conc_{n+1}^{k+1} */
2045 const cs_real_t temp = t_bulk[c_id]; /* temp_{n+1}^{k+1} */
2046 const cs_real_t conc_pre = c_bulk_pre[c_id];
2047 const cs_real_t temp_pre = t_bulk_pre[c_id];
2048 const cs_real_t gliq_pre = g_l_pre[c_id];
2049
2050 cs_real_t dgldC, dgldT, gliq, eta_new, t_liquidus, t_solidus;
2051 cs_real_t c_star, t_star, dh, dgl;
2052 cs_solidification_state_t state, state_pre;
2053
2054 /* gliq, temp and conc iterates may not be related with the gliq(temp, conc)
2055 * function until convergence is reached. So one needs to be careful. */
2056 gliq = gliq_pre; /* default initialization to avoid a warning */
2057
2058 state = _which_state(alloy, temp, conc);
2059 state_pre = _which_state(alloy, temp_pre, conc_pre);
2060 eta_new = alloy->eta_coef_array[c_id];
2061
2062 /* Knowing in which part of the phase diagram we are and we then update
2063 * the value of the liquid fraction: g_l and the concentration of the
2064 * liquid "solute" */
2065
2066 switch (state) {
2067
2068 case CS_SOLIDIFICATION_STATE_SOLID:
2069 /* ============================= */
2070
2071 switch (state_pre) {
2072 case CS_SOLIDIFICATION_STATE_LIQUID: /* Liquid --> Solid transition */
2073
2074 t_liquidus = _get_t_liquidus(alloy, conc_pre);
2075
2076 _get_dgl_mushy(alloy, t_liquidus, conc_pre, &dgldT, &dgldC);
2077
2078 t_star = ( cpovL*temp + 1 + dgldT*t_liquidus + dgldC*(conc_pre-conc) ) /
2079 ( cpovL + dgldT );
2080
2081 gliq = 1 + (dgldT*(t_star - t_liquidus) + dgldC*(conc-conc_pre));
2082
2083 /* Make sure that the liquid fraction remains inside physical bounds */
2084
2085 gliq = fmin(fmax(0, gliq), 1.);
2086
2087 if (gliq > 0) {
2088
2089 t_solidus = _get_t_solidus(alloy, conc);
2090 if (t_star > t_solidus) /* Mushy or liquid */
2091 eta_new = 1/( gliq * (1-alloy->kp) + alloy->kp );
2092
2093 else {
2094
2095 /* Remain on the solidus line and redefine a new state */
2096
2097 t_star = t_solidus;
2098 eta_new = _get_eta(alloy, conc);
2099 }
2100 }
2101 else
2102 eta_new = _get_eta(alloy, conc);
2103
2104 t_bulk[c_id] = t_star;
2105 break;
2106
2107 case CS_SOLIDIFICATION_STATE_MUSHY: /* Mushy --> Solid transition */
2108 t_solidus = _get_t_solidus(alloy, conc);
2109 _get_dgl_mushy(alloy, temp_pre, conc_pre, &dgldT, &dgldC);
2110
2111 /* Variation of enthalpy when considering a mushy zone */
2112
2113 dh = cp0*(temp - temp_pre) +
2114 L*(dgldC*(conc-conc_pre) + dgldT*(temp-temp_pre));
2115
2116 if (conc < alloy->cs1) { /* without eutectic */
2117
2118 /* Take into account the fact that the variation of gliq is only in
2119 the mushy zone */
2120
2121 c_star = conc_pre +
2122 (dh - cp0*(temp-temp_pre) - dgldT*(t_solidus-temp_pre) )
2123 / (L*dgldC);
2124
2125 gliq = gliq_pre + dgldT*(temp-temp_pre) + dgldC*(c_star-conc_pre);
2126
2127 /* Make sure that the gliq remains inside physical bounds */
2128
2129 gliq = fmin(fmax(0, gliq), 1.);
2130 if (gliq > 0) { /* still in the mushy zone */
2131 eta_new = 1/( gliq * (1-alloy->kp) + alloy->kp );
2132 t_bulk[c_id] = t_solidus + 1e-6;
2133 }
2134 else
2135 eta_new = _get_eta(alloy, conc);
2136
2137 }
2138 else { /* with eutectic */
2139
2140 c_star = conc +
2141 (dh - cp0*(t_solidus-temp_pre)
2142 - L*(dgldC*(conc-conc_pre) + dgldT*(t_solidus-temp_pre)) )
2143 / (L*alloy->dgldC_eut);
2144
2145 if (c_star < alloy->cs1 || c_star > alloy->c_eut) {
2146 gliq = 0;
2147 eta_new = _get_eta(alloy, conc);
2148 }
2149 else {
2150
2151 gliq = gliq_pre +
2152 dgldC*(conc-conc_pre) + dgldT*(t_solidus-temp_pre) +
2153 alloy->dgldC_eut * (c_star - conc);
2154
2155 /* Make sure that the gliq remains inside physical bounds */
2156
2157 gliq = fmin(fmax(0, gliq), 1.);
2158 if (gliq > 0) /* remains on the eutectic plateau */
2159 t_bulk[c_id] = t_solidus;
2160
2161 eta_new = _get_eta(alloy, c_star);
2162
2163 } /* Invalid c_star */
2164
2165 } /* Eutectic transition taken into account */
2166 break;
2167
2168 case CS_SOLIDIFICATION_STATE_EUTECTIC: /* Eutectic --> Solid transition */
2169 _get_dgl_mushy(alloy, alloy->t_eut, conc_pre, &dgldT, &dgldC);
2170
2171 /* Variation of gl when considering how is implemented the eutectic
2172 zone */
2173
2174 dgl = dgldT*(temp-temp_pre) + alloy->dgldC_eut*(conc-conc_pre);
2175 dh = cp0*(temp -temp_pre) + dgl*L;
2176
2177 /* If one remains on the eutectic plateau, then the concentration should
2178 be c_star w.r.t. dh = dgldC_eut * (C* - Cn) since Tk+1 = Tn = Teut */
2179
2180 c_star = conc_pre + dh/(L*alloy->dgldC_eut);
2181
2182 if (c_star < alloy->cs1 || c_star > alloy->c_eut) {
2183
2184 /* In fact the final state is below the eutectic plateau */
2185
2186 gliq = 0;
2187 eta_new = _get_eta(alloy, conc);
2188
2189 }
2190 else {
2191
2192 gliq = gliq_pre + alloy->dgldC_eut*(c_star-conc_pre);
2193 gliq = fmin(fmax(0, gliq), 1.);
2194 eta_new = _get_eta(alloy, c_star);
2195 if (gliq > 0)
2196 t_bulk[c_id] = alloy->t_eut;
2197
2198 }
2199 break;
2200
2201 default: /* Solid --> solid */
2202 gliq = 0;
2203 if (gliq_pre > 0) /* Otherwise keep the same value for eta */
2204 eta_new = _get_eta(alloy, conc);
2205 break;
2206
2207 } /* Switch on the previous state */
2208 break;
2209
2210
2211 case CS_SOLIDIFICATION_STATE_MUSHY:
2212 /* ============================= */
2213
2214 switch (state_pre) {
2215
2216 case CS_SOLIDIFICATION_STATE_LIQUID: /* Liquid --> Mushy transition */
2217 t_liquidus = _get_t_liquidus(alloy, conc_pre);
2218
2219 _get_dgl_mushy(alloy, t_liquidus, conc_pre, &dgldT, &dgldC);
2220
2221 t_star = ( cpovL*temp + dgldT*t_liquidus + dgldC*(conc_pre - conc) ) /
2222 ( cpovL + dgldT );
2223
2224 gliq = 1 + (dgldT*(t_star-t_liquidus) + dgldC*(conc-conc_pre));
2225
2226 t_bulk[c_id] = t_star;
2227 break;
2228
2229 case CS_SOLIDIFICATION_STATE_MUSHY:
2230 _get_dgl_mushy(alloy, temp_pre, conc_pre, &dgldT, &dgldC);
2231
2232 gliq = gliq_pre +
2233 (dgldT*(temp-temp_pre) + dgldC*(conc-conc_pre));
2234 break;
2235
2236 default:
2237 gliq = alloy->inv_kpm1 *
2238 ( alloy->kp - alloy->ml*conc / (temp - alloy->t_melt) );
2239
2240 } /* End of switch on the previous state */
2241
2242 /* Make sure that the liquid fraction remains inside physical bounds */
2243 gliq = fmin(fmax(0, gliq), 1.);
2244
2245 eta_new = 1/( gliq * (1-alloy->kp) + alloy->kp );
2246 break;
2247
2248
2249 case CS_SOLIDIFICATION_STATE_LIQUID:
2250 /* ============================== */
2251 gliq = 1;
2252 eta_new = 1;
2253 break;
2254
2255
2256 case CS_SOLIDIFICATION_STATE_EUTECTIC:
2257 /* ================================ */
2258
2259 switch (state_pre) {
2260
2261 case CS_SOLIDIFICATION_STATE_LIQUID: /* Liquid --> Eutectic transition */
2262 t_liquidus = _get_t_liquidus(alloy, conc_pre);
2263
2264 _get_dgl_mushy(alloy, t_liquidus, conc_pre, &dgldT, &dgldC);
2265
2266 t_star = ( cpovL*temp + dgldT*t_liquidus + dgldC*(conc_pre - conc) ) /
2267 ( cpovL + dgldT );
2268
2269 t_bulk[c_id] = t_star;
2270
2271 gliq = 1 + (dgldT*(t_star - t_liquidus) + dgldC*(conc-conc_pre));
2272
2273 /* Make sure that the liquid fraction remains inside physical bounds */
2274 gliq = fmin(fmax(0, gliq), 1.);
2275
2276 if (t_star > alloy->t_eut_inf)
2277 eta_new = 1/( gliq * (1-alloy->kp) + alloy->kp );
2278 else
2279 eta_new = _get_eta(alloy, conc);
2280 break;
2281
2282 case CS_SOLIDIFICATION_STATE_MUSHY: /* Mushy --> Eutectic transition */
2283 assert(conc > alloy->cs1);
2284
2285 /* First part of the path in the mushy zone */
2286 _get_dgl_mushy(alloy, temp_pre, conc_pre, &dgldT, &dgldC);
2287
2288 gliq = g_l_pre[c_id] +
2289 alloy->dgldC_eut*(conc - conc_pre) + dgldT*(alloy->t_eut - temp_pre);
2290
2291 /* Make sure that the liquid fraction remains inside physical bounds */
2292 gliq = fmin(fmax(0, gliq), 1.);
2293
2294 eta_new = _get_eta(alloy, conc);
2295 break;
2296
2297 default: /* eutectic --> eutectic or solid --> eutectic */
2298 _get_dgl_mushy(alloy, alloy->t_eut, conc_pre, &dgldT, &dgldC);
2299
2300 /* Variation of gl when considering how is implemented the eutectic
2301 zone */
2302 dgl = dgldT*(temp-temp_pre) + alloy->dgldC_eut*(conc-conc_pre);
2303 dh = cp0*(temp -temp_pre) + dgl*L;
2304
2305 /* If one remains on the eutectic plateau, then the concentration should
2306 be c_star w.r.t. dh = dgldC_eut * (C* - Cn) since Tk+1 = Tn = Teut */
2307 c_star = conc_pre + dh/(L*alloy->dgldC_eut);
2308
2309 if (c_star < alloy->cs1 || c_star > alloy->c_eut) {
2310
2311 gliq = (conc - alloy->cs1)*alloy->dgldC_eut;
2312
2313 /* In this case Cl = C_eut = eta * Cbulk--> eta = C_eut/Cbulk */
2314 eta_new = _get_eta(alloy, conc);
2315
2316 }
2317 else {
2318
2319 gliq = gliq_pre + alloy->dgldC_eut*(c_star-conc_pre);
2320 if (gliq > 0) /* Remains on the eutectic plateau */
2321 t_bulk[c_id] = alloy->t_eut;
2322
2323 eta_new = _get_eta(alloy, c_star);
2324
2325 }
2326
2327 /* Make sure that the liquid fraction remains inside physical bounds */
2328 gliq = fmin(fmax(0, gliq), 1.);
2329
2330 break;
2331
2332 }
2333 break;
2334
2335 default:
2336 bft_error(__FILE__, __LINE__, 0, " %s: Invalid state for cell %ld\n",
2337 __func__, (long)c_id);
2338 break;
2339
2340 } /* Switch on cell state */
2341
2342 /* Update the liquid fraction and apply if needed a relaxation */
2343 if (alloy->gliq_relax > 0)
2344 g_l[c_id] = (1 - alloy->gliq_relax)*gliq + alloy->gliq_relax*g_l[c_id];
2345 else
2346 g_l[c_id] = gliq;
2347
2348 /* Update eta and apply if needed a relaxation */
2349 if (alloy->eta_relax > 0) {
2350 cs_real_t eta_old = alloy->eta_coef_array[c_id];
2351 alloy->eta_coef_array[c_id] =
2352 (1-alloy->eta_relax)*eta_new + alloy->eta_relax*eta_old;
2353 }
2354 else
2355 alloy->eta_coef_array[c_id] = eta_new;
2356
2357 } /* Loop on cells */
2358
2359 }
2360
2361 /*----------------------------------------------------------------------------*/
2362 /*!
2363 * \brief Update the source term for the thermal equation.
2364 * This corresponds to the binary alloy model.
2365 *
2366 * \param[in] mesh pointer to a cs_mesh_t structure
2367 * \param[in] connect pointer to a cs_cdo_connect_t structure
2368 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2369 * \param[in] ts pointer to a cs_time_step_t structure
2370 */
2371 /*----------------------------------------------------------------------------*/
2372
2373 static void
_update_thm_binary_path(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)2374 _update_thm_binary_path(const cs_mesh_t *mesh,
2375 const cs_cdo_connect_t *connect,
2376 const cs_cdo_quantities_t *quant,
2377 const cs_time_step_t *ts)
2378 {
2379 CS_UNUSED(mesh);
2380 cs_solidification_t *solid = cs_solidification_structure;
2381 cs_solidification_binary_alloy_t *alloy
2382 = (cs_solidification_binary_alloy_t *)solid->model_context;
2383
2384 const cs_real_t *c_bulk = alloy->c_bulk->val;
2385 const cs_real_t *c_bulk_pre = alloy->c_bulk->val_pre;
2386 const cs_real_t *t_bulk = solid->temperature->val;
2387 const cs_real_t *t_bulk_pre = solid->temperature->val_pre;
2388
2389 const cs_real_t rhoL = solid->mass_density->ref_value * solid->latent_heat;
2390 const cs_real_t rhoLovdt = rhoL/ts->dt[0];
2391
2392 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
2393
2394 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL)
2395 continue; /* No update */
2396
2397 const cs_real_t conc_kp1 = c_bulk[c_id]; /* Solute transport solved */
2398 const cs_real_t conc_k = alloy->ck_bulk[c_id];
2399 const cs_real_t temp_k = t_bulk[c_id];
2400
2401 const cs_real_t conc_pre = c_bulk_pre[c_id];
2402 const cs_real_t temp_pre = t_bulk_pre[c_id];
2403
2404 const cs_real_t rhocvolLovdt = quant->cell_vol[c_id] * rhoLovdt;
2405 cs_real_t dgldC, dgldT, t_solidus, t_liquidus;
2406
2407 cs_solidification_state_t state_k = _which_state(alloy, temp_k, conc_k);
2408
2409 /* Knowing in which part of the phase diagram we are, then we update
2410 * the value of the concentration of the liquid "solute" */
2411 switch (_which_state(alloy, temp_pre, conc_pre)) {
2412
2413 case CS_SOLIDIFICATION_STATE_LIQUID:
2414 /* ==============================
2415 * From the knowledge of the previous iteration, try something smarter...
2416 */
2417 switch (state_k) {
2418 case CS_SOLIDIFICATION_STATE_MUSHY: /* Liquid --> Mushy */
2419 t_liquidus = _get_t_liquidus(alloy, conc_pre);
2420
2421 _get_dgl_mushy(alloy, t_liquidus, conc_pre, &dgldT, &dgldC);
2422
2423 solid->thermal_reaction_coef_array[c_id] = dgldT * rhoLovdt;
2424 solid->thermal_source_term_array[c_id] = rhocvolLovdt *
2425 ( dgldT * t_liquidus + dgldC * (conc_pre-conc_kp1) );
2426 break;
2427
2428 case CS_SOLIDIFICATION_STATE_EUTECTIC: /* Liquid --> Eutectic */
2429 case CS_SOLIDIFICATION_STATE_SOLID: /* Liquid --> Solid */
2430 t_liquidus = _get_t_liquidus(alloy, conc_pre);
2431 t_solidus = _get_t_solidus(alloy, conc_kp1);
2432
2433 _get_dgl_mushy(alloy, t_liquidus, conc_pre, &dgldT, &dgldC);
2434
2435 solid->thermal_reaction_coef_array[c_id] = 0;
2436 solid->thermal_source_term_array[c_id] = rhocvolLovdt *
2437 ( dgldT * (t_liquidus-t_solidus) + dgldC * (conc_pre-conc_kp1) );
2438 break;
2439
2440 default: /* Liquid */
2441 solid->thermal_reaction_coef_array[c_id] = 0;
2442 solid->thermal_source_term_array[c_id] = 0;
2443 break;
2444
2445 } /* End of switch on the state k */
2446 break;
2447
2448 case CS_SOLIDIFICATION_STATE_MUSHY:
2449 /* ============================= */
2450 {
2451 _get_dgl_mushy(alloy, temp_pre, conc_pre, &dgldT, &dgldC);
2452
2453 switch (state_k) {
2454
2455 case CS_SOLIDIFICATION_STATE_SOLID: /* Mushy --> Solid transition */
2456 solid->thermal_reaction_coef_array[c_id] = dgldT * rhoLovdt;
2457 if (conc_kp1 < alloy->cs1) /* Part without eutectic */
2458 solid->thermal_source_term_array[c_id] = rhocvolLovdt *
2459 ( dgldT*temp_pre + dgldC*(conc_pre-conc_kp1) );
2460 else /* Part with eutectic */
2461 solid->thermal_source_term_array[c_id] = rhocvolLovdt *
2462 ( dgldT*temp_pre + alloy->dgldC_eut*(conc_pre-conc_kp1) );
2463 break;
2464
2465 case CS_SOLIDIFICATION_STATE_EUTECTIC: /* Mushy --> Eutectic */
2466 assert(conc_kp1 > alloy->cs1);
2467
2468 /* First part in the mushy zone but not the final part */
2469 solid->thermal_reaction_coef_array[c_id] = dgldT * rhoLovdt;
2470 solid->thermal_source_term_array[c_id] = rhocvolLovdt *
2471 ( dgldT*temp_pre + alloy->dgldC_eut*(conc_pre-conc_kp1) );
2472 break;
2473
2474 default:
2475 solid->thermal_reaction_coef_array[c_id] = dgldT * rhoLovdt;
2476 solid->thermal_source_term_array[c_id] = rhocvolLovdt *
2477 ( dgldT * temp_pre + dgldC * (conc_pre - conc_kp1) );
2478 break;
2479
2480 }
2481
2482 }
2483 break;
2484
2485 case CS_SOLIDIFICATION_STATE_EUTECTIC:
2486 /* ================================ */
2487 {
2488 cs_real_t r_coef = 0;
2489 cs_real_t s_coef = alloy->dgldC_eut * (conc_pre - conc_kp1);
2490
2491 if (solid->options & CS_SOLIDIFICATION_WITH_PENALIZED_EUTECTIC) {
2492 if (state_k == CS_SOLIDIFICATION_STATE_EUTECTIC ||
2493 state_k == CS_SOLIDIFICATION_STATE_SOLID) {
2494 if (conc_kp1 > alloy->cs1 && conc_kp1 < alloy->c_eut) {
2495 _get_dgl_mushy(alloy, temp_pre, conc_pre, &dgldT, &dgldC);
2496 r_coef = dgldT * rhoLovdt;
2497 s_coef += dgldT*alloy->t_eut;
2498 }
2499 }
2500 }
2501
2502 solid->thermal_reaction_coef_array[c_id] = r_coef;
2503 solid->thermal_source_term_array[c_id] = rhocvolLovdt * s_coef;
2504
2505 }
2506 break;
2507
2508 case CS_SOLIDIFICATION_STATE_SOLID:
2509 /* ============================= */
2510 solid->thermal_reaction_coef_array[c_id] = 0;
2511 solid->thermal_source_term_array[c_id] = 0;
2512 break;
2513
2514 default:
2515 bft_error(__FILE__, __LINE__, 0,
2516 " %s: Invalid state for cell %ld\n", __func__, (long)c_id);
2517 break;
2518
2519 } /* Switch on cell state */
2520
2521 } /* Loop on cells */
2522
2523 }
2524
2525 /*----------------------------------------------------------------------------*/
2526 /*!
2527 * \brief Update the source term for the thermal equation.
2528 * This function is related to the Stefan problem with a liquid fraction
2529 * being a step function w.r.t. the temperature
2530 *
2531 * \param[in] mesh pointer to a cs_mesh_t structure
2532 * \param[in] connect pointer to a cs_cdo_connect_t structure
2533 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2534 * \param[in] ts pointer to a cs_time_step_t structure
2535 */
2536 /*----------------------------------------------------------------------------*/
2537
2538 static void
_update_thm_stefan(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)2539 _update_thm_stefan(const cs_mesh_t *mesh,
2540 const cs_cdo_connect_t *connect,
2541 const cs_cdo_quantities_t *quant,
2542 const cs_time_step_t *ts)
2543 {
2544 if (mesh->n_cells < 1)
2545 return;
2546
2547 cs_real_t rho_c, rhoLovdt;
2548
2549 cs_solidification_t *solid = cs_solidification_structure;
2550
2551 const cs_real_t Lovdt = solid->latent_heat/ts->dt[0];
2552 const cs_real_t *g_l = solid->g_l_field->val;
2553 const cs_real_t *g_l_pre = solid->g_l_field->val_pre;
2554 const cs_real_t *vol = quant->cell_vol;
2555
2556 bool rho_is_uniform = cs_property_is_uniform(solid->mass_density);
2557
2558 /* Use the first cell to set the value */
2559
2560 if (rho_is_uniform) {
2561 rho_c = cs_property_get_cell_value(0, ts->t_cur, solid->mass_density);
2562 rhoLovdt = rho_c * Lovdt;
2563 }
2564
2565 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
2566 for (cs_lnum_t c = 0; c < quant->n_cells; c++) {
2567
2568 /* Retrieve the value of the properties */
2569
2570 if (!rho_is_uniform) {
2571 rho_c = cs_property_get_cell_value(c, ts->t_cur, solid->mass_density);
2572 rhoLovdt = rho_c * Lovdt;
2573 }
2574
2575 if (connect->cell_flag[c] & CS_FLAG_SOLID_CELL) /* Tag as solid for all the
2576 computation */
2577 continue; /* No update: 0 by default */
2578
2579 /* reaction_coef_array is set to zero. Only the source term is updated */
2580
2581 if (fabs(g_l[c] - g_l_pre[c]) > 0)
2582 solid->thermal_source_term_array[c] = rhoLovdt*vol[c]*(g_l_pre[c]-g_l[c]);
2583 else
2584 solid->thermal_source_term_array[c] = 0;
2585
2586 } /* Loop on cells */
2587
2588 }
2589
2590 /*----------------------------------------------------------------------------*/
2591 /*!
2592 * \brief Update the liquid fraction in each cell and the temperature if
2593 * needed. Case of the Stefan model.
2594 *
2595 * \param[in] mesh pointer to a cs_mesh_t structure
2596 * \param[in] connect pointer to a cs_cdo_connect_t structure
2597 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2598 * \param[in] ts pointer to a cs_time_step_t structure
2599 */
2600 /*----------------------------------------------------------------------------*/
2601
2602 static void
_update_gl_stefan(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)2603 _update_gl_stefan(const cs_mesh_t *mesh,
2604 const cs_cdo_connect_t *connect,
2605 const cs_cdo_quantities_t *quant,
2606 const cs_time_step_t *ts)
2607 {
2608 CS_UNUSED(ts);
2609
2610 if (mesh->n_cells < 1)
2611 return;
2612
2613 cs_real_t cp_c, cpovL;
2614
2615 cs_solidification_t *solid = cs_solidification_structure;
2616 cs_solidification_stefan_t *model =
2617 (cs_solidification_stefan_t *)solid->model_context;
2618
2619 bool cp_is_uniform = cs_property_is_uniform(solid->cp);
2620
2621 cs_real_t *temp = solid->temperature->val;
2622 cs_real_t *g_l = solid->g_l_field->val;
2623
2624 /* Use the first cell to set the value */
2625
2626 if (cp_is_uniform) {
2627 cp_c = cs_property_get_cell_value(0, ts->t_cur, solid->cp);
2628 cpovL = cp_c/solid->latent_heat;
2629 }
2630
2631 /* Update g_l values in each cell as well as the cell state and the related
2632 count */
2633
2634 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
2635 for (cs_lnum_t c = 0; c < quant->n_cells; c++) {
2636
2637 /* Retrieve the value of the property */
2638
2639 if (!cp_is_uniform) {
2640 cp_c = cs_property_get_cell_value(c, ts->t_cur, solid->cp);
2641 cpovL = cp_c/solid->latent_heat;
2642 }
2643
2644 if (connect->cell_flag[c] & CS_FLAG_SOLID_CELL)
2645 continue; /* Tag as solid during all the computation
2646 => No update: 0 by default */
2647
2648 if (temp[c] > model->t_change) {
2649
2650 if (g_l[c] < 1) { /* Not in a stable state */
2651
2652 /* Compute a new g_l */
2653
2654 g_l[c] += cpovL * (temp[c] - model->t_change);
2655 if (g_l[c] < 1) {
2656
2657 temp[c] = model->t_change;
2658 solid->cell_state[c] = CS_SOLIDIFICATION_STATE_MUSHY;
2659
2660 }
2661 else { /* Overshoot of the liquid fraction */
2662
2663 solid->cell_state[c] = CS_SOLIDIFICATION_STATE_LIQUID;
2664 temp[c] = model->t_change + 1./cpovL * (g_l[c] - 1);
2665 g_l[c] = 1.;
2666
2667 }
2668
2669 }
2670 else { /* g_l = 1, stable state */
2671
2672 solid->cell_state[c] = CS_SOLIDIFICATION_STATE_LIQUID;
2673
2674 }
2675
2676 }
2677 else { /* T < T_ch */
2678
2679 if (g_l[c] > 0) { /* Not in a stable state */
2680
2681 /* Compute a new g_l */
2682
2683 g_l[c] += cpovL * (temp[c] - model->t_change);
2684
2685 if (g_l[c] < 0) { /* Undershoot of the liquid fraction */
2686
2687 solid->cell_state[c] = CS_SOLIDIFICATION_STATE_SOLID;
2688 temp[c] = model->t_change + 1./cpovL * g_l[c];
2689 g_l[c] = 0.;
2690
2691 }
2692 else {
2693
2694 temp[c] = model->t_change;
2695 solid->cell_state[c] = CS_SOLIDIFICATION_STATE_MUSHY;
2696
2697 }
2698
2699 }
2700 else { /* g_l = 0, stable state */
2701
2702 solid->cell_state[c] = CS_SOLIDIFICATION_STATE_SOLID;
2703
2704 }
2705
2706 }
2707
2708 } /* Loop on cells */
2709
2710 }
2711
2712 /*----------------------------------------------------------------------------*/
2713 /*!
2714 * \brief Function aims at computing the new couple
2715 * (temperature,liquid fraction) defining the new state
2716 *
2717 * \param[in] mesh pointer to a cs_mesh_t structure
2718 * \param[in] connect pointer to a cs_cdo_connect_t structure
2719 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2720 * \param[in] time_step pointer to a cs_time_step_t structure
2721 */
2722 /*----------------------------------------------------------------------------*/
2723
2724 static void
_stefan_thermal_non_linearities(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)2725 _stefan_thermal_non_linearities(const cs_mesh_t *mesh,
2726 const cs_cdo_connect_t *connect,
2727 const cs_cdo_quantities_t *quant,
2728 const cs_time_step_t *time_step)
2729 {
2730 cs_solidification_t *solid = cs_solidification_structure;
2731 cs_solidification_stefan_t *s_model
2732 = (cs_solidification_stefan_t *)solid->model_context;
2733
2734 const size_t csize = quant->n_cells*sizeof(cs_real_t);
2735 const cs_equation_t *t_eq = solid->thermal_sys->thermal_eq;
2736
2737 /* Retrieve the current values */
2738 cs_real_t *temp = cs_equation_get_cell_values(t_eq, false);
2739 cs_real_t *g_l = solid->g_l_field->val;
2740 cs_real_t *enthalpy = solid->enthalpy->val;
2741
2742 cs_real_t *hk = NULL; /* enthalpy h^{n+1,k} */
2743 BFT_MALLOC(hk, quant->n_cells, cs_real_t);
2744 memcpy(hk, enthalpy, csize);
2745
2746 /* Non-linear iterations (k) are performed to converge on the relation
2747 * h^{n+1,k+1} = h^{n+1,k} + eps with eps a user-defined tolerance
2748 *
2749 * h = rho.cp (Temp - Tref) + rho.L.gliq
2750 *
2751 * One sets:
2752 * T^{n+1,0} = T^n and gl^{n+1,0} = gl^n
2753 */
2754
2755 cs_equation_current_to_previous(t_eq);
2756 cs_field_current_to_previous(solid->g_l_field);
2757 cs_field_current_to_previous(solid->enthalpy);
2758
2759 /* Initialize loop stopping criteria */
2760
2761 cs_real_t delta_h = 1 + s_model->max_delta_h;
2762 int iter = 0;
2763
2764 while ( delta_h > s_model->max_delta_h && iter < s_model->n_iter_max) {
2765
2766 /* Compute the new thermal source term */
2767
2768 s_model->update_thm_st(mesh, connect, quant, time_step);
2769
2770 /* Solve the thermal system */
2771
2772 cs_thermal_system_compute(false, /* No cur2prev inside a non-linear
2773 iterative process */
2774 mesh, connect, quant, time_step);
2775
2776 /* Compute the new liquid fraction (and update the temperature if needed) */
2777
2778 s_model->update_gl(mesh, connect, quant, time_step);
2779
2780 /* Now compute the enthalpy knowing the temperature and the liquid
2781 * fraction.
2782 * enthalpy stores k+1,n+1 and hk stores k,n+1
2783 */
2784
2785 _compute_enthalpy(quant,
2786 time_step->t_cur, /* t_eval */
2787 temp, /* temperature */
2788 g_l, /* liquid fraction */
2789 s_model->t_change, /* temp_ref */
2790 solid->latent_heat, /* latent heat coeff. */
2791 solid->mass_density, /* rho */
2792 solid->cp, /* cp */
2793 enthalpy); /* computed enthalpy */
2794
2795 delta_h = -1;
2796 for (cs_lnum_t c = 0; c < quant->n_cells; c++) {
2797
2798 cs_real_t dh = fabs(enthalpy[c] - hk[c]);
2799 hk[c] = enthalpy[c];
2800
2801 if (dh > delta_h)
2802 delta_h = dh;
2803
2804 } /* Loop on cells */
2805
2806 iter++;
2807 if (solid->verbosity > 1)
2808 cs_log_printf(CS_LOG_DEFAULT,
2809 "### Solidification.NL: k= %d | delta_enthalpy= %5.3e\n",
2810 iter, delta_h);
2811
2812 } /* Until convergence */
2813
2814 BFT_FREE(hk);
2815
2816 /* Monitoring */
2817
2818 if (solid->verbosity > 1)
2819 cs_log_printf(CS_LOG_DEFAULT,
2820 "## Solidification: Stop after %d iters, delta = %5.3e\n",
2821 iter, delta_h);
2822
2823 _monitor_cell_state(connect, quant, solid);
2824
2825 /* Parallel synchronization of the number of cells in each state */
2826
2827 cs_parall_sum(CS_SOLIDIFICATION_N_STATES, CS_GNUM_TYPE, solid->n_g_cells);
2828 }
2829
2830 /*----------------------------------------------------------------------------*/
2831 /*!
2832 * \brief Compute the new system state (temperature, liquid fraction) using
2833 * the methodology defined in Voller & Prakash (87)
2834 *
2835 * \param[in] mesh pointer to a cs_mesh_t structure
2836 * \param[in] connect pointer to a cs_cdo_connect_t structure
2837 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2838 * \param[in] time_step pointer to a cs_time_step_t structure
2839 */
2840 /*----------------------------------------------------------------------------*/
2841
2842 static void
_voller_prakash_87(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)2843 _voller_prakash_87(const cs_mesh_t *mesh,
2844 const cs_cdo_connect_t *connect,
2845 const cs_cdo_quantities_t *quant,
2846 const cs_time_step_t *time_step)
2847 {
2848 cs_solidification_t *solid = cs_solidification_structure;
2849
2850 /* Solidification process with a pure component without segregation */
2851
2852 cs_solidification_voller_t *v_model =
2853 (cs_solidification_voller_t *)solid->model_context;
2854
2855 /* Solve the thermal system */
2856
2857 cs_thermal_system_compute(true, /* operate a cur2prev operation inside */
2858 mesh, connect, quant, time_step);
2859
2860 /* Update fields and properties which are related to solved variables */
2861
2862 cs_field_current_to_previous(solid->g_l_field);
2863
2864 v_model->update_gl(mesh, connect, quant, time_step);
2865
2866 v_model->update_thm_st(mesh, connect, quant, time_step);
2867
2868 /* Post-processing */
2869
2870 if (solid->post_flag & CS_SOLIDIFICATION_POST_ENTHALPY)
2871 _compute_enthalpy(quant,
2872 time_step->t_cur, /* t_eval */
2873 solid->temperature->val, /* temperature */
2874 solid->g_l_field->val, /* liquid fraction */
2875 v_model->t_solidus, /* temp_ref */
2876 solid->latent_heat, /* latent heat coeff. */
2877 solid->mass_density, /* rho */
2878 solid->cp, /* cp */
2879 solid->enthalpy->val); /* computed enthalpy */
2880
2881 /* Monitoring */
2882
2883 _monitor_cell_state(connect, quant, solid);
2884
2885 /* At this stage, the number of solid cells is a local count
2886 * Set the enforcement of the velocity for solid cells */
2887
2888 _enforce_solid_cells(connect, quant);
2889
2890 /* Parallel synchronization of the number of cells in each state (It should
2891 be done after _enforce_solid_cells() */
2892
2893 cs_parall_sum(CS_SOLIDIFICATION_N_STATES, CS_GNUM_TYPE, solid->n_g_cells);
2894 }
2895
2896 /*----------------------------------------------------------------------------*/
2897 /*!
2898 * \brief Compute the new system state (temperature, liquid fraction) using
2899 * the methodology defined in Voller & Prakash (87) but also taking
2900 * into account the non-linearities stemming from the thermal source
2901 * term
2902 *
2903 * \param[in] mesh pointer to a cs_mesh_t structure
2904 * \param[in] connect pointer to a cs_cdo_connect_t structure
2905 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2906 * \param[in] time_step pointer to a cs_time_step_t structure
2907 */
2908 /*----------------------------------------------------------------------------*/
2909
2910 static void
_voller_non_linearities(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)2911 _voller_non_linearities(const cs_mesh_t *mesh,
2912 const cs_cdo_connect_t *connect,
2913 const cs_cdo_quantities_t *quant,
2914 const cs_time_step_t *time_step)
2915 {
2916 cs_solidification_t *solid = cs_solidification_structure;
2917
2918 /* Solidification process with a pure component without segregation */
2919
2920 cs_solidification_voller_t *v_model =
2921 (cs_solidification_voller_t *)solid->model_context;
2922
2923 cs_iter_algo_t *algo = v_model->nl_algo;
2924
2925 assert(algo != NULL);
2926
2927 const size_t csize = quant->n_cells*sizeof(cs_real_t);
2928
2929 /* Retrieve the current values */
2930
2931 cs_real_t *hkp1 = solid->enthalpy->val;
2932 cs_real_t *hk = NULL; /* enthalpy ^{n+1,k} */
2933 BFT_MALLOC(hk, quant->n_cells, cs_real_t);
2934
2935 /* Initialize the stopping criteria */
2936
2937 cs_iter_algo_reset(algo);
2938
2939 algo->normalization = sqrt(cs_cdo_blas_square_norm_pcsp(hkp1));
2940 if (algo->normalization < cs_math_zero_threshold)
2941 algo->normalization = 1.0;
2942
2943 /* Non-linear iterations (k) are performed to converge on the relation
2944 * h^{n+1,k+1} = h^{n+1,k} + eps with eps a user-defined tolerance
2945 *
2946 * h = rho.cp (Temp - Tref) + rho.L.gliq
2947 *
2948 * One sets:
2949 * T^{n+1,0} = T^n and gl^{n+1,0} = gl^n
2950 */
2951
2952 cs_equation_current_to_previous(solid->thermal_sys->thermal_eq);
2953 cs_field_current_to_previous(solid->g_l_field);
2954 cs_field_current_to_previous(solid->enthalpy);
2955
2956 do {
2957
2958 /* Compute the new thermal source term */
2959
2960 v_model->update_thm_st(mesh, connect, quant, time_step);
2961
2962 /* Solve the thermal system */
2963
2964 cs_thermal_system_compute(false, /* No cur2prev inside a non-linear
2965 iterative process */
2966 mesh, connect, quant, time_step);
2967
2968 /* Compute the new liquid fraction (and update the temperature if needed) */
2969
2970 v_model->update_gl(mesh, connect, quant, time_step);
2971
2972 /* Now compute the enthalpy knowing the temperature and the liquid
2973 * fraction.
2974 * enthalpy stores k+1,n+1 and hk stores k,n+1
2975 */
2976
2977 memcpy(hk, hkp1, csize);
2978
2979 _compute_enthalpy(quant,
2980 time_step->t_cur, /* t_eval */
2981 solid->temperature->val, /* temperature */
2982 solid->g_l_field->val, /* liquid fraction */
2983 v_model->t_solidus, /* temp_ref */
2984 solid->latent_heat, /* latent heat coeff. */
2985 solid->mass_density, /* rho */
2986 solid->cp, /* cp */
2987 hkp1); /* computed enthalpy */
2988
2989 } /* Until convergence */
2990 while (_check_nl_cvg(v_model->nl_algo_type,
2991 hk, hkp1, algo) == CS_SLES_ITERATING);
2992
2993 BFT_FREE(hk);
2994
2995 /* Monitoring */
2996
2997 if (solid->verbosity > 0)
2998 cs_log_printf(CS_LOG_DEFAULT,
2999 "## Solidification: Stop non-linear algo. after %d iters,"
3000 " residual = %5.3e\n",
3001 algo->n_algo_iter, algo->res);
3002
3003 /* Monitoring */
3004
3005 _monitor_cell_state(connect, quant, solid);
3006
3007 /* At this stage, the number of solid cells is a local count
3008 * Set the enforcement of the velocity for solid cells */
3009
3010 _enforce_solid_cells(connect, quant);
3011
3012 /* Parallel synchronization of the number of cells in each state (It should
3013 be done after _enforce_solid_cells() */
3014
3015 cs_parall_sum(CS_SOLIDIFICATION_N_STATES, CS_GNUM_TYPE, solid->n_g_cells);
3016 }
3017
3018 /*----------------------------------------------------------------------------*/
3019 /*!
3020 * \brief Function aims at computing the new temperature/bulk concentration
3021 * state for the next iteration as well as updating all related
3022 * quantities
3023 *
3024 * \param[in] mesh pointer to a cs_mesh_t structure
3025 * \param[in] connect pointer to a cs_cdo_connect_t structure
3026 * \param[in] quant pointer to a cs_cdo_quantities_t structure
3027 * \param[in] time_step pointer to a cs_time_step_t structure
3028 */
3029 /*----------------------------------------------------------------------------*/
3030
3031 static void
_default_binary_coupling(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)3032 _default_binary_coupling(const cs_mesh_t *mesh,
3033 const cs_cdo_connect_t *connect,
3034 const cs_cdo_quantities_t *quant,
3035 const cs_time_step_t *time_step)
3036 {
3037 cs_solidification_t *solid = cs_solidification_structure;
3038 cs_solidification_binary_alloy_t *alloy
3039 = (cs_solidification_binary_alloy_t *)solid->model_context;
3040
3041 const size_t csize = quant->n_cells*sizeof(cs_real_t);
3042 const cs_equation_t *c_eq = alloy->solute_equation;
3043 const cs_equation_t *t_eq = solid->thermal_sys->thermal_eq;
3044 const cs_real_t rho0 = solid->mass_density->ref_value;
3045
3046 cs_real_t *temp = cs_equation_get_cell_values(t_eq, false);
3047 cs_real_t *conc = cs_equation_get_cell_values(c_eq, false);
3048 cs_real_t *g_l = solid->g_l_field->val;
3049
3050 /* Compute the state at t^(n+1) knowing that at state t^(n) */
3051 if (solid->options & CS_SOLIDIFICATION_USE_EXTRAPOLATION) {
3052
3053 /* At this stage (i.e. before previous to current: val = n, val_pre = n-1 */
3054 cs_real_t *temp_pre = cs_equation_get_cell_values(t_eq, true);
3055 cs_real_t *conc_pre = cs_equation_get_cell_values(c_eq, true);
3056
3057 /* extrapolation at f_{n+1} = 2*f_n - f_{n-1} */
3058 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
3059 alloy->tx_bulk[c_id] = 2*temp[c_id] - temp_pre[c_id];
3060 alloy->cx_bulk[c_id] = 2*conc[c_id] - conc_pre[c_id];
3061 }
3062
3063 } /* Extrapolation is requested */
3064
3065 /* Non-linear iterations (k) are also performed to converge on the relation
3066 * gliq^{k+1} = gliq(temp^{k+1}, conc^{k+1})
3067 *
3068 * Cbulk^{0}_{n+1} = Cbulk_{n}
3069 * Tbulk^{0}_{n+1} = Tbulk_{n}
3070 * gl^{0}_{n+1} = gl_{n}
3071 */
3072
3073 cs_equation_current_to_previous(c_eq);
3074 cs_equation_current_to_previous(t_eq);
3075 cs_field_current_to_previous(solid->g_l_field);
3076
3077 /* At the beginning, field_{n+1}^{k=0} = field_n */
3078
3079 memcpy(alloy->tk_bulk, temp, csize);
3080 memcpy(alloy->ck_bulk, conc, csize);
3081
3082 cs_real_t delta_temp = 1 + alloy->delta_tolerance;
3083 cs_real_t delta_cbulk = 1 + alloy->delta_tolerance;
3084
3085 alloy->iter = 0;
3086 while ( ( delta_temp > alloy->delta_tolerance ||
3087 delta_cbulk > alloy->delta_tolerance ) &&
3088 alloy->iter < alloy->n_iter_max) {
3089
3090 /* Solve Cbulk^(k+1)_{n+1} knowing Cbulk^{k}_{n+1} */
3091
3092 cs_equation_solve(false, /* No cur2prev inside a non-linear iterative
3093 process */
3094 mesh, alloy->solute_equation);
3095
3096 /* Update the source term for the thermal equation */
3097
3098 alloy->update_thm_st(mesh, connect, quant, time_step);
3099
3100 /* Solve the thermal system */
3101
3102 cs_thermal_system_compute(false, /* No cur2prev inside a non-linear
3103 iterative process */
3104 mesh, connect, quant, time_step);
3105
3106 /* Update fields and properties which are related to solved variables
3107 * g_l, state */
3108
3109 alloy->update_gl(mesh, connect, quant, time_step);
3110
3111 /* Update the diffusion property related to the solute */
3112
3113 if (alloy->diff_coef > cs_solidification_diffusion_eps) {
3114
3115 const double rho_D = rho0 * alloy->diff_coef;
3116
3117 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
3118 for (cs_lnum_t i = 0; i < quant->n_cells; i++)
3119 alloy->diff_pty_array[i] = (g_l[i] > 0) ?
3120 rho_D * g_l[i] : cs_solidification_diffusion_eps;
3121
3122 }
3123
3124 /* Evolution of the temperature and the bulk concentration during this
3125 iteration */
3126
3127 delta_temp = -1, delta_cbulk = -1;
3128 cs_lnum_t cid_maxt = -1, cid_maxc = -1;
3129 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
3130
3131 cs_real_t dtemp = fabs(temp[c_id] - alloy->tk_bulk[c_id]);
3132 cs_real_t dconc = fabs(conc[c_id] - alloy->ck_bulk[c_id]);
3133
3134 alloy->tk_bulk[c_id] = temp[c_id];
3135 alloy->ck_bulk[c_id] = conc[c_id];
3136
3137 if (dtemp > delta_temp)
3138 delta_temp = dtemp, cid_maxt = c_id;
3139 if (dconc > delta_cbulk)
3140 delta_cbulk = dconc, cid_maxc = c_id;
3141
3142 } /* Loop on cells */
3143
3144 alloy->iter += 1;
3145 if (solid->verbosity > 0) {
3146 cs_log_printf(CS_LOG_DEFAULT,
3147 "### Solidification.NL: "
3148 " k= %d | delta_temp= %5.3e | delta_cbulk= %5.3e\n",
3149 alloy->iter, delta_temp, delta_cbulk);
3150 if (solid->verbosity > 1)
3151 cs_log_printf(CS_LOG_DEFAULT,
3152 "### Solidification.NL: "
3153 " k= %d | delta_temp= %7ld | delta_cbulk= %7ld\n",
3154 alloy->iter, (long)cid_maxt, (long)cid_maxc);
3155 }
3156
3157 } /* while iterating */
3158
3159 /* Update the liquid concentration of the solute (c_l) */
3160
3161 alloy->update_clc(mesh, connect, quant, time_step);
3162
3163 /* The cell state is now updated at this stage. This will be useful for
3164 the monitoring */
3165
3166 _update_binary_alloy_final_state(connect, quant, time_step);
3167
3168 /* Update the forcing term in the momentum equation */
3169
3170 alloy->update_velocity_forcing(mesh, connect, quant, time_step);
3171
3172 if (solid->post_flag & CS_SOLIDIFICATION_POST_ENTHALPY)
3173 _compute_enthalpy(quant,
3174 time_step->t_cur, /* t_eval */
3175 solid->temperature->val, /* temperature */
3176 solid->g_l_field->val, /* liquid fraction */
3177 alloy->t_eut, /* temp_ref */
3178 solid->latent_heat, /* latent heat coeff. */
3179 solid->mass_density, /* rho */
3180 solid->cp, /* cp */
3181 solid->enthalpy->val); /* computed enthalpy */
3182 }
3183
3184 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
3185
3186 /*============================================================================
3187 * Public function prototypes
3188 *============================================================================*/
3189
3190 /*----------------------------------------------------------------------------*/
3191 /*!
3192 * \brief Test if solidification module is activated
3193 */
3194 /*----------------------------------------------------------------------------*/
3195
3196 bool
cs_solidification_is_activated(void)3197 cs_solidification_is_activated(void)
3198 {
3199 if (cs_solidification_structure == NULL)
3200 return false;
3201 else
3202 return true;
3203 }
3204
3205 /*----------------------------------------------------------------------------*/
3206 /*!
3207 * \brief Retrieve the main structure to deal with solidification process
3208 *
3209 * \return a pointer to a new allocated solidification structure
3210 */
3211 /*----------------------------------------------------------------------------*/
3212
3213 cs_solidification_t *
cs_solidification_get_structure(void)3214 cs_solidification_get_structure(void)
3215 {
3216 return cs_solidification_structure;
3217 }
3218
3219 /*----------------------------------------------------------------------------*/
3220 /*!
3221 * \brief Set the level of verbosity for the solidification module
3222 *
3223 * \param[in] verbosity level of verbosity to set
3224 */
3225 /*----------------------------------------------------------------------------*/
3226
3227 void
cs_solidification_set_verbosity(int verbosity)3228 cs_solidification_set_verbosity(int verbosity)
3229 {
3230 cs_solidification_t *solid = cs_solidification_structure;
3231 if (solid == NULL)
3232 return;
3233
3234 solid->verbosity = verbosity;
3235 }
3236
3237 /*----------------------------------------------------------------------------*/
3238 /*!
3239 * \brief Set the value of the epsilon parameter used in the forcing term
3240 * of the momentum equation
3241 *
3242 * \param[in] forcing_eps epsilon used in the penalization term to avoid a
3243 * division by zero
3244 */
3245 /*----------------------------------------------------------------------------*/
3246
3247 void
cs_solidification_set_forcing_eps(cs_real_t forcing_eps)3248 cs_solidification_set_forcing_eps(cs_real_t forcing_eps)
3249 {
3250 assert(forcing_eps > 0);
3251 cs_solidification_forcing_eps = forcing_eps;
3252 }
3253
3254 /*----------------------------------------------------------------------------*/
3255 /*!
3256 * \brief Activate the solidification module
3257 *
3258 * \param[in] model type of modelling
3259 * \param[in] options flag to handle optional parameters
3260 * \param[in] post_flag predefined post-processings
3261 * \param[in] boundaries pointer to the domain boundaries
3262 * \param[in] ns_model model equations for the NavSto system
3263 * \param[in] ns_model_flag option flag for the Navier-Stokes system
3264 * \param[in] algo_coupling algorithm used for solving the NavSto system
3265 * \param[in] ns_post_flag predefined post-processings for Navier-Stokes
3266 *
3267 * \return a pointer to a new allocated solidification structure
3268 */
3269 /*----------------------------------------------------------------------------*/
3270
3271 cs_solidification_t *
cs_solidification_activate(cs_solidification_model_t model,cs_flag_t options,cs_flag_t post_flag,const cs_boundary_t * boundaries,cs_navsto_param_model_t ns_model,cs_navsto_param_model_flag_t ns_model_flag,cs_navsto_param_coupling_t algo_coupling,cs_navsto_param_post_flag_t ns_post_flag)3272 cs_solidification_activate(cs_solidification_model_t model,
3273 cs_flag_t options,
3274 cs_flag_t post_flag,
3275 const cs_boundary_t *boundaries,
3276 cs_navsto_param_model_t ns_model,
3277 cs_navsto_param_model_flag_t ns_model_flag,
3278 cs_navsto_param_coupling_t algo_coupling,
3279 cs_navsto_param_post_flag_t ns_post_flag)
3280 {
3281 cs_flag_t thm_num = 0, thm_post = 0, thm_model = 0;
3282
3283 /* Allocate an empty structure */
3284
3285 cs_solidification_t *solid = _solidification_create();
3286
3287 /* Set members of the structure according to the given settings */
3288
3289 solid->model = model;
3290
3291 /* By default, Stefan model is with a frozen field */
3292
3293 if (model == CS_SOLIDIFICATION_MODEL_STEFAN)
3294 options |= CS_SOLIDIFICATION_NO_VELOCITY_FIELD;
3295 solid->options = options;
3296
3297 if (post_flag & CS_SOLIDIFICATION_ADVANCED_ANALYSIS)
3298 post_flag |= CS_SOLIDIFICATION_POST_LIQUIDUS_TEMPERATURE;
3299 solid->post_flag = post_flag;
3300
3301 /* The Navier-Stokes is not solved when the frozen field is set */
3302
3303 if (solid->options & CS_SOLIDIFICATION_NO_VELOCITY_FIELD) {
3304
3305 solid->forcing_mom = NULL;
3306 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY)
3307 bft_error(__FILE__, __LINE__, 0,
3308 "%s: Incompatible set of options: "
3309 "no velocity and binary alloy model.\n"
3310 "Please check your settings.\n", __func__);
3311
3312 }
3313 else {
3314
3315 ns_model_flag |=
3316 CS_NAVSTO_MODEL_BOUSSINESQ | CS_NAVSTO_MODEL_WITH_SOLIDIFICATION;
3317 thm_model |= CS_THERMAL_MODEL_NAVSTO_ADVECTION;
3318
3319 /* Add a property taking into account the head losses induced by the
3320 solidification process */
3321
3322 solid->forcing_mom = cs_property_add("forcing_momentum_coef",
3323 CS_PROPERTY_ISO);
3324
3325 /* If liquid, this coefficient is equal to zero */
3326
3327 cs_property_set_reference_value(solid->forcing_mom, 0);
3328
3329 }
3330
3331 /* Activate the Navier-Stokes module */
3332 /* --------------------------------- */
3333
3334 cs_navsto_system_t *ns = cs_navsto_system_activate(boundaries,
3335 ns_model,
3336 ns_model_flag,
3337 algo_coupling,
3338 ns_post_flag);
3339
3340 /* Activate and default settings for the thermal module */
3341 /* ---------------------------------------------------- */
3342
3343 if (options & CS_SOLIDIFICATION_USE_ENTHALPY_VARIABLE)
3344 thm_model |= CS_THERMAL_MODEL_USE_ENTHALPY;
3345 else
3346 thm_model |= CS_THERMAL_MODEL_USE_TEMPERATURE;
3347
3348 solid->thermal_sys = cs_thermal_system_activate(thm_model, thm_num, thm_post);
3349
3350 cs_equation_param_t *th_eqp =
3351 cs_equation_get_param(solid->thermal_sys->thermal_eq);
3352
3353 /* Be sure that the space discretization is a CDO-Fb scheme */
3354
3355 cs_equation_param_set(th_eqp, CS_EQKEY_SPACE_SCHEME, "cdofb");
3356
3357 if (thm_model & CS_THERMAL_MODEL_USE_TEMPERATURE) {
3358
3359 /* Add reaction property for the temperature equation */
3360
3361 solid->thermal_reaction_coef = cs_property_add("thermal_reaction_coef",
3362 CS_PROPERTY_ISO);
3363
3364 cs_property_set_reference_value(solid->thermal_reaction_coef, 0);
3365
3366 cs_equation_add_reaction(th_eqp, solid->thermal_reaction_coef);
3367
3368 }
3369
3370 /* Retrieve or add properties related to this module */
3371
3372 solid->mass_density = ns->param->mass_density;
3373 assert(solid->mass_density != NULL);
3374
3375 solid->viscosity = ns->param->tot_viscosity;
3376 assert(solid->viscosity != NULL);
3377
3378 solid->cp = cs_property_by_name(CS_THERMAL_CP_NAME);
3379 assert(solid->cp != NULL);
3380
3381 solid->g_l = cs_property_add("liquid_fraction", CS_PROPERTY_ISO);
3382 cs_property_set_reference_value(solid->g_l, 1.0);
3383
3384 /* Initialize other members */
3385
3386 solid->enthalpy = NULL; /* Will be created later */
3387 solid->latent_heat = 1.0; /* Default dummy value */
3388
3389 /* Allocate the structure storing the modelling context/settings */
3390
3391 switch (solid->model) {
3392
3393 case CS_SOLIDIFICATION_MODEL_STEFAN:
3394 {
3395 cs_solidification_stefan_t *s_model = NULL;
3396 BFT_MALLOC(s_model, 1, cs_solidification_stefan_t);
3397
3398 /* Default initialization of this model */
3399
3400 s_model->t_change = 0.;
3401 s_model->n_iter_max = 15;
3402 s_model->max_delta_h = 1e-2;
3403
3404 /* Function pointers */
3405
3406 solid->strategy = CS_SOLIDIFICATION_STRATEGY_PATH;
3407 s_model->update_gl = _update_gl_stefan;
3408 s_model->update_thm_st = _update_thm_stefan;
3409
3410 /* Set the context */
3411
3412 solid->model_context = (void *)s_model;
3413
3414 }
3415 break;
3416
3417 case CS_SOLIDIFICATION_MODEL_VOLLER_PRAKASH_87:
3418 case CS_SOLIDIFICATION_MODEL_VOLLER_NL:
3419 {
3420 cs_solidification_voller_t *v_model = NULL;
3421 BFT_MALLOC(v_model, 1, cs_solidification_voller_t);
3422
3423 /* Default initialization of this model */
3424
3425 v_model->s_das = 0.33541;
3426 v_model->t_solidus = 0.;
3427 v_model->t_liquidus = 1.0;
3428
3429 /* Non-linear algorithm */
3430
3431 cs_iter_algo_param_t nl_param = {
3432 .verbosity = 0, /* level of display to output */
3433 .n_max_algo_iter = 15, /* n_max iter. */
3434 .atol = 1e-6, /* absolute tolerance */
3435 .rtol = 1e-2, /* relative tolerance */
3436 .dtol = 1e3 }; /* divergence tolerance */
3437
3438 if (solid->model == CS_SOLIDIFICATION_MODEL_VOLLER_NL) {
3439
3440 v_model->nl_algo_type = CS_PARAM_NL_ALGO_PICARD;
3441 v_model->nl_algo = cs_iter_algo_create(nl_param);
3442
3443 }
3444 else {
3445
3446 v_model->nl_algo_type = CS_PARAM_N_NL_ALGOS;
3447 v_model->nl_algo = NULL;
3448
3449 }
3450
3451 /* Function pointers */
3452
3453 solid->strategy = CS_SOLIDIFICATION_STRATEGY_LEGACY;
3454 if (solid->options & CS_SOLIDIFICATION_NO_VELOCITY_FIELD)
3455 v_model->update_gl = _update_gl_voller_legacy_no_velocity;
3456 else
3457 v_model->update_gl = _update_gl_voller_legacy;
3458 v_model->update_thm_st = _update_thm_voller_legacy;
3459
3460 /* If the non-linear model is used, then the default strategy is
3461 modified */
3462
3463 if (solid->model == CS_SOLIDIFICATION_MODEL_VOLLER_NL) {
3464 solid->strategy = CS_SOLIDIFICATION_STRATEGY_PATH;
3465 v_model->update_thm_st = _update_thm_voller_path;
3466 }
3467
3468 /* Set the context */
3469 solid->model_context = (void *)v_model;
3470
3471 }
3472 break;
3473
3474 case CS_SOLIDIFICATION_MODEL_BINARY_ALLOY:
3475 {
3476 cs_solidification_binary_alloy_t *b_model = NULL;
3477 BFT_MALLOC(b_model, 1, cs_solidification_binary_alloy_t);
3478
3479 /* Set a default value to the model parameters */
3480
3481 b_model->ref_concentration = 1.;
3482 b_model->s_das = 0.33541;
3483 b_model->kp = 0.5;
3484 b_model->inv_kp = 2.;
3485 b_model->inv_kpm1 = -2;
3486 b_model->ml = -1;
3487 b_model->inv_ml = -1;
3488 b_model->t_melt = 0.;
3489 b_model->t_eut = b_model->t_eut_inf = b_model->t_eut_sup = 0.;
3490 b_model->c_eut = b_model->cs1 = b_model->dgldC_eut = 0.;
3491
3492 b_model->diff_coef = 1e-6;
3493
3494 /* Monitoring and criteria to drive the convergence of the coupled system
3495 (solute transport and thermal equation) */
3496
3497 b_model->iter = 0;
3498 b_model->n_iter_max = 10;
3499 b_model->delta_tolerance = 1e-3;
3500 b_model->eta_relax = 0.;
3501 b_model->gliq_relax = 0.;
3502
3503 /* Strategy to perform the main steps of the simulation of a binary alloy
3504 * Default strategy: Taylor which corresponds to the Legacy one with
3505 * improvements thanks to some Taylor expansions */
3506
3507 solid->strategy = CS_SOLIDIFICATION_STRATEGY_TAYLOR;
3508
3509 /* Functions which are specific to a strategy */
3510
3511 b_model->update_gl = _update_gl_taylor;
3512 b_model->update_thm_st = _update_thm_taylor;
3513
3514 /* Functions which are common to all strategies */
3515
3516 b_model->update_velocity_forcing = _update_velocity_forcing;
3517 b_model->update_clc = _update_clc;
3518
3519 /* Initialize pointers */
3520
3521 b_model->tk_bulk = NULL;
3522 b_model->ck_bulk = NULL;
3523 b_model->tx_bulk = NULL;
3524 b_model->cx_bulk = NULL;
3525
3526 /* alloy->temp_faces is shared and set when defined */
3527
3528 b_model->c_l_cells = NULL;
3529 b_model->c_l_faces = NULL;
3530
3531 b_model->solute_equation = NULL;
3532 b_model->c_bulk = NULL;
3533
3534 b_model->diff_pty = NULL;
3535 b_model->diff_pty_array = NULL;
3536
3537 /* eta_coef is activated with advanced options */
3538
3539 b_model->eta_coef_pty = NULL;
3540 b_model->eta_coef_array = NULL;
3541
3542 /* Optional post-processing arrays */
3543
3544 b_model->t_liquidus = NULL;
3545 b_model->tbulk_minus_tliq = NULL;
3546 b_model->cliq_minus_cbulk = NULL;
3547
3548 /* Set the context */
3549
3550 solid->model_context = (void *)b_model;
3551 }
3552 break;
3553
3554 default:
3555 bft_error(__FILE__, __LINE__, 0,
3556 " %s: Invalid model for the solidification module.\n"
3557 " Please check your setup.", __func__);
3558
3559 } /* Switch on the solidification model */
3560
3561 /* Set the global pointer */
3562
3563 cs_solidification_structure = solid;
3564
3565 return solid;
3566 }
3567
3568 /*----------------------------------------------------------------------------*/
3569 /*!
3570 * \brief Get the structure defining the Stefan model
3571 *
3572 * \return a pointer to the structure
3573 */
3574 /*----------------------------------------------------------------------------*/
3575
3576 cs_solidification_stefan_t *
cs_solidification_get_stefan_struct(void)3577 cs_solidification_get_stefan_struct(void)
3578 {
3579 cs_solidification_t *solid = cs_solidification_structure;
3580
3581 /* Sanity checks */
3582 if (solid == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_module));
3583
3584 if (solid->model != CS_SOLIDIFICATION_MODEL_STEFAN)
3585 bft_error(__FILE__, __LINE__, 0,
3586 " %s: Stefan model not declared during the"
3587 " activation of the solidification module.\n"
3588 " Please check your settings.", __func__);
3589
3590 cs_solidification_stefan_t *s_model
3591 = (cs_solidification_stefan_t *)solid->model_context;
3592 assert(s_model != NULL);
3593
3594 return s_model;
3595 }
3596
3597 /*----------------------------------------------------------------------------*/
3598 /*!
3599 * \brief Sanity checks on the consistency of the Stefan's model settings
3600 *
3601 * \return a pointer to the structure
3602 */
3603 /*----------------------------------------------------------------------------*/
3604
3605 cs_solidification_stefan_t *
cs_solidification_check_stefan_model(void)3606 cs_solidification_check_stefan_model(void)
3607 {
3608 cs_solidification_stefan_t *s_model = cs_solidification_get_stefan_struct();
3609
3610 if (s_model->n_iter_max < 1)
3611 bft_error(__FILE__, __LINE__, 0,
3612 "%s: Invalid value for n_iter_max (= %d)\n",
3613 __func__, s_model->n_iter_max);
3614 if (s_model->max_delta_h < FLT_MIN)
3615 bft_error(__FILE__, __LINE__, 0,
3616 "%s: Invalid value for max_delta_h (= %5.3e)\n",
3617 __func__, s_model->max_delta_h);
3618
3619 return s_model;
3620 }
3621
3622 /*----------------------------------------------------------------------------*/
3623 /*!
3624 * \brief Set the main physical parameters which describe the Stefan model
3625 *
3626 * \param[in] t_change liquidus/solidus temperature (in K)
3627 * \param[in] latent_heat latent heat
3628 */
3629 /*----------------------------------------------------------------------------*/
3630
3631 void
cs_solidification_set_stefan_model(cs_real_t t_change,cs_real_t latent_heat)3632 cs_solidification_set_stefan_model(cs_real_t t_change,
3633 cs_real_t latent_heat)
3634 {
3635 cs_solidification_t *solid = cs_solidification_structure;
3636 cs_solidification_stefan_t *s_model = cs_solidification_get_stefan_struct();
3637
3638 /* Model parameters */
3639 s_model->t_change = t_change;
3640 solid->latent_heat = latent_heat;
3641 }
3642
3643 /*----------------------------------------------------------------------------*/
3644 /*!
3645 * \brief Get the structure defining the Voller model
3646 *
3647 * \return a pointer to the structure
3648 */
3649 /*----------------------------------------------------------------------------*/
3650
3651 cs_solidification_voller_t *
cs_solidification_get_voller_struct(void)3652 cs_solidification_get_voller_struct(void)
3653 {
3654 cs_solidification_t *solid = cs_solidification_structure;
3655
3656 /* Sanity checks */
3657 if (solid == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_module));
3658
3659 if (solid->model != CS_SOLIDIFICATION_MODEL_VOLLER_PRAKASH_87 &&
3660 solid->model != CS_SOLIDIFICATION_MODEL_VOLLER_NL )
3661 bft_error(__FILE__, __LINE__, 0,
3662 " %s: Voller model not declared during the"
3663 " activation of the solidification module.\n"
3664 " Please check your settings.", __func__);
3665
3666 cs_solidification_voller_t *v_model
3667 = (cs_solidification_voller_t *)solid->model_context;
3668 assert(v_model != NULL);
3669
3670 return v_model;
3671 }
3672
3673 /*----------------------------------------------------------------------------*/
3674 /*!
3675 * \brief Sanity checks on the consistency of the Voller's model settings
3676 *
3677 * \return a pointer to the structure
3678 */
3679 /*----------------------------------------------------------------------------*/
3680
3681 cs_solidification_voller_t *
cs_solidification_check_voller_model(void)3682 cs_solidification_check_voller_model(void)
3683 {
3684 cs_solidification_voller_t *v_model = cs_solidification_get_voller_struct();
3685
3686 if (v_model->t_liquidus - v_model->t_solidus < 0.)
3687 bft_error(__FILE__, __LINE__, 0,
3688 " %s: The liquidus and solidus temperatures are not"
3689 " consistent.\n"
3690 " Please check your settings.", __func__);
3691 if (v_model->s_das < FLT_MIN)
3692 bft_error(__FILE__, __LINE__, 0,
3693 " %s: Invalid value %g for the secondary dendrite arms spacing",
3694 __func__, v_model->s_das);
3695
3696 return v_model;
3697 }
3698
3699 /*----------------------------------------------------------------------------*/
3700 /*!
3701 * \brief Set the main physical parameters which describe the Voller and
3702 * Prakash modelling
3703 *
3704 * \param[in] beta thermal dilatation coefficient
3705 * \param[in] t_ref reference temperature (for the Boussinesq approx)
3706 * \param[in] t_solidus solidus temperature (in K)
3707 * \param[in] t_liquidus liquidus temperature (in K)
3708 * \param[in] latent_heat latent heat
3709 * \param[in] s_das secondary dendrite space arms
3710 */
3711 /*----------------------------------------------------------------------------*/
3712
3713 void
cs_solidification_set_voller_model(cs_real_t beta,cs_real_t t_ref,cs_real_t t_solidus,cs_real_t t_liquidus,cs_real_t latent_heat,cs_real_t s_das)3714 cs_solidification_set_voller_model(cs_real_t beta,
3715 cs_real_t t_ref,
3716 cs_real_t t_solidus,
3717 cs_real_t t_liquidus,
3718 cs_real_t latent_heat,
3719 cs_real_t s_das)
3720 {
3721 cs_solidification_t *solid = cs_solidification_structure;
3722 cs_solidification_voller_t *v_model = cs_solidification_get_voller_struct();
3723
3724 /* Model parameters */
3725 v_model->t_solidus = t_solidus;
3726 v_model->t_liquidus = t_liquidus;
3727 v_model->s_das = s_das;
3728
3729 solid->latent_heat = latent_heat;
3730
3731 /* Add the Boussinesq term */
3732 cs_navsto_param_t *nsp = cs_navsto_system_get_param();
3733
3734 cs_navsto_param_add_boussinesq_term(nsp, beta, t_ref);
3735
3736 /* Set the reference temperature in the thermal module */
3737 cs_thermal_system_set_reference_temperature(t_ref);
3738 }
3739
3740 /*----------------------------------------------------------------------------*/
3741 /*!
3742 * \brief Set the main physical parameters which describe the Voller and
3743 * Prakash modelling
3744 *
3745 * \param[in] beta thermal dilatation coefficient
3746 * \param[in] t_ref reference temperature (for the Boussinesq approx)
3747 * \param[in] t_solidus solidus temperature (in K)
3748 * \param[in] t_liquidus liquidus temperature (in K)
3749 * \param[in] latent_heat latent heat
3750 */
3751 /*----------------------------------------------------------------------------*/
3752
3753 void
cs_solidification_set_voller_model_no_velocity(cs_real_t t_solidus,cs_real_t t_liquidus,cs_real_t latent_heat)3754 cs_solidification_set_voller_model_no_velocity(cs_real_t t_solidus,
3755 cs_real_t t_liquidus,
3756 cs_real_t latent_heat)
3757 {
3758 cs_solidification_t *solid = cs_solidification_structure;
3759 cs_solidification_voller_t *v_model = cs_solidification_get_voller_struct();
3760
3761 if ((solid->options & CS_SOLIDIFICATION_NO_VELOCITY_FIELD) == 0)
3762 bft_error(__FILE__, __LINE__, 0,
3763 "%s: CS_SOLIDIFICATION_NO_VELOCITY_FIELD has not been set with"
3764 " the Voller model.\n"
3765 "Please check your settings.\n", __func__);
3766
3767 /* Model parameters (those which are useful in this case) */
3768
3769 v_model->t_solidus = t_solidus;
3770 v_model->t_liquidus = t_liquidus;
3771 solid->latent_heat = latent_heat;
3772
3773 }
3774
3775 /*----------------------------------------------------------------------------*/
3776 /*!
3777 * \brief Get the structure defining the binary alloy model
3778 *
3779 * \return a pointer to the structure
3780 */
3781 /*----------------------------------------------------------------------------*/
3782
3783 cs_solidification_binary_alloy_t *
cs_solidification_get_binary_alloy_struct(void)3784 cs_solidification_get_binary_alloy_struct(void)
3785 {
3786 cs_solidification_t *solid = cs_solidification_structure;
3787
3788 /* Sanity checks */
3789 if (solid == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_module));
3790
3791 if (solid->model != CS_SOLIDIFICATION_MODEL_BINARY_ALLOY)
3792 bft_error(__FILE__, __LINE__, 0,
3793 " %s: The binary alloy model was not declared during the"
3794 " activation of the solidification module.\n"
3795 " Please check your settings.", __func__);
3796
3797 cs_solidification_binary_alloy_t *b_model
3798 = (cs_solidification_binary_alloy_t *)solid->model_context;
3799 assert(b_model != NULL);
3800
3801 return b_model;
3802 }
3803
3804 /*----------------------------------------------------------------------------*/
3805 /*!
3806 * \brief Sanity checks on the consistency of the settings of the binary alloy
3807 * model
3808 *
3809 * \return a pointer to the structure
3810 */
3811 /*----------------------------------------------------------------------------*/
3812
3813 cs_solidification_binary_alloy_t *
cs_solidification_check_binary_alloy_model(void)3814 cs_solidification_check_binary_alloy_model(void)
3815 {
3816 cs_solidification_binary_alloy_t
3817 *b_model = cs_solidification_get_binary_alloy_struct();
3818
3819 if (b_model->s_das < FLT_MIN)
3820 bft_error(__FILE__, __LINE__, 0,
3821 " %s: Invalid value %g for the secondary dendrite arms spacing",
3822 __func__, b_model->s_das);
3823 if (b_model->kp < FLT_MIN || b_model->kp > 1 - FLT_MIN)
3824 bft_error(__FILE__, __LINE__, 0,
3825 " %s: Invalid value %g for partition coefficient",
3826 __func__, b_model->kp);
3827 if (fabs(b_model->ml) < FLT_MIN)
3828 bft_error(__FILE__, __LINE__, 0,
3829 " %s: Invalid value %g for the liquidus slope",
3830 __func__, b_model->ml);
3831 if (b_model->n_iter_max < 1)
3832 bft_error(__FILE__, __LINE__, 0,
3833 "%s: Invalid value for n_iter_max (current: %d).\n"
3834 " Should be strictly greater than 0.\n",
3835 __func__, b_model->n_iter_max);
3836 if (b_model->delta_tolerance < FLT_MIN)
3837 bft_error(__FILE__, __LINE__, 0,
3838 "%s: Invalid value for \"tolerance\" (current: %5.3e).\n",
3839 __func__, b_model->delta_tolerance);
3840
3841 return b_model;
3842 }
3843
3844 /*----------------------------------------------------------------------------*/
3845 /*!
3846 * \brief Set the main physical parameters which describe a solidification
3847 * process with a binary alloy (with components A and B)
3848 * Add a transport equation for the solute concentration to simulate
3849 * the conv/diffusion of the alloy ratio between the two components of
3850 * the alloy
3851 *
3852 * \param[in] name name of the binary alloy
3853 * \param[in] varname name of the unknown related to the tracer eq.
3854 * \param[in] beta_t thermal dilatation coefficient
3855 * \param[in] temp0 reference temperature (Boussinesq term)
3856 * \param[in] beta_c solutal dilatation coefficient
3857 * \param[in] conc0 reference mixture concentration (Boussinesq term)
3858 * \param[in] kp value of the distribution coefficient
3859 * \param[in] mliq liquidus slope for the solute concentration
3860 * \param[in] t_eutec temperature at the eutectic point
3861 * \param[in] t_melt phase-change temperature for the pure material (A)
3862 * \param[in] solute_diff solutal diffusion coefficient in the liquid
3863 * \param[in] latent_heat latent heat
3864 * \param[in] s_das secondary dendrite arm spacing
3865 */
3866 /*----------------------------------------------------------------------------*/
3867
3868 void
cs_solidification_set_binary_alloy_model(const char * name,const char * varname,cs_real_t beta_t,cs_real_t temp0,cs_real_t beta_c,cs_real_t conc0,cs_real_t kp,cs_real_t mliq,cs_real_t t_eutec,cs_real_t t_melt,cs_real_t solute_diff,cs_real_t latent_heat,cs_real_t s_das)3869 cs_solidification_set_binary_alloy_model(const char *name,
3870 const char *varname,
3871 cs_real_t beta_t,
3872 cs_real_t temp0,
3873 cs_real_t beta_c,
3874 cs_real_t conc0,
3875 cs_real_t kp,
3876 cs_real_t mliq,
3877 cs_real_t t_eutec,
3878 cs_real_t t_melt,
3879 cs_real_t solute_diff,
3880 cs_real_t latent_heat,
3881 cs_real_t s_das)
3882 {
3883 /* Check the validity of some parameters */
3884 if (kp < FLT_MIN || kp > 1 - FLT_MIN)
3885 bft_error(__FILE__, __LINE__, 0,
3886 " %s: Invalid value %g for partition coefficient", __func__, kp);
3887 if (fabs(mliq) < FLT_MIN)
3888 bft_error(__FILE__, __LINE__, 0,
3889 " %s: Invalid value %g for the liquidus slope", __func__, mliq);
3890 if (s_das < FLT_MIN)
3891 bft_error(__FILE__, __LINE__, 0,
3892 " %s: Invalid value %g for the secondary dendrite arms spacing",
3893 __func__, s_das);
3894
3895 /* Retrieve and set the binary alloy structures */
3896 cs_solidification_t *solid = cs_solidification_structure;
3897 cs_solidification_binary_alloy_t
3898 *alloy = cs_solidification_get_binary_alloy_struct();
3899
3900 solid->latent_heat = latent_heat;
3901
3902 /* Set the main physical parameters/constants */
3903 alloy->ref_concentration = conc0;
3904 alloy->s_das = s_das;
3905
3906 /* Phase diagram parameters and related quantities */
3907 alloy->kp = kp; /* partition coeff. */
3908 alloy->ml = mliq; /* liquidus slope */
3909 alloy->t_melt = t_melt; /* melting temperature */
3910 alloy->t_eut = t_eutec; /* eutectic temperature */
3911
3912 /* Derived quantities from the phase diagram */
3913 alloy->inv_kp = 1./kp;
3914 alloy->inv_kpm1 = 1./(alloy->kp - 1.);
3915 alloy->inv_ml = 1./mliq;
3916 alloy->c_eut = (t_eutec - t_melt)*alloy->inv_ml;
3917 alloy->cs1 = alloy->c_eut * kp; /* Apply the lever rule */
3918
3919 assert(fabs(alloy->c_eut - alloy->cs1) > FLT_MIN);
3920 alloy->dgldC_eut = 1./(alloy->c_eut - alloy->cs1);
3921
3922 /* Define a small range of temperature around the eutectic temperature
3923 * in which one assumes an eutectic transformation */
3924 alloy->t_eut_inf =
3925 alloy->t_eut - cs_solidification_eutectic_threshold;
3926 alloy->t_eut_sup =
3927 alloy->t_eut + cs_solidification_eutectic_threshold;
3928
3929 /* Alloy equation and variable field */
3930 assert(name != NULL && varname != NULL);
3931 cs_equation_t *eq = cs_equation_add(name, varname,
3932 CS_EQUATION_TYPE_SOLIDIFICATION,
3933 1,
3934 CS_PARAM_BC_HMG_NEUMANN);
3935
3936 /* Set an upwind scheme by default since it could be a pure advection eq. */
3937 cs_equation_param_t *eqp = cs_equation_get_param(eq);
3938
3939 /* Set the default numerical options that should be used */
3940 cs_equation_param_set(eqp, CS_EQKEY_SPACE_SCHEME, "cdo_fb");
3941 cs_equation_param_set(eqp, CS_EQKEY_HODGE_DIFF_ALGO, "cost");
3942 cs_equation_param_set(eqp, CS_EQKEY_HODGE_DIFF_COEF, "sushi");
3943 cs_equation_param_set(eqp, CS_EQKEY_ADV_SCHEME, "upwind");
3944 cs_equation_param_set(eqp, CS_EQKEY_ADV_FORMULATION, "conservative");
3945
3946 alloy->solute_equation = eq;
3947 alloy->c_bulk = NULL; /* Variable field related to this equation. This will
3948 be set later (after the equation initialization) */
3949
3950 /* Always add a diffusion term (to avoid a zero block face-face when there
3951 is no more convection */
3952 if (solute_diff > 0)
3953 alloy->diff_coef = solute_diff;
3954 else
3955 alloy->diff_coef = cs_solidification_diffusion_eps;
3956
3957 char *pty_name = NULL;
3958 size_t len = strlen(varname) + strlen("_diff_pty");
3959 BFT_MALLOC(pty_name, len + 1, char);
3960 sprintf(pty_name, "%s_diff_pty", varname);
3961 pty_name[len] = '\0';
3962 alloy->diff_pty = cs_property_add(pty_name, CS_PROPERTY_ISO);
3963 BFT_FREE(pty_name);
3964
3965 cs_equation_add_diffusion(eqp, alloy->diff_pty);
3966
3967 /* Add Boussinesq terms */
3968 cs_navsto_param_t *nsp = cs_navsto_system_get_param();
3969
3970 /* Thermal effect */
3971 cs_navsto_param_add_boussinesq_term(nsp, beta_t, temp0);
3972
3973 /* Solutal effect */
3974 cs_navsto_param_add_boussinesq_term(nsp, beta_c, conc0);
3975
3976 /* Set the reference temperature in the thermal module */
3977 cs_thermal_system_set_reference_temperature(temp0);
3978
3979 }
3980
3981 /*----------------------------------------------------------------------------*/
3982 /*!
3983 * \brief Set the strategy to update quantitiess (liquid fraction and
3984 * the thermal source term for the two main quantities)
3985 *
3986 * \param[in] strategy strategy to perform the update of quantities
3987 */
3988 /*----------------------------------------------------------------------------*/
3989
3990 void
cs_solidification_set_strategy(cs_solidification_strategy_t strategy)3991 cs_solidification_set_strategy(cs_solidification_strategy_t strategy)
3992 {
3993 cs_solidification_t *solid = cs_solidification_structure;
3994
3995 switch (solid->model) {
3996
3997 case CS_SOLIDIFICATION_MODEL_STEFAN:
3998 cs_base_warn(__FILE__, __LINE__);
3999 bft_printf("%s: Only one strategy is available with the Stefan model.\n",
4000 __func__);
4001 break;
4002
4003 case CS_SOLIDIFICATION_MODEL_VOLLER_PRAKASH_87:
4004 case CS_SOLIDIFICATION_MODEL_VOLLER_NL:
4005 {
4006 cs_solidification_voller_t *v_model =
4007 cs_solidification_get_voller_struct();
4008
4009 switch (strategy) {
4010
4011 case CS_SOLIDIFICATION_STRATEGY_LEGACY:
4012 v_model->update_thm_st = _update_thm_voller_legacy;
4013 break;
4014
4015 case CS_SOLIDIFICATION_STRATEGY_PATH:
4016 v_model->update_thm_st = _update_thm_voller_path;
4017 break;
4018
4019 default:
4020 if (solid->model == CS_SOLIDIFICATION_MODEL_VOLLER_PRAKASH_87)
4021 v_model->update_thm_st = _update_thm_voller_legacy;
4022 else
4023 v_model->update_thm_st = _update_thm_voller_path;
4024 break;
4025
4026 } /* Switch on the strategy */
4027
4028 }
4029 break; /* Voller-like models */
4030
4031 case CS_SOLIDIFICATION_MODEL_BINARY_ALLOY:
4032 {
4033 cs_solidification_binary_alloy_t *alloy =
4034 cs_solidification_get_binary_alloy_struct();
4035
4036 switch (strategy) {
4037
4038 case CS_SOLIDIFICATION_STRATEGY_LEGACY:
4039 if (solid->options & CS_SOLIDIFICATION_WITH_SOLUTE_SOURCE_TERM)
4040 alloy->update_gl = _update_gl_legacy_ast;
4041 else
4042 alloy->update_gl = _update_gl_legacy;
4043 alloy->update_thm_st = _update_thm_legacy;
4044 break;
4045
4046 case CS_SOLIDIFICATION_STRATEGY_TAYLOR:
4047 if (solid->options & CS_SOLIDIFICATION_WITH_SOLUTE_SOURCE_TERM)
4048 bft_error(__FILE__, __LINE__, 0,
4049 "%s: Adding an advective source term is incompatible with"
4050 " the Taylor strategy.\n", __func__);
4051 else
4052 alloy->update_gl = _update_gl_taylor;
4053 alloy->update_thm_st = _update_thm_taylor;
4054 break;
4055
4056 case CS_SOLIDIFICATION_STRATEGY_PATH:
4057 if (solid->options & CS_SOLIDIFICATION_WITH_SOLUTE_SOURCE_TERM)
4058 bft_error(__FILE__, __LINE__, 0,
4059 "%s: Adding an advective source term is incompatible with"
4060 " the Path strategy.\n", __func__);
4061 else
4062 alloy->update_gl = _update_gl_binary_path;
4063 alloy->update_thm_st = _update_thm_binary_path;
4064 break;
4065
4066 default:
4067 bft_error(__FILE__, __LINE__, 0, "%s: Invalid strategy.\n", __func__);
4068 break;
4069
4070 } /* Switch on strategies */
4071
4072 }
4073 break;
4074
4075 default:
4076 bft_error(__FILE__, __LINE__, 0,
4077 "%s: Invalid solidification model.\n", __func__);
4078
4079 } /* Switch on the solidification model */
4080
4081 solid->strategy = strategy;
4082 }
4083
4084 /*----------------------------------------------------------------------------*/
4085 /*!
4086 * \brief Set the functions to perform the update of physical properties
4087 * and/or the computation of the thermal source term or quantities
4088 * and/or the way to perform the coupling between the thermal equation
4089 * and the bulk concentration computation. All this setting defines
4090 * the way to compute the solidification process of a binary alloy.
4091 * If a function is set to NULL then the automatic settings are kept.
4092 *
4093 * --Advanced usage-- This enables to finely control the numerical or
4094 * physical modelling aspects.
4095 *
4096 * \param[in] vel_forcing pointer to update the velocity forcing
4097 * \param[in] cliq_update pointer to update the liquid concentration
4098 * \param[in] gliq_update pointer to update the liquid fraction
4099 * \param[in] thm_st_update pointer to update thermal source terms
4100 */
4101 /*----------------------------------------------------------------------------*/
4102
4103 void
cs_solidification_set_segr_functions(cs_solidification_func_t * vel_forcing,cs_solidification_func_t * cliq_update,cs_solidification_func_t * gliq_update,cs_solidification_func_t * thm_st_update)4104 cs_solidification_set_segr_functions(cs_solidification_func_t *vel_forcing,
4105 cs_solidification_func_t *cliq_update,
4106 cs_solidification_func_t *gliq_update,
4107 cs_solidification_func_t *thm_st_update)
4108 {
4109 cs_solidification_t *solid = cs_solidification_structure;
4110 if (solid == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_module));
4111
4112 cs_solidification_binary_alloy_t *alloy
4113 = (cs_solidification_binary_alloy_t *)solid->model_context;
4114
4115 assert(solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY);
4116 assert(alloy != NULL);
4117
4118 if (vel_forcing != NULL) {
4119 alloy->update_velocity_forcing = vel_forcing;
4120 solid->options |= CS_SOLIDIFICATION_BINARY_ALLOY_M_FUNC;
4121 }
4122
4123 if (cliq_update != NULL) {
4124 alloy->update_clc = cliq_update;
4125 solid->options |= CS_SOLIDIFICATION_BINARY_ALLOY_C_FUNC;
4126 }
4127
4128 if (gliq_update != NULL) {
4129 alloy->update_gl = gliq_update;
4130 solid->options |= CS_SOLIDIFICATION_BINARY_ALLOY_G_FUNC;
4131 }
4132
4133 if (thm_st_update != NULL) {
4134 alloy->update_thm_st = thm_st_update;
4135 solid->options |= CS_SOLIDIFICATION_BINARY_ALLOY_T_FUNC;
4136 }
4137
4138 }
4139
4140 /*----------------------------------------------------------------------------*/
4141 /*!
4142 * \brief Free the main structure related to the solidification module
4143 *
4144 * \return a NULL pointer
4145 */
4146 /*----------------------------------------------------------------------------*/
4147
4148 cs_solidification_t *
cs_solidification_destroy_all(void)4149 cs_solidification_destroy_all(void)
4150 {
4151 if (cs_solidification_structure == NULL)
4152 return NULL;
4153
4154 cs_solidification_t *solid = cs_solidification_structure;
4155
4156 /* The lifecycle of properties, equations and fields is not managed by
4157 * the current structure and sub-structures.
4158 * Free only what is owned by this structure */
4159
4160 switch (solid->model) {
4161
4162 case CS_SOLIDIFICATION_MODEL_STEFAN:
4163 {
4164 cs_solidification_stefan_t *s_model
4165 = (cs_solidification_stefan_t *)solid->model_context;
4166
4167 BFT_FREE(s_model);
4168
4169 } /* Stefan modelling */
4170 break;
4171
4172 case CS_SOLIDIFICATION_MODEL_VOLLER_PRAKASH_87:
4173 case CS_SOLIDIFICATION_MODEL_VOLLER_NL:
4174 {
4175 cs_solidification_voller_t *v_model
4176 = (cs_solidification_voller_t *)solid->model_context;
4177
4178 if (CS_SOLIDIFICATION_MODEL_VOLLER_NL)
4179 BFT_FREE(v_model->nl_algo);
4180
4181 BFT_FREE(v_model);
4182
4183 } /* Voller and Prakash modelling */
4184 break;
4185
4186 case CS_SOLIDIFICATION_MODEL_BINARY_ALLOY:
4187 {
4188 cs_solidification_binary_alloy_t *alloy
4189 = (cs_solidification_binary_alloy_t *)solid->model_context;
4190
4191 BFT_FREE(alloy->diff_pty_array);
4192 BFT_FREE(alloy->c_l_cells);
4193 BFT_FREE(alloy->eta_coef_array);
4194 BFT_FREE(alloy->tk_bulk);
4195 BFT_FREE(alloy->ck_bulk);
4196
4197 if (solid->options & CS_SOLIDIFICATION_USE_EXTRAPOLATION) {
4198 BFT_FREE(alloy->tx_bulk);
4199 BFT_FREE(alloy->cx_bulk);
4200 }
4201
4202 if (solid->options & CS_SOLIDIFICATION_WITH_SOLUTE_SOURCE_TERM)
4203 BFT_FREE(alloy->c_l_faces);
4204
4205 if (solid->post_flag & CS_SOLIDIFICATION_POST_LIQUIDUS_TEMPERATURE)
4206 BFT_FREE(alloy->t_liquidus);
4207
4208 if (solid->post_flag & CS_SOLIDIFICATION_ADVANCED_ANALYSIS) {
4209 BFT_FREE(alloy->tbulk_minus_tliq);
4210 BFT_FREE(alloy->cliq_minus_cbulk);
4211 }
4212
4213 BFT_FREE(alloy);
4214
4215 } /* Binary alloy modelling */
4216 break;
4217
4218 default:
4219 break; /* Nothing to do */
4220
4221 } /* Switch on solidification model */
4222
4223 BFT_FREE(solid->thermal_reaction_coef_array);
4224 BFT_FREE(solid->thermal_source_term_array);
4225 BFT_FREE(solid->forcing_mom_array);
4226
4227 BFT_FREE(solid->cell_state);
4228
4229 if (solid->plot_state != NULL)
4230 cs_time_plot_finalize(&solid->plot_state);
4231
4232 BFT_FREE(solid);
4233
4234 return NULL;
4235 }
4236
4237 /*----------------------------------------------------------------------------*/
4238 /*!
4239 * \brief Setup equations/properties related to the solidification module
4240 */
4241 /*----------------------------------------------------------------------------*/
4242
4243 void
cs_solidification_init_setup(void)4244 cs_solidification_init_setup(void)
4245 {
4246 cs_solidification_t *solid = cs_solidification_structure;
4247
4248 /* Sanity checks */
4249
4250 if (solid == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_module));
4251
4252 const int field_mask = CS_FIELD_INTENSIVE | CS_FIELD_CDO;
4253 const int log_key = cs_field_key_id("log");
4254 const int post_key = cs_field_key_id("post_vis");
4255 const int c_loc_id = cs_mesh_location_get_id_by_name("cells");
4256
4257 /* Add a field for the liquid fraction */
4258
4259 solid->g_l_field = cs_field_create("liquid_fraction",
4260 field_mask,
4261 c_loc_id,
4262 1,
4263 true); /* has_previous */
4264
4265 cs_field_set_key_int(solid->g_l_field, log_key, 1);
4266 cs_field_set_key_int(solid->g_l_field, post_key, 1);
4267
4268 /* Add the enthalpy field if not already created */
4269
4270 solid->enthalpy = cs_field_by_name_try("enthalpy");
4271 if (solid->enthalpy == NULL) {
4272
4273 bool add_enthalpy = false;
4274
4275 if (solid->post_flag & CS_SOLIDIFICATION_POST_ENTHALPY)
4276 add_enthalpy = true;
4277 if (solid->model == CS_SOLIDIFICATION_MODEL_VOLLER_NL ||
4278 solid->model == CS_SOLIDIFICATION_MODEL_STEFAN)
4279 add_enthalpy = true;
4280
4281 if (add_enthalpy) {
4282 solid->enthalpy = cs_field_create("enthalpy",
4283 field_mask,
4284 c_loc_id,
4285 1,
4286 true); /* has_previous */
4287
4288 cs_field_set_key_int(solid->enthalpy, log_key, 1);
4289 }
4290
4291 }
4292
4293 if (solid->post_flag & CS_SOLIDIFICATION_POST_ENTHALPY)
4294 cs_field_set_key_int(solid->enthalpy, post_key, 1);
4295
4296 /* Add a reaction term to the momentum equation */
4297
4298 if (solid->forcing_mom != NULL) {
4299
4300 cs_equation_t *mom_eq = cs_navsto_system_get_momentum_eq();
4301 cs_equation_param_t *mom_eqp = cs_equation_get_param(mom_eq);
4302 assert(mom_eqp != NULL);
4303
4304 cs_equation_add_reaction(mom_eqp, solid->forcing_mom);
4305
4306 }
4307
4308 /* Add default post-processing related to the solidification module */
4309
4310 cs_post_add_time_mesh_dep_output(cs_solidification_extra_post, solid);
4311
4312 /* Model-specific part */
4313
4314 switch (solid->model) {
4315
4316 case CS_SOLIDIFICATION_MODEL_VOLLER_PRAKASH_87:
4317 case CS_SOLIDIFICATION_MODEL_VOLLER_NL:
4318 {
4319 /* Check the sanity of the model parameters and retrieve the structure */
4320
4321 cs_solidification_voller_t
4322 *v_model = cs_solidification_check_voller_model();
4323
4324 solid->forcing_coef = 180./(v_model->s_das*v_model->s_das);
4325 }
4326 break;
4327
4328 case CS_SOLIDIFICATION_MODEL_BINARY_ALLOY:
4329 {
4330 /* Check the sanity of the model parameters and retrieve the structure */
4331
4332 cs_solidification_binary_alloy_t
4333 *alloy = cs_solidification_check_binary_alloy_model();
4334
4335 cs_equation_param_t *eqp = cs_equation_get_param(alloy->solute_equation);
4336
4337 /* Add the unsteady term */
4338
4339 cs_equation_add_time(eqp, solid->mass_density);
4340
4341 /* Add an advection term to the solute concentration equation */
4342
4343 cs_equation_add_advection(eqp, cs_navsto_get_adv_field());
4344
4345 if ((solid->options & CS_SOLIDIFICATION_WITH_SOLUTE_SOURCE_TERM) == 0) {
4346
4347 alloy->eta_coef_pty = cs_property_add("alloy_adv_coef",
4348 CS_PROPERTY_ISO);
4349
4350 cs_equation_add_advection_scaling_property(eqp, alloy->eta_coef_pty);
4351
4352 }
4353
4354 solid->forcing_coef = 180./(alloy->s_das*alloy->s_das);
4355 }
4356 break; /* Binary alloy model */
4357
4358 default: /* Stefan: There is nothing else to do */
4359 break;
4360
4361 } /* Switch on model */
4362
4363 if (cs_glob_rank_id < 1) {
4364
4365 int n_output_states = CS_SOLIDIFICATION_N_STATES - 1;
4366 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY)
4367 n_output_states += 1;
4368
4369 int n_output_values = n_output_states;
4370 if (solid->post_flag & CS_SOLIDIFICATION_POST_SOLIDIFICATION_RATE)
4371 n_output_values += 1;
4372
4373 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY) {
4374 if (solid->post_flag & CS_SOLIDIFICATION_POST_SEGREGATION_INDEX)
4375 n_output_values += 1;
4376 }
4377
4378 const char **labels;
4379 BFT_MALLOC(labels, n_output_values, const char *);
4380 for (int i = 0; i < n_output_states; i++)
4381 labels[i] = _state_names[i];
4382
4383 n_output_values = n_output_states;
4384 if (solid->post_flag & CS_SOLIDIFICATION_POST_SOLIDIFICATION_RATE)
4385 labels[n_output_values++] = "SolidRate";
4386
4387 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY)
4388 if (solid->post_flag & CS_SOLIDIFICATION_POST_SEGREGATION_INDEX)
4389 labels[n_output_values++] = "SegrIndex";
4390
4391 /* Use the physical time rather than the number of iterations */
4392
4393 if (n_output_values > 0)
4394 solid->plot_state = cs_time_plot_init_probe("solidification",
4395 "",
4396 CS_TIME_PLOT_DAT,
4397 false,
4398 180, /* flush time */
4399 -1,
4400 n_output_values,
4401 NULL,
4402 NULL,
4403 labels);
4404
4405 BFT_FREE(labels);
4406
4407 } /* rank 0 */
4408
4409 }
4410
4411 /*----------------------------------------------------------------------------*/
4412 /*!
4413 * \brief Finalize the setup stage for equations related to the solidification
4414 * module
4415 *
4416 * \param[in] connect pointer to a cs_cdo_connect_t structure
4417 * \param[in] quant pointer to a cs_cdo_quantities_t structure
4418 */
4419 /*----------------------------------------------------------------------------*/
4420
4421 void
cs_solidification_finalize_setup(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant)4422 cs_solidification_finalize_setup(const cs_cdo_connect_t *connect,
4423 const cs_cdo_quantities_t *quant)
4424 {
4425 cs_solidification_t *solid = cs_solidification_structure;
4426
4427 /* Sanity checks */
4428 if (solid == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_module));
4429
4430 const cs_lnum_t n_cells = quant->n_cells;
4431 const size_t size_c = n_cells*sizeof(cs_real_t);
4432
4433 /* Retrieve the field associated to the temperature */
4434 solid->temperature = cs_field_by_name("temperature");
4435
4436 /* Define the liquid fraction */
4437 cs_property_def_by_field(solid->g_l, solid->g_l_field);
4438
4439 /* Initially one assumes that all is liquid except for cells in a
4440 * predefined solid zone for all the computation */
4441 BFT_MALLOC(solid->cell_state, n_cells, cs_solidification_state_t);
4442
4443 cs_field_set_values(solid->g_l_field, 1.);
4444
4445 # pragma omp parallel for if (n_cells > CS_THR_MIN)
4446 for (cs_lnum_t i = 0; i < n_cells; i++) {
4447
4448 if (connect->cell_flag[i] & CS_FLAG_SOLID_CELL) {
4449 solid->g_l_field->val[i] = 0;
4450 solid->g_l_field->val_pre[i] = 0;
4451 solid->cell_state[i] = CS_SOLIDIFICATION_STATE_SOLID;
4452 }
4453 else {
4454 solid->g_l_field->val_pre[i] = 1.;
4455 solid->cell_state[i] = CS_SOLIDIFICATION_STATE_LIQUID;
4456 }
4457
4458 } /* Loop on cells */
4459
4460 if (cs_flag_test(solid->options,
4461 CS_SOLIDIFICATION_NO_VELOCITY_FIELD) == false) {
4462
4463 /* Define the forcing term acting as a reaction term in the momentum
4464 equation. This term is related to the liquid fraction */
4465 BFT_MALLOC(solid->forcing_mom_array, n_cells, cs_real_t);
4466 memset(solid->forcing_mom_array, 0, size_c);
4467
4468 cs_property_def_by_array(solid->forcing_mom,
4469 cs_flag_primal_cell,
4470 solid->forcing_mom_array,
4471 false, /* definition is owner ? */
4472 NULL); /* no index */
4473
4474 /* Add the temperature array for the Boussinesq term (thermal effect) */
4475 cs_navsto_param_t *nsp = cs_navsto_system_get_param();
4476
4477 assert(nsp->n_boussinesq_terms > 0);
4478 cs_navsto_param_boussinesq_t *bp = nsp->boussinesq_param;
4479
4480 cs_navsto_param_set_boussinesq_array(bp, solid->temperature->val);
4481
4482 }
4483
4484 /* Define the reaction coefficient and the source term for the temperature
4485 equation */
4486 if (solid->thermal_reaction_coef != NULL) {
4487
4488 BFT_MALLOC(solid->thermal_reaction_coef_array, n_cells, cs_real_t);
4489 memset(solid->thermal_reaction_coef_array, 0, size_c);
4490
4491 cs_property_def_by_array(solid->thermal_reaction_coef,
4492 cs_flag_primal_cell,
4493 solid->thermal_reaction_coef_array,
4494 false, /* definition is owner ? */
4495 NULL); /* no index */
4496
4497 BFT_MALLOC(solid->thermal_source_term_array, n_cells, cs_real_t);
4498 memset(solid->thermal_source_term_array, 0, size_c);
4499
4500 cs_equation_param_t *thm_eqp = cs_equation_param_by_name(CS_THERMAL_EQNAME);
4501 cs_equation_add_source_term_by_array(thm_eqp,
4502 NULL, /* all cells selected */
4503 cs_flag_primal_cell,
4504 solid->thermal_source_term_array,
4505 false, /* definition is owner ? */
4506 NULL); /* no index */
4507
4508 }
4509
4510 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY) {
4511 /* ==================================== */
4512
4513 cs_solidification_binary_alloy_t *alloy
4514 = (cs_solidification_binary_alloy_t *)solid->model_context;
4515
4516 /* Get a shortcut to the c_bulk field */
4517 alloy->c_bulk = cs_equation_get_field(alloy->solute_equation);
4518
4519 /* Allocate an array to store the liquid concentration */
4520 BFT_MALLOC(alloy->c_l_cells, n_cells, cs_real_t);
4521
4522 # pragma omp parallel for if (n_cells > CS_THR_MIN)
4523 for (cs_lnum_t i = 0; i < n_cells; i++)
4524 alloy->c_l_cells[i] = alloy->ref_concentration;
4525
4526 /* Add the c_l_cells array for the Boussinesq term (solutal effect) */
4527 cs_navsto_param_t *nsp = cs_navsto_system_get_param();
4528
4529 assert(nsp->n_boussinesq_terms == 2);
4530 cs_navsto_param_boussinesq_t *bp = nsp->boussinesq_param + 1;
4531
4532 cs_navsto_param_set_boussinesq_array(bp, alloy->c_l_cells);
4533
4534 /* Allocate arrays storing the intermediate states during the
4535 sub-iterations to solve the non-linearity */
4536 BFT_MALLOC(alloy->tk_bulk, n_cells, cs_real_t);
4537 BFT_MALLOC(alloy->ck_bulk, n_cells, cs_real_t);
4538
4539 if (solid->options & CS_SOLIDIFICATION_USE_EXTRAPOLATION) {
4540 BFT_MALLOC(alloy->tx_bulk, n_cells, cs_real_t);
4541 BFT_MALLOC(alloy->cx_bulk, n_cells, cs_real_t);
4542 }
4543
4544 /* Allocate eta even if SOLUTE_WITH_SOURCE_TERM is activated */
4545 const cs_real_t eta_ref_value = 1.;
4546 BFT_MALLOC(alloy->eta_coef_array, n_cells, cs_real_t);
4547 # pragma omp parallel for if (n_cells > CS_THR_MIN)
4548 for (cs_lnum_t i = 0; i < n_cells; i++)
4549 alloy->eta_coef_array[i] = eta_ref_value;
4550
4551 if (solid->options & CS_SOLIDIFICATION_WITH_SOLUTE_SOURCE_TERM) {
4552
4553 BFT_MALLOC(alloy->c_l_faces, quant->n_faces, cs_real_t);
4554 memset(alloy->c_l_faces, 0, sizeof(cs_real_t)*quant->n_faces);
4555
4556 }
4557 else { /* Estimate the reference value for the solutal diffusion property
4558 * One assumes that g_l (the liquid fraction is equal to 1) */
4559
4560 cs_property_set_reference_value(alloy->eta_coef_pty, eta_ref_value);
4561
4562 cs_property_def_by_array(alloy->eta_coef_pty,
4563 cs_flag_primal_cell,
4564 alloy->eta_coef_array,
4565 false,
4566 NULL);
4567 }
4568
4569 /* Estimate the reference value for the solutal diffusion property
4570 * One assumes that g_l (the liquid fraction) is equal to 1. */
4571 const cs_real_t pty_ref_value =
4572 solid->mass_density->ref_value*alloy->diff_coef;
4573
4574 cs_property_set_reference_value(alloy->diff_pty, pty_ref_value);
4575
4576 BFT_MALLOC(alloy->diff_pty_array, n_cells, cs_real_t);
4577
4578 # pragma omp parallel for if (n_cells > CS_THR_MIN)
4579 for (cs_lnum_t i = 0; i < n_cells; i++)
4580 alloy->diff_pty_array[i] = pty_ref_value;
4581
4582 cs_property_def_by_array(alloy->diff_pty,
4583 cs_flag_primal_cell,
4584 alloy->diff_pty_array,
4585 false,
4586 NULL);
4587
4588 if (solid->post_flag & CS_SOLIDIFICATION_ADVANCED_ANALYSIS) {
4589
4590 BFT_MALLOC(alloy->tbulk_minus_tliq, n_cells, cs_real_t);
4591 memset(alloy->tbulk_minus_tliq, 0, size_c);
4592 BFT_MALLOC(alloy->cliq_minus_cbulk, n_cells, cs_real_t);
4593 memset(alloy->cliq_minus_cbulk, 0, size_c);
4594
4595 }
4596
4597 if (solid->post_flag & CS_SOLIDIFICATION_POST_LIQUIDUS_TEMPERATURE)
4598 BFT_MALLOC(alloy->t_liquidus, n_cells, cs_real_t);
4599
4600 } /* Binary alloy model */
4601
4602 }
4603
4604 /*----------------------------------------------------------------------------*/
4605 /*!
4606 * \brief Summarize the solidification module in the log file dedicated to
4607 * the setup
4608 */
4609 /*----------------------------------------------------------------------------*/
4610
4611 void
cs_solidification_log_setup(void)4612 cs_solidification_log_setup(void)
4613 {
4614 cs_solidification_t *solid = cs_solidification_structure;
4615
4616 if (solid == NULL)
4617 return;
4618
4619 const char *module = "Solidification";
4620
4621 cs_log_printf(CS_LOG_SETUP, "\nSummary of the solidification module\n");
4622 cs_log_printf(CS_LOG_SETUP, "%s\n", cs_sep_h1);
4623
4624 cs_log_printf(CS_LOG_SETUP, " * %s | Verbosity: %d\n",
4625 module, solid->verbosity);
4626
4627 /* Display options */
4628 cs_log_printf(CS_LOG_SETUP, " * %s | Strategy:", module);
4629 switch (solid->strategy) {
4630
4631 case CS_SOLIDIFICATION_STRATEGY_LEGACY:
4632 cs_log_printf(CS_LOG_SETUP, " Legacy\n");
4633 break;
4634 case CS_SOLIDIFICATION_STRATEGY_TAYLOR:
4635 cs_log_printf(CS_LOG_SETUP, " Legacy + Taylor-based updates\n");
4636 break;
4637 case CS_SOLIDIFICATION_STRATEGY_PATH:
4638 cs_log_printf(CS_LOG_SETUP, " Rely on the solidification path\n");
4639 break;
4640
4641 default:
4642 bft_error(__FILE__, __LINE__, 0, "%s: Invalid strategy\n", __func__);
4643 }
4644
4645 switch (solid->model) {
4646 case CS_SOLIDIFICATION_MODEL_STEFAN:
4647 {
4648 cs_solidification_stefan_t *s_model =
4649 (cs_solidification_stefan_t *)solid->model_context;
4650
4651 cs_log_printf(CS_LOG_SETUP, " * %s | Model: Stefan\n", module);
4652 cs_log_printf(CS_LOG_SETUP,
4653 " * %s | Tliq/sol: %5.3e\n"
4654 " * %s | Latent heat: %5.3e\n"
4655 " * %s | Max. iter: %d; Max. delta enthalpy: %5.3e\n",
4656 module, s_model->t_change, module, solid->latent_heat,
4657 module, s_model->n_iter_max, s_model->max_delta_h);
4658 }
4659 break;
4660
4661 case CS_SOLIDIFICATION_MODEL_VOLLER_PRAKASH_87:
4662 {
4663 cs_solidification_voller_t *v_model
4664 = (cs_solidification_voller_t *)solid->model_context;
4665
4666 cs_log_printf(CS_LOG_SETUP, " * %s | Model: Voller-Prakash (1987)\n",
4667 module);
4668 if (solid->options & CS_SOLIDIFICATION_NO_VELOCITY_FIELD)
4669 cs_log_printf(CS_LOG_SETUP,
4670 " * %s | Tliq: %5.3e; Tsol: %5.3e\n"
4671 " * %s | Latent heat: %5.3e\n",
4672 module, v_model->t_liquidus, v_model->t_solidus,
4673 module, solid->latent_heat);
4674 else
4675 cs_log_printf(CS_LOG_SETUP,
4676 " * %s | Tliq: %5.3e; Tsol: %5.3e\n"
4677 " * %s | Latent heat: %5.3e\n"
4678 " * %s | Forcing coef: %5.3e s_das: %5.3e\n",
4679 module, v_model->t_liquidus, v_model->t_solidus,
4680 module, solid->latent_heat,
4681 module, solid->forcing_coef, v_model->s_das);
4682
4683 }
4684 break;
4685
4686 case CS_SOLIDIFICATION_MODEL_VOLLER_NL:
4687 {
4688 cs_solidification_voller_t *v_model
4689 = (cs_solidification_voller_t *)solid->model_context;
4690
4691 cs_iter_algo_t *ia = v_model->nl_algo;
4692 assert(ia != NULL);
4693
4694 cs_log_printf(CS_LOG_SETUP, " * %s |"
4695 " Model: Voller-Prakash (1987) with non-linearities\n",
4696 module);
4697
4698 if (solid->options & CS_SOLIDIFICATION_NO_VELOCITY_FIELD)
4699 cs_log_printf(CS_LOG_SETUP,
4700 " * %s | Tliq: %5.3e; Tsol: %5.3e\n"
4701 " * %s | Latent heat: %5.3e\n"
4702 " * %s | NL Algo: max. iter: %d; rtol: %5.3e,"
4703 " atol: %5.3e, dtol: %5.3e\n",
4704 module, v_model->t_liquidus, v_model->t_solidus,
4705 module, solid->latent_heat,
4706 module, ia->param.n_max_algo_iter, ia->param.rtol,
4707 ia->param.atol, ia->param.dtol);
4708 else
4709 cs_log_printf(CS_LOG_SETUP,
4710 " * %s | Tliq: %5.3e; Tsol: %5.3e\n"
4711 " * %s | Latent heat: %5.3e\n"
4712 " * %s | Forcing coef: %5.3e s_das: %5.3e\n"
4713 " * %s | NL Algo: max. iter: %d; rtol: %5.3e,"
4714 " atol: %5.3e, dtol: %5.3e\n",
4715 module, v_model->t_liquidus, v_model->t_solidus,
4716 module, solid->latent_heat,
4717 module, solid->forcing_coef, v_model->s_das,
4718 module, ia->param.n_max_algo_iter, ia->param.rtol,
4719 ia->param.atol, ia->param.dtol);
4720 }
4721 break;
4722
4723 case CS_SOLIDIFICATION_MODEL_BINARY_ALLOY:
4724 {
4725 cs_solidification_binary_alloy_t *alloy
4726 = (cs_solidification_binary_alloy_t *)solid->model_context;
4727
4728 cs_log_printf(CS_LOG_SETUP, " * %s | Model: Binary alloy\n",
4729 module);
4730 cs_log_printf(CS_LOG_SETUP, " * %s | Alloy: %s\n",
4731 module, cs_equation_get_name(alloy->solute_equation));
4732 cs_log_printf(CS_LOG_SETUP,
4733 " * %s | Distribution coef.: %5.3e\n"
4734 " * %s | Liquidus slope: %5.3e\n"
4735 " * %s | Phase change temp.: %5.3e\n"
4736 " * %s | Eutectic conc.: %5.3e\n"
4737 " * %s | Latent heat: %5.3e\n",
4738 module, alloy->kp,
4739 module, alloy->ml, module, alloy->t_melt,
4740 module, alloy->c_eut,
4741 module, solid->latent_heat);
4742 cs_log_printf(CS_LOG_SETUP,
4743 " * %s | Forcing coef: %5.3e; s_das: %5.3e\n",
4744 module, solid->forcing_coef, alloy->s_das);
4745
4746 cs_log_printf(CS_LOG_SETUP, " * %s | Options:", module);
4747 if (solid->options & CS_SOLIDIFICATION_BINARY_ALLOY_C_FUNC)
4748 cs_log_printf(CS_LOG_SETUP,
4749 " User-defined function for the concentration eq.");
4750 else {
4751
4752 if (cs_flag_test(solid->options,
4753 CS_SOLIDIFICATION_WITH_SOLUTE_SOURCE_TERM))
4754 cs_log_printf(CS_LOG_SETUP,
4755 " Solute concentration with an advective source term");
4756 else
4757 cs_log_printf(CS_LOG_SETUP,
4758 " Solute concentration with an advective coefficient");
4759
4760 } /* Not user-defined */
4761 cs_log_printf(CS_LOG_SETUP, "\n");
4762
4763 if (solid->options & CS_SOLIDIFICATION_BINARY_ALLOY_T_FUNC)
4764 cs_log_printf(CS_LOG_SETUP,
4765 " * %s | Options: %s\n", module,
4766 " User-defined function for the thermal equation");
4767
4768 if (solid->options & CS_SOLIDIFICATION_BINARY_ALLOY_G_FUNC)
4769 cs_log_printf(CS_LOG_SETUP,
4770 " * %s | Options: %s\n", module,
4771 " User-defined function for the liquid fraction/state");
4772
4773 if (cs_flag_test(solid->options, CS_SOLIDIFICATION_USE_EXTRAPOLATION))
4774 cs_log_printf(CS_LOG_SETUP,
4775 " * %s | Options: %s\n", module,
4776 " Update using a second-order in time extrapolation");
4777
4778 if (solid->options & CS_SOLIDIFICATION_WITH_PENALIZED_EUTECTIC) {
4779 if (solid->strategy == CS_SOLIDIFICATION_STRATEGY_PATH)
4780 cs_log_printf(CS_LOG_SETUP, " * %s | Options: %s\n", module,
4781 " Penalized eutectic temperature");
4782 else
4783 cs_log_printf(CS_LOG_SETUP, " * %s | Options: %s\n", module,
4784 " Penalized eutectic temperature (unused)");
4785 }
4786
4787 if (alloy->n_iter_max > 1)
4788 cs_log_printf(CS_LOG_SETUP, " * %s | Options: Use sub-iterations"
4789 " n_iter_max %d; tolerance: %.3e\n",
4790 module, alloy->n_iter_max, alloy->delta_tolerance);
4791
4792 } /* Binary alloy */
4793
4794 default:
4795 break;
4796
4797 } /* Switch on model type */
4798
4799 cs_log_printf(CS_LOG_SETUP, "\n");
4800 }
4801
4802 /*----------------------------------------------------------------------------*/
4803 /*!
4804 * \brief Initialize the context structure used to build the algebraic system
4805 * This is done after the setup step.
4806 *
4807 * \param[in] mesh pointer to a cs_mesh_t structure
4808 * \param[in] connect pointer to a cs_cdo_connect_t structure
4809 * \param[in] quant pointer to a cs_cdo_quantities_t structure
4810 * \param[in] time_step pointer to a cs_time_step_t structure
4811 */
4812 /*----------------------------------------------------------------------------*/
4813
4814 void
cs_solidification_initialize(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)4815 cs_solidification_initialize(const cs_mesh_t *mesh,
4816 const cs_cdo_connect_t *connect,
4817 const cs_cdo_quantities_t *quant,
4818 const cs_time_step_t *time_step)
4819 {
4820 CS_UNUSED(mesh);
4821 CS_UNUSED(connect);
4822
4823 cs_solidification_t *solid = cs_solidification_structure;
4824
4825 /* Sanity checks */
4826 if (solid == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_module));
4827
4828 /* Set the first fluid/solid cell and sanity check for the mass density in the
4829 fluid/solid zone */
4830 const cs_real_t cp0 = solid->cp->ref_value;
4831 const cs_real_t rho0 = solid->mass_density->ref_value;
4832
4833 for (int i = 0; i < cs_volume_zone_n_zones(); i++) {
4834
4835 const cs_zone_t *z = cs_volume_zone_by_id(i);
4836
4837 if (z->type & CS_VOLUME_ZONE_SOLID) /* permanent solid zone */
4838 continue;
4839
4840 else { /* fluid/solid zone according to thermodynamics conditions */
4841
4842 if (z->n_elts == 0)
4843 continue;
4844
4845 if (solid->first_cell < 0)
4846 solid->first_cell = z->elt_ids[0];
4847
4848 else {
4849
4850 cs_real_t rho = cs_property_get_cell_value(z->elt_ids[0],
4851 time_step->t_cur,
4852 solid->mass_density);
4853
4854 if (fabs(rho - rho0) > FLT_MIN)
4855 bft_error(__FILE__, __LINE__, 0,
4856 "%s: A uniform value of the mass density in the"
4857 " solidification/melting area is assumed.\n"
4858 " Please check your settings.\n"
4859 " rho0= %5.3e and rho= %5.3e in zone %s\n",
4860 __func__, rho0, rho, z->name);
4861
4862 cs_real_t cp = cs_property_get_cell_value(z->elt_ids[0],
4863 time_step->t_cur,
4864 solid->cp);
4865
4866 if (fabs(cp - cp0) > FLT_MIN)
4867 bft_error(__FILE__, __LINE__, 0,
4868 "%s: A uniform value of the Cp property in the"
4869 " solidification/melting area is assumed.\n"
4870 " Please check your settings.\n"
4871 " cp0= %5.3e and cp= %5.3e in zone %s\n",
4872 __func__, cp0, cp, z->name);
4873
4874 }
4875
4876 } /* solidification/melting zone */
4877
4878 } /* Loop on volume zones */
4879
4880 /* End of sanity checks */
4881
4882 switch (solid->model) {
4883
4884 case CS_SOLIDIFICATION_MODEL_BINARY_ALLOY:
4885 {
4886 cs_solidification_binary_alloy_t *alloy
4887 = (cs_solidification_binary_alloy_t *)solid->model_context;
4888
4889 if (solid->options & CS_SOLIDIFICATION_WITH_SOLUTE_SOURCE_TERM) {
4890
4891 if (cs_equation_get_space_scheme(alloy->solute_equation) !=
4892 CS_SPACE_SCHEME_CDOFB)
4893 bft_error(__FILE__, __LINE__, 0,
4894 " %s: Invalid space scheme for equation %s\n",
4895 __func__, cs_equation_get_name(alloy->solute_equation));
4896
4897 cs_equation_add_user_hook(alloy->solute_equation,
4898 NULL, /* hook context */
4899 _fb_solute_source_term); /* hook function */
4900
4901 /* Store the pointer to the current face temperature values */
4902 alloy->temp_faces =
4903 cs_equation_get_face_values(solid->thermal_sys->thermal_eq, false);
4904
4905 } /* CS_SOLIDIFICATION_WITH_SOLUTE_SOURCE_TERM */
4906
4907 /* One assumes that all the alloy mixture is liquid thus C_l = C_bulk */
4908 const cs_lnum_t n_cells = quant->n_cells;
4909 memcpy(alloy->c_l_cells, alloy->c_bulk->val, n_cells*sizeof(cs_real_t));
4910
4911 /* Set the previous iterate before calling update functions */
4912 memcpy(alloy->tk_bulk, solid->temperature->val, n_cells*sizeof(cs_real_t));
4913 memcpy(alloy->ck_bulk, alloy->c_bulk->val, n_cells*sizeof(cs_real_t));
4914
4915 if (alloy->c_l_faces != NULL) {
4916 cs_real_t *c_bulk_faces =
4917 cs_equation_get_face_values(alloy->solute_equation, false);
4918 memcpy(alloy->c_l_faces, c_bulk_faces, quant->n_faces*sizeof(cs_real_t));
4919 }
4920
4921 if (solid->post_flag & CS_SOLIDIFICATION_POST_ENTHALPY)
4922 _compute_enthalpy(quant,
4923 time_step->t_cur, /* t_eval */
4924 solid->temperature->val, /* temperature */
4925 solid->g_l_field->val, /* liquid fraction */
4926 alloy->t_eut, /* temp_ref */
4927 solid->latent_heat, /* latent heat coeff. */
4928 solid->mass_density, /* rho */
4929 solid->cp, /* cp */
4930 solid->enthalpy->val); /* computed enthalpy */
4931
4932 } /* CS_SOLIDIFICATION_MODEL_BINARY_ALLOY */
4933 break;
4934
4935 case CS_SOLIDIFICATION_MODEL_VOLLER_PRAKASH_87:
4936 case CS_SOLIDIFICATION_MODEL_VOLLER_NL:
4937 {
4938 cs_solidification_voller_t *v_model
4939 = (cs_solidification_voller_t *)solid->model_context;
4940
4941 v_model->update_gl(mesh, connect, quant, time_step);
4942
4943 if ( (solid->post_flag & CS_SOLIDIFICATION_POST_ENTHALPY) ||
4944 (solid->model == CS_SOLIDIFICATION_MODEL_VOLLER_NL) )
4945 _compute_enthalpy(quant,
4946 time_step->t_cur, /* t_eval */
4947 solid->temperature->val, /* temperature */
4948 solid->g_l_field->val, /* liquid fraction */
4949 v_model->t_solidus, /* temp_ref */
4950 solid->latent_heat, /* latent heat coeff. */
4951 solid->mass_density, /* rho */
4952 solid->cp, /* cp */
4953 solid->enthalpy->val); /* computed enthalpy */
4954
4955 }
4956 break;
4957
4958 case CS_SOLIDIFICATION_MODEL_STEFAN:
4959 {
4960 cs_solidification_stefan_t *s_model
4961 = (cs_solidification_stefan_t *)solid->model_context;
4962
4963 /* Temperature has been initialized.
4964 * Compute a first guess for the liquid fraction knowing that the liquid
4965 * fractions is a step function w.r.t. the temperature. Thus, one assumes
4966 * that the liquid fraction is either 0 or 1 after this step. If the
4967 * temperature is equal to T_ch, one assumes that g_l = 1
4968 *
4969 * Initialize source term and reaction term
4970 */
4971
4972 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
4973 for (cs_lnum_t c = 0; c < quant->n_cells; c++) {
4974
4975 if (solid->temperature->val[c] < s_model->t_change) {
4976 /* In the solid part */
4977 solid->g_l_field->val[c] = 0;
4978 solid->cell_state[c] = CS_SOLIDIFICATION_STATE_SOLID;
4979 }
4980 else { /* In the liquid part */
4981 solid->g_l_field->val[c] = 1;
4982 solid->cell_state[c] = CS_SOLIDIFICATION_STATE_LIQUID;
4983 }
4984
4985 /* No source term and reaction term at the begining. One assumes that
4986 there is no phase change at the first sub-iteration. */
4987 solid->thermal_reaction_coef_array[c] = 0;
4988 solid->thermal_source_term_array[c] = 0;
4989
4990 } /* Loop on cells */
4991
4992 /* Now compute the enthalpy knowing the temperature and the liquid
4993 fraction */
4994 _compute_enthalpy(quant,
4995 time_step->t_cur, /* t_eval */
4996 solid->temperature->val, /* temperature */
4997 solid->g_l_field->val, /* liquid fraction */
4998 s_model->t_change, /* temp_ref */
4999 solid->latent_heat, /* latent heat coeff. */
5000 solid->mass_density, /* rho */
5001 solid->cp, /* cp */
5002 solid->enthalpy->val); /* computed enthalpy */
5003
5004 } /* Stefan model */
5005 break;
5006
5007 default:
5008 break; /* Nothing to do */
5009
5010 } /* Switch on model */
5011 }
5012
5013 /*----------------------------------------------------------------------------*/
5014 /*!
5015 * \brief Solve equations related to the solidification module
5016 *
5017 * \param[in] mesh pointer to a cs_mesh_t structure
5018 * \param[in] connect pointer to a cs_cdo_connect_t structure
5019 * \param[in] quant pointer to a cs_cdo_quantities_t structure
5020 * \param[in] time_step pointer to a cs_time_step_t structure
5021 */
5022 /*----------------------------------------------------------------------------*/
5023
5024 void
cs_solidification_compute(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step)5025 cs_solidification_compute(const cs_mesh_t *mesh,
5026 const cs_cdo_connect_t *connect,
5027 const cs_cdo_quantities_t *quant,
5028 const cs_time_step_t *time_step)
5029 {
5030 cs_solidification_t *solid = cs_solidification_structure;
5031
5032 /* Sanity checks */
5033 if (solid == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_module));
5034
5035 switch (solid->model) {
5036
5037 case CS_SOLIDIFICATION_MODEL_BINARY_ALLOY:
5038 _default_binary_coupling(mesh, connect, quant, time_step);
5039 break;
5040
5041 case CS_SOLIDIFICATION_MODEL_VOLLER_PRAKASH_87:
5042 _voller_prakash_87(mesh, connect, quant, time_step);
5043 break;
5044
5045 case CS_SOLIDIFICATION_MODEL_VOLLER_NL:
5046 _voller_non_linearities(mesh, connect, quant, time_step);
5047 break;
5048
5049 case CS_SOLIDIFICATION_MODEL_STEFAN:
5050 _stefan_thermal_non_linearities(mesh, connect, quant, time_step);
5051 break;
5052
5053 default:
5054 break; /* Nothing else to do */
5055
5056 } /* Switch on model */
5057
5058 /* Solve the Navier-Stokes system */
5059
5060 if ((solid->options & CS_SOLIDIFICATION_NO_VELOCITY_FIELD) == 0)
5061 /* The Navier-Stokes is not solved when the frozen field is set */
5062 cs_navsto_system_compute(mesh, connect, quant, time_step);
5063
5064 /* Perform the monitoring */
5065 if (solid->verbosity > 0)
5066 _do_monitoring(quant);
5067 }
5068
5069 /*----------------------------------------------------------------------------*/
5070 /*!
5071 * \brief Predefined extra-operations for the solidification module
5072 *
5073 * \param[in] connect pointer to a cs_cdo_connect_t structure
5074 * \param[in] quant pointer to a cs_cdo_quantities_t structure
5075 * \param[in] ts pointer to a cs_time_step_t structure
5076 */
5077 /*----------------------------------------------------------------------------*/
5078
5079 void
cs_solidification_extra_op(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts)5080 cs_solidification_extra_op(const cs_cdo_connect_t *connect,
5081 const cs_cdo_quantities_t *quant,
5082 const cs_time_step_t *ts)
5083 {
5084 cs_solidification_t *solid = cs_solidification_structure;
5085
5086 if (solid == NULL)
5087 return;
5088
5089 /* Estimate the number of values to output */
5090
5091 int n_output_values = CS_SOLIDIFICATION_N_STATES - 1;
5092 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY) {
5093 n_output_values += 1;
5094
5095 if (solid->post_flag & CS_SOLIDIFICATION_POST_SEGREGATION_INDEX)
5096 n_output_values += 1;
5097
5098 }
5099
5100 if (solid->post_flag & CS_SOLIDIFICATION_POST_SOLIDIFICATION_RATE)
5101 n_output_values += 1;
5102
5103 /* Compute the output values */
5104
5105 cs_real_t *output_values = NULL;
5106 BFT_MALLOC(output_values, n_output_values, cs_real_t);
5107 memset(output_values, 0, n_output_values*sizeof(cs_real_t));
5108
5109 int n_output_states = (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY) ?
5110 CS_SOLIDIFICATION_N_STATES : CS_SOLIDIFICATION_N_STATES - 1;
5111 for (int i = 0; i < n_output_states; i++)
5112 output_values[i] = solid->state_ratio[i];
5113
5114 n_output_values = n_output_states;
5115
5116 if (solid->post_flag & CS_SOLIDIFICATION_POST_SOLIDIFICATION_RATE) {
5117
5118 const cs_real_t *gl = solid->g_l_field->val;
5119
5120 cs_real_t integr = 0;
5121 for (cs_lnum_t i = 0; i < quant->n_cells; i++) {
5122 if (connect->cell_flag[i] & CS_FLAG_SOLID_CELL)
5123 continue;
5124 integr += (1 - gl[i])*quant->cell_vol[i];
5125 }
5126
5127 /* Parallel reduction */
5128
5129 cs_parall_sum(1, CS_REAL_TYPE, &integr);
5130
5131 output_values[n_output_values] = integr/quant->vol_tot;
5132 n_output_values++;
5133
5134 }
5135
5136 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY) {
5137
5138 cs_solidification_binary_alloy_t *alloy
5139 = (cs_solidification_binary_alloy_t *)solid->model_context;
5140 assert(alloy != NULL);
5141
5142 const cs_real_t *c_bulk = alloy->c_bulk->val;
5143
5144 if (solid->post_flag & CS_SOLIDIFICATION_POST_SEGREGATION_INDEX) {
5145
5146 const cs_real_t inv_cref = 1./alloy->ref_concentration;
5147
5148 cs_real_t si = 0;
5149 for (cs_lnum_t i = 0; i < quant->n_cells; i++) {
5150 if (connect->cell_flag[i] & CS_FLAG_SOLID_CELL)
5151 continue;
5152 double c = (c_bulk[i] - alloy->ref_concentration)*inv_cref;
5153 si += c*c*quant->cell_vol[i];
5154 }
5155
5156 /* Parallel reduction */
5157
5158 cs_parall_sum(1, CS_REAL_TYPE, &si);
5159
5160 output_values[n_output_values] = sqrt(si/quant->vol_tot);
5161 n_output_values++;
5162
5163 }
5164
5165 if (solid->post_flag & CS_SOLIDIFICATION_POST_LIQUIDUS_TEMPERATURE) {
5166 assert(alloy->t_liquidus != NULL);
5167
5168 /* Compute the value to be sure that it corresponds to the current
5169 state */
5170
5171 for (cs_lnum_t i = 0; i < quant->n_cells; i++) {
5172 if (connect->cell_flag[i] & CS_FLAG_SOLID_CELL)
5173 alloy->t_liquidus[i] = -999.99; /* no physical meaning */
5174 else
5175 alloy->t_liquidus[i] = _get_t_liquidus(alloy, alloy->c_bulk->val[i]);
5176 }
5177
5178 }
5179
5180 if ((solid->post_flag & CS_SOLIDIFICATION_ADVANCED_ANALYSIS) > 0) {
5181
5182 assert(alloy->t_liquidus != NULL &&
5183 alloy->cliq_minus_cbulk != NULL &&
5184 alloy->tbulk_minus_tliq != NULL);
5185
5186 const cs_real_t *c_l = alloy->c_l_cells;
5187 const cs_real_t *t_bulk = solid->temperature->val;
5188
5189 /* Compute Cbulk - Cliq */
5190
5191 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
5192
5193 if (connect->cell_flag[c_id] & CS_FLAG_SOLID_CELL)
5194 continue; /* = 0 by default */
5195
5196 const cs_real_t conc = c_bulk[c_id];
5197 const cs_real_t temp = t_bulk[c_id];
5198
5199 alloy->cliq_minus_cbulk[c_id] = c_l[c_id] - conc;
5200 alloy->tbulk_minus_tliq[c_id] = temp - alloy->t_liquidus[c_id];
5201
5202 } /* Loop on cells */
5203
5204 } /* Advanced analysis */
5205
5206 } /* Binary alloy modelling */
5207
5208 if (cs_glob_rank_id < 1 && solid->plot_state != NULL)
5209 cs_time_plot_vals_write(solid->plot_state,
5210 ts->nt_cur,
5211 ts->t_cur,
5212 n_output_values,
5213 output_values);
5214
5215 BFT_FREE(output_values);
5216 }
5217
5218 /*----------------------------------------------------------------------------*/
5219 /*!
5220 * \brief Predefined post-processing output for the solidification module.
5221 * Prototype of this function is fixed since it is a function pointer
5222 * defined in cs_post.h (\ref cs_post_time_mesh_dep_output_t)
5223 *
5224 * \param[in, out] input pointer to a optional structure (here a
5225 * cs_gwf_t structure)
5226 * \param[in] mesh_id id of the output mesh for the current call
5227 * \param[in] cat_id category id of the output mesh for this call
5228 * \param[in] ent_flag indicate global presence of cells (ent_flag[0]),
5229 * interior faces (ent_flag[1]), boundary faces
5230 * (ent_flag[2]), particles (ent_flag[3]) or probes
5231 * (ent_flag[4])
5232 * \param[in] n_cells local number of cells of post_mesh
5233 * \param[in] n_i_faces local number of interior faces of post_mesh
5234 * \param[in] n_b_faces local number of boundary faces of post_mesh
5235 * \param[in] cell_ids list of cells (0 to n-1)
5236 * \param[in] i_face_ids list of interior faces (0 to n-1)
5237 * \param[in] b_face_ids list of boundary faces (0 to n-1)
5238 * \param[in] time_step pointer to a cs_time_step_t struct.
5239 */
5240 /*----------------------------------------------------------------------------*/
5241
5242 void
cs_solidification_extra_post(void * input,int mesh_id,int cat_id,int ent_flag[5],cs_lnum_t n_cells,cs_lnum_t n_i_faces,cs_lnum_t n_b_faces,const cs_lnum_t cell_ids[],const cs_lnum_t i_face_ids[],const cs_lnum_t b_face_ids[],const cs_time_step_t * time_step)5243 cs_solidification_extra_post(void *input,
5244 int mesh_id,
5245 int cat_id,
5246 int ent_flag[5],
5247 cs_lnum_t n_cells,
5248 cs_lnum_t n_i_faces,
5249 cs_lnum_t n_b_faces,
5250 const cs_lnum_t cell_ids[],
5251 const cs_lnum_t i_face_ids[],
5252 const cs_lnum_t b_face_ids[],
5253 const cs_time_step_t *time_step)
5254 {
5255 CS_UNUSED(n_i_faces);
5256 CS_UNUSED(n_b_faces);
5257 CS_UNUSED(cell_ids);
5258 CS_UNUSED(i_face_ids);
5259 CS_UNUSED(b_face_ids);
5260
5261 if (input == NULL)
5262 return;
5263
5264 cs_solidification_t *solid = (cs_solidification_t *)input;
5265
5266 if (cat_id == CS_POST_MESH_PROBES) {
5267
5268 cs_field_t *fld = cs_field_by_name_try("liquid_fraction");
5269 assert(fld != NULL);
5270
5271 cs_post_write_probe_values(mesh_id,
5272 CS_POST_WRITER_ALL_ASSOCIATED,
5273 "liquid_fraction",
5274 fld->dim,
5275 CS_POST_TYPE_cs_real_t,
5276 CS_MESH_LOCATION_CELLS,
5277 NULL,
5278 NULL,
5279 fld->val,
5280 time_step);
5281
5282 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY) {
5283
5284 cs_solidification_binary_alloy_t *alloy
5285 = (cs_solidification_binary_alloy_t *)solid->model_context;
5286
5287 cs_post_write_probe_values(mesh_id,
5288 CS_POST_WRITER_ALL_ASSOCIATED,
5289 "C_l",
5290 1,
5291 CS_POST_TYPE_cs_real_t,
5292 CS_MESH_LOCATION_CELLS,
5293 NULL,
5294 NULL,
5295 alloy->c_l_cells,
5296 time_step);
5297
5298 if (solid->post_flag & CS_SOLIDIFICATION_POST_LIQUIDUS_TEMPERATURE) {
5299 assert(alloy->t_liquidus != NULL);
5300 cs_post_write_probe_values(mesh_id,
5301 CS_POST_WRITER_ALL_ASSOCIATED,
5302 "Tliquidus",
5303 1,
5304 CS_POST_TYPE_cs_real_t,
5305 CS_MESH_LOCATION_CELLS,
5306 NULL,
5307 NULL,
5308 alloy->t_liquidus,
5309 time_step);
5310 }
5311
5312 if (solid->post_flag & CS_SOLIDIFICATION_ADVANCED_ANALYSIS) {
5313
5314 cs_post_write_probe_values(mesh_id,
5315 CS_POST_WRITER_ALL_ASSOCIATED,
5316 "delta_cliq_minus_cbulk",
5317 1,
5318 CS_POST_TYPE_cs_real_t,
5319 CS_MESH_LOCATION_CELLS,
5320 NULL,
5321 NULL,
5322 alloy->cliq_minus_cbulk,
5323 time_step);
5324
5325 cs_post_write_probe_values(mesh_id,
5326 CS_POST_WRITER_ALL_ASSOCIATED,
5327 "delta_tbulk_minus_tliq",
5328 1,
5329 CS_POST_TYPE_cs_real_t,
5330 CS_MESH_LOCATION_CELLS,
5331 NULL,
5332 NULL,
5333 alloy->tbulk_minus_tliq,
5334 time_step);
5335
5336 if (alloy->eta_coef_array != NULL)
5337 cs_post_write_probe_values(mesh_id,
5338 CS_POST_WRITER_ALL_ASSOCIATED,
5339 "Cbulk_advection_scaling",
5340 1,
5341 CS_POST_TYPE_cs_real_t,
5342 CS_MESH_LOCATION_CELLS,
5343 NULL,
5344 NULL,
5345 alloy->eta_coef_array,
5346 time_step);
5347
5348 } /* Advanced analysis */
5349
5350 } /* Binary alloy model */
5351
5352 } /* Probes */
5353
5354 if ((cat_id == CS_POST_MESH_VOLUME) &&
5355 (ent_flag[0] == 1)) { /* ent_flag == 1 --> on cells */
5356
5357 if (solid->cell_state != NULL &&
5358 (solid->post_flag & CS_SOLIDIFICATION_POST_CELL_STATE)) {
5359
5360 cs_post_write_var(CS_POST_MESH_VOLUME,
5361 CS_POST_WRITER_DEFAULT,
5362 "cell_state",
5363 1,
5364 false, /* interlace */
5365 true, /* true = original mesh */
5366 CS_POST_TYPE_int,
5367 solid->cell_state, NULL, NULL,
5368 time_step);
5369
5370 }
5371
5372 if (solid->model == CS_SOLIDIFICATION_MODEL_BINARY_ALLOY) {
5373
5374 cs_solidification_binary_alloy_t *alloy
5375 = (cs_solidification_binary_alloy_t *)solid->model_context;
5376
5377 cs_real_t *wb = cs_equation_get_tmpbuf();
5378
5379 if (solid->post_flag & CS_SOLIDIFICATION_ADVANCED_ANALYSIS) {
5380
5381 if (alloy->cliq_minus_cbulk != NULL)
5382 cs_post_write_var(CS_POST_MESH_VOLUME,
5383 CS_POST_WRITER_DEFAULT,
5384 "delta_cliq_minus_cbulk",
5385 1,
5386 false, /* interlace */
5387 true, /* true = original mesh */
5388 CS_POST_TYPE_cs_real_t,
5389 alloy->cliq_minus_cbulk, NULL, NULL,
5390 time_step);
5391
5392 if (alloy->tbulk_minus_tliq != NULL)
5393 cs_post_write_var(CS_POST_MESH_VOLUME,
5394 CS_POST_WRITER_DEFAULT,
5395 "delta_tbulk_minus_tliq",
5396 1,
5397 false, /* interlace */
5398 true, /* true = original mesh */
5399 CS_POST_TYPE_cs_real_t,
5400 alloy->tbulk_minus_tliq, NULL, NULL,
5401 time_step);
5402
5403 if (alloy->eta_coef_array != NULL)
5404 cs_post_write_var(CS_POST_MESH_VOLUME,
5405 CS_POST_WRITER_DEFAULT,
5406 "Cbulk_advection_scaling",
5407 1,
5408 false, /* interlace */
5409 true, /* true = original mesh */
5410 CS_POST_TYPE_cs_real_t,
5411 alloy->eta_coef_array, NULL, NULL,
5412 time_step);
5413
5414 } /* Advanced analysis */
5415
5416 if (solid->post_flag & CS_SOLIDIFICATION_POST_LIQUIDUS_TEMPERATURE) {
5417
5418 if (alloy->t_liquidus != NULL)
5419 cs_post_write_var(CS_POST_MESH_VOLUME,
5420 CS_POST_WRITER_DEFAULT,
5421 "T_liquidus",
5422 1,
5423 false, /* interlace */
5424 true, /* true = original mesh */
5425 CS_POST_TYPE_cs_real_t,
5426 alloy->t_liquidus, NULL, NULL,
5427 time_step);
5428
5429 }
5430
5431 if (solid->post_flag & CS_SOLIDIFICATION_POST_CBULK_ADIM) {
5432
5433 const cs_real_t inv_cref = 1./alloy->ref_concentration;
5434 const cs_real_t *c_bulk = alloy->c_bulk->val;
5435
5436 for (cs_lnum_t i = 0; i < n_cells; i++)
5437 wb[i] = (c_bulk[i] - alloy->ref_concentration)*inv_cref;
5438
5439 cs_post_write_var(CS_POST_MESH_VOLUME,
5440 CS_POST_WRITER_DEFAULT,
5441 "C_bulk_adim",
5442 1,
5443 false, /* interlace */
5444 true, /* true = original mesh */
5445 CS_POST_TYPE_cs_real_t,
5446 wb, NULL, NULL,
5447 time_step);
5448
5449 } /* CS_SOLIDIFICATION_POST_CBULK_ADIM */
5450
5451 if (solid->post_flag & CS_SOLIDIFICATION_POST_CLIQ)
5452 cs_post_write_var(CS_POST_MESH_VOLUME,
5453 CS_POST_WRITER_DEFAULT,
5454 "C_l",
5455 1,
5456 false, /* interlace */
5457 true, /* true = original mesh */
5458 CS_POST_TYPE_cs_real_t,
5459 alloy->c_l_cells, NULL, NULL,
5460 time_step);
5461
5462 } /* Binary alloy model */
5463
5464 } /* volume_mesh + on cells */
5465 }
5466
5467 /*----------------------------------------------------------------------------*/
5468
5469 END_C_DECLS
5470