← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/mati/mttssc.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/mati/mttssc.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 mttssc &
29
 
!================
30
 
 
31
 
 ( idbia0 , idbra0 ,                                              &
32
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
33
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
34
 
   nvar   , nscal  , nphas  , ncepdp , ncesmp ,                   &
35
 
   nideve , nrdeve , nituse , nrtuse , iscal  ,                   &
36
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
37
 
   ipnfac , nodfac , ipnfbr , nodfbr , icepdc , icetsm , itypsm , &
38
 
   iconra ,                                                       &
39
 
   idevel , ituser , ia     ,                                     &
40
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
41
 
   dt     , rtpa   , rtp    , propce , propfa , propfb ,          &
42
 
   coefa  , coefb  , ckupdc , smacel ,                            &
43
 
   crvexp , crvimp ,                                              &
44
 
   viscf  , viscb  , xam    ,                                     &
45
 
   w1     , w2     , w3     , w4     , w5     ,                   &
46
 
   w6     , w7     , w8     , w9     , w10    , w11    ,          &
47
 
   rdevel , rtuser , ra     )
48
 
 
49
 
!===============================================================================
50
 
! FONCTION :
51
 
! ----------
52
 
 
53
 
 
54
 
! CALCUL DES TERMES SOURCES POUR LE SCALAIRE ISCAL
55
 
 
56
 
!    SOUS-PROGRAMME SPECIFIQUE A MATISSE (COPIE DE USTSSC)
57
 
 
58
 
!    LES SEULS SCALAIRES TRAITES SONT
59
 
!      ISCA(ITAAMT) TEMPERATURE DE L'AIR
60
 
!      ISCA(ITPCMT) TEMPERATURE DE PEAU DES COLIS
61
 
!      ISCA(ITPPMT) TEMPERATURE DE PEAU DES MURS
62
 
 
63
 
 
64
 
! ON RESOUT RHO*VOLUME*D(VAR)/DT = CRVIMP*VAR + CRVEXP
65
 
 
66
 
! ON FOURNIT ICI CRVIMP ET CRVEXP (ILS CONTIENNENT RHO*VOLUME)
67
 
!    CRVEXP en kg variable/s :
68
 
!     ex : pour les temperatures      kg degres/s
69
 
!          pour les enthalpies        Joules/s
70
 
!    CRVIMP en kg /s :
71
 
 
72
 
! VEILLER A UTILISER UN CRVIMP NEGATIF
73
 
! (ON IMPLICITERA CRVIMP
74
 
!  IE SUR LA DIAGONALE DE LA MATRICE, LE CODE AJOUTERA :
75
 
!   MAX(-CRVIMP,0) EN SCHEMA STANDARD EN TEMPS
76
 
!       -CRVIMP    SI LES TERMES SOURCES SONT A L'ORDRE 2
77
 
 
78
 
! CES TABLEAUX SONT INITIALISES A ZERO AVANT APPEL A CE SOUS
79
 
!   PROGRAMME ET AJOUTES ENSUITE AUX TABLEAUX PRIS EN COMPTE
80
 
!   POUR LA RESOLUTION
81
 
 
82
 
! EN CAS D'ORDRE 2 DEMANDE SUR LES TERMES SOURCES, ON DOIT
83
 
!   FOURNIR CRVEXP A L'INSTANT N     (IL SERA EXTRAPOLE) ET
84
 
!           CRVIMP A L'INSTANT N+1/2 (IL EST  DANS LA MATRICE,
85
 
!                                     ON LE SUPPOSE NEGATIF)
86
 
 
87
 
!-------------------------------------------------------------------------------
88
 
!ARGU                             ARGUMENTS
89
 
!__________________.____._____.________________________________________________.
90
 
! name             !type!mode ! role                                           !
91
 
!__________________!____!_____!________________________________________________!
92
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
93
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
94
 
! ndim             ! i  ! <-- ! spatial dimension                              !
95
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
96
 
! ncel             ! i  ! <-- ! number of cells                                !
97
 
! nfac             ! i  ! <-- ! number of interior faces                       !
98
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
99
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
100
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
101
 
! nnod             ! i  ! <-- ! number of vertices                             !
102
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
103
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
104
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
105
 
! nvar             ! i  ! <-- ! total number of variables                      !
106
 
! nscal            ! i  ! <-- ! total number of scalars                        !
107
 
! nphas            ! i  ! <-- ! number of phases                               !
108
 
! ncepdp           ! i  ! <-- ! number of cells with head loss                 !
109
 
! ncesmp           ! i  ! <-- ! number of cells with mass source term          !
110
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
111
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
112
 
! iscal            ! i  ! <-- ! scalar number                                  !
113
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
114
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
115
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
116
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
117
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
118
 
!  (nfml, nprfml)  !    !     !                                                !
119
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
120
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
121
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
122
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
123
 
