基于數字單元法的三維正交織物微觀幾何結構建模

馬 瑩，何田田，陳 翔,3，祿 盛,3，王友棋

(1. 重慶郵電大學 先進制造工程學院，重慶 400065； 2. 堪薩斯州立大學 復合材料實驗室，堪薩斯 66506； 3. 西安交通大學 機械結構強度與振動國家重點實驗室，陜西 西安 710049)

高性能纖維增強復合材料具有高比剛度、高比強度、抗沖擊和耐疲勞等多種優異的力學性能，在航空航天、軍事、醫療等諸多領域[1]得到廣泛的應用。相比其他復合材料增強體而言，三維機織物具有多取向性、可設計性以及結構效能優異等特點，尤其在沖擊吸能阻尼結構方面優勢巨大，被廣泛用于人體防彈衣等柔性防護裝備。

目前已有眾多學者應用多種建模軟件對不同三維織物進行了細觀結構建模，并從理論上分析其細觀結構。王旭[2]應用動畫渲染和制作軟件3ds MAX(3D Studio Max)基于曲線控制點的紗線軸線生成方法，結合截面曲線放樣技術建立了三維織物的細觀模型；陳振等[3]將紗線假設為連續實體，討論了織物幾何結構建模軟件TexGen的仿真建模方法及其結果的準確性。雖然該軟件所建的特定幾何模型與真實織物圖像對比有一定的擬合度，但并不能很好反映所有類型織物的紗線形態。余育苗等[4]以三維正交機織復合材料為研究對象，假設紗線橫截面為矩形截面，結合三維機織復合材料的結構特點建立了單胞模型。

上述建模方法較為真實地反映了三維織物細觀幾何結構，然而所建數值模型大都以紗線為最小單位，將紗線橫截面理想化地假設為橢圓、跑道、凸透鏡等形狀；但實際織造過程中織物橫截面是動態變化的，與假設的截面形狀存在較大差異，因此，為了建立更加接近三維織物真實形態的數值模型，有學者在原有理想模型基礎上提出了多種改進方法，Green等[5]建立了一種用于預測三維織物在織造和壓實過程中變形情況的精確模型；Fredrik等[6]將紗線模擬為鏈式結構，在細觀尺度上提出了一種管狀編織物幾何建模方法，該方法展現了光滑的紗線路徑和截面形狀變化；Isart等[7]分別通過將紗線截面形狀理想化的建模方法、數字單元法和理論分析方法建立了織物幾何模型，并分析評價了3種建模方法的優缺點。

三維機織物的幾何結構與力學性能高度關聯，研究其表征技術，對該類復合材料的有效利用和深層次開發具有重要工程實際意義。目前已有研究大都是參考織物截面顯微圖像假定其理想幾何形狀，然后在建模軟件中進行參數設置并一步成型，不能有效反應織物內部復雜微觀幾何結構。因此，針對現有研究中假設紗線截面為理想形狀建模的問題，本文采用由堪薩斯州立大學復合材料團隊研發的紡織建模軟件數字織物力學性能分析器(Digital Fabric Mechanics Analyzer，DFMATM)。該軟件以數字單元法為理論基礎，可實現三維機織物織造過程動態仿真及其微觀幾何結構數值模擬[8]，現已被業界廣泛用于構建二維平紋、三維正交[9-10]、角聯鎖等復雜織物組織。與其他紡織建模軟件相比，數字單元法在近似于纖維尺度建立了織物幾何結構數值模型，規避了對紗線截面形狀和材料彈性常數等的簡化及假設，在時域中模擬了織造過程，能夠真實反映紗線截面的動態變化。

本文以數字單元法為理論基礎，通過分析該方法在織物織造行為模擬中的作用機制，提出了一種計算纖維間摩擦力的方法。進而研究了三維正交織物組織結構關鍵點位置，建立其拓撲結構，并在時域中對該織物的織造行為進行模擬，獲得了5個精度遞進的單胞數值模型。通過實驗對比，分析了纖維間摩擦力對織物內部節點力和勢能的影響，揭示了紗線纖維化離散程度對仿真時間、織物厚度、纖維體積分數和紗線空間構型的影響規律。

1 數字單元法建模原理

1.1 基本要素

