des_thermo_cond_mod.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: DES_CONDUCTION                                         !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!  Author: J.Musser                                   Date: 25-Jun-10  !
!                                                                      !
!  Comments:                                                           !
!                                                                      !
!  REF: Batchelor and O'Brien, "Thermal or electrical conduction       !
!       through a granular material," Proceedings of the Royal Society !
!       of London. Series A, Mathematical and Physical Sciences,       !
!       Vol. 355, no 1682, pp. 313-333, July 1977.                     !
!                                                                      !
!  REF: Rong and Horio, "DEM simulation of char combustion in a        !
!       fluidized bed," in Second International Conference on CFD in   !
!       the materials and Process Industries, Melbourne, 1999, pp.     !
!       65-70.                                                         !
!                                                                      !
!  REF: Xavier and Davidson, "Heat transfer to surfaces immersed in    !
!       fluidised beds, particularly tub arrays," in Fluidization:     !
!       Proceedings of the Second Engineering Foundation Conference,   !
!       Cambridge, England, 2-6 April, 1978, pp. 333-338.              !
!                                                                      !
!  REF: Zhou, Yu, and Zulli, "Particle scale study of heat transfer in !
!       packed and bubbling fluidized beds," AIChE Journal, Vol. 55,   !
!       no 4, pp 868-884, 2009.                                        !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      MODULE des_thermo_cond
 

      IMPLICIT NONE
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: des_qw_cond
      LOGICAL :: do_area_correction
      CONTAINS



      FUNCTION des_conduction(I, J, CENTER_DIST, iM, iIJK)

      use constant
      use des_thermo
      use discretelement
      use funits
      use param1, only: zero, undefined
      use physprop
      use run, only: units
      IMPLICIT NONE

! Passed variables
!---------------------------------------------------------------------//
! Index of particle being looped over
      INTEGER, INTENT(IN) :: I,J
! Index of fluid cell containing particle I
      INTEGER, INTENT(IN) :: iIJK
! Solid phase indices for the given particles
      INTEGER, INTENT(IN) :: iM
! Distance between the centers of particle I and particle J (magnitude)
      DOUBLE PRECISION, INTENT(IN) :: CENTER_DIST
      DOUBLE PRECISION :: DES_CONDUCTION

! Local variables
!---------------------------------------------------------------------//
! Solid phase indices for the given particles
      INTEGER jM
! Radius of smaller particle
      DOUBLE PRECISION MIN_RAD
! Radius of larger particle
      DOUBLE PRECISION MAX_RAD
! Rate of particle-particle conduction
      DOUBLE PRECISION Q_pp
! Rate of particle-fluid-particle conduction
      DOUBLE PRECISION Q_pfp
! Outer radius of region delineating particle-fluid-particle conduction
      DOUBLE PRECISION RD_OUT
! Inner radius of region delineating particle-fluid-particle conduction
      DOUBLE PRECISION RD_IN
! The radius of the fluid lens containing the larger particle
      DOUBLE PRECISION LENS_RAD
! Temperature difference between two particles
      DOUBLE PRECISION DeltaTp
! Effective thermal conductivity
      DOUBLE PRECISION lK_eff
! Radius of contact area
      DOUBLE PRECISION lRadius
! Particle overlap (simulated and corrected)
      DOUBLE PRECISION OLAP_Sim, OLAP_actual
! Corrected center distance
      DOUBLE PRECISION CENTER_DIST_CORR
! Gas thermal conductivity
      DOUBLE PRECISION :: K_gas
! Overlap distance between each particle and contact plane
      DOUBLE PRECISION :: OLAP_1, OLAP_2
