莊黎偉，戴干策

?

中空纖維膜外壓全量過(guò)濾動(dòng)態(tài)過(guò)程的數(shù)值模擬

莊黎偉，戴干策

（華東理工大學(xué)化學(xué)工程聯(lián)合國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200237）

摘要:建立了中空纖維膜外壓式全量過(guò)濾的CFD模型，模擬膜絲長(zhǎng)度、直徑、滲透系數(shù)、裝填密度、污染指數(shù)以及跨膜壓差不同條件，得到通量分布和產(chǎn)水量的動(dòng)態(tài)演變過(guò)程。研究結(jié)果表明:通量分布會(huì)隨著過(guò)濾的進(jìn)行而逐漸變得均勻，這種通量分布的自我調(diào)節(jié)作用在膜絲較長(zhǎng)、較細(xì)，滲透性較好，裝填密度較高，污染指數(shù)較高以及跨膜壓差較高時(shí)更為明顯；產(chǎn)水流量的倒數(shù)與累積產(chǎn)水量呈線性關(guān)系，但由于通量分布不均勻并且均勻性演變，這種線性關(guān)系區(qū)別于傳統(tǒng)濾餅過(guò)濾模型；通過(guò)數(shù)據(jù)擬合得到了適用于中空纖維外壓式全量操作的濾餅過(guò)濾關(guān)聯(lián)式，可用于預(yù)測(cè)組件的性能和指導(dǎo)組件的設(shè)計(jì)。

關(guān)鍵詞:中空纖維膜；計(jì)算流體力學(xué)；模型；數(shù)值模擬；通量分布；動(dòng)態(tài)演變

引 言

中空纖維超濾膜組件憑借其自支撐、高裝填密度、高比表面積[1]的優(yōu)點(diǎn)，被廣泛用于污水處理[2]和海水脫鹽[3]等領(lǐng)域。然而，中空纖維膜組件的特殊結(jié)構(gòu)會(huì)造成膜表面通量分布不均[4-7]，從而減小膜面積利用率，縮短組件壽命。因此，有諸多學(xué)者進(jìn)行了中空纖維膜組件內(nèi)通量分布的研究。

早期研究主要采用的是模型外加輔助實(shí)驗(yàn)的方法:Chang等[8]通過(guò)建立浸沒(méi)式中空纖維膜的滲透模型，發(fā)現(xiàn)通量沿膜絲長(zhǎng)度的分布取決于膜絲內(nèi)徑、長(zhǎng)度以及滲透性；在此基礎(chǔ)上，他們先后研究了操作通量低于或高于臨界通量下浸沒(méi)式中空纖維膜組件的通量分布[9]以及全量操作下通量分布在過(guò)濾過(guò)程中的演變情況[10]。隨著計(jì)算流體力學(xué)方法的興起，越來(lái)越多膜研究者將其用于膜組件的設(shè)計(jì)和過(guò)程的優(yōu)化[11]:Günther等以Happel的自由表面模型[12]為依據(jù)，將整個(gè)膜組件簡(jiǎn)化成由單根膜絲和包裹其外圍的環(huán)狀流體層，在此基礎(chǔ)上分別研究了純水過(guò)濾時(shí)裝填密度對(duì)通量分布的影響[13]以及耦合濾餅過(guò)濾模型時(shí)裝填密度對(duì)通量和濾餅分布的影響[14]。由此可知，中空纖維膜組件內(nèi)通量分布不僅與膜絲本身性質(zhì)有關(guān)，還會(huì)受膜絲裝填以及污垢阻力的影響。由于中空纖維膜組件管程和殼程存在質(zhì)量交換，交換過(guò)程由膜多孔介質(zhì)層的流體滲透實(shí)現(xiàn)，滲透通量又與污垢阻力相互耦合，所以中空纖維膜組件內(nèi)局部通量既存在空間上的分布，又隨時(shí)間變化。中空纖維膜組件內(nèi)通量分布的動(dòng)態(tài)演變必然影響組件的產(chǎn)水性能，因而研究這種動(dòng)態(tài)過(guò)程有助于深化理解中空纖維膜組件的過(guò)濾機(jī)理。雖然文獻(xiàn)中存在中空纖維膜組件內(nèi)通量分布動(dòng)態(tài)演變的研究，但是缺乏對(duì)幾何和操作參數(shù)影響的系統(tǒng)考量，從而無(wú)法建立通量分布動(dòng)態(tài)演變與組件過(guò)濾性能的關(guān)系，所以本文的研究目的便在于此。

