王之千， 毛保全， 朱銳， 白向華， 韓小平

(陸軍裝甲兵學(xué)院 兵器與控制系, 北京 100072)

0 引言

廣義上振蕩流指流體在周期振蕩邊界條件下的流動特性。國內(nèi)外眾多學(xué)者根據(jù)不同的研究背景，對振蕩流進(jìn)行了相應(yīng)研究[1-2]。Mentzoni等[3]基于Navier-Stokes方程提出了一種與海洋應(yīng)用有關(guān)的振蕩流動條件下多孔板水動力二維黏性流場模型，求得了解析解，并通過實驗驗證了模型的準(zhǔn)確性。Fu等[4]采用實驗方法測量了通過多孔通道的速度、測試段的壓降和沿?zé)岜砻娴臏囟?，研究了穩(wěn)態(tài)和振蕩氣流作用下多孔通道的傳熱特性。Ni等[5]將實驗和仿真相結(jié)合提出了一種定量描述含多孔介質(zhì)振蕩流動特性的方法，從而指導(dǎo)基于振蕩流器件的設(shè)計。Pamuk[6]通過多組實驗研究了振蕩流動條件下流體的傳熱問題，得到了一種新的振蕩流努塞爾數(shù)相關(guān)關(guān)系，為估計振蕩流傳熱提供可靠方法。

黏彈性膠體材料屬于非牛頓流體，外力作用下在阻尼緩沖器中作周期往復(fù)運動時黏彈性膠體的流動過程可看作振蕩流。然而上述振蕩流的研究背景、材料及邊界條件與黏彈性膠體阻尼緩沖器均不相同，且主要以多孔、微孔流動或傳熱導(dǎo)熱為主，加之黏彈性流體獨特的非線性，應(yīng)用于黏彈性膠體阻尼緩沖器中較為困難。

隨著分?jǐn)?shù)階導(dǎo)數(shù)的提出，分?jǐn)?shù)Maxwell模型被證明可以較好地模擬黏彈性流體的流動過程[7-9]。Yan等[10]等將分?jǐn)?shù)Maxwell方法用于研究滾動運動中的管道振蕩流動，以模擬船上管道受到海水振蕩時的情況，得到了滾動運動中管內(nèi)速度梯度及變化情況。Tan等[11-12]利用連續(xù)分?jǐn)?shù)導(dǎo)數(shù)離散拉普拉斯逆變換，研究了黏彈性流體在平面突然運動及兩平行板間的分?jǐn)?shù)Maxwell流體模型，并求得了速度和應(yīng)力的解析解。Zhang等[13]研究了黏彈性納米流體在多孔介質(zhì)中的二維非定常邊界層流動和傳熱傳質(zhì)問題，提出了一種含布朗擴(kuò)散和熱遷移的時間- 空間分?jǐn)?shù)Maxwell導(dǎo)熱模型,分析了參數(shù)對速度、溫度和濃度分布的影響。上述研究從數(shù)學(xué)角度很好地描述了黏彈性流體的流動特性，但值得注意地是，這些表達(dá)式模型參數(shù)往往通過數(shù)值方法得出，并不始終在物理和工程上保持合理性[14]。

準(zhǔn)態(tài)特性的提出解決了上述問題，它將數(shù)學(xué)與物理意義相結(jié)合，使模型參數(shù)與實際研究對象聯(lián)系更為緊密[15]。雖然準(zhǔn)態(tài)特性已提出半個多世紀(jì)，但直到近幾年才在分?jǐn)?shù)Maxwell模型中得到廣泛應(yīng)用[16]。Holder等[17]通過研究膠原凝膠的黏彈性力學(xué)性能，證明了分?jǐn)?shù)Maxwell模型可以很好地描述成熟膠原凝膠力學(xué)性能的時間依賴性。Jaishankar等[18-19]構(gòu)建了含準(zhǔn)態(tài)特性的剪切變形分?jǐn)?shù)Maxwell模型，通過與實驗數(shù)據(jù)對比，表明在一定時間尺度范圍內(nèi)，模型預(yù)測的長時間冪律響應(yīng)與實驗數(shù)據(jù)具有很好的一致性，并證明了準(zhǔn)態(tài)特性可準(zhǔn)確描述復(fù)雜流變界面，有助于確定材料對其他變形模式的響應(yīng)。另外，文獻(xiàn)[20-22]中對準(zhǔn)態(tài)特性引入分?jǐn)?shù)Maxwell模型等分?jǐn)?shù)階本構(gòu)模型進(jìn)行了詳細(xì)研究。