! Effective Lens for particles 1 and 2 (req'd for analytic calculation)
      DOUBLE PRECISION :: RLENS_1, RLENS_2
! Effective minimum conduction distance for particles 1 and 2
      DOUBLE PRECISION :: S_1, S_2
! Effective thermal conductance for particles 1 and 2
      DOUBLE PRECISION :: H_1, H_2
! Analytic thermal conductance between particles
      DOUBLE PRECISION :: H
! Dummy variable for calculations 
      DOUBLE PRECISION :: RATIO
! Functions
!---------------------------------------------------------------------//
 !     DOUBLE PRECISION :: EVAL_H_PFP

! Identify the solid phases of the neighbor particle
         jm = pijk(j,5)

! Determine the radius of the larger and smaller particle
         min_rad = min(des_radius(i), des_radius(j))
         max_rad = max(des_radius(i), des_radius(j))

         deltatp = des_t_s(j) - des_t_s(i)

! Calculate the particle-particle conduction
! REF: Batchelor and O'Brien, 1977 (MODIFIED)
!---------------------------------------------------------------------//
         q_pp = 0.0

! Initialize corrected center-center distance
         center_dist_corr = center_dist
         IF(center_dist < (max_rad + min_rad)) THEN
! Correct particle overlap to account for artificial softening
            olap_sim = max_rad + min_rad - center_dist
            IF(do_area_correction)THEN
               IF(im.ge.jm)THEN
                  olap_actual = correct_olap(olap_sim,jm,im)
               ELSE
                  olap_actual = correct_olap(olap_sim,im,jm)
               ENDIF
               center_dist_corr = max_rad + min_rad - olap_actual
            ENDIF
! Effective thermal conductivity
            lk_eff = k_eff(k_s0(im),k_s0(jm))
! Effective contact area's radius
            lradius = radius(max_rad, min_rad)
! Compute time-correction term (not done yet)
            ! CALL CALC_TIME_CORRECTION (.... )
! Inter-particle heat transfer
            q_pp = 2.0d0 * lk_eff * lradius * deltatp

! Assign the inter-particle heat transfer to both particles.
         ENDIF

!!         IF(.TRUE.)THEN ! NEW WAY
!---------------------------------------------------------------------//
! Calculate the particle-fluid-particle conduction
! REF: Rong and Horio, 1999 (MODIFIED)
!---------------------------------------------------------------------//
         lens_rad = max_rad * (1.0d0 + flpc)
         olap_actual = max_rad + min_rad - center_dist_corr
! check to see if particles are within lens thickness from eachother
! a negative overlap indicates particles are separated 
         IF(olap_actual > (-max_rad*flpc))THEN
            rd_out = radius(lens_rad, min_rad)
            ! get overlap between each particle and contact plane
            ratio = olap_actual / max_rad
            olap_1 = max_rad*max_rad*(ratio-0.5d0*ratio*ratio)/center_dist_corr
            olap_2 = olap_actual - olap_1
            ! get effective minimum conduction distance for each particle
            ratio = (olap_actual + des_min_cond_dist) / max_rad
            s_1 = (max_rad*max_rad*(ratio-0.5d0*ratio*ratio)/&
            &(center_dist_corr-des_min_cond_dist)) - olap_1
            s_2 = des_min_cond_dist - s_1
            ! get effective lens radius for particle-contact plane
            rlens_1 = sqrt(rd_out*rd_out+(min_rad-olap_1)**2.0)
            rlens_2 = sqrt(rd_out*rd_out+(max_rad-olap_2)**2.0)
            

! GET GAS THERMAL CONDUCTIVITY (default calculation is done in calc_k_g and is for air)
            if(k_g0.eq.undefined)then
                  ! Compute gas conductivity as is done in calc_k_g
                  ! But use average of particle and wall temperature to be gas temperature
               k_gas = 6.02d-5*sqrt(0.5d0*(des_t_s(i)+des_t_s(j))/300.d0) ! cal/(s.cm.K)
                  ! 1 cal = 4.183925D0 J
               IF (units == 'SI') k_gas = 418.3925d0*k_gas !J/s.m.K
            else
               k_gas=k_g0
            endif
   
            h_1 = eval_h_pfp(rlens_1, s_1, olap_1,min_rad)*min_rad*k_gas
            h_2 = eval_h_pfp(rlens_2, s_2, olap_2,max_rad)*max_rad*k_gas
            IF(h_1.eq.zero.OR.h_2.eq.zero)THEN
               h = zero
            ELSE
               h = h_1*h_2/(h_1+h_2)
            ENDIF
     
            q_pfp = h *deltatp
! Particle-fluid-particle is analytically computed using conductance for each
! particle to the contact plane.  The effective lens radius and minimum conduction
! distance are first calculated for each particle-fluid-contact_plane conduction. 
            
         ELSE
            q_pfp = zero
         ENDIF
         des_conduction = q_pp + q_pfp
               
            
 !!        ENDIF ! NEW WAY
!-------------------------------------------------------------------------------------
! OLD WAY
!-------------------------------------------------------------------------------------
!!         IF(.FALSE.)THEN
! Calculate the particle-fluid-particle conduction
! REF: Rong and Horio, 1999 (MODIFIED)
!---------------------------------------------------------------------//
! Calculate the radius of the fluid lens surrounding the larger particle
! Default FLPC = 0.2
!!            LENS_RAD = MAX_RAD * (1.0D0 + FLPC)

! Calculate the outer radial distance of the region for particle-fluid-
! particle heat conduction.
!!            RD_OUT = RADIUS( LENS_RAD, MIN_RAD)

! If the value returned is less than zero, then the fluid lens
! surrounding the larger particle does not intersect with the surface
! of the smaller particle. In this case, particle-fluild-particle
! conduction does not occur.
!!            Q_pfp = 0.0
!!            IF(RD_OUT .GT. ZERO)THEN
! Calculate the distance from the line connecting the particles' centers
! to the point of contact between the two particles. This value is
! zero if the particles are not touching and is the radius of the
! shared contact area otherwise.
!!               RD_IN = ZERO
!!               IF(CENTER_DIST_CORR < (MAX_RAD + MIN_RAD) ) &
!!               RD_IN = RADIUS(MAX_RAD, MIN_RAD)
! Calculate the rate of heat transfer between the particles through the
! fluid using adaptive Simpson's rule to manage the integral.
!!               Q_pfp = K_g(iIJK) * DeltaTp * &
!!               ADPT_SIMPSON(RD_IN,RD_OUT)

!!            ENDIF               ! RD_OUT > ZERO
!!         ENDIF ! OLD WAY

         ! Return total heat flux
!!         DES_CONDUCTION = Q_pp + Q_pfp
!! ------------------------------------------------------------------//
!  END OLD WAY
!! ------------------------------------------------------------------//

      RETURN

      CONTAINS

!......................................................................!
! Subroutine Procedure: FUNCTION RADIUS                                !
!                                                                      !
! Propose: Calculate the center line radius between two particles.     !
! Depending on what is given as input, this function calculates:       !
!   1) radius of contact area between two particles                    !
!   2) radius delineating particle-fluid-particle region.              !
!......................................................................!
      DOUBLE PRECISION FUNCTION radius(R1, R2)

      IMPLICIT NONE