在前期的研究[4]中，分別采用解析模型和CFD模擬研究了中空纖維膜組件在外壓式全量過(guò)濾時(shí)幾何因素對(duì)通量分布以及能耗利用率的影響。本文在此基礎(chǔ)上，考慮了濾餅阻力和局部通量的耦合，建立了新的中空纖維膜過(guò)濾CFD模型，描述在不同幾何參數(shù)和操作條件下外壓式全量過(guò)濾通量分布和產(chǎn)水流量的動(dòng)態(tài)過(guò)程，并根據(jù)CFD模擬結(jié)果擬合出產(chǎn)水流量與累積產(chǎn)水量的關(guān)聯(lián)式。

圖1　中空纖維膜組件Fig.1 Hollow fiber membrane module

1　模型建立

1.1 物理模型

計(jì)算所用物理模型結(jié)構(gòu)如圖1(a)所示。4根長(zhǎng)度為0.5 m，內(nèi)外徑分別為0.7和1.3 mm，滲透系數(shù)為4.9×10-16m2的中空纖維膜絲均勻裝填于長(zhǎng)為0.5 m、直徑為6 mm的圓柱形殼體內(nèi)，裝填密度約為18.8%。原水由入口進(jìn)入到殼程，在沿著軸向流動(dòng)的同時(shí)，不斷通過(guò)膜表面滲透至管程，由于是全量操作，原水全部轉(zhuǎn)化為滲透液，并從管程出口離開。污垢全部截留在膜外表面，在滲透階段，滲透阻力不斷增大，直至進(jìn)入洗滌階段?？缒翰顬?.1 MPa，恒壓操作。原水濃度為0.2 g·L-1，濾餅比阻為3×1014m·kg-1，忽略濾餅的可壓縮性。在以上結(jié)構(gòu)參數(shù)和操作參數(shù)的基礎(chǔ)上，改變其中某一參數(shù)值，研究不同參數(shù)下過(guò)濾的動(dòng)態(tài)特征。共考慮6個(gè)參數(shù)，分別為幾何參數(shù)（膜絲長(zhǎng)度、直徑、滲透系數(shù)、裝填密度）和操作參數(shù)（跨膜壓差、污染指數(shù)），每種參數(shù)取7種情況，所有參數(shù)值列于表1。其中污染指數(shù)FI（fouling index）為原水濃度和濾餅比阻的乘積。在改變直徑時(shí)，同時(shí)改變滲透系數(shù)，以保證不同直徑膜絲的點(diǎn)跨膜阻力均為6.12×1011m-1，消除滲透系數(shù)的影響。在研究裝填密度時(shí)，采用文獻(xiàn)[4,13-14]的方法，保持膜絲尺寸不變，通過(guò)改變膜絲外圍環(huán)形殼程空間大小來(lái)調(diào)節(jié)裝填密度。

1.2 控制方程

在管程與殼程可認(rèn)為是非定態(tài)層流不可壓縮流動(dòng)，因而連續(xù)性方程和動(dòng)量方程分別為

式中，v為流體速度；ρ為流體密度；p為壓強(qiáng)；μ為流體運(yùn)動(dòng)黏度。對(duì)于膜多孔介質(zhì)區(qū)域，存在膜阻力和污垢阻力，對(duì)動(dòng)量方程進(jìn)行修正，得到多孔介質(zhì)區(qū)域的連續(xù)性方程和動(dòng)量方程

式中，λ為膜絲實(shí)體孔隙度，設(shè)定為 0.7；Km為多孔介質(zhì)的滲透系數(shù)，且隨著過(guò)濾的進(jìn)行而不斷降低。本文將污垢阻力等效成膜多孔介質(zhì)阻力的增幅，該方法簡(jiǎn)化了污垢層的計(jì)算，推導(dǎo)過(guò)程如下。

根據(jù)阻力串聯(lián)定律[14-15]，任何一點(diǎn)通量J可以表示為

式中，Rm和 Rf分別為膜阻力和污垢阻力。對(duì)于污垢阻力，主要來(lái)源于污垢造成孔道阻塞、孔道收縮以及外層形成的濾餅。本文假設(shè)污垢顆粒大于膜孔徑，所以污垢阻力僅為濾餅阻力。因此，污垢阻力Rf可以表示為[16]

式中，m為單位面積濾餅的質(zhì)量；α為濾餅比阻。任意一點(diǎn)m隨時(shí)間的變化關(guān)系可以表示為[16]

式中，cs為原水污垢質(zhì)量濃度；φf(shuō)和φs為濾餅層和原水中固體體積分?jǐn)?shù)。由于全量過(guò)濾多用于低濃度污水處理，φf(shuō)?φs，所以有式（7）的簡(jiǎn)化。將微分形式的式（7）進(jìn)行時(shí)間離散得到式（7）的差分形式，同時(shí)聯(lián)立式（6），可以得到濾餅阻力在時(shí)間步長(zhǎng)Δt內(nèi)的增幅正比于t～(t+Δt)時(shí)段內(nèi)該處的局部通量（由于時(shí)間步長(zhǎng)很小，可認(rèn)為該時(shí)段通量不發(fā)生變化），因而可以得到任意時(shí)刻任意點(diǎn)濾餅阻力與該處局部通量所有歷史值之間的關(guān)系，最終任意點(diǎn)任意時(shí)刻通量可以表示為

