體動力學網格劃分技術在骨骼有限元建模中的應用

體動力學網格劃分技術在骨骼有限元建模中的應用

FEA有道 FEA有道生物力學 前天

點擊下方「藍字」關注我們吧

鄧達人1，孟春玲1，馮敏山2，張剛1

(1．北京工商大學材料與機械工程學院，北京100048;

2．中國中醫研究院望京醫院，中醫正骨技術北京市重點實驗室，北京100102)

摘要:

目的：研究一種高質量、高效率的骨骼分網方式。

方法：結合骨結構的特點，將流體動力學(computational fluid dynamics，CFD)的分網方法引入到骨生物力學的網格劃分中，用六面體單元模擬皮質骨，四面體單元模擬松質骨。

結果：採用CFD網格劃分技術能夠劃分出高質量的六面體單元，並且較好地模擬骨骼的結構特徵，可以實現計算機自動分網，分網花費時間約為傳統分網方法的1/5，有限元模型計算結果基本符合屍體實驗測量值。

結論：CFD網格劃分技術可應用於骨生物力學領域中，為複雜人體骨骼的重構探索了一條有效途徑。

關鍵詞:骨生物力學;皮質骨;有限元模型;流體動力學;網格劃分

在人體生物力學分析中，由於人體的解剖結構十分複雜，且構成人體的三維實體表面往往是不規則的自由曲面，故如何建立合理有效的有限元模型一直是一項非常重要而煩瑣的工作［1］。網格單元質量的好壞將直接決定有限元的計算結果，網格質量可能導致應力計算結果偏差達到7%～100%［2］，低質量的網格甚至會導致有限元計算無法完成。

自1972年現代生物力學之父馮元楨首次將有限元法應用到生物力學研究以來，各國研究者先後採用了多種網格劃分方法，其中主要有四面體網格自動劃分法和六面體映射網格劃分法［3-5］。田衛軍等［6］基於逆向技術製作骨盆環四面體有限元模型;曹立波等［7］建立及驗證兒童假人頭部有限元六面體模型。作為最簡單的三維單元，四面體單元是常應變、常應力單元，計算精度較低，而且對於複雜的接觸條件常伴隨計算不收斂，無法滿足複雜載荷條件下的力學分析;當三維實體表面是十分複雜的自由曲面時，若全部採用六面體單元，人工分區比較困難，不得不對外觀做較多簡化，網格劃分耗時也較長，甚至在曲率變化複雜的畸形表面難以實現。

骨結構分為松質骨和皮質骨，皮質骨硬度和密度均很高，是骨的主要受力部分。同時，骨骼之間的非線性接觸問題一直是骨生物力學研究的重點，皮質骨表面應力梯度較大，故皮質骨需要採用精度比較高的網格。本研究結合骨結構的特點，將流體動力學(computer fluid dynamics，CFD)流體邊界層的網格劃分技術［8］，引入到骨生物力學的骨骼網格劃分中。把皮質骨的外表面劃分為二維四邊形網格，以該網格做參考網格向模型內部映射出若干層三維六面體單元，用六面體單元模擬皮質骨;對於受力較小的松質骨採用四面體單元。二維四邊形網格劃分曲面時不用考慮人工分區和網格映射，受到曲面曲率變化的影響較小，無需對模型外觀做太多簡化［9］;CFD網格劃分法吸收四面體單元分網簡單快速和六面體單元計算精度高的特點。研究工作為複雜人體骨骼的重構探索了一條有效途徑，為人體骨骼的生物力學分析奠定基礎。

1 骨模型傳統網格劃分方法

1． 1 四面體網格自動劃分法

四面體網格劃分方法是通過在所選區域內進行布點，然後將它們聯接成四面體，該方法是目前最流行的實體網格劃分方法之一，幾乎在所有的有限元前處理軟體中都可以實現該方法的網格劃分［10］。但是在網格三角化過程中，使用該方法對所生成的單元形狀難以控制，網格計算精度相對較低。如果通過增加網格密度提高精度又會導致計算量大大增加，對研究骨骼之間的接觸等複雜問題難以完成計算［11］。例如:田衛軍等［6］基於逆向技術的骨盆環三維重構與結構分析，其有限元模型採用四面體單元，椎體與椎間盤、恥骨聯合處以及小關節之間採用GLUE連接，將各個部分的模型粘連成一個整體，無法研究骨骼之間的接觸問題［12］。

1． 2 六面體映射網格劃分法

六面體網格劃分一般採用映射網格劃分方法。該方法通過適當的映射函數將待剖分物理域映射到參數空間中形成規則參數區域，然後對規則參數區域進行網格剖分，最後將參數域的網格反向映射回物理空間，從而得到物理域的有限元網格［9］。目前常用的網格劃分方法中，六面體映射網格劃分法的算法簡單、計算速度快、單元質量好、單元密度可控，因此在各個領域得到了廣泛的應用。

