混流式水輪機轉(zhuǎn)輪設(shè)計變量耦合強度分析

朱國俊，羅興锜，馮建軍，盧金玲，吳廣寬

?

混流式水輪機轉(zhuǎn)輪設(shè)計變量耦合強度分析

朱國俊，羅興锜※，馮建軍，盧金玲，吳廣寬

（西安理工大學(xué)水利水電學(xué)院，西安 710048）

在開展混流式水輪機轉(zhuǎn)輪的多學(xué)科優(yōu)化時通常需要進行設(shè)計變量的學(xué)科耦合強度分析，以便為建立簡潔、高效的多學(xué)科優(yōu)化求解策略提供依據(jù)。該文提出了一種基于全局相對靈敏度的轉(zhuǎn)輪設(shè)計變量耦合強度分析方法，該方法基于參數(shù)化全三維反問題設(shè)計理論實現(xiàn)對轉(zhuǎn)輪幾何的參數(shù)化控制，并通過改進后的Morris OAT法進行各學(xué)科目標函數(shù)對轉(zhuǎn)輪設(shè)計變量的全局相對靈敏度計算，然后以各設(shè)計變量的全局相對靈敏度集合為論域，引入模糊隸屬度函數(shù)量化計算設(shè)計變量對各學(xué)科目標函數(shù)的隸屬度。最后，以各設(shè)計變量的隸屬度值為基礎(chǔ)，提出了設(shè)計變量的耦合強度判定準則，為設(shè)計變量的耦合程度屬性確定提供了參考。采用該文所提出的分析方法對某混流式水輪機模型轉(zhuǎn)輪的設(shè)計變量開展了耦合強度分析，并根據(jù)變量的耦合強度分析結(jié)果開展了轉(zhuǎn)輪的優(yōu)化設(shè)計，優(yōu)化后的轉(zhuǎn)輪不僅使得水輪機在3個優(yōu)化工況下的水力效率分別提高0.2%、0.82%和1.2%，同時葉片的空化和強度性能也得到改善，該結(jié)果表明該文所提出的混流式水輪機轉(zhuǎn)輪設(shè)計變量耦合強度分析方法能可靠有效的界定各設(shè)計變量的耦合程度屬性，從而驗證了方法的可行性。

優(yōu)化；設(shè)計；計算機仿真；混流式轉(zhuǎn)輪；全局相對靈敏度；耦合強度分析；改進的Morris OAT 法

0 引 言

混流式水輪機轉(zhuǎn)輪的優(yōu)化設(shè)計方法由于受到傳統(tǒng)的串行設(shè)計開發(fā)模式的桎梏，設(shè)計人員通常先針對轉(zhuǎn)輪的流體力學(xué)性能開展優(yōu)化設(shè)計[1-6]，然后再進行轉(zhuǎn)輪的固體力學(xué)性能分析作為校核。這種方法曾經(jīng)取得一定的成果，但隨著如今國際水電市場的競爭環(huán)境越來越激烈，傳統(tǒng)方法已難以兼顧水輪機多層次的有時甚至是矛盾的設(shè)計需求，該方法的一些弊端（如開發(fā)時間長、不同學(xué)科的設(shè)計人員對任務(wù)目標的評價標準有差異而導(dǎo)致設(shè)計工作擱置等等）也在設(shè)計過程中暴露。為消除傳統(tǒng)方法中的弊端，轉(zhuǎn)輪多學(xué)科優(yōu)化設(shè)計方法逐漸獲得設(shè)計人員的關(guān)注并應(yīng)用于工程設(shè)計。

在混流式轉(zhuǎn)輪的多學(xué)科優(yōu)化設(shè)計過程中，多學(xué)科優(yōu)化求解策略的具體結(jié)構(gòu)決定了優(yōu)化問題求解的復(fù)雜性及難度，而轉(zhuǎn)輪設(shè)計變量的學(xué)科耦合強度屬性則是確定優(yōu)化求解策略具體結(jié)構(gòu)的關(guān)鍵。目前國內(nèi)外學(xué)者采用的分析設(shè)計變量耦合程度屬性的方法有決策分析法[7]、基于聚類思想的方法[8-9]以及基于靈敏度分析的方法[10-11]等。在這些方法中，類似于決策分析法的新方法由于應(yīng)用較少，其可靠性尚待研究；基于聚類思想的方法由于計算過程相對復(fù)雜，因此其應(yīng)用缺乏普適性；而基于靈敏度分析的方法由于出現(xiàn)時間早[12]且理論計算過程簡潔，所以獲得了相對廣泛的應(yīng)用。但是，混流式轉(zhuǎn)輪多學(xué)科優(yōu)化問題中各設(shè)計變量的物理意義不盡相同，導(dǎo)致了目標函數(shù)對各設(shè)計變量的靈敏度存在量綱及量級上的差異，不利于對比分析。此外，在建立混流式轉(zhuǎn)輪多學(xué)科優(yōu)化求解策略時，所需的設(shè)計變量耦合強度屬性的類別數(shù)通常少于設(shè)計變量數(shù)目，因此只依靠各設(shè)計變量的靈敏度屬性難以準確界定設(shè)計變量的耦合強度。

為解決基于靈敏度信息判定轉(zhuǎn)輪設(shè)計變量耦合強度屬性時出現(xiàn)的問題，本文提出了一種實用的設(shè)計變量耦合強度分析方法開展混流式轉(zhuǎn)輪多學(xué)科優(yōu)化設(shè)計中的變量耦合強度分析。以混流式轉(zhuǎn)輪的參數(shù)化全三維反問題設(shè)計方法為轉(zhuǎn)輪幾何控制策略，該方法首先通過改進原Morris OAT法中的抽樣技術(shù)提升其計算效率，然后基于改進后的Morris OAT法計算獲得各學(xué)科目標函數(shù)對轉(zhuǎn)輪設(shè)計變量的全局相對靈敏度信息，最后以設(shè)計變量的全局相對靈敏度數(shù)據(jù)集合為論域，引入模糊隸屬度函數(shù)計算獲得設(shè)計變量對各學(xué)科目標函數(shù)的隸屬度，并以此為基礎(chǔ)建立耦合強度判定準則為混流式轉(zhuǎn)輪設(shè)計變量的耦合強度屬性判定提供參考。

