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!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
40!|| multi_fvm_mod ../common_source/modules/ale/multi_fvm_mod.F90
41!|| multi_solve_eint_mod ../engine/source/multifluid/multi_solve_eint.F90
42!|| multi_submatlaw_mod ../engine/source/multifluid/multi_submatlaw.F
43!|| th_surf_mod ../common_source/modules/interfaces/th_surf_mod.F
44!||====================================================================
45 SUBROUTINE multi_inlet_ebcs(ITASK, EBCS_ID, MULTI_FVM, NELEM, ELEM_LIST, FACE_LIST,
46 . FVM_INLET_DATA, IXS, IXQ, IXTG, XGRID, WGRID, IPM, PM,
47 . FUNC_VALUE, ID_SURF, NPF,TF,FSAVSURF, TIMESTEP, MATPARAM)
48C-----------------------------------------------
49C M o d u l e s
50C-----------------------------------------------
51 USE multi_fvm_mod
52 USE th_surf_mod , only : th_surf_num_channel
53 USE matparam_def_mod , ONLY : matparam_struct_
54 USE multi_solve_eint_mod , ONLY : multi_solve_eint
55 USE multi_submatlaw_mod , ONLY : multi_submatlaw
56C-----------------------------------------------
57C I m p l i c i t T y p e s
58C-----------------------------------------------
59#include "implicit_f.inc"
60C-----------------------------------------------
61C C o m m o n B l o c k s
62C-----------------------------------------------
63! NIXS
64#include "param_c.inc"
65! ISPMD
66#include "task_c.inc"
67! ALE
68! NUMELS, NUMELQ, NUMELTG
69#include "com04_c.inc"
70! SNPC,STF
71 COMMON /tablesizf/ stf,snpc
72 INTEGER STF,SNPC
73C-----------------------------------------------
74C D u m m y A r g u m e n t s
75C-----------------------------------------------
76 INTEGER, INTENT(IN) :: EBCS_ID
77 TYPE(MULTI_FVM_STRUCT), INTENT(INOUT) :: MULTI_FVM
78 INTEGER, INTENT(IN) :: ITASK, NELEM, ELEM_LIST(NELEM), FACE_LIST(NELEM)
79 INTEGER, INTENT(IN) :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
80 my_real, INTENT(IN) :: xgrid(3, *), wgrid(3, *)
81 INTEGER, INTENT(IN) :: IPM(NPROPMI, NUMMAT)
82 my_real, INTENT(IN) :: pm(npropm, nummat), func_value(*)
83 TYPE(fvm_inlet_data_struct), INTENT(IN) :: FVM_INLET_DATA
84 INTEGER, INTENT(IN) :: NPF(SNPC), ID_SURF
85 my_real, INTENT(IN) :: tf(stf), timestep
86 my_real, INTENT(INOUT) :: fsavsurf(th_surf_num_channel,nsurf)
87 TYPE(matparam_struct_), DIMENSION(NUMMAT),INTENT(IN) :: MATPARAM !material data structure
88C-----------------------------------------------
89C L o c a l V a r i a b l e s
90C-----------------------------------------------
91 INTEGER :: IELEM, ELEMID, NODE_ID
92 INTEGER :: KFACE, NB_NOD_FOUND, KK, JJ
93 my_real :: x1(3), x2(3), x3(3), x4(3), p, normvel, vx, vy, vz, ssp, surf, nx, ny, nz
94 my_real :: rho, sstar, pstar, pfac, sl, sr, wfac(3), vii(5), vjj(5), vel2, normalw
95 my_real :: fii(5), fjj(5), viistar(5), vjjstar(5), fiistar(5), fjjstar(5), pp(5)
96 INTEGER :: NODE1, NODE2, 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), phase_sspii(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 . pjj, normal_veljj, vxjj, vyjj, vzjj, velii2, alphaii, sub_rhoii, sub_rhoeintii, sub_viistar(3),
105 . sub_fiistar(3), alphastar, sub_rhostar, sub_pii, veljj2, sub_estar, eintjj, psurf
106 INTEGER :: IFUNC, IELEM_START, IELEM_END
107 my_real :: ABURN(1) !< afterburning
108 INTEGER :: VARTMP_EOS(1,128) !nel=1
109 INTEGER :: NVARTMP_EOS !upper bound
110C-----------------------------------------------
111C B o d y
112C-----------------------------------------------
113 nvartmp_eos = 128
114 vartmp_eos(1,:) = 1
115 nbmat = multi_fvm%NBMAT
116 ielem_start = 1 + itask * nelem / nthread
117 ielem_end = (1 + itask) * nelem / nthread
118 dummy(1:6) = zero
119 dummy2(1) = one
120 aburn(:) = zero
121 DO ielem = ielem_start, ielem_end
122 elemid = elem_list(ielem)
123 IF (elemid <= numels) THEN
124 DO imat = 1, nbmat
125 local_matid(imat) = ipm(20 + imat, ixs(1, elemid))
126 matlaw(imat) = ipm(2, local_matid(imat))
127 ENDDO
128 ELSE IF (elemid <= numels + numelq) THEN
129 DO imat = 1, nbmat
130 local_matid(imat) = ipm(20 + imat, ixq(1, elemid))
131 matlaw(imat) = ipm(2, local_matid(imat))
132 ENDDO
133 ELSE
134 DO imat = 1, nbmat
135 local_matid(imat) = ipm(20 + imat, ixtg(1, elemid))
136 matlaw(imat) = ipm(2, local_matid(imat))
137 ENDDO
138 ENDIF
139
140C Face of the element
141 kface = face_list(ielem)
142
143 nx = multi_fvm%FACE_DATA%NORMAL(1, kface, elemid)
144 ny = multi_fvm%FACE_DATA%NORMAL(2, kface, elemid)
145 nz = multi_fvm%FACE_DATA%NORMAL(3, kface, elemid)
146 surf = multi_fvm%FACE_DATA%SURF(kface, elemid)
147 wfac(1:3) = multi_fvm%FACE_DATA%WFAC(1:3, kface, elemid)
148C Normal grid velocity
149 normalw = wfac(1) * nx + wfac(2) * ny + wfac(3) * nz
150
151C Current element
152 rhoii = multi_fvm%RHO(elemid)
153 pii = multi_fvm%PRES(elemid)
154 eintii = multi_fvm%EINT(elemid)
155 vxii = multi_fvm%VEL(1, elemid)
156 vyii = multi_fvm%VEL(2, elemid)
157 vzii = multi_fvm%VEL(3, elemid)
158 velii2 = vxii**2 + vyii**2 + vzii**2
159 IF (nbmat > 1) THEN
160 DO imat = 1, nbmat
161 phase_alphaii(imat) = multi_fvm%PHASE_ALPHA(imat, elemid)
162 phase_eintii(imat) = multi_fvm%PHASE_EINT(imat, elemid)
163 phase_rhoii(imat) = multi_fvm%PHASE_RHO(imat, elemid)
164 phase_presii(imat) = multi_fvm%PHASE_PRES(imat, elemid)
165 ENDDO
166 ELSE
167 phase_alphaii(1) = one
168 phase_eintii(1) = eintii
169 phase_rhoii(1) = rhoii
170 phase_presii(1) = pii
171 ENDIF
172 sspii = multi_fvm%SOUND_SPEED(elemid)
173 normal_velii = vxii * nx + vyii * ny + vzii * nz
174
175C Boundary "GHOST" element
176 rhojj = zero
177 DO imat = 1, nbmat
178 phase_rhojj(imat) = fvm_inlet_data%VAL_RHO(imat)
179 ifunc = fvm_inlet_data%FUNC_RHO(imat)
180 IF (ifunc > 0) THEN
181 phase_rhojj(imat) = phase_rhojj(imat) * func_value(ifunc)
182 ELSEIF (ifunc == 0) THEN
183 phase_rhojj(imat) = phase_rhojj(imat)
184 ELSEIF (ifunc == -1) THEN
185 phase_rhojj(imat) = phase_rhoii(imat)
186 ENDIF
187 phase_alphajj(imat) = fvm_inlet_data%VAL_ALPHA(imat)
188 ifunc = fvm_inlet_data%FUNC_ALPHA(imat)
189 IF (ifunc > 0) THEN
190 phase_alphajj(imat) = phase_alphajj(imat) * func_value(ifunc)
191 ELSEIF (ifunc == 0) THEN
192 phase_alphajj(imat) = phase_alphajj(imat)
193 ELSEIF (ifunc == -1) THEN
194 phase_alphajj(imat) = phase_alphaii(imat)
195 ENDIF
196 rhojj = rhojj + phase_rhojj(imat) * phase_alphajj(imat)
197 ENDDO
198 IF (fvm_inlet_data%FORMULATION == 2) THEN
199C VE formulation
200 sspjj = zero
201 pjj = zero
202 eintjj = zero
203 DO imat = 1, nbmat
204 phase_eintjj(imat) = fvm_inlet_data%VAL_PRES(imat)
205 ifunc = fvm_inlet_data%FUNC_PRES(imat)
206 IF (ifunc > 0) THEN
207 phase_eintjj(imat) = phase_eintjj(imat) * func_value(ifunc)
208 ELSEIF (ifunc == 0) THEN
209 phase_eintjj(imat) = phase_eintjj(imat)
210 ELSEIF (ifunc == -1) THEN
211 phase_eintjj(imat) = phase_eintii(imat)
212 ENDIF
213 IF (phase_alphajj(imat) > zero) THEN
214 CALL multi_submatlaw(
215 1 0, matlaw(imat), local_matid(imat), 1,
216 2 phase_eintjj(imat), phase_presjj(imat), phase_rhojj(imat), phase_sspjj(imat),
217 3 dummy2, dummy, pm, ipm,
218 4 npropm, npropmi, dummy, dummy2,
219 5 dummy, multi_fvm%BFRAC(imat,elemid), multi_fvm%TBURN(elemid), dummy ,
220 6 dummy(1), dummy, snpc , stf ,
221 7 npf, tf, dummy(1) , 1 ,
222 8 matparam(local_matid(imat)),nvartmp_eos, vartmp_eos ,nummat ,
223 9 aburn)
224 sspjj = sspjj + phase_alphajj(imat) * phase_rhojj(imat) * max(em20, phase_sspjj(imat))
225 pjj = pjj + phase_presjj(imat) * phase_alphajj(imat)
226 eintjj = eintjj + phase_alphajj(imat) * phase_eintjj(imat)
227 ENDIF
228 ENDDO
229
230 ELSEIF (fvm_inlet_data%FORMULATION == 1) THEN
231C VP formulation
232 sspjj = zero
233 pjj = zero
234 eintjj = zero
235 DO imat = 1, nbmat
236 phase_presjj(imat) = fvm_inlet_data%VAL_PRES(imat)
237 ifunc = fvm_inlet_data%FUNC_PRES(imat)
238 IF (ifunc > 0) THEN
239 phase_presjj(imat) = phase_presjj(imat) * func_value(ifunc)
240 ELSEIF (ifunc == 0) THEN
241 phase_presjj(imat) = phase_presjj(imat)
242 ELSEIF (ifunc == -1) THEN
243 phase_presjj(imat) = phase_presii(imat)
244 ENDIF
245 pjj = pjj + phase_presjj(imat) * phase_alphajj(imat)
246C Need to solve p(RHO, EINT) = PHASE_PRESJJ(IMAT) for EINT
247C INITIALIZE INTERNAL ENERGY
248 phase_eintjj(imat) = pm(23, local_matid(imat))
249 CALL multi_solve_eint(matlaw(imat), local_matid(imat), pm, ipm, npropm, npropmi,
250 . phase_eintjj(imat), phase_rhojj(imat), phase_presjj(imat), phase_sspjj(imat),
251 . multi_fvm%BFRAC(imat, elemid), multi_fvm%TBURN(elemid), dummy(1), dummy(1),
252 . dummy, dummy2, snpc, stf, npf, tf, dummy, 1, matparam(local_matid(imat)),
253 ? nvartmp_eos, vartmp_eos, nummat, aburn)
254
255 sspjj = sspjj + phase_alphajj(imat) * phase_rhojj(imat) * max(em20, phase_sspjj(imat))
256 eintjj = eintjj + phase_alphajj(imat) * phase_eintjj(imat)
257 ENDDO
258 ENDIF
259 IF (sspjj / rhojj > zero) THEN
260 sspjj = sqrt(sspjj / rhojj)
261 ELSE
262 sspjj = multi_fvm%SOUND_SPEED(elemid)
263 ENDIF
264 IF (fvm_inlet_data%VECTOR_VELOCITY == 0) THEN
265C Normal velocity imposed
266 normal_veljj = -fvm_inlet_data%VAL_VEL(1)
267 ifunc = fvm_inlet_data%FUNC_VEL(1)
268 IF (ifunc > 0) THEN
269 normal_veljj = normal_veljj * func_value(ifunc)
270 ELSEIF (ifunc == 0) THEN
271 normal_veljj = normal_veljj
272 ELSEIF (ifunc == -1) THEN
273 normal_veljj = normal_velii
274 ENDIF
275C NORMAL_VELJJ = NORMAL_VELII
276
277 vxjj = normal_veljj * nx
278 vyjj = normal_veljj * ny
279 vzjj = normal_veljj * nz
280 ELSE
281 vxjj = fvm_inlet_data%VAL_VEL(1)
282 ifunc = fvm_inlet_data%FUNC_VEL(1)
283 IF (ifunc > 0) THEN
284 vxjj = vxjj * func_value(ifunc)
285 ELSEIF (ifunc == -1) THEN
286 vxjj = vxii
287 ENDIF
288 vyjj = fvm_inlet_data%VAL_VEL(2)
289 ifunc = fvm_inlet_data%FUNC_VEL(2)
290 IF (ifunc > 0) THEN
291 vyjj = vyjj * func_value(ifunc)
292 ELSEIF (ifunc == -1) THEN
293 vyjj = vyii
294 ENDIF
295 vzjj = fvm_inlet_data%VAL_VEL(3)
296 ifunc = fvm_inlet_data%FUNC_VEL(3)
297 IF (ifunc > 0) THEN
298 vzjj = vzjj * func_value(ifunc)
299 ELSEIF (ifunc == -1) THEN
300 vzjj = vzii
301 ENDIF
302 normal_veljj = vxjj * nx + vyjj * ny + vzjj * nz
303 ENDIF
304 veljj2 = vxjj**2 + vyjj**2 + vzjj**2
305C HLL wave speed estimates
306 sl = min(normal_velii - sspii, normal_veljj - sspjj)
307 sr = max(normal_velii + sspii, normal_veljj + sspjj)
308
309C Intermediate wave speed
310 sstar = pjj - pii + rhoii * normal_velii * (sl - normal_velii) - rhojj * normal_veljj * (sr - normal_veljj)
311 sstar = sstar / (rhoii * (sl - normal_velii) - rhojj * (sr - normal_veljj))
312
313 pstar = pii + rhoii * (sstar - normal_velii) * (sl - normal_velii)
314 pp(1) = zero
315 pp(2) = pstar * nx
316 pp(3) = pstar * ny
317 pp(4) = pstar * nz
318 pp(5) = sstar * pstar
319
320 IF (sl > normalw) THEN
321 vii(1) = rhoii
322 vii(2) = rhoii * vxii
323 vii(3) = rhoii * vyii
324 vii(4) = rhoii * vzii
325 vii(5) = eintii + half * rhoii * velii2
326! Normal physical flux current element
327 fii(1) = vii(1) * normal_velii
328 fii(2) = vii(2) * normal_velii + pii * nx
329 fii(3) = vii(3) * normal_velii + pii * ny
330 fii(4) = vii(4) * normal_velii + pii * nz
331 fii(5) = (vii(5) + pii) * normal_velii
332! /th/surf : pressure
333 psurf = pii
334C Take the fluxes of cell II
335C ===
336C Global fluxes
337 multi_fvm%FLUXES(1:5, kface, elemid) = (fii(1:5) - normalw * vii(1:5)) * surf
338 multi_fvm%FLUXES(6, kface, elemid) = normal_velii * surf
339C ===
340C Submaterial fluxes
341 IF (nbmat > 1) THEN
342 DO imat = 1, nbmat
343 alphaii = phase_alphaii(imat)
344 sub_rhoii = phase_rhoii(imat) ! ALPHA_RHO
345 sub_rhoeintii = phase_eintii(imat)
346 sub_viistar(1) = alphaii
347 sub_viistar(2) = alphaii * sub_rhoii ! ALPHA_RHO
348 sub_viistar(3) = alphaii * sub_rhoeintii
349 sub_fiistar(1:3) = sub_viistar(1:3) * normal_velii
350 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
351 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
352 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
353 ENDDO
354 ENDIF
355C
356 ELSEIF (sl <= normalw .AND. normalw <= sstar) THEN
357 vii(1) = rhoii
358 vii(2) = rhoii * vxii
359 vii(3) = rhoii * vyii
360 vii(4) = rhoii * vzii
361 vii(5) = eintii + half * rhoii * velii2
362! Normal physical flux current element
363 fii(1) = vii(1) * normal_velii
364 fii(2) = vii(2) * normal_velii + pii * nx
365 fii(3) = vii(3) * normal_velii + pii * ny
366 fii(4) = vii(4) * normal_velii + pii * nz
367 fii(5) = (vii(5) + pii) * normal_velii
368! /th/surf : pressure
369 psurf = pii
370C Take intermediate state flux (HLLC scheme)
371C ===
372C Global fluxes
373 viistar(1:5) = fii(1:5) - (sl) * vii(1:5) - pp(1:5)
374 viistar(1:5) = viistar(1:5) / (sstar - sl)
375 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
376 multi_fvm%FLUXES(1:5, kface, elemid) = (fiistar(1:5) - normalw * viistar(1:5)) * surf
377 multi_fvm%FLUXES(6, kface, elemid) = sstar * surf
378C ===
379C Submaterial fluxes
380 IF (nbmat > 1) THEN
381 DO imat = 1, nbmat
382 matlaw(imat) = ipm(2, local_matid(imat))
383 alphastar = phase_alphaii(imat)
384 sub_rhostar = phase_rhoii(imat) * (normal_velii - sl) / (sstar - sl)
385 IF (alphastar > zero) THEN
386 sub_rhoii = phase_rhoii(imat)
387 sub_rhoeintii = phase_eintii(imat)
388 sub_pii = phase_presii(imat)
389 sub_estar = phase_eintii(imat) / phase_rhoii(imat) -
390 . phase_presii(imat) * (one / sub_rhostar - one / phase_rhoii(imat))
391
392 IF (sub_estar < zero) THEN
393 sub_estar = zero
394 ENDIF
395 ELSE
396 sub_estar = zero
397 ENDIF
398 sub_viistar(1) = alphastar
399 sub_viistar(2) = alphastar * sub_rhostar
400 sub_viistar(3) = alphastar * sub_rhostar * sub_estar
401
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
409 ELSE IF (sstar < normalw .AND. normalw <= sr) THEN
410 vii(1) = rhojj
411 vii(2) = rhojj * vxjj
412 vii(3) = rhojj * vyjj
413 vii(4) = rhojj * vzjj
414 vii(5) = eintjj + half * rhojj * veljj2
415! Normal physical flux current element
416 fii(1) = vii(1) * normal_veljj
417 fii(2) = vii(2) * normal_veljj + pjj * nx
418 fii(3) = vii(3) * normal_veljj + pjj * ny
419 fii(4) = vii(4) * normal_veljj + pjj * nz
420 fii(5) = (vii(5) + pjj) * normal_veljj
421! /th/surf : pressure
422 psurf = pjj
423C Take intermediate state flux (HLLC scheme)
424C ===
425C Global fluxes
426 viistar(1:5) = fii(1:5) - (sr) * vii(1:5) - pp(1:5)
427 viistar(1:5) = viistar(1:5) / (sstar - sr)
428 fiistar(1:5) = viistar(1:5) * sstar + pp(1:5)
429 multi_fvm%FLUXES(1:5, kface, elemid) = (fiistar(1:5) - normalw * viistar(1:5)) * surf
430 multi_fvm%FLUXES(6, kface, elemid) = sstar * surf
431C ===
432C Submaterial fluxes
433 IF (nbmat > 1) THEN
434 DO imat = 1, nbmat
435 matlaw(imat) = ipm(2, local_matid(imat))
436 alphastar = phase_alphajj(imat)
437 sub_rhostar = phase_rhojj(imat) * (normal_veljj - sr) / (sstar - sr)
438 IF (alphastar > zero) THEN
439 sub_rhoii = phase_rhojj(imat)
440 sub_rhoeintii = phase_eintjj(imat)
441 sub_pii = phase_presjj(imat)
442 sub_estar = phase_eintjj(imat) / phase_rhojj(imat) -
443 . phase_presjj(imat) * (one / sub_rhostar - one / phase_rhojj(imat))
444 IF (sub_estar < zero) THEN
445 sub_estar = zero
446 ENDIF
447 ELSE
448 sub_estar = zero
449 ENDIF
450
451 sub_viistar(1) = phase_alphajj(imat)
452 sub_viistar(2) = phase_alphajj(imat) * sub_rhostar
453 sub_viistar(3) = phase_alphajj(imat) * sub_rhostar * sub_estar
454 sub_fiistar(1:3) = sub_viistar(1:3) * sstar
455 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
456 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
457 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
458 ENDDO
459 ENDIF
460 ELSE IF (sr < normalw) THEN
461 vii(1) = rhojj
462 vii(2) = rhojj * vxjj
463 vii(3) = rhojj * vyjj
464 vii(4) = rhojj * vzjj
465 vii(5) = eintjj + half * rhojj * veljj2
466! Normal physical flux current element
467 fii(1) = vii(1) * normal_veljj
468 fii(2) = vii(2) * normal_veljj + pjj * nx
469 fii(3) = vii(3) * normal_veljj + pjj * ny
470 fii(4) = vii(4) * normal_veljj + pjj * nz
471 fii(5) = (vii(5) + pjj) * normal_veljj
472! /th/surf : pressure
473 psurf = pjj
474C Take the fluxes of cell II
475C ===
476C Global fluxes
477 multi_fvm%FLUXES(1:5, kface, elemid) = (fii(1:5) - normalw * vii(1:5)) * surf
478 multi_fvm%FLUXES(6, kface, elemid) = normal_veljj * surf
479C ===
480C Submaterial fluxes
481 IF (nbmat > 1) THEN
482 DO imat = 1, nbmat
483 alphaii = phase_alphajj(imat)
484 sub_rhoii = phase_rhojj(imat) ! ALPHA_RHO
485 sub_rhoeintii = phase_eintjj(imat)
486 sub_viistar(1) = alphaii
487 sub_viistar(2) = alphaii * sub_rhoii ! ALPHA_RHO
488 sub_viistar(3) = alphaii * sub_rhoeintii
489 sub_fiistar(1:3) = sub_viistar(1:3) * normal_veljj
490 multi_fvm%SUBVOL_FLUXES(imat, kface, elemid) = (sub_fiistar(1) - normalw * sub_viistar(1)) * surf
491 multi_fvm%SUBMASS_FLUXES(imat, kface, elemid) = (sub_fiistar(2) - normalw * sub_viistar(2)) * surf
492 multi_fvm%SUBENER_FLUXES(imat, kface, elemid) = (sub_fiistar(3) - normalw * sub_viistar(3)) * surf
493 ENDDO
494 ENDIF
495 ELSE
496 CALL arret(2)
497C
498 ENDIF
499
500 ! /TH/SURF (massflow,velocity,pressure,acceleration)
501 fsavsurf(2,id_surf) = fsavsurf(2,id_surf) + multi_fvm%FLUXES(1, kface, elemid) ! += rho.S.un
502 fsavsurf(3,id_surf) = fsavsurf(3,id_surf) + multi_fvm%FLUXES(1, kface, elemid) / vii(1) ! += S.un
503 fsavsurf(4,id_surf) = fsavsurf(4,id_surf) + surf*psurf ! += S.P
504 fsavsurf(6,id_surf) = fsavsurf(6,id_surf) + multi_fvm%FLUXES(1, kface, elemid) * timestep ! m <- m+dm (cumulative value)
505 ENDDO
506C-----------------------------------------------
507C E n d o f s u b r o u t i n e
508C-----------------------------------------------
509 END SUBROUTINE multi_inlet_ebcs
510 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)
Definition multi_inlet_ebcs.F:48
multi_submatlaw_mod
Definition multi_submatlaw.F:34
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:63
th_surf_mod
OPTION /TH/SURF outputs of Pressure and Area needed Tabs.
Definition th_surf_mod.F:60
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:99
arret
subroutine arret(nn)
Definition arret.F:87