2.5D 機織淺交彎聯復合材料數字單元建模分析

2023-06-28 00:48:36王宏越王兵方國東孟松鶴

航空學報 2023年9期

王宏越，王兵，方國東，孟松鶴

哈爾濱工業大學特種環境復合材料技術國家級重點實驗室，哈爾濱 150001

陶瓷基編織復合材料由于高比模量、高比強度、耐高溫等特性，廣泛應用于航空航天領域結構部件中［1-4］。但由于制備工藝的影響，材料在微細觀結構上往往具有一定的隨機性，例如纖維單絲的排列分布［5］、孔隙和裂紋的空間位置分布［6］、纖維束方向的偏轉扭曲及截面的形態變化［7］等。這些微細觀的隨機特性導致了材料宏觀力學性能表征困難。通過建立反映真實微結構特征的單胞模型，利用仿真分析方法可以有效評估材料性能，揭示織物變形機制，進而指導編織復合材料的工藝設計，減少織造缺陷。

目前，比較常見的編織復合材料細觀模型有理想模型和精細化模型。理想模型一般包含如下假設：紗線橫截面形狀保持不變，紗線一般遵循可用基本函數表示的理想軌跡。理想模型簡潔易用，并可利用一些通用軟件實現，如Tex-Gen［8-9］、WiseTex［10］等。然而，建立的模型忽略了真實復雜的編織結構內部不規則的紗線波動、紗線擠壓、截面扭曲等特征，利用這種理想模型分析結構力學性能通常產生較大的誤差，因此精細化建模的方法備受關注?；贛icro-CT（Micro Computed Tomography）掃描切片重構的方法可捕捉編織紗線的真實變形特征［11-15］。然而，該方法需要大量微觀高精度切片，且圖像降噪、組分材料的輪廓識別與分割等處理技術極其復雜，時間和經濟成本高昂，方法的普適性較低。尤其當組分材料的密度相近時，如C/C 復合材料，圖像中的voxel 點的灰度區分并不明顯，需要在對圖像進行降噪的同時，采用特殊方法進行輪廓切割，降低了圖像處理的效率［16］。近些年，另外一種基于數值模擬動力學的虛擬建模方法可有效實現復雜織物結構高保真模型的重構，其核心難點在于纖維束力學行為的描述，比較典型的模型有珠鏈模型［17］、彈簧-質點模型［18］、彈性細桿模型［19］等。Wang 和Sun［20］將纖維束看成數字單元，模擬纖維束變形，提出了數字單元方法。由于該方法假設纖維束截面在纖維運動過程中保持不變，Zhou 等［21］進一步將該方法拓展，通過多鏈描述纖維束，以考慮紗線截面變形。Miao［22］、Huang［23］等分別改進了該方法的模擬精度和計算效率。Green 等［24-25］將該方法拓展到商業有限元軟件LS-dyna 中，拓寬了方法的適用性。El Said等［26-27］利用殼單元模擬纖維束變形，大大提高了計算效率。除此之外，Durville［28］使用不同的初始構型，對織物成型進行了模擬。Xie 等［29］在不重構模型的情況下，使用嵌套單元的方法，對成型后的復合材料進行了力學分析。Bergou 等［30］在計算機圖形領域，利用離散桿方法，通過考慮拉伸旋轉等物理行為，對彈性繩力學行為進行了模擬。實驗觀測結果表明，這種虛擬建模方法預測編織復合材料織物結構形態是高效精確的。然而，數字單元建模的精度與效率強烈依賴于單元類型的選取、單元的彈性性能的賦予、單元數量、數字單元間的摩擦系數設置。截至目前，對這些建模參數的合理選擇原則及其對模擬結果的影響尚無明確定論。此外，數字單元法分析的模型輸出結果中纖維束是離散的桿單元，不能直接應用于結構力學性能分析。因此，需要對數字單元法建模參數進行定性及定量的分析，進而建立模擬參數的選擇原則，并給出數字單元三維實體幾何重構方法，本文對研究織物變形機制及其對材料力學性能的影響具有重要意義。

