1、第一章 COMSOL MULTIPHYSICS及數(shù)值分析基礎(chǔ)W. B. J. ZIMMERMAN , B. N. HEWAKANDAMBYDepartment of Chemical and Process Engineering, University of Sheffield,Newcastle Street, Sheffield S1 3JD United KingdomE-mail: w.本章主要介紹COMSOL Multiphysics在零維和一維模型數(shù)值分析方面的幾 個關(guān)鍵內(nèi)容。這些內(nèi)容包括求根、步進式數(shù)值積分、常微分方程數(shù)值積分和線性 系統(tǒng)分析。這幾乎是所有的化工過程數(shù)學(xué)分析方法

2、。下面通過COMSOLMultiphysics中的一些常見化工過程使用實例來介紹這些方法，包括：閃蒸、管 式反應(yīng)器設(shè)計、擴散反應(yīng)系統(tǒng)和固體中熱傳導(dǎo)。1 .簡介本章內(nèi)容很多，可以分為幾個不同的目標。首先介紹了 COMSOL Multiphysics的主要工作特性；其次介紹了如何使用這些特性來分析一些簡單的， 位于零維空間、一維空間或“空間一時間”系統(tǒng)中的化工問題。本章還希望通過 展示COMSOL Multiphysics和MATLAB工具在化工過程分析中的強大功能，激 發(fā)讀者對使用COMSOL Multiphysics進行建模和仿真的興趣。由于COMSOL Multiphysics不是一個通用的

3、問題求解工具，所以一些目標需 要迂回實現(xiàn)。作者在使用FORTRAN、Mathematica和MATLAB解決化工問題方 面有著豐富的教學(xué)經(jīng)驗，并用這些工具實現(xiàn)過這里所有的例子。而且，擴展化工問題的數(shù)值分析也已經(jīng)在POLYMATH1中實現(xiàn)，這似乎只在化工委員會的CACHE項目中使用過。本書前一版已經(jīng)介紹過在零維空間中求解非線性代數(shù)方程和和時間有關(guān)的 常微分方程的內(nèi)容。從概念上講，零維域就是一個簡單的有限元。 通過研究某一 特定有限元中的變化對理解有限元方法非常有用。但是，COMSOL Multiphysics通過獨立對話框設(shè)置，使得零維幾何方程和和時間相關(guān)的常微分方程求解變得非 常簡單。所以本章

4、將同時采用這兩種方法求解這些例子。2 .方法1:求根典型的數(shù)值分析課程會講解多種求根方法，但是從實際經(jīng)驗來看，只有兩種算法非常有用一一二分法和牛頓法。 我們這里沒有列出所有方法，而是重點考慮 為什么求根是最有效的數(shù)值分析工具。 在線性系統(tǒng)中求根非常簡單，但是對于非 線性系統(tǒng)這就是一個挑戰(zhàn)，而所有感興趣的動力學(xué)問題幾乎都是非線性系統(tǒng)。對 非線性系統(tǒng)的求根起源于對反函數(shù)的描述。 為什么呢？因為對于大多數(shù)非線性函 數(shù)，“正向" y=f(u)很好表示，但是它的反函數(shù)u=-1(y)可能不能顯式表示、多值 (無意義)或根本不存在。如果反函數(shù)存在的話，求解反函數(shù)其實就是求根的過 程一一求解滿足F(

5、u)=0的u等價于求解F(u)=f(u)-y=0。因為大多數(shù)數(shù)值分析的 目標是在系統(tǒng)約束下計算求解，所以這也等價于對所有的約束取反。COMSOLMultiphysics擁有求解非線性問題的核心函數(shù)femnlin ,本節(jié)主要介紹用它求解零維非線性問題。femnlin函數(shù)使用牛頓方法求解，由于只有一個變量u,牛頓法通過對一階倒數(shù)F'(u)迭代來求根。該方法首先估計函數(shù)的斜率范圍，然后再逼近根。該斜率可以通過理論分析（牛頓-拉夫遜方法）和數(shù)值（正割法）方法求得。如果能 用任何一種方法求得斜率，就可以用泰勒定律來逼近根。其基本思想就是使用目 前猜測值U0的泰勒展開式：(2)f (u) = f

6、(U0)(u -U0)f '(U0)IH該公式可以化簡，忽略（U-U0）的高階項，計算根如下:f（U0）f（U0）這個方法可以快速地擴展到多維求解空間，例如將 U看作未知矢量，“被 f（U0）除”看作“乘以f的雅克比矩陣的逆"。下一節(jié)介紹COMSOL MUltiphysics 中的求根過程。2.1 求根：COMSOL Multiphysics非線性求解器的使用實例如上節(jié)所述，求根本身是一個“零維”活動，至少對于“空間 -時間”系統(tǒng) 多維未知矢量U來說是這樣的。COMSOL多物理場沒有零維模式，所以我們臨 時采用一維模式。這在方面增加了我們不需要的冗余功能。 但是由于問題規(guī)模較

7、 小，COMSOL MUltiphysics編碼效率高，且現(xiàn)代微處理器的運算速度快，這點就 不成為問題了。啟動MATLAB并在命令窗口鍵入 COMSOL MUltiphysics。屏幕閃爍幾秒后， 會出現(xiàn)一個模型導(dǎo)航窗口。按照表1所示步驟，建立一個零維使用模式來解決非 線性多項式方程：u3 u2 -4u 2 = 0通過在“ Physics'菜單、"Subdomain settings”選項中的設(shè)定，使得表1中的每 個子域都滿足該方程。注意左上角，這是以矢量符號給出的方程。在一維模式下， 這個方程可以簡化為：cuc cuv u CU,、da - ,c+suf ：+au+P =f

