2024中國可再生能源大會(huì)：風(fēng)機(jī)關(guān)鍵部件動(dòng)力學(xué)分析的精細(xì)積分法研究

風(fēng)機(jī)關(guān)鍵部件動(dòng)力學(xué)分析的精細(xì)積分法研究風(fēng)機(jī)關(guān)鍵部件動(dòng)力學(xué)分析的精細(xì)積分法研究2研究背景與現(xiàn)狀風(fēng)機(jī)基礎(chǔ)結(jié)構(gòu)的基本計(jì)算模型線性與非線性動(dòng)力學(xué)問題的精細(xì)積分法350050新增累計(jì)400裝機(jī)容量(GW)裝機(jī)容量(GW)300?我國風(fēng)電裝機(jī)容量新增、累計(jì)呈快速發(fā)展趨勢(shì)裝機(jī)容量(GW)裝機(jī)容量(GW)300200?我國海上風(fēng)電發(fā)展迅猛，近幾年裝機(jī)規(guī)模龐大20010000201420152016201720182019202020212022202320142015201620172018201920202021202220232014201520162017201820192020202120222023?風(fēng)電關(guān)鍵結(jié)構(gòu)尺寸規(guī)模逐年增大，建造成本變高?風(fēng)電機(jī)組荷載會(huì)對(duì)海上風(fēng)機(jī)的基礎(chǔ)結(jié)構(gòu)產(chǎn)生影響?海上風(fēng)電機(jī)組關(guān)鍵結(jié)構(gòu)常用梁模型進(jìn)行計(jì)算量縮減4風(fēng)機(jī)基礎(chǔ)結(jié)構(gòu)的基本計(jì)算以及算法改進(jìn)基于SubDyn展開單樁三角架單樁三角架SubDyn兼容兩類梁模型：Euler-Bernoulli梁單元以及Timoshenko梁單元，在通過單元?jiǎng)偠刃纬扇謩偠鹊臅r(shí)候，局部與全局坐標(biāo)系的示意圖如上圖；最終形成的動(dòng)力學(xué)方程可以形成如下的線性形式5!?Newmark法!?Newmark法?Wilson法?廣義u法?精細(xì)積分方法風(fēng)機(jī)柔性葉片基本計(jì)算模型風(fēng)力發(fā)電機(jī)計(jì)算仿真?國內(nèi)外開發(fā)了Bladed、OpenFAST、ENFAST、HAWC2等的各類風(fēng)力發(fā)電機(jī)的仿真軟件?OpenFAST中有專門進(jìn)行葉片動(dòng)力學(xué)仿真的動(dòng)力學(xué)模塊BeamDyn?葉片在空間中的形態(tài)具有較大的幾何非線性。對(duì)于柔性葉片而言，其模型具有大位移、大轉(zhuǎn)角甚至大應(yīng)變的特點(diǎn)?Euler-Bernoulli梁?Timoshenko梁?共旋坐標(biāo)法?絕對(duì)坐標(biāo)法?幾何精確梁模型6風(fēng)機(jī)柔性葉片基本計(jì)算模型風(fēng)機(jī)柔性葉片基本計(jì)算模型?將風(fēng)機(jī)葉片視為梁模型?荷載以外力荷載形式施加在葉片上?葉片通過錨點(diǎn)固結(jié)在中心轉(zhuǎn)子上ω錨點(diǎn)錨點(diǎn)!!轉(zhuǎn)子?在錨點(diǎn)處建立局部參考系，與錨點(diǎn)運(yùn)動(dòng)狀態(tài)一致，局部參考系中定義了梁的參考構(gòu)型?梁在局部坐標(biāo)系內(nèi)柔性變形，隨局部坐標(biāo)系在整體坐標(biāo)系中旋轉(zhuǎn)轉(zhuǎn)子隨葉片一起運(yùn)動(dòng)的坐標(biāo)系ω錨點(diǎn)隨葉片一起運(yùn)動(dòng)的坐標(biāo)系整體坐標(biāo)系隨葉片一起運(yùn)動(dòng)的坐標(biāo)系7?每個(gè)單元定義了一個(gè)共旋構(gòu)型?每個(gè)單元定義了一個(gè)共旋構(gòu)型?單元變形與受力分析是基于這個(gè)共旋構(gòu)型來計(jì)算的?共旋構(gòu)型與參考構(gòu)型之間的變換是剛體旋轉(zhuǎn)變換?適用于大旋轉(zhuǎn)大變形問題，但是單元還是小應(yīng)變的2EI0--2EI0--M1M2L02EI--M1M2L02EI--內(nèi)力通過旋轉(zhuǎn)變換至整體坐標(biāo)系中，有了整體坐標(biāo)下的KM=TTKM0T?位移變換一部分是全局位移轉(zhuǎn)換到單元坐標(biāo)系后的值，另一部分是由于單元坐標(biāo)系發(fā)生旋轉(zhuǎn)所產(chǎn)生的位移?整體坐標(biāo)中的剛度陣滿足與內(nèi)力的切線關(guān)系?pK=?u8設(shè)置時(shí)間步長度η,可以得到按時(shí)間步更新的設(shè)置時(shí)間步長度η,可以得到按時(shí)間步更新的迭代格式，在第k時(shí)間步有：Vtk+η=eAηVt+∫eAη一ξg(tk+ξ)dξ若在時(shí)間步內(nèi)，外荷載可以表現(xiàn)為線性的則非齊次項(xiàng)的積分可以積分得出，在第k時(shí)間步的迭代格式如下：r1線性與非線性動(dòng)力學(xué)問題的精細(xì)積分法精細(xì)積分方法（PIM）上式是一個(gè)離散后的動(dòng)力學(xué)方程，以如下方式引入狀態(tài)矢量V，表示為矩陣形式通過這種方法，可以將原方程化成如上的等價(jià)的一階微分方程形式V(t)=eAtV0+∫eA(t一ξ)g(ξ)dξ如此形式的微分方程，其解可以表示為：齊次方程的通解+非齊次項(xiàng)的特解?！襡A(t一ξ)f(ξ)dξ是Duhamel積分的形式9在時(shí)間步開始階段，在時(shí)間步開始階段，?fi+kiu的值為。1.首先假定時(shí)間步內(nèi)?fi+kiu項(xiàng)不發(fā)生改變，即認(rèn)為時(shí)間步內(nèi)的?fi+kiu為。2.計(jì)算r0和r1，通過迭代格式從vi計(jì)算vi+13.而后通過vi+1修正ki，再計(jì)算r0和r1，從vi重新計(jì)算vi+14.如此往復(fù)直到計(jì)算獲得的vi+1和上一循環(huán)的vi+1的差別在精度范圍內(nèi)。如此完成非線性問題的預(yù)估-矯正格式線性與非線性動(dòng)力學(xué)問題的精細(xì)積分法非線性動(dòng)力學(xué)問題的精細(xì)積分法若在精細(xì)積分的計(jì)算中，于時(shí)間步中考慮全局剛度與質(zhì)量的變化，考慮其非線性效應(yīng)：內(nèi)力項(xiàng)fi在時(shí)間步內(nèi)并不一定時(shí)刻滿足fi=kiu，全局剛度與內(nèi)力僅滿足切線關(guān)系ki=?fi/?ui。則可以參考這樣的思路采用類似的預(yù)估-矯正方法構(gòu)造求解方法：回到最初的動(dòng)力學(xué)方程，可以構(gòu)造出時(shí)間步內(nèi)的線性系統(tǒng)Mi+ciuǘ+kiu=Qt?fiu+kiu即在前文推導(dǎo)的基礎(chǔ)上，外荷載取得g(t)=tt?fiu+kiu)。M?1(Q()()通過多次積分的預(yù)估-矯正方法來對(duì)f進(jìn)行矯正，來構(gòu)造非線性的迭代格式?；赟ubDyn軟件，修改并發(fā)展SubDyn中的時(shí)程積分方法。計(jì)算左圖所示的基礎(chǔ)結(jié)構(gòu)：?單樁底部與海底固結(jié)?頂部接口處輸入為零?考慮單樁自重的影響分別采用龍格庫塔與精細(xì)積分方法在不同的時(shí)間步長下進(jìn)行計(jì)算，對(duì)于不同時(shí)間步長下各類方法的結(jié)果比較，以4階龍格庫塔（RK4）方法在0.001s時(shí)的計(jì)算結(jié)果作為參考解使用MATLAB進(jìn)行數(shù)據(jù)的可視化?在時(shí)間步長為0.001s時(shí)，PIM的計(jì)算?在時(shí)間步長為0.0015s時(shí)，這是RK4所能計(jì)算的最大時(shí)間步長，P?在時(shí)間步長為0.003s時(shí)，PIM的計(jì)算三腳架模型以及離散計(jì)算簡圖基于SubDyn軟件，修改并發(fā)展SubDyn中的時(shí)程積分方法。計(jì)算左圖所示的基礎(chǔ)結(jié)構(gòu)：?底部固結(jié)?頂部輸入x=0.1m?考慮單樁自重影響分別采用龍格庫塔與精細(xì)積分方法在不同的時(shí)間步長下進(jìn)行計(jì)算，并對(duì)比結(jié)果；對(duì)于不同時(shí)間步長下各類方法的結(jié)果比較，以4階龍格庫塔（RK4）方法在0.005s時(shí)的計(jì)算結(jié)果作為參考解?在時(shí)間步長為0.005s時(shí)，PIM的計(jì)?在時(shí)間步長為0.1s時(shí)，PIM的計(jì)三腳架模型以及離散計(jì)算簡圖基于SubDyn軟件，修改并發(fā)展SubDyn中的時(shí)程積分方法。計(jì)算左圖所示的基礎(chǔ)結(jié)構(gòu)：?底部樁腿固結(jié)?頂部輸入外力為5sin(10t)?并考慮單樁自重分別采用龍格庫塔與精細(xì)積分方法在不同的時(shí)間步長下進(jìn)行計(jì)算，并對(duì)比結(jié)果；對(duì)于不同時(shí)間步長下各類方法的結(jié)果比較，以4階龍格庫塔（RK4）方法在0.005s時(shí)的計(jì)算結(jié)果作為參考解?在時(shí)間步長為0.005s時(shí)，PIM的計(jì)算結(jié)果與參考解一致，驗(yàn)證了精細(xì)積分法的準(zhǔn)確性?在時(shí)間步長為0.05s時(shí)，為RK4能計(jì)算的最大時(shí)間步長，RK4的計(jì)算結(jié)果略微偏差，PIM計(jì)算結(jié)果準(zhǔn)確?在時(shí)間步長為0.1s時(shí)，PIM的計(jì)算結(jié)果仍與參考解一致，說明了精細(xì)積分法的穩(wěn)定性?基于SubDyn的精細(xì)積分法能計(jì)算更大時(shí)間步長的數(shù)值積分，在求解大規(guī)模多自由度問題有顯著優(yōu)勢(shì)若先假設(shè)u(t)的解析形式，則可得出外力Qt的對(duì)應(yīng)形式假設(shè)u(t)具有如下解析形式u(t)=Lsinwt代入動(dòng)力學(xué)微分方程，可得Q(t)的形式如下=—MLw2sin(wt)+K1Lsinwt+K2L3sin3wtQt=MLw2sin(wt)+K1Lsinwt+K2L3sin3wt在后續(xù)具體計(jì)算模擬過程中取得參數(shù)設(shè)置Y=設(shè)置Y=0.5,β=0.25，產(chǎn)生了較大的數(shù)值阻尼其計(jì)算的角速度較解析解明顯變慢，同時(shí)振幅也很快衰竭使用廣義u法對(duì)此問題進(jìn)行計(jì)算使用與Newmark法中相同的Y,β參數(shù)，并取廣義α參數(shù)P=0，其計(jì)算結(jié)果與Newmark法類似使用廣義u法對(duì)此問題進(jìn)行計(jì)算使用與Newmark使用廣義u法對(duì)此問題進(jìn)行計(jì)算使用與Newmark法中相同的y,β參數(shù)，并取廣義α參數(shù)P=0.5，其計(jì)算結(jié)果有所改善效果較Newmark法與廣義α法好沒有出現(xiàn)角速度與振幅的偏移共旋坐標(biāo)梁的靜力學(xué)算例設(shè)置共旋坐標(biāo)梁的靜力學(xué)算例設(shè)置?設(shè)置梁長度為10m，采用φ20Q235圓鋼的材料參數(shù)；?使梁的最左側(cè)節(jié)點(diǎn)固結(jié)于空間，約束其全部的自由度；?施加彎矩于最右側(cè)節(jié)點(diǎn)，彎矩大小EIz，為其彎成圓周的理論彎矩大??；?設(shè)置30荷載步，通過NR迭代方法求解共旋坐標(biāo)梁的靜力在靜力學(xué)與動(dòng)力學(xué)求解過程中，由于模型的幾何非線性，在形成全局坐標(biāo)下的剛度陣時(shí)，需要記錄單元的旋轉(zhuǎn)歷史。?形成10節(jié)點(diǎn)的有限元梁模型：采用共旋歐拉梁的單元?jiǎng)偠?，?×6的阻尼陣C為對(duì)角陣，于對(duì)角線取得阻尼值?在最右側(cè)節(jié)點(diǎn)施加大小為彎矩大小EIz，為梁右端節(jié)點(diǎn)旋轉(zhuǎn)角度為90°的理論彎矩；?時(shí)間步長為0.01s，總時(shí)間步長為2000步；外力矩在500時(shí)間步內(nèi)線性加載至約定的外力矩M。πMM計(jì)算結(jié)構(gòu)簡圖M0彎矩施加a)梁變形后的構(gòu)型圖b)精細(xì)積分法的計(jì)算過程a)梁變形后的構(gòu)型圖b)精細(xì)積分法的計(jì)算過程?通過MATLAB進(jìn)行編程，實(shí)現(xiàn)共旋坐標(biāo)梁模型，并采用精細(xì)積