圓盤反應(yīng)器流場數(shù)值模擬

鄧 斌，戴干策

華東理工大學(xué)化學(xué)工程聯(lián)合國家重點實驗室，上海 200237

圓盤反應(yīng)器起源于1900年Weigand 提出的生物轉(zhuǎn)盤，早期主要應(yīng)用于生物廢水處理，20世紀(jì)初用于合成纖維聚酯。近年來，圓盤接觸器在中低粘物系中的應(yīng)用頗受關(guān)注，包括光催化[1]、檸檬酸生產(chǎn)[2]和精細(xì)化工[3]等領(lǐng)域。旋轉(zhuǎn)圓盤接觸器的基本結(jié)構(gòu)是光滑平面圓盤，部分浸沒于液體中，通過電機(jī)驅(qū)動水平軸帶動垂直圓盤旋轉(zhuǎn)，以此區(qū)別于垂直軸帶動完全浸沒于液體中的水平圓盤。水平圓盤反應(yīng)器物料從中心加入，流量是預(yù)先給定的，這類流動已經(jīng)有比較詳細(xì)的研究。部分浸沒的垂直圓盤，盤面膜厚及圓盤持液量則取決于流體物性、圓盤尺寸及旋轉(zhuǎn)速度等因素，稱為自計量流。從液膜流體動力學(xué)角度而言，其顯著的特點是物體從液體中穿過，形成自由面，帶動流體運動，而非物體作潛體運動。這類研究最早可以追溯到上世紀(jì)40年代Landau 等[4]平板抽出模型的經(jīng)典工作，其后為開發(fā)涂布技術(shù)，有相當(dāng)廣泛的涂布流研究，主要探索抽出速度和膜厚的關(guān)系。

Landau 等[4]通過簡化Navier-Stokes 方程最早提出了平板抽出模型，得到了平板帶動液膜的厚度與流速和物性的關(guān)系，但當(dāng)抽出速度較大時，關(guān)系式明顯偏離實驗結(jié)果。White 等[5]考慮了重力因素對膜厚的影響。Groenveld[6]系統(tǒng)研究了影響平板膜厚的各個因素，指出平板液膜厚度主要取決于物系粘度和平板抽出速度。Bintanja 等[7]研究圓盤反應(yīng)器氧氣傳遞特性時受到Groenveld 的啟發(fā)，用線速度代替平板速度得到了圓盤反應(yīng)器平均膜厚公式，但當(dāng)轉(zhuǎn)速較高時，膜厚與線速度的關(guān)系并非文中提出的線性關(guān)系。Zhevalkink 等[8]根據(jù)平板抽出模型，結(jié)合圓盤反應(yīng)器忽略離心力和毛細(xì)管力的影響，并考察了液位對膜厚的影響，簡化Navier-Stokes 方程并求解，得到了圓盤表面平均膜厚關(guān)系式。Zhevalkink的結(jié)果在低轉(zhuǎn)速（6～33 r/min）下與實驗吻合的較好，因而得到了較廣泛的應(yīng)用，但是當(dāng)轉(zhuǎn)速較高時誤差較大。Vijayraghvan 等[9]利用膜厚測量儀研究了不同轉(zhuǎn)速、不同物系下圓盤表面的膜厚分布，指出圓盤上液膜厚度的不均勻性，膜厚隨著半徑的增加而增加，實驗還發(fā)現(xiàn)各點膜厚隨轉(zhuǎn)速的增加先增加再減小。Afanasiev 等[10]針對圓盤反應(yīng)器從槽內(nèi)帶上液體的三維現(xiàn)象，通過有限元方法求解四階非線性偏微分方程得到了圓盤表面液膜的分布規(guī)律，建立了二維模型，這種模型一定程度上反映了圓盤表面液膜厚度的變化規(guī)律。楊曉宇等[11]采用流體體積函數(shù)模型，對垂直旋轉(zhuǎn)圓盤反應(yīng)器的流動和成膜特性進(jìn)行數(shù)值模擬，并考察了圓盤轉(zhuǎn)速和反應(yīng)器內(nèi)液位高度對膜厚分布規(guī)律的影響。