! Radius values
      DOUBLE PRECISION, INTENT(IN) :: R1, R2
! Distance between particle centers
      DOUBLE PRECISION BETA
! Generic value
      DOUBLE PRECISION VALUE

! Calculate
      VALUE = (r1**2 - r2**2 + center_dist_corr**2)/(2.0d0 * r1 * center_dist_corr)
! Check to ensure that VALUE is less than or equal to one. If VALUE is
! greater than one, the triangle inequality has been violated. Therefore
! there is no intersection between the fluid lens surrounding the larger
! particle and the surface of the smaller particle.
! Thus, there is no particle-fluid-particle heat transfer.
      IF( VALUE .GT. 1.0d0) THEN
         radius = -1.0d0
      ELSE
! Calculate beta (Law of cosines)
         beta = acos( VALUE )
! Calculate the radius
         radius = r1 * sin(beta)
      ENDIF

      RETURN
      END FUNCTION radius

!......................................................................!
! Subroutine Procedure: FUNCTION K_eff                                 !
!                                                                      !
! Propose: Calculate the effective thermal conductivity of two         !
! particles with conductivities K1 and K2                              !
!......................................................................!
      DOUBLE PRECISION FUNCTION k_eff(K1, K2)

      IMPLICIT NONE

! Thermal conductivities
      DOUBLE PRECISION, INTENT(IN) :: K1, K2

      k_eff = 2.0d0 * (k1*k2)/(k1 + k2)

      RETURN
      END FUNCTION k_eff

