41
42
43
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 INTEGER NEL,NPARAM,NUVAR,IPT,NFUNC,IXFEM,IXEL,ILAY,NPTOT
64 INTEGER NGL(NEL),NOFF(NEL),IFUNC(NFUNC),
65 . ELCRKINI(NXLAYMAX,NEL)
66 my_real time,timestep(nel),uparam(*),dpla(nel),epsp(nel),
67 . tstar(nel),aldt(nel)
68
69
70
73 . uvar(nel,nuvar),off(nel),offl(nel),
74 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
75 . tens(nel,5),dfmax(nel),tdel(nel)
76
77
78
79 INTEGER NPF(*)
81 EXTERNAL finter
82
83
84
85 INTEGER :: I,J,J1,J2,K,L,NINDX,ISHELL,LAYXFEM,NRATE,DMG_FLAG,
86 . IFUN_EL,IFUN_TEMP,NF_LOC
87 INTEGER, DIMENSION(NEL) :: INDX,RFLAG
88 INTEGER ,DIMENSION(NFUNC) :: IFUN_STR
89
90 my_real :: dp,p,sigm,sigm_ps,svm,ef1,ef2,df,fac,depsf,lambda,
91 . cc,bb,cr,orm,ss1,ss2,yy,yy_n,dadv,
92 . dcrit,el_ref,sc_el,sc_temp,dd,dn
93 my_real,
DIMENSION(NEL) :: epsf,dmg_scale
94 my_real,
DIMENSION(NFUNC) :: yfac,rate
95 CHARACTER (LEN=3) :: XCHAR
96
97
98
99
100
101
102 dmg_flag = 1
103 indx = 0
104 sigm_ps = one/sqrt(three)
105
106 IF (uvar(1,2) == zero) THEN
107 uvar(1:nel,2) = aldt(1:nel)
108 ENDIF
109
110 ishell = int(uparam(2))
111 dcrit = uparam(4)
112 dd = uparam(5)
113 dn = uparam(6)
114 sc_temp = uparam(7)
115 sc_el = uparam(8)
116 el_ref = uparam(9)
117 dadv = uparam(11)
118 nrate = nfunc - 2
119 yfac(1:nrate) = uparam(11+1 :11+nrate)
120 rate(1:nrate) = uparam(11+1+nrate:11+nrate*2)
121
122 layxfem = ixfem
123 IF (layxfem == 1 .and. ishell == 1) ishell=2
124
125 nindx = 0
126 rflag = 0
127
128 DO i=1,nel
129 tens(i,1) = signxx(i)
130 tens(i,2) = signyy(i)
131 tens(i,3) = signxy(i)
132 tens(i,4) = signyz(i)
133 tens(i,5) = signzx(i)
134 END DO
135
136 IF (ixel > 0) THEN
137 IF (ixel == 1) THEN
138 xchar = '1st'
139 ELSEIF (ixel == 2) THEN
140 xchar = '2nd'
141 ELSEIF (ixel == 3) THEN
142 xchar = '3rd'
143 ENDIF
144 ELSE
145 xchar = ' '
146 ENDIF
147
148
149
150 ifun_str(1:nrate) = ifunc(1:nrate)
151
152 ifun_el = ifunc(nrate+1)
153 ifun_temp = ifunc(nrate+2)
154
155
156
157 DO i=1,nel
158 j1 = 1
159 DO k=2, nrate-1
160 IF (epsp(i) > rate(j1)) j1 = k
161 ENDDO
162 p = third*(signxx(i) + signyy(i))
163 svm = sqrt(signxx(i)*signxx(i) + signyy(i)*signyy(i)
164 . - signxx(i)*signyy(i) + three*signxy(i)*signxy(i))
165 sigm = p /
max(em20,svm)
166
167 IF (nrate > 1) THEN
168 j2 = j1+1
169 ef1 = yfac(j1)*finter(ifun_str(j1),sigm,npf,tf,df)
170 ef2 = yfac(j2)*finter(ifun_str(j2),sigm,npf,tf,df)
171 fac = (epsp(i) - rate(j1))/(rate(j2) - rate(j1))
173 epsf(i) =
max(ef1 + fac*(ef2 - ef1), em20)
174 ELSE
175 epsf(i) = yfac(j1)*finter(ifun_str(j1),sigm,npf,tf,df)
176 ENDIF
177 ENDDO
178
179 DO i=1,nel
180
181 IF (ifun_el > 0) THEN
182 lambda = uvar(i,2) / el_ref
183 fac = sc_el*finter(ifun_el,lambda,npf,tf,df)
184 epsf(i) = epsf(i)* fac
185 ENDIF
186
187 IF (ifun_temp > 0) THEN
188 fac = sc_temp*finter(ifun_temp,tstar(i),npf,tf,df)
189 epsf(i) = epsf(i)* fac
190 ENDIF
191 ENDDO
192
193 IF (ishell == 1) THEN
194 IF (ixfem == 1 .OR. ixfem == 2) THEN
195 DO i=1,nel
196 IF (ishell == 1 .AND. off(i)==one) THEN
197 IF(uvar(i,1) == zero) THEN
198 dp = one
199 ELSE
200
201 dp = dn*uvar(i,1)**(one-one/dn)
202 ENDIF
203 IF (epsf(i) > zero) uvar(i,1)=
204 . uvar(i,1)+dp*dpla(i)/epsf(i)
205 IF (ixel == 0) THEN
206 IF (elcrkini(ilay,i)==0) THEN
207 IF (uvar(i,1) >= dcrit) THEN
208 elcrkini(ilay,i) = -1
209 off(i) = four_over_5
210 nindx=nindx+1
211 indx(nindx)=i
212 rflag(i) = 1
213 tdel(i)= time
214 ENDIF
215 ELSEIF (elcrkini(ilay,i) == 2) THEN
216 IF (uvar(i,1) >= dadv) THEN
217 elcrkini(ilay,i) = 1
218 off(i) = four_over_5
219 nindx=nindx+1
220 indx(nindx)=i
221 rflag(i) = -1
222 tdel(i)= time
223 ENDIF
224 ENDIF
225 ELSEIF (uvar(i,1 )>= dcrit) THEN
226 off(i) = four_over_5
227 nindx=nindx+1
228 indx(nindx)=i
229 rflag(i) = 2
230 ENDIF
231 ENDIF
232 ENDDO
233 ENDIF
234
235 IF (nindx > 0) THEN
236 DO j=1,nindx
237 i=indx(j)
238#include "lockon.inc"
239
240 IF (rflag(i)>0.AND.rflag(i)<2)
241 . WRITE(iout, 3800) ngl(i)
242 IF (rflag(i)>0.AND.rflag(i)<2)
243 . WRITE(istdo,3900) ngl(i),time
244
245 IF (rflag(i) < 0) WRITE(iout, 4000) ngl(i)
246 IF (rflag(i) < 0) WRITE(istdo,4100) ngl(i),time
247
248 IF (rflag(i) > 1) WRITE(iout, 4200)xchar,ngl(i)
249 IF (rflag(i) > 1) WRITE(istdo,4300)xchar,ngl(i),time
250#include "lockoff.inc"
251 ENDDO
252 ENDIF
253 ENDIF
254
255 IF (ishell > 1) THEN
256 IF (ixfem == 1) THEN
257 DO i=1,nel
258 IF (off(i) == one)THEN
259 IF (uvar(i,1) < dcrit) THEN
260 IF(uvar(i,1) == zero) THEN
261 dp = one
262 ELSE
263
264 dp = dn*uvar(i,1)**(one-one/dn)
265 ENDIF
266 IF (epsf(i) > zero) uvar(i,1)=
267 . uvar(i,1)+dp*dpla(i)/epsf(i)
268 IF (ixel == 0) THEN
269 IF (elcrkini(ilay,i) == 0 .AND.
270 . uvar(i,1) >= dcrit) THEN
271 IF (ishell == 2) THEN
272 signxx(i) = zero
273 signyy(i) = zero
274 signxy(i) = zero
275 signyz(i) = zero
276 signzx(i) = zero
277 ENDIF
278 nindx=nindx+1
279 indx(nindx)=i
280 elcrkini(ilay,i) = -1
281 noff(i) = noff(i) + 1
282 IF (noff(i) == nptot) THEN
283 off(i) = four_over_5
284 tdel(i)= time
285 ENDIF
286 rflag(i) = 1
287 ELSEIF (elcrkini(ilay,i) == 2 .AND.
288 . uvar(i,1) >= dadv) THEN
289 IF (ishell == 2) THEN
290 signxx(i) = zero
291 signyy(i) = zero
292 signxy(i) = zero
293 signyz(i) = zero
294 signzx(i) = zero
295 ENDIF
296 nindx=nindx+1
297 indx(nindx)=i
298 elcrkini(ilay,i) = 1
299 noff(i) = noff(i) + 1
300 IF(dadv < dcrit) uvar(i,1) = dcrit
301 IF (noff(i) == nptot) THEN
302 off(i) = four_over_5
303 tdel(i)= time
304 ENDIF
305 rflag(i) = -1
306 ENDIF
307 ELSEIF (uvar(i,1) >= dcrit) THEN
308 IF (ishell == 2) THEN
309 signxx(i) = zero
310 signyy(i) = zero
311 signxy(i) = zero
312 signyz(i) = zero
313 signzx(i) = zero
314 ENDIF
315 nindx=nindx+1
316 indx(nindx)=i
317 noff(i) = noff(i) + 1
318 rflag(i) = 3
319 IF (noff(i) == nptot) THEN
320 off(i) = four_over_5
321 rflag(i) = 4
322 ENDIF
323 ENDIF
324 ELSEIF (ishell == 2) THEN
325 signxx(i) = zero
326 signyy(i) = zero
327 signxy(i) = zero
328 signyz(i) = zero
329 signzx(i) = zero
330 ENDIF
331 ENDIF
332 ENDDO
333 ELSEIF (ixfem == 2) THEN
334 DO i=1,nel
335 IF (off(i)==one .AND. (ishell==2 .OR. ishell==3))THEN
336 IF (uvar(i,1) < dcrit) THEN
337 IF(uvar(i,1) == zero) THEN
338 dp = one
339 ELSE
340
341 dp = dn*uvar(i,1)**(one-one/dn)
342 ENDIF
343
344 IF (epsf(i) > zero) uvar(i,1)=
345 . uvar(i,1)+dp*dpla(i)/epsf(i)
346 IF (ixel == 0) THEN
347 IF (elcrkini(ilay,i) == 0 .AND.
348 . uvar(i,1) >= dcrit) THEN
349 IF (ishell == 2) THEN
350 signxx(i) = zero
351 signyy(i) = zero
352 signxy(i) = zero
353 signyz(i) = zero
354 signzx(i) = zero
355 ENDIF
356 nindx=nindx+1
357 indx(nindx)=i
358 noff(i) = noff(i) + 1
359 IF (noff(i) == nptot) THEN
360 off(i) = four_over_5
361 elcrkini(ilay,i) = -1
362 rflag(i) = 1
363 tdel(i)= time
364 ENDIF
365 ELSEIF (elcrkini(ilay,i) == 2 .AND.
366 . uvar(i,1) >= dadv) THEN
367 IF (ishell == 2) THEN
368 signxx(i) = zero
369 signyy(i) = zero
370 signxy(i) = zero
371 signyz(i) = zero
372 signzx(i) = zero
373 ENDIF
374 nindx=nindx+1
375 indx(nindx)=i
376 noff(i) = noff(i) + 1
377 IF(dadv < dcrit) uvar(i,1) = dcrit
378 IF (noff(i) == nptot) THEN
379 off(i) = four_over_5
380 elcrkini(ilay,i) = 1
381 rflag(i) = -1
382 tdel(i)= time
383 ENDIF
384 ENDIF
385 ELSEIF (uvar(i,1) >= dcrit) THEN
386 IF (ishell == 2) THEN
387 signxx(i) = zero
388 signyy(i) = zero
389 signxy(i) = zero
390 signyz(i) = zero
391 signzx(i) = zero
392 ENDIF
393 nindx=nindx+1
394 indx(nindx)=i
395 noff(i) = noff(i) + 1
396 IF (noff(i) == nptot) THEN
397 off(i) = four_over_5
398 rflag(i) = 4
399 ENDIF
400 ENDIF
401 ELSEIF (ishell == 2) THEN
402 signxx(i) = zero
403 signyy(i) = zero
404 signxy(i) = zero
405 signyz(i) = zero
406 signzx(i) = zero
407 ENDIF
408 ENDIF
409 ENDDO
410 ENDIF
411
412 IF (nindx > 0) THEN
413 DO j=1,nindx
414 i = indx(j)
415#include "lockon.inc"
416 IF(ixfem ==1)THEN
417
418 IF (rflag(i)>0.AND.rflag(i)<3)WRITE(iout,4600)ngl(i),ipt
419 IF (rflag(i)>0.AND.rflag(i)<3)WRITE(istdo,4700)
420 . ngl(i),ipt,time
421
422 IF (rflag(i) < 0) WRITE(iout, 4800) ngl(i),ipt
423 IF (rflag(i) < 0) WRITE(istdo,4900) ngl(i),ipt,time
424
425 IF (rflag(i) > 2) WRITE(iout, 4400)xchar,ngl(i),ipt
426 IF (rflag(i) > 2) WRITE(istdo,4500)xchar,ngl(i),ipt,time
427
428 IF (rflag(i) /= 0 .AND. ixel == 0)
429 . WRITE(iout, 2000) ngl(i),ipt
430 IF (rflag(i) /= 0.AND. ixel == 0)
431 . WRITE(istdo,2100) ngl(i),ipt,time
432 ELSEIF(ixfem ==2)THEN
433
434 IF (rflag(i)>0.AND.rflag(i)<3)WRITE(iout,3800)ngl(i)
435 IF (rflag(i)>0.AND.rflag(i)<3)WRITE(istdo,3900)
436 . ngl(i),time
437
438 IF (rflag(i) < 0) WRITE(iout, 4000) ngl(i)
439 IF (rflag(i) < 0) WRITE(istdo,4100) ngl(i),time
440
441 IF (rflag(i) > 2) WRITE(iout, 4200)xchar,ngl(i)
442 IF (rflag(i) > 2) WRITE(istdo,4300)xchar,ngl(i),time
443 ENDIF
444#include "lockoff.inc"
445 ENDDO
446 ENDIF
447 ENDIF
448
449
450 DO i=1,nel
451 dfmax(i)=
min(one,
max(dfmax(i),uvar(i,1)/dcrit))
452 ENDDO
453
454 2000 FORMAT(1x,'FAILURE OF SHELL ELEMENT (TAB)',i10,1x,
455 .'LAYER',i10)
456 2100 FORMAT(1x,'FAILURE OF SHELL ELEMENT (TAB)',i10,1x,
457 .'LAYER',i10,':',/,'AT TIME :',1pe12.4)
458 2200 FORMAT(1x,'STRESS TENSOR SET TO ZERO IN THE LAYER')
459 2400 FORMAT(1x,1pg20.13,' % OF THICKNESS OF SHELL BROKEN ')
460 2500 FORMAT(1x,' LOWER SKIN -> UPPER SKIN ')
461 2600 FORMAT(1x,' UPPER SKIN -> LOWER SKIN ')
462 3700 FORMAT(1x,'STRESS TENSOR SET TO ZERO, LAYER',i10)
463
464 2410 FORMAT(1x,1pg20.13,' % OF THICKNESS OF SHELL ',i10,' BROKEN ')
465 3800 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT (TAB)',i10)
466 3900 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT (TAB)',i10,
467 . 1x,':',/,' AT TIME :',1pe12.4)
468 4000 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT (TAB) ',i10)
469 4100 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT (TAB) ',i10,
470 . 1x,':',/,' AT TIME :',1pe12.4)
471 4200 FORMAT(1x,'DELETE OF ',a5,' CRACKED PHANTOM ELEMENT'/
472 . 1x,'OF THE ORIGINAL SHELL ELEMENT (TAB) ',
473 . i10)
474 4300 FORMAT(1x,'DELETE OF ',a5,' CRACKED PHANTOM ELEMENT'/
475 . 1x,'OF THE ORIGINAL SHELL ELEMENT (TAB) ',
476 . i10,':',/1x,'AT TIME :',1pe20.13)
477 4400 FORMAT(1x,'DELETE OF ',a5,' CRACKED PHANTOM ELEMENT'/
478 . 1x,'OF THE ORIGINAL SHELL ELEMENT (TAB) ',
479 . i10,' LAYER',i10)
480 4500 FORMAT(1x,'DELETE OF ',a5,' CRACKED PHANTOM ELEMENT'/
481 . 1x,'OF THE ORIGINAL SHELL ELEMENT (TAB) ',
482 . i10,' LAYER',i10,':',/1x,'AT TIME :',1pe20.13)
483 4600 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT (TAB)',i10,
484 . 1x,'LAYER',i10)
485 4700 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT (TAB)',i10,
486 . 1x,'LAYER',i10,':',/,' AT TIME :',1pe12.4)
487 4800 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT (TAB) ',i10,
488 . 1x,'LAYER',i10)
489 4900 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT (TAB) ',i10,
490 . 1x,'LAYER',i10,':',/,' AT TIME :',1pe12.4)
491 5010 FORMAT(1x,'SHELL ELEMENT FAILURE DUE TO THINNING (TAB)',i10)
492 5020 FORMAT(1x,'SHELL ELEMENT FAILURE DUE TO THINNING (TAB)',i10,
493 . 1x,':',/1x,'AT TIME :',1pe12.4)
494
495
496 RETURN