OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_wilkins_s.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!|| fail_wilkins_s ../engine/source/materials/fail/wilkins/fail_wilkins_s.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!|| mmain8 ../engine/source/materials/mat_share/mmain8.F
28!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
29!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
30!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
31!||====================================================================
32 SUBROUTINE fail_wilkins_s(
33 1 NEL ,NUPARAM,NUVAR ,
34 2 TIME ,TIMESTEP ,UPARAM ,NGL ,
35 4 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
36 5 DPLA ,UVAR ,OFF ,DFMAX ,TDELE )
37C-----------------------------------------------
38C 3D Wilkins Failure model
39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C---------+---------+---+---+--------------------------------------------
44C VAR | SIZE |TYP| RW| DEFINITION
45C---------+---------+---+---+--------------------------------------------
46C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
47C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
48C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
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 EPSPXX | NEL | F | R | STRAIN RATE XX
54C EPSPYY | NEL | F | R | STRAIN RATE YY
55C ... | | | |
56C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
57C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
58C ... | | | |
59C EPSXX | NEL | F | R | STRAIN XX
60C EPSYY | NEL | F | R | STRAIN 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 SIGVXX | NEL | F | W | VISCOUS STRESS XX
70C SIGVYY | NEL | F | W | VISCOUS STRESS YY
71C ... | | | |
72C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
73C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
74C---------+---------+---+---+--------------------------------------------
75C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
76C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
77C---------+---------+---+---+--------------------------------------------
78#include "mvsiz_p.inc"
79#include "scr17_c.inc"
80#include "units_c.inc"
81#include "comlock.inc"
82#include "param_c.inc"
83#include "impl1_c.inc"
84C-----------------------------------------------
85C I N P U T A r g u m e n t s
86C-----------------------------------------------
87C
88 INTEGER NEL, NUPARAM, NUVAR,NGL(NEL)
89 my_real TIME,TIMESTEP,UPARAM(*),
90 . SIGNXX(NEL),SIGNYY(NEL),SIGNZZ(NEL),
91 . SIGNXY(NEL),SIGNYZ(NEL),SIGNZX(NEL),
92 . dpla(nel)
93C-----------------------------------------------
94C O U T P U T A r g u m e n t s
95C-----------------------------------------------
96cc my_real
97
98C-----------------------------------------------
99C I N P U T O U T P U T A r g u m e n t s
100C-----------------------------------------------
101 my_real uvar(nel,nuvar), off(nel),dfmax(nel),tdele(nel)
102C-----------------------------------------------
103C L o c a l V a r i a b l e s
104C-----------------------------------------------
105 INTEGER I,J,IDEL,IDEV,IFLAG,INDX(MVSIZ),IADBUF,NINDX,
106 . NINDEX,INDEX(MVSIZ),JJ,IFAIL
107 my_real
108 . AL,BE,PC,DC,
109 . p,scale,s1,s2,s3,ss,e1,e2,e3,e4,e5,e6,
110 . epst2, a, cc1, b, y ,yp, d, e42, e52, c, e62,w1,w2
111C-----------------------------------------------
112 al = uparam(1)
113 be = uparam(2)
114 pc = uparam(3)
115 dc = uparam(4)
116 iflag = int(uparam(6))
117C-----------------------------------------------
118 idel=0
119 idev=0
120 scale = zero
121 IF(iflag==1)THEN
122 idel=1
123 ELSEIF(iflag==2)THEN
124 idev =1
125 END IF
126C...
127 IF(idel==1)THEN
128 DO i=1,nel
129 IF(off(i)<0.1) off(i)=0.0
130 IF(off(i)<1.0) off(i)=off(i)*0.8
131 END DO
132 END IF
133C
134 IF(idel==1)THEN
135 nindx=0
136C-------------------------------
137C RUPTURE DUCTILE
138C-------------------------------
139C Tuler Butcher
140 DO i=1,nel
141 IF(iflag==1.AND.off(i)==1.)THEN
142C-------------------
143C STREES principal 1, 4 newton iterations
144C-------------------
145 p = third*(signxx(i) + signyy(i) + signzz(i))
146 e1 = signxx(i) - p
147 e2 = signyy(i) - p
148 e3 = signzz(i) - p
149 e4 = signxy(i)
150 e5 = signyz(i)
151 e6 = signzx(i)
152C -y = (e1-x)(e2-x)(e3-x)
153C - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
154C + 2e4 e5 e6
155C e1 + e2 + e3 = 0 => terme en x^2 = 0
156C y = x^3 + c x + d
157c yp= 3 x^2 + c
158 e42 = e4*e4
159 e52 = e5*e5
160 e62 = e6*e6
161 c = - half * (e1*e1 + e2*e2 + e3*e3) - e42 - e52 - e62
162 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
163 & - two*e4*e5*e6
164 cc1 = c*third
165 s1 = sqrt(-cc1)
166 epst2 = s1 * s1
167 y = (epst2 + c)* s1 + d
168 IF(abs(y)>em8)THEN
169 s1 = 1.75 * s1
170 epst2 = s1 * s1
171 y = (epst2 + c)* s1 + d
172 yp = three*epst2 + c
173 IF(yp/=zero)s1 = s1 - y/yp
174 epst2 = s1 * s1
175 y = (epst2 + c)* s1 + d
176 yp = three*epst2 + c
177 IF(yp/=zero)s1 = s1 - y/yp
178 epst2 = s1 * s1
179 y = (epst2 + c)* s1 + d
180 yp = three*epst2 + c
181 IF(yp/=zero)s1 = s1 - y/yp
182 epst2 = s1 * s1
183 y = (epst2 + c)* s1 + d
184 yp = three*epst2 + c
185 IF(yp/=zero)s1 = s1 - y/yp
186 ENDIF
187 b = s1
188 c = c + s1**2
189 s2 = half*(-s1 + sqrt(max(zero,(s1**2 - four*c))))
190 ss = half*(-s1 - sqrt(max(zero,(s1**2 - four*c))))
191 s3 = ss
192 IF(ss>s2)THEN
193 s3 = s2
194 s2 = ss
195 ENDIF
196 a=one
197 IF(s1/=zero.AND.s3==zero) a = s2/s1
198 IF(s1==zero.AND.s3/=zero) a = s2/s3
199 IF(s1/=zero.AND.s3/=zero)
200 . a = max(s2/s3, s2/s1)
201 w2 = max(em20,(two - a))**be
202 w1 = one - p/pc
203 w1 = (max(em20,one/w1))**al
204 uvar(i,1) = uvar(i,1) + w1*w2*dpla(i)
205 IF(uvar(i,1)>=dc) THEN
206 off(i) = four_over_5
207 nindx=nindx+1
208 indx(nindx)=i
209 tdele(i) = time
210 ENDIF
211 ENDIF
212 ENDDO
213 IF(nindx>0.AND.imconv==1)THEN
214 DO j=1,nindx
215 i=indx(j)
216#include "lockon.inc"
217 WRITE(istdo,1000)ngl(i)
218 WRITE(iout,1100)ngl(i),time
219#include "lockoff.inc"
220 END DO
221 END IF
222C end Tuler Butcher
223 END IF
224
225Cc deviatoric will be vanished
226 IF(idev==1)THEN
227 nindx=0
228 nindex = 0
229 DO i=1,nel
230 IF(iflag==2.AND.off(i)==one)THEN
231 IF(uvar(i,1)<dc)THEN
232 p = third*(signxx(i) + signyy(i) + signzz(i))
233 e1 = signxx(i) - p
234 e2 = signyy(i) - p
235 e3 = signzz(i) - p
236 e4 = signxy(i)
237 e5 = signyz(i)
238 e6 = signzx(i)
239C -y = (e1-x)(e2-x)(e3-x)
240C - e5^2(e1-x) - e6^2(e2-x) - e4^2(e3-x)
241C + 2e4 e5 e6
242C e1 + e2 + e3 = 0 => terme en x^2 = 0
243C y = x^3 + c x + d
244c yp= 3 x^2 + c
245 e42 = e4*e4
246 e52 = e5*e5
247 e62 = e6*e6
248 c = - half * (e1*e1 + e2*e2 + e3*e3) - e42 - e52 - e62
249 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42
250 & - two*e4*e5*e6
251 cc1 = c*third
252 s1 = sqrt(-cc1)
253 epst2 = s1 * s1
254 y = (epst2 + c)* s1 + d
255 IF(abs(y)>em8)THEN
256 s1 = 1.75 * s1
257 epst2 = s1 * s1
258 y = (epst2 + c)* s1 + d
259 yp = three*epst2 + c
260 IF(yp/=zero)s1 = s1 - y/yp
261 epst2 = s1 * s1
262 y = (epst2 + c)* s1 + d
263 yp = three*epst2 + c
264 IF(yp/=zero)s1 = s1 - y/yp
265 epst2 = s1 * s1
266 y = (epst2 + c)* s1 + d
267 yp = three*epst2 + c
268 IF(yp/=zero)s1 = s1 - y/yp
269 epst2 = s1 * s1
270 y = (epst2 + c)* s1 + d
271 yp = three*epst2 + c
272 IF(yp/=zero)s1 = s1 - y/yp
273 ENDIF
274 b = s1
275 c = c + s1**2
276 s2 = half*(-s1 + sqrt(max(zero,(s1**2 - four*c))))
277 ss = half*(-s1 - sqrt(max(zero,(s1**2 - four*c))))
278 s3 = ss
279 IF(ss>s2)THEN
280 s3 = s2
281 s2 = ss
282 ENDIF
283 a = one
284 IF(s1/=zero.AND.s3==zero) a = s2/s1
285 IF(s1==zero.AND.s3/=zero) a = s2/s3
286 IF(s1/=zero.AND.s3/=zero)
287 . a = max(s2/s3, s2/s1)
288 w2 = max(em20,(two - a))**be
289 w1 = one - one*p/pc
290 w1 = (max(em20,one/w1))**al
291 uvar(i,1) = uvar(i,1) + w1*w2*dpla(i)
292 IF(uvar(i,1)>=dc) THEN
293 nindx=nindx+1
294 indx(nindx)=i
295 signxx(i) = p
296 signyy(i) = p
297 signzz(i) = p
298 signxy(i) = zero
299 signyz(i) = zero
300 signzx(i) = zero
301 ENDIF
302C uvar> DC
303 ELSE
304 p = third*(signxx(i) + signyy(i) + signzz(i))
305 signxx(i) = p
306 signyy(i) = p
307 signzz(i) = p
308 signxy(i) = zero
309 signyz(i) = zero
310 signzx(i) = zero
311 ENDIF
312 ENDIF
313 ENDDO
314 IF(nindx>0.AND.imconv==1)THEN
315 DO j=1,nindx
316 i = indx(j)
317#include "lockon.inc"
318 WRITE(iout, 2000) ngl(i)
319 WRITE(istdo,2100) ngl(i),time
320#include "lockoff.inc"
321 END DO
322 END IF
323 ENDIF
324C-------------Maximum Damage storing for output : 0 < DFMAX < 1--------------
325 DO i=1,nel
326 dfmax(i)= min(one,max(dfmax(i),uvar(i,1)/dc))
327 ENDDO
328C-----------------------------------------------
329 1000 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER ',i10)
330 1100 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER ',i10,
331 . ' AT TIME :',1pe20.13)
332CC
333 2000 FORMAT(1x,' DEVIATORIC STRESS WILL BE VANISHED',i10)
334 2100 FORMAT(1x,' DEVIATORIC STRESS WILL BE VANISHED',i10,
335 . ' AT TIME :',1pe20.13)
336 RETURN
337 END
fail_wilkins_s
subroutine fail_wilkins_s(nel, nuparam, nuvar, time, timestep, uparam, ngl, signxx, signyy, signzz, signxy, signyz, signzx, dpla, uvar, off, dfmax, tdele)
Definition fail_wilkins_s.F:37
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21