危銀濤 羅亦文 繆一鳴 柴德龍 馮希金

1.清華大學(xué)汽車安全與節(jié)能國家重點實驗室，北京，100084

2.貝卡爾特鋼絲簾線亞洲技術(shù)中心，江陰，214434

3.杭州朝陽橡膠有限公司，杭州，310008

0 引言

有限元方法在輪胎研究中的應(yīng)用始于20世紀(jì)70年代［1－7］。目前已發(fā)展到用有限元方法分析輪胎 滾 動 損 失［8－10］、斷 裂［11－15］和 輪 胎 動 力 學(xué)［16－23］等方面。為了研究子午線輪胎鋼絲簾線的受力與變形規(guī)律，在建立輪胎整體模型的過程中，采用rebar（加強筋）單元［24］模擬鋼簾線。rebar單元本質(zhì)上是一維桿單元，只能承受軸向拉壓變形，而實際鋼簾線是由若干單絲捻合成股后，再由股捻合而成的三維空間結(jié)構(gòu)，在輪胎結(jié)構(gòu)中除受拉外，還有彎曲變形，而現(xiàn)有的整體模型不足以全面細(xì)致地揭示簾線的力學(xué)性能。

筆者基于已有的輪胎整體建模理論，提出一種研究簾線內(nèi)部應(yīng)力和應(yīng)變的局部分層建模新方法。在局部尺度建模過程中，將簾線的設(shè)計變量（簾線捻角、捻距、簾線中心半徑）作為建模輸入?yún)?shù)建立含單絲信息的有限元模型，將從整體模型中獲得的鋼絲簾線典型變形和應(yīng)力作為載荷邊界條件。

1 全鋼載重子午線輪胎的整體建模

基于rebar技術(shù)［24］，在充氣壓力為840kPa、負(fù)荷為3675kg的條件下，對某全鋼載重子午線輪胎進(jìn)行整體建模。圖1a顯示的是輪胎整體模型中的材料分布狀況；圖1b是鋼絲簾線的分布示意圖。

1.1 簾線張力

對于載重子午線輪胎，鋼絲簾線是承載的主要部件。鋼絲的分布不僅影響自身受力，還影響膠料的剪應(yīng)力和剪應(yīng)變，最終影響到輪胎的耐磨性、剛性、生熱性和滾動阻力等。充氣壓力條件下，第二帶束層簾線的張力大小沿簾線呈拋物線形狀，簾線中心位置（冠部）為拋物線頂點，此處的張力最大，為372N；簾線兩端張力最小，為4N。對于胎體簾線，中間（冠部）區(qū)域的受力要小于兩胎側(cè)區(qū)域的受力，從上胎側(cè)到下胎側(cè)，簾線受力幾乎不變，胎圈部位簾線張力有所增加，但到達(dá)峰值之后簾線張力開始下降，到簾線端點張力下降到零。

圖2 第二帶束層簾線受力分布

圖3 胎體簾線受力分布

在充氣載荷和垂直載荷的共同作用下，輪胎簾線的受力狀況不同于純充氣狀態(tài)的受力狀況，主要表現(xiàn)為對接地區(qū)域的影響。圖2、圖3分別為第二帶束層簾線與胎體簾線的受力分布三維圖，底面上兩個互相垂直的坐標(biāo)軸分別表示沿子午方向的寬度（沿胎體簾線長度）和輪胎的周向角度（輪胎最高點為0°，最低點為180°），縱軸表示輪胎簾線張力。從圖2可知，第二層帶束層簾線為主要承載層，承受大部分的充氣壓力和外部負(fù)荷，帶束簾線受力具有3個特點：①簾線的受力與圓周的胎面中心線呈反對稱；②在接觸區(qū)域，簾線受力在冠部區(qū)域減小而在帶束層的末端增大；③簾線受力的最大值轉(zhuǎn)移到胎肩區(qū)域的帶束末端，最小值轉(zhuǎn)移到冠部中心區(qū)域，即接觸使帶束簾線受力在冠部區(qū)域松弛，在胎肩區(qū)域增加。從圖3可知，胎體簾線應(yīng)力分布的一個重要特點是：在輪胎與地面的接觸區(qū)域，簾線應(yīng)力有所減小；在接觸區(qū)域以外，胎體簾線應(yīng)力基本不變。

