← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/comb/cs_coal_physprop.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/comb/cs_coal_physprop.f90
 
1
!-------------------------------------------------------------------------------
 
2
 
 
3
! This file is part of Code_Saturne, a general-purpose CFD tool.
 
4
!
 
5
! Copyright (C) 1998-2011 EDF S.A.
 
6
!
 
7
! This program is free software; you can redistribute it and/or modify it under
 
8
! the terms of the GNU General Public License as published by the Free Software
 
9
! Foundation; either version 2 of the License, or (at your option) any later
 
10
! version.
 
11
!
 
12
! This program is distributed in the hope that it will be useful, but WITHOUT
 
13
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 
14
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 
15
! details.
 
16
!
 
17
! You should have received a copy of the GNU General Public License along with
 
18
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
 
19
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
 
20
 
 
21
!-------------------------------------------------------------------------------
 
22
 
 
23
subroutine cs_coal_physprop &
 
24
!==========================
 
25
 
 
26
 ( nvar   , nscal  ,                                              &
 
27
   ibrom  , izfppp ,                                              &
 
28
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
 
29
   coefa  , coefb   )
 
30
 
 
31
!===============================================================================
 
32
! FONCTION :
 
33
! --------
 
34
! ROUTINE PHYSIQUE PARTICULIERE : COMBUSTION CHARBON PULVERISE
 
35
! Calcul de RHO du melange
 
36
!
 
37
! Arguments
 
38
!__________________.____._____.________________________________________________.
 
39
! name             !type!mode ! role                                           !
 
40
!__________________!____!_____!________________________________________________!
 
41
! nvar             ! i  ! <-- ! total number of variables                      !
 
42
! nscal            ! i  ! <-- ! total number of scalars                        !
 
43
! ibrom            ! te ! <-- ! indicateur de remplissage de romb              !
 
44
! izfppp(nfabor)   ! te ! <-- ! numero de zone de la face de bord              !
 
45
!                  !    !     !  pour le module phys. part.                    !
 
46
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
 
47
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
 
48
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
 
49
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
 
50
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
 
51
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
 
52
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
 
53
!  (nfabor, *)     !    !     !                                                !
 
54
!__________________!____!_____!________________________________________________!
 
55
 
 
56
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
 
57
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
 
58
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
 
59
!            --- tableau de travail
 
60
!===============================================================================
 
61
 
 
62
!===============================================================================
 
63
! Module files
 
64
!===============================================================================
 
65
 
 
66
use paramx
 
67
use numvar
 
68
use optcal
 
69
use cstphy
 
70
use cstnum
 
71
use entsor
 
72
use parall
 
73
use ppppar
 
74
use ppthch
 
75
use coincl
 
76
use cpincl
 
77
use ppincl
 
78
use ppcpfu
 
79
use cs_coal_incl
 
80
use mesh
 
81
 
 
82
!===============================================================================
 
83
 
 
84
implicit none
 
85
 
 
86
! Arguments
 
87
 
 
88
integer          nvar   , nscal
 
89
 
 
90
integer          ibrom
 
91
integer          izfppp(nfabor)
 
92
 
 
93
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
 
94
double precision propce(ncelet,*)
 
95
double precision propfa(nfac,*), propfb(nfabor,*)
 
96
double precision coefa(nfabor,*), coefb(nfabor,*)
 
97
 
 
98
! Local variables
 
99
 
 
100
integer          iel, icha, icla, ipcrom, ipcro2
 
101
integer          izone, ifac
 
102
integer          ipbrom, ipcx2c, iromf , ioxy , nbclip1,nbclip2
 
103
 
 
104
double precision x1sro1, x2sro2, srrom1, uns1pw
 
105
double precision x2tot, wmolme, unsro1
 
106
double precision ff3min,ff3max,fvp2max,valmin,valmax
 
107
 
 
108
integer          ipass
 
109
data             ipass /0/
 
110
save             ipass
 
111
 
 
112
integer          iok1,iok2,iok3
 
113
double precision , dimension ( : )     , allocatable :: f1m,f2m,f3m,f4m,f5m
 
