李新磊 吳 坤 趙林英 范學軍,

* (中國科學院力學研究所高溫氣體動力學國家重點實驗室,北京 100190)

? (中國科學院大學工程科學學院,北京 100049)

** (西北工業(yè)大學機電學院,西安 710072)

?? (合肥中科重明科技有限公司,合肥 230601)

引言

超燃沖壓發(fā)動機作為一種高超聲速巡航飛行器,受氣動加熱和燃料的燃燒釋熱影響因而承受著巨大的熱載荷[1-3].在發(fā)動機的燃燒室壁中,通常嵌入矩形直通道構成再生冷卻系統(tǒng),以保證發(fā)動機能夠長時間穩(wěn)定運行.然而,隨著飛行馬赫數的提升,傳統(tǒng)的再生冷卻結構面臨著冷卻能力不足、結構超溫等問題[4-6].因此,亟需開展新型熱防護方案的設計及優(yōu)化.

國內外學者對此進行了大量研究,其中最為常見的一類做法是添置粗糙元件.如Sunden 等[7]通過在圓管內設置非對稱的波浪形翅片,使得管道的換熱系數顯著提高,壁面溫度顯著降低.Li 等[8]對帶有截斷肋的矩形通道進行了研究,分析了截斷肋周圍的強化傳熱機制,并探究了肋高和肋間距等參數對通道流動換熱性能的影響規(guī)律.楊澤南等[9]對由Kagome 型網格陣列填充的通道進行了研究,發(fā)現在相同工況下,錯排網格填充的通道,不僅減重效果明顯,網格后復雜的渦系結構也使得流體努塞爾數明顯提高,壁面溫度顯著下降.此外,部分研究還通過改變再生冷卻通道的流動模式來達到增強換熱的目的.如Li 等[10]基于變分法生成了能夠自適應幾何和熱邊界條件的變管徑冷卻通道,在犧牲部分壓降的條件下使得壁面平均溫度和溫度不均勻度均有所下降.Li 等[11]對嵌套式冷卻通道進行了研究,發(fā)現同向平行通道有著更好的傳熱性能,且內通道的嵌入有助于克服熱分層現象.Zhang 等[12]提出一種雙向流動模式,該模式使得壁面溫度分布更加均勻,正癸烷的轉化率提升以及較大的化學熱沉,但壁面最高溫度和壓降均有不同程度的增加.

除了以上研究,近年來以拓撲優(yōu)化為代表的新興優(yōu)化方法,逐漸成為散熱系統(tǒng)設計的新途徑[13-16].該方法是以微結構的材料性質作為設計變量,在固定的設計空間內自動尋找最優(yōu)的材料分布方案[17],因此在理論上擁有最高的設計自由度,被認為是最具前途的優(yōu)化方法之一[18].Tang 等[14]曾對考慮材料變物性的導熱熱沉結構進行拓撲優(yōu)化,得到了如同樹枝枝干般的熱沉結構,同時發(fā)現固體材料的變物性顯著改變了拓撲優(yōu)化構型,進而影響了其在真實熱環(huán)境中的性能表現.Zhang 等[19]基于偽密度法對沖壓發(fā)動機的再生冷卻結構進行了拓撲優(yōu)化,發(fā)現其設計域內生成了眾多垂直于主流方向的微通道,顯著改善了整體換熱,但通道內局部的逆流和滯留區(qū)導致出現了局部高溫.本文認為,針對發(fā)動機再生冷卻通道的拓撲優(yōu)化是結合了湍流流動、燃料物性劇烈變化等特征的復雜流熱耦合問題[19],其計算量大且極易出現網格依賴性布局、棋盤格單元等構型缺陷,進而影響結構性能.針對上述問題,本文將開發(fā)一套考慮變熱物性的流熱耦合拓撲優(yōu)化求解器,采用建表插值法求解熱物性,用以緩解熱物性求解計算量大的問題.此外,將完整再生冷卻設計域簡化為一個換熱設計單元,通過適當的濾波技術和松弛求解策略,以盡量避免出現構型缺陷.

本文主要內容如下,首先建立拓撲優(yōu)化模型,基于連續(xù)伴隨法對考慮變物性的伴隨方程和靈敏度進行推導,并基于開源CFD 平臺OpenFOAM[20]構建拓撲優(yōu)化求解器.針對流熱耦合問題開展拓撲優(yōu)化設計,重點考察能量耗散約束對優(yōu)化構型的影響.為了評估拓撲優(yōu)化構型的流動換熱性能,本文還將對拓撲優(yōu)化構型輪廓進行提取,開展三維共軛換熱數值模擬分析.

