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