CN105260499A - 一种三维柔性梁系统力学的运动仿真方法及系统 - Google Patents

Abstract

本发明提出一种三维柔性梁系统力学的运动仿真方法及系统。首先设定建立动力学方程的条件，然后根据设定的条件建动力学方程，将每个物体的位置和梯度作为动力学方程待求解的未知数；最后求解动力学方程，获得三维柔性梁系统中每个物体在每个时刻的位置，根据每个物体在每个时刻的位置信息构建三维柔性梁系统动画图像并显示出来。本发明可以对三维柔性梁进行动力学仿真并方便、快捷地设计三维柔性梁系统。

Description

一种三维柔性梁系统力学的运动仿真方法及系统

技术领域

本发明属于力学与机械技术领域，具体涉及一种三维柔性梁系统力学的运动仿真方法。

背景技术

从20世纪70－80年代起，传统意义上的CAD/CAE/CAM技术开始进入实用阶段，它们主要关注产品零部件质量和性能，通过采用结构设计、工程分析和制造过程控制的软件或工具，以达到设计和制造高质量零部件的目的。具体地说，传统的CAD技术基于三维实体几何造型技术，支持产品零部件的详细结构设计和形态分析。传统的CAE技术主要指应用有限元软件，完成产品零部件的结构分析、热分析、振动特性等功能分析问题。传统的CAM技术旨在提高产品零部件的可制造性，提供对机床、机器人、铸造过程、冲压过程、锻造加工等方面更好的控制。在过去的几十年里，传统的CAD/CAE/CAM技术在主要的工业领域(汽车、航空、通用机械、机械电子等)得到了广泛的应用，并且取得了巨大的成效。以汽车工业来说，在1995－1999的五年里，零部件故障率降低了40％，与之相伴的，是产品开发和制造成本的相应降低。

机械系统力学的运动仿真技术作为一种虚拟样机技术与CAD/CAE/CAM技术，向系统的设计/分析/制造、以提高产品整体质量和性能并降低开发与制造成本，提高零部件的质量和性能产生很好的作用。

现阶段国内的产品主要还是依靠于实验以及经验，对于使用仿真方法进行设计还是比较少，当前，安全、稳定、节能成为机械制造业发展的主题，为了更好的制造，对于系统进行动力学仿真是非常重要的，这有助于对于机械产品进行设计以及修改，进而进行生产销售。国内对于动力学仿真方法急需自主发展。如何进行自有仿真方法的开发，并把它投入到工业界进行学习以及发展，是当前的一个研究重点，以及急需发展的一个重要环节。

发明内容

本发明的目的在于提出一种三维柔性梁系统的动力学仿真方法，可以对三维柔性梁系统的运动进行仿真。

为了解决上述技术问题，本发明提供一种三维柔性梁系统力学的运动仿真方法，包括以下步骤：

第一步，设定建立动力学方程的条件，包括

(1)设定三维柔性梁系统中各物体的材料参数(弹性模量，密度等)、尺寸参数(长度，宽度，高度)以及初始位置参数；

(2)选择三维柔性梁系统中各个物体间的约束方式建立三维的柔性梁系统模型，约束方式包括固定约束，移动约束等；

(3)给定三维柔性梁系统模型所受的外力；

(4)根据需要设定三维柔性梁系统动力学仿真的时间，以及计算的时间步长。

第二步，构建动力学方程，将每个物体的位置和梯度作为动力学方程待求解的未知数；

第三步，求解动力学方程，获得三维柔性梁系统中每个物体在每个时刻(即每个步长)的位置；根据每个物体在每个时刻的位置信息构建三维柔性梁系统动画图像并显示出来。这样让用户有一个清晰的认识，判断该设计是否满足自己的需要，用户可以根据需要修改设计。

进一步，体三维柔性梁系统动力学方程如式(1)所示：

$\left[\begin{array}{cc}M\stackrel{\·\·}{q}& {\mathrm{\Φ}}_q^T\\ {\mathrm{\Φ}}_q& 0\end{array}\right]\left[\begin{array}{c}\stackrel{\·\·}{q}\\ \mathrm{\λ}\end{array}\right]=\left[\begin{array}{c}Q\\ \mathrm{\γ}\end{array}\right]---\left(1\right)$

式(1)中，Φ_q为Φ关于q的导数阵，q为系统的广义坐标，Φ＝0为三维柔性梁系统的约束方程，λ和γ分别为拉格朗日乘子矩阵和加速度约束方程的右项；为单体三维柔性梁系统的动力学方程；M＝diag(M₁M₂……M_N)为总的质量阵；Q＝diag(Q₁Q₂……Q_N)为总的外力阵。

