劉軒廷，孫博華

(1. 西安建筑科技大學理學院，西安 710055；2. 力學技術研究院，西安 710055；3. 西安建筑科技大學土木工程學院，西安 710055)

3D 打印的概念是由Charles Hull 于20 世紀80 年代初提出，問世至今30 余年間，3D 打印技術發(fā)展尤為迅速，并取得了巨大的成功。但當時盛行的立體光刻、熔融沉積建模方式并不適用于建筑行業(yè)，直至KHOSHNEVIS 等[1]開創(chuàng)性地提出了“輪廓工藝”，展現(xiàn)了增材制造技術在建筑領域的潛力。建筑業(yè)的增材制造方法主要是基于擠壓的3D 混凝土打印(3DCP)，其中混凝土材料不斷從噴嘴中擠出，以實現(xiàn)分層打印[2]。3DCP 技術由于材料的局限性、經(jīng)濟性等原因，尚未挖掘出其真正潛力，但是在某些極端條件下，有效地彌補了傳統(tǒng)建筑的弱項，如該技術在防疫方艙醫(yī)院建造方面的應用，于其他星球建造棲息地等[3?8]。

MOHAN 等 [9]針對基于擠出的3DCP 中混凝土的早期材料行為進行了討論，將其分為了兩個階段：1)沉積前，包括將新拌混凝土運送至噴頭，以及新拌混凝土由噴頭擠出時，混凝土打印條的完備性，層截面幾何形狀等[10?14]；2)沉積后，包括3DCP 結構的建立和強度的增加。沉積后，混凝土需具備抵抗自重以及后續(xù)層重量的能力，其次3DCP 結構也需要具有足夠的穩(wěn)定性，因此，用于考慮3DCP 結構塑性坍塌與彈性屈曲的數(shù)值模型被討論[15?22]。此類針對打印進程中力學性能的研究可優(yōu)化打印參數(shù)的設定，同時避免大量試驗。

混凝土沉積后，針對其打印進程中的力學行為，SUIKER[20]首次給出了直壁結構相應的力學模型，在引入混凝土固化性能與打印參數(shù)下，將參數(shù)模型結果與試驗數(shù)據(jù)、有限元模擬進行對比分析，發(fā)現(xiàn)可較好地分析結構在打印進程中的破壞行為，預測直壁結構打印層數(shù)，并為找尋最佳打印參數(shù)集給予了理論指導。WOLFS 等[17]對3DCP圓柱殼結構的早期力學行為進行了數(shù)值建模與試驗測試，但結果顯示，有限元模擬在預測失效變形模式方面是較為準確的，但在定量分析失效層數(shù)方面仍存在提升的空間。WOLFS 等[17]將結果的差異歸因于材料性能的高估以及對幾何、材料缺陷的忽視。OOMS 等[22]提出并討論了一種用于模擬3DCP 結構行為的數(shù)值模型，可用于自動生成數(shù)值分析輸入文件而不管打印幾何的復雜性，但在預測3DCP 圓柱殼失效層數(shù)方面仍存在較大差距。OOMS 等[22]認為結果的差異可能由材料性能的高估與有限元模擬所選的網(wǎng)格、單元類型導致，但其相關方面的研究并未得到較好的結果。因此，針對3DCP 進程中圓柱殼結構的力學性能進行分析是十分重要的[23]。

本文在相關殼體屈曲理論基礎上[24?44]，結合打印進程參數(shù)，構建了用于描述3DCP 圓柱殼彈性屈曲與塑性坍塌的參數(shù)模型，并對二者的競爭關系進行了評判。將得出的理論、模擬結果與WOLFS 等[17]試驗結果進行對比驗證，發(fā)現(xiàn)該文提出的模型可以較好地預測3DCP 圓柱殼的失效高度與失效形式。

1 非均勻內(nèi)力作用下的非均質(zhì)柱殼

在3DCP 進程中，柱殼結構會因混凝土固化行為影響，材料在垂直方向(自重方向)呈非均質(zhì)性。此外，柱殼結構承受的載荷為自重載荷，其應力在垂直方向同樣是不均勻的。與之對應的是，需基于平衡方程與邊界條件，對該非均勻內(nèi)力作用下的非均質(zhì)柱殼的力學模型進行描述。因此，本節(jié)通過柱殼的勢能最小化原理將其壓屈微分方程進行推導。

