← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/fuel/fucyno.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/fucyno.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 fucyno &
29
 
!================
30
 
 
31
 
 ( ncelet , ncel   ,                                              &
32
 
   indpdf ,                                                       &
33
 
   rtp    , propce ,                                              &
34
 
   f1m    , f3m    , f4m   , f1cl  , f3cl , f4cl ,                &
35
 
   f4m1   , f4m2   , d4cl  , d4f4  , hrec , f4s3 ,                &
36
 
   fuel1  , fuel3  , oxyd  , prod1 , prod2  ,                     &
37
 
   xiner  , xh2s   , xso2   )
38
 
 
39
 
!===============================================================================
40
 
! FONCTION :
41
 
! --------
42
 
 
43
 
! CALCUL DE :
44
 
 
45
 
!          K1 exp(-E1/RT)   (conversion HCN en NO)
46
 
!          K2 exp(-E2/RT)   (conversion HCN en N2)
47
 
!          K3 exp(-E3/RT)   (NO thermique)
48
 
 
49
 
!  VALEURS CELLULES
50
 
!  ----------------
51
 
 
52
 
! Arguments
53
 
!__________________.____._____.________________________________________________.
54
 
! name             !type!mode ! role                                           !
55
 
!__________________!____!_____!________________________________________________!
56
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
57
 
! ncel             ! i  ! <-- ! number of cells                                !
58
 
! rtp              ! tr ! <-- ! variables de calcul au centre des              !
59
 
! (ncelet,*)       !    !     !    cellules (instant courant)                  !
60
 
! fvap             ! tr ! <-- ! moyenne du traceur 1 mvl [fovm+co]             !
61
 
! fhet             ! tr !  ->           ! moyenne du traceur 3 c h�terog�ne    !
62
 
! f4p2m            ! tr ! <-- ! variance du traceur 4 (air)                    !
63
 
! indpdf           ! te ! <-- ! passage par les pdf                            !
64
 
! fuel1            ! tr !  <- ! fraction massique fovm                         !
65
 
! fuel3            ! tr !  <- ! fraction massique co                           !
66
 
! oxyd             ! tr !  <- ! fraction massique o2                           !
67
 
! prod1            ! tr !  <- ! fraction massique co2                          !
68
 
! prod2            ! tr !  <- ! fraction massique h2o                          !
69
 
! xiner            ! tr !  <- ! fraction massique n2                           !
70
 
! xh2s             ! tr !  <- ! fraction massique h2s                          !
71
 
! xso2             ! tr !  <- ! fraction massique so2                          !
72
 
! f4m1             ! tr ! <-- ! borne minimum                                  !
73
 
! f4m2             ! tr ! <-- ! borne max                                      !
74
 
! d4cl             ! tr ! <-- ! amplitude du pic de dirac en f4cl              !
75
 
! d4f4             ! tr ! <-- ! amplitude du pic de dirac en 1                 !
76
 
!                  !    !     ! (air pur)                                      !
77
 
!__________________!____!_____!________________________________________________!
78
 
 
79
 
!     TYPE : E (ENTIER), R (REEL), A (ALPHAMNUMERIQUE), T (TABLEAU)
80
 
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
81
 
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
82
 
!            --- tableau de travail
83
 
!===============================================================================
84
 
 
85
 
implicit none
86
 
 
87
 
!==============================================================================
88
 
! Common blocks
89
 
!==============================================================================
90
 
 
91
 
include "paramx.h"
92
 
include "numvar.h"
93
 
include "optcal.h"
94
 
include "cstphy.h"
95
 
include "cstnum.h"
96
 
include "entsor.h"
97
 
include "parall.h"
98
 
include "pointe.h"
99
 
include "ppppar.h"
100
 
include "ppthch.h"
101
 
include "coincl.h"
102
 
include "cpincl.h"
103
 
include "fuincl.h"
104
 
include "ppincl.h"
105
 
include "ppcpfu.h"
106
 
 
107
 
!===============================================================================
108
 
 
109
 
! Arguments
110
 
 
111
 
integer          ncelet , ncel
112
 
integer          indpdf(ncelet)
113
 
 
114
 
double precision rtp(ncelet,*), propce(ncelet,*)
115
 
