Ansys 计算刚柔耦合矩阵系数

Ansys 计算刚柔耦合系数矩阵

卫星的刚柔耦合动力学模型

• 柔性航天器的刚柔耦合动力学模型可以表示为
  $$\begin{array}{rl}\mathbf{m}\dot{\mathbf{v}}+{\mathbf{B}}_{tran}\ddot{\mathbf{\eta }}& =\mathbf{F}\\ \mathbf{J}\dot{\mathbf{\omega }}+{\mathbf{\omega }}^{\times }\mathbf{J}\mathbf{\omega }+{\mathbf{B}}_{rot}\ddot{\mathbf{\eta }}+{\mathbf{\omega }}^{\times }{\mathbf{B}}_{rot}\dot{\mathbf{\eta }}& =\mathbf{T}\\ \ddot{\mathbf{\eta }}+2\mathbf{\xi }\mathbf{\Omega }\dot{\mathbf{\eta }}+{\mathbf{\Omega }}^2\mathbf{\eta }+{\mathbf{B}}_{rot}^T\dot{\mathbf{\omega }}+{\mathbf{B}}_{tran}^T\dot{\mathbf{v}}& =0\end{array}$$
  式中：

  • $\mathbf{m}$，$J$ 分别表示柔性航天器的质量和转动惯量矩阵；
  • $\mathbf{v}$，$\mathbf{\omega }$ 分别表示航天器本体相对于惯性系的线速度和角速度矢量；
  • $\mathbf{F}$，$\mathbf{T}$ 分别表示航天受到的三轴控制力和控制力矩；
  • $\mathbf{\eta }$，$\mathbf{\xi }$，$\Omega$ 分别表示航天器的模态位移向量，柔性振动的阻尼比和频率矩阵；
  • ${\mathbf{B}}_{tran}=A\cdot {\mathbf{B}}_{tran}^b$ 为柔性附件在航天器坐标系中相对于航天器本体的平动耦合系数矩阵，$\mathbf{A}$为柔性附件坐标系到卫星坐标系的转移矩阵，${\mathbf{B}}_{tran}^b$ 为柔性附件在自身坐标系中相对于航天器本体的平动耦合系数矩阵；
  • ${\mathbf{B}}_{rot}={\mathbf{l}}_p^{\times }{\mathbf{B}}_{tran}+\mathbf{A}{\mathbf{B}}_{rot}^b$ 为柔性附件在航天器坐标系中相对于坐标原点的转动耦合矩阵，${\mathbf{l}}_p$ 为由航天器质心指向柔性附件坐标原点的向量，${}^{\times }$ 表示反对称矩阵，${\mathbf{B}}_{rot}^b$ 为柔性附件在自身坐标系中相对于坐标原点的转动耦合系数。

• ${\mathbf{B}}_{tran}$ 和 ${\mathbf{B}}_{rot}$ 的表达式如式(4)和式(5)所示
  $$\begin{array}{rl}{\mathbf{B}}_{tran}& ={\left[\begin{array}{cccccc}{\mathbf{B}}_{tran}^1& {\mathbf{B}}_{tran}^2& \cdots & {\mathbf{B}}_{tran}^i& \cdots & {\mathbf{B}}_{tran}^m\end{array}\right]}^T\\ {\mathbf{B}}_{rot}& =\left[\begin{array}{cccccc}{\mathbf{B}}_{rot}^1& {\mathbf{B}}_{rot}^2& \cdots & {\mathbf{B}}_{rot}^i& \cdots & {\mathbf{B}}_{rot}^m\end{array}\right]\end{array}$$
  其中，$m$ 为模态截断阶数；每个元素 ${\mathbf{B}}_{tran/rot}^i$ 为 $3\times 1$ 的向量。

• ${\mathbf{B}}_{tran}^i$ 和 ${\mathbf{B}}_{rot}^i$ 的表达式如式(6)和式(7)所示
  $$\begin{array}{rl}{\mathbf{B}}_{tran}^i& =\sum _{k=1}^n{m}_k{\Phi }_k^i=\sum _{k=1}^n{m}_k\left[\begin{array}{c}{u}_{kx}^i\\ u_{ky}^i\\ u_{kz}^i\end{array}\right]\\ {\mathbf{B}}_{rot}^i& =\sum _{k=1}^m{m}_k{\mathbf{r}}_k^{\times }{\Phi }_k^i=\sum _{k=1}^n{m}_k\left[\begin{array}{ccc}0& -r_{kz}& r_{ky}\\ r_{kz}& 0& -r_k\\ -r_{ky}& r_{kx}& 0\end{array}\right]\left[\begin{array}{c}{u}_{kx}^i\\ u_{ky}^i\\ u_{kz}^i\end{array}\right]\end{array}$$
  其中，$n$ 为柔性结构离散后单元总数，$m_k$ 为第 $k$ 个单元的质量，$u_k$ 为第 $k$ 个单元质心的位移，$r_k$ 为第 $k$ 个单元质心的位置。上述参量中，单元的质量和质心位置可直接从有限元模型中读取；单元的质心处位移需要在每阶模态的振型结果中读取。

采用 ANSYS 的 APDL 语言的计算方法

• 程序主要分为四部分

  • 第一部分是建立结构有限元模型进行模态分析，不再赘述。
  • 第二部分是采用 *get 命令提取所有单元的质量、质心坐标。
  • 第三部分是在每阶振型下分别计算所有单元的位移，每个单元的位移取为单元所有节点的位移平均值。
  • 第四部分是把所有单元质量、质心坐标和位移代入公式(4)和式(5)计算${\mathbf{B}}_{tran}$ 和 ${\mathbf{B}}_{rot}$ 。

• 参考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.