免费可用于 LS DYNA 的开源高精度人体肌肉模型

这是用于商业有限元分析软件 LS-DYNA 的开源肌肉模型，该材料模型可用于科学目的的研究。

通用有限元（FE）仿真软件 LS-DYNA 中将扩展的具有连续阻尼，偏心力 - 速度关系，集成的 Ca2 + 依赖性激活动力学和路由功能的四元素 Hill 型肌肉模型用作桁架的用户材料元素。从文献中选取了三组不同的哺乳动物实验数据，对这种物质模型进行了验证和验证。

将它与 LS-DYNA 中已经存在的 * MAT_MUSCLE（* MAT_156）希尔型肌肉模型进行了比较，该模型目前在有限元人体模型（HBM）中使用。给出了一个从 FE ViVA OpenHBM 提取手臂模型的应用示例，其中考虑了生理肌肉路径。仿真结果表明材料模型精度更高，与 * MAT_156 相比，计算鲁棒性更高，并且具有更好的肌肉路由能力。

用户资料子例程的 Fortran 源代码 dyn21.f 这项研究中进行的所有模拟的肌肉参数均作为开源材料下的补充材料提供。这样可以在 LS-DYNA 中快速应用建议的材料模型，尤其是在主动人体模型（AHBM）中。

1. 版本支持

该代码是为 LS-DYNA R8.1.0 编写的，更准确地说是为 Usermat 软件包编写的：

LS-DYNA_smp_d_R810_x64_redhat511_ifort131

被使用了。您还可以检查 dyn21.f-file 的开头，该文件应为

> c
c $Id: dyn21.F 96707 2015-03-18 22:29:10Z ubasu $
> c

经过修改后，它也应该可以与其他版本的 LS-DYNA 一起编译。

2. 项目内容

该项目的程序包主要包含肌肉模型源代码和相关的使用案例。

源代码：

• 包含注释肌肉模型：dyn21.f
• 补丁文件：patchfile

版本 v2.0 的示例即将推出：

• 使用扩展山丘型肌肉模型的颈部肌肉 ViVA HBM 的跌落测试装置。运行 main.k 文件。

源代码的注释版本是完整 dyn21.f 文件的剥离版本。您必须在 dyn21.f LS-DYNA 发行商提供的文件的相应部分复制源代码。或者，您可以使用 patchfile 和著名的 unix 命令 patch 自动执行此操作。

在这两种情况下，您都必须手动添加根查找算法，请参见源代码的注释版本的第 505 行。我们使用 GE Forsythe 文中介绍的 Sec 中的 ZEROIN 函数：

GE Forsythe；马萨诸塞州马尔科姆；Moler，CB：用于数学计算的计算机方法。Prentice Hall 专业技术参考，1977 年。

3. 计算结构

模型采用 41 号材料，定义具有可收缩元素的扩展 Hill 型肌肉模型，平行弹性元件和串联弹性元件以及系列阻尼元件。

• 材料模型的定义：

subroutine umat41 (lft,cm,eps,sig,epsp,hsv,dt1,capa,etype,tt,
1 temper,failel,crv,cma,qmat,elsiz,idele,x,v,a,i,d1,d2,d3,NHISVAR,
2 no_hsvs)

• 历史变量

hsv(1)    = sig(1)
hsv(2)    = STIM
hsv(3)    = q
hsv(4)    = F_MTC
hsv(5)    = F_CE
hsv(6)    = F_PEE
hsv(7)    = F_SEE
hsv(8)    = F_SDE
hsv(9)    = elleng (l_MTC)
hsv(10)   = l_CE
hsv(11)   = dot_l_MTC
hsv(12)   = dot_l_CE
hsv(13)   = counter_output
hsv(14)   = F_isom
hsv(15)   = gam_rel
hsv(16)   = l_MTC_0 (initial element length)
hsv(20)   = lCEdelay
hsv(21)   = dotlCEdelay
hsv(22)   = l_CE_ref (for reflex controller)
hsv(23)   = strain (for reflex controller)
hsv(24)   = STIM_reflex_prev (STIM value of reflex controller in previous timestep)
hsv(30)   = buffersize
hsv(31)   = idx_begin_lCEbuffer (hsv(142:145) seems to be used internally by LSDYNA, suggestion: use hsv(31)=150)
hsv(32)   = indexr1……index1 of lce-ringbuffer (not eq index of hsv)
hsv(33)   = indexr1……index2 of lce-ringbuffer (not eq index of hsv)
hsv(34)   = lasttime
hsv(36)   = begindotlCE
hsv(37)   = indexdotr1……index1 of dotlce-ringbuffer (not eq index of hsv)    
hsv(38)   = indexdotr2……index2 of dotlce-ringbuffer (not eq index of hsv)  
hsv(39)   = dotlasttime   
hsv(hsv(31):hsv(31)+hsv(30)-1) = ringbuffer_l_CE
hsv(hsv(36):hsv(36)+hsv(30)-1) = ringbuffer_dot_l_CE

• 初始化变量

dimension cm(*),eps(*),sig(*),hsv(*),crv(lq1,2,*),cma(*)
 dimension x(3,*),v(3,*),a(3,*),qmat(3,3)
 dimension d1(*),d2(*),d3(*)
 logical failel
 character*5 etype
 real*8 l_PEE0, K_PEE, d_SE_max, l_SEE_nll, v_SEE, K_SEE_nl
 real*8 F_isom, F_PEE, F_SEE, l_SEE, L_A_rel, A_rel, q,dq
 real*8 B_rel, L_B_rel, Q_Brel, D0, C2, C1, C0, dot_l_CE
 real*8 O1, O2, O3, F_CE, F_CE_init, F_SDE, F_MTC
 real*8 dot_l_MTC,Q_Arel, K_SEE_l, F_SUM, AREA,l_CE,q_old
 real*8 tol,l_MTC_0,STIM,dSTIM,tau,dtau,epsilon
 real*8 delay, gam_rel, rho_act, lCEdelay, dotlCEdelay
 real*8 lambda,dlambda
 integer curve_id, var,k,j,idpart
 common /coeff/ l_PEE0,K_PEE,d_SE_max,l_SEE_nll,v_SEE,K_SEE_nl
1 ,K_SEE_l,epsilon
 common/soundloc/sndspd(nlq),sndsp(nlq),diagm(nlq),sarea(nlq)
 common/aux33loc/
1 ix1(nlq),ix2(nlq),ix3(nlq),ix4(nlq),ix5(nlq),mxt(nlq)
 common/aux14loc/
1 sig1(nlq),sig2(nlq),sig3(nlq),sig4(nlq),
2 sig5(nlq),sig6(nlq),epx1(nlq),epx2(nlq),aux(nlq,14),dig1(nlq),
3 sig_pass(nlq),sig_actv(nlq),act_levl(nlq),out_leng(nlq),
4 eps_rate(nlq),sig_svs(nlq),sig_sde(nlq)

• 计算肌肉参数

l_PEE0=cm(19)*cm(10)
d_SE_max=cm(27)*(cm(9)*cm(15))/(cm(10)*cm(16))
l_SEE_nll=(1.0+cm(23))*cm(22)
v_SEE=cm(23)/cm(24)
K_SEE_nl=cm(25)/(cm(23)*cm(22))**v_SEE
K_SEE_l=cm(25)/(cm(24)*cm(22))
K_PEE=cm(21)*(cm(9)/(cm(10)*(cm(11)+1.0-cm(19)))**cm(20))

后续完成“计算初始肌肉力量平衡”、“更新肌肉力 / 长度 / 速度”等过程，感兴趣可以参照第 4 小节的内容。