TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Coalescence_bulles_2groupes_PolyMAC_MPFA.cpp
1/****************************************************************************
2* Copyright (c) 2021, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#include <Coalescence_bulles_2groupes_PolyMAC_MPFA.h>
17#include <Pb_Multiphase.h>
18#include <Milieu_composite.h>
19#include <Op_Diff_Turbulent_PolyMAC_MPFA_Face.h>
20#include <Viscosite_turbulente_base.h>
21#include <Matrix_tools.h>
22#include <Array_tools.h>
23#include <Coalescence_bulles_2groupes_base.h>
24#include <math.h>
25#include <Champ_Elem_PolyMAC_MPFA.h>
26#include <Champ_Face_base.h>
27#include <Domaine_PolyMAC_MPFA.h>
28
29Implemente_instanciable(Coalescence_bulles_2groupes_PolyMAC_MPFA, "Coalescence_bulles_2groupes_elem_PolyMAC_MPFA", Source_base);
30
32{
33 return os;
34}
35
37{
38 Param param(que_suis_je());
39 param.ajouter("beta_k", &beta_k_);
40 param.ajouter("Dh", &dh_);
41 param.lire_avec_accolades_depuis(is);
42
43
44 const Pb_Multiphase *pbm = sub_type(Pb_Multiphase, equation().probleme()) ? &ref_cast(Pb_Multiphase, equation().probleme()) : nullptr;
45
46 if (!pbm || pbm->nb_phases() == 1) Process::exit(que_suis_je() + " : not needed for single-phase flow!");
47 for (int n = 0; n < pbm->nb_phases(); n++) //recherche de n_l, n_g : phase {liquide,gaz}_continu en priorite
48 {
49 if (pbm->nom_phase(n).debute_par("liquide") && (n_l < 0 || pbm->nom_phase(n).finit_par("continu"))) n_l = n;
50 if (( pbm->nom_phase(n).finit_par("group1"))) n_g1 = n;
51 if (( pbm->nom_phase(n).finit_par("group2"))) n_g2 = n;
52 }
53 if (n_l < 0) Process::exit(que_suis_je() + " : liquid phase not found!");
54 if (n_g1 < 0) Process::exit(que_suis_je() + " : group 1 not found!");
55 if (n_g2 < 0) Process::exit(que_suis_je() + " : group 2 not found!");
56
57 if (pbm->has_correlation("Coalescence_bulles_2groupes")) correlation_ = pbm->get_correlation("Coalescence_bulles_2groupes"); //correlation fournie par le bloc correlation
58 else Correlation_base::typer_lire_correlation(correlation_, *pbm, "Coalescence_bulles_2groupes", is); //sinon -> on la lit
59
60 return is;
61}
62
63void Coalescence_bulles_2groupes_PolyMAC_MPFA::dimensionner_blocs(matrices_t matrices, const tabs_t& semi_impl) const
64{
65 const Domaine_PolyMAC_MPFA& domaine = ref_cast(Domaine_PolyMAC_MPFA, equation().domaine_dis());
66 const int ne = domaine.nb_elem(), ne_tot = domaine.nb_elem_tot(), N = equation().inconnue().valeurs().line_size();
67
68 for (auto &&n_m : matrices)
69 if (n_m.first == "alpha" || n_m.first == "k" || n_m.first == "tau" || n_m.first == "omega" || n_m.first == "interfacial_area")
70 {
71 Matrice_Morse& mat = *n_m.second, mat2;
72 const DoubleTab& dep = equation().probleme().get_champ(n_m.first.c_str()).valeurs();
73 int nc = dep.dimension_tot(0),
74 M = dep.line_size();
75 Stencil sten(0, 2);
76 if (n_m.first == "alpha")
77 for (int e = 0; e < ne; e++)
78 for (int n = 0; n < N; n++) sten.append_line(N * e + n, N * e + n);
79 if (n_m.first == "interfacial_area" ) // N <= M
80 for (int e = 0; e < ne; e++)
81 for (int n = 0; n < N; n++) sten.append_line(N * e + n, M * e + n);
82 if (n_m.first == "k" || n_m.first == "tau" || n_m.first == "omega") // N <= M
83 for (int e = 0; e < ne; e++)
84 for (int n = 0; n < N; n++)
85 for (int k = 0; k < M; k++) sten.append_line(N * e + n, M * e + k);
86 //tableau_trier_retirer_doublons(sten);
87 Matrix_tools::allocate_morse_matrix(N * ne_tot, M * nc, sten, mat2);
88 mat.nb_colonnes() ? mat += mat2 : mat = mat2;
89 }
90}
91
92void Coalescence_bulles_2groupes_PolyMAC_MPFA::ajouter_blocs(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl) const
93{
94 const Pb_Multiphase& pbm = ref_cast(Pb_Multiphase, equation().probleme());
95 const Domaine_PolyMAC_MPFA& domaine = ref_cast(Domaine_PolyMAC_MPFA, equation().domaine_dis());
96 const DoubleVect& pe = equation().milieu().porosite_elem(), &ve = domaine.volumes();
97
98
99 const DoubleTab& inco = equation().inconnue().valeurs(),
100 &inco_p = equation().inconnue().passe(),
101 &d_b = equation().probleme().get_champ("diametre_bulles").passe(),
102 &alpha = pbm.equation_masse().inconnue().valeurs(),
103 &alpha_p = pbm.equation_masse().inconnue().passe(),
104 &press_p = ref_cast(QDM_Multiphase,pbm.equation_qdm()).pression().passe(),
105 &temp_p = pbm.equation_energie().inconnue().passe(),
106 &rho_p = equation().milieu().masse_volumique().passe(),
107 &nu_p = equation().probleme().get_champ("viscosite_cinematique").passe(),
108 *tab_k_p = equation().probleme().has_champ("k") ? &equation().probleme().get_champ("k").passe() : nullptr,
109 *tab_k = equation().probleme().has_champ("k") ? &equation().probleme().get_champ("k").valeurs() : nullptr,
110 *tau = equation().probleme().has_champ("tau") ? &equation().probleme().get_champ("tau").valeurs() : nullptr,
111 *omega = equation().probleme().has_champ("omega") ? &equation().probleme().get_champ("omega").valeurs() : nullptr ;
112
113 const Milieu_composite& milc = ref_cast(Milieu_composite, equation().milieu());
114 int N = pbm.nb_phases(), Nk = (tab_k_p) ? (*tab_k_p).line_size() : -1, Np = equation().probleme().get_champ("pression").valeurs().line_size();
115
116 // Models use epsilon but with omega and tau we induce new variations of tau/omega and k
117 std::string Type_diss = "other"; // omega, tau or other dissipation
118 if (tau) Type_diss = "tau";
119 else if (omega) Type_diss = "omega";
120
121 DoubleTrav epsilon(alpha);
122 const Op_Diff_Turbulent_PolyMAC_MPFA_Face& op_diff = ref_cast(Op_Diff_Turbulent_PolyMAC_MPFA_Face, equation().probleme().equation(0).operateur(0).l_op_base());
123 const Viscosite_turbulente_base& visc_turb = ref_cast(Viscosite_turbulente_base, op_diff.correlation());
124 visc_turb.eps(epsilon); // Epsilon is in the past
125 double limiter = visc_turb.limiteur();
126 double dh = dh_;
127
128
129 Matrice_Morse *Ma = matrices.count("alpha") ? matrices.at("alpha") : nullptr,
130 *Mk = matrices.count("k") ? matrices.at("k") : nullptr,
131 *Mtau = matrices.count("tau") ? matrices.at("tau") : nullptr,
132 *Momega = matrices.count("omega") ? matrices.at("omega") : nullptr,
133 *Mai = matrices.count("interfacial_area") ? matrices.at("interfacial_area") : nullptr;
134
135 int cR = (rho_p.dimension_tot(0) == 1), cM = (nu_p.dimension_tot(0) == 1), n, k, e, d, D = dimension;
136 DoubleTrav a_l(N), a_l_p(N), p_l(N), T_l(N), rho_l(N), nu_l(N), sigma_l(N,N), dv(N, N), d_bulles(N), eps_l(Nk), k_l(Nk), coeff_RC(N, N), coeff_WE(N, N); //arguments pour coeff
137 const Coalescence_bulles_2groupes_base& correlation_coal = ref_cast(Coalescence_bulles_2groupes_base, correlation_.valeur());
138
139 // fill velocity at elem tab
140 DoubleTab pvit_elem(0, N * D);
141 domaine.domaine().creer_tableau_elements(pvit_elem);
142 const Champ_Face_base& ch_vit = ref_cast(Champ_Face_base,ref_cast(Pb_Multiphase, equation().probleme()).equation_qdm().inconnue());
143 ch_vit.get_elem_vector_field(pvit_elem);
144
145 const double fac_sec = 1.e4 ; // numerical security
146 const double alpha_sec = 0.509 ; // cf C.Bazin thesis
147 const double alpha_min = 1.e-3 ; // to avoid numerical problems
148
149 /* elements */
150 for (e = 0; e < domaine.nb_elem(); e++)
151 {
152 // Get field values for correlations-------------------------------------------------------------------------------------------------------------------------------------------------------
153 for (n = 0; n < N; n++)
154 {
155 a_l(n) = alpha(e, n);
156 a_l_p(n) = alpha_p(e, n);
157 p_l(n) = press_p(e, n * (Np > 1));
158 T_l(n) = temp_p(e, n);
159 rho_l(n) = rho_p(!cR * e, n);
160 nu_l(n) = nu_p(!cM * e, n);
161 for (k = 0; k < N; k++)
162 if(milc.has_interface(n, k))
163 {
164 Interface_base& sat = milc.get_interface(n, k);
165 sigma_l(n,k) = sat.sigma(temp_p(e,n), press_p(e,n * (Np > 1)));
166 }
167 else if (milc.has_saturation(n, k))
168 {
169 Saturation_base& z_sat = milc.get_saturation(n, k);
170
171 DoubleTab& sig = z_sat.get_sigma_tab();
172 sigma_l(n,k) = sig(e);
173
174 }
175
176 d_bulles(n) = d_b(e,n);
177
178 }
179
180 for (n = 0; n < Nk; n++)
181 {
182 eps_l(n) =epsilon(e, n) ;
183 k_l(n) = (tab_k_p) ? (*tab_k_p)(e,n) : 0;
184 }
185
186 for (dv =0, d = 0; d < D; d++)
187 for (n = 0; n < N; n++)
188 for (k = 0 ; k<N ; k++) dv(n, k) += (pvit_elem(e, N * d + n) - ((n!=k) ? pvit_elem(e, N * d + k) : 0) ) * (pvit_elem(e, N * d + n) - ((n!=k) ? pvit_elem(e, N * d + k) : 0) ); // nv(n,n) = ||v(n)||, nv(n, k!=n) = ||v(n)-v(k)||
189 for (n = 0; n < N; n++)
190 for (k = 0 ; k<N ; k++) dv(n, k) = sqrt(dv(n, k)) ;
191
192 // Get correlations-------------------------------------------------------------------------------------------------------------------------------------------------------
193 correlation_coal.coefficient_RC(a_l, a_l_p, p_l, T_l, rho_l, nu_l, sigma_l, dh, dv, d_bulles, eps_l, k_l, n_l, n_g1, n_g2, coeff_RC); // Semi-Explicit coeff for Random Collisions : alpha is implicit
194 correlation_coal.coefficient_WE(a_l_p, p_l, T_l, rho_l, nu_l, sigma_l, dh, dv, d_bulles, eps_l, k_l, n_l, n_g1, n_g2, coeff_WE); // Semi-Explicit coeff for Wake entrainment : alpha is implicit
195
196 // Get epsilon according to model
197 double eps_valeurs = epsilon(e, n_l) ;
198 if (Type_diss == "tau") eps_valeurs = beta_k_ * ((*tab_k)(e, n_l)>1.e-8 ? (*tab_k)(e, n_l)*(*tab_k)(e, n_l)/ std::max((*tab_k)(e, n_l) * (*tau)(e, n_l), limiter * nu_p(e, n_l)) : 0 );
199 else if (Type_diss == "omega") eps_valeurs = beta_k_ * ((*tab_k)(e, n_l)*(*omega)(e, n_l)) ;
200 else eps_valeurs = epsilon(e, n_l);
201
202
203 // Recuring functions for source terms-------------------------------------------------------------------------------------------------------------------------------------------------------
204 // fac = prefactor, da = void derivative
205
206 const double fac_RC1 = (alpha(e, n_g1)>alpha_min) ? pe(e) * ve(e) * coeff_RC(n_g1, n_l) : 0.;
207 const double dfac_RC1da = (alpha(e, n_g1)>alpha_min) ? pe(e) * ve(e) * coeff_RC(n_l, n_g1) : 0. ;
208 const double fac_WE1 = (alpha(e, n_g1)>alpha_min) ? pe(e) * ve(e) * coeff_WE(n_g1, n_l) : 0. ;
209 const double dfac_WE1da = (alpha(e, n_g1)>alpha_min) ? pe(e) * ve(e) * coeff_WE(n_l, n_g1) : 0.;
210 const double fac_RC2 = (alpha(e, n_g2)>alpha_min) ? pe(e) * ve(e) * coeff_RC(n_g2, n_l) : 0.;
211 const double dfac_RC2da = (alpha(e, n_g2)>alpha_min) ? pe(e) * ve(e) * coeff_RC(n_l, n_g2) : 0. ;
212 const double fac_WE2 = (alpha(e, n_g2)>alpha_min) ? pe(e) * ve(e) * coeff_WE(n_g2, n_l) : 0. ;
213 const double fac_RC112 = (alpha(e, n_g1)>alpha_min) ? pe(e) * ve(e) * coeff_RC(n_g1, n_g1) : 0.;
214 const double fac_RC122_1 =(alpha(e, n_g2)>alpha_min) ? pe(e) * ve(e) * coeff_RC(n_g1, n_g2) : 0.;
215 const double dfac_RC122_1da = (alpha(e, n_g2)>alpha_min) ? pe(e) * ve(e) * coeff_RC(n_l, n_l) : 0. ;
216 const double fac_RC122_2 = (alpha(e, n_g2)>alpha_min) ? pe(e) * ve(e) * coeff_RC(n_g2, n_g1) : 0. ;
217 const double fac_WE112 = (alpha(e, n_g1)>alpha_min) ? pe(e) * ve(e) * coeff_WE(n_g1, n_g1) : 0.;
218 const double fac_WE122_1 = (alpha(e, n_g2)>alpha_min) ? pe(e) * ve(e) * coeff_WE(n_g1, n_g2) : 0.;
219 const double fac_WE122_2 = (alpha(e, n_g2)>alpha_min) ? pe(e) * ve(e) * coeff_WE(n_g2, n_g1) : 0.;
220
221
222 const double ai1_p = std::max(inco_p(e, n_g1), 0.) ; // security inco negative
223 const double ai1 = std::max(inco(e, n_g1), 0.) ; // security inco negative
224
225 const double ai2_p = std::max(inco_p(e, n_g2), 0.) ; // security inco negative
226 const double ai2 = std::max(inco(e, n_g2), 0.) ; // security inco negative
227 const double eps_1_over3 = std::cbrt(eps_valeurs) ;
228
229 const double alphag1_1_over3 = std::cbrt(std::min(alpha(e, n_g1), alpha_sec)) ;
230 const double alphag1_p_1_over3 = std::cbrt(std::min(alpha_p(e, n_g1), alpha_sec)) ;
231
232 const double alphag2_1_over3 = std::cbrt(alpha(e, n_g2)) ;
233 const double alphag2p_1_over3 = std::cbrt(alpha_p(e, n_g2)) ;
234
235 const double ai2_1_over3 = std::cbrt(ai2) ;
236 const double ai2p_1_over3 = std::cbrt(ai2_p) ;
237
238 const double ai1_5_over3 = std::cbrt(ai1) * std::cbrt(ai1) * std::cbrt(ai1) * std::cbrt(ai1) * std::cbrt(ai1) ;
239 const double ai1p_5_over3 = std::cbrt(ai1_p) * std::cbrt(ai1_p) * std::cbrt(ai1_p) * std::cbrt(ai1_p) * std::cbrt(ai1_p) ;
240 const double ai2_5_over3 = std::cbrt(ai2) * std::cbrt(ai2) * std::cbrt(ai2) * std::cbrt(ai2) * std::cbrt(ai2) ;
241
242 // Fill the matrix--------------------------------------------------------------------------------------------------------------------------------------------------
243
244 // RC (1)--------------------------------------------------------------------------------------------------------------------------
245
246 secmem(e , n_g1) += (fac_RC1 > 0. ) ? fac_RC1 * std::min(alpha(e, n_g1), alpha_sec) * ai1_5_over3 * eps_1_over3 : 0. ; // (alpha1, ai1, epsilon) implicit dependance
247
248 // RC (2)---------------------------------------------------------------------------------------------------------------------------
249
250 secmem(e , n_g2) += (fac_RC2 > 0. ) ? fac_RC2 * (alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3) /std::min(ai2_1_over3,fac_sec) * eps_1_over3 : 0. ; // (alpha2, ai2, epsilon) implicit dependance
251
252 // RC (112)----------------------------------------------------------------------------------------------------------------------------
253
254 secmem(e , n_g2) += (fac_RC112 > 0. ) ? fac_RC112 * eps_1_over3 * std::min(alpha_p(e, n_g1), alpha_sec) * ai1p_5_over3: 0. ; // (epsilon) implicit dependance
255
256 // RC (122) for g1----------------------------------------------------------------------------------------------------------------------------
257
258 secmem(e , n_g1) += (fac_RC122_1 > 0. ) ? fac_RC122_1 * eps_1_over3 * (alphag1_1_over3 * alphag1_1_over3) * (alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3) * ai1 * (ai2p_1_over3 * ai2p_1_over3) : 0.; // (alpha1, ai1, epsilon) implicit dependance
259
260 // RC (122) for g2----------------------------------------------------------------------------------------------------------------------------
261
262 secmem(e , n_g2) += (fac_RC122_2 > 0. ) ? fac_RC122_2 * eps_1_over3 * (alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3) * (alphag2_1_over3) * (ai2_5_over3) : 0. ; // (alpha2, ai2, epsilon) implicit dependance
263
264 // WE (1)----------------------------------------------------------------------------------------------------------------------------
265
266 secmem(e , n_g1) += (fac_WE1 > 0. ) ? fac_WE1 * ai1 * ai1 : 0. ; // (ai1) implicit dependance
267
268 // WE (2)----------------------------------------------------------------------------------------------------------------------------
269
270 secmem(e , n_g2) += (fac_WE2 > 0. ) ? fac_WE2 * ai2 * ai2 / std::min(alpha(e, n_g2) ,fac_sec): 0.; // (alpha2, ai2) implicit dependance
271
272 // WE (112)----------------------------------------------------------------------------------------------------------------------------
273
274 secmem(e , n_g2) += (fac_WE112 > 0. ) ? fac_WE112 * ai1_p * ai1_p: 0. ;// no implicit dependance
275
276 // WE (122) for g1----------------------------------------------------------------------------------------------------------------------------
277
278 secmem(e , n_g1) += (fac_WE122_1 > 0. ) ? fac_WE122_1 * ai1 * ai2_p: 0. ; // (ai1) implicit dependance
279
280 // WE (122) for g2----------------------------------------------------------------------------------------------------------------------------
281
282 secmem(e , n_g2) += (fac_WE122_2 > 0. ) ? fac_WE122_2 * ai2 * std::min(alpha_p(e, n_g1), alpha_sec) / std::min(alpha(e, n_g2),fac_sec): 0.; // (alpha2, ai2) implicit dependance
283
284
285 if (Ma)//-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
286 {
287 // RC (1)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
288
289 (*Ma)(N * e + n_g1 , N * e + n_g1) -= (fac_RC1 > 0. ) ? fac_RC1 * ai1_5_over3 * eps_1_over3 + dfac_RC1da * std::min(alpha(e, n_g1), alpha_sec) * ai1_5_over3 * eps_1_over3: 0. ;
290
291 // RC (2)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
292
293 (*Ma)(N * e + n_g2 , N * e + n_g2) -= (fac_RC2 > 0. ) ? fac_RC2 * 7./3. * (alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3) / std::min(ai2_1_over3,fac_sec) * eps_1_over3 + dfac_RC2da * (alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3) / std::min(ai2_1_over3,fac_sec) * eps_1_over3 : 0. ;
294
295 // RC (122) for g1---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
296
297 (*Ma)(N * e + n_g1 , N * e + n_g1) -= (fac_RC122_1 > 0. ) ? fac_RC122_1 * eps_1_over3 * 2./3. / std::min(alphag1_1_over3,fac_sec) * (ai2p_1_over3 * ai2p_1_over3 * ai2p_1_over3 * ai2p_1_over3) * ai1 * (ai2p_1_over3 * ai2p_1_over3) + dfac_RC122_1da * eps_1_over3 * (alphag1_1_over3 * alphag1_1_over3) * (alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3) * ai1 * (ai2p_1_over3 * ai2p_1_over3) : 0. ;
298
299 // RC (122) for g2---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
300
301 (*Ma)(N * e + n_g2 , N * e + n_g2) -= (fac_RC122_2 > 0. ) ? fac_RC122_2 * eps_1_over3 / std::min(alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 ,fac_sec) * 1./3./ std::min( alphag2_1_over3 * alphag2_1_over3,fac_sec) * ai2_5_over3 : 0. ;//alpha
302
303 // WE (1)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
304
305 (*Ma)(N * e + n_g1 , N * e + n_g1) -=(dfac_WE1da > 0. ) ? dfac_WE1da * ai1 * ai1 : 0. ;
306
307 // WE (2)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
308
309 (*Ma)(N * e + n_g2 , N * e + n_g2) -=(fac_WE2 > 0. ) ? -fac_WE2 * (ai2 * ai2) / std::min(alpha(e, n_g2) * alpha(e, n_g2), fac_sec) : 0. ;
310
311
312 //WE (122) for g2---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
313
314 (*Ma)(N * e + n_g2 , N * e + n_g2) -=(fac_WE122_2 > 0. ) ? - fac_WE122_2 * ai2 * alpha_p(e, n_g1) / std::min(alpha(e, n_g2) * alpha(e, n_g2),fac_sec) : 0. ;
315
316 }
317 if (Mai)//-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
318 {
319 // RC (1)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
320
321 (*Mai)(N * e + n_g1 , N * e + n_g1) -=(fac_RC1 > 0. ) ? fac_RC1 * std::min(alpha(e, n_g1), alpha_sec) * 5./3. * std::cbrt(ai1) * std::cbrt(ai1) * eps_1_over3 : 0. ;
322
323 // RC (2)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
324
325 (*Mai)(N * e + n_g2 , N * e + n_g2) -=(fac_RC2 > 0. ) ? - fac_RC2 * (alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3) * 1./3. / std::min(ai2_1_over3 * ai2_1_over3 * ai2_1_over3 * ai2_1_over3,fac_sec) * eps_1_over3 : 0. ;
326
327 // RC (122) for g1---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
328
329 (*Mai)(N * e + n_g1 , N * e + n_g1) -=(fac_RC122_1 > 0. ) ? fac_RC122_1 * eps_1_over3 * (alphag1_1_over3 * alphag1_1_over3) * (alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3) * (ai2p_1_over3 * ai2p_1_over3) : 0. ;
330
331 // RC (122) for g2---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
332
333 (*Mai)(N * e + n_g2 , N * e + n_g2) -=(fac_RC122_2 > 0. ) ? fac_RC122_2 * eps_1_over3 * (alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3) * alphag2_1_over3 * 5./3. * (ai2_1_over3 * ai2_1_over3) : 0. ;
334
335 // WE (1)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
336
337 (*Mai)(N * e + n_g1 , N * e + n_g1) -=(fac_WE1 > 0. ) ? fac_WE1 * 2. * ai1 : 0. ;
338 // WE (2)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
339
340 (*Mai)(N * e + n_g2 , N * e + n_g2) -=(fac_WE2 > 0. ) ? fac_WE2 * 2. * ai2 / std::max(alpha(e, n_g2),1./fac_sec) : 0. ;
341
342 // WE (122) for g1---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
343
344 (*Mai)(N * e + n_g1 , N * e + n_g1) -=(fac_WE122_1 > 0. ) ? fac_WE122_1 * ai2_p : 0. ;
345
346 // WE (122) for g2---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
347
348 (*Mai)(N * e + n_g2 , N * e + n_g2) -=(fac_WE122_2 > 0. ) ? fac_WE122_2 * 2. * ai2 * alpha_p(e, n_g1) / std::max(alpha(e, n_g2),1./fac_sec) : 0. ;
349 }
350 if (Type_diss == "tau")//-------------------------------------------------------------------------------------------------------------------------------------------------------
351 {
352 if ((*tab_k)(e, n_l) * (*tau)(e, n_l) > limiter * nu_p(e, n_l)) // derivee en k ; depend de l'activation ou non du limiteur
353 {
354 if (Mk)//--------------------------------------------------------------------------------------------------------------------------------------------------------------------
355 {
356 const double deps = 1./3. * std::cbrt(beta_k_) / std::min(std::cbrt((*tab_k)(e, n_l)) * std::cbrt((*tab_k)(e, n_l)),fac_sec) / std::min(std::cbrt((*tau)(e, n_l)),fac_sec) ;
357
358 // RC (1)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
359
360 (*Mk)(N * e + n_g1, Nk * e + n_l) -=(fac_RC1 > 0. ) ? fac_RC1 * std::min(alpha(e, n_g1), alpha_sec) * ai1_5_over3 * 1./3. * deps : 0. ;
361
362 // RC (2)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
363
364 (*Mk)(N * e + n_g2, Nk * e + n_l) -=(fac_RC2 > 0. ) ? fac_RC2 * (alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3) / std::min(ai2_1_over3,fac_sec) * deps : 0.;
365
366 // RC (112)--------------------------------------------------------------------------------------------------------------------------------------------------------------
367
368 (*Mk)(N * e + n_g2, Nk * e + n_l) -=(fac_RC112 > 0. ) ? fac_RC112 * std::min(alpha_p(e, n_g1), alpha_sec) * ai1p_5_over3 * deps: 0. ;
369
370 // RC (122) for g1----------------------------------------------------------------------------------------------------------------------------------------------------------------
371
372 (*Mk)(N * e + n_g1, Nk * e + n_l) -=(fac_RC122_1 > 0. ) ? fac_RC122_1 * (alphag1_1_over3 * alphag1_1_over3) * (alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3) * ai1 * (ai2p_1_over3 * ai2p_1_over3) * deps : 0.;
373
374 // RC (122) for g2-------------------------------------------------------------------------------------------------------------------------------------------------------------
375
376 (*Mk)(N * e + n_g2, Nk * e + n_l) -=(fac_RC122_2 > 0. ) ? fac_RC122_2 * (alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3) * alphag2_1_over3 * (ai2_5_over3) * deps : 0.;
377
378 }
379 if (Mtau)//--------------------------------------------------------------------------------------------------------------------------------------------------------------------
380 {
381 const double deps = -1./3. * std::cbrt(beta_k_) * std::cbrt((*tab_k)(e, n_l)) / std::min(std::cbrt((*tau)(e, n_l)) * std::cbrt((*tau)(e, n_l)) * std::cbrt((*tau)(e, n_l)) * std::cbrt((*tau)(e, n_l)),fac_sec) ;
382
383 // RC (1)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
384
385 (*Mtau)(N * e + n_g1, Nk * e + n_l)-=(fac_RC1 > 0. ) ? fac_RC1 * std::min(alpha(e, n_g1), alpha_sec) * ai1_5_over3 * deps : 0.;
386
387 // RC (2)-------------------------------------------------------------------------------------------------------------------------------------------------------------------------
388
389 (*Mtau)(N * e + n_g2, Nk * e + n_l)-=(fac_RC2 > 0. ) ? fac_RC2 * (alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3) / std::min(ai2_1_over3,fac_sec) * deps : 0.;
390
391 // RC (112)----------------------------------------------------------------------------------------------------------------------------------------------------------------------
392
393 (*Mtau)(N * e + n_g2, Nk * e + n_l)-=(fac_RC112 > 0. ) ? fac_RC112 * std::min(alpha_p(e, n_g1), alpha_sec) * ai1p_5_over3 * deps : 0.;
394
395 // RC (122) for g1------------------------------------------------------------------------------------------------------------------------------------------------------------------
396
397 (*Mtau)(N * e + n_g1, Nk * e + n_l) -=(fac_RC122_1 > 0. ) ? fac_RC122_1 * (alphag1_1_over3 * alphag1_1_over3) * (alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3) * ai1 * (ai2p_1_over3 * ai2p_1_over3) * deps : 0.;
398
399 // RC (122) for g2------------------------------------------------------------------------------------------------------------------------------------------------------------------
400
401 (*Mtau)(N * e + n_g2, Nk * e + n_l) -=(fac_RC122_2 > 0. ) ? fac_RC122_2 * (alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3) * alphag2_1_over3 * (ai2_5_over3) * deps : 0.;
402
403 }
404 }
405 else if (Mk)//---------------------------------------------------------------------------------------------------------------------------------------------------------------------
406 {
407 const double deps = 1./3. * std::cbrt(beta_k_) / std::min(std::cbrt((*tab_k)(e, n_l)),fac_sec) / std::min(std::cbrt(limiter * nu_p(e, n_l)),fac_sec) ;
408
409 // RC (1)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
410
411 (*Mk)(N * e + n_g1, Nk * e + n_l) -=(fac_RC1 > 0. ) ? fac_RC1 * std::min(alpha(e, n_g1), alpha_sec) * ai1_5_over3 * deps : 0.;
412
413 // RC (2)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
414
415 (*Mk)(N * e + n_g2, Nk * e + n_l) -= (fac_RC2 > 0. ) ? fac_RC2 * (alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3) / std::min(ai2_1_over3,fac_sec) * deps : 0.;
416
417 // RC (112)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
418
419 (*Mk)(N * e + n_g2, Nk * e + n_l) -=(fac_RC112 > 0. ) ? fac_RC112 * std::min(alpha_p(e, n_g1), alpha_sec) * ai1p_5_over3 * deps : 0.;
420
421 // RC (122) for g1----------------------------------------------------------------------------------------------------------------------------------------------------------------
422
423 (*Mk)(N * e + n_g1, Nk * e + n_l) -=(fac_RC122_1 > 0. ) ? fac_RC122_1 * (alphag1_1_over3 * alphag1_1_over3) * (alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3) * ai1 * (ai2p_1_over3 * ai2p_1_over3) * deps : 0.;
424
425 //RC (122) for g2--------------------------------------------------------------------------------------------------------------------------------------------------------------------
426
427 (*Mk)(N * e + n_g2, Nk * e + n_l) -=(fac_RC122_2 > 0. ) ? fac_RC122_2 * (alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3) * alphag2_1_over3 * (ai2_5_over3) * deps : 0.;
428
429 }
430 }
431 if (Type_diss == "omega")//-------------------------------------------------------------------------------------------------------------------------------------------------------
432 {
433 if (Momega)//-----------------------------------------------------------------------------------------------------------------------------------------------------------------
434 {
435 const double deps = 1./3. * std::cbrt(beta_k_) * std::cbrt((*tab_k)(e, n_l)) / std::min(std::cbrt((*omega)(e, n_l)) * std::cbrt((*omega)(e, n_l)),fac_sec) ;
436
437 // RC (1)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
438
439 (*Momega)(N * e + n_g1 , Nk * e + n_l) -=(fac_RC1 > 0. ) ? fac_RC1 * std::min(alpha(e, n_g1), alpha_sec) * ai1_5_over3 * deps : 0.;
440
441 // RC (2)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
442
443 (*Momega)(N * e + n_g2 , Nk * e + n_l) -=(fac_RC2 > 0. ) ? fac_RC2 * (alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3) / std::min(ai2_1_over3,fac_sec) * deps : 0.;
444
445 // RC (112)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
446
447 (*Momega)(N * e + n_g2 , Nk * e + n_l) -=(fac_RC112 > 0. ) ? fac_RC112 * std::min(alpha_p(e, n_g1), alpha_sec) * ai1p_5_over3 * deps : 0.;
448
449 // RC (122) for g1----------------------------------------------------------------------------------------------------------------------------------------------------------------------
450
451 (*Momega)(N * e + n_g1, Nk * e + n_l) -=(fac_RC122_1 > 0. ) ? fac_RC122_1 * (alphag1_1_over3 * alphag1_1_over3) * (alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3) * ai1 * (ai2p_1_over3 * ai2p_1_over3) * deps : 0.;
452
453 // RC (122) for g2----------------------------------------------------------------------------------------------------------------------------------------------------------------
454
455 (*Momega)(N * e + n_g2, Nk * e + n_l) -=(fac_RC1 > 0. ) ? fac_RC122_2 * (alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3) * alphag2_1_over3 * (ai2_5_over3) * deps : 0.;
456
457 }
458 if (Mk)//---------------------------------------------------------------------------------------------------------------------------------------------------------------------
459 {
460 const double deps = 1./3. * std::cbrt(beta_k_) / std::min(std::cbrt((*tab_k)(e, n_l)) * std::cbrt((*tab_k)(e, n_l)) ,fac_sec) * std::cbrt((*omega)(e, n_l)) ;
461
462 // RC (1)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
463
464 (*Mk)(N * e + n_g1 , Nk * e + n_l) -=(fac_RC1 > 0. ) ? fac_RC1 * std::min(alpha(e, n_g1), alpha_sec) * ai1_5_over3 * deps : 0.;
465
466 // RC (2)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
467
468 (*Mk)(N * e + n_g2 , Nk * e + n_l) -=(fac_RC2 > 0. ) ? fac_RC2 * (alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3 * alphag2_1_over3) / std::min(ai2_1_over3,fac_sec) * deps : 0.;
469
470 // RC (112)---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
471
472 (*Mk)(N * e + n_g2 , Nk * e + n_l) -=(fac_RC112 > 0. ) ? fac_RC112 * std::min(alpha_p(e, n_g1), alpha_sec) * ai1p_5_over3 * deps : 0.;
473
474 // RC (122) for g1---------------------------------------------------------------------------------------------------------------------------------------------------------------------
475
476 (*Mk)(N * e + n_g1, Nk * e + n_l) -=(fac_RC122_1 > 0. ) ? fac_RC122_1 * (alphag1_1_over3 * alphag1_1_over3) * (alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3 * alphag2p_1_over3) * ai1 * (ai2p_1_over3 * ai2p_1_over3) * deps : 0.;
477
478 // RC (122) for g2--------------------------------------------------------------------------------------------------------------------------------------------------------------------
479
480 (*Mk)(N * e + n_g2, Nk * e + n_l) -=(fac_RC122_2 > 0. ) ? fac_RC122_2 * (alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3 * alphag1_p_1_over3) * alphag2_1_over3 * (ai2_5_over3) * deps : 0.;
481
482 }
483 }
484 }
485}
virtual DoubleTab & get_elem_vector_field(DoubleTab &, bool passe=false) const
DoubleTab & passe(int i=1) override
Renvoie les valeurs du champs a l'instant t-i.
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
virtual DoubleTab & valeurs()=0
virtual DoubleTab & passe(int i=1)
Definition Champ_Proto.h:50
classe Coalescence_bulles_2groupes_PolyMAC_MPFA
void ajouter_blocs(matrices_t matrices, DoubleTab &secmem, const tabs_t &semi_impl={}) const override
void dimensionner_blocs(matrices_t matrices, const tabs_t &semi_impl={}) const override
virtual void coefficient_WE(const DoubleTab &alpha, const DoubleTab &p, const DoubleTab &T, const DoubleTab &rho, const DoubleTab &nu, const DoubleTab &sigma, const double Dh, const DoubleTab &ndv, const DoubleTab &d_bulles, const DoubleTab &eps, const DoubleTab &k_turb, const int n_l, const int n_g1, const int n_g2, DoubleTab &coeff) const =0
virtual void coefficient_RC(const DoubleTab &alpha, const DoubleTab &alpha_p, const DoubleTab &p, const DoubleTab &T, const DoubleTab &rho, const DoubleTab &nu, const DoubleTab &sigma, double Dh, const DoubleTab &ndv, const DoubleTab &d_bulles, const DoubleTab &eps, const DoubleTab &k_turb, const int n_l, const int n_g1, const int n_g2, DoubleTab &coeff) const =0
static void typer_lire_correlation(OWN_PTR(Correlation_base)&, const Probleme_base &, const Nom &, Entree &)
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
virtual const Milieu_base & milieu() const =0
virtual const Champ_Inc_base & inconnue() const =0
Probleme_base & probleme()
Renvoie le probleme associe a l'equation.
DoubleTab & get_sigma_tab()
double sigma(const double T, const double P) const
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
int nb_colonnes() const override
Return local number of columns (=size on the current proc).
static void allocate_morse_matrix(const int nb_lines, const int nb_columns, const Stencil &stencil, Matrice_Morse &matrix, const bool &attach_stencil_to_matrix=false)
DoubleVect & porosite_elem()
Definition Milieu_base.h:58
virtual const Champ_base & masse_volumique() const
Renvoie la masse volumique du milieu.
Classe Milieu_composite Cette classe represente un fluide reel ainsi que.
bool has_interface(int k, int l) const
bool has_saturation(int k, int l) const
Interface_base & get_interface(int k, int l) const
Saturation_base & get_saturation(int k, int l) const
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
virtual int finit_par(const char *const n) const
Definition Nom.cpp:324
virtual int debute_par(const char *const n) const
Definition Nom.cpp:319
static int dimension
Definition Objet_U.h:99
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
: class Op_Diff_Turbulent_PolyMAC_MPFA_Face
const Correlation_base & correlation() const
classe Pb_Multiphase Cette classe represente un probleme de thermohydraulique multiphase de type "3*N...
virtual Equation_base & equation_qdm()
virtual Equation_base & equation_energie()
const Nom & nom_phase(int i) const
int nb_phases() const
virtual Equation_base & equation_masse()
bool has_champ(const Motcle &nom, OBS_PTR(Champ_base) &ref_champ) const override
const Champ_base & get_champ(const Motcle &nom) const override
int has_correlation(std::string nom_correlation) const
const Correlation_base & get_correlation(std::string nom_correlation) const
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
classe QDM_Multiphase Cette classe porte les termes de l'equation de la dynamique
Classe de base des flux de sortie.
Definition Sortie.h:52
classe Source_base Un objet Source_base est un terme apparaissant au second membre d'une
Definition Source_base.h:42
_SIZE_ dimension_tot(int) const override
Definition TRUSTTab.tpp:160
void append_line(_TYPE_)
Definition TRUSTTab.tpp:213
int line_size() const
Definition TRUSTVect.tpp:67
classe Viscosite_turbulente_base correlations de viscosite turbulente decrivant le tenseur de Reynold...
virtual void eps(DoubleTab &eps) const =0