黏彈性膠體材料在阻尼緩沖器中運動時具有較強(qiáng)的非線性和黏彈性，呈現(xiàn)非牛頓流體特性，而采用以準(zhǔn)態(tài)特性為參數(shù)的分?jǐn)?shù)Maxwell模型模擬黏彈性流體在圓管內(nèi)和平板間的振蕩流鮮有報道，同時應(yīng)用其預(yù)測滯回曲線以指導(dǎo)阻尼緩沖器設(shè)計的報道也較少。因此，本文針對黏彈性膠體在阻尼緩沖器阻尼孔和節(jié)流間隙中振蕩流問題，設(shè)定振蕩邊界條件，構(gòu)建以準(zhǔn)態(tài)特性為參數(shù)的分?jǐn)?shù)Maxwell模型，并與牛頓流體模型進(jìn)行對比，研究黏彈性流體在圓管內(nèi)和兩平板間的振蕩流動特性。通過實驗驗證分?jǐn)?shù)Maxwell模型模擬黏彈性膠體阻尼緩沖器振蕩流的準(zhǔn)確性。

1 分?jǐn)?shù)Maxwell振蕩流模型構(gòu)建

1.1 分?jǐn)?shù)Maxwell模型

廣義黏彈性模型是將彈簧與黏壺串聯(lián)看作一個機(jī)構(gòu)元件，其本構(gòu)關(guān)系可寫為

(1)

采用Maxwell模型串聯(lián)兩個彈簧壺機(jī)構(gòu)元件(κ,α)和(ψ,β)(見圖1)，其中ψ為準(zhǔn)態(tài)特性(Pa·sβ)，β為分?jǐn)?shù)階指數(shù)。由Maxwell模型關(guān)系可得

σa=σb=σ，

(2)

εa+εb=ε,

(3)

式中：σa和εa分別為彈簧壺元件(κ,α)的應(yīng)力與應(yīng)變；σb和εb分別為彈簧壺元件(ψ,β)的應(yīng)力與應(yīng)變；σ和ε分別為Maxwell模型的總應(yīng)力與總應(yīng)變。

圖1 分?jǐn)?shù)Maxwell模型Fig.1 Fractional Maxwell model

(1)式代入(2)式、(3)式，得

(4)

即分?jǐn)?shù)Maxwell模型為

(5)

1.2 結(jié)構(gòu)及簡化

某型黏彈性膠體阻尼緩沖器結(jié)構(gòu)如圖2所示。阻尼緩沖器支座用于連接阻尼緩沖器和安裝機(jī)構(gòu)，起一定承力作用。將黏彈性膠體材料通過黏彈性膠體注入器注入阻尼緩沖器腔室中。根據(jù)需要調(diào)整黏彈性膠體材料的填充量，使阻尼緩沖器具有一定的預(yù)壓力，稱為初始預(yù)緊力。當(dāng)阻尼緩沖器受到外界壓力或沖擊大于初始預(yù)緊力時，活塞桿推動活塞壓縮阻尼緩沖器腔室內(nèi)的黏彈性膠體材料。膠體材料被迫流過節(jié)流間隙和活塞上的阻尼孔，產(chǎn)生黏滯力，阻礙活塞前行，這一過程中外力轉(zhuǎn)化為熱能和勢能，消耗和存儲。當(dāng)外力被撤銷時，沖擊力小于初始預(yù)緊力，黏彈性膠體材料自行膨脹，釋放壓縮過程中儲存的勢能，將活塞推回到初始位置，而通過阻尼孔和節(jié)流間隙被擠入前腔體的膠體材料流回后腔體，等待下次沖擊或壓力。

圖2 某型黏彈性膠體阻尼緩沖器示意圖Fig.2 Schematic diagram of a viscoelastic elastomer shock absorber

