紗筒殘余氨的擴散過程建模與數值模擬

張 東, 孟 婥

(東華大學 機械工程學院， 上海 201620 )

紗筒殘余氨的擴散過程建模與數值模擬

張 東, 孟 婥

(東華大學 機械工程學院， 上海 201620 )

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

紗筒； 殘余氨回收； 液氨整理； 傳質數學模型； 數值模擬

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

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

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

1 紗筒傳質模型

1.1模型的提出

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

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

1.2殘余氨揮發模型的建立

由建模假設條件，根據多相傳質基本理論并結合紗筒的物理結構形式，建立紗線中殘余氨揮發模型，如圖1所示。

圖1 紗筒傳質模型Fig.1 Model of mass transfer

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

1.3數學模型的建立

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

(1)

式中：ρ為氨的質量濃度，kg/L；r為散發方向的坐標，m；t為傳質經歷的時間，s。

由假設條件5可知，紗層擴散初始條件為

ρ|t=0=ρ0

(2)

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

(3)

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

1.4歸一化處理

由于1.3中建立的數學模型難以用有限差分法進行數值計算，并且不易確定算法的精度與穩定性，因此，為了能與平板類傳質方程相類比、分析問題的本質，必須要轉換成直角坐標的數學方程，令

(4)

(5)

(6)

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

(7)

(8)

(9)

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

(10)

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

C|T=0=1

(11)

內邊界條件為

(12)

外邊界條件為

(13)

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

2 傳質方程離散化

傳質方程(1)的解析解求解比較復雜[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，下標i代表第i空間層，上標n代表第n時層，在數值計算中Cin即是節點(i,n)的數值。為了使式(10)對時間和空間都具有二階精度，采用對網格點(n,n+0.5)進行中心差分計算，得到式(10)的相容性隱式差分方程，如式(14)所示。

(14)

(15)

(16)

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

(17)

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

(18)

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

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

(19)

式中：

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

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

即g≤1

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

3 傳質參數確定

3.1擴散系數與分配系數

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

3.2對流擴散系數

由對流傳質Sherwood(Ss)準數[10]知：

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

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

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

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

式中：D0、P0、T0和Da、Pair、Tair分別是氨氣在標準狀態下和實驗條件下的擴散系數、壓強和溫度；u是空氣流流速，m/s；v是運動黏度，m2/s。由于Re<500 000，從而對流擴散系數為

當風速u=0.1 m/s時，k=2.3×10-3m/s；當u=0.5 m/s時，k=5.2×10-3m/s；當u=1 m/s時，k=7.3×10-3m/s。

4 數值計算結果

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

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

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

圖3示出不同風速不同時間時紗筒內的歸一化殘余氨質量濃度，圖中縱坐標C是無量綱質量濃度。可看出，提高風速可加快殘余氨的揮發速率，但隨著時間推移，揮發速率不斷降低，這時單純依靠提高風速來增加揮發速率已不可行，必須要提高環境溫度等其他因素才能達到工程需求。故本文可為優化風速提供借鑒增加達到節能減排的效果。

圖3 不同風速和不同時間時紗筒內的歸一化殘余氨質量濃度Fig.3 Normalized concentration distributions under conditions of different air flow velocities and different time

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

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

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

5 結 論

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

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

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

FZXB

[1] 雋海艷, 孟婥, 孫以澤, 等. 改性處理對亞麻纖維形態結構和結晶度的影響[J].棉紡織技術, 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] 宋偉, 孔慶媛, 李洪枚. 木家具表面揮發性污染物散發傳質特性[J].化工學報, 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散發過程模擬與傳質參數測定新方法[J].化工學報, 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] 林庚，倪文，于昕. 納米孔真空絕熱板一維穩定導熱數學推導及其數值分析[J].中國礦業, 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] 淮秀蘭，姜任秋. 快速瞬態傳質過程中非菲克效應的實驗與理論研究[J]. 工程熱物理學報, 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

長江學者與創新團隊發展計劃滾動支持項目(IRT_16R12);東華大學博士創新基金項目(16D310303)

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