[發明專利]一種單軸進給系統銑削過程的動力學集成建模方法

有效

• 申請號: 202010956712.6
• 申請日: 2020-09-12
• 公開（公告）號: CN112207331B
• 公開（公告）日: 2021-10-19
• 發明（設計）人: 張星;張偉;趙釗;趙萬華;尹佳
• 申請（專利權）人: 西安交通大學
• 主分類號: B23C1/00
• 分類號: B23C1/00;B23Q5/34
• 代理公司: 西安智大知識產權代理事務所 61215
• 代理人: 賀建斌

1.一種單軸進給系統銑削過程的動力學集成建模方法，其特征在于，包括以下步驟：

步驟1)控制系統模塊的等效建模：

控制系統模塊采用包含位置環、速度環和電流環的閉環控制方式，位置環控制器采用比例控制器，速度環、Q軸和D軸電流環控制器均采用比例-積分控制器；以數控系統輸出的位移指令為輸入，分別計算得到位置環、速度環、電流環的輸入與輸出值；

步驟2)PWM與逆變器模塊的等效建模：

基于SVPWM工作原理，根據控制系統模塊輸出的電壓矢量幅值建模和相位，確定6個絕緣柵雙極型晶體管(IGBT)的通斷狀態、通斷時間以及施加在電機定子繞組上的三相電壓值；

步驟3)伺服電機模塊的等效建模：

根據電機內部的電—磁—力—運動之間的耦合關系，分別對電學物理量定子三相電流、定子三相反電動勢、D軸與Q軸電流；磁場物理量轉子磁動勢、定子磁動勢、合成磁動勢；力學物理量電機輸出扭矩、扭矩角；以及運動學物理量工作臺實際位移、速度、加速度、電機轉子位移、速度進行計算；

步驟4)機械系統模塊的動力學等效建模：

分別對滾珠絲杠進給系統機械結構所包含的電機轉子21、工作臺22、滑塊23、導軌24、軸承25、螺母26、絲杠27和聯軸器28的進行等效，由牛頓—歐拉方程得到機械系統的動力學方程并進一步將其轉換為狀態空間表示形式，用Stribeck模型對機械系統摩擦力進行表征；

步驟5)銑削模塊的等效建模：

將零件視為剛性體，將刀具視為柔性體，建立兩自由度動力學模型，對刀具的振動狀態進行計算；考慮工作臺的實際位移和刀具的振動，對刀具-零件的實際嚙合狀態進行求解；根據刀具-零件的瞬時嚙合狀態，計算獲得瞬時切削力；

步驟6)模型的集成：

將控制系統模塊的輸出作為PWM與逆變器模塊的輸入，將PWM與逆變器模塊的輸出作為伺服電機模塊的輸入，將伺服電機模塊輸出的D、Q軸電流反饋到控制系統模塊，其輸出的電機扭矩作為機械系統模塊的輸入，機械系統模塊輸出各運動部件的運動狀態物理量，并將工作臺位移和電機轉子的速度反饋到控制系統模塊，將電機轉子的電角度反饋到控制系統模塊和伺服電機模塊，將工作臺位移輸出到銑削模塊；將銑削模塊計算獲得的切削力反饋到機械系統模塊，作用于工作臺；控制系統模塊、PWM與逆變器模塊、伺服電機模塊、機械系統模塊和銑削模塊集成為一個耦合的集成模型；

步驟7)對集成模型進行離散化：

將連續的時間離散成等間隔的時間步，將集成模型中關于時間的連續變量用當前仿真時刻對應的值進行代替，并將集成模型中所有的微分方程轉化為差分方程，將狀態空間的連續形式轉化為離散形式；

所述的控制系統模塊的位置環控制器采用比例控制器，其增益用K_p表示，設t為時間，x_r(t)為數控系統輸出的位移指令，θ_wta(t)為工作臺的實際位移輸出轉換到電機轉子端對應的轉角，則位置環控制器的輸出為

ω_r(t)＝K_p(x_r(t)η-θ_wta(t)) (1)

式中，η為螺母絲杠副的傳動比，ω_r(t)為速度指令值；

速度環控制器采用比例-積分控制器，其增益用K_v表示，時間常數用T_v表示，則速度環控制器的輸出為：

式中，ω_ra(t)為電機轉子的實際角速度值，i_r(t)為電流指令值；

