1 : C TLAP13 SOURCE CHAT 05/01/13 03:41:41 5004
2 : SUBROUTINE TLAP13(NESP,IMUC,IMTC,TSIGK,TSIGE,IGRKEP,
3 : & MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM,
4 : & ICHFLU,DT)
5 : C***********************************************************************
6 : C NOM : TLAP13
7 : C DESCRIPTION : Calcul des flux diffusifs pour k-\eps equations de
8 : C turbulence (3D)
9 : C
10 : C
11 : C LANGAGE : ESOPE
12 : C AUTEUR : S. Kudriakov (CEA/DEN/DM2S/SFME/LTMF)
13 : C mél : skudriakov@cea.fr
14 : C***********************************************************************
15 : C APPELES (UTIL) : KRIPAD : MELEME -> (num. globale->locale)
16 : C LICHT : Lecture des pointeurs (maillages, valeurs)
17 : C d'un objet de type MCHPOI.
18 : C ERREUR : Gestion des erreurs par GIBI.
19 : C APPELE PAR : ZLAP11 : Chapeau de l'opérateur Gibiane 'LAPN'
20 : C option 'VF'.
21 : C***********************************************************************
22 : C ENTREES : NESP (type ENTIER) : number of species explicitly
23 : C treated in the NS equations
24 : C IMUC (type MCHPOI) : laminar viscosity
25 : C IMTC (type MCHPOI) : turbulent viscosity
26 : C TSIGK (type REEL) : turbulent constant \sigma_k
27 : C TSIGE (type REEL) : turbulent constant
28 : C \sigma_{\eps}
29 : C IGRKEP (type MCHPOI) : gradients of k and \epsilon
30 : C
31 : C MELEMC (type MELEME) : maillage des centres des
32 : C éléments.
33 : C MELEMF (type MELEME) : maillage des faces des
34 : C éléments.
35 : C MELEFL (type MELEME) : connectivités face-(centre
36 : C gauche, centre droit).
37 : C ISURF (type MCHPOI) : surface des faces.
38 : C INORM (type MCHPOI) : normale aux faces.
39 : C IDIAM (type MCHPOI) : diamètre des éléments.
40 : C ENTREES/SORTIES : -
41 : C SORTIES : ICHFLU (type MCHPOI) : flux diffusif aux
42 : C interfaces.
43 : C DT (type REAL*8) : pas de temps de stabilité
44 : C (Fourier)
45 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
46 : C***********************************************************************
47 : C VERSION : v1, 3/12/2003, version initiale
48 : C HISTORIQUE :
49 : C***********************************************************************
50 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
51 : C en cas de modification de ce sous-programme afin de faciliter
52 : C la maintenance !
53 : C***********************************************************************
54 : IMPLICIT INTEGER(I-N)
55 : -INC CCOPTIO
56 : -INC CCREEL
57 : -INC SMCHPOI
58 : -INC SMELEME
59 : -INC SMCOORD
60 : -INC SMLENTI
61 : -INC SMLMOTS
62 : C
63 : POINTEUR IMUC.MCHPOI,IMTC.MCHPOI,IGRKEP.MCHPOI
64 : POINTEUR MPMUC.MPOVAL,MPMTC.MPOVAL,MPGRKE.MPOVAL,
65 : & MPSURF.MPOVAL, MPNORM.MPOVAL, MPDIAM.MPOVAL,
66 : & MPFLUX.MPOVAL
67 : C
68 : POINTEUR MELEMC.MELEME,MELEMF.MELEME,MELEFL.MELEME
69 : POINTEUR MLCENT.MLENTI,MLETIM.MLENTI,MLERIM.MLENTI
70 : C
71 : INTEGER ITIMP
72 : & ,ISURF,INORM,IDIAM,ICHFLU
73 : & ,NFAC, NLCF, NGCF, NGCF1, NGCEG, NGCED
74 : & ,NLCEG,NLCED,NLFTI,NLFRI
75 : & , ICOORX, IGEOM
76 :
77 : REAL*8 DT, UNSDT
78 : & ,XG,YG,ZG,XFMXG,YFMYG,ZFMZG,DRG
79 : & ,XD,YD,ZD,XFMXD,YFMYD,ZFMZD,DRD,ALPHA,UMALPH
80 : & ,XF,YF,ZF
81 : & ,CNX, CNY, CNZ, ORIENT, DIAM, DIAM2, CELL, SURF
82 : INTEGER NESP
83 : REAL*8 DKADX,DKADY,DKADZ,DEPDX,DEPDY,DEPDZ
84 : REAL*8 MU,MUT,TSIGK,TSIGE
85 : C
86 : CHARACTER*8 TYPE
87 : C
88 : C**** Initialisation de 1/DT
89 : C
90 : UNSDT = 0.0D0
91 : C
92 : C**** KRIPAD pour la correspondance global/local de centre
93 : C
93 : CALL KRIPAD(MELEMC,MLCENT)
94 : C
95 : C EN KRIPAD
96 : C SEGACT MELEMC
97 : C SEGACT MLCENT
98 : C
98 : CALL LICHT(IMUC,MPMUC,TYPE,IGEOM)
98 : CALL LICHT(IMTC,MPMTC,TYPE,IGEOM)
98 : CALL LICHT(IGRKEP,MPGRKE,TYPE,IGEOM)
98 : CALL LICHT(ISURF,MPSURF,TYPE,IGEOM)
98 : CALL LICHT(INORM,MPNORM,TYPE,IGEOM)
98 : CALL LICHT(IDIAM,MPDIAM,TYPE,IGEOM)
98 : CALL LICHT(ICHFLU,MPFLUX,TYPE,IGEOM)
99 : C
100 : C EN LICHT
101 : C SEGACT*MOD MPROC
102 : C SEGACT*MOD MPTEMC
103 : C SEGACT*MOD MPSURF
104 : C SEGACT*MOD MPNORM
105 : C SEGACT*MOD MPDIAM
106 : C SEGACT*MOD MPFLUX
107 : C
108 : C---------------------------------------------------------
109 : SEGACT MELEFL
110 : SEGACT MELEMF
111 : NFAC = MELEMF.NUM(/2)
112 : C
113 : C**** Boucle sur les faces
114 : C
115 : DO NLCF = 1, NFAC, 1
116 : C
117 : C******* NLCF = numero local du centre de facel
118 : C NGCF = numero global du centre de facel
119 : C NLCF1 = numero local du centre de face
120 : C NGCEG = numero global du centre ELT "gauche"
121 : C NLCEG = numero local du centre ELT "gauche"
122 : C NGCED = numero global du centre ELT "droite"
123 : C NLCED = numero local du centre ELT "droite"
124 : C
125 : NGCF = MELEMF.NUM(1,NLCF)
126 : NGCF1 = MELEFL.NUM(2,NLCF)
127 : IF(NGCF .NE. NGCF1)THEN
128 : MOTERR(1:40)= 'FACEL et FACE = ? '
128 : CALL ERREUR(5)
129 : GOTO 9999
130 : ENDIF
131 : C
132 : NGCEG = MELEFL.NUM(1,NLCF)
133 : NGCED = MELEFL.NUM(3,NLCF)
134 : NLCEG = MLCENT.LECT(NGCEG)
135 : NLCED = MLCENT.LECT(NGCED)
136 : C-------------------------------------------------------
137 : IF(NGCEG .NE. NGCED)THEN
138 : C
139 : C********** Parametres geometriques
140 : C
141 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1
142 : XF = MCOORD.XCOOR(ICOORX)
143 : YF = MCOORD.XCOOR(ICOORX+1)
144 : ZF = MCOORD.XCOOR(ICOORX+2)
145 : C
146 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
147 : XG = MCOORD.XCOOR(ICOORX)
148 : YG = MCOORD.XCOOR(ICOORX+1)
149 : ZG = MCOORD.XCOOR(ICOORX+2)
150 : XFMXG = XF - XG
151 : YFMYG = YF - YG
152 : ZFMZG = ZF - ZG
153 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
154 : C
155 : ICOORX = ((IDIM + 1) * (NGCED - 1))+1
156 : XD = MCOORD.XCOOR(ICOORX)
157 : YD = MCOORD.XCOOR(ICOORX+1)
158 : ZD = MCOORD.XCOOR(ICOORX+2)
159 : XFMXD = XF - XD
160 : YFMYD = YF - YD
161 : ZFMZD = ZF - ZD
162 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
163 : C
164 : C********** F=G -> DRG = 0 -> ALPHA = 0
165 : C
166 : ALPHA=DRG/(DRG+DRD)
167 : UMALPH= 1.0D0 - ALPHA
168 : C
169 : C********** Les valeurs à l'interface
170 : C
171 : C DRG=0 -> F=G
172 : C
173 : ELSE
174 : C
175 : C********** MURS
176 : C
177 : C Etat a gauche = Etat droite
178 : C
179 : ALPHA=0.0D0
180 : UMALPH=1.0D0
181 : C
182 : C********** Parametres geometriques
183 : C
184 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1
185 : XF = MCOORD.XCOOR(ICOORX)
186 : YF = MCOORD.XCOOR(ICOORX+1)
187 : ZF = MCOORD.XCOOR(ICOORX+2)
188 : C
189 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
190 : XG = MCOORD.XCOOR(ICOORX)
191 : YG = MCOORD.XCOOR(ICOORX+1)
192 : ZG = MCOORD.XCOOR(ICOORX+2)
193 : XFMXG = XF - XG
194 : YFMYG = YF - YG
195 : ZFMZG = ZF - ZG
196 : ENDIF
197 : C
198 : C******* On calcule le sign du pruduit scalare
199 : C (Normales de Castem) * (vecteur "gauche" -> "centre")
200 : C
201 : CNX = MPNORM.VPOCHA(NLCF,1)
202 : CNY = MPNORM.VPOCHA(NLCF,2)
203 : CNZ = MPNORM.VPOCHA(NLCF,3)
204 : ORIENT=(CNX * XFMXG)+(CNY * YFMYG)+(CNZ*ZFMZG)
205 : ORIENT = SIGN(1.0D0,ORIENT)
206 : IF(ORIENT .NE. 1.0D0)THEN
207 : MOTERR(1:40)=
208 : & 'LAPN , subroutine zlap12.eso. '
209 : WRITE(IOIMP,*) MOTERR(1:40)
209 : CALL ERREUR(5)
210 : GOTO 9999
211 : ENDIF
212 : C
213 : C******* Les flux aux interfaces
214 : C
215 : SURF = MPSURF.VPOCHA(NLCF,1)
216 : DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) +
217 : & ALPHA*MPDIAM.VPOCHA(NLCED,1)
218 : DIAM2=DIAM*DIAM
219 : MU = UMALPH*MPMUC.VPOCHA(NLCEG,1) +
220 : & ALPHA*MPMUC.VPOCHA(NLCED,1)
221 : MUT = UMALPH*MPMTC.VPOCHA(NLCEG,1) +
222 : & ALPHA*MPMTC.VPOCHA(NLCED,1)
223 : c------------------------------------------
224 : DKADX=MPGRKE.VPOCHA(NLCF,1)
225 : DKADY=MPGRKE.VPOCHA(NLCF,2)
226 : DKADZ=MPGRKE.VPOCHA(NLCF,3)
227 : DEPDX=MPGRKE.VPOCHA(NLCF,4)
228 : DEPDY=MPGRKE.VPOCHA(NLCF,5)
229 : DEPDZ=MPGRKE.VPOCHA(NLCF,6)
230 : * Contribution pour l'espèce IESP
231 : MPFLUX.VPOCHA(NLCF,IDIM+2+NESP)=(MU+(MUT/TSIGK))*
232 : $ ((DKADX*CNX)+(DKADY*CNY)+(DKADZ*CNZ))*SURF*(-1.D0)
233 : * Contribution pour l'énergie totale
234 : MPFLUX.VPOCHA(NLCF,IDIM+3+NESP)=(MU+(MUT/TSIGE))*
235 : $ ((DEPDX*CNX)+(DEPDY*CNY)+(DEPDZ*CNZ))*SURF*(-1.D0)
236 : * Le pas de temps
237 : CELL=(4.0D0*(MU+MUT))/DIAM2
238 : IF(CELL .GT. UNSDT)THEN
239 : UNSDT=CELL
240 : ENDIF
241 : C---------------------------------------------------
242 : ENDDO
243 : C
244 : C
245 : DT = 1.0D0 / (UNSDT+XPETIT)
246 : C
247 : SEGDES MELEFL
248 : SEGDES MELEMF
249 : SEGDES MELEMC
250 : SEGDES MPSURF
251 : SEGDES MPNORM
252 : SEGDES MPDIAM
253 : SEGDES MLCENT
254 : C
255 : SEGDES MPMUC
256 : SEGDES MPMTC
257 : SEGDES MPGRKE
258 : SEGDES MPFLUX
259 : C
260 : C
261 : 9999 CONTINUE
262 : RETURN
263 : END
264 :
265 :
266 :
267 :
268 :
269 :
270 :
271 :
272 :
© Cast3M 2003 - All rights reserved.
Disclaimer