← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/base/vorvit.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/base/vorvit.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 vorvit &
29
 
!================
30
 
 
31
 
 ( ncevor , nvor   , ient   , ivorce , visco  ,                   &
32
 
   yzcel  , xu     , xv     , xw     ,                            &
33
 
   yzvor  , signv  , sigma  , gamma  , temps )
34
 
 
35
 
!===============================================================================
36
 
!  FONCTION  :
37
 
!  ---------
38
 
 
39
 
! METHODE DES VORTEX POUR LES CONDITIONS AUX LIMITES D'ENTREE EN
40
 
!   L.E.S. : CALCUL DES LA VITESSE DES VORTEX
41
 
 
42
 
!-------------------------------------------------------------------------------
43
 
! Arguments
44
 
!__________________.____._____.________________________________________________.
45
 
! name             !type!mode ! role                                           !
46
 
!__________________!____!_____!________________________________________________!
47
 
! ncevor           ! e  ! <-- ! nombre de face a l'entree ou est               !
48
 
!                  !    !     ! utilise la methode                             !
49
 
! nvor             ! e  !     ! nombre de vortex a l'entree ou est             !
50
 
!                  !    !     ! utilisee la methode                            !
51
 
! ient             ! e  ! <-- ! numero de l'entree ou est utilisee             !
52
 
!                  !    !     ! la methode                                     !
53
 
! ivorce           ! te ! <-- ! numero du vortex le plus proche d'une          !
54
 
!     (nvomax)     !    !     ! face donnee                                    !
55
 
! visco            ! tr ! <-- ! viscosite cinematique sur les faces            !
56
 
!(icvmax,nnent)    !    !     ! d'entree                                       !
57
 
! yzcel            ! tr ! <-- ! coordonnees des faces d'entree dans            !
58
 
!   (icvmax ,2)    !    !     ! le referentiel local                           !
59
 
! xu(icvmax)       ! tr ! <-- ! composante de vitesse principale               !
60
 
! xv(icvmax)       ! tr ! <-- ! composantes de vitesse transverses             !
61
 
! xw(icvmax)       ! tr ! <-- !                                                !
62
 
! yzvor            ! tr ! <-- ! coordonnees du centre des vortex               !
63
 
!   (nvomax,2)     !    !     !                                                !
64
 
! signv(nvomax)    ! tr ! <-- ! sens de rotation des vortex                    !
65
 
! sigma            ! tr ! <-- ! taille des vortex                              !
66
 
!(nvomax,nnent)    !    !     !                                                !
67
 
! gamma            ! tr ! <-- ! intensite (circulation) des vortex             !
68
 
