蘇祥龍 許文祥 陳 文

(河海大學(xué)力學(xué)與材料學(xué)院軟物質(zhì)力學(xué)研究所，南京211100)

流體力學(xué)

基于分形導(dǎo)數(shù)對(duì)非牛頓流體層流的數(shù)值研究1)

蘇祥龍 許文祥 陳 文2)

(河海大學(xué)力學(xué)與材料學(xué)院軟物質(zhì)力學(xué)研究所，南京211100)

非牛頓流體具有復(fù)雜的流變特性，揭示該流變特性可以更加合理地指導(dǎo)非牛頓流體在工農(nóng)業(yè)生產(chǎn)中的應(yīng)用.經(jīng)典的非牛頓流體本構(gòu)模型往往形式復(fù)雜，僅能應(yīng)用于某些特定的情況.分?jǐn)?shù)階導(dǎo)數(shù)模型具有參數(shù)少和形式簡(jiǎn)單的特點(diǎn)，已成功地應(yīng)用于描述非牛頓流體的運(yùn)動(dòng).Hausdor ff分形導(dǎo)數(shù)作為一個(gè)備選的建模方法，相比分?jǐn)?shù)階導(dǎo)數(shù)具有更簡(jiǎn)單的形式以及更高的計(jì)算效率.本文基于Hausdor ff分形導(dǎo)數(shù)改進(jìn)現(xiàn)有牛頓黏性模型，提出分形黏壺模型.通過研究分形黏壺在常應(yīng)變率下表觀黏度的變化情況，以及在加、卸載條件下的蠕變及恢復(fù)特性，發(fā)現(xiàn)分形黏壺模型適合于描述具有黏彈性的非牛頓流體(本文稱之為分形流體).結(jié)合連續(xù)性方程及運(yùn)動(dòng)微分方程，推導(dǎo)出分形流體在平行板間層流的基本方程.按是否拖動(dòng)上板和是否存在水平的壓力梯度分為3種工況，分別用數(shù)值方法計(jì)算這3種工況下流速在板間的分布及其隨時(shí)間變化的情況.通過分析不同工況下的流速分布，發(fā)現(xiàn)水平的壓力梯度會(huì)改變流速隨時(shí)間變化的形狀，且會(huì)推遲流速到達(dá)穩(wěn)定的時(shí)間.在水平壓力梯度不存在的情況下，不同階數(shù)的分形流體具有相同的流速分布或是演變過程.另外，在水平壓力梯度存在的情況下，上板速度不影響不同階數(shù)分形流體間穩(wěn)定速度的差值.

分形導(dǎo)數(shù)，非牛頓流體，分形流體，層流，表觀黏度

引言

非牛頓流體廣泛地存在于食品工業(yè)，醫(yī)學(xué)，石油工業(yè)等領(lǐng)域.不同于牛頓流體，非牛頓流體展現(xiàn)出復(fù)雜的流變現(xiàn)象，包括蠕變、時(shí)間依賴和剪切依賴的黏性.黏性依賴于剪切率的現(xiàn)象在實(shí)際生活中非常常見，比如玉米淀粉的懸浮液，奶油，番茄醬，冰和血液.參照牛頓流體，定義非牛頓流體的切應(yīng)力與應(yīng)變率之比為其表觀黏度[1]，它通常是應(yīng)變率的函數(shù).這種模型被稱為廣義牛頓模型[2]，其表達(dá)式為[3]

式中，T是偏應(yīng)力張量，μ是表觀黏度函數(shù);D是應(yīng)變率張量，D=0.5h▽V+(▽V)Ti,其中，V 是速度向量，▽V是速度梯度，上標(biāo)“T”表示轉(zhuǎn)置；Π是應(yīng)變率張量D的不變量，Π=2tr(D2)，其中符號(hào)“tr”表示張量的跡，即主對(duì)角元素求和.冪律模型作為廣義牛頓模型的一個(gè)特例，其表觀黏度函數(shù)的表達(dá)式為[3]

其中，K為流體的稠度系數(shù)，n是材料參數(shù).當(dāng)n=1時(shí)冪律模型退化為牛頓黏性模型，而n＞1,n＜1分別可以描述剪切增稠和剪切變稀的流體[3].另外，Oldroyd模型和Carreau模型[3]是根據(jù)表觀黏度在應(yīng)變率趨于零和趨于無(wú)窮時(shí)的極限值而寫出的本構(gòu)關(guān)系的表達(dá)式.如果流體存在屈服應(yīng)力，且超過屈服應(yīng)力之后應(yīng)力與應(yīng)變呈線性關(guān)系，可以用Bingham模型[4]來(lái)描述.而Herschel-Bulkley模型[5]和Casson模型[6]可以描述屈服應(yīng)力之上應(yīng)力與應(yīng)變呈非線性關(guān)系的情況.另一方面，時(shí)間依賴的黏性也普遍存在.比如酸奶，水基氧化鐵凝膠和一些泥土，它們的黏性隨時(shí)間減小，稱為觸變性[7].而另外一些流體，像打印機(jī)油墨，滑膜液黏性卻是隨時(shí)間增加，稱為震凝性[8].有關(guān)非牛頓流體的觸變性和震凝性的本構(gòu)模型研究見參考文獻(xiàn)[9-11].

