曹占雪，王佳棟，溫國(guó)正

（江蘇大學(xué) 土木工程與力學(xué)學(xué)院，江蘇 鎮(zhèn)江 212013）

流體晃動(dòng)是指部分充液的儲(chǔ)液罐在外部激勵(lì)作用下內(nèi)部流體產(chǎn)生運(yùn)動(dòng)的現(xiàn)象。在實(shí)際工程中流體晃動(dòng)產(chǎn)生的液動(dòng)壓力作用于結(jié)構(gòu)，對(duì)其產(chǎn)生不可忽略的影響，例如唐山大地震中流體晃動(dòng)使罐壁與底板之間的角焊處發(fā)生斷裂破壞從而導(dǎo)致儲(chǔ)油的泄露[1]，Tokachi-oki地震中流體晃動(dòng)使儲(chǔ)液罐浮頂下沉從而引發(fā)火災(zāi)[2]等。在工程上通常采用隔板降低流體晃動(dòng)的幅度，從而減小了流體晃動(dòng)對(duì)結(jié)構(gòu)的影響；在理論上為獲得流體晃動(dòng)的各階模態(tài)并服務(wù)于后續(xù)的動(dòng)力學(xué)分析，有必要進(jìn)行流體晃動(dòng)的特性研究，因此研究?jī)?chǔ)液罐中的流體晃動(dòng)特性具有一定的工程應(yīng)用和理論價(jià)值。

通過(guò)模態(tài)分析主要可獲得儲(chǔ)液系統(tǒng)的固有頻率和各階陣型，通過(guò)模態(tài)分析可以避免激勵(lì)頻率和儲(chǔ)液系統(tǒng)的固有頻率發(fā)生共振，從而為在工程上的應(yīng)用提供依據(jù)。目前已有許多文獻(xiàn)研究了隔板對(duì)流體晃動(dòng)的影響，其中就包括模態(tài)分析，并已取得了豐碩的理論成果。Chu 等[3]研究了多個(gè)相同隔板條件下的流體晃動(dòng)問(wèn)題，發(fā)現(xiàn)系統(tǒng)的固有頻率隨著隔板上方的有效水深的變化而變化。Ebrahimian 等[4]用邊界元法獲得了軸對(duì)稱儲(chǔ)液系統(tǒng)的固有頻率和振型，其研究表明隔板的位置對(duì)固有頻率的影響更大。Wang等[5]對(duì)帶有單個(gè)或多個(gè)相同隔板的環(huán)形儲(chǔ)液系統(tǒng)進(jìn)行了研究，發(fā)現(xiàn)隔板的位置主要對(duì)第1 階模態(tài)的影響比較大。針對(duì)帶單層隔板的渡槽，房忠潔等[6]用半解析方法說(shuō)明隔板對(duì)速度勢(shì)分布的影響比較大。有些學(xué)者對(duì)圓柱形儲(chǔ)液系統(tǒng)進(jìn)行了大量的研究，如楊唱等[7]通過(guò)對(duì)其研究發(fā)現(xiàn)隔板參數(shù)與防晃效果之間的關(guān)系。同樣地，Wang 等[8]分析了多層相同隔板對(duì)流體晃動(dòng)效應(yīng)的影響。隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，數(shù)值模擬也經(jīng)常被用來(lái)解決帶隔板的儲(chǔ)液系統(tǒng)中流體晃動(dòng)問(wèn)題。陳更等[9]對(duì)帶隔板的液艙進(jìn)行了數(shù)值模擬，其研究表明縱隔板可以更好抑制液艙中流體的晃動(dòng)。Jin 等[10]通過(guò)數(shù)值模擬研究了水平多孔隔板附近的流體晃動(dòng)特性問(wèn)題，為理論研究提供了一個(gè)依據(jù)。

綜上所述，這些文獻(xiàn)主要研究了單個(gè)或多個(gè)相同隔板對(duì)流體晃動(dòng)的影響，鮮有研究?jī)蓚€(gè)不同半徑環(huán)形剛性隔板的控制機(jī)理，而文獻(xiàn)[11]中的流體子域法是一種半解析法，能夠同時(shí)滿足求解問(wèn)題的高效性和高精度，因此本文基于流體子域法研究了兩個(gè)不同半徑環(huán)形剛性隔板對(duì)流體晃動(dòng)特性的影響。

