calc_k_cp.f

Go to the documentation of this file.

! TO DO:
! 1. Kcp needs to be defined for each solids phase (?).
! 2. The part of Kcp from P_star should be based on the
!    sum of EP_s of close-packed solids.
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine:  CALC_K_cp                                              C
!  Purpose: Calculate and store dPodEp_s                               C
!                                                                      C
!  Notes: MCP must be defined to call this routine.                    C
!                                                                      C
!  Author: M. Syamlal                                 Date: 5-FEB-97   C
!  Reviewer:                                          Date:            C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!                                                                      C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE calc_k_cp(Kcp)

!-----------------------------------------------
! Modules
!----------------------------------------------
      USE compar
      USE constant
      USE fldvar
      USE functions
      USE geometry
      USE indices
      USE param
      USE param1
      USE physprop
      USE pscor
      USE rdf
      USE run
      USE sendrecv
      USE solids_pressure
      USE trace
      USE visc_s
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! dPodEP_s
      DOUBLE PRECISION :: Kcp(dimension_3)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! indices
      INTEGER :: IJK, M
! Other variables
      DOUBLE PRECISION :: Pc, DPcoDEPs, Mu, Mu_b, Mu_zeta, ZETA
      DOUBLE PRECISION :: F2, DF2oDEPs, Pf, Pfmax, N_Pff
! Blend Factor
      Double Precision :: blend
!-----------------------------------------------
! External functions
!-----------------------------------------------
      DOUBLE PRECISION :: DZETAoDEPs
!-----------------------------------------------

! initializing
      kcp(:) = zero

      IF (mcp == undefined_i) THEN
! this error should be caught earlier in the routines so that this
! branch should never be entered
         RETURN
      ELSE
! the lowest solids phase index of those solids phases that can close
! pack (i.e. close_packed=T) and the index of the solids phase that is
! used to form the solids correction equation.
         m = mcp
      ENDIF


! by definition M must be close_packed (this is a redundant check)
      IF(close_packed(m)) THEN

!!$omp    parallel do default(shared)                  &
!!$omp    private( IJK, DPcoDEPS, Pc,                  &
!!$omp             Mu, Mu_b, Mu_zeta, ZETA, N_Pff,     &
!!$omp             F2, DF2oDEPs, Pf, Pfmax, blend )

         DO ijk = ijkstart3, ijkend3
            IF(.NOT.wall_at(ijk))THEN

               IF (friction) THEN

                  IF ((one-ep_g(ijk)).GT.eps_f_min) THEN
! if friction and sufficiently packed to invoke friction
! ---------------------------------------------------------------->>>

                     IF ((one-ep_g(ijk)).GT.(one-ep_star_array(ijk))) THEN

! Linearized form of Pc; this is more stable and provides continuous
! function.
                       dpcodeps = (to_si*fr)*((delta**5)*&
                          (2d0*(one-ep_star_array(ijk)-delta) - &
                          2d0*eps_f_min)+&
                          ((one-ep_star_array(ijk)-delta)-eps_f_min)*&
                          (5*delta**4))/(delta**10)

                       pc = (to_si*fr)*&
                          ( ((one-ep_star_array(ijk)-delta)-&
                          eps_f_min)**n_pc )/(delta**d_pc)
                       pc = pc + dpcodeps*( (one-ep_g(ijk))+delta-&
                          (one-ep_star_array(ijk)))