1 拓撲優(yōu)化模型

拓撲優(yōu)化方法種類眾多,包括偽密度法、水平集法、相場法以及進化結構優(yōu)化法等[18].在各方法中,以偽密度法的應用最為廣泛和成熟.因此,本節(jié)將基于偽密度法建立拓撲優(yōu)化模型,并利用連續(xù)伴隨法對考慮變物性的伴隨方程及靈敏度進行推導.

1.1 控制方程

針對流熱耦合的拓撲優(yōu)化問題,其物理模型示意圖如圖1 所示,設計域承受來自壁面的熱載荷q或由內部產生的體積熱Q,流體在流經該設計域時與固體實現熱交換,進而達到熱平衡.該優(yōu)化問題的核心是在滿足壓降或材料體積比等約束下,找到最優(yōu)的流固材料單元的空間分布,從而使區(qū)域的平均溫度或最高溫度等目標降到最低.

圖1 拓撲優(yōu)化設計示意圖Fig.1 Schematic of the topology optimization design

首先建立流熱耦合拓撲優(yōu)化模型,該模型的控制方程包含如下的質量方程、動量方程和能量方程

其中,R(w,γ) 表示等式約束向量,γ 為設計變量,取值范圍為0～1,w=(p,u,T) 為隱式依賴于設計變量的狀態(tài)向量,在該問題中分別對應壓力、速度和溫度;ρ,Cp,μ,λ 和h分別為密度、定壓比熱容、動力黏度、熱導率和比焓,均考慮其溫度的相關性.τ為黏性應力張量,其定義式為τ=2μeffS(u)=(μ+μt)·為有效動力黏度,表示分子動力黏度μ與湍流黏度μt之和,Prt為湍流普朗特數,默認值為0.85,Q為體積熱源.

與常規(guī)形式N-S 方程不同的是,拓撲優(yōu)化模型的動量方程中引入了體積力項 -α(γ)u,α (γ) 是與多孔介質滲透率的倒數相關的插值函數.引入設計變量 γ 對 α (γ) 進行插值,用以區(qū)分流體與固體材料單元,從而實現相似多孔介質流動的描述,因此該變量又稱為“偽密度”.插值函數采用如下形式[21]

式中,θ 為形狀參數,取值范圍為0.005～0.1,該函數形狀如圖2 所示.αmax的取值則與達西數Da相關

圖2 不同形狀參數下的插值函數Fig.2 Interpolation function with different shape parameterθ

式中，l為幾何特征尺寸,Da表征了多孔介質的滲透能力,其取值小于 10-5較為合適[15].

相應地,在能量方程中對熱導率進行插值處理,λmin和 λmax分別為流體和固體材料的熱導率.當設計變量 γ=0 時,α 取值為一個較大數 αmax,表明該處單元為滲透率較低的固體單元;當設計變量 γ=1 時,α 取值為0,表明該處為流體單元.

本文采用了比焓形式的能量方程,與定壓比熱容存在如下熱力學關系

因此,在求解溫度時,可利用如下形式的牛頓迭代法計算得到

式中,pold,Told,Cp(pold,Told) 和h(pold,Told) 分別為上個迭代步中的壓力、溫度、比熱和比焓,hnew為利用式(1)計算得到的更新后的比焓,該迭代在相對誤差小于10-4時終止.

除了上述控制方程,本文定義了如下形式的Pnorm 函數[22]作為該優(yōu)化問題的目標函數

該目標有利于避免設計域內局部區(qū)域溫度過高,式中n取值為30.此外,引入如下的不等式約束函數

其中,G(w,γ) 為不等式約束向量,分別用來約束設計域的材料體積比例以及進出口能量耗散.Vmax為材料體積比例的期望值,φmax為進出口能量耗散的期望值.

綜上所述,一個完整的基于密度法的拓撲優(yōu)化問題可以由以下數學表達式來概述.

1.2 伴隨方程與靈敏度推導

本節(jié)將基于連續(xù)伴隨法對流熱耦合問題的伴隨關系式及靈敏度公式進行推導,首先構造如下形式的拉格朗日函數

式中,Ψ 為目標函數,λ=(pa,ua,Ta) 為拉格朗日乘子,又稱伴隨向量.由于R(w,γ)=0,因此原目標函數關于設計變量的靈敏度等價于如下拉格朗日函數的全導數