式中，Δp為局部跨膜壓差。此外，式（4）中膜的滲透系數(shù)可以表示為

式中，δ為膜多孔區(qū)域厚度。通過(guò)式（8）實(shí)現(xiàn)點(diǎn)通量和滲透阻力的耦合。理論上只要時(shí)間步長(zhǎng)夠小，聯(lián)立式（1）～式（4）和式（10）就可以描述整個(gè)系統(tǒng)滲透的動(dòng)態(tài)過(guò)程。

上述模型還基于以下假設(shè):全量過(guò)濾時(shí)軸向速度較小，忽略平行于膜表面流動(dòng)剪切對(duì)污垢的軸向輸運(yùn)作用，即污垢一旦與膜或?yàn)V餅接觸便黏附；污垢與流體密度相近，跟隨性好，忽略布朗擴(kuò)散和慣性作用[14]；忽略濾餅層對(duì)流動(dòng)通道尺寸的影響，該假設(shè)在裝填密度較低或者滲透周期較小時(shí)可認(rèn)為合理。

表1　幾何參數(shù)和操作參數(shù)Table 1 Geometric and operating parameters

1.3 邊界條件和初始條件

如圖1(a)所示，原水進(jìn)口和滲透液出口分別設(shè)為壓強(qiáng)入口和壓強(qiáng)出口，其中出口處壓強(qiáng)保持 0 MPa，進(jìn)口處壓強(qiáng)根據(jù)表1的跨膜壓差而進(jìn)行調(diào)整；計(jì)算域外圍其他區(qū)域均設(shè)為無(wú)滑移壁面；計(jì)算域內(nèi)部膜實(shí)體區(qū)設(shè)為多孔介質(zhì)，僅考慮黏性阻力，阻力系數(shù)根據(jù)式（10）確定；多孔介質(zhì)區(qū)域與管殼程交界面壓強(qiáng)和速度連續(xù)，該種處理方法能定性描述交界面處滑移速度的軸向分布[4]，優(yōu)于文獻(xiàn)慣用的速度無(wú)滑移處理[13-14]。初始時(shí)，內(nèi)部速度和壓強(qiáng)均設(shè)為0，原水在跨膜壓差驅(qū)動(dòng)下進(jìn)入設(shè)備。

2　數(shù)值求解

2.1 方程離散與求解

計(jì)算域內(nèi)整個(gè)流動(dòng)過(guò)程，通過(guò)基于有限體積法的計(jì)算流體力學(xué)軟件Fluent 6.3.26進(jìn)行離散求解，壓強(qiáng)速度以SIMPLE算法進(jìn)行耦合，壓強(qiáng)離散為二階格式，動(dòng)量離散為二階迎風(fēng)格式，求解過(guò)程中，壓強(qiáng)、動(dòng)量的亞松弛因子分別為 0.3、0.7。時(shí)間步長(zhǎng)設(shè)為 10-3s，每個(gè)時(shí)間步內(nèi)連續(xù)性和速度殘差均低于10-5時(shí)進(jìn)入下一時(shí)間步。先在純水過(guò)濾過(guò)程中計(jì)算至流量穩(wěn)定，1 s時(shí)開啟污水過(guò)濾模式，計(jì)算至101 s，同時(shí)記錄實(shí)時(shí)流量。多孔介質(zhì)黏性阻力與流動(dòng)的耦合通過(guò)自編程序（UDF）實(shí)現(xiàn)。求解過(guò)程通過(guò)兩個(gè)Intel Xeon CPU和64 GB的DDR4型號(hào)內(nèi)存完成。

2.2 網(wǎng)格與時(shí)間步長(zhǎng)無(wú)關(guān)性

整個(gè)計(jì)算域均為結(jié)構(gòu)化網(wǎng)格，以1.1節(jié)給出的特定尺寸的計(jì)算域?yàn)槔ざ嗫捉橘|(zhì)區(qū)軸向均勻排布50個(gè)網(wǎng)格層，與管程和殼程一致，徑向上有兩個(gè)網(wǎng)格層，多孔介質(zhì)區(qū)和殼程網(wǎng)格為六面體，管程均為五面體，如圖1(b)所示。由圖2可知，時(shí)間步長(zhǎng)為10-3s時(shí)，網(wǎng)格數(shù)由48850增加至403100，1～100 s的產(chǎn)水量隨時(shí)間的衰減曲線最大偏差在1%以內(nèi)；同時(shí)，由圖3可知，兩種網(wǎng)格尺寸下1 s和100 s兩個(gè)時(shí)刻點(diǎn)通量的軸向分布最大偏差也在 6%以內(nèi)，可認(rèn)為48850的網(wǎng)格數(shù)或相同網(wǎng)格尺寸下所得計(jì)算結(jié)果的數(shù)值誤差在允許范圍內(nèi)。根據(jù)圖2、圖3可知，時(shí)間步長(zhǎng)由10-3s減小至5×10-4s時(shí)，計(jì)算結(jié)果幾乎無(wú)差別，可認(rèn)為10-3s的時(shí)間步長(zhǎng)足夠捕捉過(guò)濾過(guò)程的動(dòng)態(tài)特征。

