基于離散元與有限元耦合的充氣輪胎沙土路面行駛性能仿真方法研究

鄭祖美， 臧孟炎， 曾海洋

(華南理工大學 機械與汽車工程學院, 廣東 廣州 510640)

基于離散元與有限元耦合的充氣輪胎沙土路面行駛性能仿真方法研究

鄭祖美， 臧孟炎， 曾海洋

(華南理工大學 機械與汽車工程學院, 廣東 廣州 510640)

為研究充氣輪胎在沙土路面的行駛性能，提出用于模擬充氣輪胎- 沙土路面行駛性能的離散元與有限元耦合模型，建立205/55R16輪胎的三維有限元模型和沙土路面離散元模型，并以離散元與有限元接觸算法模擬輪胎與路面之間的相互作用。依據被動滑轉原理，采用基于PID控制理論的加載方式為土槽裝置輪胎加載。在自主研發的越野車沙地行駛行為仿真軟件ORV-SAND中實現以上功能，仿真分析兩種胎面對應輪胎的總牽引力、掛鉤牽引力、行駛阻力和輪輞下陷量等行駛參數，以及不同滑轉率對這些參數的影響。結果表明，仿真分析與試驗結果在趨勢上具有較好的一致性，說明提出的充氣輪胎- 沙土路面離散元與有限元耦合模型及仿真方法在研究輪胎與沙土相互作用關系方面具有很大的潛力。

公路運輸； 充氣輪胎； 沙土路面； 行駛性能； 離散元與有限元耦合方法

Abstract: A combined discrete-finite element method is proposed to simulate the interaction between the pneumatic tire and the granular sand, where the continuum based finite element method is used for the pneumatic tire and the discrete element method is used for the free-flowing particles. Interactions between the pneumatic tire and the sand are described by using a contact detection algorithm. A PID-controller model based on automatic control theory is applied on the base of forced-slip principle to ascertain the loading conditions of tire. A corresponding numerical analysis code is implemented into the in-house developed code, named ORV-SAND, to simulate the tire travel performance. Travel parameters of pneumatic tires with different tread patterns (i.e. smooth and block), such as gross tractive force, drawbar pull, running resistance and sinkage, are investigated under different slip conditions. The simulated results are basically corresponded with the experimental results. It is demonstrated that the proposed discrete-finite element method and the numerical method have enormous potential in studying the interaction between the pneumatic tire and the granular sand.

Key words: highway transportation; pneumatic tire; granular sand; travel performance; discrete-finite element method

0 引言

充氣輪胎與沙土路面的相互作用對越野汽車沙土路面的行駛性能有著重要影響。截止目前，試驗方法是越野汽車行駛性能評價的主要方式。然而，試驗方法周期長、成本高，不利于提高設計開發效率。近年來，隨著計算機技術的發展，越來越多的學者采用數值仿真方法研究輪胎與路面的相互作用。常用的數值仿真方法包括有限元法、離散元法和離散元與有限元耦合的方法等。ABAQUS作為專業的有限元分析軟件，廣泛應用在車輪與路面相互作用問題中[1-4]。此時松軟的沙土路面亦采用有限元方法模擬，并遵守相應的Druck-Prager準則或者Mohr-Coulomb準則。盡管有限元方法能夠準確地模擬輪胎變形且具有較高的計算效率，但是它不能模擬車輪行駛過程中顆粒的飛濺。為了描述輪胎- 路面相互作用時沙土路面的流動和變形，很多學者采用了離散元法[5-8]。該方法中，輪胎由具有確定位置關系的離散單元組成，且通過改變離散單元之間的位置關系來模擬不同結構的輪胎，路面由離散元顆粒構成。張銳等[9-10]采用離散元商業軟件PFC3D仿真了越沙步行輪在松軟路面的行駛性能，并研究了輪面曲率半徑對沙地剛性輪沉陷性能的影響。然而，離散單元雖然能夠再現車輪行駛時的顆粒流動，卻難以準確描述輪胎結構、反映輪胎的變形情況。

