多级仿生建模及优化技术

喻长江 戴 宁 李大伟 程筱胜

南京航空航天大学机电学院，南京,210016

0 引言

晶格结构材料是极具潜力的轻质多功能结构材料，与传统材料相比，晶格结构材料拥有更加优良的力学性能。晶格结构材料在生物医疗、航空航天等领域应用潜力巨大，增材制造的发展为晶格结构的制造提供了保证。自然界中的生物经过数万年进化，拥有了独特的结构及力学特性，利用其复杂的多尺度和多相结构可实现超出人造系统的功能[1]。这些多层生物结构材料往往比单一结构材料具有更合理的力学特性。以常见的多级生物体牙齿、骨骼、珍珠等[2]为例，这些生物结构材料拥有卓越的力学特性，它们拥有坚硬的表面层，以抵抗磨损或者穿透，并具有坚韧的内在结构以适应变形的需求，这些特性与生物体内部的多级结构有非常紧密的联系。但如何参数化地设计多级结构使其拥有优良的力学性能仍具有一定的挑战性。

SIGMUND等[3]提出模型微观尺度结构优化的概念，初步探索了单晶格阵列对模型宏观弹性模量的控制效果；ARMILLOTTA等[4]利用晶格结构设计并制造出了超轻且坚硬的结构；GARCIA等[5]通过调整晶格单元结构的孔隙率对结构进行拓扑优化，进而调整了晶格的弹性力学性能,但由于单晶格结构功能单一,因此无法设计具有复合功能的结构材料;SCHUMACHER等[6]将多晶格在三维模型空间中优化分布，实现了对模型弹性模量的准确控制，但是该方法计算量大，只能够实现介观尺度的建模；XU等[7]在文献[6]的基础上，通过在单元晶体中组装两个不同的晶格单元，实现在不同空间方向上对刚度的调控，生成具有可控各向异性的晶格结构单元，但这类设计方法只能够实现性能的微调。

本文以多级生物体作为主要研究对象，从微观层面探究其坚韧的原理，构建了以晶格孔隙率驱动力学性能调控的晶格模型族，对目标对象进行填充优化，并基于有限元映射的方法使晶格呈现梯度变化，从而最大化地接近结构梯度变化特征。此外，还基于隐式曲面方法优化了晶格间的过渡连接，保证了结构间的有效连接及力的稳定传递。通过迭代优化生成具有形状结构稳健且满足坚韧生物力学特性的模型。最后将物理实验与仿真实验相结合，评估了本文多级晶格力学性能的控制方法，并将该方法应用在生物模型上。

1 技术路线

本文目标是将多级体空间根据给定的力学参数用匹配的晶格进行填充，并依据仿真结果进行优化设计，从而生成期望的复合功能结构。方法主要分为两个阶段，在预处理阶段，通过有限元仿真运用插值的方法，建立由弹性模量和抗冲击韧性等参数进行检索的晶格单元模型族。在设计过程中，首先将获得的生物体μCT(computed tomography)图片进行重建，获得多级体空间，根据前人对生物体各级力学性能的研究结果，输入各级体空间的生物力学参数，匹配晶格模型族中的单元结构，根据给定的力学参数对模型各级空间进行填充，通过全局优化的方法对模型进行初步的优化设计，并保证连接的有效性及稳定性。随后将模型导入有限元分析软件ABAQUS中，根据实际的受力情况对模型进行分析，并将模型的力学结果导出。通过三维映射的方法，将分析结果与晶格填充模型进行匹配，从而调整晶格结构的孔隙率、尺寸等参数。将优化后结果重新导入ABAQUS中进行分析，验证是否满足生物体变形的需求，如果不满足则返回重新完成优化过程，循环往复，直到生成具有形状结构稳健和满足生物力学特性的模型，其流程如图1所示。

图1 多级微结构力学性能建模流程图Fig.1 The mechanical properties of multilevel microstructures modeling flowchart

2 晶格结构设计

2.1 均质化

目前人们已经对周期性中尺度单元结构设计提出了各种优化解决方法，但是设计人员必须考虑均质化问题，即单胞尺寸必须比所有方向上的设计空间小得多。假设在所有方向上都是周期性的晶格结构，可以通过均质化的方法计算出晶格单元的弹性模量[8]。由晶格单元构成的刚度矩阵Μ可以表示为

(1)

式中,n为晶格单元体网格的个数;R为应力应变矩阵;Ci为与单元密度相关用于表示晶格的矩阵;Ωi为第i个单元域。

我们定义晶格单元Cell为c，通常为了计算出晶格单元的结构特性，需要求解6个方向的载荷：

MX(i)=F(i)i=1，2，…，6

(2)

其中，X(i)为每个晶格单元c在i方向上的位移向量；F(i)为施加的载荷，可以通过应力场计算得出：

