汪怡平，張子一

（1.武漢理工大學現(xiàn)代汽車零部件技術湖北省重點實驗室，湖北 武漢 430070；2.武漢理工大學汽車零部件技術湖北省協(xié)同創(chuàng)新中心，湖北 武漢 430070；3.甘肅建投重工科技有限公司，甘肅省專用汽車工程研究中心，甘肅 蘭州 730050）

1 引言

側風環(huán)境下車輛高速通過橋塔時，橋塔處的繞流和脫落渦會影響車輛的行駛穩(wěn)定性[1]。車輛周圍不斷變化的流場產生的瞬態(tài)氣動力作用在車輛上，會引起車輛產生橫向位移和橫擺角，車輛的穩(wěn)定性惡化，乘員的舒適性與安全性受到嚴重影響[2]。隨著高級駕駛輔助技術的發(fā)展，主動校正和自動轉向受到越來越多的關注。研究側風環(huán)境下車輛氣動力以及動態(tài)響應的變化規(guī)律，將為高級駕駛輔助系統(tǒng)的研發(fā)提供理論依據(jù)[3]。

因此，深入研究側風環(huán)境下車輛通過橋塔時的流場變化、氣動性能和動態(tài)響應，具有重要意義。

對于風-橋路況下車輛的氣動特性與穩(wěn)定性，國內外學者通過道路、風洞試驗和數(shù)值模擬開展了廣泛研究。道路試驗能夠真實地獲取車輛在側風下的動態(tài)響應。但道路測試無法用于早期無原型樣車的評估分析。文獻[4]利用風洞試驗測量橋塔周圍不同位置靜止車輛的氣動力和表面壓力。文獻[5]利用風洞中移動轎車模型模擬轎車通過橋塔尾跡獲得瞬態(tài)氣動力。但風洞試驗成本高且不能提供足夠的流場信息。計算流體力學（CFD）能提高瞬態(tài)三維流場信息成為克服實驗局限性的一個有效選擇。文獻[6]利用CFD 方法獲得側風環(huán)境下通過橋塔的重型卡車的氣動特性。然而，以往對于車輛在側風環(huán)境下通過橋塔的氣動特性研究，車輛運動軌跡是預設的，即將CFD 得到的氣動力加載到動力學模型，從而得到動力響應。但實際行車時，當車輛在側風環(huán)境下通過橋塔時，流場突變引起的氣動力變化將引起橫向位移、橫擺角速度等動態(tài)響應，進而引起氣動力的大小和方向發(fā)生變化，最終再次影響車輛的動力學響應。

因此，側風環(huán)境下車輛通過橋塔區(qū)域時，必須考慮動力學響應和氣動響應之間的相互作用。文獻[7]耦合空氣動力學與6自由度的動力學模型對汽車的氣動特性進行了分析。研究表明車輛姿態(tài)變化會產生不穩(wěn)定的氣動力和力矩，動態(tài)耦合法具有廣泛的適應性，能夠有效研究側風環(huán)境等復雜條件下的車輛氣動特性。

采用CFD與多體系統(tǒng)動力學（MBD）耦合的動態(tài)耦合法。通過改變風壓中心位置，驗證汽車動力學模型的魯棒性。通過對比風洞試驗中不同偏航角下的氣動阻力系數(shù)，驗證CFD的可靠性。深入分析動態(tài)耦合法中流場相關信息，基于氣動力與車輛運動之間相互作用機理研究車輛側風穩(wěn)定性。

2 耦合方法

車輛側風行駛狀況可以視作多個物理量之間的耦合作用：周圍的流場、車輛運動等。按照耦合變量的求解方式，可以將耦合方法分為強耦合、弱耦合。強耦合方法通過單元矩陣或荷載向量把耦合作用構造到控制方程中，然后對控制方程直接求解。弱耦合方法是在每一步內分別對CFD和MBD求解，通過把第一個物理場的結果作為外荷載加于第二個物理場來實現(xiàn)兩個場的耦合，本文的研究均采用弱耦合方式。針對汽車側風穩(wěn)定性的研究，目前主要采用靜態(tài)耦合法，其原理，如圖1所示。

圖1 靜態(tài)耦合法示意圖Fig.1 Static Coupling Method

