變量耦合反問題線積分積分方程和積分微分方程

1、第七章 變量耦合：反問題、線積分、積分方程和積分微分方程W. B. J. ZIMMERMANDepartment of Chemical and Process Engineering, University of Sheffield,Newcastle Street, Sheffield S1 3JD United KingdomE-mail:本章進一步深入探索耦合變量，考慮如何利用耦合變量求解反演方程和各種 類型的積分方程。本章選擇了四種重要的應用實例一一使用激光探測位置，污染云層的傳播，電容層析成像反演問題，纖維媒介中的非局部傳熱以及顆粒處理中 的聚集平衡方程。我們將用到之前在COMSO

2、L Multiphysics中沒有用到的新特征 本一通過MATLAB與優(yōu)化工具箱耦合，擴展網格，使用瞬態(tài)求解器作為穩(wěn)態(tài)非 線性模型的迭代工具，以及在耦合模式下選擇激活/非激活多物理場模式的能力。 后者對于單一耦合問題非常有用（如第3章中介紹的球形顆粒表面催化劑周圍的 水動力學問題）。為了處理本章中的積分方程，設置輔助域上的虛擬因變量。耦 合變量在該域進行各種運算，但是因變量本身并不需要。因此，可以更好的調節(jié) 對積分方程的有限元近似。1 .引言我們已經對邊界和域積分較為熟悉了可以在COMSOL Multiphysics的后處理菜單中選擇。邊界積分對計算各種表面量非常有用：例如靜電體表面電荷， 流

3、體中的曳力等。域積分常用于平均和在域中定義的較高程度自由度的結合。這些特征是可信的，并給出了通過積分形式表達的有限元方法，伽遼金方法（見第2章），COMSOL Multiphysics包含了有效和精確的積分算法。因此，如果讀者對 任意積分階數(shù)的數(shù)值積分或者常微分方程組 （ODES）感興趣，可以查閱MATLAB 的內置函數(shù)（見第1章），這里不再贅述。在本章中將討論積分方程及其理論的 更多復雜應用，均采用COMSOL Multiphysics求解。線積分、積分方程以及積分 微分方程是應用數(shù)學研究的對象，也是COMSOL Multiphysics之所以具備變量 耦合能力的基礎。這也是為什么本章要取名

4、為“多物理場耦合問題的伸展” 。盡 管題目選定了，但這只是說明性的，我們仍然以解決化工過程中的問題為目標， 對COMSOL Multiphysics的特征進行深入了解。這里最重要的是使用激光雷達定 位以及氣體污染云層的擴散，聚集平衡方程（通常為微分方程的初值問題） ，電 容層析的反演問題，以及纖維媒介中的非局部傳熱問題的計算。擴展多物理場在第一個方法中，我必須承認自己曾對我們要研究的擴展的多物理場問題保 持懷疑的態(tài)度。最終，標量耦合變量吸引了我，并提供了撰寫第4章的素材。它也促使我的課題組對微流網絡這樣一個新的領域產生了興趣。目前， COMSOL Multiphysics提供了其他兩種耦合變量

5、的概念拉伸和投影耦合變量。關于使用 他們的例子均與后處理有關，例如，用于在域的截面函數(shù)性的表達求解結果以分 析特定特征。盡管如此，在一些較少出現(xiàn)的情況下，耦合變量（拉伸和投影）會應用到模 型中與獨立域變量同時求解。多個域、多尺度以及多過程模型對于工程數(shù)學和數(shù) 學物理并不常見。通常情況下，模型是局部的-包含于一個特定的微分或偏微分方程組中，而方程組的邊界條件和初始條件均較好的給出了定義。這些是局部連 續(xù)性模型。從歷史觀點來說，已經從分析技術使用的角度分析了這方面的發(fā)展， 在一些領域上已經涉及閉式模型的種類。即使是能夠采用連續(xù)方法處理的模型， 仍然需要根據(jù)離散方法進行近似。連續(xù)性顆粒水動力學1和離

6、散方法2越來越受到關注，但是較老的方法，如分子動力學模擬3,蒙特卡羅方法4,微觀水動力學5,分區(qū)自動化6和氣體/等離子體的精確數(shù)值解 口為連續(xù)系統(tǒng)/離散系統(tǒng) 建立了橋梁。另一個技術策略基于滿足偏微分方程約束的優(yōu)化理論-罰不可行自由度?；旌戏e分非線性規(guī)劃8,遺傳算法9以及遺傳規(guī)劃10都對混合離散/連 續(xù)系統(tǒng)較為適用。COMSOL Multiphysics傾向于以PDE約束下的連續(xù)性模型進 行方程描述。因此，從概念上講，物理場伸展并不是處理較粗糙時的彌補性思考。 它可以允許從基礎上處理離散系統(tǒng)， 而且不僅僅是偏微分方程約束的問題， 而且 也可以處理積分約束問題。實際上，耦合變量允許非局部和離散模擬