速度環控制器的輸出將作為Q軸電流環的參考輸入，Q軸和D軸的電流環控制器均采用比例-積分控制器，其增益分別用K_iq和K_id表示，時間常數則分別用T_iq和T_id表示，則Q軸電流環控制器的輸出為

式中，i_qa(t)為電機Q軸電流的實際值，V_q(t)為Q軸電壓指令值；

D軸采用零輸入，因此D軸電流環控制器的輸出為：

式中，i_da(t)為電機D軸電流的實際值，V_d(t)為Q軸電壓指令值；

電流環控制器輸出值V_q(t)和V_d(t)均需經過限幅環節，其限定的幅值由實際控制伺服系統性能決定；設其經過限幅環節后的電壓值分別為V_qr(t)和V_dr(t)，該值為在兩相旋轉坐標系中表示的電壓值，對其進行如下坐標變換

得到在兩相靜止坐標系下的電壓值V_α(t)和V_β(t)，式中，θ_re(t)為電機的實際電角度值；進一步將兩相靜止坐標系下的電壓值用極坐標表示，得到控制系統模塊的最終輸出值為：

所述的PWM與逆變器模塊包含以下建模步驟：

控制系統模塊輸出的電壓矢量幅值V_m(t)及其相位θ_V(t)確定6個絕緣柵雙極型晶體管(IGBT)的通斷狀態，其中，T₂、T₄、T₆分別始終與T₁、T₃、T₅的通斷狀態相反，其具體關系如下：

3.1)當時：

狀態一：

通斷狀態：T₁接通，T₃、T₅斷開；

持續時間：

三相電壓：

其中，T_s為PWM開關周期，U_dc為直流電壓；

狀態二：

通斷狀態：T₁、T₃接通，T₅斷開；

持續時間：

三相電壓：

狀態三：

通斷狀態：T₁、T₃、T₅全部斷開；

持續時間：t₃＝T_s-t₁-t₂；

三相電壓：V_A(t)＝V_b(t)＝V_c(t)＝0；

3.2)當時：

狀態一：

通斷狀態：T₃接通，T₁、T₅斷開；

持續時間：

三相電壓：

狀態二：

通斷狀態：T₁、T₃接通，T₅斷開；

持續時間：

三相電壓：

狀態三：

通斷狀態：T₁、T₃、T₅全部斷開；

持續時間：t₃＝T_s-t₁-t₂；

三相電壓：V_A(t)＝V_b(t)＝V_c(t)＝0；

3.3)當時：

狀態一：

通斷狀態：T₃接通，T₁、T₅斷開；

持續時間：

三相電壓：

狀態二：

通斷狀態：T₃、T₅接通，T₁斷開；

持續時間：

三相電壓：

狀態三：

通斷狀態：T₁、T₃、T₅全部斷開；

持續時間：t₃＝T_s-t₁-t₂；

三相電壓：V_A(t)＝V_b(t)＝V_c(t)＝0；

3.4)當時：

狀態一：

通斷狀態：T₅接通，T₁、T₃斷開；

持續時間：

三相電壓：

狀態二：

通斷狀態：T₃、T₅接通，T₁斷開；

持續時間：

三相電壓：

狀態三：

通斷狀態：T₁、T₃、T₅全部斷開；

持續時間：t₃＝T_s-t₁-t₂；

三相電壓：V_A(t)＝V_b(t)＝V_c(t)＝0；

3.5)當時：

狀態一：

通斷狀態：T₅接通，T₁、T₃斷開；

持續時間：

三相電壓：

狀態二：

通斷狀態：T₁、T₅接通，T₃斷開；

持續時間：

三相電壓：

狀態三：

通斷狀態：T₁、T₃、T₅全部斷開；

持續時間：t₃＝T_s-t₁-t₂；

三相電壓：V_A(t)＝V_b(t)＝V_c(t)＝0；

3.6)當時：

狀態一：

通斷狀態：T₁接通，T₃、T₅斷開；

持續時間：

三相電壓：

狀態二：

通斷狀態：T₁、T₅接通，T₃斷開；

持續時間：

三相電壓：

狀態三：

通斷狀態：T₁、T₃、T₅全部斷開；

持續時間：t₃＝T_s-t₁-t₂；

三相電壓：V_A(t)＝V_b(t)＝V_c(t)＝0；

所述的伺服電機模塊包含以下建模步驟：