!                        Pc = 1d25*( ((ONE-EP_G(IJK)) - &
!                           (ONE-ep_star_array(ijk)))**10d0)
!                        DPcoDEPS = 1d26*(((ONE-EP_G(IJK)) - &
!                           (ONE-ep_star_array(ijk)))**9d0)
                     ELSE
                        pc = fr*(((one-ep_g(ijk)) - eps_f_min)**n_pc)/&
                           (((one-ep_star_array(ijk)) - (one-ep_g(ijk))&
                           + small_number)**d_pc)
                        dpcodeps =&
                           fr*(((one-ep_g(ijk)) - eps_f_min)**(n_pc - one))&
                           *(n_pc*((one-ep_star_array(ijk)) - (one-ep_g(ijk)))&
                           +d_pc*((one-ep_g(ijk)) - eps_f_min))&
                           / (((one-ep_star_array(ijk)) - (one-ep_g(ijk)) + &
                           small_number)**(d_pc + one))
                     ENDIF

                     mu = (5d0*dsqrt(pi*theta_m(ijk,m))&
!QX
                        *d_p(ijk,m)*ro_s(ijk,m))/96d0

                     mu_b = (256d0*mu*ep_s(ijk,m)*ep_s(ijk,m)&
                        *g_0(ijk,m,m))/(5d0*pi)

                     IF (savage.EQ.1) THEN
                        mu_zeta =&
                           ((2d0+alpha)/3d0)*((mu/(eta*(2d0-eta)*&
                           g_0(ijk,m,m)))*(1d0+1.6d0*eta*ep_s(ijk,m)*&
                           g_0(ijk,m,m))*(1d0+1.6d0*eta*(3d0*eta-2d0)*&
                           ep_s(ijk,m)*g_0(ijk,m,m))+(0.6d0*mu_b*eta))
!QX
                        zeta = ((48d0*eta*(1d0-eta)*ro_s(ijk,m)*ep_s(ijk,m)*&
                            ep_s(ijk,m)*g_0(ijk,m,m)*&
                            (theta_m(ijk,m)**1.5d0))/&
                            (sqrt_pi*d_p(ijk,m)*2d0*mu_zeta))**0.5d0

                     ELSEIF (savage.EQ.0) THEN
                        zeta = (small_number +&
                             trd_s2(ijk,m) - ((trd_s_c(ijk,m)*&
                             trd_s_c(ijk,m))/3.d0))**0.5d0

                     ELSE
                        zeta = ((theta_m(ijk,m)/(d_p(ijk,m)*d_p(ijk,m))) +&
                             (trd_s2(ijk,m) - ((trd_s_c(ijk,m)*&
                             trd_s_c(ijk,m))/3.d0)))**0.5d0

                     ENDIF

                     IF (trd_s_c(ijk,m) .GE. zero) THEN
                        n_pff = dsqrt(3d0)/(2d0*sin_phi) !dilatation
                     ELSE
                        n_pff = n_pf !compaction
                     ENDIF

                     IF ((trd_s_c(ijk,m)/(zeta*n_pff*dsqrt(2d0)*&
                          sin_phi)) .GT. 1d0) THEN
                        f2 = 0d0
                        df2odeps = zero
                     ELSEIF(trd_s_c(ijk,m) == zero) THEN
                        f2 = one
                        df2odeps = zero
                     ELSE
                        f2 = (1d0 - (trd_s_c(ijk,m)/(zeta*n_pff*&
                           dsqrt(2d0)*sin_phi)))**(n_pff-1d0)
                        IF (savage.EQ.1) THEN
                           df2odeps = (n_pff-1d0)*(f2**(n_pff-2d0))*&
                              trd_s_c(ijk,m)&
                              *dzetaodeps(ep_s(ijk,m), ijk, m)&
                              / (zeta*zeta*&
                              n_pff*dsqrt(2d0)*sin_phi)
                        ELSE
                           df2odeps=zero
                        ENDIF

                        pf = pc*f2
                        pfmax = pc*((n_pf/(n_pf-1d0))**(n_pf-1d0))

                        IF (pf> pfmax) THEN
                           f2 = (n_pf/(n_pf-1d0))**(n_pf-1d0)
                           df2odeps = zero
                        ENDIF
                     ENDIF

! Contributions to Kcp(IJK) from kinetic theory have been left out in
! the expressions below as they cause convergence problems at low solids
! volume fraction
                     kcp(ijk) = f2*dpcodeps + pc*df2odeps

                  ELSE
! the solids are not sufficiently packed to invoke friction model
                     kcp(ijk) = zero

                  ENDIF   ! end if/else branch (one-ep_g(ijk)>eps_f_min)