圖2　產(chǎn)水量隨時(shí)間的變化Fig.2 Volumetric flow rate of permeate against time

3　結(jié)果與討論

3.1 模型驗(yàn)證

首先，通過(guò)實(shí)驗(yàn)驗(yàn)證CFD模型在模擬純水滲透方面的準(zhǔn)確性，實(shí)驗(yàn)裝置和操作細(xì)節(jié)見(jiàn)文獻(xiàn)[4]。采用 5種長(zhǎng)度中空纖維膜絲（其余膜絲細(xì)節(jié)與 1.1節(jié)一致）得出單根膜絲產(chǎn)水量隨膜絲長(zhǎng)度的變化關(guān)系，如圖4所示，模擬結(jié)果與實(shí)驗(yàn)值相差在10%以內(nèi)。

圖3　點(diǎn)通量的軸向分布Fig.3 Flux distribution in axial direction

圖4　單根膜絲產(chǎn)水速率隨長(zhǎng)度變化的實(shí)驗(yàn)值和模擬值Fig.4 Volumetric flow rate of single fiber with different length based on experiment and simulation

圖5　組件跨膜壓差隨時(shí)間變化的模擬結(jié)果與文獻(xiàn)結(jié)果[16]對(duì)比Fig.5 Comparison of trans-membrane pressure against time in present study and Ref. [16]

然后，參照文獻(xiàn)[16]的中空纖維膜組件尺寸和操作參數(shù)，模擬了恒流量下跨膜壓差隨時(shí)間的變化關(guān)系，并與文獻(xiàn)的模擬結(jié)果進(jìn)行對(duì)比（圖5），最大偏差在8%以內(nèi)。需要指出的是，文獻(xiàn)[16]考慮了濾餅厚度對(duì)殼程通道的影響，但是裝填密度較低，可認(rèn)為通道尺寸未發(fā)生明顯變化，所以與本文的模擬細(xì)節(jié)一致。

因此，實(shí)驗(yàn)結(jié)果和文獻(xiàn)結(jié)果驗(yàn)證了本文CFD模型的合理性。

3.2 壓強(qiáng)場(chǎng)與速度場(chǎng)

圖6、圖7分別給出了組件內(nèi)部的壓強(qiáng)和速度分布，該圖基于1.1節(jié)中幾何和操作參數(shù)，過(guò)濾時(shí)間為101 s。由圖6可知，殼程壓強(qiáng)遠(yuǎn)高于管程，這是因?yàn)橥鈮菏竭^(guò)濾需要跨膜壓差作為推動(dòng)力；由圖7可知，殼程速度逐漸降低，而管程速度逐漸升高，這由全量過(guò)濾殼程原水不斷滲透至管程造成。由此可知，模擬所得壓強(qiáng)和速度分布符合外壓全量過(guò)濾流場(chǎng)特征。

圖6　壓強(qiáng)分布Fig.6 Pressure distribution

圖7　速度分布Fig.7 Velocity distribution

3.3 通量分布均勻性演變

圖8(a)～(f)分別給出了不同膜絲長(zhǎng)度、直徑、滲透系數(shù)、裝填密度、污染指數(shù)和跨膜壓差（具體數(shù)值見(jiàn)表1）下，1 s、52 s和101 s等3個(gè)時(shí)刻通量分布的均勻性。6幅圖中框出的數(shù)據(jù)均為基于 1.1節(jié)中幾何參數(shù)和操作參數(shù)下模擬所得結(jié)果，除了圖8（d）中裝填密度為9.1%，低于其他圖中的18.8%。通量分布均勻性采用軸向上量綱1局部通量的標(biāo)準(zhǔn)差SD表示

式中，J為平均通量。SD越高，通量分布均勻性越差。

由圖8可知，隨著過(guò)濾的進(jìn)行，通量分布的均勻性有所提高，這種自我調(diào)節(jié)作用在許多文獻(xiàn)[9,17]中都有提及。原因在于，初始局部通量較高的位置，濾餅形成速率高于初始局部通量較低區(qū)域，而濾餅形成速率對(duì)應(yīng)著過(guò)濾阻力增加的速率，從而在一定的點(diǎn)跨膜壓差分布下，局部通量高的區(qū)域通量降低速率更快，通量分布逐漸均勻。

