OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps48c.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps48c (nel0, nuparam, nuvar, mfunc, kfunc, npf, npt0, ipt, iflag, tf, time, timestep, uparam, rho0, area, eint, thkly, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thk, pla, uvar, off, ngl, ipm, mat, etse, gs, yld, epsp, dpla_i, israte, asrate, epsd_pg, inloc, dplanl, loff)

Function/Subroutine Documentation

◆ sigeps48c()

subroutine sigeps48c ( integer nel0,
integer nuparam,
integer nuvar,
integer mfunc,
integer, dimension(mfunc) kfunc,
integer, dimension(*) npf,
integer npt0,
integer ipt,
integer, dimension(*) iflag,
tf,
time,
timestep,
uparam,
rho0,
area,
eint,
thkly,
epspxx,
epspyy,
epspxy,
epspyz,
epspzx,
depsxx,
depsyy,
depsxy,
depsyz,
depszx,
epsxx,
epsyy,
epsxy,
epsyz,
epszx,
sigoxx,
sigoyy,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signxy,
signyz,
signzx,
sigvxx,
sigvyy,
sigvxy,
sigvyz,
sigvzx,
soundsp,
viscmax,
thk,
pla,
uvar,
off,
integer, dimension(nel0) ngl,
integer, dimension(npropmi,*) ipm,
integer, dimension(nel0) mat,
etse,
gs,
yld,
intent(inout) epsp,
dpla_i,
integer israte,
intent(in) asrate,
intent(in) epsd_pg,
integer inloc,
dplanl,
intent(in) loff )

Definition at line 28 of file sigeps48c.F.