另外，關(guān)于非牛頓流體流動(dòng)特性的研究也很多.范椿[12]研究了對(duì)于二維擾動(dòng)非牛頓冪律流體在重力作用下沿傾斜面流動(dòng)的穩(wěn)定性,得到臨界雷諾數(shù)作為冪律指數(shù) n和平板斜率 β的函數(shù).韓式方等[13]研究了上隨體Maxwell流體在圓管內(nèi)非定常流動(dòng)規(guī)律.董波等[14]基于插值補(bǔ)充格子波爾茲曼方法模擬了冪律流體的圓柱繞流問題.Hayat等[15]推導(dǎo)了Oldroyd-B流體的一些簡(jiǎn)單流動(dòng)的精確解，包括在無(wú)限大平行板間的流動(dòng).近年來(lái)，分?jǐn)?shù)階導(dǎo)數(shù)作為一種優(yōu)秀的數(shù)學(xué)建模工具，成功地應(yīng)用到粘彈性建模領(lǐng)域[1619].另外，分?jǐn)?shù)階本構(gòu)方程也被用于描述非牛頓流體的運(yùn)動(dòng)[2022]，包括求解非牛頓流體流動(dòng)的精確解[2326].然而，在數(shù)值計(jì)算時(shí)分?jǐn)?shù)階微積分算子的全局性特征要求相當(dāng)大的計(jì)算成本和存儲(chǔ)空間[27].

陳文[28]基于時(shí)空變換首次提出了Hausdor ff分形導(dǎo)數(shù)，與分?jǐn)?shù)階導(dǎo)數(shù)模型相比，Hausdor ff分形導(dǎo)數(shù)簡(jiǎn)單的數(shù)學(xué)定義不僅給理論推導(dǎo)帶來(lái)了便利，也為計(jì)算模擬節(jié)約了大量成本.目前Hausdor ff分形導(dǎo)數(shù)已成功地應(yīng)用于反常擴(kuò)散[2830]、松弛振動(dòng)[31]、熱控制[32]等領(lǐng)域.蔡偉等[27]首次將Hausdor ff分形導(dǎo)數(shù)應(yīng)用到黏彈性建模中，推導(dǎo)了分形Maxwell和分形Kelvin模型的蠕變?nèi)崃亢退沙谀Ａ浚瑪?shù)值結(jié)果表明這些模型可以很好地描述黏彈性行為.蘇祥龍等[33]將Hausdor ff分形導(dǎo)數(shù)用于描述非牛頓流體的流變行為，研究表明分形導(dǎo)數(shù)模型適宜描述非牛頓流體的剪切依賴和時(shí)間依賴的表觀黏度以及蠕變等流變現(xiàn)象.

本文用Hausdor ff分形導(dǎo)數(shù)替換牛頓黏性模型里的常規(guī)導(dǎo)數(shù)，得到一種描述非牛頓流體的分形黏壺模型.研究了分形黏壺在常應(yīng)變率下的應(yīng)力響應(yīng)以及蠕變和回復(fù)特性.并在此基礎(chǔ)上推導(dǎo)了分形流體在平行板層流情況下的控制方程.用數(shù)值方法計(jì)算了多種工況下板間流速分布以及演變過程.并探究了上板的速度、水平壓力梯度以及分形導(dǎo)數(shù)的階數(shù)對(duì)流速分布及演變的影響.

1 Hausdor ff分形導(dǎo)數(shù)及其性質(zhì)

基于時(shí)空尺度變換，陳文[28]提出了一種新的導(dǎo)數(shù)——Hausdor ff分形導(dǎo)數(shù).函數(shù) f(t)在時(shí)間方向上的Hausdor ff分形導(dǎo)數(shù)定義為

式中，d/dtα是時(shí)間分形導(dǎo)數(shù)算子，α是分形導(dǎo)數(shù)的階數(shù)，它代表一種分形量度[28]；t和t′分別代表不同的時(shí)刻.在分形不變量及分形平衡的假設(shè)下，運(yùn)用變量替換=tα，可將分形導(dǎo)數(shù)轉(zhuǎn)化為常規(guī)導(dǎo)數(shù)[28].