對(duì)于t=1 s，即初始通量分布，由圖8(a)、(b)可知，隨著膜絲長(zhǎng)度的增加和膜絲直徑的減小，通量分布的初始均勻性逐漸降低，這是因?yàn)楣艹虊航翟斐闪司植靠缒翰钛啬そz長(zhǎng)度逐漸增大，使靠近出口端的局部通量最高，靠近入口端的最低，這種分布從圖3可以看出，并且這種不均勻性在較長(zhǎng)[4]和較細(xì)[8]的膜絲中尤為明顯。由圖 8(c)可知，膜滲透系數(shù)越低，通量分布的初始均勻性越好，這是由于滲透系數(shù)越低意味著阻力越大，通道阻力的提高會(huì)一定程度削弱造成流動(dòng)不均勻因素的影響[18-20]，從而使流動(dòng)分布更均勻，所以超濾、微濾膜中通量分布均勻性普遍差于納濾以及反滲透。由圖8(d)可知，通量分布的初始均勻性隨著裝填密度的增大先漸漸變好，然后急劇惡化。根據(jù)文獻(xiàn)[4]的理論模型可知，中空纖維膜外壓式全量過(guò)濾時(shí)局部通量沿軸向呈U形分布，當(dāng)裝填密度較?。ǖ陀诩s50%）時(shí)，局部通量最大值靠近出口端，而當(dāng)裝填密度較大時(shí)，膜絲進(jìn)口端的局部通量會(huì)接近甚至超過(guò)出口端數(shù)值，這就解釋了圖8(d)的現(xiàn)象。但是可以看出，當(dāng)裝填密度低于68.1%時(shí)，通量分布的初始均勻性均在 0.01～0.02，差別不大。尤其是當(dāng)裝填密度低于40%時(shí)，可認(rèn)為膜絲相互獨(dú)立，此時(shí)裝填密度的變化對(duì)組件通量分布[4]和產(chǎn)水量[13]沒(méi)有明顯影響。雖然低裝填密度能防止通量分布均勻性的惡化，但是高裝填密度卻能提高膜組件在橫截面的流動(dòng)分布均勻性[19-20]。所以在組件膜絲裝填密度的取舍上，應(yīng)該選擇合適的范圍，例如文獻(xiàn)[20]建議裝填密度在30%左右以兼顧熱質(zhì)傳遞系數(shù)和壓降損耗，而本文的建議范圍是30%～68%以保證通量分布和組件橫截面流動(dòng)分布的均勻性。由圖 8(f)可以看出，跨膜壓差或者說(shuō)是操作通量的增大會(huì)降低通量分布的初始均勻性，這與文獻(xiàn)[7]所得結(jié)果一致，原因是高操作通量增大了管程的壓降。

圖8　幾何和操作參數(shù)對(duì)瞬態(tài)通量分布均勻性的影響Fig.8 Effect of geometrical and operating parameters on uniformity of transient flux distribution

對(duì)于t=52 s和101 s，通過(guò)6幅圖可以發(fā)現(xiàn)，通量分布初始均勻性差的提高幅度最大，比如圖8(a)～(d)中的L7、D7、P1、PD7，或者是最大的污染指數(shù)或最大的跨膜壓差會(huì)最大幅度提升通量分布均勻性，分別如圖8(e)、(f)中的FI7和TMP7，由此說(shuō)明，通量分布自我調(diào)節(jié)作用在初始均勻性較差或者較大污染指數(shù)或者較大跨膜壓差（對(duì)應(yīng)較大操作通量）下才更為明顯。此外，開始的51 s（t=1～52 s）通量分布均勻性提高幅度遠(yuǎn)高于后期的49 s （t=52～101 s）。原因在于，初始通量分布十分不均勻，并且前期通量較高，濾餅生成速率越高，所以自我調(diào)節(jié)作用較為明顯；而后期，由于通量分布均勻性已有明顯提高，并且由于恒壓操作，通量衰減至較小值，所以自我調(diào)節(jié)作用沒(méi)有前期明顯。

3.4 產(chǎn)水量演變

經(jīng)典的濾餅阻力模型可以表示為[5,21]

式中，A為膜表面積；V為累積產(chǎn)水體積。然而，該模型只適用于平面式膜[3]，并且，從模型表達(dá)式也可以看出，該模型只能描述局部濾餅形成過(guò)程或者通量和濾餅阻力分布完全均勻的情況，無(wú)法應(yīng)用于本文的中空纖維膜系統(tǒng)。所以，本文將給出中空纖維膜系統(tǒng)的濾餅阻力模型，適用于存在通量和濾餅阻力分布及其均勻性演變的過(guò)程。

圖9　不同幾何和操作參數(shù)下dt/dV隨V的變化Fig.9 dt/dV versus V under various geometrical and operating conditions

