赵士明 李文雷 赵静一③ 郭 锐 高 远* 卢子帅

(* 燕山大学河北省重型机械流体动力传输与控制实验室 秦皇岛 066004)(**先进锻压成形技术与科学教育部重点实验室(燕山大学) 秦皇岛 066004)(***唐山工业职业技术学院机械工程系 唐山 063299)(****河北省唐山市第二人民医院创伤外科研究所 唐山 063001)

0 引 言

组织工程(tissue engineering)是利用工程学和生命科学的基本原理，开发能恢复、维持或改善受损组织或器官功能的生物代替物。组织工程支架作为细胞和新生组织的载体要保证营养液、生长因子等能够从多孔或微管进入支架内部，与细胞进行物质交换排出代谢物[1，2]。组织工程支架为细胞接种、繁殖和保持各自功能提供了保障，其可以引导组织再生，控制组织结构。研究表明，流体剪应力等各种外载力的刺激，对骨细胞的增殖分化具有重大影响。支架内部微孔的尺寸、空间走向、结构及布局都将影响其内部微流体流场的分布。营养液在管道里究竟能给多大范围内的细胞提供足够的养分,与支架的内部空间结构、孔径和微流体的流速等有密切关系[3]。

现有的骨组织工程支架在实验过程中存在局部新生组织生长较快的问题，以致于营养液、生长因子和代谢物的流动被阻塞，其主要原因是未对骨支架内部流场的作用进行研究。目前，针对支架内部的速度场和压力场分布的研究已有报道，文献[4]通过实验和计算流体动力学数值模拟方法研究了骨组织工程支架内部流场分布，2种方法所得结果相近，能够相互验证，为量化骨组织工程支架内部流体力学环境和细胞生长、分化的关系提供了基础。文献[5]提出了将曲率驱动细胞增长与流体动力学结合的模型。燕山大学赵静一团队提出了组织工程中的微流体现象[3]，对近几年骨组织工程支架内微孔结构[6]、微流体流场研究[7]和微流体驱动控制研究[8]进行了综述，并对骨组织支架内部微孔结构进行了设计和仿真研究[9]。文献[10]研究了3种微管道中细胞浓度分布，并分析了微管道系统中细胞可能沉积的部位。生物体内流体剪应力、静水压力、机械拉伸力等机械力共同作用是骨组织生长和维持的关键，其中流体剪应力对细胞、组织的刺激更为突出[11]。在组织体外培养中各种力学刺激是必不可少的，得到越来越多的重视[12，13]。张寅星等人[14]研究了不同剪应力水平下NO、ALP的含量，发现较低剪应力条件下，成骨细胞的分化与剪应力成一定的比例关系；而在高剪应力水平下，两者之间没有明显的比例关系。文献[15]分析了不同流场下的剪应力分布，在震荡壁面边界下，剪应力水平明显提高。然而，目前对骨组织工程多孔支架结构内部流场及剪应力分布的研究较少。

本文根据工程学和生物医学原理，设计了4种不同内部结构的骨组织工程支架，分析其在不同边界条件作用下，不同支架内部速度、压力、流体剪应力的大小及分布情况，通过理论分析和数值仿真发现，针对人体各部分组织结构的不同，设计不同结构形式的支架，能够调整营养液的输送速度，通过调节流体剪应力的范围，可以增加骨细胞的生长速度，该研究具有一定的理论和实际意义。

1 建立骨组织工程支架物理模型和数学模型

1.1 物理模型

为了确定合适的组织工程支架内部结构，对不同结构内部微流体流场,使用CREO3.0建立了4种三维模型，并使用ICEM CFD进行网格划分如图1～图4所示。在人体骨组织中具有水平的浮克曼管和竖直的哈佛管，中心主管通过横向的浮克曼管与竖直的哈佛管相连，依据这种特征设计了3种仿生结构支架和一种螺旋支架。4种支架中图1是最基础的仿生模型，主管道直径为2 mm，长度为20 mm，支管直径为1.2 mm。图2和图3是在支架1基础上的改进型，而图4是从工程流体力学角度出发设计的内部结构。图1、图2、图3所示3种支架内部主管和支管尺寸相同，图1、图3支管与主管道垂直，图2浮克曼管与主管道夹角呈30 °，图3将哈佛管移至浮克曼管端部并且将每层浮克曼管间距离增大1.5倍，增加了连接哈佛管的辅助支管。图4中螺旋管道支管直径为1.2 mm。