首先由CFD模型進行瞬態(tài)計算，得到不同時刻的氣動力，然后將數(shù)據(jù)輸入MBD模型，計算得到氣動力作用下的車輛運動參數(shù)，并評價車輛的側風穩(wěn)定性。然而，靜態(tài)耦合法中，車輛的運動是預設的，忽略了氣動力對汽車姿態(tài)的影響。

所以，該方法進行車輛側風穩(wěn)定性分析所得結論與真實情況有一定的差異。針對上述問題，采用動態(tài)耦合法，其原理，如圖2所示。該方法在CFD模型和MBD模型中分別加載側風模型和道路模型，并實時交換氣動力與動力學響應數(shù)據(jù)。整個動態(tài)耦合流程，如圖3所示。

圖2 動態(tài)耦合法示意圖Fig.2 Dynamic Coupling Method

圖3 動態(tài)耦合法仿真流程Fig.3 Flow Chart of the Dynamic Coupling Process

在耦合前，CFD和MBD模型分別單獨計算0.5s，以獲得穩(wěn)定的氣動力和動力學響應。然后開始耦合，實時交換氣動力和運動參數(shù)。計算過程如下：將CFD模型計算得到的側向力和橫擺力矩輸入MBD模型；MBD模型計算車輛動態(tài)響應并將橫擺角速度和側向速度反饋給CFD模型；利用CFD計算車輛在新位置的側向力和橫擺力矩。通過動態(tài)耦合法，CFD模型中車輛氣動力變化與MBD模型中車身姿態(tài)變化得到統(tǒng)一，計算結果更加真實。

2.1 湍流模型

車輛周圍的流場是三維、瞬態(tài)、湍流、不可壓縮的。選取RNG k-ε湍流模型進行瞬態(tài)計算，其控制方程如下：

式中：ρ—流體密度；xi、xj—坐標軸分量；μ—流體粘度；t—時間；ut—時刻的速度；ui—速度分量；k—湍流動能系數(shù)；σ—渦黏性系數(shù)。σk=σ?=1.39是k與ε的湍流普朗特施密特數(shù)，Gk是關于k的經驗公式。E、C1、C2、η為常數(shù)，由實驗確定。

2.2 重疊網(wǎng)格法

采用重疊網(wǎng)格技術實現(xiàn)車輛在橋梁上的運動[8]，其連續(xù)性方程和動量方程如下所示：

式中：t—時間；ρ—流體密度；p—壓強；μ—動力粘度；dS、dV—表面積和控制體積；n—控制體表面指向外面的法向量；u—流速矢量；ug—網(wǎng)格運動速度矢量即MBD 模型中的車輛速度。網(wǎng)格運動速度矢量與網(wǎng)格體積變化應該滿足空間守恒律：

2.3 動力學響應計算

基于拉格朗日方程建立MBD模型，車輛模型由各子系統(tǒng)和約束組成。在全局坐標系中，每個約束方程可表示為：

式中：q—廣義坐標系{q}T={x y zψθφ}。歐拉角以z-x-z形式表示；Ck—部件連接處的約束函數(shù)；λk—拉格朗日乘數(shù)；m—作用部件之間節(jié)點的約束數(shù)量。每個部件與相應約束節(jié)點由約束方程聯(lián)結在一起，動力學方程如下：

式中：Φq—約束方程矩陣；Fq—作用力矩陣。

3 數(shù)值計算

3.1 幾何建模

采用一款全尺寸轎車模型，如圖4所示。車輛長L=5.013m，寬W=1.978m，高H=1.458m。CFD 模型與CBD 模型采用相同的坐標系。坐標原點位于重心處（GG）。為提高計算效率，對車身附件進行了簡化。

圖4 車輛三維模型Fig.4 Three-Dimensional Model of the Sedan

選取帶有橋塔的部分橋梁路段作為車輛行駛環(huán)境，如圖5所示。為提高計算效率，對橋梁模型進行簡化，忽略拉索、護欄等部件。橋梁模型長La=170m。橋面的寬度Wa=28m，橋面高度Hb=25m。橋塔寬度Wb=10m，高度Hb=50m。

圖5 橋梁三維模型Fig.5 Three-Dimensional Model of the Bridge

