TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Op_Dift_VDF_Face_Axi_base.cpp
1/****************************************************************************
2* Copyright (c) 2026, 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 <Op_Dift_VDF_Face_Axi_base.h>
17#include <Modele_turbulence_hyd_base.h>
18
19Implemente_base(Op_Dift_VDF_Face_Axi_base,"Op_Dift_VDF_Face_Axi_base",Op_Dift_VDF_Face_base);
20
21Sortie& Op_Dift_VDF_Face_Axi_base::printOn(Sortie& s ) const { return s << que_suis_je() ; }
23
24void Op_Dift_VDF_Face_Axi_base::associer(const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_cl_dis, const Champ_Inc_base& ch_transporte)
25{
26 const Domaine_VDF& zvdf = ref_cast(Domaine_VDF,domaine_dis);
27 const Domaine_Cl_VDF& zclvdf = ref_cast(Domaine_Cl_VDF,domaine_cl_dis);
28 const Champ_Face_VDF& inco = ref_cast(Champ_Face_VDF,ch_transporte);
29 le_dom_vdf = zvdf;
30 la_zcl_vdf = zclvdf;
31 inconnue = inco;
32 surface.ref(zvdf.face_surfaces());
34 orientation.ref(zvdf.orientation());
35 porosite.ref(la_zcl_vdf->equation().milieu().porosite_face());
36 xp.ref(zvdf.xp());
37 xv.ref(zvdf.xv());
38 Qdm.ref(zvdf.Qdm());
39 face_voisins.ref(zvdf.face_voisins());
40 elem_faces.ref(zvdf.elem_faces());
41 type_arete_bord.ref(zclvdf.type_arete_bord());
42}
43
45{
47 Equation_base& eqn_hydr = equation();
48 Champ_Face_VDF& vitesse = ref_cast(Champ_Face_VDF,eqn_hydr.inconnue());
50 const RefObjU& modele_turbulence = eqn_hydr.get_modele(TURBULENCE);
51 const Modele_turbulence_hyd_base& mod_turb = ref_cast(Modele_turbulence_hyd_base,modele_turbulence.valeur());
53}
54
56{
57 return Op_Dift_VDF_Face_base::calculer_dt_stab(le_dom_vdf.valeur()) ;
58}
59
61{
62 if (le_modele_turbulence->has_loi_paroi_hyd()) tau_tan.ref(le_modele_turbulence->loi_paroi().Cisaillement_paroi());
63}
64
65// XXX E Saikali : j'ai fait comme ca sinon nu_t est pas initialiser dans le cas var
67{
68 le_modele_turbulence = mod;
69 associer_diffusivite_turbulente(le_modele_turbulence->viscosite_turbulente());
70 associer_loipar(le_modele_turbulence->loi_paroi()); /* on fait rien pour le moment ... */
71}
72
73void Op_Dift_VDF_Face_Axi_base::ajouter_elem(const DoubleVect& visco_turb, const DoubleTab& tau_diag, DoubleTab& resu) const
74{
75 for (int num_elem = 0; num_elem < le_dom_vdf->nb_elem(); num_elem++)
76 {
77 const int fx0 = elem_faces(num_elem,0), fx1 = elem_faces(num_elem,dimension), fy0 = elem_faces(num_elem,1), fy1 = elem_faces(num_elem,1+dimension);
78 const double visc_elem = nu_(num_elem) + 2*visco_turb(num_elem);
79 const double flux_X = (visc_elem*tau_diag(num_elem,0))*0.5*(surface(fx0)+surface(fx1)), flux_Y = (visc_elem*tau_diag(num_elem,1))*0.5*(surface(fy0)+surface(fy1));
80 resu(fx0) += flux_X;
81 resu(fx1) += flux_X; // TODO FIXME EUH .?.?. Yannick help :/
82 resu(fy0) += flux_Y;
83 resu(fy1) += flux_Y; // TODO FIXME EUH .?.?. Yannick help :/
84
85 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
86 const double coef_laplacien_axi = -0.5*(tau_diag(num_elem,1)*visc_elem);
87 resu[fx0] += coef_laplacien_axi*volumes_entrelaces(fx0)*porosite(fx0)/xv(fx0,0);
88 resu[fx1] += coef_laplacien_axi*volumes_entrelaces(fx1)*porosite(fx1)/xv(fx1,0);
89 }
90}
91
92void Op_Dift_VDF_Face_Axi_base::ajouter_elem_3D(const DoubleVect& visco_turb, const DoubleTab& tau_diag, DoubleTab& resu) const
93{
94 for (int num_elem = 0; num_elem < le_dom_vdf->nb_elem(); num_elem++)
95 {
96 const int fz0 = elem_faces(num_elem,2), fz1 = elem_faces(num_elem,2+dimension);
97 const double visc_elem = nu_(num_elem) + 2*visco_turb(num_elem), flux_Z = (visc_elem*tau_diag(num_elem,2))*0.5*(surface(fz0)+surface(fz1));
98 resu[fz0] += flux_Z;
99 resu[fz1] -= flux_Z;
100 }
101}
102
103void Op_Dift_VDF_Face_Axi_base::ajouter_aretes_bords(const DoubleVect& visco_turb, const DoubleTab& tau_croises, DoubleTab& resu) const
104{
105 const int ndeb = le_dom_vdf->premiere_arete_bord(), nfin = ndeb + le_dom_vdf->nb_aretes_bord();
106 double d_visco_turb, d_visco_lam;
107 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
108 {
109 const int n_type = type_arete_bord(n_arete-ndeb);
110 switch(n_type)
111 {
112 case TypeAreteBordVDF::PAROI_PAROI: // paroi-paroi
113 {
114 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), ori3 = orientation(fac3);
115 const int rang1 = (fac1 - le_dom_vdf->premiere_face_bord()), rang2 = (fac2-le_dom_vdf->premiere_face_bord());
116 double coef;
117 if (is_var()) // XXX : E Saikali : sais pas quoi faire sinon ecarts ...
118 {
119 // Calcul du frottement identique a celui de TRIOVF : On calcule la moyenne des u_star et on l'eleve au carre. On calcule la moyenne des surfaces
120 const double tau_tan_1 = tau_tan(rang1,ori3), tau_tan_2 = tau_tan(rang2,ori3) ;
121 double tau = 0.5*(tau_tan_1 + tau_tan_2 ), surf = 0.5*(surface(fac1)+surface(fac2));
122 coef = tau*tau*surf;
123 }
124 else // Autre solution pour le calcul du frottement : On calcule u_star*u_star*surf sur chaque partie de la facette de Qdm
125 {
126 const double tau1 = tau_tan(rang1,ori3)*0.5*surface(fac1), tau2 = tau_tan(rang2,ori3)*0.5*surface(fac2);
127 coef = tau1+tau2;
128 }
129
130 resu[fac3] += coef;
131 break;
132 }
133 case TypeAreteBordVDF::FLUIDE_FLUIDE: /* fall through */
135 {
136 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), signe = Qdm(n_arete,3), ori1b = orientation(fac1), ori3b = orientation(fac3);
137 d_visco_turb = 0.5*(visco_turb(face_voisins(fac3,0)) + visco_turb(face_voisins(fac3,1)));
138 d_visco_lam = nu_mean_2_pts_(face_voisins(fac3,0),face_voisins(fac3,1));
139
140 if (ori1b == 0) // bord d'equation R = cte
141 {
142 double flux1;
143 if (ori3b == 1)
144 {
145 const double tau12 = tau_croises(n_arete,0), tau21 = tau_croises(n_arete,1);
146 flux1 = (d_visco_lam*tau12 + d_visco_turb*(tau12+tau21))*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
147 }
148 else //if (ori3 == 2)
149 {
150 assert (ori3b == 2) ;
151 const double tau13 = tau_croises(n_arete,0), tau31 = tau_croises(n_arete,1);
152 flux1 = (d_visco_lam*tau13 + d_visco_turb*(tau13+tau31))*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
153 }
154
155 resu[fac3] += signe*flux1;
156 }
157 else if (ori1b == 1) // bord d'equation teta = cte
158 {
159 double flux2;
160 if (ori3b == 0)
161 {
162 const double tau21 = tau_croises(n_arete,0), tau12 = tau_croises(n_arete,1);
163 flux2 = (d_visco_lam*tau21 + d_visco_turb*(tau12+tau21))*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
164
165 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
166 const double coef_laplacien_axi = 0.5*(d_visco_lam*tau21 + d_visco_turb*(tau12+tau21));
167 resu(fac1) += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
168 resu(fac2) += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
169 }
170 else // if (ori3b == 2) flux de tau23 a travers le bord
171 {
172 assert (ori3b == 2);
173 const double tau23 = tau_croises(n_arete,0), tau32 = tau_croises(n_arete,1);
174 flux2 = (d_visco_lam*tau23 + d_visco_turb*(tau23+tau32))*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
175 }
176
177 resu[fac3] += signe*flux2;
178 }
179 else // (ori1 == 2) bord d'equation Z = cte
180 {
181 double flux3;
182 if (ori3b == 0)
183 {
184 const double tau31 = tau_croises(n_arete,0), tau13 = tau_croises(n_arete,1);
185 flux3 = (d_visco_lam*tau31 + d_visco_turb*(tau13+tau31))*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
186 }
187 else //if (ori3 == 1) flux de tau32 a travers le bord
188 {
189 assert(ori3b == 1);
190 const double tau32 = tau_croises(n_arete,0), tau23 = tau_croises(n_arete,1);
191 flux3 = (d_visco_lam*tau32 + d_visco_turb*(tau23+tau32))*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
192 }
193
194 resu[fac3] += signe*flux3;
195 }
196 break;
197 }
199 break;
200 default :
201 {
202 Cerr << "On a rencontre un type d'arete non prevu : [ num arete : " << n_arete << " ], [ type : " << n_type << " ]" << finl;
204 break;
205 }
206 }
207 }
208}
209
210void Op_Dift_VDF_Face_Axi_base::fill_resu_aretes_mixtes(const int i , const int j , const int k , const int l , const double d_visco_lam , const double tau , DoubleTab& resu) const
211{
212 // flux de mu_lam*tau sur la facette a cheval sur les faces i et j
213 const double flux = d_visco_lam*tau*0.25*(surface(i)+surface(j))*(porosite(i)+porosite(j));
214 resu(k) += flux;
215 resu(l) -= flux;
216}
217
218void Op_Dift_VDF_Face_Axi_base::ajouter_aretes_mixtes(const DoubleTab& tau_croises, DoubleTab& resu) const
219{
220 // Sur les aretes mixtes les termes croises du tenseur de Reynolds sont nuls: il ne reste donc que la diffusion laminaire
221 const int ndeb = le_dom_vdf->premiere_arete_mixte(), nfin = le_dom_vdf->premiere_arete_interne();
222 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
223 {
224 const int fac1=Qdm(n_arete,0), fac2=Qdm(n_arete,1), fac3=Qdm(n_arete,2), fac4=Qdm(n_arete,3), ori1 = orientation(fac1), ori3 = orientation(fac3);
225 const double d_visco_lam = nu_mean_4_pts_(fac3,fac4); // XXX : BUG corrige dans cette fonction si ecart ...
226
227 if (ori1 == 1) // (seule possibilite : ori3 =0) Arete XY
228 {
229 const double tau12 = tau_croises(n_arete,0), tau21 = tau_croises(n_arete,1);
230 fill_resu_aretes_mixtes(fac1,fac2,fac3,fac4,d_visco_lam,tau21,resu);
231
232 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
233 const double coef_laplacien_axi = 0.5*d_visco_lam*tau21;
234 resu(fac1) += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
235 resu(fac2) += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
236
237 fill_resu_aretes_mixtes(fac3,fac4,fac1,fac2,d_visco_lam,tau12,resu);
238 }
239 else if (ori3 == 1) // (seule possibilite ori1 = 2) arete YZ
240 {
241 const double tau23 = tau_croises(n_arete,0), tau32 = tau_croises(n_arete,1);
242 fill_resu_aretes_mixtes(fac1,fac2,fac3,fac4,d_visco_lam,tau32,resu);
243 fill_resu_aretes_mixtes(fac3,fac4,fac1,fac2,d_visco_lam,tau23,resu);
244 }
245 else // seule possibilite ori1 = 2 et ori3 = 0: arete XZ
246 {
247 const double tau13 = tau_croises(n_arete,0), tau31 = tau_croises(n_arete,1);
248 fill_resu_aretes_mixtes(fac1,fac2,fac3,fac4,d_visco_lam,tau31,resu);
249 fill_resu_aretes_mixtes(fac3,fac4,fac1,fac2,d_visco_lam,tau13,resu);
250 }
251 }
252}
253
254void Op_Dift_VDF_Face_Axi_base::fill_resu_aretes_internes(const int i, const int j, const int k, const int l, const double v_lam, const double v_turb, const double tau1, const double tau2, DoubleTab& resu) const
255{
256 // flux de v_lam*tau2 + v_turb*(tau2+tau1) sur la facette a cheval sur les faces i et j
257 const double flux = (v_lam*tau2 + v_turb*(tau2+tau1))*0.25*(surface(i)+surface(j))*(porosite(i)+porosite(j));
258 resu(k) += flux;
259 resu(l) -= flux;
260}
261
262void Op_Dift_VDF_Face_Axi_base::ajouter_aretes_internes(const DoubleVect& visco_turb, const DoubleTab& tau_croises, DoubleTab& resu) const
263{
264 const int ndeb = le_dom_vdf->premiere_arete_interne(), nfin = le_dom_vdf->nb_aretes();
265 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
266 {
267 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), fac4 = Qdm(n_arete,3), ori1 = orientation(fac1), ori3 = orientation(fac3);
268
269 // Calcul de la viscosite turbulente au milieu de l'arete
270 const double d_visco_turb = 0.25*(visco_turb(face_voisins(fac3,0)) + visco_turb(face_voisins(fac3,1)) + visco_turb(face_voisins(fac4,0)) + visco_turb(face_voisins(fac4,1)));
271 const double d_visco_lam = nu_mean_4_pts_(face_voisins(fac3,0), face_voisins(fac3,1), face_voisins(fac4,0), face_voisins(fac4,1));
272
273 if (ori1 == 1) // (seule possibilite : ori3 =0) Arete XY
274 {
275 const double tau12 = tau_croises(n_arete,0), tau21 = tau_croises(n_arete,1);
276 fill_resu_aretes_internes(fac1,fac2,fac3,fac4,d_visco_lam,d_visco_turb,tau12,tau21,resu);
277
278 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
279 const double coef_laplacien_axi = 0.5*(d_visco_lam*tau21 + d_visco_turb*(tau21+tau12));
280 resu(fac1) += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
281 resu(fac2) += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
282
283 fill_resu_aretes_internes(fac3,fac4,fac1,fac2,d_visco_lam,d_visco_turb,tau21,tau12,resu);
284 }
285 else if (ori3 == 1) // (seule possibilite ori1 = 2) arete YZ
286 {
287 const double tau23 = tau_croises(n_arete,0), tau32 = tau_croises(n_arete,1);
288 fill_resu_aretes_internes(fac1,fac2,fac3,fac4,d_visco_lam,d_visco_turb,tau23,tau32,resu);
289 fill_resu_aretes_internes(fac3,fac4,fac1,fac2,d_visco_lam,d_visco_turb,tau32,tau23,resu);
290 }
291 else // seule possibilite ori1 = 2 et ori3 = 0: arete XZ
292 {
293 const double tau13 = tau_croises(n_arete,0), tau31 = tau_croises(n_arete,1);
294 fill_resu_aretes_internes(fac1,fac2,fac3,fac4,d_visco_lam,d_visco_turb,tau13,tau31,resu);
295 fill_resu_aretes_internes(fac3,fac4,fac1,fac2,d_visco_lam,d_visco_turb,tau31,tau13,resu);
296 }
297 }
298}
299
300DoubleTab& Op_Dift_VDF_Face_Axi_base::ajouter(const DoubleTab& inco, DoubleTab& resu) const
301{
302 if (inco.line_size() > 1) not_implemented(__func__);
303 const double temps = equation().schema_temps().temps_courant();
304 mettre_a_jour_var(temps); // seulement pour var_axi !
305
306 const Domaine_Cl_VDF& zclvdf = la_zcl_vdf.valeur();
307 const DoubleVect& visco_turb = diffusivite_turbulente().valeurs();
308 const DoubleTab& tau_diag = inconnue->tau_diag(), &tau_croises = inconnue->tau_croises();
309 ref_cast_non_const(Champ_Face_VDF,inconnue.valeur()).calculer_dercov_axi(zclvdf);
310
311 ajouter_elem(visco_turb,tau_diag,resu); // Boucle sur les elements pour traiter les facettes situees a l'interieur des elements
312 if (dimension == 3) ajouter_elem_3D(visco_turb,tau_diag,resu);
313 ajouter_aretes_bords(visco_turb,tau_croises,resu); // Boucle sur les aretes bord
314 ajouter_aretes_mixtes(tau_croises,resu); // Boucle sur les aretes mixtes
315 ajouter_aretes_internes(visco_turb,tau_croises,resu); // Boucle sur les aretes internes
316
317 return resu;
318}
319
320DoubleTab& Op_Dift_VDF_Face_Axi_base::calculer(const DoubleTab& inco, DoubleTab& resu) const
321{
322 resu = 0;
323 return ajouter(inco,resu);
324}
325
326void Op_Dift_VDF_Face_Axi_base::fill_coeff_matrice_morse(const int fac1, const int fac2, const double flux, Matrice_Morse& matrice) const
327{
328 const auto& tab1 = matrice.get_set_tab1();
329 const auto& tab2 = matrice.get_set_tab2();
330 auto& coeff = matrice.get_set_coeff();
331 for (auto k = tab1[fac1]-1; k < tab1[fac1+1]-1; k++)
332 {
333 if (tab2[k]-1 == fac1) coeff[k] += flux;
334 if (tab2[k]-1 == fac2) coeff[k] -= flux;
335 }
336 for (auto k = tab1[fac2]-1; k < tab1[fac2+1]-1; k++)
337 {
338 if (tab2[k]-1 == fac1) coeff[k] -= flux;
339 if (tab2[k]-1 == fac2) coeff[k] += flux;
340 }
341}
342
343void Op_Dift_VDF_Face_Axi_base::ajouter_contribution_elem(const DoubleVect& visco_turb, const DoubleTab& tau_diag, Matrice_Morse& matrice) const
344{
345 auto& tab1 = matrice.get_set_tab1();
346 auto& tab2 = matrice.get_set_tab2();
347 auto& coeff = matrice.get_set_coeff();
348 for (int num_elem = 0; num_elem < le_dom_vdf->nb_elem(); num_elem++)
349 {
350 const int fx0 = elem_faces(num_elem,0), fx1 = elem_faces(num_elem,dimension), fy0 = elem_faces(num_elem,1), fy1 = elem_faces(num_elem,1+dimension);
351 const double visc_elem = nu_(num_elem) + 2*visco_turb(num_elem);
352 double flux_X, flux_Y;
353
354 if (is_var())
355 {
356 flux_X = (visc_elem*tau_diag(num_elem,0))*0.5*(surface(fx0)+surface(fx1));
357 flux_Y = (visc_elem*tau_diag(num_elem,1))*0.5*(surface(fy0)+surface(fy1));
358 }
359 else
360 {
361 flux_X = visc_elem*0.5*(surface(fx0)+surface(fx1));
362 flux_Y = visc_elem*0.5*(surface(fy0)+surface(fy1));
363 }
364
365 fill_coeff_matrice_morse(fx0,fx1,flux_X,matrice);
366 fill_coeff_matrice_morse(fy0,fy1,flux_Y,matrice);
367
368 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
369 double coef_laplacien_axi;
370 if (is_var()) coef_laplacien_axi = -0.5*(tau_diag(num_elem,1)*visc_elem);
371 else coef_laplacien_axi = -0.5*visc_elem; // XXX : comprends rien la
372
373 for (auto l = tab1[fx0]-1; l < tab1[fx0+1]-1; l++)
374 if (tab2[l]-1 == fx0) coeff[l] += coef_laplacien_axi*volumes_entrelaces(fx0)*porosite(fx0)/xv(fx0,0);
375
376 for (auto l = tab1[fx1]-1; l < tab1[fx1+1]-1; l++)
377 if (tab2[l]-1 == fx1) coeff[l] += coef_laplacien_axi*volumes_entrelaces(fx1)*porosite(fx1)/xv(fx1,0);
378 }
379}
380
381void Op_Dift_VDF_Face_Axi_base::ajouter_contribution_elem_3D(const DoubleVect& visco_turb, const DoubleTab& tau_diag, Matrice_Morse& matrice) const
382{
383 for (int num_elem = 0; num_elem < le_dom_vdf->nb_elem(); num_elem++)
384 {
385 const int fz0 = elem_faces(num_elem,2), fz1 = elem_faces(num_elem,2+dimension);
386 const double visc_elem = nu_(num_elem) + 2*visco_turb(num_elem);
387 double flux_Z;
388
389 if (is_var()) flux_Z = (visc_elem*tau_diag(num_elem,2))*0.5*(surface(fz0)+surface(fz1));
390 else flux_Z = visc_elem*0.5*(surface(fz0)+surface(fz1));
391
392 fill_coeff_matrice_morse(fz0,fz1,flux_Z,matrice);
393 }
394}
395
396void Op_Dift_VDF_Face_Axi_base::ajouter_contribution_aretes_bords(const DoubleVect& visco_turb, const DoubleTab& tau_diag, Matrice_Morse& matrice) const
397{
398 auto& tab1 = matrice.get_set_tab1();
399 auto& tab2 = matrice.get_set_tab2();
400 auto& coeff = matrice.get_set_coeff();
401 int ndeb = le_dom_vdf->premiere_arete_bord(), nfin = ndeb + le_dom_vdf->nb_aretes_bord();
402 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
403 {
404 const int n_type = type_arete_bord(n_arete-ndeb);
405 switch(n_type)
406 {
408 break;
409 case TypeAreteBordVDF::FLUIDE_FLUIDE: /* fall through */
411 {
412 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), signe = Qdm(n_arete,3), ori1b = orientation(fac1), ori3b = orientation(fac3);
413 const double d_visco_turb = 0.5*(visco_turb(face_voisins(fac3,0))+ visco_turb(face_voisins(fac3,1)));
414 const double d_visco_lam = nu_mean_2_pts_(face_voisins(fac3,0),face_voisins(fac3,1));
415
416 if (ori1b == 0) // bord d'equation R = cte
417 {
418 // XXX j'ai supprime le if
419 const double flux1 = (d_visco_lam+ d_visco_turb)*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
420 for (auto l = tab1[fac3]-1; l < tab1[fac3+1]-1; l++)
421 if (tab2[l]-1 == fac3) coeff[l] += signe*flux1;
422 }
423 else if (ori1b == 1) // bord d'equation teta = cte
424 {
425 if (ori3b == 0)
426 {
427 const double flux2 = (d_visco_lam + d_visco_turb)*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
428 for (auto l = tab1[fac3]-1; l < tab1[fac3+1]-1; l++)
429 if (tab2[l]-1 == fac3) coeff[l] += signe*flux2;
430
431 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
432 const double coef_laplacien_axi = 0.5*(d_visco_lam + d_visco_turb);
433
434 for (auto l = tab1[fac1]-1; l < tab1[fac1+1]-1; l++)
435 if (tab2[l]-1 == fac1) coeff[l] += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
436
437 for (auto l = tab1[fac2]-1; l < tab1[fac2+1]-1; l++)
438 if (tab2[l]-1 == fac2) coeff[l] += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
439 }
440 else if (ori3b == 2) // flux de tau23 a travers le bord
441 {
442 const double flux3 = (d_visco_lam + d_visco_turb)*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
443 for (auto l = tab1[fac3]-1; l < tab1[fac3+1]-1; l++)
444 if (tab2[l]-1 == fac3) coeff[l] += signe*flux3;
445 }
446 }
447 else // (ori1 == 2) bord d'equation Z = cte
448 {
449 const double flux4 = (d_visco_lam + d_visco_turb)*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
450 for (auto l = tab1[fac3]-1; l < tab1[fac3+1]-1; l++)
451 if (tab2[l]-1 == fac3) coeff[l] += signe*flux4;
452 }
453 break;
454 }
456 break;
457 default :
458 {
459 Cerr << "On a rencontre un type d'arete non prevu : [ num arete : " << n_arete << " ], [ type : " << n_type << " ]" << finl;
461 break;
462 }
463 }
464 }
465}
466
467void Op_Dift_VDF_Face_Axi_base::ajouter_contribution_aretes_mixtes(Matrice_Morse& matrice) const
468{
469 auto& tab1 = matrice.get_set_tab1();
470 auto& tab2 = matrice.get_set_tab2();
471 auto& coeff = matrice.get_set_coeff();
472 // Sur les aretes mixtes les termes croises du tenseur de Reynolds sont nuls: il ne reste donc que la diffusion laminaire
473 const int ndeb = le_dom_vdf->premiere_arete_mixte(), nfin = le_dom_vdf->premiere_arete_interne();
474 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
475 {
476 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), fac4 = Qdm(n_arete,3), ori1 = orientation(fac1);
477 const double d_visco_lam = nu_mean_4_pts_(fac3,fac4); // XXX : BUG corrige dans cette fonction si ecart ...
478
479 if (ori1 == 1) // (seule possibilite : ori3 =0) Arete XY
480 {
481 double flux1;
482 // flux de mu_lam*tau21 sur la facette a cheval sur les faces fac1 et fac2
483 flux1 = d_visco_lam*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
484 fill_coeff_matrice_morse(fac3,fac4,flux1,matrice);
485
486 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
487 const double coef_laplacien_axi = 0.5*d_visco_lam;
488
489 for (auto l = tab1[fac1]-1; l < tab1[fac1+1]-1; l++)
490 if (tab2[l]-1 == fac1) coeff[l] += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
491
492 for (auto l = tab1[fac2]-1; l < tab1[fac2+1]-1; l++)
493 if (tab2[l]-1 == fac2) coeff[l] += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
494
495 // flux de mu_lam*tau12 sur la facette a cheval sur les faces fac3 et fac4
496 flux1 = d_visco_lam*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
497 fill_coeff_matrice_morse(fac1,fac2,flux1,matrice);
498 }
499 else // 2 possibilites : ori1 = 2 et ori3 == 1: arete YZ OU ori1 = 2 et ori3 = 0: arete XZ
500 {
501 // flux de mu_lam*tau32 sur la facette a cheval sur les faces fac1 et fac2 ou
502 // flux de mu_lam*tau31 sur la facette a cheval sur les faces fac1 et fac2
503
504 double flux2;
505 flux2 = d_visco_lam*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
506 fill_coeff_matrice_morse(fac3,fac4,flux2,matrice);
507
508 // flux de mu_lam*tau23 sur la facette a cheval sur les faces fac3 et fac4 ou
509 // flux de mu_lam*tau13 sur la facette a cheval sur les faces fac3 et fac4
510 flux2 = d_visco_lam*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
511 fill_coeff_matrice_morse(fac1,fac2,flux2,matrice);
512 }
513 }
514}
515
516void Op_Dift_VDF_Face_Axi_base::ajouter_contribution_aretes_internes(const DoubleVect& visco_turb, Matrice_Morse& matrice) const
517{
518 auto& tab1 = matrice.get_set_tab1();
519 auto& tab2 = matrice.get_set_tab2();
520 auto& coeff = matrice.get_set_coeff();
521 const int ndeb = le_dom_vdf->premiere_arete_interne(), nfin = le_dom_vdf->nb_aretes();
522 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
523 {
524 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), fac4 = Qdm(n_arete,3), ori1 = orientation(fac1);
525 // Calcul de la viscosite turbulente au milieu de l'arete
526 const double d_visco_turb = 0.25*(visco_turb(face_voisins(fac3,0)) + visco_turb(face_voisins(fac3,1)) + visco_turb(face_voisins(fac4,0)) + visco_turb(face_voisins(fac4,1)));
527 const double d_visco_lam = nu_mean_4_pts_(face_voisins(fac3,0), face_voisins(fac3,1), face_voisins(fac4,0), face_voisins(fac4,1));
528
529 if (ori1 == 1) // (seule possibilite : ori3 =0) Arete XY
530 {
531 // flux de mu_lam*tau21 + mu_turb*(tau21+tau12) sur la facette a cheval sur les faces fac1 et fac2
532 double flux1;
533 flux1 = (d_visco_lam + d_visco_turb)*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
534 fill_coeff_matrice_morse(fac3,fac4,flux1,matrice);
535
536 // Termes supplementaires dans le laplacien en axi : Ils sont integres comme des termes sources
537 const double coef_laplacien_axi = 0.5*(d_visco_lam + d_visco_turb);
538
539 for (auto l = tab1[fac1]-1; l < tab1[fac1+1]-1; l++)
540 if (tab2[l]-1 == fac1) coeff[l] += coef_laplacien_axi*volumes_entrelaces(fac1)*porosite(fac1)/xv(fac1,0);
541
542 for (auto l = tab1[fac2]-1; l < tab1[fac2+1]-1; l++)
543 if (tab2[l]-1 == fac2) coeff[l] += coef_laplacien_axi*volumes_entrelaces(fac2)*porosite(fac2)/xv(fac2,0);
544
545 // flux de mu_lam*tau12 + mu_turb*(tau21+tau12) sur la facette a cheval sur les faces fac3 et fac4
546 flux1 = (d_visco_lam + d_visco_turb)*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
547 fill_coeff_matrice_morse(fac1,fac2,flux1,matrice);
548 }
549 else // 2 possibilites : ori1 = 2 et ori3 == 1: arete YZ OU ori1 = 2 et ori3 = 0: arete XZ
550 {
551 double flux2;
552 // flux de mu_lam*tau31 + mu_turb*(tau31+tau13) sur la facette a cheval sur les faces fac1 et fac2
553 flux2 = (d_visco_lam+ d_visco_turb)*0.25*(surface(fac1)+surface(fac2))*(porosite(fac1)+porosite(fac2));
554 fill_coeff_matrice_morse(fac3,fac4,flux2,matrice);
555
556 // flux de mu_lam*tau13 + mu_turb*(tau13+tau31) sur la facette a cheval sur les faces fac3 et fac4
557 flux2 = (d_visco_lam + d_visco_turb)*0.25*(surface(fac3)+surface(fac4))*(porosite(fac3)+porosite(fac4));
558 fill_coeff_matrice_morse(fac1,fac2,flux2,matrice);
559 }
560 }
561}
562
563void Op_Dift_VDF_Face_Axi_base::ajouter_contribution(const DoubleTab& inco, Matrice_Morse& matrice ) const
564{
565 if (inco.line_size() > 1) not_implemented(__func__);
566 const double temps = equation().schema_temps().temps_courant();
567 mettre_a_jour_var(temps); // seulement pour var_axi !
568
569 const Domaine_Cl_VDF& zclvdf = la_zcl_vdf.valeur();
570 const DoubleVect& visco_turb = diffusivite_turbulente().valeurs();
571 const DoubleTab& tau_diag = inconnue->tau_diag();
572 ref_cast_non_const(Champ_Face_VDF,inconnue.valeur()).calculer_dercov_axi(zclvdf);
573
574 ajouter_contribution_elem(visco_turb,tau_diag,matrice); // Boucle sur les elements pour traiter les facettes situees a l'interieur des elements
575 if (dimension == 3) ajouter_contribution_elem_3D(visco_turb,tau_diag,matrice);
576 ajouter_contribution_aretes_bords(visco_turb,tau_diag,matrice); // Boucle sur les aretes bord
577 ajouter_contribution_aretes_mixtes(matrice); // Boucle sur les aretes mixtes
578 ajouter_contribution_aretes_internes(visco_turb,matrice); // Boucle sur les aretes internes
579}
580
582{
583 if (resu.line_size() > 1) not_implemented(__func__);
584 const double temps = equation().schema_temps().temps_courant();
585 mettre_a_jour_var(temps); // seulement pour var_axi !
586
587 const int ndeb = le_dom_vdf->premiere_arete_bord(), nfin = ndeb + le_dom_vdf->nb_aretes_bord();
588 for (int n_arete = ndeb; n_arete < nfin; n_arete++)
589 {
590 const int n_type = type_arete_bord(n_arete-ndeb);
591 switch(n_type)
592 {
594 {
595 const int fac1 = Qdm(n_arete,0), fac2 = Qdm(n_arete,1), fac3 = Qdm(n_arete,2), ori3 = orientation(fac3);
596 const int rang1 = (fac1 - le_dom_vdf->premiere_face_bord()), rang2 = (fac2 - le_dom_vdf->premiere_face_bord());
597 double coef;
598 if (is_var())
599 {
600 // Calcul du frottement identique a celui de TRIOVF : On calcule la moyenne des u_star et on l'eleve au carre. On calcule la moyenne des surfaces
601 const double tau = 0.5*(sqrt(tau_tan(rang1,ori3)) + sqrt(tau_tan(rang2,ori3))), surf = 0.5*(surface(fac1)+surface(fac2));
602 coef = tau*tau*surf;
603 }
604 else // Autre solution pour le calcul du frottement : On calcule u_star*u_star*surf sur chaque partie de la facette de Qdm
605 {
606 const double tau1 = tau_tan(rang1,ori3)*0.5*surface(fac1), tau2 = tau_tan(rang2,ori3)*0.5*surface(fac2);
607 coef = tau1+tau2;
608 }
609
610 resu[fac3] += coef;
611 break;
612 }
613 case TypeAreteBordVDF::FLUIDE_FLUIDE: /* fall through */
616 break;
617 default :
618 {
619 Cerr << "On a rencontre un type d'arete non prevu : [ num arete : " << n_arete << " ], [ type : " << n_type << " ]" << finl;
621 break;
622 }
623 }
624 }
625}
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
class Champ_Face_VDF Cette classe sert a representer un champ vectoriel dont on ne calcule
void dimensionner_tenseur_Grad()
Classe Champ_Inc_base.
class Domaine_Cl_VDF
int type_arete_bord(int num_arete) const
classe Domaine_Cl_dis_base Les objets Domaine_Cl_dis_base representent les conditions aux limites
class Domaine_VDF
Definition Domaine_VDF.h:64
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
int Qdm(int num_arete, int) const
virtual const DoubleVect & face_surfaces() const
Definition Domaine_VF.h:51
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
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
double xp(int num_elem, int k) const
Definition Domaine_VF.h:77
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
virtual const RefObjU & get_modele(Type_modele type) const
virtual const Champ_Inc_base & inconnue() const =0
Schema_Temps_base & schema_temps()
Renvoie le schema en temps associe a l'equation.
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
auto & get_set_tab2()
auto & get_set_coeff()
auto & get_set_tab1()
Classe Modele_turbulence_hyd_base Cette classe sert de base a la hierarchie des classes.
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
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
const Champ_Fonc_base & diffusivite_turbulente() const
void associer_modele_turbulence(const Modele_turbulence_hyd_base &)
virtual double nu_mean_4_pts_(const int, const int) const =0
void associer(const Domaine_dis_base &, const Domaine_Cl_dis_base &, const Champ_Inc_base &) override
void completer() override
Associe l'operateur au domaine_dis, le domaine_Cl_dis, et a l'inconnue de son equation.
virtual double nu_(const int) const =0
double calculer_dt_stab() const override
Calcul dt_stab.
void associer_diffusivite_turbulente(const Champ_Fonc_base &visc_turb)
void mettre_a_jour(double) override
DOES NOTHING - to override in derived classes.
void associer_loipar(const Turbulence_paroi_base &)
DoubleTab & calculer(const DoubleTab &, DoubleTab &) const override
calcule la contribution de la diffusion, la range dans resu
virtual double nu_mean_2_pts_(const int, const int) const =0
virtual bool is_var() const =0
void contribue_au_second_membre(DoubleTab &) const
virtual void mettre_a_jour_var(double) const =0
double calculer_dt_stab() const override
Calcul dt_stab.
virtual void completer()
Associe l'operateur au domaine_dis, le domaine_Cl_dis, et a l'inconnue de son equation.
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
double temps_courant() const
Renvoie le temps courant.
Classe de base des flux de sortie.
Definition Sortie.h:52
int line_size() const
Definition TRUSTVect.tpp:67
const Objet_U & valeur() const
Definition TRUST_Ref.h:134