基于應力分布的3D打印殼模型構造與優化①

徐洪飛,黃章進,宋 鵬,2

1(中國科學技術大學 計算機科學與技術學院,合肥 233027)

2(洛桑聯邦理工學院 計算機通信學院,洛桑 瑞士 1015)

從上個世紀80年代第一家3D打印設備公司誕生到現在,隨著材料與工藝技術的不斷進步,3D打印技術愈發成熟并逐漸流行起來.方便個人定制、適用于任意形狀、快速成型等突出優點使3D打印有力的推動了個性化產品的設計制造[1]、機械結構設計[2]、三維重建[3]等領域的發展.

盡管3D打印機器越來越普及,打印材料的價格依然遠高于傳統制造方式所用的原料,如常用的光敏樹脂材料價格高達2000元/千克.這一定程度上制約了3D打印技術的發展和推廣.因此,如何有效降低模型的實際打印體積,節約打印成本是3D打印技術面臨的一個重要挑戰[4,5].針對節材打印這一問題,國內外研究者提出了許多解決方案.Wang等[6]受到建筑中的桁架結構啟發提出一種“蒙皮-鋼架”的輕質結構,并基于該結構進行模型的優化.這種方法能在有效降低打印成本同時,保證打印結果滿足各種約束條件,如物理強度、平衡性等.Lu等[7]受到自然界中的一種介觀結構——蜂窩結構的啟發,提出一種內部掏空優化方法.該方案基于應力分布來對模型空間進行自適應Voronoi分割得到大小不同的Voronoi單元,通過掏空Voronoi單元構造出類似于蜂窩的內部結構.其優化結果也能夠一定程度上降低打印體積,同時有著很好的物理強度.Zhang等[8]受到物體的中軸和骨骼結構的啟發,提出一種中軸樹結構作為物體內部結構型.這種結構型主要由中軸結構、邊界框架以及一組連接桿件這三部分組成.徐文鵬等[9]借鑒傳統漸進結構優化方法,給出一種拓撲優化算法.Zhang與徐文鵬的方法都很好地反應了模型受力的傳遞路線并且大幅節省了材料的使用.

另外,考慮到打印出的模型在現實場景下可能受到的各種外力,在節省材料的同時,也要保證模型結構的力學穩定性,使得打印結果滿足一定的強度要求.這一方面可以通過采用具有更好力學性能的材料進行打印,如具有點陣結構的新型復合材料[10,11];另一方面則可以對模型做進行幾何結構優化.本文主要討論后者.在打印模型之前對模型進行準確的力學穩定性判斷,并對存在強度問題的部位進行合理的結構優化可以歸結為模型的強度優化問題[4,5].對此類問題 Stava等[12]給出一種解決方案,該方案通過自動檢測薄弱區域并對模型進行內部挖洞、局部加厚以及加支撐這三種方式來提高模型的強度,調整后的模型可以承受用戶指定的載荷.Zhou等[13]依據模型的幾何特征和材料屬性直接計算出模型的最不利載荷,同時識別出模型上最容易被破壞的區域.Zhao等[14]提出了一種使用厚度參數來構造滿足給定力學穩定性約束的殼模型的方法.該方法依據模型的邊界面和用戶指定的外力,根據應力約束來優化厚度參數以擠出內外表面,最終構造出滿足約束的殼模型.

上述研究成果在節材打印、強度優化方面取得了很好的效果,為本文提供了寶貴的經驗和參考.同時,我們注意到,節材打印的相關研究[6-8]都采用了一些特定的結構型式來生成內部稀疏的空間結構,以此來減小打印體積.這樣帶來的一個問題是模型中會存在著大量的懸空結構,如Lu等[7]的蜂窩結構、Zhang等[8]的中軸樹結構.這種模型需要使用無需支撐的打印技術如3DP、SLS或者使用可溶解的支撐材料才能打印完成,一定程度上提高了打印的硬件門檻和操作難度.而關于強度優化的相關研究[12-14]均是通過對模型做變形來解決的,雖然能有效的提高模型強度,但會導致模型外觀的變化,影響使用體驗.

