多孔壁面對高速邊界層最優增長條帶二次失穩的影響規律

王宇天，劉建新，王曉坤，李曉明，

1.空軍航空大學，長春 130022

2.天津大學 機械工程學院，天津 300072

高超聲速飛行器具有重要的軍事和經濟價值，是當今各航空航天大國的研究熱點之一。在臨近空間飛行雷諾數范圍內，邊界層通常要經歷由層流到湍流的轉捩過程。由于飛行器摩擦阻力、熱防護、發動機性能等均與邊界層流動狀態密切相關，因此邊界層轉捩是飛行器設計中必須要解決的基礎理論問題，具有重要的工程應用價值與學術意義［1-4］。

邊界層轉捩描述的是擾動產生與失穩演化直至突變為湍流（Breakdown）的全過程，其中擾動的產生由感受性機制決定，失穩演化則是流動穩定性理論的范疇。線性穩定性理論（Linear Stability Theory，LST）、拋物化穩定性方程（Linear Parabolized Stability Equations，LPSE）等可以預測模態擾動的線性發展過程，但人們在Poiseuille、Couette、Blasius 邊界層等流動實驗中發現，擾動在中性曲線以外的穩定區域依然出現了較大增長，甚至轉捩。由于擾動間的非線性作用只是擾動能量的重新分配，并不能產生擾動能量的凈增加，所以必然存在一種不同于模態增長的線性增長機制，通常稱為非模態或瞬態增長。通過數學分析可知，瞬態增長由線性擾動方程算子矩陣的非正規性所致［5］。由于非正規矩陣的特征向量并不彼此正交，因此即使所有特征向量均為衰減模態（主要為連續模態），也可通過一定的線性組合產生增長擾動。而在物理上，瞬態增長通常由Lift-up 機制［6］進行解釋：在流向渦的作用下，壁面附近的低速流體被拋向邊界層外，并將邊界層外的高速流體帶向壁面，進而形成了流向速度在展向呈周期性分布的條帶結構。

瞬態增長中各衰減模態的線性組合系數，即幅值與相位信息，同樣由感受性過程決定，所以瞬態增長率受外部擾動情況（如高幅值自由流擾動、粗糙元等）的影響。但在特定的流場參數條件下，可以通過最優擾動理論進行求解，得到在有限時間或空間范圍內實現最大瞬態增長的最優擾動。已有的研究表明［7-8］，從不可壓到高超聲速邊界層，最優擾動通常為定常流向渦的形式，其通過非線性發展增長至較大幅值的條帶結構時，新的平均流剖面在展向上將出現拐點，從而發生二次失穩觸發轉捩。因此從轉捩預測的角度上看，這是一種更為保守的亞臨界轉捩觸發閾值。

瞬態增長理論完善了傳統的線性穩定性理論框架，可以對超越了模態增長理論的旁路轉捩中一些亞臨界轉捩現象進行定性解釋，如經典的Poiseuille 管道流動［9］。除此之外，瞬態增長理論的一個重要應用是再入返回艙等大鈍度飛行器的“鈍頭體悖論”問題，即發生在弓形激波后亞聲速區的轉捩現象。根據線性穩定理論，此處的順壓梯度作用將抑制Tollmien-Schlichting（T-S）波的增長，而凸面曲率也將抑制G?rtler 渦，同時橫流速度較小使得橫流模態可以忽略，因此人們在很長一段時間內無法找到實驗中出現轉捩的原因，并且近期在李強等［10］的頓頭平板轉捩實驗中也出現了這一現象。近幾年來，NASA 課題組Paredes 等［11］對最優擾動進行了細致的研究，發現在駐點附近存在較強的瞬態增長，并認為壁面微小粗糙度便可能成為誘導轉捩的主要原因。同時瞬態增長理論也同樣能夠對鈍頭體“轉捩反轉”現象［12］進行定性解釋［13-14］，計算發現在鈍頭與錐身連接處附近存在較強的瞬態增長，并且增長率隨著鈍度增大而增大，因此能夠在Mack 第一、第二和熵層模態增長率較小的情況下觸發轉捩。盡管連接處的微小粗糙元不一定能夠恰好激發最優擾動，但也為這一轉捩過程預測提供了一種科學的閾值指導，并且還需要進一步的實驗驗證。

