多級仿生建模及優化技術

喻長江 戴 寧 李大偉 程筱勝

南京航空航天大學機電學院，南京,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)多孔結構可以有效減小模型的質量，實現輕量化，且能夠更好地保證與周邊的骨組織相互連接，利于生長，使其生物相容性更強。