本文利用數字單元法模擬2.5D 機織淺交彎聯SiO2f/SiO2復合材料的成型過程，在ABAQUS中利用通用接觸及3D 梁和桿單元建立包含隨機結構的機織體數字化離散模型。施加溫度載荷及邊界約束條件，通過分析單元類型、彈性性能、單元數量、單元間摩擦系數對模擬結果的影響，確定數字單元模擬參數設定的主要原則。提取材料變形至目標構型后數字單元的特征信息，編寫參數化建模腳本，實現3D 實體幾何模型的精確重構。

1 紗線動力學模擬

1. 1 紗線特征參數確定

研究材料為2.5D 機織淺交彎聯法向增強SiO2f/SiO2復合材料，其理想的細觀結構如圖1所示。

圖1 法向增強2.5D 機織細觀結構示意圖Fig. 1 Schematic diagram of 2.5D woven composites with normal reinforcement

為更了合理地進行動力學模擬，需要獲取較為準確的紗線特征數據。利用Diondo d2 微納米焦點Micro-CT 對SiO2f/SiO2復合材料柱體樣品進行掃描檢測，得到的切片剖面以及整體樣品如圖2 所示?？梢园l現，SiO2f/SiO2復合材料經紗較粗，截面形狀近似橢圓形，左右間距有明顯因為機織工藝引起的錯動，同時引起了上下2 束纖維的間距也發生了改變；由于經紗的錯動，使得緯紗線相互擠壓交疊，在擠壓稀疏區域近似橢圓或矩形，在擠壓密集區域近似梯形。利用Image-J 對圖片進行標注處理，獲取了不同紗線的間距、長度等特征信息，得到的數據可近似認為服從正態分布，統計信息的均值及標準差數據如表1 所示。由于緯紗線間相互擠壓，并沒有預留空間，相互間距可通過紗線截面寬度近似計算獲得。

圖2 Micro-CT 掃描2.5D SiO2f/SiO2復合材料Fig. 2 2.5D SiO2f/SiO2 woven composites with Micro-CT scanning

1. 2數字單元法建模

織物內部的紗線由數千根纖維單絲構成，將纖維抽象為離散的“數字單元”，纖維束可以簡化為多條“數字單元”聚集的結構（見圖3）。其中代表纖維束的單元通過無摩擦的點鉸接在一起。為了保留紗線的“柔性”特征，“數字單元”的長度盡可能小于單元截面直徑。分別采用了Abaqus/Explicit 中的桿單元（T3D2）和梁單元（B31）模擬“數字單元”。由于纖維的運動，相鄰纖維單鏈之間會相互接觸，因此利用通用接觸（General Contact）計算單元之間的接觸問題，取代“數字單元”中的接觸單元，并加入阻尼設置提高計算穩定性。由于單元長度小于截面直徑，在設置中應排除自相交，避免單元間的穿透距離過大。

圖3 數字單元法建模Fig. 3 Digital element modeling

1. 3 織物細觀結構隨機性

為了最終得到纖維織物的細觀結構，初始構型依據織物機織結構進行構建。相關研究表明［22］，紗線排布的初始形狀對最終計算結果并沒有影響，桿單元根數會影響最終纖維束的形態，20 根“數字單元”足以對纖維束形態進行精準描述?，F將經紗離散為61 根桿單元，緯紗離散為63 根桿單元，增強紗離散為10 根桿單元，初始位置如圖4 所示。為將不同紗線表示清楚，圖中利用不同的顏色對紗線進行區分。

圖4 2.5D 法向增強淺交彎聯織物結構離散Fig. 4 Discretization of 2.5D normal reinforced shallow cross bending fabric structure

基于Micro-CT 觀察結果，除緯紗外，其余紗線初始中心線方向均設置為直線。緯紗使用分段曲線擬合如圖5 所示，曲線組合為正弦函數與余弦函數，即

圖5 彎曲紗線離散Fig. 5 Discretization of warp yarn

式中：A、B、C、D為方程系數，可通過實際結構測量數據計算得出。

由于纖維單絲在在機織過程中僅是位置發生了變化，因此認為構成的纖維束的面積不發生改變。利用統計數據近似計算得到纖維束面積，并計算桿單元的面積。由于初始階段的離散分布特征會影響模型計算效率，故選擇一種計算量小的離散構造，如表2 所示。