若函數(shù) f(t)在定義域內(nèi)一階可導(dǎo)，且t1?α/α不為零，那么Hausdor ff分形導(dǎo)數(shù)與常規(guī)導(dǎo)數(shù)有如下關(guān)系[27]

設(shè)函數(shù) f(t)和g(t)在定義域內(nèi)一階可導(dǎo)，g(t)≠0，且各自的α階分形導(dǎo)數(shù)存在，在這種情況下我們可以推導(dǎo)出Hausdor ff分形導(dǎo)數(shù)的四則運(yùn)算法則：

(1)函數(shù)和差的分形導(dǎo)數(shù)

(2)函數(shù)乘積的分形導(dǎo)數(shù)

(3)函數(shù)商的分形導(dǎo)數(shù)

可以看到，Hausdor ff分形導(dǎo)數(shù)的四則運(yùn)算法則與常規(guī)導(dǎo)數(shù)的四則運(yùn)算法則一致.

2 分形黏壺模型

本文用Hausdor ff分形導(dǎo)數(shù)取代牛頓黏性模型里的常規(guī)導(dǎo)數(shù),得到一種新的本構(gòu)模型，稱其為分形黏壺模型，其本構(gòu)關(guān)系為

其中，τ和ε分別表示剪切應(yīng)力與應(yīng)變，η是材料參數(shù).α為無(wú)量綱數(shù)，文獻(xiàn)[27]給出了基于分形黏壺與彈簧的串并聯(lián)模型得到的分形Maxwell模型，分形Kelvin模型和分形Zener模型，這些模型成功地描述了軟土、凍土、鹽巖以及混凝土等材料的蠕變現(xiàn)象.

為了更清楚地認(rèn)識(shí)分形黏壺的性質(zhì)，本文研究：(1)分形黏壺的表觀黏度隨時(shí)間的變化；(2)分形黏壺的蠕變及恢復(fù)特性.

2.1 表觀黏度隨時(shí)間變化

令應(yīng)變率為常數(shù)c，即ε(t)=ct，代入式(8)得到應(yīng)力隨時(shí)間的變化關(guān)系為

那么這種情況下分形黏壺的表觀黏度ηb可表示為

令 η =1kg·m?1·sα?2，在 α 分別取 0.3,0.5 和 0.7 時(shí)畫出分形黏壺的表觀黏度隨時(shí)間的變化情況，見圖1.從圖1可以看到，分形黏壺的表觀黏度隨時(shí)間逐漸增大；階數(shù)越大，則表觀黏度越小，且其增長(zhǎng)速率越小.當(dāng)α趨近于1時(shí)，表觀黏度幾乎維持穩(wěn)定.由分形導(dǎo)數(shù)的定義式(3)，令α→1，那么分形導(dǎo)數(shù)即退化為常規(guī)導(dǎo)數(shù)，應(yīng)變函數(shù)ε(t)的α階分形導(dǎo)數(shù)則退化為ε(t)的一階導(dǎo)數(shù)，即應(yīng)變率函數(shù)˙ε(t)

考慮上式的變換，則當(dāng)α→1時(shí)，分形黏壺模型退化為牛頓黏性模型

圖1 分形黏壺在常應(yīng)變率下的表觀黏度Fig.1 The apparent viscosity of the fractal dashpot under constant strain rate

2.2 加、卸載時(shí)的蠕變及恢復(fù)特性

令τ(t)=τ0，代入式(8)中得到分形黏壺的蠕變響應(yīng)為

選擇如圖2所示的加載路徑，在t1=1s時(shí)刻加上應(yīng)力τ0，在t2=2s時(shí)刻將應(yīng)力撤掉，得到分形黏壺蠕變及恢復(fù)曲線如圖3.

圖2 分形黏壺加載及卸載圖Fig.2 The loading and unloading for the fractal dashpot

圖3 分形黏壺在不同階數(shù)下的蠕變及恢復(fù)曲線Fig.3 The creep and recovery curves of the fractal dashpot

從圖3可以看到，分形黏壺具有明顯的黏彈性.其中，當(dāng)α較小時(shí)，分形黏壺蠕變及恢復(fù)變形發(fā)展較快，展現(xiàn)出高彈性.另外，在α比較大的情況下，蠕變變形幾乎與時(shí)間呈線性關(guān)系而恢復(fù)變形發(fā)展很慢，這種性質(zhì)表現(xiàn)跟黏性流體很類似.由分析可知，分形黏壺的黏彈性與階數(shù)α的取值密切相關(guān).

另一方面，從非牛頓流體的角度來(lái)說(shuō)，分形黏壺可以描述表觀黏度隨時(shí)間增大的流體.而且分形黏壺的階數(shù)α取值越大，表觀黏度越小.

3 分形流體平行板層流的基本方程

