紗筒殘余氨的擴(kuò)散過程建模與數(shù)值模擬

張 東, 孟 婥

(東華大學(xué) 機(jī)械工程學(xué)院， 上海 201620 )

紗筒殘余氨的擴(kuò)散過程建模與數(shù)值模擬

張 東, 孟 婥

(東華大學(xué) 機(jī)械工程學(xué)院， 上海 201620 )

為準(zhǔn)確把握液氨整理后紗筒中殘余氨的傳質(zhì)特性，提高氨的回收率與回收速率，首先建立了經(jīng)綜合考慮擴(kuò)散傳質(zhì)和對流傳質(zhì)機(jī)制的傳質(zhì)數(shù)學(xué)模型，對柱坐標(biāo)系下的傳質(zhì)方程進(jìn)行歸一化處理，將其轉(zhuǎn)換成易于分析的一維平板類傳質(zhì)方程。然后利用Crank-Nicolson隱式差分法對歸一化的傳質(zhì)方程進(jìn)行差分近似求解，證明了算法的穩(wěn)定性與收斂性。最后將數(shù)值計算結(jié)果與實驗數(shù)據(jù)對比，驗證模型的正確性。結(jié)果表明，提高風(fēng)速以及減小紗筒外徑都可加快氨的揮發(fā)，小風(fēng)速下氨揮發(fā)過程是一個近似穩(wěn)態(tài)過程。

紗筒； 殘余氨回收； 液氨整理； 傳質(zhì)數(shù)學(xué)模型； 數(shù)值模擬

對氨和天然纖維材料之間的改性機(jī)制和工藝過程的大量深入研究證實了氨對天然纖維的改性作用具有顯著效果[1]。紗筒從改性釜中改性取出后仍然有殘余氨存在，因此，需對殘余氨的擴(kuò)散過程進(jìn)行研究。

棉麻纖維紗線通常纏繞在筒芯上形成紗筒，纖維或紗線作為固體骨架，纏繞時的物理間隙形成紗筒內(nèi)部孔隙，因此，紗筒是一種多孔材料，殘留的液氨和氣氨為流體填充在紗筒中。目前對多孔材料中流體傳質(zhì)過程特性的研究主要集中在一維平板類傳熱傳質(zhì)方面[2-4]，鮮有對紗筒這種柱狀多孔材料傳質(zhì)特性的研究，并且柱坐標(biāo)系下傳質(zhì)方程的解析解極其復(fù)雜，這給數(shù)學(xué)計算及工程應(yīng)用帶來極大困擾。

本文首先建立了基于傳質(zhì)學(xué)原理的紗筒傳質(zhì)數(shù)學(xué)模型，能夠模擬紗筒中殘余氨的散發(fā)行為，明確反映其物理過程。然后通過坐標(biāo)變換使圓柱坐標(biāo)系下的問題轉(zhuǎn)換為平板問題，為圓柱狀材料的傳熱傳質(zhì)提供計算方法。再利用隱式差分法近似求解傳質(zhì)方程，并證明算法的穩(wěn)定性與收斂性。最后將數(shù)值模擬結(jié)果與實驗數(shù)據(jù)作對比，對其傳質(zhì)特性進(jìn)行分析，指出影響殘余氨揮發(fā)速率的因素，以期為工程實踐中準(zhǔn)確預(yù)測紗筒殘余氨揮發(fā)過程提供模型與計算方法。

1 紗筒傳質(zhì)模型

1.1模型的提出

多孔材料中流體揮發(fā)具有相當(dāng)復(fù)雜的過程特性，現(xiàn)有研究多是根據(jù)流體力學(xué)基本理論和傳質(zhì)傳熱原理，對實際揮發(fā)過程進(jìn)行簡化，基于合理的假設(shè)條件建立數(shù)學(xué)模型，因此，對本文研究的改性后紗筒中殘余氨的揮發(fā)過程做以下建模條件假設(shè)。

