1 /*============================================================================
2 * Physical properties management for groundwater flow module.
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 /*----------------------------------------------------------------------------*/
31
32 /*----------------------------------------------------------------------------
33 * Standard C library headers
34 *----------------------------------------------------------------------------*/
35
36 #include <assert.h>
37 #include <stdio.h>
38 #include <stdarg.h>
39 #include <stdlib.h>
40 #include <string.h>
41 #include <math.h>
42
43 /*----------------------------------------------------------------------------
44 * Local headers
45 *----------------------------------------------------------------------------*/
46
47 #include "bft_mem.h"
48 #include "bft_error.h"
49 #include "bft_printf.h"
50
51 #include "cs_gwf_parameters.h"
52 #include "cs_math.h"
53 #include "cs_field.h"
54 #include "cs_field_pointer.h"
55 #include "cs_mesh.h"
56 #include "cs_mesh_quantities.h"
57
58 /*----------------------------------------------------------------------------
59 * Header for the current file
60 *----------------------------------------------------------------------------*/
61
62 #include "cs_gwf_physical_properties.h"
63
64 /*----------------------------------------------------------------------------*/
65
66 BEGIN_C_DECLS
67
68 /*=============================================================================
69 * Public function definitions
70 *============================================================================*/
71
72 /*----------------------------------------------------------------------------*/
73 /*!
74 * \brief Update delay of all transported species (user scalars)
75 *
76 * Species transport is delayed by retention in solid phase.
77 * This delay is computed as follows:
78 * \f[
79 * R = 1 + \rho \dfrac{K_d}{\theta} ;
80 * \f]
81 * where \f$ R \f$ is the delay factor, \f$ \rho \f$ the soil (bulk) density,
82 * \f$ K_d \f$ the contaminant distribution coefficient and \f$ \theta \f$
83 * the moisture content.
84 *
85 * Please refer to the
86 * <a href="../../theory.pdf#gwf_sp_transp"><b>dedicated section</b></a>
87 * of the theory guide for more informations.
88 *
89 */
90 /*---------------------------------------------------------------------------*/
91
cs_gwf_delay_update(void)92 void cs_gwf_delay_update(void)
93 {
94 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
95
96 cs_field_t *sca, *delay, *sat, *rosoil, *kd;
97 cs_gwf_soilwater_partition_t sorption_scal;
98
99 int key_part = cs_field_key_id("gwf_soilwater_partition");
100
101 /* Get soil property fields */
102 sat = cs_field_by_name("saturation");
103 rosoil = cs_field_by_name("soil_density");
104
105 /* Loop on user scalars */
106 for (int f_id = 0; f_id < cs_field_n_fields(); f_id++) {
107 sca = cs_field_by_id(f_id);
108 if ((sca->type & CS_FIELD_VARIABLE) && (sca->type & CS_FIELD_USER)) {
109 /* Sorption properties are needed to update delay */
110 cs_field_get_key_struct(sca, key_part, &sorption_scal);
111
112 /* Update of delay */
113 kd = cs_field_by_id(sorption_scal.ikd);
114 delay = cs_field_by_id(sorption_scal.idel);
115
116 for (int c_id = 0; c_id < n_cells; c_id++) {
117 delay->val[c_id] = 1. + rosoil->val[c_id]
118 * kd->val[c_id] / sat->val[c_id];
119 }
120 }
121 }
122 }
123
124 /*----------------------------------------------------------------------------*/
125 /*!
126 * \brief Add first-order decay to implicit part of source term array
127 *
128 * Corresponding EDO for decay phenomenon is:
129 * \f[
130 * \dfrac{dc}{dt} = - decay_rate c
131 * \f]
132 *
133 * \param[in] f_id field index of scalar on which decay is set
134 * \param[in,out] ts_imp source term implicit part for scalar of index f_id
135 */
136 /*---------------------------------------------------------------------------*/
137
cs_gwf_decay_rate(const int f_id,cs_real_t * ts_imp)138 void cs_gwf_decay_rate(const int f_id,
139 cs_real_t *ts_imp)
140 {
141 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
142 const cs_mesh_quantities_t *mesh_quantities = cs_glob_mesh_quantities;
143 const cs_real_t *vol = mesh_quantities->cell_vol;
144
145 cs_field_t *sca = cs_field_by_id(f_id);
146 const int dr_id = cs_field_key_id("fo_decay_rate");
147 cs_real_t decay_rate = cs_field_get_key_double(sca, dr_id);
148
149 /* First-order decay_rate in implict term */
150 for (int c_id = 0; c_id < n_cells; c_id++)
151 ts_imp[c_id] -= decay_rate * vol[c_id];
152 }
153
154 /*----------------------------------------------------------------------------*/
155 /*!
156 * \brief Update sorbed concentration for scalars with kinetic sorption.
157 *
158 * It is estimated by the following analytical expression :
159 * \f[
160 * S^{n+1} = S^n \exp(- (k^{-} + decay_rate) * \Delta t)
161 * - C^n \dfrac{k^{+}}{k^{-}}
162 * \left(\exp(- (k^{-} + decay_rate) * \Delta t) - 1 \right)
163 * \f]
164 *
165 * \param[in] f_id field index of scalar which properties are updated
166 */
167 /*---------------------------------------------------------------------------*/
168
cs_gwf_sorbed_concentration_update(const int f_id)169 void cs_gwf_sorbed_concentration_update(const int f_id)
170 {
171 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
172 cs_real_t *dt = CS_F_(dt)->val;
173
174 /* Get scalar, sorbed concentration fields */
175 cs_field_t *sca = cs_field_by_id(f_id);
176 int key_sorb = cs_field_key_id("gwf_sorbed_concentration_id");
177 cs_field_t *sorb = cs_field_by_id(cs_field_get_key_int(sca, key_sorb));
178
179 /* Get first-order decay rate */
180 const int dr_id = cs_field_key_id("fo_decay_rate");
181 cs_real_t decay_rate = cs_field_get_key_double(sca, dr_id);
182
183 /* Get soil-water partition structure */
184 cs_gwf_soilwater_partition_t sorption_scal;
185 const int key_part = cs_field_key_id("gwf_soilwater_partition");
186 cs_field_get_key_struct(sca, key_part, &sorption_scal);
187
188 /* Get k+ and k- fields */
189 cs_field_t *kp = cs_field_by_id(sorption_scal.ikp);
190 cs_field_t *km = cs_field_by_id(sorption_scal.ikm);
191
192 /* Update sorbed concentration */
193
194 /* First choice : analytical resolution */
195 if (sorption_scal.anai){
196 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++)
197 /* Case of reversible sorption or decay rate term */
198 if (km->val[c_id] + decay_rate > cs_math_epzero) {
199 cs_real_t expkdt = exp(-(km->val[c_id] + decay_rate) * dt[c_id]);
200 cs_real_t kpskm = kp->val[c_id] / (km->val[c_id] + decay_rate);
201 sorb->val[c_id] = expkdt * sorb->val[c_id]
202 - (expkdt-1.) * kpskm * sca->val[c_id];
203 }
204 else { /* Irreversible sorption and no decay rate */
205 sorb->val[c_id] = sorb->val[c_id]
206 + dt[c_id]*kp->val[c_id]*sca->val[c_id];
207 }
208 }
209 else { /* Second choice : explicit */
210 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
211 sorb->val[c_id] += dt[c_id] * (kp->val[c_id] * sca->val[c_id]
212 - (km->val[c_id] + decay_rate) * sorb->val[c_id]);
213 }
214 }
215 }
216
217 /*----------------------------------------------------------------------------*/
218 /*!
219 * \brief Update liquid concentration according to precipitation phenomenon.
220 *
221 * If liquid concentration exceeds solubility index, then it is clipped and
222 * transferred in precipitated concentration.
223 *
224 * Note that the decay rate is applied before the clipping.
225 *
226 * \param[in] f_id field index of scalar which properties are updated
227 */
228 /*---------------------------------------------------------------------------*/
229
cs_gwf_precipitation(const int f_id)230 void cs_gwf_precipitation(const int f_id)
231 {
232 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
233
234 cs_real_t *dt = CS_F_(dt)->val;
235 cs_real_t ctot_tmp = 0.;
236 cs_field_t *sca, *precip, *solub;
237 cs_gwf_soilwater_partition_t sorption_scal;
238
239 int key_part = cs_field_key_id("gwf_soilwater_partition");
240 int key_pre = cs_field_key_id("gwf_precip_concentration_id");
241
242 /* Get scalar, precipitation and solubility index fields */
243 sca = cs_field_by_id(f_id);
244 cs_field_get_key_struct(sca, key_part, &sorption_scal);
245 precip = cs_field_by_id(cs_field_get_key_int(sca, key_pre));
246 solub = cs_field_by_id(sorption_scal.imxsol);
247
248 /* Get first-order decay rate */
249 const int dr_id = cs_field_key_id("fo_decay_rate");
250 cs_real_t decay_rate = cs_field_get_key_double(sca, dr_id);
251
252 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
253 /* Apply the radioactive decay rate to precipitation concentration before
254 clipping. To be consistent with decay treatment of liquid concentration,
255 an implicit scheme is applied. */
256 precip->val[c_id] *= 1. / (1. + decay_rate * dt[c_id]);
257
258 /* Clipping of solute concentration, update of precip. concentration */
259 ctot_tmp = sca->val[c_id] + precip->val[c_id];
260 sca->val[c_id] = CS_MIN(ctot_tmp, solub->val[c_id]);
261 precip->val[c_id] = CS_MAX(0., ctot_tmp - solub->val[c_id]);
262 }
263 }
264
265 /*----------------------------------------------------------------------------*/
266 /*!
267 * \brief Take into account kinetic chemical reaction in evolution equation
268 * of total liquid concentration.
269 *
270 * \f[
271 * S^{n+1} = S^n \exp(- (k^{-} + decay_rate) * \Delta t)
272 * - C^n \dfrac{k^{+}}{k^{-}}
273 * \left(\exp(- (k^{-} + decay_rate) * \Delta t) - 1 \right)
274 * \f]
275 *
276 * \param[in] f_id field index of scalar which properties are updated
277 */
278 /*---------------------------------------------------------------------------*/
279
cs_gwf_kinetic_reaction(const int f_id,cs_real_t * ts_imp,cs_real_t * ts_exp)280 void cs_gwf_kinetic_reaction(const int f_id,
281 cs_real_t *ts_imp,
282 cs_real_t *ts_exp)
283 {
284 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
285 const cs_mesh_quantities_t *mesh_quantities = cs_glob_mesh_quantities;
286 const cs_real_t *vol = mesh_quantities->cell_vol;
287
288 cs_real_t *dt = CS_F_(dt)->val;
289
290 /* Get soil density */
291 cs_field_t *rosoil = cs_field_by_name("soil_density");
292
293 /* Get scalar and sorbed concentration fields */
294 cs_field_t *sca = cs_field_by_id(f_id);
295 int key_sorb = cs_field_key_id("gwf_sorbed_concentration_id");
296 cs_field_t *sorb = cs_field_by_id(cs_field_get_key_int(sca, key_sorb));
297
298 /* Get first-order decay rate */
299 const int dr_id = cs_field_key_id("fo_decay_rate");
300 cs_real_t decay_rate = cs_field_get_key_double(sca, dr_id);
301
302 /* Soil-water partition structure */
303 int key_part = cs_field_key_id("gwf_soilwater_partition");
304 cs_gwf_soilwater_partition_t sorption_scal;
305 cs_field_get_key_struct(sca, key_part, &sorption_scal);
306
307 /* Get k+ and k- fields */
308 cs_field_t *kp = cs_field_by_id(sorption_scal.ikp);
309 cs_field_t *km = cs_field_by_id(sorption_scal.ikm);
310
311 /* First choice : analytical resolution */
312 if (sorption_scal.anai) {
313 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++)
314 /* General case (reversible sorption or presence of decay rate) */
315 if (km->val[c_id] + decay_rate > cs_math_epzero) {
316 cs_real_t expkdt = exp(-(km->val[c_id]+decay_rate) * dt[c_id]);
317 cs_real_t kpskm = kp->val[c_id] / (km->val[c_id] + decay_rate);
318 ts_exp[c_id] += - vol[c_id]
319 *(decay_rate * rosoil->val[c_id] * sorb->val[c_id]
320 + rosoil->val[c_id]/dt[c_id] * (1-expkdt)
321 *(kpskm*sca->val[c_id] - sorb->val[c_id]));
322 ts_imp[c_id] += vol[c_id] * rosoil->val[c_id] / dt[c_id]
323 *(1-expkdt)*kpskm;
324 }
325 else { /* case of irreversible sorption without decay rate */
326 cs_real_t rokpl = rosoil->val[c_id] * kp->val[c_id];
327 ts_exp[c_id] += - vol[c_id] * rokpl * sca->val[c_id];
328 ts_imp[c_id] += + vol[c_id] * rokpl;
329 }
330 }
331 else { /* Second choice : explicit */
332 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
333 ts_exp[c_id] += vol[c_id] * rosoil->val[c_id]
334 *( km->val[c_id] * sorb->val[c_id]
335 - kp->val[c_id] * sca->val[c_id]);
336 ts_imp[c_id] += vol[c_id] * rosoil->val[c_id] * kp->val[c_id];
337 }
338 }
339 }
340
341 /*----------------------------------------------------------------------------*/
342
343 END_C_DECLS
344