兩種變厚度空心儲能飛輪的應力特性

蘭 晨，李文艷

（華北電力大學能源動力與機械工程學院，北京102206）

近年來，我國積極應對全球氣候變化，踐行《巴黎協(xié)定》，為實現(xiàn)“碳中和”而大力發(fā)展可再生能源，特別是太陽能、風能裝機功率在2030 年將達到12 億千瓦。但是可再生能源受各類因素所限制而有著較大的隨機性與波動性，其并網(wǎng)過程會給電網(wǎng)帶來一定的沖擊[1]。應用儲能裝置可以很好應對可再生能源給電網(wǎng)帶來的影響，其中飛輪儲能由于充能速度快、響應時間短、占地小、無污染和使用壽命長等優(yōu)點在新能源并網(wǎng)、UPS、石油石化、軌道交通等領域得到了廣泛的運用[2-5]。

儲能飛輪根據(jù)其拓撲結構主要可以分為4 種：內(nèi)飛輪內(nèi)轉(zhuǎn)子結構、分體式結構、內(nèi)轉(zhuǎn)子外飛輪結構、外轉(zhuǎn)子外飛輪結構[6]。其中內(nèi)轉(zhuǎn)子外飛輪結構可以視為輪輻加輪緣的變厚度空心飛輪結構，由于其大部分質(zhì)量分布在飛輪外緣，使得其具有較大的儲能密度，因此在商業(yè)上得到了廣泛的應用[7]。目前，已有眾多學者對此類飛輪結構進行了分析研究與優(yōu)化設計。蘇芳等[8]基于有限元軟件分析研究了空心飛輪轉(zhuǎn)子使用不同材料時飛輪徑向、環(huán)向應力的變化規(guī)律。任正義等[9]應用Ansys Workbench 軟件對3種不同形式的空心鋁合金飛輪轉(zhuǎn)子模型進行有限元分析，研究了3 種形式空心飛輪轉(zhuǎn)子的應力、變形分布情況，并對曲線輪輻飛輪進行了優(yōu)化。閆曉磊等[10]采用最優(yōu)控制理論，得到空心飛輪轉(zhuǎn)子的最優(yōu)形狀解析表達式，并對空心飛輪轉(zhuǎn)子在低速、中速和高速情況下進行了最優(yōu)化設計。以上飛輪轉(zhuǎn)子應力分析都建立在平面應力狀態(tài)假設的基礎之上，并未考慮飛輪輪緣對應力分布，特別是軸向應力的影響，在飛輪輪緣高度較大時不能有效確定飛輪應力分布。因此本文建立了兩種變厚度空心飛輪模型，在空間應力狀態(tài)假設基礎上，利用Ansys Workbench有限元分析軟件分析了輪緣高度變化對兩種飛輪模型徑向、環(huán)向、軸向應力最大值及最大變形量的影響，并給出了沿給定路徑的各項應力分布。

1 飛輪應力理論

目前的飛輪結構優(yōu)化與應力研究大多建立在盤狀飛輪的基礎上[8-10]，當飛輪外徑R1與飛輪厚度h滿足式(1)[11]時

飛輪應力狀態(tài)可以視為平面應力狀態(tài)，然而當飛輪結構不滿足式(1)時，例如功率型儲能飛輪大多采用半徑小、高度大、轉(zhuǎn)動慣量小，近似于空心長圓筒型的結構[4]，此時不能采用平面應力理論計算飛輪的應力，將飛輪應力分解時需要考慮軸向應力，根據(jù)彈性力學理論，飛輪高速旋轉(zhuǎn)產(chǎn)生的離心力可沿圓柱坐標系3個方向分解為徑向應力、環(huán)向應力與軸向應力，對于各項同性的等厚度空心圓筒，其在半徑為r 處的徑向應力，環(huán)向應力與軸向應力分別為

式中，σr為飛輪徑向應力，Pa；σθ為飛輪環(huán)向應力，Pa；σz為飛輪軸向應力，Pa；εz為飛輪軸向應變；ρ為飛輪材料密度，kg/m3；ω為飛輪旋轉(zhuǎn)角速度，rad/s；v 為材料泊松比；R0為飛輪內(nèi)徑，m；R1為飛輪外徑，m。

且明顯可知飛輪環(huán)向應力始終大于徑向應力與軸向應力，采用工程設計Tresca 屈服準則時，僅需滿足

式中，[σ]為材料許用應力，Pa。

變厚度結構可將大部分質(zhì)量集中在飛輪外徑處。變厚度結構較于等厚度結構，其優(yōu)勢在于采用變厚度結構的飛輪有著較大的有效回轉(zhuǎn)半徑，等厚度飛輪與變厚度結飛輪結構如圖1所示，其有效回轉(zhuǎn)半徑分別為

當h=h2時，式(9)減去式(10)有

