File: N:\mfix\model\set_constprop.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: set_constprop                                           C
4     !  Purpose: This routine serves two purposes:                          C
5     !    1) initializes various variables everywhere in the domain with    C
6     !       a zero value. the real need for this is unclear. undefined     C
7     !       may be a better approach...                                    C
8     !    2) if defined, sets physical properties to their specified        C
9     !       constant value in the fluid domain. cannot set in flow         C
10     !       boundaries or later checks will also complain (may be          C
11     !       overly strict check)                                           C
12     !                                                                      C
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14           SUBROUTINE SET_CONSTPROP
15     
16     ! Modules
17     !-----------------------------------------------
18           USE param1, only: zero, half, one, undefined
19     
20           use discretelement, only: des_mmax
21     
22           USE fldvar, only: ro_g
23           USE fldvar, only: ro_s, d_p
24           use fldvar, only: p_s
25     
26           USE visc_g, only: mu_gt, epmu_gt, lambda_gt, eplambda_gt
27           USE visc_s, only: mu_s, epmu_s, lambda_s, eplambda_s, lambda_s_c
28           USE visc_s, only: ep_star_array
29           USE visc_s, only: ep_g_blend_start, ep_g_blend_end
30     
31           USE physprop, only: nmax, mmax, smax
32           USE physprop, only: ro_s0, ro_g0, d_p0
33           USE physprop, only: mu_g0, mu_g, mu_s0
34           USE physprop, only: mw_avg, mw_mix_g
35           USE physprop, only: c_pg0, k_g0, c_pg, k_g, dif_g0, dif_g
36           USE physprop, only: c_ps0, k_s0, c_ps, k_s, dif_s0, dif_s
37           USE physprop, only: cv
38     
39           USE constant, only: ep_s_max_ratio, d_p_ratio, ep_s_max, m_max
40           use constant, only: ep_star
41     
42           USE run, only: call_dqmom
43           USE run, only: yu_standish, fedors_landel
44           USE run, only: kt_type_enum, ia_2005, gd_1999, gtsh_2012
45           USE run, only: blending_stress, sigm_blend, tanh_blend
46     
47           USE drag, only: f_gs, f_ss
48     
49           use kintheory, only: mu_sm_ip, mu_sl_ip, xi_sm_ip, xi_sl_ip
50           use kintheory, only: fnu_s_ip, ft_sm_ip, ft_sl_ip
51           use kintheory, only: kth_sl_ip, knu_sm_ip, knu_sl_ip, kvel_s_ip
52           use kintheory, only: ed_ss_ip, edvel_sl_ip, edt_s_ip, edvel_sm_ip
53           use kintheory, only: a2_gtsh, xsi_gtsh
54     
55           use mms, only: use_mms
56     
57           USE compar, only: ijkstart3, ijkend3
58           use functions, only: wall_at, fluid_at
59     
60           IMPLICIT NONE
61     !-----------------------------------------------
62     ! Local variables
63     !-----------------------------------------------
64     ! indices
65           INTEGER :: IJK, M, I, J
66           DOUBLE PRECISION :: old_value, DP_TMP(SMAX)
67     !-----------------------------------------------
68     
69     ! First, initialize certain transport coefficients, physical
70     ! properties, and other variable types everywhere in the
71     ! domain to zero. Then, set these to their specified value
72     ! (if defined) in fluid cells. Some are also set in flow cells.
73     ! NOTE: DO NOT simply zero existing field variables.
74           RO_g = ZERO
75           MU_g = ZERO
76           MU_gt = ZERO
77           EPMU_GT = ZERO
78           LAMBDA_GT = ZERO
79           EPLAMBDA_GT = ZERO
80           K_g = ZERO
81           C_PG = ZERO
82           MW_MIX_G = zero
83           DIF_g = ZERO
84           F_GS = ZERO
85     
86           RO_S = ZERO
87           D_P = zero     ! this is done in init_fvars
88           MU_s = ZERO
89           EPMU_s = ZERO
90           LAMBDA_s_c = ZERO
91           LAMBDA_s = ZERO
92           EPLAMBDA_s = ZERO
93     ! not defining p_s in flow cells can become a problem for certain
94     ! cases as gradients in P_s are used (e.g., when a PO is west of a
95     ! MI...?). Simply assigning it a zero value may not be the best
96     ! approach but it is currently used....
97           P_S = ZERO
98           K_s = ZERO
99           C_PS = ZERO
100           DIF_S = ZERO
101           F_SS = ZERO
102     
103     ! Set the flag for recalculating gas viscosity.
104     !      RECALC_VISC_G = (MU_g0==UNDEFINED .OR. L_SCALE0/=ZERO .OR.&
105     !                       K_EPSILON .OR. ISHII)
106     
107     ! Set default value for virtual mass coefficient
108           Cv = HALF
109     
110     ! Variables specific to various kinetic theory models
111           IF (KT_TYPE_ENUM == IA_2005) THEN
112              MU_sM_ip = ZERO
113              MU_sL_ip = ZERO
114              XI_sM_ip = ZERO
115              XI_sL_ip = ZERO
116              Fnu_s_ip = ZERO
117              FT_sM_ip = ZERO
118              FT_sL_ip = ZERO
119              Kth_sL_ip = ZERO
120              Knu_sM_ip = ZERO
121              Knu_sL_ip = ZERO
122              Kvel_s_ip = ZERO
123              ED_ss_ip = ZERO
124              EDvel_sL_ip = ZERO
125           ENDIF
126     
127           IF (KT_TYPE_ENUM == IA_2005 .OR. &
128               KT_TYPE_ENUM == GD_1999 .OR.  &
129               KT_TYPE_ENUM == GTSH_2012) THEN
130              EDT_s_ip = ZERO
131              EDvel_sM_ip = ZERO
132           ENDIF
133     
134           IF(KT_TYPE_ENUM == GTSH_2012) THEN
135              A2_gtsh = ZERO
136              xsi_gtsh = zero
137           ENDIF
138     
139     ! Set specified constant physical properties values
140           DO IJK = ijkstart3, ijkend3
141     
142              IF (.NOT.WALL_AT(IJK)) THEN
143     ! Fluid and inflow/outflow cells: FLAG < 100
144                 IF (RO_G0 /= UNDEFINED) RO_G(IJK) = RO_G0
145                 IF (C_PG0 /= UNDEFINED) C_PG(IJK) = C_PG0
146                 IF (MW_AVG /= UNDEFINED) MW_MIX_G(IJK) = MW_AVG
147              ENDIF
148     
149              IF (FLUID_AT(IJK)) THEN
150     ! Strictly Fluid cells: FLAG = 1
151                 IF (MU_G0 /= UNDEFINED) THEN
152                    MU_G(IJK) = MU_G0
153                    MU_GT(IJK) = MU_G0
154                    EPMU_GT(IJK) = MU_G0
155                    LAMBDA_GT(IJK) = -(2.0d0/3.0d0)*MU_G0
156                    EPLAMBDA_GT(IJK) = -(2.0d0/3.0d0)*MU_G0
157                 ENDIF
158                 IF (K_G0 /= UNDEFINED) K_G(IJK) = K_G0
159                 IF (DIF_G0 /= UNDEFINED) DIF_G(IJK,:NMAX(0)) = DIF_G0
160              ENDIF
161     
162              IF (USE_MMS) THEN
163                 IF (MU_G0 /= UNDEFINED) THEN
164                    MU_G(IJK) = MU_G0
165                    MU_GT(IJK) = MU_G0
166                    EPMU_GT(IJK) = MU_G0
167                    LAMBDA_GT(IJK) = -(2.0d0/3.0d0)*MU_G0
168                    EPLAMBDA_GT(IJK) = -(2.0d0/3.0d0)*MU_G0
169                 ENDIF
170                 IF (K_G0 /= UNDEFINED) K_G(IJK) = K_G0
171                 IF (DIF_G0 /= UNDEFINED) DIF_G(IJK,:NMAX(0)) = DIF_G0
172              ENDIF
173     
174           ENDDO
175     
176     ! ensure ro_s(ijk,m) is assigned to ro_s0(m) or ro_s(np) so that the
177     ! function ep_s works for discrete phases. might be able to move this
178     ! to set_ic_dem and set_ic_mppic. also required d_p(ijk,m) for hybrid
179     ! use.  note check_solids_common_all ensures d_p0 is set for all
180     ! solids defined also either ro_s0 must be set or ro_s0
181           DO M = 1, MMAX+DES_MMAX
182              DO IJK = ijkstart3, ijkend3
183                 IF(.NOT.WALL_AT(IJK)) THEN
184     ! Fluid and inflow/outflow cells: FLAG < 100
185                    IF (RO_S0(M) /= UNDEFINED) RO_S(IJK,M) = RO_S0(M)
186                    IF (C_PS0(M) /= UNDEFINED) C_PS(IJK,M) = C_PS0(M)
187                    IF (D_P0(M) /= UNDEFINED) D_P(IJK,M) = D_P0(M)
188                 ENDIF
189     
190                 IF (FLUID_AT(IJK)) THEN
191     ! Strictly fluid cells: FLAG = 1
192                    IF (MU_S0(M) /= UNDEFINED) THEN
193                       MU_S(IJK,M) = MU_S0(M)
194                       EPMU_S(IJK,M) = MU_S0(M)
195                       LAMBDA_S(IJK,M) = (-2./3.)*MU_S(IJK,M)
196                       EPLAMBDA_S(IJK,M) = (-2./3.)*MU_S(IJK,M)
197                    ENDIF
198                    IF (K_S0(M) /= UNDEFINED) K_S(IJK,M) = K_S0(M)
199                    IF (DIF_S0(M) /= UNDEFINED) DIF_S(IJK,M,:NMAX(M)) = DIF_S0(M)
200                 ENDIF
201     
202                 IF (USE_MMS) THEN
203                    IF (MU_S0(M) /= UNDEFINED) THEN
204                       MU_S(IJK,M) = MU_S0(M)
205                       EPMU_S(IJK,M) = MU_S0(M)
206                       LAMBDA_S(IJK,M) = (-2./3.)*MU_S(IJK,M)
207                       EPLAMBDA_S(IJK,M) = (-2./3.)*MU_S(IJK,M)
208                    ENDIF
209                    IF (K_S0(M) /= UNDEFINED) K_S(IJK,M) = K_S0(M)
210                    IF (DIF_S0(M) /= UNDEFINED) DIF_S(IJK,M,:NMAX(M)) = DIF_S0(M)
211                 ENDIF
212     
213     ! set ep_star_array to user input ep_star in all cells.
214                 EP_star_array(ijk) = ep_star
215     ! initializing blending stress parameters
216                 IF(BLENDING_STRESS.AND.TANH_BLEND) THEN
217                    ep_g_blend_start(ijk) = ep_star_array(ijk) * 0.99d0
218                    ep_g_blend_end(ijk)   = ep_star_array(ijk) * 1.01d0
219                 ELSE IF(BLENDING_STRESS.AND.SIGM_BLEND) THEN
220                    ep_g_blend_start(ijk) = ep_star * 0.97d0
221                    ep_g_blend_end(ijk) = ep_star * 1.01d0
222                 ELSE
223                    ep_g_blend_start(ijk) = ep_star_array(ijk)
224                    ep_g_blend_end(ijk)   = ep_star_array(ijk)
225                 ENDIF
226     
227              ENDDO   ! end loop over ijk
228           ENDDO   ! end loop over MMAX
229     
230     
231     ! Initializing parameters needed if a correlation is used to compute
232     ! ep_star: initializing the indexing system.
233           IF(YU_STANDISH .OR. FEDORS_LANDEL) THEN
234              DO M = 1, SMAX
235                 IF(EP_S_MAX(M) == UNDEFINED) EP_S_MAX(M) = ONE-EP_STAR
236              ENDDO
237     
238              IF (.NOT.CALL_DQMOM) THEN
239     
240     ! refer to Syam's dissertation
241                 IF (SMAX == 2) THEN
242                    ep_s_max_ratio(1,2) = ep_s_max(1)/ &
243                       (ep_s_max(1)+(1.-ep_s_max(1))*ep_s_max(2))
244                 ENDIF
245     
246     ! initialize local variables
247                 DO I = 1, SMAX
248                    DP_TMP(I) = D_P0(I)
249                    M_MAX(I) = I
250                 ENDDO
251     
252     ! Rearrange the indices from coarsest particles to finest to be
253     ! used in CALC_ep_star. Done here because it may need to be done
254     ! for auto_restart
255                 DO I = 1, SMAX
256                    DO J = I, SMAX
257                       IF(DP_TMP(I) < DP_TMP(J)) THEN
258                          old_value = DP_TMP(I)
259                          DP_TMP(I) = DP_TMP(J)
260                          DP_TMP(J) = old_value
261                       ENDIF
262                    ENDDO
263                 ENDDO
264     
265                 DO I = 1, SMAX
266                    DO J = 1, SMAX
267                       IF(DP_TMP(I) == D_P0(J) .AND. D_P0(I) .NE. D_P0(J)) THEN
268                          M_MAX(I) = J
269                       ENDIF
270                    ENDDO
271                 ENDDO
272              ENDIF    ! if .not. call_dqmom
273           ELSE   ! if .not. Yu-standish or Fedors-Landel
274              EP_S_MAX(:) = ZERO
275              EP_S_MAX_RATIO(:,:) = ZERO
276              D_P_RATIO(:,:) = ZERO
277              M_MAX(:) = 0
278           ENDIF
279     
280           RETURN
281           END SUBROUTINE SET_CONSTPROP
282