!``````````````````````````````````````````````````````````````````````!
! Function: F                                                          !
!                                                                      !
! Purpose: This function defines the region of particle-fluid-particle !
!          heat transfer. Note that this function develops a           !
!          singularity at the boundary of the contact region when two  !
!          particles are touching.  This singularity is resolved by    !
!          making the assumption that the surfaces of two particles    !
!          never directly touch (c.f. Rong and Horio, 1999). By        !
!          default, it is assumed that the particles are separated by  !
!          a minimum length of 4.0x10^(-10) meters. This value may be  !
!          modified by specifying a different value in the mfix.dat    !
!          file under the variable name DES_MIN_COND_DIST.             !
!                                                                      !
!          Note, the closer this value is to zero, the harder it will  !
!          be for the quadrature routine to converge.                  !
!......................................................................!
      DOUBLE PRECISION FUNCTION f(R)
        USE utilities
        USE des_thermo
      IMPLICIT NONE

      DOUBLE PRECISION, INTENT(IN) :: R

      f = (2.0d0*pi*r)/max((center_dist_corr - sqrt(max_rad**2-r**2) - &
         sqrt(min_rad**2-r**2)), des_min_cond_dist)

      IF( mfix_isnan(f))print*,'F IS NAN'

      END FUNCTION f

!``````````````````````````````````````````````````````````````````````!
! Function: ADPT_SIMPSON                                               !
!                                                                      !
! Purpose: This function applies adaptive Simpson's Rule to solve the  !
!          definite integral of the function F.                        !
!                                                                      !
!         *NOTE: This route will return the result of the method even  !
!          if convergence has not been achieved. If this occurs, a     !
!          message will be written to the log file to notify the user. !
!......................................................................!
      DOUBLE PRECISION FUNCTION adpt_simpson(A, B)

      IMPLICIT NONE

! Bounds of integration
      DOUBLE PRECISION, INTENT(IN) :: A, B
! Maximum recursive depth
      INTEGER, PARAMETER :: Kmax = 25
! Current depth of recursion
      INTEGER :: K
! Array storage for recursion
      DOUBLE PRECISION :: V(kmax,6)
! Step size
      DOUBLE PRECISION :: H
! Simpson's Rule evaluated on the left and right intervals
      DOUBLE PRECISION :: lS, rS
! Value of the function F evaluate at the midpoints of the left and
! right intervals
      DOUBLE PRECISION :: rF, lF
! Local error bound
      DOUBLE PRECISION :: Err_BND
! Convergence bound
      DOUBLE PRECISION :: EPS = 10.0**(-2)

! Error indicating that an error has been specified to the user.
      LOGICAL, SAVE :: ADPT_SIMPSON_ERR = .false.

! Initialize variables
      v(:,:) = 0.0d0   !
      adpt_simpson = 0.0d0   ! Integral value
      h = (b-a)/2.0d0 ! Dynamic interval length
! Calculate Simpson's Rule over the interval [A,B]
      ls = (f(a) + 4.0d0*f((a+b)/2.0d0) + f(b))*(h/3.0d0)
! Initialize the storage vector for the first pass
      v(1,:) = (/ a, h, f(a), f((a+b)/2.0d0), f(b), ls/)
      k = 1
! Enter recursion calculations
      DO
! Establish the new interval length
         h = v(k,2)/2.0d0
! Evaluate the function on the left interval
         lf = f(v(k,1) + h)
! Calculate Simpson's Rule over the left interval
         ls = (v(k,3) + 4.0d0*lf + v(k,4))*(h/3.0d0)
! Evaluate the function on the right interval
         rf = f(v(k,1) + 3.0d0*h)
! Calculate Simpson's Rule over the right interval
         rs = (v(k,4) + 4.0d0*rf + v(k,5))*(h/3.0d0)