依據偏微分方程約束的KKT (Karush–Kuhn–Tucker)條件[23],構造如下的伴隨變量方程組

利用積分變換對上述方程組進行化簡,便可得到求解伴隨變量所需的控制方程,其具體形式及推導過程參見附錄A.能夠發(fā)現,由于考慮了溫度相關的變熱物性和輸運性質,伴隨變量的控制方程形式相比于常物性模型[15]復雜很多.而伴隨變量控制方程與原始狀態(tài)變量的控制方程存在形式上的相似性,因此可采用相似的數值算法對其進行求解.

此時拉格朗日函數的全導數,即目標函數的靈敏度,經推導后轉化為如下形式

在求解得到原始狀態(tài)變量與伴隨變量解后,便可計算得到該靈敏度數值.

2 數值方法

本節(jié)將對拓撲優(yōu)化過程中遇到的關鍵數值問題以及所采取的數值方法進行介紹,包括濾波與投影函數的使用、熱物性的求解、拓撲優(yōu)化求解器的構建以及共軛換熱數值求解器的驗證等.

2.1 濾波與投影

在拓撲優(yōu)化模型中,由于對材料物性采用了特殊的懲罰策略,導致優(yōu)化過程中容易出現一些不切實際的布局,如網格依賴性布局或棋盤格單元等[24-26],這嚴重影響了材料布局的重構和收斂.研究表明,一些正則化技術,如密度濾波器[27]可以有效地解決上述問題.本文采用亥姆霍茲型偏微分方程作為濾波器,其方程形式如下

式中 γ0為濾波前的偽密度,γf為濾波后的偽密度,Rmin為濾波半徑.當網格單元尺寸為 ?x時,濾波半徑Rmin越大,則 γf的形狀特征越容易被抹平,因此,本文控制各算例的最大濾波半徑不超過 2 ?x.

此外,采用濾波器通常會產生灰色單元,即處于流固材料單元之間的模糊單元.因此,可以配合使用投影函數[28]來獲得清晰的流固邊界.本文采用如下的雙曲正切函數來對 γf進行投影

式中 γp為投影后的偽密度,β 和 η 為用來控制投影函數形狀的形狀因子.投影函數的 η 取值為Vmax,用來保證材料的體積比例,而控制斜率的 β 因子隨優(yōu)化過程進行而逐漸變大,本文中 β 初始值為1,最大值為32,投影后的優(yōu)化變量 γp近乎收斂到接近于{0,1}的集合,可以得到較為清晰的流固邊界.

在對設計變量 γ 進行濾波和投影操作后,此時的靈敏度形式為

為了獲得原始靈敏度,需要進行如下的鏈式求導[15]

2.2 熱物性求解

在發(fā)動機冷卻通道內,碳氫燃料的熱物性隨溫度與壓力在很大的范圍內變化[29].以正癸烷為例,作為一種典型的碳氫燃料,在再生冷卻的流動吸熱過程中,其自身溫升高達數百攝氏度,熱力學性質和輸運性質會在臨界點附近發(fā)生陡峭變化,容易引發(fā)方程求解過程中的數值積分剛性問題.此外,在利用連續(xù)伴隨方法推導得到的伴隨方程中,存在熱力學性質及輸運性質關于溫度的偏導項,包括 ?h/?T,?ρ/?T,?Cp/?T,? μ/?T和 ? λ/?T,當采用解析形式的物性求解方程時,偏導項的求解將變得十分困難.

因此,本文對物性參數及各偏導項的計算采取了建表-插值[30]的方法,首先基于CoolProp 軟件[31]建立數據表,在不考慮燃料高溫裂解的前提下,該物性數據表以壓力和溫度作為基礎變量,變量間隔分別為1 kPa 和2 K,模擬過程中采用線性插值計算得到相關的參數.圖3 所示為正癸烷在3 MPa 壓力下(Tpc=648.4 K)的密度、比熱以及二者關于溫度的熱力學偏導性質,其中利用建表-插值得到的數據與利用CoolProp 直接計算得到的數據十分接近,最大相對誤差均不超過0.2%,表明本文所采用的物性計算方法是準確的.

圖3 熱物理性質及熱力學偏導Fig.3 Thermo-physical properties and derivatives

2.3 求解器構建

