File: N:\mfix\model\des\des_thermo_cond_mod.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module name: DES_CONDUCTION                                         !
4     !                                                                      !
5     !  Purpose:                                                            !
6     !                                                                      !
7     !  Author: J.Musser                                   Date: 25-Jun-10  !
8     !                                                                      !
9     !  Comments:                                                           !
10     !                                                                      !
11     !  REF: Batchelor and O'Brien, "Thermal or electrical conduction       !
12     !       through a granular material," Proceedings of the Royal Society !
13     !       of London. Series A, Mathematical and Physical Sciences,       !
14     !       Vol. 355, no 1682, pp. 313-333, July 1977.                     !
15     !                                                                      !
16     !  REF: Rong and Horio, "DEM simulation of char combustion in a        !
17     !       fluidized bed," in Second International Conference on CFD in   !
18     !       the materials and Process Industries, Melbourne, 1999, pp.     !
19     !       65-70.                                                         !
20     !                                                                      !
21     !  REF: Xavier and Davidson, "Heat transfer to surfaces immersed in    !
22     !       fluidised beds, particularly tub arrays," in Fluidization:     !
23     !       Proceedings of the Second Engineering Foundation Conference,   !
24     !       Cambridge, England, 2-6 April, 1978, pp. 333-338.              !
25     !                                                                      !
26     !  REF: Zhou, Yu, and Zulli, "Particle scale study of heat transfer in !
27     !       packed and bubbling fluidized beds," AIChE Journal, Vol. 55,   !
28     !       no 4, pp 868-884, 2009.                                        !
29     !                                                                      !
30     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
31           MODULE DES_THERMO_COND
32      
33     
34           IMPLICIT NONE
35           DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: DES_Qw_cond
36           LOGICAL :: DO_AREA_CORRECTION
37           CONTAINS
38     
39     
40     
41           FUNCTION DES_CONDUCTION(I, J, CENTER_DIST, iM, iIJK)
42     
43           use constant
44           use des_thermo
45           use discretelement
46           use funits
47           use param1, only: zero, UNDEFINED
48           use physprop
49           use run, only: UNITS
50           IMPLICIT NONE
51     
52     ! Passed variables
53     !---------------------------------------------------------------------//
54     ! Index of particle being looped over
55           INTEGER, INTENT(IN) :: I,J
56     ! Index of fluid cell containing particle I
57           INTEGER, INTENT(IN) :: iIJK
58     ! Solid phase indices for the given particles
59           INTEGER, INTENT(IN) :: iM
60     ! Distance between the centers of particle I and particle J (magnitude)
61           DOUBLE PRECISION, INTENT(IN) :: CENTER_DIST
62           DOUBLE PRECISION :: DES_CONDUCTION
63     
64     ! Local variables
65     !---------------------------------------------------------------------//
66     ! Solid phase indices for the given particles
67           INTEGER jM
68     ! Radius of smaller particle
69           DOUBLE PRECISION MIN_RAD
70     ! Radius of larger particle
71           DOUBLE PRECISION MAX_RAD
72     ! Rate of particle-particle conduction
73           DOUBLE PRECISION Q_pp
74     ! Rate of particle-fluid-particle conduction
75           DOUBLE PRECISION Q_pfp
76     ! Outer radius of region delineating particle-fluid-particle conduction
77           DOUBLE PRECISION RD_OUT
78     ! Inner radius of region delineating particle-fluid-particle conduction
79           DOUBLE PRECISION RD_IN
80     ! The radius of the fluid lens containing the larger particle
81           DOUBLE PRECISION LENS_RAD
82     ! Temperature difference between two particles
83           DOUBLE PRECISION DeltaTp
84     ! Effective thermal conductivity
85           DOUBLE PRECISION lK_eff
86     ! Radius of contact area
87           DOUBLE PRECISION lRadius
88     ! Particle overlap (simulated and corrected)
89           DOUBLE PRECISION OLAP_Sim, OLAP_actual
90     ! Corrected center distance
91           DOUBLE PRECISION CENTER_DIST_CORR
92     ! Gas thermal conductivity
93           DOUBLE PRECISION :: K_gas
94     ! Overlap distance between each particle and contact plane
95           DOUBLE PRECISION :: OLAP_1, OLAP_2
96     ! Effective Lens for particles 1 and 2 (req'd for analytic calculation)
97           DOUBLE PRECISION :: RLENS_1, RLENS_2
98     ! Effective minimum conduction distance for particles 1 and 2
99           DOUBLE PRECISION :: S_1, S_2
100     ! Effective thermal conductance for particles 1 and 2
101           DOUBLE PRECISION :: H_1, H_2
102     ! Analytic thermal conductance between particles
103           DOUBLE PRECISION :: H
104     ! Dummy variable for calculations 
105           DOUBLE PRECISION :: RATIO
106     ! Functions
107     !---------------------------------------------------------------------//
108      !     DOUBLE PRECISION :: EVAL_H_PFP
109     
110     ! Identify the solid phases of the neighbor particle
111              jM = PIJK(J,5)
112     
113     ! Determine the radius of the larger and smaller particle
114              MIN_RAD = MIN(DES_RADIUS(I), DES_RADIUS(J))
115              MAX_RAD = MAX(DES_RADIUS(I), DES_RADIUS(J))
116     
117              DeltaTp = DES_T_s(J) - DES_T_s(I)
118     
119     ! Calculate the particle-particle conduction
120     ! REF: Batchelor and O'Brien, 1977 (MODIFIED)
121     !---------------------------------------------------------------------//
122              Q_pp = 0.0
123     
124     ! Initialize corrected center-center distance
125              CENTER_DIST_CORR = CENTER_DIST
126              IF(CENTER_DIST < (MAX_RAD + MIN_RAD)) THEN
127     ! Correct particle overlap to account for artificial softening
128                 OLAP_Sim = MAX_RAD + MIN_RAD - CENTER_DIST
129                 IF(DO_AREA_CORRECTION)THEN
130                    IF(Im.ge.Jm)THEN
131                       OLAP_Actual = CORRECT_OLAP(OLAP_Sim,Jm,Im)
132                    ELSE
133                       OLAP_Actual = CORRECT_OLAP(OLAP_Sim,Im,Jm)
134                    ENDIF
135                    CENTER_DIST_CORR = MAX_RAD + MIN_RAD - OLAP_Actual
136                 ENDIF
137     ! Effective thermal conductivity
138                 lK_eff = K_eff(K_s0(iM),K_s0(jM))
139     ! Effective contact area's radius
140                 lRadius = RADIUS(MAX_RAD, MIN_RAD)
141     ! Compute time-correction term (not done yet)
142                 ! CALL CALC_TIME_CORRECTION (.... )
143     ! Inter-particle heat transfer
144                 Q_pp = 2.0d0 * lK_eff * lRadius * DeltaTp
145     
146     ! Assign the inter-particle heat transfer to both particles.
147              ENDIF
148     
149     !!         IF(.TRUE.)THEN ! NEW WAY
150     !---------------------------------------------------------------------//
151     ! Calculate the particle-fluid-particle conduction
152     ! REF: Rong and Horio, 1999 (MODIFIED)
153     !---------------------------------------------------------------------//
154              LENS_RAD = MAX_RAD * (1.0D0 + FLPC)
155              OLAP_Actual = MAX_RAD + MIN_RAD - CENTER_DIST_CORR
156     ! check to see if particles are within lens thickness from eachother
157     ! a negative overlap indicates particles are separated 
158              IF(OLAP_actual > (-MAX_RAD*FLPC))THEN
159                 RD_OUT = RADIUS(LENS_RAD, MIN_RAD)
160                 ! get overlap between each particle and contact plane
161                 RATIO = OLAP_actual / MAX_RAD
162                 OLAP_1 = MAX_RAD*MAX_RAD*(RATIO-0.5D0*RATIO*RATIO)/CENTER_DIST_CORR
163                 OLAP_2 = OLAP_actual - OLAP_1
164                 ! get effective minimum conduction distance for each particle
165                 RATIO = (OLAP_actual + DES_MIN_COND_DIST) / MAX_RAD
166                 S_1 = (MAX_RAD*MAX_RAD*(RATIO-0.5D0*RATIO*RATIO)/&
167                 &(CENTER_DIST_CORR-DES_MIN_COND_DIST)) - OLAP_1
168                 S_2 = DES_MIN_COND_DIST - S_1
169                 ! get effective lens radius for particle-contact plane
170                 RLENS_1 = sqrt(RD_OUT*RD_OUT+(MIN_RAD-OLAP_1)**2.0)
171                 RLENS_2 = sqrt(RD_OUT*RD_OUT+(MAX_RAD-OLAP_2)**2.0)
172                 
173     
174     ! GET GAS THERMAL CONDUCTIVITY (default calculation is done in calc_k_g and is for air)
175                 if(k_g0.eq.UNDEFINED)then
176                       ! Compute gas conductivity as is done in calc_k_g
177                       ! But use average of particle and wall temperature to be gas temperature
178                    K_Gas = 6.02D-5*SQRT(0.5d0*(DES_T_S(I)+DES_T_S(J))/300.D0) ! cal/(s.cm.K)
179                       ! 1 cal = 4.183925D0 J
180                    IF (UNITS == 'SI') K_Gas = 418.3925D0*K_Gas !J/s.m.K
181                 else
182                    K_Gas=k_g0
183                 endif
184        
185                 H_1 = EVAL_H_PFP(RLENS_1, S_1, OLAP_1,MIN_RAD)*MIN_RAD*k_gas
186                 H_2 = EVAL_H_PFP(RLENS_2, S_2, OLAP_2,MAX_RAD)*MAX_RAD*k_gas
187                 IF(H_1.eq.ZERO.OR.H_2.eq.ZERO)THEN
188                    H = ZERO
189                 ELSE
190                    H = H_1*H_2/(H_1+H_2)
191                 ENDIF
192          
193                 Q_pfp = H *DeltaTp
194     ! Particle-fluid-particle is analytically computed using conductance for each
195     ! particle to the contact plane.  The effective lens radius and minimum conduction
196     ! distance are first calculated for each particle-fluid-contact_plane conduction. 
197                 
198              ELSE
199                 Q_pfp = ZERO
200              ENDIF
201              DES_CONDUCTION = Q_pp + Q_pfp
202                    
203                 
204      !!        ENDIF ! NEW WAY
205     !-------------------------------------------------------------------------------------
206     ! OLD WAY
207     !-------------------------------------------------------------------------------------
208     !!         IF(.FALSE.)THEN
209     ! Calculate the particle-fluid-particle conduction
210     ! REF: Rong and Horio, 1999 (MODIFIED)
211     !---------------------------------------------------------------------//
212     ! Calculate the radius of the fluid lens surrounding the larger particle
213     ! Default FLPC = 0.2
214     !!            LENS_RAD = MAX_RAD * (1.0D0 + FLPC)
215     
216     ! Calculate the outer radial distance of the region for particle-fluid-
217     ! particle heat conduction.
218     !!            RD_OUT = RADIUS( LENS_RAD, MIN_RAD)
219     
220     ! If the value returned is less than zero, then the fluid lens
221     ! surrounding the larger particle does not intersect with the surface
222     ! of the smaller particle. In this case, particle-fluild-particle
223     ! conduction does not occur.
224     !!            Q_pfp = 0.0
225     !!            IF(RD_OUT .GT. ZERO)THEN
226     ! Calculate the distance from the line connecting the particles' centers
227     ! to the point of contact between the two particles. This value is
228     ! zero if the particles are not touching and is the radius of the
229     ! shared contact area otherwise.
230     !!               RD_IN = ZERO
231     !!               IF(CENTER_DIST_CORR < (MAX_RAD + MIN_RAD) ) &
232     !!               RD_IN = RADIUS(MAX_RAD, MIN_RAD)
233     ! Calculate the rate of heat transfer between the particles through the
234     ! fluid using adaptive Simpson's rule to manage the integral.
235     !!               Q_pfp = K_g(iIJK) * DeltaTp * &
236     !!               ADPT_SIMPSON(RD_IN,RD_OUT)
237     
238     !!            ENDIF               ! RD_OUT > ZERO
239     !!         ENDIF ! OLD WAY
240     
241              ! Return total heat flux
242     !!         DES_CONDUCTION = Q_pp + Q_pfp
243     !! ------------------------------------------------------------------//
244     !  END OLD WAY
245     !! ------------------------------------------------------------------//
246     
247           RETURN
248     
249           CONTAINS
250     
251     !......................................................................!
252     ! Subroutine Procedure: FUNCTION RADIUS                                !
253     !                                                                      !
254     ! Propose: Calculate the center line radius between two particles.     !
255     ! Depending on what is given as input, this function calculates:       !
256     !   1) radius of contact area between two particles                    !
257     !   2) radius delineating particle-fluid-particle region.              !
258     !......................................................................!
259           DOUBLE PRECISION FUNCTION RADIUS(R1, R2)
260     
261           IMPLICIT NONE
262     
263     ! Radius values
264           DOUBLE PRECISION, INTENT(IN) :: R1, R2
265     ! Distance between particle centers
266           DOUBLE PRECISION BETA
267     ! Generic value
268           DOUBLE PRECISION VALUE
269     
270     ! Calculate
271           VALUE = (R1**2 - R2**2 + CENTER_DIST_CORR**2)/(2.0d0 * R1 * CENTER_DIST_CORR)
272     ! Check to ensure that VALUE is less than or equal to one. If VALUE is
273     ! greater than one, the triangle inequality has been violated. Therefore
274     ! there is no intersection between the fluid lens surrounding the larger
275     ! particle and the surface of the smaller particle.
276     ! Thus, there is no particle-fluid-particle heat transfer.
277           IF( VALUE .GT. 1.0d0) THEN
278              RADIUS = -1.0d0
279           ELSE
280     ! Calculate beta (Law of cosines)
281              BETA = ACOS( VALUE )
282     ! Calculate the radius
283              RADIUS = R1 * SIN(BETA)
284           ENDIF
285     
286           RETURN
287           END FUNCTION RADIUS
288     
289     !......................................................................!
290     ! Subroutine Procedure: FUNCTION K_eff                                 !
291     !                                                                      !
292     ! Propose: Calculate the effective thermal conductivity of two         !
293     ! particles with conductivities K1 and K2                              !
294     !......................................................................!
295           DOUBLE PRECISION FUNCTION K_eff(K1, K2)
296     
297           IMPLICIT NONE
298     
299     ! Thermal conductivities
300           DOUBLE PRECISION, INTENT(IN) :: K1, K2
301     
302           K_eff = 2.0d0 * (K1*K2)/(K1 + K2)
303     
304           RETURN
305           END FUNCTION K_eff
306     
307     !``````````````````````````````````````````````````````````````````````!
308     ! Function: F                                                          !
309     !                                                                      !
310     ! Purpose: This function defines the region of particle-fluid-particle !
311     !          heat transfer. Note that this function develops a           !
312     !          singularity at the boundary of the contact region when two  !
313     !          particles are touching.  This singularity is resolved by    !
314     !          making the assumption that the surfaces of two particles    !
315     !          never directly touch (c.f. Rong and Horio, 1999). By        !
316     !          default, it is assumed that the particles are separated by  !
317     !          a minimum length of 4.0x10^(-10) meters. This value may be  !
318     !          modified by specifying a different value in the mfix.dat    !
319     !          file under the variable name DES_MIN_COND_DIST.             !
320     !                                                                      !
321     !          Note, the closer this value is to zero, the harder it will  !
322     !          be for the quadrature routine to converge.                  !
323     !......................................................................!
324           DOUBLE PRECISION FUNCTION F(R)
325             USE utilities
326             USE des_thermo
327           IMPLICIT NONE
328     
329           DOUBLE PRECISION, INTENT(IN) :: R
330     
331           F = (2.0d0*Pi*R)/MAX((CENTER_DIST_CORR - SQRT(MAX_RAD**2-R**2) - &
332              SQRT(MIN_RAD**2-R**2)), DES_MIN_COND_DIST)
333     
334           IF( mfix_isnan(F))PRINT*,'F IS NAN'
335     
336           END FUNCTION F
337     
338     !``````````````````````````````````````````````````````````````````````!
339     ! Function: ADPT_SIMPSON                                               !
340     !                                                                      !
341     ! Purpose: This function applies adaptive Simpson's Rule to solve the  !
342     !          definite integral of the function F.                        !
343     !                                                                      !
344     !         *NOTE: This route will return the result of the method even  !
345     !          if convergence has not been achieved. If this occurs, a     !
346     !          message will be written to the log file to notify the user. !
347     !......................................................................!
348           DOUBLE PRECISION FUNCTION ADPT_SIMPSON(A, B)
349     
350           IMPLICIT NONE
351     
352     ! Bounds of integration
353           DOUBLE PRECISION, INTENT(IN) :: A, B
354     ! Maximum recursive depth
355           INTEGER, PARAMETER :: Kmax = 25
356     ! Current depth of recursion
357           INTEGER :: K
358     ! Array storage for recursion
359           DOUBLE PRECISION :: V(Kmax,6)
360     ! Step size
361           DOUBLE PRECISION :: H
362     ! Simpson's Rule evaluated on the left and right intervals
363           DOUBLE PRECISION :: lS, rS
364     ! Value of the function F evaluate at the midpoints of the left and
365     ! right intervals
366           DOUBLE PRECISION :: rF, lF
367     ! Local error bound
368           DOUBLE PRECISION :: Err_BND
369     ! Convergence bound
370           DOUBLE PRECISION :: EPS = 10.0**(-2)
371     
372     ! Error indicating that an error has been specified to the user.
373           LOGICAL, SAVE :: ADPT_SIMPSON_ERR = .FALSE.
374     
375     ! Initialize variables
376           V(:,:) = 0.0d0   !
377           ADPT_SIMPSON = 0.0d0   ! Integral value
378           H = (B-A)/2.0d0 ! Dynamic interval length
379     ! Calculate Simpson's Rule over the interval [A,B]
380           lS = (F(A) + 4.0d0*F((A+B)/2.0d0) + F(B))*(H/3.0d0)
381     ! Initialize the storage vector for the first pass
382           V(1,:) = (/ A, H, F(A), F((A+B)/2.0d0), F(B), lS/)
383           K = 1
384     ! Enter recursion calculations
385           DO
386     ! Establish the new interval length
387              H = V(K,2)/2.0d0
388     ! Evaluate the function on the left interval
389              lF = F(V(K,1) + H)
390     ! Calculate Simpson's Rule over the left interval
391              lS = (V(K,3) + 4.0d0*lF + V(K,4))*(H/3.0d0)
392     ! Evaluate the function on the right interval
393              rF = F(V(K,1) + 3.0d0*H)
394     ! Calculate Simpson's Rule over the right interval
395              rS = (V(K,4) + 4.0d0*rF + V(K,5))*(H/3.0d0)
396     ! Evaluate the error bound
397              Err_BND = (30.0d0*EPS*H)/(B-A)
398     ! Check to see if conversion on the interval has been met
399              IF( ABS(lS + rS - V(K,6)) .LT. Err_BND)THEN
400     ! Update the integral value
401                 ADPT_SIMPSON = ADPT_SIMPSON + lS + rS + &
402                    (1.0d0/15.0d0)*(lS + rS - V(K,6))
403     ! Decrement the recursion counter
404                 K = K - 1
405     ! If K=0, then integration has been successfully completed over the
406     ! entire interval [A,B].
407                 IF( K == 0) RETURN
408              ELSEIF( (K .GE. Kmax) .OR. &
409                 (H == (B-A)*(1.0d0/2.0d0)**(Kmax+3))) THEN
410     ! Flag that the method did not converge.
411                 IF(.NOT.ADPT_SIMPSON_ERR)THEN
412                    WRITE(*,1000)
413                    WRITE(UNIT_LOG,1000)
414                    ADPT_SIMPSON_ERR = .TRUE.
415                 ENDIF
416     ! Update the integral value
417                 ADPT_SIMPSON = ADPT_SIMPSON + lS + rS + &
418                    (1.0d0/15.0d0)*(lS + rS - V(K,6))
419     ! Decrement the recursion counter
420                 K = K - 1
421              ELSE
422     ! Refine the subintervals through recursive splitting of the intervals
423                 V(K+1,:) = (/V(K,1) + 2.0d0*H, H, V(K,4), rF, V(K,5), rS/)
424                 V(K,:) = (/V(K,1), H, V(K,3), lF, V(K,4), lS/)
425     ! Increment level counter
426                 K = K+1
427              ENDIF
428           ENDDO
429     
430      1000 FORMAT(/1X,70('*'),/' From: DES_COND_EQ',/, ' Message: ',        &
431              'Integration of the particle-fluid-particle equation did ',   &
432              'not',/' converge! No definite bound can be placed on the ',  &
433              'error.',/' Future convergence messages will be suppressed!', &
434              /1X,70('*'))
435     
436           END FUNCTION ADPT_SIMPSON
437     
438          DOUBLE PRECISION FUNCTION EVAL_H_PFP(RLENS_dim,S,OLAP_dim,RP)
439           USE CONSTANT
440           USE PARAM1
441          
442           IMPLICIT NONE
443           ! Note: Function inputs dimensional quantities
444           DOUBLE PRECISION, intent(in) :: RLENS_dim, S, OLAP_dim, RP
445           ! BELOW VARIABLES ARE NONDIMENSIONALIZED BY RP
446           DOUBLE PRECISION :: RLENS, OLAP, KN
447           DOUBLE PRECISION :: TERM1,TERM2,TERM3
448           DOUBLE PRECISION :: Rout,Rkn
449           DOUBLE PRECISION, PARAMETER :: TWO = 2.0D0
450           
451           RLENS = RLENS_dim/RP
452           KN = S/RP
453           OLAP = OLAP_dim/RP
454     
455           IF(-OLAP.ge.(RLENS-ONE))THEN
456              Rout = ZERO
457           ELSEIF(OLAP.le.TWO)THEN
458              Rout = sqrt(RLENS**2-(ONE-OLAP)**2)
459           ELSE
460              WRITE(*,*)'ERROR: Extremeley excessive overlap (Olap > Diam)'
461              WRITE(*,*)'OLAP = ',OLAP_dim,OLAP, RP
462              WRITE(*,*)'RLENS_dim',RLENS_dim, RLENS
463              write(*,*)'S, Kn', S,KN
464              STOP
465           ENDIF
466           Rout=MIN(Rout,ONE)
467           
468           IF(OLAP.ge.ZERO)THEN
469          !     Particle in contact (below code verified)
470              TERM1 = PI*((ONE-OLAP)**2-(ONE-OLAP-KN)**2)/KN
471              TERM2 = TWO*PI*(sqrt(ONE-Rout**2)-(ONE-OLAP-KN))
472              TERM3 = TWO*PI*(ONE-OLAP)*log((ONE-OLAP-sqrt(ONE-Rout**2))/Kn)
473              EVAL_H_PFP = TERM1+TERM2+TERM3
474           ELSE
475              IF(-OLAP.ge.KN)THEN
476                 Rkn = ZERO
477              ELSE
478                 Rkn=sqrt(ONE-(ONE-OLAP-Kn)**2)
479              ENDIF
480         
481              TERM1 = (Rkn**2/(TWO*KN))+sqrt(ONE-Rout**2)-sqrt(ONE-Rkn**2)
482              TERM2 = (ONE-OLAP)*log((ONE-OLAP-sqrt(ONE-Rout**2))/(ONE-OLAP-sqrt(ONE-Rkn**2)))
483              EVAL_H_PFP = TWO*PI*(TERM1+TERM2)
484              
485           ENDIF
486           RETURN
487           END FUNCTION EVAL_H_PFP
488     
489           END FUNCTION DES_CONDUCTION
490     
491     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
492     !                                                                      !
493     !  Module name: DES_CONDUCTION                                         !
494     !                                                                      !
495     !  Purpose: Compute conductive heat transfer with wall                 !
496     !                                                                      !
497     !  Author: A.Morris                                  Date: 07-Jan-16   !
498     !                                                                      !
499     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
500           FUNCTION DES_CONDUCTION_WALL(OLAP, K_sol, K_wall, K_gas, TWall, &
501           TPart, Rpart, RLens, M)
502           USE param1, only: zero
503           USE DES_THERMO, only: des_min_cond_dist
504     
505           DOUBLE PRECISION :: DES_CONDUCTION_WALL 
506           DOUBLE PRECISION, intent(in) :: OLAP, K_sol, K_wall, K_gas, TWall&
507           &                               ,TPart,RPart, RLens
508           DOUBLE PRECISION :: Rin, Rout
509           DOUBLE PRECISION :: OLAP_ACTUAL, TAUC_CORRECTION
510           DOUBLE PRECISION :: KEff
511           DOUBLE PRECISION :: Q_pw, Q_pfw
512           INTEGER :: M ! Solids phase index
513     
514           ! Initialize variables
515           Q_pw = ZERO
516           Q_pfw= ZERO
517           OLAP_ACTUAL = OLAP
518     
519           ! Check to see if particle is in contact (P-W) conduction
520           IF(OLAP.gt.ZERO)THEN
521           ! Compute effective solids conductivity
522              KEff = 2.0D0*K_sol*K_wall/(K_sol+K_wall)
523           ! Correct for overlap using "area" correction terms
524              IF(DO_AREA_CORRECTION)THEN
525                 OLAP_ACTUAL = CORRECT_OLAP(OLAP,M,-1) ! (-1 indicates wall)
526              ENDIF
527           ! Compute inner Rad
528              Rin = sqrt(2.0D0*OLAP_ACTUAL*Rpart-OLAP_ACTUAL*OLAP_ACTUAL)
529           ! Compute time correction term (commented out for now)
530             ! CALL CALC_TIME_CORRECTION(TAUC_CORRECTION,M,-1)
531           ! Compute heat transfer (particle-wall)
532              Q_pw = 2.0D0*Rin*Keff*(TWall-TPart)
533           ENDIF
534           ! Compute heat transfer (particle-fluid-wall)
535           Q_pfw=EVAL_H_PFW(RLens,DES_MIN_COND_DIST, OLAP_ACTUAL, RPart) * &
536           &     K_gas*RPart*(TWall-TPart)
537     
538           DES_CONDUCTION_WALL=Q_pw+Q_pfw
539          
540           RETURN 
541           END FUNCTION DES_CONDUCTION_WALL
542     
543     
544     !     THIS FUNCTION COMPUTES THE OVERLAP THAT WOULD OCCUR (for a given force)
545     !     IF ACTUAL MATERIAL PROPERTIES WERE USED. 
546           FUNCTION CORRECT_OLAP(OLAP,M,L)
547           use discretelement
548           IMPLICIT NONE
549           DOUBLE PRECISION :: CORRECT_OLAP
550           DOUBLE PRECISION, INTENT (IN) :: OLAP
551           INTEGER, INTENT (IN) :: M, L
552           DOUBLE PRECISION :: KN_ACTUAL, KN_SIM
553           ! L=-1 corresponds to wall
554           IF(L.eq.-1)THEN
555              ! WALL CONTACT
556              IF(DES_COLL_MODEL_ENUM .EQ. HERTZIAN)THEN
557                 CORRECT_OLAP = (HERT_KWN(M)/HERT_KWN_ACTUAL(M))**(2.0D0/3.0D0)*OLAP
558              ELSE
559                 CORRECT_OLAP = (KN_W*OLAP/HERT_KWN_ACTUAL(M))**(2.0D0/3.0D0)
560              ENDIF
561           ELSE
562              ! PARTICLE-PARTICLE
563              IF(DES_COLL_MODEL_ENUM .EQ. HERTZIAN)THEN
564                 CORRECT_OLAP = (HERT_KN(M,L)/HERT_KN_ACTUAL(M,L))**(2.0D0/3.0D0)*OLAP
565              ELSE
566                 CORRECT_OLAP = (KN*OLAP/HERT_KN_ACTUAL(M,L))**(2.0D0/3.0D0)
567              ENDIF
568           ENDIF
569           RETURN
570           END FUNCTION CORRECT_OLAP
571     
572     
573           DOUBLE PRECISION FUNCTION EVAL_H_PFW(RLENS_dim,S,OLAP_dim,RP)
574           USE CONSTANT
575           USE PARAM1
576          
577           IMPLICIT NONE
578           ! Note: Function inputs dimensional quantities
579           DOUBLE PRECISION, intent(in) :: RLENS_dim, S, OLAP_dim, RP
580           ! BELOW VARIABLES ARE NONDIMENSIONALIZED BY RP
581           DOUBLE PRECISION :: RLENS, OLAP, KN
582           DOUBLE PRECISION :: TERM1,TERM2,TERM3
583           DOUBLE PRECISION :: Rout,Rkn
584           DOUBLE PRECISION, PARAMETER :: TWO = 2.0D0
585           
586           RLENS = RLENS_dim/RP
587           KN = S/RP
588           OLAP = OLAP_dim/RP
589     
590           IF(-OLAP.ge.(RLENS-ONE))THEN
591              Rout = ZERO
592           ELSEIF(OLAP.le.TWO)THEN
593              Rout = sqrt(RLENS**2-(ONE-OLAP)**2)
594           ELSE
595              WRITE(*,*)'ERROR: Extremeley excessive overlap (Olap > Diam)'
596              WRITE(*,*)'OLAP = ',OLAP_dim,OLAP, RP
597              WRITE(*,*)'RLENS_dim',RLENS_dim, RLENS
598              write(*,*)'S, Kn', S,KN
599              STOP
600           ENDIF
601           Rout=MIN(Rout,ONE)
602           
603           IF(OLAP.ge.ZERO)THEN
604          !     Particle in contact (below code verified)
605              TERM1 = PI*((ONE-OLAP)**2-(ONE-OLAP-KN)**2)/KN
606              TERM2 = TWO*PI*(sqrt(ONE-Rout**2)-(ONE-OLAP-KN))
607              TERM3 = TWO*PI*(ONE-OLAP)*log((ONE-OLAP-sqrt(ONE-Rout**2))/Kn)
608              EVAL_H_PFW = TERM1+TERM2+TERM3
609           ELSE
610              IF(-OLAP.ge.KN)THEN
611                 Rkn = ZERO
612              ELSE
613                 Rkn=sqrt(ONE-(ONE-OLAP-Kn)**2)
614              ENDIF
615         
616              TERM1 = (Rkn**2/(TWO*KN))+sqrt(ONE-Rout**2)-sqrt(ONE-Rkn**2)
617              TERM2 = (ONE-OLAP)*log((ONE-OLAP-sqrt(ONE-Rout**2))/(ONE-OLAP-sqrt(ONE-Rkn**2)))
618              EVAL_H_PFW = TWO*PI*(TERM1+TERM2)
619              
620           ENDIF
621           RETURN
622           END FUNCTION EVAL_H_PFW
623     
624     
625           SUBROUTINE CALC_TIME_CORRECTION(time_corr, phaseI, phaseJ)
626           use discretelement, only: tau_c_base_actual, tauw_c_base_actual
627           use discretelement, only: tau_c_base_sim, tauw_c_base_sim
628           use discretelement, only: des_coll_model_enum, HERTZIAN
629           IMPLICIT NONE
630           INTEGER, intent (in) :: phaseI, phaseJ
631           DOUBLE PRECISION :: time_corr, tau_actual, tau_sim
632           DOUBLE PRECISION :: vimp ! impact veloctiy
633           INTEGER :: M, N
634           vimp = 1.0D0
635           time_corr = 1.0D0
636           if (phaseJ == -1) then
637              ! Wall contact
638              M = phaseI
639              IF (DES_COLL_MODEL_ENUM .EQ. HERTZIAN) THEN
640                 time_corr = tauw_c_base_actual(M) / tauw_c_base_sim(M)
641              ELSE
642                 vimp = 1.0D0
643                 if (vimp .le. 0.0D0)then
644                    time_corr = 1.0D0
645                 else
646                    time_corr = tauw_c_base_actual(M)*vimp**(-0.2D0) / tauw_c_base_sim(M)
647                 endif
648              ENDIF
649     
650           ELSE
651              ! particle-particle contact
652              if (phaseI .le. phaseJ)then
653                 M = phaseI
654                 N = phaseJ
655              else
656                 M = phaseJ
657                 N = phaseI
658              endif
659              
660              IF (DES_COLL_MODEL_ENUM .EQ. HERTZIAN) THEN
661                 time_corr = tau_c_base_actual(M,N) / tau_c_base_sim(M,N)
662              ELSE
663                 vimp = 1.0D0
664                 if (vimp .le. 0.0D0)then
665                    time_corr = 1.0D0
666                 else
667                    time_corr = tau_c_base_actual(M,N)*vimp**(-0.2D0) / tau_c_base_sim(M,N)
668                 endif
669              ENDIF
670           ENDIF
671           time_corr = time_corr **(2.0D0/3.0D0)
672           ! TEMPORARY 
673           time_corr = 1.0D0
674           RETURN
675           END SUBROUTINE CALC_TIME_CORRECTION
676     
677     
678           
679       END MODULE DES_THERMO_COND
680