李 偉，孫鳳芝，王金銘

(1.沈陽工業(yè)大學 理學院，遼寧 沈陽110870;2.大慶師范學院 數(shù)學科學學院，黑龍江 大慶163712)

0 引言

目前，國內(nèi)對風力發(fā)電機組本身的氣動特性分析研究主要體現(xiàn)在兩個方面:一是葉片優(yōu)化設計［1］;二是風力發(fā)電機系統(tǒng)控制、能量轉(zhuǎn)換過程［2-3］。而對風力機氣動流場模擬［4-9］主要是先對風輪機數(shù)學模型做一定程度的簡化，再采用間接耦合方法［10-11］進行數(shù)值求解。由于間接耦合方法存在收斂性難以保證等缺點，因此采用直接耦合方法［12-14］數(shù)值求解風力機流場問題，具有重要的理論意義和應用價值。

本文的研究從葉片動態(tài)旋轉(zhuǎn)入手，將風力機流場的求解區(qū)域看成圓柱體型區(qū)域，建立圓柱坐標系下的數(shù)學模型;將有限差分方法與牛頓法的結(jié)合應用到風力機流場數(shù)值模擬中，從理論上分析牛頓迭代格式局部收斂性，給出局部收斂性條件。

1 數(shù)學模型

1.1 基本控制方程

風力機流場是一個可壓縮、有黏性的非定常流場，其控制方程屬于流體力學方程。由于求解區(qū)域為圓柱體型區(qū)域，因此用圓柱坐標系(x，r，θ)較為方便。在圓柱坐標系下的控制方程可寫為:

式中:u、v、w 分別為風速在x、r、θ 方向上的分量，ρ 為空氣密度，R 為氣體常數(shù)，T 是溫度，s1= λ +2μ =4μ/3、s2= λ + μ = μ/3，λ 為流體的第二粘性系數(shù)，λ = -2μ/3，μ 為流體的動力粘度，fx= 0，fr= mg，fθ=3π/2。

1.2 定解條件、離散化格式及牛頓迭代格式

如圖1所示的求解區(qū)域Ω 包括入口區(qū)域和出口區(qū)域兩個部分，其中:邊界Γ1為區(qū)域Ω 的外表面部分，Γ2為Ω1與Ω2的交界面。

圖1 求解區(qū)域

②邊界條件:在Γ1上邊界條件:v = v0，ρ = ρ0，在Γ2| Γ3上滿足控制方程(1)～(4)，在Γ3滿足邊界條件:

其中:v0為無窮遠來流風速，ρ0為初始的空氣密度，A 為風輪掃掠面積，a 為軸向誘導因子，ω 為角速度矢量的大小?？傻糜煞匠?1)～(4)、初始條件、邊界條件組成的偏微分方程組的定解問題。

將如圖1所示的求解區(qū)域Ω 進行如下剖分:

其中:hx=取時間步長為τ，時刻tm= mτ(m = 1，2，…)。

將方程(1)在時刻tm、在空間節(jié)點(xi，rj，θk)處對時間導數(shù)項采用向后差分、對空間導數(shù)項采用中心差分，得到方程(1)的兩種情形的離散格式:

在x 軸上的空間節(jié)點(xi，r0，θk)處，控制方程(1)的離散格式如下:

在任意空間節(jié)點(xi，r0，θk)，j ＞0 處，控制方程(1)的離散格式如下:

其中:X(m)=

利用牛頓迭代法求解式(5)，牛頓迭代格式為:

其中，l 為迭代步數(shù)，初始向量［X(m)］(0)取上一時間層的近似解向量Xm-1。

2 局部收斂分析

根據(jù)文獻［13］，局部收斂分析是判斷數(shù)值方法優(yōu)劣的重要標準之一。下面的定理進一步說明了風力機流場數(shù)值模擬中的局部二階收斂性。

定理1:若時間步長滿足τ ≤min{τ1，τ2，τ3，}，則牛頓迭代格式(6)是局部二階收斂的。其中:

C1= ‖X*‖∞(X*是(5)的解)，C2表示密度ρ 在tm-1時刻求解區(qū)域Ω 內(nèi)的最小值。