離散元與有限元耦合方法作為上述兩種方法的結合，采用有限元法模擬輪胎變形，采用離散元法模擬沙土路面流動，采用離散元與有限元耦合接觸算法模擬輪胎和沙地的相互作用，充分利用了兩種方法的優點。Nakashima等[11-14]首先采用二維離散元與有限元耦合方法(DEM-FEM)對剛性車輪沙地行駛行為進行了仿真分析，并開發了相應的計算程序，隨后又改進了加載方式，仿真分析了星球探測車用彈性車輪在火星土壤上的行駛性能。但是二維仿真不能處理輪胎的側向力等問題。Michael等[15]開發了高效的離散元和有限元接觸算法，最先實現了基于三維DEM-FEM的充氣輪胎沙土路面行駛性能仿真分析。Zhao等[16]基于自主開發的三維DEM-FEM軟件ORV-SAND，分析了不同滑轉率下剛性車輪的沙地行駛行為。ORV-SAND是一款基于FORTRAN 90/95語言、具有完全自主知識產權、離散元與有限元耦合的求解工具，已獲得中國版權保護中心頒發的軟件著作權證書。該求解器包含了適用于處理沙地流動問題的離散元計算模塊，適用于計算充氣輪胎變形問題的有限元模塊和適用于處理輪胎與沙地顆粒接觸問題的離散元與有限元耦合模塊。本文在文獻[16]的基礎上開發了輪胎橡膠用Mooney-Rivlin材料模型，依據室內土槽試驗，建立了充氣輪胎沙土路面行駛的DEM-FEM模型，采用基于PID控制的加載方式，分析不同滑轉率下不同胎面結構的行駛性能，并將仿真結果與試驗結果進行對比，以驗證仿真分析方法的有效性。

1 輪胎- 沙土路面DEM-FEM模型

1.1 沙土路面離散元模型

圖1 沙土路面離散元模型Fig.1 Discrete element model of granular sand

土槽的尺寸設置為1 500 mm×500 mm×200 mm. 本文采用球形離散單元模擬沙粒。受計算條件限制，沙粒的半徑定義在5～7 mm之間，沙地離散元模型通過前處理軟件生成，并在重力作用下壓實[17]。圖1為重力壓實后的沙地離散元模型，沙粒數目為90 905.

1.2 充氣輪胎有限元模型

如圖2所示，將205/55R16輪胎模型簡化為只包含胎面、胎側、胎體和輪輞4部分的簡易輪胎模型。其中，胎面、胎側和胎體部分使用Mooney-Rivlin材料模擬，輪輞設為剛體。本文考慮兩種胎面結構的輪胎，分別為光面輪胎(見圖2(a))和塊狀花紋輪胎(見圖2(b))。輪胎使用線性六面體有限單元描述。

圖2 充氣輪胎有限元模型Fig.2 Finite element model of pneumatic tire

該款輪胎在0.2 MPa內壓的自由滾動半徑約為317 mm，斷面寬度為200 mm. 將充氣輪胎放置到壓實之后的沙土路面上，最終生成的DEM-FEM模型(以塊狀花紋輪胎為例)如圖3所示。

圖3 充氣輪胎沙土路面DEM-FEM模型Fig.3 DEM-FEM model of pneumatic tire and granular sand

2 基于PID理論的加載方式

2.1 土槽試驗

本文土槽試驗采用Shinone等[18]的裝置。圖4和圖5分別為土槽試驗裝置的照片和示意圖，該裝置采用被動滑轉方法測試輪胎在沙土路面的行駛性能。

圖4 試驗裝置照片Fig.4 Experimental set-up

圖5 單輪測試系統示意圖Fig.5 Schematic diagram of single wheel tester

如圖4所示，試驗裝置主要包括土槽和單輪測試系統，單輪測試系統用于測試輪胎的沙地行駛性能。