參照式（12），圖9(a)～(f)分別給出了不同膜絲長(zhǎng)度、直徑、滲透系數(shù)、裝填密度、污染指數(shù)和跨膜壓差（具體數(shù)值見(jiàn)表1）條件下，6～101 s（隔5s取一次數(shù)值）中dt/dV隨著滲透液累積體積V的變化。擬合結(jié)果表明，所有曲線都呈高度線性，說(shuō)明通量和濾餅阻力分布及其均勻性演變并未改變dt/dV與 V的線性關(guān)系。所以，中空纖維膜系統(tǒng)內(nèi)濾餅阻力模型可假設(shè)為

式中，a為斜率，b為截距。

圖9說(shuō)明，隨著過(guò)濾的進(jìn)行，產(chǎn)水累積總量越來(lái)越大，但產(chǎn)水流量會(huì)逐漸減小，圖9中任意曲線上相鄰兩點(diǎn)間隔時(shí)間都為5 s，但點(diǎn)的距離卻不斷減小，也說(shuō)明單位時(shí)間的產(chǎn)水量不斷減小。這是因?yàn)楹愣缒翰钕?，濾餅阻力的增大必然使產(chǎn)水流量逐漸減小。

圖9中6幅圖可以分為3類。

第1類是各曲線斜率不同，截距不同，比如圖9(a)、(b)、(f)。截距等于初始或者純水過(guò)濾體積流量的倒數(shù)，反映的是跨膜壓差（推動(dòng)力）和中空纖維膜組件結(jié)構(gòu)（設(shè)備阻力）兩者的作用結(jié)果；而斜率反映了跨膜壓差（推動(dòng)力）、組件結(jié)構(gòu)（設(shè)備阻力）和污垢（過(guò)程阻力）三者的作用結(jié)果。在圖9(a)、(b)中，隨著膜絲增長(zhǎng)和增粗，斜率和截距均不斷減小。參照式（12）可知，膜表面積的增大均會(huì)降低曲線的斜率和截距，而膜絲增長(zhǎng)和增粗正好對(duì)應(yīng)面積的增大，說(shuō)明中空纖維膜過(guò)濾產(chǎn)生濾餅的過(guò)程也有類似的規(guī)律。同樣，圖 9(f)中跨膜壓差的增大也同時(shí)減小了斜率和截距，與式（12）給出的規(guī)律相同。

第2類是各曲線斜率相同，截距不同，比如圖9(c)、(d)。在圖 9(c)中，隨著膜滲透系數(shù)的減小，或者說(shuō)膜阻力的增大，截距不斷增大，初始產(chǎn)水流量不斷減小，也與式（12）規(guī)律相同。在圖9(d)中，隨裝填密度的增大，截距逐漸增大，并且增大速率加快，說(shuō)明過(guò)高裝填度會(huì)顯著降低產(chǎn)水流量，這與文獻(xiàn)[13]所得結(jié)果一致。

第3類是各曲線截距相同，而斜率不同，比如圖9(e)。在該圖中，污染指數(shù)的增大提高了斜率，說(shuō)明進(jìn)料濃度或者濾餅比阻的增大都會(huì)加快產(chǎn)水量流量的衰減。

為了得到式（13）中a和b的表達(dá)式，首先取出圖9中42條曲線的截距，通過(guò)數(shù)據(jù)擬合并參照文獻(xiàn)[16]可知，截距b的表達(dá)式可以在該文獻(xiàn)中純水滲透表達(dá)式的基礎(chǔ)上進(jìn)行簡(jiǎn)單修正而得到；然后，取出圖9中所有曲線的斜率，參照式（12）斜率的表達(dá)式并進(jìn)行修正可以得到斜率a的表達(dá)式。最終

式（13）可以表示為

其中

式中，n為膜絲數(shù)量；Ai和Ae分別為單根膜絲內(nèi)外表面積；di和de分別為內(nèi)外徑；g（ε）=，ε為殼程孔隙度，與裝填密度之和為 1。通過(guò)與曲線（12）斜率比較可以發(fā)現(xiàn)，曲線（14）斜率跟膜的內(nèi)外表面積以及內(nèi)外徑之比都有關(guān)系，說(shuō)明中空纖維膜的圓環(huán)形結(jié)構(gòu)以及管殼程流場(chǎng)都會(huì)使濾餅形成過(guò)程變得復(fù)雜，而無(wú)法采用模型（12）進(jìn)行描述。圖10給出了表1中所有幾何和操作參數(shù)下任意V所對(duì)應(yīng)dt/dV的模擬值與擬合式（14）所得數(shù)值的相對(duì)誤差的平均值，吻合較好，說(shuō)明關(guān)聯(lián)式（14）能描述各種幾何和操作參數(shù)下產(chǎn)水流量和累積產(chǎn)水量的關(guān)系。