114
double precision , dimension ( : )     , allocatable :: f6m,f7m,f8m,f9m
 
115
double precision , dimension ( : )     , allocatable :: enth1 , x2 ,fvp2m
 
116
double precision , dimension ( : )     , allocatable :: xoxyd,enthox
 
117
 
 
118
!===============================================================================
 
119
!
 
120
!===============================================================================
 
121
! 0. ON COMPTE LES PASSAGES
 
122
!===============================================================================
 
123
 
 
124
ipass = ipass + 1
 
125
 
 
126
!===============================================================================
 
127
! 1. INITIALISATIONS A CONSERVER
 
128
!===============================================================================
 
129
 
 
130
!===============================================================================
 
131
! Deallocation dynamic arrays
 
132
!----
 
133
allocate(f1m(1:ncelet), f2m(1:ncelet), f3m(1:ncelet), STAT=iok1)
 
134
allocate(f4m(1:ncelet), f5m(1:ncelet), STAT=iok1)
 
135
allocate(f6m(1:ncelet), f7m(1:ncelet), f8m(1:ncelet), f9m(1:ncelet), STAT=iok2)
 
136
allocate(enth1(1:ncel), x2(1:ncel)   , fvp2m(1:ncel), STAT=iok3)
 
137
!----
 
138
if ( iok1 > 0 .or. iok2 > 0 .or. iok3 > 0) THEN
 
139
  write(nfecra,*) ' Memory allocation error inside: '
 
140
  write(nfecra,*) '     cs_coal_physprop             '
 
141
  call csexit(1)
 
142
ENDIF
 
143
if ( ieqnox .eq. 1 ) then
 
144
  !----
 
145
  allocate(xoxyd(1:ncelet),enthox(1:ncelet),STAT=iok1)
 
146
  !----
 
147
  if ( iok1 > 0 ) then
 
148
    write(nfecra,*) ' Memory allocation error inside:         '
 
149
    write(nfecra,*) '   cs_coal_physprop for xoxyd and enthox '
 
150
  endif
 
151
endif
 
152
!===============================================================================
 
153
!     Pointeur sur masse volumique du gaz aux cellules
 
154
iromf = ipproc(irom1)
 
155
!
 
156
!===============================================================================
 
157
! 2. CALCUL DES PROPRIETES PHYSIQUES DE LA PHASE DISPERSEE
 
158
!                    VALEURS CELLULES
 
159
!                    ----------------
 
160
!    FRACTION MASSIQUE DE SOLIDE
 
161
!    DIAMETRE
 
162
!    MASSE VOLUMIQUE
 
163
!===============================================================================
 
164
 
 
165
call cs_coal_physprop2 ( ncelet , ncel , rtp , propce )
 
166
!=====================
 
167
 
 
168
!===============================================================================
 
169
! 3. CALCUL DES PROPRIETES PHYSIQUES DE LA PHASE GAZEUSE
 
170
!                    VALEURS CELLULES
 
171
!                    ----------------
 
172
!    TEMPERATURE
 
173
!    MASSE VOLUMIQUE
 
174
!    CONCENTRATIONS DES ESPECES GAZEUSES
 
175
!===============================================================================
 
176
 
 
177
! --- Calcul de l'enthalpie du gaz     enth1
 
178
!        de F1M
 
179
!        de F2M
 
180
!        de F3M                    dans W3=1-F1M-F2M-F4M-F5M-F6M-F7M-F8M-F9M
 
181
!        de F4M
 
182
!        de F5M
 
183
!        de F6M
 
184
!        de F7M
 
185
!        de F8M
 
186
!        de F9M
 
187
!        de FVP2M
 
188
!
 
189
! Initialisation des fm et de x2 a 0
 
190
f1m( : ) = 0.d0
 
191
f2m( : ) = 0.d0
 
192
f3m( : ) = 0.d0
 
193
f4m( : ) = 0.d0
 
194
f5m( : ) = 0.d0
 
195
f6m( : ) = 0.d0
 
196
f7m( : ) = 0.d0
 
197
f8m( : ) = 0.d0
 
198
f9m( : ) = 0.d0
 
199
x2 ( : ) = 0.d0
 
200
!
 
