OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_orthbiquad_c.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ fail_orthbiquad_c()

subroutine fail_orthbiquad_c ( integer nel,
integer nuvar,
time,
uparam,
integer, dimension(nel) ngl,
integer ipt,
integer nptot,
signxx,
signyy,
signxy,
signyz,
signzx,
dpla,
epsp,
uvar,
uel1,
off,
offl,
dfmax,
tdel,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
aldt,
integer, dimension(nel), intent(inout) foff,
integer ipg )

Definition at line 31 of file fail_orthbiquad_c.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
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