(3)

其中，单元应变矩阵为

ε(1)=[1 0 0 0 0 0]T
ε(2)=[0 1 0 0 0 0]T
ε(3)=[0 0 1 0 0 0]T
ε(4)=[0 0 0 1 0 0]T
ε(5)=[0 0 0 0 1 0]T
ε(6)=[0 0 0 0 0 1]T

2.2 晶格模型族构建

在均值化理论基础之上进行晶格单元的设计，图2所示为几种常见的晶格单元结构。为高效准确地检索到具有特定参数的微结构，采取了预计算的策略。考虑到晶格单元间的连接问题，将晶格单元结构划分为三类，分别为面心立方(face center cubic，FCC)、体心立方(body center cubic，BCC)以及面体心立方(face and body center cubic，FBCC)，以保证相邻结构可以稳定地连接。随后对每一个晶格结构都进行有限元分析，以获得不同晶格结构所对应的力学参数，并为构建的单元结构赋予对应的参数信息，如晶格类别、孔隙率、抗冲击韧性、弹性模量等，建立晶格与各参数之间的映射关系。通过对晶格孔隙率的调控，可以获得一组由孔隙率驱动的晶格单元所对应的弹性模量，对这组数据进行插值，可以获得一组连续的晶格模型系。首先，利用距离场转换函数：

f(x)=sgn(x)ln(|x|+δ)

(4)

(a)十字立方 (b)面心立方

(c)十字面心立方 (d)体心立方

(e)立方体 (f)正十二面体

(g)面体心立方 (h)阵结构 图2 典型的晶格单元结构图Fig.2 Typical lattice unit structures

将模型空间进行变换。其中，δ为调控拟合精度的阈值，δ的数值越小，拟合的精度越高。然后根据给定的孔隙率φgoal和变换距离场f(x)，通过最小二乘法线性插值求解ξ，当ξ最小时可以求得φgoal以基本单元cbase插值计算得到的对应单元体结构cgoal：

(5)

式中,d为空间维数;m为多项式的次数;xi为x的第i个取值;Π为所有晶格点阵的所在域。

通过这种方法可以获得比较准确的晶格族系，后期可以通过继续增加晶格种类来提高方法的适用性及准确性，为多级参数力学性能调控的选择提供丰富的初始样本。图3展示了三类典型晶格结构BCC、FCC以及FBCC，由孔隙率驱动的五种典型结构弹性模量的变化趋势。由力学试验可知，不同孔隙率的晶格结构所对应的弹性模量呈线性分布。当单元模型族数据库越完善，则在建模时对晶格单元的选择就越多，便于设计出更具个性化需求的多级模型。

(a)BCC晶格

(b)FCC晶格

(c)FBCC晶格图3 三类常见的三维晶格单元结构，不同孔隙率所对应的弹性模量Fig.3 Three common three-dimensional lattice unit structures, different porosity corresponding to the elastic modulus

2.3 多级晶格单元力学调控方法

在单元晶格的理论研究基础之上，笔者提出了多级晶格结构的力学调控方法。从根源上可以追溯到线性弹性定律：σ=Eε，其中σ={σ1,σ2,σ3,σ4,σ5,σ6},ε={ε1,ε1,ε1,γ4,γ5,γ6}。在上文中定义了晶格单元的均值化计算方法，首先将晶格单元拓展到单级晶格模型中，定义m为单级晶格层数，n为m层晶格的个数，利用均质化理论[9]将单级晶格弹性模量表示为

(6)

ξ(i,j),(m,n)=exp(d(i,j),(m,n)·g(i,j),(m,n))

(7)

描述第i层第j个晶格单元与第m层第n个晶格单元间的匹配程度，d(i,j),(m,n)为选定结构与实际微结构间的连接误差,g(i,j),(m,n)为相邻结构间的相似度。随后对多级结构进行定义：

(8)

式中,N为模型的层数。

3 多级微结构填充实现

3.1 多级晶格填充

在对多级生物体进行仿生设计之前，必须先获得多级生物体空间的离散网格数据模型，然后将μCT扫描数据进行三角网格数据重建。由于生物体各层级结构之间是梯度变化的，没有明确的分界线，为便于后期设计，在根据实际情况对模型重建后，需进行相应的优化设计，运用等值面绘制方法，将各层结构用明确的界限划分出来，用几何参数Ωi(i∈N)表示各级体空间，其中i为生物体层级数。输入各层级需要设定的力学参数，在晶格族数据库中检索与之匹配的晶格单元，随后进行平铺填充。

3.2 基于水平集方法的梯度建模