隨著對轉捩機理的深入認識，各種延遲轉捩技術應運而生［15］，其中多孔壁面能夠有效抑制Mack 第二模態，可作為一種推遲轉捩的被動控制方式［16］。目前，微孔結構可在高超聲速飛行器表面涂層材料的制備工藝中自然形成，因此在工程領域受到了廣泛關注。但已有的研究均針對二維邊界層的局部穩定性問題，對于三維邊界層二次失穩模態的控制效果還不明確。本文通過建立最優擾動的求解方法，并以最優增長條帶形成的三維邊界層為研究對象，研究多孔壁面對全局穩定性（Bi-Global）的影響規律，評估其在亞臨界轉捩中的控制效果，以期為多孔涂層的布置方案提供一定的科學指導。

1 數值方法及驗證

本文以超聲速絕熱平板邊界層流動為研究對象，基本控制方程為直角坐標系下的可壓縮Navier-Stokes（N-S）方程，選取以下特征量與特征尺度對物理與時空變量進行無量綱化：

式中：x、y、z分別為流向、法向與展向坐標；u、v、w分別為流向、法向與展向速度；t為時間；ρ為密度；T為溫度；p為壓力；上標“*”代表有量綱物理量；下標“∞”代表自由來流值。特征長度取0.001 m，故參考雷諾數定義為。黏性系數μ由Sutherland 定律求得，傳熱系數κ則可通過普朗特數Pr得到：

流動穩定性分析采用原始變量，在非守恒型N-S 方程的基礎上，通過理想氣體狀態方程p=消去壓力項p。將瞬時物理量q=[ρ，u，v，w，T]T表示為基本流與擾動量之和：q=+q′，將其代入N-S 方程后減去基本流滿足的方程部分，并忽略非線性作用項，便可得到如下緊致形式的線性擾動方程：

式中：系數矩陣Γ、A、B、C、D、Vxx、Vyy、Vzz、Vxy、Vxz、Vyz的具體表達式可參見文獻［17］，并以可壓縮邊界層相似解作為基本流動［18］。為了理論方法介紹的完整性，LPSE 及其伴隨方程的表達式詳見附錄A［19］。

1.1 最優擾動方法

傳統的LST 通過求解常微分方程特征值問題，得到的是擾動在時間與空間趨近于無窮時的漸進演化行為。而瞬態增長則通過求解常微分方程的初值問題，得到擾動在有限時間或空間內的代數增長。目前，最優擾動有2 種求解方法：一種是在平行流假設條件下，基于LST 方程的奇異值分析［20-21］；另一種是考慮邊界層弱非平行性的，基于伴隨方程的優化方法。基于LST 的平行流最優擾動分析只能定性地給出最優擾動的展向波數范圍，而對于擾動放大倍數G的計算可能存在較大誤差，因此本文將在考慮邊界層弱非平行性的LPSE 方程的基礎上，利用Lagrange 乘數法求解擾動增長的數學最優化問題［22-23］。

最優增長條帶在流向上的尺度較大，因此擾動幅值與相位信息的變化可以全部反映在形狀函數中，即可使式（A1）中α=0。采用Hanifi等［20］建議的擾動能量范數E的計算準則，以內積的形式給出：

設入口x0處最優初始擾動為，以某一指定位置x1處的擾動的能量范數的增長倍數作為目標函數：

其中：?代表取梯度。為使式（8）恒滿足δF=0，式（8）中兩項內積需同時為0。第1 項自動滿足，即求解齊次LPSE 方程（式（A2））：

式（8）中第2 項內積可利用內積定義（式（A5））與分步積分將其變換為

式中：L?=0即為齊次伴隨LPSE方程，此時分步積分產生的邊界項僅保留式（10）中的第二、第三項。式（10）進一步變換為

要使式（11）恒等式成立，等號左邊兩項需同時為0，這樣便可以得到在x0與x1處初始擾動與伴隨向量的關系式為

對于線性問題，c0、c1為歸一化系數使得。在壁面處為奇異矩陣，因此可由法向動量方程求得。在入口位置x0處，式（12）中的流向擾動速度優化條件為

由于伴隨密度擾動在壁面處無預設齊次條件，將導致擾動速度在壁面處非零，本文采用Tempelmann 等［23］的建議方法，直接假設，下標“w”代表在壁面處的值，該方法與Zuccher 等［24］的插值方法具有相同的效果。基于以上優化條件，便可通過迭代求解得到最優擾動，具體步驟為