基于上述的模型與數值方法,本文利用開源平臺OpenFOAM 搭建了拓撲優(yōu)化求解器.圖4 所示為求解器的主要求解框架,其主程序包含以下4 部分:(1)對原始狀態(tài)變量方程求解;(2)對伴隨變量方程求解;(3)預估目標函數與約束值;(4)對靈敏度和偽密度進行濾波和投影.其中對原始狀態(tài)變量和伴隨變量控制方程的求解模塊均基于內置的標準可壓縮流動求解器rhoSimpleFoam 進行改造.對熱物性及熱力學偏導的求解模塊則按上節(jié)所述的建表-插值法編譯為動態(tài)鏈接庫.耦合了MMA (method of moving asymptotes)[32]算法作為優(yōu)化工具,同樣將其編譯成庫供主程序調用,對設計變量進行更新.求解器迭代優(yōu)化的終止條件以目標函數的相對變化小于0.1%或固定迭代步數為準則.

圖4 考慮變物性的拓撲優(yōu)化求解流程圖Fig.4 Flowchart of the topology optimization approach considering variable physical properties

2.4 共軛換熱數值驗證

本文所開展的拓撲優(yōu)化設計均基于二維模型,為了對其實際三維構型的流動換熱性能進行評估,將對拓撲優(yōu)化構型進行提取,并利用OpenFOAM 內置的標準共軛換熱求解器chtMultiRegionFoam 進行數值模擬.為了驗證該共軛換熱求解器的準確性,首先與Zhu 等[33]進行的流動換熱實驗進行對比分析.該實驗中固體材料為不銹鋼,燃料為正癸烷.按照對應的實驗條件,入口給定質量流量0.6083 g/s,入口溫度為625.93 K,出口給定背壓4.19 MPa,測試段體積熱源大小為1.758×108W/m3,連接段熱流邊界的熱流密度大小為11550 W/m2.采用k-? 湍流模型進行驗證,其結果如圖5 所示,顯然,該數值計算結果與實驗結果較為吻合,外壁溫及冷卻劑溫度的最大相對誤差均小于5%,證明該共軛換熱求解器可以較好地預測流動換熱性能.

圖5 冷卻劑及外壁面溫度的數值與實驗結果對比Fig.5 Comparison of predicted coolant and out wall temperatures against the experimental data

3 結果與分析

3.1 流熱耦合拓撲優(yōu)化設計

本節(jié)將針對發(fā)動機再生冷卻通道開展流熱耦合拓撲優(yōu)化設計,為了節(jié)省計算成本,以圖6 所示的二維換熱單元作為設計域,即不考慮構型在z方向的優(yōu)化.其幾何及邊界條件設置見圖6,計算域進出口尺寸d=5 mm,厚度為2 mm,進、出口段為長度10d的非設計域,中間設計域長度為36d,流體材料為正癸烷,固體材料為GH3128,其熱導率隨溫度的變化如下

圖6 流-熱耦合結構計算域示意圖Fig.6 Schematic of the computational domain of the coupled thermal–fluid structure

設定入口速度為1 m/s,初始溫度為300 K,對應的入口雷諾數為2500,計算域出口壓力設為3 MPa,其余邊界為對稱邊界,對設計域施加體積熱源Q=5×108W/m3.利用均勻結構化網格對區(qū)域進行離散化,考慮到要盡可能保留小尺度特征,同時還要保證加工可制造性,本文將網格尺寸設定為0.1 mm×0.1 mm,網格總數為140000,采用k-ε湍流模型進行模擬.

該優(yōu)化模型以P-norm 函數 Ψ 作為目標,施加材料體積約束gvol,Vmax取值為0.68,αmax取初始值為103,每步增長1.05 倍,最大值設為108.為了考察不同能量耗散約束值對拓撲優(yōu)化構型的影響,本文首先對 γ0=0.5,αmax=103的初始分布構型進行計算,得到其能量耗散值約為0.015 k g·m/s2,隨后將該數值作為參考值,對 φmax的取值范圍在0.02～0.25 內的算例進行了優(yōu)化,并從中提取出 φmax為0.03,0.05,0.14,0.19 和0.21 的5 個算例,分別命名為Case 1～Case 5.

圖7 為優(yōu)化得到的偽密度場分布及對應的速度分布與溫度分布.為了方便觀察,這里將計算結果按關于x軸對稱形式顯示.能夠發(fā)現,在材料比例基本一致的條件下,Case 1～Case 5 的固體域均發(fā)生了不同程度的分裂,且各分裂的固體胞元塊呈現出不同的排列方式.在Case 1 中,大部分固體域分布在通道兩側,流體流通面積較大,流速較低.而隨著能量耗散約束值增大,固體域分裂成的塊數逐漸增多,流體域內的微細通道數量增多,流體流通面積減小,流速增大,與固體域的接觸面積明顯增加.從Case 1～Case 5,流體的流動分離與再混合行為逐漸增多,多處區(qū)域存在局部的流動加速,同樣促進了流體與固體接觸面的熱量交換.