1.1 勢能

如圖1 所示，柱殼半徑為R，高度為l，均勻厚度h。軸對稱情況下，根據(jù)Goldenveizer-Novozhilov殼體理論，應變-位移關系可以表示為：

圖1 柱殼在自重作用下示意圖Fig. 1 Cylinder under self-weight

式中：x為縱向坐標；y為環(huán)向坐標；εx、εy、εxy為殼體中面線性應變分量；χx、χy、χxy為曲率變化；u、v、w分別為縱向、切向、徑向位移。

柱殼總應變能U，由薄膜應變能Us，與中面拉伸彎曲應變能Ub組成：

其中：

式中：星號“*” 代表材料參數(shù)在空間中是非均勻的(為方便起見，取ν為常數(shù))；E?為非均勻彈性模量，非均勻彎曲剛度由下式給出：

1.2 平衡方程與邊界條件

將式(7)進行勢能最小化處理：

根據(jù)式(8)可以得到非均勻內(nèi)力作用下的非均質(zhì)柱殼的軸對稱壓屈微分方程：

2 3D 打印進程中柱殼的失效機制

如圖2 所示，在3DCP 進程中可能發(fā)生兩種失效機制，分別為彈性屈曲和塑性坍塌。其中當材料硬化速率與打印速率均較慢時，在結構底部容易出現(xiàn)較大變形，從而發(fā)生塑性坍塌。

圖2 3D 打印進程中失效機制示意圖Fig. 2 Failure mechanism of cylinder in 3D concrete printing(3DCP) processes

在本節(jié)中建立了這兩種失效機制的控制方程。彈性屈曲的控制方程可以由第1 節(jié)中所導出的壓屈微分方程和邊界條件所得。在引入打印材料的固化過程和印刷速度相關的時間效應影響，本節(jié)最后給出了3DCP 圓柱殼彈性屈曲與塑性破壞的競爭標準。

2.1 打印時間因素對柱殼的屈曲影響

圖3 柱殼打印截面示意圖Fig. 3 3D printed section of a cylinder

一般情況下可以認為打印速度是固定的并設時間t=0時，l(0)=0，則存在關系l=l˙t。 這樣式(12)可以簡化成：

式中，t為時間。如果打印速率是變化的問題就比較復雜了， 將在今后的研究中考慮。

考慮材料在點x=0處，固化過程中彈性剛度隨時間演化可表示為：

將式(17)代入式(9)、式(10)可得歐拉坐標系下無量綱屈曲方程：

其中，無量綱參數(shù)如表1 所示。

表1 無量綱參數(shù)與其物理含義Table 1 Dimensionless parameter and physical significance

2.2 柱殼的塑性破壞

除彈性屈曲外，柱殼也會因自重因素影響，達到材料塑性屈服強度 σp而發(fā)生破壞，其中的下標“ p”代表“塑性破壞”。通常來說，塑性破壞的關鍵位置為柱殼底部，即x=0處，此處自重產(chǎn)生的軸向應力達到最大值。柱殼底部塑性破壞的屈服極限強度通?？杀硎緸椋?/p>

式中：|σp|為屈服強度絕對值；l為柱殼的長度。由于固化進程因素影響，屈服強度在圓柱殼長度增長方向是不均勻的。在x軸方向，可將 σp由 σp?代替，圓柱殼底部屈服強度隨時間的變化規(guī)律可以給定為：

式中：σp,0為打印進程中t=0時刻的屈服強度；h?(t)為屈服強度隨時間演化的固化特征函數(shù)。

由于塑性屈服極限由材料所決定，即可得出與文獻[20]完全相同的無量綱塑性坍塌長度隨固化速率變化函數(shù)：

圖4 塑性坍塌長度 lp 隨固化速率 ξσ變化曲線圖Fig. 4 Plastic collapse length lp vs. curing rate ξσ

2.3 彈性屈曲與塑性坍塌的競爭判據(jù)

當塑性坍塌長度小于彈性屈曲長度，即lplcr時，圓柱殼會發(fā)生彈性屈曲。這一用于分析可能失效機制的標準可以通過調(diào)用幾何、材料和打印參數(shù)來方便地表示：

3 柱殼失效問題數(shù)值解