进一步，使用龙格库塔法对动力学方程进行求解。

本发明还提供一种三维柔性梁系统力学的运动仿真系统，包含以下模块：

模型参数输入模块，用于输入三维柔性梁系统各物体的尺寸参数，材料参数以及各物体的初始位置参数；

约束方式选定模块，用于在预先存储的多种约束方式中选定各物体间需要的约束方式；

外力输入模块，用于输入模型所受的外力；

计算时间以及计算步长输入模块，用于输入三维柔性梁系统力学仿真的时间，以及计算的时间步长；

动力学计算模块，用于进行动力学分析，其包括动力学方程构建子模块和计算子模块；动力学方程构建子模块用于根据模型参数输入、连接方式选定、外力模块输入，建立动力学方程，将每个物体的位置和梯度作为待求解的未知数；计算子模块用于求解三维的柔性梁系统运动的力学方程，获得三维的柔性梁系统中每个物体的在每个时刻(即时间步长)的位置；

显示模块，用于根据每个时刻系统中每个物体的位置，显示出三维柔性梁系统的图像画面。

本发明与现有技术相比，其显著优点在于，本发明提供了输入环节，计算环节，以及结果显示功能，能全面地对三维柔性梁系统进行动力学模拟仿真；本发明的功能在国内属于先进水平，很好的填补国内空白。

附图说明

图1是三维梁单元广义坐标表示示意图。

图2是柔性单摆示意图。

图3是本发明仿真试验一中对柔性梁自由下落的动力学问题进行仿真的结果图，其中(a)为使用本发明方法获得的仿真结果图，(b)为使用现有技术获得仿真结果图。

图4是三维柔性梁系统自由下落示意图。

图5是本发明仿真试验二中末端点A和B的位置示意图，在梁的刚性比较大的情况，梁的弹性模量为E＝2.0e10N/m2，其中，(a)为第一根梁末端A点的位置示意图，(b)为第二根梁末端B点的位置示意图。

图6是本发明仿真试验二中对三维柔性梁系统自由下落的动力学问题进行仿真的结果示意图，其中E＝2.0e10N/m2。

图7是本发明仿真试验二中第一根梁末端A点的位置仿真结果图，其中，E＝2.0e6N/m2。

图8是本发明仿真试验二中对三维柔性梁系统自由下落的动力学问题进行仿真的结果图，其中E＝2.0e6N/m2。

图9是本发明仿真试验二中第一根梁末端A的位置仿真结果图，其中，E＝2.0e5N/m2。

图10是本发明仿真试验二中对三维柔性梁系统自由下落的动力学问题进行仿真的结果图，其中E＝2.0e5N/m2。

具体实施方式

一、基于绝对节点坐标法的三维梁单元

1.1梁单元的表达式

结合图1，在绝对节点坐标法中，三维柔性梁单元的广义坐标是由梁单元的两个节点的全局位置坐标和梯度组成的，梁单元和两个节点的表达式如式(1)和(2)所示：

e＝[e_Ae_B]^T(1)

$\begin{array}{cc}{e}_j={\left[\begin{array}{cccc}{r}^T& r_x^T& r_y^T& r_z^T\end{array}\right]}_j^T& j=A,B\end{array}---\left(2\right)$

式(1)和(2)中，e为梁单元的广义坐标，e_j分别为两节点的广义坐标。

其中(2)式中的r为节点A和B的全局位置坐标，表示梁单元的节点位置坐标，r_xr_yr_z为节点A和B的全局位置坐标对x,y,z方向的梯度，A和B表示梁单元的A点和B点，这样e_j为12×1的向量。

而梁单元中任意点P的位置坐标可以通过单元坐标以及式(3)所示形函数来表示：

r＝Se(3)

式(3)中S为形函数，S的表达式为：

S＝[S₁I₃S₂I₃S₃I₃…S₇I₃S₈I₃]

S₁＝1-3ξ²+2ξ³,S₂＝l(1-2ξ²+ξ³),

S₃＝l(η-ξη),S₄＝l(ζ-ξζ),

S₅＝3ξ²-2ξ³,S₆＝l(-ξ²+ξ³),

S₇＝lξη,S₈＝lξζ