double precision f1m(ncelet) , f3m(ncelet)  , f4m(ncelet)
116
 
double precision f1cl(ncelet), f3cl(ncelet) , f4cl(ncelet)
117
 
 
118
 
double precision f4m1(ncelet) , f4m2(ncelet) , d4cl(ncelet)
119
 
double precision d4f4(ncelet) , hrec(ncelet) , f4s3(ncel)
120
 
 
121
 
double precision fuel1(ncelet), fuel3(ncelet)
122
 
double precision oxyd(ncelet) , prod1(ncelet), prod2(ncelet)
123
 
double precision xiner(ncelet)
124
 
double precision xh2s(ncelet) , xso2(ncelet)
125
 
 
126
 
! Local variables
127
 
 
128
 
integer          iel , icla , i
129
 
 
130
 
integer          n1 , n2 , n3 , n4 , n5 , n6 , n7 , n8 , n9
131
 
integer          n10, n11
132
 
integer          nbrint
133
 
integer          ipdf1 , ipdf2 , ipdf3
134
 
integer          iexp1 , iexp2 , iexp3
135
 
parameter        (nbrint = 11)
136
 
integer          inttmp(nbrint)
137
 
integer          npart
138
 
parameter        (npart = 200 )
139
 
 
140
 
double precision wmo2
141
 
 
142
 
double precision ee1,ee2,ee3,kk1,kk2,kk3
143
 
double precision ts3 , ts3den , ts3num , u , t , ts3s , dirac
144
 
double precision q , r , s , lpa , lri , tmpgaz , tfuel , xmx2
145
 
double precision bb1 , bb2 , bb3 , bb4 , wmco , xo2 , bb
146
 
double precision gh1, gh2 ,gh3 ,gh4 ,gh5,gh6, gh7 , tg
147
 
double precision gh8, gh9 ,gh10,gh11
148
 
double precision dgs , rsf4 , rscl , bmoy , p , spdf
149
 
double precision val(npart+1),tt(npart+1) , gs(npart+1)
150
 
 
151
 
!===============================================================================
152
 
 
153
 
 
154
 
!===============================================================================
155
 
! 0. INITIALISATIONS
156
 
!===============================================================================
157
 
 
158
 
 
159
 
!===============================================================================
160
 
! 1. CALCULS PRELIMINAIRES
161
 
!===============================================================================
162
 
 
163
 
! --- Masses molaires
164
 
 
165
 
wmo2   = wmole(io2)
166
 
 
167
 
! --- pointeurs
168
 
 
169
 
iexp1 = ipproc(ighcn1)
170
 
iexp2 = ipproc(ighcn2)
171
 
iexp3 = ipproc(ignoth)
172
 
 
173
 
! Parametres des lois d'Arrhenius
174
 
 
175
 
kk1 = 3.0e12
176
 
ee1 = 3.0e4
177
 
kk2 = 1.2e10
178
 
ee2 = 3.35e4
179
 
kk3 = 3.4e12
180
 
ee3 = 6.69e4
181
 
 
182
 
! Pour les termes, indicateur de calcul par PDF ou non
183
 
!       = 1 --> passage par pdf
184
 
!       = 0 --> on ne passe pas par les pdf
185
 
 
186
 
ipdf1 = 0
187
 
ipdf2 = 0
188
 
ipdf3 = 0
189
 
 
190
 
! ---- Initialisation
191
 
 
192
 
do iel = 1, ncel
193
 
!        PROPCE(IEL,IEXP1)  = ZERO
194
 
!        PROPCE(IEL,IEXP2)  = ZERO
195
 
!        PROPCE(IEL,IEXP3)  = ZERO
196
 
enddo
197
 
 
198
 
!===============================================================================
199
 
! 3. CALCUL SANS LES PDF
200
 
!===============================================================================
201
 
 
202
 
do iel = 1, ncel
203
 
 
204
 
  tg  = propce(iel,ipproc(itemp1))
205
 
  propce(iel,iexp1)  = kk1*exp(-ee1/tg)
206
 
 
207
 
enddo
208
 
 
209
 
do iel = 1, ncel
210
 
 
211
 
  tg  = propce(iel,ipproc(itemp1))
