Download kfmm1.eso
Back to the list
1 : C KFMM1 SOURCE CHAT 05/01/13 00:55:42 5004
2 : SUBROUTINE KFMM1(IRN,IGN,IRETN,IGAMN,
3 : & ICHPSU,ICHPDI,ICHPVO,INORM,
4 : & MELEMC,MELEMF,MELEFE,DCFL,
5 : & MELLIM,ICHLIM,
6 : & ICHRES)
7 : C************************************************************************
8 : C
9 : C PROJET : CASTEM 2000
10 : C
11 : C NOM : KFMM1
12 : C
13 : C DESCRIPTION : Voir KFMM
14 : C
14 : C Cas deux dimensions, gaz "calorically perfect"
15 : C
16 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
17 : C
18 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF
19 : C
20 : C************************************************************************
21 : C ENTREES
22 : C
23 : C
24 : C 1) Pointeurs des CHPOINTs
25 : C
26 : C IRN : CHPOINT 'CENTRE' contenant la masse volumique
27 : C
28 : C IGN : CHPOINT 'CENTRE' contenant la q.d.m.
29 : C
30 : C IRETN : CHPOINT 'CENTRE' contenant l'energie totale
31 : C
32 : C IGAMN : CHPOINT 'CENTRE' contenant le gamma
33 : C
34 : C ICHLIM : CHPOINT conditions aux bords
35 : C
36 : C 3) Pointeurs de CHPOINTs de la table DOMAINE
37 : C
38 : C ICHPSU : CHPOINT "FACE" contenant la surface des faces
39 : C
40 : C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum
41 : C de chaque element
42 : C
43 : C ICHPVO : CHPOINT "CENTRE" contenant le volume
44 : C de chaque element
45 : C
46 : C INORM : CHPOINT "CENTRE" contenant le normales aux faces
47 : C
48 : C
49 : C 4) Pointeurs de MELEME de la table DOMAINE
50 : C
51 : C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
52 : C
53 : C MELEMF : MELEME 'FACE' du SPG des FACES
54 : C
55 : C MELEFE : MELEME 'FACEL' du connectivité FACES -> ELEM
56 : C
57 : C MELLIM : MELEME SPG conditions aux bords
58 : C
59 : C 5)
60 : C
61 : C DCFL = le double de la CFL
62 : C
63 : C
64 : C SORTIES
65 : C
66 : C ICHRES : resultat
67 : C
68 : C************************************************************************
69 : C
70 : C HISTORIQUE (Anomalies et modifications éventuelles)
71 : C
72 : C HISTORIQUE : Avril 2002 : creation
73 : C Janvier 2003: implementation de condition aux limites
74 : C
75 : C************************************************************************
76 : C
77 : C
78 : C N.B.: On suppose qu'on a déjà controllé RO, P > 0
79 : C GAMMA \in (1,3)
80 : C Y \in (0,1)
81 : C Si non il faut le faire!!!
82 : C
83 : C************************************************************************
84 : C
85 : C
86 : C**** Variables de COOPTIO
87 : C
88 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
89 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
90 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
91 : C & ,IECHO, IIMPI, IOSPI
92 : C & ,IDIM
93 : C & ,MCOORD
94 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
95 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
96 : C & ,NORINC,NORVAL,NORIND,NORVAD
97 : C & ,NUCROU, IPSAUV, IFICLE, IPREFI
98 : C
99 : IMPLICIT INTEGER(I-N)
100 : INTEGER IRN,IGN,IRETN,IGAMN,ICHPSU,ICHPDI,ICHPVO,INORM
101 : & ,MELEMC,MELEMF,MELEFE,ICHRES,ISPG1,ISPG2,NFAC
102 : & ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCF1, NLCEG, NLCED
103 : & ,ICEN,NCEN,ICHLIM,MELLIM,NLLIM
104 : REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, REC, PG, VOLG, DIAMG
105 : & , ROD, RUXD, RUYD, UND, RETD, GAMD, PD, VOLD, DIAMD
106 : & , CNX, CNY, DCFL
107 : & , SURF, SIGMAF
108 : & , UNSDT, UNSDTF, UXD, UYD
109 : CHARACTER*8 TYPE
110 : C
111 : C**** LES INCLUDES
112 : C
113 : -INC CCOPTIO
114 : -INC SMCHPOI
115 : POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL
116 : & , MPORES.MPOVAL, MPODTI.MPOVAL, MPOVOL.MPOVAL
117 : & , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL
118 : & , MPNORM.MPOVAL, MPLIM.MPOVAL
119 : -INC SMELEME
120 : -INC SMLMOTS
121 : -INC SMLENTI
122 : POINTEUR MLELIM.MLENTI
123 : C
124 : C**** Initialisation des MLENTI des conditions aux limites
125 : C
126 : C
126 : CALL KRIPAD(MELLIM,MLELIM)
127 : C SEGINI MLELIM
128 : C
129 : C**** Initialisation des MELEMEs
130 : C
131 : C 'CENTRE', 'FACEL'
132 : C
133 : IPT2 = MELEFE
134 : SEGACT IPT2
135 : NFAC = IPT2.NUM(/2)
136 : C
137 : C**** KRIPAD pour la correspondance global/local de centre
138 : C
138 : CALL KRIPAD(MELEMC,MLENT1)
139 : C
140 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
141 : C
142 : C Si i est le numero global d'un noeud de ICEN,
143 : C MLENT1.LECT(i) contient sa position, i.e.
144 : C
145 : C I = numero global du noeud centre
146 : C MLENT1.LECT(i) = numero local du noeud centre
147 : C
148 : C MLENT1 déjà activé, i.e.
149 : C
150 : C SEGINI MLENT1
151 : C
152 : C
153 : C**** KRIPAD pour la correspondance global/local de 'FACE'
154 : C
154 : CALL KRIPAD(MELEMF,MLENT2)
155 : C SEGINI MLENT2
156 : C
157 : C
158 : C**** CHPOINTs de la table DOMAINE
159 : C
159 : CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
159 : CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
159 : CALL LICHT(ICHPVO,MPOVOL,TYPE,IGEOMC)
159 : CALL LICHT(INORM,MPNORM,TYPE,IGEOMC)
160 : C
161 : C**** LICHT active les MPOVALs en *MOD
162 : C
163 : C i.e.
164 : C
165 : C SEGACT MPOVSU*MOD
166 : C SEGACT MPOVOL*MOD
167 : C SEGACT MPOVDI*MOD
168 : C SEGACT MPNORM*MOD
169 : C
169 : CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC)
170 : C SEGACT MPORES*MOD
171 : C
172 : C MPODTI initialisé a zero; MPODTI = 1 / DT
173 : C
174 : SEGINI, MPODTI=MPORES
175 : C
175 : CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
175 : CALL LICHT(IGN,MPGN,TYPE,IGEOMC)
175 : CALL LICHT(IRETN,MPRETN,TYPE,IGEOMC)
175 : CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC)
175 : CALL LICHT(ICHLIM,MPLIM,TYPE,MELLIM)
176 : C
177 : C SEGACT MPRN*MOD
178 : C SEGACT MPGN*MOD
179 : C SEGACT MPRETN*MOD
180 : C SEGACT MPGAMN*MOD
181 : C SEGACT MPLIM*MOD
182 : C
183 : IF(IGEOMF .NE. MELEMF)THEN
184 : WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso'
185 : WRITE(IOIMP,*) 'Probleme de SPG'
186 : C 21 2
187 : C Données incompatibles
187 : CALL ERREUR(21)
188 : GOTO 9999
189 : ENDIF
190 : IF(IGEOMC .NE. MELEMC)THEN
191 : WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso'
192 : WRITE(IOIMP,*) 'Probleme de SPG'
193 : C 21 2
194 : C Données incompatibles
194 : CALL ERREUR(21)
195 : GOTO 9999
196 : ENDIF
197 : C
198 : C
199 : C**** BOUCLE SUR FACEL pour le calcul du FLUX
200 : C
201 : DO NLCF = 1, NFAC
202 : C
203 : C******* NLCF = numero local du centre de facel
204 : C NGCF = numero global du centre de facel
205 : C NLCF1 = numero local du centre de face
206 : C NGCEG = numero global du centre ELT "gauche"
207 : C NLCEG = numero local du centre ELT "gauche"
208 : C NGCED = numero global du centre ELT "droite"
209 : C NLCED = numero local du centre ELT "droite"
210 : C
211 : NGCEG = IPT2.NUM(1,NLCF)
212 : NGCED = IPT2.NUM(3,NLCF)
213 : NGCF = IPT2.NUM(2,NLCF)
214 : NLCF1 = MLENT2.LECT(NGCF)
215 : NLCEG = MLENT1.LECT(NGCEG)
216 : NLCED = MLENT1.LECT(NGCED)
217 : C
218 : C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
219 : C
220 : IF(NLCF .NE. NLCF1)THEN
221 : WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso'
222 : WRITE(IOIMP,*) 'Probleme de SPG'
223 : C 21 2
224 : C Données incompatibles
224 : CALL ERREUR(21)
225 : GOTO 9999
226 : ENDIF
227 : C
228 : CNX = MPNORM.VPOCHA(NLCF,1)
229 : CNY = MPNORM.VPOCHA(NLCF,2)
230 : C
231 : C******* Recuperation des Etats "gauche" et "droite"
232 : C
233 : C
234 : ROG = MPRN.VPOCHA(NLCEG,1)
235 : RUXG = MPGN.VPOCHA(NLCEG,1)
236 : RUYG = MPGN.VPOCHA(NLCEG,2)
237 : UNG = (RUXG * CNX) + (RUYG * CNY)
238 : UNG = UNG / ROG
239 : RETG = MPRETN.VPOCHA(NLCEG,1)
240 : GAMG = MPGAMN.VPOCHA(NLCEG,1)
241 : REC = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
242 : REC = REC / ROG
243 : PG = (GAMG - 1.0D0) * (RETG - REC)
244 : VOLG = MPOVOL.VPOCHA(NLCEG,1)
245 : DIAMG = MPOVDI.VPOCHA(NLCEG,1)
246 : C
247 : ROD = MPRN.VPOCHA(NLCED,1)
248 : RUXD = MPGN.VPOCHA(NLCED,1)
249 : RUYD = MPGN.VPOCHA(NLCED,2)
250 : RETD = MPRETN.VPOCHA(NLCED,1)
251 : GAMD = MPGAMN.VPOCHA(NLCED,1)
252 : REC = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD))
253 : REC = REC / ROD
254 : PD = (GAMD - 1.0D0) * (RETD - REC)
255 : VOLD = MPOVOL.VPOCHA(NLCED,1)
256 : DIAMD = MPOVDI.VPOCHA(NLCED,1)
257 : IF(NLCEG .NE. NLCED)THEN
258 : UND = (RUXD * CNX) + (RUYD * CNY)
259 : UND = UND / ROD
260 : ELSE
261 : C Murs au condition aux limite
262 : NLLIM=MLELIM.LECT(NGCF)
263 : IF(NLLIM .EQ.0)THEN
264 : C Mur
265 : UND = -1.0D0 * UNG
266 : ELSE
267 : ROD = MPLIM.VPOCHA(NLLIM,1)
268 : UXD = MPLIM.VPOCHA(NLLIM,2)
269 : UYD = MPLIM.VPOCHA(NLLIM,3)
270 : PD = MPLIM.VPOCHA(NLLIM,4)
271 : GAMD = GAMG
272 : UND = (UXD * CNX) + (UYD * CNY)
273 : RETD = ((1.0D0/(GAMD - 1.0D0))*PD)+
274 : & (0.5D0*ROD*((UXD*UXD)+(UYD*UYD)))
275 : VOLD = VOLG
276 : ENDIF
277 : ENDIF
278 : C
279 : SURF = MPOVSU.VPOCHA(NLCF,1)
280 : SIGMAF = (0.5D0 * (GAMG + GAMD)) * (PG + PD) / (ROG + ROD)
281 : SIGMAF = SIGMAF ** 0.5D0
282 : SIGMAF = (0.5D0 * (ABS(UNG) + ABS(UND))) + SIGMAF
283 : C
284 : C******* On a defini (ROg,ROUNg,ROUTg,Pg,(Yg)), (ROd,ROUNd,ROUTd,Pd,(Yd))
285 : C et on a déjà verifié ROg, ROd, Pg, Pd > 0 et 0