← Back to branch summary

~ubuntu-branches/ubuntu/precise/code-saturne/precise

« back to all changes in this revision

Viewing changes to src/fuel/fuflux.f90

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-24 00:00:08 UTC
  • mfrom: (6.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20111124000008-2vo99e38267942q5
Tags: 2.1.0-3
Install a missing file

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/fuel/fuflux.f90
1
 
!-------------------------------------------------------------------------------
2
 
 
3
 
!     This file is part of the Code_Saturne Kernel, element of the
4
 
!     Code_Saturne CFD tool.
5
 
 
6
 
!     Copyright (C) 1998-2009 EDF S.A., France
7
 
 
8
 
!     contact: saturne-support@edf.fr
9
 
 
10
 
!     The Code_Saturne Kernel is free software; you can redistribute it
11
 
!     and/or modify it under the terms of the GNU General Public License
12
 
!     as published by the Free Software Foundation; either version 2 of
13
 
!     the License, or (at your option) any later version.
14
 
 
15
 
!     The Code_Saturne Kernel is distributed in the hope that it will be
16
 
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
17
 
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18
 
!     GNU General Public License for more details.
19
 
 
20
 
!     You should have received a copy of the GNU General Public License
21
 
!     along with the Code_Saturne Kernel; if not, write to the
22
 
!     Free Software Foundation, Inc.,
23
 
!     51 Franklin St, Fifth Floor,
24
 
!     Boston, MA  02110-1301  USA
25
 
 
26
 
!-------------------------------------------------------------------------------
27
 
 
28
 
subroutine fuflux &
29
 
!================
30
 
 
31
 
 ( idbia0 , idbra0 ,                                              &
32
 
   ncelet , ncel   ,                                              &
33
 
   rtpa   , propce , volume ,                                     &
34
 
   w1     , w2     , w3     ,                                     &
35
 
   ra     )
36
 
 
37
 
!===============================================================================
38
 
! FONCTION :
39
 
! --------
40
 
 
41
 
! CALCUL DES TERMES DE TRANSFERT DE MASSE ENTRE LA PHASE CONTINUE
42
 
! ET LA PHASE DISPERSEE
43
 
 
44
 
 
45
 
! Arguments
46
 
!__________________.____._____.________________________________________________.
47
 
! name             !type!mode ! role                                           !
48
 
!__________________!____!_____!________________________________________________!
49
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
50
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
51
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
52
 
! ncel             ! i  ! <-- ! number of cells                                !
53
 
! rtpa             ! tr ! <-- ! variables de calcul au centre des              !
54
 
! (ncelet,*)       !    !     !    cellules (instant precedent)                !
55
 
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
56
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
57
 
! w1, w2, w3       ! tr ! --- ! tableaux de travail                            !
58
 
! ra(*)            ! ra ! --- ! main real work array                           !
59
 
!__________________!____!_____!________________________________________________!
60
 
 
61
 
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
62
 
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
63
 
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
64
 
!            --- tableau de travail
65
 
!===============================================================================
66
 
 
67
 
implicit none
68
 
 
69
 
!===============================================================================
70
 
! Common blocks
71
 
!===============================================================================
72
 
 
73
 
include "paramx.h"
74
 
include "numvar.h"
75
 
include "optcal.h"
76
 
include "cstphy.h"
77
 
include "cstnum.h"
78
 
include "entsor.h"
79
 
include "parall.h"
80
 
include "ppppar.h"
81
 
include "ppthch.h"
82
 
include "coincl.h"
83
 
include "cpincl.h"
84
 
include "fuincl.h"
85
 
include "ppincl.h"
86
 
 
87
 
!===============================================================================
88
 
 
89
 
! Arguments
90
 
 
91
 
integer          idbia0 , idbra0
92
 
integer          ncelet , ncel
93
 
 
94
 
double precision rtpa(ncelet,*), propce(ncelet,*)
95
 
double precision w1(ncelet), w2(ncelet), w3(ncelet)
96
 
double precision volume(ncelet)
97
 
double precision ra(*)
98
 
 
99
 
! Local variables
100
 
 
101
 
integer          idebia , idebra
102
 
integer          iel    , iphas  , icla
103
 
integer          ipcrom , ipcte1 , ipcte2 , ipcro2 , ipcdia
104
 
integer          ipcgev , ipcght , ipcyox
105
 
integer          ipcvst,ipcvsl,ipccp,ipchgl
106
 
 
107
 
double precision xng,xnuss
108
 
double precision pparo2 , xdffli , xdfext , xdftot0 , xdftot1
109
 
double precision diacka, xuash
110
 
double precision dcoke , surf
111
 
 
112
 
!===============================================================================
113
 
! 1. INITIALISATIONS ET CALCULS PRELIMINAIRES
114
 
!===============================================================================
115
 
 
116
 
! --- Initialisation memoire
117
 
 
118
 
idebia = idbia0
119
 
idebra = idbra0
120
 
 
121
 
! --- Initialisation des termes de transfert de masse
122
 
 
123
 
do icla = 1, nclafu
124
 
  ipcgev = ipproc(igmeva(icla))
125
 
  ipcght = ipproc(igmhtf(icla))
126
 
  do iel = 1, ncel
127
 
    propce(iel,ipcgev) = zero
128
 
    propce(iel,ipcght) = zero
129
 
  enddo
130
 
enddo
131
 
 
132
 
! --- Initialisation des tableaux de travail
133
 
 
134
 
do iel = 1, ncel
135
 
  w1(iel) = zero
136
 
  w2(iel) = zero
137
 
  w3(iel) = zero
138
 
enddo
139
 
 
140
 
! --- Calcul de la masse volumique du melange gazeux
141
 
 
142
 
iphas  = 1
143
 
ipcrom = ipproc(irom(iphas))
144
 
ipcte1 = ipproc(itemp1)
145
 
ipcyox = ipproc(iym1(io2))
146
 
 
147
 
! --- Numero des grandeurs physiques (voir usclim)
148
 
ipcrom = ipproc(irom(iphas))
149
 
ipcvst = ipproc(ivisct(iphas))
150
 
 
151
 
! --> Terme source pour l'enthalpie du liquide
152
 
 
153
 
do icla = 1, nclafu
154
 
 
155
 
  ipcro2 = ipproc(irom3 (icla))
156
 
  ipcdia = ipproc(idiam3(icla))
157
 
  ipcte2 = ipproc(itemp3(icla))
158
 
  ipcght = ipproc(igmhtf(icla))
159
 
  ipchgl = ipproc(ih1hlf(icla))
160
 
 
161
 
 
162
 
! ---- Contribution aux bilans explicite et implicite
163
 
!        des echanges par diffusion moleculaire
164
 
!        6 Lambda Nu / diam**2 / Rho2 * Rho * (T1-T2)
165
 
! ------ Calcul de lambda dans W1
166
 
 
167
 
  xnuss = 2.d0
168
 
  do iel = 1, ncel
169
 
    if ( ivisls(ihm).gt.0 ) then
170
 
      ipcvsl = ipproc(ivisls(ihm))
171
 
      if ( icp(iphas).gt.0 ) then
172
 
        ipccp   = ipproc(icp(iphas))
173
 
        w1(iel) = propce(iel,ipcvsl) * propce(iel,ipccp)
174
 
      else
175
 
        w1(iel) = propce(iel,ipcvsl) * cp0(iphas)
176
 
      endif
177
 
    else
178
 
      if ( icp(iphas).gt.0 ) then
179
 
        ipccp   = ipproc(icp(iphas))
180
 
        w1(iel) = visls0(ihm) * propce(iel,ipccp)
181
 
      else
182
 
        w1(iel) = visls0(ihm) * cp0(iphas)
183
 
      endif
184
 
    endif
185
 
  enddo
186
 
 
187
 
!----Contribution aux bilans explicite et implicite des
188
 
!    echanges par diffusion moleculaire
189
 
!        6 Lambda Nu / diam**2 / Rho2 * Rho
190
 
!    le diametre est en mm  donc on multiplie par 1.D-3
191
 
!    pour l'avoir en m
192
 
 
193
 
  do iel = 1, ncel
194
 
 
195
 
 
196
 
    if ( rtpa(iel,isca(ing(icla))) .gt. epsifl .and.              &
197
 
           propce(iel,ipcte1).gt. propce(iel,ipcte2) ) then
198
 
 
199
 
 
200
 
!          PROPCE(IEL,IPCHGL) = 6.D0 * W1(IEL) * XNUSS
201
 
!     &                         /((PROPCE(IEL,IPCDIA)*1.D-3)**2)
202
 
!     &                         /PROPCE(IEL,IPCRO2)
203
 
 
204
 
      propce(iel,ipchgl) = w1(iel)*xnuss*rtpa(iel,isca(ing(icla)))&
205
 
                          *pi*propce(iel,ipcdia)*1.d6
206
 
 
207
 
    else
208
 
      propce(iel,ipchgl) =0.d0
209
 
    endif
210
 
 
211
 
  enddo
212
 
 
213
 
enddo
214
 
 
215
 
!===============================================================================
216
 
! 2. TRANSFERTS DE MASSE PAR EVAPORATION
217
 
!===============================================================================
218
 
 
219
 
do icla = 1, nclafu
220
 
 
221
 
  ipcro2 = ipproc(irom3 (icla))
222
 
  ipcdia = ipproc(idiam3(icla))
223
 
  ipcte2 = ipproc(itemp3(icla))
224
 
  ipcgev = ipproc(igmeva(icla))
225
 
  ipchgl = ipproc(ih1hlf(icla))
226
 
 
227
 
! FO & PPl 16/09/05  Version avec parametres en dur, en attendant les
228
 
!                    lectures et inclusion en include
229
 
  tevap1 = 150.d0 + tkelvi
230
 
  tevap2 = 450.d0 + tkelvi
231
 
 
232
 
 
233
 
! PPl 161205 On r�partit le flux interfacial d'enthalpie entre
234
 
!            l'�chauffement de la goutte et l'�vaporation
235
 
 
236
 
  do iel = 1, ncel
237
 
 
238
 
! --- Transfert de masse du a l'evaporation
239
 
 
240
 
    propce(iel,ipcgev) = zero
241
 
 
242
 
!         IF (PROPCE(IEL,IPCTE2)    .GE. TEVAP1 .AND.
243
 
!    &       PROPCE(IEL,IPCTE2)    .LE. TEVAP2 .AND.
244
 
!    &       RTPA(IEL,ISCA(IYFOL(ICLA))) .GE. EPSIFL      ) THEN
245
 
!          PROPCE(IEL,IPCGEV) = PROPCE(IEL,IPPROC(IH1HLF)) /
246
 
!     &    ( RTPA(IEL,ISCA(IYFOL)) * CP2FOL * (TEVAP2-TEVAP1)  /
247
 
!     &      ( ( RTPA(IEL,ISCA(IYFOL)) + RTPA(IEL,ISCA(IFVAP)) )
248
 
!     &    * (1    .d0-FKC))                                     + HRFVAP )
249
 
! PPl 161205
250
 
 
251
 
! PPl 090106 On teste plut�t les diam�tres car d�sormais les
252
 
!        gouttes finissent trop fines
253
 
 
254
 
    if ( propce(iel,ipcte2)    .gt. tevap1 .and.                  &
255
 
         propce(iel,ipcdia)    .gt. dinikf(icla) .and.            &
256
 
         rtpa(iel,isca(iyfol(icla))) .gt. epsifl       ) then
257
 
      propce(iel,ipcgev) = propce(iel,ipchgl)                     &
258
 
                        /( hrfvap + cp2fol*(tevap2-tevap1) )
259
 
    endif
260
 
 
261
 
  enddo
262
 
 
263
 
enddo
264
 
 
265
 
!===============================================================================
266
 
! 3. CALCUL DE RHO_COKE MOYEN
267
 
!    On suppose pour le calcul de la masse volumique du coke que
268
 
!    l'evaporation a lieu a volume de coke constant
269
 
!    Par la suite, on suppose (pour commencer ?) que la
270
 
!    combustion h�t�rog�ne a lieu � volume constant =>
271
 
!    � masse volumique d�croissante
272
 
!===============================================================================
273
 
 
274
 
! --- Initialisation
275
 
 
276
 
rhokf  = rho0fl
277
 
 
278
 
! --- Calcul de la masse volumique moyenne du coke
279
 
 
280
 
!===============================================================================
281
 
! 4. TRANSFERTS DE MASSE PAR COMBUSTION HETEROGENE
282
 
!===============================================================================
283
 
 
284
 
do icla = 1, nclafu
285
 
 
286
 
  ipcro2 = ipproc(irom3 (icla))
287
 
  ipcdia = ipproc(idiam3(icla))
288
 
  ipcte2 = ipproc(itemp3(icla))
289
 
  ipcgev = ipproc(igmeva(icla))
290
 
  ipcght = ipproc(igmhtf(icla))
291
 
 
292
 
  do iel = 1, ncel
293
 
 
294
 
    if ( propce(iel,ipcdia) .le. dinikf(icla) .and.               &
295
 
         propce(iel,ipcdia) .gt. diniin(icla) .and.               &
296
 
         rtpa(iel,isca(iyfol(icla))) .gt. epsifl ) then
297
 
 
298
 
      xng   = rtpa(iel,isca(ing(icla)))*1.d9
299
 
 
300
 
      xuash = xng*(1.d0-xinfol)                                   &
301
 
             *(rho0fl*pi*((dinifl(icla)*1.d-3)**2))/6.d0
302
 
 
303
 
! --- Calcul de la pression partielle en oxygene (atm)
304
 
!                                                 ---
305
 
!       PO2 = RHO1*RR*T*YO2/MO2
306
 
 
307
 
      pparo2 = propce(iel,ipproc(irom1))*rr*propce(iel,ipcte1)    &
308
 
              *propce(iel,ipcyox)/wmole(io2)
309
 
      pparo2 = pparo2 / prefth
310
 
 
311
 
! --- Calcul de Dcoke en metres
312
 
 
313
 
      dcoke = ( ( rtpa(iel,isca(iyfol(icla)))                     &
314
 
              /(rtpa(iel,isca(ing(icla)))*rho0fl)                 &
315
 
              -pi*(dinikf(icla)**3)*xinkf/6.d0  )                 &
316
 
               *6.d0/(pi*(1.d0-xinkf)) )**(1.d0/3.d0)             &
317
 
            *1.d-3
318
 
     if ( dcoke .lt. 0.d0 ) then
319
 
       WRITE(NFECRA,*) 'erreur Dcoke = ',Dcoke,IEL
320
 
       call csexit(1)
321
 
     endif
322
 
 
323
 
! --- Coefficient de cinetique chimique de formation de CO
324
 
!       en (kg.m-2.s-1.atm(-n))
325
 
 
326
 
      xdffli = ahetfl*exp(-ehetfl*4185.d0                         &
327
 
              /(rr*propce(iel,ipcte1)))
328
 
 
329
 
! --- Coefficient de diffusion en  (Kg/m2/s/atm) : XDFEXT
330
 
!     Coefficient global pour n=0.5 en (kg/m2/s) : XDFTOT0
331
 
!     Coefficient global pour n=1   en (Kg/m2/s) : XDFTOT1
332
 
 
333
 
      diacka = dcoke/(dinikf(icla)*1.d-3)
334
 
      if ( diacka .gt. epsifl ) then
335
 
        xdfext = 2.53d-7*((propce(iel,ipcte1))**0.75d0)           &
336
 
                / dcoke*2.d0
337
 
        xdftot1 = pparo2 / ( 1.d0/xdffli + 1.d0/xdfext )
338
 
        xdftot0 = -(xdffli**2)/(2.d0*xdfext**2)+(pparo2*xdffli**2 &
339
 
                  +(xdffli**4)/(2.d0*xdfext**2))**0.5d0
340
 
      else
341
 
        xdftot1 = xdffli*pparo2
342
 
        xdftot0 = xdffli*pparo2**0.5d0
343
 
      endif
344
 
 
345
 
!     Surface
346
 
 
347
 
      surf = pi*(dcoke**2)
348
 
 
349
 
! --- Calcul de PROPCE(IEL,IPCGHT) = - COXCK*XDFTOT0*PPARO2*XNP < 0
350
 
! --- ou        PROPCE(IEL,IPCGHT) = - COXCK*XDFTOT1*PPARO2*XNP < 0
351
 
 
352
 
      if (iofhet.eq.1) then
353
 
!             PROPCE(IEL,IPCGHT) = - XDFTOT1*COXCK*XNG
354
 
        propce(iel,ipcght) = - xdftot1*surf*xng
355
 
      else
356
 
!             PROPCE(IEL,IPCGHT) = - XDFTOT0*COXCK*XNG
357
 
        propce(iel,ipcght) = - xdftot0*surf*xng
358
 
      endif
359
 
 
360
 
    else
361
 
      propce(iel,ipcght) = 0.d0
362
 
    endif
363
 
 
364
 
  enddo
365
 
 
366
 
enddo
367
 
 
368
 
!===============================================================================
369
 
! FORMATS
370
 
!--------
371
 
!----
372
 
! FIN
373
 
!----
374
 
 
375
 
return
376
 
end subroutine