OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_hashin_s.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "units_c.inc"
#include "comlock.inc"
#include "param_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine fail_hashin_s (nel, nuvar, ilay, npt0, time, timestep, uparam, ngl, off, noff, signxx, signyy, signzz, signxy, signyz, signzx, uvar, nuparam, dfmax, tdele, epsp, lf_dammx)

Function/Subroutine Documentation

◆ fail_hashin_s()

subroutine fail_hashin_s ( integer nel,
integer nuvar,
integer ilay,
integer npt0,
time,
timestep,
uparam,
integer, dimension(nel) ngl,
off,
integer, dimension(nel) noff,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
uvar,
integer nuparam,
intent(inout) dfmax,
tdele,
epsp,
integer, intent(in) lf_dammx )

Definition at line 30 of file fail_hashin_s.F.

37C-----------------------------------------------
38C Hashin model ------
39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C---------+---------+---+---+--------------------------------------------
44C VAR | SIZE |TYP| RW| DEFINITION
45C---------+---------+---+---+--------------------------------------------
46C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
47C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
48C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
49C---------+---------+---+---+--------------------------------------------
50C TIME | 1 | F | R | CURRENT TIME
51C TIMESTEP| 1 | F | R | CURRENT TIME STEP
52C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
53C EPSPXX | NEL | F | R | STRAIN RATE XX
54C EPSPYY | NEL | F | R | STRAIN RATE YY
55C ... | | | |
56C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
57C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
58C ... | | | |
59C EPSXX | NEL | F | R | STRAIN XX
60C EPSYY | NEL | F | R | STRAIN YY
61C ... | | | |
62C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
63C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
64C ... | | | |
65C---------+---------+---+---+--------------------------------------------
66C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
67C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
68C ... | | | |
69C SIGVXX | NEL | F | W | VISCOUS STRESS XX
70C SIGVYY | NEL | F | W | VISCOUS STRESS YY
71C ... | | | |
72C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
73C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
74C---------+---------+---+---+--------------------------------------------
75C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
76C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
77C---------+---------+---+---+--------------------------------------------
78#include "mvsiz_p.inc"
79#include "com01_c.inc"
80#include "units_c.inc"
81#include "comlock.inc"
82#include "param_c.inc"
83C-----------------------------------------------
84C I N P U T A r g u m e n t s
85C-----------------------------------------------
86 INTEGER NEL,NUPARAM,ILAY,NPT0,NUVAR
87 INTEGER NGL(NEL)
88 INTEGER, INTENT(IN) :: LF_DAMMX
89 my_real
90 . time,timestep(nel),uparam(nuparam),
91 . signxx(nel),signyy(nel),signzz(nel),
92 . signxy(nel),signyz(nel),signzx(nel),
93 . epsp(nel)
94C-----------------------------------------------
95C O U T P U T A r g u m e n t s
96C-----------------------------------------------
97cc my_real
98
99C-----------------------------------------------
100C I N P U T O U T P U T A r g u m e n t s
101C-----------------------------------------------
102 INTEGER NOFF(NEL)
103 my_real uvar(nel,nuvar), off(nel), tdele(nel)
104 my_real ,DIMENSION(NEL,LF_DAMMX) ,INTENT(INOUT) :: dfmax
105C-----------------------------------------------
106C L o c a l V a r i a b l e s
107C-----------------------------------------------
108 integer
109 . i,j,idel,idel_l,iflag,indx(mvsiz),iadbuf,nindx,
110 . nindex,index(mvsiz),ifail,jj,kk,ir,imatly,
111 . imodel,iunidir,ifabric,indx0(mvsiz),nindx0,
112 . mode(mvsiz,7),imode,failnpt,tmod
113 my_real
114 . sigt1,sigt2,sigt3,
115 . sigc1,sigc2,csig,fsig12,
116 . msig12(mvsiz),msig23(mvsiz),msig13(mvsiz),angle,
117 . sdel,tmax,ratio,fac,s12,s23,
118 . f1,f2,f3,f4,f5,sig,p,f6,f7,epspref,filt,alpha,fs(nel),
119 . telem,k2m,instf,k0,k1,k2,k3,zep669741,dtinv,km,
120 . fsig12up
121C--------------------------------------------------------------
122CC Constant
123 zep669741 = six*em01 + six*em02 + nine*em3 + seven*em04 + four*em05 + em06
124C
125 imodel = int(uparam(1))
126 sigt1 = uparam(2)
127 sigt2 = uparam(3)
128 sigt3 = uparam(4)
129 sigc1 = uparam(5)
130 sigc2 = uparam(6)
131 csig = uparam(7)
132 fsig12 = uparam(8)
133 msig12(1:mvsiz) = uparam(9)
134 msig13(1:mvsiz) = uparam(10)
135 msig23(1:mvsiz) = uparam(11)
136 angle = uparam(12)
137 sdel = uparam(13)
138 tmax = uparam(14)
139 iflag = int(uparam(16))
140 ratio = uparam(17)
141 tmod = int(uparam(19))
142 epspref = uparam(20)
143 filt = uparam(21)
144 km = uparam(22)
145 nindx0 = 0
146 IF (iflag == 3) ratio = one - one/ npt0
147 DO i=1,nel
148 indx(i) = 0
149 indx0(i) = 0
150 ENDDO
151C-----------------------------------------------
152C USER VARIABLES INITIALIZATION
153C-----------------------------------------------
154 IF(isigi == 0) THEN
155 IF(time == zero)THEN
156 DO i=1,nel
157 uvar(i,8) = one
158 ENDDO
159 ENDIF
160 ENDIF
161C-----------------------------------------------
162 iunidir = 0
163 ifabric = 0
164 IF(imodel == 1)THEN
165 iunidir = 1
166 ELSEIF(imodel == 2)THEN
167 ifabric = 1
168 ENDIF
169c ELSEIF(IFLAG(I) == 3)THEN
170c
171 DO i=1,nel
172 IF(off(i) < em01) off(i)=zero
173 IF(off(i) < one) off(i)=off(i)*four_over_5
174C
175 mode(i,1) = 0
176 mode(i,2) = 0
177 mode(i,3) = 0
178 mode(i,4) = 0
179 mode(i,5) = 0
180 mode(i,6) = 0
181 mode(i,7) = 0
182 END DO
183C
184C initialisation
185
186C-------------------------------
187C
188C OFF = 0. si la matrice ou la fibre a rompu
189C-------------------------------
190C
191 IF(iunidir > 0)THEN
192 nindx=0
193 indx = 0
194 indx0=0
195 nindx0 = 0
196c
197 DO i=1,nel
198 f1 =zero
199 f2 =zero
200 f3 =zero
201 f4 = zero
202 f5 = zero
203c alpha = number of cycles
204 dtinv = timestep(i)/max(timestep(i)**2,em20)
205 alpha = int(max(one,filt*dtinv))
206c alpha coef de moyenne mobile exponentielle
207 alpha = two/(alpha + one)
208 fs(i)=uvar(i,10)
209c smoothed strain rate
210 epsp(i)=(one - alpha)*uvar(i,11) + alpha*epsp(i)
211 uvar(i,11)=epsp(i)
212 IF(off(i) == one .AND. imodel == 1)THEN
213C
214 IF(iflag == 1) THEN
215C-------------------------------
216C OFF = 0 when one layer fiber or matrix criteria is reached
217C-------------------------------
218 IF (uvar(i,8) < one)THEN
219 IF(tmod == 1) THEN
220 telem = max((timestep(i)*1.0),tmax*(epspref/epsp(i))**2)
221 uvar(i,8)=uvar(i,8)*exp(-timestep(i)/telem)
222 ELSEIF(tmod == 2) THEN
223 telem = uvar(i,12)
224 uvar(i,8) = exp(-(time - uvar(i,7))/telem)
225 ELSE
226 uvar(i,8)= exp(-(time - uvar(i,7))/tmax)
227 ENDIF
228 IF (uvar(i,8) < em02) uvar(i,8) = zero
229 signxx(i) = uvar(i,1)*uvar(i,8)
230 signyy(i) = uvar(i,2)*uvar(i,8)
231 signzz(i) = uvar(i,3)*uvar(i,8)
232 signxy(i) = uvar(i,4)*uvar(i,8)
233 signyz(i) = uvar(i,5)*uvar(i,8)
234 signzx(i) = uvar(i,6)*uvar(i,8)
235 IF (uvar(i,8) == zero ) THEN
236 off(i)=four_over_5
237 nindx=nindx+1
238 indx(nindx)=i
239 tdele(i) = time
240C KK = 4711
241 ENDIF
242 ELSE
243C
244C fiber criteria
245C
246 sig = half*(signxx(i) + abs(signxx(i)))
247 f1=(sig/sigt1)**2
248 . + (signxy(i)**2 + signzx(i)**2)/fsig12**2
249 dfmax(i,2) = max(dfmax(i,2),f1)
250 dfmax(i,2) = min(dfmax(i,2),one)
251 sig = -half*(signyy(i) + signzz(i))
252 sig = - signxx(i) + half*(sig + abs(sig))
253 sig = half*(sig + abs(sig))
254 f2 = (sig /sigc1)**2
255 dfmax(i,3) = max(dfmax(i,3),f2)
256 dfmax(i,3) = min(dfmax(i,3),one)
257 p = -third*(signxx(i) + signyy(i) + signzz(i))
258 IF(p > 0) f3 = (p/csig)**2
259 dfmax(i,4) = max(dfmax(i,4),f3)
260 dfmax(i,4) = min(dfmax(i,4),one)
261C
262C matrice criteria
263C
264 fac=tan(angle)
265 s12 = msig12(i)
266 s23 = msig23(i)
267 IF(signyy(i) < zero) THEN
268 s12 = msig12(i) - signyy(i)*fac
269 s23 = msig23(i) - signyy(i)*fac
270 ENDIF
271C
272 sig = half*(signyy(i) + abs(signyy(i)))
273 f4=(sig/sigt2)**2
274 . + (signyz(i)/s23)**2 + (signxy(i)/s12)**2
275 dfmax(i,5) = max(dfmax(i,5),f4)
276 dfmax(i,5) = min(dfmax(i,5),one)
277C
278 IF(signzz(i) < zero) THEN
279 msig13(i) = msig13(i) - signzz(i)*fac
280 msig23(i) = msig23(i) - signzz(i)*fac
281 ENDIF
282 sig = half*(signzz(i) + abs(signzz(i)))
283
284 f5= (sig/sigt3)**2 +
285 . (signyz(i)/msig23(i))**2 + (signzx(i)/msig13(i))**2
286 f5 = sdel*sdel*f5
287 dfmax(i,6) = max(dfmax(i,6),f5)
288 dfmax(i,6) = min(dfmax(i,6),one)
289C
290 dfmax(i,1) = min(one,max(dfmax(i,1),f1,f2,f3,f4,f5))
291 uvar(i,10) = max(f1,f2,f3,f4,f5)
292C
293 IF(f1 >= one .OR. f2 >= one .OR. f3 >= one .OR.
294 . f4 >= one .OR. f5 >= one )THEN
295 uvar(i,1) = signxx(i)
296 uvar(i,2) = signyy(i)
297 uvar(i,3) = signzz(i)
298 uvar(i,4) = signxy(i)
299 uvar(i,5) = signyz(i)
300 uvar(i,6) = signzx(i)
301 uvar(i,7) = time
302 IF(tmod == 0) THEN
303 uvar(i,8) = four_over_5
304 ELSE
305 uvar(i,8) = zep99
306 ENDIF
307 fac = hundred*max(abs(signxx(i)),abs(signyy(i)),
308 . abs(signzz(i)),abs(signxy(i)),abs(signyz(i)),
309 . abs(signzx(i)))
310 k0 = fac*epspref**2*zep669741*tmax**2 + zep9*epspref*tmax
311 k1 = zep9*epsp(i)
312 k2 = zep669741*fac*epsp(i)**2
313 k2m = zep669741*fac*km*epspref**2
314 k2 = max(k2m,k2,em20)
315 telem = (sqrt(k1**2+4*k2*k0)-k1)/2/k2
316 uvar(i,12) = telem
317 IF (f1>=one) THEN
318 mode(i,1) = 1
319 kk = 1
320 ENDIF
321 IF (f2>=one) THEN
322 mode(i,2) = 2
323 kk = 2
324 ENDIF
325 IF (f3>=one) THEN
326 mode(i,3) = 3
327 kk = 3
328 ENDIF
329 IF (f4>=one) THEN
330 mode(i,4) = 4
331 kk = 4
332 ENDIF
333 IF (f5>=one) THEN
334 mode(i,5) = 5
335 kk = 5
336 ENDIF
337 uvar(i,9) = kk
338 nindx0=nindx0+1
339 indx0(nindx0)=i
340 ENDIF
341 ENDIF
342C
343 ELSE ! iflag/= 1
344C-------------------------------
345C OFF = 0 when all layer fiber or matrix criteria is reached
346C-------------------------------
347 IF (uvar(i,8) == zero )THEN
348 signxx(i) = zero
349 signyy(i) = zero
350 signzz(i) = zero
351 signxy(i) = zero
352 signzx(i) = zero
353 signyz(i) = zero
354 ELSEIF (uvar(i,8) < one) THEN
355 IF(tmod == 1) THEN
356 telem = max((timestep(i)*1.0),tmax*(epspref/epsp(i))**2)
357 uvar(i,8)=uvar(i,8)*exp(-timestep(i)/telem)
358 ELSEIF(tmod == 2) THEN
359 telem = uvar(i,12)
360 uvar(i,8)= exp(-(time - uvar(i,7))/telem )
361 ELSE
362 uvar(i,8)= exp(-(time - uvar(i,7))/tmax)
363 ENDIF
364 IF(uvar(i,8) < em02)uvar(i,8) = zero
365 signxx(i) = uvar(i,1)*uvar(i,8)
366 signyy(i) = uvar(i,2)*uvar(i,8)
367 signzz(i) = uvar(i,3)*uvar(i,8)
368 signxy(i) = uvar(i,4)*uvar(i,8)
369 signyz(i) = uvar(i,5)*uvar(i,8)
370 signzx(i) = uvar(i,6)*uvar(i,8)
371 IF (uvar(i,8) == zero )THEN
372 noff(i) = noff(i) + 1
373 failnpt = int(npt0*ratio)
374 IF (noff(i) >= failnpt .OR. npt0 == 1) THEN
375 nindx=nindx+1
376 indx(nindx)=i
377 off(i) = four_over_5
378 tdele(i) = time
379 ENDIF
380 ENDIF
381 ELSE
382C
383C fiber criteria
384C
385 sig = half*(signxx(i) + abs(signxx(i)))
386 f1=(sig/sigt1)**2
387 . + (signxy(i)**2 + signzx(i)**2)/fsig12**2
388 dfmax(i,2) = max(dfmax(i,2),f1)
389 dfmax(i,2) = min(dfmax(i,2),one)
390 sig = -half*(signyy(i) + signzz(i))
391 sig = - signxx(i) + half*(sig + abs(sig))
392 sig = half*(sig + abs(sig))
393 f2 = (sig /sigc1)**2
394 dfmax(i,3) = max(dfmax(i,3),f2)
395 dfmax(i,3) = min(dfmax(i,3),one)
396 p = -third*(signxx(i) + signyy(i) + signzz(i))
397 IF(p > 0) f3 = (p/csig)**2
398 dfmax(i,4) = max(dfmax(i,4),f3)
399 dfmax(i,4) = min(dfmax(i,4),one)
400C
401C matrice criteria
402C
403 fac=tan(angle)
404 s12 = msig12(i)
405 s23 = msig23(i)
406 IF(signyy(i) < zero) THEN
407 s12 = msig12(i) - signyy(i)*fac
408 s23 = msig23(i) - signyy(i)*fac
409 ENDIF
410C
411 sig = half*(signyy(i) + abs(signyy(i)))
412 f4=(sig/sigt2)**2
413 . + (signyz(i)/s23)**2 + (signxy(i)/s12)**2
414 dfmax(i,5) = max(dfmax(i,5),f4)
415 dfmax(i,5) = min(dfmax(i,5),one)
416C
417 IF(signzz(i) < zero) THEN
418 msig13(i) = msig13(i) - signzz(i)*fac
419 msig23(i) = msig23(i) - signzz(i)*fac
420 ENDIF
421 sig = half*(signzz(i) + abs(signzz(i)))
422 f5= (sig/sigt3)**2 +
423 . (signyz(i)/msig23(i))**2 + (signzx(i)/msig13(i))**2
424 f5 = sdel*sdel*f5
425 dfmax(i,6) = max(dfmax(i,6),f5)
426 dfmax(i,6) = min(dfmax(i,6),one)
427C
428 dfmax(i,1) = min(one,max(dfmax(i,1),f1,f2,f3,f4,f5))
429C
430 uvar(i,10) = max(f1,f2,f3,f4,f5)
431 IF(f1 >= one .OR. f2 >= one .OR. f3 >= one .OR.
432 . f4 >= one .OR. f5 >= one )THEN
433 uvar(i,1) = signxx(i)
434 uvar(i,2) = signyy(i)
435 uvar(i,3) = signzz(i)
436 uvar(i,4) = signxy(i)
437 uvar(i,5) = signyz(i)
438 uvar(i,6) = signzx(i)
439 uvar(i,7) = time
440 IF(tmod == 0) THEN
441 uvar(i,8) = four_over_5
442 ELSE
443 uvar(i,8) = zep99
444 ENDIF
445 fac = hundred*max(abs(signxx(i)),abs(signyy(i)),
446 . abs(signzz(i)),abs(signxy(i)),abs(signyz(i)),
447 . abs(signzx(i)))
448 k0 = fac*epspref**2*zep669741*tmax**2 + zep9*epspref*tmax
449 k1 = zep9*epsp(i)
450 k2 = zep669741*fac*epsp(i)**2
451 k2m = zep669741*fac*km*epspref**2
452 k2 = max(k2m,k2,em20)
453 telem = (sqrt(k1**2+4*k2*k0)-k1)/2/k2
454 uvar(i,12) = telem
455 nindx0= nindx0+1
456 indx0(nindx0)=i
457 IF (f1>=one) THEN
458 mode(i,1) = 1
459 kk =1
460 ENDIF
461 IF (f2>=one) THEN
462 mode(i,2) = 1
463 kk =2
464 ENDIF
465 IF (f3>=one) THEN
466 mode(i,3) = 1
467 kk =3
468 ENDIF
469 IF (f4>=one) THEN
470 mode(i,4) = 1
471 kk =4
472 ENDIF
473 IF (f5>=one) THEN
474 mode(i,5) = 1
475 kk =5
476 ENDIF
477 uvar(i,9) = kk
478 ENDIF
479 ENDIF
480c
481 ENDIF ! iflag choice
482 ENDIF
483
484 ENDDO
485 ENDIF
486C ----------------------------------------------
487C fabric layer model
488C OFF = 0. si la matrice ou la fibre a rompu
489C-------------------------------
490 IF(ifabric > 0)THEN
491 indx = 0
492 nindx = 0
493 indx0 = 0
494 nindx0 = 0
495 DO i=1,nel
496 f1 = zero
497 f2 = zero
498 f3 = zero
499 f4 = zero
500 f5 = zero
501 f6 = zero
502 f7 = zero
503c alpha = number of cycles
504 dtinv = timestep(i)/max(timestep(i)**2,em20)
505 alpha = int(max(one,filt*dtinv))
506c alpha coef de moyenne mobile exponentielle
507 alpha = two/(alpha + one)
508 fs(i)=uvar(i,10)
509c smoothed strain rate
510 epsp(i)=(one - alpha)*uvar(i,11) + alpha*epsp(i)
511 uvar(i,11)=epsp(i)
512 IF(off(i) == one .AND. imodel /= 1)THEN
513 IF(iflag == 1) THEN
514 IF(uvar(i,8) < one)THEN
515 IF(tmod == 1) THEN
516 telem = max((timestep(i)*1.0),tmax*(epspref/epsp(i))**2)
517 uvar(i,8)=uvar(i,8)*exp(-timestep(i)/telem)
518 ELSEIF(tmod == 2) THEN
519 telem = uvar(i,12)
520 uvar(i,8)= exp(-(time - uvar(i,7))/telem)
521 ELSE
522 uvar(i,8)= exp(-(time - uvar(i,7))/tmax)
523 ENDIF
524 IF(uvar(i,8) < em02)uvar(i,8) = zero
525 signxx(i) = uvar(i,1)*uvar(i,8)
526 signyy(i) = uvar(i,2)*uvar(i,8)
527 signzz(i) = uvar(i,3)*uvar(i,8)
528 signxy(i) = uvar(i,4)*uvar(i,8)
529 signyz(i) = uvar(i,5)*uvar(i,8)
530 signzx(i) = uvar(i,6)*uvar(i,8)
531 IF(uvar(i,8) == zero )THEN
532 off(i)=four_over_5
533 nindx=nindx+1
534 indx(nindx)=i
535 tdele(i) = time
536 ENDIF
537 ELSE
538C
539C Fill and warp directions
540C
541 sig = half*(signxx(i) + abs(signxx(i)))
542 f1=(sig/sigt1)**2
543 . + (signxy(i)**2 + signzx(i)**2)/fsig12**2
544 dfmax(i,2) = max(dfmax(i,2),f1)
545 dfmax(i,2) = min(dfmax(i,2),one)
546C
547 sig = half*(signyy(i) + abs(signyy(i)))
548 fsig12up = fsig12*sigt2/sigt1
549 f2=(sig/sigt2)**2
550 . + (signxy(i)**2 + signyz(i)**2)/fsig12up**2
551 dfmax(i,3) = max(dfmax(i,3),f2)
552 dfmax(i,3) = min(dfmax(i,3),one)
553C
554 sig = half*(-signzz(i) + abs(signzz(i)))
555 IF(-signxx(i) + sig > zero)
556 . f3 = ((-signxx(i) + sig)/sigc1)**2
557 dfmax(i,4) = max(dfmax(i,4),f3)
558 dfmax(i,4) = min(dfmax(i,4),one)
559C
560 IF(-signyy(i) + sig > zero)
561 . f4 = ((-signyy(i) + sig)/sigc2)**2
562 dfmax(i,5) = max(dfmax(i,5),f4)
563 dfmax(i,5) = min(dfmax(i,5),one)
564C
565 p = -third*(signxx(i) + signyy(i) + signzz(i))
566 IF(p > zero) f5 = (p/csig)**2
567 dfmax(i,6) = max(dfmax(i,6),f5)
568 dfmax(i,6) = min(dfmax(i,6),one)
569C F5 = (SIGNXY(I)/MSIG12(I))**2
570 f6 = (signxy(i)/msig12(i))**2
571 dfmax(i,7) = max(dfmax(i,7),f6)
572 dfmax(i,7) = min(dfmax(i,7),one)
573C
574 sig = half*(signzz(i) + abs(signzz(i)))
575 IF(signzz(i) < zero) THEN
576 fac=tan(angle)
577 msig13(i) = msig13(i) - signzz(i)*fac
578 msig23(i) = msig23(i) - signzz(i)*fac
579 ENDIF
580C F6 = (SIG/SIGT3(I))**2 +
581C . (SIGNYZ(I)/MSIG23(I))**2 + (SIGNZX(I)/MSIG13(I))**2
582CC
583C F6 = SDEL(I)*SDEL(I)*F6
584 f7 = (sig/sigt3)**2 +
585 . (signyz(i)/msig23(i))**2 + (signzx(i)/msig13(i))**2
586C
587 f7 = sdel*sdel*f7
588 dfmax(i,8) = max(dfmax(i,8),f7)
589 dfmax(i,8) = min(dfmax(i,8),one)
590C
591 dfmax(i,1) = min(one,max(dfmax(i,1),f1,f2,f3,f4,f5,f6,f7))
592C
593 uvar(i,10) = max(f1,f2,f3,f4,f5,f6,f7)
594 IF(f1 >= one .OR. f2 >= one .OR. f3 >= one .OR.
595 . f4 >= one .OR. f5 >= one .OR. f6 >= one .OR.
596 . f7 >= one )THEN
597 uvar(i,1) = signxx(i)
598 uvar(i,2) = signyy(i)
599 uvar(i,3) = signzz(i)
600 uvar(i,4) = signxy(i)
601 uvar(i,5) = signyz(i)
602 uvar(i,6) = signzx(i)
603 uvar(i,7) = time
604 IF(tmod == 0) THEN
605 uvar(i,8) = four_over_5
606 ELSE
607 uvar(i,8) = zep99
608 ENDIF
609 fac = hundred*max(abs(signxx(i)),abs(signyy(i)),
610 . abs(signzz(i)),abs(signxy(i)),abs(signyz(i)),
611 . abs(signzx(i)))
612 k0 = fac*epspref**2*zep669741*tmax**2 + zep9*epspref*tmax
613 k1 = zep9*epsp(i)
614 k2 = zep669741*fac*epsp(i)**2
615 k2m = zep669741*fac*km*epspref**2
616 k2 = max(k2m,k2,em20)
617 telem = (sqrt(k1**2+4*k2*k0)-k1)/2/k2
618 uvar(i,12) = telem
619 IF (f1>=one) THEN
620 mode(i,1) = 1
621 kk = 1
622 ENDIF
623
624 IF (f2>=one) THEN
625 mode(i,2) = 1
626 kk = 2
627 ENDIF
628 IF (f3>=one) THEN
629 mode(i,3) = 1
630 kk = 3
631 ENDIF
632 IF (f4>=one) THEN
633 mode(i,4) = 1
634 kk = 4
635 ENDIF
636 IF (f5>=one) THEN
637 mode(i,5) = 1
638 kk = 5
639 ENDIF
640 IF (f6>=one) THEN
641 mode(i,6) = 1
642 kk = 5
643 ENDIF
644 IF (f7>=one) THEN
645 mode(i,7) = 1
646 kk = 7
647 ENDIF
648 uvar(i,9) = kk
649C
650 nindx0=nindx0+1
651 indx0(nindx0)=i
652 ENDIF
653 ENDIF
654 ELSE ! iflag/= 1
655C-------------------------------
656C OFF = 0. all layer fiber or matrix criteria is reateched
657C-------------------------------
658 IF(uvar(i,8) == zero) THEN
659 signxx(i) = zero
660 signyy(i) = zero
661 signzz(i) = zero
662 signxy(i) = zero
663 signyz(i) = zero
664 signzx(i) = zero
665 ELSEIF(uvar(i,8) < one)THEN
666 IF(tmod == 1) THEN
667 telem = max((timestep(i)*1.0),tmax*(epspref/epsp(i))**2)
668 uvar(i,8)=uvar(i,8)*exp(-timestep(i)/telem)
669 ELSEIF(tmod == 2) THEN
670 telem = uvar(i,12)
671 uvar(i,8)=exp(-(time - uvar(i,7))/telem)
672 ELSE
673 uvar(i,8)= exp(-(time - uvar(i,7))/tmax)
674 ENDIF
675 IF(uvar(i,8) < em02)uvar(i,8) = zero
676 signxx(i) = uvar(i,1)*uvar(i,8)
677 signyy(i) = uvar(i,2)*uvar(i,8)
678 signzz(i) = uvar(i,3)*uvar(i,8)
679 signxy(i) = uvar(i,4)*uvar(i,8)
680 signyz(i) = uvar(i,5)*uvar(i,8)
681 signzx(i) = uvar(i,6)*uvar(i,8)
682 IF(uvar(i,8) == zero )THEN
683 noff(i) = noff(i) + 1
684 failnpt = int((npt0*ratio))
685 IF( npt0 == 1 .OR. noff(i) >= failnpt ) THEN
686 off(i) = four_over_5
687 nindx=nindx+1
688 indx(nindx)=i
689 tdele(i) = time
690 ENDIF
691 ENDIF
692 ELSE
693C
694C Fill and warp directions
695C
696 sig = half*(signxx(i) + abs(signxx(i)))
697 f1=(sig/sigt1)**2
698 . + (signxy(i)**2 + signzx(i)**2)/fsig12**2
699 dfmax(i,2) = max(dfmax(i,2),f1)
700 dfmax(i,2) = min(dfmax(i,2),one)
701
702 sig = half*(signyy(i) + abs(signyy(i)))
703 fsig12up = fsig12*sigt2/sigt1
704 f2=(sig/sigt2)**2
705 . + (signxy(i)**2 + signyz(i)**2)/fsig12up**2
706 dfmax(i,3) = max(dfmax(i,3),f2)
707 dfmax(i,3) = min(dfmax(i,3),one)
708C
709 sig = half*(-signzz(i) + abs(signzz(i)))
710 IF(-signxx(i) + sig > zero)
711 . f3 = ((-signxx(i) + sig)/sigc1)**2
712 dfmax(i,4) = max(dfmax(i,4),f3)
713 dfmax(i,4) = min(dfmax(i,4),one)
714C
715 IF(-signyy(i) + sig > zero)
716 . f4 = ((-signyy(i) + sig)/sigc2)**2
717 dfmax(i,5) = max(dfmax(i,5),f4)
718 dfmax(i,5) = min(dfmax(i,5),one)
719C
720 p = -third*(signxx(i) + signyy(i) + signzz(i))
721 IF(p > zero) f5 = (p/csig)**2
722 dfmax(i,6) = max(dfmax(i,6),f5)
723 dfmax(i,6) = min(dfmax(i,6),one)
724C
725 f6 = (signxy(i)/msig12(i))**2
726 dfmax(i,7) = max(dfmax(i,7),f6)
727 dfmax(i,7) = min(dfmax(i,7),one)
728C
729 sig = half*(signzz(i) + abs(signzz(i)))
730 IF(signzz(i) < zero) THEN
731 fac=tan(angle)
732 msig13(i) = msig13(i) - signzz(i)*fac
733 msig23(i) = msig23(i) - signzz(i)*fac
734 ENDIF
735C F6 = (SIG/SIGT3(I))**2 +
736 f7 = (sig/sigt3)**2 +
737 . (signyz(i)/msig23(i))**2 + (signzx(i)/msig13(i))**2
738C F6 = SDEL(I)*SDEL(I)*F6
739 f7 = sdel*sdel*f7
740 dfmax(i,8) = max(dfmax(i,8),f7)
741 dfmax(i,8) = min(dfmax(i,8),one)
742C
743 dfmax(i,1) = min(one,max(dfmax(i,1),f1,f2,f3,f4,f5,f6,f7))
744 uvar(i,10) = max(f1,f2,f3,f4,f5,f6,f7)
745C
746 IF(f1 >= one .OR. f2 >= one .OR. f3 >= one .OR.
747 . f4 >= one .OR. f5 >= one .OR. f6 >= one .OR.
748 . f7 >= one )THEN
749 uvar(i,1) = signxx(i)
750 uvar(i,2) = signyy(i)
751 uvar(i,3) = signzz(i)
752 uvar(i,4) = signxy(i)
753 uvar(i,5) = signyz(i)
754 uvar(i,6) = signzx(i)
755 uvar(i,7) = time
756 IF(tmod == 0) THEN
757 uvar(i,8) = four_over_5
758 ELSE
759 uvar(i,8) = zep99
760 ENDIF
761C
762 fac = hundred*max(abs(signxx(i)),abs(signyy(i)),
763 . abs(signzz(i)),abs(signxy(i)),abs(signyz(i)),
764 . abs(signzx(i)))
765 k0 = fac*epspref**2*zep669741*tmax**2 + zep9*epspref*tmax
766 k1 = zep9*epsp(i)
767 k2 = zep669741*fac*epsp(i)**2
768 k2m = zep669741*fac*km*epspref**2
769 k2 = max(k2m,k2,em20)
770 telem = (sqrt(k1**2+4*k2*k0)-k1)/2/k2
771 uvar(i,12) = telem
772C
773 nindx0= nindx0+1
774 indx0(nindx0)=i
775 IF (f1>=one) THEN
776 mode(i,1) = 1
777 kk = 1
778 ENDIF
779 IF (f2>=one) THEN
780 mode(i,2) = 1
781 kk = 2
782 ENDIF
783 IF (f3>=one) THEN
784 mode(i,3) = 1
785 kk = 3
786 ENDIF
787 IF (f4>=one) THEN
788 mode(i,4) = 1
789 kk = 4
790 ENDIF
791 IF (f5>=one) THEN
792 mode(i,5) = 1
793 kk = 5
794 ENDIF
795 IF (f6>=one) THEN
796 mode(i,6) = 1
797 kk = 6
798 ENDIF
799 IF (f7>=one) THEN
800 mode(i,7) = 1
801 kk = 7
802 ENDIF
803 uvar(i,9) = kk
804 ENDIF
805 ENDIF
806 ENDIF ! iflag
807 ENDIF
808 ENDDO
809
810 ENDIF
811
812 IF(nindx0 > 0)THEN
813 DO j=1,nindx0
814 i = indx0(j)
815 DO jj=1,7
816 imode = mode(i,jj)
817 IF(imode == 1) THEN
818 imode = jj
819#include "lockon.inc"
820 WRITE(iout, 1100) ilay,ngl(i),imode,time
821 WRITE(iout, 1110) uvar(i,1),uvar(i,2),uvar(i,3)
822 WRITE(iout, 1115) uvar(i,4),uvar(i,6),uvar(i,5)
823 WRITE(istdo,1100) ilay,ngl(i),imode,time
824 WRITE(istdo,1110) uvar(i,1),uvar(i,2),uvar(i,3)
825 WRITE(istdo,1115) uvar(i,4),uvar(i,6),uvar(i,5)
826#include "lockoff.inc"
827 ENDIF
828 END DO
829 END DO
830 ENDIF
831C
832 IF(nindx > 0)THEN
833 DO j=1,nindx
834 i = indx(j)
835C JJ = 0
836C DO JJ=1,7
837C IMODE = MODE(JJ,I)
838#include "lockon.inc"
839C IF(IMODE == 1) THEN
840C IMODE = JJ
841 WRITE(iout, 1200) ngl(i),ilay,int(uvar(i,9)),time
842 WRITE(iout, 1210) uvar(i,1),uvar(i,2),uvar(i,3)
843 WRITE(iout, 1215) uvar(i,4),uvar(i,6),uvar(i,5)
844 WRITE(istdo,1200) ngl(i),ilay,int(uvar(i,9)),time
845 WRITE(istdo,1210) uvar(i,1),uvar(i,2),uvar(i,3)
846 WRITE(istdo,1215) uvar(i,4),uvar(i,6),uvar(i,5)
847#include "lockoff.inc"
848C ENDIF
849C END DO
850 END DO
851 ENDIF
852C--------------------------------------------
853 1000 FORMAT(1x,'FAILURE ELEMENT-1 #',i10,1x,
854 .'LAYER # ',i10,1x,'MODE #',i10)
855C . 1X, 'FAILURE NUMBER LAYER ', I10)
856 1100 FORMAT(1x,'FAILURE LAYER #',i10,1x,
857 .'ELEMENT #',i10,1x,'HASHIN MODE #',i10,1x, 'AT TIME #:',1pe12.4)
858 1110 FORMAT(1x,'RESPONSIBLE STRESS: SIG_11= ',1pe12.4,1x,
859 .'SIG_22= ',1pe12.4,1x,'SIG_33= ',1pe12.4)
860 1115 FORMAT(1x,'RESPONSIBLE STRESS: SIG_12= ',1pe12.4,1x,
861 .'SIG_23= ',1pe12.4,1x,'SIG_13= ',1pe12.4)
862 2000 FORMAT(1x,'FAILURE ELEMENT #',i10,1x,
863 .'LAYER # ',i10,1x,'MODE #',i10)
864C . 1X, 'FAILURE NUMBER LAYER ', I10)
865 2100 FORMAT(1x,'FAILURE ELEMENT-2 #',i10,1x,
866 .'LAYER #',i10,1x, 'AT TIME #:',1pe12.4,1x,'MODE #',i10)
867 1200 FORMAT(1x,'DELETE SOLID ELEMENT # ',i10,1x,'Failed layer # ',i10,
868 . 1x,' HASHIN mode # ',i10,1x,'AT TIME # ',1pe12.4)
869 1210 FORMAT(1x,'RESPONSIBLE STRESS: SIG_11= ',1pe12.4,1x,
870 .'SIG_22= ',1pe12.4,1x,'SIG_33= ',1pe12.4)
871 1215 FORMAT(1x,'RESPONSIBLE STRESS: SIG_12= ',1pe12.4,1x,
872 .'SIG_13= ',1pe12.4,1x,'SIG_23= ',1pe12.4,1x)
873C--------------------------------------------
874 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21