單輪測試系統包括圖5所示的行駛單元和支撐單元兩部分。輪胎安裝在行駛單元的驅動軸上，驅動軸通過正時皮帶與電機輸出軸(見圖4中的電機1)相連，該電機控制輪胎的轉動，輪胎轉矩由安裝在驅動軸上的傳感器測得。當輪胎與地面接觸時，會使行駛單元沿水平方向向前行駛。支撐單元由土槽端部的電機(見圖4中的電機2)帶動，可沿水平方向行駛。支撐單元與行駛單元之間由兩個平行安裝的垂直滑塊連接，使行駛單元可通過滑塊自由上下移動，滑塊處裝有兩個環形稱重傳感器，用來測量它們之間的水平作用力(即掛鉤牽引力)。

試驗時，支撐單元的速度由電機2確定，是固定值。行駛單元的速度由電機1的輸出轉矩、掛鉤牽引力和地面行駛阻力共同決定。地面行駛阻力是隨時間變化的，故行駛單元的速度是隨時間變化的。行駛單元速度的變化會改變支撐單元與行駛單元之間掛鉤牽引力的大小，而掛鉤牽引力反過來又會影響行駛單元的速度。當行駛單元速度大于支撐單元速度時，掛鉤牽引力會降低行駛單元的速度；當行駛單元速度小于支撐單元速度時，掛鉤牽引力會增加行駛單元的速度。故行駛單元的速度不僅隨時間變化，還圍繞著支撐單元的速度上下波動。保持輪胎的轉動角速度不變，改變支撐單元的水平行駛速度，即可實現滑轉率控制。

2.2 加載方式

為了真實地反映試驗過程中輪胎與沙地的相互作用，本文采用PID控制加載方式[14]再現單個輪胎在土槽中的實際運動過程。

如圖6所示，假設在固定滑轉率s下，車輪的轉動角速度為ω，沿x軸正方向的水平行駛速度為v. 需對輪胎加載的水平方向的力F可表示為

(1)

該式進一步可寫為

(2)

圖6 掛鉤牽引力計算示意圖Fig.6 Schematic diagram of drawbar pull calculation

同樣地，如圖7所示，需對輪胎進行轉矩T的加載，可表示為

(3)

圖7 轉矩計算示意圖Fig.7 Schematic diagram of torque calculation

此時，車輪的總牽引力Fw、掛鉤牽引力Fd和行駛阻力Fr可由(4)式～(6)式表示：

(4)

Fd=-F,

(5)

Fr=-(Fw-Fd).

(6)

在此過程中，水平行駛速度與轉動角速度必須滿足(7)式：

(7)

式中：r為滾動半徑。

3 仿真模型定義

3.1 仿真模型相關物性參數定義

仿真模型使用的各部分材料參數見表1，此時沙土路面考慮為干沙。沙土路面與充氣輪胎之間的接觸參數見表2. 表3給出了計算掛鉤牽引力和轉矩過程中用到的比例、積分和微分系數。仿真參數見文獻[14,16,18]。

表1 充氣輪胎- 沙土路面材料參數

表2 充氣輪胎- 沙土路面接觸參數

表3 比例、積分和微分系數

3.2 仿真分析條件定義

分別計算光面輪胎和塊狀花紋輪胎在不同滑轉率下的行駛性能，步驟如下：

1)對輪胎進行充氣，充氣氣壓為0.20 MPa，充氣時間為0～10 ms；

2)對輪胎進行豎直方向載荷加載，載荷大小為980 N，加載時間為10～12 ms；

3)對輪胎進行速度加載，使其轉動角速度和水平行駛速度達到預定值。輪胎的轉動角速度預設值固定為0.005 rad/ms，水平行駛速度隨滑轉率而變化。本文中滑轉率取5%、15%、25%和35%。速度加載時間為12～14 ms；

