骨小梁必須能夠承受日常行為和受傷時引起的載荷。由於骨小梁的高度多孔性和複雜結構，並且這種多孔性和複雜結構在不同的解剖部位和不同的個人之間，差異極大，因而研究骨小梁的機械特性非常具有挑戰性。雖然細觀有限元分析 (μFE) 是分析骨小梁機械特性的最常用的方法，但由於這些模型的尺寸很大，迫使研究人員使用自定義程式碼和線性分析方法。ABAQUS 的非線性功能可以對這些模型進行有效分析，為重要的研究課題提供答案。

ABAQUS 的主要功能和優點

建立模擬骨骼組織機械特性的本構模型

輕鬆求解幾何的和材料的非線性模型

支援並行求解

背景知識

　　骨小梁位於長骨（如股骨）的末端和立方骨（如脊骨）中，是人類骨骼中承載生物組織的主要生物組織。它的機械特性有很高的臨床價值和研究價值。增進對骨小梁機械特性的瞭解，有助於深入研究骨骼的斷裂機理，也有助於評估年齡、疾病和藥物治療的影響。骨小梁是一個充滿孔洞的組織——脊骨中 85% 以上是孔洞，並有著複雜的結構，而且這兩者都因不同的人和不同的解剖部位而存在很大差異（參見圖 1）。因此，要用統計學的方法確定骨小梁的機械特性，就需要多個樣本的機械特性資料。

圖1 人類脊骨圓柱型小梁樣本的透檢視

　　細觀有限元法 (μFE) 廣泛應用於骨小梁機械特性的研究，包括在光譜水平和微結構水平兩個方面的研究。這些模型是透過對骨小梁樣本進行高解析度成像得到的，樣本被自動劃分成元素為六面體的有限元網格（參見圖 2）。網格中所有的單元都完全一樣，一般尺寸在 50 微米。劃分網格後，一個邊長 5 毫米的立方體樣本的 μFE 模型一般具有 50 萬個自由度。與試驗用樣本（8 毫米直徑和 15 毫米長）類似的骨骼樣本的 μFE 模型則有幾百萬個自由度。

圖 2：含有 44μm 單元的骨骼樣本中一塊邊長2.5 毫米立方體的 μFE 網格圖

　　在過去，這些大量的問題會使許多研究人員不得不利用自定義程式碼，一個單元一個單元地迭代求解。由於非線性有限元模型非常複雜，這些自定義程式碼只限於線彈性分析。雖然線彈性有限元模型不能模擬骨骼受損情況，但是研究人員經常利用它與試驗資料校核，確定骨骼組織的彈性特性。然而，關於骨小梁非線性機械特性的許多問題還有待解決.因為 ABAQUS/Standard 能夠利用並行處理能力解決大型問題，包括複雜材料模型問題，所以它非常適合這類分析。在本技術簡報中，我們利用 ABAQUS/Standard 研究了幾何非線性在骨小梁機械特性中的作用。我們對一個具有四百多萬個自由度的模型進行線彈性分析。透過檢驗此分析的並行處理能力（也就是可伸縮性），我們展示了它求解大型問題的可行性。

有限元分析方法

　　利用顯微X線斷層攝影技術 (μCT 20,Scanco Medical AG, Bassersdorf, Switzerland)，以 22 微米的解析度，對容積率為 9% 的人類脊椎骨小梁樣本進行成像（參見圖 1）。建立了兩個 μFE 模型。首先，整個圓柱形樣本被劃分為大小44 微米的六面體單元網格（參見圖2）。然後，從圓柱體中心劃出一個邊長為5毫米的立體子區，建立另一個具有相同單元大小的模型。兩個模型的網格數量見表 1。

表 1：μFE 模型的網格數量

模型

單元數量

節點數

自由度數

圓柱體

828,853

1,380,834

4,142,502

立方體

131,322

216,027

648,081

　　圓柱體模型被用來評定直接稀疏求解器的並行處理能力。在沒有摩擦的情況下，在頂面和底面應用位移邊界條件，模擬 1% 壓縮應變。分別用 HP rx8620 計算機中 1、2和4 個 CPU 進行線彈性分析。

