CN108763818A - 一种适用于周期动载荷的骨结构预测方法 - Google Patents
一种适用于周期动载荷的骨结构预测方法 Download PDFInfo
- Publication number
- CN108763818A CN108763818A CN201810599983.3A CN201810599983A CN108763818A CN 108763818 A CN108763818 A CN 108763818A CN 201810599983 A CN201810599983 A CN 201810599983A CN 108763818 A CN108763818 A CN 108763818A
- Authority
- CN
- China
- Prior art keywords
- bone
- unit
- bone structure
- stress
- apparent density
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/06—Power analysis or power optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Prostheses (AREA)
Abstract
本发明公开了一种适用于周期动载荷的骨结构预测方法,包括:步骤一确定要进行预测的骨结构的拓扑范围;步骤二建模;步骤三有限元网格划分;步骤四计算骨的材料性能参数弹性模量,并对每个单元定义材料属性;步骤五:对骨结构有限元模型加载力学边界条件和载荷;步骤六:对骨结构进行有限元分析求解;步骤七:提取骨结构有限元模型各单元的主应力、米塞斯应力;步骤八:根据各单元的应力值,计算疲劳寿命循环次数;步骤九:计算各单元的力学激励;步骤十:应用骨重建控制方程计算下一步的表观密度:本申请使骨结构预测既能体现出骨作为一种有生命器官的骨功能适应性原理,又能体现出骨主要承受周期动载荷这一特点。
Description
技术领域
本发明涉及一种骨结构的预测方法,具体是一种适用于周期动载荷的骨结构预测方法。
背景技术
骨结构预测(也称为骨重建数值模拟)是研究骨重建理论的重要方法,为临床上预测骨密度变化或假体优化设计提供理论依据,也有助于各种骨代谢疾病的研究和新的骨代谢疾病治疗方案的提出。骨结构预测方程是“需要就生长,不需要就吸收”的骨功能适应性原理的数值实现算法。从工程观点看,是把骨骼作为一种工程材料的自适应优化设计。因此,骨功能适应性原理也称为骨结构自优化理论。
骨结构预测控制方程一般为:式中ρ为骨的表观密度,单位体积骨的质量,用来表征骨内部结构特性。B为骨重建速率。S为力学激励,用来表征力学效应的幅度,可以取为应变、应力、应变能密度等。K为骨重建平衡态的力学激励。雷周激欣等采用应变能密度、米塞斯等效应力、米塞斯等效应变3种不同力学激励对股骨近端骨结构进行了预测。(参见:雷周激欣,王冬梅,王春慧,等.不同力学激励对骨重建数值模拟的影响[J].医用生物力学,2015,30(4):299-303.)骨承受的载荷非常复杂,主要是日常活动以及运动引起的动载荷。正常情况下,人的运动等活动均是有规律的。因此,骨承受的主要载荷是规律性的周期动载荷。根据力学理论,材料在动载荷作用下主要失效形式是疲劳。当应力水平高于材料的持久疲劳极限时,每一水平的应力都有相应的疲劳寿命(循环次数),每一应力循环都将使材料受到损伤。因此,损伤是承受周期动载荷骨的力学效应的适当衡量形式。当前,骨结构预测方法存在的问题和缺点是力学激励不论采用应变能密度、米塞斯应力,还是米塞斯应变等衡量力学效应的表征量均未考虑承受周期动载荷骨结构疲劳形式的失效。骨疲劳损伤是临床上一种常见的骨损伤方式。因此,当前的骨结构预测算法从力学原理来看是存在一些缺陷的。
发明内容
由于现有的骨结构预测方法未充分考虑骨主要承受周期动载荷这一特点,本发明的目的是提供一个适用于周期动载荷的骨结构预测方法,使骨结构预测既能体现出骨作为一种有生命器官的骨功能适应性原理,又能体现出骨主要承受周期动载荷这一特点。
为实现上述目的,本申请的技术方案为:一种适用于周期动载荷的骨结构预测方法,包括:
步骤一:确定要进行预测的骨结构的几何拓扑范围,并根据力学原理抽象出力学边界条件和载荷;
步骤二:对步骤一确定的骨结构进行建模,建模可采用CAD软件或有限元分析软件完成;
步骤三:对步骤二得到的骨结构模型进行有限元网格划分,得到骨结构有限元网格模型;
步骤四:定义骨结构的初始表观密度,并根据弹性模量与表观密度的关系式计算骨的材料性能参数弹性模量,泊松比取为0.3,并对每个单元定义材料属性;
E=Cργ-------------------(1)
式中E为弹性模量,ρ为骨的表观密度,单位体积骨的质量,用来表征骨内部结构特性。C、γ为常数。
步骤五:对骨结构有限元模型加载力学边界条件和载荷;
步骤六:对骨结构进行有限元分析求解;
步骤七:提取骨结构有限元模型各单元的主应力、米塞斯应力等各种应力,包括但不限于;
步骤八:根据各单元的应力值,应用下面的公式计算疲劳寿命循环次数;
log Nσ=H logσ+JT+Kρ+M-------------------(2)
式中,Nσ为疲劳寿命循环次数;σ为应力(单位:MPa);T为温度,取为37℃;ρ为表观密度(g/cm3);H,J,K,M为经验常数。
步骤九:各单元力学激励的计算;
以积分步长Δt为一个应力谱周期,有σ1,σ2,…σk等k级应力,N1,N2,...NK为对应应力的疲劳寿命循环次数,n1,n2,...nK为对应应力的循环次数,S为力学激励,取一个应力谱周期的损伤。
步骤十:应用骨重建控制方程计算下一步的表观密度:
ρt+Δt=ρt+B(S-K)Δt---------------------(5)
式(4)为微分方程,式(5)为数值计算迭代公式,式中ρ为骨的表观密度,单位体积骨的质量,用来表征骨内部结构特性。B为骨重建速率。S为力学激励,用来表征力学效应的幅度,取损伤作为力学激励。K为骨重建平衡态的力学激励。ρt为t时刻骨的表观密度,ρt+Δt为t+Δt时刻骨的表观密度,Δt为积分步长。
步骤十一:根据各单元的新的表观密度计算各单元新的材料性能参数弹性模量:
步骤十二:对各单元重新定义材料参数;
步骤十三:骨量收敛判断,骨量趋于稳定,退出骨重建迭代循环,结束骨结构预测计算。否则转至步骤六开始下一步的迭代计算;
步骤十四:进入有限元分析后处理,进行各单元骨的表观密度显示,即得到预测的骨结构。
本发明由于采用以上技术方案,能够取得如下的技术效果:
1、把骨重建控制方程与力学积累损伤理论深度结合,既能体现骨功能适应性的原理,又体现出骨主要承受周期动载荷这一特点,使骨结构预测方法在理论上更为合理。
2、取一个应力谱周期的损伤作为骨重建的力学激励,根据积累损伤公式进行计算,把力学积累损伤理论引入骨生物力学。
本方法在理论上更为合理,体现了疲劳损伤是临床上一种常见的骨损伤方式,使预测结果更为准确。
附图说明
图1为方板模型图;
图2为预测过程中结构质量与迭代步数的关系曲线图。
具体实施方式
为了使本发明的目的、技术方案和优点更加清楚,下面结合附图和具体实施例对本发明进行详细描述。
实施例1
本实施例提供一种适用于周期动载荷的骨结构预测方法,具体包括如下步骤:
步骤一、选取2mm×2mm、厚度0.1mm方形骨片作为预测骨结构的初始拓扑范围,骨片底部一端施加固定铰支约束,一端施加活动铰支约束,以一天作为载荷谱周期,确定载荷大小等级数以及每级载荷的循环次数,方形骨片四边施加载荷,方形骨片划分为40×40四节点四边形单元,见图1;
步骤二:定义骨结构的初始表观密度为0.8g/cm3,根据公式1计算弹性模量,取泊松比为0.3,定义各单元的材料参数;
步骤三:有限元分析求解;
步骤四:提取各单元的不同大小载荷级别的米塞斯应力;
步骤五:根据公式2计算各大小等级应力的疲劳寿命循环次数;
步骤六:根据公式3计算各单元的力学激励;
步骤七:根据公式5计算各单元新的表观密度;
步骤八:根据公式1计算各单元新的弹性模量,并对各单元进行定义材料;
步骤九:骨量收敛判断,骨量趋于稳定(见图2),退出骨重建迭代循环,结束骨结构预测计算。否则转至步骤三开始下一步的迭代计算;
步骤十:进入有限元分析后处理,进行各单元骨的表观密度显示,即得到预测的骨结构。
以上所述,仅为本发明创造较佳的具体实施方式,但本发明创造的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明创造披露的技术范围内,根据本发明创造的技术方案及其发明构思加以等同替换或改变,都应涵盖在本发明创造的保护范围之内。
Claims (5)
1.一种适用于周期动载荷的骨结构预测方法,其特征在于,包括:
步骤一:确定要进行预测的骨结构的几何拓扑范围,并根据力学原理抽象出力学边界条件和载荷;
步骤二:对步骤一确定的骨结构进行建模;
步骤三:对步骤二得到的模型进行有限元网格划分,得到骨结构有限元网格模型;
步骤四:定义骨结构的初始表观密度,计算骨的材料性能参数弹性模量,并对每个单元定义材料属性;
步骤五:对骨结构有限元模型加载力学边界条件和载荷;
步骤六:对骨结构进行有限元分析求解;
步骤七:提取骨结构有限元模型各单元的主应力、米塞斯应力;
步骤八:根据各单元的应力值,计算疲劳寿命循环次数;
步骤九:计算一个应力谱周期各单元的损伤,即为力学激励;
步骤十:应用骨重建控制方程计算下一步的表观密度:
步骤十一:根据各单元新的表观密度计算得到各单元新的材料性能参数弹性模量:
步骤十二:对各单元重新定义材料参数;
步骤十三:骨量收敛判断,骨量趋于稳定,退出骨重建迭代循环,结束骨结构预测计算;否则转至步骤六开始下一步的迭代计算;
步骤十四:进入有限元分析后处理,进行各单元骨的表观密度显示,即得到预测的骨结构。
2.根据权利要求1所述一种适用于周期动载荷的骨结构预测方法,其特征在于,步骤四中计算骨的材料性能参数弹性模量,具体为;
E=Cργ-------------------(1)
式中E为弹性模量,ρ为骨的表观密度,单位体积骨的质量,用来表征骨内部结构特性,C、γ为常数。
3.根据权利要求1所述一种适用于周期动载荷的骨结构预测方法,其特征在于,步骤八中计算疲劳寿命循环次数,具体为:
logNσ=Hlogσ+JT+Kρ+M------------------(2)
式中,Nσ为疲劳寿命循环次数;σ为应力,T为温度;ρ为表观密度;H,J,K,M为经验常数。
4.根据权利要求1所述一种适用于周期动载荷的骨结构预测方法,其特征在于,步骤九中计算各单元的力学激励,具体为;
以积分步长Δt为一个应力谱周期,有σ1,σ2,…σk等k级应力,N1,N2,...NK为对应应力的疲劳寿命循环次数,n1,n2,...nK为对应应力的循环次数,S为力学激励,取一个应力谱周期的损伤。
5.根据权利要求1所述一种适用于周期动载荷的骨结构预测方法,其特征在于,步骤十中应用骨重建控制方程计算下一步的表观密度,具体为:
ρt+Δt=ρt+B(S-K)Δt---------------------(5)
式(4)为微分方程,式(5)为数值计算迭代公式,式中ρ为骨的表观密度,单位体积骨的质量,用来表征骨内部结构特性;B为骨重建速率;S为力学激励,用来表征力学效应的幅度,取损伤作为力学激励;K为骨重建平衡态的力学激励;ρt为t时刻骨的表观密度,ρt+Δt为t+Δt时刻骨的表观密度,Δt为积分步长。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810599983.3A CN108763818A (zh) | 2018-06-12 | 2018-06-12 | 一种适用于周期动载荷的骨结构预测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810599983.3A CN108763818A (zh) | 2018-06-12 | 2018-06-12 | 一种适用于周期动载荷的骨结构预测方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN108763818A true CN108763818A (zh) | 2018-11-06 |
Family
ID=64022066
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810599983.3A Pending CN108763818A (zh) | 2018-06-12 | 2018-06-12 | 一种适用于周期动载荷的骨结构预测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108763818A (zh) |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103313682A (zh) * | 2010-07-08 | 2013-09-18 | 因特脉管有限公司 | 用于对多个管腔内外科卡钉进行放置的部署装置 |
CN104840265A (zh) * | 2015-06-11 | 2015-08-19 | 大连大学 | 能够即可负载的一体式种植体 |
-
2018
- 2018-06-12 CN CN201810599983.3A patent/CN108763818A/zh active Pending
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103313682A (zh) * | 2010-07-08 | 2013-09-18 | 因特脉管有限公司 | 用于对多个管腔内外科卡钉进行放置的部署装置 |
CN104840265A (zh) * | 2015-06-11 | 2015-08-19 | 大连大学 | 能够即可负载的一体式种植体 |
Non-Patent Citations (5)
Title |
---|
李昊等: "皮质骨中骨单元应力集中效应的有限元分析", 《医用生物力学》 * |
马宗民等: "动载荷作用下骨重建力学调控机制", 《医用生物力学》 * |
马宗民等: "老龄引起骨质疏松的模拟", 《生物医学工程研究》 * |
马宗民等: "考虑载荷特性的骨重建力学调控机制", 《力学与实践》 * |
马超等: "基于自组织控制过程的骨重建过程的有限元仿真", 《航天医学与医学工程》 * |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Genet et al. | Modeling pathologies of diastolic and systolic heart failure | |
Jaatinen et al. | Extended phase diagram of the three-dimensional phase field crystal model | |
Prakash et al. | Simulation of micromechanical behavior of polycrystals: finite elements versus fast Fourier transforms | |
Benedetti et al. | A grain boundary formulation for crystal plasticity | |
Xiao et al. | Modeling and simulation of curled dry leaves | |
WO2023078380A1 (zh) | 基于sn曲线的程序载荷谱编制方法、系统和存储介质 | |
Anderson et al. | Virtual experiments, physical validation: dental morphology at the intersection of experiment and theory | |
Rocha et al. | A macro finite-element formulation for cardiac electrophysiology simulations using hybrid unstructured grids | |
CN109002630A (zh) | 一种超弹性材料的快速仿真方法 | |
Gong et al. | A nonlinear elimination preconditioned inexact Newton method for heterogeneous hyperelasticity | |
CN108763818A (zh) | 一种适用于周期动载荷的骨结构预测方法 | |
Yan et al. | FEM modeling method of damage structures for structural damage detection | |
CN108595885A (zh) | 一种考虑应力状态的骨结构预测方法 | |
CN105843984B (zh) | 一种基于h细化的电阻层析成像有限元建模方法 | |
Rohan et al. | Finite element modelling of nearly incompressible materials and volumetric locking: a case study | |
Thurieau et al. | A simple solution procedure to 3D-piezoelectric problems: Isotropic BEM coupled with a point collocation method | |
CN108647466A (zh) | 一种适用于动静耦合载荷的骨结构预测方法 | |
Viswanath et al. | Deep learning for topology optimization of triply periodic minimal surface based gyroid-like structures | |
Washio et al. | Analysis of spontaneous oscillations for a three-state power-stroke model | |
Li et al. | Development of an efficient wetting and drying treatment for shallow‐water modeling using the quadrature‐free Runge‐Kutta discontinuous Galerkin method | |
Álamo et al. | Adapting the simp model for topology optimization of biomechanical structures | |
Saito | Conservative numerical schemes for the Keller-Segel system and numerical results (Mathematical analysis on the self-organization and self-similarity) | |
Qu et al. | Damage model of bone under mechanical and electromagnetic loadings | |
Duong et al. | Modeling and simulation of a growing mass by the smoothed finite element method (SFEM) | |
Gao et al. | Backreaction and order reduction in initially contracting models of the Universe |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
RJ01 | Rejection of invention patent application after publication |
Application publication date: 20181106 |
|
RJ01 | Rejection of invention patent application after publication |