OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
r23l113def3.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!|| r23l113def3 ../engine/source/elements/spring/r23l113def3.F
25!||--- called by ------------------------------------------------------
26!|| r23law113 ../engine/source/elements/spring/r23law113.F
27!||--- calls -----------------------------------------------------
28!|| redef3 ../engine/source/elements/spring/redef3.F90
29!|| repla3 ../engine/source/elements/spring/repla3.F
30!||--- uses -----------------------------------------------------
31!|| python_funct_mod ../common_source/modules/python_mod.F90
32!|| redef3_mod ../engine/source/elements/spring/redef3.F90
33!||====================================================================
34 SUBROUTINE r23l113def3(PYTHON,
35 1 SKEW, IPM, IGEO, MID,
36 2 PID, GEO, UPARAM, FX,
37 3 FY, FZ, E, DX,
38 4 DY, DZ, NPF, TF,
39 5 OFF, DPX, DPY, DPZ,
40 6 DPX2, DPY2, DPZ2, FXEP,
41 7 FYEP, FZEP, X0, Y0,
42 8 Z0, XMOM, YMOM, ZMOM,
43 9 RX, RY, RZ, RPX,
44 A RPY, RPZ, XMEP, YMEP,
45 B ZMEP, RPX2, RPY2, RPZ2,
46 C ANIM, POSX, POSY, POSZ,
47 D POSXX, POSYY, POSZZ, FR_WAVE,
48 E E6, NEL, EXX2, EYX2,
49 F EZX2, EXY2, EYY2, EZY2,
50 G EXZ2, EYZ2, EZZ2, AL2DP,
51 H NGL, CRIT_NEW,X0_ERR, ALDP,
52 I YIELDX, YIELDY, YIELDZ, YIELDX2,
53 J YIELDY2, YIELDZ2, EXX, EYX,
54 K EZX, EXY, EYY, EZY,
55 L EXZ, EYZ, EZZ, XCR,
56 M RX1, RY1, RZ1, RX2,
57 N RY2, RZ2, XIN, AK,
58 O XM, XKM, XCM, XKR,
59 P VX1, VX2, VY1, VY2,
60 Q VZ1, VZ2, NUVAR, UVAR,
61 R MASS, DX0, DY0, DZ0,
62 S RX0, RY0, RZ0, NFT,
63 T STF, SANIN, IRESP, SNPC,
64 U SZYIELD_COMP,SZXXOLD_COMP,YIELDXC,YIELDYC,
65 V YIELDZC, YIELDRXC,YIELDRYC,YIELDRZC,
66 W DXOLDC, DYOLDC, DZOLDC, DRXOLDC,
67 W DRYOLDC, DRZOLDC)
68C
69 USE redef3_mod
70 USE python_funct_mod
71C-----------------------------------------------
72C I m p l i c i t T y p e s
73C-----------------------------------------------
74#include "implicit_f.inc"
75#include "comlock.inc"
76C-----------------------------------------------
77C G l o b a l P a r a m e t e r s
78C-----------------------------------------------
79#include "mvsiz_p.inc"
80C-----------------------------------------------
81C C o m m o n B l o c k s
82C-----------------------------------------------
83#include "param_c.inc"
84#include "com04_c.inc"
85#include "com08_c.inc"
86#include "scr14_c.inc"
87#include "scr17_c.inc"
88#include "units_c.inc"
89#include "com01_c.inc"
90#include "impl1_c.inc"
91C-----------------------------------------------
92C D u m m y A r g u m e n t s
93C-----------------------------------------------
94 type(python_), intent(inout) :: PYTHON
95 INTEGER, INTENT(IN) :: NFT
96 INTEGER NPF(*),IGEO(NPROPGI,*),NEL,NGL(*),PID(*),MID(*),NUVAR,
97 . IPM(NPROPMI,*)
98 INTEGER, INTENT(IN) :: STF !< Size of TF
99 INTEGER, INTENT(IN) :: SANIN !< Size of ANIM
100 INTEGER, INTENT(IN) :: IRESP !< Single precision flag
101 INTEGER, INTENT(IN) :: SNPC !< Size of NPF
102 INTEGER, INTENT(IN) :: SZYIELD_COMP
103 INTEGER, INTENT(IN) :: SZXXOLD_COMP
104C REAL
105 my_real
106 . SKEW(LSKEW,*), GEO(NPROPG,*), FX(*), FY(*), FZ(*), E(*), DX(*),
107 . DY(*), DZ(*), TF(*), OFF(*), DPX(*), DPY(*), DPZ(*), FXEP(*),
108 . FYEP(*), FZEP(*), X0(*), Y0(*), Z0(*), XMOM(*), YMOM(*),
109 . ZMOM(*), RX(*), RY(*), RZ(*), RPX(*), RPY(*), RPZ(*), XMEP(*),
110 . YMEP(*), ZMEP(*), DPX2(*), DPY2(*), DPZ2(*),RPX2(*), RPY2(*),
111 . RPZ2(*),ANIM(*),FR_WAVE(*),E6(NEL,6),
112 . EXX2(MVSIZ), EYX2(MVSIZ), EZX2(MVSIZ),
113 . EXY2(MVSIZ), EYY2(MVSIZ), EZY2(MVSIZ),
114 . EXZ2(MVSIZ), EYZ2(MVSIZ), EZZ2(MVSIZ),
115 . CRIT_NEW(*), X0_ERR(MVSIZ),YIELDX(*),YIELDY(*),
116 . YIELDZ(*),YIELDX2(*),YIELDY2(*),YIELDZ2(*),
117 . EXX(MVSIZ), EYX(MVSIZ), EZX(MVSIZ), EXY(MVSIZ),
118 . EYY(MVSIZ), EZY(MVSIZ), EXZ(MVSIZ), EYZ(MVSIZ),
119 . EZZ(MVSIZ), XCR(MVSIZ), RX1(MVSIZ), RX2(MVSIZ),
120 . RY1(MVSIZ), RY2(MVSIZ), RZ1(MVSIZ), RZ2(MVSIZ),
121 . XIN(MVSIZ),AK(MVSIZ),XM(MVSIZ),XKM(MVSIZ),XCM(MVSIZ),
122 . XKR(MVSIZ),VX1(MVSIZ),VX2(MVSIZ),VY1(MVSIZ),VY2(MVSIZ),
123 . VZ1(MVSIZ),VZ2(MVSIZ),UVAR(NUVAR,*),UPARAM(*),MASS(*),
124 . DX0(*),DY0(*),DZ0(*),RX0(*),RY0(*),RZ0(*)
125 my_real, DIMENSION(6,NEL), INTENT(INOUT) :: POSX,POSY,POSZ,POSXX,POSYY,POSZZ
126 my_real, INTENT(INOUT) :: YIELDXC(SZYIELD_COMP),YIELDYC(SZYIELD_COMP),
127 . yieldzc(szyield_comp),yieldrxc(szyield_comp),yieldryc(szyield_comp),
128 . yieldrzc(szyield_comp)
129 my_real, INTENT(INOUT) :: dxoldc(szxxold_comp),dyoldc(szxxold_comp),
130 . dzoldc(szxxold_comp),drxoldc(szxxold_comp),dryoldc(szxxold_comp),
131 . drzoldc(szxxold_comp)
132 TARGET :: uvar
133 DOUBLE PRECISION ALDP(MVSIZ),AL2DP(MVSIZ)
134C-----------------------------------------------
135C L o c a l V a r i a b l e s
136C-----------------------------------------------
137 INTEGER INDX(MVSIZ),
138 . iecrou(mvsiz), ifunc(mvsiz), ifv(mvsiz), ifunc2(mvsiz),
139 . i, ileng, j, kk, ifail(mvsiz),ifail2(mvsiz),
140 . nindx,israte, ifunc3(mvsiz),i1,i2,i3,i4,i5,i6,i7,i8,
141 . i9,i10,i11,i12,i13,i14,i15,if1,if2,if3,if4,iadbuf,nupar,
142 . if5,ifunc4(mvsiz)
143C REAL
144 my_real
145 . xk(mvsiz), yk(mvsiz),
146 . zk(mvsiz),xc(mvsiz), yc(mvsiz), zc(mvsiz),xh(mvsiz),
147 . xhr(mvsiz),dxold(mvsiz), dyold(mvsiz), dzold(mvsiz),
148 . b(mvsiz), d(mvsiz), epla(mvsiz),
149 . dv(mvsiz),vrt(mvsiz),vrr(mvsiz),
150 . dmn(mvsiz),dmx(mvsiz),xl0(mvsiz),crit(mvsiz),
151 . xn(mvsiz),ff(mvsiz),lscale(mvsiz),ee(mvsiz),gf3(mvsiz),
152 . hx(mvsiz), hy(mvsiz), hz(mvsiz)
153 my_real
154 . at,dt05,xka,yka,zka,cc,cn,xa,xb,dlim,vfail,
155 . x21, y21, z21, epxy, epxz,
156 . vx21, vy21, vz21, ryavp, rzavp,eyzp,exzp,
157 . ryav, rzav,den, c, cp, exyp,
158 . x21phi, y21phi, z21phi, vx21phi, vy21phi, vz21phi,
159 . ryav1, rzav1, ryav1p, rzav1p,asrate,not_used(2)
160 my_real :: max_slope(mvsiz)
161 DOUBLE PRECISION X0DP(MVSIZ)
162 my_real ,DIMENSION(:), POINTER :: coord_old
163 LOGICAL :: ANY_PYTHON_FUNCTION
164 TARGET :: NOT_USED
165C-----------------------------------------------
166C
167 not_used = zero
168C
169 nupar = 4
170 i1 = nupar
171 i2 = i1 + 6
172 i3 = i2 + 6
173 i4 = i3 + 6
174 i5 = i4 + 6
175 i6 = i5 + 6
176 i7 = i6 + 6
177 i8 = i7 + 6
178 i9 = i8 + 6
179 i10 = i9 + 6
180 i11 = i10 + 6
181 i12 = i11 + 6
182 i13 = i12 + 6
183 i14 = i13 + 6
184 nupar = nupar + 84
185 i15 = i14 + 6 + 28 ! 4 + 14*6 + 28
186 DO i=1,nel
187 iadbuf= ipm(7,mid(i)) - 1
188 epla(i)=zero
189 xm(i)=mass(i)
190 xk(i)=uparam(iadbuf + i11 + 1)
191 xc(i)=uparam(iadbuf + i12 + 1)
192 yk(i)=uparam(iadbuf + i11 + 2)
193 yc(i)=uparam(iadbuf + i12 + 2)
194 zk(i)=uparam(iadbuf + i11 + 3)
195 zc(i)=uparam(iadbuf + i12 + 3)
196 xka=xk(i)*uparam(iadbuf + i1 + 1)
197 yka=yk(i)*uparam(iadbuf + i1 + 2)
198 zka=zk(i)*uparam(iadbuf + i1 + 3)
199 xkm(i)= max(xka,yka,zka)
200 hx(i) = uparam(iadbuf + i14 + 1)
201 hy(i) = uparam(iadbuf + i14 + 2)
202 hz(i) = uparam(iadbuf + i14 + 3)
203
204 xh(i)= max(hx(i),hy(i),hz(i))
205 xcm(i)= max(xc(i),yc(i),zc(i))
206 xcm(i)= xcm(i)+xh(i)
207
208 xkr(i)= max(yka,zka) * aldp(i) * aldp(i)
209 xcr(i)= (max(yc(i),zc(i)) + max(hy(i),hz(i))) * aldp(i) * aldp(i)
210 vrt(i) = uparam(iadbuf + nupar + 1)
211 vrr(i) = uparam(iadbuf + nupar + 2)
212 ifail(i) = nint(uparam(iadbuf + 1 ))
213 ifail2(i)= nint(uparam(iadbuf + 3 ))
214 ENDDO
215C
216 IF (inispri /= 0 .and. tt == zero) THEN
217 DO i=1,nel
218 xl0(i)= x0(i)
219! if not initialized length
220 IF (xl0(i) == zero) xl0(i) = aldp(i)
221 ENDDO
222 ENDIF
223C
224 IF (tt == zero) THEN
225 DO i=1,nel
226 x0(i)= aldp(i) ! cast double to my_real
227 ENDDO
228 ENDIF
229C
230 IF (scodver >= 101) THEN
231 IF (tt == zero) THEN
232 DO i=1,nel
233 x0_err(i)= aldp(i)-x0(i) ! difference between double and my_real
234 ENDDO
235 ENDIF
236 ENDIF
237C
238 IF ( inispri /= 0 .and. tt == zero) THEN
239 DO i=1,nel
240 x0(i)= xl0(i)
241 ENDDO
242 ENDIF
243C
244 DO i=1,nel
245 x0dp(i)= x0(i) ! cast double to My_real
246 ENDDO
247C
248 IF (scodver >= 101) THEN
249 DO i=1,nel
250 x0dp(i)= x0(i) + x0_err(i) ! difference between double and my_real
251 ENDDO
252 ENDIF
253C---------------------
254C TRANSLATIONS
255C---------------------
256 DO i=1,nel
257 dxold(i)=dx(i)
258 dyold(i)=dy(i)
259 dzold(i)=dz(i)
260 ENDDO
261!
262 IF (inispri /= 0 .and. tt == zero) THEN
263 DO i=1,nel
264 dxold(i)=dx0(i)
265 dyold(i)=dy0(i)
266 dzold(i)=dz0(i)
267 ENDDO
268 ENDIF
269C
270 dt05=half*dt1
271 IF (ismdisp > 0) THEN
272 DO i=1,nel
273 vx21 = vx2(i)-vx1(i)
274 vy21 = vy2(i)-vy1(i)
275 vz21 = vz2(i)-vz1(i)
276 dx(i) = dxold(i)+(vx21*exx(i)+vy21*eyx(i)+vz21*ezx(i))*dt1
277 dy(i) = dyold(i)+(vx21*exy(i)+vy21*eyy(i)+vz21*ezy(i))*dt1
278 dz(i) = dzold(i)+(vx21*exz(i)+vy21*eyz(i)+vz21*ezz(i))*dt1
279C
280 x21 = (rx2(i)+rx1(i))
281 y21 = (ry2(i)+ry1(i))
282 z21 = (rz2(i)+rz1(i))
283C
284 ryav1 = (x21*exy2(i)+y21*eyy2(i)+z21*ezy2(i))
285 rzav1 = (x21*exz2(i)+y21*eyz2(i)+z21*ezz2(i))
286C
287 ryav = dt05 * ryav1
288 rzav = dt05 * rzav1
289C
290 dy(i) = dy(i) - rzav * al2dp(i)
291 dz(i) = dz(i) + ryav * al2dp(i)
292C
293 crit(i) = zero
294 ENDDO
295 ELSE
296 DO i=1,nel
297 vx21 = vx2(i)-vx1(i)
298 vy21 = vy2(i)-vy1(i)
299 vz21 = vz2(i)-vz1(i)
300C
301 epxy = (vx21*exy2(i)+vy21*eyy2(i)+vz21*ezy2(i))*dt05
302 epxz = (vx21*exz2(i)+vy21*eyz2(i)+vz21*ezz2(i))*dt05
303C
304 x21 = (rx2(i)+rx1(i))
305 y21 = (ry2(i)+ry1(i))
306 z21 = (rz2(i)+rz1(i))
307C
308 ryav1 = (x21*exy2(i)+y21*eyy2(i)+z21*ezy2(i))
309 rzav1 = (x21*exz2(i)+y21*eyz2(i)+z21*ezz2(i))
310C
311 at=epxz/max(al2dp(i),em30)
312 at=atan(at)
313 ryav = dt05 * (ryav1) + two * at
314 at=epxy/max(al2dp(i),em30)
315 at=atan(at)
316 rzav = dt05 * (rzav1) - two * at
317C
318 dx(i) = aldp(i) - x0dp(i)
319 dy(i) = dyold(i) - rzav * al2dp(i)
320 dz(i) = dzold(i) + ryav * al2dp(i)
321C
322 crit(i) = zero
323 ENDDO
324 ENDIF !(ISMDISP > 0) THEN
325C
326 DO i=1,nel
327 iadbuf = ipm(7,mid(i)) - 1
328 ileng = nint(uparam(iadbuf + 2))
329 IF (ileng /= 0) THEN
330 xl0(i)=x0dp(i)
331 ELSE
332 xl0(i)=one
333 ENDIF
334 ENDDO
335C-------------------------------
336 nindx = 0
337 if1 = 0
338 if2 = 6
339 if3 = 12
340 if4 = 18
341 if5 = 24
342 DO i=1,nel
343 iadbuf = ipm(7,mid(i)) - 1
344 ifunc(i) = ipm(10 + if1 + 1,mid(i))
345 ifv(i) = ipm(10 + if2 + 1,mid(i))
346 ifunc2(i)= ipm(10 + if3 + 1,mid(i))
347 ifunc3(i)= ipm(10 + if4 + 1,mid(i))
348 ifunc4(i)= ipm(10 + if5 + 1,mid(i))
349 iecrou(i)= nint(uparam(iadbuf + i13 + 1))
350 ak(i) = uparam(iadbuf + i1 + 1)
351 b(i) = uparam(iadbuf + i2 + 1)
352 d(i) = uparam(iadbuf + i3 + 1)
353 ee(i) = uparam(iadbuf + i4 + 1)
354 gf3(i) = uparam(iadbuf + i5 + 1)
355 ff(i) = uparam(iadbuf + i6 + 1)
356 lscale(i)= uparam(iadbuf + i7 + 1)
357 dmn(i) = uparam(iadbuf + i8 + 1)
358 dmx(i) = uparam(iadbuf + i9 + 1)
359 ENDDO
360 IF (nuvar >=1) THEN
361 coord_old => uvar(1,1:nel)
362 ELSE
363 coord_old => not_used
364 ENDIF
365 CALL redef3(python,
366 . fx, xk, dx, fxep,
367 . dxold, dpx, tf, npf,
368 . xc, off, e6(1,1), dpx2,
369 . anim, anim_fe(11),posx,
370 . xl0, dmn, dmx, dv,
371 . ff, lscale, ee, gf3,
372 . ifunc3, yieldx, aldp, ak,
373 . b, d, iecrou, ifunc,
374 . ifv, ifunc2, epla, coord_old,
375 . nel, nft, stf, sanin,
376 . dt1, iresp, impl_s, idyna,
377 . snpc, yieldc=yieldxc, xx_oldc=dxoldc,
378 . fx0=uparam(iadbuf+i15+1),ifunc4=ifunc4,pos6=posx(6,1))
379C
380 DO i=1,nel
381 cc = uparam(iadbuf + nupar + 3)
382 cn = uparam(iadbuf + nupar + 9)
383 xa = uparam(iadbuf + nupar + 15)
384 xb = uparam(iadbuf + nupar + 21)
385 IF (off(i) == one .AND. dmx(i) /= zero .AND. dmn(i)/= zero) THEN
386 IF (ifail2(i) == 0) THEN
387 xa = one
388 xb = two
389 IF (dx(i) > zero) THEN
390 dlim = dx(i) / (xl0(i)*dmx(i))
391 ELSE
392 dlim = dx(i) / (xl0(i)*dmn(i))
393 ENDIF
394 ELSE
395 vfail = cc * (abs(dv(i)/vrt(i)))**cn
396 IF (ifail2(i) == 1) THEN
397 IF (dx(i) > zero) THEN
398 dlim = dx(i) / (xl0(i)*dmx(i) + vfail)
399 ELSE
400 dlim = dx(i) / (xl0(i)*dmn(i) - vfail)
401 ENDIF
402 ELSEIF (ifail2(i) == 2) THEN
403 IF (fx(i) > zero) THEN
404 dlim = fx(i) / (dmx(i) + vfail)
405 ELSE
406 dlim = fx(i) / (dmn(i) - vfail)
407 ENDIF
408 ELSEIF (ifail2(i) == 3) THEN
409 dlim = max(zero,e6(i,1)) / (dmx(i) + vfail)
410 ENDIF
411 ENDIF
412 IF (ifail(i) == 0) THEN
413! Uniaxial failure
414 crit(i) = max(crit(i),xa*dlim)
415 ELSE
416! Multiaxial failure
417 crit(i)= crit(i) + xa*dlim**xb
418 ENDIF
419 ENDIF
420 ENDDO
421C
422 DO i=1,nel
423 iadbuf = ipm(7,mid(i)) - 1
424 ifunc(i) = ipm(10 + if1 + 2,mid(i))
425 ifv(i) = ipm(10 + if2 + 2,mid(i))
426 ifunc2(i)= ipm(10 + if3 + 2,mid(i))
427 ifunc3(i)= ipm(10 + if4 + 2,mid(i))
428 ifunc4(i)= ipm(10 + if5 + 2,mid(i))
429 iecrou(i)= nint(uparam(iadbuf + i13 + 2))
430 ak(i) = uparam(iadbuf + i1 + 2)
431 b(i) = uparam(iadbuf + i2 + 2)
432 d(i) = uparam(iadbuf + i3 + 2)
433 ee(i) = uparam(iadbuf + i4 + 2)
434 gf3(i) = uparam(iadbuf + i5 + 2)
435 ff(i) = uparam(iadbuf + i6 + 2)
436 lscale(i)= uparam(iadbuf + i7 + 2)
437 dmn(i) = uparam(iadbuf + i8 + 2)
438 dmx(i) = uparam(iadbuf + i9 + 2)
439 ENDDO
440C
441 kk = 1 + numelr * anim_fe(11)
442 IF (nuvar >= 2) THEN
443 coord_old => uvar(2,1:nel)
444 ELSE
445 coord_old => not_used
446 ENDIF
447 CALL redef3(python,
448 . fy, yk, dy, fyep,
449 . dyold, dpy, tf, npf,
450 . yc, off, e6(1,2), dpy2,
451 . anim(kk), anim_fe(12),posy,
452 . xl0, dmn, dmx, dv,
453 . ff, lscale, ee, gf3,
454 . ifunc3, yieldy, aldp, ak,
455 . b, d, iecrou, ifunc,
456 . ifv, ifunc2, epla, coord_old,
457 . nel, nft, stf, sanin,
458 . dt1, iresp, impl_s, idyna,
459 . snpc, yieldc=yieldyc, xx_oldc=dyoldc,
460 . fx0=uparam(iadbuf+i15+2),ifunc4=ifunc4, pos6=posy(6,1))
461C
462 DO i=1,nel
463 iadbuf = ipm(7,mid(i)) - 1
464 cc = uparam(iadbuf + nupar + 4)
465 cn = uparam(iadbuf + nupar + 10)
466 xa = uparam(iadbuf + nupar + 16)
467 xb = uparam(iadbuf + nupar + 22)
468 IF (off(i) == one .AND. dmx(i) /= zero .AND. dmn(i) /= zero) THEN
469 IF (ifail2(i) == 0) THEN
470 xa = one
471 xb = two
472 IF (dy(i) > zero) THEN
473 dlim = dy(i) / (xl0(i)*dmx(i))
474 ELSE
475 dlim = dy(i) / (xl0(i)*dmn(i))
476 ENDIF
477 ELSE
478 vfail = cc * (abs(dv(i)/vrt(i)))**cn
479 IF (ifail2(i) == 1) THEN
480 IF (dy(i) > zero) THEN
481 dlim = dy(i) / (xl0(i)*dmx(i) + vfail)
482 ELSE
483 dlim = dy(i) / (xl0(i)*dmn(i) - vfail)
484 ENDIF
485 ELSEIF (ifail2(i) == 2) THEN
486 IF (fy(i) > zero) THEN
487 dlim = fy(i) / (dmx(i) + vfail)
488 ELSE
489 dlim = fy(i) / (dmn(i) - vfail)
490 ENDIF
491 ELSEIF (ifail2(i) == 3) THEN
492 dlim = max(zero,e6(i,2)) / (dmx(i) + vfail)
493 ENDIF
494 ENDIF
495 IF (ifail(i) == 0) THEN
496! Uniaxial failure
497 crit(i) = max(crit(i),xa*dlim)
498 ELSE
499! Multiaxial failure
500 crit(i)= crit(i) + xa*dlim**xb
501 ENDIF
502 ENDIF
503 ENDDO
504C
505 DO i=1,nel
506 iadbuf = ipm(7,mid(i)) - 1
507 ifunc(i) = ipm(10 + if1 + 3,mid(i))
508 ifv(i) = ipm(10 + if2 + 3,mid(i))
509 ifunc2(i)= ipm(10 + if3 + 3,mid(i))
510 ifunc3(i)= ipm(10 + if4 + 3,mid(i))
511 ifunc4(i)= ipm(10 + if5 + 3,mid(i))
512 iecrou(i)= nint(uparam(iadbuf + i13 + 3))
513 ak(i) = uparam(iadbuf + i1 + 3)
514 b(i) = uparam(iadbuf + i2 + 3)
515 d(i) = uparam(iadbuf + i3 + 3)
516 ee(i) = uparam(iadbuf + i4 + 3)
517 gf3(i) = uparam(iadbuf + i5 + 3)
518 ff(i) = uparam(iadbuf + i6 + 3)
519 lscale(i)= uparam(iadbuf + i7 + 3)
520 dmn(i) = uparam(iadbuf + i8 + 3)
521 dmx(i) = uparam(iadbuf + i9 + 3)
522 ENDDO
523C
524 kk = 1 + numelr * (anim_fe(11)+anim_fe(12))
525 IF (nuvar >= 3) THEN
526 coord_old => uvar(3,1:nel)
527 ELSE
528 coord_old => not_used
529 ENDIF
530 CALL redef3(python,
531 . fz, zk, dz, fzep,
532 . dzold, dpz, tf, npf,
533 . zc, off, e6(1,3), dpz2,
534 . anim(kk), anim_fe(13),posz,
535 . xl0, dmn, dmx, dv,
536 . ff, lscale, ee, gf3,
537 . ifunc3, yieldz, aldp, ak,
538 . b, d, iecrou, ifunc,
539 . ifv, ifunc2, epla, coord_old,
540 . nel, nft, stf, sanin,
541 . dt1, iresp, impl_s, idyna,
542 . snpc, yieldc=yieldzc, xx_oldc=dzoldc,
543 . fx0=uparam(iadbuf+i15+3),ifunc4=ifunc4, pos6=posz(6,1))
544C
545 DO i=1,nel
546 iadbuf = ipm(7,mid(i)) - 1
547 cc = uparam(iadbuf + nupar + 5)
548 cn = uparam(iadbuf + nupar + 11)
549 xa = uparam(iadbuf + nupar + 17)
550 xb = uparam(iadbuf + nupar + 23)
551 IF (off(i) == one .AND. dmx(i) /= zero .AND. dmn(i) /= zero) THEN
552 IF (ifail2(i) == 0) THEN
553 xa = one
554 xb = two
555 IF (dz(i) > zero)THEN
556 dlim = dz(i) / (xl0(i)*dmx(i))
557 ELSE
558 dlim = dz(i) / (xl0(i)*dmn(i))
559 ENDIF
560 ELSE
561 vfail = cc * (abs(dv(i)/vrt(i)))**cn
562 IF (ifail2(i) == 1) THEN
563 IF (dz(i) > zero) THEN
564 dlim = dz(i) / (xl0(i)*dmx(i) + vfail)
565 ELSE
566 dlim = dz(i) / (xl0(i)*dmn(i) - vfail)
567 ENDIF
568 ELSEIF (ifail2(i) == 2) THEN
569 IF (fz(i) > zero) THEN
570 dlim = fz(i) / (dmx(i) + vfail)
571 ELSE
572 dlim = fz(i) / (dmn(i) - vfail)
573 ENDIF
574 ELSEIF (ifail2(i) == 3) THEN
575 dlim = max(zero,e6(i,3)) / (dmx(i) + vfail)
576 ENDIF
577 ENDIF
578 IF (ifail(i) == 0) THEN
579! Uniaxial failure
580 crit(i) = max(crit(i),xa*dlim)
581 ELSE
582! Multiaxial failure
583 crit(i)= crit(i) + xa*dlim**xb
584 ENDIF
585 ENDIF
586 ENDDO
587C---------------------
588C ROTATIONS
589C---------------------
590 DO i=1,nel
591 iadbuf= ipm(7,mid(i)) - 1
592 xin(i)= geo(2,pid(i))
593 xk(i) = uparam(iadbuf + i11 + 4)
594 xc(i) = uparam(iadbuf + i12 + 4)
595 yk(i) = uparam(iadbuf + i11 + 5)
596 yc(i) = uparam(iadbuf + i12 + 5)
597 zk(i) = uparam(iadbuf + i11 + 6)
598 zc(i) = uparam(iadbuf + i12 + 6)
599 hx(i) = uparam(iadbuf + i14 + 4)
600 hy(i) = uparam(iadbuf + i14 + 5)
601 hz(i) = uparam(iadbuf + i14 + 6)
602
603 xhr(i)= max(hx(i),hy(i),hz(i))
604
605 xkr(i)= max(xk(i)*uparam(iadbuf + i1 + 4),
606 . yk(i)*uparam(iadbuf + i1 + 5),
607 . zk(i)*uparam(iadbuf + i1 + 6))+xkr(i)
608 xcr(i)= max(xc(i),yc(i),zc(i)) + xhr(i) +xcr(i)+xh(i)
609 ENDDO
610C
611 DO i=1,nel
612 dxold(i)=rx(i)
613 dyold(i)=ry(i)
614 dzold(i)=rz(i)
615 ENDDO
616!
617 IF ( inispri /= 0 .AND. tt == zero) THEN
618 DO i=1,nel
619 dxold(i)=rx0(i)
620 dyold(i)=ry0(i)
621 dzold(i)=rz0(i)
622 ENDDO
623 ENDIF
624C
625 DO i=1,nel
626 x21 = (rx2(i)-rx1(i))*dt1
627 y21 = (ry2(i)-ry1(i))*dt1
628 z21 = (rz2(i)-rz1(i))*dt1
629 rx(i) = dxold(i) + x21*exx2(i)+y21*eyx2(i)+z21*ezx2(i)
630 ry(i) = dyold(i) + x21*exy2(i)+y21*eyy2(i)+z21*ezy2(i)
631 rz(i) = dzold(i) + x21*exz2(i)+y21*eyz2(i)+z21*ezz2(i)
632 ENDDO
633C-------------------------------
634 DO i=1,nel
635 iadbuf = ipm(7,mid(i)) - 1
636 ifunc(i) = ipm(10 + if1 + 4,mid(i))
637 ifv(i) = ipm(10 + if2 + 4,mid(i))
638 ifunc2(i)= ipm(10 + if3 + 4,mid(i))
639 ifunc3(i)= ipm(10 + if4 + 4,mid(i))
640 ifunc4(i)= ipm(10 + if5 + 4,mid(i))
641 iecrou(i)= nint(uparam(iadbuf + i13 + 4))
642 ak(i) = uparam(iadbuf + i1 + 4)
643 b(i) = uparam(iadbuf + i2 + 4)
644 d(i) = uparam(iadbuf + i3 + 4)
645 ee(i) = uparam(iadbuf + i4 + 4)
646 gf3(i) = uparam(iadbuf + i5 + 4)
647 ff(i) = uparam(iadbuf + i6 + 4)
648 lscale(i)= uparam(iadbuf + i7 + 4)
649 dmn(i) = uparam(iadbuf + i8 + 4)
650 dmx(i) = uparam(iadbuf + i9 + 4)
651 ENDDO
652 IF (nuvar >= 4) THEN
653 coord_old => uvar(4,1:nel)
654 ELSE
655 coord_old => not_used
656 ENDIF
657 CALL redef3(python,
658 . xmom, xk, rx, xmep,
659 . dxold, rpx, tf, npf,
660 . xc, off, e6(1,4), rpx2,
661 . anim, 0, posxx,
662 . xl0, dmn, dmx, dv,
663 . ff, lscale, ee, gf3,
664 . ifunc3, yieldx2, aldp, ak,
665 . b, d, iecrou, ifunc,
666 . ifv, ifunc2, epla, coord_old,
667 . nel, nft, stf, sanin,
668 . dt1, iresp, impl_s, idyna,
669 . snpc, yieldc=yieldrxc, xx_oldc=drxoldc,
670 . fx0=uparam(iadbuf+i15+4),ifunc4=ifunc4, pos6=posxx(6,1))
671C
672 DO i=1,nel
673 iadbuf= ipm(7,mid(i)) - 1
674 cc = uparam(iadbuf + nupar + 6)
675 cn = uparam(iadbuf + nupar + 12)
676 xa = uparam(iadbuf + nupar + 18)
677 xb = uparam(iadbuf + nupar + 24)
678 IF (off(i) == one .AND. dmx(i) /= zero .AND. dmn(i) /= zero) THEN
679 IF (ifail2(i) == 0) THEN
680 xa = one
681 xb = two
682 IF (rx(i) > zero) THEN
683 dlim = rx(i) / (xl0(i)*dmx(i))
684 ELSE
685 dlim = rx(i) / (xl0(i)*dmn(i))
686 ENDIF
687 ELSE
688 vfail = cc * (abs(dv(i)/vrr(i)))**cn
689 IF (ifail2(i) == 1) THEN
690 IF (rx(i) > zero) THEN
691 dlim = rx(i) / (xl0(i)*dmx(i) + vfail)
692 ELSE
693 dlim = rx(i) / (xl0(i)*dmn(i) - vfail)
694 ENDIF
695 ELSEIF (ifail2(i) == 2) THEN
696 IF(xmom(i)>0.)THEN
697 dlim = xmom(i)/(dmx(i) + vfail)
698 ELSE
699 dlim = xmom(i)/(dmn(i) - vfail)
700 ENDIF
701 ELSEIF (ifail2(i) == 3) THEN
702 dlim = max(zero,e6(i,4)) / (dmx(i) + vfail)
703 ENDIF
704 ENDIF
705 IF (ifail(i) == 0) THEN
706! Uniaxial failure
707 crit(i) = max(crit(i),xa*dlim)
708 ELSE
709! Multiaxial failure
710 crit(i)= crit(i) + xa*dlim**xb
711 ENDIF
712 ENDIF
713 ENDDO
714C
715 DO i=1,nel
716 iadbuf = ipm(7,mid(i)) - 1
717 ifunc(i) = ipm(10 + if1 + 5,mid(i))
718 ifv(i) = ipm(10 + if2 + 5,mid(i))
719 ifunc2(i)= ipm(10 + if3 + 5,mid(i))
720 ifunc3(i)= ipm(10 + if4 + 5,mid(i))
721 ifunc4(i)= ipm(10 + if5 + 5,mid(i))
722 iecrou(i)= nint(uparam(iadbuf + i13 + 5))
723 ak(i) = uparam(iadbuf + i1 + 5)
724 b(i) = uparam(iadbuf + i2 + 5)
725 d(i) = uparam(iadbuf + i3 + 5)
726 ee(i) = uparam(iadbuf + i4 + 5)
727 gf3(i) = uparam(iadbuf + i5 + 5)
728 ff(i) = uparam(iadbuf + i6 + 5)
729 lscale(i)= uparam(iadbuf + i7 + 5)
730 dmn(i) = uparam(iadbuf + i8 + 5)
731 dmx(i) = uparam(iadbuf + i9 + 5)
732 ENDDO
733 IF (nuvar >= 5) THEN
734 coord_old => uvar(5,1:nel)
735 ELSE
736 coord_old => not_used
737 ENDIF
738 CALL redef3(python,
739 . ymom, yk, ry, ymep,
740 . dyold, rpy, tf, npf,
741 . yc, off, e6(1,5), rpy2,
742 . anim, 0, posyy,
743 . xl0, dmn, dmx, dv,
744 . ff, lscale, ee, gf3,
745 . ifunc3, yieldy2, aldp, ak,
746 . b, d, iecrou, ifunc,
747 . ifv, ifunc2, epla, coord_old,
748 . nel, nft, stf, sanin,
749 . dt1, iresp, impl_s, idyna,
750 . snpc, yieldc=yieldryc, xx_oldc=dryoldc,
751 . fx0=uparam(iadbuf+i15+5),ifunc4=ifunc4, pos6=posyy(6,1))
752C
753 DO i=1,nel
754 iadbuf= ipm(7,mid(i)) - 1
755 cc = uparam(iadbuf + nupar + 7)
756 cn = uparam(iadbuf + nupar + 13)
757 xa = uparam(iadbuf + nupar + 19)
758 xb = uparam(iadbuf + nupar + 25)
759 IF (off(i) == one .AND. dmx(i) /= zero .AND. dmn(i) /= zero) THEN
760 IF (ifail2(i) == 0) THEN
761 xa = one
762 xb = two
763 IF (ry(i) > zero) THEN
764 dlim = ry(i) / (xl0(i)*dmx(i))
765 ELSE
766 dlim = ry(i) / (xl0(i)*dmn(i))
767 ENDIF
768 ELSE
769 vfail = cc * (abs(dv(i)/vrr(i)))**cn
770 IF (ifail2(i) == 1) THEN
771 IF (ry(i) > zero) THEN
772 dlim = ry(i) / (xl0(i)*dmx(i) + vfail)
773 ELSE
774 dlim = ry(i) / (xl0(i)*dmn(i) - vfail)
775 ENDIF
776 ELSEIF (ifail2(i) == 2) THEN
777 IF (ymom(i) > zero)THEN
778 dlim = ymom(i)/(dmx(i) + vfail)
779 ELSE
780 dlim = ymom(i)/(dmn(i) - vfail)
781 ENDIF
782 ELSEIF (ifail2(i) == 3) THEN
783 dlim = max(zero,e6(i,5)) / (dmx(i) + vfail)
784 ENDIF
785 ENDIF
786 IF (ifail(i) == 0) THEN
787! Uniaxial failure
788 crit(i) = max(crit(i),xa*dlim)
789 ELSE
790! Multiaxial failure
791 crit(i)= crit(i) + xa*dlim**xb
792 ENDIF
793 ENDIF
794 ENDDO
795C
796 DO i=1,nel
797 iadbuf = ipm(7,mid(i)) - 1
798 ifunc(i) = ipm(10 + if1 + 6,mid(i))
799 ifv(i) = ipm(10 + if2 + 6,mid(i))
800 ifunc2(i)= ipm(10 + if3 + 6,mid(i))
801 ifunc3(i)= ipm(10 + if4 + 6,mid(i))
802 ifunc4(i)= ipm(10 + if5 + 6,mid(i))
803 iecrou(i)= nint(uparam(iadbuf + i13 + 6))
804 ak(i) = uparam(iadbuf + i1 + 6)
805 b(i) = uparam(iadbuf + i2 + 6)
806 d(i) = uparam(iadbuf + i3 + 6)
807 ee(i) = uparam(iadbuf + i4 + 6)
808 gf3(i) = uparam(iadbuf + i5 + 6)
809 ff(i) = uparam(iadbuf + i6 + 6)
810 lscale(i)= uparam(iadbuf + i7 + 6)
811 dmn(i) = uparam(iadbuf + i8 + 6)
812 dmx(i) = uparam(iadbuf + i9 + 6)
813 ENDDO
814 IF (nuvar >= 6) THEN
815 coord_old => uvar(6,1:nel)
816 ELSE
817 coord_old => not_used
818 ENDIF
819 CALL redef3(python,
820 . zmom, zk, rz, zmep,
821 . dzold, rpz, tf, npf,
822 . zc, off, e6(1,6), rpz2,
823 . anim, 0, poszz,
824 . xl0, dmn, dmx, dv,
825 . ff, lscale, ee, gf3,
826 . ifunc3, yieldz2, aldp, ak,
827 . b, d, iecrou, ifunc,
828 . ifv, ifunc2, epla, coord_old,
829 . nel, nft, stf, sanin,
830 . dt1, iresp, impl_s, idyna,
831 . snpc, yieldc=yieldrzc, xx_oldc=drzoldc,
832 . fx0=uparam(iadbuf+i15+6),ifunc4=ifunc4,pos6=poszz(6,1))
833C
834 DO i=1,nel
835 iadbuf= ipm(7,mid(i)) - 1
836 cc = uparam(iadbuf + nupar + 8)
837 cn = uparam(iadbuf + nupar + 14)
838 xa = uparam(iadbuf + nupar + 20)
839 xb = uparam(iadbuf + nupar + 26)
840 IF (off(i) == one .AND. dmx(i) /= zero .AND. dmn(i) /= zero) THEN
841 IF (ifail2(i) == 0) THEN
842 xa = one
843 xb = two
844 IF (rz(i) > zero) THEN
845 dlim = rz(i) / (xl0(i)*dmx(i))
846 ELSE
847 dlim = rz(i) / (xl0(i)*dmn(i))
848 ENDIF
849 ELSE
850 vfail = cc * (abs(dv(i)/vrr(i)))**cn
851 IF (ifail2(i) == 1) THEN
852 IF (rz(i) > zero)THEN
853 dlim = rz(i) / (xl0(i)*dmx(i) + vfail)
854 ELSE
855 dlim = rz(i) / (xl0(i)*dmn(i) - vfail)
856 ENDIF
857 ELSEIF (ifail2(i) == 2) THEN
858 IF (zmom(i) > zero)THEN
859 dlim = zmom(i)/(dmx(i) + vfail)
860 ELSE
861 dlim = zmom(i)/(dmn(i) - vfail)
862 ENDIF
863 ELSEIF (ifail2(i) == 3) THEN
864 dlim = max(zero,e6(i,6)) / (dmx(i) + vfail)
865 ENDIF
866 ENDIF
867 IF (ifail(i) == 0) THEN
868! Uniaxial failure
869 crit(i) = max(crit(i),xa*dlim)
870 ELSE
871! Multiaxial failure
872 crit(i)= crit(i) + xa *dlim**xb
873 ENDIF
874 ENDIF
875 ENDDO
876C
877 DO i=1,nel
878 e(i) = e6(i,1)+e6(i,2)+e6(i,3)+e6(i,4)+e6(i,5)+e6(i,6)
879 ENDDO
880C-------------------------------
881C COUPLED FAILURE
882C-------------------------------
883 DO i=1,nel
884 iadbuf = ipm(7,mid(i)) - 1
885 israte = nint(uparam(iadbuf + nupar + 27))
886C---- smoothing factor alpha = 2PI*fc*dt/(2PI*fc*dt+1) ---
887 asrate = (2*pi*uparam(iadbuf + nupar + 28)*dt1)/(one+2*pi*uparam(iadbuf + nupar + 28)*dt1)
888 IF (israte /= 0) THEN
889 IF (crit_new(i) < one) THEN
890 crit(i) = min(crit(i),one+em3)
891 crit(i) = asrate*crit(i) + (one - asrate)*crit_new(i)
892 crit_new(i) = min(crit(i),one)
893 ELSE
894 crit_new(i) = one
895 ENDIF
896 ELSE
897 IF (crit_new(i) < one) THEN
898 crit_new(i) = min(crit(i),one)
899 ELSE
900 crit_new(i) = one
901 ENDIF
902 ENDIF
903 IF (off(i) == one) THEN
904 IF (crit(i) >= one) THEN
905 off(i)=zero
906 nindx = nindx + 1
907 indx(nindx) = i
908 idel7nok = 1
909 ENDIF
910 ENDIF
911 ENDDO
912C
913 DO j=1,nindx
914 i = indx(j)
915#include "lockon.inc"
916 WRITE(iout, 1000) ngl(i)
917 WRITE(istdo,1100) ngl(i),tt
918#include "lockoff.inc"
919 ENDDO
920C-------------------------------
921C COUPLED PLASTICITY
922C-------------------------------
923 CALL repla3(
924 1 xk, rpx, tf, npf,
925 2 iecrou, ifunc, ifv, epla,
926 3 nel)
927 CALL repla3(
928 1 yk, rpy, tf, npf,
929 2 iecrou, ifunc, ifv, epla,
930 3 nel)
931 CALL repla3(
932 1 zk, rpz, tf, npf,
933 2 iecrou, ifunc, ifv, epla,
934 3 nel)
935C
936 DO i=1,nel
937 iadbuf= ipm(7,mid(i)) - 1
938 xk(i)=uparam(iadbuf + i11 + 1)
939 yk(i)=uparam(iadbuf + i11 + 2)
940 zk(i)=uparam(iadbuf + i11 + 3)
941 ENDDO
942C
943 CALL repla3(
944 1 xk, dpx, tf, npf,
945 2 iecrou, ifunc, ifv, epla,
946 3 nel)
947 CALL repla3(
948 1 yk, dpy, tf, npf,
949 2 iecrou, ifunc, ifv, epla,
950 3 nel)
951 CALL repla3(
952 1 zk, dpz, tf, npf,
953 2 iecrou, ifunc, ifv, epla,
954 3 nel)
955 DO i=1,nel
956 xkm(i)=xkm(i)/xl0(i)
957 xcm(i)=xcm(i)/xl0(i)
958 xin(i)=xin(i)*xl0(i)
959 xkr(i)=xkr(i)/xl0(i)
960 xcr(i)=xcr(i)/xl0(i)
961 ENDDO
962C---
963 1000 FORMAT(1x,'-- RUPTURE OF SPRING ELEMENT NUMBER ',i10)
964 1100 FORMAT(1x,'-- RUPTURE OF SPRING ELEMENT :',i10,' at time :',G11.4)
965C---
966 RETURN
967 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
r23l113def3
subroutine r23l113def3(python, skew, ipm, igeo, mid, pid, geo, uparam, fx, fy, fz, e, dx, dy, dz, npf, tf, off, dpx, dpy, dpz, dpx2, dpy2, dpz2, fxep, fyep, fzep, x0, y0, z0, xmom, ymom, zmom, rx, ry, rz, rpx, rpy, rpz, xmep, ymep, zmep, rpx2, rpy2, rpz2, anim, posx, posy, posz, posxx, posyy, poszz, fr_wave, e6, nel, exx2, eyx2, ezx2, exy2, eyy2, ezy2, exz2, eyz2, ezz2, al2dp, ngl, crit_new, x0_err, aldp, yieldx, yieldy, yieldz, yieldx2, yieldy2, yieldz2, exx, eyx, ezx, exy, eyy, ezy, exz, eyz, ezz, xcr, rx1, ry1, rz1, rx2, ry2, rz2, xin, ak, xm, xkm, xcm, xkr, vx1, vx2, vy1, vy2, vz1, vz2, nuvar, uvar, mass, dx0, dy0, dz0, rx0, ry0, rz0, nft, stf, sanin, iresp, snpc, szyield_comp, szxxold_comp, yieldxc, yieldyc, yieldzc, yieldrxc, yieldryc, yieldrzc, dxoldc, dyoldc, dzoldc, drxoldc, dryoldc, drzoldc)
Definition r23l113def3.F:68
repla3
subroutine repla3(xk, dpx, tf, npf, iecrou, ifunc, ifv, epla, nel)
Definition repla3.F:39