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