摘"要：為實(shí)現(xiàn)高階高度非線性強(qiáng)耦合偏微分方程組（Partial Differential Equations， PDEs）的數(shù)值求解，本文以滲蝕強(qiáng)耦合PDEs為典型案例，剖析了PDEs的高階高度非線性，結(jié)合空間映射，基于弱形式建模與分離式算法，實(shí)現(xiàn)了滲蝕強(qiáng)耦合PDEs的數(shù)值求解，并驗(yàn)證了求解方法的可行性與可靠性。研究表明：強(qiáng)耦合是導(dǎo)致PDEs非線性特性的充分條件；非弱形式建模難以妥善解決高階高度非線性強(qiáng)耦合PDEs數(shù)值收斂性問題；分離式求解算法對(duì)于高階高度非線性強(qiáng)耦合PDEs的初始條件更具包容性。

關(guān)鍵詞：多場強(qiáng)耦合PDEs；高階高度非線性；弱形式；求解方法

中圖分類號(hào)：O175.29"""""文獻(xiàn)標(biāo)識(shí)碼：A """"""文章編號(hào)：20959699（2024）03000106

物理問題的突破與數(shù)學(xué)方法的革新密切相關(guān)。目前，系列關(guān)鍵物理現(xiàn)象已被相關(guān)數(shù)理方程所描述。其中，研究者開發(fā)了眾多極具價(jià)值的偏微分方程（Partial Differential Equation， PDE）以預(yù)測相應(yīng)關(guān)鍵物理過程。早期，受限于研究手段，從業(yè)人員依據(jù)需求進(jìn)行必要簡化，通過單個(gè)PDE詳細(xì)地剖析了相關(guān)物理規(guī)律，構(gòu)建了單物理場認(rèn)知體系。顯然，自然界本身，客觀上存在著極為復(fù)雜的事實(shí)——多物理場相互耦合。隨后，研究熱點(diǎn)轉(zhuǎn)向多物理場耦合。期間，研究人員在多個(gè)方面取得了進(jìn)步，包括：聯(lián)立PDEs從兩場耦合到多場耦合、從弱耦合到雙向強(qiáng)耦合以及從雙向強(qiáng)耦合到多向強(qiáng)耦合等[1]。然而，隨著研究的進(jìn)一步深入，多向強(qiáng)耦合PDEs逐漸高階高度非線性化，使得常規(guī)求解策略難以妥善解決非線性問題而無法滿足收斂標(biāo)準(zhǔn)[2]。

為實(shí)現(xiàn)高階高度非線性強(qiáng)耦合PDEs的數(shù)值求解，本文將以突水突泥潛伏期滲蝕四場強(qiáng)耦合PDEs為典型案例，重點(diǎn)論述PDEs的高階高度非線性強(qiáng)耦合特性，分析空間映射在滲蝕強(qiáng)耦合PDEs求解過程中的必要性，通過對(duì)比討論，闡述高階高度非線性強(qiáng)耦合PDEs求解條件的苛刻性，論證最優(yōu)求解策略的可行性，最后，結(jié)合試驗(yàn)數(shù)據(jù)，驗(yàn)證求解方法的可靠性，僅供類似問題參考。

1"強(qiáng)耦合PDEs及其高階高度非線性特性

1.1"強(qiáng)耦合PDEs

表 1所示的四場多參數(shù)強(qiáng)耦合PDEs具有典型的高階高度非線性特性，其考慮了結(jié)構(gòu)變形與空間映射，通過孔隙率場（φ）、壓力場（PF）、顆粒濃度場（cv）以及位移場（u）描述了隧道突水突泥地質(zhì)災(zāi)害潛伏期的滲蝕致災(zāi)過程，重釋了滲蝕與孔隙率的時(shí)空關(guān)聯(lián)性，揭示了滲蝕作用與力學(xué)行為的因果關(guān)系，可預(yù)測隧道突水突泥地質(zhì)災(zāi)害的最先失穩(wěn)區(qū)域與臨界撤離時(shí)刻。對(duì)于滲蝕強(qiáng)耦合PDEs的數(shù)值求解，本文將以LT隧道ZK72＋177.8里程突水突泥潛伏期的滲蝕致災(zāi)過程為案例背景來進(jìn)行。為便于計(jì)算，將致災(zāi)結(jié)構(gòu)的幾何形狀簡化為塌體和洞渣，并對(duì)其進(jìn)行了受力分析，見圖1。根據(jù)研究重點(diǎn)，此處將不再贅述PDEs的推導(dǎo)過程及工程背景，相關(guān)內(nèi)容見文獻(xiàn)[2]。

1.2"高階高度非線性強(qiáng)耦合特性

由PDE所描述的數(shù)理模型存在某個(gè)或者某些輸入條件依賴于模型的輸出的情況時(shí)，則該P(yáng)DE具有非線性。這些條件或是PDE的系數(shù)（通常為物性參數(shù)），抑或是初始條件和邊界條件，一般以代數(shù)、微分或積分方程的形式存在。此外，多物理場耦合問題的“強(qiáng)”與“弱”，描述的是PDEs之間的相互依賴程度。強(qiáng)耦合PDEs中，任意一個(gè)PDE的輸出依賴于剩余的所有PDEs的輸出。對(duì)于滲蝕強(qiáng)耦合PDEs，在孔隙率場的源項(xiàng)ξ×（2κφ/φ‖SymbolQC@PF+ρFgZ‖－τc）+SymbolQC@u·中，2κφ/φ依賴于孔隙率場的解，即孔隙率場具有直觀的源項(xiàng)非線性，與此同時(shí)，SymbolQC@PF+ρFgZ則建立了孔隙率場與壓力場、濃度場（ρF=（1－cv）ρw+cvρS）的耦合，而SymbolQC@u·則表示位移場的輸出參與了孔隙率場的求解；壓力場的系數(shù)kφ，cv=κφ/μ（cv）具有來自孔隙率場和濃度場的貢獻(xiàn)，同時(shí)，φ與cv依賴于PF，即kφ，cv取決于PF，故此處耦合與非線性共存，此外，壓力場源項(xiàng)通過kφ，cvρFg再次將壓力場與孔隙率場和濃度場進(jìn)行了耦合，同樣地，－SymbolQC@u·表示壓力場的輸出受制于位移場的輸出。從濃度場系數(shù)－kφ，cv（SymbolQC@PF－ρFg）中可以看出，耦合（與壓力場、孔隙率場）與非線性共存，另外，在濃度場PDE中，φ·su=φ·－ε·v既是系數(shù)也是源項(xiàng)，是濃度場與孔隙率場和位移場強(qiáng)耦合的雙重體現(xiàn)；無獨(dú)有偶，位移場等號(hào)左側(cè)耦合了孔隙率場，而φ的決定因素包括位移，故此處同為耦合與非線性共存，位移場與剩余兩場的耦合則體現(xiàn)在源項(xiàng)SymbolQC@PF－ρ-g中，其中，ρ-=（1－φ+φcv）ρS+ρw。不難看出，對(duì)于強(qiáng)耦合PDEs而言，只要其中一個(gè)PDE的某個(gè)局部具有非線性，則整個(gè)PDEs具有非線性。需要補(bǔ)充的是，PDEs具有非線性并不能說明其為強(qiáng)耦合，弱耦合PDEs同樣可具有非線性。顯然，強(qiáng)耦合是導(dǎo)致PDEs非線性特性的充分條件。

2"數(shù)值解實(shí)現(xiàn)方法

2.1"空間映射

滲蝕強(qiáng)耦合PDEs屬于典型的流固耦合模型，涉及絕對(duì)空間、材料空間、幾何空間和網(wǎng)格空間。起初，四者相互重合，由于流體流動(dòng)產(chǎn)生的力與物理場的誘導(dǎo)，引起結(jié)構(gòu)在材料坐標(biāo)系中產(chǎn)生了變形，導(dǎo)致其他三者與絕對(duì)空間產(chǎn)生分離。即變形之后，原本存在材料的區(qū)域，現(xiàn)已無材料但有網(wǎng)格，反之則反。由于自由度被定義在網(wǎng)格節(jié)點(diǎn)之上，因此，必須確保網(wǎng)格在每個(gè)時(shí)步緊隨材料邊界。

再者，PF屬于流場問題，研究流場中流體的物理變化，緊跟流體質(zhì)點(diǎn)并不現(xiàn)實(shí)，須基于絕對(duì)空間。而φ、cv和u則屬于材料空間，需緊跟物質(zhì)質(zhì)點(diǎn)研究，而同時(shí)包含絕對(duì)空間和材料空間的模型最終自然需在絕對(duì)空間中進(jìn)行運(yùn)算。

由此可見，必須更新網(wǎng)格節(jié)點(diǎn)的坐標(biāo)，通過ALE理論引導(dǎo)網(wǎng)格在四種空間中相互映射以追隨材料邊界。這并非表示在每個(gè)時(shí)步都要重新剖分網(wǎng)格，而是讓網(wǎng)格節(jié)點(diǎn)同樣具有自由度——空間網(wǎng)格位移，通過網(wǎng)格平滑實(shí)時(shí)匹配移動(dòng)的材料邊界。

2.2"基于PDEs弱形式求解

對(duì)于滲蝕強(qiáng)耦合PDEs，在目前的商業(yè)軟件中并無對(duì)應(yīng)或類似的內(nèi)置模塊，須基于方程建模，才能實(shí)現(xiàn)數(shù)值求解。Comsol Mutiphysics為此類問題的研究提供了契機(jī)，其內(nèi)嵌系數(shù)型、一般型以及弱形式型PDE三種接口，三者靈活性逐一增強(qiáng)，與此同時(shí)，使用難度也逐步增大。特別是弱形式PDE接口，其開放性強(qiáng)，為復(fù)雜PDEs的求解提供了更多可能。滲蝕四場多參數(shù)強(qiáng)耦合PDEs屬于復(fù)雜高階非線性偏微分方程系統(tǒng)。系數(shù)型與一般型PDE接口因無法妥善解決非線性問題所以其并不適用。因此，本文將通過弱化來實(shí)現(xiàn)高階高度非線性PDEs的數(shù)值求解。

弱化過程靈活多變，雖有基本的規(guī)則，但無統(tǒng)一的方法。從其結(jié)果來看（表 2），弱形式降低了解變量的導(dǎo)數(shù)階，從而降低了對(duì)解變量連續(xù)性的要求，但針對(duì)實(shí)際問題而言，弱形式相較于原PDE更加逼近解析解，而這個(gè)功勞要?dú)w于試函數(shù)。事實(shí)上，弱化的第一步就是用試函數(shù)與整個(gè)PDE做內(nèi)積，即弱化就是正交過程，而正交則意味著希爾伯特空間距離最短，誤差最小。毋庸置疑，弱化后的方程不僅與原PDE等價(jià)，而且降低了計(jì)算成本。更重要的是，弱形式PDE更適合求解非線性多場耦合問題。

需要指出的是，對(duì)于強(qiáng)耦合模型的邊界條件，需要厘清其屬性。若某類邊界上自由度的一階變分為0，則此類邊界屬于Dirichlet邊界，用非弱形式直接設(shè)定這類邊界即可。對(duì)于Neumann邊界，可能有來自其他場的自由度，這些自由度并非事先已知，而是來自其他場的解。例如，從壓力場方程的弱形式中可以看出：壓力場邊界弱形式中有來自位移場的速度解。因此，需要將已知條件代入邊界弱表達(dá)式，派生出各自邊界的具體弱表達(dá)式。LT隧道ZK72＋177.8里程突水突泥潛伏期滲蝕強(qiáng)耦合PDEs的邊界條件見表3。

2.3"數(shù)值精度分析與網(wǎng)格獨(dú)立性驗(yàn)證

在正式計(jì)算前，有必要進(jìn)行數(shù)值精度與網(wǎng)格獨(dú)立性分析。一般情況下，模型的網(wǎng)格剖分策略應(yīng)綜合考慮模型幾何尺度與物理場變化，需根據(jù)試算結(jié)果，對(duì)耦合作用較為劇烈的區(qū)域進(jìn)行網(wǎng)格細(xì)化。理論上講，只要不產(chǎn)生畸變，網(wǎng)格越細(xì)化，梯度解析越準(zhǔn)確，所求得的數(shù)值解便越逼近解析解。此外，在網(wǎng)格細(xì)化的基礎(chǔ)上，還需考慮整體網(wǎng)格單元的協(xié)調(diào)性與物理場的匹配性，以確保計(jì)算結(jié)果切合實(shí)際。

根據(jù)模型的幾何尺寸與試算結(jié)果，本文以最大單元大小為變量，共制定5種不同密度的網(wǎng)格剖分方案（由粗到細(xì)）。在相同工況下，對(duì)5種網(wǎng)格剖分方案進(jìn)行了計(jì)算。設(shè)置求解步長為1d，以便于整數(shù)天數(shù)據(jù)（非插值結(jié)果）的提取。

圖 2顯示了5種方案t=200d的計(jì)算結(jié)果?？梢钥闯?，除Mesh2外，其余4個(gè)方案的濃度場計(jì)算結(jié)果中，塌體底邊并未遭受顯著滲蝕。此外，Mesh2在該時(shí)步各未知量的表面平均值顯然并不符合其余4種網(wǎng)格的收斂趨勢。其原因?yàn)镸esh2的整體網(wǎng)格協(xié)調(diào)性與耦合物理場之間并不匹配。Mesh2計(jì)算結(jié)果的時(shí)間超前效應(yīng)是瞬態(tài)響應(yīng)被不協(xié)調(diào)網(wǎng)格放大的后果，是一種假象，嚴(yán)重偏離實(shí)際。由此可見，網(wǎng)格剖分應(yīng)綜合考慮網(wǎng)格密度、整體網(wǎng)格協(xié)調(diào)性與耦合物理場之間的匹配性。

事實(shí)上，對(duì)于剩余4種網(wǎng)格方案而言，最粗糙的Mesh1，其數(shù)值精度已經(jīng)收斂，但由于本文的重點(diǎn)在于數(shù)值解實(shí)現(xiàn)方法，因此選數(shù)值精度最高的網(wǎng)格剖分方案Mesh5用于討論。

3"結(jié)果與討論

為充分論證高階高度非線性強(qiáng)耦合PDEs數(shù)值求解的苛刻性與本文實(shí)現(xiàn)方法的可行性，本節(jié)將滲蝕強(qiáng)耦合PDEs的不同求解方法進(jìn)行了比較（網(wǎng)格方案均為Mesh5）。圖 3顯示了各求解方法的斂散情況及其求解日志的關(guān)鍵參數(shù)。對(duì)于滲蝕四場高階高度非線性強(qiáng)耦合PDEs，以非弱形式建模，利用全耦合算法求解，如圖3（a）所示，在Step=72時(shí)，收斂曲線首次發(fā)散。因無法滿足局域誤差估計(jì)與測試標(biāo)準(zhǔn)，導(dǎo)致求解計(jì)算在局域時(shí)間內(nèi)反復(fù)回溯，非線性求解器累計(jì)計(jì)算失敗次數(shù)（NLfail）持續(xù)上升。類似地，相同情形在分離式求解非弱形式模型過程當(dāng)中重現(xiàn)。在Step=96時(shí)，收斂圖走勢明顯反彈，如圖3（b）所示，由于非線性求解失敗記錄的累積，導(dǎo)致計(jì)算最終陷入惡性循環(huán)。根據(jù)求解日志及斂散情況，判斷非弱模型不會(huì)收斂，因此無需繼續(xù)運(yùn)算，在Step=373、404時(shí)，分別終止了非弱形式模型的兩種求解過程。顯然，對(duì)于高階高度非線性強(qiáng)耦合PDEs，非弱形式建模求解是無力的。

為此，本文將滲蝕強(qiáng)耦合PDEs進(jìn)行了弱化，采用全耦合求解與分離式求解對(duì)弱化模型進(jìn)行了運(yùn)算。由收斂曲線初期走勢可知，弱形式對(duì)于非線性的處理更為妥善。不難看出，弱形式建模結(jié)合分離式求解順利解決了滲蝕強(qiáng)耦合PDEs的非線性難題，模型收斂曲線最后基本穩(wěn)定在最大時(shí)間步長附近，且全過程非線性求解失敗累計(jì)次數(shù)為0，如圖3（e）所示。而未能遂愿的是，對(duì)于弱形式與全耦合的求解組合，在Step=781時(shí)收斂曲線明顯上揚(yáng)，如圖3（d），且計(jì)算結(jié)果已經(jīng)失真（φ＞1.0），可見，非線性強(qiáng)耦合PDEs求解條件要求苛刻，同時(shí)也意味著求解器選擇策略值得重視。

事實(shí)上，可將各種物理場進(jìn)行任意耦合的Comsol Mutiphysics在運(yùn)算非弱模型時(shí)，同樣對(duì)PDE進(jìn)行了自動(dòng)弱化處理，但不足之處在于其僅弱化了域方程，并未協(xié)調(diào)弱化耦合邊界條件。此外，在求解算法的選擇方面，全耦合求解法會(huì)根據(jù)初值條件對(duì)所有PDE同步執(zhí)行NewtonRaphson迭代，直到收斂，如圖4所示。該算法常被用于普通問題的收斂誤差控制，而滲蝕強(qiáng)耦合PDEs屬于瞬態(tài)研究，無法精確給定初始猜測，進(jìn)而導(dǎo)致了全耦合求解難以收斂。分離算法是在任意時(shí)步循序逐個(gè)求解PDE，如圖5所示。相較于前者，該算法既加快了解的收斂又節(jié)約了求解所需資源。每個(gè)分離步驟本身均可視為一個(gè)非線性問題，而該算法允許每個(gè)步驟中的每個(gè)PDE自主選擇最優(yōu)求解策略。盡管需要更多次數(shù)的迭代，但其每個(gè)環(huán)節(jié)均快于NewtonRaphson步驟。同時(shí)，該法直觀地顯示了各分離步驟的斂散情況，便于排除物理場的錯(cuò)誤設(shè)定。

需要指出的是，瞬態(tài)求解累計(jì)計(jì)算失敗次數(shù)（Tfail）僅影響求解效率，可通過適當(dāng)細(xì)化時(shí)間離散來改善。換言之，最終的斂散情況與結(jié)果的可信程度并不取決于Tfail。為說明求解方法的可靠性，本文擬合了累計(jì)細(xì)粒流失百分比隨時(shí)間的變化曲線，如圖6所示。顯然，累計(jì)細(xì)粒流失百分比隨時(shí)間符合具有時(shí)間常數(shù)參數(shù)的單相指數(shù)衰減函數(shù)，擬合結(jié)果在函數(shù)型式和表達(dá)形式上均高度吻合于以Sterpi[34]為代表的滲蝕本構(gòu)模型，驗(yàn)證了求解方法的可靠性。

4"結(jié)論

（1）強(qiáng)耦合是導(dǎo)致PDEs非線性特性的充分條件；

（2）非弱形式建模難以妥善解決高階高度非線性強(qiáng)耦合PDEs數(shù)值收斂性問題；

（3）分離式求解算法對(duì)于高階高度非線性強(qiáng)耦合PDEs的初始條件更具包容性；

（4）數(shù)值試驗(yàn)結(jié)果表明，文章所述求解方法可行性與可靠性兼得。

參考文獻(xiàn)：

[1]"孫培德，楊東全，陳奕柏.多物理場耦合模型及數(shù)值模擬導(dǎo)論[M].北京：中國科學(xué)技術(shù)出版社，2007.

[2]"魏海江.隧洞（道）充填介質(zhì)突水突泥滲蝕致災(zāi)機(jī)制研究[D].西安：西安理工大學(xué)，2024.

[3]"Sterpi D. Effects of the erosion and transport of fine particles due to seepage flow[J].International Journal of Geomechanics.2003，3（01）：111122.

[4]"魏海江，許增光，耿凱強(qiáng)，等.滲蝕過程中孔結(jié)構(gòu)參數(shù)對(duì)細(xì)顆粒流失的影響[J].排灌機(jī)械工程學(xué)報(bào)，2023，41（04）：417424，432.

責(zé)任編輯：肖祖銘

Numerical Solution Method for Highorder Highly NonlinearStrongly Coupled Partial Differential Equations：Taking the Strongly Coupled Partial Differential Equations Describing Suffusion as an Example

WEI Haijiang1， XUE Rui1， LIANG Gang2， CAO Cheng3， YANG Tian1， ZHANG Xinwei4

（1. Northwest Electric Power Design Institute Co.， Ltd. of China Power Engineering Consulting Group， Xi'an 710075， China;

2. School of Mechanical and Electronic Engineering， Jingdezhen University， Jingdezhen 333400， China;

3. State Key Laboratory of EcoHydraulics in Northwest Arid Region of China， Xi'an University of Technology， Xi'an 710048， China;

4. College of Safety Science and Engineering， Xi'an University of Science and Technology， Xi'an 710054， China）

Abstract：To solve the highorder highly nonlinear strongly coupled Partial Differential Equations （PDEs）， this paper takes the PDEs describing suffusion as a typical case， analyzes the highorder highly nonlinear of PDEs， combines spatial mapping， based on weak form modeling and segregated approach， realizes the numerical solution of PDEs describing suffusion， and verifies the feasibility and reliability of the solution method. The results show that strong coupling is a sufficient condition for PDEs to have nonlinearity. Without the help of weak form of PDEs， it will be difficult to solve the nonlinearity of strongly coupled models; Segregated approach is more inclusive to the initial guess of transient strongly coupled PDEs.

Keywords：multifield strongly coupled PDEs; highorder highly nonlinear; weak form; solution method