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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_Gw_Hw_Cw                                           C
4     !                                                                      C
5     !  Purpose: Calculate Gw, Hw, and Cw                                   C
6     !                                                                      C
7     !  Author: A. Srivastava & K. Agrawal, Princeton Univ. Date: 10-APR-98 C
8     !  Reviewer:                                           Date:           C
9     !                                                                      C
10     !  Modified: Sofiane Benyahia, Fluent Inc.             Date: 03-FEB-05 C
11     !  Purpose: Include conductivity defined by Simonin and Ahmadi         C
12     !           Also included Jenkins small frictional limit               C
13     !                                                                      C
14     !  Literature/Document References:                                     C
15     !     Johnson, P. C., and Jackson, R., Frictional-collisional          C
16     !        constitutive relations for granular materials, with           C
17     !        application to plane shearing, Journal of Fluid Mechanics,    C
18     !        1987, 176, pp. 67-93                                          C
19     !     Jenkins, J. T., and Louge, M. Y., On the flux of fluctuating     C
20     !        energy in a collisional grain flow at a flat frictional wall, C
21     !        Physics of Fluids, 1997, 9(10), pp. 2835-2840                 C
22     !                                                                      C
23     !     See calc_mu_s.f for references on frictional theory models       C
24     !     See calc_mu_s.f for references on Ahmadi and Simonin models      C
25     !                                                                      C
26     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
27     
28           SUBROUTINE CALC_Gw_Hw_Cw(g0, EPS, EPG, ep_star_avg, &
29                 g0EP_avg, TH, Mu_g_avg, RO_g_avg, ROs_avg, &
30                 DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, VREL, VSLIP,&
31                 DEL_U, S_S, S_dd, VEL, W_VEL, M, gw, hw, cw)
32     
33     !-----------------------------------------------
34     ! Modules
35     !-----------------------------------------------
36           USE param
37           USE param1
38           USE constant
39           USE physprop
40           USE fldvar
41           USE bc
42           USE run
43           USE mpi_utility
44           IMPLICIT NONE
45     !-----------------------------------------------
46     ! Dummy Arguments
47     !-----------------------------------------------
48     ! Radial distribution function of solids phase M with each
49     ! other solids phase
50           DOUBLE PRECISION, INTENT(IN) :: g0
51     ! Average solids volume fraction of each solids phase
52           DOUBLE PRECISION, INTENT(IN) :: EPS
53     ! Average solids and gas volume fraction
54           DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
55     ! Sum of eps*G_0
56           DOUBLE PRECISION, INTENT(IN) :: g0EP_avg
57     ! Average theta_m
58           DOUBLE PRECISION, INTENT(INOUT) :: TH
59     ! Average gas viscosity
60           DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
61     ! Average gas density
62           DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
63     ! Average solids density
64           DOUBLE PRECISION, INTENT(IN) :: ROS_avg
65     ! Average particle diameter of each solids phase
66           DOUBLE PRECISION, INTENT(IN) :: DP_avg
67     ! Average cross-correlation K_12 and lagrangian integral time-scale
68           DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg, Tau_1_avg
69     ! Magnitude of slip velocity between two phases
70           DOUBLE PRECISION, INTENT(IN) :: VREL
71     ! Slip velocity between wall and particles
72           DOUBLE PRECISION, INTENT(IN) :: VSLIP
73     ! Relevant solids velocity at wall
74           DOUBLE PRECISION, INTENT(IN) :: VEL
75     ! Relevant wall velocity
76           DOUBLE PRECISION, INTENT(IN) :: W_VEL
77     ! del.u
78           DOUBLE PRECISION, INTENT(IN) :: DEL_U
79     ! S:S
80           DOUBLE PRECISION, INTENT(IN) :: S_S
81     ! S_dd (d can be x,y or z)
82           DOUBLE PRECISION, INTENT(IN) :: S_dd
83     ! Solids phase index
84           INTEGER, INTENT(IN) :: M
85     ! Wall momentum coefficients:
86     ! 1st term on LHS
87           DOUBLE PRECISION, INTENT(INOUT) :: Gw
88     ! 2nd term on LHS
89           DOUBLE PRECISION, INTENT(INOUT) :: Hw
90     ! all terms appearing on RHS
91           DOUBLE PRECISION, INTENT(INOUT) :: Cw
92     !-----------------------------------------------
93     ! Local Variables
94     !-----------------------------------------------
95     ! Part of wall momentum coefficient in 1st term of LHS
96           DOUBLE PRECISION :: Mu_s
97     ! Term appearing in wall momenutm coefficient of 2nd term LHS
98           DOUBLE PRECISION :: F_2
99     ! Frictional Pressure Pf
100           DOUBLE PRECISION :: Pf
101     ! Critical pressure Pc
102           DOUBLE PRECISION :: Pc
103     ! Chi (appears in frictional boundary condition)
104           DOUBLE PRECISION :: Chi, N_Pff
105     ! Viscosity
106           DOUBLE PRECISION :: Mu
107     ! Bulk viscosity
108           DOUBLE PRECISION :: Mu_b
109     ! Viscosity corrected for interstitial fluid effects
110           DOUBLE PRECISION :: Mu_star
111     ! Reynolds number based on slip velocity
112           DOUBLE PRECISION :: Re_g
113     ! Friction Factor in drag coefficient
114           DOUBLE PRECISION :: C_d
115     ! Drag Coefficient
116           DOUBLE PRECISION :: Beta, DgA
117     ! Square root of S:S or the form suggested by Savage
118           DOUBLE PRECISION :: ZETA
119     ! Constants in Simonin model
120           DOUBLE PRECISION :: Sigma_c, Tau_2_c, Tau_12_st, Nu_t
121           DOUBLE PRECISION :: Tau_2, zeta_c_2, MU_2_T_Kin, Mu_2_Col
122           DOUBLE PRECISION :: Tmp_Ahmadi_Const
123     ! Other local terms
124           DOUBLE PRECISION :: dpc_dphi
125     ! Local variable for rdf
126           DOUBLE PRECISION :: g0_loc
127     !-----------------------------------------------
128     
129     ! This is done here similar to bc_theta to avoid small negative values of
130     ! Theta coming most probably from linear solver
131           IF(TH .LE. ZERO)THEN
132             TH = 1D-8
133             IF (myPE.eq.PE_IO) THEN
134                WRITE(*,*) &
135                   'Warning: Negative granular temp at wall set to 1e-8'
136     !          CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
137             ENDIF
138           ENDIF
139     
140           g0_loc = g0
141     
142     ! modify F_2 if Jenkins BC is used (sof)
143           IF(JENKINS) THEN
144     
145             IF (VSLIP == ZERO) THEN
146     ! if solids velocity field is initialized to zero, use free slip bc
147                F_2 = zero
148     
149             ELSE
150                IF(AHMADI) THEN
151     ! Ahmadi model uses different solids pressure model
152     ! the coefficient mu in Jenkins paper is defined as tan_Phi_w, that's how
153     ! I understand it from soil mechanic papers, i.e., G.I. Tardos, powder
154     ! Tech. 92 (1997), 61-74. See his equation (1). Define Phi_w in mfix.dat!
155     ! here F_2 divided by VSLIP to use the same bc as Johnson&Jackson
156     
157                   F_2 = tan_Phi_w*ROS_avg*EPS* &
158                      ((ONE + 4.0D0*g0EP_avg) + HALF*(ONE -C_e*C_e))*TH/VSLIP
159     
160                ELSE
161     ! Simonin or granular models use same solids pressure
162                   F_2 = tan_Phi_w*ROS_avg*EPS*(1d0+ 4.D0 * Eta *g0EP_avg)*TH/VSLIP
163                ENDIF   ! end if(Ahmadi)
164     
165             ENDIF ! endif(vslip==0)
166     
167           ELSE    ! if(.not.jenkins)
168     
169              F_2 = (PHIP*DSQRT(3d0*TH)*Pi*ROS_avg*EPS*g0_loc)&
170                   /(6d0*(ONE-ep_star_avg))
171     
172           ENDIF   ! end if(Jenkins)/else
173     
174     
175           Mu = (5d0*DSQRT(Pi*TH)*Dp_avg*ROS_avg)/96d0
176           Mu_b = (256d0*Mu*EPS*g0EP_avg)/(5d0*Pi)
177     
178     ! This is from Wen-Yu correlation, you can put here your own single particle drag
179           Re_g = EPG*RO_g_avg*Dp_avg*VREL/Mu_g_avg
180           IF (Re_g.lt.1000d0) THEN
181              C_d = (24.d0/(Re_g+SMALL_NUMBER))*(ONE + 0.15d0 * Re_g**0.687d0)
182           ELSE
183              C_d = 0.44d0
184           ENDIF
185           DgA = 0.75d0*C_d*Ro_g_avg*EPG*VREL/(Dp_avg*EPG**(2.65d0))
186           IF(VREL == ZERO) DgA = LARGE_NUMBER
187           Beta = SWITCH*EPS*DgA
188     
189     ! SWITCH enables us to turn on/off the modification to the
190     ! particulate phase viscosity. If we want to simulate gas-particle
191     ! flow then SWITCH=1 to incorporate the effect of drag on the
192     ! particle viscosity. If we want to simulate granular flow
193     ! without the effects of an interstitial gas, SWITCH=0.
194           IF(SWITCH == ZERO .OR. Ro_g_avg == ZERO)THEN
195              Mu_star = Mu
196           ELSEIF(TH .LT. SMALL_NUMBER)THEN
197              MU_star = ZERO
198           ELSE
199              Mu_star = ROS_avg*EPS* g0_loc*TH* Mu/ &
200                 (ROS_avg*g0EP_avg*TH + 2.0d0*SWITCH*DgA/ROS_avg* Mu)
201           ENDIF
202     
203           Mu_s = ((2d0+ALPHA)/3d0)*((Mu_star/(Eta*(2d0-Eta)*&
204                        g0_loc))*(ONE+1.6d0*Eta*g0EP_avg&
205                        )*(ONE+1.6d0*Eta*(3d0*Eta-2d0)*&
206                        g0EP_avg)+(0.6d0*Mu_b*Eta))
207     
208     ! particle relaxation time
209           Tau_12_st = ROS_avg/(DgA+small_number)
210     
211           IF(SIMONIN) THEN !see calc_mu_s for explanation of these definitions
212     
213              Sigma_c = (ONE+ C_e)*(3.d0-C_e)/5.d0
214              Tau_2_c = DP_avg/(6.d0*EPS*g0_loc*DSQRT(16.d0*(TH+Small_number)/PI))
215              zeta_c_2= 2.D0/5.D0*(ONE+ C_e)*(3.d0*C_e-ONE)
216              Nu_t =  Tau_12_avg/Tau_12_st
217              Tau_2 = ONE/(2.D0/Tau_12_st+Sigma_c/Tau_2_c)
218              MU_2_T_Kin = (2.0D0/3.0D0*K_12_avg*Nu_t + TH * &
219                          (ONE+ zeta_c_2*EPS*g0_loc))*Tau_2
220              Mu_2_Col = 8.D0/5.D0*EPS*g0_loc*Eta* (MU_2_T_Kin+ &
221                        Dp_avg*DSQRT(TH/PI))
222              Mu_s = EPS*ROS_avg*(MU_2_T_Kin + Mu_2_Col)
223           ELSEIF(AHMADI) THEN
224              IF(EPS < (ONE-ep_star_avg)) THEN
225                 Tmp_Ahmadi_Const = &
226                 ONE/(ONE+ Tau_1_avg/Tau_12_st * (ONE-EPS/(ONE-ep_star_avg))**3)
227              ELSE
228                 Tmp_Ahmadi_Const = ONE
229              ENDIF
230              Mu_s = Tmp_Ahmadi_Const &
231                  *0.1045D0*(ONE/g0_loc+3.2D0*EPS+12.1824D0*g0_loc*EPS*EPS)  &
232                  *Dp_avg*ROS_avg* DSQRT(TH)
233           ENDIF
234     
235     ! Calculating frictional terms
236           IF ((ONE-EPG)<= EPS_f_min) THEN
237              Pf = ZERO
238              Chi = ZERO
239              ZETA = 1d0
240           ELSE
241              IF (SAVAGE.EQ.1) THEN    !form of Savage
242                 ZETA = ((48d0*Eta*(1d0-Eta)*ROS_avg*EPS*EPS*g0_loc*&
243                         (TH**1.5d0))/&
244                         (SQRT_Pi*Dp_avg*2d0*Mu_s))**0.5d0
245              ELSEIF (SAVAGE.EQ.0)  THEN !S:S form
246                 ZETA = DSQRT(S_S)
247              ELSE
248                 ZETA = DSQRT(S_S + (TH/(Dp_avg*Dp_avg)))
249              ENDIF
250     
251              IF (EPG < ep_star_avg) THEN
252                 dpc_dphi = (to_SI*Fr)*((delta**5)*(2d0*(ONE-ep_star_avg-delta) - &
253                    2d0*eps_f_min)+((ONE-ep_star_avg-delta)-eps_f_min)&
254                    *(5*delta**4))/(delta**10)
255     
256                 Pc = (to_SI*Fr)*(((ONE-ep_star_avg-delta) - EPS_f_min)**N_Pc)/(delta**D_Pc)
257     !            Pc=  1d25*(((ONE-EPG)-(ONE-ep_star_avg))**10d0)  ! this is old Pc
258              ELSE
259                 Pc = Fr*(((ONE-EPG) - EPS_f_min)**N_Pc)/ &
260                    (((ONE-ep_star_avg) - (ONE-EPG) + SMALL_NUMBER)**D_Pc)
261              ENDIF
262     
263              IF (DEL_U .GE. ZERO) THEN
264                 N_Pff = DSQRT(3d0)/(2d0*Sin_Phi) !dilatation
265              ELSE
266                 N_Pff = N_Pf !compaction
267              ENDIF
268     
269              IF ((DEL_U/(ZETA*N_Pff*DSQRT(2d0)*Sin_Phi)) .GT. 1d0) THEN
270                 Pf = ZERO
271              ELSEIF( DEL_U == ZERO ) THEN
272                 Pf = Pc
273              ELSE
274                 Pf = Pc*(1d0 - (DEL_U/(ZETA*N_Pff*DSQRT(2d0)*Sin_Phi)))**&
275                    (N_Pff-1d0)
276              ENDIF
277     
278              Chi =  DSQRT(2d0)*Pf*Sin_Phi*(N_Pff - (N_Pff-1d0)*&
279                 ((Pf/Pc)**(1d0/(N_Pff-1d0))))
280     
281              IF (Chi< ZERO) THEN
282                 Pf = Pc*((N_Pff/(N_Pff-1d0))**(N_Pff-1d0))
283                 Chi = ZERO
284              ENDIF
285     
286     ! by writing Pf & Chi in the following form, we ensure distribution
287     ! of stresses amoung all solids phases (sof, Oct 24 2005)
288     
289              Pf = Pf * EPS/(ONE-EPG)
290              Chi = Chi * EPS/(ONE-EPG)
291     
292           ENDIF
293     
294     ! Calculating gw, hw, cw
295     
296           Gw = ONE   ! we write this in the exact same form as non-frictional JJ BC
297     
298           IF(VSLIP == ZERO .OR. ZETA == ZERO) THEN
299             Hw = ZERO
300           ELSE
301             Hw = F_2 + Pf*tan_Phi_w - Chi*S_dd*tan_Phi_w/ZETA
302             Hw = Hw / (MU_s + Chi/(2d0*ZETA))  ! this is because Gw is set to one.
303           ENDIF
304           Cw = hw * W_VEL
305     
306           RETURN
307           END SUBROUTINE CALC_Gw_Hw_Cw
308