OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps52c.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!|| sigeps52c ../engine/source/materials/mat/mat052/sigeps52c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| table_interp ../engine/source/tools/curve/table_tools.F
30!|| table_interp_law76 ../engine/source/tools/curve/table_tools.F
31!||--- uses -----------------------------------------------------
32!|| interface_table_mod ../engine/share/modules/table_mod.F
33!|| table_mod ../engine/share/modules/table_mod.F
34!||====================================================================
35 SUBROUTINE sigeps52c(
36 1 NEL0 ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,
37 2 NPF ,NPT0 ,IPT ,IFLAG ,ASRATE ,
38 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,
39 3 AREA ,EINT ,THKLY ,L_DMG ,DMG ,
40 4 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
41 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
42 6 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
43 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
44 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
45 9 SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
46 A SOUNDSP,VISCMAX,THK ,PLA ,UVAR ,
47 B OFF ,NGL ,IPM ,MAT ,ETSE ,
48 C GS ,YLD ,EPSD_PG ,TABLE ,EPSD )
49C-----------------------------------------------
50 USE table_mod
51 USE interface_table_mod
52C-----------------------------------------------
53C I m p l i c i t T y p e s
54C-----------------------------------------------
55#include "implicit_f.inc"
56C-----------------------------------------------
57C G l o b a l P a r a m e t e r s
58C-----------------------------------------------
59#include "mvsiz_p.inc"
60C---------+---------+---+---+--------------------------------------------
61C VAR | SIZE |TYP| RW| DEFINITION
62C---------+---------+---+---+--------------------------------------------
63C NEL0 | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL0
64C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
65C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
66C---------+---------+---+---+--------------------------------------------
67C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
68C IFUNC | NFUNC | I | R | FUNCTION INDEX
69C NPF | * | I | R | FUNCTION ARRAY
70C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
71C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
72C IFLAG | * | I | R | GEOMETRICAL FLAGS
73C TF | * | F | R | FUNCTION ARRAY
74C---------+---------+---+---+--------------------------------------------
75C TIME | 1 | F | R | CURRENT TIME
76C TIMESTEP| 1 | F | R | CURRENT TIME STEP
77C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAYSIGY
78C RHO0 | NEL0 | F | R | INITIAL DENSITY
79C AREA | NEL0 | F | R | AREA
80C EINT | 2*NEL0 | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
81C THKLY | NEL0 | F | R | LAYER THICKNESS
82C EPSPXX | NEL0 | F | R | STRAIN RATE XX
83C EPSPYY | NEL0 | F | R | STRAIN RATE YY
84C ... | | | |
85C DEPSXX | NEL0 | F | R | STRAIN INCREMENT XX
86C DEPSYY | NEL0 | F | R | STRAIN INCREMENT YY
87C ... | | | |
88C EPSXX | NEL0 | F | R | STRAIN XX
89C EPSYY | NEL0 | F | R | STRAIN YY
90C ... | | | |
91C SIGOXX | NEL0 | F | R | OLD ELASTO PLASTIC STRESS XX
92C SIGOYY | NEL0 | F | R | OLD ELASTO PLASTIC STRESS YY
93C ... | | | |
94C---------+---------+---+---+--------------------------------------------
95C SIGNXX | NEL0 | F | W | NEW ELASTO PLASTIC STRESS XX
96C SIGNYY | NEL0 | F | W | NEW ELASTO PLASTIC STRESS YY
97C ... | | | |
98C SIGVXX | NEL0 | F | W | VISCOUS STRESS XX
99C SIGVYY | NEL0 | F | W | VISCOUS STRESS YY
100C ... | | | |
101C SOUNDSP | NEL0 | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
102C VISCMAX | NEL0 | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
103C---------+---------+---+---+--------------------------------------------
104C THK | NEL0 | F |R/W| THICKNESS
105C PLA | NEL0 | F |R/W| PLASTIC STRAIN
106C UVAR |NEL0*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
107C OFF | NEL0 | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
108C---------+---------+---+---+--------------------------------------------
109C C o m m o n B l o c k s
110C-----------------------------------------------
111#include "param_c.inc"
112C tmp +1
113C-----------------------------------------------
114C I N P U T A r g u m e n t s
115C-----------------------------------------------
116C
117 INTEGER NEL0, NUPARAM, NUVAR, NPT0, IPT,IFLAG(*),
118 . NGL(NEL0),MAT(NEL0),IPM(NPROPMI,*),NSG
119 my_real TIME,TIMESTEP(NEL0),UPARAM(*),
120 . AREA(NEL0),RHO0(NEL0),EINT(NEL0,2),
121 . THKLY(NEL0),PLA(NEL0),
122 . EPSPXX(NEL0),EPSPYY(NEL0),
123 . EPSPXY(NEL0),EPSPYZ(NEL0),EPSPZX(NEL0),
124 . DEPSXX(NEL0),DEPSYY(NEL0),
125 . DEPSXY(NEL0),DEPSYZ(NEL0),DEPSZX(NEL0),
126 . EPSXX(NEL0) ,EPSYY(NEL0) ,
127 . EPSXY(NEL0) ,EPSYZ(NEL0) ,EPSZX(NEL0) ,
128 . SIGOXX(NEL0),SIGOYY(NEL0),
129 . SIGOXY(NEL0),SIGOYZ(NEL0),SIGOZX(NEL0),
130 . gs(*)
131 TYPE(ttable) TABLE(*)
132 INTEGER, INTENT(IN) :: L_DMG
133C-----------------------------------------------
134C O U T P U T A r g u m e n t s
135C-----------------------------------------------
136 my_real
137 . SIGNXX(NEL0),SIGNYY(NEL0),
138 . SIGNXY(NEL0),SIGNYZ(NEL0),SIGNZX(NEL0),
139 . SIGVXX(NEL0),SIGVYY(NEL0),
140 . SIGVXY(NEL0),SIGVYZ(NEL0),SIGVZX(NEL0),
141 . SOUNDSP(NEL0),VISCMAX(NEL0),ETSE(NEL0)
142C-----------------------------------------------
143C I N P U T O U T P U T A r g u m e n t s
144C-----------------------------------------------
145 my_real uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
146 my_real, DIMENSION(NEL0,L_DMG), INTENT(INOUT) :: dmg
147 my_real ,INTENT(IN) :: asrate
148 my_real ,DIMENSION(NEL0) ,INTENT(IN) :: epsd_pg
149 my_real ,DIMENSION(NEL0) ,INTENT(INOUT) :: epsd
150C-----------------------------------------------
151C VARIABLES FOR FUNCTION INTERPOLATION
152C-----------------------------------------------
153 INTEGER NPF(*), MFUNC, KFUNC(MFUNC),ITER, IFLA
154 my_real FINTER ,TF(*)
155 EXTERNAL FINTER
156C EXTERNAL FINTER
157C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
158C Y : y = f(x)
159C X : x
160C DYDX : f'(x) = dy/dx
161C IFUNC(J): FUNCTION INDEX
162C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
163C NPF,TF : FUNCTION PARAMETER
164C-----------------------------------------------
165C L o c a l V a r i a b l e s
166C-----------------------------------------------
167 INTEGER I,J,NRATE(MVSIZ),J1,J2,NINDX,NMAX,IADBUF,NFUNC,
168 . iad1(mvsiz),ipos1(mvsiz),ilen1(mvsiz),
169 . iad2(mvsiz),ipos2(mvsiz),ilen2(mvsiz),
170 . jj(mvsiz),index(mvsiz),nratem,
171 . nrate1,ifunc(mfunc),iadbufv,
172 . nfuncm,israte,iyeild_tab,
173 . itab(100),ntable,nxh,ndim,ipos,nxk,mx
174 my_real
175 . e,a1,a2,g,g3,
176 . nu
177 my_real
178 . f(mvsiz), sigm(mvsiz), epsp1(mvsiz),
179 . lamda(mvsiz), fstar(mvsiz),p0(mvsiz),pn(mvsiz),
180 . et,yeild0,csd,
181 . visp, q1, q2, q3,
182 . sn, epsn, fi,fn,
183 . fc,ff,fu,
184 . sigxx(mvsiz),sigyy(mvsiz),sigxy(mvsiz),
185 . sigzx(mvsiz),sigyz(mvsiz),en
186 my_real
187 . p,vm,r,fac,
188 . dp11,dp22,dp33,dp12,dp13,dp23,a22,a11,
189 . d11,d22,d33,d12,d13,d23,dcrf,dcrm, dsepp,dcd,lam1,
190 . c2,df,coh,sih,va,crit,conv,fg(mvsiz),fn1(mvsiz),yldc,
191 . c11,p1, var,c1,nn1,dezz,
192 . dtinv,sqr22, va1, va2, vm1, sigm1, sigm2, a21,yp1,yp2,
193 . df1,df2, xfac,yfac,dx2,xx(2),dft,yy
194C
195c
196C-----------------------------------------------
197C USER VARIABLES INITIALIZATION
198C-----------------------------------------------
199
200 ntable = ipm(226,mat(1))
201 iyeild_tab = 0
202 IF(ntable > 0) THEN
203 iyeild_tab = 1
204 itab(1)=ipm(227,mat(1))
205 iadbuf = ipm(7,mat(1))-1
206 xfac = uparam(iadbuf + 22 )
207 yfac = uparam(iadbuf + 23 )
208 ENDIF
209 mx = mat(1)
210 iadbufv = ipm(7,mx)-1
211 e = uparam(iadbufv+1)
212 nu = uparam(iadbufv+2)
213 yeild0 = uparam(iadbufv+3)
214 et = uparam(iadbufv+4)
215 en = uparam(iadbufv+5)
216 csd = one/uparam(iadbufv+6)
217 visp = uparam(iadbufv+7)
218 q1 = uparam(iadbufv+8)
219 q2 = uparam(iadbufv+9)
220 q3 = uparam(iadbufv+10)
221 sn = uparam(iadbufv+11)
222 IF(sn==zero) sn = ep20
223 epsn = uparam(iadbufv+12)
224 fi = uparam(iadbufv+13)
225 fn = uparam(iadbufv+14)
226 fc = uparam(iadbufv+15)
227 ff = uparam(iadbufv+16)
228 fu = uparam(iadbufv+17)
229 ifla = nint(uparam(iadbufv+18))
230 israte = nint(uparam(iadbufv+19))
231C
232 a1 = e/(one -nu**2)
233 a2 = nu*a1
234 g = half*e/(one + nu)
235 g3 = three*g
236 c1 = e*(one -nu) /((one + nu)*(one - two*nu))
237 c2 = c1*nu/(one - nu)
238 nn1 = nu/(1. - nu)
239C
240 IF(time==zero)THEN
241 DO i=1,nel0
242 epsp1(i) = zero
243 pla(i) = zero
244 IF (en<one) pla(i) = em20
245 dmg(i,4) = fi
246 dmg(i,5) = fi
247 yld(i) = yeild0
248 uvar(i,1) = yeild0
249 ENDDO
250 IF(iyeild_tab > 0) THEN
251 DO i=1,nel0
252 xx(1)=zero
253 xx(2)=zero
254 CALL table_interp (table(itab(1)),xx,yy)
255 yld(i) = yy *yfac
256 uvar(i,1) = yld(i)
257 ENDDO
258 ENDIF
259 ENDIF
260 sqr22 = one/sqrt(two*pi)
261C-----------------------------------------------
262C
263 DO i=1,nel0
264C
265 fg(i) = dmg(i,2)
266 fn1(i) = dmg(i,3)
267 f(i) = dmg(i,4)
268 fstar(i) = dmg(i,5)
269 sigm(i) = uvar(i,1)
270 epsp1(i) = zero
271C
272 signxx(i) = sigoxx(i) + a1 * depsxx(i) + a2 * depsyy(i)
273 signyy(i) = sigoyy(i) + a2 * depsxx(i) + a1 * depsyy(i)
274 signxy(i) = sigoxy(i) + g * depsxy(i)
275 signyz(i) = sigoyz(i) + gs(i) *depsyz(i)
276 signzx(i) = sigozx(i) + gs(i) *depszx(i)
277C
278 soundsp(i) = sqrt(a1/rho0(i))
279 viscmax(i) = zero
280 etse( i )= one
281C-------------------
282C STRAIN RATE
283C-------------------
284 IF (israte == 0) THEN
285 epsd(i) = half*( abs(epspxx(i)+epspyy(i))
286 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
287 . + epspxy(i)*epspxy(i) ) )
288 ELSE
289 epsd(i) = asrate*epsd_pg(i) + (one-asrate)*epsd(i)
290 END IF
291
292 dezz = -nn1*(depsxx(i) + depsyy(i))
293 thk(i) = thk(i) + dezz*thkly(i)*off(i)
294 ENDDO
295C-------------------
296C CRITERE
297C-----
298C
299C VISCOPLASTIC
300 IF(ifla==0)THEN
301C-----------
302C IFLAG=0,
303C-------------
304 DO i=1,nel0
305C
306 IF(off(i)==one)THEN
307C
308 IF(f(i)<=fc)THEN
309 fstar(i) = f(i)
310 df = one
311 ELSE
312 df = (fu-fc)/(ff-fc)
313 fstar(i) = fc + df*(f(i)-fc)
314 ENDIF
315 vm = (signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i)
316 . + three * (signxy(i)**2 + signzx(i)**2 + signyz(i)**2))
317 vm = sqrt(vm)
318 pn(i) = (signxx(i)+signyy(i)) * third
319 sigm1 = one/max(em20, sigm(i))
320 var = three_half * q2* pn(i) * sigm1
321 var = exp(var)
322 coh = half * (var + one/max(em20,var))
323 sih = half * (var - one/max(em20,var))
324C
325 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
326 va = sqrt(max(zero,va))
327 yldc = sigm(i) * va
328 yld(i ) = yldc
329 crit = vm - yldc
330C
331 IF(vm < yldc .OR. yldc == zero )THEN
332 lamda(i)=zero
333 ELSE
334C.. stockage des contraintes ....
335 lamda(i) = zero
336 sigxx(i) = signxx(i)
337 sigyy(i) = signyy(i)
338 sigxy(i) = signxy(i)
339 sigzx(i) = signzx(i)
340 sigyz(i) = signyz(i)
341 dp11 = zero
342 dp22 = zero
343 dp33 = zero
344 dp12 = zero
345 dp13 = zero
346 dp23 = zero
347 dtinv =timestep(i)/max(timestep(i)**2,em20)
348 IF(f(i)==one)THEN
349 a21 = ep20
350 ELSE
351 a21 = sigm1 /(one - f(i))
352 ENDIF
353 a11 = fn*exp(-half*((pla(i)-epsn)/sn)**2)/sn
354 a11 = a11*sqr22
355 IF(en/=1)THEN
356 DO iter=1,5
357 vm1 = one/max(em20, vm)
358 va1 = one/max(em20, va)
359 va2 = half * q1*q2*fstar(i)*sih*va1
360 d11 = half * (two*signxx(i)-signyy(i))*vm1 + va2
361 d22 = half * (two*signyy(i)-signxx(i))*vm1 + va2
362 d33 = half * (-signxx(i)-signyy(i))*vm1 + va2
363 d12 = three * signxy(i)*vm1
364 d13 = three * signzx(i)*vm1
365 d23 = three * signyz(i)*vm1
366 a22 = d11 * signxx(i) + d22 * signyy(i) +
367 . two*(d12 * signxy(i) + d13 * signzx(i) +
368 . d23 * signyz(i))
369C
370 a22 = a22 * a21
371 dcrf = - sigm(i)*(q3*fstar(i)*df - q1*coh*df)*va1
372 dcrm = - va - three*va2*pn(i)*sigm1
373 dsepp = et * en * pla(i) ** (en - one)*
374 . (one + ( epsd(i)*csd)**visp)
375 dcd = c1*(d11**2 + d22**2 + d33**2) +
376 . two*c2*(d11*d22 + d11*d33 + d22*d33) +
377 . two*g * d12**2 + 2*gs(i) * ( d13**2 + d23**2 )
378C
379 lam1= dcd - dcrm * dsepp * a22 -
380 . dcrf *((one -f(i)) * (d11 + d22 + d33) + a11*a22)
381 IF(lam1/=zero)lamda(i)= max(zero,( vm - yldc ))/lam1
382 dp11 = dp11 + lamda(i) * d11
383 dp22 = dp22 + lamda(i) * d22
384 dp33 = dp33 + lamda(i) * d33
385 dp12 = dp12 + lamda(i) * d12
386 dp13 = dp13 + lamda(i) * d13
387 dp23 = dp23 + lamda(i) * d23
388C NEW STRE SS
389 signxx(i)= sigxx(i) - a1 * dp11 - a2 * dp22
390 signyy(i)= sigyy(i) - a1 * dp22 - a2 * dp11
391 signxy(i)= sigxy(i) - two*g * dp12
392 signzx(i)= sigzx(i) - two*gs(i) * dp13
393 signyz(i)= sigyz(i) - two*gs(i) * dp23
394C....... EFF ECTIVE PLASTIC DEFORMATION.
395 epsp1(i)= signxx(i) * dp11 + signyy(i) * dp22 +
396 . two*(signxy(i) * dp12 + signzx(i) * dp13 +
397 . signyz(i) * dp23 )
398 epsp1(i)= epsp1(i)* a21
399 epsp1(i) = max(zero, epsp1(i))
400 vm = signxx(i)**2 + signyy(i)**2
401 . + three*(signxy(i)**2 + signzx(i)**2 +signyz(i)**2)
402 . - signxx(i)*signyy(i)
403 vm = sqrt(vm)
404C
405 pn(i) = (signxx(i)+signyy(i))*third
406 var = three_half * q2 * pn(i) * sigm1
407 var = exp(var)
408 coh = half * ( var + one/max(em20,var))
409 sih = half * ( var - one/max(em20,var))
410 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
411 va= sqrt(max(zero,va))
412 yldc = sigm(i) * va
413 ENDDO ! ite
414 ELSE
415 DO iter=1,5
416 vm1 = one/max(em20, vm)
417 va1 = one/max(em20, va)
418 va2 = half * q1*q2*fstar(i)*sih*va1
419 d11 = half * (two*signxx(i)-signyy(i))*vm1 + va2
420 d22 = half * (two*signyy(i)-signxx(i))*vm1 + va2
421 d33 = half*( - signxx(i)-signyy(i))*vm1 + va2
422 d12 = three * signxy(i) * vm1
423 d13 = three * signzx(i) * vm1
424 d23 = three * signyz(i) * vm1
425 a22 = d11 * signxx(i) + d22 * signyy(i) +
426 . two*(d12 * signxy(i) + d13 * signzx(i) +
427 . d23 * signyz(i))
428 a22 = a22 * a21
429 dcrf = - sigm(i)*(q3*fstar(i)*df - q1*coh*df)*va1
430 dcrm = - va - three*va2*pn(i)*sigm1
431 dsepp = et*(one+ ( epsd(i)*csd)**visp)
432 dcd = c1*(d11**2 + d22**2 + d33**2) +
433 . two*c2*(d11*d22 + d11*d33 + d22*d33) +
434 . two*g*d12**2 + two*gs(i)*( d13**2 + d23**2)
435C
436 lam1= dcd - dcrm * dsepp * a22 -
437 . dcrf *((one - f(i)) * (d11 + d22+d33) + a11*a22)
438 IF(lam1/=zero) lamda(i)= max(zero,( vm - yldc ))/lam1
439 dp11 = dp11 + lamda(i) * d11
440 dp22 = dp22 + lamda(i) * d22
441 dp33 = dp33 + lamda(i) * d33
442 dp12 = dp12 + lamda(i) * d12
443 dp13 = dp13 + lamda(i) * d13
444 dp23 = dp23 + lamda(i) * d23
445C NEW ST RESS
446 signxx(i)= sigxx(i) - a1 * dp11 - a2 * dp22
447 signyy(i)= sigyy(i) - a1 * dp22 - a2 * dp11
448 signxy(i)= sigxy(i) - two*g * dp12
449 signzx(i)= sigzx(i) - two*gs(i) * dp13
450 signyz(i)= sigyz(i) - two*gs(i) * dp23
451C....... E FFECTIVE PLASTIC DEFORMATION.
452 epsp1(i)= signxx(i) * dp11 + signyy(i) * dp22 +
453 . two*(signxy(i) * dp12 + signzx(i) * dp13 +
454 . signyz(i) * dp23 )
455 epsp1(i)= epsp1(i)*a21
456 epsp1(i) = max(zero, epsp1(i))
457 vm = signxx(i)**2 + signyy(i)**2
458 . + three*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
459 . - signxx(i)*signyy(i)
460 vm = sqrt(vm)
461C
462 pn(i) = (signxx(i)+signyy(i))*third
463 var = three_half * q2 * pn(i) * sigm1
464 var = exp(var)
465 coh = half * ( var + one/max(em20,var))
466 sih = half * ( var - one/max(em20,var))
467 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
468 va = sqrt(max(zero,va))
469 yldc= sigm(i) * va
470 ENDDO ! ite
471 ENDIF !LAM
472 yld(i ) = yldc
473 etse( i) = dsepp/(dsepp+e)
474 dezz = nn1*(dp11 +dp22) + dp33
475 thk(i) = thk(i) + dezz*thkly(i)*off(i)
476 fg(i) = fg(i) + (one - f(i)) * (dp11 + dp22 +dp33)
477 fn1(i) = fn1(i) + a11 * epsp1(i)
478 pla(i) = pla(i) + epsp1(i)
479 f(i) = fi + fg(i) + fn1(i)
480 IF(q1==zero.AND.
481 . q2==zero.AND.
482 . q3==zero)f(i)=zero
483 IF(f(i)<zero)f(i)=zero
484 sigm(i)= (yeild0 + et*pla(i)**en)* (one+(epsd(i)*csd)**visp)
485C....COMPUTE F*.....
486 IF(f(i)<=fc)THEN
487 fstar(i)=f(i)
488 ELSE
489 IF(ff==fc)THEN
490 fstar(i) = ep20
491 ELSE
492 fstar(i) =fc+(fu-fc)*(f(i)-fc)/(ff-fc)
493 ENDIF
494 ENDIF
495 uvar(i,1) = sigm(i)
496c..... DEVIATORIC STRESS
497 ENDIF
498 IF(off(i)==one.AND.fstar(i)>=ff) off(i)=four_over_5
499 ENDIF
500 ENDDO ! NEL0
501 ELSEIF(ifla==1)THEN
502C----
503CIFLA set to 1
504C-----
505 DO i=1,nel0
506 IF(off(i)==one)THEN
507 IF(f(i)<=fc)THEN
508 fstar(i) = f(i)
509 df = one
510 ELSE
511 df =(fu - fc)/(ff - fc)
512 fstar(i) = fc + df*(f(i)-fc)
513 ENDIF
514C
515 vm = (signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i)
516 . + three*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2))
517 pn(i) = (signxx(i)+signyy(i))*third
518 sigm1 = one/max(em20,sigm(i))
519 sigm2 = sigm1**2
520 var = three_half * q2 * pn(i) * sigm1
521 var = exp(var)
522 coh = half *(var + one/max(em20,var))
523 sih = half *(var - one/max(em20,var))
524C
525 IF(pn(i)<=zero) THEN
526 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
527 ELSE
528 va= one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
529 ENDIF
530C
531 lamda(i) = zero
532 yldc = sigm(i)
533 crit = vm * sigm2 - va
534 yld( i ) = sigm(i) * sqrt(max(zero,va))
535 IF(crit < zero .OR. yld(i) == zero )THEN
536 lamda(i)=zero
537 ELSE
538C.. stockage des contraintes ....
539 sigxx(i) = signxx(i)
540 sigyy(i) = signyy(i)
541 sigxy(i) = signxy(i)
542 sigzx(i) = signzx(i)
543 sigyz(i) = signyz(i)
544 dp11 = zero
545 dp22 = zero
546 dp33 = zero
547 dp12 = zero
548 dp13 = zero
549 dp23 = zero
550 dtinv =timestep(i)/max(timestep(i)**2,em20)
551C
552 IF(f(i)==one)THEN
553 a21 = ep20
554 ELSE
555 a21 = sigm1 /(one - f(i))
556 ENDIF
557 a11 = fn*exp(-half*((pla(i)-epsn)/sn)**2)/sn
558 a11 = a11*sqr22
559C
560 IF(en/=one) THEN
561 DO iter=1,5
562 IF(pn(i)<=zero)THEN
563 d11 = (two*signxx(i)-signyy(i))*sigm2
564 d22 = (two*signyy(i)-signxx(i))*sigm2
565 d33 = ( - signxx(i) - signyy(i))*sigm2
566 d12 = six*signxy(i)*sigm2
567 d13 = six*signzx(i)*sigm2
568 d23 = six*signyz(i)*sigm2
569 dcrm = -two*vm*sigm1*sigm2
570 dcrf = two*q1*df - two*q3*df*fstar(i)
571 ELSE
572 va2 = q1 * q2 * fstar(i) * sih * sigm1
573 d11 = (two*signxx(i) - signyy(i))*sigm2 + va2
574 d22 = (two*signyy(i) - signxx(i))*sigm2 + va2
575 d33 = ( - signxx(i) - signyy(i))*sigm2 + va2
576 d12 = six*signxy(i)*sigm2
577 d13 = six*signzx(i)*sigm2
578 d23 = six*signyz(i)*sigm2
579 dcrm = -two*vm*sigm1*sigm2 - three*va2*pn(i)*sigm1
580 dcrf = two * q1 * df * coh - two*q3*df*fstar(i)
581 ENDIF
582C
583 a22 = d11*signxx(i) + d22*signyy(i) +
584 . three*(d12*signxy(i)+ d13*signzx(i)+ d23*signyz(i))
585 a22 = a22*a21
586 dsepp = et * en * pla(i) ** (en - 1)*
587 . (one + ( epsd(i)*csd)**visp)
588 dcd = c1*(d11**2 + d22**2 + d33**2) +
589 . two*c2*(d11*d22 + d11*d33 + d22*d33) +
590 . two*g*d12**2 + two*gs(i)*( d13**2 + d23**2)
591C
592 lam1= dcd - dcrm * dsepp*a22 -
593 . dcrf*((one-f(i))*(d11 + d22 +d33) + a11*a22)
594 IF(lam1/=zero) lamda(i) = max(zero,(vm*sigm2 - va))/lam1
595C.... PLASTIC DEFORMATION
596 dp11 = dp11 + lamda(i) * d11
597 dp22 = dp22 + lamda(i) * d22
598 dp33 = dp33 + lamda(i) * d33
599 dp12 = dp12 + lamda(i) * d12
600 dp13 = dp13 + lamda(i) * d13
601 dp23 = dp23 + lamda(i) * d23
602C NEW STRESS
603 signxx(i) = sigxx(i) - a1*dp11 - a2*dp22
604 signyy(i) = sigyy(i) - a1*dp22 - a2*dp11
605 signxy(i) = sigxy(i) - two*g*dp12
606 signzx(i)= sigzx(i) - two*gs(i)*dp13
607 signyz(i)= sigyz(i) - two*gs(i)*dp23
608C....... EFFECTIVE PLASTIC DEFORMATION.
609 epsp1(i)= signxx(i)*dp11 + signyy(i)*dp22 +
610 . two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
611 epsp1(i) = epsp1(i)*a21
612 epsp1(i) = max(zero, epsp1(i))
613 vm= signxx(i)**2 + signyy(i)**2
614 . + three*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
615 . - signxx(i)*signyy(i)
616C
617 pn(i) = (signxx(i) + signyy(i)) * third
618 var = three_half * q2*pn(i)*sigm1
619 var= exp(var)
620 coh = half * (var + one/max(em20,var))
621 sih = half * (var - one/max(em20,var))
622C
623 IF(pn(i)<=zero) THEN
624 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
625 ELSE
626 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
627 ENDIF
628C
629 ENDDO ! it
630C
631 ELSE
632 DO iter=1,5
633 IF(pn(i)<=zero)THEN
634 d11 = (two*signxx(i)-signyy(i))*sigm2
635 d22 = (two*signyy(i)-signxx(i))*sigm2
636 d33 = ( - signxx(i) - signyy(i))*sigm2
637 d12 = six*signxy(i)*sigm2
638 d13 = six*signzx(i)*sigm2
639 d23 = six*signyz(i)*sigm2
640 dcrm = -two*vm*sigm1*sigm2
641 dcrf = two*q1*df - two*q3*df*fstar(i)
642 ELSE
643 va2 = q1 * q2 * fstar(i) * sih * sigm1
644 d11 = (two*signxx(i) - signyy(i))*sigm2 + va2
645 d22 = (two*signyy(i) - signxx(i))*sigm2 + va2
646 d33 = ( - signxx(i) - signyy(i))*sigm2 + va2
647 d12 = six*signxy(i)*sigm2
648 d13 = six*signzx(i)*sigm2
649 d23 = six*signyz(i)*sigm2
650 dcrm = -two*vm*sigm1*sigm2 - three*va2*pn(i)*sigm1
651 dcrf = two * q1 * df * coh - two*q3*df*fstar(i)
652 ENDIF
653C
654 a22 = d11*signxx(i) + d22*signyy(i) +
655 . three*(d12*signxy(i)+ d13*signzx(i)+ d23*signyz(i))
656 a22 = a22*a21
657 dsepp = et*(one+ ( epsd(i)*csd)**visp)
658 dcd = c1*(d11**2 + d22**2 + d33**2) +
659 . two*c2*(d11*d22 + d11*d33 + d22*d33) +
660 . two*g*d12**2 + 2*gs(i)*( d13**2 + d23**2)
661C
662 lam1= dcd - dcrm * dsepp*a22 -
663 . dcrf*((one -f(i))*(d11 + d22 + d33) + a11*a22)
664 IF(lam1/=zero)lamda(i) = max(zero,(vm*sigm2 - va))/lam1
665C.... PLASTIC DEFORMATION
666 dp11 = dp11 + lamda(i) * d11
667 dp22 = dp22 + lamda(i) * d22
668 dp33 = dp33 + lamda(i) * d33
669 dp12 = dp12 + lamda(i) * d12
670 dp13 = dp13 + lamda(i) * d13
671 dp23 = dp23 + lamda(i) * d23
672C NEW STRESS
673 signxx(i) = sigxx(i) - a1*dp11 - a2*dp22
674 signyy(i) = sigyy(i) - a1*dp22 - a2*dp11
675 signxy(i) = sigxy(i) - two*g*dp12
676 signzx(i)= sigzx(i) - two*gs(i)*dp13
677 signyz(i)= sigyz(i) - two*gs(i)*dp23
678C....... EFFECTIVE PLASTIC DEFORMATION.
679 epsp1(i)= signxx(i)*dp11 + signyy(i)*dp22 +
680 . two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
681 epsp1(i) = epsp1(i)*a21
682C
683 epsp1(i) =max(zero, epsp1(i))
684 vm= signxx(i)**2 + signyy(i)**2
685 . + three*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
686 . - signxx(i)*signyy(i)
687C
688 pn(i) = (signxx(i) + signyy(i)) * third
689 var = three_half * q2*pn(i)*sigm1
690 var = exp(var)
691 coh = half * (var + one/max(em20,var))
692 sih = half * (var - one/max(em20,var))
693C
694 IF(pn(i)<=zero) THEN
695 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
696 ELSE
697 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
698 ENDIF
699C
700 ENDDO ! it
701 ENDIF
702 yld(i ) = sigm(i)*sqrt(max(zero,va))
703 etse( i) = dsepp/(dsepp+e)
704 dezz = nn1*(dp11 +dp22) + dp33
705 thk(i) = thk(i) + dezz*thkly(i)*off(i)
706 fg(i) = fg(i) + (one - f(i))*(dp11 + dp22 + dp33)
707 fn1(i) = fn1(i) + a11*epsp1(i)
708 pla(i) = pla(i) + epsp1(i)
709 f(i) = fi + fg(i) + fn1(i)
710 IF(q1==zero.and.
711 . q2==zero.and.
712 . q3==zero)f(i)=zero
713 IF(f(i)<zero)f(i)=zero
714 sigm(i)= (yeild0 + et*pla(i)**en)*(one+(epsd(i)*csd)**visp)
715C....COMPUTE F*.....
716 IF(f(i)<=fc)THEN
717 fstar(i) = f(i)
718 ELSE
719 IF(ff==fc)THEN
720 fstar(i) = ep20
721 ELSE
722 fstar(i) = fc+ (fu-fc)*(f(i)-fc)/(ff-fc)
723 ENDIF
724 ENDIF
725C
726 dmg(i,1) = f(i)/ff
727 dmg(i,2) = fg(i)
728 dmg(i,3) = fn1(i)
729 dmg(i,4) = f(i)
730 dmg(i,5) = fstar(i)
731 uvar(i,1) = sigm(i)
732c..... DEVIATORIC STRESS
733 ENDIF
734 IF(off(i)==one.AND.fstar(i)>=ff) off(i)=four_over_5
735 ENDIF
736 ENDDO
737 ELSEIF(ifla==2)THEN
738C----
739CIFLA set to 2
740C-----
741 DO i=1,nel0
742 IF(off(i)==one)THEN
743C
744 IF(f(i)<=fc)THEN
745 fstar(i) = f(i)
746 df = one
747 ELSE
748 df =(fu - fc)/(ff - fc)
749 fstar(i) = fc + df*(f(i)-fc)
750 ENDIF
751C
752 vm = (signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i)
753 . + three*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2))
754 pn(i) = (signxx(i)+signyy(i))*third
755C
756 sigm1 = one/max(em20,sigm(i))
757 sigm2 = sigm1**2
758 var = three_half * q2 * pn(i) * sigm1
759 var = exp(var)
760 coh = half *(var + one/max(em20,var))
761 sih = half *(var - one/max(em20,var))
762C
763 IF(pn(i)<=zero) THEN
764 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
765 ELSE
766 va= one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
767 ENDIF
768C
769 lamda(i) = zero
770 yldc = sigm(i)
771 crit = vm * sigm2 - va
772C
773 yld( i ) = sigm(i) * sqrt(max(zero,va))
774 IF(crit < zero .OR. yld(i) == zero)THEN
775 lamda(i)=zero
776 ELSE
777C.. stockage des contraintes ....
778 sigxx(i) = signxx(i)
779 sigyy(i) = signyy(i)
780 sigxy(i) = signxy(i)
781 sigzx(i) = signzx(i)
782 sigyz(i) = signyz(i)
783C
784 dp11 = zero
785 dp22 = zero
786 dp33 = zero
787 dp12 = zero
788 dp13 = zero
789 dp23 = zero
790C
791 dtinv =timestep(i)/max(timestep(i)**2,em20)
792C
793 IF(f(i)==one)THEN
794 a21 = ep20
795 ELSE
796 a21 = sigm1 /(one - f(i))
797 ENDIF
798 a11 = fn*exp(-half*((pla(i)-epsn)/sn)**2)/sn
799 a11 = a11*sqr22
800 IF(en/=one) THEN
801 DO iter=1,5
802 IF(pn(i)<=zero)THEN
803 d11 = (two*signxx(i)-signyy(i))*sigm2
804 d22 = (two*signyy(i)-signxx(i))*sigm2
805 d33 = ( - signxx(i) - signyy(i))*sigm2
806 d12 = six*signxy(i)*sigm2
807 d13 = six*signzx(i)*sigm2
808 d23 = six*signyz(i)*sigm2
809 dcrm = -two*vm*sigm1*sigm2
810 dcrf = two*q1*df - two*q3*df*fstar(i)
811 ELSE
812 va2 = q1 * q2 * fstar(i) * sih * sigm1
813 d11 = (two*signxx(i) - signyy(i))*sigm2 + va2
814 d22 = (two*signyy(i) - signxx(i))*sigm2 + va2
815 d33 = ( - signxx(i) - signyy(i))*sigm2 + va2
816 d12 = six*signxy(i)*sigm2
817 d13 = six*signzx(i)*sigm2
818 d23 = six*signyz(i)*sigm2
819 dcrm = -two*vm*sigm1*sigm2 - three*va2*pn(i)*sigm1
820 dcrf = two * q1 * df * coh - two*q3*df*fstar(i)
821 ENDIF
822C
823 a22 = d11*signxx(i) + d22*signyy(i) +
824 . three*(d12*signxy(i)+ d13*signzx(i)+ d23*signyz(i))
825 a22 = a22*a21
826C
827 dsepp = et * en * pla(i) ** (en - one)*
828 . (one + ( epsd(i)*csd)**visp)
829C
830 dcd = c1*(d11**2 + d22**2 + d33**2) +
831 . two*c2*(d11*d22 + d11*d33 + d22*d33) +
832 . two*g*d12**2 + two*gs(i)*( d13**2 + d23**2)
833C
834 lam1= dcd - dcrm * dsepp*a22 -
835 . dcrf*((1-f(i))*(d11 + d22 +d33) + a11*a22)
836 IF(lam1/=zero) lamda(i) = max(zero,(vm*sigm2 - va))/lam1
837C.... PLASTIC DEFORMATION
838 dp11 = dp11 + lamda(i) * d11
839 dp22 = dp22 + lamda(i) * d22
840 dp33 = dp33 + lamda(i) * d33
841 dp12 = dp12 + lamda(i) * d12
842 dp13 = dp13 + lamda(i) * d13
843 dp23 = dp23 + lamda(i) * d23
844C NEW STRESS
845 signxx(i) = sigxx(i) - a1*dp11 - a2*dp22
846 signyy(i) = sigyy(i) - a1*dp22 - a2*dp11
847 signxy(i) = sigxy(i) - two*g*dp12
848 signzx(i)= sigzx(i) - two*gs(i)*dp13
849 signyz(i)= sigyz(i) - two*gs(i)*dp23
850C....... EFFECTIVE PLASTIC DEFORMATION.
851 epsp1(i)= signxx(i)*dp11 + signyy(i)*dp22 +
852 . two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
853 epsp1(i) = epsp1(i)*a21
854 epsp1(i) = max(zero, epsp1(i))
855C
856 vm= signxx(i)**2 + signyy(i)**2
857 . + three*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
858 . - signxx(i)*signyy(i)
859C
860 pn(i) = (signxx(i) + signyy(i))*third
861 var = three_half * q2*pn(i)*sigm1
862 var= exp(var)
863 coh = half * (var + one/max(em20,var))
864 sih = half * (var - one/max(em20,var))
865C
866 IF(pn(i)<=zero) THEN
867 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
868 ELSE
869 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
870 ENDIF
871 ENDDO ! it
872C
873 ELSE
874 DO iter=1,5
875 IF(pn(i)<=zero)THEN
876 d11 = (two*signxx(i)-signyy(i))*sigm2
877 d22 = (two*signyy(i)-signxx(i))*sigm2
878 d33 = ( - signxx(i) - signyy(i))*sigm2
879 d12 = six*signxy(i)*sigm2
880 d13 = six*signzx(i)*sigm2
881 d23 = six*signyz(i)*sigm2
882 dcrm = -two*vm*sigm1*sigm2
883 dcrf = two*q1*df - two*q3*df*fstar(i)
884 ELSE
885 va2 = q1 * q2 * fstar(i) * sih * sigm1
886 d11 = (two*signxx(i) - signyy(i))*sigm2 + va2
887 d22 = (two*signyy(i) - signxx(i))*sigm2 + va2
888 d33 = ( - signxx(i) - signyy(i))*sigm2 + va2
889 d12 = six*signxy(i)*sigm2
890 d13 = six*signzx(i)*sigm2
891 d23 = six*signyz(i)*sigm2
892 dcrm = -two*vm*sigm1*sigm2 - three*va2*pn(i)*sigm1
893 dcrf = two * q1 * df * coh - two*q3*df*fstar(i)
894 ENDIF
895C
896 a22 = d11*signxx(i) + d22*signyy(i) +
897 . three*(d12*signxy(i)+ d13*signzx(i)+ d23*signyz(i))
898 a22 = a22*a21
899C
900 dsepp = et*(1+ ( epsd(i)*csd)**visp)
901 dcd = c1*(d11**2 + d22**2 + d33**2) +
902 . two*c2*(d11*d22 + d11*d33 + d22*d33) +
903 . two*g*d12**2 + two*gs(i)*( d13**2 + d23**2)
904C
905 lam1= dcd - dcrm * dsepp*a22 -
906 . dcrf*((one -f(i))*(d11 + d22 + d33) + a11*a22)
907 IF(lam1/=zero)lamda(i) = max(zero,(vm*sigm2 - va))/lam1
908C.... PLASTIC DEFORMATION
909 dp11 = dp11 + lamda(i) * d11
910 dp22 = dp22 + lamda(i) * d22
911 dp33 = dp33 + lamda(i) * d33
912 dp12 = dp12 + lamda(i) * d12
913 dp13 = dp13 + lamda(i) * d13
914 dp23 = dp23 + lamda(i) * d23
915C NEW STRESS
916 signxx(i) = sigxx(i) - a1*dp11 - a2*dp22
917 signyy(i) = sigyy(i) - a1*dp22 - a2*dp11
918 signxy(i) = sigxy(i) - two*g*dp12
919 signzx(i) = sigzx(i) - two*gs(i)*dp13
920 signyz(i) = sigyz(i) - two*gs(i)*dp23
921C....... EFFECTIVE PLASTIC DEFORMATION.
922 epsp1(i)= signxx(i)*dp11 + signyy(i)*dp22 +
923 . two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
924 epsp1(i) = epsp1(i)*a21
925C
926 epsp1(i) =max(zero, epsp1(i))
927 vm= signxx(i)**2 + signyy(i)**2
928 . + three*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
929 . - signxx(i)*signyy(i)
930C
931 pn(i) = (signxx(i) + signyy(i))*third
932 var = three_half * q2*pn(i)*sigm1
933 var = exp(var)
934 coh = half * (var + one/max(em20,var))
935 sih = half * (var - one/max(em20,var))
936C
937 IF(pn(i)<=zero) THEN
938 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
939 ELSE
940 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
941 ENDIF
942 ENDDO ! it
943 ENDIF
944 yld(i) = sigm(i) * sqrt(max(zero, va))
945 etse( i) = dsepp/(dsepp+e)
946 dezz = nn1*(dp11 +dp22) + dp33
947 thk(i) = thk(i) + dezz*thkly(i)*off(i)
948C
949 fg(i) = fg(i) + (one - f(i))*(dp11 + dp22 + dp33)
950 fn1(i) = zero
951 pn(i) = (signxx(i) + signyy(i))*third
952 IF (pn(i)>=zero) fn1(i) = fn1(i) + a11*epsp1(i)
953 pla(i) = pla(i) + epsp1(i)
954 f(i) = fi + fg(i) + fn1(i)
955 IF(q1==zero.and.
956 . q2==zero.and.
957 . q3==zero)f(i)=zero
958 IF(f(i)<zero)f(i)=zero
959 sigm(i)= (yeild0 + et*pla(i)**en)*(one+(epsd(i)*csd)**visp)
960C....COMPUTE F*.....
961 IF(f(i)<=fc)THEN
962 fstar(i) = f(i)
963 ELSE
964 IF(ff==fc)THEN
965 fstar(i) = ep20
966 ELSE
967 fstar(i) = fc+ (fu-fc)*(f(i)-fc)/(ff-fc)
968 ENDIF
969 ENDIF
970C
971 dmg(i,1) = f(i)/ff
972 dmg(i,2) = fg(i)
973 dmg(i,3) = fn1(i)
974 dmg(i,4) = f(i)
975 dmg(i,5) = fstar(i)
976 uvar(i,1) = sigm(i)
977c..... DEVIATORIC STRESS
978 ENDIF
979C
980 IF(off(i)==one.AND.fstar(i)>=ff) off(i)=four_over_5
981 ENDIF
982C
983 ENDDO
984 ELSE
985C---------
986C iflag = 3
987c----------
988 DO i=1,nel0
989 IF(off(i)==one)THEN
990C
991 IF(f(i)<=fc)THEN
992 fstar(i) = f(i)
993 df = one
994 ELSE
995 df = (fu-fc)/(ff-fc)
996 fstar(i) = fc + df*(f(i)-fc)
997 ENDIF
998 vm = (signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i)
999 . + three * (signxy(i)**2 + signzx(i)**2 + signyz(i)**2))
1000 vm = sqrt(vm)
1001 pn(i) = (signxx(i)+signyy(i))*third
1002C
1003 sigm1 = one/max(em20, sigm(i))
1004 var = three_half * q2 * pn(i) * sigm1
1005 var = exp(var)
1006 coh = half * (var + one/max(em20,var))
1007 sih = half * (var - one/max(em20,var))
1008C
1009 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
1010 va = sqrt(max(zero,va))
1011 yldc = sigm(i) * va
1012C
1013 yld(i ) = yldc
1014 crit = vm - yldc
1015C
1016 IF(vm < yldc .OR. yldc == zero)THEN
1017 lamda(i)=zero
1018 ELSE
1019C.. stockage des contraintes ....
1020 lamda(i) = zero
1021 sigxx(i) = signxx(i)
1022 sigyy(i) = signyy(i)
1023 sigxy(i) = signxy(i)
1024 sigzx(i) = signzx(i)
1025 sigyz(i) = signyz(i)
1026C
1027 dp11 = zero
1028 dp22 = zero
1029 dp33 = zero
1030 dp12 = zero
1031 dp13 = zero
1032 dp23 = zero
1033 dtinv =timestep(i)/max(timestep(i)**2,em20)
1034C
1035 IF(f(i)==one)THEN
1036 a21 = ep20
1037 ELSE
1038 a21 = sigm1 /(one - f(i))
1039 ENDIF
1040 a11 = fn*exp(-half*((pla(i)-epsn)/sn)**2)/sn
1041 a11 = a11*sqr22
1042C
1043 IF(en/=1)THEN
1044 DO iter=1,5
1045 vm1 = one/max(em20, vm)
1046 va1 = one/max(em20, va)
1047 va2 = half * q1*q2*fstar(i)*sih*va1
1048 d11 = half * (two*signxx(i)-signyy(i))*vm1 + va2
1049 d22 = half * (two*signyy(i)-signxx(i))*vm1 + va2
1050 d33 = half* (-signxx(i)-signyy(i))*vm1 + va2
1051 d12 = three * signxy(i)*vm1
1052 d13 = three * signzx(i)*vm1
1053 d23 = three * signyz(i)*vm1
1054 a22 = d11 * signxx(i) + d22 * signyy(i) +
1055 . two*(d12 * signxy(i) + d13 * signzx(i) +
1056 . d23 * signyz(i))
1057C
1058 a22 = a22 * a21
1059 dcrf = - sigm(i)*(q3*fstar(i)*df - q1*coh*df)*va1
1060 dcrm = - va - three*va2*pn(i)*sigm1
1061 dsepp = et * en * pla(i) ** (en - 1)*
1062 . (one+ ( epsd(i)*csd)**visp)
1063 dcd = c1*(d11**2 + d22**2 + d33**2) +
1064 . two*c2*(d11*d22 + d11*d33 + d22*d33) +
1065 . two*g * d12**2 + two*gs(i) * ( d13**2 + d23**2 )
1066C
1067 lam1= dcd - dcrm * dsepp * a22 -
1068 . dcrf *((one - f(i)) * (d11 + d22 + d33) + a11*a22)
1069
1070 IF(lam1/=zero)lamda(i) = max(zero,( vm - yldc )) / lam1
1071C PLASTIC DEFORMATION
1072 dp11 = dp11 + lamda(i) * d11
1073 dp22 = dp22 + lamda(i) * d22
1074 dp33 = dp33 + lamda(i) * d33
1075 dp12 = dp12 + lamda(i) * d12
1076 dp13 = dp13 + lamda(i) * d13
1077 dp23 = dp23 + lamda(i) * d23
1078C NEW STRESS
1079 signxx(i)= sigxx(i) - a1 * dp11 - a2 * dp22
1080 signyy(i)= sigyy(i) - a1 * dp22 - a2 * dp11
1081 signxy(i)= sigxy(i) - two*g * dp12
1082 signzx(i)= sigzx(i) - two*gs(i) * dp13
1083 signyz(i)= sigyz(i) - two*gs(i) * dp23
1084C....... EFFECTIVE PLASTIC DEFORMATION.
1085 epsp1(i)= signxx(i) * dp11 + signyy(i) * dp22 +
1086 . two*(signxy(i) * dp12 + signzx(i) * dp13 +
1087 . signyz(i) * dp23 )
1088 epsp1(i)= epsp1(i)* a21
1089 epsp1(i) = max(zero, epsp1(i))
1090 vm = signxx(i)**2 + signyy(i)**2
1091 . + three*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
1092 . - signxx(i)*signyy(i)
1093 vm = sqrt(vm)
1094C
1095 pn(i) = (signxx(i)+signyy(i))* third
1096 var = three_half * q2 * pn(i) * sigm1
1097 var = exp(var)
1098 coh = half * ( var + one/max(em20,var))
1099 sih = half * ( var - one/max(em20,var))
1100 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
1101 va= sqrt(max(zero,va))
1102 yldc = sigm(i) * va
1103C IF(CONV<1E-6) GOTO 100
1104
1105 ENDDO ! ite
1106 ELSE
1107 DO iter=1,5
1108 vm1 = one/max(em20, vm)
1109 va1 = one/max(em20, va)
1110 va2 = half * q1*q2*fstar(i)*sih*va1
1111 d11 = half * (two*signxx(i)-signyy(i))*vm1 + va2
1112 d22 = half * (two*signyy(i)-signxx(i))*vm1 + va2
1113 d33 = half*( - signxx(i)-signyy(i))*vm1 + va2
1114 d12 = three * signxy(i) * vm1
1115 d13 = three * signzx(i) * vm1
1116 d23 = three * signyz(i) * vm1
1117 a22 = d11 * signxx(i) + d22 * signyy(i) +
1118 . two *(d12 * signxy(i) + d13 * signzx(i) +
1119 . d23 * signyz(i))
1120C
1121 a22 = a22 * a21
1122 dcrf = - sigm(i)*(q3*fstar(i)*df - q1*coh*df)*va1
1123 dcrm = - va - three*va2*pn(i)*sigm1
1124 dsepp = et*(one + ( epsd(i)*csd)**visp)
1125 dcd = c1*(d11**2 + d22**2 + d33**2) +
1126 . two*c2*(d11*d22 + d11*d33 + d22*d33) +
1127 . two*g * d12**2 + two *gs(i) * ( d13**2 + d23**2)
1128C
1129 lam1= dcd - dcrm * dsepp * a22 -
1130 . dcrf *((one -f(i)) * (d11 + d22+d33) + a11*a22)
1131 IF(lam1/=zero) lamda(i) = max(zero,( vm - yldc )) / lam1
1132
1133C PLASTIC DEFORMATION
1134 dp11 = dp11 + lamda(i) * d11
1135 dp22 = dp22 + lamda(i) * d22
1136 dp33 = dp33 + lamda(i) * d33
1137 dp12 = dp12 + lamda(i) * d12
1138 dp13 = dp13 + lamda(i) * d13
1139 dp23 = dp23 + lamda(i) * d23
1140C NEW STRESS
1141 signxx(i)= sigxx(i) - a1 * dp11 - a2 * dp22
1142 signyy(i)= sigyy(i) - a1 * dp22 - a2 * dp11
1143 signxy(i)= sigxy(i) - 2*g * dp12
1144 signzx(i)= sigzx(i) - 2*gs(i) * dp13
1145 signyz(i)= sigyz(i) - 2*gs(i) * dp23
1146C....... EFFECTIVE PLASTIC DEFORMATION.
1147 epsp1(i)= signxx(i) * dp11 + signyy(i) * dp22 +
1148 . two*(signxy(i) * dp12 + signzx(i) * dp13 +
1149 . signyz(i) * dp23 )
1150 epsp1(i)= epsp1(i)*a21
1151 epsp1(i) = max(zero, epsp1(i))
1152 vm = signxx(i)**2 + signyy(i)**2
1153 . + three*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
1154 . - signxx(i)*signyy(i)
1155 vm = sqrt(vm)
1156C
1157 pn(i) = (signxx(i)+signyy(i))* third
1158 var = three_half * q2 * pn(i) * sigm1
1159 var = exp(var)
1160 coh = half * ( var + one/max(em20,var))
1161 sih = half * ( var - one/max(em20,var))
1162 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
1163 va = sqrt(max(zero,va))
1164 yldc= sigm(i) * va
1165C IF(CONV<1E-6) GOTO 100
1166 ENDDO ! ite
1167 ENDIF
1168C
1169 yld(i ) = yldc
1170 etse( i) = dsepp/(dsepp+e)
1171 dezz = nn1*(dp11 +dp22) + dp33
1172 thk(i) = thk(i) + dezz*thkly(i)*off(i)
1173C
1174 pn(i) = (sigxx(i) + sigyy(i))*third
1175 fg(i) = fg(i) + (one - f(i)) * (dp11 + dp22 +dp33)
1176 fn1(i) = zero
1177 IF (pn(i)>=zero) fn1(i) = fn1(i) + a11*epsp1(i)
1178 pla(i) = pla(i) + epsp1(i)
1179 f(i) = fi + fg(i) + fn1(i)
1180 IF(q1==zero.and.
1181 . q2==zero.and.
1182 . q3==zero)f(i)=zero
1183 IF(f(i)<zero)f(i)=zero
1184 sigm(i)= (yeild0 + et*pla(i)**en)* (one+(epsd(i)*csd)**visp)
1185C....COMPUTE F*.....
1186 IF(f(i)<=fc)THEN
1187 fstar(i)=f(i)
1188 ELSE
1189 IF(ff==fc)THEN
1190 fstar(i) = ep20
1191 ELSE
1192 fstar(i) = fc+ (fu-fc)*(f(i)-fc)/(ff-fc)
1193 ENDIF
1194 ENDIF
1195 dmg(i,1) = f(i)/ff
1196 dmg(i,2) = fg(i)
1197 dmg(i,3) = fn1(i)
1198 dmg(i,4) = f(i)
1199 dmg(i,5) = fstar(i)
1200 uvar(i,1) = sigm(i)
1201c..... DEVIATORIC STRESS
1202 ENDIF
1203 IF(off(i)==one.AND.fstar(i)>=ff) off(i)=four_over_5
1204 ENDIF
1205 ENDDO
1206C
1207 ENDIF
1208C
1209C tabulated law
1210 IF(iyeild_tab > 0) THEN
1211 iadbuf = ipm(7,mat(1))-1
1212 xfac = uparam(iadbuf + 22 )
1213 yfac = uparam(iadbuf + 23 )
1214 DO i=1,nel0
1215 ndim= table(itab(1))%NDIM
1216 IF(ndim==2)THEN
1217 nxk=SIZE(table(itab(1))%X(2)%VALUES)
1218 DO j=2,nxk
1219 dx2 = table(itab(1))%X(2)%VALUES(j)*xfac- epsd(i)
1220 IF(dx2 >= zero .OR. j == nxk)THEN
1221 ipos=j-1
1222 r=(table(itab(1))%X(2)%VALUES(j)*xfac-epsd(i))/
1223 . (table(itab(1))%X(2)%VALUES(j)*xfac
1224 . -table(itab(1))%X(2)%VALUES(j-1)*xfac)
1225 EXIT
1226 ENDIF
1227 END DO ! NXK
1228 ELSE
1229 r = one
1230 ENDIF ! NDIM
1231C interpolation on the curve ---
1232 xx(1) = pla(i)
1233 xx(2) = epsd(i)*xfac
1234 CALL table_interp_law76(table(itab(1)),ipos,xx,
1235 . r,dft,yy)
1236 yld(i) = yfac*yy
1237 uvar(i,1) = yld(i)
1238 ENDDO ! NEL0
1239 ENDIF ! IYEILD_TAB
1240C-----------
1241 RETURN
1242 END
interface_table_mod::table_interp_law76
Definition table_mod.F:272
interface_table_mod::table_interp
Definition table_mod.F:242
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
sigeps52c
subroutine sigeps52c(nel0, nuparam, nuvar, mfunc, kfunc, npf, npt0, ipt, iflag, asrate, tf, time, timestep, uparam, rho0, area, eint, thkly, l_dmg, dmg, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thk, pla, uvar, off, ngl, ipm, mat, etse, gs, yld, epsd_pg, table, epsd)
Definition sigeps52c.F:49
table_mod::ttable
Definition table_mod.F:125