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