! Evaluate the error bound
         err_bnd = (30.0d0*eps*h)/(b-a)
! Check to see if conversion on the interval has been met
         IF( abs(ls + rs - v(k,6)) .LT. err_bnd)THEN
! Update the integral value
            adpt_simpson = adpt_simpson + ls + rs + &
               (1.0d0/15.0d0)*(ls + rs - v(k,6))
! Decrement the recursion counter
            k = k - 1
! If K=0, then integration has been successfully completed over the
! entire interval [A,B].
            IF( k == 0) RETURN
         ELSEIF( (k .GE. kmax) .OR. &
            (h == (b-a)*(1.0d0/2.0d0)**(kmax+3))) THEN
! Flag that the method did not converge.
            IF(.NOT.adpt_simpson_err)THEN
               WRITE(*,1000)
               WRITE(unit_log,1000)
               adpt_simpson_err = .true.
            ENDIF
! Update the integral value
            adpt_simpson = adpt_simpson + ls + rs + &
               (1.0d0/15.0d0)*(ls + rs - v(k,6))
! Decrement the recursion counter
            k = k - 1
         ELSE
! Refine the subintervals through recursive splitting of the intervals
            v(k+1,:) = (/v(k,1) + 2.0d0*h, h, v(k,4), rf, v(k,5), rs/)
            v(k,:) = (/v(k,1), h, v(k,3), lf, v(k,4), ls/)
! Increment level counter
            k = k+1
         ENDIF
      ENDDO

 1000 FORMAT(/1x,70('*'),/' From: DES_COND_EQ',/, ' Message: ',        &
         'Integration of the particle-fluid-particle equation did ',   &
         'not',/' converge! No definite bound can be placed on the ',  &
         'error.',/' Future convergence messages will be suppressed!', &
         /1x,70('*'))

      END FUNCTION adpt_simpson

     DOUBLE PRECISION FUNCTION eval_h_pfp(RLENS_dim,S,OLAP_dim,RP)
      USE constant
      USE param1
     
      IMPLICIT NONE
      ! Note: Function inputs dimensional quantities
      DOUBLE PRECISION, intent(in) :: RLENS_dim, S, OLAP_dim, RP
      ! BELOW VARIABLES ARE NONDIMENSIONALIZED BY RP
      DOUBLE PRECISION :: RLENS, OLAP, KN
      DOUBLE PRECISION :: TERM1,TERM2,TERM3
      DOUBLE PRECISION :: Rout,Rkn
      DOUBLE PRECISION, PARAMETER :: TWO = 2.0d0
      
      rlens = rlens_dim/rp
      kn = s/rp
      olap = olap_dim/rp

      IF(-olap.ge.(rlens-one))THEN
         rout = zero
      ELSEIF(olap.le.two)THEN
         rout = sqrt(rlens**2-(one-olap)**2)
      ELSE
         WRITE(*,*)'ERROR: Extremeley excessive overlap (Olap > Diam)'
         WRITE(*,*)'OLAP = ',olap_dim,olap, rp
         WRITE(*,*)'RLENS_dim',rlens_dim, rlens
         write(*,*)'S, Kn', s,kn
         stop
      ENDIF
      rout=min(rout,one)
      
      IF(olap.ge.zero)THEN
     !     particle in contact(below code verified)
         term1 = pi*((one-olap)**2-(one-olap-kn)**2)/kn
         term2 = two*pi*(sqrt(one-rout**2)-(one-olap-kn))
         term3 = two*pi*(one-olap)*log((one-olap-sqrt(one-rout**2))/kn)
         eval_h_pfp = term1+term2+term3
      ELSE
         IF(-olap.ge.kn)THEN
            rkn = zero
         ELSE
            rkn=sqrt(one-(one-olap-kn)**2)
         ENDIF
    
         term1 = (rkn**2/(two*kn))+sqrt(one-rout**2)-sqrt(one-rkn**2)
         term2 = (one-olap)*log((one-olap-sqrt(one-rout**2))/(one-olap-sqrt(one-rkn**2)))
         eval_h_pfp = two*pi*(term1+term2)
         
      ENDIF
      RETURN
      END FUNCTION eval_h_pfp

      END FUNCTION des_conduction

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: DES_CONDUCTION                                         !
!                                                                      !
!  Purpose: Compute conductive heat transfer with wall                 !
!                                                                      !
!  Author: A.Morris                                  Date: 07-Jan-16   !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      FUNCTION des_conduction_wall(OLAP, K_sol, K_wall, K_gas, TWall, &
      tpart, rpart, rlens, m)
      USE param1, only: zero
      USE des_thermo, only: des_min_cond_dist

      DOUBLE PRECISION :: DES_CONDUCTION_WALL 
      DOUBLE PRECISION, intent(in) :: OLAP, K_sol, K_wall, K_gas, TWall&
      &                               ,TPart,RPart, RLens
      DOUBLE PRECISION :: Rin, Rout
      DOUBLE PRECISION :: OLAP_ACTUAL, TAUC_CORRECTION
      DOUBLE PRECISION :: KEff
      DOUBLE PRECISION :: Q_pw, Q_pfw
      INTEGER :: M ! Solids phase index

      ! Initialize variables
      q_pw = zero
      q_pfw= zero
      olap_actual = olap

      ! Check to see if particle is in contact (P-W) conduction
      IF(olap.gt.zero)THEN
      ! Compute effective solids conductivity
         keff = 2.0d0*k_sol*k_wall/(k_sol+k_wall)
      ! Correct for overlap using "area" correction terms
         IF(do_area_correction)THEN
            olap_actual = correct_olap(olap,m,-1) ! (-1 indicates wall)
         ENDIF
      ! Compute inner Rad
         rin = sqrt(2.0d0*olap_actual*rpart-olap_actual*olap_actual)
      ! Compute time correction term (commented out for now)
        ! CALL CALC_TIME_CORRECTION(TAUC_CORRECTION,M,-1)
      ! Compute heat transfer (particle-wall)
         q_pw = 2.0d0*rin*keff*(twall-tpart)
      ENDIF
      ! Compute heat transfer (particle-fluid-wall)
      q_pfw=eval_h_pfw(rlens,des_min_cond_dist, olap_actual, rpart) * &
      &     k_gas*rpart*(twall-tpart)

      des_conduction_wall=q_pw+q_pfw
     
      RETURN 
      END FUNCTION des_conduction_wall


