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