OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps96.F File Reference
#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps96 (nel, ngl, nuparam, nuvar, nfunc, ifunc, npf, tf, uparam, uvar, jthe, rho, temp, defp, soundsp, off, epsd, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx)

Function/Subroutine Documentation

◆ sigeps96()

subroutine sigeps96 ( integer nel,
integer, dimension(nel) ngl,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
uparam,
target uvar,
integer jthe,
intent(in) rho,
intent(in) temp,
target defp,
intent(out) soundsp,
intent(in) off,
intent(inout) epsd,
intent(in) epspxx,
intent(in) epspyy,
intent(in) epspzz,
intent(in) epspxy,
intent(in) epspyz,
intent(in) epspzx,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signzz,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx )

Definition at line 28 of file sigeps96.F.

36C-----------------------------------------------
37C I m p l i c i t T y p e s
38C-----------------------------------------------
39#include "implicit_f.inc"
40C-----------------------------------------------
41C C O M M O N
42C-----------------------------------------------
43C---------+---------+---+---+--------------------------------------------
44C VAR | SIZE |TYP| RW| DEFINITION
45C---------+---------+---+---+--------------------------------------------
46C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
47C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
48C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
49C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
50C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
51C RHO | NEL | F | R | DENSITY
52C TEMP | NEL | F | R | ELEMENT TEMPERATURE (TSTAR)
53C DEFP | NEL | F | R | PLASTIC STRAIN
54C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
55C---------+---------+---+---+--------------------------------------------
56C EPSPXX | NEL | F | R | STRAIN RATE XX
57C EPSPYY | NEL | F | R | STRAIN RATE YY
58C ... | | | |
59C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
60C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
61C ... | | | |
62C---------+---------+---+---+--------------------------------------------
63C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
64C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
65C ... | | | |
66C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
67C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
68C ... | | | |
69C---------+---------+---+---+--------------------------------------------
70C I N P U T A r g u m e n t s
71C-----------------------------------------------
72 INTEGER NEL,NUPARAM,NUVAR,NFUNC,JTHE
73 INTEGER NPF(*),NGL(NEL),IFUNC(NFUNC)
74 my_real tf(*),uparam(nuparam)
75 my_real,DIMENSION(NEL), INTENT(IN) :: rho,temp,off,
76 . epspxx,epspyy,epspzz,epspxy,epspyz,epspzx,
77 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
78 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
79C-----------------------------------------------
80C O U T P U T A r g u m e n t s
81C-----------------------------------------------
82 my_real ,DIMENSION(NEL), INTENT(OUT) :: soundsp,
83 . signxx,signyy,signzz,signxy,signyz,signzx
84C-----------------------------------------------
85C I N P U T O U T P U T A r g u m e n t s
86C-----------------------------------------------
87 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT), TARGET :: uvar
88 my_real ,DIMENSION(NEL,2), INTENT(INOUT), TARGET :: defp
89 my_real ,DIMENSION(NEL), INTENT(INOUT) :: epsd
90C-----------------------------------------------
91C L o c a l V a r i a b l e s
92C-----------------------------------------------
93 INTEGER I,ITER,NITER,IVARF,IFUNC_E,IFUNC_YLD
94 INTEGER, DIMENSION(NEL) :: IPOS1,ILEN1,IAD1,IPOS2,ILEN2,IAD2
95c
96 my_real tanb,tanp,eini,kini,gini,nu,yxi,kep,h,ho,hm,res,dres,
97 . sigy,siga,sigb,sigr,ra1,ra2,rb,rc,na,fac,q2,pp1,pp2,pla2,pla,dpla,
98 . dfdp,dfdc,dgdp,denom,ldot,dwpx2,qx2,depslv,depspv,depspd,
99 . e1,e2,e3,e4,e5,e6,tr,sy1,sy2,sy3,sy4,sy5,sy6,
100 . dfds1,dfds2,dfds3,dfds4,dfds5,dfds6,depsp1,depsp2,depsp3,depsp4,
101 . depsp5,depsp6,depsl1,depsl2,depsl3,depsl4,depsl5,depsl6,epsp,
102 . qy,fy,py,f1,f2,x,x1,x2,dx,alpha
103 my_real, DIMENSION(NEL) :: k,g,g2,ra,depsv,dav,d1,d2,d3,
104 . ds1,ds2,ds3,ds4,ds5,ds6,dp,de1,de2,de3,de4,de5,de6,dydx,
105 . fo,po,so1,so2,so3,so4,so5,so6,s1,s2,s3,s4,s5,s6,p,yld,
106 . s1n,s2n,s3n,s4n,s5n,s6n,pn,qn,q,f,fac_e,frate
107 my_real small,big,tol
108 my_real ,DIMENSION(:), POINTER :: epspd,epspv
109C=======================================================================
110c UPARAM(1) E : Young modulus
111c UPARAM(2) G : Shear modulus
112c UPARAM(3) K : Bulk modulus
113c UPARAM(4) NU : Poisson coefficient
114c UPARAM(5) TANB : Tangent beta (yield friction angle)
115c UPARAM(6) TANP : Tangent Psi (plastic flow angle)
116c UPARAM(7) YXI : sqrt(1 + TANB^2)
117c UPARAM(8) SIGY : initial yield value in pure shear (p=0)
118c UPARAM(9) SIGA : Max yield gain during early hardening peak
119c UPARAM(10) SIGB : Max Yield drop during softening phase
120c UPARAM(11) SIGR : Rubbery modulus (hardening modulus in rubbery state)
121c UPARAM(12) RA1 : 1st exponent coefficient for SIGA evolution
122c UPARAM(13) RA2 : 2nd exponent coefficient for SIGA evolution
123c UPARAM(14) NA : exponent coefficient for temperature dependency of SIGA
124c UPARAM(15) RB : exponent coefficient for SIGB evolution
125c UPARAM(16) RC : coefficient for SIGR evolution
126c UPARAM(20) NUP : Plastic deformation Poisson coefficient
127c---------------------
128c UVAR(1) YIELD
129c UVAR(2) EPSD
130C=======================================================================
131 niter = 3
132 small = em10
133 big = ep20
134 tol = em20
135c---------------------
136 epspd => defp(1:nel,1) ! Von Mises Equivalent Plastic Strain
137 epspv => defp(1:nel,2) ! Volumetric Plastic Strain
138c---------------------
139 eini = uparam(1)
140 gini = uparam(2)
141 kini = uparam(3)
142 nu = uparam(4)
143 tanb = uparam(5)
144 tanp = uparam(6)
145 yxi = uparam(7)
146 sigy = uparam(8)
147 siga = uparam(9)
148 sigb = uparam(10)
149 sigr = uparam(11)
150 ra1 = uparam(12)
151 ra2 = uparam(13)
152 na = uparam(14)
153 rb = uparam(15)
154 rc = uparam(16)
155 kep = uparam(17) ! 1 / (1 + nup^2) : d_epsp = kep*ldot
156 alpha= uparam(19) ! strain rate smoothing coefficient
157c
158 ifunc_e = ifunc(1)
159 ifunc_yld = ifunc(2)
160c-------------------------------------------
161c TEMPERATURE AND STRAIN RATE DEPENDENCY (YOUNG MODULUS, YIELD)
162c-------------------------------------------
163 ivarf = 2
164 IF (jthe == 1. and. ifunc_e > 0) THEN
165 DO i=1,nel
166 ipos1(i) = nint(uvar(i, ivarf + 1))
167 iad1(i) = npf(ifunc_e) / 2 + 1
168 ilen1(i) = npf(ifunc_e + 1) / 2 - iad1(i) - ipos1(i)
169 ENDDO
170c
171 CALL vinter2(tf,iad1,ipos1,ilen1,nel,temp,dydx,fac_e)
172 uvar(1:nel,ivarf + 1) = ipos1(1:nel)
173c
174 k(1:nel) = kini*fac_e(1:nel)
175 g(1:nel) = gini*fac_e(1:nel)
176 g2(1:nel) = g(1:nel)*two
177 ivarf = ivarf + 1
178 ELSE
179 k(1:nel) = kini
180 g(1:nel) = gini
181 g2(1:nel) = gini*two
182 ENDIF
183c---
184 IF (ifunc_yld > 0) THEN
185 DO i=1,nel
186 tr = (epspxx(i) + epspyy(i) + epspzz(i))*third
187 e1 = epspxx(i) - tr
188 e2 = epspyy(i) - tr
189 e3 = epspzz(i) - tr
190 e4 = half*epspxy(i)
191 e5 = half*epspyz(i)
192 e6 = half*epspzx(i)
193 epsp = half*(e1**2 + e2**2 + e3**2) + e4**2 + e5**2 + e6**2
194 epsp = sqrt(three*epsp)/three_half
195 epsp = alpha*epsp + (one -alpha)*epsd(i)
196 epsd(i) = epsp
197 ENDDO
198c
199 DO i=1,nel
200 ipos2(i) = nint(uvar(i, ivarf + 1))
201 iad2(i) = npf(ifunc_yld) / 2 + 1
202 ilen2(i) = npf(ifunc_yld + 1) / 2 - iad2(i) - ipos2(i)
203 ENDDO
204c
205 CALL vinter2(tf,iad2,ipos2,ilen2,nel,epsd,dydx,frate)
206 uvar(1:nel,ivarf + 1) = ipos2(1:nel)
207 ELSE
208 frate(1:nel) = one
209 ENDIF
210 uvar(1:nel,2) = epsd(1:nel)
211c---
212 IF (jthe == 1) THEN
213 ra(1:nel) = ra1*temp(1:nel)**na + ra2
214 ELSE
215 ra(1:nel) = ra1 + ra2
216 ENDIF
217c------------------------------------------
218c ELASTIC TRIAL INCREMENT
219c------------------------------------------
220 DO i=1,nel
221 pla = epspd(i)
222 pla2 = pla**2
223 pp1 = one - rc*pla2
224 pp2 = pp1 + two
225 yld(i) = sigy*frate(i) + siga*(one - exp(-ra(i)*pla))
226 . - sigb *(one - exp(-rb*pla))
227 . + sigr*pla*third*pp2 / pp1
228c
229 depsv(i) = (depsxx(i) + depsyy(i) + depszz(i))
230 dav(i) = depsv(i)*third
231 d1(i) = depsxx(i)-dav(i)
232 d2(i) = depsyy(i)-dav(i)
233 d3(i) = depszz(i)-dav(i)
234c Increments elastiques pression et deviateur
235 dp(i) =-k(i)* depsv(i)
236 ds1(i)= g2(i)* d1(i)
237 ds2(i)= g2(i)* d2(i)
238 ds3(i)= g2(i)* d3(i)
239 ds4(i)= g(i) * depsxy(i)
240 ds5(i)= g(i) * depsyz(i)
241 ds6(i)= g(i) * depszx(i)
242c Valeurs elastiques pression et deviateur
243 po(i) =-(sigoxx(i)+sigoyy(i)+sigozz(i))*third
244 so1(i) = sigoxx(i) + po(i)
245 so2(i) = sigoyy(i) + po(i)
246 so3(i) = sigozz(i) + po(i)
247 so4(i) = sigoxy(i)
248 so5(i) = sigoyz(i)
249 so6(i) = sigozx(i)
250
251 s1(i) = so1(i) + ds1(i)
252 s2(i) = so2(i) + ds2(i)
253 s3(i) = so3(i) + ds3(i)
254 s4(i) = so4(i) + ds4(i)
255 s5(i) = so5(i) + ds5(i)
256 s6(i) = so6(i) + ds6(i)
257 p(i) = po(i) + dp(i)
258c Contraintes elastiques
259 signxx(i) = s1(i) - p(i)
260 signyy(i) = s2(i) - p(i)
261 signzz(i) = s3(i) - p(i)
262 signxy(i) = s4(i)
263 signyz(i) = s5(i)
264 signzx(i) = s6(i)
265 ENDDO
266c------------------------------------------
267c EVALUATION FONCTION CRITERE DE PLASTICITE
268c------------------------------------------
269c old criterion
270 DO i=1,nel
271 q2 = three_half*(so1(i)**2+so2(i)**2+so3(i)**2)
272 . + three*(so4(i)**2+so5(i)**2+so6(i)**2)
273 q(i) = sqrt(q2)
274 fo(i)= q(i) - max(zero, po(i)*tanb + yld(i))
275 ENDDO
276c------------------------------------------
277 DO i=1,nel
278 q2 = three_half*(s1(i)**2+s2(i)**2+s3(i)**2)
279 . + three*(s4(i)**2+s5(i)**2+s6(i)**2)
280 q(i)= sqrt(q2)
281 f(i)= q(i) - max(zero, p(i)*tanb + yld(i))
282 ENDDO
283c------------------------------------------
284c ITERATION POUR TROUVER LE PT D'INTERSECTION AVEC LE CRITERE
285c------------------------------------------
286 DO i=1,nel
287 IF (f(i) > tol) THEN
288c
289 IF (fo(i) < zero) THEN
290 x1 = one
291 x2 = zero
292 f1 = fo(i)
293 f2 = f(i)
294 DO iter=1,niter
295 x = (x1*f2 - x2*f1) / (f2 - f1) ! regula falsi
296 sy1 = s1(i) - x * ds1(i)
297 sy2 = s2(i) - x * ds2(i)
298 sy3 = s3(i) - x * ds3(i)
299 sy4 = s4(i) - x * ds4(i)
300 sy5 = s5(i) - x * ds5(i)
301 sy6 = s6(i) - x * ds6(i)
302 py = p(i) - x * dp(i)
303 qy = three_half*(sy1**2+sy2**2+sy3**2) + three*(sy4**2+sy5**2+sy6**2)
304 qy = sqrt(qy)
305 fy = qy - max(zero, py*tanb + yld(i))
306c
307 IF (fy*f1 > 0) THEN
308 x1 = x
309 f1 = fy
310 ELSEIF (fy*f2 > 0) THEN
311 x2 = x
312 f2 = fy
313 ELSE
314 EXIT
315 ENDIF
316 ENDDO
317c
318c constraints
319 ds1(i) = x*ds1(i)
320 ds2(i) = x*ds2(i)
321 ds3(i) = x*ds3(i)
322 ds4(i) = x*ds4(i)
323 ds5(i) = x*ds5(i)
324 ds6(i) = x*ds6(i)
325 dp(i) = x*dp(i)
326c
327 s1(i) = sy1
328 s2(i) = sy2
329 s3(i) = sy3
330 s4(i) = sy4
331 s5(i) = sy5
332 s6(i) = sy6
333 p(i) = py
334 q(i) = qy
335c deformations
336 de1(i) = x*d1(i)
337 de2(i) = x*d2(i)
338 de3(i) = x*d3(i)
339 de4(i) = x*depsxy(i)
340 de5(i) = x*depsyz(i)
341 de6(i) = x*depszx(i)
342 depsv(i)=x*depsv(i)
343 ELSE
344 de1(i) = d1(i)
345 de2(i) = d2(i)
346 de3(i) = d3(i)
347 de4(i) = depsxy(i)
348 de5(i) = depsyz(i)
349 de6(i) = depszx(i)
350 ENDIF
351 ENDIF
352 ENDDO
353c---------------------------------
354c PLASTIC FLOW : FULL NEWTON ITERATIONS
355c---------------------------------
356 DO i=1,nel
357 IF (f(i) > zero) THEN
358 dfdp =-tanb
359 dgdp =-tanp
360 dfdc =-one
361c
362 IF (q(i) > em20) THEN
363 fac = three_half/q(i)
364 dfds1 = s1(i)*fac
365 dfds2 = s2(i)*fac
366 dfds3 = s3(i)*fac
367 dfds4 = s4(i)*fac
368 dfds5 = s5(i)*fac
369 dfds6 = s6(i)*fac
370 ELSE
371 dfds1 =zero
372 dfds2 =zero
373 dfds3 =zero
374 dfds4 =zero
375 dfds5 =zero
376 dfds6 =zero
377 ENDIF
378c
379 ho = siga*ra(i)*exp(-ra(i)*epspd(i)) - sigb*rb*exp(-rb*epspd(i))
380 . + sigr*(one + two_third*rc*pla2 * pp2 / pp1**2)
381 denom = three*g(i) + k(i)*dfdp*dgdp - dfdc*ho
382 IF (denom /= zero) THEN
383 ldot = f(i) / denom
384 ELSE
385 ldot = zero
386 ENDIF
387c-----------------------------------------------------------------------
388 niter = 4
389 DO iter=1,niter
390c
391 depspd = ldot
392 depspv = ldot*dgdp*third
393c
394 depsp1 = ldot*dfds1
395 depsp2 = ldot*dfds2
396 depsp3 = ldot*dfds3
397 depsp4 = ldot*dfds4 * two
398 depsp5 = ldot*dfds5 * two
399 depsp6 = ldot*dfds6 * two
400c
401c Deviatoric elastic deformations
402c
403 depsl1 = de1(i) - depsp1
404 depsl2 = de2(i) - depsp2
405 depsl3 = de3(i) - depsp3
406 depsl4 = de4(i) - depsp4
407 depsl5 = de5(i) - depsp5
408 depsl6 = de6(i) - depsp6
409 depslv = depsv(i) + depspv
410c
411c Stress increment starting from yield surface
412c
413 ds1(i)= g2(i)*depsl1
414 ds2(i)= g2(i)*depsl2
415 ds3(i)= g2(i)*depsl3
416 ds4(i)= g(i) *depsl4
417 ds5(i)= g(i) *depsl5
418 ds6(i)= g(i) *depsl6
419 dp(i) =-k(i) *depslv
420c
421c New stress
422c
423 s1n(i) = s1(i) + ds1(i)
424 s2n(i) = s2(i) + ds2(i)
425 s3n(i) = s3(i) + ds3(i)
426 s4n(i) = s4(i) + ds4(i)
427 s5n(i) = s5(i) + ds5(i)
428 s6n(i) = s6(i) + ds6(i)
429 pn(i) = p(i) + dp(i)
430c
431c Update yield and plastic increment
432c
433 q2 = three_half*(s1n(i)**2 + s2n(i)**2 + s3n(i)**2)
434 . + three*(s4n(i)**2 + s5n(i)**2 + s6n(i)**2)
435 qn(i)= sqrt(q2)
436c
437 dpla = ldot*kep
438 pla = epspd(i) + dpla
439 pla2 = pla**2
440 pp1 = one - rc*pla2
441 pp2 = pp1 + two
442 yld(i) = sigy*frate(i) + siga*(one - exp(-ra(i)*pla))
443 . - sigb *(one - exp(-rb*pla))
444 . + sigr*pla*third*pp2 / pp1
445 yld(i) = max(yld(i), zero)
446
447c H = SIGA*RA(I)*EXP(-RA(I)*EPSPD(I)) - SIGB*RB*EXP(-RB*EPSPD(I))
448c . + SIGR*(ONE + TWO_THIRD*RC*PLA2 * PP2 / PP1**2)
449c HM = (HO + H) * HALF
450 hm = ho
451
452 res = qn(i) - max(zero, pn(i)*tanb + yld(i))
453 dres = -(three*g(i) + k(i)*dfdp*dgdp - dfdc*hm)
454
455 IF (dres /= zero) THEN
456 ldot = ldot - res / dres
457 ENDIF
458c
459 ENDDO ! ITER=1,NITER
460c-----------------------------------------------------------------------
461c CHECK REPROJECTION
462c
463 fy = qn(i) - max(zero, pn(i)*tanb + yld(i))
464c
465 IF (fy > small) THEN
466 IF (tanb*pn(i) + yld(i) <= zero) THEN
467 pn(i) =-yld(i)/tanb
468 s1n(i) = zero
469 s2n(i) = zero
470 s3n(i) = zero
471 s4n(i) = zero
472 s5n(i) = zero
473 s6n(i) = zero
474c WRITE(6,*)'TRI-TRACTION FAILURE'
475 ELSE
476 IF (qn(i) > small) THEN
477 x = (pn(i)*tanb + yld(i)) / qn(i)
478 IF (x < one-em02 .OR. x > one) THEN
479 WRITE(7,*)'REPROJ Q',x,fy,pn(i),qn(i)
480 WRITE(6,*)'REPROJ Q',x,fy,pn(i),qn(i)
481 ENDIF
482 s1n(i) = x*s1n(i)
483 s2n(i) = x*s2n(i)
484 s3n(i) = x*s3n(i)
485 s4n(i) = x*s4n(i)
486 s5n(i) = x*s5n(i)
487 s6n(i) = x*s6n(i)
488 ENDIF
489 ENDIF
490 ENDIF
491 epspd(i) = epspd(i) + ldot*kep
492 epspv(i) = epspv(i) + ldot*dgdp
493c
494c Final stress update
495c
496 signxx(i) = s1n(i) - pn(i)
497 signyy(i) = s2n(i) - pn(i)
498 signzz(i) = s3n(i) - pn(i)
499 signxy(i) = s4n(i)
500 signyz(i) = s5n(i)
501 signzx(i) = s6n(i)
502c
503 ENDIF ! F > 0
504c
505 uvar(i,1)= yld(i)
506c
507 ENDDO ! I=1,NEL
508c----------------------------------------------
509c END OF PLASTICITY PROJECTION LOOP
510c----------------------------------------------
511 DO i=1,nel
512 soundsp(i) = sqrt((k(i) + four_over_3*g(i))/rho(i))
513 ENDDO
514c-----------
515 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143