MULTI-SCALE MODELING AND SIMULATION OF SKELETAL MUSCLE BIOMECHANICAL PROPERTIES
-
摘要: 针对肌纤维微观结构模型与显微镜下观察的图像存在一定差异、微观组分生物力学模型无法有效捕获骨骼肌剪切变形时的力学行为、骨骼肌多尺度数值模型计算成本高的问题, 本文从实验、多尺度建模和仿真的角度研究骨骼肌被动力学特性, 提出以曲边泰森多边形作为肌纤维截面形状, 并建立骨骼肌微观尺度代表性体元RVE; 提出新的生物力学模型(MMA模型), 并将MMA模型作为肌纤维和结缔组织生物力学模型, MMA模型采用完全的应变不变量
${I}_{4},\;{I}_{5},\;{I}_{6}和{I}_{7}$ , 使骨骼肌的剪切行为从材料属性层面得以体现; 结合骨骼肌力学实验结果、RVE模型、肌纤维和结缔组织的生物力学模型, 建立骨骼肌多尺度数值模型. 根据实验结果确定生物力学模型参数, 通过多尺度均匀化方法实现微观尺度和宏观尺度之间的连接, 最终获得骨骼肌宏观力学行为, 通过纵向拉伸、横向拉伸、平面外纵向剪切和平面内剪切4种变形形式的仿真结果验证骨骼肌多尺度数值模型的收敛性. 研究了生物力学模型中模型参数、肌纤维体积分数和肌纤维结构对骨骼肌宏观力学行为的影响, 最后通过实验验证多尺度数值模型的有效性. 本文骨骼肌多尺度数值模型不仅能够用于研究骨骼肌微观因素对宏观力学行为的影响, 也可用于研究疾病对骨骼肌生物力学特性的影响及模拟骨骼肌重塑和再生.Abstract: Aiming at the problems that there is a certain difference between the muscle fiber microstructure model and the image observed under the microscope, the microscopic component biomechanical model cannot effectively capture the mechanical behavior of skeletal muscle during shear deformation, and the high calculation cost of multi-scale numerical models of skeletal muscle. In this thesis, the mechanical properties of skeletal muscle are studied from the perspectives of experiment, multiscale modeling and simulation. Curved-edge Voronoi polygons are proposed as the cross-section of muscle fibers, and the corresponding representative volume element (RVE) is established at the microscale. A new biomechanical model (MMA model) is proposed, and the MMA model is used as the biomechanical model of muscle fibers and connective tissue, the MMA model adopts complete strain invariants$ {I}_{4}、{I}_{5}、{I}_{6}、{I}_{7} $ , so that the shear behavior of skeletal muscle is reflected at the level of material properties. Combine the experimental results of skeletal muscle, the RVE models, the biomechanical models of muscle fibers and connective tissue to establish a multiscale numerical model of skeletal muscle. According to the experimental results, the parameters of the biomechanical model are determined, the multiscale homogenization method are used to realize the connection between the microscale and the macro-scale, and the macroscopic mechanical behavior of skeletal muscle is finally obtained, four deformation forms of Longitudinal stretch, stretch laterally, out-of-plane longitudinal shear and in-plane shear are performed to verify the convergence of the model. This thesis research the effects of model parameters, muscle fiber volume fraction and muscle fiber structure on skeletal muscle on macroscopic mechanical behavior. Combined with experimental data, the effectiveness of the multiscale numerical model is verified. In this paper, the multi-scale numerical model of skeletal muscle can not only be used to study the influence of microscopic factors on the macroscopic mechanical behavior of skeletal muscle, but also to study the influence of diseases on the biomechanical properties of skeletal muscle and to simulate skeletal muscle remodeling and regeneration. -
引 言
人体大约包含647块骨骼肌[1], 占人体重量的
$30\% \sim 40\%$ [2], 在人体内分布广泛, 是人体主要组成器官之一. 根据产生来源的不同, 骨骼肌力行为可以分为主动力行为和被动力行为. 主动力行为来自大脑和脊柱的神经控制, 由大脑和脊柱神经发出的电脉冲刺激肌纤维, 触发肌节缩短, 从而在骨骼肌中产生积极的主动收缩力行为.20世纪60年代, 生物力学作为独立的学科分支[3], 成为了世界范围内研究的热点. 著名美籍华裔学者冯元桢先生致力于生物力学领域的研究[4], 被誉为生物力学之父, 是生物力学的开创者和生物医学工程的奠基人. 自1972年Spencer[5]建立的纤维增强生物力学模型以来, 一系列模型被提出, 并广泛应用于模拟生物软组织、橡胶类材料等. Gerard等[6]在一名74岁女性尸体上取出的舌肌和面部组织, 采用两参数的Mooney-Rivlin模型表征舌肌和面部组织的生物力学行为, 成功地模拟了语言表达期间观察到的舌头运动; 陈伟等[7]通过三维重建技术获得了老年健康女性骨盆与肛提肌群的骨骼-骨骼肌系统, 采用两参数的生物力学Mooney-Rivlin模型, 该模型逼真的反映了盆底肛提肌群的立体空间位置关系; Schiavone等[8]通过两参数Yeoh模型表征了舌肌, 通过吸气实验获得最佳模型参数, 并在仿真中获得满意的仿真精度; Nazari等[9-10]利用5参数Mooney-Rivlin模型解释面部肌肉的非线性行为, 通过有限元仿真得出该模型能逼真模拟语音和非语言交流中的面部运动; Stavness等[11]建立5参数Mooney-Rivlin模型舌肌有限元模型, 基于提出的跟踪算法大大提高了计算效率, 并将该模型应用于模拟舌肌的前伸和弯曲运动, 预测的肌肉活动与舌张力和肌电图的实验记录一致; Gindre等[12]基于Voigt近似方法[13]建立骨骼肌多尺度解析模型, 在纤维方向拉伸时, 肌联蛋白发挥作用, 在压缩时肌联蛋白作用消失, 成功预测了骨骼肌在拉伸和压缩外载荷下的高度不对称性; Ogden模型[14]已经被证明能够模拟生物软组织复杂的力学行为; Al-Dirini等[15]在设计新型座椅坐垫时, 在对肌肉和皮下脂肪组织定义材料属性时才采用一阶Ogden模型; Roux等[16]在骨骼肌-肌腱组件模型中采用Ogden模型表征骨骼肌力学行为; Bosboom等[17]采用不可压缩Ogden模型用于描述骨骼肌被动力学特性; 粟思橙[18]采用Ogden模型模拟骨骼肌超弹性力学行为; Sengeh等[19]采用二阶Ogden模型表征骨骼肌和皮肤组织的非线性弹性行为; Moerman等[20]提出的耦合的幂函数形式生物力学模型用于捕获骨骼肌的非线性和各向异性力学行为; Holzapfel等[21]在2000年提出了一种用于模拟动脉组织的指数形式生物力学模型; Hernández等[22]在Holzapfel等[21]的基础上建立了一种用于模拟骨骼肌的指数形式生物力学模型; Böl等[23]在建立力-电耦合的骨骼肌生物力学模型时, 骨骼肌生物力学模型采用指数函数形式; Wu等[24]采用分段函数重现了人类骨骼肌力-拉伸关系. Calvo等[25]提出改进分段函数形式的生物力学模型, 采用Levenberg-Marquardt优化算法拟合大鼠胫骨前肌和肌腱组织的实验数据.
Spyrou等[26]假设肌纤维的截面为规则的正六边形, 他们根据肌纤维不规则截面设置肌纤维截面为Voronoi多边形(泰森多边形)[27]; Sharafi等[28-29]认为肌纤维束截面是多边形; Jiménez等[30]采用周期性分布的圆形作为肌纤维横截面; Virgilio等[31]采用基于代理的模型作为肌纤维截面形状; Kuravi等[32]通过对猪骨骼肌组织的染色切片进行图像配准, 随后进行图像分割, 确定每个切片中的肌纤维和结缔组织区域, 结合三维重建技术最终建立三维骨骼肌微观结构模型.
Valentin等[33]假设肌纤维为长方体建立一个基于骨骼肌微结构的模型通过数值均匀化方法预测骨骼肌被动力学响应. Kurav等[34]通过数值均匀化方式建立骨骼肌多尺度框架, 相较于他们之前的研究[32], 改进后的模型预测到了骨骼肌的各向异性和拉伸-压缩不对称性等典型力学行为. Spyrou等[27]以泰森多边形为肌纤维截面形状建立骨骼肌单包模型, 并对数值均匀化方法和解析法的仿真结果进行对比, 发现两者仿真结果具有较好的一致性. Bleiler等[35]在Voigt近似方法的基础上建立骨骼肌多尺度解析模型.
针对抽象化的肌纤维的截面形状的几何模型与实验观察所得肌纤维横截面形状有一定的差异、肌纤维和结缔组织的生物力学模型忽略了应变不变量
${I_{\text{5}}}$ 和${I_{\text{7}}}$ 的影响, 使生物力学模型不能有效模拟骨骼肌在剪切变形时的力学特性.本文建立更加接近真实骨骼肌微观结构的RVE模型、建立肌纤维和结缔组织的生物力学模型. 结合RVE模型和骨骼肌组分的生物力学模型实现多尺度数值模型的有限元仿真, 实现对骨骼肌的宏观生物力学行为的预测.
1. 骨骼肌准静态力学实验
为了获得骨骼肌的生物力学特性, 以新鲜猪后腿骨骼肌为实验对象, 对骨骼肌分别进行如图1所示实验, 实验采用Zwickz010电子万能材料试验机完成, 如图2所示. 为了使黏弹性的影响降至最低, 生物组织拉伸实验的应变率范围在
${10^{ - 5}} \sim {10^{ - 1}}\;{{\text{s}}^{ - 1}}$ [36], 所以本实验应变率设置为$5 \times {10^{ - 3}}\;{{\text{s}}^{ - 1}}$ .![]()
图 1 骨骼肌力学实验示意图Figure 1. Schematic diagram of the mechanics experiments of skeletal muscle1.1 试样的拉伸实验
骨骼肌拉伸实验分为纵向拉伸和横向拉伸两组, 每组5个试样, 拉伸速度为3 mm/min, 预载荷为0.1 N, 拉伸直到骨骼肌产生破坏. 使用电子游标卡尺对试样的长、宽和厚分别测量上、中、下3个不同位置, 试样平均尺寸如表1所示.
表 1 拉伸实验试样几何参数Table 1. Geometric parameters of specimens of stretch experimentExperiment No. L/mm W/mm H/mm skeletal muscle longitudinal
stretch experiment1-1 28.01 10.07 4.20 1-2 33.40 8.24 4.48 1-3 33.59 7.77 4.51 1-4 29.99 9.45 4.78 1-5 30.01 10.11 4.53 skeletal muscle lateral stretch experiment 2-1 31.98 11.94 6.26 2-2 36.25 7.27 5.63 2-3 28.71 7.48 5.77 2-4 38.10 13.13 5.21 2-5 33.74 12.03 5.23 利用软件Origin中Savitzky-Golay方法, 对实验数据做平滑预处理, “points of window”的数值设置为50. 当试样被拉伸时, 试样在纵轴方向所受的名义应力
$P$ 为$$ P{\text{ = }}\frac{F}{{{A_0}}} $$ (1) 式中,
$ F $ 为拉伸载荷;$ {A_0} $ 为试样的初始横截面积. 试样的名义应变${{\varepsilon }}$ 表达式为$$ {{\varepsilon }}{\text{ = }}\frac{{\Delta L}}{{{L_0}}} $$ (2) 式中,
$ {L_{\text{0}}} $ 为试样的初始有效长度, 本实验设置上下夹具的距离为10 mm, 即$ {L_{\text{0}}} = 10{\text{mm}} $ ;$\Delta L$ 为被拉伸位移, 实验数据经平滑处理, 结果如图3所示, 仅列举1-1试样纵向拉伸数据和2-1试样横向拉伸数据.骨骼肌拉伸实验名义应力的均值和标准差, 结果如图4所示.
![]()
图 4 拉伸实验名义应力-名义应变曲线Figure 4. Nominal stress-nominal strain curve of stretch experiment1.2 试样的剪切实验
剪切实验包括平面外纵向剪切和平面内剪切, 每种剪切实验包含5个试样, 使用电子游标卡尺对试样的长、宽和厚分别测量上、中、下3个不同位置, 几何参数如表2所示.
表 2 剪切实验试样几何参数Table 2. Geometric parameters of specimens of shear experimentExperiment No. L/mm W/mm H/mm skeletal muscle in-plane shear
experiment3-1 19.30 21.85 5.30 3-2 18.46 21.79 5.69 3-3 18.62 18.60 4.02 3-4 18.07 21.14 4.40 3-5 18.65 21.40 4.50 skeletal muscle out-of-plane
longitudinal shear experiment4-1 22.84 16.50 4.11 4-2 20.75 20.16 5.18 4-3 26.37 15.81 4.65 4-4 22.87 23.25 3.31 4-5 17.98 17.74 2.70 剪切实验数据采用Adjacent-Averaging平滑方法. 在平滑处理时, 为了尽可能增加平滑效果, “points of window”的数值设置为100, 平滑后结果如图5中红色曲线所示.
剪切应变
${\boldsymbol{\gamma }}$ 的计算公式为$$ {\boldsymbol{\gamma }}{\text{ = tan}}\frac{{\Delta L}}{n} $$ (3) 式中,
$ \Delta L $ 为剪切位移, 单位为mm;$n$ 为剪切区域宽度,$n$ = 6.5 mm.名义剪切应力
$ {\boldsymbol{\tau }} $ 的计算公式为$$ {\boldsymbol{\tau }} = \frac{F}{{L h}} $$ (4) 式中,
$ F $ 为实验中施加的力, 单位为N;$ L $ 和$h$ 分别为剪切变形区域的长度和厚度.对最终数据进行了数据拟合, 得到平滑的剪切应力-剪切应变曲线, 结果如图6和图7所示.
![]()
图 6 平面内剪切的名义剪切应力−剪切应变曲线Figure 6. Nominal shear stress-shear strain curve of in-plane shear![]()
图 7 平面外纵向剪切的名义剪切应力−剪切应变曲线Figure 7. Nominal shear stress-shear strain curve of out-of-plane longitudinal shear2. 骨骼肌微观尺度RVE模型的建立
2.1 骨骼肌RVE模型结构的确定
骨骼肌是微观结构复杂的多相介质材料, 在几何结构上, 骨骼肌宏观几何模型被认为是由有限数量的微观RVE模型周期排列而成, 如图8所示.
研究表明, 肌纤维半径在10 ~ 100 μm (假设肌纤维为圆柱形的情况下)[38], 骨骼肌的宏观尺度为1 mm或1 cm, 因此本文选取RVE的尺寸为350 μm × 350 μm × 350 μm, 通过单根肌纤维横截面面积
$ {S_f} $ 和RVE模型横截面面积$ {S_{{\rm{RVE}}}} $ 估算每个RVE模型中包含肌纤维数量${n_f}$ $$ {n_f} = \frac{{{S_{{\rm{RVE}}}}}}{{{S_f}}} $$ (5) 近似认为肌纤维半径为
${{20}}$ μm, 则每个RVE模型中包含约为49根肌纤维, 这里取整, 认为RVE模型中包含50根肌纤维. 因此在生成RVE模型时, 设置肌纤维的种子数量为50.泰森多边形又被称为冯洛诺伊图(Voronoi diagram), 是由相邻种子点的中垂线组成的多边形, 如图9所示. 泰森多边形经常用于模拟多晶体微结构[39], 2011年, Sharafi等[29]首次将泰森多边形引入到肌纤维微观结构建模领域.
2.2 泰森多边形RVE模型的建立
ABAQUS借助插件Homtools生成具有泰森多边形的二维RVE模型, 按照2.1节中确定的RVE几何模型的尺寸和种子数量, 参数设置如图10所示, 得到的二维泰森多边形及其种子分布如图11所示.
多尺度计算方法要求RVE几何模型满足对称结构, 为了满足对称结构, 在二维泰森多边形外包一层结缔组织, 实现上下边界和左右边界完全是结缔组织, 结缔组织层厚根据肌纤维体积分数而设置, 处理后的结果如图12所示.
由二维泰森多边形RVE模型生成三维泰森多边形RVE模型, 对图12所示的二维模型划分网格, 如图13(a)所示(红线所占网格为结缔组织), 创建孤立网格部件, 对孤立网格部件的网格进行网格偏置, 设置要生成三维网格模型的厚度, 设置层数, 最终即可生成三维网格模型, 如图13(b)所示.
2.3 曲边泰森多边形RVE模型的建立
曲边泰森多边形是在泰森多边形的基础上进行曲边优化处理得到的. 曲边优化处理的原理为: 标记每个二维RVE模型中泰森多边形各边的中点, 用B样条曲线连接泰森多边形各边中点, 最终由B样条曲线构成曲边优化后的曲边泰森多边形的边. 曲边优化过程如图14所示, 获得曲边泰森多边形的二维RVE模型如图15所示. 曲边优化后的二维RVE模型生成三维RVE模型, 采用2.2节中介绍的三维泰森多边形RVE模型的生成方式, 最终获得三维RVE模型的网格模型如图16所示.
生成肌纤维和结缔组织的三维网格模型如图17和图18所示, 生成的肌纤维和结缔组织为共节点连接, 满足了肌纤维和结缔组织完全绑定连接要求.
2.4 不同体积分数RVE模型的建立
为了研究肌纤维体积分数和肌纤维截面形状对骨骼肌力学特性的影响, 分别建立了肌纤维截面为泰森多边形和曲边泰森多边形, 体积分数分别为0.7, 0.8, 0.9的RVE模型, 建立的RVE模型如表3所示, 根据肌纤维体积分数, 经过多次尝试得对应结缔组织层厚(ep)数值如表4所示. 肌纤维的体积分数是通过二维模型中肌纤维与RVE面积的比值获得, 计算公式为
表 4 ep设置值Table 4. Setting value of epVolume fraction Voronoi polygons/μm Curved-edges Voronoi
polygons/μm0.7 0.0080 0.005229 0.8 0.0052 0.003000 0.9 0.0025 0.000540 $$ {v_f} = \frac{{{S_{{\rm{RVE}} - f}}}}{{{S_{{\rm{RVE}}}}}} $$ (6) 式中
${S_{{\rm{RVE}} - f}}$ 是RVE二维模型中肌纤维的面积;${S_{{\rm{RVE}}}}$ 是RVE二维模型的面积. 其中${S_{{\rm{RVE}} - f}}$ 由ABAQUS中查询.表 3 森多边形和曲边泰森多边形RVE模型Table 3. The RVE model of Voronoi polygons and curved-edge Voronoi polygonsVolume fraction Voronoi polygons RVE model Curved-edges Voronoi polygons RVE model 0.7 
0.8 
0.9 
2.5 周期性边界条件的添加
基于EasyPBC的灵活性, 本节采用EasyPBC对RVE网格模型添加周期性边界条件[40-43], 其流程如图19所示.
3. 基本模型
3.1 模型参数
由于肌纤维为横观各向同性超弹性材料, 用应变能方程表示材料的力学特性, 肌纤维应变能方程
$W$ 的一般形式为$$ W = W\left( {{I_1},{I_2},{I_3},{I_4},{I_5}} \right) $$ (7) 结缔组织为各向异性材料, 其应变能方程
$W$ 的一般形式为$$ W = W\left( {{I_1},{I_2},{I_3},{I_4},{I_5},{I_6},{I_7}} \right) $$ (8) 式中,
${I}_{1},\;{I}_{2},\;{I}_{3},\;{I}_{4},\;{I}_{5},\;{I}_{6},\;{I}_{7}$ 为应变不变量, 计算公式为$$ \left.\begin{split} & {I_1} = {\text{tr}}({\boldsymbol{C}}),{I_2} = \frac{1}{2}\left[ {{{\rm{tr}}} {{({\boldsymbol{C}})}^2} - {{\rm{tr}}} ({{\boldsymbol{C}}^2})} \right],{I_3} = \det ({\boldsymbol{C}}) \\ & J = \det ({\boldsymbol{F}}) = \sqrt {{I_3}} ,{I_4} = {{\boldsymbol{M}}_{\text{1}}} \cdot {\boldsymbol{C}}{{\boldsymbol{M}}_{\text{1}}},{I_5} = {{\boldsymbol{M}}_{\text{1}}} \cdot {{\boldsymbol{C}}^2}{{\boldsymbol{M}}_{\text{1}}} \\ & {I_6} = {{\boldsymbol{M}}_{\text{2}}} \cdot {\boldsymbol{C}}{{\boldsymbol{M}}_{\text{2}}},{I_7} = {{\boldsymbol{M}}_{\text{2}}} \cdot {{\boldsymbol{C}}^2}{{\boldsymbol{M}}_{\text{2}}} \end{split}\right\} $$ (9) 式中,
$J$ 为体积比, 对于不可压缩材料$J = 1$ . 右柯西格林变形张量${\boldsymbol{C = }}{{\boldsymbol{F}}^{\rm{T}}}{\boldsymbol{F}}$ , 左柯西格林变形张量${\boldsymbol{B}}= {\boldsymbol{F}} \cdot {{\boldsymbol{F}}^{\rm{T}}}$ [44].$ {\boldsymbol{M}} $ 表示肌纤维的最优方向,$ {\boldsymbol{M}} $ 为单位矢量. 在笛卡尔坐标中, 设定$ {\boldsymbol{M}} $ 方向平行于Y轴, 如图20所示,$ {\boldsymbol{M}} $ 的向量表示为![]()
图 20 骨骼肌肌纤维和胶原纤维的分布Figure 20. Distribution of skeletal muscle fibers and collagen fibers$$ {\boldsymbol{M}} = \left( {\begin{array}{*{20}{c}} 0&1&0 \end{array}} \right) $$ (10) 分别用
${{\boldsymbol{M}}_1}$ 和${{\boldsymbol{M}}_2}$ 表示纤维的两个方向, 假设I型胶原纤维分布在笛卡尔坐标内的X-Y平面, 胶原纤维与X轴夹角为${\theta ^m}$ , 如图20所示, 则${{\boldsymbol{M}}_1}$ 和${{\boldsymbol{M}}_2}$ 可表示为$$ \left. \begin{gathered} {{\boldsymbol{M}}_{\text{1}}} = \left( {\begin{array}{*{20}{c}} {\cos {\theta ^m}}&{\sin {\theta ^m}}&{{\text{ }}0} \end{array}} \right) \\ {{\boldsymbol{M}}_2} = \left( {\begin{array}{*{20}{c}} {\cos {\theta ^m}}&{ - \sin {\theta ^m}}&0 \end{array}} \right) \\ \end{gathered} \right\} $$ (11) 为了表征骨骼肌的可压缩性, 应变能方程被解耦为体积-等体积的和, 这种解耦基于Flory[45]提出的变形梯度的乘法分解, 即
$$ {\boldsymbol{F}} = {J^{1/3}}{{\boldsymbol{F}}^*} $$ (12) 式中,
${{\boldsymbol{F}}^*}$ 是${\boldsymbol{F}}$ 的等体积部分, 相应的应变不变量的${I}_{1},\;{I}_{2},\;{I}_{3},\; {I}_{4},\;{I}_{5},\;{I}_{6},\;{I}_{7}$ 等体积表示为${I}_{1}{}^{*},\;{I}_{2}{}^{*},\;{I}_{3}{}^{*},\; {I}_{4}{}^{*},\; {I}_{5}{}^{*},\;{I}_{6}{}^{*},\;{I}_{7}{}^{*}$ , 被称为等体积应变不变量或伪应变不变量, 两种应变不变量的关系如下$$ \left. \begin{gathered} I_a^* = {J^{ - 2/3}}{I_{(a)}},a \in \{ 1,4,6,8\} \\ I_b^* = {J^{ - 4/3}}{I_b},b \in \{ 2,5,7\} \\ \end{gathered} \right\} $$ (13) 根据应变能方程的解耦, 应变能方程式(7)和式(8)被表示为
$$ W = f\left( J \right) + {W_{^{{\rm{iso}}}}} + {W_{^{{\rm{aniso}}}}} $$ (14) 式中,
$f(J)$ 为应变能方程的体积部分, 表征模型的体积变化;${W_{{\rm{iso}}}}$ 为各向同性部分, 表征模型基质的力学特性;${W_{{\rm{aniso}}}}$ 为各向异性部分, 表征纤维的增强作用.由能量守恒和角动量守恒的原理, 可以得到应变能对时间的导数的表达式
$$ \dot W = {\boldsymbol{S}} \cdot {{\dot {\boldsymbol{E}}}} = \frac{1}{2}{\boldsymbol{S}} \cdot {{\dot {\boldsymbol{C}}}} $$ (15) 式中S为第二类皮奥拉-基尔霍夫(Piola-Kirchhoff, P-K)应力张量, E是格林应变张量, E与C的关系为
$$ {\boldsymbol{E}} = \frac{1}{2}({\boldsymbol{C}} - {\boldsymbol{I}}) $$ (16) 通过链式规则, 式(15)转化为
$$ \left( {{\boldsymbol{S}} - 2\frac{{\partial W}}{{\partial {\boldsymbol{C}}}}} \right) \cdot \frac{1}{2}\mathop {\boldsymbol{C}}\limits^. = 0 $$ (17) 进而可得
$$ {\boldsymbol{S}} = 2\frac{{\partial W}}{{\partial {\boldsymbol{C}}}} $$ (18) 3.2 MMA模型
本文在MA模型[46-47]的基础上提出MMA模型, 采用的方式为: 在MA模型的各向异性部分添加应变不变量
$ {I_5} $ 和$ {I_7} $ 的影响, 最终建立的MMA模型中各向异性部分的表达式为$$ \begin{split} & {W_{{\rm{aniso}}}} = \frac{{{c_2}}}{{2k}}\sum\limits_{i = 4,6} {\left\{ {\exp \left[ {k{{\left( {{I_i} - 1} \right)}^2}} \right] - 1} \right\}} + \\ &\qquad {c_3}\left( {{I_5} + {I_7} - 2{I_4} - 2{I_6} + 2} \right) \end{split} $$ (19) 所以MMA模型模拟可压缩性、超弹性和各向异性的表达式为
$$ \begin{split} & W = f(J) + {W_{{\rm{iso}}}} + {W_{{\rm{aniso}}}} = \;\frac{{{\kappa _0}}}{2}{(J - 1)^2} + \\ &\qquad \frac{{{c_1}}}{2}\left( {I_1^* - 3} \right) + {\text{ }}\frac{{{c_2}}}{{2k}}\sum\limits_{i = 4,6} {\left\{ {\exp \left[ {k{{\left( {{I_i} - 1} \right)}^2}} \right] - 1} \right\}} + \\ &\qquad {c_3}\left( {{I_5} + {I_7} - 2{I_4} - 2{I_6} + 2} \right) \end{split} $$ (20) 根据MA模型柯西应力所得过程, 求得MMA模型的柯西应力的表达式为
$$ \begin{split} & \sigma {\text{ = }}2{W_J}{\boldsymbol{I}} + 2{J^{ - 5/3}}{W_1}\left( {{\boldsymbol{B}} - \frac{1}{3}{I_1}{\boldsymbol{I}}} \right)\; + \\ & \qquad\frac{2}{J}{W_4}{\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1} + \frac{2}{J}{W_6}{\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}\; + \\ &\qquad \frac{2}{J}{W_5}({\boldsymbol{BF}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1} + {\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{BF}}{{\boldsymbol{M}}_1})\; + \\ & \qquad\frac{2}{J}{W_7}({\boldsymbol{BF}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2} + {\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{BF}}{{\boldsymbol{M}}_2}) \end{split} $$ (21) 同样采用MA模型中的数学推导过程, 对MMA模型的柯西应力求矩阵的迹, 得
$$ {\text{tr}}({\boldsymbol{\sigma }}) = 3{f^\prime }(J) + \frac{2}{J}\left( {{W_4} + {W_6}} \right)\left( {{C_{11}}{{\cos^2 \theta }} + {C_{22}}{{\sin^2 \theta }}} \right) $$ (22) 根据式(21), MMA模型在静水压力作用下的应力与体积变化和纤维增强作用有关, 因此MMA模型能够模拟可压缩性和各向异性.
对于剪切变形, 同样以X-Y平面的剪切变形为例, 根据MA模型的变形梯度F, 计算获得柯西应力分量
${\sigma _{12}}$ 的表达式为$$ \begin{split} & {\sigma _{12}}{\text{ = }}2{J^{ - 5/3}}{W_1}{{\text{B}}_{12}} +\\ &\qquad\frac{2}{J}{W_5}{({\boldsymbol{BF}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1} + {\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{BF}}{{\boldsymbol{M}}_1})_{12}} + \\ & \qquad \frac{2}{J}{W_7}{({\boldsymbol{BF}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2} + {\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{BF}}{{\boldsymbol{M}}_2})_{12}}= \\ & \qquad 2{J^{ - 5/3}}{W_1}\gamma + \frac{2}{J}\left( {{W_5} + {W_7}} \right)\cdot\\ &\qquad\left[ {\gamma {{\left( {\sin \theta + \gamma \cos \theta } \right)}^2} + \gamma {{ {\cos^2 \theta }}}} \right] \end{split} $$ (23) 对比MA模型和MMA模型的柯西应力分量
${\sigma _{12}}$ , MA模型的柯西应力分量${\sigma _{12}}$ 仅包含剪切应变$\gamma $ 的影响, 并没有纤维相关量, 因此在剪切变形时表现为各向同性. 而MMA模型柯西应力分量${\sigma _{12}}$ , 包含纤维增强的相关量${W_5}$ ,${W_7}$ 和$\theta $ . 可知, MMA模型在模拟剪切变形时, 柯西应力分量${\sigma _{12}}$ 受到纤维增强作用的影响. 因此剪切变形时, 剪切平面X-Y平面受到纤维增强作用, 使同X-Y平面的剪切模量包含纤维的增强作用, 从而在剪切变形时体现各向异性.3.3 MMA模型条件验证
(1) 初始条件验证
骨骼肌在初始状态(参考构型)下, 假设骨骼肌未发生变形, 在不考虑预紧力和重力的作用下, 骨骼肌应力为0, 相应的应变不变量的初始值分别为
$$ I_1^* = 3,\;\;{J^2} = {I_3} = {I_4} = {I_5} = {I_6} = {I_7} = 1 $$ (24) 在初始条件下, MMA需要满足以下要求
$$\qquad\qquad {W^0} = 0 $$ (25) $$\qquad\qquad W_1^0 + 2W_2^0 + W_3^0 = 0 $$ (26) $$\qquad\qquad W_4^0 + 2W_5^0 = 0 $$ (27) $$\qquad\qquad W_6^0 + 2W_7^0 = 0 $$ (28) 式中上标0表示参考构型状态下的值.
根据MMA模型分别计算式(25) ~ 式(28)可得
$$ \begin{split} & {W^0} = \frac{{{\kappa _0}}}{2}{(1 - 1)^2} + \frac{{{c_1}}}{2}\left( {3 - 3} \right) + \\ & \qquad\frac{{{c_2}}}{{2k}}\sum\limits_{i = 4,6} {\left\{ {\exp \left[ {k{{\left( {1 - 1} \right)}^2}} \right] - 1} \right\}} + \\ & \qquad{c_3}\left( {1 + 1 - 2 - 2 + 2} \right) = 0 \end{split} $$ (29) $$ W_1^0 = 2W_2^0 = W_3^0 = 0 $$ (30) $$ \begin{split} & W_4^0 + 2W_5^0 = W_6^0 + 2W_7^0= \\ &\qquad \frac{{{c_2}}}{{2k}}\left[ {k{{\left( {1 - 1} \right)}^2}} \right] \times 2 \times (1 - 1)\exp \left[ {k{{\left( {1 - 1} \right)}^2}} \right] - \\ &\qquad2{c_3} + 2{c_3} = 0 \end{split} $$ (31) 根据式(29) ~ 式(31)的计算结果可得, MMA模型满足式(25) ~ 式(28)初始条件要求.
(2) 连续性条件验证
由于MMA模型是
$ f(J) $ ,${W_{{\rm{iso}}}}$ ,${W_{{\rm{aniso}}}}$ 的叠加, 当$ f(J) $ ,${W_{{\rm{iso}}}}$ ,${W_{{\rm{aniso}}}}$ 分别满足连续性要求, 则MMA模型满足连续性要求.对于
$ f(J) $ $$ \frac{{\partial f(J)}}{{\partial {\boldsymbol{C}}}} = \frac{{\partial f(J)}}{{\partial J}} \times \frac{{\partial J}}{{\partial {\boldsymbol{C}}}} = \frac{1}{2}J{{\boldsymbol{C}}^{ - 1}}\left( {J - 1} \right) $$ (32) 对于
$ {W_{{\rm{iso}}}} $ $$ \frac{{\partial {W_{{\rm{iso}}}}}}{{\partial {\boldsymbol{C}}}} = \frac{{\partial {W_{{\rm{iso}}}}}}{{\partial I_1^*}} \times \frac{{\partial I_1^*}}{{\partial {\boldsymbol{C}}}} = \frac{{{c_1}}}{2}\left( {I_3^{ - 1/3}I - \frac{1}{3}I_1^ * {C^{ - 1}}} \right) $$ (33) 对于
$ {W_{{\rm{aniso}}}} $ $$ \begin{split} & \frac{{\partial {W_{{\rm{aniso}}}}}}{{\partial {\boldsymbol{C}}}} = \frac{{\partial {W_{{\rm{aniso}}}}}}{{\partial {I_4}}} \times \frac{{\partial {I_4}}}{{\partial {\boldsymbol{C}}}} + \frac{{\partial {W_{{\rm{aniso}}}}}}{{\partial {I_5}}} \times \frac{{\partial {I_5}}}{{\partial {\boldsymbol{C}}}} + \frac{{\partial {W_{{\rm{aniso}}}}}}{{\partial {I_6}}} \times \\ &\qquad\frac{{\partial {I_6}}}{{\partial {\boldsymbol{C}}}} + \frac{{\partial {W_{{\rm{aniso}}}}}}{{\partial {I_7}}} \times \frac{{\partial {I_7}}}{{\partial {\boldsymbol{C}}}} = {W_4}{\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1} + \\ &\qquad {W_5}\left( {{\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{CF}}{{\boldsymbol{M}}_1} + {\boldsymbol{CF}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}} \right) + \\ & \qquad{W_6}{\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2} + {W_7}\left( {{\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{CF}}{{\boldsymbol{M}}_2} + {\boldsymbol{CF}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}} \right) \end{split} $$ (34) 通过计算可得
$ f(J) $ ,${W_{{\rm{iso}}}}$ 和${W_{{\rm{aniso}}}}$ 的导数都是存在的, 因此MMA模型是连续的, 可得MMA模型在求解域内是连续的.(3) 稳定性条件验证
稳定性要求保证了应变能函数局部极小值的存在. 为了满足稳定性标准, 应变能方程必须是凸函数[48], 如果3部分都是凸函数, 则应变能方程满足凸函数要求[49], 即体积部分、各向同性部分和各向异性部分的二阶偏导数大于等于0.
$$ \;\;f{(J)_{(JJ)}} = \frac{{{\partial ^2}f(J)}}{{\partial {J^2}}} = {\kappa _0} $$ (35) $$ {W_{11}} = \frac{{{\partial ^2}W}}{{\partial {{ {{I^2_1}}}}}} = 0 $$ (36) $$ \begin{split} &{W_{44}} = \frac{{{\partial ^2}W}}{{\partial {{ {{I^2_4}}}}}} = [ 3{c_{\text{2}}}{\left( {{I_4} - 1} \right)^2} + 2{c_{\text{2}}}k{\left( {{I_4} - 1} \right)^4}{{] }}\cdot\\ &\qquad \exp \left[ {k{{\left( {{I_4} - 1} \right)}^2}} \right]\end{split} $$ (37) $$ \begin{split} &{W_{66}} = \frac{{{\partial ^2}W}}{{\partial {{{{I^2_6}} }}}} = [ 3{c_{\text{2}}}{\left( {{I_6} - 1} \right)^2} + 2{c_{\text{2}}}k{\left( {{I_6} - 1} \right)^4}{{] }}\cdot\\ &\qquad\exp \left[ {k{{\left( {{I_6} - 1} \right)}^2}} \right]\end{split} $$ (38) $$ {W_{55}} = \frac{{{\partial ^2}W}}{{\partial {{ {{I^2_5}} }}}} = 0 $$ (39) $$ {W_{77}} = \frac{{{\partial ^2}W}}{{\partial {{ {{I^2_7}} }}}} = 0 $$ (40) 显然, 当
$ {\kappa _0} $ ,$ {c_{\text{2}}} $ 和$ k $ 均大于等于0时, MMA模型满足稳定性要求. 由于骨骼肌各向异性和可压缩性要求,$ {\kappa _0} $ ,$ {c_{\text{2}}} $ 和$ k $ 不能为0, 因此要求参数$ {\kappa _0} $ ,$ {c_{\text{2}}} $ 和$ k $ 必须大于0.经过分析, MMA模型可以表征可压缩性、各向异性、超弹性和剪切变形力学行为, 并且满足初始条件要求、连续性要求和稳定性要求. MMA模型适用于各向异性材料, 能直接应用于结缔组织的模拟, 但是对于肌纤维(横观各向同性材料), MMA模型的各向异性部分应该改写为
$$ {W_{{\rm{aniso}}}} = \frac{{{c_2}}}{{2k}}\exp \left[ {k{{\left( {{I_{\text{4}}} - 1} \right)}^2}} \right] - 1 + {c_3}\left( {{I_5} - 2{I_4} + {\text{1}}} \right) $$ (41) 经过验证, 式(41)表征的横观各向同性材料模型同样需要满足初始条件要求、稳定性要求和连续性要求.
3.4 用户材料子程序的开发
本文涉及的仿真运算属于静力学范畴, 因此选用ABAQUS/Standard(隐式求解器), 开发对象为UMAT. 在Visual Studio2015中创建.for工程, UMAT子程序的一般形式为
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,
DSTRAN, TIME,
2 DTIME,TEMP,DTEMP,PREDEF,DPRED,
MATERL,NDI,NSHR, NTENS,
3 NSTATV,PROPS,NPROPS,COORDS,DROT,
PNEWDT,CELENT,
4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,
KSPT,KSTEP,KINC)
INCLUDE 'ABA_PARAM.INC
$\vdots $ C求解柯西应力STREE
C求解雅可比矩阵DDSDDE
C 柯西应力STREE和雅可比矩阵DDSDDE的转化
$\vdots $ RETURN
END SUBROUTINE UMAT
通过FORTRAN语言编辑完成的用户材料子程序UMAT在每个积分点的运算中都会被ABAQUS/Standard中被调用. 在有限元计算过程中, UMAT被调用过程如图21所示.
在3.2节中已经计算获得了MMA模型的柯西应力的最终表达式, 但是UMAT中还需要计算雅可比矩阵DDSDDE, 用于计算应力更新.
根据式(15)、式(18)和式(19)可得MMA模型的第二个皮奥拉-基尔霍夫应力张量的表达式为
$$ \begin{split} & {\boldsymbol{S}}/2 = W_1^* \times \frac{{\partial I_{\text{1}}^{\text{*}}}}{{\delta {\boldsymbol{C}}}} + f{\left( J \right)_J} \times \frac{1}{2}J{{\boldsymbol{C}}^{ - 1}} + \\ & \qquad{W_4}\frac{{\partial {I_4}}}{{\delta {\boldsymbol{C}}}} + {W_6}\frac{{\partial {I_6}}}{{\delta {\boldsymbol{C}}}} + {W_5}\frac{{\partial {I_5}}}{{\delta {\boldsymbol{C}}}} + {W_7}\frac{{\partial {I_7}}}{{\delta {\boldsymbol{C}}}} \end{split} $$ (42) 对式中应变不变量偏导数求解, 分别设
${{\boldsymbol{N}}}_{1} = {{\boldsymbol{M}}}_{1}\otimes {\boldsymbol{C}}{{\boldsymbol{M}}}_{1} + {\boldsymbol{C}}{{\boldsymbol{M}}}_{1}\otimes {{\boldsymbol{M}}}_{1},\;{{\boldsymbol{N}}}_{2} = {{\boldsymbol{M}}}_{2}\otimes {\boldsymbol{C}}{{\boldsymbol{M}}}_{2} + {\boldsymbol{C}}{{\boldsymbol{M}}}_{2}\otimes {{\boldsymbol{M}}}_{2}$ 得$$ \begin{split} & {\boldsymbol{S}}/2 = W_1^* \times \left( {I_3^{ - 1/3}{\boldsymbol{I}} - \frac{1}{3}I_{\text{1}}^{\text{*}}{{\boldsymbol{C}}^{ - 1}}} \right) + f{\left( J \right)_J} \times \frac{1}{2}J{{\boldsymbol{C}}^{ - 1}} + \\ & \qquad{W_4} \times {{\boldsymbol{M}}_1} \otimes {{\boldsymbol{M}}_1} + {W_6} \times {{\boldsymbol{M}}_2} \otimes {{\boldsymbol{M}}_2} + {W_5} \times \\ &\qquad{{{{\boldsymbol{N}}}}_1} + {W_7} \times {{\boldsymbol{N}}_2}\\[-12pt] \end{split} $$ (43) 根据式(43)定义的非线性本构方程可转化为如下增量形式
$$ \Delta {\boldsymbol{S}} = {\boldsymbol{\varPi}} \cdot \frac{1}{2}\Delta {\boldsymbol{C}} $$ (44) $\Delta {\boldsymbol{S}}$ 与$\Delta {\boldsymbol{C}}$ 之间是线性关系, 所以式(44)通常称为线性化本构方程,${\boldsymbol{\varPi}}$ 为四阶弹性张量, 这个数值近似是基于基尔霍夫应力的Jaumann率的线性化增量形式[50], 其被定义为$${\boldsymbol{\varPi}} = 2\frac{{\partial {\boldsymbol{S}}}}{{\partial {\boldsymbol{C}}}} = 4\frac{{{\partial ^2}W}}{{\partial {{\boldsymbol{C}}^2}}} $$ (45) 可得MMA模型的4阶弹性张量
${\boldsymbol{\varPi}}$ 表示为$$ \begin{split} & \frac{1}{4}{\boldsymbol{\varPi}} = J f{\left( J \right)_J} \left( {{{\boldsymbol{C}}^{ - 1}} \otimes {{\boldsymbol{C}}^{ - 1}} - 2{{\boldsymbol{I}}_{{{{C}}^{ - 1}}}}} \right) + \\ &\qquad{J^2} f{\left( J \right)_{JJ}} \left( {{{\boldsymbol{C}}^{ - 1}} \otimes {{\boldsymbol{C}}^{ - 1}}} \right) - \frac{4}{3}{J^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. } 4}}}\cdot\\ &\qquad W_1^*\left[ {{\boldsymbol{I}} \otimes {{\boldsymbol{C}}^{ - 1}} + {{\boldsymbol{C}}^{ - 1}} \otimes {\boldsymbol{I}} - {I_1}\left( {{{\boldsymbol{I}}_{{{\boldsymbol{C}}^{ - 1}}}} + \frac{1}{3}{{\boldsymbol{C}}^{ - 1}} \otimes {{\boldsymbol{C}}^{ - 1}}} \right)} \right]+ \\ &\qquad \;{W_{44}}{{\boldsymbol{M}}_1} \otimes {{\boldsymbol{M}}_1} \otimes {{\boldsymbol{M}}_1} \otimes {{\boldsymbol{M}}_1} + {W_{66}}{{\boldsymbol{M}}_2} \otimes {{\boldsymbol{M}}_2} \otimes {{\boldsymbol{M}}_2} \otimes {{\boldsymbol{M}}_2}+ \\ &\qquad {W_5}\frac{{\partial {{\boldsymbol{N}}_1}}}{{\partial {\boldsymbol{C}}}} + {W_7}\frac{{\partial {{\boldsymbol{N}}_2}}}{{\partial {\boldsymbol{C}}}}\\[-12pt] \end{split} $$ (46) 将参考构型中所建立的增量本构方程(46)可转化为现时构型
$$ {\mathcal{L}_v}{\boldsymbol{\tau }}= \Delta {\boldsymbol{\tau}} - \left( {\Delta {\boldsymbol{F}}{{\boldsymbol{F}}^{ - 1}}} \right){\boldsymbol{\tau}} -{\boldsymbol{ \tau}} {\left( {\Delta {\boldsymbol{F}}{{\boldsymbol{F}}^{ - 1}}} \right)^{\rm{T}}} = {{\boldsymbol{\varPi}}^{\tau c}} \cdot \Delta {\boldsymbol{D}} $$ (47) 式中
${\mathcal{L}_v}{\boldsymbol{\tau }}$ 为基尔霍夫应力${\boldsymbol{\tau }}$ 的奥尔德罗伊德率增量形式;${{\boldsymbol{\varPi}}^{\tau c}}$ 是一种变换后的弹性张量形式, 其分量为$$ {\boldsymbol{\varPi}}_{ijkl}^{\tau c} = \frac{1}{J}{F_{iI}}{F_{jJ}}{F_{kK}}{F_{IL}}{\mathbb{C}_{ILKL}} $$ (48) 将式(46)代入式(48), 可得MMA模型推导出4阶弹性张量
${{\boldsymbol{\varPi}}^{\tau c}}$ 的分量表达式为$$ \begin{split} & \frac{J}{4}{\boldsymbol{\varPi}}_{ijkl}^{\tau c} = J \times f{\left( J \right)_J} \times \left[ {{{\left( {{\boldsymbol{I}} \otimes {\boldsymbol{I}}} \right)}_{ijkl}} - {\delta _{ik}}{\delta _{jl}} - {\delta _{il}}{\delta _{jk}}} \right] + \\ & \qquad {J^2} \times f{\left( J \right)_{JJ}} \times {\left( {{\boldsymbol{I}} \otimes {\boldsymbol{I}}} \right)_{ijkl}} - \frac{4}{3}{J^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. } 4}}}W_1^*\Biggr\{ {{{\left( {{\boldsymbol{B}} \otimes {\boldsymbol{I}} + {\boldsymbol{I}} \otimes {\boldsymbol{B}}} \right)}_{ijkl}}} - \\ &\qquad {{I_1}\left[ { - \frac{1}{2}\left( {{\delta _{ik}}{\delta _{jl}} + {\delta _{il}}{\delta _{jk}}} \right) + \frac{1}{3}{{\left( {{\boldsymbol{I}} \otimes {\boldsymbol{I}}} \right)}_{ijkl}}} \right]} \Biggr\}\; + \\ & \qquad {W_{44}}{\left( {{\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}} \right)_{ijkl}} + \\ & \qquad {W_{66}}{\left( {{\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}} \right)_{_{ijkl}}} + \\ & \qquad {W_5}\frac{1}{2}\left[ {{\delta _{ik}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}{\text{)}}}_{jl}}{\text{ + }}{\delta _{il}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}{\text{)}}}_{jk}}} \right. + \\ & \qquad \left. {{\delta _{jk}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}{\text{)}}}_{il}}{\text{ + }}{\delta _{jl}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}{\text{)}}}_{ik}}} \right] + \\ & \qquad {W_7}\frac{1}{2}\left[ {{\delta _{ik}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}{\text{)}}}_{jl}}{\text{ + }}{\delta _{il}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}{\text{)}}}_{jk}}} \right. + \\ & \qquad \left. {{\delta _{jk}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}{\text{)}}}_{il}}{\text{ + }}{\delta _{jl}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}{\text{)}}}_{ik}}} \right]\\[-12pt] \end{split} $$ (49) ABAQUS利用了zaremba-jaumann率的增量形式, 结合式(49)和
${\boldsymbol{\varPi}}_{ijkl}^{\tau Z - J}$ [51], 可得MMA模型的雅可比矩阵的最终表达式(50), 将该式在UMAT通过FORTRAN语言实现代码编程$$ \begin{split} & \frac{J}{4}{\boldsymbol{\varPi}}_{ijkl}^{\tau Z - J} = J \times f{\left( J \right)_J} \times \left[ {{{\left( {{\boldsymbol{I}} \otimes {\boldsymbol{I}}} \right)}_{ijkl}} - {\delta _{ik}}{\delta _{jl}} - {\delta _{il}}{\delta _{jk}}} \right] + \\ & \qquad{J^2} \times f{\left( J \right)_{JJ}} \times {\left( {{\boldsymbol{I}} \otimes {\boldsymbol{I}}} \right)_{ijkl}} - \frac{4}{3}{J^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. } 4}}}W_1^*\Biggr\{ {{{\left( {{\boldsymbol{B}} \otimes {\boldsymbol{I}} + {\boldsymbol{I}} \otimes {\boldsymbol{B}}} \right)}_{ijkl}}} - \\ & \qquad {{I_1}\left[ { - \frac{1}{2}\left( {{\delta _{ik}}{\delta _{jl}} + {\delta _{il}}{\delta _{jk}}} \right) + \frac{1}{3}{{\left( {{\boldsymbol{I}} \otimes {\boldsymbol{I}}} \right)}_{ijkl}}} \right]}\Biggr\}\; + \\ &\qquad \;{W_{44}}{\left( {{\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}} \right)_{ijkl}} + \\ & \qquad{W_{66}}{\left( {{\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}} \right)_{_{ijkl}}} + \\ & \qquad{W_5}\frac{1}{2}\left[ {{\delta _{ik}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}{\text{)}}}_{jl}}{\text{ + }}{\delta _{il}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}{\text{)}}}_{jk}}} \right.{\text{ + }} \\ & \qquad\left. {{\delta _{jk}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}{\text{)}}}_{il}}{\text{ + }}{\delta _{jl}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_1} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_1}{\text{)}}}_{ik}}} \right] + \\ & \qquad{W_7}\frac{1}{2}\left[ {{\delta _{ik}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}{\text{)}}}_{jl}}{\text{ + }}{\delta _{il}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}{\text{)}}}_{jk}}} \right.{\text{ + }} \\ & \qquad\left. {{\delta _{jk}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}{\text{)}}}_{il}}{\text{ + }}{\delta _{jl}}{{({\boldsymbol{F}}{{\boldsymbol{M}}_2} \otimes {\boldsymbol{F}}{{\boldsymbol{M}}_2}{\text{)}}}_{ik}}} \right] + \\ &\qquad \frac{2}{J}\left( {{\sigma _{il}}{\delta _{jk}} + {\sigma _{jk}}{\delta _{il}} + {\sigma _{ik}}{\delta _{jl}} + {\sigma _{jl}}{\delta _{ik}}} \right) - \frac{4}{J}{\sigma _{ij}}{\delta _{kl}} \\[-12pt]\end{split} $$ (50) 式中
$\delta $ 为克罗内克函数,$ i,j,k,l,I,J,K,L = 1,2,3 $ .通过式(21)计算获得柯西应力为3 × 3的矩阵, 式(50)计算得到的雅可比矩阵是3 × 3 × 3 × 3的4阶张量. 在UMAT子程序中柯西应力STREE被保存为1 × 6的矢量, 雅可比矩阵DDSDDE为6 × 6矩阵, 因此需要做进一步变换. 其中柯西应力的对应方式为
$$ {\text{STREES}}[6] = \left[ {\begin{array}{*{20}{c}} 1&4&5 \\ {}&2&6 \\ {sym}&{}&3 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {11}&{12}&{13} \\ {}&{22}&{23} \\ {sym}&{}&{33} \end{array}} \right] $$ (51) 雅可比矩阵的对应关系为
DDSDDE(1,1) =
${\boldsymbol{\varPi}}_{1,1,1,1}^{\tau Z - J}$ , DDSDDE(2,2) =${\boldsymbol{\varPi}}_{2,2,2,2}^{\tau Z - J}$ DDSDDE(3,3) =
${\boldsymbol{\varPi}}_{3,3,3,3}^{\tau Z - J}$ , DDSDDE(1,2) =${\boldsymbol{\varPi}}_{1,1,2,2}^{\tau Z - J}$ DDSDDE(1,3) =
${\boldsymbol{\varPi}}_{1,1,3,3}^{\tau Z - J}$ , DDSDDE(2,3) =${\boldsymbol{\varPi}}_{2,2,3,3}^{\tau Z - J}$ DDSDDE(1,4) =
${\boldsymbol{\varPi}}_{1,1,1,2}^{\tau Z - J}$ , DDSDDE(2,4) =${\boldsymbol{\varPi}}_{2,2,1,2}^{\tau Z - J}$ DDSDDE(3,4) =
${\boldsymbol{\varPi}}_{3,3,1,2}^{\tau Z - J}$ , DDSDDE(4,4) =${\boldsymbol{\varPi}}_{1,2,1,2}^{\tau Z - J}$ DDSDDE(4,4) =
${\boldsymbol{\varPi}}_{1,2,1,2}^{\tau Z - J}$ , DDSDDE(1,5) =${\boldsymbol{\varPi}}_{1,1,1,3}^{\tau Z - J}$ DDSDDE(2,5) =
${\boldsymbol{\varPi}}_{2,2,1,3}^{\tau Z - J}$ , DDSDDE(3,5) =${\boldsymbol{\varPi}}_{3,3,1,3}^{\tau Z - J}$ DDSDDE(1,6) =
${\boldsymbol{\varPi}}_{1,1,2,3}^{\tau Z - J}$ , DDSDDE(2,6) =${\boldsymbol{\varPi}}_{2,2,2,3}^{\tau Z - J}$ DDSDDE(3,6) =
${\boldsymbol{\varPi}}_{3,3,2,3}^{\tau Z - J}$ , DDSDDE(5,5) =${\boldsymbol{\varPi}}_{1,3,1,3}^{\tau Z - J}$ DDSDDE(6,6) =
${\boldsymbol{\varPi}}_{2,3,2,3}^{\tau Z - J}$ , DDSDDE(4,5) =${\boldsymbol{\varPi}}_{1,2,1,3}^{\tau Z - J}$ DDSDDE(4,6) =
${\boldsymbol{\varPi}}_{1,2,2,3}^{\tau Z - J}$ , DDSDDE(5,6) =${\boldsymbol{\varPi}}_{1,3,2,3}^{\tau Z - J}$ DDSDDE下标对应满足表5所示关系.
表 5 雅可比矩阵下标对应关系Table 5. Correspondence between the subscripts of the Jacobian matrixa b 1 2 3 4 5 6 i k 1 2 3 1 1 2 j l 1 2 3 2 3 3 通过以上分析研究, 用户材料子程序UMAT中的雅可比矩阵DDSDDE下标通过数组的嵌套实现. 本文通过中间数组A实现雅可比矩阵DDSDDE的转化. 数组A为1 × 6的一维数组, 其中每个元素为一个含有两个元素的子数组, 如
${{A}}(4) = \left\{ {1,2} \right\}$ , 则${\text{DDSDD}}{{\text{E}}_{4,i}} = {\boldsymbol{\varPi}}_{1,2,k,l}^{\tau Z - J}$ , 通过这样的数据存储方式实现雅可比矩阵DDSDDE从3 × 3 × 3 × 3 4阶张量到6 × 6矩阵的转化.4. 骨骼肌多尺度数值模型的仿真分析和验证
4.1 骨骼肌多尺度数值模型材料参数的确定
肌纤维和结缔组织生物力学模型代表了肌纤维和结缔组织微观力学特性. 结缔组织的体积比与拉伸比拟合的线性关系式为
$$ {J^m} = 1 - 0.436\;09({\lambda ^m} - 1) $$ (52) 式中,
${J^m}$ 是结缔组织体积比;${\lambda ^m}$ 是结缔组织拉伸比. 肌纤维的体积比与拉伸比拟合的线性关系为$$ {J^f} = 1 - 0.025\;3({\lambda ^f} - 1) $$ (53) 式中,
${J^f}$ 是肌纤维体积比;${\lambda ^f}$ 是肌纤维拉伸比.为了区分肌纤维和结缔组织的生物力学模型的模型参数, 这里将第3节的生物力学模型改写为以下形式
$$ \begin{split} & {W^f} = \;\frac{{\kappa _0^f}}{2}{({J^f} - 1)^2} + \frac{{c_1^f}}{2}\left( {I_1^* - 3} \right) + \\ & \qquad \frac{{c_2^f}}{{2{k^f}}}\sum\limits_{i = 4,6} {\left\{ {\exp \left[ {{k^f}{{\left( {{I_i} - 1} \right)}^2}} \right] - 1} \right\}} + \\ &\qquad c_3^f\left( {{I_5} + {I_7} - 2{I_4} - 2{I_6} + 2} \right) \end{split} $$ (54) $$ \begin{split} & {W^m} = \;\frac{{\kappa _0^m}}{2}{({J^m} - 1)^2} + \frac{{c_1^m}}{2}\left( {I_1^* - 3} \right) + \\ & \qquad \frac{{c_2^m}}{{2{k^m}}}\sum\limits_{i = 4,6} {\left\{ {\exp \left[ {{k^m}{{\left( {{I_i} - 1} \right)}^2}} \right] - 1} \right\}} + \\ & \qquad c_3^m\left( {{I_5} + {I_7} - 2{I_4} - 2{I_6} + 2} \right) \end{split} $$ (55) 式中, 上标f表示肌纤维对应参数; 上标m表示结缔组织对应参数. 文中结合Flory[45]和Nolan等[46-47]中应变不变量与拉伸比的关系可得:
${I}_{1}^{*} = {J}^{-2/3}{I}_{1},$ $ {I_1} = \lambda _1^2 + \lambda _2^2 + \lambda _3^2 ,$ ${I_4} = {I_6} = \lambda _1^2{\cos^2 \theta } + \lambda _2^2{\sin^2 \theta },$ ${I_5} = {I_7} = \lambda _1^4{\cos^2 \theta } + \lambda _2^4{\sin^2 \theta}$ .本文选用Böl等[52]实验测量获得的猪后腿肌纤维名义应力-名义应变数据, 以及Smith等[25]实验测量的肌纤维体积比-拉伸比数据. 根据选取的实验数据, 利用Levenberg-Marquardt最小化算法与式(54)拟合获得肌纤维模型参数. 对于结缔组织生物力学模型的模型参数的确定, 选用Kohn等[4]通过实验获得猪后腿骨骼肌结缔组织的应力-应变数据, 以及Smith等[25]测量的肌纤维束体积比-应变数据. 同样选用Levenberg-Marquardt最小化算法与式(55)拟合获得结缔组织模型参数.
拟合过程在软件Origin2018实现, 在自定义函数中编辑式(54)整理后的C语言程序, 如下所示:
y = 1.0/2.0*k0*(x1-1.0)^2 + 1.0/2.0*c1*((x1^(-2.0/3.0)*(c2^2 + c3^2 + c4^2))-3)* + c2/k*exp(k*(x2^2*(cos(A))^2 + c3^2*(sin(A))^2)) + c3*(2*(c2^4(cos(A))^2 + C3^4*(sin(A))^2-4*(c2^2*(cos(A))^2 + c3^2*(sin(A))^2-2))).
在C语言程序中, 采用x1代表
${J^f}$ ; x2代表${\lambda _1}$ ; x3代表${\lambda _2}$ ; A代表角度${\theta ^f}$ ; k0表示${\kappa _0}^f$ ; c1, c2, c3分别表示${c_1}^f$ ,${c_2}^f$ 和${c_3}^f$ . 同理可得结缔组织C语言程序.拟合所得肌纤维和结缔组织的材料参数如表6所示, 其中肌纤维数据拟合的相关性系数为
${R^2} = 0.998$ , 结缔组织数据拟合的相关性系数为${R^2} = 0.999$ .表 6 肌纤维和结缔组织材料参数Table 6. Parameters of fibers and connective tissue$\kappa _0^f/ { {\text{MPa} } }$ $c_1^f/ { {\text{kPa} } }$ $c_2^f/ { {\text{kPa} } }$ 0.31444 7 160.44 $\kappa _0^m/ { {\text{MPa} } }$ $c_1^m/ { {\text{kPa} } }$ $c_2^m/ { {\text{kPa} } }$ 0.17074 3.84 106.39 $c_3^f/{\text{kPa} }$ ${k^f}$ ${\theta ^f}$ 5.6 0.0002643 1.5708 $c_3^m/{\text{kPa} }$ ${k^m}$ ${\theta ^m}$ 5.1 0.00132877 0.96 在结缔组织中, 胶原纤维与X轴夹角为0.96(弧度制), 由于胶原纤维分布于X-Y平面, 所以计算可得胶原纤维方向与Y轴夹角为34.97°, 此结果与Bleiler等[35]和Purslow等[53]获得的35°极为接近, 证明了数据拟合结果的有效性. 表6中的模型参数为参考值, 结合第1节中骨骼肌纵向拉伸实验测得数据, 进一步调整模型参数, 模型参数的最终结果如表7所示.
表 7 调整后的肌纤维和结缔组织材料参数Table 7. Adjusted material parameters of muscle fiber phase and connective tissue phase$\kappa _0^f/{ {\text{MPa} } }$ $c_1^f/ { {\text{kPa} } }$ $c_2^f/ { {\text{kPa} } }$ 0.31444 7 180.44 $\kappa _0^m/{ {\text{MPa} } }$ $c_1^m/ { {\text{kPa} } }$ $c_2^m/ { {\text{kPa} } }$ 0.17074 3.84 106.39 $c_3^f/{\text{kPa} }$ ${k^f}$ ${\theta ^f}(1)$ 5.6 1.5043 1.5708 $c_3^m/{\text{kPa} }$ ${k^m}$ ${\theta ^m}(1)$ 5.1 1.33 0.96 4.2 宏观-微观尺度之间的连接
在多尺度数值建模中, 宏观应力可以通过周期性边界条件和有限运动学中的均匀定理进行估计[32], RVE的体积平均柯西应力张量即为宏观柯西应力张量. 可得RVE的体积平均应力张量
$ {\boldsymbol{\hat \sigma }} $ 的计算公式为$$ {\boldsymbol{\hat \sigma }} = \frac{1}{V}\int_V {{\boldsymbol{\sigma }}{\rm{d}}V} $$ (56) 式中
$ {\boldsymbol{\sigma }} $ 为微观应力;$ V $ 为RVE体积. 通过平均柯西应力张量计算平均名义应力张量${\boldsymbol{\hat P}}$ 为$$ {{\hat {\boldsymbol{P}}}} = \hat J{{{\hat {\boldsymbol{F}}}}^{ - 1}}{{\hat {\boldsymbol{\sigma}} }} $$ (57) 式中
$\hat J = \det {{\hat {\boldsymbol{F}}}}$ ;${{\hat {\boldsymbol{F}}}}$ 为宏观变形梯度. 通过式(56)和式(57)建立宏观应力和微观应力之间的关系, 进而从微观响应获得宏观响应. 宏观名义应力张量的计算公式进一步表示为$$ {{\hat {\boldsymbol{P}}}} = \hat J{{{\hat {\boldsymbol{F}}}}^{ - 1}}\left( {\frac{1}{{{V_1}}}\int_{{V_1}} {{\sigma _1}} {\rm{d}}{V_1} + \frac{1}{{{V_2}}}\int_{{V_2}} {{\sigma _2}} {\rm{d}}{V_2}} \right) $$ (58) 式中, 下标1, 2分别代表肌纤维和结缔组织相关物理量. 为了书写方便, 下面采用
$ {\boldsymbol{\sigma }} $ 和$ {\boldsymbol{P}} $ 代表宏观柯西应力和名义应力.4.3 骨骼肌多尺度数值模型有限元仿真
本节使用第2节建立的RVE模型和第3节建立的生物力学模型, 结合4.1节确定的模型参数和4.2节中的宏观-微观尺度之间的转化方法共同完成多尺度数值模型的建立.
4.3.1 模型参数对多尺度数值模型仿真结果的影响分析
本小节针对肌纤维生物力学模型中的参数进行了灵敏度分析. 在对生物力学模型参数进行灵敏度分析时, 选取肌纤维体积分数为0.7的泰森多边形RVE模型为研究对象.
(1) 纵向拉伸
纵向拉伸中, 加载方向沿着肌纤维的分布方向, 也就是Y轴. 通过有限元仿真获得肌纤维的模型参数对纵向拉伸仿真结果的影响, 结果如图22所示.
![]()
图 22 模型参数对纵向拉伸仿真结果的影响Figure 22. The influence of model parameters on the simulation results of longitudinal stretch经过对图22(a) ~ 图22(e)的分析, 每个参数对名义应力-名义应变的影响范围和影响趋势清晰明了, 对每个模型参数对骨骼肌宏观生物力学特性的影响进行综合, 结果如图22(f)所示.
$ \kappa _0^f $ 的影响并不明显;$ c_1^f $ 的影响大于$ \kappa _0^f $ , 影响范围为整个变形区域;$ c_2^f $ 和${k^f}$ 的影响范围为$\varepsilon > 0.35$ 时;$ c_3^f $ 的影响范围与$ c_1^f $ 相同, 但是影响程度远大于$ c_1^f $ , 在踝区最为明显.(2) 横向拉伸
横向拉伸测试中, 在垂直于肌纤维方向添加拉伸位移载荷. 由于结缔组织中的胶原纤维分布在X-Y平面, 因此胶原纤维在横向拉伸过程中会抵抗拉伸位移载荷, 增加X轴方向的刚度. 仿真结果如图23所示.
![]()
图 23 模型参数对垂直于纤维方向拉伸仿真结果的影响Figure 23. The influence of model parameters on the simulation results of stretching longitudinal to fiber direction从仿真结果可得: 在垂直于肌纤维方向的位移载荷下, 肌纤维的模型参数
$ c_2^f $ ,${k^f}$ 和$ c_3^f $ 的改变不影响名义应力大小, 如图23(b)所示, 证明了仿真结果与理论分析结果一致.(3) 平面外纵向剪切
在平面外纵向剪切仿真测试中, 沿肌纤维方向添加剪切位移载荷. 仿真结果如图24(a)和图24(c)所示, 模型参数
$ c_1^f $ 和$ c_3^f $ 对名义剪切应力${\tau _{12}}$ 有明显影响; 如图24(b)所示, 模型参数$ c_2^f $ 不影响纵向剪切的力学行为, 改变${k^f}$ 得到了与改变$ c_2^f $ 相同的结论.仿真结果说明在平面外纵向剪切变形中, 肌纤维各向同性部分的参数影响剪切应力的大小. 但是应变不变量
$ I_4^f $ 的系数并不能影响剪切应力, 可见应变能方程中仅包含$I_1^*$ ,$ I_4^f $ ,$J$ 不能有效捕获材料的剪切变形行为, 也说明了在应变能方程中添加$ I_5^f $ 和$ I_7^f $ 是很有必要的. 模型参数对平面外纵向剪切变形的影响也证明了本文提出的生物力学模型的有效性.![]()
图 24 模型参数对平面外纵向剪切仿真结果的影响Figure 24. The influence of model parameters on the simulation results of out-of-plane longitudinal shear(4) 平面内剪切
对于平面内剪切, 剪切位移载荷方向垂直于肌纤维方向. 肌纤维模型参数对平面内剪切变形仿真结果的影响如图25所示.
![]()
图 25 模型参数对平面内剪切仿真结果的影响Figure 25. The influence of model parameters on the simulation results of in-plane shear从图25(a)中可以得出, 模型参数
$ c_1^f $ 对名义剪切应力${\tau _{13}}$ 有影响, 并且随着剪切应变的增加影响变得更加明显; 从图25(b)中可以得出, 随着$ c_2^f $ 变化, 名义剪切应力${\tau _{13}}$ 的仿真结果没有变化. 由于$ c_3^f $ 和${k^f}$ 对名义剪切应力${\tau _{13}}$ 的影响与$ c_2^f $ 相同, 因此$ c_3^f $ 和${k^f}$ 对平面内剪切仿真结果的影响曲线图不再展示.在剪切变形仿真测试中, 模型参数
$ c_3^f $ 和$ c_3^m $ 的值直接影响了应变不变量$ I_5^f $ 和$ I_5^m $ 贡献的大小. 此处设置$ c_3^f $ 和$ c_3^m $ 分别为0, 并将仿真结果与$ c_3^f = 0.005\;6 $ 和$ c_3^m = 0.005\;1 $ 的仿真结果进行对比, 结果如图26所示.![]()
图 26 应变不变量${I_5}$ 和${I_7}$ 对剪切变形仿真结果的影响Figure 26. Effects of strain invariants${I_5}$ and${I_7}$ on simulation results of shear deformation当
$ c_3^f = 0 $ 和$ c_3^m = 0 $ 时, 平面外纵向剪切的名义剪切应力${\tau _{12}}$ 曲线与平面内剪切的名义剪切应力${\tau _{13}}$ 曲线重合, 说明在生物力学模型中不含应变不变量${I_5}$ 和${I_7}$ 时, X-Y平面的剪切模量等于X-Z平面的剪切模量. 在肌纤维增强并且肌纤维体积分数为0.7的骨骼肌中, 这样的结果显然是不合理的. 因为在平面外纵向剪切中, 剪切力沿着肌纤维方向, 肌纤维也必然产生变形进而引起不同的剪切模量. 当$c_3^f = 0.005\;6$ 和$c_3^m = 0.005\;1$ 时, 平面外纵向剪切的名义剪切应力${\tau _{12}}$ 大于平面内剪切的名义剪切应力${\tau _{13}}$ . 这样的仿真结果与肌纤维增强的骨骼肌力学属性相匹配, 证明了本文提出的生物力学模型的有效性, 也进一步说明: 为了有效捕获纤维增强复合材料的剪切变形时的力学行为, 生物力学模型中必须包含应变不变量${I_5}$ 和${I_7}$ .4.3.2 体积分数对多尺度数值模型仿真结果影响的分析
本小节分别对肌纤维体积分数为0.7, 0.8和0.9的泰森多边形RVE模型进行仿真分析, 仿真结果如图27所示.
经过对比图27(a)和图27(b), 在拉伸中, 肌纤维体积分数的增加导致肌纤维方向的刚度增加, 分析原因: 肌纤维刚度大于结缔组织刚度, 肌纤维体积分数增加, 使骨骼肌总体刚度增加. 在横向拉伸中, 主要是结缔组织抵抗横向载荷, 肌纤维体积分数增加, 结缔组织的体积分数随之降低, 导致结缔组织抵抗作用减弱, 因此名义应力
${P_{11}}$ 随着体积分数的增加而减小.在平面外纵向剪切中, 名义剪切应变
${\tau _{12}}$ 与纵向拉伸的名义应力${P_{22}}$ 相同, 随着肌纤维体积分数的增加而增加, 如图27(c)所示. 产生这样的现象是由于肌纤维体积分数增加, 肌纤维增强作用增大, 引起名义剪切应变${\tau _{12}}$ 的增加.在平面内剪切中, 名义剪切应变
${\tau _{13}}$ 与肌纤维体积分数变化的关系并不显著, 如图27(d)所示. 体积分数为0.7和0.8的曲线近乎重合, 仅体积分数为0.9的曲线略有差别, 但是能够看出: 体积分数越大, 名义剪切应变${\tau _{13}}$ 越小. 肌纤维体积分数增加, 在X-Z平面内的结缔组织含量降低, 结缔组织中的胶原纤维含量随之降低, 对X-Z平面的增强作用减小, 引起名义剪切应变${\tau _{13}}$ 越小.4.3.3 肌纤维结构对多尺度数值模型仿真结果影响的分析
本小节分别对体积分数为0.7, 0.8和0.9的肌纤维结构分别为泰森多边形和曲边泰森多边形的RVE模型进行有限元仿真, 获得相应名义应力-名义应变关系如图28所示.
![]()
图 28 肌纤维结构对仿真结果影响Figure 28. The influence of muscle fiber structure on simulation results![]()
28 肌纤维结构对仿真结果影响(续)28. The influence of muscle fiber structure on simulation results (continued)经过对比图28(a)和图28(b), 图28(c)和图28(d), 发现无论是在拉伸测试中还是在剪切测试, 肌纤维的结构对垂直于纤维方向力学行为的影响大于对沿纤维方向力学行为的影响. 通过对比图28(a)和图28(c), 图28(b)和图28(d), 发现肌纤维的结构对拉伸变形的影响大于对剪切变形的影响. 在理论上, 肌纤维结构并不会影响纤维方向变形的仿真结果, 因为结缔组织和肌纤维沿着纤维方向以相同的应变发生变形. 但是, 在垂直于纤维方向, 肌纤维结构影响结缔组织分布, 影响名义应力大小分布, 因此肌纤维结构更多的影响垂直于纤维方向的力学行为. 最终可得仿真结果与理论分析结果一致.
经过以上分析, 肌纤维结构对垂直于肌纤维方向的力学行为影响较大, 为了进一步阐明曲边泰森多边形和泰森多边形对垂直于肌纤维方向的力学行为仿真精度的影响, 接下来将二者的横向拉伸和平面内剪切的仿真结果分别于实验数据进行对比. 选取体积分数为0.7的曲边泰森多边形RVE模型和泰森多边形RVE模型进行仿真计算, 获得对比结果如图29所示.
![]()
图 29 曲边泰森多边形和泰森多边形对仿真精度对比Figure 29. Comparison of simulation accuracy between curved-edge Voronoi polygons and Voronoi polygons在横向拉伸变形中, 对比结果如图29(a)所示. 当
$0 < \varepsilon < 0.32$ 时, 肌纤维截面为泰森多边形的仿真结果更接近实验曲线; 当$0.32 < \varepsilon < 0.6$ 时, 肌纤维截面为曲边泰森多边形的仿真结果更接近实验曲线. 为了对比二者精度, 计算仿真结果与实验数据的均方根误差(root mean squared error, RMSE)$$ {{RMSE}} = \sqrt {\frac{1}{N}\sum\limits_1^N {{{\left( {{X_{{\text{obs}},i}} - {X_{{\text{model}},i}}} \right)}^2}} } $$ (59) 式中N为数据点的个数;
$ {X_{{\text{obs}},i}} $ 为骨骼肌实验数据值;$ {X_{{\text{model}},i}} $ 为骨骼肌多尺度数值模型仿真值.肌纤维截面为泰森多边形的RMSE为0.00599 MPa, 肌纤维截面为曲边泰森多边形的RMSE为0.00525 MPa. 因此从整体计算结果可得曲边泰森多边形的仿真结果更加贴近实验数据, 说明肌纤维截面形状选择曲边泰森多边形时, 不仅使截面结构更加接近实验观察图像, 而且能够获得更高的仿真精度.
在平面内剪切变形中, 对比结果如图29(b)所示. 肌纤维截面结构为曲边泰森多边形的RVE仿真结果介于实验数据和肌纤维截面结构泰森多边形的RVE仿真结果之间, 说明肌纤维截面结构为曲边泰森多边形时的仿真精度更高.
4.4 骨骼肌多尺度数值模型的有效性验证
为了验证模型的有效性, 将纵向拉伸、横向拉伸、平面外纵向剪切和平面内剪切四种变形的仿真结果与第1节的骨骼肌准静态生物力学实验结果进行对比, 对比结果如图30所示.
![]()
图 30 多尺度数值模型仿真结果与实验数据的对比Figure 30. Comparison of simulation results of multiscale numerical model and experimental data在纵向拉伸中, 仿真结果曲线几乎与实验结果曲线重合. 但是在横向拉伸中, 仿真结果与实验数据存在一定的偏差. 由于实验样本中的结缔组织未完全剔除, 导致在平面外纵向剪切变形中, 两条曲线未完全重合; 而在平面内剪切变形中, 两曲线则完全重合.
仿真结果与实验数据的对比证明了本文建立的多尺度数值模型能够有效预测骨骼肌拉伸和剪切变形的力学行为. 从实验的角度和仿真的角度阐明了, 平面外纵向剪切和平面内剪切具有不同的剪切力学行为, 也就是在相同的剪切应变条件下具有不同的剪切应力.
为了进一步量化仿真结果和实验数据之间的误差, 分别计算纵向拉伸、横向拉伸、平面外纵向剪切和平面内剪切的仿真结果与实验结果的均方根误差. 纵向拉伸名义应力的RMSE为0.0035 MPa; 横向拉伸名义应力的RMSE为0.00599 MPa; 平面外纵向剪切的名义剪切应力的RMSE为0.00085 MPa; 平面内剪切的名义剪切应力的RMSE为0.000602 MPa. 经过对比分析, 在纵向拉伸中, RMSE是最大名义应力的2.1%; 在横向拉伸中, RMSE是最大名义应力的26.4%; 在平面外纵向剪切中, RMSE是最大名义剪切应力的7.1%; 在平面内剪切中, RMSE是最大名义剪切应力的8.1%. 纵向拉伸、平面外纵向剪切和平面内剪切的计算结果均在10%以内, 说明多尺度数值模型在这3种变形中能够精确预测骨骼肌生物力学行为.
4.5 骨骼肌多尺度数值模型对比验证
为了进一步验证本文提出的多尺度数值模型更加精确地预测骨骼肌宏观力学行为, 本小节将现有多尺度模型与本文建立的多尺度数值模型进行对比. 对比模型为Spyrou等[27]的多尺度模型. 由于Spyrou等[27]将他们的仿真结果与Morrow等[51]的实验进行了比较, 为统一标准, 本文也同Morrow等[54]的实验数据进行对比. 依照4.3.1节中多尺度数值模型参数的灵敏度分析, 调整获得的模型参数如表8所示. 本文多尺度数值模型与Spyrou等[27]的多尺度模型的对比结果如图31所示.
表 8 拟合Morrow等[54]的实验数据获得的模型参数Table 8. Model parameters obtained by fitting the experimental data of Morrow et al.[54]$\kappa _0^f/ { {\text{MPa} } }$ $c_1^f/{ {\text{kPa} } }$ $c_2^f/{ {\text{kPa} } }$ 11 9.16 3.6 $\kappa _0^m/{ {\text{MPa} } }$ $c_1^m/{ {\text{kPa} } }$ $c_2^m/ { {\text{kPa} } }$ 0.17074 3.84 1.0639 $c_3^f/ { {\text{kPa} } }$ ${k^f}\left( 1 \right)$ ${\theta ^f}(1)$ 15.15 3.2044 1.57 $c_3^m/{ {\text{kPa} } }$ ${k^m}\left( 1 \right)$ ${\theta ^m}(1)$ 5.1 1.33 0.436 ![]()
骨骼肌多尺度数值模型阐明了骨骼肌微观结构和组分生物力学模型对骨骼肌宏观力学行为的影响规律, 具有精度高和反应客观微观结构的优势. 以泰森多边形的RVE模型为例, 统计获得了泰森多边形RVE模型的网格数量和仿真所需时间, 如图32所示. 图中n表示RVE模型包含的网格数量.
![]()
图 32 泰森多边形RVE模型的网格数量和仿真时间Figure 32. Mesh number and simulation time of the Voronoi polygons RVE model从图32可知, 在相同网格数量、相同体积分数和相同计算机的条件下, 不同变形类型对应的计算时间存在明显差异, 且具有一定的大小关系, 即纵向拉伸 > 横向拉伸 > 平面外纵向剪切 > 平面内剪切. 说明多尺度数值模型在剪切变形时收敛速度快. 这样的大小关系在体积分数为0.9时表现最为明显. 经过分析, 当肌纤维体积分数达到0.9时, 结缔组织的厚度减小, 为了仿真能够顺利收敛, 划分网格数量增加, 因此时间差异更加明显. 而对于同一变形类型, 不同体积分数没有明显的计算时间大小关系.
5. 结 论
本文针对肌纤维微观结构模型与显微镜下观察的图像存在一定差异、微观组分生物力学模型无法有效捕获骨骼肌剪切变形时的力学行为、骨骼肌多尺度数值模型计算成本高的问题, 经过深入研究骨骼肌微观组分和结构, 确定影响骨骼肌宏观力学特性的内在因素. 在此基础上, 本文建立更加真实的骨骼肌微观RVE模型、建立肌纤维和结缔组织的生物力学模型. 结合骨骼肌准生物力学实验数据、RVE模型和生物力学模型实现多尺度数值模型的仿真, 得到以下结论.
(1)得出了骨骼肌具有各向异性、超弹性和不同平面具有不同剪切模量的生物力学特性. 完成了骨骼肌纵向拉伸, 横向拉伸, 平面外纵向剪切和平面内剪切实验.
(2)得出了肌纤维结构对骨骼肌宏观力学行为的影响. 分别建立了肌纤维横截面为泰森多边形和曲边泰森多边形的RVE模型. 从仿真结果分析可得, 在横向拉伸测试和平面内剪切测试中, 曲边泰森多边形结构的仿真结果均大于泰森多边形结构的仿真结果.
(3)得出了肌纤维体积分数对骨骼肌宏观力学行为的影响规律. 在多尺度数值模型中, 分别建立肌纤维的体积分数分别为0.7, 0.8和0.9的RVE模型, 并进行有限元仿真.
(4)建立了骨骼肌微观组分(肌纤维和结缔组织)生物力学模型, 有效预测了骨骼肌组分的生物力学行为. 通过仿真分析, 建立的生物力学模型能够捕获骨骼肌的各向异性、可压缩性和超弹性, 且在不同剪切平面获得不同剪切模量.
-
![]()
图 1 骨骼肌力学实验示意图
Figure 1. Schematic diagram of the mechanics experiments of skeletal muscle
![]()
图 4 拉伸实验名义应力-名义应变曲线
Figure 4. Nominal stress-nominal strain curve of stretch experiment
![]()
图 6 平面内剪切的名义剪切应力−剪切应变曲线
Figure 6. Nominal shear stress-shear strain curve of in-plane shear
![]()
图 7 平面外纵向剪切的名义剪切应力−剪切应变曲线
Figure 7. Nominal shear stress-shear strain curve of out-of-plane longitudinal shear
![]()
图 20 骨骼肌肌纤维和胶原纤维的分布
Figure 20. Distribution of skeletal muscle fibers and collagen fibers
![]()
图 22 模型参数对纵向拉伸仿真结果的影响
Figure 22. The influence of model parameters on the simulation results of longitudinal stretch
![]()
图 23 模型参数对垂直于纤维方向拉伸仿真结果的影响
Figure 23. The influence of model parameters on the simulation results of stretching longitudinal to fiber direction
![]()
图 24 模型参数对平面外纵向剪切仿真结果的影响
Figure 24. The influence of model parameters on the simulation results of out-of-plane longitudinal shear
![]()
图 25 模型参数对平面内剪切仿真结果的影响
Figure 25. The influence of model parameters on the simulation results of in-plane shear
![]()
图 26 应变不变量
${I_5}$ 和${I_7}$ 对剪切变形仿真结果的影响Figure 26. Effects of strain invariants
${I_5}$ and${I_7}$ on simulation results of shear deformation![]()
图 28 肌纤维结构对仿真结果影响
Figure 28. The influence of muscle fiber structure on simulation results
![]()
28 肌纤维结构对仿真结果影响(续)
28. The influence of muscle fiber structure on simulation results (continued)
![]()
图 29 曲边泰森多边形和泰森多边形对仿真精度对比
Figure 29. Comparison of simulation accuracy between curved-edge Voronoi polygons and Voronoi polygons
![]()
图 30 多尺度数值模型仿真结果与实验数据的对比
Figure 30. Comparison of simulation results of multiscale numerical model and experimental data
![]()
![]()
图 32 泰森多边形RVE模型的网格数量和仿真时间
Figure 32. Mesh number and simulation time of the Voronoi polygons RVE model
表 1 拉伸实验试样几何参数
Table 1 Geometric parameters of specimens of stretch experiment
Experiment No. L/mm W/mm H/mm skeletal muscle longitudinal
stretch experiment1-1 28.01 10.07 4.20 1-2 33.40 8.24 4.48 1-3 33.59 7.77 4.51 1-4 29.99 9.45 4.78 1-5 30.01 10.11 4.53 skeletal muscle lateral stretch experiment 2-1 31.98 11.94 6.26 2-2 36.25 7.27 5.63 2-3 28.71 7.48 5.77 2-4 38.10 13.13 5.21 2-5 33.74 12.03 5.23 
下载: 导出CSV
表 2 剪切实验试样几何参数
Table 2 Geometric parameters of specimens of shear experiment
Experiment No. L/mm W/mm H/mm skeletal muscle in-plane shear
experiment3-1 19.30 21.85 5.30 3-2 18.46 21.79 5.69 3-3 18.62 18.60 4.02 3-4 18.07 21.14 4.40 3-5 18.65 21.40 4.50 skeletal muscle out-of-plane
longitudinal shear experiment4-1 22.84 16.50 4.11 4-2 20.75 20.16 5.18 4-3 26.37 15.81 4.65 4-4 22.87 23.25 3.31 4-5 17.98 17.74 2.70
下载: 导出CSV
表 4 ep设置值
Table 4 Setting value of ep
Volume fraction Voronoi polygons/μm Curved-edges Voronoi
polygons/μm0.7 0.0080 0.005229 0.8 0.0052 0.003000 0.9 0.0025 0.000540
下载: 导出CSV
表 3 森多边形和曲边泰森多边形RVE模型
Table 3 The RVE model of Voronoi polygons and curved-edge Voronoi polygons
Volume fraction Voronoi polygons RVE model Curved-edges Voronoi polygons RVE model 0.7 0.8 0.9
下载: 导出CSV
表 5 雅可比矩阵下标对应关系
Table 5 Correspondence between the subscripts of the Jacobian matrix
a b 1 2 3 4 5 6 i k 1 2 3 1 1 2 j l 1 2 3 2 3 3
下载: 导出CSV
表 6 肌纤维和结缔组织材料参数
Table 6 Parameters of fibers and connective tissue
$\kappa _0^f/ { {\text{MPa} } }$ $c_1^f/ { {\text{kPa} } }$ $c_2^f/ { {\text{kPa} } }$ 0.31444 7 160.44 $\kappa _0^m/ { {\text{MPa} } }$ $c_1^m/ { {\text{kPa} } }$ $c_2^m/ { {\text{kPa} } }$ 0.17074 3.84 106.39 $c_3^f/{\text{kPa} }$ ${k^f}$ ${\theta ^f}$ 5.6 0.0002643 1.5708 $c_3^m/{\text{kPa} }$ ${k^m}$ ${\theta ^m}$ 5.1 0.00132877 0.96
下载: 导出CSV
表 7 调整后的肌纤维和结缔组织材料参数
Table 7 Adjusted material parameters of muscle fiber phase and connective tissue phase
$\kappa _0^f/{ {\text{MPa} } }$ $c_1^f/ { {\text{kPa} } }$ $c_2^f/ { {\text{kPa} } }$ 0.31444 7 180.44 $\kappa _0^m/{ {\text{MPa} } }$ $c_1^m/ { {\text{kPa} } }$ $c_2^m/ { {\text{kPa} } }$ 0.17074 3.84 106.39 $c_3^f/{\text{kPa} }$ ${k^f}$ ${\theta ^f}(1)$ 5.6 1.5043 1.5708 $c_3^m/{\text{kPa} }$ ${k^m}$ ${\theta ^m}(1)$ 5.1 1.33 0.96
下载: 导出CSV
表 8 拟合Morrow等[54]的实验数据获得的模型参数
Table 8 Model parameters obtained by fitting the experimental data of Morrow et al.[54]
$\kappa _0^f/ { {\text{MPa} } }$ $c_1^f/{ {\text{kPa} } }$ $c_2^f/{ {\text{kPa} } }$ 11 9.16 3.6 $\kappa _0^m/{ {\text{MPa} } }$ $c_1^m/{ {\text{kPa} } }$ $c_2^m/ { {\text{kPa} } }$ 0.17074 3.84 1.0639 $c_3^f/ { {\text{kPa} } }$ ${k^f}\left( 1 \right)$ ${\theta ^f}(1)$ 15.15 3.2044 1.57 $c_3^m/{ {\text{kPa} } }$ ${k^m}\left( 1 \right)$ ${\theta ^m}(1)$ 5.1 1.33 0.436
下载: 导出CSV
-
[1] 陈胜国, 汪华侨. 人体骨骼肌的分布和变异. 解剖学研究, 2013, 35(3): 237-240 (Chen Shengguo, Wang Huaqiao. Distribution and variation of human skeletal muscle. Anatomy Research, 2013, 35(3): 237-240 (in Chinese) [2] Janssen I, Heymsfield S, Wang Z, et al. Skeletal muscle mass and distribution in 468 men and women aged 18-88 Yr. Journal of Applied Physiology, 2000, 89(1): 81-88 doi: 10.1152/jappl.2000.89.1.81
[3] 姜宗来. 从生物力学到力学生物学的进展. 力学进展, 2017, 47: 313-336 (Jiang Zonglai. Advances from biomechanics to mechanical biology. Advances in Mechanics, 2017, 47: 313-336 (in Chinese) doi: 10.6052/1000-0992-16-023 [4] Fung YC. Biomechanics: mechanical properties of living tissues. Journal of Biomechanical Engineering, 1981, 103(4): 231-298 doi: 10.1115/1.3138285
[5] Spencer A. Deformations of Fibre-reinforced Materials. New York: Oxford University Press, 1972
[6] Gerard JM, Ohayon J, Luboz V, et al. Non-Linear elastic properties of the lingual and facial tissues assessed by indentation technique application to the biomechanics of speech production. Medical Engineering & Physics, 2005, 27(10): 884-892
[7] 陈伟, 吴立军, 严志汉等. 老年健康女性盆底肛提肌有限元模型的建立及意义. 生物医学工程学杂志, 2011, 28(5): 927-931 (Chen Wei, Wu Lijun, Yan Zhihan, et al. Establishment and significance of finite element model of levator anus pelvic floor muscle in elderly healthy women. Journal of Biomedical Engineering, 2011, 28(5): 927-931 (in Chinese) [8] Schiavone P, Boudou T, Promayon E, et al. A light sterilizable pipette device for the in vivo estimation of human soft tissues constitutive laws//International Conference of the IEEE Engineering in Medicine & Biology Society, Vancouver, 2008: 4298-4310
[9] Nazari MA, Perrier P, Chabanas M, et al. Simulation of dynamic orofacial movements using a constitutive law varying with muscle activation. Computer Methods in Biomechanics and Biomedical Engineering, 2010, 13(4): 469-482 doi: 10.1080/10255840903505147
[10] Nazari MA, Perrier P, Chabanas M, et al. Shaping by stiffening: a modeling study for lips. Motor Control, 2011, 15(1): 141-168 doi: 10.1123/mcj.15.1.141
[11] Stavness I, Lloyd JE, Fels S. Automatic prediction of tongue muscle activations using a finite element model. Journal of Biomechanics, 2012, 45(16): 2841-2848 doi: 10.1016/j.jbiomech.2012.08.031
[12] Gindre J, Takaza M, Moerman KM, et al. A structural model of passive skeletal muscle shows two reinforcement processes in resisting deformation. Journal of the Mechanical Behavior of Biomedical Materials, 2013, 22: 84-94 doi: 10.1016/j.jmbbm.2013.02.007
[13] Voigt W. Ueber die beziehung zwischen den beiden elasticitätsconstanten isotroper körper. Annalen der Physik (Leipzig)
, 1889, 274(12): 573-587 doi: 10.1002/andp.18892741206 [14] Ogden RW. Nonlinear elasticity, anisotropy, material stability and residual stresses in soft tissue. Biomechanics of Soft Tissue in Cardiovascular Systems, 2003, 441: 65-108
[15] Al-Dirini RM, Reed MP, Hu J, et al. Development and validation of a high anatomical fidelity fe model for the buttock and thigh of a seated individual. Annals of Biomedical Engineering, 2016, 44(9): 2805-2816 doi: 10.1007/s10439-016-1560-3
[16] Roux A, Laporte S, Lecompte J, et al. Influence of muscle-tendon complex geometrical parameters on modeling passive stretch behavior with the discrete element method. Journal of Biomechanics, 2016, 49(2): 252-258 doi: 10.1016/j.jbiomech.2015.12.006
[17] Bosboom EMH, Hesselink MKC, Oomens CWJ, et al. Passive lateral mechanical properties of skeletal muscle under in vivo compression. Journal of Biomechanics, 2001, 34(10): 1365-1368 doi: 10.1016/S0021-9290(01)00083-5
[18] 粟思橙. 基于肌肉主动力的颈部有限元建模研究. [硕士论文]. 长沙: 湖南大学, 2014 Su Sicheng. Research on Finite Element Modeling of Neck Based on Muscle Initiative. [Master Thesis]. Changsha: Hunan University, 2014 (in Chinese))
[19] Sengeh DM, Moerman KM, Petron A, et al. Multi-material 3-D viscoelastic model of a transtibial residuum from in-vivo indentation and mri data. Journal of the Mechanical Behavior of Biomedical Materials, 2016, 59: 379-392 doi: 10.1016/j.jmbbm.2016.02.020
[20] Moerman KM, Simms CK, Nagel T. Control of stretch-compression asymmetry in ogden hyperelasticity with application to soft tissue modelling. Journal of the Mechanical Behavior of Biomedical Materials, 2016, 56: 218-228 doi: 10.1016/j.jmbbm.2015.11.027
[21] Holzapfel GA, Gasser TC, Ogden RW. A new constitutive framework for arterial wall mechanics and a comparative study of material models. Journal of Elasticity, 2000, 61(1/3): 1-48 doi: 10.1023/A:1010835316564
[22] Hernández B, Pena E, Pascual G, et al. Mechanical and histological characterization of the abdominal muscle. A Previous Step to Modelling Hernia Surgery. J Mech Behav Biomed Mater, 2011, 4(3): 392-404
[23] Böl M, Weikert R, Weichert C. A coupled electromechanical model for the excitation-dependent contraction of skeletal muscle. Journal of the Mechanical Behavior of Biomedical Materials, 2011, 4(7): 1299-1310 doi: 10.1016/j.jmbbm.2011.04.017
[24] Wu T, Hung APL, Hunter P, et al. Modelling facial expressions: a framework for simulating nonlinear soft tissue deformations using embedded 3D muscles. Finite Elements in Analysis and Design, 2013, 76: 63-70 doi: 10.1016/j.finel.2013.08.002
[25] Calvo B, Ramírez A, Alonso A, et al. Passive nonlinear elastic behaviour of skeletal muscle: experimental results and model formulation. Journal of Biomechanics, 2010, 43(2): 318-325 doi: 10.1016/j.jbiomech.2009.08.032
[26] Spyrou LA, Agoras M, Danas K. A homogenization model of the voigt type for skeletal muscle. Journal of Theoretical Biology, 2017, 414: 50-61 doi: 10.1016/j.jtbi.2016.11.018
[27] Spyrou LA, Brisard S, Danas K. Multiscale modeling of skeletal muscle tissues based on analytical and numerical homogenization. Journal of the Mechanical Behavior of Biomedical Materials, 2019, 92: 97-117 doi: 10.1016/j.jmbbm.2018.12.030
[28] Sharafi B, Blemker SS. A micromechanical model of skeletal muscle to explore the effects of fiber and fascicle geometry. Journal of Biomechanics, 2010, 43: 3207-3213 doi: 10.1016/j.jbiomech.2010.07.020
[29] Sharafi B, Blemker SS. A mathematical model of force transmission from intrafascicularly terminating muscle fibers. Journal of Biomechanics, 2011, 44(11): 2031-2039 doi: 10.1016/j.jbiomech.2011.04.038
[30] Jiménez FL. Modeling of soft composites under three-dimensional loading. Composites Part B-Engineering, 2014, 59: 173-180 doi: 10.1016/j.compositesb.2013.11.020
[31] Virgilio KM, Martin KS, Peirce SM, et al. Multiscale models of skeletal muscle reveal the complex effects of muscular dystrophy on tissue mechanics and damage susceptibility. Interface Focus, 2015, 5(2): 20140080 doi: 10.1098/rsfs.2014.0080
[32] Kuravi R, Leichsenring K, Böl M, et al. 3D finite element models from serial section histology of skeletal muscle tissue-the role of micro-architecture on mechanical behaviour. Journal of the Mechanical Behavior of Biomedical Materials, 2021, 113: 104109 doi: 10.1016/j.jmbbm.2020.104109
[33] Valentin T, Simms C. An inverse model of the mechanical response of passive skeletal muscle: implications for microstructure. Journal of Biomechanics, 2020, 99: 109483 doi: 10.1016/j.jbiomech.2019.109483
[34] Kuravi R, Leichsenring K, Trostorf R, et al. Predicting muscle tissue response from calibrated component models and histology-based finite element models. Journal of the Mechanical Behavior of Biomedical Materials, 2021, 117: 104375 doi: 10.1016/j.jmbbm.2021.104375
[35] Bleiler C, Ponte CP, Rohrle O. A microstructurally-based, multi-Scale, continuum-mechanical model for the passive behaviour of skeletal muscle tissue. Journal of the Mechanical Behavior of Biomedical Materials, 2019, 97: 171-186 doi: 10.1016/j.jmbbm.2019.05.012
[36] 王礼立. 高应变率下材料动态力学性能. 力学与实践, 1982, 1: 9-19 (Wang Lili. Dynamic mechanical properties of materials under high strain rate. Mechanics in Engineering, 1982, 1: 9-19 (in Chinese) [37] Nemat-Nasser S, Lori MSKD. Micromechanics: overall properties of heterogeneous materials. Journal of Applied Mechanics, 1996, 63(2): 561 doi: 10.1115/1.2788912
[38] Wisdom KM, Delp SL, Kuhl E. Use it or lose it: multiscale skeletal muscle adaptation to mechanical stimuli. Biomech Model Mechanobiol, 2015, 14(2): 195-215 doi: 10.1007/s10237-014-0607-3
[39] Kanit T, Forest S, Galliet I, et al. Determination of the size of the representative volume element for random composites: statistical and numerical approach. International Journal of Solids and Structures, 2003, 40(13-14): 3647-3679 doi: 10.1016/S0020-7683(03)00143-4
[40] 武鹏伟. 非均质材料微—宏观非线性分析的多尺度研究. [硕士论文]. 杭州: 浙江工业大学, 2016 Wu Pengwei. Multi-scale study of micro-macroscopic nonlinear analysis of heterogeneous materials. [Master Thesis]. Hangzhou: Zhejiang University of Technology, 2016 (in Chinese))
[41] 李庆, 杨晓翔. 周期性边界条件下炭黑增强橡胶基复合材料有效弹性性能数值模拟. 福州大学学报(自然科学版), 2013, 41(1): 97-103 (Li Qing, Yang Xiaoxiang. Numerical simulation of effective elastic properties of carbon black reinforced rubber matrix composites under periodic boundary conditions. Journal of Fuzhou University (Natural Science Edition, 2013, 41(1): 97-103 (in Chinese) [42] Xia Z, Zhang Y, Ellyin F. A unified periodical boundary conditions for representative volume elements of composites and applications. International Journal of Solids and Structures, 2003, 40(8): 1907-1921 doi: 10.1016/S0020-7683(03)00024-6
[43] Xu R, Bouby C, Zahrouni H, et al. 3D modeling of shape memory alloy fiber reinforced composites by multiscale finite element method. Composite Structures, 2018, 200: 408-419 doi: 10.1016/j.compstruct.2018.05.108
[44] Holzapfel GA. Nonlinear solid mechanics: a continuum approach for engineering. Chichester: John Wiley & Sons, 2000, 37: 489-490
[45] Flory PJ. Thermodynamic relations for high elastic materials. Transactions of the Faraday Society, 1961, 57(5): 829-838
[46] Nolan DR, Gower AL, Destrade M, et al. A robust anisotropic hyperelastic formulation for the modelling of soft tissue. Journal of the Mechanical Behavior of Biomedical Materials, 2014, 39: 48-60 doi: 10.1016/j.jmbbm.2014.06.016
[47] Nolan DR, Mcgarry JP. On the compressibility of arterial tissue. Annals of Biomedical Engineering, 2016, 44(4): 993-1007 doi: 10.1007/s10439-015-1417-1
[48] Klisch SM. A bimodular polyconvex anisotropic strain energy function for articular cartilage. Journal of Biomechanical Engineering, 2007, 129(2): 250-258 doi: 10.1115/1.2486225
[49] Schrder J, Neff P. Invariant formulation of hyperelastic lateral isotropy based on polyconvex free energy functions. International Journal of Solids & Structure, 2003, 40(2): 401-445
[50] Sun W, Chaikof EL, Levenston ME. Numerical approximation of tangent Moduli for finite element implementations of nonlinear hyperelastic material models. Journal of Biomechanical Engineering, 2008, 130(6): 061003
[51] Bazant ZP, Gattu M, Vorel J. Work conjugacy error in commercial finite-element codes: its magnitude and how to compensate for it. Proceedings of the Royal Society A-Mathematical Physical and Engineering Sciences, 2012, 468(2146): 3047-3058 doi: 10.1098/rspa.2012.0167
[52] Böl M, Iyer R, Dittmann J, et al. Investigating the passive mechanical behaviour of skeletal muscle fibres: micromechanical experiments and bayesian hierarchical modelling. Acta Biomaterialia, 2019, 92: 277-289 doi: 10.1016/j.actbio.2019.05.015
[53] Purslow PP. Muscle fascia and force transmission. Journal of Bodywork and Movement Therapies, 2010, 14(4): 411-417 doi: 10.1016/j.jbmt.2010.01.005
[54] Morrow DA, Haut Donahue TL, Odegard GM, et al. Laterally isotropic stretch material properties of skeletal muscle tissue. Journal of the Mechanical Behavior of Biomedical Materials, 2010, 3(1): 124-129 doi: 10.1016/j.jmbbm.2009.03.004
计量
- 文章访问数: 826
- HTML全文浏览量: 382
- PDF下载量: 181