表2 桿單元截面面積及其初始構型Table 2 Cross sectional area of rod element and initial configuration

為保證織物能夠打緊收縮，需設置合理的邊界條件。由于桿單元無法傳遞彎矩，故經紗端點位置在位移上采用了周期邊界條件；緯紗及增強紗端點沿紗線方向位移為0。若采用梁單元，由于其可傳遞彎矩，經紗約束條件可放松為x方向位移為0。

在計算時采用顯示的準靜態模擬分析，中心差分計算方法是條件收斂的，收斂極限與整個系統的最高頻率相關，計算公式為

式中：ωmax為系統中最高頻率；ξmax為最高頻率下模型中臨界阻尼的百分比。

收斂極限為膨脹波穿過所有單元中的最小時間，即

式中：Lmin為網格中最小單元尺寸；cd為用λ?、μ?計算得到的膨脹波速度，計算公式為

式中：ρ為材料密度；在各向同性彈性材料中，λ?、μ?計分別為拉梅第一常數和拉梅第二常數，可由楊氏模量E和泊松比ν計算得出，計算公式為

增量時間可通過每個單元的收斂時間預估，即

式中：Le為單元特征長度。

完整的機織過程的模擬分為2 個分析步：第1分析步將纖維束進行隨機移動，以達到隨機位置；第2 分析步將機織結構的整體溫度降低，邊界節點處設置位移邊界條件。結果如圖6 所示。

圖6 機織過程模擬Fig. 6 Simulation of woven process

觀察計算結果可知，對織物纖維束合理的離散，利用Abaqus/Explicit 可以完成纖維束隨機移動及織物動態收緊過程的模擬。其中，纖維束的隨機移動可以看作在編織機機織過程中，隨機擾動造成的纖維束錯位，動態收緊可以看作機織過程中的纖維束受力過程。

1. 4 數字單元法關鍵參數及其影響

1）初始截面形狀

在數字單元方法模擬中，初始數字單元離散所形成的截面構型并不影響最終織物成型后的截面輪廓［22］。雖然圓形截面具有通用性，但矩形截面與最終紗線截面形狀較為相似，因此其計算時間大大降低。由于本文研究的2.5D 編織結構較為復雜，纖維間距較為緊湊，因此初始離散構型不能過于分散，否則無法預估材料屬性及纖維束移動的距離。當初始離散構型較為緊湊時，為了保證纖維束截面面積不發生變化，可采用與真實材料纖維束截面形狀相似的形狀。由于考慮了織物中纖維束位置的離散性，纖維束在模擬過程中將移動到隨機指定的位置，運動距離將受初始纖維束離散截面影響。對于特定的織物結構，初始纖維束離散截面形狀可與最終形狀相似，不僅可以較好地保留纖維束彎曲及表面形貌的真實特征，同時也可以節約計算時間并降低計算復雜度。

2）材料屬性

利用桿單元進行模擬，通過施加合理的邊界條件，可以發現模量的大小并不會影響最終的結果，如圖7 所示。這是由于桿單元無法傳遞彎矩，模量僅影響桿單元的收縮力與接觸力，模量在較大范圍內變化，降低的溫度足夠使各纖維束在最終位置處于收緊狀態。但在模量較小的時候，纖維束較其他情況更松散，這是由于接觸力較小使得計算出的摩擦力較小，從而導致纖維束之間的約束力降低。

圖7 不同模量桿單元模擬結果對比Fig. 7 Comparison of simulation results obtained by different modulus of rod element

利用梁單元進行模擬時，為保證單元能夠模擬纖維柔軟的性質，需要額外定義剪切剛度。若材料模量較高，則導致單元變形困難，模擬得到的纖維束彎曲程度更接近初始形狀，直纖維不易彎折，彎曲纖維傾向于保持彎曲，模量過大會導致模擬過程中織物無法變形。模量越高，在相同時刻，整個織物較離散前的形態變化就越小。不同梁單元特征尺寸對模擬結果的影響如圖8 所示，由式（7）可得，使用梁單元進行模擬時，其模量范圍與其單元尺寸具有明顯相關性，合理的模量不僅可以更好地模擬真實纖維柔軟的特性，更能增加計算收斂速度，模擬出的織物形態更加真實?？筛鶕剑?）對劃分好的網格進行計算，得到材料的模量與密度。本文所采用的材料密度為2 200 kg/m3，經紗模量為1.65 MPa，緯紗模量為3 MPa。