7、。在第3節(jié)(Scalar),第4節(jié)(projection)以及第5節(jié)(extrusion)中，我們 重新回顧了耦合變量來研究如何使用 COMSOL Multiphysics處理反演問題，線積 分以及積分/積分-微分方程問題。標量耦合變量毫無疑問，標量耦合變量從概念上最容易理解。 在第4章中，標量耦合變量 用來連接非均相反應器中的冉循環(huán)流體-反應器的輸出，通過合理改變比例，重新進入原料流中。建立了零維單元和第二個幾何體用于模擬緩沖容器，該容器獲得了循環(huán)流體與反應器出口之間的代數(shù)聯(lián)系。一個標量耦合向目標域、子域、壁面或邊線傳遞同一個量，在模型中任何一個地方都不發(fā)生改變。標量耦合變量由 一個源域中的

8、積分創(chuàng)建。在我們的例子中，源是零維的(末點或者單個單元結構), 積分與被積函數(shù)相同。更進一步說，緩沖容器模型是人為建立的，其再循環(huán)關系 能夠由弱邊界條件約束來處理而不需要考慮第二個域。因此，我們已經看到標量耦合變量源積分并不常見。下一章中，我們會著手解決反演問題，這里耦合將是 必要的而且錯綜復雜。反演問題意味著存在一個描述較為清楚的正向問題，但是反演問題卻沒有描述清楚。我們選擇的反演問題是一個電容層析成像問題。電容過程層析成像目前已經成為一種較成熟的工程科學技術。最常見的結構即電容層析成像，經常用于圓柱管道內的多相流圖像處理。測定多相流體中分散相的 細節(jié)信息對于化工單元的過程控制非常重要。兩相

9、流的非嵌入和非侵入測量非常 難以獲得。而油氣生產行業(yè)內經常碰到的高時間變動性流體會進一步增加測量的 困難程度。對流體瞬態(tài)特性的精確測量以及相分布的實時分布情況難以獲得。一個可能的解決方法是使用電容層析成像技術測量流體的介電常數(shù)分布，測得數(shù)據(jù)可以給出管道截面上相分布的情況。層析成像設備能夠提供容器或者管道中的非嵌入性測試的組分分布圖像.電容層析成像技術提供了管道中組分電介質分布的二維圖像.通過電極實現(xiàn)非嵌入式電容測試時，電極通常由充-放電理論11進行激勵，然后使用數(shù)學重構算法來 建立具有不同介電常數(shù)的材料圖像。這個過程可以確定不同相.至今為止，許多關于過程工程的研究采用了 ECT方法，但是一些應

10、用涉及了流化床和氣動運輸。 McKee等12報道了電容層析成像技術在兩個工業(yè)裝置尺度上氣動傳輸圖像方 面的應用。該項工作開創(chuàng)了 ECT技術在密相氣動傳輸方面的應用，并展示了電 容層析成像技術在在線過程控制方面的發(fā)展?jié)摿?。文獻 13中對該領域的工作做 了較好的回顧。層析成像裝置包括三個主要部件：傳感器陣列（通常采用12電極；66獨立測試），數(shù)據(jù)采集系統(tǒng)和圖像重構系統(tǒng)。電容測量從電極中獲得，但電極有可能 被污染。對于每個電極對，均進行如下的充-放電操作：工作電極在給定電壓下 充電（15V）,測試電極接地，同時連接電流測試裝置。該測試裝置給出振蕩電 流的平均數(shù)值，獲得與未知電容值成一定比例的電壓?；?/p>

11、本的電容數(shù)據(jù)采集系統(tǒng)基于電荷傳遞定律。 放電電流產生了一個負的輸出 電壓。典型的充/放電循環(huán)以1MHZ頻率重復，連續(xù)的充電放電電流脈沖在兩個 電極之間進行平均，產生兩個直流輸出電壓。在電極布置之前需進行儀器校準，傳感器裝置由較低介電常數(shù)物質填充。該 步驟可以提供介電常數(shù)的參考值。然后當管道中充滿較高介電常數(shù)物質時改變電 流測試的敏感性。校正步驟對于所研究的任何一類材料均是必要的。通過使用幻燈算法進行圖像重構，可獲得管道中濃度分布的圖像?，F(xiàn)有的ECI算法技術能夠以每秒100幀的速度生成圖像，幾乎可以獲得過程中的實時信 息。盡管如此，現(xiàn)有的ECT系統(tǒng)分辨率較低（1/10管道半徑）。造成該局限的主

12、要原因是每個電極的表面積較大， 對于所有的實際目的，充放電時電極之間的電 場線是平行的，這樣可使圖像的重構過程變得較為簡單。如果使用更多或者更小 的電極，則有希望獲得較高的分辨率，但是會使重構算法變得更為復雜。 這種情 況下將不得不求解泊松方程以獲得內部介電常數(shù)場。這部分內容里，我們將給出在稀疏系統(tǒng)、大面積電極以及圖像重構反演問題， 數(shù)據(jù)采集為12通道（網格劃分如圖1所示）。（a） 網格劃分（b）場力線的正向解圖1左側：管道包含四個桿狀物時的矩形橫截面網格劃分情況， 每根桿的介 電常數(shù)為ex=gy=e z= £ 4=0.005 ,中間介質具有單位介電常數(shù)e。右側：當西側 電極保持單位

