1 : C YLAP2B SOURCE MAGN 09/07/30 21:22:29 6458
2 : SUBROUTINE YLAP2B(ICOGRT,MPROC,MPVITC,MPTEMC,
3 : $ MPVOLU,MPNORM,MPSURF,MELEFL,
4 : $ KRFACE,KRCENT,LCLIMT,KRTIMP,LCLIMQ,KRQIMP,
5 : $ NOMINC,
6 : $ KAPPA,CV,
7 : $ IJACO,
8 : $ IMPR,IRET)
9 : IMPLICIT INTEGER(I-N)
10 : IMPLICIT REAL*8 (A-H,O-Z)
11 : C***********************************************************************
12 : C NOM : YLAP2B
13 : C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
14 : C VF 3D.
15 : C Ici, on ne calcule que les contributions à la matrice
16 : C jacobienne faisant intervenir les coefficients pour le
17 : C calcul des gradients de température (ICOGRT).
18 : C (contributions à d Res_{\rho e_t} / d var
19 : C var prenant successivement les valeurs :
20 : C \rho, \rho u, \rho v, \rho w, \rho e_t)
21 : C
22 : C
23 : C
24 : C
25 : C LANGAGE : ESOPE
26 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
27 : C mél : gounand@semt2.smts.cea.fr
28 : C***********************************************************************
29 : C APPELES (UTIL) : AJMTK : ajoute un objet de type MATSIM (non
30 : C standard) à un objet de type MATRIK.
31 : C APPELE PAR : YLAP2A : Calcul de la matrice jacobienne du
32 : C résidu du laplacien VF 3D.
33 : C***********************************************************************
34 : C ENTREES : ICOGRT (type MCHELM) : coefficients pour le
35 : C calcul du gradient de la température aux
36 : C interfaces.
37 : C MPROC (type MPOVAL) : masse volumique par
38 : C élément.
39 : C MPVITC (type MPOVAL) : vitesse par élément.
40 : C MPTEMC (type MPOVAL) : température par élément.
41 : C MPVOLU (type MPOVAL) : volume des éléments.
42 : C MPNORM (type MPOVAL) : normale aux faces.
43 : C MPSURF (type MPOVAL) : surface des faces.
44 : C MELEFL (type MELEME) : connectivités face-(centre
45 : C gauche, centre droit).
46 : C KRFACE (type MLENTI) : tableau de repérage dans
47 : C le maillage des faces des éléments.
48 : C KRCENT (type MLENTI) : tableau de repérage dans
49 : C le maillage des centres des éléments.
50 : C LCLIMT (type logique) : .TRUE. => CL de Dirichlet
51 : C sur la température.
52 : C KRTIMP (type MLENTI) : tableau de repérage dans
53 : C maillage des CL de Dirichlet sur la
54 : C température.
55 : C LCLIMQ (type logique) : .TRUE. => CL de Dirichlet
56 : C sur le flux de chaleur.
57 : C KRTIMP (type MLENTI) : tableau de repérage dans
58 : C maillage des CL de Dirichlet sur la
59 : C température.
60 : C KRQIMP (type MLENTI) : tableau de repérage dans
61 : C maillage des CL de Dirichlet sur le flux de
62 : C chaleur.
63 : C NOMINC (type MLMOTS) : noms des inconnues.
64 : C KAPPA (type réel) : conductivité thermique (SI)
65 : C CV (type réel) : chaleur spécifique à volume
66 : C constant (SI).
67 : C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du
68 : C résidu du laplacien VF 3D.
69 : C SORTIES : -
70 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
71 : C***********************************************************************
72 : C VERSION : v1, 28/08/2001, version initiale
73 : C HISTORIQUE : v1, 28/08/2001, création
74 : C HISTORIQUE :
75 : C HISTORIQUE :
76 : C***********************************************************************
77 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
78 : C en cas de modification de ce sous-programme afin de faciliter
79 : C la maintenance !
80 : C***********************************************************************
81 : -INC CCOPTIO
82 : -INC SMCHPOI
83 : POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL,MPTEMC.MPOVAL
84 : POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
85 : -INC SMCHAML
86 : POINTEUR ICOGRT.MCHELM,JCOGRT.MCHAML
87 : POINTEUR KCDTDX.MELVAL,KCDTDY.MELVAL,KCDTDZ.MELVAL
88 : -INC SMELEME
89 : POINTEUR MELEFL.MELEME
90 : POINTEUR MCOGRT.MELEME
91 : -INC SMLENTI
92 : POINTEUR KRTIMP.MLENTI,KRQIMP.MLENTI
93 : POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
94 : -INC SMLMOTS
95 : POINTEUR NOMINC.MLMOTS
96 : -INC SMMATRIK
97 : POINTEUR IJACO.MATRIK
98 : *
99 : * Objet matrice élémentaire simplifié
100 : *
101 : SEGMENT GMATSI
102 : INTEGER POIPR1(NPP1,NEL1)
103 : INTEGER POIDU1(1,NEL1)
104 : INTEGER POIPR2(NPP2,NEL2)
105 : INTEGER POIDU2(2,NEL2)
106 : POINTEUR LMATSI(0).MATSIM
107 : ENDSEGMENT
108 : * Contributions de la part du gradient de température (CTGRT)
109 : POINTEUR CTGRT.GMATSI
110 : SEGMENT MATSIM
111 : CHARACTER*8 NOMPRI,NOMDUA
112 : REAL*8 VALMA1(1,NPP1,NEL1)
113 : REAL*8 VALMA2(2,NPP2,NEL2)
114 : ENDSEGMENT
115 : POINTEUR RETRHO.MATSIM
116 : POINTEUR RETROU.MATSIM
117 : POINTEUR RETROV.MATSIM
118 : POINTEUR RETROW.MATSIM
119 : POINTEUR RETRET.MATSIM
120 : *
121 : REAL*8 KAPPA,CV
122 : *
123 : INTEGER IMPR,IRET
124 : *
125 : LOGICAL LCLIMT,LCLIMQ
126 : LOGICAL LMUR
127 : LOGICAL LCTRB1
128 : *
129 : INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
130 : INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
131 : INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
132 : INTEGER NLCENP,NLCEND,NLFACE,NLCLQ,NLCLT
133 : INTEGER NPTEL
134 : *
135 : REAL*8 BETAX,BETAY,BETAZ,CNX,CNY,CNZ
136 : REAL*8 SIGNOR,SURFFA,VOLUEL
137 : REAL*8 RHOP,UP,VP,WP,TP
138 : REAL*8 FACTOR
139 : REAL*8 DTDRHO,DTDROU,DTDROV,DTDROW,DTDRET
140 : C
141 : *
142 : * Executable statements
143 : *
144 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap2b.eso'
145 : * On calcule la partie de d Res_{\rho e_t} / d var
146 : * (var prend successivement les valeurs :
147 : * \rho, \rho u, \rho v, \rho w, \rho e_t)
148 : * faisant intervenir les coefficients pour le calcul des gradients de
149 : * température (ICOGRT).
150 : * C'est la partie contenant le terme : - \vect{q} \pscal \vect{n}
151 : * Les noms de matrices élémentaires (type MATSIM) associées sont :
152 : * RETRHO, RETROU, RETROV, RETROW, RETRET
153 : IF (LCLIMT) THEN
154 : SEGACT KRTIMP
155 : ENDIF
156 : IF (LCLIMQ) THEN
157 : SEGACT KRQIMP
158 : ENDIF
159 : SEGACT NOMINC
160 : SEGACT KRCENT
161 : SEGACT KRFACE
162 : SEGACT MELEFL
163 : SEGACT MPSURF
164 : SEGACT MPNORM
165 : SEGACT MPVOLU
166 : SEGACT MPROC
167 : SEGACT MPVITC
168 : SEGACT MPTEMC
169 : SEGACT ICOGRT
170 : NSOUCH=ICOGRT.IMACHE(/1)
171 : DO 1 ISOUCH=1,NSOUCH
172 : MCOGRT=ICOGRT.IMACHE(ISOUCH)
173 : JCOGRT=ICOGRT.ICHAML(ISOUCH)
174 : SEGACT JCOGRT
175 : KCDTDX=JCOGRT.IELVAL(1)
176 : KCDTDY=JCOGRT.IELVAL(2)
177 : KCDTDZ=JCOGRT.IELVAL(3)
178 : SEGDES JCOGRT
179 : SEGACT KCDTDX
180 : SEGACT KCDTDY
181 : SEGACT KCDTDZ
182 : SEGACT MCOGRT
183 : NELEM=MCOGRT.NUM(/2)
184 : NPTEL=MCOGRT.NUM(/1)
185 : NPP1=NPTEL-1
186 : NPP2=NPTEL-1
187 : NEL1=NELEM
188 : NEL2=NELEM
189 : IEL1=1
190 : IEL2=1
191 : SEGINI RETRHO
192 : SEGINI RETROU
193 : SEGINI RETROV
194 : SEGINI RETROW
195 : SEGINI RETRET
196 : SEGINI CTGRT
197 : RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
198 : RETRHO.NOMPRI(5:8)=' '
199 : RETRHO.NOMDUA(1:4)=NOMINC.MOTS(5)
200 : RETRHO.NOMDUA(5:8)=' '
201 : RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
202 : RETROU.NOMPRI(5:8)=' '
203 : RETROU.NOMDUA(1:4)=NOMINC.MOTS(5)
204 : RETROU.NOMDUA(5:8)=' '
205 : RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
206 : RETROV.NOMPRI(5:8)=' '
207 : RETROV.NOMDUA(1:4)=NOMINC.MOTS(5)
208 : RETROV.NOMDUA(5:8)=' '
209 : RETROW.NOMPRI(1:4)=NOMINC.MOTS(4)
210 : RETROW.NOMPRI(5:8)=' '
211 : RETROW.NOMDUA(1:4)=NOMINC.MOTS(5)
212 : RETROW.NOMDUA(5:8)=' '
213 : RETRET.NOMPRI(1:4)=NOMINC.MOTS(5)
214 : RETRET.NOMPRI(5:8)=' '
215 : RETRET.NOMDUA(1:4)=NOMINC.MOTS(5)
216 : RETRET.NOMDUA(5:8)=' '
217 : DO 12 IELEM=1,NELEM
218 : * Le premier point du support de ICOGRT est un point FACE
219 : NGFACE=MCOGRT.NUM(1,IELEM)
220 : NLFACE=KRFACE.LECT(NGFACE)
221 : IF (NLFACE.EQ.0) THEN
222 : WRITE(IOIMP,*) 'Erreur de programmation n°1'
223 : GOTO 9999
224 : ENDIF
225 : * On calcule la contribution à la matrice jacobienne IJACO de la face
226 : * NGFAC (points duaux : centres à gauche et à droite de la face)
227 : * (points primaux : une partie (bicoz conditions aux limites)
228 : * de ceux du stencil pour le calcul du gradient
229 : * à la face, ils doivent être des points centres)
230 : * Si le flux de chaleur sur la face est imposé par les conditions
231 : * aux limites, la contribution de la face à IJACO est nulle.
232 : LCTRB1=.TRUE.
233 : IF (LCLIMQ) THEN
234 : NLCLQ=KRQIMP.LECT(NGFACE)
235 : IF (NLCLQ.NE.0) THEN
236 : LCTRB1=.FALSE.
237 : ENDIF
238 : ENDIF
239 : IF (LCTRB1) THEN
240 : NGCGAU=MELEFL.NUM(1,NLFACE)
241 : NGCDRO=MELEFL.NUM(3,NLFACE)
242 : LMUR=(NGCGAU.EQ.NGCDRO)
243 : * On distingue le cas où la face est un bord du maillage (mur)
244 : * du cas où la face est interne au maillage
245 : IF (.NOT.LMUR) THEN
246 : NPD=2
247 : ELSE
248 : NPD=1
249 : ENDIF
250 : NPP=NPTEL-1
251 : * IPD=1 : point à gauche du point NGFACE
252 : * IPD=2 : point à droite du point NGFACE
253 : DO 122 IPD=1,NPD
254 : NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
255 : IF (.NOT.LMUR) THEN
256 : CTGRT.POIDU2(IPD,IEL2)=NPDUAL
257 : ELSE
258 : CTGRT.POIDU1(IPD,IEL1)=NPDUAL
259 : ENDIF
260 : NLCEND=KRCENT.LECT(NPDUAL)
261 : IF (NLCEND.EQ.0) THEN
262 : WRITE(IOIMP,*) 'Erreur grave n°1'
263 : GOTO 9999
264 : ENDIF
265 : DO 124 IPP=1,NPP
266 : NPPRIM=MCOGRT.NUM(IPP+1,IELEM)
267 : C
268 : C AB : we do not check the BC anymore
269 : C
270 : NLCENP=KRCENT.LECT(NPPRIM)
271 : IF (NLCENP .EQ. 0) THEN
272 : * Lorsque une contribution est nulle, on fixe artificiellement le
273 : * point primal égal au point dual.
274 : IF (.NOT.LMUR) THEN
275 : CTGRT.POIPR2(IPP,IEL2)=NPDUAL
276 : RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
277 : RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
278 : RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
279 : RETROW.VALMA2(IPD,IPP,IEL2)=0.D0
280 : RETRET.VALMA2(IPD,IPP,IEL2)=0.D0
281 : ELSE
282 : CTGRT.POIPR1(IPP,IEL1)=NPDUAL
283 : RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
284 : RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
285 : RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
286 : RETROW.VALMA1(IPD,IPP,IEL1)=0.D0
287 : RETRET.VALMA1(IPD,IPP,IEL1)=0.D0
288 : ENDIF
289 : ELSE
290 : * Les contributions valent :
291 : * (d Res_{\rho e_t})_d / (d var)_p =
292 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \kappa
293 : * * [ ((n_x * \beta_x) + (n_y * \beta_y) + (n_z * \beta_z)) *
294 : * ((dT)_p / (d var)_p)]
295 : * avec :
296 : * (dT)_p / (d \rho)_p = (1 / \rho_p) *
297 : * ( (((u_p)^2 + (v_p)^2 + (w_p)^2) / (2 * c_v)) - T )
298 : * (dT)_p / (d \rho u)_p = - u_p / (\rho_p * c_v)
299 : * (dT)_p / (d \rho v)_p = - v_p / (\rho_p * c_v)
300 : * (dT)_p / (d \rho w)_p = - w_p / (\rho_p * c_v)
301 : * (dT)_p / (d \rho e_t)_p = 1 / (\rho_p * c_v)
302 : * \beta_x : coefficients pour le calcul de dT/dx
303 : * \beta_y : coefficients pour le calcul de dT/dy
304 : * \beta_z : coefficients pour le calcul de dT/dz
305 : *
306 : * normale sortante pour IPD=1, rentrante pour IPD=2
307 : SIGNOR=(-1.D0)**(IPD+1)
308 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
309 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
310 : CNX =MPNORM.VPOCHA(NLFACE,1)
311 : CNY =MPNORM.VPOCHA(NLFACE,2)
312 : CNZ =MPNORM.VPOCHA(NLFACE,3)
313 : BETAX =KCDTDX.VELCHE(IPP+1,IELEM)
314 : BETAY =KCDTDY.VELCHE(IPP+1,IELEM)
315 : BETAZ =KCDTDZ.VELCHE(IPP+1,IELEM)
316 : RHOP =MPROC.VPOCHA(NLCENP,1)
317 : UP =MPVITC.VPOCHA(NLCENP,1)
318 : VP =MPVITC.VPOCHA(NLCENP,2)
319 : WP =MPVITC.VPOCHA(NLCENP,3)
320 : TP =MPTEMC.VPOCHA(NLCENP,1)
321 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*KAPPA
322 : $ *((CNX*BETAX)+(CNY*BETAY)+(CNZ*BETAZ))
323 : DTDRHO=((((UP*UP)+(VP*VP)+(WP*WP))
324 : $ /(2.D0*CV))-TP)/RHOP
325 : DTDROU=-UP/(RHOP*CV)
326 : DTDROV=-VP/(RHOP*CV)
327 : DTDROW=-WP/(RHOP*CV)
328 : DTDRET=1.D0/(RHOP*CV)
329 : IF (.NOT.LMUR) THEN
330 : CTGRT.POIPR2(IPP,IEL2)=NPPRIM
331 : RETRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRHO
332 : RETROU.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROU
333 : RETROV.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROV
334 : RETROW.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROW
335 : RETRET.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRET
336 : ELSE
337 : CTGRT.POIPR1(IPP,IEL1)=NPPRIM
338 : RETRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRHO
339 : RETROU.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROU
340 : RETROV.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROV
341 : RETROW.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROW
342 : RETRET.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRET
343 : ENDIF
344 : ENDIF
345 : 124 CONTINUE
346 : 122 CONTINUE
347 : IF (.NOT.LMUR) THEN
348 : IEL2=IEL2+1
349 : ELSE
350 : IEL1=IEL1+1
351 : ENDIF
352 : ENDIF
353 : 12 CONTINUE
354 : NPP1=NPTEL-1
355 : NPP2=NPTEL-1
356 : NEL1=IEL1-1
357 : NEL2=IEL2-1
358 : SEGADJ RETRHO
359 : SEGADJ RETROU
360 : SEGADJ RETROV
361 : SEGADJ RETROW
362 : SEGADJ RETRET
363 : SEGADJ CTGRT
364 : CTGRT.LMATSI(**)=RETRHO
365 : CTGRT.LMATSI(**)=RETROU
366 : CTGRT.LMATSI(**)=RETROV
367 : CTGRT.LMATSI(**)=RETROW
368 : CTGRT.LMATSI(**)=RETRET
369 : * On accumule les matrices résultantes dans IJACO
369 : CALL AJMTK(CTGRT,IJACO,IMPR,IRET)
370 : IF (IRET.NE.0) GOTO 9999
371 : SEGSUP RETRHO
372 : SEGSUP RETROU
373 : SEGSUP RETROV
374 : SEGSUP RETROW
375 : SEGSUP RETRET
376 : SEGSUP CTGRT
377 : *
378 : SEGDES MCOGRT
379 : SEGDES KCDTDZ
380 : SEGDES KCDTDY
381 : SEGDES KCDTDX
382 : 1 CONTINUE
383 : SEGDES ICOGRT
384 : SEGDES MPTEMC
385 : SEGDES MPVITC
386 : SEGDES MPROC
387 : SEGDES MPVOLU
388 : SEGDES MPNORM
389 : SEGDES MPSURF
390 : SEGDES MELEFL
391 : SEGDES KRFACE
392 : SEGDES KRCENT
393 : SEGDES NOMINC
394 : IF (LCLIMQ) THEN
395 : SEGDES KRQIMP
396 : ENDIF
397 : IF (LCLIMT) THEN
398 : SEGDES KRTIMP
399 : ENDIF
400 : *
401 : * Normal termination
402 : *
403 : IRET=0
404 : RETURN
405 : *
406 : * Format handling
407 : *
408 : *
409 : * Error handling
410 : *
411 : 9999 CONTINUE
412 : IRET=1
413 : WRITE(IOIMP,*) 'An error was detected in subroutine ylap2b'
414 : RETURN
415 : *
416 : * End of subroutine YLAP2B
417 : *
418 : END
419 :
420 :
421 :
422 :
423 :
424 :
425 :
426 :
© Cast3M 2003 - All rights reserved.
Disclaimer