骨骼之間的接觸問題一直是骨生物力學研究的重點，由於該接觸問題為非線性問題，皮質骨表面應力梯度較大，故皮質骨需要採用精度比較高的網格;目前而言，六面體單元精度比四面體單元高，因此，皮質骨需要劃分為六面體單元，但是六面體分網無法實現類似於四面體的自動分網，當三維實體的表面曲率變化十分複雜時，六面體映射網格劃分法需要創建大量的輔助線，剖分模型網格映射的區域十分困難，故不得不對模型外觀做較多的簡化，劃分高質量的六面體單元需耗費大量的時間［13］，甚至在曲率變化複雜的畸形表面難以實現。特別在研究畸形骨骼模型時，需要對網格模型進行反覆的更改，導致分網的工作量非常大。

2 CFD 網格劃分法

2． 1 CFD 流體邊界層分網

骨是一種具有獨特結構的高密度結締組織，大體結構分為松質骨(骨松質)和皮質骨(骨密質)，前者分布於骨的內部佔全身骨量20%，密度低於皮質骨且富有彈性;後者位於骨外側，佔全身骨量80%，其硬度和密度均很高，主要用於發揮骨的支撐作用，是骨生物力學研究的重點之一［14］。

本研究根據骨結構組成及特點，將CFD中流體邊界層的網格劃分技術引入到骨生物力學中。CFD網格劃分法是利用大型有限元前處理軟體Hyper Mesh的CFD-Tetramesh功能，將流體管道模型的外表面劃分為二維四邊形網格，以該二維網格做參考網格，向模型內部映射出若干層三維六面體單元作為流體邊界，再將模型內部劃分為四面體單元作為流體域［9］。因此，在對骨骼進行有限元分網時，可以使用CFD網格劃分法，將皮質骨的部分採用精度高、不易分網的六面體單元，將松質骨的部分採用精度較低、容易分網的四面體單元。

2． 2 CFD 網格劃分法的分網過程

以某醫院提供的枕寰樞複合體模型為研究對象，採集CT斷層圖片，應用Mimics醫學影像處理軟體完成枕寰樞複合體的三維重建並獲取點雲數據［見圖1(a)］。

將枕寰樞複合體的點雲數據導入三維建模軟體SolidWorks，進行逆向建模，得出骨骼三維實體模型［見圖1(b)］。

使用大型有限元前處理軟體HyperMesh對骨骼的三維模型進行CFD網格劃分，將骨骼模型的幾何面劃分為二維四邊形網格見圖1(c)］;

以二維四邊形網格為參考，向模型內部映射出3層高質量的六面體單元作為皮質骨［見圖1(d)］;

再將模型其餘部分劃為四面體單元作為松質骨［見圖1(e)］，

完成骨骼模型的網格劃分［見圖1(f)］。

CFD網格劃分法極大減小了分區的難度，保證了模型內部的網格質量;該方法根據模型的幾何面生成二維四邊形網格，保留了原模型的幾何特徵，降低簡化對模型精度的影響，由於二維四邊形網格劃分對象是曲面，劃分曲面時無需考慮網格映射，曲面的複雜程度對二維四邊形網格劃分造成的影響較小，故採用二維四邊形網格劃分幾何面比採用六面體單元劃分幾何實體的難度小;該方法能夠在保持模型幾何面形狀特徵不變的前提下，調整任意位置六面體單元層的厚度和層數，實現模擬皮質骨厚度不均的目的。圖2所示為樞椎網格模型的上關節面剖視圖，該關節麵皮質骨的最大厚度被調整為正常厚度的1．5倍;此外，通過控制外部六面體層厚度、層數和網格密度以及內部四面體的網格生長率，能夠再次調整骨骼模型的網格質量和網格密度。

採用CFD網格劃分法對枕寰樞複合體模型的枕骨、寰椎、樞椎等複雜的三維實體進行網格劃分，建立的最終模型如圖3所示。

圖3枕寰樞複合體的網格模型

該模型保留了骨骼的幾何特徵，骨骼接觸面網格尺寸為0.5mm，其餘部分網格尺寸為1.5mm，皮質骨厚度為1mm，由3層尺寸為0.33mm的六面體單元組成;模型包含103497個六面體單元和617784個四面體單元，單元類型包括C3D20、C3D8、C3D6、C3D4等。模型的單元質量見表1。

3 網格模型驗證

3． 1 有限元模型的建立

將HyperMesh製作的枕寰樞複合體網格模型，通過大型有限元軟體ABAQUS進行韌帶添加、材料參數賦值、接觸定義、設置載荷以及約束條件等，從而完成有限元模型的建立。

韌帶採用非線性彈簧單元模擬，彈簧加載點的選取參考各韌帶的生理解剖部位;將椎骨材料視為各向同性彈性材料，分別把骨性結構、軟骨以及韌帶賦予不同的材料參數，相關材料參數來源於已發表的文獻［15］。根據枕寰樞複合體的生理結構特徵，本模型的接觸發生在以下幾個部位:寰枕關節、寰樞關節、寰齒關節以及齒突與橫韌帶關節處的軟骨之間，且寰枕關節、寰樞關節分別具有兩個側方關節［14］，故本模型共建立12個接觸面、6對接觸對，用於模擬各關節的接觸問題;樞椎椎體下緣採取完全約束，力矩加載點位於枕骨中上部，加載點與枕骨做耦合連接。

