不同流固耦合方法進行輪胎滑水性能仿真研究的對比分析

苑 陽，劉從臻*，錢 浩

（1.山東理工大學 交通與車輛工程學院，山東 淄博 255049；2.江蘇大學 汽車與交通工程學院，江蘇 鎮江 212013）

流固耦合方法涉及到流體與固體間的相互影響，其對流體域和固體域分別求解的特性可更加真實地模擬滑水現象的規律。近年來其在模擬輪胎滑水及深入解析滑水性能方面得到廣泛應用。董斌等[1-3]使用混合花紋輪胎仿真模型，采用大型有限元仿真軟件Fluent計算獲得了不同工況下胎面所受的水流動壓力及胎面接觸區域的水流速度分布。王國林等[4-6]采用計算流體力學（CFD）的氣-液二相流數值模型將輪胎滑水過程簡化為流體動壓力作用使輪胎產生附加變形，并將臨界滑水速度仿真預測值與已有的經驗公式計算值進行對比分析。

從計算原理上看，單向耦合算法僅建立在水流壓力場對輪胎結構場影響的情況下，與實際滑水現象有差異，從本質上存在誤差。

隨著多場耦合技術和計算機軟、硬件技術的發展，雙向流固耦合算法的應用逐漸受到各行專業人士和學者的青睞。在輪胎滑水雙向流固耦合計算方面，臧孟炎等[7]建立了基于Ls-Dyna軟件的輪胎滑水模型，通過分析輪胎接地前緣形成的楔狀水層逐漸被轉動的車輪擠壓進胎面排水溝槽的狀態及其沿縱向排水溝槽噴出的狀態等研究不同水膜厚度下車輛臨界滑水速度的變化。H.Grogger等[8]以光面輪胎和縱向花紋輪胎為研究對象，在考慮結構、大變形和與路面相互作用等因素的基礎上，提出一種基于弱耦合流體-結構相互作用的集成軟件工具，建立了有效且準確傳遞基本變量的耦合模塊和處理流體域網格以匹配輪胎外表面的變形模塊，將流體壓力分布與輪胎仿真模型相結合，重復迭代運算，獲得準確的流體壓力分布和上浮力。

目前部分仿真分析方法未考慮輪胎與水層流固耦合的雙向作用，而且很少通過流體域分析輪胎的滑水特性。因此，通過單向、雙向流固耦合方法進行仿真對比，并基于流體域分析輪胎的滑水特性，從而明晰單向耦合算法的不足是很有必要的。

本研究以某廠家提供的185/60R15輪胎為例，基于Workbench平臺的輪胎滑水流固耦合仿真分析方法，建立精確的輪胎滑水水流模型，在相同工況下進行單向、雙向流固耦合算法的仿真對比。

1 理論背景

1.1 流固耦合研究簡介

隨著計算科學以及數值分析方法的不斷發展，流固耦合或流固交互作用研究從20世紀80年代以來受到了廣泛關注，從數據傳遞角度出發，其可以分為單向流固耦合和雙向流固耦合[9]。一般來說，對于大多數耦合作用現象，如果只考慮靜態結構性能，單向流固耦合分析已經足夠，但是如果要考慮結構的動力學特性，雙向流固耦合分析就必不可少。雙向流固耦合同時考慮流體與固體，計算精度高，但是設置復雜，計算周期長。相比之下，單向流固耦合計算所需時間短、資源較少，僅適用于固體結構形變量不足以顯著影響流體域的情況。單向流固耦合方法將流體域計算所得結果作為負荷施加到結構域中，而沒有將結構域分析結果傳遞給流體域，設置相對簡單，計算周期較短。

1.2 水流控制方程

本研究輪胎滑水流固耦合模型中的水流流動滿足質量守恒、能量守恒和動量守恒定律，這3個物理量守恒構成了水流的基本控制方程[10]。

質量守恒方程：

能量守恒方程：

動量守恒方程：

式中ρw——流體密度；

t——時間；

ui，uj——速度張量；

xi，xj——坐標張量；

E——動能；

σw，εw——流體柯西應力、應變；

p——流體靜壓；

τij——應力張量；

ρwgi——重力體積力；

Fi——其他體積力。

1.3 湍流模型

適合湍動能（k）-耗散率（ε）模型的流動類型比較廣泛，包括旋均勻剪切流、自由流（射流和混合層）、腔道流動和邊界層流動。該模型k和ε的方程[10]為

其中C1＝max[0.43，η/（η＋5）]，η＝Sk/ε

式中μ——動力粘性系數；