為了探究分形流體的流動(dòng)特點(diǎn)，設(shè)置一個(gè)簡(jiǎn)單流動(dòng).如圖4所示，兩無(wú)限寬的平行剛性板水平放置在分形流體中，已知平行板間距為a，分形流體的密度為ρ，下板固定.

圖4 兩平行板間分形流體的層流運(yùn)動(dòng)Fig.4 The laminar fl ow of the fractal fl uid between two parallel plates

將滿足式(8)的具有黏彈性的非牛頓流體稱為分形流體，考慮分形流體是不可壓縮、各向同性的，那么分形流體的本構(gòu)關(guān)系可以寫成如下形式

其中，σ是柯西應(yīng)力張量，T是偏應(yīng)力張量，D是應(yīng)變率張量，I是單位張量，p表示壓力，μ是表觀黏度函數(shù)，具體寫作

考慮分形流體是不可壓縮的，其連續(xù)性方程為

其中f是單位質(zhì)量力矢量，ρ是密度.

對(duì)于分形流體在如圖4所示的平行板間的層流運(yùn)動(dòng)，速度矢量、應(yīng)力張量只在x方向上的分量不為零，且是z和t的函數(shù)；單位質(zhì)量力矢量?jī)H有z方向上的重力項(xiàng)，記作

其中，u表示速度矢量在x方向上的分量；σ表示應(yīng)力張量在x方向上的分量；g表示重力加速度，負(fù)號(hào)表示與圖4中指示的坐標(biāo)軸正方向相反.根據(jù)上式，結(jié)合式(12)～式(14)，簡(jiǎn)化分形流體的運(yùn)動(dòng)微分方程式(15)得

對(duì)式(18)積分得

其中，f(x)是關(guān)于x的任意函數(shù).可知?p/?x僅是x的函數(shù)，而u不是x的函數(shù)，所以求解u時(shí)?p/?x可做常數(shù)處理.將式(13)代入式(17)中，即得到分形流體平行板層流的基本方程

這是一個(gè)變系數(shù)的二階線性非齊次偏微分方程，本文采用軟件MATLAB中的pdetool工具箱求解.另外，牛頓流體在平行板間層流的穩(wěn)定精確解為[34]

4 計(jì)算結(jié)果及討論

根據(jù)上板是否存在初始速度以及是否存在沿x方向的壓力梯度，分別在3種工況下研究分形流體在平行板間層流的流速分布及變化規(guī)律.

采用pdetool計(jì)算時(shí)，為計(jì)算方便，取a=1m,ρ =1kg/m3,η =1kg·m?1·sα?1，并取計(jì)算時(shí)間為 3s，時(shí)間和空間上分別用30個(gè)網(wǎng)格點(diǎn)劃分.

對(duì)比式(13)和式(10)可知，平行板層流情況下分形流體的表觀黏度函數(shù)與常應(yīng)變率下的表觀黏度函數(shù)一致.初始時(shí)刻分形流體的表觀黏度為零，且隨時(shí)間以冪函數(shù)形式增長(zhǎng).對(duì)于不同的階數(shù)α，階數(shù)越大，表觀黏度函數(shù)反而越小.具體變化情況見圖1.

工況1：存在水平壓力梯度且只拖動(dòng)上板

在此工況下，可以給出分形流體流動(dòng)的初邊值條件為

令U=1m/s,?p/?x= ?1Pa/m，當(dāng) α 取 0.3,0.5,0.7時(shí)分別計(jì)算流速.當(dāng)α=0.3時(shí)求解得到的流速在板間的分布以及隨時(shí)間的變化，見圖5.

由圖5可知，從下板(z=0)到上板(z=1m)流速由0變化到1m/s，在板拖動(dòng)的初始時(shí)刻，板間流速分布變化較為劇烈，隨著時(shí)間的推移，流速分布趨于平穩(wěn).

圖5 當(dāng)U=1m/s,?p/?x=?1Pa/m時(shí)分形流體的流速分布及隨時(shí)間變化(α=0.3)Fig.5 The distribution and evolution of velocity of fractal fl uid for U=1m/s,?p/?x= ?1Pa/m(α =0.3)

由圖1或是式(13)可知，分形流體的表觀黏度隨階數(shù)α增大而減小，且隨時(shí)間延長(zhǎng).若分形流體的表觀黏度影響其流速的分布，那么兩板間的流速分布永遠(yuǎn)不可能達(dá)到真正的穩(wěn)定狀態(tài).但表觀黏度的變化率卻是逐漸減小的

