分數Maxwell黏彈性膠體阻尼緩沖器振蕩流模型研究

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

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

0 引言

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

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

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

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

黏彈性膠體材料在阻尼緩沖器中運動時具有較強的非線性和黏彈性，呈現非牛頓流體特性，而采用以準態特性為參數的分數Maxwell模型模擬黏彈性流體在圓管內和平板間的振蕩流鮮有報道，同時應用其預測滯回曲線以指導阻尼緩沖器設計的報道也較少。因此，本文針對黏彈性膠體在阻尼緩沖器阻尼孔和節流間隙中振蕩流問題，設定振蕩邊界條件，構建以準態特性為參數的分數Maxwell模型，并與牛頓流體模型進行對比，研究黏彈性流體在圓管內和兩平板間的振蕩流動特性。通過實驗驗證分數Maxwell模型模擬黏彈性膠體阻尼緩沖器振蕩流的準確性。

1 分數Maxwell振蕩流模型構建

1.1 分數Maxwell模型

廣義黏彈性模型是將彈簧與黏壺串聯看作一個機構元件，其本構關系可寫為

(1)

采用Maxwell模型串聯兩個彈簧壺機構元件(κ,α)和(ψ,β)(見圖1)，其中ψ為準態特性(Pa·sβ)，β為分數階指數。由Maxwell模型關系可得

σa=σb=σ，

(2)

εa+εb=ε,

(3)

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

圖1 分數Maxwell模型Fig.1 Fractional Maxwell model

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

(4)

即分數Maxwell模型為

(5)

1.2 結構及簡化

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

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

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

圖3 阻尼緩沖器結構與黏彈性膠體的運動關系簡化Fig.3 Simplified kinematic relation between shock absorber structure and viscoelastic elastomer

1.3 控制方程

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

ma=∑F,

(6)

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

與坐標系無關的運動方程矢量表達式為

(7)

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

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

(8)

即振蕩流運動方程可寫為

(9)

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

(10)

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

(11)

(12)

由Riemann-Liouville分數階微積分公式可知:

(13)

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

(14)

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

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

1)阻尼孔為

(15)

2)節流間隙為

(16)

1.4 無量綱分析

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

(17)

式中：η為準態特性無量綱系數。

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

(18)

(19)

對應的阻尼孔和節流間隙的無量綱初始和邊界條件為

(20)

(21)

1.5 數值求解

(22)

(23)

o()為高階無窮小，則

(24)

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

則

(25)

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

(26)

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

(27)

根據分數階積分定義:

(28)

(29)

(30)

(31)

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

(32)

(33)

2 數值結果及討論

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

2.1 節流間隙數值模擬

圖4 ω=2時分數Maxwell模型節流間隙流速分布(α=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)

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

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

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

2.2 阻尼孔數值模擬

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

圖8 ω=2時分數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)流速分布三維圖，可以看出分數Maxwell模型流速分布形狀呈柱塞形，而牛頓流體模型流速分布呈拋物線狀。對比圖10(b)和圖11(b)兩種模型二維流速分布圖可知，牛頓流體模型流速分布振蕩傳遞性明顯較分數Maxwell模型弱，因此其非線性也較弱。結合圖8～圖11可以判斷分數Maxwell模型對頻率的依賴性更強，振蕩頻率不同流速分布也不同，而牛頓流體模型則對頻率依賴性較弱，該結果與節流間隙類似。

圖10 ω=6時分數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 參數變化影響分析

圖12 ω=6時各參數變化對分數Maxwell模型阻尼孔y=0.5處應力- 應變率分布影響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處應力- 應變率分布影響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 實驗驗證

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

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

表1 黏彈性膠體阻尼緩沖器結構尺寸

Tab.1 Basic parameters of viscoelastic elastomer shock absorber

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

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

表2 膠體材料及模型參數[27]

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

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

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

表3 實驗與數值滯回曲線結果及誤差對比

Tab.3 Test and simulated results errors of hysteretic curves

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

4 結論

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

1)分數Maxwell模型流速分布較牛頓流體模型非線性強，且具有較強的頻率依賴性；不同參數下，分數Maxwell模型應力- 應變率曲線均呈橢圓形，而牛頓流體模型應力- 應變率曲線呈直線。

2)隨著參數α、β和η的增加，分數Maxwell模型的應力- 應變率曲線橢圓長軸沿逆時針旋轉。

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

因此，分數Maxwell模型可模擬黏彈性膠體在阻尼緩沖器阻尼孔和節流間隙中振蕩流的流動特性，用以預測黏彈性膠體阻尼緩沖器滯回曲線的形狀及變化趨勢，從而指導黏彈性膠體阻尼緩沖器的設計。