μt——湍流粘性系數；

Gk——由于平均速度梯度引起的湍動能；

Gb——由于浮力引起的湍動能；

YM——可壓縮湍流脈動膨脹對總耗散率的影響；

η——用戶修正系數；

S——長度；

C1ε，C2，C3ε——在Fluent軟件中為默認值常數，C1ε＝1.44，C2＝1.90，C3ε＝0.09；

σk，σε——k和ε的湍流普朗特數，σk＝1，σε＝1.20。

1.4 滑水經驗公式

一般而言，有限元模型建立及計算結果有效性最合適的驗證方法是與試驗結果進行對比，然而滑水試驗對測試場地條件和設備要求很高，因此本仿真采用滑水經驗公式對輪胎滑水模型的合理性進行驗證。W.B.Horne等[11]通過分析大量輪胎滑水的試驗結果，提出了著名的NASA滑水經驗公式：

式中vp——臨界滑水速度；

P——輪胎充氣壓力。

NASA滑水經驗公式只適合于積水較薄且被水覆蓋的路面狀況，為此，D.Dunlap等[12]在公式（6）的基礎上，提出了滑水速度與胎面花紋深度、胎面寬度和水膜厚度之間關系的滑水預測公式：

DT——胎面花紋深度，一般情況下取值范圍為1.6～10.0 mm；

DW——水膜厚度；

T——胎面寬度。

本仿真分析中，P＝220 kPa，DT＝7 mm，DW＝10 mm，T＝185 mm，代入公式（7）中，計算pv′ 為102.59 km·h-1（28.50 m·s-1）。

2 滑水物理模型

以某廠家提供的185/60R15輪胎為例，通過參考文獻[13]獲得花紋溝排列方向及花紋類型等對滑水性能的影響，進而建立垂直花紋輪胎和路面三維殼體模型，并賦予輪胎和路面材料屬性，定義輪胎-路面接觸算法，進行網格劃分及邊界條件定義；在CFD模塊中建立流體域模型并進行網格劃分和邊界條件定義[3]。

2.1 結構模型與邊界條件設置

結構場計算的對象為輪胎和剛性路面。輪胎胎面與胎側分開建模、劃分網格，共產生64 783個單元和36 483個節點，如圖1所示。

圖1 輪胎網格模型

網格劃分完成后對模型設置邊界條件，邊界條件包括結構負荷和約束，如圖2所示。

圖2 輪胎邊界條件

2.2 流體域模型與邊界條件設置

對分析水層建立六面體模型并與輪胎進行布爾減操作，剔除兩者相交部分，而剩余部分即為流體域模型。然后進行網格劃分和邊界條件定義。在Fluent軟件中對流體域采用補片協調方法的四面體網格劃分，共產生196 350個四面體單元和40 333個節點，如圖3所示。流體域模型的邊界條件定義為速度入口、壓力出口以及壁面3種形式，并且在速度入口施加Profile文件定義水的運動，如圖4所示。

圖3 流體域模型

圖4 流體域邊界條件

2.3 求解器設置

為了使單向、雙向兩種流固耦合算法的計算結果具有可比性，不僅兩種算法的邊界條件設置、網格劃分完全一致，流體域、結構域對應的求解器設置也完全一致。其中，流體域求解器設置項目較多，主要設置項如下：基于壓力傳遞求解，采用可實現的k-ε湍流模型、平順可重塑的動網格和半隱式方法求解壓力耦合方程組。

3 結果與分析

3.1 流體域

t＝1.5 s時的單向、雙向流固耦合流體域的動水壓力分布如圖5所示。

圖5 流體域動水壓力云圖

由圖5可見，兩種方法的壓力峰值位置分布情況基本一致，但在單向流固耦合流體域中接地前端胎面最大動水壓力作用區域比較集中，說明該處水流動能大部分轉化為作用于胎面的壓力能，容易對輪胎產生明顯的徑向抬升作用，而雙向流固耦合模擬輪胎變形改變了周圍水流的壓力場，正、負壓力均降低，使壓力分布相對較為均勻，難以對胎面形成集中有效的徑向抬升作用，從而降低了輪胎發生滑水現象的幾率。

將流體域流固耦合面處的動水壓力等于或大于輪胎所受負荷時定義為輪胎發生完全滑水的時刻。在相同工況下輪胎-干路面法向接觸力及單向、雙向流固耦合動水壓力隨時間的變化曲線如圖6所示。

圖6 在相同工況下輪胎-干路面法向接觸力及單向、雙向流固耦合動水壓力隨時間的變化曲線