圖8 不同梁單元特征尺寸對模型的影響Fig. 8 Effect of different beam element characteristic sizes on the model

通過圖7、圖8 的結果對比可知，由于梁單元具有彎曲剛度，采用梁單元較易模擬結構松散的織物；梁單元變形能力極易受模量影響，若織物紗線成型后位置距離散前位置移動過大，采用梁單元進行計算將降低收斂性。故梁單元適用于模擬紗線體積分數小，位置變化小的織物，桿單元適用于模擬紗線體積分數大、擠壓變形大、位置變化大的織物。

3）虛擬纖維束中數字單元數

若將纖維束中每根纖維單絲都建模為梁單元或桿單元顯然是不合適的。計算成本極大限制了建模時所使用代表纖維束的單元數量。文獻［22］提出對于模擬2D 平紋機織織物，為準確描述每根纖維束截面需要使用52 條數字單元，3D 編織織物需要19 條。對于本文所構建的2.5D 機織復合材料，其緯紗與經紗具有完全不同的特征，緯紗較為緊湊，截面大多呈梯形；經紗橫縱比較大，呈橢圓形?，F采用不同數量的單元對機織結構進行模擬，如圖9 所示。

圖9 不同單元數量模擬結果對比Fig. 9 Comparison of simulation results obtained by different element numbers

通過觀察可以發現，盡管單元數量不同，最終模擬所得纖維輪廓較為相似。單元數量越多，模擬所得到的纖維束輪廓越清晰，纖維束的形態越自然。對緯紗而言，當單元數量減少，即單元半徑增加時，由于接觸判定范圍增加，導致單元彎曲形態不平滑。本文認為，纖維束中的單元數量初始可試取20～50，后可根據計算結果與真實形貌進行調整。由于單元數量會影響接觸判定范圍及接觸搜索，在考慮計算成本的情況下，不建議選擇使用過多數量。

4）摩擦系數

摩擦系數是織物成型過程起重要作用的參數。纖維束與纖維束之間，纖維束內的纖維單絲之間均存在摩擦。在模擬中，分別考慮摩擦系數為0.1、0.3、0.6 的情況，若在模型中不考慮隨機性，由于纖維束內部存在較大摩擦，纖維束橫截面更加緊湊，即摩擦力起到了阻止纖維束截面變形的目的；若考慮隨機性，當纖維束移動到位置時，低摩擦系數將會導致纖維束位置錯動，當摩擦系數高于0.3 時，纖維束最終位置與第1 分析步最終位置較為一致，不會隨著緯紗的拉緊產生位置錯動，如圖10 所示。

圖10 不同摩擦系數織物形態Fig. 10 Fabric morphology with different friction coefficients

觀察織物最終形態可以發現，增加摩擦系數可有效增加織物成型厚度。由于摩擦系數增加，導致單元間摩擦力增加，單元更容易攢聚在一起，形成纖維束截面長軸變短、短軸變長，使得整體織物高度增加。通過調整摩擦系數，將有效調整纖維束截面形態。但摩擦系數不易過高，隨著摩擦系數的增加，計算時間將成倍增加，計算成本大大增加。設定在溫度320 ℃時摩擦系數為0.1，隨著溫度降低至300 ℃，摩擦系數線性增加為0.3。

利用建立的基于數字單元法的動力學建模方法模擬得到的2.5D 機織淺交彎聯SiO2f/SiO2復合材料的織物形態如圖10 所示，采用的計算參數如表3 所示。

表3 關鍵計算參數Table 3 Simulation parameters

通過圖11 可以發現，模擬結果與Micro-CT觀測結果吻合一致，驗證了建模方法的正確性。

圖11 模擬結果與Micro-CT 觀測結果對比Fig. 11 Comparison of simulation and Micro-CT observation

2 機織復合材料高保真模型重構

通過數字單元方法模擬得到了機織復合材料細觀纖維束結構的形態，還需要進一步提取每根纖維束的信息，才能實現對機織復合材料細觀結構的高保真重構。

2. 1 纖維束三維空間描述