1) 紗筒中紗線層視為均質(zhì)的等效圓柱體多孔材料，且只考慮該等效圓柱體內(nèi)部擴(kuò)散和表面對流擴(kuò)散[5]。2) 圓柱體內(nèi)的擴(kuò)散是由紗線層內(nèi)不同位置處液氨濃度差引起的，由于回收殘余氨過程中紗筒是首尾相接逐個堆砌的，故只考慮圓柱體內(nèi)徑向一維擴(kuò)散，忽略軸向擴(kuò)散。采用等效擴(kuò)散系數(shù)代替徑向擴(kuò)散系數(shù)。3) 假定紗層圓柱體表面存在一個界限層，該界限層靠近紗層一側(cè)，以液相狀態(tài)存在；靠近空氣流一側(cè)，殘余氨以氣相狀態(tài)存在。在界限層處液相氨和氣相氨發(fā)生相變，并處于平衡狀態(tài)，液相濃度和氣相濃度分配關(guān)系服從亨利定律[5]。4) 初始時刻，紗線層內(nèi)氨濃度是均勻的。5) 揮發(fā)過程中紗線層溫度和空氣流溫度均視為不變，整個過程中氨的質(zhì)量傳輸是連續(xù)的。

1.2殘余氨揮發(fā)模型的建立

由建模假設(shè)條件，根據(jù)多相傳質(zhì)基本理論并結(jié)合紗筒的物理結(jié)構(gòu)形式，建立紗線中殘余氨揮發(fā)模型，如圖1所示。

圖1 紗筒傳質(zhì)模型Fig.1 Model of mass transfer

紗筒簡化為空心圓柱體，以紗筒底面中心為原點，沿軸向和徑向建立柱坐標(biāo)系，紗層位于內(nèi)外半徑為a和b的同心環(huán)形體內(nèi)，紗筒外部為空氣流，其流動方向與紗層表面平行，且紗筒圓柱體表面各位置處的空氣流流速及方向相同，由于是數(shù)個紗筒堆砌放在改性釜中，可忽略上表面及下表面空氣流作用。紗層表面與外界空氣流之間存在邊界層，紗層內(nèi)液相氨以擴(kuò)散方式向邊界層傳輸，在邊界層內(nèi)液相氨轉(zhuǎn)化為氣相氨，并在邊界層外側(cè)通過對流擴(kuò)散方式擴(kuò)散到空氣流中。如圖1所示，殘余氨揮發(fā)模型將紗筒中氨揮發(fā)過程描述為由紗層到邊界層到空氣流層的3層結(jié)構(gòu)傳質(zhì)體系，其中邊界層內(nèi)通過動態(tài)吸附脫附作用在紗層表面形成液相氨與氣相氨的濃度分配，實現(xiàn)內(nèi)部擴(kuò)散與外部對流擴(kuò)散之間物質(zhì)傳輸過程的銜接。

1.3數(shù)學(xué)模型的建立

擴(kuò)散傳質(zhì)區(qū)內(nèi)，r=a處為紗層與筒芯交界面，在筒芯側(cè)液氨濃度梯度為零，紗層側(cè)液氨濃度為ρ(r,t)，簡記為ρ。其擴(kuò)散由濃度梯度與等效擴(kuò)散系數(shù)D決定，故由菲克第二定律知，紗層內(nèi)傳質(zhì)控制方程[6]為

(1)

式中：ρ為氨的質(zhì)量濃度，kg/L；r為散發(fā)方向的坐標(biāo)，m；t為傳質(zhì)經(jīng)歷的時間，s。

由假設(shè)條件5可知，紗層擴(kuò)散初始條件為

ρ|t=0=ρ0

(2)

邊界層指紗筒和空氣的接觸面，厚度可忽略，層內(nèi)氨濃度由紗層到空氣發(fā)生驟降，存在動態(tài)吸附平衡，控制方程邊界條件為

(3)

式中：a為紗筒外徑，m；b為紗筒的內(nèi)徑，m；ρ∞為空氣流中氨的質(zhì)量濃度，kg/L；D為等效擴(kuò)散系數(shù)，m2/s；k為空氣中氨的對流擴(kuò)散系數(shù)，m/s；β為氨在紗筒和空氣交界面處的濃度分配系數(shù)，為無量綱數(shù)。

1.4歸一化處理

由于1.3中建立的數(shù)學(xué)模型難以用有限差分法進(jìn)行數(shù)值計算，并且不易確定算法的精度與穩(wěn)定性，因此，為了能與平板類傳質(zhì)方程相類比、分析問題的本質(zhì)，必須要轉(zhuǎn)換成直角坐標(biāo)的數(shù)學(xué)方程，令

(4)

(5)

(6)

式中：T為無量綱時間；X為無量綱坐標(biāo)；C為無量綱質(zhì)量濃度。且T≥0，0≤X≤1，0≤C≤1，對式(4)～(6)進(jìn)行坐標(biāo)變換得

(7)

(8)

(9)

將式(7)～(9)代入式(1)，得到歸一化的紗筒內(nèi)擴(kuò)散控制方程

(10)

