蔣德飛，張竹林，鄒彥冉

山東交通學(xué)院 汽車工程學(xué)院，濟(jì)南 250357

0 引言

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

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

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

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

1 計(jì)算模型與計(jì)算方法

1.1 罐體模型和網(wǎng)格劃分

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

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

a)罐體結(jié)構(gòu)簡(jiǎn)圖 b)防波板結(jié)構(gòu)示意圖 圖1 罐體結(jié)構(gòu)尺寸 圖2 罐體有限元模型

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

圖3 罐體結(jié)構(gòu)域網(wǎng)格模型 圖4 罐體流體域網(wǎng)格模型

1.2 雙向流固耦合

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

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

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

1.3 兩相流計(jì)算模型

本文對(duì)于流體域中流體運(yùn)動(dòng)采用ANSYS-Fluent中的流體體積模型(volume of fluid model,VOF)模擬，用來(lái)處理無(wú)相互穿插的多相流問(wèn)題。

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

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

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

式中：u、v、w為流體質(zhì)點(diǎn)在x、y、z軸方向上的速度分量，ρ為流體質(zhì)點(diǎn)處的密度，t為時(shí)間。

1.3.2 動(dòng)量守恒方程

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

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

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

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

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

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

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

1.5 邊界條件

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

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

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

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

2 極限工況下罐體應(yīng)力計(jì)算結(jié)果及分析

2.1 極限制動(dòng)工況

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

圖6 極限制動(dòng)工況下不同時(shí)刻罐體內(nèi)的氣液2相分布

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

圖7 極限制動(dòng)工況下罐體應(yīng)力云圖 圖8 極限制動(dòng)工況下防波板應(yīng)力云圖

2.2 極限轉(zhuǎn)向工況

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

圖9 極限轉(zhuǎn)向工況下不同時(shí)刻罐體內(nèi)的氣液2相分布

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

圖10 極限轉(zhuǎn)向工況下罐體應(yīng)力云圖 圖11 極限轉(zhuǎn)向工況下防波板應(yīng)力云圖

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

2.3.1 極限制動(dòng)工況

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

圖12 極限制動(dòng)工況下罐體最大應(yīng)力時(shí)間歷程 圖13 極限制動(dòng)工況下防波板最大應(yīng)力時(shí)間歷程

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

2.3.2 極限轉(zhuǎn)向工況

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

圖14 極限轉(zhuǎn)向工況下罐體最大應(yīng)力時(shí)間歷程 圖15 極限轉(zhuǎn)向工況下防波板最大應(yīng)力時(shí)間歷程

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

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

3 結(jié)論

1)本文以某型號(hào)液體運(yùn)輸車輛罐體為研究對(duì)象，建立罐體結(jié)構(gòu)域與流體域模型，采用雙向流固耦合方法利用有限元軟件ANSYS對(duì)罐車在制動(dòng)和轉(zhuǎn)向2種極限運(yùn)行工況下罐體內(nèi)液體的沖擊晃動(dòng)進(jìn)行仿真分析，得到充液比為50%時(shí)2種極限工況下罐體和防波板的應(yīng)力時(shí)間歷程和氣液2相分布圖，從而得到2種工況下罐體與防波板的應(yīng)力。

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