數字單元法主要包括3個基本要素：數字纖維、數字紗線和接觸單元，如圖1所示。可以看出，數字纖維由節點和桿單元鏈接而成，節點與桿單元間無摩擦作用。當桿單元長度趨近于零時，數字纖維可靈活彎曲，可模擬出實際織物中纖維單絲的真實形態，并賦予其真實材料屬性。數字紗線由多根數字纖維組成，數字纖維的數量由真實紗線所包含的纖維數量和組成符合真實紗線截面形狀所需數字纖維數量共同決定。數字紗線通常包含10～100根數字纖維，數字纖維的排列方式決定了紗線的截面形狀。在仿真過程中，當相鄰纖維之間的間距小于數字纖維直徑時，建立接觸單元。

圖1 數字單元法三要素Fig.1 Three elements of digital element approach. (a) Digital fiber; (b) Digital yarn; (c) Contact element

數字單元法仿真結果的準確性主要由網格劃分精細程度決定，使用該方法進行織造過程的仿真，通過網格劃分的紗線纖維化離散和纖維離散實現。在紗線纖維化離散中，1根紗線被離散為與原有紗線路徑平行的多根數字纖維，且離散前后紗線截面積和所有數字纖維的截面積總和保持不變。在纖維離散中，組成數字纖維的桿單元長度減小，節點數量增加，且桿單元長度與數字纖維直徑呈定比。在張力的作用下，因紗線纖維化離散產生的數字纖維間發生相對運動，紗線截面形狀沿軸線方向發生變化。纖維離散保證了計算纖維間作用力的準確性，2種離散機制的共同作用實現了三維機織物微觀尺度幾何結構預測，因此，紗線中的數字纖維數量和桿單元長度是影響仿真結果準確性的重要因素。

1.2 纖維間的相互作用力

為了真實反映織物織造過程中紗線的受力情況和成型過程，數字單元法不斷優化改進，在力學模型中先后加入了紗線張力、接觸力[11]、阻尼力[12]與摩擦力。

纖維間的摩擦存在于整個織造過程中，摩擦力的計算方法影響了紗線的幾何形狀和織物勢能，因此，依據纖維間的相互作用機制可建立其摩擦力計算模型。纖維間的摩擦力由節點的運動狀態和接觸單元間的相互作用力大小共同決定。摩擦力計算方法示意圖如圖2所示。

圖2 摩擦力計算方法示意圖Fig.2 Illustration of friction force calculation. (a) Before translation; (b) After translation

節點i和節點j組成1對接觸單元，分別屬于2根不同的纖維。平移前該接觸單元所處的相對位置如圖2(a)所示，從t時刻到t+1時刻，節點i和j分別移動到i′和j′位置。令ij為向量v0，i′j′為對應向量v′。令節點i和j在t時刻首次發生接觸，并在t+1時刻仍然保持接觸，則摩擦從t+1時刻開始計算。將向量ij平移到i′j′處，使點i與點i′重合，平移后的接觸單元如圖2(b)所示。令節點j相對節點i的相位移矢量為u；相對切向位移矢量為us；相對法向位移矢量為un。摩擦力可由式(1)、(2)計算得出。

Fs=ksus

(1)

ks=μkn

(2)

式中：Fs為纖維間摩擦力，N；ks為接觸單元摩擦剛度，N/m；μ為摩擦因數；kn為接觸單元法向摩擦剛度，N/m。

當|Fs|≤μFn，節點i，j無相對運動，Fs=μknus；當|Fs|>μFn，節點i、j發生相對運動，摩擦力大小為Fs=μFnus。其中，Fn為接觸單元法向相互作用力。

如圖2(b)所示，從t時刻到t+1時刻，節點i、j的切向位移矢量us可由式(3)、(4)、(5)計算得出。

u=v′-v0

(3)

un=uv′2/|v′|2

(4)

us=u-un

(5)

2 建模方法

2.1 織物拓撲結構的建立

本文建立的三維正交織物模型由低結晶碳化硅纖維組成，紗線截面積為8.66×10-8m2，纖維軸向彈性模量為190 GPa，纖維橫向彈性模量約為軸向的1/10，纖維密度為2 500 kg/m3。三維正交織物組織結構如圖3所示，可以看出虛線框內為1個代表性體積單元。該單元內2根結構對稱的接結經紗捆綁1列緯紗，浮長為1個組織點。由真實織物樣本測得該結構代表性體積單元長度、寬度及厚度，分別為0.002 903、0.001 587 5和0.003 82 m。織物拓撲結構建立步驟：1)將經紗、緯紗和接結經紗分別定義為3種不同的紗線類型。緯向包含緯紗，經向包含經紗、接結經紗。2)依次將緯紗、經紗和接結紗結構轉換為矩陣表征。