201
do icla = 1, nclacp
 
202
  ipcx2c = ipproc(ix2(icla))
 
203
  do iel = 1, ncel
 
204
    x2(iel) =  x2(iel) + propce(iel,ipcx2c)
 
205
  enddo
 
206
enddo
 
207
!
 
208
do icha = 1, ncharb
 
209
  do iel = 1, ncel
 
210
    f1m(iel) =  f1m(iel) + rtp(iel,isca(if1m(icha)))
 
211
    f2m(iel) =  f2m(iel) + rtp(iel,isca(if2m(icha)))
 
212
  enddo
 
213
enddo
 
214
 
 
215
if ( ieqnox .eq. 1 ) then
 
216
  do iel = 1, ncel
 
217
    xoxyd(iel)= (1.d0-x2(iel))-f1m(iel)-f2m(iel)
 
218
  enddo
 
219
endif
 
220
 
 
221
ff3min = 1.d+20
 
222
ff3max =-1.d+20
 
223
nbclip1= 0
 
224
nbclip2= 0
 
225
valmin = 1.d+20
 
226
valmax =-1.d+20
 
227
do iel = 1, ncel
 
228
  uns1pw = 1.d0/(1.d0-x2(iel))
 
229
  if ( noxyd .ge. 2 ) then
 
230
    f4m(iel) = rtp(iel,isca(if4m))
 
231
    if ( noxyd .eq. 3 ) then
 
232
      f5m(iel) = rtp(iel,isca(if5m))
 
233
    endif
 
234
  endif
 
235
 
 
236
  if ( ippmod(iccoal) .ge. 1 ) then
 
237
    f6m(iel) = rtp(iel,isca(if6m))
 
238
  endif
 
239
 
 
240
  f7m(iel) =  rtp(iel,isca(if7m))
 
241
 
 
242
  if ( ihtco2 .eq. 1 ) then
 
243
    f8m(iel) = rtp(iel,isca(if8m))
 
244
  endif
 
245
 
 
246
  if ( ihth2o .eq. 1 ) then
 
247
    f9m(iel) = rtp(iel,isca(if9m))
 
248
  endif
 
249
 
 
250
  fvp2m(iel) = rtp(iel,isca(ifvp2m))
 
251
 
 
252
  f1m(iel)  = f1m(iel)    *uns1pw
 
253
  f2m(iel)  = f2m(iel)    *uns1pw
 
254
  f4m(iel)  = f4m(iel)    *uns1pw
 
255
  f5m(iel)  = f5m(iel)    *uns1pw
 
256
  f6m(iel)  = f6m(iel)    *uns1pw
 
257
  f7m(iel)  = f7m(iel)    *uns1pw
 
258
  f8m(iel)  = f8m(iel)    *uns1pw
 
259
  f9m(iel)  = f9m(iel)    *uns1pw
 
260
!
 
261
!  fvp2m(iel)= fvp2m(iel)  *uns1pw
 
262
!
 
263
! correction variance : uniquement phase 1
 
264
!
 
265
  fvp2m(iel)= fvp2m(iel)*uns1pw - x2(iel)*((f1m(iel)+f2m(iel))**2.d0)
 
266
  if ( fvp2m(iel) .lt. 0.d0 ) then
 
267
    nbclip1 = nbclip1 + 1
 
268
    valmin = min(valmin,fvp2m(iel))
 
269
    fvp2m(iel)= 0.d0
 
270
  endif
 
271
  fvp2max =(f1m(iel)+f2m(iel))*(1.d0-(f1m(iel)+f2m(iel)))
 
272
  if ( fvp2m(iel) .gt. fvp2max  ) then
 
273
    nbclip2 = nbclip2 + 1
 
274
    valmax = max(valmax,fvp2m(iel))
 
275
    fvp2m(iel)= fvp2max
 
276
  endif
 
277
!
 
278
! fin correction
 
279
!
 
280
  f3m(iel) = 1.d0                                        &
 
281
           -( f1m(iel)+f2m(iel)+f4m(iel)+f5m(iel)        &
 
282
             +f6m(iel)+f7m(iel)+f8m(iel)+f9m(iel) )
 