車輛在橋面上的初始位置，如圖6 所示。車輛距離橋塔48m。設vi為車輛行駛速度，全程行駛3.5s。側風垂直于橋梁軸線，設vr為側風速度。

圖6 車輛行駛過程Fig.6 The Process of Driving on the Cross-Sea Bridge Deck

3.2 計算域與邊界條件

對橋梁距離速度入口、壓力出口的凈寬、橋梁上方的凈高、橋梁離地間隙進行不同設置以考察計算域尺度對計算結果的影響，如圖7所示。不同的計算域方案，如表1所示。因橋梁上各橋塔相距較遠，彼此之間氣流干擾較小，為減少計算量，選取一段橋塔，計算域的長度均為170m，設橋寬為特征長度B。

圖7 計算域Fig.7 Computational Domain

側風入口為速度入口，出口邊界為壓力出口，側面邊界和頂部邊界為對稱邊界，橋梁的表面、地面、車輛模型為壁面邊界。車輛以30m/s速度行駛，風速為10m/s，不同尺度計算域計算得到的側向力，如表1所示。當入口和出口的凈寬等于1.0B，車輛上方的凈高等于2.0B，橋下的凈空等于1.0B時，計算結果對計算域大小的依賴性消失。因此最終選取的計算域長度為170m，寬度為84m，高度為88m。

表1 不同計算域的車輛側向力Tab.1 Aerodynamic Force with Different Computational Domain

流場計算中，車輛運動姿態(tài)的變化通過重疊網(wǎng)格技術實現(xiàn)。在重疊網(wǎng)格技術中，網(wǎng)格分為兩部分，背景網(wǎng)格區(qū)域（靜止固定的環(huán)境域）與重疊網(wǎng)格區(qū)域（包含車輛車體的運動域），如圖8所示。

圖8 車輛周圍網(wǎng)格橫向截面圖Fig.8 Mesh Around Sedan

采用STAR-CCM在計算域內生成六面體網(wǎng)格。為保證計算精度，對車輛和橋梁進行網(wǎng)格加密。車身表面和重疊網(wǎng)格區(qū)域的面網(wǎng)格尺寸范圍為（20～180）mm。為模擬附面層流動，在車身表面拉伸出4層棱柱網(wǎng)格，如圖9所示。

圖9 車身表面網(wǎng)格Fig.9 Surface Mesh of Sedan

重疊網(wǎng)格區(qū)域與背景網(wǎng)格區(qū)域之間交界面的網(wǎng)格尺寸大小設置為相同的數(shù)量級，以確保瞬態(tài)數(shù)值模擬計算的穩(wěn)定性。瞬態(tài)計算時間時間步設置為0.001s。

3.3 多體系統(tǒng)動力學模型

汽車是一個多自由度的非線性連續(xù)體振動系統(tǒng)。采用的MBD模型由10個子系統(tǒng)構成，包括車身剛體，前后懸架，前后輪胎，制動，動力和轉向。前懸架采用麥弗遜獨立懸架，轉向采用液壓助力齒輪系轉向系統(tǒng)。應用Magic Formula輪胎模型模擬輪胎動力學模型。整車動力學模型共計85個自由度。根據(jù)實際車型的關鍵參數(shù)調整模型的前后輪距、軸距、四輪定位參數(shù)等。整車的主參數(shù)，如表2所示。整車虛擬樣機模型，如圖10所示。計算仿真中，汽車方向盤鎖死。CFD計算得到的側向力作用于MBD模型重心處。輪胎采用基于Pacejka89模型，參數(shù)為185/65R15。MBD計算的時間步長為0.001s。

表2 車輛整車參數(shù)Tab.2 Vehicle Parameters for Multi-Body Dynamic Mode

圖10 汽車多體系統(tǒng)動力學模型Fig.10 Vehicle Multi-Body Dynamic Model

4 模型驗證

4.1 CFD仿真驗證

在某大學HD-2風洞實驗室進行1：3的車輛模型風洞試驗，驗證CFD 仿真的準確性，如圖11所示。風洞風速為30m/s。測量車輛不同橫擺角下的氣動力和橫擺力矩。橫擺角在-15°至+15°之間變化，每隔3°記錄一次數(shù)據(jù)。