其中，I₃为3阶单位矩阵，ξ＝x/l,η＝y/l,ζ＝z/l，l表示梁单元的长度，x,y,z表示P点在梁单元坐标中的位置。

1.2梁单元的质量阵

将式(3)对时间求导，可以得到梁单元中任一点的绝对速度公式，如式(4)所示，

$\stackrel{\·}{r}=S\stackrel{\·}{e}---\left(4\right)$

式(4)中，表示梁单元中任一点的绝对速度，表示广义坐标的导数。利用(4)式可以求得梁单元的动能，动能如式(5)所示，

$\begin{array}{c}{T}^e=\frac{1}{2}{\∫}_V\mathrm{\ρ}{\stackrel{\·}{r}}^T\stackrel{\·}{r}dV\\ =\frac{1}{2}{\stackrel{\·}{e}}^T{\∫}_V{\mathrm{\ρS}}^{T}SdV\stackrel{\·}{e}=\frac{1}{2}{\stackrel{\·}{e}}^T{m}^e\stackrel{\·}{e}\end{array}---\left(5\right)$

式(5)中，T^e为梁单元的动能，m^e＝∫_VρS^TSdV为梁单元的质量阵，ρ为材料的质量密度，V为梁单元的体积。由(5)可以看出，质量矩阵是一个常量矩阵。

1.3梁单元的刚度阵

梁单元对x,y,z方向的位移梯度J如式(6)所示，

$J=\frac{\∂r}{\∂x}=\left(\begin{array}{ccc}{S}_{1x}e& S_{1y}e& S_{1z}e\\ S_{2x}e& S_{2y}e& S_{2z}e\\ S_{3x}e& S_{3y}e& S_{3z}e\end{array}\right)---\left(6\right)$

式(6)中，x＝[xyz]^T， $\begin{array}{ccc}{S}_{ix}=\∂S_i/\∂x,& S_{iy}=\∂S_i/\∂y,& S_{iz}=\∂S_i/\∂z,\end{array}$ S_i为单元形函数的第i行。

利用位移梯度，可以获得梁单元的格林-柯西应变张量ε_G如式(7)所示，

$\begin{array}{c}{\mathrm{\ϵ}}_G=\frac{1}{2}(J^{T}J-I)\\ =\frac{1}{2}\left(\begin{array}{ccc}(e^T{a}_{1}e-1)& e^T{a}_{4}e& e^T{a}_{5}e\\ e^T{a}_{4}e& (e^T{a}_{2}e-1)& e^T{a}_{6}e\\ e^T{a}_{5}e& e^T{a}_{6}e& (e^T{a}_{3}e-1)\end{array}\right)\end{array}---\left(7\right)$

式(7)中，