根據(jù)黏彈性膠體阻尼緩沖器振蕩流實際工況，黏彈性膠體在阻尼緩沖器阻尼孔中的運動可簡化為流體中心以ucos(ωt)(u為一個恒定的振幅，表示速度；ω為振蕩頻率)的流速在圓管內(nèi)振蕩流動，而靠近管壁的流體流速為0 m/s(見圖3(a)所示，即流體在兩平行板間以ucos(ωt)的速度左右振蕩流動，而兩平板始終保持不動)；黏彈性膠體在節(jié)流間隙中振蕩流的運動可簡化為一板不動、一板以ucos(ωt)的速度振蕩運動(見圖3(b)所示，即上板以ucos(ωt)的速度左右振蕩流動，下板始終保持不動，流體由于粘性在上板的帶動下左右振蕩流動)。圖3中：L為阻尼孔或節(jié)流間隙間距；p為一個恒定的振幅，表示壓力。

圖3 阻尼緩沖器結(jié)構(gòu)與黏彈性膠體的運動關(guān)系簡化Fig.3 Simplified kinematic relation between shock absorber structure and viscoelastic elastomer

1.3 控制方程

流體運動中總動量的時間變化率等于其自身體積力和相互作用表面力的總和。將加速度a看作單位質(zhì)量流體動量隨時間的變化率，根據(jù)動量平衡率的牛頓運動定律可知：

ma=∑F,

(6)

式中：m為質(zhì)量；F為體積力和表面力總和。

與坐標(biāo)系無關(guān)的運動方程矢量表達(dá)式為

(7)

式中：f為單位體積力；ρ為黏彈性流體密度；σ為應(yīng)力張量。

考慮阻尼緩沖器中黏彈性流體介質(zhì)運動時黏性較大、總質(zhì)量較小，體積力遠(yuǎn)小于黏性力，因此將體積力忽略不計。則有

(8)

即振蕩流運動方程可寫為

(9)

其中考慮沿x軸方向外加振蕩的壓力梯度[23-24]，如圖3所示，其形式為

(10)

式中：u(t,y)和σ(t,y)分別表示在時刻t空間位置y處的流速和應(yīng)力。

(11)

(12)

由Riemann-Liouville分?jǐn)?shù)階微積分公式可知:

(13)

式中：t>0，θ>0，?=「θ?表示不小于θ的最小整數(shù)，且

(14)

式中：Γ(z)是伽馬函數(shù)。

假定邊界無滑移，由圖3可得相應(yīng)的初始和邊界條件：

1)阻尼孔為

(15)

2)節(jié)流間隙為

(16)

1.4 無量綱分析

采用Π定理引入如下無量綱變量：

(17)

式中：η為準(zhǔn)態(tài)特性無量綱系數(shù)。

(17)式的無量綱變量分別代入(11)式、(12)式中，得到分?jǐn)?shù)Maxwell振蕩流無量綱方程(以下變量均為無量綱變量，為簡便起見，省略無量綱標(biāo)記“*”)：

(18)

(19)

對應(yīng)的阻尼孔和節(jié)流間隙的無量綱初始和邊界條件為

(20)

(21)

1.5 數(shù)值求解

(22)

(23)

o()為高階無窮小，則

(24)

式中：r為高階項系數(shù)；t0=0；h0=1；hi=(i+1)1-h-i1-h,i=0,1,…,k.

則

(25)

式中：高階項系數(shù)r=2-h；hj=(j+1)1-h-j1-h,j=0,1,…,k. 則

(26)

同理，利用中心差分定理：

(27)

根據(jù)分?jǐn)?shù)階積分定義:

(28)

(29)

(30)

(31)

另外，牛頓流體模型根據(jù)(1)式可寫為

(32)

(33)

2 數(shù)值結(jié)果及討論

將模型和邊界條件代入數(shù)學(xué)仿真軟件MATLAB中，為保證結(jié)果收斂，選取時間步長Δt=0.001、空間步長Δy=0.025. 通過數(shù)值計算，對比分?jǐn)?shù)Maxwell模型和牛頓流體模型的流動特性，選取基本模型參數(shù)為α=0.8,β=0.3,η=0.2.

2.1 節(jié)流間隙數(shù)值模擬

圖4 ω=2時分?jǐn)?shù)Maxwell模型節(jié)流間隙流速分布(α=0.8,β=0.3,η=0.2)Fig.4 Flow velocity distributions of the fractional Maxwell model in the gap for ω=2 (α=0.8,β=0.3,η=0.2)

