42
43
44
45 USE multi_fvm_mod
46 USE ebcs_mod
48 USE matparam_def_mod , ONLY : matparam_struct_
50 USE output_mod , ONLY : output_
51 use element_mod , only : nixs,nixq,nixtg
52
53
54
55#include "implicit_f.inc"
56
57
58
59
60#include "param_c.inc"
61
62#include "task_c.inc"
63
64#include "com01_c.inc"
65! numels, numelq, numeltg
66#include "com04_c.inc"
67
68 COMMON /tablesizf/ stf,snpc
69 INTEGER STF,SNPC
70
71
72
74 INTEGER, INTENT(IN) :: EBCS_ID
75 TYPE(MULTI_FVM_STRUCT), INTENT(INOUT) :: MULTI_FVM
76 INTEGER, INTENT(IN) :: ITASK, NELEM, ELEM_LIST(NELEM), FACE_LIST(NELEM)
77 INTEGER, INTENT(IN) :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
78 my_real,
INTENT(IN) :: xgrid(3, *), wgrid(3, *)
79 INTEGER, INTENT(IN) :: IPM(NPROPMI, *)
80 my_real,
INTENT(IN) :: pm(npropm, *), func_value(*)
81 TYPE(FVM_INLET_DATA_STRUCT), INTENT(IN) :: FVM_INLET_DATA
82 TYPE(t_ebcs_nrf), INTENT(INOUT) :: EBCS
83 INTEGER, INTENT(IN) :: NPF(SNPC)
84 my_real,
INTENT(IN) :: tf(stf), timestep
85 TYPE(MATPARAM_STRUCT_),DIMENSION(NUMMAT), INTENT(IN) :: MATPARAM
86 TYPE(OUTPUT_),INTENT(INOUT) :: OUTPUT
87 LOGICAL,INTENT(IN) :: PRED
88
89
90
91 INTEGER :: IELEM, ELEMID
92 INTEGER :: KFACE
94 my_real :: sstar, pstar, sl, sr, wfac(3), vii(5), normalw
95 my_real :: fii(5), viistar(5), fiistar(5), pp(5)
96 INTEGER :: IMAT, NBMAT
97 INTEGER :: MATLAW(MULTI_FVM%NBMAT), LOCAL_MATID(MULTI_FVM%NBMAT)
98 my_real :: phase_rhoii(multi_fvm%NBMAT), phase_presii(multi_fvm%NBMAT),
99 . phase_eintii(multi_fvm%NBMAT),
100 . phase_alphaii(multi_fvm%NBMAT), phase_rhojj(multi_fvm%NBMAT),
101 . phase_presjj(multi_fvm%NBMAT), phase_eintjj(multi_fvm%NBMAT),
102 . phase_sspjj(multi_fvm%NBMAT), phase_alphajj(multi_fvm%NBMAT)
103 my_real :: dummy(6),dummy2(1), rhoii, pii, eintii, vxii, vyii, vzii, sspii, normal_velii, rhojj, sspjj,
104 . p_jj, normal_veljj, vxjj, vyjj, vzjj, velii2, alphaii, sub_rhoii, sub_rhoeintii
105 . sub_fiistar(3), alphastar, sub_rhostar, sub_pii, veljj2, sub_estar, eintjj
106 INTEGER :: IELEM_START, IELEM_END, ID_SURF
107 my_real :: tcar_p, tcar_vf,
alpha,beta,dp0,vnew,vold,pvois,pold,mach,rhoc2,roc,psurf
108 INTEGER :: VARTMP_EOS(1,128)
109 INTEGER :: NVARTMP_EOS
111 my_real :: m_in, m_out, e_in, e_out
112
113
114
115 nvartmp_eos = 128
116 vartmp_eos(1,:) = 1
117
118 m_in = zero
119 m_out = zero
120 e_in = zero
121 e_out = zero
122
123 tcar_p = ebcs%TCAR_P
124 tcar_vf = ebcs%TCAR_VF
125 id_surf = ebcs%SURF_ID
126 aburn(:) = zero
127
128 alpha = timestep/tcar_p
129 beta = timestep/
max(timestep,tcar_vf)
130 IF(tcar_vf>=ep20)beta=zero
131 IF(timestep == zero)THEN
133 beta = one
134 ENDIF
135
136 nbmat = multi_fvm%NBMAT
137 ielem_start = 1 + itask * nelem / nthread
138 ielem_end = (1 + itask) * nelem / nthread
139 dummy(1:6) = zero
140 dummy2(:) = one
141 DO ielem = ielem_start, ielem_end
142 elemid = elem_list(ielem)
143 IF (elemid <= numels) THEN
144 DO imat = 1, nbmat
145 local_matid(imat) = ipm(20 + imat, ixs(1, elemid))
146 matlaw(imat) = ipm(2, local_matid(imat))
147 ENDDO
148 ELSE IF (elemid <= numels + numelq) THEN
149 DO imat = 1, nbmat
150 local_matid(imat) = ipm(20 + imat, ixq(1, elemid))
151 matlaw(imat) = ipm(2, local_matid(imat))
152 ENDDO
153 ELSE
154 DO imat = 1, nbmat
155 local_matid(imat) = ipm(20 + imat, ixtg(1, elemid))
156 matlaw(imat) = ipm(2, local_matid(imat))
157 ENDDO
158 ENDIF
159
160
161 kface = face_list(ielem)
162
163 nx = multi_fvm%FACE_DATA%NORMAL(1, kface, elemid)
164 ny = multi_fvm%FACE_DATA%NORMAL(2, kface, elemid)
165 nz = multi_fvm%FACE_DATA%NORMAL(3, kface, elemid)
166 surf = multi_fvm%FACE_DATA%SURF(kface, elemid)
167 wfac(1:3) = multi_fvm%FACE_DATA%WFAC(1:3, kface, elemid)
168
169 normalw = wfac(1) * nx + wfac(2) * ny + wfac(3) * nz
170
171
172 rhoii = multi_fvm%RHO(elemid)
173 pii = multi_fvm%PRES(elemid)
174 eintii = multi_fvm%EINT(elemid)
175 vxii = multi_fvm%VEL(1, elemid)
176 vyii = multi_fvm%VEL(2, elemid)
177 vzii = multi_fvm%VEL(3, elemid)
178 velii2 = vxii**2 + vyii**2 + vzii**2
179 IF (nbmat > 1) THEN
180 DO imat = 1, nbmat
181 phase_alphaii(imat) = multi_fvm%PHASE_ALPHA(imat, elemid)
182 phase_eintii(imat) = multi_fvm%PHASE_EINT(imat, elemid)
183 phase_rhoii(imat) = multi_fvm%PHASE_RHO(imat, elemid)
184 phase_presii(imat) = multi_fvm%PHASE_PRES(imat, elemid)
185 ENDDO
186 ELSE
187 phase_alphaii(1) = one
188 phase_eintii(1) = eintii
189 phase_rhoii(1) = rhoii
190 phase_presii(1) = pii
191 ENDIF
192 sspii = multi_fvm%SOUND_SPEED(elemid)
193 normal_velii = vxii * nx + vyii * ny + vzii * nz
194
195
196 rhojj = zero
197 DO imat = 1, nbmat
198 phase_rhojj(imat) = fvm_inlet_data%VAL_RHO(imat)
199 phase_rhojj(imat) = phase_rhoii(imat)
200 phase_alphajj(imat) = fvm_inlet_data%VAL_ALPHA(imat)
201 phase_alphajj(imat) = phase_alphaii(imat)
202 rhojj = rhojj + phase_rhojj(imat) * phase_alphajj(imat)
203 ENDDO
204
205
206 sspjj = zero
207 p_jj = zero
208 eintjj = zero
209 DO imat = 1, nbmat
210 phase_eintjj(imat) = fvm_inlet_data%VAL_PRES(imat)
211 phase_eintjj(imat) = phase_eintii(imat)
212 IF (phase_alphajj(imat) > zero) THEN
214 1 0, matlaw(imat), local_matid(imat), 1,
215 2 phase_eintjj(imat), phase_presjj(imat), phase_rhojj(imat), phase_sspjj(imat),
216 3 dummy2, dummy, pm, ipm,
217 4 npropm, npropmi, dummy, dummy2,
218 5 dummy, multi_fvm%BFRAC(imat,elemid), multi_fvm%TBURN(elemid), dummy,
219 6 dummy(1), dummy, snpc, stf,
220 7 npf, tf, dummy, 1,
221 7 matparam(local_matid(imat)), nvartmp_eos, vartmp_eos, nummat, aburn)
222 sspjj = sspjj + phase_alphajj(imat) * phase_rhojj(imat) *
max(em20, phase_sspjj(imat))
223 p_jj = p_jj + phase_presjj(imat) * phase_alphajj(imat)
224 eintjj = eintjj + phase_alphajj(imat) * phase_eintjj(imat)
225 ENDIF
226 ENDDO
227
228 vnew = normal_velii
229 pvois = p_jj
230 mach = abs(normal_velii / sspii)
231 pold = ebcs%pold(ielem)
232 dp0 = ebcs%DP0(ielem)
233 IF(ncycle == 1)pold=pvois+dp0
234 vold = ebcs%vold(ielem)
235 IF(ncycle == 1) vold=normal_velii
236 ebcs%vold(ielem) = vnew
237
238 IF(mach >= one .AND. vnew > zero)THEN
239
240 pp = pvois
241 ELSE
242 rhoc2 = rhojj*sspii*sspii
243 roc = sqrt(rhojj*rhoc2)
244 p_jj = one/(one+
alpha)*(pold+roc*(vnew-vold))+
alpha*(pvois+dp0)/(
alpha + one)
245 ENDIF
246
247 ebcs%pold(ielem) = p_jj
248
249 IF (sspjj / rhojj > zero) THEN
250 sspjj = sqrt(sspjj / rhojj)
251 ELSE
252 sspjj = multi_fvm%SOUND_SPEED(elemid)
253 ENDIF
254
255
256 vxjj = fvm_inlet_data%VAL_VEL(1)
257 vxjj = vxii
258 vyjj = fvm_inlet_data%VAL_VEL(2)
259 vyjj = vyii
260 vzjj = fvm_inlet_data%VAL_VEL(3)
261 vzjj = vzii
262 normal_veljj = vxjj * nx + vyjj * ny + vzjj * nz
263
264 veljj2 = vxjj**2 + vyjj**2 + vzjj**2
265
266 sl =
min(normal_velii - sspii, normal_veljj - sspjj)
267 sr =
max(normal_velii + sspii, normal_veljj + sspjj)
268
269
270 sstar = p_jj - pii + rhoii * normal_velii * (sl - normal_velii) - rhojj * normal_veljj * (sr - normal_veljj)
271 sstar = sstar / (rhoii * (sl - normal_velii) - rhojj * (sr - normal_veljj))
272
273 pstar = pii + rhoii * (sstar - normal_velii) * (sl - normal_velii)
274 pp(1) = zero
275 pp(2) = pstar * nx
276 pp(3) = pstar * ny
277 pp(4) = pstar * nz
278 pp(5) = sstar * pstar
279
280 IF (sl > normalw) THEN
281 vii(1) = rhoii
282 vii(2) = rhoii * vxii
283 vii(3) = rhoii * vyii
284 vii(4) = rhoii * vzii
285 vii(5) = eintii + half * rhoii * velii2
286
287 fii(1) = vii(1) * normal_velii
288 fii(2) = vii(2) * normal_velii + pii * nx
289 fii(3) = vii(3) * normal_velii + pii * ny
290 fii(4) = vii(4) * normal_velii + pii * nz
291 fii(5) = (vii(5) + pii) * normal_velii
292
293 psurf = pii
294
295
296
297 multi_fvm%FLUXES(1:5, kface, elemid) = (fii(1:5) - normalw * vii(1:5)) * surf
298 multi_fvm%FLUXES(6, kface, elemid) = normal_velii * surf
299
300
301 IF (nbmat > 1) THEN
302 DO imat = 1, nbmat
303 alphaii = phase_alphaii(imat)
304 sub_rhoii = phase_rhoii(imat)
305 sub_rhoeintii = phase_eintii(imat)
306 sub_viistar(1) = alphaii
307 sub_viistar(2) = alphaii * sub_rhoii
308 sub_viistar(3) = alphaii * sub_rhoeintii
309 sub_fiistar(1:3) = sub_viistar(1:3) * normal_velii
310 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
311 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
312 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
313 ENDDO
314 ENDIF
315
316 ELSEIF (sl <= normalw .AND. normalw <= sstar) THEN
317 vii(1) = rhoii
318 vii(2) = rhoii * vxii
319 vii(3) = rhoii * vyii
320 vii(4) = rhoii * vzii
321 vii(5) = eintii + half * rhoii * velii2
322
323 fii(1) = vii(1) * normal_velii
324 fii(2) = vii(2) * normal_velii + pii * nx
325 fii(3) = vii(3) * normal_velii + pii * ny
326 fii(4) = vii(4) * normal_velii + pii * nz
327 fii(5) = (vii(5) + pii) * normal_velii
328
329 psurf = pii
330
331
332
333 viistar(1:5) = fii(1:5) - (sl) * vii(1:5) - pp(1:5)
334 viistar(1:5) = viistar(1:5) / (sstar - sl)
335 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
336 multi_fvm%FLUXES(1:5, kface, elemid) = (fiistar(1:5) - normalw * viistar(1:5)) * surf
337 multi_fvm%FLUXES(6, kface, elemid) = sstar * surf
338
339
340 IF (nbmat > 1) THEN
341 DO imat = 1, nbmat
342 matlaw(imat) = ipm(2, local_matid(imat))
343 alphastar = phase_alphaii(imat)
344 sub_rhostar = phase_rhoii(imat) * (normal_velii - sl) / (sstar - sl)
345 IF (alphastar > zero) THEN
346 sub_rhoii = phase_rhoii(imat)
347 sub_rhoeintii = phase_eintii(imat)
348 sub_pii = phase_presii(imat)
349 sub_estar = phase_eintii(imat) / phase_rhoii(imat) -
350 . phase_presii(imat) * (one / sub_rhostar - one / phase_rhoii(imat))
351
352 IF (sub_estar < zero) THEN
353 sub_estar = zero
354 ENDIF
355 ELSE
356 sub_estar = zero
357 ENDIF
358 sub_viistar(1) = alphastar
359 sub_viistar(2) = alphastar * sub_rhostar
360 sub_viistar(3) = alphastar * sub_rhostar * sub_estar
361
362 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
363 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
364 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
365 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
366 ENDDO
367 ENDIF
368
369 ELSE IF (sstar < normalw .AND. normalw <= sr) THEN
370 vii(1) = rhojj
371 vii(2) = rhojj * vxjj
372 vii(3) = rhojj * vyjj
373 vii(4) = rhojj * vzjj
374 vii(5) = eintjj + half * rhojj * veljj2
375
376 fii(1) = vii(1) * normal_veljj
377 fii(2) = vii(2) * normal_veljj + p_jj * nx
378 fii(3) = vii(3) * normal_veljj + p_jj * ny
379 fii(4) = vii(4) * normal_veljj + p_jj * nz
380 fii(5) = (vii(5) + p_jj) * normal_veljj
381
382 psurf = p_jj
383
384
385
386 viistar(1:5) = fii(1:5) - (sr) * vii(1:5) - pp(1:5)
387 viistar(1:5) = viistar(1:5) / (sstar - sr)
388 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
389 multi_fvm%FLUXES(1:5, kface, elemid) = (fiistar(1:5) - normalw * viistar(1:5)) * surf
390 multi_fvm%FLUXES(6, kface, elemid) = sstar * surf
391
392
393 IF (nbmat > 1) THEN
394 DO imat = 1, nbmat
395 matlaw(imat) = ipm(2, local_matid(imat))
396 alphastar = phase_alphajj(imat)
397 sub_rhostar = phase_rhojj(imat) * (normal_veljj - sr) / (sstar - sr)
398 IF (alphastar > zero) THEN
399 sub_rhoii = phase_rhojj(imat)
400 sub_rhoeintii = phase_eintjj(imat)
401 sub_pii = phase_presjj(imat)
402 sub_estar = phase_eintjj(imat) / phase_rhojj(imat) -
403 . phase_presjj(imat) * (one / sub_rhostar - one / phase_rhojj(imat))
404 IF (sub_estar < zero) THEN
405 sub_estar = zero
406 ENDIF
407 ELSE
408 sub_estar = zero
409 ENDIF
410
411 sub_viistar(1) = phase_alphajj(imat)
412 sub_viistar(2) = phase_alphajj(imat) * sub_rhostar
413 sub_viistar(3) = phase_alphajj(imat) * sub_rhostar * sub_estar
414 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
415 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
416 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
417 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
418 ENDDO
419 ENDIF
420 ELSE IF (sr < normalw) THEN
421 vii(1) = rhojj
422 vii(2) = rhojj * vxjj
423 vii(3) = rhojj * vyjj
424 vii(4) = rhojj * vzjj
425 vii(5) = eintjj + half * rhojj * veljj2
426
427 fii(1) = vii(1) * normal_veljj
428 fii(2) = vii(2) * normal_veljj + p_jj * nx
429 fii(3) = vii(3) * normal_veljj + p_jj * ny
430 fii(4) = vii(4) * normal_veljj + p_jj * nz
431 fii(5) = (vii(5) + p_jj) * normal_veljj
432
433 psurf = p_jj
434
435
436
437 multi_fvm%FLUXES(1:5, kface, elemid) = (fii(1:5) - normalw * vii(1:5)) * surf
438 multi_fvm%FLUXES(6, kface, elemid) = normal_veljj * surf
439
440
441 IF (nbmat > 1) THEN
442 DO imat = 1, nbmat
443 alphaii = phase_alphajj(imat)
444 sub_rhoii = phase_rhojj(imat)
445 sub_rhoeintii = phase_eintjj(imat)
446 sub_viistar(1) = alphaii
447 sub_viistar(2) = alphaii * sub_rhoii
448 sub_viistar(3) = alphaii * sub_rhoeintii
449 sub_fiistar(1:3) = sub_viistar(1:3) * normal_veljj
450 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
451 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
452 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
453 ENDDO
454 ENDIF
455
456 ELSE
458
459 ENDIF
460
461 IF(.NOT. pred)THEN
462
463
464
465 fsavsurf(2,id_surf) = fsavsurf(2,id_surf) + multi_fvm%FLUXES(1, kface, elemid)
466 fsavsurf(3,id_surf) = fsavsurf(3,id_surf) + multi_fvm%FLUXES(1, kface, elemid) / vii(1)
467 fsavsurf(4,id_surf) = fsavsurf(4,id_surf) + surf*psurf
468 fsavsurf(6,id_surf) = fsavsurf(6,id_surf) + multi_fvm%FLUXES(1, kface, elemid)
469
470
471 m_in = m_in - timestep *
min(zero, multi_fvm%FLUXES(1, kface, elemid) )
472 m_out = m_out - timestep *
max(zero, multi_fvm%FLUXES(1, kface, elemid) )
473 e_in = e_in - timestep *
min(zero, multi_fvm%FLUXES(5, kface, elemid) )
474 e_out = e_out - timestep *
max(zero, multi_fvm%FLUXES(5, kface, elemid) )
475 ENDIF
476 ENDDO
477
478 IF(.NOT. pred)THEN
479
480 output%DATA%INOUT%DM_IN = output%DATA%INOUT%DM_IN + m_in
481 output%DATA%INOUT%DM_OUT = output%DATA%INOUT%DM_OUT + m_out
482 output%DATA%INOUT%DE_IN = output%DATA%INOUT%DE_IN + e_in
483 output%DATA%INOUT%DE_OUT = output%DATA%INOUT%DE_OUT + e_out
484
485 ENDIF
486
487
488
subroutine multi_submatlaw(iflag, matlaw, local_matid, nel, eint, pres, rho, ssp, vol, grun, pm, ipm, npropm, npropmi, bufmat, off, theta, burnfrac, burntime, deltax, current_time, sigold, snpc, stf, npf, tf, vareos, nvareos, mat_param, nvartmp_eos, vartmp_eos, nummat, aburn)
OPTION /TH/SURF outputs of Pressure and Area needed Tabs.
integer, parameter th_surf_num_channel
number of /TH/SURF channels : AREA, VELOCITY, MASSFLOW, P A, MASS