1.2 胎體彎曲應(yīng)力

如果胎體曲率已知，那么簾線的應(yīng)力［25－27］：

式中，r為單絲半徑，r＝0.11mm；E為簾線材料的彈性模量，E≈200MPa；Δκ為胎體曲率。

式（1）表明簾線中的單絲均獨立彎曲，且曲率相同。因此，一根單絲的曲率可以從的胎體整體變形中獲得

式中，κu為胎體變形后的曲率；κ0為原始曲率；uy為y方向的位移，由整體分析得到。

變形前后y的二階偏導(dǎo)數(shù)可用差分求得

由式（2）、式（3）分別可以得到彎曲應(yīng)力和彎曲曲率，為進(jìn)一步分析簾線的局部變形提供了載荷邊界條件。簾線的局部變形包括絲與絲之間接觸應(yīng)力需要更精細(xì)的局部模型。

2 復(fù)雜多股鋼絲簾線的有限元建模

為了全面細(xì)致地揭示簾線的變形機(jī)制，創(chuàng)建簾線局部分層有限元模型。在局部建模中采用的胎體彎曲曲率為25m－1（或曲率半徑為40mm），同時，對簾線施加從50～300N的張力。

2.1 復(fù)雜鋼簾線雙螺旋結(jié)構(gòu)的幾何描述

為了給復(fù)雜的多股鋼絲簾線劃分網(wǎng)格，首先要精確地描述其幾何變形特征。圖4所示是一種3＋9＋15×0.22＋0.15的三層鋼絲簾線結(jié)構(gòu)，其內(nèi)層為3根單絲，中間層為9根單絲，外層為15根單絲，單絲直徑均為0.22mm，外繞絲直徑為0.15mm。圖5所示簾線為一種雙螺旋結(jié)構(gòu)鋼簾線，簾線由3股鋼絲螺旋捻合組成，每股又由4根單絲螺旋捻合組成，單絲直徑均為0.22mm，其表述形式為3×4×0.22。鋼絲簾線指定層的單螺旋中心線均可由以下參數(shù)方程描述：

式中，Rs為單螺旋中心線的半徑；αs為捻角；θs為旋轉(zhuǎn)角；θs0為初相角。

對于三層鋼絲簾線，單一螺旋線結(jié)構(gòu)足以描述其幾何形狀。對于雙螺旋結(jié)構(gòu)簾線，以局部坐標(biāo)軸為基礎(chǔ)的數(shù)學(xué)方程能更方便描述鋼絲的幾何特性。對于圖6所示的螺旋結(jié)構(gòu)，定義局部Frenet－Serret軸t、n、b（分別為沿著螺旋線的切線、法線和副法線的單位向量）［28-29］：

圖4 三層結(jié)構(gòu)鋼絲簾線

圖5 雙螺旋結(jié)構(gòu)多股鋼絲簾線

圖6 螺旋結(jié)構(gòu)的幾何特性描述

整體坐標(biāo)系到局部Frenet－Serret坐標(biāo)系的坐標(biāo)轉(zhuǎn)換矩陣可以根據(jù)簾線參數(shù)寫為

則局部Frenet－Serret坐標(biāo)系到整體坐標(biāo)系的坐標(biāo)轉(zhuǎn)換矩陣對a求逆即可：

局部坐標(biāo)系中，矢量r轉(zhuǎn)化為整體坐標(biāo)系：

不同層旋轉(zhuǎn)角度關(guān)系如下：

式中，αs（i）為第i層的捻角；θs（i）為第i層的旋轉(zhuǎn)角度；Rs（i）為第i層的單螺旋曲線的中心線半徑。

雙螺旋中，股的旋轉(zhuǎn)角度和單絲旋轉(zhuǎn)角度之間的關(guān)系為

式中，θw（i，j）為第i層的第j子層的旋轉(zhuǎn)角度；αwt（i，j）為第i層的第j子層的捻角。