3.1 特征值問題的求解過程

微分方程式(18)可利用伽遼金(Galerkin)方法進行求解：

3.2 固化進程影響下數(shù)值求解分析

采用2 次、3 次、4 次、5 次多項式基函數(shù)對圓柱體分岔屈曲數(shù)值解的精度進行了評價。從圖5可以看出，對于N=3 的多項式，臨界屈曲長度的估計明顯過高，精確性不足。對于N=4 時，精確性有所提高，但仍存在偏差。但對于N=5 和N=6，數(shù)值解在小數(shù)點后四位的結果完全相同，因此似乎已經(jīng)收斂到精確解，且滿足針對3DCP 圓柱殼屈曲長度分析需求。

圖5 無量綱屈曲長度 隨無量綱固化速率 變化曲線Fig. 5 Critical buckling length vs. curing rate E

圖6 3DCP 圓柱殼臨界屈曲長度 lcr隨固化速率 變化曲線Fig. 6 Critical buckling length vs. curing ratefor 3DCP cylinder

3.3 彈性屈曲與塑性破壞數(shù)值求解比較分析

針對柱殼彈性屈曲與塑性坍塌進行數(shù)值求解，用于分析線性固化進程中彈性屈曲與塑性坍塌的競爭關系。

圖7 3DCP 圓柱殼失效機制Fig. 7 The failure mechanism of 3DCP cylinder

3.4 實例驗證

采用WOLFS 等[17]完全相同的打印參數(shù)與柱殼半徑r=250 mm，打印層高tl=10 mm，打印噴頭移動速度vn=5000 mm/min ?；炷疗骄芏圈?2070 kg/m3。

彈性模量方程E/kPa(t/min) ：

抗壓屈服強度方程σp/kPa(t/min) ：

本文首先利用商用有限元軟件ABAQUS 并采用與WOLFS 等[17]相同的建模方法進行模擬，得出3DCP 圓柱殼失效層數(shù)為49 層，與相關研究中的數(shù)值模擬結果相差不大。在考慮初始缺陷的影響下，采用Riks 方法進行非線性屈曲分析，缺陷幅值設為厚度的1%，得出3DCP 圓柱殼失效層數(shù)為48 層，直到設置缺陷幅值為殼厚的27.5%時，得到與試驗[17]相同的失效層數(shù)，即29 層。這表明3DCP 圓柱殼存在較大的初始缺陷，但該方法很難定量分析3DCP 結構的屈曲失效行為。

在考慮到實際打印過程中層截面幾何形狀并非矩形[12]的情況下，采取WOLFS 等[17]相關工作中的圖片目視所得結果，即圓角梯形進行模擬，層截面幾何形狀如圖8 所示。

圖8 3DCP 圓柱殼圓角梯形層截面示意圖Fig. 8 Rounded trapezoid schematic diagram of 3DCP cylinder model

使用ABAQUS 中的RESTART 功能模擬3DCP圓柱殼的打印過程，其中分析步將采用具有更好收斂性的Dynamic、 Implicit 類型。RESTART 分析功能允許重新啟動模擬并對已有的載荷歷史響應進行計算，即結構的應變/位移、應力均被保留。其中每一個job 的生成代表一段物料的添加，直至模擬到結構失效。相比于目前模擬3DCP 過程的主流建模方式(MODEL CHANGE)[22]，該方法不易造成單元扭曲但計算量較大。層間接觸將被建模為基于庫侖-摩擦模型的法向和切向行為的組合，并假設摩擦系數(shù)為0.6[22]。選用三維實體八節(jié)點縮減積分單元(C3D8R)構建FE 網(wǎng)格，單層單元大約為1600 個，整個模型單元約為46 000 個。整個建模過程將利用自定義開發(fā)的參數(shù)化Python腳本進行控制。同時，利用Python 腳本對每一段打印體的材料屬性進行刷新，用以模擬打印材料的時間依賴性。