式中，h 為等厚度飛輪厚度；h1為飛輪輪盤厚度，m；h2為飛輪輪緣厚度，m；a 為飛輪內(nèi)徑，m；b為飛輪外徑，m；c為飛輪輪緣內(nèi)徑。

圖1 等厚度飛輪結構圖(a)與變厚度飛輪結構圖(b)Fig.1 Structure diagram of equal thickness flywheel(a)and variable thickness flywheel(b)

而飛輪儲能密度通常由飛輪有效回轉(zhuǎn)半徑?jīng)Q定，通常情況下飛輪儲能密度指質(zhì)量儲能密度，飛輪質(zhì)量儲能密度可表示為

式中，E為飛輪儲能量，J；I為飛輪轉(zhuǎn)動慣量，kg/m2；ω為飛輪角速度，rad/s；R為飛輪有效回轉(zhuǎn)半徑，m；K為飛輪形狀因子；[σ]為許用應力，Pa。

圖2 給出了飛輪形狀因子與飛輪結構的關系，形狀因子K越大代表飛輪儲能性能越好。

圖2 不同飛輪結構的形狀因子Fig.2 Shape factors of different flywheel structures

由圖2 可以得知等厚度結構飛輪形狀因子為0.305，而變厚度結構飛輪形狀因子介于0.305～0.5 之間。由上述分析可知變厚度結構有著比等厚度結構更好的儲能特性。

2 變厚度空心飛輪模型構建

在飛輪角速度與材料都確定時，飛輪應力僅與其結構有關，建立如圖3的轉(zhuǎn)子模型，其中X1代表軸孔半徑、X2代表輪輻半徑、X3代表輪緣厚度、X4代表輪緣高度、X5代表輪輻厚度，兩種轉(zhuǎn)子模型的結構參數(shù)見表1。將兩模型導入有限元軟件Ansys Workbench 中并施加10000 r/min 的旋轉(zhuǎn)速度與合適的約束條件，為了同時考慮計算精度與計算速度，選取MultiZone 劃分方法，網(wǎng)格長度取2 mm，劃分網(wǎng)格后的飛輪轉(zhuǎn)子模型如圖4所示；為了分析輪緣高度對飛輪輪輻上應力分布的影響，建立從飛輪內(nèi)徑到外徑的路徑如圖3 中A-B。飛輪材料使用7075鋁合金，其參數(shù)見表2。

圖3 飛輪模型一(a)與飛輪模型二結構圖(b)Fig.3 Structure diagram of flywheel model 1(a)and flywheel model 2(b)

圖4 X4=100 mm時飛輪模型一的有限元網(wǎng)格劃分Fig.4 Meshing of flywheel model 1 when X4=100 mm

表1 飛輪轉(zhuǎn)子模型參數(shù)Table 1 Model parameters of flywheel rotor

表2 飛輪材料參數(shù)Table 2 flywheel material parameters

3 飛輪轉(zhuǎn)子模型應力分析

由表1可知，兩模型除了輪緣高度其他參數(shù)都相同。由理論力學可知，當模型二輪緣高度為模型一輪緣高度的兩倍，且兩飛輪模型以相同角速度繞中心旋轉(zhuǎn)軸旋轉(zhuǎn)時，兩者轉(zhuǎn)動慣量相等。將模型二輪緣高度乘0.5 得到模型二的折算輪緣高度，為了方便說明，之后也稱輪緣高度。

由圖3可知，輪緣的存在會使得飛輪厚度突變并必然導致應力集中。為了分析應力集中對飛輪整體應力的影響，通過逐步增大輪緣高度得到在不同輪緣高度下兩種飛輪模型的最大徑向、最大環(huán)向、最大軸向應力值與最大變形量。最大應力值可表征飛輪模型應對應力集中的能力，最大變形量則反映了飛輪模型應對變形的能力。兩飛輪模型各項最大應力值隨輪緣高度變化趨勢見圖5～7，比值為模型一應力值與模型二對應應力值之比。由圖5 可知，隨著輪緣高度的增加，模型一最大徑向應力值單調(diào)增長，隨后在輪緣高度為60 mm 時達到穩(wěn)定值49.1 MPa；模型二最大徑向應力值先增后減，分別在輪緣高度為10 mm 和30 mm 時達到最大值32 MPa 和極小值28.86 MPa，并在輪緣高度為60 mm時達到穩(wěn)定值29.5 MPa。比值先減小并于輪緣高度為3.75 mm 時取得最小值0.698，隨后增長至輪緣高度為60 mm時達到穩(wěn)定值1.65，輪緣高度為17.5 mm時比值為1，此時兩模型應力值相等。

圖5 飛輪模型的最大徑向應力與應力比值Fig.5 Maximum radial stress and stress ratio of flywheel model