8、(4)din 1nrct exexJtx顯然，a、丫和B在一維簡化里是多余的。既然我們是求零維的根，所有公式4左側(cè)的系數(shù)都可以設(shè)為00通過重新組合多項式，我們發(fā)現(xiàn) 2=4和£=口3+ U2+2。 注意我們是通過設(shè)定最大基元尺寸為1、將域離散化為只有一個基元的域，從而得到零維空間的！表1系數(shù)模式下的求根實例。文件名稱：roofinder.mph.Model Navigator選擇 1D 空間維數(shù)，COMSOL Multiphysics: PDE modes: PDE, coefficient 形式Draw菜單Specify objects： Line。Coordinates彈出菜單。x

9、： 0 1 名稱：interval,完成Physics 菜單：Boundary settings選擇域：1和2 （按住Ctrl鍵同時選中） 選擇Neumann邊界條件采用默認設(shè)置q=0 g=0,完成Physics 菜單：Subdomain settings選擇域：1設(shè)定：c=0; a=4;f=uA3+uA2+2; da =0使用，選擇Init選項卡：設(shè)定u(t0)=-2,完成Mesh菜單：Mesh parameters設(shè)定最大基元尺寸1 點擊remesh,完成Solve菜單：Solver parameters求解器選擇 Stationary nonlinear,求解，完成Post-proces

10、sing Point Evalution選擇邊界1,表達式：u,完成設(shè)定初始猜測值為u(t0)=-2,下面我們來尋找最接近這個值的根。你可能對為什么設(shè)置a=4而不是把所有有關(guān)參數(shù)放進f中感到奇怪，那是因為對公式(4)右側(cè)有 限元離散后得不到一個奇異剛度矩陣。后處理階段在輸出窗口中顯示出結(jié)果：Value: -2.732051, Expression: u, Boundary: 1.分析解得出的根為-1-V3,數(shù)值解能夠很好的和其吻合。根據(jù)代數(shù)中二次公式的解知道，另外一個根是-1 + T3,第三個根是1?；氐阶佑蛟O(shè)置中來，如果最初設(shè)置的猜測值為 u(t0)=-0.5,則 COMSOL Multip

11、hysics 解出 u=0.762051,又一 個很好的近似數(shù)值解。如果u(t0)=1.2作為最初猜想時，得到u=10該練習(xí)體現(xiàn)了非線性求解器和問題的兩個特性一一(1)非線性問題可以有多 個解；(2)最初猜測值對于求解很關(guān)鍵。牛頓法通常收斂到最接近的解，但是在 高階非線性問題中，這點可能不成立。這些特性在多維求解空間和“空間-時間”相關(guān)系統(tǒng)中同樣適用。COMSOL Multiphysics 模型 mph 文件(rootfinder.mph)中包含了所有的 MATLAB/COMSOL Script源代碼和FEMLAB的擴展功能，以便再次恢復(fù)到 FEMLAB當(dāng)前圖形用戶界面。該文件在網(wǎng)頁 http

12、:eyrie.shef.ac.uk/femlab中可以 得到。只需打開“Menu”菜單，選擇“Open model m-file”選項，在彈出的“Open file”對話框中選中它，你就可以在“ Subdomain settings”中快速設(shè)定非線性函 數(shù)，給定初始猜測值，然后用非線性求解器得到收斂解。但是如果函數(shù)中沒有線 性項可以放在公式(4)左側(cè)怎么辦？例如，F(xiàn)EMLAB將tan(u)-u2+5=0放在公式(4) 左側(cè)時，會生成奇異剛度矩陣。建議在子域中設(shè)定一個u的二次派生系數(shù)，c=1o 通過和Neumann邊界條件相結(jié)合，這個人為擴散因素不能改變基元中解為常數(shù) 的事實，但是它可以避免剛度

13、矩陣奇異。 在一般模式下求根tan(u)-u2+5=0中奇異剛度矩陣的問題可以通過一般模式來避免，它可以求解下 式：-2二 u .ea . 2 d a二 t-:t二F-X其中Ru,ux)和(4)式中的系數(shù)項功能類似，但是被求解器處理方式不同。在系數(shù) 模式中，認為系數(shù)和u無關(guān)，除非使用數(shù)值雅可比矩陣，這會產(chǎn)生一些非線性依賴關(guān)系一一迭代次數(shù)增加。一般模式下構(gòu)建剛度矩陣時，精確雅可比矩陣會將r和f對u進行微分。通常一般模式比使用數(shù)值雅可比矩陣的系數(shù)模式需要的迭 代次數(shù)少。下面使用的雅可比矩陣，不需要任何特殊處理來避免剛度矩陣奇異?？偟膩碚f，一般模式在處理非線性問題時比系數(shù)模式更加實用。從我的個人觀點

14、 看來，系數(shù)模式是COMSOL Multiphysics-一個遺傳特性 MATLAB的偏微 分方程工具箱，它在很多方面是 FEMLAB和COMSOL Multiphysics的先驅(qū)，它 廣泛使用了系數(shù)形式。止匕外，數(shù)值雅克比矩陣系數(shù)公式是一個長期存在的標準有 限元方法，所以相對于其它代碼，它是一種非常有用的公式。表2給出了使用一般模式的步驟一一對我們剛才的步驟做了小小改動。表2 一般模式下的求根實例。文件名：rfingen.mphModel Navigator1D, COMSOL Multiphysics: PDE modes, general 模式Options設(shè)定 Axes/Grid 為0

