diff --git a/Components/Elements/Extended-Hill-Type-Muscle-1D-user-spring b/Components/Elements/Extended-Hill-Type-Muscle-1D-user-spring new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/Components/Elements/Extended-Hill-Type-Muscle-1D-user-spring @@ -0,0 +1 @@ + diff --git a/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Code/lecg30.F b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Code/lecg30.F new file mode 100644 index 0000000..52e041f --- /dev/null +++ b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Code/lecg30.F @@ -0,0 +1,524 @@ +C MIT License +C +C Copyright (c) 2025 RM412 +C +C Permission is hereby granted, free of charge, to any person obtaining a copy +C of this software and associated documentation files (the "Software"), to deal +C in the Software without restriction, including without limitation the rights +C to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +C copies of the Software, and to permit persons to whom the Software is +C furnished to do so, subject to the following conditions: +C +C The above copyright notice and this permission notice shall be included in all +C copies or substantial portions of the Software. +C +C THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +C IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +C FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +C AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +C LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +C OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +C SOFTWARE. +C +C Suggestion : +C +C BSD 3-Clause License +C Copyright (c) 2014 +C D. Haeufle, M. Günther, A. Bayer, S. Schmitt +C Copyright (c) 2017 +C C. Kleinbach, O. Martynenko, J. Promies, D. F. B. Haeufle, J. Fehr, S. Schmitt +C Copyright (c) 2019 +C O. Martynenko, F. Kempter, C. Kleinbach, S. Schmitt, J. Fehr +C Copyright (c) 2020 +C P. Lerge, O. Martynenko, L. Nölle, S. Schmitt, I. Wochner +C Copyright (c) 2022 +C O. Martynenko, F. Kempter, C. Kleinbach, L. Nölle, P. Lerge, S. Schmitt, J. Fehr +C Copyright (c) 2026 +C A. Ardisson, E. Wagnac, M.-H. Beauséjour/ÉTS +C Adaptation and port to OpenRadioss +C All rights reserved. +C Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following C C conditions are met: +C Redistributions of source code must retain the above copyright notice, +C this list of conditions and the following disclaimer. +C +C +C Redistributions in binary form must reproduce the above copyright notice, +C this list of conditions and the following disclaimer in the documentation +C and/or other materials provided with the distribution. +C +C +C Neither the name of the copyright holder nor the names of its contributors +C may be used to endorse or promote products derived from this software +C without specific prior written permission. +C +C +C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” +C AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +C IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +C ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +C LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +C POSSIBILITY OF SUCH DAMAGE. + +Chd|==================================================================== +Chd| LECG30 +Chd|-- called by ----------- +Chd|-- calls --------------- +Chd| SET_U_GEO +Chd| SET_U_PNU +Chd|==================================================================== + SUBROUTINE LECG30(IIN ,IOUT ,NUVAR ,PARGEO) +C----------------------------------------------- +C I m p l i c i t T y p e s +C----------------------------------------------- + IMPLICIT NONE +C----------+---------+---+---+-------------------------------------------- +C VAR | SIZE |TYP| RW| DEFINITION +C----------+---------+---+---+-------------------------------------------- +C IIN | 1 | I | R | INPUT FILE UNIT (D00 file) +C IOUT | 1 | I | R | OUTPUT FILE UNIT (L00 file) +C NUVAR | 1 | I | W | NUMBER OF USER ELEMENT VARIABLES +C----------+---------+---+---+-------------------------------------------- +C PARGEO | * | F | W | 1)SKEW NUMBER +C | | | | 2)STIFNESS FOR INTERFACE +C | | | | 3)FRONT WAVE OPTION +C | | | | 4)... not yet used +C----------+---------+---+---+-------------------------------------------- +C +C This subroutine read the user geometry parameters. +C +C The geometry datas has to bee stored in radioss storage +C with the function SET_U_GEO(value_index,value). +C +C If some standard radioss functions (time function or +C x,y function) are used, this function IDs has to +C bee stored with the function SET_U_PNU(func_index,func_id,KFUNC). +C +C If this property refers to a user material, this +C material IDs has to bee stored with the function +C SET_U_PNU(mat_index,mat_id,KMAT). +C +C If this property refers to a user property, this +C sub-property IDs has to bee stored with the function +C SET_U_PNU(sub_prop_index,sub_prop_id,KMAT). +C +C SET_U_GEO and SET_U_PNU return 0 if no error +C SET_U_GEO and SET_U_PNU return the maximum allowed index +C if index is larger than this maximum +C----------------------------------------------- +C D u m m y A r g u m e n t s +C----------------------------------------------- + INTEGER IIN,IOUT,NUVAR + DOUBLE PRECISION + . PARGEO(*) + INTEGER SET_U_PNU,SET_U_GEO, + . KFUNC,KMAT,KPROP + EXTERNAL SET_U_PNU,SET_U_GEO + PARAMETER (KFUNC=29) + PARAMETER (KMAT=31) + PARAMETER (KPROP=30) +C================================================================= +C----------------------------------------------- +C L o c a l V a r i a b l e s +C----------------------------------------------- + INTEGER SENS_ID,IERROR,ActOpt,output_MTD,IFUNC1 + INTEGER IFUNC2,IFUNC3 + DOUBLE PRECISION + . AMAS,DIAMETER,ELASTIF,XK,SENS_IDf,q0,tauq,betaq,k,m,muscle, + . Fmax,t_lCEopt,dWdes, + . nuCEdes,dWasc,nuCEasc,Arel0,Brel0, + . Secc,Fecc,LPEE0,nuPEE,FPEE,dUSEEnll,duSEEl, + . dFSEE0,Damping,Param1,Param2, + . timeste,Activation,target_l_CE,kp,kd, + . dotlCEdelay,t_PreSim,threshold,ActOpt_f,SENS_ID_f,output_MTD_f, + . STIM_ID +C +C================================================================= +C + WRITE(IOUT,*) ' MUSCLE User routine : Debug Start USER ' + WRITE(IOUT,*) ' by C. Kleinbach, O. Martynenko, J. Promies, D.F.B. Haeufle, J. Fehr, S. Schmitt' + WRITE(IOUT,*) ' Extended Hill-type Muscle Model ' + WRITE(IOUT,*) ' *** Springs *** ' + WRITE(IOUT,*) ' M.Bulla, A.Ardisson 2025.03.10 ' + WRITE(IOUT,*) ' ' + + NUVAR = 39 + +c! hsv(1) to hsv(39) +c! ------------------------------------------------------------ +c! ------ history variables - overview ------------------------ +c! ------------------------------------------------------------ +c! hsv(1) = sig(1) +c! hsv(2) = STIM +c! hsv(3) = q +c! hsv(4) = F_MTC +c! hsv(5) = F_CE +c! hsv(6) = F_PEE +c! hsv(7) = F_SEE +c! hsv(8) = F_SDE +c! hsv(9) = elleng (l_MTC) +c! hsv(10) = l_CE +c! hsv(11) = dot_l_MTC +c! hsv(12) = dot_l_CE +c! hsv(13) = counter_output +c! hsv(14) = F_isom +c! hsv(15) = gam_rel +c! hsv(16) = l_MTC_0 (initial element length) +c! ------------------------------------------------------------ +c! ------ only needed if controller are used ------------------ +c! ------------------------------------------------------------ +c! hsv(20) = lCEdelay +c! hsv(21) = dotlCEdelay +c! hsv(22) = l_CE_ref (for reflex controller) +c! hsv(23) = strain (for reflex controller) +c! hsv(24) = STIM_reflex_prev (STIM value of reflex controller in previous timestep) +c! ------------------------------------------------------------ +c! ------ delay buffer----------------------------------------- +c! ------------------------------------------------------------ +c! hsv(30) = buffersize +c! hsv(31) = idx_begin_lCEbuffer (hsv(142:145) seems to be used internally by lsdyna, suggestion: use hsv(31)=150) +c! hsv(32) = indexr1...index1 of lce-ringbuffer (not eq index of hsv) +c! hsv(33) = indexr1...index2 of lce-ringbuffer (not eq index of hsv) +c! hsv(34) = lasttime +c! hsv(36) = begindotlCE +c! hsv(37) = indexdotr1...index1 of dotlce-ringbuffer (not eq index of hsv) +c! hsv(38) = indexdotr2...index2 of dotlce-ringbuffer (not eq index of hsv) +c! hsv(39) = dotlasttime +c! hsv(hsv(31):hsv(31)+hsv(30)-1) = ringbuffer_l_CE +c! hsv(hsv(36):hsv(36)+hsv(30)-1) = ringbuffer_dot_l_CE +C +C +c!*MAT_USER_DEFINED_MATERIAL_MODELS +c!$ Hill-type Muscle for a cat soleus moerl12 +c!$# mid ro mt lmc nhv iortho ibulk ig +c! 1 1.000E-6 41 28 14 0 3 4 +c!$# ivect ifail itherm ihyper ieos lmca unused unused +c! 0 0 0 0 0 0 +c!$# cm(1) cm(2) cm(3) cm(4) cm(5) cm(6) cm(7) cm(8) +c!$# ActOpt STIM_ID q0 tauq/c betaq/eta k m muscle +c! 2 3 1.0E-4 1.373E-4 5.27E4 2.9 22.54 +c!$# cm(9) cm(10) cm(11) cm(12) cm(13) cm(14) cm(15) cm(16) +c!$# Fmax lCEopt dWdes nuCEdes dWasc nuCEasc Arel0 Brel0 +c! 10.0 0.053 0.35 1.5 0.35 3.0 0.07 0.2 +c!$# cm(17) cm(18) cm(19) cm(20) cm(21) cm(22) cm(23) cm(24) +c!$# Secc Fecc LPEE0 nuPEE FPEE lSEE0 dUSEEnll duSEEl +c! 2.0 1.5 0.9 2.5 2.0 0.060 0.0425 0.017 +c!$# cm(25) cm(26) cm(27) cm(28) cm(29) cm(30) cm(31) cm(32) +c!$# dFSEE0 Damping Param1 Param2 output me o.timeste +c! 4.0 3.0 0.3 0.01 +c! Variables (activation dynamic) +c! cm(1)=Activation Option (EQ.0. Activation Values see STIM ID, EQ.1 Zajac, EQ.2 Hatze) +c! cm(2)=GT0: LE.0.0 constant values for STIM or Activation GT.0.0 curve id for STIM or Activation +c! cm(3)=q0 minimum value of q +c! cm(4)=tau_q time constant of rising activation LT0 curve id of tau_q over time +c! cm(4)=c +c! cm(5)=beta_q ratio between tau_q and time constant of falling activation +c! cm(5)=eta +c! cm(6)=k +c! cm(7)=m +c! cm(8)=muscle length offset (not needed, use instead *Part_Averaged for routing) +c! Variables (activation dynamic) +c! +c! Variables (isometric force) +c! +c! cm(9)=F_max maximum isometric force +c! cm(10)=l_CE_opt optimal fibre length +c! cm(11)=dW_des +c! cm(12)=nu_CE_des +c! cm(13)=dW_asc +c! cm(14)=nu_CE_asc +c! +c! Variables (Hill Parameter, concentric) +c! +c! cm(15)=A_rel_0 +c! cm(16)=B_rel_0 +c! +c! Variables (van Soest Parameter, eccentric) +c! +c! cm(17)=S_ecc +c! cm(18)=S_ecc +c! +c! Variables (Parallel elastic element) +c! +c! cm(19)=L_PEE_0 +c! cm(20)=nu_PEE +c! cm(21)=F_PEE +c! +c! Variables (Seriell elastic element) +c! +c! cm(22)=l_SEE_0 +c! cm(23)=dU_SEE_nllfindkParams +c! cm(24)=dU_SEE_l +c! cm(25)=dF_SEE_0 +c! +c! Variables (Damping element) +c! +c! cm(26)= --- (former damping method) +c! cm(27)=D_SE +c! cm(28)=R_SE +c! +c! Variables (Output definition; musout.(partid)) +c! +c! cm(29)=output method (EQ.0. no output +c! EQ.1. basic output (idpart, tt, hsv(2:10)) +c! EQ.2. advanced output (basic output plus dot_l_CE, dot_l_MTC, lCEdelay, dotlCEdelay)) +c! +c! cm(30)=timestep of outputfile +c! +c! Variables for Controller +c! +c! cm(33)=Activation Method (EQ.1. lambda_controller +c! EQ.2. hybrid_controller +c! EQ.3. reflexive controller) +c! +c! cm(34)=target l_CE !!! POSSIBLE CURVE/FUNCTION +c! cm(35)=kp +c! cm(36)=kd +c! cm(37)=delay of lCEdelay / dotlCEdelay +c! cm(38)=time till swap from alpha to lambda / t_PreSim for reflexive controller +c! cm(39)=threshold for reflex controller (e.g. 0.10 for a 10% strain threshold) +c! +c! [...] +c! +c! d1 - strain rate/increment in x direction, local x for shells +c! d2 - strain rate/increment in y direction, local y for shells +c! d3 - strain rate/increment in z direction, local z for shells +c! ------------------------------------------------------------ +c! ------ history variables - overview ------------------------ +c! ------------------------------------------------------------ +c! hsv(1) = sig(1) +c! hsv(2) = STIM +c! hsv(3) = q +c! hsv(4) = F_MTC +c! hsv(5) = F_CE +c! hsv(6) = F_PEE +c! hsv(7) = F_SEE +c! hsv(8) = F_SDE +c! hsv(9) = elleng (l_MTC) +c! hsv(10) = l_CE +c! hsv(11) = dot_l_MTC +c! hsv(12) = dot_l_CE +c! hsv(13) = counter_output +c! hsv(14) = F_isom +c! hsv(15) = gam_rel +c! hsv(16) = l_MTC_0 (initial element length) +c! ------------------------------------------------------------ +c! ------ only needed if controller are used ------------------ +c! ------------------------------------------------------------ +c! hsv(20) = lCEdelay +c! hsv(21) = dotlCEdelay +c! hsv(22) = l_CE_ref (for reflex controller) +c! hsv(23) = strain (for reflex controller) +c! hsv(24) = STIM_reflex_prev (STIM value of reflex controller in previous timestep) +c! ------------------------------------------------------------ +c! ------ delay buffer----------------------------------------- +c! ------------------------------------------------------------ +c! hsv(30) = buffersize +c! hsv(31) = idx_begin_lCEbuffer (hsv(142:145) seems to be used internally by lsdyna, suggestion: use hsv(31)=150) +c! hsv(32) = indexr1...index1 of lce-ringbuffer (not eq index of hsv) +c! hsv(33) = indexr1...index2 of lce-ringbuffer (not eq index of hsv) +c! hsv(34) = lasttime +c! hsv(36) = begindotlCE +c! hsv(37) = indexdotr1...index1 of dotlce-ringbuffer (not eq index of hsv) +c! hsv(38) = indexdotr2...index2 of dotlce-ringbuffer (not eq index of hsv) +c! hsv(39) = dotlasttime +c! hsv(hsv(31):hsv(31)+hsv(30)-1) = ringbuffer_l_CE +c! hsv(hsv(36):hsv(36)+hsv(30)-1) = ringbuffer_dot_l_CE +c! + READ(IIN,ERR=999,END=999,FMT='(3F20.0,2I10)') + . AMAS,DIAMETER,XK,ActOpt,SENS_ID + IF (SENS_ID .NE.0) SENS_IDf = real(SENS_ID)+0.00001 + + READ(IIN,ERR=999,END=999,FMT='(5F20.0)') + . STIM_ID,q0,tauq,betaq,k + READ(IIN,ERR=999,END=999,FMT='(5F20.0)') + . m,muscle,Fmax,t_lCEopt,dWdes + READ(IIN,ERR=999,END=999,FMT='(5F20.0)') + . nuCEdes,dWasc,nuCEasc,Arel0,Brel0 + READ(IIN,ERR=999,END=999,FMT='(5F20.0)') + . Secc,Fecc,LPEE0,nuPEE,FPEE + READ(IIN,ERR=999,END=999,FMT='(3F20.0)') + . dUSEEnll, duSEEl, dFSEE0 + READ(IIN,ERR=999,END=999,FMT='(3F20.0,I10)') + . Damping,Param1,Param2,output_MTD + READ(IIN,ERR=999,END=999,FMT='(5F20.0)') + . timeste,Activation,target_l_CE,kp,kd + READ(IIN,ERR=999,END=999,FMT='(3F20.0)') + . dotlCEdelay,t_PreSim,threshold +c Variables (isometric force) + + PARGEO(1) = 0 + PARGEO(2) = XK + PARGEO(3) = 1 +C +c! WRITE(IOUT,1000) +c! . AMAS,ELASTIF,XK,SENS_ID +c! WRITE(IOUT,1000) +c! . AMAS,DIAMETER,XK,ActOpt,STIM_ID,q0,tauq,betaq,k,m,muscle,Fmax,lCEopt,dWdes, +c! . nuCEdes,dWasc,nuCEasc,Arel0,Brel0,Secc,Fecc,LPEE0,nuPEE,FPEE, +c! . lSEE0,dUSEEnll,duSEEl,dFSEE0,Damping,Param1,Param2,output_MTD, +c! . timeste,Activation,target_l_CE,kp,kd,dotlCEdelay,t_PreSim,threshold + WRITE(IOUT,1000) + . AMAS,DIAMETER,XK,ActOpt,STIM_ID,q0,tauq,betaq,k,m,muscle,Fmax,t_lCEopt,dWdes, + . nuCEdes,dWasc,nuCEasc,Arel0,Brel0,Secc,Fecc,LPEE0,nuPEE,FPEE, + . dUSEEnll,duSEEl,dFSEE0,Damping,Param1,Param2,output_MTD, + . timeste,Activation,target_l_CE,kp,kd,dotlCEdelay,t_PreSim,threshold +c!C + IERROR = SET_U_GEO(1,AMAS) +c! IERROR = SET_U_GEO(2,ELASTIF) + IERROR = SET_U_GEO(2,XK) + IERROR = SET_U_GEO(3,DIAMETER) + + ActOpt_f= REAL(ActOpt)+0.00001 + print *,'ActOpt = ',ActOpt + print *,'REAL(ActOpt) = ',REAL(ActOpt) +c! SENS_ID_f = REAL(SENS_ID)+0.00001 + IERROR = SET_U_GEO(4, ActOpt_f) + SENS_ID_f = DBLE(SENS_ID)+0.00001 + SENS_ID_f = SENS_ID+0.00001 + IERROR = SET_U_GEO(5,SENS_ID_f) +c! IERROR = SET_U_GEO(5,REAL(SENS_ID)+0.00001) + IERROR = SET_U_GEO(6,STIM_ID) + print*, STIM_ID +C----------------------------------------------------------------------- + IF (STIM_ID .GT. 0.0) THEN + IFUNC1 = INT(STIM_ID+0.00001) + print*, 'IFUNC1_lecg = ', IFUNC1 + IERROR = SET_U_PNU(1,IFUNC1,KFUNC) + ENDIF +C----------------------------------------------------------------------- + IERROR = SET_U_GEO(7,q0) + IERROR = SET_U_GEO(8,tauq) +C----------------------------------------------------------------------- + IF (tauq .LT. 0.0) THEN !not sure of this + IFUNC2 = INT(tauq+0.00001) + print*, 'IFUNC2_lecg', IFUNC2 + IERROR = SET_U_PNU(2,IFUNC2,KFUNC) + ENDIF +C----------------------------------------------------------------------- + IERROR = SET_U_GEO(9,betaq) + IERROR = SET_U_GEO(10,k) + + IERROR = SET_U_GEO(11,m) + IERROR = SET_U_GEO(12,muscle) + IERROR = SET_U_GEO(13,Fmax) + IERROR = SET_U_GEO(14,t_lCEopt) + IERROR = SET_U_GEO(15,dWdes) + + IERROR = SET_U_GEO(16,nuCEdes) + IERROR = SET_U_GEO(17,dWasc) + IERROR = SET_U_GEO(18,nuCEasc) + IERROR = SET_U_GEO(19,Arel0) + IERROR = SET_U_GEO(20,Brel0) + + IERROR = SET_U_GEO(21,Secc) + IERROR = SET_U_GEO(22,Fecc) + IERROR = SET_U_GEO(23,LPEE0) + IERROR = SET_U_GEO(24,nuPEE) + IERROR = SET_U_GEO(25,FPEE) + + IERROR = SET_U_GEO(27,dUSEEnll) + IERROR = SET_U_GEO(28,duSEEl) + IERROR = SET_U_GEO(29,dFSEE0) + + IERROR = SET_U_GEO(30,Damping) + IERROR = SET_U_GEO(31,Param1) + IERROR = SET_U_GEO(32,Param2) + output_MTD_f = REAL(output_MTD)+0.00001 + IERROR = SET_U_GEO(33,output_MTD_f) +c! IERROR = SET_U_GEO(33,REAL(output_MTD)+0.00001) + + IERROR = SET_U_GEO(34,timeste) + IERROR = SET_U_GEO(35,Activation) + IERROR = SET_U_GEO(36,target_l_CE) +C----------------------------------------------------------------------- + !IF (target_l_CE .GT. 0.0) THEN + IFUNC3 = INT(target_l_CE+0.00001) + print*, 'IFUNC3_lecg', IFUNC3 + IERROR = SET_U_PNU(3,IFUNC3,KFUNC) + !ENDIF +C----------------------------------------------------------------------- + IERROR = SET_U_GEO(37,kp) + IERROR = SET_U_GEO(38,kd) + + IERROR = SET_U_GEO(39,dotlCEdelay) + IERROR = SET_U_GEO(40,t_PreSim) + IERROR = SET_U_GEO(41,threshold) + +c! IERROR = SET_U_GEO(4,SENS_IDf) + +C + RETURN + 999 CONTINUE + WRITE(IOUT,*)' **ERROR IN USER PROPERTY INPUT !!!' + RETURN + 1000 FORMAT( + & 5X,' Extended Hill-type Muscle Model:',/, + & 5X,' Version 25.06.2024 a ',/, + & 5X,'',/, + & 5X,'Density. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'Crossection area . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'STIFFNESS FOR INTERFACE. . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'Variables (activation dynamic):',/, + & 5X,'Activation Option. . . . . . . . . . . . . . . . . . . . . . . . . . . .=',I10/, + & 5X,'(EQ.0. Activation Values see STIM ID, EQ.1 Zajac, EQ.2 Hatze)',/, + & 5X,'STIM_ID . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'LE.0.0 constant values for STIM or Activation GT.0.0 curve id for STIM or Activation',/, + & 5X,'q0 minimum value of q . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'tau_q time constant of rising activation LT0 curve id of tau_q over time=',E12.4/, + & 5X,'beta_q ratio between tau_q and time constant of falling activation eta .=',E12.4/, + & 5X,'k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'muscle length offset (not needed, use instead *Part_Averaged for routing)',E12.4/, + & 5X,'Variables (isometric force):',/, + & 5X,'F_max maximum isometric force. . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'t_l_CE_opt optimal fibre length proportion . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'dW_des . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'nu_CE_des. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'dW_asc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'nu_CE_asc. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'Variables (Hill Parameter, concentric):',/, + & 5X,'A_rel_0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'B_rel_0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'Variables (van Soest Parameter, eccentric):',/, + & 5X,'S_ecc. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'F_ecc. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'Variables (Parallel elastic element):',/, + & 5X,'L_PEE_0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'nu_PEE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'F_PEE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'Variables (Seriell elastic element):',/, + & 5X,'l_SEE_0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'dU_SEE_nllfindkParams. . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'dU_SEE_l . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'dF_SEE_0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'Variables (Damping element):',/, + & 5X,'former damping method. . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'D_SE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'R_SE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'Variables (Output definition; musout.(partid)):',/, + & 5X,'output method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',I10/, + & 5X,'EQ.0. no output ',/, + & 5X,'EQ.1. basic output (idpart, tt, hsv(2:10)) ',/, + & 5X,'EQ.2. advanced output (basic output plus dot_l_CE, dot_l_MTC, lCEdelay, dotlCEdelay) ',/, + & 5X,'timestep of outputfile . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'Variables (for Controller):',/, + & 5X,'Activation Method. . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'EQ.1. lambda_controller',/, + & 5X,'EQ.2. hybrid_controller',/, + & 5X,'EQ.3. reflexive controller',/, + & 5X,'target l_CE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'kp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'kd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'delay of lCEdelay / dotlCEdelay. . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'time till swap from alpha to lambda / t_PreSim for reflexive controller.=',E12.4/, + & 5X,'threshold for reflex controller (e.g. 0.10 for a 10% strain threshold) .=',E12.4//) + + +c! & 5X,'STIFFNESS FOR INTERFACE . . . . . . . . .=',E12.4/, +c! & 5X,'SENSOR ID . . . . . . . . . . . . . . . .=',I10//) + END diff --git a/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Code/rini30.F b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Code/rini30.F new file mode 100644 index 0000000..194650a --- /dev/null +++ b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Code/rini30.F @@ -0,0 +1,228 @@ +C MIT License +C +C Copyright (c) 2025 RM412 +C +C Permission is hereby granted, free of charge, to any person obtaining a copy +C of this software and associated documentation files (the "Software"), to deal +C in the Software without restriction, including without limitation the rights +C to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +C copies of the Software, and to permit persons to whom the Software is +C furnished to do so, subject to the following conditions: +C +C The above copyright notice and this permission notice shall be included in all +C copies or substantial portions of the Software. +C +C THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +C IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +C FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +C AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +C LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +C OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +C SOFTWARE. +C +C Suggestion : +C +C BSD 3-Clause License +C Copyright (c) 2014 +C D. Haeufle, M. Günther, A. Bayer, S. Schmitt +C Copyright (c) 2017 +C C. Kleinbach, O. Martynenko, J. Promies, D. F. B. Haeufle, J. Fehr, S. Schmitt +C Copyright (c) 2019 +C O. Martynenko, F. Kempter, C. Kleinbach, S. Schmitt, J. Fehr +C Copyright (c) 2020 +C P. Lerge, O. Martynenko, L. Nölle, S. Schmitt, I. Wochner +C Copyright (c) 2022 +C O. Martynenko, F. Kempter, C. Kleinbach, L. Nölle, P. Lerge, S. Schmitt, J. Fehr +C Copyright (c) 2026 +C A. Ardisson, E. Wagnac, M.-H. Beauséjour/ÉTS +C Adaptation and port to OpenRadioss +C All rights reserved. +C Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following C C conditions are met: +C Redistributions of source code must retain the above copyright notice, +C this list of conditions and the following disclaimer. +C +C +C Redistributions in binary form must reproduce the above copyright notice, +C this list of conditions and the following disclaimer in the documentation +C and/or other materials provided with the distribution. +C +C +C Neither the name of the copyright holder nor the names of its contributors +C may be used to endorse or promote products derived from this software +C without specific prior written permission. +C +C +C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” +C AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +C IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +C ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +C LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +C POSSIBILITY OF SUCH DAMAGE. + +Cgw41n4 source deplace de lecg30.F vers rini30.F +Chd|==================================================================== +Chd| RINI30 +Chd|-- called by ----------- +Chd|-- calls --------------- +Chd|==================================================================== + SUBROUTINE RINI30(NEL ,IOUT ,IPROP , + 3 IX ,XL ,MASS ,XINER ,STIFM , + 4 STIFR ,VISCM ,VISCR ,UVAR ,NUVAR ) +C------------------------------------------------------------------------- +C This subroutine initialize springs using user properties. +C------------------------------------------------------------------------- +C----------+---------+---+---+-------------------------------------------- +C VAR | SIZE |TYP| RW| DEFINITION +C----------+---------+---+---+-------------------------------------------- +C IOUT | 1 | I | R | OUTPUT FILE UNIT (L00 file) +C IPROP | 1 | I | R | PROPERTY NUMBER +C----------+---------+---+---+-------------------------------------------- +C IX | 3*NEL | I | R | SPRING CONNECTIVITY +C | IX(1,I) NODE 1 ID +C | IX(2,I) NODE 2 ID +C | IX(3,I) OPTIONNAL NODE 3 ID +C | IX(4,I) SPRING ID +C XL | NEL | F | R | ELEMENT LENGTH +C----------+---------+---+---+-------------------------------------------- +C MASS | NEL | F | W | ELEMENT MASS +C XINER | NEL | F | W | ELEMENT INERTIA (SPHERICAL) +C STIFM | NEL | F | W | ELEMENT STIFNESS (TIME STEP) +C STIFR | NEL | F | W | ELEMENT ROTATION STIFNESS (TIME STEP) +C VISCM | NEL | F | W | ELEMENT VISCOSITY (TIME STEP) +C VISCR | NEL | F | W | ELEMENT ROTATION VISCOSITY (TIME STEP) +C----------+---------+---+---+-------------------------------------------- +C UVAR |NUVAR*NEL| F | W | USER ELEMENT VARIABLES +C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES +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,KPROP) +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) NUMBERS. +C SEE LECG29 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 NUMBERS(function +C refered by users materials). +C SEE LECM29 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------------------------------------------------------------------------- +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 MID = 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 PID = 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----------------------------------------------- +C I m p l i c i t T y p e s +C----------------------------------------------- + IMPLICIT NONE +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,NEL,IPROP, + . IX(4,NEL) , + . GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU, + . KFUNC,KMAT,KPROP + double precision + . XL(NEL) ,MASS(NEL) ,XINER(NEL) ,STIFM(NEL) , + . STIFR(NEL),VISCM(NEL) ,VISCR(NEL),UVAR(NUVAR,*), + . GET_U_MAT,GET_U_GEO + EXTERNAL GET_U_PNU,GET_U_MNU,GET_U_MAT,GET_U_GEO,GET_U_PID, + . GET_U_MID + 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----------------------------------------------- + double precision + . AMAS,ELASTIF,DIAMETER,ActOpt_f,q0,STIM_ID + INTEGER I + DOUBLE PRECISION UN,ZERO +C----------------------------------------------- + UN=1.0 + ZERO = 0.0 + print *,' *** 1: rini START' + + AMAS = GET_U_GEO(1,IPROP) + ELASTIF = GET_U_GEO(2,IPROP) + DIAMETER = GET_U_GEO(3,IPROP) + ActOpt_f = GET_U_GEO(4,IPROP) + STIM_ID = GET_U_GEO(6,IPROP) + q0 = GET_U_GEO(7,IPROP) +C-------------------------------------- +C ELEMENT CHECK +C-------------------------------------- + DO I=1,NEL + IF(XL(I).EQ.0.0)THEN + WRITE(IOUT,*)' **ERROR: ZERO LENGTH SPRING :' + ENDIF + ENDDO +C-------------------------------------- +C ELEMENT INITIALIZATION +C-------------------------------------- + DO I=1,NEL + MASS(I) = AMAS*3.14*((DIAMETER/2)**2)*XL(I) !here the real mass has to be computed density * volume + XINER(I) = ZERO +C +C +C FOR NODAL AND ELEMENT TIME STEP COMPUTATION + STIFM(I) = ELASTIF + STIFR(I) = ZERO + VISCM(I) = ZERO + VISCR(I) = ZERO + print *,' *** 1: rini values: I =',I + print *,' *** 1: rini values: AMAS =',AMAS + print *,' *** 1: rini values: DIAMETER =',DIAMETER + print *,' *** 1: rini values: XL(I) =',XL(I) + print *,' *** 1: rini values: MASS =',MASS(I) + print *,' *** 1: rini values: ActOpt_f =',ActOpt_f + print *,' *** 1: rini values: q0 =',q0 + print *,' *** 1: rini values: STIM_ID =',STIM_ID + print *,' *** 1: rini values: STIFM(I) =',STIFM(I) + + WRITE(IOUT,1000) + . AMAS,MASS(I) + + ENDDO +C + print *,' *** rini END' + + RETURN + 1000 FORMAT( + & 5X,' Extended Hill-type Muscle Model INITAILIZATION :',/, + & 5X,'AMAS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/, + & 5X,'ELEMENT MASS. . . . . . . . . . . . . . . . . . . . . . . . . . . . .=',E12.4/) + + END diff --git a/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Code/ruser30.F b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Code/ruser30.F new file mode 100644 index 0000000..d89fc58 --- /dev/null +++ b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Code/ruser30.F @@ -0,0 +1,1405 @@ +C MIT License +C +C Copyright (c) 2025 RM412 +C +C Permission is hereby granted, free of charge, to any person obtaining a copy +C of this software and associated documentation files (the "Software"), to deal +C in the Software without restriction, including without limitation the rights +C to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +C copies of the Software, and to permit persons to whom the Software is +C furnished to do so, subject to the following conditions: +C +C The above copyright notice and this permission notice shall be included in all +C copies or substantial portions of the Software. +C +C THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +C IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +C FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +C AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +C LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +C OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +C SOFTWARE. +C +C Suggestion : +C +C BSD 3-Clause License +C Copyright (c) 2014 +C D. Haeufle, M. Günther, A. Bayer, S. Schmitt +C Copyright (c) 2017 +C C. Kleinbach, O. Martynenko, J. Promies, D. F. B. Haeufle, J. Fehr, S. Schmitt +C Copyright (c) 2019 +C O. Martynenko, F. Kempter, C. Kleinbach, S. Schmitt, J. Fehr +C Copyright (c) 2020 +C P. Lerge, O. Martynenko, L. Nölle, S. Schmitt, I. Wochner +C Copyright (c) 2022 +C O. Martynenko, F. Kempter, C. Kleinbach, L. Nölle, P. Lerge, S. Schmitt, J. Fehr +C Copyright (c) 2026 +C A. Ardisson, E. Wagnac, M.-H. Beauséjour/ÉTS +C Adaptation and port to OpenRadioss +C All rights reserved. +C Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following C C conditions are met: +C Redistributions of source code must retain the above copyright notice, +C this list of conditions and the following disclaimer. +C +C +C Redistributions in binary form must reproduce the above copyright notice, +C this list of conditions and the following disclaimer in the documentation +C and/or other materials provided with the distribution. +C +C +C Neither the name of the copyright holder nor the names of its contributors +C may be used to endorse or promote products derived from this software +C without specific prior written permission. +C +C +C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” +C AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +C IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +C ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +C LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +C POSSIBILITY OF SUCH DAMAGE. + +! C----------------------------------------------- +! C calc_F_isom FUNCTION !output calc_F_ISOM +! C---------------------------------------------- + double precision FUNCTION calc_F_isom(l_CE,l_CE_opt,cm) +! ! C! Computation of the isometric force developped by the muscle + double precision :: cm(*) + double precision :: l_CE , f_isom_temp,l_CE_opt + IF (l_CE.GE.l_CE_opt) THEN + f_isom_temp = exp(-(abs(((l_CE/l_CE_opt)-1.0)/cm(11)))**cm(12)) + ELSE + f_isom_temp = exp(-(abs(((l_CE/l_CE_opt)-1.0)/cm(13)))**cm(14)) + END IF +C Additional code to detect values less or equal than zero + IF (f_isom_temp.LE.1E-13) THEN + calc_F_isom = 1E-13 + ELSE + calc_F_isom = f_isom_temp + END IF + RETURN + END FUNCTION calc_F_isom +C +C----------------------------------------------- +C calc_F_PEE FUNCTION ! output calc_F_PEE +C---------------------------------------------- +C! Force of the parallel elastic element PEE + double precision FUNCTION calc_F_PEE(l_CE, cm, l_PEE0, K_PEE) + double precision :: cm(*) + double precision :: l_CE, l_PEE0,K_PEE + IF (l_CE.GE.l_PEE0) THEN + calc_F_PEE = K_PEE*((l_CE-l_PEE0)**cm(20)) + ELSE + calc_F_PEE = 0.0d0 + END IF + RETURN + END FUNCTION calc_F_PEE +C +C----------------------------------------------- +C calc_F_SEE FUNCTION !output calc_F_SEE +C---------------------------------------------- +C! Force of the serial elastic element SEE + double precision function calc_F_SEE(l_CE,l_SEE_0,cm,elleng,l_SEE_nll,v_SEE,K_SEE_nl,K_SEE_l) + double precision :: cm(*) + double precision :: l_CE,K_SEE_nl,K_SEE_l,l_SEE,l_SEE_nll + double precision :: v_SEE,elleng,l_SEE_0 + l_SEE = elleng - l_CE + IF ((l_SEE.LT.l_SEE_nll).and.(l_SEE.GT.l_SEE_0)) THEN + calc_F_SEE = K_SEE_nl*((l_SEE-l_SEE_0)**v_SEE) + ELSE IF ((l_SEE.GE.l_SEE_nll).and.(l_SEE.GT.l_SEE_0)) THEN + calc_F_SEE = cm(25) + K_SEE_l*(l_SEE - l_SEE_nll) + ELSE + calc_F_SEE = 0.0d0 + ENDIF + RETURN + END FUNCTION calc_F_SEE +C +C----------------------------------------------- +C calc_F_SUM FUNCTION +C----------------------------------------------- +C + DOUBLE PRECISION FUNCTION calc_F_SUM(l_CE,l_CE_opt,l_SEE_0,elleng,q,cm, l_PEE0,K_PEE,l_SEE_nll,v_SEE,K_SEE_nl, K_SEE_l) + DOUBLE PRECISION, dimension(39) :: cm + DOUBLE PRECISION l_PEE0, K_PEE, v_SEE, K_SEE_nl, K_SEE_l,l_SEE_nll,l_CE + DOUBLE PRECISION F_isom, F_PEE, F_SEE, F_CE_init, elleng, q + DOUBLE PRECISION calc_F_isom, calc_F_PEE, calc_F_SEE,l_SEE_0,l_CE_opt + F_isom = calc_F_isom(l_CE,l_CE_opt,cm) + F_PEE = calc_F_PEE(l_CE, cm, l_PEE0, K_PEE) + F_SEE = calc_F_SEE(l_CE,l_SEE_0,cm,elleng,l_SEE_nll,v_SEE,K_SEE_nl,K_SEE_l) + F_CE_init = cm(9)*q*F_isom + calc_F_SUM = F_SEE - F_CE_init - F_PEE + RETURN + END FUNCTION calc_F_SUM +C +C----------------------------------------------- +C calc_STIM FUNCTION +C----------------------------------------------- +C + DOUBLE PRECISION FUNCTION calc_STIM(cm,UVAR,tt,lCEdelay,dotlCEdelay,ncycle,NUVAR,I,IFUNC1,IFUNC2,IFUNC3,GET_U_FUNC) + INTEGER NUVAR,I,ncycle,IFUNC1,IFUNC2,IFUNC3 + DOUBLE PRECISION cm(*),UVAR(NUVAR,I) + DOUBLE PRECISION lCEdelay, dotlCEdelay + DOUBLE PRECISION STIM, lambda, dlambda + DOUBLE PRECISION STIM_lambda, STIM_hybrid, STIM_reflex + DOUBLE PRECISION STIM_closed,STIM_open,tt,DXDY + DOUBLE PRECISION GET_U_FUNC +c +c get STIM + SELECT CASE (int(cm(33))) + CASE (0) +c NO CONTROLLER + IF (cm(2).GT.0.0) THEN +C STIM level from curve + STIM = GET_U_FUNC(IFUNC1,tt,DXDY) + ELSE +c constant STIM level + STIM = abs(cm(2)) + END IF + IF (STIM.GT.1.0) then + STIM = 1.0 + ELSE IF (STIM.LT.0.0) THEN + STIM = 0.0 + ENDIF + CASE (1) +c LAMBDA CONTROLLER + IF (cm(34).LT.0.0) THEN +c lambda from curve + lambda = GET_U_FUNC(IFUNC3,tt,DXDY) + dlambda=0.0 + ELSE +c lambda = constant + lambda = abs(cm(34)) + dlambda=0.0 + ENDIF + IF (lambda.LT.0.0) THEN + lambda = 0.0 + ENDIF + UVAR(50,I)=lambda + IF (tt.LT.cm(38)) THEN + IF (cm(2).GT.0.0) THEN +c STIM from curve + STIM = GET_U_FUNC(IFUNC1,tt,DXDY) + ELSE +c constant STIM level + STIM = abs(cm(2)) + ENDIF + IF (STIM.GT.1.0) THEN + STIM = 1.0 + ELSE IF (STIM.LT.0.0) THEN + STIM = 0.0 + ENDIF + ELSE + STIM = STIM_lambda(cm ,l_CE_opt,UVAR,lCEdelay,dotlCEdelay,tt,lambda,dlambda,NUVAR,I) + ENDIF + CASE (2,3) +c HYBRID CONTROLLER STANDARD and KISTEMAKER + IF (cm(2).GT.0.0) THEN +c STIM from curve + STIM = GET_U_FUNC(IFUNC1,tt,DXDY) + ELSE +c constant STIM level + STIM = abs(cm(2)) + ENDIF + IF (STIM.GT.1.0) THEN + STIM = 1.0 + ELSE IF (STIM.LT.0.0) THEN + STIM = 0.0 + ENDIF + IF (cm(34).LT.0.0) THEN +c lambda from curve + lambda = GET_U_FUNC(IFUNC3,tt,DXDY) + dlambda=0.0 + ELSE +c lambda = constant + lambda = abs(cm(34)) + dlambda=0.0 + ENDIF + IF (lambda.LT.0.0) THEN + lambda = 0.0 + ENDIF + UVAR(50,I) = lambda + IF (tt.GE.cm(38)) THEN + STIM_closed = STIM_lambda(cm ,l_CE_opt,UVAR,lCEdelay,dotlCEdelay,tt,lambda,dlambda,NUVAR,I) + STIM = STIM_hybrid(cm,l_CE_opt,UVAR,STIM_open,lambda,dlambda,lCEdelay,DOtlCEdelay,tt,NUVAR,I) + ENDIF + CASE (4) +c REFLEX CONTROLLER + IF (tt.LT.cm(38)) THEN + IF (cm(2).GT.0.0) THEN +c STIM from curve + STIM = GET_U_FUNC(IFUNC1,tt,DXDY) + ELSE +c constant STIM level + STIM = abs(cm(2)) + ENDIF + IF (STIM.GT.1.0) THEN + STIM = 1.0 + ELSE IF (STIM.LT.0.0) THEN + STIM = 0.0 + END IF + ELSE + STIM = STIM_reflex(ncycle,cm,UVAR,tt,I,NUVAR) + ENDIF + CASE DEFAULT + STIM = 0 + END SELECT + calc_STIM = STIM + UVAR(2,I) = STIM + RETURN + END FUNCTION calc_STIM +C +C----------------------------------------------- +C calc_dgam FUNCTION +C----------------------------------------------- +C + DOUBLE PRECISION FUNCTION calc_dgam(cm,STIM,UVAR,NUVAR,I) + DOUBLE PRECISION cm(*), UVAR(NUVAR,I) + DOUBLE PRECISION :: STIM + INTEGER I +C + calc_dgam = cm(7)*(STIM-UVAR(15,I)) +C + RETURN + END FUNCTION calc_dgam +C +C----------------------------------------------- +C calc_dq FUNCTION +C----------------------------------------------- +c + DOUBLE PRECISION FUNCTION calc_dq(UVAR,cm,STIM,tau_q,NUVAR,I) +c calculates activation with the differential equation +c + INTEGER NUVAR,I + DOUBLE PRECISION cm(*),UVAR(NUVAR,I) + DOUBLE PRECISION STIM,tau_q +c + calc_dq = 1/tau_q*(STIM-STIM*(1.0-cm(5))*(UVAR(3,I)-cm(3))-cm(5) + 1 *(UVAR(3,I)-cm(3))) + RETURN + END FUNCTION calc_dq +C----------------------------------------------- +C calc_q FUNCTION +C----------------------------------------------- + + DOUBLE PRECISION FUNCTION calc_q(cm,l_CE_opt,UVAR,tt,STIM,DT,NUVAR,I,GET_U_FUNC,IFUNC2) + INTEGER NUVAR,I,IFUNC2 + DOUBLE PRECISION cm(*), UVAR(NUVAR,I) + DOUBLE PRECISION STIM, q, gam_rel, rho_act, tau + DOUBLE PRECISION calc_dq, calc_dgam,tt, DT,GET_U_FUNC,l_CE_opt + q = cm(3) + IF (cm(1).EQ.1.0) THEN +c ZAJAC activation dynamics +c get tau + IF (cm(4).LT.0.0) THEN +c tau from curve + tau = GET_U_FUNC(IFUNC2,tt,DXDY) + ELSE + tau = cm(4) + ENDIF +c calc new q + q = UVAR(3,I) + calc_dq(UVAR,cm,STIM,tau,NUVAR,I)*DT + ELSE IF (cm(1).EQ.2.0) THEN +c HATZE activation dynamics +c calc new gamma + gam_rel = UVAR(15,I) + calc_dgam(cm,STIM,UVAR,NUVAR,I)*DT + UVAR(15,I) = gam_rel +c calc new rho + rho_act = cm(4)*cm(5)*(cm(6)-1.0)/(cm(6)-(UVAR(10,I)/l_CE_opt)) + 1 *(UVAR(10,I)/l_CE_opt) +c calc new q + q = (cm(3)+(rho_act*gam_rel)**3)/(1.0+(rho_act*gam_rel)**3) + ELSE IF (cm(1).EQ.3.0) THEN +c Modified HATZE activation dynamics by Rockenfeller, Günther 2018 +c calc new gamma + gam_rel = UVAR(15,I) + calc_dgam(cm,STIM,UVAR,NUVAR,I)*DT + UVAR(15,I) = gam_rel +c calc new rho + rho_act = cm(4)*cm(5)*UVAR(10,I)/l_CE_opt +c calc new q + q = (cm(3)+(rho_act*gam_rel)**3)/(1.0+(rho_act*gam_rel)**3) + ENDIF + IF (q.GT.1.0) THEN + q = 1.0 + ELSE IF (q.LT.cm(3)) THEN + q = cm(3) + ENDIF + UVAR(3,I) = q + calc_q = q + RETURN + END FUNCTION calc_q +C----------------------------------------------- +C delaylCE FUNCTION +C----------------------------------------------- + DOUBLE PRECISION function delaylCE(ncycle,UVAR,cm,tt,NUVAR,I) + DOUBLE PRECISION lasttime, DT1,tt + INTEGER NUVAR,I,ncycle + DOUBLE PRECISION cm(*),UVAR(NUVAR,I) + INTEGER index1,index2,indexr1, indexr2, iter +c param declaration + lasttime = UVAR(34,I); +c dt = delay/buffersize + DT1 = cm(37)/UVAR(30,I); + indexr1 = UVAR(32,I); + indexr2 = UVAR(33,I); + IF (ncycle.LE.3) THEN +c initialisation +c set ring buffer values to 0. + DO iter =1 , (UVAR(30,I)) + UVAR((UVAR(31,I)+ iter - 1),I) = 0; + ENDDO + ELSEIF (tt.gt.(lasttime + DT1)) THEN + lasttime = lasttime + DT1 +c mod ensures jump from end to the beginning of the ringbuffer + indexr1 = mod((UVAR(32,I)+1),(UVAR(30,I))) + indexr2 = mod((UVAR(33,I)+1),(UVAR(30,I))) +c save current lce + index1 = UVAR(31,I) + indexr1 + UVAR(index1,I) = UVAR(10,I); + ENDIF +c define absolute index position of current delayedlCE value + index2 = UVAR(31,I) + indexr2 +c return delayedlCE + delaylCE = UVAR(index2,I) + UVAR(32,I) = indexr1 + UVAR(33,I) = indexr2 + UVAR(34,I) = lasttime + return + end +c +C----------------------------------------------- +C delaydotlCE FUNCTION +C----------------------------------------------- + DOUBLE PRECISION FUNCTION delaydotlCE(ncycle,UVAR,cm,tt,NUVAR,I) + DOUBLE PRECISION lasttimed,DT1,tt + integer ncycle,NUVAR,I + DOUBLE PRECISION cm(*),UVAR(NUVAR,I) + integer indexd1,indexd2,indexdr1,indexdr2,iterd +c ---------initialisation ------ +c param declaration + lasttimed = UVAR(39,I) + DT1 = cm(37)/UVAR(30,I) + indexdr1 = UVAR(37,I) + indexdr2 = UVAR(38,I) + IF (ncycle.LE.3) THEN +c set ring buffer values to 0. + DO iterd=1,(UVAR(30,I)) + UVAR((UVAR(36,I) + iterd-1),I) = 0 + ENDDO + ELSEIF (tt.gt.(lasttimed + DT1)) THEN + lasttimed = lasttimed + DT1 +c mod ensures jump from end to the beginning of the ringbuffer + indexdr1 = mod((UVAR(37,I)+1),(UVAR(30,I))) + indexdr2 = mod((UVAR(38,I)+1),(UVAR(30,I))) + indexd1 = UVAR(36,I) + indexdr1 +c save current dot_lce + UVAR(indexd1,I) = UVAR(12,I) + ENDIF +c define absolute index position of current delayeddotlCE value + indexd2 = UVAR(36,I) + indexdr2 +c return delayeddotlCE + delaydotlCE = UVAR(indexd2,I) + UVAR(39,I) = lasttimed + UVAR(37,I) = indexdr1 + UVAR(38,I) = indexdr2 + return + end +C----------------------------------------------- +C ZEROIN FUNCTION +C----------------------------------------------- +C + DOUBLE PRECISION function ZEROIN(elleng,l_CE_opt,l_SEE_0,act,cm,l_SEE_nll,K_SEE_nl,v_SEE,l_PEE0,K_PEE,K_SEE_l) + IMPLICIT NONE +C!----------------------------------------------- +c! Root-finding algorithm using the ZEROIN function +c! from Sec. 7.2 of: +c! +c! Forsythe, G.E.; Malcolm, M.A.; Moler, C.B.: Computer Methods for +c! Mathematical Computations. Prentice Hall Professional Technical +c! Reference, 1977 +c! +c! as interval we used [0,elleng], the function to evaluate F(X) +c! is calc_F_sum(X,elleng,act,cm). + DOUBLE PRECISION elleng, act,l_SEE_nll,K_SEE_nl + DOUBLE PRECISION cm(*) +c! a zero of the function F(X) is computed in the Interval AX,BX +c! INPUT: +c! AX LEFT ENDPOINT OF INITIAL INTERVAL +c! BX RIGHT ENDPOINT OF INITIAL INTERVAL +c! F FUNCTION SUBPROGRAM WHICH EVALUATES F(X) +c! ANY X IN THE INTERVAL AX,BX +c! TOL DESIRED LENGTH OF THE INTERVAL OF UNCERTAINTY +c! OF THE FINAL RESULT (.GE. 0.0) +c! +c! OUTPUT: +c! ZEROIN ABCISSA APPROXIMATING A ZERO OF F IN THE INTERVAL AX,BX + DOUBLE PRECISION A,B,C,D,E,EPS,FA,FB,FC,TOL1,XM,P,Q,R,S,TOL + DOUBLE PRECISION AX,BX, TOL_CALC + DOUBLE PRECISION calc_F_SUM,v_SEE,l_PEE0,K_PEE,K_SEE_l,l_CE_opt,l_SEE_0 +c! +c! Get the numerical machine tolerance +c! + TOL = 2.2E-16 + EPS = 1.0 + 10 EPS = EPS/2.0 + TOL_CALC = 1.0 + EPS + IF (TOL_CALC.GT.1.0) GO TO 10 +c! FIND INTERVAL + AX = 0 + BX = elleng +c INITIALIZATION + A = AX + B = BX + FA = calc_F_SUM(A,l_CE_opt,l_SEE_0,elleng,act,cm, l_PEE0,K_PEE, + 1 l_SEE_nll,v_SEE,K_SEE_nl, K_SEE_l) +C + FB = calc_F_SUM(B ,l_CE_opt,l_SEE_0,elleng,act,cm, l_PEE0, + 2 K_PEE,l_SEE_nll,v_SEE,K_SEE_nl, K_SEE_l) +c BEGIN STEP +20 C = A + FC = FA + D = B-A + E = D +30 IF (ABS(FC).GE.ABS(FB)) GO TO 40 + A = B + B = C + C = A + FA = FB + FB = FC + FC = FA +c CONVERGENCE TEST +40 TOL1= 2.0*EPS*ABS(B)+0.5*TOL + XM = 0.5*(C-B) + IF (ABS(XM).LE.TOL1) GO TO 90 + IF (FB.EQ.0.0) GO TO 90 +c IS BISECTION NECESSARY + IF (ABS(E).LT.TOL1) GO TO 70 + IF (ABS(FA).LE.ABS(FB)) GO TO 70 +c IS QUADRATIC INTERPOLATION POSSIBLE + IF (A.NE.C) GO TO 50 +c LINEAR INTERPOLATION + S = FB/FA + P = 2.0d0*XM*S + Q = 1.0d0-S + GO TO 60 +c INVERSE QUADRATIC INTERPOLATION +50 Q = FA/FC + R = FB/FC + S = FB/FA + P = S*(2.0*XM*Q*(Q-R)-(B-A)*(R-1.0)) + Q = (Q-1.0)*(R-1.0)*(S-1.0) +c ADJUST SIGNS +60 IF (P.GT.0.0) Q=-Q + P = ABS(P) +c IS INTERPOLATION ACCEPTABLE + IF ((2.0*P) .GE. (3.0*XM*Q-ABS(TOL1*Q))) GO TO 70 + IF (P .GE. ABS(0.5*E*Q)) GO TO 70 + E = D + D = P/Q + GO TO 80 +c BISECTION +70 D = XM + E = D +c COMPLETE STEP +80 A = B + FA = FB + IF (ABS(D) .GT. TOL1) B = B+D + IF (ABS(D) .LE. TOL1) B = B+SIGN(TOL1,XM) + FB = calc_F_SUM(B,l_CE_opt,l_SEE_0,elleng,act,cm, l_PEE0, + 1 K_PEE,l_SEE_nll,v_SEE,K_SEE_nl, K_SEE_l) + IF ((FB*(FC/ABS(FC))) .GT. 0.0) GO TO 20 + GO TO 30 +c DONE +90 ZEROIN=B + RETURN + END FUNCTION ZEROIN +C +C----------------------------------------------- +C STIM_lambda FUNCTION ! output STIM_lambda +C---------------------------------------------- +C + DOUBLE PRECISION FUNCTION STIM_lambda(cm ,l_CE_opt,UVAR,lCEdelay,dotlCEdelay,tt,lambda,dlambda,NUVAR,I) + DOUBLE PRECISION cm(*),UVAR(NUVAR,I) + DOUBLE PRECISION lCEdelay,dotlCEdelay,lambda, tt, dlambda,l_CE_opt + INTEGER I,NUVAR +C + STIM_lambda = (cm(35)*(lCEdelay-lambda)+cm(36)*(dotlCEdelay-dlambda))/l_CE_opt + IF (cm(33).EQ.3) THEN + IF (STIM_lambda.GT.1.0) THEN + STIM_lambda = 1.0d0 + ELSE IF (STIM_lambda.LT.-1.0) THEN + STIM_lambda = -1.0d0 + ELSE + CONTINUE + END IF + ELSE + IF (STIM_lambda.GT.1.0) THEN + STIM_lambda = 1.0d0 + ELSE IF (STIM_lambda.LT.0.0) THEN + STIM_lambda = 0.0d0 + ELSE + CONTINUE + END IF + END IF + RETURN + END FUNCTION STIM_lambda +C +C----------------------------------------------- +C Hybrid controller FUNCTION ! output STIM_hybrid +C---------------------------------------------- +C + DOUBLE PRECISION function STIM_hybrid(cm,l_CE_opt,UVAR,STIM_open,lambda,dlambda,lCEdelay,DOtlCEdelay,tt,NUVAR,I) + IMPLICIT NONE + INTEGER NUVAR, I + DOUBLE PRECISION cm(*),UVAR(NUVAR,I) + DOUBLE PRECISION STIM_open, lambda, dlambda, tt + DOUBLE PRECISION lCEdelay, DOtlCEdelay,l_CE_opt + IF (lCEdelay.GT.lambda) THEN + STIM_hybrid=(STIM_open+(cm(35)*(lCEdelay-lambda)+cm(36)*(DOtlCEdelay-dlambda))/l_CE_opt) + ELSE + STIM_hybrid=STIM_open + END IF + IF (STIM_hybrid.GT.1.0) THEN + STIM_hybrid = 1.0 + ELSE IF (STIM_hybrid.LT.0.0) THEN + STIM_hybrid = 0.0 + END IF + RETURN + END FUNCTION STIM_hybrid +C +C----------------------------------------------- +C Reflexive Controller function ! output STIM_reflex +C---------------------------------------------- +C + DOUBLE PRECISION function STIM_reflex(ncycle,cm,UVAR,tt,I,NUVAR) + IMPLICIT NONE + DOUBLE PRECISION cm(*),UVAR(NUVAR,*) + DOUBLE PRECISION tt + INTEGER ncycle,I,NUVAR +c! initialization of strain (UVAR(23)) + UVAR(23,I) = 0.0d0 +c! initialization of l_CE_ref = l_CE(t=0) OR +c! IF time < t_PreSim THEN l_CE_ref = l_CE -- l_CE_ref = l_CE(t=t_PreSim) + IF (ncycle.LE. 3 .OR. tt .LE. cm(38).OR. UVAR(22,I).EQ.0.) THEN + UVAR(22,I) = UVAR(10,I) + END IF +c! UVAR(22) (=l_CE_ref) = l_CE(t=t_PreSim) +c! IF time t_PreSim + delay: calc strain (strain = 0 for t UVAR +C----------------------------------------------- + . STIM,q,F_MTC, F_CE,F_PEE,F_SEE,F_SDE,l_MTC,l_CE, DOt_l_MTC, + . DOt_l_CE, counter_output, F_isom, gam_rel, l_MTC_0, + . lCEdelay, l_CE_ref, strain, STIM_reflex_prev, + . buffersize, idx_begin_lCEbuffer, lasttime, calc_F_isom, calc_F_PEE, calc_F_SEE, + . calc_F_SUM, calc_dq,calc_q, calc_dgam, epsilon,tau, lambda,dlambda, + . rho_act, STIM_lambda, STIM_hybrid, STIM_reflex,dq,l_CE_opt,l_SEE_0 +C----------------------------------------------- +C----------------------------------------------- + DOUBLE PRECISION + . cm(39) +C ----------------------------------------------- + CHARACTER*100 FILENAME, COMMAND ! USED FOR VARIABLES CHECK + INTEGER J, IO_STATUS + DOUBLE PRECISION :: PLOT_X, PLOT_Y, test1, test2, test3 +C----------------------------------------------- +C----------------------------------------------- + UN = 1.0d0 + ZERO = 0.0d0 +c! ------------------------------------------------------------ +c! ------ history variables - overview ------------------------ +c! ------------------------------------------------------------ +c! UVAR(1) = sig(1) +c! uvar(2) = STIM +c! uvar(3) = q +c! uvar(4) = F_MTC +c! uvar(5) = F_CE +c! uvar(6) = F_PEE +c! uvar(7) = F_SEE +c! uvar(8) = F_SDE +c! uvar(9) = elleng (l_MTC) +c! uvar(10) = l_CE +c! uvar(11) = dot_l_MTC +c! uvar(12) = dot_l_CE +c! uvar(13) = counter_output +c! uvar(14) = F_isom +c! uvar(15) = gam_rel +c! uvar(16) = l_MTC_0 (initial element length) +c! ------------------------------------------------------------ +c! ------ only needed if controller are used ------------------ +c! ------------------------------------------------------------ +c! uvar(20) = lCEdelay +c! uvar(21) = dotlCEdelay +c! uvar(22) = l_CE_ref (for reflex controller) +c! uvar(23) = strain (for reflex controller) +c! uvar(24) = STIM_reflex_prev (STIM value of reflex controller in previous timestep) +c! ------------------------------------------------------------ +c! ------ delay buffer----------------------------------------- +c! ------------------------------------------------------------ +c! uvar(30) = buffersize +c! uvar(31) = idx_begin_lCEbuffer (uvar(142:145) seems to be used internally by lsdyna, suggestion: use uvar(31)=150) +c! uvar(32) = indexr1...index1 of lce-ringbuffer (not eq index of uvar) +c! uvar(33) = indexr1...index2 of lce-ringbuffer (not eq index of uvar) +c! uvar(34) = lasttime +c! uvar(36) = begindotlCE +c! uvar(37) = indexdotr1...index1 of dotlce-ringbuffer (not eq index of uvar) +c! uvar(38) = indexdotr2...index2 of dotlce-ringbuffer (not eq index of uvar) +c! uvar(39) = dotlasttime +c! uvar(uvar(31):uvar(31)+uvar(30)-1) = ringbuffer_l_CE +c! uvar(uvar(36):uvar(36)+uvar(30)-1) = ringbuffer_dot_l_CE +C---------------------------------------------------------- +! FILENAME = "C:/Work/Tools/userlib_sdk/examples/USER_dump/"// ! PATH FOR CHEKING PLOT FILES +! 1 "Springs/EXCEL_TESTS/" // +! 2 "EXCEL_TEST_DATA.csv" + FILENAME = "C:/Memoire/Etude_sensibilite/MISE_A_JOUR/OBSERVATIONS"// + 1 "/OBSERVATIONS.csv" +C---------------------------------------------------------- + !print *,' ****** ruser START' + ncycle = GET_U_CYCLE() + tt = GET_U_TIME() + + ELASTIF = GET_U_GEO(2,IPROP) + DIAMETER = GET_U_GEO(3,IPROP) + + AMAS = GET_U_GEO(1,IPROP) +c! IERROR = SET_U_GEO(2,ELASTIF) +c! Variables (activation dynamic) + ActOpt = INT(GET_U_GEO(4,IPROP)+0.0001) + ActOpt_f = GET_U_GEO(4,IPROP) + cm(1) = INT(ActOpt_f) + SENS_ID = INT(GET_U_GEO(5,IPROP)+0.0001) + + STIM_ID = INT(GET_U_GEO(6,IPROP)+0.0001) + STIM_ID_f = GET_U_GEO(6,IPROP) + cm(2) = STIM_ID_f ! *** is/can be evtl. a function // NOT USED YET + +c! IFUNC1 = GET_U_PNU(1,IPROP,KFUNC) +c! +c! IF (IFUNC1.EQ.ZERO)THEN +c! SCALE_D1 = 1.0 +c! ELSE +c! SCALE_D1 = GET_U_FUNC(IFUNC1,TT,DXDY) +c! ENDIF +C---------------------------------------------------------- + +C---------------------------------------------------------- + q0 = GET_U_GEO(7,IPROP) + cm(3) = q0 + tauq = GET_U_GEO(8,IPROP) + cm(4) = tauq + betaq = GET_U_GEO(9,IPROP) + cm(5) = betaq + k = GET_U_GEO(10,IPROP) + cm(6) = k + + m = GET_U_GEO(11,IPROP) + cm(7) = m + muscle = GET_U_GEO(12,IPROP) + cm(8) = muscle +c! +c! Variables (isometric force) +c! + Fmax = GET_U_GEO(13,IPROP) + cm(9) = Fmax + t_lCEopt = GET_U_GEO(14,IPROP) + cm(10) = t_lCEopt + dWdes = GET_U_GEO(15,IPROP) + cm(11) = dWdes + + nuCEdes = GET_U_GEO(16,IPROP) + cm(12) = nuCEdes + dWasc = GET_U_GEO(17,IPROP) + cm(13) = dWasc + nuCEasc = GET_U_GEO(18,IPROP) + cm(14) = nuCEasc +c! +c! Variables (Hill Parameter, concentric) +c! + + Arel0 = GET_U_GEO(19,IPROP) + cm(15) = Arel0 + Brel0 = GET_U_GEO(20,IPROP) + cm(16) = Brel0 +c! +c! Variables (van Soest Parameter, eccentric) +c! + Secc = GET_U_GEO(21,IPROP) + cm(17) = Secc + Fecc = GET_U_GEO(22,IPROP) + cm(18) = Fecc +c! +c! Variables (Parallel elastic element) +c! + + LPEE0 = GET_U_GEO(23,IPROP) + cm(19) = LPEE0 + nuPEE = GET_U_GEO(24,IPROP) + cm(20) = nuPEE + FPEE = GET_U_GEO(25,IPROP) + cm(21) = FPEE +c! +c! Variables (Seriell elastic element) +c! + + dUSEEnll = GET_U_GEO(27,IPROP) + cm(23) = dUSEEnll + duSEEl = GET_U_GEO(28,IPROP) + cm(24) = duSEEl + dFSEE0 = GET_U_GEO(29,IPROP) + cm(25) = dFSEE0 +c! +c! Variables (Damping element) +c! + + Damping = GET_U_GEO(30,IPROP) + cm(26) = Damping + Param1 = GET_U_GEO(31,IPROP) + cm(27) = Param1 + Param2 = GET_U_GEO(32,IPROP) + cm(28) = Param2 +c! +c! Variables (Output definition; musout.(partid)) +c! +c! cm(29)=output method (EQ.0. no output +c! EQ.1. advanced output (basic output plus dot_l_CE, dot_l_MTC, lCEdelay, dotlCEdelay)) +c! + output_MTD = INT(GET_U_GEO(33,IPROP)+0.0001) + cm(29) = GET_U_GEO(33,IPROP) + + timeste = GET_U_GEO(34,IPROP) + cm(30) = timeste +c! +c! Variables for Controller +c! +c! cm(33) = Activation Method (EQ.1. lambda_controller +c! EQ.2. hybrid_controller +c! EQ.3. reflexive controller) +c! + Activation = GET_U_GEO(35,IPROP) + cm(33) = Activation +C + target_l_CE = GET_U_GEO(36,IPROP) + cm(34) = target_l_CE +C + kp = GET_U_GEO(37,IPROP) + cm(35) = kp + kd = GET_U_GEO(38,IPROP) + cm(36) = kd + + dotlCEdelay = GET_U_GEO(39,IPROP) + cm(37) = dotlCEdelay + t_PreSim = GET_U_GEO(40,IPROP) + cm(38) = t_PreSim + threshold = GET_U_GEO(41,IPROP) + cm(39) = threshold + + IFUNC1 = GET_U_PNU(1,IPROP,KFUNC) + IFUNC2 = GET_U_PNU(2,IPROP,KFUNC) + IFUNC3 = GET_U_PNU(3,IPROP,KFUNC) +C +C initialization of initial laength + IF (tt.EQ.0.0) THEN + DO I = 1,NEL + !UVAR(1,I) = XL(I) + UVAR(16,I) = XL(I) + cm(8) + print *,'' + print *,'*************************************' + print *,'* VERSION from 27.06.2024 - ver. a1 *' + print *,'*************************************' + print *,'' + + print *,'Spring #= ',GET_SPRING_ELNUM(I),' ; MASS =',MASS(I),' Initialization:' + print *,'AMAS =',AMAS + print *,'ELASTIF ',ELASTIF + print *,'DIAMETER ',DIAMETER + print *,'ActOpt = ',ActOpt + print *,'ActOpt_f = ',ActOpt_f , cm(1) + print *,'SENS_ID = ',SENS_ID + + print *,'STIM_ID = ',STIM_ID , cm(2) + print *,'q0 = ',q0 , cm(3) + print *,'tauq = ',tauq , cm(4) + print *,'betaq = ',betaq , cm(5) + print *,'k = ',k , cm(6) + + print *,'m = ',m , cm(7) + print *,'muscle = ',muscle , cm(8) + print *,'Fmax = ',Fmax , cm(9) + print *,'t_lCEopt = ',t_lCEopt , cm(10) + print *,'dWdes = ',dWdes , cm(11) + + print *,'nuCEdes = ',nuCEdes , cm(12) + print *,'dWasc = ',dWasc , cm(13) + print *,'nuCEasc = ',nuCEasc , cm(14) + print *,'Arel0 = ',Arel0 , cm(15) + print *,'Brel0 = ',Brel0 , cm(16) + + print *,'Secc = ',Secc , cm(17) + print *,'Fecc = ',Fecc , cm(18) + print *,'LPEE0 = ',LPEE0 , cm(19) + print *,'nuPEE = ',nuPEE , cm(20) + print *,'FPEE = ',FPEE , cm(21) + + print *,'dUSEEnll = ',dUSEEnll , cm(23) + print *,'duSEEl = ',duSEEl , cm(24) + print *,'dFSEE0 = ',dFSEE0 , cm(25) + + print *,'Damping = ',Damping , cm(26) + print *,'Param1 = ',Param1 , cm(27) + print *,'Param2 = ',Param2 , cm(28) + print *,'output_MTD = ',output_MTD , cm(29) + + print *,'timeste = ',timeste , cm(30) +c! cm(31) and cm(32) not used ? + print *,'Activation = ',Activation , cm(33) + print *,'target_l_CE = ',target_l_CE, cm(34) + print *,'kp = ',kp , cm(35) + print *,'kd = ',kd , cm(36) + + print *,'dotlCEdelay = ',dotlCEdelay, cm(37) + print *,'t_PreSim = ',t_PreSim , cm(38) + print *,'threshold = ',threshold , cm(39) + ENDDO + ENDIF +C----------------------------------------------- +C Muscle constant parameters calculation +C----------------------------------------------- +C +C ----------------------------------- +C Loop over elements +C ----------------------------------- + DO I = 1,NEL +C----------------------------------------------- +C Calc length element +C----------------------------------------------- +C + l_CE_opt = cm(10)*UVAR(16,I) !relative elongation of contractile element + l_SEE_0 = (1-cm(10))*UVAR(16,I) + l_PEE0 = cm(19)*l_CE_opt + d_SE_max = cm(27)*(cm(9)*cm(15))/(l_CE_opt*cm(16)) + l_SEE_nll = (1.0+cm(23))*l_SEE_0 + v_SEE = cm(23)/cm(24) + K_SEE_nl = cm(25)/(cm(23)*l_SEE_0)**v_SEE + K_SEE_l = cm(25)/(cm(24)*l_SEE_0) + K_PEE = cm(21)*(cm(9)/(l_CE_opt*(cm(11) + 1.0-cm(19)))**cm(20)) +! elleng = sqrt((x(1,ix1(i))-x(1,ix2(i)))**2 +! 1 +(x(2,ix1(i))-x(2,ix2(i)))**2 +! 2 +(x(3,ix1(i))-x(3,ix2(i)))**2)+cm(8) + elleng = XL(I) + cm(8) +c muscle length is element length plus offset + IF ((DT.EQ.0).OR.(ncycle.LE.1)) then + q = cm(3) + UVAR(3,I) = q + UVAR(13,I) = 0.0 + UVAR(15,I) = 0.0 + UVAR(16,I) = elleng +C RETURN + END IF +C----------------------------------------------- +C Calculation initial muscle force eqilibrium +C----------------------------------------------- + IF ((ncycle.LT.5).and.(UVAR(13,I).EQ.0)) THEN +C! UVAR(13,I)=anint(cm(37)/dt1) + dot_l_MTC = 0.0 + q = cm(3) + l_CE = ZEROIN(elleng ,l_CE_opt,l_SEE_0,q,cm, + 1 l_SEE_nll,K_SEE_nl,v_SEE,l_PEE0,K_PEE,K_SEE_l) + + CALL hill_model(l_CE,l_CE_opt,l_SEE_0,elleng, dot_l_MTC, q, cm, + 1 F_MTC, F_SEE, F_SDE, F_CE, F_PEE, F_isom, DT, UVAR, + 2 dot_l_CE, tt, i, l_PEE0, K_PEE, d_SE_max, l_SEE_nll, v_SEE, K_SEE_nl, K_SEE_l, epsilon, NUVAR) + END IF +C! +C----------------------------------------------- +C Calculation of delayed lCE/DOtlCE +C----------------------------------------------- +c ncycle.EQ.3, because for ncycle<=2: dt1==0 in case of part_averaged +c cm(37) - delay +c cm(33) - controller higher than Alpha + IF (cm(33).ge.1.0) THEN +c all controllers + IF (cm(37).EQ.0.0) THEN +c all controllers without delay: +c lCEdelay=lCE dotlCEdelay=dot_l_CE + UVAR(20,I) = UVAR(10,I) + UVAR(21,I) = UVAR(12,I) + ELSE +c controller with delay used --> ringbuffer needed + CALL ringbuffer(cm,UVAR,ncycle,tt,NUVAR,I) + ENDIF + lCEdelay = UVAR(20,I) + dotlCEdelay = UVAR(20,I) + ENDIF +C----------------------------------------------- +C Calculation of activation +C---------------------------------------------- +C + IF (cm(1).EQ.0) THEN +c! Use Activation Values directly + IF (cm(2).LT.0.0) THEN + q = GET_U_FUNC(IFUNC1,tt,DXDY) + ELSE + q = abs(cm(2)) + ENDIF + IF (q.GT.1.0) THEN + q = 1.0 + ELSE IF (q.LT.0.0) THEN + q = 0.0 + END IF + IF (int(cm(33)).EQ.0) THEN + UVAR(3,I) = q + ELSE + IF (tt.LT.cm(38)) THEN + UVAR(3,I) = q + ELSE + cm(1) = 3 + STIM = calc_STIM(cm,UVAR,tt,lCEdelay,dotlCEdelay,ncycle,NUVAR,I,IFUNC1,IFUNC2,IFUNC3,GET_U_FUNC) + q = calc_q(cm,l_CE_opt,UVAR,tt,STIM,DT,NUVAR,I,GET_U_FUNC,IFUNC2) + END IF + ENDIF + ELSE + STIM = calc_STIM(cm,UVAR,tt,lCEdelay,dotlCEdelay,ncycle,NUVAR,I,IFUNC1,IFUNC2,IFUNC3,GET_U_FUNC) + q = calc_q(cm,l_CE_opt,UVAR,tt,STIM,DT,NUVAR,I,GET_U_FUNC,IFUNC2) + ENDIF +C----------------------calc dot_l_MTC------------ +C IF (ncycle.GT.1) THEN + l_CE = UVAR(10,I) + dot_l_MTC = VX(I) + UVAR(11,I) = dot_l_MTC +C------------------------------------------------ + CALL hill_model(l_CE,l_CE_opt,l_SEE_0,elleng, dot_l_MTC, q, cm, + 1 F_MTC, F_SEE, F_SDE, F_CE, F_PEE, F_isom, DT, UVAR, + 2 dot_l_CE, tt, i, l_PEE0, K_PEE, d_SE_max, l_SEE_nll, v_SEE, K_SEE_nl, K_SEE_l, epsilon, NUVAR) +C ENDIF + ! Problably have to allocate new value to FX(I) + DX = DT * dot_l_MTC + elleng = UVAR(9,I) + DX +C + FX(I) = F_MTC +C----------------------------------------------- + UVAR(4,I)= F_MTC + UVAR(7,I)= F_SEE + UVAR(8,I)= F_SDE + UVAR(5,I)= F_CE + UVAR(6,I)= F_PEE + UVAR(14,I)= F_isom + UVAR(9,I)= elleng + UVAR(10,I)= l_CE + UVAR(12,I) = dot_l_CE +C---------------------------------------------- + IF (ncycle/100 == (REAL(ncycle)/100)) THEN +C print*, "F_MTC = ", F_MTC +C print*, "F_SEE = ", F_SEE +C print*, "F_SDE = ", F_SDE +C print*, "F_CE = ", F_CE +C print*, "F_PEE = ", F_PEE +C print*, "F_isom = ", F_isom +C print*, "elleng = ", elleng +C print*, "l_CE = ", l_CE +C print*, "dot_l_CE = ", dot_l_CE +C print*, "dot_l_MTC = ", dot_l_MTC +C print*, "q = ", q + ENDIF +C----------------------------------------------- +C IF (int(cm(29)).EQ.1.0) THEN +C IF (ncycle/100 == (REAL(ncycle)/100)) THEN +C OPEN(UNIT=10, FILE=FILENAME, STATUS='OLD', ACTION='WRITE', +C 1 POSITION='APPEND') +C IF (IO_STATUS /= 0) THEN +C PRINT *, 'Error: Cannot open file ', FILENAME +C STOP +C END IF +C WRITE(10,'(F10.4," ",F10.4," ",F10.4," ",F10.4," ",F10.4," ",F10.4," ",F10.4," ",F10.4," ",F10.4" ",F10.4)') tt, +C 1 UVAR(10,I)-UVAR(9,I), UVAR(10,I), UVAR(14,I), +C 2 UVAR(5,I), UVAR(7,I), UVAR(8,I),q,STIM,UVAR(4,I) +C CLOSE(UNIT=10) +C ENDIF +C ENDIF +C----------------------------------------------- +C *** TIMESTEP *** +C============================================ +C TIME STEP +C============================================ +C + STIFM(I) = ELASTIF + STIFR(I) = ZERO + VISCM(I) = ZERO + VISCR(I) = ZERO + XINER(I) = ZERO + + END DO +C----------------------------------------------- + RETURN + END subroutine ruser30 +C----------------------------------------------- +C----------------------------------------------- +C hill_model SUBROUTINE +C----------------------------------------------- + SUBROUTINE hill_model(l_CE,l_CE_opt,l_SEE_0,elleng,dot_l_MTC,q,cm, + 1 F_MTC,F_SEE, F_SDE ,F_CE,F_PEE,F_isom,DT,UVAR,dot_l_CE,tt,I,l_PEE0, + 2 K_PEE,d_SE_max,l_SEE_nll,v_SEE,K_SEE_nl,K_SEE_l,epsilon,NUVAR) +c + DOUBLE PRECISION, dimension(39) :: cm + DOUBLE PRECISION :: UVAR(NUVAR,*) + integer I,NUVAR + DOUBLE PRECISION A1,B1,dot_h,dot_h_0,ent_rate,mech_eff + DOUBLE PRECISION l_PEE0, K_PEE, d_SE_max, l_SEE_nll, v_SEE, K_SEE_nl + DOUBLE PRECISION l_CE, F_isom, F_PEE, F_SEE, l_SEE, L_A_rel, A_rel + DOUBLE PRECISION B_rel,Q_Brel,D0,C2,C1,C0,dot_l_CE + DOUBLE PRECISION v_max, P, O1, O2, O3, F_CE, F_CE_init, F_SDE, F_MTC + DOUBLE PRECISION dot_l_MTC, elleng, K_SEE_l, Q_Arel, F_SUM, AREA + DOUBLE PRECISION l_MTC_0,DT,tt,epsilon, Power, dot_l_SDE,q + DOUBLE precision :: calc_F_isom, calc_F_PEE, calc_F_SEE,l_CE_opt,l_SEE_0 + !common /coeff/ l_PEE0,K_PEE,d_SE_max,l_SEE_nll,v_SEE,K_SEE_nl + !1 ,K_SEE_l,epsilon +c Isometric Force +c + F_isom = calc_F_isom(l_CE,l_CE_opt,cm) +c +c Force of the parallel elastic element PEE + F_PEE = calc_F_PEE(l_CE, cm, l_PEE0, K_PEE) +c +c Force of the serial elastic element SEE +c + F_SEE = calc_F_SEE(l_CE,l_SEE_0,cm,elleng,l_SEE_nll,v_SEE,K_SEE_nl,K_SEE_l) +c +c Hill Parameters concentric contraction +c + if (l_CE.LT.l_CE_opt) then + L_A_rel = 1.0 + else + L_A_rel = F_isom + end if + Q_Arel = (0.25)*(1.0+3.0*q) + A_rel = cm(15)*L_A_rel*Q_Arel + Q_Brel = (1.0/7.0)*(3.0+4.0*q) + B_rel = cm(16)*Q_Brel +C CE CONTRACTION VELOCITY CALCULATION + call calc_damping(q,l_CE_opt,cm,A_rel,B_rel,F_SEE,F_PEE,F_isom,dot_l_MTC,D0,C2,C1, + 1 C0,l_PEE0,K_PEE,d_SE_max) +c + DD = C1**2.0-4.0*C2*C0 + IF ((abs(C2).GT.0.0).and.(DD.GE.0.0)) THEN + dot_l_CE = (-C1-sqrt(DD))/(2.0*C2) + ELSE IF (abs(C1).GT.0.0) THEN + dot_l_CE = -C0/C1 + ELSE + dot_l_CE = 0.0 + END IF +C ECCENTRIC CONTRACTION CASE + IF (dot_l_CE.GT.0.0) THEN +C + B_rel = (q*F_isom*(1.0-cm(18))/(q*F_isom+A_rel)*B_rel/cm(17)) + A_rel = -cm(18)*q*F_isom +C +C CE ECCENTRIC CONTRACTION VELOCITY CALCULATION + call calc_damping(q,l_CE_opt,cm,A_rel,B_rel,F_SEE,F_PEE,F_isom,dot_l_MTC,D0,C2,C1, + 1 C0,l_PEE0,K_PEE,d_SE_max) +c + DD = C1**2.0-4.0*C2*C0 +c + IF ((abs(C2).GT.0.0).and.(DD.GE.0.0)) THEN + dot_l_CE = (-C1+sqrt(DD))/(2.0*C2) + ELSE IF (abs(C1).GT.0.0) THEN + dot_l_CE = -C0/C1 + ELSE + dot_l_CE = 0.0 + ENDIF +C +C IF (dot_l_CE.LT.0.0) THEN +C PRINT*, 'eccentric case error in muscle' +C ENDIF +c + ENDIF +C NUMERICAL RESTRICTION CHECK +C + IF (l_CE.LT.(0.001*l_CE_opt)) THEN + dot_l_CE = dot_l_MTC + ELSE IF (l_CE.GT.(1.999*l_CE_opt)) THEN + l_SEE = elleng-l_CE + IF (l_SEE.GT.l_SEE_nll) THEN + dot_l_CE = dot_l_MTC/(1.0+(K_PEE*cm(20)*(l_CE-l_PEE0)** + 1 (cm(20)-1.0))/K_SEE_l) + ELSE + dot_l_CE = dot_l_MTC/(1.0+(K_PEE*cm(20)*(l_CE-l_PEE0)** + 1 (cm(20)-1.0))/(K_SEE_nl*v_SEE*(l_SEE-l_SEE_0)** + 2 (v_SEE-1.0))) + ENDIF + ENDIF +c +c contractile element force +c +c overwrite existing and (possibly) modyfied A_rel and B_rel + A_rel = cm(15)*L_A_rel*Q_Arel + B_rel = cm(16)*Q_Brel +c + IF (dot_l_CE.GT.0.0) THEN + B_rel = (q*F_isom*(1.0-cm(18))/(q*F_isom+A_rel)*B_rel/cm(17)) + A_rel = -cm(18)*q*F_isom + ENDIF + F_CE = cm(9)*(((q*F_isom+A_rel)/(1.0-dot_l_CE/(l_CE_opt*B_rel)))-A_rel) +c + IF (F_CE.LT.0.0) THEN + F_CE = 0.0 + ENDIF +c +c +c force of the serial damping element +c + F_SDE = d_SE_max*((1.0-cm(28))*((F_CE+F_PEE)/cm(9))+cm(28)) + 1 *(dot_l_MTC-dot_l_CE) +c +c F_SDE force limit start + SELECT CASE (int(cm(26))) + CASE (1) + IF (F_SDE.LT.-10*abs(F_SEE)) THEN + F_SDE = -10*abs(F_SEE) + ENDIF + END SELECT +c F_SDE force limit stop +c + F_MTC = F_SEE+F_SDE +c +c calc l_CE (integrate dot_l_CE with finite difference methode) +cs + l_CE = l_CE + dot_l_CE * DT + UVAR(10,I) = l_CE + UVAR(12,I) = dot_l_CE +c + RETURN + END SUBROUTINE hill_model +C----------------------------------------------- +C calc_damping SUBROUTINE +C---------------------------------------------- +C + SUBROUTINE calc_damping(q,l_CE_opt,cm,A_rel,B_rel,F_SEE,F_PEE,F_isom,dot_l_MTC,D0,C2,C1, + 1 C0,l_PEE0,K_PEE,d_SE_max) + DOUBLE PRECISION,dimension(39) :: cm + DOUBLE PRECISION :: q,A_rel,B_rel,F_SEE,F_PEE,F_isom + DOUBLE PRECISION :: dot_l_MTC,d_SE_max,D0,C2,C1,C0,l_PEE0,K_PEE,l_CE_opt +C +C! 3. Force-dependent serial (SE)-Damping + D0 = l_CE_opt*B_rel*d_SE_max*(cm(28)+(1.0-cm(28))*(q*F_isom+F_PEE/cm(9))) + C2 = d_SE_max*(cm(28)-(A_rel-F_PEE/cm(9))*(1.0-cm(28))) + C1 = -C2*dot_l_MTC-D0-F_SEE+F_PEE-cm(9)*A_rel + C0 = D0*dot_l_MTC+l_CE_opt*B_rel*(F_SEE-F_PEE-cm(9)*q*F_isom) + RETURN + END +C +C----------------------------------------------- +C ring buffer subroutine +C---------------------------------------------- +C + SUBROUTINE ringbuffer(cm,UVAR,ncycle,tt,NUVAR,I) +c + DOUBLE PRECISION cm(*),UVAR(NUVAR,I) + INTEGER ncycle,I,NUVAR + DOUBLE PRECISION lCEdelay,dotlCEdelay,tt + DOUBLE PRECISION delaylCE, delaydotlCE +c +c + if (UVAR(31,I).EQ.0) then +c ----------- ringbuffer initialisation ---------------- +c introducing Ring-Buffer with starting index - beginlCE 150 +c because of lsdyna messing with our data at lower index + UVAR(31,I) = 150 +c total ring-buffer size based on NHISVAR (FLOOR = "abrunden") + UVAR(30,I) = floor((NUVAR-UVAR(31,I)-4)/2) +c NHISVAR is the amount of maximum hsv (defined in compile process) +c no_hsvs is the user defined amount of hsvs (!<=NHISVAR) +c +c indexr1 + UVAR(32,I) = 0 +c indexr2 + UVAR(33,I) = 1 +c lasttime + UVAR(34,I) = 0 +c begindotlCE + UVAR(36,I) = UVAR(31,I)+UVAR(30,I)+2 +c indexdotr1 + UVAR(37,I) = 0 +c indexdotr2 + UVAR(38,I) = 1 +c dotlasttime + UVAR(39,I) = 0 +c +c init lCEdelay and dotlCEdelay + lCEdelay = 0 + dotlCEdelay = 0 +c + ELSE +c ring buffer already initialised + SELECT CASE (int(cm(33))) + CASE (1:3) +c lambda or hybrid controller + lCEdelay = delaylCE(ncycle,UVAR,cm,tt,NUVAR,I) + dotlCEdelay = delaydotlCE(ncycle,UVAR,cm,tt,NUVAR,I) + CASE (4) +c reflex controller only delaylCE needed + lCEdelay = delaylCE(ncycle,UVAR,cm,tt,NUVAR,I) + dotlCEdelay = UVAR(12,I) + CASE DEFAULT + END SELECT + ENDIF + UVAR(30,I) = lCEdelay + UVAR(21,I) = dotlCEdelay + RETURN + END SUBROUTINE ringbuffer +c \ No newline at end of file diff --git a/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Example/Gunther_1D_0000.rad b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Example/Gunther_1D_0000.rad new file mode 100644 index 0000000..e061117 --- /dev/null +++ b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Example/Gunther_1D_0000.rad @@ -0,0 +1,208 @@ +#RADIOSS STARTER +##======================================================================================== +## +# MIT License +# +# Copyright (c) 2025 RM412 +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# Suggestion : +# +# BSD 3-Clause License +# Copyright (c) 2014 +# D. Haeufle, M. Günther, A. Bayer, S. Schmitt +# Copyright (c) 2017 +# C. Kleinbach, O. Martynenko, J. Promies, D. F. B. Haeufle, J. Fehr, S. Schmitt +# Copyright (c) 2019 +# O. Martynenko, F. Kempter, C. Kleinbach, S. Schmitt, J. Fehr +# Copyright (c) 2020 +# P. Lerge, O. Martynenko, L. Nölle, S. Schmitt, I. Wochner +# Copyright (c) 2022 +# O. Martynenko, F. Kempter, C. Kleinbach, L. Nölle, P. Lerge, S. Schmitt, J. Fehr +# Copyright (c) 2026 +# A. Ardisson, E. Wagnac, M.-H. Beauséjour/ÉTS +# Adaptation and port to OpenRadioss +# All rights reserved. +# Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following C # conditions are met: +# Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# +# Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# +# Neither the name of the copyright holder nor the names of its contributors +# may be used to endorse or promote products derived from this software +# without specific prior written permission. +# +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +##======================================================================================== +## +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/PARAMETER/GLOBAL/REAL/1 + +param1 6.0 +## +/BEGIN +#Runname +Gunther_100g_V5 +# Invers Irun + 2022 0 +# Input_f_unit Input_length_unit Input_time_unit + g mm ms +# Work_mass_unit Work_length_unit Work_time_unit + g mm ms +## +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/NODE +# node_ID Xc Yc Zc + 1 0.0 60.0 0.0 + 2 0.0 0.0 0.0 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/SPRING/1 +# sprg_ID node_ID1 node_ID2 skew_ID + 1 1 2 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/PART/1 +#part_title +hill_element +# prop_ID mat_ID Thick + 2 0 0 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +#HWCOLOR properties 260003151 10 +/PROP/USER2/2 +USER_SPRING +# Rho Diameter Contact Stif ActOpt SENS_ID + 0.001 12.0 8.0 2 17 +# STIM_ID q0 tauq/c betaq/eta k + 3.0 5.0E-3 1.373E-4 5.27E+4 2.9 +# m muscle len offset Fmax lCEopt dWdes + 11.3E-3 0.0 30 0.25 0.14 +# nuCEdes dWasc nuCEasc Arel0 Brel0 + 3.0 0.57 4.0 0.1 1.0E-3 +# Secc Fecc LPEE0 nuPEE FPEE + 2.0 1.8 0.9 2.5 1.0 +# dUSEEnll duSEEl dFSEE0 + 0.1825 0.073 60.0 +# Damping Param1 Param2output_MTD + 3.0 0.3 0.01 1 +# timestep Activation target_l_CE kp kd + 0.0 0.75 6 +# dotlCEdelay t_PreSim threshold + +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +##HWCOLOR solvermasses 1 6 +/ADMAS/0/1 +100g_MASS +# MASS grnd_ID + 100 2 +##HMNAME LOADCOLS 1EMBEDDED +##HWCOLOR loadcollectors 1 21 +/BCS/1 +Loads +# Tra rot skew_ID grnod_ID + 111 111 0 1 +/BCS/2 +Loads +#FREE TRANSLATION IN Y +# Tra rot skew_ID grnod_ID + 101 111 0 2 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +##HMNAME LOADCOLS 2GRAV +##HWCOLOR loadcollectors 2 55 +/GRAV/2 +champ_GRAV +#funct_IDT DIR skew_ID sensor_ID grnod_ID Ascale_x Fscale_Y + 19 Y 0 0 3 -1.0 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +##HWCOLOR curves 19 5 +/FUNCT/19 +#title +CHMP_GRAV +# X Y + 0.0 0.00981 + 2000.0 0.00981 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/FUNCT/1 +aa + 0.0 0.0 + 10.0 0.0 + 100.0 0.0 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/FUNCT/2 +tauq_vs_time + 0.0 0.0 + 10.0 20.0 + 100.0 20.0 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/FUNCT/3 +Activaton_Gunther +# X Y + 0.0 0.0 + 1399.9 0.0 + 1400.0 1.0 + 2499.0 1.0 + 2500.0 0.0 + 3000.0 0.0 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/GRNOD/NODE/1 +#title +EMBEDDED_NODE +# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10 + 1 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/GRNOD/NODE/2 +#title +LOADING_NODE +# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10 + 2 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/GRNOD/NODE/3 +#title +MASS_NODE +# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10 + 2 +#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----| +/TH/NODE/1 +#thgroup_name +point_suivi +# var_ID1 var_ID2 var_ID3 var_ID4 var_ID5 var_ID6 var_ID7 var_ID8 var_ID9 var_ID10 +DY V DY VY +# node_ID skew_ID node_name + 2 0 +##-------------------------------------------------------------------------------------------------- +## End Of Radioss Block Deck +##-------------------------------------------------------------------------------------------------- +/END + diff --git a/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Example/Gunther_1D_0001.rad b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Example/Gunther_1D_0001.rad new file mode 100644 index 0000000..595a0eb --- /dev/null +++ b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Example/Gunther_1D_0001.rad @@ -0,0 +1,22 @@ +/RUN/Gunther_1D/1/ +2000.00000000000 +/STOP +0.000000000000000 0.000000000000000 0.000000000000000 1 1 0 +/MON/FULL +/PRINT/-10/0 +/VERS/2023 +/DT/NODA/CST/0 +0.500000000000000 0.000000000000000 0.000000000000000 +/TFILE/0 +20.00000000000000 +/ANIM/DT +0.000000000000000 200.0000000000000 +/ANIM/ELEM/EPSP +/ANIM/ELEM/VONM +/ANIM/NODA/DMAS +/ANIM/SHELL/TENS/STRESS/UPPER +/ANIM/SHELL/TENS/STRESS/LOWER +/ANIM/SHELL/EPSP/UPPER +/ANIM/SHELL/EPSP/LOWER +/ANIM/VECT/VEL +/ANIM/SPRING/FORCE \ No newline at end of file diff --git a/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Git_figure.jpg b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Git_figure.jpg new file mode 100644 index 0000000..1c8ac3d Binary files /dev/null and b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Git_figure.jpg differ diff --git a/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/LICENSE b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/LICENSE new file mode 100644 index 0000000..3bde7ef --- /dev/null +++ b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/LICENSE @@ -0,0 +1,66 @@ +MIT License + +Copyright (c) 2025 RM412 + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + +Suggestion : + +BSD 3-Clause License +Copyright (c) 2014 + D. Haeufle, M. Günther, A. Bayer, S. Schmitt +Copyright (c) 2017 + C. Kleinbach, O. Martynenko, J. Promies, D. F. B. Haeufle, J. Fehr, S. Schmitt +Copyright (c) 2019 + O. Martynenko, F. Kempter, C. Kleinbach, S. Schmitt, J. Fehr +Copyright (c) 2020 + P. Lerge, O. Martynenko, L. Nölle, S. Schmitt, I. Wochner +Copyright (c) 2022 + O. Martynenko, F. Kempter, C. Kleinbach, L. Nölle, P. Lerge, S. Schmitt, J. Fehr +Copyright (c) 2026 +A. Ardisson, E. Wagnac, M.-H. Beauséjour/ÉTS + Adaptation and port to OpenRadioss +All rights reserved. +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + +Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + +Neither the name of the copyright holder nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + POSSIBILITY OF SUCH DAMAGE. + diff --git a/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Muscle setup.jpeg b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Muscle setup.jpeg new file mode 100644 index 0000000..b9e36d6 Binary files /dev/null and b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/Muscle setup.jpeg differ diff --git a/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/README.md b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/README.md new file mode 100644 index 0000000..54e84a0 --- /dev/null +++ b/Components/Materials/Extended-Hill-Type-Muscle-1D-user-spring/README.md @@ -0,0 +1,247 @@ +# Extended-Hill-Type-Muscle--1D-user-spring-for-OpenRadioss +# Introduction +This project is a derivative work of the original LS-DYNA implementation by Haeufle (2014) et al. and Martynenko et al. (2023), redistributed under the BSD 3-Clause License. The model is based on the Enhanced Hill-Type Model (EHTM) and has been translated, refactored, and validated to operate within the OpenRadioss material framework. +The aim is to give FEM users a ready-to-use muscle contraction material law that can be integrated into larger biomechanical simulations, such as head–neck models for high-impact accident studies. + +The model represents a single muscle–tendon unit and captures the main mechanical and physiological behaviors of active muscle tissue: +- Contractile element (CE) for active force generation +- Parallel elastic element (PEE) for passive stiffness +- Series elastic element (SEE) for tendon elasticity +- Serial damping element (SDE) adds a velocity-dependent force term in series with the CE and SEE, improving stability and making the model better at reproducing rapid loading/unloading responses + +All parameters are taken from experimental data and literature Günther et al. (2007), (Kleinbach et al., 2017), (Martynenko et al., 2023), ensuring that the model is both physiologically realistic and numerically stable. +Example input files, Fortran implementation, and plots of the force–length and force–velocity relationships are provided for reference. + +# Options and Keywords used +# Model description +The EHTM encapsulates the physiomechanical behavior of a muscle–tendon unit via four sequential components: + +### 1. Contractile Element CE +- Generates active force based on excitation and activation dynamics, modulated by force–length and force–velocity relationships. +### 2. Series Elastic Element SEE +- Models tendon elasticity using a piecewise law: a non-linear region followed by a linear response at higher elongations. +### 3. Parallel Elastic Element PEE +- Models passive force in parallel with CE +- Activated only when CE length exceeds LPEE,0 +- uses an exponential stiffness law +### 4. Serial Damping Element SDE +- Adds velocity-dependent damping to stabilize the system under dynamic loading. + +## Muscle Activation Dynamics + +This model simulates the activation of a single muscle fiber using a step-type stimulation (STIM) signal as the input. The stimulation represents the neural trigger from the experiment and is converted into muscle activation through a chemo-physiological process. + +The model includes two calcium ion diffusion methods: Zajac (1989) and Hatze (1978). In this implementation, the Hatze method is used, following Kleinbach et al. (2017). The Hatze approach calculates the muscle activation level q based on: +- The muscle activation $q$ depends on the contractile element length $l_{CE}$ and calcium ion concentration $\gamma_{\mathrm{rel}}$. + + +$$ +q_{\text{rel}} = q_0 + \frac{(\rho \, \gamma_{\text{rel}})^3}{1 + (\rho \, \gamma_{\text{rel}})^3} +$$ + + +$$ +\frac{d\gamma_{\mathrm{rel}}}{dt} = m\,(\mathrm{STIM} - \gamma_{\mathrm{rel}}) +$$ + + + + +$$ +\rho = c \, \eta^{(k - l_{CE,\text{rel}})} +$$ + + +#### Muscle Model Setup +![EHTM Muscle Model Setup](Git_figure.jpg) + +## Mechanical Contraction Model + +The active force of the contractile unit depends on the activation level *q*, the unit’s length $l_{CE}$, and its contraction velocity (dlCE/dt). +Fmax and the Hill constants $A_{rel}$ and $B_{rel}$, which shape the force–velocity and force–length relationships. + +$$ +F_{CE}\left(l_{CE}, \dot{i}_{CE}, q\right) = F_{\max}\left(\frac{q F_{\mathrm{isom}} + A_{\mathrm{rel}}}{1 - \dfrac{\dot{i}_{CE}}{B_{\mathrm{rel}}\, l_{CE,\mathrm{opt}}}}- A_{\mathrm{rel}}\right) +$$ + +The isometric force depends solely on 𝑙𝐶𝐸, following a bell-shaped curve whose slope and width are controlled by ∆𝑊𝑙𝑖𝑚𝑏 and 𝜈𝐶𝐸,𝑙𝑖𝑚𝑏. + +$$ +F_{isom}(l_{CE}) = \exp \left( - \left| \frac{\frac{l_{CE}}{l_{CE,opt}} - 1}{\Delta W_{limb}} \right|^{\nu_{CE,limb}} \right) +$$ + + +The passive mechanical property of muscle fibres where K and v are stiffness components: + +$$ +F_{PEE} = +\begin{cases} +0, & l_{CE} < l_{PEE,0} \\ +K_{PEE} \left( l_{CE,opt} - l_{PEE,0} \right)^{\nu_{PEE}}, & l_{CE} \ge l_{PEE,0} +\end{cases} +$$ + +The second part of the EHTM models tendons, using a nonlinear elastic element and an adaptive damping component to reflect their fibrous nature. This approach captures the initial nonlinear response caused by collagen fiber reorganization, followed by a linear stiffness phase. + +$$ +F_{SEE}(l_{SEE}) = +\begin{cases} +0, & \text{if } l_{SEE} < l_{SEE,0} \\ +K_{SEE,nl} \left( l_{SEE} - l_{SEE,0} \right)^{\nu_{SEE}}, & \text{if } l_{SEE} \ge l_{SEE,0} \text{ and } l_{SEE} < l_{SEE,nll} \\ +\Delta F_{SEE,0} + K_{SEE,l} \left( l_{SEE} - l_{SEE,nll} \right), & \text{if } l_{SEE} \ge l_{SEE,nll} +\end{cases} +$$ + + +In the EHTM model, the SDE damping element is adaptive, using a minimum damping value $R_{SDE}$ and a non-dimensional scale factor $r_{SDE}$ to reduce high-frequency instability. + +$$ +F_{SDE}\left(l_{CE}, \dot{l}_{SDE}, q\right)=d_{SDE,\max}\left[\left(1 - R_{SDE}\right)\frac{F_{CE} + F_{PEE}}{F_{\max}}+ R_{SDE}\right]\dot{l}_{SDE} +$$ + +$$ +d_{SDE,\max} = D_{SDE}\frac{F_{\max}\,A_{\mathrm{rel},0}}{l_{CE,\mathrm{opt}}\,B_{\mathrm{rel},0}} +$$ + +# Boundary Conditions + +This 1D muscle–tendon unit is exercised in pure axial loading: + +- **Nodes & element** + - `NODE 1` (origin/proximal end) — fixed support + - `NODE 2` (distal end/loading point) — moves only along the muscle axis (Y) + - One user spring (EHTM) connects `NODE 1` → `NODE 2` + +- **Kinematic constraints** + - `NODE 1`: fully constrained (translations + rotations) + - `NODE 2`: X and Z fixed; Y free (axial motion) + +- **Loading** + - A **point mass of (100 g)** is attached at the loading node + - **Gravity** acts in −Y using a constant time function (9.81 m/s² in the working unit system) + +- **What to monitor** + - Displacement **DY** and velocity **VY** of `NODE 2` via `/TH/NODE` + + +## OpenRadioss keywords + +```text +# --- Node groups (example) --------------------------------- +# Group 1: fixed node (NODE 1) +# Group 2: loading node (NODE 2) +# Group 3: gravity mass node group (usually same as Group 2) + +# --- Boundary conditions ----------------------------------- +# Fix NODE 1 (all translations + rotations) +#/BCS/1 +# Tra Rot skew_ID grnod_ID + 111 111 0 1 + +# Lock X,Z and free Y on NODE 2 → translation mask 101 (X,Z fixed; Y free) +#/BCS/2 +# Tra Rot skew_ID grnod_ID + 101 111 0 2 + +# --- Added mass at loading node (0.1 kg = 100 g) ----------- +#/ADMAS/0/1 +# MASS(g) grnod_ID + 100 2 + +# --- Gravity in −Y ----------------------------------------- +# Constant gravity defined by function 19 +#/GRAV/2 +#funct_IDT DIR skew_ID sensor_ID grnod_ID Ascale_X Fscale_Y + 19 Y 0 0 3 -1.0 + +# Function 19: constant acceleration +#/FUNCT/19 +# X(time) Y(accel) + 0.0 0.00981 + 2000 0.00981 + +# --- Time histories at NODE 2 ------------------------------- +#/TH/NODE/1 +# var_ID1 var_ID2 + DY VY +# node_ID skew_ID + 2 0 + +``` +# Material characterization +The parameters below define the single-fiber Extended Hill-Type Muscle (EHTM) used in this repo. +They come from the muscle-scale validation and literature synthesis. + +### Activation (Hatze-type) +| Parameter | Symbol | Value | Units | Notes | +| ------------------ | ------ | -------: | ----- | -------------------------------------------------------------------- | +| Minimum activation | $q_0$ | 5.0e-3 | – | Baseline activation level | +| Calcium gain | $c$ | 1.373e-4 | mol/L | Hatze parameter | +| Calcium scaling | $\eta$ | 5.27e+4 | L/mol | Hatze parameter | +| Exponent | $k$ | 2.9 | – | Hatze parameter | +| Relaxation rate | $m$ | 11.3e-3 | 1/ms | $\frac{d\gamma_{\mathrm{rel}}}{dt} = m\(\mathrm{STIM} - \gamma_{\mathrm{rel}})$ | + +### Isometric Force - Length +| Parameter | Symbol | Value | Units | Notes | +| --------------------- | ------------------------- | ----: | ----- | ---------------------- | +| Max. isometric force | $F_{\max}$ | 30 | N | Peak active force | +| Optimal fiber length | $l_{CE,\mathrm{opt}}$ | 15 | mm | At $F_{\max}$ | +| Width (descending) | $\Delta W_{\mathrm{des}}$ | 0.14 | – | F-L curve width (desc.) | +| Exponent (descending) | $\nu_{CE,\mathrm{des}}$ | 3.0 | – | F-L shape (desc.) | +| Width (ascending) | $\Delta W_{\mathrm{asc}}$ | 0.57 | – | F-L curve width (asc.) | +| Exponent (ascending) | $\nu_{CE,\mathrm{asc}}$ | 4.0 | – | F-L shape (asc.) | + +### Force - Velocity, Hill Parabola +| Parameter | Symbol | Value | Units | Notes | +| --------------- | -------------------- | -----: | ----- | ----- | +| Hill constant A | $A_{\mathrm{rel},0}$ | 0.1 | – | | +| Hill constant B | $B_{\mathrm{rel},0}$ | 1.0e-3 | – | | + +### Parallel Elastic Element +| Parameter | Symbol | Value | Units | Notes | +| ----------------------- | ----------- | ----: | ----- | ---------------------------------------------------------------- | +| Activation length ratio | $L_{PEE,0}$ | 0.9 | – | PEE active when $l_{CE} \ge L_{PEE,0} \cdot l_{CE,\mathrm{opt}}$ | +| Stiffness exponent | $\nu_{PEE}$ | 2.5 | – | $F_{PEE}\propto(l_{CE}-l_{PEE,0})^{\nu_{PEE}}$ | +| PEE scale | $F_{PEE}$ | 1.0 | – | Dimensionless scale (maps to $K_{PEE}$ in code) | + +### Series Elastic Element +| Parameter | Symbol | Value | Units | Notes | +| -------------- | -------------------- | -----: | ----- | ----------------------- | +| Slack length | $l_{SEE,0}$ | 45 | mm | Zero-force length | +| Nonlinear span | $\Delta U_{SEE,nll}$ | 0.1825 | – | To transition point | +| Linear span | $\Delta U_{SEE,l}$ | 0.073 | – | Linear region scale | +| Force at SEE | $\Delta F_{SEE,0}$ | 60 | N | Offset at $l_{SEE,nll}$ | + +### Serial Damping Element +| Parameter | Symbol | Value | Units | Notes | +| ------------------ | -------------- | -------------------------------------------------------------------------------------: | ------- | ------------------- | +| Damping scale | $D_{SDE}$ | 0.3 | – | Sets $d_{SDE,\max}$ | +| Min. damping ratio | $R_{SDE}$ | 0.01 | – | | +| Max. damping | $d_{SDE,\max}$ | $D_{SDE}\,\dfrac{F_{\max} A_{\mathrm{rel},0}}{l_{CE,\mathrm{opt}} B_{\mathrm{rel},0}}$ | N·ms/mm | Derived from above | + +##### Units: This model uses mm and ms (with N for force). + +# References +The following publications provide the theoretical background and validation for the Extended Hill-Type Muscle (EHTM) model and related muscle material models used in this repository. + +1. **Kleinbach, C., Martynenko, O., Promies, J., Haeufle, D. F. B., Fehr, J., & Schmitt, S.** (2017). + *Implementation and validation of the extended Hill-type muscle model with robust routing capabilities in LS-DYNA for active human body models.* + BioMedical Engineering OnLine, 16(1), 109. + [https://doi.org/10.1186/s12938-017-0399-7](https://doi.org/10.1186/s12938-017-0399-7) + +2. **Martynenko, O., Kleinbach, C., Promies, J., Fehr, J., Haeufle, D. F. B., & Schmitt, S.** (2023). + *Development and verification of a physiologically based Hill-type muscle material model in LS-DYNA for Active Human Body Models.* + Computer Methods in Biomechanics and Biomedical Engineering, 26(14), 1630–1648. + [https://doi.org/10.1080/10255842.2023.2251766](https://doi.org/10.1080/10255842.2023.2251766) + + +3. **Günther, M., Röhrle, O., Schmitt, S., & Blickhan, R.** (2007). + *High-frequency oscillations as a consequence of neglected serial damping in Hill-type muscle models.* + Biological Cybernetics, 97(1), 63–79. + [https://doi.org/10.1007/s00422-007-0154-3](https://doi.org/10.1007/s00422-007-0154-3) + +4. **Haeufle, D. F. B., Günther, M., Bayer, A., & Schmitt, S.** (2014). + *Hill-type muscle model with serial damping and eccentric force–velocity relation.* + Journal of Biomechanics, 47(6), 1531–1536. + [https://doi.org/10.1016/j.jbiomech.2014.02.009](https://doi.org/10.1016/j.jbiomech.2014.02.009) +