分數Maxwell黏彈性膠體阻尼緩沖器啟動流模型

王之千，毛保全，馮帥，楊雨迎，朱銳

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

0 引言

黏彈性膠體材料屬于非牛頓流體，是一種半液半固黏稠狀高分子聚合物，具有良好的彈性、可壓縮性及流動性等優異的特點。將其裝入特殊的密閉空腔機械結構中，起到減振、緩沖、吸能等效果，該裝置即為黏彈性膠體阻尼緩沖器。黏彈性膠體在阻尼緩沖器中的啟動流形式根據緩沖器機械結構可分為間隙流和孔隙流兩種。

由于黏彈性流體獨特的非線性，國內外眾多學者根據不同的研究背景和黏彈性材料，選用不同的數學模型和本構關系研究阻尼緩沖器中黏彈性流體的流動特性。研究主要分為兩類：一類為經驗或半經驗公式，這類研究公式推導基于大量試驗數據，所推出的經驗公式根據研究對象的不同而不同，使用限定條件較多，適用范圍較為局限。另一類為經典模型及其衍生模型，其中以Maxwell模型最為典型。尤其隨著分數階導數的提出，分數Maxwell模型被證明可以較好地模擬黏彈性流體的流動過程[1-5]。Yang等[6]采用分數Maxwell模型研究了黏彈性流體在無限長管道中恒定壓力梯度驅動下從靜止開始的單向啟動流動情況，求出了精確解。Yao[7]在逆行階梯Maxwell模型基礎上通過張量和不同分數階導數定義將分數階黏壺模型用于模擬大變形非線性黏彈性流體的流動。上述研究從數學角度很好地描述了黏彈性流體的流動特性，但值得注意的是，這些表達式模型參數往往通過數值方法得出，并不始終在物理和工程上保持合理性[8]。

圖1 彈簧壺元件原理圖Fig.1 Schematic diagram of a spring-pot as an element

黏彈性膠體材料在阻尼緩沖器中運動時呈現非牛頓流體特性，具有較強的黏彈性，可采用分數Maxwell模型進行研究。但上述文獻中均求得模型的解析解，然而分數階微積分方程求解析解難度較大、過程較為復雜，計算機模擬時運算速度較慢；而數值解仿真模擬易于實現。并且采用以準態特性為參數的分數Maxwell模型模擬黏彈性流體在圓管內和平板間加速流動的啟動流鮮有報道，同時在工程上用于研究阻尼緩沖器節流間隙和阻尼孔的分數Maxwell啟動流模型也報道較少。

因此，本文針對黏彈性膠體在阻尼緩沖器阻尼孔和節流間隙中的啟動流問題，構建以準態特性為參數的分數Maxwell模型，采用有限差分法求得模型的數值解，并設計啟動流測試試驗系統，驗證分數Maxwell模型模擬黏彈性膠體阻尼緩沖器啟動流的準確性。

1 分數Maxwell啟動流模型構建

1.1 分數Maxwell模型

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

(1)

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

σa=σb=σ,

(2)

εa+εb=ε，

(3)

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

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

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

(4)

即分數Maxwell模型為

(5)

1.2 結構及簡化

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

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

為符合黏彈性膠體阻尼緩沖器啟動流時的實際工況，將黏彈性膠體在阻尼緩沖器阻尼孔和節流間隙中啟動流的運動關系簡化至二維空間中，分別看作流體在圓管內加速流動(見圖4(a)，即流體在兩平行板間以U=at的速度流動，而兩平行板始終保持不動)和兩平行板間一板不動、一板加速運動這兩種形式(見圖4(b)，即下板以U=at的速度向右移動，上板始終保持不動，流體由于黏性在下板的帶動下向右流動)，圖中L為阻尼孔或節流間隙間距，U=at表示流速，a為加速度。

圖4 阻尼緩沖器結構與黏彈性膠體的運動關系Fig.4 Motion relationship between shock absorberstructure and viscoelastic elastomer

1.3 控制方程

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

Ma=∑F,

(6)

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

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

(7)

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

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

(8)