步驟1選取任意滿足齊次邊界條件的擾動作為初始擾動，可利用LST 算子進行當地求解。

步驟2利用LPSE 方程（式（A2））從x0→x1推進求解得到。

步驟3由式（13）得到伴隨向量，利用伴隨LPSE 方程從x1→x0推進求解得到。

步驟4由式（12）與式（14）得到新的初始擾動0。

重復步驟2～步驟4，直至目標函數J()=G收斂。計算結果表明該優化系統具有吸引子，通常迭代3～4 步便可以收斂至10-4。

以x1=1 m 為優化目標位置，在Ma=3.0 工況下，當展向波數β=0.25 時取不同的初始位置x0，或當x0=0.36 m 時取不同展向波數β，最優擾動增長倍數的計算結果如圖1 所示，其中縱坐標以來流單位雷諾數Re∞做歸一化處理。通過與Paredes 等［8］的結果對比驗證了程序的正確性，并且可以看出最優擾動增長倍數與優化初始位置的選取密切相關。

圖1 最優擾動增長倍數的驗證Fig.1 Validation of optimal disturbance amplification

1.2 全局穩定性分析

不同于局部LST 分析，對于基本流在展向z同樣是快變，而流向x仍是慢變的三維邊界層，需采用全局穩定性分析（Bi-Global）。此時擾動可以寫為如下行進波的形式：

方程在法向上同樣采用10 階中心差分格式，展向上采用適用于周期性邊界條件的Fourier 譜方法離散，結合齊次邊界條件（式（A4）），便可通過求解廣義特征值問題得到二次失穩擾動的空間增長率。對于大型稀疏矩陣，本文采用隱式自啟動Arnoldi 方法。

對于多孔壁面，由于孔隙直徑一般為幾十微米量級，孔隙間的相互作用較小，因此其誘導的流向與橫向擾動速度可忽略不計，但對于法向速度與溫度擾動具有半透射作用。Fedorov 等［25］通過引入聲導納系數將其與壓力擾動相關聯，擾動在壁面處的邊界條件變為

式中：Ac、Bc為復導納系數，其與涂層材料、平均流以及擾動特征相關。由于涂層厚度（孔隙深度）h通常遠大于孔隙直徑2r，可將其看作細長管。基于聲波在細長管內的傳播理論，推導出復導納系數的計算公式為

式中：?為開孔率，當孔隙截面為圓形時，易知最大開孔率?→π/4；下標“w”代表基本流動在壁面處的值；Jn則為n階第一類Bessel 函數。

已有的直接數值模擬（Direct Numerical Simulation，DNS）與實驗結果均證明該理論邊界條件的合理性［26］。本文選取Fedorov 等［25］針對二維邊界層的計算結果進行驗證，流場條件為：Ma∞=6.0，=243.9 K，=1.4，計算結果如圖2 所示，驗證了程序的正確性。

圖2 多孔壁面條件下第二模態增長率的驗證Fig.2 Validation of growth rate of the second modes over porous wall

1.3 直接數值模擬

對于最優增長擾動的非線性演化過程，本文采用高精度有限差分方法求解N-S 方程得到。計算程序在中國科學院李新亮等的開源程序OpenCFD-SC［27］的基礎上發展而來。當采用層流邊界層相似解作為基本流時，由于相似解并不是N-S 方程的穩態解，因此需要在N-S 方程中引入人工源項以維持基本流保持不變［28］。在對無黏通量進行空間離散時，采用截斷誤差最小的Steger-Warming 通量分裂格式得到正負通量，再采用7 階迎風格式進行離散，黏性通量采用8 階中心差分格式進行離散。時間推進采用具有TVD（Total Variation Diminishing）性質的3 步3階Runge-Kutta 方法。在壁面處，速度與溫度采用無滑移等溫條件，并假設?p/?y=0 作為壁面壓力條件。

2 計算結果與分析

已有的研究表明，多孔壁面通常只對無黏的Mack 第二模態起穩定作用。但根據3 層結構理論［29］，滿足的Mack 第一模態同樣為無黏模態，但多孔壁面卻起促進作用。為了研究多孔壁面在不同中高馬赫數范圍內對二次失穩模態的影響規律，本文將研究來流馬赫數為3.0 與6.0 這2 種情況，來流總溫=333 K，Pr=0.7，單位雷諾數Re∞=106m-1。