　　利用邊長為 5 毫米的立方體模型進行非線性分析。這個尺寸的立方體已大得足夠確定平面特性，同時又小得足夠確保非線性分析的可行性。骨骼組織模型是用鑄鐵塑性材料製造的。鑄鐵塑性材料在受到拉伸和壓縮時，其彈塑性狀態會有不同的屈服強度和硬化，因此會產生一個非對稱的單元剛度矩陣。因此，需要使用非對稱儲存的並行稀疏直接求解器。樣本是一個彈性模量為 13.4Gpa，泊松比為 0.3 的組織。（參見參考文獻 2）根據人類股骨骨小梁組織的屈服應變（參見參考文獻 3），鑄鐵塑性模型組織拉伸時的屈服應力為 55.2MPa，壓縮時為 110.6MPa。在拉伸和壓縮時，使用的是相當於彈性模量 5% 的硬化斜度。在無摩擦的位移邊界條件下，拉伸和壓縮採用了 2% 的公稱應變。在這樣低的公稱應變條件下，骨骼微結構的自力接觸可以忽略。此外，每個模型都進行了考慮和不考慮幾何非線性變形的模擬。總共進行了四個非線性分析，為了進行比較，還計算了平面屈服應變。所有對立方體的分析都是在一臺 IBM Power4 計算機上進行的，使用了兩個CPU。

結果和結論

　　利用 4 個 CPU 對圓柱體模型進行線性分析，用時不到 16 分鐘，佔用記憶體不到 11 GB（參見表 2）。表 2 還包括了平行直接求解器的計數結果；加速因數是根據求解時間得到的。對具有幾何非線性的立方體 μFE 模型進行非線性分析，用時不到 7.4 小時，佔用記憶體 4.1 GB。每個非線性分析需要大約 100 個線性方程的解，這就強調了求解器可伸縮性的重要性。骨骼結構中初始屈服點的定位使得非線性分析的收斂變得更加具有挑戰性（參見圖 3）。

表 2：6.4-3 版直接稀疏求解器效能

CPU 數求解時間（秒） 加速 總時間（秒） 1 5541.00 1348 2 2951.88 1107 4 1713.24 945

圖 3：在 2% 壓縮應變情況下邊長為 2.5 毫米立方體的μFE 模型的骨骼結構區域性應力分佈圖

圖 4 是根據表觀應變（樣本長度的變化/原始樣本長度）作出的表觀應力（外力/橫截面面積（25 平方毫米））圖。初始屈服點定義為偏移量達到 0.2% 的點。與試驗資料（參見參考文獻 4）類似，壓縮時的屈服應變比拉伸時的大。

表 3：不同模型組合的屈服應變

幾何非線性拉伸 壓縮 考慮 0.61 0.78 不考慮 0.590.86

圖 4：四個非線性分析的應力應變關係圖幾何非線性在壓縮時引起軟化，拉伸時引起硬化標記顯示由 0.2% 偏移量方法（點劃線）確定的初始屈服點

　　雖然組織材料在硬化，但當考慮幾何非線性時，很明顯地觀察到有軟化發生（參見圖 4）。另外屈服應變與試驗測量結果類似，特別是壓縮的情況下（參見參考文獻 4）。這些結果表明在拉伸和壓縮時骨小梁組織有著不同的屈服特性，應該把 μFE 模型和幾何非線性結合起來，精確地模擬骨小梁的平面屈服特性。

致謝

ABAQUS, Inc. 衷心感謝加州大學伯克利分校的 Tony M. Keaveny 教授提供骨骼樣本成像資料和有限元網格。

參考文獻 van Rietbergen, B.; H. Weinans; R. Huiskes; A. Odgaard, “A New Method to Determine the Trabecular Bone Elastic Properties and Loading Using Micromechanical Finite Element Models,” Journal of Biomechanics, vol. 28, pp. 69–81, 1995.Rho, J. Y.; T. Y. Tsui; G. M. Pharr, “Elastic Properties of Human Cortical and Trabecular Lamellar Bone Measured by Nanoindentation,” Biomaterials, vol. 18, pp. 1325–1330, 1997.Bayraktar, H. H.; E. F. Morgan; G. L. Niebur; G. E. Morris; E. K. Wong; T. M. Keaveny, “Comparison of the Elastic and Yield Properties of Human Femoral Trabecular and Cortical Bone Tissue,” Journal of Biomechanics, vol. 37, pp. 27–35, 2004.Morgan, E. F.; and T. M. Keaveny, “Dependence of Yield Strain of Human Trabecular Bone on Anatomic Site,” Journal of Biomechanics, vol. 34, pp. 569–577, 2001.