極限運行工況下液罐雙向流固耦合分析

蔣德飛，張竹林，鄒彥冉

山東交通學院 汽車工程學院，濟南 250357

0 引言

液罐車輛在道路運輸過程中由于行駛狀態(tài)的改變，罐體受到慣性力的作用引起液體晃動，表現(xiàn)出罐體與充裝介質復雜的耦合效應[1]。液體的晃動導致罐車行駛不平穩(wěn)，影響罐車的操縱穩(wěn)定性及行車安全，同時晃動的液體產生瞬態(tài)載荷易導致罐體應力集中，造成罐體及防波板損壞。利用雙向流固耦合方法對罐體內液體的晃動進行仿真分析，研究罐體的結構受力對保證罐車的安全穩(wěn)定運行具有理論價值和實際意義。

近年來，國內外相關學者針對罐體內液體晃動對罐體的作用力進行了大量研究。何華等[2]采用Fluent軟件仿真證明增加防波板數(shù)量可減小液體晃動，防波板數(shù)量增多時液體晃動的峰值載荷明顯降低，同時對罐車不同制動初速度下罐體的瞬態(tài)應力進行對比分析，證明峰值應力存在時間與制動初速度有關；陳志偉[3]以充裝液體介質的罐體為基礎，建立液體晃動模型，分析剛性壁面應力分布和罐體內充裝系數(shù)的大小及其對防波板和封頭作用力的影響；劉忠亮[4]提出可通過對罐體V型支座的結構和位置的優(yōu)化設計，使罐體封頭和筒體中的應力最??；張韜等[5]對罐體制動工況液體的晃動進行流固耦合分析，得出罐車在達到極限制動工況0.52 s后液體晃動對罐體的壓力達到最大；蔡昌全等[6]建立液化氣體運輸半掛車的整體模型，進行應力分析，對應力集中位置進行應力分類并進行強度校核，以確保滿足運輸安全要求；王若平等[7]按照文獻[8]中的相關規(guī)定分別對縱向工況、橫向工況、垂直向上工況、垂直向下工況進行應力、應變分析，指出罐體結構設計中的薄弱環(huán)節(jié)，為鋁合金罐車的設計制造提供了可靠依據。

上述研究大多集中于采用單向流固耦合方式將仿真分析過程中的數(shù)據從流體傳遞給固體，該方法不易得到耦合過程中罐體的應力變化及分布規(guī)律，同時也忽略了罐體的結構變形對流體的影響，仿真結果相比雙向流固耦合法求解精度較差。雙向流固耦合方法流體與結構進行雙向數(shù)據交換，固體以邊界位移為載荷傳遞給流體，流體域將計算的流場數(shù)據返回到固體邊界上，該方法能夠得到結構域在流場提供的邊界條件下的結構形變、應力分布等。

相對于其他行駛工況，罐車在制動和轉向時罐體受力較復雜。根據文獻[9]的要求，本文基于罐車在附著系數(shù)最大的瀝青混凝土路面上制動時沿行駛方向不發(fā)生抱死拖滑可達到的最大制動減速度，以及轉向時罐車沿行駛方向不發(fā)生拖滑的最大減速度和不發(fā)生側翻的最大側向減速度，采用雙向流固耦合方法系統(tǒng)分析罐體在制動和轉向2種極限工況下，不同充液比下的罐體所受應力隨時間的變化趨勢以及防波板的應力分布情況，對罐體的強度進行檢算。

1 計算模型與計算方法

1.1 罐體模型和網格劃分

某罐車罐體總長度為10 350 mm，寬度為2410 mm，罐體內部包含有6個防波板，將罐體內腔室分隔為7個艙室，相鄰2個防波板上的過人孔交錯布置，罐體結構如圖1所示(圖中單位mm)。

本文側重于研究極限運行工況下罐內液體晃動產生的動載荷在液罐車罐體上的應力分布，采用三維軟件SolidWorks構建罐體物理模型，將物理模型導入ANSYS-Workbench中，建立罐體有限元模型，采用ANSYS-Workbench中基于直接建模思想的重要幾何前處理軟件ANSYS-SpaceClaim對罐體在仿真過程中可能存在的尖角、小邊及縫隙等進行檢查和修補[10]，得到較精確的罐體有限元模型如圖2所示。

a)罐體結構簡圖 b)防波板結構示意圖 圖1 罐體結構尺寸 圖2 罐體有限元模型

