1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv! 2 ! ! 3 ! Module name: DES_REACTION_MODEL ! 4 ! ! 5 ! Purpose: ! 6 ! ! 7 ! ! 8 ! Author: J.Musser Date: 16-Jun-10 ! 9 ! ! 10 ! Comments: ! 11 ! ! 12 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! 13 SUBROUTINE DES_REACTION_MODEL 14 15 USE compar 16 Use constant 17 Use des_rxns 18 Use des_thermo 19 Use discretelement 20 USE geometry 21 USE indices 22 Use param1 23 use run, only: ANY_SPECIES_EQ, SPECIES_EQ 24 use physprop, only: SMAX, NMAX 25 use run, only: SOLVE_ROs 26 use functions 27 28 IMPLICIT NONE 29 30 ! Passed variables 31 !----------------------------------------------- 32 ! None 33 34 ! Local variables 35 !----------------------------------------------- 36 ! Index of neighbor particle of particle I such that I < J 37 INTEGER IJK 38 ! Index value of particle 39 INTEGER lNP, NP 40 ! index of solids phase and species 41 INTEGER M, N 42 ! total rate of consumption/production of species (g/sec) 43 DOUBLE PRECISION SUM_DES_R_sc, SUM_DES_Rs 44 ! masses of species comprising the particle (g) 45 DOUBLE PRECISION S_MASS( DIMENSION_N_s ) 46 47 LOGICAL ALL_GONE( DIMENSION_N_s ) 48 49 DOUBLE PRECISION, PARAMETER :: P43 = 4.0d0/3.0d0 50 51 DOUBLE PRECISION dMdt, dXdt, dXMdt, dXMdt_DTs 52 53 ! Logical for Adams-Bashfort integration. 54 LOGICAL,SAVE:: FIRST_PASS = .TRUE. 55 56 !---------------------------------------------------------------------// 57 58 IF(.NOT.ANY_SPECIES_EQ) RETURN 59 60 ! Loop over fluid cells. 61 !---------------------------------------------------------------------// 62 IJK_LP: DO IJK = IJKSTART3, IJKEND3 63 IF(.NOT.FLUID_AT(IJK)) CYCLE IJK_LP 64 IF(PINC(IJK) == 0) CYCLE IJK_LP 65 66 ! Loop over all particles in cell IJK. 67 !---------------------------------------------------------------------// 68 lNP_LP: DO lNP = 1, PINC(IJK) 69 NP = PIC(IJK)%p(lNP) 70 ! Skip indices that do not represent particles 71 IF(IS_NONEXISTENT(NP)) CYCLE lNP_LP 72 ! Skip indices that represent ghost particles 73 IF(IS_GHOST(NP) .OR. IS_ENTERING_GHOST(NP) .OR. IS_EXITING_GHOST(NP)) CYCLE lNP_LP 74 75 ! Set the particle phase index 76 M = PIJK(NP,5) + SMAX 77 78 ! Skip particles when not solving species equations 79 IF(.NOT.SPECIES_EQ(M)) CYCLE lNP_LP 80 81 ! Check to see if any of the reaction rate will consume more than the 82 ! available species in the particle. 83 !---------------------------------------------------------------------// 84 ! Initialize local variables 85 SUM_DES_Rs = ZERO 86 SUM_DES_R_sc = ZERO 87 ALL_GONE(:) = .FALSE. 88 89 ! Reset the flag 90 DO N=1,NMAX(M) 91 ! Calculate the current mass of species N in the particle 92 S_MASS(N) = DES_X_s(NP,N) * PMASS(NP) 93 ! Calculte the amount of species N mass that is produced/consumed. 94 ! dXMdt_DTs > 0 :: Net production in mass units 95 ! dXMdt_DTs < 0 :: Net consumption in mass units 96 dXMdt_DTs = (DES_R_sp(NP,N)-DES_R_sc(NP,N))*DTSOLID 97 ! Check to see if the amount of solids phase species N consumed will 98 ! be larger than the amount of available solids. 99 IF((S_MASS(N) + dXMdt_DTs) < ZERO) THEN 100 ! Indicate that all of species is consumed. 101 ALL_GONE(N) = .TRUE. 102 ! Limit the consumption rate so that only the amount of species mass 103 ! that the particle contains is consumed. 104 dXMdt = S_MASS(N)/DTSOLID 105 ! Calculate the total rate of change in particle mass 106 SUM_DES_Rs = SUM_DES_Rs + dXMdt 107 ! Calculate the total rate of consumption 108 SUM_DES_R_sc = SUM_DES_R_sc + dXMdt 109 ELSE 110 ! Calculate the total particle mass rate of change (mass/time) 111 SUM_DES_Rs = SUM_DES_Rs + (DES_R_sp(NP,N)-DES_R_sc(NP,N)) 112 ! Calculate the total rate of consumption (mass/time) 113 SUM_DES_R_sc = SUM_DES_R_sc + DES_R_sc(NP,N) 114 ENDIF 115 ENDDO 116 117 ! Update the particle's mass. 118 ! The mass of the particle is updated first so that it can be used in 119 ! updating the species mass percent of the particle. 120 !---------------------------------------------------------------------// 121 dMdt = SUM_DES_Rs 122 ! First-order method: Euler 123 IF(INTG_EULER) THEN 124 PMASS(NP) = PMASS(NP) + DTSOLID*dMdt 125 ! Second-order Adams-Bashforth scheme 126 ELSEIF(INTG_ADAMS_BASHFORTH) THEN 127 IF(FIRST_PASS) dMdt_OLD(NP) = dMdt 128 PMASS(NP) = PMASS(NP) + DTSOLID * & 129 (1.50d0*dMdt - 0.5d0*dMdt_OLD(NP)) 130 ! Store the current value as the old value. 131 dMdt_OLD(NP) = dMdt 132 ENDIF 133 134 ! Update the species mass percent 135 !---------------------------------------------------------------------// 136 DO N=1,NMAX(M) 137 IF(ALL_GONE(N))THEN 138 DES_X_s(NP,N) = ZERO 139 ELSE 140 dXdt = ((DES_R_sp(NP,N) - DES_R_sc(NP,N)) - & 141 DES_X_s(NP,N) * SUM_DES_Rs)/PMASS(NP) 142 ! First-order method: Euler 143 IF(INTG_EULER) THEN 144 DES_X_s(NP,N) = DES_X_s(NP,N) + DTSOLID*dXdt 145 ! Second-order Adams-Bashforth scheme 146 ELSEIF(INTG_ADAMS_BASHFORTH) THEN 147 IF(FIRST_PASS) dXdt_OLD(NP,N) = dXdt 148 DES_X_s(NP,N) = DES_X_s(NP,N) + DTSOLID * & 149 (1.50d0*dXdt - 0.5d0*dXdt_OLD(NP,N)) 150 ! Store the current value as the old value. 151 dXdt_OLD(NP,N) = dXdt 152 ENDIF 153 ENDIF 154 ENDDO 155 156 ! Variable density 157 !---------------------------------------------------------------------// 158 IF(SOLVE_ROs(M)) THEN 159 ! update variables 160 RO_Sol(NP) = PMASS(NP)/PVOL(NP) 161 162 ! Shrinking particle 163 !---------------------------------------------------------------------// 164 ELSE 165 DES_RADIUS(NP) = & 166 (PMASS(NP)/(P43*Pi*Ro_Sol(NP)))**(1.0d0/3.0d0) 167 ! update variables 168 PVOL(NP) = PMASS(NP)/Ro_Sol(NP) 169 ENDIF 170 171 ! Update one over the particle's moment of inertia 172 OMOI(NP) = 5.0d0 / (2.0d0 * PMASS(NP) * DES_RADIUS(NP)**2) 173 174 ENDDO lNP_LP ! End loop over all particles 175 ENDDO IJK_LP ! End loop over fluid cells 176 177 ! Clear the necessary variables. 178 DES_R_sp(:,:) = ZERO 179 DES_R_sc(:,:) = ZERO 180 181 ! Flag that the first pass is over 182 FIRST_PASS = .FALSE. 183 184 RETURN 185 END SUBROUTINE DES_REACTION_MODEL 186