File: /nfs/home/0/users/jenkins/mfix.git/model/adjust_theta.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: ADJUST_THETA                                            C
4     !  Purpose: Remove small negative values of theta caused by linear     C
5     !           solvers                                                    C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 02-APR-98  C
8     !                                                                      C
9     !  Modified: S. Benyahia                              Date: 02-AUG-06  C
10     !  Purpose: check for small negative numbers at walls                  C
11     !           (not just fluid cells)                                     C
12     !                                                                      C
13     !                                                                      C
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
15     
16           SUBROUTINE ADJUST_THETA(M, IER)
17     
18     !-----------------------------------------------
19     ! Modules
20     !-----------------------------------------------
21           USE param1, only: zero
22           USE toleranc, only: zero_ep_s
23           USE constant, only: pi, to_si
24     ! granular temperature of solids phase m
25           USE fldvar, only: theta_m
26     ! material density of solids phase m
27           USE fldvar, only: ro_s
28     ! particle diameter of solids phase m
29           USE fldvar, only: d_p
30     ! number of solids phases
31           USE physprop, only: smax
32     ! kt types
33           USE run, only: kt_type
34           USE run, only: kt_type_enum
35           USE run, only: lun_1984
36           USE run, only: simonin_1996
37           USE run, only: ahmadi_1995
38           USE run, only: gd_1999
39           USE run, only: gtsh_2012
40           USE run, only: ia_2005
41           USE run, only: ghd_2007
42     ! needed for function.inc
43           USE compar
44           USE geometry
45           USE indices
46           USE functions
47     
48           IMPLICIT NONE
49     !-----------------------------------------------
50     ! Dummy arguments
51     !-----------------------------------------------
52     ! solids phase
53           INTEGER, INTENT(IN) :: M
54     ! error indicator
55           INTEGER, INTENT(INOUT) :: IER
56     !-----------------------------------------------
57     ! Local variables
58     !-----------------------------------------------
59     ! Indices
60           INTEGER :: IJK
61     ! Solids phase index
62           INTEGER :: L
63     ! Particle mass and diameter for use with those kinetic theories
64     ! that include mass of particle in definition of theta
65           DOUBLE PRECISION :: M_PM, D_PM
66     ! small value of theta_m
67           DOUBLE PRECISION :: smallTheta
68     !-----------------------------------------------
69     
70           IER = 0
71           smallTheta = (to_SI)**4 * ZERO_EP_S
72     
73           DO IJK = IJKSTART3, IJKEND3
74             IF ( FLUID_AT(IJK) ) THEN
75     
76               SELECT CASE(KT_TYPE_ENUM)
77                 CASE (LUN_1984, SIMONIN_1996, AHMADI_1995, GD_1999, &
78                       GTSH_2012)
79                   IF (THETA_M(IJK,M) < smallTheta) &
80                      THETA_M(IJK,M) = smallTheta
81     
82                 CASE (IA_2005)
83                   D_PM = D_P(IJK,M)
84                   M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
85                   IF (THETA_M(IJK,M) < smallTheta*M_PM) &
86                     THETA_M(IJK,M) = smallTheta*M_PM
87     
88                 CASE (GHD_2007)
89                   M_PM = ZERO
90                   DO L = 1,SMAX
91                     D_PM = D_P(IJK,L)
92                     M_PM = M_PM +(PI/6.d0)*(D_PM**3)*RO_S(IJK,L)
93                   ENDDO
94                   M_PM = M_PM/DBLE(SMAX)
95                   IF (THETA_M(IJK,M) < smallTheta*M_PM) &
96                     THETA_M(IJK,M) = smallTheta*M_PM
97     
98                 CASE DEFAULT
99     ! should never hit this
100                    WRITE (*, '(A)') 'ADJUST_THETA'
101                    WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
102                    call mfix_exit(myPE)
103               END SELECT   ! end selection of kt_type_enum
104             ENDIF   ! end if (fluid_at)
105           ENDDO  ! end do ijk
106     
107           RETURN
108           END SUBROUTINE ADJUST_THETA
109     
110