圖11 1：3的轎車模型風洞試驗型Fig.11 A 1:3 Scale Model Wind Tunnel Tests

實驗和仿真數(shù)據(jù)對比表明，如圖12所示。數(shù)值模擬與實驗結果具有良好的一致性。氣動阻力系數(shù)Cd隨橫擺角增大而增大，當橫擺角超過9°過時，阻力系數(shù)Cd增量逐漸減小，氣動側向力系數(shù)隨橫擺角呈線性變化。

圖12 風洞實驗結果與仿真計算結果的比較Fig.12 The Comparisons Between Experiment and Computational Results

4.2 MBD模型魯棒性驗證

作用在車身上的側向力會使車輛發(fā)生偏航和側滑。在研究中，首先在車輛的重心（CG）處施加階躍型側向力（SLF）。然后改變SLF的作用點，計算車輛相應的動力響應。如果動態(tài)響應符合預期，則MBD模型具有良好的魯棒性。

側向力作用點位于重心的后方時，可以有效提高車輛側風穩(wěn)定性。魯棒性驗證中，SLF 的作用點位置如下（:1）向前移動150mm 和300mm（分別用SLFfwd150 和SLFfwd300 表示）；（2）向后移動150mm 和300mm（分別用SLFbwd150 和SLFbwd300 表示）；（3）向上移動200mm和400mm（分別用SLFuwd200和SLFuwd400 表示）；向下移動200mm 和400mm（分別用SLFdwd200 和SLFdwd400表示）。側向力作用在不同位置的側向位移，如圖13（a）所示。當側向力施力點向重心后移動時，車輛的側向位移減小，這與理論結果一致。這表明MBD模型能較好地捕捉XOY平面的運動趨勢。側向力作用在不同作用點下的側傾角，如圖13（b）所示。側向力作用點上移產生負的側傾力矩，風壓中心越高，側傾力矩越大，側向力與側傾力矩共同作用下的側傾角越大，將增加車輛的側傾運動。這表明MBD模型能較好地捕捉YOZ平面的運動趨勢。因此，該車輛多體系統(tǒng)動力學模型具有較強的魯棒性。

圖13 風壓中心移動對車輛運動的影響Fig.13 Effect of Center of Wind Pressure Movement on Vehicle Movement

5 結果分析

采用靜態(tài)和動態(tài)耦合法計算得到的側風環(huán)境下汽車通過橋塔區(qū)域的氣動力變化，如圖14所示。結果表明，靜態(tài)和動態(tài)耦合得到的側向力與橫擺力矩的變化規(guī)律基本一致。在t=1.75s時，車輛駛入橋塔區(qū)域，側風被橋塔遮擋，車輛受到的側向力與橫擺力矩大幅降低。在t=2s時，車輛駛離橋塔區(qū)域，橋塔尾流對側向力和橫擺力矩的影響減小。2s以后，車輛只受到側風作用，氣動力變化趨于穩(wěn)定。動態(tài)耦合法獲得的側向力和橫擺力矩較靜態(tài)耦合法更大。原因在于，車輛此前受到側風與橋塔尾流作用，氣動力變化引起車輛橫向位移、橫擺角速度、車身姿態(tài)等動態(tài)響應。進而又使氣動力的大小和方向不斷變化。因車輛姿態(tài)發(fā)生變化，動態(tài)耦合法較靜態(tài)耦合法，車輛受到的側向力與橫擺力矩更大。

圖14 側向力與橫擺力矩隨時間變化曲線Fig.14 The Lateral Force and Yaw Moment

采用靜態(tài)耦合和動態(tài)耦合法計算得到的側風環(huán)境下車輛通過橋塔區(qū)域的動力學響應，如圖15所示。

圖15 橫擺角與側向位移隨時間變化曲線Fig.15 The Yaw Angle and Lateral Displacement

