1 : C FCOLCP SOURCE CHAT 05/01/12 23:57:19 5004
2 : SUBROUTINE FCOLCP(NESP,
3 : & GAMG,ROG,PG,UNG,UTG,
4 : & GAMD,ROD,PD,UND,UTD,
5 : & YG,YD,FLU,
6 : & CELLT,
7 : & LOGNC,LOGAN,MESERR)
8 : C
9 : C*** N.B. Cela est une modification non optimisée de la subroutine FCOLTE:
10 : C
11 : C OU
12 : C
13 : C NSCA = 0
14 : C NORDP1 = 1
15 : C XST = 0.0D0
16 : C RTOTG = RTOTD = 1.0D0
17 : C ACVTOG(1) = RTOTG / (GAMG - 1.0D0)
18 : C ACVTOD(1) = RTOTG / (GAMD - 1.0D0)
19 : C TG = PG / ROG
20 : C TD = PD / ROD
21 : C ETHERG = ACVTOG(1) * RTOTG
22 : C ETHERD = ACVTOD(1) * RTOTD
23 : C
24 : C
25 : C SUBROUTINE FCOLTE(NESP,NSCA,NORDP1,
26 : C & XST,
27 : C & GAMG,RTOTG,ACVTOG,ROG,PG,TG,UNG,UTG,ETHERG,
28 : C & GAMD,RTOTD,ACVTOD,ROD,PD,TD,UND,UTD,ETHERD,
29 : C & YG,YD,SCAG,SCAD,
30 : C & FLU,CELLT,
31 : C & LOGNC,MESERR,LOGAN)
32 : C************************************************************************
33 : C
34 : C PROJET : CASTEM 2000
35 : C
36 : C NOM : FCOLTE
37 : C
38 : C DESCRIPTION : Formulation Volumes Finis pour les Equations
39 : C d'Euler Multi-Especes relatives à un melange
40 : C de gaz parfaits.
41 : C
42 : C Calcul du flux aux interfaces avec la methode
43 : C de Colella-Glaz, avec les chocs qui remplacent
44 : C les ondes de rarefaction + entropy fix.
45 : C
46 : C (voir:
47 : C 1) BECCANTINI, Thèse de doctorat)
48 : C
49 : C LANGUAGE : FORTRAN 77
50 : C
51 : C AUTEUR : A. BECCANTINI DRN/DMT/SEMT/LTMF
52 : C
53 : C************************************************************************
54 : C
55 : C APPELES : RACCOL
56 : C
57 : C************************************************************************
58 : C
59 : C**** Entrées:
60 : C
61 : C NESP = nombre d'especes considérées dans les Equations
62 : C d'Euler
63 : C
64 : C NSCA = nombre de scalaires passifs a transporter
65 : C
66 : C NORDP1 = (ordre des polynoms cv(T)) + 1
67 : C
68 : C XST = droite characteristique ou je cherche la solution
69 : C
70 : C GAMG, GAMD = les "gamma" du gaz (gauche et droite)
71 : C
72 : C RTOTG, RTOTD = les constantes des gaz
73 : C
74 : C ACVTOG, ACVTOD = les sommes des coefficients des cv,i
75 : C ACVTOT(j) = \sum_{i=1,nesp} Y_i*ACV_{i,j}
76 : C j = 1 , NORDP1
77 : C
78 : C ROG, ROD = les densités
79 : C
80 : C PG, PD = les pressions
81 : C
82 : C TG, TD = les temperatures
83 : C
84 : C UNG, UND = vitesses normales
85 : C
86 : C UTG, UTD = vitesses tangentielles
87 : C
88 : C ETHERG, ETHERD = les energies
89 : C
90 : C YG, YD = tables des fractiones massiques
91 : C
92 : C SCAG, SCAD = tables des scalaires passifs
93 : C
94 : C**** Sorties:
95 : C
96 : C FLU = table du flux a l'interface dans le repaire
97 : C (n,t), i.e.
98 : C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
99 : C rho*un*y1, ...)
100 : C
101 : C CELLT = condition de stabilité, i.e.
102 : C
103 : C dT/diamax < cellt
104 : C
105 : C LOGNC = si TRUE, la methode de newton pour le calcul
106 : C des etats intermediares ne converges pas
107 : C
108 : C MESERR = message d'erreuer
109 : C
110 : C LOGAN = si TRUE, une anomalie a été detectée
111 : C
112 : C************************************************************************
113 : C
114 : C HISTORIQUE (Anomalies et modifications éventuelles)
115 : C
116 : C HISTORIQUE : Créé le 08.02.00
117 : C
118 : C 21.02.00 transport de scalaires passifs
119 : C
120 : C************************************************************************
121 : C
122 : C N.B.: Toutes les variables sont DECLAREES
123 : C
124 : C
125 : IMPLICIT INTEGER(I-N)
126 : INTEGER NESP, NORDP1, I1, NSCA
127 : REAL*8 GAMG,RTOTG,ACVTOG(1),ROG,PG,TG,UNG,UTG,ETHERG
128 : & ,GAMD,RTOTD,ACVTOD(1),ROD,PD,TD,UND,UTD,ETHERD
129 : & ,YG(*),YD(*),SCAG(1),SCAD(1)
130 : & ,FLU(*),CELLT
131 : & ,PSROG,PSROD,ASONG,ASOND,PVAL,UVAL
132 : & ,HTHEG,HTHED
133 : & ,ROGS,PGS,TGS,ETHEGS,UGS,GAMGS
134 : & ,RODS,PDS,TDS,ETHEDS,UDS,GAMDS
135 : & ,XST,UNGXST,UNDXST,UGSXST,UDSXST,DGSXST,DDSXST,DGS,DDS
136 : & ,ROU,ECFLUX,XLIMG,XLIMD,ULIMG,ULIMD
137 : & ,UNGMCG,UNDPCD,PFLUX,ROFLUX,UNFLUX,UTFLUX,HTFLUX
138 : & ,FC1,FC2,FPOIS1,FPOIS2,FPOIS3,FPOIS4,SIG12,SIG45,SIG3
139 : & ,FPOISG,FPOISD,EPSU
140 : & ,UGMCGS,UDPCDS,VITG,VITGS,VITD,VITDS,SIGCHO,SIG123,SIG456
141 : & ,SIG56
142 : & ,FUNT, TFLUX, RTOTFL, ETHFLU
143 : CHARACTER*(40) MESERR
144 : LOGICAL LOGAN, LOGNC
145 : C
146 : C**** Definition pour le passage "thermally-perfect"
146 : C -> "calorically perfect"
147 : C
148 : NSCA = 0
149 : NORDP1 = 1
150 : XST = 0.0D0
151 : RTOTG = 1.0D0
152 : RTOTD = 1.0D0
153 : ACVTOG(1) = RTOTG / (GAMG - 1.0D0)
154 : ACVTOD(1) = RTOTG / (GAMD - 1.0D0)
155 : TG = PG / (ROG * RTOTG)
156 : TD = PD / (ROD * RTOTD)
157 : ETHERG = ACVTOG(1) * TG
158 : ETHERD = ACVTOD(1) * TD
159 : SCAG(1) = 0.0D0
160 : SCAD(1) = 0.0D0
161 : C
162 : C**** YG, YD, FLU
163 : C
164 : C Dans le cas Euler monoespece, on doit
165 : C avoir :
166 : C YG(1) = YD(1) = 0.0D0
167 : C
168 : PSROG = PG / ROG
169 : PSROD = PD / ROD
170 : ASONG = SQRT(GAMG * PSROG)
171 : ASOND = SQRT(GAMD * PSROD)
172 : C
173 : C**** Valeurs de référence
174 : C
175 : PVAL = MAX(PG,PD)
176 : UVAL = MAX(ASONG,ASOND)
177 : EPSU = UVAL*1.0D-6
178 : C
179 : C**** États intermédiaires
180 : C
181 : HTHEG = ETHERG + PSROG
182 : HTHED = ETHERD + PSROD
183 : XLIMG = PSROG / (2.0D0 * HTHEG - PSROG)
184 : XLIMD = PSROD / (2.0D0 * HTHED - PSROD)
185 : ULIMG = UNG + SQRT( 2.0D0 * HTHEG *
186 : & (1.0D0 - XLIMG) / (1.0D0 + XLIMG))
187 : ULIMD = UND - SQRT( 2.0D0 * HTHED *
188 : & (1.0D0 - XLIMD) / (1.0D0 + XLIMD))
189 : UNGMCG = UNG - ASONG
190 : UNDPCD = UND + ASOND
191 : IF(ULIMD .GE. ULIMG)THEN
192 : C
193 : C******* Vide
194 : C
195 : C N.B.: la surface de contact n'existe pas: on mit le vide
196 : C
197 : IF((ULIMG .LE. 0.0D0) .AND. (ULIMD .GE. 0.0D0))THEN
198 : FLU(1) = 0.0D0
199 : FLU(2) = 0.0D0
200 : FLU(3) = 0.0D0
201 : FLU(4) = 0.0D0
202 : DO I1 = 1, NESP, 1
203 : FLU(I1+4) = 0.0D0
204 : ENDDO
205 : DO I1 = 1, NSCA, 1
206 : FLU(I1+4+NESP) = 0.0D0
207 : ENDDO
208 : C
209 : ELSE
210 : ROGS = XLIMG * ROG
211 : RODS = XLIMD * ROD
212 : C
213 : C********** Les functions pois
214 : C
215 : FC1 = (ULIMG - XST)/(ULIMG - UNGMCG)
216 : SIG12 = 0.5D0 * (SIGN(1.0D0,FC1)+1.0D0)
217 : FC2 = 0.5D0*(FC1+1.0D0-ABS(FC1-1.0D0))
218 : FPOIS1 = SIG12*FC2
219 : C
220 : FC1 = (XST - ULIMD)/(UNDPCD - ULIMD)
221 : SIG45 = 0.5D0 * (SIGN(1.0D0,FC1)+1.0D0)
222 : FC2 = 0.5D0*(FC1+1.0D0-ABS(FC1-1.0D0))
223 : FPOIS4 = SIG45*FC2
224 : C
225 : SIG3 = 1.0D0 - (SIG12+SIG45)
226 : FC1 = (ULIMD - XST)/(ULIMD - ULIMG + EPSU)
227 : FC1 = 0.5D0*(FC1+1.0D0-ABS(FC1-1.0D0))
228 : FC2 = 1.0D0 - FC1
229 : FPOIS2 = SIG12 * (1.0D0 - FPOIS1) + SIG3 * FC1
230 : FPOIS3 = SIG45 * (1.0D0 - FPOIS4) + SIG3 * FC2
231 : C
232 : C********** La solution pour le calcul du flux
233 : C
234 : ROFLUX = ROG * FPOIS1 + ROGS * FPOIS2
235 : & + RODS * FPOIS3 + ROD * FPOIS4
236 : TFLUX = TG * FPOIS1 + TD * FPOIS4
237 : UNFLUX = UNG * FPOIS1 + ULIMG * FPOIS2
238 : & + ULIMD * FPOIS3 + UND * FPOIS4
239 : C
240 : FPOISG = FPOIS1 + FPOIS2
241 : FPOISD = FPOIS3 + FPOIS4
242 : RTOTFL = RTOTG * FPOISG + RTOTD * FPOISD
243 : PFLUX = ROFLUX * RTOTFL * TFLUX
244 : UTFLUX = UTG * FPOISG + UTD * FPOISD
245 : ECFLUX = 0.5D0 * (UNFLUX * UNFLUX + UTFLUX * UTFLUX)
246 : C
247 : FUNT = 1.0D0
248 : ETHFLU = 0.0D0
249 : DO I1 = 1, NORDP1, 1
250 : FUNT = FUNT * TFLUX
251 : ETHFLU = ETHFLU + (FPOISG * ACVTOG(I1) * FUNT / I1) +
252 : & (FPOISD * ACVTOD(I1) * FUNT / I1)
253 : ENDDO
254 : HTFLUX = ETHFLU + PFLUX / ROFLUX
255 : C
256 : ROU = ROFLUX * UNFLUX
257 : FLU(1) = ROU
258 : FLU(2) = ROU * UNFLUX + PFLUX
259 : FLU(3) = ROU * UTFLUX
260 : FLU(4) = ROU * (HTFLUX + ECFLUX)
261 : DO I1 = 1, NESP, 1
262 : FLU(I1+4) = ROU *
263 : & (FPOISG * YG(I1) + FPOISD * YD(I1))
264 : ENDDO
265 : DO I1 = 1, NSCA, 1
266 : FLU(I1+4+NESP) = ROU *
267 : & (FPOISG * SCAG(I1) + FPOISD * SCAD(I1))
268 : ENDDO
269 : ENDIF
270 : CELLT = 1.0D0/MAX(ABS(UNGMCG),ABS(UNDPCD))
271 : C
272 : C
273 : C******* Test EF
274 : C
275 : C write(4,'(5D14.6)') XST, FPOIS1, FPOIS2, FPOIS3, FPOIS4
276 : C
277 : ELSE
278 : UNGXST=UNG-XST
279 : UNDXST=UND-XST
279 : CALL RACCOL(PVAL,UVAL,NORDP1,
280 : & RTOTG,ACVTOG,RTOTD,ACVTOD,
281 : & ROG,TG,ETHERG,UNGXST,PSROG,GAMG,ASONG,
282 : & ROD,TD,ETHERD,UNDXST,PSROD,GAMD,ASOND,
283 : & ROGS,PGS,TGS,ETHEGS,UGSXST,DGSXST,GAMGS,
284 : & RODS,PDS,TDS,ETHEDS,UDSXST,DDSXST,GAMDS,
285 : & LOGNC)
286 : C
287 : IF(LOGNC)THEN
288 : MESERR = 'Subroutine RACCOL.ESO, Newton-Raphson.'
289 : C
290 : C******* Message de warning; mais on s'arrete pas
291 : C
292 : ELSEIF((DGSXST .GT. UGSXST) .OR. (UDSXST .GT. DDSXST))THEN
293 : LOGAN = .TRUE.
294 : MESERR = 'Subroutine RACCOL.ESO'
295 : GOTO 9999
296 : ENDIF
297 : UGS=UGSXST+XST
298 : UDS=UDSXST+XST
299 : DGS=DGSXST+XST
300 : DDS=DDSXST+XST
301 : UGMCGS=UGS-SQRT(GAMGS*PGS/ROGS)
302 : UDPCDS=UDS+SQRT(GAMDS*PDS/RODS)
303 :
304 : C
305 : C******* Tests
306 : C
307 : C write(*,*) ROGS,PGS,TGS,ETHEGS,UGS,GAMGS
308 : C write(*,*) RODS,PDS,TDS,ETHEDS,UDS,GAMDS
309 : C
310 : C
311 : C******* Les functions pois
312 : C
313 : C Les vitesses characteristiques des ondes
314 : C
315 : VITG = 0.5D0*(DGS+UNGMCG-ABS(DGS-UNGMCG))
316 : VITGS = 0.5D0*(DGS+UGMCGS+ABS(DGS-UGMCGS))
317 : VITD = 0.5D0*(DDS+UNDPCD+ABS(DDS-UNDPCD))
318 : VITDS = 0.5D0*(DDS+UDPCDS-ABS(DDS-UDPCDS))
319 : C
320 : C******* ZONES 1,2,3
321 : C
322 : FC1 = (VITGS-XST)/(ABS(VITGS - VITG)+EPSU)
323 : SIG12 = 0.5D0*(SIGN(1.0D0,FC1)+1.0D0)
324 : FC2 = 0.5D0*(FC1 + 1.0D0 - ABS(FC1 - 1.0D0))
325 : SIGCHO = 0.5D0 *(SIGN(1.0D0,(VITG -VITGS+EPSU))+1.0D0)
326 : SIG123 = 0.5D0 *(SIGN(1.0D0,((0.5D0*(UGS+UDS))-XST))+1.0D0)
327 : C
328 : FPOIS2 = (1.0D0 - SIGCHO)*(1.0D0 - FC2*SIG12)*SIG123 +
329 : & SIGCHO * (SIG123 - SIG12)
330 : FPOIS1 = (1.0D0 -FPOIS2)*SIG123
331 : C
332 : C******** ZONES 4,5,6
333 : C
334 : SIG456 = 1.0D0 - SIG123
335 : FC1 = (XST - VITDS)/(ABS(VITD-VITDS)+EPSU)
336 : SIG56 = 0.5D0*(SIGN(1.0D0,FC1)+1.0D0)
337 : FC2 = 0.5D0 * (FC1 + 1.0D0 - ABS(FC1 - 1.0D0))
338 : SIGCHO=0.5D0*(SIGN(1.0D0,VITDS-VITD+EPSU)+1.0D0)
339 : C
340 : FPOIS3= (1.0D0 - SIGCHO)*(1.0D0 - FC2*SIG56)*SIG456 +
341 : & SIGCHO * (SIG456 - SIG56)
342 : FPOIS4 = (1.0D0 - FPOIS3) * SIG456
343 : C
344 : C******* La solution pour le calcul du flux
345 : C (ro, T, w lineaire sur les detantes
346 : C p = ro R T
347 : C h = sum_i ACVTO(i)/i T^i
348 : C
349 : C
350 : ROFLUX = ROG * FPOIS1 + ROGS * FPOIS2
351 : & + RODS * FPOIS3 + ROD * FPOIS4
352 : TFLUX = TG * FPOIS1 + TGS * FPOIS2
353 : & + TDS * FPOIS3 + TD * FPOIS4
354 : UNFLUX = UNG * FPOIS1 + UGS * FPOIS2
355 : & + UDS * FPOIS3 + UND * FPOIS4
356 : C
357 : FPOISG = FPOIS1 + FPOIS2
358 : FPOISD = FPOIS3 + FPOIS4
359 : RTOTFL = RTOTG * FPOISG + RTOTD * FPOISD
360 : PFLUX = ROFLUX * RTOTFL * TFLUX
361 : UTFLUX = UTG * FPOISG + UTD * FPOISD
362 : ECFLUX = 0.5D0 * (UNFLUX * UNFLUX + UTFLUX * UTFLUX)
363 : C
364 : FUNT = 1.0D0
365 : ETHFLU = 0.0D0
366 : DO I1 = 1, NORDP1, 1
367 : FUNT = FUNT * TFLUX
368 : ETHFLU = ETHFLU + (FPOISG * ACVTOG(I1) * FUNT / I1) +
369 : & (FPOISD * ACVTOD(I1) * FUNT / I1)
370 : ENDDO
371 : HTFLUX = ETHFLU + PFLUX / ROFLUX
372 : C
373 : ROU = ROFLUX * UNFLUX
374 : FLU(1) = ROU
375 : FLU(2) = ROU * UNFLUX + PFLUX
376 : FLU(3) = ROU * UTFLUX
377 : FLU(4) = ROU * (HTFLUX + ECFLUX)
378 : DO I1 = 1, NESP, 1
379 : FLU(I1+4) = ROU *
380 : & (FPOISG * YG(I1) + FPOISD * YD(I1))
381 : ENDDO
382 : DO I1 = 1, NSCA, 1
383 : FLU(I1+4+NESP) = ROU *
384 : & (FPOISG * SCAG(I1) + FPOISD * SCAD(I1))
385 : ENDDO
386 : C
387 : CELLT = 1.0/MAX(ABS(VITG),ABS(VITD))
388 : C
389 : C******* Tests: les vitesses
390 : C
391 : C write(3,'(A8,D16.8)') 'VITG =', VITG
392 : C write(3,'(A8,D16.8)') 'VITGS =', VITGS
393 : C write(3,'(A8,2D16.8)') 'VITS =', UGS, UDS
394 : C write(3,'(A8,D16.8)') 'VITDS =', VITDS
395 : C write(3,'(A8,D16.8)') 'VITD =', VITD
396 : C
397 : C
398 : C******* Test EF
399 : C
400 : C write(4,'(5D14.6)') XST, FPOIS1, FPOIS2, FPOIS3, FPOIS4
401 : C
402 : ENDIF
403 : C
404 : 9999 CONTINUE
405 : RETURN
406 : END
407 :
408 :
409 :
410 :
411 :
412 :
413 :
414 :
415 :
© Cast3M 2003 - All rights reserved.
Disclaimer