趙金虎

（阜陽師范大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，安徽 阜陽 236037）

0 引言

非牛頓流體廣泛存在于日常生活和自然界中，比如油漆、牛奶、泥漿、高分子溶液、聚合體、生物體的血液、細(xì)胞液等［1］.隨著現(xiàn)代科技的發(fā)展，非牛頓流體力學(xué)的研究理論在許多工業(yè)生產(chǎn)和應(yīng)用科學(xué)領(lǐng)域都有重要應(yīng)用，并對工業(yè)的發(fā)展具有重大的現(xiàn)實(shí)意義［2］.粘彈性流體是非牛頓流體的一種，在剪切流場中，這種流體在外力作用下發(fā)生形變或流動，當(dāng)外力消除后，它的形變會隨時間恢復(fù)或部分恢復(fù)［3］.近年來，分?jǐn)?shù)階模型憑借靈活多變的非局部性質(zhì)在非牛頓流體的粘彈性動力學(xué)表征方面獲得廣泛的應(yīng)用［4-7］.Feng等［8］提出平行平板中非牛頓流體流動的隱式有限差分格式，驗(yàn)證分?jǐn)?shù)階線性方程組的穩(wěn)定性和收斂性.Yang等［9］研究矩形微通道中分?jǐn)?shù)階Maxwell流體的電滲流動，使用譜方法有效解決該模型的控制方程組.Zhang 等［10］考慮矩形管中分?jǐn)?shù)階Maxwell 磁流體在變壓強(qiáng)梯度下的二維流動，求得分析解和數(shù)值解，并作對比.Bai等［11］分析雙向伸長板上的分?jǐn)?shù)階Maxwell流體的流動與傳熱傳質(zhì)，采用有限差分方法求解非線性方程組，收斂精度與數(shù)值特例進(jìn)行驗(yàn)證.Rasheed等［12］討論分?jǐn)?shù)階粘彈性流體流動過程中的化學(xué)反應(yīng)作用，從數(shù)值結(jié)果分析反應(yīng)物的存在對流動的顯著影響.Ahmed［13］研究多孔介質(zhì)方腔中納米流體的對流作用，引用分?jǐn)?shù)階Darcy模型，模擬納米混合物的波形流動.更多有關(guān)分?jǐn)?shù)階模型的應(yīng)用，詳見文獻(xiàn)［14-19］.

在關(guān)于分?jǐn)?shù)階粘彈性流體流動與傳熱的數(shù)值理論研究中，很少有文獻(xiàn)使用有限體積方法求解非線性的分?jǐn)?shù)階控制方程.鑒于此，本文擬結(jié)合有限體積方法與分?jǐn)?shù)階L1格式，用于求解分?jǐn)?shù)階粘彈性流體的自然對流與傳熱傳質(zhì)問題.分?jǐn)?shù)階雙參數(shù)Maxwell模型與流動的動量方程相結(jié)合，建立強(qiáng)耦合非線性且含有多項(xiàng)時間分?jǐn)?shù)階導(dǎo)數(shù)的邊界層控制方程組，最后獲得數(shù)值解并進(jìn)行分析.

1 控制方程

考慮豎直平板上粘彈性流體的自然對流與傳熱傳質(zhì)（受到熱和濃度的浮升力作用），溫度方程加入粘性耗散和放射吸熱，濃度方程也包含熱泳遷移和化學(xué)反應(yīng).除密度變化符合Boussinesq′s假設(shè)，流體的其他性質(zhì)均設(shè)定為常數(shù).平板分別控制在常數(shù)溫度Tw和常數(shù)濃度Cw，且遠(yuǎn)離平板足夠遠(yuǎn)處的溫度T∞和濃度C∞非常小.在邊界層流動中，粘彈性流體的本構(gòu)關(guān)系采用分?jǐn)?shù)階雙參數(shù)Maxwell模型［20］：

其中τxy是剪切力分量，ε是剪切應(yīng)變，ε?=dε/dt是剪切速率，G是剪切模量，λ=μ/G是松弛時間，α和β分別是剪切力和剪切應(yīng)變的分?jǐn)?shù)階導(dǎo)數(shù)參數(shù)，?α/?tα是Caputo型分?jǐn)?shù)階導(dǎo)數(shù)算子［21］：

其中Γ(·)是Gamma函數(shù).

強(qiáng)耦合的動量、溫度和濃度的無量綱控制方程分別為

其中：u和v分別是x軸方向（沿著豎直平板）和y軸方向（垂直于平板）的速度分量，T是溫度，C是濃度，Gr是Grashof 數(shù)，Gm是濃度的Grashof 數(shù)，Nr是浮升力比值，Pr是Prandtl 數(shù)，αf是熱擴(kuò)散系數(shù)，Ec是Eckert數(shù)，Ql是放射吸收系數(shù)，Sc是Schmidt數(shù)，ψ是熱泳系數(shù)，γ是化學(xué)反應(yīng)參數(shù).本文由于篇幅有限，只建立分?jǐn)?shù)階的控制方程組，各無量綱參數(shù)的物理意義參見文獻(xiàn)［22］.

無量綱的初始條件和邊界條件為：

2 數(shù)值方法

