OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps65.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!|| sigeps65 ../engine/source/materials/mat/mat065/sigeps65.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| eqsolv_2 ../engine/source/materials/mat/mat065/sigeps65c.F
29!|| finter ../engine/source/tools/curve/finter.F
30!|| finter2 ../engine/source/tools/curve/vinter.F
31!|| interstar ../engine/source/materials/mat/mat065/sigeps65c.F
32!||====================================================================
33 SUBROUTINE sigeps65(
34 1 NEL ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,NPF ,
35 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,RHO ,
36 3 VOLUME ,EINT ,
37 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
38 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
39 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
40 7 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
41 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 9 SIGVXX ,SIGVYY ,SIGVZZ ,SIGVXY ,SIGVYZ ,SIGVZX ,
43 A SOUNDSP,VISCMAX,UVAR ,OFF ,NGL ,IPM ,
44 B MAT ,EPSP ,IPLA ,YLD ,PLA ,ETSE ,
45 C DPLA ,AMU )
46
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51#include "comlock.inc"
52C-----------------------------------------------
53C G l o b a l P a r a m e t e r s
54C-----------------------------------------------
55C---------+---------+---+---+--------------------------------------------
56C VAR | SIZE |TYP| RW| DEFINITION
57C---------+---------+---+---+--------------------------------------------
58C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
59C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
60C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
61C---------+---------+---+---+--------------------------------------------
62C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
63C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
64C NPF | * | I | R | FUNCTION ARRAY
65C TF | * | F | R | FUNCTION ARRAY
66C---------+---------+---+---+--------------------------------------------
67C TIME | 1 | F | R | CURRENT TIME
68C TIMESTEP| 1 | F | R | CURRENT TIME STEP
69C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
70C RHO0 | NEL | F | R | INITIAL DENSITY
71C RHO | NEL | F | R | DENSITY
72C VOLUME | NEL | F | R | VOLUME
73C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
74C EPSPXX | NEL | F | R | STRAIN RATE XX
75C EPSPYY | NEL | F | R | STRAIN RATE YY
76C ... | | | |
77C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
78C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
79C ... | | | |
80C EPSXX | NEL | F | R | STRAIN XX
81C EPSYY | NEL | F | R | STRAIN YY
82C ... | | | |
83C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
84C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
85C ... | | | |
86C---------+---------+---+---+--------------------------------------------
87C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
88C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
89C ... | | | |
90C SIGVXX | NEL | F | W | VISCOUS STRESS XX
91C SIGVYY | NEL | F | W | VISCOUS STRESS YY
92C ... | | | |
93C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
94C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
95C---------+---------+---+---+--------------------------------------------
96C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
97C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
98C---------+---------+---+---+--------------------------------------------
99#include "com01_c.inc"
100#include "param_c.inc"
101C=======================================================================
102C-----------------------------------------------
103C I N P U T A r g u m e n t s
104C-----------------------------------------------
105 INTEGER NEL, NUPARAM, NUVAR,IPT,IPLA, JTHE,
106 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*)
107 my_real
108 . TIME,TIMESTEP,UPARAM(*),
109 . RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
110 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
111 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
112 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
113 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
114 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
115 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
116 . SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
117 . sigoxy(nel),sigoyz(nel),sigozx(nel),
118 . epsp(nel),etse(nel), tempel(nel),amu(nel)
119C-----------------------------------------------
120C O U T P U T A r g u m e n t s
121C-----------------------------------------------
122 my_real
123 . signxx(nel),signyy(nel),signzz(nel),
124 . signxy(nel),signyz(nel),signzx(nel),
125 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
126 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
127 . soundsp(nel),viscmax(nel),dpla(nel)
128C-----------------------------------------------
129C I N P U T O U T P U T A r g u m e n t s
130C-----------------------------------------------
131 my_real
132 . uvar(nel,nuvar), off(nel), yld(nel), pla(nel)
133C-----------------------------------------------
134C VARIABLES FOR FUNCTION INTERPOLATION
135C-----------------------------------------------
136 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
137 my_real
138 . FINTER2, TF(*),FINTER
139 EXTERNAL FINTER2,FINTER
140C-----------------------------------------------
141C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
142C Y : y = f(x)
143C X : x
144C DYDX : f'(x) = dy/dx
145C IFUNC(J): FUNCTION INDEX
146C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
147C NPF,TF : FUNCTION PARAMETER
148C-----------------------------------------------
149C L o c a l V a r i a b l e s
150C-----------------------------------------------
151 INTEGER I,I1,I2,J,IADBUF,J1,J2,NFUNC,FUNC,FUND,
152 . NP1,NP2,K1,K,NRATE,N0,N1,N2,N3,N4
153 INTEGER JJ(NEL),MX,
154 . IPOSC1(NEL),IADC1(NEL),ILENC1(NEL),
155 . IPOSC2(NEL),IADC2(NEL),ILENC2(NEL),
156 . IPOSD1(NEL),IADD1(NEL),ILEND1(NEL),FLAG(NEL),
157 . IPOSD2(NEL),IADD2(NEL),ILEND2(NEL),IFUNC(100)
158
159C
160 my_real
161 . S1,S2,T1,T2,X1,X2,Y1,Y2,SX,TY,DYDX,DTDS,RATE1,
162 . RATE2,SIGY1,SIGY2,SVTEST,TEST,HH,DESEL
163 my_real
164 . e(nel),nu,g3,g,eps(nel),c1,
165 . epsmax,trsigo(nel),nuu(nel),epss(nel),
166 . eps0(nel),depseq(nel),dsigxx(nel),dsigyy(nel),
167 . dsigzz(nel),sv2(nel),sigo(nel),svm(nel),
168 . strxx(nel),stryy(nel),strzz(nel),sigm(nel),
169 . strxy(nel),stryz(nel),strzx(nel),trsign(nel),
170 . rbulk,g2,
171 . dydxc1(nel),dydxc2(nel),yc1(nel),yc2(nel),
172 . dydxd1(nel),dydxd2(nel),yd1(nel),yd2(nel),
173 . sigyld(nel),hc(nel),hd(nel),sigfc(nel), sigfd(nel),
174 . fac(nel),epse(nel),yfac(nel,2),depss(nel),h(nel),
175 . depse(nel),deq(nel),yield(nel),unload(nel),
176 . epsc1(nel),epsc2(nel),epsd1(nel),epsd2(nel),work(nel),
177 . dydxc(nel),yc(nel),yfac1(nel),yfac2(nel),coefb(nel),
178 . dexx(nel),deyy(nel),dezz(nel),p0(nel),coefa(nel),
179 . soxx(nel),soyy(nel),sozz(nel) ,svmo2(nel),dsig(nel),
180 . a,b,cc,x11,x22,dav,alpha,sfc,sfd,p
181C
182C-----------------------------------------------
183C USER VARIABLES INITIALIZATION
184C-----------------------------------------------
185 mx = mat(1)
186 nfunc = ipm(10,mx)
187 DO j=1,nfunc
188 ifunc(j)=ipm(10+j,mx)
189 ENDDO
190 iadbuf = ipm(7,mx) - 1
191 nrate = nint(uparam(iadbuf+1))
192 g = uparam(iadbuf+3)
193 nu = uparam(iadbuf+4)
194 g2 = uparam(iadbuf+7)
195 g3 = uparam(iadbuf+8)
196 epsmax= uparam(iadbuf+11)
197 c1=uparam(iadbuf+2)/three/(one - two*nu)
198C
199 rbulk = uparam(iadbuf+2)*third/(one-two*nu)
200 DO i=1,nel
201 e(i) = uparam(iadbuf+2)
202 ENDDO
203 nrate = nint(uparam(iadbuf+1))
204 n0 = 5
205 n1 = nrate + n0
206 n2 = nrate + n1
207 n3 = nrate + n2
208 n4 = nrate + n3 + 1 !! = 5 + 4*NRATE +1
209c------------------------------------------
210 IF (time == zero)THEN
211 IF (isigi == 0) THEN
212 DO i=1,nel
213 DO j=1,10
214 uvar(i,j)=zero
215 ENDDO
216C --- OLD EFFECTIVE STRESS
217 trsigo(i)=-(sigoxx(i)+sigoyy(i)+sigozz(i))*third
218 uvar(i,3) = sqrt(three_half*(
219 . (sigoxx(i)+trsigo(i))*(signxx(i)+trsigo(i))
220 . +(sigoyy(i)+trsigo(i))*(signyy(i)+trsigo(i))
221 . +(sigozz(i)+trsigo(i))*(signzz(i)+trsigo(i))
222 . +two*(sigoxy(i)*sigoxy(i)+sigoyz(i)*sigoyz(i)
223 . +sigozx(i)*sigozx(i))))
224 ENDDO
225 ENDIF
226c---- find yield values
227 DO i=1,nel
228 DO i1=1,nrate
229 func = ifunc(i1)
230 fund = ifunc(i1+nrate)
231 np1 = (npf(func+1)-npf(func))/2
232 np2 = (npf(fund+1)-npf(fund))/2
233c---
234 DO j=3,np1
235 j1=2*(j-2)
236 s1=tf(npf(func)+j1)
237 s2=tf(npf(func)+j1+2)
238 t1=tf(npf(func)+j1+1)
239 t2=tf(npf(func)+j1+3)
240 DO k=3,np2
241 k1=2*(k-2)
242 x1=tf(npf(fund)+k1)
243 x2=tf(npf(fund)+k1+2)
244 y1=tf(npf(fund)+k1+1)
245 y2=tf(npf(fund)+k1+3)
246 IF (y2>=t1 .AND. y1<=t2 .AND. x2>=s1 .AND. x1<=s2) THEN
247 dydx = (y2-y1) / (x2-x1)
248 dtds = (t2-t1) / (s2-s1)
249 IF (dydx > dtds) THEN
250 sx = (t1-y1-dtds*s1+dydx*x1) / (dydx-dtds)
251 ty = t1 + dtds*(sx - s1)
252 IF (ty>=y1 .AND. ty<=y2 .AND. sx>=x1 .AND. sx<=x2)THEN
253 iadbuf = ipm(7,mx) + 13
254 uparam(iadbuf+i1+nrate*2) = ty
255 GOTO 150
256 ENDIF
257 ENDIF
258 ENDIF
259 ENDDO
260 ENDDO
261 150 CONTINUE
262c---
263 ENDDO
264 ENDDO
265 ENDIF !TIME == ZERO
266C-------------------
267C STRAIN
268C-------------------
269 DO i=1,nel
270 epss(i) = uvar(i,1)
271 epse(i) = uvar(i,2)
272 ENDDO
273
274 DO i=1,nel
275 eps(i) = (one/(one+nu))*sqrt((epsxx(i)*epsxx(i)+epsyy(i)*epsyy(i)+
276 . epszz(i)*epszz(i)+two*(epsxy(i)*epsxy(i)+
277 . epsyz(i)*epsyz(i)+epszx(i)*epszx(i))))
278 depseq(i) = (one/(one+nu))*sqrt((depsxx(i)*depsxx(i) + depsyy(i)*
279 . depsyy(i)+depszz(i)*depszz(i)+two*(depsxy(i)*depsxy(i)+
280 . depsyz(i)*depsyz(i) +depszx(i)*depszx(i))))
281 ENDDO
282C-------------------
283C ELASTIC STRESS
284C-------------------
285 DO i=1,nel
286C---
287c ..strain increment deviator ..
288 dav = (depsxx(i)+depsyy(i)+depszz(i))*third
289c ..pressure..
290 p0(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third
291c ..strain increment deviator tensor ..
292 dexx(i) = depsxx(i)-dav
293 deyy(i) = depsyy(i)-dav
294 dezz(i) = depszz(i)-dav
295c ..stress increment deviator tensor ..
296 soxx(i) =sigoxx(i)+p0(i)
297 soyy(i) =sigoyy(i)+p0(i)
298 sozz(i) =sigozz(i)+p0(i)
299
300 svmo2(i) = two_third*uvar(i,3)*uvar(i,3)
301c ..deviatoric elastic predictor ..
302 signxx(i)=soxx(i)+g2*dexx(i)
303 signyy(i)=soyy(i)+g2*deyy(i)
304 signzz(i)=sozz(i)+g2*dezz(i)
305 signxy(i)=sigoxy(i)+g *depsxy(i)
306 signyz(i)=sigoyz(i)+g *depsyz(i)
307 signzx(i)=sigozx(i)+g *depszx(i)
308
309C
310 sv2(i) = three_half*(signxx(i)*signxx(i)+signyy(i)*signyy(i)
311 . +signzz(i) *signzz(i)
312 . +two*(signxy(i)*signxy(i)
313 . +signyz(i)*signyz(i)
314 . +signzx(i)*signzx(i)))
315 svm(i) = sqrt(sv2(i))
316
317
318 p = c1 * (rho(i)/rho0(i)-one)
319 work(i) = (signxx(i)-p)*depsxx(i) +
320 . (signyy(i)-p)*depsyy(i) +
321 . (signzz(i)-p)*depszz(i) +
322 . two*signxy(i) *depsxy(i) +
323 . two*signyz(i) *depsyz(i) +
324 . two*signzx(i) *depszx(i)
325
326 dsig(i) = sqrt(three_half* ( (g2*dexx(i))**2+
327 . (g2*deyy(i))**2+
328 . (g2*dezz(i))**2+
329 . two*((g *depsxy(i))**2+
330 . (g *depsyz(i))**2+
331 . (g *depszx(i))**2) ))
332
333C-------------------
334 coefa(i) = two* g*g *(two*(dexx(i)*dexx(i)+
335 . deyy(i)*deyy(i) + dezz(i)*dezz(i)) +
336 . depsxy(i)*depsxy(i)+depsyz(i)*depsyz(i)+
337 . depszx(i)*depszx(i))
338
339 coefb(i) = four *g*(soxx(i)*dexx(i)
340 . +soyy(i)*deyy(i)+sozz(i)*dezz(i)+sigoxy(i)*depsxy(i)
341 . +sigoyz(i)*depsyz(i)+sigozx(i)*depszx(i))
342
343C-------------------
344C SOUND SPEED
345C-------------------
346 soundsp(i) = sqrt((rbulk+four_over_3*g)/rho0(i))
347 viscmax(i) = zero
348 ENDDO
349C-------------------
350C Interpolation
351C-------------------
352 iadbuf = ipm(7,mat(1))-1+14
353 mx = mat(1)
354 DO i=1,nel
355 jj(i) = 1
356 DO j=2,nrate-1
357 IF (epsp(i) > uparam(iadbuf+j)) jj(i) = j
358 ENDDO
359C
360 i1 = iadbuf+jj(i)
361 j1 = jj(i)
362 rate1 = uparam(i1)
363 yfac1(i) = uparam(i1+nrate)
364 sigy1 = uparam(i1+nrate*2)*yfac1(i)
365C
366 IF (nrate > 1) THEN
367 i2 = i1 + 1
368 j2 = j1 + 1
369c
370 rate2 = uparam(i2)
371 yfac2(i) = uparam(i2+nrate)
372 fac(i) = (epsp(i) - rate1) / (rate2 - rate1)
373 sigy2 = uparam(i2+nrate*2)*yfac2(i)
374 sigyld(i) = sigy1 + fac(i)*(sigy2-sigy1)
375c loading curves
376 iposc1(i) = nint(uvar(i,n0+j1))
377 iposc2(i) = nint(uvar(i,n0+j2))
378 epsc1(i) = uvar(i,n1+j1)
379 epsc2(i) = uvar(i,n1+j2)
380 iadc1(i) = npf(ipm(10+j1,mx))/2 + 1
381 iadc2(i) = npf(ipm(10+j2,mx))/2 + 1
382 ilenc1(i) = npf(ipm(10+j1,mx)+1)/2-iadc1(i)-iposc1(i)
383 ilenc2(i) = npf(ipm(10+j2,mx)+1)/2-iadc2(i)-iposc2(i)
384c unloading curves
385 iposd1(i) = nint(uvar(i,n2+j1))
386 iposd2(i) = nint(uvar(i,n2+j2))
387 epsd1(i) = uvar(i,n3+j1)
388 epsd2(i) = uvar(i,n3+j2)
389 iadd1(i) = npf(ipm(10+nrate+j1,mx))/2 + 1
390 iadd2(i) = npf(ipm(10+nrate+j2,mx))/2 + 1
391 ilend1(i) = npf(ipm(10+nrate+j1,mx)+1)/2
392 . - iadd1(i) - iposd1(i)
393 ilend2(i) = npf(ipm(10+nrate+j2,mx)+1)/2
394 . - iadd2(i) - iposd2(i)
395 ELSE
396 sigyld(i) = sigy1
397 iposc1(i) = nint(uvar(i,n0+j1))
398 iposd1(i) = nint(uvar(i,n2+j1))
399 epsc1(i) = uvar(i,n1+j1)
400 epsd1(i) = uvar(i,n3+j1)
401
402 iadc1(i) = npf(ipm(10+j1,mx))/2 + 1
403 iadd1(i) = npf(ipm(10+nrate+j1,mx))/2 + 1
404 ilenc1(i) = npf(ipm(10+j1,mx)+1)/2-iadc1(i)-iposc1(i)
405 ilend1(i) = npf(ipm(10+nrate+j1,mx)+1)/2
406 . - iadd1(i) - iposd1(i)
407 ENDIF
408 ENDDO
409 IF (nrate > 1) THEN
410c Sigma load = f(total strain)
411 CALL interstar(tf ,iadc1,iposc1,ilenc1,nel ,
412 . epss ,epsc1,dydxc1,e ,yc1 ,yfac1 )
413 CALL interstar(tf ,iadc2,iposc2,ilenc2,nel ,
414 . epss ,epsc2,dydxc2,e ,yc2 ,yfac2 )
415c Sigma unload = f(elastic strain)
416 CALL interstar(tf ,iadd1,iposd1,ilend1,nel ,
417 . epse ,epsd1,dydxd1,e ,yd1 ,yfac1 )
418 CALL interstar(tf ,iadd2,iposd2,ilend2,nel ,
419 . epse ,epsd2,dydxd2,e ,yd2 ,yfac2 )
420 ELSE
421 CALL interstar(tf ,iadc1,iposc1,ilenc1,nel ,
422 . epss ,epsc1,dydxc1,e ,yc1 ,yfac1 )
423 CALL interstar(tf ,iadd1,iposd1,ilend1,nel ,
424 . epse ,epsd1,dydxd1,e ,yd1 ,yfac1 )
425 ENDIF
426c-------------------
427 IF (nrate > 1) THEN
428 DO i=1,nel
429 j1 = jj(i)
430 j2 = j1 + 1
431 uvar(i,n0+j1) = iposc1(i)
432 uvar(i,n0+j2) = iposc2(i)
433 uvar(i,n1+j1) = epsc1(i)
434 uvar(i,n1+j2) = epsc2(i)
435 uvar(i,n2+j1) = iposd1(i)
436 uvar(i,n2+j2) = iposd2(i)
437 uvar(i,n3+j1) = epsd1(i)
438 uvar(i,n3+j2) = epsd2(i)
439c
440 yc1(i) = yc1(i)
441 yd1(i) = yd1(i)
442 yc2(i) = yc2(i)
443 yd2(i) = yd2(i)
444 dydxc1(i)= dydxc1(i)
445 dydxd1(i)= dydxd1(i)
446 dydxc2(i)= dydxc2(i)
447 dydxd2(i)= dydxd2(i)
448c
449 hc(i) = dydxc1(i) + (dydxc2(i)-dydxc1(i))*fac(i)
450 hd(i) = dydxd1(i) + (dydxd2(i)-dydxd1(i))*fac(i)
451 sigfc(i) = yc1(i) + (yc2(i)-yc1(i))*fac(i)
452 sigfd(i) = yd1(i) + (yd2(i)-yd1(i))*fac(i)
453c
454 ENDDO
455 ELSE
456 DO i=1,nel
457 j1 = jj(i)
458 uvar(i,n0+j1) = iposc1(i)
459 uvar(i,n1+j1) = epsc1(i)
460 uvar(i,n2+j1) = iposd1(i)
461 uvar(i,n3+j1) = epsd1(i)
462c
463 sigfc(i) = yc1(i)
464 sigfd(i) = yd1(i)
465 hc(i) = dydxc1(i)
466 hd(i) = dydxd1(i)
467 ENDDO
468 ENDIF
469c-------------------
470c PROJECTION
471c-------------------
472 DO i=1,nel
473 depss(i) = zero
474 depse(i) = zero
475 deq(i) = eps(i) -uvar(i,n4+1)
476 flag(i) = int(uvar(i,n4+2) )
477C
478 IF (svm(i) < sigfd(i))THEN !decharge
479 flag(i) = -1
480 depss(i) = (svm(i) - sigfd(i)) / (e(i) + hd(i))
481 sfd = sigfd(i)
482 sigfd(i) = sigfd(i) + hd(i)*depss(i)
483 desel = ( sigfd(i) - sfd ) / e(i)
484 test = desel + depss(i)
485 IF((depseq(i)-abs(test))>em15.AND. dsig(i)> sqrt(svmo2(i)))THEN
486 depss(i) = -(svm(i) + sfd) / (e(i) + hd(i))
487 sfd = sigfd(i)
488 sigfd(i) = sigfd(i) + hd(i)*depss(i)
489 ENDIF
490 a = coefa(i)
491 b = coefb(i)
492 cc = svmo2(i) - two_third * sigfd(i)*sigfd(i)
493 CALL eqsolv_2(a,b,cc,x11,x22)
494 x11 = max(x11,zero)
495 x22 = max(x22,zero)
496 alpha = min(x11,x22)
497c
498 signxx(i)=soxx(i)+alpha*g2*dexx(i)
499 signyy(i)=soyy(i)+alpha*g2*deyy(i)
500 signzz(i)=sozz(i)+alpha*g2*dezz(i)
501 signxy(i)=sigoxy(i)+alpha*g *depsxy(i)
502 signyz(i)=sigoyz(i)+alpha*g *depsyz(i)
503 signzx(i)=sigozx(i)+alpha*g *depszx(i)
504C
505 ELSEIF (svm(i) >= sigfc(i)) THEN ! .AND.WORK(I)>ZERO
506 flag(i) = 1
507 depss(i) = (svm(i) - sigfc(i)) / (e(i) + hc(i))
508 sfc = sigfc(i) + hc(i)*depss(i)
509 sfd = sigfd(i) + hd(i)*depss(i)
510 dpla(i)= zero
511 a = coefa(i)
512 b = coefb(i)
513 cc = svmo2(i) - two_third * sfc*sfc
514 CALL eqsolv_2(a,b,cc,x11,x22)
515 alpha = max(x11,x22)
516 signxx(i)=soxx(i)+alpha*g2*dexx(i)
517 signyy(i)=soyy(i)+alpha*g2*deyy(i)
518 signzz(i)=sozz(i)+alpha*g2*dezz(i)
519 signxy(i)=sigoxy(i)+alpha*g *depsxy(i)
520 signyz(i)=sigoyz(i)+alpha*g *depsyz(i)
521 signzx(i)=sigozx(i)+alpha*g *depszx(i)
522 !-----------------------------------------
523 IF (sfd > sfc) THEN !plastic
524 depss(i) = (svm(i) - sigfc(i)) / (g3 + hc(i))
525 sigfc(i) = sigfc(i) + hc(i)*depss(i)
526 sigfd(i) = sigfd(i) + hd(i)*depss(i)
527 a = coefa(i)
528 b = coefb(i)
529 cc = svmo2(i) - two_third * sigfc(i)*sigfc(i)
530 CALL eqsolv_2(a,b,cc,x11,x22)
531 alpha = max(x11,x22)
532 signxx(i)=soxx(i)+alpha*g2*dexx(i)
533 signyy(i)=soyy(i)+alpha*g2*deyy(i)
534 signzz(i)=sozz(i)+alpha*g2*dezz(i)
535 signxy(i)=sigoxy(i)+alpha*g *depsxy(i)
536 signyz(i)=sigoyz(i)+alpha*g *depsyz(i)
537 signzx(i)=sigozx(i)+alpha*g *depszx(i)
538 dpla(i) = depss(i)
539 pla(i) = pla(i) + off(i)*dpla(i)
540 ELSE
541 sigfc(i) = sfc
542 sigfd(i) = sfd
543 dpla(i)= zero
544 ENDIF
545 ELSEIF ( svm(i) < sigfc(i))THEN !cas elastique,pente E => depss = 0
546 dpla(i) = zero
547 IF (flag(i) == -1 .AND. work(i) < zero.AND.epss(i)>pla(i))THEN !passage cote negatif,faux svm
548 flag(i) = -1
549 yield(i)=svm(i)+sigfd(i)
550 depss(i) = -(yield(i)) / (e(i) + hd(i))
551 sigfd(i) = sigfd(i) + hd(i)*depss(i)
552 a = coefa(i)
553 b = coefb(i)
554 cc = svmo2(i) - two_third * sigfd(i)*sigfd(i)
555 CALL eqsolv_2(a,b,cc,x11,x22)
556 x11 = max(x11,zero)
557 x22 = max(x22,zero)
558 alpha = min(x11,x22)
559 signxx(i)=soxx(i)+alpha*g2*dexx(i)
560 signyy(i)=soyy(i)+alpha*g2*deyy(i)
561 signzz(i)=sozz(i)+alpha*g2*dezz(i)
562 signxy(i)=sigoxy(i)+alpha*g *depsxy(i)
563 signyz(i)=sigoyz(i)+alpha*g *depsyz(i)
564 signzx(i)=sigozx(i)+alpha*g *depszx(i)
565 ENDIF
566C
567 ENDIF
568 ENDDO
569c----------------------------------------------
570 DO i=1,nel
571 svm(i) = sqrt(three_half*(signxx(i)*signxx(i)+signyy(i)*signyy(i)
572 . +signzz(i) *signzz(i)
573 . +two*(signxy(i)*signxy(i)
574 . +signyz(i)*signyz(i)
575 . +signzx(i)*signzx(i))))
576c P = C1 * (RHO(I)/RHO0(I)-ONE)
577 p = c1 * amu(i)
578 signxx(i)=signxx(i)-p
579 signyy(i)=signyy(i)-p
580 signzz(i)=signzz(i)-p
581
582 IF( dpla(i) > zero)THEN
583 h(i)=(svm(i)-uvar(i,3))/dpla(i)
584 uvar(i,n4)=h(i)/g2
585 ENDIF
586
587 epss(i) = epss(i) + depss(i)
588 epse(i) = epss(i) - pla(i)
589
590 uvar(i,1) = epss(i)
591 uvar(i,2) = epse(i)
592 uvar(i,3) = svm(i)
593 uvar(i,4) = pla(i)
594 uvar(i,5) = epsp(i)
595 uvar(i,n4+1) = eps(i)
596 uvar(i,n4+2) = flag(i)
597
598
599 etse(i) = uvar(i,n4)
600 IF(pla(i) > epsmax .AND. off(i) == one) off(i)=four_over_5
601 ENDDO
602 RETURN
603 END
604C
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
sigeps65
subroutine sigeps65(nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, ngl, ipm, mat, epsp, ipla, yld, pla, etse, dpla, amu)
Definition sigeps65.F:46
eqsolv_2
subroutine eqsolv_2(a, b, c, x1, x2)
Definition sigeps65c.F:539
interstar
subroutine interstar(tf, iad, ipos, ilen, nel, x, epss_i, dydx, e, y, yfac)
Definition sigeps65c.F:579