3． 2 模型驗證

通過計算枕寰樞複合體有限元模型在4種生理載荷下(前屈、後伸、側屈、軸向旋轉)枕寰、寰樞之間的相對運動角度，並與Panjabi等［16］1．5N·m載荷下屍體試驗數據進行比較，驗證有限元模型的有效性。模型在4種生理載荷下計算的枕寰(C0～1)、寰樞(C1～2)間相對運動角度與屍體試驗結果見表2。

在1．5N·m載荷作用下，側屈運動中寰枕間偏轉角度比屍體實驗數據大10．3°;軸向旋轉中，寰枕間比屍體實驗數據大0．42°，本模型其他運動結果均位於屍體試驗數據範圍內。當1．5N·m外載荷全部施加於模型時，模型整體的最大旋轉角度發生在軸向旋轉運動。

通過比較有限元模型的計算結果與屍體試驗數據，證明採用CFD網格劃分法製作的骨骼有限元模型，具有較高的計算精度，能夠模擬複雜生理運動下的骨骼運動，對上位頸椎的臨床研究具有指導作用。

3． 3 網格收斂性分析

在有限元模擬運算中，不同的網格尺寸會對有限元計算精度產生較大的影響，網格尺寸越小，計算精度越高;但是過小的網格會大大增加模型的單元數目和模型解算時間步，故在建模中需要綜合考慮模擬精度和計算時間成本等要素［10］。

通過比較4組不同網格尺寸的枕寰樞複合體有限元模型在相同生理運動條件下(前屈、後伸)的計算結果，驗證不同網格尺寸對計算結果的影響。有限元模型的網格尺寸見表 3。

將4組有限元模型的計算結果與Panjabi等［16］在1．5N·m載荷條件下的屍體試驗數據及某醫院提供的1．5N·m載荷下屍體試驗數據進行比較，具體數據見表4。1～4號模型的收斂時間分別為27、6、8、3h。

通過表 4 可知:

(1)模型的網格尺寸越小，計算結果越接近實驗值。通過1、2號模型可以發現，隨著網格尺寸的減小，模型計算精度逐漸提高，但是隨著網格的減小其尺寸變化對計算精度的影響逐漸降低。同時，隨著網格尺寸的減小，模型計算的時間成倍增加，大大提高了計算的時間成本。

(2)由1、3、4號模型可以發現，模型計算結果受到接觸面網格尺寸的影響較大。較小的接觸面網格尺寸，能夠顯著提升模型計算精度［12］。

4 結論

根據骨結構的特點，把CFD流體邊界層的網格劃分技術引入到骨生物力學的骨骼網格劃分中。網格模型外觀逼真，較好模擬了骨骼的生理特徵;該分網方法綜合了四面體網格自動劃分法和六面體映射網格劃分方法的優點，降低了分網的難度，提高了模型的網格質量，保證了模型的計算精度;同時受到模型幾何面曲率變化的影響較小，保留了原始模型的幾何特徵，降低了簡化對模型精度的影響;此外，可以調節任意位置邊界層厚度，滿足皮質骨在骨骼不同部位厚度不一致的特點;分網後的模型通過控制外部六面體層厚度、層數和網格密度以及內部四面體的網格生長率，能夠再次調整骨骼模型的網格質量和網格密度。研究工作為複雜人體骨骼的重構探索了一條有效途徑，為人體骨骼的生物力學分析奠定基礎。

FEA有道生物力學服務介紹

BUSINESS PRODUCT

01

承接醫學有限元分析項目委託

盆骨、腰椎、頸椎、肩關節、髖關節、肘關節、膝關節、踝關節、義齒、種植體、上下頜骨、黏膜、牙冠 等分析。

02

醫學有限元分析培訓課程

1. ABAQUS醫學有限元分析及SolidWorks醫學3D建模培訓
2. ANSYS醫學有限元分析及SolidWorks醫學3D建模培訓

一對一教學，隨到隨學！ 可選擇線上或線下培訓。

我們服務過的客戶:（部分）

醫院：

中山大學附屬第一 / 第二 / 第三醫院、中山大學附屬第三醫院、南方醫院、廣州中醫藥大學附屬第一醫院、南部戰區總醫院、北部戰區總醫院、廣東省人民醫院、復旦大學附屬華山醫院、第三軍區醫院、廣東省中醫院珠海醫院、中山大學光華口腔醫學院附屬口腔醫院、廣東省工傷康復醫院、廣州市正骨醫院、佛山市南海區人民醫院、廣東藥科大學附屬第一醫院、雲南省中醫院、山東省臨沂市人民醫院、山東整骨醫院、浙江中醫藥大學附屬第三醫院、海南省人民醫院 等醫院。

科研院校：

南方醫科大學、中山大學、廣州中醫藥大學、廣州醫科大學、廣西醫科大學、復旦大學、同濟大學、華南理工大學、暨南大學、蘭州大學、南昌大學、雲南中醫學院、工業信息化部電子第五研究所、電子五所、核動力四所、中航光電設備研究所、國科軍工、中科華核電技術研究院、等。