OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps83.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!|| sigeps83 ../engine/source/materials/mat/mat083/sigeps83.F
25!||--- called by ------------------------------------------------------
26!|| suser43 ../engine/source/elements/solid/sconnect/suser43.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||====================================================================
30 SUBROUTINE sigeps83(
31 1 NEL ,TIME ,TIMESTEP,UPARAM ,OFF ,
32 2 EPSD ,STIFM ,IFUNC ,MAXFUNC ,NPF ,TF ,
33 3 AREA ,DEPSZZ ,DEPSYZ ,DEPSZX ,NUPARAM ,EPSZZ ,
34 4 SIGOZZ ,SIGOYZ ,SIGOZX ,SIGNZZ ,SIGNYZ ,SIGNZX ,
35 5 PLA ,JSMS ,DMELS ,SYM, UVAR ,NUVAR ,
36 6 DMG ,ASRATE )
37C-----------------------------------------------
38C I m p l i c i t T y p e s
39C-----------------------------------------------
40#include "implicit_f.inc"
41#include "com01_c.inc"
42C---------+---------+---+---+--------------------------------------------
43C VAR | SIZE |TYP| RW| DEFINITION
44C---------+---------+---+---+--------------------------------------------
45C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
46C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
47C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
48C JSMS | 1 | I | R | 0/1 (=1 IF /DT/AMS APPLIES TO THIS ELEMENT GROUP)
49C---------+---------+---+---+--------------------------------------------
50C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
51C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
52C NPF | * | I | R | FUNCTION ARRAY
53C TF | * | F | R | FUNCTION ARRAY
54C---------+---------+---+---+--------------------------------------------
55C TIME | 1 | F | R | CURRENT TIME
56C TIMESTEP| 1 | F | R | CURRENT TIME STEP
57C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
58C ... | | | |
59C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
60C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
61C ... | | | |
62C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
63C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
64C ... | | | |
65C---------+---------+---+---+--------------------------------------------
66C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
67C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
68C ... | | | |
69C DMELS | NEL | F | W | NON DIAGONAL TERM FOR AMS
70C---------+---------+---+---+--------------------------------------------
71C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
72C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
73C---------+---------+---+---+--------------------------------------------
74#include "sms_c.inc"
75C----------------------------------------------------------
76C D u m m y A r g u m e n t s
77C----------------------------------------------------------
78 INTEGER NEL,NUPARAM,NUVAR,MAXFUNC,JSMS
79 INTEGER IFUNC(*),NPF(*)
80 my_real
81 . TIME,TIMESTEP,ASRATE
82 my_real
83 . UPARAM(NUPARAM),OFF(NEL),TF(*),PLA(NEL),AREA(NEL),
84 . epsd(nel),depszz(nel),depsyz(nel),depszx(nel),
85 . sigozz(nel),sigoyz(nel),sigozx(nel),epszz(nel),
86 . stifm(nel),signzz(nel),signyz(nel),signzx(nel),
87 . dmels(*),uvar(nel,nuvar), sym(nel),dmg(nel)
88C-----------------------------------------------
89C L o c a l V a r i a b l e s
90C-----------------------------------------------
91 INTEGER :: II,IEL,ITER,NITER,IRATE,NRATE,IDYIELD,IPLAS,IFUNN,IFUNT,
92 . ICOMP,NTRAC,NCOMP,VP
93 INTEGER ,DIMENSION(NEL) :: ELTRAC,ELCOMP
94 my_real :: YOUNGT,YOUNGC,SHEAR,DYDX,
95 . SVMN,SVMT,DTB,ALPHA,BETA,SVM,TS,TN,XSCALE,RNC,RSC,XFAC,YFAC,
96 . AA,AN,AS,NORMEF,NZZ,NYZ,NZX,FACN,FACT,SZZ,SYZ,SZX,
97 . fpx,fpy,fprim,sign,sigt,sig_eff,phi,dtinv,dszz,dsyz,dszx,dyld,yld0
98 my_real ,DIMENSION(NEL) :: dsigx,dsigy,dsigz,dep,depst,stf,epsp,ht,
99 . fyield,hyield, rn,hn,rs,sigtrzz,sigtryz,sigtrzx,young
100C----------------------------------------------------------
101C E x t e r n a l F u n c t i o n
102C----------------------------------------------------------
103 my_real :: finter
104c-----------------------------------------------------
105c UVAR(1) = strain rate
106C=======================================================================
107C INPUT PARAMETERS INITIALIZATION
108C-----------------------------------------------
109 ifunn = ifunc(1)
110 ifunt = ifunc(2)
111 idyield = ifunc(3)
112 youngt = uparam(1)
113 alpha = uparam(2)
114 beta = uparam(3)
115 yfac = uparam(4)
116 xscale = uparam(5)
117 rnc = uparam(6)
118 rsc = uparam(7)
119 xfac = uparam(8)
120 iplas = uparam(10)
121 shear = uparam(11)
122 icomp = nint(uparam(12))
123 youngc = uparam(13)
124 vp = nint(uparam(14))
125 dtinv = one / max(em20,timestep)
126 niter = 3
127c
128 IF (vp == 0) THEN ! total displacement rate
129 epsp(1:nel) = sqrt(depszz(1:nel)**2 + depsyz(1:nel)**2 + depszx(1:nel)**2) *dtinv
130 uvar(1:nel,1) = asrate*epsp(1:nel) + (one - asrate)*uvar(:nel,1)
131 epsd(1:nel) = uvar(1:nel,1)
132 ENDIF
133 epsp(:nel) = uvar(:nel,1)
134 dep(:nel) = zero
135c----------------------------------
136 ntrac = 0
137 ncomp = 0
138 IF (icomp == 0) THEN ! elasto-plastic in tension & compression
139 IF (youngc == youngt) THEN ! common Young modulus
140 young(:nel) = youngt
141 ELSE
142 DO iel=1,nel
143 IF (sigozz(iel) > zero) THEN ! element in tension
144 young(iel) = youngt
145 ELSE ! element in compression
146 young(iel) = youngc
147 ENDIF
148 ENDDO
149 END IF
150 ELSE ! ICOMP = 1 (linear in compression, elasto-plastic in tension)
151 eltrac(1:nel) = 0
152 elcomp(1:nel) = 0
153 DO iel=1,nel
154 IF (sigozz(iel) > zero) THEN ! element in tension
155 young(iel) = youngt
156 ntrac = ntrac + 1
157 eltrac(ntrac) = iel
158 ELSE ! element in compression
159 young(iel) = youngc
160 ncomp = ncomp + 1
161 elcomp(ncomp) = iel
162 ENDIF
163 ENDDO
164 ENDIF
165c
166 stf(1:nel) = young(1:nel) * area(1:nel)
167 dsigz(1:nel) = young(1:nel) * depszz(1:nel)
168 dsigy(1:nel) = shear*depsyz(1:nel)
169 dsigx(1:nel) = shear*depszx(1:nel)
170 signzz(1:nel) = sigozz(1:nel) + dsigz(1:nel) * off(1:nel)
171 signyz(1:nel) = sigoyz(1:nel) + dsigy(1:nel) * off(1:nel)
172 signzx(1:nel) = sigozx(1:nel) + dsigx(1:nel) * off(1:nel)
173 stifm(1:nel) = stifm(1:nel) + stf(1:nel)*off(1:nel)
174 sigtrzz(1:nel) = signzz(1:nel)
175 sigtryz(1:nel) = signyz(1:nel)
176 sigtrzx(1:nel) = signzx(1:nel)
177c-----------------------------------------------------
178c omega = sqrt(2k/2*dmels), dt=2/omega, 2*dmels=dt**2 * 2k / 4
179 IF (idtmins==2 .AND. jsms/=0) THEN
180 dtb = (dtmins/dtfacs)**2
181 DO iel=1,nel
182 dmels(iel)=max(dmels(iel),half*dtb*stf(iel)*off(iel))
183 ENDDO
184 END IF
185c
186c----------------------Interpolation
187c
188 IF (ifunn /= 0) THEN
189 DO iel=1,nel
190 rn(iel) = finter(ifunn,epsp(iel)*xscale,npf,tf,dydx) * rnc !/ area(iel)
191 ENDDO
192 ELSE
193 DO iel=1,nel
194 rn(iel) = rnc !/ AREA(IEL)
195 ENDDO
196 ENDIF
197
198 IF (ifunt /=0)THEN
199 DO iel=1,nel
200 rs(iel) = finter(ifunt,epsp(iel)*xscale,npf,tf,dydx) * rsc !/ AREA(IEL)
201 ENDDO
202 ELSE
203 DO iel=1,nel
204 rs(iel) = rsc ! / AREA(IEL)
205 ENDDO
206 ENDIF
207 DO iel=1,nel
208 fyield(iel) = max(zero, finter(idyield,pla(iel)*xfac,npf,tf,dydx))
209 fyield(iel) = fyield(iel)*(one-dmg(iel))*yfac
210 hyield(iel) = dydx*yfac
211 ENDDO
212c-----------------------------------------------------------------------------
213c Plasticity projection (radial return from sigma trial)
214c-----------------------------------------------------------------------------
215 IF (beta == two) THEN
216 IF (icomp == 0) THEN ! elasto-plastic in tension and compression
217 DO iel = 1,nel
218 aa = rn(iel)*(one-alpha*sin(sym(iel)))
219 an = one/max(em20,aa)**2
220 as = one/max(em20,rs(iel))**2
221 sig_eff = an * sigtrzz(iel)** 2 + as * (sigtryz(iel)**2 + sigtrzx(iel)**2)
222 yld0 = fyield(iel)
223 phi = sig_eff - yld0**2
224 IF (sig_eff > zero .and. phi > zero) THEN ! plastic
225 normef = sqrt(sigtrzz(iel)**2 + sigtryz(iel)**2 + sigtrzx(iel)**2)
226 nzz = sigtrzz(iel) / normef
227 nyz = sigtryz(iel) / normef
228 nzx = sigtrzx(iel) / normef
229 facn = two*an*young(iel)
230 fact = two*as*shear
231 DO iter= 1,niter
232
233 dszz = -facn*signzz(iel)*nzz
234 dsyz = -fact*signyz(iel)*nyz
235 dszx = -fact*signzx(iel)*nzx
236 dyld = -two*fyield(iel)*hyield(iel)
237 fprim = dszz + dsyz + dszx + dyld
238 IF(fprim .NE. zero) THEN
239 dep(iel) = dep(iel) - phi / fprim
240 signzz(iel) = sigtrzz(iel) - young(iel)*dep(iel) * nzz
241 signyz(iel) = sigtryz(iel) - shear*dep(iel) * nyz
242 signzx(iel) = sigtrzx(iel) - shear*dep(iel) * nzx
243 fyield(iel) = yld0 + hyield(iel)*dep(iel)
244 sig_eff = an * signzz(iel)** 2 + as * (signyz(iel)**2 + signzx(iel)**2)
245 phi = sig_eff - fyield(iel)**2 ! yield criterion function
246 ENDIF
247 ENDDO
248 dep(iel) = max(zero, dep(iel))*off(iel)
249 pla(iel) = pla(iel) + dep(iel)
250 ENDIF
251 ENDDO
252c--------------
253 ELSE ! ICOMP = 1 (linear in compression, elasto-plastic in tension)
254c--------------
255 DO ii = 1,ntrac
256 iel= eltrac(ii)
257 aa = rn(iel)*(one-alpha*sin(sym(iel)))
258 an = one/max(em20,aa)**2
259 as = one/max(em20,rs(iel))**2
260 sig_eff = an * sigtrzz(iel)** 2 + as * (sigtryz(iel)**2 + sigtrzx(iel)**2)
261 yld0 = fyield(iel)
262 phi = sig_eff - yld0**2
263 IF (sig_eff > zero .and. phi > zero) THEN ! plastic
264 normef = sqrt(sigtrzz(iel)**2 + sigtryz(iel)**2 + sigtrzx(iel)**2)
265 nzz = sigtrzz(iel) / normef
266 nyz = sigtryz(iel) / normef
267 nzx = sigtrzx(iel) / normef
268 facn = two*an*young(iel)
269 fact = two*as*shear
270 DO iter= 1,niter
271 dszz = -facn*signzz(iel)*nzz
272 dsyz = -fact*signyz(iel)*nyz
273 dszx = -fact*signzx(iel)*nzx
274 dyld = -two*fyield(iel)*hyield(iel)
275 fprim = dszz + dsyz + dszx + dyld
276
277 IF(fprim .NE. zero) THEN
278 dep(iel) = dep(iel) - phi / fprim
279 signzz(iel) = sigtrzz(iel) - young(iel)*dep(iel) * nzz
280 signyz(iel) = sigtryz(iel) - shear*dep(iel) * nyz
281 signzx(iel) = sigtrzx(iel) - shear*dep(iel) * nzx
282 fyield(iel) = yld0 + hyield(iel)*dep(iel)
283 sig_eff = an * signzz(iel)** 2 + as * (signyz(iel)**2 + signzx(iel)**2)
284 phi = sig_eff - fyield(iel)**2 ! yield criterion function
285 ENDIF
286 ENDDO
287 dep(iel) = max(zero, dep(iel))*off(iel)
288 pla(iel) = pla(iel) + dep(iel)
289 ENDIF
290 ENDDO
291c
292 DO ii = 1,ncomp
293 iel= elcomp(ii)
294 as = one/max(em20,rs(iel))**2
295 sig_eff = as * (sigtryz(iel)**2 + sigtrzx(iel)**2)
296 yld0 = fyield(iel)
297 phi = sig_eff - yld0**2
298 IF (sig_eff > zero .and. phi > zero) THEN ! plastic
299 normef = sqrt(sigtryz(iel)**2 + sigtrzx(iel)**2)
300 nyz = sigtryz(iel) / normef
301 nzx = sigtrzx(iel) / normef
302 fact = two*as*shear
303 DO iter= 1,niter
304 dsyz = -fact*signyz(iel)*nyz
305 dszx = -fact*signzx(iel)*nzx
306 dyld = -two*fyield(iel)*hyield(iel)
307 fprim = dsyz + dszx + dyld
308 IF(fprim .NE. zero) THEN
309 dep(iel) = dep(iel) - phi / fprim
310 signyz(iel) = sigtryz(iel) - shear*dep(iel) * nyz
311 signzx(iel) = sigtrzx(iel) - shear*dep(iel) * nzx
312 fyield(iel) = yld0 + hyield(iel)*dep(iel)
313 sig_eff = as * (signyz(iel)**2 + signzx(iel)**2)
314 phi = sig_eff - fyield(iel)**2 ! yield criterion function
315 ENDIF
316 ENDDO
317 dep(iel) = max(zero, dep(iel))*off(iel)
318 pla(iel) = pla(iel) + dep(iel)
319 ENDIF
320 ENDDO
321c
322 END IF
323c------------------------
324 ELSE ! BETA /= 2
325c------------------------
326 IF (icomp == 0) THEN ! elasto-plastic in tension and compression
327 DO iel=1,nel
328 aa = rn(iel)*(one-alpha*sin(sym(iel)))
329 an = one/max(em20,aa)**beta
330 as = one/max(em20,rs(iel))**beta
331 svmn = abs(signzz(iel)) ! plast - normal stress
332 svmt = sqrt(signyz(iel)**2 + signzx(iel)**2) ! plast - shear stress
333 sig_eff = an * svmn**beta + as * svmt**beta
334 yld0 = fyield(iel)
335 phi = sig_eff - yld0**beta
336 IF (sig_eff > zero .and. phi > zero) THEN ! plastic
337 normef = sqrt(sigtrzz(iel)**2+sigtryz(iel)**2 + sigtrzx(iel)**2)
338 nzz = sigtrzz(iel) / normef
339 nyz = sigtryz(iel) / normef
340 nzx = sigtrzx(iel) / normef
341 facn = beta*an*young(iel)
342 fact = beta*as*shear
343
344 DO iter=1,niter
345 svmt = max(svmt,em20)
346 szz = signzz(iel)
347 syz = signyz(iel)
348 szx = signzx(iel)
349 dszz = -facn*nzz*abs(szz)**(beta-one) * sign(one,szz)
350 dsyz = -fact*nyz*syz*svmt**(beta-two)
351 dszx = -fact*nzx*szx*svmt**(beta-two)
352 dyld = -beta*hyield(iel)*fyield(iel)**(beta-one)
353 fprim = dszz + dsyz + dszx + dyld
354 IF(fprim .NE. zero) THEN
355 dep(iel) = dep(iel) - phi / fprim
356 signzz(iel) = sigtrzz(iel) - young(iel)*dep(iel) * nzz
357 signyz(iel) = sigtryz(iel) - shear*dep(iel) * nyz
358 signzx(iel) = sigtrzx(iel) - shear*dep(iel) * nzx
359 fyield(iel) = max(zero, yld0 + hyield(iel)*dep(iel))
360 svmt = sqrt(signyz(iel)**2 + signzx(iel)**2)
361 sig_eff = an * abs(signzz(iel))**beta + as * svmt**beta
362 phi = sig_eff - fyield(iel)**beta
363 ENDIF
364 ENDDO
365 dep(iel) = max(zero, dep(iel))*off(iel)
366 pla(iel) = pla(iel) + dep(iel)
367 ENDIF
368 ENDDO
369c
370 ELSE ! ICOMP = 1 ! elastic in compression
371c
372 DO ii = 1,ntrac
373 iel= eltrac(ii)
374 aa = rn(iel)*(one-alpha*sin(sym(iel)))
375 an = one/max(em20,aa)**beta
376 as = one/max(em20,rs(iel))**beta
377 svmn = abs(signzz(iel)) ! plast - normal stress
378 svmt = sqrt(signyz(iel)**2 + signzx(iel)**2) ! plast - shear stress
379 sig_eff = an * svmn**beta + as * svmt**beta
380 yld0 = fyield(iel)
381 phi = sig_eff - yld0**beta
382 IF (sig_eff > zero .and. phi > zero) THEN ! plastic
383 normef = sqrt(sigtrzz(iel)**2+sigtryz(iel)**2 + sigtrzx(iel)**2)
384 nzz = sigtrzz(iel) / normef
385 nyz = sigtryz(iel) / normef
386 nzx = sigtrzx(iel) / normef
387 facn = beta*an*young(iel)
388 fact = beta*as*shear
389
390 DO iter=1,niter
391 svmt = max(svmt,em20)
392 szz = signzz(iel)
393 syz = signyz(iel)
394 szx = signzx(iel)
395 dszz = -facn*nzz*abs(szz)**(beta-one) * sign(one,szz)
396 dsyz = -fact*nyz*syz*svmt**(beta-two)
397 dszx = -fact*nzx*szx*svmt**(beta-two)
398 dyld = -beta*hyield(iel)*fyield(iel)**(beta-one)
399 fprim = dszz + dsyz + dszx + dyld
400 IF(fprim .NE. zero) THEN
401 dep(iel) = dep(iel) - phi / fprim
402 signzz(iel) = sigtrzz(iel) - young(iel)*dep(iel) * nzz
403 signyz(iel) = sigtryz(iel) - shear*dep(iel) * nyz
404 signzx(iel) = sigtrzx(iel) - shear*dep(iel) * nzx
405 fyield(iel) = max(zero, yld0 + hyield(iel)*dep(iel))
406 svmt = sqrt(signyz(iel)**2 + signzx(iel)**2)
407 sig_eff = an * abs(signzz(iel))**beta + as * svmt**beta
408 phi = sig_eff - fyield(iel)**beta
409 ENDIF
410 ENDDO
411 dep(iel) = max(zero, dep(iel))*off(iel)
412 pla(iel) = pla(iel) + dep(iel)
413 ENDIF
414 ENDDO
415c
416 DO ii = 1,ncomp
417 iel= elcomp(ii)
418 as = one/max(em20,rs(iel))**beta
419 svmt = sqrt(signyz(iel)**2 + signzx(iel)**2) ! plast - shear stress
420 sig_eff = as * svmt**beta
421 yld0 = fyield(iel)
422 phi = sig_eff - yld0**beta
423 IF (sig_eff > zero .and. phi > zero) THEN ! plastic
424 normef = svmt
425 nyz = sigtryz(iel) / normef
426 nzx = sigtrzx(iel) / normef
427 fact = beta*as*shear
428
429 DO iter=1,niter
430 svmt = max(svmt,em20)
431 syz = signyz(iel)
432 szx = signzx(iel)
433 dsyz = -fact*nyz*syz*svmt**(beta-two)
434 dszx = -fact*nzx*szx*svmt**(beta-two)
435 dyld = -beta*hyield(iel)*fyield(iel)**(beta-one)
436 fprim = dsyz + dszx + dyld
437 IF(fprim .NE. zero) THEN
438 dep(iel) = dep(iel) - phi / fprim
439 signyz(iel) = sigtryz(iel) - shear*dep(iel) * nyz
440 signzx(iel) = sigtrzx(iel) - shear*dep(iel) * nzx
441 fyield(iel) = yld0 + hyield(iel)*dep(iel)
442 svmt = sqrt(signyz(iel)**2 + signzx(iel)**2)
443 sig_eff = as * svmt**beta
444 phi = sig_eff - fyield(iel)**beta
445 ENDIF
446 ENDDO
447 dep(iel) = max(zero, dep(iel))*off(iel)
448 pla(iel) = pla(iel) + dep(iel)
449 ENDIF
450 ENDDO
451c
452 END IF ! ICOMP
453c
454 ENDIF ! BETA
455c----------------------------end iplas 2--------------------------------------
456c normal projection option (hidden: iplas = 1) deleted in august 2021
457c------------------
458 IF (vp == 1) THEN
459 DO iel=1,nel
460 epsp(iel) = dep(iel) * dtinv ! plastic strain rate
461 uvar(iel,1) = asrate*epsp(iel) + (one - asrate)*uvar(iel,1)
462 epsd(iel) = uvar(iel,1)
463 ENDDO
464 ENDIF
465C-----------
466 RETURN
467 END
max
#define max(a, b)
Definition macros.h:21
sigeps83
subroutine sigeps83(nel, time, timestep, uparam, off, epsd, stifm, ifunc, maxfunc, npf, tf, area, depszz, depsyz, depszx, nuparam, epszz, sigozz, sigoyz, sigozx, signzz, signyz, signzx, pla, jsms, dmels, sym, uvar, nuvar, dmg, asrate)
Definition sigeps83.F:37