關(guān)于圓盤反應(yīng)器液相流型，有兩個經(jīng)典模型：Batchelor[12]流型與Stewartson[13]流型。王良生等[14]通過實驗研究發(fā)現(xiàn)，在一定的粘度和轉(zhuǎn)速范圍內(nèi)，盤間距對流型有顯著的影響。Fujie 等[15]通過類比生物轉(zhuǎn)盤與連續(xù)攪拌釜得到了圓盤反應(yīng)器功率的相關(guān)影響因素，并將功率數(shù)與Re數(shù)相關(guān)聯(lián)。

本研究嘗試用計算流體力學(xué)方法，通過流體體積函數(shù)（VOF）模型考察流場，揭示流型、速度分布特征，在此基礎(chǔ)上分析宏觀特性如軸功率和平均膜厚，探討盤間距、盤與壁面距離等因素的影響。

1 計算方法和模型驗證

根據(jù)圓盤反應(yīng)器內(nèi)氣液兩相的流動特點和接觸方式，采用Hirt 等[16]提出的VOF 模型進(jìn)行計算。計算中將圓盤和軸表面的邊界設(shè)置為旋轉(zhuǎn)壁面，壁面邊界條件為無滑移，自由面通過求解VOF 模型方程得到。采用Gambit 進(jìn)行網(wǎng)格劃分，ANSYS FLUENT 商業(yè)軟件進(jìn)行求解。圓盤反應(yīng)器網(wǎng)格劃分和計算區(qū)域均為六面體網(wǎng)格，為了減少網(wǎng)格數(shù)，節(jié)約計算時間，根據(jù)圓盤兩側(cè)流場關(guān)于圓盤對稱的特點，只計算圓盤一側(cè)的流場區(qū)域，另一側(cè)流場采用對稱邊界條件獲得。具體計算方法和模型驗證參見作者另一篇關(guān)于圓盤反應(yīng)器表面更新數(shù)值模擬的研究[17]。

2 結(jié)果與討論

2.1 流場分析

2.1.1 圓盤表面液膜速度分布

對半徑R為150 mm 的圓盤，μ為1 000 mPa·s，圓盤轉(zhuǎn)速N為90 r/min 進(jìn)行數(shù)值模擬，研究圓盤表面液膜的速度分布。以網(wǎng)格內(nèi)液相體積分?jǐn)?shù)為0.99 的等值面作為液膜自由面，即氣液接觸面。速度矢量分布如圖1所示，根據(jù)速度差異可以將圓盤表面分為三個區(qū)：起始區(qū)，加速區(qū)，勻速區(qū)。液體從槽體內(nèi)的靜止?fàn)顟B(tài)到隨圓盤轉(zhuǎn)動必然存在著加速過程。在膜厚方向（z方向）壁面邊界無滑移，圓盤轉(zhuǎn)動帶動內(nèi)層液體運動；由于粘性力作用，內(nèi)層液體通過動量交換使外層流體加速，形成加速區(qū)；隨著外層流體速度增加，z方向速度梯度減小，流體內(nèi)摩擦力減小，加速度減小，最后勻速，形成勻速區(qū)。（在本研究所有的流場圖中，x方向表示水平方向，-y方向表示重力方向，z方向表示膜厚方向。）

圖1 表面液膜速度矢量圖譜Fig.1 Liquid film velocity vector of the surface

圖2 起始區(qū)和加速區(qū)速度矢量圖譜Fig.2 The velocity vector of initial region and accelerating region

類似于平板抽出模型，在起始區(qū)，速度矢量如圖2中的下部分，液膜較厚，膜厚方向速度梯度小，部分液體有向下的速度回流到槽里，沒有隨圓盤轉(zhuǎn)動。如圖3，起始階段速度先減小再增加也是因為部分液體回流，（以VL/VD比值作為縱坐標(biāo)，VL為液膜表面速度，VD為r為150 mm 處圓盤線速度）。

在加速區(qū)，如圖2上部分，速度隨著角度增加而增加，根據(jù)質(zhì)量守恒，液膜加速，里面的液體暴露出來，表面得以更新。如圖3，當(dāng)θ大于π/2 后，加速度減小，進(jìn)入勻速區(qū)域，勻速區(qū)的速度變化較小，液膜速度與半徑成正比，與角度無關(guān)。

2.1.2 圓盤表面液膜厚度分布

