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