OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ebcs10_nrf.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!|| ebcs10_nrf_mod ../engine/source/boundary_conditions/ebcs/ebcs10_nrf.f
25!||--- called by ------------------------------------------------------
26!|| ebcs_main ../engine/source/boundary_conditions/ebcs/ebcs_main.f
27!||====================================================================
28 MODULE ebcs10_nrf_mod
29 IMPLICIT NONE
30 CONTAINS
31!||====================================================================
32!|| ebcs10_nrf ../engine/source/boundary_conditions/ebcs/ebcs10_nrf.F
33!||--- called by ------------------------------------------------------
34!|| ebcs_main ../engine/source/boundary_conditions/ebcs/ebcs_main.F
35!||--- uses -----------------------------------------------------
36!|| ale_mod ../common_source/modules/ale/ale_mod.F
37!|| ebcs_mod ../common_source/modules/boundary_conditions/ebcs_mod.F90
38!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
39!|| element_mod ../common_source/modules/elements/element_mod.F90
40!|| multi_fvm_mod ../common_source/modules/ale/multi_fvm_mod.F90
41!|| multimat_param_mod ../common_source/modules/multimat_param_mod.F90
42!|| output_mod ../common_source/modules/output/output_mod.F90
43!|| segvar_mod ../engine/share/modules/segvar_mod.F
44!|| th_surf_mod ../common_source/modules/interfaces/th_surf_mod.F
45!||====================================================================
46 SUBROUTINE ebcs10_nrf(NSEG,ISEG,SEGVAR,
47 . A,V,W,X,
48 . LISTE,NOD,IRECT,IELEM,IFACE,
49 . LA,EBCS,IPARG,ELBUF_TAB,MULTI_FVM,IXQ,IXS,IXTG,
50 . ELEM_ADRESS,FSKY,FSAVSURF,OUTPUT,DT1,TIME)
51C-----------------------------------------------
52C M o d u l e s
53C-----------------------------------------------
54 USE ebcs_mod , ONLY : t_ebcs_nrf
55 USE elbufdef_mod , ONLY : elbuf_struct_, l_bufel_, g_bufel_, buf_mat_
56 USE multi_fvm_mod , ONLY : multi_fvm_struct
57 USE segvar_mod , ONLY : t_segvar
58 USE ale_mod , ONLY : ale
59 USE th_surf_mod , ONLY : th_surf_num_channel
60 USE multimat_param_mod , ONLY : m51_n0phas, m51_nvphas
61 USE output_mod , ONLY : output_
62 use element_mod , only : nixs,nixq,nixtg
63C-----------------------------------------------
64C I m p l i c i t T y p e s
65C-----------------------------------------------
66#include "implicit_f.inc"
67C-----------------------------------------------
68C C o m m o n B l o c k s
69C-----------------------------------------------
70#include "param_c.inc"
71#include "com01_c.inc"
72#include "com04_c.inc"
73#include "parit_c.inc"
74C-----------------------------------------------
75C D u m m y A r g u m e n t s
76C-----------------------------------------------
77 my_real, INTENT(INOUT) :: fsavsurf(th_surf_num_channel,nsurf)
78 INTEGER,INTENT(IN) :: NSEG,NOD,ISEG(NSEG),LISTE(NOD),IRECT(4,NSEG),IELEM(NSEG),IFACE(NSEG)
79 INTEGER,INTENT(IN) :: IXQ(NIXQ,NUMELQ),IXS(NIXS,NUMELS),IXTG(NIXTG,NUMELTG)
80 my_real,INTENT(INOUT) :: A(3,NUMNOD)
81 my_real v(3,numnod),w(3,numnod),x(3,numnod),la(3,nod)
82 TYPE(t_ebcs_nrf), INTENT(INOUT) :: EBCS
83 INTEGER :: IPARG(NPARG,NGROUP)
84 TYPE(elbuf_struct_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
85 TYPE(multi_fvm_struct),INTENT(IN) :: MULTI_FVM
86 TYPE(t_segvar),INTENT(INOUT) :: SEGVAR
87 INTEGER, DIMENSION(4,NSEG), INTENT(IN) :: ELEM_ADRESS ! adress for fsky array (only used with parith/on)
88 my_real, DIMENSION(8,LSKY), INTENT(INOUT) :: fsky ! acceleration array for parith/on option
89 my_real,INTENT(IN) :: dt1 !< time step
90 my_real,INTENT(IN) :: time !< current time
91 TYPE(OUTPUT_), INTENT(INOUT) :: OUTPUT !< output structure
92C-----------------------------------------------
93C L o c a l V a r i a b l e s
94C-----------------------------------------------
95 TYPE(g_bufel_), POINTER :: GBUF
96 TYPE(L_BUFEL_) ,POINTER :: LBUF
97 INTEGER :: II,IS,KK,KSEG,NN(4),NNG(4),NUM,KTY,KLT,MFT,NGRP,ILOC,NPT,IVOI,IDX(6),IX(4)
98 INTEGER :: ICF_2d(2,4), ICF_3d(4,6), JJ, ISUBMAT, IPOS, NBMAT, MTN
99 my_real :: orient,rho,roc,fac1,fac2,vol,mass,
100 . x13,y13,z13,x24,y24,z24,
101 . alpha, beta, v0(3,nod),
102 . xn, yn, zn, tmp(3), vold, vnew, pold, pvois,
103 . dp0,mach,pp,ssp,rhoc2,
104 . tcar_p, tcar_vf, surf, eint,phase_alpha(21),phase_rho(21), phase_eint(21)
105 my_real :: de_in, de_out, dm_in, dm_out
106 my_real :: roou, enou, fluxo, fluxi
107 my_real :: ym
108 LOGICAL bFOUND
109 TYPE(buf_mat_) ,POINTER :: MBUF
110 INTEGER :: ADRESS ! adress for parit/on array
111 INTEGER :: SURF_ID ! EBCS Surface Identify
112C-----------------------------------------------
113 DATA icf_2d /1,2,2,3,3,4,4,1/
114 DATA icf_3d /1,4,3,2,3,4,8,7,5,6,7,8,1,2,6,5,2,3,7,6,1,5,8,4/
115C-----------------------------------------------
116C S o u r c e L i n e s
117C-----------------------------------------------
118 tcar_p = ebcs%TCAR_P
119 tcar_vf = ebcs%TCAR_VF
120 surf_id = ebcs%SURF_ID
121 de_in = zero
122 de_out = zero
123 dm_in = zero
124 dm_out = zero
125 ym = one
126
127 alpha = dt1/tcar_p
128 beta = dt1/max(dt1,tcar_vf)
129 IF(tcar_vf>=ep20)beta=zero
130 IF(dt1 == zero)THEN
131 alpha = one !zero
132 beta = one
133 ENDIF
134
135 IF(iale == 1)THEN
136
137 DO ii=1,nod
138 num=liste(ii)
139 v0(1,ii)=v(1,num)-w(1,num)
140 v0(2,ii)=v(2,num)-w(2,num)
141 v0(3,ii)=v(3,num)-w(3,num)
142 ENDDO
143 ELSEIF(ieuler == 1)THEN
144 DO ii=1,nod
145 num=liste(ii)
146 v0(1,ii)=v(1,num)
147 v0(2,ii)=v(2,num)
148 v0(3,ii)=v(3,num)
149 ENDDO
150 ENDIF
151
152 DO ii=1,nod
153 num=liste(ii)
154 la(1,ii)=zero
155 la(2,ii)=zero
156 la(3,ii)=zero
157 ENDDO
158
159 DO is=1,nseg
160
161 kseg=abs(iseg(is))
162 orient=float(iseg(is)/kseg)
163
164 !---OUTWARD NORMAL
165 IF(n2d == 0)THEN
166 jj = iface(is)
167 ix(1)=ixs(icf_3d(1,jj)+1,ielem(is))
168 ix(2)=ixs(icf_3d(2,jj)+1,ielem(is))
169 ix(3)=ixs(icf_3d(3,jj)+1,ielem(is))
170 ix(4)=ixs(icf_3d(4,jj)+1,ielem(is))
171 x13=x(1,ix(3))-x(1,ix(1))
172 y13=x(2,ix(3))-x(2,ix(1))
173 z13=x(3,ix(3))-x(3,ix(1))
174 x24=x(1,ix(4))-x(1,ix(2))
175 y24=x(2,ix(4))-x(2,ix(2))
176 z24=x(3,ix(4))-x(3,ix(2))
177 xn=y13*z24-z13*y24
178 yn=z13*x24-x13*z24
179 zn=x13*y24-y13*x24
180 fac2=one/sqrt(xn**2+yn**2+zn**2) * orient
181
182 xn = xn*fac2
183 yn = yn*fac2
184 zn = zn*fac2
185 surf = half/fac2
186 IF(ix(4) == ix(3))THEN ; npt=3;fac1=third; else; npt=4;fac1=fourth; ENDIF
187 ELSE
188 fac1=half
189 npt=2
190 jj = iface(is)
191 IF(numeltg > 0)THEN
192 ix(1) = ixtg(icf_2d(1,jj)+1,ielem(is))
193 ix(2) = ixtg(icf_2d(2,jj)+1,ielem(is))
194 ELSE
195 ix(1) = ixq(icf_2d(1,jj)+1,ielem(is))
196 ix(2) = ixq(icf_2d(2,jj)+1,ielem(is))
197 ENDIF
198 xn = zero
199 yn = -(-x(3,ix(2))+x(3,ix(1)))
200 zn = (-x(2,ix(2))+x(2,ix(1)))
201 fac2 = one/sqrt(yn*yn+zn*zn)*orient
202 IF(n2d==1)ym = half*(x(2,ix(1))+x(2,ix(2)))
203 yn=yn*fac2
204 zn=zn*fac2
205 surf = one/fac2
206 ENDIF
207
208 !-- VELOCITY
209 nn(1)=irect(1,is)
210 nn(2)=irect(2,is)
211 nn(3)=irect(3,is)
212 nn(4)=irect(4,is)
213 nng(1)=liste(nn(1))
214 nng(2)=liste(nn(2))
215 nng(3)=liste(nn(3))
216 nng(4)=liste(nn(4))
217
218 tmp(1:3) = zero
219 DO kk=1,npt
220 tmp(1) = tmp(1) + v0(1,nn(kk))
221 tmp(2) = tmp(2) + v0(2,nn(kk))
222 tmp(3) = tmp(3) + v0(3,nn(kk))
223 ENDDO
224 vold = ebcs%vold(is)
225 vnew = fac1 * (tmp(1)*xn + tmp(2)*yn + tmp(3)*zn)
226 IF(time == zero) THEN
227 IF(ale%GRID%NWALE == 7)THEN
228 vold = zero ! 0.0 with /ALE/GRID/FLOW-TRACKING
229 ENDIF
230 vold = vnew
231 ENDIF
232 !-- storage for next cycle
233 ebcs%vold(is) = vnew
234 !Static-PRessure (increment)
235 dp0 = ebcs%DP0(is)
236
237 !-- ADJACENT STATE
238 bfound = .false.
239 ivoi = ielem(is)
240 DO ngrp=1,ngroup
241 kty = iparg(5,ngrp)
242 klt = iparg(2,ngrp)
243 mft = iparg(3,ngrp)
244 IF(n2d == 0)THEN
245 IF(kty /= 1)cycle
246 ELSE
247 IF(kty /= 2 .AND. kty /= 7)cycle
248 ENDIF
249 IF (ivoi <= klt+mft)THEN
250 bfound = .true.
251 EXIT
252 ENDIF
253 ENDDO
254 IF(.NOT.bfound)cycle !next I
255 gbuf => elbuf_tab(ngrp)%GBUF
256 lbuf => elbuf_tab(ngrp)%BUFLY(1)%LBUF(1,1,1)
257 mtn = iparg(1,ngrp)
258
259 !PRESSURE
260 iloc = ivoi-mft-1
261 DO kk=1,6
262 idx(kk) = klt*(kk-1)
263 ENDDO
264 pvois = -third*(gbuf%SIG(idx(1)+iloc+1) + gbuf%SIG(idx(2)+iloc+1) + gbuf%SIG(idx(3)+iloc+1))
265 pold = ebcs%pold(is)
266 IF(time == zero)pold=pvois+dp0
267
268 !DENSITY
269 rho = gbuf%RHO(iloc+1)
270
271 !VOLUME
272 vol=gbuf%VOL(iloc+1)
273 mass = rho*vol
274
275 !SOUND SPEED
276 IF(multi_fvm%IS_USED)THEN
277 ssp = multi_fvm%SOUND_SPEED(ivoi)
278 ELSE
279 ssp = lbuf%SSP(iloc+1)
280 ENDIF
281
282 !ENERGY
283 eint = gbuf%EINT(iloc+1)
284
285 !VOL FRAC
286 IF(mtn == 51)THEN
287 mbuf => elbuf_tab(ngrp)%BUFLY(1)%MAT(1,1,1)
288 DO isubmat=1,4
289 ipos = 1
290 kk = (m51_n0phas + (isubmat-1)*m51_nvphas +ipos-1) * klt + iloc+1
291 phase_alpha(isubmat) = mbuf%VAR(kk)
292 ipos = 9
293 kk = (m51_n0phas + (isubmat-1)*m51_nvphas +ipos-1) * klt + iloc+1
294 phase_rho(isubmat) = mbuf%VAR(kk)
295 ipos = 8
296 kk = (m51_n0phas + (isubmat-1)*m51_nvphas +ipos-1) * klt + iloc+1
297 phase_eint(isubmat) = mbuf%VAR(kk)
298 enddo!next ITRIMAT
299 ELSEIF(mtn == 151)THEN
300 DO isubmat=1,multi_fvm%NBMAT
301 phase_alpha(isubmat) = multi_fvm%PHASE_ALPHA(isubmat,iloc+1+klt)
302 phase_rho(isubmat) = multi_fvm%PHASE_RHO(isubmat,iloc+1+klt)
303 phase_eint(isubmat) = multi_fvm%PHASE_EINT(isubmat,iloc+1+klt)
304 ENDDO
305 ENDIF
306 nbmat=ebcs%NBMAT
307
308 !-- ghost cell data for advection step ( inlet/outlet )
309 segvar%RHO(kseg)=rho
310 segvar%EINT(kseg)=eint
311 IF(segvar%has_phase_alpha)segvar%PHASE_ALPHA(1:4,kseg)=phase_alpha(1:4)
312 IF(segvar%has_phase_rho)segvar%PHASE_RHO(1:4,kseg)=phase_rho(1:4)
313 IF(segvar%has_phase_eint)segvar%PHASE_EINT(1:4,kseg)=phase_eint(1:4)
314
315 !-- post-treatment for mass and energy balance
316 roou = segvar%RHO(kseg)
317 enou = segvar%EINT(kseg)
318 fluxo=vnew*surf*dt1*ym
319 fluxi=min(fluxo,zero)
320 fluxo=max(fluxo,zero)
321 dm_out=dm_out-fluxo*rho
322 dm_in=dm_in-fluxi*rho
323 de_out=de_out-fluxo*eint
324 de_in=de_in-fluxi*eint
325 !-- FORMULATION
326 mach = one
327 IF(ssp /= zero)mach = abs(vnew / ssp)
328
329 IF(mach >= one .AND. vnew > zero)THEN
330 !outgoing supersonic velocity: State = Neighboring State
331 pp = pvois
332 ELSE
333 rhoc2 = rho*ssp*ssp
334 roc =sqrt(rho*rhoc2)
335C PP = Pold + (ROC*(Vnew-Vold)) !+ALPHA*(Pvois+DP0)/(ALPHA + ONE)
336 pp = one/(one+alpha)*(pold+roc*(vnew-vold))+alpha*(pvois+dp0)/(alpha + one)
337
338 ENDIF
339 ebcs%pold(is) = pp
340
341 IF(segvar%nbmat > 0)THEN
342 !Volume Fractions update
343 IF(vnew < zero)THEN !Returning flow, relaxation of the volume percentage
344 IF(mtn == 51)THEN
345 DO isubmat=1,4
346 ! (1-beta) vf_current + beta vf_0
347 phase_alpha(isubmat)=(one-beta)*phase_alpha(isubmat)
348 . + beta*mbuf%VAR((m51_n0phas+(isubmat-1)*m51_nvphas +23-1)* klt+iloc+1)
349 ENDDO
350 ELSEIF(mtn == 151)THEN
351
352 ENDIF
353 ENDIF
354 segvar%PHASE_ALPHA(1:4,kseg) = phase_alpha(1:4)
355 ENDIF
356
357 !expand pressure loading to segment nodes
358 DO kk=1,npt
359 la(1,nn(kk)) = la(1,nn(kk)) - pp*surf*xn*fac1
360 la(2,nn(kk)) = la(2,nn(kk)) - pp*surf*yn*fac1
361 la(3,nn(kk)) = la(3,nn(kk)) - pp*surf*zn*fac1
362 ENDDO
363
364 !/th/surf(massflow, velocity, pressure)
365 fsavsurf(2,surf_id) = fsavsurf(2,surf_id) + rho*surf*vnew !rho.S.u = dm/dt
366 fsavsurf(3,surf_id) = fsavsurf(3,surf_id) + surf*vnew !S.u
367 fsavsurf(4,surf_id) = fsavsurf(4,surf_id) + surf*pp !S.P
368 fsavsurf(6,surf_id) = fsavsurf(6,surf_id) + rho*surf*vnew*dt1 ! m<-m+dm (cumulative)
369
370 ! -------------
371 ! for parith/on option : need to update with FSKY array
372 IF(iparit>0) THEN
373 DO kk=1,npt
374 adress = elem_adress(kk,is) ! adress of FSKY array for element IS and node KK
375 fsky(1,adress) = -pp*surf*xn*fac1
376 fsky(2,adress) = -pp*surf*yn*fac1
377 fsky(3,adress) = -pp*surf*zn*fac1
378 fsky(4:8,adress) = zero
379 ENDDO
380 ENDIF
381 ! -------------
382 ENDDO
383
384 ! -------------
385 ! for parith/off option : update directly the acceleration array A()
386 IF(iparit==0) THEN
387 DO ii=1,nod
388 num=liste(ii)
389 a(1,num)=a(1,num)+la(1,ii)
390 a(2,num)=a(2,num)+la(2,ii)
391 a(3,num)=a(3,num)+la(3,ii)
392 ENDDO
393 ENDIF
394 ! -------------
395
396!$OMP CRITICAL
397 output%DATA%INOUT%DM_IN = output%DATA%INOUT%DM_IN + dm_in
398 output%DATA%INOUT%DM_OUT = output%DATA%INOUT%DM_OUT + dm_out
399 output%DATA%INOUT%DE_IN = output%DATA%INOUT%DE_IN + de_in
400 output%DATA%INOUT%DE_OUT = output%DATA%INOUT%DE_OUT + de_out
401!$OMP END CRITICAL
402
403 RETURN
404 END SUBROUTINE ebcs10_nrf
405 END MODULE ebcs10_nrf_mod
my_real
#define my_real
Definition cppsort.cpp:32
ebcs_main
subroutine ebcs_main(igrsurf, segvar, volmon, a, v, w, x, fsav, fv, ms, stifn, iparg, elbuf_tab, ebcs_tab, multi_fvm, ixq, ixs, ixtg, fsky, fsavsurf, time, dt1, nsensor, sensor_tab, python, npc, tf, snpc, stf, output, matparam, pm, bufmat)
Definition ebcs_main.F:68
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
ale_mod
Definition ale_mod.F:150
ale_mod::ale
type(ale_) ale
Definition ale_mod.F:253
ebcs10_nrf_mod
Definition ebcs10_nrf.F:28
ebcs10_nrf_mod::ebcs10_nrf
subroutine ebcs10_nrf(nseg, iseg, segvar, a, v, w, x, liste, nod, irect, ielem, iface, la, ebcs, iparg, elbuf_tab, multi_fvm, ixq, ixs, ixtg, elem_adress, fsky, fsavsurf, output, dt1, time)
Definition ebcs10_nrf.F:51
segvar_mod
Definition segvar_mod.F:63
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
segvar_mod::t_segvar
Definition segvar_mod.F:71
velocity
subroutine velocity(a, ar, v, vr, fzero, itab, nale)
Definition velocity.F:29