283
 
 
284
  ff3max = max(ff3max,f3m(iel))
 
285
  ff3min = min(ff3min,f3m(iel))
 
286
!
 
287
  if ( ieqnox .eq. 1 ) then
 
288
    enthox(iel) = rtp(iel,isca(ihox))/xoxyd(iel)
 
289
  endif
 
290
!
 
291
enddo
 
292
 
 
293
if ( irangp .ge. 0 ) then
 
294
  !--
 
295
  call parmin(ff3min)
 
296
  call parmax(ff3max)
 
297
  call parcpt(nbclip1)
 
298
  call parcpt(nbclip2)
 
299
  call parmin(valmin)
 
300
  call parmax(valmax)
 
301
  !--
 
302
endif
 
303
WRITE(NFECRA,*) ' Values of F3 min and max: ',FF3MIN,FF3MAX
 
304
if ( nbclip1 .gt. 0 ) then
 
305
  write(nfecra,*) ' Clipping phase gas variance in min :',nbclip1,valmin
 
306
endif
 
307
if ( nbclip2 .gt. 0 ) then
 
308
  write(nfecra,*) ' Clipping phase gas variance in max :',nbclip2,valmax
 
309
endif
 
310
 
 
311
! ---- Enthalpie du gaz H1
 
312
enth1( : ) =0.D0
 
313
do icla = 1, nclacp
 
314
  do iel = 1, ncel
 
315
    enth1(iel) =  enth1(iel) + rtp(iel,isca(ih2(icla)))
 
316
  enddo
 
317
enddo
 
318
do iel = 1, ncel
 
319
  enth1(iel) = (rtp(iel,isca(ihm))-enth1(iel))/ ( 1.d0-x2(iel) )
 
320
enddo
 
321
 
 
322
call cs_coal_physprop1 &
 
323
!=====================
 
324
 ( ncelet , ncel   ,                                      &
 
325
   f1m    , f2m    , f3m    , f4m    , f5m    ,           &
 
326
   f6m    , f7m    , f8m    , f9m    , fvp2m  ,           &
 
327
   enth1  , enthox ,                                      &
 
328
   rtp    , propce , propce(1,iromf)   )
 
329
 
 
330
!===============================================================================
 
331
! 4. CALCUL DES PROPRIETES PHYSIQUES DE LA PHASE DISPERSEE
 
332
!                    VALEURS CELLULES
 
333
!                    ----------------
 
334
!    TEMPERATURE
 
335
!===============================================================================
 
336
 
 
337
! --- Transport d'H2
 
338
 
 
339
call  cs_coal_thfieldconv2  ( ncelet , ncel , rtp , propce )
 
340
!=========================
 
341
 
 
342
!===============================================================================
 
343
! 5. CALCUL DES PROPRIETES PHYSIQUES DU MELANGE
 
344
!                    VALEURS CELLULES
 
345
!                    ----------------
 
346
!    MASSE VOLUMIQUE
 
347
!===============================================================================
 
348
! --- Calcul de Rho du melange : 1/Rho = X1/Rho1 + Somme(X2/Rho2)
 
349
!     On sous relaxe quand on a un rho n a disposition, ie
 
350
!       a partir du deuxieme passage ou
 
351
!       a partir du premier passage si on est en suite de calcul et
 
352
!         qu'on a relu la masse volumique dans le fichier suite.
 
353
 
 
354
ipcrom = ipproc(irom)
 
355
 
 
356
if (ipass.gt.1.or.(isuite.eq.1.and.initro.eq.1)) then
 
357
  srrom1 = srrom
 
358
else
 
359
  srrom1 = 0.d0
 
360
endif
 
361
 
 
362
do iel = 1, ncel
 
363
  x2sro2 = zero
 
364
  do icla = 1, nclacp
 
365
    ipcro2 = ipproc(irom2(icla))
 
366
    ipcx2c = ipproc(ix2(icla))
 
367
    x2sro2 = x2sro2 + propce(iel,ipcx2c) / propce(iel,ipcro2)
 
368
  enddo
 
369
  x1sro1 = (1.d0-x2(iel)) / propce(iel,iromf)
 
