OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps48.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"
#include "scr17_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps48 (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, ipt, ipm, mat, epsp, sigy, defp, dpla1, amu)

Function/Subroutine Documentation

◆ sigeps48()

subroutine sigeps48 ( integer nel,
integer nuparam,
integer nuvar,
integer mfunc,
integer, dimension(mfunc) kfunc,
integer, dimension(*) 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,
integer, dimension(nel) ngl,
integer ipt,
integer, dimension(npropmi,*) ipm,
integer, dimension(nel) mat,
epsp,
sigy,
defp,
dpla1,
amu )

Definition at line 29 of file sigeps48.F.

42C-----------------------------------------------
43C I m p l i c i t T y p e s
44C-----------------------------------------------
45#include "implicit_f.inc"
46C-----------------------------------------------
47C G l o b a l P a r a m e t e r s
48C-----------------------------------------------
49#include "mvsiz_p.inc"
50C---------+---------+---+---+--------------------------------------------
51C VAR | SIZE |TYP| RW| DEFINITION
52C---------+---------+---+---+--------------------------------------------
53C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
54C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
55C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
56C---------+---------+---+---+--------------------------------------------
57C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
58C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
59C NPF | * | I | R | FUNCTION ARRAY
60C TF | * | F | R | FUNCTION ARRAY
61C---------+---------+---+---+--------------------------------------------
62C TIME | 1 | F | R | CURRENT TIME
63C TIMESTEP| 1 | F | R | CURRENT TIME STEP
64C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
65C RHO0 | NEL | F | R | INITIAL DENSITY
66C RHO | NEL | F | R | DENSITY
67C VOLUME | NEL | F | R | VOLUME
68C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
69C EPSPXX | NEL | F | R | STRAIN RATE XX
70C EPSPYY | NEL | F | R | STRAIN RATE YY
71C ... | | | |
72C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
73C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
74C ... | | | |
75C EPSXX | NEL | F | R | STRAIN XX
76C EPSYY | NEL | F | R | STRAIN YY
77C ... | | | |
78C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
79C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
80C ... | | | |
81C---------+---------+---+---+--------------------------------------------
82C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
83C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
84C ... | | | |
85C SIGVXX | NEL | F | W | VISCOUS STRESS XX
86C SIGVYY | NEL | F | W | VISCOUS STRESS YY
87C ... | | | |
88C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
89C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
90C---------+---------+---+---+--------------------------------------------
91C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
92C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
93C---------+---------+---+---+--------------------------------------------
94#include "param_c.inc"
95#include "scr17_c.inc"
96C
97 INTEGER NEL, NUPARAM, NUVAR,IPT,
98 . IPM(NPROPMI,*),NGL(NEL),MAT(NEL)
99 my_real time,timestep,uparam(*),
100 . rho(nel),rho0(nel),volume(nel),eint(nel),
101 . epspxx(nel),epspyy(nel),epspzz(nel),
102 . epspxy(nel),epspyz(nel),epspzx(nel),
103 . depsxx(nel),depsyy(nel),depszz(nel),
104 . depsxy(nel),depsyz(nel),depszx(nel),
105 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,
106 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
107 . sigoxx(nel),sigoyy(nel),sigozz(nel),
108 . sigoxy(nel),sigoyz(nel),sigozx(nel),
109 . epsp(nel),amu(nel)
110C-----------------------------------------------
111C O U T P U T A r g u m e n t s
112C-----------------------------------------------
113 my_real
114 . signxx(nel),signyy(nel),signzz(nel),
115 . signxy(nel),signyz(nel),signzx(nel),
116 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
117 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
118 . soundsp(nel),viscmax(nel),sigy(nel),defp(nel),
119 . dpla1(mvsiz)
120C-----------------------------------------------
121C I N P U T O U T P U T A r g u m e n t s
122C-----------------------------------------------
123 my_real
124 . uvar(nel,nuvar), off(nel)
125C-------------------------
126C Variables not used (user functions)
127C-------------------------
128 INTEGER NPF(*),MFUNC,KFUNC(MFUNC)
129 my_real
130 . tf(*)
131C EXTERNAL FINTER
132C-----------------------------------------------
133C L o c a l V a r i a b l e s
134C-----------------------------------------------
135 INTEGER I,IADBUF,MX
136 my_real
137 . e,nu,p,dav,vm,r,dpla,epst,e1,e2,e3,e4,e5,e6,c,cd,d,
138 . y,yp,e42,e52,e62,epst2,pa,pb,pc,pda,pdb,yy,
139 . c1(mvsiz),c2(mvsiz),c3(mvsiz),
140 . g(mvsiz),g2(mvsiz),g3(mvsiz),epsgm(mvsiz),
141 . epmax(mvsiz),smax(mvsiz),pla(mvsiz),fail(mvsiz),
142 . epsr1(mvsiz),epsr2(mvsiz),yld(mvsiz),h(mvsiz),
143 . ca(mvsiz),cb(mvsiz),cn(mvsiz),cc(mvsiz),co(mvsiz),
144 . cm(mvsiz),ce(mvsiz),ck(mvsiz),eps0(mvsiz)
145C==================================================================|
146 IF(time==zero)THEN
147 DO i=1,nel
148 uvar(i,1)=zero
149 ENDDO
150 ENDIF
151C
152C---
153 mx = mat(1)
154 iadbuf = ipm(7,mx)-1
155 DO i=1,nel
156 e = uparam(iadbuf+1)
157 nu = uparam(iadbuf+2)
158 ca(i) = uparam(iadbuf+3)
159 smax(i) = uparam(iadbuf+4)
160 epmax(i)= uparam(iadbuf+5)
161 epsr1(i)= uparam(iadbuf+6)
162 epsr2(i)= uparam(iadbuf+7)
163 cb(i) = uparam(iadbuf+8)
164 cn(i) = uparam(iadbuf+9)
165 cc(i) = uparam(iadbuf+10)
166 co(i) = uparam(iadbuf+11)
167 cm(i) = uparam(iadbuf+12)
168 g(i) = uparam(iadbuf+16)
169 g2(i) = uparam(iadbuf+17)
170 g3(i) = uparam(iadbuf+18)
171 c1(i) = uparam(iadbuf+19)
172 c2(i) = uparam(iadbuf+20)
173 c3(i) = uparam(iadbuf+21)
174 epsgm(i)= uparam(iadbuf+22)
175 eps0(i) = uparam(iadbuf+23)
176 ce(i) = uparam(iadbuf+24)
177 ck(i) = uparam(iadbuf+25)
178 ENDDO
179C
180C--- DEVIATORIC STRESS ESTIMATE
181C
182 DO i=1,nel
183 pla(i) = uvar(i,1)
184 p = -(sigoxx(i)+sigoyy(i)+sigozz(i))* third
185 dav = (depsxx(i)+depsyy(i)+depszz(i))* third
186 signxx(i)=sigoxx(i) + p + g2(i)*(depsxx(i)-dav)
187 signyy(i)=sigoyy(i) + p + g2(i)*(depsyy(i)-dav)
188 signzz(i)=sigozz(i) + p + g2(i)*(depszz(i)-dav)
189 signxy(i)=sigoxy(i) + g(i)*depsxy(i)
190 signyz(i)=sigoyz(i) + g(i)*depsyz(i)
191 signzx(i)=sigozx(i) + g(i)*depszx(i)
192C
193 soundsp(i) = sqrt((c1(i)+four*g(i)/three)/rho0(i))
194 viscmax(i) = zero
195 dpla1(i) = zero
196 ENDDO
197C__________
198C
199C--- STRAIN principal 1, 4 newton iterations
200C
201 DO i=1,nel
202C
203 dav = (epsxx(i)+epsyy(i)+epszz(i))*third
204 e1 = epsxx(i) - dav
205 e2 = epsyy(i) - dav
206 e3 = epszz(i) - dav
207 e4 = half*epsxy(i)
208 e5 = half*epsyz(i)
209 e6 = half*epszx(i)
210C -y = (e1-x)(e2-x)(e3-x)
211C - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
212C + 2e4 e5 e6
213C e1 + e2 + e3 = 0 => terme en x^2 = 0
214C y = x^3 + c x + d
215c yp= 3 x^2 + c
216 e42 = e4*e4
217 e52 = e5*e5
218 e62 = e6*e6
219 c = - e1*e1 - e2*e2 - e3*e3 - e42 - e52 - e62
220 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
221 & - two*e4*e5*e6
222 cd = c*third
223 epst = sqrt(-cd)
224 epst2 = epst * epst
225 y = (epst2 + c)* epst + d
226 IF(abs(y)>em8)THEN
227 epst = onep75 * epst
228 epst2 = epst * epst
229 y = (epst2 + c)* epst + d
230 yp = three*epst2 + c
231 IF(yp/=zero)epst = epst - y/yp
232 epst2 = epst * epst
233 y = (epst2 + c)* epst + d
234 yp = three*epst2 + c
235 IF(yp/=zero)epst = epst - y/yp
236 epst2 = epst * epst
237 y = (epst2 + c)* epst + d
238 yp = three*epst2 + c
239 IF(yp/=zero)epst = epst - y/yp
240 epst2 = epst * epst
241 y = (epst2 + c)* epst + d
242 yp = 3*epst2 + c
243 IF(yp/=zero)epst = epst - y/yp
244 epst = epst + dav
245 ENDIF
246C
247C--- tension failure
248C
249 fail(i) = max(zero,min(one,
250 . (epsr2(i)-epst)/(epsr2(i)-epsr1(i)) ))
251C
252 ENDDO
253C
254C--- STRAIN RATE EFFECT, CURRENT YIELD & HARDENING
255C---------------------
256C Law formulation:
257C sig = ca*SigY + cb*epla^cn + (cc-co*epla^cm)*ln(epsp/eps0) + ce*epsp^ck
258C---------------------
259 DO i=1,nel
260 IF(pla(i)<=zero) THEN
261 pa=ca(i)
262 ELSEIF(pla(i)>epmax(i)) THEN
263 pa=ca(i)+cb(i)*epmax(i)**cn(i)
264 ELSE
265 pa=ca(i)+cb(i)*pla(i)**cn(i)
266 ENDIF
267 IF(epsp(i)<=eps0(i)) THEN
268 pb=zero
269 ELSEIF(pla(i)<=zero) THEN
270 pb=cc(i)*log(epsp(i)/eps0(i))
271 ELSE
272 pb=(cc(i)-co(i)*pla(i)**cm(i))*log(epsp(i)/eps0(i))
273 ENDIF
274 IF(epsp(i)<=zero) THEN
275 pc=zero
276 ELSE
277 pc=ce(i)*epsp(i)**ck(i)
278 ENDIF
279C-----
280 IF(pla(i)>zero. and .cn(i)>=one) THEN
281 pda = cb(i)*cn(i)*pla(i)**(cn(i)- one)
282 ELSEIF(pla(i)>zero. and .cn(i)<one)THEN
283 pda = cb(i)*cn(i)*pla(i)**(one-cn(i))
284 ELSE
285 pda = e
286 ENDIF
287 IF(pla(i)<=zero. or .epsp(i)<=eps0(i)) THEN
288 pdb = zero
289 ELSEIF(cm(i)>=one) THEN
290 pdb = co(i)*cm(i)*pla(i)**(cm(i)- one)*log(epsp(i)/eps0(i))
291 ELSE
292 pdb = co(i)*cm(i)*pla(i)**(one - cm(i))*log(epsp(i)/eps0(i))
293 ENDIF
294C-----
295 yy = pa + pb + pc
296 yld(i)= min(smax(i)+pc, yy)
297 h(i) = pda + pdb
298 IF (yld(i)<yy) h(i) = zero
299 yld(i)= fail(i)*yld(i)
300 h(i) = fail(i)*h(i)
301 sigy(i)=yld(i)
302 IF(pla(i)>epmax(i)) yld(i)=zero
303 ENDDO
304C
305C-------------------
306C VON MISES & RADIAL RETURN
307C-------------------
308 DO i=1,nel
309 vm =half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
310 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
311 vm =sqrt(three*vm)
312 r = min(one,yld(i)/ max(vm,em20))
313C plastic strain increment.
314 dpla=(one -r)*vm/max(g3(i)+h(i),em20)
315C actual yield stress.
316 yld(i)=yld(i)+dpla*h(i)
317 IF(pla(i)>epmax(i)) yld(i)=zero
318 r = min(one,yld(i)/ max(vm,em20))
319C
320c P = C1(I) * (RHO(I)/RHO0(I)- ONE)
321 p = c1(i) * amu(i)
322 signxx(i)=signxx(i)*r-p
323 signyy(i)=signyy(i)*r-p
324 signzz(i)=signzz(i)*r-p
325 signxy(i)=signxy(i)*r
326 signyz(i)=signyz(i)*r
327 signzx(i)=signzx(i)*r
328 pla(i)=pla(i) + dpla
329 uvar(i,1)=pla(i)
330 dpla1(i) = dpla
331 ENDDO
332C
333 DO i=1,nel
334 sigy(i)=max(sigy(i),yld(i))
335 defp(i)=uvar(i,1)
336 ENDDO
337C________________________________________________________
338 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21