13、電壓，其余電極和邊界保持接地時的穩(wěn)態(tài)等電勢（電壓）線。一日T二二二二二二日管道中的電荷密度通過對無磁場耦合的麥克斯韋方程的合理簡化與電勢進 行關聯(lián)14:(1)其中，為單位體積中的電荷量，在體相流體和內含物中為零，但是在電極表面不為零，e為介電常數(shù)或者是介質的介電常數(shù)，由縮放因子的選擇而決定 為電勢（電壓）。使用方程（1）,應用較小的控制體積考慮電極和體相流體截面,(2)可得到電流通量邊界條件:n n其中，s為固體組元的介電常數(shù)，w為體相流體介質介電常數(shù)。LHS為由電極一側流出的電流通量，穩(wěn)態(tài)下的電位差由電極積累的電荷平衡。整理方程（2）可得邊界條件為：一q（3）n 0 s其中，我們令q為電極電

14、荷量。有了這些控制方程，我們能夠定義相關的層析成像數(shù)學問題。1.1 正問題假設電極電壓保持在穩(wěn)定電壓（見圖1）,傳感電極保持接地（0V）。方程（1） 的解計算了電極上總的電荷為：q；i d（4）n在內含物介電常數(shù)已知的情況下，該問題是一個正問題。圖1 （右）為該問題的解，我們將在 COMSOL Multiphysics中表示此問題。反演問題我們模擬的最終目的在于解決一個反演問題-根據(jù)測量得到的輸出結果計算未知的參數(shù)。為了事先這個目標，我們首先設置正問題，然后使用正問題來尋 找反演算法。對于未知參數(shù)（如同在高中代數(shù)中所學），我們很少能夠使用正向 模型來求解反演問題。更為復雜的是如何定義一個待求的

15、反演問題。例如，假設 進行同樣的實驗，但是管道中的電介質場是未知的。 電荷q可在電極上測試，介 電常數(shù)e與通過（4）測得的值保持一致。這個問題可劃分為反演問題，下一小節(jié)中， 我們將使用正問題模型來獲得反演問題的解。在COMSOL Multiphysics模擬ECT的正問題啟動COMSOL Multiphysics以及模型導航器，按照表 4中的說明設置ECT 的正問題。方程（1）描述的偏微分方程組，以及相對應的邊界條件，這個問題是 線性的、穩(wěn)定的嗎？我們是否應該使用穩(wěn)態(tài)線性求解器呢？而且我們可能需要一 個瞬態(tài)求解器，因此，確定da系數(shù)現(xiàn)在較為重要。非線性求解器考慮了兩個牛 頓迭代步驟，在實際應用

16、中可獲得比線性求解器更準確的解。如果你不相信我， 可以檢查求解器記錄。在這個問題中，第一個牛頓迭代步達到的誤差估計為 10-14, 第二個牛頓迭代步則可達到 10-16。但是，使用非線性求解器還有另外一個較為 細微的原因-為了保證方程組所有的解都得到了識別。如果一些相互關聯(lián)關系隱藏于通用PDE方程的“F”項中，那么隨后的牛頓迭代可以將這些因素分離出來。 這與線性求解器并無較大差別，盡管如此，有的時候線性求解器可以得出一結果， 但非線性求解器可能不收斂。表1,靜電場線模型。文件名：模型導航欄選擇二維模型選擇 COMSOL Multiphysics| PDE Modes| General Mode

17、 因變量名稱：phiOptions 菜單 | Axes/Settings設定x軸為一1到1;設定x和y的柵格間距均為0.1Draw菜單設定矩形R,左下角(-0.25,1),高0.1,寬0.5;設定矩形R2,左下角(1,-0.25),高0.5,寬0.1;設定矩形R3,左下角(-0.25,-1.1),高0.1,寬0.5;設定矩形R4,左下角(-1.1,-0.25),高0.5,寬0.1;設定止方形SQ1,左下角(-1,-1),寬2;設定圓形，圓心(0.6,0.2),半徑0.1;設定圓形,圓心(-0.4,0.7),半徑0.1;設定圓形，圓心(0.4,-0.5),半徑0.1;設定圓形，圓心(-0.8,-

18、0.2),半徑0.1;選擇 Create Composite object 輸入公式SQ1+R1+R2+R3+R4；取消選中 “ keep internal boundary"Options 菜單 |Constants 選項e0e1e2e3e410.050.050.050.05phi1phi2phi3phi41000Physics 菜單 |Boundary settings設定邊界 4, 5, 6, 7, 13, 15, 16, 18保持 Dirichelt 邊界 條件（R= phi）；設定邊界1, 2, 3為Dirichlet邊界條件，R=phi1-phi設定邊界8, 9, 12為

19、Dirichlet邊界條件,R=phi2-phi設定邊界 10, 11, 14為 Dirichlet 邊界條件,R=phi3-phi設定邊界 17, 19, 20 為 Dirichlet 邊界條件,R=phi4-phiPhysics 菜單 |Subdomain settings子域_rFdaea1e0*phix e0*phi y0102e2*phix e2*phi y0103e3*phix e2*phi y0104e4*phix e2*phi y0105e1*phix e1*phi y010Mesh標準網格加密一次Solve菜單設定穩(wěn)態(tài)非線性求解器 點擊工具欄計算按鈕(=)表2.每個電極的邊界

