1 : C TLAP12 SOURCE CHAT 05/01/13 03:41:38 5004
2 : SUBROUTINE TLAP12(NESP,IMUC,IMTC,TSIGK,TSIGE,IGRKEP,
3 : & MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM,
4 : & ICHFLU,DT)
5 : C***********************************************************************
6 : C NOM : TLAP12
7 : C DESCRIPTION : Calcul des flux diffusifs pour k-\eps equations de
8 : C turbulence (2D)
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,XFMXG,YFMYG,DRG
79 : & ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH
80 : & ,XF,YF
81 : & ,CNX, CNY, ORIENT, DIAM, DIAM2, CELL, SURF
82 : INTEGER NESP
83 : REAL*8 DKADX,DKADY,DEPDX,DEPDY
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 : C
145 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
146 : XG = MCOORD.XCOOR(ICOORX)
147 : YG = MCOORD.XCOOR(ICOORX+1)
148 : XFMXG = XF - XG
149 : YFMYG = YF - YG
150 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG))
151 : C
152 : ICOORX = ((IDIM + 1) * (NGCED - 1))+1
153 : XD = MCOORD.XCOOR(ICOORX)
154 : YD = MCOORD.XCOOR(ICOORX+1)
155 : XFMXD = XF - XD
156 : YFMYD = YF - YD
157 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD))
158 : C
159 : C********** F=G -> DRG = 0 -> ALPHA = 0
160 : C
161 : ALPHA=DRG/(DRG+DRD)
162 : UMALPH= 1.0D0 - ALPHA
163 : C
164 : C********** Les valeurs à l'interface
165 : C
166 : C DRG=0 -> F=G
167 : C
168 : ELSE
169 : C
170 : C********** MURS
171 : C
172 : C Etat a gauche = Etat droite
173 : C
174 : ALPHA=0.0D0
175 : UMALPH=1.0D0
176 : C
177 : C********** Parametres geometriques
178 : C
179 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1
180 : XF = MCOORD.XCOOR(ICOORX)
181 : YF = MCOORD.XCOOR(ICOORX+1)
182 : C
183 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
184 : XG = MCOORD.XCOOR(ICOORX)
185 : YG = MCOORD.XCOOR(ICOORX+1)
186 : XFMXG = XF - XG
187 : YFMYG = YF - YG
188 : ENDIF
189 : C
190 : C******* On calcule le sign du pruduit scalare
191 : C (Normales de Castem) * (vecteur "gauche" -> "centre")
192 : C
193 : CNX = MPNORM.VPOCHA(NLCF,1)
194 : CNY = MPNORM.VPOCHA(NLCF,2)
195 : ORIENT = CNX * XFMXG + CNY * YFMYG
196 : ORIENT = SIGN(1.0D0,ORIENT)
197 : IF(ORIENT .NE. 1.0D0)THEN
198 : MOTERR(1:40)=
199 : & 'LAPN , subroutine zlap12.eso. '
200 : WRITE(IOIMP,*) MOTERR(1:40)
200 : CALL ERREUR(5)
201 : GOTO 9999
202 : ENDIF
203 : C
204 : C******* Les flux aux interfaces
205 : C
206 : SURF = MPSURF.VPOCHA(NLCF,1)
207 : DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) +
208 : & ALPHA*MPDIAM.VPOCHA(NLCED,1)
209 : DIAM2=DIAM*DIAM
210 : MU = UMALPH*MPMUC.VPOCHA(NLCEG,1) +
211 : & ALPHA*MPMUC.VPOCHA(NLCED,1)
212 : MUT = UMALPH*MPMTC.VPOCHA(NLCEG,1) +
213 : & ALPHA*MPMTC.VPOCHA(NLCED,1)
214 : c------------------------------------------
215 : DKADX=MPGRKE.VPOCHA(NLCF,1)
216 : DKADY=MPGRKE.VPOCHA(NLCF,2)
217 : DEPDX=MPGRKE.VPOCHA(NLCF,3)
218 : DEPDY=MPGRKE.VPOCHA(NLCF,4)
219 : * Contribution pour l'espèce IESP
220 : MPFLUX.VPOCHA(NLCF,IDIM+2+NESP)=
221 : $ (MU+(MUT/TSIGK))*((DKADX*CNX)+(DKADY*CNY))*SURF*(-1.D0)
222 : * Contribution pour l'énergie totale
223 : MPFLUX.VPOCHA(NLCF,IDIM+3+NESP)=
224 : $ (MU+(MUT/TSIGE))*((DEPDX*CNX)+(DEPDY*CNY))*SURF*(-1.D0)
225 : * Le pas de temps (proportional to \rho) - not implemented
226 : CELL=(4.0D0*(MU+MUT))/DIAM2
227 : IF(CELL .GT. UNSDT)THEN
228 : UNSDT=CELL
229 : ENDIF
230 : C---------------------------------------------------
231 : ENDDO
232 : C
233 : C-----------------------------------------------
234 : DT = 1.0D0 / (UNSDT+XPETIT)
235 : C
236 : SEGDES MELEFL
237 : SEGDES MELEMF
238 : SEGDES MELEMC
239 : SEGDES MPSURF
240 : SEGDES MPNORM
241 : SEGDES MPDIAM
242 : SEGDES MLCENT
243 : C
244 : SEGDES MPMUC
245 : SEGDES MPMTC
246 : SEGDES MPGRKE
247 : SEGDES MPFLUX
248 : C
249 : C
250 : 9999 CONTINUE
251 : RETURN
252 : END
253 :
254 :
255 :
256 :
257 :
258 :
259 :
260 :
261 :
© Cast3M 2003 - All rights reserved.
Disclaimer