OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_hashin_s.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| fail_hashin_s ../engine/source/materials/fail/hashin/fail_hashin_s.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
28!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
29!||====================================================================
30 SUBROUTINE fail_hashin_s(
31 1 NEL ,NUVAR ,ILAY ,NPT0 ,
32 2 TIME ,TIMESTEP,UPARAM ,
33 3 NGL ,OFF ,NOFF ,SIGNXX ,
34 4 SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
35 5 UVAR ,NUPARAM ,DFMAX ,TDELE ,EPSP ,
36 6 LF_DAMMX)
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
875 END
alpha
#define alpha
Definition eval.h:35
fail_hashin_s
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)
Definition fail_hashin_s.F:37
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21