佟振鳴,張笑丹,張慧,袁竹林,劉成,張文良,金圣毅,陳超*
1 上海新型煙草制品研究院有限公司,上海 200080;2 東南大學能源與環境學院,江蘇南京 210096;3 上海煙草集團有限責任公司,上海 200082
電子煙[1]作為一種新型煙草產品在近些年來得到快速發展。其工作原理大致如下:電子煙霧化器通過電熱絲加熱具有毛細作用[2]的導油繩[3],導油繩的兩端與煙油[4]相接觸,煙油在毛細力的作用下輸送到電加熱區被加熱汽化而產生煙氣提供給消費者。電子煙的質量影響因素涉及兩個重要方面:一是煙油的成分組成,煙油是電子煙液的俗稱,通常由煙堿、水、丙二醇、丙三醇和其他添加劑組成[5],其品質對感官有重要影響;二是電熱霧化器[6]的工作性能,電熱功率、導油繩的毛細作用和傳熱特性都會對煙油的汽化速率產生影響,若電熱絲與導油繩接觸面溫度過高,將會導致與電熱絲周圍的煙油發生碳化,會嚴重影響煙氣的品質和吸食口感。因此了解與掌握導油繩內的溫度場、煙油濃度場以及煙油汽化量在不同加熱功率、加熱時間和導油繩擴散系數等參數條件下的變化規律對電子煙霧化器的設計與改進具有重要的意義。
由于電子煙霧化器的體積較小,采用實驗的方法精確測量導油繩內溫度和煙油濃度隨時間的變化規律相對困難,國內外部分學者嘗試通過各種實驗[7-10]方法測定電子煙霧化器的溫度特性,雖然獲得了一些初步結果,但由于受到測量手段的限制,迄今對導油繩內的傳熱傳質的研究仍不是很充分,尤其缺乏導油繩內溫度場、煙油濃度場和汽化量隨不同電加熱功率、加熱時間和導油繩擴散系數的定量變化關系。此外,由于抽煙是間歇性的行為,因此霧化器的工作過程是一非穩態的傳熱傳質過程,這更增加了導油繩內各物理場變化的復雜性?;谝陨犀F狀,本文采用數值模擬的方法,首先建立電熱絲與導油繩的傳熱模型、煙油在導油繩內的輸運模型以及煙油達到相變點后的汽化模型,而后對這些數學模型進行編程求解,獲得了相應電子煙霧化器工作條件下的導油繩內溫度場和煙油濃度場的非定常變化規律,煙油液體汽化量與加熱時間的關系,以及干燒時間點隨加熱功率和導油繩擴散系數的變化關系。本研究旨在為電子煙霧化器的設計和避免煙油過熱碳化提供依據。
導油繩是電子煙霧化器的核心部件,市場上多款主流電子煙(表1)均采用導油繩供油,導油繩的主要作用是基于內部微孔結構產生的毛細力將煙油傳輸至電熱絲以加熱霧化。對導油繩的數值模擬是基于霧化器的幾何結構及工作過程相關的物理定律和數學方程的一種算法模擬。圖1是本文所研究的電子煙霧化器的物理模型,主要包括煙油、導油繩等。

圖1 電子煙霧化器物理模型Fig. 1 Physical model of the electronic cigarette nebulizer

表1 通過導油繩供油的電子煙型號Tab. 1 Models of electronic cigarettes with guide rope
電子煙霧化器的工作過程包括以下三個環節:(1)電熱絲產生的熱量通過熱傳導方式傳遞給導油繩,并使其溫度升高;(2)當導油繩內的溫度達到煙油沸點時,煙油發生汽化并帶走汽化潛熱;(3)導油繩兩端的煙油在毛細力的作用下向電加熱區域進行質量傳遞,對發生汽化的煙油形成補充。
1.2.1 導油繩內熱量的傳遞
導油繩為圓柱形狀,由于溫度場、煙油濃度場沿圓周方向分布一致(圖2),因此建模時取導油繩軸線部位二維薄片進行建模。