43C-----------------------------------------------
44C I m p l i c i t T y p e s
45C-----------------------------------------------
46#include "implicit_f.inc"
47C-----------------------------------------------
48C G l o b a l P a r a m e t e r s
49C-----------------------------------------------
50#include "mvsiz_p.inc"
51C---------+---------+---+---+--------------------------------------------
52C VAR | SIZE |TYP| RW| DEFINITION
53C---------+---------+---+---+--------------------------------------------
54C NEL0 | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
55C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
56C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
57C---------+---------+---+---+--------------------------------------------
58C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
59C IFUNC | NFUNC | I | R | FUNCTION INDEX
60C NPF | * | I | R | FUNCTION ARRAY
61C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
62C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
63C IFLAG | * | I | R | GEOMETRICAL FLAGS
64C TF | * | F | R | FUNCTION ARRAY
65C---------+---------+---+---+--------------------------------------------
66C TIME | 1 | F | R | CURRENT TIME
67C TIMESTEP| 1 | F | R | CURRENT TIME STEP
68C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
69C RHO0 | NEL0 | F | R | INITIAL DENSITY
70C AREA | NEL0 | F | R | AREA
71C EINT | 2*NEL0 | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
72C THKLY | NEL0 | F | R | LAYER THICKNESS
73C EPSPXX | NEL0 | F | R | STRAIN RATE XX
74C EPSPYY | NEL0 | F | R | STRAIN RATE YY
75C ... | | | |
76C DEPSXX | NEL0 | F | R | STRAIN INCREMENT XX
77C DEPSYY | NEL0 | F | R | STRAIN INCREMENT YY
78C ... | | | |
79C EPSXX | NEL0 | F | R | STRAIN XX
80C EPSYY | NEL0 | F | R | STRAIN YY
81C ... | | | |
82C SIGOXX | NEL0 | F | R | OLD ELASTO PLASTIC STRESS XX
83C SIGOYY | NEL0 | F | R | OLD ELASTO PLASTIC STRESS YY
84C ... | | | |
85C---------+---------+---+---+--------------------------------------------
86C SIGNXX | NEL0 | F | W | NEW ELASTO PLASTIC STRESS XX
87C SIGNYY | NEL0 | F | W | NEW ELASTO PLASTIC STRESS YY
88C ... | | | |
89C SIGVXX | NEL0 | F | W | VISCOUS STRESS XX
90C SIGVYY | NEL0 | F | W | VISCOUS STRESS YY
91C ... | | | |
92C SOUNDSP | NEL0 | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
93C VISCMAX | NEL0 | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
94C---------+---------+---+---+--------------------------------------------
95C THK | NEL0 | F |R/W| THICKNESS
96C PLA | NEL0 | F |R/W| PLASTIC STRAIN
97C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
98C OFF | NEL0 | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
99C---------+---------+---+---+--------------------------------------------
100C C o m m o n B l o c k s
101C-----------------------------------------------
102#include "param_c.inc"
103C-----------------------------------------------
104C I N P U T A r g u m e n t s
105C-----------------------------------------------
106C
107 INTEGER NEL0, NUPARAM, NUVAR, NPT0,ISRATE, IPT,IFLAG(*),
108 . NGL(NEL0),MAT(NEL0),IPM(NPROPMI,*),INLOC
109 my_real ,INTENT(IN) :: asrate
110 my_real ,DIMENSION(NEL0) ,INTENT(IN) :: epsd_pg
111 my_real ,DIMENSION(NEL0) ,INTENT(INOUT) :: epsp
112 my_real
113 . time,timestep,uparam(*),
114 . area(nel0),rho0(nel0),eint(nel0,2),
115 . thkly(nel0),pla(nel0),
116 . epspxx(nel0),epspyy(nel0),
117 . epspxy(nel0),epspyz(nel0),epspzx(nel0),
118 . depsxx(nel0),depsyy(nel0),
119 . depsxy(nel0),depsyz(nel0),depszx(nel0),
120 . epsxx(nel0) ,epsyy(nel0) ,
121 . epsxy(nel0) ,epsyz(nel0) ,epszx(nel0) ,
122 . sigoxx(nel0),sigoyy(nel0),
123 . sigoxy(nel0),sigoyz(nel0),sigozx(nel0),
124 . gs(*),dplanl(nel0)
125 my_real, DIMENSION(NEL0), INTENT(IN) :: loff
126C-----------------------------------------------
127C O U T P U T A r g u m e n t s
128C-----------------------------------------------
129 my_real
130 . signxx(nel0),signyy(nel0),
131 . signxy(nel0),signyz(nel0),signzx(nel0),
132 . sigvxx(nel0),sigvyy(nel0),
133 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
134 . soundsp(nel0),viscmax(nel0),etse(nel0),
135 . dpla_i(nel0)
136C-----------------------------------------------
137C I N P U T O U T P U T A r g u m e n t s
138C-----------------------------------------------
139 my_real
140 . uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
141C-------------------------
142C variables non utilisees (Fonctions utilisateur)
143C-------------------------
144 INTEGER NPF(*),MFUNC,KFUNC(MFUNC)
145 my_real
146 . tf(*)
147C-----------------------------------------------
148C L o c a l V a r i a b l e s
149C-----------------------------------------------
150 INTEGER I,J,N,NINDX,NMAX,IADBUF,INDEX(MVSIZ)
151 my_real
152 . svm(mvsiz),dr(mvsiz),
153 . aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
154 . pp(mvsiz),qq(mvsiz),fail(mvsiz),h(mvsiz),hs(mvsiz),
155 . epsm(mvsiz),
156 . ylo(mvsiz),epsgm(mvsiz),
157 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz)
158 my_real
159 . r,umr,nux,a,b,c,s11,s22,s12,p,p2,dezz,sigz,s1,s2,s3,
160 . vm2,epst,nnu2,f,df,q2,yld_i,sigpxx,e1,a11,a21,g1,g31,
161 . nnu11,nu11,nu21,nu31,nu41,nu51,nu61,sigpyy,sigpxy,sigm1,
162 . epsm1,epsr11,epsr21,yd1,hc1,cn1,c11,c21,cm1,c31,cl1,
163 . fisokin1,epsgm1,eps01,pa,pb,pc,pda,pdb,yy,
164 . dsxx,dsyy,dsxy,dexx,deyy,dexy, alpha
165C
166 DATA nmax/3/
167C=======================================================================
168 IF(time==zero)THEN
169 DO i=1,nel0
170 uvar(i,1)=zero
171 uvar(i,2)=zero
172 uvar(i,3)=zero
173 ENDDO
174 ENDIF
175 iadbuf = ipm(7,mat(1))-1
176 e1 = uparam(iadbuf+1)
177 nux = uparam(iadbuf+2)
178 sigm1= uparam(iadbuf+4)
179 epsm1= uparam(iadbuf+5)
180 epsr11=uparam(iadbuf+6)
181 epsr21=uparam(iadbuf+7)
182 yd1 = uparam(iadbuf+3)
183 hc1 = uparam(iadbuf+8)
184 cn1 = uparam(iadbuf+9)
185 c11 = uparam(iadbuf+10)
186 c21 = uparam(iadbuf+11)
187 cm1 = uparam(iadbuf+12)
188 fisokin1=uparam(iadbuf+15)
189 g1 = uparam(iadbuf+16)
190 g31 = uparam(iadbuf+18)
191 a11 = uparam(iadbuf+20)
192 a21 = uparam(iadbuf+21)
193c EPSGM1=UPARAM(IADBUF+22)
194 eps01= uparam(iadbuf+23)
195 c31 = uparam(iadbuf+24)
196 cl1 = uparam(iadbuf+25)
197C
198 nnu11 = nux / (one - nux)
199 nnu2 = nnu11*nnu11
200 nu11 = one/(one-nux)
201 nu21 = one/(one+nux)
202 nu31 = one-nnu11
203 nu41 = one + nnu2 + nnu11
204 nu51 = one + nnu2 - two*nnu11
205 nu61 = half- nnu2 + half*nnu11
206
207C
208C------------------------------------------
209C KINEMATIC HARDENING
210C------------------------------------------
211 DO i=1,nel0
212! sigoxx(i) = sigoxx(i) - uvar(i,1)
213! SIGOYY(I) = SIGOYY(I) - UVAR(I,2)
214! SIGOXY(I) = SIGOXY(I) - UVAR(I,3)
215 ENDDO
216C
217C--- ELASTIC STRESS ESTIMATE
218C
219 DO i=1,nel0
220C
221 signxx(i)=sigoxx(i) - uvar(i,1) +a11*depsxx(i)+a21*depsyy(i)
222 signyy(i)=sigoyy(i) - uvar(i,2) +a21*depsxx(i)+a11*depsyy(i)
223 signxy(i)=sigoxy(i) - uvar(i,3) +g1 *depsxy(i)
224 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
225 signzx(i)=sigozx(i)+gs(i) *depszx(i)
226 sigexx(i) = signxx(i)
227 sigeyy(i) = signyy(i)
228 sigexy(i) = signxy(i)
229 soundsp(i) = sqrt(a11/rho0(i))
230 viscmax(i) = zero
231 etse(i) = one
232C-------------------
233C STRAIN RATE
234C-------------------
235 IF (israte == 0) THEN
236 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
237 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
238 . + epspxy(i)*epspxy(i) ) )
239 ELSE
240 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
241 END IF
242C-------------------
243C STRAIN
244C-------------------
245 epst = half*( epsxx(i)+epsyy(i)
246 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
247 . + epsxy(i)*epsxy(i) ) )
248 fail(i) = max(zero,min(one,(epsr21-epst)/(epsr21-epsr11)))
249C----
250 ENDDO
251C-------------------
252C CURRENT YIELD AND HARDENING
253C-------------------
254C
255 DO i=1,nel0
256 IF(pla(i)<=zero) THEN
257 pa=yd1
258 ELSE
259 pa=yd1+hc1*pla(i)**cn1
260 ENDIF
261 IF(epsp(i)<=eps01) THEN
262 pb=0.
263 ELSEIF(pla(i)<=zero) THEN
264 pb=c11*log(epsp(i)/eps01)
265 ELSE
266 pb=(c11-c21*pla(i)**cm1)*log(epsp(i)/eps01)
267 ENDIF
268 IF(epsp(i)<=zero) THEN
269 pc=zero
270 ELSE
271 pc=c31*epsp(i)**cl1
272 ENDIF
273C-----
274 IF(pla(i)>zero. and .cn1>=one) THEN
275 pda = hc1*cn1*pla(i)**(cn1-one)
276 ELSEIF(pla(i)>zero. and .cn1<one)THEN
277 pda = hc1*cn1*pla(i)**(one -cn1)
278 ELSE
279 pda = e1
280 ENDIF
281 IF(pla(i)<=0. or .epsp(i)<=eps01) THEN
282 pdb = zero
283 ELSEIF(cm1>=one) THEN
284 pdb = c21*cm1*pla(i)**(cm1 - one)*log(epsp(i)/eps01)
285 ELSE
286 pdb = c21*cm1*pla(i)**(one-cm1)*log(epsp(i)/eps01)
287 ENDIF
288C
289c IF (ICC1/=1.AND.CN1/=ZERO)
290c & EPSGM1=((SIGM1/RQ-YD1)/HC1)**(1./CN1)
291c IF (PLA(I)>=EPSGM1) THEN
292c YLD(I) = SIGM1
293c HS(I) = ZERO
294c ENDIF
295C----
296 ylo(i)= yd1 + pc
297 yy = pa + pb + pc
298 hs(i) = pda + pdb
299 yld(i)= min(sigm1+pc, yy)
300 IF (yld(i)<yy) hs(i) = zero
301 hs(i) = fail(i)*hs(i)
302 yld(i) = fail(i)*yld(i)
303 ENDDO
304C
305C--- KINEMATIC HARDENING
306C
307 IF (fisokin1==one) THEN
308 DO i=1,nel0
309 yld(i) = fail(i)*ylo(i)
310 ENDDO
311 ELSEIF (fisokin1>zero) THEN
312 DO i=1,nel0
313 yld(i) = (one-fisokin1)*yld(i) + fisokin1*fail(i)*ylo(i)
314 ENDDO
315 ENDIF
316C
317C-------------------------
318C PROJECTION IFLAG = 0 (projection radiale)
319C-------------------------
320C
321 IF(iflag(1)==0) THEN
322 DO i=1,nel0
323 svm(i)=sqrt(signxx(i)*signxx(i)
324 . +signyy(i)*signyy(i)
325 . -signxx(i)*signyy(i)
326 . +three*signxy(i)*signxy(i))
327 r = min(one,yld(i)/max(em20,svm(i)))
328 signxx(i)=signxx(i)*r
329 signyy(i)=signyy(i)*r
330 signxy(i)=signxy(i)*r
331 umr = one - r
332 dpla_i(i) = off(i)*svm(i)*umr/(g31+hs(i))
333 pla(i) = pla(i) + dpla_i(i)
334 s1=half*(signxx(i)+signyy(i))
335 IF (inloc == 0) THEN
336 dezz = dpla_i(i) * s1 / yld(i)
337 dezz = -(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
338 thk(i) = thk(i) + dezz*thkly(i)*off(i)
339 ENDIF
340 ENDDO
341C
342C-------------------------
343C PROJECTION IFLAG = 1
344C-------------------------
345C
346 ELSEIF(iflag(1)==1)THEN
347C
348C--- CRITERE DE VON MISES
349C
350 DO i=1,nel0
351 s1=signxx(i)+signyy(i)
352 s2=signxx(i)-signyy(i)
353 s3=signxy(i)
354 aa(i)=fourth*s1*s1
355 bb(i)=three_over_4*s2*s2+3.*s3*s3
356 svm(i)=sqrt(aa(i)+bb(i))
357 IF (inloc == 0) THEN
358 dezz = -(depsxx(i)+depsyy(i))*nnu11
359 thk(i) = thk(i) + dezz*thkly(i)*off(i)
360 ENDIF
361 ENDDO
362C
363C--- GATHER PLASTIC FLOW
364C
365 nindx=0
366 DO i=1,nel0
367 IF(svm(i)>yld(i).AND.off(i)==one) THEN
368 nindx=nindx+1
369 index(nindx)=i
370 ENDIF
371 ENDDO
372C
373C--- DEP EN CONTRAINTE PLANE
374C
375 IF(nindx/=0) THEN
376 DO j=1,nindx
377 i=index(j)
378 hs(i) = max(zero,hs(i))
379 etse(i)= hs(i)/(hs(i)+e1)
380 h(i) = (one-fisokin1)*hs(i)
381 dpla_j(i)=(svm(i)-yld(i))/(g31+hs(i))
382 ENDDO
383C
384 DO n=1,nmax
385#include "vectorize.inc"
386 DO j=1,nindx
387 i =index(j)
388 dpla_i(i)=dpla_j(i)
389 yld_i = yld(i)+h(i)*dpla_i(i)
390 dr(i) = half*e1*dpla_i(i)/yld_i
391 pp(i) = one/(one+dr(i)*nu11)
392 qq(i) = one/(one + three*dr(i)*nu21)
393 p2 = pp(i)*pp(i)
394 q2 = qq(i)*qq(i)
395 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
396 df =-(aa(i)*nu11*p2*pp(i)+three*bb(i)*nu21*q2*qq(i))
397 . *(e1-two*dr(i)*h(i))/yld_i
398 . -two*h(i)*yld_i
399 IF(dpla_i(i)>zero) THEN
400 dpla_j(i) = max(zero,dpla_i(i)-f/df)
401 ELSE
402 dpla_j(i) = zero
403 ENDIF
404 ENDDO
405 ENDDO
406C
407C--- CONTRAINTES PLASTIQUEMENT ADMISSIBLES
408C
409#include "vectorize.inc"
410 DO j=1,nindx
411 i=index(j)
412 pla(i) = pla(i) + dpla_i(i)
413 s1=(signxx(i)+signyy(i))*pp(i)
414 s2=(signxx(i)-signyy(i))*qq(i)
415 signxx(i)=half*(s1+s2)
416 signyy(i)=half*(s1-s2)
417 signxy(i)=signxy(i)*qq(i)
418 IF (inloc == 0) THEN
419 dezz =-nu31*dr(i)*s1/e1
420 thk(i) = thk(i) + dezz*thkly(i)*off(i)
421 ENDIF
422 ENDDO
423 ENDIF
424C
425C-------------------------
426C PROJECTION IFLAG = 2
427C-------------------------
428C projection radial sur le deviateur sur un critere reduit
429C projection elastique en z => sig33 = 0
430C le coef. de reduction du critere est tel que
431C l'on se trouve sur le critere apres les 2 projections
432C
433 ELSEIF(iflag(1)==2)THEN
434C
435 DO i=1,nel0
436 p = -(signxx(i)+signyy(i))*third
437 s11 = signxx(i)+p
438 s22 = signyy(i)+p
439C s33 = p = -(S11 + S22)
440 s12 = signxy(i)
441C
442 p2 = p*p
443 vm2= three*(s12*s12 - s11*s22)
444 a = p2*nu41 + vm2
445 vm2= 3.*p2 + vm2
446 b = p2*nu61
447 c = p2*nu51 - yld(i)*yld(i)
448 r = min(one,(-b + sqrt(max(zero,b*b-a*c)))/max(a ,em20))
449 signxx(i) = s11*r - p
450 signyy(i) = s22*r - p
451 signxy(i) = s12*r
452C signzz = p*r - p
453C proj. signzz = 0.
454 umr = one-r
455 sigz = nnu11*p*umr
456 signxx(i) = signxx(i) + sigz
457 signyy(i) = signyy(i) + sigz
458 svm(i)=sqrt(vm2)
459 dpla_i(i) = off(i)*svm(i)*umr/(g31+hs(i))
460 pla(i) = pla(i) + dpla_i(i)
461 IF (inloc == 0) THEN
462 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
463 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
464 thk(i) = thk(i) + dezz*thkly(i)*off(i)
465 ENDIF
466 IF(r<one) etse(i)= hs(i)/(hs(i)+e1)
467 ENDDO
468 ENDIF
469C
470 DO i=1,nel0
471 IF(pla(i)>epsm1.AND.off(i)==one)off(i)=four_over_5
472 ENDDO
473C
474C--- KINEMATIC HARDENING
475C
476 IF (fisokin1/=zero) THEN
477 DO i=1,nel0
478 dsxx = sigexx(i) - signxx(i)
479 dsyy = sigeyy(i) - signyy(i)
480 dsxy = sigexy(i) - signxy(i)
481 dexx = (dsxx - nux*dsyy)
482 deyy = (dsyy - nux*dsxx)
483 dexy = two*(one+nux)*dsxy
484 alpha = fisokin1*hs(i)/(e1+hs(i))/three
485 sigpxx = alpha*(four*dexx+two*deyy)
486 sigpyy = alpha*(four*deyy+two*dexx)
487 sigpxy = alpha*dexy
488C
489 signxx(i) = signxx(i) + uvar(i,1)
490 signyy(i) = signyy(i) + uvar(i,2)
491 signxy(i) = signxy(i) + uvar(i,3)
492 uvar(i,1) = uvar(i,1) + sigpxx
493 uvar(i,2) = uvar(i,2) + sigpyy
494 uvar(i,3) = uvar(i,3) + sigpxy
495 ENDDO
496 ENDIF
497C
498C--- NON-LOCAL THICKNESS VARIATION
499C
500 IF (inloc > 0) THEN
501 DO i = 1,nel0
502 IF (loff(i) == one) THEN
503 svm(i) = sqrt(signxx(i)*signxx(i)
504 . + signyy(i)*signyy(i)
505 . - signxx(i)*signyy(i)
506 . + three*signxy(i)*signxy(i))
507 dezz = max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(svm(i),em20)
508 dezz = -nux*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e1) - dezz
509 thk(i) = thk(i) + dezz*thkly(i)*off(i)
510 ENDIF
511 ENDDO
512 ENDIF
513C
514 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21