xforc28.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    xforc28       ../engine/source/elements/xelem/xforc28.F
!||--- called by ------------------------------------------------------
!||    xforc3        ../engine/source/elements/xelem/xforc3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg        ../engine/source/output/message/message.F
!||    arret         ../engine/source/system/arret.F
!||    get_u_func    ../engine/source/user_interface/ufunc.f
!||    get_u_mid     ../engine/source/user_interface/upidmid.F
!||    get_u_mnu     ../engine/source/user_interface/upidmid.F
!||    get_u_pid     ../engine/source/user_interface/upidmid.F
!||    get_u_pnu     ../engine/source/user_interface/upidmid.F
!||--- uses       -----------------------------------------------------
!||    message_mod   ../engine/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE xforc28(NX     ,
     2             XEL     ,VEL     ,VREL   ,UIX    ,UID    ,
     3             IOUT    ,IPROP   ,IMAT   ,OFF    ,KEINT  ,
     4             EINT   ,MASS    ,XINER   ,STIFM  ,STIFR  ,
     5             VISCM  ,VISCR   ,FORC    ,TORQ   ,
     6             NUVAR   ,UVAR    ,NUVARN ,UVARN  ,
     7             DT      ,DTE     )
C-------------------------------------------------------------------------
C     This subroutine computes multi-purpose element forces and moments
C     when element uses property TYPE28==NSTRAND.
C-------------------------------------------------------------------------
C----------+---------+---+---+--------------------------------------------
C VAR      | SIZE    |TYP| RW| DEFINITION
C----------+---------+---+---+--------------------------------------------
C NX       |    1    | I | R | NUMBER OF NODES (CONSTANT IN THE GROUP)
C----------+---------+---+---+--------------------------------------------
C XEL      | 3*NX    | F | R | NODES COORDINATES
C VEL      | 3*NX    | F | R | NODES VELOCITIES
C VREL     | 3*NX    | F | R | NODES ROTATIONAL VELOCITIES
C----------+---------+---+---+--------------------------------------------
C UIX      |   NX    | I | R | ELEMENT CONNECTIVITY
C                            | IX(J) (1<=J<=NX) : NODE J USER ID
C UID      |    1    | I | R | ELEMENT USER IDENTIFIER
C----------+---------+---+---+--------------------------------------------
C IOUT     |  1      | I | R | OUTPUT FILE UNIT (L00 file)
C IPROP    |  1      | I | R | PROPERTY NUMBER
C IMAT     |  1      | I | R | MATERIAL NUMBER
C----------+---------+---+---+--------------------------------------------
C OFF      |  1      | F | W | ELEMENT STATE (ON=1./OFF=0.)
C KEINT    |  1      | I | W | ELEMENT INTERNAL ENERGY FLAG
C                        | 0 | HAS TO BE COMPUTED BY RADIOSS OUT OF USER ROUTINE.
C                        | 1 | IS RETURNED BY THIS USER ROUTINE.
C EINT     |  1      | F | W | ELEMENT INTERNAL ENERGY IF KEINT=1
C----------+---------+---+---+--------------------------------------------
C MASS     |   NX    | F | W | NODAL MASS
C XINER    |   NX    | F | W | NODAL INERTIA (SPHERICAL)
C STIFM    |   NX    | F | W | NODAL STIFNESS (TIME STEP)
C STIFR    |   NX    | F | W | NODAL ROTATION STIFNESS (TIME STEP)
C VISCM    |   NX    | F | W | NODAL VISCOSITY (TIME STEP)
C VISCR    |   NX    | F | W | NODAL ROTATION VISCOSITY (TIME STEP)
C----------+---------+---+---+--------------------------------------------
C FORC     | 3*NX    | F | W | NODAL FORCES
C TORQ     | 3*NX    | F | W | NODAL TORQUS
C----------+---------+---+---+--------------------------------------------
C NUVAR    |  1      | I | R | NUMBER OF USER ELEMENT VARIABLES
C UVAR     |NUVAR    | F |R/W| USER ELEMENT VARIABLES 
C                            |               (FIX SIZE ZONE)
C NUVARN   |  1      | I | R | NUMBER OF USER ELEMENT VARIABLES PER NODE
C UVARN    |NUVARN*NX| F |R/W| USER ELEMENT VARIABLES PER NODE 
C                            |               (NX DEPENDENT SIZE ZONE)
C----------+---------+---+---+--------------------------------------------
C DT       |  1      | F | R | PREVIOUS TIME INCREMENT :
C                            |    CURRENT TIME POINT IS T = TOLD+DT
C----------+---------+---+---+--------------------------------------------
C DTE      |  1      | F | W | ELEMENT TIME STEP
C----------+---------+---+---+--------------------------------------------
C-------------------------------------------------------------------------
C FUNCTION 
C-------------------------------------------------------------------------
C INTEGER II = GET_U_PNU(I,IP,KK)
C         IFUNCI = GET_U_PNU(I,IP,KFUNC)
C         IPROPI = GET_U_PNU(I,IP,KFUNC)
C         IMATI  = GET_U_PNU(I,IP,KMAT)
C         I     :     VARIABLE INDEX(1 for first variable,...)
C         IP    :     PROPERTY NUMBER
C         KK    :     PARAMETER KFUNC,KMAT,KPROP
C         THIS FUNCTION RETURN THE USER STORED FUNCTION(IF KK=KFUNC), 
C         MATERIAL(IF KK=KMAT) OR PROPERTY(IF KK=KPROP) NUMBER. 
C         SEE LECG28 FOR CORRESPONDING ID STORAGE.
C-------------------------------------------------------------------------
C INTEGER IFUNCI = GET_U_MNU(I,IM,KFUNC)
C         I     :     VARIABLE INDEX(1 for first function)
C         IM    :     MATERIAL NUMBER
C         KFUNC :     ONLY FUNCTION ARE YET AVAILABLE.
C         THIS FUNCTION RETURN THE USER STORED FUNCTION NUMBER(function 
C         referred by users materials).
C         SEE LECM28 FOR CORRESPONDING ID STORAGE.
C-------------------------------------------------------------------------
C my_real PARAMI = GET_U_GEO(I,IP)
C         I     :     PARAMETER INDEX(1 for first parameter,...)
C         IP    :     PROPERTY NUMBER
C         THIS FUNCTION RETURN THE USER GEOMETRY PARAMETERS 
C         NOTE: IF(IP==IPROP) UPARAG(I) == GET_U_GEO(I,IPROP)
C-------------------------------------------------------------------------
C my_real PARAMI = GET_U_MAT(I,IM)
C         I     :     PARAMETER INDEX(1 for first parameter,...)
C         IM    :     MATERIAL NUMBER
C         THIS FUNCTION RETURN THE USER MATERIAL PARAMETERS 
C         NOTE: GET_U_MAT(0,IMAT) RETURN THE DENSITY
C-------------------------------------------------------------------------
C INTEGER PID = GET_U_PID(IP)
C         IP    :     PROPERTY NUMBER
C         THIS FUNCTION RETURN THE USER PROPERTY ID CORRESPONDING TO
C         USER PROPERTY NUMBER IP. 
C-------------------------------------------------------------------------
C INTEGER MID = GET_U_MID(IM)
C         IM   :     MATERIAL NUMBER
C         THIS FUNCTION RETURN THE USER MATERIAL ID CORRESPONDING TO
C         USER MATERIAL NUMBER IM. 
C-------------------------------------------------------------------------
C my_real Y = GET_U_FUNC(IFUNC,X,DXDY)
C         IFUNC :     function number obtained by 
C                     IFUNC = GET_U_MNU(I,IM,KFUNC) or IFUNC = GET_U_PNU(I,IP,KFUNC)
C         X     :     X value
C         DXDY  :     slope dX/dY
C         THIS FUNCTION RETURN Y(X)
C-------------------------------------------------------------------------
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C----------------------------------------------------------
C   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
C----------------------------------------------------------
      INTEGER IOUT,NUVAR,NUVARN,IPROP,IMAT,
     .        NX ,UIX(NX) ,UID, KEINT,
     .        GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
     .        KFUNC,KMAT,KPROP
      my_real
     .   OFF, EINT, XEL(3,NX),  VEL(3,NX) ,VREL(3,NX),
     .   mass(nx) ,  xiner(nx) ,stifm(nx) ,
     .   stifr(nx),  viscm(nx) ,viscr(nx) ,
     .   forc(3,nx), torq(3,nx),
     .   uvar(nuvar),uvarn(nuvarn*nx),dt ,dte ,alpha,
     .   get_u_mat, get_u_geo, get_u_func
      EXTERNAL get_u_mnu,get_u_pnu,get_u_mid,get_u_pid,
     .         get_u_mat,get_u_geo, get_u_func
      parameter(kfunc=29)
      parameter(kmat=31)
      parameter(kprop=33)
