1
!-------------------------------------------------------------------------------
5
! This file is part of Code_Saturne, a general-purpose CFD tool.
7
! Copyright (C) 1998-2011 EDF S.A.
9
! This program is free software; you can redistribute it and/or modify it under
10
! the terms of the GNU General Public License as published by the Free Software
11
! Foundation; either version 2 of the License, or (at your option) any later
14
! This program is distributed in the hope that it will be useful, but WITHOUT
15
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19
! You should have received a copy of the GNU General Public License along with
20
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
21
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
23
!-------------------------------------------------------------------------------
28
( nvar , nscal , ncepdp , ncesmp , &
29
icepdc , icetsm , itypsm , &
30
dt , rtpa , propce , propfa , propfb , &
31
coefa , coefb , ckupdc , smacel , tinssa , divu , &
34
!===============================================================================
40
! Additional right-hand side source terms for the viscosity-like variable
41
! when using the Spalart-Allmaras model (ITURB=70)
46
! The additional source term is decomposed into an explicit part (crvexp) and
47
! an implicit part (crvimp) that must be provided here.
48
! The resulting equations solved by the code are:
50
! rho*volume*d(nusa)/dt + .... = crvimp*nusa + crvexp
52
! Note that crvexp, crvimp are defined after the Finite Volume
53
! integration over the cells, so they include the "volume" term. More precisely:
54
! - crvexp is expressed in kg.m2/s2
55
! - crvimp is expressed in kg/s
57
! The crvexp, crvimp arrays are already initialized to 0 before
58
! entering the routine. It is not needed to do it in the routine (waste of CPU time).
60
! For stability reasons, Code_Saturne will not add -crvimp directly to the
61
! diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp term is
62
! treated implicitely only if it strengthens the diagonal of the matrix.
63
! However, when using the second-order in time scheme, this limitation cannot
64
! be done anymore and -crvimp is added directly. The user should therefore test
65
! the negativity of crvimp by himself.
67
! When using the second-order in time scheme, one should supply:
69
! - crvimp at time n+1/2
71
! When entering the routine, two additional work arrays are already set for
72
! potential user need:
73
! tinssa = 2 (S11)**2 + 2 (S22)**2 + 2 (S33)**2
74
! + (2 S12)**2 + (2 S13)**2 + (2 S23)**2
76
! where Sij = (dUi/dxj+dUj/dxi)/2
78
! divu = du/dx + dv/dy + dw/dz
83
! The selection of cells where to apply the source terms is based on a getcel
84
! command. For more info on the syntax of the getcel command, refer to the
85
! user manual or to the comments on the similar command getfbr in the routine
88
!-------------------------------------------------------------------------------
90
!__________________.____._____.________________________________________________.
91
! name !type!mode ! role !
92
!__________________!____!_____!________________________________________________!
93
! nvar ! i ! <-- ! total number of variables !
94
! nscal ! i ! <-- ! total number of scalars !
95
! ncepdp ! i ! <-- ! number of cells with head loss terms !
96
! ncesmp ! i ! <-- ! number of cells with mass source term !
97
! icepdc(ncepdp) ! ia ! <-- ! index number of cells with head loss terms !
98
! icetsm(ncesmp) ! ia ! <-- ! index number of cells with mass source terms !
99
! itypsm ! ia ! <-- ! type of mass source term for each variable !
100
! (ncesmp,nvar) ! ! ! (see ustsma) !
101
! dt(ncelet) ! ra ! <-- ! time step (per cell) !
102
! rtpa ! ra ! <-- ! calculated variables at cell centers !
103
! (ncelet, *) ! ! ! (preceding time steps) !
104
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers !
105
! propfa(nfac, *) ! ra ! <-- ! physical properties at interior face centers !
106
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers !
107
! coefa, coefb ! ra ! <-- ! boundary conditions !
108
! (nfabor, *) ! ! ! !
109
! ckupdc(ncepdp,6) ! ra ! <-- ! head loss coefficient !
110
! smacel ! ra ! <-- ! value associated to each variable in the mass !
111
! (ncesmp,nvar) ! ! ! source terms or mass rate (see ustsma) !
112
! tinssa ! ra ! <-- ! tubulent production term (see comment above) !
113
! divu ! ra ! <-- ! velocity divergence (see comment above) !
114
! crvexp ! ra ! --> ! explicit part of the source term for k !
115
! crvimp ! ra ! --> ! implicit part of the source term for k !
116
!__________________!____!_____!________________________________________________!
118
! Type: i (integer), r (real), s (string), a (array), l (logical),
119
! and composite types (ex: ra real array)
120
! mode: <-- input, --> output, <-> modifies data, --- work array
121
!===============================================================================
123
!===============================================================================
125
!===============================================================================
136
!===============================================================================
143
integer ncepdp , ncesmp
145
integer icepdc(ncepdp)
146
integer icetsm(ncesmp), itypsm(ncesmp,nvar)
148
double precision dt(ncelet), rtpa(ncelet,*)
149
double precision propce(ncelet,*)
150
double precision propfa(nfac,*), propfb(nfabor,*)
151
double precision coefa(nfabor,*), coefb(nfabor,*)
152
double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar)
153
double precision tinssa(ncelet), divu(ncelet)
154
double precision crvexp(ncelet), crvimp(ncelet)
159
double precision ff, tau, xx
161
integer, allocatable, dimension(:) :: lstelt
163
!===============================================================================
165
! TEST_TO_REMOVE_FOR_USE_OF_SUBROUTINE_START
166
!===============================================================================
170
!===============================================================================
171
! TEST_TO_REMOVE_FOR_USE_OF_SUBROUTINE_END
174
!===============================================================================
176
!===============================================================================
178
! Allocate a temporary array for cells selection
179
allocate(lstelt(ncel))
182
! --- Index number of the density in the propce array
183
ipcrom = ipproc(irom)
185
if(iwarni(inusa).ge.1) then
189
!===============================================================================
190
! 2. Example of arbitrary additional source term for nusa
192
! Source term for k :
193
! rho volume d(k)/dt = ...
194
! ... - rho*volume*ff - rho*volume*nusa/tau
196
! With xx = 2.d0, ff=3.d0 and tau = 4.d0
198
!===============================================================================
200
! --- Explicit source terms
207
crvexp(iel) = -propce(iel,ipcrom)*volume(iel)*ff
210
! --- Implicit source terms
211
! crvimp is already initialized to 0, no need to set it here
214
crvimp(iel) = -propce(iel,ipcrom)*volume(iel)/tau
221
1000 format(' User source terms for Spalart-Allmaras',/)
227
! Deallocate the temporary array
added
removed
users/base/ustssa.f90