從圖5的流速分布圖也可以看出，相比于初始時(shí)刻，流速在t=3s時(shí)刻變化已經(jīng)非常小了.以圖6為例，在t為0～0.5s階段，流速?gòu)?增加到0.56m/s，平均增長(zhǎng)率為1.12m/s2；而在t為2.5～3s階段，流速?gòu)?.5198m/s降到0.5174m/s，平均減少速率為4.8×10?3m/s2.前后兩時(shí)間段流速的平均變化率相差了近200倍！因此，將t=3s時(shí)刻近似穩(wěn)定的流速作為考察量，稱其為“穩(wěn)定流速”.

如圖7給出了當(dāng)α分別取0.3,0.5和0.7時(shí)分形流體板間穩(wěn)定流速的分布以及由式(21)得到的牛頓流體的經(jīng)典的穩(wěn)定流速分布.從圖7可知，流速在板間以類拋物線形式分布，且下板速度為0，上板速度為1m/s.分形流體的階數(shù)α越大，即表觀黏度越小，其穩(wěn)定流速也越大，且越接近于(低于)經(jīng)典流速分布.

圖6給出了z=0.5m(板間中間點(diǎn))處的分形流體的流速隨時(shí)間的變化曲線.由圖可知，流速?gòu)牧阍龃蟮揭粋€(gè)極大值，之后衰減并趨于穩(wěn)定，清晰地再現(xiàn)了分形流體在平行板間流動(dòng)的啟動(dòng)及穩(wěn)定過程.且分形流體的表觀黏度越小，對(duì)應(yīng)的流速值越大.

圖6 當(dāng)U=1m/s,?p/?x=?1Pa/m時(shí)分形流體在z=0.5m處的流速隨時(shí)間變化圖Fig.6 The evolution of velocity of fractal fl uid on z=0.5m for U=1m/s,?p/?x= ?1Pa/m

圖7 當(dāng)U=1m/s,?p/?x=?1Pa/m時(shí)分形流體的穩(wěn)定流速分布(t=3s)Fig.7 The stable velocity distribution(t=3s)of fractal fl uid for U=1m/s,?p/?x= ?1Pa/m

工況2：存在水平壓力梯度且上、下板固定

在此工況下，可以給出分形流體流動(dòng)的初邊值條件為

令?p/?x=?1Pa/m，在α取0.3,0.5,0.7時(shí)分別計(jì)算分形流體的流速.圖8給出了α取0.3時(shí)分形流體的流速分布及變化情況，可以發(fā)現(xiàn)，流速值明顯比工況1下的小.

圖8 當(dāng)U=0,?p/?x=?1Pa/m時(shí)分形流體的流速分布及隨時(shí)間變化圖(α=0.3)Fig.8 The distribution and evolution of velocity of fractal fl uid for U=0,?p/?x= ?1Pa/m(α =0.3)

圖9給出了當(dāng)α分別取0.3,0.5和0.7時(shí)分形流體的穩(wěn)定流速分布以及牛頓流體的經(jīng)典穩(wěn)定流速分布.可以看到，流速在板間是以類拋物線形式分布的，板兩端(z=0和z=1m)速度為零，最大速度在板中間取得.另外也可以看到，表觀黏度越小，分形流體達(dá)到的穩(wěn)定流速越大.

圖9 當(dāng)U=0,?p/?x=?1Pa/m時(shí)分形流體的穩(wěn)定流速分布(t=3s)Fig.9 The stable velocity distribution(t=3s)of fractal fl uid for U=0,?p/?x= ?1Pa/m

另外，圖10給出了z=0.5m(板中間點(diǎn))處的流速隨時(shí)間的變化曲線.從圖中可以看出，流速在一個(gè)短時(shí)間內(nèi)增加到極大值，隨后逐漸衰減并趨于穩(wěn)定.

對(duì)比工況1，雖然兩種情況下流速值相差較大，但是對(duì)于不同階數(shù)的分形流體它們的穩(wěn)定流速的差值沒有太大變化.比如在工況1的板中間點(diǎn)處,α=0.3和α=0.7兩種分形流體穩(wěn)定流速相差大概0.05m/s(圖6或圖7)，工況2下兩種分形流體穩(wěn)定流速同樣相差0.05m/s左右(圖9或圖10).說(shuō)明上板的速度對(duì)不同分形流體間相對(duì)速度大小沒有影響.另外，對(duì)比圖10與圖6，發(fā)現(xiàn)拖動(dòng)上板會(huì)縮短流速達(dá)到穩(wěn)定的時(shí)間.

圖10 當(dāng)U=0,?p/?x=?1Pa/m時(shí)分形流體在z=0.5m處的流速隨時(shí)間變化圖Fig.10 The evolution of velocity of fractal fl uid on z=0.5m for U=0,?p/?x= ?1Pa/m

工況3：不存在水平壓力梯度且只拖動(dòng)上板

