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

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps93 (nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, pla, uvar, rho0, off, et, yld, seq, epsp, asrate, nvartmp, vartmp, dpla)

Function/Subroutine Documentation

◆ sigeps93()

subroutine sigeps93 ( integer nel,
integer nuparam,
integer nuvar,
integer mfunc,
integer, dimension(mfunc) kfunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
epspxx,
epspyy,
epspzz,
epspxy,
epspyz,
epspzx,
depsxx,
depsyy,
depszz,
depsxy,
depsyz,
depszx,
sigoxx,
sigoyy,
sigozz,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
soundsp,
pla,
uvar,
rho0,
off,
et,
yld,
seq,
epsp,
asrate,
integer nvartmp,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
dpla )

Definition at line 32 of file sigeps93.F.

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