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

Go to the source code of this file.

Functions/Subroutines

subroutine fail_orthbiquad_s (nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, ngl, dpla, epsp, uvar, off, signxx, signyy, signzz, signxy, signyz, signzx, dfmax, tdele, aldt)

Function/Subroutine Documentation

◆ fail_orthbiquad_s()

subroutine fail_orthbiquad_s ( integer nel,
integer nuparam,
integer nuvar,
integer mfunc,
integer, dimension(mfunc) kfunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
integer, dimension(nel) ngl,
dpla,
epsp,
uvar,
off,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
dfmax,
tdele,
aldt )

Definition at line 34 of file fail_orthbiquad_s.F.

40C!-----------------------------------------------
41C! I m p l i c i t T y p e s
42C!-----------------------------------------------
43#include "implicit_f.inc"
44C!---------+---------+---+---+--------------------------------------------
45C! VAR | SIZE |TYP| RW| DEFINITION
46C!---------+---------+---+---+--------------------------------------------
47C! NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
48C! NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
49C! NUVAR | 1 | I | R | NUMBER OF FAILURE ELEMENT VARIABLES
50C!---------+---------+---+---+--------------------------------------------
51C! MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
52C! KFUNC | NFUNC | I | R | FUNCTION INDEX not used
53C! NPF | * | I | R | FUNCTION ARRAY
54C! TF | * | F | R | FUNCTION ARRAY
55C!---------+---------+---+---+--------------------------------------------
56C! TIME | 1 | F | R | CURRENT TIME
57C! TIMESTEP| 1 | F | R | CURRENT TIME STEP
58C! UPARAM | NUPARAM | F | R | USER FAILURE PARAMETER ARRAY
59C!---------+---------+---+---+--------------------------------------------
60C! SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
61C! SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
62C! ... | | | |
63C! ... | | | |
64C!---------+---------+---+---+--------------------------------------------
65C! UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
66C! OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
67C!---------+---------+---+---+--------------------------------------------
68#include "mvsiz_p.inc"
69#include "scr17_c.inc"
70#include "units_c.inc"
71#include "comlock.inc"
72C!-----------------------------------------------
73 INTEGER NEL, NUPARAM, NUVAR,NGL(NEL)
74 my_real time,timestep,uparam(*),
75 . signxx(nel),signyy(nel),signzz(nel),
76 . signxy(nel),signyz(nel),signzx(nel),uvar(nel,nuvar),
77 . dpla(nel),epsp(nel),off(nel),dfmax(nel),tdele(nel),
78 . aldt(nel)
79C!-----------------------------------------------
80C! VARIABLES FOR FUNCTION INTERPOLATION
81C!-----------------------------------------------
82 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
83 my_real finter ,tf(*)
84 EXTERNAL finter
85C! Y = FINTER(IFUNC(J),X,NPF,TF,DF)
86C! Y : y = f(x)
87C! X : x
88C! DF : f'(x) = dy/dx
89C! IFUNC(J): FUNCTION INDEX
90C! J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
91C! NPF,TF : FUNCTION PARAMETER
92C!-----------------------------------------------
93C! L o c a l V a r i a b l e s
94C!-----------------------------------------------
95 INTEGER I,K,J,INDX(MVSIZ),NINDX,SEL,NANGLE,IREG,IRATE
96 my_real
97 . df,p,triaxs,svm,sxx,syy,szz,eps_fail,
98 . p1x,p1y,s1x,s1y,s2y,a1,b1,c1,ref_el_len,lambda,fac,
99 . x_1(3),x_2(3),cos2(10,10),epsd0,cjc,rate_scale,frate,
100 . mohr_radius,cos2theta(nel)
101 my_real sig1(nel),sig2(nel),mohr_center
102 my_real, DIMENSION(:), ALLOCATABLE :: q_x11,q_x12,q_x13,q_x21,q_x22,q_x23,q_inst
103C
104 DATA cos2/
105 1 1. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
106 2 0. ,1. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
107 3 -1. ,0. ,2. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
108 4 0. ,-3. ,0. ,4. ,0. ,0. ,0. ,0. ,0. ,0. ,
109 5 1. ,0. ,-8. ,0. ,8. ,0. ,0. ,0. ,0. ,0. ,
110 6 0. ,5. ,0. ,-20. ,0. ,16. ,0. ,0. ,0. ,0. ,
111 7 -1. ,0. ,18. ,0. ,-48. ,0. ,32. ,0. ,0. ,0. ,
112 8 0. ,-7. ,0. ,56. ,0. ,-112.,0. ,64. ,0. ,0. ,
113 9 1. ,0. ,-32. ,0. ,160. ,0. ,-256.,0. ,128. ,0. ,
114 a 0. ,9. ,0. ,-120.,0. ,432. ,0. ,-576 ,0. ,256. /
115C!--------------------------------------------------------------
116 !=======================================================================
117 ! - INITIALISATION OF COMPUTATION ON TIME STEP
118 !=======================================================================
119 ! Recovering model parameters
120 sel = int(uparam(2))
121 ref_el_len = uparam(3)
122 epsd0 = uparam(4)
123 cjc = uparam(5)
124 rate_scale = uparam(6)
125 nangle = int(uparam(7))
126c
127 ! -> Allocation of factors
128 ALLOCATE(q_x11(nangle),q_x12(nangle),q_x13(nangle),
129 . q_x21(nangle),q_x22(nangle),q_x23(nangle))
130 ! Recovering factors
131 IF (sel == 3) THEN
132 ALLOCATE(q_inst(nangle))
133 DO j = 1,nangle
134 q_x11(j) = uparam(8 + 7*(j-1))
135 q_x12(j) = uparam(9 + 7*(j-1))
136 q_x13(j) = uparam(10 + 7*(j-1))
137 q_x21(j) = uparam(11 + 7*(j-1))
138 q_x22(j) = uparam(12 + 7*(j-1))
139 q_x23(j) = uparam(13 + 7*(j-1))
140 q_inst(j) = uparam(14 + 7*(j-1))
141 ENDDO
142 ELSE
143 DO j = 1,nangle
144 q_x11(j) = uparam(8 + 6*(j-1))
145 q_x12(j) = uparam(9 + 6*(j-1))
146 q_x13(j) = uparam(10 + 6*(j-1))
147 q_x21(j) = uparam(11 + 6*(j-1))
148 q_x22(j) = uparam(12 + 6*(j-1))
149 q_x23(j) = uparam(13 + 6*(j-1))
150 ENDDO
151 ENDIF
152 IF (sel == 3) sel = 2
153c
154 ! Recovering functions
155 irate = 0
156 IF (rate_scale /= zero) THEN
157 irate = 1
158 ENDIF
159 ireg = 0
160 IF (ref_el_len /= zero) THEN
161 ireg = irate + 1
162 ENDIF
163c
164 ! At initial time, compute the element size regularization factor
165 IF (uvar(1,3) == zero .AND. ireg > 0) THEN
166 DO i=1,nel
167 lambda = aldt(i)/ref_el_len
168 fac = finter(kfunc(ireg),lambda,npf,tf,df)
169 uvar(i,3) = fac
170 ENDDO
171 ENDIF
172c
173 ! Checking element failure and recovering user variable
174 DO i=1,nel
175 IF (off(i) < em01) off(i) = zero
176 IF (off(i) < one .AND. off(i) > zero) off(i) = off(i)*four_over_5
177 END DO
178C
179 ! Initialization of variable
180 nindx = 0
181c
182 !====================================================================
183 ! - LOOP OVER THE ELEMENT TO COMPUTE THE DAMAGE VARIABLE
184 !====================================================================
185 DO i=1,nel
186c
187 ! Failure strain initialization
188 eps_fail = zero
189c
190 ! If the element is not broken
191 IF (off(i) == one .AND. dpla(i) /= zero) THEN
192c
193 ! Computation of loading angle
194 mohr_radius = sqrt(((signxx(i)-signyy(i))/two)**2 + signxy(i)**2)
195 mohr_center = (signxx(i)+signyy(i))/two
196 IF (mohr_radius > em20) THEN
197 cos2theta(i) = ((signxx(i)-signyy(i))/two)/mohr_radius
198 ELSE
199 cos2theta(i) = one
200 ENDIF
201 sig1(i) = mohr_center + mohr_radius
202 sig2(i) = mohr_center - mohr_radius
203 IF (sig1(i)<zero.OR.((sig2(i)<zero).AND.(sig2(i)<-sig1(i)))) THEN
204 cos2theta(i) = -cos2theta(i)
205 ENDIF
206c
207 ! Computation of the BIQUAD parameters
208 x_1(1:3) = zero
209 x_2(1:3) = zero
210 DO j = 1,nangle
211 DO k = 1,j
212 x_1(1) = x_1(1) + q_x11(j)*cos2(k,j)*(cos2theta(i))**(k-1)
213 x_1(2) = x_1(2) + q_x12(j)*cos2(k,j)*(cos2theta(i))**(k-1)
214 x_1(3) = x_1(3) + q_x13(j)*cos2(k,j)*(cos2theta(i))**(k-1)
215 x_2(1) = x_2(1) + q_x21(j)*cos2(k,j)*(cos2theta(i))**(k-1)
216 x_2(2) = x_2(2) + q_x22(j)*cos2(k,j)*(cos2theta(i))**(k-1)
217 x_2(3) = x_2(3) + q_x23(j)*cos2(k,j)*(cos2theta(i))**(k-1)
218 ENDDO
219 ENDDO
220c
221 ! Computation of hydrostatic stress, Von Mises stress, and stress triaxiality
222 p = third*(signxx(i) + signyy(i) + signzz(i))
223 sxx = signxx(i) - p
224 syy = signyy(i) - p
225 szz = signzz(i) - p
226 svm = half*(sxx**2 + syy**2 + szz**2)
227 . + signxy(i)**2 + signzx(i)**2 + signyz(i)**2
228 svm = sqrt(three*svm)
229 triaxs = p/max(em20,svm)
230 IF (triaxs < -two_third) triaxs = -two_third
231 IF (triaxs > two_third) triaxs = two_third
232c
233 ! Computing the plastic strain at failure according to stress triaxiality
234 IF (triaxs <= third) THEN
235 eps_fail = x_1(1) + x_1(2)*triaxs + x_1(3)*triaxs**2
236 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
237 ELSE
238 SELECT CASE (sel)
239 CASE(1)
240 eps_fail = x_2(1) + x_2(2)*triaxs + x_2(3)*triaxs**2
241 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
242 CASE(2)
243 IF (triaxs <= one/sqr3) THEN ! triax < 0.57735
244 p1x = third
245 p1y = x_1(1) + x_1(2)*p1x + x_1(3)*p1x**2
246 s1x = one/sqr3
247 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
248 a1 = (p1y - s1y)/(p1x - s1x)**2
249 b1 = -two*a1*s1x
250 c1 = a1*s1x**2 + s1y
251 eps_fail = c1 + b1*triaxs + a1*triaxs**2
252 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
253 ELSE ! triax > 0.57735
254 p1x = two*third
255 p1y = x_2(1) + x_2(2)*p1x + x_2(3)*p1x**2
256 s1x = one/sqr3
257 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
258 a1 = (p1y - s1y)/(p1x - s1x)**2
259 b1 = -two*a1*s1x
260 c1 = a1*s1x**2 + s1y
261 eps_fail = c1 + b1*triaxs + a1*triaxs**2
262 IF (ireg > 0) eps_fail = eps_fail*uvar(i,3)
263 ENDIF
264 END SELECT
265 ENDIF
266c
267 ! Strain-rate effects
268 IF (cjc /= zero .OR. irate /= 0) THEN
269 IF (cjc /= zero) THEN
270 frate = one
271 IF (epsp(i) > epsd0) frate = frate + cjc*log(epsp(i)/epsd0)
272 ELSEIF (irate /= 0) THEN
273 frate = finter(kfunc(irate),epsp(i)/rate_scale,npf,tf,df)
274 ENDIF
275 eps_fail = eps_fail*frate
276 ENDIF
277c
278 ! Computation of damage variables
279 dfmax(i) = dfmax(i) + dpla(i)/max(eps_fail,em6)
280 dfmax(i) = min(one,dfmax(i))
281c
282 ! Checking element failure using global damage
283 IF (dfmax(i) >= one .AND. off(i) == one) THEN
284 off(i) = four_over_5
285 nindx = nindx + 1
286 indx(nindx) = i
287 tdele(i) = time
288 ENDIF
289 ENDIF
290 ENDDO
291c------------------------
292c------------------------
293 ! Deallocation of table
294 IF (ALLOCATED(q_x11)) DEALLOCATE(q_x11)
295 IF (ALLOCATED(q_x12)) DEALLOCATE(q_x12)
296 IF (ALLOCATED(q_x13)) DEALLOCATE(q_x13)
297 IF (ALLOCATED(q_x21)) DEALLOCATE(q_x21)
298 IF (ALLOCATED(q_x22)) DEALLOCATE(q_x22)
299 IF (ALLOCATED(q_x23)) DEALLOCATE(q_x23)
300 IF (ALLOCATED(q_inst)) DEALLOCATE(q_inst)
301c------------------------
302c------------------------
303 IF (nindx > 0) THEN
304 DO j=1,nindx
305 i = indx(j)
306#include "lockon.inc"
307 WRITE(iout, 1000) ngl(i),time
308 WRITE(istdo,1100) ngl(i),time
309#include "lockoff.inc"
310 END DO
311 END IF
312c------------------------
313 1000 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER (ORTHBIQUAD) el#',i10,
314 . ' AT TIME :',1pe12.4)
315 1100 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER (ORTHBIQUAD) el#',i10,
316 . ' AT TIME :',1pe12.4)
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