File: N:\mfix\model\calc_k_cp.f

1     ! TO DO:
2     ! 1. Kcp needs to be defined for each solids phase (?).
3     ! 2. The part of Kcp from P_star should be based on the
4     !    sum of EP_s of close-packed solids.
5     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
6     !                                                                      C
7     !  Subroutine:  CALC_K_cp                                              C
8     !  Purpose: Calculate and store dPodEp_s                               C
9     !                                                                      C
10     !  Notes: MCP must be defined to call this routine.                    C
11     !                                                                      C
12     !  Author: M. Syamlal                                 Date: 5-FEB-97   C
13     !  Reviewer:                                          Date:            C
14     !                                                                      C
15     !                                                                      C
16     !  Literature/Document References:                                     C
17     !                                                                      C
18     !  Variables referenced:                                               C
19     !  Variables modified:                                                 C
20     !                                                                      C
21     !  Local variables:                                                    C
22     !                                                                      C
23     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
24     
25           SUBROUTINE CALC_K_cp(Kcp)
26     
27     !-----------------------------------------------
28     ! Modules
29     !----------------------------------------------
30           USE compar
31           USE constant
32           USE fldvar
33           USE functions
34           USE geometry
35           USE indices
36           USE param
37           USE param1
38           USE physprop
39           USE pscor
40           USE rdf
41           USE run
42           USE sendrecv
43           USE solids_pressure
44           USE trace
45           USE visc_s
46           IMPLICIT NONE
47     !-----------------------------------------------
48     ! Dummy arguments
49     !-----------------------------------------------
50     ! dPodEP_s
51           DOUBLE PRECISION :: Kcp(DIMENSION_3)
52     !-----------------------------------------------
53     ! Local variables
54     !-----------------------------------------------
55     ! indices
56           INTEGER :: IJK, M
57     ! Other variables
58           DOUBLE PRECISION :: Pc, DPcoDEPs, Mu, Mu_b, Mu_zeta, ZETA
59           DOUBLE PRECISION :: F2, DF2oDEPs, Pf, Pfmax, N_Pff
60     ! Blend Factor
61           Double Precision :: blend
62     !-----------------------------------------------
63     ! External functions
64     !-----------------------------------------------
65           DOUBLE PRECISION :: DZETAoDEPs
66     !-----------------------------------------------
67     
68     ! initializing
69           KCP(:) = ZERO
70     
71           IF (MCP == UNDEFINED_I) THEN
72     ! this error should be caught earlier in the routines so that this
73     ! branch should never be entered
74              RETURN
75           ELSE
76     ! the lowest solids phase index of those solids phases that can close
77     ! pack (i.e. close_packed=T) and the index of the solids phase that is
78     ! used to form the solids correction equation.
79              M = MCP
80           ENDIF
81     
82     
83     ! by definition M must be close_packed (this is a redundant check)
84           IF(CLOSE_PACKED(M)) THEN
85     
86     !!$omp    parallel do default(shared)                  &
87     !!$omp    private( IJK, DPcoDEPS, Pc,                  &
88     !!$omp             Mu, Mu_b, Mu_zeta, ZETA, N_Pff,     &
89     !!$omp             F2, DF2oDEPs, Pf, Pfmax, blend )
90     
91              DO IJK = ijkstart3, ijkend3
92                 IF(.NOT.WALL_AT(IJK))THEN
93     
94                    IF (FRICTION) THEN
95     
96                       IF ((ONE-EP_G(IJK)).GT.EPS_f_min) THEN
97     ! if friction and sufficiently packed to invoke friction
98     ! ---------------------------------------------------------------->>>
99     
100                          IF ((ONE-EP_G(IJK)).GT.(ONE-ep_star_array(ijk))) THEN
101     
102     ! Linearized form of Pc; this is more stable and provides continuous
103     ! function.
104                            DPcoDEPS = (to_SI*Fr)*((delta**5)*&
105                               (2d0*(ONE-ep_star_array(IJK)-delta) - &
106                               2d0*eps_f_min)+&
107                               ((ONE-ep_star_array(ijk)-delta)-eps_f_min)*&
108                               (5*delta**4))/(delta**10)
109     
110                            Pc = (to_SI*Fr)*&
111                               ( ((ONE-ep_star_array(IJK)-delta)-&
112                               EPS_f_min)**N_Pc )/(delta**D_Pc)
113                            Pc = Pc + DPcoDEPS*( (ONE-EP_G(IJK))+delta-&
114                               (ONE-ep_star_array(IJK)))
115     
116     !                        Pc = 1d25*( ((ONE-EP_G(IJK)) - &
117     !                           (ONE-ep_star_array(ijk)))**10d0)
118     !                        DPcoDEPS = 1d26*(((ONE-EP_G(IJK)) - &
119     !                           (ONE-ep_star_array(ijk)))**9d0)
120                          ELSE
121                             Pc = Fr*(((ONE-EP_G(IJK)) - EPS_f_min)**N_Pc)/&
122                                (((ONE-ep_star_array(ijk)) - (ONE-EP_G(IJK))&
123                                + SMALL_NUMBER)**D_Pc)
124                             DPcoDEPs =&
125                                Fr*(((ONE-EP_G(IJK)) - EPS_f_min)**(N_Pc - ONE))&
126                                *(N_Pc*((ONE-ep_star_array(ijk)) - (ONE-EP_G(IJK)))&
127                                +D_Pc*((ONE-EP_G(IJK)) - EPS_f_min))&
128                                / (((ONE-ep_star_array(ijk)) - (ONE-EP_G(IJK)) + &
129                                SMALL_NUMBER)**(D_Pc + ONE))
130                          ENDIF
131     
132                          Mu = (5d0*DSQRT(Pi*Theta_m(IJK,M))&
133     !QX
134                             *D_p(IJK,M)*RO_S(IJK,M))/96d0
135     
136                          Mu_b = (256d0*Mu*EP_s(IJK,M)*EP_s(IJK,M)&
137                             *G_0(IJK,M,M))/(5d0*Pi)
138     
139                          IF (SAVAGE.EQ.1) THEN
140                             Mu_zeta =&
141                                ((2d0+ALPHA)/3d0)*((Mu/(Eta*(2d0-Eta)*&
142                                G_0(IJK,M,M)))*(1d0+1.6d0*Eta*EP_s(IJK,M)*&
143                                G_0(IJK,M,M))*(1d0+1.6d0*Eta*(3d0*Eta-2d0)*&
144                                EP_s(IJK,M)*G_0(IJK,M,M))+(0.6d0*Mu_b*Eta))
145     !QX
146                             ZETA = ((48d0*Eta*(1d0-Eta)*RO_S(IJK,M)*EP_s(IJK,M)*&
147                                 EP_s(IJK,M)*G_0(IJK,M,M)*&
148                                 (Theta_m(IJK,M)**1.5d0))/&
149                                 (SQRT_Pi*D_p(IJK,M)*2d0*Mu_zeta))**0.5d0
150     
151                          ELSEIF (SAVAGE.EQ.0) THEN
152                             ZETA = (SMALL_NUMBER +&
153                                  trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
154                                  trD_s_C(IJK,M))/3.d0))**0.5d0
155     
156                          ELSE
157                             ZETA = ((Theta_m(IJK,M)/(D_p(IJK,M)*D_p(IJK,M))) +&
158                                  (trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
159                                  trD_s_C(IJK,M))/3.d0)))**0.5d0
160     
161                          ENDIF
162     
163                          IF (trD_s_C(IJK,M) .GE. ZERO) THEN
164                             N_Pff = DSQRT(3d0)/(2d0*Sin_Phi) !dilatation
165                          ELSE
166                             N_Pff = N_Pf !compaction
167                          ENDIF
168     
169                          IF ((trD_s_C(IJK,M)/(ZETA*N_Pff*DSQRT(2d0)*&
170                               Sin_Phi)) .GT. 1d0) THEN
171                             F2 = 0d0
172                             DF2oDEPs = ZERO
173                          ELSEIF(trD_s_C(IJK,M) == ZERO) THEN
174                             F2 = ONE
175                             DF2oDEPs = ZERO
176                          ELSE
177                             F2 = (1d0 - (trD_s_C(IJK,M)/(ZETA*N_Pff*&
178                                DSQRT(2d0)*Sin_Phi)))**(N_Pff-1d0)
179                             IF (SAVAGE.EQ.1) THEN
180                                DF2oDEPs = (N_Pff-1d0)*(F2**(N_Pff-2d0))*&
181                                   trD_s_C(IJK,M)&
182                                   *DZETAoDEPs(EP_s(IJK,M), IJK, M)&
183                                   / (ZETA*ZETA*&
184                                   N_Pff*DSQRT(2d0)*Sin_Phi)
185                             ELSE
186                                DF2oDEPs=ZERO
187                             ENDIF
188     
189                             Pf = Pc*F2
190                             Pfmax = Pc*((N_Pf/(N_Pf-1d0))**(N_Pf-1d0))
191     
192                             IF (Pf> Pfmax) THEN
193                                F2 = (N_Pf/(N_Pf-1d0))**(N_Pf-1d0)
194                                DF2oDEPS = ZERO
195                             ENDIF
196                          ENDIF
197     
198     ! Contributions to Kcp(IJK) from kinetic theory have been left out in
199     ! the expressions below as they cause convergence problems at low solids
200     ! volume fraction
201                          Kcp(IJK) = F2*DPcoDEPS + Pc*DF2oDEPS
202     
203                       ELSE
204     ! the solids are not sufficiently packed to invoke friction model
205                          Kcp(IJK) = ZERO
206     
207                       ENDIF   ! end if/else branch (one-ep_g(ijk)>eps_f_min)
208     ! end if friction and sufficiently packed to invoke friction
209     ! ----------------------------------------------------------------<<<
210     
211                    ELSE ! FRICTION = .FALSE.
212     
213     
214                       IF(EP_g(IJK) .LT. ep_g_blend_end(ijk)) THEN
215     ! not friction but the solids are packed so that the plastic pressure
216     ! model is invoked
217     ! ---------------------------------------------------------------->>>
218     
219                          Kcp(IJK) = dPodEP_s(EP_s(IJK, M),ep_g_blend_end(ijk))
220     
221                          IF(BLENDING_STRESS) THEN
222                             blend =  blend_function(IJK)
223                             Kcp(IJK) = (1.0d0-blend) * Kcp(IJK)
224                          ENDIF
225                       ELSE
226     ! the solids are not sufficiently packed to invoke the plastic stress
227     ! model
228                          Kcp(IJK) = ZERO
229                       ENDIF
230     ! end if not friction but sufficiently packed to invoke plastic pressure
231     ! ----------------------------------------------------------------<<<
232     
233                    ENDIF   ! end if/else branch if(friction)
234     
235                 ELSE    ! else branch of if (.not.wall_at(ijk))
236                    Kcp(IJK) = ZERO
237                 ENDIF   ! end if/else (.not.wall_at(ijk))
238     
239              ENDDO  ! do ijk=ijkstart3,ijkend3
240           ENDIF   ! end if (close_packed(m))
241     
242           CALL send_recv(Kcp, 2)
243     
244           RETURN
245           END
246     
247     
248     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
249     !                                                                      C
250     !  Function: DZETAoDEPs (EPs,IJK,M)                                    C
251     !  Purpose: Calculate derivative of zeta                               C
252     !           w.r.t granular volume fraction                             C
253     !                                                                      C
254     !  Author: A. Srivastava                              Date: 8-JUNE-98  C
255     !  Reviewer:                                          Date:            C
256     !                                                                      C
257     !  Revision Number:                                                    C
258     !  Purpose:                                                            C
259     !  Author:                                            Date: dd-mmm-yy  C
260     !  Reviewer:                                          Date: dd-mmm-yy  C
261     !                                                                      C
262     !  Literature/Document References:                                     C
263     !                                                                      C
264     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
265     
266           DOUBLE PRECISION FUNCTION DZETAoDEPs(EPs, IJK, M)
267     
268     !-----------------------------------------------
269     ! Modules
270     !-----------------------------------------------
271           USE compar
272           USE constant
273           USE fldvar
274           USE functions
275           USE geometry
276           USE indices
277           USE param
278           USE param1
279           USE physprop
280           USE pscor
281           USE rdf
282           USE run
283           USE trace
284           USE visc_s
285           IMPLICIT NONE
286     !-----------------------------------------------
287     ! Dummy Arguments
288     !-----------------------------------------------
289     ! solids volume fraction
290           DOUBLE PRECISION, INTENT(IN) :: EPs
291     ! indices
292           INTEGER, INTENT(IN) :: IJK,M
293     !-----------------------------------------------
294     ! Local variables
295     !-----------------------------------------------
296     ! local radial distribution function
297           DOUBLE PRECISION :: g0
298     ! Other variables
299           DOUBLE PRECISION :: Mu, Mu_b,  DEPs2G_0oDEPs, F1, DF1oDEPs
300     !-----------------------------------------------
301     
302           g0 = G_0(IJK, M, M)
303     !QX
304           Mu = (5d0*DSQRT(Pi*Theta_m(IJK,M))*D_p(IJK,M)*RO_S(IJK,M))/96d0
305     
306           Mu_b = (256d0*Mu*EPs*EPs*g0/(5d0*Pi))
307     
308           DEPs2G_0oDEPs = EPs*EPs*DG_0DNU(EPs) + 2d0*EPs*g0
309     
310           F1 = ((2d0+ALPHA)/3d0)*((2*Mu/(Eta*(2d0-Eta)*&
311                g0))*(1d0+1.6d0*Eta*EPs*g0)*(1d0+1.6d0*Eta*(3d0*Eta-2d0)&
312                *EPs*g0)+(1.2d0*Mu_b*Eta))
313     
314           DF1oDEPs = ((2d0+ALPHA)/3d0)*((2*Mu/(Eta*(2d0-Eta))*&
315              ((-DG_0DNU(EPs)/(g0*g0)) + (1.6d0*Eta*(3d0*Eta-1d0))&
316             + (64d0*Eta*Eta*(3d0*Eta-2d0)*DEPs2G_0oDEPs/25d0))) +&
317             3.2d0*Eta*RO_S(IJK,M)*D_p(IJK,M)*((Theta_m(IJK,M)/Pi)**0.5d0)&
318               *DEPs2G_0oDEPs)
319     
320           DZETAoDEPs = 0.5d0*((48d0*Eta*(1d0-Eta)*RO_S(IJK,M)*F1*&
321                   (Theta_m(IJK,M)**1.5d0)/&
322                (SQRT_Pi*D_p(IJK,M)*EPs*EPs*g0))**0.5d0)*&
323                (F1*DEPs2G_0oDEPs - EPs*Eps*g0*DF1oDEPs)/(F1*F1)
324     
325     
326           RETURN
327           END
328     
329