4)速度達到預定值后解除速度定義(相當于輪胎獲得了初速度)，對輪胎進行掛鉤牽引力和轉矩加載，此后輪胎正常行駛。

輪胎與路面之間的接觸算法使用自主開發的適用于處理顆粒與復雜結構接觸問題的DZcell算法[19]，輪胎與沙粒之間的作用力按照Hertz-Mindlin理論[20]計算。輪胎從圖3所示位置開始沿x軸正方向行駛，記錄行駛速度、法向力、總牽引力、掛鉤牽引力和輪輞下陷量等參數。

4 仿真結果及分析

圖8給出了滑轉率為5%時的輪胎轉動角速度和水平行駛速度。由圖8可以看出，除了加載初始時刻的波動外，輪胎的速度基本保持在預設值附近(轉動角速度為0.005 rad/ms，水平行駛速度為1.5 mm/ms)，誤差控制在5%以內。由此可見，本文使用的基于PID控制理論的加載方式能夠使輪胎穩定行駛。

圖8 滑轉率為5%時輪胎的轉動角速度和水平行駛速度Fig.8 Rotational and translational velocities of tire under the slip ratio of 5%

圖9 充氣輪胎滑轉率為5%時行駛至470 ms的車轍情況Fig.9 Wheel ruts of pneumatic tires at 470 ms under the slip ratio of 5%

圖9給出了兩款輪胎滑轉率為5%、行駛至470 ms時的行駛軌跡，云圖顯示了沙土路面顆粒的豎直方向位移。由圖9可知：輪胎駛過的區域，沙土路面顆粒發生下陷及流動；輪胎兩側及后側路面稍有抬升，前方路面抬升較多，這是因為輪胎行駛接近土槽邊界，該區域的路面受到來自輪胎和剛性墻的共同作用，堆積現象比較明顯；塊狀花紋輪胎駛過的區域留下了較清晰的車轍，路面下陷深度更大，這是因為帶花紋的輪胎對路面的剪切作用更大。

圖10 充氣輪胎滑轉率為35%時的法向力、輪輞下陷量、掛鉤牽引力和總牽引力Fig.10 Normal force, felloe sinkage, drawbar pull and gross tractive force of pneumatic tires under the slip of 35%

圖10給出了兩款輪胎在滑轉率為35%時的法向力Fn、輪輞下陷量df、掛鉤牽引力Fd和總牽引力Fw隨時間變化的曲線圖。圖10中，B表示塊狀花紋輪胎，S表示光面輪胎，Fw,B即為塊狀花紋輪胎的總牽引力，其余以此類推。由圖10可知：輪胎受力在初始時刻有較大的波動，隨后逐漸穩定；兩款輪胎的法向力隨時間變化情況基本一致，都在980 N附近波動；花紋輪胎的總牽引力、掛鉤牽引力和輪輞下陷量均大于光面輪胎。

圖11和圖12給出了仿真和試驗[18]中不同花紋結構的輪胎行駛參數與滑轉率的關系曲線。其中，圖11為總牽引力、掛鉤牽引力和行駛阻力隨滑轉率的變化趨勢。圖12為輪輞下陷量隨滑轉率的變化趨勢。仿真結果取輪胎行駛到300～600 mm之間的平均值。

圖11 總牽引力、掛鉤牽引力和行駛阻力仿真結果與試驗結果對比Fig.11 Comparison of simulated and experimental results: gross tractive force, drawbar pull and running resistance

圖12 輪輞下陷量仿真結果與試驗結果對比Fig.12 Comparison of simulated and experimental results of felloe sinkage

由圖11和圖12可知：

1) 輪胎的總牽引力隨著滑轉率的增大單調增加。

2) 隨著滑轉率的增大，輪胎的掛鉤牽引力增加；當滑轉率大于20%時，掛鉤牽引力的增長速率減小；當滑轉率大于20%后，塊狀花紋輪胎的掛鉤牽引力比光面輪胎大。

