36
37
38
39#include "implicit_f.inc"
40
41
42
43#include "units_c.inc"
44#include "comlock.inc"
45
46
47
48 INTEGER, INTENT(IN) :: NEL,NUPARAM,NUVAR,IPG,NPG,ILAY
49 INTEGER, DIMENSION(NEL), INTENT(IN) :: NGL
51 my_real,
DIMENSION(NEL),
INTENT(IN) :: depsxx,depsyy,
52 . depszz,depsxy,depsyz,depszx,signxx,signyy,signzz,
53 . signxy,signyz,signzx,aldt
54 my_real,
DIMENSION(NUPARAM),
INTENT(IN) :: uparam
55
56
57
58 INTEGER, DIMENSION(NEL), INTENT(INOUT) :: NOFF
59 my_real,
DIMENSION(NEL),
INTENT(INOUT) :: dfmax,tdel,off,loff
60 my_real,
DIMENSION(NEL,6),
INTENT(INOUT) :: dmg_scale
61 my_real,
DIMENSION(NEL,NUVAR),
INTENT(INOUT) :: uvar
62
63
64
65 INTEGER I,J,NINDX,NINDX2,FAILMOD,MODE,ISHAP11T,ISHAP11C,
66 . ISHAP22T,ISHAP22C,ISHAP33T,ISHAP33C,ISHAP12T,ISHAP12C,
67 . ISHAP23T,ISHAP23C,ISHAP31T,ISHAP31C,FAILIP,NMOD
68 INTEGER ,DIMENSION(NEL) :: INDX, INDX2
69 INTEGER ,DIMENSION(NEL,12) :: FMODE
70 my_real dam(nel,12),ener(nel,12),le(nel)
71 my_real sigma_11t,sigma_11c,sigma_22t,sigma_22c,sigma_33t,
72 . sigma_33c,sigma_12t,sigma_12c,sigma_23t,sigma_23c,
73 . sigma_31t,sigma_31c,g_11t,g_11c,g_22t,g_22c,g_33t,g_33c,
74 . g_12t,g_12c,g_23t,g_23c,g_31t,g_31c
75
76
77
78
79
80 failip = nint(uparam(2))
81 failip =
min(failip,npg)
82 sigma_11t = uparam(3)
83 sigma_11c = uparam(4)
84 sigma_22t = uparam(5)
85 sigma_22c = uparam(6)
86 sigma_33t = uparam(7)
87 sigma_33c = uparam(8)
88 sigma_12t = uparam(9)
89 sigma_12c = uparam(10)
90 sigma_23t = uparam(11)
91 sigma_23c = uparam(12)
92 sigma_31t = uparam(13)
93 sigma_31c = uparam(14)
94 g_11t = uparam(15)
95 g_11c = uparam(16)
96 g_22t = uparam(17)
97 g_22c = uparam(18)
98 g_33t = uparam(19)
99 g_33c = uparam(20)
100 g_12t = uparam(21)
101 g_12c = uparam(22)
102 g_23t = uparam(23)
103 g_23c = uparam(24)
104 g_31t = uparam(25)
105 g_31c = uparam(26)
106 ishap11t = nint(uparam(27))
107 ishap11c = nint(uparam(28))
108 ishap22t = nint(uparam(29))
109 ishap22c = nint(uparam(30))
110 ishap33t = nint(uparam(31))
111 ishap33c = nint(uparam(32))
112 ishap12t = nint(uparam(33))
113 ishap12c = nint(uparam(34))
114 ishap23t = nint(uparam(35))
115 ishap23c = nint(uparam(36))
116 ishap31t = nint(uparam(37))
117 ishap31c = nint(uparam(38))
118 nmod = nint(uparam(39))
119
120
121 IF (uvar(1,25) == zero) THEN
122 uvar(1:nel,25) = aldt(1:nel)
123 ENDIF
124 le(1:nel) = uvar(1:nel,25)
125
126
127 DO j = 1,12
128 DO i = 1,nel
129
130 dam(i,j) = uvar(i,j)
131
132 ener(i,j) = uvar(i,j+12)
133 ENDDO
134 ENDDO
135
136
137
138
139
140 nindx = 0
141 nindx2 = 0
142 indx(1:nel) = 0
143 indx2(1:nel) = 0
144
145
146 DO i = 1, nel
147
148 IF ((loff(i) == one).AND.(off(i) == one)) THEN
149
150
151 fmode(i,1:12) = 0
152 failmod = 0
153
154
155 mode = 1
156 IF (signxx(i) > sigma_11t .AND. signxx(i) > zero .AND. dam(i,mode) < one) THEN
157
158 IF (ishap11t == 1) THEN
159 dam(i,mode) = dam(i,mode) + le(i)*depsxx(i)*sigma_11t/(two*g_11t)
160 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
161
162 ELSEIF (ishap11t == 2) THEN
163 ener(i,mode) = ener(i,mode) + signxx(i)*le(i)*depsxx(i)
164 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
165 dam(i,mode) = (one - exp(-ener(i,mode)/g_11t))
166 IF (dam(i,mode) > 0.999d0) dam(i,mode)
167 ENDIF
168 dam(i,mode) =
min(dam(i,mode),one)
169
170 dfmax(i) =
max(dfmax(i),dam(i,mode))
171 ENDIF
172
173 IF (dam(i,mode) >= one) THEN
174 failmod = failmod + 1
175 fmode(i,mode) = 1
176 ENDIF
177
178
179 mode = 2
180 IF (signyy(i) > sigma_22t .AND. signyy(i) > zero .AND. dam(i,mode) < one) THEN
181
182 IF (ishap22t == 1) THEN
183 dam(i,mode) = dam(i,mode) + le(i)*depsyy(i)*sigma_22t/(two*g_22t)
184 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
185
186 ELSEIF (ishap22t == 2) THEN
187 ener(i,mode) = ener(i,mode) + signyy(i)*le(i)*depsyy(i)
188 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
189 dam(i,mode) = (one - exp(-ener(i,mode)/g_22t))
190 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
191 ENDIF
192 dam(i,mode) =
min(dam(i,mode),one)
193
194 dfmax(i) =
max(dfmax(i),dam(i,mode))
195 ENDIF
196
197 IF (dam(i,mode) >= one) THEN
198 failmod = failmod + 1
199 fmode(i,mode) = 1
200 ENDIF
201
202
203 mode = 3
204 IF (signxy(i) > sigma_12t .AND. signxy(i) > zero .AND. dam(i,mode) < one) THEN
205
206 IF (ishap12t == 1) THEN
207 dam(i,mode) = dam(i,mode) + le(i)*depsxy(i)*sigma_12t/(four*g_12t)
208 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
209
210 ELSEIF (ishap12t == 2) THEN
211 ener(i,mode) = ener(i,mode) + signxy(i)*le(i)*half*depsxy(i)
212 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
213 dam(i,mode) = (one - exp(-ener(i,mode)/g_12t))
214 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
215 ENDIF
216 dam(i,mode) =
min(dam(i,mode),one)
217
218 dfmax(i) =
max(dfmax(i),dam(i,mode))
219 ENDIF
220
221 IF (dam(i,mode) >= one) THEN
222 failmod = failmod + 1
223 fmode(i,mode) = 1
224 ENDIF
225
226
227 mode = 4
228 IF (-signxx(i) > sigma_11c .AND. signxx(i) < zero .AND. dam(i,mode) < one) THEN
229
230 IF (ishap11c == 1) THEN
231 dam(i,mode) = dam(i,mode) + le(i)*abs(depsxx(i))*sigma_11c/(two*g_11c)
232 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
233
234 ELSEIF (ishap11c == 2) THEN
235 ener(i,mode) = ener(i,mode) + abs(signxx(i))*le(i)*abs(depsxx(i))
236 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
237 dam(i,mode) = (one - exp(-ener(i,mode)/g_11c))
238 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
239 ENDIF
240 dam(i,mode) =
min(dam(i,mode),one)
241
242 dfmax(i) =
max(dfmax(i),dam(i,mode))
243 ENDIF
244
245 IF (dam(i,mode) >= one) THEN
246 failmod = failmod + 1
247 fmode(i,mode) = 1
248 ENDIF
249
250
251 mode = 5
252 IF (-signyy(i) > sigma_22c .AND. signyy(i) < zero .AND. dam(i,mode) < one) THEN
253
254 IF (ishap22c == 1) THEN
255 dam(i,mode) = dam(i,mode) + le(i)*abs(depsyy(i))*sigma_22c/(two*g_22c)
256 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
257 ! -> exponential stress softening
258 ELSEIF (ishap22c == 2) THEN
259 ener(i,mode) = ener(i,mode) + abs(signyy(i))*le(i)*abs(depsyy(i))
260 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
261 dam(i,mode) = (one - exp(-ener(i,mode)/g_22c))
262 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
263 ENDIF
264 dam(i,mode) =
min(dam(i,mode),one)
265
266 dfmax(i) =
max(dfmax(i),dam(i,mode))
267 ENDIF
268
269 IF (dam(i,mode) >= one) THEN
270 failmod = failmod + 1
271 fmode(i,mode) = 1
272 ENDIF
273
274
275 mode = 6
276 IF (-signxy(i) > sigma_12c .AND. signxy(i) < zero .AND. dam(i,mode) < one) THEN
277
278 IF (ishap12c == 1) THEN
279 dam(i,mode) = dam(i,mode) + le(i)*abs(depsxy(i))*sigma_12c/(four*g_12c)
280 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
281
282 ELSEIF (ishap12c == 2) THEN
283 ener(i,mode) = ener(i,mode) + abs(signxy(i))*le(i)*half*abs(depsxy(i))
284 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
285 dam(i,mode) = (one - exp(-ener(i,mode)/g_12c))
286 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
287 ENDIF
288 dam(i,mode) =
min(dam(i,mode),one)
289
290 dfmax(i) =
max(dfmax(i),dam(i,mode))
291 ENDIF
292
293 IF (dam(i,mode) >= one) THEN
294 failmod = failmod + 1
295 fmode(i,mode) = 1
296 ENDIF
297
298
299 mode = 7
300 IF (signzz(i) > sigma_33t .AND. signzz(i) > zero .AND. dam(i,mode) < one) THEN
301
302 IF (ishap33t == 1) THEN
303 dam(i,mode) = dam(i,mode) + le(i)*depszz(i)*sigma_33t/(two*g_33t)
304 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
305
306 ELSEIF (ishap33t == 2) THEN
307 ener(i,mode) = ener(i,mode) + signzz(i)*le(i)*depszz(i)
308 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
309 dam(i,mode) = (one - exp(-ener(i,mode)/g_33t))
310 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
311 ENDIF
312 dam(i,mode) =
min(dam(i,mode),one)
313
314 dfmax(i) =
max(dfmax(i),dam(i,mode))
315 ENDIF
316
317 IF (dam(i,mode) >= one) THEN
318 failmod = failmod + 1
319 fmode(i,mode) = 1
320 ENDIF
321
322
323 mode = 8
324 IF (-signzz(i) > sigma_33c .AND. signzz(i) < zero .AND. dam(i,mode) < one) THEN
325
326 IF (ishap33c == 1) THEN
327 dam(i,mode) = dam(i,mode) + le(i)*abs(depszz(i))*sigma_33c/(two*g_33c)
328 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
329
330 ELSEIF (ishap33c == 2) THEN
331 ener(i,mode) = ener(i,mode) + abs(signzz(i))*le(i)*abs(depszz(i))
332 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
333 dam(i,mode) = (one - exp(-ener(i,mode)/g_33c))
334 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
335 ENDIF
336 dam(i,mode) =
min(dam(i,mode),one)
337
338 dfmax(i) =
max(dfmax(i),dam(i,mode))
339 ENDIF
340
341 IF (dam(i,mode) >= one) THEN
342 failmod = failmod + 1
343 fmode(i,mode) = 1
344 ENDIF
345
346
347 mode = 9
348 IF (signyz(i) > sigma_23t .AND. signyz(i) > zero .AND. dam(i,mode) < one) THEN
349
350 IF (ishap23t == 1) THEN
351 dam(i,mode) = dam(i,mode) + le(i)*depsyz(i)*sigma_23t/(four*g_23t)
352 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
353
354 ELSEIF (ishap23t == 2) THEN
355 ener(i,mode) = ener(i,mode) + signyz(i)*le(i)*half*depsyz(i)
356 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
357 dam(i,mode) = (one - exp(-ener(i,mode)/g_23t))
358 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
359 ENDIF
360 dam(i,mode) =
min(dam(i,mode),one)
361
362 dfmax(i) =
max(dfmax(i),dam(i,mode))
363 ENDIF
364
365 IF (dam(i,mode) >= one) THEN
366 failmod = failmod + 1
367 fmode(i,mode) = 1
368 ENDIF
369
370
371 mode = 10
372 IF (-signyz(i) > sigma_23c .AND. signyz(i) < zero .AND. dam(i,mode) < one) THEN
373
374 IF (ishap23c == 1) THEN
375 dam(i,mode) = dam(i,mode) + le(i)*abs(depsyz(i))*sigma_23c/(four*g_23c)
376 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
377
378 ELSEIF (ishap23c == 2) THEN
379 ener(i,mode) = ener(i,mode) + abs(signyz(i))*le(i)*half*abs(depsyz(i))
380 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
381 dam(i,mode) = (one - exp(-ener(i,mode)/g_23c))
382 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
383 ENDIF
384 dam(i,mode) =
min(dam(i,mode),one)
385
386 dfmax(i) =
max(dfmax(i),dam(i,mode))
387 ENDIF
388
389 IF (dam(i,mode) >= one) THEN
390 failmod = failmod + 1
391 fmode(i,mode) = 1
392 ENDIF
393
394
395 mode = 11
396 IF (signzx(i) > sigma_31t .AND. signzx(i) > zero .AND. dam(i,mode) < one) THEN
397
398 IF (ishap31t == 1) THEN
399 dam(i,mode) = dam(i,mode) + le(i)*depszx(i)*sigma_31t/(four*g_31t)
400 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
401
402 ELSEIF (ishap31t == 2) THEN
403 ener(i,mode) = ener(i,mode) + signzx(i)*le(i)*half*depszx(i)
404 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
405 dam(i,mode) = (one - exp(-ener(i,mode)/g_31t))
406 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
407 ENDIF
408 dam(i,mode) =
min(dam(i,mode),one)
409
410 dfmax(i) =
max(dfmax(i),dam(i,mode))
411 ENDIF
412
413 IF (dam(i,mode) >= one) THEN
414 failmod = failmod + 1
415 fmode(i,mode) = 1
416 ENDIF
417
418
419 mode = 12
420 IF (-signzx(i) > sigma_31c .AND. signzx(i) < zero .AND. dam(i,mode) < one) THEN
421
422 IF (ishap31c == 1) THEN
423 dam(i,mode) = dam(i,mode) + le(i)*abs(depszx(i))*sigma_31c/(four*g_31c)
424 dam(i,mode) =
max(dam(i,mode),uvar(i,mode))
425
426 ELSEIF (ishap31c == 2) THEN
427 ener(i,mode) = ener(i,mode) + abs(signzx(i))*le(i)*half*abs(depszx(i))
428 ener(i,mode) =
max(ener(i,mode),uvar(i,mode+12))
429 dam(i,mode) = (one - exp(-ener(i,mode)/g_31c))
430 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
431 ENDIF
432 dam(i,mode) =
min(dam(i,mode),one)
433
434 dfmax(i) =
max(dfmax(i),dam(i,mode))
435 ENDIF
436
437 IF (dam(i,mode) >= one) THEN
438 failmod = failmod + 1
439 fmode(i,mode) = 1
440 ENDIF
441
442
443 IF (failmod >= nmod) THEN
444 loff(i) = zero
445 nindx = nindx + 1
446 indx(nindx) = i
447 noff(i) = noff(i) + 1
448 IF (noff(i) >= failip) THEN
449 off(i) = zero
450 tdel(i) = time
451 nindx2 = nindx2 + 1
452 indx2(nindx2) = i
453 ENDIF
454 ENDIF
455 ENDIF
456 ENDDO
457
458
459
460
461 DO i = 1, nel
462 IF ((loff(i) == one) .AND. (off(i) == one)) THEN
463
464 IF (signxx(i) >= zero) THEN
465 dmg_scale(i,1) = one - dam(i,1)
466 ELSE
467 dmg_scale(i,1) = one - dam(i,4)
468 ENDIF
469
470 IF (signyy(i) >= zero) THEN
471 dmg_scale(i,2) = one - dam(i,2)
472 ELSE
473 dmg_scale(i,2) = one - dam(i,5)
474 ENDIF
475
476 IF (signzz(i) >= zero) THEN
477 dmg_scale(i,3) = one - dam(i,7)
478 ELSE
479 dmg_scale(i,3) = one - dam(i,8)
480 ENDIF
481
482 IF (signxy(i) >= zero) THEN
483 dmg_scale(i,4) = one - dam(i,3)
484 ELSE
485 dmg_scale(i,4) = one - dam(i,6)
486 ENDIF
487
488 IF (signyz(i) >= zero) THEN
489 dmg_scale(i,5) = one - dam(i,9)
490 ELSE
491 dmg_scale(i,5) = one - dam(i,10)
492 ENDIF
493
494 IF (signzx(i) >= zero) THEN
495 dmg_scale(i,6) = one - dam(i,11)
496 ELSE
497 dmg_scale(i,6) = one - dam(i,12)
498 ENDIF
499 ENDIF
500 ENDDO
501
502
503
504
505 DO j = 1,12
506 DO i = 1,nel
507
508 uvar(i,j) = dam(i,j)
509
510 uvar(i,j+12) = ener(i,j)
511 ENDDO
512 ENDDO
513
514
515
516
517 IF (nindx > 0) THEN
518 DO j=1,nindx
519 i = indx(j)
520#include "lockon.inc"
521 WRITE(iout, 1000) ngl(i),ipg,ilay
522 WRITE(istdo,1000) ngl(i),ipg,ilay
523 IF (fmode(i,1) == 1) WRITE(iout,2000) '1 - TRACTION XX'
524 IF (fmode(i,2) == 1) WRITE(iout,2000) '2 - TRACTION YY'
525 IF (fmode(i,3) == 1) WRITE(iout,2000) '3 - POSITIVE SHEAR XY'
526 IF (fmode(i,4) == 1) WRITE(iout,2000) '4 - COMPRESSION XX'
527 IF (fmode(i,5) == 1) WRITE(iout,2000) '5 - COMPRESSION YY'
528 IF (fmode(i,6) == 1) WRITE(iout,2000) '6 - NEGATIVE SHEAR XY'
529 IF (fmode(i,7) == 1) WRITE(iout,2000) '7 - TRACTION ZZ'
530 IF (fmode(i,8) == 1) WRITE(iout,2000) '8 - COMPRESSION ZZ'
531 IF (fmode(i,9) == 1) WRITE(iout,2000) '9 - POSITIVE SHEAR YZ'
532 IF (fmode(i,10) == 1) WRITE(iout,2000) '10 - NEGATIVE SHEAR YZ'
533 IF (fmode(i,11) == 1) WRITE(iout,2000) '11 - POSITIVE SHEAR ZX'
534 IF (fmode(i,12) == 1) WRITE(iout,2000) '12 - NEGATIVE SHEAR ZX'
535#include "lockoff.inc"
536 ENDDO
537 ENDIF
538
539 IF (nindx2 > 0) THEN
540 DO j=1,nindx2
541 i = indx2(j)
542#include "lockon.inc"
543 WRITE(iout, 3000) ngl(i),time
544 WRITE(istdo,3000) ngl(i),time
545#include "lockoff.inc"
546 ENDDO
547 ENDIF
548
549 1000 FORMAT(1x,'FAILURE (ORTHENERG) OF SOLID ELEMENT ',i10,1x,',GAUSS PT',
550 . i5,1x,',LAYER',i3)
551 2000 FORMAT(1x,'MODE',1x,a)
552 3000 FORMAT(1x,'-- RUPTURE OF SOLID ELEMENT : ',i10,1x,'AT TIME : ',1pe12.4)
553