建立電機局部坐標系X_mY_mZ_m，原點位于伺服電機轉軸中心，Z_m軸與轉子軸線重合，X_m軸與電機定子A相繞組軸線方向一致，Y_m軸同時垂直于X_m軸和Z_m軸；

轉子永磁體在氣隙產生的磁通密度為

式中，B_r為永磁體剩余磁通密度，g為氣隙長度，h為永磁體厚度，μ_rm為永磁體的相對磁導率，則電機轉子每極的氣隙磁通為：

式中，D_r表示轉子直徑，L為轉子鐵心的軸向長度，P_n為電機轉子磁極對數；

電機轉子每極的磁動勢為

式中，μ₀為空氣磁導率，A_m為永磁體每極對應的氣隙面積；磁路總磁勢為：

F_r＝2F_p (10)

對于定子，其相電流由下式求得

式中，e_a(t)、e_b(t)、e_c(t)分別為定子三相繞組反電動勢，R為定子相電阻；

各相繞組每極產生的磁動勢為：

式中，k_w為繞組系數，用于考慮繞組分布效應；T_ph為每相繞組匝數；分別將各相磁動勢沿X_m軸和Y_m軸投影，得到定子合成磁場每極磁動勢在X_m軸、Y_m軸的分量為：

求得定子合成磁通沿X_m軸、Y_m軸的分量分別為：

將轉子磁通分別沿X_m軸和Y_m軸進行投影，得其分量為：

從而獲得轉子與定子合成的總磁通量在x軸和y軸的分量分別為：

則電機內部總磁通向量幅值及其相位為：

這里，θ_Φsr(t)的取值范圍為[-π,π]，其具體所在象限根據Φ_srx(t)和Φ_sry(t)的符號進行判斷；

進一步，求得定子三相繞組處的總磁通分別為：

則定子繞組的反電動勢為：

根據式(13)，計算得到定子合成磁動勢為：

則定子與轉子的總合成磁動勢通過下式求得：

其中，δ_sr(t)表示定子磁動勢軸線與轉子磁動勢軸線夾角且fix(·)表示向零取整；

電機輸出轉矩根據下式求得：

對于電機的d、q軸電流，通過對電機三相電流進行Park變換得到，即

電機轉子的機械角速度通過以下機械運動方程計算得到：

所述的機械系統模塊包含以下建模步驟：

滾珠絲杠進給系統動力學模型所考慮的結構包括電機轉子(21)、工作臺(22)、滑塊(23)、導軌(24)、軸承(25)、螺母(26)、絲杠(27)和聯軸器(28)；電機轉子(21)等效為轉動慣量J_r；工作臺(22)等效為質量塊m_wt；將絲杠(27)等分為三段并分別等效為三個質量m_s1、m_s2、m_s3和三個轉動慣量J_s1、J_s2、J_s3，并分別用絲杠拉壓剛度k_ssL、k_ssR和絲杠扭轉剛度k_θsL、k_θsR連接；將聯軸器(28)等效為扭轉彈簧單元，其扭轉剛度用k_θc表示；將螺母(26)等效為拉壓彈簧，其剛度用k_sn表示；將絲杠兩端軸承(25)分別等效為剛度為k_sbL和k_sbR的拉壓彈簧；考慮的阻尼環節包括導軌(24)與滑塊(23)之間的滑動阻尼c_wt，螺母(26)的平動阻尼c_sn和轉動阻尼c_θn，軸承(25)的平動阻尼c_sbL、c_sbR和轉動阻尼c_θbL、c_θbR，絲杠(27)的平動阻尼c_ssL、c_ssR和轉動阻尼c_θsL、c_θsR，聯軸器(28)的轉動阻尼c_θc以及電機轉子軸承阻尼c_θr；

由牛頓—歐拉方程，得到機械系統的動力學方程為

式中，M為機械系統的質量矩陣，且

K為機械系統的剛度矩陣，

C為機械系統的阻尼矩陣，

q(t)為機械系統的廣義坐標向量，

q(t)＝[x_a(t),x_s1(t),x_s2(t),x_s3(t),θ_s1(t),θ_s2(t),θ_s3(t),θ_r(t)]^T (29)

其中，x_wt為工作臺位移，x_s1、x_s2、x_s3分別為三段絲杠的平動位移，θ_s1、θ_s2、θ_s3分別為三段絲杠的扭轉角度，θ_r為電機轉子的扭轉角度；

F(t)為廣義力向量，

