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