OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps44p.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!|| sigeps44p ../engine/source/materials/mat/mat044/sigeps44p.F
25!||--- called by ------------------------------------------------------
26!|| main_beam3 ../engine/source/elements/beam/main_beam3.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||====================================================================
30 SUBROUTINE sigeps44p(
31 . NEL ,NGL ,MAT ,PID ,NUPARAM ,UPARAM ,
32 . GEO ,OFF ,PLA ,AL ,
33 . EXX ,EXY ,EXZ ,KXX ,KYY ,KZZ ,
34 . FA1 ,FA2 ,FA3 ,MA1 ,MA2 ,MA3 ,
35 . FOR ,MOM ,PM ,NUVAR ,UVAR ,NFUNC ,
36 . IFUNC ,TF ,NPF ,SIGY)
37C-----------------------------------------------
38C I m p l i c i t T y p e s
39C-----------------------------------------------
40#include "implicit_f.inc"
41#include "comlock.inc"
42C-----------------------------------------------
43C C o m m o n B l o c k s
44C-----------------------------------------------
45#include "param_c.inc"
46#include "scr17_c.inc"
47#include "com04_c.inc"
48#include "com08_c.inc"
49#include "units_c.inc"
50#include "impl1_c.inc"
51C-----------------------------------------------
52C D u m m y A r g u m e n t s
53C-----------------------------------------------
54 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,NFUNC,IFUNC(NFUNC),
55 . NPF(*)
56 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: MAT,PID,NGL
57 my_real ,DIMENSION(NPROPM ,NUMMAT) ,INTENT(IN) :: PM
58 my_real ,DIMENSION(NPROPG ,NUMGEO) ,INTENT(IN) :: GEO
59 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) :: UPARAM
60 my_real ,DIMENSION(NEL,NUVAR) :: uvar
61 my_real ,DIMENSION(NEL,3) :: for,mom
62 my_real ,DIMENSION(NEL) :: off,pla,al,exx,exy,exz,kxx,kyy,kzz,
63 . fa1,fa2,fa3,ma1,ma2,ma3
64 my_real
65 . tf(*),finter
66 EXTERNAL finter
67 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: sigy
68C-----------------------------------------------
69C L o c a l V a r i a b l e s
70C-----------------------------------------------
71 INTEGER ,DIMENSION(NEL) :: INDX,ICC,ISRATE,VFLAG,IPOS,IAD,ILEN
72 INTEGER :: I,J,IPID,NINDX
73 my_real :: dmg,fact,epif,frate,alpha,epsdot
74 my_real ,DIMENSION(NEL) :: e,g,g3,ca,cb,cn,cp,a1,b1,b2,b3,shf,
75 . degmb,degfx,dmpm,dmpf,f1,m1,m2,m3,degsh,yeq,yld,yma2,yldmax,
76 . dwpla,epmax,epst,epsr1,epsr2,epdr,epsp,rho,epmx,dwelm,dwelf,rr,
77 . sh,sh10,sh20,sh0,sh1,sh2,asrate,yscale,dpla,dfdpla
78C=======================================================================
79 epif = zero
80C
81 IF (impl_s == 0 .OR. idyna > 0) THEN
82 DO i=1,nel
83 dmpm(i)=geo(16,pid(i))*al(i)
84 dmpf(i)=geo(17,pid(i))*al(i)
85 ENDDO
86 ELSE
87 DO i=1,nel
88 dmpm(i)=zero
89 dmpf(i)=zero
90 ENDDO
91 ENDIF
92C
93 DO i=1,nel
94 ipid = pid(i)
95 e(i) = uparam(1)
96 ca(i) = uparam(3)
97 yldmax(i)= uparam(4)
98 epmax(i) = uparam(5)
99 epsr1(i) = uparam(6)
100 epsr2(i) = uparam(7)
101 cb(i) = uparam(8)
102 cn(i) = uparam(9)
103 icc(i) = nint(uparam(10))
104 epdr(i) = uparam(11)
105 epif = max(epif,epdr(i))
106 cp(i) = uparam(12)
107 g(i) = uparam(16)
108 g3(i) = uparam(18)
109 israte(i)= nint(uparam(13))
110 asrate(i)= uparam(14)
111 vflag(i) = nint(uparam(23))
112 yscale(i)= uparam(24)
113c
114 rho(i) = pm(1,mat(i))
115c
116 a1(i) = geo(1 ,ipid)
117 b1(i) = geo(2 ,ipid)
118 b2(i) = geo(18,ipid)
119 b3(i) = geo(4 ,ipid)
120 shf(i) = geo(37,ipid)
121 dpla(i) = zero
122 ENDDO
123C-----------------------------
124C
125 DO i=1,nel
126 degmb(i) = for(i,1)*exx(i)
127 degsh(i) = for(i,2)*exy(i)+for(i,3)*exz(i)
128 degfx(i) = mom(i,1)*kxx(i)+mom(i,2)*kyy(i)+mom(i,3)*kzz(i)
129 ENDDO
130C
131C CISAILLEMENT TRANSVERSAL CALCULE AVEC K1=12EI/L**2 K2=5/6GA
132C
133 DO i=1,nel
134 sh(i) = five_over_6*g(i)*a1(i)
135 yma2(i)= twelve*e(i)/al(i)**2
136 sh10(i)= yma2(i)*b1(i)
137 sh20(i)= yma2(i)*b2(i)
138 sh0(i) = (one - shf(i))*sh(i)
139 sh1(i) = sh0(i)*sh10(i)/(sh(i)+sh10(i)) + shf(i)*sh10(i)
140 sh2(i) = sh0(i)*sh20(i)/(sh(i)+sh20(i)) + shf(i)*sh20(i)
141C
142 f1(i) = for(i,1) + exx(i)*a1(i)*e(i)
143 for(i,2)= for(i,2) + exy(i)*sh2(i)
144 for(i,3)= for(i,3) + exz(i)*sh1(i)
145 m1(i) = mom(i,1) + kxx(i)*g(i) *b3(i)
146 m2(i) = mom(i,2) + kyy(i)*e(i)*b1(i)
147 m3(i) = mom(i,3) + kzz(i)*e(i)*b2(i)
148 epst(i) = f1(i)/a1(i)/e(i)
149 ENDDO
150C-------------
151C CRITERE
152C-------------
153 DO i=1,nel
154 yeq(i) = f1(i)*f1(i) + three * a1(i) *
155 . (m1(i)*m1(i) / max(b3(i),em20)
156 . + m2(i)*m2(i) / max(b1(i),em20)
157 . + m3(i)*m3(i) / max(b2(i),em20) )
158 yeq(i) = max(em20, sqrt(yeq(i)) / a1(i) )
159 ENDDO
160C-----------------------------------
161C YIELD
162C-----------------------------------
163 IF (nfunc > 0) THEN
164 !IPOS(1:NEL) = VARTMP(1:NEL,1)
165 !IAD (1:NEL) = NPF(IFUNC(1)) / 2 + 1
166 !ILEN(1:NEL) = NPF(IFUNC(1)+1) / 2 - IAD(1:NEL) - IPOS(1:NEL)
167 !CALL VINTER(TF,IAD,IPOS,ILEN,NEL,PLA,DFDPLA,YLD)
168 DO i = 1, nel
169 yld(i) = yscale(i)*finter(ifunc(1),pla(i),npf,tf,dfdpla(i))
170 ENDDO
171 !VARTMP(1:NEL,1) = IPOS(1:NEL)
172 ENDIF
173c
174c
175c
176 DO i = 1,nel
177 IF (nfunc > 0) THEN
178 yld(i) = yscale(i)*yld(i)
179 ELSE
180 yld(i) = ca(i)
181 ENDIF
182 ENDDO
183c
184c
185c
186 DO i = 1,nel
187 IF (pla(i) > zero) THEN
188 IF (nfunc > 0) THEN
189 yld(i) = yscale(i)*yld(i)
190 ELSE
191 yld(i) = ca(i) + cb(i)*exp(cn(i)*log(pla(i)))
192 ENDIF
193 ENDIF
194 ENDDO
195C-------------
196C STRAIN RATE EFFECT
197C-------------
198 IF (epif > zero) THEN
199 DO i = 1,nel
200 IF (epdr(i) > zero) THEN
201 epsp(i) = abs(degmb(i) + degfx(i))/yeq(i)/a1(i)
202 IF (vflag(i) /= 1) THEN
203 IF (israte(i) == 1) THEN
204 alpha = min(one, asrate(i)*dt1)
205 epsdot = alpha*abs(epsp(i)/max(em20,dt1)) + (one-alpha)*uvar(i,1)
206 uvar(i,1) = epsdot
207 ELSE
208 epsdot = abs(epsp(i)/max(em20,dt1))
209 ENDIF
210 ELSE
211 epsdot = uvar(i,1)
212 ENDIF
213 frate = one + (epsdot*epdr(i))**cp(i)
214 IF (icc(i)== 1) yldmax(i) = yldmax(i) * frate
215 IF ((nfunc > 0) .AND. (ca(i) /= zero)) THEN
216 IF (pla(i)>zero) THEN
217 yld(i) = yld(i) + (ca(i) + cb(i)*exp(cn(i)*log(pla(i))))*(frate-one)
218 ELSE
219 yld(i) = yld(i) + ca(i)*(frate-one)
220 ENDIF
221 ELSE
222 yld(i) = yld(i) * frate
223 ENDIF
224 ENDIF
225 ENDDO
226 ELSE
227 epsp(1:nel )= one
228 ENDIF
229c-------------------
230c PROJECTION
231c-------------------
232 DO i=1,nel
233 yld(i) = min(yld(i),yldmax(i))
234 sigy(i) = yld(i)
235 rr(i) = min(one, yld(i) / yeq(i))
236 f1(i) = f1(i)*rr(i)
237 dwelm(i) =(f1(i) + for(i,1))*(f1(i)-for(i,1)) / e(i)/a1(i)
238 degmb(i) = degmb(i)+ f1(i)*exx(i)
239 ENDDO
240C
241 DO i=1,nel
242 m1(i) = m1(i)*rr(i)
243 m2(i) = m2(i)*rr(i)
244 m3(i) = m3(i)*rr(i)
245 dwelf(i) =(m1(i)+mom(i,1))*(m1(i)-mom(i,1))/ g(i)/b3(i)+
246 . (m2(i)+mom(i,2))*(m2(i)-mom(i,2))/e(i)/b1(i)+
247 . (m3(i)+mom(i,3))*(m3(i)-mom(i,3))/e(i)/b2(i)
248 degfx(i) = degfx(i) + m1(i)*kxx(i) + m2(i)*kyy(i) + m3(i)*kzz(i)
249 ENDDO
250C
251 DO i=1,nel
252 dwpla(i) = degmb(i) + degfx(i) - dwelm(i) - dwelf(i)
253 degsh(i) = degsh(i) + for(i,2)*exy(i) + for(i,3)*exz(i)
254 fact = half*off(i)*al(i)
255 ENDDO
256C-----------------------
257C EPS PLASTIQUE
258C-----------------------
259 DO i=1,nel
260 fact = dwpla(i)/yld(i)
261 dpla(i) = off(i)*max(zero,half*fact/a1(i))
262 pla(i) = pla(i) + off(i)*max(zero,half*fact/a1(i))
263 ENDDO
264c--------------------------------
265c DUCTILE RUPTURE
266c--------------------------------
267 DO i=1,nel
268 IF (off(i) < em01) off(i) = zero
269 IF (off(i) < one ) off(i) = off(i)*four_over_5
270 ENDDO
271c--------------------------------
272c AXIAL TENSION OR PLASTIC STRAIN FAILURE
273c--------------------------------
274 nindx = 0
275 DO i = 1,nel
276 IF (off(i) == one) THEN
277 dmg = one
278 IF (epst(i) > epsr1(i)) THEN
279 dmg = (epsr2(i) - epst(i)) / (epsr2(i) - epsr1(i))
280 dmg = max(dmg, zero)
281 for(i,1) = f1(i)*dmg
282 for(i,2) = for(i,2)*dmg
283 for(i,3) = for(i,3)*dmg
284 mom(i,1) = m1(i)*dmg
285 mom(i,2) = m2(i)*dmg
286 mom(i,3) = m3(i)*dmg
287 ENDIF
288c test strain failure
289 IF (dmg == zero .or. pla(i) >= epmax(i)) THEN
290 off(i) = four_over_5
291 idel7nok = 1
292 nindx = nindx+1
293 indx(nindx) = i
294 ENDIF
295c
296 IF (vflag(i) == 1) THEN
297 alpha = min(one, asrate(i)*dt1)
298 epsdot = dpla(i)/max(em20,dt1)
299 uvar(i,1) = alpha*epsdot + (one - alpha)*uvar(i,1)
300 ENDIF
301c
302 ENDIF
303 ENDDO
304c--------------------------------
305 IF (nindx > 0 .AND. imconv == 1) THEN
306 DO j=1,nindx
307 i = indx(j)
308#include "lockon.inc"
309 WRITE(iout, 1000) ngl(i)
310 WRITE(istdo,1100) ngl(i),tt
311#include "lockoff.inc"
312 ENDDO
313 ENDIF
314C
315 DO i=1,nel
316 fa1(i) = f1(i)
317 fa2(i) = for(i,2)
318 fa3(i) = for(i,3)
319 ma1(i) = m1(i)
320 ma2(i) = m2(i)
321 ma3(i) = m3(i)
322 ENDDO
323c-----------------------------------------------
324 1000 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT NUMBER ',i10)
325 1100 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT :',i10,' AT TIME :',g11.4)
326c-----------------------------------------------
327 RETURN
328 END
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
sigeps44p
subroutine sigeps44p(nel, ngl, mat, pid, nuparam, uparam, geo, off, pla, al, exx, exy, exz, kxx, kyy, kzz, fa1, fa2, fa3, ma1, ma2, ma3, for, mom, pm, nuvar, uvar, nfunc, ifunc, tf, npf, sigy)
Definition sigeps44p.F:37