212
 
  xo2 = oxyd(iel)*propce(iel,ipproc(immel))                       &
213
 
       /wmo2
214
 
  if ( xo2 .gt. 0.018d0 ) then
215
 
    bb = 0.d0
216
 
  else if ( xo2 .lt. 0.0025d0 ) then
217
 
    bb= 1.d0
218
 
  else
219
 
    bb=(0.018d0-xo2)/(0.018d0-0.0025d0)
220
 
  endif
221
 
 
222
 
  propce(iel,iexp2)  = kk2*exp(-ee2/tg)*(xo2**bb)
223
 
 
224
 
enddo
225
 
 
226
 
do iel = 1, ncel
227
 
 
228
 
  tg  = propce(iel,ipproc(itemp1))
229
 
  xo2 = oxyd(iel)*propce(iel,ipproc(immel))                       &
230
 
       /wmo2
231
 
 
232
 
  propce(iel,iexp3)  = kk3*exp(-ee3/tg)*(xo2**0.5d0)
233
 
 
234
 
enddo
235
 
 
236
 
!===============================================================================
237
 
! 3. CALCUL AVEC LES PDF
238
 
!===============================================================================
239
 
 
240
 
if ( ipdf1 .eq. 1 .or. ipdf2 .eq. 1 .or. ipdf3 .eq. 1 ) then
241
 
 
242
 
  do iel=1,ncel
243
 
 
244
 
    if ( indpdf(iel).eq.1 ) then
245
 
 
246
 
!         Calcul de Tfioul moyen
247
 
 
248
 
      xmx2  = 0.d0
249
 
      tfuel = 0.d0
250
 
 
251
 
      do icla=1,nclafu
252
 
        xmx2 = xmx2 + rtp(iel,isca(iyfol(icla)))
253
 
      enddo
254
 
 
255
 
      if ( xmx2 .gt. 0.d0 ) then
256
 
        do icla=1,nclafu
257
 
          tfuel = tfuel + rtp(iel,isca(iyfol(icla)))              &
258
 
                         *propce(iel,ipproc(itemp3(icla)))
259
 
        enddo
260
 
        tfuel = tfuel/xmx2
261
 
      else
262
 
        tfuel = propce(iel,ipproc(itemp1))
263
 
      endif
264
 
 
265
 
!  On recupere la valeur de TAIR
266
 
 
267
 
      taire = rtp(iel,isca(itaire))
268
 
 
269
 
! On initialise par les temperatures Tair (en f4) et Tcl (TFUEL) aux extr��mit��s
270
 
 
271
 
      dirac  = d4cl(iel)*tfuel + d4f4(iel)*taire
272
 
 
273
 
!1 On recupere la valeur de la temperature moyenne
274
 
 
275
 
      tmpgaz = propce(iel,ipproc(itemp1))
276
 
 
277
 
!1 Integration premier intervalle, entre CL et S3 (stoechio derni��re r��action)
278
 
 
279
 
!1 Parametres de la droite entre CL et S3
280
 
 
281
 
!1 Detection des bornes d'integration (celles de la pdf) du domaine pauvre
282
 
 
283
 
      if(tfuel.gt.epsicp .and. tmpgaz.gt.epsicp) then
284
 
 
285
 
        if(f4s3(iel).le.f4m1(iel)) then
286
 
          p=((f4m1(iel)+f4m2(iel))/2.-f4s3(iel))/(1.-f4s3(iel))
287
 
          ts3s=((tmpgaz-dirac)/hrec(iel)                          &
288
 
                              /(f4m2(iel)-f4m1(iel))-taire*p)     &
289
 
                              /(1.d0-p)
290
 
        else if(f4s3(iel).ge.f4m2(iel)) then
291
 
          ts3s=((tmpgaz-dirac)/hrec(iel)                          &
292
 
               /(f4m2(iel)-f4m1(iel))-tfuel)                      &
293
 
               *(f4s3(iel)-f4cl(iel))                             &
294
 
               /((f4m1(iel)+f4m2(iel))/2.d0-f4cl(iel))+tfuel
295
 
        else
