TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Calcul_Production_K_VDF.cpp
1/****************************************************************************
2* Copyright (c) 2015 - 2016, 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 <Calcul_Production_K_VDF.h>
17#include <Transport_K_Eps.h>
18#include <Convection_Diffusion_Temperature.h>
19#include <Convection_Diffusion_Concentration.h>
20#include <Modele_turbulence_scal_base.h>
21#include <Probleme_base.h>
22#include <TRUSTTrav.h>
23#include <Dirichlet_entree_fluide_leaves.h>
24#include <Entree_fluide_concentration_imposee.h>
25#include <Champ_Uniforme.h>
26#include <Domaine_VDF.h>
27#include <Champ_Face_VDF.h>
28#include <Domaine_Cl_VDF.h>
29#include <Fluide_Quasi_Compressible.h>
30#include <Debog.h>
31#include <TRUSTTrav.h>
32#include <K_Omega_Eps_constants.h>
33
34////////////////////////////////////////////////////////////////////////////
35//
36// Implementation des fonctions de la classe Calcul_Production_K_VDF
37//
38/////////////////////////////////////////////////////////////////////////////
39DoubleVect& Calcul_Production_K_VDF::
40calculer_terme_production_K(const Domaine_VDF& domaine_VDF, const Domaine_Cl_VDF& domaine_Cl_VDF , DoubleVect& S, const DoubleTab& K_eps,
41 const DoubleTab& vitesse,const Champ_Face_VDF& vit, const DoubleTab& visco_turb ) const
42{
43 int nb_elem = domaine_VDF.domaine().nb_elem();
44 const IntTab& elem_faces = domaine_VDF.elem_faces();
45 const IntVect& orientation = domaine_VDF.orientation();
46
47 int elem;
48 //IntVect element(4);
49
50 S = 0. ;
51
52 // Apres la factorisation du calcul de Sij*Sij (v1.4.6) il est apparu que pour
53 // le terme de production de K la contribution des parois ne doit pas etre
54 // pris en compte, car K et Eps sont ensuite imposes par la loi de paroi,
55 // Si l'on prend en compte la contribution de la paroi les viscosites
56 // turbulentes a la paroi sont tres fortes, et les pas de temps tres faibles
57 // Pour la v1.4.8, on revient a la modelisation de la 1.4.5 qui ne
58 // n'ajoutait pas a Sij*Sij la contribution des parois
59 //Dans la methode calcul_S_barre_sans_contrib_paroi, contribution_paroi
60 //est fixe a 0 par defaut
61
62 vit.calcul_S_barre_sans_contrib_paroi(vitesse,S,domaine_Cl_VDF);
63
64 double coef;
65
66 Debog::verifier("Source_Transport_K_Eps_VDF_P0_VDF:: calculer_terme_production_K ",S);
67 for ( elem=0; elem<nb_elem; elem++)
68 {
69
70 //P= 2*nut*Sij*Sij
71 S(elem) *= visco_turb(elem) ;
72 coef = 0. ;
73 for (int i=0; i<Objet_U::dimension; i++)
74 {
75 coef += (vitesse[elem_faces(elem,i)] - vitesse[elem_faces(elem,i+Objet_U::dimension)]) /domaine_VDF.dim_elem(elem,orientation(elem_faces(elem,i)));
76 }
77 //Corrections pour prendre en compte la divergence de u
78 //non nulle en Quasi-Compressible
79
80 S(elem) += -(2./3.)*visco_turb(elem)*coef * coef;
81 S(elem) += -(2./3.)*K_eps(elem,0) * coef ;
82
83 }
84
85 Debog::verifier("Source_Transport_K_Eps_VDF_P0_VDF:: calculer_terme_production_K ",S);
86 return S ;
87}
88
89DoubleVect& Calcul_Production_K_VDF::
90calculer_terme_production_K_for_komega(const Domaine_VDF& domaine_VDF, const Domaine_Cl_VDF& domaine_Cl_VDF,
91 DoubleVect& S, const DoubleTab& K_Omega,
92 const DoubleTab& vitesse, const Champ_Face_VDF& vit,
93 const DoubleTab& visco_turb,const bool deactivate_production_limiter,const double cst_production_limiter) const
94{
95 const IntTab& elem_faces = domaine_VDF.elem_faces();
96 const IntVect& orientation = domaine_VDF.orientation();
97
98 // cAlan, 23/01/2023. Pour le modele k-omega, on prend un
99 // calcul_S_barre car on n a pas de loi de paroi. La fonction sera a
100 // factoriser avec la precedente.
101 S = 0.;
102 vit.calcul_S_barre(vitesse, S, domaine_Cl_VDF);
103
104 Debog::verifier("Source_Transport_K_Omega_VDF_P0_VDF::calculer_terme_production_K",
105 S);
106
107 int nb_elem = domaine_VDF.nb_elem();
108 for (int elem = 0; elem < nb_elem; elem++)
109 {
110 // P = 2*nut*Sij*Sij
111 S(elem) *= visco_turb(elem);
112 double coef = 0.;
113 for (int i = 0; i < Objet_U::dimension; i++)
114 {
115 coef += (vitesse[elem_faces(elem, i)] - vitesse[elem_faces(elem, i+Objet_U::dimension)])
116 / domaine_VDF.dim_elem(elem, orientation(elem_faces(elem, i)));
117 }
118
119 //Corrections pour prendre en compte la divergence de u
120 //non nulle en Quasi-Compressible
121 S(elem) += -(2./3.)*visco_turb(elem)*coef*coef;
122 S(elem) += -(2./3.)*K_Omega(elem, 0)*coef;
123
124 if (!deactivate_production_limiter)
125 S(elem)=std::min(S(elem),cst_production_limiter*BETA_K*K_Omega(elem,0)*K_Omega(elem,1));
126 }
127
128 Debog::verifier("Source_Transport_K_Omega_VDF_P0_VDF::calculer_terme_production_K_for_komega",
129 S);
130 return S ;
131}
132
133DoubleVect& Calcul_Production_K_VDF::
134calculer_terme_production_K_Axi(const Domaine_VDF& domaine_VDF,
135 const Champ_Face_VDF& vitesse,
136 DoubleVect& P,
137 const DoubleTab& K_Eps,
138 const DoubleTab& visco_turb) const
139{
140 // P est discretise comme K et Eps i.e au centre des elements
141 //
142 // P(elem) = - ( Rey11*tau11 + Rey22*tau22 + Rey33*tau33
143 // + Rey12*(tau12+tau21) + Rey13*(tau13 + tau31)
144 // + Rey23*(tau23+tau32) )
145 //
146 //
147
148 P= 0;
149
150 int nb_elem = domaine_VDF.nb_elem();
151 int nb_aretes = domaine_VDF.nb_aretes();
152 const IntTab& face_voisins = domaine_VDF.face_voisins();
153 const IntTab& Qdm = domaine_VDF.Qdm();
154 const DoubleTab& tau_diag = vitesse.tau_diag();
155 const DoubleTab& tau_croises = vitesse.tau_croises();
156 double d_visco_turb,k_elem,P_arete;
157
158 // Calcul des composantes du tenseur de Reynolds a partir du tenseur GradU
159
160 DoubleTrav Rey_diag(nb_elem,Objet_U::dimension);
161 DoubleTrav Rey_croises(nb_aretes);
162
163 Rey_croises = 0;
164 int num_elem;
165 for (num_elem=0; num_elem<nb_elem; num_elem++)
166 {
167 k_elem = K_Eps(num_elem,0);
168 for (int i=0; i<Objet_U::dimension; i++)
169 Rey_diag(num_elem,i) = -2*visco_turb(num_elem)*tau_diag(num_elem,i)
170 + 2./3.*k_elem;
171 }
172
173 // Boucle sur les aretes bord
174
175 int ndeb,nfin,n_arete;
176 ndeb = domaine_VDF.premiere_arete_bord();
177 nfin = ndeb + domaine_VDF.nb_aretes_bord();
178 int fac1,fac2,fac3,fac4;
179
180 for (n_arete=ndeb; n_arete<nfin; n_arete++)
181 {
182 fac3 = Qdm(n_arete,2);
183
184 d_visco_turb = 0.5*(visco_turb(face_voisins(fac3,0))
185 + visco_turb(face_voisins(fac3,1)));
186 Rey_croises(n_arete) = - d_visco_turb*(tau_croises(n_arete,0) + tau_croises(n_arete,1));
187
188 }
189
190 // Pas de boucle sur les aretes mixtes: les composantes croisees
191 // du tenseur de Reynolds sont prises nulles sur les aretes mixtes
192
193 // Boucle sur les aretes internes
194
195 ndeb = domaine_VDF.premiere_arete_interne();
196 nfin = domaine_VDF.nb_aretes();
197
198 for (n_arete=ndeb; n_arete<nfin; n_arete++)
199 {
200 fac1=Qdm(n_arete,0);
201 fac2=Qdm(n_arete,1);
202 fac3=Qdm(n_arete,2);
203 fac4=Qdm(n_arete,3);
204
205 d_visco_turb = 0.25*(visco_turb(face_voisins(fac3,0)) + visco_turb(face_voisins(fac3,1))
206 + visco_turb(face_voisins(fac4,0)) + visco_turb(face_voisins(fac4,1)));
207 Rey_croises(n_arete) = - d_visco_turb*(tau_croises(n_arete,0) + tau_croises(n_arete,1));
208 }
209
210 // Calcul du terme de production P
211
212 for (num_elem=0; num_elem<nb_elem; num_elem++)
213 {
214 for (int i=0; i<Objet_U::dimension; i++)
215 P(num_elem) -= Rey_diag(num_elem,i)*tau_diag(num_elem,i);
216 }
217
218 ndeb = domaine_VDF.premiere_arete_bord();
219 nfin = ndeb + domaine_VDF.nb_aretes_bord();
220
221 for (n_arete=ndeb; n_arete<nfin; n_arete++)
222 {
223 fac3=Qdm(n_arete,2);
224
225 P_arete = - Rey_croises(n_arete)*(tau_croises(n_arete,0)+tau_croises(n_arete,1));
226 P[face_voisins(fac3,0)] += 0.25*P_arete;
227 P[face_voisins(fac3,1)] += 0.25*P_arete;
228 }
229
230 ndeb = domaine_VDF.premiere_arete_interne();
231 nfin = domaine_VDF.nb_aretes();
232
233 for (n_arete=ndeb; n_arete<nfin; n_arete++)
234 {
235 fac1=Qdm(n_arete,0);
236 fac2=Qdm(n_arete,1);
237
238 P_arete = - Rey_croises(n_arete)*(tau_croises(n_arete,0)+tau_croises(n_arete,1));
239 P[face_voisins(fac1,0)] += 0.25*P_arete;
240 P[face_voisins(fac1,1)] += 0.25*P_arete;
241 P[face_voisins(fac2,0)] += 0.25*P_arete;
242 P[face_voisins(fac2,1)] += 0.25*P_arete;
243 }
244 P.echange_espace_virtuel();
245 return P;
246}
247
248
249DoubleTab& Calcul_Production_K_VDF::calculer_u_teta(const Domaine_VDF& domaine_VDF,
250 const Domaine_Cl_VDF& zcl_VDF,
251 const DoubleTab& temper,
252 const DoubleTab& alpha_turb,
253 DoubleTab& u_teta) const
254{
255 // ---->
256 // On note u_teta le vecteur alpha_turb.gradT
257 //
258 // Sur chaque face on calcule la composante de u_teta normale a la face
259
260 int nb_faces= domaine_VDF.nb_faces();
261 int n0,n1;
262 double alpha,dist;
263 int face;
264 const IntTab& face_voisins = domaine_VDF.face_voisins();
265 u_teta = 0;
266
267 // Traitement des faces internes
268
269
270 int premiere_face_int=domaine_VDF.premiere_face_int();
271 nb_faces=domaine_VDF.nb_faces();
272
273 if (Objet_U::axi)
274
275 for (face=premiere_face_int; face<nb_faces; face++)
276 {
277 n0 = face_voisins(face,0);
278 n1 = face_voisins(face,1);
279 dist = domaine_VDF.dist_norm_axi(face);
280 alpha = 0.5*(alpha_turb(n0)+alpha_turb(n1));
281 u_teta[face] = alpha*(temper[n1] - temper[n0])/dist;
282 }
283 else
284
285 for (face=premiere_face_int; face<nb_faces; face++)
286 {
287 n0 = face_voisins(face,0);
288 n1 = face_voisins(face,1);
289 dist = domaine_VDF.dist_norm(face);
290 alpha = 0.5*(alpha_turb(n0)+alpha_turb(n1));
291 u_teta[face] = alpha*(temper[n1] - temper[n0])/dist;
292 }
293
294 // Traitement des conditions limites de type Entree_fluide_K_Eps_impose :
295
296 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
297 {
298
299 const Cond_lim& la_cl = zcl_VDF.les_conditions_limites(n_bord);
300
301 if (sub_type(Entree_fluide_temperature_imposee,la_cl.valeur()) )
302 {
303 const Entree_fluide_temperature_imposee& la_cl_diri=ref_cast(Entree_fluide_temperature_imposee,la_cl.valeur());
304 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
305 int ndeb = le_bord.num_premiere_face();
306 int nfin = ndeb + le_bord.nb_faces();
307 for (face=ndeb; face<nfin; face++)
308 {
309 double T_imp = la_cl_diri.val_imp(face-ndeb);
310 n0 = face_voisins(face,0);
311 n1 = face_voisins(face,1);
312 if (Objet_U::axi)
313 dist = 2*domaine_VDF.dist_norm_bord_axi(face);
314 else
315 dist = 2*domaine_VDF.dist_norm_bord(face);
316 if (n0 != -1)
317 {
318 alpha = alpha_turb(n0);
319 u_teta[face] = alpha*(T_imp-temper[n0])/dist;
320 }
321 else // n1 != -1
322 {
323 alpha = alpha_turb(n1);
324 u_teta[face] = alpha*(temper[n1]-T_imp)/dist;
325 }
326 }
327 }
328 }
329 return u_teta;
330}
331
332DoubleTab& Calcul_Production_K_VDF::calculer_u_conc(const Domaine_VDF& domaine_VDF,
333 const Domaine_Cl_VDF& zcl_VDF,
334 const DoubleTab& conc,
335 const DoubleTab& alpha_turb,
336 DoubleTab& u_conc,int nb_consti) const
337
338{
339 // ---->
340 // On note u_conc le vecteur alpha_turb.gradC
341 //
342 // Sur chaque face on calcule la composante de u_conc normale a la face
343
344 int nb_faces= domaine_VDF.nb_faces();
345 int n0,n1;
346 double alpha,dist;
347 int face;
348 const IntTab& face_voisins = domaine_VDF.face_voisins();
349 u_conc = 0;
350
351 // Traitement des faces internes
352
353
354 int premiere_face_int=domaine_VDF.premiere_face_int();
355 nb_faces=domaine_VDF.nb_faces();
356 int k;
357
358 if (Objet_U::axi)
359
360 for (face=premiere_face_int; face<nb_faces; face++)
361 {
362 n0 = face_voisins(face,0);
363 n1 = face_voisins(face,1);
364 dist = domaine_VDF.dist_norm_axi(face);
365 alpha = 0.5*(alpha_turb(n0)+alpha_turb(n1));
366 for (k=0; k<nb_consti; k++)
367 u_conc(face,k) = alpha*(conc(n1,k) - conc(n0,k))/dist;
368 }
369 else
370
371 for (face=premiere_face_int; face<nb_faces; face++)
372 {
373 n0 = face_voisins(face,0);
374 n1 = face_voisins(face,1);
375 dist = domaine_VDF.dist_norm(face);
376 alpha = 0.5*(alpha_turb(n0)+alpha_turb(n1));
377 for (k=0; k<nb_consti; k++)
378 u_conc(face,k) = alpha*(conc(n1,k) - conc(n0,k))/dist;
379 }
380
381 // Traitement des conditions limites de type Entree_fluide_K_Eps_impose :
382
383 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
384 {
385
386 const Cond_lim& la_cl = zcl_VDF.les_conditions_limites(n_bord);
387
388 if ( sub_type(Entree_fluide_concentration_imposee,la_cl.valeur()) )
389 {
390 const Entree_fluide_concentration_imposee& la_cl_diri= ref_cast(Entree_fluide_concentration_imposee,la_cl.valeur());
391 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
392 int ndeb = le_bord.num_premiere_face();
393 int nfin = ndeb + le_bord.nb_faces();
394 for (face=ndeb; face<nfin; face++)
395 {
396 n0 = face_voisins(face,0);
397 n1 = face_voisins(face,1);
398 if (Objet_U::axi)
399 dist = 2*domaine_VDF.dist_norm_bord_axi(face);
400 else
401 dist = 2*domaine_VDF.dist_norm_bord(face);
402 if (n0 != -1)
403 {
404 alpha = alpha_turb(n0);
405 for (k=0; k<nb_consti; k++)
406 u_conc(face,k) = alpha*(la_cl_diri.val_imp(face-ndeb,k)-conc(n0,k))/dist;
407 }
408 else // n1 != -1
409 {
410 alpha = alpha_turb(n1);
411 for (k=0; k<nb_consti; k++)
412 u_conc(face,k) = alpha*(conc(n1,k)-la_cl_diri.val_imp(face-ndeb,k))/dist;
413 }
414 }
415 }
416 }
417 return u_conc;
418}
419
420
421DoubleVect& Calcul_Production_K_VDF::calculer_terme_destruction_K(const Domaine_VDF& domaine_VDF,
422 const Domaine_Cl_VDF& zcl_VDF,
423 DoubleVect& G,const DoubleTab& temper,
424 const DoubleTab& alpha_turb,
425 double beta,const DoubleVect& gravite) const
426{
427 //
428 // G est discretise comme K et Eps i.e au centre des elements
429 //
430 // --> ----->
431 // G(elem) = beta alpha_t(elem) G . gradT(elem)
432 //
433
434 int nb_elem = domaine_VDF.nb_elem();
435 int nb_faces= domaine_VDF.nb_faces();
436 DoubleTrav u_teta(nb_faces);
437 const DoubleVect& porosite_face = zcl_VDF.equation().milieu().porosite_face();
438
439 // ---->
440 // On note u_teta le vecteur alpha_turb.gradT
441 //
442 // Appel a la fonction qui calcule sur chaque face la composante
443 // de u_teta normale a la face
444
445 calculer_u_teta(domaine_VDF,zcl_VDF,temper,alpha_turb,u_teta);
446
447 // ------> ----->
448 // On calcule ensuite une valeur moyenne de gravite.u_teta sur chaque
449 // element.
450
451 G = 0;
452
453 // -------> ------>
454 // Calcul de beta.gravite . u_teta
455
456 const Domaine& le_dom=domaine_VDF.domaine();
457 int nb_faces_elem = le_dom.nb_faces_elem();
458
459 IntTrav numfa(nb_faces_elem);
460 DoubleVect coef(Objet_U::dimension);
461 const IntTab& les_elem_faces = domaine_VDF.elem_faces();
462
463 for (int elem=0; elem<nb_elem; elem++)
464 {
465 for (int i=0; i<nb_faces_elem; i++)
466 numfa[i] = les_elem_faces(elem,i);
467
468 coef(0) = 0.5*(u_teta(numfa[0])*porosite_face(numfa[0])
469 + u_teta(numfa[Objet_U::dimension])*porosite_face(numfa[Objet_U::dimension]));
470 coef(1) = 0.5*(u_teta(numfa[1])*porosite_face(numfa[1])
471 + u_teta(numfa[Objet_U::dimension+1])*porosite_face(numfa[Objet_U::dimension+1]));
472
473 if (Objet_U::dimension ==2)
474 G[elem] = beta*(gravite(0)*coef(0) + gravite(1)*coef(1));
475
476 else if (Objet_U::dimension == 3)
477 {
478 coef(2) = 0.5*(u_teta(numfa[2])*porosite_face(numfa[2])
479 + u_teta(numfa[5])*porosite_face(numfa[5]));
480 G[elem] = beta*(gravite(0)*coef(0) + gravite(1)*coef(1) + gravite(2)*coef(2));
481 }
482
483 }
484 G.echange_espace_virtuel();
485 return G;
486}
487
488DoubleVect& Calcul_Production_K_VDF::calculer_terme_destruction_K(const Domaine_VDF& domaine_VDF,
489 const Domaine_Cl_VDF& zcl_VDF,
490 DoubleVect& G,const DoubleTab& temper,
491 const DoubleTab& alpha_turb,
492 const DoubleTab& beta,const DoubleVect& gravite) const
493{
494 //
495 // G est discretise comme K et Eps i.e au centre des elements
496 //
497 // --> ----->
498 // G(elem) = beta(elem) alpha_t(elem) G . gradT(elem)
499 //
500
501 int nb_elem = domaine_VDF.nb_elem();
502 int nb_faces= domaine_VDF.nb_faces();
503 DoubleTrav u_teta(nb_faces);
504 const DoubleVect& porosite_face = zcl_VDF.equation().milieu().porosite_face();
505
506 // ---->
507 // On note u_teta le vecteur alpha_turb.gradT
508 //
509 // Appel a la fonction qui calcule sur chaque face la composante
510 // de u_teta normale a la face
511
512 calculer_u_teta(domaine_VDF,zcl_VDF,temper,alpha_turb,u_teta);
513
514 // ------> ----->
515 // On calcule ensuite une valeur moyenne de gravite.u_teta sur chaque
516 // element.
517
518 G = 0;
519
520 // -------> ------>
521 // Calcul de beta.gravite . u_teta
522
523 const Domaine& le_dom=domaine_VDF.domaine();
524 int nb_faces_elem = le_dom.nb_faces_elem();
525 IntTrav numfa(nb_faces_elem);
526 const IntTab& les_elem_faces = domaine_VDF.elem_faces();
527 DoubleVect coef(Objet_U::dimension);
528
529 for (int elem=0; elem<nb_elem; elem++)
530 {
531 for (int i=0; i<nb_faces_elem; i++)
532 numfa[i] = les_elem_faces(elem,i);
533
534 coef(0) = 0.5*(u_teta(numfa[0])*porosite_face(numfa[0])
535 + u_teta(numfa[Objet_U::dimension])*porosite_face(numfa[Objet_U::dimension]));
536 coef(1) = 0.5*(u_teta(numfa[1])*porosite_face(numfa[1])
537 + u_teta(numfa[Objet_U::dimension+1])*porosite_face(numfa[Objet_U::dimension+1]));
538
539 if (Objet_U::dimension ==2)
540 G[elem] = beta(elem)*(gravite(0)*coef(0) + gravite(1)*coef(1));
541
542 else if (Objet_U::dimension == 3)
543 {
544 coef(2) = 0.5*(u_teta(numfa[2])*porosite_face(numfa[2])
545 + u_teta(numfa[5])*porosite_face(numfa[5]));
546 G[elem] = beta(elem)*(gravite(0)*coef(0) + gravite(1)*coef(1) + gravite(2)*coef(2));
547 }
548
549 }
550 G.echange_espace_virtuel();
551 return G;
552}
553
554DoubleVect& Calcul_Production_K_VDF::calculer_terme_destruction_K(const Domaine_VDF& domaine_VDF,
555 const Domaine_Cl_VDF& zcl_VDF,
556 DoubleVect& G,const DoubleTab& temper,
557 const DoubleTab& alpha_turb,
558 const DoubleVect& beta,
559 const DoubleVect& gravite,
560 int nb_consti) const
561{
562 //
563 // G est discretise comme K et Eps i.e au centre des elements
564 //
565 // G(elem) est la somme sur les nb_consti de Gk(elem,k)
566 //
567 // avec
568 // -> --->
569 // Gk(i,k) = beta_c(k).alpha_t(i).G.gradC(i,k)
570 //
571
572 int nb_elem = domaine_VDF.nb_elem();
573 int nb_faces= domaine_VDF.nb_faces();
574 DoubleTrav u_conc(nb_faces,nb_consti);
575 const DoubleVect& porosite_face = zcl_VDF.equation().milieu().porosite_face();
576
577 // ------> ---->
578 // On note u_conc le vecteur alpha_turb.gradC
579 //
580 // Appel a la fonction qui calcule sur chaque face la composante
581 // de u_conc normale a la face
582
583 calculer_u_conc(domaine_VDF,zcl_VDF,temper,alpha_turb,u_conc,nb_consti);
584
585 // ------> ----->
586 // On calcule ensuite une valeur moyenne de gravite.u_conc sur chaque
587 // element.
588
589 G = 0;
590
591 // -------> ----->
592 // Calcul de beta.gravite . u_conc
593
594 const Domaine& le_dom=domaine_VDF.domaine();
595 int nb_faces_elem = le_dom.nb_faces_elem();
596 IntTrav numfa(nb_faces_elem);
597 const IntTab& les_elem_faces = domaine_VDF.elem_faces();
598 DoubleVect coef(Objet_U::dimension);
599 int k;
600
601 for (int elem=0; elem<nb_elem; elem++)
602 {
603 for (int i=0; i<nb_faces_elem; i++)
604 numfa[i] = les_elem_faces(elem,i);
605
606 for (k=0; k<nb_consti; k++)
607 {
608 coef(0) = 0.5*(u_conc(numfa[0],k)*porosite_face(numfa[0])
609 + u_conc(numfa[Objet_U::dimension],k)*porosite_face(numfa[Objet_U::dimension]));
610 coef(1) = 0.5*(u_conc(numfa[1],k)*porosite_face(numfa[1])
611 + u_conc(numfa[Objet_U::dimension+1],k)*porosite_face(numfa[Objet_U::dimension+1]));
612
613 if (Objet_U::dimension ==2)
614 G[elem] += beta(k)*(gravite(0)*coef(0) + gravite(1)*coef(1));
615 else if (Objet_U::dimension == 3)
616 {
617 coef(2) = 0.5*(u_conc(numfa[2],k)*porosite_face(numfa[2])
618 + u_conc(numfa[5],k)*porosite_face(numfa[5]));
619 G[elem] += beta(k)*(gravite(0)*coef(0) + gravite(1)*coef(1) + gravite(2)*coef(2));
620 }
621 }
622 }
623 G.echange_espace_virtuel();
624 return G;
625}
626
627DoubleVect& Calcul_Production_K_VDF::
628calculer_terme_production_K_BiK(const Domaine_VDF& domaine_VDF, const Domaine_Cl_VDF& domaine_Cl_VDF , DoubleVect& S, const DoubleTab& K,
629 const DoubleTab& vitesse,const Champ_Face_VDF& vit, const DoubleTab& visco_turb ) const
630{
631 int nb_elem = domaine_VDF.domaine().nb_elem();
632 const IntTab& elem_faces = domaine_VDF.elem_faces();
633 const IntVect& orientation = domaine_VDF.orientation();
634
635 int elem;
636 //IntVect element(4);
637
638 S = 0. ;
639
640 // Apres la factorisation du calcul de Sij*Sij (v1.4.6) il est apparu que pour
641 // le terme de production de K la contribution des parois ne doit pas etre
642 // pris en compte, car K et Eps sont ensuite imposes par la loi de paroi,
643 // Si l'on prend en compte la contribution de la paroi les viscosites
644 // turbulentes a la paroi sont tres fortes, et les pas de temps tres faibles
645 // Pour la v1.4.8, on revient a la modelisation de la 1.4.5 qui ne
646 // n'ajoutait pas a Sij*Sij la contribution des parois
647 //Dans la methode calcul_S_barre_sans_contrib_paroi, contribution_paroi
648 //est fixe a 0 par defaut
649
650 vit.calcul_S_barre_sans_contrib_paroi(vitesse,S,domaine_Cl_VDF);
651
652 double coef;
653
654 Debog::verifier("Source_Transport_K_Eps_VDF_P0_VDF:: calculer_terme_production_K_BiK ",S);
655 for ( elem=0; elem<nb_elem; elem++)
656 {
657
658 //P= 2*nut*Sij*Sij
659 S(elem) *= visco_turb(elem) ;
660 coef = 0. ;
661 for (int i=0; i<Objet_U::dimension; i++)
662 {
663 coef += (vitesse[elem_faces(elem,i)] - vitesse[elem_faces(elem,i+Objet_U::dimension)]) /domaine_VDF.dim_elem(elem,orientation(elem_faces(elem,i)));
664 }
665 //Corrections pour prendre en compte la divergence de u
666 //non nulle en Quasi-Compressible
667
668 S(elem) += -(2./3.)*visco_turb(elem)*coef * coef;
669 S(elem) += -(2./3.)*K(elem) * coef ;
670
671 }
672
673 Debog::verifier("Source_Transport_K_Eps_VDF_P0_VDF:: calculer_terme_production_K_BiK ",S);
674 return S ;
675}
676
677DoubleVect& Calcul_Production_K_VDF::
678calculer_terme_production_K_BiK_Axi(const Domaine_VDF& domaine_VDF,
679 const Champ_Face_VDF& vitesse,
680 DoubleVect& P,
681 const DoubleTab& K,
682 const DoubleTab& visco_turb) const
683{
684 // P est discretise comme K et Eps i.e au centre des elements
685 //
686 // P(elem) = - ( Rey11*tau11 + Rey22*tau22 + Rey33*tau33
687 // + Rey12*(tau12+tau21) + Rey13*(tau13 + tau31)
688 // + Rey23*(tau23+tau32) )
689 //
690 //
691
692 P= 0;
693
694 int nb_elem = domaine_VDF.nb_elem();
695 int nb_aretes = domaine_VDF.nb_aretes();
696 const IntTab& face_voisins = domaine_VDF.face_voisins();
697 const IntTab& Qdm = domaine_VDF.Qdm();
698 const DoubleTab& tau_diag = vitesse.tau_diag();
699 const DoubleTab& tau_croises = vitesse.tau_croises();
700 double d_visco_turb,k_elem,P_arete;
701
702 // Calcul des composantes du tenseur de Reynolds a partir du tenseur GradU
703
704 DoubleTrav Rey_diag(nb_elem,Objet_U::dimension);
705 DoubleTrav Rey_croises(nb_aretes);
706
707 Rey_croises = 0;
708 int num_elem;
709 for (num_elem=0; num_elem<nb_elem; num_elem++)
710 {
711 k_elem = K(num_elem);
712 for (int i=0; i<Objet_U::dimension; i++)
713 Rey_diag(num_elem,i) = -2*visco_turb(num_elem)*tau_diag(num_elem,i)
714 + 2./3.*k_elem;
715 }
716
717 // Boucle sur les aretes bord
718
719 int ndeb,nfin,n_arete;
720 ndeb = domaine_VDF.premiere_arete_bord();
721 nfin = ndeb + domaine_VDF.nb_aretes_bord();
722 int fac1,fac2,fac3,fac4;
723
724 for (n_arete=ndeb; n_arete<nfin; n_arete++)
725 {
726 fac3 = Qdm(n_arete,2);
727
728 d_visco_turb = 0.5*(visco_turb(face_voisins(fac3,0))
729 + visco_turb(face_voisins(fac3,1)));
730 Rey_croises(n_arete) = - d_visco_turb*(tau_croises(n_arete,0) + tau_croises(n_arete,1));
731
732 }
733
734 // Pas de boucle sur les aretes mixtes: les composantes croisees
735 // du tenseur de Reynolds sont prises nulles sur les aretes mixtes
736
737 // Boucle sur les aretes internes
738
739 ndeb = domaine_VDF.premiere_arete_interne();
740 nfin = domaine_VDF.nb_aretes();
741
742 for (n_arete=ndeb; n_arete<nfin; n_arete++)
743 {
744 fac1=Qdm(n_arete,0);
745 fac2=Qdm(n_arete,1);
746 fac3=Qdm(n_arete,2);
747 fac4=Qdm(n_arete,3);
748
749 d_visco_turb = 0.25*(visco_turb(face_voisins(fac3,0)) + visco_turb(face_voisins(fac3,1))
750 + visco_turb(face_voisins(fac4,0)) + visco_turb(face_voisins(fac4,1)));
751 Rey_croises(n_arete) = - d_visco_turb*(tau_croises(n_arete,0) + tau_croises(n_arete,1));
752 }
753
754 // Calcul du terme de production P
755
756 for (num_elem=0; num_elem<nb_elem; num_elem++)
757 {
758 for (int i=0; i<Objet_U::dimension; i++)
759 P(num_elem) -= Rey_diag(num_elem,i)*tau_diag(num_elem,i);
760 }
761
762 ndeb = domaine_VDF.premiere_arete_bord();
763 nfin = ndeb + domaine_VDF.nb_aretes_bord();
764
765 for (n_arete=ndeb; n_arete<nfin; n_arete++)
766 {
767 fac3=Qdm(n_arete,2);
768
769 P_arete = - Rey_croises(n_arete)*(tau_croises(n_arete,0)+tau_croises(n_arete,1));
770 P[face_voisins(fac3,0)] += 0.25*P_arete;
771 P[face_voisins(fac3,1)] += 0.25*P_arete;
772 }
773
774 ndeb = domaine_VDF.premiere_arete_interne();
775 nfin = domaine_VDF.nb_aretes();
776
777 for (n_arete=ndeb; n_arete<nfin; n_arete++)
778 {
779 fac1=Qdm(n_arete,0);
780 fac2=Qdm(n_arete,1);
781
782 P_arete = - Rey_croises(n_arete)*(tau_croises(n_arete,0)+tau_croises(n_arete,1));
783 P[face_voisins(fac1,0)] += 0.25*P_arete;
784 P[face_voisins(fac1,1)] += 0.25*P_arete;
785 P[face_voisins(fac2,0)] += 0.25*P_arete;
786 P[face_voisins(fac2,1)] += 0.25*P_arete;
787 }
788 P.echange_espace_virtuel();
789 return P;
790}
Calcul_Production_K_VDF::calculer_u_teta
DoubleTab & calculer_u_teta(const Domaine_VDF &, const Domaine_Cl_VDF &, const DoubleTab &, const DoubleTab &, DoubleTab &) const
Definition Calcul_Production_K_VDF.cpp:249
Calcul_Production_K_VDF::calculer_terme_production_K
DoubleVect & calculer_terme_production_K(const Domaine_VDF &, const Domaine_Cl_VDF &, DoubleVect &, const DoubleTab &, const DoubleTab &, const Champ_Face_VDF &, const DoubleTab &) const
Definition Calcul_Production_K_VDF.cpp:40
Calcul_Production_K_VDF::calculer_terme_production_K_BiK
DoubleVect & calculer_terme_production_K_BiK(const Domaine_VDF &, const Domaine_Cl_VDF &, DoubleVect &, const DoubleTab &, const DoubleTab &, const Champ_Face_VDF &, const DoubleTab &) const
Definition Calcul_Production_K_VDF.cpp:628
Calcul_Production_K_VDF::calculer_terme_production_K_BiK_Axi
DoubleVect & calculer_terme_production_K_BiK_Axi(const Domaine_VDF &, const Champ_Face_VDF &, DoubleVect &, const DoubleTab &, const DoubleTab &) const
Definition Calcul_Production_K_VDF.cpp:678
Calcul_Production_K_VDF::calculer_terme_production_K_for_komega
DoubleVect & calculer_terme_production_K_for_komega(const Domaine_VDF &, const Domaine_Cl_VDF &, DoubleVect &, const DoubleTab &, const DoubleTab &, const Champ_Face_VDF &, const DoubleTab &, const bool deactivate_production_limiter, const double cst_production_limiter) const
Definition Calcul_Production_K_VDF.cpp:90
Calcul_Production_K_VDF::calculer_u_conc
DoubleTab & calculer_u_conc(const Domaine_VDF &, const Domaine_Cl_VDF &, const DoubleTab &, const DoubleTab &, DoubleTab &, int) const
Definition Calcul_Production_K_VDF.cpp:332
Calcul_Production_K_VDF::calculer_terme_destruction_K
DoubleVect & calculer_terme_destruction_K(const Domaine_VDF &, const Domaine_Cl_VDF &, DoubleVect &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleVect &) const
Definition Calcul_Production_K_VDF.cpp:488
Calcul_Production_K_VDF::calculer_terme_production_K_Axi
DoubleVect & calculer_terme_production_K_Axi(const Domaine_VDF &, const Champ_Face_VDF &, DoubleVect &, const DoubleTab &, const DoubleTab &) const
Definition Calcul_Production_K_VDF.cpp:134
Champ_Face_VDF
class Champ_Face_VDF Cette classe sert a representer un champ vectoriel dont on ne calcule
Definition Champ_Face_VDF.h:42
Champ_Face_VDF::calcul_S_barre
DoubleVect & calcul_S_barre(const DoubleTab &, DoubleVect &, const Domaine_Cl_VDF &) const
Definition Champ_Face_VDF.cpp:1025
Champ_Face_VDF::calcul_S_barre_sans_contrib_paroi
DoubleVect & calcul_S_barre_sans_contrib_paroi(const DoubleTab &, DoubleVect &, const Domaine_Cl_VDF &) const
Methode qui renvoie SMA_barre aux elements a partir de la vitesse aux faces.
Definition Champ_Face_VDF.cpp:834
Champ_Face_VDF::tau_croises
const DoubleTab & tau_croises() const
Definition Champ_Face_VDF.h:56
Champ_Face_VDF::tau_diag
const DoubleTab & tau_diag() const
Definition Champ_Face_VDF.h:55
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Debog::verifier
static void verifier(const char *const msg, double)
Definition Debog.cpp:21
Dirichlet::val_imp
virtual double val_imp(int i) const
Renvoie la valeur imposee sur la i-eme composante du champ a la frontiere au temps par defaut du cham...
Definition Dirichlet.cpp:35
Domaine_32_64::nb_elem
int_t nb_elem() const
Definition Domaine.h:131
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Renvoie le nombre de face de type i des elements geometriques constituants le domaine.
Definition Domaine.h:484
Domaine_Cl_VDF
class Domaine_Cl_VDF
Definition Domaine_Cl_VDF.h:63
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Renvoie la i-ieme condition aux limites.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VDF
class Domaine_VDF
Definition Domaine_VDF.h:64
Domaine_VDF::dist_norm
double dist_norm(int num_face) const override
Definition Domaine_VDF.h:342
Domaine_VDF::dist_norm_bord_axi
double dist_norm_bord_axi(int num_face) const
Definition Domaine_VDF.h:400
Domaine_VDF::nb_aretes
int nb_aretes() const
Definition Domaine_VDF.h:570
Domaine_VDF::dim_elem
double dim_elem(int, int) const
Definition Domaine_VDF.h:627
Domaine_VDF::orientation
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
Definition Domaine_VDF.h:297
Domaine_VDF::premiere_arete_bord
int premiere_arete_bord() const
Definition Domaine_VDF.h:538
Domaine_VDF::premiere_arete_interne
int premiere_arete_interne() const
Definition Domaine_VDF.h:578
Domaine_VDF::nb_aretes_bord
int nb_aretes_bord() const
Definition Domaine_VDF.h:530
Domaine_VDF::dist_norm_axi
double dist_norm_axi(int num_face) const
Definition Domaine_VDF.h:356
Domaine_VDF::Qdm
int Qdm(int num_arete, int) const
Definition Domaine_VDF.h:231
Domaine_VDF::dist_norm_bord
double dist_norm_bord(int num_face) const override
Definition Domaine_VDF.h:382
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::elem_faces
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
Definition Domaine_VF.h:543
Domaine_VF::premiere_face_int
int premiere_face_int() const
une face est interne ssi elle separe deux elements.
Definition Domaine_VF.h:463
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_dis_base::nb_elem
int nb_elem() const
Definition Domaine_dis_base.h:49
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Entree_fluide_concentration_imposee
Definition Entree_fluide_concentration_imposee.h:22
Entree_fluide_temperature_imposee
classe Entree_fluide_temperature_imposee Cas particulier de la classe Dirichlet_entree_fluide pour la...
Definition Dirichlet_entree_fluide_leaves.h:86
Equation_base::milieu
virtual const Milieu_base & milieu() const =0
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_premiere_face
int num_premiere_face() const
Definition Front_VF.h:63
Milieu_base::porosite_face
DoubleVect & porosite_face()
Definition Milieu_base.h:62
MorEqn::equation
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::axi
static int axi
Definition Objet_U.h:101
TRUSTVect::echange_espace_virtuel
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
Definition TRUSTVect.tpp:282