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