296
 
          ts3s=(2.*(tmpgaz-dirac)/hrec(iel)                       &
297
 
               - tfuel*(f4s3(iel)-f4m1(iel))**2                   &
298
 
                /(f4s3(iel)-f4cl(iel))                            &
299
 
               - taire*(f4m2(iel)-f4s3(iel))**2                   &
300
 
                /(1.d0-f4s3(iel)))                                &
301
 
      /((2.*f4m2(iel)-f4m1(iel)-f4s3(iel))                        &
302
 
               +(f4m1(iel)-f4cl(iel))*(f4s3(iel)-f4m1(iel))       &
303
 
                                     /(f4s3(iel)-f4cl(iel))       &
304
 
               -(f4m2(iel)-f4s3(iel))**2/(1.d0-f4s3(iel)))
305
 
        endif
306
 
 
307
 
        bb1 = max( f4m1(iel) , f4cl(iel) )
308
 
        bb2 = min( f4m2(iel) , f4s3(iel) )
309
 
        bb3 = max( f4m1(iel) , f4s3(iel) )
310
 
        bb4 = min( f4m2(iel) , 1.d0 )
311
 
 
312
 
! On definit quelques parametres intermediaires pour simplifier le code
313
 
 
314
 
        if (bb2 .gt. bb1 ) then
315
 
          lri = hrec(iel)
316
 
        else
317
 
          lri = 0.d0
318
 
        endif
319
 
        if (bb4 .gt. bb3 ) then
320
 
          lpa = hrec(iel)
321
 
        else
322
 
          lpa = 0.d0
323
 
        endif
324
 
 
325
 
        p = f4s3(iel) - f4cl(iel)
326
 
        q = bb2**2 - bb1**2
327
 
        r = bb2 - bb1
328
 
        s = 1.d0 - f4s3(iel)
329
 
        t = bb4**2 - bb3**2
330
 
        u = bb4 - bb3
331
 
 
332
 
! On calcul de TS3
333
 
 
334
 
!           TS3=(TMPGAZ-DIRAC+TFUEL*(((LRI*Q)/(2.D0*P))-((LRI*F4S3*R)/P))
335
 
!     &         - TAIRE*( ((LPA*T)/(2.D0*S)) +(LPA*U) -((LPA*U)/S)))
336
 
!     &         /
337
 
!     &         (((LRI*Q)/(2.D0*P))+(LRI*R)-((LRI*F4S3*R)/P)-((LPA*T)
338
 
!     &         /(2.D0*S)) +((LPA*U)/S) )
339
 
 
340
 
        gh1 = -lri*q/(p*2.d0)
341
 
        gh2 =  lri*r*f4s3(iel)/p
342
 
        gh3 =  lpa*t/(2.d0*s)
343
 
        gh4 =  lpa*u
344
 
        gh5 =  lpa*u/s
345
 
        gh6 =  lri*q/(p*2.d0)
346
 
        gh7 =  lri*r
347
 
        gh8 =  f4s3(iel)*lri*r/p
348
 
        gh9 = lpa*t/(s*2.d0)
349
 
        gh10= lpa*u/s
350
 
        gh11 = tmpgaz - dirac
351
 
 
352
 
        ts3num = gh11 - tfuel*(gh1+gh2) - taire*(gh3+gh4-gh5)
353
 
        ts3den = gh6 + gh7 -gh8 - gh9 + gh10
354
 
 
355
 
 
356
 
        ts3 = ts3num / ts3den
357
 
        if(abs(ts3-ts3s).gt.1.d0) then
358
 
          WRITE(NFECRA,*) 'TS3 paul-TS3 sandro ', IEL,TS3,TS3S
359
 
        endif
360
 
      endif
361
 
 
362
 
      spdf = hrec(iel) * (f4m2(iel) - f4m1(iel))
363
 
 
364
 
! Integration
365
 
 
366
 
      xo2 = oxyd(iel)*propce(iel,ipproc(immel))                   &
367
 
           /wmo2
368
 
 
369
 
      do i = 1, npart+1
370
 
        gs(i) = f4m1(iel) + dble(i-1)/dble(npart)                 &
371
 
                           *(f4m2(iel)-f4m1(iel))
372
 
 
373
 
        if( gs(i) .le. f4s3(iel) ) then
374
 
          tt(i) = (ts3-tfuel)/(f4s3(iel)-f4cl(iel))* gs(i)        &