圖2 導油繩二維簡化模型示意圖Fig. 2 Schematic diagram of two-dimensional simplified model of the guide rope
電子煙霧化器接通電源后,電熱絲通過導熱方式加熱浸有煙油的導油繩,煙油達到一定溫度后發生相變汽化。導油繩中的熱量傳遞屬于二維非穩態導熱過程,可用公式(1)進行描述[11]。

式中t為任意時刻液體的溫度,K;τ為時間,s;a=λ/ρc為熱擴散系數,m2/s,其中λ為導熱系數,W/(m·K),ρ為密度,kg/m3。
1.2.2 導油繩內煙油的汽化
導油繩內煙油液體的加熱過程可分為兩個階段:(1)當煙油液體溫度Tf低于其沸點Tfbp時,煙油吸收熱量后升溫,此階段內的煙油不發生蒸發相變;(2)當煙油液體溫度Tf達到其沸點Tfbp時,煙油液體的溫度保持不變,所吸收的熱量用于煙油液體的蒸發。蒸發前液態煙油的吸熱量為:

式中cf為液態煙油的比熱容,kJ/(kg·K);mf為液體的質量,kg;tf1、tf2分別為煙油液體在不同時刻的溫度,℃。
煙油蒸發的熱平衡方程為:

式中γf為煙油的蒸發潛熱,kJ/kg;m′f為煙油蒸發量,kg;Δτ為時間步長。
1.2.3 導油繩內的傳質
煙油液體在導油繩的毛細作用下從兩端傳入中間蒸發區,屬于二維非穩態傳質過程,可用如下微分方程表示[11]:

式中C為任意時刻液體的濃度,kg/m3;D為導油繩擴散系數(或稱傳質系數),m2/s,與導油繩毛細性能有關。電子煙霧化器中的液體主要在導油繩的毛細作用下傳輸,而導油繩屬于多孔介質,因此液體的擴散過程受導油繩孔隙率、孔徑大小等因素的影響,為了獲得液體在導油繩中的擴散系數,根據菲克定律采用如圖3所示實驗裝置(配合精密天平)進行估測。

圖3 導油繩擴散系數估測裝置圖Fig. 3 Schematic of the diffusion coefficient estimation device of the guide rope
為了直觀地看出煙油液體的擴散過程,在液體中加入色素。實驗前把導油繩從低到高分成相等的若干段,將導油繩下端浸入盛有煙油液體的燒杯中,在不同時間測量導油繩內的液體質量,根據導油繩內液體質量與時間的關系可擬合得到液體在導油繩中的擴散系數。
數值模擬采用Visual Basic(VB)編程實現,對物理模型進行了簡化并采用結構網格,因為結構比較簡單而省略了網格獨立性驗證。研究主要關注液體在導油繩內的非穩態傳輸過程,不包含相變后氣體的冷凝過程,因此在圖中沒有注明出口區。根據以上數學模型及相關求解方式,基于VB語言設計了模擬仿真軟件,軟件界面如圖4所示,分為參數輸入模塊與功能操作模塊。輸入參數包含主要包括導油繩尺寸、加熱功率、過熱度等,操作模塊可以執行既定參數的模擬仿真以及圖片保存等功能。

圖4 模擬仿真軟件Fig. 4 Simulation software
由于煙油液體在導油繩中的升溫、汽化和傳輸是非穩態過程,因此在計算溫度和濃度時要先離散傳熱傳質方程,利用SIMPLE[12]算法計算液體溫度和濃度。首先設置邊界條件和初始計算條件并計算溫度場,通過離散傳熱微分方程獲得不同網格的液體溫度;然后判斷每個網格的液體是否發生蒸發。隨后求解濃度場,液體蒸發使不同位置產生濃度梯度,煙油在導油繩產生的毛細力作用下進行傳遞,通過離散傳質微分方程計算獲得不同網格的液體濃度(圖5)。

