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

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps44c (nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, thkly, off, etse, epsd_pg, 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, soundsp, viscmax, thk, pla, uvar, gs, yld, epsp, dpla_i, asrate, nvartmp, vartmp, sigp, inloc, dplanl, loff)

Function/Subroutine Documentation

◆ sigeps44c()

subroutine sigeps44c ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
integer, intent(in) mfunc,
integer, dimension(mfunc), intent(in) kfunc,
integer, dimension(*), intent(in) npf,
dimension(*), intent(in) tf,
intent(in) time,
intent(in) timestep,
intent(in) uparam,
intent(in) rho0,
intent(in) thkly,
intent(inout) off,
intent(out) etse,
intent(in) epsd_pg,
intent(in) epspxx,
intent(in) epspyy,
intent(in) epspxy,
intent(in) epspyz,
intent(in) epspzx,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) epsxx,
intent(in) epsyy,
intent(in) epsxy,
intent(in) epsyz,
intent(in) epszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) soundsp,
intent(out) viscmax,
intent(inout) thk,
intent(inout) pla,
intent(inout) uvar,
intent(in) gs,
intent(inout) yld,
intent(inout) epsp,
intent(out) dpla_i,
intent(in) asrate,
integer, intent(in) nvartmp,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
intent(inout) sigp,
integer, intent(in) inloc,
intent(in) dplanl,
intent(in) loff )

Definition at line 30 of file sigeps44c.F.