根據圖2的罐體有限元模型，通過ANSYS-SpaceClaim填充得到罐體的結構域與流體域模型[11]。通過ANSYS-Mechanical模塊(該模塊是ANSYS重要的前處理軟件之一，主要功能包括網格智能生成和各種結果的數(shù)據處理)對罐體結構域和流體域模型進行網格劃分，需要保證結構域和流體域網格尺寸與形狀的一致性，網格過渡部分保證平滑且節(jié)點相互對應，可減少網格總數(shù)，提高耦合數(shù)據交互的便利性，避免在軟件分析中產生模型的網格不匹配，影響仿真結果，提高計算效率。網格劃分后，罐體結構域有限元模型包括277 264個單元，553 300個節(jié)點，罐體流體域模型包括352 429個單元，66 644個節(jié)點，結構域和流體域有限元模型如圖3、4所示。

圖3 罐體結構域網格模型 圖4 罐體流體域網格模型

1.2 雙向流固耦合

流體域模型可用于模擬罐體內部一定量的空氣和水，根據相關約束和工況求解罐體結構域與流體域網格模型，流體域與結構域的數(shù)據相互交換，且在各個時間步都相互作用。采用雙向流固耦合方法能夠得到結構域在流場提供的邊界條件下的時間歷程，如結構形變、應力等變量的時間歷程[12]。

雙向流固耦合方法仿真精度相比單向流固耦合更符合實際運行工況，但對罐體的結構域和流體域模型進行仿真分析計算量較大，雙向耦合分析具體流程如圖5所示。進行瞬態(tài)計算時，罐體的結構域模型在ANSYS-Transient Structural中進行計算得到罐體的位移變形量；罐體的流體域模型在ANSYS-Fluent中仿真計算得到罐體在液體晃動下的壓力、速度分布；結構域和流體域通過ANSYS-System Coupling完成耦合數(shù)據的交互傳遞。在耦合分析過程中，將總時間劃分為n個時間步，在每個時間步內都要進行流體域與固體域的求解和數(shù)據傳遞，同時該時間步的結果作為下一時間步的初值進行計算，繼而完成所有時間步求解[13-16]。

圖5 雙向流固耦合分析流程

1.3 兩相流計算模型

本文對于流體域中流體運動采用ANSYS-Fluent中的流體體積模型(volume of fluid model,VOF)模擬，用來處理無相互穿插的多相流問題。

采用體積率函數(shù)表示流體自由面的位置和流體所占的體積，用每個單元內含有流體的體積和單元體積的比β表示流體的填充程度，當β=0時，該單元不含有流體；當β=1時，該單元充滿流體；當0<β<1時，該單元含有自由液面。

VOF模型通過求解運動方程組和處理穿過區(qū)域的流體可以精確模擬氣液分界面的變化情況和液體晃動過程[17-18]。多相流VOF模型的基本控制方程包括流體連續(xù)性方程與動量守恒方程。

1.3.1 流體連續(xù)性方程

式中：u、v、w為流體質點在x、y、z軸方向上的速度分量，ρ為流體質點處的密度，t為時間。

1.3.2 動量守恒方程

式中：Fx、Fy、Fz分別為流體質點在x、y、z方向上力的分量，p為流體質點處的壓力。

1.4 求解器參數(shù)設置

罐體的流體域模型在ANSYS-Fluent中采用壓力基(pressure-based)進行求解，鑒于罐體中氣液兩相的交界面隨著運算時間發(fā)生變化，因此采用瞬態(tài)求解。在ANSYS-Fluent中湍流模型采用標準k-epsilon數(shù)值模型[19]。

流體域求解采用基于非迭代算法(pressure-implicit with splitting of operators，PISO)，該算法以壓力基為理論基礎，為瞬態(tài)不可壓縮流動設計的一種求解方法，這種算法在求解瞬態(tài)問題時有明顯優(yōu)勢，瞬態(tài)計算時，PISO算法能夠在大時間步、動量及壓力松弛因子(取值范圍為0～1，通過改變壓力松弛因子可控制收斂的速度和改善收斂情況)為1時仍能保持計算穩(wěn)定，提高收斂速度[20-21]。

