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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ sigeps56()

subroutine sigeps56 ( integer nel,
integer nuparam,
integer nuvar,
integer mfunc,
integer, dimension(mfunc) kfunc,
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,
integer, dimension(nel) ngl,
integer ipt,
integer, dimension(npropmi,*) ipm,
integer, dimension(nel) mat,
epsp,
integer ipla,
yld,
pla,
dpla,
amu )

Definition at line 33 of file sigeps56.F.

46CGW087 B PM ,MAT )
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51#include "comlock.inc"
52C-----------------------------------------------
53C G l o b a l P a r a m e t e r s
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56C---------+---------+---+---+--------------------------------------------
57C VAR | SIZE |TYP| RW| DEFINITION
58C---------+---------+---+---+--------------------------------------------
59C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
60C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
61C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
62C---------+---------+---+---+--------------------------------------------
63C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
64C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
65C NPF | * | I | R | FUNCTION ARRAY
66C TF | * | F | R | FUNCTION ARRAY
67C---------+---------+---+---+--------------------------------------------
68C TIME | 1 | F | R | CURRENT TIME
69C TIMESTEP| 1 | F | R | CURRENT TIME STEP
70C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
71C RHO0 | NEL | F | R | INITIAL DENSITY
72C RHO | NEL | F | R | DENSITY
73C VOLUME | NEL | F | R | VOLUME
74C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
75C EPSPXX | NEL | F | R | STRAIN RATE XX
76C EPSPYY | NEL | F | R | STRAIN RATE YY
77C ... | | | |
78C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
79C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
80C ... | | | |
81C EPSXX | NEL | F | R | STRAIN XX
82C EPSYY | NEL | F | R | STRAIN YY
83C ... | | | |
84C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
85C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
86C ... | | | |
87C---------+---------+---+---+--------------------------------------------
88C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
89C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
90C ... | | | |
91C SIGVXX | NEL | F | W | VISCOUS STRESS XX
92C SIGVYY | NEL | F | W | VISCOUS STRESS YY
93C ... | | | |
94C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
95C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
96C---------+---------+---+---+--------------------------------------------
97C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
98C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
99C---------+---------+---+---+--------------------------------------------
100#include "param_c.inc"
101#include "parit_c.inc"
102#include "scr17_c.inc"
103#include "com01_c.inc"
104#include "com08_c.inc"
105#include "units_c.inc"
106C
107 INTEGER NEL, NUPARAM, NUVAR,IPT,IPLA,
108 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*)
109 my_real
110 . time,timestep,uparam(*),
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 . epsp(nel),amu(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),dpla(nel)
130C-----------------------------------------------
131C I N P U T O U T P U T A r g u m e n t s
132C-----------------------------------------------
133 my_real
134 . uvar(nel,nuvar), off(nel), yld(nel),pla(nel)
135C-----------------------------------------------
136C VARIABLES FOR FUNCTION INTERPOLATION
137C-----------------------------------------------
138 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
139 my_real
140 . finter2, tf(*)
141 EXTERNAL finter2
142C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
143C Y : y = f(x)
144C X : x
145C DYDX : f'(x) = dy/dx
146C IFUNC(J): FUNCTION INDEX
147C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
148C NPF,TF : FUNCTION PARAMETER
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,IPOS1(MVSIZ),IPOS2(MVSIZ),IAD1(MVSIZ),
154 . ILEN1(MVSIZ),IAD2(MVSIZ),ILEN2(MVSIZ),
155 . IFUNC(100),NFUNCV,PFUN,
156 . IPOSP(MVSIZ),IADP(MVSIZ),ILENP(MVSIZ),NFUNCM,NRATEM,
157 . IPFLAG,IPARAM,NPAR,
158 . NINDX,INDX(MVSIZ),MX
159 my_real
160 . e,nu,p,dav,vm,r,fac,epst,ep1,ep2,ep3,ep4,ep5,ep6,
161 . e1,e2,e3,e4,e5,e6,c,cc,d,y,yp,e42,e52,e62,epst2,
162 . c1,g,g2,g3,
163 . epsmax,rate(mvsiz,2),yfac(mvsiz,2),
164 . y1(mvsiz),y2(mvsiz),h(mvsiz),dydx1(mvsiz),
165 . dydx2(mvsiz),fail(mvsiz),epsr1,
166 . epsr2,p0(mvsiz),pfac(mvsiz),
167 . dfdp(mvsiz)
168c . ONE,ZERO,EM20
169c DATA ZERO/0.0/,ONE/1.0/,EM20/1.E-20/
170C-----------------------------------------------
171C USER VARIABLES INITIALIZATION
172C-----------------------------------------------
173 IF (ivector/=1) THEN
174 mx = mat(1)
175 nfunc = ipm(10,mx)
176 DO j=1,nfunc
177 ifunc(j)=ipm(10+j,mx)
178 ENDDO
179 ELSE
180 mx = mat(1)
181 nfuncm = 0
182 nfuncv = ipm(10,mx)
183 nfuncm = max(nfuncm,nfuncv)
184
185 DO j=1,nfuncm
186 IF(nfuncv>=j) THEN
187 ifunc(j)=ipm(10+j,mx)
188 ENDIF
189 ENDDO
190 ENDIF
191C
192 nratem = 0
193C
194 mx = mat(1)
195 iadbufv = ipm(7,mx)-1
196 e = uparam(iadbufv+2)
197 g = uparam(iadbufv+5)
198 g2 = two*g
199 g3 = three*g
200 nu = uparam(iadbufv+6)
201 c1 = e/three/(one - two*nu)
202C
203 nrate = nint(uparam(iadbufv+1))
204 nratem = max(nratem,nrate)
205 epsmax=uparam(iadbufv+6+2*nrate+1)
206 IF(epsmax==zero)THEN
207 IF(tf(npf(ifunc(1)+1)-1)==zero)THEN
208 epsmax=tf(npf(ifunc(1)+1)-2)
209 ELSE
210 epsmax= ep30
211 ENDIF
212 ENDIF
213 epsr1 =uparam(iadbufv+6+2*nrate+2)
214 IF(epsr1==zero)epsr1=ep30
215 epsr2 =uparam(iadbufv+6+2*nrate+3)
216 IF(epsr2==zero)epsr2=twoep30
217 ipflag = 0
218 mx = mat(1)
219 npar = ipm(9,mx)
220 iparam = 15 + 2*nrate
221 DO i=1,nel
222 pfac(i) = one
223 pfun = 0
224 IF (npar>=iparam) THEN
225 j = nrate+1
226 ifunc(j)=ipm(10+j,mx)
227 pfun=nint(uparam(iadbufv+iparam))
228 IF (pfun>0) ipflag = ipflag + 1
229 ENDIF
230 ENDDO
231C
232 IF (isigi==0) THEN
233 IF(time==zero)THEN
234 DO i=1,nel
235 uvar(i,1)=zero
236 uvar(i,2)=zero
237 DO j=1,nrate
238 uvar(i,j+2)=zero
239 ENDDO
240 IF (pfun>0) uvar(i,nrate+5)=zero
241 ENDDO
242 ENDIF
243 ENDIF
244C-----------------------------------------------
245C
246 DO i=1,nel
247C
248 pla(i) = uvar(i,1)
249 p0(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third
250 dav = (depsxx(i)+depsyy(i)+depszz(i))*third
251 signxx(i)=sigoxx(i)+p0(i)+g2*(depsxx(i)-dav)
252 signyy(i)=sigoyy(i)+p0(i)+g2*(depsyy(i)-dav)
253 signzz(i)=sigozz(i)+p0(i)+g2*(depszz(i)-dav)
254 signxy(i)=sigoxy(i)+g *depsxy(i)
255 signyz(i)=sigoyz(i)+g *depsyz(i)
256 signzx(i)=sigozx(i)+g *depszx(i)
257C
258 soundsp(i) = sqrt((c1+four*g/three)/rho0(i))
259 viscmax(i) = zero
260 dpla(i) = zero
261 ENDDO
262C-------------------
263C STRAIN RATE
264C-------------------
265 DO i=1,nel
266C 087 DAV = (EPSPXX(I)+EPSPYY(I)+EPSPZZ(I))/3.
267C 087 E1 = EPSPXX(I) - DAV
268C 087 E2 = EPSPYY(I) - DAV
269C 087 E3 = EPSPZZ(I) - DAV
270C 087 E4 = 0.5*EPSPXY(I)
271C 087 E5 = 0.5*EPSPYZ(I)
272C 087 E6 = 0.5*EPSPZX(I)
273C 087 EPSP(I) =.5*(E1**2+E2**2+E3**2) +E4**2+E5**2+E6**2
274C 087 EPSP(I) = SQRT(3.*EPSP(I))/1.5
275C-------------------
276C STRAIN principal 1, 4 newton iterations
277C-------------------
278 dav = (epsxx(i)+epsyy(i)+epszz(i)) * third
279 e1 = epsxx(i) - dav
280 e2 = epsyy(i) - dav
281 e3 = epszz(i) - dav
282 e4 = half*epsxy(i)
283 e5 = half*epsyz(i)
284 e6 = half*epszx(i)
285C -y = (e1-x)(e2-x)(e3-x)
286C - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
287C + 2e4 e5 e6
288C e1 + e2 + e3 = 0 => terme en x^2 = 0
289C y = x^3 + c x + d
290c yp= 3 x^2 + c
291 e42 = e4*e4
292 e52 = e5*e5
293 e62 = e6*e6
294 c = - e1*e1 - e2*e2 - e3*e3 - e42 - e52 - e62
295 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
296 & - two*e4*e5*e6
297 cc = c*third
298 epst = sqrt(-cc)
299 epst2 = epst * epst
300 y = (epst2 + c)* epst + d
301 IF(abs(y)>em8)THEN
302 epst = onep75 * epst
303 epst2 = epst * epst
304 y = (epst2 + c)* epst + d
305 yp = three*epst2 + c
306 IF(yp/=zero)epst = epst - y/yp
307 epst2 = epst * epst
308 y = (epst2 + c)* epst + d
309 yp = three*epst2 + c
310 IF(yp/=zero)epst = epst - y/yp
311 epst2 = epst * epst
312 y = (epst2 + c)* epst + d
313 yp = three*epst2 + c
314 IF(yp/=zero)epst = epst - y/yp
315 epst2 = epst * epst
316 y = (epst2 + c)* epst + d
317 yp = three*epst2 + c
318 IF(yp/=zero)epst = epst - y/yp
319 epst = epst + dav
320 ENDIF
321C-------------------
322C tension failure
323C-------------------
324 fail(i) = max(zero,min(one,
325 . (epsr2-epst)/(epsr2-epsr1) ))
326C
327c jbm023
328c
329c FAIL(I)=1.
330 ENDDO
331C-------------------
332C CRITERE
333C-------------------
334 DO i=1,nel
335 jj(i) = 1
336 ENDDO
337 IF (ivector/=1) THEN
338 DO i=1,nel
339 DO j=2,nrate-1
340 IF(epsp(i)>=uparam(iadbufv+6+j)) jj(i) = j
341 ENDDO
342 ENDDO
343 ELSE
344 DO j=2,nratem-1
345 DO i=1,nel
346 IF(nrate-1>=j) THEN
347 IF(epsp(i)>=uparam(iadbufv+6+j)) jj(i) = j
348 ENDIF
349 ENDDO
350 ENDDO
351 ENDIF
352C
353 DO i=1,nel
354 rate(i,1)=uparam(iadbufv+6+jj(i))
355 rate(i,2)=uparam(iadbufv+6+jj(i)+1)
356 yfac(i,1)=uparam(iadbufv+6+nrate+jj(i))
357 yfac(i,2)=uparam(iadbufv+6+nrate+jj(i)+1)
358 ENDDO
359C
360 DO i=1,nel
361 j1 = jj(i)
362 j2 = j1+1
363C 037 IPOS1(I) = NINT(UVAR(I,J1))
364 ipos1(i) = nint(uvar(i,j1+2))
365 iad1(i) = npf(ifunc(j1)) / 2 + 1
366 ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i) - ipos1(i)
367C 037 IPOS2(I) = NINT(UVAR(I,J2))
368 ipos2(i) = nint(uvar(i,j2+2))
369 iad2(i) = npf(ifunc(j2)) / 2 + 1
370 ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i) - ipos2(i)
371 ENDDO
372C
373 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
374 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
375C
376C--- PRESSURE DEPENDENT YIELD FUNCTION FACTOR
377 IF (ipflag==nel.AND.(iparit==0)) THEN
378C------- to optimize in SPMD & PARIT/ON case add a global flag
379C in a case of homogeneous element groups (to do if needed)
380 DO i=1,nel
381 iposp(i) = nint(uvar(i,nrate+5))
382 iadp(i) = npf(ifunc(pfun)) / 2 + 1
383 ilenp(i) = npf(ifunc(pfun)+1) / 2 - iadp(i) - iposp(i)
384 ENDDO
385 CALL vinter2(tf,iadp,iposp,ilenp,nel,p0,dfdp,pfac)
386 DO i=1,nel
387 uvar(i,nrate+5) = iposp(i)
388 ENDDO
389 ELSEIF (ipflag>0) THEN
390 DO i=1,nel
391 IF (pfun>0) THEN
392 iposp(i) = nint(uvar(i,nrate+5))
393 iadp(i) = npf(ifunc(pfun)) / 2 + 1
394 ilenp(i) = npf(ifunc(pfun)+1) / 2
395 . - iadp(i) - iposp(i)
396 pfac(i) = finter2(tf ,iadp(i),iposp(i),ilenp(i),
397 . p0(i),dfdp(i))
398 uvar(i,nrate+5) = iposp(i)
399 ENDIF
400 ENDDO
401 ENDIF
402 DO i=1,nel
403 j1 = jj(i)
404 j2 = j1+1
405 y1(i)=y1(i)*yfac(i,1)
406 y2(i)=y2(i)*yfac(i,2)
407 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
408 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
409 dydx1(i)=dydx1(i)*yfac(i,1)
410 dydx2(i)=dydx2(i)*yfac(i,2)
411 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
412 h(i) = h(i) * max(zero, pfac(i))
413 yld(i) = yld(i) * max(zero, pfac(i))
414
415C 037 UVAR(I,J1) = IPOS1(I)
416C 037 UVAR(I,J2) = IPOS2(I)
417 uvar(i,j1+2) = ipos1(i)
418 uvar(i,j2+2) = ipos2(i)
419 ENDDO
420C-------------------
421C PROJECTION
422C-------------------
423 IF(ipla==0)THEN
424 DO i=1,nel
425 vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
426 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
427 vm =sqrt(three*vm)
428 r = min(one,yld(i)/ max(vm,em20))
429c P = C1 * (RHO(I)/RHO0(I)- ONE)
430 p = c1 * amu(i)
431 signxx(i)=signxx(i)*r-p
432 signyy(i)=signyy(i)*r-p
433 signzz(i)=signzz(i)*r-p
434 signxy(i)=signxy(i)*r
435 signyz(i)=signyz(i)*r
436 signzx(i)=signzx(i)*r
437 pla(i)=pla(i) + (one - r)*vm/max(g3+h(i),em20)
438 uvar(i,1)=pla(i)
439 dpla(i) = (one - r)*vm/max(g3+h(i),em20)
440 ENDDO
441 ELSEIF(ipla==2)THEN
442 DO i=1,nel
443 vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
444 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
445 vm =sqrt(three*vm)
446 r = min(one,yld(i)/ max(vm,em20))
447c P = C1 * (RHO(I)/RHO0(I)-ONE)
448 p = c1 * amu(i)
449 signxx(i)=signxx(i)*r-p
450 signyy(i)=signyy(i)*r-p
451 signzz(i)=signzz(i)*r-p
452 signxy(i)=signxy(i)*r
453 signyz(i)=signyz(i)*r
454 signzx(i)=signzx(i)*r
455 pla(i)=pla(i) + (one-r)*vm/max(g3,em20)
456 uvar(i,1)=pla(i)
457 dpla(i) = (one-r)*vm/max(g3,em20)
458 ENDDO
459 ELSEIF(ipla==1)THEN
460 DO i=1,nel
461 vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
462 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
463 vm =sqrt(three*vm)
464 r = min(one,yld(i)/ max(vm,em20))
465C plastic strain increment.
466 dpla(i)=(one - r)*vm/max(g3+h(i),em20)
467C actual yield stress.
468 yld(i)=max(yld(i)+dpla(i)*h(i),zero)
469 r = min(one,yld(i)/ max(vm,em20))
470C
471c P = C1 * (RHO(I)/RHO0(I)-ONE)
472 p = c1 * amu(i)
473 signxx(i)=signxx(i)*r-p
474 signyy(i)=signyy(i)*r-p
475 signzz(i)=signzz(i)*r-p
476 signxy(i)=signxy(i)*r
477 signyz(i)=signyz(i)*r
478 signzx(i)=signzx(i)*r
479 pla(i)=pla(i) + dpla(i)
480 uvar(i,1)=pla(i)
481 ENDDO
482 ENDIF
483C-------------------
484 DO i=1,nel
485 IF(off(i)<em01) off(i)=zero
486 IF(off(i)<one) off(i)=off(i)*four_over_5
487 ENDDO
488C
489 nindx=0
490 DO i=1,nel
491 IF(pla(i)>epsmax.AND.off(i)==one) THEN
492 off(i)=four_over_5
493 nindx=nindx+1
494 indx(nindx)=i
495 ENDIF
496 ENDDO
497 IF(nindx>0)THEN
498 DO j=1,nindx
499#include "lockon.inc"
500 WRITE(iout, 1000) ngl(indx(j))
501 WRITE(istdo,1100) ngl(indx(j)),tt
502#include "lockoff.inc"
503 ENDDO
504 ENDIF
505C
506 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
507 1100 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10,
508 . ' AT TIME :',g11.4)
509C
510 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72