下面將首先利用伴隨方法得到線性最優擾動，然后通過DNS 計算不同初始幅值最優擾動的非線增長過程得到最優增長條帶，再以某一流向位置的y-z截面作為新的三維邊界層剖面進行二次失穩分析。

由于條帶基本流在展向上具有周期性，對于基頻（Fundamental）或亞諧型（Subharmonic）二次失穩模態，可以根據其對稱或反對稱性質分為奇（Sinuous 型）、偶（Varicose 型）模態，二者具有不同的增長率。偶模態的擾動分量在一個周期內呈對稱分布，為反對稱分布，而奇模態的對稱性與之相反。本文在實際計算中，為更好地區分奇/偶模態，展向上采用對稱離散方式，即只需保留一半的展向網格即可。通過網格分辨率驗證，Ny×Nz=201×64 已滿足計算需求。

2.1 最優擾動

圖3 分別為Ma∞=3.0 與Ma∞=6.0 流場條件下，不同起始位置的最優擾動增長倍數隨展向波數β的變化規律。可以看出在高超聲速條件下，最大瞬態增長倍數相對降低，并且產生最大瞬態增長的展向波數β減小。

圖3 初始位置x0 與展向波數β 對最優擾動增長倍數的影響（x1=1 m）Fig.3 Effect of initial position x0 and spanwise wavenumber β on optimal disturbance amplification（x1=1 m）

由圖3 可以看出，選取不同的初始位置將改變擾動的最優增長倍數。從感受性的角度看，由于前緣區域的邊界層較薄且具有較強的非平行性，因此自由來流渦擾動或壁面粗糙度具有更大的感受性系數，因此由較高幅值環境擾動所引起的瞬態增長通常發生在中性曲線下支前。為不失研究對象的一般性，本文選取靠近前緣的初始位置，即x0=0.09 m，同時在Ma∞=3.0 與Ma∞=6.0 流場條件下，分別選取展向波數β=0.25 與β=0.10 的最優擾動作為初始擾動，初始與終點位置處擾動的形狀函數如圖4 所示。可以看出初始最優擾動具有流向渦的形式，也就是展向與法向擾動速度較大，其增長由Lift-up 機制主導。并且，隨著馬赫數增大至高超聲速，最優初始擾動的溫度和密度擾動分量與橫向和法向擾動速度具有相同量級。當流向渦發展為條帶時，溫度與密度擾動將超過流向速度擾動成為擾動的主要分量，這與定常G?rtler 模態相似［30］。

圖4 初始與終點位置處的最優擾動形狀函數（x0=0.09 m，x1=1 m）Fig.4 Shape function of optimal disturbances at initial and final locations（x0=0.09 m，x1=1 m）

2.2 馬赫數3.0 工況

當最優擾動初始能量（式（4））E0=0.000 4時，x=0.2 m 處的截面條帶如圖5 所示。圖中深灰色線為流向速度等值線，其大小從0.1 依次遞增0.1 至1.0，可以看出此時展向流場已具有一定的非平行。

圖5 x=0.2 m 處的流向速度等值線與梯度云圖（E0=0.000 4）Fig.5 Contours of isolines and gradient of streamwise velocity at x=0.2 m（E0=0.000 4）

多孔壁面條件設為：開孔率?=0.5，孔隙半徑r=0.1 mm，約為當地邊界層厚度的2%。研究表明［26］：當孔隙深度大于0.8 倍的當地邊界層厚度時，控制效果基本不再隨孔隙深度而改變。設置孔隙深度h=2 mm，在本文中的超聲速條件下|Λh|＞3，根據式（19）中雙曲正切函數tanh 的性質可知tanh(Λh)→1，因此其與Λh→∞時等效。圖6 為基頻與亞諧頻失穩模態的增長率-αi隨角頻率ω的變化規律，以及一維基本流的LST增長率。當最優擾動條帶的幅值較小時，二次失穩模態的增長率未超過第一模態，并且基頻偶模態（FV）增長率高于基頻奇模態（FS），而對于亞諧型失穩則是奇模態（SS）高于偶模態（SV）。與第一模態類似，多孔壁面促使二次失穩模態增長率變大，并且對亞諧失穩的作用更強。

圖6 有無多孔壁面作用下的二次失穩模態增長率（E0=0.000 4）Fig.6 Growth rates of secondary instability modes with and without porous wall（E0=0.000 4）