F(t)＝[F_f(t)+F_x(t) 0 0 0 0 0 0 T_e(t)]^T (30)

其中，F_x(t)為X方向的銑削力，F_f(t)為機械系統所受到的摩擦力，采用Stribeck模型計算如下：

其中，v_wt(t)為工作臺速度，F_c、F_s、V_s可通過實驗辨識得到，δ＝2；

進一步將機械系統動力學方程用狀態空間形式表示，選擇狀態變量x₁(t)＝q(t)，得機械系統的狀態空間表達式為：

其中，狀態向量：系統矩陣：常系數矩陣：輸出矩陣：C＝[I O]，輸入向量：u(t)＝F(t)，y(t)為輸出向量；

所述的銑削模塊的等效建模包括以下建模步驟：

(1)銑削過程動力學建模：

建立刀具局部坐標系O_tX_tY_tZ_t，坐標系原點O_t位于初始狀態下刀具底部中心處，坐標軸X_t與機床進給方向相同，Z_t為刀具軸線方向，Y_t垂直于X_tO_tZ_t平面；將零件視為剛體，將刀具簡化為兩自由度動力學系統，其振動方向分別與X_t軸、Y_t軸方向相同；

設某一時刻沿X_t軸、Y_t軸方向產生的銑削力分別為F_x(t)和F_y(t)，則刀具的振動狀態通過下式求得：

式中，m_xt、m_yt分別為刀具在X_t軸、Y_t軸方向的質量，c_xt、c_yt分別為X_t軸、Y_t軸方向的阻尼系數，k_xt、k_yt分別為X_t軸、Y_t軸方向的剛度，以上各物理量均可通過對刀具進行模態測試實驗獲得；x_t，y_t分別為刀具中心位置在X_t軸、Y_t軸的坐標；

(2)刀具-零件瞬時嚙合狀態與零件已加工表面的仿真計算：

設螺旋立銑刀的齒數為Z，螺旋角為β，半徑為R，對于刀齒均勻分布的銑刀來說，其齒間角φ_p表示為：

將刀具切削刃沿軸線方向分割成M個厚度為d_z的切削微元，并設初始時刻時，第1個刀齒端部切削微元的瞬時徑向接觸角為0°，則求得某一時刻t時，第j個刀齒上的第l個切削刃微元處的瞬時徑向接觸角為：

式中，S為主軸轉速；

對于工作臺做進給運動、刀具靜止的機床，在某一瞬時t，由步驟4)求得工作臺的位移為x_wt(t)，由式(33)求得刀具的瞬時振動位置為(x_t(t),y_t(t))，則可知在時刻t，在考慮工作臺實際位移和刀具振動的情況下，刀具中心在零件坐標系內的實際位置為：

進一步求得在時刻t，第j個刀齒上的第l個切削刃微元的坐標為

上式即為銑削加工過程中，螺旋銑刀各切削刃掃略過的空間軌跡曲面表達式，由其形成的包絡面即為零件的已加工表面；

在考慮工作臺實際位移和刀具振動的情況下，實際的未變形切屑厚度h_jl(t)為當前切削刃微元位置S_jl(t)與過該點和刀具軸線的連線與上一切削刃軌跡交點之間的距離；

(3)銑削力計算：

銑削力計算采用常用的瞬時剛性力計算模型，

式中：dF_t,jl，dF_r,jl，dF_a,jl分別表示切向、徑向和軸向切削力微元；ds表示切削刃微元長度；K_tc，K_rc，K_ac分別表示切向、徑向和軸向切削力系數；K_te，K_re，K_ae分別表示切向、徑向和軸向刃口力系數；切削力系數和刃口力系數通過實驗辨識得到；

考慮到切削刃微元高度較小，認為ds≈dz，則作用在刀齒j上的厚度為dz的第l個切削刃微元上的切向、徑向和軸向切削力微元表示為：

式中：g_jl(φ_jl(t))為單位階躍函數，用于表示當前切削刃微元是否參與切削，其定義為：

式中：φ_st，φ_ex——分別表示切入角和切出角；

對于逆銑，切入角和切出角的計算表達式分別為：

式中，a_e表示切削寬度；

對于順銑，切入角和切出角的計算表達式分別為：

通過坐標變換，得到作用在直角坐標系中的三個軸的切削力分量如下：

通過沿軸向積分和對每個刀齒求和，得到作用于整個銑刀上在進給、法線和軸向方向上的瞬時切削力分別為：

