石懷龍 王勇 鄔平波

(西南交通大學(xué) 牽引動(dòng)力國家重點(diǎn)實(shí)驗(yàn)室,成都 610031)

引言

液體晃動(dòng)會(huì)影響貨物運(yùn)輸安全,如公路、鐵路和航空運(yùn)輸中的液罐車以及海底運(yùn)輸管道等,主要是由于液體自由表面發(fā)生非線性運(yùn)動(dòng),對(duì)罐體施加動(dòng)態(tài)載荷.液體晃動(dòng)行為研究起源于航空航天領(lǐng)域,主要難點(diǎn)為容器和液體接觸面之間的運(yùn)動(dòng)邊界條件模擬、液體自由表面的非線性運(yùn)動(dòng)模擬等.

Abramson等人基于線性彈簧振子模型研究了航天器內(nèi)儲(chǔ)油箱的液體振動(dòng),分析具有不同罐體形狀內(nèi)的液體晃動(dòng)頻率和振型,并給出接觸作用力解析表達(dá)式,指出液體第一階晃動(dòng)頻率對(duì)系統(tǒng)影響最大[1].Ranganathan等人建立了單擺-固定質(zhì)量塊模型,用單擺模擬液體一階晃動(dòng)效應(yīng),用固定質(zhì)量塊模擬未參與晃動(dòng)的部分液體,并應(yīng)用于公路槽罐車制動(dòng)過程中的液體晃動(dòng)行為[2,3].文獻(xiàn)[4,5]詳細(xì)闡述了等效模型研究進(jìn)展,分析了不同建模方法和充比液時(shí)重載槽罐車的傾覆性能.然而,這些典型的等效物理模型過于簡化,不能描述液體慣量和自由表面的連續(xù)性變化.Celebi基于Navier-Stokes方程建立了連續(xù)非線性模型,研究帶有隔板的矩形貯箱內(nèi)的液體晃動(dòng)[6].Ibrahim[7,8]總結(jié)了二十一世紀(jì)初之前的液體晃動(dòng)研究進(jìn)展,明確了液體晃動(dòng)的被動(dòng)控制方法以及其在結(jié)構(gòu)振動(dòng)抑制中的應(yīng)用,討論了液體的非線性屬性和結(jié)構(gòu)非線性振動(dòng)等編程實(shí)現(xiàn)的困難.Aliabadi[9]對(duì)比了連續(xù)流體模型與彈簧振子和單擺系統(tǒng)構(gòu)成的等效模型的結(jié)果差異,分析了液體自由液面的變化和動(dòng)態(tài)慣量特性,這些都是等效模型所不能實(shí)現(xiàn)的.文獻(xiàn)[10]總結(jié)了近期液固耦合的建模手段和數(shù)值方法,并與試驗(yàn)數(shù)據(jù)對(duì)比分析各自優(yōu)缺點(diǎn).Bogomaz[11]和Vera等人[12]采用單擺等效模型模擬液體,建立考慮液體晃動(dòng)的三維車輛系統(tǒng)動(dòng)力學(xué)模型,分析車輛在啟動(dòng)、制動(dòng)和連掛以及過曲線時(shí)的車鉤振動(dòng)、車體加速度變化情況,指出液體晃動(dòng)對(duì)車輛承載和振動(dòng)響應(yīng)產(chǎn)生顯著影響. Wang等[13]和Wei等[14]較早地采用連續(xù)介質(zhì)力學(xué)方法模擬矩形貯箱內(nèi)液體的晃動(dòng)效應(yīng),可模擬液體表面的連續(xù)性變化.

虛擬質(zhì)量法將液體考慮為附加質(zhì)量代入系統(tǒng)質(zhì)量陣,借助有限元軟件計(jì)算當(dāng)前液體質(zhì)量,然后通過插值得到任意時(shí)刻的液體質(zhì)量,即做了線性疊加假設(shè)[15].劉富等[16]通過光滑粒子流體動(dòng)力方法研究矩形容器內(nèi)隔板設(shè)置問題,合理設(shè)置隔板可有效抑制液體晃動(dòng).等效力學(xué)模型需借助流體動(dòng)力學(xué)和有限元仿真等手段推算或擬合模型參數(shù),并基于線性假設(shè)組建車輛系統(tǒng)三維晃動(dòng)模型,可提高計(jì)算效率,但僅適用于基本規(guī)律性研究,不能描述液體自由表面和慣量的非線性變化[17,18].