!     THIS FUNCTION COMPUTES THE OVERLAP THAT WOULD OCCUR (for a given force)
!     IF ACTUAL MATERIAL PROPERTIES WERE USED. 
      FUNCTION correct_olap(OLAP,M,L)
      use discretelement
      IMPLICIT NONE
      DOUBLE PRECISION :: CORRECT_OLAP
      DOUBLE PRECISION, INTENT (IN) :: OLAP
      INTEGER, INTENT (IN) :: M, L
      DOUBLE PRECISION :: KN_ACTUAL, KN_SIM
      ! L=-1 corresponds to wall
      IF(l.eq.-1)THEN
         ! WALL CONTACT
         IF(des_coll_model_enum .EQ. hertzian)THEN
            correct_olap = (hert_kwn(m)/hert_kwn_actual(m))**(2.0d0/3.0d0)*olap
         ELSE
            correct_olap = (kn_w*olap/hert_kwn_actual(m))**(2.0d0/3.0d0)
         ENDIF
      ELSE
         ! PARTICLE-PARTICLE
         IF(des_coll_model_enum .EQ. hertzian)THEN
            correct_olap = (hert_kn(m,l)/hert_kn_actual(m,l))**(2.0d0/3.0d0)*olap
         ELSE
            correct_olap = (kn*olap/hert_kn_actual(m,l))**(2.0d0/3.0d0)
         ENDIF
      ENDIF
      RETURN
      END FUNCTION correct_olap


      DOUBLE PRECISION FUNCTION eval_h_pfw(RLENS_dim,S,OLAP_dim,RP)
      USE constant
      USE param1
     
      IMPLICIT NONE
      ! Note: Function inputs dimensional quantities
      DOUBLE PRECISION, intent(in) :: RLENS_dim, S, OLAP_dim, RP
      ! BELOW VARIABLES ARE NONDIMENSIONALIZED BY RP
      DOUBLE PRECISION :: RLENS, OLAP, KN
      DOUBLE PRECISION :: TERM1,TERM2,TERM3
      DOUBLE PRECISION :: Rout,Rkn
      DOUBLE PRECISION, PARAMETER :: TWO = 2.0d0
      
      rlens = rlens_dim/rp
      kn = s/rp
      olap = olap_dim/rp

      IF(-olap.ge.(rlens-one))THEN
         rout = zero
      ELSEIF(olap.le.two)THEN
         rout = sqrt(rlens**2-(one-olap)**2)
      ELSE
         WRITE(*,*)'ERROR: Extremeley excessive overlap (Olap > Diam)'
         WRITE(*,*)'OLAP = ',olap_dim,olap, rp
         WRITE(*,*)'RLENS_dim',rlens_dim, rlens
         write(*,*)'S, Kn', s,kn
         stop
      ENDIF
      rout=min(rout,one)
      
      IF(olap.ge.zero)THEN
     !     particle in contact(below code verified)
         term1 = pi*((one-olap)**2-(one-olap-kn)**2)/kn
         term2 = two*pi*(sqrt(one-rout**2)-(one-olap-kn))
         term3 = two*pi*(one-olap)*log((one-olap-sqrt(one-rout**2))/kn)
         eval_h_pfw = term1+term2+term3
      ELSE
         IF(-olap.ge.kn)THEN
            rkn = zero
         ELSE
            rkn=sqrt(one-(one-olap-kn)**2)
         ENDIF
    
         term1 = (rkn**2/(two*kn))+sqrt(one-rout**2)-sqrt(one-rkn**2)
         term2 = (one-olap)*log((one-olap-sqrt(one-rout**2))/(one-olap-sqrt(one-rkn**2)))
         eval_h_pfw = two*pi*(term1+term2)
         
      ENDIF
      RETURN
      END FUNCTION eval_h_pfw


      SUBROUTINE calc_time_correction(time_corr, phaseI, phaseJ)
      use discretelement, only: tau_c_base_actual, tauw_c_base_actual
      use discretelement, only: tau_c_base_sim, tauw_c_base_sim
      use discretelement, only: des_coll_model_enum, hertzian
      IMPLICIT NONE
      INTEGER, intent (in) :: phaseI, phaseJ
      DOUBLE PRECISION :: time_corr, tau_actual, tau_sim
      DOUBLE PRECISION :: vimp ! impact veloctiy
      INTEGER :: M, N
      vimp = 1.0d0
      time_corr = 1.0d0
      if (phasej == -1) then
         ! Wall contact
         m = phasei
         IF (des_coll_model_enum .EQ. hertzian) THEN
            time_corr = tauw_c_base_actual(m) / tauw_c_base_sim(m)
         ELSE
            vimp = 1.0d0
            if (vimp .le. 0.0d0)then
               time_corr = 1.0d0
            else
               time_corr = tauw_c_base_actual(m)*vimp**(-0.2d0) / tauw_c_base_sim(m)
            endif
         ENDIF

      ELSE
         ! particle-particle contact
         if (phasei .le. phasej)then
            m = phasei
            n = phasej
         else
            m = phasej
            n = phasei
         endif
         
         IF (des_coll_model_enum .EQ. hertzian) THEN
            time_corr = tau_c_base_actual(m,n) / tau_c_base_sim(m,n)
         ELSE
            vimp = 1.0d0
            if (vimp .le. 0.0d0)then
               time_corr = 1.0d0
            else
               time_corr = tau_c_base_actual(m,n)*vimp**(-0.2d0) / tau_c_base_sim(m,n)
            endif
         ENDIF
      ENDIF
      time_corr = time_corr **(2.0d0/3.0d0)
      ! TEMPORARY 
      time_corr = 1.0d0
      RETURN
      END SUBROUTINE calc_time_correction


      
  END MODULE des_thermo_cond