圖7 拓撲優(yōu)化構型及狀態(tài)變量云圖Fig.7 Contours of the optimized layouts and state variables

圖7 拓撲優(yōu)化構型及狀態(tài)變量云圖 (續(xù))Fig.7 Contours of the optimized layouts and state variables (continued)

通過觀察各偽密度法場對應的溫度云圖能夠發(fā)現,盡管采用了投影等數值策略,在優(yōu)化構型中仍存在部分灰色單元區(qū)域,如在Case 5 的后段,由于灰色單元的存在,導致該區(qū)域中出現流動死區(qū),其溫度呈現出非連續(xù)變化,這在一定程度上干擾了對流動換熱性能的評估.因此,為了評估真實的拓撲優(yōu)化通道的流動換熱性能,本文將進一步對拓撲優(yōu)化構型進行提取,將其對應的三維結構進行共軛換熱數值模擬.

3.2 拓撲冷卻通道的流動換熱性能

本節(jié)將對三維拓撲優(yōu)化通道開展共軛換熱數值模擬,并對其流動換熱性能進行分析.首先將5 種拓撲優(yōu)化構型中 γ <0.5 的區(qū)域進行提取,利用樣條曲線對拓撲胞元進行建模,圖8 所示為5 種拓撲通道(Case 1～Case 5)的三維模型,此外,加入了傳統(tǒng)直通道構型Case 0 作為對比算例,其固體壁肋寬約為1.70 mm,高度為2 mm,冷熱壁厚度均為0.75 mm.

圖8 拓撲優(yōu)化通道的三維幾何視圖Fig.8 Three-dimensional layouts of topology optimization channels

隨后對上述模型進行計算,其入口速度仍為1 m/s,初始溫度為300 K,出口壓力為3 MPa.按能量守恒原則,將體積熱源Q換算為在固體域底部面施加的面熱流q,其大小約為1 MW/m2,上下層固體域左右面為對稱邊界,其余壁面為絕熱邊界,采用kε湍流模型并結合壁面函數,對通道內各固體壁面添加邊界層網格,第一層網格高度為0.004 mm,以保證無量綱距離y+<10.

圖9 所示為計算得到的各通道中心截面和底部被加熱面的溫度分布云圖.能夠發(fā)現,由于固體胞元的不規(guī)則排列,5 個通道的熱壁面溫度呈現非均勻分布,在某些較大的胞元間隙位置,如Case 3 的190 mm 和210 mm 位置處,由于流體流通面積較大,流速較低,導致底面加熱位置的壁溫較高,但從Case 1到Case 5,各通道內的高溫區(qū)域顯著減少.

圖9 xy 截面的溫度分布云圖Fig.9 Contours of temperature on xy plane

圖10 為中心xy截面上的速度云圖、壓力云圖以及濕周底部面的努塞爾數云圖.在Case 1～Case 5中,冷卻劑平均流速逐漸增大,壓力損失也逐漸增大.在Case 3～Case 5 中,通道前端位置的固體胞元首先起到了節(jié)流的作用,使冷卻劑流速加快,隨后在流動過程中發(fā)生了大量的流動分離和再混合,帶來了較大的能量耗散和壓力損失.但由于冷卻劑流動邊界層和熱邊界層的分離和再附著,使得傳熱熱阻降低,增強了流體與固體間的換熱,局部努塞爾數增大.此外,在多個固體胞元的鈍前緣位置,局部努塞爾數也有明顯提高,這表明流體對固體的沖擊效應,同樣使得胞元前緣位置的傳熱增強.

表1 所示為各拓撲冷卻通道的流動和換熱性能數據,圖11 為熱壁面平均溫度和最高溫度與通道壓降的性能關系曲線.能夠發(fā)現,從Case 1～Case 5,各通道的壓降逐漸增加,且熱壁面平均溫度和最高溫度均隨壓降升高而顯著降低,其中Case 5 熱壁面平均溫度和最高溫度分別降低了91.1 K 和55.1 K.各通道的平均努塞爾數逐漸增加,表明各通道整體的換熱能力逐漸增強.但相比于Case 0,僅Case 3～Case 5 起到了強化換熱的效果,其平均努塞爾數增益百分比分別為12.6%,16.0%和23.4%.