圖7 ω=0.06 時二次失穩模態流向擾動速度模值Fig.7 Modulus of streamwise disturbance velocity of secondary instability at ω=0.06

當提高最優擾動的初始能量至E0=0.01時，x=0.2 m 處截面條帶如圖8 所示。此時條帶基本流進一步扭曲，展向剪切增大了10 倍以上。并且可以注意到等值線在中心處出現下凹區，這在G?rtler 條帶［30］中并未出現，將產生新的二次失穩模態。

圖8 x=0.2 m 處的流向速度等值線與梯度云圖（E0=0.01）Fig.8 Contours of isolines and gradient of streamwise velocity at x=0.2 m（E0=0.01）

圖9 為二次失穩模態的增長率-αi隨角頻率ω的變化規律，此時二次失穩模態的增長率遠高于第一模態，將成為促發轉捩的主要因素。與不可壓Blasius 邊界層［31］最優擾動條帶及G?rtler渦［30］的二次失穩相同，亞諧模態與基頻模態具有相近的增長率，且低于基頻模態，因此在轉捩預測時通常更關心基頻失穩。與圖6 相似，基頻奇模態（FS）增長率高于基頻偶模態（FV），而亞諧型失穩模態則與之相反。

圖9 有無多孔壁面作用下的二次失穩模態增長率（E0=0.01）Fig.9 Growth rates of secondary instability modes with and without porous wall（E0=0.01）

圖10 ω=0.085 時二次失穩模態流向速度擾動實部Fig.10 Real parts of streamwise disturbance velocity of secondary instability modes at ω=0.085

2.3 馬赫數6.0 工況

在進行Bi-Global 分析之前，首先對高超聲速二維基本流動進行穩定性分析。圖11 為x=0.2 m 處LST 分析得到的Mack 第二模態增長率-αi隨角頻率ω的變化規律。多孔壁面條件設為：開孔率?=0.5，孔隙半徑r=0.15 mm，約為當地邊界層厚度的1.5%，孔隙深度h=2 mm。從圖中可以看出：多孔壁面對低頻模態（即第一模態）有促進作用，而對高頻模態（即第二模態）有明顯的抑制作用，并且對于三維模態具有相同的控制效果。在該工況下，第二模態由慢模態［32］（即第一模態）發展而來，圖12 為二維快/慢模態相速度c=ω/αr隨角頻率ω的變化規律，多孔壁面對相速度的影響整體較小，并且對慢模態的影響主要體現在第一模態頻率區，對第二模態無明顯影響。快/慢模態的同步點頻率［32］約為0.237，而圖11 中控制效果的轉折點頻率約為0.22，二者較為接近，并且對于三維模態具有相同的規律。

圖11 x=0.2 m 處LST 分析得到的Mack 模態增長率Fig.11 Growth rates of Mack mode analyzed by LST at x=0.2 m

圖12 x=0.2 m 處LST 分析得到的快/慢模態相速度Fig.12 Phase velocity of fast/slow modes analyzed by LST at x=0.2 m

當擾動初始能量E0=0.004 時，x=0.2 m處截面條帶如圖13 所示，圖中深灰色線為流向速度等值線，其大小從0.1 依次遞增0.1 至1.0，流向速度法向梯度與展向剪切與Ma∞=3.0 工況具有相似的分布（圖5）。

圖13 x=0.2 m 處的流向速度等值線與梯度云圖(E0=0.004)Fig.13 Contours of isolines and gradient of streamwise velocity at x=0.2 m (E0=0.004)

高超聲速條帶的二次失穩情況更加復雜，由于基頻失穩通常在轉捩中占主導地位，下面將重點研究基頻失穩。該流場條件下，共出現2 個奇模態（S1，S2），3 個偶模態（V1，V2，V3）。光滑壁面條件下，以流向擾動速度模值||的最大值做歸一化處理，圖14 為二次失穩模態的流向擾動速度模值||及溫度擾動模值||的分布云圖。對于高超聲速流動，溫度擾動幅值遠高于速度擾動，并且速度擾動峰值靠近壁面區，而溫度擾動峰值位于靠近邊界層外緣的臨界層。

圖14 ω=0.28 時二次失穩模態流向速度與溫度擾動模值Fig.14 Modulus of streamwise velocity and temperature disturbances of secondary instability modes at ω=0.28

