Ansys 计算刚柔耦合系数矩阵
卫星的刚柔耦合动力学模型
柔性航天器的刚柔耦合动力学模型可以表示为
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ω˙+BtranTv˙=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} BtranBrot=[Btran1Btran2⋯Btrani⋯Btranm]T=[Brot1Brot2⋯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} BtraniBroti=k=1∑nmkΦki=k=1∑nmk ukxiukyiukzi =k=1∑mmkrk×Φki=k=1∑nmk 0rkz−rky−rkz0rkxrky−rk0 ukxiukyiukzi
其中， 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,Z Amat2(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,1 TableA(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 = 0 elem_y = 0 elem_z = 0 ESEL,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_x elem_y = elem_y+node_y elem_z = elem_z+node_z cnt_node = NDNEXT(cnt_node) ! 获取下一个节点编号 *ENDDO Bmat3(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,1 TableB1(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 = 0 tran_y = 0 tran_z = 0 ! 重置转动标量 rot_x = 0 rot_y = 0 rot_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_x rot_y = rot_y+Amat2(NE,2)*tmp_y rot_z = rot_z+Amat2(NE,2)*tmp_z *ENDDO Btran(1,NM) = tran_x Btran(2,NM) = tran_y Btran(3,NM) = tran_z Brot(1,NM) = rot_x Brot(2,NM) = rot_y Brot(3,NM) = rot_z *ENDDO !——————————! 输出平动与转动耦合系数矩阵 !——————————! *PRINT, Btran, Btran.txt *PRINT,Brot,Brot.txt FINISH
系统转动惯量的求解方法
采用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_x cent_y = Amat2(i, 4) - centriod_y cent_z = Amat2(i, 5) - centriod_z Ic_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.