膜厚分布如圖4，分別取圓盤表面半徑r為50，75，100 mm 局部各點的膜厚，如圖5，膜厚在切向隨角度的增加而單調(diào)減小，表明圓盤抽出時液膜彎月面的存在，與平板抽出模型膜厚隨平板高度分布規(guī)律一致。另外，膜厚隨著半徑的增加而增加，當(dāng)圓盤從液體中抽出時，盤上各點的切線速度與半徑成正比，速度越大，圓盤帶起的液體越多。起始區(qū)和加速區(qū)膜厚隨半徑增加而明顯增加，在勻速區(qū)膜厚與半徑和角度關(guān)系較小，基本趨于恒定。

圖3 液膜表面速度大小徑向和切向分布規(guī)律Fig.3 The radial and tangential distribution of liquid film velocity

圖4 圓盤表面膜厚分布Fig.4 The distribution of film thickness

圖5 膜厚徑向和切向分布規(guī)律Fig.5 The radial and tangential distribution of film thickness

2.1.3 圓盤間距對槽內(nèi)液相流型的影響

模擬了μ為1 000 mPa·s，N為60 r/min，盤間距L分別為25，35，45，60，80 mm，半徑R為75 mm 的兩個圓盤之間的流場。

模擬發(fā)現(xiàn)，當(dāng)L為80 mm，槽內(nèi)流體出現(xiàn)Stewartson 流型(A 流型)，如圖6所示。Stewartson認(rèn)為盤兩側(cè)均為離心機(jī)，流體只在邊界層內(nèi)隨盤面旋轉(zhuǎn)，而軸向流動和徑向流動均發(fā)生在邊界層外。從俯視圖可以看出，兩盤之間的液相區(qū)域形成兩個渦，渦中心沿豎直方向，稱為“垂直渦”，渦中間存在明顯的“死區(qū)”，流體只在其中做局部循環(huán)。圖中紅色代表液相，藍(lán)色代表氣相，如圖7，在z為10 mm 平面處的速度矢量圖上氣液界面水平，兩相均存在中心沿水平方向的渦，稱為“水平渦”，“水平渦”在液相內(nèi)形成一個總體循環(huán)，更新液相與氣液的接觸面。這兩個漩渦主導(dǎo)槽內(nèi)流體的運動。因此可以推斷在靠近圓盤區(qū)域10～20 mm 內(nèi)液相存在著三維的循環(huán)。

