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