圖5 導油繩的網格劃分Fig. 5 Meshing of the guide rope
為了能夠直觀地觀察導油繩內的溫度場和煙油濃度場,本研究利用計算機多媒體技術,把計算所得不同時刻各網格的溫度和濃度數據轉化為不同顏色及深度表示,以便清楚地觀察導油繩內溫度場和濃度場的變化規律。
煙油的主要成分是丙三醇,為了使模擬能夠貼合實際煙油的霧化過程,因此模擬計算選擇丙三醇代表煙油,設置時間步長為0.01 s,網格尺度為0.05 mm,液體初始溫度為20℃,初始濃度為1000 kg/m3。圖6是在表2所示計算條件下所獲得的導油繩內液體溫度場和濃度場隨時間的變化過程,從上至下為時間遞增方向。其中紅色表示導油繩內的溫度場,藍色表示導油繩內的煙油濃度場,顏色越深表明溫度、濃度值越高。

表2 模型的輸入參數Tab. 2 Input parameters of the model
由圖6可得:當電熱絲通電后,最接近電熱絲的導油繩溫度迅速提高(紅色不斷加深),升溫過程逐漸向軸線蔓延。當電熱絲附近溫度達到液體的沸點后,液體開始汽化,該區域導油繩溫度不再繼續升高。液體開始汽化時,導油繩內的液體濃度場開始變化,最靠近電熱絲區域的液體發生相變導致其濃度下降(藍色不斷變淺),該趨勢逐漸向導油繩的中軸線蔓延,說明液體汽化量逐漸增大。

圖6 不同時刻導油繩內的(a)溫度場及(b)濃度場Fig. 6 Field in the guide rope at different timing: (a) Temperature,(b) Concentration
為了驗證以上數值模擬結果的有效性和準確性,采用相同條件的實驗結果對現有結果加以驗證。謝國勇等[9]利用配備超細熱電偶的溫度檢測系統,建立了檢測電子煙霧化器核心區溫度的實驗系統,獲得了霧化器中煙油在霧化過程中的溫度變化。在所建立的霧化器的數值模擬平臺設置與實驗相同的輸入條件,將計算獲得的模擬結果與實驗結果對比從而驗證所建立的數學模型的準確性。
圖7是當加熱功率分別為6 W和10 W時實驗和模擬獲得的丙三醇溫度隨時間的變化圖,從圖中可以看出:實驗值和模擬值隨時間變化的總體趨勢基本一致,液體的溫度變化都分為兩個階段:升溫段和穩定段,與實驗結果對比可以發現,模擬中液體的溫度要比實驗略高,這可能是因為在模擬計算時沒有考慮熱量損失,但是兩項結果平均相對誤差不超過5%,因此可以證明所建立模型是準確和有效的。

圖7 實驗和模擬獲得丙三醇的溫度變化Fig. 7 Variation of the glycerol temperature obtained by experiment and simulation
2.2.1 導油繩內的非定常傳熱傳質特性
從圖6的數值模擬結果可見,隨著電加熱時間的延長,導油繩徑向的溫度梯度和液體的濃度梯度均發生相應的變化,圖8給出了導油繩徑向最大溫度梯度及最大濃度梯度隨時間的變化關系,從圖8(a)可以看出徑向最大溫度梯度先迅速增大,然后平穩一段時間后出現減小,這是由于熱量不斷從電熱絲的位置向導油繩中心位置擴散,處于邊緣的液體比中心液體升溫快,因此開始階段的溫度梯度越來越大,但是在邊緣液體溫度到達到沸點后發生汽化的過程中其溫度基本不變或出現輕度過熱,中心液體溫度升高速度較慢,所以梯度出現相對平穩段,而隨著中心液體溫度不斷升高且蒸汽釋放導致其過熱度不再增加,最大溫度梯度又逐漸下降。擬合得到徑向最大溫度梯度隨時間變化曲線可得其數學關系式為


圖8 導油繩徑向最大梯度隨時間的變化Fig. 8 Variation of maximum radial gradient of guide rope with time
從圖8(b)可以看出徑向最大濃度梯度在前3 s內基本保持不變,說明該階段內液體尚未發生蒸發相變,而隨后徑向最大濃度梯度隨時間快速增大,說明該階段內液體開始發生蒸發相變,隨后濃度梯度的平穩反映出導油繩供油量與液體蒸發量大致平衡,最后階段梯度的下降說明導油繩供液量略顯不足。最大濃度梯度隨時間變化與最大溫度梯度隨時間變化各階段匹配性較好,擬合得到徑向最大液體濃度梯度隨時間的變化關系式為

