薩克遜碗沉沒耗時(shí)的數(shù)值求解

張廷鈞，石昕宇，王徐升

(南京大學(xué) 匡亞明學(xué)院，江蘇 南京 210046)

盎格魯-薩克遜人這一概念最早出現(xiàn)于公元5世紀(jì)[1]，屆時(shí)西方缺乏計(jì)時(shí)器具，而薩克遜人發(fā)明了一種能夠計(jì)時(shí)的碗[2]，薩克遜碗因而得名.薩克遜碗的底部中心處有一個(gè)孔，如果放在水面上，隨著水從小孔流入，碗會(huì)逐漸下沉.碗的運(yùn)動(dòng)過程由流體力學(xué)原理決定，其沉沒耗時(shí)與碗的質(zhì)量、底部小孔孔徑、高度、外徑以及液體密度等參量有關(guān)[3].流體力學(xué)方程形式復(fù)雜，往往難以找到解析解，因此我們進(jìn)行理論分析建模，并編寫程序數(shù)值求解了薩克遜碗的沉沒耗時(shí)，并與實(shí)驗(yàn)結(jié)果進(jìn)行對比.

本文的研究內(nèi)容來源于2020年中國大學(xué)生物理學(xué)術(shù)競賽賽題[4]，主要工作如下：1) 分析薩克遜碗的運(yùn)動(dòng)過程，列出動(dòng)力學(xué)方程，編寫數(shù)值求解程序；2) 改變參量，研究碗沉沒耗時(shí)的變化情況，并與實(shí)驗(yàn)結(jié)果進(jìn)行對比；3) 對本文建立的模型進(jìn)行評價(jià).

作為一道競賽賽題，一些已經(jīng)發(fā)表的學(xué)術(shù)論文和論壇網(wǎng)站也對該問題進(jìn)行了討論[5-7].相比之下，本文的工作具有如下優(yōu)點(diǎn)和創(chuàng)新：1) 較為系統(tǒng)地從理論計(jì)算和實(shí)驗(yàn)兩方面，研究了孔徑、質(zhì)量和液體密度等多個(gè)參量對薩克遜碗沉沒耗時(shí)的影響.先以實(shí)驗(yàn)數(shù)據(jù)檢驗(yàn)理論模型的正確性，再通過數(shù)值計(jì)算探索更為普遍的情況；2) 討論了本文理論模型的適用范圍，在該范圍內(nèi)計(jì)算結(jié)果能精確地描述實(shí)驗(yàn)數(shù)據(jù)；3) 通過對連續(xù)性方程進(jìn)行簡單的修正，即取得較好的效果；4) 定量考慮了形狀阻力和摩擦阻力的影響.

1 模型建立、條件設(shè)定以及說明

1.1 模型建立與條件設(shè)定

圖1為圓柱形薩克遜碗在容器中下沉的剖面示意圖.其中最外側(cè)的黑色豎線表示容器壁，中部對稱的兩個(gè)L形部分為碗壁和碗底，顏色較淺的部分表示液體，實(shí)際情況下容器的尺寸比碗大得多.圖中h0為碗的總高度，H表示碗沉入液體中的深度，h表示碗內(nèi)液面高度，r為孔徑，R為外徑，T指碗的壁厚而T′指碗的底厚.設(shè)定x軸正方向向下.將碗和碗中的液體視為一個(gè)整體，整體的質(zhì)心坐標(biāo)為

(1)

其中mb是碗的質(zhì)量，ρ是液體密度.xb為碗的質(zhì)心坐標(biāo)，它與H只相差一個(gè)常量(該常量可以由碗的形狀計(jì)算得到).

圖1 薩克遜碗下沉的剖面示意圖與各符號含義

接下來對液體從碗底小孔涌入碗這一過程進(jìn)行分析.在實(shí)際應(yīng)用的情況下，薩克遜碗所在液體的表面積遠(yuǎn)大于碗的橫截面積，因此可以認(rèn)為液面高度不變，且除了碗下方的部分液體之外都是靜止的.實(shí)驗(yàn)中，絕大多數(shù)碗下沉速度較慢，不超過3 mm/s，且沒有運(yùn)動(dòng)狀態(tài)突變，因此可認(rèn)為流體定常流動(dòng)，且沿流線流動(dòng).此外實(shí)際情況下顯然可以認(rèn)為液體不可壓縮.因此對碗內(nèi)與碗外液面處的液體列出伯努利方程，再考慮小孔處碗與液體之間的摩擦，對方程進(jìn)行修正：

(2)

