OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps77.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!|| sigeps77 ../engine/source/materials/mat/mat077/sigeps77.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||====================================================================
30 SUBROUTINE sigeps77(
31 1 NEL ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,NPF ,
32 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,RHO ,
33 3 VOLUME ,EINT ,
34 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
35 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
36 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
37 7 SIG0XX ,SIG0YY ,SIG0ZZ ,SIG0XY ,SIG0YZ ,SIG0ZX ,
38 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
39 9 SIGVXX ,SIGVYY ,SIGVZZ ,SIGVXY ,SIGVYZ ,SIGVZX ,
40 A SOUNDSP,FSSP ,VISCMAX,UVAR ,OFF ,NGL ,
41 B PM ,IPM , MAT ,EPSP ,P_AIR )
42C-----------------------------------------------
43C I m p l i c i t T y p e s
44C-----------------------------------------------
45#include "implicit_f.inc"
46#include "comlock.inc"
47C-----------------------------------------------
48C G l o b a l P a r a m e t e r s
49C-----------------------------------------------
50#include "mvsiz_p.inc"
51C---------+---------+---+---+--------------------------------------------
52C VAR | SIZE |TYP| RW| DEFINITION
53C---------+---------+---+---+--------------------------------------------
54C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
55C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
56C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
57C---------+---------+---+---+--------------------------------------------
58C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
59C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
60C NPF | * | I | R | FUNCTION ARRAY
61C TF | * | F | R | FUNCTION ARRAY
62C---------+---------+---+---+--------------------------------------------
63C TIME | 1 | F | R | CURRENT TIME
64C TIMESTEP| 1 | F | R | CURRENT TIME STEP
65C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
66C RHO0 | NEL | F | R | INITIAL DENSITY
67C RHO | NEL | F | R | DENSITY
68C VOLUME | NEL | F | R | VOLUME
69C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
70C EPSPXX | NEL | F | R | STRAIN RATE XX
71C EPSPYY | NEL | F | R | STRAIN RATE YY
72C ... | | | |
73C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
74C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
75C ... | | | |
76C EPSXX | NEL | F | R | STRAIN XX
77C EPSYY | NEL | F | R | STRAIN YY
78C ... | | | |
79C SIG0XX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
80C SIG0YY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
81C ... | | | |
82C---------+---------+---+---+--------------------------------------------
83C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
84C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
85C ... | | | |
86C SIGVXX | NEL | F | W | VISCOUS STRESS XX
87C SIGVYY | NEL | F | W | VISCOUS STRESS YY
88C ... | | | |
89C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
90C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
91C---------+---------+---+---+--------------------------------------------
92C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
93C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
94C---------+---------+---+---+--------------------------------------------
95C---------+---------+---+---+--------------------------------------------
96C-----------------------------------------------
97C air UVAR(1) = rho
98C UVAR(2) = flow_vol
99C UVAR(3) = E_int
100C UVAR(4) = flow_dvol
101C UVAR(5) = SIGXX
102C UVAR(6 ) = SIGYY
103C UVAR(7 ) = SIGZZ
104C UVAR(8 ) = SIGXY
105C UVAR(9 ) = SIGYZ
106C UVAR(10 ) = SIGZX
107C foam variable
108C UVAR(10 + ..)
109#include "param_c.inc"
110C
111 INTEGER NEL, NUPARAM, NUVAR,IPT,
112 . NGL(NEL),MAT(NEL),IPLA,IPM(NPROPMI,*)
113 my_real
114 . TIME,TIMESTEP,UPARAM(*),
115 . RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
116 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
117 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
118 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
119 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
120 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
121 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
122 . sig0xx(nel),sig0yy(nel),sig0zz(nel),
123 . sig0xy(nel),sig0yz(nel),sig0zx(nel),
124 . pm(npropm,*),epsp(nel),fssp(nel)
125C-----------------------------------------------
126C O U T P U T A r g u m e n t s
127C-----------------------------------------------
128 my_real
129 . signxx(nel),signyy(nel),signzz(nel),
130 . signxy(nel),signyz(nel),signzx(nel),
131 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
132 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
133 . soundsp(nel),viscmax(nel)
134C-----------------------------------------------
135C I N P U T O U T P U T A r g u m e n t s
136C-----------------------------------------------
137 my_real
138 . uvar(nel,nuvar), off(nel), pla(nel)
139C-----------------------------------------------
140C VARIABLES FOR FUNCTION INTERPOLATION
141C-----------------------------------------------
142 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
143 my_real
144 . TF(*),FINTER
145 EXTERNAL FINTER
146C EXTERNAL FINTER
147C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
148C Y : y = f(x)
149C X : x
150C DYDX : f'(x) = dy/dx
151C IFUNC(J): FUNCTION INDEX
152C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
153C NPF,TF : FUNCTION PARAMETER
154C-----------------------------------------------
155C L o c a l V a r i a b l e s
156C-----------------------------------------------
157 INTEGER I,J,J1,J2,I1,I2,IADBUF,NC,NFUNC,ILOAD(MVSIZ),
158 . IFUNC(MVSIZ,100),NLOAD,NUNLOAD,
159 . IE_CST(MVSIZ),ILOAD0,IDAMAGE,NPAR_FOAM,IFUNCR,
160 . IFUNCK
161 my_real
162 . R,FAC,YP1,YP2,YN1,YN2,COEFF,RMIN,RMAX,
163 . E,E1,E2,E3,E4,E5,E6,DF,BB,CC,DELTA,ALPHA,X1,X2,SVM2,
164 . SVM,RHO_AIR,E_AIR,V0,DVOL,MU ,GAMA,RHO_AIR0,P0,
165 . E0,G(MVSIZ),NU,AA1(MVSIZ),AA2(MVSIZ),PFOAM,
166 . FAIL(MVSIZ),
167 . EPST(MVSIZ),AA,YLDMIN(MVSIZ),YLDMAX(MVSIZ),
168 . yld(mvsiz),rate(mvsiz,2),ef(mvsiz),
169 . yfac(mvsiz,2),eps0(mvsiz),epss(mvsiz),
170 . epssmax,df1,df2,dav,emax,pext,p_air(mvsiz),
171 . eps_max,dsig ,yldelas(mvsiz),p,svm1(mvsiz),expo,hys,
172 . de,pgaz,vnew,el(mvsiz),den,pgaz0,e_air0,espe,pair0,
173 . var,kk,alpha0,et
174C-----------------------------------------------
175C USER VARIABLES INITIALIZATION
176C-----------------------------------------------
177C
178 DO i= 1, nel
179 nfunc = ipm(10,mat(i))
180 DO j=1,nfunc
181 ifunc(i,j)=ipm(10+j,mat(i))
182 ENDDO
183 ENDDO
184C
185 e0 = uparam(2)
186 aa = uparam(3)
187 epssmax = uparam(4)
188cc G = UPARAM(5)
189 nu = uparam(6)
190 nload = uparam(7)
191 nunload = uparam(8)
192 idamage = uparam(9 + 2*nfunc )
193 expo = uparam(10 + 2*nfunc )
194 hys = uparam(11 + 2*nfunc)
195 emax = uparam(12 + 2*nfunc )
196C air parameter
197 rho_air0 = uparam(14 + 2*nfunc)
198 p0 = uparam(15 + 2*nfunc)
199 gama = uparam(16 + 2*nfunc)
200 alpha0 = uparam(17 + 2*nfunc)
201 pext = uparam(18 + 2*nfunc)
202C initial foam pressure
203 pfoam = uparam(19 + 2*nfunc)
204 kk = uparam(21 + 2*nfunc)
205C function for porosity & for darcy K
206cc IFUNCR = IFUNC(1, NFUNC )
207cc IFUNCK = IFUNC(1, NFUNC -1 )
208C old stress with out gaz pressure
209 DO i=1,nel
210C
211 pair0 = uvar(i,19)
212C
213 alpha = uvar(i,21)
214 sig0xx(i) = sig0xx(i) + alpha * pair0
215 sig0yy(i) = sig0yy(i) + alpha * pair0
216 sig0zz(i) = sig0zz(i) + alpha * pair0
217C
218 rho_air = uvar(i,1)
219 e_air0 = uvar(i,2)
220 vnew = uvar(i,3)
221 dvol = uvar(i,4)
222 mu = rho_air/rho_air0
223 v0 = vnew*mu
224Cccc PGAZ0 = UVAR(I,20)
225 espe = e_air0/max(em15,v0)
226 pgaz0= (gama - one)*mu*espe
227 e_air = e_air0 - half*pgaz0*dvol
228 bb = one + half*(gama - one)*mu*dvol/max(em15,v0)
229 e_air = e_air/max(em20,bb)
230 espe = e_air/max(em15,v0)
231 pgaz =(gama - one)*mu*espe
232C
233 p_air(i) = pgaz - pext
234 p_air(i)= max(p_air(i),-pext)
235
236cc pour tester ALe P_AIR = P0*MU**GAMA
237 uvar(i,5) = -p_air(i)
238 uvar(i,6) = -p_air(i)
239 uvar(i,7) = -p_air(i)
240 uvar(i,8) = zero
241 uvar(i,9) = zero
242 uvar(i,10)= zero
243 uvar(i,20)= pgaz
244 uvar(i,2) = e_air/max(em15,vnew)
245 ef(i) = p0*gama*mu**(gama - 1)
246 fssp(i) = sqrt(ef(i)/rho_air0)
247C
248C compute of alpha & v/v0 for foam element
249C
250 ifuncr = ifunc(i,nfunc )
251 var = rho0(i)/rho(i)
252 IF( ifuncr > 0 )
253 . uvar(i,21) = alpha0*finter(ifuncr,var,npf,tf,df)
254
255 ifunck = ifunc(i,nfunc - 1)
256 IF( ifunck > 0 )
257 . uvar(i,22) = kk*finter(ifunck,var,npf,tf,df)
258C
259 uvar(i,23)= var
260 ENDDO
261C-----------------------------------------------
262 DO i=1,nel
263C-------------------
264C epst_spherique
265C-------------------
266 epst(i) = epsxx(i)**2+epsyy(i)**2 + epszz(i)**2 +
267 . half*(epsxy(i)**2+epsyz(i)**2+epszx(i)**2)
268 epst(i) = sqrt(epst(i))
269C-------------------
270C estimation elastique
271C-------------------
272 eps0(i) = uvar(i,11)
273 ENDDO
274
275C-------------------
276C CRITERE
277C-------------------
278 DO i=1,nel
279C YLD_elast
280 rate(i,1)=uparam(9)
281 yfac(i,1)=uparam(9 + nfunc )
282 IF(epst(i) >= epssmax) THEN
283 yldelas(i)=yfac(i,1)*finter(ifunc(i,1),epssmax,npf,tf,df)
284 yldelas(i)=emax*(epst(i) - epssmax) + yldelas(i)
285 ELSE
286 yldelas(i) = yfac(i,1)*finter(ifunc(i,1),epst(i),npf,tf,df)
287 ENDIF
288C yld_max
289 nc = nload
290 j1 = 1
291 DO j=2,nc-1
292 IF(abs(epsp(i)) >= abs(uparam(8 + j )))THEN
293 j1 = j
294 ENDIF
295 ENDDO
296 rate(i,1)=uparam( 8 + j1)
297 yfac(i,1)=uparam( 8 + nfunc + j1)
298 IF(epst(i) >= epssmax) THEN
299 IF(nc > 1)THEN
300 j2 = j1+1
301 rate(i,2)=uparam(8 + j2 )
302 yfac(i,2)=uparam(8 + nfunc + j2 )
303C
304 yp1 = yfac(i,1)*finter(ifunc(i,j1),epssmax,npf,tf,df1)
305 yp2 = yfac(i,2)*finter(ifunc(i,j2),epssmax,npf,tf,df2)
306C
307 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
308 yldmax(i) = max(yp1 + fac*(yp2 - yp1), em20)
309 yldmax(i) = emax*(epst(i) - epssmax) + yldmax(i)
310C
311 ELSE
312 yldmax(i) =
313 . yfac(i,1)*finter(ifunc(i,j1),epssmax,npf,tf,df)
314 yldmax(i) = emax*(epst(i) - epssmax) + yldmax(i)
315 ENDIF
316 ELSE
317 IF(nc > 1)THEN
318 j2 = j1+1
319 rate(i,2)=uparam( 8 + j2 )
320 yfac(i,2)=uparam( 8 + nfunc + j2 )
321C
322 yp1 = yfac(i,1)*finter(ifunc(i,j1),epst(i),npf,tf,df1)
323 yp2 = yfac(i,2)*finter(ifunc(i,j2),epst(i),npf,tf,df2)
324C
325 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
326 yldmax(i) = max(yp1 + fac*(yp2 - yp1), em20)
327 ELSE
328 yldmax(i) = yfac(i,1)*finter(ifunc(i,j1),epst(i),npf,tf,df)
329 ENDIF
330 ENDIF
331C yld_min
332 nc = nunload
333 j1 = 1 + nload
334 yldmin(i) = zero
335 IF(nc > 0 ) THEN
336 DO j=2,nc-1
337 IF(abs(epsp(i)) >=
338 . abs(uparam(nload + 8 + j )))THEN
339 j1 = nload + j
340 ENDIF
341 ENDDO
342C
343 rate(i,1)=uparam(8 + j1)
344 yfac(i,1)=uparam(8 + nfunc + j1)
345
346 IF(epst(i) >= epssmax) THEN
347 IF(nc > 1)THEN
348 j2 = j1+1
349 rate(i,2)=uparam(8 + j2 )
350 yfac(i,2)=uparam(8 + nfunc + j2 )
351C
352 yp1 =yfac(i,1)*finter(ifunc(i,j1),epssmax,npf,tf,df1)
353 yp2 =yfac(i,2)*finter(ifunc(i,j2),epssmax,npf,tf,df2)
354C
355 IF(yp2 < yp1 ) THEN
356 fac = (rate(i,2) - epsp(i))/(rate(i,2) - rate(i,1))
357 yldmin(i) = max(yp2 + fac*(yp1-yp2), em20)
358 ELSE
359 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
360 yldmin(i) = max(yp1 + fac*(yp2 - yp1), em20)
361 ENDIF
362 yldmin(i) = emax*(epst(i) - epssmax) + yldmin(i)
363 ELSE
364 yldmin(i)=
365 . yfac(i,1)*finter(ifunc(i,j1),epssmax,npf,tf,df)
366 yldmin(i) = emax*(epst(i) - epssmax) + yldmin(i)
367 ENDIF
368 ELSE
369 IF(nc > 1)THEN
370 j2 = j1+1
371 rate(i,2)=uparam( 8 + j2 )
372 yfac(i,2)=uparam( 8 + nfunc + j2 )
373C
374 yp1 = yfac(i,1)*finter(ifunc(i,j1),epst(i),npf,tf,df1)
375 yp2 = yfac(i,2)*finter(ifunc(i,j2),epst(i),npf,tf,df2)
376C
377 IF(yp2 < yp1 ) THEN
378 fac = (rate(i,2) - epsp(i))/(rate(i,2) - rate(i,1))
379 yldmin(i) = max(yp2 + fac*(yp1-yp2), em20)
380 ELSE
381 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
382 yldmin(i) = max(yp1 + fac*(yp2 - yp1), em20)
383 ENDIF
384 ELSE
385 yldmin(i) =
386 . yfac(i,1)*finter(ifunc(i,j1),epst(i),npf,tf,df)
387 ENDIF
388 ENDIF
389 ENDIF
390 ENDDO
391C =====================================================
392 DO i = 1,nel
393 ie_cst(i)= 0
394 delta = epst(i) - uvar(i,13)
395 IF(delta >= zero)THEN
396 yld(i) = yldmax(i)
397 iload(i) = 1
398 ELSEIF(delta < zero)THEN
399 yld(i) = yldmin(i)
400 iload(i) = -1
401 IF(idamage /= 1 )yld(i) = yldmax(i)
402 ENDIF
403C
404 e = uvar(i,12)
405 iload0 = int(uvar(i,14))
406 epss(i) = epst(i) - yld(i)/ e
407 epss(i) = max(zero, epss(i))
408 de = aa*(epss(i) - uvar(i,11))
409 IF(iload(i) == 1) THEN
410 e = e + max(de, zero)
411 IF(iload0 == -1) e= uvar(i,12)
412 uvar(i,11) = max(uvar(i,11), epss(i))
413 ELSE
414 e = e + min(de ,zero)
415 IF(iload0 == 1) e= uvar(i,12)
416 uvar(i,11) = min(epss(i),uvar(i,11))
417 ENDIF
418 e = min(e, emax)
419 e = max(e, e0)
420 uvar(i,12) = e
421C==================================================
422 aa1(i) = e*(one-nu)/(one + nu)/(one - two*nu)
423 aa2(i) = aa1(i)*nu/(one - nu)
424 g(i) =half*e/(one + nu)
425C ----
426 signxx(i)= aa1(i)*depsxx(i) + aa2(i)*(depsyy(i) + depszz(i))
427 signyy(i)= aa1(i)*depsyy(i) + aa2(i)*(depsxx(i) + depszz(i))
428 signzz(i)= aa1(i)*depszz(i) + aa2(i)*(depsxx(i) + depsyy(i))
429 signxy(i)= g(i) *depsxy(i)
430 signyz(i)= g(i) *depsyz(i)
431 signzx(i)= g(i) *depszx(i)
432 dsig = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
433 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
434 dsig =sqrt(dsig)
435C estimate stress
436 signxx(i)=sig0xx(i) + aa1(i)*depsxx(i)
437 . + aa2(i)*(depsyy(i) + depszz(i))
438 signyy(i)=sig0yy(i) + aa1(i)*depsyy(i)
439 . + aa2(i)*(depsxx(i) + depszz(i))
440 signzz(i)=sig0zz(i) + aa1(i)*depszz(i)
441 . + aa2(i)*(depsxx(i) + depsyy(i))
442 signxy(i)=sig0xy(i) + g(i) *depsxy(i)
443 signyz(i)=sig0yz(i) + g(i) *depsyz(i)
444 signzx(i)=sig0zx(i) + g(i) *depszx(i)
445c
446
447 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
448 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
449 svm =sqrt(svm2)
450C sound speed
451 et = aa1(i)
452 soundsp(i) = sqrt((et + ef(i))/rho0(i))
453 viscmax(i) = zero
454C
455 IF(idamage == 1 ) THEN
456 IF(svm >= yldmax(i) )THEN
457 yld(i) = yldmax(i)
458cc IF(DELTA < ZERO .AND. DSIG > YLDMIN(I).AND. DSIG > SVM )
459cc . YLD(I) = YLDMIN(I)
460 IF(delta < zero ) yld(i) = yldmin(i)
461 ELSEIF(svm <= yldmin(i)) THEN
462 yld(i) = yldmin(i)
463C
464 ELSE
465 yld(i) = svm
466 ie_cst(i) = 1
467 IF(delta < zero .AND. dsig > yldmin(i) .AND. dsig > svm)THEN
468 yld(i) = yldmin(i)
469 ie_cst(i) = 0
470 ENDIF
471 ENDIF
472 ELSE
473 yld(i) = yldmax(i)
474 IF(delta > zero .AND. svm < yldmax(i))yld(i)=svm
475cc IF(DELTA < ZERO .AND. SVM > YLDMAX(I)
476cc . .AND. DSIG < SVM )YLD(I) = SVM
477 uvar(i,17) = uvar(i,17) +
478 . half*(yld(i ) + uvar(i,15))*delta
479 uvar(i,17) = max(zero, uvar(i,17))
480 uvar(i,18) = max(uvar(i,18), uvar(i,17))
481 ENDIF
482C
483 ENDDO
484C-------------------
485C projection spherique
486C-------------------
487 IF(idamage == 1 ) THEN
488 DO i=1,nel
489 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
490 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
491 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
492 signxy(i)= g(i) *epsxy(i)
493 signyz(i)= g(i) *epsyz(i)
494 signzx(i)= g(i) *epszx(i)
495C
496 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
497 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
498 svm =sqrt(svm2)
499 r = yld(i)/max(em20,svm)
500 signxx(i)=signxx(i)*r
501 signyy(i)=signyy(i)*r
502 signzz(i)=signzz(i)*r
503 signxy(i)=signxy(i)*r
504 signyz(i)=signyz(i)*r
505 signzx(i)=signzx(i)*r
506C
507 IF(ie_cst(i) == 1) THEN
508 iload0 = int(uvar(i,14))
509 IF(iload0 /= iload(i))THEN
510 iload(i) = iload0
511 uvar(i,11) = eps0(i)
512 ENDIF
513 ENDIF
514 uvar(i,13) = epst(i)
515 uvar(i,14) = iload(i)
516 uvar(i,15) = yld(i)
517 uvar(i,16) = epsp(i)
518 ENDDO
519 ELSEIF(idamage == 2 ) THEN
520 DO i=1,nel
521C
522 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
523 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
524 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
525 signxy(i)= g(i) *epsxy(i)
526 signyz(i)= g(i) *epsyz(i)
527 signzx(i)= g(i) *epszx(i)
528C
529 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
530 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
531 svm =sqrt(svm2)
532 r = yld(i)/max(em20,svm)
533 signxx(i)=signxx(i)*r
534 signyy(i)=signyy(i)*r
535 signzz(i)=signzz(i)*r
536 signxy(i)=signxy(i)*r
537 signyz(i)=signyz(i)*r
538 signzx(i)=signzx(i)*r
539C
540 IF(iload(i) == -1) THEN
541 r = yldmin(i)/max(em20,yldelas(i))
542 signxx(i)=signxx(i)*r
543 signyy(i)=signyy(i)*r
544 signzz(i)=signzz(i)*r
545 signxy(i)=signxy(i)*r
546 signyz(i)=signyz(i)*r
547 signzx(i)=signzx(i)*r
548 ENDIF
549C
550 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
551 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
552 svm =sqrt(svm2)
553 uvar(i,13) = epst(i)
554 uvar(i,14) = iload(i)
555 uvar(i,15) = yld(i)
556 uvar(i,16) = epsp(i)
557 ENDDO
558 ELSEIF(idamage == 3) THEN
559 DO i=1,nel
560C
561 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
562 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
563 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
564 signxy(i)= g(i) *epsxy(i)
565 signyz(i)= g(i) *epsyz(i)
566 signzx(i)= g(i) *epszx(i)
567C
568C
569 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
570 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
571 svm =sqrt(svm2)
572 r = (yld(i)/max(em20,svm))
573 signxx(i)=signxx(i)*r
574 signyy(i)=signyy(i)*r
575 signzz(i)=signzz(i)*r
576 signxy(i)=signxy(i)*r
577 signyz(i)=signyz(i)*r
578 signzx(i)=signzx(i)*r
579C
580 IF(iload(i) == -1) THEN
581 r = one - (uvar(i,17)/max(em20,uvar(i,18)))**expo
582 r = one - (one - hys)*r
583 signxx(i)=signxx(i)*r
584 signyy(i)=signyy(i)*r
585 signzz(i)=signzz(i)*r
586 signxy(i)=signxy(i)*r
587 signyz(i)=signyz(i)*r
588 signzx(i)=signzx(i)*r
589 ENDIF
590C
591 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
592 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
593 svm =sqrt(svm2)
594 uvar(i,13) = epst(i)
595 uvar(i,14) = iload(i)
596 uvar(i,15) = yld(i)
597 uvar(i,16) = epsp(i)
598 ENDDO
599 ENDIF
600C------------------------------------
601 IF(time == 0 .AND. pfoam /= zero)THEN
602 DO i=1,nel
603 signxx(i)= -pfoam
604 signyy(i)= -pfoam
605 signzz(i)= -pfoam
606 ENDDO
607 ELSE
608
609 IF(pfoam /= zero) THEN
610 DO i=1,nel
611 IF( epst(i) == zero ) THEN
612 signxx(i)= -pfoam
613 signyy(i)= -pfoam
614 signzz(i)= -pfoam
615 ENDIF
616 ENDDO
617 ENDIF
618 ENDIF
619C
620 DO i=1,nel
621 alpha = uvar(i,21)
622 signxx(i) = signxx(i) - alpha * p_air(i)
623 signyy(i) = signyy(i) - alpha * p_air(i)
624 signzz(i) = signzz(i) - alpha * p_air(i)
625 uvar(i,19) = p_air(i)
626 ENDDO
627c-----------
628 RETURN
629 END
630C
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps77
subroutine sigeps77(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, sig0xx, sig0yy, sig0zz, sig0xy, sig0yz, sig0zx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, fssp, viscmax, uvar, off, ngl, pm, ipm, mat, epsp, p_air)
Definition sigeps77.F:42