文獻(xiàn)[12]指出存在多層相同隔板的情況下，下層隔板對(duì)流體晃動(dòng)的影響比較小，故本文在此主要考慮兩層隔板情況，并且上隔板的內(nèi)徑大于下隔板，其目的是為了凸顯下隔板對(duì)于晃動(dòng)特性的影響。為了使所有的子域具有連續(xù)的邊界條件，引入8 個(gè)人工界面將儲(chǔ)液罐中的流體區(qū)域劃分為8 個(gè)子域，利用疊加原理和分離變量法求出各子域的速度勢(shì)的試函數(shù)，將其代入子域間的連續(xù)條件與自由液面的動(dòng)力學(xué)邊界條件，得到關(guān)于待定系數(shù)與晃動(dòng)頻率的級(jí)數(shù)方程，通過(guò)廣義傅里葉展開(kāi)將級(jí)數(shù)方程轉(zhuǎn)換成線性方程組，從而得到關(guān)于頻率與模態(tài)的特征方程。利用求解廣義特征值的方法確定待定系數(shù)，從而使問(wèn)題得到解決。

1 基本方程

圖1 為帶兩層隔板的圓柱形儲(chǔ)液罐，儲(chǔ)液罐內(nèi)半徑為R3，上隔板內(nèi)半徑為R2，下層隔板內(nèi)半徑為R1。如圖1 所示。R1、R2、R3滿足關(guān)系：R1＜R2＜R3。本文主要研究?jī)?chǔ)液罐中流體的晃動(dòng)特性，因此罐體與隔板均假設(shè)為剛性。以儲(chǔ)液罐底板的中心為坐標(biāo)原點(diǎn)，建立如圖1所示的柱體坐標(biāo)系Orθz，上層隔板的位置設(shè)為z2，下層隔板的位置設(shè)為z1，液面高度設(shè)為H。

圖1 帶兩層隔板的圓柱形儲(chǔ)液罐

按照?qǐng)D2 所示的子域劃分方法，將儲(chǔ)液罐中的流體區(qū)域劃分成8 個(gè)子域Ωi(i=1,2,…,8)，Sk(k=1,2…,8)為相鄰子域間的人工界面。將子域Ωi的速度勢(shì)函數(shù)設(shè)為φi，由于儲(chǔ)液罐中的流體為無(wú)旋、無(wú)黏、不可壓縮的理想流體，因此任意子域的速度勢(shì)函數(shù)滿足如下拉普拉斯方程：

在柱坐標(biāo)系下，流體域Ωi中每一點(diǎn)的速度為[12]：

如圖2所示。流體子域在流-固耦合界面上需滿足固壁邊界條件，因此各子域的速度勢(shì)函數(shù)在罐壁上應(yīng)滿足如下邊界條件：

圖2 流體子域的劃分

在罐底應(yīng)滿足如下邊界條件：

在隔板上應(yīng)滿足如下邊界條件：

柱形子域的速度勢(shì)函數(shù)應(yīng)滿足如下邊界條件：

本文考慮流體做微幅線性晃動(dòng)，根據(jù)Bernoulli方程，流體子域Ω6、Ω7以及Ω8在自由液面上應(yīng)滿足如下邊界條件：

相鄰子域在其界面上應(yīng)該滿足壓力及速度連續(xù)條件，即有：

式中：n為人工界面Sk上的法向單位矢量；i、q為相鄰的子域。

2 子域速度勢(shì)函數(shù)求解

當(dāng)儲(chǔ)液罐中流體發(fā)生自由晃動(dòng)時(shí)，可以將流體速度勢(shì)函數(shù)設(shè)為如下形式[13]：

在罐底應(yīng)滿足如下邊界條件：

在隔板上應(yīng)滿足如下邊界條件：

位于中心區(qū)域的速度勢(shì)函數(shù)應(yīng)滿足如下邊界條件：

Φ2|r=0=有限值,Φ5|r=0=有限值。

流體子域Ω6、Ω7以及Ω8的振型函數(shù)在自由液面上應(yīng)滿足如下邊界條件：

相鄰子域在其界面上應(yīng)該滿足壓力及速度連續(xù)條件，即有：

式中：n為人工界面Sk上的法向單位矢量；i、q為相鄰的子域。

