1 : C YLAP1B SOURCE MAGN 09/07/30 21:22:15 6458
2 : SUBROUTINE YLAP1B(ICOGRT,MPROC,MPVITC,MPTEMC,
3 : $ MPVOLU,MPNORM,MPSURF,MELEFL,
4 : $ KRFACE,KRCENT,LCLIMT,KRTIMP,LCLIMQ,KRQIMP,
5 : $ LMIXT,KRMIXT,NOMINC,ICLAU,
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 : YLAP1B
13 : C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
14 : C VF 2D.
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 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 : YLAP1A : Calcul de la matrice jacobienne du
32 : C résidu du laplacien VF 2D.
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 LMIXT (type logique) : .TRUE. => CL mixtes
64 : C KRMIXT (type MLENTI) : tableau de repérage dans
65 : C maillage des CL mixtes.
66 : C NOMINC (type MLMOTS) : noms des inconnues.
67 : C ICLAU : option pour ne calculer
68 : C que la thermique
69 : C KAPPA (type réel) : conductivité thermique (SI)
70 : C CV (type réel) : chaleur spécifique à volume
71 : C constant (SI).
72 : C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du
73 : C résidu du laplacien VF 2D.
74 : C SORTIES : -
75 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
76 : C***********************************************************************
77 : C VERSION : v1, 08/08/2001, version initiale
78 : C HISTORIQUE : v1, 08/08/2001, création
79 : C HISTORIQUE : 11/02/2003, ajout de l'option 'MIXT'
80 : C HISTORIQUE :
81 : C***********************************************************************
82 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
83 : C en cas de modification de ce sous-programme afin de faciliter
84 : C la maintenance !
85 : C***********************************************************************
86 : -INC CCOPTIO
87 : -INC SMCHPOI
88 : POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL,MPTEMC.MPOVAL
89 : POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
90 : -INC SMCHAML
91 : POINTEUR ICOGRT.MCHELM,JCOGRT.MCHAML
92 : POINTEUR KCDTDX.MELVAL,KCDTDY.MELVAL
93 : -INC SMELEME
94 : POINTEUR MELEFL.MELEME
95 : POINTEUR MCOGRT.MELEME
96 : -INC SMLENTI
97 : POINTEUR KRTIMP.MLENTI,KRQIMP.MLENTI,KRMIXT.MLENTI
98 : POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
99 : -INC SMLMOTS
100 : POINTEUR NOMINC.MLMOTS
101 : -INC SMMATRIK
102 : POINTEUR IJACO.MATRIK
103 : *
104 : * Objet matrice élémentaire simplifié
105 : *
106 : SEGMENT GMATSI
107 : INTEGER POIPR1(NPP1,NEL1)
108 : INTEGER POIDU1(1,NEL1)
109 : INTEGER POIPR2(NPP2,NEL2)
110 : INTEGER POIDU2(2,NEL2)
111 : POINTEUR LMATSI(0).MATSIM
112 : ENDSEGMENT
113 : * Contributions de la part du gradient de température (CTGRT)
114 : POINTEUR CTGRT.GMATSI
115 : SEGMENT MATSIM
116 : CHARACTER*8 NOMPRI,NOMDUA
117 : REAL*8 VALMA1(1,NPP1,NEL1)
118 : REAL*8 VALMA2(2,NPP2,NEL2)
119 : ENDSEGMENT
120 : POINTEUR RETRHO.MATSIM
121 : POINTEUR RETROU.MATSIM
122 : POINTEUR RETROV.MATSIM
123 : POINTEUR RETRET.MATSIM
124 : *
125 : REAL*8 KAPPA,CV
126 : *
127 : INTEGER IMPR,IRET,ICLAU
128 : *
129 : LOGICAL LCLIMT,LCLIMQ,LMIXT
130 : LOGICAL LMUR
131 : LOGICAL LCTRB1,LCTRB2
132 : *
133 : INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
134 : INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
135 : INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
136 : INTEGER NLCENP,NLCEND,NLFACE,NLCLQ,NLCLT
137 : INTEGER NPTEL
138 : *
139 : REAL*8 BETAX,BETAY,CNX,CNY
140 : REAL*8 SIGNOR,SURFFA,VOLUEL
141 : REAL*8 RHOP,UP,VP,TP
142 : REAL*8 FACTOR
143 : REAL*8 DTDRHO,DTDROU,DTDROV,DTDRET
144 : C
145 : *
146 : * Executable statements
147 : *
148 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap1b.eso'
149 : * On calcule la partie de d Res_{\rho e_t} / d var
150 : * (var prend successivement les valeurs :
151 : * \rho, \rho u, \rho v, \rho e_t)
152 : * faisant intervenir les coefficients pour le calcul des gradients de
153 : * température (ICOGRT).
154 : * C'est la partie contenant le terme : - \vect{q} \pscal \vect{n}
155 : * Les noms de matrices élémentaires (type MATSIM) associées sont :
156 : * RETRHO, RETROU, RETROV, RETRET
157 : IF (LCLIMT) THEN
158 : SEGACT KRTIMP
159 : ENDIF
160 : IF (LCLIMQ) THEN
161 : SEGACT KRQIMP
162 : ENDIF
163 : IF (LMIXT) THEN
164 : SEGACT KRMIXT
165 : ENDIF
166 : SEGACT NOMINC
167 : SEGACT KRCENT
168 : SEGACT KRFACE
169 : SEGACT MELEFL
170 : SEGACT MPSURF
171 : SEGACT MPNORM
172 : SEGACT MPVOLU
173 : SEGACT MPROC
174 : SEGACT MPVITC
175 : SEGACT MPTEMC
176 : SEGACT ICOGRT
177 : NSOUCH=ICOGRT.IMACHE(/1)
178 : DO 1 ISOUCH=1,NSOUCH
179 : MCOGRT=ICOGRT.IMACHE(ISOUCH)
180 : JCOGRT=ICOGRT.ICHAML(ISOUCH)
181 : SEGACT JCOGRT
182 : KCDTDX=JCOGRT.IELVAL(1)
183 : KCDTDY=JCOGRT.IELVAL(2)
184 : SEGDES JCOGRT
185 : SEGACT KCDTDX
186 : SEGACT KCDTDY
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 RETRET
200 : SEGINI CTGRT
201 : RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
202 : RETRHO.NOMPRI(5:8)=' '
203 : RETRHO.NOMDUA(1:4)=NOMINC.MOTS(4)
204 : RETRHO.NOMDUA(5:8)=' '
205 : RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
206 : RETROU.NOMPRI(5:8)=' '
207 : RETROU.NOMDUA(1:4)=NOMINC.MOTS(4)
208 : RETROU.NOMDUA(5:8)=' '
209 : RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
210 : RETROV.NOMPRI(5:8)=' '
211 : RETROV.NOMDUA(1:4)=NOMINC.MOTS(4)
212 : RETROV.NOMDUA(5:8)=' '
213 : RETRET.NOMPRI(1:4)=NOMINC.MOTS(4)
214 : RETRET.NOMPRI(5:8)=' '
215 : RETRET.NOMDUA(1:4)=NOMINC.MOTS(4)
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 : LCTRB2=.TRUE.
268 : IF (LCLIMT) THEN
269 : NLCLT=KRTIMP.LECT(NPPRIM)
270 : IF (NLCLT.NE.0) THEN
271 : LCTRB2=.FALSE.
272 : ENDIF
273 : ENDIF
274 : IF (LCLIMQ) THEN
275 : NLCLQ=KRQIMP.LECT(NPPRIM)
276 : IF (NLCLQ.NE.0) THEN
277 : LCTRB2=.FALSE.
278 : ENDIF
279 : ENDIF
280 : IF (LMIXT) THEN
281 : NLCLM=KRMIXT.LECT(NPPRIM)
282 : IF (NLCLM.NE.0) THEN
283 : LCTRB2=.FALSE.
284 : ENDIF
285 : ENDIF
286 : IF (.NOT.LCTRB2) THEN
287 : * Lorsque une contribution est nulle, on fixe artificiellement le
288 : * point primal égal au point dual.
289 : IF (.NOT.LMUR) THEN
290 : CTGRT.POIPR2(IPP,IEL2)=NPDUAL
291 : RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
292 : RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
293 : RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
294 : RETRET.VALMA2(IPD,IPP,IEL2)=0.D0
295 : ELSE
296 : CTGRT.POIPR1(IPP,IEL1)=NPDUAL
297 : RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
298 : RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
299 : RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
300 : RETRET.VALMA1(IPD,IPP,IEL1)=0.D0
301 : ENDIF
302 : ELSE
303 : * Les contributions valent :
304 : * (d Res_{\rho e_t})_d / (d var)_p =
305 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \kappa
306 : * * [ ((n_x * \beta_x) +(n_y * \beta_y)) *
307 : * ((dT)_p / (d var)_p)]
308 : * avec :
309 : * (dT)_p / (d \rho)_p = (1 / \rho_p) *
310 : * ( (((u_p)^2 + (v_p)^2) / (2 * c_v)) - T )
311 : * (dT)_p / (d \rho u)_p = - u_p / (\rho_p * c_v)
312 : * (dT)_p / (d \rho v)_p = - v_p / (\rho_p * c_v)
313 : * (dT)_p / (d \rho e_t)_p = 1 / (\rho_p * c_v)
314 : * \beta_x : coefficients pour le calcul de dT/dx
315 : * \beta_y : coefficients pour le calcul de dT/dy
316 : *
317 : NLCENP=KRCENT.LECT(NPPRIM)
318 : IF (NLCENP.EQ.0) THEN
319 : WRITE(IOIMP,*) 'Erreur grave n°2'
320 : GOTO 9999
321 : ENDIF
322 : * normale sortante pour IPD=1, rentrante pour IPD=2
323 : SIGNOR=(-1.D0)**(IPD+1)
324 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
325 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
326 : CNX =MPNORM.VPOCHA(NLFACE,1)
327 : CNY =MPNORM.VPOCHA(NLFACE,2)
328 : BETAX =KCDTDX.VELCHE(IPP+1,IELEM)
329 : BETAY =KCDTDY.VELCHE(IPP+1,IELEM)
330 : RHOP =MPROC.VPOCHA(NLCENP,1)
331 : UP =MPVITC.VPOCHA(NLCENP,1)
332 : VP =MPVITC.VPOCHA(NLCENP,2)
333 : TP =MPTEMC.VPOCHA(NLCENP,1)
334 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*KAPPA
335 : $ *((CNX*BETAX)+(CNY*BETAY))
336 : DTDRHO=((((UP*UP)+(VP*VP))/(2.D0*CV))-TP)/RHOP
337 : DTDROU=-UP/(RHOP*CV)
338 : DTDROV=-VP/(RHOP*CV)
339 : DTDRET=1.D0/(RHOP*CV)
340 : IF (.NOT.LMUR) THEN
341 : CTGRT.POIPR2(IPP,IEL2)=NPPRIM
342 : RETRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRHO
343 : RETROU.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROU
344 : RETROV.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROV
345 : RETRET.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRET
346 : ELSE
347 : CTGRT.POIPR1(IPP,IEL1)=NPPRIM
348 : RETRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRHO
349 : RETROU.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROU
350 : RETROV.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROV
351 : RETRET.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRET
352 : ENDIF
353 : ENDIF
354 : 124 CONTINUE
355 : 122 CONTINUE
356 : IF (.NOT.LMUR) THEN
357 : IEL2=IEL2+1
358 : ELSE
359 : IEL1=IEL1+1
360 : ENDIF
361 : ENDIF
362 : 12 CONTINUE
363 : NPP1=NPTEL-1
364 : NPP2=NPTEL-1
365 : NEL1=IEL1-1
366 : NEL2=IEL2-1
367 : SEGADJ RETRHO
368 : SEGADJ RETROU
369 : SEGADJ RETROV
370 : SEGADJ RETRET
371 : SEGADJ CTGRT
372 : IF (ICLAU.EQ.0) THEN
373 : CTGRT.LMATSI(**)=RETRHO
374 : CTGRT.LMATSI(**)=RETROU
375 : CTGRT.LMATSI(**)=RETROV
376 : ENDIF
377 : CTGRT.LMATSI(**)=RETRET
378 : * On accumule les matrices résultantes dans IJACO
378 : CALL AJMTK(CTGRT,IJACO,IMPR,IRET)
379 : IF (IRET.NE.0) GOTO 9999
380 : SEGSUP RETRHO
381 : SEGSUP RETROU
382 : SEGSUP RETROV
383 : SEGSUP RETRET
384 : SEGSUP CTGRT
385 : *
386 : SEGDES MCOGRT
387 : SEGDES KCDTDY
388 : SEGDES KCDTDX
389 : 1 CONTINUE
390 : SEGDES ICOGRT
391 : SEGDES MPTEMC
392 : SEGDES MPVITC
393 : SEGDES MPROC
394 : SEGDES MPVOLU
395 : SEGDES MPNORM
396 : SEGDES MPSURF
397 : SEGDES MELEFL
398 : SEGDES KRFACE
399 : SEGDES KRCENT
400 : SEGDES NOMINC
401 : IF (LCLIMQ) THEN
402 : SEGDES KRQIMP
403 : ENDIF
404 : IF (LCLIMT) THEN
405 : SEGDES KRTIMP
406 : ENDIF
407 : IF (LMIXT) THEN
408 : SEGDES KRMIXT
409 : ENDIF
410 : *
411 : * Normal termination
412 : *
413 : IRET=0
414 : RETURN
415 : *
416 : * Format handling
417 : *
418 : *
419 : * Error handling
420 : *
421 : 9999 CONTINUE
422 : IRET=1
423 : WRITE(IOIMP,*) 'An error was detected in subroutine ylap1b'
424 : RETURN
425 : *
426 : * End of subroutine YLAP1B
427 : *
428 : END
429 :
430 :
431 :
432 :
433 :
434 :
435 :
436 :
© Cast3M 2003 - All rights reserved.
Disclaimer