針對當前研究中存在的上述兩個問題,本文提出一種基于應力分布的殼模型構造優化方案.本方案首先構造出均勻厚度的初始殼模型,然后進行基于應力分布的殼模型優化迭代,得到具有極小體積的殼模型.經過多組實驗驗證,優化后模型的體積為原模型的17.2%～24.4%,節材效果顯著.此外,本文方案生成的殼模型還具有如下優點:不改變模型的外觀;使用普通的桌面級打印設備即可打印,如最常見的FDM打印機;添加的支撐材料方便去除.

1 優化方案概述

給定一個待打印的實心三維模型M,使用打印材料m,指定載荷F,屈服強度σy,最小打印厚度tmin.本文所要解決的問題是如何通過優化M的幾何形狀來減小模型的打印體積.同時,需要滿足如下約束條件,模型的外觀不變、在載荷的作用下不發生屈服,以及模型的最小厚度tm不低于最小打印厚度.所以優化問題可以描述為:

其中,第一個約束條件為屈服強度約束條件,用來保證模型可以承載指定的載荷而不發生破壞.依據第四強度理論設置,即當構件各部位的von Mises等效應力(見第2節)的最大值σmax達到材料的屈服強度σy時,材料就會開始屈服,發生塑性形變,外力撤銷后無法恢復原來形狀.

針對上述優化問題,本文給出一種基于應力分布的殼模型構造優化方案,概述如下:首先構造與M對應的均勻厚度空殼模型M0,對M0進行應力分析,如果不滿足屈服強度約束,則均勻增加M0的厚度,直到滿足約束為止;然后基于M0的應力分布情況,自適應地調整Mi各位置的厚度得到Mi+1,直到迭代結果Mk達到某個約束條件時停止,Mk-1即是具有極小打印體積的優化后的殼模型.由于模型外觀不變,即模型的外表面R不變 (注:外表面R圍成實心模型M),因此優化殼的厚度也就等價于優化內表面S的形狀.

2 應力分析

應力描述了連續介質內部各部分之間相互作用力的強度,在連續介質力學中將其定義為單位面積所承受的作用力.本文的殼模型構造優化方案一方面借助應力分析結果指導迭代過程,調整模型厚度,另一方面也使用應力分布信息驗證優化結果,判斷是否還可以繼續優化.

由于絕大多數桌面級3D打印機僅支持塑料(PLA、ABS、尼龍等)作為打印材料,因此本文只考慮塑性材料作為打印材料的情況.對塑性材料的屈服情況的判斷,一般使用材料力學的第四強度理論[15].該理論認為引起材料屈服破壞的主要因素是畸變能密度,其強度條件見公式 (2),其中σ1、σ2、σ3為主應力,σy為材料的屈服強度.

von Mises是一種屈服準則,表示在一定的變形條件下,當受力物體內一點的等效應力達到某一值時,該點就開始進入塑性狀態,產生屈服.屈服準則的值即為von Mises等效應力,計算公式如下:

由式(2)、式(3)可知,von Mises等效應力遵循第四強度理論,所以本文使用von Mises等效應力來度量模型各個部位的應力水平,并通過與材料的屈服強度進行比較來判斷模型是否發生屈服.

為了給構件必要的強度儲備,材料力學中引入了許用應力[15]的概念,即將材料的屈服強度除以一個大于1的安全因數n,作為構件應力不可超過的值,通常對于塑性材料n取1.1～2,此時屈服強度條件約束為:

同理,本文取不大于1的安全系數s1作為屈服強度的系數,則引入安全系數的屈服強度約束條件為:

本文中使用PLA材質時將s1設置為0.7,使用PC-ABS材質時將s1設置為1(為了與文獻[7]中Lu方案進行對比).若想給模型更高的強度儲備,可以將s1設置為更小的值(一般不小于0.5).此外,為了消除打印機的打印誤差帶來的影響,本文還引入了最小打印厚度對應的安全系數s2,實驗中將s2設為1.2.至此,引入了安全系數的優化問題表示為:

為了實現對模型應力分布的高效準確計算,本文借助機械設計軟件SolidWorks的Simulation工具包進行基于有限元的靜力分析.依次進行模型的四面體化,物理環境的配置(包括所受外力F,材料屬性m,固定位置等),以及基于有限元的靜力分析.最后將四面體化的結果和應力分布結果導出,作為模型優化程序的輸入.圖1展示了此應力分析步驟的工作流程.需要說明的是,雖然此工具能夠得到精確程度較高的應力分布結果,但計算效率不高,本文中所有模型的應力分析時間在 3～10 min/次.

