許曉勤, 陳淑梅, 林紹義

(1. 福建船政交通職業(yè)學(xué)院汽車(chē)運(yùn)用工程系， 福建 福州 350007； 2. 福州大學(xué)機(jī)械工程及自動(dòng)化學(xué)院， 福建 福州 350108)

0 引言

非牛頓流體在工程和工業(yè)中有著廣泛的應(yīng)用， 例如聚合物溶液、 鉆井泥漿、 油漆、 番茄醬、 洗發(fā)水等流體均不服從牛頓粘滯定律， 其應(yīng)力和應(yīng)變率之間的關(guān)系是非線性的. 不同場(chǎng)合的非牛頓流體控制方程更加復(fù)雜且不盡相同， 近年來(lái)比較流行的非牛頓流體模型主要有Maxwell流體、 Williamson流體、 粘彈性流體、 切線雙曲流體等. 其中Casson流體是最常見(jiàn)的非牛頓流體流變模型之一， 能成功描述血液、 果汁、 湯、 醬油、 巧克力等非線性流變行為. 如Sankar等[1]在分析通過(guò)狹窄動(dòng)脈的血流脈動(dòng)情況時(shí)就將血液看作Casson流體， 而將周?chē)难獫{看作牛頓流體處理； Cikrikci等[2]在研究低熱量巧克力配方的物理特性時(shí)也認(rèn)為Casson模型最能擬合不同巧克力配方的流變數(shù)據(jù).

由拉伸板運(yùn)動(dòng)產(chǎn)生的非牛頓流體流動(dòng)在實(shí)際生產(chǎn)中也廣泛存在， 如通過(guò)狹縫擠出熔融聚合物、 塑料片生產(chǎn)、 食品加工、 造紙、 電線和纖維涂層等. 在這些過(guò)程中， 產(chǎn)品的最終質(zhì)量在很大程度上取決于傳熱過(guò)程中的冷卻速率， 而磁流體動(dòng)力學(xué)(magnetohydrodynamic， MHD)參數(shù)是控制冷卻速率， 獲得所需產(chǎn)品質(zhì)量的重要參數(shù)之一[3]. 這也是廣大學(xué)者研究的熱點(diǎn)問(wèn)題之一， 如Nadeem等[4]研究在熱輻射存在下， 指數(shù)收縮板引起的Casson流體MHD流動(dòng)與傳熱； Mukhopadyay等[5]研究在給定熱通量條件下， 指數(shù)拉伸板引起的Casson流體穩(wěn)態(tài)MHD邊界層流動(dòng)與熱傳遞. 結(jié)果均表明磁場(chǎng)能夠抑制速度場(chǎng)而增大溫度場(chǎng).

當(dāng)兩個(gè)接觸物體或同一物體溫度不同時(shí)， 就會(huì)產(chǎn)生熱交換， 過(guò)去主要采用傅立葉定律來(lái)描述. 傅立葉熱傳導(dǎo)定律給出的是拋物線方程， 意味著任何初始溫度變化均會(huì)瞬間被整個(gè)物體感受到， 而這與實(shí)際并不相符[6]. 為了解決這個(gè)問(wèn)題， Cattaneo[7]在傅立葉熱傳導(dǎo)定律基礎(chǔ)上引入熱松弛時(shí)間， 認(rèn)為熱是以低速熱波的形式傳遞. 隨后， Christov[8]利用熱松弛時(shí)間及Oldroyd上對(duì)流導(dǎo)數(shù)修正了Cattaneo定律(即后來(lái)學(xué)者們經(jīng)常提到的Cattaneo-Christov熱通量模型)， 以保持模型的材料不變性. 近年來(lái)， Cattaneo-Christov熱通量模型在各種物理場(chǎng)合的應(yīng)用研究引起了越來(lái)越多學(xué)者的關(guān)注. Hayat等[9]研究了在Cattaneo-Christov熱通量存在時(shí)， 雙擴(kuò)散效應(yīng)對(duì)3D納米流體流動(dòng)的影響； Mustafa[10]探討拉伸板引起的粘彈性流體旋轉(zhuǎn)流動(dòng)與傳熱， 采用Cattaneo-Christov熱通量模型； Mahmood等[11]分析Casson納米流體邊界層流動(dòng)的熵產(chǎn)與傳熱現(xiàn)象， 熱傳導(dǎo)模型采用Cattaneo-Christov熱通量形式.

受上述文獻(xiàn)啟發(fā)， 本文旨在研究由線性拉伸板引起的Casson流體作磁流體流動(dòng)時(shí)的邊界層流動(dòng)與傳熱現(xiàn)象， 熱傳導(dǎo)方程采用Cattaneo-Christov熱通量模型, 能量方程采用Cattaneo-Christov熱通量模型. 利用打靶法， 結(jié)合龍格庫(kù)塔法和牛頓迭代法數(shù)值計(jì)算, 用圖形詳細(xì)分析各控制參數(shù)對(duì)速度、 溫度及表面摩擦系數(shù)、 局部Nusselt數(shù)的影響. 并與經(jīng)典的傅立葉熱傳導(dǎo)定律進(jìn)行對(duì)比分析. 最后文章還給出了Casson參數(shù)和磁場(chǎng)參數(shù)對(duì)流體流動(dòng)和傳熱的影響.

1 數(shù)學(xué)模型

圖1 Casson流體流動(dòng)示意圖 Fig.1 Schematic of Casson fluid flow

研究不可壓縮Casson流體在拉伸表面上的磁流體流動(dòng)， 考慮了Cattaneo-Christov熱通量形式， 施加的磁場(chǎng)為恒定磁場(chǎng)B0. 采用笛卡爾坐標(biāo)系，x軸平行于拉伸表面,y軸垂直于拉伸表面， 如圖1所示. Casson流體作各向同性、 不可壓縮流動(dòng)時(shí)的流變狀態(tài)方程為[12]：

(1)

式中：μB為非牛頓流體的塑性動(dòng)力學(xué)粘度；Py為流體的屈服應(yīng)力；eij為(i，j)變形率分量；π=eijeij為變形率分量的自乘量；πc為基于非牛頓模型變形率分量自乘量的臨界值.

由質(zhì)量守恒定律、 動(dòng)量守恒定理及能量守恒定律可得如下邊界層方程：

(2)

(3)

(4)

邊界條件：

u=uw(x)=cx,v=0,T=Tw(y=0)

(5)

u→0,T→T∞(y→∞)

(6)

引入下列相似變量：

(7)

將式(7)代入方程(2)～(4)， 可知方程(2)自動(dòng)滿(mǎn)足， 方程(3)～(4)簡(jiǎn)化為：

(8)

θ″+Prfθ′-λ0Pr(f2θ″+ff′θ′)=0

(9)

邊界條件(5)～(6)轉(zhuǎn)化為：

f=0,f′=1,θ=1 (η=0)

(10)

f′→0,θ→0 (η→∞)

(11)

(12)

此處沿拉伸板表面摩擦力τw、 熱通量qw分別定義為：

(13)

將式(13)代入式(12)得：

(14)

其中:Rex=uwx/υ為局部Reynolds數(shù). 因此(1+1/β)f″(0)可代表表面摩擦系數(shù)， 而θ′(0)可代表局部Nusselt數(shù).

2 數(shù)值方法

圖2 程序框圖Fig.2 Block diagram

上述非線性微分方程的求解采用龍格庫(kù)塔法， 在Matlab程序包中結(jié)合ode45算法求解. 首先通過(guò)打靶法將邊值問(wèn)題轉(zhuǎn)化成初值問(wèn)題， 假設(shè)初始值的未知參數(shù)， 通過(guò)參數(shù)賦值后采用龍格庫(kù)塔法求解， 通過(guò)牛頓迭代法不斷地調(diào)整未知參數(shù)的值， 直到方程的解滿(mǎn)足精度要求.

龍格-庫(kù)塔法是求解常微分初值問(wèn)題的一種重要方法， 其特點(diǎn)是精度高、 收斂性及穩(wěn)定性好、 計(jì)算過(guò)程可改變步長(zhǎng)等， 是目前求解初值問(wèn)題的一種有效方法. 為進(jìn)一步提高計(jì)算效率和計(jì)算精度， 同時(shí)采用打靶法， 其實(shí)質(zhì)是將邊值問(wèn)題轉(zhuǎn)化成初值問(wèn)題來(lái)求解， 基本思想是將邊值問(wèn)題中未知的初始值設(shè)為未知參數(shù)， 從而將此問(wèn)題看作是帶有未知參數(shù)的初值問(wèn)題. 文章設(shè)置誤差極限為10-5， 邊界層厚度η∞取有限值10， 整個(gè)過(guò)程運(yùn)用Matlab編寫(xiě)程序代碼進(jìn)行迭代求解， 程序框圖如圖2所示. 為驗(yàn)證該數(shù)值方法的準(zhǔn)確性， 將表面摩擦力-f″(0)和局部Nusselt數(shù)-θ′(0)與現(xiàn)有文獻(xiàn)結(jié)果進(jìn)行比較， 結(jié)果很吻合， 如表1～2所示.

表1 M取不同值時(shí)的-f ″(0)

表2 Pr取不同值時(shí)的-θ′(0)

3 結(jié)果討論

圖3給出Casson參數(shù)β對(duì)速度和溫度分布的影響， 由圖可知， 速度輪廓隨著Casson參數(shù)的增大而下降， 而溫度變化趨勢(shì)卻相反. 這是因?yàn)棣略酱螅?塑性動(dòng)力學(xué)粘度越大， 阻礙流體流動(dòng)， 同時(shí)粘度增加使傳熱加快， 從而使溫度上升.