如圖9 所示，F(xiàn)EM 結果下的圓柱殼變形圖將被繪制，并與WOLFS 等[17]試驗結果進行對比，其中彩色點線圖為模擬結果，黑色實線圖為試驗結果。從結果對比可以看出，模擬與試驗的最大徑向位移十分接近，其中FEM 結果為15.44 mm，試驗結果為12.36 mm 與17.78 mm。模擬得出的破壞前結構高度為266 mm，試驗結果為270.24 mm與266.42 mm。但模擬結果得出的最大徑向位移點相對較低，且變形圖相對較陡。這些差異可能是由于目視所得層幾何形狀與實際具有一定的差異，將導致屈曲模態(tài)發(fā)生變化?？傮w來說，更改層幾何形狀能得到與試驗較為一致的模擬結果。同樣，從徑向位移云圖可以看出，圓柱殼在破壞前的變形并未呈現(xiàn)完美軸對稱情況，在結構中上部出現(xiàn)了不同程度的內(nèi)凹現(xiàn)象。

圖9 FEM 與實驗對比圖Fig. 9 Comparison of FEM and experiment

在保證層截面幾何面積一致即自重載荷不變的情況下，選取了形狀為圓角矩形的情況進行有限元模擬分析，幾何尺寸如圖10 所示。其中，圓角半徑R=12 mm ，有效殼厚為h=20.33 mm 。從結果對比可以看出，F(xiàn)EM 得出的最大徑向位移為17.68 mm ，破壞前結構高度為283.43 mm ，失效層數(shù)為25 層。模擬與試驗結果符合較好，并通過有限元對比分析可以看出，層截面的有效接觸寬度為控制3DCP 圓柱殼屈曲高度的重要因素，而非測量所得的實際殼體厚度。關于層截面幾何形狀對3DCP 圓柱殼失效形式的影響，將在后續(xù)的工作中開展。

圖10 3DCP 圓柱殼圓角矩形層截面示意圖Fig. 10 Rounded rectangle schematic diagram of 3DCP cylinder model

如圖11 所示，將FEM 所得圓角梯形截面的非線性屈曲分析結果與理論結果進行對比發(fā)現(xiàn)，理論結果相較模擬，其屈曲長度lcr較低。理論結果為25.8 層，有限元結果為29 層。二者結果之間的差異可能是因為，理論分析所得結果建立在軸對稱假設下，與實際情況是存在差距的，但模擬顯示3DCP 圓柱殼的變形形式將由軸對稱開展至非軸對稱，直至失效。但此現(xiàn)象開展非常迅速，在3 層～5 層內(nèi)便開展完畢。因此，本文認為軸對稱模型在分析3DCP 圓柱殼穩(wěn)定性時，存在一定的不足，但可作為評估屈曲破壞的下限。

圖11 最大徑向位移wmax 隨層高變化曲線Fig. 11 Max. radial deflection wmax vs. layer number

表2 相關工作失效層數(shù)對比Table 2 Comparison of failure layer number with related work

3.5 失效機制分析

圖12 3DCP 圓柱殼無量綱失效長度lˉ隨徑厚比 r/h變化圖Fig. 12 The dimensionless failure length lˉ of 3DCP cylinders vs. the radius-thickness ratio r/h

4 結論

在本文研究中，描述了3D 混凝土打印(3DCP)圓柱殼在打印進程中的力學行為。建立了兩種失效機制(彈性屈曲和塑性坍塌)的控制方程，并對這兩種失效機制間的競爭關系進行了描述。最后與Wolfs 等[17]試驗結果進行了對比驗證。本文得出的主要結論與展望如下：

(1)建立了用于描述3DCP 圓柱殼在自重作用下的結構失效行為的數(shù)值模型，該模型可準確預測打印結構的失效高度與失效形式，這是在相關工作中均未達到的。

(2) 描述了3DCP 圓柱殼彈性屈曲與塑性坍塌之間的競爭關系，其中彈、塑性固化速率之比決定了彈性屈曲與塑性坍塌這兩種失效機制的主導地位。并且在相關打印參數(shù)固定時，可以確定3DCP 結構的彈性屈曲區(qū)域與塑性坍塌區(qū)域，這為提供最優(yōu)打印參數(shù)集，針對失效機制提供特定結構增強方案，減少試錯性試驗提供了幫助。

(3) 通過試驗、模擬、理論三者結果對比分析可以看出，層截面的有效接觸寬度為控制3DCP圓柱殼屈曲高度的重要因素，而非測量所得的實際殼體厚度。但仍需要對層截面幾何形狀以及層間界面進行更廣泛的研究。