3 殼模型的構造與優化

3.1 均勻厚度殼模型的構造

利用三維建模軟件創建的待打印模型往往只有一個封閉的外表面R,該外表面一般由一組三角面片連接組成.因此如果不對模型進一步處理,直接打印會得到一個實心的模型,材料耗費會非常大.所以一般在打印前需要對模型進行空殼處理,即計算出一層內表面用來與外表面共同構成封閉的空殼模型.對于凸多面體,可以直接將外表面頂點沿著法線方向向內偏置一段距離來得到內表面.但對于一般非凸多面體,簡單的向內偏置外表面會導致三角面片的自交.Jung等[16]提出的自交檢查和去除方法可以消除三角面片自交的情況,但殼模型卻不再厚度可控.為此,本文采用了基于模型體素化表示的等值面提取算法來構造內表面(見算法1).

圖1 應力分析流程

算法1.基于體素化的均勻厚度殼模型構造算法輸入:模型外表面R,指定厚度為t.輸出:模型內表面S (與R圍成均勻厚度殼模型(S,R)).1.以a為邊長,體素化R所在包圍盒,得到一組體素,記為Voxel.2.FOR每個在Voxel中的體素vox.3.判斷vox相對R的位置(分為內部,外部,相交).4.以R為邊界,Voxel中體素為離散元素構建有向距離場Dc.5.以0.25a為邊長,構建細粒度有向距離場Df.6.取Dc中具有距離值t/a的一層體素Voxelc.7.FOR每個在Voxelc中的體素voxc.8.獲取8個頂點在Df中的距離.9.對Voxelc應用Marching Cubes算法[17]提取內表面S.10.RETURN S.

圖2展示了算法1的相關流程.此算法不但可以避免三角面片自交,還能夠提取出過渡非常平滑的內表面,見圖2(c).需要說明的是,在算法1第2、3行的實現上,本文采用了由粗到細的體素化方式進行加速,大幅提升了算法效率.為了得到足夠精細的內表面,需要保證體素尺寸足夠小,因此就需要處理數目龐大的體素,以尺寸為10×20×15 cm3的Kitten模型為例,若體素邊長為1 mm,則需要計算 3e9個體素的位置,這涉及到大量的相交檢測判斷.在 Intel Core i5-4590@ 3.30 GHz 平臺下,直接的對細粒度的體素進行位置判斷的平均耗時為1.9分鐘.為了減少時間消耗,我們采用由粗到細的體素化方式進行加速.具體做法是,先進行粗粒度的體素化,體素邊長設為 4 mm.對這些大體素進行位置判斷,刪除在模型外部的大體素,由模型內部的大體素細分出的小體素(邊長1 mm)一定在模型內部,所以只需對那些由模型邊界上的大體素細分出的小體素進行位置判斷,從而大大減少了相交檢測判斷次數.且所有在模型內部或邊界上的邊長1 mm的小體素的相對位置均被正確判斷,不會影響到有效信息的完整性.使用這個新的方式對該模型進行處理,平均耗時縮短為15秒,為優化前的13.2%,效率提升明顯.

3.2 基于應力分布的殼模型構造與優化

為了保證模型打印完成之后具有一定的力學穩定性(即可以承受給定的載荷),需要模型上各單元的von Mises等效應力不超過材料的屈服強度.所以當模型簡單處理為均勻厚度的殼模型時,為保證所有單元的應力值在安全范圍內,對于那些本身應力值較低的區域需要和應力值高的區域具有相同厚度.這會導致材料的浪費.這一問題在應力分布不均的模型身上會突顯出來.故而對于算法1構造出的均勻厚度模型,需要進一步去調整內表面形狀,最大限度的減小模型體積.

本文提出一種基于應力分布自適應調整內表面形狀的方法.具體做法是讓頂點沿著法線方向外擴,應力大的頂點具有較小外擴速度,應力小的頂點具有較大外擴速度.圖3給出了基于此方法進行一次內表面頂點外擴的示意.其中頂點a,b具有較小應力,頂點c,d具有較大應力,外擴后的形狀如圖3(c)所示.

