OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps52.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"
#include "scr17_c.inc"
#include "units_c.inc"
#include "com08_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps52 (nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, l_dmg, dmg, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, ngl, ipt, ipm, mat, epsp, ipla, sigy, defp, table)

Function/Subroutine Documentation

◆ sigeps52()

subroutine sigeps52 ( integer nel,
integer nuparam,
integer nuvar,
integer mfunc,
integer, dimension(mfunc) kfunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
rho0,
rho,
volume,
eint,
integer, intent(in) l_dmg,
intent(inout) dmg,
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,
integer, dimension(nel) ngl,
integer ipt,
integer, dimension(npropmi,*) ipm,
integer, dimension(nel) mat,
epsp,
integer ipla,
sigy,
defp,
type(ttable), dimension(*) table )

Definition at line 36 of file sigeps52.F.

49C-----------------------------------------------
50 USE table_mod
51 USE interface_table_mod
52C-----------------------------------------------
53C I m p l i c i t T y p e s
54C-----------------------------------------------
55#include "implicit_f.inc"
56C-----------------------------------------------
57C G l o b a l P a r a m e t e r s
58C-----------------------------------------------
59#include "mvsiz_p.inc"
60C---------+---------+---+---+--------------------------------------------
61C VAR | SIZE |TYP| RW| DEFINITION
62C---------+---------+---+---+--------------------------------------------
63C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
64C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
65C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
66C---------+---------+---+---+--------------------------------------------
67C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
68C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
69C NPF | * | I | R | FUNCTION ARRAY
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 RHO | NEL | F | R | DENSITY
77C VOLUME | NEL | F | R | VOLUME
78C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
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 UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
102C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
103C---------+---------+---+---+--------------------------------------------
104#include "param_c.inc"
105#include "scr17_c.inc"
106#include "units_c.inc"
107#include "com08_c.inc"
108C
109 INTEGER NEL, NUPARAM, NUVAR,IPT,
110 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*),IPLA, ITER, IFLAG
111 my_real
112 . time,timestep,uparam(*),
113 . rho(nel),rho0(nel),volume(nel),eint(nel),
114 . epspxx(nel),epspyy(nel),epspzz(nel),
115 . epspxy(nel),epspyz(nel),epspzx(nel),
116 . depsxx(nel),depsyy(nel),depszz(nel),
117 . depsxy(nel),depsyz(nel),depszx(nel),
118 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,
119 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
120 . sigoxx(nel),sigoyy(nel),sigozz(nel),
121 . sigoxy(nel),sigoyz(nel),sigozx(nel),
122 . epsp(nel)
123C
124 TYPE(TTABLE) TABLE(*)
125 INTEGER, INTENT(IN) :: L_DMG
126C-----------------------------------------------
127C O U T P U T A r g u m e n t s
128C-----------------------------------------------
129 my_real
130 . signxx(nel),signyy(nel),signzz(nel),
131 . signxy(nel),signyz(nel),signzx(nel),
132 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
133 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
134 . soundsp(nel),viscmax(nel),sigy(nel),defp(nel)
135C-----------------------------------------------
136C I N P U T O U T P U T A r g u m e n t s
137C-----------------------------------------------
138 my_real
139 . uvar(nel,nuvar), off(nel)
140 my_real, DIMENSION(NEL,L_DMG), INTENT(INOUT) :: dmg
141C-----------------------------------------------
142C VARIABLES FOR FUNCTION INTERPOLATION
143C-----------------------------------------------
144 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
145 my_real
146 . finter2, tf(*)
147 EXTERNAL finter2
148
149C-----------------------------------------------
150C L o c a l V a r i a b l e s
151C-----------------------------------------------
152 INTEGER I,J,IADBUFV,J1,J2,JJ(MVSIZ),NFUNC,
153 . NRATE(MVSIZ),IPOS1(MVSIZ),IPOS2(MVSIZ),IAD1(MVSIZ),
154 . ILEN1(MVSIZ),IAD2(MVSIZ),ILEN2(MVSIZ),
155 . NFUNCV(MVSIZ),PFUN(MVSIZ),
156 . IPOSP(MVSIZ),IADP(MVSIZ),ILENP(MVSIZ),NFUNCM,NRATEM,
157 . IPFLAG,IPARAM,NPAR,IYEILD_TAB,IADBUF,
158 , ITAB(100),NTABLE,NXH,NDIM,IPOS,NXK
159 my_real
160 . f(mvsiz), sigm(mvsiz), epsm(mvsiz), epsp1(mvsiz),
161 . lamda(mvsiz), fstar(mvsiz),p0(mvsiz),pn(mvsiz),
162 . et,yeild0(mvsiz),csd,
163 . visp(mvsiz), q1, q2, q3,
164 . sn, epsn, fi,fn,
165 . fc,ff,fu, n,
166 . sigxx(mvsiz),sigyy(mvsiz),sigxy(mvsiz),
167 . sigzx(mvsiz),sigyz(mvsiz), sigzz(mvsiz),
168 . sganxx(mvsiz), sganyy(mvsiz), sganzz(mvsiz),
169 . sganxy(mvsiz), sganxz(mvsiz), sganyz(mvsiz)
170 my_real
171 . e,nu,p,dav,vm,r,fac,epst,ep1,ep2,ep3,ep4,ep5,ep6,
172 . e1,e2,e3,e4,e5,e6,c,cc,d,y,yp,e42,e52,e62,epst2,
173 . dp11,dp22,dp33,dp12,dp13,dp23,a1,a2,
174 . d11,d22,d33,d12,d13,d23,dcrf,dcrm, dsepp,dcd,lam1,
175 . c1(mvsiz),c2(mvsiz) ,g(mvsiz),df,coh,sih,va,crit,conv,
176 . fg(mvsiz),fn1(mvsiz),yld,c11,p1, var, dtinv, sqr22,
177 . sxx,syy,szz, sigm1, sigm2, va1, vm1, va11,a21,yp1,yp2,
178 . df1, df2, xfac,yfac,dx2,xx(2),dft,yy
179C-----------------------------------------------
180C USER VARIABLES INITIALIZATION
181C-----------------------------------------------
182 ntable = ipm(226,mat(1))
183 iyeild_tab = 0
184 IF(ntable > 0) THEN
185 iyeild_tab = 1
186 itab(1)=ipm(226+1,mat(1))
187 iadbuf = ipm(7,mat(1))-1
188 xfac = uparam(iadbuf + 22 )
189 yfac = uparam(iadbuf + 23 )
190 ENDIF
191C
192 iadbufv = ipm(7,mat(1))-1
193 e = uparam(iadbufv+1)
194 nu = uparam(iadbufv+2)
195 iflag = uparam(iadbufv+18)
196 et = uparam(iadbufv+4)
197 n = uparam(iadbufv+5)
198 csd = one/uparam(iadbufv+6)
199 q1 = uparam(iadbufv+8)
200 q2 = uparam(iadbufv+9)
201 q3 = uparam(iadbufv+10)
202 sn = uparam(iadbufv+11)
203 IF(sn==zero)sn = ep20
204 epsn = uparam(iadbufv+12)
205 fi = uparam(iadbufv+13)
206 fn = uparam(iadbufv+14)
207 fc = uparam(iadbufv+15)
208 ff = uparam(iadbufv+16)
209 fu = uparam(iadbufv+17)
210 DO i=1,nel
211
212 yeild0(i)= uparam(iadbufv+3)
213 visp(i) = uparam(iadbufv+7)
214
215 g(i) = half*e/(one + nu)
216 c11 = e/three/(one- two*nu)
217 soundsp(i) = sqrt((c11 + four_over_3*g(i))/rho0(i))
218 viscmax(i) = zero
219 c1(i) = e*(one-nu) /((one + nu)*(one - two*nu))
220 c2(i) = c1(i)*nu/(one - nu)
221 ENDDO
222C
223 IF(time==0.0)THEN
224 DO i=1,nel
225 epsp1(i)=zero
226 epsm(i)=zero
227 IF(n<one)epsm(i)=em20
228 uvar(i,1)=epsm(i)
229 uvar(i,2)=yeild0(i)
230 dmg(i,4) = fi
231 dmg(i,5) = fi
232 sigy(i) =yeild0(i)
233 ENDDO
234C
235 IF(iyeild_tab > 0) THEN
236 DO i=1,nel
237 xx(1)=zero
238 xx(2)=zero
239 CALL table_interp (table(itab(1)),xx,yy)
240 yeild0(i) = yy *yfac
241 uvar(i,2) = yeild0(i)
242 sigy(i) = yeild0(i)
243 ENDDO
244 ENDIF
245 ENDIF
246C
247 sqr22= one/sqrt(two*pi)
248C
249C-----------------------------------------------
250C
251 DO i=1,nel
252 epsm(i) = uvar(i,1)
253 sigm(i) = uvar(i,2)
254 fg(i) = dmg(i,2)
255 fn1(i) = dmg(i,3)
256 f(i) = dmg(i,4)
257 fstar(i) = dmg(i,5)
258 epsp1(i)= zero
259C
260 signxx(i)=sigoxx(i) + c1(i)*depsxx(i)
261 . + c2(i)*(depsyy(i) + depszz(i))
262 signyy(i)=sigoyy(i) + c1(i)*depsyy(i)
263 . + c2(i)*(depsxx(i) + depszz(i))
264 signzz(i)=sigozz(i) + c1(i)*depszz(i)
265 . + c2(i)*(depsxx(i) + depsyy(i))
266 signxy(i)=sigoxy(i) + g(i)*depsxy(i)
267 signyz(i)=sigoyz(i) + g(i)*depsyz(i)
268 signzx(i)=sigozx(i) + g(i)*depszx(i)
269 IF(off(i)==zero) THEN
270 signxx(i)=zero
271 signyy(i)=zero
272 signzz(i)=zero
273 signxy(i)=zero
274 signzx(i)=zero
275 signyz(i)=zero
276 ENDIF
277 ENDDO
278C
279 dtinv=timestep/max(timestep**2,em20)
280C-------------------
281C STRAIN RATE
282C-------------------
283C VISCOPLASTIC
284 IF(iflag==0)THEN
285C-------------
286C IFLAG = 0,
287C-----------
288 DO i=1,nel
289 IF(off(i)==one)THEN
290 IF(f(i)<=fc)THEN
291 fstar(i) = f(i)
292 df = one
293 ELSE
294 df = (fu- fc)/(ff-fc)
295 fstar(i)= fc + df*(f(i)-fc)
296 ENDIF
297C
298 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
299 sxx = signxx(i)-pn(i)
300 syy = signyy(i)-pn(i)
301 szz = signzz(i)-pn(i)
302 vm = half*(sxx**2 + syy**2 + szz**2)
303 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
304 vm = sqrt(three*vm)
305 sigm1 = one/max(sigm(i),em20)
306 var = three_half * q2 * pn(i)* sigm1
307 var = exp(var)
308 coh = half * (var + one/max(em20,var))
309 sih = half * (var - one/max(em20,var))
310 va = one + q3 * fstar(i)**2 - two * q1 * fstar(i) * coh
311 va= sqrt(max(zero,va))
312 yld = sigm(i) * va
313 sigy(i) = yld
314 crit = vm - yld
315 IF(vm < yld .OR. sigy(i) == zero )THEN
316 lamda(i)=zero
317 ELSE
318 lamda(i) = zero
319 sigxx(i) = signxx(i)
320 sigyy(i) = signyy(i)
321 sigzz(i) = signzz(i)
322 sigxy(i) = signxy(i)
323 sigzx(i) = signzx(i)
324 sigyz(i) = signyz(i)
325 dp11 = zero
326 dp22 = zero
327 dp33 = zero
328 dp12 = zero
329 dp13 = zero
330 dp23 = zero
331 IF(f(i)==one)THEN
332 a21 = ep20
333 ELSE
334 a21 = sigm1 /(one - f(i))
335 ENDIF
336 a1 = fn*exp(-half*((epsm(i)-epsn)/sn)**2)/sn
337 a1 = a1*sqr22
338 DO iter=1,5
339 vm1 = one/max(vm,em20)
340 va1 = one/max(va,em20)
341 va11 =half*q1*q2*fstar(i)*sih*va1
342 d11 = half*(two*signxx(i) - signyy(i)-signzz(i))*vm1
343 . + va11
344 d22 = half*(two*signyy(i) - signxx(i)-signzz(i))*vm1
345 . + va11
346 d33 = half*(two*signzz(i) - signxx(i)-signyy(i))*vm1
347 . + va11
348 d12 = three*signxy(i)*vm1
349 d13 = three*signzx(i)*vm1
350 d23 = three*signyz(i)*vm1
351C
352 a2 = d11*signxx(i) + d22*signyy(i) + d33*signzz(i) +
353 . two*(d12*signxy(i) + d13*signzx(i) + d23*signyz(i))
354 a2 = a2*a21
355C
356 dcrf=-sigm(i)*(q3*fstar(i)*df-q1*coh*df)*va1
357 dcrm= - va - three*va11*pn(i)*sigm1
358 IF(n==1)THEN
359 dsepp = et*(one+ ( epsp(i)*csd)**visp(i))
360 ELSE
361 dsepp = et * n * epsm(i) ** (n - 1)*
362 . (one + ( epsp(i)*csd)**visp(i))
363 ENDIF
364C
365 dcd = c1(i)*(d11**2 + d22**2 + d33**2) +
366 . two*c2(i)*(d11*d22 + d11*d33 + d22*d33) +
367 . two*g(i)*(d12**2 + d13**2 + d23**2)
368C
369 lam1= dcd - dcrm* dsepp*a2 -
370 . dcrf*((one - f(i))*(d11+d22+d33) + a1*a2)
371 IF(lam1/=zero) lamda(i) = max(zero,(vm - yld))/lam1
372C PLASTIC DEFORMATION
373 dp11 = dp11 + lamda(i) * d11
374 dp22 = dp22 + lamda(i) * d22
375 dp33 = dp33 + lamda(i) * d33
376 dp12 = dp12 + lamda(i) * d12
377 dp13 = dp13 + lamda(i) * d13
378 dp23 = dp23 + lamda(i) * d23
379C NEW STRESS
380 signxx(i)= sigxx(i) - c1(i)*dp11 - c2(i)*(dp22 + dp33)
381 signyy(i)= sigyy(i) - c1(i)*dp22 - c2(i)*(dp11 + dp33)
382 signzz(i)= sigzz(i) - c1(i)*dp33 - c2(i)*(dp22 + dp11)
383 signxy(i)= sigxy(i) - two*g(i)*dp12
384 signzx(i)= sigzx(i) - two*g(i)*dp13
385 signyz(i)= sigyz(i) - two*g(i)*dp23
386C....... AMISSIBLE PLASTIC DEFORMATION.
387
388 epsp1(i)= signxx(i)*dp11 + signyy(i)*dp22 + signzz(i)*dp33
389 . + two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
390 epsp1(i)=epsp1(i)*a21
391 epsp1(i) = max(zero, epsp1(i))
392 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
393 sxx = signxx(i)-pn(i)
394 syy = signyy(i)-pn(i)
395 szz = signzz(i)-pn(i)
396 vm = half*(sxx**2 + syy**2 + szz**2)
397 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
398 vm = sqrt(three*vm)
399 var = three_half * q2*pn(i) *sigm1
400 var = exp( var)
401 coh = half * ( var + one/max(em20,var))
402 sih = half * ( var - one/max(em20,var))
403 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
404 va = sqrt(max(zero,va))
405 yld = sigm(i) * va
406C IF(CONV<1E-3)GO TO 100
407 ENDDO ! ITER
408C
409 fg(i) = fg(i) + (one-f(i))*(dp11+dp22+dp33)
410 sigy(i) = yld
411 fn1(i) = fn1(i) + a1*epsp1(i)
412 epsm(i) = epsm(i) + epsp1(i)
413 f(i) = fi + fg(i) + fn1(i)
414 IF(q1==zero.and.q2==zero.and.q3==zero)f(i)=zero
415 IF(f(i)<zero) f(i)=zero
416 sigm(i) = (yeild0(i) + et * epsm(i) ** n)*
417 . (one+ (epsp(i)*csd) ** visp(i))
418C....COMPUTE F*.....
419 IF(f(i)<=fc)THEN
420 fstar(i) = f(i)
421 ELSE
422 IF(ff==fc)THEN
423 fstar(i) = ep20
424 ELSE
425 fstar(i) = fc+ (fu-fc)/(ff-fc)*(f(i)-fc)
426 ENDIF
427 ENDIF
428 uvar(i,1) = epsm(i)
429 uvar(i,2) = sigm(i)
430 dmg(i,1) = f(i)/ff
431 dmg(i,2) = fg(i)
432 dmg(i,3) = fn1(i)
433 dmg(i,4) = f(i)
434 dmg(i,5) = fstar(i)
435 ENDIF
436 IF(off(i)==one.AND.fstar(i)>=ff)THEN
437 off(i)=four_over_5
438 WRITE(iout, 1000) ngl(i)
439 WRITE(istdo,1100) ngl(i),tt
440 ENDIF
441 ELSE
442 IF(off(i)<em01)off(i)=zero
443 off(i)=off(i)*four_over_5
444 ENDIF
445 ENDDO
446 ELSEIF(iflag==1)THEN
447C----------
448C IFLAG set to 1
449C..........
450 DO i=1,nel
451 IF(off(i)==one)THEN
452C
453 IF(f(i)<=fc)THEN
454 fstar(i) = f(i)
455 df = one
456 ELSE
457 fstar(i) = fc + (fu - fc)/(ff-fc)*(f(i)-fc)
458 df = (fu - fc)/(ff - fc)
459 fstar(i) = fc + df*(f(i)-fc)
460 ENDIF
461 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
462 sxx = signxx(i)-pn(i)
463 syy = signyy(i)-pn(i)
464 szz = signzz(i)-pn(i)
465 vm = half*(sxx**2 + syy**2 + szz**2)
466 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
467 vm = three*vm
468 sigm1 = one/max(em20,sigm(i))
469 sigm2 = sigm1**2
470 var = three_half * q2 * pn(i) * sigm1
471 var = exp(var)
472 coh = half*(var + one/max(em20,var))
473 sih = half*(var - one/max(em20,var))
474 IF(pn(i)<=zero) THEN
475 va = one +q3*fstar(i)**2 - two*q1*fstar(i)
476 ELSE
477 va= one +q3*fstar(i)**2 - two*q1*fstar(i)*coh
478 ENDIF
479C
480 yld = sigm(i)
481 crit = vm * sigm2 - va
482C
483 sigy(i) = yld*sqrt(max(zero, va))
484 IF(crit < zero .OR. sigy(i) == zero)THEN
485 lamda(i) = zero
486 ELSE
487C
488 lamda(i) = zero
489 sigxx(i) = signxx(i)
490 sigyy(i) = signyy(i)
491 sigzz(i) = signzz(i)
492 sigxy(i) = signxy(i)
493 sigzx(i) = signzx(i)
494 sigyz(i) = signyz(i)
495C
496 IF(f(i)==one)THEN
497 a21 = ep20
498 ELSE
499 a21 = sigm1 /(one - f(i))
500 ENDIF
501 a1 = fn*exp(-half*((epsm(i)-epsn)/sn)**2)/sn
502 a1 = a1*sqr22
503 dp11 = zero
504 dp22 = zero
505 dp33 = zero
506 dp12 = zero
507 dp13 = zero
508 dp23 = zero
509 DO iter=1,5
510 IF(pn(i)<=zero)THEN
511 d11 = (two*signxx(i) - signyy(i) - signzz(i))*sigm2
512 d22 = (two*signyy(i) - signxx(i) - signzz(i))*sigm2
513 d33 = (two*signzz(i) - signxx(i) - signyy(i))*sigm2
514 d12 = six*signxy(i)*sigm2
515 d13 = six*signzx(i)*sigm2
516 d23 = six*signyz(i)*sigm2
517 dcrm = -two * vm * sigm2 * sigm1
518 dcrf = two * q1 * df - two * q3 * df * fstar(i)
519 ELSE
520 va1 = q1 * q2 * fstar(i) * sih*sigm1
521 d11 = (two*signxx(i) - signyy(i) - signzz(i))*sigm2 + va1
522 d22 = (two*signyy(i) - signxx(i) - signzz(i))*sigm2 + va1
523 d33 = (two*signzz(i) - signxx(i) - signyy(i))*sigm2 + va1
524 d12 = six * signxy(i)*sigm2
525 d13 = six * signzx(i)*sigm2
526 d23 = six * signyz(i)*sigm2
527 dcrm =-two* vm*sigm1*sigm2 - three * va1 * pn(i)*sigm1
528 dcrf = two * q1*df*coh - two*q3*df*fstar(i)
529 ENDIF
530C
531 a2 = d11 * signxx(i) + d22 * signyy(i) + d33 * signzz(i) +
532 . two*(d12 * signxy(i) + d13 * signzx(i) + d23 * signyz(i) )
533 a2 = a2 * a21
534C
535 IF(n==1)THEN
536 dsepp = et*(one+ ( epsp(i)*csd)**visp(i))
537 ELSE
538 dsepp = et * n * epsm(i) ** (n - 1)*
539 . (one+ ( epsp(i)*csd)**visp(i))
540 ENDIF
541 dcd = c1(i) * (d11**2 + d22**2 + d33**2) +
542 . two*c2(i)*(d11*d22 + d11*d33 + d22*d33) +
543 . two*g(i) *(d12**2 + d13**2 + d23**2)
544 lam1 = dcd - dcrm * dsepp * a2 -
545 . dcrf * ((one - f(i)) * (d11 + d22 + d33) + a1*a2)
546 IF(lam1/=zero) lamda(i) = max(zero,(vm * sigm2 - va))/lam1
547C
548C PLASTIC DEFORMATION
549C
550 dp11 = dp11 + lamda(i) * d11
551 dp22 = dp22 + lamda(i) * d22
552 dp33 = dp33 + lamda(i) * d33
553 dp12 = dp12 + lamda(i) * d12
554 dp13 = dp13 + lamda(i) * d13
555 dp23 = dp23 + lamda(i) * d23
556C NEW STRESS
557 signxx(i)= sigxx(i) - c1(i) * dp11 - c2(i) * (dp22 + dp33)
558 signyy(i)= sigyy(i) - c1(i) * dp22 - c2(i) * (dp11 + dp33)
559 signzz(i)= sigzz(i) - c1(i) * dp33 - c2(i) * (dp22 + dp11)
560 signxy(i)= sigxy(i) - two*g(i) * dp12
561 signzx(i)= sigzx(i) - two*g(i) * dp13
562 signyz(i)= sigyz(i) - two*g(i) * dp23
563C
564C....... EFFECTIVE PLASTIC DEFORMATION.
565C
566 epsp1(i) = signxx(i)*dp11 + signyy(i)*dp22 + signzz(i)*dp33+
567 . two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
568 epsp1(i) = epsp1(i)*a21
569 epsp1(i) = max(zero, epsp1(i))
570 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
571 sxx = signxx(i)-pn(i)
572 syy = signyy(i)-pn(i)
573 szz = signzz(i)-pn(i)
574 vm = half*(sxx**2 + syy**2 + szz**2)
575 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
576 vm = three*vm
577 var = three_half * q2 * pn(i)*sigm1
578 var = exp(var)
579 coh = half * (var + one/max(em20,var))
580 sih = half * (var - one/max(em20,var))
581 IF(pn(i)<=zero) THEN
582 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
583 ELSE
584 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
585 ENDIF
586C
587 yld = va
588 ENDDO ! ITER
589 fg(i) = fg(i) + (one-f(i))*(dp11 + dp22 +dp33)
590 sigy(i) = sigm(i)*sqrt(max(zero,va))
591 fn1(i) = fn1(i) + a1 * epsp1(i)
592 epsm(i) = epsm(i) + epsp1(i)
593 f(i) = fi + fg(i) + fn1(i)
594 IF(f(i)<zero) f(i)=zero
595 sigm(i)= (yeild0(i) + et * epsm(i)**n)*(one +
596 . (epsp(i)*csd)**visp(i))
597C....COMPUTE F*.....
598 IF(f(i)<=fc)THEN
599 fstar(i) = f(i)
600 ELSE
601 IF(ff==fc)THEN
602 fstar(i) = ep20
603 ELSE
604 fstar(i) = fc+ (fu-fc)*(f(i)-fc)/(ff-fc)
605 ENDIF
606 ENDIF
607C
608 uvar(i,1) = epsm(i)
609 uvar(i,2) = sigm(i)
610 dmg(i,1) = f(i)/ff
611 dmg(i,2) = fg(i)
612 dmg(i,3) = fn1(i)
613 dmg(i,4) = f(i)
614 dmg(i,5) = fstar(i)
615C
616 ENDIF
617 IF(off(i)==one.AND.fstar(i)>=ff)THEN
618 off(i)=four_over_5
619 WRITE(iout, 1000) ngl(i)
620 WRITE(istdo,1100) ngl(i),tt
621 ENDIF
622 ELSE
623 IF(off(i)<em01)off(i)=zero
624 off(i)=off(i)*four_over_5
625 ENDIF
626 ENDDO
627 ELSEIF(iflag==2) THEN
628C----------
629C IFLAG set to 2
630C..........
631 DO i=1,nel
632 IF(off(i)==one)THEN
633C
634 IF(f(i)<=fc)THEN
635 fstar(i) = f(i)
636 df = one
637 ELSE
638 df = (fu - fc)/(ff - fc)
639 fstar(i) = fc + df*(f(i)-fc)
640 ENDIF
641C
642 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
643 sxx = signxx(i)-pn(i)
644 syy = signyy(i)-pn(i)
645 szz = signzz(i)-pn(i)
646 vm = half*(sxx**2 + syy**2 + szz**2)
647 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
648 vm = three*vm
649C
650 sigm1 = one/max(em20,sigm(i))
651 sigm2 = sigm1**2
652 var = half * q2 * pn(i) * sigm1
653 var = exp(var)
654 coh = half*(var + one/max(em20,var))
655 sih = half*(var - one/max(em20,var))
656 IF(pn(i)<=zero) THEN
657 va = one +q3*fstar(i)**2 - two*q1*fstar(i)
658 ELSE
659 va= one+q3*fstar(i)**2 - two*q1*fstar(i)*coh
660 ENDIF
661 yld = sigm(i)
662 crit = vm * sigm2 - va
663 sigy(i) = yld*sqrt(max(zero, va))
664 IF(crit < zero .OR. sigy(i) == zero)THEN
665 lamda(i) = zero
666 ELSE
667 lamda(i) = zero
668 sigxx(i) = signxx(i)
669 sigyy(i) = signyy(i)
670 sigzz(i) = signzz(i)
671 sigxy(i) = signxy(i)
672 sigzx(i) = signzx(i)
673 sigyz(i) = signyz(i)
674C
675 IF(f(i)==one)THEN
676 a21 = ep20
677 ELSE
678 a21 = sigm1 /(one - f(i))
679 ENDIF
680 a1 = fn*exp(-half*((epsm(i)-epsn)/sn)**2)/sn
681 a1 = a1*sqr22
682C
683 dp11 = zero
684 dp22 = zero
685 dp33 = zero
686 dp12 = zero
687 dp13 = zero
688 dp23 = zero
689C
690 DO iter=1,5
691 IF(pn(i)<=zero)THEN
692 d11 = (two*signxx(i) - signyy(i) - signzz(i))*sigm2
693 d22 = (two*signyy(i) - signxx(i) - signzz(i))*sigm2
694 d33 = (two*signzz(i) - signxx(i) - signyy(i))*sigm2
695 d12 = six*signxy(i)*sigm2
696 d13 = six*signzx(i)*sigm2
697 d23 = six*signyz(i)*sigm2
698 dcrm = -two * vm * sigm2 * sigm1
699 dcrf = two * q1 * df - two * q3 * df * fstar(i)
700 ELSE
701 va1 = q1 * q2 * fstar(i) * sih*sigm1
702 d11 = (two*signxx(i) - signyy(i) - signzz(i))*sigm2 + va1
703 d22 = (two*signyy(i) - signxx(i) - signzz(i))*sigm2 + va1
704 d33 = (two*signzz(i) - signxx(i) - signyy(i))*sigm2 + va1
705 d12 = six * signxy(i)*sigm2
706 d13 = six * signzx(i)*sigm2
707 d23 = six * signyz(i)*sigm2
708 dcrm =-two * vm*sigm1*sigm2 - three * va1 * pn(i)*sigm1
709 dcrf = two * q1*df*coh - two*q3*df*fstar(i)
710 ENDIF
711C
712 a2 = d11 * signxx(i) + d22 * signyy(i) + d33 * signzz(i) +
713 . two*(d12 * signxy(i) + d13 * signzx(i) + d23 * signyz(i) )
714 a2 = a2 * a21
715C
716 IF(n==1)THEN
717 dsepp = et*(one + ( epsp(i)*csd)**visp(i))
718 ELSE
719 dsepp = et * n * epsm(i) ** (n - one)*
720 . (one + ( epsp(i)*csd)**visp(i))
721 ENDIF
722 dcd = c1(i) * (d11**2 + d22**2 + d33**2) +
723 . two*c2(i)*(d11*d22 + d11*d33 + d22*d33) +
724 . two*g(i) *(d12**2 + d13**2 + d23**2)
725 lam1 = dcd - dcrm * dsepp * a2 -
726 . dcrf * ((one - f(i)) * (d11 + d22 + d33) + a1*a2)
727 IF(lam1/=zero) lamda(i) = max(zero,(vm * sigm2 - va))/lam1
728C
729C PLASTIC DEFORMATION
730C
731 dp11 = dp11 + lamda(i) * d11
732 dp22 = dp22 + lamda(i) * d22
733 dp33 = dp33 + lamda(i) * d33
734 dp12 = dp12 + lamda(i) * d12
735 dp13 = dp13 + lamda(i) * d13
736 dp23 = dp23 + lamda(i) * d23
737C NEW STRESS
738 signxx(i)= sigxx(i) - c1(i) * dp11 - c2(i) * (dp22 + dp33)
739 signyy(i)= sigyy(i) - c1(i) * dp22 - c2(i) * (dp11 + dp33)
740 signzz(i)= sigzz(i) - c1(i) * dp33 - c2(i) * (dp22 + dp11)
741 signxy(i)= sigxy(i) - two*g(i) * dp12
742 signzx(i)= sigzx(i) - two*g(i) * dp13
743 signyz(i)= sigyz(i) - two*g(i) * dp23
744C
745C....... EFFECTIVE PLASTIC DEFORMATION.
746C
747 epsp1(i) = signxx(i)*dp11 + signyy(i)*dp22 + signzz(i)*dp33+
748 . two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
749 epsp1(i) = epsp1(i)*a21
750 epsp1(i) = max(zero, epsp1(i))
751
752 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
753 sxx = signxx(i)-pn(i)
754 syy = signyy(i)-pn(i)
755 szz = signzz(i)-pn(i)
756 vm = half*(sxx**2 + syy**2 + szz**2)
757 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
758 vm = three*vm
759 var = three_half * q2 * pn(i)*sigm1
760 var = exp(var)
761 coh = half * (var + one/max(em20,var))
762 sih = half * (var - one/max(em20,var))
763 IF(pn(i)<=zero) THEN
764 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
765 ELSE
766 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
767 ENDIF
768C
769 yld = va
770 ENDDO ! ITER
771 fg(i) = fg(i) + (one -f(i))*(dp11 + dp22 +dp33)
772 sigy(i) = sigm(i)*sqrt(max(zero,va))
773 fn1(i) = zero
774 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
775 IF(pn(i)>=zero) fn1(i)= fn1(i) + a1 * epsp1(i)
776 epsm(i) = epsm(i) + epsp1(i)
777 f(i) = fi + fg(i) + fn1(i)
778 IF(f(i)<zero) f(i)=zero
779 sigm(i)= (yeild0(i) + et * epsm(i)**n)*(one+
780 . (epsp(i)*csd)**visp(i))
781C....COMPUTE F*.....
782 IF(f(i)<=fc)THEN
783 fstar(i) = f(i)
784 ELSE
785 IF(ff==fc)THEN
786 fstar(i) = ep20
787 ELSE
788 fstar(i) = fc+ (fu-fc)*(f(i)-fc)/(ff-fc)
789 ENDIF
790 ENDIF
791C
792 uvar(i,1) = epsm(i)
793 uvar(i,2) = sigm(i)
794 dmg(i,1) = f(i)/ff
795 dmg(i,2) = fg(i)
796 dmg(i,3) = fn1(i)
797 dmg(i,4) = f(i)
798 dmg(i,5) = fstar(i)
799C
800 ENDIF
801CCsm41w11 +++
802 IF(off(i)==one.AND.fstar(i)>=ff)THEN
803 off(i)=four_over_5
804 WRITE(iout, 1000) ngl(i)
805 WRITE(istdo,1100) ngl(i),tt
806 ENDIF
807 ELSE
808
809 IF(off(i)<em01)off(i)=zero
810 off(i)=off(i)*four_over_5
811 ENDIF
812 ENDDO
813 ELSE
814C-----
815C IFLAG = 3
816C------------
817C-------------
818C
819C-----------
820 DO i=1,nel
821 IF(off(i)==one)THEN
822 IF(f(i)<=fc)THEN
823 fstar(i) = f(i)
824 df = one
825 ELSE
826 df = (fu- fc)/(ff-fc)
827 fstar(i)= fc + df*(f(i)-fc)
828 ENDIF
829C
830 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
831 sxx = signxx(i)-pn(i)
832 syy = signyy(i)-pn(i)
833 szz = signzz(i)-pn(i)
834 vm = half*(sxx**2 + syy**2 + szz**2)
835 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
836 vm = sqrt(three*vm)
837 sigm1 = one/max(sigm(i),em20)
838 var = three_half * q2 * pn(i)* sigm1
839 var = exp(var)
840 coh = half * (var + one/max(em20,var))
841 sih = half * (var - one/max(em20,var))
842 va = one + q3 * fstar(i)**2 - two *q1*fstar(i)*coh
843 va= sqrt(max(zero,va))
844 yld = sigm(i) * va
845C
846 sigy(i) = yld
847 crit = vm - yld
848 IF(vm < yld .OR. sigy(i) == zero)THEN
849 lamda(i)=zero
850 ELSE
851 lamda(i) = zero
852 sigxx(i) = signxx(i)
853 sigyy(i) = signyy(i)
854 sigzz(i) = signzz(i)
855 sigxy(i) = signxy(i)
856 sigzx(i) = signzx(i)
857 sigyz(i) = signyz(i)
858C
859 dp11 = zero
860 dp22 = zero
861 dp33 = zero
862 dp12 = zero
863 dp13 = zero
864 dp23 = zero
865C
866 IF(f(i)==one)THEN
867 a21 = ep20
868 ELSE
869 a21 = sigm1 /(one - f(i))
870 ENDIF
871 a1 = fn*exp(-half*((epsm(i)-epsn)/sn)**2)/sn
872 a1 = a1*sqr22
873c
874 DO iter=1,5
875 vm1 = one/max(vm,em20)
876 va1 = one/max(va,em20)
877 va11 =half*q1*q2*fstar(i)*sih*va1
878 d11 = half*(two*signxx(i)
879 . - signyy(i)-signzz(i))*vm1 + va11
880 d22 = half*(two*signyy(i)
881 . - signxx(i)-signzz(i))*vm1 + va11
882 d33 = half*(two*signzz(i)
883 . - signxx(i)-signyy(i))*vm1 + va11
884 d12 = three*signxy(i)*vm1
885 d13 = three*signzx(i)*vm1
886 d23 = three*signyz(i)*vm1
887C
888 a2 = d11*signxx(i) + d22*signyy(i) + d33*signzz(i) +
889 . two*(d12*signxy(i) + d13*signzx(i) + d23*signyz(i))
890 a2 = a2*a21
891C
892 dcrf=-sigm(i)*(q3*fstar(i)*df-q1*coh*df)*va1
893 dcrm= - va - three*va11*pn(i)*sigm1
894 IF(n==1)THEN
895 dsepp = et*(one + ( epsp(i)*csd)**visp(i))
896 ELSE
897 dsepp = et * n * epsm(i) ** (n - one)*
898 . (one + ( epsp(i)*csd)**visp(i))
899 ENDIF
900 dcd = c1(i)*(d11**2 + d22**2 + d33**2) +
901 . two*c2(i)*(d11*d22 + d11*d33 + d22*d33) +
902 . two*g(i)*(d12**2 + d13**2 + d23**2)
903
904 lam1= dcd - dcrm* dsepp*a2 -
905 . dcrf*((one-f(i))*(d11+d22+d33) + a1*a2)
906 IF(lam1/=zero) lamda(i) = max(zero,(vm - yld))/lam1
907C PLASTIC DEFORMATION
908 dp11 = dp11 + lamda(i) * d11
909 dp22 = dp22 + lamda(i) * d22
910 dp33 = dp33 + lamda(i) * d33
911 dp12 = dp12 + lamda(i) * d12
912 dp13 = dp13 + lamda(i) * d13
913 dp23 = dp23 + lamda(i) * d23
914C NEW STRESS
915 signxx(i)= sigxx(i) - c1(i)*dp11 - c2(i)*(dp22 + dp33)
916 signyy(i)= sigyy(i) - c1(i)*dp22 - c2(i)*(dp11 + dp33)
917 signzz(i)= sigzz(i) - c1(i)*dp33 - c2(i)*(dp22 + dp11)
918 signxy(i)= sigxy(i) - two*g(i)*dp12
919 signzx(i)= sigzx(i) - two*g(i)*dp13
920 signyz(i)= sigyz(i) - two*g(i)*dp23
921C....... AMISSIBLE PLASTIC DEFORMATION.
922 epsp1(i)= signxx(i)*dp11 + signyy(i)*dp22 + signzz(i)*dp33
923 . + two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
924 epsp1(i)=epsp1(i)*a21
925 epsp1(i) = max(zero, epsp1(i))
926 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
927 sxx = signxx(i)-pn(i)
928 syy = signyy(i)-pn(i)
929 szz = signzz(i)-pn(i)
930 vm = half*(sxx**2 + syy**2 + szz**2)
931 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
932 vm = sqrt(three*vm)
933 var = three_half * q2*pn(i) *sigm1
934 var = exp( var)
935 coh = half * ( var + one/max(em20,var))
936 sih = half * ( var - one/max(em20,var))
937 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
938 va = sqrt(max(zero,va))
939 yld = sigm(i) * va
940 ENDDO ! ITER
941C
942 fg(i) = fg(i) + (one-f(i))*(dp11+dp22+dp33)
943 pn(i) = (sigxx(i) + sigyy(i) + sigzz(i))*third
944 sigy(i) = yld
945 fn1(i) = zero
946 IF(pn(i)>=zero) fn1(i) = fn1(i) + a1*epsp1(i)
947 epsm(i) = epsm(i) + epsp1(i)
948 f(i) = fi + fg(i) + fn1(i)
949 IF(q1==zero.and.
950 . q2==zero.and.q3==zero) f(i)= zero
951 IF(f(i)<zero) f(i)=zero
952 sigm(i) = (yeild0(i) + et * epsm(i) ** n)*
953 . (one + (epsp(i)*csd) ** visp(i))
954C....COMPUTE F*.....
955 IF(f(i)<=fc)THEN
956 fstar(i) = f(i)
957 ELSE
958 IF(ff==fc)THEN
959 fstar(i) = ep20
960 ELSE
961 fstar(i) = fc+ (fu-fc)/(ff-fc)*(f(i)-fc)
962 ENDIF
963 ENDIF
964 uvar(i,1) = epsm(i)
965 uvar(i,2) = sigm(i)
966 dmg(i,1) = f(i)/ff
967 dmg(i,2) = fg(i)
968 dmg(i,3) = fn1(i)
969 dmg(i,4) = f(i)
970 dmg(i,5) = fstar(i)
971 ENDIF
972 IF(off(i)==one.AND.fstar(i)>=ff)THEN
973 off(i)=four_over_5
974 WRITE(iout, 1000) ngl(i)
975 WRITE(istdo,1100) ngl(i),tt
976 ENDIF
977 ELSE
978 IF(off(i)<em01)off(i)=zero
979 off(i)=off(i)*four_over_5
980 ENDIF
981 ENDDO
982 ENDIF
983 DO i=1,nel
984 defp(i)=uvar(i,1)
985 ENDDO
986C
987C tabulated law
988 IF(iyeild_tab > 0) THEN
989 iadbuf = ipm(7,mat(1))-1
990 xfac = uparam(iadbuf + 22 )
991 yfac = uparam(iadbuf + 23 )
992 DO i=1,nel
993 ndim= table(itab(1))%NDIM
994 IF(ndim==2)THEN
995 nxk=SIZE(table(itab(1))%X(2)%VALUES)
996 DO j=2,nxk
997 dx2 = table(itab(1))%X(2)%VALUES(j)*xfac- epsp(i)
998 IF(dx2 >= zero .OR. j == nxk)THEN
999 ipos=j-1
1000 r=(table(itab(1))%X(2)%VALUES(j)*xfac-epsp(i))/
1001 . (table(itab(1))%X(2)%VALUES(j)*xfac
1002 . -table(itab(1))%X(2)%VALUES(j-1)*xfac)
1003 EXIT
1004 ENDIF
1005 END DO ! NXK
1006 ELSE
1007 r = one
1008 ENDIF ! NDIM
1009C interpolation on the curve ---
1010 xx(1) = epsm(i)
1011 xx(2) = epsp(i)*xfac
1012 CALL table_interp_law76(table(itab(1)),ipos,xx,
1013 . r,dft,yy)
1014 uvar(i,2)= yfac*yy
1015 ENDDO ! NEL0
1016 ENDIF ! IYEILD_TAB
1017C ---------------------------------------------------------------------
1018 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
1019 1100 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10,
1020 . ' AT TIME :',g20.12)
1021 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
interface_table_mod::table_interp_law76
Definition table_mod.F:272
interface_table_mod::table_interp
Definition table_mod.F:242
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112