圖10　不同幾何和操作參數(shù)下模擬與關(guān)聯(lián)式（14）所得dt/dV的相對(duì)誤差Fig.10 Relative error of dt/dV obtained from simulation to the one based on correlation Eq. (14) under various geometrical and operating conditions

4　結(jié) 論

本文建立了中空纖維膜外壓式全量過(guò)濾 CFD模型，模擬展示了通量分布和產(chǎn)水體積流量的動(dòng)態(tài)演變過(guò)程，并通過(guò)數(shù)據(jù)擬合得到了適用于中空纖維膜過(guò)濾的關(guān)聯(lián)式，結(jié)果表明:

（1）中空纖維膜通量分布會(huì)隨著過(guò)濾的進(jìn)行逐漸變得均勻，這種自我調(diào)節(jié)作用在通量均勻性較差、產(chǎn)水流量較大以及進(jìn)料濃度或?yàn)V餅比阻較大時(shí)更加明顯；

（2）中空纖維膜外壓式全量產(chǎn)水體積流量隨過(guò)濾的進(jìn)行逐漸減小，當(dāng)污垢阻力主要為不可壓縮濾餅阻力時(shí)，產(chǎn)水體積流量的倒數(shù)與累積產(chǎn)水量呈線性關(guān)系，但這種線性關(guān)系是基于通量分布不均勻以及均勻性存在演變的特性，區(qū)別于基于通量分布完全均勻假設(shè)的傳統(tǒng)濾餅阻力模型，因而更能反映中空纖維膜組件過(guò)濾的產(chǎn)水特性；

（3）關(guān)聯(lián)式（14）包含中空纖維膜組件的結(jié)構(gòu)和操作參數(shù)，且適用于較寬的參數(shù)范圍，因而可用于預(yù)測(cè)特定結(jié)構(gòu)下組件的產(chǎn)水性能，并指導(dǎo)組件的優(yōu)化設(shè)計(jì)；同時(shí)，本文的CFD模型在模擬計(jì)算時(shí)耗費(fèi)計(jì)算資源較少，有望用于大型中空纖維膜組件過(guò)濾過(guò)程的模擬，這部分工作將在后期的研究中開展。

符 號(hào) 說(shuō) 明

A ——膜面積，m2

Ae——膜外表面積，m2

Ai——膜內(nèi)表面積，m2

cs——原水污垢質(zhì)量濃度，kg·m-3

de——膜絲外徑，mm

di——膜絲內(nèi)徑，mm

J ——滲透通量，m·s-1

Km——膜滲透系數(shù)，m2

L ——膜絲長(zhǎng)度，m

m ——單位面積濾餅的質(zhì)量，kg·m-2

n ——膜絲數(shù)量

p ——靜壓強(qiáng)，MPa

Δp ——局部跨膜壓差，MPa

QV——體積流量，m3·s-1

Rf——污垢阻力，m-1

Rm——膜阻力，m-1

TMP ——總跨膜壓差，MPa

t ——時(shí)間，s

α ——濾餅比阻，m·kg-1

δ ——膜厚度，mm

ε ——膜絲束殼程空隙度

λ ——膜絲孔隙度

μ ——?jiǎng)恿︷ざ?，Pa·s

ρ——密度，kg·m-3

?f——原水固體體積分?jǐn)?shù)

?s——濾餅層固體體積分?jǐn)?shù)

References

[1] YANG X, WANG R, FANE A G, et al. Membrane module design and dynamic shear-induced techniques to enhance liquid separation by hollow fiber modules: a review [J]. Desalination and Water Treatment, 2013, 51 (16/17/18): 3604-3627.

[2] DELGADO S, DIAZ F, VERA L, et al. Modelling hollow-fibre ultrafiltration of biologically treated wastewater with and without gas sparging [J]. Journal of Membrane Science, 2004, 228 (1): 55-63.

[3] LI N N, FANE A G, HO W W, et al. Advanced Membrane Technology and Applications [M]. New Jersey: John Wiley & Sons, 2011: 47-242.

[4] ZHUANG L W, GUO H F, WANG P H, et al. Study on the flux distribution in a dead-end outside-in hollow fiber membrane module [J]. Journal of Membrane Science, 2015, 495: 372-383.

[5] KIM J, DIGIANO F A. Defining critical flux in submerged membranes: influence of length-distributed flux [J]. Journal of Membrane Science, 2006, 280 (1): 752-761.

[6] LI X H, LI J X, WANG J, et al. Experimental investigation of local flux distribution and fouling behavior in double-end and dead-end submerged hollow fiber membrane modules [J]. Journal of Membrane Science, 2014, 453: 18-26.

