OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps36c.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!|| sigeps36c ../engine/source/materials/mat/mat036/sigeps36c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| vinter ../engine/source/tools/curve/vinter.F
29!|| vinter2 ../engine/source/tools/curve/vinter.F
30!||====================================================================
31 SUBROUTINE sigeps36c(
32 1 NEL ,NUVAR ,NVARTMP ,MFUNC ,
33 2 KFUNC ,NPF ,IPLAS ,
34 2 TF ,TIMESTEP,UPARAM ,RHO0 ,
35 3 THKLY ,ISRATE ,ASRATE ,EPSD_PG ,EPSD ,
36 4 EPSPXX ,EPSPYY ,EPSPXY ,
37 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 6 EPSXX ,EPSYY ,EPSXY ,
39 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 9 SOUNDSP,VISCMAX,THK ,PLA ,UVAR ,
42 A VARTMP ,OFF ,IPM ,IMAT ,
43 B ETSE ,GS ,YLD ,
44 C DPLA_I ,GAMA_IMP,SIGNOR ,SHF ,HARDM ,
45 D FACYLDI,INLOC ,DPLANL ,DMG ,PLANL ,
46 E L_PLANL,SIGBXX,SIGBYY,SIGBXY ,LOFF)
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.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 NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
63C IFUNC | NFUNC | I | R | FUNCTION INDEX
64C NPF | * | I | R | FUNCTION ARRAY
65C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
66C IPLAS | * | I | R | GEOMETRICAL FLAGS
67C TF | * | F | R | FUNCTION ARRAY
68C---------+---------+---+---+--------------------------------------------
69C TIME | 1 | F | R | CURRENT TIME
70C TIMESTEP| 1 | F | R | CURRENT TIME STEP
71C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
72C RHO0 | NEL | F | R | INITIAL DENSITY
73C THKLY | NEL | F | R | LAYER THICKNESS
74C EPSPXX | NEL | F | R | STRAIN RATE XX
75C EPSPYY | NEL | F | R | STRAIN RATE YY
76C ... | | | |
77C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
78C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
79C ... | | | |
80C EPSXX | NEL | F | R | STRAIN XX
81C EPSYY | NEL | F | R | STRAIN YY
82C ... | | | |
83C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
84C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
85C ... | | | |
86C---------+---------+---+---+--------------------------------------------
87C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
88C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
89C ... | | | |
90C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
91C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
92C---------+---------+---+---+--------------------------------------------
93C THK | NEL | F |R/W| THICKNESS
94C PLA | NEL | F |R/W| PLASTIC STRAIN
95C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
96C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
97C---------+---------+---+---+--------------------------------------------
98C C o m m o n B l o c k s
99C-----------------------------------------------
100#include "param_c.inc"
101#include "impl1_c.inc"
102C-----------------------------------------------
103C I N P U T A r g u m e n t s
104C-----------------------------------------------
105 INTEGER NEL, NUVAR, NVARTMP,IMAT,ISRATE,INLOC
106 INTEGER IPLAS(*),IPM(NPROPMI,*)
107 INTEGER, INTENT(IN) :: L_PLANL
108 my_real :: TIMESTEP,ASRATE
109 my_real :: UPARAM(*),RHO0(NEL),
110 . THKLY(NEL),PLA(NEL),FACYLDI(NEL),
111 . EPSPXX(NEL),EPSPYY(NEL),EPSPXY(NEL),
112 . DEPSXX(NEL),DEPSYY(NEL),DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
113 . EPSXX(NEL) ,EPSYY(NEL),EPSXY(NEL),
114 . SIGOXX(NEL),SIGOYY(NEL),SIGOXY(NEL),SIGOYZ(NEL),SIGOZX(NEL),
115 . GS(NEL),SHF(NEL),HARDM(NEL),DPLANL(NEL),DMG(NEL),
116 . PLANL(NEL*L_PLANL)
117 my_real, DIMENSION(NEL), INTENT(IN) :: LOFF
118 my_real, DIMENSION(NEL), INTENT(IN) :: EPSD_PG ! global element strain rate
119 my_real, DIMENSION(NEL), INTENT(INOUT) :: EPSD ! lbuf%epsd
120C-----------------------------------------------
121C O U T P U T A r g u m e n t s
122C-----------------------------------------------
123 my_real
124 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
125 . soundsp(nel),viscmax(nel),etse(nel),
126 . dpla_i(nel),gama_imp(nel),signor(mvsiz,5)
127C-----------------------------------------------
128C I N P U T O U T P U T A r g u m e n t s
129C-----------------------------------------------
130 my_real
131 . uvar(nel,nuvar), off(nel),thk(nel),yld(nel)
132 INTEGER, INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
133 my_real, DIMENSION(NEL), INTENT(INOUT) :: sigbxx,sigbyy,sigbxy
134C-----------------------------------------------
135C VARIABLES FOR FUNCTION INTERPOLATION
136C-----------------------------------------------
137 INTEGER :: NPF(*), MFUNC, KFUNC(MFUNC)
138 my_real :: tf(*)
139C-----------------------------------------------
140C L o c a l V a r i a b l e s
141C-----------------------------------------------
142 INTEGER :: I,J,N,II,J1,J2,NINDX,VP,IFAIL,OPTE,NINDEX_PLA,
143 . NITER,NRATE,IPFUN,NFUNC,PFUN,IFUNCE,FUN1,FUN2,YLDCHECK,
144 . ISMOOTH
145 INTEGER IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
146 . IAD2(MVSIZ),IPOS2(MVSIZ),ILEN2(MVSIZ),
147 . jj(mvsiz),index(mvsiz),iposp(mvsiz),ipospe(mvsiz),
148 . ifunc(100),iadp(mvsiz),ilenp(mvsiz)
149 INTEGER, DIMENSION(MVSIZ) :: INDEX_PLA
150 my_real :: A,P2,R,FAC,DEZZ,SVM,
151 . HKIN,DTINV,AAA,SOUNDSP1,ALPHA,CE,EINF,EPSP1,EPSP2,
152 . f,df,q2,yld_i,sigpxx,sigpyy,sigpxy,
153 . e1,a11,a21,g1,g31,nu11,nu21,nu_mnu,u_mnu,t_pnu,nux,nu3,
154 . epsmax,epsr1,epsr2,epsf,fisokin,s1,s2,s3,yld0
155 my_real :: fact(nel),e(mvsiz),a1(mvsiz),a2(mvsiz),g(mvsiz),
156 . dydx1(mvsiz),dydx2(mvsiz),rfac(mvsiz),
157 . y1(mvsiz),y2(mvsiz),dr(mvsiz),yfac(mvsiz,2),escale(mvsiz),
158 . aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
159 . pp(mvsiz),qq(mvsiz),fail(mvsiz),h(mvsiz),hk(mvsiz),
160 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz),
161 . svm2(mvsiz),yld2(mvsiz),hi(mvsiz),g3(mvsiz),epst(mvsiz),
162 . pfac(mvsiz),p0(mvsiz),dfdp(mvsiz),pscale(mvsiz),
163 . dydxe(mvsiz),plap(mvsiz)
164 my_real, DIMENSION(MVSIZ) :: tab_loc
165c------------------
166 DATA niter/3/
167C=======================================================================
168C USER VARIABLES INITIALIZATION
169C-----------------------------------------------
170 nfunc = ipm(10,imat)
171 DO j=1,nfunc
172 ifunc(j) = ipm(10+j,imat)
173 ENDDO
174 ipfun = ifunc(nfunc-1)
175C
176 e1 = uparam(2)
177 a11 = uparam(3)
178 a21 = uparam(4)
179 g1 = uparam(5)
180 nux = uparam(6)
181 nrate = nint(uparam(1))
182 epsmax = uparam(2*nrate + 7)
183 epsr1 = uparam(2*nrate + 8)
184 epsr2 = uparam(2*nrate + 9)
185 g31 = uparam(2*nrate + 11)
186 fisokin= uparam(2*nrate + 14)
187 epsf = uparam(2*nrate + 15)
188 pfun = nint(uparam(2*nrate + 16))
189 soundsp1 = uparam(2*nrate + 18)
190 nu_mnu = uparam(2*nrate + 19) ! NU / (1 - NU)
191 t_pnu = uparam(2*nrate + 20) ! 3 / (1 + NU)
192 u_mnu = uparam(2*nrate + 21) ! 1 / (1 - NU)
193 opte = uparam(2*nrate + 23)
194 einf = uparam(2*nrate + 24)
195 ce = uparam(2*nrate + 25)
196 vp = nint(uparam(2*nrate + 26)) ! flag for plastic strain dependency
197 ifail = nint(uparam(2*nrate + 27)) ! flag for failure
198 yldcheck = nint(uparam(2*nrate + 28))
199 ismooth = nint(uparam(2*nrate + 29)) ! strain rate interpolation flag
200c-----------------------------
201 nindex_pla = 0
202 viscmax(1:nel) = zero
203 etse(1:nel) = one
204 e(1:nel) = e1
205 a1(1:nel) = a11
206 a2(1:nel) = a21
207 g(1:nel) = g1
208 g3(1:nel) = g31
209 soundsp(1:nel) = soundsp1
210c-----------------------
211 IF (opte == 1) THEN ! Variable Young modulus defined by tabulated function
212 ifunce = uparam(2*nrate + 22)
213 IF (ifunce > 0) THEN
214 DO i=1,nel
215 IF(pla(i) > zero)THEN
216 nindex_pla = nindex_pla + 1
217 index_pla(nindex_pla) = i
218 ipospe(nindex_pla) = vartmp(i,1)
219 iadp(nindex_pla) = npf(kfunc(ifunce)) / 2 + 1
220 ilenp(nindex_pla) = npf(kfunc(ifunce)+1) / 2 - iadp(nindex_pla) - ipospe(nindex_pla)
221 tab_loc(nindex_pla) = pla(i)
222 ENDIF
223 ENDDO
224 CALL vinter2(tf,iadp,ipospe,ilenp,nindex_pla,tab_loc,dydxe,escale)
225 vartmp(index_pla(1:nindex_pla),1) = ipospe(1:nindex_pla)
226 ENDIF
227
228#include "vectorize.inc"
229 DO ii=1,nindex_pla
230 i = index_pla(ii)
231 e(i) = escale(ii)* e1
232 a1(i) = e(i)/(one - nux*nux)
233 a2(i) = nux*a1(i)
234 g(i) = half*e(i)/(one+nux)
235 g3(i) = three*g(i)
236 gs(i) = g(i) * shf(i)
237 soundsp(i) = sqrt(e(i)/(one - nux*nux)/rho0(i))
238 ENDDO
239 ELSEIF (ce /= zero) THEN ! Variable Young modulus defined analytically
240 DO i=1,nel
241 IF (pla(i) > zero) THEN
242 e(i) = e1-(e1-einf)*(one-exp(-ce*pla(i)))
243 a1(i) = e(i)/(one - nux*nux)
244 a2(i) = nux*a1(i)
245 g(i) = half*e(i)/(one+nux)
246 g3(i) = three*g(i)
247 gs(i) = g(i) * shf(i)
248 soundsp(i) = sqrt(e(i)/(one - nux*nux)/rho0(i))
249 ENDIF
250 ENDDO
251 ENDIF ! OPTE
252c-------------------
253c damage factor based on equivalent strain
254c-------------------
255 IF (ifail == 2) THEN
256 epst(1:nel) = half*( epsxx(1:nel)+epsyy(1:nel)
257 . + sqrt( (epsxx(1:nel)-epsyy(1:nel))*(epsxx(1:nel)-epsyy(1:nel))
258 . + epsxy(1:nel)*epsxy(1:nel) ) )
259 fail(1:nel) = max(em20,min(one,(epsr2-epst(1:nel))/(epsr2-epsr1)))
260 dmg(1:nel) = one - max(zero,min(one,(epsr2-epst(1:nel))/(epsr2-epsr1)))
261 ELSE
262 fail(1:nel) = one
263 ENDIF
264C=======================================================================
265c Split code in 2 independent part :
266c VP = 1 dependent on plastic strain rate
267c VP = 0 dependent on total strain rate
268C=======================================================================
269 IF (vp == 0) THEN
270c------------------------------------------
271 sigoxx(1:nel) = sigoxx(1:nel) - sigbxx(1:nel)
272 sigoyy(1:nel) = sigoyy(1:nel) - sigbyy(1:nel)
273 sigoxy(1:nel) = sigoxy(1:nel) - sigbxy(1:nel)
274C
275 p0(1:nel) = -(sigoxx(1:nel) + sigoyy(1:nel))*third
276 signxx(1:nel) = sigoxx(1:nel)+a1(1:nel)*depsxx(1:nel)+a2(1:nel)*depsyy(1:nel)
277 signyy(1:nel) = sigoyy(1:nel)+a2(1:nel)*depsxx(1:nel)+a1(1:nel)*depsyy(1:nel)
278 signxy(1:nel) = sigoxy(1:nel)+g(1:nel) *depsxy(1:nel)
279 signyz(1:nel) = sigoyz(1:nel)+gs(1:nel) *depsyz(1:nel)
280 signzx(1:nel) = sigozx(1:nel)+gs(1:nel) *depszx(1:nel)
281 sigexx(1:nel) = signxx(1:nel)
282 sigeyy(1:nel) = signyy(1:nel)
283 sigexy(1:nel) = signxy(1:nel)
284C-------------------
285C STRAIN RATE
286C-------------------
287 IF (israte == 0) THEN
288 DO i=1,nel
289 epsd(i) = half*( abs(epspxx(i)+epspyy(i))
290 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
291 . + epspxy(i)*epspxy(i) ) )
292 ENDDO
293 ELSE
294 epsd(1:nel) = asrate*epsd_pg(1:nel) + (one-asrate)*epsd(1:nel)
295 ENDIF
296
297c-----------------------------------------------
298c PRESSURE DEPENDENT YIELD FUNCTION FACTOR
299c-----------------------------------------------
300 IF (pfun > 0) THEN
301 pscale(1:nel) = uparam(17 + 2*nrate )*p0(1:nel)
302 DO i=1,nel
303 iposp(i) = vartmp(i,2)
304 iadp(i) = npf(ipfun) / 2 + 1
305 ilenp(i) = npf(ipfun) / 2 - iadp(i) - iposp(i)
306 ENDDO
307 CALL vinter2(tf,iadp,iposp,ilenp,nel,pscale,dfdp,pfac)
308 pfac(1:nel) = max(zero, pfac(1:nel))
309 vartmp(1:nel,2) = iposp(1:nel)
310 ELSE
311 pfac(1:nel) = one
312 ENDIF
313C-------------------
314C CRITERE
315C-------------------
316 IF (nrate == 1) THEN ! only static curve => no strain rate interpolation
317 ipos1(1:nel) = vartmp(1:nel,3)
318 iad1(1:nel) = npf(ifunc(1)) / 2 + 1
319 ilen1(1:nel) = npf(ifunc(1)+1) / 2 - iad1(1:nel) - ipos1(1:nel)
320c
321 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
322c
323 yfac(1:nel,1) = uparam(6+nrate+1)*facyldi(1:nel)
324 vartmp(1:nel,3) = ipos1(1:nel)
325 fact(1:nel) = fail(1:nel) * pfac(1:nel) * yfac(1:nel,1)
326 h(1:nel) = dydx1(1:nel) * fact(1:nel)
327c
328 IF (fisokin == zero) THEN
329 yld(1:nel) = y1(1:nel) * fact(1:nel)
330 ELSE IF (fisokin == one) THEN
331 yld0 = tf(npf(ifunc(1))+1)
332 yld(1:nel) = yld0 * fact(1:nel)
333 ELSE IF (fisokin > zero) THEN
334 yld0 = tf(npf(ifunc(1))+1)
335 yld(1:nel) = ((one-fisokin)*y1(1:nel) + fisokin *yld0)*fact(1:nel)
336 END IF
337c
338 ELSE ! strain rate interpolation
339 jj(1:nel) = 1
340 DO j=2,nrate-1
341 DO i=1,nel
342 IF (epsd(i) >= uparam(6+j)) jj(i) = j
343 ENDDO
344 ENDDO
345 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
346#include "vectorize.inc"
347 DO i=1,nel
348 epsp1 = max(uparam(6+jj(i)), em20)
349 epsp2 = uparam(7+jj(i))
350 rfac(i) = log(max(epsd(i),em20)/epsp1) / log(epsp2/epsp1)
351 yfac(i,1) = uparam(6+nrate+jj(i)) * facyldi(i)
352 yfac(i,2) = uparam(7+nrate+jj(i)) * facyldi(i)
353 ENDDO
354 ELSE ! linear interpolation of strain rate
355#include "vectorize.inc"
356 DO i=1,nel
357 epsp1 = uparam(6+jj(i))
358 epsp2 = uparam(7+jj(i))
359 rfac(i) = (epsd(i) - epsp1) / (epsp2 - epsp1)
360 yfac(i,1) = uparam(6+nrate+jj(i)) * facyldi(i)
361 yfac(i,2) = uparam(7+nrate+jj(i)) * facyldi(i)
362 ENDDO
363 END IF
364
365#include "vectorize.inc"
366 DO i=1,nel
367 j1 = jj(i)
368 j2 = j1+1
369 fun1 = ifunc(j1)
370 fun2 = ifunc(j2)
371 ipos1(i) = vartmp(i,2+j1)
372 ipos2(i) = vartmp(i,2+j2)
373 iad1(i) = npf(fun1) / 2 + 1
374 ilen1(i) = npf(fun1+1) / 2 - iad1(i)-ipos1(i)
375 iad2(i) = npf(fun2) / 2 + 1
376 ilen2(i) = npf(fun2+1) / 2 - iad2(i)-ipos2(i)
377 END DO
378C
379 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
380 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
381c
382 IF (fisokin == zero) THEN
383#include "vectorize.inc"
384 DO i=1,nel
385 j1 = jj(i)
386 j2 = j1+1
387 y1(i) = y1(i)*yfac(i,1)
388 y2(i) = y2(i)*yfac(i,2)
389 fac = rfac(i)
390 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
391 yld(i) = max(yld(i),em20)
392 dydx1(i)=dydx1(i)*yfac(i,1)
393 dydx2(i)=dydx2(i)*yfac(i,2)
394 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
395 yld(i) = yld(i)*max(zero,pfac(i))
396 h(i) = h(i)*max(zero,pfac(i))
397 ENDDO
398 DO i=1,nel
399 j1 = jj(i)
400 j2 = j1+1
401 vartmp(i,2+j1) = ipos1(i)
402 vartmp(i,2+j2) = ipos2(i)
403 ENDDO
404 ELSEIF (fisokin == one) THEN
405#include "vectorize.inc"
406 DO i=1,nel
407 j1 = jj(i)
408 j2 = j1+1
409 fun1 = ifunc(j1)
410 fun2 = ifunc(j2)
411 fac = rfac(i)
412 dydx1(i)=dydx1(i)*yfac(i,1)
413 dydx2(i)=dydx2(i)*yfac(i,2)
414 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
415C ECROUISSAGE CINEMATIQUE
416 y1(i)=tf(npf(fun1)+1)
417 y2(i)=tf(npf(fun2)+1)
418 y1(i)=y1(i)*yfac(i,1)
419 y2(i)=y2(i)*yfac(i,2)
420 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
421 yld(i) = yld(i)*max(zero,pfac(i))
422 h(i) = h(i)*max(zero,pfac(i))
423 ENDDO
424 DO i=1,nel
425 j1 = jj(i)
426 j2 = j1+1
427 vartmp(i,2+j1) = ipos1(i)
428 vartmp(i,2+j2) = ipos2(i)
429 ENDDO
430 ELSE ! Mixed hardening
431#include "vectorize.inc"
432 DO i=1,nel
433 j1 = jj(i)
434 j2 = j1+1
435 fun1 = ifunc(j1)
436 fun2 = ifunc(j2)
437 y1(i)=y1(i)*yfac(i,1)
438 y2(i)=y2(i)*yfac(i,2)
439 fac = rfac(i)
440 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
441 yld(i) = max(yld(i),em20)
442 dydx1(i)=dydx1(i)*yfac(i,1)
443 dydx2(i)=dydx2(i)*yfac(i,2)
444 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
445 y1(i)=tf(npf(fun1)+1)
446 y2(i)=tf(npf(fun2)+1)
447 y1(i)=y1(i)*yfac(i,1)
448 y2(i)=y2(i)*yfac(i,2)
449 yld(i) = (one - fisokin) * yld(i) +
450 . fisokin * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
451 yld(i) = yld(i)*max(zero,pfac(i))
452 h(i) = h(i)*max(zero,pfac(i))
453 ENDDO
454 DO i=1,nel
455 j1 = jj(i)
456 j2 = j1+1
457 vartmp(i,2+j1) = ipos1(i)
458 vartmp(i,2+j2) = ipos2(i)
459 ENDDO
460 ENDIF ! FISOKIN
461 END IF ! NRATE > 0
462c
463 IF (yldcheck == 1) THEN
464 DO i=1,nel
465 yld(i) = max(yld(i), em20)
466 END DO
467 END IF
468c------------------------------
469C PROJECTION
470C-------------------
471 IF (iplas(1) == 0) THEN
472c projection radiale
473 nu3 = one-nu_mnu
474 DO i=1,nel
475 svm2(i)= signxx(i)*signxx(i)
476 . + signyy(i)*signyy(i)
477 . - signxx(i)*signyy(i)
478 . + three*signxy(i)*signxy(i)
479 IF (svm2(i)>yld(i)*yld(i)) THEN
480 svm =sqrt(svm2(i))
481 r = yld(i)/svm
482 signxx(i)=signxx(i)*r
483 signyy(i)=signyy(i)*r
484 signxy(i)=signxy(i)*r
485 dpla_i(i) = off(i)*svm*(one-r)/(g3(i)+h(i))
486 pla(i) = pla(i) + dpla_i(i)
487 IF (inloc == 0) THEN
488 IF(yld(i) .NE. 0) THEN
489 !YLD may be 0 when SIGNXX+SIGNYY == 0, that will cause floating point
490 !exception in debug mode
491 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
492 ELSE
493 dezz = zero
494 ENDIF
495 dezz=-(depsxx(i)+depsyy(i))*nu_mnu-nu3*dezz
496 thk(i) = thk(i) + dezz*thkly(i)*off(i)
497 ENDIF
498 etse(i)= h(i)/(h(i)+e(i)) !
499 ENDIF
500 ENDDO
501C
502 ELSEIF (iplas(1) == 1) THEN ! implicit scheme, 3 Newton iterations
503C-------------------------
504C CRITERE DE VON MISES
505C-------------------------
506 DO i=1,nel
507 h(i) = max(zero,h(i))
508 s1=signxx(i)+signyy(i)
509 s2=signxx(i)-signyy(i)
510 s3=signxy(i)
511 aa(i)= fourth*s1*s1
512 bb(i)= three_over_4*s2*s2 + three*s3*s3
513 svm2(i)= aa(i)+bb(i)
514 ENDDO
515 IF (inloc == 0) THEN
516 DO i=1,nel
517 dezz = -(depsxx(i)+depsyy(i))*nu_mnu
518 thk(i) = thk(i) + dezz*thkly(i)*off(i)
519 ENDDO
520 ENDIF
521C-------------------------
522C GATHER PLASTIC FLOW
523C-------------------------
524 nindx=0
525 DO i=1,nel
526 IF (svm2(i) > yld(i)*yld(i) .AND. off(i) == one) THEN
527 nindx=nindx+1
528 index(nindx)=i
529 ENDIF
530 ENDDO
531C
532 IF (nindx > 0) THEN
533C---------------------------
534C DEP EN CONTRAINTE PLANE
535C---------------------------
536#include "vectorize.inc"
537 DO j=1,nindx
538 i=index(j)
539 svm=sqrt(svm2(i))
540 dpla_j(i)=(svm-yld(i))/(g3(i)+h(i))
541 etse(i)= h(i)/(h(i)+e(i))
542 hi(i) = h(i)*(one-fisokin)
543 hk(i) = two_third*h(i)*fisokin
544 ENDDO
545C
546 nu3 = one-nu_mnu
547 DO n=1,niter
548#include "vectorize.inc"
549 DO j=1,nindx
550 i=index(j)
551 dpla_i(i)=dpla_j(i)
552 yld_i = yld(i) + hi(i)*dpla_i(i)
553 dr(i) = half*e(i)*dpla_i(i)/yld_i
554 aaa = three*hk(i)/e(i)
555 nu11 = u_mnu+aaa
556 nu21 = t_pnu+aaa
557 pp(i) = one/(one+dr(i)*nu11)
558 qq(i) = one/(one+dr(i)*nu21)
559 p2 = pp(i)*pp(i)
560 q2 = qq(i)*qq(i)
561 f = aa(i)*p2+bb(i)*q2-yld_i*yld_i
562 df =-(aa(i)*nu11*p2*pp(i)+bb(i)*nu21*q2*qq(i))
563 . *(e(i)-two*dr(i)*hi(i))/yld_i
564 . -two*hi(i)*yld_i
565 df = sign(max(abs(df),em20),df)
566 IF(dpla_i(i) > zero) THEN
567 dpla_j(i)=max(zero,dpla_i(i)-f/df)
568 ELSE
569 dpla_j(i)=zero
570 ENDIF
571 ENDDO
572 ENDDO
573C------------------------------------------
574C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
575C------------------------------------------
576#include "vectorize.inc"
577 DO j=1,nindx
578 i=index(j)
579 pla(i) = pla(i) + dpla_i(i)
580 s1=(signxx(i)+signyy(i))*pp(i)
581 s2=(signxx(i)-signyy(i))*qq(i)
582 signxx(i)=half*(s1+s2)
583 signyy(i)=half*(s1-s2)
584 signxy(i)=signxy(i)*qq(i)
585 IF (inloc == 0) THEN
586 dezz = - nu3*dr(i)*s1/e(i)
587 thk(i) = thk(i) + dezz*thkly(i)*off(i)
588 ENDIF
589C----- for kin hardening
590 yld(i) =yld(i)+hi(i)*dpla_i(i)
591 ENDDO
592 ENDIF
593c-------------------------------------------
594 ELSEIF (iplas(1) == 2) THEN
595c-------------------------------------------
596C CRITERE DE VON MISES
597C-------------------------
598 DO i=1,nel
599 h(i) = max(zero,h(i))
600 svm2(i) = signxx(i)*signxx(i) + signyy(i)*signyy(i)
601 . - signxx(i)*signyy(i) + three*signxy(i)*signxy(i)
602 IF (inloc == 0) THEN
603 dezz = -(depsxx(i)+depsyy(i))*nu_mnu
604 thk(i) = thk(i) + dezz*thkly(i)*off(i)
605 ENDIF
606 ENDDO
607C-------------------------
608C GATHER PLASTIC FLOW
609C-------------------------
610 nindx=0
611 DO i=1,nel
612 yld2(i)=yld(i)*yld(i)
613 IF (svm2(i) > yld2(i) .AND. off(i) == one) THEN
614 nindx=nindx+1
615 index(nindx)=i
616 ENDIF
617 ENDDO
618C
619 IF (nindx > 0) THEN
620C-------------
621C PROJ NORMALE AU CRITERE AVEC CALCUL APPROCHE DE LA NORMALE + RETOUR RADIAL
622C-------------
623 nu3 = one-nu_mnu
624#include "vectorize.inc"
625 DO j=1,nindx
626 i=index(j)
627 a = (svm2(i)-yld2(i))
628 . / (five*svm2(i)+three*(-signxx(i)*signyy(i)+signxy(i)*signxy(i)))
629 s1= (one-two*a)*signxx(i)+ a*signyy(i)
630 s2= a*signxx(i)+(one-two*a)*signyy(i)
631 s3=(one-three*a)*signxy(i)
632 signxx(i)=s1
633 signyy(i)=s2
634 signxy(i)=s3
635 svm=sqrt(svm2(i))
636 dpla_i(i) = off(i)*(svm-yld(i))/(g3(i)+h(i))
637C
638 hk(i) = h(i)*(one-fisokin)
639 yld(i) =yld(i)+hk(i)*dpla_i(i)
640 END DO
641C
642 DO j=1,nindx
643 i=index(j)
644 svm = sqrt( signxx(i)*signxx(i)
645 . +signyy(i)*signyy(i)
646 . -signxx(i)*signyy(i)
647 . +three*signxy(i)*signxy(i))
648 r = min(one,yld(i)/max(em20,svm))
649 signxx(i)=signxx(i)*r
650 signyy(i)=signyy(i)*r
651 signxy(i)=signxy(i)*r
652 pla(i) = pla(i) + dpla_i(i)
653 IF (inloc == 0) THEN
654 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
655 dezz = -nu3*dezz
656 thk(i) = thk(i) + dezz*thkly(i)*off(i)
657 ENDIF
658 etse(i)= h(i)/(h(i)+e(i))
659 END DO
660 END IF
661C=======================================================================
662 ENDIF ! IPLAS
663C=======================================================================
664 ELSE ! VP = 1 : plastic strain rate dependency
665C=======================================================================
666c
667 sigoxx(1:nel) = sigoxx(1:nel) - sigbxx(1:nel)
668 sigoyy(1:nel) = sigoyy(1:nel) - sigbyy(1:nel)
669 sigoxy(1:nel) = sigoxy(1:nel) - sigbxy(1:nel)
670C
671 p0(1:nel) = -(sigoxx(1:nel) + sigoyy(1:nel))*third
672 signxx(1:nel) = sigoxx(1:nel)+a1(1:nel)*depsxx(1:nel)+a2(1:nel)*depsyy(1:nel)
673 signyy(1:nel) = sigoyy(1:nel)+a2(1:nel)*depsxx(1:nel)+a1(1:nel)*depsyy(1:nel)
674 signxy(1:nel) = sigoxy(1:nel)+g(1:nel) *depsxy(1:nel)
675 signyz(1:nel) = sigoyz(1:nel)+gs(1:nel) *depsyz(1:nel)
676 signzx(1:nel) = sigozx(1:nel)+gs(1:nel) *depszx(1:nel)
677 sigexx(1:nel) = signxx(1:nel)
678 sigeyy(1:nel) = signyy(1:nel)
679 sigexy(1:nel) = signxy(1:nel)
680c-----------------------------------------------
681c PRESSURE DEPENDENT YIELD FUNCTION FACTOR
682c-----------------------------------------------
683 IF (pfun > 0) THEN
684 pscale(1:nel) = uparam(17 + 2*nrate )*p0(1:nel)
685 DO i=1,nel
686 iposp(i) = vartmp(i,2)
687 iadp(i) = npf(ipfun) / 2 + 1
688 ilenp(i) = npf(ipfun) / 2 - iadp(i) - iposp(i)
689 ENDDO
690 CALL vinter2(tf,iadp,iposp,ilenp,nel,pscale,dfdp,pfac)
691 pfac(1:nel) = max(zero, pfac(1:nel))
692 vartmp(1:nel,2) = iposp(1:nel)
693 ELSE
694 pfac(1:nel) = one
695 ENDIF
696C-------------------
697C PLASTIC STRAIN RATE
698C-------------------
699 plap(1:nel) = uvar(1:nel,2)
700C-------------------
701C CRITERE
702C-------------------
703 jj(1:nel) = 1
704 DO j=2,nrate-1
705 DO i=1,nel
706 IF (plap(i) >= uparam(6+j)) jj(i) = j
707 ENDDO
708 ENDDO
709c
710 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
711#include "vectorize.inc"
712 DO i=1,nel
713 epsp1 = max(uparam(6+jj(i)), em20)
714 epsp2 = uparam(7+jj(i))
715 rfac(i) = log(max(plap(i),em20)/epsp1) / log(epsp2/epsp1)
716 ENDDO
717 ELSE ! linear interpolation of strain rate
718#include "vectorize.inc"
719 DO i=1,nel
720 epsp1 = uparam(6+jj(i))
721 epsp2 = uparam(7+jj(i))
722 rfac(i) = (plap(i) - epsp1) / (epsp2 - epsp1)
723 ENDDO
724 END IF
725#include "vectorize.inc"
726 DO i=1,nel
727 yfac(i,1) = uparam(6+nrate+jj(i)) * facyldi(i)
728 yfac(i,2) = uparam(7+nrate+jj(i)) * facyldi(i)
729 ENDDO
730
731#include"vectorize.inc"
732 DO i=1,nel
733 j1 = jj(i)
734 j2 = j1+1
735 fun1 = ifunc(j1)
736 fun2 = ifunc(j2)
737 ipos1(i) = vartmp(i,2+j1)
738 ipos2(i) = vartmp(i,2+j2)
739 iad1(i) = npf(fun1) / 2 + 1
740 ilen1(i) = npf(fun1+1) / 2 - iad1(i)-ipos1(i)
741 iad2(i) = npf(fun2) / 2 + 1
742 ilen2(i) = npf(fun2+1) / 2 - iad2(i)-ipos2(i)
743 END DO
744C
745 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
746 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
747c-------------------------------
748 IF (fisokin == zero) THEN
749c-------------------------------
750#include "vectorize.inc"
751 DO i=1,nel
752 j1 = jj(i)
753 j2 = j1+1
754 y1(i)= y1(i)*yfac(i,1)
755 y2(i)= y2(i)*yfac(i,2)
756 fac = rfac(i)
757 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
758 yld(i) = max(yld(i),em20)
759 dydx1(i)= dydx1(i)*yfac(i,1)
760 dydx2(i)= dydx2(i)*yfac(i,2)
761 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
762 yld(i) = yld(i)*max(zero,pfac(i))
763 h(i) = h(i)*max(zero,pfac(i))
764 ENDDO
765 DO i=1,nel
766 j1 = jj(i)
767 j2 = j1+1
768 vartmp(i,2+j1) = ipos1(i)
769 vartmp(i,2+j2) = ipos2(i)
770 ENDDO
771c-------------------------------
772 ELSEIF (fisokin == one) THEN ! ecrouissage cinematique
773c-------------------------------
774#include "vectorize.inc"
775 DO i=1,nel
776 j1 = jj(i)
777 j2 = j1+1
778 fun1 = ifunc(j1)
779 fun2 = ifunc(j2)
780 fac = rfac(i)
781 dydx1(i)=dydx1(i)*yfac(i,1)
782 dydx2(i)=dydx2(i)*yfac(i,2)
783 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
784 y1(i) = tf(npf(fun1)+1)
785 y2(i) = tf(npf(fun2)+1)
786 y1(i) = y1(i)*yfac(i,1)
787 y2(i) = y2(i)*yfac(i,2)
788 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
789 yld(i) = yld(i)*max(zero,pfac(i))
790 h(i) = h(i)*max(zero,pfac(i))
791 ENDDO
792 DO i=1,nel
793 j1 = jj(i)
794 j2 = j1+1
795 vartmp(i,2+j1) = ipos1(i)
796 vartmp(i,2+j2) = ipos2(i)
797 ENDDO
798c-------------------------------
799 ELSE ! ecrouissage mixte
800c-------------------------------
801#include "vectorize.inc"
802 DO i=1,nel
803 j1 = jj(i)
804 j2 = j1+1
805 fun1 = ifunc(j1)
806 fun2 = ifunc(j2)
807 y1(i)= y1(i)*yfac(i,1)
808 y2(i)= y2(i)*yfac(i,2)
809 fac = rfac(i)
810 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
811 yld(i) = max(yld(i),em20)
812 dydx1(i)= dydx1(i)*yfac(i,1)
813 dydx2(i)= dydx2(i)*yfac(i,2)
814 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
815 y1(i) = tf(npf(fun1)+1)
816 y2(i) = tf(npf(fun2)+1)
817 y1(i) = y1(i)*yfac(i,1)
818 y2(i) = y2(i)*yfac(i,2)
819 yld(i) = (one - fisokin) * yld(i) +
820 . fisokin * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
821 yld(i) = yld(i)*max(zero,pfac(i))
822 h(i) = h(i)*max(zero,pfac(i))
823 ENDDO
824 DO i=1,nel
825 j1 = jj(i)
826 j2 = j1+1
827 vartmp(i,2+j1) = ipos1(i)
828 vartmp(i,2+j2) = ipos2(i)
829 ENDDO
830 ENDIF ! ISOKIN
831C-------------------
832C PROJECTION
833C-------------------
834C CRITERE DE VON MISES
835C-------------------------
836 DO i=1,nel
837 h(i) = max(zero,h(i))
838 s1=signxx(i)+signyy(i)
839 s2=signxx(i)-signyy(i)
840 s3=signxy(i)
841 aa(i)=fourth*s1*s1
842 bb(i)=three_over_4*s2*s2+3.*s3*s3
843 svm2(i)= aa(i)+bb(i)
844 IF (inloc == 0) THEN
845 dezz = -(depsxx(i)+depsyy(i))*nu_mnu
846 thk(i) = thk(i) + dezz*thkly(i)*off(i)
847 ENDIF
848 ENDDO
849C-------------------------
850C GATHER PLASTIC FLOW
851C-------------------------
852 nindx=0
853 DO i=1,nel
854 IF(svm2(i) > yld(i)*yld(i) .AND. off(i) == one) THEN
855 nindx=nindx+1
856 index(nindx)=i
857 ENDIF
858 ENDDO
859C
860 IF (nindx > 0) THEN
861 !---------------------------
862 !DEP EN CONTRAINTE PLANE
863 !---------------------------
864#include "vectorize.inc"
865 DO j=1,nindx
866 i=index(j)
867 svm = sqrt(svm2(i))
868 dpla_j(i)=(svm-yld(i))/(g3(i)+h(i))
869 etse(i)= h(i)/(h(i)+e(i))
870 hi(i) = h(i)*(one-fisokin)
871 hk(i) = two_third*h(i)*fisokin
872 ENDDO
873C
874 nu3 = one-nu_mnu
875 DO n=1,niter
876#include "vectorize.inc"
877 DO j=1,nindx
878 i=index(j)
879 dpla_i(i)=dpla_j(i)
880 yld_i =yld(i)+hi(i)*dpla_i(i)
881 dr(i) =half*e(i)*dpla_i(i)/yld_i
882 aaa = three*hk(i)/e(i)
883 nu11 = u_mnu+aaa
884 nu21 = t_pnu+aaa
885 pp(i) = one/(one+dr(i)*nu11)
886 qq(i) = one/(one+dr(i)*nu21)
887 p2 = pp(i)*pp(i)
888 q2 = qq(i)*qq(i)
889 f = aa(i)*p2+bb(i)*q2-yld_i*yld_i
890 df =-(aa(i)*nu11*p2*pp(i)+bb(i)*nu21*q2*qq(i))
891 . *(e(i)-two*dr(i)*hi(i))/yld_i
892 . -two*hi(i)*yld_i
893 df = sign(max(abs(df),em20),df)
894 IF(dpla_i(i) > zero) THEN
895 dpla_j(i)=max(zero,dpla_i(i)-f/df)
896 ELSE
897 dpla_j(i)=zero
898 ENDIF
899 ENDDO ! NINDX
900 ENDDO ! NITER
901C------------------------------------------
902C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
903C------------------------------------------
904#include "vectorize.inc"
905 DO j=1,nindx
906 i=index(j)
907 pla(i) = pla(i) + dpla_i(i)
908 s1=(signxx(i)+signyy(i))*pp(i)
909 s2=(signxx(i)-signyy(i))*qq(i)
910 signxx(i)=half*(s1+s2)
911 signyy(i)=half*(s1-s2)
912 signxy(i)=signxy(i)*qq(i)
913 IF (inloc == 0) THEN
914 dezz = - nu3*dr(i)*s1/e(i)
915 thk(i) = thk(i) + dezz*thkly(i)*off(i)
916 ENDIF
917 !-----for kin hardening
918 yld(i) =yld(i)+hi(i)*dpla_i(i)
919 ENDDO
920 ENDIF
921C=======================================================================
922 ENDIF ! VP flag
923C=======================================================================
924c----------------------------------------------------------------------
925c Element failure
926c----------------------------------------------------------------------
927 IF (ifail == 1) THEN
928 IF (inloc > 0) THEN
929 DO i=1,nel
930 IF (epsmax < ep20) dmg(i) = planl(i)/epsmax
931 IF (off(i) == one .AND. planl(i) > epsmax) off(i) = four_over_5
932 ENDDO
933 ELSE
934 DO i=1,nel
935 IF (epsmax < ep20) dmg(i) = pla(i)/epsmax
936 IF (off(i) == one .AND. pla(i) > epsmax) off(i) = four_over_5
937 ENDDO
938 ENDIF
939 ELSEIF (ifail == 2) THEN
940 IF (inloc > 0) THEN
941 DO i=1,nel
942 IF (epsmax < ep20) dmg(i) = max(dmg(i),planl(i)/epsmax)
943 IF (off(i) == one .AND. (planl(i) > epsmax .OR. epst(i) > epsf)) off(i) = four_over_5
944 ENDDO
945 ELSE
946 DO i=1,nel
947 IF (epsmax < ep20) dmg(i) = max(dmg(i),pla(i)/epsmax)
948 IF (off(i) == one .AND. (pla(i) > epsmax .OR. epst(i) > epsf)) off(i) = four_over_5
949 ENDDO
950 ENDIF
951 ENDIF
952c----------------------------------------------------------------------
953 IF (impl_s > 0) THEN
954 IF (ikt > 0 ) THEN
955 DO i = 1,nel
956 IF (dpla_i(i) > zero) THEN
957c ...... Parameter d(gama)
958 gama_imp(i)= three_half*dpla_i(i)/yld(i)
959c ...... HK,HH---
960 signor(i,4)=fisokin*h(i)
961 signor(i,5)=(one-fisokin)*h(i)
962c ...... Deviatoric stresses shifted by modified back stress ->ksi
963 signor(i,1)=third*(two*signxx(i)-signyy(i))
964 signor(i,2)=third*(two*signyy(i)-signxx(i))
965 signor(i,3)=two*signxy(i)
966 ELSE
967 gama_imp(i) = zero
968 END IF
969 END DO
970 END IF !(IKT > 0) THEN
971 END IF
972C------------------------------------------
973C PLASTIC STRAIN FILTERING
974C------------------------------------------
975 IF (vp == 1) THEN
976 dtinv = one/max(timestep,em20)
977 DO i=1,nel
978 plap(i) = asrate*dpla_i(i)*dtinv + (one - asrate)*uvar(i,2)
979 uvar(i,2) = plap(i)
980 ENDDO
981 ENDIF
982C------------------------------------------
983C ECROUISSAGE CINE
984C------------------------------------------
985 IF (fisokin > zero) THEN
986 DO i=1,nel
987 hkin = fisokin*h(i)
988 alpha = hkin*dpla_i(i)/yld(i)
989 sigpxx = alpha*signxx(i)
990 sigpyy = alpha*signyy(i)
991 sigpxy = alpha*signxy(i)
992
993 sigbxx(i) = sigbxx(i) + sigpxx
994 sigbyy(i) = sigbyy(i) + sigpyy
995 sigbxy(i) = sigbxy(i) + sigpxy
996C
997 signxx(i) = signxx(i) + sigbxx(i)
998 signyy(i) = signyy(i) + sigbyy(i)
999 signxy(i) = signxy(i) + sigbxy(i)
1000 ENDDO
1001 ENDIF
1002C--------------------------------
1003C HARDENING MODULUS
1004C--------------------------------
1005 DO i=1,nel
1006 hardm(i) = h(i)
1007 ENDDO
1008C--------------------------------
1009C NON-LOCAL THICKNESS VARIATION
1010C--------------------------------
1011 IF (inloc > 0) THEN
1012 DO i = 1,nel
1013 IF (loff(i) == one) THEN
1014 svm = sqrt(signxx(i)*signxx(i)
1015 . + signyy(i)*signyy(i)
1016 . - signxx(i)*signyy(i)
1017 . + three*signxy(i)*signxy(i))
1018 dezz = max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(svm,em20)
1019 dezz = -nux*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e(i)) - dezz
1020 thk(i) = thk(i) + dezz*thkly(i)*off(i)
1021 ENDIF
1022 ENDDO
1023 ENDIF
1024c-----------
1025 RETURN
1026 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps36c
subroutine sigeps36c(nel, nuvar, nvartmp, mfunc, kfunc, npf, iplas, tf, timestep, uparam, rho0, thkly, israte, asrate, epsd_pg, epsd, epspxx, epspyy, epspxy, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, soundsp, viscmax, thk, pla, uvar, vartmp, off, ipm, imat, etse, gs, yld, dpla_i, gama_imp, signor, shf, hardm, facyldi, inloc, dplanl, dmg, planl, l_planl, sigbxx, sigbyy, sigbxy, loff)
Definition sigeps36c.F:47
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