OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m36iter_imp.F File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine m36iter_imp (signxx, signyy, signzz, signxy, signyz, signzx, pla, yld, g3, yfac, dpla1, h, tf, iad1, ilen1, nel, fisokin, vartmp, pla0, plam, y1, dydx1, ipflag, pfun, ipfun, iposp, nrate, npf, iadp, ilenp, pfac, pscale1, dfdp, fail, nvartmp, sigbxx, sigbyy, sigbzz, sigbxy, sigbyz, sigbzx)

Function/Subroutine Documentation

◆ m36iter_imp()

subroutine m36iter_imp ( signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
pla,
yld,
g3,
yfac,
dpla1,
h,
tf,
integer, dimension(mvsiz) iad1,
integer, dimension(mvsiz) ilen1,
integer nel,
fisokin,
integer, dimension(nel,nvartmp) vartmp,
pla0,
plam,
y1,
dydx1,
integer ipflag,
integer pfun,
integer ipfun,
integer, dimension(mvsiz) iposp,
integer nrate,
integer, dimension(*) npf,
integer, dimension(mvsiz) iadp,
integer, dimension(mvsiz) ilenp,
pfac,
pscale1,
dfdp,
fail,
integer nvartmp,
sigbxx,
sigbyy,
sigbzz,
sigbxy,
sigbyz,
sigbzx )

Definition at line 33 of file m36iter_imp.F.