[7] LI X H, LI J X, WANG H, et al. A filtration model for prediction of local flux distribution and optimization of submerged hollow fiber membrane module [J]. AIChE Journal, 2015, 61 (12): 4377-4386.

[8] CHANG S, FANE A G. The effect of fibre diameter on filtration and flux distribution—relevance to submerged hollow fibre modules [J]. Journal of Membrane Science, 2001, 184 (2): 221-231.

[9] CHANG S, FANE A G, VIGNESWARAN S. Modeling and optimizing submerged hollow fiber membrane modules [J]. AIChE Journal, 2002, 48 (10): 2203-2212.

[10] CHANG S, FANE A G, WAITE T D. Analysis of constant permeate flow filtration using dead-end hollow fiber membranes [J]. Journal of Membrane Science, 2006, 268 (2): 132-141.

[11] GHIDOSSI R, VEYRET D, MOULIN P. Computational fluid dynamics applied to membranes: state of the art and opportunities [J]. Chemical Engineering and Processing: Process Intensification, 2006, 45 (6): 437-454.

[12] HAPPEL J. Viscous flow relative to arrays of cylinders [J]. AIChE Journal, 1959, 5 (2): 174-177.

[13] GüNTHER J, SCHMITZ P, ALBASI C, et al. A numerical approach to study the impact of packing density on fluid flow distribution in hollow fiber module [J]. Journal of Membrane Science, 2010, 348 (1): 277-286.

[14] GüNTHER J, HOBBS D, ALBASI C, et al. Modeling the effect of packing density on filtration performances in hollow fiber microfiltration module: a spatial study of cake growth [J]. Journal of Membrane Science, 2012, 389: 126-136.

[15] TANSEL B, BAO W Y, TANSEL I N. Characterization of fouling kinetics in ultrafiltration systems by resistances in series model [J]. Desalination, 2000, 129 (1): 7-14.

[16] SERRA C, CLIFTON M J, MOULIN P, et al. Dead-endultrafiltration in hollow fiber modules: module design and process simulation [J]. Journal of Membrane Science, 1998, 145 (2): 159-172.

[17] MENDRET J, GUIGUI C, CABASSUD C, et al. Numerical investigations of the effect of non-uniform membrane permeability on deposit formation and filtration process [J]. Desalination, 2010, 263 (1/2/3): 122-132.

[18] ZHANG L Z. Flow maldistribution and performance deteriorations in membrane-based heat and mass exchangers [J]. Journal of Heat Transfer, 2009, 131 (11): 111801.

[19] ZHANG L Z, LI Z X, ZHONG T S, et al. Flow maldistribution and performance deteriorations in a cross flow hollow fiber membrane module for air humidification [J]. Journal of Membrane Science, 2013, 427: 1-9.

[20] LI Z X, ZHANG L Z. Flow maldistribution and performance deteriorations in a counter flow hollow fiber membrane module for air humidification/dehumidification [J]. International Journal of Heat and Mass Transfer, 2014, 74: 421-430.

[21] CHELLAM S, JACANGELO J G, BONACQUISTI T P. Modeling and experimental verification of pilot-scale hollow fiber, direct flow microfiltration with periodic backwashing [J]. Environmental Science & Technology, 1998, 32 (1): 75-81.

2016-01-19收到初稿，2016-03-09收到修改稿。

聯(lián)系人:戴干策。第一作者:莊黎偉（1988—），男，博士研究生。

Received date: 2016-01-19.

中圖分類號(hào):TP 273

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):0438—1157（2016）07—2841—10

DOI:10.11949/j.issn.0438-1157.20160088

Corresponding author:Prof. DAI Gance, gcdai@ecust.edu.cn

Numerical simulation of dynamic process during outside-in dead-end filtration in hollow fiber membrane module

ZHUANG Liwei, DAI Gance
(State Key Laboratory of Chemical Engineering, East China University of Science and Technology, Shanghai 200237, China)

Abstract:A CFD model was developed based on the filtration in the dead-end outside-in hollow fiber membrane module. Various fiber length, diameter, permeability, packing density, fouling index and transmembrane pressure were chosen during numerical simulation of the dynamic evolution of flux distribution and permeate volumetric flow rate. The simulation revealed that the uniformity of flux distribution improved as the filtration processes. The self-adjustment of the flux distribution was more pronounced with longer, narrower, more permeable fibers, higher packing density, fouling index and trans-membrane pressure. The inverse of the water volumetric flow rate increased linearly with the accumulated volume of the permeate. Due to the non-uniformity of the flux distribution and its dynamic evolution, the linear relationship differed from the one presented in the classic cake filtration model. A correlation equation to characterize the dead-end outside-in cake filtration in the hollow fiber membrane module was obtained through curve fitting of the simulation data. The equation enabled the prediction of the module performance and better design of the module.

Key words:hollow fiber membrane; CFD; model; numerical simulation; flux distribution; dynamic evolution