1 參數(shù)化全三維反問題設(shè)計方法

本文以混流式轉(zhuǎn)輪全三維有旋流動反問題設(shè)計方法為基礎(chǔ)，引入貝塞爾曲線參數(shù)化技術(shù)控制全三維反問題設(shè)計所需的邊界條件，形成了混流式轉(zhuǎn)輪的幾何控制策略?；炝魇睫D(zhuǎn)輪全三維反問題設(shè)計理論假設(shè)轉(zhuǎn)輪內(nèi)部的流動沿圓周方向呈周期性變化，因此其內(nèi)部流場可采用周向平均流場與周期性脈動流場的疊加進行模擬。周向平均流場通過由周向平均相對流動流函數(shù)和置于葉片中面的源匯所產(chǎn)生的勢函數(shù)構(gòu)建流動控制方程進行求解，周期性脈動流場則通過引入周期性鋸齒函數(shù)構(gòu)建周期性脈動流動控制方程進行求解，具體的方程形式參見文獻[13]。轉(zhuǎn)輪葉片方程根據(jù)葉片表面相對速度與葉片表面法向量垂直的物理邊界條件構(gòu)建，并在柱坐標系下最終化簡如式（1）。

1.1 葉片區(qū)速度矩的參數(shù)化控制

注：H和S分別為上冠、下環(huán)處速度矩控制點，=1,2,3,4。

Note:HandSare velocity momentcontrol point in hub and shroud,respectively,=1,2,3,4.

圖1 速度矩參數(shù)化控制示意圖

Fig.1 Schematic diagram of parameterized control of velocity moment

為了實現(xiàn)對葉片區(qū)速度矩分布的控制，如圖1中所示，采用四階三次貝塞爾曲線分別對上冠和下環(huán)處的速度矩分布實現(xiàn)參數(shù)化控制，上冠和下環(huán)間其余各軸面流線上的速度矩分布則根據(jù)線性分布插值獲得。轉(zhuǎn)輪葉片進口處速度矩只受導(dǎo)葉開度影響，所以在優(yōu)化過程中保持不變。此外，為了保證機組運行的穩(wěn)定性，維持轉(zhuǎn)輪出口速度矩從上冠至下環(huán)等于0，因此轉(zhuǎn)輪葉片區(qū)進口和出口處速度矩都固定不變，總結(jié)后即得速度矩控制點1、4、1、4不變。剩余可變化的控制點只有2、3、2、3，為了減少設(shè)計變量數(shù)目，只選擇它們的縱坐標y2、y3、y2、y3作為控制變量，最終，葉片區(qū)速度矩控制變量數(shù)量為4個。

1.2 葉片軸面流道的參數(shù)化控制

轉(zhuǎn)輪的軸面流道幾何也是全三維反問題設(shè)計所需的邊界條件，但是由于軸面流道中的下環(huán)幾何型線的改變會導(dǎo)致標稱直徑1的變化，同時會使下環(huán)出口錐角與尾水管錐管尺寸不匹配，因此本文只將上冠型線的幾何作為優(yōu)化設(shè)計的對象，通過控制轉(zhuǎn)輪上冠型線的形狀可以改善葉片的空化性能以及調(diào)節(jié)轉(zhuǎn)輪的過流能力。上冠型線是直接參與加工制造的關(guān)鍵幾何形狀之一，本文采用五階貝塞爾曲線對其開展參數(shù)化控制，如圖2所示。從圖2可知，1～5是轉(zhuǎn)輪上冠型線的貝塞爾控制頂點。為了保證轉(zhuǎn)輪上冠與活動導(dǎo)葉流道的光滑連接，優(yōu)化過程中保持控制點1坐標不變，同時只允許改變控制點2橫坐標x2。葉片出口過流斷面的面積是影響轉(zhuǎn)輪過流量和葉片空化性能的關(guān)鍵因素，因此控制點4的橫坐標x4也被選擇為優(yōu)化變量之一。除此之外，根據(jù)精簡優(yōu)化變量的原則，保持控制點3和5不變，則最終的上冠型線控制變量只有x2、x4。

注：P1~P5是轉(zhuǎn)輪上冠型線的貝塞爾控制頂點。

1.3 葉片厚度分布的參數(shù)化控制

轉(zhuǎn)輪葉片的厚度分布規(guī)律對葉片流體力學(xué)性能和強度性能都有直接影響，三維葉片的厚度分布規(guī)律主要由各流面上的翼型厚度分布規(guī)律組合形成，因此可通過控制各流面上的翼型厚度分布實現(xiàn)對其控制。流面上的翼型厚度分布規(guī)律可分解為最大厚度值和相對厚度分布規(guī)律兩部分，為確保最終三維葉片的正背面光滑，各流面翼型的最大厚度值從上冠到下環(huán)按線性進行分布[14]，同時參考相近水頭段優(yōu)秀轉(zhuǎn)輪葉片的最大厚度值分布，確定上冠與下環(huán)處翼型最大厚度的關(guān)系如式（3）。

注：T為控制點，=1,2,3,4。

Note:Tiscontrol point,=1,2,3,4.

圖3 翼型相對厚度控制示意圖

Fig.3 Schematic diagram of foil's relative thickness controlled

綜上所述，混流式轉(zhuǎn)輪葉片的參數(shù)化控制策略所需的設(shè)計變量參數(shù)集合如式（4）。