! icepdc(ncelet    ! te ! <-- ! numero des ncepdp cellules avec pdc            !
124
 
! icetsm(ncesmp    ! te ! <-- ! numero des cellules a source de masse          !
125
 
! itypsm           ! te ! <-- ! type de source de masse pour les               !
126
 
! (ncesmp,nvar)    !    !     !  variables (cf. ustsma)                        !
127
 
! iconra           ! te ! <-- ! tab de connectivite pour                       !
128
 
! (ncelet+1)       !    !     !  le rayonnement et les panaches                !
129
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
130
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
131
 
! ia(*)            ! ia ! --- ! main integer work array                        !
132
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
133
 
!  (ndim, ncelet)  !    !     !                                                !
134
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
135
 
!  (ndim, nfac)    !    !     !                                                !
136
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
137
 
!  (ndim, nfabor)  !    !     !                                                !
138
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
139
 
!  (ndim, nfac)    !    !     !                                                !
140
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
141
 
!  (ndim, nfabor)  !    !     !                                                !
142
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
143
 
!  (ndim, nnod)    !    !     !                                                !
144
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
145
 
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
146
 
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
147
 
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
148
 
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
149
 
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
150
 
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
151
 
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
152
 
!  (nfabor, *)     !    !     !                                                !
153
 
! ckupdc           ! tr ! <-- ! tableau de travail pour pdc                    !
154
 
!  (ncepdp,6)      !    !     !                                                !
155
 
! smacel           ! tr ! <-- ! valeur des variables associee a la             !
156
 
! (ncesmp,*   )    !    !     !  source de masse                               !
157
 
!                  !    !     !  pour ivar=ipr, smacel=flux de masse           !
158
 
! crvexp(ncelet    ! tr ! --> ! tableau de travail pour part explicit          !
159
 
! crvimp(ncelet    ! tr ! --> ! tableau de travail pour terme instat           !
160
 
! viscf(nfac)      ! tr ! --- ! tableau de travail    faces internes           !
161
 
! viscb(nfabor     ! tr ! --- ! tableau de travail    faces de bord            !
162
 
! xam(nfac,2)      ! tr ! --- ! tableau de travail    faces de bord            !
163
 
! w1..11(ncelet    ! tr ! --- ! tableau de travail    cellules                 !
164
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
165
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
166
 
! ra(*)            ! ra ! --- ! main real work array                           !
167
 
!__________________!____!_____!________________________________________________!
168
 
 
169
 
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
170
 
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
171
 
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
172
 
!            --- tableau de travail
173
 
!-------------------------------------------------------------------------------
174
 
!===============================================================================
175
 
 
176
 
implicit none
177
 
 
178
 
!===============================================================================
179
 
! Common blocks
180
 
!===============================================================================
181
 
 
182
 
include "cstnum.h"
183
 
include "paramx.h"
184
 
include "pointe.h"
185
 
include "numvar.h"
186
 
include "entsor.h"
187
 
include "optcal.h"
188
 
include "cstphy.h"
189
 
include "parall.h"
190
 
include "period.h"
191
 
include "matiss.h"
192
 
 
193
 
 
194
 
!===============================================================================
195
 
 
196
 
! Arguments
197
 
 
198
 
integer          idbia0 , idbra0
199
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
200
 
integer          nfml   , nprfml
201
 
integer          nnod   , lndfac , lndfbr , ncelbr
202
 
integer          nvar   , nscal  , nphas
203
 
integer          ncepdp , ncesmp
204
 
integer          nideve , nrdeve , nituse , nrtuse
205
 
integer          iscal
206
 
 
207
 
integer          ifacel(2,nfac) , ifabor(nfabor)
208
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
209
 
integer          iprfml(nfml,nprfml)
210
 
integer          ipnfac(nfac+1), nodfac(lndfac)
211
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
212
 
integer          icepdc(ncepdp)
213
 
integer          icetsm(ncesmp), itypsm(ncesmp,nvar)
214
 
integer          iconra(ncelet)
215
 
integer          idevel(nideve), ituser(nituse), ia(*)
216
 
 
217
 
double precision xyzcen(ndim,ncelet)
218
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
219
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
220
 
double precision xyznod(ndim,nnod), volume(ncelet)
221
 
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
222
 
double precision propce(ncelet,*)
223
 
double precision propfa(nfac,*), propfb(nfabor,*)
224
 
double precision coefa(nfabor,*), coefb(nfabor,*)
225
 
double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar)
226
 
double precision crvexp(ncelet), crvimp(ncelet)
227
 
double precision viscf(nfac), viscb(nfabor)
228
 
double precision xam(nfac,2)
229
 
double precision w1(ncelet), w2(ncelet), w3(ncelet)
230
 
double precision w4(ncelet), w5(ncelet), w6(ncelet)
231
 
double precision w7(ncelet), w8(ncelet), w9(ncelet)
232
 
double precision w10(ncelet), w11(ncelet)
233
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
234
 
 
235
 
! Local variables
236
 
 
237
 
character*80     chaine
238
 
integer          idebia, idebra
239
 
integer          ivar  , iiscvr, ipcrom, iphas
240
 
integer          icoul , ifml  , iel   , jel
241
 
integer          modntl
242
 
integer          nzones, izone
243
 
double precision xnu0  , xlamb0, prand0, phi0
244
 
double precision gmat  , untier
245
 
double precision fform , emis
246
 
double precision hplus , hmoin , hmax  , difx  , dify  , dhmur
247
 
double precision fragh , xxav  , tbulk , un0   , ptot
248
 
double precision raleiz, reynol
249
 
double precision xnsraz, xnsrey, xnsalv
250
 
double precision hraz  , hrey  , h0    , hray
251
 
double precision hmrey , hmraz , hm0   , hmray , hnb
252
 
double precision cpcont, cpparo
253
 
double precision xlimin, xlimax, yrgmin, yrgmax, zalmin, zalmax
254
 
double precision sqz
255
 
double precision dhe   , dhs   , ptcn
256
 
double precision fptsto, pvccsc, pccsc
257
 
 
258
 
!===============================================================================
259
 
!===============================================================================
260
 
! 1. INITIALISATION
261
 
!===============================================================================
262
 
 
263
 
! --- Gestion memoire
264
 
idebia = idbia0
265
 
idebra = idbra0
266
 
 
267
 
 
268
 
! --- Numero du scalaire a traiter : ISCAL (argument)
269
 
 
270
 
! --- Numero de la variable associee au scalaire a traiter ISCAL
271
 
ivar = isca(iscal)
272
 
 
273
 
! --- Nom de la variable associee au scalaire a traiter ISCAL
274
 
chaine = nomvar(ipprtp(ivar))
275
 
 
276
 
! --- Indicateur de variance
277
 
!         Si ISCAVR = 0 :
278
 
!           le scalaire ISCAL n'est pas une variance
279
 
!         Si ISCAVR > 0 et ISCAVR < NSCAL + 1 :
280
 
!           le scalaire ISCAL est une variance associee
281
 
!           au scalaire ISCAVR
282
 
iiscvr = iscavr(iscal)
283
 
 
284
 
! --- Numero de phase associee au scalaire ISCAL
285
 
iphas = iphsca(iscal)
286
 
 
287
 
! --- Masse volumique
288
 
ipcrom = ipproc(irom(iphas))
289
 
 
290
 
 
291
 
! --- Reperage du pas de temps pour les impressions listing
292
 
if(ntlist.gt.0) then
293
 
  modntl = mod(ntcabs,ntlist)
294
 
elseif(ntlist.eq.-1.and.ntcabs.eq.ntmabs) then
295
 
  modntl = 0
296
 
else
297
 
  modntl = 1
298
 
endif
299
 
 
300
 
! --- Impression
301
 
if(iwarni(ivar).ge.1) then
302
 
  write(nfecra,1000) chaine(1:8)
303
 
endif
304
 
if((irangp.le.0).and.(modntl.eq.0)) then
305
 
  write(nfecra,2001) chaine(1:8)
306
 
endif
307
 
 
308
 
! --- Numerique
309
 
untier = 1.d0/3.d0
310
 
 
311
 
!===============================================================================
312
 
! 1. DIMENSIONS GEOMETRIQUES
313
 
!===============================================================================
314
 
 
315
 
 
316
 
!===============================================================================
317
 
! 2. PROPRIETES PHYSIQUES
318
 
!===============================================================================
319
 
 
320
 
! --- Proprietes de l'air
321
 
 
322
 
!     Nombre de Prandtl moleculaire
323
 
prand0 = 0.7d0
324
 
!     Viscosite cinematique moleculaire (nu en m2.s-1)
325
 
xnu0   = xmumat/ro0(iphas)
326
 
!     Condutivite thermique moleculaire (lambda en W.m-1.K-1)
327
 
xlamb0 = (xnu0/prand0)*cp0(iphas)*ro0(iphas)
328
 
 
329
 
 
330
 
! --- Capacite thermique fictive pour le solide
331
 
 
332
 
!     On cherche un etat stationnaire comme limite d'un transitoire
333
 
!       Les equations portant sur la temperature T des solides
334
 
!         comportent un terme en CP dT/dt qui tend vers zero a
335
 
!         convergence. La valeur de CP est utilisee pour accelerer
336
 
!         la convergence : elle correspond a l'utilisation d'un
337
 
!         pas de temps adapte pour la temperature des solides (comme
338
 
!         on pourrait le faire avec CDTVAR)
339
 
 
340
 
cpcont = 1.d0
341
 
cpparo = 1.d0
342
 
 
343
 
 
344
 
!===============================================================================
345
 
! 2. EROSION DES PANACHES DE CONVECTION NATURELLE
346
 
!===============================================================================
347
 
 
348
 
!     Le but est de reporter dans la zone de ciel une partie de la
349
 
!       puissance transmise a l'air dans la zone de stockage.
350
 
!       Ceci permet de simuler l'effet des panaches, le debit
351
 
!       enthalpique de ceux-ci etant prescrit par l'utilisateur.
352
 
 
353
 
!     Calcul de DHS, DHE
354
 
 
355
 
! --- Debit enthalpique de convection naturelle calcule (initialisation)
356
 
dhs = 0.d0
357
 
 
358
 
! --- Debit enthalpique de convection naturelle impose pour les panaches
359
 
 
360
 
!     Nul par defaut (pas de modelisation des panaches)
361
 
dhe = 0.d0
362
 
!     Calcul, si modelisation des panaches
363
 
if(imdcnt.eq.1) then
364
 
  dhe = dhpcnt/frdtra
365
 
endif
366
 
 
367
 
!===============================================================================
368
 
! 2. RAYONNEMENT
369
 
!===============================================================================
370
 
 
371
 
! --- Flux d'un conteneur (puissance ramenee a la surface en W.m-2)
372
 
 
373
 
!       En cas de stockage en alveole, il y a un plenum inferieur libre
374
 
if(ialveo.eq.1) then
375
 
  phi0 = puicon/((hreso-hplen)*dmcont*pi)
376
 
else
377
 
  phi0 = puicon/(hreso*dmcont*pi)
378
 
endif
379
 
 
380
 
! --- Facteur de forme (alveole ou mur) / conteneur
381
 
 
382
 
!     En alveole, le conteneur voit l'alveole entiere et c'est tout
383
 
if(ialveo.eq.1) then
384
 
 
385
 
  fform = 1.d0
386
 
 
387
 
!     Sans alveole, le facteur de forme depend de la structure du reseau
388
 
else
389
 
 
390
 
!       Pour un reseau en ligne, considerer une rangee suffit
391
 
!         (valable jusqu'a PTRRES/DMCONT=4)
392
 
  fform =  sqrt( sqrt(dmcont/ptrres) + (dmcont/ptrres) )
393
 
 
394
 
!       Pour un reseau en quinconce (reseau a pas triangulaire), on
395
 
!         complete par la contribution d'une seconde rangee
396
 
  if(iconlg.eq.0) then
397
 
    fform = max( (fform + 0.22d0), 1.d0 )
398
 
  endif
399
 
 
400
 
endif
401
 
 
402
 
 
403
 
!  --- Emissivite equivalente (conteneur / mur)
404
 
emis = 1.d0/(1.d0/emicon + 1.d0/emimur)
405
 
 
406
 
!  --- Densite de surface d'echange (par unite de hauteur)
407
 
xxav = pi*dmcont/(ptrres*plgres)
408
 
 
409
 
 
410
 
!  --- Connectivite pour le rayonnement colis / plafond
411
 
!                et pour la modelisation des panaches
412
 
!        Une connectivite est necessaire pour relier la couche
413
 
!          superieure de la zone occupee par des colis et les
414
 
!          mailles situees directement sous le plafond
415
 
 
416
 
!     Cote du haut des mailles du haut des conteneurs
417
 
hplus = hreso
418
 
!     Cote du bas des mailles du haut des conteneurs
419
 
hmoin = hreso - epchel
420
 
!     Cote du bas des mailles touchant le plafond
421
 
hmax  = epchel*nchest - epchel
422
 
 
423
 
!     Calcul de la connectivite
424
 
!       Attention, c'est quadratique ...
425
 
!       mais on ne le fait qu'au premier passage
426
 
 
427
 
!     Calcul uniquement si pas deja fait et si utile
428
 
!       (utile si les colis ne touchent pas le plafond)
429
 
if( (icnrok.eq.0).and.(hreso.lt.epchel*nchest) ) then
430
 
!     On repere les mailles de plafond dans la zone de stockage
431
 
  do iel = 1, ncel
432
 
    ifml  = ifmcel(iel   )
433
 
    icoul = iprfml(ifml,1)
434
 
    if(icoul.eq.icmtst.and.xyzcen(3,iel).gt.hmax)then
435
 
!     On repere la maille du haut des colis situee dessous
436
 
      do jel = 1, ncel
437
 
        difx = abs(xyzcen(1,jel)-xyzcen(1,iel))
438
 
        dify = abs(xyzcen(2,jel)-xyzcen(2,iel))
439
 
        if(difx.lt.ptrres*1.d-2.and.dify.lt.plgres*1.d-2) then
440
 
          if ( xyzcen(3,jel).lt.hplus.and.                        &
441
 
               xyzcen(3,jel).gt.hmoin ) then
442
 
            iconra(iel) = jel
443
 
          endif
444
 
        endif
445
 
      enddo
446
 
    endif
447
 
  enddo
448
 
!     On indique que la connectivite a ete calculee
449
 
  icnrok = 1
450
 
endif
451
 
 
452
 
!===============================================================================
453
 
 
454
 
!  ON VA COMPLETER CI-DESSOUS LES TERMES SOURCES
455
 
!     - PUREMENT EXPLICITES (PUISSANCE)
456
 
!     - PARTIELLEMENT EXPLICITES ET IMPLICITES (ECHANGES)
457
 
 
458
 
 
459
 
!  ON UTILISE LE MODELE DE TERME SOURCE SUIVANT, POUR UN SCALAIRE F :
460
 
 
461
 
!                             S = A * F + B
462
 
 
463
 
!     APPARAISSANT DANS LES EQUATIONS DE TRANSPORT DE F SOUS LA FORME :
464
 
 
465
 
!                       RHO VOLUME D(F)/Dt = VOLUME*S
466
 
 
467
 
 
468
 
!   CE TERME A UNE PARTIE QU'ON VEUT IMPLICITER         : A
469
 
!           ET UNE PARTIE QU'ON VA TRAITER EN EXPLICITE : B
470
 
 
471
 
 
472
 
!   PAR EXEMPLE SI ON A :
473
 
!     A = - RHO / TAUF
474
 
!     B =   RHO * PRODF
475
 
!        AVEC
476
 
!     TAUF   = 10.D0  [secondes  ] (TEMPS DE DISSIPATION DE F)
477
 
!     PRODF  = 100.D0 [variable/s] (PRODUCTION DE F PAR UNITE DE TEMPS)
478
 
 
479
 
!   ON A ALORS
480
 
!     CRVIMP(IEL) = VOLUME(IEL)* A = - VOLUME(IEL) (RHO / TAUF )
481
 
!     CRVEXP(IEL) = VOLUME(IEL)* B =   VOLUME(IEL) (RHO * PRODF)
482
 
 
483
 
!===============================================================================
484
 
 
485
 
! === Initialisation
486
 
! ===========================================================
487
 
 
488
 
! --- Compteurs
489
 
!       pour moyenne et affichage, mais pas seulement
490
 
 
491
 
!     Coefficient d'echange moyen de convection forcee
492
 
hmrey = 0.d0
493
 
!     Coefficient d'echange moyen de convection naturelle
494
 
hmraz = 0.d0
495
 
!     Coefficient d'echange convectif moyen efficace
496
 
hm0   = 0.d0
497
 
!     Coefficient d'echange radiatif moyen pour les parois
498
 
hmray = 0.d0
499
 
!     Compteur pour moyenne
500
 
hnb   = 0.d0
501
 
!     Puissance totale sans panache
502
 
ptot  = 0.d0
503
 
!     Puissance totale avec panache
504
 
ptcn  = 0.d0
505
 
 
506
 
!     Gravite
507
 
gmat = sqrt(gx**2+gy**2+gz**2)
508
 
 
509
 
 
510
 
! === Boucle sur les cellules et selection selon couleur
511
 
! ===========================================================
512
 
 
513
 
do iel = 1, ncel
514
 
 
515
 
!     Vitesse horizontale
516
 
  un0 = sqrt(rtp(iel,iu(iphas))**2+rtp(iel,iv(iphas))**2)
517
 
 
518
 
!     Temperature de reference Kelvin
519
 
  tbulk = (tinit+rtp(iel,isca(itaamt)))*0.5d0 + tkelvi
520
 
 
521
 
!     Raleigh/DeltaT = FRAGH = g*H**3*beta/(nu a)
522
 
!       . g    = GMAT
523
 
!       . H    = XYZCEN(3,IEL) : on suppose que le sol est a z=0
524
 
!       . beta = 1/TBULK
525
 
!       . nu   = XNU0
526
 
!       . a    = nu/Pr = XNU0/PRAND0
527
 
 
528
 
  fragh = gmat*xyzcen(3,iel)**3*prand0/(tbulk*xnu0**2)
529
 
 
530
 
!     Couleur pour selection de la zone de stockage
531
 
  ifml  = ifmcel(iel   )
532
 
  icoul = iprfml(ifml,1)
533
 
 
534
 
  if(icoul.eq.icmtst) then
535
 
 
536
 
 
537
 
! === Calcul du coefficient d'echange par rayonnement
538
 
! ===========================================================
539
 
 
540
 
!     Coeff d'echange = flux / (S (Tcolis-Tparoi))
541
 
!       = Sigma Epsilon FactForme (Tcolis+Tparoi)(Tcolis**2+Tparoi**2)
542
 
 
543
 
!     La paroi consideree est le plafond ou les murs ou les alveoles
544
 
 
545
 
! --- Zone de ciel
546
 
!       On utilise la connectivite calculee precedemment
547
 
    if(icnrok.eq.1.and.xyzcen(3,iel).gt.hmax) then
548
 
 
549
 
!     En presence d'alveoles,
550
 
!       le plafond echange avec les alveoles face a face
551
 
!       (seule l'alveole directement en face est vue,
552
 
!        l'alveole voit un disque de plafond)
553
 
!       La temperature du plafond et des alveoles est representee
554
 
!       par le meme scalaire
555
 
!       Ce coefficient est utilise pour le calcul des termes sources
556
 
!       permettant de determiner la temperature de peau des parois
557
 
!       (plafond, mur, alveoles)
558
 
!       On le verra aussi intervenir pour un calcul de temperature
559
 
!       de colis, mais dans la zone de ciel ou la temperature de colis
560
 
!       n'a pas de signification (et comme il n'y a pas de diffusion
561
 
!       en alveoles, c'est ok)
562
 
      if(ialveo.eq.1) then
563
 
        hray = stephn*emimur                                      &
564
 
             * (0.25d0*pi*dmcont**2)/(ptrres*plgres)              &
565
 
             * ( rtp(iconra(iel),isca(itppmt))+tkelvi             &
566
 
                +rtp(       iel ,isca(itppmt))+tkelvi)            &
567
 
             * ((rtp(iconra(iel),isca(itppmt))+tkelvi)**2         &
568
 
               +(rtp(       iel ,isca(itppmt))+tkelvi)**2)
569
 
 
570
 
!     Sans alveoles,
571
 
!       le plafond echange avec les colis face a face
572
 
!       Ce coefficient sera utilise pour le calcul des termes sources
573
 
!       permettant de determiner la temperature de peau du plafond
574
 
      else
575
 
        hray = stephn*emis                                        &
576
 
             * (0.25d0*pi*dmcont**2)/(ptrres*plgres)              &
577
 
             * ( rtp(iconra(iel),isca(itpcmt))+tkelvi             &
578
 
                +rtp(       iel ,isca(itppmt))+tkelvi)            &
579
 
             * ((rtp(iconra(iel),isca(itpcmt))+tkelvi)**2         &
580
 
               +(rtp(       iel ,isca(itppmt))+tkelvi)**2)
581
 
      endif
582
 
 
583
 
! --- Zones de murs ou d'alveoles
584
 
!       les murs et les alveoles echangent avec les colis
585
 
!       le facteur de forme calcule precedemment s'applique
586
 
!       il prend en compte le cas avec ou sans alveoles
587
 
!       Ce coefficient sera utilise pour le calcul des termes sources
588
 
!       permettant de determiner
589
 
!         - la temperature de peau des parois
590
 
!         - la temperature de peau des colis, en alveole (FFORM=1)
591
 
 
592
 
    else
593
 
      hray = stephn*emis                                          &
594
 
           * fform                                                &
595
 
           * ( rtp(iel,isca(itpcmt))+tkelvi                       &
596
 
              +rtp(iel,isca(itppmt))+tkelvi)                      &
597
 
           * ((rtp(iel,isca(itpcmt))+tkelvi)**2                   &
598
 
             +(rtp(iel,isca(itppmt))+tkelvi)**2)
599
 
    endif
600
 
 
601
 
 
602
 
! === Puissance transmise a l'air et a la peau des colis par conduction
603
 
! ===========================================================
604
 
 
605
 
!     Pour l'air :
606
 
!       en stationnaire, toute la puissance des colis se retrouve tot ou
607
 
!       tard dans l'air (elle peut transiter par la peau des colis,
608
 
!       par les alveoles, par les murs, mais elle est extraite par
609
 
!       l'air).
610
 
 
611
 
!     Pour la peau des colis :
612
 
!       toute la puissance des colis est transmise par conduction a
613
 
!       la peau des colis.
614
 
 
615
 
 
616
 
    if (iscal.eq.itaamt.or.iscal.eq.itpcmt) then
617
 
 
618
 
! --- Puissance des colis par unite de surface au sol
619
 
!       On prend en compte la carte de remplissage ensuite
620
 
 
621
 
!     On doit fournit un CRVEXP en Watt/CP, soit de la dimension de
622
 
!       PUICON/CP.
623
 
 
624
 
!     Comme la puissance integree sur la hauteur du colis doit etre
625
 
!       la puissance totale du colis, PUICON, on doit imposer dans
626
 
!       chaque maille une fraction de PUICON.
627
 
!     Pour une repartition uniforme en hauteur, cette fraction est
628
 
!       la hauteur de la maille VOLUME/(PTRRES*PLGRES) divisee par
629
 
!       la hauteur totale du colis.
630
 
!     Les repartitions non uniformes en hauteur sont traitees par
631
 
!       les cartes de repartition en altitude, plus bas.
632
 
!     On impose donc CRVEXP = PUICON/CP * VOLUME/(PTRRES*PLGRES)
633
 
!       puis on divise par la hauteur de colis (plus bas, lors du
634
 
!       traitement des cartes de repartition en altitude).
635
 
 
636
 
!     Noter que comme on multiplie par VOLUME/(PTRRES*PLGRES) ici
637
 
!       et qu'on divise plus bas par EPCHEL, on pourrait simplifier,
638
 
!       puisque, a priori, EPCHEL = VOLUME/(PTRRES*PLGRES)
639
 
 
640
 
 
641
 
!     Air ambiant
642
 
      if (iscal.eq.itaamt) then
643
 
        crvexp(iel) = puicon/cp0(iphas)                           &
644
 
             * volume(iel)/(ptrres*plgres)
645
 
!     Peau des colis
646
 
      elseif(iscal.eq.itpcmt) then
647
 
        crvexp(iel) = puicon/cpcont                               &
648
 
             * volume(iel)/(ptrres*plgres)
649
 
      endif
650
 
 
651
 
! --- Prise en compte des cartes de remplissage (adimensionnelles)
652
 
!     Les entrepots ne sont pas uniformement remplis.
653
 
 
654
 
!     Pour repr�senter la carte d'encombrement des lignes de colis, on
655
 
!       repere les lignes occupees (coordonnee X, indicateur ILIGNE)
656
 
!       par leurs bornes XLIMIN et XLIMAX qui sont stockees dans VIZCAR
657
 
!       (avec l'indicateur ICPUIS faisant reference a la puissance)
658
 
      nzones = nzocar(iligne,icpuis)
659
 
      do izone = 1, nzones
660
 
        xlimin = ptrres*vizcar(1,izone,iligne,icpuis)
661
 
        xlimax = ptrres*vizcar(2,izone,iligne,icpuis)
662
 
        if ((xyzcen(1,iel).ge.xlimin).and.                        &
663
 
             (xyzcen(1,iel).lt.xlimax)) then
664
 
          crvexp(iel) = crvexp(iel) * vcarth(izone,iligne)
665
 
        endif
666
 
      enddo
667
 
 
668
 
!     Pour repr�senter la carte d'encombrement des rangees de colis, on
669
 
!       repere les rangees occupees (coordonnee Y, indicateur IRANGE)
670
 
!       par leurs bornes YRGMIN et YRGMAX qui sont stockees dans VIZCAR
671
 
!       (avec l'indicateur ICPUIS faisant reference a la puissance)
672
 
      nzones = nzocar(irange,icpuis)
673
 
      do izone = 1, nzones
674
 
        yrgmin = plgres*vizcar(1,izone,irange,icpuis)
675
 
        yrgmax = plgres*vizcar(2,izone,irange,icpuis)
676
 
        if ((xyzcen(2,iel).ge.yrgmin).and.                        &
677
 
             (xyzcen(2,iel).lt.yrgmax)) then
678
 
          crvexp(iel) = crvexp(iel) * vcarth(izone,irange)
679
 
        endif
680
 
      enddo
681
 
 
682
 
!     Pour repr�senter la carte d'encombrement en hauteur de colis, on
683
 
!       repere les couches occupees (coordonnee Z, indicateur IALTIT)
684
 
!       par leurs bornes ZALMIN et ZALMAX qui sont stockees dans VIZCAR
685
 
!       (avec l'indicateur ICPUIS faisant reference a la puissance)
686
 
!     SQZ permet d'integrer VCARTH sur les zones afin que
687
 
!       de normer VCARTH de telle sorte que la multiplication par
688
 
!       VCARTH/SQZ permette de retrouver, en integrant sur la hauteur,
689
 
!       toute la puissance du colis.
690
 
      sqz = 0.d0
691
 
      nzones = nzocar(ialtit,icpuis)
692
 
      do izone = 1, nzones
693
 
        zalmin = epchel*vizcar(1,izone,ialtit,icpuis)
694
 
        zalmax = epchel*vizcar(2,izone,ialtit,icpuis)
695
 
        sqz = sqz + vcarth(izone,ialtit) *                        &
696
 
             ( vizcar(2,izone,ialtit,icpuis) -                    &
697
 
             vizcar(1,izone,ialtit,icpuis) )
698
 
      enddo
699
 
      do izone = 1, nzones
700
 
        zalmin = epchel*vizcar(1,izone,ialtit,icpuis)
701
 
        zalmax = epchel*vizcar(2,izone,ialtit,icpuis)
702
 
        if ((xyzcen(3,iel).ge.zalmin).and.                        &
703
 
             (xyzcen(3,iel).lt.zalmax)) then
704
 
          crvexp(iel) = crvexp(iel)*                              &
705
 
               vcarth(izone,ialtit)/(sqz*epchel)
706
 
        endif
707
 
      enddo
708
 
 
709
 
!       On calcule la puissance totale injectee (Watt)
710
 
 
711
 
!     Air ambiant
712
 
      if (iscal.eq.itaamt) then
713
 
        ptot = ptot + crvexp(iel)*cp0(iphas)
714
 
      endif
715
 
 
716
 
    endif
717
 
 
718
 
 
719
 
! === Echange convectif conteneur / air et radiatif conteneur / mur
720
 
! ===========================================================
721
 
 
722
 
!     Pour la determination de la temperature de peau des colis
723
 
 
724
 
!     On calcule le coefficient d'echange convectif efficace selon
725
 
!       la configuration (convection naturelle, forcee, mixte).
726
 
!     On ajoute le terme de rayonnement a partir du coefficient
727
 
!       equivalent HRAY calcule precedemment.
728
 
 
729
 
    if (iscal.eq.itpcmt) then
730
 
 
731
 
! --- Convection naturelle verticale le long des colis
732
 
!       Van Vliet&Ross et Vliet&Liu
733
 
 
734
 
!     Rayleigh a partir du flux d'un conteneur
735
 
      raleiz = fragh*phi0*xyzcen(3,iel)/xlamb0
736
 
!     Nusselt associe
737
 
      xnsraz = 0.17d0 * raleiz**0.25d0
738
 
      xnsraz = max(xnsraz,(0.6d0 * raleiz**0.2d0))
739
 
!     Coefficient d'echange de convection naturelle       conteneur / air
740
 
!       Hypothese que le zero est a la base des colis
741
 
      hraz   = xnsraz*xlamb0/xyzcen(3,iel)
742
 
!     Coefficient d'echange moyen de convection naturelle
743
 
      hmraz  = hmraz + hraz
744
 
 
745
 
 
746
 
! --- Convection forcee transversale dans le reseau de colis
747
 
!       Zukauskas
748
 
 
749
 
!     Reynolds (vitesse horizontale gap, diametre
750
 
      reynol = (un0*ptrres/(ptrres-dmcont)) * dmcont / xnu0
751
 
!     Nusselt associe
752
 
      xnsrey = 0.35d0*(ptrres/plgres)**0.2d0                      &
753
 
           * 0.89d0*reynol**0.6d0
754
 
!     Coefficient d'echange de convection forcee          conteneur / air
755
 
      hrey   = xnsrey*xlamb0/dmcont
756
 
!     Coefficient d'echange moyen de convection forcee
757
 
      hmrey = hmrey + hrey
758
 
 
759
 
! --- Compteur pour les moyennes
760
 
 
761
 
      hnb = hnb + 1.d0
762
 
 
763
 
 
764
 
! --- Prise en compte du type de reseau (en ligne ou en quinconce)
765
 
!       et de la presence d'alveoles
766
 
 
767
 
!   - Alveoles
768
 
!       convection naturelle uniquement
769
 
 
770
 
      if(ialveo.eq.1) then
771
 
 
772
 
!     Nusselt associe
773
 
        xnsalv = 0.17d0 * raleiz**0.25d0
774
 
!     Coefficient d'echange de convection naturelle       conteneur / air
775
 
!       Hypothese que le zero est a la base des colis
776
 
        h0 = xnsalv*xlamb0/xyzcen(3,iel)
777
 
 
778
 
 
779
 
!   - Reseau en ligne (pas carre) sans alveoles
780
 
!       2/3 conv.force (face+cotes) + 1/3 conv.nat (sillage arriere)
781
 
 
782
 
      elseif(iconlg.eq.1)then
783
 
        h0 = (hraz+2.d0*hrey)/3.d0
784
 
 
785
 
!   - Reseau en quinconce (pas triangulaire) sans alveoles
786
 
!       conv.force (pas de protection du colis amont)
787
 
 
788
 
      else
789
 
        h0 = hrey
790
 
      endif
791
 
 
792
 
!   - Sans alveoles
793
 
!       le coefficient d'echange est toujours au moins egal a celui
794
 
!       de la convection naturelle
795
 
      if(ialveo.ne.1) then
796
 
        h0 = max(hraz,h0)
797
 
      endif
798
 
 
799
 
!    - Coefficient d'echange convectif moyen efficace
800
 
      hm0 = hm0 + h0
801
 
 
802
 
 
803
 
! --- Ajout du terme conductif et du terme radiatif
804
 
 
805
 
!       CRVEXP contient deja la puissance des colis transmise
806
 
!         par conduction
807
 
!       On veut ajouter dans CRVEXP des Watt/CP soit H0*Surface*T/CP.
808
 
!         On a ici XXAV*VOLUME = PI*DMCONT*EPCHEL, qui est la surface
809
 
!         d'�change dans la cellule.
810
 
 
811
 
 
812
 
!   - Alveoles
813
 
!       echange par convection naturelle et
814
 
!       les colis rayonnent avec le plafond et les alveoles
815
 
 
816
 
      if(ialveo.eq.1) then
817
 
 
818
 
        crvexp(iel) = crvexp(iel)                                 &
819
 
             +h0  *(xxav*volume(iel))                             &
820
 
                  *rtp(iel,isca(itaamt))/cpcont                   &
821
 
             +hray*(xxav*volume(iel))                             &
822
 
                  *rtp(iel,isca(itppmt))/cpcont
823
 
        crvimp(iel) = -(h0+hray)*(xxav*volume(iel))/cpcont
824
 
 
825
 
 
826
 
!   - Sans alveole
827
 
!       echange par convection mixte seul
828
 
!       les colis ne rayonnent qu'entre eux, et ceci est pris en compte
829
 
!         par diffusion
830
 
 
831
 
      else
832
 
 
833
 
        crvexp(iel) = crvexp(iel)                                 &
834
 
             +h0  *(xxav*volume(iel))                             &
835
 
                  *rtp(iel,isca(itaamt))/cpcont
836
 
        crvimp(iel) = - h0      *(xxav*volume(iel))/cpcont
837
 
 
838
 
      endif
839
 
 
840
 
    endif
841
 
 
842
 
 
843
 
 
844
 
! === Echange convectif mur / air et radiatif mur / conteneur
845
 
! ===========================================================
846
 
 
847
 
!     Pour la determination de la temperature de peau des parois
848
 
!       (murs et alveoles)
849
 
 
850
 
    if (iscal.eq.itppmt) then
851
 
 
852
 
! --- Echange par convection paroi / air
853
 
 
854
 
!   - Alveoles
855
 
!       convection naturelle uniquement
856
 
 
857
 
      if(ialveo.eq.1) then
858
 
 
859
 
!     Rayleigh a partir du flux d'un conteneur
860
 
        raleiz = fragh*phi0*xyzcen(3,iel)/xlamb0
861
 
!     Nusselt associe
862
 
        xnsalv = 0.17d0 * raleiz**0.25d0
863
 
!     Coefficient d'echange de convection naturelle       alveole / air
864
 
!       Hypothese que le zero est a la base des colis
865
 
        hraz = xnsalv*xlamb0/xyzcen(3,iel)
866
 
!     Coefficient d'echange efficace alveole / air
867
 
!       Convection naturelle seule
868
 
        h0 = hraz
869
 
 
870
 
 
871
 
!   - Sans alveole
872
 
!       convection naturelle (Mac Adams) ou convection forcee
873
 
 
874
 
      else
875
 
 
876
 
!    .Convection naturelle
877
 
 
878
 
!     Rayleigh a partir de l'ecart de temperature paroi - air
879
 
        raleiz =                                                  &
880
 
             fragh*(rtp(iel,isca(itppmt))-rtp(iel,isca(itaamt)))
881
 
!     Nusselt associe (Mac Adams)
882
 
        xnsraz = 0.12d0*(abs(raleiz))**untier
883
 
!     Coefficient d'echange de convection naturelle       paroi / air
884
 
        hraz   = xnsraz*xlamb0/xyzcen(3,iel)
885
 
!     Coefficient d'echange moyen de convection naturelle
886
 
        hmraz  = hmraz + hraz
887
 
 
888
 
 
889
 
!    .Convection forcee
890
 
 
891
 
!     Espace mur reseau
892
 
        dhmur  = ptrres-dmcont
893
 
!     Reynolds associ� (vitesse horizontale hors reseau, pres des murs)
894
 
        reynol = dhmur*un0/xnu0
895
 
!     Nusselt associe (Colburn)
896
 
        xnsrey = 0.023d0*reynol**0.8d0*prand0**untier
897
 
!     Coefficient d'echange de convection forcee          paroi / air
898
 
        hrey   = xnsrey*xlamb0/dhmur
899
 
!     Coeficient d'echange moyen de convection forcee
900
 
        hmrey  = hmrey + hrey
901
 
 
902
 
 
903
 
!     .Coefficient d'echange efficace (le plus efficace est dominant)
904
 
        h0     = max(hraz,hrey)
905
 
 
906
 
      endif
907
 
 
908
 
 
909
 
!   - Coefficient d'echange convectif moyen efficace
910
 
      hm0   = hm0 + h0
911
 
!   - Coefficient d'echange moyen de convection naturelle
912
 
      hmraz = hmraz + hraz
913
 
!   - Compteur pour les moyennes
914
 
      hnb   = hnb + 1.d0
915
 
 
916
 
 
917
 
! --- Ajout du terme conductif et du terme radiatif
918
 
 
919
 
!       On veut ajouter dans CRVEXP des Watt/CP soit H0*Surface*T/CP.
920
 
!         On calcule H0*VOLUME*T/CP et on divise par une distance
921
 
!         plus bas. On pourrait faire les choses plus directement.
922
 
 
923
 
!       Les echanges radiatifs sont les echanges avec les colis.
924
 
 
925
 
 
926
 
!   - Coefficient d'echange radiatif moyen pour les parois
927
 
      hmray = hmray + hray
928
 
 
929
 
 
930
 
!   - Echanges convectifs
931
 
      crvexp(iel) = h0*volume(iel)*rtp(iel,isca(itaamt))/cpparo
932
 
      crvimp(iel) =-h0*volume(iel)/cpparo
933
 
 
934
 
 
935
 
!   - Echange radiatif dans le ciel
936
 
      if(icnrok.eq.1.and.xyzcen(3,iel).gt.hmax) then
937
 
 
938
 
!    .Alveoles : le plafond echange avec les alveoles
939
 
        if(ialveo.eq.1) then
940
 
          crvexp(iel) = crvexp(iel)                               &
941
 
               +hray*volume(iel)                                  &
942
 
               *rtp(iconra(iel),isca(itppmt))/cpparo
943
 
 
944
 
!    .Sans alveole : le plafond echange avec le dessus des colis
945
 
        else
946
 
          crvexp(iel) = crvexp(iel)                               &
947
 
               +hray*volume(iel)                                  &
948
 
               *rtp(iconra(iel),isca(itpcmt))/cpparo
949
 
        endif
950
 
 
951
 
 
952
 
!   - Echange radiatif avec les colis pour les autres parois
953
 
!       (murs et alveoles)
954
 
      else
955
 
        crvexp(iel) = crvexp(iel)                                 &
956
 
             +hray*volume(iel)*rtp(iel,isca(itpcmt))/cpparo
957
 
 
958
 
      endif
959
 
 
960
 
 
961
 
!   - La partie implicite est la meme pour tous
962
 
      crvimp(iel) = crvimp(iel) - hray*volume(iel)/cpparo
963
 
 
964
 
 
965
 
!   - On divise par une distance
966
 
!       Plus exactement, on multiplie par la surface d'echange et
967
 
!       on divise par le volume (on aurait donc pu eviter de multiplier
968
 
!       par le volume).
969
 
 
970
 
 
971
 
!    .Au dessus ou au dessous des colis, le plafond ou le sol echangent
972
 
!       sur une surface de PTRRES*PLGRES
973
 
!       dans des cellules de volume EPCHEL*PTRRES*PLGRES.
974
 
!       Le rapport Surface/Volume est donc 1/EPCHEL
975
 
 
976
 
      if(xyzcen(3,iel).lt.epchel.or.xyzcen(3,iel).gt.hmax) then
977
 
        crvexp(iel) = crvexp(iel)/epchel
978
 
        crvimp(iel) = crvimp(iel)/epchel
979
 
 
980
 
 
981
 
!    .En stockage en alveoles, les alveoles echangent
982
 
!       sur une surface de PI*DHALVE*EPCHEL
983
 
!       dans des cellules de volume EPCHEL*PTRRES*PLGRES.
984
 
!       Le rapport Surface/Volume est donc PI*DHALVE/(PTRRES*PLGRES)
985
 
!     A priori DHALVE est positif, mais on laisse le ABS (pour le
986
 
!       cas, ou il serait defini comme la difference entre le
987
 
!       diametre des colis et des alveoles)
988
 
 
989
 
      elseif(ialveo.eq.1) then
990
 
        crvexp(iel) =                                             &
991
 
             crvexp(iel)*pi*abs(dhalve)/(ptrres*plgres)
992
 
        crvimp(iel) =                                             &
993
 
             crvimp(iel)*pi*abs(dhalve)/(ptrres*plgres)
994
 
 
995
 
 
996
 
!    .Sans alveoles, les murs echangent
997
 
!       sur une surface de PLGRES*EPCHEL
998
 
!       dans des cellules de volume EPCHEL*PTRRES*PLGRES.
999
 
!       Le rapport Surface/Volume est donc 1/PTRRES
1000
 
!       La ou il n'y a pas de mur, la temperature de mur ne signifie
1001
 
!       rien, mais on peut quand meme l'y calculer (il n'y a pas de
1002
 
!       phenomene de diffusion, tout est local)
1003
 
 
1004
 
      else
1005
 
        crvexp(iel) = crvexp(iel)/ptrres
1006
 
        crvimp(iel) = crvimp(iel)/ptrres
1007
 
 
1008
 
      endif
1009
 
 
1010
 
    endif
1011
 
 
1012
 
 
1013
 
!     Fin du test sur la zone de stockage
1014
 
  endif
1015
 
enddo
1016
 
 
1017
 
 
1018
 
! === Pour l'air, sans modelisation des panaches de convection naturelle
1019
 
! ===========================================================
1020
 
 
1021
 
if(iscal.eq.itaamt) then
1022
 
 
1023
 
  if(imdcnt.eq.0)then
1024
 
 
1025
 
! --- Calcul du debit enthalpique de convection naturelle
1026
 
!       sortant en haut de la zone des colis : rho*W*S*Cp*DeltaT
1027
 
!       (il aurait ete plus correct d'utiliser FLUMAS)
1028
 
 
1029
 
    do iel = 1, ncel
1030
 
 
1031
 
      ifml  = ifmcel(iel   )
1032
 
      icoul = iprfml(ifml,1)
1033
 
      if(icoul.eq.icmtst) then
1034
 
 
1035
 
        if(icnrok.eq.1.and.xyzcen(3,iel).gt.hmax) then
1036
 
          jel = iconra(iel)
1037
 
          dhs = dhs +                                             &
1038
 
               propce(jel,ipcrom)*rtp(jel,iw(iphas))              &
1039
 
               *(volume(jel)/epchel)                              &
1040
 
               *cp0(iphas)*(rtp(jel,isca(itaamt))-tinit)
1041
 
        endif
1042
 
 
1043
 
      endif
1044
 
    enddo
1045
 
 
1046
 
 
1047
 
! --- Calcul de la puissance en kW,
1048
 
!       avec correction si representation 3D en 2D
1049
 
 
1050
 
    puitot = ptot*frdtra*1.d-3
1051
 
 
1052
 
 
1053
 
! --- Calcul du debit enthalpique de convection naturelle en kW,
1054
 
!       avec correction si representation 3D en 2D
1055
 
 
1056
 
    debcon = dhs*frdtra*1.d-3
1057
 
 
1058
 
! --- Impression de PUITOT et DEBCON
1059
 
    if (irangp.le.0.and.modntl.eq.0) then
1060
 
      write(nfecra,2002) puitot, debcon
1061
 
    endif
1062
 
 
1063
 
 
1064
 
! === Pour l'air, modelisation des panaches de convection naturelle
1065
 
! ===========================================================
1066
 
 
1067
 
  elseif(imdcnt.eq.1) then
1068
 
 
1069
 
 
1070
 
! --- Le debit enthalpique des panaches ne peut pas etre plus grand
1071
 
!       que la puissance totale disponible ; sinon, c'est qu'il
1072
 
!       a ete mal predit ou que la puissance a ete mal imposee :
1073
 
!       les cartes ou la puissance d'un conteneur peuvent etre en
1074
 
!       cause
1075
 
!     On fait ce test ici, dans la mesure ou PTOT est calcule juste
1076
 
!       au dessus.
1077
 
 
1078
 
    if(dhpcnt.ge.ptot*frdtra)then
1079
 
      write(nfecra,9001) imdcnt, dhpcnt, ptot*frdtra,             &
1080
 
           dhpcnt, ptot*frdtra
1081
 
      call csexit (1)
1082
 
      !==========
1083
 
    endif
1084
 
 
1085
 
! --- Modification du terme source pour la modelisation des panaches
1086
 
 
1087
 
!       On redistribue la puissance transmise a l'air uniquement :
1088
 
!         une fraction de la puissance totale est retiree de la zone
1089
 
!         de stockage et injectee directement dans le ciel au dessus
1090
 
!         des conteneurs. La repartition se fait a priori de maniere
1091
 
!         homogene (le panache de tous les conteneurs presents recoit
1092
 
!         la meme fraction de la puissance totale), puis on applique
1093
 
!         la carte de repartition thermique XY.
1094
 
 
1095
 
!   - Calcul des grandeurs independantes de la maille
1096
 
 
1097
 
!     Fraction de la puissance totale a conserver
1098
 
!       dans la zone de stockage
1099
 
    fptsto = 1.d0-dhe/ptot
1100
 
 
1101
 
!     Puissance Volumique a ajouter dans la zone de Ciel au dessus
1102
 
!       d'un Conteneur de puissance PUICON, divise par CP
1103
 
    pvccsc = ( (dhe/ptot)*puicon / (hercnt*ptrres*plgres) )       &
1104
 
         / cp0(iphas)
1105
 
 
1106
 
    do iel = 1, ncel
1107
 
 
1108
 
      ifml  = ifmcel(iel   )
1109
 
      icoul = iprfml(ifml,1)
1110
 
 
1111
 
      if(icoul.eq.icmtst) then
1112
 
 
1113
 
!   - Reduction dans la zone de stockage
1114
 
 
1115
 
!     La puissance est reduite pour tous les conteneurs de maniere
1116
 
!       homogene, en fonction du rapport entre le debit enthalpique
1117
 
!       des panaches impose par l'utilisateur et la puissance totale
1118
 
!       presente dans l'entrepot. Les cartes de repartition thermique
1119
 
!       ayant deja ete appliquees a CRVEXP, c'est bien le rapport de
1120
 
!       reduction (et non pas la valeur absolue de la reduction) qui
1121
 
!       est homogene.
1122
 
 
1123
 
        crvexp(iel) = crvexp(iel)*fptsto
1124
 
 
1125
 
 
1126
 
        if ( xyzcen(3,iel).lt.(hercnt+hreso).and.                 &
1127
 
             xyzcen(3,iel).gt.hreso) then
1128
 
 
1129
 
!   - Report dans la zone des panaches (1/2)
1130
 
 
1131
 
!     Calcul de la puissance a ajouter dans la zone de ciel, au dessus
1132
 
!       d'un conteneur de puissance PUICON (deja divise par CP)
1133
 
 
1134
 
          pccsc = pvccsc * volume(iel)
1135
 
 
1136
 
!     Les entrepots ne sont pas uniformement remplis.
1137
 
!     Application de la carte thermique horizontale XY pour
1138
 
!       selectionner les positions auxquelles un conteneur (ou
1139
 
!       une fraction de conteneur est presente, puisque l'on
1140
 
!       peut avoir des demi-conteneurs, par exemple sur un plan
1141
 
!       de symetrie)
1142
 
 
1143
 
!     Pour repr�senter la carte d'encombrement des lignes de colis, on
1144
 
!       repere les lignes occupees (coordonnee X, indicateur ILIGNE)
1145
 
!       par leurs bornes XLIMIN et XLIMAX qui sont stockees dans VIZCAR
1146
 
!       (avec l'indicateur ICPUIS faisant reference a la puissance)
1147
 
          nzones = nzocar(iligne,icpuis)
1148
 
          do izone = 1, nzones
1149
 
            xlimin = ptrres*vizcar(1,izone,iligne,icpuis)
1150
 
            xlimax = ptrres*vizcar(2,izone,iligne,icpuis)
1151
 
            if ((xyzcen(1,iel).ge.xlimin).and.                    &
1152
 
                 (xyzcen(1,iel).lt.xlimax)) then
1153
 
              pccsc = pccsc *                                     &
1154
 
                   vcarth(izone,iligne)
1155
 
            endif
1156
 
          enddo
1157
 
 
1158
 
!     Pour repr�senter la carte d'encombrement des rangees de colis, on
1159
 
!       repere les rangees occupees (coordonnee Y, indicateur IRANGE)
1160
 
!       par leurs bornes YRGMIN et YRGMAX qui sont stockees dans VIZCAR
1161
 
!       (avec l'indicateur ICPUIS faisant reference a la puissance)
1162
 
          nzones = nzocar(irange,icpuis)
1163
 
          do izone = 1, nzones
1164
 
            yrgmin = plgres*vizcar(1,izone,irange,icpuis)
1165
 
            yrgmax = plgres*vizcar(2,izone,irange,icpuis)
1166
 
            if ((xyzcen(2,iel).ge.yrgmin).and.                    &
1167
 
                 (xyzcen(2,iel).lt.yrgmax)) then
1168
 
              pccsc = pccsc *                                     &
1169
 
                   vcarth(izone,irange)
1170
 
            endif
1171
 
          enddo
1172
 
 
1173
 
!     Ajout de la puissance des panaches a la puissance deja existante
1174
 
 
1175
 
          crvexp(iel) = crvexp(iel) + pccsc
1176
 
 
1177
 
!     Fin si zone de panache
1178
 
        endif
1179
 
 
1180
 
!     Fin si zone de stockage
1181
 
      endif
1182
 
 
1183
 
 
1184
 
! --- Calcul de la puissance totale incluant la modelisation
1185
 
!       des panaches de convection natuelle
1186
 
!       (normalement la puissance totale n'a pas ete modifiee)
1187
 
 
1188
 
      ptcn = ptcn + crvexp(iel)*cp0(iphas)
1189
 
 
1190
 
!     Fin boucle NCEL
1191
 
    enddo
1192
 
 
1193
 
! --- Calcul et impression de la puissance totale incluant la modelisation
1194
 
!       des panaches de convection naturelle
1195
 
!       en kW, avec correction si representation 3D en 2D
1196
 
 
1197
 
    puitot = ptcn*frdtra*1.d-3
1198
 
    if ((irangp.le.0).and.(modntl.eq.0)) then
1199
 
      write(nfecra,2003) puitot
1200
 
    endif
1201
 
 
1202
 
!     Fin si on modelise les panaches
1203
 
  endif
1204
 
!     Fin si scalaire = air ambiant
1205
 
endif
1206
 
 
1207
 
 
1208
 
 
1209
 
! === Impressions pour la temperature de peau des colis et des parois
1210
 
! ===========================================================
1211
 
 
1212
 
! --- Valeur fictive : pas de convection forcee en alveole
1213
 
!     (valeur bizarre, mais uniquement pour affichage)
1214
 
if (ialveo.eq.1) then
1215
 
  hmrey = -1.00001d10
1216
 
endif
1217
 
 
1218
 
! --- Temperature de peau des colis
1219
 
 
1220
 
if(iscal.eq.itpcmt) then
1221
 
 
1222
 
!     Si rien a moyenner, on prend une valeur negative
1223
 
!       (mais les numerateurs seront nuls, normalement)
1224
 
  if (hnb.lt.epzero) then
1225
 
    hnb = -epzero
1226
 
  endif
1227
 
 
1228
 
!     Stockage en common pour impression ulterieure
1229
 
  cfecca = hm0/hnb
1230
 
 
1231
 
!     Impressions des coefficients d'echange moyens
1232
 
  if ((irangp.le.0).and.(modntl.eq.0)) then
1233
 
    write(nfecra,2004) hmrey/hnb, hmraz/hnb, hm0/hnb
1234
 
  endif
1235
 
 
1236
 
endif
1237
 
 
1238
 
! --- Temperature de peau des parois
1239
 
 
1240
 
if(iscal.eq.itppmt) then
1241
 
 
1242
 
!     Si rien a moyenner, on prend une valeur negative
1243
 
!       (mais les numerateurs seront nuls, normalement)
1244
 
  if (hnb.lt.epzero) then
1245
 
    hnb = -epzero
1246
 
  endif
1247
 
 
1248
 
!     Stockage en common pour impression ulterieure
1249
 
  cfecma = hm0/hnb
1250
 
 
1251
 
!     Impressions des coefficients d'echange moyens
1252
 
  if ((irangp.le.0).and.(modntl.eq.0)) then
1253
 
    write(nfecra,2005)                                            &
1254
 
         hmrey/hnb, hmraz/hnb, hmray/hnb, hm0/hnb
1255
 
  endif
1256
 
 
1257
 
endif
1258
 
 
1259
 
 
1260
 
!--------
1261
 
! FORMATS
1262
 
!--------
1263
 
 
1264
 
 1000 format(' TERMES SOURCES MATISSE POUR LA VARIABLE ',A8,/)
1265
 
 
1266
 
 2001 FORMAT (/,3X,'** INFORMATIONS SUR MATISSE, VARIABLE ',A8,/, &
1267
 
          3X,'   -------------------------------------------')
1268
 
 2002 format (                                                          &
1269
 
'mati --------------------------------------------------------',/,&
1270
 
'mati Puissance (sans modelisation panaches)                  ',/,&
1271
 
'mati --------------------------------------                  ',/,&
1272
 
'mati Puissance totale                      ', E12.5  ,'    kW',/,&
1273
 
'mati Debit enthalpique de conv. naturelle  ', E12.5  ,'    kW',/,&
1274
 
'mati --------------------------------------------------------',/)
1275
 
 2003 format (                                                          &
1276
 
'mati --------------------------------------------------------',/,&
1277
 
'mati Puissance (avec modelisation panaches)                  ',/,&
1278
 
'mati --------------------------------------                  ',/,&
1279
 
'mati Puissance totale avec erosion panaches', E12.5  ,'    kW',/,&
1280
 
'mati --------------------------------------------------------',/)
1281
 
 2004 format (                                                          &
1282
 
'mati --------------------------------------------------------',/,&
1283
 
'mati Coeff. d echange moyens (calcul T� de peau des colis)   ',/,&
1284
 
'mati -----------------------------------------------------   ',/,&
1285
 
'mati Convection forcee                   ', E12.5  ,' W/m2/�C',/,&
1286
 
'mati Convection naturelle                ', E12.5  ,' W/m2/�C',/,&
1287
 
'mati Efficace                            ', E12.5  ,' W/m2/�C',/,&
1288
 
'mati --------------------------------------------------------',/)
1289
 
 2005 format (                                                          &
1290
 
'mati --------------------------------------------------------',/,&
1291
 
'mati Coeff. d echange moyens (calcul T� de peau des parois)  ',/,&
1292
 
'mati ------------------------------------------------------  ',/,&
1293
 
'mati Convection forcee                   ', E12.5  ,' W/m2/�C',/,&
1294
 
'mati Convection naturelle                ', E12.5  ,' W/m2/�C',/,&
1295
 
'mati Rayonnement                         ', E12.5  ,' W/m2/�C',/,&
1296
 
'mati Efficace                            ', E12.5  ,' W/m2/�C',/,&
1297
 
'mati --------------------------------------------------------',/)
1298
 
 
1299
 
 9001 format(                                                           &
1300
 
'@                                                            ',/,&
1301
 
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
1302
 
'@                                                            ',/,&
1303
 
'@ @@ ATTENTION : ARRET MATISSE (MTTSSC)                      ',/,&
1304
 
'@    =========                                               ',/,&
1305
 
'@      DEBIT ENTHALPIQUE DES PANACHES INADAPTE.              ',/,&
1306
 
'@                                                            ',/,&
1307
 
'@    La modelisation des panaches est activee                ',/,&
1308
 
'@      avec IMDCNT = ',I10   ,'                              ',/,&
1309
 
'@    Le debit enthalpique des panaches de convection         ',/,&
1310
 
'@      naturelle doit etre strictement inferieur a la        ',/,&
1311
 
'@      puissance totale des colis presents dans l''entrepot  ',/,&
1312
 
'@      puisqu''il en represente une partie.                  ',/,&
1313
 
'@                                                            ',/,&
1314
 
'@    Le debit enthalpique prescrit est  DHPCNT  =',E12.5      ,/,&
1315
 
'@    La puissance totale reelle est PTOT*FRDTRA =',E12.5      ,/,&
1316
 
'@                                                            ',/,&
1317
 
'@    L''inegalite suivante n''est pas verifiee :             ',/,&
1318
 
'@          DHPCNT    <    PTOT*FRDTRA                        ',/,&
1319
 
'@       ',E12.5   ,'     ',E12.5   ,'                        ',/,&
1320
 
'@                                                            ',/,&
1321
 
'@    Le calcul ne sera pas execute.                          ',/,&
1322
 
'@                                                            ',/,&
1323
 
'@    Desactiver la modelisation des panaches ou              ',/,&
1324
 
'@      modifier la valeur du debit enthalpique prescrit ou   ',/,&
1325
 
'@      modifier la puissance totale (carte de remplissage ou ',/,&
1326
 
'@      puissance des colis).                                 ',/,&
1327
 
'@                                                            ',/,&
1328
 
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
1329
 
'@                                                            ',/)
1330
 
 
1331
 
!----
1332
 
! FIN
1333
 
!----
1334
 
 
1335
 
return
1336
 
 
1337
 
end subroutine