圖6 L 為80 mm 速度矢量俯視圖譜 (y=-10 mm)Fig.6 The top view of velocity vector (y=-10 mm）

圖7 L 為80 mm 速度矢量側(cè)視圖譜 (z=10 mm)Fig.7 The side view of velocity vector (z=10 mm)

圖8 L 為45 mm 時速度矢量圖譜 (y=-10 mm)Fig.8 The top view of velocity vector (y=-10 mm)

當(dāng)L為45 mm，槽內(nèi)流體流型如圖8所示(B 流型)，隨著圓盤間距的減小，“垂直渦”被“擠壓”，漩渦的中心分別向圓盤靠攏，這是個過渡區(qū)域，“垂直渦”依然存在，但是“死區(qū)”消失，最有利于混合。

當(dāng)L為25 mm，隨著圓盤間距進(jìn)一步減小出現(xiàn)如圖9所示C 流型，漩渦中心消失，“垂直渦”消失，從z為10 mm 的速度矢量圖上可以看出，“水平漩渦”依然存在。如圖10，氣液界面開始傾斜，隨著盤間距的減小，部分液體被夾帶到氣相，部分氣體被夾帶到液相?？梢酝茢啵P間距過小，會使兩個盤之間液體相連，兩個盤合并成一個盤，盤間液體隨圓盤一起轉(zhuǎn)動，無法形成有效地更新。

圖9 L 為25 mm 時速度矢量圖譜 (y=-10 mm)Fig.9 The top view of velocity vector (y=-10 mm)

圖10 L 為25 mm 時 速度矢量圖譜 (z=10mm)Fig.10 The side view of velocity vector (z=10mm)

總體而言，以上流型都屬于Stewartson 流型，只是在不同條件下有一定差別。根據(jù)各種流型隨盤間距的演變過程分析，認(rèn)為B 流型可能是最有利于槽體內(nèi)液相混合。

2.1.4 圓盤與槽體之間的物質(zhì)交換

Yamane等[18]在1972年最早研究了生物轉(zhuǎn)盤的圓盤表面與槽體之間的物質(zhì)交換。假設(shè)圓盤轉(zhuǎn)動時，盤上液體與釜內(nèi)液體完全混合，即認(rèn)為交換更新率為100%，并得到了圓盤的傳質(zhì)速率表達(dá)式。Bintanja等[7]對這個假設(shè)進(jìn)行了檢驗發(fā)現(xiàn)，這一假設(shè)未必成立。Zeevalkink 等[19]研究發(fā)現(xiàn)，一定轉(zhuǎn)速范圍內(nèi)實驗測得的傳質(zhì)系數(shù)比理論值小50%以上。認(rèn)為圓盤每旋轉(zhuǎn)一周，僅有25%～50%的液膜被更新，遠(yuǎn)遠(yuǎn)低于Yamane 等假設(shè)的100%。

本研究先通過三維非穩(wěn)態(tài)VOF 兩相流模擬半徑均為75 mm 的兩個圓盤之間的流場。液體粘度μ為1 000 mPa·s，N為60 r/min，L分別為25，35，45，60，80 mm。流場穩(wěn)定后，在槽體內(nèi)均勻加入少量具有同種性質(zhì)的第三類互溶組分，待圓盤旋轉(zhuǎn)一周后統(tǒng)計圓盤表面第三類組分含量，根據(jù)其比例獲得液膜旋轉(zhuǎn)一周后被更新的百分率。不同間距圓盤液膜更新率如圖11所示。

由圖可知，在CFD 模擬的圓盤間距與圓盤直徑之比（L/D）范圍內(nèi)，液膜更新率也只有25%～45%，與Zeevalkink 等的研究結(jié)果一致，遠(yuǎn)小于100%。當(dāng)L/D為0.3 時圓盤表面液膜獲得最大的更新百分率。此時處于B 流型，這也與流場分析結(jié)果一致。

圖11 圓盤液膜更新率與間距關(guān)系Fig.11 The relationship between liquid renewal and distance of discs

2.2 性能預(yù)測

在圓盤反應(yīng)器的不同應(yīng)用中，功率特性和持液量都是其基本性能參數(shù)，通過CFD 對其進(jìn)行研究，給出了其估算方法，并探討了相關(guān)影響因素。功率與力矩和轉(zhuǎn)速的關(guān)系為：P=Mω

已知力矩即可得知功率，因而獲得功率的問題轉(zhuǎn)化為CFD 計算力矩。對液面以上網(wǎng)格進(jìn)行液相體積積分，能夠直接獲得圓盤持液量，根據(jù)圓盤面積轉(zhuǎn)換成相應(yīng)的平均膜厚。

2.2.1 圓盤間距對性能的影響

以L/D為橫坐標(biāo)，膜厚和功率隨盤間距變化如圖12和圖13所示。由圖可知，膜厚和功率都隨著圓盤間距的增加而減小，當(dāng)L/D大于0.4 時，膜厚和功率基本不變。

圖12 平均膜厚隨盤間距的變化規(guī)律Fig.12 The law of average film thickness changes with disc spacing

圖13 功率隨圓盤間距的變化規(guī)律Fig.13 The law of power changes with disc spacing

2.2.2 槽體壁效應(yīng)對性能的影響

改變圓盤與槽體間距d，進(jìn)行模擬計算，圓盤半徑為75 mm。壁效應(yīng)是筒體對圓盤性能的影響，與圓盤和筒體間距有關(guān)，也與圓盤性能有關(guān)，膜厚是性能的重要體現(xiàn)，當(dāng)間距一定時，膜厚越厚，壁效應(yīng)也會越明顯。因此以d/h為橫坐標(biāo)無量綱化。模擬結(jié)果以間距d/h為橫坐標(biāo)，示于圖14和圖15。由圖可得，膜厚隨著圓盤與筒體間距的增加而增加，功率則減小。當(dāng)圓盤與筒體間距d/h大于10 時，圓筒的壁效應(yīng)可以忽略?？梢姳谛?yīng)對于膜厚和功率的影響是同步的，但效應(yīng)相反。

圖14 壁效應(yīng)對膜厚的影響Fig.14 The influence of wall effect on film thickness

圖15 壁效應(yīng)對功率的影響Fig.15 The influence of wall effect on power

2.2.3 圓盤功率特性