結果表明，車輛在橫擺力矩的作用下產生了橫擺運動。橫擺角持續(xù)增大直至側風消失。在1.75s 時，車輛駛入橋塔區(qū)域，側風被橋塔遮擋，車輛受到的側向力與橫擺力矩大幅減小。受側風影響，車身姿態(tài)持續(xù)發(fā)生變化。動態(tài)耦合法因考慮車身姿態(tài)與氣動力的相互作用，車輛受到較大的側向力與橫擺力矩。因此，動態(tài)耦合法較靜態(tài)耦合法，車輛橫擺角度更大。側向位移是評價車輛側風穩(wěn)定性的重要參數(shù)。隨著橫擺角的增大，汽車前進速度在側向的分速度也逐漸增大。此外，車輛在側向力的作用下產生側向速度，車輛側向位移持續(xù)增大。動態(tài)耦合法較靜態(tài)耦合法，車輛側向位移更大。并且隨著時間累積，兩者的差距持續(xù)增大。在t=3.5s 時，動態(tài)耦合法下的側向位移達到1.73m，嚴重危及車輛和駕駛員的安全。動態(tài)耦合法中，計算域水平截面不同時刻的壓力云圖以及汽車迎風側與背風側的壓力云圖，如圖16、圖17所示。

圖16 動態(tài)耦合法各個時刻x-y方向壓力云圖（z=0.8）Fig.16 Pressure Distribution at the x-y Section（z=0.8m）

圖17 動態(tài)耦合法各個時刻車輛表面壓力云圖Fig.17 The Surface Pressure Distribution on the Sedan

受側風作用，車輛在流場的姿態(tài)是不斷變化的。從圖16、圖17中可以看出，在t=0.6s，側風速度為0m/s時，流場壓力云圖幾乎是對稱的，而車身表面迎風側和背風側的壓力分布幾乎相同。此時車輛的側向位移為0m。在t=1.2s和t=1.6s時，側風速度為10m/s時，汽車在側風作用下，受到相對穩(wěn)定的側向力和橫擺力矩作用，具有較大的側向速度和橫擺角速度，車輛表面壓力分布不對稱，并且氣流在迎風側C柱和背風側的A柱上發(fā)生分離。迎風側和背風側的壓力差導致側向力和橫擺力矩增加。車輛開始沿側風方向滑行，橫擺角與側向位移逐漸增大。在t＝1.6s時，當車輛接近橋塔時，車輛迎風側的正壓區(qū)減小。車輛背風側A柱附近的氣流分離區(qū)也減小。側向力和橫擺力矩逐漸減小。在t=1.75s時，側風被橋塔遮擋，車輛表面壓力分布幾乎是對稱的，此時反向橫擺力矩使橫擺角減小。

在t=1.9s時，車輛開始離開橋塔尾流區(qū)，車輛背風側的A柱和C柱處發(fā)生氣流分離。車輛迎風側壓力增加，側向力和橫擺力矩再次增加。在t=2.3s時，當車輛完全駛離橋塔并再次駛入側風區(qū)域時，車輛表面壓力分布與進入橋塔區(qū)類似。不同之處在于，由于車輛的橫擺角變大，車身姿態(tài)發(fā)生變化，迎風側C柱的氣流分離面積減小。在t=3.5s時側風消失，車輛的表明壓力分布由不對稱逐漸趨于對稱。整個過程，車身姿態(tài)的變化影響了流場及氣動力。

6 結論

采用動靜態(tài)耦合法對側風環(huán)境下車輛通過橋塔區(qū)域時的側風穩(wěn)定性進行分析，得出結論如下：

（1）側風環(huán)境中，橋塔周圍會產生繞流和脫落渦。車輛通過該區(qū)域時，側向力和橫擺力矩的數(shù)值先減少后增加，車輛受橋塔遮蔽時出現(xiàn)反方向側向力。橫擺力矩減小至36N·m，惡劣的側風環(huán)境使車輛的行駛穩(wěn)定性變差。

（2）側風環(huán)境下車輛通過橋塔區(qū)域，動態(tài)耦合法考慮了車輛受側風作用而產生的運動狀態(tài)的變化，獲得動力學響應更符合實際情況，計算得到的動力學響應如側向位移為1.75m，橫擺角為2.12°，相比于靜態(tài)耦合法分別增加了12.91%、17.12%。

（3）車輛受側風影響，橫擺角持續(xù)增加，車身姿態(tài)的變化影響外部流場，進一步影響氣動力。因此，研究車輛側風穩(wěn)定性應考慮氣動力和車輛動態(tài)響應之間的相互作用。