OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_orthstrain_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_orthstrain ../engine/source/materials/fail/orthstrain/fail_orthstrain_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!||--- calls -----------------------------------------------------
30!|| finter ../engine/source/tools/curve/finter.F
31!|| jacobiew ../engine/source/materials/mat/mat033/sigeps33.F
32!|| vinter2 ../engine/source/tools/curve/vinter.F
33!||====================================================================
34 SUBROUTINE fail_orthstrain (
35 1 NEL ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,
36 2 NPF ,TF ,TIME ,TIMESTEP,UPARAM ,ISMSTR ,
37 3 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
38 4 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
39 5 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 6 UVAR ,OFF ,IPT ,NGL ,DFMAX ,TDEL ,
41 7 UELR ,NPT ,DELTAX ,LF_DAMMX)
42C-----------------------------------------------
43c Orthotropic strain failure model
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C-----------------------------------------------
49C G l o b a l P a r a m e t e r s
50C-----------------------------------------------
51#include "units_c.inc"
52#include "comlock.inc"
53C---------+---------+---+---+--------------------------------------------
54C VAR | SIZE |TYP| RW| DEFINITION
55C---------+---------+---+---+--------------------------------------------
56C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
57C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
58C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
59C---------+---------+---+---+--------------------------------------------
60C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
61C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
62C NPF | * | I | R | FUNCTION ARRAY
63C TF | * | F | R | FUNCTION ARRAY
64C---------+---------+---+---+--------------------------------------------
65C TIME | 1 | F | R | CURRENT TIME
66C TIMESTEP| 1 | F | R | CURRENT TIME STEP
67C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
68C EPSPXX | NEL | F | R | STRAIN RATE XX
69C EPSPYY | NEL | F | R | STRAIN RATE YY
70C ... | | | |
71C EPSXX | NEL | F | R | STRAIN XX
72C EPSYY | NEL | F | R | STRAIN YY
73C ... | | | |
74C---------+---------+---+---+--------------------------------------------
75C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
76C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
77C ... | | | |
78C---------+---------+---+---+--------------------------------------------
79C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
80C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
81C---------+---------+---+---+--------------------------------------------
82C I N P U T A r g u m e n t s
83C--------------------------------------------
84 INTEGER :: NEL,NUPARAM,NUVAR,IPT,NPT,ISMSTR,LF_DAMMX
85 INTEGER ,DIMENSION(NEL) :: NGL
86 my_real :: TIME,TIMESTEP
87 my_real ,DIMENSION(NUPARAM) :: UPARAM
88 my_real ,DIMENSION(NEL) ,INTENT(IN) :: EPSXX,EPSYY,EPSZZ,
89 . EPSXY,EPSYZ,EPSZX,EPSPXX,EPSPYY,EPSPZZ,EPSPXY,EPSPYZ,EPSPZX
90 my_real ,DIMENSION(NEL) :: UELR,DELTAX,
91 . signxx,signyy,signzz,signxy,signyz,signzx
92C-----------------------------------------------
93C I N P U T O U T P U T A r g u m e n t s
94C-----------------------------------------------
95 my_real uvar(nel,nuvar),off(nel),tdel(nel)
96 my_real, DIMENSION(NEL,LF_DAMMX), INTENT(INOUT) :: dfmax
97C-----------------------------------------------
98C VARIABLES FOR FUNCTION INTERPOLATION
99C-----------------------------------------------
100 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
101 my_real FINTER ,TF(*)
102 EXTERNAL finter
103C-----------------------------------------------
104C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
105C Y : y = f(x)
106C X : x
107C DYDX : f'(x) = dy/dx
108C IFUNC(J): FUNCTION INDEX
109C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
110C NPF,TF : FUNCTION PARAMETER
111C-----------------------------------------------
112C L o c a l V a r i a b l e s
113C-----------------------------------------------
114 INTEGER I,J,NINDX,FAILI,MODE,XX,XY,YY,ZZ,YZ,ZX,IFUNC_SIZ,STRDEF,STRFLAG,
115 . M11T, M11C,M22T,M22C,M33T,M33C,M12T,M12C,M23T,M23C,M31T,M31C
116 INTEGER ,DIMENSION(NEL) :: INDX,IPOSP,IADP,ILENP
117 INTEGER ,DIMENSION(NEL,12) :: FMODE
118 my_real ODAM(NEL,12),EPSP(NEL,6),SIGO(NEL,6),FSIZE(NEL),
119 . ESCAL(NEL),DYDX(NEL)
120 my_real ,DIMENSION(NEL) :: EPS11,EPS22,EPS33,EPS12,EPS23,EPS31,
121 . epsp11,epsp22,epsp33,epsp12,epsp23,epsp31,
122 . epsm1,epsm2,epsm3,epsm4,epsm5,epsm6,
123 . epsm7,epsm8,epsm9,epsm10,epsm11,epsm12
124 my_real frate,epspref,epdot,alpha,depsm,deps,fact,p,i1,i2,i3,q,r,r_inter,
125 . phi,et1,et2,et12,ec1,ec2,ec12,et3,ec3,ec23,et23,ec31,et31,
126 . et1m,et2m,et12m,ec1m,ec2m,ec12m,et3m,ec3m,ec23m,et23m,ec31m,et31m,
127 . di,damxx,damyy,damxy,damzz,damyz,damzx,
128 . odamxx,odamyy,odamxy,odamzz,odamyz,odamzx,fscale_siz,ref_siz
129 my_real epstens(3,3),epstens2(3,3),eps_pr(3),vect_pr(3,3),
130 . epsptens(3,3),epsptens2(3,3),epsp_pr(3)
131 INTEGER NROT,K,L,M
132C=======================================================================
133 ALPHA = min(two*pi*uparam(16)*max(timestep,em20),one)
134 epspref = uparam(17)
135 fscale_siz = uparam(26)
136 ref_siz = uparam(27)
137 strdef = nint(uparam(28))
138 ifunc_siz = ifunc(13)
139c----------------------------------------------
140c strain transformation following input definition
141c-------------------
142 strflag = 0
143 IF (strdef == 2) THEN ! failure defined as engineering strain
144 IF (ismstr /=1 .AND. ismstr /=3 .AND. ismstr /= 11) THEN
145 strflag = 2
146 END IF
147 ELSE IF (strdef == 3) THEN ! failure defined as true strain
148 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN
149 strflag = 3
150 END IF
151 END IF
152c--------------------------
153 SELECT CASE (strflag)
154c
155 CASE (2) ! transform true strain to engineering
156 DO i=1,nel
157 IF (off(i) == one ) THEN
158 ! -> Fill the strain tensor 3x3
159 epstens(1,1) = epsxx(i)
160 epstens(2,2) = epsyy(i)
161 epstens(3,3) = epszz(i)
162 epstens(1,2) = half*epsxy(i)
163 epstens(2,3) = half*epsyz(i)
164 epstens(1,3) = half*epszx(i)
165 epstens(2,1) = epstens(1,2)
166 epstens(3,1) = epstens(1,3)
167 epstens(3,2) = epstens(2,3)
168 ! -> Compute principal strains and vectors
169 CALL jacobiew(epstens,3,eps_pr,vect_pr,nrot)
170 ! -> Transform true principal strains to engineering principal strains
171 eps_pr(1) = exp(eps_pr(1)) - one
172 eps_pr(2) = exp(eps_pr(2)) - one
173 eps_pr(3) = exp(eps_pr(3)) - one
174 ! -> Recover engineering strain tensor
175 epstens(1:3,1:3) = zero
176 epstens(1,1) = eps_pr(1)
177 epstens(2,2) = eps_pr(2)
178 epstens(3,3) = eps_pr(3)
179 ! --> Principal vector*Principal strains*Transpose(Principal vector)
180 DO k = 1,3
181 DO l = 1,3
182 epstens2(k,l) = zero
183 DO m = 1,3
184 epstens2(k,l) = epstens2(k,l) + epstens(k,m)*vect_pr(l,m)
185 ENDDO
186 ENDDO
187 ENDDO
188 DO k = 1,3
189 DO l = 1,3
190 epstens(k,l) = zero
191 DO m = 1,3
192 epstens(k,l) = epstens(k,l) + vect_pr(k,m)*epstens2(m,l)
193 ENDDO
194 ENDDO
195 ENDDO
196 ! -> Save engineering strains
197 eps11(i) = epstens(1,1)
198 eps22(i) = epstens(2,2)
199 eps33(i) = epstens(3,3)
200 eps12(i) = epstens(1,2)
201 eps23(i) = epstens(2,3)
202 eps31(i) = epstens(3,1)
203 ! -> Fill the strain rate tensor
204 epsptens(1,1) = epspxx(i)
205 epsptens(2,2) = epspyy(i)
206 epsptens(3,3) = epspzz(i)
207 epsptens(1,2) = half*epspxy(i)
208 epsptens(2,3) = half*epspyz(i)
209 epsptens(1,3) = half*epspzx(i)
210 epsptens(2,1) = epsptens(1,2)
211 epsptens(3,1) = epsptens(1,3)
212 epsptens(3,2) = epsptens(2,3)
213 ! -> Compute principal strain rates and vectors
214 CALL jacobiew(epsptens,3,epsp_pr,vect_pr,nrot)
215 ! -> Transform true principal strain rates to engineering principal strain rates
216 epsp_pr(1) = (eps_pr(1)+one)*epsp_pr(1)
217 epsp_pr(2) = (eps_pr(2)+one)*epsp_pr(2)
218 epsp_pr(3) = (eps_pr(3)+one)*epsp_pr(3)
219 ! -> Recover engineering strain rate tensor
220 epsptens(1:3,1:3) = zero
221 epsptens(1,1) = epsp_pr(1)
222 epsptens(2,2) = epsp_pr(2)
223 epsptens(3,3) = epsp_pr(3)
224 ! --> Principal vector*Principal strains*Transpose(Principal vector)
225 DO k = 1,3
226 DO l = 1,3
227 epsptens2(k,l) = zero
228 DO m = 1,3
229 epsptens2(k,l) = epsptens2(k,l) + epsptens(k,m)*vect_pr(l,m)
230 ENDDO
231 ENDDO
232 ENDDO
233 DO k = 1,3
234 DO l = 1,3
235 epsptens(k,l) = zero
236 DO m = 1,3
237 epsptens(k,l) = epsptens(k,l) + vect_pr(k,m)*epsptens2(m,l)
238 ENDDO
239 ENDDO
240 ENDDO
241 ! -> Save engineering strain rates
242 epsp11(i) = epsptens(1,1)
243 epsp22(i) = epsptens(2,2)
244 epsp33(i) = epsptens(3,3)
245 epsp12(i) = epsptens(1,2)
246 epsp23(i) = epsptens(2,3)
247 epsp31(i) = epsptens(3,1)
248 END IF
249 ENDDO
250c
251 CASE (3) ! transform engineering to true strain
252 DO i=1,nel
253 IF (off(i) == one ) THEN
254 ! -> Fill the strain tensor 3x3
255 epstens(1,1) = epsxx(i)
256 epstens(2,2) = epsyy(i)
257 epstens(3,3) = epszz(i)
258 epstens(1,2) = half*epsxy(i)
259 epstens(2,3) = half*epsyz(i)
260 epstens(1,3) = half*epszx(i)
261 epstens(2,1) = epstens(1,2)
262 epstens(3,1) = epstens(1,3)
263 epstens(3,2) = epstens(2,3)
264 ! -> Compute principal strains and vectors
265 CALL jacobiew(epstens,3,eps_pr,vect_pr,nrot)
266 ! -> Transform engineering principal strains to true principal strains
267 eps_pr(1) = log(max(one + eps_pr(1),em20))
268 eps_pr(2) = log(max(one + eps_pr(2),em20))
269 eps_pr(3) = log(max(one + eps_pr(3),em20))
270 ! -> Recover true strain tensor
271 epstens(1:3,1:3) = zero
272 epstens(1,1) = eps_pr(1)
273 epstens(2,2) = eps_pr(2)
274 epstens(3,3) = eps_pr(3)
275 ! --> Principal vector*Principal strains*Transpose(Principal vector)
276 DO k = 1,3
277 DO l = 1,3
278 epstens2(k,l) = zero
279 DO m = 1,3
280 epstens2(k,l) = epstens2(k,l) + epstens(k,m)*vect_pr(l,m)
281 ENDDO
282 ENDDO
283 ENDDO
284 DO k = 1,3
285 DO l = 1,3
286 epstens(k,l) = zero
287 DO m = 1,3
288 epstens(k,l) = epstens(k,l) + vect_pr(k,m)*epstens2(m,l)
289 ENDDO
290 ENDDO
291 ENDDO
292 ! -> Save true strains
293 eps11(i) = epstens(1,1)
294 eps22(i) = epstens(2,2)
295 eps33(i) = epstens(3,3)
296 eps12(i) = epstens(1,2)
297 eps23(i) = epstens(2,3)
298 eps31(i) = epstens(3,1)
299 ! -> Fill the strain rate tensor
300 epsptens(1,1) = epspxx(i)
301 epsptens(2,2) = epspyy(i)
302 epsptens(3,3) = epspzz(i)
303 epsptens(1,2) = half*epspxy(i)
304 epsptens(2,3) = half*epspyz(i)
305 epsptens(1,3) = half*epspzx(i)
306 epsptens(2,1) = epsptens(1,2)
307 epsptens(3,1) = epsptens(1,3)
308 epsptens(3,2) = epsptens(2,3)
309 ! -> Compute principal strain rates and vectors
310 CALL jacobiew(epsptens,3,epsp_pr,vect_pr,nrot)
311 ! -> Transform engineering principal strain rates to true principal strain rates
312 epsp_pr(1) = (one/(exp(eps_pr(1))))*epsp_pr(1)
313 epsp_pr(2) = (one/(exp(eps_pr(2))))*epsp_pr(2)
314 epsp_pr(3) = (one/(exp(eps_pr(3))))*epsp_pr(3)
315 ! -> Recover true strain rates tensor
316 epsptens(1:3,1:3) = zero
317 epsptens(1,1) = epsp_pr(1)
318 epsptens(2,2) = epsp_pr(2)
319 epsptens(3,3) = epsp_pr(3)
320 ! --> Principal vector*Principal strains*Transpose(Principal vector)
321 DO k = 1,3
322 DO l = 1,3
323 epsptens2(k,l) = zero
324 DO m = 1,3
325 epsptens2(k,l) = epsptens2(k,l) + epsptens(k,m)*vect_pr(l,m)
326 ENDDO
327 ENDDO
328 ENDDO
329 DO k = 1,3
330 DO l = 1,3
331 epsptens(k,l) = zero
332 DO m = 1,3
333 epsptens(k,l) = epsptens(k,l) + vect_pr(k,m)*epsptens2(m,l)
334 ENDDO
335 ENDDO
336 ENDDO
337 ! -> Save true strain rates
338 epsp11(i) = epsptens(1,1)
339 epsp22(i) = epsptens(2,2)
340 epsp33(i) = epsptens(3,3)
341 epsp12(i) = epsptens(1,2)
342 epsp23(i) = epsptens(2,3)
343 epsp31(i) = epsptens(3,1)
344 END IF
345 ENDDO
346c
347 CASE DEFAULT
348 ! no transformation : failure strain measure is defined by Ismstr
349 eps11(1:nel) = epsxx(1:nel)
350 eps22(1:nel) = epsyy(1:nel)
351 eps33(1:nel) = epszz(1:nel)
352 eps12(1:nel) = half*epsxy(1:nel)
353 eps23(1:nel) = half*epsyz(1:nel)
354 eps31(1:nel) = half*epszx(1:nel)
355 epsp11(1:nel) = epspxx(1:nel)
356 epsp22(1:nel) = epspyy(1:nel)
357 epsp33(1:nel) = epspzz(1:nel)
358 epsp12(1:nel) = half*epspxy(1:nel)
359 epsp23(1:nel) = half*epspyz(1:nel)
360 epsp31(1:nel) = half*epspzx(1:nel)
361 END SELECT
362c------------------------------
363c Element size scale factor
364 IF (ifunc_siz > 0) THEN
365 IF (uvar(1,13)==zero) THEN
366 escal(1:nel) = deltax(1:nel) / ref_siz
367 iposp(1:nel) = 0
368 iadp(1:nel) = npf(ifunc_siz) / 2 + 1
369 ilenp(1:nel) = npf(ifunc_siz + 1) / 2 - iadp(1:nel) - iposp(1:nel)
370 CALL vinter2(tf,iadp,iposp,ilenp,nel,escal,dydx,fsize)
371 fsize(1:nel) = fscale_siz * fsize(1:nel)
372 uvar(1:nel,13) = fsize(1:nel)
373 ELSE
374 fsize(1:nel) = uvar(1:nel,13)
375 ENDIF
376 ELSE
377 fsize(1:nel) = one
378 ENDIF
379c------------------------------
380c
381 xx=1
382 yy=2
383 zz=3
384 xy=4
385 yz=5
386 zx=6
387 m11t = 1
388 m11c = 4
389 m22t = 2
390 m22c = 5
391 m33t = 7
392 m33c = 8
393 m12t = 3
394 m12c = 6
395 m23t = 9
396 m23c = 10
397 m31t = 11
398 m31c = 12
399c current variable initialisation
400 DO i=1, nel
401 IF (off(i) == one) THEN
402 DO j=1 ,12
403 odam(i,j) = min(dfmax(i,j+1),one-em3)
404 ENDDO
405 DO j=1 ,6
406 sigo(i,j) = uvar(i,j)
407 epsp(i,j) = uvar(i,j+6)
408 ENDDO
409 ENDIF
410 ENDDO
411c
412c strain rate filtering
413 DO i=1,nel
414 IF (off(i) == one) THEN
415 epsp(i,xx) = alpha * abs(epsp11(i)) + (one-alpha)*epsp(i,xx)
416 epsp(i,yy) = alpha * abs(epsp22(i)) + (one-alpha)*epsp(i,yy)
417 epsp(i,zz) = alpha * abs(epsp33(i)) + (one-alpha)*epsp(i,zz)
418 epsp(i,xy) = alpha * abs(epsp12(i)) + (one-alpha)*epsp(i,xy)
419 epsp(i,yz) = alpha * abs(epsp23(i)) + (one-alpha)*epsp(i,yz)
420 epsp(i,zx) = alpha * abs(epsp31(i)) + (one-alpha)*epsp(i,zx)
421 ELSE
422 epsp(i,xx) = zero
423 epsp(i,yy) = zero
424 epsp(i,zz) = zero
425 epsp(i,xy) = zero
426 epsp(i,yz) = zero
427 epsp(i,zx) = zero
428 ENDIF
429 ENDDO
430c
431c damage and failure calculation in each mode (dir)
432c
433 nindx = 0
434 indx(1:nel) = 0
435 fmode(1:nel,1:12) = 0
436 DO i=1, nel
437 IF (off(i) == one) THEN
438 et1 = uparam(4)
439 et1m = uparam(5)
440 et2 = uparam(6)
441 et2m = uparam(7)
442 ec1 = uparam(8)
443 ec1m = uparam(9)
444 ec2 = uparam(10)
445 ec2m = uparam(11)
446 et12 = uparam(12)
447 et12m = uparam(13)
448 ec12 = uparam(14)
449 ec12m = uparam(15)
450 et3 = uparam(18)
451 et3m = uparam(19)
452 ec3 = uparam(20)
453 ec3m = uparam(21)
454 et23 = uparam(22)
455 et23m = uparam(23)
456 ec23 = uparam(24)
457 ec23m = uparam(25)
458 et31 = uparam(22)
459 et31m = uparam(23)
460 ec31 = uparam(24)
461 ec31m = uparam(25)
462c
463 faili = 0
464c
465c MODE XX TRACTION
466c
467 mode = 1
468 epdot = abs(epsp(i,xx))/epspref
469 IF (ifunc(mode) > 0) THEN
470 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
471 ELSE
472 frate = one
473 ENDIF
474 frate = frate*fsize(i)
475 et1m = frate * et1m
476 et1 = frate * et1
477 deps = max(eps11(i) - et1, zero)
478 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
479 depsm = (et1m - et1)
480 fact = et1m / max(eps11(i),em20)
481 di = fact*deps/max(depsm,em20)
482 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
483 IF (dfmax(i,1+mode) >= one) THEN
484 faili = 1
485 fmode(i,mode) = 1
486 dfmax(i,1+mode) = one
487 epsm1(i) = et1m
488 ENDIF
489 ENDIF
490c
491c MODE YY TRACTION
492c
493 mode = 2
494 epdot = abs(epsp(i,yy))/epspref
495 IF (ifunc(mode) > 0) THEN
496 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
497 ELSE
498 frate = one
499 ENDIF
500 frate = frate*fsize(i)
501 et2m = frate * et2m
502 et2 = frate * et2
503 deps = max(eps22(i) - et2, zero)
504 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
505 depsm = (et2m - et2)
506 fact = et2m / max(eps22(i),em20)
507 di = fact * deps / max(depsm,em20)
508 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
509 IF (dfmax(i,1+mode) >= one) THEN
510 faili = 1
511 fmode(i,mode) = 1
512 dfmax(i,1+mode) = one
513 epsm2(i) = et2m
514 ENDIF
515 ENDIF
516c
517c MODE XY TRACTION
518c
519 mode = 3
520 epdot = abs(epsp(i,xy))/epspref
521 IF (ifunc(mode) > 0) THEN
522 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
523 ELSE
524 frate = one
525 ENDIF
526 frate = frate*fsize(i)
527 et12m = frate * et12m
528 et12 = frate * et12
529 deps = max(eps12(i) - et12, zero)
530 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
531 depsm = (et12m - et12)
532 fact = et12m / max(eps12(i),em20)
533 di = fact * deps / max(depsm,em20)
534 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
535 IF (dfmax(i,1+mode) >= one) THEN
536 faili = 1
537 fmode(i,mode) = 1
538 dfmax(i,1+mode) = one
539 epsm3(i) = et12m
540 ENDIF
541 ENDIF
542c
543c MODE XX COMPRESSION
544c
545 mode = 4
546 epdot = abs(epsp(i,xx))/epspref
547 IF (ifunc(mode) > 0) THEN
548 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
549 ELSE
550 frate = one
551 ENDIF
552 frate = frate*fsize(i)
553 ec1m = frate * ec1m
554 ec1 = frate * ec1
555 deps = max(-eps11(i) - ec1, zero)
556 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
557 depsm = (ec1m - ec1)
558 fact = ec1m / max(abs(eps11(i)),em20)
559 di = fact * deps / max(depsm,em20)
560 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
561 IF (dfmax(i,1+mode) >= one) THEN
562 faili = 1
563 fmode(i,mode) = 1
564 dfmax(i,1+mode) = one
565 epsm4(i) = -ec1m
566 ENDIF
567 ENDIF
568c
569c MODE YY COMPRESSION
570c
571 mode = 5
572 epdot = abs(epsp(i,yy))/epspref
573 IF (ifunc(mode) > 0) THEN
574 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
575 ELSE
576 frate = one
577 ENDIF
578 frate = frate*fsize(i)
579 ec2m = frate * ec2m
580 ec2 = frate * ec2
581 deps = max(-eps22(i) - ec2, zero)
582 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
583 depsm = (ec2m - ec2)
584 fact = ec2m / max(abs(eps22(i)),em20)
585 di = fact * deps / max(depsm,em20)
586 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
587 IF (dfmax(i,1+mode) >= one) THEN
588 faili = 1
589 fmode(i,mode) = 1
590 dfmax(i,1+mode) = one
591 epsm5(i) = -ec2m
592 ENDIF
593 ENDIF
594c
595c MODE XY COMPRESSION
596c
597 mode = 6
598 epdot = abs(epsp(i,xy))/epspref
599 IF (ifunc(mode) > 0) THEN
600 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
601 ELSE
602 frate = one
603 ENDIF
604 frate = frate*fsize(i)
605 ec12m = frate * ec12m
606 ec12 = frate * ec12
607 deps = max(-eps12(i) - ec12, zero)
608 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
609 depsm = (ec12m - ec12)
610 fact = ec12m / max(abs(eps12(i)),em20)
611 di = fact * deps / max(depsm,em20)
612 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
613 IF (dfmax(i,1+mode) >= one) THEN
614 faili = 1
615 fmode(i,mode) = 1
616 dfmax(i,1+mode) = one
617 epsm6(i) = -ec12m
618 ENDIF
619 ENDIF
620c
621c MODE ZZ Traction
622c
623 mode = m33t
624 epdot = abs(epsp(i,zz))/epspref
625 IF (ifunc(mode) > 0) THEN
626 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
627 ELSE
628 frate = one
629 ENDIF
630 frate = frate*fsize(i)
631 et3m = frate * et3m
632 et3 = frate * et3
633 deps = max(eps33(i) - et3, zero)
634 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
635 depsm = (et3m - et3)
636 fact = et3m / max(eps33(i),em20)
637 di = fact * deps / max(depsm,em20)
638 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
639 IF (dfmax(i,1+mode) >= one) THEN
640 faili = 1
641 fmode(i,mode) = 1
642 dfmax(i,1+mode) = one
643 epsm7(i) = et3m
644 ENDIF
645 ENDIF
646c
647c MODE ZZ Compression
648c
649 mode = m33c
650 epdot = abs(epsp(i,zz))/epspref
651 IF (ifunc(mode) > 0) THEN
652 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
653 ELSE
654 frate = one
655 ENDIF
656 frate = frate*fsize(i)
657 ec3m = frate * ec3m
658 ec3 = frate * ec3
659 deps = max(-eps33(i) - ec3, zero)
660 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
661 depsm = (ec3m - ec3)
662 fact = ec3m / max(abs(eps33(i)),em20)
663 di = fact * deps / max(depsm,em20)
664 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
665 IF (dfmax(i,1+mode) >= one) THEN
666 faili = 1
667 fmode(i,mode) = 1
668 dfmax(i,1+mode) = one
669 epsm8(i) = -ec3m
670 ENDIF
671 ENDIF
672c
673c MODE YZ TRACT
674c
675 mode = m23t
676 epdot = abs(epsp(i,yz))/epspref
677 IF (ifunc(mode) > 0) THEN
678 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
679 ELSE
680 frate = one
681 ENDIF
682 frate = frate*fsize(i)
683 et23m = frate * et23m
684 et23 = frate * et23
685 deps = max(eps23(i) - et23, zero)
686 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
687 depsm = (et23m - et23)
688 fact = et23m / max(eps23(i),em20)
689 di = fact * deps / max(depsm,em20)
690 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
691 IF (dfmax(i,1+mode) >= one) THEN
692 faili = 1
693 fmode(i,mode) = 1
694 dfmax(i,1+mode) = one
695 epsm9(i) = et23m
696 ENDIF
697 ENDIF
698c
699c MODE YZ COMP
700c
701 mode = m23c
702 epdot = abs(epsp(i,yz))/epspref
703 IF (ifunc(mode) > 0) THEN
704 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
705 ELSE
706 frate = one
707 ENDIF
708 frate = frate*fsize(i)
709 ec23m = frate * ec23m
710 ec23 = frate * ec23
711 deps = max(-eps23(i) - ec23, zero)
712 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
713 depsm = (ec23m - ec23)
714 fact = ec23m / max(abs(eps23(i)),em20)
715 di = fact * deps / max(depsm,em20)
716 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
717 IF (dfmax(i,1+mode) >= one) THEN
718 faili = 1
719 fmode(i,mode) = 1
720 dfmax(i,1+mode) = one
721 epsm10(i) = -ec23m
722 ENDIF
723 ENDIF
724c
725c MODE ZX Traction
726c
727 mode = m31t
728 epdot = abs(epsp(i,zx))/epspref
729 IF (ifunc(mode) > 0) THEN
730 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
731 ELSE
732 frate = one
733 ENDIF
734 frate = frate*fsize(i)
735 et31m = frate * et31m
736 et31 = frate * et31
737 deps = max(eps31(i) - et31, zero)
738 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
739 depsm = (et31m - et31)
740 fact = et31m / max(eps31(i),em20)
741 di = fact * deps / max(depsm,em20)
742 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
743 IF (dfmax(i,1+mode) >= one) THEN
744 faili = 1
745 fmode(i,mode) = 1
746 dfmax(i,1+mode) = one
747 epsm11(i) = et31m
748 ENDIF
749 ENDIF
750c
751c mode ZX comp
752c
753 mode = m31c ! Mode 12
754 epdot = abs(epsp(i,zx))/epspref
755 IF (ifunc(mode) > 0) THEN
756 frate = finter(ifunc(mode),epdot,npf,tf,dydx)
757 ELSE
758 frate = one
759 ENDIF
760 frate = frate*fsize(i)
761 ec31m = frate * ec31m
762 ec31 = frate * ec31
763 deps = max(-eps31(i) - ec31, zero)
764 IF (deps > zero .AND. dfmax(i,1+mode) < one) THEN
765 depsm = (ec31m - ec31)
766 fact = ec31m / max(abs(eps31(i)),em20)
767 di = fact * deps / max(depsm,em20)
768 dfmax(i,1+mode) = max(dfmax(i,1+mode), di)
769 IF (dfmax(i,1+mode) >= one) THEN
770 faili = 1
771 fmode(i,mode) = 1
772 dfmax(i,1+mode) = one
773 epsm12(i) = -ec31m
774 ENDIF
775 ENDIF
776c
777c
778 IF (faili == 1) THEN
779 faili = 0
780 uelr(i) = uelr(i)+1
781 nindx = nindx + 1
782 indx(nindx) = i
783 IF (uelr(i) >= npt) THEN
784 tdel(i) = time
785 off(i) = zero
786 ENDIF
787 ENDIF
788 ENDIF
789 ENDDO
790c------------------------
791c Apply damage to stress
792c------------------------
793 DO i=1, nel
794 IF (off(i) == one) THEN
795c damage / direction is equal to the max of damages par mode in the same direction
796 damxx = max(dfmax(i,2) ,dfmax(i,5) )
797 damyy = max(dfmax(i,3) ,dfmax(i,6) )
798 damzz = max(dfmax(i,8) ,dfmax(i,9) )
799 damxy = max(dfmax(i,4) ,dfmax(i,7) )
800 damyz = max(dfmax(i,10),dfmax(i,11))
801 damzx = max(dfmax(i,12),dfmax(i,13))
802c
803 odamxx = max(odam(i,1) ,odam(i,4) )
804 odamyy = max(odam(i,2) ,odam(i,5) )
805 odamzz = max(odam(i,7) ,odam(i,8) )
806 odamxy = max(odam(i,3) ,odam(i,6) )
807 odamyz = max(odam(i,9) ,odam(i,10))
808 odamzx = max(odam(i,11),odam(i,12))
809c
810 dfmax(i,1) = max(damxx,damyy)
811 dfmax(i,1) = max(damzz,dfmax(i,1))
812 dfmax(i,1) = max(damxy,dfmax(i,1))
813 dfmax(i,1) = max(damyz,dfmax(i,1))
814 dfmax(i,1) = max(damzx,dfmax(i,1))
815c Undamaged sign estimate
816 signxx(i) = sigo(i,xx)/max(one-odamxx,em20) + (signxx(i) - sigo(i,xx))
817 signyy(i) = sigo(i,yy)/max(one-odamyy,em20) + (signyy(i) - sigo(i,yy))
818 signzz(i) = sigo(i,zz)/max(one-odamzz,em20) + (signzz(i) - sigo(i,zz))
819 signxy(i) = sigo(i,xy)/max(one-odamxy,em20) + (signxy(i) - sigo(i,xy))
820 signyz(i) = sigo(i,yz)/max(one-odamyz,em20) + (signyz(i) - sigo(i,yz))
821 signzx(i) = sigo(i,zx)/max(one-odamzx,em20) + (signzx(i) - sigo(i,zx))
822c Damaged sign
823 signxx(i) = signxx(i) * (one-damxx)*off(i)
824 signyy(i) = signyy(i) * (one-damyy)*off(i)
825 signzz(i) = signzz(i) * (one-damzz)*off(i)
826 signxy(i) = signxy(i) * (one-damxy)*off(i)
827 signyz(i) = signyz(i) * (one-damyz)*off(i)
828 signzx(i) = signzx(i) * (one-damzx)*off(i)
829 ELSE
830 signxx(i) = zero
831 signyy(i) = zero
832 signzz(i) = zero
833 signxy(i) = zero
834 signyz(i) = zero
835 signzx(i) = zero
836 ENDIF
837 ENDDO
838c------------------------
839c Update UVAR
840c------------------------
841 DO i=1, nel
842 IF (off(i) == one) THEN
843 DO j=1,6
844 uvar(i,j+6) = epsp(i,j)
845 ENDDO
846 ENDIF
847 ENDDO
848 DO i =1, nel
849 IF (off(i) == one) THEN
850 uvar(i,xx) = signxx(i)
851 uvar(i,yy) = signyy(i)
852 uvar(i,zz) = signzz(i)
853 uvar(i,xy) = signxy(i)
854 uvar(i,yz) = signyz(i)
855 uvar(i,zx) = signzx(i)
856 ENDIF
857 ENDDO
858c------------------------
859c Output
860c------------------------
861 IF (nindx > 0) THEN
862 DO j=1,nindx
863 i = indx(j)
864#include "lockon.inc"
865 WRITE(istdo,1000) ngl(i),ipt
866 WRITE(iout, 2000) ngl(i),ipt
867 IF (fmode(i,1) == 1) WRITE(iout,4000) '1 - TRACTION XX',eps11(i),epsm1(i),epsp(i,xx)
868 IF (fmode(i,2) == 1) WRITE(iout,4000) '2 - TRACTION YY',eps22(i),epsm2(i),epsp(i,yy)
869 IF (fmode(i,3) == 1) WRITE(iout,4000) '3 - POSITIVE SHEAR XY',eps12(i),epsm3(i),epsp(i,xy)
870 IF (fmode(i,4) == 1) WRITE(iout,4000) '4 - COMPRESSION XX',eps11(i),epsm4(i),epsp(i,xx)
871 IF (fmode(i,5) == 1) WRITE(iout,4000) '5 - COMPRESSION YY',eps22(i),epsm5(i),epsp(i,yy)
872 IF (fmode(i,6) == 1) WRITE(iout,4000) '6 - NEGATIVE SHEAR XY',eps12(i),epsm6(i),epsp(i,xy)
873 IF (fmode(i,7) == 1) WRITE(iout,4000) '7 - TRACTION ZZ',eps33(i),epsm7(i),epsp(i,zz)
874 IF (fmode(i,8) == 1) WRITE(iout,4000) '8 - COMPRESSION ZZ',eps33(i),epsm8(i),epsp(i,zz)
875 IF (fmode(i,9) == 1) WRITE(iout,4000) '9 - POSITIVE SHEAR YZ',eps23(i),epsm9(i),epsp(i,yz)
876 IF (fmode(i,10) == 1) WRITE(iout,4000) '10 - NEGATIVE SHEAR YZ',eps23(i),epsm10(i),epsp(i,yz)
877 IF (fmode(i,11) == 1) WRITE(iout,4000) '11 - POSITIVE SHEAR ZX',eps31(i),epsm11(i),epsp(i,zx)
878 IF (fmode(i,12) == 1) WRITE(iout,4000) '12 - NEGATIVE SHEAR ZX',eps31(i),epsm12(i),epsp(i,zx)
879 IF (off(i) == zero) WRITE(iout, 3000) ngl(i),time
880#include "lockoff.inc"
881 END DO
882 END IF
883c-----------------------------------------------------------------------
884 1000 FORMAT(1x,'FAILURE (ORTHSTR) OF SOLID ELEMENT ',i10,1x,',GAUSS PT ',i5)
885 2000 FORMAT(1x,'FAILURE (ORTHSTR) OF SOLID ELEMENT ',i10,1x,',GAUSS PT ',i5)
886 3000 FORMAT(1x,'RUPTURE OF ELEMENT #',i10,1x,'AT TIME # ',1pe12.4)
887 4000 FORMAT(1x,'MODE',1x,a,', STRAIN',1pe12.4,', CRITICAL VALUE',1pe12.4,
888 . ', STRAIN RATE',1pe12.4)
889c-----------------------------------------------------------------------
890 RETURN
891 END
alpha
#define alpha
Definition eval.h:35
fail_orthstrain
subroutine fail_orthstrain(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, ismstr, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, signxx, signyy, signzz, signxy, signyz, signzx, uvar, off, ipt, ngl, dfmax, tdel, uelr, npt, deltax, lf_dammx)
Definition fail_orthstrain_s.F:42
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
jacobiew
subroutine jacobiew(a, n, ew, ev, nrot)
Definition matrix.F:29
mmain
subroutine mmain(pm, elbuf_str, ix, nix, x, geo, iparg, nel, skew, bufmat, ipart, ipartel, nummat, matparam, imat, ipm, ngl, pid, npf, tf, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, rx, ry, rz, sx, sy, sz, gama, voln, dvol, s1, s2, s3, s4, s5, s6, dxx, dyy, dzz, d4, d5, d6, wxx, wyy, wzz)
Definition mmain.F:43
mulaw
subroutine mulaw(lft, llt, nft, mtn, jcvt, pm, off, sig, eint, rho, vol, strain, gama, uvar, bufmat, tf, npf, imat, ngl, nuvar, nvartmp, vartmp, geo, pid, epsd, wxx, wyy, wzz, jsph, ssp, voln, vis, d1, d2, d3, d4, d5, d6, dvol, sold1, sold2, sold3, sold4, sold5, sold6, rx, ry, rz, sx, sy, sz, tx, ty, tz, ismstr, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, ipm, isorth, nel, matparam)
Definition mulaw.F:54
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143