1 : C YLA1BT SOURCE GOUNAND 10/10/13 21:16:06 6779
2 : C YLAP1B SOURCE LEPOTIER 02/12/17 21:19:20 4510
3 : SUBROUTINE YLA1BT(ICOGRT,MPTEMC,
4 : $ MPVOLU,MPNORM,MPSURF,MELEFL,
5 : $ KRFACE,KRCENT,LCLIMT,KRTIMP,LCLIMQ,KRQIMP,
6 : $ LMIXT,KRMIXT,IJACO,
7 : $ IMPR,IRET)
8 : IMPLICIT INTEGER(I-N)
9 : IMPLICIT REAL*8 (A-H,O-Z)
10 : C***********************************************************************
11 : C NOM : YLA1BT
12 : C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
13 : C VF 2D.
14 : C Ici, on ne calcule que les contributions à la matrice
15 : C jacobienne faisant intervenir les coefficients pour le
16 : C calcul des gradients de température (ICOGRT).
17 : C (contributions à d Res_{\rho e_t} / d var
18 : C var prenant successivement les valeurs :
19 : C \rho, \rho u, \rho v, \rho e_t)
20 : C
21 : C
22 : C
23 : C
24 : C LANGAGE : ESOPE
25 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
26 : C mél : gounand@semt2.smts.cea.fr
27 : C***********************************************************************
28 : C APPELES (UTIL) : AJMTK : ajoute un objet de type MATSIM (non
29 : C standard) à un objet de type MATRIK.
30 : C APPELE PAR : YLAP1A : Calcul de la matrice jacobienne du
31 : C résidu du laplacien VF 2D.
32 : C***********************************************************************
33 : C ENTREES : ICOGRT (type MCHELM) : coefficients pour le
34 : C calcul du gradient de la température aux
35 : C interfaces.
36 : C MPROC (type MPOVAL) : masse volumique par
37 : C élément.
38 : C MPVITC (type MPOVAL) : vitesse par élément.
39 : C MPTEMC (type MPOVAL) : température par élément.
40 : C MPVOLU (type MPOVAL) : volume des éléments.
41 : C MPNORM (type MPOVAL) : normale aux faces.
42 : C MPSURF (type MPOVAL) : surface des faces.
43 : C MELEFL (type MELEME) : connectivités face-(centre
44 : C gauche, centre droit).
45 : C KRFACE (type MLENTI) : tableau de repérage dans
46 : C le maillage des faces des éléments.
47 : C KRCENT (type MLENTI) : tableau de repérage dans
48 : C le maillage des centres des éléments.
49 : C LCLIMT (type logique) : .TRUE. => CL de Dirichlet
50 : C sur la température.
51 : C KRTIMP (type MLENTI) : tableau de repérage dans
52 : C maillage des CL de Dirichlet sur la
53 : C température.
54 : C LCLIMQ (type logique) : .TRUE. => CL de Dirichlet
55 : C sur le flux de chaleur.
56 : C KRTIMP (type MLENTI) : tableau de repérage dans
57 : C maillage des CL de Dirichlet sur la
58 : C température.
59 : C KRQIMP (type MLENTI) : tableau de repérage dans
60 : C maillage des CL de Dirichlet sur le flux de
61 : C chaleur.
62 : C LMIXT (type logique) : .TRUE. => CL mixtes
63 : C KRMIXT (type MLENTI) : tableau de repérage dans
64 : C maillage des CL mixtes.
65 : C NOMINC (type MLMOTS) : noms des inconnues.
66 : C ICLAU : option pour ne calculer
67 : C que la thermique
68 : C KAPPA (type réel) : conductivité thermique (SI)
69 : C CV (type réel) : chaleur spécifique à volume
70 : C constant (SI).
71 : C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du
72 : C résidu du laplacien VF 2D.
73 : C SORTIES : -
74 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
75 : C***********************************************************************
76 : C VERSION : v1, 08/08/2001, version initiale
77 : C HISTORIQUE : v1, 08/08/2001, création
78 : C HISTORIQUE : 11/02/2003, ajout de l'option 'MIXT'
79 : C HISTORIQUE :
80 : C***********************************************************************
81 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
82 : C en cas de modification de ce sous-programme afin de faciliter
83 : C la maintenance !
84 : C***********************************************************************
85 : -INC CCOPTIO
86 : -INC SMCHPOI
87 : POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL,MPTEMC.MPOVAL
88 : POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
89 : -INC SMCHAML
90 : POINTEUR ICOGRT.MCHELM,JCOGRT.MCHAML
91 : POINTEUR KCDTDX.MELVAL,KCDTDY.MELVAL
92 : -INC SMELEME
93 : POINTEUR MELEFL.MELEME
94 : POINTEUR MCOGRT.MELEME
95 : -INC SMLENTI
96 : POINTEUR KRTIMP.MLENTI,KRQIMP.MLENTI,KRMIXT.MLENTI
97 : POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
98 : -INC SMLMOTS
99 : POINTEUR NOMINC.MLMOTS
100 : -INC SMMATRIK
101 : POINTEUR IJACO.MATRIK
102 : *
103 : * Objet matrice élémentaire simplifié
104 : *
105 : SEGMENT GMATSI
106 : INTEGER POIPR1(NPP1,NEL1)
107 : INTEGER POIDU1(1,NEL1)
108 : INTEGER POIPR2(NPP2,NEL2)
109 : INTEGER POIDU2(2,NEL2)
110 : POINTEUR LMATSI(0).MATSIM
111 : ENDSEGMENT
112 : * Contributions de la part du gradient de température (CTGRT)
113 : POINTEUR CTGRT.GMATSI
114 : SEGMENT MATSIM
115 : CHARACTER*8 NOMPRI,NOMDUA
116 : REAL*8 VALMA1(1,NPP1,NEL1)
117 : REAL*8 VALMA2(2,NPP2,NEL2)
118 : ENDSEGMENT
119 : POINTEUR RETRET.MATSIM
120 : *
121 : REAL*8 KAPPA,CV
122 : *
123 : INTEGER IMPR,IRET,ICLAU
124 : *
125 : LOGICAL LCLIMT,LCLIMQ,LMIXT
126 : LOGICAL LMUR
127 : LOGICAL LCTRB1,LCTRB2
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,CNX,CNY
136 : REAL*8 SIGNOR,SURFFA,VOLUEL
137 : REAL*8 RHOP,UP,VP,TP
138 : REAL*8 FACTOR
139 : REAL*8 DTDRHO,DTDROU,DTDROV,DTDRET
140 : C
141 : *
142 : * Executable statements
143 : *
144 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap1b.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 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, RETRET
153 : IF (LCLIMT) THEN
154 : SEGACT KRTIMP
155 : ENDIF
156 : IF (LCLIMQ) THEN
157 : SEGACT KRQIMP
158 : ENDIF
159 : IF (LMIXT) THEN
160 : SEGACT KRMIXT
161 : ENDIF
162 : SEGACT KRCENT
163 : SEGACT KRFACE
164 : SEGACT MELEFL
165 : SEGACT MPSURF
166 : SEGACT MPNORM
167 : SEGACT MPVOLU
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 : SEGDES JCOGRT
177 : SEGACT KCDTDX
178 : SEGACT MCOGRT
179 : NELEM=MCOGRT.NUM(/2)
180 : NPTEL=MCOGRT.NUM(/1)
181 : NPP1=NPTEL-1
182 : NPP2=NPTEL-1
183 : NEL1=NELEM
184 : NEL2=NELEM
185 : IEL1=1
186 : IEL2=1
187 :
188 :
189 : c BOUCLE POUR EVALUER LE NOMBRE DE FACES FRONTIERES
190 : c POUR CHAQUE SOUS DOMAINE
191 : NFRON = 1
192 : DO 1200 IELEM=1,NELEM
193 : NGFACE=MCOGRT.NUM(1,IELEM)
194 : NLFACE=KRFACE.LECT(NGFACE)
195 : IF (NLFACE.EQ.0) THEN
196 : WRITE(IOIMP,*) 'Erreur de programmation n1'
197 : GOTO 9999
198 : ENDIF
199 : NGCGAU=MELEFL.NUM(1,NLFACE)
200 : NGCDRO=MELEFL.NUM(3,NLFACE)
201 : IF ( NGCGAU.EQ.NGCDRO) THEN
202 : NFRON = NFRON + 1
203 : ENDIF
204 : 1200 CONTINUE
205 : c WRITE(6,*) 'NFRON= ',NFRON
206 : NEL1 = NFRON
207 : c WRITE(6,*) 'NPP1= ',NPP1,'NPP2= ',NPP2
208 : c WRITE(6,*) 'NEL1= ',NEL1,'NEL2= ',NEL2
209 :
210 : SEGINI RETRET
211 : SEGINI CTGRT
212 : RETRET.NOMPRI(1:4)='RETN'
213 : RETRET.NOMPRI(5:8)=' '
214 : RETRET.NOMDUA(1:4)='RETN'
215 : RETRET.NOMDUA(5:8)=' '
216 : DO 12 IELEM=1,NELEM
217 : * Le premier point du support de ICOGRT est un point FACE
218 : NGFACE=MCOGRT.NUM(1,IELEM)
219 : NLFACE=KRFACE.LECT(NGFACE)
220 : IF (NLFACE.EQ.0) THEN
221 : WRITE(IOIMP,*) 'Erreur de programmation n1'
222 : GOTO 9999
223 : ENDIF
224 : * On calcule la contribution à la matrice jacobienne IJACO de la face
225 : * NGFAC (points duaux : centres à gauche et à droite de la face)
226 : * (points primaux : une partie (bicoz conditions aux limites)
227 : * de ceux du stencil pour le calcul du gradient
228 : * à la face, ils doivent être des points centres)
229 : * Si le flux de chaleur sur la face est imposé par les conditions
230 : * aux limites, la contribution de la face à IJACO est nulle.
231 : LCTRB1=.TRUE.
232 : c IF (LCLIMQ) THEN
233 : c NLCLQ=KRQIMP.LECT(NGFACE)
234 : c IF (NLCLQ.NE.0) THEN
235 : c LCTRB1=.FALSE.
236 : c ENDIF
237 : c ENDIF
238 : * WRITE(6,*) 'NGFACE= ',NGFACE,'LCTRB1= ', LCTRB1
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 n1'
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 : RETRET.VALMA2(IPD,IPP,IEL2)=0.D0
292 : ELSE
293 : CTGRT.POIPR1(IPP,IEL1)=NPDUAL
294 : RETRET.VALMA1(IPD,IPP,IEL1)=0.D0
295 : ENDIF
296 : ELSE
297 : * Les contributions valent :
298 : * (d Res_{\rho e_t})_d / (d var)_p =
299 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \kappa
300 : * * [ ((n_x * \beta_x) +(n_y * \beta_y)) *
301 : * ((dT)_p / (d var)_p)]
302 : * avec :
303 : * (dT)_p / (d \rho)_p = (1 / \rho_p) *
304 : * ( (((u_p)^2 + (v_p)^2) / (2 * c_v)) - T )
305 : * (dT)_p / (d \rho u)_p = - u_p / (\rho_p * c_v)
306 : * (dT)_p / (d \rho v)_p = - v_p / (\rho_p * c_v)
307 : * (dT)_p / (d \rho e_t)_p = 1 / (\rho_p * c_v)
308 : * \beta_x : coefficients pour le calcul de dT/dx
309 : * \beta_y : coefficients pour le calcul de dT/dy
310 : *
311 : NLCENP=KRCENT.LECT(NPPRIM)
312 : IF (NLCENP.EQ.0) THEN
313 : WRITE(IOIMP,*) 'Erreur grave n2'
314 : GOTO 9999
315 : ENDIF
316 : * normale sortante pour IPD=1, rentrante pour IPD=2
317 : SIGNOR=(-1.D0)**(IPD+1)
318 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
319 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
320 : CNX =MPNORM.VPOCHA(NLFACE,1)
321 : CNY =MPNORM.VPOCHA(NLFACE,2)
322 : BETAX =KCDTDX.VELCHE(IPP+1,IELEM)
323 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA
324 : $ *(BETAX)
325 : DTDRET=1.D0
326 : IF (.NOT.LMUR) THEN
327 : CTGRT.POIPR2(IPP,IEL2)=NPPRIM
328 : RETRET.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRET
329 : c WRITE(6,*) 'IPD=',IPD,'IPP=',IPP,'IEL2=',IEL2
330 : c WRITE(6,*) 'FACT=',FACTOR*DTDRET
331 : c WRITE(6,*) 'SIGNOR=',SIGNOR,'VOLUEL=',VOLUEL,
332 : c & 'SURFFA=',SURFFA,'BETAX= ',BETAX
333 : ELSE
334 : CTGRT.POIPR1(IPP,IEL1)=NPPRIM
335 : RETRET.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRET
336 : c WRITE(6,*) 'IPD=',IPD,'IPP=',IPP,'IEL1=',IEL1
337 : c WRITE(6,*) 'FACT=',FACTOR*DTDRET
338 : c WRITE(6,*) 'SIGNOR=',SIGNOR,'VOLUEL=',VOLUEL,
339 : c & 'SURFFA=',SURFFA,'BETAX= ',BETAX
340 : ENDIF
341 :
342 : ENDIF
343 : 124 CONTINUE
344 : 122 CONTINUE
345 : IF (.NOT.LMUR) THEN
346 : IEL2=IEL2+1
347 : ELSE
348 : IEL1=IEL1+1
349 : ENDIF
350 : ENDIF
351 : 12 CONTINUE
352 : NPP1=NPTEL-1
353 : NPP2=NPTEL-1
354 : NEL1=IEL1-1
355 : NEL2=IEL2-1
356 :
357 : c WRITE(6,*) 'APRES SEGADJ'
358 : c WRITE(6,*) 'NPP1= ',NPP1,'NPP2= ',NPP2
359 : c WRITE(6,*) 'NEL1= ',NEL1,'NEL2= ',NEL2
360 : SEGDES MCOGRT
361 : SEGDES KCDTDX
362 :
363 : SEGADJ RETRET
364 : SEGADJ CTGRT
365 : CTGRT.LMATSI(**)=RETRET
366 : * On accumule les matrices résultantes dans IJACO
366 : CALL AJMTK(CTGRT,IJACO,IMPR,IRET)
367 : IF (IRET.NE.0) GOTO 9999
368 : SEGSUP RETRET
369 : SEGSUP CTGRT
370 : *
371 : 1 CONTINUE
372 : SEGDES ICOGRT
373 : SEGDES MPTEMC
374 : SEGDES MPVOLU
375 : SEGDES MPNORM
376 : SEGDES MPSURF
377 : SEGDES MELEFL
378 : SEGDES KRFACE
379 : SEGDES KRCENT
380 : IF (LCLIMQ) THEN
381 : SEGDES KRQIMP
382 : ENDIF
383 : IF (LCLIMT) THEN
384 : SEGDES KRTIMP
385 : ENDIF
386 : IF (LMIXT) THEN
387 : SEGDES KRMIXT
388 : ENDIF
389 : *
390 : * Normal termination
391 : *
392 : IRET=0
393 : RETURN
394 : *
395 : * Format handling
396 : *
397 : *
398 : * Error handling
399 : *
400 : 9999 CONTINUE
401 : IRET=1
402 : WRITE(IOIMP,*) 'An error was detected in subroutine ylap1b'
403 : RETURN
404 : *
405 : * End of subroutine YLAP1B
406 : *
407 : END
408 :
409 :
410 :
411 :
412 :
413 :
414 :
415 :
© Cast3M 2003 - All rights reserved.
Disclaimer