基于上述雙螺旋簾線結(jié)構(gòu)的數(shù)學(xué)描述，設(shè)計一個遞歸程序來獲得復(fù)雜簾線的幾何形狀，其算法如下：

（1）每個股層中股的初相角由下式?jīng)Q定：

式中，θs0（k）為第i層第k（k＝1，2，…，ks）根單絲的初相角；ks是第i層中的單絲總數(shù)；θw0（K）為在第i層的第j子層的第K（K＝1，2，…，kw）線的初相角；kw為在第i層的第j子層的鋼絲總數(shù)。

（2）層和子層的旋轉(zhuǎn)角度增量。第k層股的增量旋轉(zhuǎn)角度為

式中，ls（1）為第1層股的長度；lcord為模型中簾線的長度；nrope為沿簾線中心線方向的節(jié)點數(shù)量。

單絲旋轉(zhuǎn)角度為

式中，Δθw（i，j）為第i層的第j子層的增量旋轉(zhuǎn)角度。

（3）單螺旋中心線上點的坐標(biāo)可以用式（4）進(jìn)行計算。

（4）單螺旋表面點的計算方法如下：

為了推導(dǎo)式（15），首先將從單螺旋中心線到表面點矢量寫成

將式（8）的坐標(biāo)轉(zhuǎn)換關(guān)系應(yīng)用到式（16）就可以推導(dǎo)出式（15），得到單螺旋表面點坐標(biāo)的解析描述。

（5）對雙螺旋結(jié)構(gòu)，坐標(biāo)變換公式必須使用兩次。

雙螺旋結(jié)構(gòu)的中心線可描述為

式中，Rwt為股中心到單絲中心的距離。

雙螺旋結(jié)構(gòu)的表面節(jié)點的計算方法如下：

重復(fù)上述步驟（1）～步驟（5）（式（11）～式（19）），可以得到任意復(fù)雜的多股簾線精確數(shù)學(xué)描述。

2.2 有限元網(wǎng)格劃分

基于上述工作，筆者開發(fā)了專用程序Cablemesh來建立多股鋼簾線有限元模型，輸入?yún)?shù)包括簾線設(shè)計變量（捻距、捻角、螺旋中心線半徑）和網(wǎng)格參數(shù)（簾線模型長度，沿軸向、徑向和圓周方向的網(wǎng)格密度）。建模過程中，所有的單絲表面將產(chǎn)生接觸。圖7所示為典型的雙螺旋鋼簾線有限元模型。

圖7 雙螺旋鋼絲簾線有限元模型

3 簾線拉力仿真

研究了上述兩種型號鋼簾線（簾線長度均為20mm）。仿真中，鋼絲本體材料采用彈塑性本構(gòu)模型，材料模型參數(shù)用拉伸試驗方法得到，試驗在Zwick材料試驗機(jī)上進(jìn)行。為正確仿真單絲的受力，有限元模型的單元必須保證精確，一個基本的準(zhǔn)則是在單絲的徑向至少劃分3個單元，而在單絲的周向至少劃分20個單元。根據(jù)這一準(zhǔn)則三層鋼簾線模型被劃分為168 000個單元，而雙螺旋模型被劃分為72 000個單元，模型考慮了所有可能接觸的單絲對，計算在32核CPU并行機(jī)上進(jìn)行。

圖8、圖9分別給出了兩種簾線的測試與有限元仿真應(yīng)力－應(yīng)變曲線。曲線表明：有限元仿真和試驗數(shù)據(jù)吻合得很好；雙螺旋簾線是一種高伸長簾線，仿真和試驗都可以清楚地觀察到在拉伸過程中分別產(chǎn)生結(jié)構(gòu)變形、彈性變形和彈塑性變形（圖8）。三層鋼簾線是一種高強度鋼簾線，其結(jié)構(gòu)變形階段不明顯（圖9）。產(chǎn)生差異的主要原因在于兩種簾線的結(jié)構(gòu)設(shè)計不同，相比于三層鋼簾線，雙螺旋簾線有相對較小的捻距和捻角，在拉伸過程中要經(jīng)歷較大的結(jié)構(gòu)變形后，才能在絲與絲之間形成穩(wěn)定的接觸平衡；兩種型號簾線的基本特征都能夠用有限元來模擬，從而證明了模型的有效性和合理性。