46c--------------------------------------------------------------
47c
48c This subroutine computes plastic strain rate multiplier
49c "delta_lambda" for the von Mises material characterized by
50c piecewise-linear hardening.
51c In addition, it computes the deviatoric stresses
52c via radial return.
53c
54c -------------------
55c -- parameters in --
56c -------------------
57c
58c SIGNXX = Elastic Predictor (on input)
59c (XX-Component Deviatoric Stress Tensor)
60c (...)
61c SIGNZX = Elastic Predictor (on input)
62c (ZX-Component Deviatoric Stress Tensor)
63c
64c PLA = previous time-step equivalent plastic strain (on input)
65c YLD = previous time-step yield stress (on input)
66c G3 = 3 x Kirchhoff's modulus (elastic shear modulus))
67c YFAC = scaling factor of sigma-eps_plast relation
68c
69c TF = array containing discrete sigma-eps_plast relation
70c IAD1, ILEN1 = pre-computed data required by interpolation
71c routine "VINTER"
72c NEL = size of the element group being processed
73c FISOKIN = indicator for isotropic/kinematic/moixed hardening
74c UVAR = array containing various state variables
75c IPFLAG = flag for existence of pressure scaling (scalar)
76c PFUN = flag for existence of pressure scaling (array)
77c PFAC = scaling factor for pressure at time t (not needed)
78c PSCALE1 = pressure at time t + delta_t
79c IPOSP, ILENP, NRATE, NPF, IADP = arrays and parameters required
80c by pressure interpolating routine (for pressure scaling)
81c
82c
83c --------------------
84c -- parameters out --
85c --------------------
86c
87c DPLA1 = equivalent plastic strain increment
88c PLA = updated equivalent plastic strain (on output)
89c H = updated hardening parameter
90c SIGNXX = updated XX-component of deviatoric stress tensor
91c (...)
92c SIGNZX = updated ZX-component of deviatoric stress tensor,
93c UVAR = array containing various state variables
94c
95c -----------------
96c -- work arrays --
97c -----------------
98c
99c Y1 = yield stress
100c DYDX1 = derivative in sig0-eps_plast curve
101c
102C-----------------------------------------------
103C I m p l i c i t T y p e s
104C-----------------------------------------------
105#include "implicit_f.inc"
106C-----------------------------------------------
107C G l o b a l P a r a m e t e r s
108C-----------------------------------------------
109#include "mvsiz_p.inc"
110C-----------------------------------------------
111C C o m m o n B l o c k s
112C-----------------------------------------------
113#include "parit_c.inc"
114#include "scr05_c.inc"
115#include "units_c.inc"
116C-----------------------------------------------
117C D u m m y A r g u m e n t s
118C-----------------------------------------------
119C
120 INTEGER NEL, NVARTMP
121 INTEGER IAD1(MVSIZ), ILEN1(MVSIZ),
122 . IPFLAG, PFUN, IPFUN, IPOSP(MVSIZ),
123 . IADP(MVSIZ), ILENP(MVSIZ), NRATE, NPF(*)
124 INTEGER :: VARTMP(NEL,NVARTMP)
125C
126C REAL
127 my_real
128 . signxx(*), signyy(*), signzz(*),
129 . signxy(*), signyz(*), signzx(*),
130 . pla(*), g3(*), yld(*),
131 . dpla1(*), h(*),
132 . tf(*), y1(*), dydx1(*),
133 . yfac(mvsiz,*),
134 . pla0(mvsiz), plam(mvsiz),
135 . fisokin,
136 . pfac(mvsiz), pscale1(mvsiz), dfdp(mvsiz),
137 . fail(mvsiz)
138 my_real
139 . sigbxx(nel),sigbyy(nel),sigbzz(nel),sigbxy(nel),sigbyz(nel),sigbzx(nel)
140CC
141C-----------------------------------------------
142C L o c a l V a r i a b l e s
143C-----------------------------------------------
144 INTEGER I, II, ITER, NITER, I_BILIN, IPOS
145C
146C REAL
147 my_real
148 . xsi, dxsi, lhs, rhs, alpha_radial, vm, hep, r, dpla,
149 . yld_curr, yld_prev, yld_m, h_m, p_scal, y_scal,
150 . total_scal, yld_m_prev, h_m_prev,
151 . sxx, syy, szz, sxy, syz, szx
152
153C
154C-----------------------------------------------
155C External Functions
156C-----------------------------------------------
157C
158 my_real
159 . finter2
160 EXTERNAL finter2
161C
162C=======================================================================
163C
164C
165c ---- hardening law flag: i_bilin = 1 -- bilinear hardening law
166c ---- i_bilin = 0 -- general nonlinear hardening
167c
168 i_bilin = 0
169c
170c ---- max number of iterations ----
171 niter = 20
172c
173c -- scaling if pressure dependent yield function --
174c -- compute the scaling factor for pressure + delta-pressure
175c
176 DO i=1,nel
177 pfac(i) = one
178 ENDDO
179c
180 IF (ipflag == nel.AND.(iparit == 0)) THEN
181 DO i=1,nel
182 iposp(i) = vartmp(i,2)
183 iadp(i) = npf(ipfun) / 2 + 1
184 ilenp(i) = npf(ipfun + 1) / 2 - iadp(i) - iposp(i)
185 ENDDO
186 IF (iresp == 1) THEN
187 CALL vinter2dp(tf,iadp,iposp,ilenp,nel,pscale1,dfdp,pfac)
188 ELSE
189 CALL vinter2(tf,iadp,iposp,ilenp,nel,pscale1,dfdp,pfac)
190 ENDIF
191 pfac(i) = max(zero, pfac(i))
192 ELSEIF (ipflag > 0) THEN
193 IF (pfun > 0) THEN
194 DO i=1,nel
195 iposp(i) = vartmp(i,2)
196 iadp(i) = npf(ipfun) / 2 + 1
197 ilenp(i) = npf(ipfun + 1) / 2 - iadp(i) - iposp(i)
198 pfac(i) = finter2(tf ,iadp(i),iposp(i),ilenp(i),
199 . pscale1(i),dfdp(i))
200 pfac(i) = max(zero, pfac(i))
201 ENDDO
202 ENDIF
203 ENDIF ! .. IPFLAG ...
204c -- end of computation of pressure-scaling factor --
205c
206c
207c -- store equivalent plastic strain at the beginning of time-step
208c -- store "fractional" equivalent plastic strain
209c
210 DO i=1,nel
211 pla0(i) = pla(i)
212 plam(i) = (one - fisokin) * pla(i)
213 ENDDO
214c
215c -- computes the yield stress --
216c
217 DO i=1,nel
218C 2Andrzej--------!!!IPOS is not initialized-----
219 ipos = 0
220 yld(i) = finter2( tf, iad1(i), ipos, ilen1(i),
221 $ plam(i),dydx1(i) )
222C------------------same than explicit
223 yld(i) = max(yld(i),em20)
224 ENDDO
225c
226 DO 100 i = 1, nel ! .. loop on a group of elements ..
227c
228c ... temporarily stores elastic predictor
229 sxx = signxx(i)
230 syy = signyy(i)
231 szz = signzz(i)
232 sxy = signxy(i)
233 syz = signyz(i)
234 szx = signzx(i)
235c
236 IF(i_bilin/=1 )THEN !.. general, nonlinear hardening
237c
238 xsi = zero ! .. plastic strain rate multiplier ..
239c
240c -- PREVIOUS TIME-STEP EQUIVALENT PLASTIC STRAIN
241 pla(i) = max(zero, pla(i))
242c
243c -- von Mises stress at elastic predictor --
244 vm =three*(half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
245 $ +signxy(i)**2+signyz(i)**2+signzx(i)**2)
246 vm =sqrt(vm)
247c
248c
249c -- COMPARE VM with PREVIOUS TIME-STEP YIELD STRESS --
250c -- YLD is already multiplied by PFAC (at pressure)
251c
252 IF( vm > yld(i) )then ! .. plastically active process
253c
254c
255c -- scaling for sigma-eps_plast function --
256 y_scal = yfac(i,1)
257c
258c -- scaling if pressure dependent yield function --
259 p_scal = pfac(i) * fail(i)
260c
261c -- total scaling
262 total_scal = y_scal * p_scal
263c
264c ------------------------------------------------------
265c yield stress and hardening parameter at m*(eps)
266c m =: (1 - FISOKIN)
267c ------------------------------------------------------
268c -- compute yield stress and hardening parameter DYDX1
269c at equivalent plastic strain PLAM, via interpolation
270c
271 ipos = 0 !.. FINTER2 will search the whole table
272c
273c -- yield stress at m*(eps) --
274 yld_m_prev = finter2( tf, iad1(i), ipos, ilen1(i),
275 $ plam(i),dydx1(i) )
276c
277c -- hardening parameter --
278 h_m_prev = dydx1(i)
279c
280c -- apply scaling --
281 yld_m_prev = yld_m_prev * total_scal
282 h_m_prev = h_m_prev * total_scal
283 yld_m_prev = max(yld_m_prev,em20)
284c -- check if H is nonnegative and y0 positive
285 if( .not.( h_m_prev >= zero ).or.
286 $ .not.( yld_m_prev > zero ) )then
287 write(istdo,1000)plam(i),h_m_prev,yld_m_prev
288 write(iout,1000)plam(i),h_m_prev,yld_m_prev
289 CALL imp_stop(-1)
290 endif
291c
292c -------------------
293c yield stress at eps
294c -------------------
295c
296 ipos = 0 !.. FINTER2 will search the whole table
297c
298c -- yield stress at eps --
299c
300 yld_prev = finter2( tf, iad1(i), ipos, ilen1(i),
301 $ pla0(i),dydx1(i) )
302c
303c -- apply scaling --
304 yld_prev = yld_prev * total_scal
305 yld_prev = max(yld_prev,em20)
306c
307c -- check if H is nonnegative and y0 positive
308 if( .not.( dydx1(i) >= zero ).or.
309 $ .not.( yld_prev > zero ) )then
310 write(istdo,1000)pla0(i),dydx1(i),yld_prev
311 write(iout,1000)pla0(i),dydx1(i),yld_prev
312 CALL imp_stop(-1)
313 endif
314c
315c
316c
317c =================
318c NR iterations
319c =================
320c
321 DO 80 iter = 1,niter
322c
323c ---------------------------------------------------
324c yield stress and hardening parameter at eps + deps
325c ---------------------------------------------------
326C -- compute yield stress and hardening parameter DYDX1
327c at equivalent plastic strain PLA, via interpolation
328c
329 ipos = 0 !.. FINTER2 will search the whole table
330c
331c -- yield stress at eps + deps --
332c
333 yld_curr = finter2( tf,iad1(i), ipos, ilen1(i),
334 $ pla(i),dydx1(i) )
335c
336c -- hardening parameter --
337 h(i) = dydx1(i)
338c
339c -- apply scaling --
340 yld_curr = yld_curr * total_scal
341 h(i) = h(i) * total_scal
342 yld_curr = max(yld_curr,em20)
343c
344c -- check if H is nonnegative and y0 positive
345 if( .not.( h(i) >= zero ).or.
346 $ .not.( yld_curr > zero ) )then
347 write(istdo,1000)pla(i),h(i),yld_curr
348 write(iout,1000)pla(i),h(i),yld_curr
349 CALL imp_stop(-1)
350 endif
351c
352c
353c -----------------------------------------
354c RHS and LHS for Newton at the Gauss point
355c -----------------------------------------
356c
357 IF(fisokin == zero)then
358 rhs = vm - g3(i) * xsi - yld_curr
359 lhs = g3(i) + h(i)
360 ELSEIF(fisokin == one)then
361 rhs = vm - g3(i)*xsi + yld_prev - yld_curr - yld_m_prev
362 lhs = g3(i) + h(i)
363 ELSE
364 rhs = vm - g3(i)*xsi + yld_prev - yld_curr - yld_m_prev
365 lhs = g3(i) + h(i)
366 ENDIF
367c
368c
369c --------------------------------------
370c Solution for Newton at the Gauss point
371c --------------------------------------
372 dxsi = rhs/lhs
373C -- update plastic strain rate multiplier
374 xsi = xsi + dxsi
375C -- update equivalent plastic strain --
376 pla(i) = pla(i) + dxsi
377 pla(i) = max(zero,pla(i)) !..to prevent negative value
378c
379c
380c -- convergence criterion should be stringent --
381 IF( abs(dxsi)<em10 ) GO TO 90
382c
383 80 CONTINUE ! -- ITER --
384c -- we need a warning here --
385ccc WRITE(*,*)
386ccc $ 'M36ITER_IMP--NON-CONVERGED ITERATION', ABS(DXSI),ABS(RHS)
387c
388 90 CONTINUE
389c
390c ... equivalent plastic strain increment
391 dpla = max(zero,xsi) !..to prevent negative value
392c
393c
394c
395c --------------------------------------
396c update of stresses and back stresses
397c --------------------------------------
398c
399 IF(fisokin == zero)then
400 alpha_radial = one - ( g3(i)/max(vm,em20) )*dpla
401c
402 signxx(i) = signxx(i) * alpha_radial
403 signyy(i) = signyy(i) * alpha_radial
404 signzz(i) = signzz(i) * alpha_radial
405 signxy(i) = signxy(i) * alpha_radial
406 signyz(i) = signyz(i) * alpha_radial
407 signzx(i) = signzx(i) * alpha_radial
408c
409 ELSEIF(fisokin > zero.and.fisokin<=one)then
410c
411c -- update "fractional" equivalent plastic strain --
412 plam(i) = (one - fisokin)*pla(i) !. for mixed hardeneing
413c
414c
415c ------------------------------------------------------
416c yield stress and hardening parameter at m*(eps+deps)
417c m =: (1 - FISOKIN)
418c ------------------------------------------------------
419c -- compute yield stress and hardening parameter DYDX1
420c at equivalent plastic strain PLAM, via interpolation
421c
422 ipos = 0 !.. FINTER2 will search the whole table
423c
424c -- yield stress at m*(eps+deps) --
425 yld_m = finter2( tf, iad1(i), ipos, ilen1(i),
426 $ plam(i),dydx1(i) )
427c
428c -- hardening parameter --
429 h_m = dydx1(i)
430c -- apply scaling --
431 yld_m = yld_m * total_scal
432 h_m = h_m * total_scal
433 yld_m = max(yld_m,em20)
434c
435c -- check if H is nonnegative and y0 positive
436 if( .not.( h_m >= zero ).or.
437 $ .not.( yld_m > zero ) )then
438 write(istdo,1000)plam(i),h_m,yld_m
439 write(iout,1000)plam(i),h_m,yld_m
440 CALL imp_stop(-1)
441 endif
442c
443c
444c
445c ..updates (shifted) stresses
446c
447c -- ALPHA_RADIAL := PARAMETER ALPHA OF THE RADIAL RETURN
448c
449 alpha_radial = one - ( g3(i)/max(vm,em20) )*dpla
450 $ - (yld_curr - yld_prev -yld_m + yld_m_prev )/max(vm,em20)
451c
452 signxx(i) = signxx(i) * alpha_radial
453 signyy(i) = signyy(i) * alpha_radial
454 signzz(i) = signzz(i) * alpha_radial
455 signxy(i) = signxy(i) * alpha_radial
456 signyz(i) = signyz(i) * alpha_radial
457 signzx(i) = signzx(i) * alpha_radial
458c
459c ..updates back stresses
460 hep = ( yld_curr - yld_prev
461 $ -yld_m + yld_m_prev )/max(vm,em20)
462c
463 sigbxx(i) = sigbxx(i) + sxx *hep
464 sigbyy(i) = sigbyy(i) + syy *hep
465 sigbzz(i) = sigbzz(i) + szz *hep
466 sigbxy(i) = sigbxy(i) + sxy *hep
467 sigbyz(i) = sigbyz(i) + syz *hep
468 sigbzx(i) = sigbzx(i) + szx *hep
469C
470c --adds the back stresses to shifted stresses
471 signxx(i) = signxx(i) + sigbxx(i)
472 signyy(i) = signyy(i) + sigbyy(i)
473 signzz(i) = signzz(i) + sigbzz(i)
474 signxy(i) = signxy(i) + sigbxy(i)
475 signyz(i) = signyz(i) + sigbyz(i)
476 signzx(i) = signzx(i) + sigbzx(i)
477c
478 ENDIF ! ..FISOKIN > ZERO.and.FISOKIN<=ONE ..
479c
480c
481c .. updates yield stress.
482 IF(fisokin == zero)then
483 yld(i) = yld_curr
484 ELSEIF(fisokin == one)then
485cccc YLD(I) = YLD(I)
486 ELSE
487 yld(i) = yld_m
488 ENDIF
489c
490 dpla1(i) = dpla !--PLASTIC STRAIN MULTIPLIER "delta lambda"
491
492c
493c
494c --- end of plastically active process ---
495c
496 ELSE !.. elasticity
497
498c
499c .. gets stresses from shifted stresses and back stresses
500 IF(fisokin > zero)then
501 signxx(i) = signxx(i) + sigbxx(i)
502 signyy(i) = signyy(i) + sigbyy(i)
503 signzz(i) = signzz(i) + sigbzz(i)
504 signxy(i) = signxy(i) + sigbxy(i)
505 signyz(i) = signyz(i) + sigbyz(i)
506 signzx(i) = signzx(i) + sigbzx(i)
507 ENDIF
508
509 ENDIF ! .. VM > YLD(I)
510c
511c-------------------------------------------------------------------
512 ELSE ! -- i_bilin == 1
513c --- bilinear plastic hardening (isotropic, kinematic and mixed)
514c-------------------------------------------------------------------
515c
516c .. elastic predictor is temporarily stored in sxx, ... , szx
517c
518c .. von Mises stress**2 at the elastic predictor ..
519 vm =three*(half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
520 $ +signxy(i)**2+signyz(i)**2+signzx(i)**2)
521c
522 IF(vm > yld(i)*yld(i))THEN !.. plastically active process
523 vm =sqrt(vm)
524 r = yld(i)/ max(vm,em20)
525c .. plastic strain increment.
526 dpla=(one - r)*vm/max(g3(i)+h(i),em20)
527c
528c ... updates stresses ..
529 alpha_radial = one - g3(i)*( (one - r)/( g3(i) + h(i) ) )
530c
531 signxx(i) = signxx(i) * alpha_radial
532 signyy(i) = signyy(i) * alpha_radial
533 signzz(i) = signzz(i) * alpha_radial
534 signxy(i) = signxy(i) * alpha_radial
535 signyz(i) = signyz(i) * alpha_radial
536 signzx(i) = signzx(i) * alpha_radial
537c
538 IF(fisokin > zero)then
539c
540c --adds the back stresses at the beginning of the increment
541 signxx(i) = signxx(i) + sigbxx(i)
542 signyy(i) = signyy(i) + sigbyy(i)
543 signzz(i) = signzz(i) + sigbzz(i)
544 signxy(i) = signxy(i) + sigbxy(i)
545 signyz(i) = signyz(i) + sigbyz(i)
546 signzx(i) = signzx(i) + sigbzx(i)
547c
548c ..updates back stresses
549 hep = fisokin*h(i)* dpla/vm
550c
551 sigbxx(i) = sigbxx(i) + sxx *hep
552 sigbyy(i) = sigbyy(i) + syy *hep
553 sigbzz(i) = sigbzz(i) + szz *hep
554 sigbxy(i) = sigbxy(i) + sxy *hep
555 sigbyz(i) = sigbyz(i) + syz *hep
556 sigbzx(i) = sigbzx(i) + szx *hep
557c
558 ENDIF ! .. FISOKIN > ZERO ..
559c
560c .. updated yield stress.
561 yld(i) = max(yld(i)+(one - fisokin)*dpla*h(i),zero)
562c
563c .. updates equivalent plastic strain
564 pla(i) = pla(i) + dpla
565c
566 dpla1(i) = dpla
567c
568 ELSE !.. elasticity
569c
570 IF(fisokin > zero)then
571c .. gets stresses from shifted stresses and back stresses
572 signxx(i) = signxx(i) + sigbxx(i)
573 signyy(i) = signyy(i) + sigbyy(i)
574 signzz(i) = signzz(i) + sigbzz(i)
575 signxy(i) = signxy(i) + sigbxy(i)
576 signyz(i) = signyz(i) + sigbyz(i)
577 signzx(i) = signzx(i) + sigbzx(i)
578 ENDIF
579c
580 ENDIF ! .. VM > YLD(I)*YLD(I)
581c
582c
583 ENDIF ! .. i_bilin ..
584c
585 100 CONTINUE !..LOOP ON A GROUP OF ELEMENTS
586c
587 1000 FORMAT(3x,'Hardening parameters in law36 is not positive'/,
588 $ 2x,'Eps_plast =',e11.4,2x,'H =',e11.4,
589 $ 2x,'Sy0 =',e11.4)
590c-----------
591 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
imp_stop
subroutine imp_stop(istop)
Definition imp_solv.F:1992
max
#define max(a, b)
Definition macros.h:21
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143
vinter2dp
subroutine vinter2dp(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:214