OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
xforc28.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine xforc28 (nx, xel, vel, vrel, uix, uid, iout, iprop, imat, off, keint, eint, mass, xiner, stifm, stifr, viscm, viscr, forc, torq, nuvar, uvar, nuvarn, uvarn, dt, dte)

Function/Subroutine Documentation

◆ xforc28()

subroutine xforc28 ( integer nx,
xel,
vel,
vrel,
integer, dimension(nx) uix,
integer uid,
integer iout,
integer iprop,
integer imat,
off,
integer keint,
eint,
mass,
xiner,
stifm,
stifr,
viscm,
viscr,
forc,
torq,
integer nuvar,
uvar,
integer nuvarn,
uvarn,
dt,
dte )

Definition at line 38 of file xforc28.F.

45C-------------------------------------------------------------------------
46C This subroutine computes multi-purpose element forces and moments
47C when element uses property TYPE28==NSTRAND.
48C-------------------------------------------------------------------------
49C----------+---------+---+---+--------------------------------------------
50C VAR | SIZE |TYP| RW| DEFINITION
51C----------+---------+---+---+--------------------------------------------
52C NX | 1 | I | R | NUMBER OF NODES (CONSTANT IN THE GROUP)
53C----------+---------+---+---+--------------------------------------------
54C XEL | 3*NX | F | R | NODES COORDINATES
55C VEL | 3*NX | F | R | NODES VELOCITIES
56C VREL | 3*NX | F | R | NODES ROTATIONAL VELOCITIES
57C----------+---------+---+---+--------------------------------------------
58C UIX | NX | I | R | ELEMENT CONNECTIVITY
59C | IX(J) (1<=J<=NX) : NODE J USER ID
60C UID | 1 | I | R | ELEMENT USER IDENTIFIER
61C----------+---------+---+---+--------------------------------------------
62C IOUT | 1 | I | R | OUTPUT FILE UNIT (L00 file)
63C IPROP | 1 | I | R | PROPERTY NUMBER
64C IMAT | 1 | I | R | MATERIAL NUMBER
65C----------+---------+---+---+--------------------------------------------
66C OFF | 1 | F | W | ELEMENT STATE (ON=1./OFF=0.)
67C KEINT | 1 | I | W | ELEMENT INTERNAL ENERGY FLAG
68C | 0 | HAS TO BE COMPUTED BY RADIOSS OUT OF USER ROUTINE.
69C | 1 | IS RETURNED BY THIS USER ROUTINE.
70C EINT | 1 | F | W | ELEMENT INTERNAL ENERGY IF KEINT=1
71C----------+---------+---+---+--------------------------------------------
72C MASS | NX | F | W | NODAL MASS
73C XINER | NX | F | W | NODAL INERTIA (SPHERICAL)
74C STIFM | NX | F | W | NODAL STIFNESS (TIME STEP)
75C STIFR | NX | F | W | NODAL ROTATION STIFNESS (TIME STEP)
76C VISCM | NX | F | W | NODAL VISCOSITY (TIME STEP)
77C VISCR | NX | F | W | NODAL ROTATION VISCOSITY (TIME STEP)
78C----------+---------+---+---+--------------------------------------------
79C FORC | 3*NX | F | W | NODAL FORCES
80C TORQ | 3*NX | F | W | NODAL TORQUS
81C----------+---------+---+---+--------------------------------------------
82C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
83C UVAR |NUVAR | F |R/W| USER ELEMENT VARIABLES
84C | (FIX SIZE ZONE)
85C NUVARN | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES PER NODE
86C UVARN |NUVARN*NX| F |R/W| USER ELEMENT VARIABLES PER NODE
87C | (NX DEPENDENT SIZE ZONE)
88C----------+---------+---+---+--------------------------------------------
89C DT | 1 | F | R | PREVIOUS TIME INCREMENT :
90C | CURRENT TIME POINT IS T = TOLD+DT
91C----------+---------+---+---+--------------------------------------------
92C DTE | 1 | F | W | ELEMENT TIME STEP
93C----------+---------+---+---+--------------------------------------------
94C-------------------------------------------------------------------------
95C FUNCTION
96C-------------------------------------------------------------------------
97C INTEGER II = GET_U_PNU(I,IP,KK)
98C IFUNCI = GET_U_PNU(I,IP,KFUNC)
99C IPROPI = GET_U_PNU(I,IP,KFUNC)
100C IMATI = GET_U_PNU(I,IP,KMAT)
101C I : VARIABLE INDEX(1 for first variable,...)
102C IP : PROPERTY NUMBER
103C KK : PARAMETER KFUNC,KMAT,KPROP
104C THIS FUNCTION RETURN THE USER STORED FUNCTION(IF KK=KFUNC),
105C MATERIAL(IF KK=KMAT) OR PROPERTY(IF KK=KPROP) NUMBER.
106C SEE LECG28 FOR CORRESPONDING ID STORAGE.
107C-------------------------------------------------------------------------
108C INTEGER IFUNCI = GET_U_MNU(I,IM,KFUNC)
109C I : VARIABLE INDEX(1 for first function)
110C IM : MATERIAL NUMBER
111C KFUNC : ONLY FUNCTION ARE YET AVAILABLE.
112C THIS FUNCTION RETURN THE USER STORED FUNCTION NUMBER(function
113C referred by users materials).
114C SEE LECM28 FOR CORRESPONDING ID STORAGE.
115C-------------------------------------------------------------------------
116C my_real PARAMI = GET_U_GEO(I,IP)
117C I : PARAMETER INDEX(1 for first parameter,...)
118C IP : PROPERTY NUMBER
119C THIS FUNCTION RETURN THE USER GEOMETRY PARAMETERS
120C NOTE: IF(IP==IPROP) UPARAG(I) == GET_U_GEO(I,IPROP)
121C-------------------------------------------------------------------------
122C my_real PARAMI = GET_U_MAT(I,IM)
123C I : PARAMETER INDEX(1 for first parameter,...)
124C IM : MATERIAL NUMBER
125C THIS FUNCTION RETURN THE USER MATERIAL PARAMETERS
126C NOTE: GET_U_MAT(0,IMAT) RETURN THE DENSITY
127C-------------------------------------------------------------------------
128C INTEGER PID = GET_U_PID(IP)
129C IP : PROPERTY NUMBER
130C THIS FUNCTION RETURN THE USER PROPERTY ID CORRESPONDING TO
131C USER PROPERTY NUMBER IP.
132C-------------------------------------------------------------------------
133C INTEGER MID = GET_U_MID(IM)
134C IM : MATERIAL NUMBER
135C THIS FUNCTION RETURN THE USER MATERIAL ID CORRESPONDING TO
136C USER MATERIAL NUMBER IM.
137C-------------------------------------------------------------------------
138C my_real Y = GET_U_FUNC(IFUNC,X,DXDY)
139C IFUNC : function number obtained by
140C IFUNC = GET_U_MNU(I,IM,KFUNC) or IFUNC = GET_U_PNU(I,IP,KFUNC)
141C X : X value
142C DXDY : slope dX/dY
143C THIS FUNCTION RETURN Y(X)
144C-------------------------------------------------------------------------
145C-----------------------------------------------
146C M o d u l e s
147C-----------------------------------------------
148 USE message_mod
149C-----------------------------------------------
150C I m p l i c i t T y p e s
151C-----------------------------------------------
152#include "implicit_f.inc"
153#include "comlock.inc"
154C----------------------------------------------------------
155C 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
156C----------------------------------------------------------
157 INTEGER IOUT,NUVAR,NUVARN,IPROP,IMAT,
158 . NX ,UIX(NX) ,UID, KEINT,
159 . GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
160 . KFUNC,KMAT,KPROP
161 my_real
162 . off, eint, xel(3,nx), vel(3,nx) ,vrel(3,nx),
163 . mass(nx) , xiner(nx) ,stifm(nx) ,
164 . stifr(nx), viscm(nx) ,viscr(nx) ,
165 . forc(3,nx), torq(3,nx),
166 . uvar(nuvar),uvarn(nuvarn*nx),dt ,dte ,alpha,
167 . get_u_mat, get_u_geo, get_u_func
168 EXTERNAL get_u_mnu,get_u_pnu,get_u_mid,get_u_pid,
169 . get_u_mat,get_u_geo, get_u_func
170 parameter(kfunc=29)
171 parameter(kmat=31)
172 parameter(kprop=33)
173C=======================================================================
174C-----------------------------------------------
175C L o c a l V a r i a b l e s
176C-----------------------------------------------
177 INTEGER NB1, NB2, NB3, MB1, MB2, MB3, MB4, MB5,
178 . K, IFUNCT, IFV, NIT, ECHANGE
179 my_real
180 . lprev, lnext, vprev(3), vnext(3), vv,xk, xc,
181 . l0, l, dx, dxold, ddx, dt11, dvx, fx, g, y, dfdx, dgdx, f,
182 . stif, mu1, mu2, fric, mu, beta, fmax, dfs, df, ff,
183 . epsmoy, df1, df2, xn, fxold, epstot, epsmin, epsmax, deps,
184 . rho, xm, xkm, xks, xcm, dtc, dtk, fact, ffac,xfac
185C-----------------------------------------------
186C initial total length UVAR(NB1:NB1)
187 nb1=1
188C previous elongation UVAR(NB2:NB2)
189 nb2=nb1+1
190C previous force average UVAR(NB3:NB3)
191 nb3=nb2+1
192C initial nodes masses UVARN(MB1:MB1+NX-1)
193 mb1=1
194C forces into strands UVARN(MB2:MB2+NX-1)
195C using UVARN(MB2:MB2+NX-2) only.
196 mb2=mb1+nx
197C strands initial length UVARN(MB3:MB3+NX-1)
198C using UVARN(MB3:MB3+NX-2) only.
199 mb3=mb2+nx
200C strands elongations UVARN(MB4:MB4+NX-1)
201C using UVARN(MB4:MB4+NX-2) only.
202 mb4=mb3+nx
203C strands internal energy UVARN(MB5:MB5+NX-1)
204C using UVARN(MB5:MB5+NX-2) only.
205 mb5=mb4+nx
206C-----------------------------------------------
207C already broken only returns mass & null stiffness.
208 IF (off==zero) GOTO 200
209C-----------------------------------------------
210C FORCE & PLASTIC DISPLACEMENT,
211C-----------------------------------------------
212C ACTUAL LENGTH.
213 l=0.0
214 DO k=1,nx-1
215 lnext=
216 . sqrt((xel(1,k+1)-xel(1,k))*(xel(1,k+1)-xel(1,k))
217 . +(xel(2,k+1)-xel(2,k))*(xel(2,k+1)-xel(2,k))
218 . +(xel(3,k+1)-xel(3,k))*(xel(3,k+1)-xel(3,k)))
219 l=l+lnext
220 mass(k) =lnext
221 ENDDO
222C-------
223 xk =get_u_geo(4,iprop)
224 ifunct =get_u_pnu(1,iprop,kfunc)
225 ifv =get_u_pnu(2,iprop,kfunc)
226 ffac=get_u_geo(12,iprop)
227 xfac=get_u_geo(13,iprop)
228C-------
229C STIF (tangent stiffness) is considered for friction ..
230 IF (ifunct==0.AND.ifv==0) THEN
231C------------------------------------------
232C LINEAR ELASTIC
233 l0 =uvar(nb1)
234 dx =l-l0
235C raideur tgte par unite de longueur.
236 stif=xk
237 f =xk*dx/l0
238 epstot=dx/l0
239 ELSEIF (ifunct==0.AND.ifv/=0) THEN
240C------------------------------------------
241C G(Deps) only.
242 l0 =uvar(nb1)
243 dx =l-l0
244C no friction when only viscous ..
245 stif=0.
246 f = ffac
247 epstot=dx/l0
248 ELSE
249C------------------------------------------
250C NON LINEAR ELASTIC
251 l0 =uvar(nb1)
252 dx =l-l0
253 epstot=dx/l0
254 f =get_u_func(ifunct,epstot,dfdx)
255C raideur tgte par unite de longueur.
256 f = ffac*f
257 stif=ffac*dfdx
258 ENDIF
259C-------------------------------------
260C LINEAR DAMPING.
261 dt11 = dt
262 IF(dt11==zero)dt11 = ep30
263C
264 dxold=uvar(nb2)
265 ddx=dx-dxold
266 dvx=ddx/dt11
267C
268 deps =(epstot-dxold/l0)/dt11
269C
270 uvar(nb2)=dx
271C
272 xc =get_u_geo(5,iprop)
273 dgdx=zero
274 IF (ifv>0) THEN
275 g = get_u_func(ifv,deps*xfac,dgdx)
276 dgdx = dgdx*xfac
277 fx = f*g+xc*deps
278 stif=stif*g
279 ELSE
280 fx= f+xc*deps
281 g = one
282 ENDIF
283C-------
284 stif=stif*l/l0
285C-----------------------------------------------
286C FRICTION PULLEYS.
287C STRAND STIFFNESS == K*/l0 is assumed to be K*/l * L/L0
288C-----------------------------------------------
289C TORQ is used for temporary storage (computation optimization).
290 DO k=1,nx-1
291 vnext(1)=xel(1,k+1)-xel(1,k)
292 vnext(2)=xel(2,k+1)-xel(2,k)
293 vnext(3)=xel(3,k+1)-xel(3,k)
294 vv=1./max(em15,mass(k))
295 vnext(1)=vv*vnext(1)
296 vnext(2)=vv*vnext(2)
297 vnext(3)=vv*vnext(3)
298 torq(1,k)=vnext(1)
299 torq(2,k)=vnext(2)
300 torq(3,k)=vnext(3)
301 ENDDO
302C-------
303C temporary storage for internal energy computation.
304 DO k=1,nx-1
305 viscr(k)=uvarn(mb2+k-1)
306 xiner(k)=uvarn(mb4+k-1)
307 ENDDO
308C-------------------------------------
309C rupture criteria is based upon average of total deformation
310C into nstrand element.
311 epsmin=get_u_geo(8,iprop)
312 epsmax=get_u_geo(9,iprop)
313 IF (epstot<epsmin.OR.epstot>epsmax) THEN
314 off=zero
315 fx =zero
316#include "lockon.inc"
317 WRITE(iout, 1000) uid
318#include "lockoff.inc"
319 GOTO 100
320 ENDIF
321C-------
322 mu1 = get_u_geo(10,iprop)
323 mu2 = get_u_geo(11,iprop)
324C also available at time 0.
325 epsmoy=(dx-dxold)/max(em15,l)
326C-------
327 xn=nx
328 DO k=1,nx-1
329 vprev(1)=torq(1,k)
330 vprev(2)=torq(2,k)
331 vprev(3)=torq(3,k)
332 ddx = vprev(1) * (vel(1,k+1) - vel(1,k))
333 . + vprev(2) * (vel(2,k+1) - vel(2,k))
334 . + vprev(3) * (vel(3,k+1) - vel(3,k))
335 uvarn(mb2+k-1)=uvarn(mb2+k-1)
336C tous les brins ont la meme longueur correspond a :
337 . +stif*(dt*ddx*(xn-one)/max(em15,l)-epsmoy)
338 ENDDO
339C-------
340C
341 echange=1
342 nit=0
343 DO WHILE(echange==1.AND.nit<10)
344 echange=0
345 nit=nit+1
346 DO k=1,nx-1
347C stifr utilise de maniere temporaire.
348 stifr(k)=uvarn(mb2+k-1)
349 END DO
350 DO k=2,nx-1
351C default pulley friction.
352 fric=mu1
353 mu = get_u_geo(50+k,iprop)
354 IF (mu>=zero) fric=mu
355 mu = get_u_geo(525+k-1,iprop)
356 IF (mu>=zero) THEN
357 fric=fric + half*mu
358 ELSE
359 fric=fric + half*mu2
360 ENDIF
361 mu = get_u_geo(525+k,iprop)
362 IF (mu>=zero) THEN
363 fric=fric+half*mu
364 ELSE
365 fric=fric + half*mu2
366 ENDIF
367C
368 vprev(1)=torq(1,k-1)
369 vprev(2)=torq(2,k-1)
370 vprev(3)=torq(3,k-1)
371 vnext(1)=-torq(1,k)
372 vnext(2)=-torq(2,k)
373 vnext(3)=-torq(3,k)
374 alpha=vprev(1)*vnext(1)+vprev(2)*vnext(2)+vprev(3)*vnext(3)
375 alpha = min(alpha,one)
376 alpha = max(alpha,-one)
377 beta = pi-acos(alpha)
378C
379 ff = two*fx+stifr(k-1)+stifr(k)
380 fmax = max(zero,ff*tanh(half*fric*beta))
381 dfs = stifr(k-1)-stifr(k)
382 IF(abs(dfs)>fmax)THEN
383 df =sign(abs(dfs)-fmax,dfs)
384C tous les brins ont la meme longueur correspond a :
385 uvarn(mb2+k-2)=uvarn(mb2+k-2)-half*df
386 uvarn(mb2+k-1)=uvarn(mb2+k-1)+half*df
387 echange=1
388 END IF
389 END DO
390 END DO
391C-----------------------------------------------
392C RETURN FORCES.
393C-----------------------------------------------
394 vprev(1)=torq(1,1)
395 vprev(2)=torq(2,1)
396 vprev(3)=torq(3,1)
397 forc(1,1)=fx*vprev(1)
398 forc(2,1)=fx*vprev(2)
399 forc(3,1)=fx*vprev(3)
400 DO k=2,nx-1
401 vnext(1)=torq(1,k)
402 vnext(2)=torq(2,k)
403 vnext(3)=torq(3,k)
404 forc(1,k)=fx*(vnext(1)-vprev(1))
405 forc(2,k)=fx*(vnext(2)-vprev(2))
406 forc(3,k)=fx*(vnext(3)-vprev(3))
407 vprev(1)=vnext(1)
408 vprev(2)=vnext(2)
409 vprev(3)=vnext(3)
410 ENDDO
411 forc(1,nx)=-fx*vprev(1)
412 forc(2,nx)=-fx*vprev(2)
413 forc(3,nx)=-fx*vprev(3)
414C-------
415 DO k=1,nx-1
416 dfs=uvarn(mb2+k-1)
417 vprev(1)=dfs*torq(1,k)
418 vprev(2)=dfs*torq(2,k)
419 vprev(3)=dfs*torq(3,k)
420 forc(1,k) =forc(1,k) +vprev(1)
421 forc(2,k) =forc(2,k) +vprev(2)
422 forc(3,k) =forc(3,k) +vprev(3)
423 forc(1,k+1)=forc(1,k+1)-vprev(1)
424 forc(2,k+1)=forc(2,k+1)-vprev(2)
425 forc(3,k+1)=forc(3,k+1)-vprev(3)
426 ENDDO
427C-------
428 DO k=1,nx
429 torq(1,k)=zero
430 torq(2,k)=zero
431 torq(3,k)=zero
432 ENDDO
433C-----------------------------------------------
434C SAVE ELONGATION AND INTERNAL ENERGY OF EACH STRAND (for TH).
435C-----------------------------------------------
436100 CONTINUE
437 fxold=uvar(nb3)
438 DO k=1,nx-1
439 uvarn(mb4+k-1)=mass(k)-uvarn(mb3+k-1)
440 uvarn(mb5+k-1)=uvarn(mb5+k-1)
441 . + half*(uvarn(mb4+k-1)-xiner(k))
442 . *(fx+fxold+uvarn(mb2+k-1)+viscr(k))
443 ENDDO
444C return internal energy of element.
445 keint=1
446 eint =zero
447 DO k=1,nx-1
448 eint=eint+uvarn(mb5+k-1)
449 ENDDO
450 uvar(nb3)=fx
451C-------
452 IF (off==zero) THEN
453C vanish forces into strands.
454 DO k=1,nx-1
455 uvarn(mb2+k-1)=zero
456 ENDDO
457 GOTO 200
458 ENDIF
459C-----------------------------------------------
460C ELEMENT TIME STEP == MINIMUM TIME STEP OF ALL STRANDS.
461C-----------------------------------------------
462 rho = get_u_geo(3,iprop)
463 dte = ep20
464 DO k=1,nx-1
465C-------
466C CHECK NSTRAND STIFNESS AND DAMPING..
467 IF(mass(k)<=em15)THEN
468#include "lockon.inc"
469 WRITE(iout,*)
470 .' ** ERROR NSTRAND : NULL STRAND LENGTH, ELEMENT=',
471 . uid
472#include "lockoff.inc"
473 GOTO 999
474 ENDIF
475C-------
476 xm = rho*uvarn(mb3+k-1)
477C rigidite tangente
478 xkm =stif*(nx-1)/max(em15,l)
479 xcm = (f*dgdx+xc)*l/(max(em15,mass(k))*l0)
480C-------
481C CHECK NSTRAND STIFNESS AND DAMPING..
482 IF(xkm<=em15.AND.xcm<=em15)THEN
483#include "lockon.inc"
484 WRITE(iout,*)
485 .' ** ERROR NSTRAND : NULL STRAND STIFFNESS & DAMPING, ELEMENT=',
486 . uid
487#include "lockoff.inc"
488 GOTO 999
489 ENDIF
490C-------
491C DTK=(SQRT(0.25*XCM*XCM+XM*XKM)-0.5*XCM)/MAX(EM15,XKM)
492C DTC=0.5 * 2.*XM/MAX(EM15,XCM)
493 dtk=(sqrt(xcm*xcm+xm*xkm)-xcm)/max(em15,xkm)
494 dtc=xm/max(em15,xcm)
495 IF (dtk==zero) THEN
496C XKM==0.
497 dtk=dtc
498 ELSE
499 dtk=min(dtk,dtc)
500 ENDIF
501 dte=min(dte,dtk)
502 ENDDO
503C-----------------------------------------------
504C FOR NODAL TIME STEP COMPUTATION :
505C-----------------------------------------------
506200 CONTINUE
507 DO k=1,nx
508 xiner(k) =zero
509 stifr(k) =zero
510 viscr(k) =zero
511 ENDDO
512C
513 IF (off==zero) THEN
514 DO k=1,nx
515 viscm(k) =zero
516 stifm(k) =zero
517 ENDDO
518 ELSE
519 viscm(1) =(f*dgdx+xc)*l/(max(em15,mass(1))*l0)
520 stifm(1) =stif*(nx-1)/max(em15,l)
521 DO k=2,nx-1
522 fact =one/max(em15,mass(k-1))+one/max(em15,mass(k))
523 viscm(k) =fact*(f*dgdx+xc)*l/l0
524 stifm(k) =2.*stif*(nx-1)/max(em15,l)
525 ENDDO
526 viscm(nx) =(f*dgdx+xc)*l/(max(em15,mass(nx-1))*l0)
527 stifm(nx) =stif*(nx-1)/max(em15,l)
528 ENDIF
529C
530 DO k=1,nx
531 mass(k)=uvarn(mb1+k-1)
532 ENDDO
533C-----------------------------------------------
534 1000 FORMAT(1x,'-- RUPTURE OF NSTRAND ELEMENT NUMBER ',i10)
535 RETURN
536C-----
537 999 CONTINUE
538 CALL ancmsg(msgid=217,anmode=aninfo)
539 CALL arret(2)
my_real
#define my_real
Definition cppsort.cpp:32
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
message_mod
Definition message_mod.F:1249
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87
get_u_pid
integer function get_u_pid(ip)
Definition uaccess.F:626
get_u_pnu
integer function get_u_pnu(ivar, ip, k)
Definition uaccess.F:482
get_u_mid
integer function get_u_mid(im)
Definition uaccess.F:668
get_u_mnu
integer function get_u_mnu(ivar, im, k)
Definition uaccess.F:565