同理可得歸一化的初始條件為

C|T=0=1

(11)

內(nèi)邊界條件為

(12)

外邊界條件為

(13)

自此，已將模型中的物理方程進(jìn)行了歸一化處理，為后面驗證算法穩(wěn)定性與收斂性奠定了基礎(chǔ)。

2 傳質(zhì)方程離散化

傳質(zhì)方程(1)的解析解求解比較復(fù)雜[7]，在工程實踐中不易計算，本文采用有限差分法求解式(10)，使其更具實用性。式(10)是拋物線型偏微分方程，其影響域為半域，宜采用中心差分計算。將0≤X≤1等分成N份，從而△X=N-1，Xi=i△X，Tn=n△T，C(Xi,Tn)=C(i△X,n△T)簡記為Cin，下標(biāo)i代表第i空間層，上標(biāo)n代表第n時層，在數(shù)值計算中Cin即是節(jié)點(i,n)的數(shù)值。為了使式(10)對時間和空間都具有二階精度，采用對網(wǎng)格點(n,n+0.5)進(jìn)行中心差分計算，得到式(10)的相容性隱式差分方程，如式(14)所示。

(14)

(15)

(16)

即得到內(nèi)邊界條件。同理假設(shè)外虛擬點CN+1，由式(13)得：

(17)

將其代入式(15)中得到外邊界條件：

(18)

綜合式(15)、(16)、(18)得到時間推進(jìn)方式的代數(shù)方程組

C(n)=A-1BnC(0)

(19)

式中：

A、B為N+1階三對角方陣；C(n+1)、C(n)為N+1維向量。

由Von Neumann分析法即傅里葉分析法知，差分格式的誤差放大系數(shù)[8]為

即g≤1

式中，θ為相鄰2個節(jié)點上的諧波相位角之差，rad。從而可知式(14)的差分格式是穩(wěn)定的，由Lax定理可知差分格式也是收斂的[8]。三對角陣A滿足優(yōu)對角條件，代數(shù)方程組(19)的解亦恒穩(wěn)定，可取較大時間步長，計算省時，精度也滿足要求。通過解方程組(19)便可得出氨的濃度隨位置與時間的變化規(guī)律，再通過高精度數(shù)值積分，便可獲得任意時刻紗筒內(nèi)剩余氨的質(zhì)量，可與實驗數(shù)據(jù)作對比，驗證模型和算法的合理性。

3 傳質(zhì)參數(shù)確定

3.1擴(kuò)散系數(shù)與分配系數(shù)

擴(kuò)散系數(shù)D與紗筒內(nèi)紗線的組織結(jié)構(gòu)、纏繞方式以及松緊程度等因素相關(guān)。分配系數(shù)β描述了在紗筒邊界層液相氨與氣相氨的比例關(guān)系，其中D和β的取值可參考文獻(xiàn)[9]，D=1.777 8×10-6m2/s，β=8 340。

3.2對流擴(kuò)散系數(shù)

由對流傳質(zhì)Sherwood(Ss)準(zhǔn)數(shù)[10]知：

層流時Ss=0.664Sc1/3Re1/2(Re<500 000)

湍流時Ss=0.037Sc1/3Re4/5(Re>500 000)

式中：Sc=kzDa-1；雷諾數(shù)Re=uzv-1；z為紗筒高度，m；k為對流擴(kuò)散系數(shù)；Da為氨在空氣中的傳質(zhì)系數(shù)。

Da=D0P0Pair-1Tair3/2T0-3/2

式中：D0、P0、T0和Da、Pair、Tair分別是氨氣在標(biāo)準(zhǔn)狀態(tài)下和實驗條件下的擴(kuò)散系數(shù)、壓強(qiáng)和溫度；u是空氣流流速，m/s；v是運(yùn)動黏度，m2/s。由于Re<500 000，從而對流擴(kuò)散系數(shù)為

當(dāng)風(fēng)速u=0.1 m/s時，k=2.3×10-3m/s；當(dāng)u=0.5 m/s時，k=5.2×10-3m/s；當(dāng)u=1 m/s時，k=7.3×10-3m/s。

4 數(shù)值計算結(jié)果

