OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fv_up_switch.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!|| fv_up_switch ../engine/source/airbag/fv_up_switch.F
25!||--- called by ------------------------------------------------------
26!|| fvbag0 ../engine/source/airbag/fvbag0.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../engine/source/output/message/message.F
29!|| arret ../engine/source/system/arret.F
30!|| fvinjt6 ../engine/source/airbag/fvinjt6.F
31!|| fvinjt8 ../engine/source/airbag/fvinjt8.F
32!|| fvinjt8_1 ../engine/source/airbag/fvinjt8_1.F
33!|| fvtemp ../engine/source/airbag/fvtemp.F
34!|| fvvent0 ../engine/source/airbag/fvvent0.f
35!|| get_u_func ../engine/source/user_interface/ufunc.F
36!|| spmd_fvb_gath ../engine/source/mpi/airbags/spmd_fvb_gath.F
37!|| spmd_fvb_igath ../engine/source/mpi/airbags/spmd_fvb_igath.F
38!||--- uses -----------------------------------------------------
39!|| anim_mod ../common_source/modules/output/anim_mod.F
40!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
41!|| finter_mixed_mod ../engine/source/tools/finter_mixed.F90
42!|| fvbag_mod ../engine/share/modules/fvbag_mod.F
43!|| fvmbag_meshcontrol_mod ../common_source/modules/airbag/fvmbag_meshcontrol_mod.F
44!|| message_mod ../engine/share/message_module/message_mod.F
45!|| python_funct_mod ../common_source/modules/python_mod.F90
46!|| sensor_mod ../common_source/modules/sensor_mod.F90
47!||====================================================================
48 SUBROUTINE fv_up_switch(NN, NEL, ELEM, NJET, NPOLY,LENH,NBA,
49 2 IBAGJET, RBAGJET, NVENT , IBAGHOL, RBAGHOL ,
50 3 P , RHO , TK , U , SSPK ,
51 4 X , V , A , NSENSOR, SENSOR_TAB,
52 5 FSAV , NPC , TF , IVOLU , RVOLU ,
53 6 MPOLH , QPOLH , EPOLH ,
54 7 PPOLH , RPOLH , GPOLH ,
55 8 NPOLH , IFVNOD , RFVNOD , IFVTRI ,
56 9 IFVPOLY, IFVTADR, IFVPOLH ,
57 A IFVPADR, INFO , NNS , NNTR , IFV ,
58 B NPOLHA , DLH , CPAPOLH ,
59 C CPBPOLH, CPCPOLH, RMWPOLH ,
60 D ITAGEL , ELSINI , ICONTACT , IDPOLH ,
61 E ELFMASS, ELFVEL , IBUFA , ELEMA , TAGELA ,
62 F PA , RHOA , TKA , UA , BRNA ,
63 G NNA , NTGA , IBPOLH , DTPOLH , NNT ,
64 H NELT , XXXA , VVVA , NCONA , POROSITY,
65 I ITYP , IGEO , SSPKA ,
66 J GEO , PM , IPM , TPOLH , ELFEHPY ,
67 K CPDPOLH, CPEPOLH, CPFPOLH ,
68 L ELTG , IPARG , MATTG ,
69 M IGROUPTG,IGROUPC, ELBUF_TAB, CFL_COEF,
70 N PDISP_OLD, PDISP, WFEXT, PYTHON )
71C-----------------------------------------------
72C M O D U L E S
73C-----------------------------------------------
74 USE fvbag_mod
75 USE message_mod
76 USE elbufdef_mod
77 USE fvmbag_meshcontrol_mod
78! USE OMP_LIB
79 USE sensor_mod
80 USE anim_mod
81 USE python_funct_mod
82 use finter_mixed_mod
83C-----------------------------------------------
84C I M P L I C I T T Y P E S
85C-----------------------------------------------
86#include "implicit_f.inc"
87C-----------------------------------------------
88C C O M M O N B L O C K S
89C-----------------------------------------------
90#include "com01_c.inc"
91#include "com04_c.inc"
92#include "com06_c.inc"
93#include "com08_c.inc"
94#include "scr02_c.inc"
95#include "scr07_c.inc"
96#include "param_c.inc"
97#include "units_c.inc"
98#include "task_c.inc"
99#include "warn_c.inc"
100#include "tabsiz_c.inc"
101#include "mvsiz_p.inc"
102C-----------------------------------------------
103C D U M M Y A R G U M E N T S
104C-----------------------------------------------
105 INTEGER, INTENT(IN) :: NPOLY, LENH, NPOLH, NBA,NSENSOR
106 INTEGER, INTENT(IN) :: NN, NEL, NELT, ELEM(3, NELT), NJET, IBAGJET(NIBJET, NJET),
107 . NVENT, NPC(SNPC), IFVNOD(3,NNS), IFVTRI(6,NNTR),
108 . IFVPOLY(NNTR),IFVTADR(NPOLY+1), IFVPOLH(LENH), IFVPADR(NPOLH + 1),
109 . NNS, NNTR, IFV, NPOLHA, ITAGEL(NELT), ICONTACT(SICONTACT),
110 . IDPOLH(NPOLH), IBUFA(NNA), ELEMA(3,NTGA), TAGELA(NTGA), BRNA(8,NBA),
111 . NNA, NTGA, IBPOLH(NPOLH), NNT, NCONA(16,NNA),
112 . ITYP, IGEO(NPROPGI,*), IPM(NPROPMI,NUMMAT), IPARG(NPARG,NGROUP)
113 INTEGER, INTENT(INOUT) :: IVOLU(NIMV), INFO, IBAGHOL(NIBHOL, NVENT)
114 INTEGER, INTENT(IN) :: ELTG(NELT), MATTG(NELT), IGROUPTG(NUMELTG), IGROUPC(NUMELC)
115
116 MY_REAL, INTENT(IN) ::
117 . X(3,NUMNOD), V(3,NUMNOD), A(3,NUMNOD),
118 . TF(STF), XXXA(3,NNA), VVVA(3,NNA),
119 . GEO(NPROPG,NUMGEO), PM(NPROPM,NUMMAT), CFL_COEF, PDISP_OLD
120 MY_REAL, INTENT(INOUT) :: POROSITY(NELT - NEL)
121 MY_REAL, INTENT(INOUT) :: RBAGJET(NRBJET,NJET)
122 MY_REAL, INTENT(INOUT) ::
123 . P(NNT), RHO(NNT), TK(NNT), U(3,NNT), SSPK(NNT), MPOLH(NPOLH), QPOLH(3,NPOLH),
124 . EPOLH(NPOLH), PPOLH(NPOLH), RPOLH(NPOLH), GPOLH(NPOLH), RFVNOD(2,NNS), DLH,
125 . CPAPOLH(NPOLH), CPBPOLH(NPOLH), CPCPOLH(NPOLH), RMWPOLH(NPOLH), ELSINI(NELT),
126 . ELFMASS(NELT), ELFVEL(NELT), PA(NNA), RHOA(NNA), TKA(NNA),SSPKA(NNA), UA(3,NNA),
127 . dtpolh(npolh), fsav(nthvki), tpolh(npolh), cpdpolh(npolh), cpepolh(npolh),
128 . cpfpolh(npolh), rbaghol(nrbhol,nvent), pdisp, elfehpy(nelt), rvolu(nrvolu)
129 TYPE(elbuf_struct_), DIMENSION(NGROUP), INTENT(IN) :: ELBUF_TAB
130 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
131 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
132 TYPE(PYTHON_),INTENT(IN) :: PYTHON
133C-----------------------------------------------
134C L O C A L V A R I A B L E S
135C-----------------------------------------------
136 INTEGER I, II, IINJ, IEL, N1, N2, N3, JEL,
137 . NN1, NN2, NN3, J, JJ, K, KK, ISENS, IVEL,
138 . itagp(npolh), nv, ilvout, ipri,
139 . ivent, idef, iporp,
140 . itvent, phdt, i1, i2, icont(nnt),ivdp,ittf, idpdef,
141 . idtpdef,itypl,
142 . prepare_anim,iventyp
143 INTEGER TITREVENT(20)
144
145 my_real
146 . nodarea(nnt), x1, y1, z1, x2, y2, z2, x3, y3, z3,
147 . x12, y12, z12, x13, y13, z13, nrx, nry, nrz, area2,
148 . elarea(nelt), norm(3,nelt), ksi, eta,
149 . pnorm(3,nntr), pvolu(npolh), area, nx, ny,
150 . nz, gamai, cpai, cpbi, cpci, rmwi, ti,
151 . cpdi, cpei, cpfi,
152 . rhoi, datainj(6,njet), scalt,
153 . fvel, tstart, tsg, injvel,
154 . cpa, cpb, cpc, efac,
155 . pu(3,npolh), gama, fel(nelt), dm(npolh),
156 . dq(3,npolh), de(npolh), dmi, rho1, re1, ux1,
157 . uy1, uz1, rho2, re2, ux2, uy2, uz2, vfx, vfy, vfz, vrel(nelt), ss(nelt), ss_,
158 . alpha, rhom, rem, ruxm, ruym, ruzm, p1, p2,
159 . uel(3,nelt), area1, area3, ff, temp, pext, volg, volph,
160 . areap, areael, qa, qb, dtx, al, dd, ad, ssp,
161 . qx, vv, dmini(npolh), dmcpa(npolh), dmcpb(npolh),
162 . dmcpc(npolh), dmrmw(npolh), msini,
163 . rmw, cp, cv, rhoel(nelt), tkel(nelt), sspel(nelt), gama1, gama2,
164 . qvisc(npolh), amtot, pmean, pmean2, elsout(nelt),tmp(nel),
165 . area_vent(nvent), pm_vent(nvent), scalp, scals,
166 . aout, deri, pdef, pvoltmp,
167 . dtpdefi, dtpdefc, tvent, poro(nelt), vola(nna),
168 . x23, y23, z23, x31, y31, z31, l12, l23, l31, h1,
169 . fac1, fac2, fac3, area_old, pp, dls(npolh), vmax, uu,
170 . pcrit, vx1, vy1, vz1, vx2, vy2, vz2,
171 . vx3, vy3, vz3, vvx, vvy, vvz, vel(nelt), fac,
172 . aoutot, svent(nvent), sventot, xxx(3,nnt),
173 . vvv(3,nnt), fvdp, rbid, ttf,
174 . dmcpd(npolh), dmcpe(npolh), dmcpf(npolh),
175 . cpd, cpe, cpf, temp0, pel(nelt), pstag, volno,tstope,
176 . enint, massflow, wfext0, eldmass(nel), eldehpy(nel), dmout, dhout,
177 . xx1, xx2, xx3, yy1, yy2, yy3, zz1, zz2, zz3, nnx, nny, nnz,
178 . pres_av, cp_av, cv_av, rgas_av, rho_av
179
180
181 LOGICAL :: EXIT_NEG_VOL,INJECTION_STARTED
182 LOGICAL :: UP_SWITCH
183C
184 MY_REAL
185 . GET_U_FUNC
186 EXTERNAL GET_U_FUNC
187 CHARACTER*20 VENTTITLE
188 INTEGER :: NODE_ID(5), INODE
189 MY_REAL :: TAB_PVOL(NELT), DOT_PROD
190 MY_REAL :: MOMENTUM_FLUX_X(NELT), MOMENTUM_FLUX_Y(NELT), MOMENTUM_FLUX_Z(NELT),
191 . MASS_FLUX(NELT), ENERGY_FLUX(NELT),
192 . CPA_FLUX(NELT), CPB_FLUX(NELT), CPC_FLUX(NELT), CPD_FLUX(NELT),
193 . CPE_FLUX(NELT), CPF_FLUX(NELT), RGAS_FLUX(NELT)
194 INTEGER(8) :: VEC_PTR_PLUS, VEC_PTR_MINUS
195 INTEGER :: COUNT, IDT_FVMBAG, PARAM_L_TYPE
196 MY_REAL :: CPAM, CPBM, CPCM, RGASM, CPDM, CPEM, CPFM
197 MY_REAL :: DTOLD, LAMBDA
198
199C-----------------------------------------------
200C C O M M E N T S
201C-----------------------------------------------
202C NEL : number of triangles (external surface only)
203C NELT : number of triangles (both internal & external surfaces)
204C NNS : number of nodes of TETRA mesh
205C NNA : number of 'additional' nodes
206C NNT : total number of nodes
207C NNTR : number of triangles (all polyhedra)
208C NPOLH : number of polyhedron
209C NPOLY : number of polygons
210
211C-----------------------------------------------
212C S O U R C E L I N E S
213C-----------------------------------------------
214 volph = zero
215 temp = zero
216
217 IF(nbgauge > 0 ) THEN
218 ALLOCATE(fvspmd(ifv)%GGG(3,nnt))
219 ALLOCATE(fvspmd(ifv)%GGA(3,nna))
220 ENDIF
221 ALLOCATE(fvspmd(ifv)%AAA(3,nnt))
222 fvspmd(ifv)%AAA(1:3,1:nnt) = zero
223
224! =================== !
225! PARITH/ON TREATMENT !
226! =================== !
227 IF (nspmd == 1) THEN
228 DO i=1,fvspmd(ifv)%NN_L+fvspmd(ifv)%NNI_L
229 i1=fvspmd(ifv)%IBUF_L(1,i)
230 i2=fvspmd(ifv)%IBUF_L(2,i)
231 xxx(1,i1)=x(1,i2)
232 xxx(2,i1)=x(2,i2)
233 xxx(3,i1)=x(3,i2)
234 ENDDO
235
236 DO i=1,fvspmd(ifv)%NN_L+fvspmd(ifv)%NNI_L
237 i1=fvspmd(ifv)%IBUF_L(1,i)
238 i2=fvspmd(ifv)%IBUF_L(2,i)
239 vvv(1,i1)=v(1,i2)
240 vvv(2,i1)=v(2,i2)
241 vvv(3,i1)=v(3,i2)
242 ENDDO
243
244 IF (intbag/=0.AND.nvent>0.AND.ivolu(39)/=0) THEN
245 DO i=1,fvspmd(ifv)%NN_L+fvspmd(ifv)%NNI_L
246 i1=fvspmd(ifv)%IBUF_L(1,i)
247 i2=fvspmd(ifv)%IBUF_L(2,i)
248 icont(i1)=icontact(i2)
249 ENDDO
250 ENDIF
251
252 ELSE
253 CALL spmd_fvb_gath(ifv, x, xxx, rbid, rbid, 4)
254
255 CALL spmd_fvb_gath(ifv, v, vvv, rbid, rbid, 4)
256
257 IF (intbag/=0.AND.nvent>0.AND.ivolu(39)/=0)
258 . CALL spmd_fvb_igath(ifv, icontact, icont)
259
260
261 IF (ispmd/=fvspmd(ifv)%PMAIN-1) RETURN
262 END IF
263
264! ================================= !
265! REBUILD CONNECTIVITY AFTER SWITCH
266! ================================= !
267
268! INITIALIZATION PHASE, DONE ONLY ONCE
269! ------------------------------------
270
271! NB, FOR AN ELEMENT IEL, THE 'MINUS SIDE' IS BEHIND THE NORMAL, WHILE 'PLUS SIDE' IS IN THE DIRECTION OF
272! THE NORMAL
273
274 IF (ivolu(74) == -2) THEN
275 ALLOCATE(fvdata(ifv)%TRI_TO_ELEM(2, nelt))
276 fvdata(ifv)%TRI_TO_ELEM(:, :) = 0
277 ALLOCATE(fvdata(ifv)%ELEM_TO_TRI(npolh))
278 CALL intvector_create(vec_ptr_plus)
279 CALL intvector_create(vec_ptr_minus)
280 DO i = 1, npolh
281! CLEAR INT VECTORS
282 CALL intvector_clear(vec_ptr_plus)
283 CALL intvector_clear(vec_ptr_minus)
284 count = 0
285! LOOP ON POLYHEDRAL FACES
286 DO j = ifvpadr(i), ifvpadr(i+1)-1
287 jj = ifvpolh(j)
288! LOOP ON TRIANGLES OF THE POLYHEDRAL FACE
289 DO k = ifvtadr(jj), ifvtadr(jj+1)-1
290 kk = ifvpoly(k)
291 iel = ifvtri(4, kk)
292 n1 = ifvtri(1, kk)
293 IF (ifvnod(1, n1) /= 2) THEN
294 CALL abort()
295 ENDIF
296 n1 = ifvnod(3, n1)
297 x1 = xxxa(1, n1)
298 y1 = xxxa(2, n1)
299 z1 = xxxa(3, n1)
300 n2 = ifvtri(2, kk)
301 IF (ifvnod(1, n2) /= 2) THEN
302 CALL abort()
303 ENDIF
304 n2 = ifvnod(3, n2)
305 x2 = xxxa(1, n2)
306 y2 = xxxa(2, n2)
307 z2 = xxxa(3, n2)
308 n3 = ifvtri(3, kk)
309 IF (ifvnod(1, n3) /= 2) THEN
310 CALL abort()
311 ENDIF
312 n3 = ifvnod(3, n3)
313 x3 = xxxa(1, n3)
314 y3 = xxxa(2, n3)
315 z3 = xxxa(3, n3)
316 nx = (y2 - y1) * (z3 - z1) - (z2 - z1) * (y3 - y1)
317 ny = (z2 - z1) * (x3 - x1) - (x2 - x1) * (z3 - z1)
318 nz = (x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)
319 IF (iel /= 0) THEN
320 count = count + 1
321 nn1 = elem(1, abs(iel))
322 nn2 = elem(2, abs(iel))
323 nn3 = elem(3, abs(iel))
324 xx1 = xxx(1, nn1)
325 yy1 = xxx(2, nn1)
326 zz1 = xxx(3, nn1)
327 xx2 = xxx(1, nn2)
328 yy2 = xxx(2, nn2)
329 zz2 = xxx(3, nn2)
330 xx3 = xxx(1, nn3)
331 yy3 = xxx(2, nn3)
332 zz3 = xxx(3, nn3)
333 nnx = (yy2 - yy1) * (zz3 - zz1) - (zz2 - zz1) * (yy3 - yy1)
334 nny = (zz2 - zz1) * (xx3 - xx1) - (xx2 - xx1) * (zz3 - zz1)
335 nnz = (xx2 - xx1) * (yy3 - yy1) - (yy2 - yy1) * (xx3 - xx1)
336 dot_prod = nx * nnx + ny * nny + nz * nnz
337! EXTERNAL (IEL > 0) OR INTERNAL (IEL < 0) SURFACE ELEMENT
338 IF (iel < 0) THEN
339! INTERNAL, THERE IS A POTENTIAL NEIGHBOR
340 iel = -iel
341 i1 = ifvtri(5, kk)
342 i2 = ifvtri(6, kk)
343 IF (i1 == i2) THEN
344 CALL abort()
345 ENDIF
346 fvdata(ifv)%TRI_TO_ELEM(1, iel) = i1
347 fvdata(ifv)%TRI_TO_ELEM(2, iel) = i2
348 IF (i1 == i) THEN
349 IF (dot_prod >= zero) THEN
350 CALL intvector_push_back(vec_ptr_minus, iel)
351 ELSE
352 CALL intvector_push_back(vec_ptr_plus, iel)
353 fvdata(ifv)%TRI_TO_ELEM(1, iel) = i2
354 fvdata(ifv)%TRI_TO_ELEM(2, iel) = i1
355 ENDIF
356 ELSE
357 IF (dot_prod >= zero) THEN
358 CALL intvector_push_back(vec_ptr_plus, iel)
359 ELSE
360 CALL intvector_push_back(vec_ptr_minus, iel)
361 fvdata(ifv)%TRI_TO_ELEM(1, iel) = i2
362 fvdata(ifv)%TRI_TO_ELEM(2, iel) = i1
363 ENDIF
364 ENDIF
365 ELSE
366! EXTERNAL SURFACE : NO NEIGHBOR
367 IF (ifvtri(5, kk) /= ifvtri(6, kk)) THEN
368 CALL abort()
369 ENDIF
370 fvdata(ifv)%TRI_TO_ELEM(1, iel) = i
371 fvdata(ifv)%TRI_TO_ELEM(2, iel) = i
372 IF (dot_prod >= zero) THEN
373 CALL intvector_push_back(vec_ptr_minus, iel)
374 ELSE
375 CALL intvector_push_back(vec_ptr_plus, iel)
376 ENDIF
377 ENDIF
378 ELSE
379 CALL abort()
380 ENDIF
381 ENDDO
382 ENDDO
383 CALL intvector_get_size(vec_ptr_minus, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE)
384 CALL intvector_get_size(vec_ptr_plus, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE)
385 ALLOCATE(fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE))
386 CALL intvector_copy_to(vec_ptr_minus, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(1))
387 ALLOCATE(fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE))
388 CALL intvector_copy_to(vec_ptr_plus, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(1))
389 ENDDO
390 ivolu(74) = -20
391 CALL intvector_delete(vec_ptr_plus)
392 CALL intvector_delete(vec_ptr_minus)
393 ENDIF
394
395
396!$OMP PARALLEL PRIVATE(I,J)
397!$OMP+ PRIVATE(I1,I2,KSI,ETA,N1,N2,N3,IEL,X1,Y1,Z1,FAC,NN1,NN2,NN3)
398!$OMP+ PRIVATE(X2,Y2,Z2,X3,Y3,Z3,VX1,VX2,VX3,VY1,VY2,VY3,VZ1,VZ2,VZ3)
399!$OMP+ PRIVATE(X12,Y12,Z12,X13,Y13,Z13,AREA2,NRX,NRY,NRZ,JJ,AREA,KK,K)
400!$OMP+ PRIVATE(NX,NY,NZ,PVOLTMP)
401
402! ================================================= !
403! SURFACE AND NORMAL COMPUTATION FOR FABRIC ELEMENT !
404! ================================================= !
405
406!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
407 DO i=1,nnt
408 nodarea(i)=zero
409 ENDDO
410!$OMP END DO
411!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
412 DO iel=1,nelt
413 n1=elem(1,iel)
414 n2=elem(2,iel)
415 n3=elem(3,iel)
416 x1=xxx(1,n1)
417 x2=xxx(1,n2)
418 x3=xxx(1,n3)
419 y1=xxx(2,n1)
420 y2=xxx(2,n2)
421 y3=xxx(2,n3)
422 z1=xxx(3,n1)
423 z2=xxx(3,n2)
424 z3=xxx(3,n3)
425 x12=x2-x1
426 y12=y2-y1
427 z12=z2-z1
428 x13=x3-x1
429 y13=y3-y1
430 z13=z3-z1
431 nrx=y12*z13-z12*y13
432 nry=z12*x13-x12*z13
433 nrz=x12*y13-y12*x13
434 area2=sqrt(nrx**2+nry**2+nrz**2)
435 elarea(iel)=area2
436 norm(1,iel)=nrx/area2
437 norm(2,iel)=nry/area2
438 norm(3,iel)=nrz/area2
439 IF(iel<=nel) THEN
440 tmp(iel) = one_over_6 * (x1*nrx+y1*nry+z1*nrz)
441 ENDIF
442 pel(iel)=zero
443 ENDDO
444!$OMP END DO
445
446!$OMP SINGLE
447 volg=zero
448 DO iel=1,nelt
449 n1=elem(1,iel)
450 n2=elem(2,iel)
451 n3=elem(3,iel)
452 area2 = elarea(iel)
453 elarea(iel) = half * area2
454 nodarea(n1)=nodarea(n1)+one_over_6*area2
455 nodarea(n2)=nodarea(n2)+one_over_6*area2
456 nodarea(n3)=nodarea(n3)+one_over_6*area2
457 IF(iel<=nel) THEN
458 volg=volg+tmp(iel)
459 ENDIF
460 ENDDO
461C
462 IF (dt1==zero.OR.ivolu(39)==0) THEN
463 DO iel=1,nelt
464 elfmass(iel)=zero
465 elfehpy(iel)=zero
466 ENDDO
467 fsav(1:nthvki)=zero
468 ENDIF
469!$OMP END SINGLE
470
471! ============================ !
472! VOLUME OF THE FINITE VOLUMES !
473! ============================ !
474
475!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
476 DO iel = 1, nelt
477 nx = norm(1, iel)
478 ny = norm(2, iel)
479 nz = norm(3, iel)
480 n1 = elem(1, iel)
481 x1 = xxx(1, n1)
482 x2 = xxx(2, n1)
483 x3 = xxx(3, n1)
484 area = elarea(iel)
485 pvoltmp = third * area * (nx * x1 + ny * x2 + nz * x3)
486 tab_pvol(iel) = pvoltmp
487 ENDDO
488!$OMP END DO
489
490!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
491 DO i = 1, npolh
492 pvolu(i) = zero
493 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
494 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
495! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
496 pvolu(i) = pvolu(i) + tab_pvol(iel)
497 ENDDO
498 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
499 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
500! FINITE VOLUME I IS ON THE RIGHT SIDE OF ELEMENT IEL, HENCE A - SIGN
501 pvolu(i) = pvolu(i) - tab_pvol(iel)
502 ENDDO
503 ENDDO
504!$OMP END DO
505!$OMP END PARALLEL
506
507
508C---------------------------------------------------
509C VERIFICATIONS
510C---------------------------------------------------
511 areael=zero
512 DO iel=1,nel
513 areael=areael+elarea(iel)
514 ENDDO
515 areap=areael
516 rvolu(18)=areael
517 ilvout=ivolu(44)
518 volph=zero
519 exit_neg_vol = .false.
520 DO i=1,npolh
521 volph=volph+pvolu(i)
522 IF(pvolu(i) <= 0) exit_neg_vol = .true.
523 ENDDO
524
525 ipri=mod(ncycle,iabs(ncpri))
526 IF(exit_neg_vol) THEN
527 DO i=1,npolh
528 IF (pvolu(i)<=zero) THEN
529 info=1
530 ierr=ierr+1
531 CALL ancmsg(msgid=185,anmode=aninfo,
532 . i1=idpolh(i),r1=pvolu(i),i3=i,i4=npolh)
533 CALL arret(2)
534 ENDIF
535 ENDDO
536 ENDIF
537
538
539C---------------------------------------------------
540C COMPUTATION OF QUANTITIES IN POLYEDRA
541C---------------------------------------------------
542 IF (dt1==zero) THEN
543 gamai=rvolu(1)
544 gpolh(1:npolh)=gamai
545 cpai =rvolu(7)
546 cpapolh(1:npolh)=cpai
547 cpbi =rvolu(8)
548 cpbpolh(1:npolh)=cpbi
549 cpci =rvolu(9)
550 cpcpolh(1:npolh)=cpci
551 rmwi =rvolu(10)
552 rmwpolh(1:npolh)=rmwi
553 ti =rvolu(13)
554 tpolh(1:npolh)=ti
555 cpdi =rvolu(56)
556 cpdpolh(1:npolh)=cpdi
557 cpei =rvolu(57)
558 cpepolh(1:npolh)=cpei
559 cpfi =rvolu(58)
560 cpfpolh(1:npolh)=cpfi
561 rhoi =rvolu(62)
562 rpolh(1:npolh)=rhoi
563 ENDIF
564 IF (dt1==zero.OR.ivolu(39)==0) THEN
565 rhoi =rvolu(62)
566 efac =rvolu(66)
567 mpolh(1:npolh)=rhoi*pvolu(1:npolh)
568 epolh(1:npolh)=mpolh(1:npolh)*efac
569 pu(1,1:npolh) = rvolu(67)
570 pu(2,1:npolh) = rvolu(68)
571 pu(3,1:npolh) = rvolu(69)
572 qpolh(1,1:npolh) = mpolh(1:npolh)*pu(1,1:npolh)
573 qpolh(2,1:npolh) = mpolh(1:npolh)*pu(2,1:npolh)
574 qpolh(3,1:npolh) = mpolh(1:npolh)*pu(3,1:npolh)
575 ENDIF
576
577 pext =rvolu(3)
578 IF(ivolu(39) /= 0) THEN
579
580C---------------------------------------------------
581C INJECTORS
582C---------------------------------------------------
583 datainj(1:6,1:njet)=zero
584 DO iel=1,nelt
585 IF (itagel(iel)>0) THEN
586 iinj=itagel(iel)
587 datainj(1,iinj)=datainj(1,iinj)+elarea(iel)
588 ENDIF
589 ENDDO
590C
591 scalt =rvolu(26)
592 IF(ityp==6) THEN
593 CALL fvinjt6(njet, ibagjet, rbagjet, npc, tf,nsensor,
594 . sensor_tab,scalt, datainj, python )
595 ELSEIF(ityp==8) THEN
596 CALL fvinjt8(njet, ibagjet, rbagjet, npc, tf,nsensor,
597 . sensor_tab,scalt, igeo, geo, pm, ivolu, datainj,python)
598 ENDIF
599C
600 DO iinj=1,njet
601 isens=ibagjet(4,iinj)
602 fvel=rbagjet(15,iinj)
603 ivel=ibagjet(11,iinj)
604 IF(isens==0)THEN
605 tstart=zero
606 ELSE
607 tstart=sensor_tab(isens)%TSTART
608 ENDIF
609C ITTF : 1 SHIFT EVENT TIME (CAN BE REMOVED IN A LATER VERSION)
610C ITTF : 2 SHIFT EVENT TIME + EVENT FUNCTION
611C ITTF : 1 OR 2 +10 1 FLAG HAS BEEN ACTIVATED AND RVOLU(60) IS FILLED
612 ittf=ivolu(17)
613 IF (tt>=tstart.AND.(ittf==1.OR.ittf==2.OR.ittf==3))THEN
614 ittf=ittf+10
615 rvolu(60)=tstart
616 ivolu(17)=ittf
617 END IF
618 IF (tt>=tstart.AND.dt1>zero)THEN
619 tsg=(tt-tstart)*scalt
620 IF (ivel>0) THEN
621 injvel=fvel*finter_mixed(python,nfunct,ivel,tsg,npc,tf)
622 ELSE
623 injvel = fvel
624 ENDIF
625 ELSE
626 injvel=zero
627 ENDIF
628 datainj(3,iinj)=injvel
629 ENDDO
630C---------------------------------------------------
631C VENT HOLES
632C---------------------------------------------------
633 IF (intbag==0) THEN
634 poro(1:nelt)=zero
635 ELSE
636!$OMP PARALLEL PRIVATE(IEL,N1,N2,N3)
637!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
638 DO iel=1,nelt
639 n1=elem(1,iel)
640 n2=elem(2,iel)
641 n3=elem(3,iel)
642 poro(iel)=zero
643 IF (icont(n1)/=0) poro(iel)=poro(iel)+third
644 IF (icont(n2)/=0) poro(iel)=poro(iel)+third
645 IF (icont(n3)/=0) poro(iel)=poro(iel)+third
646 ENDDO
647!$OMP END DO
648!$OMP END PARALLEL
649 ENDIF
650C
651 pext =rvolu(3)
652 scalt=rvolu(26)
653 scalp=rvolu(27)
654 scals=rvolu(28)
655 ttf =rvolu(60)
656C-----------------------------------------------------------
657C EFFECTIVE LEAKAGE SURFACE OF TRIANGLES
658C CORRECTION OF SCALE FACTORS ON THE VENT HOLES SURFACE
659C-----------------------------------------------------------
660 CALL fvvent0(
661 1 elsout ,aoutot ,nvent ,nelt ,ittf ,
662 2 elarea ,elsini ,elem ,itagel ,svent ,
663 3 ibaghol ,rvolu ,rbaghol ,poro ,p ,
664 4 eltg ,iparg ,mattg ,nel ,porosity,
665 5 ipm ,pm ,elbuf_tab,igroupc ,igrouptg)
666C---------------------------------------------------
667C FINITE VOLUME BALANCES
668C---------------------------------------------------
669!$OMP PARALLEL PRIVATE(I,J,JJ,K,KK,IEL,JEL,AREA,N1,N2,N3)
670!$OMP+ PRIVATE(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X12,Y12,Z12,X23,Y23,Z23,Z31,Y31,X31)
671!$OMP+ PRIVATE(L12,L23,L31,H1,NX,NY,NZ,IINJ,DMI,IVENT,IDEF,ITVENT,AOUT)
672!$OMP+ PRIVATE(P1,RHO1,RE1,GAMA1,UX1,UY1,UZ1,UU,RHO2,VV,ETA,KSI,PSTAG)
673!$OMP+ PRIVATE(PCRIT,P2,VMAX,IVDP,FVDP,GAMAI,RHOI,DMOUT,DHOUT,CPAI,CPBI,CPCI)
674!$OMP+ PRIVATE(RMWI,CPDI,CPEI,CPFI,VX1,VY1,VZ1,VX2,VY2,VZ2,VX3,VY3,VZ3,VVX,VVY,VVZ,I1,I2,GAMA2)
675!$OMP+ PRIVATE(VFX,VFY,VFZ,SS_,ALPHA,RHOM,RUXM,RUYM,RUZM,REM,MASSFLOW,CPAM,CPBM,CPCM,RGASM,CPDM,CPEM,CPFM)
676!$OMP SECTIONS
677!$OMP SECTION
678 fel(1:nelt) = zero
679 elfvel(1:nelt)=zero
680!$OMP SECTION
681 eldmass(1:nel)=zero
682 eldehpy(1:nel)=zero
683!$OMP SECTION
684 itagp(1:npolh)=0
685 dm(1:npolh)=zero
686!$OMP SECTION
687 dq(1,1:npolh)=zero
688 dq(2,1:npolh)=zero
689!$OMP SECTION
690 dq(3,1:npolh)=zero
691 de(1:npolh)=zero
692!$OMP SECTION
693 dmini(1:npolh)=zero
694 dmcpa(1:npolh)=zero
695!$OMP SECTION
696 dmcpb(1:npolh)=zero
697 dmcpc(1:npolh)=zero
698!$OMP SECTION
699 dmrmw(1:npolh)=zero
700 dmcpd(1:npolh)=zero
701!$OMP SECTION
702 dmcpe(1:npolh)=zero
703 dmcpf(1:npolh)=zero
704!$OMP SECTION
705C----------------------
706C TRANSPORTS
707C----------------------
708 pu(1,1:npolh)=qpolh(1,1:npolh)/mpolh(1:npolh)
709 pu(2,1:npolh)=qpolh(2,1:npolh)/mpolh(1:npolh)
710 pu(3,1:npolh)=qpolh(3,1:npolh)/mpolh(1:npolh)
711!$OMP END SECTIONS
712
713
714! ================== !
715! FLUXES COMPUTATION !
716! ================== !
717!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
718 DO iel = 1, nelt
719 mass_flux(iel) = zero
720 momentum_flux_x(iel) = zero
721 momentum_flux_y(iel) = zero
722 momentum_flux_z(iel) = zero
723 energy_flux(iel) = zero
724 cpa_flux(iel) = zero
725 cpb_flux(iel) = zero
726 cpc_flux(iel) = zero
727 rgas_flux(iel) = zero
728 IF (ityp == 8) THEN
729 cpd_flux(iel) = zero
730 cpe_flux(iel) = zero
731 cpf_flux(iel) = zero
732 ENDIF
733 elfmass(iel) = zero
734 elfehpy(iel) = zero
735 elfvel(iel) = zero
736 vel(iel) = zero
737! SURFACE
738 area = elarea(iel)
739 ss(iel) = zero
740 vrel(iel) = zero
741! NORMAL VECTOR
742 nx = norm(1, iel)
743 ny = norm(2, iel)
744 nz = norm(3, iel)
745 IF (itagel(iel) > 0) THEN
746! INJECTOR
747 iinj = itagel(iel)
748 dmi = datainj(2, iinj) * area / datainj(1, iinj)
749 mass_flux(iel) = dmi
750 momentum_flux_x(iel) = -dmi * nx * datainj(3, iinj)
751 momentum_flux_y(iel) = -dmi * ny * datainj(3, iinj)
752 momentum_flux_z(iel) = -dmi * nz * datainj(3, iinj)
753 energy_flux(iel) = dmi * datainj(4, iinj)
754 cpa_flux(iel) = dmi * rbagjet(2, iinj)
755 cpb_flux(iel) = dmi * rbagjet(3, iinj)
756 cpc_flux(iel) = dmi * rbagjet(4, iinj)
757 rgas_flux(iel) = dmi * rbagjet(1, iinj)
758 IF (ityp == 8) THEN
759 cpd_flux(iel) = dmi * rbagjet(16, iinj)
760 cpe_flux(iel) = dmi * rbagjet(17, iinj)
761 cpf_flux(iel) = dmi * rbagjet(18, iinj)
762 ENDIF
763 elfmass(iel) = elfmass(iel) + dmi * dt1
764 elfehpy(iel) = elfehpy(iel) + dmi * datainj(4, iinj) * dt1
765 elfvel(iel) = -datainj(3, iinj)
766 ELSE IF (itagel(iel) < 0) THEN
767! VENT HOLE
768 ivent = -itagel(iel)
769 idef = ibaghol(1, ivent)
770 itvent = ibaghol(10, ivent)
771 aout = elsout(iel)
772 IF (fvdata(ifv)%TRI_TO_ELEM(1, iel) /= fvdata(ifv)%TRI_TO_ELEM(2, iel)) THEN
773 CALL abort()
774 ENDIF
775 i = fvdata(ifv)%TRI_TO_ELEM(1, iel)
776 p1 = ppolh(i)
777 rho1 = rpolh(i)
778 re1 = epolh(i) / pvolu(i)
779 gama1 = gpolh(i)
780 ux1 = pu(1,i)
781 uy1 = pu(2,i)
782 uz1 = pu(3,i)
783 uu = zero
784! LOCAL VELOCITY
785 IF (itvent == 1) THEN
786 uu = nx * ux1 + ny * uy1 + nz * uz1
787 rho2 = rho1
788 IF (p1 < pext) uu = zero
789! BERNOULLI
790 ELSEIF ((itvent == 2 .AND. p1 > zero) .OR. (itvent == 4 .AND. p1 >= pext)) THEN
791 vv = nx * ux1 + ny * uy1 + nz * uz1
792 vv = max(vv, zero)
793 eta = (gama1 - one) / gama1
794 ksi = one / eta
795 pstag = p1 * (one + half * eta * rho1 * vv * vv / p1)**ksi
796 pcrit = pstag * (two / (gama1+one))**ksi
797 p2 = max(pext, pcrit)
798 rho2 = rho1 * (p2 / p1)**(one / gama1)
799 uu = two * p1 * (one - (p2 / p1)**eta) / (rho1 * eta)
800 uu = max(uu, zero)
801 uu = sqrt(uu)
802 vmax = half * ((p1 - pext) * pvolu(i) / (gama1 - one))
803 . / max(em20, rho2 * aout * dt1 * gama1 * re1 / rho1)
804 vmax = max(vmax, zero)
805 uu = min(uu, vmax)
806! CHEMKIN
807 ELSE IF (itvent == 3) THEN
808 p2 = p1 - pext
809 ivdp = ibaghol(9, ivent)
810 fvdp = rbaghol(13, ivent)
811 uu = fvdp * get_u_func(ivdp, p2 * scalp, deri)
812 IF (p2 < zero) THEN
813 rho2 = rvolu(62)
814 ELSE
815 rho2 = rho1
816 ENDIF
817! BERNOULLI WITH FLOW-IN
818 ELSEIF (itvent == 4 .AND. p1 < pext) THEN
819 gamai = rvolu(1)
820 rhoi = rvolu(62)
821 eta = (gamai - one) / gamai
822 pcrit = pext * (two / (gamai + one))**(one / eta)
823 p2 = max(p1, pcrit)
824 rho2 = rhoi * (p2 / pext)**(one / gamai)
825 uu = two * pext * (one - (p2 / pext)**eta)/(rhoi * eta)
826 uu = max(uu, zero)
827 uu = sqrt(uu)
828 vmax = half*((pext - p1) * pvolu(i) / (gama1 - one))
829 . / max(em20, rho2 * aout * dt1 * rvolu(63))
830 vmax = max(vmax, zero)
831 uu = min(uu, vmax)
832 uu = -uu
833! GRAEFE
834 ELSE IF (itvent == 5) THEN
835 p2 = p1 - pext
836 rho2 = rho1
837 uu = two * p2 / rho2
838 uu = max(uu, zero)
839 uu = sqrt(uu)
840 ENDIF
841
842 IF (uu > zero .AND. idef == 1) THEN
843 dmout = uu * rho2 * aout
844 mass_flux(iel) = -dmout
845 momentum_flux_x(iel) = -dmout * ux1
846 momentum_flux_y(iel) = -dmout * uy1
847 momentum_flux_z(iel) = -dmout * uz1
848 dhout = dmout * gama1 * re1 / rho1
849 energy_flux(iel) = -dhout
850 cpai = cpapolh(i)
851 cpbi = cpbpolh(i)
852 cpci = cpcpolh(i)
853 rmwi = rmwpolh(i)
854 cpa_flux(iel) = -dmout * cpai
855 cpb_flux(iel) = -dmout * cpbi
856 cpc_flux(iel) = -dmout * cpci
857 rgas_flux(iel) = -dmout * rmwi
858 IF(ityp == 8) THEN
859 cpdi = cpdpolh(i)
860 cpei = cpepolh(i)
861 cpfi = cpfpolh(i)
862 cpd_flux(iel) = -dmout * cpdi
863 cpe_flux(iel) = -dmout * cpei
864 cpf_flux(iel) = -dmout * cpfi
865 ENDIF
866 !/TODO
867 !DMINI(I)=DMINI(I)-DMOUT
868 elfmass(iel)=elfmass(iel)-dmout*dt1
869 elfehpy(iel)=elfehpy(iel)-dhout*dt1
870 elfvel(iel)=elfvel(iel)+uu*area/elarea(iel)
871 eldmass(iel)=-dmout
872 eldehpy(iel)=-dhout
873 ENDIF
874C
875 IF (uu < zero .AND. idef == 1) THEN
876 dmi = -uu * rho2 * aout
877 mass_flux(iel) = dmi
878 momentum_flux_x(iel) = dmi * ux1
879 momentum_flux_y(iel) = dmi * uy1
880 momentum_flux_z(iel) = dmi * uz1
881 cpai = rvolu(7)
882 cpbi = rvolu(8)
883 cpci = rvolu(9)
884 rmwi = rvolu(10)
885 cpa_flux(iel) = dmi * cpai
886 cpb_flux(iel) = dmi * cpbi
887 cpc_flux(iel) = dmi * cpci
888 rgas_flux(iel) = dmi * rmwi
889 IF (ityp == 8) THEN
890 cpdi = rvolu(56)
891 cpei = rvolu(57)
892 cpfi = rvolu(58)
893 cpd_flux(iel) = dmi * cpdi
894 cpe_flux(iel) = dmi * cpei
895 cpf_flux(iel) = dmi * cpfi
896 ENDIF
897 energy_flux(iel) = dmi * rvolu(63)
898 elfmass(iel) = elfmass(iel) + dmi * dt1
899 elfehpy(iel) = elfehpy(iel) + dmi * dt1 * rvolu(63)
900 elfvel(iel) = elfvel(iel) + uu * area / elarea(iel)
901 eldmass(iel) = dmi
902 eldehpy(iel) = dmi * rvolu(63)
903 ENDIF
904 ELSE IF (iel > nel) THEN
905! TRIANGLE ON INTERNAL SURFACE TAKE POROSITY INTO ACCOUNT
906 IF (porosity(iel - nel) /= zero) THEN
907 n1=elem(1,iel)
908 n2=elem(2,iel)
909 n3=elem(3,iel)
910 IF (tagela(iel) > 0) THEN
911 vx1=vvv(1,n1)
912 vy1=vvv(2,n1)
913 vz1=vvv(3,n1)
914 vx2=vvv(1,n2)
915 vy2=vvv(2,n2)
916 vz2=vvv(3,n2)
917 vx3=vvv(1,n3)
918 vy3=vvv(2,n3)
919 vz3=vvv(3,n3)
920 ELSE
921 vx1=vvva(1,n1)
922 vy1=vvva(2,n1)
923 vz1=vvva(3,n1)
924 vx2=vvva(1,n2)
925 vy2=vvva(2,n2)
926 vz2=vvva(3,n2)
927 vx3=vvva(1,n3)
928 vy3=vvva(2,n3)
929 vz3=vvva(3,n3)
930 ENDIF
931 !grid velocity
932 vvx=third*(vx1+vx2+vx3)
933 vvy=third*(vy1+vy2+vy3)
934 vvz=third*(vz1+vz2+vz3)
935 !normal vector
936 nx=norm(1,iel)
937 ny=norm(2,iel)
938 nz=norm(3,iel)
939 !normal velocity
940 vel(iel)=nx*vvx+ny*vvy+nz*vvz
941 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
942 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
943 IF (i1 == i2) THEN
944 CALL abort()
945 ENDIF
946 area = elarea(iel) * porosity(iel - nel)
947 rho1 = rpolh(i1)
948 re1 = epolh(i1) / pvolu(i1)
949 ux1 = pu(1, i1)
950 uy1 = pu(2, i1)
951 uz1 = pu(3, i1)
952 gama1 = gpolh(i1)
953 rho2 = rpolh(i2)
954 re2 = epolh(i2) / pvolu(i2)
955 ux2 = pu(1, i2)
956 uy2 = pu(2, i2)
957 uz2 = pu(3, i2)
958 gama2 = gpolh(i2)
959 vfx = half * (ux1 + ux2)
960 vfy = half * (uy1 + uy2)
961 vfz = half * (uz1 + uz2)
962 vrel(iel) = (vfx - vvx)*(vfx - vvx) + (vfy - vvy)*(vfy - vvy) + (vfz - vvz)*(vfz - vvz)
963 ss_= nx * (vfx - vvx) + ny * (vfy - vvy) + nz * (vfz - vvz)
964 ss(iel) = elarea(iel)*abs(ss_)
965 IF (ss_ <= zero) THEN
966 alpha = zero
967 ELSE
968 alpha = one
969 ENDIF
970 rhom = alpha * rho1 + (one - alpha) * rho2
971 ruxm = alpha * rho1 * ux1 + (one - alpha) * rho2 * ux2
972 ruym = alpha * rho1 * uy1 + (one - alpha) * rho2 * uy2
973 ruzm = alpha * rho1 * uz1 + (one - alpha) * rho2 * uz2
974 rem = alpha * gama1 * re1 + (one - alpha) * gama2 * re2
975
976 massflow = rhom * ss_ * area
977 mass_flux(iel) = -massflow
978 momentum_flux_x(iel) = -ruxm * ss_ * area
979 momentum_flux_y(iel) = -ruym * ss_ * area
980 momentum_flux_z(iel) = -ruzm * ss_ * area
981 energy_flux(iel) = -rem * ss_ * area
982
983 elfmass(iel) = elfmass(iel) + massflow * dt1
984 elfvel(iel) = elfvel(iel) + ss_ * area
985
986 cpam = alpha * cpapolh(i1) + (one - alpha) * cpapolh(i2)
987 cpbm = alpha * cpbpolh(i1) + (one - alpha) * cpbpolh(i2)
988 cpcm = alpha * cpcpolh(i1) + (one - alpha) * cpcpolh(i2)
989 rgasm = alpha * rmwpolh(i1) + (one - alpha) * rmwpolh(i2)
990 cpa_flux(iel) = -massflow * cpam
991 cpb_flux(iel) = -massflow * cpbm
992 cpc_flux(iel) = -massflow * cpcm
993 rgas_flux(iel) = -massflow * rgasm
994 IF (ityp == 8) THEN
995 cpdm = alpha * cpdpolh(i1) + (one - alpha) * cpdpolh(i2)
996 cpem = alpha * cpepolh(i1) + (one - alpha) * cpepolh(i2)
997 cpfm = alpha * cpfpolh(i1) + (one - alpha) * cpfpolh(i2)
998 cpd_flux(iel) = -massflow * cpdm
999 cpe_flux(iel) = -massflow * cpem
1000 cpf_flux(iel) = -massflow * cpfm
1001 ENDIF
1002 ENDIF
1003 ENDIF
1004 ENDDO
1005!$OMP END DO
1006
1007!$OMP END PARALLEL
1008
1009!$OMP PARALLEL PRIVATE(I,MSINI,CPA,CPB,CPC,CPD,CPE,CPF,EFAC,TEMP0)
1010!$OMP+ PRIVATE(TEMP,CP,CV,RMW)
1011!$OMP+ PRIVATE(AL,DD,AD,GAMA,SSP,QX,VV,SS_,AREA)
1012! --- PRIVATE(AREAMAX)
1013!$OMP+ PRIVATE(ITYPL,J,IEL,COUNT)
1014
1015!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1016 DO i = 1, npolh
1017 msini = mpolh(i)
1018 cpapolh(i) = msini * cpapolh(i)
1019 cpbpolh(i) = msini * cpbpolh(i)
1020 cpcpolh(i) = msini * cpcpolh(i)
1021 rmwpolh(i) = msini * rmwpolh(i)
1022 IF (ityp == 8) THEN
1023 cpdpolh(i) = msini * cpdpolh(i)
1024 cpepolh(i) = msini * cpepolh(i)
1025 cpfpolh(i) = msini * cpfpolh(i)
1026 ENDIF
1027 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
1028 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
1029! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
1030 mpolh(i) = mpolh(i) + mass_flux(iel) * dt1
1031 qpolh(1, i) = qpolh(1, i) + momentum_flux_x(iel) * dt1
1032 qpolh(2, i) = qpolh(2, i) + momentum_flux_y(iel) * dt1
1033 qpolh(3, i) = qpolh(3, i) + momentum_flux_z(iel) * dt1
1034 epolh(i) = epolh(i) + energy_flux(iel) * dt1
1035 IF (mpolh(i) <= zero .OR. epolh(i) <= zero) THEN
1036 tpolh(i) = zero
1037 ELSE
1038 cpapolh(i) = cpapolh(i) + cpa_flux(iel) * dt1
1039 cpbpolh(i) = cpbpolh(i) + cpb_flux(iel) * dt1
1040 cpcpolh(i) = cpcpolh(i) + cpc_flux(iel) * dt1
1041 rmwpolh(i) = rmwpolh(i) + rgas_flux(iel) * dt1
1042 IF (ityp == 8) THEN
1043 cpdpolh(i) = cpdpolh(i) + cpd_flux(iel) * dt1
1044 cpepolh(i) = cpepolh(i) + cpe_flux(iel) * dt1
1045 cpfpolh(i) = cpfpolh(i) + cpf_flux(iel) * dt1
1046 ENDIF
1047 ENDIF
1048 ENDDO
1049 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
1050 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
1051! FINITE VOLUME I IS ON THE RIGHT SIDE OF ELEMENT IEL, HENCE A - SIGN
1052 mpolh(i) = mpolh(i) - mass_flux(iel) * dt1
1053 qpolh(1, i) = qpolh(1, i) - momentum_flux_x(iel) * dt1
1054 qpolh(2, i) = qpolh(2, i) - momentum_flux_y(iel) * dt1
1055 qpolh(3, i) = qpolh(3, i) - momentum_flux_z(iel) * dt1
1056 epolh(i) = epolh(i) - energy_flux(iel) * dt1
1057 IF (mpolh(i) <= zero .OR. epolh(i) <= zero) THEN
1058 tpolh(i) = zero
1059 ELSE
1060 cpapolh(i) = cpapolh(i) - cpa_flux(iel) * dt1
1061 cpbpolh(i) = cpbpolh(i) - cpb_flux(iel) * dt1
1062 cpcpolh(i) = cpcpolh(i) - cpc_flux(iel) * dt1
1063 rmwpolh(i) = rmwpolh(i) - rgas_flux(iel) * dt1
1064 IF (ityp == 8) THEN
1065 cpdpolh(i) = cpdpolh(i) - cpd_flux(iel) * dt1
1066 cpepolh(i) = cpepolh(i) - cpe_flux(iel) * dt1
1067 cpfpolh(i) = cpfpolh(i) - cpf_flux(iel) * dt1
1068 ENDIF
1069 ENDIF
1070 ENDDO
1071! TEMPERATURE COMPUTATION
1072 efac = epolh(i) / mpolh(i)
1073 cpapolh(i) = cpapolh(i) / mpolh(i)
1074 cpbpolh(i) = cpbpolh(i) / mpolh(i)
1075 cpcpolh(i) = cpcpolh(i) / mpolh(i)
1076 rmwpolh(i) = rmwpolh(i) / mpolh(i)
1077 IF (ityp == 8) THEN
1078 cpdpolh(i) = cpdpolh(i) / mpolh(i)
1079 cpepolh(i) = cpepolh(i) / mpolh(i)
1080 cpfpolh(i) = cpfpolh(i) / mpolh(i)
1081 ENDIF
1082 cpa = cpapolh(i)
1083 cpb = cpbpolh(i)
1084 cpc = cpcpolh(i)
1085 rmw = rmwpolh(i)
1086 IF (ityp == 8) THEN
1087 cpd = cpdpolh(i)
1088 cpe = cpepolh(i)
1089 cpf = cpfpolh(i)
1090 ENDIF
1091 temp0 = rvolu(25)
1092 CALL fvtemp(ityp , efac , cpa , cpb , cpc ,
1093 . cpd , cpe , cpf , rmw , temp0,
1094 . temp )
1095 cp = cpa + cpb * temp + cpc * temp * temp
1096 IF (ityp == 8) THEN
1097 cp = cp + cpd * temp**3 + cpf * temp**4
1098 IF (cpe /= zero .AND. temp > zero) THEN
1099 cp = cp + cpe / (temp * temp)
1100 ENDIF
1101 ENDIF
1102 cv = cp - rmw
1103 rmwpolh(i) = rmw
1104 gpolh(i) = cp / cv
1105 tpolh(i) = temp
1106! PRESSURE AND MASS DENSITY
1107 rpolh(i) = mpolh(i) / pvolu(i)
1108 ppolh(i) = rpolh(i) * rmw * temp
1109 ENDDO
1110!$OMP END DO
1111
1112C---------------------------------------------------
1113C QA, QB AND STABILITY
1114C---------------------------------------------------
1115!$OMP SINGLE
1116
1117 ! /DT/FVMBAG/[ID_DT_OPTION]
1118 ! ID_DT_OPTION = 0 obsolete
1119 ! ID_DT_OPTION = 1 (legacy method), apply to all FVMBAG1
1120 ! ID_DT_OPTION = 2 (2022.3), apply to given FVMBAG1 or FVMBAG2.
1121 ! in this case characteristic length is calculated depending on L_TYPE flag
1122 ! and smoothing is applied : dt <- min(dt, dt_old*(1+lambda))
1123 !--------!-------------------------------------------------------------------!
1124 ! VALUE ! DESCRIPTION !
1125 !------------!--------!-------------------------------------------------------------------!
1126 ! ! 0,1 ! Constant (Starter minimum characteristic length) !
1127 ! L_TYPE ! 2 ! VOLUME^1/3 at current time !
1128 ! ! 3 ! VOLUME / SMAX at current time !
1129 !------------!--------!-------------------------------------------------------------------!
1130
1131
1132 !------------------------------!
1133 idt_fvmbag = fvdata(ifv)%ID_DT_OPTION !IDTMIN(52) ! /DT/FVMBAG/[IDT_FVMBAG]
1134
1135 qa=rvolu(46)
1136 qb=rvolu(47)
1137 dtx=ep30
1138 phdt=0
1139
1140!$OMP END SINGLE
1141
1142 !------------------------------!
1143 ! --- CHARACTERISTIC LENGTH ---!
1144 !------------------------------!
1145 idt_fvmbag = fvdata(ifv)%ID_DT_OPTION ! /DT/FVMBAG/[IDT_FVMBAG]
1146
1147 SELECT CASE(idt_fvmbag)
1148
1149 ! /DT/FVMBAG/1
1150 CASE(0,1)
1151 ! ------------
1152 param_l_type = 0
1153!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1154 DO i=1,npolh
1155 IF (ibpolh(i) == 0) THEN
1156 dls(i)=dlh !POLYHEDRON CONSTANT CHARACTERISTIC LENGTH COMPUTED IN THE STARTER
1157 ELSE
1158 dls(i)=pvolu(i)**third !ADDITIONAL BRICKS
1159 ENDIF
1160 ENDDO
1161!$OMP END DO
1162
1163 ! /DT/FVMBAG/2
1164 CASE(2)
1165 ! ------------
1166 param_l_type = fvdata(ifv)%L_TYPE
1167
1168 SELECT CASE(param_l_type)
1169 CASE(0,1) !dlh = starter constant
1170 ! /DT/FVMBAG/2 + L_TYPE=0,1
1171!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1172 DO i=1,npolh
1173 dls(i)=dlh
1174 ENDDO
1175!$OMP END DO
1176 CASE(2) ! dl=VOL**1/3
1177 ! /DT/FVMBAG/2 + L_TYPE=2
1178!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1179 DO i=1,npolh
1180 dls(i)=pvolu(i)**third
1181 ENDDO
1182!$OMP END DO
1183 CASE(3) ! dl=VOL**1/3
1184 ! /DT/FVMBAG/2 + L_TYPE=3
1185!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1186 DO i=1,npolh
1187 dls(i)=pvolu(i)
1188 ENDDO
1189!$OMP END DO
1190 END select!PARAM_L_TYPE
1191
1192 END select!IDT_FVMBAG
1193
1194
1195 !------------------------------!
1196 ! --- TIME STEP ---!
1197 !------------------------------!
1198 SELECT CASE (idt_fvmbag)
1199
1200 CASE(0,1)
1201!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1202 DO i=1,npolh
1203 IF(mpolh(i)<=zero.OR.epolh(i)<=zero.OR.gpolh(i)<=one) THEN
1204 dtpolh(i)=ep30
1205 qvisc(i) =zero
1206 ELSE
1207C
1208 al=pvolu(i)**third
1209 dd=dm(i)/mpolh(i)
1210 ad=max(zero,dd)
1211 gama=gpolh(i)
1212 ssp=sqrt((gama-one)*gama*epolh(i)/mpolh(i))
1213 qvisc(i)=rpolh(i)*ad*al*(qa*qa*al*ad+qb*ssp)
1214C
1215 qx=qb*ssp+al*qa*qa*ad
1216 ssp=qx+sqrt(qx*qx+ssp*ssp)
1217 vv=sqrt(pu(1,i)**2+pu(2,i)**2+pu(3,i)**2)
1218 dtpolh(i)=cfl_coef*dls(i)/(ssp+vv)
1219 ENDIF
1220 ENDDO
1221!$OMP END DO
1222
1223 CASE(2) !time step computed from a characteristic length
1224
1225 param_l_type = fvdata(ifv)%L_TYPE !characteristic length formulation
1226
1227 IF(param_l_type /= 3)THEN !L_TYPE=0/1 & 2
1228!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1229 DO i=1,npolh
1230 IF(mpolh(i) <= zero.OR.epolh(i) <= zero.OR.gpolh(i) <= one) THEN
1231 dtpolh(i)=ep30
1232 qvisc(i) =zero
1233 ELSE
1234 ! 1.calculating Grid Velocity for each polyedron
1235 count = 0
1236 vv = zero
1237 !loop over polygons
1238 !AREAMAX=ZERO
1239 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
1240 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
1241 area = elarea(iel)
1242! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
1243 !AREAMAX=MAX(AREAMAX,AREA)
1244 vv = vv+sqrt(vrel(iel))
1245 count=count+1
1246 ENDDO
1247 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
1248 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
1249! FINITE VOLUME I IS ON THE RIGHT SIDE OF ELEMENT IEL, HENCE A - SIGN
1250 vv = vv+sqrt(vrel(iel))
1251 count=count+1
1252 ENDDO
1253 al=pvolu(i)**third
1254 dd=dm(i)/mpolh(i)
1255 ad=max(zero,dd)
1256 gama=gpolh(i)
1257 ssp=sqrt((gama-one)*gama*epolh(i)/mpolh(i))
1258 qvisc(i)=rpolh(i)*ad*al*(qa*qa*al*ad+qb*ssp)
1259 qx=qb*ssp+al*qa*qa*ad
1260 ssp=qx+sqrt(qx*qx+ssp*ssp)
1261 ! 3.relative velocty (PU:material velocity : PW:grid velocity)
1262 IF(count > 0)THEN
1263 vv = vv / float(count)
1264 ENDIF
1265 ! 4.time step
1266 dtpolh(i)=cfl_coef*dls(i)/(ssp+vv)
1267 ENDIF
1268 ENDDO
1269!$OMP END DO
1270
1271 ELSE !L_TYPE=3 ! check all faces and compute normal relative velocities
1272!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1273 DO i = 1, npolh
1274 ss_=zero
1275 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
1276 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
1277! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
1278 ss_=max(ss_,ss(iel)) ! <|u-w|,n>
1279 ENDDO
1280 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
1281 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
1282! FINITE VOLUME I IS ON THE RIGHT SIDE OF ELEMENT IEL, HENCE A - SIGN
1283 ss_=max(ss_,ss(iel)) ! <|u-w|,n>
1284 ENDDO
1285
1286 al=pvolu(i)**third
1287 dd=dm(i)/mpolh(i)
1288 ad=max(zero,dd)
1289 gama=gpolh(i)
1290 ssp=sqrt((gama-one)*gama*epolh(i)/mpolh(i))
1291 qvisc(i)=rpolh(i)*ad*al*(qa*qa*al*ad+qb*ssp)
1292 qx=qb*ssp+al*qa*qa*ad
1293 ssp=qx+sqrt(qx*qx+ssp*ssp)
1294 ! 3.relative velocty (PU:material velocity : PW:grid velocity)
1295 ! SS_ <- max( area*(u-w).n )
1296 ! u-w = 0 unless in case of porosity
1297 ! 4.time step
1298 dtpolh(i)=cfl_coef*dls(i)/(ssp+ss_)
1299 ENDDO
1300!$OMP END DO
1301 ENDIF
1302
1303 END SELECT
1304
1305
1306!$OMP END PARALLEL
1307
1308 DO i=1,npolh
1309 IF (dtpolh(i)<dtx) THEN
1310 dtx=dtpolh(i)
1311 phdt=i
1312 ENDIF
1313 ENDDO
1314
1315 !optional time step growth limitor (/DT/FVMBAG/2, 3rd parameter)
1316 IF(ncycle > 0 .AND. idt_fvmbag == 2)THEN
1317 lambda = fvdata(ifv)%LAMBDA
1318 IF(lambda > zero)THEN
1319 dtold = fvdata(ifv)%DTOLD
1320 IF(dtold > zero)THEN
1321 dtx = min( dtx, dtold*(one+lambda))
1322 ENDIF
1323 ENDIF
1324 ENDIF
1325 fvdata(ifv)%DTOLD = dtx
1326
1327 IF (ilvout >= 1.AND.ipri==0) THEN
1328 WRITE(iout,'(A,I10,A,E12.4,A,I10)')' ** MONITORED VOLUME ID: ',ivolu(1),' TIME STEP:',dtx,' FINITE VOLUME:',idpolh(phdt)
1329 ENDIF
1330
1331C
1332 IF (dtx<dt2) THEN
1333 dt2=dtx
1334 nelts=ivolu(1)
1335 itypts=52
1336 ENDIF
1337C---------------------------------------------------
1338C PRESSURE FORCES AND WORK COMPUTATION
1339C---------------------------------------------------
1340!$OMP PARALLEL PRIVATE(N1,N2,N3,IEL,NX,NY,NZ)
1341!$OMP+ PRIVATE(VX1,VX2,VX3,VY1,VY2,VY3,VZ1,VZ2,VZ3)
1342!$OMP+ PRIVATE(VVX,VVY,VVZ,I1,I2,AREA,IVENT,IDEF,IVENTYP,AOUT,PP,P1,P2,PMEAN,AREA_OLD,AREA1,SS_)
1343!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1344 DO iel=1,nelt
1345 n1=elem(1,iel)
1346 n2=elem(2,iel)
1347 n3=elem(3,iel)
1348 vx1=vvv(1,n1)
1349 vy1=vvv(2,n1)
1350 vz1=vvv(3,n1)
1351 vx2=vvv(1,n2)
1352 vy2=vvv(2,n2)
1353 vz2=vvv(3,n2)
1354 vx3=vvv(1,n3)
1355 vy3=vvv(2,n3)
1356 vz3=vvv(3,n3)
1357 vvx=third*(vx1+vx2+vx3)
1358 vvy=third*(vy1+vy2+vy3)
1359 vvz=third*(vz1+vz2+vz3)
1360 nx=norm(1,iel)
1361 ny=norm(2,iel)
1362 nz=norm(3,iel)
1363 vel(iel)=nx*vvx+ny*vvy+nz*vvz
1364 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
1365 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
1366 area = elarea(iel)
1367 momentum_flux_x(iel) = zero
1368 momentum_flux_y(iel) = zero
1369 momentum_flux_z(iel) = zero
1370 energy_flux(iel) = zero
1371 fel(iel) = zero
1372 IF (i1 == i2) THEN
1373! EXTERNAL TRIANGLE
1374 IF (itagel(iel) < 0) THEN
1375! VENT HOLE
1376 ivent = -itagel(iel)
1377 idef = ibaghol(1, ivent)
1378 iventyp = ibaghol(13, ivent)
1379 IF (idef == 1 .AND. iventyp == 0) THEN
1380 aout = elsout(iel)
1381 area = area - aout
1382 pp = pext
1383 fel(iel) = pp * aout
1384 momentum_flux_x(iel) = -pp * aout * nx
1385 momentum_flux_y(iel) = -pp * aout * ny
1386 momentum_flux_z(iel) = -pp * aout * nz
1387 energy_flux(iel) = -pp * vel(iel) * aout
1388 ENDIF
1389 ENDIF
1390 pp = ppolh(i1) + qvisc(i1)
1391 fel(iel) = fel(iel) + pp * area
1392 momentum_flux_x(iel) = momentum_flux_x(iel) - pp * area * nx
1393 momentum_flux_y(iel) = momentum_flux_y(iel) - pp * area * ny
1394 momentum_flux_z(iel) = momentum_flux_z(iel) - pp * area * nz
1395 energy_flux(iel) = energy_flux(iel) - pp * vel(iel) * area
1396 energy_flux(iel) = energy_flux(iel) - rvolu(19) * area * (tpolh(i1) - rvolu(25))
1397 ELSE
1398! INTERNAL TRIANGLE, SEPARATES 2 FINITE VOLUMES
1399 p1 = ppolh(i1) + qvisc(i1)
1400 p2 = ppolh(i2) + qvisc(i2)
1401 pmean = half * (p1 + p2)
1402 area_old = area
1403 area = area * porosity(iel - nel)
1404 area1 = area_old - area
1405 ss_ = nx * vvx + ny * vvy + nz * vvz
1406 IF (itagel(iel) > 0) THEN
1407! INJECTOR
1408 fel(iel) = (p1 - p2) * area1
1409 momentum_flux_x(iel) = -pmean * nx * area
1410 momentum_flux_y(iel) = -pmean * ny * area
1411 momentum_flux_z(iel) = -pmean * nz * area
1412 energy_flux(iel) = -pmean * ss_ * area
1413 ELSE
1414 fel(iel) = (p1 - p2) * area1
1415 momentum_flux_x(iel) = -pmean * nx * area
1416 momentum_flux_y(iel) = -pmean * ny * area
1417 momentum_flux_z(iel) = -pmean * nz * area
1418 energy_flux(iel) = -pmean * ss_ * area
1419 ENDIF
1420 ENDIF
1421 ENDDO
1422!$OMP END DO
1423!$OMP END PARALLEL
1424
1425C---------------------------------------------------
1426C NEW QUANTITIES IN POLYEDRA
1427C---------------------------------------------------
1428!$OMP PARALLEL PRIVATE(I, J, IEL)
1429!$OMP+ PRIVATE(NX, NY, NZ, N1, N2, N3, VX1, VX2, VX3, VY1, VY2, VY3, VZ1, VZ2, VZ3)
1430!$OMP+ PRIVATE(VVX, VVY, VVZ, I1, I2, PP, AREA_OLD, AREA, AREA1, SS_)
1431!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1432 DO i = 1, npolh
1433 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%MINUS_SIZE
1434 iel = fvdata(ifv)%ELEM_TO_TRI(i)%MINUS(j)
1435 nx=norm(1,iel)
1436 ny=norm(2,iel)
1437 nz=norm(3,iel)
1438 n1=elem(1,iel)
1439 n2=elem(2,iel)
1440 n3=elem(3,iel)
1441 vx1=vvv(1,n1)
1442 vy1=vvv(2,n1)
1443 vz1=vvv(3,n1)
1444 vx2=vvv(1,n2)
1445 vy2=vvv(2,n2)
1446 vz2=vvv(3,n2)
1447 vx3=vvv(1,n3)
1448 vy3=vvv(2,n3)
1449 vz3=vvv(3,n3)
1450 vvx=third*(vx1+vx2+vx3)
1451 vvy=third*(vy1+vy2+vy3)
1452 vvz=third*(vz1+vz2+vz3)
1453 area = elarea(iel)
1454 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
1455 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
1456 IF (i1 /= i) THEN
1457 CALL abort()
1458 ENDIF
1459 pp = ppolh(i) + qvisc(i)
1460 IF (iel > nel) THEN
1461 area_old = area
1462 area = area * porosity(iel - nel)
1463 area1 = area_old - area
1464 ELSE
1465 area1 = zero
1466 ENDIF
1467
1468 ss_ = nx * vvx + ny * vvy + nz * vvz
1469
1470! FINITE VOLUME I IS ON THE LEFT SIDE OF ELEMENT IEL, HENCE A + SIGN
1471 qpolh(1, i) = qpolh(1, i) + (momentum_flux_x(iel) - pp * nx * area1) * dt1
1472 qpolh(2, i) = qpolh(2, i) + (momentum_flux_y(iel) - pp * ny * area1) * dt1
1473 qpolh(3, i) = qpolh(3, i) + (momentum_flux_z(iel) - pp * nz * area1) * dt1
1474 IF (itagel(iel) > 0) THEN
1475 epolh(i) = epolh(i) + (energy_flux(iel) + pp * vel(iel) * area1) * dt1
1476 ELSE
1477 epolh(i) = epolh(i) + (energy_flux(iel) - pp * ss_ * area1) * dt1
1478 ENDIF
1479 ENDDO
1480 DO j = 1, fvdata(ifv)%ELEM_TO_TRI(i)%PLUS_SIZE
1481 iel = fvdata(ifv)%ELEM_TO_TRI(i)%PLUS(j)
1482 nx=norm(1,iel)
1483 ny=norm(2,iel)
1484 nz=norm(3,iel)
1485 n1=elem(1,iel)
1486 n2=elem(2,iel)
1487 n3=elem(3,iel)
1488 vx1=vvv(1,n1)
1489 vy1=vvv(2,n1)
1490 vz1=vvv(3,n1)
1491 vx2=vvv(1,n2)
1492 vy2=vvv(2,n2)
1493 vz2=vvv(3,n2)
1494 vx3=vvv(1,n3)
1495 vy3=vvv(2,n3)
1496 vz3=vvv(3,n3)
1497 vvx=third*(vx1+vx2+vx3)
1498 vvy=third*(vy1+vy2+vy3)
1499 vvz=third*(vz1+vz2+vz3)
1500 area = elarea(iel)
1501 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
1502 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
1503 IF (i2 /= i) THEN
1504 CALL abort()
1505 ENDIF
1506 pp = ppolh(i) + qvisc(i)
1507 IF (iel > nel) THEN
1508 area_old = area
1509 area = area * porosity(iel - nel)
1510 area1 = area_old - area
1511 ELSE
1512 area1 = zero
1513 ENDIF
1514
1515 ss_ = nx * vvx + ny * vvy + nz * vvz
1516
1517
1518! FINITE VOLUME I IS ON THE RIGHT SIDE OF ELEMENT IEL, HENCE A - SIGN
1519 qpolh(1, i) = qpolh(1, i) - (momentum_flux_x(iel) - pp * nx * area1) * dt1
1520 qpolh(2, i) = qpolh(2, i) - (momentum_flux_y(iel) - pp * ny * area1) * dt1
1521 qpolh(3, i) = qpolh(3, i) - (momentum_flux_z(iel) - pp * nz * area1) * dt1
1522 IF (itagel(iel) > 0) THEN
1523 epolh(i) = epolh(i) - (energy_flux(iel) - pp * vel(iel) * area1) * dt1
1524 ELSE
1525 epolh(i) = epolh(i) - (energy_flux(iel) - pp * ss_ * area1) * dt1
1526 ENDIF
1527 ENDDO
1528 ENDDO
1529!$OMP END DO
1530!$OMP END PARALLEL
1531
1532 ENDIF
1533
1534C---------------------------------------------------
1535C PRINT OUTS
1536C---------------------------------------------------
1537 250 IF (ilvout >= 1 .AND. ipri == 0) THEN
1538 WRITE(iout,*)
1539 WRITE(iout,'(A25,I10,A31)')' ** MONITORED VOLUME ID: ',ivolu(1),' - FINITE VOLUME MESH UPDATE **'
1540 WRITE(iout,'(A,I8)')' NUMBER OF FINITE VOLUMES : ',npolh
1541 WRITE(iout,'(A,G16.9,A17,G16.9,A1)')' SUM VOLUME FINITE VOLUMES : ',volph,' (VOLUME AIRBAG: ',volg,')'
1542 WRITE(iout,'(A,G16.9,A17,G16.9,A1)')' SUM AREA SURFACE POLYGONS : ',areap,' (AREA AIRBAG : ',areael,')'
1543 WRITE(iout,*)
1544 ENDIF
1545C---------------------------------------------------
1546C TIME HISTORY
1547C---------------------------------------------------
1548 amtot=zero
1549 pmean=zero
1550 enint=zero
1551 fsav(1:17)=zero
1552 DO i=1,npolh
1553 amtot=amtot+mpolh(i)
1554 enint=enint+epolh(i)
1555 pmean=pmean+ppolh(i)*pvolu(i)
1556 gama=gpolh(i)
1557 rmw=rmwpolh(i)
1558 temp=tpolh(i)
1559 cpa=cpapolh(i)
1560 cpb=cpbpolh(i)
1561 cpc=cpcpolh(i)
1562 cp=cpa+cpb*temp+cpc*temp*temp
1563 IF(ityp==8) THEN
1564 cpd=cpdpolh(i)
1565 cpe=cpepolh(i)
1566 cpf=cpfpolh(i)
1567 cp=cp+cpd*temp**3+cpf*temp**4
1568 IF(temp > zero) cp=cp+cpe/(temp*temp)
1569 ENDIF
1570 cv=cp-rmw
1571 fsav(5) =fsav(5) +mpolh(i)*temp
1572 fsav(10)=fsav(10)+mpolh(i)*cp
1573 fsav(11)=fsav(11)+mpolh(i)*cv
1574 fsav(12)=fsav(12)+mpolh(i)*gama
1575 ENDDO
1576C PRESSURE DISPERSION
1577 pdisp = zero
1578 pmean2 = pmean / volph
1579 DO i = 1, npolh
1580 pdisp = pdisp + pvolu(i) * (ppolh(i) - pmean2)**2
1581 ENDDO
1582 pdisp = sqrt(pdisp / volph) / pmean2
1583 fsav(19) = pdisp
1584 injection_started = .false.
1585 DO iinj = 1, njet
1586 injection_started = injection_started .OR. (datainj(2, iinj) > zero)
1587 ENDDO
1588 IF (.NOT. injection_started) THEN
1589 pdisp = zero
1590 ENDIF
1591
1592C
1593 fsav(1)=amtot
1594 fsav(2)=volph
1595 pres_av = pmean / volph
1596 fsav(3) = pres_av
1597 fsav(4)=areap
1598 fsav(7)=zero
1599 dmout =zero
1600 dhout =zero
1601 cp_av = fsav(10) / amtot
1602 fsav(10) = cp_av
1603 cv_av = fsav(11) / amtot
1604 fsav(11) = cv_av
1605 fsav(12) = cp_av / cv_av
1606 fsav(14) = npolh
1607 rgas_av = cp_av - cv_av
1608 rho_av = amtot / volph
1609 fsav(5) = fsav(5) / amtot ! sum(mass[i]*T[i])/total_mass
1610
1611 IF(ivolu(39) == 0) THEN
1612 fsav(6) =zero
1613 fsav(13)=zero
1614 ELSE
1615 fsav(6) =aoutot
1616 DO iel=1,nel
1617 IF (itagel(iel)<0) THEN
1618 ivent=-itagel(iel)
1619 rbaghol(18,ivent)=rbaghol(18,ivent)+elfvel(iel)*elsout(iel)
1620 rbaghol(19,ivent)=rbaghol(19,ivent)-elfmass(iel)
1621 rbaghol(20,ivent)=rbaghol(20,ivent)-elfehpy(iel)
1622 rbaghol(21,ivent)=rbaghol(21,ivent)+eldmass(iel)
1623 rbaghol(22,ivent)=rbaghol(22,ivent)+eldehpy(iel)
1624 ENDIF
1625 ENDDO
1626 DO ivent=1,nvent
1627 sventot=rbaghol(16,ivent)+rbaghol(17,ivent)
1628 IF(sventot <= zero) cycle
1629 rbaghol(18,ivent)=rbaghol(18,ivent)/sventot
1630 ENDDO
1631 DO ivent=1,nvent
1632 fsav(7)=fsav(7)+(rbaghol(16,ivent)+rbaghol(17,ivent))*rbaghol(18,ivent)
1633 dmout =dmout+rbaghol(21,ivent)
1634 dhout =dhout+rbaghol(22,ivent)
1635 ENDDO
1636 fsav(7) =fsav(7) /max(aoutot,em20)
1637 fsav(13)=dtx ! not DTPOLH(PHDT) since DTX may be updated with smoothing factor
1638
1639 IF(ityp == 8) THEN
1640 CALL fvinjt8_1(njet, ibagjet , rbagjet , igeo, geo, pm,
1641 . ivolu, rvolu, dmout, dhout)
1642 ttf =rvolu(60)
1643 IF (ivolu(74) == 0) THEN
1644 up_switch = tt-ttf >= rvolu(70) .OR. npolh <= ivolu(37)
1645 ENDIF
1646 IF (ivolu(74) == 1) THEN
1647C AUTOMATIC SWITCH TO UNIFORM PRESSURE WHEN DISPERSION OF PRESSURE IS LOW
1648 up_switch = ((pdisp < pdisp_old) .AND. (pdisp < rvolu(73)))
1649 up_switch = up_switch .OR. tt-ttf >= rvolu(70)
1650 ENDIF
1651 IF (ivolu(74) == 2) THEN
1652C AUTOMATIC SWITCH TO UNIFORM PRESSURE WHEN DISPERSION OF PRESSURE IS LOW
1653 up_switch = ((pdisp < pdisp_old) .AND. (pdisp < rvolu(73)))
1654 up_switch = up_switch .AND. tt-ttf >= rvolu(70)
1655 ENDIF
1656
1657 IF (up_switch) THEN
1658 rvolu(12)=fsav(3) ! PRESSION
1659 rvolu(13)=fsav(5) ! TEMPERATURE
1660 rvolu(16)=fsav(2) ! VOLUME
1661 rvolu(20)=fsav(1) ! MASSE
1662 ENDIF
1663 ENDIF
1664
1665 ENDIF
1666
1667 DO iel=1,nelt
1668 IF (itagel(iel) > 0) THEN
1669 fsav(15)=fsav(15)+elfmass(iel)
1670 fsav(16)=fsav(16)+elfehpy(iel)
1671 ENDIF
1672 ENDDO
1673 fsav(17)=enint
1674
1675
1676C---------------------------------------------------
1677C NODAL FLUID VELOCITIES
1678C---------------------------------------------------
1679 IF( (tt>=tanim .AND. tt<=tanim_stop) .OR.tt>=toutp.OR. nbgauge > 0 .OR.
1680 . (manim>=4.AND.manim<=15) ) THEN
1681C COMPUTATIONS FOR ANIM, ALSO NEEDED IF THERE ARE GAUGES.
1682 prepare_anim = 1
1683 uel(1:3,1:nelt)=zero
1684 rhoel(1:nelt)=zero
1685 tkel(1:nelt) =zero
1686 sspel(1:nelt)=zero
1687 ELSE
1688 prepare_anim = 0
1689 ENDIF
1690
1691C PEL(1:NELT)=ZERO
1692
1693!$OMP PARALLEL PRIVATE(I,J,JJ,K,KK,IEL,JEL,II,NV)
1694!$OMP+ PRIVATE(NX,NY,NZ,NRX,NRY,NRZ)
1695!$OMP+ PRIVATE(FAC,AREA,I1,I2)
1696!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1697 DO iel = 1, nelt
1698 area = elarea(iel)
1699 i1 = fvdata(ifv)%TRI_TO_ELEM(1, iel)
1700 i2 = fvdata(ifv)%TRI_TO_ELEM(2, iel)
1701 pel(iel) = half * (ppolh(i1) + ppolh(i2))
1702 IF (prepare_anim == 1) THEN
1703 uel(1, iel) = half * (pu(1, i1) + pu(1, i2))
1704 uel(2, iel) = half * (pu(2, i1) + pu(2, i2))
1705 uel(3, iel) = half * (pu(3, i1) + pu(3, i2))
1706 rhoel(iel) = half * (rpolh(i1) + rpolh(i2))
1707 tkel(iel) = half * (tpolh(i1) + tpolh(i2))
1708 sspel(iel) = half * (sqrt(gpolh(i1)*ppolh(i1)/rpolh(i1)) + sqrt(gpolh(i2)*ppolh(i2)/rpolh(i2)))
1709 ENDIF
1710 ENDDO
1711!$OMP END DO
1712
1713C------------------
1714C VENT HOLES
1715C------------------
1716 IF(prepare_anim == 1) THEN
1717!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1718 DO iel=1,nel
1719 IF(itagel(iel)<0) THEN
1720 uel(1,iel)=elfvel(iel)*norm(1,iel)
1721 uel(2,iel)=elfvel(iel)*norm(2,iel)
1722 uel(3,iel)=elfvel(iel)*norm(3,iel)
1723 ENDIF
1724 ENDDO
1725!$OMP END DO NOWAIT
1726!$OMP SINGLE
1727 u(1:3,1:nnt)=zero
1728 rho(1:nnt) =zero
1729 tk(1:nnt) =zero
1730 sspk(1:nnt) =zero
1731 vola(1:nna)=zero
1732 pa(1:nna)=zero
1733 rhoa(1:nna)=zero
1734 tka(1:nna)=zero
1735 sspka(1:nna)=zero
1736 ua(1:3,1:nna)=zero
1737!$OMP END SINGLE
1738 ENDIF
1739
1740!$OMP DO SCHEDULE(DYNAMIC,MVSIZ)
1741 DO i=1,nnt
1742 p(i)=zero
1743 ENDDO
1744!$OMP END DO
1745!$OMP END PARALLEL
1746
1747 DO iel=1,nelt
1748 n1=elem(1,iel)
1749 n2=elem(2,iel)
1750 n3=elem(3,iel)
1751 area1=nodarea(n1)
1752 area2=nodarea(n2)
1753 area3=nodarea(n3)
1754 fac1=third*elarea(iel)/area1
1755 fac2=third*elarea(iel)/area2
1756 fac3=third*elarea(iel)/area3
1757 p(n1) =p(n1) +fac1*pel(iel)
1758 p(n2) =p(n2) +fac2*pel(iel)
1759 p(n3) =p(n3) +fac3*pel(iel)
1760 ENDDO
1761
1762 IF(prepare_anim == 1) THEN
1763 DO iel=1,nelt
1764 n1=elem(1,iel)
1765 n2=elem(2,iel)
1766 n3=elem(3,iel)
1767 area1=nodarea(n1)
1768 area2=nodarea(n2)
1769 area3=nodarea(n3)
1770 fac1=third*elarea(iel)/area1
1771 fac2=third*elarea(iel)/area2
1772 fac3=third*elarea(iel)/area3
1773 u(1,n1)=u(1,n1)+fac1*uel(1,iel)
1774 u(2,n1)=u(2,n1)+fac1*uel(2,iel)
1775 u(3,n1)=u(3,n1)+fac1*uel(3,iel)
1776 u(1,n2)=u(1,n2)+fac2*uel(1,iel)
1777 u(2,n2)=u(2,n2)+fac2*uel(2,iel)
1778 u(3,n2)=u(3,n2)+fac2*uel(3,iel)
1779 u(1,n3)=u(1,n3)+fac3*uel(1,iel)
1780 u(2,n3)=u(2,n3)+fac3*uel(2,iel)
1781 u(3,n3)=u(3,n3)+fac3*uel(3,iel)
1782 rho(n1)=rho(n1)+fac1*rhoel(iel)
1783 rho(n2)=rho(n2)+fac2*rhoel(iel)
1784 rho(n3)=rho(n3)+fac3*rhoel(iel)
1785 tk(n1)=tk(n1)+fac1*tkel(iel)
1786 tk(n2)=tk(n2)+fac2*tkel(iel)
1787 tk(n3)=tk(n3)+fac3*tkel(iel)
1788 sspk(n1)=sspk(n1)+fac1*sspel(iel)
1789 sspk(n2)=sspk(n2)+fac2*sspel(iel)
1790 sspk(n3)=sspk(n3)+fac3*sspel(iel)
1791 ENDDO
1792
1793 DO i=1,npolh
1794 IF (ibpolh(i)/=0) THEN
1795 temp=tpolh(i)
1796 ii=iabs(ibpolh(i))
1797 IF(brna(1,ii)==brna(4,ii).AND.brna(5,ii)==brna(8,ii))THEN
1798C PENTAHEDRON
1799 DO k=1,6
1800 j=k
1801 IF(k>=4) j=j+1
1802 jj=brna(j,ii)
1803 vola(jj)=vola(jj)+one_over_6*pvolu(i)
1804 pa(jj)=pa(jj)+one_over_6*ppolh(i)*pvolu(i)
1805 rhoa(jj)=rhoa(jj)+one_over_6*rpolh(i)*pvolu(i)
1806 tka(jj)=tka(jj)+one_over_6*temp*pvolu(i)
1807 sspka(jj)=sspka(jj)+one_over_6*sqrt(gpolh(i)*ppolh(i)/rpolh(i))*pvolu(i)
1808 ua(1,jj)=ua(1,jj)+one_over_6*pu(1,i)*pvolu(i)
1809 ua(2,jj)=ua(2,jj)+one_over_6*pu(2,i)*pvolu(i)
1810 ua(3,jj)=ua(3,jj)+one_over_6*pu(3,i)*pvolu(i)
1811 ENDDO
1812 ELSEIF(brna(5,ii)==brna(8,ii).AND.
1813 . brna(6,ii)==brna(8,ii).AND.
1814 . brna(7,ii)==brna(8,ii))THEN
1815C PYRAMID
1816 DO j=1,5
1817 jj=brna(j,ii)
1818 vola(jj)=vola(jj)+one_fifth*pvolu(i)
1819 pa(jj)=pa(jj)+one_fifth*ppolh(i)*pvolu(i)
1820 rhoa(jj)=rhoa(jj)+one_fifth*rpolh(i)*pvolu(i)
1821 tka(jj)=tka(jj)+one_fifth*temp*pvolu(i)
1822 sspka(jj)=sspka(jj)+one_fifth*sqrt(gpolh(i)*ppolh(i)/rpolh(i))*pvolu(i)
1823 ua(1,jj)=ua(1,jj)+one_fifth*pu(1,i)*pvolu(i)
1824 ua(2,jj)=ua(2,jj)+one_fifth*pu(2,i)*pvolu(i)
1825 ua(3,jj)=ua(3,jj)+one_fifth*pu(3,i)*pvolu(i)
1826 ENDDO
1827 ELSE
1828C HEXAHEDRON OR TETRAEDRON
1829 DO j=1,8
1830 jj=brna(j,ii)
1831 vola(jj)=vola(jj)+one_over_8*pvolu(i)
1832 pa(jj)=pa(jj)+one_over_8*ppolh(i)*pvolu(i)
1833 rhoa(jj)=rhoa(jj)+one_over_8*rpolh(i)*pvolu(i)
1834 tka(jj)=tka(jj)+one_over_8*temp*pvolu(i)
1835 sspka(jj)=sspka(jj)+one_over_8*sqrt(gpolh(i)*ppolh(i)/rpolh(i))*pvolu(i)
1836 ua(1,jj)=ua(1,jj)+one_over_8*pu(1,i)*pvolu(i)
1837 ua(2,jj)=ua(2,jj)+one_over_8*pu(2,i)*pvolu(i)
1838 ua(3,jj)=ua(3,jj)+one_over_8*pu(3,i)*pvolu(i)
1839 ENDDO
1840 ENDIF
1841 ENDIF
1842 ENDDO
1843
1844 DO i=1,nna
1845 volno=vola(i)
1846 IF(volno > zero) THEN
1847 pa(i) =pa(i)/volno
1848 rhoa(i)=rhoa(i)/volno
1849 tka(i) =tka(i)/volno
1850 sspka(i) =sspka(i)/volno
1851 ua(1,i)=ua(1,i)/volno
1852 ua(2,i)=ua(2,i)/volno
1853 ua(3,i)=ua(3,i)/volno
1854 ENDIF
1855 ENDDO
1856 ENDIF ! PREPARE_ANIM
1857
1858
1859C---------------------------------------------------
1860C OPENING OF VENTHOLES
1861C---------------------------------------------------
1862 DO ivent=1,nvent
1863 area_vent(ivent)=zero
1864 pm_vent(ivent)=zero
1865 ENDDO
1866
1867C
1868 DO iel=1,nel
1869 IF (itagel(iel)<0) THEN
1870 ivent=-itagel(iel)
1871 area_vent(ivent)=area_vent(ivent)+elarea(iel)
1872 pm_vent(ivent)=pm_vent(ivent)+fel(iel)
1873 ENDIF
1874 ENDDO
1875
1876
1877
1878 DO ivent=1,nvent
1879 IF (area_vent(ivent)>zero) THEN
1880 pm_vent(ivent)=pm_vent(ivent)/area_vent(ivent)
1881 ELSE
1882 pm_vent(ivent)=zero
1883 ENDIF
1884C
1885 idef =ibaghol(1,ivent)
1886 idtpdef= ibaghol(11,ivent)
1887 idpdef = ibaghol(12,ivent)
1888C
1889 IF(idef==2) cycle
1890C
1891 pdef =rbaghol(1,ivent)
1892 dtpdefi=rbaghol(4,ivent)
1893 dtpdefc=rbaghol(5,ivent)
1894 tvent =rbaghol(3,ivent)
1895 tstope =rbaghol(14,ivent)
1896 IF (idtpdef==0) THEN
1897 IF(pm_vent(ivent)>pdef+pext) THEN
1898 dtpdefc = dtpdefc+dt1
1899 END IF
1900 ELSE IF (idtpdef==1) THEN
1901 IF (pm_vent(ivent)>pdef+pext) THEN
1902 idpdef = 1
1903 END IF
1904 IF (idpdef==1) THEN
1905 dtpdefc = dtpdefc+dt1
1906 END IF
1907 END IF
1908 ibaghol(12,ivent) = idpdef
1909 rbaghol(5,ivent) = dtpdefc
1910C
1911 ittf = ivolu(17)
1912 ttf =rvolu(60)
1913 DO k=1,20
1914 titrevent(k)=ibaghol(14+k,ivent)
1915 venttitle(k:k) = achar(titrevent(k))
1916 ENDDO
1917 IF(ittf==11.OR.ittf==12.OR.ittf==13) THEN
1918 IF (idef==0.AND.pm_vent(ivent)>pdef+pext
1919 . .AND.dtpdefc>dtpdefi
1920 . .AND.volph>em3*areap**three_half
1921 . .AND.tt<tstope+ttf) THEN
1922 idef=1
1923 WRITE(iout,'(A)')
1924 . ' ** AIRBAG VENT HOLE MEMBRANE IS DEFLATED **'
1925 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1926 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1927 WRITE(istdo,'(2A)')
1928 . ' ** VENT HOLE MEMBRANE IS DEFLATED ',venttitle
1929 ENDIF
1930 IF(idef==0 .AND. tt>tvent+ttf
1931 . .AND. tt<tstope+ttf) THEN
1932 idef=1
1933 WRITE(iout,'(A)') ' ** AIRBAG VENTING STARTS **'
1934 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1935 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1936 WRITE(istdo,'(2A)') ' ** VENTING STARTS ',venttitle
1937 ENDIF
1938 IF(idef==1 .AND. tt>=tstope+ttf) THEN
1939 idef=2
1940 WRITE(iout,'(A)') ' ** AIRBAG VENTING STOPS **'
1941 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1942 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1943 WRITE(istdo,'(2A)') ' ** VENTING STOPS ',venttitle
1944 ENDIF
1945C
1946 ELSE IF(ittf==0) THEN
1947 IF (idef==0.AND.pm_vent(ivent)>pdef+pext
1948 . .AND.dtpdefc>dtpdefi
1949 . .AND.volph>em3*areap**three_half
1950 . .AND.tt<tstope) THEN
1951 idef=1
1952 WRITE(iout,'(A)')
1953 . ' ** AIRBAG VENT HOLE MEMBRANE IS DEFLATED **'
1954 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1955 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1956 WRITE(istdo,'(2A)')
1957 . ' ** VENT HOLE MEMBRANE IS DEFLATED ',venttitle
1958 ENDIF
1959 IF(idef==0 .AND. tt>tvent
1960 . .AND. tt<tstope) THEN
1961 idef=1
1962 WRITE(iout,'(A)') ' ** AIRBAG VENTING STARTS **'
1963 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1964 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1965 WRITE(istdo,'(2A)') ' ** VENTING STARTS ',venttitle
1966 ENDIF
1967 IF(idef==1 .AND. tt>=tstope) THEN
1968 idef=2
1969 WRITE(iout,'(A)') ' ** AIRBAG VENTING STOPS **'
1970 WRITE(iout,'(3X,2(A,I10),2A)')' MONITORED VOLUME ',ivolu(1),
1971 . ' VENT HOLE NUMBER',ivent,' ',venttitle
1972 WRITE(istdo,'(2A)') ' ** VENTING STOPS ',venttitle
1973 ENDIF
1974 END IF
1975C
1976 ibaghol(1,ivent)=idef
1977 ENDDO
1978C---------------------------------------------------
1979C FORCES ON AIRBAGS AND OUTPUTS
1980C---------------------------------------------------
1981C FVSPMD(IFV)%AAA(1:3,1:NNT) = ZERO
1982C
1983 IF(ivolu(39)==1) THEN
1984 wfext0=wfext
1985 DO iel=1,nelt
1986 n1=elem(1,iel)
1987 n2=elem(2,iel)
1988 n3=elem(3,iel)
1989 nrx=norm(1,iel)
1990 nry=norm(2,iel)
1991 nrz=norm(3,iel)
1992 IF(iel<=nel) THEN
1993 ff=fel(iel)-pext*elarea(iel)
1994 ELSE
1995 ff=fel(iel)
1996 ENDIF
1997 fvspmd(ifv)%AAA(1,n1)=fvspmd(ifv)%AAA(1,n1)+third*ff*nrx
1998 fvspmd(ifv)%AAA(2,n1)=fvspmd(ifv)%AAA(2,n1)+third*ff*nry
1999 fvspmd(ifv)%AAA(3,n1)=fvspmd(ifv)%AAA(3,n1)+third*ff*nrz
2000 fvspmd(ifv)%AAA(1,n2)=fvspmd(ifv)%AAA(1,n2)+third*ff*nrx
2001 fvspmd(ifv)%AAA(2,n2)=fvspmd(ifv)%AAA(2,n2)+third*ff*nry
2002 fvspmd(ifv)%AAA(3,n2)=fvspmd(ifv)%AAA(3,n2)+third*ff*nrz
2003 fvspmd(ifv)%AAA(1,n3)=fvspmd(ifv)%AAA(1,n3)+third*ff*nrx
2004 fvspmd(ifv)%AAA(2,n3)=fvspmd(ifv)%AAA(2,n3)+third*ff*nry
2005 fvspmd(ifv)%AAA(3,n3)=fvspmd(ifv)%AAA(3,n3)+third*ff*nrz
2006C
2007 wfext=wfext+ff*vel(iel)*dt1
2008 ENDDO
2009 fsav(18)=fsav(18)+wfext-wfext0
2010 ENDIF
2011C----------------------
2012C OUTPUT TH GAUGES
2013C----------------------
2014 IF (nbgauge > 0) THEN
2015 DO i=1,nnt
2016 fvspmd(ifv)%GGG(1,i)=p(i)
2017 fvspmd(ifv)%GGG(2,i)=rho(i)
2018 fvspmd(ifv)%GGG(3,i)=tk(i)
2019 ENDDO
2020
2021 DO i=1,nna
2022 fvspmd(ifv)%GGA(1,i)=pa(i)
2023 fvspmd(ifv)%GGA(2,i)=rhoa(i)
2024 fvspmd(ifv)%GGA(3,i)=tka(i)
2025 ENDDO
2026 ENDIF
2027
2028 RETURN
2029 END
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alpha
#define alpha
Definition eval.h:35
fv_up_switch
subroutine fv_up_switch(nn, nel, elem, njet, npoly, lenh, nba, ibagjet, rbagjet, nvent, ibaghol, rbaghol, p, rho, tk, u, sspk, x, v, a, nsensor, sensor_tab, fsav, npc, tf, ivolu, rvolu, mpolh, qpolh, epolh, ppolh, rpolh, gpolh, npolh, ifvnod, rfvnod, ifvtri, ifvpoly, ifvtadr, ifvpolh, ifvpadr, info, nns, nntr, ifv, npolha, dlh, cpapolh, cpbpolh, cpcpolh, rmwpolh, itagel, elsini, icontact, idpolh, elfmass, elfvel, ibufa, elema, tagela, pa, rhoa, tka, ua, brna, nna, ntga, ibpolh, dtpolh, nnt, nelt, xxxa, vvva, ncona, porosity, ityp, igeo, sspka, geo, pm, ipm, tpolh, elfehpy, cpdpolh, cpepolh, cpfpolh, eltg, iparg, mattg, igrouptg, igroupc, elbuf_tab, cfl_coef, pdisp_old, pdisp, wfext, python)
Definition fv_up_switch.F:71
fvinjt6
subroutine fvinjt6(njet, ibagjet, rbagjet, npc, tf, nsensor, sensor_tab, scalt, datainj, python)
Definition fvinjt6.F:35
fvinjt8
subroutine fvinjt8(njet, ibagjet, rbagjet, npc, tf, nsensor, sensor_tab, scalt, igeo, geo, pm, ivolu, datainj, python)
Definition fvinjt8.F:36
fvinjt8_1
subroutine fvinjt8_1(njet, ibagjet, rbagjet, igeo, geo, pm, ivolu, rvolu, dmout, dhout)
Definition fvinjt8_1.F:32
fvtemp
subroutine fvtemp(ityp, efac_in, cpa_in, cpb_in, cpc_in, cpd_in, cpe_in, cpf_in, rmw_in, temp0_in, temp_in)
Definition fvtemp.F:35
fvvent0
subroutine fvvent0(elsout, aoutot, nvent, nelt, ittf, elarea, elsini, elem, itagel, svent, ibaghol, rvolu, rbaghol, poro, p, eltg, iparg, mattg, nel, porosity, ipm, pm, elbuf_tab, igroupc, igrouptg)
Definition fvvent0.F:45
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
anim_mod
Definition anim_mod.F:96
fvbag_mod
Definition fvbag_mod.F:44
fvbag_mod::fvspmd
type(fvbag_spmd), dimension(:), allocatable fvspmd
Definition fvbag_mod.F:129
fvbag_mod::fvdata
type(fvbag_data), dimension(:), allocatable fvdata
Definition fvbag_mod.F:128
fvmbag_meshcontrol_mod
Definition fvmbag_meshcontrol_mod.F:46
message_mod
Definition message_mod.F:1249
poro
subroutine poro(geo, nodpor, ms, x, v, w, af, am, skew, weight, nporgeo)
Definition poro.F:40
spmd_fvb_gath
subroutine spmd_fvb_gath(ifv, x, xxx, xxxa, xxxsa, ido)
Definition spmd_fvb_gath.F:40
spmd_fvb_igath
subroutine spmd_fvb_igath(ifv, itabs, itabg)
Definition spmd_fvb_igath.F:34
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87