40
41
42
43
44
45#include "implicit_f.inc"
46
47
48
49#include "units_c.inc"
50#include "comlock.inc"
51
52
53
54
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
83
84
85
86
87 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,IPG,ILAY,IPT,ISMSTR
88 INTEGER ,INTENT(IN) :: LF_DAMMX
89 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
90 my_real ,
INTENT(IN) :: time,timestep
91 my_real ,
DIMENSION(NEL) ,
INTENT(IN) :: off,epsxx,epsyy,epsxy,
92 . epspxx,epspyy,epspxy,aldt
93 my_real ,
DIMENSION(NUPARAM) ,
INTENT(IN) :: uparam
94 TARGET :: uvar
95
96
97
98 INTEGER ,INTENT(OUT) :: DMG_FLAG
99 INTEGER ,DIMENSION(NEL) ,INTENT(INOUT) :: OFFLY,FOFF
100 my_real ,
DIMENSION(NEL) ,
INTENT(INOUT) :: signxx,signyy,signxy
101 my_real ,
DIMENSION(NEL,LF_DAMMX),
INTENT(INOUT) :: dfmax
102 my_real ,
DIMENSION(NEL) ,
INTENT(OUT) :: tdel,dmg_scale
103 my_real ,
DIMENSION(NEL,NUVAR) ,
INTENT(INOUT) :: uvar
104
105
106
107 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
109 EXTERNAL finter
110
111
112
113
114
115
116
117
118
119
120
121 INTEGER I,J,NINDX,FAILI,MODE,XX,XY,YY,IFUNC_SIZ,STRDEF,STRFLAG
122 INTEGER ,DIMENSION(NEL) :: INDX,IPOSP,IADP,ILENP
123 INTEGER ,DIMENSION(NEL,6) :: FMODE
124 my_real odam(nel,6),epsp(nel,6),
125 . sigoxx(nel),sigoyy(nel),sigoxy(nel),escal(nel),dydxe(nel)
126 my_real dydx,frate,et1,et2,et12,ec1,ec2,ec12,et3,ec3,ec23,et23,ec31,
127 . et31,et1m,et2m,et12m,ec1m,ec2m,ec12m,et3m,ec3m,ec23m,et23m,ec31m,
128 . et31m,di,
alpha,damxx,damyy,damxy,odamxx,odamyy,odamxy,theta,c,s,
129 . epspref,epdot,deps,depsm,fact,fscale_siz,ref_siz,epsd1,epsd2,
130 . epst1,epst2,epst_a,epst_b
131 my_real ,
DIMENSION(NEL) :: eps11,eps22,eps12,epsp11,epsp22,epsp12,
132 . epsm1,epsm2,epsm3,epsm4,epsm5,epsm6
133 my_real ,
DIMENSION(:),
POINTER :: fsize
134
135
136
137 dmg_flag = int(uparam(1))
138 alpha =
min(two*pi*uparam(16)*timestep,one)
139 epspref = uparam(17)
140 fscale_siz = uparam(26)
141 ref_siz = uparam(27)
142 strdef = nint(uparam(28))
143 ifunc_siz = ifunc(13)
144
145
146
147
148
149
150
151
152
153
154
155 strflag = 0
156 IF (strdef == 2) THEN
157 IF (ismstr /= 1 .AND. ismstr /= 3 .AND. ismstr /= 11) THEN
158 strflag = 2
159 END IF
160 ELSE IF (strdef == 3) THEN
161 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN
162 strflag = 3
163 END IF
164 END IF
165
166 SELECT CASE (strflag)
167
168 CASE (2)
169 DO i=1,nel
170 IF (off(i) == one ) THEN
171 theta = half*atan(epsxy(i) /
max((epsxx(i)-epsyy(i)),em20))
172 c = cos(theta)
173 s = sin(theta)
174
175 epst_a = half*(epsxx(i)+epsyy(i))
176 epst_b = sqrt( (half*(epsxx(i)-epsyy(i)))**2 + (half*epsxy(i))**2)
177 epst1 = epst_a + epst_b
178 epst2 = epst_a - epst_b
179 epst1 = exp(epst1) - one
180 epst2 = exp(epst2) - one
181 eps11(i) = c*c*epst1 + s*s*epst2
182 eps22(i) = s*s*epst1 + c*c*epst2
183 eps12(i) = s*c*(epst2 - epst1)
184
185 epst_a = half*(epspxx(i)+epspyy(i))
186 epst_b = sqrt( (half*(epspxx(i)-epspyy(i)))**2 + (half*epspxy(i))**2)
187 epsd1 = epst_a + epst_b
188 epsd2 = epst_a - epst_b
189 epsd1 = (one+epst1)*epsd1
190 epsd2 = (one+epst2)*epsd2
191 epsp11(i) = c*c*epsd1 + s*s*epsd2
192 epsp22(i) = s*s*epsd1 + c*c*epsd2
193 epsp12(i) = s*c*(epsd2 - epsd1)
194 END IF
195 ENDDO
196
197 CASE (3)
198 DO i=1,nel
199 IF (off(i) == one ) THEN
200 theta = half*atan(epsxy(i) /
max((epsxx(i)-epsyy(i)),em20))
201 c = cos(theta)
202 s = sin(theta)
203
204 epst_a = half*(epsxx(i)+epsyy(i))
205 epst_b = sqrt( (half*(epsxx(i)-epsyy(i)))**2 + (half*epsxy(i))**2)
206 epst1 = epst_a + epst_b
207 epst2 = epst_a - epst_b
208 epst1 = log(one + epst1)
209 epst2 = log(one + epst2)
210 eps11(i) = c*c*epst1 + s*s*epst2
211 eps22(i) = s*s*epst1 + c*c*epst2
212 eps12(i) = s*c*(epst2 - epst1)
213
214 epst_a = half*(epspxx(i)+epspyy(i))
215 epst_b = sqrt( (half*(epspxx(i)-epspyy(i)))**2 + (half*epspxy(i))**2)
216 epsd1 = epst_a + epst_b
217 epsd2 = epst_a - epst_b
218 epsd1 = (one/exp(epst1))*epsd1
219 epsd2 = (one/exp(epst2))*epsd2
220 epsp11(i) = c*c*epsd1 + s*s*epsd2
221 epsp22(i) = s*s*epsd1 + c*c*epsd2
222 epsp12(i) = s*c*(epsd2 - epsd1)
223 END IF
224 ENDDO
225
226 CASE DEFAULT
227
228 eps11(1:nel) = epsxx(1:nel)
229 eps22(1:nel) = epsyy(1:nel)
230 eps12(1:nel) = half*epsxy(1:nel)
231 epsp11(1:nel) = epspxx(1:nel)
232 epsp22(1:nel) = epspyy(1:nel)
233 epsp12(1:nel) = half*epspxy(1:nel)
234 END SELECT
235
236
237 fsize => uvar(1:nel,13)
238
239 IF (fsize(1)==zero) THEN
240 IF (ifunc_siz > 0) THEN
241 escal(1:nel) = aldt(1:nel) / ref_siz
242 iposp(1:nel) = 0
243 iadp(1:nel) = npf(ifunc_siz) / 2 + 1
244 ilenp(1:nel) = npf(ifunc_siz + 1) / 2 - iadp(1:nel)
245 CALL vinter2(tf,iadp,iposp,ilenp,nel,escal,dydxe,fsize)
246 fsize(1:nel) = fscale_siz * fsize(1:nel)
247 ELSE
248 fsize(1:nel) = one
249 ENDIF
250 ENDIF
251
252 xx=1
253 yy=2
254 xy=3
255
256 DO j=1 ,6
257 DO i=1, nel
258 odam(i,j) =
min(dfmax(i,1+j),one-em3)
259 epsp(i,j) = uvar(i,j)
260 ENDDO
261 ENDDO
262 DO i =1, nel
263 sigoxx(i) = uvar(i,6+xx)
264 sigoyy(i) = uvar(i,6+yy)
265 sigoxy(i) = uvar(i,6+xy)
266 ENDDO
267
268
269 DO i=1,nel
270 epsp(i,xx) =
alpha * abs(epsp11(i)) + (one-
alpha)*epsp(i,xx)
271 epsp(i,yy) =
alpha * abs(epsp22(i)) + (one-
alpha)*epsp(i,yy)
272 epsp(i,xy) =
alpha * abs(epsp12(i)) + (one-
alpha)*epsp(i,xy)
273 ENDDO
274
275
276
277 nindx = 0
278 DO i=1, nel
279 IF (off(i) == one .AND. offly(i) == 1 .AND. foff(i) == 1) THEN
280 et1 = uparam(4)
281 et1m = uparam(5)
282 et2 = uparam(6)
283 et2m = uparam(7)
284 ec1 = uparam(8)
285 ec1m = uparam(9)
286 ec2 = uparam(10)
287 ec2m = uparam(11)
288 et12 = uparam(12)
289 et12m = uparam(13)
290 ec12 = uparam(14)
291 ec12m = uparam(15)
292 et3 = uparam(18)
293 et3m = uparam(19)
294 ec3 = uparam(20)
295 ec3m = uparam(21)
296 et23 = uparam(22)
297 et23m = uparam(23)
298 ec23 = uparam(24)
299 ec23m = uparam(25)
300 et31 = uparam(22)
301 et31m = uparam(23)
302 ec31 = uparam(24)
303 ec31m = uparam(25)
304
305 fmode(i,1:6) = 0
306 faili = 0
307
308 mode = 1
309 epdot = abs(epsp(i,xx))/epspref
310 IF (ifunc(mode) > 0) THEN
311 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
312 ELSE
313 frate = one
314 ENDIF
315 frate = frate*fsize(i)
316 et1m = frate * et1m
317 et1 = frate * et1
318 deps = eps11(i) - et1
319 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
320 depsm = et1m - et1
321 fact = et1m / eps11(i)
322 di = fact * deps / depsm
323 dfmax(i,1+mode) =
max(dfmax(i,1+mode), di)
324 IF (dfmax(i,1+mode) >= one) THEN
325 faili = 1
326 fmode(i,mode) = 1
327 dfmax(i,1+mode) = one
328 epsm1(i) = et1m
329 ENDIF
330 ENDIF
331
332 mode = 2
333 epdot = abs(epsp(i,yy))/epspref
334 IF (ifunc(mode) > 0) THEN
335 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
336 ELSE
337 frate = one
338 ENDIF
339 frate = frate*fsize(i)
340 et2m = frate * et2m
341 et2 = frate * et2
342 deps =
max(eps22(i) - et2, zero)
343 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
344 depsm = (et2m - et2)
345 fact = et2m / eps22(i)
346 di = fact * deps / depsm
347 dfmax(i,1+mode) =
max(dfmax(i,1+mode), di)
348 IF (dfmax(i,1+mode) >= one) THEN
349 faili = 1
350 fmode(i,mode) = 1
351 dfmax(i,1+mode) = one
352 epsm2(i) = et2m
353 ENDIF
354 ENDIF
355
356 mode = 3
357 epdot = abs(epsp(i,xy))/epspref
358 IF (ifunc(mode) > 0) THEN
359 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
360 ELSE
361 frate = one
362 ENDIF
363 frate = frate*fsize(i)
364 et12m = frate * et12m
365 et12 = frate * et12
366 deps =
max(eps12(i) - et12, zero)
367 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
368 depsm = (et12m - et12)
369 fact = et12m / eps12(i)
370 di = fact * deps / depsm
371 dfmax(i,1+mode) =
max(dfmax(i,1+mode), di)
372 IF (dfmax(i,1+mode) >= one) THEN
373 faili = 1
374 fmode(i,mode) = 1
375 dfmax(i,1+mode) = one
376 epsm3(i) = et12m
377 ENDIF
378 ENDIF
379
380 mode = 4
381 epdot = abs(epsp(i,xx))/epspref
382 IF (ifunc(mode) > 0) THEN
383 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
384 ELSE
385 frate = one
386 ENDIF
387 frate = frate*fsize(i)
388 ec1m = frate * ec1m
389 ec1 = frate * ec1
390 deps = -eps11(i) - ec1
391 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
392 depsm = (ec1m - ec1)
393 fact = ec1m / abs(eps11(i))
394 di = fact * deps / depsm
395 dfmax(i,1+mode) =
max(dfmax(i,1+mode), di)
396 IF (dfmax(i,1+mode) >= one) THEN
397 faili = 1
398 fmode(i,mode) = 1
399 dfmax(i,1+mode) = one
400 epsm4(i) = -ec1m
401 ENDIF
402 ENDIF
403
404 mode = 5
405 epdot = abs(epsp(i,yy))/epspref
406 IF (ifunc(mode) > 0) THEN
407 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
408 ELSE
409 frate = one
410 ENDIF
411 frate = frate*fsize(i)
412 ec2m = frate * ec2m
413 ec2 = frate * ec2
414 deps =
max(-eps22(i) - ec2, zero)
415 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
416 depsm = (ec2m - ec2)
417 fact = ec2m / abs(eps22(i))
418 di = fact * deps / depsm
419 dfmax(i,1+mode) =
max(dfmax(i,1+mode), di)
420 IF (dfmax(i,1+mode) >= one) THEN
421 faili = 1
422 fmode(i,mode) = 1
423 dfmax(i,1+mode) = one
424 epsm5(i) = -ec2m
425 ENDIF
426 ENDIF
427
428 mode = 6
429 epdot = abs(epsp(i,xy))/epspref
430 IF (ifunc(mode) > 0) THEN
431 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
432 ELSE
433 frate = one
434 ENDIF
435 frate = frate*fsize(i)
436 ec12m = frate * ec12m
437 ec12 = frate * ec12
438 deps =
max(-eps12(i) - ec12, zero)
439 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
440 depsm = (ec12m - ec12)
441 fact = ec12m / abs(eps12(i))
442 di = fact * deps / depsm
443 dfmax(i,1+mode) =
max(dfmax(i,1+mode), di)
444 IF (dfmax(i,1+mode) >= one) THEN
445 faili = 1
446 fmode(i,mode) = 1
447 dfmax(i,1+mode) = one
448 epsm6(i) = -ec12m
449 ENDIF
450 ENDIF
451
452 IF (faili == 1) THEN
453 foff(i) = 0
454 tdel(i) = time
455 nindx = nindx + 1
456 indx(nindx) = i
457 ENDIF
458 ENDIF
459 ENDDO
460
461 IF (nindx > 0) THEN
462 DO j=1,nindx
463 i = indx(j)
464#include "lockon.inc"
465 WRITE(iout, 1000) ngl(i),ipg,ilay,ipt
466 WRITE(istdo,1000) ngl(i),ipg,ilay,ipt
467 IF (fmode(i,1)==1) WRITE(iout, 2000) '1 - TRACTION XX',eps11(i),epsm1(i),epsp(i,xx)
468 IF (fmode(i,2)==1) WRITE(iout, 2000) '2 - TRACTION YY',eps22(i),epsm2(i),epsp(i,yy)
469 IF (fmode(i,3)==1) WRITE(iout, 2000) '3 - POSITIVE SHEAR XY',eps12(i),epsm3(i),epsp(i,xy)
470 IF (fmode(i,4)==1) WRITE(iout, 2000) '4 - COMPRESSION XX',eps11(i),epsm4(i),epsp(i,xx)
471 IF (fmode(i,5)==1) WRITE(iout, 2000) '5 - COMPRESSION YY',eps22(i),epsm5(i),epsp(i,yy)
472 IF (fmode(i,6)==1) WRITE(iout, 2000) '6 - NEGATIVE SHEAR XY',eps12(i),epsm6(i),epsp(i,xy)
473#include "lockoff.inc"
474 END DO
475 END IF
476
477
478
479
480 IF (dmg_flag == 1) THEN
481 DO i=1, nel
482 damxx =
max(dfmax(i,2),dfmax(i,5))
483 damyy =
max(dfmax(i,3),dfmax(i,6))
484 damxy =
max(dfmax(i,4),dfmax(i,7))
485 dfmax(i,1) =
max(dfmax(i,1), damxx)
486 dfmax(i,1) =
max(dfmax(i,1), damyy)
487 dfmax(i,1) =
max(dfmax(i,1), damxy)
488 dfmax(i,1) =
min(one,dfmax(i,1))
489 dmg_scale(i) = one - dfmax(i,1)
490 ENDDO
491 ELSE
492 DO i=1, nel
493 damxx =
max(dfmax(i,2),dfmax(i,5))
494 damyy =
max(dfmax(i,3),dfmax(i,6))
495 damxy =
max(dfmax(i,4),dfmax(i,7))
496 dfmax(i,1) =
max(dfmax(i,1), damxx)
497 dfmax(i,1) =
max(dfmax(i,1), damyy)
498 dfmax(i,1) =
max(dfmax(i,1), damxy)
499 dfmax(i,1) =
min(one,dfmax(i,1))
500 odamxx =
max(odam(i,1),odam(i,4))
501 odamyy =
max(odam(i,2),odam(i,5))
502 odamxy =
max(odam(i,3),odam(i,6))
503
504 signxx(i) = sigoxx(i)/(one-odamxx) + (signxx(i) - sigoxx(i))
505 signyy(i) = sigoyy(i)/(one-odamyy) + (signyy(i) - sigoyy(i))
506 signxy(i) = sigoxy(i)/(one-odamxy) + (signxy(i) - sigoxy(i))
507
508 signxx(i) = signxx(i)* (one-damxx)
509 signyy(i) = signyy(i)* (one-damyy)
510 signxy(i) = signxy(i)* (one-damxy)
511 uvar(i,7) = signxx(i)
512 uvar(i,8) = signyy(i)
513 uvar(i,9) = signxy(i)
514 ENDDO
515 ENDIF
516
517 DO j=1 ,6
518 DO i=1, nel
519 uvar(i,j) = epsp(i,j)
520 ENDDO
521 ENDDO
522
523 1000 FORMAT(1x,'FAILURE (ORTHSTR) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
524 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3)
525 2000 FORMAT(1x,'MODE',1x,a,', STRAIN',1pe12.4,', CRITICAL VALUE',1pe12.4,
526 . ', STRAIN RATE',1pe12.4)
527
528 RETURN
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)