綜上,有必要提出可適用于描述液體慣量連續(xù)變化的非線性模型,并易于與車輛多體系統(tǒng)動(dòng)力學(xué)模型集成.Shabana提出絕對(duì)節(jié)點(diǎn)坐標(biāo)方法(Absolute Nodal Coordinate Formulation,ANCF),在全局坐標(biāo)系下用節(jié)點(diǎn)位失的梯度代替?zhèn)鹘y(tǒng)有限元中的轉(zhuǎn)角坐標(biāo),形成的柔性體質(zhì)量矩陣為常數(shù)陣,在大變形問題研究中具有高精度等特點(diǎn)[19,20].

本文基于采用絕對(duì)節(jié)點(diǎn)坐標(biāo)方法的柔性多體系統(tǒng)動(dòng)力學(xué)理論,提出一種液體晃動(dòng)模擬新方法,并應(yīng)用于液罐車內(nèi)的液體晃動(dòng)問題.首先,基于流體力學(xué)基礎(chǔ)理論推導(dǎo)液體粘性方程和體積不可壓條件方程,然后,采用絕對(duì)節(jié)點(diǎn)坐標(biāo)方法描述的實(shí)體單元對(duì)液體進(jìn)行網(wǎng)格劃分,并通過罰函數(shù)描述液體與罐體的接觸作用,最后組建罐體-液體耦合多體系統(tǒng)動(dòng)力學(xué)模型,仿真研究液罐車內(nèi)液體的橫向和縱向晃動(dòng)行為.

1 基于絕對(duì)節(jié)點(diǎn)坐標(biāo)方法的柔性多體系統(tǒng)

1.1 運(yùn)動(dòng)學(xué)描述

1.2 慣性力表達(dá)

(1)

1.3 彈性力表達(dá)

(2)

其中,單元?jiǎng)偠汝嘖為節(jié)點(diǎn)坐標(biāo)的強(qiáng)非線性陣；格林-拉格朗日應(yīng)變張量ε定義為：

(3)

圖1 單元各構(gòu)型定義及相互關(guān)系Fig.1 ANCF deformable body configurations

1.4 廣義力表達(dá)

假設(shè)單元任意一點(diǎn)處受外力矢量作用f,則外力對(duì)應(yīng)的虛功為fTδr,其中r為單元上的力作用點(diǎn)處的位置矢量,δr為r的虛變量.為獲得與系統(tǒng)廣義坐標(biāo)相關(guān)的廣義外力,將廣義力寫成節(jié)點(diǎn)坐標(biāo)虛變量形式為：

(4)

則Qf=STf對(duì)應(yīng)外力f的廣義力.

1.5 運(yùn)動(dòng)微分代數(shù)方程組

由于絕對(duì)節(jié)點(diǎn)坐標(biāo)方法的位置矢量定義在全局坐標(biāo)系下,因此其運(yùn)動(dòng)約束方程描述和雅克比矩陣形式較簡單.因此,這里僅給出含非線性約束方程的增廣式多體系統(tǒng)運(yùn)動(dòng)微分代數(shù)方程組：

(5)

式中,Qe為系統(tǒng)廣義外力,Qd為與約束方程相關(guān)的廣義約束力,λ為拉格朗日乘子,Cq為與約束方程相關(guān)的雅克比矩陣.

2 流體力學(xué)基礎(chǔ)理論

2.1 Navier-Stokes方程

流體力學(xué)研究方法通常可分為兩類,歐拉法關(guān)注空間網(wǎng)格點(diǎn)處的流體流動(dòng)情況,而拉格朗日法則追蹤流體物質(zhì)點(diǎn)運(yùn)動(dòng)軌跡.當(dāng)流體粘性屬性只與溫度場(chǎng)相關(guān),且內(nèi)部應(yīng)力與應(yīng)變線性相關(guān),稱該類流體滿足牛頓型特征,適用牛頓第二定律描述其運(yùn)動(dòng).假設(shè)流體為不可壓的牛頓類,連續(xù)性方程寫為：

(6)

式中,ρ為流體密度(kg/m3),t為時(shí)間(s),r和v分別為位置(m)和速度矢量(m/s),為散度算子.假設(shè)流體不可壓縮,即流體密度保持不變,則連續(xù)性方程式(6)簡化為·v=0,即不可壓條件.流體運(yùn)動(dòng)偏微分方程為：

(7)

σ={-p+λtr(D)}I+2μD

(8)

將公式(8)代入(7)得流體單元狀態(tài)平衡方程：

(9)

(10)

式中,J=|J|；Cr=JTJ為右柯西-格林應(yīng)變張量.

2.2 流體不可壓條件