根據(jù)疊加原理，對(duì)各子域的振型函數(shù)進(jìn)行分解，結(jié)果如圖3 所示，圖3(a)所示為環(huán)柱形子域，圖3(b)所示為柱形子域，圖3(a)中的細(xì)點(diǎn)畫(huà)線表示子域的對(duì)稱軸。子域6、7、8 中的上實(shí)線為零壓力邊界條件，其余所有實(shí)線均為剛性邊界條件，虛線代表非齊次性邊界條件，即有：

圖3 振型函數(shù)分解示意圖

根據(jù)式(16)，對(duì)于子域Ω1，振型函數(shù)的分量Φ11應(yīng)滿足如下邊界條件：

對(duì)于子域Ω2，振型函數(shù)的分量Φ12、Φ22應(yīng)滿足如下邊界條件：

對(duì)于子域Ω3，振型函數(shù)的分量Φ13應(yīng)滿足如下邊界條件：

對(duì)于子域Ω4，振型函數(shù)的分量Φ14、Φ24、Φ34應(yīng)滿足如下邊界條件：

對(duì)于子域Ω5，振型函數(shù)的分量Φ15、Φ25、Φ35應(yīng)滿足如下邊界條件：

對(duì)于子域Ω6，振型函數(shù)的分量Φ16、Φ26應(yīng)滿足如下邊界條件：

對(duì)于子域Ω7，振型函數(shù)的分量Φ17、Φ27、Φ37、Φ47應(yīng)滿足如下邊界條件：

對(duì)于子域Ω8，振型函數(shù)Φ18、Φ28、Φ38應(yīng)滿足如下邊界條件：

根據(jù)式(17)至式(24)，利用分離變量法可求得各個(gè)分量的試函數(shù)，再根據(jù)式(16)可得各個(gè)子域振型函數(shù)的形式解。將振型函數(shù)代入式(14)～式(15)中可以得到19 個(gè)級(jí)數(shù)線性方程組。根據(jù)三角函數(shù)和貝塞爾函數(shù)的正交性，消除級(jí)數(shù)方程中的空間坐標(biāo)，并將方程截?cái)嘀罭，從而得到有關(guān)待定系數(shù)的線性方程組，即有：

式中，矩陣V、D中的非零元素為級(jí)數(shù)方程展開(kāi)后的系數(shù)；Λ為無(wú)量綱化后的晃動(dòng)頻率，滿足：

級(jí)數(shù)方程(25)中的Λ2定義為廣義特征值，通過(guò)求解廣義特征值的方法確定和的數(shù)值。將代入振型函數(shù)就可以確定各個(gè)子域的速度勢(shì)函數(shù)。

3 收斂性和比較研究

首先進(jìn)行收斂性分析，表1和表2分別對(duì)應(yīng)于上層隔板參數(shù)變化時(shí)的前3 階晃動(dòng)頻率，從表中可以發(fā)現(xiàn)，隨著截?cái)囗?xiàng)數(shù)的增加，晃動(dòng)頻率會(huì)逐漸趨于精確解。雖然隔板參數(shù)對(duì)晃動(dòng)頻率的收斂性會(huì)有一定的影響，但這種影響會(huì)隨著截?cái)囗?xiàng)數(shù)的增加而逐漸減小，當(dāng)截?cái)囗?xiàng)數(shù)足夠大時(shí)，該影響可忽略不計(jì)。當(dāng)N=26時(shí)，表1和表2中的晃動(dòng)頻率均可保證4位有效數(shù)字，因此在下面的分析中，截?cái)囗?xiàng)數(shù)均取26。

表1 (l=0,1,n=1,2,3)隨截?cái)囗?xiàng)數(shù)N變化的收斂性（R1=0.2/m;R3=1/m;z1=0.4/m;z2=0.8/m;H=1/m)

表1 (l=0,1,n=1,2,3)隨截?cái)囗?xiàng)數(shù)N變化的收斂性（R1=0.2/m;R3=1/m;z1=0.4/m;z2=0.8/m;H=1/m)

