[发明专利]一种基于微观结构参数个体化关节软骨仿真方法

摘要

一种基于微观结构参数个体化关节软骨仿真方法，本发明涉及基于微观结构参数个体化关节软骨仿真方法。本发明是为了解决现有获取专门患者关节软骨参数的手段为生理切片等破坏性手段的问题。本发明步骤为步骤一关节软骨微观结构参数赋值；步骤二根据步骤一建立个体化线弹性关节软骨微观模型；步骤三将粘弹性材料参数整合到步骤二建立的个体化线弹性关节软骨微观模型中，建立粘弹性关节软骨微观模型；步骤四对步骤三建立的粘弹性关节软骨微观模型进行有限元仿真。本发明应用于生物医学工程领域。

主权项

一种基于微观结构参数个体化关节软骨仿真方法，其特征在于，所述基于微观结构参数个体化关节软骨仿真方法包括以下步骤：步骤一：关节软骨微观结构参数赋值；步骤二：根据步骤一建立个体化线弹性关节软骨微观模型；步骤三：将粘弹性材料参数整合到步骤二建立的个体化线弹性关节软骨微观模型中，建立粘弹性关节软骨微观模型；步骤四：对步骤三建立的粘弹性关节软骨微观模型进行有限元仿真；所述步骤一中关节软骨微观结构参数赋值具体为：基于微观结构参数个体化关节软骨模型参数分为两类：一类是组分材料参数，关节软骨的组分分为纤维和基质两部分，与粘弹性纤维相关联的材料参数是λF,μF，与弹性基质相关联的材料参数为λM,μM；λF为纤维拉梅常数，μF为纤维剪切模量，λM为基质拉梅常数，μM为基质剪切模量；第二类参数是专门患者结构参数，包括特征体元的纤维体积分数VF和特征体元在整个关节软骨中纤维取向分布函数Φ(θ,φ)；第一类组分材料参数通过力学试验获取；第二类结构参数获取包括以下步骤：步骤一一：关节软骨微观结构参数的选取；选择纤维体积分数和纤维取向分布函数作为建立个体化关节软骨模型依据的微观结构参数，所述纤维体积分数为纤维所占面积与观察面积之比，纤维取向分布函数是纤维三维空间分布的一种方位表达形式；步骤一二：关节软骨微观结构参数的计算；基于纤维占光学图像的比例确定纤维体积分数，分为两个步骤；步骤一二一：应用全局算法和局部算法，选择阈值使光学成像分为纤维和背景两个区域；步骤一二二：纤维占据图像区域的计算；VF=Σk=1k=NAFKNAP]]>式中：为第k层图像中纤维所占的面积；N为断层图像的总数；AP为图像的总面积；步骤一三：采用基于正交滤波器的方法，通过滤波器定义获取滤波输出及基于取向张量获得纤维取向分布，确定纤维取向分布函数；具体确定方法如下：方向函数为：式中：是常规化的频率向量，为定义滤波器k方向的单位向量；取向张量的建立需要滤波器的输出qk，即图像与每个滤波器卷积后的结果；滤波器Fk定义在频域上，通过卷积定理得到滤波器的输出qk为：式中：为傅里叶逆变换，I′(ω)为与图像相关的像素强度函数；R(ω)为滤波器的径向函数；得到滤波器输出qk后，每个像素处的取向张量Τ由下式得到：T=Σk||qk||Mk]]>式中：I为单位矩阵；Τ有3个特征值e1≥e2≥e3和3个相应的特征向量通过对特征值大小的比较，得到确定方向的方法：a)e1≈e2＞＞e3表示线性特征；b)e1＞＞e2≈e3表示平面特征；c)e1≈e2≈e3表示没有确定取向的各向同性区域；在球面坐标系中特征向量的表达式为：式中：和分别为在X，Y和Z上的分量；从每个像素处估计的取向张量中转换出在球面坐标系中表示的纤维取向的直方图，直方图常规化为纤维取向分布函数；所述步骤二中建立个体化线弹性关节软骨微观模型的具体过程为：步骤二一：基于微观结构参数线弹性关节软骨微观模型的建立；关节软骨微观力学模型特征体元为同心圆柱模型，在特征体元上建立直角坐标系(1,2,3)，1沿着纤维轴向方向，称为特征体元的径向，2和3位于垂直于纤维的平面内，称为特征体元的横向；由1、2和3建立的坐标系称为复合材料的局部坐标系，全局坐标系放置于复合材料的中心位置；每一个特征体元由一对同轴的圆柱体组成，内部圆柱为内部纤维，外部圆柱为外部基质；设内部纤维和外部基质是线弹性的，线弹性模型本构方程表述为：σ＝λeΙ+2με＝Cε式中，σ为应力，ε为应变，e为体积应变，λ,μ是弹性模量E和泊松比ν的函数；λ=Ev(1+v)(1-2v)]]>μ=E2(1+v)]]>纤维的线弹性特性由λF,μF确定，基质的线弹性特性由λM,μM确定；步骤二二：计算刚度矩阵；对于每一个特征体元，横观各向同性刚度矩阵为：Cijkl=C[11]C[12]C[12]000C[22]C[23]000C[22]000C[22]-C[23]200SymC[55]0C[55]]]>式中，Sym代表对称矩阵中的对称项；将上式中的各系数替换为工程常数的表达式如下：Cijkl=E11+4v122K232K23v122K23v12μ23+K23K23-μ23μ23+K23μ23Symμ12μ12]]>式中：E11为纵向的杨氏模量；v12为纵向泊松比；K23为平面应变体积模量；μ12为面内剪切模量；μ23为特征体元的横观剪切模量；表征的工程常数与各组分材料参数之间的关系式，其表达形式如下：E11=VFμF(3λF+2μF)ξ1+(1-VF)μM(3λM+2μM)ξ2+(1-VF)VFμMξ32ξ4]]>v12=VMμM+VFμF+VFVM(vF-vM)(33KM+μM-33KF+μF)3VM3KF+μF+3VF3KM+μM+1μM]]>K23=ξ2+VF1ξ1-ξ2+1-VFλM+2μM]]>μ12=μM[(1+VF)μF+(1-VF)μM](1-VF)μF+(1+VF)μM]]>μ23=μM+VFμMμMμF-μM+(1-VF)(λM+3μM)2(λM+2μM)]]>式中：ξ1＝λF+μFξ2＝λM+μMξ3=λFξ1-λMξ2]]>ξ4=1+VFμFξ2+(1-VF)μMξ1]]>Vα，Eα，να，μα和Kα分别为α相的体积分数，弹性模量，泊松比，剪切模量和体积模量；则刚度矩阵中的各刚度系数表达式为：C[11]=(λM+2μM)(1-VF)+VF[λF+2μF-(1-VF)(λF-λM)2(1-VF)ξ1+VFξ2+μM]]]>C[12]=(1-VF)λM(μF+μM)+λF(λM+2VFμM)(1-VF)ξ1+VFξ2+μM]]>C[22]＝λM+2μM+ξ5+ξ6C[23]＝λM+ξ5‑ξ6C[55]=μ12=μM[(1+VF)μF+(1-VF)μM](1-VF)μF+(1-VF)μM]]>式中：ξ5=VF1ξ1-ξ2+1-VFλM+2μM]]>ξ6=VFμMμMμF-μM+(1-VF)(λM+3μM)2(λM+2μM)]]>本构方程在全局坐标系表达为下式的：σ′ij＝C′ijklε'klC′ijkl张量在全局坐标系中通过θ和φ角度描述单个特征体元的刚度；设整个组织是由多个特征体元以随机分布的形式组成的，组织的整体刚度通过对θ和φ所有角度范围内的C′ijkl积分得到；则组织作为整体的等效刚度矩阵由下式给出：C‾ijkl=∫0π∫0πΦ(θ,φ)Cijkl′sinφdφdθ]]>所述步骤三中建立粘弹性关节软骨微观模型的具体过程为：步骤三一：确定基于个体化微观结构参数粘弹性关节软骨粘弹性参数；将粘弹性材料参数整合进个体化线弹性关节软骨微观模型中，采用Prony级数形式替换个体化线弹性关节软骨微观模型中的弹性参数；在个体化线弹性关节软骨微观模型中，纤维相的参数是λF,μF，利用粘弹性Prony级数参数替换μF；粘弹性Prony级数表示为：Ω(t)=Ω∞+ΣK=1NΩKe-t/τK]]>其中Ω∞为平衡模量，ΩK是松弛模量，τK是松弛时间常数，t为时间；步骤三二：确定应力应变关系矩阵；σij(t)=C‾ijkl0ϵkl-ΣK=1NDijK]]>式中，将Cijklεkl划分为时间相关项和时间无关项两个部分，为时间无关项的系数，为时间相关项；时间相关相定义：DijK=∫0tC‾ijklK(1-e(t′-t)/xK)×dϵkldt′dt′]]>式中，为全局等效刚度矩阵中与时间相关的项；所述步骤四中对步骤三建立的粘弹性关节软骨微观模型进行有限元仿真的输入变量、输出变量和计算公式为：输入变量：材料参数：λF、μF、λM和μM，松弛时间常数τK；结构参数：纤维体积分数VF和纤维取向分布函数Φ(θ,φ)；输出变量：用一致切线算子表示关节软骨力学特性，一致切线算子为计算公式：计算全局等效刚度矩阵：C‾ijkl=∫0x∫0xΦ(θ,φ)Cijkl′sinφdφdθ]]>更新时间相关项：Dijn+1K=DijnKΔt+C‾ijklKϵkln+1τK1Δt+1τK]]>更新包含时间相关项和时间无关项的应力应变等式：σijn+1(t)=C‾ijkl0ϵkln+1-ΣK=1NDijn+1K]]>更新一致切线算子为：∂σijn+1∂Δϵkl=C‾ijkl0-ΣK=1NΔtΔt+τKC‾ijklK.]]>