流速只在y軸方向有分布，則

(9)

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

(10)

(11)

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

(12)

(13)

Γ(z)是伽馬函數。

1.4 初始值和邊界條件

當t≤0時，u(t)=0；當t>0時，阻尼孔啟動流中心點的流速以加速度a逐漸增大，而靠近管壁的流速始終為0 m/s. 在節流間隙中，靠近運動板的流體流速隨運動板以加速度a逐漸增大，而靜止板附近的流速始終為0 m/s. 因此，初始條件和邊界條件為

1)阻尼孔：

(14)

2)節流間隙：

(15)

1.5 無量綱分析

采用Π定理分析分數Maxwell模型的無量綱變量，

(16)

將(16)式的無量綱變量分別代入(10)式、(11)式中，得到分數Maxwell啟動流無量綱方程(為簡便起見，省略無量綱標記“*”)：

(17)

(18)

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

(19)

(20)

1.6 數值求解

(21)

(22)

式中：o(Δt)為皮亞諾余項。

則

(23)

式中：t0=0；h0=1；hi=(i+1)1-h-i1-h，i=0,1,…,k.

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

(24)

根據分數階積分定義：

(25)

(26)

(27)

(28)

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

(29)

(30)

2 數值計算結果及分析

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

2.1 阻尼孔數值模擬

圖5 阻尼孔流速分布Fig.5 Flow velocity distribution in orifice at different times

圖5為分數Maxwell模型和牛頓流體模型隨時間增大時阻尼孔流速的變化情況。從圖5中可以看出，隨著時間的不斷增加，分數Maxwell模型流速分布曲線由扁平狀逐漸發展為拋物線狀，符合實際黏彈性流體的流動特性，且流速分布曲線保持著良好的非線性特點，曲線曲率較牛頓流體大。這意味著分數Maxwell模型模擬的流體始終保持著較好的黏彈性。

圖6為阻尼孔中心點的流體應力變化情況。隨著時間的增加，分數Maxwell模型和牛頓流體模型的應力均呈線性增長。由于分數Maxwell模型是通過廣義黏彈性模型推導衍生而得，而牛頓流體模型則是廣義黏彈性模型的一種特殊形式，因此二者應力變化趨勢基本相同。但分數Maxwell模型應力小于牛頓流體模型。

取t=0.8，比較參數變化對分數Maxwell模型流速的影響，如圖7～圖9所示。參數α和β均為分數Maxwell模型的階數，二者的變化情況對流場的影響基本一致，隨著α和β的增加，分數Maxwell模型流體曲率逐漸下降，但變化不大。而ζ相比α和β對曲線曲率影響較大。從圖9中可以看出，當ζ=0.1時流速曲率最大，而隨著ζ的增加，流速曲率逐漸減小，并向牛頓流體發展。ζ與準態特性有關，而準態特性與松弛時間有關，松弛時間越長，材料越接近彈性體。這意味著ζ越大，黏度越小。

圖7 不同α時分數Maxwell模型阻尼孔流速分布Fig.7 Flow velocity distribution of fractional Maxwell model in orifice for different α

圖8 不同β時分數Maxwell模型阻尼孔流速分布Fig.8 Flow velocity distribution of fractional Maxwell model in orifice for different β

圖9 不同ζ時分數Maxwell模型阻尼孔流速分布Fig.9 Flow velocity distribution of fractional Maxwell model in orifice for different ζ

2.2 節流間隙數值模擬

阻尼緩沖器節流間隙啟動流流速分布隨時間的變化情況如圖10所示。從圖10中可以看出，隨著時間的增加，牛頓流體模型流速分布逐漸向直線形狀發展，而分數Maxwell模型始終保持著良好的非線性。進一步證明了阻尼孔流速模擬的結果，分數Maxwell模型非線性良好，具有較強的黏彈性。

圖10 節流間隙流速分布Fig.10 Flow velocity distributions in gap at different times

