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

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps187 (nel, nuparam, nuvar, nfunc, ifunc, 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, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, yld, pla, dpla, etse, ipm, mat, israte, facyldi, epsp)

Function/Subroutine Documentation

◆ sigeps187()

subroutine sigeps187 ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) 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,
sigvxx,
sigvyy,
sigvzz,
sigvxy,
sigvyz,
sigvzx,
soundsp,
viscmax,
uvar,
off,
yld,
pla,
dpla,
etse,
integer, dimension(npropmi,*) ipm,
integer, dimension(nel) mat,
integer, dimension(nel) israte,
facyldi,
epsp )

Definition at line 31 of file sigeps187.F.

44C=======================================================================
45C BARLAT YLD2000 material model
46C=======================================================================
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C G L O B A L P A R A M E T E R S
53C-----------------------------------------------
54#include "com01_c.inc"
55#include "param_c.inc"
56C-----------------------------------------------
57C C O M M O N
58C-----------------------------------------------
59C---------+---------+---+---+--------------------------------------------
60C VAR | SIZE |TYP| RW| DEFINITION
61C---------+---------+---+---+--------------------------------------------
62C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
63C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
64C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
65C---------+---------+---+---+--------------------------------------------
66C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
67C IFUNC | NFUNC | I | R | FUNCTION INDEX
68C NPF | * | I | R | FUNCTION ARRAY
69C TF | * | F | R | FUNCTION ARRAY
70C---------+---------+---+---+--------------------------------------------
71C TIME | 1 | F | R | CURRENT TIME
72C TIMESTEP| 1 | F | R | CURRENT TIME STEP
73C UVAR | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
74C RHO0 | NEL | F | R | INITIAL DENSITY
75C RHO | NEL | F | R | DENSITY
76C VOLUME | NEL | F | R | VOLUME
77C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
78C EPSPXX | NEL | F | R | STRAIN RATE XX
79C EPSPYY | NEL | F | R | STRAIN RATE YY
80C ... | | | |
81C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
82C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
83C ... | | | |
84C EPSXX | NEL | F | R | STRAIN XX
85C EPSYY | NEL | F | R | STRAIN YY
86C ... | | | |
87C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
88C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
89C ... | | | |
90C---------+---------+---+---+--------------------------------------------
91C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
92C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
93C ... | | | |
94C SIGVXX | NEL | F | W | VISCOUS STRESS XX
95C SIGVYY | NEL | F | W | VISCOUS STRESS YY
96C ... | | | |
97C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
98C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
99C---------+---------+---+---+--------------------------------------------
100C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
101C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
102C---------+---------+---+---+--------------------------------------------
103C-----------------------------------------------
104C I N P U T A r g u m e n t s
105C-----------------------------------------------
106C
107 INTEGER :: NEL, NUPARAM, NUVAR
108 INTEGER :: IPM(NPROPMI,*),ISRATE(NEL),MAT(NEL)
109 my_real
110 . time,timestep,uparam(nuparam),
111 . rho(nel),rho0(nel),volume(nel),eint(nel),
112 . epspxx(nel),epspyy(nel),epspzz(nel),
113 . epspxy(nel),epspyz(nel),epspzx(nel),
114 . depsxx(nel),depsyy(nel),depszz(nel),
115 . depsxy(nel),depsyz(nel),depszx(nel),
116 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,
117 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
118 . sigoxx(nel),sigoyy(nel),sigozz(nel),
119 . sigoxy(nel),sigoyz(nel),sigozx(nel),
120 . facyldi(nel),epsp(nel)
121C-----------------------------------------------
122C O U T P U T A r g u m e n t s
123C-----------------------------------------------
124 my_real
125 . signxx(nel),signyy(nel),signzz(nel),
126 . signxy(nel),signyz(nel),signzx(nel),
127 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
128 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
129 . soundsp(nel),viscmax(nel),yld(nel) ,
130 . pla(nel),dpla(nel),etse(nel)
131C-----------------------------------------------
132C I N P U T O U T P U T A r g u m e n t s
133C-----------------------------------------------
134 my_real
135 . uvar(nel,nuvar), off(nel)
136C-----------------------------------------------
137C VARIABLES FOR FUNCTION INTERPOLATION
138C-----------------------------------------------
139 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
140 my_real
141 . finter ,tf(*)
142 EXTERNAL finter
143c Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
144c Y : y = f(x)
145c X : x
146c DYDX : f'(x) = dy/dx
147c IFUNC(J): FUNCTION INDEX
148c J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
149c NPF,TF : FUNCTION PARAMETER
150C-----------------------------------------------
151C L o c a l V a r i a b l e s
152C-----------------------------------------------
153 INTEGER I, II, J, K, NRATE, JJ(NEL), J1, J2, NINDX, IFLAG, IFLAGSR, EXPA, EXPAM2,
154 . NMAX,INDX(NEL), IPOS1(NEL),ILEN1(NEL),IAD1(NEL),
155 . IPOS2(NEL),ILEN2(NEL),IAD2(NEL)
156 my_real
157 . y1(nel),y2(nel),dydx(nel),dydx1(nel),dydx2(nel),
158 . xpxx(nel),xpyy(nel),xpxy(nel),xpyz(nel),xpzx(nel),
159 . xppxx(nel),xppyy(nel),xppxy(nel),xppyz(nel),xppzx(nel),
160 . phip(nel),phipp(nel),
161 . deltap(nel), deltapp1(nel), deltapp2(nel), deltapp(nel),
162 . deplxx(nel),deplyy(nel),deplxy(nel),deplzz(nel),
163 . sig(nel),h(nel),sigtrxx(nel),sigtryy(nel),
164 . sigtrxy(nel),rate(nel),sigtryz(nel),sigtrzz(nel),
165 . sigtrzx(nel),deplyz(nel),deplzx(nel),ddeplyz(nel),ddeplzx(nel),
166 . ddeplxx(nel),ddeplyy(nel),ddeplxy(nel),ddeplzz(nel)
167 my_real
168 . e,nu,bulk,g,unsa,du,g2,lamda,
169 . al1, al2 , al3 , al4 ,
170 . al5, al6 , al7 , al8 ,
171 . al9, al10, al11, al12,
172 . dav,
173 . lpp11,lpp12,lpp21,lpp22,lpp23,
174 . lpp34,lpp45,lpp56,lpp13,
175 . fct,fp1,fp2,fp3,fp4,fp5,
176 . fpp1,fpp2,fpp3,dfp1,dfp2,dfp3,fpp4,fpp5,
177 . dfpp1,dfpp2,dfpp3,df1,df2,df3,df,dfp4,dfp5,dfp6,
178 . epspi,dfpp4,dfpp5,dfpp6,
179 . aswift ,epso,qvoce,beta,ko,alpha,nexp,df4,df5,df6,
180 . expv,kswift,kvoce,unsp,unsc,pla_i,yld_i,
181 . dswift,dvoce,dfepsp,yldp,unspt,unsct,
182 . tpp1, tpp2,tpp3,tp1
183 my_real
184 . plap(nel),f_epsp(nel),
185 . f_epspt(nel)
186 my_real
187 . yfac(nel,2),pui_1
188 DATA nmax/5/
189c-----------------------------------------------
190c USER VARIABLES INITIALIZATION
191c-----------------------------------------------
192c Swift/Voce
193c k=alpha*k1(eps,p)+(1-alpha)*k2(eps,p)
194c k1=A*(eps,p+e0)**n
195c k2=Q*(1-exp(-b*eps,p))+sig0
196c-----------------------------------------------
197c UVAR1 =PLA
198c UVAR2 =YLD
199c UVAR3 =DPLA
200c UVAR4 =ETSE
201c UVAR5 =IPOS
202c UVAR6 =H
203c-----------------------------------------------
204C
205 yldp = zero
206 e = uparam(1)
207 nu = uparam(2)
208 g = uparam(23)
209 g2 = two *g
210 bulk = e/three/(one-two*nu)
211
212 al1 = uparam(3)
213 al2 = uparam(4)
214 al3 = uparam(5)
215 al4 = uparam(6)
216 al5 = uparam(7)
217 al6 = uparam(8)
218 al7 = uparam(9)
219 al8 = uparam(10)
220 al9 = uparam(11)
221 al10 = uparam(12)
222 al11 = uparam(13)
223 al12 = uparam(14)
224
225
226c
227 expa = nint(uparam(15) )
228 unsa = one/expa
229 expam2= expa - 2
230c
231 alpha = uparam(16)
232 nexp = uparam(17)
233 epso = uparam(18)
234 aswift = uparam(19)
235
236 nrate = nint(uparam(25))
237 lamda = uparam(25+2*nfunc+9)
238
239
240 qvoce = uparam(25+2*nrate+1)
241 beta = uparam(25+2*nrate+2)
242 ko = uparam(25+2*nrate+3)
243
244 unsp = uparam(25+2*nrate+4)
245 unsc = uparam(25+2*nrate+5)
246
247
248 unspt = uparam(25+2*nrate+6)
249 unsct = uparam(25+2*nrate+7)
250
251 iflagsr = nint(uparam(25+2*nrate+8))
252 !DEPLACER AU STARTER
253 lpp11=(-two* al3 +two*al4 + eight*al5- two*al6)/nine
254 lpp12=(-four*al4 +four*al6+ al3 - four*al5)/nine
255 lpp13=(al3 +two*al4 -four*al5- two*al6 )/nine
256 lpp21=(four*al3 -four*al4-four*al5+ al6)/nine
257 lpp22=(-two* al3+ eight*al4 + two*al5- two*al6)/nine
258 lpp23=(-two* al3 -four*al4 + two*al5+ al6)/nine
259 lpp34=al8
260 lpp45=al11
261 lpp56=al12
262
263
264 iflag = int(uparam(24))
265
266c ---
267 DO i=1,nel
268 dav = depsxx(i)+depsyy(i)+depszz(i)
269 signxx(i) = sigoxx(i) + depsxx(i)*g2 + lamda*dav
270 signyy(i) = sigoyy(i) + depsyy(i)*g2 + lamda*dav
271 signzz(i) = sigozz(i) + depszz(i)*g2 + lamda*dav
272 signxy(i) = sigoxy(i) + depsxy(i)*g
273 signyz(i) = sigoyz(i) + depsyz(i)*g
274 signzx(i) = sigozx(i) + depszx(i)*g
275
276C------------------
277C SIGMA TRIAL
278C------------------
279 sigtrxx(i)=signxx(i)
280 sigtryy(i)=signyy(i)
281 sigtrzz(i)=signzz(i)
282 sigtrxy(i)=signxy(i)
283 sigtryz(i)=signyz(i)
284 sigtrzx(i)=signzx(i)
285c-------------------
286c SOUND SPEED, TANGENT MODULUS
287c-------------------
288 soundsp(i) = sqrt((bulk+four_over_3*g)/rho0(i))
289 viscmax(i) = zero
290 etse(i) = one
291 dpla(i) = zero
292 ENDDO
293C-------------------
294c flag formulation ************************************************
295C-------------------
296
297 IF (iflag == 1) THEN ! ANALYTIC SWIFT-VOCE YIELD
298c ----- Yield ---------------
299 DO i=1,nel
300 expv = exp(-beta*pla(i))
301 IF((pla(i) + epso)>zero) THEN
302 pui_1 = exp(nexp*log((pla(i) + epso)))
303 ELSE
304 pui_1 = zero
305 ENDIF
306! KSWIFT = ASWIFT*(PLA(I) + EPSO)**NEXP
307 kswift = aswift*pui_1
308 kvoce = ko + qvoce*(one - expv)
309 yld(i) = alpha*kswift + (one-alpha)*kvoce
310 f_epspt(i) = one
311 f_epsp(i) = one
312 !!print*, 'yld ', YLD(I)
313 ENDDO
314 IF (unspt /= zero )THEN
315 DO i=1,nel
316 IF(epsp(i)/=zero)
317 . f_epspt(i) =one + exp(unspt * log(unsct*epsp(i)))
318 ENDDO
319 ENDIF
320 DO i=1,nel
321 !--------------------
322 !BARLAT CRITERION
323 !--------------------
324 !XP(5) = LP(5*6) * SIGMA(6)
325 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
326 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
327 xpxy(i) = al7*signxy(i)
328 xpyz(i) = al9*signyz(i)
329 xpzx(i) = al10*signzx(i)
330
331 !XPP(5) = LPP(5*6) * SIGMA(6)
332 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
333 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
334 xppxy(i) = lpp34*signxy(i)
335 xppyz(i) = lpp45*signyz(i)
336 xppzx(i) = lpp56*signzx(i)
337
338 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
339 phip(i) = sqrt(max(zero, deltap(i)))**expa
340
341
342
343 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
344 deltapp(i) = sqrt(max(zero, deltapp(i)))
345 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
346 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
347
348 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
349
350 !SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
351 IF(((phip(i)+phipp(i)))>zero) THEN
352 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
353 ELSE
354 sig(i) = zero
355 ENDIF
356
357 ENDDO
358
359C Check Yield Condition
360C ----------------------------------------------------
361 nindx = 0
362 DO i=1,nel
363 IF (sig(i)>yld(i).AND.off(i) == one)THEN ! Plastic Loading
364 nindx = nindx + 1
365 indx(nindx) = i
366 ENDIF
367 ENDDO
368 !--------------------
369 !PLASTIC FLOW
370 !--------------------
371 DO ii=1,nindx
372 i = indx(ii)
373 !IF (SIG(I)>YLD(I).AND.OFF(I) == ONE)THEN
374 dpla(i)= uvar(i,3)+ em09
375
376 !dpla(i)= (sig(i)-yld(i))/(three*g+dydx(i))
377 deplxx(i)=zero
378 deplyy(i)=zero
379 deplzz(i)=zero
380 deplxy(i)=zero
381 deplyz(i)=zero
382 deplzx(i)=zero
383
384 !YLD_I=YLD(I)
385 pla_i=pla(i)
386
387 DO k=1,4
388 !algo newton
389 !calcul de dsigma/ddpla=dphi/dsigma*dsigma/ddpla
390 !phi = phip + phipp
391 !dphip/dsigma = dphip/dxp*dxp/dsigma where dxp/dsigma = Lp
392 pla(i) = pla_i + dpla(i)
393 plap(i)= dpla(i)/timestep
394C
395 expv = exp(-beta*pla(i))
396! KSWIFT = ASWIFT*(PLA(I) + EPSO)**NEXP
397 IF((pla(i) + epso)>zero) THEN
398 pui_1 = exp(nexp*log(pla(i) + epso))
399 ELSE
400 pui_1 = zero
401 ENDIF
402 kswift = aswift*pui_1
403 kvoce = ko + qvoce*(one - expv)
404 yld(i) = alpha*kswift + (one-alpha)*kvoce
405
406 !===========ADD PL SR DEPENDENCY THROUGH COWPER SEYMONDS
407 !Compute plastic strain rate Cowper seymonds
408 IF( unsp/=zero)THEN
409 IF(plap(i)>zero)THEN
410 f_epsp(i) = one+exp(unsp * log(unsc*plap(i)))
411 ELSE
412 f_epsp(i) = one
413 ENDIF
414 ELSE
415 f_epsp(i) = one
416 ENDIF
417 !DERIVATIVE OF YLD
418 dswift = nexp*kswift/(pla(i) + epso)
419 dvoce = qvoce*beta*expv
420 !DFEPSP = UNSP *F_EPSP(I)/ (DPLA(I)+TIMESTEP/MAX(EM20,UNSC))
421 dfepsp = unsp *(f_epsp(i)-one)/ max(em20,dpla(i))
422 yldp = alpha*dswift + (one-alpha)*dvoce
423 yldp = yldp*f_epsp(i) + yld(i)*dfepsp
424 yld(i) = yld(i) * f_epsp(i) *f_epspt(i)
425 !===========END PL SR DEPENDENCY
426 fct = sig(i)- yld(i) ! residu
427
428 tp1 = expa*(sqrt(max(zero,deltap(i))))**(expa-2)
429 !D_phip/D_XP
430 fp1 = tp1 * (xpxx(i) - xpyy(i))
431 fp2 = -fp1
432 fp3 = three * tp1 * xpxy(i)
433 fp4 = three * tp1 * xpyz(i)
434 fp5 = three * tp1 * xpzx(i)
435
436
437 tpp1 = expa * deltapp1(i) ** (expa - 1)
438 tpp2 = expa * deltapp2(i) ** (expa - 1)
439 tpp3 = two *(tpp1 - tpp2) /deltapp(i)
440
441 !D_phipp/D_XPP
442 fpp1 = tpp1 *(three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
443 . + tpp2 *(three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
444 fpp2 = tpp1 *(-three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
445 . + tpp2 *(-three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
446 fpp3 = tpp3 * xppxy(i)
447 fpp4 = tpp3 * xppyz(i)
448 fpp5 = tpp3 * xppzx(i)
449
450 dfp1=third*(two*al1*fp1-al2*fp2)
451 dfp2=third*(two*al2*fp2-al1*fp1)
452 dfp3=third*(-al1*fp1-al2*fp2)
453 dfp4=al7*fp3
454 dfp5=al9*fp4
455 dfp6=al10*fp5
456
457 dfpp1=fpp1*lpp11+fpp2*lpp21
458 dfpp2=fpp1*lpp12+fpp2*lpp22
459 dfpp3=fpp1*lpp13+fpp2*lpp23
460 dfpp4=fpp3*lpp34
461 dfpp5=fpp4*lpp45
462 dfpp6=fpp5*lpp56
463
464 du = unsa * sig(i) / (max(em20, phip(i)+phipp(i)) )
465
466
467 df1 = dfp1+dfpp1
468 df2 = dfp2+dfpp2
469 df3 = dfp3+dfpp3
470 df4 = dfp4+dfpp4
471 df5 = dfp5+dfpp5
472 df6 = dfp6+dfpp6
473
474 !DF = DF1*DSDEPL1+DF2*DSDEPL2+DSDEPL3*DF3
475 ddeplxx(i)=df1*du
476 ddeplyy(i)=df2*du
477 ddeplzz(i)=df3*du
478 ddeplxy(i)=df4*du
479 ddeplyz(i)=df5*du
480 ddeplzx(i)=df6*du
481
482 df = - g * du**2
483 . * (two*(df1*df1 + df2*df2 + df3*df3 )
484 . + df4*df4 + df5*df5 + df6*df6 )
485
486c DF = -TWO*G * (DU * DF1) **2
487c . -TWO*G * (DU * DF2) **2
488c . -TWO*G * (DU * DF3) **2
489c . -g * (du * df4) **2
490c . -g * (du * df5) **2
491c . -g * (du * df6) **2
492
493 df = df - yldp
494 !CALCUL PAS PLASTIQUE
495 dpla(i)=dpla(i)-fct/df
496 dpla(i)=max(zero,dpla(i))
497
498 !TENSEUR PLASTIQUE
499 deplxx(i)=dpla(i)*df1*du
500 deplyy(i)=dpla(i)*df2*du
501 deplzz(i)=dpla(i)*df3*du
502 deplxy(i)=dpla(i)*df4*du
503 deplyz(i)=dpla(i)*df5*du
504 deplzx(i)=dpla(i)*df6*du
505
506 !update of the stresses
507 signxx(i)=sigtrxx(i)- g2*deplxx(i)
508 signyy(i)=sigtryy(i)- g2*deplyy(i)
509 signzz(i)=sigtrzz(i)- g2*deplzz(i)
510 signxy(i)=sigtrxy(i)- g*deplxy(i)
511 signyz(i)=sigtryz(i)- g*deplyz(i)
512 signzx(i)=sigtrzx(i)- g*deplzx(i)
513 !CALCUL NOUVEAU SIGMA EFFECTIF
514 !--------------------
515 !BARLAT CRITERION
516 !--------------------
517 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
518 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
519 xpxy(i) = al7*signxy(i)
520 xpyz(i) = al9*signyz(i)
521 xpzx(i) = al10*signzx(i)
522
523 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
524 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
525 xppxy(i) = lpp34*signxy(i)
526 xppyz(i) = lpp45*signyz(i)
527 xppzx(i) = lpp56*signzx(i)
528
529 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
530
531 phip(i) = sqrt(max(zero, deltap(i)))**expa
532
533
534 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
535 deltapp(i) = sqrt(max(em20, deltapp(i)))
536
537 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
538 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
539
540 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
541
542! SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
543 IF(((phip(i)+phipp(i)))>zero) THEN
544 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
545 ELSE
546 sig(i) = zero
547 ENDIF
548
549 ENDDO ! iter
550 uvar(i,3)=dpla(i)
551 ENDDO ! NINDX
552C-------------------
553C-------------------
554 ELSE ! IFLAG *TABULATED YIELD
555C-------------------
556c yield tabulated in below
557C-------------------
558 IF(iflagsr == 0)THEN
559C-------------------
560C-------------------
561 DO i=1,nel
562 jj(i) = 1
563 DO j=2,nrate-1
564 IF (epsp(i) > uparam(25+j)) jj(i) = j
565 ENDDO
566 ENDDO
567 DO i=1,nel
568 epspi = uparam(25+jj(i)) !
569 rate(i)=(epsp(i) - epspi)/(uparam(26+jj(i)) - epspi)
570 yfac(i,1)=uparam(25+nrate+jj(i))*facyldi(i)
571 yfac(i,2)=uparam(25+nrate+jj(i)+1)*facyldi(i)
572 ENDDO
573 DO i=1,nel
574 j1 = jj(i)
575 j2 = j1+1
576 ipos1(i) = nint(uvar(i,j1+5))
577 iad1(i) = npf(ipm(10+j1,mat(1))) / 2 + 1
578 ilen1(i) = npf(ipm(10+j1,mat(1))+1) / 2 - iad1(i)-ipos1(i)
579 ipos2(i) = nint(uvar(i,j2+5))
580 iad2(i) = npf(ipm(10+j2,mat(1))) / 2 + 1
581 ilen2(i) = npf(ipm(10+j2,mat(1))+1) / 2 - iad2(i)-ipos2(i)
582 END DO
583 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
584 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
585 DO i=1,nel
586 j1 = jj(i)
587 j2 = j1+1
588 y1(i)=y1(i)*yfac(i,1)
589 y2(i)=y2(i)*yfac(i,2)
590 yld(i) = (y1(i) + rate(i)*(y2(i)-y1(i)))
591 yld(i) = max(yld(i),em20)
592 dydx1(i)=dydx1(i)*yfac(i,1)
593 dydx2(i)=dydx2(i)*yfac(i,2)
594 dydx(i) = (dydx1(i) + rate(i)*(dydx2(i)-dydx1(i)))
595 uvar(i,j1+5) = ipos1(i)
596 uvar(i,j2+5) = ipos2(i)
597 ENDDO
598c------------------------------
599 ELSE !flagSR
600c------------------------------
601 IF (isigi == 0) THEN
602 IF(time == zero)THEN
603 DO i=1,nel
604 yfac(i,1)= uparam(25+nrate+1)*facyldi(i)
605 yld(i) = yfac(i,1)*finter(ipm(11,mat(1)) ,zero,npf,tf,dydx(i))
606 uvar(i,2)= yld(i)
607 ENDDO
608 ENDIF
609 ENDIF
610 DO i=1,nel
611 pla(i) = uvar(i,1)
612 yld(i) = uvar(i,2)
613 dydx(i) = uvar(i,5)
614 plap(i) = uvar(i,6)
615 jj(i) = 1
616 ENDDO
617 ENDIF !flagSR
618C-------------------
619C-------------------
620c------------------------------
621 DO i=1,nel
622 !--------------------
623 !BARLAT CRITERION
624 !--------------------
625 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
626 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
627 xpxy(i) = al7*signxy(i)
628 xpyz(i) = al9*signyz(i)
629 xpzx(i) = al10*signzx(i)
630
631 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
632 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
633 xppxy(i) = lpp34*signxy(i)
634 xppyz(i) = lpp45*signyz(i)
635 xppzx(i) = lpp56*signzx(i)
636
637 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
638 phip(i) = sqrt(max(zero, deltap(i)))**expa
639
640
641 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
642 deltapp(i) = sqrt(max(em20, deltapp(i)))
643 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
644 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
645
646 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
647
648! SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
649 IF((phip(i)+phipp(i))>zero) THEN
650 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
651 ELSE
652 sig(i) = zero
653 ENDIF
654 ENDDO
655C ----------------------------------------------------
656C Check Yield Condition
657C ----------------------------------------------------
658 nindx = 0
659 DO i=1,nel
660 IF (sig(i)>yld(i).AND.off(i) == one)THEN ! Plastic Loading
661 nindx = nindx + 1
662 indx(nindx) = i
663 ENDIF
664 ENDDO
665C ----------------------------------------------------
666
667C ----------------------------------------------------
668c flag strain rate option *****************************************
669C ----------------------------------------------------
670 IF(iflagsr == 0)THEN
671 DO ii=1,nindx
672 i = indx(ii)
673 !--------------------
674 !PLASTIC FLOW
675 !--------------------
676 !DPLA(I)= UVAR(I,3)
677 !DPLA(I)= (SIG(I)-YLD(I))/(THREE*G+DYDX(I))
678 dpla(i)= uvar(i,3)+ em09
679 deplxx(i)=zero
680 deplyy(i)=zero
681 deplzz(i)=zero
682 deplxy(i)=zero
683 deplyz(i)=zero
684 deplzx(i)=zero
685
686 yld_i=yld(i)
687 DO k=1,5
688 !algo newton
689 !calcul de dsigma/ddpla=dphi/dsigma*dsigma/ddpla
690 !phi = phip + phipp
691 !dphip/dsigma = dphip/dxp*dxp/dsigma where dxp/dsigma = Lp
692 yld(i) =yld_i+dydx(i)*dpla(i)
693 fct = sig(i)- yld(i) ! residu
694
695 tp1 = expa*(sqrt(max(zero,deltap(i))))**(expa-2)
696
697 !D_phip/D_XP
698 fp1 = tp1 * (xpxx(i) - xpyy(i))
699 fp2 = -fp1
700 fp3 = three * tp1 * xpxy(i)
701 fp4 = three * tp1 * xpyz(i)
702 fp5 = three * tp1 * xpzx(i)
703
704
705 tpp1 = expa * deltapp1(i) ** (expa - 1)
706 tpp2 = expa * deltapp2(i) ** (expa - 1)
707 tpp3 = two *(tpp1 - tpp2) /deltapp(i)
708
709 !D_phipp/D_XPP
710 fpp1 = tpp1 *(three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
711 . +tpp2 *(three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
712 fpp2 = tpp1 *(-three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
713 . +tpp2 *(-three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
714 fpp3 = tpp3 * xppxy(i)
715 fpp4 = tpp3 * xppyz(i)
716 fpp5 = tpp3 * xppzx(i)
717
718 dfp1=third*(two*al1*fp1-al2*fp2)
719 dfp2=third*(two*al2*fp2-al1*fp1)
720 dfp3=third*(-al1*fp1-al2*fp2)
721 dfp4=al7*fp3
722 dfp5=al9*fp4
723 dfp6=al10*fp5
724
725
726 dfpp1=fpp1*lpp11+fpp2*lpp21
727 dfpp2=fpp1*lpp12+fpp2*lpp22
728 dfpp3=fpp1*lpp13+fpp2*lpp23
729 dfpp4=fpp3*lpp34
730 dfpp5=fpp4*lpp45
731 dfpp6=fpp5*lpp56
732
733 du = unsa * sig(i) / (max(em20, phip(i)+phipp(i)) )
734
735 df1 = dfp1+dfpp1
736 df2 = dfp2+dfpp2
737 df3 = dfp3+dfpp3
738 df4 = dfp4+dfpp4
739 df5 = dfp5+dfpp5
740 df6 = dfp6+dfpp6
741
742 df = -two*g * (du * df1) **2
743 . -two*g * (du * df2) **2
744 . -two*g * (du * df3) **2
745 . -g * (du * df4) **2
746 . -g * (du * df5) **2
747 . -g * (du * df6) **2
748
749 df = df - dydx(i)
750
751 !CALCUL PAS PLASTIQUE
752 dpla(i)=dpla(i)-fct/df
753 dpla(i)=max(zero,dpla(i))
754 !TENSEUR PLASTIQUE
755 deplxx(i)=dpla(i)*df1*du
756 deplyy(i)=dpla(i)*df2*du
757 deplzz(i)=dpla(i)*df3*du
758 deplxy(i)=dpla(i)*df4*du
759 deplyz(i)=dpla(i)*df5*du
760 deplzx(i)=dpla(i)*df6*du
761
762 !update of the stresses
763 signxx(i)=sigtrxx(i)- g2*deplxx(i)
764 signyy(i)=sigtryy(i)- g2*deplyy(i)
765 signzz(i)=sigtrzz(i)- g2*deplzz(i)
766 signxy(i)=sigtrxy(i)- g*deplxy(i)
767 signyz(i)=sigtryz(i)- g*deplyz(i)
768 signzx(i)=sigtrzx(i)- g*deplzx(i)
769
770 !CALCUL NOUVEAU SIGMA EFFECTIF
771 !--------------------
772 !BARLAT CRITERION
773 !--------------------
774 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
775 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
776 xpxy(i) = al7*signxy(i)
777 xpyz(i) = al9*signyz(i)
778 xpzx(i) = al10*signzx(i)
779
780 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
781 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
782 xppxy(i) = lpp34*signxy(i)
783 xppyz(i) = lpp45*signyz(i)
784 xppzx(i) = lpp56*signzx(i)
785
786 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
787 phip(i) = sqrt(max(zero, deltap(i)))**expa
788
789
790 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
791 deltapp(i) = sqrt(max(em20, deltapp(i)))
792 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
793 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
794
795 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
796
797! !SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
798 IF((half*(phip(i)+phipp(i)))>zero) THEN
799 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
800 ELSE
801 sig(i) = zero
802 ENDIF
803
804 ENDDO !iter
805 uvar(i,3)=dpla(i)
806 ENDDO ! NINDX
807C=======================================================================
808 ELSE !IFLAGSR PLASTIC STRAIN RATE ******************************
809C=======================================================================
810 DO ii=1,nindx
811 i = indx(ii)
812 !--------------------
813 !PLASTIC FLOW
814 !--------------------
815 dpla(i) = uvar(i,3) + em09
816 plap(i) = dpla(i)/timestep
817 !--------------
818 !update yield :
819 !--------------
820 jj(i) = 1
821 DO j=2,nrate-1
822 IF(plap(i) >= uparam(25+j)) jj(i) = j
823 ENDDO
824 j1 = jj(i)
825 j2 = j1+1
826 rate(i)=(plap(i) - uparam(25+j1))/(uparam(26+j1) - uparam(25+j1))
827 yfac(i,1)=uparam(25+nrate+j1)*facyldi(i)
828 yfac(i,2)=uparam(25+nrate+j2)*facyldi(i)
829 !Y1 and Y2 for Pl SR dependency
830 y1(i) = yfac(i,1)*finter(ipm(10+j1,mat(1)),pla(i),npf,tf,dydx1(i))
831 y2(i) = yfac(i,2)*finter(ipm(10+j2,mat(1)),pla(i),npf,tf,dydx2(i))
832 yld(i) = y1(i) + rate(i)*(y2(i)-y1(i))
833 yld(i) = max(yld(i),em20)
834 dydx1(i)= dydx1(i)*yfac(i,1)
835 dydx2(i)= dydx2(i)*yfac(i,2)
836 dydx(i) = dydx1(i) + rate(i)*(dydx2(i)-dydx1(i))
837 yldp = dydx(i) + (y2(i)-y1(i))/
838 . (uparam(26+j1)-uparam(25+j1)) /timestep
839
840
841
842 !DPLA(I)= UVAR(I,3)
843 !DPLA(I)= (SIG(I)-YLD(I))/(THREE*G+DYDX(I))
844 deplxx(i)=zero
845 deplyy(i)=zero
846 deplzz(i)=zero
847 deplxy(i)=zero
848 deplyz(i)=zero
849 deplzx(i)=zero
850
851
852 DO k=1,5
853 !algo newton
854 !calcul de dsigma/ddpla=dphi/dsigma*dsigma/ddpla
855 !phi = phip + phipp
856 !dphip/dsigma = dphip/dxp*dxp/dsigma where dxp/dsigma = Lp
857 !---------------------------------------
858 !update yield :
859 !---------------------------------------
860 pla(i) = uvar(i,1) + dpla(i)
861 plap(i) = dpla(i) / timestep
862 jj(i) = 1
863 DO j=2,nrate-1
864 IF(plap(i) >= uparam(25+j)) jj(i) = j
865 ENDDO
866 j1 = jj(i)
867 j2 = j1+1
868 rate(i)=(plap(i) - uparam(25+j1))/(uparam(26+j1) - uparam(25+j1))
869 yfac(i,1)=uparam(25+nrate+j1)*facyldi(i)
870 yfac(i,2)=uparam(25+nrate+j2)*facyldi(i)
871 !Y1 and Y2 for Pl SR dependency
872 y1(i) = yfac(i,1)*finter(ipm(10+j1,mat(1)),pla(i),npf,tf,dydx1(i))
873 y2(i) = yfac(i,2)*finter(ipm(10+j2,mat(1)),pla(i),npf,tf,dydx2(i))
874 yld(i) = y1(i) + rate(i)*(y2(i)-y1(i))
875 yld(i) = max(yld(i),em20)
876 dydx1(i)=dydx1(i)*yfac(i,1)
877 dydx2(i)=dydx2(i)*yfac(i,2)
878 dydx(i) = (dydx1(i) + rate(i)*(dydx2(i)-dydx1(i)))
879 !---------------------------------------
880 yldp = dydx(i)+ (y2(i)-y1(i))
881 . /(uparam(26+j1)-uparam(25+j1)) /timestep
882
883 !---------------------------------------
884
885 !YLD(I) =YLD_I+DYDX(I)*DPLA(I)
886
887 fct = sig(i)- yld(i)
888
889 tp1 = expa*(sqrt(max(zero,deltap(i))))**(expa-2)
890
891 !D_phip/D_XP
892 fp1 = tp1 * (xpxx(i) - xpyy(i))
893 fp2 = -fp1
894 fp3 = three * tp1 * xpxy(i)
895 fp4 = three * tp1 * xpyz(i)
896 fp5 = three * tp1 * xpzx(i)
897
898
899 tpp1 = expa * deltapp1(i) ** (expa - 1)
900 tpp2 = expa * deltapp2(i) ** (expa - 1)
901 tpp3 = two *(tpp1 - tpp2) /deltapp(i)
902
903 !D_phipp/D_XPP
904 fpp1 = tpp1 *(three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
905 . +tpp2 *(three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
906 fpp2 = tpp1 *(-three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
907 . +tpp2 *(-three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
908 fpp3 = tpp3 * xppxy(i)
909 fpp4 = tpp3 * xppyz(i)
910 fpp5 = tpp3 * xppzx(i)
911
912 dfp1=third*(two*al1*fp1-al2*fp2)
913 dfp2=third*(two*al2*fp2-al1*fp1)
914 dfp3=third*(-al1*fp1-al2*fp2)
915 dfp4=al7*fp3
916 dfp5=al9*fp4
917 dfp6=al10*fp5
918
919
920 dfpp1=fpp1*lpp11+fpp2*lpp21
921 dfpp2=fpp1*lpp12+fpp2*lpp22
922 dfpp3=fpp1*lpp13+fpp2*lpp23
923 dfpp4=fpp3*lpp34
924 dfpp5=fpp4*lpp45
925 dfpp6=fpp5*lpp56
926
927 du = unsa * sig(i) / (max(em20, phip(i)+phipp(i)) )
928
929 df1 = dfp1+dfpp1
930 df2 = dfp2+dfpp2
931 df3 = dfp3+dfpp3
932 df4 = dfp4+dfpp4
933 df5 = dfp5+dfpp5
934 df6 = dfp6+dfpp6
935
936 df = -two*g * (du * df1) **2
937 . -two*g * (du * df2) **2
938 . -two*g * (du * df3) **2
939 . -g * (du * df4) **2
940 . -g * (du * df5) **2
941 . -g * (du * df6) **2
942
943 df = df - yldp
944
945 !CALCUL PAS PLASTIQUE
946 dpla(i)=dpla(i)-fct/df
947 dpla(i)=max(zero,dpla(i))
948 !TENSEUR PLASTIQUE
949 deplxx(i)=dpla(i)*df1*du
950 deplyy(i)=dpla(i)*df2*du
951 deplzz(i)=dpla(i)*df3*du
952 deplxy(i)=dpla(i)*df4*du
953 deplyz(i)=dpla(i)*df5*du
954 deplzx(i)=dpla(i)*df6*du
955
956 !update of the stresses
957 signxx(i)=sigtrxx(i)- g2*deplxx(i)
958 signyy(i)=sigtryy(i)- g2*deplyy(i)
959 signzz(i)=sigtrzz(i)- g2*deplzz(i)
960 signxy(i)=sigtrxy(i)- g*deplxy(i)
961 signyz(i)=sigtryz(i)- g*deplyz(i)
962 signzx(i)=sigtrzx(i)- g*deplzx(i)
963
964 !CALCUL NOUVEAU SIGMA EFFECTIF
965 !--------------------
966 !BARLAT CRITERION
967 !--------------------
968 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
969 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
970 xpxy(i) = al7*signxy(i)
971 xpyz(i) = al9*signyz(i)
972 xpzx(i) = al10*signzx(i)
973
974 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
975 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
976 xppxy(i) = lpp34*signxy(i)
977 xppyz(i) = lpp45*signyz(i)
978 xppzx(i) = lpp56*signzx(i)
979
980 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
981 phip(i) = sqrt(max(zero, deltap(i)))**expa
982
983
984 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
985 deltapp(i) = sqrt(max(em20, deltapp(i)))
986 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
987 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
988
989 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
990
991! !SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
992 IF((half*(phip(i)+phipp(i)))>zero) THEN
993 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
994 ELSE
995 sig(i) = zero
996 ENDIF
997 ENDDO ! iter k
998 uvar(i,3) =dpla(i)
999 ENDDO !ndx
1000c----------------------------
1001 ENDIF !flag formulation
1002 ENDIF
1003C=======================================================================
1004 DO i=1,nel
1005C
1006 IF( dpla(i) > zero)THEN
1007 h(i)=(yld(i)-uvar(i,2))/dpla(i)
1008 h(i)=max(em10,h(i))
1009 uvar(i,2)=yld(i)
1010 uvar(i,4)=h(i)/g2!(H(I)+E)
1011 uvar(i,5)=h(i)
1012 ENDIF
1013 pla(i)=uvar(i,1) + dpla(i)
1014 uvar(i,1) = pla(i)
1015 etse(i)= uvar(i,4)
1016 yld(i) = uvar(i,2)
1017 h(i) = uvar(i,5)
1018 ENDDO
1019 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
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:73