R2 0.4/m 0.6/m l 0 1 0 1 n 1 2 3 1 2 3 1 2 3 1 2 3 N=21 2.543 6.478 9.981 0.783 4.722 8.219 2.846 6.744 10.029 1.028 4.864 8.349 N=22 2.546 6.480 9.983 0.783 4.722 8.219 2.848 6.748 10.031 1.028 4.864 8.349 N=23 2.546 6.478 9.982 0.783 4.722 8.219 2.846 6.744 10.030 1.029 4.864 8.349 N=24 2.546 6.480 9.984 0.783 4.722 8.219 2.848 6.748 10.031 1.029 4.864 8.349 N=25 2.546 6.478 9.983 0.783 4.722 8.219 2.847 6.744 10.031 1.029 4.864 8.349 N=26 2.546 6.480 9.984 0.783 4.722 8.219 2.848 6.748 10.032 1.029 4.864 8.349 N=27 2.546 6.479 9.983 0.783 4.722 8.219 2.847 6.745 10.031 1.029 4.864 8.349

表2 (l=0,1,n=1,2,3)隨截?cái)囗?xiàng)數(shù)N變化的收斂性（R1=0.3/m;R2=0.6/m;R3=1/m;z1=0.2/m;H=1/m）

表2 (l=0,1,n=1,2,3)隨截?cái)囗?xiàng)數(shù)N變化的收斂性（R1=0.3/m;R2=0.6/m;R3=1/m;z1=0.2/m;H=1/m）

z2 0.4/m 0.8/m l 0 1 0 1 n 1 2 3 1 2 3 1 2 3 1 2 3 N=21 3.772 6.996 10.137 1.598 5.315 8.529 2.851 6.749 10.030 1.063 4.918 8.349 N=22 3.774 6.998 10.141 1.598 5.316 8.529 2.853 6.753 10.031 1.063 4.918 8.350 N=23 3.773 6.999 10.144 1.598 5.317 8.530 2.851 6.750 10.031 1.063 4.918 8.350 N=24 3.775 7.001 10.147 1.598 5.317 8.531 2.853 6.754 10.031 1.063 4.918 8.350 N=25 3.774 7.002 10.149 1.598 5.318 8.531 2.852 6.751 10.031 1.063 4.918 8.350 N=26 3.775 7.003 10.151 1.598 5.318 8.531 2.853 6.754 10.032 1.064 4.918 8.350 N=27 3.775 7.004 10.153 1.598 5.318 8.532 2.852 6.751 10.032 1.064 4.918 8.350

為驗(yàn)證本文方法的正確性，使用通用有限元軟件ADINA對(duì)本文模型進(jìn)行模態(tài)分析，得到了流體晃動(dòng)的頻率和模態(tài)。表3 對(duì)本文方法解和ADINA 解進(jìn)行了比較。結(jié)果表明，本文計(jì)算的結(jié)果和ADINA有限元軟件模擬的結(jié)果具有較好的一致性，其最大相對(duì)誤差為0.65%，最小相對(duì)誤差為0.02%。因此本文方法滿足求解問(wèn)題的高精度要求。

表3 (l=0,1,n=1,2,3)的比較研究（R1=0.35/m;R3=1/m;z1=0.5/m;z2=0.8/m;H=1/m）

表3 (l=0,1,n=1,2,3)的比較研究（R1=0.35/m;R3=1/m;z1=0.5/m;z2=0.8/m;H=1/m）

R2 0.7/m 0.8/m n 1 2 3 1 2 3 l=0本文解3.088 4 6.765 7 10.089 1 3.367 1 3.783 7 10.091 6 ADINA解3.107 0 6.746 2 10.092 1 3.389 0 3.797 8 10.098 2相對(duì)誤差/（%）0.60 0.16 0.03 0.65 0.21 0.07 l=1本文解1.155 0 4.880 0 8.394 7 1.280 7 5.013 6 8.402 1 ADINA解1.162 3 4.885 0 8.396 2 1.286 7 5.023 0 8.405 5相對(duì)誤差/（%）0.63 0.10 0.02 0.47 0.19 0.04

4 參數(shù)研究