圖2 均勻厚度殼模型的構造流程

圖3 基于應力外擴頂點

基于以上思路,我們給出每次迭代時頂點v的外擴距離計算公式,見式(7).其中,do為頂點v到模型外表面的距離,σv為v的應力值,s1、s2分別為屈服強度、最小打印厚度的安全系數(見第2節).公式第一項用來保證最小打印厚度約束,第二項用來計算頂點在當前應力情況下的外擴距離,其中p用來控制每次外擴的最大距離.

基于式(7),我們提出基于應力分布的殼模型構造與優化算法,見算法2.

算法2.基于應力分布的殼模型構造與優化算法輸入:外表面R,屈服強度σy,最小打印厚度tmin.輸出:殼模型(S,R),S為構造的內表面.1.對R進行體素化并構建距離場Dc、Df.2.利用Marching Cubes算法按厚度t提取S0.3.WHILE (S0,R)最大應力 σmax＞s1σy.4.t←t+t0.5.利用Marching Cubes算法按厚度t提取S0.6.WHILE (Si,R)最大應力 σmax≤s1σy.

7.FOR 每個Si上的頂點v.8.利用式(7)計算外擴距離dv.9.IF dv==0.10.BREAK.11.FOR 每個Si上的頂點v 12.沿頂點法線方向向外移動dv距離.13.對Si應用算法1按厚度0提取Si+1.14.i←i+1.images/BZ_111_1362_1506_1382_1531.png15.IFdv==0.16.S←Si.17.ELSE 18.S←Si-1.19.RETURN (S,R).images/BZ_111_1378_1812_1399_1837.png

算法2分為兩個階段:前5行進行均勻厚度殼模型的迭代構造;后面則是進行基于應力的殼模型的迭代優化.此算法可以構造出具有極小體積的非均勻厚度殼模型,從而大幅減少打印材料的成本.

注意這里得到的不是具有全局最小體積值的模型.因為在每次調整內表面時,各頂點移動方向均是向模型外部,模型體積單調減小,而載荷保持不變,各單元的應力會逐漸增大,直到某一單元的應力達到材料的屈服強度時為止.若此時允許那些達到屈服強度約束的頂點向內移動,使得其應力回到屈服強度約束之下,則可以固定這些頂點,繼續向外擴展其他頂點,這會進一步減少模型體積.然而這種增加頂點移動自由度的做法會極大地增加算法的復雜程度,考慮到每次迭代過程中應力分析耗時較大的特點,如此進行下去會非常耗時,可行性不強.因此,本文優化算法不允許頂點向內移動,而是讓應力較小的頂點具有較大的外擴速度,好處是可以很快收斂生成具有極小體積的殼模型,保證方案的可用性.

3.3 算法實例分析

這里以Kitten模型為例,展示算法2的迭代過程與優化效果.初始Kitten模型包圍盒尺寸為6.5×10×6 cm3,最小打印厚度為1 mm,施加載荷為 1000 N,這里使用的打印材料是PLA,其相關物理屬性見表1.另外,式(4)中相關系數為:s1=0.7,s2=1.2,p=0.002.

表1 打印材料屬性參數

首先構造均勻厚度殼模型,初始厚度為5 mm.然后對其進行應力分析(見圖4).上方箭頭表示載荷,底部箭頭表示模型的固定位置.此時最大應力為22.7 MPa,小于屈服強度約束s1σy(=49 MPa).

圖4 均勻厚度殼模型的von Mises 應力分布

接著對初始殼模型進行優化迭代,得到最終結果.圖5展示模型優化過程中剖面的變化.結合圖4中的應力分布信息,可以發現優化后模型在初始應力較高的區域(方框標識),厚度較大,在應力較低的區域(圓圈標識),厚度較小.圖6展示了優化過程中體積、應力的變化,可見隨著模型體積的降低,模型的平均應力在逐漸增加,直至模型的最大應力超過屈服強度約束.經計算,優化后Kitten模型的體積為輸入模型的22.4%,節材效果顯著.

圖5 模型Kitten的優化過程

4 實驗與討論