改性前稱量紗筒，改性后再次稱量紗筒，二者相差即為紗筒內(nèi)殘余氨的初始質(zhì)量M0。將改性后的紗筒首尾相接堆砌放置在烘干爐中，在恒定溫度(293 K，使傳質(zhì)系數(shù)都是常量)下，再調(diào)控風(fēng)機(jī)風(fēng)速，不同風(fēng)速(u分別為0.1、0.5、1.0 m/s)條件下進(jìn)行殘余氨的回收，揮發(fā)出的氨氣通入水中，再經(jīng)過后處理回收氨。由于在持續(xù)通風(fēng)條件下ρ∞≈0，每隔1 h稱量1次，共稱量20次，結(jié)果取平均值，獲得了紗筒殘余氨揮發(fā)速率曲線，將所得實驗數(shù)據(jù)與數(shù)值計算結(jié)果作對比，結(jié)果如圖2所示。

圖2 實驗結(jié)果與模擬結(jié)果對比Fig.2 Comparison between experimental data and calculated data

圖2中M0是初始時刻紗筒中殘余氨的質(zhì)量，Mt是任意時刻紗筒中殘余氨的剩余質(zhì)量，縱坐標(biāo)Mt/M0即是紗筒內(nèi)殘余氨與總傳質(zhì)質(zhì)量之比。在開始一段時間內(nèi)，仿真數(shù)據(jù)與實驗數(shù)據(jù)有些偏差，主要原因是沒有考慮傳質(zhì)速度的有限性[11]，以及氨揮發(fā)時大量吸熱，導(dǎo)致溫度降低，揮發(fā)速率降低。在傳質(zhì)穩(wěn)定后，傳質(zhì)速度有限性對擴(kuò)散速率的影響可以不計，從而能很好地表達(dá)實際情況，故本文模型可作為工程實踐的計算模型。從圖還可看出，在開始的一段時間氨的剩余量下降得最快，表明揮發(fā)速率最快。

圖3示出不同風(fēng)速不同時間時紗筒內(nèi)的歸一化殘余氨質(zhì)量濃度，圖中縱坐標(biāo)C是無量綱質(zhì)量濃度?？煽闯?，提高風(fēng)速可加快殘余氨的揮發(fā)速率，但隨著時間推移，揮發(fā)速率不斷降低，這時單純依靠提高風(fēng)速來增加揮發(fā)速率已不可行，必須要提高環(huán)境溫度等其他因素才能達(dá)到工程需求。故本文可為優(yōu)化風(fēng)速提供借鑒增加達(dá)到節(jié)能減排的效果。

圖3 不同風(fēng)速和不同時間時紗筒內(nèi)的歸一化殘余氨質(zhì)量濃度Fig.3 Normalized concentration distributions under conditions of different air flow velocities and different time

圖4示出不同時間不同位置處的濃度曲線。由圖可知，在同一時刻紗筒內(nèi)外的濃度差別是很小的，其原因是殘余氨揮發(fā)過程是緩慢的，是一近似穩(wěn)態(tài)過程，但當(dāng)邊界層揮發(fā)速度較快即風(fēng)速極快時便不再是準(zhǔn)穩(wěn)態(tài)過程，必須要考慮到傳質(zhì)松弛時間的問題[11]。

圖4 不同時間不同位置處的濃度曲線Fig.4 Concentration curves under conditions of different time and positions

圖5 不同紗筒外徑下的傳質(zhì)速率曲線Fig.5 Mass transfer rate curves of different diameters

5 結(jié) 論

1)在綜合考慮了二大傳質(zhì)機(jī)制條件下，建立了紗筒殘氨揮發(fā)模型。通過對方程進(jìn)行歸一化處理使得柱狀傳質(zhì)問題轉(zhuǎn)變成了一維平板類傳質(zhì)問題，為柱坐標(biāo)系下的傳質(zhì)問題提供了一種新的求解方法。

2)利用Crank-Nicolson隱式差分法對歸一化的傳質(zhì)方程進(jìn)行差分近似求解，證明了算法的穩(wěn)定性與收斂性，并求得了與實驗數(shù)據(jù)吻合很好的結(jié)果，驗證了模型的正確性，可為工程實踐提供借鑒。

3)通過對計算結(jié)果的分析可知，提高風(fēng)速可加快氨的揮發(fā)，并指出了紗筒外徑對氨揮發(fā)時間和風(fēng)機(jī)空氣流利用率的影響，為工程實踐中優(yōu)化紗筒外徑提供了模型與計算基礎(chǔ)。

FZXB

[1] 雋海艷, 孟婥, 孫以澤, 等. 改性處理對亞麻纖維形態(tài)結(jié)構(gòu)和結(jié)晶度的影響[J].棉紡織技術(shù), 2013, 41(7):499-501. JUAN Haiyan, MENG Zhuo, SUN Yize, et al. Effect modifiction on morphology and crystallinity of flax fiber[J]. Cotton Textile Technology, 2013, 41(7): 499-501.

