1 : C XLAP1B SOURCE MAGN 09/07/30 21:21:18 6458
2 : SUBROUTINE XLAP1B(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 : XLAP1B
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 : XLAP1A : 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 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 2D.
66 : C SORTIES : -
67 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
68 : C***********************************************************************
69 : C VERSION : v1, 10/12/2001, version initiale
70 : C HISTORIQUE : v1, 10/12/2001, 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
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 RETRET.MATSIM
118 : *
119 : REAL*8 KAPPA,CV
120 : *
121 : INTEGER IMPR,IRET
122 : *
123 : LOGICAL LCLIMT,LCLIMQ
124 : LOGICAL LMUR
125 : LOGICAL LCTRB1,LCTRB2
126 : *
127 : INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
128 : INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
129 : INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
130 : INTEGER NLCENP,NLCEND,NLFACE,NLCLQ,NLCLT
131 : INTEGER NPTEL
132 : *
133 : REAL*8 BETAX,BETAY,CNX,CNY
134 : REAL*8 SIGNOR,SURFFA,VOLUEL
135 : REAL*8 RHOP,UP,VP,TP
136 : REAL*8 FACTOR
137 : REAL*8 DTDRHO,DTDROU,DTDROV,DTDRET
138 : *
139 : INTEGER ICOORX,NLCGAU,NLCDRO
140 : REAL*8 XF,YF,XG,YG,XFMXG,YFMYG,DRG
141 : & ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH
142 : C
143 : *
144 : * Executable statements
145 : *
146 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans xlap1b.eso'
147 : * On calcule la partie de d Res_{\rho e_t} / d var
148 : * (var prend successivement les valeurs :
149 : * \rho, \rho u, \rho v, \rho e_t)
150 : * faisant intervenir les coefficients pour le calcul des gradients de
151 : * température (ICOGRT).
152 : * C'est la partie contenant le terme : - \vect{q} \pscal \vect{n}
153 : * Les noms de matrices élémentaires (type MATSIM) associées sont :
154 : * RETRHO, RETROU, RETROV, RETRET
155 : IF (LCLIMT) THEN
156 : SEGACT KRTIMP
157 : ENDIF
158 : IF (LCLIMQ) THEN
159 : SEGACT KRQIMP
160 : ENDIF
161 : SEGACT NOMINC
162 : SEGACT KRCENT
163 : SEGACT KRFACE
164 : SEGACT MELEFL
165 : SEGACT MPSURF
166 : SEGACT MPNORM
167 : SEGACT MPVOLU
168 : SEGACT MPKAPC
169 : SEGACT MPCVC
170 : SEGACT MPROC
171 : SEGACT MPVITC
172 : SEGACT MPTEMC
173 : SEGACT ICOGRT
174 : NSOUCH=ICOGRT.IMACHE(/1)
175 : DO 1 ISOUCH=1,NSOUCH
176 : MCOGRT=ICOGRT.IMACHE(ISOUCH)
177 : JCOGRT=ICOGRT.ICHAML(ISOUCH)
178 : SEGACT JCOGRT
179 : KCDTDX=JCOGRT.IELVAL(1)
180 : KCDTDY=JCOGRT.IELVAL(2)
181 : SEGDES JCOGRT
182 : SEGACT KCDTDX
183 : SEGACT KCDTDY
184 : SEGACT MCOGRT
185 : NELEM=MCOGRT.NUM(/2)
186 : NPTEL=MCOGRT.NUM(/1)
187 : NPP1=NPTEL-1
188 : NPP2=NPTEL-1
189 : NEL1=NELEM
190 : NEL2=NELEM
191 : IEL1=1
192 : IEL2=1
193 : SEGINI RETRHO
194 : SEGINI RETROU
195 : SEGINI RETROV
196 : SEGINI RETRET
197 : SEGINI CTGRT
198 : RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
199 : RETRHO.NOMPRI(5:8)=' '
200 : RETRHO.NOMDUA(1:4)=NOMINC.MOTS(4)
201 : RETRHO.NOMDUA(5:8)=' '
202 : RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
203 : RETROU.NOMPRI(5:8)=' '
204 : RETROU.NOMDUA(1:4)=NOMINC.MOTS(4)
205 : RETROU.NOMDUA(5:8)=' '
206 : RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
207 : RETROV.NOMPRI(5:8)=' '
208 : RETROV.NOMDUA(1:4)=NOMINC.MOTS(4)
209 : RETROV.NOMDUA(5:8)=' '
210 : RETRET.NOMPRI(1:4)=NOMINC.MOTS(4)
211 : RETRET.NOMPRI(5:8)=' '
212 : RETRET.NOMDUA(1:4)=NOMINC.MOTS(4)
213 : RETRET.NOMDUA(5:8)=' '
214 : DO 12 IELEM=1,NELEM
215 : * Le premier point du support de ICOGRT est un point FACE
216 : NGFACE=MCOGRT.NUM(1,IELEM)
217 : NLFACE=KRFACE.LECT(NGFACE)
218 : IF (NLFACE.EQ.0) THEN
219 : WRITE(IOIMP,*) 'Erreur de programmation n°1'
220 : GOTO 9999
221 : ENDIF
222 : * On calcule la contribution à la matrice jacobienne IJACO de la face
223 : * NGFAC (points duaux : centres à gauche et à droite de la face)
224 : * (points primaux : une partie (bicoz conditions aux limites)
225 : * de ceux du stencil pour le calcul du gradient
226 : * à la face, ils doivent être des points centres)
227 : * Si le flux de chaleur sur la face est imposé par les conditions
228 : * aux limites, la contribution de la face à IJACO est nulle.
229 : LCTRB1=.TRUE.
230 : IF (LCLIMQ) THEN
231 : NLCLQ=KRQIMP.LECT(NGFACE)
232 : IF (NLCLQ.NE.0) THEN
233 : LCTRB1=.FALSE.
234 : ENDIF
235 : ENDIF
236 : IF (LCTRB1) THEN
237 : NGCGAU=MELEFL.NUM(1,NLFACE)
238 : NGCDRO=MELEFL.NUM(3,NLFACE)
239 : NLCGAU=KRCENT.LECT(NGCGAU)
240 : NLCDRO=KRCENT.LECT(NGCDRO)
241 : LMUR=(NGCGAU.EQ.NGCDRO)
242 : * On distingue le cas où la face est un bord du maillage (mur)
243 : * du cas où la face est interne au maillage
244 : IF (.NOT.LMUR) THEN
245 : NPD=2
246 : ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
247 : XF = MCOORD.XCOOR(ICOORX)
248 : YF = MCOORD.XCOOR(ICOORX+1)
249 : ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
250 : XG = MCOORD.XCOOR(ICOORX)
251 : YG = MCOORD.XCOOR(ICOORX+1)
252 : XFMXG = XF - XG
253 : YFMYG = YF - YG
254 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG))
255 : ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
256 : XD = MCOORD.XCOOR(ICOORX)
257 : YD = MCOORD.XCOOR(ICOORX+1)
258 : XFMXD = XF - XD
259 : YFMYD = YF - YD
260 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD))
261 : ALPHA=DRG/(DRG+DRD)
262 : UMALPH= 1.0D0 - ALPHA
263 : ELSE
264 : NPD=1
265 : ALPHA=0.0D0
266 : UMALPH=1.0D0
267 : ENDIF
268 : KAPPA=UMALPH*MPKAPC.VPOCHA(NLCGAU,1) +
269 : & ALPHA*MPKAPC.VPOCHA(NLCDRO,1)
270 : CV=UMALPH*MPCVC.VPOCHA(NLCGAU,1) +
271 : & ALPHA*MPCVC.VPOCHA(NLCDRO,1)
272 : NPP=NPTEL-1
273 : * IPD=1 : point à gauche du point NGFACE
274 : * IPD=2 : point à droite du point NGFACE
275 : DO 122 IPD=1,NPD
276 : NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
277 : IF (.NOT.LMUR) THEN
278 : CTGRT.POIDU2(IPD,IEL2)=NPDUAL
279 : ELSE
280 : CTGRT.POIDU1(IPD,IEL1)=NPDUAL
281 : ENDIF
282 : NLCEND=KRCENT.LECT(NPDUAL)
283 : IF (NLCEND.EQ.0) THEN
284 : WRITE(IOIMP,*) 'Erreur grave n°1'
285 : GOTO 9999
286 : ENDIF
287 : DO 124 IPP=1,NPP
288 : NPPRIM=MCOGRT.NUM(IPP+1,IELEM)
289 : LCTRB2=.TRUE.
290 : IF (LCLIMT) THEN
291 : NLCLT=KRTIMP.LECT(NPPRIM)
292 : IF (NLCLT.NE.0) THEN
293 : LCTRB2=.FALSE.
294 : ENDIF
295 : ENDIF
296 : IF (.NOT.LCTRB2) THEN
297 : * Lorsque une contribution est nulle, on fixe artificiellement le
298 : * point primal égal au point dual.
299 : IF (.NOT.LMUR) THEN
300 : CTGRT.POIPR2(IPP,IEL2)=NPDUAL
301 : RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
302 : RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
303 : RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
304 : RETRET.VALMA2(IPD,IPP,IEL2)=0.D0
305 : ELSE
306 : CTGRT.POIPR1(IPP,IEL1)=NPDUAL
307 : RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
308 : RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
309 : RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
310 : RETRET.VALMA1(IPD,IPP,IEL1)=0.D0
311 : ENDIF
312 : ELSE
313 : * Les contributions valent :
314 : * (d Res_{\rho e_t})_d / (d var)_p =
315 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \kappa
316 : * * [ ((n_x * \beta_x) +(n_y * \beta_y)) *
317 : * ((dT)_p / (d var)_p)]
318 : * avec :
319 : * (dT)_p / (d \rho)_p = (1 / \rho_p) *
320 : * ( (((u_p)^2 + (v_p)^2) / (2 * c_v)) - T )
321 : * (dT)_p / (d \rho u)_p = - u_p / (\rho_p * c_v)
322 : * (dT)_p / (d \rho v)_p = - v_p / (\rho_p * c_v)
323 : * (dT)_p / (d \rho e_t)_p = 1 / (\rho_p * c_v)
324 : * \beta_x : coefficients pour le calcul de dT/dx
325 : * \beta_y : coefficients pour le calcul de dT/dy
326 : *
327 : NLCENP=KRCENT.LECT(NPPRIM)
328 : IF (NLCENP.EQ.0) THEN
329 : WRITE(IOIMP,*) 'Erreur grave n°2'
330 : GOTO 9999
331 : ENDIF
332 : * normale sortante pour IPD=1, rentrante pour IPD=2
333 : SIGNOR=(-1.D0)**(IPD+1)
334 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
335 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
336 : CNX =MPNORM.VPOCHA(NLFACE,1)
337 : CNY =MPNORM.VPOCHA(NLFACE,2)
338 : BETAX =KCDTDX.VELCHE(IPP+1,IELEM)
339 : BETAY =KCDTDY.VELCHE(IPP+1,IELEM)
340 : RHOP =MPROC.VPOCHA(NLCENP,1)
341 : UP =MPVITC.VPOCHA(NLCENP,1)
342 : VP =MPVITC.VPOCHA(NLCENP,2)
343 : TP =MPTEMC.VPOCHA(NLCENP,1)
344 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*KAPPA
345 : $ *((CNX*BETAX)+(CNY*BETAY))
346 : DTDRHO=((((UP*UP)+(VP*VP))/(2.D0*CV))-TP)/RHOP
347 : DTDROU=-UP/(RHOP*CV)
348 : DTDROV=-VP/(RHOP*CV)
349 : DTDRET=1.D0/(RHOP*CV)
350 : IF (.NOT.LMUR) THEN
351 : CTGRT.POIPR2(IPP,IEL2)=NPPRIM
352 : RETRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRHO
353 : RETROU.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROU
354 : RETROV.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROV
355 : RETRET.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRET
356 : ELSE
357 : CTGRT.POIPR1(IPP,IEL1)=NPPRIM
358 : RETRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRHO
359 : RETROU.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROU
360 : RETROV.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROV
361 : RETRET.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRET
362 : ENDIF
363 : ENDIF
364 : 124 CONTINUE
365 : 122 CONTINUE
366 : IF (.NOT.LMUR) THEN
367 : IEL2=IEL2+1
368 : ELSE
369 : IEL1=IEL1+1
370 : ENDIF
371 : ENDIF
372 : 12 CONTINUE
373 : NPP1=NPTEL-1
374 : NPP2=NPTEL-1
375 : NEL1=IEL1-1
376 : NEL2=IEL2-1
377 : SEGADJ RETRHO
378 : SEGADJ RETROU
379 : SEGADJ RETROV
380 : SEGADJ RETRET
381 : SEGADJ CTGRT
382 : CTGRT.LMATSI(**)=RETRHO
383 : CTGRT.LMATSI(**)=RETROU
384 : CTGRT.LMATSI(**)=RETROV
385 : CTGRT.LMATSI(**)=RETRET
386 : * On accumule les matrices résultantes dans IJACO
386 : CALL AJMTK(CTGRT,IJACO,IMPR,IRET)
387 : IF (IRET.NE.0) GOTO 9999
388 : SEGSUP RETRHO
389 : SEGSUP RETROU
390 : SEGSUP RETROV
391 : SEGSUP RETRET
392 : SEGSUP CTGRT
393 : *
394 : SEGDES MCOGRT
395 : SEGDES KCDTDY
396 : SEGDES KCDTDX
397 : 1 CONTINUE
398 : SEGDES ICOGRT
399 : SEGDES MPTEMC
400 : SEGDES MPVITC
401 : SEGDES MPROC
402 : SEGDES MPCVC
403 : SEGDES MPKAPC
404 : SEGDES MPVOLU
405 : SEGDES MPNORM
406 : SEGDES MPSURF
407 : SEGDES MELEFL
408 : SEGDES KRFACE
409 : SEGDES KRCENT
410 : SEGDES NOMINC
411 : IF (LCLIMQ) THEN
412 : SEGDES KRQIMP
413 : ENDIF
414 : IF (LCLIMT) THEN
415 : SEGDES KRTIMP
416 : ENDIF
417 : *
418 : * Normal termination
419 : *
420 : IRET=0
421 : RETURN
422 : *
423 : * Format handling
424 : *
425 : *
426 : * Error handling
427 : *
428 : 9999 CONTINUE
429 : IRET=1
430 : WRITE(IOIMP,*) 'An error was detected in subroutine xlap1b'
431 : RETURN
432 : *
433 : * End of subroutine XLAP1B
434 : *
435 : END
436 :
437 :
438 :
439 :
440 :
441 :
442 :
443 :
© Cast3M 2003 - All rights reserved.
Disclaimer