$\left\{\begin{array}{c}{a}_1=S_{1x}^T{S}_{1x}+S_{2x}^T{S}_{2x}+S_{3x}^T{S}_{3x},\\ a_2=S_{1y}^T{S}_{1y}+S_{2y}^T{S}_{2y}+S_{3y}^T{S}_{3y},\\ a_3=S_{1z}^T{S}_{1z}+S_{2z}^T{S}_{2z}+S_{3z}^T{S}_{3z},\\ a_4=S_{1x}^T{S}_{1y}+S_{2x}^T{S}_{2y}+S_{3x}^T{S}_{3y},\\ a_5=S_{1x}^T{S}_{1z}+S_{2x}^T{S}_{2z}+S_{3x}^T{S}_{3z},\\ a_6=S_{1y}^T{S}_{1z}+S_{2y}^T{S}_{2z}+S_{3y}^T{S}_{3z}\end{array}\right.$

其中，a₁,a₂,..,a₆均为24×24的矩阵。

应变张量的列矢量ε的形式如式(8)所示，

ε＝[ε_G11ε_G22ε_G332ε_G122ε_G132ε_G23]^T(8)

式(8)中，ε_G11,...,ε_G23表示式子(7)中ε_G的第1行第1列,...,第2行第三列等。

根据材料的本构关系，单元应力σ如式(9)所示，

σ＝Dε(9)

式(9)中，D为弹性模量矩阵，其表达式为：

$D=\left[\begin{array}{cccccc}\mathrm{\λ}+2\mathrm{\μ}& & & & & \\ \mathrm{\λ}& \mathrm{\λ}+2\mathrm{\μ}& & symmetric& & \\ \mathrm{\λ}& \mathrm{\λ}& \mathrm{\λ}+2\mathrm{\μ}& & & \\ 0& 0& 0& \mathrm{\μ}& & \\ 0& 0& 0& 0& \mathrm{\μ}& \\ 0& 0& 0& 0& 0& \mathrm{\μ}\end{array}\right]$

其中，γ为梁的泊松比，其表达式如式(10)所示，

$\mathrm{\λ}=\frac{\mathrm{\γ}E}{(1+\mathrm{\γ})(1-2\mathrm{\γ})},\mathrm{\μ}=\frac{E}{2(1+\mathrm{\γ})}---\left(10\right)$

式(10)中，E为梁的弹性模量。

梁单元的弹性力所作的虚功δU^e的表达式如式(11)所示：

$\begin{array}{c}{\mathrm{\δU}}^e=\mathrm{\δ}(-\frac{1}{2}{\∫}_V{\mathrm{\σ}}^T\mathrm{\ϵ}dV)\\ =-\frac{1}{2}\mathrm{\δ}{\∫}_V{\mathrm{\ϵ}}^{T}D\mathrm{\ϵ}dV=-{\∫}_V{\left(D\mathrm{\ϵ}\right)}^T\mathrm{\δ}\mathrm{\ϵ}dV={\left(Q_{\mathrm{\ϵ}}^e\right)}^T\mathrm{\δ}e\end{array}---\left(11\right)$

其中，为梁的广义弹性力，其表达式如式(12)所示：

$Q_{\mathrm{\ϵ}}^e=-{\∫}_V\left[\begin{array}{cccccc}{a}_{1}e& a_{2}e& a_{3}e& 2a_{4}e& 2a_{5}e& 2a_{6}e\end{array}\right]\left(D\mathrm{\ϵ}\right)dV---\left(12\right)$

1.4单元的广义外力阵

设F_P是作用于梁单元上P点的外力，则外力所做的虚功如式(13)所示：

${\mathrm{\δP}}_F^e={\mathrm{\δr}}^T{F}_P={\mathrm{\δe}}^T{Q}_F^e---\left(13\right)$

其中，为梁单元的集中力的广义阵。

重力的虚功如式(14)所示：

${\mathrm{\δP}}_g^e={\mathrm{\δr}}^T{\∫}_V\mathrm{\ρ}gdV={\mathrm{\δe}}^T{Q}_g^e---\left(14\right)$

其中为梁单元的重力的广义力阵，g为重力加速度矢量，令总的广义外力 $Q^e=Q_F^e+Q_g^e+Q_{\mathrm{\ϵ}}^e.$

2梁系统的动力学方程

设B_e为单元节点坐标和总体节点坐标的转换阵，则有q_e＝B_eq_i，其中，q_i为梁Bi总体节点坐标列，q_e为单元的节点坐标列阵，则Bi的质量阵M_i和外力阵Q_i如式(15)所示：

$\left\{\begin{array}{c}{M}_i=\underset{j=1}{\overset{n}{\Σ}}{B}_i^T{m}_j^e{B}_i\\ Q_i=\underset{j=1}{\overset{n}{\Σ}}{B}_i^T{Q}_j^e\end{array}\right.---\left(15\right)$

设柔性梁系统由梁Bi(i＝1,…,N)组成，则单体系统的动力学方程如式(16)所示，

$M\stackrel{\·\·}{q}=Q---\left(16\right)$

式(16)中：

M＝diag(M₁M₂……M_N)为总的质量阵；

Q＝diag(Q₁Q₂……Q_N)为总的外力阵；

q为系统的广义坐标。

设Φ＝0为系统的约束方程，则多体系统的动力学方程如式(17)所示：

$\left[\begin{array}{cc}M\stackrel{\·\·}{q}& {\mathrm{\Φ}}_q^T\\ {\mathrm{\Φ}}_q& 0\end{array}\right]\left[\begin{array}{c}\stackrel{\·\·}{q}\\ \mathrm{\λ}\end{array}\right]=\left[\begin{array}{c}Q\\ \mathrm{\γ}\end{array}\right]---\left(17\right)$

式(17)中，Φ_q为Φ关于q的导数阵；λ和γ分别为拉格朗日乘子矩阵和加速度约束方程的右项。

3动力学求解与图像表示

式(17)是一组微分方程，使用龙格库塔法对式(17)进行求解，这样就可以计算得到系统中各个物体的每个时刻的位置q。

通过所得到的位置q的结果，对于某一时刻，画出图形。根据梁节点的位置，求得梁中每个点的位置，得到每个梁的图形。从而得到每个时刻系统的图形。

为了说明本发明，做了一下仿真实验。

仿真实验一

对有代表意义的柔性单摆的自由下落问题进行了计算，如图2，柔性单摆的参数为：L＝1m，宽度和高度b＝h＝0.01m，泊松比v＝0.3。弹性模量为E＝2.0e6N/m2。

本仿真实验把仿真的结果画在x-y平面上并与现有技术]进行了比较，如图3。两者的结果一致，这样验证了本发明方法的有效性。

仿真实验二

对三维柔性梁系统的自由下落的例子进行了计算，如图4，其中，第一个梁与地面，第一个梁与第二个梁采用移动约束，两个梁刚开始成直角。其中梁的参数为L＝1m，宽度和高度b＝h＝0.02m，泊松比v＝0.3。

通过改变梁的弹性模量E可以计算不同结果，仿真中E分别取为2.0e10N/m2，2.0e6N/m2直到2.0e5N/m2。计算结果如图5至图10中所示，图中ADAMS结果取自多体动力学软件MSC.ADAMS与ANSYS联合使用得到的计算结果。从图5中可以发现，在刚性较大时，两者结果一致。

而当梁的柔性比较大时，弹性模量达到2.0e6N/m2时(图7)，ADAMS计算结果出现异常。随着弹性模量的继续减小，当弹性模量达到2.0e5N/m2时(图9)，ADAMS计算失败。这个结果提醒我们这样一个事实，在使用ADAMS软件进行动力学计算时，一定要谨慎看待所得的结果，或许这个结果有较大误差，甚至是错误的结论。

图8与图10可以明显的表明本发明方法可以很好的计算三维柔性梁系统的动力学问题。

Claims (4)

1.一种三维柔性梁系统的动力学仿真方法，其特征在于，包括以下步骤：

第一步，设定建立动力学方程的条件，包括：

设定三维的柔性梁系统中各物体的材料参数、尺寸参数以及初始位置参数、选择三维柔性梁系统中各个物体间的约束方式；

给定三维柔性梁系统模型所受的外力；

根据需要设定三维的柔性梁系统动力学分析的时间，以及计算的时间步长；

第二步，构建动力学方程，将每个物体的位置和梯度作为动力学方程待求解的未知数；

第三步，求解动力学方程，获得三维柔性梁系统中每个物体在每个时刻的位置；根据每个物体在每个时刻的位置信息构建三维柔性梁系统动画图像并显示出来。

2.如权利要求1所述三维柔性梁系统力学的运动仿真方法，其特征在于，多体三维柔性梁系统动力学方程如式(1)所示：

$\left[\begin{array}{cc}M\stackrel{\·\·}{q}& {\mathrm{\Φ}}_q^T\\ {\mathrm{\Φ}}_q& 0\end{array}\right]\left[\begin{array}{c}\stackrel{\·\·}{q}\\ \mathrm{\λ}\end{array}\right]=\left[\begin{array}{c}Q\\ \mathrm{\γ}\end{array}\right]---\left(1\right)$

式(1)中，Φ_q为Φ关于q的导数阵，q为系统的广义坐标，Φ＝0为三维柔性梁系统的约束方程，λ和γ分别为拉格朗日乘子矩阵和加速度约束方程的右项；为单体三维柔性梁系统的动力学方程；M＝diag(M₁M₂……M_N)为总的质量阵；Q＝diag(Q₁Q₂……Q_N)为总的外力阵。

3.如权利要求2所述三维柔性梁系统力学的运动仿真方法，其特征在于，使用龙格库塔法对动力学方程进行求解。

4.一种三维柔性梁系统力学的运动仿真系统，其特征在于，包含以下模块：

模型参数输入模块，用于输入三维柔性梁系统各物体的尺寸参数，材料参数以及各物体的初始位置参数；

约束方式选定模块，用于在预先存储的多种约束方式中选定各物体间需要的约束方式；

外力输入模块，用于输入模型所受的外力值；

计算时间以及计算步长输入模块；

动力学计算模块，用于进行动力学分析，其包括动力学方程构建子模块和计算子模块；动力学方程构建子模块用于构建三维柔性梁系统运动的力学方程，将每个物体的位置和梯度作为待求解的未知数；计算子模块用于求解三维柔性梁系统运动的动力学方程，获得三维柔性梁系统运动中每个物体的在每个时刻的位置；

显示模块，用于根据每个时刻系统中每个物体的位置，显示出三维柔性梁系统的图像画面。