OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps36.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| sigeps36 ../engine/source/materials/mat/mat036/sigeps36.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
28!||--- calls -----------------------------------------------------
29!|| finter ../engine/source/tools/curve/finter.F
30!|| m36iter_imp ../engine/source/materials/mat/mat036/m36iter_imp.F
31!|| vinter ../engine/source/tools/curve/vinter.f
32!|| vinter2 ../engine/source/tools/curve/vinter.F
33!||====================================================================
34 SUBROUTINE sigeps36(
35 1 NEL ,NUVAR ,MFUNC ,KFUNC ,NPF ,
36 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,
37 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 5 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
39 6 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 7 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 8 SOUNDSP,VISCMAX,UVAR ,OFF ,NGL ,IEOS ,
42 9 IPM ,MAT ,EPSP ,IPLA ,YLD ,PLA ,
43 A DPLA1 ,ETSE ,AL_IMP ,SIGNOR ,AMU ,DPDM ,
44 B FACYLDI,NVARTMP,VARTMP ,DMG ,INLOC ,PLANL ,
45 C SIGBXX ,SIGBYY ,SIGBZZ ,SIGBXY ,SIGBYZ ,SIGBZX )
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50#include "comlock.inc"
51C-----------------------------------------------
52C G l o b a l P a r a m e t e r s
53C-----------------------------------------------
54#include "mvsiz_p.inc"
55C---------+---------+---+---+--------------------------------------------
56C VAR | SIZE |TYP| RW| DEFINITION
57C---------+---------+---+---+--------------------------------------------
58C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
59C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
60C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
61C---------+---------+---+---+--------------------------------------------
62C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
63C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
64C NPF | * | I | R | FUNCTION ARRAY
65C TF | * | F | R | FUNCTION ARRAY
66C---------+---------+---+---+--------------------------------------------
67C TIME | 1 | F | R | CURRENT TIME
68C TIMESTEP| 1 | F | R | CURRENT TIME STEP
69C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
70C RHO0 | NEL | F | R | INITIAL DENSITY
71C RHO | NEL | F | R | DENSITY
72C ... | | | |
73C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
74C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
75C ... | | | |
76C EPSXX | NEL | F | R | STRAIN XX
77C EPSYY | NEL | F | R | STRAIN YY
78C ... | | | |
79C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
80C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
81C ... | | | |
82C---------+---------+---+---+--------------------------------------------
83C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
84C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
85C ... | | | |
86C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
87C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
88C---------+---------+---+---+--------------------------------------------
89C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
90C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
91C---------+---------+---+---+--------------------------------------------
92#include "param_c.inc"
93#include "scr17_c.inc"
94#include "com01_c.inc"
95#include "com04_c.inc"
96#include "com08_c.inc"
97#include "units_c.inc"
98#include "impl1_c.inc"
99C=======================================================================
100 INTEGER, INTENT(IN) :: NEL, NUVAR,IPLA,NVARTMP,INLOC
101 INTEGER, INTENT(IN) :: IEOS
102 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL,MAT
103 INTEGER ,DIMENSION(NPROPMI,NUMMAT), INTENT(IN) :: IPM
104 my_real
105 . TIME,TIMESTEP,UPARAM(*),
106 . RHO0(NEL),
107 . EPSXX(NEL),EPSYY(NEL),EPSZZ(NEL),
108 . EPSXY(NEL),EPSYZ(NEL),EPSZX(NEL),
109 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
110 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
111 . sigoxx(nel),sigoyy(nel),sigozz(nel),
112 . sigoxy(nel),sigoyz(nel),sigozx(nel),
113 . epsp(nel),etse(nel),amu(nel),facyldi(nel),
114 . dmg(nel),planl(nel)
115 my_real, DIMENSION(MVSIZ) ,INTENT(IN) :: dpdm
116C-----------------------------------------------
117C O U T P U T A r g u m e n t s
118C-----------------------------------------------
119 my_real
120 . signxx(nel),signyy(nel),signzz(nel),
121 . signxy(nel),signyz(nel),signzx(nel),
122 . soundsp(nel),viscmax(nel),dpla1(nel),
123 . al_imp(nel),signor(mvsiz,6)
124C-----------------------------------------------
125C I N P U T O U T P U T A r g u m e n t s
126C-----------------------------------------------
127 my_real
128 . uvar(nel,nuvar), off(nel), yld(nel),
129 . pla(nel)
130 INTEGER, INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
131 my_real, DIMENSION(NEL), INTENT(INOUT) :: SIGBXX,SIGBYY,SIGBZZ,SIGBXY,SIGBYZ,SIGBZX
132C-----------------------------------------------
133C VARIABLES FOR FUNCTION INTERPOLATION
134C-----------------------------------------------
135 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
136 my_real :: TF(*),FINTER
137 EXTERNAL FINTER
138C-----------------------------------------------
139C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
140C Y : y = f(x)
141C X : x
142C DYDX : f'(x) = dy/dx
143C IFUNC(J): FUNCTION INDEX
144C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
145C NPF,TF : FUNCTION PARAMETER
146C-----------------------------------------------
147C L o c a l V a r i a b l e s
148C-----------------------------------------------
149 INTEGER I,J,K,J1,J2,NFUNC,VP,II,NITER,IFAIL,NRATE,PFUN,IPFUN,
150 . IPFLAG,ISMOOTH,YLDCHECK,OPTE,IFUNCE,NINDX
151 INTEGER ILEN1(MVSIZ),IAD2(MVSIZ),ILEN2(MVSIZ),INDEX(MVSIZ),
152 . IFUNC(100),IAD1(MVSIZ), IADP(MVSIZ),ILENP(MVSIZ),
153 . INDX(MVSIZ),JJ(MVSIZ), IADP2(MVSIZ),ILENP2(MVSIZ)
154 INTEGER, DIMENSION(MVSIZ) :: IPOS1,IPOS2,IPOSPE,IPOSP,IPOSP2
155 my_real
156 . E11,NU,DAV,VM,R(MVSIZ),FAC,EPST,P,EPSP1,EPSP2,
157 . E1,E2,E3,E4,E5,E6,C,CC,D,Y,YP,E42,E52,E62,EPST2,
158 . C11,G11,G21,G31,EPSR1,DPLA,EPSF,HKIN,FISOKIN,
159 . DSXX,DSYY,DSZZ,DSXY,DSYZ,
160 . dszx,sigpxx,sigpyy,sigpxy,sigpyz,sigpzx,sigpzz,alpha,
161 . epsr1dav,epsr1f,vm_1,g3h,norm_1,epsmax,ssp,
162 . einf,ce1,epsr2,svm,df,f
163 my_real :: fact(nel),
164 . bulk(mvsiz),g(mvsiz),g2(mvsiz),g3(mvsiz),
165 . y1(mvsiz),y2(mvsiz),h(mvsiz),dydx1(mvsiz),
166 . dydx2(mvsiz),fail(mvsiz),e(mvsiz),
167 . p0(mvsiz),pfac(mvsiz),pscale(mvsiz),rfac(mvsiz),
168 . pscale1(mvsiz),pla0(mvsiz), plam(mvsiz),
169 . dfdp(mvsiz),epstt(mvsiz),
170 . sigexx(mvsiz),sigeyy(mvsiz),sigezz(mvsiz),
171 . sigexy(mvsiz),sigeyz(mvsiz),sigezx(mvsiz),
172 . dydxe(mvsiz),plaold(mvsiz),
173 . yfac(mvsiz,2),
174 . coef(mvsiz),
175 . dfsr(mvsiz)
176 my_real
177 . plap(mvsiz),svm2(mvsiz),dpla_j(mvsiz),dpla_i(mvsiz),
178 . hi(mvsiz),escale12(mvsiz),
179 . svm_tab(mvsiz)
180 INTEGER :: NINDEX_PLA,NINDEX_VINTER
181 INTEGER, DIMENSION(MVSIZ) :: INDEX_PLA,INDEX_VINTER
182 my_real, DIMENSION(MVSIZ) :: TAB_LOC,Y1_LOC,DYDX1_LOC, Y2_LOC,DYDX2_LOC
183C=======================================================================
184C USER VARIABLES INITIALIZATION
185C-----------------------------------------------
186 NITER = 5
187 nindex_pla = 0
188 nfunc = ipm(10,mat(1))
189 DO j=1,nfunc
190 ifunc(j)= ipm(10+j,mat(1))
191 ENDDO
192 ipfun = ifunc(nfunc-1)
193C
194 e11 = uparam(2)
195 g11 = uparam(5)
196 nu = uparam(6)
197 nrate = nint(uparam(1))
198 epsmax = uparam(6+2*nrate+1)
199 epsr1 = uparam(6+2*nrate+2)
200 epsr2 = uparam(6+2*nrate+3)
201 g21 = uparam(6+2*nrate+4)
202 g31 = uparam(6+2*nrate+5)
203 c11 = uparam(6+2*nrate+6) ! ..bulk modulus
204 ssp = uparam(6+2*nrate+7)
205 fisokin = uparam(6+2*nrate+8)
206 epsf = uparam(6+2*nrate+9)
207 pfun = nint(uparam(16+2*nrate))
208 opte = uparam(2*nrate+23)
209 einf = uparam(2*nrate+24)
210 ce1 = uparam(2*nrate+25)
211 vp = nint(uparam(2*nrate + 26)) ! flag for plastic strain dependency
212 ifail = nint(uparam(2*nrate + 27)) ! flag rupture
213 yldcheck = nint(uparam(2*nrate + 28))
214 ismooth = nint(uparam(2*nrate + 29)) ! strain rate interpolation flag
215C
216#include "vectorize.inc"
217 DO i=1,nel
218 pfac(i) = one
219 pscale(i) = uparam(17+2*nrate)
220 e(i) = e11
221 g(i) = g11
222 g2(i) = g21
223 g3(i) = g31
224 bulk(i) = c11
225 soundsp(i) = ssp
226 ENDDO
227c-----------------------------
228
229C------------------------------------------
230C ECROUISSAGE CINE
231C------------------------------------------
232 IF (fisokin > zero) THEN
233c ..subtracts back stresses for kinmatic or mixed hardening
234c from the old stresses, thus, SIGOxx is a shifted old stress
235#include "vectorize.inc"
236 DO i=1,nel
237 sigoxx(i) = sigoxx(i) - sigbxx(i)
238 sigoyy(i) = sigoyy(i) - sigbyy(i)
239 sigozz(i) = sigozz(i) - sigbzz(i)
240 sigoxy(i) = sigoxy(i) - sigbxy(i)
241 sigoyz(i) = sigoyz(i) - sigbyz(i)
242 sigozx(i) = sigozx(i) - sigbzx(i)
243 ENDDO
244 ENDIF
245c------------------------------------------
246c element length scale function
247c------------------------------------------
248 IF (opte == 1) THEN
249 ifunce = uparam( 2*nrate + 22)
250 DO i=1,nel
251 IF(pla(i) > zero)THEN
252 nindex_pla = nindex_pla + 1
253 index_pla(nindex_pla) = i
254 ipospe(nindex_pla) = vartmp(i,1)
255 iadp(nindex_pla) = npf(kfunc(ifunce)) / 2 + 1
256 ilenp(nindex_pla) = npf(kfunc(ifunce)+1) / 2 - iadp(nindex_pla) - ipospe(nindex_pla)
257 tab_loc(nindex_pla) = pla(i)
258 ENDIF
259 ENDDO
260 CALL vinter2(tf,iadp,ipospe,ilenp,nindex_pla,tab_loc,dydxe,escale12)
261 vartmp(index_pla(1:nindex_pla),1) = ipospe(1:nindex_pla)
262#include "vectorize.inc"
263 DO ii=1,nindex_pla
264 i = index_pla(ii)
265 e(i) = escale12(ii)* e(i)
266 g(i) = half*e(i)/(one+nu)
267 g2(i) = two*g(i)
268 g3(i) = three*g(i)
269 bulk(i) = e(i)/three/(one - two*nu)
270 soundsp(i) = sqrt((bulk(i) + four_over_3*g(i))/rho0(i))
271 ENDDO
272 ELSEIF ( ce1 /= zero) THEN
273#include "vectorize.inc"
274 DO i=1,nel
275 IF(pla(i) > zero)THEN
276 e(i) = e11-(e11-einf)*(one-exp(-ce1*pla(i)))
277 g(i) = half*e(i)/(one+nu)
278 g2(i) = two*g(i)
279 g3(i) = three*g(i)
280 bulk(i) = e(i)/three/(one - two*nu)
281 soundsp(i) = sqrt((bulk(i) + four_over_3*g(i))/rho0(i))
282 ENDIF
283 ENDDO
284 ENDIF
285c------------------------------------------
286c PRESSURE DEPENDENT YIELD FUNCTION FACTOR
287c------------------------------------------
288 IF (pfun > 0) THEN
289 DO i=1,nel
290 iadp(i) = npf(ipfun) / 2 + 1
291 ilenp(i) = npf(ipfun + 1) / 2 - iadp(i) - vartmp(i,2)
292 ENDDO
293 CALL vinter2(tf,iadp,vartmp(1:nel,2),ilenp,nel,pscale,dfdp,pfac)
294 pfac(1:nel) = max(zero, pfac(1:nel))
295 ENDIF
296c
297c ..pressure + delta_pressure (scaled) ..
298 IF (impl_s > 0) pscale1(1:nel)=pscale(1:nel) * ( bulk(1:nel) * amu(1:nel) )
299c
300c------------------------------------------
301c PRESSURE DEPENDENT YIELD FUNCTION FACTOR
302c------------------------------------------
303#include "vectorize.inc"
304 DO i=1,nel
305C
306c ..strain increment axiator ..
307 dav = (depsxx(i)+depsyy(i)+depszz(i))*third
308c ..pressure..
309 p0(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third
310C
311c ..pressure (scaled, typically the scale factor = 1)
312 pscale(i) = pscale(i)*p0(i)
313c
314c ..deviatoric elastic predictor ..
315 signxx(i)=sigoxx(i)+p0(i)+g2(i)*(depsxx(i)-dav)
316 signyy(i)=sigoyy(i)+p0(i)+g2(i)*(depsyy(i)-dav)
317 signzz(i)=sigozz(i)+p0(i)+g2(i)*(depszz(i)-dav)
318 ENDDO
319C
320 signxy(1:nel)=sigoxy(1:nel)+g(1:nel) *depsxy(1:nel)
321 signyz(1:nel)=sigoyz(1:nel)+g(1:nel) *depsyz(1:nel)
322 signzx(1:nel)=sigozx(1:nel)+g(1:nel) *depszx(1:nel)
323 sigexx(1:nel) = signxx(1:nel)
324 sigeyy(1:nel) = signyy(1:nel)
325 sigezz(1:nel) = signzz(1:nel)
326 sigexy(1:nel) = signxy(1:nel)
327 sigeyz(1:nel) = signyz(1:nel)
328 sigezx(1:nel) = signzx(1:nel)
329 viscmax(1:nel) = zero
330 dpla1(1:nel) =zero
331 epstt(1:nel) = zero
332 fail(1:nel) = one
333C-------------------
334C STRAIN principal 1, 4 newton iterations (only when rupture)
335C-------------------
336 IF (ifail > 1) THEN
337 epsr1f = min(epsr1,epsf)
338 DO i=1,nel
339c IF(EPSR1F == EP30)CYCLE
340 dav = (epsxx(i)+epsyy(i)+epszz(i)) * third
341 e1 = epsxx(i) - dav
342 e2 = epsyy(i) - dav
343 e3 = epszz(i) - dav
344 e4 = half*epsxy(i)
345 e5 = half*epsyz(i)
346 e6 = half*epszx(i)
347C -y = (e1-x)(e2-x)(e3-x)
348C - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
349C + 2e4 e5 e6
350C e1 + e2 + e3 = 0 => terme en x^2 = 0
351C y = x^3 + c x + d
352c yp= 3 x^2 + c
353 e42 = e4*e4
354 e52 = e5*e5
355 e62 = e6*e6
356 c = - half *(e1*e1 + e2*e2 + e3*e3) - e42 - e52 - e62
357 epst = sqrt(-c*third)
358c 2*EPST est un majorant de la solution
359 epsr1dav = epsr1f - dav
360 IF(epst+epst < epsr1dav)cycle
361 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
362 & - two*e4*e5*e6
363 epst2 = epst * epst
364 y = (epst2 + c)* epst + d
365 IF(abs(y) > em8)THEN
366 epst = onep75 * epst
367 epst2 = epst * epst
368 y = (epst2 + c)* epst + d
369 yp = three*epst2 + c
370 epst = epst - y/yp
371c EPST est un majorant de la solution
372 IF(epst < epsr1dav)cycle
373 epst2 = epst * epst
374 y = (epst2 + c)* epst + d
375 yp = three*epst2 + c
376 epst = epst - y/yp
377 IF(epst < epsr1dav)cycle
378 epst2 = epst * epst
379 y = (epst2 + c)* epst + d
380 yp = three*epst2 + c
381 epst = epst - y/yp
382 IF(epst < epsr1dav)cycle
383 epst2 = epst * epst
384 y = (epst2 + c)* epst + d
385 yp = three*epst2 + c
386 epst = epst - y/yp
387 ENDIF
388 epst = epst + dav
389 epstt(i)= epst
390C-------------------
391C tension failure
392C-------------------
393 fail(i) = max(em20, min(one, (epsr2-epst)/(epsr2-epsr1) ))
394 dmg(i) = one - max(zero, min(one,(epsr2-epst)/(epsr2-epsr1)))
395 ENDDO
396 ENDIF
397C=======================================================================
398c Split code in 2 independent part :
399c VP = 1 dependent on plastic strain rate
400c VP = 0 dependent on total strain rate
401C=======================================================================
402 IF (vp == 0) THEN
403C=======================================================================
404 IF (nrate == 1) THEN ! only static curve => no strain rate interpolation
405 iad1(1:nel) = npf(ifunc(1)) / 2 + 1
406 ilen1(1:nel) = npf(ifunc(1)+1) / 2 - iad1(1:nel) - vartmp(1:nel,3)
407c
408 CALL vinter(tf,iad1,vartmp(1:nel,3),ilen1,nel,pla,dydx1,y1)
409c
410 yfac(1:nel,1) = uparam(6+nrate+1)*facyldi(1:nel)
411 fact(1:nel) = fail(1:nel) * pfac(1:nel) * yfac(1:nel,1)
412 h(1:nel) = dydx1(1:nel) * fact(1:nel)
413c
414 IF (fisokin == zero) THEN
415 yld(1:nel) = y1(1:nel) * fact(1:nel)
416 ELSE IF (fisokin == one) THEN
417 yld(1:nel) = tf(npf(ifunc(1))+1) * fact(1:nel)
418 ELSE IF (fisokin > zero) THEN
419 yld(1:nel) = ((one-fisokin)*y1(1:nel) + fisokin *tf(npf(ifunc(1))+1))*fact(1:nel)
420 END IF
421c-----------
422 ELSE ! strain rate interpolation
423c-----------
424 jj(1:nel) = 1
425 DO j=2,nrate-1
426#include "vectorize.inc"
427 DO i=1,nel
428 IF (epsp(i) >= uparam(6+j)) jj(i) = j
429 ENDDO
430 ENDDO
431C
432 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
433#include "vectorize.inc"
434 DO i=1,nel
435 epsp1 = max(uparam(6+jj(i)), em20)
436 epsp2 = uparam(7+jj(i))
437 rfac(i) = log(max(epsp(i),em20)/epsp1) / log(epsp2/epsp1)
438 ENDDO
439 ELSE ! linear interpolation of strain rate
440#include "vectorize.inc"
441 DO i=1,nel
442 epsp1 = uparam(6+jj(i))
443 epsp2 = uparam(7+jj(i))
444 rfac(i) = (epsp(i) - epsp1) / (epsp2 - epsp1)
445 ENDDO
446 END IF
447C
448#include "vectorize.inc"
449 DO i=1,nel
450 j1 = jj(i)
451 j2 = j1+1
452 ipos1(i) = vartmp(i,j1+2)
453 ipos2(i) = vartmp(i,j2+2)
454 yfac(i,1)= uparam(6+nrate+j1) * facyldi(i)
455 yfac(i,2)= uparam(6+nrate+j2) * facyldi(i)
456 ENDDO
457 iad1(1:nel) = npf(ifunc(jj(1:nel))) / 2 + 1
458 ilen1(1:nel) = npf(ifunc(jj(1:nel))+1) / 2 - iad1(1:nel) - ipos1(1:nel)
459 iad2(1:nel) = npf(ifunc(jj(1:nel)+1)) / 2 + 1
460 ilen2(1:nel) = npf(ifunc(jj(1:nel)+1)+1) / 2 - iad2(1:nel) - ipos2(1:nel)
461C
462 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
463 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
464c
465 DO i=1,nel
466 j1 = jj(i)
467 j2 = j1+1
468 vartmp(i,j1+2) = ipos1(i)
469 vartmp(i,j2+2) = ipos2(i)
470 ENDDO
471C
472 IF (fisokin == zero) THEN ! ecrouissage isotrope
473#include "vectorize.inc"
474 DO i=1,nel
475 j1 = jj(i)
476 j2 = j1+1
477 y1(i)= y1(i)*yfac(i,1)
478 y2(i)= y2(i)*yfac(i,2)
479 fac = rfac(i)
480 cc = fail(i)* pfac(i)
481 yld(i) = (y1(i) + fac*(y2(i)-y1(i))) * cc
482 dydx1(i)= dydx1(i)*yfac(i,1)
483 dydx2(i)= dydx2(i)*yfac(i,2)
484 h(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i))) * cc
485 ENDDO
486c
487 ELSEIF (fisokin == one) THEN ! ecrouissage cinematique
488c
489#include "vectorize.inc"
490 DO i=1,nel
491 j1 = jj(i)
492 j2 = j1+1
493 fac = rfac(i)
494 cc = fail(i) * pfac(i)
495 dydx1(i)=dydx1(i)*yfac(i,1)
496 dydx2(i)=dydx2(i)*yfac(i,2)
497 h(i) = cc*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
498 y1(i)= tf(npf(ifunc(j1))+1)
499 y2(i)= tf(npf(ifunc(j2))+1)
500 y1(i)= y1(i)*yfac(i,1)
501 y2(i)= y2(i)*yfac(i,2)
502 yld(i) = cc*(y1(i) + fac*(y2(i)-y1(i)))
503 ENDDO
504c
505 ELSE ! 0 < FISOKIN < 1 ecrouissage mixte
506c
507#include "vectorize.inc"
508 DO i=1,nel
509 j1 = jj(i)
510 j2 = j1+1
511 y1(i) = y1(i)*yfac(i,1)
512 y2(i) = y2(i)*yfac(i,2)
513 fac = rfac(i)
514 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
515 yld(i) = max(yld(i),em20)
516 dydx1(i)= dydx1(i)*yfac(i,1)
517 dydx2(i)= dydx2(i)*yfac(i,2)
518 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
519 h(i) = h(i) * pfac(i)
520
521 y1(i)=tf(npf(ifunc(j1))+1)
522 y2(i)=tf(npf(ifunc(j2))+1)
523 y1(i)=y1(i)*yfac(i,1)
524 y2(i)=y2(i)*yfac(i,2)
525 yld(i) = (one - fisokin) * yld(i)
526 . + fisokin * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
527 yld(i) = yld(i) * pfac(i)
528 ENDDO
529 ENDIF ! FISOKIN
530 END IF ! NRATE > 0
531c
532 IF (yldcheck == 1) THEN
533 DO i=1,nel
534 yld(i) = max(yld(i), em20)
535 END DO
536 END IF
537C-------------------
538C PROJECTION
539C-------------------
540 IF (ipla == 0) THEN
541 DO i=1,nel
542 vm =three*(half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
543 . + signxy(i)**2+signyz(i)**2+signzx(i)**2)
544 IF(vm > yld(i)*yld(i))THEN
545 vm = sqrt(vm)
546 r(i) = yld(i)/ max(vm,em20)
547 signxx(i)=signxx(i)*r(i)
548 signyy(i)=signyy(i)*r(i)
549 signzz(i)=signzz(i)*r(i)
550 signxy(i)=signxy(i)*r(i)
551 signyz(i)=signyz(i)*r(i)
552 signzx(i)=signzx(i)*r(i)
553 pla(i) = pla(i) + (one - r(i))*vm/max(g3(i)+h(i),em20)
554 dpla1(i) = (one - r(i))*vm/max(g3(i)+h(i),em20)
555 ENDIF
556 ENDDO
557c
558 ELSEIF(ipla == 2)THEN
559c
560 DO i=1,nel
561 vm = three*(half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
562 . + signxy(i)**2+signyz(i)**2+signzx(i)**2)
563 IF (vm > yld(i)*yld(i)) THEN
564 vm = sqrt(vm)
565 r(i) = yld(i)/ max(vm,em20)
566 signxx(i)=signxx(i)*r(i)
567 signyy(i)=signyy(i)*r(i)
568 signzz(i)=signzz(i)*r(i)
569 signxy(i)=signxy(i)*r(i)
570 signyz(i)=signyz(i)*r(i)
571 signzx(i)=signzx(i)*r(i)
572 pla(i)=pla(i) + (one-r(i))*vm/max(g3(i),em20)
573 dpla1(i) = (one-r(i))*vm/max(g3(i),em20)
574 ENDIF
575 ENDDO
576c
577 ELSEIF(ipla == 1)THEN
578c
579 IF (impl_s == 0) THEN
580 DO i=1,nel
581 vm = three*(half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
582 . + signxy(i)**2+signyz(i)**2+signzx(i)**2)
583 IF (vm > yld(i)*yld(i)) THEN
584 vm = sqrt(vm)
585 r(i) = yld(i)/ max(vm,em20)
586C plastic strain increment.
587 dpla =(one - r(i)) * vm/max(g3(i)+h(i),em20)
588C updated yield stress.
589 yld(i) = max(yld(i)+(one - fisokin)*dpla*h(i),zero)
590 r(i) = min(one,yld(i) / max(vm,em20))
591C
592c . . here the updated stresses are shifted stresses
593c (for kinematic or mixed hardening)
594 signxx(i)=signxx(i)*r(i)
595 signyy(i)=signyy(i)*r(i)
596 signzz(i)=signzz(i)*r(i)
597 signxy(i)=signxy(i)*r(i)
598 signyz(i)=signyz(i)*r(i)
599 signzx(i)=signzx(i)*r(i)
600 pla(i)=pla(i) + dpla
601 dpla1(i) = dpla
602 ENDIF
603 ENDDO
604C
605 ELSE ! -- IMPL_S > 0
606C
607c --- nonlinear hardening requires iterations in radial return --
608 ipflag = 0
609 jj(1:nel) = 1
610 IF (pfun > 0) ipflag = 1
611 DO i=1,nel
612C inside M36ITER_IMP, IPOS1=0
613 j1 = jj(i)
614 iad1(i) = npf(ifunc(j1)) / 2 + 1
615 ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i)
616 ENDDO
617 CALL m36iter_imp( signxx, signyy, signzz,
618 . signxy, signyz, signzx,
619 . pla, yld, g3,
620 . yfac, dpla1, h,
621 . tf, iad1,
622 . ilen1, nel,
623 . fisokin, vartmp,pla0,
624 . plam,y1, dydx1,
625 . ipflag, pfun, ipfun, iposp,
626 . nrate, npf, iadp, ilenp,
627 . pfac, pscale1, dfdp, fail,
628 . nvartmp ,
629 . sigbxx,sigbyy,sigbzz,sigbxy,sigbyz,sigbzx)
630C
631 ENDIF ! -- IMPL_S<=0
632c
633 ENDIF ! -- IPLA --
634C=======================================================================
635C=======================================================================
636 ELSE ! VP=1
637C=======================================================================
638C=======================================================================
639 IF (isigi == 0) THEN
640 IF(time == zero)THEN
641 DO i=1,nel
642 yfac(i,1) = uparam(6+nrate+1)*facyldi(i)
643 yld(i) = yfac(i,1)*finter(ifunc(1) ,zero,npf,tf,dydx1(i))
644 h(i) = dydx1(i)
645 yld(i) = yld(i)*max(zero,pfac(i))
646 h(i) = h(i) *max(zero,pfac(i))
647 uvar(i,3) = yld(i)
648 ENDDO
649 ENDIF
650 ENDIF
651C
652 DO i=1,nel
653 plaold(i) = pla(i)
654 plap(i) = uvar(i,2)
655 yld(i) = uvar(i,3)
656 jj(i) = 1
657 ENDDO
658C-------------------------
659C CRITERE DE VON MISES
660C-------------------------
661 DO i=1,nel
662 svm2(i) = three*(half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
663 . + signxy(i)**2+signyz(i)**2+signzx(i)**2)
664 ENDDO
665 nindx=0
666 DO i=1,nel
667 IF(svm2(i)>yld(i)*yld(i).AND.off(i)== one) THEN
668 nindx=nindx+1
669 index(nindx)=i
670 ENDIF
671 ENDDO
672C=======================================================================
673c PLASTICITY
674C=======================================================================
675 IF (nindx > 0 ) THEN
676c
677 IF (fisokin == zero) THEN !===No kinematic hardening
678 !-------------------
679 !Plastic increment :
680 !-------------------
681#include "vectorize.inc"
682 DO ii=1,nindx
683 i = index(ii)
684 svm_tab(i) = sqrt(svm2(i))
685 dpla_j(i) = uvar(i,1) + em09
686 !--------------
687 !update yield :
688 !--------------
689 jj(i) = 1
690 !Y1 and Y2 for Pl SR dependency
691 ENDDO
692 DO j=2,nrate-1
693 DO ii=1,nindx
694 i = index(ii)
695 IF (plap(i) >= uparam(6+j)) jj(i) = j
696 ENDDO
697 ENDDO
698 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
699#include "vectorize.inc"
700 DO ii=1,nindx
701 i = index(ii)
702 epsp1 = max(uparam(6+jj(i)), em20)
703 epsp2 = uparam(7+jj(i))
704 rfac(i) = log(max(plap(i),em20)/epsp1) / log(epsp2/epsp1)
705 ENDDO
706 ELSE ! linear interpolation of strain rate
707#include "vectorize.inc"
708 DO ii=1,nindx
709 i = index(ii)
710 epsp1 = uparam(6+jj(i))
711 epsp2 = uparam(7+jj(i))
712 rfac(i) = (plap(i) - epsp1) / (epsp2 - epsp1)
713 ENDDO
714 END IF ! linear interpolation of strain rate
715! ------------------------
716 nindex_vinter = 0
717#include "vectorize.inc"
718 DO ii=1,nindx
719 i = index(ii)
720 nindex_vinter = nindex_vinter + 1
721 index_vinter(nindex_vinter) = i
722 j1=jj(i)
723 j2=jj(i)+1
724 iposp(nindex_vinter) = vartmp(i,j1+2)
725 iadp(nindex_vinter) = npf(ifunc(jj(i)) ) / 2 + 1
726 ilenp(nindex_vinter) = npf(ifunc(jj(i)) + 1) / 2 - iadp(nindex_vinter) - iposp(nindex_vinter)
727 iposp2(nindex_vinter) = vartmp(i,j2+2)
728 iadp2(nindex_vinter) = npf(ifunc(jj(i)+1)) / 2 + 1
729 ilenp2(nindex_vinter) = npf(ifunc(jj(i)+1) + 1) / 2 - iadp2(nindex_vinter) - iposp2(nindex_vinter)
730 tab_loc(nindex_vinter) = pla(i)
731 yfac(i,1) = uparam(6+nrate+j1) * facyldi(i)
732 yfac(i,2) = uparam(6+nrate+j2) * facyldi(i)
733 ENDDO
734
735 CALL vinter2(tf,iadp,iposp,ilenp,nindex_vinter,tab_loc,dydx1_loc,y1_loc)
736 CALL vinter2(tf,iadp2,iposp2,ilenp2,nindex_vinter,tab_loc,dydx2_loc,y2_loc)
737
738 DO ii=1,nindex_vinter
739 i = index_vinter(ii)
740 j1 = jj(i)
741 j2 = j1+1
742 vartmp(i,j1+2) = iposp(ii)
743 vartmp(i,j2+2) = iposp2(ii)
744 ENDDO
745#include "vectorize.inc"
746 DO ii=1,nindex_vinter
747 i = index_vinter(ii)
748 y1(i) = y1_loc(ii)
749 y2(i) = y2_loc(ii)
750 dydx1(i)=dydx1_loc(ii)
751 dydx2(i)=dydx2_loc(ii)
752 ENDDO
753! ------------------------
754#include "vectorize.inc"
755 DO ii=1,nindx
756 i = index(ii)
757 y1(i) = y1(i) * yfac(i,1)
758 y2(i) = y2(i) * yfac(i,2)
759 fac = rfac(i)
760 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
761 yld(i) = max(yld(i),em20)
762 dydx1(i)=dydx1(i)*yfac(i,1)
763 dydx2(i)=dydx2(i)*yfac(i,2)
764 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
765 yld(i) = yld(i)*max(zero,pfac(i))
766 h(i) = h(i) *max(zero,pfac(i))
767 dfsr(i) = h(i) + max(zero,pfac(i))*fail(i)*(y2(i)-y1(i))/
768 . (uparam(7+jj(i))-uparam(6+jj(i))) /timestep
769 ENDDO
770c------------------------
771 DO k=1,niter
772 !update svm
773#include "vectorize.inc"
774 DO ii=1,nindx
775 i = index(ii)
776 dpla_i(i) = dpla_j(i)
777 pla(i) = plaold(i) + dpla_i(i)
778 plap(i) = dpla_i(i) / timestep
779 r(i) = yld(i)/(yld(i) +g3(i)*dpla_i(i))
780 svm = r(i)* svm_tab(i)
781
782 f = svm - yld(i) - g3(i) *dpla_i(i)
783 df = - g3(i) -dfsr(i)
784 df = sign(max(abs(df),em20),df)
785 IF(dpla_i(i) > zero) THEN
786 dpla_j(i)=max( em10 ,dpla_i(i)-f/df)
787 ELSE
788 dpla_j(i)=em10
789 ENDIF
790 !--------------
791 !update yield :
792 !--------------
793 pla(i) = plaold(i) + dpla_j(i)
794 plap(i) = dpla_j(i) / timestep
795 jj(i) = 1
796 ENDDO
797 DO j=2,nrate-1
798 DO ii=1,nindx
799 i = index(ii)
800 IF(plap(i) >= uparam(6+j)) jj(i) = j
801 ENDDO
802 ENDDO
803 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
804#include "vectorize.inc"
805 DO i=1,nel
806 epsp1 = max(uparam(6+jj(i)), em20)
807 epsp2 = uparam(7+jj(i))
808 rfac(i) = log(max(plap(i),em20)/epsp1) / log(epsp2/epsp1)
809 ENDDO
810 ELSE ! linear interpolation of strain rate
811#include "vectorize.inc"
812 DO ii=1,nindx
813 i = index(ii)
814 epsp1 = uparam(6+jj(i))
815 epsp2 = uparam(7+jj(i))
816 rfac(i) = (plap(i) - epsp1) / (epsp2 - epsp1)
817 ENDDO
818 END IF
819! ------------------------
820 nindex_vinter = 0
821#include "vectorize.inc"
822 DO ii=1,nindx
823 i = index(ii)
824 nindex_vinter = nindex_vinter + 1
825 index_vinter(nindex_vinter) = i
826 j1=jj(i)
827 j2=jj(i)+1
828 iposp(nindex_vinter) = vartmp(i,j1+2)
829 iadp(nindex_vinter) = npf(ifunc(jj(i)) ) / 2 + 1
830 ilenp(nindex_vinter) = npf(ifunc(jj(i)) + 1) / 2 - iadp(nindex_vinter) - iposp(nindex_vinter)
831 iposp2(nindex_vinter) = vartmp(i,j2+2)
832 iadp2(nindex_vinter) = npf(ifunc(jj(i)+1)) / 2 + 1
833 ilenp2(nindex_vinter) = npf(ifunc(jj(i)+1) + 1) / 2 - iadp2(nindex_vinter) - iposp2(nindex_vinter)
834 tab_loc(nindex_vinter) = pla(i)
835 yfac(i,1)=uparam(6+nrate+j1)*facyldi(i)
836 yfac(i,2)=uparam(6+nrate+j2)*facyldi(i)
837 ENDDO
838
839 CALL vinter2(tf,iadp,iposp,ilenp,nindex_vinter,tab_loc,dydx1_loc,y1_loc)
840 CALL vinter2(tf,iadp2,iposp2,ilenp2,nindex_vinter,tab_loc,dydx2_loc,y2_loc)
841
842 DO ii=1,nindex_vinter
843 i = index_vinter(ii)
844 j1 = jj(i)
845 j2 = j1+1
846 vartmp(i,j1+2) = iposp(ii)
847 vartmp(i,j2+2) = iposp2(ii)
848 ENDDO
849#include "vectorize.inc"
850 DO ii=1,nindex_vinter
851 i = index_vinter(ii)
852 y1(i) = y1_loc(ii)
853 y2(i) = y2_loc(ii)
854 dydx1(i)=dydx1_loc(ii)
855 dydx2(i)=dydx2_loc(ii)
856 ENDDO
857! ------------------------
858#include "vectorize.inc"
859 DO ii=1,nindx
860 i = index(ii)
861 !Y1 and Y2 for Pl SR dependency
862 ! PFAC => scale factor pression
863 y1(i) = yfac(i,1) * y1(i)
864 y2(i) = yfac(i,2) * y2(i)
865 fac = rfac(i)
866 yld(i)= fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
867 yld(i) = max(yld(i),em20)
868 yld(i) = yld(i)*max(zero,pfac(i))
869 dydx1(i)=dydx1(i)*yfac(i,1)
870 dydx2(i)=dydx2(i)*yfac(i,2)
871 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
872 h(i) = h(i) *max(zero,pfac(i))
873 dfsr(i)= h(i)+ max(zero,pfac(i))*fail(i)*(y2(i)-y1(i))
874 . /(uparam(7+jj(i))-uparam(6+jj(i))) /timestep
875 ENDDO
876 ENDDO !iter
877! ------------------------
878#include "vectorize.inc"
879 DO ii=1,nindx
880 i = index(ii)
881 pla(i) = plaold(i) + dpla_i(i)
882 plap(i) = dpla_i(i) / timestep
883 signxx(i) = signxx(i)*r(i)
884 signyy(i) = signyy(i)*r(i)
885 signzz(i) = signzz(i)*r(i)
886 signxy(i) = signxy(i)*r(i)
887 signyz(i) = signyz(i)*r(i)
888 signzx(i) = signzx(i)*r(i)
889 dpla1(i) = dpla_i(i)
890 uvar(i,1) = dpla_i(i)
891 uvar(i,2) = plap(i)
892 uvar(i,3) = yld(i)
893 ENDDO !NINDX
894c------------------------------------------------
895 ELSEIF (fisokin == one ) THEN ! PURE kinematic hardening
896c------------------------------------------------
897 !Plastic increment :
898 !-------------------
899#include "vectorize.inc"
900 DO ii=1,nindx
901 i = index(ii)
902 svm_tab(i) = sqrt(svm2(i))
903 dpla_j(i) = uvar(i,1) + em09
904 !--------------
905 !update yield :
906 !--------------
907 jj(i) = 1
908 ENDDO
909 DO j=2,nrate-1
910#include "vectorize.inc"
911 DO ii=1,nindx
912 i = index(ii)
913 IF (plap(i) >= uparam(6+j)) jj(i) = j
914 ENDDO
915 ENDDO
916c
917 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
918#include "vectorize.inc"
919 DO ii=1,nindx
920 i = index(ii)
921 epsp1 = max(uparam(6+jj(i)), em20)
922 epsp2 = uparam(7+jj(i))
923 rfac(i) = log(max(plap(i),em20)/epsp1) / log(epsp2/epsp1)
924 ENDDO
925 ELSE ! linear interpolation of strain rate
926#include "vectorize.inc"
927 DO ii=1,nindx
928 i = index(ii)
929 epsp1 = uparam(6+jj(i))
930 epsp2 = uparam(7+jj(i))
931 rfac(i) = (plap(i) - epsp1) / (epsp2 - epsp1)
932 ENDDO
933 END IF
934! ------------------------
935 nindex_vinter = 0
936#include "vectorize.inc"
937 DO ii=1,nindx
938 i = index(ii)
939 nindex_vinter = nindex_vinter + 1
940 index_vinter(nindex_vinter) = i
941 j1=jj(i)
942 j2=jj(i)+1
943 iposp(nindex_vinter) = vartmp(i,j1+2)
944 iadp(nindex_vinter) = npf(ifunc(jj(i)) ) / 2 + 1
945 ilenp(nindex_vinter) = npf(ifunc(jj(i)) + 1) / 2 - iadp(nindex_vinter) - iposp(nindex_vinter)
946 iposp2(nindex_vinter) = vartmp(i,j2+2)
947 iadp2(nindex_vinter) = npf(ifunc(jj(i)+1)) / 2 + 1
948 ilenp2(nindex_vinter) = npf(ifunc(jj(i)+1) + 1) / 2 - iadp2(nindex_vinter) - iposp2(nindex_vinter)
949 tab_loc(nindex_vinter) = pla(i)
950 yfac(i,1)=uparam(6+nrate+j1)*facyldi(i)
951 yfac(i,2)=uparam(6+nrate+j2)*facyldi(i)
952 ENDDO
953
954 CALL vinter2(tf,iadp,iposp,ilenp,nindex_vinter,tab_loc,dydx1_loc,y1_loc)
955 CALL vinter2(tf,iadp2,iposp2,ilenp2,nindex_vinter,tab_loc,dydx2_loc,y2_loc)
956
957 DO ii=1,nindex_vinter
958 i = index_vinter(ii)
959 j1 = jj(i)
960 j2 = jj(i)+1
961 vartmp(i,j1+2) = iposp(ii)
962 vartmp(i,j2+2) = iposp2(ii)
963 ENDDO
964#include "vectorize.inc"
965 DO ii=1,nindex_vinter
966 i = index_vinter(ii)
967 y1(i) = y1_loc(ii)
968 y2(i) = y2_loc(ii)
969 dydx1(i)=dydx1_loc(ii)
970 dydx2(i)=dydx2_loc(ii)
971 ENDDO
972! ------------------------
973#include "vectorize.inc"
974 DO ii=1,nindx
975 i = index(ii)
976 !Y1 and Y2 for Pl SR dependency
977 j1 = jj(i)
978 j2 = jj(i)+1
979 y1(i) = yfac(i,1)*y1(i)
980 y2(i) = yfac(i,2)*y2(i)
981 fac = rfac(i)
982 yld(i)= fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
983
984 dydx1(i)=dydx1(i)*yfac(i,1)
985 dydx2(i)=dydx2(i)*yfac(i,2)
986 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
987 h(i) = h(i) *max(zero,pfac(i))
988 hi(i) = h(i)*(one-fisokin)
989
990 coef(i) = (y2(i)-y1(i))/(uparam(6+j2)-uparam(6+j1)) /timestep
991 dfsr(i) = h(i)+ coef(i)
992 dfsr(i) = dfsr(i) * (one-fisokin)*max(zero,pfac(i))
993 !Hardening
994 y1(i)=tf(npf(ifunc(j1))+1)
995 y2(i)=tf(npf(ifunc(j2))+1)
996 y1(i)= y1(i)*yfac(i,1)
997 y2(i)= y2(i)*yfac(i,2)
998 coef(i) = max(zero,pfac(i))*fisokin * fail(i)*(y2(i)-y1(i))
999 . /(uparam(6+j2)-uparam(6+j1)) /timestep
1000 yld(i) = (one-fisokin) * yld(i) +
1001 . fisokin * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
1002 yld(i) = max(yld(i),em20)
1003 yld(i) = yld(i) *max(zero,pfac(i))
1004 ENDDO
1005 !----------------------
1006 !-- NEWTON ITERATIONS :
1007 !----------------------
1008 DO k=1,niter
1009#include "vectorize.inc"
1010 DO ii=1,nindx
1011 i = index(ii)
1012 !update svm
1013 dpla_i(i) = dpla_j(i)
1014 pla(i) = plaold(i) + dpla_i(i)
1015 plap(i) = dpla_i(i) / timestep
1016 r(i) = yld(i)/(yld(i) +(g3(i)+fisokin*h(i))*dpla_i(i))
1017 svm = r(i)* svm_tab(i)
1018 f = svm - yld(i) - (g3(i)+fisokin*h(i)) *dpla_i(i)
1019 df =-(g3(i)+fisokin*h(i)) - (dfsr(i) + coef(i) )
1020 df = sign(max(abs(df),em20),df)
1021 IF(dpla_i(i) > zero) THEN
1022 dpla_j(i)=max( em10 ,dpla_i(i)-f/df)
1023 ELSE
1024 dpla_j(i)=em10
1025 ENDIF
1026 !--------------
1027 !update yield :
1028 !--------------
1029 pla(i) = plaold(i) + dpla_j(i)
1030 plap(i) = dpla_j(i) / timestep
1031 jj(i) = 1
1032 ENDDO
1033 DO j=2,nrate-1
1034#include "vectorize.inc"
1035 DO ii=1,nindx
1036 i = index(ii)
1037 IF(plap(i) >= uparam(6+j)) jj(i) = j
1038 ENDDO
1039 ENDDO
1040c
1041 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
1042#include "vectorize.inc"
1043 DO ii=1,nindx
1044 i = index(ii)
1045 epsp1 = max(uparam(6+jj(i)), em20)
1046 epsp2 = uparam(7+jj(i))
1047 rfac(i) = log(max(plap(i),em20)/epsp1) / log(epsp2/epsp1)
1048 ENDDO
1049 ELSE ! linear interpolation of strain rate
1050#include "vectorize.inc"
1051 DO ii=1,nindx
1052 i = index(ii)
1053 epsp1 = uparam(6+jj(i))
1054 epsp2 = uparam(7+jj(i))
1055 rfac(i) = (plap(i) - epsp1) / (epsp2 - epsp1)
1056 ENDDO
1057 END IF
1058! ------------------------
1059 nindex_vinter = 0
1060#include "vectorize.inc"
1061 DO ii=1,nindx
1062 i = index(ii)
1063 nindex_vinter = nindex_vinter + 1
1064 index_vinter(nindex_vinter) = i
1065 j1=jj(i)
1066 j2=jj(i)+1
1067 iposp(nindex_vinter) = vartmp(i,j1+2)
1068 iadp(nindex_vinter) = npf(ifunc(jj(i)) ) / 2 + 1
1069 ilenp(nindex_vinter) = npf(ifunc(jj(i)) + 1) / 2 - iadp(nindex_vinter) - iposp(nindex_vinter)
1070 iposp2(nindex_vinter) = vartmp(i,j2+2)
1071 iadp2(nindex_vinter) = npf(ifunc(jj(i)+1)) / 2 + 1
1072 ilenp2(nindex_vinter) = npf(ifunc(jj(i)+1) + 1) / 2 - iadp2(nindex_vinter) - iposp2(nindex_vinter)
1073 tab_loc(nindex_vinter) = pla(i)
1074 yfac(i,1)=uparam(6+nrate+j1)*facyldi(i)
1075 yfac(i,2)=uparam(6+nrate+j2)*facyldi(i)
1076 ENDDO
1077
1078 CALL vinter2(tf,iadp,iposp,ilenp,nindex_vinter,tab_loc,dydx1_loc,y1_loc)
1079 CALL vinter2(tf,iadp2,iposp2,ilenp2,nindex_vinter,tab_loc,dydx2_loc,y2_loc)
1080
1081 DO ii=1,nindex_vinter
1082 i = index_vinter(ii)
1083 j1 = jj(i)
1084 j2 = j1+1
1085 vartmp(i,j1+2) = iposp(ii)
1086 vartmp(i,j2+2) = iposp2(ii)
1087 ENDDO
1088#include "vectorize.inc"
1089 DO ii=1,nindex_vinter
1090 i = index_vinter(ii)
1091 y1(i) = y1_loc(ii)
1092 y2(i) = y2_loc(ii)
1093 dydx1(i)=dydx1_loc(ii)
1094 dydx2(i)=dydx2_loc(ii)
1095 ENDDO
1096! ------------------------
1097#include "vectorize.inc"
1098 DO ii=1,nindx
1099 i = index(ii)
1100 !Y1 and Y2 for Pl SR dependency
1101 j1 = jj(i)
1102 j2 = jj(i)+1
1103 y1(i) = yfac(i,1) * y1(i)
1104 y2(i) = yfac(i,2) * y2(i)
1105 fac = rfac(i)
1106 yld(i)= fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
1107 dydx1(i)=dydx1(i)*yfac(i,1)
1108 dydx2(i)=dydx2(i)*yfac(i,2)
1109 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
1110 h(i) = h(i) *max(zero,pfac(i))
1111 hi(i) = h(i)*(one-fisokin)
1112 dfsr(i) = hi(i)+ max(zero,pfac(i))* (one-fisokin)*(y2(i)-y1(i))
1113 . / (uparam(7+jj(i))-uparam(6+jj(i))) /timestep
1114 y1(i)=tf(npf(ifunc(j1))+1)
1115 y2(i)=tf(npf(ifunc(j2))+1)
1116 y1(i)=y1(i)*yfac(i,1)
1117 y2(i)=y2(i)*yfac(i,2)
1118 coef(i) = max(zero,pfac(i))*fisokin * fail(i)*(y2(i)-y1(i))
1119 . / (uparam(6+j2)-uparam(6+j1)) /timestep
1120 yld(i) = (one-fisokin) * yld(i)
1121 . + fisokin * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
1122 yld(i) = max(yld(i),em20)
1123 yld(i) = yld(i) * max(zero,pfac(i))
1124 ENDDO
1125 ENDDO !iter
1126! ------------------------
1127#include "vectorize.inc"
1128 DO ii=1,nindx
1129 i = index(ii)
1130 pla(i) = plaold(i) + dpla_i(i)
1131 plap(i) = dpla_i(i) / timestep
1132 signxx(i) = signxx(i)*r(i)
1133 signyy(i) = signyy(i)*r(i)
1134 signzz(i) = signzz(i)*r(i)
1135 signxy(i) = signxy(i)*r(i)
1136 signyz(i) = signyz(i)*r(i)
1137 signzx(i) = signzx(i)*r(i)
1138 dpla1(i) = dpla_i(i)
1139 uvar(i,1) = dpla_i(i)
1140 uvar(i,2) = plap(i)
1141 uvar(i,3) = yld(i)
1142 ENDDO !NINDX
1143c------------------------------------------------
1144 ELSE ! (1 > FISOKIN > 0) : Mixed ISO-KINE HARDENING
1145c------------------------------------------------
1146#include "vectorize.inc"
1147 DO ii=1,nindx
1148 i = index(ii)
1149 svm_tab(i) = sqrt(svm2(i))
1150 dpla_j(i) = uvar(i,1) + em09
1151 !DPLA_J(I) = (SVM-YLD(I))/(G3(I)+H(I))!converge moins bien avec estimation radial
1152 !--------------
1153 !update yield :
1154 !--------------
1155 jj(i) = 1
1156 ENDDO
1157 DO j=2,nrate-1
1158#include "vectorize.inc"
1159 DO ii=1,nindx
1160 i = index(ii)
1161 IF(plap(i) >= uparam(6+j)) jj(i) = j
1162 ENDDO
1163 ENDDO
1164c
1165 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
1166#include "vectorize.inc"
1167 DO ii=1,nindx
1168 i = index(ii)
1169 epsp1 = max(uparam(6+jj(i)), em20)
1170 epsp2 = uparam(7+jj(i))
1171 rfac(i) = log(max(plap(i),em20)/epsp1) / log(epsp2/epsp1)
1172 ENDDO
1173 ELSE ! linear interpolation of strain rate
1174#include "vectorize.inc"
1175 DO ii=1,nindx
1176 i = index(ii)
1177 epsp1 = uparam(6+jj(i))
1178 epsp2 = uparam(7+jj(i))
1179 rfac(i) = (plap(i) - epsp1) / (epsp2 - epsp1)
1180 ENDDO
1181 END IF
1182! ------------------------
1183 nindex_vinter = 0
1184#include "vectorize.inc"
1185 DO ii=1,nindx
1186 i = index(ii)
1187 nindex_vinter = nindex_vinter + 1
1188 index_vinter(nindex_vinter) = i
1189 j1=jj(i)
1190 j2=jj(i)+1
1191 iposp(nindex_vinter) = vartmp(i,j1+2)
1192 iadp(nindex_vinter) = npf(ifunc(jj(i)) ) / 2 + 1
1193 ilenp(nindex_vinter) = npf(ifunc(jj(i)) + 1) / 2 - iadp(nindex_vinter) - iposp(nindex_vinter)
1194 iposp2(nindex_vinter) = vartmp(i,j2+2)
1195 iadp2(nindex_vinter) = npf(ifunc(jj(i)+1)) / 2 + 1
1196 ilenp2(nindex_vinter) = npf(ifunc(jj(i)+1) + 1) / 2 - iadp2(nindex_vinter) - iposp2(nindex_vinter)
1197 tab_loc(nindex_vinter) = pla(i)
1198 yfac(i,1)=uparam(6+nrate+j1)*facyldi(i)
1199 yfac(i,2)=uparam(6+nrate+j2)*facyldi(i)
1200 ENDDO
1201
1202 CALL vinter2(tf,iadp,iposp,ilenp,nindex_vinter,tab_loc,dydx1_loc,y1_loc)
1203 CALL vinter2(tf,iadp2,iposp2,ilenp2,nindex_vinter,tab_loc,dydx2_loc,y2_loc)
1204
1205 DO ii=1,nindex_vinter
1206 i = index_vinter(ii)
1207 j1 = jj(i)
1208 j2 = j1+1
1209 vartmp(i,j1+2) = iposp(ii)
1210 vartmp(i,j2+2) = iposp2(ii)
1211 ENDDO
1212#include "vectorize.inc"
1213 DO ii=1,nindex_vinter
1214 i = index_vinter(ii)
1215 y1(i) = y1_loc(ii)
1216 y2(i) = y2_loc(ii)
1217 dydx1(i)=dydx1_loc(ii)
1218 dydx2(i)=dydx2_loc(ii)
1219 ENDDO
1220! ------------------------
1221#include "vectorize.inc"
1222 DO ii=1,nindx
1223 i = index(ii)
1224 !Y1 and Y2 for Pl SR dependency
1225 j1 = jj(i)
1226 j2 = jj(i)+1
1227 y1(i) = yfac(i,1) * y1(i)
1228 y2(i) = yfac(i,2) * y2(i)
1229 fac = rfac(i)
1230 yld(i)= fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
1231
1232 dydx1(i)=dydx1(i)*yfac(i,1)
1233 dydx2(i)=dydx2(i)*yfac(i,2)
1234 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
1235 h(i) = h(i) *max(zero,pfac(i))
1236 hi(i) = h(i)*(one-fisokin)
1237
1238 coef(i) = (y2(i)-y1(i))/(uparam(6+j2)-uparam(6+j1)) /timestep
1239 dfsr(i) = h(i)+ coef(i)
1240 dfsr(i) = dfsr(i) * (one-fisokin)*max(zero,pfac(i))
1241 !Hardening
1242 y1(i)=tf(npf(ifunc(j1))+1)
1243 y2(i)=tf(npf(ifunc(j2))+1)
1244 y1(i)= y1(i)*yfac(i,1)
1245 y2(i)= y2(i)*yfac(i,2)
1246 coef(i) = max(zero,pfac(i))*fisokin * fail(i)*(y2(i)-y1(i))
1247 . /(uparam(6+j2)-uparam(6+j1)) /timestep
1248 yld(i) = (one-fisokin) * yld(i) +
1249 . fisokin * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
1250 yld(i) = max(yld(i),em20)
1251 yld(i) = yld(i)*max(zero,pfac(i))
1252 ENDDO
1253 !----------------------
1254 !-- NEWTON ITERATIONS :
1255 !----------------------
1256 DO k=1,niter
1257#include "vectorize.inc"
1258 DO ii=1,nindx
1259 i = index(ii)
1260 !update svm
1261 dpla_i(i) = dpla_j(i)
1262 pla(i) = plaold(i) + dpla_i(i)
1263 plap(i) = dpla_i(i) / timestep
1264 r(i) = yld(i)/(yld(i) +(g3(i)+fisokin*h(i))*dpla_i(i))
1265 svm = r(i)* svm_tab(i)
1266 f = svm - yld(i) - (g3(i)+fisokin*h(i)) *dpla_i(i)
1267 df =-(g3(i)+fisokin*h(i)) - (dfsr(i) + coef(i))
1268 df = sign(max(abs(df),em20),df)
1269 IF(dpla_i(i) > zero) THEN
1270 dpla_j(i)=max( em10 ,dpla_i(i)-f/df)
1271 ELSE
1272 dpla_j(i)=em10
1273 ENDIF
1274 !--------------
1275 !update yield :
1276 !--------------
1277 pla(i) = plaold(i) + dpla_j(i)
1278 plap(i) = dpla_j(i) / timestep
1279 jj(i) = 1
1280 ENDDO
1281 DO j=2,nrate-1
1282#include "vectorize.inc"
1283 DO ii=1,nindx
1284 i = index(ii)
1285 IF(plap(i) >= uparam(6+j)) jj(i) = j
1286 ENDDO
1287 ENDDO
1288c
1289 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
1290#include "vectorize.inc"
1291 DO ii=1,nindx
1292 i = index(ii)
1293 epsp1 = max(uparam(6+jj(i)), em20)
1294 epsp2 = uparam(7+jj(i))
1295 rfac(i) = log(max(plap(i),em20)/epsp1) / log(epsp2/epsp1)
1296 ENDDO
1297 ELSE ! linear interpolation of strain rate
1298#include "vectorize.inc"
1299 DO ii=1,nindx
1300 i = index(ii)
1301 epsp1 = uparam(6+jj(i))
1302 epsp2 = uparam(7+jj(i))
1303 rfac(i) = (plap(i) - epsp1) / (epsp2 - epsp1)
1304 ENDDO
1305 END IF
1306! ------------------------
1307 nindex_vinter = 0
1308#include "vectorize.inc"
1309 DO ii=1,nindx
1310 i = index(ii)
1311 nindex_vinter = nindex_vinter + 1
1312 index_vinter(nindex_vinter) = i
1313 j1=jj(i)
1314 j2=jj(i)+1
1315 iposp(nindex_vinter) = vartmp(i,j1+2)
1316 iadp(nindex_vinter) = npf(ifunc(jj(i)) ) / 2 + 1
1317 ilenp(nindex_vinter) = npf(ifunc(jj(i)) + 1) / 2 - iadp(nindex_vinter) - iposp(nindex_vinter)
1318 iposp2(nindex_vinter) = vartmp(i,j2+2)
1319 iadp2(nindex_vinter) = npf(ifunc(jj(i)+1)) / 2 + 1
1320 ilenp2(nindex_vinter) = npf(ifunc(jj(i)+1) + 1) / 2 - iadp2(nindex_vinter) - iposp2(nindex_vinter)
1321 tab_loc(nindex_vinter) = pla(i)
1322 yfac(i,1)=uparam(6+nrate+j1)*facyldi(i)
1323 yfac(i,2)=uparam(6+nrate+j2)*facyldi(i)
1324 ENDDO
1325
1326 CALL vinter2(tf,iadp,iposp,ilenp,nindex_vinter,tab_loc,dydx1_loc,y1_loc)
1327 CALL vinter2(tf,iadp2,iposp2,ilenp2,nindex_vinter,tab_loc,dydx2_loc,y2_loc)
1328
1329 DO ii=1,nindex_vinter
1330 i = index_vinter(ii)
1331 j1 = jj(i)
1332 j2 = j1+1
1333 vartmp(i,j1+2) = iposp(ii)
1334 vartmp(i,j2+2) = iposp2(ii)
1335 ENDDO
1336#include "vectorize.inc"
1337 DO ii=1,nindex_vinter
1338 i = index_vinter(ii)
1339 y1(i) = y1_loc(ii)
1340 y2(i) = y2_loc(ii)
1341 dydx1(i)=dydx1_loc(ii)
1342 dydx2(i)=dydx2_loc(ii)
1343 ENDDO
1344! ------------------------
1345#include "vectorize.inc"
1346 DO ii=1,nindx
1347 i = index(ii)
1348 !Y1 and Y2 for Pl SR dependency
1349 j1 = jj(i)
1350 j2 = jj(i)+1
1351 y1(i) = yfac(i,1) * y1(i)
1352 y2(i) = yfac(i,2) * y2(i)
1353 fac = rfac(i)
1354 yld(i)= fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
1355 dydx1(i)=dydx1(i)*yfac(i,1)
1356 dydx2(i)=dydx2(i)*yfac(i,2)
1357 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
1358 h(i) = h(i) *max(zero,pfac(i))
1359 hi(i)= h(i)*(one-fisokin)
1360 dfsr(i) = hi(i)+ max(zero,pfac(i))* (one-fisokin)*(y2(i)-y1(i))
1361 . /(uparam(7+jj(i))-uparam(6+jj(i))) /timestep
1362 y1(i)=tf(npf(ifunc(j1))+1)
1363 y2(i)=tf(npf(ifunc(j2))+1)
1364 y1(i)=y1(i)*yfac(i,1)
1365 y2(i)=y2(i)*yfac(i,2)
1366 coef(i) = max(zero,pfac(i))*fisokin * fail(i)*(y2(i)-y1(i))
1367 . /(uparam(6+j2)-uparam(6+j1)) /timestep
1368 yld(i) = (one-fisokin) * yld(i) +
1369 . fisokin * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
1370 yld(i) = max(yld(i),em20)
1371 yld(i) = yld(i) *max(zero,pfac(i))
1372 ENDDO
1373 ENDDO !iter
1374#include "vectorize.inc"
1375 DO ii=1,nindx
1376 i = index(ii)
1377 pla(i) = plaold(i) + dpla_i(i)
1378 plap(i) = dpla_i(i) / timestep
1379 signxx(i) = signxx(i)*r(i)
1380 signyy(i) = signyy(i)*r(i)
1381 signzz(i) = signzz(i)*r(i)
1382 signxy(i) = signxy(i)*r(i)
1383 signyz(i) = signyz(i)*r(i)
1384 signzx(i) = signzx(i)*r(i)
1385 dpla1(i) = dpla_i(i)
1386 uvar(i,1) = dpla_i(i)
1387 uvar(i,2) = plap(i)
1388 uvar(i,3) = yld(i)
1389 ENDDO
1390 ENDIF ! FISOKIN
1391 ENDIF ! NINDX > 0 (plasticity)
1392c------------------------------------------
1393 ENDIF ! VP flag
1394C=======================================================================
1395 IF (ipla /= 1 .OR. impl_s <= 0) THEN
1396C------------------------------------------
1397 IF (fisokin == one ) THEN ! ECROUISSAGE CINE
1398#include "vectorize.inc"
1399 DO i=1,nel
1400 dsxx = sigexx(i) - signxx(i)
1401 dsyy = sigeyy(i) - signyy(i)
1402 dszz = sigezz(i) - signzz(i)
1403 dsxy = sigexy(i) - signxy(i)
1404 dsyz = sigeyz(i) - signyz(i)
1405 dszx = sigezx(i) - signzx(i)
1406
1407 hkin = two_third*fisokin*h(i)
1408 alpha = hkin/(g2(i)+hkin)
1409 sigpxx = alpha*dsxx
1410 sigpyy = alpha*dsyy
1411 sigpzz = alpha*dszz
1412 sigpxy = alpha*dsxy
1413 sigpyz = alpha*dsyz
1414 sigpzx = alpha*dszx
1415
1416c ..updates back stresses
1417 sigbxx(i) = sigbxx(i) + sigpxx
1418 sigbyy(i) = sigbyy(i) + sigpyy
1419 sigbzz(i) = sigbzz(i) + sigpzz
1420 sigbxy(i) = sigbxy(i) + sigpxy
1421 sigbyz(i) = sigbyz(i) + sigpyz
1422 sigbzx(i) = sigbzx(i) + sigpzx
1423
1424c ..gets stresses from shifted stresses and back stresses
1425 signxx(i) = signxx(i) + sigbxx(i)
1426 signyy(i) = signyy(i) + sigbyy(i)
1427 signzz(i) = signzz(i) + sigbzz(i)
1428 signxy(i) = signxy(i) + sigbxy(i)
1429 signyz(i) = signyz(i) + sigbyz(i)
1430 signzx(i) = signzx(i) + sigbzx(i)
1431 ENDDO
1432
1433 ELSEIF (fisokin > zero) THEN ! ecrouissage mixte
1434#include "vectorize.inc"
1435 DO i=1,nel
1436 dsxx = sigexx(i) - signxx(i)
1437 dsyy = sigeyy(i) - signyy(i)
1438 dszz = sigezz(i) - signzz(i)
1439 dsxy = sigexy(i) - signxy(i)
1440 dsyz = sigeyz(i) - signyz(i)
1441 dszx = sigezx(i) - signzx(i)
1442
1443 hkin = two_third*fisokin*h(i)
1444 alpha = hkin/(g2(i)+hkin)
1445 sigpxx = alpha*dsxx
1446 sigpyy = alpha*dsyy
1447 sigpzz = alpha*dszz
1448 sigpxy = alpha*dsxy
1449 sigpyz = alpha*dsyz
1450 sigpzx = alpha*dszx
1451
1452c ..updates back stresses
1453 sigbxx(i) = sigbxx(i) + sigpxx
1454 sigbyy(i) = sigbyy(i) + sigpyy
1455 sigbzz(i) = sigbzz(i) + sigpzz
1456 sigbxy(i) = sigbxy(i) + sigpxy
1457 sigbyz(i) = sigbyz(i) + sigpyz
1458 sigbzx(i) = sigbzx(i) + sigpzx
1459
1460c ..gets stresses from shifted stresses and back stresses
1461 signxx(i) = signxx(i) + sigbxx(i)
1462 signyy(i) = signyy(i) + sigbyy(i)
1463 signzz(i) = signzz(i) + sigbzz(i)
1464 signxy(i) = signxy(i) + sigbxy(i)
1465 signyz(i) = signyz(i) + sigbyz(i)
1466 signzx(i) = signzx(i) + sigbzx(i)
1467 ENDDO
1468 ENDIF ! ..FISOKIN
1469 END IF !(IPLA/=1.OR.IMPL_S<=0 ) THEN
1470
1471 IF (ieos == 0) THEN ! add pressure to the deviatoric stress
1472 DO i=1,nel
1473! P = BULK * (RHO/RHO0-ONE)
1474 p = bulk(i) * amu(i)
1475 signxx(i) = signxx(i) - p
1476 signyy(i) = signyy(i) - p
1477 signzz(i) = signzz(i) - p
1478 ENDDO
1479 ELSE
1480 ! if EOS is used, material law calculates only deviatoric stress tensor
1481 ! sound speed depends on pressure derivative over volume change
1482 ! calculated in EOS
1483 DO i = 1, nel
1484 soundsp(i) = sqrt((dpdm(i) + four*g(i)/three)/rho0(i))
1485 ENDDO
1486 END IF
1487c
1488c-------------------------
1489 IF (impl_s > 0) THEN
1490c ..saves the shifted elastic predictors, plastic strain multiplier,
1491c the Kirchhoff's modulus and the plastic hardening modulus
1492c for the sake of computation of elasto-plastic stiffness matrix
1493c
1494 DO i=1,nel
1495 IF(dpla1(i) > 0) etse(i)= h(i)/g2(i)
1496 ENDDO
1497 DO i = 1,nel
1498 IF (dpla1(i) > zero) THEN
1499
1500c ...... Von Mises stress at the elastic predictor (point B)
1501c ...... SIGEXX, etc. are deviatoric stresses
1502 vm =half*(sigexx(i)**2+sigeyy(i)**2+sigezz(i)**2)
1503 . +sigexy(i)**2+sigeyz(i)**2+sigezx(i)**2
1504 vm_1 =one/sqrt(three*vm)
1505 g3h = g3(i)+h(i)
1506 r(i) = max(zero,one-g3h*dpla1(i)*vm_1)
1507c ...... NORM_1 normalizes deviatoric stresses, includes consistent
1508c ...... stiffness matrix parameter beta, von Mises at B, and two_pmi
1509 norm_1=g3(i)*vm_1*sqrt(r(i)/g3h)
1510c ...... Deviatoric stresses "normalized"
1511 signor(i,1)=sigexx(i)*norm_1
1512 signor(i,2)=sigeyy(i)*norm_1
1513 signor(i,3)=sigezz(i)*norm_1
1514 signor(i,4)=sigexy(i)*norm_1
1515 signor(i,5)=sigeyz(i)*norm_1
1516 signor(i,6)=sigezx(i)*norm_1
1517
1518c ...... Parameter alpha of consistent matrix
1519 al_imp(i)= one - g3(i)*dpla1(i)*vm_1
1520 ELSE
1521 al_imp(i)=one
1522 ENDIF
1523 ENDDO
1524 ENDIF ! IMPL_S > 0
1525c-------------------------
1526 DO i=1,nel
1527 IF (off(i) < em01) off(i) = zero
1528 IF (off(i) < one) off(i) = off(i)*four_over_5
1529 ENDDO
1530c
1531 IF (ifail > 0) THEN
1532 nindx = 0
1533 IF (ifail == 2) THEN
1534 IF (inloc > 0) THEN
1535 DO i=1,nel
1536 IF (epsmax < ep20) dmg(i) = max(dmg(i),planl(i)/epsmax)
1537 IF ((planl(i) > epsmax .OR. epstt(i) > epsf) .AND. off(i) == one) THEN
1538 off(i)= four_over_5
1539 nindx = nindx+1
1540 indx(nindx) = i
1541 ENDIF
1542 ENDDO
1543 ELSE
1544 DO i=1,nel
1545 IF (epsmax < ep20) dmg(i) = max(dmg(i),pla(i)/epsmax)
1546 IF ((pla(i) > epsmax .OR. epstt(i) > epsf) .AND. off(i) == one) THEN
1547 off(i)= four_over_5
1548 nindx = nindx+1
1549 indx(nindx) = i
1550 ENDIF
1551 ENDDO
1552 ENDIF
1553 ELSE
1554 IF (inloc > 0) THEN
1555 DO i=1,nel
1556 IF (epsmax < ep20) dmg(i) = planl(i)/epsmax
1557 IF (planl(i) > epsmax .AND. off(i) == one) THEN
1558 off(i)= four_over_5
1559 nindx = nindx+1
1560 indx(nindx) = i
1561 ENDIF
1562 ENDDO
1563 ELSE
1564 DO i=1,nel
1565 IF (epsmax < ep20) dmg(i) = pla(i)/epsmax
1566 IF (pla(i) > epsmax .AND. off(i) == one) THEN
1567 off(i)= four_over_5
1568 nindx = nindx+1
1569 indx(nindx) = i
1570 ENDIF
1571 ENDDO
1572 ENDIF
1573 END IF
1574 IF (nindx > 0 .AND. imconv == 1) THEN
1575 DO j=1,nindx
1576#include "lockon.inc"
1577 WRITE(iout, 1000) ngl(indx(j))
1578 WRITE(istdo,1100) ngl(indx(j)),tt
1579#include "lockoff.inc"
1580 ENDDO
1581 ENDIF
1582 ENDIF ! IFAIL
1583C-----------
1584 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
1585 1100 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10, ' AT TIME :',g11.4)
1586C-----------
1587 RETURN
1588 END
alpha
#define alpha
Definition eval.h:35
m36iter_imp
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)
Definition m36iter_imp.F:46
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps36
subroutine sigeps36(nel, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, viscmax, uvar, off, ngl, ieos, ipm, mat, epsp, ipla, yld, pla, dpla1, etse, al_imp, signor, amu, dpdm, facyldi, nvartmp, vartmp, dmg, inloc, planl, sigbxx, sigbyy, sigbzz, sigbxy, sigbyz, sigbzx)
Definition sigeps36.F:46
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72