圖8 雙螺旋簾線的拉力－應(yīng)變曲線

圖9 三層簾線的拉力－應(yīng)變曲線

4 拉伸與彎曲組合分析

對長為20mm的三層鋼絲簾線進(jìn)行拉彎組合變形仿真。在簾線末端截面中心施加Δθ轉(zhuǎn)角和張力，產(chǎn)生拉伸－彎曲組合變形，使簾線緊貼在半徑為40mm滾筒上，張力大小分別為50、100、200、300N。

旋轉(zhuǎn)角度Δθ為圓筒上變形簾線的角度跨度，可表示為

式中，L為簾線長度，L＝20mm；ε為等效拉伸應(yīng)變；r為圓筒半徑，r＝40mm；Rs為鋼簾線半徑，Rs＝0.11mm。

假設(shè)單絲獨立彎曲，簾線中鋼絲的應(yīng)力水平為［25-27］

式中，Ew為鋼絲的彈性模量；Rw為單絲半徑。

根據(jù)式（21），計算得出三層中單絲的彎曲應(yīng)力水平為523MPa；采用有限元仿真，三層模型中單絲最大拉伸應(yīng)力為470MPa，最大壓縮應(yīng)力為484MPa。圖10給出了離加載末端10mm處截面上的Von Mises應(yīng)力分布，從中可以得出以下結(jié)論：①單絲在拉伸與彎曲組合載荷作用下均獨立彎曲。②彎曲中性軸的幾何位置取決于簾線所在層的位置和拉力水平。最內(nèi)層單絲中性軸偏移到壓縮區(qū)，最外層單絲中性軸幾乎仍在幾何中軸。③單絲之間的接觸顯著地改變了局部的應(yīng)力狀態(tài)，單絲之間接觸面的增加極大地提高了單絲局部壓力。

5 輪胎鋼簾線應(yīng)變的測量

5.1 應(yīng)變計系統(tǒng)設(shè)計

考慮到 1／4橋式傳感器法［30－31］與鋼條法［31］均有不足，筆者提出一種測量全鋼載重子午線輪胎簾線受力的新方法——焊接法。如圖11所示，該方法在簾線特定的位置上焊錫，再將焊料磨成方形柱狀，然后在方柱的表面粘合應(yīng)變計。簾線直徑約為1.7mm，輪胎制造過程中的最高溫度約為185℃，故選用QFLK－1－11－1LE應(yīng)變計及60通道靜態(tài)應(yīng)變儀。

5.2 拉伸曲線標(biāo)定

鋼簾線埋膠硫化后的試樣如圖12所示。比較用Zwick試驗機(jī)測量埋膠硫化試樣得到的拉伸力－應(yīng)變曲線和從應(yīng)變計得到的拉伸力－應(yīng)變曲線，可以發(fā)現(xiàn)二者之間有一定的差別（圖13）。所以在應(yīng)用應(yīng)變計的數(shù)據(jù)之前，要對其進(jìn)行標(biāo)定，即用Zwick測量數(shù)據(jù)對應(yīng)變計測量得到的曲線特性進(jìn)行標(biāo)定。

圖10 不同拉力下三層簾線離末端10mm處截面的Mises應(yīng)力分布

圖11 簾線應(yīng)變測量實驗中試樣的焊接方法

圖12 鋼簾線埋膠硫化后的試樣

5.3 樣胎制備與簾線受力測量

用粘貼好應(yīng)變片的簾線替換帶束層、胎體層等骨架材料中選定位置的簾線，在成形鼓上完成輪胎制造并硫化成形，得到測試用樣胎。測量輪胎內(nèi)簾線的受力時，充氣壓力為840kPa，將埋于輪胎淺表層的導(dǎo)線引出，與試驗臺、應(yīng)變儀一起搭建簾線受力測量平臺，見圖14。

圖13 埋膠硫化后試樣的拉力－應(yīng)變曲線