首先對纖維束的三維空間幾何進行嚴格的數學描述，表示為P，利用纖維束的中心路徑S和截面C定義纖維束，公式為

式中：X′和Y′為纖維束路徑的局部坐標。

式中：U為纖維束“上”向量，如圖12 所示，該向量方向垂直纖維束中心線。假設織物在XY平面上，則默認其為（0，0，1）。

圖12 纖維束路徑坐標計算方法Fig. 12 Calculation method for yarn path

2. 2 纖維束中心路徑插值

纖維束的中心線可利用線上離散點插值來描述。選取適當的插值函數可利用少量插值點足夠準確地描述纖維束形態。為保證纖維束在切線方向上光滑并且纖維束與纖維束之間沒有空隙，插值函數至少應該是C1連續的，最常用的方法是樣條插值。由于高階樣條插值會出現Runge 現象，低階樣條插值在實際應用中更有效。

式中：ti（i=0，1，…，k-1）為節點。常見的插值類型為三次Bezier 樣條插值、三次自然樣條插值、三次周期樣條插值，這些插值方法都是三階的，后兩者為C2連續的。本文采用三次Bezier 樣條插值。

一般可利用4 個節點P1、P2、P3、P4構造Bezier樣條插值曲線，如圖13 所示。曲線從P1開始朝向P2方向，最終經由P3方向到達P4。在P1點處，切線方向從P1指向P2；在P4點的切向方向由P3指向P4。三次Bezier樣條插值的參數化方程為

圖13 三次Bezier 樣條曲線Fig. 13 Cubic Bezier spline curve

為保證C0連續，點P1i、P4i可用分段的樣條曲線Si定義

為保證C1連續，點P2i、P3i的選取需滿足

在t=a和t=b處的斜率并沒有被定義。在大部分情況下，可以認為纖維是周期性的，因此，可以假設

如果纖維是非周期的，可假設為

這些假設可以保證生成光滑的曲線。若在特殊情況下，已有參數不足以描述纖維束，可以在特定節點處調整切線方向向量以使曲線光滑。

2. 3 纖維束截面的插值與離散

纖維束的截面是由垂直纖維束路徑的平面與纖維束相交切出的2D 圖形。一般情況下，認為纖維束是實體，纖維束截面被認為是所有纖維壓縮在一起形成的最小區域。截面的邊緣可以通過二維參數方程進行描述。假設纖維束的2 個截面A(t)、B(t)，兩者之間的截面C(t)可以定義為

在截面A(t)、B(t)之間線性插值是保證纖維束C0連續最簡單的方式。此外，截面A(t)和截面B(t)的定義方式必須是相似的，截面輪廓的點需要在相似的位置。

三維實體表面最簡單的顯示方式是將其表面離散為插值曲線并利用插值曲線表示，可以采用OpenGL 或者Direct3D 對圖形進行渲染。獲取表面插值曲線的方法就是表面離散，如圖14 所示。除了顯示纖維束，表面離散同時可以用于計算體積或者邊界單元分析。將式（8）離散，假設a、b分別為0、1，將其用k和n個點平均分，則每個線段長度可表示為1/(k-1)和1/n，即

圖14 參數化表面的離散Fig. 14 Discretization of parameterized surfaces

因為在v=1 時，和v=0 時是同一個點，所以式（24）的分母為n而不是n-1。

離散點的密度由k、n控制，一般需要讓所有的分段長度近似。邊長Lˉ的計算公式為

式中：[ ·]為取整計算。

2. 4 幾何模型重構及分析

實體幾何模型重構最主要的過程就是纖維束輪廓的提取與纖維束中心線方向的計算。兩者是相互迭代的過程：提取纖維束輪廓，必須首先知道纖維束中心線方向；計算纖維束中心線方向，必須通過纖維束輪廓計算得到線上的點（該點為纖維束輪廓的中心）。本文首先假設已知纖維束中心線方向，即從一束桿單元中，選擇靠近中心位置的桿的方向作為初始纖維束中心方向，從而計算得到纖維束輪廓，利用獲取的纖維束輪廓，計算得到中心點，再次修正纖維束中心線方向，重復該過程直到前后兩次的纖維束中心線方向誤差不超過1%，一般三次迭代可以達到最終目標，如圖15 所示［14］，其中坐標軸為纖維束空間位置坐標。