上述變量參數(shù)集合中的9個變量涵蓋了葉片區(qū)載荷控制參數(shù)、轉(zhuǎn)輪軸面流道控制參數(shù)以及葉片厚度控制參數(shù)，基本實現(xiàn)了對整個轉(zhuǎn)輪幾何的掌控。

2 目標函數(shù)選取及計算方法

轉(zhuǎn)輪優(yōu)化設(shè)計問題的主要目標是獲取一個具有良好能量、空化及強度性能的轉(zhuǎn)輪。為了考慮轉(zhuǎn)輪與上下游過流部件的匹配性能，轉(zhuǎn)輪的能量性能采用水輪機的水力效率作為衡量指標。轉(zhuǎn)輪的空化性能采用轉(zhuǎn)輪葉片上的最低壓力系數(shù)min作為衡量指標，而轉(zhuǎn)輪葉片的強度性能則采用葉片上的最大等效靜應(yīng)力max作為指標。因此，本文中混流式轉(zhuǎn)輪的多學(xué)科優(yōu)化目標函數(shù)集合為

式中水力效率和最低壓力系數(shù)min為無量綱數(shù)，max為最大等效靜應(yīng)力，MPa。

上述目標函數(shù)全部通過數(shù)值計算方法獲取。水力效率和葉片表面最低壓力系數(shù)通過采用CFX對水輪機進行全三維黏性數(shù)值計算獲得，葉片最大等效靜應(yīng)力則通過ANSYS有限元分析計算獲取。

3 基于改進Morris OAT法全局靈敏度分析

對設(shè)計變量參數(shù)集合中各設(shè)計變量開展耦合強度分析的前提是獲得各設(shè)計變量的全局靈敏度信息。常用的全局靈敏度計算方法有：Sobol法[19]、傅里葉靈敏度計算法[20]、多元回歸分析法[21]以及Morris OAT法[22-23]。Sobol法和傅里葉靈敏度計算法的計算公式復(fù)雜，在實際工程中應(yīng)用較少。多元回歸法和Morris OAT法都是基于全局抽樣的分析法，其分析計算過程簡便所以應(yīng)用更為廣泛，而兩者中Morris OAT法以其計算量少的優(yōu)勢又更勝一籌，所以本文采用其作為全局靈敏度計算（global sensitive analysis, GSA）的方法。

常規(guī)Morris OAT法中的抽樣技術(shù)是傳統(tǒng)拉丁超立方抽樣法，該方法存在所抽取的樣本點在設(shè)計空間內(nèi)分布不均勻的缺點，因此通常需要開展3～4次重復(fù)抽樣計算以確保抽樣點在設(shè)計空間內(nèi)分布的均勻性，極大的增加了全局靈敏度分析的計算量。為了改進該方法的缺陷，本文引入優(yōu)化拉丁超立方抽樣技術(shù)替代Morris OAT法中的抽樣方法。優(yōu)化拉丁超立方技術(shù)是以中心離差度量最小為準則開展設(shè)計空間樣本點抽取的試驗設(shè)計方法[24]，其抽取的樣本點在設(shè)計空間內(nèi)分布均勻且確保每個樣本的每個水平都被覆蓋。采用該抽樣技術(shù)與Morris OAT法相結(jié)合開展全局靈敏度計算時，只需進行一次樣本點抽取和靈敏度計算即可。通過抽樣技術(shù)獲取所有的設(shè)計樣本點后，根據(jù)以下方法開展各設(shè)計變量的全局相對靈敏度計算。

假如設(shè)計變量參數(shù)集合的維數(shù)為，全局靈敏度計算所需抽取的樣本點總數(shù)為，則對于個設(shè)計樣本點中的某個樣本點X，其表達式如下

X=(x1,x2,…,x,…,x) （6）

式中x代表設(shè)計變量參數(shù)集合中的第個設(shè)計變量x在樣本點X處的取值，=1, 2, …,，=1, 2, …,。

則根據(jù)Morris OAT法可知[25-26]：在樣本點X處，目標函數(shù)()對設(shè)計變量參數(shù)集合中的第個變量x的相對靈敏度為

式中(…)為目標函數(shù)值，為設(shè)計變量x在樣本點X處的擾動，=1, 2, …,，=1, 2, …,。樣本點X在全局靈敏度的計算過程中也被稱為基值點。根據(jù)文獻[27]，式（7）中各設(shè)計變量的擾動可取為5%，然后可通過式（7）計算獲得x在所有抽樣點處的相對靈敏度。最后，根據(jù)Morris OAT法，依據(jù)式（8）計算x的全局相對靈敏度。

4 耦合強弱程度判定準則的建立

設(shè)計變量的耦合強弱程度由其對各學(xué)科目標函數(shù)的影響程度大小體現(xiàn)。同時對各個學(xué)科的目標函數(shù)影響程度大的變量即為強耦合變量，只對其中某個學(xué)科的目標函數(shù)有影響，而對其他學(xué)科目標函數(shù)影響程度很小或沒有影響的變量為弱耦合變量，由此可見，設(shè)計變量耦合程度的“強”和“弱”是一個模糊概念。為了對其進行定量劃分，本文引入模糊數(shù)學(xué)中隸屬度函數(shù)的概念，以各設(shè)計變量的全局相對靈敏度數(shù)值為基礎(chǔ)，通過歸一化計算獲得設(shè)計變量對目標函數(shù)的隸屬度，如此即可在統(tǒng)一的量化標準下制定設(shè)計變量的強弱耦合程度判定準則。

在本文所關(guān)注的3個目標函數(shù)中，水力效率和最低壓力系數(shù)min為流體力學(xué)學(xué)科的目標函數(shù)，而葉片的最大等效靜應(yīng)力max則為固體力學(xué)學(xué)科的目標函數(shù)。各目標函數(shù)對設(shè)計變量的全局相對靈敏度可寫成式（9）所示的集合模式。