C=======================================================================
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NB1, NB2, NB3, MB1, MB2, MB3, MB4, MB5,
     .        k, ifunct, ifv, nit, echange
      my_real
     .  lprev, lnext, vprev(3), vnext(3), vv,xk, xc,
     .  l0, l, dx, dxold, ddx, dt11, dvx, fx, g, y, dfdx, dgdx, f,
     .  stif, mu1, mu2, fric, mu, beta, fmax, dfs, df,   ff,
     .  epsmoy, df1, df2, xn, fxold, epstot, epsmin, epsmax, deps,
     .  rho, xm, xkm, xks, xcm, dtc, dtk, fact, ffac,xfac  
C-----------------------------------------------
C       initial total length UVAR(NB1:NB1)
        nb1=1
C       previous elongation  UVAR(NB2:NB2)
        nb2=nb1+1
C       previous force average UVAR(NB3:NB3)
        nb3=nb2+1
C       initial nodes masses UVARN(MB1:MB1+NX-1)
        mb1=1
C       forces into strands UVARN(MB2:MB2+NX-1)
C       using UVARN(MB2:MB2+NX-2) only.
        mb2=mb1+nx
C       strands initial length UVARN(MB3:MB3+NX-1)
C       using UVARN(MB3:MB3+NX-2) only.
        mb3=mb2+nx
