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 31 of file sigeps48.F.

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