圖15 為二次失穩基頻奇/偶模態的增長率-αi隨角頻率ω的變化規律。V1 模態增長率高于S1 模態，是二次失穩的主導模態，而S2 與V3增長率較小可以忽略。V2 模態是一個較為特殊的模態，當ω＜0.27 時，其增長率明顯小于V1 模態與S1 模態，但其高頻區的增長率較高，可能成為轉捩中的最危險的模態。當采用多孔壁面邊界層條件時，其影響規律與Mack 模態相似，即只對高頻二次失穩模態起抑制作用，而對低頻二次失穩模態起促進作用。S1 模態與V1 模態的轉折頻率分別為0.215 與0.19，這與在圖11 中Mack模態的轉折頻率相近。V2 模態的轉折頻率為0.28，遠高于其他模態，但V2 的快速增長只能出現在高頻區域，相對的低頻區域即使有促進作用，由于其增長率遠小于S1 模態與V1 模態，因此不會對轉捩產生影響。

圖15 有無多孔壁面作用下的二次失穩模態增長率（Ma∞=6.0）Fig.15 Growth rates of secondary instability modes with and without porous wall（Ma∞=6.0）

綜合以上研究結果可以看出，雖然二次失穩模態在本質上都屬于無黏失穩，特別是對于奇模態，其應是由展向剪切主導的失穩模態，但在不同的流場條件下，多孔壁面對其影響規律卻不盡相同，并與無條帶時的局部穩定性特征密切相關。Song 等［33］針對G?rtler 渦的二次失穩問題，證明了二次失穩模態可以由上游的Mack 第一、第二模態演化而來，由此可以推斷二次失穩模態也將保留二維邊界層中Mack 模態的穩定性特征。但在物理空間中，上游某一特定頻率的模態擾動演化到下游也應為一特定擾動模態，但在Bi-Global 分析中卻會得到多個二次失穩模態，這其中的關聯機制還尚不清楚，需要借助模態分解技術［28］做進一步詳細研究。

3 結論

本文針對瞬態增長理論，利用Lagrange 乘數法與變分法，推導出以伴隨拋物化穩定性方程為約束條件的優化系統，建立并實現了最優擾動的數值求解方法。以可壓縮平板邊界層流動為研究對象，針對Mack 第一、第二模態分別主導轉捩的超聲速（Ma∞=3.0）與高超聲速（Ma∞=6.0）2 種工況開展研究。計算表明：最優擾動具有流向渦的形式，并通過Lift-up［6］增長機制發展為條帶結構，并且在高超聲速條件下，其與定常G?rtler 模態相似［30］，溫度與密度擾動將超過速度擾動成為擾動的主要分量。

以有限幅值最優擾動非線性發展形成的三維邊界層為研究對象，并利用Fedorov 等［25］建立的等效邊界條件，研究了多孔壁面對最優增長條帶二次失穩模態的影響規律。結果表明：與第一模態相同，多孔壁面對超聲速條帶的二次失穩擾動只起促進作用。對于高超聲速條帶，多孔壁面對第一模態頻率范圍內的二次失穩擾動為促進作用，對第二模態頻率范圍內的二次失穩擾動起抑制作用，并且轉折頻率接近局部快/慢模態的同步頻率。雖然目前還無法對其中的物理機制做出解釋，但對工程應用中多孔涂層的布置方案具有一定的指導意義。

附錄A：

基于邊界層慢變假設，將擾動分解為快變的波數函數和慢變的形狀函數兩部分：

擾動在邊界處滿足Dirichlet 邊界條件：

在法向上采用10 階中心差分離散后，結合齊次邊界層條件（式（A4）），此時便將擾動方程求解由LST 的特征值問題轉變為偏微分方程的初邊值問題，在流向上采用隱式歐拉格式推進求解即可得到擾動的空間發展過程，詳細計算方法與殘余橢圓性問題可參見文獻［19］。

伴隨LPSE 方程可以通過Green-Lagrange恒等式推導而來，首先定義向量內積為

式中：上標“H”表示共軛轉置。本文以上標“?”表示相應的伴隨變量，將伴隨向量與齊次LPSE 方程（式（A2））作內積，通過分步積分可以得到如下Green-Lagrange 恒等式：

式中：L?=0 即為伴隨LPSE 方程；為分步積分產生的邊界項，伴隨向量同樣滿足齊次邊界條件（式（A4））。伴隨算子L?為

伴隨LPSE 方程求解采用與LPSE 相同的數值方法，只是改為由下游向上游推進求解。