20、條件邊界被積函數(shù)Q1, 2, 3nx*phix+ny*phiy2.0500658, 9, 12nx*phix+ny*phiy-0.03068910, 11, 14nx*phix+ny*phiy:-0.033562117, 19, 20nx*phix+ny*phiy-0.015282表2包含了正交通量的邊界積分注意（nx,ny）是COMSOL中的保留幾何變量。 這些是外部點在邊界上的正交矢量分量，可根據(jù)標準方程計算：一?（5）n在靜電問題中，可解釋為正交電場分量。值得注意的是在靜電問題和穩(wěn)態(tài)傳 熱問題之間存在著精確的比擬關系。 可以理解為溫度，通量q為熱流通量，介 質為棒狀。因此，正問題模型描述

21、了由于在特定邊界施加電壓而引起的通量，而通過電極的場泄漏是可測的。因此，反演問題即通過測量場泄漏來確定域中的介 質性質。例如，我們能夠知道體相介質的介電常數(shù)以及棒的位置嗎？主要的問題 在于測量，比如，通量q對于參數(shù)ej非常敏感。一些情況下，反演問題可以由 COMSOL Multiphysics圖形用戶界面求解。耦 合變量對于GUI反演方法起著重要作用。最簡單的問題是通過電極1通量q,邊 界1,2,3被測量，并用于推斷子域2中的介電常數(shù)及。我們可以想象通量q1對于 e2的變化較為敏感，我們能夠將這樣的敏感性顯示出來嗎？這個問題可以使用參 數(shù)求解器來處理。按照表3中的步驟創(chuàng)建一個積分變量來計算 q

22、1,使用參數(shù)求解 器顯示e2的敏感性。圖2顯示了敏感性研究的結果。首先，好的消息是 q1隨e2 的變化是單調的，因此通過反演問題來確定改的值的解是單一的。但不好的消息是這樣的敏感性水平（2從0.05變化到0.2時，q1改變1%）使得鑒別實驗噪 聲變得較為困難。而且也會使其他參數(shù)的確定變得較為困難。 注意，模型的自由 度為21625。表3前面ECT模型的變量敏感性修正。文件名：ectrect模型導航欄打開Options 菜單 | Integration coupling Variables| Boundary variables選擇邊界1, 2, 3名稱：q1表式：nx*phi x+ny*phi

23、y保持全局目標Solve菜單選擇非線性變量求解器，General選項卡 變量名稱：e2變量化 0.01:0.01:0.2點擊工具欄求解按鈕（=）Post-processin球單選擇橫截面變量繪制Point選項卡，設定表達式 q1/2.050065General選項卡，選擇 point plot圖2常規(guī)通量qi隨e2變化的變量敏感性圖1.2 在GUI中使用反演方法現(xiàn)在讓我們使用GUI來求解反演問題。假設我們測得 qi=2.07,想要知道介 電常數(shù)e2我們已經知道參數(shù)求解器可處理這個問題一一我們能夠插入圖2中，并且能夠找到必需的值。我們我們有一個低維參數(shù)空間和低維可測空間。我們能 夠使用圖形的方法

24、或者表格查詢的方法。由于我們有一個測量值和一個未知值， 因此可以歸類為尋根問題。一個策略是將 COMSOL Multiphysics正問題，例如 qi=f(e2)打包，使用MATLAB內置的尋根函數(shù)。但正如我們在第1章中所看到的， COMSOL Multiphysics的非線性求解器采用牛頓法，所以我們該如何確認 GUI 能夠為我們找到根呢？最簡單的答案是我們需要另外一個自由度-待求變量， E2,以及另外一個約束如q1=q1,m,后者為測得的通量，前者為模型計算得到的通 量。增進另外的自由度可以通過設置一個零維域來實現(xiàn)，如我們在第1章中所作的。另一個辦法是使用 ODE設置。狀態(tài)變量是全局的，能

25、夠用于模型的任意位 置?；貞浀?章中的內容，緩沖容器中的濃度 uc和vc用作狀態(tài)變量并與進口、 出口的邊界條件耦合。第三個方法是在弱模式下設置額外的自由度。零維結構對于此問題來說較為冗繁。ODE中對于狀態(tài)變量的設置較為簡單，但是變量的向后耦合不如上述兩 種方法靈活。弱結點模式的設置較為容易，但是需要對COMSOL Multiphysics平臺下弱形式的設置方法有較深入的理解。只有訓練有素的有限元數(shù)值計算師對 于弱形式和弱形式編程更為合適。我們將在這里使用弱形式，也要提請讀者注意 可以使用增加額外獨立變量和自由度的方法。表4給出了設置反演問題在 GUI中的設置。注意，由于有36個頂點，弱結 點模

26、式創(chuàng)建了 36個新的自由度。將除頂點1以外的所有頂點設置為非激活狀態(tài)， 可以得到額外的自由度?？赡苄枰獙θ鯃霰磉_式E2_test*(q1-2.07)稍加解釋，E2_test為與獨立變量E2共腕的“測試函數(shù)”，而該獨立變量在第2章中已加以定 義。在這個表達式中提供了約束要求的格式，如 q1=2.07。由于E2僅在頂點1處定義，為了使變量E2成為狀態(tài)變量，需要新定義一個積分變量(這里命名為P2)表4 ECT模型求解反問題時圖形用戶界面的修改內容。文件名 ：ectrect模型導航欄打開 ectrect_點擊 multiphysics 按鈕。COMSOL Multiphysics| PDE modes

