趙春來(lái)，臧孟炎
（華南理工大學(xué)機(jī)械與汽車工程學(xué)院 廣州）
摘 要：采用離散元與有限元耦合方法（FEM/DEM）模擬輪胎沙石路面相互作用，提出 “路面更替法”，以控制仿真模型規(guī)模，提高計(jì)算效率。建立剛性有限元輪胎在離散元沙石路面上行駛的三維仿真模型，使用自主開發(fā)軟件CDFP對(duì)剛性車輪砂石路面行駛行為進(jìn)行仿真分析。獲取了輪胎掛鉤牽引力、下陷量、沙石顆粒流動(dòng)趨勢(shì)等行駛行為參數(shù)。通過(guò)與相關(guān)實(shí)驗(yàn)結(jié)果對(duì)比，驗(yàn)證了該方法的可行性。
關(guān)鍵詞：越野車輛；通過(guò)性能；沙石路面；FEM/DEM；路面更替
中圖分類號(hào)：U461.5+1 文獻(xiàn)標(biāo)志碼：A
深入研究輪胎-沙粒間的相互作用對(duì)高水平越野車的設(shè)計(jì)、性能評(píng)價(jià)等具有重要的指導(dǎo)意義。然而，單純使用有限元方法（FEM）雖然能夠很好地描述輪胎特性[1-5]（包括輪胎的復(fù)雜花紋形狀、充氣輪胎的變形等），但在研究越野車輛沙地行駛行為評(píng)價(jià)方面遇到了挑戰(zhàn)：由于沙石路面具有典型的非連續(xù)介質(zhì)特性，用傳統(tǒng)的連續(xù)介質(zhì)理論（如：FEM）不能合理描述其運(yùn)動(dòng)行為。離散元方法（DEM）在描述散體介質(zhì)的非連續(xù)特性方面具有明顯優(yōu)勢(shì)[6-9]。
為了兼顧FEM和DEM的優(yōu)點(diǎn)，Nakashima[10-12]等采用二維離散元與有限元耦合方法（FEM/DEM）研究輪胎在沙石路面上的行駛行為，并開發(fā)了相應(yīng)的計(jì)算程序。但是，二維方法在研究輪胎行駛行為時(shí)存在很多局限性，例如：無(wú)法描述輪胎花紋、沙石顆粒的側(cè)向流動(dòng)及輪胎側(cè)向力等。本文介紹了一種我們自主開發(fā)的、適用于越野車輪胎行駛行為仿真分析的三維FEM/DEM計(jì)算方法，然后針對(duì)越野車輛沙地行駛過(guò)程的特點(diǎn)，提出了“路面更替法”以控制仿真模型規(guī)模，并基于FORTRAN 95語(yǔ)言開發(fā)了相應(yīng)的專用分析軟件。
1 方法介紹
1.1 FEM/DEM算法介紹
本文所述FEM/DEM方法在車輪-沙地相互作用分析的本質(zhì)是使用FEM來(lái)模擬輪胎，用DEM方法來(lái)描述沙粒。離散單元及與沙石接觸輪胎有限單元節(jié)點(diǎn)的運(yùn)動(dòng)受牛頓第二定律控制。算法中包含圖1（a）、（b）所示的兩種接觸類型：（1）離散單元間的接觸；（2）離散單元與有限單元間的接觸。這兩種接觸類型分別用來(lái)描述沙石顆粒間的相互作用和輪胎與沙石顆粒間的相互作用。其中前者采用Williams等提出的‘C-grid’方法[13]進(jìn)行接觸判斷，后者采用作者開發(fā)的適用于球形離散單元和六面體有限單元的NTS算法[14]。
（a） 離散單元間相互作用 （b） 離散元與有限元相互作用 （c） 單元間作用力
圖1 單元接觸模型
圖1中，hij為單元間的重疊量，vi、vj、ωi、ωj分別為單元i和j的速度和角速度。Fn和Fs分別為單元間的法向力和切向力，如式（1）、（2）所示，其中切向力Fs滿足庫(kù)倫摩擦定律。
（1）
（2）
式中，Fn,e、Fn,v、Fs,e、Fs,v分別為單元間的法向彈簧力、法向阻尼力、切向彈簧力和切向阻尼力，通過(guò)Hertz-Mindlin理論求解[15]。μ為單元間的摩擦系數(shù)。求解有限單元與離散單元間作用力時(shí)，將有限單元視為半徑無(wú)窮大的離散單元[16]。
輪胎與沙石顆粒相互作用的二維示意圖如圖2所示，車輪用有限元描述，沙石路面用離散元描述。輪胎沿
正方向行駛。輪胎的掛鉤牽引力N、輪胎所受的路面豎直方向反作用力P及輪胎滑轉(zhuǎn)率s分別為
（3）
（4）
（5）
其中，G=Σfx+；R=Σfx-；f為輪胎與沙石顆粒間的相互作用力；V為輪胎質(zhì)心水平速度； r為輪胎的滾動(dòng)半徑，對(duì)剛性輪胎，滾動(dòng)半徑等于輪胎自然半徑；
為輪胎的角速度。
圖2 輪胎與路面間的相互作用示意圖
2.2 路面置換法
離散元方法在描述非連續(xù)的沙石路面時(shí)具有優(yōu)勢(shì)，但該法計(jì)算效率較低，而影響計(jì)算效率的主要因素是參與計(jì)算的離散單元的數(shù)量。實(shí)際當(dāng)中，車輪后方一定范圍外，被碾壓過(guò)的沙石顆粒對(duì)車輪行駛行為的影響可以忽略，同時(shí)，車輪前方一定距范圍外，沙石顆粒對(duì)車輪行駛行為影響也較小。基于此，我們提出一種“路面更替法”，以控制參與計(jì)算的離散單元數(shù)目，提高計(jì)算效率。下面通過(guò)二維示意圖描述該方法的主要步驟：
（1）首先，在給定區(qū)域內(nèi)生成離散單元集，作為初始沙石路面樣本，記為DES A0，長(zhǎng)度為a，如圖3（a）所示；
（2）第二，將上一步中的離散單元集DES A0復(fù)制兩份，分別記為DES A1和DES A2，并按順序排列，組成長(zhǎng)度為2a的初始路面，如圖3（b）所示；
（3）第三，將車輪放置在路面DES A1中間位置，車輪在自重及外載荷作用下自由下沉，直到輪胎所受的路面法向反作用力等于車輪自重及外載荷之和。然后給車輪施加恒定的角速度和相應(yīng)的水平速度，使車輪在某一滑轉(zhuǎn)率下行駛，如圖3（c）所示；
（4）第四，當(dāng)車輪行駛到DES A2的中間位置時(shí)，DES A1對(duì)車輪行駛行為影響可以忽略，且DES A2前端沒(méi)有發(fā)生明顯破壞，此時(shí)開始進(jìn)行第一次路面更替：將DES A1刪除，同時(shí)重新備份路面樣本，記為DES A3，放置在DES A2前端，與DES A2組成新的路面，如圖3（d）所示。
反復(fù)執(zhí)行步驟（4），參與計(jì)算的沙石路面長(zhǎng)度始終等于2a，可以實(shí)現(xiàn)較少數(shù)量的離散單元來(lái)模擬任意長(zhǎng)度的沙石路面，大幅提高計(jì)算效率的目的。對(duì)于由多種結(jié)構(gòu)路面構(gòu)成的越野車輛沙地通過(guò)性能評(píng)價(jià)，只需對(duì)每種結(jié)構(gòu)路面生成樣本路面，然后根據(jù)車輪即將行駛的路面結(jié)構(gòu)調(diào)入對(duì)應(yīng)結(jié)構(gòu)樣本路面即可。
（a） 沙石區(qū)域樣本 （b） 組成初始沙石區(qū)域
（c） 放置輪胎 （d） 路面更替
圖3 路面置換方法的二維示意圖
同種路面結(jié)構(gòu)相應(yīng)的程序流程圖如圖4所示，其中Tter為計(jì)算終止時(shí)間，Δt為顯式算法的時(shí)間步長(zhǎng)，T為當(dāng)前計(jì)算時(shí)間。
圖4 執(zhí)行流程圖
3 數(shù)值算例
為驗(yàn)證方法的有效性，基于室內(nèi)土槽實(shí)驗(yàn)[17]，建立了輪胎在沙石路面行駛的三維仿真模型，分析了不同滑轉(zhuǎn)率下輪胎的行駛行為。
3.1 離散元沙石路面
首先，采用作者之前開發(fā)的分層生成方法[18]，在指定區(qū)域內(nèi)，隨機(jī)生成給定半徑范圍的球形離散單元集，區(qū)域邊界采用剛性墻約束。然后，讓初始生成的離散單元在自重作用下達(dá)到穩(wěn)定狀態(tài)，以模擬自然沙地。最終生成的離散元在X: （0,735）、Y: （0,480）、Z: （0,280）范圍，離散單元集孔隙率約為0.32，以此作為沙石路面初始樣本，如圖5所示。本文側(cè)重于對(duì)開發(fā)方法的驗(yàn)證，因此，離散單元半徑范圍取6~7mm，大于實(shí)際沙粒半徑。自重壓實(shí)過(guò)程離散單元總勢(shì)能的時(shí)間歷程如圖6所示，在1.1s以后，其趨勢(shì)基本平穩(wěn)，此時(shí)可以認(rèn)為壓實(shí)過(guò)程結(jié)束。
圖5 初始沙石路面樣本 圖6 離散元沙石顆粒重力勢(shì)能時(shí)間歷程
楊氏模量：75000 MPa，泊松比：0.3，密度：2400 kg/m3，單元數(shù)目：45551，摩擦系數(shù)：0.3
3.2 初始路面及車輪模型
初始沙石路面區(qū)域由兩段經(jīng)自重壓實(shí)后的“路面樣本”組成，單元總數(shù)為91102個(gè)，長(zhǎng)度等于路面樣本長(zhǎng)度二倍的1470mm。有限元車輪采用如圖7所示的簡(jiǎn)化模型，暫不考慮輪胎變形。輪胎與沙路間的摩擦系數(shù)為0.4，輪胎放置在如圖所示位置。給輪胎施加自重及外載荷（計(jì)1295N）同時(shí)加載恒定的角速度0.5rad/s和相應(yīng)的水平速度，使輪胎在給定滑轉(zhuǎn)率下沿X軸正方向行駛。
（a） 主視圖 （b） 左視圖
圖7 初始路面和有限元輪胎模型
輪胎參數(shù)：楊氏模量：2 MPa，泊松比：0.49，密度：1800 kg/m3，單元數(shù)目：1344
3.3 仿真結(jié)果
當(dāng)車輪行駛距離為k×735 mm（k為路面更替次數(shù)）時(shí)執(zhí)行路面更替。圖8給出輪胎在30%滑轉(zhuǎn)率下，離散單元Z方向的位移云圖：圖8（a）為初始狀態(tài)；圖8（b）為輪胎行駛到第二段樣本路面中間位置時(shí)的車轍；圖（c）為執(zhí)行第一次沙石路面更替后的車轍；圖8（d）給出了第二路面更替前的車轍狀況；圖8（e）執(zhí)行了第二次路面更替；圖8（f）為計(jì)算終止?fàn)顟B(tài)，車輪行駛距離為1560mm。
（a） （b）
（c） （d）
（e） （f）
圖8 車輪滑轉(zhuǎn)率為30%時(shí)沙石顆粒Z方向位移云圖
圖9為輪胎在滑轉(zhuǎn)率為30%時(shí)的輪胎掛鉤牽引力時(shí)間歷程。掛鉤牽引力在初始階段有一個(gè)較大峰值，主要原因是初始階段，輪胎與路面間豎直方向有相對(duì)沖擊速度，產(chǎn)生一個(gè)沖擊力，引起路面法向反作用力出現(xiàn)峰值，與此相應(yīng)，掛鉤牽引力初始階段也出現(xiàn)峰值。隨后掛鉤牽引力快速減小，但仍有較大波動(dòng)，主要原因是離散單元顆粒半徑取值較大，對(duì)接觸力計(jì)算結(jié)果影響所致。
圖9 滑轉(zhuǎn)率為30%時(shí)車輪行駛過(guò)程 圖10 滑轉(zhuǎn)率為30%時(shí)車輪行駛過(guò)程中的
中的輪胎掛鉤牽引力時(shí)間歷程 輪胎下陷量時(shí)間歷程
圖10為輪胎行駛過(guò)程中的下陷量時(shí)間歷程，由圖可知：在初始階段，輪胎下陷量快速增加，主要原因是輪胎豎直方向載荷未達(dá)到平衡狀態(tài)。在行駛200mm之后，下陷量基本上趨于一個(gè)穩(wěn)定的狀態(tài)，均值約為45mm，且在路面更替時(shí)未出現(xiàn)明顯波動(dòng)。
圖11為輪胎在30%滑轉(zhuǎn)率下行駛時(shí)，沙石顆粒的流動(dòng)趨勢(shì)，圖中箭頭表示離散單元在X-Z平面內(nèi)的速度矢量。由圖可知：沙石顆粒的流動(dòng)趨勢(shì)可分為兩個(gè)區(qū)域：前區(qū)，繞順時(shí)針?lè)较颍缓髤^(qū)，繞逆時(shí)針?lè)较颉Ｔ摤F(xiàn)象與莊繼德[19]等人的實(shí)驗(yàn)研究結(jié)果基本一致。
（a） 仿真結(jié)果 （b） 實(shí)驗(yàn)結(jié)果
圖11 車輪行駛過(guò)程中沙石顆粒的運(yùn)動(dòng)趨勢(shì)
滑轉(zhuǎn)率是影響輪胎行駛行為的主要因素，為研究不同滑轉(zhuǎn)率下輪胎的行駛行為，圖12和圖13分別給出了輪胎下陷量和掛鉤牽引力隨滑轉(zhuǎn)率的變化關(guān)系，其中下陷量和掛鉤牽引力值取各工況穩(wěn)定狀態(tài)下的均值。由圖可知：輪胎的下陷量隨滑轉(zhuǎn)率的增加而逐漸增加，而且其增加趨勢(shì)逐漸“陡峭”，該趨勢(shì)與文獻(xiàn)[17]研究結(jié)果基本一致。
掛鉤牽引力隨滑轉(zhuǎn)率的增加逐漸增加，當(dāng)滑轉(zhuǎn)率小于25%時(shí)，增加趨勢(shì)較陡峭，當(dāng)滑轉(zhuǎn)率大于25%時(shí)，增加趨勢(shì)緩慢，該現(xiàn)象也與文獻(xiàn)[17]實(shí)驗(yàn)結(jié)果整體趨勢(shì)一致。
圖12 輪胎下陷量隨滑轉(zhuǎn)率變化關(guān)系 圖13 掛鉤牽引力隨滑轉(zhuǎn)率變化關(guān)系
4 結(jié) 論
將FEM/DEM方法應(yīng)用于車輪在沙石路面的行駛行為仿真分析，并提出了一種“路面更替法”。路面更替法可以實(shí)現(xiàn)一定數(shù)量的離散單元來(lái)模擬任意長(zhǎng)度的沙石路面，達(dá)到控制模型規(guī)模，提高計(jì)算效率的目的。
輪胎行駛行為，包括路面豎直方向反作用力，掛鉤牽引力，下陷量及輪胎行駛過(guò)程中沙石顆粒的流動(dòng)狀況等都可很方便地通過(guò)計(jì)算求得，數(shù)值試驗(yàn)結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比驗(yàn)證了該方法的有效性。
參考文獻(xiàn)
[1] Xia K. Finite element modeling of tire/terrain interaction: Application to predicting soil compaction and tire mobility[J]. Journal of Terramechanics. 2011;48:113-23.
[2] Xia K, Yang YM. Three-dimensional finite element modeling of tire/ground interaction[J]. International Journal For Numerical and Analytical Methods in Geomechanics. 2012;36:498-516.
[3] Li H, Schindler C. Analysis of soil compaction and tire mobility with finite element method[J]. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics. 2013;227:275-91.
[4] González Cueto O, Iglesias Coronel CE, Recarey Morfa CA, Urriolagoitia Sosa G, Hernández Gómez LH, Urriolagoitia Calderón G, et al. Three dimensional finite element model of soil compaction caused by agricultural tire traffic[J]. Computers and Electronics in Agriculture. 2013;99:146-52.
[5] 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:42-60.
[6] Smith W, Peng H. Modeling of wheel-soil interaction over rough terrain using the discrete element method[J]. Journal of Terramechanics. 2013;50:277-87.
[7] Li W, Huang Y, Cui Y, Dong S, Wang J. Trafficability analysis of lunar mare terrain by means of the discrete element method for wheeled rover locomotion[J]. Journal of Terramechanics. 2010;47:161-72.
[8] 張銳, 劉芳, 曾桂銀, 李建橋. 非規(guī)則車輪 散體模擬月壤相互作用仿真系統(tǒng)研究[C]. 顆粒材料計(jì)算力學(xué)研究進(jìn)展2012.
[9] Jiang M, Liu F, Shen Z, Zheng M. Distinct element simulation of lugged wheel performance under extraterrestrial environmental effects[J]. Acta Astronautica. 2014;99:37-51.
[10] Nakashima H, Takatsu Y, Shinone H. Analysis of tire tractive performance on deformable terrain by finite element-discrete element method[J]. Journal of Computational Science and Technology. 2008;4:423-34.
[11] Nakashima H, Takatsu Y, Shinone H. 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 and Logistics. 2009;2:55-65.
[12] 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:127-37.
[13] Williams J, Perkins E, Cook B. A contact algorithm for partitioning N arbitrary sized objects[J]. Engineering Computations. 2004;21:235-48.
[14] Zang MY, Gao W, Lei Z. A contact algorithm for 3D discrete and finite element contact problems based on penalty function method[J]. Computational Mechanics. 2011;48:541-50.
[15] Balevi?ius R, D?iugys A, Ka?ianauskas R. Discrete element method and its application to the analysis of penetration into granular media[J]. Journal of Civil Engineering and Management. 2004;1:3-14.
[16] Han K, Peric D, Crook A, Owen D. A combined finite/discrete element simulation of shot peening processes–Part I: studies on 2D interaction laws[J]. Engineering Computations. 2000;17:593-620.
[17] Shinone H, Nakashima H, Takatsu Y. 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 and Food. 2010;3:61-6.
[18] Zang M, Zhao C. Numerical Simulation of Rigid Wheel Running Behavior on Sand Terrain[C]. 2013.
[19] 莊繼德. 計(jì)算汽車地面力學(xué)[M]: 機(jī)械工業(yè)出版社; 2002.
