OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps44.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!|| sigeps44 ../engine/source/materials/mat/mat044/sigeps44.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
28!||--- calls -----------------------------------------------------
29!|| mstrain_rate ../engine/source/materials/mat_share/mstrain_rate.F
30!|| vinter ../engine/source/tools/curve/vinter.F
31!||====================================================================
32 SUBROUTINE sigeps44(
33 1 NEL ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,NPF ,
34 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,RHO ,
35 3 VOLUME ,EINT ,IEOS ,DPDM ,
36 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
37 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
39 7 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 9 SIGVXX ,SIGVYY ,SIGVZZ ,SIGVXY ,SIGVYZ ,SIGVZX ,
42 A SOUNDSP,VISCMAX,UVAR ,OFF ,NGL ,IPT ,
43 B IPM ,MAT ,EPSP ,IPLA ,SIGY ,PLA ,
44 C DPLA1 ,AMU ,ISRATE ,ASRATE ,NVARTMP ,VARTMP ,
45 D ET )
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50C-----------------------------------------------
51C G l o b a l P a r a m e t e r s
52C-----------------------------------------------
53#include "mvsiz_p.inc"
54#include "param_c.inc"
55#include "scr17_c.inc"
56C---------+---------+---+---+--------------------------------------------
57C VAR | SIZE |TYP| RW| DEFINITION
58C---------+---------+---+---+--------------------------------------------
59C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
60C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
61C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
62C---------+---------+---+---+--------------------------------------------
63C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
64C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
65C NPF | * | I | R | FUNCTION ARRAY
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 | NEL | F | R | INITIAL DENSITY
72C RHO | NEL | F | R | DENSITY
73C VOLUME | NEL | F | R | VOLUME
74C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
75C EPSPXX | NEL | F | R | STRAIN RATE XX
76C EPSPYY | NEL | F | R | STRAIN RATE YY
77C ... | | | |
78C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
79C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
80C ... | | | |
81C EPSXX | NEL | F | R | STRAIN XX
82C EPSYY | NEL | F | R | STRAIN YY
83C ... | | | |
84C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
85C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
86C ... | | | |
87C---------+---------+---+---+--------------------------------------------
88C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
89C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
90C ... | | | |
91C SIGVXX | NEL | F | W | VISCOUS STRESS XX
92C SIGVYY | NEL | F | W | VISCOUS STRESS YY
93C ... | | | |
94C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
95C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
96C---------+---------+---+---+--------------------------------------------
97C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
98C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
99C---------+---------+---+---+--------------------------------------------
100 INTEGER ,INTENT(IN) :: IEOS
101 INTEGER NEL, NUPARAM, NUVAR,IPT,NVARTMP,
102 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*), IPLA,IEXPAN,
103 . ISRATE
104 my_real
105 . TIME,TIMESTEP,UPARAM(NUPARAM),
106 . RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
107 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
108 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
109 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
110 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
111 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
112 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
113 . sigoxx(nel),sigoyy(nel),sigozz(nel),
114 . sigoxy(nel),sigoyz(nel),sigozx(nel),
115 . epsp(nel),amu(nel), asrate
116 my_real, DIMENSION(MVSIZ) ,INTENT(IN) :: dpdm
117C-----------------------------------------------
118C O U T P U T A r g u m e n t s
119C-----------------------------------------------
120 my_real
121 . signxx(nel),signyy(nel),signzz(nel),
122 . signxy(nel),signyz(nel),signzx(nel),
123 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
124 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
125 . soundsp(nel),viscmax(nel),sigy(nel),
126 . pla(nel),dpla1(mvsiz)
127 INTEGER :: VARTMP(NEL,NVARTMP)
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)
133 my_real, DIMENSION(NEL), INTENT(INOUT) ::
134 . et
135C-------------------------
136C variables non utilisees (fonctions utilisateur)
137C-------------------------
138 INTEGER NPF(*),MFUNC,KFUNC(MFUNC)
139 my_real
140 . tf(*)
141C-----------------------------------------------
142C L o c a l V a r i a b l e s
143C-----------------------------------------------
144 INTEGER I,IADBUF,ICC,VFLAG,IPOS(MVSIZ),IAD(MVSIZ),ILEN(MVSIZ),IDEV
145 my_real
146 . e,nu,p,dav,vm,r,epst,e1,e2,e3,e4,e5,e6,c,cd,d,
147 . y,yp,e42,e52,e62,epst2,rq,bulk,
148 . g,g2,g3,epsgm,epmax,smax(mvsiz),fail(mvsiz),
149 . epsr1,epsr2,yld(mvsiz),h(mvsiz),dfdpla(mvsiz),
150 . ca,cb,cn,cc,cp,dpla,yscale,dpdt
151C=======================================================================
152C-----------------------------------------------
153C MATERIAL PARAMETERS
154C-----------------------------------------------
155 e = uparam(1)
156 nu = uparam(2)
157 ca = uparam(3)
158 smax(1:nel) = uparam(4)
159 epmax = uparam(5)
160 epsr1 = uparam(6)
161 epsr2 = uparam(7)
162 cb = uparam(8)
163 cn = uparam(9)
164 icc = nint(uparam(10))
165 cc = uparam(11)
166 cp = uparam(12)
167 g = uparam(16)
168 g2 = uparam(17)
169 g3 = uparam(18)
170 bulk = uparam(19)
171 epsgm = uparam(22)
172 vflag = nint(uparam(23))
173 yscale = uparam(24)
174C-----------------------------------------------
175C TOTAL OR DEVIATORIC STRAIN-RATE COMPUTATION
176C-----------------------------------------------
177 idev = vflag-2
178 IF ((vflag == 2) .OR. (vflag == 3)) THEN
179 CALL mstrain_rate(nel ,israte ,asrate ,epsp ,idev ,
180 . epspxx ,epspyy ,epspzz ,epspxy ,epspyz ,
181 . epspzx )
182 ENDIF
183C-----------------------------------------------
184C DEVIATORIC STRESS ESTIMATION
185C-----------------------------------------------
186 DO i=1,nel
187 p = -(sigoxx(i)+sigoyy(i)+sigozz(i))* third
188 dav = (depsxx(i)+depsyy(i)+depszz(i)) * third
189 signxx(i) = sigoxx(i)+p+g2*(depsxx(i)-dav)
190 signyy(i) = sigoyy(i)+p+g2*(depsyy(i)-dav)
191 signzz(i) = sigozz(i)+p+g2*(depszz(i)-dav)
192 signxy(i) = sigoxy(i)+g *depsxy(i)
193 signyz(i) = sigoyz(i)+g *depsyz(i)
194 signzx(i) = sigozx(i)+g *depszx(i)
195 viscmax(i) = zero
196 dpla1(i) = zero
197 ENDDO
198C
199C-----------------------------------------------
200C STRAIN principal 1, 4 newton iterations
201C-----------------------------------------------
202 DO i=1,nel
203C
204 dav = (epsxx(i)+epsyy(i)+epszz(i)) * third
205 e1 = epsxx(i) - dav
206 e2 = epsyy(i) - dav
207 e3 = epszz(i) - dav
208 e4 = half*epsxy(i)
209 e5 = half*epsyz(i)
210 e6 = half*epszx(i)
211C -y = (e1-x)(e2-x)(e3-x)
212C - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
213C + 2e4 e5 e6
214C e1 + e2 + e3 = 0 => terme en x^2 = 0
215C y = x^3 + c x + d
216c yp= 3 x^2 + c
217 e42 = e4*e4
218 e52 = e5*e5
219 e62 = e6*e6
220 c = - e1*e1 - e2*e2 - e3*e3 - e42 - e52 - e62
221 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
222 & - two*e4*e5*e6
223 cd = c*third
224 epst = sqrt(-cd)
225 epst2 = epst * epst
226 y = (epst2 + c)* epst + d
227 IF(abs(y)>em8)THEN
228 epst = onep75 * epst
229 epst2 = epst * epst
230 y = (epst2 + c)* epst + d
231 yp = three*epst2 + c
232 IF(yp/=zero)epst = epst - y/yp
233 epst2 = epst * epst
234 y = (epst2 + c)* epst + d
235 yp = three*epst2 + c
236 IF(yp/=zero)epst = epst - y/yp
237 epst2 = epst * epst
238 y = (epst2 + c)* epst + d
239 yp = three*epst2 + c
240 IF(yp/=zero)epst = epst - y/yp
241 epst2 = epst * epst
242 y = (epst2 + c)* epst + d
243 yp = three*epst2 + c
244 IF(yp/=zero)epst = epst - y/yp
245 epst = epst + dav
246 ENDIF
247C
248C--- tension failure
249C
250 fail(i) = max(zero,min(one,
251 . (epsr2-epst)/(epsr2-epsr1) ))
252C
253 ENDDO
254C-----------------------------------------------
255C STRAIN RATE EFFECT, CURRENT YIELD & HARDENING
256C-----------------------------------------------
257 IF (mfunc > 0) THEN
258 ipos(1:nel) = vartmp(1:nel,1)
259 iad(1:nel) = npf(kfunc(1)) / 2 + 1
260 ilen(1:nel) = npf(kfunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
261 CALL vinter(tf,iad,ipos,ilen,nel,pla,dfdpla,yld)
262 vartmp(1:nel,1) = ipos(1:nel)
263 ENDIF
264
265 IF (ipla /= 2) THEN
266 DO i=1,nel
267 rq = one
268 IF (cc /= zero) rq = one + (cc*epsp(i))**cp
269 IF (icc == 1) smax(i) = smax(i)*rq
270 IF (pla(i) > zero) THEN
271 IF ((mfunc > 0) .AND. (ca == zero)) THEN
272 yld(i) = yscale * yld(i) * rq
273 h(i) = fail(i) * (yscale*dfdpla(i)*rq)
274 ELSEIF ((mfunc > 0) .AND. (ca > zero)) THEN
275 yld(i) = yscale * yld(i)
276 . + (ca + cb*pla(i)**cn) * (rq-one)
277 h(i) = fail(i) * (yscale*dfdpla(i) + cn*cb*(rq-one)/(pla(i)**(one-cn)))
278 ELSE
279 yld(i) = (ca + cb*pla(i)**cn) * rq
280 h(i) = fail(i)*cn*cb*rq/(pla(i)**(one-cn))
281 ENDIF
282 yld(i) = fail(i) * min(smax(i),yld(i))
283 ELSE
284 IF ((mfunc > 0) .AND. (ca == zero)) THEN
285 yld(i) = fail(i) * yscale * yld(i) * rq
286 ELSEIF ((mfunc > 0) .AND. (ca > zero)) THEN
287 yld(i) = fail(i) * (yscale * yld(i) + ca * (rq-one))
288 ELSE
289 yld(i) = fail(i) * ca * rq
290 ENDIF
291 h(i) = e
292 ENDIF
293 IF (pla(i) >= epsgm) THEN
294 yld(i) = fail(i) * smax(i)
295 h(i) = zero
296 ENDIF
297 sigy(i) = yld(i)
298 IF (pla(i) >= epmax) yld(i) = zero
299 ENDDO
300 ELSE
301 DO i=1,nel
302 rq = one
303 IF (cc /= zero) rq = one + (cc*epsp(i))**cp
304 IF (icc == 1) smax(i) = smax(i)*rq
305 IF (pla(i) > zero) THEN
306 IF ((mfunc > 0) .AND. (ca == zero)) THEN
307 yld(i) = yscale * yld(i) * rq
308 h(i) = fail(i)*(yscale*dfdpla(i)*rq)
309 ELSEIF ((mfunc > 0) .AND. (ca > zero)) THEN
310 yld(i) = yscale * yld(i)
311 . + (ca + cb*pla(i)**cn) * (rq-one)
312 h(i) = fail(i)* (yscale*dfdpla(i) + cn*cb*(rq-one)/(pla(i)**(one-cn)))
313 ELSE
314 yld(i) = (ca + cb*pla(i)**cn) * rq
315 h(i) = fail(i)*cn*cb*rq/(pla(i)**(one-cn))
316 ENDIF
317 yld(i) = fail(i) * min(smax(i),yld(i))
318 ELSE
319 IF ((mfunc > 0) .AND. (ca == zero)) THEN
320 yld(i) = fail(i) * yscale * yld(i) * rq
321 ELSEIF ((mfunc > 0) .AND. (ca > zero)) THEN
322 yld(i) = fail(i) * (yscale * yld(i) + ca * (rq-one))
323 ELSE
324 yld(i) = fail(i) * ca * rq
325 ENDIF
326 h(i) = e
327 ENDIF
328 IF (pla(i) >= epsgm) THEN
329 yld(i) = fail(i) * smax(i)
330 yld(i) = zero
331 ENDIF
332 sigy(i) = yld(i)
333 IF (pla(i) >= epmax) yld(i)=zero
334 ENDDO
335 ENDIF
336C
337C-----------------------------------------------
338C VON MISES & RADIAL RETURN
339C-----------------------------------------------
340 IF (ipla == 0) THEN
341 DO i = 1,nel
342 vm = half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
343 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
344 vm = sqrt(three*vm)
345 r = min(one,yld(i)/ max(vm,em20))
346 signxx(i) = signxx(i)*r
347 signyy(i) = signyy(i)*r
348 signzz(i) = signzz(i)*r
349 signxy(i) = signxy(i)*r
350 signyz(i) = signyz(i)*r
351 signzx(i) = signzx(i)*r
352 pla(i) = pla(i) + (one -r)*vm/max(g3+h(i),em20)
353 dpla1(i) = (one -r)*vm/max(g3+h(i),em20)
354 ENDDO
355 ELSEIF (ipla == 2) THEN
356 DO i = 1,nel
357 vm = half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
358 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
359 vm = sqrt(three*vm)
360 r = min(one,yld(i)/ max(vm,em20))
361 signxx(i) = signxx(i)*r
362 signyy(i) = signyy(i)*r
363 signzz(i) = signzz(i)*r
364 signxy(i) = signxy(i)*r
365 signyz(i) = signyz(i)*r
366 signzx(i) = signzx(i)*r
367 pla(i) = pla(i) + (one -r)*vm/max(g3,em20)
368 dpla1(i) = (one -r)*vm/max(g3,em20)
369 ENDDO
370 ELSEIF (ipla == 1) THEN
371 DO i = 1,nel
372 vm = half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
373 1 +signxy(i)**2+signyz(i)**2+signzx(i)**2
374 vm = sqrt(three*vm)
375 r = min(one,yld(i)/ max(vm,em20))
376C plastic strain increment.
377 dpla = (one - r)*vm/max(g3+h(i),em20)
378C actual yield stress
379 yld(i) = yld(i) + dpla*h(i)
380 IF (pla(i) >= epmax) yld(i)=zero
381 r = min(one,yld(i)/ max(vm,em20))
382 signxx(i) = signxx(i)*r
383 signyy(i) = signyy(i)*r
384 signzz(i) = signzz(i)*r
385 signxy(i) = signxy(i)*r
386 signyz(i) = signyz(i)*r
387 signzx(i) = signzx(i)*r
388 pla(i) = pla(i) + dpla
389 dpla1(i) = dpla
390 ENDDO
391 ENDIF
392C----------------------------------
393 ! update total stress tensor
394 IF (ieos == 0) THEN
395 DO i = 1, nel
396 ! P = BULK(I) * (RHO(I)/RHO0(I) - ONE)
397 p = bulk * amu(i) ! AMU is computed by taking into account thermal strain
398 signxx(i) = signxx(i) - p
399 signyy(i) = signyy(i) - p
400 signzz(i) = signzz(i) - p
401 signxy(i) = signxy(i)
402 signyz(i) = signyz(i)
403 signzx(i) = signzx(i)
404 soundsp(i) = sqrt((bulk + four*g/three)/rho0(i))
405 ENDDO
406 ELSE
407 ! if EOS is used, material law calculates only deviatoric stress tensor
408 ! sound speed depends on pressure derivative over volume change
409 ! calculated in EOS
410 DO i = 1, nel
411 soundsp(i) = sqrt((dpdm(i) + four*g/three)/rho0(i))
412 ENDDO
413 END IF
414C----------------------------------
415 IF (vflag == 1) THEN ! plastic strain rate filtering
416 DO i=1,nel
417 dpdt = dpla1(i) / max(em20,timestep)
418 epsp(i) = asrate * dpdt + (one - asrate) * epsp(i)
419 ENDDO
420 ENDIF
421C
422 DO i=1,nel
423 sigy(i) = max(sigy(i),yld(i))
424 IF (dpla1(i) > zero) THEN
425 et(i) = h(i) / (h(i) + e)
426 ELSE
427 et(i) = one
428 ENDIF
429 ENDDO
430C-----------
431 RETURN
432 END
433C
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
mstrain_rate
subroutine mstrain_rate(nel, israte, asrate, epsd, idev, ep1, ep2, ep3, ep4, ep5, ep6)
Definition mstrain_rate.F:33
sigeps44
subroutine sigeps44(nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, ieos, dpdm, 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, ipla, sigy, pla, dpla1, amu, israte, asrate, nvartmp, vartmp, et)
Definition sigeps44.F:46
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72