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