蘭 晨，李文艷

（華北電力大學(xué)能源動(dòng)力與機(jī)械工程學(xué)院，北京102206）

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

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

1 飛輪應(yīng)力理論

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

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

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

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

式中，[σ]為材料許用應(yīng)力，Pa。

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

當(dāng)h=h2時(shí)，式(9)減去式(10)有

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

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

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

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

圖2 給出了飛輪形狀因子與飛輪結(jié)構(gòu)的關(guān)系，形狀因子K越大代表飛輪儲(chǔ)能性能越好。

圖2 不同飛輪結(jié)構(gòu)的形狀因子Fig.2 Shape factors of different flywheel structures

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

2 變厚度空心飛輪模型構(gòu)建

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

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

圖4 X4=100 mm時(shí)飛輪模型一的有限元網(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)子模型應(yīng)力分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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