Ansys 計算剛柔耦合系數矩陣
文章目錄
- Ansys 計算剛柔耦合系數矩陣
- 衛星的剛柔耦合動力學模型
- 采用 ANSYS 的 APDL 語言的計算方法
- 系統轉動慣量的求解方法
- 參考文獻
衛星的剛柔耦合動力學模型
-
柔性航天器的剛柔耦合動力學模型可以表示為
m v ˙ + B t r a n η ¨ = F J ω ˙ + ω × J ω + B r o t η ¨ + ω × B r o t η ˙ = T η ¨ + 2 ξ Ω η ˙ + Ω 2 η + B r o t T ω ˙ + B t r a n T v ˙ = 0 \begin{align} \boldsymbol{m}\dot{\boldsymbol{v}} + \boldsymbol{B}_{tran}\ddot{\boldsymbol{\eta}} &= \boldsymbol{F}\\ \boldsymbol{J}\dot{\boldsymbol{\omega}}+\boldsymbol{\omega}^\times\boldsymbol{J}\boldsymbol{\omega}+\boldsymbol{B}_{rot}\ddot{\boldsymbol{\eta}}+\boldsymbol{\omega}^\times\boldsymbol{B}_{rot}\dot{\boldsymbol{\eta}}&=\boldsymbol{T}\\ \ddot{\boldsymbol{\eta}}+2\boldsymbol{\xi}\boldsymbol{\Omega}\dot{\boldsymbol{\eta}}+\boldsymbol{\Omega}^2\boldsymbol{\eta}+\boldsymbol{B}_{rot}^T\dot{\boldsymbol{\omega}} + \boldsymbol{B}_{tran}^T\dot{\boldsymbol{v}}&=\boldsymbol{0} \end{align} mv˙+Btran?η¨?Jω˙+ω×Jω+Brot?η¨?+ω×Brot?η˙?η¨?+2ξΩη˙?+Ω2η+BrotT?ω˙+BtranT?v˙?=F=T=0??
式中：- m \boldsymbol{m} m， J J J 分別表示柔性航天器的質量和轉動慣量矩陣；
- v \boldsymbol{v} v， ω \boldsymbol{\omega} ω 分別表示航天器本體相對于慣性系的線速度和角速度矢量；
- F \boldsymbol{F} F， T \boldsymbol{T} T 分別表示航天受到的三軸控制力和控制力矩；
- η \boldsymbol{\eta} η， ξ \boldsymbol{\xi} ξ， Ω \boldsymbol{\varOmega} Ω 分別表示航天器的模態位移向量，柔性振動的阻尼比和頻率矩陣；
- B t r a n = A ? B t r a n b \boldsymbol{B}_{tran}=A\cdot \boldsymbol{B}_{tran}^b Btran?=A?Btranb? 為柔性附件在航天器坐標系中相對于航天器本體的平動耦合系數矩陣， A \boldsymbol{A} A為柔性附件坐標系到衛星坐標系的轉移矩陣， B t r a n b \boldsymbol{B}_{tran}^b Btranb? 為柔性附件在自身坐標系中相對于航天器本體的平動耦合系數矩陣；
- B r o t = l p × B t r a n + A B r o t b \boldsymbol{B}_{rot} = \boldsymbol{l}_p^\times \boldsymbol{B}_{tran} + \boldsymbol{A}\boldsymbol{B}_{rot}^b Brot?=lp×?Btran?+ABrotb? 為柔性附件在航天器坐標系中相對于坐標原點的轉動耦合矩陣， l p \boldsymbol{l}_p lp? 為由航天器質心指向柔性附件坐標原點的向量， × ^\times × 表示反對稱矩陣， B r o t b \boldsymbol{B}_{rot}^b Brotb? 為柔性附件在自身坐標系中相對于坐標原點的轉動耦合系數。
-
B t r a n \boldsymbol{B}_{tran} Btran? 和 B r o t \boldsymbol{B}_{rot} Brot? 的表達式如式(4)和式(5)所示
B t r a n = [ B t r a n 1 B t r a n 2 ? B t r a n i ? B t r a n m ] T B r o t = [ B r o t 1 B r o t 2 ? B r o t i ? B r o t m ] \begin{align} \boldsymbol{B}_{tran} &= \begin{bmatrix} \boldsymbol{B}_{tran}^1 & \boldsymbol{B}_{tran}^2 & \cdots & \boldsymbol{B}_{tran}^i & \cdots & \boldsymbol{B}_{tran}^m \end{bmatrix}^T \\ \boldsymbol{B}_{rot} &= \begin{bmatrix} \boldsymbol{B}_{rot}^1 & \boldsymbol{B}_{rot}^2 & \cdots & \boldsymbol{B}_{rot}^i & \cdots & \boldsymbol{B}_{rot}^m \end{bmatrix} \end{align} Btran?Brot??=[Btran1??Btran2????Btrani????Btranm??]T=[Brot1??Brot2????Broti????Brotm??]??
其中， m m m 為模態截斷階數；每個元素 B t r a n / r o t i \boldsymbol{B}_{tran/rot}^i Btran/roti? 為 3 × 1 3\times1 3×1 的向量。 -
B t r a n i \boldsymbol{B}_{tran}^i Btrani? 和 B r o t i \boldsymbol{B}_{rot}^i Broti? 的表達式如式(6)和式(7)所示
B t r a n i = ∑ k = 1 n m k Φ k i = ∑ k = 1 n m k [ u k x i u k y i u k z i ] B r o t i = ∑ k = 1 m m k r k × Φ k i = ∑ k = 1 n m k [ 0 ? r k z r k y r k z 0 ? r k ? r k y r k x 0 ] [ u k x i u k y i u k z i ] \begin{align} \boldsymbol{B}_{tran}^i &= \sum_{k=1}^nm_k\boldsymbol{\varPhi}_k^i = \sum_{k=1} ^n m_k\begin{bmatrix} u_{kx}^i \\ u_{ky}^i \\ u_{kz}^i\end{bmatrix} \\ \boldsymbol{B}_{rot}^i &= \sum_{k=1}^m m_k \boldsymbol{r}_k^\times\boldsymbol{\varPhi}_k^i = \sum_{k=1}^n m_k \begin{bmatrix} 0 & -r_{kz} & r_{ky} \\ r_{kz} & 0 & -r_k \\ -r_{ky} & r_{kx} & 0 \end{bmatrix} \begin{bmatrix} u_{kx}^i \\ u_{ky}^i \\ u_{kz}^i \end{bmatrix} \end{align} Btrani?Broti??=k=1∑n?mk?Φki?=k=1∑n?mk? ?ukxi?ukyi?ukzi?? ?=k=1∑m?mk?rk×?Φki?=k=1∑n?mk? ?0rkz??rky???rkz?0rkx??rky??rk?0? ? ?ukxi?ukyi?ukzi?? ???
其中， n n n 為柔性結構離散后單元總數， m k m_k mk? 為第 k k k 個單元的質量， u k u_k uk? 為第 k k k 個單元質心的位移， r k r_k rk? 為第 k k k 個單元質心的位置。上述參量中，單元的質量和質心位置可直接從有限元模型中讀取；單元的質心處位移需要在每階模態的振型結果中讀取。
采用 ANSYS 的 APDL 語言的計算方法
-
程序主要分為四部分
- 第一部分是建立結構有限元模型進行模態分析，不再贅述。
- 第二部分是采用 *get 命令提取所有單元的質量、質心坐標。
- 第三部分是在每階振型下分別計算所有單元的位移，每個單元的位移取為單元所有節點的位移平均值。
- 第四部分是把所有單元質量、質心坐標和位移代入公式(4)和式(5)計算 B t r a n \boldsymbol{B}_{tran} Btran? 和 B r o t \boldsymbol{B}_{rot} Brot? 。
-
參考APDL代碼
! ———————————————————— 計算平動耦合系數矩陣和轉動耦合系數矩陣 ———————————————————— !!!!!!!!!!!!!!!!!!!! /POST1!——————————! 提取所有單元的質量、質心坐標 !——————————! *GET, Num_Element, ELEM, ,COUNT ! 獲取模型中單元總數 *GET, cnt_element, ELEM, ,NUM,MIN ! 獲取模型中第一個單元的編號 *DIM, Amat2, ,Num_Element,5 ! 定義一個二維矩陣 *DO, NE, 1,Num_Element,1 ! 遍歷所有單元Amat2(NE, 1) = cnt_element ! 第1列存儲當前單元編號*GET, volu_e, ELEM, cnt_element, VOLU ! 獲取當前單元的體積*GET, mat_e, ELEM, cnt_element, ATTR,MAT ! 獲取材料編號dens_e = mat_mp(mat_e,1) ! 密度 !! 注意，如果涉及多種密度，需要自定義一個矩陣，按材料序號順序存放密度信息，比如這里的矩陣mat_mp，其第一列是密度信息。!!Amat2(NE, 2) = volu_e*dens_e ! 第2列存儲當前單元質量! 獲取單元質心位置*GET, cent_x, ELEM, cnt_element,CENT,X*GET, cent_y, ELEM, cnt_element,CENT,Y*GET, cent_z, ELEM, cnt_element,CENT,ZAmat2(NE, 3) = cent_x ! 第3列存儲單元質心X軸位置Amat2(NE, 4) = cent_y ! 第4列存儲單元質心Y軸位置Amat2(NE, 5) = cent_z ! 第5列存儲單元質心Z軸位置*GET, cnt_element, ELEM,cnt_element,NXTH ! 獲取下一個單元編號 *ENDDO *DMAT, TableA, D, ALLOC, Num_Element,5 *DO, i,1,Num_Element,1*DO,j,1,5,1TableA(i,j) = Amat2(i, j)*ENDDO *ENDDO *PRINT, TableA,TableA.txt!——————————! 從模態分析每一階結果中提取單元的質心位移 !——————————! *DIM, Bmat3, ,Num_Element,4,Num_Modal ! 定義一個三維矩陣 !! Num_Modal為自定義 !! *DO, NM, 1,Num_Modal,1 ! 遍歷所有模態SET,,,,,,,NM ! 從結果文件中讀取指定編號的數據集*GET, cnt_element, ELEM, ,NUM,MIN ! 獲取模型中第一個單元的編號*DO, NE, 1,Num_Element,1 ! 遍歷所有單元Bmat3(NE, 1, NM) = cnt_element ! 第1列存儲當前單元編號! 重置單元位移elem_x = 0elem_y = 0elem_z = 0ESEL,S,,,cnt_element ! 選擇指定單元編號的單元NSLE,S,ALL ! 選擇附屬在單元上的所有節點*GET, Num_Node, NODE,,COUNT ! 獲取節點總數*GET, cnt_node, NODE,,NUM,MIN ! 獲取第一個節點的編號*DO,NN, 1,Num_Node,1 ! 遍歷所有節點! 獲取節點的位移*GET, node_x, NODE, cnt_node,U,X*GET, node_y, NODE, cnt_node,U,Y*GET, node_z, NODE, cnt_node,U,Z! 計算該單元中節點位移的總和elem_x = elem_x+node_xelem_y = elem_y+node_yelem_z = elem_z+node_zcnt_node = NDNEXT(cnt_node) ! 獲取下一個節點編號*ENDDOBmat3(NE,2,NM) = elem_x/Num_Node ! 第2列存儲當前單元X軸位移Bmat3(NE,3,NM) = elem_y/Num_Node ! 第3列存儲當前單元Y軸位移Bmat3(NE,4,NM) = elem_z/Num_Node ! 第4列存儲當前單元Z軸位移ESEL,ALL ! 重新選擇所有單元*GET, cnt_element, ELEM,cnt_element,NXTH! 獲取下一個單元編號*ENDDO *ENDDO *DMAT, TableB1, D, ALLOC, Num_Element,4 *DMAT, TableB2, D, ALLOC, Num_Element,4 *DO, i,1,Num_Element,1*DO,j,1,4,1TableB1(i,j) = Bmat3(i,j,1)TableB2(i,j) = Bmat3(i,j,2)*ENDDO *ENDDO *PRINT, TableB1,TableB.txt *PRINT, TableB2,TableB.txt!——————————! 計算平動與轉動耦合系數矩陣 !——————————! *DMAT, Btran, D, ALLOC, 3,Num_Modal ! 定義平動耦合系數矩陣 *DMAT, Brot, D, ALLOC, 3,Num_Modal ! 定義轉動耦合系數矩陣 !*DIM, Btran, ,3,Num_Modal !*DIM, Brot, ,3,Num_Modal *DO, NM,1,Num_Modal,1 ! 遍歷模態! 重置平動標量tran_x = 0tran_y = 0tran_z = 0! 重置轉動標量rot_x = 0rot_y = 0rot_z = 0*DO,NE,1,Num_Element,1 ! 遍歷單元! 平動tran_x = tran_x+Amat2(NE,2)*Bmat3(NE,2,NM)tran_y = tran_y+Amat2(NE,2)*Bmat3(NE,3,NM)tran_z = tran_z+Amat2(NE,2)*Bmat3(NE,4,NM)! 臨時變量tmp_x = Amat2(NE,4)*Bmat3(NE,4,NM)-Amat2(NE,5)*Bmat3(NE,3,NM)tmp_y = Amat2(NE,5)*Bmat3(NE,2,NM)-Amat2(NE,3)*Bmat3(NE,4,NM)tmp_z = Amat2(NE,3)*Bmat3(NE,3,NM)-Amat2(NE,4)*Bmat3(NE,2,NM)! 轉動rot_x = rot_x+Amat2(NE,2)*tmp_xrot_y = rot_y+Amat2(NE,2)*tmp_yrot_z = rot_z+Amat2(NE,2)*tmp_z*ENDDOBtran(1,NM) = tran_xBtran(2,NM) = tran_yBtran(3,NM) = tran_zBrot(1,NM) = rot_xBrot(2,NM) = rot_yBrot(3,NM) = rot_z *ENDDO!——————————! 輸出平動與轉動耦合系數矩陣 !——————————! *PRINT, Btran, Btran.txt *PRINT,Brot,Brot.txtFINISH
系統轉動慣量的求解方法
-
采用apdl命令
!——————————! 續：求總質量及質心坐標 !——————————! ! mass_total = 0 sum_mx = 0 sum_my = 0 sum_mz = 0 *DO, i, 1, Num_Element, 1 mass_total = mass_total + Amat2(i, 2) ! 累加總質量 sum_mx = sum_mx + Amat2(i, 2)*Amat2(i, 3) ! 累加 X 坐標的加權值 sum_my = sum_my + Amat2(i, 2)*Amat2(i, 4) ! 累加 Y 坐標的加權值 sum_mz = sum_mz + Amat2(i, 2)*Amat2(i, 5) ! 累加 Z 坐標的加權值 *ENDDO ! 計算質心坐標 centriod_x = sum_mx/mass_total ! 計算質心 X 坐標 centriod_y = sum_my/mass_total ! 計算質心 Y 坐標 centriod_z = sum_mz/mass_total ! 計算質心 Z 坐標 ! 輸出 *DMAT, MassCentriod, D, ALLOC, 1,4 MassCentriod(1,1) = mass_total MassCentriod(1,2) = centriod_x MassCentriod(1,3) = centriod_y MassCentriod(1,4) = centriod_z *PRINT, MassCentriod, MassCentriod.txt !——————————! 續：求整個模型所有單元的總轉動慣量 !——————————! ! 相對原點 Io_xx = 0 $ Io_yy = 0 $ Io_zz = 0 Io_xy = 0 $ Io_xz = 0 $ Io_yz = 0 *DO, i, 1, Num_Element, 1 mass_e = Amat2(i, 2)cent_x = Amat2(i, 3)cent_y = Amat2(i, 4)cent_z = Amat2(i, 5)Io_xx = Io_xx + mass_e*(cent_y*cent_y+cent_z*cent_z)Io_yy = Io_yy + mass_e*(cent_x*cent_x+cent_z*cent_z)Io_zz = Io_zz + mass_e*(cent_x*cent_x+cent_y*cent_y)Io_xy = Io_xy + mass_e*(cent_x*cent_y)Io_xz = Io_xz + mass_e*(cent_x*cent_z)Io_yz = Io_yz + mass_e*(cent_y*cent_z) *ENDDO ! 相對質心 Ic_xx = 0 $ Ic_yy = 0 $ Ic_zz = 0 Ic_xy = 0 $ Ic_xz = 0 $ Ic_yz = 0 *DO, i, 1, Num_Element, 1 mass_e = Amat2(i, 2)cent_x = Amat2(i, 3) - centriod_xcent_y = Amat2(i, 4) - centriod_ycent_z = Amat2(i, 5) - centriod_zIc_xx = Ic_xx + mass_e*(cent_y*cent_y+cent_z*cent_z)Ic_yy = Ic_yy + mass_e*(cent_x*cent_x+cent_z*cent_z)Ic_zz = Ic_zz + mass_e*(cent_x*cent_x+cent_y*cent_y)Ic_xy = Ic_xy + mass_e*(cent_x*cent_y)Ic_xz = Ic_xz + mass_e*(cent_x*cent_z)Ic_yz = Ic_yz + mass_e*(cent_y*cent_z) *ENDDO ! 輸出 *DMAT, Io, D, ALLOC, 3,3 Io(1,1)=Io_xx $ Io(1,2)=-Io_xy $ Io(1,3)=-Io_xz Io(2,1)=-Io_xy $ Io(2,2)=Io_yy $ Io(2,3)=-Io_yz Io(3,1)=-Io_xz $ Io(3,2)=-Io_yz $ Io(3,3)=Io_zz *PRINT, Io, Io.txt *DMAT, Ic, D, ALLOC, 3,3 Ic(1,1)=Ic_xx $ Ic(1,2)=-Ic_xy $ Ic(1,3)=-Ic_xz Ic(2,1)=-Ic_xy $ Ic(2,2)=Ic_yy $ Ic(2,3)=-Ic_yz Ic(3,1)=-Ic_xz $ Ic(3,2)=-Ic_yz $ Ic(3,3)=Ic_zz *PRINT, Ic, Ic.txt
參考文獻
[^1]: 衛曉娜,董云峰.基于ANSYS的平動和轉動耦合系數矩陣計算[C]//北京力學會,北京振動工程學會.北京力學會第21屆學術年會暨北京振動工程學會第22屆學術年會論文集.北京航空航天大學宇航學院;,2015:758-760.
[^2]: 郭江.衛星撓性附件平動轉動耦合系數的有限元分析[D].南京理工大學,2019.DOI:10.27241/d.cnki.gnjgu.2019.002186.
)
 方法)
