File: N:\mfix\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 exit, only: mfix_exit
45           USE functions
46           USE geometry
47           USE indices
48     
49           IMPLICIT NONE
50     !-----------------------------------------------
51     ! Dummy arguments
52     !-----------------------------------------------
53     ! solids phase
54           INTEGER, INTENT(IN) :: M
55     ! error indicator
56           INTEGER, INTENT(INOUT) :: IER
57     !-----------------------------------------------
58     ! Local variables
59     !-----------------------------------------------
60     ! Indices
61           INTEGER :: IJK
62     ! Solids phase index
63           INTEGER :: L
64     ! Particle mass and diameter for use with those kinetic theories
65     ! that include mass of particle in definition of theta
66           DOUBLE PRECISION :: M_PM, D_PM
67     ! small value of theta_m
68           DOUBLE PRECISION :: smallTheta
69     !-----------------------------------------------
70     
71           IER = 0
72           smallTheta = (to_SI)**4 * ZERO_EP_S
73     
74           DO IJK = IJKSTART3, IJKEND3
75             IF ( FLUID_AT(IJK) ) THEN
76     
77               SELECT CASE(KT_TYPE_ENUM)
78                 CASE (LUN_1984, SIMONIN_1996, AHMADI_1995, GD_1999, &
79                       GTSH_2012)
80                   IF (THETA_M(IJK,M) < smallTheta) &
81                      THETA_M(IJK,M) = smallTheta
82     
83                 CASE (IA_2005)
84                   D_PM = D_P(IJK,M)
85                   M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
86                   IF (THETA_M(IJK,M) < smallTheta*M_PM) &
87                     THETA_M(IJK,M) = smallTheta*M_PM
88     
89                 CASE (GHD_2007)
90                   M_PM = ZERO
91                   DO L = 1,SMAX
92                     D_PM = D_P(IJK,L)
93                     M_PM = M_PM +(PI/6.d0)*(D_PM**3)*RO_S(IJK,L)
94                   ENDDO
95                   M_PM = M_PM/DBLE(SMAX)
96                   IF (THETA_M(IJK,M) < smallTheta*M_PM) &
97                     THETA_M(IJK,M) = smallTheta*M_PM
98     
99                 CASE DEFAULT
100     ! should never hit this
101                    WRITE (*, '(A)') 'ADJUST_THETA'
102                    WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
103                    call mfix_exit(myPE)
104               END SELECT   ! end selection of kt_type_enum
105             ENDIF   ! end if (fluid_at)
106           ENDDO  ! end do ijk
107     
108           RETURN
109           END SUBROUTINE ADJUST_THETA
110     
111