本部分介紹求解分?jǐn)?shù)階控制方程組的有限體積方法.為獲得單個控制單元上控制方程的數(shù)值積分，把計(jì)算區(qū)域劃分為離散的控制系統(tǒng)，如圖1所示，其中P為中心節(jié)點(diǎn).然后，在單個控制體積以及有限時間步長Δt下對方程組進(jìn)行積分.最后可以得到關(guān)于每一個節(jié)點(diǎn)P的離散差分方程.

圖1 網(wǎng)格系統(tǒng)的控制單元

從方程（3）開始，首項(xiàng)關(guān)于時間的導(dǎo)數(shù)通過向后差分格式進(jìn)行離散，得到：

其中ΔV=Δx·Δy是控制單元的體積，Δx和Δy分別是空間步長.

對流項(xiàng)的體積積分通過高斯公式轉(zhuǎn)化為面積分，并且采用一階迎風(fēng)差分格式進(jìn)行近似：

其中A表示控制體積的表面積，Aw=Ae=Δx，An=As=Δy.

右端的擴(kuò)散項(xiàng)通過中心差分格式進(jìn)行離散：

粘性耗散的二次項(xiàng)經(jīng)過線性化處理后離散為：

單項(xiàng)濃度采用半隱格式進(jìn)行積分：

方程（2）和（4）中的整數(shù)階導(dǎo)數(shù)采用相似的過程進(jìn)行處理.特別地，方程（4）右端的熱泳作用項(xiàng)的積分為：

對于ω（0＜ω ＜1）階的時間分?jǐn)?shù)階導(dǎo)數(shù)，本文引用L1近似格式進(jìn)行離散［23］：

因此，方程（2）中的多項(xiàng)時間分?jǐn)?shù)階導(dǎo)數(shù)的積分結(jié)果為：

最終獲得內(nèi)節(jié)點(diǎn)P處的迭代差分方程：

其中φ可以代表速度、溫度或濃度.例如，關(guān)于速度差分方程的各項(xiàng)系數(shù)為：

為獲得y軸方向速度分量的數(shù)值解，對連續(xù)性方程（1）進(jìn)行積分得到：

最后，獲得內(nèi)節(jié)點(diǎn)P的迭代差分方程組.迭代截斷誤差控制在10-4以內(nèi)，計(jì)算區(qū)域的2個邊界長分別為Xmax=1 和Ymax=20，其中Ymax對應(yīng)于y→∞且位于速度邊界層之外.為平衡數(shù)值解的精度和計(jì)算耗時，網(wǎng)格數(shù)選定為M=20，N=200，時間步長控制為Δt=0.1.圖2給出當(dāng)空間網(wǎng)格數(shù)分別加倍時，同一組參數(shù)下的速度分布吻合良好.結(jié)果表明，本文選取的網(wǎng)格步長適合計(jì)算.特別指出，本文提出的有限體積方法數(shù)值格式是條件穩(wěn)定的，當(dāng)網(wǎng)格步長的比值超過一定范圍時，數(shù)值解不穩(wěn)定.

圖2 網(wǎng)格無關(guān)性驗(yàn)證

3 結(jié)果討論

圖3和圖4給出分?jǐn)?shù)階導(dǎo)數(shù)參數(shù)α對數(shù)值結(jié)果的影響.圖3表明，隨著α的增加，垂直速度u增大，且邊界層厚度也增加.α=0 表示分?jǐn)?shù)階二階流體，具有最小的速度分布.這個結(jié)果表明，流體采用分?jǐn)?shù)階模型時，隨著粘彈性的增強(qiáng)，表現(xiàn)出剪切增稠的特性.圖4表明，溫度分布隨著α的增加而略微減少，這與圖3中α對速度的作用剛好相反.

圖3 不同α 下的速度分布

圖4 不同α 下的溫度分布

圖5~7給出分?jǐn)?shù)階導(dǎo)數(shù)參數(shù)β對數(shù)值結(jié)果的影響.在圖5中，隨著β的增加，速度分布呈現(xiàn)明顯減少的趨勢，且邊界層的厚度減少.這個結(jié)果與圖3中α對速度的作用截然相反，這是因?yàn)榉謹(jǐn)?shù)階模型中剪切力和剪切速率的分?jǐn)?shù)階參數(shù)異號，當(dāng)分?jǐn)?shù)階參數(shù)小于0時表示的是分?jǐn)?shù)階積分.當(dāng)β=1時表示單參數(shù)的Maxwell模型，該模型具有最小的速度分布和邊界層厚度.圖6和圖7分別給出β對溫度和濃度分布的影響.隨著β的增加，溫度和濃度均顯著增加，且傳熱和傳質(zhì)邊界層厚度也增大.

圖5 不同β 下的速度分布

圖6 不同β 下的溫度分布

圖7 不同β 下的濃度分布

4 小結(jié)

本文結(jié)合有限體積方法和分?jǐn)?shù)階L1格式，求解分?jǐn)?shù)階雙參數(shù)Maxwell 模型的邊界層控制方程組，結(jié)果表明，剪切力和剪切應(yīng)變的分?jǐn)?shù)階導(dǎo)數(shù)參數(shù)對數(shù)值解的作用相反.對于本文建立的強(qiáng)耦合非線性的分?jǐn)?shù)階控制方程組，目前研究還很難求出解析解，所以無法給出穩(wěn)定性和收斂性的證明，這些問題將在下一步的工作中解決.