結構域與流體域有限元分析模型在仿真分析時需要定義材料屬性。結構域材料采用鋁合金，該材料屈服極限為325 MPa；為了安全起見，在對罐體進行應力校驗時用屈服應力除以安全系數(shù)得到許用應力，根據文獻[22-23]的規(guī)定，設計安全系數(shù)取1.5，計算得到罐體材料的許用應力為217 MPa。流體域中的氣液兩相模型通過ANSYS-Fluent中的Patch功能將流體域劃分為氣、液2部分，本文首先研究罐體充液比(罐體充裝液體的體積與罐體容積的比值)為50%時的液體晃動，即50%的氣體沿y軸正向分布，50%的液體沿y軸負向分布，氣體和液體材料分別設置為空氣和水；通過ANSYS-Fluent中Cell Registers建立不同充液比下的流體模型，從而實現(xiàn)多種充液比下罐體受力的仿真分析。

為防止流體域發(fā)生液面晃動后可能產生的網格變形、歪曲率過大或產生負體積網格[24]，本文動網格采取彈簧光順網格控制方法，其計算控制方程為：

式中：Fi為網格節(jié)點i所受的力，Δxi、Δxj分別為網格節(jié)點i、j的位移，ni為與i相連的網格節(jié)點數(shù)，kij為i與j之間的彈簧剛度。

1.5 邊界條件

雙向流固耦合分析過程中罐體結構域內表面和流體域外表面需要進行耦合數(shù)據交互，在ANSYS-Transient Structural中可完成流體域耦合面的設置；結構域內表面命名為Solid_Move，流體域外表面命名為Fluid_Move，方便在ANSYS-System Coupling模塊中創(chuàng)建數(shù)據交互通道，完成仿真過程中的雙向數(shù)據交換。

設定罐體z方向和x方向對應的減速度用來模擬罐體制動工況和轉向工況；流體域在Fluent中的邊界條件設置為：給定流體域y軸負方向的重力加速度g=9.81 m/s2，不考慮傳熱影響，設置壓力參考點，定義耦合面為流體區(qū)域的外表面，打開變形的相鄰邊界層，建立初始液體區(qū)域和氣體區(qū)域，在液面處設置101 kPa，進行全局初始化。為了提高罐體雙向流固耦合仿真精度，求解時間步長設置為0.005s，罐體結構域計算求解子步設置為3，總的模擬仿真時間為2 s，分析達到極限工況后2 s內流體域的運動參數(shù)和結構域的應力。

罐車行駛在不同路面時，在瀝青或干混凝土路面上的附著系數(shù)最大，其滑動附著系數(shù)為0.75，制動減速度達到7.34 m/s2；峰值附著系數(shù)為0.8，制動減速度最高可達到7.8 m/s2，因此在模擬汽車制動工況時最大制動減速度為8.0 m/s2。

由于在良好路面上，車輛輪胎的附著系數(shù)最高可達0.8，即側向減速度達到0.8g時，汽車開始側滑；而罐車的側翻閾值為0.4～0.6g，當側向減速度達到0.4g，罐車尚未發(fā)生側滑時，就已發(fā)生側翻[25]，因此在模擬轉向工況時，罐體z軸方向的最大減速度為8.0 m/s2，x軸方向的最大減速度為4.0 m/s2。

2 極限工況下罐體應力計算結果及分析

2.1 極限制動工況

在極限制動工況下對罐體進行瞬態(tài)分析，仿真分析完成后可在ANSYS-Fluent后處理軟件CFD-Post中查看不同時刻的罐內氣液2相分布圖，觀察仿真過程中的流體運動狀況。氣液2相圖選取的平面為流體域沿y軸方向的縱切面；罐體內的液體晃動如圖6所示，罐體后端的液體不斷向罐體前端流動，同時空氣通過防波板向罐體后端流動，最終罐體前端液體增多，空氣與液體的交界面形成傾斜的近似斜面。

圖6 極限制動工況下不同時刻罐體內的氣液2相分布

罐體與防波板的應力分布分別如圖7、8所示。由圖7可知：罐體在晃動過程中最大應力為55.709 MPa，出現(xiàn)在罐體前封頭與筒體接觸的部位。由圖8可以看出：防波板每處過人孔所受的應力沿汽車行駛方向逐漸增大，最大應力為41.596 MPa，出現(xiàn)在罐體沿行駛方向最前方的防波板直徑最大的通液口處。

圖7 極限制動工況下罐體應力云圖 圖8 極限制動工況下防波板應力云圖

2.2 極限轉向工況