將各目標函數(shù)的全局相對靈敏度集合U、U、U分別作為論域，定義“對目標函數(shù)影響程度大的設(shè)計變量”為論域上的模糊集，則通過該模糊集的隸屬度函數(shù)即可完成論域中各元素對模糊集的隸屬度確定，從而實現(xiàn)論域到隸屬度集合的變換。

不同的研究人員會根據(jù)自身的習慣和偏好建立或選擇不同型式的隸屬度函數(shù)，盡管隸屬度函數(shù)在型式上有區(qū)別，但在反映、解決實際模糊問題時殊途同歸。本文中的模糊集定義是“對目標函數(shù)影響程度大的設(shè)計變量”，屬于尋找論域中“偏大型”元素的模糊問題，因此本文選擇了半梯形分布的偏大型函數(shù)()作為隸屬度函數(shù)。

式中取所研究論域中的最小值，取所研究論域中的最大值。通過式（9）和（10），設(shè)計變量對各目標函數(shù)的隸屬度被統(tǒng)一映射到[0,1]的數(shù)值區(qū)間，并可寫成如下形式。

A={a1,a2,…,a9}

A={a1,a2,…,a9} （11）

A={a1,a2,…,a9}

式中a、a、a分別表示設(shè)計變量對水力效率、最低壓力系數(shù)min以及葉片上最大等效靜應(yīng)力max的隸屬度值，=1, 2, …, 9。

在獲得每個設(shè)計變量對各目標函數(shù)的隸屬度以后，為了使設(shè)計變量的耦合強度分類更為合理，本文根據(jù)設(shè)計變量的隸屬度數(shù)值將設(shè)計變量的耦合強度分為強耦合、中等耦合和弱耦合3種。對于設(shè)計變量參數(shù)集中的任意一個設(shè)計變量，本文將其耦合強度屬性分類的原則定義如下：

1）若設(shè)計變量對各個學(xué)科目標函數(shù)的隸屬度數(shù)值都為1，則該設(shè)計變量為強耦合設(shè)計變量。

2）若設(shè)計變量對各個學(xué)科目標函數(shù)的隸屬度函數(shù)值都屬于（0,1）的數(shù)值區(qū)間，或只對其中某個學(xué)科目標函數(shù)的隸屬度函數(shù)值為1而對其他學(xué)科目標函數(shù)的隸屬度函數(shù)值屬于（0,1）的數(shù)值區(qū)間，則稱該設(shè)計變量為中等耦合設(shè)計變量。

3）若設(shè)計變量只對某個學(xué)科目標函數(shù)的隸屬度函數(shù)值不為0而對其他學(xué)科目標函數(shù)的隸屬度函數(shù)值都為0，則表明該設(shè)計變量為弱耦合設(shè)計變量。

根據(jù)上述設(shè)計變量的耦合強度屬性定義，對本文所研究的混流式轉(zhuǎn)輪葉片的任一設(shè)計變量，可制定其耦合強度判定的隸屬度取值原則如表1所示。

表1 設(shè)計變量耦合強度判定的隸屬度取值原則

注：a、a、a分別代表任一設(shè)計變量對目標函數(shù)、min及max的隸屬度。

Note:,aandarepresent the membership of the design variables to the objective functions,minandmaxrespectively.

針對混流式轉(zhuǎn)輪葉片，綜合式（11）中各目標函數(shù)隸屬度集合中的數(shù)據(jù)及表1中的判定原則即可判定各設(shè)計變量的耦合強度屬性，然后再根據(jù)設(shè)計變量集合中的強、弱耦合變量的數(shù)量即可確定多學(xué)科優(yōu)化設(shè)計求解策略，為混流式轉(zhuǎn)輪的多學(xué)科優(yōu)化奠定基礎(chǔ)。

5 應(yīng)用算例

5.1 算例相關(guān)參數(shù)

表2 設(shè)計變量初始值及變化范圍

注：當變量下標中含有.ori時，表示其為對應(yīng)變量的初始值。

Note: The variable which subscript contains .ori represents the initial value of corresponding variable.

確定設(shè)計變量的變化范圍后，選取每個設(shè)計變量的變化水平數(shù)為10，并使水平取值在各變量的變化范圍內(nèi)均布，然后將每個設(shè)計變量的10個水平取值分別以1～10的整數(shù)進行編碼，則根據(jù)優(yōu)化拉丁超立方抽樣方法可以產(chǎn)生10個基值點的水平取值矩陣，在確定了基值點以后即可按照本文第3節(jié)所述方法開展各設(shè)計變量的全局相對靈敏度計算。

5.2 目標函數(shù)的計算方法

在相對靈敏度計算過程中，各基值點對應(yīng)的目標函數(shù)值以及設(shè)計變量發(fā)生擾動后的目標函數(shù)值全部通過數(shù)值計算方法獲取，根據(jù)Morris OAT法的計算原理，總共需要開展數(shù)值分析的次數(shù)為：基值點數(shù)×設(shè)計變量數(shù)=90。

水力效率和最低壓力系數(shù)min為流體力學(xué)學(xué)科的目標函數(shù)，本文擬通過對水輪機內(nèi)部流動進行流體動力學(xué)分析獲得其數(shù)值。盡管目前的計算機運算速度已經(jīng)獲得了巨大提升，但是開展90次混流式水輪機全流道計算仍然具有很大的計算工作量，因此，出于減少靈敏度分析過程中計算工作量的目的，本文選擇了單周期的活動導(dǎo)葉、轉(zhuǎn)輪與尾水管作為流體動力學(xué)計算域開展數(shù)值分析，采用六面體結(jié)構(gòu)化網(wǎng)格離散的計算域如圖4所示。

圖4 流體力學(xué)計算域