圖6 飛輪模型的最大環(huán)向應力與應力比值Fig.6 Maximum circumferential stress and stress ratio of flywheel model

由圖6可知，最大環(huán)向應力與最大徑向應力變化趨勢相似。隨著輪緣高度增加，模型一最大環(huán)向應力值增大，并在輪緣高度為60 mm 時達到穩(wěn)定值67 MPa；模型二最大環(huán)向應力值先增后減，最終在輪緣高度為60 mm 時達到穩(wěn)定值59.1 MPa，其中最大值63.1 MPa 和極小值58.32 MPa 分別在輪緣高度為10 mm 和30 mm 時取得。比值在輪緣高度為3.75 mm 時取得最小值0.745，隨后升高并在輪緣高度為22.55 mm 時達到1，最終在輪緣高度達到60 mm時達到穩(wěn)定值1.133。

圖7 飛輪模型的最大軸向應力與應力比值Fig.7 Maximum axial stress and stress ratio of flywheel model

如圖7所示，兩模型的最大軸向應力值相差較大，隨著輪緣高度增加，模型一最大軸向應力明顯增大，其穩(wěn)定值47.5 MPa 出現(xiàn)在輪緣高度為80 mm時；模型二的最大軸向應力值在輪緣高度為10 mm 處取得最大值16.3 MPa，之后應力值隨輪緣高度增加而減小，并在輪緣高度為60 mm時達到穩(wěn)定值8.96 MPa。比值在輪緣高度為2.5 mm時取得最小值0.52，隨后快速增長至5.3并達到穩(wěn)定。

對于最大變形量，如圖8所示，隨著輪緣高度增加，模型一最大變形量不斷增大但增速逐漸變緩，當輪緣高度為150 mm 時，最大變形量為50.447 μm；模型二最大變形量變化呈增-減-增趨勢，在輪緣高度為10 mm 與45 mm 時分別取得極大值89.902 μm 與極小值68.471 μm，兩種飛輪模型最大變形量最終都不斷增大。應力比值分別于輪緣高度為7.5 mm和90 mm時取得極小值0.2254和極大值0.6077，且始終不大于1，表明模型二最大變形量始終大于模型一最大變形量。

由上述分析可得，在相同載荷和材料條件下，輪緣高度很小時，模型一應力性能略優(yōu)于模型二，但輪緣高度較大時，模型二的應力性能遠優(yōu)于模型一。相對的，由于模型二的輪緣集中在一側，使其最大變形量始終大于模型一的最大變形量，變形量比值還將隨著輪緣高度增大而不斷降低。并且注意到輪緣高度超過某一臨界值后，兩種飛輪的最大應力值都將達到穩(wěn)定。

圖8 飛輪模型的最大變形量與變形量比值Fig.8 Maximum deformation and deformation ratio of flywheel model

圖9 模型一沿路徑方向的徑向應力分布Fig.9 Radial stress distribution along path of model 1

圖10 模型二沿路徑方向的徑向應力分布Fig.10 Radial stress distribution along path of model 2

為了進一步研究在輪緣高度對應力在飛輪半徑方向上變化的影響，結合在不同輪緣高度下飛輪模型最大應力值的變化情況，選取7個合適的輪緣高度值，分析兩種飛輪模型的各項應力值在給定的輪緣高度下沿圖3 中路徑A-B 的變化情況。如圖9 和圖10 所示，兩飛輪模型沿路徑方向的徑向應力都呈先增后減的趨勢，并且在路徑長度為0 mm 與100 mm，即飛輪模型的內(nèi)徑與外徑處，徑向應力值都接近0，這與式(2)相符。如圖11和圖12所示，徑向應力最大值出現(xiàn)在飛輪內(nèi)徑，環(huán)向應力沿路徑方向單調(diào)遞減。如圖13，路徑長度小于40 mm時，模型一軸向應力接近0；路徑長度為60 mm時，軸向應力取得最大值。由圖14 可知，模型二軸向應力在路徑長度小于35 mm 時接近0，路徑長度為55 mm 時取得最大值，路徑長度為70 mm 時取得負的極小值，代表此處取得反向的應力極值。由上述分析可知，隨著輪緣高度的不斷增大，徑向應力與環(huán)向應力隨之增大，但變化趨勢沒有發(fā)生明顯改變，僅在厚度突變處(60 mm)應力曲線出現(xiàn)了略微的彎曲，代表此處應力降低速度變快；軸向應力受輪緣影響較大，應力曲線在輪緣附近出現(xiàn)顯著的凸起，代表軸向應力在此處突增，說明在輪緣高度較大的情況下，飛輪轉(zhuǎn)子不能簡單地利用平面應力理論進行分析，否則可能會導致較大的誤差。

圖11 模型一沿路徑方向的環(huán)向應力分布Fig.11 Circumferential stress distribution along path of model 1