15、,1Draw名稱：Interval ； Start= 0; Stop=1Physics 菜單：Boundary Settings設(shè)定兩個端點(域)為Neumann邊界條件Physics 菜單：Subdomain Settings設(shè)定 r=0; da=0; F=u3+u2-4*u+2初值 u(t0)=-0.5Mesh模式設(shè)定最大基元尺寸，general： 1;點擊 RemeshSolve使用默認設(shè)置(nonlinear solver, exact Jacobiar)Post-Processing經(jīng)過5次迭代，結(jié)果收斂。點擊圖片讀出 u = 0.73205%改 變初始條件找到另外兩個根雖然這里給出

16、了單變量簡單函數(shù)求根的模板(rfindgen.mph),但是實際上MATLAB有更簡單的計算子程序，通過內(nèi)置函數(shù)fzero和函數(shù)聲明來求根，現(xiàn)在COMSOL Multiphysics圖形用戶界面可以使用輔助方程中的輔助變量來求解幾何約束，該變量稱為狀態(tài)變量。作為對一般模式求根模型的補充， 使用狀態(tài)變量 v按照表3的步驟來求解一個非線性方程。表3在OD殷置中求根Physics 菜單：ODE settings名稱：v。方程：tanh(v)-vA2+5,完成Solve菜單：Solver parameters選擇 Stationary nonlinear求解器，求解，完成Post-processing

17、選擇Point evaluation,邊界2,表遼式：vReport 窗口值： 2.008819,表達式：v,邊界：2下一節(jié)我們將非線性求根方法使用到一個常見的化工實例中一一閃蒸，它用到了 COMSOL Multiphysics圖形用戶界面的更多新特性。2.2 求根：閃蒸實例化學(xué)熱力學(xué)中廣泛存在求根法的使用實例，由于化學(xué)平衡和質(zhì)量守衡的約束 條件很充分，再加上狀態(tài)方程，就產(chǎn)生了和問題中未知變量個數(shù)相同的約束條件。 在本節(jié)中，我們以閃蒸為例來介紹簡單的求根方法，構(gòu)建只有一個自由度的系統(tǒng), 這里以相分率力作為未知變量。表4閃蒸單元的初始組分以及平衡時的分配系數(shù)K組分Xi65 C, 3.4bar下的

18、 Ki乙烷0.007916.2丙烷0.12815.2i-丁烷0.08492.6n-丁烷0.26901.98i-正戊烷0.05890.91n-正戊烷0.13610.72己烷0.31510.28液相碳氫化合物混合物通過閃蒸到 65C、3.4 bar狀態(tài)。液相組成和每種組 分閃蒸“ K”值見表4。我們希望知道閃蒸過程中氣相和液相產(chǎn)物的成分和液相 的比例。表4給出了原始成分。組分i的質(zhì)量平衡滿足以下關(guān)系(6)Xi =(1 -)yi 為Xi是閃蒸前液體的莫爾分數(shù)，Xi是閃蒸后液體莫爾分數(shù)，yi是蒸發(fā)產(chǎn)物莫爾分數(shù)， ©是液相產(chǎn)物和總供給摩爾流率的比值。平衡系數(shù)的定義為： Ki= yi/ X。將該