C       strands elongations UVARN(MB4:MB4+NX-1)
C       using UVARN(MB4:MB4+NX-2) only.
        mb4=mb3+nx
C       strands internal energy UVARN(MB5:MB5+NX-1)
C       using UVARN(MB5:MB5+NX-2) only.
        mb5=mb4+nx
C-----------------------------------------------
C       already broken only returns mass & null stiffness.
        IF (off==zero) GOTO 200
C-----------------------------------------------
C       FORCE & PLASTIC DISPLACEMENT,
C-----------------------------------------------
C       ACTUAL LENGTH.
        l=0.0
        DO k=1,nx-1
          lnext=
     .       sqrt((xel(1,k+1)-xel(1,k))*(xel(1,k+1)-xel(1,k))
     .           +(xel(2,k+1)-xel(2,k))*(xel(2,k+1)-xel(2,k))
     .           +(xel(3,k+1)-xel(3,k))*(xel(3,k+1)-xel(3,k)))
          l=l+lnext
          mass(k) =lnext
        ENDDO
C-------
        xk  =get_u_geo(4,iprop)
        ifunct =get_u_pnu(1,iprop,kfunc)
        ifv    =get_u_pnu(2,iprop,kfunc)
        ffac=get_u_geo(12,iprop)
        xfac=get_u_geo(13,iprop)
C-------
C       STIF (tangent stiffness) is considered for friction ..
        IF (ifunct==0.AND.ifv==0) THEN