圖14 TBR簾線應(yīng)變的測量平臺

5.4 測量結(jié)果與討論

如圖15所示，簾線上的字母a～k表示在該點的應(yīng)變計能夠測到實驗數(shù)據(jù)，有少量的應(yīng)變計無法測到試驗數(shù)據(jù)。圖16顯示的是已測的試驗數(shù)據(jù)和有限元分析結(jié)果，可以發(fā)現(xiàn)試驗和仿真較一致。

圖15 有效應(yīng)變計的位置

圖16 充氣狀態(tài)下輪胎簾線受力試驗與有限元仿真的結(jié)果

6 結(jié)語

結(jié)合已有的整體尺度模型，筆者提出了一種全面細(xì)致分析載重子午輪胎鋼絲簾線變形和應(yīng)力的雙尺度建模方法。開發(fā)的局部尺度模型須以整體尺度模型得到的簾線典型彎曲變形和拉力分布為輸入。研究結(jié)果表明：對于雙螺旋（高伸長率）簾線，在拉伸仿真和實驗中能清楚地觀察到結(jié)構(gòu)變形、彈性變形及彈塑性變形階段，而對三層（高強度）簾線，在拉伸仿真和實驗中發(fā)現(xiàn)其結(jié)構(gòu)變形階段不明顯、彈性階段的剛度更大、塑性變形比較?。焕瓘澖M合仿真中，單絲均獨立彎曲，最內(nèi)層單絲中性軸偏移到壓縮區(qū)，最外層單絲中性軸幾乎仍在幾何中軸，單絲之間接觸面的增大極大地提高了單絲局部壓力；與試驗結(jié)果的比較驗證了有限元分析方法的正確性。

［1］Zorowski C F.Mathematical Prediction of Dynamic Tire Behavior［J］.Tire Science and Technology，1973，1（1）：99－117.

［2］De Eskinazi J，Soedel W，Yang T Y.Contact of an Inflated Toroidal Membrane with a Flat Surface as an Approach to the Tire Deflection Problem［J］.Tire Science and Technology，1975，3（1）：43－61.

［3］Kaga H，Okamoto K，Tozawa Y.Stress Analysis of a Tire under Vertical Load by a Finite Element Method［J］.Tire Science and Technology，1977，5（2）：102－118.

［4］De Eskinazi J，Yang T Y，Soedel W.Displacements and Stresses Resulting from Contact of a Steel Belted Radial Tire with a Flat Surface［J］.Tire Science and Technology，1978，6（1）：48－70.

［5］Watanabe Y.Finite Element Model for Analysis of Deformations of Bias－ply Motorcycle Tires Subject to Inflation Pressure［J］.Vehicle System Dynamics，1984，13（3）：113－128.

［6］Kung L E，Soedel W，Yang T Y.Free Vibration of a Pneumatic Tire－wheel Unit Using a Ring on an Elastic Foundation and a Finite Element Model［J］.Journal of Sound and Vibration，1986，107（2）：181－194.

［7］Gall R，Tabaddor F，Robbins D，et al.Some Notes on the Finite Element Analysis of Tires［J］.Tire Science and Technology，1995，23（3）：175－188.

［8］Tan H F，Du X W，Wei Y T，et al.Mechanical Properties of Cord－rubber Composites and Tire Finite Element Analysis［J］.Vehicle System Dynamics，2004，40（S）：161－174.

［9］Shida Z，Koishi M，Kogure T，et al.Rolling Resistance Simulation of Tires Using Static Finite Element Analysis［J］.Tire Science and Technology，1999，27（2）：84－105.

［10］Wei Y T，Nasdala L，Rothert H.Analysis of Tire Rolling Contact Response by REF Model［J］.Tire Science and Technology，2004，32（4）：214－235.

［11］HanY H，Becker E B，F(xiàn)ahrenthold E P，et al.Fatigue Life Prediction for Cord－rubber Composite Tires Using a Global－lcal Finite Element Method［J］.Tire Science and Technology，2004，32（1）：23－40.

［12］Wei Y T，Tian Z H，Du X W.Finite Element Model for The Rolling Loss Prediction and Fracture A－nalysis of Radial Tires［J］.Tire Science and Technology，1999，27（4）：250－276.

