CSTR相平面圖的繪制大作業(yè)教學(xué)文案

1、 化工過程分析與合成大作業(yè)CSTR 相平面圖的繪制指導(dǎo)老師：孫巍 張衛(wèi)東化工0703 楊福紳 200711084化工0706 丁志毅 200711159化工0706 沈楨琪 200711178化工0708 曹華瑜 200711246化工0708 馬 超 200711224目錄 TOC o 1-3 h z u HYPERLINK l _Toc266649844 摘要 PAGEREF _Toc266649844 h 1 HYPERLINK l _Toc266649845 關(guān)鍵詞 PAGEREF _Toc266649845 h 1 HYPERLINK l _Toc266649846 一、問題綜述 P

2、AGEREF _Toc266649846 h 1 HYPERLINK l _Toc266649847 二、反應(yīng)動力學(xué)模型及其反應(yīng)器設(shè)計 PAGEREF _Toc266649847 h 3 HYPERLINK l _Toc266649848 2.1反應(yīng)本征動力學(xué)模型 PAGEREF _Toc266649848 h 3 HYPERLINK l _Toc266649849 2.1.1模型假設(shè) PAGEREF _Toc266649849 h 3 HYPERLINK l _Toc266649850 2.1.2反應(yīng)動力學(xué)模型 PAGEREF _Toc266649850 h 3 HYPERLINK l _T

3、oc266649851 2.2反應(yīng)器設(shè)計 PAGEREF _Toc266649851 h 6 HYPERLINK l _Toc266649852 三、算法 PAGEREF _Toc266649852 h 9 HYPERLINK l _Toc266649853 四、相圖的繪制 PAGEREF _Toc266649853 h 10 HYPERLINK l _Toc266649854 五、雅克比矩陣判定模型局部穩(wěn)定性 PAGEREF _Toc266649854 h 11 HYPERLINK l _Toc266649855 六、結(jié)果分析 PAGEREF _Toc266649855 h 12 HYPER

4、LINK l _Toc266649856 6.1穩(wěn)態(tài)點的個數(shù)討論 PAGEREF _Toc266649856 h 12 HYPERLINK l _Toc266649857 6 .2進料狀態(tài)影對穩(wěn)定點與消耗時間的影響 PAGEREF _Toc266649857 h 13 HYPERLINK l _Toc266649859 七、參考文獻 PAGEREF _Toc266649859 h 14 HYPERLINK l _Toc266649860 八、附錄 PAGEREF _Toc266649860 h 15 HYPERLINK l _Toc266649861 附錄一 主要符號說明 PAGEREF _T

5、oc266649861 h 15 HYPERLINK l _Toc266649862 附錄二 程序 PAGEREF _Toc266649862 h 16摘要：本文對CSTR 中進行的放熱二級反應(yīng)進行了初步研究，建立了簡單模型，對CSTR 的操作進行了動態(tài)模擬，并采用龍格-庫塔法對此模型進行了求解，把不同時間下溫度和濃度的瞬時變化值繪出CSTR的相平面圖，并對模型進行了評價和推廣。關(guān)鍵詞：CSTR ，相平面圖 ，局部穩(wěn)態(tài)化一、問題綜述CSTR(Constant Str Tank Reactor )即全混流反應(yīng)器，又稱為全混釜或連續(xù)流動充分攪拌槽式反應(yīng)器，是一類在工業(yè)生產(chǎn)中廣泛應(yīng)用的反應(yīng)器。CST

6、R操作方式為連續(xù)進料、連續(xù)反應(yīng)、連續(xù)出料，為帶有攪拌槳葉的槽式反應(yīng)設(shè)備，攪拌的目的在于使物料體系達到均勻狀態(tài)，以有利于反應(yīng)的均勻和傳熱。在穩(wěn)態(tài)操作時，反應(yīng)器同一部位的操作參數(shù)不隨時間而變，因此可視為集中參數(shù)系統(tǒng)，有利于產(chǎn)品質(zhì)量控制和過程自動控制。與間歇反應(yīng)器操作方式不同，沒有裝料、卸料、升溫等不發(fā)生化學(xué)反應(yīng)的輔助時間，因而生產(chǎn)能力較大，輔助勞動少。原料加入后立即與釜反應(yīng)敏感的化學(xué)反應(yīng)，不會引起副反應(yīng)，由于釜內(nèi)物料容量大，當(dāng)進料條件發(fā)生一定程度波動時，釜內(nèi)反應(yīng)條件不會發(fā)生明顯變化，穩(wěn)定性好，操作安全。反應(yīng)物流入反應(yīng)器在瞬間與反應(yīng)器內(nèi)物料混合均勻，即在反映其中各處物料的溫度、濃度都相同。CSTR的