网格划分是计算流体力学数值模拟中重要的一环，生成一套高质量网格将显著提高计算精度和收敛速度。ICEM CFD因其友好的操作界面、完善的几何功能、灵活的拓扑创建、先进的O型网格技术得到快速发展。图1、图2、图3三种支架因其结构相对复杂、相贯界面较多，划分结构网格将非常复杂，消耗大量时间，所以采用非结构网格。图4结构相对简单，使用结构网格，可以大幅度提高网格质量和计算时间。

图1 仿生支架1

图2 仿生支架2

图3 仿生支架3

图4 螺旋结构支架4

1.2 数学模型

骨组织工程支架内的微流孔道虽然达到了微米级别，在其流动特性上呈现了与宏观的不同，但是根据前人的大量研究来看，支架孔道内的流动仍然满足宏观的物理守恒定律，包括质量守恒定律和动量守恒定律和能量守恒定律。孔内流体流动的质量守恒方程即连续方程可用式(1)表示；动量方程的本质是满足牛顿第二定律，可用式(2)表示；能量守恒定律本质是热力学第一定律，可用式(3)[7]表示。营养液粘度较高，流动速度低，由雷诺公式(4)可知Re≈12，其液在支架内的流动状态为层流，其他管道与主管道相比直径小，所以雷诺数更低。假设营养液为不可压缩流体，其运动微分方程可以简化为式(5)的形式[7]。骨组织工程支架内的微流孔道达到了微米级别，在其流动特性上呈现了与宏观的不同。

(1)

式中，ρ为流体密度，单位为kg/m3。

(2)

式(2)中，τxx、τxy、τxz等是因粘性作用而产生的作用在微元体表面上的粘性应力τ的分量；fx、fy、fz为x、y、z方向上的单位质量力。

(3)

其中，E为流体微团的总能，h为焓，keff为有效热传导系数，Sh为体积源项。

(4)

式(4)中，u为平均流速，单位为m/s；d为管道内径，单位为m；μ为动力粘度系数，单位为kg/m·s。

(5)

式中，p为流体微元表面压力，v为流体运动粘度。

流体在流动过程中因受到管壁的摩擦阻力和管道直径变化的影响导致流动方向改变而产生流动阻力，主要有沿程阻力和局部阻力。由沿程阻力引起的能量损失与流程长度成正比，简称沿程损失，用式(6)计算；过流断面的大小、形状、方位等发生变化时流体所受到的阻力是各种各样的，集中在很短的流段内，这种阻力为局部阻力，所引起的能量损失可用式(7)表示。

(6)

(7)

式中，λ为沿程阻力系数，l为流程，d为管道内径，v为流速，ζ为局部阻力系数。

2 通过有限元进行数值模拟与分析

根据所建立的物理模型和数值模型，对骨组织工程支架内部流动的压力与流速进行了数值模拟分析。

2.1 压力和速度分析

把ICEM CFD中划分好的网格文件导入到Fluent中进行仿真研究，营养液的密度ρ为1×103kg/m3，粘度μ为0.0018 Pa·s为不可压缩牛顿流体。假设流体与壁面的接触边界是无滑移壁面，选取管道入口端面为速度入口边界条件，营养液流速恒定为0.01 m/s，出口边界压力为0 Pa。仿生骨支架速度分布云图，如图5～图8所示。

图5 仿生支架1内流体速度分布云图

图6 仿生支架2内流体速度分布云图

图7(a)、7(b)为支架3在主视面和左侧面的速度分布图。由图5、图6、图7速度分布云图可知，这3种结构支架沿着主管道从入口到出口速度变化明显，由于沿程阻力和局部阻力以及各个支管的分流作用，速度逐渐减小。在第一层浮克曼管处的分流现象明显，分流速度接近入口速度的一半。支架2和支架1、支架3相比较，其浮克曼管与主管夹角呈60 °，流体的动量与分流支管的夹角小，流体流动的局部阻力小，所以在支架2中支管的流速明显高于支架1、支架3。正是因为流体动量的矢量性以及主管和哈佛管之间距离不足，无法完全消除向下的动量分量，支架1出现在第一层浮克曼管处向上的速度比向下的速度低。由于支架1的浮克曼管的另一端直接出口分流作用明显，而支架2、支架3的哈佛管与浮克曼管的另一端相连，支架2的哈佛管中流体的流速高于支架1。支架3第一层浮克曼管与哈佛管连接位置向上出口的距离小于下出口，阻力比较小，大部分流体从上出口流出。由图8所示螺旋支架的速度分布云图可知，其流动速度从上到下逐渐变缓，流动状态均匀，脉动小。

(a) 主视

(b) 左视

(a) 纵截面

(b) 横截面

