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