證明:F'(X(m))于X*的領域S(X*，δ)∈Ω 內(nèi)連續(xù)G -可微，為了證明［F'(X(*))］-1存在，只需證明F'(X(*))為嚴格對角占優(yōu)矩陣。

記F'(X*)= alq，F(xiàn)'(X*)中的第l = 4［(k -1)N1N2+ (j -1)N1+ (i -1)］+1 行非零元素是:

若要使F'(X*)中的第l 行嚴格對角占優(yōu)，則只需

成立，即時間步長τ 滿足τ ≤［2A1C1］-1= τ1。

類似的，可得到時間步長滿足τ2時，F(xiàn)'(X*)中的第l + 1、l + 2 行嚴格對角占優(yōu)成立;滿足τ3時，F(xiàn)'(X*)中的第l + 3 行嚴格對角占優(yōu)成立。因此，當τ ≤min{τ1，τ2，τ3}時，F(xiàn)'(X*)嚴格對角占優(yōu)，［F'(X*)］-1存在。

對于?h ∈Rn，h ≠0，F(xiàn)″(X*)hh 的第l 個分量是:

其中m1= l -4NxNr，m2= l +4NxNr，m3= l -4Nr，m4= l +4Nr。

同理可得到F″(X*)hh 中的第l +1、l +2、l +3 個分量是非零的，從而我們有:F″(X*)hh ≠0。綜上，根據(jù)定理1 可知，牛頓迭代格式(7)是局部二階收斂的。

3 結(jié)語

1)建立了圓柱坐標系下的風力機流場的數(shù)學模型，即非線性、非定常的偏微分方程組定解問題。

2)將有限差分方法與牛頓法的結(jié)合運用到數(shù)值求解中，以其中一個方程為例給出相應的有限差分方法離散化格式與牛頓法迭代格式。

3)分析迭代格式的局部收斂性，給出局部收斂條件。

［1］張立茹，汪建文，劉冬冬，等.水平軸風力機風輪子午面流場結(jié)構的實驗研究［J］.太陽能學報，2011，32(3):323 -329.

［2］王果毅，董海濤.水平軸風力機的氣動數(shù)值模擬［J］.太陽能，2008(3):30 -33.

［3］楊從新，宋顯成.一種大型風力機葉片的氣動優(yōu)化設計方法［J］.空氣動力學學報，2011，29(02):222 -225.

［4］何玉林，任海軍，金鑫等.基于間接轉(zhuǎn)矩控制的風電機組功率測控平臺［J］.太陽能學報，2011，32(2):198 -203.

［5］Zhaoxue Cheng，Rennian Li.Criterion of aerodynamic performance of large scale offshore horizontal axis wind turbines［J］.Applied Mathematics and Mechanics，2010，31(1):13 -20.

［6］錢翼稷.空氣動力學［M］.北京:北京航空航天大學出版社，2005:1 -53.

［7］Versteeg H K，Malalasekera W.An introduction to Computational Fluid Dynamics［M］.The Finite Volume Method.Longman Prentice Hall Press，1995:85 -54.

［8］張慶麟.風力機葉片三維氣動性能的數(shù)值研究［J］.清華大學學報，2007，29(9):1172 -1174.

［9］葛永斌，田振夫，馬紅磊.三維泊松方程的高精度多重網(wǎng)格解法［J］.應用數(shù)學，2006，19(2):313 -318.

［10］傅作新，鄭雄.彈性與彈塑性問題的有限元與邊界元耦合解法［J］.計算數(shù)學，1993，10(4):375 -382.

［11］Hsu M H.Dynamic behaviour of wind turbine blades［J］.Proceedings of the Institution of Mechanical Engineers，Part C:Joumal of Mechanical Engineering Science，2008，222(8):1453 -1464.

［12］徐利治.幾個斂速階數(shù)較高的迭代程序和一個測定實根重數(shù)的方法［J］.數(shù)學學報，1962，12(4):331 -340.

［13］馮果忱.非線性方程組迭代解法［M］.上海:上?？茖W技術出版社，1989:112 -124.