晶格结构间连接的稳定性，将对多级结构的宏观力学性能产生巨大影响。如图4a所示，晶格结构间不可能完全匹配。为解决这一问题，一方面需丰富晶格模型族备选晶格，选择匹配程度系数ξ较高的晶格。在检索晶格结构时匹配类别能够相互连接的结构，比如将FBCC-BCC、FBCC-FCC、FBCC-FBCC等结构相连。值得注意的是，FCC不能与BCC直接相连，必须以FBCC作为过渡晶格结构。另一方面，基于隐式曲面的方法[6]对晶格进行融合。由于距离函数的变化率是均匀的，因此有利于实现晶格梯度均匀变化，并提高数值计算的稳定性。图4b为晶格融合后的效果，图4c为晶格梯度变化的效果。

(a)四种常见晶格的连接 (b)晶格融合后结果

(c)晶格杆径从3.0 mm到0.5 mm梯度变化图4 晶格融合及梯度变化Fig.4 Lattice fusion and gradient change

3.3 性能空间分布

首先利用OOFEM有限元库[10]的分析方法将应力映射到模型上，获得多级映射空间Ωi(i∈N)。然后，通过预计算选取不同空间Ωi对应的相似力学特性单元模型Li,j(i,j∈N)，并用求得的单元结构ci,j对模型空间Ωi进行填充，其关系为

(9)

为介绍基于有限元的纹理映射原理，以15 mm×7.5 mm×2.5 mm长方体为实验物体，进行有限元仿真实验，对长方体的中间施加压力(图5a)，获得图5b所示应力场云图。采用LU等[11]提出的基于全局对比度保持转换方法，构建如下参数化转化函数：

I=f(c,s)
c=(r,g,b)

式中，r、g、b为每个像素的颜色向量；s为该像素点的显著性度。

将彩色的应力分布云图转换为灰度，并保持良好的灰度图，如图5c所示。在上文梯度变形的原理之上，构建应力云图灰度值与模型梯度Osolution之间的映射关系为

G(I)=Osolution

(10)

根据不同的灰度值调整晶格杆径的梯度变化，以匹配不同应力值的大小，效果见图5d。

(a)受力分析

(b)应力云图

(c)灰度图

(d)图像映射结果图5 有限元结果映射实验Fig.5 FEA results in mapping experiment

4 实验分析

4.1 压力实验

为了验证多级微结构对弹性模量参数调控的有效性，我们从排列组合方式、方向以及孔隙率等多个角度，对多级晶格结构进行了物理实验以及力学仿真实验，并将实验与仿真结果进行了对比。如图6a所示，运用熔融沉积制造(fused deposition modeling，FDM)工艺使用聚乳酸(PLA)材料制造了三组对比实验样件，三组多级晶格结构单元的排列方式分别是B-FB-F(BCC-FBCC-FCC)、B-FB-F(vertical)、F-FB-B，每组设计6个单元，杆径从0.45～0.20 mm均匀变化，且每个多级单元结构由3×3×3个单级单元结构所组成。在结构测试过程中，使用Instron 3367的设备进行压力试验，如图6b所示。试验前，先对模型施加了一个预载荷，使晶格处于预紧状态,然后以0.8 mm/min的速率进行测试。

(a)多级晶格模型样件

(b)力学压缩试验设备图6 多级晶格结构物理压力试验Fig.6 Physical pressure test of multilevel lattice element structure

(a)φ=68.25%

(b)φ=55.02%图7 B-FB-F晶格单元不同孔隙率应力-应变曲线图Fig.7 A set of stress-strain curves for a group of B-FB-F lattices

将压力测试的应力-应变曲线线性部分进行线性拟合,图7为一组B-FB-F结构弹性模量的拟合效果图，其中实线表示真实应力-应变曲线数值变化,拟合出的虚线斜率为弹性模量的数值。

图8a所示为多级晶格结构，考虑到其各向异性，在三维空间中可能受到来自各个方向施加的力，为保证设计后的结果与真实受力的一致性，将有限元仿真结果与物理试验结果进行对比验证，从A、B、C三个方向对晶格进行压力试验，图8b～图8d显示了三组仿真及物理试验的结果对比。通过表1可知，仿真和试验结果的拟合结果较好，所以本文利用仿真实验的结果是具有可行性的。

(a)施加力的方向

(b)F-FB-B晶格

(c)B-FB-F晶格

(d)F-FB-B晶格垂直方向图8 三种结构仿真结果与物理压缩试验数据对比Fig.8 Results of comparison between simulation and physical tests

结构F-FB-BB-FB-FF-FB-B(V)相似度87.7190.3688.27

4.2 冲击试验

利用光固化成形(stereolithigraphy apparatus,SLA)工艺制造了4组共24个试验样件(80 mm×10 mm×6 mm)，每个长条样件由360个晶格单元组成，样件杆径在0.5～1.0 mm之间均匀变化。使用DELTA-TPO摆锤式冲击试验机评估多级晶格模型的抗冲击性能，也就是晶格模型的韧性，设置试验温度为2 ℃，摆锤的冲击能量为25 J，每次都对摆锤的角度进行调零并保证样件左右对称均匀摆放，记录下每组模型的抗冲击韧性ak值。