其中Patm為大氣壓強(qiáng).v0應(yīng)該理解如下：液體通過小孔流入碗、且在沒有徑向擴(kuò)散的情況下到達(dá)碗內(nèi)液面時(shí)的速度.或通過如下假想來理解：有一根直徑與孔徑相同的硬質(zhì)管道連接了小孔和內(nèi)液面，液體從該管道上端流出時(shí)的速度即為v0.考慮到液體近似靜止，式(2)中已經(jīng)忽略了液體初流速.

由于液體不可壓縮，單位時(shí)間通過同一流管的液體體積恒定，通過連續(xù)性方程，可以將上述假想的管道上端處液體流速與內(nèi)液面上升速度聯(lián)系起來，有

(3)

聯(lián)立式(2)和式(3)，得

(4)

(5)

其中CD是速度矯正項(xiàng)，它與碗的特征和液體黏滯性等因素有關(guān)，其物理意義將于下一小節(jié)討論.

在下沉過程中，碗除了受到重力、浮力和表面張力之外，還會(huì)受到阻力，包括形狀阻力和摩擦阻力.根據(jù)Blasius邊界層理論，阻力的表達(dá)式為[8,9]

(6)

其中A為碗的橫截面積，μ為黏度，經(jīng)驗(yàn)參數(shù)cf可以取0.82.

綜合考慮上述各類受力，對碗及碗中液體組成的整體列出動(dòng)力學(xué)方程：

(7)

其中

mc=mb+ρπ(R-T)2h

(8)

指整體的質(zhì)量.式(7)右邊的最后一項(xiàng)為液體對碗的表面張力，γ為液體的表面張力系數(shù)，θ是液體與碗的接觸角.化簡式(7)得

(9)

對式(1)、式(5)、式(9)進(jìn)行聯(lián)合求解，即可得到薩克遜碗的運(yùn)動(dòng)情況.

本文使用的初始條件為碗的位置使重力與浮力平衡(碗內(nèi)沒有液體)，且速度為0.當(dāng)碗的上沿到達(dá)水面(即H=h0)時(shí)，認(rèn)為碗沉沒.為了簡單起見，設(shè)液體與碗的接觸角為一定值(這其實(shí)是與實(shí)驗(yàn)事實(shí)不符的).

1.2 關(guān)于伯努利方程速度矯正項(xiàng)CD的說明

事實(shí)上，修正系數(shù)CD有明確物理意義，它衡量了液體流速沿徑向不均勻分布的情況.由于碗底小孔的邊緣會(huì)對向上涌入的水施加作用[3]，所以接近孔緣的液體流速低于中心處.而伯努利方程和連續(xù)性方程均默認(rèn)液體流速沿徑向均勻分布，因此需要將液體流速(或速度的平方)修改為平均速度或速度的方均值，有

(10)

以及

(11)

(12)

將上式與式(5)對比，只要令

(13)

2 碗沉沒耗時(shí)的計(jì)算及其與實(shí)驗(yàn)對比

本章節(jié)研究了各參量變化時(shí)碗的沉沒耗時(shí)，這些參量包括碗的孔徑、質(zhì)量和液體密度等.部分條件下還開展實(shí)驗(yàn)，將實(shí)驗(yàn)數(shù)據(jù)與計(jì)算結(jié)果進(jìn)行對比.

實(shí)驗(yàn)部分采用了控制變量法，使用的所有薩克遜碗為同種塑料材質(zhì)，密度相同且均勻，約為1 600 kgm3.對于相同外徑的碗，高度、壁厚和內(nèi)徑等參量均相同，僅有孔徑可能不同，這對質(zhì)量造成的影響不超過1.5%，故本文近似認(rèn)為它們的質(zhì)量相同.內(nèi)徑為90 mm、75 mm和63 mm的空碗，質(zhì)量可分別取為157.5 g、98.5 g和64.5 g(數(shù)值計(jì)算即采用這些數(shù)值)，實(shí)際情況下各種孔徑的碗質(zhì)量與此偏差小于1 g.探究沉沒耗時(shí)與孔徑和質(zhì)量關(guān)系時(shí)，使用的液體為清水，密度恒定為997 kgm3.

2.1 沉沒耗時(shí)與孔徑的關(guān)系

對內(nèi)徑為90 mm和75 mm的碗，本小節(jié)研究了孔徑從2 mm變化到8 mm過程中沉沒耗時(shí)的變化情況，并與實(shí)驗(yàn)數(shù)據(jù)作了對比，結(jié)果如圖2所示.

圖2 沉沒耗時(shí)與孔徑的關(guān)系

對內(nèi)徑為90 mm的碗，固定CD=0.68；對內(nèi)徑為75 mm的碗，固定CD=0.645，在這一組參數(shù)下，數(shù)值計(jì)算的結(jié)果與實(shí)驗(yàn)數(shù)據(jù)符合得較好.但是對于內(nèi)徑為63 mm的碗，很難有一個(gè)固定的CD取值能很好地對各個(gè)小孔半徑的實(shí)驗(yàn)結(jié)果進(jìn)行描述.這可能是模型的一個(gè)缺陷，我們將在后文對此進(jìn)行討論.

