OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_orthbiquad_c.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_orthbiquad_c ../engine/source/materials/fail/orthbiquad/fail_orthbiquad_c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!|| usermat_shell ../engine/source/materials/mat_share/usermat_shell.F
28!||--- calls -----------------------------------------------------
29!|| finter ../engine/source/tools/curve/finter.F
30!||====================================================================
31 SUBROUTINE fail_orthbiquad_c(
32 1 NEL ,NUVAR ,
33 2 TIME ,UPARAM ,NGL ,IPT ,NPTOT ,
34 3 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
35 4 DPLA ,EPSP ,UVAR ,UEL1 ,
36 5 OFF ,OFFL ,DFMAX ,TDEL ,NFUNC ,
37 6 IFUNC ,NPF ,TF ,ALDT ,FOFF ,
38 7 IPG )
39C-----------------------------------------------
40C ORTHOTROPIC BIQUADRATIC FAILURE MODEL
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C G l o b a l P a r a m e t e r s
47C-----------------------------------------------
48#include "units_c.inc"
49#include "comlock.inc"
50#include "mvsiz_p.inc"
51C---------+---------+---+---+--------------------------------------------
52C VAR | SIZE |TYP| RW| DEFINITION
53C---------+---------+---+---+--------------------------------------------
54C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
55C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
56C---------+---------+---+---+--------------------------------------------
57C NPT | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
58C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
59C NGL | NEL | I | R | ELEMENT NUMBER
60C---------+---------+---+---+--------------------------------------------
61C TIME | 1 | F | R | CURRENT TIME
62C---------+---------+---+---+--------------------------------------------
63C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
64C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
65C ... | | | |
66C---------+---------+---+---+--------------------------------------------
67C PLA | NEL | F |R/W| PLASTIC STRAIN
68C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
69C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
70C FOFF | NEL | I |R/W| DELETED INTEGRATION POINT FLAG (=1 ON, =0 OFF)
71C---------+---------+---+---+--------------------------------------------
72C IPT CURRENT INTEGRATION POINT IN ALL LAYERS
73C IPTT CURRENT INTEGRATION POINT IN THE LAYER (FOR OUTPUT ONLY)
74C NPTOT NUMBER OF INTEGRATION POINTS IN ALL LAYERS (TOTAL)
75C-----------------------------------------------
76C I N P U T A r g u m e n t s
77C-----------------------------------------------
78 INTEGER NEL,NUVAR,IPT,NPTOT,NFUNC,IPG
79 INTEGER NGL(NEL),IFUNC(NFUNC)
80 my_real TIME,UPARAM(*),DPLA(NEL),EPSP(NEL),
81 . UEL1(NEL),ALDT(NEL)
82C-----------------------------------------------
83C I N P U T O U T P U T A r g u m e n t s
84C-----------------------------------------------
85 INTEGER, DIMENSION(NEL), INTENT(INOUT) :: FOFF
86 my_real
87 . uvar(nel,nuvar),off(nel),offl(nel),
88 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
89 . dfmax(nel),tdel(nel)
90C-----------------------------------------------
91C VARIABLES FOR FUNCTION INTERPOLATION
92C-----------------------------------------------
93 INTEGER NPF(*)
94 my_real FINTER ,TF(*)
95 EXTERNAL FINTER
96C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
97C Y : y = f(x)
98C X : x
99C DYDX : f'(x) = dy/dx
100C IFUNC(J): FUNCTION INDEX
101C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
102C NPF,TF : FUNCTION PARAMETER
103C-----------------------------------------------
104C L o c a l V a r i a b l e s
105C-----------------------------------------------
106 my_real hydros,vmises,triaxs,ref_el_len
107 INTEGER I,J,K,L,NINDX1,NINDX2,SEL,NANGLE
108 INTEGER FAIL_COUNT,IT,IREG,IRATE
109 INTEGER, DIMENSION(MVSIZ) :: INDX1,INDX2
110 my_real X_1(3),X_2(3),c3,DYDX,COS2(10,10),EPSD0,CJC,RATE_SCALE,FRATE
111 my_real EPS_FAIL,EPS_FAIL2,DAMAGE,INST,FSCALE_LEN,MOHR_RADIUS
112 my_real P1X,P1Y,S1X,S1Y,S2Y,A1,A2,B1,B2,C1,C2,FAC,LAMBDA,COS2THETA(NEL)
113 my_real sig1(nel),sig2(nel),mohr_center
114 my_real, DIMENSION(:), ALLOCATABLE :: q_x11,q_x12,q_x13,q_x21,q_x22,q_x23,q_inst
115C
116 DATA cos2/
117 1 1. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
118 2 0. ,1. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
119 3 -1. ,0. ,2. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
120 4 0. ,-3. ,0. ,4. ,0. ,0. ,0. ,0. ,0. ,0. ,
121 5 1. ,0. ,-8. ,0. ,8. ,0. ,0. ,0. ,0. ,0. ,
122 6 0. ,5. ,0. ,-20. ,0. ,16. ,0. ,0. ,0. ,0. ,
123 7 -1. ,0. ,18. ,0. ,-48. ,0. ,32. ,0. ,0. ,0. ,
124 8 0. ,-7. ,0. ,56. ,0. ,-112.,0. ,64. ,0. ,0. ,
125 9 1. ,0. ,-32. ,0. ,160. ,0. ,-256.,0. ,128. ,0. ,
126 a 0. ,9. ,0. ,-120.,0. ,432. ,0. ,-576 ,0. ,256. /
127C-----------------------------------------------
128c! UVAR1 = damage due to instability (triax between 1/3 and 2/3)
129c! UVAR2 = actual integration point ON or OFF
130c! UVAR3 =
131C-----------------------------------------------
132c
133 !=======================================================================
134 ! - INITIALISATION OF COMPUTATION ON TIME STEP
135 !=======================================================================
136 ! Recovering model parameters
137 sel = int(uparam(2))
138 ref_el_len = uparam(3)
139 epsd0 = uparam(4)
140 cjc = uparam(5)
141 rate_scale = uparam(6)
142 nangle = int(uparam(7))
143 frate = -huge(frate)
144 ! -> Allocation of factors
145 ALLOCATE(q_x11(nangle),q_x12(nangle),q_x13(nangle),
146 . q_x21(nangle),q_x22(nangle),q_x23(nangle))
147 ! Recovering factors
148 IF (sel == 3) THEN
149 ALLOCATE(q_inst(nangle))
150 DO j = 1,nangle
151 q_x11(j) = uparam(8 + 7*(j-1))
152 q_x12(j) = uparam(9 + 7*(j-1))
153 q_x13(j) = uparam(10 + 7*(j-1))
154 q_x21(j) = uparam(11 + 7*(j-1))
155 q_x22(j) = uparam(12 + 7*(j-1))
156 q_x23(j) = uparam(13 + 7*(j-1))
157 q_inst(j) = uparam(14 + 7*(j-1))
158 ENDDO
159 ELSE
160 DO j = 1,nangle
161 q_x11(j) = uparam(8 + 6*(j-1))
162 q_x12(j) = uparam(9 + 6*(j-1))
163 q_x13(j) = uparam(10 + 6*(j-1))
164 q_x21(j) = uparam(11 + 6*(j-1))
165 q_x22(j) = uparam(12 + 6*(j-1))
166 q_x23(j) = uparam(13 + 6*(j-1))
167 ENDDO
168 ENDIF
169c
170 ! Recovering functions
171 irate = 0
172 IF (rate_scale /= zero) THEN
173 irate = 1
174 ENDIF
175 ireg = 0
176 IF (ref_el_len /= zero) THEN
177 ireg = irate + 1
178 ENDIF
179c
180 ! at initial time, compute the element size regularization factor
181 IF (uvar(1,3) == zero .AND. ireg > 0) THEN
182 DO i=1,nel
183 lambda = aldt(i)/ref_el_len
184 fac = finter(ifunc(ireg),lambda,npf,tf,dydx)
185 uvar(i,3) = fac
186 ENDDO
187 ENDIF
188C
189 ! Initialization of variable
190 nindx1 = 0
191 nindx2 = 0
192c
193 !====================================================================
194 ! - LOOP OVER THE ELEMENT TO COMPUTE THE DAMAGE VARIABLE
195 !====================================================================
196 DO i = 1,nel
197c
198 ! Failure strain initialization
199 eps_fail = zero
200 eps_fail2 = zero
201c
202 ! If the element is not broken
203 IF (off(i) == one .AND. dpla(i) /= zero .AND. foff(i) == 1) THEN
204c
205 ! Computation of loading angle
206 mohr_radius = sqrt(((signxx(i)-signyy(i))/two)**2 + signxy(i)**2)
207 mohr_center = (signxx(i)+signyy(i))/two
208 IF (mohr_radius > em20) THEN
209 cos2theta(i) = ((signxx(i)-signyy(i))/two)/mohr_radius
210 ELSE
211 cos2theta(i) = one
212 ENDIF
213 sig1(i) = mohr_center + mohr_radius
214 sig2(i) = mohr_center - mohr_radius
215 IF (sig1(i)<zero.OR.((sig2(i)<zero).AND.(sig2(i)<-sig1(i)))) THEN
216 cos2theta(i) = -cos2theta(i)
217 ENDIF
218c
219 ! Computation of the BIQUAD parameters
220 x_1(1:3) = zero
221 x_2(1:3) = zero
222 inst = zero
223 DO j = 1,nangle
224 DO k = 1,j
225 x_1(1) = x_1(1) + q_x11(j)*cos2(k,j)*(cos2theta(i))**(k-1)
226 x_1(2) = x_1(2) + q_x12(j)*cos2(k,j)*(cos2theta(i))**(k-1)
227 x_1(3) = x_1(3) + q_x13(j)*cos2(k,j)*(cos2theta(i))**(k-1)
228 x_2(1) = x_2(1) + q_x21(j)*cos2(k,j)*(cos2theta(i))**(k-1)
229 x_2(2) = x_2(2) + q_x22(j)*cos2(k,j)*(cos2theta(i))**(k-1)
230 x_2(3) = x_2(3) + q_x23(j)*cos2(k,j)*(cos2theta(i))**(k-1)
231 IF (sel == 3) inst = inst + q_inst(j)*cos2(k,j)*(cos2theta(i))**(k-1)
232 ENDDO
233 ENDDO
234c
235 ! Computation of hydrostatic stress, Von Mises stress, and stress triaxiality
236 hydros = (signxx(i)+ signyy(i))/three
237 vmises = sqrt((signxx(i)**2)+(signyy(i)**2)-(signxx(i)*signyy(i))+(three*signxy(i)**2))
238 triaxs = hydros / max(em20,vmises)
239c
240 ! Computing the plastic strain at failure according to stress triaxiality
241 IF (triaxs <= third) THEN ! triax < 1/3
242 eps_fail = x_1(1) + x_1(2)*triaxs + x_1(3)*triaxs**2
243 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
244 ELSE ! triax > 1/3
245 SELECT CASE (sel)
246 CASE(1)
247 eps_fail = x_2(1) + x_2(2)*triaxs + x_2(3)*triaxs**2
248 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
249 CASE(2)
250 IF (triaxs <= one/sqr3) THEN ! triax < 0.57735
251 p1x = third
252 p1y = x_1(1) + x_1(2)*p1x + x_1(3)*p1x**2
253 s1x = one/sqr3
254 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
255 a1 = (p1y - s1y)/(p1x - s1x)**2
256 b1 = -two*a1*s1x
257 c1 = a1*s1x**2 + s1y
258 eps_fail = c1 + b1*triaxs + a1*triaxs**2
259 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
260 ELSE ! triax > 0.57735
261 p1x = two*third
262 p1y = x_2(1) + x_2(2)*p1x + x_2(3)*p1x**2
263 s1x = one/sqr3
264 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
265 a1 = (p1y - s1y)/(p1x - s1x)**2
266 b1 = -two*a1*s1x
267 c1 = a1*s1x**2 + s1y
268 eps_fail = c1 + b1*triaxs + a1*triaxs**2
269 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
270 ENDIF
271 CASE(3)
272 IF (triaxs <= one/sqr3) THEN ! triax < 0.57735
273 p1x = third
274 p1y = x_1(1) + x_1(2)*p1x + x_1(3)*p1x**2
275 s1x = one/sqr3
276 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
277 s2y = inst
278 a1 = (p1y - s1y)/(p1x - s1x)**2
279 a2 = (p1y - s2y)/(p1x - s1x)**2
280 b1 = -two*a1*s1x
281 b2 = -two*a2*s1x
282 c1 = a1*s1x**2 + s1y
283 c2 = a2*s1x**2 + s2y
284 eps_fail = c1 + b1*triaxs + a1*triaxs**2
285 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
286 eps_fail2 = c2 + b2*triaxs + a2*triaxs**2 ! INSTABILITY
287 IF (ireg > 0) eps_fail2 = eps_fail2*uvar(i,3)
288 ELSE ! triax > 0.57735
289 p1x = two*third
290 p1y = x_2(1) + x_2(2)*p1x + x_2(3)*p1x**2
291 s1x = one/sqr3
292 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
293 s2y = inst
294 a1 = (p1y - s1y)/(p1x - s1x)**2
295 a2 = (p1y - s2y)/(p1x - s1x)**2
296 b1 = -two*a1*s1x
297 b2 = -two*a2*s1x
298 c1 = a1*s1x**2 + s1y
299 c2 = a2*s1x**2 + s2y
300 eps_fail = c1 + b1*triaxs + a1*triaxs**2
301 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
302 eps_fail2 = c2 + b2*triaxs + a2*triaxs**2 ! INSTABILITY
303 IF (ireg > 0) eps_fail2 = eps_fail2*uvar(i,3)
304 ENDIF
305 END SELECT
306 ENDIF
307c
308 ! Strain-rate effects
309 IF (cjc /= zero .OR. irate /= 0) THEN
310 IF (cjc /= zero) THEN
311 frate = one
312 IF (epsp(i) > epsd0) frate = frate + cjc*log(epsp(i)/epsd0)
313 ELSEIF (irate /= 0) THEN
314 frate = finter(ifunc(irate),epsp(i)/rate_scale,npf,tf,dydx)
315 ENDIF
316 eps_fail = eps_fail*frate
317 IF (sel == 3) eps_fail2 = eps_fail2*frate
318 ENDIF
319c
320 ! Computation of damage variables
321 dfmax(i) = dfmax(i) + dpla(i)/max(eps_fail,em6)
322 dfmax(i) = min(one,dfmax(i))
323 IF (sel == 3) THEN
324 damage = uvar(i,1)
325 IF (eps_fail2 > zero) THEN
326 damage = damage + dpla(i)/max(eps_fail2,em6)
327 uvar(i,1) = damage
328 ENDIF
329 ELSE
330 damage = zero
331 ENDIF
332c
333 ! Checking element failure using global damage
334 IF (offl(i) == one .AND. dfmax(i) >= one) THEN
335 foff(i) = 0
336 nindx1 = nindx1 + 1
337 indx1(nindx1) = i
338 ENDIF
339c
340 ! Checking element failure using instability damage
341 IF (damage > one .AND. uvar(i,2) == zero .AND. offl(i) == one) THEN
342 uvar(i,2) = one
343 uel1(i) = uel1(i) + one
344 IF (int(uel1(i)) == nptot) THEN
345 nindx2 = nindx2 + 1
346 indx2(nindx2) = i
347 tdel(i) = time
348 off(i) = zero
349 signxx(i) = zero
350 signyy(i) = zero
351 signxy(i) = zero
352 signyz(i) = zero
353 signzx(i) = zero
354 ENDIF
355 ENDIF
356 ENDIF
357 ENDDO
358c------------------------
359c------------------------
360 ! Deallocation of table
361 IF (ALLOCATED(q_x11)) DEALLOCATE(q_x11)
362 IF (ALLOCATED(q_x12)) DEALLOCATE(q_x12)
363 IF (ALLOCATED(q_x13)) DEALLOCATE(q_x13)
364 IF (ALLOCATED(q_x21)) DEALLOCATE(q_x21)
365 IF (ALLOCATED(q_x22)) DEALLOCATE(q_x22)
366 IF (ALLOCATED(q_x23)) DEALLOCATE(q_x23)
367 IF (ALLOCATED(q_inst)) DEALLOCATE(q_inst)
368c------------------------
369c------------------------
370 IF (nindx1 > 0) THEN
371 DO j=1,nindx1
372 i = indx1(j)
373#include "lockon.inc"
374 WRITE(iout, 2000) ngl(i),ipg,ipt,time
375 WRITE(istdo,2000) ngl(i),ipg,ipt,time
376#include "lockoff.inc"
377 END DO
378 END IF
379 IF (nindx2 > 0) THEN
380 DO j=1,nindx2
381 i = indx2(j)
382#include "lockon.inc"
383 WRITE(iout, 3000) ngl(i)
384 WRITE(iout, 2200) ngl(i), time
385 WRITE(istdo,3000) ngl(i)
386 WRITE(istdo,2200) ngl(i), time
387#include "lockoff.inc"
388 END DO
389 END IF
390c------------------------
391 2000 FORMAT(1x,'for shell element(orthbiquad)',I10,1X,'gauss point',I3,
392 . 1X,'layer',I3,':',/,
393 . 1X,'stress tensor set to zero',1X,'at time :',1PE12.4)
394 2200 FORMAT(1X,' *** rupture of shell element (orthbiquad)',I10,1X,
395 . ' at time :',1PE12.4)
396 3000 FORMAT(1X,'for shell element(orthbiquad)',I10,
397 . 1X,'instability reached.')
398c------------------------
399 RETURN
400 END
the
end diagonal values have been computed in the(sparse) matrix id.SOL
fail_orthbiquad_c
subroutine fail_orthbiquad_c(nel, nuvar, time, uparam, ngl, ipt, nptot, signxx, signyy, signxy, signyz, signzx, dpla, epsp, uvar, uel1, off, offl, dfmax, tdel, nfunc, ifunc, npf, tf, aldt, foff, ipg)
Definition fail_orthbiquad_c.F:39
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