OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
r27def3.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!|| r27def3 ../engine/source/elements/spring/r27def3.F
25!||--- called by ------------------------------------------------------
26!|| rforc3 ../engine/source/elements/spring/rforc3.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||--- uses -----------------------------------------------------
30!|| python_funct_mod ../common_source/modules/python_mod.F90
31!||====================================================================
32 SUBROUTINE r27def3(PYTHON,
33 1 F, E, DL, AL0,
34 2 IPOS, GEO, IGEO, NPF,
35 3 TF, V, OFF, ANIM,
36 4 AL0_ERR, X1DP, X2DP, NGL,
37 5 MGN, EX, EY, EZ,
38 6 XK, XM, XC, FSCALE,
39 7 NEL, NFT, CRIT)
40 USE python_funct_mod
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45#include "comlock.inc"
46C-----------------------------------------------
47C G l o b a l P a r a m e t e r s
48C-----------------------------------------------
49#include "mvsiz_p.inc"
50C-----------------------------------------------
51C C o m m o n B l o c k s
52C-----------------------------------------------
53#include "units_c.inc"
54#include "scr17_c.inc"
55#include "scr14_c.inc"
56#include "param_c.inc"
57#include "com08_c.inc"
58C-----------------------------------------------
59C D u m m y A r g u m e n t s
60C-----------------------------------------------
61 type(python_), intent(inout) :: PYTHON
62 INTEGER, INTENT(IN) :: NEL
63 INTEGER, INTENT(IN) :: NFT
64 INTEGER NPF(*),IGEO(NPROPGI,*),NGL(*),MGN(*)
65C REAL
66 my_real
67 . GEO(NPROPG,*),F(*),AL0(*),E(*),DL(*),TF(*),OFF(*),
68 . anim(*),ipos(*),v(3,*),
69 . al0_err(*),ex(mvsiz),ey(mvsiz),ez(mvsiz),xk(mvsiz),
70 . xm(mvsiz),xc(mvsiz),fscale(mvsiz)
71 DOUBLE PRECISION X1DP(3,*),X2DP(3,*)
72 my_real, DIMENSION(NEL), INTENT(INOUT) ::
73 . crit
74C-----------------------------------------------
75C L o c a l V a r i a b l e s
76C-----------------------------------------------
77 INTEGER I,J,NINDX,PID
78 INTEGER NC1(MVSIZ), NC2(MVSIZ),INDX(MVSIZ),IFUNC(MVSIZ),
79 . ifunc2(mvsiz),ifail(mvsiz),ileng(mvsiz),itens(mvsiz),
80 . fsmooth(mvsiz)
81C REAL
82 my_real
83 . dlold(mvsiz),dmn(mvsiz),dmx(mvsiz),xl0(mvsiz),df1,df2,
84 . dvl(mvsiz),fk(mvsiz),fd(mvsiz),ddl(mvsiz),gap(mvsiz),
85 . fscale2(mvsiz),ascale(mvsiz),ascale2(mvsiz),fold(mvsiz),
86 . nexp(mvsiz),fcut(mvsiz),alpha(mvsiz),fkold(mvsiz)
87 my_real
88 . bid,sum,dt11,damp,damm
89 DOUBLE PRECISION EXDP(MVSIZ),EYDP(MVSIZ),EZDP(MVSIZ),ALDP(MVSIZ),
90 . AL0DP(MVSIZ)
91 LOGICAL :: ANY_PYTHON_FUNCTION
92 INTEGER :: NFUNCT,PYID
93C-----------------------------------------------
94C E x t e r n a l
95C-----------------------------------------------
96 my_real
97 . finter
98 EXTERNAL finter
99C-----------------------------------------------
100c
101C ==================================================================
102C RECOVERING SPRING PARAMETERS
103C ==================================================================
104 any_python_function = .false.
105 nfunct = python%FUNCT_OFFSET + python%nb_functs - python%nb_sensors! OFFSET = nb of non-python functions
106
107 DO i=1,nel
108 pid = mgn(i)
109 xm(i) = geo(1,pid) ! Spring mass
110 xk(i) = geo(2,pid) ! Spring linear stiffness
111 xc(i) = geo(3,pid) ! Spring linear damping
112 fscale(i) = geo(10,pid) ! Stiffness tabulated function scale factor
113 dmn(i) = geo(15,pid) ! Negative limit for failure
114 dmx(i) = geo(16,pid) ! Positive limit for failure
115 ascale2(i) = geo(18,pid) ! Velocity scale factor for damping tabulated function
116 gap(i) = geo(19,pid) ! Compression gap for spring activation
117 ascale(i) = geo(39,pid) ! Displacement scale factor for stiffness tabulated function
118 ifail(i) = nint(geo(43,pid)) ! Flag for failure criterion
119 ileng(i) = nint(geo(93,pid)) ! Flag for unit length values
120 fscale2(i) = geo(132,pid) ! Damping tabulated function scale factor
121 itens(i) = nint(geo(133,pid)) ! Tensile behavior flag
122 nexp(i) = geo(134,pid) ! Non-linear exponent
123 fsmooth(i) = nint(geo(135,pid)) ! Spring velocity filtering flag
124 fcut(i) = geo(136,pid) ! Cutoff frequency
125 ifunc(i) = igeo(101,pid) ! Function ID for stiffness tabulated force
126 ifunc2(i) = igeo(102,pid) ! Function ID for damping tabulated force
127 IF (ifunc2(i) /= 0) xc(i) = geo(141,pid)
128 pyid = python_funct_id(nfunct, ifunc(i),npf)
129 IF(pyid > 0) THEN
130 any_python_function = .true.
131 ifunc(i) = -pyid
132 ENDIF
133 pyid = python_funct_id(nfunct, ifunc2(i),npf)
134 IF(pyid > 0) THEN
135 any_python_function = .true.
136 ifunc2(i) = -pyid
137 ENDIF
138 ENDDO
139C ==================================================================
140C
141C ==================================================================
142C COMPUTATION OF THE NEW SPRING LENGTH
143C ==================================================================
144 DO i=1,nel
145 exdp(i) = x2dp(1,i)-x1dp(1,i)
146 eydp(i) = x2dp(2,i)-x1dp(2,i)
147 ezdp(i) = x2dp(3,i)-x1dp(3,i)
148 dlold(i) = dl(i)
149 aldp(i) = sqrt(exdp(i)*exdp(i)+eydp(i)*eydp(i)+ezdp(i)*ezdp(i))
150 ENDDO
151C
152 ! Conversion DOUBLE -> MY_REAL if TT = ZERO
153 IF (tt == zero) THEN
154 DO i=1,nel
155 al0(i) = aldp(i) ! cast double to My_real
156 al0_err(i) = aldp(i)-al0(i) ! difference between double and My_real
157 ENDDO
158 ENDIF
159C
160 DO i=1,nel
161 al0dp(i) = al0(i) ! cast My_real to double
162 al0dp(i) = al0dp(i) + al0_err(i) ! AL_DP doit etre recalcule ainsi afin de garantir la coherence absolue entre AL0_DP et AL_DP
163 ENDDO
164C
165 DO i=1,nel
166 sum = max(aldp(i),em15)
167 exdp(i) = exdp(i)/sum
168 eydp(i) = eydp(i)/sum
169 ezdp(i) = ezdp(i)/sum
170 ex(i) = exdp(i)
171 ey(i) = eydp(i)
172 ez(i) = ezdp(i)
173 ENDDO
174C
175 ! Spring total elongation
176 DO i=1,nel
177 dl(i) = aldp(i) - al0dp(i)
178 ENDDO
179C
180 ! ILENG flag
181 DO i=1,nel
182 ! Value per length unit
183 IF (ileng(i) /= 0) THEN
184 xl0(i) = al0dp(i)
185 ! Classical units
186 ELSE
187 xl0(i) = one
188 ENDIF
189 ENDDO
190C ==================================================================
191C
192C ==================================================================
193C COMPUTATION OF THE SPRING FORCE
194C ==================================================================
195 ! Timestep
196 dt11 = dt1
197 IF (dt11 == zero) dt11 = ep30
198C
199 ! Recovering spring variables + Unit conversion if needed (ILENG)
200 DO i = 1,nel
201 ! Current total elongation
202 dl(i) = dl(i)/xl0(i)
203 ! Old elongation
204 dlold(i) = dlold(i)/xl0(i)
205 ! Elongation increment
206 ddl(i) = (dl(i)-dlold(i))
207 ! Elongation velocity
208 dvl(i) = ddl(i)/dt11
209 ! Spring energy
210 e(i) = e(i)/xl0(i)
211 ! Save old force value
212 fold(i) = f(i)
213 ENDDO
214C
215C ------------------------------------------------------------------
216C LOOP OVER THE ELEMENTS
217C ------------------------------------------------------------------
218 DO i=1,nel
219c
220 ! Computation if the element is not broken only
221 IF ((off(i) == one).AND.(dl(i) < gap(i) .OR. itens(i) > 0)) THEN
222 ! Stiffness part
223 ! Computation of the tabulated stiffness force
224 IF (ifunc(i) > 0) THEN
225 fk(i) = fscale(i)*finter(ifunc(i),(dl(i)-gap(i))/ascale(i),npf,tf,df1)
226 ELSEIF(ifunc(i) < 0) THEN
227 CALL python_call_funct1d(python, -ifunc(i),(dl(i)-gap(i))/ascale(i), fk(i))
228 ELSE
229 IF (abs(dl(i)-gap(i)) > zero) THEN
230 fk(i) = sign(one,(dl(i)-gap(i)))*xk(i)*exp(nexp(i)*log(abs(dl(i)-gap(i))))
231 ELSE
232 fk(i) = zero
233 ENDIF
234 ! If non-linear, the stiffness must be re-evaluated for the critical timestep
235 IF (nexp(i) > one) THEN
236 ! Old stiffness force
237 IF (abs(dlold(i)-gap(i)) > zero) THEN
238 fkold(i) = sign(one,(dlold(i)-gap(i)))*xk(i)*exp(nexp(i)*log(abs(dlold(i)-gap(i))))
239 ELSE
240 fkold(i) = zero
241 ENDIF
242 ! Non-linear stiffness slope
243 xk(i) = max(abs((fk(i)-fkold(i))/sign(max(abs(ddl(i)),em20),ddl(i))),xk(i))
244 ENDIF
245 ENDIF
246c
247 ! Damping part
248 ! Computation of the tabulated damping force
249 IF (ifunc2(i) > 0) THEN
250 fd(i) = fscale2(i)*finter(ifunc2(i),dvl(i)/ascale2(i),npf,tf,df2)
251 ! Computation of the linear damping
252 ELSEIF(ifunc2(i) < 0) THEN
253 CALL python_call_funct1d(python, -ifunc2(i),dvl(i)/ascale2(i), fd(i))
254 ELSE
255 fd(i) = xc(i)*dvl(i)
256 ENDIF
257c
258 ! Assembling forces
259 IF (abs(fk(i)) > abs(fd(i))) THEN
260 f(i) = fk(i) + fd(i)
261 ELSE
262 f(i) = two*fk(i)
263 xk(i) = two*xk(i)
264 xc(i) = zero
265 ENDIF
266 ELSE
267 f(i) = zero
268 xc(i) = zero
269 ENDIF
270c
271 ! Spring force filtering
272 IF (fsmooth(i) > 0) THEN
273 alpha(i) = (two*pi*dt11*fcut(i))/(two*pi*dt11*fcut(i) + one)
274 f(i) = alpha(i)*f(i) + (one - alpha(i))*fold(i)
275 ENDIF
276c
277 ENDDO
278C
279 ! Computation of damage variable for animations
280 IF (anim_fe(11) /= 0) THEN
281 DO i=1,nel
282 j = i+nft
283 IF (ifail(i) == 1) THEN
284 IF (itens(i) > 0) THEN
285 damp = dl(i)/max(dmx(i),em15)
286 ELSE
287 damp = zero
288 ENDIF
289 damm = dl(i)/min(dmn(i),-em15)
290 ELSEIF (ifail(i) == 2) THEN
291 IF (itens(i) > 0) THEN
292 damp = f(i)/max(dmx(i),em15)
293 ELSE
294 damp = zero
295 ENDIF
296 damm = f(i)/min(dmn(i),-em15)
297 ELSE
298 damp = zero
299 damm = zero
300 ENDIF
301 anim(j) = max(anim(j),damp,damm)
302 anim(j) = min(anim(j),one)
303 ENDDO
304 ENDIF
305c
306 ! Computation of the Spring energy + Unit conversion if needed (ILENG)
307 DO i=1,nel
308 e(i) = e(i) + (dl(i)-dlold(i))*(f(i)+fold(i)) * half
309 e(i) = e(i)*xl0(i)
310 dl(i) = dl(i)*xl0(i)
311 dlold(i) = dlold(i)*xl0(i)
312 xm(i) = xm(i)*xl0(i)
313 xk(i) = xk(i)/xl0(i)
314 xc(i) = xc(i)/xl0(i)
315 ENDDO
316C ==================================================================
317C
318C ==================================================================
319C CHECKING BROKEN SPRINGS
320C ==================================================================
321 nindx = 0
322 DO i=1,nel
323 IF (off(i) == one .AND. ifail(i) /= 0) THEN
324 IF (ifail(i) == 1) THEN
325 IF (itens(i) > 0) THEN
326 crit(i) = max(dl(i)/(dmx(i)*xl0(i)),dl(i)/(dmn(i)*xl0(i)))
327 crit(i) = min(crit(i),one)
328 crit(i) = max(crit(i),zero)
329 IF (dl(i) > dmx(i)*xl0(i) .OR. dl(i) < dmn(i)*xl0(i)) THEN
330 crit(i) = one
331 off(i) = zero
332 nindx = nindx + 1
333 indx(nindx) = i
334 idel7nok = 1
335 ENDIF
336 ELSE
337 crit(i) = dl(i)/(dmn(i)*xl0(i))
338 crit(i) = min(crit(i),one)
339 crit(i) = max(crit(i),zero)
340 IF (dl(i) < dmn(i)*xl0(i)) THEN
341 crit(i) = one
342 off(i) = zero
343 nindx = nindx + 1
344 indx(nindx) = i
345 idel7nok = 1
346 ENDIF
347 ENDIF
348 ELSEIF (ifail(i) == 2) THEN
349 IF (itens(i) > 0) THEN
350 crit(i) = max(f(i)/dmx(i),f(i)/dmn(i))
351 crit(i) = min(crit(i),one)
352 crit(i) = max(crit(i),zero)
353 IF (f(i) > dmx(i) .OR. f(i) < dmn(i)) THEN
354 crit(i) = one
355 off(i) = zero
356 nindx = nindx + 1
357 indx(nindx) = i
358 idel7nok = 1
359 ENDIF
360 ELSE
361 crit(i) = f(i)/dmn(i)
362 crit(i) = min(crit(i),one)
363 crit(i) = max(crit(i),zero)
364 IF (f(i) < dmn(i)) THEN
365 crit(i) = one
366 off(i) = zero
367 nindx = nindx + 1
368 indx(nindx) = i
369 idel7nok = 1
370 ENDIF
371 ENDIF
372 ENDIF
373 ENDIF
374 ENDDO
375 DO j=1,nindx
376 i = indx(j)
377#include "lockon.inc"
378 WRITE(iout, 1000) ngl(i)
379 WRITE(istdo,1100) ngl(i),tt
380#include "lockoff.inc"
381 ENDDO
382C-----------
383 1000 FORMAT(1x,'-- RUPTURE OF SPRING ELEMENT NUMBER ',i10)
384 1100 FORMAT(1x,'-- RUPTURE OF SPRING ELEMENT :',i10,' AT TIME :',g11.4)
385C-----------
386 END SUBROUTINE r27def3
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
r27def3
subroutine r27def3(python, f, e, dl, al0, ipos, geo, igeo, npf, tf, v, off, anim, al0_err, x1dp, x2dp, ngl, mgn, ex, ey, ez, xk, xm, xc, fscale, nel, nft, crit)
Definition r27def3.F:40