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