19、方程消去質(zhì)量平衡中的Xi,得到一個關(guān)于yi和Xi的方程:Xiyi =1、1 -* 1 -I Ki J由于各yi相加為1,所以我們有小的非線性方程:nf =1- i 1Xi1-* 1L0(8)其中n是組分數(shù)量。該函數(shù)f (電可用root(s)4求解，將計算結(jié)果回代就可以得到所有產(chǎn)物的莫爾分數(shù)。牛頓-拉夫遜方法需要知道當(dāng)前導(dǎo)數(shù)f'S)來決定計算方向，在COMSOL Multiphysics中將分析計算該值。通過下式很容易得到牛頓 -拉 夫遜迭代方法：Xi 11f' ； 1一Ki現(xiàn)在再回到COMSOL Multiphysics求根計算。作為一個練習(xí)，我們用一般 PDE 模式來計算求解

20、。我們可以只載入 roofinder.mph或rtfindgen.mph并加以修改來 計算這個問題，但是熟悉 COMSOL Multiphysics功能也是非常重要的目的。表5閃蒸實例Model Navigator選擇 1D,COMSOL Multiphysics: PDE modes:PDE, general 形式Draw菜單Specify objects: Line, Coordinates彈出菜單，x： 0 1 名稱：interval,完成Options 菜單： Constants輸入表4中數(shù)據(jù)Xi, 0.0079X2, 0.1281X3, 0.0849X4, 0.2690X5, 0.0

21、589X6, 0.1361X7, 0.3151Options 菜單：Scalar Expressions為(8)式中右側(cè)各項定義表送式t1-X1/(1-u*(1-1/K1)t2-X2/(1-u*(1-1/K2)t3-X3/(1-u*(1-1/K3)t4-X4/(1-u*(1-1/K4)t5-X5/(1-u*(1-1/K5)t6-X6/(1-u*(1-1/K6)t7-X7/(1-u*(1-1/K7)Physics 菜單：Boundary settings選擇域：1和2 (按住Ctrl鍵)選擇Neumann邊界條件 采用默認設(shè)置q=0, g=0完成Physics 菜單：Subdomain sett

22、ings選擇域：1設(shè)定 F=1+t+t2+t7; da =0選主舉Init選項卡，設(shè)定u(t0)=0.5Mesh菜單：Mesh parameters設(shè)定最大基元尺寸為1 點擊Remesh,完成Solve菜單：Solver parameters選擇Stationary nonlinear求解器求解，完成Post-processingPoint Evaluation選擇邊界：1,表遼式：u,元成Report窗口，值：0.458509,表達式：u啟動COMSOL Multiphysics,打開模型導(dǎo)航欄。如果你已經(jīng)打開 COMSOLMultiphysics對話框，那么將你的工作空間保存為 mph文件

23、或者將命令保存為 m文件，然后打開“ File”菜單選擇“New”。按照表5的步驟建立閃蒸模型。注意 本例中有兩個新增內(nèi)容" Options： Constants' 和 "Options： Expressions'。在這里可以定義常數(shù)，然后在COMSOL Multiphysics中任何可以使用純數(shù)字的輸 入框中使用定義的常數(shù)。表達式和常數(shù)類似，可以在任何COMSOL Multiphysics 數(shù)字輸入框中使用，但是不同之處在于表達式依賴于因變量。定義的常數(shù)和表達 式同樣可以在后處理過程中使用。后處理過程顯示t1和t2為：Value: -0.013865, E

24、xpression:t Boundary/Value: -0.203441, Expression: 2t, Boundary/另外求解器日志中經(jīng)常包含有用的信息。打開“solver”菜單，在最底部選擇“View Log”，會彈出求解器日志對話框。它顯示了求解器何時開始計算，執(zhí)練習(xí)行了什么命令，如何從初始狀態(tài)得到計算結(jié)果。 這里經(jīng)過三次迭代，絕對誤差達 到10 90如果你的解不是不收斂或收斂很慢的話，很容易得到這樣的信息。1.1 計算方程 f (u) = u3-3u2+5u-122的根。由于該函數(shù)是三次多項式， 以存在分析的無理數(shù)解，/1/1u =1 u =1= u=1+應(yīng)V2exp(u)-1

25、意味著不存在分析的無理1.2 計算方程f(u)=ueu -1=0的根。這是個超越方程,數(shù)解。如果你使用系數(shù)模式，設(shè) c=1來幫助收斂。3.方法2:步進法數(shù)值積分數(shù)值積分是數(shù)值分析的支柱。在有計算機之前，科學(xué)計算的首要工作是制作 特殊方程手冊，其中幾乎包括了所有特殊常微分方程的解。那么用什么方法計算 得到的？就是用一維數(shù)值積分。一維數(shù)值積分分為兩種：初值問題(IVP)和邊界值問題(BVP)。我們下一 節(jié)再介紹邊界值問題。最容易積分的是初值問題，因為所有的初始條件都在一點 定義，可以根據(jù)一階導(dǎo)數(shù)直接步進積分。顯然，如果常微分方程是一階的，例如：dy = f(t) y- = yt 十&f(t

26、)(10)dt口那么當(dāng)At-0時，(10)中第二項嚴格成立。它滿足歐拉法的條件，是積分一階常 微分方程最簡單的辦法。對于一維情況，可以簡單地根據(jù)f在(Xn,yn)點的導(dǎo)數(shù)向前步進積分，這里n是指積分過程中的第n次離散化步驟。所以，(11)yn 1 = yn hf (Xn , y)Xn 1 = Xn h這里假設(shè)導(dǎo)數(shù)變化不超過步進量 h,這只對線性函數(shù)才嚴格成立。對于任何 有曲率的函數(shù)，該假設(shè)都不成立。下面考慮在圖1中，如果步進量h取得過大會 造成什么錯誤。顯然，歐拉法需要改進的一點就是增大步進量， 因為它需要小步 進量以保證準確性。歐拉法也被稱作“一階”準確性，因為誤差只能降低到一階 h的程度。

27、圖1數(shù)值積分中的歐拉步進龍格庫塔方法所以如果我們希望增大步進量，那么就需要一個“更高階的方法”，隨著步進量的減小，誤差快速降低。k階方法誤差大約在hk量級。假設(shè)我們忽略了曲率， 可以通過計算斜率f(X)在Xn和Xn+1之間幾個中間點的值來計算y(X)的曲率?？梢?首先根據(jù)初始導(dǎo)數(shù)計算間隔中點的值，然后再用整個間隔中點的微分值來獲得二 階準確性。ki =hf(Xn，yn)./1,1, )k2 =hf Xn +2h, yn +2匕 |(12)3、yn 1 = yn k2 O(h )這樣我們通過對兩個函數(shù)的計算又得到一階準確性。所以，如果使用一階方法，N次計算得到O(1/N)誤差，而二階方法2N次計

28、算得到O(1/4N2)誤差，如果使用 一階方法的話，這需要N2次計算。更高階龍格庫塔方法我們可以做的更好嗎？顯然，如果用三個中間點可以獲得三階準確性，用四個中間點可以獲得四階準確性等等。 不過更高階的方法需要更多的編程工作，所以也要考慮我們的時間。但是從本質(zhì)上講，我們計算函數(shù)的k階導(dǎo)數(shù)不一定光滑。 可能在增加“近似準確性”的同時，高階導(dǎo)數(shù)的整體誤差會變得很大，如果是這 樣，每次連續(xù)步進都可能導(dǎo)致誤差的急劇增大。 這表明高階方法沒有低階方法穩(wěn)常微分方程通常選用四階龍格庫塔方法積分， 這樣程序比較緊湊，準確性高， 也有很好的穩(wěn)定性。其它方法在數(shù)值積分編程的時候，還需要注意兩個問題：數(shù)值穩(wěn)定性：即使

29、使用高階準確性的方法，你的積分結(jié)果和檢驗值偏差仍然 很大，這可能就是因為你的數(shù)值方法不穩(wěn)定。 你可以通過降低步進量來獲得數(shù)值 穩(wěn)定性，但是這也意味著更長的計算時間。如果你的計算量很大，計算時間太長， 可以選用半隱式方法，例如預(yù)測校正法等。剛性系統(tǒng)：物理機理起作用時，剛性系統(tǒng)通常會有兩個完全不同的長度或時 間尺度。剛性系統(tǒng)可能會出現(xiàn)上面提到的“系統(tǒng)非穩(wěn)定性”現(xiàn)象，或者非物理現(xiàn) 象的振蕩?？梢詤⒖糋ear2的書來解決這個問題。3.1 數(shù)值積分簡單實例高階常微分方程通常通過降階和步進法求解，考慮如下方程：dy q(X)dy=r(X)(13)dX dX除非q(X)和r(X)是常數(shù)，那么你很幸運能夠從大

30、多數(shù)分析方法教材中找到答案。也有一些特殊的q(X)和r(X)可以求得分析解，但是最好能計算出各種情況下的數(shù) 值解。為什么呢？因為為了搞清楚 y(X),你需要畫出曲線，得到分析解時你也需 要花費一定的精力在作圖上。那么應(yīng)該怎么做？首先將上面的二階系統(tǒng)降低為兩 個一階系統(tǒng)：dy dx= z(x)梟r(x)-q(x)z(x)上式的常微分方程能夠仿照(11)和(12)進行時間步進法數(shù)值積分。舉個簡單例子:d2udt2u =0(14)降階產(chǎn)生了兩個一階常微分方程:du1(15)=-u2 dtdu2二u1dt假設(shè)初始條件為U1 = 1, U2=0,可以建立零維空間系統(tǒng)來積分這一對常微分方程。打開COMSO

31、L Multiphysics,等待模型導(dǎo)航欄窗口。如果已經(jīng)打開一個 COMSOL Multiphysics窗口，把工作空間保存為 mph文件或者把命令保存為 m 文件，然后在“File”菜單中選擇“ New”選項。按照表6的步驟進行設(shè)置。表6時間步進積分的簡單實例Model Navigator選擇 1D, COMSOL Multiphysics: PDE modes: PDE, general 模式，插入因變重：Ui u2Draw菜單Specify objects: Line。Coordinates彈出菜單，x： 0 1 名稱：interval,完成Physics 菜單：Boundary se

32、ttings選擇域：1和2(按住Ctrl鍵) 選擇Neumann邊界條件保持默認選項q=0, g=0完成Physics 菜單：Subdomain settings選擇域：1設(shè)定 Fi = U2, F2 Ui選擇Init選項卡，設(shè)定U1 (tw)=1Mesh菜單：Mesh parameters設(shè)置最大基元尺寸為1 點擊Remesh,完成Solve菜單：Solver parameters選擇time dependent求解器，在 General選項卡中輸入 Times: linspace(0,2*pi,50)求解，完成Post-processing Cross-section plot param

33、etersPoint選項卡，接受ui的默認設(shè)置 General選項卡，接受所有時間的默認設(shè)置 完成這個例子給了我們兩個獨立變量ui、u2和空間坐標x。注意“Subdomainsettings窗口左上角方程中的矢量標志。由于是矢量變量，所有輸入選項卡都是 矢量(F)或者矩陣(da)輸入。MATLAB 命令linspace(0,2*pi,50)生成50個基元的矢量，給出從 0到2九的 數(shù)據(jù)，u1(t=24=1.004414。假設(shè)分析解是u1(t=24=1,結(jié)果不夠準確(0.4%誤差)。 之前FEMLAB允許用戶在MATLAB和FEMLAB中選擇針對時間積分的預(yù)置的求解器。也允許用戶在求解器變量對話

34、框的通用選項卡中調(diào)整容錯值 （相對和絕 對）。將相對誤差調(diào)整為0.001,絕對誤差調(diào)整為0.0001,得到一個相對較好的計 算結(jié)果ui（t=2力=0.998027。注意累積的全局誤差和求解方法準確性 0.001的量級 一致。圖2 一個周期的ui（t）u2圖3 一個周期的U2（t）圖2和圖3表明如果設(shè)定足夠小的步進時間，F(xiàn)EMLAB能夠完成高精度的 正弦和余弦函數(shù)數(shù)值積分。雖然我們一般認為正弦和余弦是“分析函數(shù)”，但是通過比較，還是能清楚地看到分析函數(shù)和需要數(shù)值積分的函數(shù)是似是而非的一一 沒有比Bessel函數(shù)和橢圓方程更有分析性的函數(shù)。練習(xí)1.3使用物理場菜單中的常微分方程設(shè)置和相同的初時條件

35、，試著積分方程（15）,假設(shè)狀態(tài)變量為V1和V2o基于時間的常微分方程和代數(shù)方程的唯一區(qū)別就是 必須使用符號代替dvi/dt和vit等。4.管式反應(yīng)器數(shù)值積分本節(jié)介紹在化工管式反應(yīng)器設(shè)計時，如何同時求解一階非線性常微分方程 組。通常管式反應(yīng)器設(shè)計的關(guān)鍵要素是反應(yīng)器的長度。一個氣態(tài)酒精脫水管式反應(yīng)器，工作在 2bar和150c下。反應(yīng)方程式為：C2H50H，C2H 4 H 2O已有試驗表明反應(yīng)速率可以表示如下，其中Ca是酒精濃度（mol/L） , R是酒精的消耗速率（mol/s/m3）：252.7CA彳 0.0131 -Ca反應(yīng)器直徑為0.05m,入口酒精流速為10g/s。問：如何確定反應(yīng)器的長

36、度， 進而獲得各種酒精轉(zhuǎn)化率？我們希望確定酒精出口摩爾分數(shù)分別為0.5, 0.4,0.3,0.2和0.1時的反應(yīng)器長度。化工設(shè)計理論假設(shè)反應(yīng)熱很小，反應(yīng)器內(nèi)為柱狀流，理想氣體，可以確定反應(yīng)流體由四個 常微分方程控制，其中獨立變量為 Ca, Cw （水濃度），V （速度）和x （反應(yīng)器 長度）:dCACa二一R 11 出.CdCW =R 1 -CW(16)dt . CdV RVdt Cdx 、/二V dt最后一個方程說明表面速度使得反應(yīng)器長度和反應(yīng)物在反應(yīng)器內(nèi)的存在時 間之間存在一個平衡。初始入口條件為：Ca（0）=C V（0）=VoCw（0）=0x（0）=0方法從初時條件和化學(xué)計量系數(shù)可以看出

37、，Cw= Ce （酒精濃度），由于溫度和壓力設(shè)為常數(shù)，所以濃度C也為常數(shù)，可以看作是理想氣體，有:pT 18314J一 kmolK(18)初始進口速度Vo可以由給定流速、入口密度（酒精分子量為 46kg/kmol）和管的 橫截面積計算出。該方程需要對時間t數(shù)值積分到需要的酒精摩爾分數(shù)。使用較 為簡單的歐拉方法或者龍格庫塔方法積分。讀者也許注意到，積分Ca可能不需要其它變量，但是 Cw, V和x和時間t 嚴格相關(guān)。但是由于我們需要求出對應(yīng)摩爾分數(shù)的 x值，所以最好同時求解四個 變量。部分結(jié)果計算得到的數(shù)值解為：CA=0.1Ct =5.65225(19)其中Ca/C如圖4所示x = 18.5435

38、COMSOL Multiphysics 計算我們希望再次建立虛擬的零維模擬空間，這次使用四個獨立變量。打開 COMSOL Multiphysics,等待模型導(dǎo)航欄窗口。根據(jù)表 7的步驟建立管式反應(yīng)器 模型。該模型給出四個獨立變量 U1、U2、U3、U4和一個空間坐標x。由于變量比 較多，建議用變量名稱進行常數(shù)設(shè)定。進一步來說，由于使用了速度定律，用矢 量表達式更為方便。表7 COMSOL Multiphysics中管式反應(yīng)器設(shè)計建模步驟Model Navigator選擇 1D 空間維數(shù)，COMSOL Multiphysics: PDE modes: PDE, general形式設(shè)止獨立變量：U

39、i u2 u3 u4選擇 Element： Lagrange-Linear,完成Draw菜單Specify objects: Line, Coordinates彈出菜單，x： 0 1名稱：interval,完成Options 菜單： Constants如下填寫變量表格： 名稱ExpressionP200000T423R8314MM46Flowrate0.01Dia0.05CP/(RT)areapi*DiaA2/4rhoMM*CvelFlowrate/rho/aera完成Options 菜單：Scalar Expressions定義 rate=52.7*U1A2/(1+0.013/u1)Phys

40、ics 菜單：Boundary settings選擇域：1和2 (按住Ctrl鍵)選擇Neumann邊界條件默認設(shè)置q = 0 g=0,完成Physics 菜單：Subdomain settings選擇域1F 選項卡；設(shè)定 F1=-rate*(1 + u1/C); F2=rate*(1-U2/C)；F3=rate*u3/C; F4=u3/CInit 選項卡；設(shè)定 U1(t0)=C; U3(t0)=vel,完成Mesh菜單：Mesh parameters設(shè)定最大基元尺寸為1 點擊Remesh,完成Solve菜單：Solver parameters選擇time dependent求解器，在Gene

41、ral選項卡中輸入Times: linspace(0,10,100)求解，完成Post-processing Cross-section plot parametersPoint選項卡，接受U1的默認設(shè)置 General選項卡，接受所有時間的默認設(shè)置 完成試著在整個時間范圍內(nèi)畫出U1, U2, U3和U4的曲線，和圖4是否保持一致?和之前計算結(jié)果是否相符？練習(xí)1.4計算以下方程組中y'(x=1)的值，并畫出y隨x變化曲線，x取0至|3區(qū)間2 _y' y' y = 0y(x =0) =1y'(x = 0) =01.5連續(xù)攪拌反應(yīng)器中的一階可逆反應(yīng)系統(tǒng)可以用線性常微

42、分方程組表示。例 如，考慮如下異構(gòu)反應(yīng)：A B- C正向反應(yīng)速率為k1和k3,相應(yīng)逆向反應(yīng)速率為k2和k40 一階動力學(xué)常微分方程組為:dcA"df = -k1cA ' k2 cB(20)dCBk1cA -'k2 cB -"k3cB k4cC dtdecdt-k3cB - k4cC你也許會感到奇怪，由于該方程組為線性，它有通用的分析解。雖然是通用 分析解，但是對動力學(xué)系統(tǒng)的反映并不多。假設(shè)剛開始時，Ca=1, ki=1Hz,k2=0Hz, k3=2Hz, k4=3Hz。針對該初值問題，畫出濃度隨時間變化的曲線。5.方法3:常微分方程數(shù)值積分前一節(jié)介紹了步進法

43、數(shù)值積分，通常也稱作“時間步進”，但是在反應(yīng)器設(shè) 計中，多是空間積分問題。步進法是順序求解未知項，而其它常見積分方法是同 時求解某一網(wǎng)格點的所有未知獨立變量。步進法只能求解初值問題( IVP),初 值個數(shù)必須嚴格等于常微分方程組的階數(shù)。 但是對于二階或者更高階的系統(tǒng)，可 能會用到第二類邊界條件一一邊界值問題， 在一維情況下，這些邊界條件只出現(xiàn) 在域的起點和終點，這就是兩點邊界值問題。步進法也能夠勉強求解邊界值問題 能一人為設(shè)定初值，通過嘗試和誤差修正找到滿足邊界值問題的初始條件。對于 更高維數(shù)的偏微分方程，邊界值問題發(fā)生在整個域的每一個邊界上。有限元方法的一個最主要優(yōu)點就是能夠很容易地求解兩點

44、邊界值問題。例 如，一維反應(yīng)和擴散方程：2=R(u)(21)D二 u22-2L二 x這里u是組分濃度，D是擴散系數(shù)，L是域的長度，R(u)是反應(yīng)消耗速率，x是 無量綱空間坐標。如果未知函數(shù) 口僅)在x=xj=jzx網(wǎng)格點能夠用不連續(xù)值Uj=U(Xj)L02近似，那么根據(jù)中心差分，方程變?yōu)椋?22)N、MMj i其中Rj=R(Uj), Mj是對角元素為-2,上、下對角元素為1的三對角矩陣：-2 i 0 0 tin 1-210 in M = 01 工工工(23)| 00工工工+FJRj =R(Uj)。通過矩陣轉(zhuǎn)置和對uin進行積分可以求解該方程，這里n指第n次猜測:nUi2 2 R2 . NMRj

45、D j I(24)Rj=R（Uin-1）。無論是初值問題還是邊界值問題，都可以將（23）式中矩陣M的行向量變得符合邊界條件。（23）式假設(shè)在乂=0和乂= 1處都有u=0。這是狄利克雷邊 界條件，也是有限微分法的自然邊界條件一一說它自然是因為如果在設(shè)定邊界條 件時沒有改變（23）式中的行向量，那么就會出現(xiàn)這樣的邊界條件。下面介紹如何用COMSOL Multiphysics在一維域上求解方程（21）,對于一階 反應(yīng)R（u）=ku,典型的無量綱變量 Damkohler數(shù)為：DakL2D23 . 13. 110 s 10 m1.2 10-9m2/s=0.833(25)邊界條件為x=0時U=1,在U=1

46、處沒有流量。下面結(jié)合MATLAB來做這個練習(xí)，進而了解 COMSOL Multiphysics如何表 示計算數(shù)據(jù)和模型輸出的結(jié)構(gòu)。在windows中，結(jié)合 MATLAB的COMSOLMultiphysics是一個可選的桌面圖標，在 UNIX/linux中，該功能可以通過在命令 欄中執(zhí)行以下語句實現(xiàn)：Comsol matlab path-ml nodesktop-ml nosplash其中“matlab”告訴femlab打開一個 matlab命令窗口，“path”根據(jù) COMSOL 命令庫建立 matlab命令窗口。首先運行COMSOL Multiphysics和模型導(dǎo)航欄窗口，然后根據(jù)表8中的

47、步驟 建立模型。該實例給了一個獨立變量 u和一個一維空間坐標x。h和r是系數(shù)模 式中狄利克雷邊界條件的兩個句柄。如果你希望邊界速度 u為常數(shù)Uo,可以通 過設(shè)定h=1, r = U。來實現(xiàn)。點擊三角形按鈕生成默認網(wǎng)格（15個），然后點擊 兩個嵌套三角形按鈕加密網(wǎng)格（30個）。表8反應(yīng)-擴散系統(tǒng)中兩點邊界值問題的常微分方程實例Model Navigator選擇 1D 空間維數(shù)，COMSOL Multiphysics: PDE modes: PDE, coefficient 形式 設(shè)置因變量：u選擇 Element: Lagrange-Linear,完成Draw菜單Specify objects:

48、 Line, Coordinates彈出菜單，x： 0 1 名稱：interval,完成Physics 菜單：Boundary settings選擇域1選擇Dirichlet邊界條件，設(shè)定h=1; r=1選擇域2選才N Neumann邊界條件，完成Physics 菜單：Subdomain settings選擇域1設(shè)定 C=-1; f=0.833*u; da =0選擇Init選項卡；設(shè)定u（t0）=1-x 完成Meshing點擊三角形按鈕進行網(wǎng)格劃分Solve菜單：Solver parametersStationary linear求解器默認設(shè)置，元成General選項卡，設(shè)止解的形式為Coef

49、ficient求解(=)Post-processingPoint evalution選擇邊界2,輸入表式u。Reports窗口顯小： 值：0.861167,表達式：u,邊界：2根據(jù)計算結(jié)果可以畫出圖5,顯然邊界條件需要滿足：x=0處u=1,在x =1處斜率減為零。但是COMSOL Multiphysics有沒有按照我們設(shè)想的求解呢？Line: u圖5穩(wěn)態(tài)下的濃度曲線現(xiàn)在用穩(wěn)態(tài)非線性求解器再計算一次。 首先查看求解器菜單中的日志，迭代 13次收斂。你注意到最終值從0.86降至I 0.69 了嗎？為什么會這樣呢？線性(系 數(shù)型)求解器只在初始時刻u(t0)=1-x時計算一次R(u),所以方程(24

50、)只需要一次 迭代。而非線性求解器在每次迭代時計算一次R(u),計算時使用上次迭代的 u值。所以非線性求解器在向收斂靠近時，經(jīng)過幾次迭代可能完全“忘記” 了初時 猜測。檢驗一下這種解釋對不對。改變初值將會改變穩(wěn)定的線性解。 返回子域設(shè)置 對話框，設(shè)定初值為u(t0)=1,最終得到的u(x=1)是多少？再試一下穩(wěn)態(tài)非線性求 解器，能得到和其它初值一樣的結(jié)果嗎？該例說明為方程選擇正確求解器的重要性。如果函數(shù)f依賴于任何非獨立變 量，就應(yīng)該選擇非線性求解器。線性求解器更快，但是它假設(shè)偏微分方程的系數(shù) 不依賴于非獨立變量 u (否則為非線性問題)。如果不確定，還是應(yīng)該用非線性 求解器。別忘了，(21)

51、式滿足R(u)=ku,是線性問題，但是COMSOL Multiphysics 使用非線性求解器才得到正確解！收斂比較慢也是模型形式導(dǎo)致的結(jié)果一一選中 雅克比求解器選項，非線性求解器只需要兩次迭代就可以得出正確結(jié)果。我們一開始認為(22)式是該問題有限微分矩陣方程，但是實際上最后用(24)式描述COMSOL Multiphysics的有限元問題。因為我們這里使用拉格朗日線性微 元，在這種特殊情況下有限元和有限微分矩陣一致。根據(jù)這一點，我們看一下MATLAB/COMSOL Script 如何計算該 COMSOL Multiphysics 問題。打開“File”菜單，選擇“Export FEM St

52、ructure as femi'。這就把現(xiàn)在的計算 結(jié)果轉(zhuǎn)化為MATLAB/COMSOL Script工作空間的數(shù)據(jù)結(jié)構(gòu)，然后用 MATLAB/COMSOL Script預(yù)置的函數(shù)和命令來計算它。在MATLAB/COMSOL Script工作空間輸入以下命令：>>x=fem.mesh.p;>>u=fem.sol.u;>>plot(x,u)將會彈出一個包含網(wǎng)格節(jié)點 u值的MATLAB/COMSOL Script圖形窗口。不 要懷疑你的圖形失真。這是因為COMSOL Multiphysics存儲網(wǎng)格點和對應(yīng)結(jié)果的 方式是為了使矩陣變得稀疏和緊湊。可以用

53、MATLAB排序命令對網(wǎng)格點進行排 序，進而使圖形有意義。在MATLAB工作空間輸入以下命令：>>xx, idx=sort(x);>>plot (xx, u(idx)最后考慮MATLAB對有限元矩陣的訪問。Fem文件結(jié)構(gòu)沒有包含有限元剛 度矩陣，但是它包含重建這個矩陣的 FEMLAB函數(shù)的必要信息。這對有限元方 法非常重要，這個FEMLAB函數(shù)是assemble，輸入以下命令：>>K,L,M,N= assemble(fem);>>K' /15這張圖包含上次COMSOL Multiphysics的計算結(jié)果，類似圖5。實際上，我 們之所以能夠

54、快速搞清楚fem求解結(jié)果的結(jié)構(gòu)，是因為這是只有一個獨立變量的 一維問題。多變量和多維產(chǎn)生的網(wǎng)格和結(jié)果的數(shù)據(jù)結(jié)構(gòu)只能用COMSOLMultiphysics的工具/函數(shù)來讀取。圖6給出了由MATLAB命令spy(K')產(chǎn)生的松 散結(jié)構(gòu)。很有必要將它和(23)式相互比較，因為能夠清楚地看到具有統(tǒng)一網(wǎng)格劃 分的一維拉格朗日線性基元和生成矩陣方程的有限微分法是類似的。D2468101214160246010121416nz = 46圖6由MATLA命令spy( K')產(chǎn)生的K'松散結(jié)構(gòu)下面看一下MATLAB中矩陣的松散結(jié)構(gòu)，所有的元素只有 1, -2和1,位 于三條對角線上。這是

55、有限元法的剛度矩陣，對于未知類型，它等于(23)式。如果返回“Subdomain settings'對話框的“element”選項卡，選擇拉格朗日二次元 素，再次求解并導(dǎo)出FEM,生成上面的K,你會發(fā)現(xiàn)雖然矩陣也很疏松，但是 和拉格朗日線性元素完全不同。練習(xí)1.6偏微分方程的系數(shù)形式有一項 刈，令a= 0.833, f=0,再次計算以上“反應(yīng) -擴散”的例子，比較穩(wěn)態(tài)線性和非線性求解器得出的結(jié)果。你能解釋為什 么這樣的表達式會導(dǎo)致這樣的結(jié)果嗎？這樣的表達式對剛度矩陣K又有什么影響？如果把Da選為某特定值，例如 DaAx2 = 2,使得對角元素幾乎為 認，你認為求解時會出現(xiàn)什么困難？6.方法4:線性系統(tǒng)分析MATLAB和COMSOL Multiphysics的核心是線性系統(tǒng)分析。本節(jié)將簡要回 顧線性算法理論一一本科生工程數(shù)學(xué)中的“矩陣方程”就是這類典型問題。好消息是不需要你自己動手編程處理矩陣，這就是MATLAB的用途：提供一個工程矩陣計算子程序庫的用戶界面。應(yīng)當(dāng)注意到，對于本章的例子，COMSOL Script也能很好地實現(xiàn)這個功能。以前科學(xué)計算主要集中在矩陣計算的高效和松散方法 上。Golub和Van Loan3的書是一本很好的專家級矩陣計算指南，但是對于 MATLAB入門水平，Hanselman和Littlefield4的書是個不錯的選擇。簡單來