OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps74.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!|| sigeps74 ../engine/source/materials/mat/mat074/sigeps74.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| finter2 ../engine/source/tools/curve/vinter.F
30!|| table_vinterp ../engine/source/tools/curve/table_tools.F
31!||--- uses -----------------------------------------------------
32!|| interface_table_mod ../engine/share/modules/table_mod.F
33!|| table_mod ../engine/share/modules/table_mod.F
34!||====================================================================
35 SUBROUTINE sigeps74(
36 1 NEL ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,NPF ,
37 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,RHO ,
38 3 VOLUME ,EINT ,FHEAT ,JLAG ,
39 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
40 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
41 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
42 7 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
43 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
44 9 SIGVXX ,SIGVYY ,SIGVZZ ,SIGVXY ,SIGVYZ ,SIGVZX ,
45 A SOUNDSP,VISCMAX,UVAR ,OFF ,NGL ,IPT ,
46 B IPM ,MAT ,EPSP ,IPLA ,YLD ,PLA ,
47 C DPLA ,ETSE ,JTHE ,TEMP ,TABLE ,SEQ_OUTPUT,
48 D AMU ,ISEQ )
49C-----------------------------------------------
50 USE table_mod
51 USE interface_table_mod
52C-----------------------------------------------
53C I m p l i c i t T y p e s
54C-----------------------------------------------
55#include "implicit_f.inc"
56#include "comlock.inc"
57C-----------------------------------------------
58C G l o b a l P a r a m e t e r s
59C-----------------------------------------------
60#include "mvsiz_p.inc"
61C---------+---------+---+---+--------------------------------------------
62C VAR | SIZE |TYP| RW| DEFINITION
63C---------+---------+---+---+--------------------------------------------
64C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
65C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
66C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
67C---------+---------+---+---+--------------------------------------------
68C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
69C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
70C NPF | * | I | R | FUNCTION ARRAY
71C TF | * | F | R | FUNCTION ARRAY
72C---------+---------+---+---+--------------------------------------------
73C TIME | 1 | F | R | CURRENT TIME
74C TIMESTEP| 1 | F | R | CURRENT TIME STEP
75C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
76C RHO0 | NEL | F | R | INITIAL DENSITY
77C RHO | NEL | F | R | DENSITY
78C VOLUME | NEL | F | R | VOLUME
79C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
80C EPSPXX | NEL | F | R | STRAIN RATE XX
81C EPSPYY | NEL | F | R | STRAIN RATE YY
82C ... | | | |
83C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
84C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
85C ... | | | |
86C EPSXX | NEL | F | R | STRAIN XX
87C EPSYY | NEL | F | R | STRAIN YY
88C ... | | | |
89C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
90C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
91C ... | | | |
92C---------+---------+---+---+--------------------------------------------
93C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
94C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
95C ... | | | |
96C SIGVXX | NEL | F | W | VISCOUS STRESS XX
97C SIGVYY | NEL | F | W | VISCOUS STRESS YY
98C ... | | | |
99C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
100C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
101C---------+---------+---+---+--------------------------------------------
102C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
103C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
104C---------+---------+---+---+--------------------------------------------
105#include "com01_c.inc"
106#include "com08_c.inc"
107#include "impl1_c.inc"
108#include "param_c.inc"
109#include "scr17_c.inc"
110#include "units_c.inc"
111C=======================================================================
112 INTEGER NEL, NUPARAM, NUVAR,IPT,IPLA, JTHE,
113 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*),ISEQ
114 INTEGER ,INTENT(IN) :: JLAG
115 my_real
116 . TIME,TIMESTEP,UPARAM(*),
117 . RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
118 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
119 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
120 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
121 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
122 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
123 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
124 . SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
125 . sigoxy(nel),sigoyz(nel),sigozx(nel),
126 . epsp(nel),etse(nel), temp(nel),seq_output(nel),
127 . amu(nel)
128 TYPE(ttable) TABLE(*)
129C-----------------------------------------------
130C O U T P U T A r g u m e n t s
131C-----------------------------------------------
132 my_real
133 . SIGNXX(NEL),SIGNYY(NEL),SIGNZZ(NEL),
134 . SIGNXY(NEL),SIGNYZ(NEL),SIGNZX(NEL),
135 . SIGVXX(NEL),SIGVYY(NEL),SIGVZZ(NEL),
136 . SIGVXY(NEL),SIGVYZ(NEL),SIGVZX(NEL),
137 . SOUNDSP(NEL),VISCMAX(NEL),DPLA(NEL)
138C-----------------------------------------------
139C I N P U T O U T P U T A r g u m e n t s
140C-----------------------------------------------
141 my_real :: uvar(nel,nuvar), off(nel), yld(nel), pla(nel)
142 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: fheat
143C-----------------------------------------------
144C VARIABLES FOR FUNCTION INTERPOLATION
145C-----------------------------------------------
146 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
147 my_real
148 . FINTER2, TF(*),FINTER
149 EXTERNAL FINTER2,FINTER
150C-----------------------------------------------
151C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
152C Y : y = f(x)
153C X : x
154C DYDX : f'(x) = dy/dx
155C IFUNC(J): FUNCTION INDEX
156C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
157C NPF,TF : FUNCTION PARAMETER
158C-----------------------------------------------
159C L o c a l V a r i a b l e s
160C-----------------------------------------------
161 INTEGER I,J,IADBUF,J1,J2,ITABLE,
162 . iparam,npar,opte,ifunce,
163 . nindx,indx(mvsiz)
164 my_real
165 . e,nu,p,dav,r,fac,epst,ep1,ep2,ep3,ep4,ep5,ep6,
166 . e1,e2,e3,e4,e5,e6,c,cc,d,y,yp,e42,e52,e62,epst2,
167 . c1,g,g2,g3,ff,gg,hh,nn,ll,mm,vol0, rhocp,
168 . epsmax,h(mvsiz),cri(mvsiz),
169 . fail(mvsiz),epsr1,epsr2,p0(mvsiz),hkin,
170 . fisokin,sigexx(mvsiz),sigeyy(mvsiz),sigezz(mvsiz),
171 . sigexy(mvsiz),sigeyz(mvsiz),sigezx(mvsiz),
172 . dexx,deyy,dexy,dezz,deyz,dezx,dsxx,dsyy,dszz,dsxy,dsyz,
173 . dszx,sigpxx,sigpyy,sigpxy,sigpyz,sigpzx,sigpzz,alpha,
174 . xfac, yfac,e0(mvsiz), g1(mvsiz), g21(mvsiz), g31(mvsiz),
175 . c11(mvsiz),escale(mvsiz),einf,dydxe(mvsiz),ce
176 my_real
177 . yk(mvsiz), dydx(mvsiz)
178 my_real,
179 . DIMENSION(NEL,3) :: xx3
180 my_real,
181 . DIMENSION(NEL,2) :: xx2
182 INTEGER,DIMENSION(NEL,3) :: IPOS3
183 INTEGER,DIMENSION(NEL,2) :: IPOS2
184C-----------------------------------------------
185C USER VARIABLES INITIALIZATION
186C-----------------------------------------------
187
188 DO I=1,nel
189 etse(i) = one
190 ENDDO
191C
192 itable = ipm(227,mat(1))
193 iadbuf = ipm(7,mat(1))-1
194 fisokin=uparam(iadbuf+1)
195 rhocp = uparam(iadbuf+23)
196 opte = uparam(iadbuf+25)
197 ce = uparam(iadbuf+27)
198 einf = uparam(iadbuf+26)
199
200 IF (opte==1.OR. ce > zero)THEN
201 DO i=1,nel
202 e0(i) = uparam(iadbuf+2)
203 g1(i) = uparam(iadbuf+5)
204 g21(i) = uparam(iadbuf+18)
205 g31(i) = uparam(iadbuf+19)
206 c11(i) = uparam(iadbuf+20)
207 ENDDO
208 nu = uparam(iadbuf+6)
209 epsmax= uparam(iadbuf+15)
210 IF(epsmax==zero)THEN
211c NXK=SIZE(TABLE(ITABLE)%X(1)%VALUES)
212c IF(TABLE(ITABLE)%Y%VALUES(NXK)==ZERO)THEN
213c EPSMAX = TABLE(ITABLE)%X(1)%VALUES(NXK)
214c ELSE
215 epsmax = ep30
216c ENDIF
217 ENDIF
218C
219 IF (opte == 1)THEN
220 DO i=1,nel
221 pla(i) = uvar(i,1)
222 IF(pla(i) > zero)THEN
223 ifunce = uparam(iadbuf+24)
224 escale(i) = finter(kfunc(ifunce),pla(i),npf,tf,dydxe(i))
225 e0(i) = escale(i)* e0(i)
226 g1(i) = half*e0(i)/(one+nu)
227 g21(i) = two*g1(i)
228 g31(i) = three*g1(i)
229 c11(i) =e0(i)/three/(one - two*nu)
230 ENDIF
231 ENDDO
232 ELSEIF ( ce /= zero) THEN
233 DO i=1,nel
234 pla(i) = uvar(i,1)
235 IF(pla(i) > zero)THEN
236 e0(i) = e0(i)-(e0(i)-einf)*(one-exp(-ce*pla(i)))
237 g1(i) = half*e0(i)/(one+nu)
238 g21(i) = two*g1(i)
239 g31(i) = three*g1(i)
240 c11(i) = e0(i)/three/(one - two*nu)
241 ENDIF
242 ENDDO
243 ENDIF
244 xfac =uparam(iadbuf+13)
245 yfac =uparam(iadbuf+14)
246C
247 epsr1 =uparam(iadbuf+16)
248 IF(epsr1==zero)epsr1=ep30
249 epsr2 =uparam(iadbuf+17)
250 IF(epsr2==zero)epsr2=twoep30
251C
252 ff = uparam(iadbuf+7)
253 gg = uparam(iadbuf+8)
254 hh = uparam(iadbuf+9)
255 ll = uparam(iadbuf+10)
256 mm = uparam(iadbuf+11)
257 nn = uparam(iadbuf+12)
258C
259 IF (time == zero) THEN
260 temp(1:nel) = uparam(iadbuf + 22)
261 IF (isigi==0) THEN
262 DO i=1,nel
263 uvar(i,1)= zero
264 uvar(i,2)= zero
265 uvar(i,3)= zero
266 uvar(i,4)= zero
267 DO j=1,6
268 uvar(i,4 + j) = zero
269 ENDDO
270 ENDDO
271 ENDIF
272 ENDIF
273C------------------------------------------
274C ECROUISSAGE CINE
275C------------------------------------------
276 IF (fisokin/=zero) THEN
277 DO i=1,nel
278 sigoxx(i) = sigoxx(i) - uvar(i,5)
279 sigoyy(i) = sigoyy(i) - uvar(i,6)
280 sigozz(i) = sigozz(i) - uvar(i,7)
281 sigoxy(i) = sigoxy(i) - uvar(i,8)
282 sigoyz(i) = sigoyz(i) - uvar(i,9)
283 sigozx(i) = sigozx(i) - uvar(i,10)
284 ENDDO
285 ENDIF
286C------------------------------------------
287 DO i=1,nel
288C
289 pla(i) = uvar(i,1)
290 p0(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third
291 dav = (depsxx(i)+depsyy(i)+depszz(i))*third
292 signxx(i)=sigoxx(i)+p0(i)+g21(i)*(depsxx(i)-dav)
293 signyy(i)=sigoyy(i)+p0(i)+g21(i)*(depsyy(i)-dav)
294 signzz(i)=sigozz(i)+p0(i)+g21(i)*(depszz(i)-dav)
295 signxy(i)=sigoxy(i)+g1(i) *depsxy(i)
296 signyz(i)=sigoyz(i)+g1(i) *depsyz(i)
297 signzx(i)=sigozx(i)+g1(i) *depszx(i)
298C
299 sigexx(i) = signxx(i)
300 sigeyy(i) = signyy(i)
301 sigezz(i) = signzz(i)
302 sigexy(i) = signxy(i)
303 sigeyz(i) = signyz(i)
304 sigezx(i) = signzx(i)
305 soundsp(i) = sqrt((c11(i)+four*g1(i)/three)/rho0(i))
306 viscmax(i) = zero
307 dpla(i) =zero
308
309 ENDDO
310C-------------------
311C STRAIN RATE
312C-------------------
313 DO i=1,nel
314C-------------------
315C STRAIN principal 1, 4 newton iterations
316C-------------------
317 dav = (epsxx(i)+epsyy(i)+epszz(i)) * third
318 e1 = epsxx(i) - dav
319 e2 = epsyy(i) - dav
320 e3 = epszz(i) - dav
321 e4 = half*epsxy(i)
322 e5 = half*epsyz(i)
323 e6 = half*epszx(i)
324C -y = (e1-x)(e2-x)(e3-x)
325C - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
326C + 2e4 e5 e6
327C e1 + e2 + e3 = 0 => terme en x^2 = 0
328C y = x^3 + c x + d
329c yp= 3 x^2 + c
330 e42 = e4*e4
331 e52 = e5*e5
332 e62 = e6*e6
333 c = - half *(e1*e1 + e2*e2 + e3*e3) - e42 - e52 - e62
334 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
335 & - two*e4*e5*e6
336 cc = c*third
337 epst = sqrt(-cc)
338 epst2 = epst * epst
339 y = (epst2 + c)* epst + d
340 IF(abs(y)>em8)THEN
341 epst = onep75 * epst
342 epst2 = epst * epst
343 y = (epst2 + c)* epst + d
344 yp = three*epst2 + c
345 IF(yp/=zero)epst = epst - y/yp
346 epst2 = epst * epst
347 y = (epst2 + c)* epst + d
348 yp = three*epst2 + c
349 IF(yp/=zero)epst = epst - y/yp
350 epst2 = epst * epst
351 y = (epst2 + c)* epst + d
352 yp = three*epst2 + c
353 IF(yp/=zero)epst = epst - y/yp
354 epst2 = epst * epst
355 y = (epst2 + c)* epst + d
356 yp = three*epst2 + c
357 IF(yp/=zero)epst = epst - y/yp
358 ENDIF
359 epst = epst + dav
360C-------------------
361C tension failure
362C-------------------
363 fail(i) = max(zero,min(one,
364 . (epsr2-epst)/(epsr2-epsr1) ))
365 ENDDO
366C-------------------
367C CRITERE
368C-------------------
369 IF(table(itable)%NDIM==2)THEN
370 DO i=1,nel
371 ipos2(i,1) = nint(uvar(i,2))
372 ipos2(i,2) = nint(uvar(i,3))
373C
374 xx2(i,1)=pla(i)
375 xx2(i,2)=epsp(i)*xfac
376 END DO
377C
378 CALL table_vinterp(table(itable),nel,nel,ipos2,xx2,yld,dydx)
379C
380 DO i=1,nel
381 yld(i) = yfac*yld(i)
382 yld(i) = fail(i)*yld(i)
383 h(i) = fail(i)*dydx(i)
384 uvar(i,2) = ipos2(i,1)
385 uvar(i,3) = ipos2(i,2)
386 END DO
387C-----
388 IF(fisokin/=zero)THEN
389 DO i=1,nel
390 ipos2(i,1) = 1
391C
392 xx2(i,1)=zero
393 xx2(i,2)=epsp(i)*xfac
394 END DO
395 CALL table_vinterp(table(itable),nel,nel,ipos2,xx2,yk,dydx)
396
397 DO i=1,nel
398C ECROUISSAGE CINEMATIQUE
399 yld(i) = (one-fisokin) * yld(i) +
400 . fisokin * fail(i) * yfac * yk(i)
401 yld(i) = max(yld(i),em20)
402 ENDDO
403 END IF
404 ELSE
405C-----------------------------------------------
406 DO i=1,nel
407 ipos3(i,1) = nint(uvar(i,2))
408 ipos3(i,2) = nint(uvar(i,3))
409 ipos3(i,3) = nint(uvar(i,4))
410C
411 xx3(i,1) = pla(i)
412 xx3(i,2) = epsp(i)*xfac
413 xx3(i,3) = temp(i)
414 END DO
415C
416 CALL table_vinterp(table(itable),nel,nel,ipos3,xx3,yld,dydx)
417C
418 DO i=1,nel
419 yld(i) = yfac*yld(i)
420 yld(i) = fail(i)*yld(i)
421 h(i) = fail(i)*dydx(i)
422 uvar(i,2) = ipos3(i,1)
423 uvar(i,3) = ipos3(i,2)
424 uvar(i,4) = ipos3(i,3)
425 END DO
426C-----
427 IF(fisokin/=zero)THEN
428 DO i=1,nel
429 ipos3(i,1) = 1
430C
431 xx3(i,1)=zero
432 xx3(i,2)=epsp(i)*xfac
433 xx3(i,3)=temp(i)
434 END DO
435 CALL table_vinterp(table(itable),nel,nel,ipos3,xx3,yk,dydx)
436
437 DO i=1,nel
438C ECROUISSAGE CINEMATIQUE
439 yld(i) = (one-fisokin) * yld(i) +
440 . fisokin * fail(i) * yfac * yk(i)
441 yld(i) = max(yld(i),em20)
442 ENDDO
443 END IF
444 END IF
445C-------------------
446C PROJECTION
447C-------------------
448 IF(ipla==0)THEN
449 DO i=1,nel
450 cri(i) =
451 . ff*(signyy(i) - signzz(i))**2 + gg*(signzz(i) - signxx(i))**2
452 . + hh*(signxx(i) - signyy(i))**2 + two*ll*signyz(i)**2
453 . + two*mm*signzx(i)**2 + two*nn*signxy(i)**2
454 cri(i) = sqrt(cri(i))
455 r = min(one,yld(i)/ max(cri(i),em20))
456c P = C11(I) * (RHO(I)/RHO0(I)- ONE)
457 p = c11(i) * amu(i)
458 signxx(i)=signxx(i)*r-p
459 signyy(i)=signyy(i)*r-p
460 signzz(i)=signzz(i)*r-p
461 signxy(i)=signxy(i)*r
462 signyz(i)=signyz(i)*r
463 signzx(i)=signzx(i)*r
464 pla(i)=pla(i) + (one - r)*cri(i)/max(g31(i)+h(i),em20)
465 uvar(i,1)=pla(i)
466 dpla(i) = (one - r)*cri(i)/max(g31(i)+h(i),em20)
467 ENDDO
468 ELSEIF(ipla==2)THEN
469 DO i=1,nel
470 cri(i) =
471 . ff*(signyy(i) - signzz(i))**2 + gg*(signzz(i) - signxx(i))**2
472 . + hh*(signxx(i) - signyy(i))**2 + two*ll*signyz(i)**2
473 . + two*mm*signzx(i)**2 + two*nn*signxy(i)**2
474 cri(i) = sqrt(cri(i))
475 r = min(one,yld(i)/ max(cri(i),em20))
476c P = C11(I) * (RHO(I)/RHO0(I)-ONE)
477 p = c11(i) * amu(i)
478 signxx(i)=signxx(i)*r-p
479 signyy(i)=signyy(i)*r-p
480 signzz(i)=signzz(i)*r-p
481 signxy(i)=signxy(i)*r
482 signyz(i)=signyz(i)*r
483 signzx(i)=signzx(i)*r
484 pla(i)=pla(i) + (one-r)*cri(i)/max(g31(i),em20)
485 uvar(i,1)=pla(i)
486 dpla(i) = (one-r)*cri(i)/max(g31(i),em20)
487 ENDDO
488 ELSEIF(ipla==1)THEN
489 DO i=1,nel
490 cri(i) =
491 . ff*(signyy(i) - signzz(i))**2 + gg*(signzz(i) - signxx(i))**2
492 . + hh*(signxx(i) - signyy(i))**2 + two*ll*signyz(i)**2
493 . + two*mm*signzx(i)**2 + two*nn*signxy(i)**2
494 cri(i) = sqrt(cri(i))
495 r = min(one,yld(i)/ max(cri(i),em20))
496C plastic strain increment.
497 dpla(i)=(one - r)*cri(i)/max(g31(i)+h(i),em20)
498C actual yield stress.
499 yld(i)=max(yld(i)+(one - fisokin)*dpla(i)*h(i),zero)
500 r = min(one,yld(i)/ max(cri(i),em20))
501C
502c P = C11(I) * (RHO(I)/RHO0(I)-ONE)
503 p = c11(i) * amu(i)
504 signxx(i)=signxx(i)*r-p
505 signyy(i)=signyy(i)*r-p
506 signzz(i)=signzz(i)*r-p
507 signxy(i)=signxy(i)*r
508 signyz(i)=signyz(i)*r
509 signzx(i)=signzx(i)*r
510 pla(i)=pla(i) + dpla(i)
511 uvar(i,1)=pla(i)
512 ENDDO
513 ENDIF
514C------------------------------------------
515C ECROUISSAGE CINE
516C------------------------------------------
517 IF (fisokin/=zero) THEN
518 DO i=1,nel
519c P = C11(I) * (RHO(I)/RHO0(I)-ONE)
520 p = c11(i) * amu(i)
521 dsxx = sigexx(i) - signxx(i) - p
522 dsyy = sigeyy(i) - signyy(i) - p
523 dszz = sigezz(i) - signzz(i) - p
524 dsxy = sigexy(i) - signxy(i)
525 dsyz = sigeyz(i) - signyz(i)
526 dszx = sigezx(i) - signzx(i)
527C
528 hkin = two_third*fisokin*h(i)
529 alpha = hkin/(g21(i)+hkin)
530 sigpxx = alpha*dsxx
531 sigpyy = alpha*dsyy
532 sigpzz = alpha*dszz
533 sigpxy = alpha*dsxy
534 sigpyz = alpha*dsyz
535 sigpzx = alpha*dszx
536C
537 uvar(i, 5) = uvar(i, 5) + sigpxx
538 uvar(i, 6) = uvar(i, 6) + sigpyy
539 uvar(i, 7) = uvar(i, 7) + sigpzz
540 uvar(i, 8) = uvar(i, 8) + sigpxy
541 uvar(i, 9) = uvar(i, 9) + sigpyz
542 uvar(i,10) = uvar(i,10) + sigpzx
543C
544 signxx(i) = signxx(i) + uvar(i, 5)
545 signyy(i) = signyy(i) + uvar(i, 6)
546 signzz(i) = signzz(i) + uvar(i, 7)
547 signxy(i) = signxy(i) + uvar(i, 8)
548 signyz(i) = signyz(i) + uvar(i, 9)
549 signzx(i) = signzx(i) + uvar(i,10)
550 ENDDO
551 ENDIF
552
553 IF (impl_s>0) THEN
554 DO i=1,nel
555 IF(dpla(i)>0) etse(i)= h(i)/g21(i)
556 ENDDO
557c IF (IKT==4)
558c . CALL PUT_SIGE0(NEL,SIGEXX,SIGEYY,SIGEZZ,SIGEXY,
559c . SIGEYZ,SIGEZX,DPLA,G ,H )
560 ENDIF
561 ELSE !opte==1 CE/=0
562 itable = ipm(227,mat(1))
563 iadbuf = ipm(7,mat(1))-1
564 fisokin=uparam(iadbuf+1)
565
566 e = uparam(iadbuf+2)
567 g = uparam(iadbuf+5)
568 g2 = uparam(iadbuf+18)
569 g3 = uparam(iadbuf+19)
570 nu = uparam(iadbuf+6)
571 c1 = uparam(iadbuf+20)
572 epsmax= uparam(iadbuf+15)
573 IF(epsmax==zero)THEN
574c NXK=SIZE(TABLE(ITABLE)%X(1)%VALUES)
575c IF(TABLE(ITABLE)%Y%VALUES(NXK)==ZERO)THEN
576c EPSMAX = TABLE(ITABLE)%X(1)%VALUES(NXK)
577c ELSE
578 epsmax = ep30
579c ENDIF
580 ENDIF
581C
582 xfac =uparam(iadbuf+13)
583 yfac =uparam(iadbuf+14)
584C
585 epsr1 =uparam(iadbuf+16)
586 IF(epsr1==zero)epsr1=ep30
587 epsr2 =uparam(iadbuf+17)
588 IF(epsr2==zero)epsr2=twoep30
589C
590 ff = uparam(iadbuf+7)
591 gg = uparam(iadbuf+8)
592 hh = uparam(iadbuf+9)
593 ll = uparam(iadbuf+10)
594 mm = uparam(iadbuf+11)
595 nn = uparam(iadbuf+12)
596C
597 IF (isigi==0) THEN
598 IF(time==zero)THEN
599 DO i=1,nel
600 uvar(i,1)=zero
601 uvar(i,2)=zero
602 uvar(i,3)=zero
603 uvar(i,4)=zero
604 DO j=1,6
605 uvar(i,4 + j ) = zero
606 ENDDO
607 ENDDO
608 ENDIF
609 ENDIF
610C------------------------------------------
611C ECROUISSAGE CINE
612C------------------------------------------
613 IF (fisokin/=zero) THEN
614 DO i=1,nel
615 sigoxx(i) = sigoxx(i) - uvar(i,5)
616 sigoyy(i) = sigoyy(i) - uvar(i,6)
617 sigozz(i) = sigozz(i) - uvar(i,7)
618 sigoxy(i) = sigoxy(i) - uvar(i,8)
619 sigoyz(i) = sigoyz(i) - uvar(i,9)
620 sigozx(i) = sigozx(i) - uvar(i,10)
621 ENDDO
622 ENDIF
623C------------------------------------------
624 DO i=1,nel
625C
626 pla(i) = uvar(i,1)
627 p0(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third
628 dav = (depsxx(i)+depsyy(i)+depszz(i))*third
629 signxx(i)=sigoxx(i)+p0(i)+g2*(depsxx(i)-dav)
630 signyy(i)=sigoyy(i)+p0(i)+g2*(depsyy(i)-dav)
631 signzz(i)=sigozz(i)+p0(i)+g2*(depszz(i)-dav)
632 signxy(i)=sigoxy(i)+g *depsxy(i)
633 signyz(i)=sigoyz(i)+g *depsyz(i)
634 signzx(i)=sigozx(i)+g *depszx(i)
635C
636 sigexx(i) = signxx(i)
637 sigeyy(i) = signyy(i)
638 sigezz(i) = signzz(i)
639 sigexy(i) = signxy(i)
640 sigeyz(i) = signyz(i)
641 sigezx(i) = signzx(i)
642 soundsp(i) = sqrt((c1+four*g/three)/rho0(i))
643 viscmax(i) = zero
644 dpla(i) =zero
645
646 ENDDO
647C-------------------
648C STRAIN RATE
649C-------------------
650 DO i=1,nel
651C-------------------
652C STRAIN principal 1, 4 newton iterations
653C-------------------
654 dav = (epsxx(i)+epsyy(i)+epszz(i)) * third
655 e1 = epsxx(i) - dav
656 e2 = epsyy(i) - dav
657 e3 = epszz(i) - dav
658 e4 = half*epsxy(i)
659 e5 = half*epsyz(i)
660 e6 = half*epszx(i)
661C -y = (e1-x)(e2-x)(e3-x)
662C - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
663C + 2e4 e5 e6
664C e1 + e2 + e3 = 0 => terme en x^2 = 0
665C y = x^3 + c x + d
666c yp= 3 x^2 + c
667 e42 = e4*e4
668 e52 = e5*e5
669 e62 = e6*e6
670 c = - half *(e1*e1 + e2*e2 + e3*e3) - e42 - e52 - e62
671 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
672 & - two*e4*e5*e6
673 cc = c*third
674 epst = sqrt(-cc)
675 epst2 = epst * epst
676 y = (epst2 + c)* epst + d
677 IF(abs(y)>em8)THEN
678 epst = onep75 * epst
679 epst2 = epst * epst
680 y = (epst2 + c)* epst + d
681 yp = three*epst2 + c
682 IF(yp/=zero)epst = epst - y/yp
683 epst2 = epst * epst
684 y = (epst2 + c)* epst + d
685 yp = three*epst2 + c
686 IF(yp/=zero)epst = epst - y/yp
687 epst2 = epst * epst
688 y = (epst2 + c)* epst + d
689 yp = three*epst2 + c
690 IF(yp/=zero)epst = epst - y/yp
691 epst2 = epst * epst
692 y = (epst2 + c)* epst + d
693 yp = three*epst2 + c
694 IF(yp/=zero)epst = epst - y/yp
695 ENDIF
696 epst = epst + dav
697C-------------------
698C tension failure
699C-------------------
700 fail(i) = max(zero,min(one,
701 . (epsr2-epst)/(epsr2-epsr1) ))
702 ENDDO
703C-------------------
704C CRITERE
705C-------------------
706 IF(table(itable)%NDIM==2)THEN
707 DO i=1,nel
708 ipos2(i,1) = nint(uvar(i,2))
709 ipos2(i,2) = nint(uvar(i,3))
710C
711 xx2(i,1)=pla(i)
712 xx2(i,2)=epsp(i)*xfac
713 END DO
714C
715 CALL table_vinterp(table(itable),nel,nel,ipos2,xx2,yld,dydx)
716C
717 DO i=1,nel
718 yld(i) = yfac*yld(i)
719 yld(i) = fail(i)*yld(i)
720 h(i) = fail(i)*dydx(i)
721 uvar(i,2) = ipos2(i,1)
722 uvar(i,3) = ipos2(i,2)
723 END DO
724C-----
725 IF(fisokin/=zero)THEN
726 DO i=1,nel
727 ipos2(i,1) = 1
728C
729 xx2(i,1)=zero
730 xx2(i,2)=epsp(i)*xfac
731 END DO
732 CALL table_vinterp(table(itable),nel,nel,ipos2,xx2,yk,dydx)
733
734 DO i=1,nel
735C ECROUISSAGE CINEMATIQUE
736 yld(i) = (one-fisokin) * yld(i) +
737 . fisokin * fail(i) * yfac * yk(i)
738 yld(i) = max(yld(i),em20)
739 ENDDO
740 END IF
741c
742 ELSE
743C
744 DO i=1,nel
745 ipos3(i,1) = nint(uvar(i,2))
746 ipos3(i,2) = nint(uvar(i,3))
747 ipos3(i,3) = nint(uvar(i,4))
748C
749 xx3(i,1)=pla(i)
750 xx3(i,2)=epsp(i)*xfac
751 xx3(i,3)=temp(i)
752 END DO
753C
754 CALL table_vinterp(table(itable),nel,nel,ipos3,xx3,yld,dydx)
755C
756 DO i=1,nel
757 yld(i) = yfac*yld(i)
758 yld(i) = fail(i)*yld(i)
759 h(i) = fail(i)*dydx(i)
760 uvar(i,2) = ipos3(i,1)
761 uvar(i,3) = ipos3(i,2)
762 uvar(i,4) = ipos3(i,3)
763 END DO
764C-----
765 IF(fisokin/=zero)THEN
766 DO i=1,nel
767 ipos3(i,1) = 1
768C
769 xx3(i,1)=zero
770 xx3(i,2)=epsp(i)*xfac
771 xx3(i,3)=temp(i)
772 END DO
773 CALL table_vinterp(table(itable),nel,nel,ipos3,xx3,yk,dydx)
774
775 DO i=1,nel
776C ECROUISSAGE CINEMATIQUE
777 yld(i) = (one-fisokin) * yld(i) +
778 . fisokin * fail(i) * yfac * yk(i)
779 yld(i) = max(yld(i),em20)
780 ENDDO
781 END IF
782 END IF
783C-------------------
784C PROJECTION
785C-------------------
786 IF(ipla==0)THEN
787 DO i=1,nel
788 cri(i) =
789 . ff*(signyy(i) - signzz(i))**2 + gg*(signzz(i) - signxx(i))**2
790 . + hh*(signxx(i) - signyy(i))**2 + two*ll*signyz(i)**2
791 . + two*mm*signzx(i)**2 + two*nn*signxy(i)**2
792 cri(i) = sqrt(cri(i))
793 r = min(one,yld(i)/ max(cri(i),em20))
794c P = C1 * (RHO(I)/RHO0(I)- ONE)
795 p = c1 * amu(i)
796 signxx(i)=signxx(i)*r-p
797 signyy(i)=signyy(i)*r-p
798 signzz(i)=signzz(i)*r-p
799 signxy(i)=signxy(i)*r
800 signyz(i)=signyz(i)*r
801 signzx(i)=signzx(i)*r
802 pla(i)=pla(i) + (one - r)*cri(i)/max(g3+h(i),em20)
803 uvar(i,1)=pla(i)
804 dpla(i) = (one - r)*cri(i)/max(g3+h(i),em20)
805 ENDDO
806 ELSEIF(ipla==2)THEN
807 DO i=1,nel
808 cri(i) =
809 . ff*(signyy(i) - signzz(i))**2 + gg*(signzz(i) - signxx(i))**2
810 . + hh*(signxx(i) - signyy(i))**2 + two*ll*signyz(i)**2
811 . + two*mm*signzx(i)**2 + two*nn*signxy(i)**2
812 cri(i) = sqrt(cri(i))
813 r = min(one,yld(i)/ max(cri(i),em20))
814c P = C1 * (RHO(I)/RHO0(I)-ONE)
815 p = c1 * amu(i)
816 signxx(i)=signxx(i)*r-p
817 signyy(i)=signyy(i)*r-p
818 signzz(i)=signzz(i)*r-p
819 signxy(i)=signxy(i)*r
820 signyz(i)=signyz(i)*r
821 signzx(i)=signzx(i)*r
822 pla(i)=pla(i) + (one-r)*cri(i)/max(g3,em20)
823 uvar(i,1)=pla(i)
824 dpla(i) = (one-r)*cri(i)/max(g3,em20)
825 ENDDO
826 ELSEIF(ipla==1)THEN
827 DO i=1,nel
828 cri(i) =
829 . ff*(signyy(i) - signzz(i))**2 + gg*(signzz(i) - signxx(i))**2
830 . + hh*(signxx(i) - signyy(i))**2 + two*ll*signyz(i)**2
831 . + two*mm*signzx(i)**2 + two*nn*signxy(i)**2
832 cri(i) = sqrt(cri(i))
833 r = min(one,yld(i)/ max(cri(i),em20))
834C plastic strain increment.
835 dpla(i)=(one - r)*cri(i)/max(g3+h(i),em20)
836C actual yield stress.
837 yld(i)=max(yld(i)+(one - fisokin)*dpla(i)*h(i),zero)
838 r = min(one,yld(i)/ max(cri(i),em20))
839C
840c P = C1 * (RHO(I)/RHO0(I)-ONE)
841 p = c1 * amu(i)
842 signxx(i)=signxx(i)*r-p
843 signyy(i)=signyy(i)*r-p
844 signzz(i)=signzz(i)*r-p
845 signxy(i)=signxy(i)*r
846 signyz(i)=signyz(i)*r
847 signzx(i)=signzx(i)*r
848 pla(i)=pla(i) + dpla(i)
849 uvar(i,1)=pla(i)
850 ENDDO
851 ENDIF
852C------------------------------------------
853C ECROUISSAGE CINE
854C------------------------------------------
855 IF (fisokin/=zero) THEN
856 DO i=1,nel
857 p = c1 * (rho(i)/rho0(i)-one)
858 dsxx = sigexx(i) - signxx(i) - p
859 dsyy = sigeyy(i) - signyy(i) - p
860 dszz = sigezz(i) - signzz(i) - p
861 dsxy = sigexy(i) - signxy(i)
862 dsyz = sigeyz(i) - signyz(i)
863 dszx = sigezx(i) - signzx(i)
864C
865 hkin = two_third*fisokin*h(i)
866 alpha = hkin/(g2+hkin)
867 sigpxx = alpha*dsxx
868 sigpyy = alpha*dsyy
869 sigpzz = alpha*dszz
870 sigpxy = alpha*dsxy
871 sigpyz = alpha*dsyz
872 sigpzx = alpha*dszx
873C
874 uvar(i, 5) = uvar(i, 5) + sigpxx
875 uvar(i, 6) = uvar(i, 6) + sigpyy
876 uvar(i, 7) = uvar(i, 7) + sigpzz
877 uvar(i, 8) = uvar(i, 8) + sigpxy
878 uvar(i, 9) = uvar(i, 9) + sigpyz
879 uvar(i,10) = uvar(i,10) + sigpzx
880C
881 signxx(i) = signxx(i) + uvar(i, 5)
882 signyy(i) = signyy(i) + uvar(i, 6)
883 signzz(i) = signzz(i) + uvar(i, 7)
884 signxy(i) = signxy(i) + uvar(i, 8)
885 signyz(i) = signyz(i) + uvar(i, 9)
886 signzx(i) = signzx(i) + uvar(i,10)
887 ENDDO
888 ENDIF
889
890 IF (impl_s>0) THEN
891 DO i=1,nel
892 IF(dpla(i)>0) etse(i)= h(i)/g2
893 ENDDO
894c IF (IKT==4)
895c . CALL PUT_SIGE0(NEL,SIGEXX,SIGEYY,SIGEZZ,SIGEXY,
896c . SIGEYZ,SIGEZX,DPLA,G ,H )
897 ENDIF
898c
899 ENDIF ! opte==1 CE/=0
900C-----------------------------------------------
901C Update temperature due to plastic work
902C--------------------
903 IF (jthe < 0 .AND. jlag /= 0) THEN
904 DO i=1,nel
905 fheat(i) = fheat(i) + yld(i)*dpla(i)*volume(i)
906 ENDDO
907 ELSE IF (rhocp > zero) THEN
908 DO i=1,nel
909 temp(i) = temp(i) + yld(i)*dpla(i)*rhocp/max(em15,volume(i))
910 ENDDO
911 ENDIF
912C
913C-------------------
914 DO i=1,nel
915 IF(off(i)<em01) off(i)=zero
916 IF(off(i)<one) off(i)=off(i)*four_over_5
917 ENDDO
918C
919 nindx=0
920 DO i=1,nel
921 IF(pla(i)>epsmax.AND.off(i)==one) THEN
922 off(i)=four_over_5
923 nindx=nindx+1
924 indx(nindx)=i
925 ENDIF
926 ENDDO
927 IF(nindx>0.AND.imconv==1)THEN
928 DO j=1,nindx
929#include "lockon.inc"
930 WRITE(iout, 1000) ngl(indx(j))
931 WRITE(istdo,1100) ngl(indx(j)),tt
932#include "lockoff.inc"
933 ENDDO
934 ENDIF
935!
936 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
937 1100 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10,
938 . ' AT TIME :',g11.4)
939!
940!!-------------------------
941!! EQUIVALENT STRESS FOR OUTPUT - THERMAL HILL ORTHOTROPIC 3D CRITERION -
942!!-------------------------
943!! IF (ISEQ > 0) THEN
944!! DO I=1,NEL
945!! CRI(I) =
946!! . FF*(SIGNYY(I) - SIGNZZ(I))**2 + GG*(SIGNZZ(I) - SIGNXX(I))**2
947!! . + HH*(SIGNXX(I) - SIGNYY(I))**2 + TWO*LL*SIGNYZ(I)**2
948!! . + TWO*MM*SIGNZX(I)**2 + TWO*NN*SIGNXY(I)**2
949!! SEQ_OUTPUT(I) = SQRT(CRI(I))
950!! ENDDO
951!! ENDIF
952!---
953 RETURN
954 END
955C
alpha
#define alpha
Definition eval.h:35
interface_table_mod::table_vinterp
Definition table_mod.F:251
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
sigeps74
subroutine sigeps74(nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, fheat, jlag, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, ngl, ipt, ipm, mat, epsp, ipla, yld, pla, dpla, etse, jthe, temp, table, seq_output, amu, iseq)
Definition sigeps74.F:49
table_mod::ttable
Definition table_mod.F:125