分?jǐn)?shù)Maxwell模型和牛頓流體模型在振蕩頻率ω為2時節(jié)流間隙振蕩流流速分布數(shù)值結(jié)果分別如圖4和圖5所示，其中，u為無量綱速度，t為無量綱時間，y為無量綱位移。結(jié)合圖4和圖5可知，兩種模型均發(fā)生周期振蕩，且靠近運動板的流體(y=0處)流速變化幅度最大，而靠近停止板的流體(y=1處)流速始終保持為0. 對比圖4(a)和圖5(a)兩種模型節(jié)流間隙振蕩流流速分布三維圖可以看出，兩種模型在平板間隨時間變化的流速分布趨勢類似，均呈平緩的梯度分布，其振蕩幅度逐漸降低，但分?jǐn)?shù)Maxwell模型流速分布較牛頓流體模型流速分布非線性更強(qiáng)。圖4(b)和圖5(b)分別為兩種模型在一個周期內(nèi)隨時間增加時的流速分布二維圖。從圖4和圖5中可以看出在一個周期內(nèi)隨著時間的增長，兩種模型在y=0處的流體流速振蕩幅度最大，并通過振蕩將流速傳遞至相鄰位置，對比可以看出分?jǐn)?shù)Maxwell模型流速分布較牛頓流體模型非線性明顯，而牛頓流體模型在流速最大時發(fā)展近似呈直線。

圖6和圖7分別為ω=6時分?jǐn)?shù)Maxwell模型和牛頓流體模型節(jié)流間隙振蕩流流速分布。對比圖6(a)和圖7(a)可以看出：分?jǐn)?shù)Maxwell模型的流速分布隨y的增加由較大的流速降至較低流速后降為0，變化幅度較大；而牛頓流體模型流速下降較為平緩，有明顯過度。圖6(b)和圖7(b)為一個周期內(nèi)兩種模型隨時間增加的二維流速分布，可以看出分?jǐn)?shù)Maxwell模型較牛頓流體模型具有較強(qiáng)非線性，振蕩傳遞性較為明顯。結(jié)合圖4～圖7可以看出：振蕩頻率ω對分?jǐn)?shù)Maxwell模型影響較大，而對牛頓流體模型影響較小，原因在于分?jǐn)?shù)Maxwell模型具有較強(qiáng)的非線性，后運動的流體流速受前一刻流體流速和頻率影響均較為明顯；而牛頓流體模型中后運動的流體流速僅受前一刻流體流速的影響，頻率對其影響較小，這也是圖7(a)中牛頓流體模型相比Maxwell模型流速分布下降平緩的原因。因此，分?jǐn)?shù)Maxwell模型模擬的流體較牛頓流體模型具有較強(qiáng)的非線性。

圖7 ω=6時牛頓流體模型節(jié)流間隙流速分布(η=0.2)Fig.7 Flow velocity distributions of the Newtonian fluid model in the gap for ω=6 (η=0.2)

2.2 阻尼孔數(shù)值模擬

圖8和圖9分別為ω=2時分?jǐn)?shù)Maxwell模型和牛頓流體模型阻尼孔振蕩流流速分布數(shù)值結(jié)果，從中可以看出兩種模型均發(fā)生周期振蕩，且在阻尼孔內(nèi)的流速分布曲線沿y=0.5的軸線對稱。對比圖8和圖9中ω=2時阻尼孔振蕩流流速分布可知，兩種模型的流速分布趨勢類似，一個周期內(nèi)同一時刻分?jǐn)?shù)Maxwell模型較牛頓流體模型的非線性稍強(qiáng)，當(dāng)振幅較大時牛頓流體模型流速分布在0

圖8 ω=2時分?jǐn)?shù)Maxwell模型阻尼孔流速分布(α=0.8,β=0.3,η=0.2)Fig.8 Flow velocity distributions of the fractional Maxwell model in the orifice for ω=2 (α=0.8,β=0.3,η=0.2)

圖9 ω=2時牛頓流體模型阻尼孔流速分布(η=0.2)Fig.9 Flow velocity distributions of the Newtonian fluid model in theorifice for ω=2 (η=0.2)