7、特性有：（1）物料在反應(yīng)器內(nèi)充分返混；（2）反應(yīng)器內(nèi)各處物料參數(shù)均一；（3）反應(yīng)器的出口組成與器內(nèi)物料組成想通；（4）連續(xù)、穩(wěn)定流動，是一定態(tài)操作。模型化（Modeling）是現(xiàn)代化工方法論的重要組成部分，尤其是過程動態(tài)學(xué)的核心。針對CSTR反應(yīng)器的特性，我們通過對所研究的反應(yīng)器、反應(yīng)系統(tǒng)及關(guān)鍵組分進行質(zhì)量、能量橫算，二級化學(xué)反應(yīng)速率表達式和化學(xué)反應(yīng)平衡常數(shù)計算式，從而建立CSTR反應(yīng)器內(nèi)的反應(yīng)的確定性集中參數(shù)模型，以此對在CSTR內(nèi)進行的反應(yīng)進行描述。為了通過確定性模型認識化工過程（乙酸乙酯合成）的內(nèi)在規(guī)律，解決與動態(tài)特性有關(guān)的工程實際問題，我們根據(jù)分析建立連續(xù)攪拌釜式反應(yīng)器的動態(tài)模型，進行

8、必須的數(shù)學(xué)處理，其中涉及到代數(shù)方程組、常微分方程組和偏微分方程組，二級放熱反應(yīng)混合方程組的求解。以對CSTR反應(yīng)器建立的確定性幾種模型的數(shù)學(xué)表達式所構(gòu)成的數(shù)學(xué)模型的正問題，求模型方程組的解析解是不可能的，不得不借助計算機求數(shù)值解。對于典型的常微分方程的初值問題，通?？梢岳谬埜?庫塔（R-K）、基爾（Gear）法等通用程序求數(shù)值解。反應(yīng)的確定性集中參數(shù)模型，應(yīng)用相關(guān)計算法則，編制計算機程序，并利用計算機輔助本篇闡述了我們通過文獻調(diào)研，數(shù)據(jù)查找及相關(guān)資料，對一個在CSTR反應(yīng)器內(nèi)進行的二級不可逆放熱反應(yīng)進行分析、建立起設(shè)計繪制了狀態(tài)空間的相平面圖，實施了對這一反應(yīng)的狀態(tài)空間分析和熱平衡分析，并對

9、所得結(jié)果進行討論和綜合的過程。對整個過程中，我們始終堅持用系統(tǒng)工程的思想、方法來解決化工過程系統(tǒng)的設(shè)計、開發(fā)、操作、控制等問題，取得了可惜的成績，使我們加深了對化工過程分析與合成的認識以及對化工過程系統(tǒng)動態(tài)模擬與分析的理解，為我們在今后的學(xué)習(xí)和工作時間奠定了思想基礎(chǔ)。模型建立后，在對過程模型的描述中，由于表現(xiàn)出強烈的非線性特性，意味著系統(tǒng)將會出現(xiàn)多重定態(tài)點。因此，對狀態(tài)空間的分析，運用了圖解的方法，既相圖。它是反映了我們所關(guān)心的狀態(tài)變量變化范圍內(nèi)，系統(tǒng)所有動態(tài)學(xué)定性特征的圖形。通過對相圖的分析，我們可以直觀的找到定態(tài)點、穩(wěn)定定態(tài)點、不穩(wěn)定定態(tài)點等信息，以便應(yīng)用于系統(tǒng)開工、調(diào)優(yōu)等方面。二、反應(yīng)動

