40
41
42
44
45
46
47
48
49
50#include "implicit_f.inc"
51
52
53
54#include "units_c.inc"
55#include "comlock.inc"
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,ILAY,IPT,NPT
77 INTEGER ,INTENT(INOUT) :: DMG_FLAG
78 INTEGER, DIMENSION(NEL) ,INTENT(IN) :: NGL,FWAVE_EL
79 INTEGER ,DIMENSION(NEL) ,INTENT(INOUT) :: OFFLY,FOFF
80
82 my_real,
DIMENSION(NUPARAM) ,
INTENT(IN) :: uparam
83 my_real,
DIMENSION(NEL) :: ssp,aldt,tdel1,tdel2,
84 . signxx,signyy,signxy
85 my_real ,
DIMENSION(NEL,2) ,
INTENT(OUT) :: crkdir
86 my_real ,
DIMENSION(NEL) ,
INTENT(INOUT) :: off,dfmax,dadv
87 my_real ,
DIMENSION(NEL,NUVAR),
TARGET,
INTENT(INOUT) :: uvar
88
89
90
91 INTEGER I,J,NINDXF1,NINDXF2,NINDXD1,NINDXD2,IDEB,IMOD,ISRATE,
92 . ISIDE,ITGLASS
93 INTEGER ,DIMENSION(NEL) :: INDXF1,INDXF2,INDXD1,INDXD2,RFLAG
94
95 my_real sigmax,aa,bb,cc,cr,cosx,sinx,s1,s2,s3,k_ic,k_th,
96 . geored,cr_foil,cr_air,cr_core,cr_edge,
alpha,alphai,beta,
97 . kres1,kres2,exp_n,exp_m,v0,vc,reflen,dmg1,dmg2,dmg3,tanphi
98 my_real ,
DIMENSION(NEL) :: tpropg,ai,formf,sigp_akt,dam1,dam2,
99 . sigp1,sigp2,sigdt1,sigdt2,sig_dtf1
100 . sp1,sp2,sp3,sin2,cos2,cosi,tfact1,tfact2,cr_len,cr_depth,cr_ang
101 my_real,
DIMENSION(:),
POINTER :: dir1,dir2
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127 rflag(1:nel) = 0
128 nindxf1 = 0
129 nindxf2 = 0
130 nindxd1 = 0
131 nindxd2 = 0
132
133 exp_n = uparam(1)
134 cr_foil = uparam(2)
135 cr_air = uparam(3)
136 cr_core = uparam(4)
137 cr_edge = uparam(5)
138 k_ic = uparam(6)
139 k_th = uparam(7)
140 v0 = uparam(8)
141 vc = uparam(9)
143 geored = uparam(11)
144 reflen = uparam(14)
145 imod = nint(uparam(15))
146 israte = nint(uparam(16))
147 ideb = nint(uparam(17))
148 iside = nint(uparam(18))
149 trelax = uparam(19)
150 kres1 = uparam(20)
151 kres2 = uparam(21)
152 itglass = nint(uparam(22))
153
154 IF (trelax > 0) dmg_flag = 1
156 exp_m = one / (one + exp_n)
157
158 sigp_min(1:nel) = uvar(1:nel,5)
159 sigp_max(1:nel) = uvar(
160 formf(1:nel) = uvar(1:nel,7)
161 ai(1:nel) = uvar(1:nel,8)
162 sigp_akt(1:nel) = uvar(1:nel,9)
163 dam1(1:nel) = uvar(1:nel,11)
164 dam2(1:nel) = uvar(1:nel,12)
165
166
167
168 IF (itglass == 1 .and. (ipt == 1 .or. ipt =THEN
170 . nel ,nuparam ,nuvar ,time ,timestep ,
171 . uparam ,ngl ,signxx ,signyy ,signxy ,
172 . uvar ,off ,ipt ,nindxf1 ,indxf1 ,
173 . tdel1 )
174
175 ELSE IF (time == zero) THEN
176
177
178 uvar(1:nel,15) = one
179
180 END IF
181
182 DO i = 1,nel
183 IF (fwave_el(i) > 0) THEN
184 dadv(i) =
min(one, geored / sqrt(aldt(i)/reflen) )
185 rflag(i) = -1
186 uvar(i,15) = one
187 ELSE IF (dadv(i) < one) THEN
188 rflag(i) = -1
189 ELSE IF (dadv(i) == one) THEN
190 rflag(i) = 1
191 ENDIF
192 tpropg(i) = aldt(i) /
min(vc,ssp(i))
193 tpropg(i) =
max(tpropg(i), timestep)
194 END DO
195
196
197
198 DO i = 1,nel
199 s1 = signxx(i)
200 s2 = signyy(i)
201 cc = signxy(i)
202 aa = (s1 + s2)*half
203 bb = (s1 - s2)*half
204 cr = sqrt(bb**2 + cc**2)
205 sigp1(i) = aa + cr
206 sigp2(i) = aa - cr
207 sp1(i) = sigp1(i)
208 sp2(i) = sigp2(i)
209 sp3(i) = zero
210 END DO
211
212
213
214 IF (israte == 0) THEN
215
216 DO i = 1,nel
217 IF (off(i) == one) THEN
218 sigdt1(i) = abs(sigp1(i) - uvar(i,1)) /
max(timestep, em20)
219 sigdt2(i) = abs(sigp2(i) - uvar(i,2)) /
max(timestep, em20)
220 sigdt1(i) =
alpha * sigdt1(i) + alphai * uvar(i,3)
221 sigdt2(i) =
alpha * sigdt2(i) + alphai * uvar(i,4)
222 ENDIF
223 END DO
224 ELSE
225
226 DO i = 1,nel
227 IF (off(i) == one) THEN
228 uvar(i,21) = abs(sigp1(i) - uvar(i,1)) /
max(timestep, em20)
231
232 uvar(i,71) = abs(sigp2(i) - uvar(i,2)) /
max(timestep, em20)
235 ELSE
236 sigdt1(i) = zero
237 sigdt2(i) = zero
238 ENDIF
239 END DO
242 ENDIF
243
244
245
246 IF (ipt == npt) THEN
247 DO i = 1,nel
248 sig_dtf1(i) = sigp_akt(i) *abs(sigdt1(i))**exp_m
249 sig_dtf2(i) = sigp_akt(i) *abs(sigdt2(i))**exp_m
250 sig_dtf1(i) =
max(sig_dtf1(i),sigp_min(i))
251 sig_dtf1(i) =
min(sig_dtf1(i),sigp_max(i))
252 sig_dtf2(i) =
max(sig_dtf2(i),sigp_min(i))
253 sig_dtf2(i) =
min(sig_dtf2(i),sigp_max(i))
254 END DO
255 ELSEIF (ipt < npt) THEN
256 DO i = 1,nel
257 IF (uvar(i,10) == one) THEN
258 sig_dtf1(i) = sigp_akt(i)*abs(sigdt1(i))**exp_m
259 sig_dtf2(i) = sigp_akt(i)*abs(sigdt2(i))**exp_m
260 sig_dtf1(i) =
max(sig_dtf1(i),sigp_min(i))
261 sig_dtf1(i) =
min(sig_dtf1(i),sigp_max(i))
262 sig_dtf2(i) =
max(sig_dtf2(i),sigp_min(i))
263 sig_dtf2(i) =
min(sig_dtf2(i),sigp_max(i))
264 ELSE IF (iside == 1) THEN
265 sig_dtf1(i) = sigp_akt(i)*abs(sigdt1(i))**exp_m
266 sig_dtf2(i) = sigp_akt(i)*abs(sigdt2(i))**exp_m
267 sig_dtf1(i) =
max(sig_dtf1(i),sigp_min(i))
268 sig_dtf1(i) =
min(sig_dtf1(i),sigp_max(i))
269 sig_dtf2(i) =
max(sig_dtf2(i),sigp_min(i))
270 sig_dtf2(i) =
min(sig_dtf2(i),sigp_max(i))
271 ELSE
272 sig_dtf1(i) = sigp_max(i)
273 sig_dtf2(i) = sigp_max(i)
274 ENDIF
275 END DO
276 ENDIF
277
278 DO i = 1,nel
279 sig_dtf1(i) = sig_dtf1(i)*dadv(i)
280 sig_dtf2(i) = sig_dtf2(i)*dadv(i)
281 END DO
282
283
284
285 DO i = 1,nel
286 IF (off(i) == one) THEN
287 IF (tdel1(i) == zero) THEN
288 IF (sigp1(i) > sig_dtf1(i) .and. uvar(i,15) == one) THEN
289 tdel1(i) = time
290 nindxf1 = nindxf1+1
291 indxf1(nindxf1) = i
292 ENDIF
293 ENDIF
294
295 IF (tdel1(i) > zero) THEN
296 tfact1(i) = (time - tdel1(i)) / tpropg(i)
297 tfact1(i) =
min(one, tfact1(i))
298 IF (uvar(i,11) > zero .and. tfact1(i) == one) THEN
299 offly(i) = -1
300 ENDIF
301 dfmax(i) =
max(dfmax(i), tfact1(i))
302 nindxd1 = nindxd1+1
303 indxd1(nindxd1) = i
304 ENDIF
305
306 ENDIF
307 ENDDO
308
309
310
311 IF (nindxf1 > 0) THEN
312#include "vectorize.inc"
313 DO j=1,nindxf1
314 i = indxf1(j)
315 s1 = signxy(i)
316 s2 = sigp1(i) - signxx(i)
317 cr = sqrt(s1**2 + s2**2)
318 IF (cr > zero) THEN
319 crkdir(i,1) = s1 / cr
320 crkdir(i,2) = s2 / cr
321 ELSEIF (s1 > s2) THEN
322 crkdir(i,1) = zero
323 crkdir(i,2) = one
324 ELSE
325 crkdir(i,1) = one
326 crkdir(i,2) = zero
327 ENDIF
328 ENDDO
329 ENDIF
330
331
332
333 IF (nindxd1 > 0) THEN
334#include "vectorize.inc"
335 DO j=1,nindxd1
336 i = indxd1(j)
337 s1 = signxx(i)
338 s2 = signyy(i)
339 s3 = signxy(i)
340 cosx = crkdir(i,1)
341 sinx = crkdir(i,2)
342 cos2(i) = cosx * cosx
343 sin2(i) = sinx * sinx
344 cosi(i) = cosx * sinx
345
346 sp1(i) = cos2(i)*s1 + sin2(i)*s2 + two*cosi(i)*s3
347 sp2(i) = sin2(i)*s1 + cos2(i)*s2 - two*cosi(i)*s3
348 sp3(i) = cosi(i)*(s2 - s1) + (cos2(i) - sin2(i))*s3
349
350 beta =
max(zero , one - tfact1(i))
351 dmg1 =
max(kres1, one - (one - kres1) * tfact1(i))
352 dmg3 =
min(one, 0.6 + 0.4 * beta)
353 dam1(i) = dmg1
354
355 sp1(i) = sp1(i) * dmg1
356 sp3(i) = sp3(i) * dmg3
357
358 signxx(i) = cos2(i)*sp1(i) + sin2(i)*sp2(i) - two*cosi(i)*sp3(i)
359 signyy(i) = sin2(i)*sp1(i) + cos2(i)*sp2(i) + two*cosi(i)*sp3(i)
360 signxy(i) = cosi(i)*(sp1(i) - sp2(i)) + (cos2(i) - sin2(i))*sp3(i)
361 ENDDO
362 ENDIF
363
364
365
366 DO i = 1,nel
367 IF (off(i) == one) THEN
368 IF (tdel1(i) > zero .and. tdel2(i) == zero) THEN
369 IF (sp2(i) > sig_dtf2(i)) THEN
370 tdel2(i) = time
371 nindxf2 = nindxf2 + 1
372 indxf2(nindxf2) = i
373 ENDIF
374 ENDIF
375
376 IF (tdel2(i) > zero) THEN
377 tfact2(i) = (time - tdel2(i)) / tpropg(i)
378 tfact2(i) =
min(one, tfact2(i))
379 IF (uvar(i,12) > zero .and. tfact2(i) == one) THEN
380 IF (offly(i) == 0) THEN
381 offly(i) = -2
382 ELSE IF (offly(i) == -1) THEN
383 offly(i) = -3
384 ENDIF
385 ENDIF
386 nindxd2 = nindxd2+1
387 indxd2(nindxd2) = i
388 ENDIF
389 ENDIF
390 ENDDO
391
392 IF (nindxd2 > 0) THEN
393#include "vectorize.inc"
394 DO j=1,nindxd2
395 i = indxd2(j)
396 beta =
max(zero , one - tfact2(i))
397 dmg2 =
max(kres2, one - (one - kres2) * tfact2(i))
398 dmg3 =
min(one, 0.6 + 0.4 * beta)
399 dam2(i) = dmg2
400
401 sp2(i) = sp2(i) * dmg2
402 sp3(i) = sp3(i) * dmg3
403
404 signxx(i) = cos2(i)*sp1(i) + sin2(i)*sp2(i) - two*cosi(i)*sp3(i)
405 signyy(i) = sin2(i)*sp1(i) + cos2(i)*sp2(i) + two*cosi(i)*sp3(i)
406 signxy(i) = cosi(i)*(sp1(i) - sp2(i)) + (cos2(i) - sin2(i))*sp3(i)
407 IF (dmg2 == zero) THEN
408 foff(i) = 0
409 ENDIF
410 ENDDO
411 ENDIF
412
413 DO i=1,nel
414 uvar(i,1) = sigp1(i)
415 uvar(i,2) = sigp2(i)
416 uvar(i,3) = sigdt1(i)
417 uvar(i,4) = sigdt2(i)
418 uvar(i,11) = dam1(i)
419 uvar(i,12) = dam2(i)
420 END DO
421
422 IF (nindxf1 > 0) THEN
423 DO j=1,nindxf1
424 i = indxf1(j)
425#include "lockon.inc"
426 IF (rflag(i) == 1) THEN
427
428 WRITE(iout, 3000) ngl(i),ilay,ipt
429 WRITE(istdo,3100) ngl(i),ilay,ipt,time
430 ELSEIF (rflag(i) == -1) THEN
431
432 WRITE(iout, 4000) ngl(i),ilay,ipt
433 WRITE(istdo,4100) ngl(i),ilay,ipt,time
434 ENDIF
435 IF (ideb == 1 .and. abs(rflag(i)) == 1) THEN
436 tanphi = atan(crkdir(i,2) / sign(
max(abs(crkdir(i,1)),em20),crkdir(i,1)))
437 WRITE(iout,3500) ngl(i),tanphi,dadv(i)
438 WRITE(iout,3700) ngl(i),sigp1(i),sigdt1(i),sig_dtf1(i)
439 ENDIF
440#include "lockoff.inc"
441 ENDDO
442 ENDIF
443
444 IF (nindxf2 > 0) THEN
445 DO j=1,nindxf2
446 i = indxf2(j)
447#include "lockon.inc"
448 WRITE(iout, 5000) ngl(i),ilay,ipt
449 WRITE(istdo,5100) ngl(i),ilay,ipt,time
450 IF (ideb == 1 .and. abs(rflag(i)) == 1) THEN
451 WRITE(iout,3600) ngl(i),dadv(i)
452 WRITE(iout,3800) ngl(i),sigp2(i),sigdt2(i),sig_dtf2(i)
453 ENDIF
454#include "lockoff.inc"
455 ENDDO
456 ENDIF
457
458 3000 FORMAT(1x, 'FAILURE INITIALIZATION IN SHELL',i10,1x,' 1ST DIR, LAYER',i2,1x,'INT POINT',i2)
459 3100 FORMAT(1x, 'FAILURE INITIALIZATION IN SHELL',i10,1x,' 1ST DIR, LAYER',i2,1x,'INT POINT',i2,
460 . 1x, 'AT TIME ',1pe12.4)
461 3500 FORMAT(10x,'SHELL ',i10,' MAJ ANGLE= ',1pe12.4,' STRESS REDUCTION= ',1pe12.4)
462 3600 FORMAT(10x,'SHELL ',i10,' STRESS REDUCTION= ',1pe12.4)
463 3700 FORMAT(10x,'SHELL ',i10,' SIGP1=',1pe12.4,' SIGDT1=',1pe12.4,' SIG_DTF1=',1pe12.4)
464 3800 FORMAT(10x,'SHELL ',i10,' SIGP2=',1pe12.4,' SIGDT2=',1pe12.4,' SIG_DTF2=',1pe12.4)
465 4000 FORMAT(1x, 'FAILURE ADVANCEMENT IN SHELL',i10,1x,' 1ST DIR, LAYER',i2,1x,'INT POINT',i2)
466 4100 FORMAT(1x, 'FAILURE ADVANCEMENT IN SHELL',i10,1x,' 1ST DIR, LAYER',i2,1x,'INT POINT',i2,
467 . 1x, 'AT TIME ',1pe12.4)
468 5000 FORMAT(1x, 'FAILURE IN SHELL',i10,1x,' 2ND DIR, LAYER',i2,1x,'INT POINT',i2)
469 5100 FORMAT(1x, 'FAILURE IN SHELL',i10,1x,' 2ND DIR, LAYER',i2,1x,'INT POINT',i2,
470 . 1x, 'AT TIME ',1pe12.4)
471
472 RETURN
subroutine fail_brokmann(nel, nuparam, nuvar, time, timestep, uparam, ngl, signxx, signyy, signxy, uvar, off, ipt, nindxf, indxf, tdel)
integer, dimension(:), allocatable index_over_50_cycles