節流間隙啟動流靠近移動平板一側(y=0)應力分布曲線如圖11所示。由圖11可見，隨著時間的增加，分數Maxwell模型與牛頓流體模型應力均呈線性增長，該結果與阻尼孔啟動流應力分布增長趨勢一致。當t=6時，分數Maxwell模型和牛頓流體模型節流間隙啟動流應力值分別為1.33和1.52，而阻尼孔啟動流應力值分別為2.17和2.54；分數Maxwell模型和牛頓流體模型阻尼孔啟動流的應力分別為節流間隙啟動流應力的1.63和1.67倍。原因在于，阻尼孔啟動流最大速度發生在平板間距中心點處(y=0.5)，而節流間隙啟動流最大速度出現在靠近移動平板一側處(y=0)，二者邊界條件不同。當時間和空間步長相同時，阻尼孔啟動流計算所得速度數據是節流間隙啟動流數據數量的一半。應力通過速度差值計算所得，而分數Maxwell模型和牛頓流體模型應力均呈線性增長，致使阻尼孔啟動流應力值遠比節流間隙啟動流大。雖然減小步長可以使流速數據數量保持一致，但步長的調整會導致差值系數發生變化，因此計算結果將保持不變。

圖11 節流間隙靠近移動平板一側應力分布Fig.11 Stress distribution in gap near the side of the moving plate

3 試驗驗證

試驗系統由伺服電機、電機控制器、玻璃管緩沖器、高速相機、圖像采集平臺和射燈光源組成，如圖12所示。為了方便捕捉黏彈性流體流動過程，采用塑料管作為活塞，玻璃管作為緩沖器腔壁，并按照緩沖器活塞和阻尼孔實際尺寸設計試驗裝置，其中腔體內壁直徑20 mm，活塞厚度15 mm，阻尼孔直徑2.5 mm. 為方便拍攝阻尼孔的流動過程，取消活塞與腔壁之間的間隙，并在塑料管活塞頭部安裝密封圈，防止運動過程中有黏彈性流體從縫隙中擠入。

圖12 啟動流試驗測試系統Fig.12 Test system for the start-up flow

采用日本Mitsubishi公司生產的交流伺服電機作為速度源推動活塞運動，通過電機控制器調整轉速，使活塞桿始終保持恒定速度或加速度。利用美國Fastec公司生產的HiSpec 5 高速相機和配套的圖像采集平臺捕捉活塞啟動時黏彈性流體在阻尼孔中的流動情況。其中，拍攝過程中采用射燈作為輔助光源以提高相機識別度。

在玻璃管緩沖器腔體中注入主體成分為甲基硅氧烷的黏彈性膠體材料，其動力黏度為5×105mm2/s，密度為1.1 kg/m3，調節電機控制器使活塞桿速度從0 mm/s啟動，加速度始終保持1.0 mm/s2，在常溫下進行試驗。高速相機采集頻率為500幀，記錄下不同時刻黏彈性流體流動情況，如圖13所示，圖13中黑色部分為黏彈性膠體材料，白色部分為塑料管活塞。試驗過程中輔助光源在玻璃管緩沖器左側，高速相機拍攝時左半部分容易出現反光點，導致試驗結果左側有部分區域出現缺失，變為白色。由于加速度較小，前0.6 s流體速度也較小，膠體變化不明顯。

圖13 高速相機記錄下阻尼孔啟動流運動情況Fig.13 Test results of start-up flow in orifice

采用表1所列模型參數，代入分數Maxwell模型中得到阻尼孔啟動流流速隨時間的變化結果。對速度進行積分，獲得分數Maxwell模型流體的液面高度分布曲線，如圖14所示。對比觀察試驗和模型結果可知：阻尼孔中心位置流速最大，而靠近阻尼孔邊緣位置流速較小；隨著時間和活塞速度的增大，試驗和數值結果的變化趨勢相同，液面形狀均從扁平狀逐漸發展成為拋物線狀，且阻尼孔啟動流過程始終保持著較強的黏彈性。

表1 分數Maxwell模型參數

圖14 分數Maxwell模型阻尼孔流體液面高度分布Fig.14 Fluid level distribution of fractional Maxwell model in orifice