10、力學(xué)模型及其反應(yīng)器設(shè)計2.1反應(yīng)本征動力學(xué)模型乙酸乙酯是一種重要的基礎(chǔ)有機化工原料，在造漆、人造革、醫(yī)藥和塑料等工業(yè)中有廣泛應(yīng)用。傳統(tǒng)的乙酸乙酯生產(chǎn)工藝是以乙醇和乙酸為原料，濃硫酸作催化劑。雖然濃硫酸對酯化反應(yīng)的催化活性較高，工藝成熟，但對生產(chǎn)設(shè)備腐蝕嚴重，易產(chǎn)生副反應(yīng)，而且產(chǎn)生大量酸性廢水。傳統(tǒng)酯化工藝中由于反應(yīng)受平衡的限制，往往采用乙醇過量的方法以提高乙酸的單程轉(zhuǎn)化率，使后續(xù)產(chǎn)品分離負荷大，能耗高。目前有關(guān)乙酸乙酯合成的新工藝研究不斷取得進展，但主要集中在用非濃硫酸催化劑?；趯嶒灁?shù)據(jù)，我們建立了本反應(yīng)的動力學(xué)模型。2.1.1模型假設(shè)a.催化劑量對反應(yīng)過程有影響，需要定義一個(x)的函數(shù)對