由圖6可見，單向流固耦合計算發生完全滑水時刻為t＝1.1 s，而雙向流固耦合計算發生完全滑水時刻為t＝1.2 s。

從圖6可直觀地看出不同耦合方法模擬整個滑水過程流固耦合面處動水壓力變化的異同，動水壓力曲線變化趨勢基本一致，動水壓力變化大致分為3個階段。

（1）充氣加載階段。該階段流速設置為零。

（2）水流沖擊階段。待輪胎充氣加載穩定后，水逐漸加速沖擊輪胎，由于垂直花紋溝槽的排水能力有限，兩種算法的胎面所受動水壓力都呈現迅速增大的趨勢，而雙向流固耦合模擬輪胎變形改變了周圍水流的壓力場，正、負壓力均降低，使相同時刻動水壓力較小，因此該階段單向流固耦合計算的動水壓力先與輪胎所受法向負荷平衡，即單向流固耦合模擬輪胎滑水發生先于雙向流固耦合模擬。

（3）完全滑水階段。該階段輪胎動水壓力大于或等于所受法向負荷，輪胎處于完全滑水狀態。

3.2 輪胎結構性能分析與對比

以接觸面中心為參考點，t＝1.5 s時輪胎的變形云圖如圖7所示。

圖7 輪胎變形云圖

由圖7可見，該時刻單向、雙向流固耦合計算的輪胎最大變形分別為17.652和15.649 mm。輪胎最大變形發生在接地區域前后兩側，并且變形量分布具有對稱性；由于輪胎接地前緣楔形區域的存在，輪胎與流體互相擠壓，使該區域水流出現滯留，從而使胎面中間區域動水壓力高于胎側區域；相同時刻下，雙向流固耦合計算的最大動水壓力和最大變形量小于單向流固耦合計算結果。

為了方便對兩種流固耦合方法的計算結果進行對比，生成相同工況下輪胎變形隨時間的變化曲線，如圖8所示。

由圖8可見，輪胎發生變形也可大致分為3個階段。

圖8 輪胎變形量隨時間的變化曲線

（1）充氣加載階段。該階段兩種算法的邊界條件設置、網格劃分完全一致，結構域對應求解器設置也完全一致，因此輪胎發生的變形量也趨于一致。

（2）水流沖擊階段。該階段是輪胎充氣加載穩定后水逐漸加速沖擊輪胎導致動水壓力升高的過程。兩種耦合方法都表現出輪胎的非線性變形。同時，相同時刻下單向流固耦合計算的變形量逐漸超過雙向流固耦合，因而提前達到變形量的穩態值，即單向流固耦合計算先發生完全滑水現象。

（3）完全滑水階段。該階段輪胎的垂直負荷完全由流體動水壓力來平衡，輪胎的變形量也趨于穩定。

3.3 誤差分析

在流體域模型的邊界條件設定速度入口，并且在速度入口導入水流運動速度與時間的關系文件。通過單向、雙向流固耦合仿真分析得到的輪胎臨界滑水速度分別為25.80和27.50 m·s-1（分別對應t＝1.1 s和t＝1.2 s時刻）。

將仿真分析得到的臨界滑水速度與經驗公式計算值（28.50 m·s-1）進行對比，單向、雙向流固耦合計算臨界滑水速度的相對誤差分別為9.47%和3.51%，由此在一定程度上說明了雙向流固耦合模擬輪胎滑水性能比單向流固耦合模擬更可靠。

4 結論

為了比較輪胎滑水單向、雙向流固耦合兩種方法計算結果的差異，本研究以185/60R15輪胎為例，基于Workbench平臺進行輪胎滑水流固耦合仿真分析，得到如下結論。

（1）采用單向、雙向流固耦合方法計算的動水壓力分布迥異，雙向流固耦合方法計算的表面正、負壓力峰值均降低，壓力分布更均勻。

（2）輪胎最大變形發生在接地區域前后兩側，并且變形量分布具有對稱性；由于輪胎接地前緣楔形區域的存在，輪胎與流體互相擠壓，使該區域水流出現滯留，從而使胎面中間區域動水壓力高于胎側區域；相同時刻下，雙向流固耦合計算的最大動水壓力和最大變形量小于單向流固耦合計算結果。

（3）相對經驗公式計算值，單向、雙向流固耦合仿真預測的臨界滑水速度的相對誤差分別為9.47%和3.51%，可見雙向流固耦合方法的計算結果更真實合理，從原理上提高了計算精度，在一定程度上可說明雙向流固耦合方法模擬輪胎滑水過程的可靠性。