圖3 三維正交織物組織結構Fig.3 Three-dimensional orthogonal fabric structure

由圖3可以看出：經紗共有9列，由上至下分別對應編號1#～9#；緯紗共有10行，由上至下依次對應編號11#～20#；經緯紗線交替排列，經紗列數比緯紗行數少1。經紗的位置矩陣由與緯紗的相對位置決定。令位于19#和20#緯紗之間的9#經紗位置矩陣為(1,1)，位于18#和19#緯紗之間的8#經紗位置矩陣為(2,2)，以此類推。

接結經紗的位置矩陣定義方式與經紗類似。0#和10#接結經紗的捆綁方向(連接緯紗方向)呈對稱分布，其位置矩陣由與緯紗的相對位置決定。例如，0#接結經紗位于第1列緯紗之下和第2列緯紗之上，其位置矩陣為(10,0)。同理，10#接結經紗位于第1列緯紗之上和第2列緯紗之下，其位置矩陣為(0,10)。受周期性邊界條件的約束，接結經紗的首尾高度應保持一致，由內置算法自動生成。

通過以上方法得到三維正交織物拓撲結構，如圖4所示。

圖4 三維正交織物拓撲結構Fig.4 Three-dimensional orthogonal woven topology

2.2 織物建模過程

數字單元法建模流程如圖5所示。數字單元法建模過程對應以下5個步驟：1)輸入材料參數。2)將織物組織結構轉換為矩陣表征，在系統中建立織物拓撲結構。3)設置迭代參數，在時域中利用多次迭代分析模擬整個織造過程。4)虛線框內為整個迭代分析過程。該過程以數字桿單元和節點為基數循環。首先判斷和創建纖維間的接觸單元對，并計算其相互作用力大小；其次通過紗線纖維化離散和纖維離散，提高模型精度和獲得紗線真實截面形狀。5)最后判斷模型是否穩定，即達到最小勢能，如果未穩定則返回步驟4)繼續分析，穩定則輸出最終模型。

圖5 數字單元法建模流程Fig.5 Flow chart of digital element approach modeling

3 織造過程模擬與結果分析

根據文獻[13-14]，選取本次仿真所用緯紗、經紗和接結經紗張力值分別為0.2、0.2、0.01 N，分別建立了精度遞進的5種單胞模型。不同仿真時間的三維正交織物織造過程模擬如圖6所示。每種模型的紗線分別由4、7、12、19、37根數字纖維組成，其織造過程大致分為3個階段。模型1如圖6(a)所示，由2次離散得到4根數字纖維模型：第1次離散是在初次運算后將紗線纖維化離散為2根數字纖維；第2次離散是在將每根數字纖維再離散為2根數字纖維。模型2如圖6(b)所示，由1次離散得到7根數字纖維模型。模型3如圖6(c)所示，由2次離散得到12根數字纖維模型，第1次將紗線分成4根數字纖維，第2次將每根數字纖維離散為3根數字纖維。模型4、5分別如圖6(d)、(e)所示，均由1次離散得到19、37根數字纖維模型。其中，圖6(b)、(d)和(e)所示的織物單胞模型，經紗線纖維化離散后其紗線截面形狀從初始的圓形變為不規則的近似真實的形狀。

圖6 不同仿真時間的三維正交織物織造過程模擬Fig.6 Simulation of three-dimensional orthogonal woven fabric weaving process with different simulation time. (a) Model 1; (b) Model 2; (c) Model 3; (d) Model 4; (e) Model 5

3.1 纖維間摩擦對節點力和勢能影響

織物單胞勢能U由纖維間的張力勢能Ut和接觸勢能Uc組成，其中接觸勢能由2個接觸單元內的節點作用力和摩擦力構成，可由式(6)計算得出。

(6)

式中：Uci為接觸勢能，J;i為接觸單元對的編號；2個相互接觸的節點編號分別為i1和i2，其接觸勢能Uci1和Uci2可由式(7)、(8)計算得出。

(7)

(8)

式中，rip為Uci1和Uci2之比，可由式(9)計算得出。

(9)

式中：Ei1為節點i1所在纖維的軸向彈性模量，Pa；Ei2為節點i1和i2所在纖維的軸向彈性模量，Pa；Li1為節點i1所在纖維的桿單元長度，m；Li2為節點i2所在纖維的桿單元長度，m。

接觸勢能Uci可由式(10)計算得出。

Uci=0.5Fi(Ri-lci)

(10)

