1 : C XLAP2D SOURCE MAGN 09/07/30 21:21:34 6458
2 : SUBROUTINE XLAP2D(ICOGRV,MPGRVF,MPROC,MPVITC,
3 : $ MPVOLU,MPNORM,MPSURF,MELEFL,
4 : $ KRFACE,KRCENT,
5 : $ LCLIMV,KRVIMP,
6 : $ LCLITO,KRTOIM,MPTOIM,
7 : $ NOMINC,
8 : $ MPMUC,
9 : $ IJACO,
10 : $ IMPR,IRET)
11 : IMPLICIT INTEGER(I-N)
12 : IMPLICIT REAL*8 (A-H,O-Z)
13 : C***********************************************************************
14 : C NOM : XLAP2D
15 : C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
16 : C VF 3D.
17 : C Ici, on ne calcule que des contributions 'compliquées' à
18 : C la matrice jacobienne faisant intervenir les
19 : C coefficients pour le calcul des gradients de vitesse
20 : C (ICOGRV)
21 : C (contributions à (d Res_{\rho e_t} / d var)
22 : C var prenant successivement les valeurs :
23 : C \rho, \rho u, \rho v, \rho w, \rho e_t )
24 : C Les contributions sont plus "compliquées" à calculer que
25 : C les simples (cf. xlap2c) car on a à dériver des produits
26 : C de fonctions de la vitesse.
27 : C xlap2e calcule aussi une partie des contributions
28 : C 'compliquées'.
29 : C
30 : C NOTE PERSO : On pourrait programmer ça plus lisiblement en stockant
31 : C les contributions dans un tableau de pointeurs (2
32 : C indices, c'est possible ?) et en effectuant des produits
33 : C matriciels (coeff. x matrices de dérivées).
34 : C On n'y coupera pas si on veut regrouper 2D et 3D...
35 : C On ne va pas le faire.
36 : C
37 : C
38 : C LANGAGE : ESOPE
39 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
40 : C mél : gounand@semt2.smts.cea.fr
41 : C***********************************************************************
42 : C APPELES (UTIL) : AJMTK : ajoute un objet de type MATSIM (non
43 : C standard) à un objet de type MATRIK.
44 : C APPELE PAR : XLAP2A : Calcul de la matrice jacobienne du
45 : C résidu du laplacien VF 3D.
46 : C***********************************************************************
47 : C ENTREES : ICOGRV (type MCHELM) : coefficients pour le
48 : C calcul du gradient de la vitesse aux
49 : C interfaces.
50 : C MPGRVF (type MPOVAL) : gradient de la vitesse
51 : C aux interfaces.
52 : C MPROC (type MPOVAL) : masse volumique par
53 : C élément.
54 : C MPVITC (type MPOVAL) : vitesse par élément.
55 : C MPVOLU (type MPOVAL) : volume des éléments.
56 : C MPNORM (type MPOVAL) : normale aux faces.
57 : C MPSURF (type MPOVAL) : surface des faces.
58 : C MELEFL (type MELEME) : connectivités face-(centre
59 : C gauche, centre droit).
60 : C KRFACE (type MLENTI) : tableau de repérage dans
61 : C le maillage des faces des éléments.
62 : C KRCENT (type MLENTI) : tableau de repérage dans
63 : C le maillage des centres des éléments.
64 : C LCLIMV (type logique) : .TRUE. => CL de Dirichlet
65 : C sur la vitesse.
66 : C KRVIMP (type MLENTI) : tableau de repérage dans
67 : C maillage des CL de Dirichlet sur la
68 : C vitesse.
69 : C LCLITO (type logique) : .TRUE. => CL de Dirichlet
70 : C sur le tenseur des contraintes.
71 : C KRTOIM (type MLENTI) : tableau de repérage dans
72 : C maillage des CL de Dirichlet sur le tenseur des
73 : C contraintes
74 : C MPTOIM (type MPOVAL) : valeurs des CL de
75 : C Dirichlet sur le tenseur des contraintes.
76 : C NOMINC (type MLMOTS) : noms des inconnues.
77 : C MPMUC (type MPOVAL) : viscosité dynamique (SI).
78 : C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du
79 : C résidu du laplacien VF 3D.
80 : C SORTIES : -
81 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
82 : C***********************************************************************
83 : C VERSION : v1, 08/03/2002, version initiale
84 : C HISTORIQUE : v1, 08/03/2002, création
85 : C HISTORIQUE :
86 : C HISTORIQUE :
87 : C***********************************************************************
88 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
89 : C en cas de modification de ce sous-programme afin de faciliter
90 : C la maintenance !
91 : C***********************************************************************
92 : -INC CCOPTIO
93 : -INC SMCOORD
94 : -INC SMCHPOI
95 : POINTEUR MPGRVF.MPOVAL
96 : POINTEUR MPMUC.MPOVAL
97 : POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL
98 : POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
99 : POINTEUR MPTOIM.MPOVAL
100 : -INC SMCHAML
101 : POINTEUR ICOGRV.MCHELM,JCOGRV.MCHAML
102 : POINTEUR KDUNDX.MELVAL,KDUNDY.MELVAL,KDUNDZ.MELVAL
103 : -INC SMELEME
104 : POINTEUR MELEFL.MELEME
105 : POINTEUR MCOGRV.MELEME
106 : -INC SMLENTI
107 : POINTEUR KRVIMP.MLENTI,KRTOIM.MLENTI
108 : POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
109 : -INC SMLMOTS
110 : POINTEUR NOMINC.MLMOTS
111 : -INC SMMATRIK
112 : POINTEUR IJACO.MATRIK
113 : *
114 : * Objet matrice élémentaire simplifié
115 : *
116 : SEGMENT GMATSI
117 : INTEGER POIPR1(NPP1,NEL1)
118 : INTEGER POIDU1(1,NEL1)
119 : INTEGER POIPR2(NPP2,NEL2)
120 : INTEGER POIDU2(2,NEL2)
121 : POINTEUR LMATSI(0).MATSIM
122 : ENDSEGMENT
123 : * Contributions compliquées 2 de la part du gradient de
124 : * vitesse (CCGRV2)
125 : POINTEUR CCGRV2.GMATSI
126 : SEGMENT MATSIM
127 : CHARACTER*8 NOMPRI,NOMDUA
128 : REAL*8 VALMA1(1,NPP1,NEL1)
129 : REAL*8 VALMA2(2,NPP2,NEL2)
130 : ENDSEGMENT
131 : POINTEUR RETRHO.MATSIM
132 : POINTEUR RETROU.MATSIM
133 : POINTEUR RETROV.MATSIM
134 : POINTEUR RETROW.MATSIM
135 : *
136 : REAL*8 MU
137 : *
138 : INTEGER IMPR,IRET
139 : *
140 : LOGICAL LCLIMV,LCLITO
141 : LOGICAL LMUR
142 : LOGICAL LCTR2A,LCTR2B
143 : *
144 : INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
145 : INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
146 : INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
147 : INTEGER NLCENP,NLCEND,NLFACE,NLCLV,NLFTOI
148 : INTEGER NPTEL
149 : INTEGER ICOORX,NLCGAU,NLCDRO
150 : *
151 : REAL*8 ALPHA,UMALPH
152 : REAL*8 DRD,DRG
153 : REAL*8 DUXXF,DUXYF,DUXZF
154 : REAL*8 DUYXF,DUYYF,DUYZF
155 : REAL*8 DUZXF,DUZYF,DUZZF
156 : REAL*8 XD,XF,XG,XFMXD,XFMXG
157 : REAL*8 YD,YF,YG,YFMYD,YFMYG
158 : REAL*8 ZD,ZF,ZG,ZFMZD,ZFMZG
159 : REAL*8 ALPHAX,ALPHAY,ALPHAZ
160 : REAL*8 CNX,CNY,CNZ
161 : REAL*8 SIGNOR,SURFFA,VOLUEL
162 : REAL*8 RHOP,UP,VP,WP
163 : REAL*8 FACTOR,FACTDU,FACTDV,FACTDW
164 : REAL*8 DUDRHO,DUDROU
165 : REAL*8 DVDRHO,DVDROV
166 : REAL*8 DWDRHO,DWDROW
167 : REAL*8 DSTDU
168 : REAL*8 TAUXX,TAUYY,TAUZZ
169 : REAL*8 TAUXY,TAUXZ,TAUYZ
170 : REAL*8 EPSIXX,EPSIXY,EPSIXZ
171 : REAL*8 EPSIYX,EPSIYY,EPSIYZ
172 : REAL*8 EPSIZX,EPSIZY,EPSIZZ
173 : REAL*8 GAMMXD,GAMMXG
174 : REAL*8 GAMMYD,GAMMYG
175 : REAL*8 GAMMZD,GAMMZG
176 : *
177 : * Executable statements
178 : *
179 : IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans xlap2d.eso'
180 : * On calcule les contributions à (d Res_{\rho e_t} / d var) ; var
181 : * prenant successivement les valeurs :
182 : * \rho, \rho u, \rho v, \rho w, \rho e_t
183 : * On dérive les termes : (\tens{\tau(\grad{u})} \prod \vect{u})
184 : * \pscal \vect{n}
185 : * ce qui donne deux contributions.
186 : * Contribution 1 :
187 : * ( (d \tens{\tau} / d var) \prod \vect{u}) \pscal \vect{n}
188 : * Contribution 2 :
189 : * ( \tens{\tau} \prod (d \vect{u} / d var)) \pscal \vect{n}
190 : * Note :
191 : * pas de contribution à (d Res_{\rho e_t} / d \rho e_t).
192 : * Les noms de matrices élémentaires (type MATSIM) associées sont :
193 : * RETRHO, RETROU, RETROV, RETROW
194 : IF (LCLIMV) THEN
195 : SEGACT KRVIMP
196 : ENDIF
197 : IF (LCLITO) THEN
198 : SEGACT KRTOIM
199 : SEGACT MPTOIM
200 : ENDIF
201 : SEGACT NOMINC
202 : SEGACT KRCENT
203 : SEGACT KRFACE
204 : SEGACT MELEFL
205 : SEGACT MPSURF
206 : SEGACT MPNORM
207 : SEGACT MPVOLU
208 : SEGACT MPMUC
209 : SEGACT MPROC
210 : SEGACT MPVITC
211 : SEGACT MPGRVF
212 : SEGACT ICOGRV
213 : NSOUCH=ICOGRV.IMACHE(/1)
214 : DO 1 ISOUCH=1,NSOUCH
215 : MCOGRV=ICOGRV.IMACHE(ISOUCH)
216 : JCOGRV=ICOGRV.ICHAML(ISOUCH)
217 : SEGACT JCOGRV
218 : KDUNDX=JCOGRV.IELVAL(1)
219 : KDUNDY=JCOGRV.IELVAL(2)
220 : KDUNDZ=JCOGRV.IELVAL(3)
221 : SEGDES JCOGRV
222 : SEGACT KDUNDX
223 : SEGACT KDUNDY
224 : SEGACT KDUNDZ
225 : SEGACT MCOGRV
226 : NELEM=MCOGRV.NUM(/2)
227 : NPTEL=MCOGRV.NUM(/1)
228 : NPP1=NPTEL
229 : NPP2=NPTEL+1
230 : NEL1=NELEM
231 : NEL2=NELEM
232 : IEL1=1
233 : IEL2=1
234 : SEGINI RETRHO
235 : SEGINI RETROU
236 : SEGINI RETROV
237 : SEGINI RETROW
238 : SEGINI CCGRV2
239 : RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
240 : RETRHO.NOMPRI(5:8)=' '
241 : RETRHO.NOMDUA(1:4)=NOMINC.MOTS(5)
242 : RETRHO.NOMDUA(5:8)=' '
243 : RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
244 : RETROU.NOMPRI(5:8)=' '
245 : RETROU.NOMDUA(1:4)=NOMINC.MOTS(5)
246 : RETROU.NOMDUA(5:8)=' '
247 : RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
248 : RETROV.NOMPRI(5:8)=' '
249 : RETROV.NOMDUA(1:4)=NOMINC.MOTS(5)
250 : RETROV.NOMDUA(5:8)=' '
251 : RETROW.NOMPRI(1:4)=NOMINC.MOTS(4)
252 : RETROW.NOMPRI(5:8)=' '
253 : RETROW.NOMDUA(1:4)=NOMINC.MOTS(5)
254 : RETROW.NOMDUA(5:8)=' '
255 : DO 12 IELEM=1,NELEM
256 : * Le premier point du support de ICOGRV est un point FACE
257 : NGFACE=MCOGRV.NUM(1,IELEM)
258 : NLFACE=KRFACE.LECT(NGFACE)
259 : IF (NLFACE.EQ.0) THEN
260 : WRITE(IOIMP,*) 'Erreur de programmation n°1'
261 : GOTO 9999
262 : ENDIF
263 : * On calcule la contribution 2 à la matrice jacobienne IJACO de la
264 : * face NGFACE
265 : * (points duaux : centres à gauche et à droite de la face)
266 : * (points primaux : une partie (bicoz conditions aux limites)
267 : * de ceux du stencil pour le calcul du gradient
268 : * à la face, ils doivent être des points centres
269 : * ET les centres à gauche et à droite qui servent pour le
270 : * calcul de la vitesse sur la face)
271 : * Si la vitesse sur la face est imposée par les
272 : * conditions aux limites, la contribution 2 de la face à IJACO est
273 : * nulle.
274 : LCTR2A=.TRUE.
275 : IF (LCLIMV) THEN
276 : NLCLV=KRVIMP.LECT(NGFACE)
277 : IF (NLCLV.NE.0) THEN
278 : LCTR2A=.FALSE.
279 : ENDIF
280 : ENDIF
281 : IF (LCTR2A) THEN
282 : NGCGAU=MELEFL.NUM(1,NLFACE)
283 : NGCDRO=MELEFL.NUM(3,NLFACE)
284 : NLCGAU=KRCENT.LECT(NGCGAU)
285 : NLCDRO=KRCENT.LECT(NGCDRO)
286 : LMUR=(NGCGAU.EQ.NGCDRO)
287 : * On calcule ensuite les valeurs des coefficients qui servent pour
288 : * interpoler la vitesse sur la face NGFACE à partir des vitesses
289 : * gauche (et droite si non mur) et des dérivées de la vitesse sur la
290 : * face NGFACE (de la même manière que dans xlap12.eso)
291 : IF (LMUR) THEN
292 : * Parametres geometriques
293 : ICOORX=((IDIM+1)*(NGFACE-1))+1
294 : XF=MCOORD.XCOOR(ICOORX)
295 : YF=MCOORD.XCOOR(ICOORX+1)
296 : ZF=MCOORD.XCOOR(ICOORX+2)
297 : ICOORX=((IDIM+1)*(NGCGAU-1))+1
298 : XG=MCOORD.XCOOR(ICOORX)
299 : YG=MCOORD.XCOOR(ICOORX+1)
300 : ZG=MCOORD.XCOOR(ICOORX+2)
301 : XFMXG=XF-XG
302 : YFMYG=YF-YG
303 : ZFMZG=ZF-ZG
304 : ALPHA=0.0D0
305 : UMALPH=1.0D0
306 : * coeffs pour calculer la vitesse
307 : GAMMXG=1.D0
308 : EPSIXX=XFMXG
309 : EPSIXY=YFMYG
310 : EPSIXZ=ZFMZG
311 : GAMMYG=1.D0
312 : EPSIYX=XFMXG
313 : EPSIYY=YFMYG
314 : EPSIYZ=ZFMZG
315 : GAMMZG=1.D0
316 : EPSIZX=XFMXG
317 : EPSIZY=YFMYG
318 : EPSIZZ=ZFMZG
319 : ELSEIF (.NOT.LMUR) THEN
320 : ICOORX=((IDIM+1)*(NGFACE-1))+1
321 : XF=MCOORD.XCOOR(ICOORX)
322 : YF=MCOORD.XCOOR(ICOORX+1)
323 : ZF=MCOORD.XCOOR(ICOORX+2)
324 : ICOORX=((IDIM+1)*(NGCGAU-1))+1
325 : XG=MCOORD.XCOOR(ICOORX)
326 : YG=MCOORD.XCOOR(ICOORX+1)
327 : ZG=MCOORD.XCOOR(ICOORX+2)
328 : XFMXG=XF-XG
329 : YFMYG=YF-YG
330 : ZFMZG=ZF-ZG
331 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
332 : ICOORX=((IDIM+1)*(NGCDRO-1))+1
333 : XD=MCOORD.XCOOR(ICOORX)
334 : YD=MCOORD.XCOOR(ICOORX+1)
335 : ZD=MCOORD.XCOOR(ICOORX+2)
336 : XFMXD=XF-XD
337 : YFMYD=YF-YD
338 : ZFMZD=ZF-ZD
339 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
340 : ALPHA=DRG/(DRG+DRD)
341 : UMALPH=1.0D0-ALPHA
342 : GAMMXG=UMALPH
343 : GAMMXD=ALPHA
344 : EPSIXX=(XFMXG*UMALPH)+(XFMXD*ALPHA)
345 : EPSIXY=(YFMYG*UMALPH)+(YFMYD*ALPHA)
346 : EPSIXZ=(ZFMZG*UMALPH)+(ZFMZD*ALPHA)
347 : GAMMYG=UMALPH
348 : GAMMYD=ALPHA
349 : EPSIYX=(XFMXG*UMALPH)+(XFMXD*ALPHA)
350 : EPSIYY=(YFMYG*UMALPH)+(YFMYD*ALPHA)
351 : EPSIYZ=(ZFMZG*UMALPH)+(ZFMZD*ALPHA)
352 : GAMMZG=UMALPH
353 : GAMMZD=ALPHA
354 : EPSIZX=(XFMXG*UMALPH)+(XFMXD*ALPHA)
355 : EPSIZY=(YFMYG*UMALPH)+(YFMYD*ALPHA)
356 : EPSIZZ=(ZFMZG*UMALPH)+(ZFMZD*ALPHA)
357 : ELSE
358 : WRITE(IOIMP,*) 'Erreur de programmation n°2'
359 : GOTO 9999
360 : ENDIF
361 : * On distingue le cas où la face est un bord du maillage (mur)
362 : * du cas où la face est interne au maillage
363 : IF (.NOT.LMUR) THEN
364 : NPD=2
365 : NPP=NPTEL+1
366 : ELSE
367 : NPD=1
368 : NPP=NPTEL
369 : ENDIF
370 : MU=UMALPH*MPMUC.VPOCHA(NLCGAU,1) +
371 : & ALPHA*MPMUC.VPOCHA(NLCDRO,1)
372 : * On calcule tout d'abord les valeurs du tenseur des contraintes sur
373 : * la face NGFACE de la même manière que dans xlap12.eso
374 : IF (LCLITO) THEN
375 : NLFTOI=KRTOIM.LECT(NGFACE)
376 : ELSE
377 : NLFTOI=0
378 : ENDIF
379 : IF (NLFTOI.NE.0) THEN
380 : TAUXX=MPTOIM.VPOCHA(NLFTOI,1)
381 : TAUYY=MPTOIM.VPOCHA(NLFTOI,2)
382 : TAUZZ=MPTOIM.VPOCHA(NLFTOI,3)
383 : TAUXY=MPTOIM.VPOCHA(NLFTOI,4)
384 : TAUXZ=MPTOIM.VPOCHA(NLFTOI,5)
385 : TAUYZ=MPTOIM.VPOCHA(NLFTOI,6)
386 : ELSE
387 : DUXXF=MPGRVF.VPOCHA(NLFACE,1)
388 : DUXYF=MPGRVF.VPOCHA(NLFACE,2)
389 : DUXZF=MPGRVF.VPOCHA(NLFACE,3)
390 : DUYXF=MPGRVF.VPOCHA(NLFACE,4)
391 : DUYYF=MPGRVF.VPOCHA(NLFACE,5)
392 : DUYZF=MPGRVF.VPOCHA(NLFACE,6)
393 : DUZXF=MPGRVF.VPOCHA(NLFACE,7)
394 : DUZYF=MPGRVF.VPOCHA(NLFACE,8)
395 : DUZZF=MPGRVF.VPOCHA(NLFACE,9)
396 : DSTDU=2.0D0/3.0D0*(DUXXF+DUYYF+DUZZF)
397 : TAUXX=MU*(2.0D0*DUXXF-DSTDU)
398 : TAUYY=MU*(2.0D0*DUYYF-DSTDU)
399 : TAUZZ=MU*(2.0D0*DUZZF-DSTDU)
400 : TAUXY=MU*(DUXYF+DUYXF)
401 : TAUXZ=MU*(DUXZF+DUZXF)
402 : TAUYZ=MU*(DUZYF+DUYZF)
403 : ENDIF
404 : * IPD=1 : point à gauche du point NGFACE
405 : * IPD=2 : point à droite du point NGFACE
406 : DO 122 IPD=1,NPD
407 : NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
408 : IF (.NOT.LMUR) THEN
409 : CCGRV2.POIDU2(IPD,IEL2)=NPDUAL
410 : ELSE
411 : CCGRV2.POIDU1(IPD,IEL1)=NPDUAL
412 : ENDIF
413 : NLCEND=KRCENT.LECT(NPDUAL)
414 : IF (NLCEND.EQ.0) THEN
415 : WRITE(IOIMP,*) 'Erreur grave n°1'
416 : GOTO 9999
417 : ENDIF
418 : DO 124 IPP=1,NPP
419 : IF (IPP.EQ.NPTEL) THEN
420 : NPPRIM=NGCGAU
421 : ELSEIF (IPP.EQ.NPTEL+1) THEN
422 : NPPRIM=NGCDRO
423 : ELSEIF (IPP.GE.1.AND.IPP.LT.NPTEL) THEN
424 : NPPRIM=MCOGRV.NUM(IPP+1,IELEM)
425 : ELSE
426 : WRITE(IOIMP,*) 'Erreur grave n°2'
427 : GOTO 9999
428 : ENDIF
429 : LCTR2B=.TRUE.
430 : IF (LCLIMV) THEN
431 : NLCLV=KRVIMP.LECT(NPPRIM)
432 : IF (NLCLV.NE.0) THEN
433 : LCTR2B=.FALSE.
434 : ENDIF
435 : ENDIF
436 : IF (.NOT.LCTR2B) THEN
437 : * Lorsque une contribution est nulle, on fixe artificiellement le
438 : * point primal égal au point dual.
439 : IF (.NOT.LMUR) THEN
440 : CCGRV2.POIPR2(IPP,IEL2)=NPDUAL
441 : RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
442 : RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
443 : RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
444 : RETROW.VALMA2(IPD,IPP,IEL2)=0.D0
445 : ELSE
446 : CCGRV2.POIPR1(IPP,IEL1)=NPDUAL
447 : RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
448 : RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
449 : RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
450 : RETROW.VALMA1(IPD,IPP,IEL1)=0.D0
451 : ENDIF
452 : ELSE
453 : * Les contributions 2 valent :
454 : * (d Res_{\rho e_t})_d / (d var)_p =
455 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f)
456 : * * [ [ [ ((n_x * \tau_xx) + (n_y * \tau_yx) + (n_z * \tau_zx))
457 : * * ((\epsilon_xx * \alpha_x) + (\epsilon_xy * \alpha_y)
458 : * + (\epsilon_xz * \alpha_z))
459 : * ]
460 : * * ((du)_p / (d var)_p)
461 : * ]
462 : * + [ [ ((n_x * \tau_xy) + (n_y * \tau_yy) + (n_z * \tau_zy))
463 : * * ((\epsilon_yx * \alpha_x) + (\epsilon_yy * \alpha_y)
464 : * + (\epsilon_yz * \alpha_z))
465 : * ]
466 : * * ((dv)_p / (d var)_p)
467 : * ]
468 : * + [ [ ((n_x * \tau_xz) + (n_y * \tau_yz) + (n_z * \tau_zz))
469 : * * ((\epsilon_zx * \alpha_x) + (\epsilon_zy * \alpha_y)
470 : * + (\epsilon_zz * \alpha_z))
471 : * ]
472 : * * ((dw)_p / (d var)_p)
473 : * ]
474 : * ]
475 : * et (m étant le centre de gauche dans le cas mur
476 : * et les centre gauche, puis droite dans le cas normal)
477 : * (d Res_{\rho e_t})_d / (d var)_m =
478 : * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f)
479 : * * [ [ [ ((n_x * \tau_xx) + (n_y * \tau_yx) + (n_z * \tau_zx))
480 : * * (\gamma_x)
481 : * ]
482 : * * ((du)_p / (d var)_p)
483 : * ]
484 : * + [ [ ((n_x * \tau_xy) + (n_y * \tau_yy) + (n_z * \tau_zy))
485 : * * (\gamma_y)
486 : * ]
487 : * * ((dv)_p / (d var)_p)
488 : * ]
489 : * + [ [ ((n_x * \tau_xz) + (n_y * \tau_yz) + (n_z * \tau_zz))
490 : * * (\gamma_z)
491 : * ]
492 : * * ((dw)_p / (d var)_p)
493 : * ]
494 : * ]
495 : * sachant que l'expression donnant l'interpolation de la
496 : * vitesse (u,v) sur une face i est de la forme :
497 : * u_i = (\gamma_x,gauche * u_gauche)
498 : * + (\gamma_x,droite * u_droite) (=0 dans le cas mur)
499 : * + (\epsilon_xx,i * (du/dx)_i)
500 : * + (\epsilon_xy,i * (du/dy)_i)
501 : * + (\epsilon_xz,i * (du/dz)_i)
502 : * v_i = (\gamma_y,gauche * v_gauche)
503 : * + (\gamma_y,droite * v_droite) (=0 dans le cas mur)
504 : * + (\epsilon_yx,i * (dv/dx)_i)
505 : * + (\epsilon_yy,i * (dv/dy)_i)
506 : * + (\epsilon_yz,i * (dv/dz)_i)
507 : * w_i = (\gamma_z,gauche * w_gauche)
508 : * + (\gamma_z,droite * w_droite) (=0 dans le cas mur)
509 : * + (\epsilon_zx,i * (dw/dx)_i)
510 : * + (\epsilon_zy,i * (dw/dy)_i)
511 : * + (\epsilon_zz,i * (dw/dz)_i)
512 : * (voir dans xlap13.eso et ci-dessus pour les valeurs des coeffs
513 : * \gamma et \epsilon)
514 : NLCENP=KRCENT.LECT(NPPRIM)
515 : IF (NLCENP.EQ.0) THEN
516 : WRITE(IOIMP,*) 'Erreur grave n°3'
517 : GOTO 9999
518 : ENDIF
519 : * normale sortante pour IPD=1, rentrante pour IPD=2
520 : SIGNOR=(-1.D0)**(IPD+1)
521 : VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
522 : SURFFA=MPSURF.VPOCHA(NLFACE,1)
523 : CNX =MPNORM.VPOCHA(NLFACE,1)
524 : CNY =MPNORM.VPOCHA(NLFACE,2)
525 : CNZ =MPNORM.VPOCHA(NLFACE,3)
526 : RHOP =MPROC.VPOCHA(NLCENP,1)
527 : UP =MPVITC.VPOCHA(NLCENP,1)
528 : VP =MPVITC.VPOCHA(NLCENP,2)
529 : WP =MPVITC.VPOCHA(NLCENP,3)
530 : FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA
531 : IF (IPP.EQ.NPTEL) THEN
532 : FACTDU= ((CNX*TAUXX)+(CNY*TAUXY)+(CNZ*TAUXZ))
533 : $ * GAMMXG
534 : FACTDV= ((CNX*TAUXY)+(CNY*TAUYY)+(CNZ*TAUYZ))
535 : $ * GAMMYG
536 : FACTDW= ((CNX*TAUXZ)+(CNY*TAUYZ)+(CNZ*TAUZZ))
537 : $ * GAMMZG
538 : ELSEIF (IPP.EQ.NPTEL+1) THEN
539 : FACTDU= ((CNX*TAUXX)+(CNY*TAUXY)+(CNZ*TAUXZ))
540 : $ * GAMMXD
541 : FACTDV= ((CNX*TAUXY)+(CNY*TAUYY)+(CNZ*TAUYZ))
542 : $ * GAMMYD
543 : FACTDW= ((CNX*TAUXZ)+(CNY*TAUYZ)+(CNZ*TAUZZ))
544 : $ * GAMMZD
545 : ELSEIF (IPP.GE.1.AND.IPP.LT.NPTEL) THEN
546 : ALPHAX=KDUNDX.VELCHE(IPP+1,IELEM)
547 : ALPHAY=KDUNDY.VELCHE(IPP+1,IELEM)
548 : ALPHAZ=KDUNDZ.VELCHE(IPP+1,IELEM)
549 : FACTDU= ((CNX*TAUXX)+(CNY*TAUXY)+(CNZ*TAUXZ))
550 : $ * ((EPSIXX*ALPHAX)+(EPSIXY*ALPHAY)
551 : $ +(EPSIXZ*ALPHAZ))
552 : FACTDV= ((CNX*TAUXY)+(CNY*TAUYY)+(CNZ*TAUYZ))
553 : $ * ((EPSIYX*ALPHAX)+(EPSIYY*ALPHAY)
554 : $ +(EPSIYZ*ALPHAZ))
555 : FACTDW= ((CNX*TAUXZ)+(CNY*TAUYZ)+(CNZ*TAUZZ))
556 : $ * ((EPSIZX*ALPHAX)+(EPSIZY*ALPHAY)
557 : $ +(EPSIZZ*ALPHAZ))
558 : ELSE
559 : WRITE(IOIMP,*) 'Erreur grave n°3'
560 : GOTO 9999
561 : ENDIF
562 : DUDRHO=-UP /RHOP
563 : DUDROU=1.D0/RHOP
564 : DVDRHO=-VP /RHOP
565 : DVDROV=1.D0/RHOP
566 : DWDRHO=-WP /RHOP
567 : DWDROW=1.D0/RHOP
568 : IF (.NOT.LMUR) THEN
569 : CCGRV2.POIPR2(IPP,IEL2)=NPPRIM
570 : RETRHO.VALMA2(IPD,IPP,IEL2)=
571 : $ FACTOR*( (FACTDU*DUDRHO)
572 : $ +(FACTDV*DVDRHO)
573 : $ +(FACTDW*DWDRHO))
574 : RETROU.VALMA2(IPD,IPP,IEL2)=
575 : $ FACTOR* (FACTDU*DUDROU)
576 : RETROV.VALMA2(IPD,IPP,IEL2)=
577 : $ FACTOR* (FACTDV*DVDROV)
578 : RETROW.VALMA2(IPD,IPP,IEL2)=
579 : $ FACTOR* (FACTDW*DWDROW)
580 : ELSE
581 : CCGRV2.POIPR1(IPP,IEL1)=NPPRIM
582 : RETRHO.VALMA1(IPD,IPP,IEL1)=
583 : $ FACTOR*( (FACTDU*DUDRHO)
584 : $ +(FACTDV*DVDRHO)
585 : $ +(FACTDW*DWDRHO))
586 : RETROU.VALMA1(IPD,IPP,IEL1)=
587 : $ FACTOR* (FACTDU*DUDROU)
588 : RETROV.VALMA1(IPD,IPP,IEL1)=
589 : $ FACTOR* (FACTDV*DVDROV)
590 : RETROW.VALMA1(IPD,IPP,IEL1)=
591 : $ FACTOR* (FACTDW*DWDROW)
592 : ENDIF
593 : ENDIF
594 : 124 CONTINUE
595 : 122 CONTINUE
596 : IF (.NOT.LMUR) THEN
597 : IEL2=IEL2+1
598 : ELSE
599 : IEL1=IEL1+1
600 : ENDIF
601 : ENDIF
602 : 12 CONTINUE
603 : NPP1=NPTEL
604 : NPP2=NPTEL+1
605 : NEL1=IEL1-1
606 : NEL2=IEL2-1
607 : SEGADJ RETRHO
608 : SEGADJ RETROU
609 : SEGADJ RETROV
610 : SEGADJ RETROW
611 : SEGADJ CCGRV2
612 : CCGRV2.LMATSI(**)=RETRHO
613 : CCGRV2.LMATSI(**)=RETROU
614 : CCGRV2.LMATSI(**)=RETROV
615 : CCGRV2.LMATSI(**)=RETROW
616 : * On accumule les matrices résultantes dans IJACO
616 : CALL AJMTK(CCGRV2,IJACO,IMPR,IRET)
617 : IF (IRET.NE.0) GOTO 9999
618 : SEGSUP RETRHO
619 : SEGSUP RETROU
620 : SEGSUP RETROV
621 : SEGSUP RETROW
622 : SEGSUP CCGRV2
623 : *
624 : SEGDES MCOGRV
625 : SEGDES KDUNDZ
626 : SEGDES KDUNDY
627 : SEGDES KDUNDX
628 : 1 CONTINUE
629 : SEGDES ICOGRV
630 : SEGDES MPGRVF
631 : SEGDES MPVITC
632 : SEGDES MPMUC
633 : SEGDES MPROC
634 : SEGDES MPVOLU
635 : SEGDES MPNORM
636 : SEGDES MPSURF
637 : SEGDES MELEFL
638 : SEGDES KRFACE
639 : SEGDES KRCENT
640 : SEGDES NOMINC
641 : IF (LCLITO) THEN
642 : SEGDES KRTOIM
643 : SEGDES MPTOIM
644 : ENDIF
645 : IF (LCLIMV) THEN
646 : SEGDES KRVIMP
647 : ENDIF
648 : *
649 : * Normal termination
650 : *
651 : IRET=0
652 : RETURN
653 : *
654 : * Format handling
655 : *
656 : *
657 : * Error handling
658 : *
659 : 9999 CONTINUE
660 : IRET=1
661 : WRITE(IOIMP,*) 'An error was detected in subroutine xlap2d'
662 : RETURN
663 : *
664 : * End of subroutine XLAP2D
665 : *
666 : END
667 :
668 :
669 :
670 :
671 :
672 :
673 :
674 :
© Cast3M 2003 - All rights reserved.
Disclaimer