OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ebcs10_nrf_mod Module Reference

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ ebcs10_nrf()

subroutine ebcs10_nrf_mod::ebcs10_nrf ( integer, intent(in) nseg,
integer, dimension(nseg), intent(in) iseg,
type(t_segvar), intent(inout) segvar,
dimension(3,numnod), intent(inout) a,
v,
w,
x,
integer, dimension(nod), intent(in) liste,
integer, intent(in) nod,
integer, dimension(4,nseg), intent(in) irect,
integer, dimension(nseg), intent(in) ielem,
integer, dimension(nseg), intent(in) iface,
la,
type(t_ebcs_nrf), intent(inout) ebcs,
integer, dimension(nparg,ngroup) iparg,
type(elbuf_struct_), dimension(ngroup), target elbuf_tab,
type(multi_fvm_struct), intent(in) multi_fvm,
integer, dimension(nixq,numelq), intent(in) ixq,
integer, dimension(nixs,numels), intent(in) ixs,
integer, dimension(nixtg,numeltg), intent(in) ixtg,
integer, dimension(4,nseg), intent(in) elem_adress,
intent(inout) fsky,
dimension(th_surf_num_channel,nsurf), intent(inout) fsavsurf,
type(output_), intent(inout) output,
intent(in) dt1,
intent(in) time )
Parameters
[in,out]outputoutput structure

Definition at line 46 of file ebcs10_nrf.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
iface
integer function iface(ip, n)
Definition iface.F:36
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
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