! end if friction and sufficiently packed to invoke friction
! ----------------------------------------------------------------<<<

               ELSE ! FRICTION = .FALSE.


                  IF(ep_g(ijk) .LT. ep_g_blend_end(ijk)) THEN
! not friction but the solids are packed so that the plastic pressure
! model is invoked
! ---------------------------------------------------------------->>>

                     kcp(ijk) = dpodep_s(ep_s(ijk, m),ep_g_blend_end(ijk))

                     IF(blending_stress) THEN
                        blend =  blend_function(ijk)
                        kcp(ijk) = (1.0d0-blend) * kcp(ijk)
                     ENDIF
                  ELSE
! the solids are not sufficiently packed to invoke the plastic stress
! model
                     kcp(ijk) = zero
                  ENDIF
! end if not friction but sufficiently packed to invoke plastic pressure
! ----------------------------------------------------------------<<<

               ENDIF   ! end if/else branch if(friction)

            ELSE    ! else branch of if (.not.wall_at(ijk))
               kcp(ijk) = zero
            ENDIF   ! end if/else (.not.wall_at(ijk))

         ENDDO  ! do ijk=ijkstart3,ijkend3
      ENDIF   ! end if (close_packed(m))

      CALL send_recv(kcp, 2)

      RETURN
      END


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Function: DZETAoDEPs (EPs,IJK,M)                                    C
!  Purpose: Calculate derivative of zeta                               C
!           w.r.t granular volume fraction                             C
!                                                                      C
!  Author: A. Srivastava                              Date: 8-JUNE-98  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number:                                                    C
!  Purpose:                                                            C
!  Author:                                            Date: dd-mmm-yy  C
!  Reviewer:                                          Date: dd-mmm-yy  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      DOUBLE PRECISION FUNCTION dzetaodeps(EPs, IJK, M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE compar
      USE constant
      USE fldvar
      USE functions
      USE geometry
      USE indices
      USE param
      USE param1
      USE physprop
      USE pscor
      USE rdf
      USE run
      USE trace
      USE visc_s
      IMPLICIT NONE
!-----------------------------------------------
! Dummy Arguments
!-----------------------------------------------
! solids volume fraction
      DOUBLE PRECISION, INTENT(IN) :: EPs
! indices
      INTEGER, INTENT(IN) :: IJK,M
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! local radial distribution function
      DOUBLE PRECISION :: g0
! Other variables
      DOUBLE PRECISION :: Mu, Mu_b,  DEPs2G_0oDEPs, F1, DF1oDEPs
!-----------------------------------------------

      g0 = g_0(ijk, m, m)
!QX
      mu = (5d0*dsqrt(pi*theta_m(ijk,m))*d_p(ijk,m)*ro_s(ijk,m))/96d0

      mu_b = (256d0*mu*eps*eps*g0/(5d0*pi))

      deps2g_0odeps = eps*eps*dg_0dnu(eps) + 2d0*eps*g0

      f1 = ((2d0+alpha)/3d0)*((2*mu/(eta*(2d0-eta)*&
           g0))*(1d0+1.6d0*eta*eps*g0)*(1d0+1.6d0*eta*(3d0*eta-2d0)&
           *eps*g0)+(1.2d0*mu_b*eta))

      df1odeps = ((2d0+alpha)/3d0)*((2*mu/(eta*(2d0-eta))*&
         ((-dg_0dnu(eps)/(g0*g0)) + (1.6d0*eta*(3d0*eta-1d0))&
        + (64d0*eta*eta*(3d0*eta-2d0)*deps2g_0odeps/25d0))) +&
        3.2d0*eta*ro_s(ijk,m)*d_p(ijk,m)*((theta_m(ijk,m)/pi)**0.5d0)&
          *deps2g_0odeps)

      dzetaodeps = 0.5d0*((48d0*eta*(1d0-eta)*ro_s(ijk,m)*f1*&
              (theta_m(ijk,m)**1.5d0)/&
           (sqrt_pi*d_p(ijk,m)*eps*eps*g0))**0.5d0)*&
           (f1*deps2g_0odeps - eps*eps*g0*df1odeps)/(f1*f1)


      RETURN
      END