[2] 宋偉, 孔慶媛, 李洪枚. 木家具表面揮發(fā)性污染物散發(fā)傳質(zhì)特性[J].化工學(xué)報, 2013, 64(5): 1550-1558.SONG Wei, KONG Qingyuan, LI Hongmei. Characteristics of VOC mass transfer from wood furniture surface[J]. CIESC Journal, 2013, 64(5): 1550-1558.

[3] 宋偉, 孔慶媛, 李洪枚. 建材VOC散發(fā)過程模擬與傳質(zhì)參數(shù)測定新方法[J].化工學(xué)報, 2013, 64 (3):913-922. SONG Wei, KONG Qingyuan, LI Hongmei. New method for simulation of VOC emission from building materials and measurement of mass transfer para-meters[J]. CIESC Journal, 2013, 64(3): 913-922.

[4] 林庚，倪文，于昕. 納米孔真空絕熱板一維穩(wěn)定導(dǎo)熱數(shù)學(xué)推導(dǎo)及其數(shù)值分析[J].中國礦業(yè), 2008, 17(3):85-88. LIN Geng, NI Wen, YU Xin.Nonporous vacuum insulation panel one dimensional stable heat conduction mathematic deduction and numerical analysis[J]. Chinese Mining Industry, 2008, 17(3):85-88.

[5] RANA P, BHARGAVA R, BEG O A. Numerical solution for mixed convection boundary layer flow of a nanofluid along an inclined plate embedded in a porous medium[J]. Comput Math Appl, 2012,64(9):2816-2832.

[6] THORSTENSON D C, POLLOCK D W. Gas transport in unsaturated zones: multi-component systems and the adequacy of Fick′s laws[J]. Water Resources Research, 1989, 25(3):477-507.

[7] CRANK J. The Mathematics of Diffusion[M]. London: Clarendon Press, 1975:138-159.

[8] ANDERSON J D. Computational Fluid Dynamics: the Basics with Applications[M]. New York: McGraw-Hill, 1995:10-50.

[9] ZHU Helei, MENG Zhuo, SUN Yize, et al. Unsteady multi-phase mass transfer modeling and analysis for the process of the ammonia volatilization from yarn bobbin[J]. International Journal of Heat and Mass Transfer, 2016(103): 231-241.

[10] WHITE FM. Heat and Mass Transfer[M]. [S.l.]:Addison-Wesley, 1988:20-30.

[11] 淮秀蘭，姜任秋. 快速瞬態(tài)傳質(zhì)過程中非菲克效應(yīng)的實驗與理論研究[J]. 工程熱物理學(xué)報, 2000, 21(5):595-599. HUAI Xiulan, JIANG Renqiu. Theoretical and experimental study on the effect of the non-Fick in fast transient mass transfer process[J]. Journal of Engineering Thermophysics, 2000, 21(5):595-599.

Modelingandnumericalsimulationforresidualammoniavolatilizationfromyarnbobbin

ZHANG Dong， MENG Zhuo

(CollegeofMechanicalEngineering,DonghuaUniversity,Shanghai201620,China)

In order to accurately grasp the mass transfer characteristics of residual ammonia in yarn bobbins after liquid ammonia treatment and increase the recovery ratio and recovery rate of ammonia, a mathematical model of mass transfer was established by considering diffusion and convection. It became easier to analyze the mechanism of the emission of ammonia when the mass transfer equations was normalized. Then the Crank-Nicolson implicit difference method was used to solve the normalized mass transfer equations, and the stability and convergence of the algorithm were proved. Finally, the numerical results were compared with the experimental data, which proved the correctness of the model. The conclusion shows that the method of increasing wind velocity and decreasing the bobbin outer diameter can accelerate the volatilization of ammonia, and the ammonia volatilization process under low wind speed condition is an approximate steady state process.

yarn bobbin； residual ammonia recovery； liquid ammonia treatment; mathematical model of mass transfer； numerical simulation

TS 19

：A

10.13475/j.fzxb.20160903706

2016-09-21

：2017-06-16

長江學(xué)者與創(chuàng)新團(tuán)隊發(fā)展計劃滾動支持項目(IRT_16R12);東華大學(xué)博士創(chuàng)新基金項目(16D310303)

張東(1989—)，男，博士生。主要研究方向為流體及熱質(zhì)傳遞學(xué)。孟婥，通信作者，E-mail:mz@dhu.edu.cn。