根據(jù)CFD 計算得到的功率進(jìn)行數(shù)據(jù)關(guān)聯(lián)，線性擬合得到：P∝R3.19，P∝ω1.97。假設(shè)：NP=βRem(Fr/Re)n，通過P與R和ω量綱分析得到：m=-0.95，n=-0.09。線性擬合如圖16所示。因此，NP=6.81Re-0.95(Fr/Re)-0.09。其中：0

圖16 功率數(shù)線性擬合Fig.16 Linear fitting of power number

圖17 平均膜厚擬合Fig.17 Linear fitting of average film thickness

2.2.4 圓盤持液量特性

根據(jù)CFD 計算得到的平均膜厚進(jìn)行數(shù)據(jù)關(guān)聯(lián)，線性擬合得到：h∝R0.441，h∝ω0.731，h∝μ0.517。假設(shè)h/R為kFraReb，其中k為比例系數(shù)，通過量綱分析得到：a=0.662，b=-0.610。將h/R與Fr0.662Re-0.610線性擬合，如圖17所示。因此k為1，即：h/R=Fr0.662Re-0.610。其中：0

圖18 不同文獻(xiàn)平均膜厚對比Fig.18 The comparison of average film thickness with different literature

與文獻(xiàn)報道的關(guān)聯(lián)式進(jìn)行對比，以CFD 計算得到的h/R值作為橫坐標(biāo)，以相應(yīng)條件下通過文獻(xiàn)公式計算得到的h/R值作為縱坐標(biāo)作圖，結(jié)果如圖18所示。由圖可知，通過CFD 計算得到的膜厚剛好分布在對角線上，不同文獻(xiàn)得到的數(shù)值與其比較。1942年Levich 平板抽出模型得到的膜厚最大，在對角線上方，由于沒有考慮重力因素，膜厚明顯增加，另外平板抽出模型與圓盤旋轉(zhuǎn)的流體受力和流動行為有一定差異。在轉(zhuǎn)速較低、膜厚較薄時，本研究得到的關(guān)聯(lián)式與Zhevalkink[8]的關(guān)聯(lián)式計算結(jié)果是一致的，但隨著轉(zhuǎn)速和粘度增加，膜厚變厚，兩者的差異增大，主要原因是Zhevalkink 實驗條件下Fr小于0.4，忽略了離心力和慣性力的作用，隨著轉(zhuǎn)速和粘度增加，這些因素的影響也趨于明顯。

3 結(jié) 論

運用CFD 方法，通過VOF 模型研究圓盤反應(yīng)器流場特性和主要性能，所得結(jié)論如下：

a）圓盤表面分為起始區(qū)，加速區(qū)和勻速區(qū)三個部分。起始區(qū)內(nèi)液膜較厚，不同位置膜厚變化大，液膜速度較低，部分液體回流，存在速度為零的駐點。加速區(qū)膜厚變化大，液膜加速，內(nèi)層的液體暴露出來，更新較快。勻速區(qū)液膜厚度最薄，速度最大，但是兩者變化較小。

b）圓盤反應(yīng)器釜內(nèi)液相主體呈現(xiàn)Stewartson 流型，由“垂直渦”和“水平渦”主導(dǎo)，受圓盤間距影響明顯，隨著盤間距的減小“垂直渦”相互擠壓，“死區(qū)”消失，當(dāng)盤間距進(jìn)一步減小時，液相和氣相夾帶嚴(yán)重，不能正常操作?！按怪睖u”是槽體與圓盤進(jìn)行物質(zhì)交換的主要方式。計算表明，圓盤適宜間距：0.2

c）通過量綱分析和數(shù)值模擬得到圓盤反應(yīng)器的功率數(shù)：

d）平均膜厚表達(dá)式：

適用于更廣泛的操作范圍和物系范圍。

符號說明

Ca——Capillary 數(shù)，Ca=μωR/σ

D——圓盤直徑，mm

Fr——Froude 數(shù)，F(xiàn)r=ω2R/g

h——液膜厚度，mm

h∞——恒定膜厚，mm

H——液位到圓盤中心的距離，mm

L——圓盤間距，mm

r——某點位置半徑，mm

R——圓盤半徑，mm

Re——雷諾數(shù)，Re=ρω2R/μ

T——無量綱膜厚

u——速度矢量，m/s