在此工況下，可以給出分形流體流動(dòng)的初邊值條件為

令 U=1m/s,?p/?x=0，在 α 取 0.3,0.5,0.7時(shí)分別計(jì)算分形流體的流速.求解得到分形流體的流速在板間的分布以及隨時(shí)間的變化，見圖11.由圖可知，從下板(z=0)到上板(z=1m)流速由0變化到1m/s，與圖5類似.

圖11 當(dāng)U=1m/s,?p/?x=0時(shí)分形流體流速分布及隨時(shí)間變化圖(α=0.3)Fig.11 The distribution and evolution of velocity of fractal fl uid for U=1m/s,?p/?x=0(α =0.3)

圖12給出了當(dāng)α分別取0.3,0.5和0.7時(shí)分形流體的穩(wěn)定流速分布以及牛頓流體的經(jīng)典穩(wěn)定流速分布.可以看到，流速在板間是線性分布的.與上兩種工況不同，這種情況下不同分形流體之間流速分布相同，且與經(jīng)典的流速分布重合.

圖12 當(dāng)U=1m/s,?p/?x=0時(shí)分形流體的穩(wěn)定流速分布(t=3s)Fig.12 The stable velocity distribution(t=3s)of fractal fl uid for U=1m/s,?p/?x=0

由圖13可以看到，板中間位置點(diǎn)的流速在很短的時(shí)間內(nèi)從零增加到極大值，而后直接進(jìn)入穩(wěn)定狀態(tài)，流速不再變化.

圖13 當(dāng)U=1m/s,?p/?x=0時(shí)分形流體在z=0.5m處的流速隨時(shí)間變化圖Fig.13 The evolution of velocity of fractal fl uid on z=0.5m for U=1m/s,?p/?x=0

對(duì)比圖6和圖13可知，水平的壓力梯度會(huì)使流速達(dá)到一個(gè)高峰然后衰減，而無(wú)水平壓力梯度時(shí)分形流體的流速?zèng)]有衰減過程，直接到達(dá)穩(wěn)定.而且水平壓力梯度的存在推遲了分形流體速度穩(wěn)定的時(shí)間.另外，水平壓力梯度存在時(shí)分形流體的階數(shù)α對(duì)流速大小影響顯著，而無(wú)水平壓力梯度時(shí)流速對(duì)階數(shù)α的變化不敏感.這說(shuō)明分形流體的表觀黏度只在水平壓力梯度存在時(shí)對(duì)流速分布有影響.

5 結(jié)論

(1)本文基于Hausdor ff分形導(dǎo)數(shù)提出了分形黏壺模型，探究了分形黏壺的力學(xué)特性.發(fā)現(xiàn)分形黏壺呈現(xiàn)出明顯的黏彈性特征，而且適合描述表觀黏度隨時(shí)間變化的非牛頓流體.結(jié)合運(yùn)動(dòng)方程和連續(xù)性方程，推導(dǎo)了分形流體在平行板層流情況下的基本方程.用數(shù)值方法計(jì)算了3種工況下分形流體在板間的流速分布以及隨時(shí)間變化的情況.

(2)計(jì)算結(jié)果顯示，流速在短時(shí)間內(nèi)增加到極大值，而后隨時(shí)間衰減并趨于穩(wěn)定或直接達(dá)到穩(wěn)定.分形流體層流的穩(wěn)定流速分布與經(jīng)典的牛頓流體流速分布有類似的形狀，但都不大于經(jīng)典流速值.

(3)對(duì)比3種工況發(fā)現(xiàn)，上板的速度會(huì)顯著影響流速的大小以及在板間的速度分布情況.而水平的壓力梯度會(huì)改變流速隨時(shí)間變化的形狀，且會(huì)推遲流速到達(dá)穩(wěn)定的時(shí)間.在不存在水平壓力梯度的情況下，不同階數(shù)的分形流體具有相同的流速分布或是演變.另外，在存在水平壓力梯度的情況下，上板速度不影響不同階數(shù)分形流體間穩(wěn)定速度的差值.

1 Fox RW,Mcdonald AT,Pritchard PJ,et al.Fluid Mechanics.[s.l.]:John Wiley&Sons Inc,2012

2 Galindo-Rosales FJ,Rubio-Hern′andez FJ,Sevilla A.An apparent viscosity function for shear thickening fl uids. Journal of Non-Newtonian Fluid Mechanics,2011,166(5):321-325

3 陳文芳,蔡扶時(shí).非牛頓流體的一些本構(gòu)方程.力學(xué)學(xué)報(bào),1983,19(1):16-26(Chen Wenfang,Cai fushi.Some constitutive equations for non-Newtonian fl uid.Acta Mechanica Sinica,1983,19(1):16-26(in Chinese))