!(nvomax,2,nnen    !    !     ! dans les deux directions du plan               !
69
 
! temps            ! tr ! <-- ! temps ecoule depuis la creation                !
70
 
!     (nvomax)     !    !     ! du vortex                                      !
71
 
!__________________!____!_____!________________________________________________!
72
 
 
73
 
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
74
 
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
75
 
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
76
 
!            --- tableau de travail
77
 
!===============================================================================
78
 
 
79
 
implicit none
80
 
 
81
 
!===============================================================================
82
 
! Common blocks
83
 
!===============================================================================
84
 
 
85
 
include "paramx.h"
86
 
include "cstnum.h"
87
 
include "cstphy.h"
88
 
include "entsor.h"
89
 
include "vortex.h"
90
 
 
91
 
!===============================================================================
92
 
 
93
 
! Arguments
94
 
 
95
 
integer          ncevor , nvor   , ient
96
 
integer          ivorce(nvomax)
97
 
 
98
 
double precision yzcel(icvmax,2) , visco(icvmax)
99
 
double precision xu(icvmax )     , xv(icvmax )     , xw(icvmax )
100
 
double precision yzvor(nvomax,2) , signv(nvomax)   , sigma(nvomax)
101
 
double precision gamma(nvomax,2) , temps(nvomax)
102
 
 
103
 
! Local variables
104
 
 
105
 
 
106
 
integer          ii, jj, iii
107
 
 
108
 
double precision yy, zz, xvisc
109
 
double precision norme, vv, ww, theta, dv, dw
110
 
double precision yparoi, zparoi, yperio, zperio, ysym, zsym
111
 
double precision phidat
112
 
double precision alpha, ek_vor, ee_vor, v_vor, w_vor
113
 
double precision lt, lk, deuxpi
114
 
 
115
 
!===============================================================================
116
 
! 1. CALCUL DE GAMMA
117
 
!===============================================================================
118
 
 
119
 
alpha = 4.d0 * sqrt(pi * xsurfv(ient)/                            &
120
 
       (3.d0 * nvor*(2.d0*log(3.d0)-3.d0*log(2.d0))))
121
 
 
122
 
do ii = 1, nvor
123
 
  yy = yzvor(ii,1)
124
 
  zz = yzvor(ii,2)
125
 
  iii = 0
126
 
  ek_vor =  phidat(nfecra,icas(ient),ndat(ient),yy,zz,            &
127
 
       ydat(1,ient),zdat(1,ient),kdat(1,ient),iii)
128
 
  gamma(ii,1) = alpha*sqrt(ek_vor)*signv(ii)
129
 
  gamma(ii,2) = gamma(ii,1)
130
 
enddo
131
 
 
132
 
!===============================================================================
133
 
! 2. CALCUL DE SIGMA
134
 
!===============================================================================
135
 
 
136
 
if(isgmvo(ient).eq.1) then
137
 
  do ii = 1, nvor
138
 
    sigma(ii) = xsgmvo(ient)
139
 
  enddo
140
 
elseif (isgmvo(ient).eq.2) then
141
 
  do ii = 1, nvor
142
 
    yy = yzvor(ii,1)
143
 
    zz = yzvor(ii,2)
144
 
    iii = 0
145
 
    ek_vor =  phidat(nfecra,icas(ient),ndat(ient),yy,zz,          &
146
 
         ydat(1,ient),zdat(1,ient),kdat(1,ient),iii)
147
 
    ee_vor =  phidat(nfecra,icas(ient),ndat(ient),yy,zz,          &
148
 
         ydat(1,ient),zdat(1,ient),epsdat(1,ient),iii)
149
 
    sigma(ii) = (cmu**(0.75d0))*(ek_vor**(1.5d0))/ee_vor
150
 
  enddo
151
 
elseif (isgmvo(ient).eq.3) then
152
 
  do ii = 1, nvor
153
 
    yy = yzvor(ii,1)
154
 
    zz = yzvor(ii,2)
155
 
    iii = 0
156
 
    ek_vor =  phidat(nfecra,icas(ient),ndat(ient),yy,zz,          &
157
 
         ydat(1,ient),zdat(1,ient),kdat(1,ient),iii)
158
 
    ee_vor =  phidat(nfecra,icas(ient),ndat(ient),yy,zz,          &
159
 
         ydat(1,ient),zdat(1,ient),epsdat(1,ient),iii)
160
 
    xvisc  =  visco(ivorce(ii))
161
 
    lt = sqrt(5.d0*xvisc*ek_vor/ee_vor)
162
 
    lk = 200.d0*(xvisc**3/ee_vor)**(0.25d0)
163
 
    sigma(ii) = max(lt,lk)
164
 
  enddo
165
 
endif
166
 
 
167
 
!===============================================================================
168
 
! 3. CALCUL DU CHAMP DE VITESSE INDUIT (AU CENTRE DES CELLULES)
169
 
!===============================================================================
170
 
 
171
 
deuxpi = 2.d0*pi
172
 
 
173
 
do ii = 1, ncevor
174
 
  vv = 0.d0
175
 
  ww = 0.d0
176
 
  do jj = 1, nvor
177
 
    yy = yzvor(jj,1)-yzcel(ii,1)
178
 
    zz = yzvor(jj,2)-yzcel(ii,2)
179
 
    norme = yy**2+zz**2
180
 
    alpha = sigma(jj)
181
 
    if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
182
 
      theta = -norme/(2.d0*alpha**2)
183
 
      theta = exp(theta)
184
 
      dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
185
 
      dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
186
 
      vv = vv - dv
187
 
      ww = ww + dw
188
 
    endif
189
 
 
190
 
!===============================================================================
191
 
! 4. TRAITEMENT DES CONDITIONS AUX LIMITES
192
 
!===============================================================================
193
 
! suite des boucles en DO
194
 
! -----------------------
195
 
! Conduite rectangulaire
196
 
! -----------------------
197
 
    if(icas(ient).eq.1) then
198
 
 
199
 
! Periodicites
200
 
 
201
 
      if(iclvor(1,ient).eq.3) then
202
 
        if(yzvor(jj,1).gt.0.d0) then
203
 
          yperio = yzvor(jj,1) - lly(ient)
204
 
          zperio = yzvor(jj,2)
205
 
        else
206
 
          yperio = yzvor(jj,1) + lly(ient)
207
 
          zperio = yzvor(jj,2)
208
 
        endif
209
 
        yy = yperio - yzcel(ii,1)
210
 
        zz = zperio - yzcel(ii,2)
211
 
        norme = yy**2+zz**2
212
 
        alpha = sigma(jj)
213
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
214
 
          theta = -norme/(2.d0*alpha**2)
215
 
          theta = exp(theta)
216
 
          dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
217
 
          dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
218
 
          vv = vv - dv
219
 
          ww = ww + dw
220
 
        endif
221
 
      endif
222
 
 
223
 
      if(iclvor(2,ient).eq.3) then
224
 
        if(yzvor(jj,2).gt.0.d0) then
225
 
          yperio = yzvor(jj,1)
226
 
          zperio = yzvor(jj,2) - llz(ient)
227
 
        else
228
 
          yperio = yzvor(jj,1)
229
 
          zperio = yzvor(jj,2) + llz(ient)
230
 
        endif
231
 
        yy = yperio - yzcel(ii,1)
232
 
        zz = zperio - yzcel(ii,2)
233
 
        norme = yy**2+zz**2
234
 
        alpha = sigma(jj)
235
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
236
 
          theta = -norme/(2.d0*alpha**2)
237
 
          theta = exp(theta)
238
 
          dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
239
 
          dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
240
 
          vv = vv - dv
241
 
          ww = ww + dw
242
 
        endif
243
 
      endif
244
 
 
245
 
! Parois
246
 
 
247
 
      if(iclvor(1,ient).eq.1) then
248
 
        yparoi = lly(ient)/2.d0
249
 
        zparoi = yzcel(ii,2)
250
 
        yparoi = 2.d0*yparoi - yzvor(jj,1)
251
 
        zparoi = 2.d0*zparoi - yzvor(jj,2)
252
 
        yy = yparoi - yzcel(ii,1)
253
 
        zz = zparoi - yzcel(ii,2)
254
 
        norme = yy**2+zz**2
255
 
        alpha = sigma(jj)
256
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
257
 
          theta = -norme/(2.d0*alpha**2)
258
 
          theta = exp(theta)
259
 
          dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
260
 
          dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
261
 
          vv = vv - dv
262
 
          ww = ww + dw
263
 
        endif
264
 
      endif
265
 
 
266
 
      if(iclvor(3,ient).eq.1) then
267
 
        yparoi = - lly(ient)/2.d0
268
 
        zparoi = yzcel(ii,2)
269
 
        yparoi = 2.d0*yparoi - yzvor(jj,1)
270
 
        zparoi = 2.d0*zparoi - yzvor(jj,2)
271
 
        yy = yparoi - yzcel(ii,1)
272
 
        zz = zparoi - yzcel(ii,2)
273
 
        norme = yy**2+zz**2
274
 
        alpha = sigma(jj)
275
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
276
 
          theta = -norme/(2.d0*alpha**2)
277
 
          theta = exp(theta)
278
 
          dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
279
 
          dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
280
 
          vv = vv - dv
281
 
          ww = ww + dw
282
 
        endif
283
 
      endif
284
 
 
285
 
      if(iclvor(2,ient).eq.1) then
286
 
        yparoi = yzcel(ii,1)
287
 
        zparoi = llz(ient)/2.d0
288
 
        yparoi = 2.d0*yparoi - yzvor(jj,1)
289
 
        zparoi = 2.d0*zparoi - yzvor(jj,2)
290
 
        yy = yparoi - yzcel(ii,1)
291
 
        zz = zparoi - yzcel(ii,2)
292
 
        norme = yy**2+zz**2
293
 
        alpha = sigma(jj)
294
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
295
 
          theta = -norme/(2.d0*alpha**2)
296
 
          theta = exp(theta)
297
 
          dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
298
 
          dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
299
 
          vv = vv - dv
300
 
          ww = ww + dw
301
 
        endif
302
 
      endif
303
 
 
304
 
      if(iclvor(4,ient).eq.1) then
305
 
        yparoi = yzcel(ii,1)
306
 
        zparoi = -llz(ient)/2.d0
307
 
        yparoi = 2.d0*yparoi - yzvor(jj,1)
308
 
        zparoi = 2.d0*zparoi - yzvor(jj,2)
309
 
        yy = yparoi - yzcel(ii,1)
310
 
        zz = zparoi - yzcel(ii,2)
311
 
        norme = yy**2+zz**2
312
 
        alpha = sigma(jj)
313
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
314
 
          theta = -norme/(2.d0*alpha**2)
315
 
          theta = exp(theta)
316
 
          dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
317
 
          dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
318
 
          vv = vv - dv
319
 
          ww = ww + dw
320
 
        endif
321
 
      endif
322
 
 
323
 
! Symetries
324
 
 
325
 
      if(iclvor(1,ient).eq.2) then
326
 
        ysym = lly(ient)/2.d0
327
 
        zsym = yzcel(ii,2)
328
 
        ysym = 2.d0*ysym - yzvor(jj,1)
329
 
        zsym = 2.d0*zsym - yzvor(jj,2)
330
 
        yy = ysym - yzcel(ii,1)
331
 
        zz = zsym - yzcel(ii,2)
332
 
        norme = yy**2+zz**2
333
 
        alpha = sigma(jj)
334
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
335
 
          theta = -norme/(2.d0*alpha**2)
336
 
          theta = exp(theta)
337
 
          dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
338
 
          vv = vv - dv
339
 
        endif
340
 
      endif
341
 
 
342
 
      if(iclvor(3,ient).eq.2) then
343
 
        ysym = - lly(ient)/2.d0
344
 
        zsym = yzcel(ii,2)
345
 
        ysym = 2.d0*ysym - yzvor(jj,1)
346
 
        zsym = 2.d0*zsym - yzvor(jj,2)
347
 
        yy = ysym - yzcel(ii,1)
348
 
        zz = zsym - yzcel(ii,2)
349
 
        norme = yy**2+zz**2
350
 
        alpha = sigma(jj)
351
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
352
 
          theta = -norme/(2.d0*alpha**2)
353
 
          theta = exp(theta)
354
 
          dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
355
 
          vv = vv - dv
356
 
        endif
357
 
      endif
358
 
 
359
 
      if(iclvor(2,ient).eq.2) then
360
 
        ysym = yzcel(ii,1)
361
 
        zsym = llz(ient)/2.d0
362
 
        ysym = 2.d0*ysym - yzvor(jj,1)
363
 
        zsym = 2.d0*zsym - yzvor(jj,2)
364
 
        yy = ysym - yzcel(ii,1)
365
 
        zz = zsym - yzcel(ii,2)
366
 
        norme = yy**2+zz**2
367
 
        alpha = sigma(jj)
368
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
369
 
          theta = -norme/(2.d0*alpha**2)
370
 
          theta = exp(theta)
371
 
          dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
372
 
          ww = ww + dw
373
 
        endif
374
 
      endif
375
 
 
376
 
      if(iclvor(4,ient).eq.2) then
377
 
        ysym = yzcel(ii,1)
378
 
        zsym = -llz(ient)/2.d0
379
 
        ysym = 2.d0*ysym - yzvor(jj,1)
380
 
        zsym = 2.d0*zsym - yzvor(jj,2)
381
 
        yy = ysym - yzcel(ii,1)
382
 
        zz = zsym - yzcel(ii,2)
383
 
        norme = yy**2+zz**2
384
 
        alpha = sigma(jj)
385
 
        if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
386
 
          theta = -norme/(2.d0*alpha**2)
387
 
          theta = exp(theta)
388
 
          dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
389
 
          ww = ww + dw
390
 
        endif
391
 
      endif
392
 
 
393
 
! --------------------
394
 
! Conduite circulaire
395
 
! --------------------
396
 
    elseif(icas(ient).eq.2) then
397
 
 
398
 
      yparoi = yzcel(ii,1)*((lld(ient)/2.0d0)                     &
399
 
           /sqrt(yzcel(ii,1)**2 + yzcel(ii,2)**2))
400
 
      zparoi = yzcel(ii,2)*((lld(ient)/2.0d0)                     &
401
 
           /sqrt(yzcel(ii,1)**2 + yzcel(ii,2)**2))
402
 
      yparoi = 2.d0*yparoi - yzvor(jj,1)
403
 
      zparoi = 2.d0*zparoi - yzvor(jj,2)
404
 
 
405
 
      yy = yparoi - yzcel(ii,1)
406
 
      zz = zparoi - yzcel(ii,2)
407
 
      norme = yy**2+zz**2
408
 
      alpha = sigma(jj)
409
 
 
410
 
      if(norme.gt.epzero.and.norme.lt.4.d0*alpha) then
411
 
        theta = -norme/(2.d0*alpha**2)
412
 
        theta = exp(theta)
413
 
        dv = zz/norme*(1.d0-theta)*gamma(jj,1)*theta/deuxpi
414
 
        dw = yy/norme*(1.d0-theta)*gamma(jj,2)*theta/deuxpi
415
 
        vv = vv - dv
416
 
        ww = ww + dw
417
 
      endif
418
 
 
419
 
    endif
420
 
 
421
 
  enddo
422
 
  xv(ii) = vv
423
 
  xw(ii) = ww
424
 
enddo
425
 
 
426
 
!===============================================================================
427
 
! 5. AJOUT DE LA VITESSE MOYENNE POUR V ET W
428
 
!===============================================================================
429
 
if(icas(ient).eq.1.or.icas(ient).eq.2.or.icas(ient).eq.3) then
430
 
  do ii = 1, ncevor
431
 
    yy = yzcel(ii,1)
432
 
    zz = yzcel(ii,2)
433
 
    iii = 0
434
 
    v_vor = phidat(nfecra,icas(ient),ndat(ient),yy,zz,            &
435
 
         ydat(1,ient),zdat(1,ient),vdat(1,ient),iii)
436
 
    w_vor = phidat(nfecra,icas(ient),ndat(ient),yy,zz,            &
437
 
         ydat(1,ient),zdat(1,ient),wdat(1,ient),iii)
438
 
    xv(ii) = xv(ii) + v_vor
439
 
    xw(ii) = xw(ii) + w_vor
440
 
  enddo
441
 
endif
442
 
 
443
 
!----
444
 
! FIN
445
 
!----
446
 
 
447
 
return
448
 
end subroutine