質(zhì)量守恒條件為dm=ρodV=ρdv,其中ρo和dV分別為流體在初始構(gòu)型時(shí)密度和體積,而ρ和dv分別為當(dāng)前構(gòu)型即發(fā)生形變后的密度和體積.流體變形前后的體積關(guān)系為[22]：

dv=rxdx·(rydy×rzdz)

=(rx·(ry×rz))dxdydz=JdV

(11)

上式兩邊同時(shí)乘以ρoρ,再利用質(zhì)量守恒關(guān)系可得出流體發(fā)生形變前后的密度關(guān)系為ρo=ρJ.若流體不可壓條件為J=1,表明流體發(fā)生形變前后的體積和密度均保持不變.

3 基于絕對(duì)節(jié)點(diǎn)坐標(biāo)方法的液體運(yùn)動(dòng)描述

3.1 六面體實(shí)體單元

(12)

式中,a、b和c分別為單元沿軸x、y和z方向的長度；ξ=x/a,η=y/b,ξ=z/c,其中,ξ、η和ζ∈[0,1],ξk、ηk和ζk為節(jié)點(diǎn)k的自然坐標(biāo).單元j上任意一點(diǎn)的位置矢量表示為：

(13)

式中,I為3×3單位陣,Sj和ej分別為單元的形函數(shù)矩陣和節(jié)點(diǎn)坐標(biāo)矢量.六面體實(shí)體單元含8個(gè)節(jié)點(diǎn),共計(jì)96個(gè)自由度,其位置矢量及梯度矢量可以表達(dá)液體形狀及其慣量變化.

3.2 體積不可壓條件

(14)

(15)

3.3 單元內(nèi)部粘性力

其中,vj為單元當(dāng)前構(gòu)型的體積,Vj為單元初始構(gòu)型的體積.結(jié)合公式(10),液體單元的內(nèi)部粘性力虛功表達(dá)式可簡化為：

(16)

(17)

3.4 液體運(yùn)動(dòng)方程

根據(jù)液體單元的質(zhì)量矩陣、外力矢量(式(4))、保證體積不可壓條件的約束力(式(14)和(15))和粘性力(式(17)),則其運(yùn)動(dòng)微分方程組為：

(18)

若液體被劃分為多個(gè)單元,通過組裝單元運(yùn)動(dòng)方程來組建系統(tǒng)運(yùn)動(dòng)方程,利用數(shù)值仿真求得系統(tǒng)速度和位置矢量.

4 基于罰函數(shù)描述的液體與罐體接觸關(guān)系

(19)

式中,K和C為罰函數(shù)剛度和阻尼系數(shù).

由于正壓力fn應(yīng)垂直于罐體上接觸點(diǎn)處切平面,因此需要計(jì)算接觸點(diǎn)處的法向量v.已知罐體結(jié)構(gòu)兩端為半橢球體,中間部分為圓柱體,如圖2所示,可根據(jù)液體空間位置計(jì)算出罐體上對(duì)應(yīng)的接觸點(diǎn)位置.摩擦力fc=μfn沿接觸點(diǎn)處切平面方向,μ為摩擦系數(shù).正壓力和摩擦力同時(shí)作用在液體和罐體上,幅值相等但方向相反.

圖2 液體與罐體接觸關(guān)系的坐標(biāo)系定義Fig.2 Interaction definition between the liquid and tank

(20)

(21)

4.1 中部圓柱體區(qū)域接觸關(guān)系

(22)

式中：r為圓柱體半徑,假設(shè)罐體為剛體.接觸點(diǎn)處穿透速度定義為：

(23)

作用在液體上的力為：

(24)

(25)

4.2 端部橢球體區(qū)域接觸關(guān)系

(26)

5 液罐車內(nèi)液體晃動(dòng)行為仿真

以液罐車裝載類似水的液體為例,密度ρ=1000kg/m3,重力系數(shù)g=-9.81m/s2,粘性系數(shù)取μ=9.3×10-4Pa[14].體積不可壓罰函數(shù)系數(shù)kIC=1.0×106N/m和cTD=1.0×104N·s/m.罐體長度12m,直徑2.8m.模擬液罐車半載情況,液面高度1.4m,則液體質(zhì)量約37.5t.

5.1 橫向晃動(dòng)模擬