圖15 纖維束輪廓及中心路徑提取流程［14］Fig. 15 Yarn outline and center path extraction process［14］

利用計算幾何中尋找凸包的方法提取纖維束輪廓，再利用插值節點和該節點處的方向向量，構建虛擬平面與其相交，通過求解方程得到截面處的空間位置坐標。然而在三維空間中，直接尋找凸包是十分困難的。若利用虛擬平面求得其與纖維束相交的交點，可在虛擬平面內尋找到這些點的凸包，利用坐標變換將二維平面上的點重新轉換為三維空間中的點，這樣就得到了纖維束截面的輪廓。這樣的處理將三維問題轉化為了二維問題，降低了求解難度。

本文采用的尋找凸包方法為Scipy 模塊中的ConvexHull 函數。利用該函數可近似得到點集的凸包，并返回控制點坐標。對這些不均勻的控制點進行重新計算，得到截面坐標均勻分布滿足截面穩定插值。通過觀察未發現纖維束加捻情況，對纖維束作平滑假設，假設其截面變化均由擠壓變形導致，自身并不發生旋轉，因此在計算中并沒有考慮旋轉角。以緯紗為例，共構造63 根單元，通過觀察，在計算過程中，單元之間并無明顯位置錯動，為加快迭代周期，選取第30 根單元的方向向量作為初始假設的纖維束方向向量。將方向向量按節點離散，在節點處取垂直于方向向量的平面，將63 根單元節點向平面投影。利用Scipy 模塊中的ConvexHull 函數獲取集合點的外輪廓點。由于所獲取的外輪廓點不均勻，現利用Geomdl 模塊中的Bezier 樣條插值命令，將輪廓點插值為閉合曲線，并返回所得到1 000 個非均勻的點，最終將整個閉合曲線平分，返回得到25 個近似均勻分布的點，這些點所圍成的輪廓即為纖維束截面在二維平面上的投影。利用平均角度的方法，將閉合的輪廓分為24 個均勻的角度覆蓋區域，取角度平分線與輪廓的交點為輪廓特征點。由圖16 所示，輪廓線整體上貼合單元外側。雖然較少的單元數即可表示較為真實的纖維束輪廓，但在輪廓提取時，較少數量的單元會使得個單元截面積增加，在單元發生較大錯動時產生較大空隙，這些空隙在一定程度上會對纖維束輪廓的提取范圍造成影響。

圖16 纖維束截面輪廓提取Fig. 16 Fiber bundle profile extraction

此投影并非最終真實纖維束空間截面輪廓，需將二維平面上的點轉換為空間坐標下的點，得到最終真實纖維素空間截面輪廓和中心線如圖17 所示。

圖17 纖維束截面提取與坐標變換Fig. 17 Fiber bundle cross section extraction and coordinate transformation

重構后的織物如圖18 所示，可以發現，本文所采用的重構方法能夠較完整地完成模型重構。在整體織物規模較大的情況下，細觀結構并非簡單的重復，模擬的數值波動會帶來細微的結構隨機性，織物越大，變化程度越復雜。同時存在紗線間相互交疊難以提取特征的問題。利用TexGen進行了穿透檢驗后發現，重構后的模型整體并無明顯穿透跡象，無論是單元數較少的增強紗，或是纖維束方向變化劇烈的緯紗，或是較為簡單的類直線經紗，本文所采用的重構方法都能夠較真實地保存其特征，即便在纖維交疊處，在紗線表面也并未觀察到嚴重的扭曲，因此該方法具有一定的普適性。

圖18 織物重構Fig. 18 Fabric reconstruction

提取模型重構后的經紗和緯紗如圖19 所示。在紗線相互搭接處，紗線由于接觸力會發生變形。由于受到壓縮，紗線截面趨向橢圓。在交疊位置之間的過渡區域，紗線由于失去約束，紗線截面重新擴大，并形成紗線褶皺。因此，本文所采用的重構方法能夠較真實的捕捉纖維束細節，能夠將真實的纖維束形貌引入后續力學性能分析中，后續研究中將繼續對紗線細觀形貌及其引起的性能偏差進行分析。