OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps84.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!|| sigeps84 ../engine/source/materials/mat/mat084/sigeps84.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||====================================================================
28 SUBROUTINE sigeps84(
29 1 NEL, NUPARAM ,NUVAR , JLAG,
30 2 TIME, TIMESTEP,UPARAM, RHO0, RHO,
31 3 VOLUME, EINT, TEMP, JTHE,
32 4 EPSPXX, EPSPYY, EPSPZZ, EPSPXY,
33 5 EPSPYZ, EPSPZX, DEPSXX, DEPSYY,
34 6 DEPSZZ, DEPSXY, DEPSYZ, DEPSZX,
35 7 EPSXX, EPSYY, EPSZZ, EPSXY,
36 8 EPSYZ, EPSZX, SIGOXX, SIGOYY,
37 9 SIGOZZ, SIGOXY, SIGOYZ, SIGOZX,
38 A SIGNXX, SIGNYY, SIGNZZ, SIGNXY,
39 B SIGNYZ, SIGNZX, SIGVXX, SIGVYY,
40 C SIGVZZ, SIGVXY, SIGVYZ, SIGVZX,
41 D SOUNDSP, VISCMAX, UVAR, OFF,
42 E YLD, PLA, DEP, ETSE ,FHEAT )
43C-----------------------------------------------
44C I m p l i c i t T y p e s
45C-----------------------------------------------
46#include "implicit_f.inc"
47C-----------------------------------------------
48C G L O B A L P A R A M E T E R S
49C-----------------------------------------------
50#include "com01_c.inc"
51#include "param_c.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 TIME | 1 | F | R | CURRENT TIME
60C TIMESTEP| 1 | F | R | CURRENT TIME STEP
61C UVAR | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
62C RHO0 | NEL | F | R | INITIAL DENSITY
63C RHO | NEL | F | R | DENSITY
64C VOLUME | NEL | F | R | VOLUME
65C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
66C EPSPXX | NEL | F | R | STRAIN RATE XX
67C EPSPYY | NEL | F | R | STRAIN RATE YY
68C ... | | | |
69C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
70C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
71C ... | | | |
72C EPSXX | NEL | F | R | STRAIN XX
73C EPSYY | NEL | F | R | STRAIN YY
74C ... | | | |
75C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
76C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
77C ... | | | |
78C---------+---------+---+---+--------------------------------------------
79C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
80C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
81C ... | | | |
82C SIGVXX | NEL | F | W | VISCOUS STRESS XX
83C SIGVYY | NEL | F | W | VISCOUS STRESS YY
84C ... | | | |
85C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
86C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
87C---------+---------+---+---+--------------------------------------------
88C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
89C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
90C---------+---------+---+---+--------------------------------------------
91C-----------------------------------------------
92C I N P U T A r g u m e n t s
93C-----------------------------------------------
94 INTEGER, INTENT(IN) :: JTHE
95 INTEGER, INTENT(IN) :: JLAG
96 INTEGER NEL, NUPARAM, NUVAR
97 my_real
98 . TIME,TIMESTEP,UPARAM(NUPARAM),
99 . RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
100 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
101 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
102 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
103 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
104 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
105 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
106 . SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
107 . SIGOXY(NEL),SIGOYZ(NEL),SIGOZX(NEL)
108C-----------------------------------------------
109C O U T P U T A r g u m e n t s
110C-----------------------------------------------
111 my_real
112 . signxx(nel),signyy(nel),signzz(nel),
113 . signxy(nel),signyz(nel),signzx(nel),
114 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
115 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
116 . soundsp(nel),viscmax(nel),yld(nel) ,
117 . pla(nel),dep(nel),etse(nel)
118C-----------------------------------------------
119C I N P U T O U T P U T A r g u m e n t s
120C-----------------------------------------------
121 my_real :: uvar(nel,nuvar), off(nel)
122 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: fheat
123 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: temp
124C-----------------------------------------------
125C L o c a l V a r i a b l e s
126C-----------------------------------------------
127 INTEGER I,II,J,K,NITER,NINDX
128 INTEGER INDX(NEL)
129 my_real
130 . TOL,LAMDA,RES,DRES,SQR2INV,SEQ,YLDP,DAV,NU,YOUNG,SHEAR,G2,
131 . BULK,QVOCE,BVOCE,KSWIFT,KVOCE,KSWIFTP,KVOCEP,EXPV,
132 . K0,ALPHA,AN,EPS0,NN,DEPS0,EPSP0,CEPSP,CEPSPN,
133 . ETA,CP,TINI,TREF,TMELT,MTEMP,DEPSAD,SEFF,GEFF,F_EPS,GSS,
134 . GGSIG,G2GSIG,P12,P22,P33,G12,G22,G33,OMEGA,PLAP0
135 my_real
136 . gsig(6),psig(6),f_epsp(nel),f_temp(nel),svm(nel),
137 . dlam(nel),dlamin(nel),plap(nel),dheat(nel),
138 . sigtrxx(nel),sigtryy(nel),sigtrzz(nel),
139 . sigtrxy(nel),sigtryz(nel),sigtrzx(nel)
140c=========================================================================
141c State Variables
142c ---------------
143c UVAR(1) = PLA ----- equiv. plastic strain
144c UVAR(2) = YLD ----- Yield
145c UVAR(3) = PLAP --- Plastic Strain rate
146c UVAR(4) = DLAM ----- Plastic Multiplier
147C-----------------------------------------------
148c Plasticity model :
149c 1- Quadratic Yield surface
150c 2- Quadratic non associated flow rule
151c 3- Isotropic strain hardening
152c a. Swift/Voce
153c k=alpha*k1(eps,p)+(1-alpha)*k2(eps,p)
154c k1=A*(eps,p+e0)**n
155c k2=Q*(1-exp(-b*eps,p))+sig0
156c
157c b. Temperature Softening Johnson/Cook :
158c Temp_factor = 1-((T-Ttrans)/(Tmelt-Ttrans))**m
159c
160c c. Strain Rate Hardening Johnson/Cook
161c Epsp_factor = 1 + C * ln(dEps/dEps0)
162c
163c Material Properties
164c --------------------------
165c Equivalent stress
166c f(sig)=SQRT((P.sig).sig)
167c Flow potential
168c g(sig)=SQRT((G.sig).sig)
169c Hardening parameters
170c
171c Calibration of P and G
172c for planar isotropic,
173c R = (R0 + 2*R45 + R90) / 4
174c P12 = -R / (1 + R)
175c P22 = 1
176c P33 = 2 * (1 + 2*R) / (1 + R)
177c
178c for orthotropic
179c P12 = -R0 / (1 + R0)
180c P22 = (R0 / R90) * (1 + R90)/(1 + R0)
181c P33 = (1 + 2*R45) / R90 * (R0 + R90)/(1 + R0)
182c=========================================================================
183C
184C Initialize material parameters
185C --------------
186 young = uparam(1)
187 nu = uparam(2)
188 shear = young*half/(one+nu)
189 g2 = two*shear
190 lamda = young*nu/(one+nu)/(one-two*nu)
191 bulk = young*third/(one-two*nu)
192
193 p12 = uparam(3)
194 p22 = uparam(4)
195 p33 = uparam(5)
196 g12 = uparam(6)
197 g22 = uparam(7)
198 g33 = uparam(8)
199 qvoce = uparam(9)
200 bvoce = uparam(10)
201 k0 = uparam(11)
202 alpha = uparam(12)
203
204 an = uparam(13)
205 eps0 = uparam(14)
206 nn = uparam(15)
207 cepsp = uparam(16)
208 deps0 = uparam(17)
209 cepspn = cepsp
210
211 eta = uparam(18)
212 cp = uparam(19)
213 tini = uparam(20)
214 tref = uparam(21)
215 tmelt = uparam(22)
216 mtemp = uparam(23)
217 depsad = uparam(24)
218
219c Initialize state variables
220C ------------------------
221 IF (isigi == 0 .and. time == zero) THEN
222 DO i=1,nel
223 f_temp(i)= min(one,max(zero,((tini-tref)/(tmelt-tref))**mtemp))
224 f_temp(i)= one - f_temp(i)
225 f_eps = alpha*(an*eps0**nn) + (one-alpha)*k0
226 yld(i) = f_eps * f_temp(i)
227 uvar(i,2) = yld(i)
228 ENDDO
229 ENDIF
230C ------------------------
231 DO i=1,nel
232 pla(i) = uvar(i,1)
233 yld(i) = uvar(i,2)
234 plap(i) = zero
235 dlam(i) = zero
236 dep(i) = zero
237 ENDDO
238c
239C Trial stress
240C ------------
241 DO i=1,nel
242 dav = depsxx(i) + depsyy(i) + depszz(i)
243 signxx(i) = sigoxx(i) + depsxx(i)*g2 + lamda*dav
244 signyy(i) = sigoyy(i) + depsyy(i)*g2 + lamda*dav
245 signzz(i) = sigozz(i) + depszz(i)*g2 + lamda*dav
246 signxy(i) = sigoxy(i) + depsxy(i)*shear
247 signyz(i) = sigoyz(i) + depsyz(i)*shear
248 signzx(i) = sigozx(i) + depszx(i)*shear
249 sigtrxx(i) = signxx(i)
250 sigtryy(i) = signyy(i)
251 sigtrzz(i) = signzz(i)
252 sigtrxy(i) = signxy(i)
253 sigtryz(i) = signyz(i)
254 sigtrzx(i) = signzx(i)
255
256 svm(i) = sqrt(signxx(i)*signxx(i)
257 . + signyy(i)*signyy(i)*p22
258 . + signzz(i)*signzz(i)*(one+two*p12+p22)
259 . + signxy(i)*signxy(i)*p33
260 . + signyz(i)*signyz(i)*three
261 . + signzx(i)*signzx(i)*three
262 . + signxx(i)*signyy(i)*two*p12
263 . - signzz(i)*signyy(i)*two*(p12+p22)
264 . - signxx(i)*signzz(i)*two*(one+p12))
265 ENDDO
266c
267C Check Yield Condition
268C ----------------------------------------------------
269 nindx = 0
270 DO i=1,nel
271 IF (svm(i) > yld(i)) THEN ! Plastic Loading
272 nindx = nindx + 1
273 indx(nindx) = i
274 ENDIF
275 ENDDO
276c
277C Plastic increment
278C ----------------------------------------------------
279 DO ii=1,nindx
280 i = indx(ii)
281
282C----- Initialisation ----------------------------------------
283C
284 tol = 0.00001
285 sqr2inv = one/sqr2
286 niter = 5
287 cepspn = cepsp
288 IF( cepspn == zero )THEN
289 dlamin(i) = uvar(i,4) + em10
290 ELSE
291 dlamin(i) = five*em4*timestep*exp((em4)/cepspn)
292 ENDIF
293 !DLAM(I) = (SVM(I) - YLD(I)) / (SHEAR*THREE + UVAR(I,6))
294 dlam(i) = uvar(i,4) + em9
295 dlam(i) = max(dlam(i), dlamin(i))
296c
297 plap(i) = dlam(i) / timestep
298 f_temp(i) = ((temp(i)-tref)/(tmelt-tref))**mtemp
299 f_temp(i) = min(one,max(zero, f_temp(i)))
300 f_temp(i) = one - f_temp(i)
301C------------------------------------------------------------
302C------------------------------------------------------------
303 psig(1) = signxx(i) + signyy(i)*p12-signzz(i)*(one+p12)
304 psig(2) = signxx(i)*p12 + signyy(i)*p22-signzz(i)*(p12+p22)
305 psig(3) =-signxx(i)*(one+p12)- signyy(i)*(p12+p22)
306 . + signzz(i)*(one+two*p12+p22)
307 psig(4) = signxy(i)*p33
308 psig(5) = signyz(i)*three
309 psig(6) = signzx(i)*three
310
311 ggsig = dlam(i)*shear/svm(i)
312 g2gsig = ggsig*two
313 signxx(i) = sigtrxx(i) - psig(1)*g2gsig
314 signyy(i) = sigtryy(i) - psig(2)*g2gsig
315 signzz(i) = sigtrzz(i) - psig(3)*g2gsig
316 signxy(i) = sigtrxy(i) - psig(4)*ggsig
317 signyz(i) = sigtryz(i) - psig(5)*ggsig
318 signzx(i) = sigtrzx(i) - psig(6)*ggsig
319C------------------------------------------------------------
320c-------- NEWTON ITERATIONS ------
321C------------------------------------------------------------
322 DO k=1,niter
323 !HILL
324 seff = sqrt(signxx(i)*signxx(i)
325 . + signyy(i)*signyy(i)*p22
326 . + signzz(i)*signzz(i)*(one+two*p12+p22)
327 . + signxy(i)*signxy(i)*p33
328 . + signyz(i)*signyz(i)*three
329 . + signzx(i)*signzx(i)*three
330 . + signxx(i)*signyy(i)*two*p12
331 . - signzz(i)*signyy(i)*two*(p12+p22)
332 . - signxx(i)*signzz(i)*two*(one+p12))
333 seff = max(em20,seff)
334 !FLOW FUNCTION G
335 geff = sqrt(signxx(i)*signxx(i)
336 . + signyy(i)*signyy(i)*g22
337 . + signzz(i)*signzz(i)*(one+two*g12+g22)
338 . + signxy(i)*signxy(i)*g33
339 . + signyz(i)*signyz(i)*three
340 . + signzx(i)*signzx(i)*three
341 . + signxx(i)*signyy(i)*two*g12
342 . - signzz(i)*signyy(i)*two*(g12+g22)
343 . - signxx(i)*signzz(i)*two*(one+g12))
344 geff = max(em20,geff)
345c
346 gss = geff/seff
347c
348 ! P * SIGMA
349 psig(1) = signxx(i) + signyy(i)*p12 - signzz(i)*(one+p12)
350 psig(2) = signxx(i)*p12 + signyy(i)*p22 - signzz(i)*(p12+p22)
351 psig(3) =-signxx(i)*(one+p12) - signyy(i)*(p12+p22)
352 . + signzz(i)*(one+two*p12+p22)
353 psig(4) = signxy(i)*p33
354 psig(5) = signyz(i)*three
355 psig(6) = signzx(i)*three
356 ! G * SIGMA
357 gsig(1) = signxx(i) + signyy(i)*g12 - signzz(i)*(one+g12)
358 gsig(2) = signxx(i)*g12 + signyy(i)*g22 - signzz(i)*(g12+g22)
359 gsig(3) =-signxx(i)*(one+g12) - signyy(i)*(g12+g22)
360 . + signzz(i)*(one+two*g12+g22)
361 gsig(4) = signxy(i)*g33
362 gsig(5) = signyz(i)*three
363 gsig(6) = signzx(i)*three
364c
365c----- Yield, Yield prime
366 expv = exp(-bvoce*pla(i))
367 kswift = an*(pla(i) + eps0)**nn
368 kvoce = k0 + qvoce*(one - expv)
369 f_eps = alpha*kswift + (one-alpha)*kvoce
370 uvar(i,2) = f_eps*f_temp(i)
371 yld(i) = uvar(i,2)
372c-----
373 f_epsp(i) = one + cepspn*log(plap(i)/deps0)
374 kswiftp = kswift*nn / (pla(i) + eps0)
375 kvocep = qvoce*bvoce*expv
376 yldp = alpha*kswiftp + (one - alpha)*kvocep
377 yldp = yldp * f_temp(i) * f_epsp(i)
378 yldp = yldp + yld(i) * cepspn / dlam(i)
379 yld(i) = uvar(i,2) * f_epsp(i)
380c f(sig)-yield
381 res = seff - yld(i)
382c
383 ! Check for convergence
384 IF ( (abs(res) < tol*uvar(i,2)) .and.
385 + (yld(i) > uvar(i,2)*(one-tol)) .and.
386 + ( plap(i) > deps0 ) .and.
387 + (abs(plap(i)-plap0) < plap0 *tol)
388 + ) EXIT
389c
390 ! Update d_lambda & plastic strain
391c
392 dres = psig(1)*gsig(1) + psig(2)*gsig(2) + psig(3)*gsig(3)
393 . + half*(psig(4)*gsig(4)+psig(5)*gsig(5)+psig(6)*gsig(6))
394 dres = dres / (seff*geff)
395 dres = -g2 * dres - yldp * gss
396 dlam(i) = dlam(i) - res / dres
397 uvar(i,4) = dlam(i)
398 ! Check for convergence and Strain Rate > Reference Strain Rate
399 IF ( (abs(res) < tol*uvar(i,2)) .and.
400 + (gss*dlam(i) / timestep < deps0) ) THEN
401 cepspn = zero
402 dlam(i) = uvar(i,4) + em10
403 dlamin(i) = em20
404 END IF
405c
406 plap0 = plap(i)
407 dlam(i) = max(dlam(i), dlamin(i))
408 dep(i) = dlam(i) * gss
409 pla(i) = uvar(i,1) + dep(i)
410 plap(i) = dep(i) / timestep
411 uvar(i,4) = dlam(i)
412C
413 ggsig = dlam(i)*shear/geff
414 g2gsig = ggsig*two
415 signxx(i) = sigtrxx(i) - gsig(1)*g2gsig
416 signyy(i) = sigtryy(i) - gsig(2)*g2gsig
417 signzz(i) = sigtrzz(i) - gsig(3)*g2gsig
418 signxy(i) = sigtrxy(i) - gsig(4)*ggsig
419 signyz(i) = sigtryz(i) - gsig(5)*ggsig
420 signzx(i) = sigtrzx(i) - gsig(6)*ggsig
421
422c
423 END DO ! K=1,NITER
424c----------- end Newton iterations
425c
426 ENDDO ! plastic increment loop
427c
428!-------------------------------------------------------------
429! Calculate Plastic Work / Temperature depending on omega
430!-------------------------------------------------------------
431 DO ii=1,nindx
432 i = indx(ii)
433 IF (plap(i) <= deps0)THEN
434 omega = zero
435 ELSEIF (plap(i) > depsad) THEN
436 omega = one
437 ELSE
438 omega = ((plap(i)-deps0)**2 )*
439 . (three*depsad - two*plap(i) - deps0)
440 . /((depsad-deps0)**3)
441 ENDIF
442 dheat(i) = eta*omega*yld(i)*dep(i)
443 ENDDO
444!
445 IF (jthe < 0 .AND. jlag /= 0) THEN
446 DO ii=1,nindx
447 i = indx(ii)
448 fheat(i) = fheat(i) + dheat(i)*volume(i)
449 ENDDO
450 ELSE
451 DO ii=1,nindx
452 i = indx(ii)
453 temp(i) = temp(i) + dheat(i) / (rho0(i)*cp)
454 ENDDO
455 ENDIF
456C-----------
457 DO i=1,nel
458 uvar(i,1) = pla(i)
459 uvar(i,2) = yld(i)
460 uvar(i,3) = plap(i)
461 uvar(i,4) = dlam(i)
462 etse(i) = one
463 soundsp(i)= sqrt((bulk+four_over_3*shear)/rho0(i))
464 viscmax(i)= zero
465 ENDDO
466c-----------
467 RETURN
468 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps84
subroutine sigeps84(nel, nuparam, nuvar, jlag, time, timestep, uparam, rho0, rho, volume, eint, temp, jthe, 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, yld, pla, dep, etse, fheat)
Definition sigeps84.F:43