由于目標函數(shù)對設(shè)計變量的敏感度實質(zhì)是目標函數(shù)對轉(zhuǎn)輪幾何形狀的敏感度，而轉(zhuǎn)輪幾何變化對水輪機性能的影響主要體現(xiàn)在2點：第一，轉(zhuǎn)輪幾何變化對其自身性能的改變；第二，轉(zhuǎn)輪幾何變化對下游尾水管水力性能的影響，因此，選擇圖4所示計算域進行流體動力學(xué)分析足以反映目標函數(shù)對轉(zhuǎn)輪設(shè)計變量的敏感度。圖4所示計算域的進口為活動導(dǎo)葉進口，給定質(zhì)量流量邊界條件，方向為固定導(dǎo)葉出流方向。計算域出口為尾水管出口，給定靜壓邊界條件。所有固壁面都采用無滑移壁面。計算域網(wǎng)格總數(shù)經(jīng)過以水力效率為基準的無關(guān)性驗證后確定為150萬。湍流模型采用SST湍流模型，該模型結(jié)合Wilcox模型和模型的優(yōu)點，并且對于逆壓梯度引起的分離流動有較高的預(yù)測精度[28-29]。為了檢驗數(shù)值計算的可靠性，在最優(yōu)單位轉(zhuǎn)速11= 67.5 r/min時，選取了單位流量510、549、588 L/s 3個工況對原始轉(zhuǎn)輪下該模型水輪機的性能進行了數(shù)值計算，將計算所得水力效率與試驗值進行對比如圖5所示，圖5中的偏差為絕對偏差，其計算公式如式（12）所示。

式中ηCFD為數(shù)值計算所得水力效率，ηEXP為模型試驗的數(shù)據(jù)值。

從圖5的數(shù)據(jù)對比中發(fā)現(xiàn)，盡管由于沒考慮蝸殼、固定導(dǎo)葉的損失使得計算值偏高于試驗值，但數(shù)值偏差小于2%，同時數(shù)值計算較好地反映了試驗數(shù)據(jù)的趨勢，因此采用本文數(shù)值計算方法分析轉(zhuǎn)輪形狀變化對水輪機流體力學(xué)性能目標函數(shù)的影響是可行的。另一個流體力學(xué)目標函數(shù)min則根據(jù)CFD計算結(jié)果以式（13）進行計算獲得其數(shù)值。

式中min為葉片上壓力最低值，Pa；為流體密度，kg/m3；U1為轉(zhuǎn)輪葉片進口標稱直徑1處的圓周速度，m/s。

固體力學(xué)目標函數(shù)max通過對葉片進行有限元分析獲取。依據(jù)文獻[30]中的研究結(jié)論，為了縮短整個靈敏度計算的時間周期及有限元計算工作量，可只選取葉片作為分析模型研究葉片的強度，因此本文只選取葉片實體域作為有限元分析對象。葉片實體域的計算載荷條件設(shè)置如下圖6a中所示。

葉片首先需要考慮重力載荷，其次，葉片還受到旋轉(zhuǎn)所引發(fā)的離心力作用，因此還需考慮旋轉(zhuǎn)效應(yīng)。葉片表面設(shè)置為流固耦合面，根據(jù)葉片的受力特點，葉片實體與上冠、下環(huán)的接觸端設(shè)為固定約束條件。葉片實體采用solid185六面體單元進行離散，劃分的葉片實體網(wǎng)格如圖6b中所示，葉片實體網(wǎng)格節(jié)點總數(shù)為238 915。有限元計算中流固耦合面上各節(jié)點的水壓力載荷按其對應(yīng)的坐標位置在CFD計算結(jié)果中插值獲得。計算完成后提取葉片上的最大等效應(yīng)力作為固體力學(xué)的目標函數(shù)max。

圖6 載荷條件與有限元網(wǎng)格

5.3 耦合強度分析結(jié)果及有效性驗證

在設(shè)計點工況下，根據(jù)上述目標函數(shù)計算方法對本文5.1節(jié)中的10個基值點以及各基值點處設(shè)計變量發(fā)生擾動后的目標函數(shù)值進行了計算，然后根據(jù)式（7）、（8）計算獲得各目標函數(shù)對各設(shè)計變量的全局相對靈敏度如表3所示。

表3 設(shè)計變量的全局相對靈敏度和隸屬度

圖7 多學(xué)科可行性求解策略結(jié)構(gòu)

根據(jù)圖7的優(yōu)化求解策略建立具體的優(yōu)化設(shè)計框架結(jié)構(gòu)開展轉(zhuǎn)輪的優(yōu)化，優(yōu)化工況點如表4。

表4 優(yōu)化工況點

每個工況包含、min和max3個目標函數(shù)，目標函數(shù)高達9個，使優(yōu)化問題成為高維優(yōu)化問題。為有效求解該優(yōu)化問題，采用作者提出的嵌入設(shè)計者偏好的多工況優(yōu)化設(shè)計方法對高維目標函數(shù)集合進行了降維處理。該方法通過引入設(shè)計者對各目標函數(shù)的期望將優(yōu)化目標轉(zhuǎn)變?yōu)椋涸跐M足所有約束條件下，尋求設(shè)計空間內(nèi)的一個解，使其在優(yōu)化工況下的性能與期望性能間的歐幾里德距離最小。該方法的可行性和有效性已經(jīng)獲得驗證[31]。采用該方法轉(zhuǎn)換后的新目標函數(shù)如式（14）～（16）所示，目標函數(shù)由9個縮減為3個，優(yōu)化問題的求解難度極大降低。

式中d、d和d為新的目標函數(shù)；η、min, j和max, j分別為工況下的3個目標函數(shù)值；η,j、min,Des,j和max,Des,j則為工況下3個目標函數(shù)的期望值，=1,2,3。