C------------------------------------------
C        LINEAR ELASTIC
         l0  =uvar(nb1)
         dx  =l-l0
C        raideur tgte par unite de longueur.
         stif=xk
         f   =xk*dx/l0
         epstot=dx/l0
        ELSEIF (ifunct==0.AND.ifv/=0) THEN
C------------------------------------------
C        G(Deps) only.
         l0  =uvar(nb1)
         dx  =l-l0
C        no friction when only viscous ..
         stif=0.
         f = ffac
         epstot=dx/l0
        ELSE
C------------------------------------------
C        NON LINEAR ELASTIC
         l0    =uvar(nb1)
         dx    =l-l0
         epstot=dx/l0
         f     =get_u_func(ifunct,epstot,dfdx)
C        raideur tgte par unite de longueur.
         f = ffac*f
         stif=ffac*dfdx
        ENDIF
C-------------------------------------
C       LINEAR DAMPING.
        dt11 = dt
        IF(dt11==zero)dt11 = ep30
C
        dxold=uvar(nb2)
        ddx=dx-dxold
        dvx=ddx/dt11
C
        deps =(epstot-dxold/l0)/dt11
C
        uvar(nb2)=dx
C
        xc  =get_u_geo(5,iprop)
        dgdx=zero
        IF (ifv>0) THEN
         g  = get_u_func(ifv,deps*xfac,dgdx)
         dgdx = dgdx*xfac
         fx = f*g+xc*deps
         stif=stif*g
        ELSE
         fx= f+xc*deps
         g = one
        ENDIF
C-------
        stif=stif*l/l0
C-----------------------------------------------
C       FRICTION PULLEYS.
C       STRAND STIFFNESS == K*/l0 is assumed to be K*/l * L/L0
C-----------------------------------------------
C       TORQ is used for temporary storage (computation optimization).
        DO k=1,nx-1
         vnext(1)=xel(1,k+1)-xel(1,k)
         vnext(2)=xel(2,k+1)-xel(2,k)
         vnext(3)=xel(3,k+1)-xel(3,k)
         vv=1./max(em15,mass(k))
         vnext(1)=vv*vnext(1)
         vnext(2)=vv*vnext(2)
         vnext(3)=vv*vnext(3)
         torq(1,k)=vnext(1)
         torq(2,k)=vnext(2)
         torq(3,k)=vnext(3)
        ENDDO
C-------
C       temporary storage for internal energy computation.
        DO k=1,nx-1
         viscr(k)=uvarn(mb2+k-1)
         xiner(k)=uvarn(mb4+k-1)
        ENDDO
C-------------------------------------
C       rupture criteria is based upon average of total deformation
C       into nstrand element.
        epsmin=get_u_geo(8,iprop)
        epsmax=get_u_geo(9,iprop)
        IF (epstot<epsmin.OR.epstot>epsmax) THEN
         off=zero
         fx =zero
#include "lockon.inc"
          WRITE(iout, 1000) uid
#include "lockoff.inc"
         GOTO 100
        ENDIF                 
C-------
        mu1    = get_u_geo(10,iprop)
        mu2    = get_u_geo(11,iprop)
C       also available at time 0.
        epsmoy=(dx-dxold)/max(em15,l)
C-------
        xn=nx
        DO k=1,nx-1
         vprev(1)=torq(1,k)
         vprev(2)=torq(2,k)
         vprev(3)=torq(3,k)
         ddx = vprev(1) * (vel(1,k+1) - vel(1,k))
     .       + vprev(2) * (vel(2,k+1) - vel(2,k))
     .       + vprev(3) * (vel(3,k+1) - vel(3,k))
         uvarn(mb2+k-1)=uvarn(mb2+k-1)
C tous les brins ont la meme longueur correspond a :
     .       +stif*(dt*ddx*(xn-one)/max(em15,l)-epsmoy)
        ENDDO
C-------
C
        echange=1
        nit=0
        DO WHILE(echange==1.AND.nit<10)
          echange=0
          nit=nit+1
          DO k=1,nx-1
