OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps66c.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!|| sigeps66c ../engine/source/materials/mat/mat066/sigeps66c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.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!||--- uses -----------------------------------------------------
32!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
33!||====================================================================
34 SUBROUTINE sigeps66c(
35 1 NEL ,NUPARAM ,NUVAR ,MFUNC ,KFUNC ,
36 2 NPF ,NPT0 ,IPT ,IFLAG ,
37 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,
38 3 AREA ,EINT ,THKLY ,ISRATE ,ASRATE ,
39 4 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
40 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
41 6 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
42 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
43 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
44 9 SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
45 A SOUNDSP,VISCMAX ,THK ,PLA ,UVAR ,
46 B OFF ,NGL ,IPM ,MAT ,ETSE ,
47 C GS ,YLD ,EPSD_PG ,EPSP ,INLOC ,
48 D DPLANL ,MAT_PARAM,NUVARV ,UVARV ,LOFF )
49C-----------------------------------------------
50C M o d u l e s
51C-----------------------------------------------
52 USE matparam_def_mod
53C-----------------------------------------------
54C I m p l i c i t T y p e s
55C-----------------------------------------------
56#include "implicit_f.inc"
57C---------+---------+---+---+--------------------------------------------
58C VAR | SIZE |TYP| RW| DEFINITION
59C---------+---------+---+---+--------------------------------------------
60C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
61C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
62C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
63C---------+---------+---+---+--------------------------------------------
64C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
65C IFUNC | NFUNC | I | R | FUNCTION INDEX
66C NPF | * | I | R | FUNCTION ARRAY
67C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
68C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
69C IFLAG | * | I | R | GEOMETRICAL FLAGS
70C TF | * | F | R | FUNCTION ARRAY
71C---------+---------+---+---+--------------------------------------------
72C TIME | 1 | F | R | CURRENT TIME
73C TIMESTEP| 1 | F | R | CURRENT TIME STEP
74C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
75C RHO0 | NEL | F | R | INITIAL DENSITY
76C AREA | NEL | F | R | AREA
77C EINT | 2*NEL | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
78C THKLY | NEL | F | R | LAYER THICKNESS
79C EPSPXX | NEL | F | R | STRAIN RATE XX
80C EPSPYY | NEL | F | R | STRAIN RATE YY
81C ... | | | |
82C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
83C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
84C ... | | | |
85C EPSXX | NEL | F | R | STRAIN XX
86C EPSYY | NEL | F | R | STRAIN YY
87C ... | | | |
88C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
89C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
90C ... | | | |
91C---------+---------+---+---+--------------------------------------------
92C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
93C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
94C ... | | | |
95C SIGVXX | NEL | F | W | VISCOUS STRESS XX
96C SIGVYY | NEL | F | W | VISCOUS STRESS YY
97C ... | | | |
98C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
99C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
100C---------+---------+---+---+--------------------------------------------
101C THK | NEL | F |R/W| THICKNESS
102C PLA | NEL | F |R/W| PLASTIC STRAIN
103C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
104C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
105C---------+---------+---+---+--------------------------------------------
106C C o m m o n B l o c k s
107C-----------------------------------------------
108#include "param_c.inc"
109#include "com01_c.inc"
110C-----------------------------------------------
111C I N P U T A r g u m e n t s
112C-----------------------------------------------
113 INTEGER ,INTENT(IN) :: NUVARV
114 INTEGER NEL, NUPARAM, NUVAR, NPT0, IPT,IFLAG(*),
115 . NGL(NEL),MAT(NEL),ISRATE,IPM(NPROPMI,*),INLOC
116 my_real ,INTENT(IN) :: ASRATE
117 my_real ,DIMENSION(NEL) ,INTENT(IN) :: EPSD_PG !< global element strain rate in Gauss pt
118 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: EPSP !< local strain rate used in equations
119 my_real
120 . TIME,TIMESTEP(NEL),UPARAM(*),
121 . AREA(NEL),RHO0(NEL),EINT(NEL,2),
122 . THKLY(NEL),PLA(NEL),
123 . EPSPXX(NEL),EPSPYY(NEL),
124 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
125 . DEPSXX(NEL),DEPSYY(NEL),
126 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
127 . epsxx(nel) ,epsyy(nel) ,
128 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
129 . sigoxx(nel),sigoyy(nel),
130 . sigoxy(nel),sigoyz(nel),sigozx(nel),
131 . gs(*),dplanl(nel)
132 TYPE(matparam_struct_) ,INTENT(IN) :: MAT_PARAM
133 my_real, DIMENSION(NEL), INTENT(IN) :: loff
134C-----------------------------------------------
135C O U T P U T A r g u m e n t s
136C-----------------------------------------------
137 my_real
138 . signxx(nel),signyy(nel),
139 . signxy(nel),signyz(nel),signzx(nel),
140 . sigvxx(nel),sigvyy(nel),
141 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
142 . soundsp(nel),viscmax(nel),etse(nel),
143 . dpla_i(nel)
144C-----------------------------------------------
145C I N P U T O U T P U T A r g u m e n t s
146C-----------------------------------------------
147 my_real ,DIMENSION(NEL,NUVARV) ,INTENT(INOUT) :: uvarv
148 my_real ,INTENT(INOUT) :: uvar(nel,nuvar), off(nel),thk(nel),yld(nel)
149C-----------------------------------------------
150C VARIABLES FOR FUNCTION INTERPOLATION
151C-----------------------------------------------
152 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
153 my_real
154 . FINTER ,TF(*),FINTER2
155 EXTERNAL FINTER,FINTER2
156C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
157C Y : y = f(x)
158C X : x
159C DYDX : f'(x) = dy/dx
160C IFUNC(J): FUNCTION INDEX
161C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
162C NPF,TF : FUNCTION PARAMETER
163C-----------------------------------------------
164C L o c a l V a r i a b l e s
165C-----------------------------------------------
166 INTEGER I,J,N,NINDX,NMAX,IRATE,NCC,NCT,NFUNC,
167 . j1,j2,nprony,iadbuf,
168 . iad1(nel),ipos1(nel),ilen1(nel),
169 . iad2(nel),ipos2(nel),ilen2(nel),
170 . index(nel),jjc(nel),jjt(nel),ifunc(mfunc)
171 my_real ec, e1t,a11t,a21t,g1t,
172 . yfacc,yfact,cp,edp,fisokin,
173 . nu,pc,pt,rpct, epspo,nnu11,epsp0,g31,dsxx,dsyy,dsxy,
174 . dexx,deyy,dexy,sigpxx,sigpyy,sigpxy,kv,vp,sigy,
175 . fac,epd,nu31,r,umr,dezz,s2,s3,nu11,s1,nu21,
176 . nu110,nu210,p2,q2,f,df,a,yld_i,yrate,hkin,alpha,
177 . dtinv
178 my_real
179 . rate0(mfunc),gv(100),beta(100),yfac(mfunc)
180 my_real ,DIMENSION(NEL) ::
181 . e, a11 , a21, g,g3,c1, sigexx ,sigeyy ,sigexy ,
182 . dydx1 ,dydx2 ,svm ,
183 . yc ,yt ,p ,h ,
184 . hk ,nnu1 ,hi ,dpla_j ,
185 . aa ,bb ,dr ,pp ,qq ,
186 . svm2 ,yld2 ,y1 ,
187 . y2 ,ht ,rate ,hc ,epspzz
188
189 DATA nmax/3/
190C-----------------------------------------------
191C USER VARIABLES INITIALIZATION
192C-----------------------------------------------
193 iadbuf = ipm(7,mat(1))
194C
195 irate = nint(uparam(iadbuf ))
196 e1t = uparam(iadbuf +1)
197 a11t = uparam(iadbuf +2)
198 a21t = uparam(iadbuf +3)
199 g1t = uparam(iadbuf +4)
200 nu = uparam(iadbuf +5)
201 pc = uparam(iadbuf +6)
202 pt = uparam(iadbuf +7)
203 epsp0 = uparam(iadbuf +8)
204 cp = uparam(iadbuf +9)
205 ncc = nint(uparam(iadbuf +10))
206 nct = nint(uparam(iadbuf +11))
207 fisokin = uparam(iadbuf + 12)
208
209 nfunc = ipm(10,mat(1))
210C
211 DO i=1,nfunc
212 ifunc(i) = ipm(10+i,mat(1))
213 yfac(i) = uparam(iadbuf +12+i)
214 rate0(i) = uparam(iadbuf +12+nfunc + i)
215 ENDDO
216C
217 sigy = uparam(iadbuf + 13 + 2*mfunc)
218 vp = uparam(iadbuf + 14 + 2*mfunc)
219c
220 ec = uparam(iadbuf + 15 + 2*mfunc)
221 rpct = uparam(iadbuf + 16 + 2*mfunc)
222C
223 IF(vp /=0 .and .iflag(1) /= 1) vp = 0
224C
225 nnu11 = nu / (one - nu)
226C
227 IF (isigi==0) THEN
228 IF(time==zero)THEN
229 DO i=1,nel
230 uvar(i,1)=zero
231 uvar(i,2)=zero
232 uvar(i,3)=zero
233 uvar(i,4)=zero
234 DO j=1,nfunc
235 uvar(i,j+4)=zero
236 ENDDO
237 ENDDO
238 ENDIF
239 ENDIF
240C------------------------------------------
241c estimation if compression mode
242C------------------------------------------
243 DO i=1,nel
244 signxx(i)=sigoxx(i) - uvar(i,2) + a11t*depsxx(i)+a21t*depsyy(i)
245 signyy(i)=sigoyy(i) - uvar(i,3) + a21t*depsxx(i)+a11t*depsyy(i)
246 p(i) = -third*(signxx(i) + signyy(i) )
247 ENDDO
248 e(1:nel) = e1t
249 g(1:nel) = g1t
250 IF(ec > zero)THEN
251 DO i=1,nel
252 IF(pc == zero .and. pt == zero .AND. abs(p(i)) < em10) THEN
253 e(i) = ec
254 ELSEIF(p(i) <= - rpct * pt) THEN
255 e(i) = e1t
256 ELSEIF(p(i) >= rpct *pc) THEN
257 e(i) = ec
258 ELSE
259 fac = rpct *(pc + pt)
260 fac = (rpct * pc - p(i))/fac
261 e(i) = fac*e1t + (one -fac)*ec
262 ENDIF
263 ENDDO
264 ENDIF
265 DO i=1,nel
266 g(i) = half*e(i)/( one + nu)
267 a11(i) = e(i)/(one - nu*nu)
268 a21(i) = nu * a11(i)
269 g3(i) = three*g(i)
270 ENDDO
271
272
273C------------------------------------------
274C estiamte stress
275C------------------------------------------
276 DO i=1,nel
277C
278C en cas d' crouissage cinematique
279C
280
281 signxx(i)=sigoxx(i) - uvar(i,2) + a11(i)*depsxx(i)+a21(i)*depsyy(i)
282 signyy(i)=sigoyy(i) - uvar(i,3) + a21(i)*depsxx(i)+a11(i)*depsyy(i)
283 signxy(i)=sigoxy(i) - uvar(i,4) + g(i) *depsxy(i)
284 signyz(i)=sigoyz(i) + gs(i) *depsyz(i)
285 signzx(i)=sigozx(i) + gs(i) *depszx(i)
286C
287 sigexx(i) = signxx(i)
288 sigeyy(i) = signyy(i)
289 sigexy(i) = signxy(i)
290
291 p(i) = -third*(signxx(i) + signyy(i) )
292C
293 soundsp(i) = sqrt(a11t/rho0(i))
294 viscmax(i) = zero
295 etse(i) = one
296C-------------------
297C STRAIN RATE
298C-------------------
299 IF (israte == 0) THEN
300 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
301 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
302 . + epspxy(i)*epspxy(i) ) )
303 ELSE
304 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
305 ENDIF
306 jjc(i) = 1
307 jjt(i) = ncc + 1
308 ENDDO
309C
310C id function
311C
312 IF(irate <=3) THEN
313 DO i=1,nel
314 ipos1(i) = nint(uvar(i,5))
315 iad1(i) = npf(ifunc(1)) / 2 + 1
316 ilen1(i) = npf(ifunc(1)+1) / 2 - iad1(i)-ipos1(i)
317 ipos2(i) = nint(uvar(i,6))
318 iad2(i) = npf(ifunc(2)) / 2 + 1
319 ilen2(i) = npf(ifunc(2) + 1) / 2 - iad2(i)-ipos2(i)
320C
321 uvar(i,5) = ipos1(i)
322 uvar(i,6) = ipos2(i)
323 END DO
324C
325 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,yc)
326 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,yt)
327
328 IF(fisokin == zero) THEN
329 DO i=1,nel
330 yc(i)=yc(i)*yfac(1)
331 yt(i)=yt(i)*yfac(2)
332 hc(i)=dydx1(i)*yfac(1)
333 ht(i)=dydx2(i)*yfac(2)
334 ENDDO
335 ELSEIF(fisokin == one ) THEN
336 DO i=1,nel
337 hc(i)=dydx1(i)*yfac(1)
338 ht(i)=dydx2(i)*yfac(2)
339 yc(i)=tf(npf(ifunc(1)) + 1)
340 yt(i)=tf(npf(ifunc(2)) + 1)
341 yc(i)=yc(i)*yfac(1)
342 yt(i)=yt(i)*yfac(2)
343 yc(i) = max(em20, yc(i))
344 yt(i) = max(em20, yt(i))
345 ENDDO
346 ELSE
347 DO i=1,nel
348 yc(i)=yc(i)*yfac(1)
349 yt(i)=yt(i)*yfac(2)
350 yc(i) = max(yc(i),em20)
351 yt(i) = max(yt(i),em20)
352 hc(i) = dydx1(i)*yfac(1)
353 ht(i) = dydx2(i)*yfac(2)
354C ECROUISSAGE CINEMATIQUE
355 y1(i)=yfac(1)*tf(npf(ifunc(1))+1)
356 y2(i)=yfac(2)*tf(npf(ifunc(2))+1)
357 yc(i) = (one - fisokin) * yc(i) + fisokin * y1(i)
358 yt(i) = (one - fisokin) * yt(i) + fisokin * y2(i)
359 ENDDO
360 ENDIF
361 ELSE
362C multiples curves
363C
364C compression
365C
366 DO j = 2,ncc - 1
367 DO i=1,nel
368 IF(epsp(i) >= rate0(j) ) jjc(i) = j
369 ENDDO
370 ENDDO
371 DO i=1,nel
372 fac=rate0(jjc(i))
373 rate(i)=(epsp(i) - fac)/(rate0(jjc(i) +1) - fac)
374 ENDDO
375 DO i=1,nel
376 j1 = jjc(i)
377 j2 = j1+1
378 ipos1(i) = nint(uvar(i,4+j1 ))
379 iad1(i) = npf(ifunc(j1)) / 2 + 1
380 ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i)-ipos1(i)
381 ipos2(i) = nint(uvar(i,4+j2))
382 iad2(i) = npf(ifunc(j2)) / 2 + 1
383 ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i)-ipos2(i)
384C
385 uvar(i,4+j1) = ipos1(i)
386 uvar(i,4+j2) = ipos2(i)
387 END DO
388 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
389 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
390C
391 IF (fisokin == zero) THEN
392 DO i=1,nel
393 j1 = jjc(i)
394 j2 = j1+1
395 y1(i)=y1(i)*yfac(j1)
396 y2(i)=y2(i)*yfac(j2)
397 fac = rate(i)
398 yc(i) =(y1(i) + fac*(y2(i)-y1(i)))
399 yc(i) = max(yc(i),em20)
400 dydx1(i)=dydx1(i)*yfac(j1)
401 dydx2(i)=dydx2(i)*yfac(j2)
402 hc(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
403 ENDDO
404 ELSEIF (fisokin == one ) THEN
405 DO i=1,nel
406 j1 = jjc(i)
407 j2 = j1+1
408 fac = rate(i)
409 dydx1(i)=dydx1(i)*yfac(j1)
410 dydx2(i)=dydx2(i)*yfac(j2)
411 hc(i) =(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
412C ECROUISSAGE CINEMATIQUE
413 y1(i)=tf(npf(ifunc(j1)) + 1)
414 y2(i)=tf(npf(ifunc(j2)) + 1)
415 y1(i)=y1(i)*yfac(j1)
416 y2(i)=y2(i)*yfac(j2)
417 yc(i) = (y1(i) + fac*(y2(i)-y1(i)))
418 yc(i) = max(em20,yc(i))
419 ENDDO
420 ELSE
421 DO i=1,nel
422 j1 = jjc(i)
423 j2 = j1 + 1
424 y1(i)=y1(i)*yfac(j1)
425 y2(i)=y2(i)*yfac(j2)
426 fac = rate(i)
427 yc(i) = (y1(i) + fac*(y2(i)-y1(i)))
428 yc(i) = max(yc(i),em20)
429 dydx1(i)=dydx1(i)*yfac(j1)
430 dydx2(i)=dydx2(i)*yfac(j2)
431 hc(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
432C ECROUISSAGE CINEMATIQUE
433 y1(i)=tf(npf(ifunc(j1))+1)
434 y2(i)=tf(npf(ifunc(j2))+1)
435 y1(i)=y1(i)*yfac(j1)
436 y2(i)=y2(i)*yfac(j2)
437 yc(i) = (one - fisokin) * yc(i) +
438 . fisokin * ((y1(i) + fac*(y2(i)-y1(i))))
439 ENDDO
440 ENDIF
441C
442C traction
443C
444 DO j = 2,nct - 1
445 DO i=1,nel
446 IF(epsp(i) >= rate0(ncc+j)) jjt(i) = ncc + j
447 ENDDO
448 ENDDO
449 DO i=1,nel
450 fac=rate0(jjt(i))
451 rate(i)=(epsp(i) - fac)/(rate0(jjt(i)+1) - fac)
452 ENDDO
453 DO i=1,nel
454 j1 = jjt(i)
455 j2 = j1+1
456 ipos1(i) = nint(uvar(i,4+j1))
457 iad1(i) = npf(ifunc(j1)) / 2 + 1
458 ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i)-ipos1(i)
459 ipos2(i) = nint(uvar(i,4+j2))
460 iad2(i) = npf(ifunc(j2)) / 2 + 1
461 ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i)-ipos2(i)
462C
463 uvar(i,4+j1) = ipos1(i)
464 uvar(i,4+j2) = ipos2(i)
465 END DO
466 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
467 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
468C
469 IF (fisokin == zero) THEN
470 DO i=1,nel
471 j1 = jjt(i)
472 j2 = j1+1
473 y1(i)=y1(i)*yfac(j1)
474 y2(i)=y2(i)*yfac(j2)
475 fac = rate(i)
476 yt(i) =(y1(i) + fac*(y2(i)-y1(i)))
477 yt(i) = max(yt(i),em20)
478 dydx1(i)=dydx1(i)*yfac(j1)
479 dydx2(i)=dydx2(i)*yfac(j2)
480 ht(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
481 ENDDO
482 ELSEIF (fisokin == one ) THEN
483 DO i=1,nel
484 j1 = jjt(i)
485 j2 = j1+1
486 fac = rate(i)
487 dydx1(i)=dydx1(i)*yfac(j1)
488 dydx2(i)=dydx2(i)*yfac(j2)
489 ht(i) =(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
490C ECROUISSAGE CINEMATIQUE
491 y1(i)=yfac(j1)*tf(npf(ifunc(j1)) + 1)
492 y2(i)=yfac(j2)*tf(npf(ifunc(j2)) + 1)
493 yt(i) = (y1(i) + fac*(y2(i)-y1(i)))
494 yt(i) = max(em20,yt(i))
495 ENDDO
496 ELSE
497 DO i=1,nel
498 j1 = jjt(i)
499 j2 = j1 + 1
500 y1(i)=y1(i)*yfac(j1)
501 y2(i)=y2(i)*yfac(j2)
502 fac = rate(i)
503 yt(i) = (y1(i) + fac*(y2(i)-y1(i)))
504 yt(i) = max(yt(i),em20)
505 dydx1(i)=dydx1(i)*yfac(j1)
506 dydx2(i)=dydx2(i)*yfac(j2)
507 ht(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
508C ECROUISSAGE CINEMATIQUE
509 y1(i)=yfac(j1)*tf(npf(ifunc(j1))+1)
510 y2(i)=yfac(j2)*tf(npf(ifunc(j2))+1)
511 yt(i) = (one - fisokin) * yt(i) +
512 . fisokin * ((y1(i) + fac*(y2(i)-y1(i))))
513 ENDDO
514 ENDIF
515 ENDIF
516c
517C strain effect on compression and traction
518 IF(irate == 3) THEN
519 DO i=1,nel
520C compression
521 yrate = yfac(3)*finter(ifunc(3),epsp(i),npf,tf,df)
522 yc(i) = yrate*yc(i)
523 hc(i) = hc(i)*yrate
524c traction
525 yrate = yfac(4)*finter(ifunc(4),epsp(i),npf,tf,df)
526 yt(i) = yrate*yt(i)
527 ht(i) = ht(i)*yrate
528 ENDDO
529 ENDIF
530C interpolation entre (tracation & compression
531 DO i=1,nel
532 IF(pc == zero .AND. pt == zero .AND. abs(p(i)) < em10) THEN
533 yld(i) = max(yc(i), em20)
534 h(i) = hc(i)
535 ELSEIF(p(i) <= -pt) THEN
536 yld(i) = max(yt(i),em20)
537 h(i) = ht(i)
538 ELSEIF(p(i) >= pc) THEN
539 yld(i) = max(yc(i), em20)
540 h(i) = hc(i)
541 ELSE
542 fac = pc + pt
543 fac = (pc - p(i))/fac
544 yld(i) = fac*yt(i) + (one -fac)*yc(i)
545 yld(i) = max(em20,yld(i))
546 h(i) = fac*ht(i) + (one -fac)*hc(i)
547 ENDIF
548 ENDDO
549C
550C strain rate effect
551c (1 + (epsp/epsp0)**(1/cp)
552 IF(vp == 0) THEN
553 IF(irate == 1) THEN
554 DO i=1,nel
555 epd = max(em20,epsp(i)/epsp0)
556 yrate = one + exp(cp*log(epd))
557 yld(i) = yld(i)*yrate
558 h(i) = h(i)*yrate
559 ENDDO
560C 1 + cp*ln(epep/epsp0)
561 ELSEIF(irate == 2) THEN
562 DO i=1,nel
563 epd = max(em20,epsp(i)/epsp0)
564 yrate = one + cp*log(epd)
565 yld(i) = yld(i)*yrate
566 h(i) = h(i)*yrate
567 ENDDO
568 ENDIF
569 ENDIF
570C-------------------
571C PROJECTION
572C-------------------
573 IF(iflag(1) == 0)THEN
574 nu31 = one-nnu11
575C projection radiale
576 DO i=1,nel
577 svm(i)=sqrt(signxx(i)*signxx(i)
578 . +signyy(i)*signyy(i)
579 . -signxx(i)*signyy(i)
580 . +three*signxy(i)*signxy(i))
581 r = min(one,yld(i)/max(em20,svm(i)))
582 signxx(i)=signxx(i)*r
583 signyy(i)=signyy(i)*r
584 signxy(i)=signxy(i)*r
585 umr = one - r
586 dpla_i(i) = off(i)*svm(i)*umr/(g3(i)+h(i))
587 pla(i) = pla(i) + dpla_i(i)
588 s1=half*(signxx(i)+signyy(i))
589 IF (inloc ==0) THEN
590 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
591 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
592 thk(i) = thk(i) + dezz*thkly(i)*off(i)
593 ENDIF
594 IF(r<one) etse(i)= h(i)/(h(i)+e(i))
595 ENDDO
596C
597 ELSEIF(iflag(1)==1)THEN
598C-------------------------
599C CRITERE DE VON MISES
600C-------------------------
601 DO i=1,nel
602 h(i) = max(zero,h(i))
603 s1=signxx(i)+signyy(i)
604 s2=signxx(i)-signyy(i)
605 s3=signxy(i)
606 aa(i)=fourth*s1*s1
607 bb(i)=three_over_4*s2*s2+3.*s3*s3
608 svm(i)=sqrt(aa(i)+bb(i))
609 IF (inloc == 0) THEN
610 dezz = -(depsxx(i)+depsyy(i))*nnu11
611 thk(i) = thk(i) + dezz*thkly(i)*off(i)
612 ENDIF
613 ENDDO
614C-------------------------
615C GATHER PLASTIC FLOW
616C-------------------------
617 nindx=0
618 DO i=1,nel
619 IF(svm(i)>yld(i).AND.off(i)==one) THEN
620 nindx=nindx+1
621 index(nindx)=i
622 ENDIF
623 ENDDO
624C
625 IF(nindx/=0) THEN
626C---------------------------
627C DEP EN CONTRAINTE PLANE
628C---------------------------
629 DO j=1,nindx
630 i=index(j)
631 dpla_j(i)=(svm(i)-yld(i))/(g3(i)+h(i))
632 etse(i)= h(i)/(h(i)+e(i))
633 hi(i) = h(i)*(one-fisokin)
634 hk(i) = two_third*h(i)*fisokin
635 ENDDO
636C
637 IF(vp == 1) THEN
638 DO j=1,nindx
639 i=index(j)
640 dtinv = timestep(i)/max(em20, timestep(i)**2)
641 epd = dpla_j(i)*dtinv
642 epd = max(em20,epd/epsp0)
643 yrate = one + exp(cp*log(epd))
644 IF(sigy == zero) THEN
645 yld(i) = yld(i)*yrate
646 hi(i) = hi(i)*yrate
647 hk(i) = hk(i)*yrate
648 ELSE
649 yld(i) = yld(i) + sigy*(yrate - one)
650 ENDIF
651 ENDDO
652 ENDIF
653C
654 nu11 = one/(one-nu)
655 nu21 = one/(one+nu)
656 nu31 = one-nnu11
657 DO n=1,nmax
658#include "vectorize.inc"
659 DO j=1,nindx
660 i=index(j)
661 dpla_i(i)=dpla_j(i)
662 yld_i =yld(i)+hi(i)*dpla_i(i)
663 dr(i) =half*e(i)*dpla_i(i)/yld_i
664 nu110 = nu11+three*hk(i)/e(i)
665 nu210 = nu21+hk(i)/e(i)
666 pp(i) =one/(one+dr(i)*nu110)
667 qq(i) =one/(one+three*dr(i)*nu210)
668 p2 =pp(i)*pp(i)
669 q2 =qq(i)*qq(i)
670 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
671 df =-(aa(i)*nu110*p2*pp(i)+three*bb(i)*nu210*q2*qq(i))
672 . *(e(i)-two*dr(i)*hi(i))/yld_i
673 . -two*hi(i)*yld_i
674 df = sign(max(abs(df),em20),df)
675 IF(dpla_i(i)>zero) THEN
676 dpla_j(i)=max(zero,dpla_i(i)-f/df)
677 ELSE
678 dpla_j(i)=zero
679 ENDIF
680 ENDDO
681 ENDDO
682C------------------------------------------
683C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
684C------------------------------------------
685#include "vectorize.inc"
686 DO j=1,nindx
687 i=index(j)
688 pla(i) = pla(i) + dpla_i(i)
689 s1=(signxx(i)+signyy(i))*pp(i)
690 s2=(signxx(i)-signyy(i))*qq(i)
691 signxx(i)=half*(s1+s2)
692 signyy(i)=half*(s1-s2)
693 signxy(i)=signxy(i)*qq(i)
694 IF (inloc == 0) THEN
695 dezz = - nu31*dr(i)*s1/e(i)
696 thk(i) = thk(i) + dezz*thkly(i)*off(i)
697 ENDIF
698 ENDDO
699 ENDIF
700C-------------------------------------------
701 ELSEIF(iflag(1)==2)THEN
702C-------------------------
703C CRITERE DE VON MISES
704C-------------------------
705 DO i=1,nel
706 h(i) = max(zero,h(i))
707 svm2(i)=signxx(i)*signxx(i)
708 . +signyy(i)*signyy(i)
709 . -signxx(i)*signyy(i)
710 . +three*signxy(i)*signxy(i)
711 svm(i)=sqrt(svm2(i))
712 IF (inloc == 0) THEN
713 dezz = -(depsxx(i)+depsyy(i))*nnu11
714 thk(i) = thk(i) + dezz*thkly(i)*off(i)
715 ENDIF
716 ENDDO
717C-------------------------
718C GATHER PLASTIC FLOW
719C-------------------------
720 nindx=0
721 DO i=1,nel
722 yld2(i)=yld(i)*yld(i)
723 IF(svm2(i)>yld2(i).AND.off(i)==one) THEN
724 nindx=nindx+1
725 index(nindx)=i
726 ENDIF
727 ENDDO
728C
729 IF(nindx/=0) THEN
730C-------------
731C PROJ NORMALE AU CRITERE AVEC CALCUL APPROCHE DE LA NORMALE + RETOUR RADIAL
732C-------------
733 nu31 = one-nnu11
734 DO j=1,nindx
735 i=index(j)
736 a=(svm2(i)-yld2(i))
737 . /(five*svm2(i)+three*(-signxx(i)*signyy(i)+signxy(i)*signxy(i)))
738 s1=(one-two*a)*signxx(i)+ a*signyy(i)
739 s2= a*signxx(i)+(one-two*a)*signyy(i)
740 s3=(one-three*a)*signxy(i)
741 signxx(i)=s1
742 signyy(i)=s2
743 signxy(i)=s3
744 dpla_i(i) = off(i)*(svm(i)-yld(i))/(g3(i)+h(i))
745C
746 hk(i) = h(i)*(one-fisokin)
747 yld(i)= yld(i)+hk(i)*dpla_i(i)
748 END DO
749C
750 DO j=1,nindx
751 i=index(j)
752 svm(i)= signxx(i)*signxx(i) + signyy(i)*signyy(i)
753 . - signxx(i)*signyy(i) + three*signxy(i)*signxy(i)
754 IF (svm(i) > yld(i)*yld(i)) THEN
755 svm(i) = sqrt(svm(i))
756 r = yld(i) / svm(i)
757 signxx(i) = signxx(i)*r
758 signyy(i) = signyy(i)*r
759 signxy(i) = signxy(i)*r
760 pla(i) = pla(i) + dpla_i(i)
761 IF (inloc == 0) THEN
762 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
763 dezz = -nu31*dezz
764 thk(i) = thk(i) + dezz*thkly(i)*off(i)
765 ENDIF
766 etse(i)= h(i)/(h(i)+e(i))
767 ENDIF
768 END DO
769 END IF
770C
771 ENDIF
772C------------------------------------------
773C ECROUISSAGE CINE
774C------------------------------------------
775
776 IF (fisokin > zero) THEN
777 DO i=1,nel
778 dsxx = sigexx(i) - signxx(i)
779 dsyy = sigeyy(i) - signyy(i)
780 dsxy = sigexy(i) - signxy(i)
781 dexx = (dsxx - nu*dsyy)
782 deyy = (dsyy - nu*dsxx)
783C
784 dexy = two*(one+nu)*dsxy
785 hkin = two_third*fisokin*h(i)
786 alpha = hkin/(e(i)+hkin)
787 sigpxx = alpha*(two*dexx+deyy)
788 sigpyy = alpha*(two*deyy+dexx)
789 sigpxy = alpha*dexy*half
790 uvar(i,2) = uvar(i,2) + sigpxx
791 uvar(i,3) = uvar(i,3) + sigpyy
792 uvar(i,4) = uvar(i,4) + sigpxy
793C
794 signxx(i) = signxx(i) + uvar(i,2)
795 signyy(i) = signyy(i) + uvar(i,3)
796 signxy(i) = signxy(i) + uvar(i,4)
797 ENDDO
798 ENDIF
799C
800C visco elastic model (prony) moved to mulawc
801C
802C--------------------------------
803C NON-LOCAL THICKNESS VARIATION
804C--------------------------------
805 IF (inloc > 0) THEN
806 DO i = 1,nel
807 IF (loff(i) == one) THEN
808 svm(i) = sqrt(signxx(i)*signxx(i) + signyy(i)*signyy(i)
809 . - signxx(i)*signyy(i) + three*signxy(i)*signxy(i))
810 dezz = max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(svm(i),em20)
811 dezz = -nu*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e(i)) - dezz
812 thk(i) = thk(i) + dezz*thkly(i)*off(i)
813 ENDIF
814 ENDDO
815 ENDIF
816C -----------------------------------------
817 RETURN
818 END
819C
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
sigeps66c
subroutine sigeps66c(nel, nuparam, nuvar, mfunc, kfunc, npf, npt0, ipt, iflag, tf, time, timestep, uparam, rho0, area, eint, thkly, israte, asrate, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thk, pla, uvar, off, ngl, ipm, mat, etse, gs, yld, epsd_pg, epsp, inloc, dplanl, mat_param, nuvarv, uvarv, loff)
Definition sigeps66c.F:49
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72