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 CONTAINS 34 35 FUNCTION DES_CONDUCTION(I, J, CENTER_DIST, iM, iIJK) 36 37 use constant 38 use des_thermo 39 use discretelement 40 use funits 41 use physprop 42 43 IMPLICIT NONE 44 45 ! Passed variables 46 !---------------------------------------------------------------------// 47 ! Index of particle being looped over 48 INTEGER, INTENT(IN) :: I,J 49 ! Index of fluid cell containing particle I 50 INTEGER, INTENT(IN) :: iIJK 51 ! Solid phase indices for the given particles 52 INTEGER, INTENT(IN) :: iM 53 ! Distance between the centers of particle I and particle J (magnitude) 54 DOUBLE PRECISION, INTENT(IN) :: CENTER_DIST 55 DOUBLE PRECISION :: DES_CONDUCTION 56 57 ! Local variables 58 !---------------------------------------------------------------------// 59 ! Solid phase indices for the given particles 60 INTEGER jM 61 ! Radius of smaller particle 62 DOUBLE PRECISION MIN_RAD 63 ! Radius of larger particle 64 DOUBLE PRECISION MAX_RAD 65 ! Rate of particle-particle conduction 66 DOUBLE PRECISION Q_pp 67 ! Rate of particle-fluid-particle conduction 68 DOUBLE PRECISION Q_pfp 69 ! Outer radius of region delineating particle-fluid-particle conduction 70 DOUBLE PRECISION RD_OUT 71 ! Inner radius of region delineating particle-fluid-particle conduction 72 DOUBLE PRECISION RD_IN 73 ! The radius of the fluid lens containing the larger particle 74 DOUBLE PRECISION LENS_RAD 75 ! Temperature difference between two particles 76 DOUBLE PRECISION DeltaTp 77 ! Effective thermal conductivity 78 DOUBLE PRECISION lK_eff 79 ! Radius of contact area 80 DOUBLE PRECISION lRadius 81 82 ! Functions 83 !---------------------------------------------------------------------// 84 85 ! Identify the solid phases of the neighbor particle 86 jM = PIJK(J,5) + SMAX 87 88 ! Determine the radius of the larger and smaller particle 89 MIN_RAD = MIN(DES_RADIUS(I), DES_RADIUS(J)) 90 MAX_RAD = MAX(DES_RADIUS(I), DES_RADIUS(J)) 91 92 DeltaTp = DES_T_s_NEW(J) - DES_T_s_NEW(I) 93 94 ! Calculate the particle-particle conduction 95 ! REF: Batchelor and O'Brien, 1977 (MODIFIED) 96 !---------------------------------------------------------------------// 97 Q_pp = 0.0 98 IF(CENTER_DIST < (MAX_RAD + MIN_RAD)) THEN 99 ! Effective thermal conductivity 100 lK_eff = K_eff(K_s0(iM),K_s0(jM)) 101 ! Effective contact area's radius 102 lRadius = RADIUS(MAX_RAD, MIN_RAD) 103 ! Inter-particle heat transfer 104 Q_pp = 2.0d0 * lK_eff * lRadius * DeltaTp 105 ! Assign the inter-particle heat transfer to both particles. 106 ENDIF 107 108 ! Calculate the particle-fluid-particle conduction 109 ! REF: Rong and Horio, 1999 (MODIFIED) 110 !---------------------------------------------------------------------// 111 ! Calculate the radius of the fluid lens surrounding the larger particle 112 ! Default FLPC = 0.2 113 LENS_RAD = MAX_RAD * (1.0D0 + FLPC) 114 115 ! Calculate the outer radial distance of the region for particle-fluid- 116 ! particle heat conduction. 117 RD_OUT = RADIUS( LENS_RAD, MIN_RAD) 118 119 ! If the value returned is less than zero, then the fluid lens 120 ! surrounding the larger particle does not intersect with the surface 121 ! of the smaller particle. In this case, particle-fluild-particle 122 ! conduction does not occur. 123 Q_pfp = 0.0 124 IF(RD_OUT .GT. ZERO)THEN 125 ! Calculate the distance from the line connecting the particles' centers 126 ! to the point of contact between the two particles. This value is 127 ! zero if the particles are not touching and is the radius of the 128 ! shared contact area otherwise. 129 RD_IN = ZERO 130 IF(CENTER_DIST < (MAX_RAD + MIN_RAD) ) & 131 RD_IN = RADIUS(MAX_RAD, MIN_RAD) 132 ! Calculate the rate of heat transfer between the particles through the 133 ! fluid using adaptive Simpson's rule to manage the integral. 134 Q_pfp = K_g(iIJK) * DeltaTp * & 135 ADPT_SIMPSON(RD_IN,RD_OUT) 136 137 ENDIF ! RD_OUT > ZERO 138 139 ! Return total heat flux 140 DES_CONDUCTION = Q_pp + Q_pfp 141 142 RETURN 143 144 CONTAINS 145 146 !......................................................................! 147 ! Subroutine Procedure: FUNCTION RADIUS ! 148 ! ! 149 ! Propose: Calculate the center line radius between two particles. ! 150 ! Depending on what is given as input, this function calculates: ! 151 ! 1) radius of contact area between two particles ! 152 ! 2) radius delineating particle-fluid-particle region. ! 153 !......................................................................! 154 DOUBLE PRECISION FUNCTION RADIUS(R1, R2) 155 156 IMPLICIT NONE 157 158 ! Radius values 159 DOUBLE PRECISION, INTENT(IN) :: R1, R2 160 ! Distance between particle centers 161 DOUBLE PRECISION BETA 162 ! Generic value 163 DOUBLE PRECISION VALUE 164 165 ! Calculate 166 VALUE = (R1**2 - R2**2 + CENTER_DIST**2)/(2.0d0 * R1 * CENTER_DIST) 167 ! Check to ensure that VALUE is less than or equal to one. If VALUE is 168 ! greater than one, the triangle inequality has been violated. Therefore 169 ! there is no intersection between the fluid lens surrounding the larger 170 ! particle and the surface of the smaller particle. 171 ! Thus, there is no particle-fluid-particle heat transfer. 172 IF( VALUE .GT. 1.0d0) THEN 173 RADIUS = -1.0d0 174 ELSE 175 ! Calculate beta (Law of cosines) 176 BETA = ACOS( VALUE ) 177 ! Calculate the radius 178 RADIUS = R1 * SIN(BETA) 179 ENDIF 180 181 RETURN 182 END FUNCTION RADIUS 183 184 !......................................................................! 185 ! Subroutine Procedure: FUNCTION K_eff ! 186 ! ! 187 ! Propose: Calculate the effective thermal conductivity of two ! 188 ! particles with conductivities K1 and K2 ! 189 !......................................................................! 190 DOUBLE PRECISION FUNCTION K_eff(K1, K2) 191 192 IMPLICIT NONE 193 194 ! Thermal conductivities 195 DOUBLE PRECISION, INTENT(IN) :: K1, K2 196 197 K_eff = 2.0d0 * (K1*K2)/(K1 + K2) 198 199 RETURN 200 END FUNCTION K_eff 201 202 !``````````````````````````````````````````````````````````````````````! 203 ! Function: F ! 204 ! ! 205 ! Purpose: This function defines the region of particle-fluid-particle ! 206 ! heat transfer. Note that this function develops a ! 207 ! singularity at the boundary of the contact region when two ! 208 ! particles are touching. This singularity is resolved by ! 209 ! making the assumption that the surfaces of two particles ! 210 ! never directly touch (c.f. Rong and Horio, 1999). By ! 211 ! default, it is assumed that the particles are separated by ! 212 ! a minimum length of 4.0x10^(-10) meters. This value may be ! 213 ! modified by specifying a different value in the mfix.dat ! 214 ! file under the variable name DES_MIN_COND_DIST. ! 215 ! ! 216 ! Note, the closer this value is to zero, the harder it will ! 217 ! be for the quadrature routine to converge. ! 218 !......................................................................! 219 DOUBLE PRECISION FUNCTION F(R) 220 USE utilities 221 IMPLICIT NONE 222 223 DOUBLE PRECISION, INTENT(IN) :: R 224 225 F = (2.0d0*Pi*R)/MAX((CENTER_DIST - SQRT(MAX_RAD**2-R**2) - & 226 SQRT(MIN_RAD**2-R**2)), DES_MIN_COND_DIST) 227 228 IF( mfix_isnan(F))PRINT*,'F IS NAN' 229 230 END FUNCTION F 231 232 !``````````````````````````````````````````````````````````````````````! 233 ! Function: ADPT_SIMPSON ! 234 ! ! 235 ! Purpose: This function applies adaptive Simpson's Rule to solve the ! 236 ! definite integral of the function F. ! 237 ! ! 238 ! *NOTE: This route will return the result of the method even ! 239 ! if convergence has not been achieved. If this occurs, a ! 240 ! message will be written to the log file to notify the user. ! 241 !......................................................................! 242 DOUBLE PRECISION FUNCTION ADPT_SIMPSON(A, B) 243 244 IMPLICIT NONE 245 246 ! Bounds of integration 247 DOUBLE PRECISION, INTENT(IN) :: A, B 248 ! Maximum recursive depth 249 INTEGER, PARAMETER :: Kmax = 25 250 ! Current depth of recursion 251 INTEGER :: K 252 ! Array storage for recursion 253 DOUBLE PRECISION :: V(Kmax,6) 254 ! Step size 255 DOUBLE PRECISION :: H 256 ! Simpson's Rule evaluated on the left and right intervals 257 DOUBLE PRECISION :: lS, rS 258 ! Value of the function F evaluate at the midpoints of the left and 259 ! right intervals 260 DOUBLE PRECISION :: rF, lF 261 ! Local error bound 262 DOUBLE PRECISION :: Err_BND 263 ! Convergence bound 264 DOUBLE PRECISION :: EPS = 10.0**(-2) 265 266 ! Error indicating that an error has been specified to the user. 267 LOGICAL, SAVE :: ADPT_SIMPSON_ERR = .FALSE. 268 269 ! Initialize variables 270 V(:,:) = 0.0d0 ! 271 ADPT_SIMPSON = 0.0d0 ! Integral value 272 H = (B-A)/2.0d0 ! Dynamic interval length 273 ! Calculate Simpson's Rule over the interval [A,B] 274 lS = (F(A) + 4.0d0*F((A+B)/2.0d0) + F(B))*(H/3.0d0) 275 ! Initialize the storage vector for the first pass 276 V(1,:) = (/ A, H, F(A), F((A+B)/2.0d0), F(B), lS/) 277 K = 1 278 ! Enter recursion calculations 279 DO 280 ! Establish the new interval length 281 H = V(K,2)/2.0d0 282 ! Evaluate the function on the left interval 283 lF = F(V(K,1) + H) 284 ! Calculate Simpson's Rule over the left interval 285 lS = (V(K,3) + 4.0d0*lF + V(K,4))*(H/3.0d0) 286 ! Evaluate the function on the right interval 287 rF = F(V(K,1) + 3.0d0*H) 288 ! Calculate Simpson's Rule over the right interval 289 rS = (V(K,4) + 4.0d0*rF + V(K,5))*(H/3.0d0) 290 ! Evaluate the error bound 291 Err_BND = (30.0d0*EPS*H)/(B-A) 292 ! Check to see if conversion on the interval has been met 293 IF( ABS(lS + rS - V(K,6)) .LT. Err_BND)THEN 294 ! Update the integral value 295 ADPT_SIMPSON = ADPT_SIMPSON + lS + rS + & 296 (1.0d0/15.0d0)*(lS + rS - V(K,6)) 297 ! Decrement the recursion counter 298 K = K - 1 299 ! If K=0, then integration has been successfully completed over the 300 ! entire interval [A,B]. 301 IF( K == 0) RETURN 302 ELSEIF( (K .GE. Kmax) .OR. & 303 (H == (B-A)*(1.0d0/2.0d0)**(Kmax+3))) THEN 304 ! Flag that the method did not converge. 305 IF(.NOT.ADPT_SIMPSON_ERR)THEN 306 WRITE(*,1000) 307 WRITE(UNIT_LOG,1000) 308 ADPT_SIMPSON_ERR = .TRUE. 309 ENDIF 310 ! Update the integral value 311 ADPT_SIMPSON = ADPT_SIMPSON + lS + rS + & 312 (1.0d0/15.0d0)*(lS + rS - V(K,6)) 313 ! Decrement the recursion counter 314 K = K - 1 315 ELSE 316 ! Refine the subintervals through recursive splitting of the intervals 317 V(K+1,:) = (/V(K,1) + 2.0d0*H, H, V(K,4), rF, V(K,5), rS/) 318 V(K,:) = (/V(K,1), H, V(K,3), lF, V(K,4), lS/) 319 ! Increment level counter 320 K = K+1 321 ENDIF 322 ENDDO 323 324 1000 FORMAT(/1X,70('*'),/' From: DES_COND_EQ',/, ' Message: ', & 325 'Integration of the particle-fluid-particle equation did ', & 326 'not',/' converge! No definite bound can be placed on the ', & 327 'error.',/' Future convergence messages will be suppressed!', & 328 /1X,70('*')) 329 330 END FUNCTION ADPT_SIMPSON 331 332 END FUNCTION DES_CONDUCTION 333 334 END MODULE DES_THERMO_COND 335