圖10和圖11分別為ω=6時兩種模型的阻尼孔流速分布。對比圖10(a)和圖11(a)流速分布三維圖，可以看出分?jǐn)?shù)Maxwell模型流速分布形狀呈柱塞形，而牛頓流體模型流速分布呈拋物線狀。對比圖10(b)和圖11(b)兩種模型二維流速分布圖可知，牛頓流體模型流速分布振蕩傳遞性明顯較分?jǐn)?shù)Maxwell模型弱，因此其非線性也較弱。結(jié)合圖8～圖11可以判斷分?jǐn)?shù)Maxwell模型對頻率的依賴性更強(qiáng)，振蕩頻率不同流速分布也不同，而牛頓流體模型則對頻率依賴性較弱，該結(jié)果與節(jié)流間隙類似。

圖10 ω=6時分?jǐn)?shù)Maxwell模型阻尼孔流速分布(α=0.8,β=0.3,η=0.2)Fig.10 Flow velocity distributions of the fractional Maxwell model in the orifice for ω=6 (α=0.8,β=0.3,η=0.2)

圖11 ω=6時牛頓流體模型阻尼孔流速分布(η=0.2)Fig.11 Flow velocity distributions of the Newtonian fluid model in the damping orifice for ω=6 (η=0.2)

2.3 參數(shù)變化影響分析

圖12 ω=6時各參數(shù)變化對分?jǐn)?shù)Maxwell模型阻尼孔y=0.5處應(yīng)力- 應(yīng)變率分布影響Fig.12 Influences of parameters on the stress-strain rate distributions of the fractional Maxwell model in y=0.5 of the damping orifice for ω=6

圖13 ω=6時不同η對牛頓流體模型阻尼孔y=0.5處應(yīng)力- 應(yīng)變率分布影響Fig.13 Influence of η on the stress-strain rate curves of the Newtonian fluid model in y=0.5 of damping orifice for ω=6

4 實驗驗證

為檢驗分?jǐn)?shù)Maxwell模型對黏彈性膠體阻尼緩沖器滯回曲線的預(yù)測性，采用MTS 810型電液伺服疲勞實驗機(jī)(如圖14所示)加載正弦波激勵模擬阻尼緩沖器的實際工況。為簡化模型運算，僅考慮阻尼孔的影響，黏彈性膠體阻尼緩沖器實驗樣機(jī)采用單出桿阻尼孔式，其結(jié)構(gòu)如圖2去除節(jié)流間隙所示，具體尺寸參數(shù)如表1所示。

圖14 MTS 810型電液伺服疲勞實驗機(jī)Fig.14 MTS 810 electro-hydraulic servo fatigue testing machine

表1 黏彈性膠體阻尼緩沖器結(jié)構(gòu)尺寸

Tab.1 Basic parameters of viscoelastic elastomer shock absorber

阻尼孔直徑/mm活塞直徑/mm最大行程/mm2.52012

為使黏彈性膠體阻尼緩沖器實際工況與模型初始及邊界條件相近，采用文獻(xiàn)[26]的實驗結(jié)果，當(dāng)位移小于等于1 mm時，邊界滑移現(xiàn)象不明顯。因此，選擇加載頻率振幅為0.3 mm，由于阻尼緩沖器無法承受拉應(yīng)力，實驗過程中采用疲勞實驗機(jī)先將阻尼緩沖器壓縮0.3 mm，后分別加載8 Hz、10 Hz的正弦波激勵頻率及0.3 mm的振蕩振幅，并采集阻尼緩沖器滯回曲線。表2膠體材料及模型參數(shù)代入分?jǐn)?shù)Maxwell模型中，得到應(yīng)力- 應(yīng)變率曲線，對應(yīng)變率進(jìn)行積分可得模型的應(yīng)力- 應(yīng)變曲線，將其看作模型所模擬的阻尼緩沖器滯回曲線。

表2 膠體材料及模型參數(shù)[27]

圖15 頻率為8 Hz時實驗滯回曲線和數(shù)值模擬滯回曲線對比Fig.15 Test and simulated hysteretic curves at 8 Hz

