OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps56c.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!|| sigeps56c ../engine/source/materials/mat/mat056/sigeps56c.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 sigeps56c(
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 -- NEC Commonblock saves the state between IPT calls
107 COMMON /nectrxi/ jst ,ic ,ifunc,minlen,nfuncv,
108 . nfuncm,imix ,imixv ,ifuncm,
109 . jend
110 COMMON /nectrxr/ ng_nrates
111C
112 INTEGER MINLEN,GLIMIT
113 PARAMETER (GLIMIT = 64)
114c
115C-----------------------------------------------
116C I N P U T A r g u m e n t s
117C-----------------------------------------------
118C
119 INTEGER NEL0, NUPARAM, NUVAR, NPT0,ISRATE, IPT,IFLAG(*),
120 . ipm(npropmi,*),ngl(nel0),mat(nel0)
121 my_real ,INTENT(IN) :: asrate
122 my_real ,DIMENSION(NEL0) ,INTENT(IN) :: epsd_pg
123 my_real ,DIMENSION(NEL0) ,INTENT(INOUT) :: epsp
124 my_real time,timestep,uparam(*),
125 . area(nel0),rho0(nel0),eint(nel0,2),
126 . thkly(nel0),pla(nel0),
127 . epspxx(nel0),epspyy(nel0),
128 . epspxy(nel0),epspyz(nel0),epspzx(nel0),
129 . depsxx(nel0),depsyy(nel0),
130 . depsxy(nel0),depsyz(nel0),depszx(nel0),
131 . epsxx(nel0) ,epsyy(nel0) ,
132 . epsxy(nel0) ,epsyz(nel0) ,epszx(nel0) ,
133 . sigoxx(nel0),sigoyy(nel0),
134 . sigoxy(nel0),sigoyz(nel0),sigozx(nel0),gs(*)
135C-----------------------------------------------
136C O U T P U T A r g u m e n t s
137C-----------------------------------------------
138 my_real
139 . signxx(nel0),signyy(nel0),
140 . signxy(nel0),signyz(nel0),signzx(nel0),
141 . sigvxx(nel0),sigvyy(nel0),
142 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
143 . soundsp(nel0),viscmax(nel0),etse(nel0),
144 . dpla_i(nel0)
145C-----------------------------------------------
146C I N P U T O U T P U T A r g u m e n t s
147C-----------------------------------------------
148 my_real uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
149C-----------------------------------------------
150C VARIABLES FOR FUNCTION INTERPOLATION
151C-----------------------------------------------
152 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
153 my_real FINTER ,TF(*)
154 EXTERNAL FINTER
155C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
156C Y : y = f(x)
157C X : x
158C DYDX : f'(x) = dy/dx
159C IFUNC(J): FUNCTION INDEX
160C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
161C NPF,TF : FUNCTION PARAMETER
162C-----------------------------------------------
163C L o c a l V a r i a b l e s
164C-----------------------------------------------
165 INTEGER I,J,J1,J2,N,NINDX,NMAX,IADBUF,NFUNC,
166 . iad1(mvsiz),ipos1(mvsiz),ilen1(mvsiz),
167 . iad2(mvsiz),ipos2(mvsiz),ilen2(mvsiz),
168 . jj(mvsiz),index(mvsiz),ifunc(mvsiz,100),nratem,
169 . nrate1,ifunc1(100),nfuncv(mvsiz),
170 . nfuncm,iperf
171 INTEGER JST(MVSIZ+1),IC,MNRATE
172 my_real
173 . MX,MX2,MX3,MX4,ME,MA1,MA2,MG,MNU,
174 . MEPSMAX,MEPSR1,MEPSR2,MFISOKIN
175 my_real
176 . DYDX1(MVSIZ),DYDX2(MVSIZ),RATE(MVSIZ,2),SVM(MVSIZ),
177 . Y1(MVSIZ),Y2(MVSIZ),DR(MVSIZ),
178 . yfac(mvsiz,2),
179 . aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
180 . pp(mvsiz),qq(mvsiz),fail(mvsiz),svmo(mvsiz),h(mvsiz),
181 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz),
182 . hk(mvsiz)
183 my_real
184 . r,umr,a,b,c,amu,s11,s22,s12,p,p2,fac,dezz,
185 . sigz,s1,s2,s3,dpla,vm2,epst,nnu2,
186 . err,f,df,pla_i,q2,yld_i,sigpxx,sigpyy,sigpxy,alpha,
187 . e1,a11,a21,g1,g31,nnu11,nu11,nu21,nu31,nu41,nu51,nu61,
188 . epsmax1,epsr11,epsr21,fisokin1,
189 . dsxx,dsyy,dsxy,dexx,deyy,dexy,nux
190C
191C -- Added variables for NEC grouping
192 my_real
193 . ng_e,ng_a1,ng_a2,ng_g,ng_nu,ng_g3,
194 . ng_espmax,ng_epsr1,ng_epsr2,ng_fisokin,
195 . ng_nrates(mvsiz,100),
196 . ng_rate1,ng_rate2,ng_yfac1,ng_yfac2,
197 . ng_ipos1,ng_ipos2,ng_iad1,ng_iad2,ng_ilen1,ng_ilen2,
198 . ng_eps1,ng_eps2,ng_y11,ng_y12,ng_epsmax,
199 . y11(mvsiz),y12(mvsiz)
200 INTEGER NG_NRATE, NG_CRVVAL,EPSMAX0LIST(MVSIZ)
201 INTEGER
202 . K,L,M,JJS(MVSIZ+1,100+1),JCNT,
203 . ISQRT,ISQRTLIST(MVSIZ),
204 . DIFF(MVSIZ),IMIX,IMIXV(MVSIZ),TEMPIC,JEND(MVSIZ),
205 . ifuncm(mvsiz,100),nfuncmv(mvsiz),jlast(mvsiz+1),
206 . iepsmax0,jjindex(mvsiz),jjstart,jjend
207C -- End NEC grouping
208 DATA nmax/3/,iperf/0/
209C-----------------------------------------------
210C USER VARIABLES INITIALIZATION
211C-----------------------------------------------
212 nfunc = ipm(10,mat(1))
213 DO j=1,nfunc
214 ifunc1(j)=ipm(10+j,mat(1))
215 ENDDO
216C
217 iadbuf = ipm(7,mat(1))-1
218 e1 = uparam(iadbuf+2)
219 a11 = uparam(iadbuf+3)
220 a21 = uparam(iadbuf+4)
221 g1 = uparam(iadbuf+5)
222 g31 = three*g1
223 nux = uparam(iadbuf+6)
224 nrate1 = nint(uparam(iadbuf+1))
225 epsmax1=uparam(iadbuf+6+2*nrate1+1)
226 IF(epsmax1==zero)THEN
227 IF(tf(npf(ifunc1(1)+1)-1)==zero)THEN
228 epsmax1=tf(npf(ifunc1(1)+1)-2)
229 ELSE
230 epsmax1= ep30
231 ENDIF
232 ENDIF
233C
234 nnu11 = nux / (one - nux)
235 nnu2 = nnu11*nnu11
236 nu11 = one/(one-nux)
237 nu21 = one/(one+nux)
238 nu31 = one-nnu11
239 nu41 = one + nnu2 + nnu11
240 nu51 = one + nnu2 - two*nnu11
241 nu61 = half - nnu2 + half*nnu11
242C
243 epsr11 =uparam(iadbuf+6+2*nrate1+2)
244 epsr21 =uparam(iadbuf+6+2*nrate1+3)
245 fisokin1=uparam(iadbuf+6+2*nrate1+8)
246C
247 IF (isigi==0) THEN
248 IF(time==zero)THEN
249 DO i=1,nel0
250 uvar(i,1)=zero
251 uvar(i,2)=zero
252 uvar(i,3)=zero
253 uvar(i,4)=zero
254 DO j=1,nrate1
255 uvar(i,j+4)=zero
256 ENDDO
257 ENDDO
258 ENDIF
259 ENDIF
260C------------------------------------------
261 DO i=1,nel0
262 dpla_i(i) =zero
263 ENDDO
264C-----------------------------------------------
265C
266 DO i=1,nel0
267C
268 signxx(i)=sigoxx(i) - uvar(i,2) +a11*depsxx(i)+a21*depsyy(i)
269 signyy(i)=sigoyy(i) - uvar(i,3) +a21*depsxx(i)+a11*depsyy(i)
270 signxy(i)=sigoxy(i) - uvar(i,4) +g1 *depsxy(i)
271 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
272 signzx(i)=sigozx(i)+gs(i) *depszx(i)
273 sigexx(i) = signxx(i)
274 sigeyy(i) = signyy(i)
275 sigexy(i) = signxy(i)
276C
277 soundsp(i) = sqrt(a11/rho0(i))
278 viscmax(i) = zero
279 etse(i) = one
280C-------------------
281C STRAIN RATE
282C-------------------
283 IF (israte==0) THEN
284 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
285 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
286 . + epspxy(i)*epspxy(i) ) )
287 ELSE
288 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
289 ENDIF
290C-------------------
291C STRAIN
292C-------------------
293 epst = half*( epsxx(i)+epsyy(i)
294 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
295 . + epsxy(i)*epsxy(i) ) )
296 fail(i) = max(zero,min(one,(epsr21-epst)/(epsr21-epsr11)))
297C
298 ENDDO
299C-------------------
300C CRITERE
301C-------------------
302 DO i=1,nel0
303 jj(i) = 1
304 ENDDO
305 iadbuf = ipm(7,mat(1))-1
306C inversion boucles
307 DO j=2,nrate1-1
308 DO i=1,nel0
309 IF(epsp(i)>=uparam(iadbuf+6+j)) jj(i) = j
310 ENDDO
311 ENDDO
312 DO i=1,nel0
313 rate(i,1)=uparam(iadbuf+6+jj(i))
314 rate(i,2)=uparam(iadbuf+6+jj(i)+1)
315 yfac(i,1)=uparam(iadbuf+6+nrate1+jj(i))
316 yfac(i,2)=uparam(iadbuf+6+nrate1+jj(i)+1)
317 ENDDO
318C
319 DO i=1,nel0
320 j1 = jj(i)
321 j2 = j1+1
322 ipos1(i) = nint(uvar(i,j1+4))
323 iad1(i) = npf(ifunc1(j1)) / 2 + 1
324 ilen1(i) = npf(ifunc1(j1)+1) / 2 - iad1(i) - ipos1(i)
325 ipos2(i) = nint(uvar(i,j2+4))
326 iad2(i) = npf(ifunc1(j2)) / 2 + 1
327 ilen2(i) = npf(ifunc1(j2)+1) / 2 - iad2(i) - ipos2(i)
328 ENDDO
329C
330 CALL vinter(tf,iad1,ipos1,ilen1,nel0,pla,dydx1,y1)
331 CALL vinter(tf,iad2,ipos2,ilen2,nel0,pla,dydx2,y2)
332C
333 IF (fisokin1==zero) THEN
334 DO i=1,nel0
335 j1 = jj(i)
336 j2 = j1+1
337 y1(i)=y1(i)*yfac(i,1)
338 y2(i)=y2(i)*yfac(i,2)
339 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
340 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
341 yld(i) = max(yld(i),em20)
342 dydx1(i)=dydx1(i)*yfac(i,1)
343 dydx2(i)=dydx2(i)*yfac(i,2)
344 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
345 uvar(i,j1+4) = ipos1(i)
346 uvar(i,j2+4) = ipos2(i)
347 ENDDO
348 ELSEIF (fisokin1==1.) THEN
349 DO i=1,nel0
350 j1 = jj(i)
351 j2 = j1+1
352 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
353 dydx1(i)=dydx1(i)*yfac(i,1)
354 dydx2(i)=dydx2(i)*yfac(i,2)
355 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
356 uvar(i,j1+4) = ipos1(i)
357 uvar(i,j2+4) = ipos2(i)
358C ECROUISSAGE CINEMATIQUE
359 y1(i)=tf(npf(ifunc1(j1))+1)
360 y2(i)=tf(npf(ifunc1(j2))+1)
361 y1(i)=y1(i)*yfac(i,1)
362 y2(i)=y2(i)*yfac(i,2)
363 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
364 ENDDO
365 ELSE
366 DO i=1,nel0
367 j1 = jj(i)
368 j2 = j1+1
369 y1(i)=y1(i)*yfac(i,1)
370 y2(i)=y2(i)*yfac(i,2)
371 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
372 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
373 yld(i) = max(yld(i),em20)
374 dydx1(i)=dydx1(i)*yfac(i,1)
375 dydx2(i)=dydx2(i)*yfac(i,2)
376 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
377 uvar(i,j1+4) = ipos1(i)
378 uvar(i,j2+4) = ipos2(i)
379C ECROUISSAGE CINEMATIQUE
380 y1(i)=tf(npf(ifunc1(j1))+1)
381 y2(i)=tf(npf(ifunc1(j2))+1)
382 y1(i)=y1(i)*yfac(i,1)
383 y2(i)=y2(i)*yfac(i,2)
384 yld(i) = (1.-fisokin1) * yld(i) +
385 . fisokin1 * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
386 ENDDO
387 ENDIF
388
389C-------------------
390C PROJECTION
391C-------------------
392 IF(iflag(1)==0)THEN
393C projection radiale
394 DO i=1,nel0
395 svm(i)=sqrt(signxx(i)*signxx(i)
396 . +signyy(i)*signyy(i)
397 . -signxx(i)*signyy(i)
398 . +three*signxy(i)*signxy(i))
399 r = min(one,yld(i)/max(em20,svm(i)))
400 signxx(i)=signxx(i)*r
401 signyy(i)=signyy(i)*r
402 signxy(i)=signxy(i)*r
403 umr = one - r
404 dpla_i(i) = off(i)*svm(i)*umr/(g31+h(i))
405 pla(i) = pla(i) + dpla_i(i)
406 s1=half*(signxx(i)+signyy(i))
407 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
408 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
409 thk(i) = thk(i) + dezz*thkly(i)*off(i)
410 IF(r<one) etse(i)= h(i)/(h(i)+e1)
411 ENDDO
412C
413 ELSEIF(iflag(1)==1)THEN
414C-------------------------
415C CRITERE DE VON MISES
416C-------------------------
417 DO i=1,nel0
418 h(i) = max(zero,h(i))
419 s1=signxx(i)+signyy(i)
420 s2=signxx(i)-signyy(i)
421 s3=signxy(i)
422 aa(i)=fourth*s1*s1
423 bb(i)=three_over_4*s2*s2+3.*s3*s3
424 svm(i)=sqrt(aa(i)+bb(i))
425 dezz = -(depsxx(i)+depsyy(i))*nnu11
426 thk(i) = thk(i) + dezz*thkly(i)*off(i)
427 ENDDO
428C-------------------------
429C GATHER PLASTIC FLOW
430C-------------------------
431 nindx=0
432 DO i=1,nel0
433 IF(svm(i)>yld(i).AND.off(i)==one) THEN
434 nindx=nindx+1
435 index(nindx)=i
436 ENDIF
437 ENDDO
438C
439 IF(nindx/=0) THEN
440C---------------------------
441C DEP EN CONTRAINTE PLANE
442C---------------------------
443 DO j=1,nindx
444 i=index(j)
445 dpla_j(i)=(svm(i)-yld(i))/(g31+h(i))
446 etse(i)= h(i)/(h(i)+e1)
447 hk(i) = h(i)*(one-fisokin1)
448 ENDDO
449C
450 DO n=1,nmax
451#include "vectorize.inc"
452 DO j=1,nindx
453 i=index(j)
454 dpla_i(i)=dpla_j(i)
455 yld_i =yld(i)+hk(i)*dpla_i(i)
456 dr(i) =half*e1*dpla_i(i)/yld_i
457 pp(i) =one/(one+dr(i)*nu11)
458 qq(i) =one/(one+three*dr(i)*nu21)
459 p2 =pp(i)*pp(i)
460 q2 =qq(i)*qq(i)
461 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
462 df =-(aa(i)*nu11*p2*pp(i)+three*bb(i)*nu21*q2*qq(i))
463 . *(e1-two*dr(i)*hk(i))/yld_i
464 . -two*hk(i)*yld_i
465 df = sign(max(abs(df),em20),df)
466 IF(dpla_i(i)>zero) THEN
467 dpla_j(i)=max(zero,dpla_i(i)-f/df)
468 ELSE
469 dpla_j(i)=zero
470 ENDIF
471 ENDDO
472 ENDDO
473C------------------------------------------
474C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
475C------------------------------------------
476#include "vectorize.inc"
477 DO j=1,nindx
478 i=index(j)
479 pla(i) = pla(i) + dpla_i(i)
480 s1=(signxx(i)+signyy(i))*pp(i)
481 s2=(signxx(i)-signyy(i))*qq(i)
482 signxx(i)=half*(s1+s2)
483 signyy(i)=half*(s1-s2)
484 signxy(i)=signxy(i)*qq(i)
485 dezz = - nu31*dr(i)*s1/e1
486 thk(i) = thk(i) + dezz*thkly(i)*off(i)
487 ENDDO
488 ENDIF
489C-------------------------------------------
490 ELSEIF(iflag(1)==2)THEN
491C projection radial sur le deviateur sur un critere reduit
492C projection elastique en z => sig33 = 0
493C le coef. de reduction du critere est tel que
494C l'on se trouve sur le critere apres les 2 projections
495 DO i=1,nel0
496 p = -(signxx(i)+signyy(i))*third
497 s11 = signxx(i)+p
498 s22 = signyy(i)+p
499C s33 = p = -(S11 + S22)
500 s12 = signxy(i)
501C
502 p2 = p*p
503 vm2= three*(s12*s12 - s11*s22)
504 a = p2*nu41 + vm2
505 vm2= three*p2 + vm2
506 b = p2*nu61
507 c = p2*nu51 - yld(i)*yld(i)
508 r = min(one,(-b + sqrt(max(zero,b*b-a*c)))/max(a ,em20))
509 signxx(i) = s11*r - p
510 signyy(i) = s22*r - p
511 signxy(i) = s12*r
512C signzz = p*r - p
513C proj. signzz = 0.
514 umr = one - r
515 sigz = nnu11*p*umr
516 signxx(i) = signxx(i) + sigz
517 signyy(i) = signyy(i) + sigz
518 svm(i)=sqrt(vm2)
519 dpla_i(i) = off(i)*svm(i)*umr/(g31+h(i))
520 pla(i) = pla(i) + dpla_i(i)
521 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
522 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
523 thk(i) = thk(i) + dezz*thkly(i)*off(i)
524 IF(r<one) etse(i)= h(i)/(h(i)+e1)
525 ENDDO
526 ENDIF
527C
528 DO i=1,nel0
529 IF(pla(i)>epsmax1.AND.off(i)==one)off(i)=four_over_5
530 ENDDO
531C------------------------------------------
532C ECROUISSAGE CINE
533C------------------------------------------
534C test bypass
535 IF (fisokin1/=zero) THEN
536 DO i=1,nel0
537 dsxx = sigexx(i) - signxx(i)
538 dsyy = sigeyy(i) - signyy(i)
539 dsxy = sigexy(i) - signxy(i)
540 dexx = (dsxx - nux*dsyy)
541 deyy = (dsyy - nux*dsxx)
542 dexy = two*(one+nux)*dsxy
543 alpha = fisokin1*h(i)/(e1+h(i))/three
544 sigpxx = alpha*(four*dexx+two*deyy)
545 sigpyy = alpha*(four*deyy+two*dexx)
546 sigpxy = alpha*dexy
547 signxx(i) = signxx(i) + uvar(i,2)
548 signyy(i) = signyy(i) + uvar(i,3)
549 signxy(i) = signxy(i) + uvar(i,4)
550 uvar(i,2) = uvar(i,2) + sigpxx
551 uvar(i,3) = uvar(i,3) + sigpyy
552 uvar(i,4) = uvar(i,4) + sigpxy
553 ENDDO
554 ENDIF
555C
556 RETURN
557 END
558C
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps56c
subroutine sigeps56c(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 sigeps56c.F:45
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72