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

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps93c (nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, soundsp, thk, pla, uvar, rho0, off, etse, thkly, shf, yld, hardm, seq, epsp, asrate, nvartmp, vartmp, dpla, inloc, dplanl, loff)

Function/Subroutine Documentation

◆ sigeps93c()

subroutine sigeps93c ( integer nel,
integer nuparam,
integer nuvar,
integer mfunc,
integer, dimension(mfunc) kfunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
epspxx,
epspyy,
epspxy,
epspyz,
epspzx,
depsxx,
depsyy,
depsxy,
depsyz,
depszx,
sigoxx,
sigoyy,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signxy,
signyz,
signzx,
soundsp,
thk,
pla,
uvar,
rho0,
off,
etse,
thkly,
shf,
yld,
hardm,
seq,
epsp,
asrate,
integer nvartmp,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
dpla,
integer inloc,
dplanl,
intent(in) loff )

Definition at line 30 of file sigeps93c.F.

41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C G l o b a l P a r a m e t e r s
47C-----------------------------------------------
48#include "mvsiz_p.inc"
49C-----------------------------------------------
50C I N P U T A r g u m e n t s
51C-----------------------------------------------
52C
53 INTEGER NEL, NUPARAM, NUVAR,NVARTMP,INLOC
54 my_real
55 . time,timestep(nel),uparam(nuparam),
56 . rho0(nel),thkly(nel),pla(nel),
57 . epspxx(nel),epspyy(nel),
58 . epspxy(nel),epspyz(nel),epspzx(nel),
59 . depsxx(nel),depsyy(nel),
60 . depsxy(nel),depsyz(nel),depszx(nel),
61 . sigoxx(nel),sigoyy(nel),
62 . sigoxy(nel),sigoyz(nel),sigozx(nel),
63 . shf(*),hardm(nel),seq(nel),dplanl(nel),
64 . epsp(nel),asrate
65 my_real, DIMENSION(NEL), INTENT(IN) :: loff
66C-----------------------------------------------
67C O U T P U T A r g u m e n t s
68C-----------------------------------------------
69 INTEGER, INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
70 my_real
71 . signxx(nel),signyy(nel),
72 . signxy(nel),signyz(nel),signzx(nel),
73 . soundsp(nel),etse(nel),dpla(nel)
74C-----------------------------------------------
75C I N P U T O U T P U T A r g u m e n t s
76C-----------------------------------------------
77 my_real
78 . uvar(nel,nuvar),off(nel),thk(nel),yld(nel)
79C-----------------------------------------------
80C VARIABLES FOR FUNCTION INTERPOLATION
81C-----------------------------------------------
82 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
83 my_real
84 . tf(*)
85C-----------------------------------------------
86C L o c a l V a r i a b l e s
87C-----------------------------------------------
88 INTEGER I,II,J,NITER,ITER,NINDX,INDEX(MVSIZ),
89 . J1,J2,ITAB,JJ(MVSIZ),VP,
90 . IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
91 . IAD2(MVSIZ),IPOS2(MVSIZ),ILEN2(MVSIZ)
92 my_real
93 . e11,e22,a11,a22,a12,g12,g13,g23,d13,d23,d33,invd33,
94 . ff,gg,hh,nn,sigy,qr1,cr1,qr2,cr2,nu12,nu13,nu23,
95 . fac,sighl(mvsiz),h(mvsiz),dydx1(mvsiz),dydx2(mvsiz),
96 . y1(mvsiz),y2(mvsiz),dezz(mvsiz),deelzz(mvsiz),
97 . dpxx(mvsiz),dpyy(mvsiz),dpzz(mvsiz),dpxy(mvsiz),
98 . normxx,normyy,normxy,dpla_dlam(mvsiz),dlam,
99 . dphi_dlam(mvsiz),phi(mvsiz),dav,deve1,deve2,
100 . deve3,deve4,dfdsig2,sig_dfdsig,ddep,dpdt
101 my_real,
102 . DIMENSION(:), ALLOCATABLE :: rate, yfac
103C-----------------------------------------------
104C USER VARIABLES INITIALIZATION
105C-----------------------------------------------
106C
107 ! Elastic parameters
108 a11 = uparam(1) ! 1st Diag. component of plane stress elastic matrix
109 a22 = uparam(2) ! 2nd Diag. component of plane stress elastic matrix
110 a12 = uparam(3) ! Non Diag. component of plane stress elastic matrix
111 d13 = uparam(6) ! Component 13 of the compliance matrix
112 d23 = uparam(8) ! Component 23 of the compliance matrix
113 d33 = uparam(9) ! Component 33 of the compliance matrix
114 g12 = uparam(10) ! Shear modulus in 12
115 g13 = uparam(11) ! Shear modulus in 13
116 g23 = uparam(12) ! Shear modulus in 23
117 e11 = uparam(13) ! Young modulus in 11
118 e22 = uparam(14) ! Young modulus in 22
119 nu12 = uparam(16) ! Poisson's ration in 12
120 nu13 = uparam(17) ! Poisson's ration in 13
121 nu23 = uparam(18) ! Poisson's ration in 23
122 ! Hill yield criterion coefficient
123 ff = uparam(19) ! First Hill coefficient
124 gg = uparam(20) ! Second Hill coefficient
125 hh = uparam(21) ! Third Hill coefficient
126 nn = uparam(24) ! Fourth Hill coefficient
127 ! Continuous hardening yield stress
128 sigy = uparam(25) ! Initial yield stress
129 qr1 = uparam(26) ! Voce 1 first parameter
130 cr1 = uparam(27) ! Voce 1 second parameter
131 qr2 = uparam(28) ! Voce 2 first parameter
132 cr2 = uparam(29) ! Voce 2 second parameter
133 ! Strain-rate computation flag
134 vp = nint(uparam(30))
135 ! Tabulated hardening yield stress
136 itab = 0
137 IF (mfunc > 0) THEN
138 itab = 1
139 ALLOCATE(rate(mfunc),yfac(mfunc))
140 DO i = 1,mfunc
141 rate(i) = uparam(30+i)
142 yfac(i) = uparam(30+mfunc+i)
143 ENDDO
144 ENDIF
145C
146 ! Recovering internal variables
147 DO i=1,nel
148 ! Checking deletion flag value
149 IF (off(i) < one) off(i) = four_over_5*off(i)
150 IF (off(i) < em01) off(i) = zero
151 ! Hourglass coefficient
152 etse(i) = one
153 h(i) = zero
154 ! Plastic strain increment
155 dpla(i) = zero
156 dpxx(i) = zero
157 dpyy(i) = zero
158 dpzz(i) = zero
159 dpxy(i) = zero
160 ! Computing strain-rate (only if more than 1 strain-rate function is indicated in the input deck)
161 IF ((itab == 1).AND.(mfunc > 1)) THEN
162 ! Plastic strain-rate
163 IF (vp == 1) THEN
164 epsp(i) = uvar(i,1)
165 ! Deviatoric strain-rate
166 ELSEIF (vp == 3) THEN
167 dav = (epspxx(i)+epspyy(i))*third
168 deve1 = epspxx(i) - dav
169 deve2 = epspyy(i) - dav
170 deve3 = - dav
171 deve4 = half*epspxy(i)
172 epsp(i) = half*(deve1**2 + deve2**2 + deve3**2) + deve4**2
173 epsp(i) = sqrt(three*epsp(i))/three_half
174 epsp(i) = asrate*epsp(i) + (one - asrate)*uvar(i,1)
175 uvar(i,1) = epsp(i)
176 ENDIF
177 ! Looking for corresponding rate in the tables
178 jj(i) = 1
179 DO j = 2,mfunc-1
180 IF (epsp(i) >= rate(j)) jj(i) = j
181 ENDDO
182 ENDIF
183 ENDDO
184C
185 ! Computing yield stress
186 ! -> Continuous yield stress
187 IF (itab == 0) THEN
188 DO i = 1,nel
189 yld(i) = sigy
190 . + qr1*(one - exp(-cr1*pla(i)))
191 . + qr2*(one - exp(-cr2*pla(i)))
192 h(i) = qr1*cr1*exp(-cr1*pla(i)) + qr2*cr2*exp(-cr2*pla(i))
193 ENDDO
194 ELSE
195 ! -> Tabulated yield stress
196 ! -> No strain-rate effect
197 IF (mfunc == 1) THEN
198 ! Recovering position tables
199 ipos1(1:nel) = vartmp(1:nel,1)
200 iad1(1:nel) = npf(kfunc(1)) / 2 + 1
201 ilen1(1:nel) = npf(kfunc(1)+1) / 2 - iad1(1:nel) - ipos1(1:nel)
202 ! Computing the plastic strain interpolation
203 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
204 ! Storing the position
205 vartmp(1:nel,1) = ipos1(1:nel)
206 ! Computing the yield stress and its derivative
207 yld(1:nel) = yfac(1)*y1(1:nel)
208 h(1:nel) = yfac(1)*dydx1(1:nel)
209 ! -> Strain-rate effect
210 ELSE
211 DO i=1,nel
212 j1 = jj(i)
213 j2 = j1 + 1
214 ! Recovering position tables for the first rate
215 ipos1(i) = vartmp(i,j1)
216 iad1(i) = npf(kfunc(j1)) / 2 + 1
217 ilen1(i) = npf(kfunc(j1)+1) / 2 - iad1(i) - ipos1(i)
218 ! Recovering position tables for the second rate
219 ipos2(i) = vartmp(i,j2)
220 iad2(i) = npf(kfunc(j2)) / 2 + 1
221 ilen2(i) = npf(kfunc(j2)+1) / 2 - iad2(i) - ipos2(i)
222 ENDDO
223 ! Computing the plastic strain interpolation
224 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
225 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
226 ! Storing new positions and computing yield stress
227 DO i=1,nel
228 j1 = jj(i)
229 j2 = j1+1
230 ! Hardening yield stress
231 y1(i) = y1(i)*yfac(j1)
232 y2(i) = y2(i)*yfac(j2)
233 fac = (epsp(i) - rate(j1))/(rate(j2) - rate(j1))
234 yld(i) = y1(i) + fac*(y2(i)-y1(i))
235 yld(i) = max(yld(i),em20)
236 ! Derivatives of yield stress
237 dydx1(i) = dydx1(i)*yfac(j1)
238 dydx2(i) = dydx2(i)*yfac(j2)
239 h(i) = dydx1(i)+fac*(dydx2(i)-dydx1(i))
240 ! New positions
241 vartmp(i,j1) = ipos1(i)
242 vartmp(i,j2) = ipos2(i)
243 ENDDO
244 ENDIF
245 ENDIF
246c
247 !========================================================================
248 ! - COMPUTATION OF TRIAL VALUES
249 !========================================================================
250 DO i=1,nel
251c
252 ! Computation of the trial stress tensor
253 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
254 signyy(i) = sigoyy(i) + a12*depsxx(i) + a22*depsyy(i)
255 signxy(i) = sigoxy(i) + g12*depsxy(i)
256 signyz(i) = sigoyz(i) + shf(i)*g23*depsyz(i)
257 signzx(i) = sigozx(i) + shf(i)*g13*depszx(i)
258C
259 ! Hill equivalent stress
260 sighl(i) = (ff + hh)*signyy(i)**2 + (gg + hh)*signxx(i)**2
261 . - two*hh*signxx(i)*signyy(i) + two*nn*signxy(i)**2
262 sighl(i) = sqrt(max(zero,sighl(i)))
263C
264 ENDDO
265c
266 !========================================================================
267 ! - COMPUTATION OF YIELD FONCTION
268 !========================================================================
269 phi(1:nel) = sighl(1:nel) - yld(1:nel)
270 ! Checking plastic behavior for all elements
271 nindx = 0
272 index = 0
273 DO i=1,nel
274 IF ((phi(i)>zero).AND.(off(i) == one)) THEN
275 nindx = nindx+1
276 index(nindx) = i
277 ENDIF
278 ENDDO
279c
280 !====================================================================
281 ! - PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM (NEWTON ITERATION)
282 !====================================================================
283c
284 ! number of maximum newton iterations
285 niter = 3
286c
287 ! Loop over the iterations
288 DO iter = 1, niter
289#include "vectorize.inc"
290 ! Loop over yielding elements
291 DO ii=1,nindx
292 i = index(ii)
293c
294 ! Note : in this part, the purpose is to compute for each iteration
295 ! a plastic multiplier allowing to update internal variables to satisfy
296 ! the consistency condition using the backward Euler implicit method
297 ! with a Newton-Raphson iterative procedure
298 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
299 ! -> PHI : current value of yield function (known)
300 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
301 ! into account of internal variables kinetic :
302 ! plasticity, temperature and damage (to compute)
303c
304 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
305 !-------------------------------------------------------------
306 normxx = (gg*signxx(i) + hh*(signxx(i)-signyy(i)))/(max(sighl(i),em20))
307 normyy = (ff*signyy(i) + hh*(signyy(i)-signxx(i)))/(max(sighl(i),em20))
308 normxy = two*nn*signxy(i)/(max(sighl(i),em20))
309c
310 ! 2 - Computation of DPHI_DLAMBDA
311 !---------------------------------------------------------
312c
313 ! a) Derivative with respect stress increments tensor DSIG
314 ! --------------------------------------------------------
315 dfdsig2 = normxx * (a11*normxx + a12*normyy)
316 . + normyy * (a12*normxx + a22*normyy)
317 . + normxy * normxy * g12
318c
319 ! b) Derivatives with respect to plastic strain P
320 ! ------------------------------------------------
321c
322 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
323 ! ----------------------------------------------------------------------------
324 ! Already computed
325c
326 ! ii) Derivative of dPLA with respect to DLAM
327 ! -------------------------------------------
328 sig_dfdsig = signxx(i) * normxx
329 . + signyy(i) * normyy
330 . + signxy(i) * normxy
331 dpla_dlam(i) = sig_dfdsig / max(yld(i),em20)
332c
333 ! 3 - Computation of plastic multiplier and variables update
334 !----------------------------------------------------------
335c
336 ! Derivative of PHI with respect to DLAM
337 dphi_dlam(i) = - dfdsig2 - h(i)*dpla_dlam(i)
338 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
339c
340 ! Computation of the plastic multiplier
341 dlam = -phi(i)/dphi_dlam(i)
342c
343 ! Plastic strains tensor update
344 dpxx(i) = dlam * normxx
345 dpyy(i) = dlam * normyy
346 dpxy(i) = dlam * normxy
347c
348 ! Elasto-plastic stresses update
349 signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
350 signyy(i) = signyy(i) - (a22*dpyy(i) + a12*dpxx(i))
351 signxy(i) = signxy(i) - dpxy(i)*g12
352c
353 ! Cumulated plastic strain and strain rate update
354 ddep = dlam*dpla_dlam(i)
355 dpla(i) = max(zero, dpla(i) + ddep)
356 pla(i) = pla(i) + ddep
357c
358 ! Update Hill equivalent stress
359 sighl(i) = (ff+hh)*signyy(i)**2 + (gg+hh)*signxx(i)**2
360 . - two*hh*signxx(i)*signyy(i) + two*nn*signxy(i)**2
361 sighl(i) = sqrt(max(zero,sighl(i)))
362c
363 ! Transverse strain update
364 dpzz(i) = dpzz(i) - (dpxx(i)+dpyy(i))
365c
366 ! If the continuous hardening yield stress is chosen
367 IF (itab == 0) THEN
368 ! Update the hardening yield stress
369 yld(i) = sigy
370 . + qr1*(one - exp(-cr1*pla(i)))
371 . + qr2*(one - exp(-cr2*pla(i)))
372 h(i) = qr1*cr1*exp(-cr1*pla(i)) + qr2*cr2*exp(-cr2*pla(i))
373c
374 ! Update yield function value
375 phi(i) = sighl(i) - yld(i)
376 ENDIF
377c
378 ENDDO
379 ! End of the loop over the yielding elements
380c
381 ! If the tabulated yield stress is chosen
382 IF (itab == 1) THEN
383 ! -> No strain-rate effect
384 IF (mfunc == 1) THEN
385 ! Recovering position tables
386 ipos1(1:nel) = vartmp(1:nel,1)
387 iad1(1:nel) = npf(kfunc(1)) / 2 + 1
388 ilen1(1:nel) = npf(kfunc(1)+1) / 2 - iad1(1:nel) - ipos1(1:nel)
389 ! Computing the plastic strain interpolation
390 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
391 ! Storing the position
392 vartmp(1:nel,1) = ipos1(1:nel)
393 ! Computing the yield stress and its derivative
394 yld(1:nel) = yfac(1)*y1(1:nel)
395 h(1:nel) = yfac(1)*dydx1(1:nel)
396 ! Update the yield function
397 phi(1:nel) = sighl(1:nel) - yld(1:nel)
398 ! -> Strain-rate effect
399 ELSE
400 DO i=1,nel
401 j1 = jj(i)
402 j2 = j1 + 1
403 ! Recovering position tables for the first rate
404 ipos1(i) = vartmp(i,j1)
405 iad1(i) = npf(kfunc(j1)) / 2 + 1
406 ilen1(i) = npf(kfunc(j1)+1) / 2 - iad1(i) - ipos1(i)
407 ! Recovering position tables for the second rate
408 ipos2(i) = vartmp(i,j2)
409 iad2(i) = npf(kfunc(j2)) / 2 + 1
410 ilen2(i) = npf(kfunc(j2)+1) / 2 - iad2(i) - ipos2(i)
411 ENDDO
412 ! Computing the plastic strain interpolation
413 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
414 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
415 ! Storing new positions and computing yield stress
416 DO i=1,nel
417 j1 = jj(i)
418 j2 = j1+1
419 ! Hardening yield stress
420 y1(i) = y1(i)*yfac(j1)
421 y2(i) = y2(i)*yfac(j2)
422 fac = (epsp(i) - rate(j1))/(rate(j2) - rate(j1))
423 yld(i) = y1(i) + fac*(y2(i)-y1(i))
424 yld(i) = max(yld(i),em20)
425 ! Derivatives of yield stress
426 dydx1(i) = dydx1(i)*yfac(j1)
427 dydx2(i) = dydx2(i)*yfac(j2)
428 h(i) = dydx1(i)+fac*(dydx2(i)-dydx1(i))
429 ! New positions
430 vartmp(i,j1) = ipos1(i)
431 vartmp(i,j2) = ipos2(i)
432 ! Update the yield function
433 phi(i) = sighl(i) - yld(i)
434 ENDDO
435 ENDIF
436 ENDIF
437c
438 ENDDO
439 ! End of the loop over the iterations
440 !===================================================================
441 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM
442 !===================================================================
443C
444 ! Update and store new internal variables
445 DO i=1,nel
446 ! Equivalent stress
447 seq(i) = sighl(i)
448 ! Sound-speed
449 soundsp(i) = sqrt(max(a11,a22)/rho0(i))
450 ! Coefficient for hourglass
451 IF (dpla(i)>zero) THEN
452 hardm(i) = h(i)
453 etse(i) = h(i) / (h(i) + max(e11,e22))
454 ELSE
455 etse(i) = one
456 hardm(i) = zero
457 ENDIF
458 ! Computation of the thickness variation
459 invd33 = one/max(em20,d33)
460 deelzz(i) = -(nu13/e11)*(signxx(i)-sigoxx(i)) -(nu23/e22)*(signyy(i)-sigoyy(i))
461 IF ((inloc > 0).AND.(loff(i) == one)) THEN
462 normxx = (gg*signxx(i) + hh*(signxx(i)-signyy(i)))/(max(sighl(i),em20))
463 normyy = (ff*signyy(i) + hh*(signyy(i)-signxx(i)))/(max(sighl(i),em20))
464 normxy = two*nn*signxy(i)/(max(sighl(i),em20))
465 dpzz(i) = -max(dplanl(i),zero)*(normxx + normyy)
466 ENDIF
467 dezz(i) = deelzz(i) + dpzz(i)
468 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
469 ! Plastic strain-rate filtering (if needed)
470 IF ((itab == 1).AND.(mfunc > 1).AND.(vp == 1)) THEN
471 dpdt = dpla(i)/max(em20,timestep(i))
472 uvar(i,1) = asrate * dpdt + (one - asrate) * uvar(i,1)
473 epsp(i) = uvar(i,1)
474 ENDIF
475 ENDDO
476 IF (ALLOCATED(rate)) DEALLOCATE(rate)
477 IF (ALLOCATED(yfac)) DEALLOCATE(yfac)
478C
my_real
#define my_real
Definition cppsort.cpp:32
the
end diagonal values have been computed in the(sparse) matrix id.SOL
max
#define max(a, b)
Definition macros.h:21
tabulated
subroutine tabulated(iflag, nel, pm, off, eint, mu, espe, dvol, df, vnew, mat, psh, pnew, dpdm, dpde, npf, tf)
Definition tabulated.F:32
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72