OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
multi_inlet_ebcs.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| multi_inlet_ebcs_mod ../engine/source/multifluid/multi_inlet_ebcs.F
25!||--- called by ------------------------------------------------------
26!|| multi_ebcs ../engine/source/multifluid/multi_ebcs.F
27!||====================================================================
28 MODULE multi_inlet_ebcs_mod
29 CONTAINS
30!||====================================================================
31!|| multi_inlet_ebcs ../engine/source/multifluid/multi_inlet_ebcs.F
32!||--- called by ------------------------------------------------------
33!|| multi_ebcs ../engine/source/multifluid/multi_ebcs.F
34!||--- calls -----------------------------------------------------
35!|| arret ../engine/source/system/arret.F
36!|| multi_solve_eint ../engine/source/multifluid/multi_solve_eint.F90
37!|| multi_submatlaw ../engine/source/multifluid/multi_submatlaw.F
38!||--- uses -----------------------------------------------------
39!|| element_mod ../common_source/modules/elements/element_mod.F90
40!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
41!|| multi_fvm_mod ../common_source/modules/ale/multi_fvm_mod.F90
42!|| multi_solve_eint_mod ../engine/source/multifluid/multi_solve_eint.F90
43!|| multi_submatlaw_mod ../engine/source/multifluid/multi_submatlaw.F
44!|| output_mod ../common_source/modules/output/output_mod.F90
45!|| th_surf_mod ../common_source/modules/interfaces/th_surf_mod.F
46!||====================================================================
47 SUBROUTINE multi_inlet_ebcs(ITASK, EBCS_ID, MULTI_FVM, NELEM, ELEM_LIST, FACE_LIST,
48 . FVM_INLET_DATA, IXS, IXQ, IXTG, XGRID, WGRID, IPM, PM,
49 . FUNC_VALUE, ID_SURF, NPF,TF,FSAVSURF, TIMESTEP, MATPARAM, OUTPUT, PRED)
50C-----------------------------------------------
51C M o d u l e s
52C-----------------------------------------------
53 USE multi_fvm_mod
54 USE th_surf_mod , only : th_surf_num_channel
55 USE matparam_def_mod , ONLY : matparam_struct_
56 USE multi_solve_eint_mod , ONLY : multi_solve_eint
57 USE multi_submatlaw_mod , ONLY : multi_submatlaw
58 USE output_mod , ONLY : output_
59 use element_mod , only : nixs,nixq,nixtg
60C-----------------------------------------------
61C I m p l i c i t T y p e s
62C-----------------------------------------------
63#include "implicit_f.inc"
64C-----------------------------------------------
65C C o m m o n B l o c k s
66C-----------------------------------------------
67! NIXS
68#include "param_c.inc"
69! ISPMD
70#include "task_c.inc"
71! ALE
72! NUMELS, NUMELQ, NUMELTG
73#include "com04_c.inc"
74! SNPC,STF
75 COMMON /tablesizf/ stf,snpc
76 INTEGER STF,SNPC
77C-----------------------------------------------
78C D u m m y A r g u m e n t s
79C-----------------------------------------------
80 INTEGER, INTENT(IN) :: EBCS_ID
81 TYPE(MULTI_FVM_STRUCT), INTENT(INOUT) :: MULTI_FVM
82 INTEGER, INTENT(IN) :: ITASK, NELEM, ELEM_LIST(NELEM), FACE_LIST(NELEM)
83 INTEGER, INTENT(IN) :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
84 my_real, INTENT(IN) :: xgrid(3, *), wgrid(3, *)
85 INTEGER, INTENT(IN) :: IPM(NPROPMI, NUMMAT)
86 my_real, INTENT(IN) :: pm(npropm, nummat), func_value(*)
87 TYPE(fvm_inlet_data_struct), INTENT(IN) :: FVM_INLET_DATA
88 INTEGER, INTENT(IN) :: NPF(SNPC), ID_SURF
89 my_real, INTENT(IN) :: tf(stf), timestep
90 my_real, INTENT(INOUT) :: fsavsurf(th_surf_num_channel,nsurf)
91 TYPE(matparam_struct_), DIMENSION(NUMMAT),INTENT(IN) :: MATPARAM !material data structure
92 TYPE(output_), INTENT(INOUT) :: OUTPUT !< output structure
93 LOGICAL,INTENT(IN) :: PRED
94C-----------------------------------------------
95C L o c a l V a r i a b l e s
96C-----------------------------------------------
97 INTEGER :: IELEM, ELEMID
98 INTEGER :: KFACE
99 my_real :: surf, nx, ny, nz
100 my_real :: sstar, pstar, sl, sr, wfac(3), vii(5), normalw
101 my_real :: fii(5), viistar(5), fiistar(5), pp(5)
102 INTEGER :: IMAT, NBMAT
103 INTEGER :: MATLAW(MULTI_FVM%NBMAT), LOCAL_MATID(MULTI_FVM%NBMAT)
104 my_real :: phase_rhoii(multi_fvm%NBMAT), phase_presii(multi_fvm%NBMAT),
105 . phase_eintii(multi_fvm%NBMAT),
106 . phase_alphaii(multi_fvm%NBMAT), phase_rhojj(multi_fvm%NBMAT),
107 . phase_presjj(multi_fvm%NBMAT), phase_eintjj(multi_fvm%NBMAT),
108 . phase_sspjj(multi_fvm%NBMAT), phase_alphajj(multi_fvm%NBMAT)
109 my_real :: dummy(6), dummy2(1), rhoii, pii, eintii, vxii, vyii, vzii, sspii, normal_velii, rhojj, sspjj,
110 . pjj, normal_veljj, vxjj, vyjj, vzjj, velii2, alphaii, sub_rhoii, sub_rhoeintii, sub_viistar(3),
111 . sub_fiistar(3), alphastar, sub_rhostar, sub_pii, veljj2, sub_estar, eintjj, psurf
112 INTEGER :: IFUNC, IELEM_START, IELEM_END
113 my_real :: ABURN(1) !< afterburning
114 INTEGER :: VARTMP_EOS(1,128) !nel=1
115 INTEGER :: NVARTMP_EOS !upper bound
116 my_real :: m_in, m_out, e_in, e_out
117C-----------------------------------------------
118C B o d y
119C-----------------------------------------------
120 nvartmp_eos = 128
121 vartmp_eos(1,:) = 1
122 nbmat = multi_fvm%NBMAT
123 ielem_start = 1 + itask * nelem / nthread
124 ielem_end = (1 + itask) * nelem / nthread
125 dummy(1:6) = zero
126 dummy2(1) = one
127 aburn(:) = zero
128 m_in = zero
129 m_out = zero
130 e_in = zero
131 e_out = zero
132 DO ielem = ielem_start, ielem_end
133 elemid = elem_list(ielem)
134 IF (elemid <= numels) THEN
135 DO imat = 1, nbmat
136 local_matid(imat) = ipm(20 + imat, ixs(1, elemid))
137 matlaw(imat) = ipm(2, local_matid(imat))
138 ENDDO
139 ELSE IF (elemid <= numels + numelq) THEN
140 DO imat = 1, nbmat
141 local_matid(imat) = ipm(20 + imat, ixq(1, elemid))
142 matlaw(imat) = ipm(2, local_matid(imat))
143 ENDDO
144 ELSE
145 DO imat = 1, nbmat
146 local_matid(imat) = ipm(20 + imat, ixtg(1, elemid))
147 matlaw(imat) = ipm(2, local_matid(imat))
148 ENDDO
149 ENDIF
150
151C Face of the element
152 kface = face_list(ielem)
153
154 nx = multi_fvm%FACE_DATA%NORMAL(1, kface, elemid)
155 ny = multi_fvm%FACE_DATA%NORMAL(2, kface, elemid)
156 nz = multi_fvm%FACE_DATA%NORMAL(3, kface, elemid)
157 surf = multi_fvm%FACE_DATA%SURF(kface, elemid)
158 wfac(1:3) = multi_fvm%FACE_DATA%WFAC(1:3, kface, elemid)
159C Normal grid velocity
160 normalw = wfac(1) * nx + wfac(2) * ny + wfac(3) * nz
161
162C Current element
163 rhoii = multi_fvm%RHO(elemid)
164 pii = multi_fvm%PRES(elemid)
165 eintii = multi_fvm%EINT(elemid)
166 vxii = multi_fvm%VEL(1, elemid)
167 vyii = multi_fvm%VEL(2, elemid)
168 vzii = multi_fvm%VEL(3, elemid)
169 velii2 = vxii**2 + vyii**2 + vzii**2
170 IF (nbmat > 1) THEN
171 DO imat = 1, nbmat
172 phase_alphaii(imat) = multi_fvm%PHASE_ALPHA(imat, elemid)
173 phase_eintii(imat) = multi_fvm%PHASE_EINT(imat, elemid)
174 phase_rhoii(imat) = multi_fvm%PHASE_RHO(imat, elemid)
175 phase_presii(imat) = multi_fvm%PHASE_PRES(imat, elemid)
176 ENDDO
177 ELSE
178 phase_alphaii(1) = one
179 phase_eintii(1) = eintii
180 phase_rhoii(1) = rhoii
181 phase_presii(1) = pii
182 ENDIF
183 sspii = multi_fvm%SOUND_SPEED(elemid)
184 normal_velii = vxii * nx + vyii * ny + vzii * nz
185
186C Boundary "GHOST" element
187 rhojj = zero
188 DO imat = 1, nbmat
189 phase_rhojj(imat) = fvm_inlet_data%VAL_RHO(imat)
190 ifunc = fvm_inlet_data%FUNC_RHO(imat)
191 IF (ifunc > 0) THEN
192 phase_rhojj(imat) = phase_rhojj(imat) * func_value(ifunc)
193 ELSEIF (ifunc == 0) THEN
194 phase_rhojj(imat) = phase_rhojj(imat)
195 ELSEIF (ifunc == -1) THEN
196 phase_rhojj(imat) = phase_rhoii(imat)
197 ENDIF
198 phase_alphajj(imat) = fvm_inlet_data%VAL_ALPHA(imat)
199 ifunc = fvm_inlet_data%FUNC_ALPHA(imat)
200 IF (ifunc > 0) THEN
201 phase_alphajj(imat) = phase_alphajj(imat) * func_value(ifunc)
202 ELSEIF (ifunc == 0) THEN
203 phase_alphajj(imat) = phase_alphajj(imat)
204 ELSEIF (ifunc == -1) THEN
205 phase_alphajj(imat) = phase_alphaii(imat)
206 ENDIF
207 rhojj = rhojj + phase_rhojj(imat) * phase_alphajj(imat)
208 ENDDO
209 IF (fvm_inlet_data%FORMULATION == 2) THEN
210C VE formulation
211 sspjj = zero
212 pjj = zero
213 eintjj = zero
214 DO imat = 1, nbmat
215 phase_eintjj(imat) = fvm_inlet_data%VAL_PRES(imat)
216 ifunc = fvm_inlet_data%FUNC_PRES(imat)
217 IF (ifunc > 0) THEN
218 phase_eintjj(imat) = phase_eintjj(imat) * func_value(ifunc)
219 ELSEIF (ifunc == 0) THEN
220 phase_eintjj(imat) = phase_eintjj(imat)
221 ELSEIF (ifunc == -1) THEN
222 phase_eintjj(imat) = phase_eintii(imat)
223 ENDIF
224 IF (phase_alphajj(imat) > zero) THEN
225 CALL multi_submatlaw(
226 1 0, matlaw(imat), local_matid(imat), 1,
227 2 phase_eintjj(imat), phase_presjj(imat), phase_rhojj(imat), phase_sspjj(imat),
228 3 dummy2, dummy, pm, ipm,
229 4 npropm, npropmi, dummy, dummy2,
230 5 dummy, multi_fvm%BFRAC(imat,elemid), multi_fvm%TBURN(elemid), dummy ,
231 6 dummy(1), dummy, snpc , stf ,
232 7 npf, tf, dummy(1) , 1 ,
233 8 matparam(local_matid(imat)),nvartmp_eos, vartmp_eos ,nummat ,
234 9 aburn)
235 sspjj = sspjj + phase_alphajj(imat) * phase_rhojj(imat) * max(em20, phase_sspjj(imat))
236 pjj = pjj + phase_presjj(imat) * phase_alphajj(imat)
237 eintjj = eintjj + phase_alphajj(imat) * phase_eintjj(imat)
238 ENDIF
239 ENDDO
240
241 ELSEIF (fvm_inlet_data%FORMULATION == 1) THEN
242C VP formulation
243 sspjj = zero
244 pjj = zero
245 eintjj = zero
246 DO imat = 1, nbmat
247 phase_presjj(imat) = fvm_inlet_data%VAL_PRES(imat)
248 ifunc = fvm_inlet_data%FUNC_PRES(imat)
249 IF (ifunc > 0) THEN
250 phase_presjj(imat) = phase_presjj(imat) * func_value(ifunc)
251 ELSEIF (ifunc == 0) THEN
252 phase_presjj(imat) = phase_presjj(imat)
253 ELSEIF (ifunc == -1) THEN
254 phase_presjj(imat) = phase_presii(imat)
255 ENDIF
256 pjj = pjj + phase_presjj(imat) * phase_alphajj(imat)
257C Need to solve p(RHO, EINT) = PHASE_PRESJJ(IMAT) for EINT
258C INITIALIZE INTERNAL ENERGY
259 phase_eintjj(imat) = pm(23, local_matid(imat))
260 CALL multi_solve_eint(matlaw(imat), local_matid(imat), pm, ipm, npropm, npropmi,
261 . phase_eintjj(imat), phase_rhojj(imat), phase_presjj(imat), phase_sspjj(imat),
262 . multi_fvm%BFRAC(imat, elemid), multi_fvm%TBURN(elemid), dummy(1), dummy(1),
263 . dummy, dummy2, snpc, stf, npf, tf, dummy, 1, matparam(local_matid(imat)),
264 ? nvartmp_eos, vartmp_eos, nummat, aburn)
265
266 sspjj = sspjj + phase_alphajj(imat) * phase_rhojj(imat) * max(em20, phase_sspjj(imat))
267 eintjj = eintjj + phase_alphajj(imat) * phase_eintjj(imat)
268 ENDDO
269 ENDIF
270 IF (sspjj / rhojj > zero) THEN
271 sspjj = sqrt(sspjj / rhojj)
272 ELSE
273 sspjj = multi_fvm%SOUND_SPEED(elemid)
274 ENDIF
275 IF (fvm_inlet_data%VECTOR_VELOCITY == 0) THEN
276C Normal velocity imposed
277 normal_veljj = -fvm_inlet_data%VAL_VEL(1)
278 ifunc = fvm_inlet_data%FUNC_VEL(1)
279 IF (ifunc > 0) THEN
280 normal_veljj = normal_veljj * func_value(ifunc)
281 ELSEIF (ifunc == 0) THEN
282 normal_veljj = normal_veljj
283 ELSEIF (ifunc == -1) THEN
284 normal_veljj = normal_velii
285 ENDIF
286C NORMAL_VELJJ = NORMAL_VELII
287
288 vxjj = normal_veljj * nx
289 vyjj = normal_veljj * ny
290 vzjj = normal_veljj * nz
291 ELSE
292 vxjj = fvm_inlet_data%VAL_VEL(1)
293 ifunc = fvm_inlet_data%FUNC_VEL(1)
294 IF (ifunc > 0) THEN
295 vxjj = vxjj * func_value(ifunc)
296 ELSEIF (ifunc == -1) THEN
297 vxjj = vxii
298 ENDIF
299 vyjj = fvm_inlet_data%VAL_VEL(2)
300 ifunc = fvm_inlet_data%FUNC_VEL(2)
301 IF (ifunc > 0) THEN
302 vyjj = vyjj * func_value(ifunc)
303 ELSEIF (ifunc == -1) THEN
304 vyjj = vyii
305 ENDIF
306 vzjj = fvm_inlet_data%VAL_VEL(3)
307 ifunc = fvm_inlet_data%FUNC_VEL(3)
308 IF (ifunc > 0) THEN
309 vzjj = vzjj * func_value(ifunc)
310 ELSEIF (ifunc == -1) THEN
311 vzjj = vzii
312 ENDIF
313 normal_veljj = vxjj * nx + vyjj * ny + vzjj * nz
314 ENDIF
315 veljj2 = vxjj**2 + vyjj**2 + vzjj**2
316C HLL wave speed estimates
317 sl = min(normal_velii - sspii, normal_veljj - sspjj)
318 sr = max(normal_velii + sspii, normal_veljj + sspjj)
319
320C Intermediate wave speed
321 sstar = pjj - pii + rhoii * normal_velii * (sl - normal_velii) - rhojj * normal_veljj * (sr - normal_veljj)
322 sstar = sstar / (rhoii * (sl - normal_velii) - rhojj * (sr - normal_veljj))
323
324 pstar = pii + rhoii * (sstar - normal_velii) * (sl - normal_velii)
325 pp(1) = zero
326 pp(2) = pstar * nx
327 pp(3) = pstar * ny
328 pp(4) = pstar * nz
329 pp(5) = sstar * pstar
330
331 IF (sl > normalw) THEN
332 vii(1) = rhoii
333 vii(2) = rhoii * vxii
334 vii(3) = rhoii * vyii
335 vii(4) = rhoii * vzii
336 vii(5) = eintii + half * rhoii * velii2
337! Normal physical flux current element
338 fii(1) = vii(1) * normal_velii
339 fii(2) = vii(2) * normal_velii + pii * nx
340 fii(3) = vii(3) * normal_velii + pii * ny
341 fii(4) = vii(4) * normal_velii + pii * nz
342 fii(5) = (vii(5) + pii) * normal_velii
343! /th/surf : pressure
344 psurf = pii
345C Take the fluxes of cell II
346C ===
347C Global fluxes
348 multi_fvm%FLUXES(1:5, kface, elemid) = (fii(1:5) - normalw * vii(1:5)) * surf
349 multi_fvm%FLUXES(6, kface, elemid) = normal_velii * surf
350C ===
351C Submaterial fluxes
352 IF (nbmat > 1) THEN
353 DO imat = 1, nbmat
354 alphaii = phase_alphaii(imat)
355 sub_rhoii = phase_rhoii(imat) ! ALPHA_RHO
356 sub_rhoeintii = phase_eintii(imat)
357 sub_viistar(1) = alphaii
358 sub_viistar(2) = alphaii * sub_rhoii ! ALPHA_RHO
359 sub_viistar(3) = alphaii * sub_rhoeintii
360 sub_fiistar(1:3) = sub_viistar(1:3) * normal_velii
361 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
362 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
363 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
364 ENDDO
365 ENDIF
366C
367 ELSEIF (sl <= normalw .AND. normalw <= sstar) THEN
368 vii(1) = rhoii
369 vii(2) = rhoii * vxii
370 vii(3) = rhoii * vyii
371 vii(4) = rhoii * vzii
372 vii(5) = eintii + half * rhoii * velii2
373! Normal physical flux current element
374 fii(1) = vii(1) * normal_velii
375 fii(2) = vii(2) * normal_velii + pii * nx
376 fii(3) = vii(3) * normal_velii + pii * ny
377 fii(4) = vii(4) * normal_velii + pii * nz
378 fii(5) = (vii(5) + pii) * normal_velii
379! /th/surf : pressure
380 psurf = pii
381C Take intermediate state flux (HLLC scheme)
382C ===
383C Global fluxes
384 viistar(1:5) = fii(1:5) - (sl) * vii(1:5) - pp(1:5)
385 viistar(1:5) = viistar(1:5) / (sstar - sl)
386 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
387 multi_fvm%FLUXES(1:5, kface, elemid) = (fiistar(1:5) - normalw * viistar(1:5)) * surf
388 multi_fvm%FLUXES(6, kface, elemid) = sstar * surf
389C ===
390C Submaterial fluxes
391 IF (nbmat > 1) THEN
392 DO imat = 1, nbmat
393 matlaw(imat) = ipm(2, local_matid(imat))
394 alphastar = phase_alphaii(imat)
395 sub_rhostar = phase_rhoii(imat) * (normal_velii - sl) / (sstar - sl)
396 IF (alphastar > zero) THEN
397 sub_rhoii = phase_rhoii(imat)
398 sub_rhoeintii = phase_eintii(imat)
399 sub_pii = phase_presii(imat)
400 sub_estar = phase_eintii(imat) / phase_rhoii(imat) -
401 . phase_presii(imat) * (one / sub_rhostar - one / phase_rhoii(imat))
402
403 IF (sub_estar < zero) THEN
404 sub_estar = zero
405 ENDIF
406 ELSE
407 sub_estar = zero
408 ENDIF
409 sub_viistar(1) = alphastar
410 sub_viistar(2) = alphastar * sub_rhostar
411 sub_viistar(3) = alphastar * sub_rhostar * sub_estar
412
413 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
414 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
415 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
416 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
417 ENDDO
418 ENDIF
419
420 ELSE IF (sstar < normalw .AND. normalw <= sr) 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! Normal physical flux current element
427 fii(1) = vii(1) * normal_veljj
428 fii(2) = vii(2) * normal_veljj + pjj * nx
429 fii(3) = vii(3) * normal_veljj + pjj * ny
430 fii(4) = vii(4) * normal_veljj + pjj * nz
431 fii(5) = (vii(5) + pjj) * normal_veljj
432! /th/surf : pressure
433 psurf = pjj
434C Take intermediate state flux (HLLC scheme)
435C ===
436C Global fluxes
437 viistar(1:5) = fii(1:5) - (sr) * vii(1:5) - pp(1:5)
438 viistar(1:5) = viistar(1:5) / (sstar - sr)
439 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
440 multi_fvm%FLUXES(1:5, kface, elemid) = (fiistar(1:5) - normalw * viistar(1:5)) * surf
441 multi_fvm%FLUXES(6, kface, elemid) = sstar * surf
442C ===
443C Submaterial fluxes
444 IF (nbmat > 1) THEN
445 DO imat = 1, nbmat
446 matlaw(imat) = ipm(2, local_matid(imat))
447 alphastar = phase_alphajj(imat)
448 sub_rhostar = phase_rhojj(imat) * (normal_veljj - sr) / (sstar - sr)
449 IF (alphastar > zero) THEN
450 sub_rhoii = phase_rhojj(imat)
451 sub_rhoeintii = phase_eintjj(imat)
452 sub_pii = phase_presjj(imat)
453 sub_estar = phase_eintjj(imat) / phase_rhojj(imat) -
454 . phase_presjj(imat) * (one / sub_rhostar - one / phase_rhojj(imat))
455 IF (sub_estar < zero) THEN
456 sub_estar = zero
457 ENDIF
458 ELSE
459 sub_estar = zero
460 ENDIF
461
462 sub_viistar(1) = phase_alphajj(imat)
463 sub_viistar(2) = phase_alphajj(imat) * sub_rhostar
464 sub_viistar(3) = phase_alphajj(imat) * sub_rhostar * sub_estar
465 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
466 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
467 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
468 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
469 ENDDO
470 ENDIF
471 ELSE IF (sr < normalw) THEN
472 vii(1) = rhojj
473 vii(2) = rhojj * vxjj
474 vii(3) = rhojj * vyjj
475 vii(4) = rhojj * vzjj
476 vii(5) = eintjj + half * rhojj * veljj2
477! Normal physical flux current element
478 fii(1) = vii(1) * normal_veljj
479 fii(2) = vii(2) * normal_veljj + pjj * nx
480 fii(3) = vii(3) * normal_veljj + pjj * ny
481 fii(4) = vii(4) * normal_veljj + pjj * nz
482 fii(5) = (vii(5) + pjj) * normal_veljj
483! /th/surf : pressure
484 psurf = pjj
485C Take the fluxes of cell II
486C ===
487C Global fluxes
488 multi_fvm%FLUXES(1:5, kface, elemid) = (fii(1:5) - normalw * vii(1:5)) * surf
489 multi_fvm%FLUXES(6, kface, elemid) = normal_veljj * surf
490C ===
491C Submaterial fluxes
492 IF (nbmat > 1) THEN
493 DO imat = 1, nbmat
494 alphaii = phase_alphajj(imat)
495 sub_rhoii = phase_rhojj(imat) ! ALPHA_RHO
496 sub_rhoeintii = phase_eintjj(imat)
497 sub_viistar(1) = alphaii
498 sub_viistar(2) = alphaii * sub_rhoii ! ALPHA_RHO
499 sub_viistar(3) = alphaii * sub_rhoeintii
500 sub_fiistar(1:3) = sub_viistar(1:3) * normal_veljj
501 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
502 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
503 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
504 ENDDO
505 ENDIF
506 ELSE
507 CALL arret(2)
508C
509 ENDIF
510
511
512 IF(.NOT.pred)THEN
513 !muscl 'off', or muscl 'on' after 2nd step of prediction/correction
514
515 ! /TH/SURF (massflow,velocity,pressure,acceleration)
516 fsavsurf(2,id_surf) = fsavsurf(2,id_surf) + multi_fvm%FLUXES(1, kface, elemid) ! += rho.S.un
517 fsavsurf(3,id_surf) = fsavsurf(3,id_surf) + multi_fvm%FLUXES(1, kface, elemid) / vii(1) ! += S.un
518 fsavsurf(4,id_surf) = fsavsurf(4,id_surf) + surf*psurf ! += S.P
519 fsavsurf(6,id_surf) = fsavsurf(6,id_surf) + multi_fvm%FLUXES(1, kface, elemid) * timestep ! m <- m+dm (cumulative value)
520
521 ! FLUX > 0 if flow is exiting (OUTWARD NORMAL)
522 m_in = m_in - timestep * min(zero, multi_fvm%FLUXES(1, kface, elemid) )
523 m_out = m_out - timestep * max(zero, multi_fvm%FLUXES(1, kface, elemid) )
524 e_in = e_in - timestep * min(zero, multi_fvm%FLUXES(5, kface, elemid) )
525 e_out = e_out - timestep * max(zero, multi_fvm%FLUXES(5, kface, elemid) )
526 ENDIF
527 ENDDO
528
529 IF(.NOT. pred)THEN
530!$OMP CRITICAL
531 output%DATA%INOUT%DM_IN = output%DATA%INOUT%DM_IN + m_in
532 output%DATA%INOUT%DM_OUT = output%DATA%INOUT%DM_OUT + m_out
533 output%DATA%INOUT%DE_IN = output%DATA%INOUT%DE_IN + e_in
534 output%DATA%INOUT%DE_OUT = output%DATA%INOUT%DE_OUT + e_out
535!$OMP END CRITICAL
536 ENDIF
537
538C-----------------------------------------------
539C E n d o f s u b r o u t i n e
540C-----------------------------------------------
541 END SUBROUTINE multi_inlet_ebcs
542 END MODULE multi_inlet_ebcs_mod
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
multi_inlet_ebcs_mod
Definition multi_inlet_ebcs.F:28
multi_inlet_ebcs_mod::multi_inlet_ebcs
subroutine multi_inlet_ebcs(itask, ebcs_id, multi_fvm, nelem, elem_list, face_list, fvm_inlet_data, ixs, ixq, ixtg, xgrid, wgrid, ipm, pm, func_value, id_surf, npf, tf, fsavsurf, timestep, matparam, output, pred)
Definition multi_inlet_ebcs.F:50
multi_submatlaw_mod
Definition multi_submatlaw.F:33
multi_submatlaw_mod::multi_submatlaw
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)
Definition multi_submatlaw.F:61
th_surf_mod
OPTION /TH/SURF outputs of Pressure and Area needed Tabs.
Definition th_surf_mod.F:61
th_surf_mod::th_surf_num_channel
integer, parameter th_surf_num_channel
number of /TH/SURF channels : AREA, VELOCITY, MASSFLOW, P A, MASS
Definition th_surf_mod.F:100
arret
subroutine arret(nn)
Definition arret.F:86