1
!-------------------------------------------------------------------------------
3
! This file is part of Code_Saturne, a general-purpose CFD tool.
5
! Copyright (C) 1998-2011 EDF S.A.
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
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
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.
21
!-------------------------------------------------------------------------------
26
( nvar , nscal , lndnod , icocel , itycel ,ifrlag , &
27
nbpmax , nvp , nvp1 , nvep , nivep , &
28
ntersl , nvlsta , nvisbr , &
31
dt , rtpa , propce , propfa , propfb , &
32
ettp , ettpa , tepa , statis , &
33
taup , tlag , piil , &
34
bx , vagaus , gradpr , gradvf , romp , &
35
brgaus , terbru , fextla )
37
!===============================================================================
41
! Subroutine of the Lagrangian particle-tracking module:
42
! ------------------------------------------------------
45
! Deposition sub-model:
47
! Main subroutine of the submodel
49
! 1/ Calculation of the normalized wall-normal distance of the boundary-cell particles
50
! 2/ Sorting of the particles with respect to their normalized wall-normal distance
51
! * If y^+ > depint : the standard Langevin model is applied
52
! * If y^+ < depint : the deposition submodel is applied
54
!-------------------------------------------------------------------------------
56
!__________________.____._____.________________________________________________.
57
! name !type!mode ! role !
58
!__________________!____!_____!________________________________________________!
59
! lndnod ! e ! <-- ! dim. connectivite cellules->faces !
60
! nvar ! i ! <-- ! total number of variables !
61
! nscal ! i ! <-- ! total number of scalars !
62
! nbpmax ! e ! <-- ! nombre max de particulies autorise !
63
! nvp ! e ! <-- ! nombre de variables particulaires !
64
! nvp1 ! e ! <-- ! nvp sans position, vfluide, vpart !
65
! nvep ! e ! <-- ! nombre info particulaires (reels) !
66
! nivep ! e ! <-- ! nombre info particulaires (entiers) !
67
! ntersl ! e ! <-- ! nbr termes sources de couplage retour !
68
! nvlsta ! e ! <-- ! nombre de var statistiques lagrangien !
69
! nvisbr ! e ! <-- ! nombre de statistiques aux frontieres !
70
! icocel ! te ! --> ! connectivite cellules -> faces !
71
! (lndnod) ! ! ! face de bord si numero negatif !
72
! itycel ! te ! --> ! connectivite cellules -> faces !
73
! (ncelet+1) ! ! ! pointeur du tableau icocel !
74
! ifrlag ! te ! --> ! numero de zone de la face de bord !
75
! (nfabor) ! ! ! pour le module lagrangien !
76
! itepa ! te ! <-- ! info particulaires (entiers) !
77
! (nbpmax,nivep ! ! ! (cellule de la particule,...) !
78
! dlgeo ! tr ! --> ! tableau contenant les donnees geometriques !
79
!(nfabor,ngeol) ! ! ! !
80
! dt(ncelet) ! ra ! <-- ! time step (per cell) !
81
! rtpa ! tr ! <-- ! variables de calcul au centre des !
82
! (ncelet,*) ! ! ! cellules (pas de temps precedent) !
83
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers !
84
! propfa(nfac, *) ! ra ! <-- ! physical properties at interior face centers !
85
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers !
86
! ettp ! tr ! --> ! tableaux des variables liees !
87
! (nbpmax,nvp) ! ! ! aux particules etape courante !
88
! ettpa ! tr ! <-- ! tableaux des variables liees !
89
! (nbpmax,nvp) ! ! ! aux particules etape precedente !
90
! tepa ! tr ! <-- ! info particulaires (reels) !
91
! (nbpmax,nvep) ! ! ! (poids statistiques,...) !
92
! statis ! tr ! <-- ! cumul des statistiques volumiques !
93
!(ncelet,nvlsta ! ! ! !
94
! taup(nbpmax) ! tr ! <-- ! temps caracteristique dynamique !
95
! tlag(nbpmax) ! tr ! <-- ! temps caracteristique fluide !
96
! piil(nbpmax,3 ! tr ! <-- ! terme dans l'integration des eds up !
97
! bx(nbpmax,3,2 ! tr ! <-- ! caracteristiques de la turbulence !
98
! vagaus ! tr ! <-- ! variables aleatoires gaussiennes !
99
!(nbpmax,nvgaus ! ! ! !
100
! gradpr(ncel,3 ! tr ! <-- ! gradient de pression !
101
! gradvf(ncel,3 ! tr ! <-- ! gradient de la vitesse du fluide !
102
! romp ! tr ! <-- ! masse volumique des particules !
103
! fextla ! tr ! <-- ! champ de forces exterieur !
104
!(ncelet,3) ! ! ! utilisateur (m/s2) !
105
!__________________!____!_____!________________________________________________!
107
! TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
108
! L (LOGIQUE) .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
109
! MODE : <-- donnee, --> resultat, <-> Donnee modifiee
110
! --- tableau de travail
111
!===============================================================================
113
!===============================================================================
115
!===============================================================================
131
!===============================================================================
138
integer nbpmax , nvp , nvp1 , nvep , nivep
139
integer ntersl , nvlsta , nvisbr , lndnod
141
integer itepa(nbpmax,nivep)
142
integer icocel(lndnod), ifrlag(nfabor), itycel(ncelet+1)
144
double precision dt(ncelet) , rtpa(ncelet,*)
145
double precision propce(ncelet,*)
146
double precision propfa(nfac,*) , propfb(nfabor,*)
147
double precision ettp(nbpmax,nvp) , ettpa(nbpmax,nvp)
148
double precision tepa(nbpmax,nvep) , statis(ncelet,*)
149
double precision taup(nbpmax) , tlag(nbpmax,3)
150
double precision piil(nbpmax,3) , bx(nbpmax,3,2)
151
double precision vagaus(nbpmax,*)
152
double precision brgaus(nbpmax,*) , terbru(nbpmax)
153
double precision gradpr(ncelet,3) , gradvf(ncelet,9)
154
double precision romp(nbpmax)
155
double precision fextla(nbpmax,3)
156
double precision dlgeo(nfabor,ngeol)
162
integer iel , ip , id , i0 , iromf , mode
165
double precision aa , bb , cc , dd , ee
166
double precision aux1 , aux2 ,aux3 , aux4 , aux5 , aux6
167
double precision aux7 , aux8 , aux9 , aux10 , aux11
168
double precision ter1f , ter2f , ter3f
169
double precision ter1p , ter2p , ter3p , ter4p , ter5p
170
double precision ter1x , ter2x , ter3x , ter4x , ter5x
171
double precision tci , force
172
double precision gama2 , omegam , omega2
173
double precision grga2 , gagam , gaome
174
double precision p11 , p21 , p22 , p31 , p32 , p33
175
double precision grav(3) , romf , vitf
176
double precision tempf, lvisq, tvisq
177
double precision d3 , vpart, vvue
178
double precision px , py , pz , distp , d1
179
double precision dismin,dismax, ustar, visccf,depint
181
integer, allocatable, dimension(:) :: ifacl
183
!===============================================================================
185
!===============================================================================
186
! 0. Memory management
187
!===============================================================================
190
! Allocae a work array
191
allocate(ifacl(nbpart))
193
!===============================================================================
195
!===============================================================================
197
! Initialize variables to avoid compiler warnings
205
! Interface location between near-wall region
206
! and core of the flow (normalized units)
210
! The density pointer according to the flow location
212
if ( ippmod(icp3pl).ge.0 .or. ippmod(icfuel).ge.0 ) then
213
iromf = ipproc(irom1)
219
!===============================================================================
220
! 2. loop on the particles
221
!===============================================================================
225
if (itepa(ip,jisor).gt.0) then
227
iel = itepa(ip,jisor)
228
romf = propce(iel,iromf)
229
visccf = propce(iel,ipproc(iviscl)) / romf
231
! Fluid temperature computation depending on the type of flow
233
if (ippmod(icp3pl).ge.0 .or. &
234
ippmod(icpl3c).ge.0 .or. &
235
ippmod(icfuel).ge.0) then
237
tempf = propce(iel,ipproc(itemp1))
239
else if (ippmod(icod3p).ge.0 .or. &
240
ippmod(icoebu).ge.0 .or. &
241
ippmod(ielarc).ge.0 .or. &
242
ippmod(ieljou).ge.0 ) then
244
tempf = propce(iel,ipproc(itemp))
246
else if (iscalt.gt.0) then
247
if (iscsth(iscalt).eq.-1) then
248
tempf = rtpa(iel,isca(iscalt)) + tkelvi
250
else if ( iscsth(iscalt).eq.1 ) then
251
tempf = rtpa(iel,isca(iscalt))
253
else if ( iscsth(iscalt).eq.2 ) then
255
call usthht(mode,rtpa(iel,isca(iscalt)),tempf)
263
!===============================================================================
264
! 2. calculation of the normalized wall-normal particle distance (y^+)
265
!===============================================================================
272
tepa(ip,jryplu) = 1.0d3
275
do il = itycel(iel),itycel(iel+1)-1
285
! Test if the particle is located in a boundary cell
288
if ( iusclb(izone) .eq. idepo1 .or. &
289
iusclb(izone) .eq. idepo2 .or. &
290
iusclb(izone) .eq. idepo3 .or. &
291
iusclb(izone) .eq. idepfa .or. &
292
iusclb(izone) .eq. irebol ) then
295
! Calculation of the wall units
298
lvisq = visccf / ustar
299
tvisq = visccf / (ustar * ustar)
305
d1 = abs(px*dlgeo(ifac,1)+py*dlgeo(ifac,2) &
306
+pz*dlgeo(ifac,3)+dlgeo(ifac,4)) &
307
/sqrt( dlgeo(ifac,1)*dlgeo(ifac,1) &
308
+dlgeo(ifac,2)*dlgeo(ifac,2) &
309
+dlgeo(ifac,3)*dlgeo(ifac,3))
311
if (d1.lt.distp) then
314
tepa(ip,jryplu)= distp/lvisq
323
!=========================================================================
324
! If y^+ is greater than the interface location,
325
! the standard model is applied
326
!=========================================================================
328
if (tepa(ip,jryplu).gt.depint) then
330
itepa(ip,jimark) = -1
335
if (id.eq.1) vitf = rtpa(iel,iu)
336
if (id.eq.2) vitf = rtpa(iel,iv)
337
if (id.eq.3) vitf = rtpa(iel,iw)
339
tci = piil(ip,id) * tlag(ip,id) + vitf
341
force = (romf * gradpr(iel,id) / romp(ip) &
342
+ grav(id) + fextla(ip,id)) * taup(ip)
344
aux1 = exp( -dtp / taup(ip))
345
aux2 = exp( -dtp / tlag(ip,id))
346
aux3 = tlag(ip,id) / (tlag(ip,id)-taup(ip))
348
aux4 = tlag(ip,id) / (tlag(ip,id)+taup(ip))
349
aux5 = tlag(ip,id) * (1.d0-aux2)
350
aux6 = bx(ip,id,nor) * bx(ip,id,nor) * tlag(ip,id)
352
aux7 = tlag(ip,id) - taup(ip)
353
aux8 = bx(ip,id,nor) * bx(ip,id,nor) * aux3**2
355
!---> trajectory terms
357
aa = taup(ip) * (1.d0 - aux1)
358
bb = (aux5 - aa) * aux3
361
ter1x = aa * ettpa(ip,jup+i0)
362
ter2x = bb * ettpa(ip,juf+i0)
364
ter4x = (dtp - aa) * force
368
ter1f = ettpa(ip,juf+i0) * aux2
369
ter2f = tci * (1.d0-aux2)
371
!---> termes pour la vitesse des particules
373
dd = aux3 * (aux2 - aux1)
376
ter1p = ettpa(ip,jup+i0) * aux1
377
ter2p = ettpa(ip,juf+i0) * dd
378
ter3p = tci * (ee-dd)
381
!---> (2.3) Coefficients computation for the stochastic integral
383
!---> Integral on the particles position
385
gama2 = 0.5d0 * (1.d0 - aux2*aux2 )
386
omegam = 0.5d0 * aux4 * ( aux5 - aux2*aa ) &
388
omegam = omegam * sqrt(aux6)
391
* (aux7*dtp - 2.d0 * (tlag(ip,id)*aux5-taup(ip)*aa)) &
392
+ 0.5d0 * tlag(ip,id) * tlag(ip,id) * aux5 &
394
+ 0.5d0 * taup(ip) * taup(ip) * aa * (1.d0+aux1) &
395
- 2.d0 * aux4 * tlag(ip,id) * taup(ip) * taup(ip) &
398
omega2 = aux8 * omega2
401
if (abs(gama2).gt.epzero) then
403
p21 = omegam / sqrt(gama2)
404
p22 = omega2 - p21**2
406
p22 = sqrt( max(zero,p22) )
413
ter5x = p21 * vagaus(ip,id) + p22 * vagaus(ip,3+id)
415
!---> integral on the vu fluid velocity
417
p11 = sqrt( gama2*aux6 )
418
ter3f = p11*vagaus(ip,id)
420
!---> integrale on the particules velocity
422
aux9 = 0.5d0 * tlag(ip,id) * (1.d0 - aux2*aux2)
423
aux10 = 0.5d0 * taup(ip) * (1.d0 - aux1*aux1)
424
aux11 = taup(ip) * tlag(ip,id) * (1.d0 - aux1*aux2) &
425
/ (taup(ip) + tlag(ip,id))
427
grga2 = (aux9 - 2.d0*aux11 + aux10) * aux8
428
gagam = (aux9 - aux11) * (aux8 / aux3)
429
gaome = ( (tlag(ip,id) - taup(ip)) * (aux5 - aa) &
430
- tlag(ip,id) * aux9 - taup(ip) * aux10 &
431
+ (tlag(ip,id) + taup(ip)) * aux11) * aux8
433
if(p11.gt.epzero) then
439
if(p22.gt.epzero) then
440
p32 = (gaome-p31*p21) / p22
445
p33 = grga2 - p31**2 - p32**2
447
p33 = sqrt( max(zero,p33) )
449
ter5p = p31 * vagaus(ip,id) &
450
+ p32 * vagaus(ip,3+id) &
451
+ p33 * vagaus(ip,6+id)
454
! Update of the particle state vector
457
ettp(ip,jxp+i0) = ettpa(ip,jxp+i0) &
458
+ ter1x + ter2x + ter3x + ter4x + ter5x
460
ettp(ip,juf+i0) = ter1f + ter2f + ter3f
462
ettp(ip,jup+i0) = ter1p + ter2p + ter3p + ter4p + ter5p
466
!===============================================================================
467
! Otherwise, the deposition submodel is applied :
468
!!===============================================================================
472
if (tepa(ip,jryplu).lt.tepa(ip,jrinpf)) then
474
if ( itepa(ip,jimark) .lt. 0 ) then
475
itepa(ip,jimark) = 10
480
if ( itepa(ip,jimark) .eq. 0 .and. &
481
itepa(ip,jdiel) .eq. iel .and. &
482
ifacl(ip) .ne. itepa(ip,jdfac)) then
483
itepa(ip,jimark) = 10
488
! if (jrinpf < y^+ < depint)
490
if ( itepa(ip,jimark) .lt. 0 ) then
491
itepa(ip,jimark) = 20
492
else if ( itepa(ip,jimark) .eq. 0 ) then
493
itepa(ip,jimark) = 30
498
itepa(ip,jdfac) = ifacl(ip)
499
itepa(ip,jdiel) = itepa(ip,jisor)
501
dismin=min(dismin,distp)
502
dismax=max(dismax,distp)
508
nbpmax , nvp , nvp1 , nvep , nivep , &
509
ntersl , nvlsta , nvisbr , &
512
dt , rtpa , propce , propfa , propfb , &
513
ettp , ettpa , tepa , &
514
statis , taup , tlag , piil , &
515
bx , vagaus , gradpr , gradvf , romp , &
516
tempf , romf , ustar , lvisq ,tvisq , depint )
added
removed
src/lagr/lagdep.f90