3) 當滑轉率大于5%時，輪胎的行駛阻力的絕對值隨著滑轉率的增大而增大。

4) 輪輞下陷量隨著滑轉率的增加而增大。

圖11、圖12中仿真結果與試驗結果相比，雖然二者存在一定的差別，但趨勢基本一致，從而定性說明了仿真分析方法的有效性。產生差異的原因主要是沙粒大小和形狀都與真實沙粒有較大的區別。后續研究將進行非球形離散單元開發，以逼真表現沙粒形狀；開展分析軟件的并行算法開發，使離散元顆粒逐步接近沙粒尺寸，實現輪胎沙地性能的定量評價；還會針對具體的應用，采取不同的沙地樣本作相應的試驗，確定仿真中需要的各項參數與實際沙地相對應。

5 結論

1)本文建立了充氣輪胎沙土路面行駛性能評價的三維DEM-FEM耦合模型，采取了基于PID控制理論的加載方式，重現了室內土槽試驗中輪胎的受力情況，仿真分析了不同花紋胎面結構的充氣輪胎在不同滑轉率下的行駛性能。

2)結果表明，隨著滑轉率增大，輪胎的總牽引力和掛鉤牽引力都增加，但當滑轉率大于20%以后，掛鉤牽引力的增長率會減小。當滑轉率大于20%以后，花紋輪胎的掛鉤牽引力比光面輪胎大。當滑轉率大于5%后，輪胎行駛阻力的絕對值隨著滑轉率的增加而增大。

3)盡管采用簡化的充氣輪胎代替真實輪胎、粒徑較大的球形離散元顆粒模擬非球形沙粒，仍然得到了與試驗在趨勢上基本一致的仿真結果，說明了本文提出的DEM-FEM模型及仿真方法在越野車輛沙地行駛性能評價方面有很大的應用潛力。

References)

[1] 王萌. 基于軟路面上車輛平順性仿真的輪胎土壤接觸模型研究[D]. 南京: 東南大學, 2009. WANG Meng. Tire-soil contact model for vehicle ride comfort simulation on soft ground[D]. Nanjing: Southeast University, 2009. (in Chinese)

[2] 任茂文, 韓卿, 張曉陽. 采用 ABAQUS/Explicit 分析滾動輪胎與變形地面相互作用[J]. 現代制造工程, 2012 (12):40-43. REN Mao-wen, HAN Qing, ZHANG Xiao-yang. The analysis on the interaction between the rolling tire and the deformed ground based on ABAQUS/Explicit[J]. Modern Manufacturing Engineering, 2012(12):40-43. (in Chinese)

[3] Li H, Schindler C. Three-dimensional finite element and analytical modelling of tyre-soil interaction[J]. Proceedings of the Institution of Mechanical Engineers Part K Journal of Multi-body Dynamics, 2013, 227(1):42-60.

[4] Cueto O G, Coronel C E I, Morfa C A R, et al. Three dimensional finite element model of soil compaction caused by agricultural tire traffic[J]. Computers & Electronics in Agriculture, 2013, 99(6):146-152.

[5] 方俊, 閆民, 許立峰. 顆粒輪胎與顆粒地面模型[J]. 科技導報, 2007, 25(24): 69-72. FANG Jun, YAN Min, XU Li-feng.Model of granule tyre and granule ground[J]. Science & Technology Review, 2007, 25(24): 69-72. (in Chinese)

[6] 葉克東. 三維 GVT 動力學模型研究[D]. 北京：北京林業大學, 2011. YE Ke-dong.Research of dynamics model of 3D GVT[D]. Beijing: Beijing Forestry University, 2011. (in Chinese)

[7] Khot L R, salokhe V M, Jayasuriya P W, et al. Experimental va-lidation of distinct element simulation for dynamic wheel-soil interaction[J]. Journal of Terramechanics, 2007, 44(6):429-437.