α——相體積分率

η——液相更新率

μ——液體粘度，mPa·s

ω——圓盤轉(zhuǎn)速，rad/s

[1]Li K, Yang C, Xu Y, et al.Effect of inorganic anions on Rhodamine B removal under visible light irradiation using Bi2O3/Ti rotating disk reactor [J].Chemical Engineering Journal, 2012, 211:208-215.

[2]Sakurai A, Imai H, Takenaka Y, et al.Simulation of citric acid production by rotating disk contactor [J].Biotechnology and Bioengineering, 1997, 56(6):689-696.

[3]Vasiliadou I A, Molina R, Martínez F, et al.Experimental and modeling study on removal of pharmaceutically active compounds in rotating biological contactors [J].Journal of Hazardous Materials, 2014, 274:473-482.

[4]Landau L, Levich B.Dragging of a liquid by a moving plate [J].Acta Physico-Chemica USSR, 1942, 17:42-54.

[5]White D A, Tallmadge J A.Theory of drag out of liquids on flat plates [J].Chemical Engineering Science, 1965, 20(1):33-37.

[6]Groenveld P.Dip-coating by withdrawal of liquid films [D].Delft:Delft University, 1970.

[7]Bintanja H H J, van der Erve J J V M, Boelhouwer C.Oxygen transfer in a rotating disc treatment plant [J].Water Research, 1975,9(12):1147-1153.

[8]Zhevalkink J A, Kelderman P, Boelhouwer C.Liquid film thickness in a rotating disc gas-liquid contactor [J].Water Research, 1978,12(8):577-581.

[9]Vijayraghvan K, Gupta J P.Thickness of the film formed on a vertically rotating disk partially immersed in a Newtonian liquid [J].Industrial &Engineering Chemistry Fundamentals, 1982, 21(4):333-336.

[10]Afanasiev K, Münch A, Wagner B.Thin film dynamics on a vertically rotating disk partially immersed in a liquid bath [J].Applied Mathematical Modelling, 2008, 32(9):1894-1911.

[11]楊曉宇, 李志鵬, 高正明.圓盤反應(yīng)器內(nèi)流動和成膜性能的數(shù)值模擬 [J].北京化工大學(xué)學(xué)報:自然科學(xué)版, 2010, 37(6):29-34.Yang Xiaoyu, Li Zhipeng, Gao Zhengming.Numerical simulation of the flow and film-forming characteristics in a horizontal rotating disk reactor [J].Journal of Beijing University of Chemical Technology, 2010, 37(6):29-34.

[12]Batchelor G K.Note on a class of solutions of the Navier-Stokes equations representing steady rotationally-symmetric flow [J].The Quarterly Journal of Mechanics and Applied Mathematics, 1951, 4(1):29-41.

[13]Stewartson K.On the flow between two rotating coaxial disks[C]//Proc.Camb.Phil.Soc.1953, 49(2):333-341.

[14]王良生, 周賢爵.圓盤反應(yīng)器的混合特性 [J].高?；瘜W(xué)工程學(xué)報, 2000, 14(1):37-42.Wang Liangsheng, Zhou Xianjue.The research of mixing characteristic in disc-ring reactor [J].Journal of Chemical Engineering of Chinese Universities, 2000, 14(1):37-42.

[15]Fujie K, Bravo H E, Kubota H.Operational design and power economy of a rotating biological contactor [J].Water Research, 1983,17(9):1153-1162.

[16]Hirt C W, Nichols B D.Volume of fluid (VOF) method for the dynamics of free boundaries [J].Journal of Computational Physics, 1981,39(1):201-225.

[17]鄧 斌, 戴干策.圓盤反應(yīng)器液膜表面更新數(shù)值模擬 [J].化工學(xué)報, 2015, 66(4):1407-1416.Deng Bin, Dai Gance.Numerical simulation of surface renewal frequency on vertically rotating disc [J].Journal of Chemical Industry and Engineering, 2015, 66(4):1407-1416.

[18]Yamane T, Yoshida F.Absorption in a rotating-disk gas-liquid contactor [J].Journal of Chemical Engineering of Japan, 1972,5(4):381-385.

[19]Zeevalkink J A, Kelderman P, Visser D C, et al.Physical mass transfer in a rotating disc gas-liquid contactor [J].Water Research, 1979,13(9):913-919.