OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps36pi.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!|| sigeps36pi ../engine/source/materials/mat/mat036/sigeps36pi.F
25!||--- called by ------------------------------------------------------
26!|| mulaw_ib ../engine/source/elements/beam/mulaw_ib.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| vinter ../engine/source/tools/curve/vinter.F
30!||====================================================================
31 SUBROUTINE sigeps36pi(
32 1 NEL ,KFUNC ,IPT ,NPF ,TF ,
33 2 NGL ,TIME ,TIMESTEP,UPARAM ,IPM ,
34 3 IMAT ,E ,G ,OFF ,DEPSXX ,
35 4 DEPSXY ,DEPSXZ ,SIGOXX ,SIGOXY ,SIGOXZ ,
36 5 SIGNXX ,SIGNXY ,SIGNXZ ,ETSE ,PLA ,
37 6 EPSP ,NVARTMP ,VARTMP ,SIGY)
38C---------+---------+---+---+--------------------------------------------
39C VAR | SIZE |TYP| RW| DEFINITION
40C---------+---------+---+---+--------------------------------------------
41C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
42C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
43C---------+---------+---+---+--------------------------------------------
44C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
45C KFUNC | NFUNC | I | R | FUNCTION INDEX BIDON
46C NPF | * | I | R | FUNCTION ARRAY
47C IPT | 1 | I | R | INTEGRATION POINT NUMBER
48C TF | * | F | R | FUNCTION ARRAY
49C---------+---------+---+---+--------------------------------------------
50C TIME | 1 | F | R | CURRENT TIME
51C TIMESTEP| 1 | F | R | CURRENT TIME STEP
52C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
53C EPSP | NEL | F | R | STRAIN RATE
54C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
55C ... | | | |
56C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
57C ... | | | |
58C---------+---------+---+---+--------------------------------------------
59C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
60C ... | | | |
61C---------+---------+---+---+--------------------------------------------
62C PLA | NEL | F |R/W| PLASTIC STRAIN
63C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
64C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
65C-----------------------------------------------
66C I m p l i c i t T y p e s
67C-----------------------------------------------
68#include "implicit_f.inc"
69C-----------------------------------------------
70C G l o b a l P a r a m e t e r s
71C-----------------------------------------------
72#include "comlock.inc"
73C---------+---------+---+---+--------------------------------------------
74C C o m m o n B l o c k s
75C-----------------------------------------------
76#include "param_c.inc"
77#include "units_c.inc"
78C-----------------------------------------------
79C I N P U T A r g u m e n t s
80C-----------------------------------------------
81 INTEGER ,INTENT(IN) :: IMAT
82 INTEGER NEL,NUPARAM,IPT,NVARTMP
83 INTEGER NGL(NEL),NPF(*),IPM(NPROPMI,*),KFUNC(*)
84 my_real
85 . TIME,TIMESTEP,UPARAM(*),PLA(NEL),EPSP(NEL),
86 . DEPSXX(NEL),DEPSXY(NEL),DEPSXZ(NEL),
87 . sigoxx(nel),sigoxy(nel),sigoxz(nel),
88 . e(*),g(*),tf(*)
89C-----------------------------------------------
90C O U T P U T A r g u m e n t s
91C-----------------------------------------------
92 my_real
93 . signxx(nel),signxy(nel),signxz(nel),etse(nel),
94 . off(nel)
95 INTEGER, INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
96 my_real, INTENT(INOUT) :: SIGY(NEL)
97C-----------------------------------------------
98C E X T E R N A L F u n c t i o n s
99C-----------------------------------------------
100 my_real
101 . finter
102 EXTERNAL finter
103C-----------------------------------------------
104C L o c a l V a r i a b l e s
105C-----------------------------------------------
106 INTEGER I,J,J1,J2,N,NINDX,IADBUF,NFUNC,IPLAS,IFAIL,MFUNCE,ISMOOTH
107 INTEGER NRATE(NEL),IAD1(NEL),IPOS1(NEL),ILEN1(NEL),
108 . iad2(nel),ipos2(nel),ilen2(nel),jj(nel),
109 . ifunc(nel,100),opte(nel)
110 my_real :: g3(nel),dydx1(nel),dydx2(nel),rate(nel,2),
111 . y1(nel),y2(nel),dr(nel),yfac(nel,2),
112 . pp(nel),qq(nel),fail(nel),svm(nel),h(nel),
113 . epsmax(nel),epsr1(nel),epsr2(nel),
114 . sigexx(nel),sigeyy(nel),sigexy(nel),
115 . hk(nel),nu(nel),nu1(nel),epsf(nel),yld(nel),
116 . aa(nel),bb(nel),cc(nel),einf(nel),ce(nel),
117 . escale(nel),dydxe(nel)
118 my_real :: r,umr,fac,gs,svm1,shfact,epsp1,epsp2,rfac
119C=======================================================================
120 iplas = 1
121 shfact = five_over_6
122 DO i=1,nel
123 nfunc = ipm(10,imat)
124 DO j=1,nfunc
125 ifunc(i,j)=ipm(10+j,imat)
126 ENDDO
127 fail(i) = one
128 ENDDO
129C
130 iadbuf = ipm(7,imat)-1
131 ifail = 0
132 DO i=1,nel
133 e(i) = uparam(iadbuf+2)
134 g(i) = uparam(iadbuf+5)
135 g3(i) = three*g(i)
136 nu(i) = uparam(iadbuf+6)
137 nrate(i) = nint(uparam(iadbuf+1))
138 epsmax(i)= uparam(iadbuf+6+2*nrate(i)+1)
139c---------------------
140 opte(i) = uparam(iadbuf+2*nrate(i) + 23)
141 einf(i) = uparam(iadbuf+2*nrate(i) + 24)
142 ce(i) = uparam(iadbuf+2*nrate(i) + 25)
143 ifail = nint(uparam(iadbuf+2*nrate(i) + 27)) ! flag for failure
144 ismooth = nint(uparam(iadbuf+2*nrate(i) + 29)) ! strain rate interpolation flag
145c--------------------
146 IF (opte(i) == 1)THEN
147 IF(pla(i) > zero)THEN
148 mfunce= uparam(iadbuf+2*nrate(i)+ 22)
149 escale(i) = finter(ifunc(i,mfunce),pla(i),npf,tf,dydxe(i))
150 e(i) = escale(i)* e(i)
151 g(i) = half*e(i)/(one+nu(i))
152 g3(i) = three*g(i)
153 ENDIF
154 ELSEIF (ce(i) /= zero) THEN
155 IF(pla(i) > zero)THEN
156 e(i) = e(i)-(e(i)-einf(i))*(one-exp(-ce(i)*pla(i)))
157 g(i) = half*e(i)/(one+nu(i))
158 g3(i) = three*g(i)
159 ENDIF
160 ENDIF
161 ENDDO
162C
163C-------------------
164c DO I=1,NEL
165C--- Module Cisaillement Transverse
166c SH = FIVE_OVER_6*G(I)
167c YMA = TWELVE*E(I)/AL(I)**2
168c SH1 = YMA*IYY(I)
169c SH2 = YMA*IZZ(I)
170C shf = ishear ( no shear = 1, shear = 0 = default)
171c SH0 = (ONE - SHF(I))*SH
172c GS1 = SH0*SH1/(SH+SH1) + SHF(I)*SH1
173c GS2 = SH0*SH2/(SH+SH2) + SHF(I)*SH2
174C--- Contraintes Elastiques
175C--- CONTRAINTES PLASTIQUEMENT ADMISSIBLES
176 DO i=1,nel
177 gs = shfact*g(i)
178 signxx(i) = sigoxx(i) + e(i)*depsxx(i)
179 signxy(i) = sigoxy(i) + gs*depsxy(i)
180 signxz(i) = sigoxz(i) + gs*depsxz(i)
181 etse(i) = one
182 ENDDO
183C-------------------
184C Yield Stress
185C-------------------
186 jj(1:nel) = 1
187 DO i=1,nel
188 iadbuf = ipm(7,imat)-1
189 DO j=2,nrate(i)-1
190 IF (epsp(i) > uparam(iadbuf+6+j)) jj(i) = j
191 ENDDO
192 ENDDO
193 DO i=1,nel
194 iadbuf = ipm(7,imat)-1
195 rate(i,1) = uparam(iadbuf+6+jj(i))
196 yfac(i,1) = uparam(iadbuf+6+nrate(i)+jj(i))
197 IF (nrate(i) > 1) THEN
198 rate(i,2) = uparam(iadbuf+6+jj(i)+1)
199 yfac(i,2) = uparam(iadbuf+6+nrate(i)+jj(i)+1)
200 ENDIF
201 ENDDO
202C
203 DO i=1,nel
204 j1 = jj(i)
205 j2 = j1+1
206 ipos2(i) = 1
207 iad2(i) = 1
208 ilen2(i) = 1
209 ipos1(i) = vartmp(i,j1)
210 iad1(i) = npf(ifunc(i,j1)) / 2 + 1
211 ilen1(i) = npf(ifunc(i,j1)+1) / 2 - iad1(i) - ipos1(i)
212 IF (nrate(i) > 1) THEN
213 ipos2(i) = vartmp(i,j2)
214 iad2(i) = npf(ifunc(i,j2)) / 2 + 1
215 ilen2(i) = npf(ifunc(i,j2)+1) / 2 - iad2(i) - ipos2(i)
216 ENDIF
217 ENDDO
218C
219 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
220 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
221C
222 DO i=1,nel
223 IF (nrate(i) == 1) THEN
224 j1 = jj(i)
225 yld(i) = fail(i)*y1(i) * yfac(i,1)
226 yld(i) = max(yld(i),em20)
227 h(i) = fail(i)*dydx1(i)*yfac(i,1)
228 vartmp(i,j1) = ipos1(i)
229 ELSE
230 j1 = jj(i)
231 j2 = j1+1
232 y1(i) = y1(i)*yfac(i,1)
233 y2(i) = y2(i)*yfac(i,2)
234 IF (ismooth == 2) THEN ! log_n based interpolation of strain rate
235 epsp1 = max(rate(i,1), em20)
236 epsp2 = rate(i,2)
237 fac = log(max(epsp(i),em20)/epsp1) / log(epsp2/epsp1)
238 ELSE ! linear interpolation of strain rate
239 epsp1 = rate(i,1)
240 epsp2 = rate(i,2)
241 fac = (epsp(i) - epsp1) / (epsp2 - epsp1)
242 END IF
243 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
244 yld(i) = max(yld(i),em20)
245 dydx1(i)= dydx1(i)*yfac(i,1)
246 dydx2(i)= dydx2(i)*yfac(i,2)
247 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
248 vartmp(i,j1) = ipos1(i)
249 vartmp(i,j2) = ipos2(i)
250 END IF
251 sigy(i) = yld(i)
252 ENDDO
253C-------------------
254C PROJECTION
255C-------------------
256c IF (IPLAS == 1) THEN
257C--- radial return
258 DO i=1,nel
259 svm1 = signxx(i)**2 + three*(signxy(i)**2 + signxz(i)**2)
260 svm1 = sqrt(svm1)
261 r = min( one, yld(i)/max(em20,svm1) )
262 signxx(i)=signxx(i)*r
263 signxy(i)=signxy(i)*r
264 signxz(i)=signxz(i)*r
265 umr = one - r
266 pla(i) = pla(i) + off(i)*svm1*umr/(g3(i)+h(i))
267 IF (r < one) etse(i)= h(i)/(h(i)+e(i))
268 ENDDO
269c ELSE
270C--- plasticite uniaxiale : 3 iterations Newton
271ctmp DO J=1,NEL
272ctmp AA(I) = SIGNXX(I)*SIGNXX(I)
273ctmp BB(I) = THREE*SIGNXY(I)*SIGNXY(I)
274ctmp CC(I) = THREE*SIGNXZ(I)*SIGNXZ(I)
275ctmp SVM(I)= AA(I) + BB(I) + CC(I)
276ctmp ENDDO
277ctmpC plastic flow
278ctmp NINDX=0
279ctmp DO I=1,NEL
280ctmp IF(SVM(I)>YLD(I).AND.OFF(I)==1.) THEN
281ctmp NINDX = NINDX+1
282ctmp INDEX(NINDX) = I
283ctmp ENDIF
284ctmp ENDDO
285ctmpC def plastiques en contrainte uniaxiale
286ctmp DO J=1,NINDX
287ctmp I=INDEX(J)
288ctmp DPLA_J(I) = (SVM(I)-YLD(I))/(G3(I)+H(I))
289ctmp ETSE(I) = H(I)/(H(I)+E(I))
290ctmp ENDDO
291ctmpC 3 iterations Newton
292ctmp DO N=1,NITER
293ctmp DO J=1,NINDX
294ctmp I = INDEX(J)
295ctmp DPLA_I(I) = DPLA_J(I)
296ctmp YLD_I = YLD(I) + H(I)*DPLA_I(I)
297ctmp DR(I) = E(I)*DPLA_I(I)/YLD_I
298ctmp PP(I) = ONE/(ONE+DR(I))
299ctmp QQ(I) = ONE/(ONE+NU1(I)*DR(I))
300ctmp P2 = PP(I)*PP(I)
301ctmp Q2 = QQ(I)*QQ(I)
302ctmp F = AA(I)*P2 + (BB(I)+CC(I))*Q2 - YLD_I*YLD_I
303ctmp DF = AA(I)*P2*PP(I) + NU*(BB(I)+CC(I))*Q2*QQ(I)
304ctmp DF = DF *(H(I)*DR(I) - E(I)) / YLD_I - H(I)*YLD_I
305ctmp DF = DF * TWO
306ctmp DF = SIGN(MAX(ABS(DF),EM20),DF)
307ctmp IF(DPLA_I(I) > ZERO) THEN
308ctmp DPLA_J(I) = MAX(ZERO,DPLA_I(I)-F/DF)
309ctmp ELSE
310ctmp DPLA_J(I) = ZERO
311ctmp ENDIF
312ctmp ENDDO
313ctmp ENDDO
314ctmpC Contraintes plastiquement admissibles
315ctmp DO J=1,NINDX
316ctmp I=INDEX(J)
317ctmp PLA(I) = PLA(I) + DPLA_I(I)
318ctmp SIGNXX(I) = SIGNXX(I)*PP(I)
319ctmp SIGNXY(I) = SIGNXY(I)*QQ(I)
320ctmp SIGNXZ(I) = SIGNXZ(I)*QQ(I)
321ctmp ENDDO
322C---
323c ENDIF
324C---
325 IF (ifail > 0) THEN
326 DO i=1,nel
327 IF (pla(i) > epsmax(i) .AND. off(i)==one) THEN
328 off(i) = four_over_5
329#include "lockon.inc"
330 WRITE(iout, 1000) ngl(i),pla(i),epsp(i)
331#include "lockoff.inc"
332 END IF
333 ENDDO
334 END IF
335c---------------------------------------------------------
336 1000 FORMAT(5x,'FAILURE BEAM ELEMENT NUMBER ',i10,
337 . ' PLASTIC STRAIN =',1pe16.9,' STRAIN RATE =',1pe16.9)
338c---------------------------------------------------------
339 RETURN
340 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps36pi
subroutine sigeps36pi(nel, kfunc, ipt, npf, tf, ngl, time, timestep, uparam, ipm, imat, e, g, off, depsxx, depsxy, depsxz, sigoxx, sigoxy, sigoxz, signxx, signxy, signxz, etse, pla, epsp, nvartmp, vartmp, sigy)
Definition sigeps36pi.F:38
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72