OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pblast_3.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!|| pblast_3 ../engine/source/loads/pblast/pblast_3.F
25!||--- called by ------------------------------------------------------
26!|| pblast_load_computation ../engine/source/loads/pblast/pblast.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!|| my_barrier ../engine/source/system/machine.f
30!|| pblast_parameters__air_burst ../common_source/modules/loads/pblast_mod.F90
31!||--- uses -----------------------------------------------------
32!|| groupdef_mod ../common_source/modules/groupdef_mod.F
33!|| h3d_inc_mod ../engine/share/modules/h3d_inc_mod.F
34!|| h3d_mod ../engine/share/modules/h3d_mod.F
35!|| names_and_titles_mod ../common_source/modules/names_and_titles_mod.F
36!|| output_mod ../common_source/modules/output/output_mod.F90
37!|| pblast_mod ../common_source/modules/loads/pblast_mod.F90
38!|| th_mod ../engine/share/modules/th_mod.f
39!|| th_surf_mod ../common_source/modules/interfaces/th_surf_mod.F
40!||====================================================================
41 SUBROUTINE pblast_3(PBLAST,
42 * ILOADP ,FAC ,A ,V ,X ,
43 1 IADC ,FSKY ,LLOADP ,FEXT ,NODA_SURF,NODA_PEXT,
44 2 ITAB ,H3D_DATA ,NL ,DTMIN_LOC ,WFEXT_LOC,
45 4 TH_SURF ,NSEGPL)
46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49 USE h3d_mod
50 USE pblast_mod
51 USE groupdef_mod
52 USE h3d_inc_mod
53 USE th_surf_mod , ONLY : th_surf_
54 USE names_and_titles_mod , ONLY : ncharline
55 USE th_mod ,ONLY : th_has_noda_pext
56 USE output_mod , ONLY : h3d_has_noda_pext, anim_has_noda_pext
57C-----------------------------------------------
58C I m p l i c i t T y p e s
59C-----------------------------------------------
60#include "implicit_f.inc"
61#include "comlock.inc"
62#include "param_c.inc"
63C-----------------------------------------------
64C C o m m o n B l o c k s
65C-----------------------------------------------
66#include "com04_c.inc"
67#include "com08_c.inc"
68#include "parit_c.inc"
69#include "scr14_c.inc"
70#include "scr16_c.inc"
71#include "mvsiz_p.inc"
72#include "units_c.inc"
73#include "sysunit.inc"
74#include "tabsiz_c.inc"
75C-----------------------------------------------
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 INTEGER,INTENT(IN) :: LLOADP(SLLOADP)
79 INTEGER,INTENT(INOUT) :: ILOADP(SIZLOADP,NLOADP)
80 INTEGER,INTENT(IN) :: IADC(*)
81 INTEGER, INTENT(IN) :: ITAB(NUMNOD),NL
82 my_real,INTENT(INOUT) :: dtmin_loc
83 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT_LOC
84 my_real,INTENT(IN) :: v(3,numnod),x(3,numnod)
85 my_real,INTENT(INOUT) :: fac(lfacload,nloadp)
86 my_real,INTENT(INOUT) :: a(3,numnod),fsky(8,sfsky/8)
87 my_real,INTENT(INOUT) :: fext(3,numnod)
88 my_real,INTENT(INOUT) :: noda_surf(numnod)
89 my_real,INTENT(INOUT) :: noda_pext(numnod)
90 TYPE(h3d_database),INTENT(IN) :: H3D_DATA
91 TYPE (TH_SURF_) , INTENT(INOUT) :: TH_SURF
92 INTEGER, INTENT(INOUT) :: NSEGPL
93 TYPE(pblast_),INTENT(INOUT) :: PBLAST
94C-----------------------------------------------
95C L o c a l V a r i a b l e s
96C-----------------------------------------------
97 CHARACTER(LEN=NCHARLINE) :: MSGOUT1,MSGOUT2
98 TYPE(FRIEDLANDER_PARAMS_) :: FRIEDLANDER_PARAMS
99 INTEGER N1, N2, N3, N4,IL,IS,IAD,I,IANIM_OR_H3D,IZ_UPDATE,ABAC_ID,ISIZ_SEG,IERR1
100 INTEGER ID, ITA_SHIFT,IMODEL
101 INTEGER :: NS,KSURF
102 INTEGER :: curve_id1, curve_id2, NN(4), NDT
103 LOGICAL :: BOOL_UNDERGROUND_CURRENT_SEG, IS_DECAY_TO_BE_COMPUTED
104 LOGICAL,SAVE :: IS_UPDATED
105 INTEGER :: NWARN !< number of segment for which it is not possible to correlate the positive impulse. It may happen that for a given Pmax and dt0 value even building a triangle shape leads to a lower impulse that the targeted one.
106
107 my_real :: zx,zy,zz,norm,nx,ny,nz ! target face centroid and normal vector
108 my_real nnx,nny,nnz,norm_nn, norm2_nn, tmp_var ! ground vector
109 my_real hz ! height of centroid
110 my_real lg,zg ! ground distance and scaled ground distance(SHDC : scaled horizontal distance from charge)
111 my_real zc ! scaled height of charge, and its value in ft/(lb**1/3)caled height of triple point
112 my_real angle_g ! angle of incidence at groundh front incident pressure and refelcted pressurese
113 my_real projz(3),projdet(3), tmp(3)
114 my_real z_struct ! scaled height for target face
115 my_real alpha_zc ! interpolation between curves Scaled Height of Charge
116
117 my_real cos_theta, alpha_inci, alpha_refl,p_inci,p_refl,z,
118 . i_inci,i_refl,dt_0,t_a,wave_refl,wave_inci, w13
119 my_real dnorm,xdet,ydet,zdet,tdet,wtnt,pmin,t_stop,dx,dy,dz,p,
120 . fac_m_bb, fac_l_bb, fac_t_bb, fac_p_bb, fac_i_bb, ta_first, tt_star
121 my_real hc ! height of charge
122 my_real z1_
123 my_real decay_inci,decay_refl,area
124 my_real npt
125 my_real :: ta_inf_loc, zeta
126 my_real :: surf_patch
127
128 my_real pi_
129 DATA pi_/3.141592653589793238462643d0/
130
131 my_real dzc
132 DATA dzc /0.07058823500000000/
133
134 my_real :: cst_255_div_ln_z1_on_zn
135 DATA cst_255_div_ln_z1_on_zn/-38.147316611455952998/
136
137 my_real ::log10_
138 DATA log10_ /2.30258509299405000000/
139
140 my_real :: cst_180
141 DATA cst_180 /180.000000000000000000/
142
143 INTEGER,EXTERNAL :: DICHOTOMIC_SEARCH_R_DESC,DICHOTOMIC_SEARCH_R_ASC
144
145 my_real fac_unit ! convert scaled distance from cm/g^1/3 to ft/lb^1/3
146 DATA fac_unit/3.966977216838196139019/
147
148C-----------------------------------------------
149C D e s c r i p t i o n
150C-----------------------------------------------
151C This subroutines is applying pressure load to a segment submitted to a blast wave (AIR BURST model).
152C Preussre time histories are built from "UFC 3-340-02, Dec. 5th 2008" tables which are hardcoded in unit system {g, cm, mus}
153C (table sampling : 256 points ; see pblast_mod.F)
154C
155C T* = T + TA_INF (shift with first arrival time for all pblast option, TA_INF=0. if ita_shift/=2)
156C If request made to update blast profile (iz_update==2) then it is made once on T*=TDET
157C
158C-----------------------------------------------
159C P r e - C o n d i t i o n
160C-----------------------------------------------
161 !--- subroutine not relevant if no /LOAD/PBLAST option
162 IF(pblast%NLOADP_B==0)RETURN
163
164 !--- time step
165 ta_first = fac(07,nl) ! arrival time for first segment to be loaded
166 il= nl-nloadp_f
167 tt_star = tt + pblast%PBLAST_DT%TA_INF
168 iz_update = iloadp(06,nl)
169 tdet = fac(01,nl)
170 ta_first = fac(07,nl) + pblast%PBLAST_DT%TA_INF ! arrival time for first segment to be loaded
171 IF(iz_update ==1)THEN
172 !if no update during engine
173 dtmin_loc = (one+em06)*(ta_first - tt)
174 dtmin_loc=max(pblast%PBLAST_TAB(il)%DTMIN, dtmin_loc) !go directly to arrival time but avoid too
175 IF(tt_star<ta_first)RETURN ! (nothing to load for now)
176 ELSE
177 IF(tdet > tt)THEN
178 dtmin_loc = (one+em06)*(tdet - tt)
179 dtmin_loc=max(pblast%PBLAST_TAB(il)%DTMIN, dtmin_loc)!go directly to arrival time but avoid too
180 ELSE
181 dtmin_loc = pblast%PBLAST_TAB(il)%DTMIN
182 ENDIF
183 IF(tt_star<tdet)RETURN ! (nothing to update for now)
184 ENDIF
185C-----------------------------------------------,
186C S o u r c e C o d e
187C-----------------------------------------------
188 wfext_loc = zero
189 ianim_or_h3d = anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT + anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT !output
190
191 !Z Range (tables)
192 z1_ = 0.500000000000000
193 !ZN = 400.000
194
195 !translation from Working unit System to {big bang} unit system
196 fac_m_bb = fac_mass*ep03
197 fac_l_bb = fac_length*ep02
198 fac_t_bb = fac_time*ep06
199 fac_p_bb = fac_m_bb/fac_l_bb/fac_t_bb/fac_t_bb
200 fac_i_bb = fac_p_bb*fac_t_bb !/FAC_M_bb**THIRD
201
202 is_updated=.false.
203 CALL my_barrier
204
205 !-----------------------------------------------,
206 ! AIR BURST
207 !-----------------------------------------------
208 il = nl-nloadp_f
209 tdet = fac(01,nl)
210 xdet = fac(02,nl)
211 ydet = fac(03,nl)
212 zdet = fac(04,nl)
213 wtnt = fac(05,nl)
214 pmin = fac(06,nl)
215 ta_first = fac(07,nl)
216 nnx = fac(08,nl)
217 nny = fac(09,nl)
218 nnz = fac(10,nl)
219 hc = fac(11,nl)
220 alpha_zc = fac(12,nl) !curve_id1+alpha_zc
221 t_stop = fac(13,nl)
222 isiz_seg = iloadp(01,nl)/4
223 is = iloadp(02,nl)
224 iz_update = iloadp(06,nl)
225 abac_id = iloadp(07,nl)
226 id = iloadp(08,nl) !user_id
227 ita_shift = iloadp(09,nl)
228 ndt = iloadp(10,nl)
229 imodel = iloadp(11,nl)
230 ta_inf_loc = ep20
231
232 is_decay_to_be_computed = .false.
233 IF(imodel == 2)is_decay_to_be_computed=.true.
234 nwarn = 0
235
236 !curve_id1+alpha_zc
237 curve_id1=int(alpha_zc)
238 curve_id2=min(10,curve_id1+1)
239 alpha_zc = alpha_zc - curve_id1
240
241 ierr1 = 0
242 w13 = (wtnt*fac_m_bb)**third ! '*FAC_M' g->work unit '/FAC_M' : WORK_UNIT -> g
243 z = zero
244 norm2_nn=nnx*nnx+nny*nny+nnz*nnz
245 norm_nn =sqrt(norm2_nn)
246
247 !---------------------------------------------
248 ! LOOP ON SEGMENTS (4N or 3N)
249 !---------------------------------------------
250!$OMP DO SCHEDULE(GUIDED,MVSIZ)
251 DO i = 1,isiz_seg
252
253 bool_underground_current_seg = .false.
254
255 n1=lloadp(iloadp(4,nl)+4*(i-1))
256 n2=lloadp(iloadp(4,nl)+4*(i-1)+1)
257 n3=lloadp(iloadp(4,nl)+4*(i-1)+2)
258 n4=lloadp(iloadp(4,nl)+4*(i-1)+3)
259
260 IF(n4 == 0 .OR. n3 == n4)THEN
261 !3-NODE-SEGMENT
262 pblast%PBLAST_TAB(il)%NPt(i) = three
263 npt = three
264 !Segment centroid
265 zx = x(1,n1)+x(1,n2)+x(1,n3)
266 zy = x(2,n1)+x(2,n2)+x(2,n3)
267 zz = x(3,n1)+x(3,n2)+x(3,n3)
268 zx = zx*third
269 zy = zy*third
270 zz = zz*third
271 !Normal vector (NX,NY,NZ) = 2*S*n where |n|=1.0
272 nx = (x(2,n3)-x(2,n1))*(x(3,n3)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n3)-x(2,n2))
273 ny = (x(3,n3)-x(3,n1))*(x(1,n3)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n3)-x(3,n2))
274 nz = (x(1,n3)-x(1,n1))*(x(2,n3)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n3)-x(1,n2))
275 !NORM = 2*S
276 norm = sqrt(nx*nx+ny*ny+nz*nz)
277 ELSE
278 !4-NODE-SEGMENT
279 pblast%PBLAST_TAB(il)%NPt(i) = four
280 npt = four
281 !Segment centroid
282 zx = x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4)
283 zy = x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4)
284 zz = x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4)
285 zx = zx*fourth
286 zy = zy*fourth
287 zz = zz*fourth
288 !Normal vector (NX,NY,NZ) = 2*S*n where |n|=1.0
289 nx = (x(2,n3)-x(2,n1))*(x(3,n4)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n4)-x(2,n2))
290 ny = (x(3,n3)-x(3,n1))*(x(1,n4)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n4)-x(3,n2))
291 nz = (x(1,n3)-x(1,n1))*(x(2,n4)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n4)-x(1,n2))
292 !NORM = 2*S
293 norm = sqrt(nx*nx+ny*ny+nz*nz)
294 ENDIF
295 nn(1)=n1
296 nn(2)=n2
297 nn(3)=n3
298 nn(4)=n4
299 pblast%PBLAST_TAB(il)%N(1,i) = n1
300 pblast%PBLAST_TAB(il)%N(2,i) = n2
301 pblast%PBLAST_TAB(il)%N(3,i) = n3
302 pblast%PBLAST_TAB(il)%N(4,i) = n4
303
304!------------------------------- ---
305
306 !--------------------------------------------!
307 ! Update Wave parameters !
308 ! (experimental) !
309 ! If requested. Otherwise use Starter param. !
310 ! Default : do not update:use Starter param. !
311 !--------------------------------------------!
312 IF(iz_update == 2)THEN
313
314 dtmin_loc = ep20
315
316 !Dist
317 dx = (xdet - zx)*fac_l_bb ! => working unit to cm
318 dy = (ydet - zy)*fac_l_bb ! => ... to cm
319 dz = (zdet - zz)*fac_l_bb ! => ... to cm
320 dnorm = sqrt(dx*dx+dy*dy+dz*dz)
321
322 !scaled distance
323 z = dnorm / w13 !in abac unit ID g,cm,mus
324
325 !Determine Height of centroid (structure face)
326 ! Basis = DETPOINT - NN
327 projdet(1)=xdet + nnx
328 projdet(2)=ydet + nny
329 projdet(3)=zdet + nnz
330 !Height is length Proj->Det. Storing Det->Proj into NN array
331 hz=-(nnx*zx + nnx*zy + nnz*zz - nnx*projdet(1)-nnx*projdet(2)-nnz*projdet(3))/hc
332
333 cos_theta = zero
334
335 IF(hz < zero)THEN
336
337 !underground segment (nothing to compute)
338 p_inci = zero
339 i_inci = zero
340 p_refl = zero
341 i_refl = zero
342 dt_0 = ep20
343 t_a = ep20
344 decay_refl = one
345 decay_inci = one
346 bool_underground_current_seg = .true.
347
348 ELSE
349
350 z_struct = hz*fac_l_bb/w13 ! scaled value in same unit system as hardcoded tables {g,cm,mus,Mbar}
351
352 !Scaled Height of Charge
353 zc = hc * fac_l_bb/w13! scaled value in same unit system as hardcoded tables {g,cm,mus,Mbar}
354 zc = zc/fac_unit ! ft/lb**0.3333 (fig 2-13 : 10 plots, one for a given ZC value)
355
356 !Horizontal Distance between Charge and Centroid : LG
357 ! ZG = scaled distance |ProjC->ProjZ|
358 projz(1) = zx + hz*nnx/hc
359 projz(2) = zy + hz*nny/hc
360 projz(3) = zz + hz*nnz/hc
361 tmp(1) = (projz(1)-projdet(1))
362 tmp(2) = (projz(2)-projdet(2))
363 tmp(3) = (projz(3)-projdet(3))
364 lg = sqrt(tmp(1)*tmp(1)+tmp(2)*tmp(2)+tmp(3)*tmp(3))
365 zg = lg*fac_l_bb/w13 !scaled horizontal distance (ground) ; same unit system as hardcoded tables {g,cm,mus,Mbar}
366
367 !Angle of structural face (mach wave is planar wave)
368 cos_theta = (projdet(1)-projz(1))*nx + (projdet(2)-projz(2))*ny + (projdet(3)-projz(3))*nz
369 cos_theta = cos_theta/max(em20,lg*norm)
370
371 !determine angle of incidence at ground (AANGLE_g) and interpolation factor (alpha_angle)
372 tmp(1)=xdet-projz(1)
373 tmp(2)=ydet-projz(2)
374 tmp(3)=zdet-projz(3)
375 tmp_var=sqrt( tmp(1)*tmp(1) + tmp(2)*tmp(2) + tmp(3)*tmp(3) )
376 angle_g = -( nnx*tmp(1) + nny*tmp(2) + nnz*tmp(3) ) /hc/tmp_var !must be between [-PI_,PI_]
377 angle_g = min(one,max(-one,angle_g)) ! bound it to expected range (epsilon)
378 angle_g = acos(angle_g)
379 angle_g = cst_180/pi_*angle_g !debug purpose
380 IF(angle_g < zero .AND. pblast%PBLAST_TAB(il)%TAGMSG(i) == 0 )THEN
381 write(iout,*) 'Warning : /LOAD/PBLAST id=',id,' NEGATIVE ANGLE,',angle_g
382 write(istdo,*) 'Warning : /LOAD/PBLAST id=',id,' NEGATIVE ANGLE,',angle_g
383 if(n4 == 0 .or. n3 == n4)then
384 write(iout,fmt= '(A,3I11)')' FACE:',itab((/n1,n2,n3/)),' SEEMS TO BE BELOW THE GROUND'
385 write(istdo,fmt='(A,3I11)')' FACE:',itab((/n1,n2,n3/)),' SEEMS TO BE BELOW THE GROUND'
386 else
387 write(iout,fmt= '(A,4I11)')' FACE:',itab((/n1,n2,n3,n4/)),' SEEMS TO BE BELOW THE GROUND'
388 write(istdo,fmt='(A,4I11)')' FACE:',itab((/n1,n2,n3,n4/)),' SEEMS TO BE BELOW THE GROUND'
389 endif
390 angle_g = zero
391 pblast%PBLAST_TAB(il)%TAGMSG(i) = 1
392 ELSEIF(angle_g > 85.00000 .AND. pblast%PBLAST_TAB(il)%TAGMSG(i) == 0)THEN
393 write(iout, fmt='(A,I0,A,E10.4)')'Warning : /LOAD/PBLAST id=',id,' ANGLE IS OVER THE UPPER BOUND,',angle_g
394 write(istdo,fmt='(A,I0,A,E10.4)')'Warning : /LOAD/PBLAST id=',id,' ANGLE IS OVER THE UPPER BOUND,',angle_g
395 if(n4 == 0 .or. n3 == n4)then
396 write(iout, fmt='(A,3I11)') ' ANGLE SET TO 85.00 FOR FACE:',itab((/n1,n2,n3/))
397 write(istdo,fmt='(A,3I11)') ' angle set to 85.00 for face:',itab((/n1,n2,n3/))
398 else
399 write(iout, FMT='(a,4i11)') ' angle set to 85.00 for face:',itab((/n1,n2,n3,n4/))
400 write(istdo,FMT='(a,4i11)') ' angle set to 85.00 for face:',itab((/n1,n2,n3,n4/))
401 endif
402
403 ANGLE_g = 85.00000
404 PBLAST%PBLAST_TAB(IL)%TAGMSG(I) = 1
405 ENDIF
406
407
408 !------------------------------------------------------------------!
409 CALL PBLAST_PARAMETERS__AIR_BURST( PBLAST,
410 + Z_struct, Zc, Zg, ANGLE_g, W13, TDET,
411 + FAC_P_bb, FAC_I_bb, FAC_T_bb,
412 + IS_DECAY_TO_BE_COMPUTED,
413 + ID,'load',.TRUE.,
414 + FRIEDLANDER_PARAMS,NWARN)
415 P_inci = FRIEDLANDER_PARAMS%P_inci
416 P_refl = FRIEDLANDER_PARAMS%P_refl
417 I_inci = FRIEDLANDER_PARAMS%I_inci
418 I_refl = FRIEDLANDER_PARAMS%I_refl
419 T_A = FRIEDLANDER_PARAMS%T_A
420 DT_0 = FRIEDLANDER_PARAMS%DT_0
421 decay_inci = FRIEDLANDER_PARAMS%decay_inci
422 decay_refl = FRIEDLANDER_PARAMS%decay_refl
423 !------------------------------------------------------------------!
424
425 TA_INF_LOC = MIN(TA_INF_LOC, T_A)
426
427 !update wave parameters
428 PBLAST%PBLAST_TAB(IL)%cos_theta(I) = cos_theta
429 PBLAST%PBLAST_TAB(IL)%P_inci(I) = P_inci
430 PBLAST%PBLAST_TAB(IL)%P_refl(I) = P_refl
431 PBLAST%PBLAST_TAB(IL)%ta(I) = T_A
432 PBLAST%PBLAST_TAB(IL)%t0(I) = DT_0
433 PBLAST%PBLAST_TAB(IL)%decay_inci(I) = decay_inci
434 PBLAST%PBLAST_TAB(IL)%decay_refl(I) = decay_refl
435
436 ENDIF !HZ
437
438 DTMIN_LOC = MIN(DTMIN_LOC,DT_0/NDT)
439 IZ_UPDATE = 1 !update done
440 ILOADP(06,NL) = IZ_UPDATE
441 IS_UPDATED=.TRUE.
442
443 ELSE
444
445 !Use wave parameters from Starter
446 cos_theta = PBLAST%PBLAST_TAB(IL)%cos_theta(I)
447 P_inci = PBLAST%PBLAST_TAB(IL)%P_inci(I)
448 P_refl = PBLAST%PBLAST_TAB(IL)%P_refl(I)
449 T_A = PBLAST%PBLAST_TAB(IL)%ta(I)
450 DT_0 = PBLAST%PBLAST_TAB(IL)%t0(I)
451 decay_inci = PBLAST%PBLAST_TAB(IL)%decay_inci(I)
452 decay_refl = PBLAST%PBLAST_TAB(IL)%decay_refl(I)
453 DTMIN_LOC = PBLAST%PBLAST_TAB(IL)%DTMIN
454
455 ENDIF ! IZ_UPDATE
456
457!------------------------------- ---
458
459 !Coefficients for wave superimposition
460 !PressureLoad = Reflected_Pressure * cos2X + IncidentPressure * (1 + cos2X -2 cosX)
461 IF(cos_theta<=ZERO)THEN
462 !Surface not facing the point of explosion
463 alpha_refl = ZERO
464 alpha_inci = ONE
465 ELSE
466 alpha_refl = cos_theta**2 ! cos**2 a
467 alpha_inci = ONE + cos_theta - TWO * alpha_refl ! 1 + cos a -2 cos**2 a
468 ENDIF
469
470 !Building pressure waves from Friedlander model. (Modified model can bu introduced later if needed)
471 WAVE_INCI = ZERO
472 WAVE_REFL = ZERO
473 IF(TT_STAR >= T_A)THEN
474 WAVE_INCI = P_inci*(ONE-(TT_STAR-T_A)/DT_0)*exp(-DECAY_inci*(TT_STAR-T_A)/DT_0)
475 WAVE_REFL = P_refl*(ONE-(TT_STAR-T_A)/DT_0)*exp(-DECAY_refl*(TT_STAR-T_A)/DT_0)
476 ELSE
477 WAVE_INCI = ZERO
478 WAVE_REFL = ZERO
479 ENDIF
480 P = alpha_refl * WAVE_REFL + alpha_inci * WAVE_INCI
481 P = MAX(P,PMIN)
482 PBLAST%PBLAST_TAB(IL)%PRES(I) = P
483
484 !Expand Pressure load to nodes
485 SURF_PATCH = HALF*SQRT(NX*NX+NY*NY+NZ*NZ) / NPT
486 PBLAST%PBLAST_TAB(IL)%FX(I)= -P * HALF*NX / NPT
487 PBLAST%PBLAST_TAB(IL)%FY(I)= -P * HALF*NY / NPT
488 PBLAST%PBLAST_TAB(IL)%FZ(I)= -P * HALF*NZ / NPT
489 PBLAST%PBLAST_TAB(IL)%SURF_PATCH(I) = SURF_PATCH
490
491 !External Force work
492 ! on a given node : DW = <F,V>*dt
493 ! for this current 4-node or 3-node face : DW = sum( <F_k,V_k>*dt k=1,NPT) where F_k=Fel/NPT
494 WFEXT_LOC=WFEXT_LOC+DT1*( PBLAST%PBLAST_TAB(IL)%FX(I) * SUM( V( 1, NN(1:NINT(NPt)) ) )
495 + + PBLAST%PBLAST_TAB(IL)%FY(I) * SUM( V( 2, NN(1:NINT(NPt)) ) )
496 + + PBLAST%PBLAST_TAB(IL)%FZ(I) * SUM( V( 3, NN(1:NINT(NPt)) ) )
497 + )
498
499C----- /TH/SURF -------
500 IF(TH_SURF%LOADP_FLAG > 0 ) THEN
501 NSEGPL = NSEGPL + 1
502 AREA = SURF_PATCH*NPT
503 DO NS=TH_SURF%LOADP_KSEGS(NSEGPL) +1,TH_SURF%LOADP_KSEGS(NSEGPL+1)
504 KSURF = TH_SURF%LOADP_SEGS(NS)
505 th_surf%channels(4,KSURF)= th_surf%channels(4,KSURF) + AREA*P ! mean pressure
506 th_surf%channels(5,KSURF)= th_surf%channels(5,KSURF) + AREA ! surface where pressure is applied
507 ENDDO
508 ENDIF
509
510
511 ENDDO!next I
512!$OMP END DO
513
514.AND. IF(IMODEL == 2 NWARN > 0)THEN
515 MSGOUT1=''
516 WRITE(MSGOUT1,FMT='(i0,a)') NWARN,
517 . ' segment(s) has excessive positive impulse regarding the peak pressure and positive duration.'
518 MSGOUT2=''
519 MSGOUT2='a triangular waveform will be used instead to maximize the impulse. defining a pmin VALUE is strongly recommended'
520 write(IOUT , FMT='(a,i10,/a,/a)') "Updated parameters for /LOAD/PBLAST id=", ID, MSGOUT1, MSGOUT2
521 write(ISTDO, FMT='(a,i10,/a,/a)') "Updated parameters for /LOAD/PBLAST id=", ID, MSGOUT1, MSGOUT2
522 ENDIF
523
524 CALL MY_BARRIER
525 IF(IS_UPDATED)THEN
526#include "lockon.inc"
527 !---arrival time
528 ZETA = FAC(07,NL)
529 FAC(07,NL) = MIN(TA_INF_LOC, FAC(07,NL)) !smp min value
530 !---time step
531 DTMIN_LOC = (ONE+EM06)*(FAC(07,NL) - TT) ! go directly to trigger time
532 DTMIN_LOC=MAX(PBLAST%PBLAST_TAB(IL)%DTMIN, DTMIN_LOC)
533 !---no update on next cycle
534 IZ_UPDATE = 1 !update done
535 ILOADP(06,NL) = IZ_UPDATE
536#include "lockoff.inc"
537!$OMP SINGLE
538 write(IOUT ,FMT='(a,i10,a,e16.8,a,e16.8)') "Updated parameters for /LOAD/PBLAST id=",
539 . ID,' previous first arrival time :',ZETA,
540 . ' is now updated to :',FAC(07,NL)
541 write(istdo,FMT='(a,i10,a,e16.8,a,e16.8)') "Updated parameters for /LOAD/PBLAST id=",
542 . ID,' previous first arrival time :',ZETA,
543 . ' is now updated to :',FAC(07,NL)
544!$OMP END SINGLE
545 ENDIF
546
547 !-------------------------------------------------------------------!
548 ! FORCE ASSEMBLY !
549 ! /PARITH/OFF : F directly added in A(1:3,1:NUMNOD). !
550 ! /PARITH/ON : F added FSKY & and automatically treated later !
551 !-------------------------------------------------------------------!
552 ! SPMD/SMP Parith/OFF
553 IF(IPARIT==0) THEN
554!$OMP SINGLE
555 DO I = 1,ISIZ_SEG
556 N1=LLOADP(ILOADP(4,NL)+4*(I-1))
557 N2=LLOADP(ILOADP(4,NL)+4*(I-1)+1)
558 N3=LLOADP(ILOADP(4,NL)+4*(I-1)+2)
559 N4=LLOADP(ILOADP(4,NL)+4*(I-1)+3)
560 A(1,N1)=A(1,N1)+PBLAST%PBLAST_TAB(IL)%FX(I)
561 A(2,N1)=A(2,N1)+PBLAST%PBLAST_TAB(IL)%FY(I)
562 A(3,N1)=A(3,N1)+PBLAST%PBLAST_TAB(IL)%FZ(I)
563 A(1,N2)=A(1,N2)+PBLAST%PBLAST_TAB(IL)%FX(I)
564 A(2,N2)=A(2,N2)+PBLAST%PBLAST_TAB(IL)%FY(I)
565 A(3,N2)=A(3,N2)+PBLAST%PBLAST_TAB(IL)%FZ(I)
566 A(1,N3)=A(1,N3)+PBLAST%PBLAST_TAB(IL)%FX(I)
567 A(2,N3)=A(2,N3)+PBLAST%PBLAST_TAB(IL)%FY(I)
568 A(3,N3)=A(3,N3)+PBLAST%PBLAST_TAB(IL)%FZ(I)
569 IF(PBLAST%PBLAST_TAB(IL)%NPt(I) == FOUR)THEN
570 A(1,N4)=A(1,N4)+PBLAST%PBLAST_TAB(IL)%FX(I)
571 A(2,N4)=A(2,N4)+PBLAST%PBLAST_TAB(IL)%FY(I)
572 A(3,N4)=A(3,N4)+PBLAST%PBLAST_TAB(IL)%FZ(I)
573 ENDIF
574 ENDDO
575!$OMP END SINGLE
576 ELSE
577!$OMP DO SCHEDULE(GUIDED,MVSIZ)
578 DO I = 1,ISIZ_SEG
579 IAD =IADC(ILOADP(4,NL)+4*(I-1))
580 FSKY(1,IAD) =PBLAST%PBLAST_TAB(IL)%FX(I)
581 FSKY(2,IAD) =PBLAST%PBLAST_TAB(IL)%FY(I)
582 FSKY(3,IAD) =PBLAST%PBLAST_TAB(IL)%FZ(I)
583 IAD =IADC(ILOADP(4,NL)+4*(I-1)+1)
584 FSKY(1,IAD) =PBLAST%PBLAST_TAB(IL)%FX(I)
585 FSKY(2,IAD) =PBLAST%PBLAST_TAB(IL)%FY(I)
586 FSKY(3,IAD) =PBLAST%PBLAST_TAB(IL)%FZ(I)
587 IAD =IADC(ILOADP(4,NL)+4*(I-1)+2)
588 FSKY(1,IAD) =PBLAST%PBLAST_TAB(IL)%FX(I)
589 FSKY(2,IAD) =PBLAST%PBLAST_TAB(IL)%FY(I)
590 FSKY(3,IAD) =PBLAST%PBLAST_TAB(IL)%FZ(I)
591 IF(PBLAST%PBLAST_TAB(IL)%NPt(I) == FOUR)THEN
592 IAD =IADC(ILOADP(4,NL)+4*(I-1)+3)
593 FSKY(1,IAD) =PBLAST%PBLAST_TAB(IL)%FX(I)
594 FSKY(2,IAD) =PBLAST%PBLAST_TAB(IL)%FY(I)
595 FSKY(3,IAD) =PBLAST%PBLAST_TAB(IL)%FZ(I)
596 ENDIF
597 ENDDO
598!$OMP END DO
599 ENDIF !IPARIT
600
601
602 !-------------------------------------------!
603 ! ANIMATION FILE /ANIM/VECT/FEXT !
604 ! H3D FILE /H3D/NODA/FEXT !
605 !-------------------------------------------!
606
607!$OMP SINGLE
608 IF(IANIM_OR_H3D > 0) THEN
609 DO I = 1,ISIZ_SEG
610 N1=PBLAST%PBLAST_TAB(IL)%N(1,I)
611 N2=PBLAST%PBLAST_TAB(IL)%N(2,I)
612 N3=PBLAST%PBLAST_TAB(IL)%N(3,I)
613 N4=PBLAST%PBLAST_TAB(IL)%N(4,I)
614 FEXT(1,N1) = FEXT(1,N1)+PBLAST%PBLAST_TAB(IL)%FX(I)
615 FEXT(2,N1) = FEXT(2,N1)+PBLAST%PBLAST_TAB(IL)%FY(I)
616 FEXT(3,N1) = FEXT(3,N1)+PBLAST%PBLAST_TAB(IL)%FZ(I)
617 FEXT(1,N2) = FEXT(1,N2)+PBLAST%PBLAST_TAB(IL)%FX(I)
618 FEXT(2,N2) = FEXT(2,N2)+PBLAST%PBLAST_TAB(IL)%FY(I)
619 FEXT(3,N2) = FEXT(3,N2)+PBLAST%PBLAST_TAB(IL)%FZ(I)
620 FEXT(1,N3) = FEXT(1,N3)+PBLAST%PBLAST_TAB(IL)%FX(I)
621 FEXT(2,N3) = FEXT(2,N3)+PBLAST%PBLAST_TAB(IL)%FY(I)
622 FEXT(3,N3) = FEXT(3,N3)+PBLAST%PBLAST_TAB(IL)%FZ(I)
623 IF(PBLAST%PBLAST_TAB(IL)%NPt(I)==FOUR)THEN
624 FEXT(1,N4) = FEXT(1,N4)+PBLAST%PBLAST_TAB(IL)%FX(I)
625 FEXT(2,N4) = FEXT(2,N4)+PBLAST%PBLAST_TAB(IL)%FY(I)
626 FEXT(3,N4) = FEXT(3,N4)+PBLAST%PBLAST_TAB(IL)%FZ(I)
627 ENDIF
628 ENDDO
629 ENDIF
630.OR..OR. IF(TH_HAS_NODA_PEXT > 0 ANIM_HAS_NODA_PEXT > 0 H3D_HAS_NODA_PEXT > 0) THEN
631 DO I = 1,ISIZ_SEG
632 N1 = PBLAST%PBLAST_TAB(IL)%N(1,I)
633 N2 = PBLAST%PBLAST_TAB(IL)%N(2,I)
634 N3 = PBLAST%PBLAST_TAB(IL)%N(3,I)
635 N4 = PBLAST%PBLAST_TAB(IL)%N(4,I)
636 SURF_PATCH = PBLAST%PBLAST_TAB(IL)%SURF_PATCH(I)
637 NODA_SURF(N1) = NODA_SURF(N1) + SURF_PATCH
638 NODA_SURF(N2) = NODA_SURF(N2) + SURF_PATCH
639 NODA_SURF(N3) = NODA_SURF(N3) + SURF_PATCH
640 P = PBLAST%PBLAST_TAB(IL)%PRES(I) * SURF_PATCH
641 NODA_PEXT(N1) = NODA_PEXT(N1) + P
642 NODA_PEXT(N2) = NODA_PEXT(N2) + P
643 NODA_PEXT(N3) = NODA_PEXT(N3) + P
644 IF(PBLAST%PBLAST_TAB(IL)%NPT(I) == FOUR)THEN
645 NODA_SURF(N4) = NODA_SURF(N4) + SURF_PATCH
646 NODA_PEXT(N4) = NODA_PEXT(N4) + P
647 ENDIF
648 ENDDO
649 ENDIF
650!$OMP END SINGLE
651
652 RETURN
653
654C-----------------------------------------------
655 IF (IERR1/=0) THEN
656 write(iout,*)' ** error in memory allocation - pblast loading'
657 write(istdo,*)' ** error in memory allocation - pblast loading'
658 CALL ARRET(2)
659 END IF
660C-----------------------------------------------
661
662 END SUBROUTINE
my_real
#define my_real
Definition cppsort.cpp:32
the
end diagonal values have been computed in the(sparse) matrix id.SOL
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
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
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
groupdef_mod
Definition groupdef_mod.F:662
h3d_inc_mod
Definition h3d_inc_mod.F:42
h3d_mod
Definition h3d_mod.F:308
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::ncharline
integer, parameter ncharline
Definition names_and_titles_mod.F:1001
th_mod
Definition th_mod.F:42
th_mod::th_has_noda_pext
integer th_has_noda_pext
Definition th_mod.F:121
th_surf_mod
OPTION /TH/SURF outputs of Pressure and Area needed Tabs.
Definition th_surf_mod.F:60
pblast_3
subroutine pblast_3(pblast, iloadp, fac, a, v, x, iadc, fsky, lloadp, fext, noda_surf, noda_pext, itab, h3d_data, nl, dtmin_loc, wfext_loc, th_surf, nsegpl)
Definition pblast_3.F:46
my_barrier
subroutine my_barrier
Definition machine.F:31
h3d_mod::h3d_database
Definition h3d_mod.F:435
th_surf_mod::th_surf_
Definition th_surf_mod.F:75