在罐體上施加橫向簡諧激勵(lì)f=0.5Hz、A=0.1m,仿真時(shí)間4s,模擬液罐車內(nèi)液體的橫向往復(fù)晃動(dòng)過程,如圖3和圖4所示.t=0.0s時(shí),半圓柱形液體在重力作用下開始與罐體內(nèi)壁發(fā)生接觸,最終構(gòu)型與罐體內(nèi)壁貼合；t=1.5s、2.5s時(shí),液體發(fā)生縱向、橫向和垂向變形以貼合罐體內(nèi)壁,但不同斷面處的液體自由表面形狀不完全相同,即液體自由表面的運(yùn)動(dòng)呈現(xiàn)非線性；t=3.0s、4.0s時(shí)可以觀察到罐體內(nèi)液體的顯著晃動(dòng)行為,液體上表面位置呈非線性變化.仿真計(jì)算時(shí)僅采用了1個(gè)單元,參考第1.3節(jié)和第4節(jié)定義單元初始外形,以盡可能保證液體初始形狀與罐體內(nèi)壁貼合,從而節(jié)省計(jì)算工時(shí).通過增加單元數(shù),可以獲得更加連續(xù)的液體自由表面變化形狀,但基本上與1個(gè)單元情況相似,因此不再列寫多個(gè)單元時(shí)的仿真結(jié)果.

5.2 縱向晃動(dòng)模擬

對(duì)罐體施加縱向簡諧激勵(lì)f=0.1Hz、A=1m,模擬液罐車內(nèi)液體的縱向晃動(dòng)過程,仿真時(shí)間20s,結(jié)果如圖5所示.可見,液體的縱向晃動(dòng)必然引起液體質(zhì)心的縱向移動(dòng).仿真計(jì)算表明其質(zhì)心最大位移可達(dá)0.5m以上,按此值估算液罐車前后轉(zhuǎn)向架軸重轉(zhuǎn)移量可達(dá)15%,從而影響液罐車的動(dòng)力學(xué)性能和車輪磨耗.液罐車采用長大編組方式運(yùn)行,制動(dòng)過程中的液體縱向晃動(dòng)必然引起車鉤/緩沖器受力的突變或制動(dòng)效能損失.

圖3 液罐車內(nèi)液體的橫向晃動(dòng)(主視圖)Fig.3 Liquid sloshing along lateral direction (axial view)

6 結(jié)論

1)本文提出采用絕對(duì)節(jié)點(diǎn)坐標(biāo)方法的液體晃動(dòng)模擬新方法,用于模擬液罐車內(nèi)液體自由表面的連續(xù)性變化.該方法描述液體自由表面形狀的連續(xù)性變化,并適用于研究具有復(fù)雜結(jié)構(gòu)或形狀容器的內(nèi)部液體晃動(dòng)問題.

2)基于流體力學(xué)牛頓類體基礎(chǔ)理論,推導(dǎo)粘性方程和滿足體積不可壓縮的條件方程；采用絕對(duì)節(jié)點(diǎn)坐標(biāo)方法描述的實(shí)體單元進(jìn)行液體網(wǎng)格劃分.

3)采用罰函數(shù)方法描述液體與罐體的接觸作用,給出液體與罐體中部圓柱體或端部橢球體的接觸關(guān)系,組建罐體-液體耦合系統(tǒng)動(dòng)力學(xué)模型.

圖4 液罐車內(nèi)液體的橫向晃動(dòng)(俯視圖)Fig.4 Liquid sloshing along lateral direction (side view)

圖5 液罐車內(nèi)液體的縱向晃動(dòng)Fig.5 Liquid sloshing along longitudinal direction in tank

4)仿真計(jì)算液罐車內(nèi)液體的橫向和縱向晃動(dòng)過程,發(fā)現(xiàn)液體自由表面形狀呈非線性變化,不同斷面處的高度和形狀不同；液體的橫向、縱向晃動(dòng)必然引起液罐車的質(zhì)心發(fā)生改變,從而影響輪重減載和車鉤力變化.

5)本文限于研究罐體和液體之間的耦合運(yùn)動(dòng),暫未考慮整車模型,也未考慮罐體內(nèi)含隔板情況.

1Abramson H N. The dynamic behavior of liquids in moving containers-with applications to space vehicle technology. Technical Report NASA SP-106, 1966. San Antonio,TX: Southwest Research Institute

2Ranganathan R,Ying Y, Miles J B. Analysis of fluid slosh in partially filled tanks and their impact on the directional response of tank vehicles. SAE Technical Paper,(No. 932942), 1993

3Ranganathan R,Ying Y, MilesJ B. Development of a mechanical analogy model to predict the dynamic behavior of liquids in partially filled tank vehicles,SAE Technical Paper,(No. 942307), 1994

4Salem M I. Rollover stability of partially filled heavy-duty elliptical tankers using trammel pendulums to simulate fluid sloshing[Ph.D Thesis]. West Virginia: West Virginia University, 2000

5Zheng X, Li X, Ren Y. Equivalent mechanical model for lateral liquid sloshing in partially filled tank vehicles.MathematicalProblemsinEngineering, 2012:123～132