4.3 冲击试验结果

(a)FCC晶格

(b)BCC晶格

(c)FBCC晶格

(d)多级晶格单元图9 单级及多级晶格单元结构冲击试验结果Fig.9 Impact test results of single-stage and multi-stage lattice cell

将不同类别的晶格单元结构、不同孔隙率下的冲击试验结果在图9中用黑色圆点表示，并对试验数据进行线性拟合。FCC、BCC单元晶格结构和多级晶格单元结构的孔隙率在0～1之间呈波浪形变化，而FBCC单元孔隙率在0～0.36之间呈线性变化(图9c)。试验结果整体呈现了很好的规律性，且冲击韧性结果虽然存在一定的误差，但却是可控的，为后期设计不同性质的物体提供了重要的设计参考。

5 多级仿生微结构设计案例

人们对牙釉质和牙本质的弹性性能、硬度和断裂力学性能进行了大量的研究，其中弹性模量[12]和硬度[13]等参数是目前研究最为广泛的力学特性。由于不同牙齿的力学性能有较大差别[14]，而第三磨牙的力学性能最为丰富，所以选取成年男性的第三磨牙作为研究对象。主要针对天然牙釉质、牙本质以及牙髓构成的多级结构弹性模量作为主要的调控目标，能够实现多级晶格结构的优势：外硬内软，可以凸显牙齿坚韧的特性。

(a)体重建(b)三维模型剖视图

(c)晶格填充(d)有限元模拟力与约束

(e)有限元变形对比(f)有限元应力云图

(g)结果(h)结果验证图10 牙齿仿生多级模型制作过程Fig.10 The primary processes of the bionic multilevel model design of tooth

主要通过以下步骤实现多级牙齿模型的设计(图10)。首先根据牙齿μCT切片重建出牙釉质层、牙本质层以及牙髓腔层三个独立的层级空间，并将其分别表示为Ω1、Ω2和Ω3。然后通过设定各级空间的力学参数范围(表2)，在晶格模型族中间按照晶格选取原则选取三种晶格，并对各级空间进行平铺填充，如图10c所示。随后固定住多级牙齿模型的根部，模拟牙周韧带和牙槽骨将自然牙约束的情况，施加垂直咬合载荷，见图10d。通过有限元分析，获得多级牙齿应力场云图，见图10f。由于在牙釉质层的右上方应力较集中，因此应降低晶格孔隙率，提升其硬度，增强其抗磨切能力；在各层级连接处也存在一定的应力集中现象，所以过渡的晶格杆径较大，以提高承受应力集中的能力。将应力场云图对多级晶格结构进行映射，获得梯度变化的多级晶格模型，并将模型各级间的连接进行增强，以保证力的有效传递。随后将力在50～450 N之间变化，以判断模型变形区间是否在最大变形范围内，并及时调整晶格单元的孔隙率，间接驱动各级晶格弹性模量的变化，确保满足实际的力学需求。整个模型在图10c～图10h之间迭代了了五次，为保证美观，将最外层超出原始边界的晶格单元进行裁剪，最终获得图10g的结果，为保证设计结果的有效性，对生成的模型进行验证，如图10h所示。结果表明生成的模型最大微动尺寸小于且接近自然牙的最大微动尺寸45.6 μm[15]，满足设计要求。

本文提出的方法可以解决种植牙与牙周骨的刚性连接问题，生成的牙齿模型不仅具有坚硬的外表层，而且拥有能够缓冲变形的内部结构，实现牙齿坚韧的多级结构特性，且多孔的结构非常利于骨组织长入，随着时间的推移，种植体与牙周骨间的连接会更加紧密，不会像传统的种植牙出现松动的问题。

表2 第三磨牙牙釉质与牙本质弹性模量测量结果[14]

6 结论

(1)本文提出了多级晶格仿生微结构建模的概念和几何设计方法，实现了多级晶格结构的力学性能调控及多种晶格连接合成。

(2)通过预构建具有不同弹性模量的微单元结构模型族，提高了建模的效率，通过有限元映射、梯度变化、隐式曲面融合等方法获得了特定的弹性模量及抗冲击韧性。

(3)将该技术应用到生物医学假体牙齿的建模之中，突出了多级晶格结构的功能性优势，生成的模型不仅拥有坚硬的外部结构以抵抗破坏，又拥有坚韧的内在结构以缓解冲击变形，实现了结构坚韧的生物力学特性。

(4)多孔结构可以有效减小模型的质量，实现轻量化，且能够更好地保证与周边的骨组织相互连接，利于生长，使其生物相容性更强。