11、此予以描述，x為催化劑質(zhì)量與乙酸質(zhì)量比；b.由于催化劑粒度很小，故整個反應(yīng)按照擬均相反應(yīng)處理；c.假設(shè)反應(yīng)釜內(nèi)料液瞬間加滿，則Ca=Caf.2.1.2反應(yīng)動力學(xué)模型乙酸與乙醇酯化反應(yīng)如下：CH3COOH +CH3CH2OH CH3COOCH2CH3+ H2O本動力學(xué)模型考慮催化劑量對反應(yīng)過程的影響，按擬均相二級可逆反應(yīng)處理。反應(yīng)速率可描述為：-r=(x)k0(CaCb CeCwK)(mol (Lmin) (1)(x)為催化劑對反應(yīng)速率的影響因子，無因次；x為催化劑質(zhì)量與乙酸質(zhì)量比；為活化能，Jmol；R為8314 J(molK)；C表示濃度，molL；下標a，b，e和w 分別表示乙酸， 乙醇，

12、乙酸乙酯和水；K 表示化學(xué)平衡常數(shù)。圖1. 反應(yīng)溫度對反應(yīng)速率的影響 圖2.催化劑用量對反應(yīng)的影響根據(jù)實驗測定的平衡濃度，我們做出圖一與圖二，擬合數(shù)據(jù)并回歸得：(x)=.= 74534584.79L/（molmin）, =40159.035J/mol。為了求出平衡常數(shù)，在催化劑用量為500 g，乙醇和乙酸各為10 mol，反應(yīng)溫度為3372 K的條件下進行實驗，每隔2O min取樣分析乙酸的濃度，直至乙酸的濃度不再發(fā)生變化為止，此時各組分的濃度作為平衡濃度，實驗結(jié)果如圖3所示。圖3.轉(zhuǎn)化率與時間的關(guān)系對于二級均相反應(yīng)，轉(zhuǎn)化率與時間關(guān)系為：根據(jù)圖3，擬合得到反應(yīng)常數(shù)=3.94。綜上所述，反應(yīng)本征

13、動力學(xué)模型可以描述為：-r= 74534584.79 (CaCb CeCw394)mol/(Lmin)在我們設(shè)計的CSTR中，每一mol乙酸配128g催化劑，催化劑粒度為0.053-0.074mm，反應(yīng)在T=347K下進行，此時反應(yīng)平衡轉(zhuǎn)化率=0.667。假設(shè)乙酸：乙醇=1:1（摩爾比）進料，則-r=1819.88 ()對CSTR 進行物料衡算和熱量衡算，得到如下三個微分式：dCadt=FVCaf-Ca-rdCedt=-FVCe+rdTdt=FVTf-T-UAVCpT-Tc+(-HCp)r2.2反應(yīng)器設(shè)計根據(jù)進料參數(shù)的設(shè)定，通過查閱得到247K下反應(yīng)體系各物質(zhì)物理性質(zhì)如下： 表1.247K反應(yīng)

14、體系各物質(zhì)物理性質(zhì)密度/kgm3黏度/Pas定壓比熱容/Jkg-1K-1熱導(dǎo)率/Wm-1K-1乙醇748.6 0.0005222613.00 0.1653乙酸984.0 0.0005842311.67 0.15728乙酸乙酯831.5 0.0002562068.47 0.12545攪拌側(cè)水971.8 0.0003554195.00 0.674表2.303K空氣的物理性質(zhì)密度/kgm3黏度/Pas定壓比熱容/Jkg-1K-1熱導(dǎo)率/Wm-1K-1空氣1.1650.000018610050.02673根據(jù)各反應(yīng)物特性及反應(yīng)機理，我們選擇了容積為內(nèi)徑1000mm，900L的鋼制CSTR反應(yīng)器，攪拌器

15、轉(zhuǎn)速為250r/min；由于反應(yīng)放熱微弱，故反應(yīng)器夾套采用空氣冷卻。反應(yīng)器與環(huán)境的傳熱量可以由以下公式進行計算：Q=UAT其中，Q總傳熱量（包括壁面熱傳導(dǎo)和兩側(cè)對流傳熱），W U總傳熱系數(shù)，W/（m2K） A總傳熱面積，m2 T傳熱溫差，K。由于反應(yīng)器內(nèi)反應(yīng)過程的的熱量變化還包含有物料輸入輸出所帶熱量和反應(yīng)熱，因此反應(yīng)器內(nèi)溫差較為復(fù)雜，故熱量選擇夾套側(cè)空氣溫差進行求算，換熱面積為反應(yīng)器有效面積，所以傳熱量的求算就變成傳熱系數(shù)的求算。換熱總傳熱系數(shù)由以下公式計算，1U=11+i=1nii+12+Rs1+RS2其中， 1反應(yīng)器內(nèi)對流傳熱系數(shù)， J/（m2K） 2夾套內(nèi)對流傳熱系數(shù)， J/（m2K）

16、 i壁厚，m i器壁材料熱導(dǎo)率，W/（mK） Rs1反應(yīng)器側(cè)污垢熱阻 Rs2夾套側(cè)污垢熱阻。器壁熱阻與反應(yīng)器材料有關(guān)，器壁兩側(cè)對流傳熱系數(shù)則跟各側(cè)物系狀態(tài)有關(guān)，因此，只需計算出1、2，即可得到出傳熱系數(shù)，為此，我們做出幾點假設(shè)：忽略反應(yīng)副反應(yīng)，忽略催化劑的影響，忽略混合熱的影響；空氣在夾套內(nèi)充分湍流；物理性質(zhì)取平均溫度下各物質(zhì)的物理性質(zhì)的加和平均數(shù)，既M=xiMi；反應(yīng)器內(nèi)衡為常壓； 圖4.反應(yīng)器模型反應(yīng)器全新，即忽略污垢熱阻；反應(yīng)器內(nèi)物系混合均勻，各點傳熱相等；反應(yīng)器材料符合國家標準。表3.反應(yīng)器設(shè)備參數(shù)反應(yīng)器內(nèi)徑D/m反應(yīng)器釜高/m反應(yīng)器容積V/m3反應(yīng)器材料11.376900鋼制材料攪拌

17、器直徑d/m攪拌器位置C/m轉(zhuǎn)速N/rmin-1夾套直徑Dj/m0.333 0.333 2501.1攪拌器側(cè)傳熱系數(shù)由經(jīng)驗公式求得1=0.33DNd223Cp13w0.14dD-0.25CH0.15d攪拌器直徑，mD反應(yīng)器內(nèi)徑，m流體密度，kg/m3N轉(zhuǎn)速，r/sCp流體的定壓熱容比， J/（kgK）流體在主體平均溫度下的熱導(dǎo)率，W/（mK）液體在主體平均溫度下的黏度，Pasw壁溫下流體黏度，Pa/s 取w0.14=1.05C/H常數(shù)，0.33。當(dāng)溫度在347K時，轉(zhuǎn)化率在0.65的情況下，帶入各參數(shù)得攪拌器側(cè)傳熱系數(shù)1=1934.981 J/（m2K）夾套側(cè)傳熱系數(shù)由經(jīng)驗公式求得2=0.02

18、3ddu0.8Cp0.4液體的熱導(dǎo)率，W/（mK）u空氣流速，m/s，取15 m/sd管子內(nèi)經(jīng)，（D2j-D2）/D，m設(shè)空氣進口溫度298K，出口溫度308K，帶入換熱介質(zhì)各參數(shù)后得夾套側(cè)換熱系數(shù)為2=43.700 J/（m2K）故得到反應(yīng)器與環(huán)境換熱總傳熱系數(shù)為U=42.835 J/（m2K）三、算法對我們所研究的二級可逆放熱反應(yīng)，在選定的反應(yīng)釜，一定的換熱狀況下，以及在進料溫度Tf，進料濃度CAf，進料流量F一定的情況下，無論釜內(nèi)初始濃度與初始溫度怎樣改變，最后均達到一個穩(wěn)定點。次狀態(tài)點及為在選定狀況下經(jīng)開工過程之后，CSTR反應(yīng)器穩(wěn)定運行的狀態(tài)點。它表現(xiàn)了集中參數(shù)模型CSTR反應(yīng)器的穩(wěn)

19、定性，即該反應(yīng)器操作受到外來干擾后的平衡能力。原微分方程組為：dCadt=FVCaf-Ca-rdCedt=-FVCe+rdTdt=FVTf-T-UAVCpT-Tc+(-HCp)r-r=1819.88 ()邊值條件為t=0時，Ca=Ca0,Ce=0 , T=T0解此微分方程我們應(yīng)用四級四階標準（古典、典則）龍格-庫塔方法我們在xn,n+1區(qū)間內(nèi)選取四個點xn, xn+12,xn+12 , xn+1的斜率值K1, K2 ,K3 ,K4加權(quán)平均生成平均斜率。在此方程組環(huán)境下，K1, K2 ,K3 ,K4的具體形式是K1=yn-2xnynK2=yn+h2K1-2xn+hyn+h2K1K3=yn+h2K

20、2-2xn+hyn+h2K2K4=yn+hK3-2(xn+h)yn+hK3其中，h為時間的步長四、相圖的繪制相圖可以直觀地反應(yīng)出定態(tài)點的位置和局部穩(wěn)定性。我們的目的是找到在一定初始條件下系統(tǒng)所能達到的穩(wěn)定狀態(tài)，因此我們可將常微分方程組，對于給定的一組初值（Ca0,Ce0,T0），可以通過龍格庫塔（RK）法,用C+求出方程組的數(shù)值解：t0 t1 t2 ti Ca0 Ca1 Ca2 Cai Ce0 Ce1 Ce2 Cei T0 T1 T2 Ti 用matlab畫出相圖（），（Ca,T）能夠收斂于一定點（Cas, Ts）， Cas=36.2636mol/m3, Ces=73.7346mol/m3,T

21、=294.563K五、雅克比矩陣判定模型局部穩(wěn)定性對于原微分方程組：dCadt=FVCaf-Ca-rdCedt=-FVCe+rdTdt=FVTf-T-UAVCpT-Tc+(-HCp)r其中r=-1819.88 ()dCadt=f(Ca,Ce,T)dCedt=z(Ca,Ce,T)dTdt=g(Ca,Ce,T)fCa=-FV-2k(Ts)CasfCe=-2k(Ts)Ce3.94fT=k(Ts）()-40159.035R-1TS2zCa=-2k(Ts)CaszCe=-2k(Ts)Ce/3.94zT=- k(Ts）()-40159.035R-1TS2gCa=2-HjCPk(Ts)CagCe=-2(-H

22、j)3.94CPk(Ts)CegT=-FV-UAVCP+-H kTsCP40159.035RTS2()k(Ts)= 1819.88列雅克比行列式，代入數(shù)值：A=fCAfCefTzCazCezTgCAgCegT=-0.016119319-0.0044541310.002176467-0.0099848740.004454131-0.002176467-1.29043E-05-2.2341E-12-0.006240503求得三個特征值1，2，3分別為1=-0.0180901， 2=-0.00642613，3=-0.00624167它們均具有負實部，證明模型局部穩(wěn)定性良好。六、結(jié)果分析6.1穩(wěn)態(tài)點的

23、個數(shù)討論通過放熱反應(yīng)研究CSTR 的穩(wěn)態(tài)問題，對于熱量衡算式，可以認為等號右邊即為反應(yīng)移熱速率G(T )與體系放熱速率R(T )之差，即dTdG(T ) R(T )G(T )= V0CPT-T0+UA(T-Tc)R(T )=VRCA0-Hr1+rtt=VR/V0當(dāng)反應(yīng)放熱速率G(T)與散熱速率R(T)相等時,可以認為反應(yīng)器已達到穩(wěn)態(tài)，對應(yīng)于G(T)線與R(T)線的交點。移熱速率為一條直線，放熱速率為曲線，在放熱量較大的反應(yīng)體系中將出現(xiàn)三個穩(wěn)態(tài)點，但由于本反應(yīng)體系放熱量極少，因此只存在一個穩(wěn)態(tài)點。6 .2進料狀態(tài)影對穩(wěn)定點與消耗時間的影響對于研究的二級可逆放熱反應(yīng)，在選定的反應(yīng)釜。一定換熱狀況下

24、，流量一定下，改變其進料濃度及進料溫度，最后均達到穩(wěn)定點：Cas=36.2636mol/m3, T=294.563K體現(xiàn)了集中參數(shù)模型CSTR反應(yīng)器的穩(wěn)定性，即改反映其操作受外界干擾的自衡能力。由下表可知，隨著進料濃度的上升，反應(yīng)需要的時間減少，這取決于一定量內(nèi)的高濃度，促進了反應(yīng)的進行；隨著進料溫度的降低，反應(yīng)需要的時間減少但不明顯，這取決于此反應(yīng)體系為放熱反應(yīng)，低溫利于反應(yīng)速率的提高及反應(yīng)時間的減少，然而，由于次反應(yīng)體系為反應(yīng)熱較少的反應(yīng)體系，因此溫度的影響效果并不明顯。表4.不同進料溫度與濃度下的反應(yīng)時間CA0/mol/m3102030405060t/hT/k3055.65395.600

25、85.53975.46865.38285.27503005.65255.59925.53815.46695.38175.27502955.65115.59755.53645.46535.38005.27392905.64945.59585.53445.46335.37815.27222855.64815.59425.53285.46115.37585.27002805.64645.59225.53065.45895.37335.2672圖6. 不同進料溫度與濃度下的相圖（添加時間軸）七、參考文獻1、陳志平，張序文等編著 攪拌預(yù)混合設(shè)備設(shè)計選用手冊 化學(xué)工業(yè)出版社 2004.5（1）2、陳敏恒，

26、叢德滋等編 化工原理（上冊）第三版 化學(xué)工業(yè)出版社 2006.5（3）3、馬沛生編著 有機化合物實驗物性數(shù)據(jù)手冊 化學(xué)工業(yè)出版社4、麻德賢，李成岳，張衛(wèi)東主編 化工過程分析與合成 化學(xué)工業(yè)出版社 2002.6（1）5、鄧建中，劉之行編 計算方法 西安交通大學(xué)出版社 2004.9(5)6、郭鍇，唐小恒等編 化學(xué)反應(yīng)工程 化學(xué)工業(yè)出版社 2008.1(2)7、廖安平，張雷等 強酸性離子交換樹脂催化合成乙酸乙酯動力學(xué) 化學(xué)反應(yīng)工程與工藝 24卷第4期 2008.8八、附錄附錄一 主要符號說明主要符號說明符號意義計量單位密度kgm-3粘度Pas催化劑對反應(yīng)速率的影響因子x催化劑質(zhì)量與乙酸質(zhì)量比Ea活化

27、能JmolR通用氣體常數(shù)J(molK)C濃度molm3Ca0初始反應(yīng)濃度molm3K化學(xué)平衡常數(shù)r化學(xué)反應(yīng)速率（molm3）2t時間sV體積m3對流傳熱系數(shù)W/(m2K)熱導(dǎo)率W/(mK)N轉(zhuǎn)速rs-1Cp比熱J(kgK)F進料量m3/sU與環(huán)境換熱的總傳熱系數(shù)W/(m2K)T反應(yīng)溫度KTf加料溫度KTc冷卻劑平均溫度KA換熱面積m2-H反應(yīng)熱效應(yīng)D反應(yīng)器內(nèi)徑mH反應(yīng)器高度m附錄二 程序程序說明：本模擬共有兩個程序源文件，其中RungeKuttaLib.h定義了龍格庫塔類，該類提供了解N維微分方程組通用方法，在數(shù)據(jù)的傳入傳出上使用了STL容器，達到了性能和穩(wěn)定性上較好的折中；main.cpp函

28、數(shù)是程序的主文件，它調(diào)用龍格庫塔類，實現(xiàn)了對CSTR反應(yīng)器模擬。另外，為方便使用程序，我們在附件中提供了編譯好的二進制程序（注：對已編譯好的二進制文件以init.txt和config.txt中控制模擬參數(shù)。），其中，RungeKuttaSimulation.exe實現(xiàn)了CSTR二維相圖的模擬（生成結(jié)果文件CSTR.txt）；RungeKuttaGetTheTime.exe實現(xiàn)了開工時間的計算（生成結(jié)果文件為time.txt）。RungeKuttaLib.h源碼如下：#ifndef INC_RUNGEKUTTALIB#define INC_RUNGEKUTTALIB#include#includ

29、e#include#includeusing namespace std;class RungeKuttapublic:RungeKutta();vectorvector Calculate();vector getState();vector getFinalState();bool IsReady();void setInit(vector);void setConfig(vector);void setf(double (*g)(vector);void clear();void setFindSteady(bool);void setFullSim(bool);static bool

30、OutSolution(vectorvector);static bool OutState(vector);static vector In(char);private:vectorvector _solution;vector vadd(vector x ,double y);vector getArray();double getElement(int index);vector init;vector now;vector steady;int dim;bool isCalcSteady;bool isSteady;bool haveF;bool isfull;double time;

31、double start;double end;double step;double (*f)(vector);void retime();RungeKutta:RungeKutta()dim =0;start = 0;end = 0;step =0;time = 0;haveF = false;isCalcSteady = false;isSteady = false;isfull = true;void RungeKutta:clear()dim =0;start = 0;end = 0;step =0;time = 0;haveF = false;isCalcSteady = false

32、;isSteady = false;isfull = true;vector RungeKutta:vadd(vector x ,double y)int i;for(i=0;ix.size();i+)xi = xi + y;return x;vectorvector RungeKutta:Calculate()vectorvector solution;if(!IsReady() return solution;vector temp;temp = this-now;temp.push_back(time);retime();solution.push_back(temp);double t

33、;for(t=this-start;tend;t+=this-step)temp = getArray();solution.push_back(temp);if(isCalcSteady)if(abs(temp0-now0)now = temp;retime();_solution = solution;return solution;vector RungeKutta:getState()return steady;vector RungeKutta:getFinalState()return now;bool RungeKutta:IsReady()return !(dim=0|end

34、= 0|step=0| haveF=false);void RungeKutta:setInit(vector vinit)this-init = vinit;this-now = this-init;this-dim = vinit.size();void RungeKutta:setConfig(vector config)this-start = config0;this-end = config1;this-step = config2;time = config0;void RungeKutta:setf(double (*g)(vector)this-f = g;this-have

35、F=true;vector RungeKutta:getArray()int i;double temp;vector sol;for(i=0;idim;i+)temp = getElement(i);sol.push_back(temp);sol.push_back(time);return sol;double RungeKutta:getElement(int i)double Xi,Xin,K1,K2,K3,K4;Xi = nowi;K1 = fi(now)*this-step;K2 = fi(vadd(now,K1/2)*this-step;K3 = fi(vadd(now,K2/2

36、)*this-step;K4 = fi(vadd(now,K3)*this-step;Xin = Xi +(K1+2*K2+2*K3+K4)/6;return Xin;void RungeKutta:retime()time = time + step;bool RungeKutta:OutSolution(vectorvector data)int iSize,jSize;ofstream fout(CSTR.txt);iSize = data.size();jSize = data0.size();for(int i=0;iiSize;i+)for(int j=0;jjSize;j+)fo

37、utdataijt;foutendl;return 0;bool RungeKutta:OutState(vector state)ofstream fout(state.txt);for(int i=0 ;istate.size();i+)foutstateit;foutendl;return 0;vector RungeKutta:In(char s)double temp;vector data;ifstream fin(s);while(1)fintemp;if(fin)data.push_back(temp);elsebreak;return data;void RungeKutta:setFindSteady(bool flag)isCalcSteady = flag;void RungeKutta:setFullSim(bool flag)isfull = flag;#endifmain.cpp源碼如下：