OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps66.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!|| sigeps66 ../engine/source/materials/mat/mat066/sigeps66.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!|| vinter ../engine/source/tools/curve/vinter.F
31!||====================================================================
32 SUBROUTINE sigeps66(
33 1 NEL ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,
34 2 NPF ,TF ,TIME ,TIMESTEP,UPARAM ,
35 3 RHO0 ,RHO ,VOLUME ,EINT ,
36 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
37 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
39 7 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 9 SOUNDSP,VISCMAX,UVAR ,OFF ,NGL ,IPM ,
42 A MAT ,EPSP ,IPLA ,YLD ,PLA , ETSE ,
43 B AMU )
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48#include "comlock.inc"
49#include "param_c.inc"
50#include "com01_c.inc"
51#include "impl1_c.inc"
52C---------+---------+---+---+--------------------------------------------
53C VAR | SIZE |TYP| RW| DEFINITION
54C---------+---------+---+---+--------------------------------------------
55C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
56C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
57C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
58C---------+---------+---+---+--------------------------------------------
59C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
60C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
61C NPF | * | I | R | FUNCTION ARRAY
62C TF | * | F | R | FUNCTION ARRAY
63C---------+---------+---+---+--------------------------------------------
64C TIME | 1 | F | R | CURRENT TIME
65C TIMESTEP| 1 | F | R | CURRENT TIME STEP
66C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
67C RHO0 | NEL | F | R | INITIAL DENSITY
68C RHO | NEL | F | R | DENSITY
69C VOLUME | NEL | F | R | VOLUME
70C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
71C EPSPXX | NEL | F | R | STRAIN RATE XX
72C EPSPYY | NEL | F | R | STRAIN RATE YY
73C ... | | | |
74C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
75C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
76C ... | | | |
77C EPSXX | NEL | F | R | STRAIN XX
78C EPSYY | NEL | F | R | STRAIN YY
79C ... | | | |
80C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
81C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
82C ... | | | |
83C---------+---------+---+---+--------------------------------------------
84C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
85C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
86C ... | | | |
87C SIGVXX | NEL | F | W | VISCOUS STRESS XX
88C SIGVYY | NEL | F | W | VISCOUS STRESS YY
89C ... | | | |
90C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
91C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
92C---------+---------+---+---+--------------------------------------------
93C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
94C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
95C=======================================================================
96 INTEGER NEL, NUPARAM, NUVAR,IPLA,
97 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*)
98 my_real
99 . TIME,TIMESTEP,UPARAM(*),
100 . RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
101 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
102 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
103 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
104 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
105 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
106 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
107 . sigoxx(nel),sigoyy(nel),sigozz(nel),
108 . sigoxy(nel),sigoyz(nel),sigozx(nel),
109 . signxx(nel),signyy(nel),signzz(nel),
110 . signxy(nel),signyz(nel),signzx(nel),
111 . epsp(nel),etse(nel),amu(nel)
112C-----------------------------------------------
113C O U T P U T A r g u m e n t s
114C-----------------------------------------------
115 my_real
116 . soundsp(nel),viscmax(nel)
117C-----------------------------------------------
118C I N P U T O U T P U T A r g u m e n t s
119C-----------------------------------------------
120 my_real
121 . uvar(nel,nuvar), off(nel), yld(nel),pla(nel)
122C-----------------------------------------------
123C VARIABLES FOR FUNCTION INTERPOLATION
124C-----------------------------------------------
125 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
126 my_real
127 . finter,finter2, tf(*)
128 EXTERNAL finter,finter2
129C-----------------------------------------------
130C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
131C Y : y = f(x)
132C X : x
133C DYDX : f'(x) = dy/dx
134C IFUNC(J): FUNCTION INDEX
135C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
136C NPF,TF : FUNCTION PARAMETER
137C-----------------------------------------------
138C L o c a l V a r i a b l e s
139C-----------------------------------------------
140 INTEGER I,J,J1,J2,IYLDC,IYLDT,IRATE,NCC,NCT,NFUNC,VP,IADBUF
141 INTEGER IPOS1(NEL),IAD1(NEL), ILEN1(NEL),
142 . IPOS2(NEL),IAD2(NEL), ILEN2(NEL),IFUNC(MFUNC),
143 . JJC(NEL),JJT(NEL)
144 my_real ET, EC, C1T, GT,NU, CP, PC,PT,RPCT,EPSP0,
145 . DAV,DPLA,FAC,EPD,VM,
146 . R,FISOKIN,YRATE,HKIN ,SIGPXX,SIGPYY,
147 . SIGPZZ,SIGPXY,SIGPYZ,SIGPZX,DSXX,DSYY,DSZZ,
148 . DSXY,DSYZ,DSZX,DF,SIGY,DTINV
149 my_real ,DIMENSION(MFUNC):: YFAC,RATE0
150 my_real ,DIMENSION(NEL) ::
151 . yc ,yt ,h ,dydx2 ,dydx1 ,dpla1 ,
152 . p0 ,p ,hc ,ht ,rate ,alpha,
153 . y1 ,y2 ,e ,c1 ,g , g2 ,g3,
154 . sigexx ,sigeyy ,sigezz ,sigexy ,sigeyz ,sigezx
155C=======================================================================
156C USER VARIABLES INITIALIZATION
157C-----------------------------------------------
158 DO i=1,nel
159 etse(i) = one
160 dpla1(i) = zero
161 ENDDO
162C
163 iadbuf = ipm(7,mat(1))
164 irate = nint(uparam(iadbuf ))
165 et = uparam(iadbuf +1)
166 gt = uparam(iadbuf +4)
167
168 nu = uparam(iadbuf +5)
169 pc = uparam(iadbuf +6)
170 pt = uparam(iadbuf +7)
171 epsp0 = uparam(iadbuf +8)
172 cp = uparam(iadbuf +9)
173 ncc = nint(uparam(iadbuf + 10))
174 nct = nint(uparam(iadbuf + 11))
175 fisokin = uparam(iadbuf + 12)
176C
177C
178 nfunc = ipm(10,mat(1))
179 DO i=1,nfunc
180 ifunc(i) = ipm(10+i,mat(1))
181 yfac(i) = uparam(iadbuf + 12 + i )
182 rate0(i) = uparam(iadbuf + 12 + nfunc + i )
183 ENDDO
184C
185 sigy = uparam(iadbuf + 13 + 2*mfunc)
186 vp = nint(uparam(iadbuf + 14 + 2*mfunc))
187 ec = uparam(iadbuf + 15 + 2*mfunc)
188 rpct = uparam(iadbuf + 16 + 2*mfunc)
189C
190 IF (isigi==0) THEN
191 IF(time==zero)THEN
192 DO i=1,nel
193 uvar(i,1)=zero
194 uvar(i,2)=zero
195 uvar(i,3)=zero
196 uvar(i,4)=zero
197 uvar(i,5)=zero
198 uvar(i,6)=zero
199 uvar(i,7)=zero
200 DO j=1,nfunc
201 uvar(i,j+7)=zero
202 ENDDO
203 ENDDO
204 ENDIF
205 ENDIF
206c------------------------------------------------
207c------------------------------------------------
208 e(1:nel) = et
209 g(1:nel) = gt
210 c1t = et/three/(one - two*nu)
211 IF(ec > zero)THEN
212 DO i=1,nel
213c P(I) = C1 * (RHO(I)/RHO0(I)- ONE)
214 p(i) = c1t * amu(i)
215 c1t = max(c1t,ec/three/(one - two*nu))
216 IF(pc == zero .and. pt == zero .AND. abs(p(i)) < em10) THEN
217 e(i) = ec
218 ELSEIF(p(i) <= - rpct * pt) THEN
219 e(i) = et
220 ELSEIF(p(i) >= rpct *pc) THEN
221 e(i) = ec
222 ELSE
223 fac = rpct *(pc + pt)
224 fac = (rpct * pc - p(i))/fac
225 e(i) = fac*et + (one -fac)*ec
226 ENDIF
227 ENDDO
228 ENDIF
229 DO i=1,nel
230 g(i) = half*e(i)/( one + nu)
231 g2(i) = two*g(i)
232 g3(i) = three*g(i)
233 c1(i) = e(i)/three/(one - two*nu)
234 ENDDO
235
236C------------------------------------------
237C estiamte stresse CINE
238C------------------------------------------
239 IF (fisokin > zero ) THEN
240 DO i=1,nel
241 sigoxx(i) = sigoxx(i) - uvar(i, 2)
242 sigoyy(i) = sigoyy(i) - uvar(i, 3)
243 sigozz(i) = sigozz(i) - uvar(i, 4)
244 sigoxy(i) = sigoxy(i) - uvar(i, 5)
245 sigoyz(i) = sigoyz(i) - uvar(i, 6)
246 sigozx(i) = sigozx(i) - uvar(i, 7)
247 ENDDO
248 ENDIF
249C
250 DO i=1,nel
251 pla(i) = uvar(i,1)
252 p0(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third
253 dav = (depsxx(i)+depsyy(i)+depszz(i))*third
254 signxx(i)=sigoxx(i) + p0(i) + g2(i)*(depsxx(i)-dav)
255 signyy(i)=sigoyy(i) + p0(i) + g2(i)*(depsyy(i)-dav)
256 signzz(i)=sigozz(i) + p0(i) + g2(i)*(depszz(i)-dav)
257 signxy(i)=sigoxy(i) + g(i) *depsxy(i)
258 signyz(i)=sigoyz(i) + g(i) *depsyz(i)
259 signzx(i)=sigozx(i) + g(i) *depszx(i)
260C
261 soundsp(i) = sqrt((c1t + four*g(i)/three)/rho0(i))
262 viscmax(i) = zero
263C
264 sigexx(i) = signxx(i)
265 sigeyy(i) = signyy(i)
266 sigezz(i) = signzz(i)
267 sigexy(i) = signxy(i)
268 sigeyz(i) = signyz(i)
269 sigezx(i) = signzx(i)
270C
271 jjc(i) = 1
272 jjt(i) = ncc + 1
273 ENDDO
274C-------------------
275C CRITERE
276C-------------------
277 IF (irate <= 3) THEN
278 DO i=1,nel
279 ipos1(i) = nint(uvar(i,8))
280 iad1(i) = npf(ifunc(1)) / 2+ 1
281 ilen1(i) = npf(ifunc(1)+ 1) / 2 - iad1(i)-ipos1(i)
282 ipos2(i) = nint(uvar(i,9))
283 iad2(i) = npf(ifunc(2)) / 2 + 1
284 ilen2(i) = npf(ifunc(2) + 1) / 2 - iad2(i)-ipos2(i)
285C
286 uvar(i,8) = ipos1(i)
287 uvar(i,9) = ipos2(i)
288 END DO
289C
290 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,yc)
291 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,yt)
292C
293 IF(fisokin == zero) THEN
294 DO i=1,nel
295 yc(i)=yc(i)*yfac(1)
296 yt(i)=yt(i)*yfac(2)
297 hc(i)=dydx1(i)*yfac(1)
298 ht(i)=dydx2(i)*yfac(2)
299 ENDDO
300 ELSEIF(fisokin == one ) THEN
301 DO i=1,nel
302 hc(i)=dydx1(i)*yfac(1)
303 ht(i)=dydx2(i)*yfac(2)
304 yc(i)=tf(npf(ifunc(1)) + 1)
305 yt(i)=tf(npf(ifunc(2)) + 1)
306 yc(i)=yc(i)*yfac(1)
307 yt(i)=yc(i)*yfac(2)
308 yc(i) = max(em20, yc(i))
309 yt(i) = max(em20, yt(i))
310 ENDDO
311 ELSE
312 DO i=1,nel
313 yc(i)=yc(i)*yfac(1)
314 yt(i)=yt(i)*yfac(2)
315 yc(i) = max(yc(i),em20)
316 yt(i) = max(yt(i),em20)
317 hc(i) = dydx1(i)*yfac(1)
318 ht(i) = dydx2(i)*yfac(2)
319C ECROUISSAGE CINEMATIQUE
320 y1(i)=yfac(1)*tf(npf(ifunc(1))+1)
321 y2(i)=yfac(2)*tf(npf(ifunc(2))+1)
322 yc(i) = (one - fisokin) * yc(i) + fisokin * y1(i)
323 yt(i) = (one - fisokin) * yt(i) + fisokin * y2(i)
324 ENDDO
325 ENDIF
326 ELSE ! IRATE =4 => multiples curves (tracation & compression)
327C----------------------------------------------------------------------
328C compression
329C------------------------------------------------------------------------
330 DO j = 2,ncc - 1
331 DO i=1,nel
332 IF(epsp(i) >= rate0(j))jjc(i) = j
333 ENDDO
334 ENDDO
335 DO i=1,nel
336 fac=rate0(jjc(i))
337 rate(i)=(epsp(i) - fac)/(rate0(jjc(i)+1 ) - fac)
338 ENDDO
339 DO i=1,nel
340 j1 = jjc(i)
341 j2 = j1+1
342 ipos1(i) = nint(uvar(i,7+j1 ))
343 iad1(i) = npf(ifunc(j1)) / 2 + 1
344 ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i)-ipos1(i)
345 ipos2(i) = nint(uvar(i,7+j2))
346 iad2(i) = npf(ifunc(j2)) / 2 + 1
347 ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i)-ipos2(i)
348C
349 uvar(i,7+j1) = ipos1(i)
350 uvar(i,7+j2) = ipos2(i)
351 END DO
352 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
353 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
354C
355 IF (fisokin == zero) THEN
356 DO i=1,nel
357 j1 = jjc(i)
358 j2 = j1+1
359 y1(i)=y1(i)*yfac(j1)
360 y2(i)=y2(i)*yfac(j2)
361 fac = rate(i)
362 yc(i) =(y1(i) + fac*(y2(i)-y1(i)))
363 yc(i) = max(yc(i),em20)
364 dydx1(i)=dydx1(i)*yfac(j1)
365 dydx2(i)=dydx2(i)*yfac(j2)
366 hc(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
367 ENDDO
368 ELSEIF (fisokin == one ) THEN
369 DO i=1,nel
370 j1 = jjc(i)
371 j2 = j1+1
372 fac = rate(i)
373 dydx1(i)=dydx1(i)*yfac(j1)
374 dydx2(i)=dydx2(i)*yfac(j2)
375 hc(i) =(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
376C ECROUISSAGE CINEMATIQUE
377 y1(i)=yfac(j1)*tf(npf(ifunc(j1)) + 1)
378 y2(i)=yfac(j2)*tf(npf(ifunc(j2)) + 1)
379 yc(i) = (y1(i) + fac*(y2(i)-y1(i)))
380 yc(i) = max(em20,yc(i))
381 ENDDO
382 ELSE
383 DO i=1,nel
384 j1 = jjc(i)
385 j2 = j1 + 1
386 y1(i)=y1(i)*yfac(j1)
387 y2(i)=y2(i)*yfac(j2)
388 fac = rate(i)
389 yc(i) = (y1(i) + fac*(y2(i)-y1(i)))
390 yc(i) = max(yc(i),em20)
391 dydx1(i)=dydx1(i)*yfac(j1)
392 dydx2(i)=dydx2(i)*yfac(j2)
393 hc(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
394C ECROUISSAGE CINEMATIQUE
395 y1(i)=yfac(j1)*tf(npf(ifunc(j1))+1)
396 y2(i)=yfac(j2)*tf(npf(ifunc(j2))+1)
397 yc(i) = (one - fisokin) * yc(i) +
398 . fisokin * (y1(i) + fac*(y2(i)-y1(i)))
399 ENDDO
400 ENDIF
401C----------------------------------------------------------------------
402C traction
403C------------------------------------------------------------------------
404 DO j = 2,nct - 1
405 DO i=1,nel
406 IF (epsp(i)>=rate0(ncc + j))jjt(i)=ncc+j
407 ENDDO
408 ENDDO
409 DO i=1,nel
410 fac=rate0(jjt(i))
411 rate(i)=(epsp(i) - fac)/( rate0(jjt(i)+1)- fac)
412 ENDDO
413 DO i=1,nel
414 j1 = jjt(i)
415 j2 = j1+1
416 ipos1(i) = nint(uvar(i,7+j1))
417 iad1(i) = npf(ifunc(j1)) / 2 + 1
418 ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i)-ipos1(i)
419 ipos2(i) = nint(uvar(i,7+j2))
420 iad2(i) = npf(ifunc(j2)) / 2 + 1
421 ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i)-ipos2(i)
422C
423 uvar(i,7+j1) = ipos1(i)
424 uvar(i,7+j2) = ipos2(i)
425 END DO
426 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
427 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
428C
429 IF (fisokin == zero) THEN
430 DO i=1,nel
431 j1 = jjt(i)
432 j2 = j1+1
433 y1(i)=y1(i)*yfac(j1)
434 y2(i)=y2(i)*yfac(j2)
435 fac = rate(i)
436 yt(i) =(y1(i) + fac*(y2(i)-y1(i)))
437 yt(i) = max(yt(i),em20)
438 dydx1(i)=dydx1(i)*yfac(j1)
439 dydx2(i)=dydx2(i)*yfac(j2)
440 ht(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
441 ENDDO
442 ELSEIF (fisokin == one ) THEN
443 DO i=1,nel
444 j1 = jjt(i)
445 j2 = j1+1
446 fac = rate(i)
447 dydx1(i)=dydx1(i)*yfac(j1)
448 dydx2(i)=dydx2(i)*yfac(j2)
449 ht(i) =(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
450C ECROUISSAGE CINEMATIQUE
451 y1(i)=yfac(j1)*tf(npf(ifunc(j1)) + 1)
452 y2(i)=yfac(j2)*tf(npf(ifunc(j2)) + 1)
453 yt(i) = (y1(i) + fac*(y2(i)-y1(i)))
454 yt(i) = max(em20,yt(i))
455 ENDDO
456 ELSE
457 DO i=1,nel
458 j1 = jjt(i)
459 j2 = j1 + 1
460 y1(i)=y1(i)*yfac(j1)
461 y2(i)=y2(i)*yfac(j2)
462 fac = rate(i)
463 yt(i) = (y1(i) + fac*(y2(i)-y1(i)))
464 yt(i) = max(yt(i),em20)
465 dydx1(i)=dydx1(i)*yfac(j1)
466 dydx2(i)=dydx2(i)*yfac(j2)
467 ht(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
468C ECROUISSAGE CINEMATIQUE
469 y1(i)=yfac(j1)*tf(npf(ifunc(j1))+1)
470 y2(i)=yfac(j2)*tf(npf(ifunc(j2))+1)
471 yt(i) = (one - fisokin) * yt(i) +
472 . fisokin * (y1(i) + fac*(y2(i)-y1(i)))
473 ENDDO
474 ENDIF
475 ENDIF
476C
477C strain effect on compression and traction
478 IF(irate == 3) THEN
479 DO i=1,nel
480C compression
481 yrate = yfac(3)*finter(ifunc(3),epsp(i),npf,tf,df)
482 yc(i) = yrate*yc(i)
483 hc(i) = hc(i)*yrate
484c traction
485 yrate = yfac(4)*finter(ifunc(4),epsp(i),npf,tf,df)
486 yt(i) = yrate*yt(i)
487 ht(i) = ht(i)*yrate
488 ENDDO
489 ENDIF
490C interpolation entre (tracation & compression
491 DO i=1,nel
492c P(I) = C1 * (RHO(I)/RHO0(I)- ONE)
493 p(i) = c1(i) * amu(i)
494 IF(pc == zero .AND. pt == zero .AND. abs(p(i)) < em10) THEN
495 yld(i) = max(yc(i), em20)
496 h(i) = hc(i)
497 ELSEIF(p(i) <= -pt) THEN
498 yld(i) = max(yt(i),em20)
499 h(i) = ht(i)
500 ELSEIF(p(i) >= pc) THEN
501 yld(i) = max(yc(i), em20)
502 h(i) = hc(i)
503 ELSE
504 fac = pc + pt
505 fac = (pc - p(i))/fac
506 yld(i) = fac*yt(i) + (one -fac)*yc(i)
507 yld(i) = max(em20,yld(i))
508 h(i) = fac*ht(i) + (one -fac)*hc(i)
509 ENDIF
510 ENDDO
511C
512C strain rate effect
513c
514C
515c analytical strain rate effect
516c (1 + (epsp/epsp0)**(1/cp)
517 IF (vp == 0) THEN
518 IF(irate == 1) THEN
519 DO i=1,nel
520 epd = max(em20,epsp(i)/epsp0)
521 yrate = one + exp(cp*log(epd))
522 yld(i) = yld(i)*yrate
523 h(i) = h(i)*yrate
524 ENDDO
525C 1 + cp*ln(epep/epsp0)
526 ELSEIF(irate == 2) THEN
527 DO i=1,nel
528 epd = max(em20,epsp(i)/epsp0)
529 yrate = one + cp*log(epd)
530 yld(i) = yld(i)*yrate
531 h(i) = h(i)*yrate
532 ENDDO
533 ENDIF
534 ENDIF
535C-------------------
536C PROJECTION
537C-------------------
538 IF(ipla==0)THEN
539C
540 IF(vp > 0 ) THEN
541 dtinv = timestep/max(em20, timestep**2)
542 DO i=1,nel
543 vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
544 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
545 vm =sqrt(three*vm)
546 r = min(one,yld(i)/ max(vm,em20))
547 epd = (one - r)*vm/max(g3(i)+h(i),em20)
548 epd = epd*dtinv
549 epd = max(em20,epd/epsp0)
550 yrate = one + exp(cp*log(epd))
551 IF(sigy == zero) THEN
552 yld(i) = yld(i)*yrate
553 h(i) = h(i)*yrate
554 ELSE
555 yld(i) = yld(i) + sigy*(yrate - one)
556 ENDIF
557 ENDDO
558 ENDIF
559C
560 DO i=1,nel
561 vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
562 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
563 vm =sqrt(three*vm)
564 r = min(one,yld(i)/ max(vm,em20))
565 signxx(i)=signxx(i)*r-p(i)
566 signyy(i)=signyy(i)*r-p(i)
567 signzz(i)=signzz(i)*r-p(i)
568 signxy(i)=signxy(i)*r
569 signyz(i)=signyz(i)*r
570 signzx(i)=signzx(i)*r
571 pla(i)=pla(i) + (one - r)*vm/max(g3(i) + h(i),em20)
572 uvar(i,1)=pla(i)
573 dpla1(i) = (one - r)*vm/max(g3(i)+h(i),em20)
574 ENDDO
575 ELSEIF(ipla==2)THEN
576C
577 IF(vp > 0 ) THEN
578 dtinv = timestep/max(em20, timestep**2)
579 DO i=1,nel
580 vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
581 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
582 vm =sqrt(three*vm)
583 r = min(one,yld(i)/ max(vm,em20))
584 epd = (one-r)*vm/max(g3(i),em20)
585 epd = epd*dtinv
586 epd = max(em20,epd/epsp0)
587 yrate = one + exp(cp*log(epd))
588 IF(sigy == zero) THEN
589 yld(i) = yld(i)*yrate
590 h(i) = h(i)*yrate
591 ELSE
592 yld(i) = yld(i) + sigy*(yrate - one)
593 ENDIF
594 ENDDO
595 ENDIF
596 DO i=1,nel
597 vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
598 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
599 vm =sqrt(three*vm)
600 r = min(one,yld(i)/ max(vm,em20))
601 signxx(i)=signxx(i)*r-p(i)
602 signyy(i)=signyy(i)*r-p(i)
603 signzz(i)=signzz(i)*r-p(i)
604 signxy(i)=signxy(i)*r
605 signyz(i)=signyz(i)*r
606 signzx(i)=signzx(i)*r
607 pla(i)=pla(i) + (one-r)*vm/max(g3(i),em20)
608 uvar(i,1)=pla(i)
609 dpla1(i) = (one-r)*vm/max(g3(i),em20)
610 ENDDO
611 ELSEIF (ipla == 1) THEN
612CC Viscoplatic
613 IF (vp > 0 ) THEN
614 dtinv = timestep/max(em20, timestep**2)
615 DO i=1,nel
616 vm = half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
617 . + signxy(i)**2+signyz(i)**2+signzx(i)**2
618 vm = sqrt(three*vm)
619 r = min(one,yld(i)/ max(vm,em20))
620C
621 epd = (one - r)*vm/max(g3(i)+h(i),em20)
622 epd = epd*dtinv
623 epd = max(em20,epd/epsp0)
624 yrate = one + exp(cp*log(epd))
625 IF (sigy == zero) THEN
626 yld(i) = yld(i)*yrate
627 h(i) = h(i)*yrate
628 ELSE
629 yld(i) = yld(i) + sigy*(yrate - one)
630 ENDIF
631 ENDDO
632 ENDIF
633C
634 DO i=1,nel
635 vm = three*(half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
636 . + signxy(i)**2+signyz(i)**2+signzx(i)**2)
637 IF (vm > yld(i)*yld(i)) THEN
638 vm = sqrt(vm)
639 r = yld(i) / vm
640c plastic strain increment
641 dpla = (one - r)*vm / max(g3(i)+h(i),em20)
642c actualized yield stress with kinematic hardening
643 yld(i) = max(yld(i) + (one-fisokin)*dpla*h(i), zero)
644 r = min(one, yld(i) / vm)
645C
646 signxx(i) = signxx(i)*r
647 signyy(i) = signyy(i)*r
648 signzz(i) = signzz(i)*r
649 signxy(i) = signxy(i)*r
650 signyz(i) = signyz(i)*r
651 signzx(i) = signzx(i)*r
652 pla(i) = pla(i) + dpla
653 uvar(i,1)= pla(i)
654 dpla1(i) = dpla
655 ENDIF
656 ENDDO
657c
658 DO i=1,nel
659 signxx(i) = signxx(i) - p(i)
660 signyy(i) = signyy(i) - p(i)
661 signzz(i) = signzz(i) - p(i)
662 ENDDO
663c
664 ENDIF
665C+
666C------------------------------------------
667C ECROUISSAGE CINE
668C------------------------------------------
669C
670 IF(fisokin > zero) THEN
671 DO i=1,nel
672 dsxx = sigexx(i) - signxx(i) - p(i)
673 dsyy = sigeyy(i) - signyy(i) - p(i)
674 dszz = sigezz(i) - signzz(i) - p(i)
675 dsxy = sigexy(i) - signxy(i)
676 dsyz = sigeyz(i) - signyz(i)
677 dszx = sigezx(i) - signzx(i)
678CC
679 hkin = two_third*fisokin*h(i)
680 alpha(i) = hkin/(g2(i)+hkin)
681 sigpxx = alpha(i)*dsxx
682 sigpyy = alpha(i)*dsyy
683 sigpzz = alpha(i)*dszz
684 sigpxy = alpha(i)*dsxy
685 sigpyz = alpha(i)*dsyz
686 sigpzx = alpha(i)*dszx
687C
688 uvar(i, 2) = uvar(i, 2) + sigpxx
689 uvar(i, 3) = uvar(i, 3) + sigpyy
690 uvar(i, 4) = uvar(i, 4) + sigpzz
691 uvar(i, 5) = uvar(i, 5) + sigpxy
692 uvar(i, 6) = uvar(i, 6) + sigpyz
693 uvar(i, 7) = uvar(i, 7) + sigpzx
694C
695 signxx(i) = signxx(i) + uvar(i, 2)
696 signyy(i) = signyy(i) + uvar(i, 3)
697 signzz(i) = signzz(i) + uvar(i, 4)
698 signxy(i) = signxy(i) + uvar(i, 5)
699 signyz(i) = signyz(i) + uvar(i, 6)
700 signzx(i) = signzx(i) + uvar(i, 7)
701 ENDDO
702 ENDIF
703C
704 IF (impl_s > zero) THEN
705 DO i=1,nel
706 IF(dpla1(i) > zero) etse(i)= h(i)/g2(i)
707 ENDDO
708 ENDIF
709C-----------
710 RETURN
711 END
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
sigeps66
subroutine sigeps66(nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, 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, soundsp, viscmax, uvar, off, ngl, ipm, mat, epsp, ipla, yld, pla, etse, amu)
Definition sigeps66.F:44
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72