6Celebi M S, Akyildiz H. Nonlinear modeling of liquid sloshing in a moving rectangular tank.OceanEngineering, 2001,29(12):1527～1553

7Ibrahim R A, Pilipchuk V N, Ikeda T. Recent advances in liquid sloshing dynamics.AppliedMechanicsReviews, 2001,54(2):133～199

8Ibrahim R A. Liquid Sloshing Dynamics: Theory and Applications. Cambridge University Press, 2005

9Aliabadi S, Johnson A, Abedi J. Comparison of finite element and pendulum models for simulation of sloshing.Computers&Fluids, 2003,32(4):535～545

10 Rebouillat S, Liksonov D. Fluid-structure interaction in partially filled liquid containers: a comparative review of numerical approaches.Computers&Fluids, 2010,39(5):739～746

11 Bogomaz G I, Markova O M, Chernomashentseva Y G. Mathematical modelling of vibrations and loading of railway tanks taking into account the liquid cargo mobility.VehicleSystemDynamics, 1998,30(3-4):285～294

12 Vera C, Paulin J, Suarez B, et al. Simulation of freight trains equipped with partially filled tank containers and related resonance phenomenon.ProceedingsoftheInstitutionofMechanicalEngineers,PartF:JournalofRailandRapidTransit, 2005,219(4):245～259

13 Wang L, Octavio J R J, Wei C, et al. Low order continuum-based liquid sloshing formulation for vehicle system dynamics.JournalofComputationalandNonlinearDynamics, 2015,10(2):021022

14 Wei C, Wang L, Shabana A A. A total Lagrangian ANCF liquid sloshing approach for multibody system applications.JournalofComputationalandNonlinearDynamics, 2016,10(5):051014

15 馬馳騁,張希農(nóng),柳征勇等. 變質(zhì)量貯箱類流固耦合系統(tǒng)的振動(dòng)響應(yīng)及時(shí)頻特性分析. 振動(dòng)與沖擊, 2014,33(21):166～171 (Ma C C, Zhang X N, Liu Z Y, et al. Dynamic responses and time-frequency feature analysis for a fluid-structure coupling system with a variable mass tank.JournalofVibrationandShock, 2014,33(21):166～171 (in Chinese))

16 劉富,童明波,陳建平. 儲(chǔ)箱內(nèi)液體晃動(dòng)動(dòng)力學(xué)分析及防晃結(jié)構(gòu)優(yōu)化. 計(jì)算機(jī)應(yīng)用與軟件, 2011,28(12):202～205 (Liu F, Tong M B, Chen J P. Inside-container liquid sloshing dynamics analysis and sloshing suppression structure optimization.ComputerApplicationsandSoftware, 2011,28(12):202～205 (in Chinese))

17 王勇. 考慮液體晃動(dòng)的三大件轉(zhuǎn)向架罐車耦合系統(tǒng)動(dòng)力學(xué)性能研究[博士學(xué)位論文]. 成都:西南交通大學(xué), 2004 (Wang Y. Study on the coupling system dynamics for three-piece bogie tank cars taking into account the liquid sloshing[Ph.D Thesis]. Chengdu, Southwest Jiaotong University, 2004 (in Chinese))

18 盧軍. 任意充液比油罐車液體晃動(dòng)及整車橫向穩(wěn)定性研究[碩士學(xué)位論文]. 成都:西南交通大學(xué), 2009 (Lu J. Liquid sloshing and transverse stability analysis of tank truck with arbitrary proportional liquid[Master Thesis]. Chengdu:Southwest Jiaotong University, 2009 (in Chinese))

19 Shabana A A, Hussien H A, Escalona J L. Application of the absolute nodal coordinate formulation to large rotation and large deformation problems.JournalofMechanicalDesign, 1998,120(2):188～195

20 田強(qiáng),劉鋮,李培等. 多柔體系統(tǒng)動(dòng)力學(xué)研究進(jìn)展與挑戰(zhàn). 動(dòng)力學(xué)與控制學(xué)報(bào), 2017,15(5):385～405 (Tian Q, Liu C, Li P, et al. Advances and challenges in dynamics of flexible multibody systems.JournalofDynamicsandControl, 2017,15(5):385～405 (in Chinese))

21 Olshevskiy A, Dmitrochenko O, Kim C W. Three-dimensional solid brick element using slopes in the absolute nodal coordinate formulation.JournalofComputationalandNonlinearDynamics, 2014,9(2):021001

22 Shabana A A. Computational continuum mechanics,3rd ed. Cambridge, UK: Cambridge University Press, 2017