4 Bingham EC.Fluidity and Plasticity.[s.l.]:McGraw-Hill Book Co.,1922

5 Herschel W,Bulkley R.Measurement of consistency as applied to rubber-benzene solutions.Am Soc Test Proc,1926,26(2):621-633

6 Casson N.Rheology of disperse systems//Proceedings of a conference organized by the British Society of Rheology,Pergamon Press,New York,1959

7 Barnes HA.Thixotropy—areview.JournalofNon-NewtonianFluid Mechanics,1997,70(1-2):1-33

8 Oates KM,Krause WE,Jones RL,et al.Rheopexy of synovial fl uid and protein aggregation.Journal of the Royal Society Interface,2006,3(6):167-174

9 Bautista F,De Santos JM,Puig JE,et al.Understanding thixotropic and antithixotropic behavior of viscoelastic micellar solutions and liquid crystalline dispersions. I.The model. Journal of Non-Newtonian Fluid Mechanics,1999,80(2):93-113

10 Mujumdar A,Beris AN,Metzner AB.Transient phenomena in thixotropic systems.Journal of Non-Newtonian Fluid Mechanics,2002,102(2):157-178

11 Blackwell BC,Ewoldt RH.A simple thixotropic–viscoelastic constitutive model produces unique signatures in large-amplitude oscillatory shear(LAOS).Journal of Nonnewtonian Fluid Mechanics,2014,s(208-209):27-41

12 范椿.非牛頓冪律流體沿傾斜面流動(dòng)的穩(wěn)定性.力學(xué)學(xué)報(bào),1982,18(2):155-160(Fan Chun.Stability of non-Newtonian power law fl uid fl owing down an inclined plane.Acta Mechanica Sinica,1982,18(2):155-160(in Chinese))

13 韓式方,伍岳慶.管內(nèi)上隨體Maxwell流體非定常流動(dòng).力學(xué)學(xué)報(bào),1990,22(5):519-525(Han Shifang,Wu Yueqing.A study on non-steady fl ow of upper-convected Maxwell fl uid in tube.Acta Mechanica Sinica,1990,22(5):519-525(in Chinese))

14 董波,李維仲,馮玉靜等.冪律流體圓柱繞流的格子波爾茲曼模擬.力學(xué)學(xué)報(bào),2014,46(1):44-53(Dong Bo,Li Weizhong,Feng Yujing,et al.Lattice boltzmann simulation of a power-law fl uid past a circular cylinder.Acta Mechanica Sinica,2014,46(1):44-53(in Chinese))

15 Hayat T,Siddiqui AM,Asghar S.Some simple fl ows of an Oldroyd-B fl uid.International Journal of Engineering Science,2001,39(2):135-147

16 Fan W,Jiang X,Qi H.Parameter estimation for the generalized fractional element network Zener model based on the Bayesian method.Physica A:Statistical Mechanics&Its Applications,2015,427:40-49

17 Bagley RL,Torvik PJ.Fractional calculus-A di ff erent approach to the analysis of viscoelastically damped structures.AIAA Journal,2012,21(5):741-748

18 Bagley RL,Torvik PJ.On the fractional calculus model of viscoelastic behavior.Journal of Rheology,1986,30(1):133-155

19 Blair GWS.Analytical and integrative aspects of the stress-straintime problem.Journal of Scienti fi c Instruments,2002,21(5):80-84

20 Hayat T,Nadeem S,Asghar S.Periodic unidirectional fl ows of a viscoelastic fl uid with the fractional Maxwell model.Applied Mathematics&Computation,2004,151(1):153-161

21 Tan W,Xu M.Plane surface suddenly set in motion in a viscoelastic fl uid with fractional Maxwell model.Acta Mechanica Sinica,2002,18(4):342-349

22 Tripathi D.Peristaltic transport of fractional Maxwell fl uids in uniform tubes:Applications in endoscopy.Computers&Mathematics with Applications,2011,62(3):1116-1126

23 Mahmood A,Parveen S,Ara A,et al.Exact analytic solutions for the unsteady fl ow of a non-Newtonian fl uid between two cylinders with fractional derivative model.Communications in Nonlinear Science&Numerical Simulation,2009,14(7):3309-3319

24 Tong D,Wang R,Yang H.Exact solutions for the fl ow of non-Newtonian fl uid with fractional derivative in an annular pipe.Science in China Series G:Physics,Mechanics and Astronomy,2005,48(4):485-495

25 Tan W,Pan W,Xu M.A note on unsteady fl ows of a viscoelastic fl uid with the fractional Maxwell model between two parallel plates.International Journal of Non-Linear Mechanics,2003,38(5):645-650

26 Tong D,Liu Y.Exact solutions for the unsteady rotational fl ow of non-Newtonian fl uid in an annular pipe.International Journal of Engineering Science,2005,43(3):281-289