C           stifr utilise de maniere temporaire.
            stifr(k)=uvarn(mb2+k-1)
          END DO
          DO k=2,nx-1
C          default pulley friction.
           fric=mu1
           mu  = get_u_geo(50+k,iprop)
           IF (mu>=zero) fric=mu
           mu  = get_u_geo(525+k-1,iprop)
           IF (mu>=zero) THEN 
             fric=fric + half*mu
           ELSE
             fric=fric + half*mu2
           ENDIF
           mu  = get_u_geo(525+k,iprop)
           IF (mu>=zero) THEN 
             fric=fric+half*mu
           ELSE
             fric=fric + half*mu2
           ENDIF
C
           vprev(1)=torq(1,k-1)
           vprev(2)=torq(2,k-1)
           vprev(3)=torq(3,k-1)
           vnext(1)=-torq(1,k)
           vnext(2)=-torq(2,k)
           vnext(3)=-torq(3,k)
           alpha=vprev(1)*vnext(1)+vprev(2)*vnext(2)+vprev(3)*vnext(3)
           alpha = min(alpha,one)
           alpha = max(alpha,-one)
           beta  = pi-acos(alpha)
C
           ff   = two*fx+stifr(k-1)+stifr(k)
           fmax = max(zero,ff*tanh(half*fric*beta))
           dfs  = stifr(k-1)-stifr(k)
           IF(abs(dfs)>fmax)THEN
             df =sign(abs(dfs)-fmax,dfs)
C            tous les brins ont la meme longueur correspond a :
             uvarn(mb2+k-2)=uvarn(mb2+k-2)-half*df
             uvarn(mb2+k-1)=uvarn(mb2+k-1)+half*df
             echange=1
           END IF
       END DO
       END DO
C-----------------------------------------------
C       RETURN FORCES.
C-----------------------------------------------
        vprev(1)=torq(1,1)
        vprev(2)=torq(2,1)
        vprev(3)=torq(3,1)
        forc(1,1)=fx*vprev(1)
        forc(2,1)=fx*vprev(2)
        forc(3,1)=fx*vprev(3)
        DO k=2,nx-1
         vnext(1)=torq(1,k)
         vnext(2)=torq(2,k)
         vnext(3)=torq(3,k)
         forc(1,k)=fx*(vnext(1)-vprev(1))
         forc(2,k)=fx*(vnext(2)-vprev(2))
         forc(3,k)=fx*(vnext(3)-vprev(3))
         vprev(1)=vnext(1)
         vprev(2)=vnext(2)
         vprev(3)=vnext(3)
        ENDDO
        forc(1,nx)=-fx*vprev(1)
        forc(2,nx)=-fx*vprev(2)
        forc(3,nx)=-fx*vprev(3)
C-------
        DO k=1,nx-1
         dfs=uvarn(mb2+k-1)
         vprev(1)=dfs*torq(1,k)
         vprev(2)=dfs*torq(2,k)
         vprev(3)=dfs*torq(3,k)
         forc(1,k)  =forc(1,k)  +vprev(1)
         forc(2,k)  =forc(2,k)  +vprev(2)
         forc(3,k)  =forc(3,k)  +vprev(3)
         forc(1,k+1)=forc(1,k+1)-vprev(1)
         forc(2,k+1)=forc(2,k+1)-vprev(2)
         forc(3,k+1)=forc(3,k+1)-vprev(3)
        ENDDO
C-------
        DO k=1,nx
         torq(1,k)=zero
         torq(2,k)=zero
         torq(3,k)=zero 
        ENDDO
C-----------------------------------------------
C       SAVE ELONGATION AND INTERNAL ENERGY OF EACH STRAND (for TH).
C-----------------------------------------------
100     CONTINUE
        fxold=uvar(nb3)
        DO k=1,nx-1
         uvarn(mb4+k-1)=mass(k)-uvarn(mb3+k-1)
         uvarn(mb5+k-1)=uvarn(mb5+k-1)
     .    + half*(uvarn(mb4+k-1)-xiner(k))
     .        *(fx+fxold+uvarn(mb2+k-1)+viscr(k))
        ENDDO
