OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps86c.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!|| sigeps86c ../engine/source/materials/mat/mat086/sigeps86c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| vinter ../engine/source/tools/curve/vinter.F
30!||====================================================================
31 SUBROUTINE sigeps86c(
32 1 NEL0 ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,
33 2 NPF ,NPT0 ,IPT ,IFLAG ,
34 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,
35 3 AREA ,EINT ,THKLY ,ISRATE ,ASRATE ,
36 4 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
37 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 6 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
39 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 9 SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
42 A SOUNDSP,VISCMAX,THK ,PLA ,UVAR ,
43 B OFF ,NGL ,IPM ,MAT ,ETSE ,
44 C GS ,YLD ,EPSD_PG ,EPSP ,DPLA_I)
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C G l o b a l P a r a m e t e r s
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53C---------+---------+---+---+--------------------------------------------
54C VAR | SIZE |TYP| RW| DEFINITION
55C---------+---------+---+---+--------------------------------------------
56C NEL0 | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL0
57C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
58C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
59C---------+---------+---+---+--------------------------------------------
60C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
61C IFUNC | NFUNC | I | R | FUNCTION INDEX
62C NPF | * | I | R | FUNCTION ARRAY
63C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
64C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
65C IFLAG | * | I | R | GEOMETRICAL FLAGS
66C TF | * | F | R | FUNCTION ARRAY
67C---------+---------+---+---+--------------------------------------------
68C TIME | 1 | F | R | CURRENT TIME
69C TIMESTEP| 1 | F | R | CURRENT TIME STEP
70C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
71C RHO0 | NEL0 | F | R | INITIAL DENSITY
72C AREA | NEL0 | F | R | AREA
73C EINT | 2*NEL0 | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
74C THKLY | NEL0 | F | R | LAYER THICKNESS
75C EPSPXX | NEL0 | F | R | STRAIN RATE XX
76C EPSPYY | NEL0 | F | R | STRAIN RATE YY
77C ... | | | |
78C DEPSXX | NEL0 | F | R | STRAIN INCREMENT XX
79C DEPSYY | NEL0 | F | R | STRAIN INCREMENT YY
80C ... | | | |
81C EPSXX | NEL0 | F | R | STRAIN XX
82C EPSYY | NEL0 | F | R | STRAIN YY
83C ... | | | |
84C SIGOXX | NEL0 | F | R | OLD ELASTO PLASTIC STRESS XX
85C SIGOYY | NEL0 | F | R | OLD ELASTO PLASTIC STRESS YY
86C ... | | | |
87C---------+---------+---+---+--------------------------------------------
88C SIGNXX | NEL0 | F | W | NEW ELASTO PLASTIC STRESS XX
89C SIGNYY | NEL0 | F | W | NEW ELASTO PLASTIC STRESS YY
90C ... | | | |
91C SIGVXX | NEL0 | F | W | VISCOUS STRESS XX
92C SIGVYY | NEL0 | F | W | VISCOUS STRESS YY
93C ... | | | |
94C SOUNDSP | NEL0 | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
95C VISCMAX | NEL0 | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
96C---------+---------+---+---+--------------------------------------------
97C THK | NEL0 | F |R/W| THICKNESS
98C PLA | NEL0 | F |R/W| PLASTIC STRAIN
99C UVAR |NEL0*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
100C OFF | NEL0 | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
101C---------+---------+---+---+--------------------------------------------
102C C o m m o n B l o c k s
103C-----------------------------------------------
104#include "param_c.inc"
105#include "com01_c.inc"
106C-----------------------------------------------
107C I N P U T A r g u m e n t s
108C-----------------------------------------------
109C
110 INTEGER NEL0, NUPARAM, NUVAR, NPT0, IPT,IFLAG(*),
111 . NGL(NEL0),MAT(NEL0),ISRATE,IPM(NPROPMI,*)
112 my_real ,INTENT(IN) :: ASRATE
113 my_real ,DIMENSION(NEL0) ,INTENT(IN) :: EPSD_PG !< global element strain rate in Gauss pt
114 my_real ,DIMENSION(NEL0) ,INTENT(INOUT) :: EPSP !< local strain rate used in equations
115 my_real
116 . TIME,TIMESTEP,
117 . AREA(NEL0),RHO0(NEL0),EINT(NEL0,2),
118 . THKLY(NEL0),PLA(NEL0),
119 . EPSPXX(NEL0),EPSPYY(NEL0),
120 . EPSPXY(NEL0),EPSPYZ(NEL0),EPSPZX(NEL0),
121 . DEPSXX(NEL0),DEPSYY(NEL0),
122 . DEPSXY(NEL0),DEPSYZ(NEL0),DEPSZX(NEL0),
123 . epsxx(nel0) ,epsyy(nel0) ,
124 . epsxy(nel0) ,epsyz(nel0) ,epszx(nel0) ,
125 . sigoxx(nel0),sigoyy(nel0),
126 . sigoxy(nel0),sigoyz(nel0),sigozx(nel0),gs(*)
127 my_real ,DIMENSION(NUPARAM) :: uparam
128C-----------------------------------------------
129C O U T P U T A r g u m e n t s
130C-----------------------------------------------
131 my_real
132 . signxx(nel0),signyy(nel0),
133 . signxy(nel0),signyz(nel0),signzx(nel0),
134 . sigvxx(nel0),sigvyy(nel0),
135 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
136 . soundsp(nel0),viscmax(nel0),etse(nel0),
137 . dpla_i(nel0)
138C-----------------------------------------------
139C I N P U T O U T P U T A r g u m e n t s
140C-----------------------------------------------
141 my_real
142 . uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
143C-----------------------------------------------
144C VARIABLES FOR FUNCTION INTERPOLATION
145C-----------------------------------------------
146 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
147 my_real
148 . FINTER ,TF(*)
149 EXTERNAL FINTER
150C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
151C Y : y = f(x)
152C X : x
153C DYDX : f'(x) = dy/dx
154C IFUNC(J): FUNCTION INDEX
155C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
156C NPF,TF : FUNCTION PARAMETER
157C-----------------------------------------------
158C L o c a l V a r i a b l e s
159C-----------------------------------------------
160 INTEGER I,J,NRATE(MVSIZ),J1,J2,N,NINDX,NMAX,
161 . IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
162 . iad2(mvsiz),ipos2(mvsiz),ilen2(mvsiz),
163 . jj(mvsiz),index(mvsiz),nratem,
164 . nrate1
165 my_real
166 . e(mvsiz),a1(mvsiz),a2(mvsiz),g(mvsiz),
167 . dydx1(mvsiz),dydx2(mvsiz),rate(mvsiz),svm(mvsiz),
168 . y1(mvsiz),y2(mvsiz),dr(mvsiz),
169 . yfac(mvsiz,2),
170 . aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
171 . pp(mvsiz),qq(mvsiz),fail(mvsiz),svmo(mvsiz),h(mvsiz),
172 . epsmax(mvsiz),epsr1(mvsiz),epsr2(mvsiz),fisokin(mvsiz),
173 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz),
174 . hk(mvsiz),nu(mvsiz),nnu1(mvsiz),epsf(mvsiz),
175 . svm2(mvsiz),yld2(mvsiz)
176 my_real
177 . r,umr,a,b,c,amu,s11,s22,s12,p,p2,fac,dezz,
178 . sigz,s1,s2,s3,dpla,vm2,epst(mvsiz),nnu2,
179 . err,f,df,pla_i,q2,yld_i,sigpxx,sigpyy,sigpxy,
180 . alpha,
181 . e1,a11,a21,g1,nnu11,nu11,nu21,nu31,nu41,nu51,nu61,
182 . epsmax1,epsr11,epsr21,fisokin1,g31,
183 . dsxx,dsyy,dsxy,dexx,deyy,dexy,nux, epsf1, mepsf
184 my_real
185 . y11(mvsiz),y21(mvsiz), me, ma1, ma2, mg, mnu,
186 . mepsmax,mepsr1,mepsr2,mfisokin
187 INTEGER JST(MVSIZ+1), IC, MNRATE
188C
189 DATA nmax/3/
190C-----------------------------------------------
191C USER VARIABLES INITIALIZATION
192C-----------------------------------------------
193 e1 = uparam(2)
194 a11 = uparam(3)
195 a21 = uparam(4)
196 g1 = uparam(5)
197 g31 = three*g1
198 nux = uparam(6)
199 nnu11 = nux / (1. - nux)
200 nrate1 = nint(uparam(1))
201 epsmax1 =uparam(7+2*nrate1)
202 epsr11 =uparam(8+2*nrate1)
203 epsr21 =uparam(9+2*nrate1)
204 fisokin1=uparam(14+2*nrate1)
205 epsf1= uparam(15+2*nrate1)
206 IF(epsmax1==zero)THEN
207 IF(tf(npf(ipm(11,mat(1))+1)-1)==zero)THEN
208 epsmax1=tf(npf(ipm(11,mat(1))+1)-2)
209 ELSE
210 epsmax1= ep30
211 ENDIF
212 ENDIF
213C
214 IF (isigi==0) THEN
215 IF(time==zero)THEN
216 DO i=1,nel0
217 uvar(i,1)=zero
218 uvar(i,2)=zero
219 uvar(i,3)=zero
220 uvar(i,4)=zero
221 DO j=1,nrate1
222 uvar(i,j+4)=zero
223 ENDDO
224 ENDDO
225 ENDIF
226 ENDIF
227C------------------------------------------
228C ECROUISSAGE CINE
229C------------------------------------------
230 DO i=1,nel0
231! SIGOXX(I) = SIGOXX(I) - UVAR(I,2)
232! SIGOYY(I) = SIGOYY(I) - UVAR(I,3)
233! SIGOXY(I) = SIGOXY(I) - UVAR(I,4)
234C
235 signxx(i)=sigoxx(i) - uvar(i,2) +a11*depsxx(i)+a21*depsyy(i)
236 signyy(i)=sigoyy(i) - uvar(i,3) +a21*depsxx(i)+a11*depsyy(i)
237 signxy(i)=sigoxy(i) - uvar(i,4) +g1 *depsxy(i)
238 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
239 signzx(i)=sigozx(i)+gs(i) *depszx(i)
240 sigexx(i) = signxx(i)
241 sigeyy(i) = signyy(i)
242 sigexy(i) = signxy(i)
243C
244 soundsp(i) = sqrt(a11/rho0(i))
245 viscmax(i) = zero
246 etse(i) = one
247C-------------------
248C STRAIN RATE
249C-------------------
250C
251 IF (israte == 0) THEN
252 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
253 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
254 . + epspxy(i)*epspxy(i) ) )
255 ELSE
256 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
257 END IF
258C-------------------
259C STRAIN
260C-------------------
261C
262 epst(i) = half*( epsxx(i)+epsyy(i)
263 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
264 . + epsxy(i)*epsxy(i) ) )
265 jj(i) = 1
266 ENDDO
267C-------------------
268C CRITERE
269C-------------------
270C inversion boucles
271 DO j=2,nrate1-1
272 DO i=1,nel0
273 IF(epsp(i)>=uparam(6+j)) jj(i) = j
274 ENDDO
275 ENDDO
276
277 DO i=1,nel0
278 fac=uparam(6+jj(i))
279 rate(i)=(epsp(i) - fac)/(uparam(7+jj(i)) - fac)
280 yfac(i,1)=uparam(6+nrate1+jj(i))
281 yfac(i,2)=uparam(6+nrate1+jj(i)+1)
282 ENDDO
283
284 DO i=1,nel0
285 j1 = jj(i)
286 j2 = j1+1
287 ipos1(i) = nint(uvar(i,j1+4))
288 iad1(i) = npf(ipm(10+j1,mat(1))) / 2 + 1
289 ilen1(i) = npf(ipm(10+j1,mat(1))+1) / 2 - iad1(i)-ipos1(i)
290 ipos2(i) = nint(uvar(i,j2+4))
291 iad2(i) = npf(ipm(10+j2,mat(1))) / 2 + 1
292 ilen2(i) = npf(ipm(10+j2,mat(1))+1) / 2 - iad2(i)-ipos2(i)
293 END DO
294C
295 CALL vinter(tf,iad1,ipos1,ilen1,nel0,pla,dydx1,y1)
296 CALL vinter(tf,iad2,ipos2,ilen2,nel0,pla,dydx2,y2)
297C
298 IF (fisokin1==zero) THEN
299 DO i=1,nel0
300 j1 = jj(i)
301 j2 = j1+1
302 y1(i)=y1(i)*yfac(i,1)
303 y2(i)=y2(i)*yfac(i,2)
304 fac = rate(i)
305 yld(i) = y1(i) + fac*(y2(i)-y1(i))
306 yld(i) = max(yld(i),em20)
307 dydx1(i)=dydx1(i)*yfac(i,1)
308 dydx2(i)=dydx2(i)*yfac(i,2)
309 h(i) =dydx1(i) + fac*(dydx2(i)-dydx1(i))
310 uvar(i,j1+4) = ipos1(i)
311 uvar(i,j2+4) = ipos2(i)
312 ENDDO
313 ELSEIF (fisokin1==1.) THEN
314C
315C ECROUISSAGE CINEMATIQUE
316C
317 DO i=1,nel0
318 j1 = jj(i)
319 j2 = j1+1
320 fac = rate(i)
321 dydx1(i)=dydx1(i)*yfac(i,1)
322 dydx2(i)=dydx2(i)*yfac(i,2)
323 h(i) = dydx1(i) + fac*(dydx2(i)-dydx1(i))
324 uvar(i,j1+4) = ipos1(i)
325 uvar(i,j2+4) = ipos2(i)
326C ECROUISSAGE CINEMATIQUE
327 y1(i)=tf(npf(ipm(10+j1,mat(1)))+1)
328 y2(i)=tf(npf(ipm(10+j2,mat(1)))+1)
329 y1(i)=y1(i)*yfac(i,1)
330 y2(i)=y2(i)*yfac(i,2)
331 yld(i) = y1(i) + fac*(y2(i)-y1(i))
332 ENDDO
333
334 ELSE
335C
336C ECROUISSAGE CINEMATIQUE
337C
338 DO i=1,nel0
339 j1 = jj(i)
340 j2 = j1+1
341 y1(i)=y1(i)*yfac(i,1)
342 y2(i)=y2(i)*yfac(i,2)
343 fac = rate(i)
344 yld(i) = y1(i) + fac*(y2(i)-y1(i))
345 yld(i) = max(yld(i),em20)
346 dydx1(i)=dydx1(i)*yfac(i,1)
347 dydx2(i)=dydx2(i)*yfac(i,2)
348 h(i) = dydx1(i) + fac*(dydx2(i)-dydx1(i))
349 uvar(i,j1+4) = ipos1(i)
350 uvar(i,j2+4) = ipos2(i)
351C ECROUISSAGE CINEMATIQUE
352 y1(i)=tf(npf(ipm(10+j1,mat(1)))+1)
353 y2(i)=tf(npf(ipm(10+j2,mat(1)))+1)
354 y1(i)=y1(i)*yfac(i,1)
355 y2(i)=y2(i)*yfac(i,2)
356 yld(i) = (1.-fisokin1) * yld(i) +
357 . fisokin1 * (y1(i) + fac*(y2(i)-y1(i)))
358 ENDDO
359 ENDIF
360
361C-------------------
362C PROJECTION
363C-------------------
364 IF(iflag(1)==0)THEN
365 nu31 = one-nnu11
366C projection radiale
367 DO i=1,nel0
368 svm(i)=sqrt(signxx(i)*signxx(i)
369 . +signyy(i)*signyy(i)
370 . -signxx(i)*signyy(i)
371 . +three*signxy(i)*signxy(i))
372 r = min(one,yld(i)/max(em20,svm(i)))
373 signxx(i)=signxx(i)*r
374 signyy(i)=signyy(i)*r
375 signxy(i)=signxy(i)*r
376 umr = one - r
377 dpla_i(i) = off(i)*svm(i)*umr/(g31+h(i))
378 pla(i) = pla(i) + dpla_i(i)
379 s1=half*(signxx(i)+signyy(i))
380 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
381 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
382 thk(i) = thk(i) + dezz*thkly(i)*off(i)
383 IF(r<one) etse(i)= h(i)/(h(i)+e1)
384 ENDDO
385C
386 ELSEIF(iflag(1)==1)THEN
387C-------------------------
388C CRITERE DE VON MISES
389C-------------------------
390 DO i=1,nel0
391 h(i) = max(zero,h(i))
392 s1=signxx(i)+signyy(i)
393 s2=signxx(i)-signyy(i)
394 s3=signxy(i)
395 aa(i)=fourth*s1*s1
396 bb(i)=three_over_4*s2*s2+3.*s3*s3
397 svm(i)=sqrt(aa(i)+bb(i))
398 dezz = -(depsxx(i)+depsyy(i))*nnu11
399 thk(i) = thk(i) + dezz*thkly(i)*off(i)
400 ENDDO
401C-------------------------
402C GATHER PLASTIC FLOW
403C-------------------------
404 nindx=0
405 DO i=1,nel0
406 IF(svm(i)>yld(i).AND.off(i)==one) THEN
407 nindx=nindx+1
408 index(nindx)=i
409 ENDIF
410 ENDDO
411C
412 IF(nindx/=0) THEN
413C---------------------------
414C DEP EN CONTRAINTE PLANE
415C---------------------------
416 DO j=1,nindx
417 i=index(j)
418 dpla_j(i)=(svm(i)-yld(i))/(g31+h(i))
419 etse(i)= h(i)/(h(i)+e1)
420 hk(i) = h(i)*(one-fisokin1)
421 ENDDO
422C
423 nu11 = one/(one-nux)
424 nu21 = one/(one+nux)
425 nu31 = one-nnu11
426 DO n=1,nmax
427#include "vectorize.inc"
428 DO j=1,nindx
429 i=index(j)
430 dpla_i(i)=dpla_j(i)
431 yld_i =yld(i)+hk(i)*dpla_i(i)
432 dr(i) =half*e1*dpla_i(i)/yld_i
433 pp(i) =one/(one+dr(i)*nu11)
434 qq(i) =one/(one+three*dr(i)*nu21)
435 p2 =pp(i)*pp(i)
436 q2 =qq(i)*qq(i)
437 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
438 df =-(aa(i)*nu11*p2*pp(i)+three*bb(i)*nu21*q2*qq(i))
439 . *(e1-two*dr(i)*hk(i))/yld_i
440 . -two*hk(i)*yld_i
441 df = sign(max(abs(df),em20),df)
442 IF(dpla_i(i)>zero) THEN
443C IF(N==NMAX) THEN
444C ERR=ABS(F/(DF*DPLA_I(I)))
445C WRITE(50,'(2I5,F8.5,4E12.3)') N,I,ERR,DPLA_I(I),F
446C ENDIF
447 dpla_j(i)=max(zero,dpla_i(i)-f/df)
448 ELSE
449 dpla_j(i)=zero
450 ENDIF
451 ENDDO
452 ENDDO
453C------------------------------------------
454C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
455C------------------------------------------
456#include "vectorize.inc"
457 DO j=1,nindx
458 i=index(j)
459 pla(i) = pla(i) + dpla_i(i)
460 s1=(signxx(i)+signyy(i))*pp(i)
461 s2=(signxx(i)-signyy(i))*qq(i)
462 signxx(i)=half*(s1+s2)
463 signyy(i)=half*(s1-s2)
464 signxy(i)=signxy(i)*qq(i)
465 dezz = - nu31*dr(i)*s1/e1
466 thk(i) = thk(i) + dezz*thkly(i)*off(i)
467 ENDDO
468 ENDIF
469C-------------------------------------------
470 ELSEIF(iflag(1)==2)THEN
471C-------------------------
472C CRITERE DE VON MISES
473C-------------------------
474 DO i=1,nel0
475 h(i) = max(zero,h(i))
476 svm2(i)=signxx(i)*signxx(i)
477 . +signyy(i)*signyy(i)
478 . -signxx(i)*signyy(i)
479 . +three*signxy(i)*signxy(i)
480 svm(i)=sqrt(svm2(i))
481 dezz = -(depsxx(i)+depsyy(i))*nnu11
482 thk(i) = thk(i) + dezz*thkly(i)*off(i)
483 ENDDO
484C-------------------------
485C GATHER PLASTIC FLOW
486C-------------------------
487 nindx=0
488 DO i=1,nel0
489 yld2(i)=yld(i)*yld(i)
490 IF(svm2(i)>yld2(i).AND.off(i)==one) THEN
491 nindx=nindx+1
492 index(nindx)=i
493 ENDIF
494 ENDDO
495C
496 IF(nindx/=0) THEN
497C-------------
498C PROJ NORMALE AU CRITERE AVEC CALCUL APPROCHE DE LA NORMALE + RETOUR RADIAL
499C-------------
500 nu31 = one-nnu11
501 DO j=1,nindx
502 i=index(j)
503 a=(svm2(i)-yld2(i))
504 . /(five*svm2(i)+three*(-signxx(i)*signyy(i)+signxy(i)*signxy(i)))
505 s1=(one-two*a)*signxx(i)+ a*signyy(i)
506 s2= a*signxx(i)+(one-two*a)*signyy(i)
507 s3=(one-three*a)*signxy(i)
508 signxx(i)=s1
509 signyy(i)=s2
510 signxy(i)=s3
511 dpla_i(i) = off(i)*(svm(i)-yld(i))/(g31+h(i))
512C
513 hk(i) = h(i)*(one-fisokin1)
514 yld(i) =yld(i)+hk(i)*dpla_i(i)
515 END DO
516C
517 DO j=1,nindx
518 i=index(j)
519 svm(i)=sqrt( signxx(i)*signxx(i)
520 . +signyy(i)*signyy(i)
521 . -signxx(i)*signyy(i)
522 . +three*signxy(i)*signxy(i))
523 r = min(one,yld(i)/max(em20,svm(i)))
524 signxx(i)=signxx(i)*r
525 signyy(i)=signyy(i)*r
526 signxy(i)=signxy(i)*r
527 pla(i) = pla(i) + dpla_i(i)
528 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
529 dezz = -nu31*dezz
530 thk(i) = thk(i) + dezz*thkly(i)*off(i)
531 etse(i)= h(i)/(h(i)+e1)
532 END DO
533 END IF
534 ENDIF
535C
536 DO i=1,nel0
537 IF(off(i)==one)THEN
538 IF(pla(i)>epsmax1.OR.epst(i)>epsf1)off(i)=four_over_5
539 END IF
540 ENDDO
541C------------------------------------------
542C ECROUISSAGE CINE
543C------------------------------------------
544C test bypass
545 IF (fisokin1/=zero) THEN
546 DO i=1,nel0
547 dsxx = sigexx(i) - signxx(i)
548 dsyy = sigeyy(i) - signyy(i)
549 dsxy = sigexy(i) - signxy(i)
550 dexx = (dsxx - nux*dsyy)
551 deyy = (dsyy - nux*dsxx)
552 dexy = two*(one+nux)*dsxy
553 alpha = fisokin1*h(i)/(e1+h(i))/three
554 sigpxx = alpha*(four*dexx+two*deyy)
555 sigpyy = alpha*(four*deyy+two*dexx)
556 sigpxy = alpha*dexy
557C
558 signxx(i) = signxx(i) + uvar(i,2)
559 signyy(i) = signyy(i) + uvar(i,3)
560 signxy(i) = signxy(i) + uvar(i,4)
561 uvar(i,2) = uvar(i,2) + sigpxx
562 uvar(i,3) = uvar(i,3) + sigpyy
563 uvar(i,4) = uvar(i,4) + sigpxy
564 ENDDO
565 ENDIF
566C
567C
568 RETURN
569 END
570C
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps86c
subroutine sigeps86c(nel0, nuparam, nuvar, mfunc, kfunc, npf, npt0, ipt, iflag, tf, time, timestep, uparam, rho0, area, eint, thkly, israte, asrate, 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, epsp, dpla_i)
Definition sigeps86c.F:45
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72