2.2 沉沒耗時(shí)與質(zhì)量的關(guān)系

在實(shí)驗(yàn)過程中，可以在碗內(nèi)均勻擺放砝碼以增加碗的質(zhì)量，同時(shí)不影響碗的平衡.對內(nèi)徑為90 mm和75 mm，孔徑分別為2 mm、4 mm、6 mm和8 mm的8個(gè)碗，取CD=0.67，計(jì)算了沉沒耗時(shí)隨質(zhì)量的變化情況.其中內(nèi)徑為90 mm、孔徑為2 mm和4 mm的兩組進(jìn)行了實(shí)驗(yàn)，質(zhì)量增量最高達(dá)到120 g.值得注意的是，實(shí)驗(yàn)中只能增加碗的質(zhì)量，而數(shù)值計(jì)算也對碗質(zhì)量減小的情形進(jìn)行了研究.結(jié)果如圖3所示.

圖3 沉沒耗時(shí)與質(zhì)量的關(guān)系

在這組參數(shù)下，數(shù)值計(jì)算的結(jié)果與實(shí)驗(yàn)數(shù)據(jù)符合得很好.由圖3可見，隨著碗質(zhì)量的減小，沉沒耗時(shí)增加得很快；而隨著碗質(zhì)量的增加，沉沒耗時(shí)逐漸減小并趨于0.這是符合預(yù)期的：當(dāng)碗的質(zhì)量足夠小時(shí)，由于浮力和表面張力的作用，碗將會(huì)漂浮在水面上而不下沉.而當(dāng)碗的質(zhì)量很大時(shí)，即使碗內(nèi)沒有液體且碗的上沿與液面持平，浮力和表面張力也不足以平衡重力，初始條件無法滿足，碗在釋放后立即沉沒.但是數(shù)值計(jì)算的結(jié)果卻沒有揭示實(shí)驗(yàn)中觀察到的一種特殊現(xiàn)象，這將在后文中進(jìn)行討論.

2.3 沉沒耗時(shí)與液體密度的關(guān)系

取若干種參量的碗，研究碗所處液體密度變化時(shí)沉沒耗時(shí)的變化情況，其中內(nèi)徑為63 mm、孔徑為2 mm(取CD=0.7)和內(nèi)徑為90 mm、孔徑為4 mm(取CD=0.68)的碗還與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對比，結(jié)果如圖4所示.

較小范圍

較大范圍圖4 沉沒耗時(shí)與液體密度的關(guān)系

數(shù)值計(jì)算的結(jié)果也能準(zhǔn)確地反映實(shí)驗(yàn)數(shù)據(jù).在圖4(a)中，沉沒耗時(shí)看似隨著液體密度線性上升，其實(shí)不然，如果液體密度繼續(xù)增加，沉沒耗時(shí)會(huì)增加得更快.當(dāng)液體密度大于1 600 kg/m3時(shí)，數(shù)值計(jì)算的結(jié)果出現(xiàn)異常.這是因?yàn)榇藭r(shí)液體密度已經(jīng)超過了碗的密度，碗將一直漂浮在液面上.而當(dāng)液體密度減小到375 kg/m3以下時(shí)，液體密度不足以平衡碗的重力，碗被釋放之后立即沉沒，數(shù)值計(jì)算會(huì)得出負(fù)的結(jié)果.

3 對模型及結(jié)果的討論

3.1 模型的精確性與局限性

從前文可見，適當(dāng)?shù)厝?shù)CD，理論模型能精確地描述碗的沉沒耗時(shí)隨碗的質(zhì)量和液體密度變化的情況.在研究沉沒耗時(shí)隨小孔半徑變化時(shí)，從結(jié)果來看，若外徑始終在孔徑的10倍以上，數(shù)值計(jì)算結(jié)果較為精確；而對于外徑較小的碗，很難找到一個(gè)合適的參數(shù)CD精確描述這一變化趨勢，甚至?xí)霈F(xiàn)沉沒時(shí)間小于0這樣的非物理結(jié)果.我們認(rèn)為這是模型內(nèi)在局限性的結(jié)果：隨著孔徑的增加或外徑的減小，液體的流動(dòng)模式向湍流轉(zhuǎn)變；而伯努利方程并不能用于描述湍流，即使作了式(2)的簡單修正效果也不理想.此外，當(dāng)孔徑較大時(shí)，碗相對液體的運(yùn)動(dòng)速度較快，而模型沒有考慮液體的初流速，這未必是合理的.即使在更大的范圍內(nèi)調(diào)整參數(shù)CD，計(jì)算結(jié)果依然與實(shí)驗(yàn)數(shù)據(jù)有偏差.

