CN116525118B - 平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法 - Google Patents

平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法 Download PDF

Info

Publication number
CN116525118B
CN116525118B CN202310283842.1A CN202310283842A CN116525118B CN 116525118 B CN116525118 B CN 116525118B CN 202310283842 A CN202310283842 A CN 202310283842A CN 116525118 B CN116525118 B CN 116525118B
Authority
CN
China
Prior art keywords
lung
respiratory
model
pressure
stage
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.)
Active
Application number
CN202310283842.1A
Other languages
English (en)
Other versions
CN116525118A (zh
Inventor
鲍玉冬
王思远
李旭
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin University of Science and Technology
Original Assignee
Harbin University of Science and Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Harbin University of Science and Technology filed Critical Harbin University of Science and Technology
Priority to CN202310283842.1A priority Critical patent/CN116525118B/zh
Publication of CN116525118A publication Critical patent/CN116525118A/zh
Application granted granted Critical
Publication of CN116525118B publication Critical patent/CN116525118B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/50—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for simulation or modelling of medical disorders
    • 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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • 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/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Medical Informatics (AREA)
  • Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Pathology (AREA)
  • Algebra (AREA)
  • General Health & Medical Sciences (AREA)
  • Epidemiology (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Biomedical Technology (AREA)
  • Primary Health Care (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)

Abstract

平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法，属于人体呼吸运动数值模拟技术领域。为解决在人体呼吸运动下的人体呼吸运动数值模拟技术，本发明根据肺部静力学特性将呼气阶段分为吸气阶段、呼气阶段，在膈肌和肋间肌上建立吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型，通过牛顿第二定律构建吸气阶段、呼气阶段的呼吸肌作用力的力学模型，进行计算整合和可靠性验证，建立呼吸系统有限元模型，进行呼吸系统动力学仿真计算分析得到呼吸系统动力学仿真计算分析结果，分别对气管未变形和最大变形时的有限元模型进行流体力学仿真计算，基于第三方虚拟平台进行平静呼吸下人体呼吸运动数值模拟。本发明用于人体呼吸运动数值模拟技术领域。

Description

平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法
技术领域
本发明属于人体呼吸运动数值模拟技术领域，具体涉及平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法。
背景技术
目前，支气管镜术依然是呼吸系统疾病治疗的主要依赖手段，借助不同匹配的介入性仪器，通过支气管镜在气管内完成药物注射、肺灌洗、针穿刺、异物钳取、支架置取等多种功能，最终达到治疗病灶的目的。然而，支气管镜作为典型的人工主导的介入性医疗器械，依然面临着人工介入带来的术中不适、术后并发症等主要问题。因此支气管镜术的发展和研究不断趋向自动化和智能化，以解决人工介入问题对支气管镜术治疗过程中带来的困扰。但这类研究由于并未考虑气管作为介入环境实际存在呼吸运动的情况，成果并不能很好的进行实践推广，亟需一种对人体呼吸运动进行数值模拟的方法对该类技术研究进行推进。
发明内容
本发明是为了解决支气管镜术自动化、智能化发展中的关键技术研究因忽略人体呼吸运动情况而不能更好的完成相关实验和理论实践应用的问题，提出涉及平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法。
为实现上述目的，本发明通过以下技术方案实现：
一种平静呼吸下人体呼吸运动数值模拟方法，包括如下步骤：
S1、根据肺部静力学特性将呼气阶段分为吸气阶段、呼气阶段，在膈肌和肋间肌上建立吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型，然后通过牛顿第二定律构建吸气阶段、呼气阶段的呼吸肌作用力的力学模型，吸气阶段呼吸肌作用力的力学模型为：
其中，Pimf为肋间肌作用力，Rup为上呼吸道产生的气阻，CL为肺顺应性，Cth为胸廓顺应性，PT为跨肺压，q为体积流量，Valv为肺容积，Vc为胸腔容积，PRfbs为无软骨气道组织的压强；
Pdf为膈肌作用力，Cab为腹部顺应性，μ1、μ2、Ap、Bp、Ae、Be、Ke、η1、η2均为常数；
VRmax取值为0.15，VRfbs为余下部分支气管系统容积，Ppl为胸膜腔内压，TLC为肺总容量，RV为余气量，K为修正系数，Pw为以压强形式表达的壁面液面表面张力；
呼气阶段呼吸肌作用力的力学模型为：
其中，P’imf为肋间肌名义作用力，P’df为膈肌名义作用力，RRfbs为无软骨气道组织的气阻，μ3为常数；
S2、将步骤S1得到的吸气阶段、呼气阶段的呼吸肌作用力的力学模型进行计算整合和可靠性验证；
S3、建立呼吸系统有限元模型，包括肺、气管、膈肌和肋间肌，以步骤S2获得的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果作为有限元仿真计算的载荷条件，分别施加在呼吸系统有限元模型的膈肌和肋间肌上，进行呼吸系统动力学仿真计算分析，得到呼吸系统动力学仿真计算分析结果；
S4、从步骤S3的呼吸系统动力学仿真计算分析结果中得到气管未变形的有限元模型、气管最大变形时的有限元模型，分别对气管未变形和最大变形时的有限元模型进行流体力学仿真计算；
S5、基于第三方虚拟平台Unity 3D，根据步骤S4得到的气管未变形和最大变形时的有限元模型进行流体力学仿真计算结果，将步骤S3得的呼吸系统动力学仿真计算分析结果与Unity 3D进行数据关联，在Unity 3D中进行平静呼吸下人体呼吸运动数值模拟。
进一步的，步骤S1的具体实现方法包括如下步骤：
S1.1、根据肺部静力学构建跨肺压模型，其中跨肺压为肺组织所受的压力差，其具体实现方法包括如下步骤：
S1.1.1、在跨肺压的作用下，胸膜壁与肺内相互联结的肺泡隔膜壁所受应力之和为0，所以肺组织的径向力平衡公式为：
其中，Pi、Po分别为壁内压强、壁外压强，ΣFo/A、ΣFi/A分别为壁面抗衡壁内压强的内力、壁面抗衡壁外压强的内力；
则对胸膜壁和肺泡壁的力平衡关系，可以得到以下公式：
其中，Ppi、Ppo分别为胸膜壁内压和胸膜壁外压，Σσpi为胸膜壁回缩应力，Σσpo为胸膜壁扩张应力，Pai、Pao分别为肺泡隔膜壁内压和肺泡隔膜壁外压，Σσai为局部回缩应力，Σσao为局部扩张应力；
S1.1.2、根据Dalton气体分压定律和帕斯卡原理，在能量损失忽略不计的情况下，认为跨肺压在肺内的传递不失真，跨肺压表示为以下公式：
PT＝(Ppi-Ppo)+∑σI-∑σO+Pw
其中，ΣσI＝ΣΣσai为肺总体回缩张力，ΣσO＝ΣΣσao为肺总体扩张张力；
肺总体回缩张力和肺总体扩张张力主要是由肺泡隔膜的弹性应力及外力作用引起的，肺总体回缩张力代表了肺组织的弹性回缩力，肺总体扩张张力代表了肺泡壁隔膜两侧的压力差，在肺部受到呼吸肌牵引而变形引发其内部应力发生改变的时候，肺本身所受内力的作用主要可以分为表面张力和弹性应力，弹性应力表示为：
其中，Nn为肺泡隔膜微元的法向张力，R为肺泡初始半径；
法向张力使用肺的拟应变能函数表达：
其中，c，a1，a2，a4为常数；E1，E2为主应变；ρo为参考状态下的肺密度，ρoW(e)为参考状态每单位肺体积的应变能；
以二维应力状态表示肺泡隔膜的应力-应变关系，肺泡隔膜微元的法向张力为拟应变能函数对应变的导数：
其中，δ为肺泡隔膜的厚度，根据各向同性得到E1＝E2；
S1.1.3、肺部组织为连续介质，依据肺应力-应变关系建立肺组织变形运动的场方程为：
其中，xj为三维笛卡尔坐标x1、x2、x3，v为肺某点的变形速度，ρ为肺组织的密度，为肺某点的位置方程，U为肺部组织的位移场函数，肺部为连续介质，肺组织变形运动的场方程满足应变协调方程；
S1.1.4、利用肺组织和胸壁的弹性回缩力表达肺组织的弹性回缩力Fe为：
胸壁的弹性回缩力Fp为：
最终得放到跨肺压：
S1.2、进行呼吸系统吸气和呼气阶段的作用力分析，然后在膈肌和肋间肌上建立吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型，具体实现方法包括如下步骤：
S1.2.1、整个呼吸系统呼气过程主要受到五个力的作用：肋间肌作用力、膈肌作用力、跨肺压PT、跨膈压Pdi和气阻R，则得到以下计算公式：
PT＝Palv-Ppl
其中，Palv为肺内压；Ppl为胸膜腔内压；
其中，C为顺应性，ΔV为容积变化量，ΔP为压力变化量；
S1.2.2、跨膈压和腹部对于呼吸系统的影响使用腹部顺应性Cab进行替代，平静呼吸时取值为1L/cmH2O，
气阻是气道阻力，属于空气在气管内流动而产生的力，由于软骨的支撑作用，气管的变形量不明显，则上呼吸道产生的气阻由下式表示：
Rup＝μ1+μ2|q|
余下支气管系统没有软骨支撑，但不包含2mm管径以下的部分，余下部分支气管系统气阻RRfbs由下式确定：
S1.2.3、吸气过程中，只考虑上呼吸道气阻，呼气过程主要受三个力的作用，跨肺压、跨膈压和上呼吸道及余下支气管系统部分的气阻，呼吸过程中还要考虑没有软骨包裹的余下支气管系统的弹性，则得的公式为：
S1.2.4、以肺部为基础，分别就肋间肌和膈肌建立吸气和呼气阶段的弹簧阻尼等效模型，肺部的弹簧阻尼模型以顺应性和跨肺压为主要部分，肺组织的粘滞特性代表了阻尼RL，它是肺的拟应变能函数的部分表示，另一部分表示为作用在肺上的弹性回缩力Fe和Fp，Pw和Ppl是跨肺压关系式中的剩余元素，肺部顺应性CL代表肺组织的弹性，作用在肺上的还有余下支气管系统的弹性PRfbs，作用在肋间肌和膈肌之上的力为自身作用力以及上呼吸道的气道阻力，胸廓顺应性Cth表达了胸廓本身的弹性行为，腹部顺应性Cab代表了腹部和跨膈压对于呼吸系统的影响，呼气过程相对于吸气过程而言，呼气过程中膈肌和肋间肌不考虑自身作用力且多考虑了余下支气管系统气阻RRfbs；
S1.3、通过牛顿第二定律，依据吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型建立呼吸肌作用力的力学模型，力学模型确定的输入变量为肺容积Valv、胸腔容积Vc、余下支气管系统容积VRfbs以及体积流量q。
进一步的，步骤S2的具体实现方法包括如下步骤：
S2.1、根据步骤S1得到的吸气阶段、呼气阶段的呼吸肌作用力的力学模型的输入变量，利用Origin从文献曲线获取模型求解的初始数据；
S2.2、根据步骤S2.1得到的初始数据使用MATLAB求解步骤S1得到的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果，求解时按照整个周期的呼吸过程进行阶段计算整合；
S2.3、使用MATLAB Cftool拟合步骤S2.2获得的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果，平滑作用力曲线；
S2.4、对比步骤S2.3拟合后的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果和文献曲线，验证吸气阶段、呼气阶段的呼吸肌作用力的力学模型的可靠性。
进一步的，步骤S2中的文献曲线从《人体呼吸系统力学属性下呼吸动力过程的仿真研究》中的图2获得。
进一步的，步骤S3的具体实现方法包括如下步骤：
S3.1、自开源数据库https://www.cancerimagingarchive.net获得CT数据，使用Mimics建立呼吸系统有限元模型，包括肺、气管、膈肌和肋间肌；
S3.2、采用Hypermesh完成有限元模型处理和网格划分，单元类型为四面体单元，气管部分为三角形面网格单元，导入到ABAQUS中赋予组织的材料参数；
S3.3、以步骤S2获得的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果作为有限元仿真计算的载荷条件，将肺门处及肺部里侧进行除z轴位移外五自由度约束，对膈肌的中央部分、上呼吸道部分及肺尖处进行完全约束，进行呼吸系统动力学仿真计算分析，得到呼吸系统动力学仿真计算分析结果。
进一步的，步骤S4的具体实现方法包括如下步骤：
S4.1、通过ABAQUS从呼吸系统动力学仿真计算分析结果中获取气管未变形和最大变形时的有限元模型，然后经Geomagic Wrap、Geomagic Design和SolidWorks实体化处理，导入到ANSYS中进行网格划分；
S4.2、在Fluent中以入口速度曲线的吸气峰值和呼气峰值为边界条件，进行2次气管未变形的有限元模型、气管最大变形时的有限元模型的流体力学仿真计算，对比2次流体力学仿真计算结果，得到平静呼吸时，流体对于呼吸运动中气管变形的影响忽略不计的结论。
进一步的，步骤S4.1中经气管未变形时及最大变形时模型在吸气峰4.02m/s和呼气峰-4.59m/s入口速度下的流体力学仿真计算分析。
进一步的，步骤S5的具体实现方法包括如下步骤：
S5.1、步骤S3得的呼吸系统动力学仿真计算分析结果与Unity 3D进行数据关联的具体实现方法包括如下步骤：
S5.1.1、将步骤S3得到的气管未变形的有限元模型导入到3DMax中按照气管级别进行分割，分割后的部分导入到Unity 3D中；
S5.1.2、针对分割后的每个部分，用Mimics测量气管未变形时自定义点到分割面处的距离，用ABAQUS查询气管变形时自定义点处的位移量；
S5.1.3、根据分割后的每个部分得到的气管未变形时自定义点到分割面处的距离及气管变形时自定义点处的位移量，计算分割后的每个部分的缩放系数，并定义缩放系数的变动范围、变动方式以及变动周期，缩放系数计算公式为：
其中，ζ为缩放系数，S1为变形后节点到局部坐标系原点的距离，S0为变形前节点到局部坐标系原点的距离，Δx为节点变形前后的位移。
S5.2、在Unity 3D中进行平静呼吸下人体呼吸运动数值模拟的方法包括如下步骤：
S5.2.1、设置脚本1获取分割面中心坐标，使用public函数公开功能选项板，通过鼠标点选的方式选取分割部分的分割面中心坐标；
S5.2.1、设置脚本2读取脚本1获得的中心坐标后，使用public函数公开分割部分的最大缩放系数和呼吸运动周期文本输入功能选项板，输入分割部分的最大缩放系数和周期，利用SetEase缓动函数实现线性缩放，SetLoops(-1,LoopType.Yoyo)实现无限循环。
一种平静呼吸下人体呼吸运动数值模拟系统，依托于所述的一种平静呼吸下人体呼吸运动数值模拟方法实现，包括力学模型计算模块、力学模型拟合模块、呼吸系统动力学仿真计算模块，所述力学模型计算模块、力学模型拟合模块、呼吸系统动力学仿真计算模块之间相互连接，所述呼吸系统动力学仿真计算模块具有Unity 3D关联启动功能。
进一步的，所述力学模型计算模块的输入变量与患者的肺功能检查报告中的潮气量产生关联，关联方式具体为在平静呼吸状态下，Valv等同于潮气量与功能残气量的和，Vc、VRfbs与Valv的关系为Vc＝Valv+VRfbs+100ml，q、VRfbs取值于文献数据，实现了平静呼吸下患者的个性化呼吸运动数值模拟的功能。
本发明的有益效果：
本发明所述的一种平静呼吸下人体呼吸运动数值模拟系统，提供了一个有效的呼吸肌作用力的力学模型，为患者呼吸力学个性化发展奠定了一定的基础，用于计算呼吸肌作用力，该作用力可作为边界条件施加到有限元仿真计算中，进而完成人体呼吸运动的数值模拟。同时，搭建了一个人体呼吸运动数值模拟系统，可完成从力学模型计算到呼吸运动数值模拟的全部过程。
本发明所述的一种平静呼吸下人体呼吸运动数值模拟方法，该数值模拟方法可用于建立数字孪生技术和虚拟手术研究中人体呼吸系统平静呼吸下的虚拟可视化模型，为呼吸系统虚拟模型的发展提供了一定的基础。该方法能够行之有效的实现患者个性化的平静呼吸下呼吸运动的数值模拟，可用于解决数字孪生及虚拟手术中对呼吸系统虚拟模型的需求，且该方法具有可开发性，能够根据不同的目的进行开发，完成相应不同类别的研究，同时可更好的推进当前支气管镜术自动化、智能化研究的进展和成果的实践应用
附图说明
图1为本发明所述的一种平静呼吸下人体呼吸运动数值模拟方法的流程图；
图2为本发明平静呼吸下吸气阶段呼吸系统动力学分析结果；
图3为本发明平静呼吸下呼气阶段呼吸系统动力学分析结果；
图4为本发明平静呼吸下吸气阶段弹簧阻尼等效模型；
图5为本发明平静呼吸下呼气阶段弹簧阻尼等效模型；
图6为本发明吸气阶段、呼气阶段的呼吸肌作用力的力学模型整合计算及验证流程图；
图7为本发明吸气阶段、呼气阶段的呼吸肌作用力的力学模型验证对比结果；
图8为本发明呼吸系统有限元模型建立结果；
图9为本发明呼吸系统动力学仿真计算位移云图；
图10为本发明呼吸系统动力学仿真计算气管节点位移数据对比；
图11为本发明呼吸系统动力学仿真计算肺部节点位移数据对比；
图12为本发明气管流体力学仿真计算的入口速度；
图13为本发明气管未变形和最大变形时的有限元模型吸气峰入口速度流体力学仿真计算结果图；
图14为本发明气管未变形和最大变形时的有限元模型呼气峰入口速度流体力学仿真计算结果图；
图15为本发明动力学仿真计算结果与系统的数据关联方法示意图；
图16为本发明气管未变形和最大变形时的有限元模型的分割结果；
图17为本发明气管的呼吸运动模拟结果；
图18为本发明一种平静呼吸下人体呼吸运动数值模拟系统的界面；
图19为本发明一种平静呼吸下人体呼吸运动数值模拟系统的工作示意图；
图20为本发明一种平静呼吸下人体呼吸运动数值模拟系统工作功能图；
图21为本发明针对患者在所述的一种平静呼吸下人体呼吸运动数值模拟系统中动力学仿真计算结果；
图22为本发明针对患者在所述的一种平静呼吸下人体呼吸运动数值模拟系统中呼吸运动模拟结果。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及具体实施方式，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅用以解释本发明，并不用于限定本发明，即所描述的具体实施方式仅仅是本发明一部分实施方式，而不是全部的具体实施方式。通常在此处附图中描述和展示的本发明具体实施方式的组件可以以各种不同的配置来布置和设计，本发明还可以具有其他实施方式。
因此，以下对在附图中提供的本发明的具体实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定具体实施方式。基于本发明的具体实施方式，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他具体实施方式，都属于本发明保护的范围。
为能进一步了解本发明的发明内容、特点及功效，兹例举以下具体实施方式，并配合附图1-附图22详细说明如下：
具体实施方式一：
一种平静呼吸下人体呼吸运动数值模拟方法，包括如下步骤：
S1、根据肺部静力学特性将呼气阶段分为吸气阶段、呼气阶段，在膈肌和肋间肌上建立吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型，然后通过牛顿第二定律构建吸气阶段、呼气阶段的呼吸肌作用力的力学模型；
进一步的，步骤S1的具体实现方法包括如下步骤：
S1.1、根据肺部静力学构建跨肺压模型，其中跨肺压为肺组织所受的压力差，其具体实现方法包括如下步骤：
呼吸系统动力学特性的分析以肺部静力学特性为基础，表现为跨肺压PT和顺应性CL，顺应性可取值为2L/cm H2O，跨肺压可推导为：肺可认为由许多肺泡隔膜构成的小单元集合体，在正常情况下，仅考虑肺组织的均匀性。跨肺压由此可演变为肺组织所受的压力差，肺内结构多数为径向对称；
S1.1.1、在跨肺压的作用下，胸膜壁与肺内相互联结的肺泡隔膜壁所受应力之和为0，所以肺组织的径向力平衡公式为：
其中，Pi、Po分别为壁内压强、壁外压强，ΣFo/A、ΣFi/A分别为壁面抗衡壁内压强的内力、壁面抗衡壁外压强的内力；
则对胸膜壁和肺泡壁的力平衡关系，可以得到以下公式：
其中，Ppi、Ppo分别为胸膜壁内压和胸膜壁外压，Σσpi为胸膜壁回缩应力，Σσpo为胸膜壁扩张应力，Pai、Pao分别为肺泡隔膜壁内压和肺泡隔膜壁外压，Σσai为局部回缩应力，Σσao为局部扩张应力；
S1.1.2、根据Dalton气体分压定律和帕斯卡原理，在能量损失忽略不计的情况下，认为跨肺压在肺内的传递不失真，跨肺压表示为以下公式：
PT＝(Ppi-Ppo)+∑σI-∑σO+Pw
其中，ΣσI＝ΣΣσai为肺总体回缩张力，ΣσO＝ΣΣσao为肺总体扩张张力；
肺总体回缩张力和肺总体扩张张力主要是由肺泡隔膜的弹性应力及外力作用引起的，肺总体回缩张力代表了肺组织的弹性回缩力，肺总体扩张张力代表了肺泡壁隔膜两侧的压力差，在肺部受到呼吸肌牵引而变形引发其内部应力发生改变的时候，肺本身所受内力的作用主要可以分为表面张力和弹性应力，弹性应力表示为：
其中，Nn为肺泡隔膜微元的法向张力，R为肺泡初始半径；
法向张力使用肺的拟应变能函数表达：
其中，c，a1，a2，a4为常数；E1，E2为主应变；ρo为参考状态下的肺密度，ρoW(e)为参考状态每单位肺体积的应变能；
以二维应力状态表示肺泡隔膜的应力-应变关系，肺泡隔膜微元的法向张力为拟应变能函数对应变的导数：
其中，δ为肺泡隔膜的厚度，根据各向同性得到E1＝E2；
S1.1.3、肺部组织为连续介质，依据肺应力-应变关系建立肺组织变形运动的场方程为：
其中，xj为三维笛卡尔坐标x1、x2、x3，v为肺某点的变形速度，ρ为肺组织的密度，为肺某点的位置方程，U为肺部组织的位移场函数，肺部为连续介质，肺组织变形运动的场方程满足应变协调方程；
进一步的，肺内弹性应力主要与肺内组织材料本身以及变形位移相关联，而肺部材料本身又以其位移场为输入。显然，肺部的位移场与肺容积的改变量是唯一相关的，肺组织的弹性应力与肺容积的变化量具有相关性，考虑以上边界条件的获取十分复杂，可利用肺组织和胸壁的弹性回缩力进行近似表达；
S1.1.4、利用肺组织和胸壁的弹性回缩力表达肺组织的弹性回缩力Fe为：
胸壁的弹性回缩力Fp为：
最终得放到跨肺压：
进一步的，弹簧阻尼等效模型的建立过程为先进行呼吸系统吸气和呼气阶段的作用力分析，然后搭建阶段内弹簧阻尼等效模型，具体为：如图2所示，整个呼吸系统呼气过程主要受到五个力的作用：肋间肌作用力、膈肌作用力、跨肺压PT、跨膈压Pdi和气阻R。跨肺压和顺应性C是描述肺组织特性的重要参量。跨肺压表示为肺内压和胸膜腔内压的压差，顺应性表示为组织变形难易的程度，即粘滞特性；
S1.2、进行呼吸系统吸气和呼气阶段的作用力分析，然后在膈肌和肋间肌上建立吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型，具体实现方法包括如下步骤：
S1.2.1、整个呼吸系统呼气过程主要受到五个力的作用：肋间肌作用力、膈肌作用力、跨肺压PT、跨膈压Pdi和气阻R，则得到以下计算公式：
PT＝Palv-Ppl
其中，Palv为肺内压；Ppl为胸膜腔内压；
其中，C为顺应性，ΔV为容积变化量，ΔP为压力变化量；
进一步的，跨膈压是腹部压力和肺内压的压差，膈肌是具有密封性质的器官，当吸气时膈肌下降，会促使腹部容积减少，腹部压力增加，从而对膈肌的作用产生影响。本研究中为动力学模型的建立，跨膈压和腹部对于呼吸系统的影响使用腹部顺应性Cab进行近似替代，平静呼吸时取值为1L/cmH2O。
气阻是气道阻力，属于空气在气管内流动而产生的力。气道阻力主要可分为三个部分：上呼吸道、余下支气管系统以及细管径支气管部分和末梢。不考虑气体交换和弥散，也不考虑气流进入呼吸道的变化所带来的能量耗散和湍流动能。气道阻力的大部分占比都在上呼吸道，余下的气管支气管体系所造成的阻力不大，尤其是细管径支气管部分和末梢阻力占比更少，可忽略不计。
在上呼吸道，由于鼻咽以及声门的作用，气流进入气道由层流变为湍流，不论是吸气还是呼气，气流在鼻咽和声门处变动最大，意味着阻力较大以及能量的耗散较多。针对于有软骨包络的气管部分，由于软骨的支撑作用，气管的变形量不明显，可认为没有变化；
S1.2.2、跨膈压和腹部对于呼吸系统的影响使用腹部顺应性Cab进行替代，平静呼吸时取值为1L/cmH2O，
气阻是气道阻力，属于空气在气管内流动而产生的力，由于软骨的支撑作用，气管的变形量不明显，则上呼吸道产生的气阻由下式表示：
Rup＝μ1+μ2|q|
余下支气管系统没有软骨支撑，但不包含2mm管径以下的部分，余下部分支气管系统气阻RRfbs由下式确定：
S1.2.3、吸气过程中，只考虑上呼吸道气阻，呼气过程主要受三个力的作用，跨肺压、跨膈压和上呼吸道及余下支气管系统部分的气阻，呼吸过程中还要考虑没有软骨包裹的余下支气管系统的弹性，则得的公式为：
/>
然后依据分析结果以肺部为基础，分别就肋间肌和膈肌建立吸气和呼气阶段的弹簧阻尼等效模型，如图4、5所示；
S1.2.4、以肺部为基础，分别就肋间肌和膈肌建立吸气和呼气阶段的弹簧阻尼等效模型，肺部的弹簧阻尼模型以顺应性和跨肺压为主要部分，肺组织的粘滞特性代表了阻尼RL，它是肺的拟应变能函数的部分表示，另一部分表示为作用在肺上的弹性回缩力Fe和Fp，Pw和Ppl是跨肺压关系式中的剩余元素，肺部顺应性CL代表肺组织的弹性，作用在肺上的还有余下支气管系统的弹性PRfbs，作用在肋间肌和膈肌之上的力为自身作用力以及上呼吸道的气道阻力，胸廓顺应性Cth表达了胸廓本身的弹性行为，腹部顺应性Cab代表了腹部和跨膈压对于呼吸系统的影响，呼气过程相对于吸气过程而言，呼气过程中膈肌和肋间肌不考虑自身作用力且多考虑了余下支气管系统气阻RRfbs；
S1.3、通过牛顿第二定律，依据吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型建立呼吸肌作用力的力学模型，力学模型确定的输入变量为肺容积Valv、胸腔容积Vc、余下支气管系统容积VRfbs以及体积流量q
吸气阶段呼吸肌作用力的力学模型为：
其中，Pimf为肋间肌作用力，Rup为上呼吸道产生的气阻，CL为肺顺应性，Cth为胸廓顺应性，PT为跨肺压，q为体积流量，Valv为肺容积，Vc为胸腔容积，PRfbs为无软骨气道组织的压强；
Pdf为膈肌作用力，Cab为腹部顺应性，μ1、μ2、Ap、Bp、Ae、Be、Ke、η1、η2均为常数；
VRmax取值为0.15，VRfbs为余下部分支气管系统容积，Ppl为胸膜腔内压，TLC为肺总容量，RV为余气量，K为修正系数，Pw为以压强形式表达的壁面液面表面张力；
呼气阶段呼吸肌作用力的力学模型为：
/>
其中，P’imf为肋间肌名义作用力，P’df为膈肌名义作用力，RRfbs为无软骨气道组织的气阻，μ3为常数；
进一步的，RRfbs为余下部分支气管系统的气阻(<2mm管径的支气管除外)，cmH2O·s/L；Cth为胸廓顺应性；VRmax取值为0.15。μ1、μ2、μ3、Ap、Bp、Ae、Be、Ke、η1、η2分别取值为0.34cmH2O·s/L、0.46cm H2O·s2/L2、0.2cm H2O·s·L、1.4cmH2O、-3.5cmH2O、0.2cmH2O、–0.5cmH2O、1L-1、2cm H2O、0.35cmH2O-1；Pw取值为15，修正系数K在呼气和吸气阶段不同，这里取呼气58.59，吸气53.11；TLC取值为5.19，RV取值为1.24；
S2、将步骤S1得到的吸气阶段、呼气阶段的呼吸肌作用力的力学模型进行计算整合和可靠性验证；
进一步的，步骤S2的具体实现方法包括如下步骤：
S2.1、根据步骤S1得到的吸气阶段、呼气阶段的呼吸肌作用力的力学模型的输入变量，利用Origin从文献曲线获取模型求解的初始数据；
S2.2、根据步骤S2.1得到的初始数据使用MATLAB求解步骤S1得到的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果，求解时按照整个周期的呼吸过程进行阶段计算整合；
S2.3、使用MATLAB Cftool拟合步骤S2.2获得的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果，平滑作用力曲线；
S2.4、对比步骤S2.3拟合后的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果和文献曲线，验证吸气阶段、呼气阶段的呼吸肌作用力的力学模型的可靠性。
步骤S2中的文献曲线从《人体呼吸系统力学属性下呼吸动力过程的仿真研究》中的图2获得，力学模型的验证结果如图7所示；
S3、建立呼吸系统有限元模型，包括肺、气管、膈肌和肋间肌，以步骤S2获得的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果作为有限元仿真计算的载荷条件，分别施加在呼吸系统有限元模型的膈肌和肋间肌上，进行呼吸系统动力学仿真计算分析，得到呼吸系统动力学仿真计算分析结果；
进一步的，步骤S3的具体实现方法包括如下步骤：
S3.1、自开源数据库https://www.cancerimagingarchive.net获得CT数据，使用Mimics建立呼吸系统有限元模型，包括肺、气管、膈肌和肋间肌；
S3.2、采用Hypermesh完成有限元模型处理和网格划分，单元类型为四面体单元，气管部分为三角形面网格单元，导入到ABAQUS中赋予组织的材料参数；
S3.3、以步骤S2获得的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果作为有限元仿真计算的载荷条件，将肺门处及肺部里侧进行除z轴位移外五自由度约束，对膈肌的中央部分、上呼吸道部分及肺尖处进行完全约束，进行呼吸系统动力学仿真计算分析，得到呼吸系统动力学仿真计算分析结果；
进一步的，通过CT数据在Mimics中完成呼吸系统有限元模型建立，如图8所示，CT数据来自开源数据库(https://www.cancerimagingarchive.net)，根据解剖学图谱，将气管建立到7级。在Hypermesh中进行气管部分的中空处理，并对模型进行网格划分，单元类型为四面体单元，气管部分为三角形面网格单元，导出到ABAQUS中。各部分的材料参数按照表1赋予。使用步骤2所得计算数据(图5中的计算数据曲线)作为模型的载荷条件，作用类型为压强载荷，根据呼吸运动过程中肺部受力的情况，将肺门处及肺部里侧进行除z轴位移外五自由度约束来近似简化人体内肺门的结缔组织及纵膈胸膜对肺的作用，保证肺运动的平衡。为了模拟横膈的中央肌腱的固定作用及上呼吸道几乎不产生运动的情况，在ABAQUS中分别对膈肌的中央部分、上呼吸道部分及肺尖处进行了完全约束。然后提交计算，得到动力学仿真计算结果，如图9所示；
表1各组织的材料参数
如图11、12所示，呼吸系统动力学仿真计算结果采用体积变化值和选定节点位移进行分析验证，短虚线为仿真计算数据，条形图为文献数据。体积变化值具体为利用ABAQUS的质量属性查询可获得肺部变形前后的体积差值与力学模型的输入肺容积Valv的变化值接近程度，图11结果肺部变形前后体积差为425ml，输入肺容积Valv的变化值490ml，偏差为13％；选定节点位移具体为文献中提供的选定节点的位移数据与仿真计算中对应节点的位移数据的匹配程度。
S4、从步骤S3的呼吸系统动力学仿真计算分析结果中得到气管未变形的有限元模型、气管最大变形时的有限元模型，分别对气管未变形和最大变形时的有限元模型进行流体力学仿真计算；
进一步的，步骤S4的具体实现方法包括如下步骤：
S4.1、通过ABAQUS从呼吸系统动力学仿真计算分析结果中获取气管未变形和最大变形时的有限元模型，然后经Geomagic Wrap、Geomagic Design和SolidWorks实体化处理，导入到ANSYS中进行网格划分；
步骤S4.1中经气管未变形时及最大变形时模型在吸气峰4.02m/s和呼气峰-4.59m/s入口速度下的流体力学仿真计算分析；
S4.2、在Fluent中以入口速度曲线的吸气峰值和呼气峰值为边界条件，进行2次气管未变形的有限元模型、气管最大变形时的有限元模型的流体力学仿真计算，对比2次流体力学仿真计算结果，得到平静呼吸时，流体对于呼吸运动中气管变形的影响忽略不计的结论。
进一步的，步骤S4中，通过ABAQUS获取的气管最大变形时及未变形时的三维模型，在Geomagic Wrap、Geomagic Design和SolidWorks进行实体化处理，导入到ANSYS中进行网格划分，单元类型为四面体网格，边界层第一层厚度为0.1，共5层，增长率1.1，膨胀层边界设置到四级支气管，其余支气管部分不使用膨胀层，进行网格加密。然后在Fluent中求解，边界条件设置具体为取图8所示的入口速度吸气峰4.02m/s和呼气峰-4.59m/s作为两气管模型第一和第二次仿真计算条件，选用压力基及瞬态分析，考虑重力，开启能量方程，选择k-ε(RNG)模型，入口速度因使用绝对速度，入口温度设置为25℃，出口为压力出口，壁面设置为静止壁面且剪切无滑移，壁面温度采用人体正常体温36.5℃，热通量设为0。求解方法中，方案选择为SIMPLE，空间离散梯度为Green-Guass Node Based，空间离散均选用了Second OrderUpwind，压力设为默认，计算时的收敛残差标准设定为小于1e-6。两次仿真计算求解后，将第一和第二次的仿真计算结果进行自身对比，如图13所示为第一次仿真计算，如图14所示为第二次仿真计算。分析结果为平静呼吸时，流体对于呼吸运动中气管变形的影响可以忽略不计，那么动力学仿真计算结果可直接用于虚拟平台上呼吸系统的呼吸运动模拟；
S5、基于第三方虚拟平台Unity 3D，根据步骤S4得到的气管未变形和最大变形时的有限元模型进行流体力学仿真计算结果，将步骤S3得的呼吸系统动力学仿真计算分析结果与Unity 3D进行数据关联，在Unity 3D中进行平静呼吸下人体呼吸运动数值模拟。
进一步的，步骤S5的具体实现方法包括如下步骤：
S5.1、步骤S3得的呼吸系统动力学仿真计算分析结果与Unity 3D进行数据关联的具体实现方法包括如下步骤：
S5.1.1、将步骤S3得到的气管未变形的有限元模型导入到3DMax中按照气管级别进行分割，分割后的部分导入到Unity 3D中；
S5.1.2、针对分割后的每个部分，用Mimics测量气管未变形时自定义点到分割面处的距离，用ABAQUS查询气管变形时自定义点处的位移量；
S5.1.3、根据分割后的每个部分得到的气管未变形时自定义点到分割面处的距离及气管变形时自定义点处的位移量，计算分割后的每个部分的缩放系数，并定义缩放系数的变动范围、变动方式以及变动周期，缩放系数计算公式为：
其中，ζ为缩放系数，S1为变形后节点到局部坐标系原点的距离，S0为变形前节点到局部坐标系原点的距离，Δx为节点变形前后的位移。
S5.2、在Unity 3D中进行平静呼吸下人体呼吸运动数值模拟的方法包括如下步骤：
S5.2.1、设置脚本1获取分割面中心坐标，使用public函数公开功能选项板，通过鼠标点选的方式选取分割部分的分割面中心坐标；
S5.2.1、设置脚本2读取脚本1获得的中心坐标后，使用public函数公开分割部分的最大缩放系数和呼吸运动周期文本输入功能选项板，输入分割部分的最大缩放系数和周期，利用SetEase缓动函数实现线性缩放，SetLoops(-1,LoopType.Yoyo)实现无限循环。
具体的，数据关联方法如图15所示，具体步骤为：
(1)导入步骤S3中气管未变形模型到3DMax中按照气管级别的1、2、3、4级进行分割，分割效果如图16所示，分割后导入到Unity 3D；
(2)针对分割的各个部分，分别测量和查询气管未变形时用户自定义点到分割面处的距离及气管变形时自定义点处的位移量；
(3)根据本步骤(2)的数据，计算每个分割部分对应的缩放系数，并定义缩放系数的变动范围、变动方式以及变动周期。缩放系数变动范围为1到各分割部分计算的缩放系数值，变动方式为x、y、z三轴同步线性变动，变动周期为呼吸周期。
具体的，Unity 3D编程过程为：利用脚本1获取分割面中心坐标，使用public函数公开该功能选项板，用户可依靠鼠标点选的方式选取分割部分的分割面中心坐标；脚本2读取脚本1获得的中心坐标后，使用public函数公开分割部分的最大缩放系数和呼吸运动周期文本输入功能选项板，用户可手动输入分割部分的最大缩放系数和周期，利用SetEase缓动函数实现线性缩放，SetLoops(-1,LoopType.Yoyo)实现无限循环。
根据步骤S3的仿真计算数据可得气管呼吸运动模拟结果，如图17所示，一个为气管未变形时，一个为气管最大变形时，以红色参考线作为对比。按图16的分割部分计算，各个部分的缩放系数见表2，呼吸周期为3.2s。
表2呼吸运动模拟的数据结果
具体实施方式二：
一种平静呼吸下人体呼吸运动数值模拟系统，依托于具体实施方式一所述的一种平静呼吸下人体呼吸运动数值模拟方法实现，包括力学模型计算模块、力学模型拟合模块、呼吸系统动力学仿真计算模块，所述力学模型计算模块、力学模型拟合模块、呼吸系统动力学仿真计算模块之间相互连接，所述呼吸系统动力学仿真计算模块具有Unity 3D关联启动功能。
进一步的，所述力学模型计算模块的输入变量与患者的肺功能检查报告中的潮气量产生关联，关联方式具体为在平静呼吸状态下，Valv等同于潮气量与功能残气量的和，Vc、VRfbs与Valv的关系为Vc＝Valv+VRfbs+100ml，q、VRfbs取值于文献数据，实现了平静呼吸下患者的个性化呼吸运动数值模拟的功能。
平静呼吸下人体呼吸运动数值模拟系统中力学模型计算模块主要负责变量输入和呼吸力学模型的数值计算；力学模型拟合模块主要负责将力计算模块的计算结果进行有效的拟合，以便于ABAQUS动力学仿真计算；呼吸系统动力学仿真计算模块主要完成ABAQUS动力学仿真计算结果的查看，并进入到Unity 3D中完成虚拟气管模型中相关参数的修改，系统界面、交互流程及功能流程如图18、19和20所示。
具体的，步骤S1的力学模型输入变量与患者肺功能检查数据存在以下关系：在平静呼吸状态下，Valv等同于潮气量与功能残气量的和，Vc、VRfbs与Valv的关系可认定为Vc＝Valv+VRfbs+100ml，q、VRfbs取用步骤S2中的文献数据。现获取三个患者的肺功能检查报告单，获得的患者平静呼吸下的潮气量(VT)见表3，并通过模拟系统完成整个流程，根据步骤3中所述的仿真计算分析结果，可得结果见表4和图21所示，缩放系数的计算结果见表5，气管呼吸运动模拟结果如图22所示，由于患者2和患者3的缩放系数差别不大，所以图中上方为患者1，下方为患者2和患者3。
表3患者数据信息
表4患者的求解结果，包括肺容积变化值对比和曲线拟合方法
表5患者个性化呼吸运动模拟的数据结果
/>
需要说明的是，术语“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
虽然在上文中已经参考具体实施方式对本申请进行了描述，然而在不脱离本申请的范围的情况下，可以对其进行各种改进并且可以用等效物替换其中的部件。尤其是，只要不存在结构冲突，本申请所披露的具体实施方式中的各项特征均可通过任意方式相互结合起来使用，在本说明书中未对这些组合的情况进行穷举性的描述仅仅是出于省略篇幅和节约资源的考虑。因此，本申请并不局限于文中公开的特定具体实施方式，而是包括落入权利要求的范围内的所有技术方案。

Claims (9)

1.一种平静呼吸下人体呼吸运动数值模拟方法，其特征在于，包括如下步骤：
S1、根据肺部静力学特性将呼气阶段分为吸气阶段、呼气阶段，在膈肌和肋间肌上建立吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型，然后通过牛顿第二定律构建吸气阶段、呼气阶段的呼吸肌作用力的力学模型，吸气阶段呼吸肌作用力的力学模型为：
其中，Pimf为肋间肌作用力，Rup为上呼吸道产生的气阻，CL为肺顺应性，Cth为胸廓顺应性，PT为跨肺压，q为体积流量，Valv为肺容积，Vc为胸腔容积，PRfbs为无软骨气道组织的压强；
Pdf为膈肌作用力，Cab为腹部顺应性，μ1、μ2、Ap、Bp、Ae、Be、Ke、η1、η2均为常数；
VRmax取值为0.15，VRfbs为余下部分支气管系统容积，Ppl为胸膜腔内压，TLC为肺总容量，RV为余气量，K为修正系数，Pw为以压强形式表达的壁面液面表面张力；
呼气阶段呼吸肌作用力的力学模型为：
其中，P’imf为肋间肌名义作用力，P’df为膈肌名义作用力，RRfbs为无软骨气道组织的气阻，μ3为常数；
S2、将步骤S1得到的吸气阶段、呼气阶段的呼吸肌作用力的力学模型进行计算整合和可靠性验证；
S3、建立呼吸系统有限元模型，包括肺、气管、膈肌和肋间肌，以步骤S2获得的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果作为有限元仿真计算的载荷条件，分别施加在呼吸系统有限元模型的膈肌和肋间肌上，进行呼吸系统动力学仿真计算分析，得到呼吸系统动力学仿真计算分析结果；
S4、从步骤S3的呼吸系统动力学仿真计算分析结果中得到气管未变形的有限元模型、气管最大变形时的有限元模型，分别对气管未变形和最大变形时的有限元模型进行流体力学仿真计算；
S5、基于第三方虚拟平台Unity 3D，根据步骤S4得到的气管未变形和最大变形时的有限元模型进行流体力学仿真计算结果，将步骤S3得的呼吸系统动力学仿真计算分析结果与Unity 3D进行数据关联，在Unity 3D中进行平静呼吸下人体呼吸运动数值模拟。
2.根据权利要求1所述的一种平静呼吸下人体呼吸运动数值模拟方法，其特征在于，步骤S1的具体实现方法包括如下步骤：
S1.1、根据肺部静力学构建跨肺压模型，其中跨肺压为肺组织所受的压力差，其具体实现方法包括如下步骤：
S1.1.1、在跨肺压的作用下，胸膜壁与肺内相互联结的肺泡隔膜壁所受应力之和为0，所以肺组织的径向力平衡公式为：
其中，Pi、Po分别为壁内压强、壁外压强，ΣFo/A、ΣFi/A分别为壁面抗衡壁内压强的内力、壁面抗衡壁外压强的内力；
则对胸膜壁和肺泡壁的力平衡关系，可以得到以下公式：
其中，Ppi、Ppo分别为胸膜壁内压和胸膜壁外压，Σσpi为胸膜壁回缩应力，Σσpo为胸膜壁扩张应力，Pai、Pao分别为肺泡隔膜壁内压和肺泡隔膜壁外压，Σσai为局部回缩应力，Σσao为局部扩张应力；
S1.1.2、根据Dalton气体分压定律和帕斯卡原理，在能量损失忽略不计的情况下，认为跨肺压在肺内的传递不失真，跨肺压表示为以下公式：
PT＝(Ppi-Ppo)+ΣσI-∑σO+Pw
其中，ΣσI＝ΣΣσai为肺总体回缩张力，ΣσO＝ΣΣσao为肺总体扩张张力；
肺总体回缩张力和肺总体扩张张力主要是由肺泡隔膜的弹性应力及外力作用引起的，肺总体回缩张力代表了肺组织的弹性回缩力，肺总体扩张张力代表了肺泡壁隔膜两侧的压力差，在肺部受到呼吸肌牵引而变形引发其内部应力发生改变的时候，肺本身所受内力的作用主要可以分为表面张力和弹性应力，弹性应力表示为：
其中，Nn为肺泡隔膜微元的法向张力，R为肺泡初始半径；
法向张力使用肺的拟应变能函数表达：
其中，c，a1，a2，a4为常数；E1，E2为主应变；ρo为参考状态下的肺密度，ρoW(e)为参考状态每单位肺体积的应变能；
以二维应力状态表示肺泡隔膜的应力-应变关系，肺泡隔膜微元的法向张力为拟应变能函数对应变的导数：
其中，δ为肺泡隔膜的厚度，根据各向同性得到E1＝E2；
S1.1.3、肺部组织为连续介质，依据肺应力-应变关系建立肺组织变形运动的场方程为：
其中，xj为三维笛卡尔坐标x1、x2、x3，v为肺某点的变形速度，ρ为肺组织的密度，为肺某点的位置方程，U为肺部组织的位移场函数，肺部为连续介质，肺组织变形运动的场方程满足应变协调方程；
S1.1.4、利用肺组织和胸壁的弹性回缩力表达肺组织的弹性回缩力Fe为：
胸壁的弹性回缩力Fp为：
最终得放到跨肺压：
S1.2、进行呼吸系统吸气和呼气阶段的作用力分析，然后在膈肌和肋间肌上建立吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型，具体实现方法包括如下步骤：
S1.2.1、整个呼吸系统呼气过程主要受到五个力的作用：肋间肌作用力、膈肌作用力、跨肺压PT、跨膈压Pdi和气阻R，则得到以下计算公式：
PT＝Palv-Ppl
其中，Palv为肺内压；Ppl为胸膜腔内压；
其中，C为顺应性，ΔV为容积变化量，ΔP为压力变化量；
S1.2.2、跨膈压和腹部对于呼吸系统的影响使用腹部顺应性Cab进行替代，平静呼吸时取值为1L/cmH2O，
气阻是气道阻力，属于空气在气管内流动而产生的力，由于软骨的支撑作用，气管的变形量不明显，则上呼吸道产生的气阻由下式表示：
Rup＝μ1+μ2|q|
余下支气管系统没有软骨支撑，但不包含2mm管径以下的部分，余下部分支气管系统气阻RRfbs由下式确定：
S1.2.3、吸气过程中，只考虑上呼吸道气阻，呼气过程主要受三个力的作用，跨肺压、跨膈压和上呼吸道及余下支气管系统部分的气阻，呼吸过程中还要考虑没有软骨包裹的余下支气管系统的弹性，则得的公式为：
S1.2.4、以肺部为基础，分别就肋间肌和膈肌建立吸气和呼气阶段的弹簧阻尼等效模型，肺部的弹簧阻尼模型以顺应性和跨肺压为主要部分，肺组织的粘滞特性代表了阻尼RL，它是肺的拟应变能函数的部分表示，另一部分表示为作用在肺上的弹性回缩力Fe和Fp，Pw和Ppl是跨肺压关系式中的剩余元素，肺部顺应性CL代表肺组织的弹性，作用在肺上的还有余下支气管系统的弹性PRfbs，作用在肋间肌和膈肌之上的力为自身作用力以及上呼吸道的气道阻力，胸廓顺应性Cth表达了胸廓本身的弹性行为，腹部顺应性Cab代表了腹部和跨膈压对于呼吸系统的影响，呼气过程相对于吸气过程而言，呼气过程中膈肌和肋间肌不考虑自身作用力且多考虑了余下支气管系统气阻RRfbs；
S1.3、通过牛顿第二定律，依据吸气阶段弹簧阻尼等效模型、呼气阶段弹簧阻尼等效模型建立呼吸肌作用力的力学模型，力学模型确定的输入变量为肺容积Valv、胸腔容积Vc、余下支气管系统容积VRfbs以及体积流量q。
3.根据权利要求1或2所述的一种平静呼吸下人体呼吸运动数值模拟方法，其特征在于，步骤S2的具体实现方法包括如下步骤：
S2.1、根据步骤S1得到的吸气阶段、呼气阶段的呼吸肌作用力的力学模型的输入变量，利用Origin从文献曲线获取模型求解的初始数据；
S2.2、根据步骤S2.1得到的初始数据使用MATLAB求解步骤S1得到的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果，求解时按照整个周期的呼吸过程进行阶段计算整合；
S2.3、使用MATLAB Cftool拟合步骤S2.2获得的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果，平滑作用力曲线；
S2.4、对比步骤S2.3拟合后的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果和文献曲线，验证吸气阶段、呼气阶段的呼吸肌作用力的力学模型的可靠性。
4.根据权利要求3所述的一种平静呼吸下人体呼吸运动数值模拟方法，其特征在于，步骤S3的具体实现方法包括如下步骤：
S3.1、获得CT数据，使用Mimics建立呼吸系统有限元模型，包括肺、气管、膈肌和肋间肌；
S3.2、采用Hypermesh完成有限元模型处理和网格划分，单元类型为四面体单元，气管部分为三角形面网格单元，导入到ABAQUS中赋予组织的材料参数；
S3.3、以步骤S2获得的吸气阶段、呼气阶段的呼吸肌作用力的力学模型计算结果作为有限元仿真计算的载荷条件，将肺门处及肺部里侧进行除z轴位移外五自由度约束，对膈肌的中央部分、上呼吸道部分及肺尖处进行完全约束，进行呼吸系统动力学仿真计算分析，得到呼吸系统动力学仿真计算分析结果。
5.根据权利要求4所述的一种平静呼吸下人体呼吸运动数值模拟方法，其特征在于，步骤S4的具体实现方法包括如下步骤：
S4.1、通过ABAQUS从呼吸系统动力学仿真计算分析结果中获取气管未变形和最大变形时的有限元模型，然后经Geomagic Wrap、Geomagic Design和SolidWorks实体化处理，导入到ANSYS中进行网格划分；
S4.2、在Fluent中以入口速度曲线的吸气峰值和呼气峰值为边界条件，进行2次气管未变形的有限元模型、气管最大变形时的有限元模型的流体力学仿真计算，对比2次流体力学仿真计算结果，得到平静呼吸时，流体对于呼吸运动中气管变形的影响忽略不计的结论。
6.根据权利要求5所述的一种平静呼吸下人体呼吸运动数值模拟方法，其特征在于，步骤S4.1中经气管未变形时及最大变形时模型在吸气峰4.02m/s和呼气峰-4.59m/s入口速度下的流体力学仿真计算分析。
7.根据权利要求6所述的一种平静呼吸下人体呼吸运动数值模拟方法，其特征在于，步骤S5的具体实现方法包括如下步骤：
S5.1、步骤S3得的呼吸系统动力学仿真计算分析结果与Unity 3D进行数据关联的具体实现方法包括如下步骤：
S5.1.1、将步骤S3得到的气管未变形的有限元模型导入到3DMax中按照气管级别进行分割，分割后的部分导入到Unity 3D中；
S5.1.2、针对分割后的每个部分，用Mimics测量气管未变形时自定义点到分割面处的距离，用ABAQUS查询气管变形时自定义点处的位移量；
S5.1.3、根据分割后的每个部分得到的气管未变形时自定义点到分割面处的距离及气管变形时自定义点处的位移量，计算分割后的每个部分的缩放系数，并定义缩放系数的变动范围、变动方式以及变动周期，缩放系数计算公式为：
其中，ζ为缩放系数，S1为变形后节点到局部坐标系原点的距离，S0为变形前节点到局部坐标系原点的距离，Δx为节点变形前后的位移；
S5.2、在Unity 3D中进行平静呼吸下人体呼吸运动数值模拟的方法包括如下步骤：
S5.2.1、设置脚本1获取分割面中心坐标，使用public函数公开功能选项板，通过鼠标点选的方式选取分割部分的分割面中心坐标；
S5.2.1、设置脚本2读取脚本1获得的中心坐标后，使用public函数公开分割部分的最大缩放系数和呼吸运动周期文本输入功能选项板，输入分割部分的最大缩放系数和周期，利用SetEase缓动函数实现线性缩放，SetLoops(-1,LoopType.Yoyo)实现无限循环。
8.一种平静呼吸下人体呼吸运动数值模拟系统，依托于权利要求1-7之一所述的一种平静呼吸下人体呼吸运动数值模拟方法实现，其特征在于，包括力学模型计算模块、力学模型拟合模块、呼吸系统动力学仿真计算模块，所述力学模型计算模块、力学模型拟合模块、呼吸系统动力学仿真计算模块之间相互连接，所述呼吸系统动力学仿真计算模块具有Unity 3D关联启动功能。
9.根据权利要求8所述的一种平静呼吸下人体呼吸运动数值模拟系统，其特征在于，所述力学模型计算模块的输入变量与患者的肺功能检查报告中的潮气量产生关联，关联方式具体为在平静呼吸状态下，Valv等同于潮气量与功能残气量的和，Vc、VRfbs与Valv的关系为Vc＝Valv+VRfbs+100ml，q、VRfbs取值于文献数据，实现了平静呼吸下患者的个性化呼吸运动数值模拟的功能。
CN202310283842.1A 2023-03-22 2023-03-22 平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法 Active CN116525118B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310283842.1A CN116525118B (zh) 2023-03-22 2023-03-22 平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310283842.1A CN116525118B (zh) 2023-03-22 2023-03-22 平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法

Publications (2)

Publication Number Publication Date
CN116525118A CN116525118A (zh) 2023-08-01
CN116525118B true CN116525118B (zh) 2023-09-26

Family

ID=87400155

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310283842.1A Active CN116525118B (zh) 2023-03-22 2023-03-22 平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法

Country Status (1)

Country Link
CN (1) CN116525118B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106446572A (zh) * 2016-09-27 2017-02-22 上海精劢医疗科技有限公司 基于边界元模型和局部区域修正的肺部呼吸运动获取方法
CN114913752A (zh) * 2022-05-26 2022-08-16 中国人民解放军陆军军医大学 一种基于集总参数的人体呼吸系统模型
WO2022264045A1 (en) * 2021-06-15 2022-12-22 Fisher & Paykel Healthcare Limited Patient simulation training system for a breathing assistance or respiratory apparatus
CN115563782A (zh) * 2022-10-11 2023-01-03 中煤科工集团重庆研究院有限公司 多因素影响下综掘面呼吸性粉尘分布模型的构建方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7544137B2 (en) * 2003-07-30 2009-06-09 Richardson Todd E Sports simulation system
FR3071398A1 (fr) * 2017-09-22 2019-03-29 Universite de Bordeaux Procede de simulation d’une dynamique respiratoire d’un poumon virtuel, simulateur virtuel, ensemble respiratoire.

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106446572A (zh) * 2016-09-27 2017-02-22 上海精劢医疗科技有限公司 基于边界元模型和局部区域修正的肺部呼吸运动获取方法
WO2022264045A1 (en) * 2021-06-15 2022-12-22 Fisher & Paykel Healthcare Limited Patient simulation training system for a breathing assistance or respiratory apparatus
CN114913752A (zh) * 2022-05-26 2022-08-16 中国人民解放军陆军军医大学 一种基于集总参数的人体呼吸系统模型
CN115563782A (zh) * 2022-10-11 2023-01-03 中煤科工集团重庆研究院有限公司 多因素影响下综掘面呼吸性粉尘分布模型的构建方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
人体呼吸系统力学属性下呼吸动力过程的仿真研究;刘国辉;肖华军;于立华;;科学技术与工程(第29期);313-318 *

Also Published As

Publication number Publication date
CN116525118A (zh) 2023-08-01

Similar Documents

Publication Publication Date Title
Wall et al. Fluid–structure interaction in lower airways of CT‐based lung geometries
Ma et al. An anatomically based hybrid computational model of the human lung and its application to low frequency oscillatory mechanics
EP1551293B1 (en) System and method for the automatic detection of the expiratory flow limitation
Xi et al. Parametric study on mouth–throat geometrical factors on deposition of orally inhaled aerosols
Bates et al. The effects of curvature and constriction on airflow and energy loss in pathological tracheas
Bertram Flow-induced oscillation of collapsed tubes and airway structures
Wang et al. Numerical analysis of respiratory flow patterns within human upper airway
Baffico et al. Multiscale modeling of the respiratory tract
Shi et al. Pressure dynamic characteristics of pressure controlled ventilation system of a lung simulator
Kannan et al. A quasi‐3D compartmental multi‐scale approach to detect and quantify diseased regional lung constriction using spirometry data
Rios et al. Computational fluid dynamics analysis of surgical approaches to bilateral vocal fold immobility
Smith et al. Influence of subglottic stenosis on the flow-induced vibration of a computational vocal fold model
CN116525118B (zh) 平静呼吸下人体呼吸运动数值模拟系统及数值模拟方法
Ionescu et al. Viscoelasticity and fractal structure in a model of human lungs
Nucci et al. A morphometric model of lung mechanics for time-domain analysis of alveolar pressures during mechanical ventilation
Choi et al. Multiscale numerical analysis of airflow in CT-based subject specific breathing human lungs
CN114913752B (zh) 一种基于集总参数的人体呼吸系统模型
Furman et al. Mathematical model of breath sound propagation in respiratory tract
Lumb et al. Computational fluid dynamic modelling of the effect of ventilation mode and tracheal tube position on air flow in the large airways
Kageyama et al. Fluid dynamic assessment of positive end-expiratory pressure in a tracheostomy tube connector during respiration
Wang et al. Analysis of airway model of tumor patients based on CFD
Dyachenko et al. Generalization of the mathematical model of lungs for describing the intensity of the tracheal sounds during forced expiration
Giamagas INTEGRATED COMPUTATIONAL SIMULATION OF THE MECHANICS OF BREATHING
Mesit et al. Simulation of lung respiration function using soft body model
Μονοκρούσος et al. Integrated computational simulation of the mechanics of breathing

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
GR01 Patent grant
GR01 Patent grant