(5)(6)兩式揭示了導油繩內徑向最大溫度梯度和濃度梯度隨時間近似成二次函數的變化規律,可為電子煙霧化器導油繩的選取、電熱絲功率的確定提供依據。
2.2.2 導油繩內液體汽化量的變化規律
針對表2的輸入條件,通過數值模擬所獲得的液體汽化量隨時間的變化關系如圖9所示,該結果與圖8所示結果各階段匹配性較好,再次反映出模型的有效性。據圖9可見電熱絲通電后的3 s時間內,丙三醇主要是升溫過程而沒有發生汽化,3 s后出現汽化且汽化量m(mg)隨時間τ(s)的延長而增加,通過擬合得到汽化量與時間關系為:

圖9 丙三醇的汽化量隨時間的變化Fig. 9 Variation of vaporization amount of glycerol with time

加熱功率也是影響煙油汽化量的一個重要因素,圖10是在不同加熱功率下計算獲得的丙三醇汽化量隨時間的變化關系,從圖中可以看出在不同加熱功率下,丙三醇的汽化量與加熱時間近似成線性關系,而隨著加熱功率的增大,液體的汽化量逐漸增加。當加熱功率相差5 W時,汽化量相差兩倍左右,因此為了提高霧化效果,可以適當增大加熱功率。

圖10 在不同加熱功率下丙三醇的汽化量隨時間的變化Fig. 10 Variation of vaporization amount of glycerol with time under different heating powers
2.2.3 導油繩發生干燒現象的影響因素研究
電子煙在工作過程中,加熱功率過高和導油繩毛細作用差均會導致電熱絲與導油繩接觸面溫度過高而發生干燒現象,嚴重影響煙氣的品質和吸食口感。根據干燒產生的特點和條件將導油繩中心高溫點處液體濃度降為0時定義為干燒發生。從發生干燒的原因分析可知影響產生干燒時間點的主要因素為加熱功率和導油繩的擴散系數,通過數值模擬計算即可得不同條件下導油繩干燒時間點的變化曲線。
圖11(a)展示了在擴散系數為0.65e-9 m2/s時的導油繩干燒時間點隨加熱功率的變化,可見干燒時間點隨加熱功率增大而逐漸減小,因為隨著加熱功率增大,液體的升溫速率和濃度的變化速率都顯著增加,液體將很快達到沸點并發生汽化使得導油繩內液體濃度降低。圖11(b)則展示了在加熱功率為10 W時的導油繩干燒時間點隨擴散系數的變化,可知干燒時間點隨導油繩擴散系數的增大而逐漸增大,這是因為擴散系數反映了導油繩輸送煙油的能力,擴散系數越大則液體傳輸越快,對應的供液量越充足,所以干燒時間點會推遲。由此可知為了防止導油繩發生干燒,在其他條件不變的情況下,應盡量選擇傳輸性能好的導油繩材料。通過對相關曲線的擬合可得干燒時間點τ(s)隨加熱功率q(W)近似成指數變化,而隨擴散系數D(m2/s)近似成對數變化,具體關系式分別為:

圖11 干燒時間點隨加熱功率和擴散系數的變化Fig. 11 Variation of dry burning timing with heating power and diffusion coefficient

上式可為避免干燒提供設計基礎,也可用于指導霧化器加熱功率和導油繩的選取。
通過對電子煙霧化器導油繩內傳熱傳質過程的建模和求解分析,可將本文主要結論歸納如下:
(1)本文針對電子煙霧化器所建立的傳熱傳質數學模型能揭示導油繩內溫度場和煙油液體濃度場隨時間的變化規律,模型具有一定的正確性和可靠性。

(3)加熱功率和擴散系數是影響導油繩干燒的兩個重要因素。相應條件下的干燒時間點隨加熱功率的增大而近似成指數規律遞減,其關系式為τ= 68.89e-0.267q,干燒時間點隨擴散系數的增加而成對數關系遞減,其數學關系式為τ= 0.8152ln(D) + 2.7221。