27 Cai W,Chen W,Xu W.Characterizing the creep of viscoelastic materials by fractal derivative models.International Journal of Non-Linear Mechanics,2016,87:58-63

28 Chen W.Time–space fabric underlying anomalous di ff usion.Chaos,Solitons&Fractals,2006,28(4):923-929

29 Chen W,Sun H,Zhang X,et al.Anomalous di ff usion modeling by fractal and fractional derivatives.Computers&Mathematics with Applications,2010,59(5):1754-1758

30 Liang Y,Ye AQ,Chen W,et al.A fractal derivative model for the characterization of anomalous di ff usion in magnetic resonance imaging.Communications in Nonlinear Science&Numerical Simulation,2016,39:529-537

31 Chen W,Zhang X,Koro?sak D.Investigation on fractional and fractal derivative relaxation-oscillation models.International Journal of Nonlinear Sciences&Numerical Simulation,2010,11(1):3-10

32 Reyesmarambio J,Moser F,Gana F,et al.A fractal time thermal model for predicting the surface temperature of air-cooled cylindrical Li-ion cells based on experimental measurements.Journal of Power Sources,2016,306:636-645

33 Su X,Chen W,Xu W.Characterizing the rheological behaviors of non-Newtonian fl uid via a viscoelastic component:Fractal dashpot.Advances in Mechanical Engineering,2017,9(3):1-12

34 王惠民.流體力學(xué)基礎(chǔ).北京：清華大學(xué)出版社,2013(Wang Huiming.Basis of fl uid mechanics.Beijing:Tsinghua University Press,2013(in Chinese))

NUMERICAL STUDY FOR LAMINAR FLOW OF NON-NEWTONIAN FLUID BASED ON FRACTAL DERIVATIVE1)

Su Xianglong Xu Wenxiang Chen Wen2)
(Institute of Soft Matter Mechanics，College of Mechanics and Materials，Hohai University，Nanjing 211100，China)

Non-Newtonian fl uid has complex rheological characteristics.It is very helpful to reveal these characteristics for the applications of non-Newtonian fl uid in industry and agriculture.The classical rheological models of non-Newtonian fl uid usually have sophisticated forms and the limitations of speci fi c materials or rheological situations.Fractional models have been successfully applied to describe the motion of non-Newtonian fl uid due to their simplicity and few parameters.As an alternative method,the Hausdor fffractal derivative possesses simpler form and higher computational efficiency compared with the fractional derivative.This paper proposes a fractal dashpot model that improves the current Newton’s Law by using the Hausdor fffractal derivative.By investigating the apparent viscosity,the creep and recovery characteristics of the fractal dashpot,it shows that the proposed fractal dashpot model is suitable to describe the non-Newtonian fl uid with viscoelasticity(the so-called fractal fl uid).Combined the fractal dashpot model with the continuity and motion equations,the basic equation for the fractal fl uid for the laminar fl ow between two parallel plates is derived.Moreover,the velocity distributions between two plates are numerically calculated in three cases,which can be obtained through whether there is horizontal pressure gradient or the initial velocity of upper plate.It is found that the horizontal pressure gradient can change the shape of velocity over time and delay the arrival of stable velocity.The fractal fl uid with di ff erent orders has the same velocity distribution and evolution when the horizontal pressure gradient doesn’t exist.In addition,the velocity of upper plate doesn’t in fl uence the di ff erence of stable velocity between di ff erent orders of fractal fl uid when the horizontal pressure gradient exists.

fractal derivative,non-Newtonian fl uid,fractal fl uid,laminar fl ow,apparent viscosity

O373

A

10.6052/0459-1879-16-318

2016–11–07收稿，2017–07–17 錄用,2017–07–19 網(wǎng)絡(luò)版發(fā)表.

1)國(guó)家自然科學(xué)基金(11372097,11402076,11572111)，江蘇省自然科學(xué)基金(BK20130841)和“111”計(jì)劃(B12032)資助項(xiàng)目.2)陳文，教授，主要研究方向：分?jǐn)?shù)階及分形導(dǎo)數(shù)對(duì)軟物質(zhì)建模,無(wú)網(wǎng)格方法研究.E-mail：chenwen@hhu.edu.cn

蘇祥龍,許文祥,陳文.基于分形導(dǎo)數(shù)對(duì)非牛頓流體層流的數(shù)值研究.力學(xué)學(xué)報(bào),2017,49(5):1020-1028

Su Xianglong,Xu Wenxiang,Chen Wen.Numerical study for laminar fl ow of non-Newtonian fl uid based on fractal derivative.Chinese Journal of Theoretical and Applied Mechanics,2017,49(5):1020-1028