27、 weak point form因變量名稱：E2應用模式名稱：errorPhysics 菜單 | Point settings 模式error選擇點2 36,取消選中“ Active in this domain”選項 選擇點 1。在 weak value/expressionM入 E2_test*(q1-2.07) Init選項卡。輸入E2(to)=e2Options 菜單 |Integration couplingVariablesPoint variables選擇點1輸入變量名：P2輸入表達式：E2 保持全局目標Physics 菜單 | Subdomain settings 模式g選擇域

28、2修改 I=-P2*phix-P2*phi y完成Solve菜單Solver ManagerInitial value選項卡：設止初值表遼式Solve for選項卡，只選擇模式gSolve菜單選擇 Get initial valueSolve菜單Solver Parameters選擇非線性求解器Advanced選項卡。設止解形式為 weak 點擊工具欄求解按鈕(=)Solve菜單Solver ManagerSolve for選項卡。同時選中g和error模式Post-processin球單Data Display Global輸入P2和q1查找值：0.167904,表達式：P2值：2.07,表

29、達式：q1注意在網格檢測狀態(tài)下，最終的模型有 21646個自由度。你可能想知道為什 么求解器管理器在誤差模式下調用關和開。實際上，這對我來說也是個迷。如果我試著從初始環(huán)境下同時求解兩個模式，則會看到“Sigular Matrix”的錯誤信息，表明一個奇異的DOF,使得線性求解器不能正常工作。由于我們已經一步一步 建立了該模型，我們精確的知道 DOF為E2,該變量沒有在任何方程中直接出現(xiàn)， 而是間接的收到誤差模式的約束。通過對解分段，求解器對E2通過耦合變量形成的關聯(lián)開始警惕。當我成為耦合變量的常用用戶后，這個反常的問題出現(xiàn)的越 來越頻繁，對解的分段也就成了必須要做的工作。練習7.1應該注意的是

30、在上面例子中，因為未知變量數(shù)目與測試變量相等，反演問題 即為尋根問題。GUI模型能夠用于多于一個的測量。假設q1 = 2.061215以及q2=-0.031372給出額外的自由度，反演問題成為誤差的最小化問題。典型的誤差 是平方根和誤差，例如誤差=(q1p-q1m)2+(q2p-q2m)2。下標p和m表示模型計算以及測量。盡管如此，模擬者可以選擇按比例縮放誤差以便獲得相對誤差，例如誤2Qlp -q1mOlm2。每個誤差項都是可比的。我們能夠通過在GUI中輸 q2m入上述誤差表達式來求解弱形式下的誤差最小值，例如輸入表達式E2_test*error。 我們讓牛頓求解器通過改變自由度來滿足上述不可

31、能實現(xiàn)的約束，從而得到能夠盡量滿足該約束的非收斂解，例如找到一個E2的值能夠使誤差最小。你可以發(fā)現(xiàn)可以先從線性求解器開始嘗試，然后以非收斂解作為初值重新運行求解器。1.3 在Matlab中對反演方法編程COMSOL Multiphysics GUI能夠通過使用一定的技巧來較容易的處理反演問 題（如例7.1）。盡管如此，隨著反演問題復雜度的提高，MATLAB編程的靈活性也需要隨之提高。如我們在參數(shù)研究中看到的，選擇對所求參數(shù)敏感的測量值 較為重要。盡管如此，選擇實驗條件會影響到敏感性。在我們的正問題中，我們 特意創(chuàng)建電極電壓常數(shù) phil, phi2, phi3, phi4,但是只觸發(fā)其中一個，

32、另外的 保持接地。在許多靜電裝置中，電極按照一定順序觸發(fā)，正交場泄漏可在每個電 極上測量，可得到按一定順序的關于電介質分布的信息。其中的一個傳感器或許對某些介電常數(shù)分布不敏感，但是另一個傳感器則可能較為敏感。 有效的使用這 些信息需要一個策略，但是如果有一個正問題模型，則反演問題能夠通過多個實 驗環(huán)境來表達，這里需要采用 COMSOL Multiphysics外圍的MATLAB編程。最容易的代碼編寫方法是將fem結構體保存到一個文件中，然后加載到MATLAB函數(shù)里。加載第一個正問題的例子，輸出fem結構體到MATLAB中c 使用Matlab命令來保存它：save fem下面的MATLAB函數(shù)計

33、算了正問題的解，觸發(fā)兩個電極，計算通過這些電 極的正交通量：function Q1,Q2=forward(E2)% Computes normal fluxes when electrodes 1 and 2 fire sequentially. load fem% Constants= eO' , ' 1' , ' e1 ' , ' 0.05 ',.'e2' ,E2, .e3' , ' 0.05 ' , ' e4' , ' 0.05 ',.'phi1 ; r

34、 ,.'phi2 ' , ' 0',.'phi3 ' , ' 0' , ' phi4 ' , ' 0' ;% Fire first electrodefem=multiphysics(fem); =meshextend(fem);=femnlin(fem, ' solcomp ' , ' phi ' , ' outcomp ' , ' phi ' , ' report ' , ' off '); % Co

35、mpute normal fluxQ1=postint(fem, ' nx*phix+ny*phiy ' ,' dl ' ,1)2,3,' edim ,% Constants= e0' , ' 1' , ' e1 ' , ' 0.05 ',.'e2' ,E2, .e3' , ' 0.05 ' , ' e4' , ' 0.05 ',.'phil ' , ' 0',.'phi2 ； r ,.&