在極限轉向工況下對罐體進行瞬態(tài)分析，不同時刻罐體內氣液2相分布情況如圖9所示，以0.4 s為時間間隔分別列出罐車達極限轉向工況后0、0.4、0.8、1.2、1.6、2.0 s時罐體的氣液2相圖，轉向時罐體具有側向加速度，氣液2相圖采用立體模型顯示。從圖9中可以看出：極限轉向工況初始時刻液體域平面呈水平狀，隨著時間延長，罐體后端艙室的液體涌入前端艙室的同時，也有沿x方向的側向運動。

圖9 極限轉向工況下不同時刻罐體內的氣液2相分布

與極限制動工況相比，封頭、防波板和筒體的峰值應力都明顯增大，罐體最大應力為126.620 MPa，如圖10所示，約為制動工況的2倍，防波板最大應力為89.073 MPa，如圖11所示。與制動工況相比，罐體與防波板的承載應力區(qū)域位置發(fā)生了變化，罐體最大應力集中在前封頭與筒體接觸的上側區(qū)域，防波板的最大應力發(fā)生在沿側向加速度方向的防波板上部區(qū)域。

圖10 極限轉向工況下罐體應力云圖 圖11 極限轉向工況下防波板應力云圖

2.3 不同充液比下的仿真數(shù)值模擬

2.3.1 極限制動工況

為了分析同一運行工況下罐體充液比不同時罐車罐體的受力情況，選取充液比分別為30%、50%、80%，研究罐體及防波板的受力隨時間的變化情況。極限制動工況下不同充液比時罐體及防波板的最大應力隨時間的變化曲線如圖12、13所示。

圖12 極限制動工況下罐體最大應力時間歷程 圖13 極限制動工況下防波板最大應力時間歷程

由圖12、13可知：1)各防波板由于受到液體載荷沖擊，應力均隨時間呈波動性變化，其應力時間歷程與罐體相似。2)隨著充液比的增大，罐體和防波板所受應力逐漸增大，且應力增長速度加快，當充液比為80%時，罐體和防波板所受的應力最大，分別為94.433、70.846 MPa；隨著充液比的增大，最大應力出現(xiàn)的時間延后，峰值應力出現(xiàn)在極限制動工況后0.07～0.12 s；隨著時間繼續(xù)增長，應力幅值逐漸減小。

2.3.2 極限轉向工況

極限轉向工況下，充液比為30%、50%、80%時罐體和防波板的最大應力變化趨勢與制動工況相似，如圖14、15所示。由圖14、15可以看出：當充液比為80%時，罐體與防波板所受的應力最大，分別為223.12、148.76 MPa,出現(xiàn)在極限轉向工況后0.05～0.10 s。

圖14 極限轉向工況下罐體最大應力時間歷程 圖15 極限轉向工況下防波板最大應力時間歷程

由極限工況下罐體應力計算結果分析可知：罐體與防波板的最大應力均發(fā)生在充液比為80%的極限轉向工況，分別為223.12、148.76 MPa，分別位于罐體筒體與前封頭的連接處與防波板沿側向加速度方向的上部區(qū)域。

由于存在邊緣應力的影響，罐體局部應力較大，屬于一次應力中的局部薄膜應力，根據文獻[22]計算局部薄膜應力強度的許用極限應力為325 MPa，因此該工況下罐體與防波板均滿足強度要求。

3 結論

1)本文以某型號液體運輸車輛罐體為研究對象，建立罐體結構域與流體域模型，采用雙向流固耦合方法利用有限元軟件ANSYS對罐車在制動和轉向2種極限運行工況下罐體內液體的沖擊晃動進行仿真分析，得到充液比為50%時2種極限工況下罐體和防波板的應力時間歷程和氣液2相分布圖，從而得到2種工況下罐體與防波板的應力。

2)分析制動、轉向2種極限運行工況充液比分別為30%、50%、80%時罐體與防波板的受力情況，得到峰值應力分別出現(xiàn)在極限制動工況后0.07～0.12 s、極限轉向工況后0.05～0.10 s。罐體與防波板的最大應力均發(fā)生在充液比為80%的極限轉向工況，分別為223.12、148.7 MPa。罐體的最大應力部位位于罐體筒體與前封頭的連接處，防波板最大應力集中在防波板沿側向加速度方向的上部區(qū)域，2處應力均小于許用極限應力325 MPa，故罐體與防波板均滿足強度要求。