其中，松弛模量Gt、松弛時間λ和分數階指數α、β通過動態模量理論公式與試驗數據擬合而得(見文獻[20])；而準態特性參數κ和ψ[21]通過上述參數近似獲得：

(31)

以t=1.0 s時黏彈性膠體液面分布試驗結果為例，提出最大液面高度H和最大曲率距離d兩個考量參數，如圖15所示。最大液面高度為黏彈性流體啟動時初始位置到當前時刻液面最高點的垂直距離，是檢驗同一時刻模型結果與試驗結果的一致性。最大曲率距離是從液面最高點到液面最低點做輔助線，平行于輔助線的最大液面切線到輔助線之間的距離，主要用來考察流體的黏彈性。根據實際標尺，測量最大液面高度和最大曲率距離，結果如表2所示，其變化趨勢如圖16和圖17所示。

圖15 黏彈性流體液面測量方法Fig.15 Viscoelastic fluid level measurement method

t/sH/mmd/mm試驗結果分數Maxwell模型結果牛頓流體模型結果試驗結果分數Maxwell模型結果牛頓流體模型結果0.60.37780.36000.36000.16090.15600.07820.80.66120.63460.63460.21020.21780.10961.00.93300.97780.97780.26370.26410.13321.21.39911.39421.39420.30880.29640.14651.41.81511.88041.88040.32500.30980.1549

圖16 最大液面高度結果對比Fig.16 Comparison of measured and simulated maximum fluid levels

圖17 最大曲率距離結果對比Fig.17 Comparison of maximum curvature distances

結合圖16、圖17及表2可知：牛頓流體模型隨著時間的增大，最大曲率距離有趨向于穩態的趨勢，模型的黏彈性很弱；而分數Maxwell模型結果與試驗結果吻合較好，其中最大液面高度最大誤差為4.79%，最大曲率距離的最大誤差4.67%，如圖18所示。由此可見，分數Maxwell模型可以很好地描述黏彈性膠體的黏彈性特點，準確地模擬黏彈性膠體在阻尼孔中啟動流的實際流動過程。

圖18 分數Maxwell模型的誤差Fig.18 Errors of fractional Maxwell model

4 結論

本文根據黏彈性膠體在阻尼緩沖器阻尼孔和節流間隙中的實際流動情況，將其簡化為黏彈性流體在圓管內加速流動和兩平行板間一板不動、一板加速運動兩種單向加速啟動流形式，根據實際工作環境設定初始和邊界條件，提出并構建了含準態特性參數用于研究黏彈性膠體阻尼緩沖器啟動流的分數Maxwell模型，采用有限差分方法求得了模型的數值解。通過與牛頓流體模型進行對比，分析了流速和應力的變化情況及各參數對流速的影響。最后，設計了試驗系統，捕捉了啟動流時黏彈性膠體材料在阻尼孔中的流動過程。通過定義、測量和對比分數Maxwell模型模擬結果與試驗結果的最大液面高度和最大曲率距離驗證了模型的準確性。得出以下結論：

1)黏彈性流體在圓管內和兩平行板間加速流動時，分數Maxwell模型流速分布始終保持著良好的黏彈性。

2)分數Maxwell模型與牛頓流體模型應力呈現線性增長趨勢。

3)隨著參數α、β和ζ的增加，分數Maxwell模型流速分布曲線曲率逐漸下降，其中ζ對曲線曲率的影響比α和β大。

4)啟動流試驗測試系統捕捉到的黏彈性膠體在阻尼孔中啟動流形狀從扁平狀逐漸發展成為拋物線狀，且中心點處速度最大，流動過程具有較強的黏彈性。該結果與分數Maxwell模型模擬結果變化趨勢相同。

5)分數Maxwell模型模擬結果的最大液面高度和最大曲率距離與試驗結果較為吻合，其最大誤差分別為4.79%和4.67%，證明了模型的準確性。

因此，分數Maxwell模型可準確地模擬黏彈性膠體在阻尼緩沖器阻尼孔和節流間隙中啟動流的流動情況，可用于指導黏彈性膠體阻尼緩沖器的設計。