36、#39;phi3 ' , ' 0' , ' phi4 ' , ' 0' ;% Fire second electrodefem=multiphysics(fem); =meshextend(fem);=femnlin(fem, ' solcomp ' , ' phi ' , ' outcomp ' , ' phi ' , ' report ' , ' off ');% Compute normal fluxQ2=postint(fem, 

37、9; nx*phix+ny*phiy ' , ' dl ' ,8,9,12,' edim ' ,1);Matlab調用較為簡單：>> Q1,Q2=forward(0.143)Q1 =2.0646Q2 =2.0035這兩個對非線性求解器的調用計算了兩個獨立模型的解，僅根據(jù)操作條件離散，下一步是創(chuàng)建誤差函數(shù)。這對于正問題較為容易，注意誤差函數(shù)必須返回一 個單值，即估計誤差，輸入向量和參數(shù)向量需要被調用：function b=errornm(v);q1,q2=forward(v);Q1=2.0646; Q2=2.0035;x=(q1-Q1)/Q1;

38、y=(q2-Q2)/Q2;v,q1,q2,sqrt(xA2+yA2)b=sqrt(xA2+yA2);如果猜測v=E2沒有給出測量的Q1=2.0646Q2=2.0035，這個函數(shù)返回了模型誤 差，例如：>>erornm(0.05)ans=0.0500 2.05012.0034 0.0070通過反復嘗試，可以使用Matlab使誤差最小化。在實際中，最好能夠使尋找 的過程自動化。Matlab中有內建的優(yōu)化函數(shù)，名稱為fminsearch,可以實現(xiàn)自動 搜尋模型最小值。這種情況下的命令為：v=0.05;fminsearch(errornm,v)符號處理了函數(shù)名稱。第二項表達式為初始條件。f

39、minsearch()提供了多變量標量函數(shù)最小化的簡便算法。它實現(xiàn)了 Nelder-Mead單純形搜索算法，能夠 通過改變輸入參數(shù)v來尋找f(v)的最小值。這個方法與其他的最小值計算方法相 比效率并不算高，尤其在涉及微分方程運算時體現(xiàn)的更明顯，但是從另一方面講， 該方法不需要采用耗費更多系統(tǒng)資源來計算函數(shù)梯度，因此對于平滑函數(shù)的魯棒 性更好。如果希望在較小的計算量情況下求方程的極小值，那么 Nelder-Mead算 法通常較為適用。匕=0.143時，使用fminsearch函數(shù)共耗費28步來得到最小誤差。每步計算 需要求解26000個自由度，這是一個計算量較大的方法。我們已經看過另外兩種不同的

40、求解方法一一 GUI為基礎以及Matlab為基礎 來計算反演問題。我知道另一個COMSOL Multiphysics中可能較為有用的反演問 題求解方法。COMSOL公司的Niklas Rom在Knowledge Base上發(fā)表了一種根據(jù) 已知濃度分布計算擴散率的方法。 該方法涉及了穩(wěn)態(tài)格式的誤差函數(shù)，要求得擴 散率使得函數(shù)誤差達到最小。這個方法對于參數(shù)關系可以顯示表達出來得情況較 為有用。以上兩個方法中得梯度均為黑箱，并沒有顯式給出。如果域中對組分得位置或者尺寸有更高的分辨率要求，則需要給出更高質量的邊界信息，可能會使包含的數(shù)據(jù)模糊化。圖像重構對于電容層析成像技術來說 是一個復雜的問題。Dya

41、kowski等5對于該類應用給出了較好的綜述，特別是基于MATLAB的軟件包IDORS是先進反演問題技術中最好的源代碼。練習7.2寫出求解練習7.1反演問題的 matlab腳本并進行求解，比較 GUI求解和 MATLAB優(yōu)化求解分別耗費的時間，哪個方法耗費的CPU時間較少？2 .投影耦合變量與線積分投影耦合變量在特定坐標系下沿二維域進行線積分，例如路徑積分。這個概 念很有用，例如，在描述量子電動力學時會涉及14。作為線積分最簡單的形式， 路徑可以選為某個坐標軸(簡單的網格線)，可以降低域或者變量的依賴性。.y2 xI x f x, y dy 在域 口2上(6)y xI(x)為耦合變量，必須在域

42、D1中定義，該域的維數(shù)應小于D2=(x1,x2) Xy1(x),y2(x)。 更復雜的映射能夠通過網格變換來實現(xiàn)，使用空間坐標(獨立變量x,y,z.)或者局部網格參數(shù)，例如，s, si或者s2,然后這些參數(shù)用于生成新的源網格，以便獲 取線段或者映射積分的計算區(qū)域。例如：I x f s , s ds(7)C x其中，曲線C由x進行參數(shù)化。因此，目標域中對應每一點 x均有這樣的一 條曲線?？傮w來說，映射耦合變量應比源域低一階，因此必須在一個新的域中定 義，可能需要顯式創(chuàng)建以便在其目標域中得到耦合變量。作為一種內在特征，映射耦合要求兩個目標域(目標可能是源域的邊界)，因此從初始狀態(tài)開始就必須 具備兩