[8] Du Y, Gao J, Jiang L, et al. Numerical analysis of lug effects on tractive performance of off-road wheel by DEM[J]. Journal of the Brazilian Society of Mechanical Sciences & Engineering, 2016,39(6):1-11.

[9] 張銳,吉巧麗,張四華,等. 越沙步行輪仿生設計及動力學性能仿真[J]. 農業工程學報，2016,32(15):26-31. ZHANG Rui, JI Qiao-li, ZHANG Si-hua, et al. Bionic design and dynamics performance simulation of walking wheel to travel on sand[J]. Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(15): 26-31. (in Chinese)

[10] 張銳,吉巧麗,張四華,等. 輪面曲率半徑對沙地剛性輪沉陷性能影響研究[J]. 農業機械學報，2016,47(11):341-349. ZHANG Rui, JI Qiao-li, ZHANG Si-hua, et al. Effect of wheel surface curvature radius on sinkage performance of sand rigid wheel[J]. Transactions of the Chinese Society for Agricultural Machinery, 2016,47(11):341-349. (in Chinese)

[11] Nakashima H, Oida A. Algorithm and implementation of soil-tire contact analysis code based on dynamic FE-DE method[J]. Journal of Terramechanics, 2004, 41(2):127-137.

[12] Nakahima H, Takatsu Y, ShinoneH, et al. FE-DEM analysis of the effect of tread pattern on the tractive performance of tires operating on sand[J]. Journal of Mechanical Systems for Transportation & Logistics, 2009, 2(1):55-65.

[13] Ono T, Nakashima H, Shimizu H, et al. Analysis of elastic wheel performance for off-road mobile robots using FE-DEM[J]. IFAC Proceedings Volumes, 2010, 43(26): 61-66.

[14] Nishiyama K, Nakashima H, Yoshida T, et al. 2D FE-DEM analysis of tractive performance of an elastic wheel for planetary rovers[J]. Journal of Terramechanics, 2016, 64:23-35.

[15] Michael M, Vogel F, Peters B. DEM-FEM coupling simulations of the interactions between a tire tread and granular terrain[J]. Computer Methods in Applied Mechanics & Engineering, 2015, 289:227-248.

[16] Zhao C, Zang M. Analysis of rigid tire traction performance on a sandy soil by 3D finite element-discrete element method[J]. Journal of Terramechanics, 2014, 55(7):29-37.

[17] 趙春來. 基于DEM/FEM的越野車沙地行駛行為仿真評價方法研究[D]. 廣州：華南理工大學, 2015. ZHAO Chun-lai. The research on the simulation evaluation method for off-road vehicle driving behavior on sand based on the DEM/FEM [D]. Guangzhou: South China University of Technology, 2015. (in Chinese)

[18] Shinone H, Nakashima H, Takatsu Y, et al. Experimental analysis of tread pattern effects on tire tractive performance on sand using an indoor traction measurement system with forced-slip mechanism[J]. Engineering in Agriculture Environment & Food, 2010, 3(2):61-66.

[19] Zheng Z, Zang M, Chen S, et al. An improved 3D DEM-FEM contact detection algorithm for the interaction simulations between particles and structures[J]. Powder Technology, 2016, 305:308-322.

AnalysisofTravelPerformanceofPneumaticTireonUnpavedRoadbyDiscrete-finiteElementMethod

ZHENG Zu-mei， ZANG Meng-yan， ZENG Hai-yang

(School of Mechanical & Automotive Engineering, South China University of Technology, Guangzhou 510640, Guangdong, China)

U461.5+1

A

1000-1093(2017)09-1822-08

10.3969/j.issn.1000-1093.2017.09.020

2017-01-11

國家自然科學基金項目(11672344)

鄭祖美(1990—), 女, 博士研究生。E-mail: z.zumei@mail.scut.edu.cn

臧孟炎(1961—), 男,教授,博士生導師。E-mail: myzang@.scut.edu.cn