圖15和圖16分別為8 Hz和10 Hz時疲勞實驗機(jī)實驗所得滯回曲線和分?jǐn)?shù)Maxwell模型模擬所得滯回曲線，實驗和模擬結(jié)果如表3所示。對比圖15和圖16不同頻率下分?jǐn)?shù)Maxwell模型模擬的黏彈性膠體緩沖器滯回曲線與實驗所得滯回曲線可知：數(shù)值模擬的滯回曲線與實驗所得滯回曲線具有相同的形狀趨勢，均呈橢圓形，其能量吸收率、最大位移、最小力和最大力與實驗結(jié)果相近，相對于實驗結(jié)果平均誤差分別為3.60%、2.55%、1.81%和3.77%. 數(shù)值模擬阻尼緩沖器滯回曲線時，能量吸收率和曲線形狀趨勢是考察模型準(zhǔn)確與否的主要指標(biāo)。分?jǐn)?shù)Maxwell模型模擬阻尼緩沖器滯回曲線形狀趨勢與實驗結(jié)果相近，且能量吸收率平均誤差為3.60%. 因此，分?jǐn)?shù)Maxwell模型可用以預(yù)測黏彈性膠體阻尼緩沖器滯回曲線的形狀及變化趨勢，從而指導(dǎo)黏彈性膠體阻尼緩沖器的設(shè)計。加之其非線性較強(qiáng)，可很好地描述黏彈性膠體的黏彈性特性。

圖16 頻率為10 Hz時實驗滯回曲線和數(shù)值模擬滯回曲線對比Fig.16 Test and simulated hysteretic curves at 10 Hz

表3 實驗與數(shù)值滯回曲線結(jié)果及誤差對比

Tab.3 Test and simulated results errors of hysteretic curves

不同激勵頻率滯回曲線結(jié)果最大力/kN最小力/kN最大位移/mm能量吸收率/%8Hz時實驗結(jié)果2.49750.19390.5879.298Hz時數(shù)值模擬結(jié)果2.45390.18870.5775.9110Hz時實驗結(jié)果2.48570.23600.5479.3610Hz時數(shù)值模擬結(jié)果2.36160.23820.5277.03數(shù)值模擬相對實驗結(jié)果誤差最大力誤差/%最小力誤差/%最大位移誤差/%能量吸收率誤差/%8Hz時1.752.681.724.2610Hz時4.990.933.372.94平均誤差3.371.812.553.60

4 結(jié)論

本文根據(jù)黏彈性膠體阻尼緩沖器受力時周期往復(fù)運動的實際工況，簡化黏彈性膠體在阻尼緩沖器阻尼孔和節(jié)流間隙中的振蕩流動過程，設(shè)定振蕩流初始和邊界條件，提出并構(gòu)建了含準(zhǔn)態(tài)特性參數(shù)用于研究黏彈性膠體阻尼緩沖器振蕩流的分?jǐn)?shù)Maxwell模型，采用有限差分法求得數(shù)值解。通過與牛頓流體模型結(jié)果進(jìn)行對比，分析了不同振蕩頻率下的流速分布及各參數(shù)對應(yīng)力- 應(yīng)變率曲線的影響。采用疲勞實驗機(jī)加載頻率模擬阻尼緩沖器實際工況，采集滯回曲線并與分?jǐn)?shù)Maxwell模型模擬結(jié)果進(jìn)行對比，驗證了模型的可預(yù)測性。數(shù)值及實驗結(jié)果表明：

1)分?jǐn)?shù)Maxwell模型流速分布較牛頓流體模型非線性強(qiáng)，且具有較強(qiáng)的頻率依賴性；不同參數(shù)下，分?jǐn)?shù)Maxwell模型應(yīng)力- 應(yīng)變率曲線均呈橢圓形，而牛頓流體模型應(yīng)力- 應(yīng)變率曲線呈直線。

2)隨著參數(shù)α、β和η的增加，分?jǐn)?shù)Maxwell模型的應(yīng)力- 應(yīng)變率曲線橢圓長軸沿逆時針旋轉(zhuǎn)。

3)分?jǐn)?shù)Maxwell模型模擬的滯回曲線與實驗所得滯回曲線具有相同的形狀趨勢，均呈橢圓形，且相對于實驗結(jié)果能量吸收率平均誤差為3.60%.

因此，分?jǐn)?shù)Maxwell模型可模擬黏彈性膠體在阻尼緩沖器阻尼孔和節(jié)流間隙中振蕩流的流動特性，用以預(yù)測黏彈性膠體阻尼緩沖器滯回曲線的形狀及變化趨勢，從而指導(dǎo)黏彈性膠體阻尼緩沖器的設(shè)計。