43、個域。耦合變量I(x)比一個線積分包含更多的信息。因此，如果你對線積分的特定 值感興趣，那么你需要點擊目標域中的點，信息窗口中將顯示在該點時耦合變量 的內插值。換一種方法，你可以輸出FEM結構，使用postinterp函數(shù)來提供一個 數(shù)值解。2.1 例：激光定位和污染云分散問題激光雷達與其他光學設備的工作原理一樣，例如光譜測定法和頻譜學.某一特 定波長的光的接收信號總會由于化學組分的吸收作用而導致光強度降低。對于稀釋的化學物質，接收信號與沿光路方向的物質積分濃度成正比，例如：proj c x s ,y s ds(8)C其中，C為曲線(x(s),y(s), s為沿弧長方向的坐標。c(x,y)為物

44、質濃度場。假設域 是準二維的，且雷達陣列沿x軸方向排布,x軸即映射到0,11的域的下邊界。然后，雷達陣列接收到離散投影耦合變量proji：1proj1 x °c x, y dy(9)這里將(8)中曲線C考慮為垂直線。這是COMSOL Multiphysics投影耦合變量在 2D域上的標準用法，投影耦合變量僅定義為 1D獨立變量的函數(shù)，局部網格變換的默認選擇對于源域為x x,y y ,對于目標域為x x ,在單位正方形 上得到得到方程(9)。如果使用局部坐標，非矩形域的選擇更為重要：2D上的(si,s2) 作為曲線邊界，1D中的s用于非線性曲線?，F(xiàn)在回頭看例子，致密的大氣污染物氣體瞬時

45、排放模型是一個2D的高斯分布。Zimmerman和Chatwin17分析了密集氣體排放的風洞數(shù)據(jù)，表明流體的瞬 時結構是高度間歇性的。然后，當云開始稀釋時，需得到2-D高斯分布的統(tǒng)計平 均值或者窗口時間平均值。因此，我們假設初始濃度分布：22x xoy y。Co x, y exp -2 -2(10)1 x1 y假設該分布所在的速度場是均勻的 向量為：U=(U0,V0)o雷達初始測得的水平和垂直投影proj 1 xC| exp2x。lxproj 2 yc2 exp2yy。(11)(12)(13)令M為總濃度:M c x, y d然后這些映射耦合變量可被認為是歸一化的條件密度函數(shù),對于 proj1

46、:1 proj1 xm* 1xdxx,10 M1 2 proj1 x mx 2x dx0 M對于c=co(x,y),初始高斯分布能夠顯示x軸上的靜矩，mx,1=x。.二次矩可用于計算標準差:sxm2 lx,2mx,1:2(14)以上結論是在y坐標軸上的x投影，所以投影能夠作為高斯分布形狀的云來獲得 位置和大小。對于非高斯云，他們至少能夠提供集中性和速度程度方面的信息?，F(xiàn)在我們在 COMSOL Multiphysics中演示這個例子。打開COMSOL Multiphysics,在模型導航器中進行如表5操作。表6和表7演示了后處理過程， 包括計算靜矩和標準差，t=0.06時，根據(jù)(13)式的計算在

47、表6中。后面的兩個積分值為sx=0.2738, sy=0.2765，幾乎是等速擴散，但是這是在擴 散的最后狀態(tài)出現(xiàn)的近似。時間為0時，表7中得出了類似的結果。后面的兩個 積分值為sx=0.0707,sy=0.0872，與lx,ly的設置保持一致。圖3給出了擴散非常迅速的初始條件。盡管模擬中取Pe=1,在這樣的網格分辨率下，數(shù)值擴散率較大。在一個細化的網格上，擴散速度可能較低。圖 4 給出了映射耦合變量示范的由人造 雷達”捕捉的近高斯分布，出現(xiàn)了由(13)和(14) 確定的中間峰值區(qū)間。表5映射耦合變量模仿雷達模型模型導航欄選擇二維空間選擇 Chemical Engineering| Mass

48、balance Convection and Diffusion application mode (chcd)接受默認設定Options 菜單 | Constants常數(shù)名稱：xo表達式：0.4常數(shù)名稱：yo表達式：0.6常數(shù)名稱：lx表達式：0.1常數(shù)名稱：ly表達式：0.12常數(shù)名稱：u0表達式：1常數(shù)名稱：V0表達式：0常數(shù)名稱：Pe 表式：sqrt(u0A2+v0A2)Options 菜單 |Axes/Grid settings設定柵格為(-0.1,1.1)x(-0.1,1.1),柵格間距 0.1, 0.1Draw菜單|Specify Objects選擇正方形，邊長為1,左下角點在原

49、點Physics 菜單 |Subdomain settings選擇子域1D=1/Pe； u=u。； v=v0Init 選項卡：c(t0)=exp(-(x-x0)A2)/lxA2-(y-y0)A2/lyA2)Physics 菜單 |Periodic Conditongs|Periodic boundaryConditions選擇邊界1。Source選項卡：輸入表式c 點擊約束名欄Destination選項卡：選擇邊界 4,輸入表達式c 源頂點：1, 2。目標頂點：3, 4。Source選項卡：選擇邊界2,輸入表達式c點擊約束名欄Destination選項卡，選擇邊界3,輸入表達式c 源頂點：1,

50、 3。目標頂點：2, 4。完成Options 菜單 | Projection CouplingVariables| Subdomain variablesSource選項卡。選擇子域1。輸入名稱proj1。表遼式：c 選擇 general transformationx-x, yy）Destination選項卡。選擇邊界2Evaluation point （xx）Source選項卡。選擇子域1。輸入名稱proj2。表達式：c 選擇 general transformation一x, yy）Destination選項卡。選擇邊界1Evaluation point （xx）Solve菜單|Solv