圖3 Casson參數(shù)β對(duì)速度輪廓和溫度分布的影響Fig.3 Effects of Casson parameters β on the velocity profiles and temperature distributions

圖4給出Hartmann數(shù)M對(duì)速度和溫度分布的影響， 由圖4(a)可知， 速度和邊界層厚度隨著M的增大而減小， 這是因?yàn)榇艌?chǎng)的上升會(huì)產(chǎn)生與流動(dòng)方向相反的力(洛倫磁力)， 阻礙流體流動(dòng)， 使速度分布減小. 由圖4(b)可知， 溫度的變化趨勢(shì)與速度相反， 這是因?yàn)槁鍌惔帕?huì)產(chǎn)生焦耳熱， 使溫度增加. 因此磁場(chǎng)可以控制冷卻速率， 改善產(chǎn)品最終質(zhì)量.

圖4 Hartmann數(shù)M對(duì)速度輪廓和溫度分布的影響Fig.4 Effects of Hartmann number M on the velocity profiles and temperature distributions

圖5給出Prandtl數(shù)Pr和熱松弛參數(shù)λ0對(duì)溫度分布的影響， 由圖5(a)可知， 溫度分布及其邊界層厚度隨Prandtl數(shù)的增大而下降， 這是因?yàn)镻randtl數(shù)越大， 熱擴(kuò)散系數(shù)越小， 熱擴(kuò)散速度越慢， 限制溫度的升高. 由圖5(b)可知， 隨著熱松弛參數(shù)λ0的增大， 溫度也下降且邊界層變薄， 從物理層面解釋?zhuān)?熱松弛參數(shù)越大， 材料顆粒將熱量傳遞到相鄰顆粒所需的時(shí)間越長(zhǎng)， 從而降低了溫度. 值得一提的是， 當(dāng)熱松弛參數(shù)λ0=0時(shí)， 熱通量表達(dá)式與經(jīng)典傅立葉熱傳導(dǎo)定律相同. 因此與采用經(jīng)典的傅立葉熱傳導(dǎo)定律相比， 采用Cattaneo-Christov熱通量模型得到的溫度和邊界層厚度均比較低.

