OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ruser44.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!|| ruser44 ../engine/source/elements/spring/ruser44.F
25!||--- called by ------------------------------------------------------
26!|| rforc3 ../engine/source/elements/spring/rforc3.F
27!||--- calls -----------------------------------------------------
28!|| get_u_func ../engine/source/user_interface/ufunc.F
29!|| get_u_mid ../engine/source/user_interface/upidmid.F
30!|| get_u_mnu ../engine/source/user_interface/upidmid.F
31!|| get_u_pid ../engine/source/user_interface/upidmid.F
32!|| get_u_pnu ../engine/source/user_interface/upidmid.F
33!|| get_v_func ../engine/source/user_interface/ufunc.F
34!||====================================================================
35 SUBROUTINE ruser44(NEL,IOUT ,IPROP ,UVAR ,NUVAR ,
36 2 FX ,FY ,FZ ,XMOM ,YMOM ,
37 3 ZMOM ,E ,OFF ,STIFM ,STIFR ,
38 4 VISCM ,VISCR ,MASS ,XINER ,DT ,
39 5 XL ,VX ,RY1 ,RZ1 ,RX ,
40 6 RY2 ,RZ2 ,FR_WAVE_E)
41C-------------------------------------------------------------------------
42C Crushable frame spring property
43C----------+---------+---+---+--------------------------------------------
44C VAR | SIZE |TYP| RW| DEFINITION
45C----------+---------+---+---+--------------------------------------------
46C IOUT | 1 | I | R | OUTPUT FILE UNIT (L00 file)
47C IPROP | 1 | I | R | PROPERTY NUMBER
48C----------+---------+---+---+--------------------------------------------
49C XL | NEL | F | R | ELEMENT LENGTH
50C----------+---------+---+---+--------------------------------------------
51C UVAR |NUVAR*NEL| F |R/W| USER ELEMENT VARIABLES
52C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
53C-------------------------------------------------------------------------
54C FUNCTION
55C-------------------------------------------------------------------------
56C INTEGER II = GET_U_PNU(I,IP,KK)
57C IFUNCI = GET_U_PNU(I,IP,KFUNC)
58C IPROPI = GET_U_PNU(I,IP,KFUNC)
59C IMATI = GET_U_PNU(I,IP,KMAT)
60C I : VARIABLE INDEX(1 for first variable,...)
61C IP : PROPERTY NUMBER
62C KK : PARAMETER KFUNC,KMAT,KPROP
63C THIS FUNCTION RETURN THE USER STORED FUNCTION(IF KK=KFUNC),
64C MATERIAL(IF KK=KMAT) OR PROPERTY(IF KK=KPROP) NUMBER.
65C SEE LECG29 FOR CORRESPONDING ID STORAGE.
66C-------------------------------------------------------------------------
67C INTEGER IFUNCI = GET_U_MNU(I,IM,KFUNC)
68C I : VARIABLE INDEX(1 for first function)
69C IM : MATERIAL NUMBER
70C KFUNC : ONLY FUNCTION ARE YET AVAILABLE.
71C THIS FUNCTION RETURN THE USER STORED FUNCTION NUMBER(function
72C referred by users materials).
73C SEE LECM29 FOR CORRESPONDING ID STORAGE.
74C-------------------------------------------------------------------------
75C my_real PARAMI = GET_U_GEO(I,IP)
76C I : PARAMETER INDEX(1 for first parameter,...)
77C IP : PROPERTY NUMBER
78C THIS FUNCTION RETURN THE USER GEOMETRY PARAMETERS
79C-------------------------------------------------------------------------
80C my_real PARAMI = GET_U_MAT(I,IM)
81C I : PARAMETER INDEX(1 for first parameter,...)
82C IM : MATERIAL NUMBER
83C THIS FUNCTION RETURN THE USER MATERIAL PARAMETERS
84C NOTE: GET_U_MAT(0,IMAT) RETURN THE DENSITY
85C-------------------------------------------------------------------------
86C INTEGER PID = GET_U_PID(IP)
87C IP : PROPERTY NUMBER
88C THIS FUNCTION RETURN THE USER PROPERTY ID CORRESPONDING TO
89C USER PROPERTY NUMBER IP.
90C-------------------------------------------------------------------------
91C INTEGER MID = GET_U_MID(IM)
92C IM : MATERIAL NUMBER
93C THIS FUNCTION RETURN THE USER MATERIAL ID CORRESPONDING TO
94C USER MATERIAL NUMBER IM.
95C-------------------------------------------------------------------------
96C my_real Y = GET_U_FUNC(IFUNC,X,DYDX)
97C IFUNC : function number obtained by
98C IFUNC = GET_U_MNU(I,IM,KFUNC) or IFUNC = GET_U_PNU(I,IP,KFUNC)
99C X : X value
100C DYDX : slope dY/dX
101C THIS FUNCTION RETURN Y(X)
102C-------------------------------------------------------------------------
103C I m p l i c i t T y p e s
104C-----------------------------------------------
105#include "implicit_f.inc"
106C----------------------------------------------------------
107C D u m m y A r g u m e n t s a n d F u n c t i o n
108C----------------------------------------------------------
109 INTEGER IOUT,NEL,NUVAR,IPROP,ICO,
110 . GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
111 . KFUNC,KMAT,KPROP
112 my_real
113 . UVAR(NUVAR,*),DT ,
114 . fx(*), fy(*), fz(*), e(*), vx(*),mass(*) ,xiner(*),
115 . ry1(*), rz1(*), off(*), xmom(*), ymom(*),
116 . zmom(*), rx(*), ry2(*), rz2(*),xl(*),
117 . stifm(*) ,stifr(*) , viscm(*) ,viscr(*) ,fr_wave_e(*) ,
118 . get_u_mat, get_u_geo, get_u_func
119 EXTERNAL get_u_mnu,get_u_pnu,get_u_mid,get_u_pid,
120 . get_u_mat,get_u_geo, get_u_func
121 parameter(kfunc=29)
122 parameter(kmat=31)
123 parameter(kprop=47)
124C-----------------------------------------------
125C L o c a l V a r i a b l e s
126C-----------------------------------------------
127 INTEGER I,ID1,ID2,ID10,ID20,
128 . ifun_xp,ifun_xmi,ifun_xxpi,ifun_xxmi,ifun_yy1pi,
129 . ifun_yy1mi,ifun_yy2pi,ifun_yy2mi,ifun_zz1pi,
130 . ifun_zz1mi,ifun_zz2pi,ifun_zz2mi,
131 . ifun_xmr,ifun_xxpr,ifun_xxmr,ifun_yy1pr,
132 . ifun_yy1mr,ifun_yy2pr,ifun_yy2mr,ifun_zz1pr,
133 . ifun_zz1mr,ifun_zz2pr,ifun_zz2mr,
134 . ifun_damp_xi,ifun_damp_yi,ifun_damp_zi,
135 . ifun_damp_xxi,ifun_damp_yyi,ifun_damp_zzi
136 INTEGER IFUNXM(NEL),IFUNXXM(NEL),IFUNYY1M(NEL),
137 . IFUNZZ1M(NEL),IFUNYY2M(NEL),IFUNZZ2M(NEL),
138 . IFUNXP(NEL),IFUNXXP(NEL),IFUNYY1P(NEL),
139 . IFUNZZ1P(NEL),IFUNYY2P(NEL),IFUNZZ2P(NEL),
140 . JPOSXM(NEL),JPOSXXM(NEL),JPOSYY1M(NEL),
141 . jposzz1m(nel),jposyy2m(nel),jposzz2m(nel),
142 . jposxp(nel),jposxxp(nel),jposyy1p(nel),
143 . jposzz1p(nel),jposyy2p(nel),jposzz2p(nel),
144 . jpos_damp_x(nel),jpos_damp_y(nel),jpos_damp_z(nel),
145 . jpos_damp_xx(nel),jpos_damp_yy(nel),jpos_damp_zz(nel),
146 . ifun_damp_x(nel),ifun_damp_y(nel),ifun_damp_z(nel),
147 . ifun_damp_xx(nel),ifun_damp_yy(nel),ifun_damp_zz(nel)
148 my_real ::
149 . xlimg,xlim,xxlim,yy1lim,yy2lim,zz1lim,zz2lim,
150 . k11,k44,k55,k66,k5b,k6c,leng2,leng,uleng0,r,ax,ay,az,
151 . fscal_x,fscal_rx,fscal_ry1,fscal_ry2,fscal_rz1,fscal_rz2,
152 . a,b,d,ff1,ff2,r1,mm,ff,fscal_damp_x,fscal_damp_y,
153 . fscal_damp_z,fscal_damp_xx,fscal_damp_yy,
154 . fscal_damp_zz,f_x,f_y,f_z,f_xx,f_yy,f_zz,alpha,ncf,
155 . idamping
156 my_real , DIMENSION(NEL) ::
157 . dydx,x,fxm,fxp,xx,fxxm,fxxp,yy1,fyy1m,fyy1p,
158 . zz1,fzz1m,fzz1p,yy2,fyy2m,fyy2p,zz2,fzz2m,fzz2p,
159 . uleng,my1,my2,mz1,mz2,fx_damp,fy_damp,fz_damp,
160 . fxx_damp,fyy_damp,fzz_damp,xdot,xxdot,yy1dot,zz1dot,yy2dot,zz2dot,
161 . xdot_old,xxdot_old,yy1dot_old,zz1dot_old,yy2dot_old,zz2dot_old,
162 . elodotx,elodoty,elodotz,elodotxx,elodotyy,elodotzz
163C=======================================================================
164 xlimg = get_u_geo(1,iprop)
165 xlim = get_u_geo(2,iprop)
166 xxlim = get_u_geo(3,iprop)
167 yy1lim = get_u_geo(4,iprop)
168 zz1lim = get_u_geo(5,iprop)
169 yy2lim = get_u_geo(6,iprop)
170 zz2lim = get_u_geo(7,iprop)
171 ico = nint(get_u_geo(16,iprop))
172 fscal_x = get_u_geo(17,iprop)
173 fscal_rx = get_u_geo(18,iprop)
174 fscal_ry1 = get_u_geo(19,iprop)
175 fscal_ry2 = get_u_geo(20,iprop)
176 fscal_rz1 = get_u_geo(21,iprop)
177 fscal_rz2 = get_u_geo(22,iprop)
178!------------------------------------------
179 ifun_xmi = get_u_pnu(1,iprop,kfunc)
180 ifun_xxmi = get_u_pnu(2,iprop,kfunc)
181 ifun_yy1mi= get_u_pnu(3,iprop,kfunc)
182 ifun_zz1mi= get_u_pnu(4,iprop,kfunc)
183 ifun_yy2mi= get_u_pnu(5,iprop,kfunc)
184 ifun_zz2mi= get_u_pnu(6,iprop,kfunc)
185 ifun_xp = get_u_pnu(7,iprop,kfunc)
186 ifun_xxpi = get_u_pnu(8,iprop,kfunc)
187 ifun_yy1pi= get_u_pnu(9,iprop,kfunc)
188 ifun_zz1pi= get_u_pnu(10,iprop,kfunc)
189 ifun_yy2pi= get_u_pnu(11,iprop,kfunc)
190 ifun_zz2pi= get_u_pnu(12,iprop,kfunc)
191 ifun_xmr = get_u_pnu(13,iprop,kfunc)
192 ifun_xxmr = get_u_pnu(14,iprop,kfunc)
193 ifun_yy1mr= get_u_pnu(15,iprop,kfunc)
194 ifun_zz1mr= get_u_pnu(16,iprop,kfunc)
195 ifun_yy2mr= get_u_pnu(17,iprop,kfunc)
196 ifun_zz2mr= get_u_pnu(18,iprop,kfunc)
197 ifun_xxpr = get_u_pnu(19,iprop,kfunc)
198 ifun_yy1pr= get_u_pnu(20,iprop,kfunc)
199 ifun_zz1pr= get_u_pnu(21,iprop,kfunc)
200 ifun_yy2pr= get_u_pnu(22,iprop,kfunc)
201 ifun_zz2pr= get_u_pnu(23,iprop,kfunc)
202!------------------------------------------
203! strain rate filtering number of cycles:
204 ncf = get_u_geo(35,iprop)
205! flag for damping activation
206 idamping = get_u_geo(36,iprop)
207!
208! add damping:
209!
210 IF (idamping > zero) THEN
211 ifun_damp_xi = get_u_pnu(24,iprop,kfunc)
212 ifun_damp_yi = get_u_pnu(25,iprop,kfunc)
213 ifun_damp_zi = get_u_pnu(26,iprop,kfunc)
214 ifun_damp_xxi = get_u_pnu(27,iprop,kfunc)
215 ifun_damp_yyi = get_u_pnu(28,iprop,kfunc)
216 ifun_damp_zzi = get_u_pnu(29,iprop,kfunc)
217!
218 fscal_damp_x = get_u_geo(23,iprop)
219 fscal_damp_y = get_u_geo(24,iprop)
220 fscal_damp_z = get_u_geo(25,iprop)
221 fscal_damp_xx = get_u_geo(26,iprop)
222 fscal_damp_yy = get_u_geo(27,iprop)
223 fscal_damp_zz = get_u_geo(28,iprop)
224! coefficients for function damping:
225 f_x = get_u_geo(29,iprop)
226 f_y = get_u_geo(30,iprop)
227 f_z = get_u_geo(31,iprop)
228 f_xx = get_u_geo(32,iprop)
229 f_yy = get_u_geo(33,iprop)
230 f_zz = get_u_geo(34,iprop)
231 ENDIF ! IF (IDAMPING > ZERO)
232!
233! add filtering:
234!
235 IF (ncf > zero) THEN
236 alpha = two/(ncf+one)
237!
238 DO i=1,nel
239 xdot_old(i) = uvar(37,i) ! old strain rate
240 xxdot_old(i) = uvar(38,i)
241 yy1dot_old(i) = uvar(39,i)
242 zz1dot_old(i) = uvar(40,i)
243 yy2dot_old(i) = uvar(41,i)
244 zz2dot_old(i) = uvar(42,i)
245 ENDDO
246!
247 DO i=1,nel
248 uleng0 = uvar(30,i)
249 xdot(i) = vx(i) * uleng0 ! current strain rate
250 xxdot(i) = rx(i)
251 yy1dot(i) = ry1(i)
252 zz1dot(i) = rz1(i)
253 yy2dot(i) = ry2(i)
254 zz2dot(i) = rz2(i)
255 ENDDO
256! filter current strain rate
257 DO i=1,nel
258 xdot(i) = alpha * xdot(i) + (one-alpha) * xdot_old(i)
259 xxdot(i) = alpha * xxdot(i) + (one-alpha) * xxdot_old(i)
260 yy1dot(i) = alpha * yy1dot(i) + (one-alpha) * yy1dot_old(i)
261 zz1dot(i) = alpha * zz1dot(i) + (one-alpha) * zz1dot_old(i)
262 yy2dot(i) = alpha * yy2dot(i) + (one-alpha) * yy2dot_old(i)
263 zz2dot(i) = alpha * zz2dot(i) + (one-alpha) * zz2dot_old(i)
264 ENDDO
265! save new filtered strain rate:
266 DO i=1,nel
267 uvar(37,i) = xdot(i) ! new filtered strain rate
268 uvar(38,i) = xxdot(i)
269 uvar(39,i) = yy1dot(i)
270 uvar(40,i) = zz1dot(i)
271 uvar(41,i) = yy2dot(i)
272 uvar(42,i) = zz2dot(i)
273 ENDDO
274! total filtered strain :
275 DO i=1,nel
276!! ULENG0 = UVAR(30,I)
277!! X(I) = UVAR(1,I) + DT * XDOT(I) * ULENG0 ! filtered strain
278 x(i) = uvar(1,i) + dt * xdot(i) ! filtered strain
279 xx(i) = uvar(2,i) + dt * xxdot(i)
280 yy1(i) = uvar(3,i) + dt * yy1dot(i)
281 zz1(i) = uvar(4,i) + dt * zz1dot(i)
282 yy2(i) = uvar(5,i) + dt * yy2dot(i)
283 zz2(i) = uvar(6,i) + dt * zz2dot(i)
284 ENDDO
285 ENDIF ! IF (NCF > ZERO)
286C=======================================================================
287C ELASTIQUE
288C=======================================================================
289 DO i=1,nel
290 my1(i) = -ymom(i)+half*xl(i)*fz(i)
291 my2(i) = ymom(i)+half*xl(i)*fz(i)
292 mz1(i) = -zmom(i)-half*xl(i)*fy(i)
293 mz2(i) = zmom(i)-half*xl(i)*fy(i)
294C
295 leng = max(xl(i),em20)
296 leng2 = leng*leng
297 uleng(i) = one/leng
298C
299 k11 = uvar(19,i)
300 k44 = uvar(20,i)
301 k55 = uvar(21,i)*leng2
302 k66 = uvar(22,i)*leng2
303 k5b = uvar(23,i)*leng2
304 k6c = uvar(24,i)*leng2
305C
306 fx(i) = fx(i) + dt * k11 * vx(i)
307 xmom(i)= xmom(i)+ dt * k44 * rx(i)
308 my1(i) = my1(i) + dt * (k55 * ry1(i) + k5b * ry2(i))
309 my2(i) = my2(i) + dt * (k5b * ry1(i) + k55 * ry2(i))
310 mz1(i) = mz1(i) + dt * (k66 * rz1(i) + k6c * rz2(i))
311 mz2(i) = mz2(i) + dt * (k6c * rz1(i) + k66 * rz2(i))
312C
313 uleng0 = uvar(30,i)
314C
315 IF (ncf == zero) THEN ! no filtering
316 x(i) = uvar(1,i) + dt * vx(i) * uleng0 ! strain
317 xx(i) = uvar(2,i) + dt * rx(i)
318 yy1(i) = uvar(3,i) + dt * ry1(i)
319 zz1(i) = uvar(4,i) + dt * rz1(i)
320 yy2(i) = uvar(5,i) + dt * ry2(i)
321 zz2(i) = uvar(6,i) + dt * rz2(i)
322 ENDIF ! IF (NCF == 0)
323C
324 jposxm(i) = nint(uvar(7,i))
325 jposxxm(i) = nint(uvar(8,i))
326 jposyy1m(i) = nint(uvar(9,i))
327 jposzz1m(i) = nint(uvar(10,i))
328 jposyy2m(i) = nint(uvar(11,i))
329 jposzz2m(i) = nint(uvar(12,i))
330 jposxp(i) = nint(uvar(13,i))
331 jposxxp(i) = nint(uvar(14,i))
332 jposyy1p(i) = nint(uvar(15,i))
333 jposzz1p(i) = nint(uvar(16,i))
334 jposyy2p(i) = nint(uvar(17,i))
335 jposzz2p(i) = nint(uvar(18,i))
336C
337 ifunxp(i) = ifun_xp
338 ifunxm(i) = ifun_xmi
339 ifunxxm(i) = ifun_xxmi
340 ifunxxp(i) = ifun_xxpi
341 ifunyy1m(i) = ifun_yy1mi
342 ifunyy1p(i) = ifun_yy1pi
343 ifunyy2m(i) = ifun_yy2mi
344 ifunyy2p(i) = ifun_yy2pi
345 ifunzz1m(i) = ifun_zz1mi
346 ifunzz1p(i) = ifun_zz1pi
347 ifunzz2m(i) = ifun_zz2mi
348 ifunzz2p(i) = ifun_zz2pi
349!
350! add function damping:
351!
352 IF (idamping > zero) THEN
353 ifun_damp_x(i) = ifun_damp_xi
354 ifun_damp_y(i) = ifun_damp_yi
355 ifun_damp_z(i) = ifun_damp_zi
356 ifun_damp_xx(i) = ifun_damp_xxi
357 ifun_damp_yy(i) = ifun_damp_yyi
358 ifun_damp_zz(i) = ifun_damp_zzi
359!
360 jpos_damp_x(i) = nint(uvar(31,i))
361 jpos_damp_y(i) = nint(uvar(32,i))
362 jpos_damp_z(i) = nint(uvar(33,i))
363 jpos_damp_xx(i) = nint(uvar(34,i))
364 jpos_damp_yy(i) = nint(uvar(35,i))
365 jpos_damp_zz(i) = nint(uvar(36,i))
366 ENDIF ! IF (IDAMPING > ZERO)
367 ENDDO
368!
369 IF (idamping > zero) THEN
370! elongation rate for damping (linear + function):
371 DO i=1,nel
372 elodotx(i) = vx(i) ! FX
373 elodoty(i) = - half * (rz2(i) + rz1(i)) * xl(i) ! FY
374 elodotz(i) = half * (ry2(i) + ry1(i)) * xl(i) ! FZ
375 elodotxx(i) = rx(i) ! XMOM
376 elodotyy(i) = ry2(i) - ry1(i) ! YMOM
377 elodotzz(i) = rz2(i) - rz1(i) ! ZMOM
378 ENDDO
379 ENDIF ! IF (IDAMPING > ZERO)
380C--------------------------------
381 IF (ico == 0) THEN ! => ELASTO PLASTIQUE - Classique
382C--------------------------------
383 DO i=1,nel
384 id1 = nint(uvar(29,i))/2
385 id2 = nint(uvar(29,i)) - 2*id1
386 id10=id1
387 id20=id2
388 IF (fr_wave_e(i) == one) THEN
389 jposxm(i) = 0
390 ifunxm(i) = ifun_xmr
391 ENDIF
392 IF (x(i) < xlimg) THEN !!! collapsed element
393 fr_wave_e(i)=one
394 ELSE
395 fr_wave_e(i)=zero
396 ENDIF
397C
398 IF (x(i) < xlim .OR. abs(xx(i)) > xxlim) THEN
399 id1 = 1
400 id2 = 1
401 ENDIF
402 IF (abs(yy1(i)) > yy1lim .OR. abs(zz1(i)) > zz1lim) id1 = 1
403 IF (abs(yy2(i)) > yy2lim .OR. abs(zz2(i)) > zz2lim) id2 = 1
404C
405 IF (id1 == 1) THEN
406 ifunxm(i) = ifun_xmr
407 ifunxxm(i) = ifun_xxmr
408 ifunxxp(i) = ifun_xxpr
409 ifunyy1m(i) = ifun_yy1mr
410 ifunyy1p(i) = ifun_yy1pr
411 ifunzz1m(i) = ifun_zz1mr
412 ifunzz1p(i) = ifun_zz1pr
413 IF (id10 == 0) THEN
414 jposxm(i) = 0
415 jposxxm(i) = 0
416 jposxxp(i) = 0
417 jposyy1m(i) = 0
418 jposzz1m(i) = 0
419 jposyy1p(i) = 0
420 jposzz1p(i) = 0
421 ENDIF
422 ENDIF
423 IF (id2 == 1) THEN
424 ifunxm(i) = ifun_xmr
425 ifunxxm(i) = ifun_xxmr
426 ifunxxp(i) = ifun_xxpr
427 ifunyy2m(i) = ifun_yy2mr
428 ifunyy2p(i) = ifun_yy2pr
429 ifunzz2m(i) = ifun_zz2mr
430 ifunzz2p(i) = ifun_zz2pr
431 IF (id20 == 0) THEN
432 jposxm(i) = 0
433 jposxxm(i) = 0
434 jposxxp(i) = 0
435 jposyy2m(i) = 0
436 jposzz2m(i) = 0
437 jposyy2p(i) = 0
438 jposzz2p(i) = 0
439 ENDIF
440 ENDIF
441 uvar(29,i) = 2*id1 + id2
442 ENDDO
443!
444 CALL get_v_func(ifunxm,nel,x,dydx,fxm,jposxm)
445 CALL get_v_func(ifunxp,nel,x,dydx,fxp,jposxp)
446 CALL get_v_func(ifunxxm,nel,xx,dydx,fxxm,jposxxm)
447 CALL get_v_func(ifunxxp,nel,xx,dydx,fxxp,jposxxp)
448 CALL get_v_func(ifunyy1m,nel,yy1,dydx,fyy1m,jposyy1m)
449 CALL get_v_func(ifunyy1p,nel,yy1,dydx,fyy1p,jposyy1p)
450 CALL get_v_func(ifunzz1m,nel,zz1,dydx,fzz1m,jposzz1m)
451 CALL get_v_func(ifunzz1p,nel,zz1,dydx,fzz1p,jposzz1p)
452 CALL get_v_func(ifunyy2m,nel,yy2,dydx,fyy2m,jposyy2m)
453 CALL get_v_func(ifunyy2p,nel,yy2,dydx,fyy2p,jposyy2p)
454 CALL get_v_func(ifunzz2m,nel,zz2,dydx,fzz2m,jposzz2m)
455 CALL get_v_func(ifunzz2p,nel,zz2,dydx,fzz2p,jposzz2p)
456!-------------------------------------------------------------
457!
458! add damping function interpolation:
459!
460 IF (idamping > zero) THEN
461 IF (ifun_damp_xi > 0) THEN
462 DO i=1,nel
463 elodotx(i) = elodotx(i) / f_x
464 ENDDO
465 CALL get_v_func(ifun_damp_x,nel,elodotx,dydx,fx_damp,jpos_damp_x)
466 ENDIF
467!
468 IF (ifun_damp_yi > 0) THEN
469 DO i=1,nel
470 elodoty(i) = elodoty(i) / f_y
471 ENDDO
472 CALL get_v_func(ifun_damp_y,nel,elodoty,dydx,fy_damp,jpos_damp_y)
473 ENDIF
474!
475 IF (ifun_damp_zi > 0) THEN
476 DO i=1,nel
477 elodotz(i) = elodotz(i) / f_z
478 ENDDO
479 CALL get_v_func(ifun_damp_z,nel,elodotz,dydx,fz_damp,jpos_damp_z)
480 ENDIF
481!
482 IF (ifun_damp_xxi > 0) THEN
483 DO i=1,nel
484 elodotxx(i) = elodotxx(i) / f_xx
485 ENDDO
486 CALL get_v_func(ifun_damp_xx,nel,elodotxx,dydx,fxx_damp,jpos_damp_xx)
487 ENDIF
488!
489 IF (ifun_damp_yyi > 0) THEN
490 DO i=1,nel
491 elodotyy(i) = elodotyy(i) / f_yy
492 ENDDO
493 CALL get_v_func(ifun_damp_yy,nel,elodotyy,dydx,fyy_damp,jpos_damp_yy)
494 ENDIF
495!
496 IF (ifun_damp_zzi > 0) THEN
497 DO i=1,nel
498 elodotzz(i) = elodotzz(i) / f_zz
499 ENDDO
500 CALL get_v_func(ifun_damp_zz,nel,elodotzz,dydx,fzz_damp,jpos_damp_zz)
501 ENDIF
502 ENDIF ! IF (IDAMPING > ZERO)
503!
504 DO i=1,nel
505 fxp(i) = fxp(i) * fscal_x
506 fxm(i) = fxm(i) * fscal_x
507 fxxp(i) = fxxp(i) * fscal_rx
508 fxxm(i) = fxxm(i) * fscal_rx
509 fyy1m(i) = fyy1m(i)* fscal_ry1
510 fyy1p(i) = fyy1p(i)* fscal_ry1
511 fzz1m(i) = fzz1m(i)* fscal_rz1
512 fzz1p(i) = fzz1p(i)* fscal_rz1
513 fyy2m(i) = fyy2m(i)* fscal_ry2
514 fyy2p(i) = fyy2p(i)* fscal_ry2
515 fzz2m(i) = fzz2m(i)* fscal_rz2
516 fzz2p(i) = fzz2p(i)* fscal_rz2
517 ENDDO
518!
519 DO i=1,nel
520 uvar(1,i) = x(i)
521 uvar(2,i) = xx(i)
522 uvar(3,i) = yy1(i)
523 uvar(4,i) = zz1(i)
524 uvar(5,i) = yy2(i)
525 uvar(6,i) = zz2(i)
526 uvar(7,i) = jposxm(i)
527 uvar(8,i) = jposxxm(i)
528 uvar(9,i) = jposyy1m(i)
529 uvar(10,i) = jposzz1m(i)
530 uvar(11,i) = jposyy2m(i)
531 uvar(12,i) = jposzz2m(i)
532 uvar(13,i) = jposxp(i)
533 uvar(14,i) = jposxxp(i)
534 uvar(15,i) = jposyy1p(i)
535 uvar(16,i) = jposzz1p(i)
536 uvar(17,i) = jposyy2p(i)
537 uvar(18,i) = jposzz2p(i)
538!---
539! add damping:
540!---
541 IF (idamping > zero) THEN
542 uvar(31,i) = jpos_damp_x(i)
543 uvar(32,i) = jpos_damp_y(i)
544 uvar(33,i) = jpos_damp_z(i)
545 uvar(34,i) = jpos_damp_xx(i)
546 uvar(35,i) = jpos_damp_yy(i)
547 uvar(36,i) = jpos_damp_zz(i)
548 ENDIF ! IF IDAMPING > ZERO)
549!---
550 fx(i) = max(min(fx(i) ,fxp(i)),fxm(i) )
551 xmom(i) = max(min(xmom(i),fxxp(i)),fxxm(i))
552 my1(i) = max(min(my1(i),fyy1p(i)),fyy1m(i))
553 mz1(i) = max(min(mz1(i),fzz1p(i)),fzz1m(i))
554 my2(i) = max(min(my2(i),fyy2p(i)),fyy2m(i))
555 mz2(i) = max(min(mz2(i),fzz2p(i)),fzz2m(i))
556 fz(i) = (my1(i)+my2(i)) * uleng(i)
557 fy(i) = -(mz1(i)+mz2(i)) * uleng(i)
558 ymom(i) = -half*(my1(i)-my2(i))
559 zmom(i) = -half*(mz1(i)-mz2(i))
560C
561 viscm(i) = zero
562 viscr(i) = zero
563 stifm(i) = uvar(25,i)
564 stifr(i) = uvar(26,i)
565 mass(i) = uvar(27,i)
566 xiner(i) = uvar(28,i)
567!
568! add damping term (linear or function):
569!
570 IF (idamping > zero) THEN
571 IF (ifun_damp_x(i) == 0) THEN ! ! linear damping
572 fx(i) = fx(i) + fscal_damp_x * elodotx(i)
573 ELSE
574 fx(i) = fx(i) + fscal_damp_x * fx_damp(i)
575 ENDIF
576!
577 IF (ifun_damp_y(i) == 0) THEN ! ! linear damping
578 fy(i) = fy(i) + fscal_damp_y * elodoty(i)
579 ELSE
580 fy(i) = fy(i) + fscal_damp_y * fy_damp(i)
581 ENDIF
582!
583 IF (ifun_damp_z(i) == 0) THEN ! ! linear damping
584 fz(i) = fz(i) + fscal_damp_z * elodotz(i)
585 ELSE
586 fz(i) = fz(i) + fscal_damp_z * fz_damp(i)
587 ENDIF
588!
589 IF (ifun_damp_xx(i) == 0) THEN ! ! linear damping
590 xmom(i) = xmom(i) + fscal_damp_xx * elodotxx(i)
591 ELSE
592 xmom(i) = xmom(i) + fscal_damp_xx * fxx_damp(i)
593 ENDIF
594!
595 IF (ifun_damp_yy(i) == 0) THEN ! ! linear damping
596 ymom(i) = ymom(i) + fscal_damp_yy * elodotyy(i)
597 ELSE
598 ymom(i) = ymom(i) + fscal_damp_yy * fyy_damp(i)
599 ENDIF
600!
601 IF (ifun_damp_zz(i) == 0) THEN ! ! linear damping
602 zmom(i) = zmom(i) + fscal_damp_zz * elodotzz(i)
603 ELSE
604 zmom(i) = zmom(i) + fscal_damp_zz * fzz_damp(i)
605 ENDIF
606 ENDIF ! IF (IDAMPING > ZERO)
607!
608 ENDDO
609C--------------------------------
610 ELSE ! => ELASTO PLASTIQUE - Directions couplees
611C--------------------------------
612 DO i=1,nel
613 id1 = nint(uvar(29,i))/2
614 id2 = nint(uvar(29,i)) - 2*id1
615 id10=id1
616 id20=id2
617 IF (fr_wave_e(i) == one) THEN !!! collapsed element
618 jposxm(i) = 0
619 ifunxm(i) = ifun_xmr
620 ENDIF
621 IF (x(i) < xlimg) THEN
622 fr_wave_e(i)=one
623 ELSE
624 fr_wave_e(i)=zero
625 ENDIF
626C
627 IF (x(i) < xlim .OR. abs(xx(i)) > xxlim) THEN
628 id1 = 1
629 id2 = 1
630 ENDIF
631 IF (abs(yy1(i)) > yy1lim .OR. abs(zz1(i)) > zz1lim) id1 = 1
632 IF (abs(yy2(i)) > yy2lim .OR. abs(zz2(i)) > zz2lim) id2 = 1
633C
634 IF (id1 == 1) THEN
635 ifunxm(i) = ifun_xmr
636 ifunxxm(i) = ifun_xxmr
637 ifunxxp(i) = ifun_xxpr
638 ifunyy1m(i) = ifun_yy1mr
639 ifunyy1p(i) = ifun_yy1pr
640 ifunzz1m(i) = ifun_zz1mr
641 ifunzz1p(i) = ifun_zz1pr
642 IF (id10 == 0) THEN
643 jposxm(i) = 0
644 jposxxm(i) = 0
645 jposxxp(i) = 0
646 jposyy1m(i) = 0
647 jposzz1m(i) = 0
648 jposyy1p(i) = 0
649 jposzz1p(i) = 0
650 ENDIF
651 ENDIF
652 IF (id2 == 1) THEN
653 ifunxm(i) = ifun_xmr
654 ifunxxm(i) = ifun_xxmr
655 ifunxxp(i) = ifun_xxpr
656 ifunyy2m(i) = ifun_yy2mr
657 ifunyy2p(i) = ifun_yy2pr
658 ifunzz2m(i) = ifun_zz2mr
659 ifunzz2p(i) = ifun_zz2pr
660 IF (id20 == 0) THEN
661 jposxm(i) = 0
662 jposxxm(i) = 0
663 jposxxp(i) = 0
664 jposyy2m(i) = 0
665 jposzz2m(i) = 0
666 jposyy2p(i) = 0
667 jposzz2p(i) = 0
668 ENDIF
669 ENDIF
670 uvar(29,i) = 2*id1 + id2
671 ENDDO
672C
673 CALL get_v_func(ifunxm,nel,x,dydx,fxm,jposxm)
674 CALL get_v_func(ifunxp,nel,x,dydx,fxp,jposxp)
675 CALL get_v_func(ifunxxm,nel,xx,dydx,fxxm,jposxxm)
676 CALL get_v_func(ifunxxp,nel,xx,dydx,fxxp,jposxxp)
677 CALL get_v_func(ifunyy1m,nel,yy1,dydx,fyy1m,jposyy1m)
678 CALL get_v_func(ifunyy1p,nel,yy1,dydx,fyy1p,jposyy1p)
679 CALL get_v_func(ifunzz1m,nel,zz1,dydx,fzz1m,jposzz1m)
680 CALL get_v_func(ifunzz1p,nel,zz1,dydx,fzz1p,jposzz1p)
681 CALL get_v_func(ifunyy2m,nel,yy2,dydx,fyy2m,jposyy2m)
682 CALL get_v_func(ifunyy2p,nel,yy2,dydx,fyy2p,jposyy2p)
683 CALL get_v_func(ifunzz2m,nel,zz2,dydx,fzz2m,jposzz2m)
684 CALL get_v_func(ifunzz2p,nel,zz2,dydx,fzz2p,jposzz2p)
685!
686! add damping function interpolation:
687!
688 IF (idamping > zero) THEN
689 IF (ifun_damp_xi > 0) THEN
690 DO i=1,nel
691 elodotx(i) = elodotx(i) / f_x
692 ENDDO
693 CALL get_v_func(ifun_damp_x,nel,elodotx,dydx,fx_damp,jpos_damp_x)
694 ENDIF
695!
696 IF (ifun_damp_yi > 0) THEN
697 DO i=1,nel
698 elodoty(i) = elodoty(i) / f_y
699 ENDDO
700 CALL get_v_func(ifun_damp_y,nel,elodoty,dydx,fy_damp,jpos_damp_y)
701 ENDIF
702!
703 IF (ifun_damp_zi > 0) THEN
704 DO i=1,nel
705 elodotz(i) = elodotz(i) / f_z
706 ENDDO
707 CALL get_v_func(ifun_damp_z,nel,elodotz,dydx,fz_damp,jpos_damp_z)
708 ENDIF
709!
710 IF (ifun_damp_xxi > 0) THEN
711 DO i=1,nel
712 elodotxx(i) = elodotxx(i) / f_xx
713 ENDDO
714 CALL get_v_func(ifun_damp_xx,nel,elodotxx,dydx,fxx_damp,jpos_damp_xx)
715 ENDIF
716!
717 IF (ifun_damp_yyi > 0) THEN
718 DO i=1,nel
719 elodotyy(i) = elodotyy(i) / f_yy
720 ENDDO
721 CALL get_v_func(ifun_damp_yy,nel,elodotyy,dydx,fyy_damp,jpos_damp_yy)
722 ENDIF
723!
724 IF (ifun_damp_zzi > 0) THEN
725 DO i=1,nel
726 elodotzz(i) = elodotzz(i) / f_zz
727 ENDDO
728 CALL get_v_func(ifun_damp_zz,nel,elodotzz,dydx,fzz_damp,jpos_damp_zz)
729 ENDIF
730 ENDIF ! IF (IDAMPING > ZERO)
731!
732 DO i=1,nel
733 fxp(i) = fxp(i) * fscal_x
734 fxm(i) = fxm(i) * fscal_x
735 fxxp(i) = fxxp(i) * fscal_rx
736 fxxm(i) = fxxm(i) * fscal_rx
737 fyy1m(i) = fyy1m(i)* fscal_ry1
738 fyy1p(i) = fyy1p(i)* fscal_ry1
739 fzz1m(i) = fzz1m(i)* fscal_rz1
740 fzz1p(i) = fzz1p(i)* fscal_rz1
741 fyy2m(i) = fyy2m(i)* fscal_ry2
742 fyy2p(i) = fyy2p(i)* fscal_ry2
743 fzz2m(i) = fzz2m(i)* fscal_rz2
744 fzz2p(i) = fzz2p(i)* fscal_rz2
745 ENDDO
746C
747 DO i=1,nel
748 uvar(1,i) = x(i)
749 uvar(2,i) = xx(i)
750 uvar(3,i) = yy1(i)
751 uvar(4,i) = zz1(i)
752 uvar(5,i) = yy2(i)
753 uvar(6,i) = zz2(i)
754 uvar(7,i) = jposxm(i)
755 uvar(8,i) = jposxxm(i)
756 uvar(9,i) = jposyy1m(i)
757 uvar(10,i) = jposzz1m(i)
758 uvar(11,i) = jposyy2m(i)
759 uvar(12,i) = jposzz2m(i)
760 uvar(13,i) = jposxp(i)
761 uvar(14,i) = jposxxp(i)
762 uvar(15,i) = jposyy1p(i)
763 uvar(16,i) = jposzz1p(i)
764 uvar(17,i) = jposyy2p(i)
765 uvar(18,i) = jposzz2p(i)
766!---
767! add damping:
768!---
769 IF (idamping > zero) THEN
770 uvar(31,i) = jpos_damp_x(i)
771 uvar(32,i) = jpos_damp_y(i)
772 uvar(33,i) = jpos_damp_z(i)
773 uvar(34,i) = jpos_damp_xx(i)
774 uvar(35,i) = jpos_damp_yy(i)
775 uvar(36,i) = jpos_damp_zz(i)
776 ENDIF ! IF IDAMPING > ZERO)
777!
778 my1(i) = max(min(my1(i),fyy1p(i)),fyy1m(i))
779 mz1(i) = max(min(mz1(i),fzz1p(i)),fzz1m(i))
780 my2(i) = max(min(my2(i),fyy2p(i)),fyy2m(i))
781 mz2(i) = max(min(mz2(i),fzz2p(i)),fzz2m(i))
782!
783 xmom(i) = max(min(xmom(i),fxxp(i)),fxxm(i))
784!
785 fyy1p(i) = max(em20, fyy1p(i))
786 fyy1m(i) = -max(em20,-fyy1m(i))
787 fyy2p(i) = max(em20, fyy2p(i))
788 fyy2m(i) = -max(em20,-fyy2m(i))
789!
790 fzz1p(i) = max(em20, fzz1p(i))
791 fzz1m(i) = -max(em20,-fzz1m(i))
792 fzz2p(i) = max(em20, fzz2p(i))
793 fzz2m(i) = -max(em20,-fzz2m(i))
794!
795 fxxp(i) = max(em20, fxxp(i))
796 fxxm(i) = -max(em20,-fxxm(i))
797!
798C---------------------------------------------------
799C Mo = R M
800C Mn = Mo [l^2 Fo^2 - l^2Fo F + Mo M ]/[l^2Fo^2+Mo^2]
801C Fn = Fo [ Mo^2 - Mo M + l^2Fo F ]/[l^2Fo^2+Mo^2]
802C Fn = Fo [ 1 - Mn/Mo ]
803C Fn = Fo [ 1 - [ l^2Fo^2 - l^2Fo F + Mo M ]/[l^2Fo^2+Mo^2]]
804C Fn = Fo [ 1 - [ l^2Fo^2 - l^2Fo F + R M^2 ]/[l^2Fo^2+R^2M^2]]
805C Fn = Fo [ 1 - [ l^2Fo(Fo - F) + R M^2 ]/[l^2Fo^2+R^2M^2]]
806C Mn = Mo [ 1 - Fn/Fo ]
807C Rn = Mn/M = Mo/M [ 1 - Fn/Fo ]
808C Rn = R [ 1 - Fn/Fo ]
809C si l = Mo/Fo
810C Fn = Fo [ 1/2 + F/2Fo - 1/2R ]
811C Rn = R [ 1 - Fn/Fo ] ! twice written
812C---------------------------------------------------
813!
814! - NODE 1 -
815!
816 ax = max(xmom(i)/fxxp(i),xmom(i)/fxxm(i))
817 ay = max(my1(i)/fyy1p(i),my1(i)/fyy1m(i))
818 az = max(mz1(i)/fzz1p(i),mz1(i)/fzz1m(i))
819!
820 mm = sqrt(ax*ax+ay*ay+az*az)
821!
822 ff = fx(i)/min(-em20,fxm(i))
823 r = (one+mm-ff)/max(em20,two*mm)
824 IF (r > zero .AND. r < one) THEN
825 ff1 = half*(one-mm+ff)*fxm(i)
826 ELSE
827 ff1 = fxm(i)
828 r = max(zero,min(one,r))
829 ENDIF
830 r1 = r
831!
832! - NODE 2 -
833!
834 ay = max(my2(i)/fyy2p(i),my2(i)/fyy2m(i))
835 az = max(mz2(i)/fzz2p(i),mz2(i)/fzz2m(i))
836 mm = sqrt(ax*ax+ay*ay+az*az)
837 r = (one+mm-ff)/max(em20,two*mm)
838 IF (r > zero .AND. r < one) THEN
839 ff2 = half*(one-mm+ff)*fxm(i)
840 ELSE
841 ff2 = fxm(i)
842 r = max(zero,min(one,r))
843 ENDIF
844 my1(i) = r1*my1(i)
845 mz1(i) = r1*mz1(i)
846 my2(i) = r*my2(i)
847 mz2(i) = r*mz2(i)
848C
849 xmom(i) = min(r,r1)*xmom(i)
850C
851 fx(i) = max(min(fx(i),fxp(i)), ff1, ff2)
852C
853 fz(i) = (my1(i)+my2(i)) * uleng(i)
854 fy(i) = -(mz1(i)+mz2(i)) * uleng(i)
855 ymom(i) = -half*(my1(i)-my2(i))
856 zmom(i) = -half*(mz1(i)-mz2(i))
857C
858 viscm(i) = zero
859 viscr(i) = zero
860 stifm(i) = uvar(25,i)
861 stifr(i) = uvar(26,i)
862 mass(i) = uvar(27,i)
863 xiner(i) = uvar(28,i)
864!
865! add damping term (linear or function):
866!
867 IF (idamping > zero) THEN
868 IF (ifun_damp_x(i) == 0) THEN ! ! linear damping
869 fx(i) = fx(i) + fscal_damp_x * elodotx(i)
870 ELSE
871 fx(i) = fx(i) + fscal_damp_x * fx_damp(i)
872 ENDIF
873!
874 IF (ifun_damp_y(i) == 0) THEN ! ! linear damping
875 fy(i) = fy(i) + fscal_damp_y * elodoty(i)
876 ELSE
877 fy(i) = fy(i) + fscal_damp_y * fy_damp(i)
878 ENDIF
879!
880 IF (ifun_damp_z(i) == 0) THEN ! ! linear damping
881 fz(i) = fz(i) + fscal_damp_z * elodotz(i)
882 ELSE
883 fz(i) = fz(i) + fscal_damp_z * fz_damp(i)
884 ENDIF
885!
886 IF (ifun_damp_xx(i) == 0) THEN ! ! linear damping
887 xmom(i) = xmom(i) + fscal_damp_xx * elodotxx(i)
888 ELSE
889 xmom(i) = xmom(i) + fscal_damp_xx * fxx_damp(i)
890 ENDIF
891!
892 IF (ifun_damp_yy(i) == 0) THEN ! ! linear damping
893 ymom(i) = ymom(i) + fscal_damp_yy * elodotyy(i)
894 ELSE
895 ymom(i) = ymom(i) + fscal_damp_yy * fyy_damp(i)
896 ENDIF
897!
898 IF (ifun_damp_zz(i) == 0) THEN ! ! linear damping
899 zmom(i) = zmom(i) + fscal_damp_zz * elodotzz(i)
900 ELSE
901 zmom(i) = zmom(i) + fscal_damp_zz * fzz_damp(i)
902 ENDIF
903 ENDIF ! IF (IDAMPING > ZERO)
904!
905 ENDDO
906C-------------------------------
907 ENDIF
908C-------------------------------
909 RETURN
910 END
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
ruser44
subroutine ruser44(nel, iout, iprop, uvar, nuvar, fx, fy, fz, xmom, ymom, zmom, e, off, stifm, stifr, viscm, viscr, mass, xiner, dt, xl, vx, ry1, rz1, rx, ry2, rz2, fr_wave_e)
Definition ruser44.F:41
get_v_func
subroutine get_v_func(ifunc, llt, xx, dydx, yy, jpos)
Definition ufunc.F:230