375
 
               + ts3 - f4s3(iel)*(ts3 - tfuel)                    &
376
 
                                / ( f4s3(iel) - f4cl(iel) )
377
 
 
378
 
        else
379
 
          tt(i) = (taire -ts3)/(1.d0-f4s3(iel))*gs(i)             &
380
 
               + taire - (taire - ts3)/(1.d0-f4s3(iel))
381
 
 
382
 
        endif
383
 
      enddo
384
 
 
385
 
      if(xo2.gt.0.018d0) then
386
 
        bmoy=0.d0
387
 
      else if(xo2 .lt. 0.0025d0) then
388
 
        bmoy=1.d0
389
 
      else
390
 
        bmoy=(0.018d0-xo2)/(0.018d0-0.0025d0)
391
 
      endif
392
 
 
393
 
!      DGS est le pas d'integration
394
 
 
395
 
      dgs = ( f4m2(iel) - f4m1(iel) ) / dble(npart)
396
 
 
397
 
! Calcul de K1*EXP(-E1/T)
398
 
 
399
 
      if ( ipdf1 .eq. 1 ) then
400
 
 
401
 
        rsf4 = kk1*exp(-ee1/taire)*d4f4(iel)
402
 
        rscl = kk1*exp(-ee1/tfuel)*d4cl(iel)
403
 
 
404
 
        do i = 1, npart+1
405
 
          val(i) = kk1 * exp(-ee1/tt(i)) * hrec(iel)
406
 
        enddo
407
 
 
408
 
        propce(iel,iexp1)= rsf4 + rscl
409
 
 
410
 
        do i = 1, npart
411
 
          propce(iel,iexp1) = propce(iel,iexp1)                   &
412
 
                           +0.5d0*dgs*(val(i)+val(i+1))
413
 
        enddo
414
 
 
415
 
      endif
416
 
 
417
 
!  Calcul de K2*EXP(-E2/T)
418
 
 
419
 
      if ( ipdf2 .eq. 1 ) then
420
 
 
421
 
        if ( xo2 .gt. 0.d0 ) then
422
 
          propce(iel,iexp2) = kk2*exp(-ee2/taire)                 &
423
 
                                 *d4f4(iel)*xo2**bmoy             &
424
 
                             +kk2*exp(-ee2/tfuel)                 &
425
 
                                 *d4cl(iel)*xo2**bmoy
426
 
 
427
 
          do i = 1, npart+1
428
 
            val(i) = kk2*exp(-ee2/tt(i))*hrec(iel)
429
 
          enddo
430
 
 
431
 
          do i = 1, npart
432
 
            propce(iel,iexp2) = propce(iel,iexp2)                 &
433
 
                               +0.5d0*dgs*(val(i)+val(i+1))       &
434
 
                                     *xo2**bmoy
435
 
          enddo
436
 
        else
437
 
          propce(iel,iexp2) = 0.d0
438
 
        endif
439
 
 
440
 
      endif
441
 
 
442
 
!  Calcul de K3*EXP(-E3/T)
443
 
 
444
 
      if ( ipdf3 .eq. 1 ) then
445
 
 
446
 
        if ( xo2 .gt. 0.d0 ) then
447
 
          propce(iel,iexp3) = kk3*exp(-ee3/taire)                 &
448
 
                                 *d4f4(iel)*xo2**(0.5d0)          &
449
 
                             +kk3*exp(-ee3/tfuel)                 &
450
 
                                 *d4cl(iel)*xo2**(0.5d0)
451
 
 
452
 
          do i = 1, npart+1
453
 
            val(i) = kk3*exp(-ee3/tt(i))*hrec(iel)
454
 
          enddo
455
 
 
456
 
          do i = 1, npart
457
 
            propce(iel,iexp3)= propce(iel,iexp3)                  &
458
 
                     + 0.5d0*dgs*(val(i)+val(i+1))*xo2**(0.5d0)
459
 
          enddo
460
 
        else
461
 
          propce(iel,iexp3)= 0.d0
462
 
        endif
463
 
      endif
464
 
 
465
 
    endif
466
 
 
467
 
  enddo
468
 
 
469
 
endif
470
 
 
471
 
return
472
 
end subroutine