表1 6 種通道的流動換熱性能Table 1 Flow and heat transfer performances of six channels

圖11 各通道熱壁面平均溫度和最高溫度隨壓降的性能關系曲線Fig.11 Relationship between the average and maximum temperature in terms of the pressure drop of the cooling channels

3.3 拓撲冷卻通道的強化換熱機理

圖12 所示為各通道在流向不同位置處yz截面上的渦量云圖.能夠觀察到,相比于Case 0～Case 2,Case 3～Case 5 通道內冷卻劑在多個位置呈現二次渦形態(tài),如在Case 4 和Case 5 的x=70～85 mm、Case 5 的105～115 mm 等位置.該渦系結構的形成是由于流體在流經固體胞元間隙時,受其他分支的流體摻混與擾動而形成,其渦系形態(tài)一直保持到胞元后緣,但沿流向渦強度逐漸減弱,直到下一個固體域分裂的區(qū)域,流體再次形成明顯的二次渦結構.該渦系結構促使近壁面流體與主流發(fā)生摻混,促進熱交換.

圖12 不同流向位置yz 截面上的渦量云圖Fig.12 Contours of vorticity on yz planes at different streamwise locations

圖13 所示為各拓撲冷卻通道在xy中心截面上冷卻劑的湍動能分布云圖.在Case 3～Case 5 通道內,多個固體胞元的側壁以及尾部回流區(qū)等位置,由于局部的流動加速或流動混合,使流體的湍動能得以激發(fā),換熱能力增強.此外,流體湍動能增強區(qū)域與二次渦形成區(qū)域基本吻合,表明二次渦結構有助于激發(fā)湍動能,進而增強換熱.

圖13 xy 截面湍動能云圖Fig.13 Contours of turbulence energy on xy plane

4 結論

本文針對發(fā)動機的再生冷卻結構,開展了考慮變物性的流熱耦合拓撲優(yōu)化設計,并對三維拓撲結構的流動換熱性能進行了數值模擬與分析,得到主要結論如下.

(1)在基于連續(xù)伴隨法的求解框架中,考慮工質的變物性會使伴隨方程的形式更加復雜,求解更加困難.利用建表-插值法對方程中熱物性、輸運性質及各偏導項進行求解,是一種有效且準確的方法.

(2)通過對流熱耦合結構進行拓撲優(yōu)化設計,發(fā)現其優(yōu)化構型隨著能量耗散約束值的增大而愈加復雜,且伴有多重的流動分離和再混合現象.

(3)通過對三維拓撲優(yōu)化構型進行共軛換熱數值模擬,結果表明,相比于直通道構型Case 0,拓撲優(yōu)化冷卻通道內的冷卻劑流態(tài)復雜,存在二次渦等渦系結構,在帶來更大壓力損失的同時,各通道局部的換熱效果得到增強.在5 個拓撲通道中(Case 1～Case 5),Case 3～Case 5 整體實現了強化換熱的效果,其平均努塞爾數增益百分比分別為12.6%,16.0%和23.4%.

綜上所述,面向發(fā)動機再生冷卻的拓撲優(yōu)化設計具有良好的應用前景,后續(xù)將繼續(xù)開展更多三維流熱耦合結構的設計工作,并逐步完善求解器的其他功能.

附錄A

如1.2 節(jié)所述,伴隨變量方程組滿足正文關系式(12),在對其進行化簡前,首先將目標函數的定義域分為邊界和內部域兩部分

則目標函數關于狀態(tài)變量的偏導數為

式中wd表示目標關于狀態(tài)變量的導數的方向,變量p,u和T對應的導數方向分別為pd,ud和Td. 下面對方程組(12)中的各積分項進行積分變換,結果如下

此時,拉格朗日函數關于狀態(tài)變量的偏導數可以推導為如下形式

若要滿足正文式(12),可得到伴隨變量控制方程

以及對應的邊界條件方程

對于方程(A12)～式(A14),本文將其分別作為伴隨變量ua,pa和Ta的顯式邊界條件,其數值隨優(yōu)化迭代而更新.方程中所涉及的目標函數關于狀態(tài)變量的各項偏導數表達式如表A1所示.此時化簡后的靈敏度表達式如下所示

對于能量耗散約束函數的靈敏度,其推導過程與換熱目標函數靈敏度相似,本文不予以贅述.

附表1 目標及約束函數關于狀態(tài)變量的偏導數Table 1 Derivatives of the objectives and constraints w.r.t state variables