本文算法由一個待打印的實心3D網格模型開始,對其進行內殼模型的構造,并基于應力分布進行內殼優化,最終得到的內表面與輸入的外表面圍成體積極小的殼模型.本文的實驗環境為Windows操作系統,處理器為Intel Core i5-4590 @ 3.30 GHz,內存為8 GB.本文算法2的殼模型優化部分需要迭代5～9次得到最終結果,而每次迭代時應力分析過程是影響算法速度的瓶頸,一般在3～10分鐘/次.本文進行的所有實驗中優化時間最多為1.5小時,平均為35分鐘.為了充分驗證算法的有效性,我們在對多個經典模型進行優化.圖7和表2展示了優化結果的信息.從表2中模型的優化數據來看,使用本文算法優化后殼模型的體積僅為輸入模型的17.2%～24.4%,大幅降低了材料的用量.

圖7 初始殼模型的應力分布與最終優化結果

表2 模型的優化數據

我們與同樣基于應力分布的Lu等在文獻[7]中提出的節材優化方案做了3組對比.兩個方案均使用PCABS作為打印材料,屈服強度約束取41 MPa.對比情況見表3和圖8.可以發現,Lu等方案的優化結果雖然具有較小的最大應力,但模型體積為原來的26.0%～43.6%,遠高于本文的19.8%～24.4%.這主要因為Lu等方案的優化結果受到蜂窩結構這種特定內部結構型式的限制,無法進一步減小模型體積.

表3 與Lu等[7]方案的優化結果對比

圖8 Lu等[7]方案結果(左)與本文結果(右)

我們使用基于FDM打印技術的3DTALK MINI打印機打印制造了3組PLA材質的優化模型及剖面(見圖9).模型高度均為9～10 cm.打印優化模型需要3～3.6小時,打印剖面模型需要1.4～2小時.在模型底部預留小孔,方便打印完成之后使用工具取出打印支撐.為了驗證優化后模型的力學穩定性,本文使用英斯特朗電子動靜態疲勞試驗機對打印出的模型進行壓力測試,見圖10.以Bimba模型為例,觀察機器數據采樣系統給出的折線圖10(b)可以發現,隨著壓力的增加,模型形變大小呈線性增長,并在載荷達到1000 N時,形變為1.2 mm.這說明模型處在彈性形變狀態,沒有發生屈服,表明優化結果模型可以承受指定載荷而不被破壞.

圖9 優化模型的打印結果

圖10 模型受力測試

5 結論與展望

本文針對3D打印領域中的材料節省問題提出了一種基于應力分布的殼模型構造優化算法.該算法能夠生成體積極小的非均勻厚度殼模型.在初始殼模型的構造過程中,先將模型體素化表示,采用由粗到細的體素化方式加速體素相對位置的判斷,接著構建距離場,利用Marching Cubes算法提取初始內殼;在內殼的優化過程中,基于應力分布的內殼外擴過程可以自適應降低模型各位置厚度,最終得到具有極小體積的殼模型.

本文實驗模型的優化結果體積為原來體積的17.2%～24.4%,顯著降低了材料消耗.此外,本文方法還具有不改變模型外觀、內部支撐容易去除、適用于多種類型打印機的突出優點.最后,我們通過壓力測試驗證了模型的力學穩定性,進一步驗證了本文算法的有效性.

本文工作還存在很多改進之處.首先,基于有限元方法的應力分析方法雖然計算準確,但效率偏低,影響了算法的執行速度.因此,提高有限元分析的計算效率、減少應力分析所需的時間是讓本文方法能夠更加高效應用的關鍵.其次,本文方法僅支持能夠在打印機中一次打印完成的模型,對于那些超出打印機尺寸的模型,需要先進行模型的分割[18,19],然后在進行體積優化.這涉及到多個部件模型組裝后的應力分析以及協同優化,很具有挑戰性,也是我們下一步要研究的主要內容.再次,本文的約束條件包括屈服應力、最小打印厚度,以及外表面保持不變,而在實際應用中,可能會有其他的要求,比如模型的平衡性,即如何讓模型以某個姿態穩定放置或穩定運動[20,21].因此,下一步考慮從這個角度對本文優化方案進行擴展.