3.2 參數(shù)的普適性與不敏感性

如前文所述，在研究不同的碗、改變不同的物理量時(shí)，參數(shù)CD僅需要微調(diào)，就能使計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)相符.這說明參數(shù)CD具有普適性，在已知條件不足且對精度要求不高時(shí)，可以考慮取某個(gè)固定的CD值(如0.675)對實(shí)驗(yàn)結(jié)果做出大致的預(yù)測.

此外，模型還是對參數(shù)CD不敏感的，在微調(diào)參數(shù)CD時(shí)，數(shù)值計(jì)算的結(jié)果也僅有微小變化，且都能很好地?cái)M合實(shí)驗(yàn)數(shù)據(jù).為了對此進(jìn)行驗(yàn)證，對內(nèi)徑為90 mm的碗，使CD從0.65以0.01為步長變化到0.7，計(jì)算沉沒耗時(shí)隨孔徑的變化如圖5所示.

圖5 模型穩(wěn)定性檢驗(yàn)

3.2 模型對特殊現(xiàn)象的描述

特殊現(xiàn)象一方面指極端情況，即參量特別大或特別小的情況.如前文所述，當(dāng)質(zhì)量取極端值時(shí)，模型的表現(xiàn)良好：當(dāng)碗的質(zhì)量增大時(shí)，沉沒耗時(shí)趨于0，質(zhì)量減小時(shí)，沉沒耗時(shí)快速地增長；而當(dāng)液體密度超過碗的密度時(shí)，沉沒耗時(shí)的異常值也反映了碗不下沉的情況.

然而當(dāng)碗底小孔半徑取極端值時(shí)，該模型無法給出合理的結(jié)果.若孔徑很大，則由于出現(xiàn)湍流等原因，碗的運(yùn)動(dòng)規(guī)律完全不同，該模型失效；若孔徑很小，則碗會(huì)由于表面張力作用不下沉，示意圖及實(shí)驗(yàn)中拍攝的照片如圖6所示.

圖6 孔徑很小時(shí)碗不下沉現(xiàn)象的圖示

對這一現(xiàn)象可以作簡單的理論解釋如下：重力和浮力平衡，即

mg=ρghπR2

(14)

在小孔處，表面張力產(chǎn)生的力大于液體壓力，有

ρghπr2≤2πrσ

(15)

由以上二式得

(16)

這表明，當(dāng)r較小時(shí)，表面張力產(chǎn)生了堵塞小孔的作用，阻止液體流入，碗不會(huì)下沉，而模型并沒有揭示這一機(jī)制.類似地，當(dāng)m較小時(shí)也有同樣的情況出現(xiàn)，模型雖然展示了碗質(zhì)量減小時(shí)不會(huì)下沉的情況，但主要基于重力與浮力關(guān)系，與該機(jī)制不同.

4 總結(jié)

基于伯努利方程和邊界層理論，本文對薩克遜碗的下沉過程進(jìn)行了物理建模以及數(shù)值求解.將數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，兩者基本相符.在此基礎(chǔ)上，本文擴(kuò)展了參量的變化范圍，從而對更廣泛的情況進(jìn)行探索，進(jìn)而得到了一些有物理意義的結(jié)果.除此以外，本文還對模型及結(jié)果進(jìn)行了討論和評價(jià)，指出模型的優(yōu)缺點(diǎn).具體總結(jié)如下：

1) 該模型主要適用于外徑與碗底孔徑比值較大的薩克遜碗(10倍以上)，此時(shí)液體的流動(dòng)模式為層流.在適用范圍內(nèi)，模型能精確描述碗的沉沒耗時(shí)隨孔徑、碗的質(zhì)量和液體密度的變化.在此范圍之外，很難通過一個(gè)固定參數(shù)CD描述各種情況.

2) 模型正確地揭示了一些極端情況下的物理結(jié)果.一是隨著碗的質(zhì)量增大，碗的沉沒耗時(shí)趨于0.二是若液體密度超過碗的密度，碗不下沉.

3) 在模型適用范圍內(nèi)，參數(shù)CD具有良好的不敏感性和普適性.一方面，對CD進(jìn)行微調(diào)幾乎不影響計(jì)算結(jié)果，且都與實(shí)驗(yàn)數(shù)據(jù)符合；另一方面，無論是改變孔徑、質(zhì)量還是液體密度，取值相近的CD都能對各種碗的沉沒耗時(shí)變化做出擬合.

4) 模型對實(shí)驗(yàn)中的一類特殊情況未能做出描述：當(dāng)碗質(zhì)量或孔徑較小時(shí)，表面張力阻止液體流入碗，進(jìn)而使得碗不下沉.