基于氧擴(kuò)散問(wèn)題的半光滑牛頓算法
2022-08-08 05:45:44曹夢(mèng)霖宇振盛
曹夢(mèng)霖, 宇振盛
(上海理工大學(xué) 理學(xué)院, 上海 200093)
0 引言
氧擴(kuò)散問(wèn)題主要出現(xiàn)在生物醫(yī)學(xué)和流體力學(xué)等領(lǐng)域,例如細(xì)胞內(nèi)的氧氣吸收、多孔介質(zhì)的原位燃燒[1-2]和晶體退火過(guò)程中的氧擴(kuò)散等。實(shí)質(zhì)上,解決氧氣擴(kuò)散問(wèn)題等價(jià)于求解氧擴(kuò)散方程,該方程是移動(dòng)邊界問(wèn)題的一類重要分支。
氧氣擴(kuò)散方程是典型的非線性問(wèn)題,其解析解很難獲得，因此,在過(guò)往的研究中一般采用數(shù)值方法進(jìn)行求解。Crank等[3]將積分法首次應(yīng)用于求解氧擴(kuò)散方程,對(duì)于移動(dòng)邊界位置的網(wǎng)格采用拉格朗日型公式,并利用泰勒級(jí)數(shù)展開(kāi)、歐拉方程和固定網(wǎng)格差分法計(jì)算出數(shù)值解,但該方法在小范圍內(nèi)是無(wú)效的。Guptar等[4]采用約束積分法,利用可變的時(shí)間步長(zhǎng)確定離散區(qū)間,避免固定時(shí)間步長(zhǎng)造成的計(jì)算量過(guò)大問(wèn)題。Chapiro等[5]基于傳統(tǒng)的隱式中心空間有限差分方案的數(shù)值方法與非線性互補(bǔ)算法相結(jié)合,與積分法相比,該方法允許更大的時(shí)間步長(zhǎng),提升了計(jì)算效率。Mitchell[6]消除邊界條件和初始條件的不一致性,并利用最優(yōu)積分法解決之前積分法中無(wú)法從t=0開(kāi)始求解的問(wèn)題。
本文仍采用隱式有限差分格式和非線性互補(bǔ)算法求解氧擴(kuò)散方程，但與文獻(xiàn)[5]不同的是,引入非線性互補(bǔ)函數(shù)(NCP-function),將非線性互補(bǔ)問(wèn)題轉(zhuǎn)化為基于NCP函數(shù)的半光滑方程組,從而使問(wèn)題轉(zhuǎn)化為無(wú)約束問(wèn)題,進(jìn)而使用廣義牛頓法進(jìn)行求解,相較于文獻(xiàn)[5]的內(nèi)點(diǎn)法,在線搜索時(shí)無(wú)需進(jìn)行可行性判別,減少了計(jì)算步驟,數(shù)值結(jié)果表明該算法有效。
1 一維氧擴(kuò)散模型
1.1 初始模型
氧擴(kuò)散模型如圖1所示。氧氣擴(kuò)散問(wèn)題分為2個(gè)階段:第一階段,介質(zhì)沒(méi)有封閉,氧氣可以自由進(jìn)出,直到達(dá)到穩(wěn)定狀態(tài)。在穩(wěn)定狀態(tài)下,介質(zhì)中再無(wú)氧氣進(jìn)出。第二階段,介質(zhì)的表面是封閉的,氧氣無(wú)法進(jìn)入介質(zhì)或擴(kuò)散到外界。介質(zhì)中氧氣吸收的持續(xù)性導(dǎo)致介質(zhì)開(kāi)始消耗內(nèi)在氧氣。此時(shí),穩(wěn)定狀態(tài)下滲透深度的邊界開(kāi)始向密封面后退。本文主要研究氧擴(kuò)散問(wèn)題的第二階段。
(a)第一階段 (b)第二階段
第二階段實(shí)質(zhì)上是一維單向移動(dòng)邊界氧氣擴(kuò)散問(wèn)題,其數(shù)學(xué)模型[3]為
式中：C(X,T)為T時(shí)刻距離介質(zhì)外表面X處自由擴(kuò)散的氧氣濃度；D為恒定的擴(kuò)散系數(shù)；r(C)為單位體積介質(zhì)內(nèi)的氧消耗速率。
表面初始位置和移動(dòng)邊界位置分別為X=0及X0(t)。在初始位置被密封后(X=0),模型為
邊界條件:
初始條件:
1.2 無(wú)量綱形式模型
(1)
(2)
(3)
(4)
式中s(t)表示x對(duì)應(yīng)于X0(T)的值,即在無(wú)量綱形式下移動(dòng)邊界的位置。當(dāng)r(c)=1時(shí),文獻(xiàn)[5]應(yīng)用FDA-NCP算法對(duì)線性情況進(jìn)行了研究。本文將考慮文獻(xiàn)[7]r(c)=1+c(1-c)的非線性情況。
2 求解氧擴(kuò)散方程算法設(shè)計(jì)
非線性互補(bǔ)問(wèn)題(NCP)的一般形式為,給定向量函數(shù)F(c):Rn→Rn求解c∈Rn,滿足
c≥0,F(c)≥0,cTF(c)=0,
依據(jù)文獻(xiàn)[8],將式(1)至式(4)轉(zhuǎn)化為變分形式,
(5)
且滿足
(6)
當(dāng)0
另證,當(dāng)x∈[0,1],式(5)成立，因此,等式(6)是有效的。
對(duì)于求解非線性互補(bǔ)問(wèn)題的一種重要途徑是將其轉(zhuǎn)化為等價(jià)的無(wú)約束問(wèn)題,主要手段是利用非線性互補(bǔ)函數(shù)函數(shù)(NCP-函數(shù)),將其轉(zhuǎn)化為等價(jià)的非線性方程組。本文使用經(jīng)典的Fischer-Burmeister函數(shù)(簡(jiǎn)稱FB互補(bǔ)函數(shù))[9],下面給出NCP函數(shù)的定義和FB函數(shù)的具體形式。
定義1 設(shè)函數(shù)φ:R2→R,如果具有
φ(a,b)=0?a≥0,b≥0,ab=0,
(7)
則稱φ為一個(gè)NCP函數(shù)。
由式(7)可得,求解非線性互補(bǔ)問(wèn)題等價(jià)于求解如下方程組Φ(c),其具體形式為
(8)
FB再生方程所得自然勢(shì)函數(shù),定義為
(9)
φFB(c,F(c))是除原點(diǎn)外處處連續(xù)的二維凸函數(shù),由凸優(yōu)化理論可得,φFB和ΦFB滿足局部Lipschitz連續(xù),因此ΦFB在Rn中的任何點(diǎn)上具有Clarke意義下的廣義Jacobian矩陣?Φ。
性質(zhì)1 函數(shù)Φ是半光滑的；函數(shù)Ψ(c)是連續(xù)可微的,且梯度可以表示為?Φ(c)Φ(c)。
性質(zhì)2[10]對(duì)?c∈Rn,有
?Φ(c)T?(A(c)-I)+F(c)(B(c)-I),
(10)
式中I是n×n階的單位矩陣,A(c)和B(c)是對(duì)角矩陣,定義為
A(c)=diag(Aii(c)),B(c)=diag(Bii(c)),
式中(ξi,ρi)滿足‖(ξi,ρi)‖≤1。
本文,采用文獻(xiàn)[10]的半光滑牛頓算法，該算法實(shí)質(zhì)是使用廣義牛頓算法求解FB再生方程組,下降方向選取牛頓方向,如果下降效果不好,則選取負(fù)梯度方向代替。
基于FB函數(shù)的算法1,步驟如下:
Step2:如果滿足‖Ψ(ck)‖≤ε,則停止,否則轉(zhuǎn)Step 3；
Step3:計(jì)算搜索方向,Hkd=-Φ(ck),其中Hk∈?BΦ(ck)。如果不可解,或Ψ(ck)Tdk≤-ρ‖dk‖p不成立,則選取dk=-Ψ(ck)；
Step4:更新c,令ck+1=ck+λkdk,其中λk=2-ik,而ik=0,1,2,…,k中使Ψ(ck+λkdk)≤Ψ(ck)+βλkΨ(ck)Tdk成立的最小者；
Step5:k=k+1,返回Step 2。
下面給出,選取Clarke意義下的廣義Jacobian矩陣H的算法2,步驟如下:
Step1:令Q={i|ci=Fi(c)=0}；
Step2:選取z∈Rn,對(duì)?i∈Q,zi≠0；
Step3:構(gòu)造矩陣H的第i行為
首先,算法1產(chǎn)生的搜索方向dk一定是下降方向。由于函數(shù)Ψ連續(xù)可微,滿足Ψ(ck)=VTΦ(ck),?V∈?BΦ(ck),并且由算法2產(chǎn)生的具有Clarke意義下的廣義Jacobian矩陣Hk滿足Ψ(ck)Tdk=Φ(ck)THkdk=-‖Φ(ck)‖2。
下面給出算法1的收斂性:
性質(zhì)3[11]如果H在點(diǎn)c∈Rn是BD-正則的,則①必存在c的一個(gè)鄰域N和一個(gè)常數(shù)l>0,使對(duì)?y∈N和V∈?BH(y)都有V非奇異且‖V-1‖≤l;②進(jìn)一步,如果H(c)=0且H在點(diǎn)c是半光滑的，則存在一個(gè)鄰域N′和一個(gè)常數(shù)θ>0,使對(duì)?y∈N′都有‖H(y)‖≥θ‖y-c‖。
定理1[10]對(duì)于算法1產(chǎn)生的無(wú)窮點(diǎn)列{ck}滿足:①序列{ck}任意的一個(gè)聚點(diǎn)都是函數(shù)Ψ的穩(wěn)定點(diǎn)；②設(shè)c*是一個(gè)極限點(diǎn),且c*是NCP的一個(gè)孤立解,則ck→c*。
定理2[10]假定算法1產(chǎn)生一個(gè)無(wú)窮點(diǎn)列{ck},若c*是它的一個(gè)極限點(diǎn)且是Φ(c)=0的一個(gè)BD-正則解,則①{ck}超線性收斂到c*；②若F局部Lipschitz連續(xù),則收斂速率是二階的。
3 有限差分格式
氧擴(kuò)散方程是一類典型拋物線類型的偏微分方程,想要得到其解析解十分困難。為了研究此類偏微分方程,需要將其進(jìn)行離散化,從而得到代數(shù)方程組,最終得到此類問(wèn)題的近似解。
將使用Crank-Nicolson有限差分格式[7]進(jìn)行離散,并使用算法1求解每個(gè)時(shí)間步內(nèi)的互補(bǔ)問(wèn)題。由于在離散格式與算法選取上的獨(dú)立性,因此解決方法的選擇更加靈活。由于柯朗-弗里德里希斯-列維條件(Courant-Friedrichs-Lewy, CFL)穩(wěn)定性的限制,顯式方案時(shí)間步長(zhǎng)較小,因此求解時(shí)間長(zhǎng)。然而,隱式方案可以允許更大的時(shí)間步長(zhǎng),并且可能是無(wú)條件穩(wěn)定的,但它們求解較為復(fù)雜。
給出具有M+1個(gè)點(diǎn)關(guān)于變量x的均勻網(wǎng)格,其中x0和xM是計(jì)算區(qū)間的邊界點(diǎn)?？臻g網(wǎng)格間距Δx=xm+1-xm=1/M。類似地,t為時(shí)間,定義時(shí)間指數(shù)為n,Δt為時(shí)間網(wǎng)格間距。
離散形式為
(11)
其中
和
下面通過(guò)結(jié)合算法1和有限差分格式,求解氧擴(kuò)散方程算法3,步驟如下:
Step2:應(yīng)用Crank-Nicolson有限差分格式,獲得離散算子FΔ(c)；
Step3:調(diào)用算法1以及離散算子FΔ(c)去求解式(8)；
Step4:在下一個(gè)時(shí)間步驟中,將Step 3中的解賦值給變量c；
返回Step 2中的算法,令n=n+1,當(dāng)n=N時(shí)算法停止。
4 數(shù)值結(jié)果
半光滑牛頓算法對(duì)應(yīng)于所有時(shí)刻的氧氣濃度數(shù)值模擬如圖2所示。給出文獻(xiàn)[5]中的FDA-NCP算法和算法1分別結(jié)合有限差分法對(duì)應(yīng)于t為0、0.06、0.10、0.15、0.18的數(shù)值模擬,如圖3所示。
圖2 半光滑牛頓算法對(duì)應(yīng)于所有時(shí)刻的氧氣濃度數(shù)值模擬
5 結(jié)語(yǔ)
結(jié)合一維氧擴(kuò)散方程模型,給出利用Fischer-Burmeister函數(shù)求解非線性互補(bǔ)問(wèn)題的半光滑牛頓算法,并在離散方法上選擇了Crank-Nicolson差分格式。通過(guò)數(shù)值模擬,與FDA-NCP算法對(duì)比,驗(yàn)證算法的有效性。對(duì)于其他移動(dòng)邊界問(wèn)題或者采用其他NCP函數(shù),如原位燃燒模型和分段NCP函數(shù),值得進(jìn)一步研究。