圖12 模型二沿路徑方向的環(huán)向應力分布Fig.12 Circumferential stress distribution along path of model 2

圖13 模型一沿路徑方向軸向應力分布Fig.13 Axial stress distribution along path of model 1

圖14 模型二沿路徑方向軸向應力分布Fig.14 Axial stress distribution along path of model 2

如圖9～14 所示，輪緣高度達到50 mm 后，再增大輪緣高度對兩種飛輪模型各項應力的影響變小，輪緣高度為50 mm和100 mm時的應力變化曲線幾乎重合，這與飛輪整體最大應力相似，超過輪緣高度臨界值后路徑方向應力也達到穩(wěn)定值。為進一步得出路徑方向達到應力穩(wěn)定時對應的臨界輪緣高度，逐漸增大輪緣高度，得到不同輪緣高度下兩種飛輪模型沿路徑方向的各項最大應力值，如圖15 所示，隨著輪緣高度增加，模型一沿路徑方向的最大徑向、最大環(huán)向、最大軸向應力分別在輪緣高度為60、60、70 mm時達到穩(wěn)定值25.9、67、12 MPa；模型二沿路徑方向的最大徑向、最大環(huán)向應力、最大軸向應力則分別在輪緣高度為60、60、40 mm 時達到穩(wěn)定值16.7、46.3、0.8 MPa。通過上述分析可以得知，模型一沿路徑的最大徑向、最大環(huán)向、最大軸向應力穩(wěn)定值比模型二對應應力穩(wěn)定值分別大54.3%、44.4%和1420%。通過對比圖15 與圖5～7 可以發(fā)現(xiàn)，隨著輪緣高度增大，無論是受應力集中影響較大的飛輪整體最大應力，還是受應力集中影響較小的沿路徑方向的飛輪應力都最終達到穩(wěn)定狀態(tài)。

圖15 飛輪模型沿路徑的最大應力Fig.15 Maximum radial stress and stress ratio along path of flywheel model

4 結 論

（1）隨著輪緣高度增大，飛輪模型一的最大應力值與沿路徑的最大應力值都將上升；飛輪模型二的最大應力值與沿路徑的最大軸向應力隨著輪緣高度增加先增大后減小，而沿路徑的最大徑向、最大環(huán)向應力值單調(diào)增長。輪緣高度超過臨界高度后，兩種飛輪模型的各項應力值都將達到穩(wěn)定值。在本文的飛輪轉(zhuǎn)子模型和載荷條件下輪緣高度臨界值約為60 mm，針對不同的飛輪轉(zhuǎn)子模型需要根據(jù)實際情況確定輪緣高度臨界值。

（2）在不同輪緣高度下兩種飛輪模型沿路徑方向的徑向應力、環(huán)向應力變化曲線與理論情況相符，徑向應力呈先增后減的趨勢；環(huán)向應力單調(diào)減少，最大環(huán)向應力出現(xiàn)在飛輪內(nèi)壁。需要注意的是，兩種飛輪模型沿路徑方向的軸向應力值在無輪緣部分都接近零，但是靠近有輪緣部分后軸向應力明顯增大，并在厚度突變處附近取得最大軸向應力值，隨后應力下降并反向，對于飛輪模型二，其軸向應力下降后還會再次上升回到正值。

（3）輪緣存在時飛輪應力不能完全視為平面應力，在工程設計中需要適當考慮軸向應力的影響，在本文中應力值穩(wěn)定后，兩個飛輪模型的軸向應力最大值分別為其環(huán)向應力最大值的70.3% 和15.16%，若忽略軸向應力則會造成較大的誤差。

（4）穩(wěn)定后，模型一的最大徑向應力、環(huán)向應力、軸向應力比模型二對應應力分別大65%、13.3%、430%，沿路徑的最大徑向、環(huán)向、軸向應力穩(wěn)定值比模型二對應應力分別大54.3%、44.4%、1420%，最大變形量比模型二的最大變形量小43.89%，這說明采用模型一在控制變形上更有優(yōu)勢，而采用模型二能更有效降低應力。

由上述結論可得，通過確定飛輪輪緣高度臨界值后可以確定其應力穩(wěn)定值。在相同的條件下，采用類似模型一的飛輪結構對材料性能要求更加嚴格，但采用高彈性模量材料滿足應力要求時可以有效減小飛輪變形，降低飛輪變形造成安全事故的可能性；采用類似模型二的飛輪結構可以選用彈性模量較低的材料以降低制造成本，或者適當增加飛輪轉(zhuǎn)速以提高飛輪儲能量，但必須注意輪緣變形帶來的影響。對于兩種模型而言，輪緣高度超出臨界高度后但變形量仍留有較大裕度的情況下，可以合理增加輪緣高度以提高轉(zhuǎn)動慣量。