OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps116.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!|| sigeps116 ../engine/source/materials/mat/mat116/sigeps116.F
25!||--- called by ------------------------------------------------------
26!|| suser43 ../engine/source/elements/solid/sconnect/suser43.F
27!||====================================================================
28 SUBROUTINE sigeps116(
29 1 NEL ,NUPARAM ,NUVAR ,JSMS ,TIME ,TIMESTEP,
30 2 UPARAM ,UVAR ,AREA ,EPSD ,OFF ,OFFL ,
31 3 EPSZZ ,EPSYZ ,EPSZX ,DEPSZZ ,DEPSYZ ,DEPSZX ,
32 4 SIGNZZ ,SIGNYZ ,SIGNZX ,STIFM ,DMELS ,DMG ,
33 5 PLA_N ,PLA_T ,IPG ,NFAIL ,NGL )
34C-----------------------------------------------
35C I m p l i c i t T y p e s
36C-----------------------------------------------
37#include "implicit_f.inc"
38C----------------------------------------------------------
39#include "units_c.inc"
40#include "comlock.inc"
41#include "sms_c.inc"
42C----------------------------------------------------------
43C D u m m y A r g u m e n t s
44C----------------------------------------------------------
45 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,JSMS,IPG
46 INTEGER ,INTENT(OUT) :: NFAIL
47 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
48 my_real ,INTENT(IN) :: TIME,TIMESTEP
49 my_real ,DIMENSION(NEL) :: OFF,OFFL,PLA_N,PLA_T,AREA,DMELS,
50 . epszz,epsyz,epszx,depszz,depsyz,depszx,signzz,signyz,signzx
51 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: dmg
52 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: epsd,stifm
53 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) :: uparam
54 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
55C-----------------------------------------------
56C L o c a l V a r i a b l e s
57C-----------------------------------------------
58 INTEGER :: II,IEL,IORDER1,IORDER2,IFAIL1,IFAIL2,ICRIT,NINDXF,NINDXD
59 INTEGER ,DIMENSION(NEL) :: INDXF,INDXD
60 my_real :: RHO0,E,G,GC1_INI,GC2_INI,GC1_INF,GC2_INF,RATG1,RATG2,
61 . FG1,FG2,SIGA1,SIGA2,SIGB1,SIGB2,RATE1,RATE2,DTINV,EPSDOT,
62 . ALPHA,ALPHAI,BETA,DM,DPN,DPT,DPT1,DPT2,DAM,DEPS,STF,ET2,
63 . bb,dtb,fact1,fact2,thick
64 my_real, DIMENSION(NEL) :: yld_n,yld_t,gc_n,gc_t,epsm1,epsm2,epsmf,
65 . cosg,sing,epsn,epst,epsm,epsn1,epsn2,epst1,epst2,epsp,
66 . pla_t1,pla_t2
67c-----------------------------------------------
68c UVAR(1) = PLA_N : plastic separation strain in normal direction
69c UVAR(2) = PLA_T1 : plastic separation strain in shear direction XZ
70c UVAR(3) = PLA_T2 : plastic separation strain in shear direction YZ
71c UVAR(4) = EPSP : strain rate
72c UVAR(5) = STRFLAG : strain rate update flag (0 => recalculate, 1 => const)
73c UVAR(6) = DAM : damage coefficient
74c check only / output variables
75c UVAR(7) = EPSN
76c UVAR(8) = EPST
77c UVAR(9) = EPSM
78c UVAR(10) = EPSM1 : mixed mode yield initiation strain
79c UVAR(11) = EPSM2 : mixed mode damage initiation strain
80c UVAR(12) = EPSMF : mixed mode damage initiation strain
81C=======================================================================
82C INPUT PARAMETERS INITIALIZATION
83C-----------------------------------------------
84 e = uparam(1)
85 g = uparam(2)
86 nfail = uparam(4)
87 gc1_ini = uparam(5)
88 gc1_inf = uparam(6)
89 ratg1 = uparam(7)
90 fg1 = uparam(8)
91 gc2_ini = uparam(9)
92 gc2_inf = uparam(10)
93 ratg2 = uparam(11)
94 fg2 = uparam(12)
95 siga1 = uparam(13)
96 sigb1 = uparam(14)
97 rate1 = uparam(15)
98 iorder1 = uparam(16)
99 ifail1 = uparam(17)
100 siga2 = uparam(18)
101 sigb2 = uparam(19)
102 rate2 = uparam(20)
103 iorder2 = uparam(21)
104 ifail2 = uparam(22)
105 icrit = uparam(23)
106 thick = uparam(24)
107 alpha = uparam(25) ! filtering coefficient of strain rate (exponential avg)
108 alphai = one - alpha
109 dtinv = one / (max(timestep, em20))
110 stf = e + g
111 nindxf = 0
112 nindxd = 0
113c-------------------------
114 pla_n(:nel) = uvar(:nel,1)
115 pla_t1(:nel) = uvar(:nel,2)
116 pla_t2(:nel) = uvar(:nel,3)
117c-------------------------
118 DO iel=1,nel
119 et2 = epsyz(iel)**2 + epszx(iel)**2
120 epsn(iel) = max(epszz(iel), zero)
121 epst(iel) = sqrt(et2)
122 epsm(iel) = sqrt(et2 + epsn(iel)**2)
123 ENDDO
124c-------------------------
125c Strain rate
126c-------------------------
127 DO iel=1,nel
128 IF (uvar(iel,5) == zero) THEN ! strain rate calculation only when epsm < EPSM1
129 epsdot = sqrt(depsyz(iel)**2 + depszx(iel)**2 + depszz(iel)**2)
130 epsdot = epsdot * dtinv / thick
131 epsp(iel) = alpha *epsdot + alphai * uvar(iel,4)
132 epsd(iel) = epsp(iel)
133 uvar(iel,4) = epsp(iel)
134 ELSE ! constant strain rate
135 epsp(iel) = uvar(iel,4)
136 END IF
137 ENDDO
138c-------------------------
139c Strain rate dependent yield values
140c-------------------------
141 yld_n(:nel) = siga1
142 yld_t(:nel) = siga2
143 IF (sigb1 > zero) THEN ! normal yld is strain rate dependent
144 IF (iorder1 == 1) THEN ! linear log dependency
145 DO iel=1,nel
146 IF (epsp(iel) > rate1) THEN
147 yld_n(iel) = siga1 + sigb1*log(epsp(iel)/rate1)
148 END IF
149 ENDDO
150 ELSE IF (iorder1 == 2) THEN ! quadratic log dependency
151 DO iel=1,nel
152 IF (epsp(iel) > rate1) THEN
153 yld_n(iel) = siga1 + sigb1*log(epsp(iel)/rate1)**2
154 END IF
155 ENDDO
156 END IF
157 END IF
158c
159 IF (sigb2 > zero) THEN ! tangent yld is strain rate dependent
160 IF (iorder2 == 1) THEN ! linear log dependency
161 DO iel=1,nel
162 IF (epsp(iel) > rate2) THEN
163 yld_t(iel) = siga2 + sigb2*log(epsp(iel)/rate2)
164 END IF
165 ENDDO
166 ELSE IF (iorder2 == 2) THEN ! quadratic log dependency
167 DO iel=1,nel
168 IF (epsp(iel) > rate2) THEN
169 yld_t(iel) = siga2 + sigb2*log(epsp(iel)/rate2)**2
170 END IF
171 ENDDO
172 END IF
173 END IF
174c-------------------------
175c Strain rate dependent fracture energies
176c-------------------------
177 gc_n(:nel) = gc1_ini
178 gc_t(:nel) = gc2_ini
179 IF (gc1_inf > zero) THEN ! normal GC is strain rate dependent
180 DO iel=1,nel
181! IF (EPSP(IEL) > RATG1) THEN
182 IF (epsp(iel) > zero) THEN
183 gc_n(iel) = gc1_ini + (gc1_inf-gc1_ini)*exp(-ratg1/epsp(iel))
184 END IF
185 ENDDO
186 END IF
187c
188 IF (gc2_inf > zero) THEN ! tangent GC is strain rate dependent
189 DO iel=1,nel
190! IF (EPSP(IEL) > RATG2) THEN
191 IF (epsp(iel) > zero) THEN
192 gc_t(iel) = gc2_ini + (gc2_inf-gc2_ini)*exp(-ratg2/epsp(iel))
193 END IF
194 ENDDO
195 END IF
196c-------------------------
197 DO iel=1,nel
198 epsn1(iel) = yld_n(iel) / e
199 epst1(iel) = yld_t(iel) / g
200 END DO
201c-------------------------
202c Check max FG values
203c-------------------------
204 IF (ifail1 == 1) THEN
205 DO iel=1,nel
206 epsn2(iel) = epsn1(iel) + fg1 * gc_n(iel)/yld_n(iel)
207 END DO
208 ELSE
209 DO iel=1,nel
210 epsn2(iel) = (epsn1(iel) + two*fg1*gc_n(iel) / yld_n(iel)) / (fg1 + one)
211 END DO
212 END IF
213c
214 IF (ifail2 == 1)THEN
215 DO iel=1,nel
216 epst2(iel) = epst1(iel) + fg2 * gc_t(iel)/yld_t(iel)
217 END DO
218 ELSE
219 DO iel=1,nel
220 epst2(iel) = (epst1(iel) + two*fg2*gc_t(iel) / yld_t(iel)) / (fg2 + one)
221 END DO
222 END IF
223c--------------------------------
224c Update failure initiation strains in tension, shear and mixed mode
225c--------------------------------
226 IF (icrit == 1) THEN
227 DO iel=1,nel
228 IF (epst(iel) == zero) THEN
229 epsm1(iel) = epsn1(iel)
230 epsm2(iel) = epsn2(iel)
231 ELSE IF (epsn(iel) == zero) THEN
232 epsm1(iel) = epst1(iel)
233 epsm2(iel) = epst2(iel)
234 ELSE ! mixed mode
235 beta = epst(iel) / epsn(iel)
236 fact1 = sqrt((one+beta**2) /(epst1(iel)**2 + (beta*epsn1(iel))**2))
237 fact2 = sqrt((one+beta**2) /(epst2(iel)**2 + (beta*epsn2(iel))**2))
238 epsm1(iel) = epsn1(iel)*epst1(iel)*fact1
239 epsm2(iel) = epsn2(iel)*epst2(iel)*fact2
240 END IF
241 IF (epsm(iel) > epsm1(iel)) uvar(iel,5) = one
242 END DO
243 ELSE
244 DO iel=1,nel
245 IF (epst(iel) == zero) THEN
246 epsm1(iel) = epsn1(iel)
247 epsm2(iel) = epsn2(iel)
248 ELSE IF (epsn(iel) == zero) THEN
249 epsm1(iel) = epst1(iel)
250 epsm2(iel) = epst2(iel)
251 ELSE ! mixed mode
252 beta = epst(iel) / epsn(iel)
253 IF (beta*epsn1(iel) > epst1(iel)) THEN
254 epsm1(iel) = epst1(iel)*sqrt(one + beta**2) / beta
255 ELSE
256 epsm1(iel) = epsn1(iel)*sqrt(one + beta**2)
257 END IF
258 IF (beta*epsn2(iel) > epst2(iel)) THEN
259 epsm2(iel) = epst2(iel)*sqrt(one + beta**2) / beta
260 ELSE
261 epsm2(iel) = epsn2(iel)*sqrt(one + beta**2)
262 END IF
263 END IF
264 IF (epsm(iel) > epsm1(iel)) uvar(iel,5) = one
265 END DO
266 END IF
267c--------------------------------
268c Power law of damage evolution
269c--------------------------------
270 DO iel=1,nel
271 IF (epsn(iel) == zero) THEN
272 cosg(iel) = zero
273 sing(iel) = one
274 ELSE IF (epst(iel) == zero) THEN
275 cosg(iel) = one
276 sing(iel) = zero
277 ELSE IF (epsm(iel) > zero) THEN
278 cosg(iel) = epszz(iel) / epsm(iel)
279 sing(iel) = epst(iel) / epsm(iel)
280 END IF
281 IF (epsm(iel) > zero) THEN
282 dm = epsm1(iel) - epsm2(iel)
283 bb = epsm1(iel)*(e*gc_t(iel)*cosg(iel)**2 + g*gc_n(iel)*sing(iel)**2)
284 epsmf(iel) = dm + two*gc_t(iel)*gc_n(iel) / bb
285 ELSE
286 epsmf(iel) = ep20
287 END IF
288 END DO
289c
290 DO iel=1,nel
291 dm = epsm(iel) - epsm2(iel)
292 IF (off(iel) > zero .and. dm > zero) THEN
293 IF (uvar(iel,6) == zero) THEN
294 nindxd = nindxd+1
295 indxd(nindxd) = iel
296 END IF
297 dam = dm / max((epsmf(iel) - epsm2(iel)), em20)
298 dam = max(dam, uvar(iel,6))
299 uvar(iel,6) = dam
300 dmg(iel) = max(dmg(iel), dam)
301 dmg(iel) = min(dmg(iel), one)
302 IF (off(iel)*offl(iel) == one .and. epsm(iel) > epsmf(iel)) THEN
303 nindxf = nindxf+1
304 indxf(nindxf) = iel
305 offl(iel) = zero
306 END IF
307 END IF
308 END DO
309c-------------------------
310c Plastic strains
311c-------------------------
312 DO iel=1,nel
313 IF (off(iel)*offl(iel) == one) THEN
314 dpn = epsn(iel) - epsm1(iel)*cosg(iel)
315 IF (dpn > zero) THEN
316 pla_n(iel) = max(dpn, pla_n(iel))
317 uvar(iel,1) = pla_n(iel)
318 END IF
319 dpt1 = epsyz(iel) - pla_t1(iel)
320 dpt2 = epszx(iel) - pla_t2(iel)
321 dpt = sqrt(dpt1**2 + dpt2**2)
322 IF (dpt > epsm1(iel)*sing(iel)) THEN
323 pla_t1(iel) = pla_t1(iel) + depsyz(iel)
324 pla_t2(iel) = pla_t2(iel) + depszx(iel)
325 uvar(iel,2) = pla_t1(iel)
326 uvar(iel,3) = pla_t2(iel)
327 pla_t(iel) = sqrt(pla_t1(iel)**2 + pla_t2(iel)**2)
328 END IF
329 END IF
330 END DO
331c-------------------------
332c Stress update
333c-------------------------
334 DO iel=1,nel
335 signzz(iel) = e * (epszz(iel) - pla_n(iel))*(one-dmg(iel))*off(iel)
336 signyz(iel) = g * (epsyz(iel) - pla_t1(iel))*(one-dmg(iel))*off(iel)
337 signzx(iel) = g * (epszx(iel) - pla_t2(iel))*(one-dmg(iel))*off(iel)
338c
339 uvar(iel,7) = epsn(iel)
340 uvar(iel,8) = epst(iel)
341 uvar(iel,9) = epsm(iel)
342 uvar(iel,10) = epsm1(iel)
343 uvar(iel,11) = epsm1(iel)
344 uvar(iel,12) = epsmf(iel)
345 ENDDO
346c-----------------------------------------------------
347c omega = sqrt(2k/2*dmels), dt=2/omega, 2*dmels=dt**2 * 2k / 4
348 IF (idtmins==2 .AND. jsms/=0) THEN
349 dtb = (dtmins/dtfacs)**2
350 DO iel=1,nel
351 dmels(iel)=max(dmels(iel),half*dtb*stf*area(iel)*off(iel))
352 ENDDO
353 END IF
354 stifm(1:nel) = stifm(1:nel) + stf*area(1:nel)*off(1:nel)
355c-----------------------------------------------------
356 IF (nindxd > 0) THEN
357 DO ii=1,nindxd
358 iel = indxd(ii)
359#include "lockon.inc"
360 WRITE(iout ,1000) ngl(iel),ipg,epsm(iel)
361#include "lockoff.inc"
362 END DO
363 END IF
364c
365 IF (nindxf > 0) THEN
366 DO ii=1,nindxf
367 iel = indxf(ii)
368#include "lockon.inc"
369 WRITE(iout ,2000) ngl(iel),ipg,epsm(iel)
370 WRITE(istdo,2100) ngl(iel),ipg,epsm(iel),time
371#include "lockoff.inc"
372 END DO
373 END IF
374c-----------------------------------------------------
375 1000 FORMAT(5x,'START DAMAGE COHESIVE ELEMENT ',i10,
376 . ' INTEGRATION POINT',i2,', MIXED MODE STRAIN=',1pe16.9)
377 2000 FORMAT(5x,'FAILURE COHESIVE ELEMENT ',i10,
378 . ' INTEGRATION POINT',i2,', MIXED MODE STRAIN=',1pe16.9)
379 2100 FORMAT(5x,'FAILURE COHESIVE ELEMENT ',i10,
380 . ' INTEGRATION POINT',i2,', MIXED MODE STRAIN=',1pe16.9,
381 . ' AT TIME ',1pe16.9)
382c-----------
383 RETURN
384 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps116
subroutine sigeps116(nel, nuparam, nuvar, jsms, time, timestep, uparam, uvar, area, epsd, off, offl, epszz, epsyz, epszx, depszz, depsyz, depszx, signzz, signyz, signzx, stifm, dmels, dmg, pla_n, pla_t, ipg, nfail, ngl)
Definition sigeps116.F:34