43C-----------------------------------------------
44C I m p l i c i t T y p e s
45C-----------------------------------------------
46#include "implicit_f.inc"
47C---------+---------+---+---+--------------------------------------------
48C VAR | SIZE |TYP| RW| DEFINITION
49C---------+---------+---+---+--------------------------------------------
50C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
51C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
52C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
53C---------+---------+---+---+--------------------------------------------
54C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
55C IFUNC | NFUNC | I | R | FUNCTION INDEX
56C NPF | * | I | R | FUNCTION ARRAY
57C TF | * | F | R | FUNCTION ARRAY
58C---------+---------+---+---+--------------------------------------------
59C TIME | 1 | F | R | CURRENT TIME
60C TIMESTEP| 1 | F | R | CURRENT TIME STEP
61C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
62C RHO0 | NEL | F | R | INITIAL DENSITY
63C THKLY | NEL | F | R | LAYER THICKNESS
64C EPSPXX | NEL | F | R | STRAIN RATE XX
65C EPSPYY | NEL | F | R | STRAIN RATE YY
66C ... | | | |
67C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
68C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
69C ... | | | |
70C EPSXX | NEL | F | R | STRAIN XX
71C EPSYY | NEL | F | R | STRAIN YY
72C ... | | | |
73C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
74C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
75C ... | | | |
76C---------+---------+---+---+--------------------------------------------
77C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
78C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
79C ... | | | |
80C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
81C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
82C---------+---------+---+---+--------------------------------------------
83C THK | NEL | F |R/W| THICKNESS
84C PLA | NEL | F |R/W| PLASTIC STRAIN
85C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
86C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
87C---------+---------+---+---+--------------------------------------------
88C C o m m o n B l o c k s
89C-----------------------------------------------
90#include "param_c.inc"
91#include "mvsiz_p.inc"
92C-----------------------------------------------
93C I N P U T A r g u m e n t s
94C-----------------------------------------------
95 INTEGER ,INTENT(IN) :: NEL,NUPARAM, NUVAR,NVARTMP,INLOC
96 my_real ,INTENT(IN) :: time,timestep,asrate
97 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) :: uparam
98 my_real ,DIMENSION(NEL) ,INTENT(IN) :: epsd_pg
99 my_real ,DIMENSION(NEL) ,INTENT(IN) :: rho0,thkly,gs,dplanl,
100 . epspxx,epspyy,epspxy,epspyz,epspzx,
101 . depsxx,depsyy,depsxy,depsyz,depszx,
102 . epsxx ,epsyy ,epsxy ,epsyz ,epszx ,
103 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx
104 INTEGER ,INTENT(IN) :: NPF(*),MFUNC,KFUNC(MFUNC)
105 my_real ,INTENT(IN) :: tf(*)
106 my_real ,DIMENSION(NEL), INTENT(IN) :: loff
107 my_real ,DIMENSION(NEL), INTENT(INOUT) :: epsp
108C-----------------------------------------------
109C O U T P U T A r g u m e n t s
110C-----------------------------------------------
111 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: soundsp,viscmax,etse,
112 . signxx,signyy,signxy,signyz,signzx,dpla_i
113C-----------------------------------------------
114C I N P U T O U T P U T A r g u m e n t s
115C-----------------------------------------------
116 INTEGER ,DIMENSION(NEL,NVARTMP) ,INTENT(INOUT) :: VARTMP
117 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: pla,off,thk,yld
118 my_real ,DIMENSION(NEL,3) ,INTENT(INOUT) :: sigp
119 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
120C-----------------------------------------------
121C L o c a l V a r i a b l e s
122C-----------------------------------------------
123 INTEGER :: I,J,N,NINDX,NMAX,IADBUF,ICC1,ISRATE,VFLAG
124 my_real :: r,umr,nux,a,b,c,s11,s22,s12,p,p2,dezz,
125 . sigz,s1,s2,s3,vm2,epst,f,df,q2,yld_i,e1,a11,a21,
126 . g1,g31,nnu11,nu11,nu21,nu31,deve1,deve2,deve3,deve4,dpdt,
127 . epsm1,epsr11,epsr21,fisokin1,ca1,cb1,cn1,cc1,cp1,
128 . dsxx,dsyy,dsxy,dexx,deyy,dexy,alpha,yscale,dav
129 INTEGER ,DIMENSION(MVSIZ) :: INDEX,IPOS,ILEN,IAD,IPOS0
130 my_real ,DIMENSION(MVSIZ) :: svm,dr,aa,bb,dpla_j,pp,qq,fail,h,hs,
131 . ylo,zeror,sigexx,sigeyy,sigexy,sigm,epsgm,dfdpla
132 my_real, DIMENSION(MVSIZ) :: rq
133C
134 DATA nmax/3/
135C-----------------------------------------------
136C MATERIAL PARAMETERS
137C-----------------------------------------------
138 e1 = uparam(1)
139 nux = uparam(2)
140 epsm1 = uparam(5)
141 epsr11 = uparam(6)
142 epsr21 = uparam(7)
143 ca1 = uparam(3)
144 cb1 = uparam(8)
145 cn1 = uparam(9)
146 icc1 = nint(uparam(10))
147 cc1 = uparam(11)
148 cp1 = uparam(12)
149 fisokin1 = uparam(15)
150 g1 = uparam(16)
151 g31 = uparam(18)
152 a11 = uparam(20)
153 a21 = uparam(21)
154 israte = nint(uparam(13))
155 vflag = nint(uparam(23))
156 yscale = uparam(24)
157C
158 nnu11 = nux / (one - nux)
159 nu11 = one/(one-nux)
160 nu21 = one/(one+nux)
161 nu31 = one - nnu11
162C
163 DO i=1,nel
164 epsgm(i) = uparam(22)
165 sigm(i) = uparam(4)
166 ENDDO
167C
168 IF (mfunc > 0) THEN
169 zeror(1:nel) = zero
170 ipos0(1:nel) = 1
171 iad(1:nel) = npf(kfunc(1)) / 2 + 1
172 ilen(1:nel) = npf(kfunc(1)+1) / 2 - iad(1:nel) - ipos0(1:nel)
173 CALL vinter(tf,iad,ipos0,ilen,nel,zeror,dfdpla,ylo)
174 ylo(1:nel) = yscale * ylo(1:nel)
175 ENDIF
176C-----------------------------------------------
177C ELASTIC STRESS ESTIMATE
178C-----------------------------------------------
179 IF (fisokin1 > zero) THEN
180 signxx(1:nel) = sigoxx(1:nel) - sigp(1:nel,1) + a11*depsxx(1:nel) + a21*depsyy(1:nel)
181 signyy(1:nel) = sigoyy(1:nel) - sigp(1:nel,2) + a21*depsxx(1:nel) + a11*depsyy(1:nel)
182 signxy(1:nel) = sigoxy(1:nel) - sigp(1:nel,3) + g1 *depsxy(1:nel)
183 ELSE
184 signxx(1:nel) = sigoxx(1:nel) + a11*depsxx(1:nel) + a21*depsyy(1:nel)
185 signyy(1:nel) = sigoyy(1:nel) + a21*depsxx(1:nel) + a11*depsyy(1:nel)
186 signxy(1:nel) = sigoxy(1:nel) + g1 *depsxy(1:nel)
187 ENDIF
188 signyz(1:nel) = sigoyz(1:nel) + gs(1:nel) * depsyz(1:nel)
189 signzx(1:nel) = sigozx(1:nel) + gs(1:nel) * depszx(1:nel)
190 sigexx(1:nel) = signxx(1:nel)
191 sigeyy(1:nel) = signyy(1:nel)
192 sigexy(1:nel) = signxy(1:nel)
193C
194 soundsp(1:nel) = sqrt(a11/rho0(1:nel))
195 viscmax(1:nel) = zero
196 etse(1:nel) = one
197C-----------------------------------------------
198C TOTAL OR DEVIATORIC STRAIN-RATE COMPUTATION
199C-----------------------------------------------
200 IF (vflag == 1) THEN
201 epsp(1:nel) = uvar(1:nel,1)
202 ENDIF
203 IF (israte > 0) THEN
204 IF (vflag == 3) THEN
205 DO i=1,nel
206 dav = (epspxx(i)+epspyy(i))*third
207 deve1 = epspxx(i) - dav
208 deve2 = epspyy(i) - dav
209 deve3 = - dav
210 deve4 = half*epspxy(i)
211 epsp(i) = half*(deve1**2 + deve2**2 + deve3**2) + deve4**2
212 epsp(i) = sqrt(three*epsp(i))/three_half
213 epsp(i) = asrate*epsp(i) + (one - asrate)*uvar(i,1)
214 uvar(i,1) = epsp(i)
215 ENDDO
216 ELSEIF (vflag == 2) THEN
217 DO i=1,nel
218 epsp(i) = asrate*epsd_pg(i) + (one - asrate)*uvar(i,1)
219 uvar(i,1) = epsp(i)
220 ENDDO
221 ENDIF
222 ELSEIF (israte == 0) THEN
223 IF (vflag == 3) THEN
224 DO i=1,nel
225 dav = (epspxx(i)+epspyy(i))*third
226 deve1 = epspxx(i) - dav
227 deve2 = epspyy(i) - dav
228 deve3 = - dav
229 deve4 = half*epspxy(i)
230 epsp(i) = half*(deve1**2 + deve2**2 + deve3**2) + deve4**2
231 epsp(i) = sqrt(three*epsp(i))/three_half
232 uvar(i,1) = epsp(i)
233 ENDDO
234 ELSEIF (vflag == 2) THEN
235 DO i=1,nel
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 ENDDO
240 ENDIF
241 ENDIF
242C-------------------
243C STRAIN & TENSION FAILURE
244C-------------------
245C
246 DO i=1,nel
247 epst = half*( epsxx(i)+epsyy(i)
248 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
249 . + epsxy(i)*epsxy(i) ) )
250 fail(i) = max(zero,min(one,(epsr21-epst)/(epsr21-epsr11)))
251 dpla_i(i) = zero
252 ENDDO
253C-------------------
254C CURRENT YIELD AND HARDENING
255C-------------------
256C
257 IF (mfunc > 0) THEN
258 ipos(1:nel) = vartmp(1:nel,1)
259 iad(1:nel) = npf(kfunc(1)) / 2 + 1
260 ilen(1:nel) = npf(kfunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
261 CALL vinter(tf,iad,ipos,ilen,nel,pla,dfdpla,yld)
262 vartmp(1:nel,1) = ipos(1:nel)
263 ENDIF
264C
265 rq(1:nel) = one
266 IF(cc1 /= zero) rq(1:nel) = one + (cc1*epsp(1:nel))**cp1
267C--------------------
268 IF((mfunc > 0) .AND. (ca1 == zero)) THEN
269 DO i=1,nel
270 ylo(i) = ylo(i) * rq(i)
271 IF (pla(i)>zero) THEN
272 yld(i) = yscale * yld(i) * rq(i)
273 hs(i) = yscale * dfdpla(i) * rq(i)
274 ELSE
275 yld(i) = yscale * yld(i) * rq(i)
276 hs(i) = e1
277 ENDIF
278 ENDDO
279C---------------------------
280 ELSEIF ((mfunc > 0) .AND. (ca1 /= zero)) THEN
281C--------------------
282 DO i=1,nel
283 ylo(i) = ylo(i) + ca1 * (rq(i)-one)
284 IF (pla(i)>zero) THEN
285 IF (cn1 == one) THEN
286 yld(i) = yscale * yld(i) + (ca1 + cb1*pla(i)) * (rq(i)-one)
287 hs(i) = yscale * dfdpla(i) + cb1 * (rq(i)-one)
288 ELSE
289 yld(i) = yscale * yld(i) + (ca1 + cb1*pla(i)**cn1) * (rq(i)-one)
290 IF (cn1>one) THEN
291 hs(i) = yscale * dfdpla(i) + cn1*cb1*(rq(i)-one) * (pla(i)**(cn1-one))
292 ELSE
293 hs(i) = yscale * dfdpla(i) + cn1*cb1*(rq(i)-one) / ((pla(i)**(one-cn1)))
294 ENDIF
295 ENDIF
296 ELSE ! PLA <= 0
297 yld(i) = yscale * yld(i) + ca1 * (rq(i)-one)
298 hs(i) = e1
299 ENDIF
300 ENDDO
301 ELSE
302C--------------------
303 DO i=1,nel
304 ylo(i) = ca1 * rq(i)
305 IF (pla(i)>zero) THEN
306 IF (cn1 == one) THEN
307 yld(i) = (ca1 + cb1*pla(i)) * rq(i)
308 hs(i) = cb1 * rq(i)
309 ELSE
310 yld(i) = (ca1 + cb1*pla(i)**cn1) * rq(i)
311 IF (cn1>one) THEN
312 hs(i) = cn1*cb1*rq(i) * (pla(i)**(cn1-one))
313 ELSE
314 hs(i) = cn1*cb1*rq(i) / ((pla(i)**(one-cn1)))
315 ENDIF
316 ENDIF
317 ELSE
318 yld(i) = ca1 * rq(i)
319 hs(i) = e1
320 ENDIF
321 ENDDO
322 ENDIF
323
324 DO i=1,nel
325 IF(icc1 == 1) sigm(i) = sigm(i) * rq(i)
326 IF (icc1 /= 1 .and. cn1 /= zero .and. cb1 /= zero)
327 & epsgm(i)=((sigm(i)/rq(i)-ca1)/cb1)**(one/cn1)
328 IF (pla(i)>=epsgm(i)) THEN
329 yld(i) = sigm(i)
330 hs(i) = zero
331 ENDIF
332 hs(i) = fail(i)*hs(i)
333C------ kinematic hardening
334 yld(i)= fail(i)*((one-fisokin1)*yld(i)+fisokin1*ylo(i))
335 yld(i) = max(yld(i),em20)
336 ENDDO
337
338C----------------------------------------
339C--- VON MISES
340C
341 DO i=1,nel
342 s1=signxx(i)+signyy(i)
343 s2=signxx(i)-signyy(i)
344 s3=signxy(i)
345 aa(i)=fourth*s1*s1
346 bb(i)=three_over_4*s2*s2+3.*s3*s3
347 svm(i)=sqrt(aa(i)+bb(i))
348 IF (inloc == 0) THEN
349 dezz = -(depsxx(i)+depsyy(i))*nnu11
350 thk(i) = thk(i) + dezz*thkly(i)*off(i)
351 ENDIF
352 ENDDO
353C
354C--- GATHER PLASTIC FLOW
355C
356 nindx = 0
357 DO i = 1,nel
358 IF ((svm(i) > yld(i)).AND.(off(i) == one)) THEN
359 nindx=nindx+1
360 index(nindx)=i
361 ENDIF
362 ENDDO
363C
364C--- DEP EN CONTRAINTE PLANE
365C
366 IF (nindx > 0) THEN
367#include "vectorize.inc"
368 DO j=1,nindx
369 i=index(j)
370 hs(i) = max(zero,hs(i))
371 dpla_j(i)=(svm(i)-yld(i))/(g31+hs(i))
372 etse(i)= hs(i)/(hs(i)+e1)
373 h(i) = hs(i)*(one-fisokin1)
374 ENDDO
375C
376 DO n=1,nmax
377#include "vectorize.inc"
378 DO j=1,nindx
379 i=index(j)
380 dpla_i(i)=dpla_j(i)
381 yld_i =yld(i)+h(i)*dpla_i(i)
382 dr(i) =half*e1*dpla_i(i)/yld_i
383 pp(i) =one/(one + dr(i)*nu11)
384 qq(i) =one/(one + three*dr(i)*nu21)
385 p2 =pp(i)*pp(i)
386 q2 =qq(i)*qq(i)
387 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
388 df=-(aa(i)*nu11*p2*pp(i)+ three*bb(i)*nu21*q2*qq(i))
389 . *(e1-two*dr(i)*h(i))/yld_i
390 . -two*h(i)*yld_i
391 df = sign(max(abs(df),em20),df)
392 IF(dpla_i(i)>zero) THEN
393 dpla_j(i)=max(zero,dpla_i(i)-f/df)
394 ELSE
395 dpla_j(i)=zero
396 ENDIF
397 ENDDO
398 ENDDO
399C
400C--- CONTRAINTES PLASTIQUEMENT ADMISSIBLES
401C
402#include "vectorize.inc"
403 DO j=1,nindx
404 i=index(j)
405 pla(i) = pla(i) + dpla_i(i)
406 s1=(signxx(i)+signyy(i))*pp(i)
407 s2=(signxx(i)-signyy(i))*qq(i)
408 signxx(i)=half*(s1+s2)
409 signyy(i)=half*(s1-s2)
410 signxy(i)=signxy(i)*qq(i)
411 IF (inloc == 0) THEN
412 dezz = - nu31*dr(i)*s1/e1
413 thk(i) = thk(i) + dezz*thkly(i)*off(i)
414 ENDIF
415 ENDDO
416
417 ENDIF ! nindx /= 0
418C
419 DO i=1,nel
420 IF ((pla(i) > epsm1).AND.(off(i) == one)) off(i) = four_over_5
421 IF (vflag == 1) THEN
422 dpdt = dpla_i(i)/max(em20,timestep)
423 uvar(i,1) = asrate * dpdt + (one - asrate) * uvar(i,1)
424 epsp(i) = uvar(i,1)
425 ENDIF
426 ENDDO
427C
428C--- KINEMATIC HARDENING
429C
430 IF (fisokin1 > zero) THEN
431 DO i=1,nel
432 dsxx = sigexx(i) - signxx(i)
433 dsyy = sigeyy(i) - signyy(i)
434 dsxy = sigexy(i) - signxy(i)
435 dexx = (dsxx - nux*dsyy)
436 deyy = (dsyy - nux*dsxx)
437 dexy = two*(one+nux)*dsxy
438 alpha = fisokin1*hs(i)/(e1+hs(i)) * third
439 signxx(i) = signxx(i) + sigp(i,1)
440 signyy(i) = signyy(i) + sigp(i,2)
441 signxy(i) = signxy(i) + sigp(i,3)
442 sigp(i,1) = sigp(i,1) + alpha*(four*dexx+two*deyy)
443 sigp(i,2) = sigp(i,2) + alpha*(four*deyy+two*dexx)
444 sigp(i,3) = sigp(i,3) + alpha*dexy
445 ENDDO
446 ENDIF
447C
448C--- NON-LOCAL THICKNESS VARIATION
449C
450 IF (inloc > 0) THEN
451 DO i = 1,nel
452 IF (loff(i) == one) THEN
453 svm(i) = sqrt(signxx(i)*signxx(i)
454 . + signyy(i)*signyy(i)
455 . - signxx(i)*signyy(i)
456 . + three*signxy(i)*signxy(i))
457 dezz = max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(svm(i),em20)
458 dezz = -nux*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e1) - dezz
459 thk(i) = thk(i) + dezz*thkly(i)*off(i)
460 ENDIF
461 ENDDO
462 ENDIF
463C-----------
464 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72
zeror
subroutine zeror(a, n)
Definition zero.F:39