圖5 Prandtl數(shù)Pr和熱松弛參數(shù)λ0對(duì)溫度分布的影響Fig.5 Effects of Prandtl number Pr and thermal relaxation parameter λ0 on the temperature distributions

圖6給出了不考慮熱松弛情況下Casson參數(shù)β和Hartmann數(shù)M對(duì)表面摩擦系數(shù)和局部Nusselt數(shù)的影響. 由圖6可知， 當(dāng)Hartmann數(shù)M不變時(shí)， 表面摩擦系數(shù)和局部Nusselt數(shù)均隨著Casson參數(shù)β的增大而減小, 這是因?yàn)殡S著β值的增加， 從壁面?zhèn)鬟f給流體的熱量減少， 局部Nusselt數(shù)降低. 而對(duì)于所有的Casson參數(shù)β， 表面摩擦系數(shù)均隨著Hartmann數(shù)M的增大而增大， 局部Nusselt數(shù)均隨著Hartmann數(shù)M的增大而減小， 這可能是由于洛倫磁力的存在增大了壁面摩擦力， 阻礙流體流動(dòng)， 使流體對(duì)流傳熱程度減弱.

圖6 Casson參數(shù)β和Hartmann數(shù)M對(duì)表面摩擦系數(shù)及局部Nusselt數(shù)的影響Fig.6 Effects of the Casson parameter β and Hartmann number M on the skin-friction coefficient and local Nusselt number

4 結(jié)語(yǔ)

1) 隨著Casson參數(shù)或Hartmann 數(shù)的增大， 速度下降而溫度上升， 通過(guò)控制磁場(chǎng)強(qiáng)度可以控制流體流動(dòng)速度和冷卻速率.

2) 溫度隨Prandtl數(shù)或熱松弛參數(shù)增大而下降， 與采用經(jīng)典的傅立葉熱傳導(dǎo)定律相比， 采用Cattaneo-Christov熱通量模型得到的溫度和邊界層厚度均比較低.

3) 表面摩擦系數(shù)和局部Nusselt數(shù)均隨著Casson參數(shù)β的增大而減小.

4) 表面摩擦系數(shù)隨著Hartmann數(shù)M的增大而增大， 而局部Nusselt數(shù)卻相反.