本发明涉及一种骨结构的预测方法，具体说是一种考虑应力状态的骨结构预测方法。
背景技术：
骨结构预测(也称为骨重建数值模拟)是研究骨重建理论的重要方法，为临床上预测骨密度变化或假体优化设计提供理论依据，也有助于各种骨代谢疾病的研究和新的骨代谢疾病治疗方案的提出。骨结构预测方程是“需要就生长，不需要就吸收”的骨功能适应性原理的数值实现算法。从工程观点看，是把骨骼作为一种工程材料的自适应优化设计。因此，骨功能适应性原理也称为骨结构自优化理论。
骨结构预测控制方程一般为：
当前骨结构预测方法存在的问题和缺点是力学激励无论采用应变能密度、米塞斯等效应力，还是米塞斯等效应变等均未考虑在骨结构预测过程中应力状态的变化，在整个计算过程中采用固定公式计算力学刺激。当采用相当应力作为力学激励进行骨结构预测时，相当应力的计算公式应随着应力状态的变化进行改变。当前的骨结构预测算法从力学原理来看是存在缺陷的。
技术实现要素：
由于现有的骨结构预测方法未考虑应力状态对力学激励的影响，本发明的目的是提供一个考虑应力状态的骨结构预测方法，使骨结构预测既能体现出骨作为一种有生命器官的骨功能适应性原理，又能体现出把骨骼作为一种工程材料的自适应优化设计。
为了实现上述目的，本发明方案的技术要点是：一种考虑应力状态的骨结构预测方法，包括：
步骤一：确定要进行预测的骨结构的拓扑范围，并根据力学原理抽象出力学边界条件和载荷；
步骤二：对步骤一确定的骨结构进行建模，建模可采用cad软件或有限元分析软件完成；
步骤三：对步骤二建立的骨结构模型进行有限元网格划分，得到骨结构有限元网格模型；
步骤四：定义骨结构的初始表观密度，并根据弹性模量与表观密度的关系式计算骨的材料性能参数弹性模量，泊松比取为0.3，并对每个单元定义材料属性；
e＝cργ------------------------------------------(1)
式中e为弹性模量，ρ为骨的表观密度，单位体积骨的质量，用来表征骨内部结构特性。c、γ为常数。
步骤五：对骨结构有限元网格模型加载力学边界条件和载荷；
步骤六：对骨结构进行有限元分析求解；
步骤七：提取骨结构有限元模型各单元的主应力；
步骤八：根据各单元的主应力判断各单元的应力状态；
步骤九：基于各单元的应力状态，根据力学原理选择各单元适用的力学强度理论；
式中s为力学激励，σ1、σ2、σ3为主应力，且σ1≥σ2≥σ3，μ为材料泊松比，[σt]为材料的抗拉许用应力，[σc]为材料的抗压许用应力。
步骤十：根据确定的力学强度理论应用对应的相当应力计算公式进行各单元力学激励的计算；
步骤十一：应用骨重建控制方程计算下一步的表观密度：
ρt＝ρt+δt+b(s-k)δt-------------------------------(4)
式(3)为微分方程，式(4)为数值计算迭代公式，式中ρ为骨的表观密度，单位体积骨的质量，用来表征骨内部结构特性。b为骨重建速率。s为力学激励，用来表征力学效应的幅度，取相当应力作为力学激励。k为骨重建平衡态的力学激励；ρt为t时刻骨的表观密度，ρt+δt为t+δt时刻骨的表观密度，δt为积分步长。
步骤十二：根据各单元的新的表观密度计算各单元新的弹性模量材料参数：
步骤十三：对各单元重新定义材料参数
步骤十四：骨量收敛判断，骨量趋于稳定(骨量与迭代步数曲线趋于平稳或者骨量不再变化)，退出骨重建迭代循环，结束骨结构预测计算。否则转至步骤六开始下一步的迭代计算；
步骤十五：进入有限元分析后处理，进行各单元骨的表观密度显示，即得到预测的骨结构。
本发明与现有技术相比有益效果在于：
1、把骨重建控制方程与力学原理深度结合，既能体现骨功能适应性的原理，又体现出在每一瞬时把骨看作是工程材料，使骨结构预测方法在理论上更为合理。
2、取相当应力作为骨重建的力学激励，根据适用的强度理论选择计算公式，而不是在整个骨结构预测迭代计算过程中力学激励计算公式保持不变，使骨结构的预测结果更为准确。
附图说明
图1是方形骨片模型图。
图2是预测过程中结构质量与迭代步数的关系曲线图。
具体实施方式
下面将结合说明书附图，对本发明作进一步说明。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
本实施例提供一种考虑应力状态的骨结构预测方法，具体为：
步骤一、选取2mm×2mm、厚度0.1mm方形骨片作为预测骨结构的初始拓扑范围，骨片底部加载约束，顶部承受三角形分布载荷，方形骨片划分为40×40四节点四边形单元，见图1；
步骤二：定义骨结构的初始表观密度为0.8g/cm3,根据公式1计算弹性模量，取泊松比为0.3，定义各单元的材料参数；
步骤三：有限元分析求解；
步骤四：提取各单元的主应力，并判断应力状态，选取合适力学激励计算公式；
步骤五：根据公式4计算各单元新的表观密度；
步骤六：根据公式1计算各单元新的弹性模量，并对各单元进行定义材料；
步骤七：骨量收敛判断，骨量趋于稳定(见图2)，退出骨重建迭代循环，结束骨结构预测计算。否则转至步骤三开始下一步的迭代计算；
步骤八：进入有限元分析后处理，进行各单元骨的表观密度显示，即得到预测的骨结构。
本申请中的一种考虑应力状态的骨结构预测方法，在理论上更为合理，预测结果更为准确。
以上所述，仅为本发明创造较佳的具体实施方式，但本发明创造的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明创造披露的技术范围内，根据本发明创造的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明创造的保护范围之内。