43
44
45
46
47
48#include "implicit_f.inc"
49
50
51
52#include "units_c.inc"
53#include "comlock.inc"
54#include "com_xfem1.inc"
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82 INTEGER NEL,NUPARAM,NUVAR,NFUNC,IPT,IPTT,NPTOT,IXFEM,IXEL,ILAY,NPTT
83 INTEGER ELCRKINI(NXLAYMAX,NEL),NOFF(NEL),IFUNC(NFUNC),NGL(NEL),FLD_IDX(NEL)
85 . time,uparam(nuparam),signxx(nel),signyy(nel),
86 . signzz(nel),signxy(nel),signyz(nel),signzx(nel),ssp(nel),
87 . epsxx(nel),epsyy(nel),epsxy(nel),epsyz(nel),epszx(nel),
88 . tens(nel,5),dfmax(nel,lf_dammx),tdel(nel),npttf(nel),dam(nel)
89 INTEGER, INTENT(IN) :: NIPARAM, LF_DAMMX
90 INTEGER, DIMENSION(NIPARAM), INTENT(IN) :: IPARAM
91 my_real,
DIMENSION(NEL),
INTENT(IN) :: pla,depsxx,depsyy,depsxy
93
94
95
96 my_real uvar(nel,nuvar),off(nel),offl(nel)
97
98
99
100 INTEGER NPF(*)
101 my_real finter , finterfld ,tf(*)
102 EXTERNAL finter
103
104
105
106
107
108
109
110
111
112
113 INTEGER I,II,J,K,L,IENG,IAD,NINDX,NINDXP,ISHELL,LAYXFEM,LENF
114 INTEGER INDX(NEL),INDXP(NEL),RFLAG(NEL),RFLAGP(NEL)
116 . rani,r1,r2,dadv,sigmax,s1,s2,ss,q,dydx,e12
117 my_real,
ALLOCATABLE,
DIMENSION(:) :: xf
119 . tba(nel),tbk(nel),sigr(nel),emaj(nel),emin(nel),
120 . em(nel),em_dcrit(nel),em_dadv(nel),demaj(nel),demin(nel),
121 . beta(nel),
alpha,fcut
122 CHARACTER (LEN=3) :: XCHAR
123
124 iad = (ipt-1)*nel
125 xchar = ' '
126 DO i=1,nel
127 rflag(i) = 0
128 rflagp(i)= 0
129 indxp(i) = 0
130 indx(i) = 0
131 ENDDO
132 nindx = 0
133 nindxp = 0
134
135
136
137
138
139
140 ishell = iparam(1)
141 ieng = iparam(3)
142
143 dadv = uparam(2)
144 rani = uparam(3)
145 fcut = uparam(5)
146 IF (uparam(6) > zero) THEN
148 ELSE
149 alpha = two*pi*fcut*dt1/(one + two*pi*fcut*dt1)
150 ENDIF
151
152 layxfem = ixfem
153 IF (ixfem == 1 .and. ishell == 1) ishell=2
154
155 DO i=1,nel
156 em_dcrit(i)= zero
157 em_dadv(i) = zero
158 END DO
159
160
161
162
163 DO i = 1,nel
164
165 e12= half*epsxy(i)
166 s1 = half*(epsxx(i) + epsyy(i))
167 s2 = half*(epsxx(i) - epsyy(i))
168 q = sqrt(s2**2 + e12**2)
169 emaj(i) = s1 + q
170 emin(i) = s1 - q
171 IF (emin(i) >= emaj(i)) THEN
172 ss = emin(i)
173 emin(i) = emaj(i)
174 emaj(i) = ss
175 ENDIF
176
177 e12 = half*depsxy(i)
178 s1 = half*(depsxx(i) + depsyy(i))
179 s2 = half*(depsxx(i) - depsyy(i))
180 q = sqrt(s2**2 + e12**2)
181 demaj(i) = s1 + q
182 demin(i) = s1 - q
183
184 demaj(i) =
alpha*demaj(i) + (one -
alpha)*uvar(i,2)
185 demin(i) =
alpha*demin(i) + (one -
alpha)*uvar(i,3)
186 beta(i) = demin(i)/sign(
max(abs(demaj(i)),em20),demaj(i))
187 uvar(i,2) = demaj(i)
188 uvar(i,3) = demin(i)
189 IF (ieng == 2) THEN
190 dfmax(i,4) = beta(i)
191 ENDIF
192 ENDDO
193
194
195
196
197
198 IF (ieng == 1) THEN
199 ii = npf(ifunc(1))
200 lenf = npf(ifunc(1)+ 1) - npf(ifunc(1))
201 ALLOCATE(xf(lenf))
202 DO i = 1,lenf
203 xf(i) = log(tf(ii + i-1) + one)
204 ENDDO
205
206 DO i = 1,nel
207 em(i) = finterfld(emin(i),lenf,xf)
208 dam(i) = emaj(i) / em(i)
209 dfmax(i,2) = dam(i)
210 dfmax(i,1) =
min(one, dam(i))
211 ENDDO
212
213 ELSE
214
215 IF (ieng == 0) THEN
216 DO i = 1,nel
217 em(i) = finter(ifunc(1),emin(i),npf,tf,dydx)
218 dam(i) = emaj(i) / em(i)
219 dfmax(i,2) = dam(i)
220 dfmax(i,1) =
min(one, dam(i))
221 ENDDO
222
223 ELSEIF (ieng == 2) THEN
224 DO i = 1,nel
225 em(i) = finter(ifunc(1),emin(i),npf,tf,dydx)
226 dam(i) =
max(pla(i) / em(i),dam(i))
227 dfmax(i,2) = dam(i)
228 dfmax(i,1) =
min(one, dam(i))
229 ENDDO
230 ENDIF
231 ENDIF
232
233
234
235
236 r1 = two*em02
237 r2 = rani/(rani+one)
238 IF (ieng < 2) THEN
239 DO i = 1,nel
240 IF (emaj(i) >= em(i)) THEN
241 fld_idx(i) = 6
242 ELSEIF (emaj(i) >= em(i) - em01) THEN
243 fld_idx(i) = 5
244 ELSEIF (emaj(i)**2 + emin(i)**2 < r1**2) THEN
245 fld_idx(i) = 1
246 ELSEIF (emaj(i) >= abs(emin(i))) THEN
247 fld_idx(i) = 4
248 ELSEIF (emaj(i) >= r2*abs(emin(i))) THEN
249 fld_idx(i) = 3
250 ELSE
251 fld_idx(i) = 2
252 ENDIF
253 dfmax(i,3) = fld_idx(i)
254 ENDDO
255 ELSE
256 DO i = 1,nel
257 IF (pla(i) >= em(i)) THEN
258 fld_idx(i) =
max(6,fld_idx(i))
259 ELSEIF (pla(i) >= em(i) - em01) THEN
260 fld_idx(i) =
max(5,fld_idx(i))
261 ELSEIF (pla(i)**2 + beta(i)**2 < r1**2) THEN
262 fld_idx(i) =
max(1,fld_idx(i))
263 ELSEIF (pla(i) >= abs(beta(i))) THEN
264 fld_idx(i) =
max(4,fld_idx(i))
265 ELSEIF (pla(i) >= r2*abs(beta(i))) THEN
266 fld_idx(i) =
max(3,fld_idx(i))
267 ELSE
268 fld_idx(i) =
max(2,fld_idx(i))
269 ENDIF
270 dfmax(i,3) = fld_idx(i)
271 ENDDO
272 ENDIF
273
274 DO i = 1,nel
275 tens(i,1) = signxx(i)
276 tens(i,2) = signyy(i)
277 tens(i,3) = signxy(i)
278 tens(i,4) = signyz(i)
279 tens(i,5) = signzx(i)
280 em_dcrit(i)= em(i)
281 IF (nfunc > 1 .and. ifunc(2) > 0) THEN
282 em_dadv(i) = dadv*finter(ifunc(2),emin(i),npf,tf,dydx)
283 ELSE
284 em_dadv(i) = dadv*em_dcrit(i)
285 ENDIF
286 END DO
287
288 IF (ixel > 0) THEN
289 IF (ixel == 1) THEN
290 xchar = '1st'
291 ELSEIF (ixel == 2) THEN
292 xchar = '2nd'
293 ELSEIF (ixel == 3) THEN
294 xchar = '3rd'
295 ENDIF
296 ELSE
297 xchar = 'standard'
298 ENDIF
299
300 SELECT CASE (layxfem)
301
302 CASE (1)
303
304 DO i = 1,nel
305 IF (ieng == 2) emaj(i) = pla(i)
306 IF (off(i) == one) THEN
307 IF (uvar(i,1) == zero) THEN
308 IF (ixel == 0) THEN
309 IF (elcrkini(ilay,i)==0 .and. emaj(i)>=em_dcrit(i)) THEN
310 nindxp = nindxp+1
311 indxp(nindxp) = i
312 uvar(i,1) = one
313 offl(i) = zero
314 noff(i) = noff(i) + 1
315 npttf(i) = npttf(i) + one
316 rflagp(i) = 1
317 IF (int(npttf(i)) == nptt) THEN
318 nindx = nindx+1
319 indx(nindx) = i
320 elcrkini(ilay,i) = -1
321 rflag(i) = 1
322 ENDIF
323 IF (noff(i) == nptot) THEN
324 off(i) = four_over_5
325 tdel(i)= time
326 ENDIF
327 ELSEIF (elcrkini(ilay,i)==2 .and. emaj(i)>=em_dadv(i)) THEN
328 nindxp = nindxp + 1
329 indxp(nindxp) = i
330 uvar(i,1) = one
331 offl(i) = zero
332 noff(i) = noff(i) + 1
333 npttf(i) = npttf(i)+ one
334 rflagp(i) = 1
335 IF (int(npttf(i)) == nptt) THEN
336 nindx = nindx+1
337 indx(nindx) = i
338 elcrkini(ilay,i) = 1
339 rflag(i) = -1
340 ENDIF
341 IF (noff(i) == nptot) THEN
342 off(i) = four_over_5
343 tdel(i)= time
344 ENDIF
345 ENDIF
346 ELSEIF (emaj(i) >= em_dcrit(i)) THEN
347 nindxp = nindxp + 1
348 indxp(nindxp) = i
349 uvar(i,1) = one
350 offl(i) = zero
351 npttf(i) = npttf(i) + one
352 rflagp(i) = 1
353 IF (int(npttf(i)) == nptt) THEN
354 nindx = nindx+1
355 indx(nindx) = i
356 off(i) = four_over_5
357 tdel(i) = time
358 rflag(i) = 3
359 ENDIF
360 ENDIF
361 ELSEIF (uvar(i,1) == one) THEN
362 signxx(i) = zero
363 signyy(i) = zero
364 signxy(i) = zero
365 signyz(i) = zero
366 signzx(i) = zero
367 ENDIF
368 ENDIF
369 ENDDO
370 IF (nindxp > 0) THEN
371 DO j=1,nindxp
372 i = indxp(j)
373 signxx(i) = zero
374 signyy(i) = zero
375 signxy(i) = zero
376 signyz(i) = zero
377 signzx(i) = zero
378#include "lockon.inc"
379 IF (rflagp(i) == 1) WRITE(iout, 3800)ngl(i),ilay,iptt
380 IF (rflagp(i) == 1) WRITE(istdo,3900)ngl(i),ilay,iptt,time
381#include "lockoff.inc"
382 ENDDO
383 ENDIF
384 IF (nindx > 0) THEN
385 DO j=1,nindx
386 i = indx(j)
387#include "lockon.inc"
388
389 IF (rflag(i) == 1)WRITE
390 IF (rflag(i) == 1) WRITE(istdo,3100) ngl(i),ilay,time
391
392 IF (rflag(i) == -1) WRITE(iout ,3200) ngl(i),ilay
393 IF (rflag(i) == -1) WRITE(istdo,3300) ngl(i),ilay,time
394
395 IF (rflag(i) == 3) WRITE(iout, 3400) xchar,ngl(i),ilay
396 IFWRITE(istdo,3500) xchar,ngl(i),ilay,time
397#include "lockoff.inc"
398 ENDDO
399 ENDIF
400
401 CASE (2)
402
403 IF (ishell == 1) THEN
404
405 IF (layxfem > 0) THEN
406 DO i = 1,nel
407 IF (ieng == 2) emaj(i) = pla(i)
408 IF (off(i) == one) THEN
409 IF (ixel == 0) THEN
410 IF (elcrkini(ilay,i) == 0 .and.
411 . emaj(i) >= em_dcrit(i)) THEN
412 elcrkini(ilay,i) = -1
413 off(i) = four_over_5
414 nindx=nindx+1
415 indx(nindx)=i
416 rflag(i) = 1
417 tdel(i)= time
418 ELSEIF (elcrkini(ilay,i) == 2 .and.
419 . emaj(i) >= em_dadv(i)) THEN
420 elcrkini(ilay,i) = 1
421 off(i) = four_over_5
422 nindx=nindx+1
423 indx(nindx)=i
424 rflag(i) = -1
425 tdel(i)= time
426 ENDIF
427 ELSE IF (emaj(i) >= em_dcrit(i)) then
428 off(i) = four_over_5
429 nindx=nindx+1
430 indx(nindx)=i
431 rflag(i) = 3
432 ENDIF
433 ENDIF
434 ENDDO
435 ENDIF
436
437 IF (nindx > 0) THEN
438 DO j=1,nindx
439 i=indx(j)
440#include "lockon.inc"
441 WRITE(iout, 4800)ngl(i),iptt
442 WRITE(istdo,4900)ngl(i),iptt,time
443
444 IF (rflag(i) == 1) WRITE(iout ,4000) ngl(i)
445 IF (rflag(i) == 1) WRITE(istdo,4100) ngl(i),time
446
447 IF (rflag(i) == -1) WRITE(iout, 4200) ngl(i)
448 IF (rflag(i) == -1) WRITE(istdo,4300) ngl(i),time
449
450 IF (rflag(i) == 3) WRITE(iout, 4400) xchar,ngl(i)
451 IF (rflag(i) == 3) WRITE(istdo,4500) xchar,ngl(i),time
452#include "lockoff.inc"
453 ENDDO
454 ENDIF
455
456 ELSEIF (ishell == 2)THEN
457
458 DO i = 1,nel
459 IF (ieng == 2) emaj(i) = pla(i)
460 IF (off(i) == one) THEN
461 IF (uvar(i,1) == zero) THEN
462 IF (ixel == 0) THEN
463 IF (elcrkini(ilay,i) == 0 .and. emaj(i) >= em_dcrit(i)) THEN
464 uvar(i,1) = one
465 nindx=nindx+1
466 indx(nindx)=i
467 nindxp = nindxp + 1
468 indxp(nindxp) = i
469 rflagp(i) = 1
470 noff(i) = noff(i) + 1
471 IF (noff(i) == nptot) THEN
472 elcrkini(ilay,i) = -1
473 off(i) = four_over_5
474 rflag(i) = 1
475 tdel(i)= time
476 ENDIF
477 ELSEIF (elcrkini(ilay,i) == 2 .and. emaj(i) >= em_dadv(i)) THEN
478 uvar(i,1) = one
479 nindx=nindx+1
480 indx(nindx)=i
481 nindxp = nindxp + 1
482 indxp(nindxp) = i
483 rflagp(i) = 1
484 noff(i) = noff(i) + 1
485 IF (noff(i) == nptot) THEN
486 off(i) = four_over_5
487 elcrkini(ilay,i) = 1
488 rflag(i) = -1
489 tdel(i)= time
490 ENDIF
491 ENDIF
492 ELSEIF (emaj(i) >= em_dcrit(i)) THEN
493 uvar(i,1) = one
494 nindx = nindx+1
495 indx(nindx) = i
496 nindxp = nindxp + 1
497 indxp(nindxp) = i
498 rflagp(i) = 1
499 noff(i) = noff(i) + 1
500 IF (noff(i) == nptot) THEN
501 off(i) = four_over_5
502 rflag(i) = 3
503 ENDIF
504 ENDIF
505 ELSEIF (uvar(i,1) == one) THEN
506 signxx(i) = zero
507 signyy(i) = zero
508 signxy(i) = zero
509 signyz(i) = zero
510 signzx(i) = zero
511 ENDIF
512 ENDIF
513 ENDDO
514
515 IF (nindxp > 0) THEN
516 DO j=1,nindxp
517 i = indxp(j)
518 signxx(i) = zero
519 signyy(i) = zero
520 signxy(i) = zero
521 signyz(i) = zero
522 signzx(i) = zero
523#include "lockon.inc"
524 IF (rflagp(i) == 1) WRITE(iout, 4800)ngl(i),iptt
525 IF (rflagp(i) == 1) WRITE(istdo,4900)ngl(i),iptt,time
526#include "lockoff.inc"
527 ENDDO
528 ENDIF
529 IF (nindx > 0) THEN
530 DO j=1,nindx
531 i = indx(j)
532#include "lockon.inc"
533
534 IF (rflag(i) == 1) WRITE(iout ,4000) ngl(i)
535 IF (rflag(i) == 1) WRITE(istdo,4100) ngl(i),time
536
537 IF (rflag(i) == -1) WRITE(iout, 4200) ngl(i)
538 IF (rflag(i) == -1) WRITE(istdo,4300) ngl(i),time
539
540 IF (rflag(i) == 3) WRITE(iout, 4400) xchar,ngl(i)
541 IF (rflag(i) == 3) WRITE(istdo,4500) xchar,ngl(i),time
542#include "lockoff.inc"
543 ENDDO
544 ENDIF
545 ENDIF
546
547 END SELECT
548
549 IF (ieng == 1) DEALLOCATE(xf)
550
551 2000 FORMAT(1x,'FOR SHELL ELEMENT (FLD)',i10,1x,'LAYER',i3,':',/
552 . 1x, 'STRESS TENSOR SET TO ZERO')
553 2100 FORMAT(1x,'FOR SHELL ELEMENT (FLD)',i10,1x,'LAYER',i3,':',/,
554 . 1x, 'STRESS TENSOR SET TO ZERO',1x,'AT TIME :',1pe20.13)
555 2400 FORMAT(1x,1pg20.13,' % OF THICKNESS OF SHELL BROKEN ')
556 2410 FORMAT(1x,1pg20.13,' % OF THICKNESS OF SHELL ',i10,' BROKEN '
557 2500 FORMAT(1x,' LOWER SKIN -> UPPER SKIN ')
558 2600 FORMAT(1x,' UPPER SKIN -> LOWER SKIN ')
559
560 3000 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT',i10,1x,'LAYER',i3)
561 3100 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT',i10,1x,'LAYER',i3,/
562 . 1x,'AT TIME :',1pe12.4)
563 3200 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT',i10,' LAYER',i3)
564 3300 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT',i10,' LAYER',i3/
565 . 1x,'AT TIME :',1pe12.4)
566 3400 FORMAT(1x,'DELETE ',a4,' PHANTOM ELEMENT, SHELL ID=',i10,' LAYER',i3)
567 3500 FORMAT(1x,'DELETE ',a4,' PHANTOM ELEMENT, SHELL ID=',i10,' LAYER',i3,/
568 . 1x,'AT TIME :',1pe12.4)
569 3800 FORMAT(1x,'fld failure in shell',I10,1X,'layer',I3,1X,'int point',I2)
570 3900 FORMAT(1X,'fld failure in shell',I10,1X,'layer',I3,1X,'int point',I2,/
571 . 1X,'at time :',1PE12.4)
572
573 4000 FORMAT(1X,'crack initialization in shell element',I10)
574 4100 FORMAT(1X,'crack initialization in shell element',I10,/
575 . 1X,'at time :',1PE12.4)
576 4200 FORMAT(1X,'crack advancement in shell element',I10)
577 4300 FORMAT(1X,'crack advancement in shell element',I10,/
578 . 1X,'at time :',1PE12.4)
579 4400 FORMAT(1X,'delete
',A4,' phantom element, shell
id=
',I10)
580 4500 FORMAT(1X,'delete
',A4,' phantom element, shell
id=
',I10,/
581 . 1X,'at time :',1PE12.4)
582 4800 FORMAT(1X,'fld failure in shell',I10,1X,'int point',I2)
583 4900 FORMAT(1X,'fld failure in shell',I10,1X,'int point',I2,1X,'at time :',1PE12.4)
584
585 RETURN