图8所示支架4的最高流速在螺旋管中间达到了0.0196 m/s，支架4内部螺旋管之间相互独立并不存在分流现象。由于螺旋管不存在交叉口没有局部阻力，但是沿程阻力会有所提高。流体在螺旋管中由于离心力的作用流速逐渐提高，阻尼随之增大，直到阻力与离心力平衡。在4种结构中支架4的流场分布更为均匀，并且可以通过减螺距、增加螺旋管的数量提高孔隙率。

仿生骨支架压力分布云图如图9～图12所示。

由于出口边界条件为零压，从整体趋势上看4种支架从入口到出口压力逐步降低到0。支架1中压力作用范围最小，只在入口附近存在较高的压力。在支架2、3中哈佛管和浮克曼管交叉汇合，并且只在上下两端存在出口，造成较高的压力以及较大的作用范围。对于支架3第一层交叉口处的压力梯度高也造成了向上出口的流速高。在支架4中流体的流程远大于其他3种结构，因此沿程阻力大，压力损失高，会在入口处形成高压，并且压力梯度大，3层螺旋管中可以明显看出从外层到内层压力逐步降低。

2.2 流体剪应力分析

体外培养过程中流体剪应力难以测量，实验复杂，数值模拟为支架内流体剪应力的分布及水平提供了良好的条件。牛顿流体对壁面的剪应力可用如下公式计算：

(8)

式中，τ为流体剪应力，μ为动力粘度系数，vx为流体在x轴方向的速度分量，vy为流体在y轴方向的速度分量。

由式(8)可以看出剪应力大小与流体动力粘度、速度梯度成正比。目前，壁面流体剪应力的调制主要有增加流体粘度、增大入口速度以及在支架内设置可动的芯片，但细胞间传质有最优速度范围，营养液有最佳粘度等限制。通过改变支架的进出口边界条件，使用Fluent UDF自定义函数功能改变速度梯度，在不增加近壁面速度的条件下进行流体剪应力的调制。

由文献[5,11]可知支架内平均剪应力0.232～0.304 Pa时能显著促进细胞分化,平均剪应力低于0.1 Pa影响了细胞分化的速度，而剪应力大于0.304 Pa时，细胞分化水平不再明显变化。图13～图16为4种支架结构在优化参数后边界条件下的剪应力分布曲线。通过曲线对比，调制后4种结构支架的流体剪应力显著提高，接近或达到最优的范围，支架1、2、3剪应力提高了0.1 Pa左右，从入口到出口逐渐降低，与速度分布趋势一致。对于支架4流体剪应力提高了0.2 Pa，在很小的范围内波动，整体上剪应力分布均匀,支架4达到了对细胞分化最优的剪应力。

图9 仿生支架1内流体压力分布云图

图10 仿生支架2内流体压力分布云图

(a) 主视

(b) 左视

(a) 纵截面

(b) 横截面

图13 仿生支架1剪应力调制曲线图

图14 仿生支架2剪应力调制曲线图

图15 仿生支架3剪应力调制曲线图

图16 螺旋支架4剪应力调制曲线图

3 结 论

使用CREO3.0通过仿生学与螺旋理论建立了4种整体结构尺寸相近、内部微孔道布局走向不同的骨组织工程支架三维模型；运用N-S方程结合流体力学管路能量损失理论建立了4种支架的数学模型；使用ICEM CFD划分网格，并利用Fluent进行了数值仿真，得到了组织工程支架内部的流动状态，并对骨支架内微流体剪应力进行了优化，所得结果如下。

(1) 支架2的速度场分布更均匀，优于支架1和3。支架2哈佛管流速高于支架1和3，提高了营养物质的交换及代谢废物的排出速度；螺旋支架4内部速度场由于离心惯性力和沿程阻力的共同作用，流速较高，整个速度场分布均匀，但是流程长。

(2) 4种支架的压力流场中支架4更均匀，由入口到出口压力逐步降低，支架1的压力作用面积最小，支架4的压力值最大，内外层压力梯度较高。

(3) 改变进出口边界条件进行流体剪应力的调制使其达到对细胞增殖、分化的最佳值。通过流体剪应力曲线对比发现，支架1、2、3的剪应力从入口到出口逐步减小，调制后剪应力提高约2倍；支架4剪应力分布均匀，调制后提高到0.3 Pa，提高了成骨细胞的生长速度。

使用ANSYS Fluent对所设计的4种骨组织工程支架内部微流体流场的模拟计算与分析，通过支架结构选择和剪应力的优化，可以使支架内部的微流体流场分布更均匀，剪应力更符合人体组织生长需要，进而更加适合骨细胞的增殖分化、营养物质的交换及代谢废物的排出，可为组织工程骨支架的设计和使用提供一定的理论依据。