［13］Pek O S，Becker E B.Finite Element Method for the Study of Belt Edge Delaminations in Truck Tires［J］.Rubber Chemistry and Technology，2005，78（4）：557－571.

［14］Feng X，Yan X，Wei Y，et al.Analysis of Extension Propagation Process of Interface Crack between Belts of a Radial Tire Using a Finite Element Method［J］.Applied Mathematical Modeling，2004，28（2）：145－162.

［15］Feng X，Yan X，Wei Y，et al.Nonlinear Finite Element Modeling of Delamination Crack Growth Process between Belts of a Radial Tire［J］.Journal of Reinforced Plastics and Composites，2004，23（4）：373－388.

［16］ZhangY，Hazard C.Effects of Tire Properties and Their Interaction with the Ground and Suspension on Vehicle Dynamic Behavior－a Finite Element Approach［J］.Tire Science and Technology，1999，27（4）：227－249.

［17］Rao K，Kumar R，Bohara P.Transient Finite Element Analysis of Tire Dynamic Behavior［J］.Tire Science and Technology，1998，31（2）：104－127.

［18］Kerchman V.Tire－suspension－chassis Dynamics in Rolling over Obstacles for Ride and Harshness Analysis［J］.Tire Science and Technology，2008，36（3）：158－191.

［19］Brinkmeier M，Nackenhorst U.An Approach for Large－scale Gyroscopic Eigenvalue Problems with Application to High－frequency Response of Rolling Tires［J］.Computational Mechanics，2008，41（4）：503－515.

［20］Brinkmeler M，Nackenhorst U，Petersen S，et al.A Finite Element Approach for the Simulation of Tire Rolling Noise［J］.Journal of Sound and Vibration，2008，309（1／2）：20－39.

［21］Mousseau C W，Hulbert G M.An Efficient Tire Model for the Analysis of Spindle Forces Produced by a Tire Impacting Large Obstacles［J］.Computer Methods in Applied Mechanics and Engineering，1996，135（1／2）：15－34.

［22］Mousseau C W，Hulbert G M.The Dynamic Response of Spindle Forces Produced by a Tire Impacting Large Obstacles in a Plane［J］.Journal of Sound and Vibration，1996，195（5）：775－796.

［23］Mousseau C W，Laursen T A，Lidberg M，et al.Vehicle Dynamics Simulations with Coupled Multibody and Finite Element Models［J］.Finite Elements in Analysis and Design，1999，31（4）：295－315.

［24］HKS Corporation.ABAQUS Theory and User Manual V.6.8［M］.Waltham，MA，USA：Dassault Systems，2008.

［25］Bourgois L.Belt Test for the Evaluation of the Fretting Fatigue and Adhesion Behavior of Steel Cord in Rubber［J］ASTM Special Technical Publication，1979（694）：103－109.

［26］Bourgois L.Survey of Mechanical Properties of Steel Cord and Related Test Methods［J］ASTM Special Technical Publication，1979（694）：19－46.

［27］Zhang Z.Estimate of Endurance Limits of Steel Cord under Different Bending Fatigue Conditions［J］.Wire Journal International，2000，33（10）：116－119.

［28］Usabiaga H，Pagalday J M.Analytical Procedure for Modeling Recursively and Wire by Wire Stranded Ropes Subjected to Traction and Torsion Loads［J］.International Journal of Solids and Structures，2008，45（21）：5503－5520.

［29］Ghoreishi S R，Davies P，Cartraud P，et al.Analytical Modeling of Synthetic Fiber Ropes.Part II：A Linear Elastic Model for 1＋6Fibrous Structures［J］.International Journal of Solids and Structures，2007，44（9）：2943－2960.

［30］Lou A Y C，Walter J D.Interlaminar Shear Strain Measurements in Cord－rubber Composites［J］.Rubber Chemistry and Technology，1979，52（4）：792－804.

［31］Cernik B M.Method of Preparing a Steel Cord for the Measurement of Stress Therein：US，3930918［P］.1976－06－01.