首先研究了上層隔板內(nèi)半徑和上層隔板位置對(duì)流體晃動(dòng)特性的影響，結(jié)果如圖4和圖5所示。從圖中發(fā)現(xiàn)，當(dāng)上層隔板的內(nèi)半徑接近儲(chǔ)液罐的內(nèi)半徑時(shí)，晃動(dòng)頻率將不再受上層隔板的影響，此時(shí)下層隔板對(duì)流體晃動(dòng)的頻率起主導(dǎo)影響。但當(dāng)隔板的位置靠近自由液面時(shí)，圖5中曲線的斜率降低很快，表明上層隔板的位置是影響流體晃動(dòng)最主要的因素。因此對(duì)上層隔板的位置進(jìn)行優(yōu)化設(shè)計(jì)能夠有效減弱流體的晃動(dòng)。

圖4 Λ2l1(l=0,1)隨上層隔板內(nèi)徑的變化曲線（R1=0.2/m;R3=1/m;z1=0.2/m;H=1/m）

圖5 (l=0,1)隨上層隔板位置的變化曲線（R1=0.2/m;R3=1/m;z1=0.2/m;H=1/m）

其次研究了下層隔板內(nèi)半徑和下層隔板位置對(duì)流體晃動(dòng)的影響，結(jié)果如圖6和圖7所示。下層隔板內(nèi)半徑對(duì)晃動(dòng)頻率的影響隨著隔板位置的降低而減小，而這主要因?yàn)楦舭褰咏诠薜讜r(shí)，隔板對(duì)流體晃動(dòng)的影響很小。但當(dāng)下層隔板的內(nèi)半徑逐漸減小時(shí)，對(duì)晃動(dòng)頻率的影響也就越來(lái)越明顯。因此改變下層隔板的內(nèi)半徑仍然能夠明顯降低流體晃動(dòng)的固有頻率，即也可達(dá)到抑制流體晃動(dòng)的目的。

圖6 (l=0,1)隨下層隔板內(nèi)徑的變化曲線（R2=0.9/m;R3=1/m;z2=0.8/m;H=1/m）

圖7 (l=0,1)隨下層隔板位置的變化曲線（R2=0.7/m;R3=1/m;z2=0.85/m;H=1/m）

將求得的流體晃動(dòng)頻率代入方程式(25)就可以確定待定系數(shù)，將各個(gè)子域的振型函數(shù)代入方程式(9)就可以確定各個(gè)子域的速度勢(shì)函數(shù)。圖8和圖9 分別為環(huán)向波數(shù)為0 和1 的前兩階振型。很明顯，當(dāng)環(huán)向波數(shù)為0時(shí)其表現(xiàn)為對(duì)稱振型，當(dāng)環(huán)向波數(shù)為1時(shí)表現(xiàn)為反對(duì)稱振型。

圖8 環(huán)向波數(shù)為0時(shí)的前兩階振型（R1=0.35/m;R2=0.7/m;R3=1/m;z1=0.5/m;z2=0.8/m;H=1/m）

圖9 環(huán)向波數(shù)為1時(shí)的前兩階振型（R1=0.35/m;R2=0.7/m;R3=1/m;z1=0.5/m;z2=0.8/m;H=1/m）

5 結(jié)語(yǔ)

本文的主要?jiǎng)?chuàng)新點(diǎn)在于研究了兩個(gè)不同半徑環(huán)形剛性隔板對(duì)流體晃動(dòng)的影響，文中隔板及罐壁均采用剛性材料，流體為理想流體。通過(guò)流體子域法求解了各個(gè)子域的速度勢(shì)函數(shù)，并與ADINA解進(jìn)行比較驗(yàn)證了本文方法的正確性，進(jìn)而對(duì)流體的晃動(dòng)特性進(jìn)行了研究。本文中的模型只適合于流體做微幅線性晃動(dòng)的情況。雖然本文的研究對(duì)象為帶兩層不同內(nèi)徑環(huán)形剛性隔板的儲(chǔ)液系統(tǒng)，但本文的研究方法可以進(jìn)一步擴(kuò)展到帶有任意個(gè)不同內(nèi)徑隔板的情況。

通過(guò)對(duì)參數(shù)的研究可以得到以下結(jié)論：存在兩層隔板的情況下，下層隔板的內(nèi)徑越小，對(duì)流體晃動(dòng)的影響越明顯，但當(dāng)下層隔板位置靠近罐底時(shí)，隔板內(nèi)徑的變化對(duì)流體晃動(dòng)的影響越來(lái)越?。幌聦痈舭鍍?nèi)徑越小，對(duì)環(huán)向波數(shù)l=0的影響越小。