限于篇幅，下文只給出轉(zhuǎn)輪的優(yōu)化結(jié)果以驗證本文方法的有效性。優(yōu)化采用的NSGA-II算法的種群數(shù)為48，進化代數(shù)為35，交叉概率為0.8。優(yōu)化求解在服務(wù)器上完成，優(yōu)化過程中的種群在各二維目標函數(shù)平面內(nèi)的投影如圖8所示。

從計算所得的Pareto解中選取了性能提升較為均衡的解作為優(yōu)化解，優(yōu)化前后的幾何對比如圖9所示。

圖8 優(yōu)化中的種群分布

圖9 葉片優(yōu)化前后的幾何對比

采用數(shù)值計算方法對優(yōu)化前后的轉(zhuǎn)輪開展全流道性能對比，網(wǎng)格采用經(jīng)無關(guān)性驗證后總數(shù)為654萬的六面體結(jié)構(gòu)化網(wǎng)格，湍流模型采用SST，蝸殼進口邊界條件給定質(zhì)量流量，尾水管出口邊界條件給定靜壓，計算后的性能對比如圖10所示。

圖10a中給出了轉(zhuǎn)輪優(yōu)化前后的水輪機效率、轉(zhuǎn)輪葉片上的最低壓力系數(shù)的對比，圖10b則給出了轉(zhuǎn)輪葉片強度性能的對比。綜合上述對比結(jié)果可知，對于優(yōu)化的3個工況，機組的水力效率分別提升了0.2%、0.82%和1.2%，同時，轉(zhuǎn)輪葉片的空化性能和強度也獲得提升，說明轉(zhuǎn)輪的流體力學(xué)和固體力學(xué)性能獲得了較為全面的改善。優(yōu)化結(jié)果證明了根據(jù)耦合強度分析結(jié)果選擇的優(yōu)化求解策略的可行性和有效性，也驗證了本文所提出的耦合強度分析方法是有效的。

圖10 優(yōu)化前后的水輪機性能對比

6 結(jié) 論

1）本文提出了綜合最優(yōu)拉丁超立方抽樣技術(shù)與Morris OAT法的全局相對靈敏度分析方法，通過該方法在混流式轉(zhuǎn)輪設(shè)計變量靈敏度分析中的應(yīng)用表明，其只需進行一次基值點的抽樣計算。與傳統(tǒng)Morris OAT法至少需要進行2～3次抽樣計算相比，計算工作量減少了50%以上。

2）建立了一種判定混流式轉(zhuǎn)輪設(shè)計變量耦合強度的方法，并應(yīng)用該方法對某混流式模型轉(zhuǎn)輪的設(shè)計變量開展了耦合強度分析。結(jié)果表明，在該混流式模型轉(zhuǎn)輪的優(yōu)化設(shè)計變量中，除了最大厚度值為弱耦合變量以外，其余變量均為中等耦合變量。同時，分析結(jié)果也表明本文所提出的耦合強度分析方法具備可行性。

3）根據(jù)耦合強度分析結(jié)果選擇了綜合計算成本低的多學(xué)科優(yōu)化策略開展某混流式模型轉(zhuǎn)輪的優(yōu)化，最終結(jié)果表明優(yōu)化設(shè)計不僅使得水輪機的水力效率在3個優(yōu)化工況下分別提高了0.2%、0.82%和1.2%，還改善了轉(zhuǎn)輪葉片的空化、強度性能，證明了所選擇的優(yōu)化求解策略的可行性，也驗證了本文所提出的耦合強度分析方法是有效的。

[1] Vishnu Prasad, Gahlot V K, Krishnamachar P. CFD approach for design optimization and validation for axial flow hydraulic turbine[J]. Indian Journal of Engineering and Materials Sciences, 2009, 16(4): 229－236.

[2] Derakhshan S H, Mostafavi A. Optimization of GAMM francis turbine runner[J]. World Academy of Science- Engineering and Technology, 2011(59): 717－723.

[3] Priyono Sutikno. Numerical simulation of the Francis Turbine and CAD used to Optimized the Runner Design[C]// The 10thAsian International Conference on Fluid Machinery, Malaysia, 2010: 816－825.

[4] Suthep Kaewnai, Somchai Wongwises. Improvement of the runner design of francis turbine using computational fluid dynamics[J]. American Journal of Engineering and Applied Science, 2011, 4(4): 540－547.

[5] 劉敏，周大慶，吳國穎，等. 超低水頭軸流式水輪機CFD優(yōu)化及流動特性研究[J]. 水電能源科學(xué)，2016(1)：150－153.

Liu Min, Zhou Daqing, Wu Guoying, et al. CFD optimization and study of flow characteristics of axial flow turbine with ultra-low head[J]. Water Resources and Power, 2016(1): 150－153. (in Chinese with English abstract)

[6] 孫見波，李穎，徐偉，等. 基于CFD和FEM的混流式水輪機數(shù)值模擬研究[J]. 電網(wǎng)與清潔能源，2016，32(9)：146－150.

Sun Jianbo，Li Ying，Xu Wei，et al. Research on numerical simulation of the Francis turbine based on CFD and FEM[J]. Power System and Clean Energy, 2016, 32(9): 146－150. (in Chinese with English abstract)

[7] 曹鵬彬，肖人彬，庫瓊. 公理設(shè)計過程中耦合設(shè)計問題的結(jié)構(gòu)化分析方法[J]. 機械工程學(xué)報，2006，42(3)：46－55.

Cao Pengbin, Xiao Renbin, Ku Qiong. Structural analytical approach to coupled design in design with axiomatic design[J]. Chinese Journal of Mechanical Engineering, 2006, 42(3): 46－55. (in Chinese with English abstract)

[8] Yao Tsung Ko. A dynamic planning method for new product development management[J]. Journal of the Chinese Institute of Industrial Engineers, 2010, 27(2): 103－120.

