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),
169 . nrate1,ifunc1(100),nfuncv(mvsiz),
170 . nfuncm,iperf
171 INTEGER JST(MVSIZ+1), IC
172 my_real
173 . MX234
174 my_real
175 . DYDX1(MVSIZ),DYDX2(MVSIZ),RATE(MVSIZ,2),SVM(MVSIZ),
176 . Y1(MVSIZ),Y2(MVSIZ),DR(MVSIZ),
177 . YFAC(MVSIZ,2),
178 . aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
179 . pp(mvsiz),qq(mvsiz),fail(mvsiz),h(mvsiz),
180 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz),
181 . hk(mvsiz)
182 my_real
183 . r,umr,a,b,c,s11,s22,s12,p,p2,fac,dezz,
184 . sigz,s1,s2,s3,vm2,epst,nnu2,
185 . f,df,q2,yld_i,sigpxx,sigpyy,sigpxy,alpha,
186 . e1,a11,a21,g1,g31,nnu11,nu11,nu21,nu31,nu41,nu51,nu61,
187 . epsmax1,epsr11,epsr21,fisokin1,
188 . dsxx,dsyy,dsxy,dexx,deyy,dexy,nux
189C
190C -- Added variables for NEC grouping
191 my_real
192 . ng_nu3,
193 . ng_nrates(mvsiz,100)
194
195 integer
196 . imix,imixv(mvsiz),jend(mvsiz),
197 . ifuncm(mvsiz,100)
198C -- End NEC grouping
199 DATA nmax/3/,iperf/0/
200C-----------------------------------------------
201C USER VARIABLES INITIALIZATION
202C-----------------------------------------------
203 nfunc = ipm(10,mat(1))
204 DO j=1,nfunc
205 ifunc1(j)=ipm(10+j,mat(1))
206 ENDDO
207C
208 iadbuf = ipm(7,mat(1))-1
209 e1 = uparam(iadbuf+2)
210 a11 = uparam(iadbuf+3)
211 a21 = uparam(iadbuf+4)
212 g1 = uparam(iadbuf+5)
213 g31 = three*g1
214 nux = uparam(iadbuf+6)
215 nrate1 = nint(uparam(iadbuf+1))
216 epsmax1=uparam(iadbuf+6+2*nrate1+1)
217 IF(epsmax1==zero)THEN
218 IF(tf(npf(ifunc1(1)+1)-1)==zero)THEN
219 epsmax1=tf(npf(ifunc1(1)+1)-2)
220 ELSE
221 epsmax1= ep30
222 ENDIF
223 ENDIF
224C
225 nnu11 = nux / (one - nux)
226 nnu2 = nnu11*nnu11
227 nu11 = one/(one-nux)
228 nu21 = one/(one+nux)
229 nu31 = one-nnu11
230 nu41 = one + nnu2 + nnu11
231 nu51 = one + nnu2 - two*nnu11
232 nu61 = half - nnu2 + half*nnu11
233C
234 epsr11 =uparam(iadbuf+6+2*nrate1+2)
235 epsr21 =uparam(iadbuf+6+2*nrate1+3)
236 fisokin1=uparam(iadbuf+6+2*nrate1+8)
237C
238 IF (isigi==0) THEN
239 IF(time==zero)THEN
240 DO i=1,nel0
241 uvar(i,1)=zero
242 uvar(i,2)=zero
243 uvar(i,3)=zero
244 uvar(i,4)=zero
245 DO j=1,nrate1
246 uvar(i,j+4)=zero
247 ENDDO
248 ENDDO
249 ENDIF
250 ENDIF
251C------------------------------------------
252 DO i=1,nel0
253 dpla_i(i) =zero
254 ENDDO
255C-----------------------------------------------
256C
257 DO i=1,nel0
258C
259 signxx(i)=sigoxx(i) - uvar(i,2) +a11*depsxx(i)+a21*depsyy(i)
260 signyy(i)=sigoyy(i) - uvar(i,3) +a21*depsxx(i)+a11*depsyy(i)
261 signxy(i)=sigoxy(i) - uvar(i,4) +g1 *depsxy(i)
262 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
263 signzx(i)=sigozx(i)+gs(i) *depszx(i)
264 sigexx(i) = signxx(i)
265 sigeyy(i) = signyy(i)
266 sigexy(i) = signxy(i)
267C
268 soundsp(i) = sqrt(a11/rho0(i))
269 viscmax(i) = zero
270 etse(i) = one
271C-------------------
272C STRAIN RATE
273C-------------------
274 IF (israte==0) THEN
275 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
276 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
277 . + epspxy(i)*epspxy(i) ) )
278 ELSE
279 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
280 ENDIF
281C-------------------
282C STRAIN
283C-------------------
284 epst = half*( epsxx(i)+epsyy(i)
285 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
286 . + epsxy(i)*epsxy(i) ) )
287 fail(i) = max(zero,min(one,(epsr21-epst)/(epsr21-epsr11)))
288C
289 ENDDO
290C-------------------
291C CRITERE
292C-------------------
293 DO i=1,nel0
294 jj(i) = 1
295 ENDDO
296 iadbuf = ipm(7,mat(1))-1
297C inversion boucles
298 DO j=2,nrate1-1
299 DO i=1,nel0
300 IF(epsp(i)>=uparam(iadbuf+6+j)) jj(i) = j
301 ENDDO
302 ENDDO
303 DO i=1,nel0
304 rate(i,1)=uparam(iadbuf+6+jj(i))
305 rate(i,2)=uparam(iadbuf+6+jj(i)+1)
306 yfac(i,1)=uparam(iadbuf+6+nrate1+jj(i))
307 yfac(i,2)=uparam(iadbuf+6+nrate1+jj(i)+1)
308 ENDDO
309C
310 DO i=1,nel0
311 j1 = jj(i)
312 j2 = j1+1
313 ipos1(i) = nint(uvar(i,j1+4))
314 iad1(i) = npf(ifunc1(j1)) / 2 + 1
315 ilen1(i) = npf(ifunc1(j1)+1) / 2 - iad1(i) - ipos1(i)
316 ipos2(i) = nint(uvar(i,j2+4))
317 iad2(i) = npf(ifunc1(j2)) / 2 + 1
318 ilen2(i) = npf(ifunc1(j2)+1) / 2 - iad2(i) - ipos2(i)
319 ENDDO
320C
321 CALL vinter(tf,iad1,ipos1,ilen1,nel0,pla,dydx1,y1)
322 CALL vinter(tf,iad2,ipos2,ilen2,nel0,pla,dydx2,y2)
323C
324 IF (fisokin1==zero) THEN
325 DO i=1,nel0
326 j1 = jj(i)
327 j2 = j1+1
328 y1(i)=y1(i)*yfac(i,1)
329 y2(i)=y2(i)*yfac(i,2)
330 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
331 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
332 yld(i) = max(yld(i),em20)
333 dydx1(i)=dydx1(i)*yfac(i,1)
334 dydx2(i)=dydx2(i)*yfac(i,2)
335 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
336 uvar(i,j1+4) = ipos1(i)
337 uvar(i,j2+4) = ipos2(i)
338 ENDDO
339 ELSEIF (fisokin1==1.) THEN
340 DO i=1,nel0
341 j1 = jj(i)
342 j2 = j1+1
343 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
344 dydx1(i)=dydx1(i)*yfac(i,1)
345 dydx2(i)=dydx2(i)*yfac(i,2)
346 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
347 uvar(i,j1+4) = ipos1(i)
348 uvar(i,j2+4) = ipos2(i)
349C ECROUISSAGE CINEMATIQUE
350 y1(i)=tf(npf(ifunc1(j1))+1)
351 y2(i)=tf(npf(ifunc1(j2))+1)
352 y1(i)=y1(i)*yfac(i,1)
353 y2(i)=y2(i)*yfac(i,2)
354 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
355 ENDDO
356 ELSE
357 DO i=1,nel0
358 j1 = jj(i)
359 j2 = j1+1
360 y1(i)=y1(i)*yfac(i,1)
361 y2(i)=y2(i)*yfac(i,2)
362 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
363 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
364 yld(i) = max(yld(i),em20)
365 dydx1(i)=dydx1(i)*yfac(i,1)
366 dydx2(i)=dydx2(i)*yfac(i,2)
367 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
368 uvar(i,j1+4) = ipos1(i)
369 uvar(i,j2+4) = ipos2(i)
370C ECROUISSAGE CINEMATIQUE
371 y1(i)=tf(npf(ifunc1(j1))+1)
372 y2(i)=tf(npf(ifunc1(j2))+1)
373 y1(i)=y1(i)*yfac(i,1)
374 y2(i)=y2(i)*yfac(i,2)
375 yld(i) = (1.-fisokin1) * yld(i) +
376 . fisokin1 * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
377 ENDDO
378 ENDIF
379
380C-------------------
381C PROJECTION
382C-------------------
383 IF(iflag(1)==0)THEN
384C radial projection
385 DO i=1,nel0
386 svm(i)=sqrt(signxx(i)*signxx(i)
387 . +signyy(i)*signyy(i)
388 . -signxx(i)*signyy(i)
389 . +three*signxy(i)*signxy(i))
390 r = min(one,yld(i)/max(em20,svm(i)))
391 signxx(i)=signxx(i)*r
392 signyy(i)=signyy(i)*r
393 signxy(i)=signxy(i)*r
394 umr = one - r
395 dpla_i(i) = off(i)*svm(i)*umr/(g31+h(i))
396 pla(i) = pla(i) + dpla_i(i)
397 s1=half*(signxx(i)+signyy(i))
398 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
399 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
400 thk(i) = thk(i) + dezz*thkly(i)*off(i)
401 IF(r<one) etse(i)= h(i)/(h(i)+e1)
402 ENDDO
403C
404 ELSEIF(iflag(1)==1)THEN
405C-------------------------
406C von mises criterion
407C-------------------------
408 DO i=1,nel0
409 h(i) = max(zero,h(i))
410 s1=signxx(i)+signyy(i)
411 s2=signxx(i)-signyy(i)
412 s3=signxy(i)
413 aa(i)=fourth*s1*s1
414 bb(i)=three_over_4*s2*s2+3.*s3*s3
415 svm(i)=sqrt(aa(i)+bb(i))
416 dezz = -(depsxx(i)+depsyy(i))*nnu11
417 thk(i) = thk(i) + dezz*thkly(i)*off(i)
418 ENDDO
419C-------------------------
420C GATHER PLASTIC FLOW
421C-------------------------
422 nindx=0
423 DO i=1,nel0
424 IF(svm(i)>yld(i).AND.off(i)==one) THEN
425 nindx=nindx+1
426 index(nindx)=i
427 ENDIF
428 ENDDO
429C
430 IF(nindx/=0) THEN
431C---------------------------
432C DEP EN CONTRAINTE PLANE
433C---------------------------
434 DO j=1,nindx
435 i=index(j)
436 dpla_j(i)=(svm(i)-yld(i))/(g31+h(i))
437 etse(i)= h(i)/(h(i)+e1)
438 hk(i) = h(i)*(one-fisokin1)
439 ENDDO
440C
441 DO n=1,nmax
442#include "vectorize.inc"
443 DO j=1,nindx
444 i=index(j)
445 dpla_i(i)=dpla_j(i)
446 yld_i =yld(i)+hk(i)*dpla_i(i)
447 dr(i) =half*e1*dpla_i(i)/yld_i
448 pp(i) =one/(one+dr(i)*nu11)
449 qq(i) =one/(one+three*dr(i)*nu21)
450 p2 =pp(i)*pp(i)
451 q2 =qq(i)*qq(i)
452 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
453 df =-(aa(i)*nu11*p2*pp(i)+three*bb(i)*nu21*q2*qq(i))
454 . *(e1-two*dr(i)*hk(i))/yld_i
455 . -two*hk(i)*yld_i
456 df = sign(max(abs(df),em20),df)
457 IF(dpla_i(i)>zero) THEN
458 dpla_j(i)=max(zero,dpla_i(i)-f/df)
459 ELSE
460 dpla_j(i)=zero
461 ENDIF
462 ENDDO
463 ENDDO
464C------------------------------------------
465C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
466C------------------------------------------
467#include "vectorize.inc"
468 DO j=1,nindx
469 i=index(j)
470 pla(i) = pla(i) + dpla_i(i)
471 s1=(signxx(i)+signyy(i))*pp(i)
472 s2=(signxx(i)-signyy(i))*qq(i)
473 signxx(i)=half*(s1+s2)
474 signyy(i)=half*(s1-s2)
475 signxy(i)=signxy(i)*qq(i)
476 dezz = - nu31*dr(i)*s1/e1
477 thk(i) = thk(i) + dezz*thkly(i)*off(i)
478 ENDDO
479 ENDIF
480C-------------------------------------------
481 ELSEIF(iflag(1)==2)THEN
482C radial projection on the deviator on a reduced criterion
483C Elastic projection in z => sig33 = 0
484C the criterion reduction coefficient is such that
485C one is on the criterion after the 2 projections
486 DO i=1,nel0
487 p = -(signxx(i)+signyy(i))*third
488 s11 = signxx(i)+p
489 s22 = signyy(i)+p
490C s33 = p = -(S11 + S22)
491 s12 = signxy(i)
492C
493 p2 = p*p
494 vm2= three*(s12*s12 - s11*s22)
495 a = p2*nu41 + vm2
496 vm2= three*p2 + vm2
497 b = p2*nu61
498 c = p2*nu51 - yld(i)*yld(i)
499 r = min(one,(-b + sqrt(max(zero,b*b-a*c)))/max(a ,em20))
500 signxx(i) = s11*r - p
501 signyy(i) = s22*r - p
502 signxy(i) = s12*r
503C signzz = p*r - p
504C proj. signzz = 0.
505 umr = one - r
506 sigz = nnu11*p*umr
507 signxx(i) = signxx(i) + sigz
508 signyy(i) = signyy(i) + sigz
509 svm(i)=sqrt(vm2)
510 dpla_i(i) = off(i)*svm(i)*umr/(g31+h(i))
511 pla(i) = pla(i) + dpla_i(i)
512 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
513 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
514 thk(i) = thk(i) + dezz*thkly(i)*off(i)
515 IF(r<one) etse(i)= h(i)/(h(i)+e1)
516 ENDDO
517 ENDIF
518C
519 DO i=1,nel0
520 IF(pla(i)>epsmax1.AND.off(i)==one)off(i)=four_over_5
521 ENDDO
522C------------------------------------------
523C ECROUISSAGE CINE
524C------------------------------------------
525C test bypass
526 IF (fisokin1/=zero) THEN
527 DO i=1,nel0
528 dsxx = sigexx(i) - signxx(i)
529 dsyy = sigeyy(i) - signyy(i)
530 dsxy = sigexy(i) - signxy(i)
531 dexx = (dsxx - nux*dsyy)
532 deyy = (dsyy - nux*dsxx)
533 dexy = two*(one+nux)*dsxy
534 alpha = fisokin1*h(i)/(e1+h(i))/three
535 sigpxx = alpha*(four*dexx+two*deyy)
536 sigpyy = alpha*(four*deyy+two*dexx)
537 sigpxy = alpha*dexy
538 signxx(i) = signxx(i) + uvar(i,2)
539 signyy(i) = signyy(i) + uvar(i,3)
540 signxy(i) = signxy(i) + uvar(i,4)
541 uvar(i,2) = uvar(i,2) + sigpxx
542 uvar(i,3) = uvar(i,3) + sigpyy
543 uvar(i,4) = uvar(i,4) + sigpxy
544 ENDDO
545 ENDIF
546C
547 RETURN
548 END
549C
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:73