式中：Ri為節點i1所在的纖維半徑和節點i2所在的纖維半徑之和,m；lci為節點間距，m；Fi為接觸力,N。Fi可由式(11)計算得出。

(11)

節點平均作用力和時間關系見圖7。可以看出：37根和7根數字纖維的節點平均作用力在計算初期，出現較大波動；20 ms以后，只有12根數字纖維模型在60 ms附近出現小幅波動，其余4個模型的節點平均作用力均逐漸平衡，并在仿真結束時趨近于0。

圖7 節點平均作用力和時間關系Fig.7 Relationship between average nodal force and time

織物單胞勢能和時間關系見圖8。可以看出，5個模型的單胞勢能在仿真過程中均呈現下降趨勢。當節點平均作用力降低，勢能也隨之降低；當勢能趨于平衡時，節點平均作用力趨近于0，此時單胞中僅剩張力勢能。因此，紗線纖維化離散程度越高，纖維間的張力勢能越小。

圖8 織物單胞勢能和時間關系Fig.8 Relationship between unit-cell potential energy and time

圖9 織物單胞模型與切片顯微圖像的正視圖對比結果Fig.9 Front view comparison results between fabric unit-cell model and microscope picture. (a) Model 1; (b) Model 2; (c) Model 3; (d) Model 4; (e) Model 5

圖10 織物單胞模型與切片顯微圖像的側視圖對比結果Fig.10 Side view comparison results between fabric unit-cell model and microscope picture. (a) Model 1; (b) Model 2; (c) Model 3; (d) Model 4; (e) Model 5

3.2 紗線纖維化離散程度對幾何結構影響

三維正交織物織造過程數據如表1所示。每個織物模型的織造過程分為3個階段，每個階段用迭代次數劃分。每次迭代后獲得的仿真結果與圖6三維正交織物織造過程模擬中展示的單胞微觀幾何結構一一對應。每個階段的紗線纖維化離散程度、纖維體積分數以及累計耗時可由表1示出。

表1 三維正交織物織造過程數據Tab.1 Simulation parameters of three-dimensional orthogonal woven fabric weaving process

織物單胞模型與切片顯微圖像的正視圖對比結果見圖9，織物單胞模型與切片顯微圖像的側視圖對比結果見圖10。

每根紗線由19根數字纖維組成的織物模型與顯微圖片擬合度最高，且隨著紗線纖維化離散程度提高，模型幾何結構更均勻，紗線空間構型更精確。5個模型厚度略低于真實織物厚度，這是因為在數值模擬中，模型的勢能達到最小時方可結束計算，而三維機織過程是個動態的過程，織物的微觀幾何形狀受織造速度和紗線張力等因素的影響，較難達到最小勢能狀態。

織物模型厚度隨時間的變化如圖11所示。可以看出，紗線所含數字纖維數量越多，結構穩定時織物的厚度越小，織物厚度下降速度越慢，達到穩定所需時間越長。其中37根數字纖維模型計算耗時約為其余模型的2～5倍。

圖11 織物模型厚度隨時間的變化Fig.11 Change of fabric thickness with time

通過對比織物模型與實驗所得織物樣本的6個因素(微觀幾何結構、節點平均作用力、勢能、仿真耗時、織物厚度、纖維體積分數)可以看出：37根、12根和7根數字纖維模型的節點平均作用力在仿真的不同時期出現波動；紗線纖維化離散程度越高，織物結構穩定時的勢能越小、厚度越小、纖維體積分數越大，紗線空間構型越接近真實紗線形態。綜合所建模型穩定性、仿真耗時和與切片顯微圖像的重合度得出，對于三維正交織物而言，將每根紗線離散為19根數字纖維的建模方法為較優。

4 結 論

為提高三維正交織物幾何結構建模精度，并為之提供一種科學有效的仿真方法，運用數字單元法理論，建立和探討了5種精度遞進的微觀幾何結構數值模型，得到如下結論。

1)紗線纖維化離散程度的提高增加了仿真時長，減緩了織物勢能、節點平均作用力和厚度的下降趨勢；當織物結構穩定時，紗線纖維化離散程度越高的模型其厚度和勢能越小。

2)綜合仿真時長、織物節點平均作用力和勢能變化情況與織物樣本厚度和顯微圖的擬合程度，將每根紗線一次性分為19根數字纖維的仿真方法較為合理。

3)采用數字單元法仿真所得的三維正交織物微觀幾何結構，較為準確地反映了真實織物樣本的紗線空間構型，為后續力學性能的研究奠定了良好的理論基礎。