51、er Parameters選擇時間依賴求解器General選項卡：設定輸出時間為 0:0.001:0.06 點擊工具欄上求解按鈕（=）Post-processing Subdomain integration通過積分c計算總質量（任何時刻）Post-processingBoundary integration分別對初始和最終情況，計算proj2*x/0.037699和 proj2*y/0.037699 （一階動量）和 proj2*xA2/0.037699 和 proj2*yA/0.037699 （二階累計量）的邊界積分表6 t=0.06s時的矩和累積量計算后處理模式t=0.06子域積分：域1

52、c (任意時間)11=0.037702邊界積分：邊界 2 proji*x/0.037702 12=0.49325邊界積分：邊界 1 proj2*y/0.037702 13=0.51522邊界積分：邊界 2 proji*xA2/0.037702 14=0.31825邊界積分：邊界 1 proj2*yA2/0.037702 15=0.34188表7 t = 0時的矩和累積量計算后處理模式t=0邊界積分：邊界 2 proj1*x/0.037702 12=0.39999邊界積分：邊界 1 proj2*y/0.037702 13=0.60012邊界積分：邊界 2 proj1*xA2/0.037702 1

53、4=0.165邊界積分：邊界 1 proj2*yA2/0.037702 15=0.367261I! rGH-H1-.!:;.!: 一 "!"二"E F 回 皿5-值|0.|?|'.».3|?人|0.5 32|»|.3口55 U-r二二二PT 二二二 一，T1«圖3時間t =0.001（左）與t =0.06（右）時根據(jù)對流擴散模型，Pe=1，均勻水 平流雙重周期邊界條件下濃度場的等密度線。(b) boundary 1型5 L 1kl I I I I a 0.1 QJ DJ DI 4 QEi DE D7 口日 0 9-y-coor

54、dinate圖4 t =0.06時，圖3所示模型水平邊界（左）和垂直邊界（右）的線性積分映射。練習7.3:假（數(shù)值）擴散在細化的網格上重復雷達的例子。 高斯云是否分散的慢了一些？你怎樣使用 這個計算例子來量化人為造成的數(shù)值擴散率？3.伸展耦合變量擴展耦合變量以其最常見的用途來命名，它將信息從一個n維域映射到n+1域。然而伸展只是它的一個可能的應用，還可以用到插值、投影、映射，具體場 合視信息而定。另外兩個耦合變量類型-標量、 投影，在整個源域或子域中進行 積分，因此能夠用于積分方程。伸展耦合變量映射了從一個域到另一個域的細節(jié) 及分布信息。目標位置由局部網格變換得到。 所以，伸展耦合變量對于多幾

55、何域 耦合的模型中間步驟較為有用。他們需要在不同的幾何域上定義。在 COMSOL Multiphysics討論會上，伸展耦合變量的一個常用例子是出于美觀的原因而選擇 的。通常，給出物理結構的對稱性，模型能夠在一部分區(qū)域上求解，甚至簡化到 低維求解，但需要采用真實的物理結構來進行后處理和可視化。8坐標和r-z坐標方向的 伸展可滿足可視化需求。假設隔板存在六邊形對稱性，則可在0,冗/3 區(qū)間上進行求解，以r和z為邊界。如果需要在整個流道中實現(xiàn)可視化，可以將 隔板伸展為其余五個隔板。所以，伸展耦合變量也可以將用于后處理的信息映射 到其他區(qū)域。3.1 積分方程積分方程的特點在于積分中包含未知函數(shù)。在微

56、分方程中，線性方程組研究 的最為充分，因此最為常見。方程組的分類較為直接。如果積分上、下限確定，方程為 Fredhom類型。如果上、下限中的一個為變 量，則為Volterra類型。如果未知函數(shù)僅僅在積分符號下出現(xiàn)，則稱為“第一類”，如果積分符號里面外面都出現(xiàn)了未知函數(shù)，則稱之為“第二類”。這里是四種混合類型的積分方程：第一類Fredhom積分方程:f第二類Fredhom積分方程：g x第一類Volterra積分方程：f第二類Volterra積分方程：g xbx a K x,t g t dtbf x a K x, t g t dtxx a K x,t g t dtxf x K x,t g t d

57、ta(15)(16)(17)(18)在上述所有類型的積分方程中，g為未知函數(shù)，K(x,t)稱為核心，并假設f(x)為已 知的。當f(x)=0時，方程為均相。讀者很可能問為什么我們要為積分方程如此費心。這是因為本章的主題-積 分方程是非局部的。一些物理現(xiàn)象的特性從本質上講即非局部的。 所以，對他們 的描述引出了積分或者積分微分方程。例如， Shaqfeh得到了混合材料的傳遞 特性，是對有效特性的非局部描述。許多方程從本質上講是 橢圓”的-邊界數(shù)據(jù) 在域中各處均有效，如穩(wěn)態(tài)傳熱或傳質，這使得一個點的解需取決于域中各點的解。這樣的非局部方程組能夠方便的以格林函數(shù)的形式來描述，從而引出對非均相系統(tǒng)的積分方程描述。最后，一些過程可以由相空間( Fourier空間，Laplac