[9] 彭翔，劉振宇，譚建榮，等. 基于多重耦合聚類的復(fù)雜產(chǎn)品多變量關(guān)聯(lián)設(shè)計模型分解[J]. 機械工程學(xué)報，2013，49(3)：111－121.

Peng Xiang, Liu Zhenyu, Tan Jianrong, et al. Multivariable correlative model decomposition for complex product design based on clustering with multiple couplings[J]. Journal of Mechanical Engineering, 2013, 49(3): 111－121. (in Chinese with English abstract)

[10] 陳柏鴻，肖人彬，劉繼紅，等. 復(fù)雜產(chǎn)品協(xié)同優(yōu)化設(shè)計中耦合因素的研究[J]. 機械工程學(xué)報，2001，37(1)：19－23.

Chen Baihong, Xiao Renbin, Liu Jihong, et al. Research on the coupled factors in MDO for complex products[J]. Chinese Journal of Mechanical Engineering, 2001, 37(1): 19－23. (in Chinese with English abstract)

[11] 袁文強，劉玉生，葉曉平，等. 結(jié)合拓撲解耦的多學(xué)科優(yōu)化方法[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報，2015，27(9)：1795－1804.

Yuan Wenqiang, Liu Yusheng, Ye Xiaoping, et al. A multidisciplinary approach combining with topological decoupling[J]. Journal of Computer-Aided Design & Computer Graphics, 2015, 27(9): 1795－1804. (in Chinese with English abstract)

[12] Sobieszczanski-Sobieski J. A linear decomposition method for optimization problems-blueprint for development[R]. NASA Technical Memorandum 83248, 1982.

[13] 羅興锜. 水力機械現(xiàn)代設(shè)計理論及應(yīng)用[M]. 西安：西安交通大學(xué)出版社，1997.

[14] 高建銘，姚志民. 水輪機的水力計算[M]. 北京：電力工業(yè)出版社，1982.

[15] 肖若富，王正偉，羅永要. 基于流固耦合的混流式水輪機轉(zhuǎn)輪靜應(yīng)力特性分析[J]. 水力發(fā)電學(xué)報，2007，26(3)：120－123.

Xiao Ruofu, Wang Zhengwei, Luo Yongyao. Stress analysis of Francis turbine runners based on FSI[J]. Journal of Hydroelectric Engineering, 2007, 26(3): 120－123. (in Chinese with English abstract)

[16] 闞闞，鄭源，趙連輝，等. 基于流固耦合的混流式水輪機轉(zhuǎn)輪強度分析[J]. 水電站機電技術(shù)，2015(2)：9－11.

Kan Kan, Zheng Yuan, Zhao Lianhui, et al. Strength analysis of Francis turbine runner based on fluid structure interaction[J]. Mechanical & Electrical Technique of Hydropower Station, 2015(2): 9－11. (in Chinese with English abstract)

[17] 崔軍玲，余波，廖翔，等. 基于流固耦合的低水頭混流式水輪機轉(zhuǎn)輪葉片靜力分析[J]. 水力發(fā)電，2017，3(3)：77－80.

Cui Junling, Yu Bo, Liao Xiang, et al. Static analysis of low head Francis turbine runner blades based on fluid-structure interaction[J]. Water Power, 2017, 3(3): 77－80. (in Chinese with English abstract)

[18] 吳廣寬，羅興锜，馮建軍，等. 基于瞬態(tài)流固耦合的混流式轉(zhuǎn)輪葉片裂紋成因分析[J]. 農(nóng)業(yè)工程學(xué)報，2015，31(8)：92－98.

Wu Guangkuan, Luo Xingqi, Feng Jianjun, et al. Cracking reason for Francis turbine blades based on transient fluid structure interaction[J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of CSAE), 2015, 31(8): 92－98. (in Chinese with English abstract)

[19] Sobol I M. Sensitivity estimates for nonlinear mathematical models[J]. Mathematical Modeling Computational Experiments, 1993, 1 (4): 407－414.

[20] Cukier R I, Fortuin C M, Shuler K E, et al. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients: I. Theory[J]. The Journal of Chemical Physics, 1973, 59(3): 3873－3878.

[21] 馮海寬，李振海，金秀良，等. 基于EFAST和PLS的蘋果葉片等效水厚度高光譜估算[J]. 農(nóng)業(yè)工程學(xué)報，2016，32(12)：165－171.

Feng Haikuan, Li Zhenhai, Jin Xiuliang, et al. Estimating equivalent water thickness of apple leaves using hypersecptral data based on EFAST and PLS[J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of CSAE), 2016, 32(12): 165－171. (in Chinese with English abstract)

[22] Loubière P, Jourdan A, Siarry P, et al. A sensitivity analysis method for driving the artificial bee colony algorithm’s search process[J]. Applied Soft Computing, 2016, 41: 515－531.

[23] Senthamaraikkannan G, Gates I, Prasad V. Development of a multiscale microbial kinetics coupled gas transport model for the simulation of biogenic coal bed methane production[J]. Fuel, 2016, 167: 188－198.

[24] Jin Ruichen, Chen Wei, Agus Sudjianto. An efficient algorithm for constructing optimal design of computer experiments[J]. Journal of Statistical Planning and Inference, 2005, 134(1): 268－287.

[25] van Griensven A, Meixner T, Grunwald S, et al. A global sensitivity analysis tool for the parameters of multi-variable catchment models[J]. Journal of Hydrology, 2006, 324(1): 10－23.

[26] 朱嘉祺，徐向陽，何爽. 基于 LH-OAT的SWMM 模型參數(shù)敏感性分析[J]. 中國農(nóng)村水利水電，2014(3)：84－87.

Zhu Jiaqi, Xu Xiangyang, He Shuang. An analysis of parameter sensitivity of SWMM model based on LH-OAT method[J]. China Rural Water and Hydropower, 2014(3): 84－87. (in Chinese with English abstract)

[27] 李軍，梅琨，張明華. 水質(zhì)模型的全局參數(shù)敏感性分析[J]. 浙江農(nóng)業(yè)學(xué)報，2012，24(2)：314－320.

Li Jun, Mei Kun, Zhang Minghua. Global parameter sensitivity analysis on water quality model[J]. Acta Agriculturae Zhejiangensis, 2012, 24(2): 314－320. (in Chinese with English abstract)

[28] Menter F R. Two-equation eddy viscosity turbulence models for engineering applications[J]. AIAA Journal, 1994, 32(8): 1598－1605.

[29] Menter F R, Kuntz M, Langtry R. Ten years of industrial experience with the SST turbulence model[C]// 4th Int. Symp. Turbulence, Heat Mass Transfer, Antalya, Turkey, 2003.

[30] 張麗霞. 混流式水輪機轉(zhuǎn)輪葉片疲勞裂紋控制研究[D]. 北京：清華大學(xué)，2010.

Zhang Lixia. Research on Fatigue Crack Control for Francis Turbine Blade[D]. Beijing: Tsinghua University. (in Chinese with English abstract)

[31] Luo Xingqi, Zhu Guojun, Feng Jianjun. Multi-point design optimization of hydrofoil for marine current turbine[J]. Journal of Hydrodynamics, Ser. B, 2014, 26(5): 807－817.

朱國俊，羅興锜，馮建軍，盧金玲，吳廣寬.混流式水輪機轉(zhuǎn)輪設(shè)計變量耦合強度分析[J]. 農(nóng)業(yè)工程學(xué)報，2017，33(22)：36－45. doi：10.11975/j.issn.1002-6819.2017.22.005 http://www.tcsae.org

Zhu Guojun, Luo Xingqi, Feng Jianjun, Lu Jinling, Wu Guangkuan. Variable coupling strength analysis of Francis turbine runner design[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(22): 36－45. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2017.20.005 http://www.tcsae.org

Variable coupling strength analysis of Francis turbine runner design

Zhu Guojun, Luo Xingqi※, Feng Jianjun, Lu Jinling, Wu Guangkuan

(,,710048,)

Hydropower is one of the most reliable source of green energy and is widely used to meet real-time electricity demand. The ever escalating energy demand and iterative development in hydropower have pushed the generation of hydropower towards a low lost-cost and more flexible state which can work in a wider hydrodynamic region without compromising its performance. So the research on multidisciplinary optimization of hydro turbine needs to be carried out. Francis-type turbine is attractive because the efficiency of the unit at its design point exceeds the other types. As mentioned above, with the needs of robust operation in wide region without compromising its performance, the Francis turbine runner has to be optimized under multidisciplinary constraints. Multidisciplinary optimization technology of Francis runners has always been a key point in the field of turbine research. The coupling strength of each variable is the basis for establishing a concise and efficient multidisciplinary optimization strategy. So in the multidisciplinary optimization of Francis runner, it is necessary to analyze the coupling strength of design variables. In this paper, based on the global relative sensitivity of design variables, a coupling strength analysis method of runner design variables has been presented. First, combing three-dimensional (3D) inverse design theory with curve parametric method, a parametric 3D inverse-design method has been developed to control the Francis runner geometry. According to this parametric 3D inverse-design method, the meridional flow passage and geometry of runner blade can be controlled by a set of discrete parameters. These parameters are treated as runner design variables. The hydraulic efficiency of turbine, the minimal pressure coefficient on runner blade and the maximum static stress of blade are selected as objective functions. Then the global relative sensitivity of each objective function to design variables is calculated by an improved Morris OAT method. Compared with the traditional Morris OAT method, an optimal Latin hypercube design technique is adopted in the improved method. The optimal Latin hypercube design technique has a better space filling property than the traditional one. So the sampling only needs to be done once, which reduces the time of global relative sensitivity analysis without compromising its precision. On the basis of the improved Morris OAT method, the speed of global relative sensitivity calculation is improved obviously. Secondly, the set which contains global relative sensitivity of design variables is treated as universe, and the semi-trapezoidal function in fuzzy mathematics is imported to quantify the membership of each variable to objective functions. After obtaining the membership of each variable, a judgment principle of variable’s coupling strength is presented in this paper. According to this judgment principle, the coupling strength property of each design variable can be defined as strong coupling, middle coupling and weak coupling. Finally, based on the coupling strength of each design variable, the multidisciplinary optimization strategy of Francis runner is established. And then the multidisciplinary optimization of Francis runner can be carried out. In order to validate the coupling strength analysis method, an application case is presented in the paper. In that case, the design variables of a Francis turbine runner with the head of 200 m are chosen as the analysis object. From the analysis result, it can be found that only the coupling strength property of maximum thickness is weak coupling, and the others are middle coupling. According to the result, a multidisciplinary optimization strategy is established and used in the optimization of Francis turbine runner. After optimization, both the hydraulic performance and strength property of the optimized runner are improved, demonstrating the effectiveness of the coupling strength analysis method. In conclusion, the method proposed in this paper can be adopted to analyze the coupling strength of design variables of Francis runner.

optimization; design; computer simulation; Francis runner; global relative sensitivity; coupling strength analyze; improved Morris OAT method

10.11975/j.issn.1002-6819.2017.22.005

TK733+.3

A

1002-6819(2017)-22-0036-10

2017-05-16

2017-08-11

國家自然科學(xué)基金項目（51379174、51679195）；國家自然科學(xué)基金重點項目（51339005）

朱國俊，博士，講師，研究方向為流體機械優(yōu)化設(shè)計理論、流動分析等。Email：guojun_zhu1984@126.com

羅興锜，教授，研究方向為流體機械流動理論及優(yōu)化設(shè)計、流體機械振動與穩(wěn)定性。Email：luoxq@xaut.edu.cn