所述的模型集成包含以下步驟：

將控制系統模塊的輸出V_m(t)和θ_V(t)作為PWM與逆變器模塊的輸入；將PWM與逆變器模塊的輸出V_A(t)、V_B(t)、V_C(t)作為伺服電機模塊的輸入；將伺服電機模塊的輸出i_qa(t)、i_da(t)反饋到控制系統模塊，其輸出T_e(t)作為機械系統模塊的輸入；機械系統模塊輸出各運動部件的運動狀態物理量，并將工作臺位移和電機轉子的速度反饋到控制系統模塊，將電機轉子的電角度反饋到控制系統模塊和伺服電機模塊；將機械系統模塊計算獲得的工作臺位移x_wt(t)輸入到銑削模塊，并將銑削模塊計算獲得的銑削力F_x(t)反饋到機械系統模塊，作用于工作臺；

所述的離散化步驟為：

在離散時域內，將連續的時間t離散成N個等間隔的時間步Δt＝t/N，則對于任一時刻t_i，對機電集成模型做如下變換：

7.1)控制系統模塊離散化：

由式(1)得，位置環控制器的輸出為：

ω_r(t_i)＝K_p(x_r(t_i)η-θ_wta(t_i-1)) (46)

速度環控制器的輸出在離散時域內的表達式為：

其中，i_r0(t_i)＝i_r0(t_i-1)+(ω_r(t_i)-ω_ra(t_i-1))Δt；

同樣，Q軸和D軸電流環控制器的輸出在離散時域內的表達式也分別轉化為如下形式：

式(5)和式(6)分別轉化為下列各式：

7.2)伺服電機模塊離散化：

建立的三相交流永磁同步電機模型在轉化為時域離散模型進行仿真時，采用如下計算方法：

定子繞組的反電動勢為

在Δt時間內，定子三相繞組處的總磁通變化量為：

從而得到當前時刻定子三相繞組處的總磁通為：

根據三相磁通的大小，分別按下列三種情況進行計算：

當max{Φ_a(t_i),Φ_b(t_i),Φ_c(t_i)}＝Φ_a(t_i)時：

當max{Φ_a(t_i),Φ_b(t_i),Φ_c(t_i)}＝Φ_b(t_i)時：

當max{Φ_a(t_i),Φ_b(t_i),Φ_c(t_i)}＝Φ_c(t_i)時：

則得到總磁通向量幅值為：

其在X_m、Y_m軸的分量分別為：

轉子磁通在X_m、Y_m軸的分量由下式計算得到：

從而得到定子磁通在X_m、Y_m軸的分量為：

進一步得到定子磁動勢在X_m、Y_m軸的分量為：

則定子三相磁動勢按下式計算得到：

其中，F_so(t_i)為引入的虛擬項，F_so(t_i)≡0；

電機定子三相電流通過磁動勢直接求得：

定子合成磁動勢幅值及其相位為：

定子磁動勢軸線與轉子磁動勢軸線夾角為：

則定子與轉子的總合成磁動勢可通過下式求得

電機輸出轉矩根據下式求得

7.3)機械系統模塊離散化：

將式所表示的機械系統狀態空間方程轉換為離散形式：

其中，G＝e^AΔt,

7.4)銑削模塊離散化：

7.4.1)銑削過程動力學模型的離散求解：

將步驟5)獲得的微分方程形式的銑削過程動力學方程變換為差分形式，根據式(33)，在t_i時刻，刀具沿X_t軸和Y_t軸的瞬時加速度通過下式求得

進一步可得到其瞬時速度為

最終，其振動位移可通過下式獲得：

7.4.2)未變形切屑厚度的搜索求解：

在時域離散求解過程中，由式(37)得到的螺旋銑刀各切削刃掃略過的空間軌跡為一系列空間離散點，將步驟5)中的未變形切屑厚度用上一切削刃軌跡上與當前切削刃微元瞬時徑向接觸角φ_jl(t_i)相同位置附近，與當前切削刃微元距離最近點之間的距離代替，由于切削刃軌跡為離散點，因此通過在上一切削刃軌跡上瞬時徑向接觸角為的點附近進行搜索，獲得距離最近的點，并進一步獲得未變形切屑厚度，即

其中，T_ω為主軸旋轉一周的時間，K為一給定整數，用于指定搜索區間；

對于銑削力的計算，只需將相應表達式中的連續時間t用離散時間t_i替代即可。