C       return internal energy of element.
        keint=1
        eint =zero
        DO k=1,nx-1
         eint=eint+uvarn(mb5+k-1)
        ENDDO
        uvar(nb3)=fx
C-------
        IF (off==zero) THEN
C        vanish forces into strands.
         DO k=1,nx-1
          uvarn(mb2+k-1)=zero
         ENDDO
         GOTO 200
        ENDIF
C-----------------------------------------------
C       ELEMENT TIME STEP == MINIMUM TIME STEP OF ALL STRANDS.
C-----------------------------------------------
        rho  = get_u_geo(3,iprop)
        dte  = ep20
        DO k=1,nx-1
C-------
C         CHECK NSTRAND STIFNESS AND DAMPING..
          IF(mass(k)<=em15)THEN
#include "lockon.inc"
             WRITE(iout,*)
     .' ** ERROR NSTRAND : NULL STRAND LENGTH, ELEMENT=',
     .        uid
#include "lockoff.inc"
            GOTO 999
          ENDIF
C-------
          xm   = rho*uvarn(mb3+k-1)
C         rigidite tangente
          xkm  =stif*(nx-1)/max(em15,l)
          xcm  = (f*dgdx+xc)*l/(max(em15,mass(k))*l0)
C-------
C         CHECK NSTRAND STIFNESS AND DAMPING..
          IF(xkm<=em15.AND.xcm<=em15)THEN
#include "lockon.inc"
             WRITE(iout,*)
     .' ** ERROR NSTRAND : NULL STRAND STIFFNESS & DAMPING, ELEMENT=',
     .        uid
#include "lockoff.inc"
            GOTO 999
          ENDIF
C-------
C           DTK=(SQRT(0.25*XCM*XCM+XM*XKM)-0.5*XCM)/MAX(EM15,XKM)
C           DTC=0.5 * 2.*XM/MAX(EM15,XCM)
          dtk=(sqrt(xcm*xcm+xm*xkm)-xcm)/max(em15,xkm)
          dtc=xm/max(em15,xcm)
          IF (dtk==zero) THEN
C          XKM==0.
           dtk=dtc
          ELSE
           dtk=min(dtk,dtc)
          ENDIF
          dte=min(dte,dtk)
        ENDDO
C-----------------------------------------------
C       FOR NODAL TIME STEP COMPUTATION :
C-----------------------------------------------
200     CONTINUE
        DO k=1,nx
         xiner(k) =zero
         stifr(k) =zero
         viscr(k) =zero 
        ENDDO
C
        IF (off==zero) THEN
         DO k=1,nx
          viscm(k) =zero
          stifm(k) =zero  
         ENDDO
        ELSE
         viscm(1) =(f*dgdx+xc)*l/(max(em15,mass(1))*l0)
         stifm(1) =stif*(nx-1)/max(em15,l)
         DO k=2,nx-1
          fact     =one/max(em15,mass(k-1))+one/max(em15,mass(k))
          viscm(k) =fact*(f*dgdx+xc)*l/l0
          stifm(k) =2.*stif*(nx-1)/max(em15,l)
         ENDDO
         viscm(nx) =(f*dgdx+xc)*l/(max(em15,mass(nx-1))*l0)
         stifm(nx) =stif*(nx-1)/max(em15,l)
        ENDIF
C
        DO k=1,nx
         mass(k)=uvarn(mb1+k-1)
        ENDDO
C-----------------------------------------------
 1000 FORMAT(1x,'-- RUPTURE OF NSTRAND ELEMENT NUMBER ',i10)
      RETURN
C-----
 999  CONTINUE
      CALL ancmsg(msgid=217,anmode=aninfo)
      CALL arret(2)
      END