370
! ---- Sous relaxation eventuelle a donner dans ppini1.F
 
371
  propce(iel,ipcrom) = srrom1*propce(iel,ipcrom)                  &
 
372
                     + (1.d0-srrom1)/(x1sro1+x2sro2)
 
373
enddo
 
374
 
 
375
 
 
376
!===============================================================================
 
377
! 6. CALCUL DE RHO DU MELANGE
 
378
!                      VALEURS FACETTES
 
379
!                      ----------------
 
380
!===============================================================================
 
381
 
 
382
ibrom = 1
 
383
ipbrom = ipprob(irom)
 
384
ipcrom = ipproc(irom)
 
385
 
 
386
! ---> Masse volumique au bord pour toutes les facettes
 
387
!      Les facettes d'entree seront recalculees.
 
388
 
 
389
do ifac = 1, nfabor
 
390
  iel = ifabor(ifac)
 
391
  propfb(ifac,ipbrom) = propce(iel,ipcrom)
 
392
enddo
 
393
 
 
394
! ---> Masse volumique au bord pour les facettes d'entree UNIQUEMENT
 
395
!     Le test sur IZONE sert pour les reprises de calcul
 
396
 
 
397
if ( ipass.gt.1 .or. isuite.eq.1 ) then
 
398
  do ifac = 1, nfabor
 
399
 
 
400
    izone = izfppp(ifac)
 
401
    if(izone.gt.0) then
 
402
      if ( ientat(izone).eq.1 .or. ientcp(izone).eq.1 ) then
 
403
        x2sro2 = zero
 
404
        x2tot  = zero
 
405
        do icla = 1, nclacp
 
406
          x2sro2 = x2sro2 + x20(izone,icla)/rho20(icla)
 
407
          x2tot  = x2tot  + x20(izone,icla)
 
408
        enddo
 
409
 
 
410
        ioxy = inmoxy(izone)
 
411
        wmolme =( oxyo2(ioxy)+oxyn2(ioxy)                         &
 
412
                 +oxyh2o(ioxy)+oxyco2(ioxy))                      &
 
413
               /( wmole(io2) *oxyo2(ioxy)                         &
 
414
                 +wmole(in2) *oxyn2(ioxy)                         &
 
415
                 +wmole(ih2o)*oxyh2o(ioxy)                        &
 
416
                 +wmole(ico2)*oxyco2(ioxy) )
 
417
 
 
418
        unsro1 = (wmolme*rr*timpat(izone)) / p0
 
419
        x1sro1 = (1.d0-x2tot) * unsro1
 
420
        propfb(ifac,ipbrom) = 1.d0 / (x1sro1+x2sro2)
 
421
      endif
 
422
    endif
 
423
 
 
424
  enddo
 
425
endif
 
426
 
 
427
!----
 
428
! Formats
 
429
!----
 
430
 
 
431
!===============================================================================
 
432
! Deallocation dynamic arrays
 
433
!----
 
434
deallocate(f1m,f2m,f3m,f4m,f5m,STAT=iok1)
 
435
deallocate(f6m,f7m,f8m,f9m,    STAT=iok2)
 
436
deallocate(enth1,x2,fvp2m,     STAT=iok3)
 
437
!----
 
438
IF ( iok1 > 0 .or. iok2 > 0 .or. iok3 > 0 ) THEN
 
439
  write(nfecra,*) ' Memory deallocation error inside: '
 
440
  write(nfecra,*) '     cs_coal_physprop              '
 
441
  call csexit(1)
 
442
ENDIF
 
443
if ( ieqnox .eq. 1 ) then
 
444
  !--
 
445
  deallocate(xoxyd,enthox)
 
446
  !--
 
447
  IF ( iok1 > 0 ) THEN
 
448
    write(nfecra,*) ' Memory deallocation error inside:       '
 
449
    write(nfecra,*) '   cs_coal_physprop for xoxyd and enthox '
 
450
    call csexit(1)
 
451
  endif
 
452
endif
 
453
!===============================================================================
 
454
!
 
455
!----
 
456
! End
 
457
!----
 
458
return
 
459
end subroutine