本发明涉及柔性体变形场离散方法技术领域，特别是一种基于B样条计算柔性矩形薄板刚柔耦合动力学响应的方法。
背景技术：
柔性体的变形场离散问题是柔性多体系统动力学中的基本问题。目前假设模态法和有限元法是两种使用较为广泛的变形场离散方法。但当工程结构并非规则的梁式或板式结构的时候，假设模态法的局限性在于对于复杂的柔性体结构很难求出振型函数，而且当工程模型发生改变的时候，为了满足精度要求需要叠加更多的模态，这样降低了计算效率。有限元法的缺点是精度分析需要耗费大量计算资源。因此对于多体系统动力学问题人们仍需关注和探索新的变形场离散方法。
文献“作大范围运动柔性梁和柔性薄板刚柔耦合动力学建模与仿真[D].南京：南京理工大学，2009.”采用假设模态法计算柔性矩形薄板的刚柔耦合动力学响应。该方法计算的柔性矩形薄板外侧角点横向变形较小，而现代工业生产中的柔性机械往往具有较大的变形，若仍然采用该方法计算精度偏低，因此需要寻找适用于计算柔性体大变形动力学响应的变形场离散方法。
技术实现要素：
本发明提出了一种基于B样条计算柔性矩形薄板刚柔耦合动力学响应的方法。
实现本发明的技术解决方案为：一种基于B样条计算柔性矩形薄板刚柔耦合动力学响应的方法，具体步骤为：
步骤1、建立薄板的动能和势能的运动学表达式；
步骤2、采用B样条确定薄板的关于w1、w2、u3的离散化变形场函数；
步骤3、将薄板的关于w1、w2、u3的离散化变形场函数代回到薄板的动能和势能表达式，并根据拉格朗日方程建立薄板的刚柔耦合动力学方程。
本发明与现有技术相比，其显著优点为：
本方法在大变形柔性体刚柔耦合动力学响应问题的计算精度上具有优势，具体为：通过减小薄板弹性模量获得了较大的外侧角点的横向变形，采用B样条法计算结果与有限元法结果吻合度高，而假设模态法的结果误差较大，证明B样条法精度高于假设模态法。本方法丰富和拓展了柔性体变形场离散方法。
下面结合附图对本发明做进一步详细的描述。
附图说明
图1是建立薄板运动学模型所依据的矩形薄板模型图。图中，坐标系O-XYZ为惯性坐标系，oxyz为连体坐标系。连体坐标系的xoy平面与未变形前板的中面重合。板的长度为a，宽度为b，厚度为h，密度为ρ，弹性模型为E，泊松比为μ。变形前板中面上一点P0变形后至P点。
图2是实施例1中仿真的绕y轴作定轴转动的悬臂薄板模型图。
图3是当Ω取0.8Hz时薄板外侧中点z方向变形量(面外或横向变形)在B样条法、假设模态法以及有限元法三种离散方法下的计算结果图。
图4是当Ω取0.8Hz且薄板弹性模量减小至3×1010N/m2时薄板外侧中点z方向变形量在B样条法、假设模态法以及有限元法三种离散方法下的计算结果图。
图5是基于MATLAB程序建立的图形用户界面的初始界面图。
图6是输入相关参量按下“画图”按钮后得到的图形用户界面图。
具体实施方式
一种基于B样条计算柔性矩形薄板刚柔耦合动力学响应的方法，具体步骤为：
步骤1、建立薄板的动能和势能的运动学表达式，具体方法为：
步骤1.1、建立变形后薄板内任意一点P在惯性坐标系下的速度矢量，所述速度矢量表示为：
VP＝Vo+ωA×(ρ0+u)+VPA (1)
式(1)中，Vo为连体坐标系相对于惯性坐标系的速度，ωA为连体坐标系相对于惯性坐标系的角速度矢量，ρ0为薄板内点P变形前对应的点在连体坐标系中的位置矢量，u为点P变形前对应的点在连体坐标系中的变形矢量，VPA为P点相对连体坐标系的速度矢量；其中，各矢量在连体坐标系的分量形式分别为：
Vo＝v1a1+v2a2+v3a3 (2)
ωA＝ω1a1+ω2a2+ω3a3 (3)
ρ0＝xa1+ya2 (4)
u＝u1a1+u2a2+u3a3(5)
式(1)-(6)中，a1、a2、a3为连体坐标系的单位方向矢量，v1、v2、v3为Vo沿a1、a2、a3方向的分量，ω1、ω2、ω3为ωA沿α1、α2、α3方向的分量，x、y为ρ0沿a1、a2、a3方向的分量，u1、u2、u3为u沿a1、a2、a3方向的分量，为VPA沿a1、a2、a3方向的分量；
将式(2)-(6)分量形式代回式(1)中，得到
步骤1.2、将任意一点P在惯性坐标系下的速度矢量代入薄板的动能公式中得到薄板动能的运动学表达式，其中，薄板的动能公式为：
式(8)中，ρ为薄板的密度，h为薄板厚度；
薄板动能的运动学表达式为：
式(9)中，薄板任意一点P纵向变形位移u1、u2具体为：
式(10)中，w1、w2和u3为P0点沿x、y和z方向的总伸长量，第二项和为耦合变形量。
步骤1.3、根据Love-Kirchhoff假设和Von-Karman原理建立薄板内任意一点P应变位移和应力应变关系式，其中，应变位移关系式为：
式(11)中，εx、εy、γxy表示薄板单元体中应变分量；
应力应变关系式为：
式(12)中，σx、σy、τxy表示薄板单元体中应力分量，E为薄板弹性模量，μ为泊松比。
步骤1.4、将式(11)和(12)代入薄板的势能公式中得到薄板势能的运动学表达式为：
其中：
式(13)-式(15)中，Um为板面内的变形能，Ub为板弯曲变形能。
步骤2、采用B样条确定薄板的关于w1、w2、u3的离散化变形场函数，具体步骤为：
步骤2.1、将薄板上分别沿x方向和y方向作均匀划分，具体为
0＝x0<x1<x2<…<xN＝a (16)
xi＝x0+ihx，hx＝xi+1-xi＝a/N (17)
0＝y0<y1<y2<…<yM＝b (18)
yj＝y0+jhy，hy＝yj+1-yj＝b/M (19)
式(16)-式(19)中，x0、x1…xi…xN分别为划分点的横坐标，y0、y1…yj…yM分别为划分点的纵坐标，hx为x方向相邻两划分点之间距离，hy为y方向相邻两划分点之间距离。a为薄板长度，b为薄板宽度，N和M分别为x方向和y方向划分的数量。
步骤2.2、将薄板上一点沿x、y和z方向总伸长量w1、w2和u3描述为：
式(20)中，表示矩阵的Kronecker乘积,φi(x)及ψj(y)分别为两组三次B样条基函数，aij、bij、cij为待定系数，式(20)中各项分别为：
[φ]＝[φ-1φ0φ1φ2…φN+1] (21)
[ψ]＝[ψ-1ψ0ψ1ψ2…ψM+1] (22)
{a}j＝[a-1,ja0,ja1,j…aN+1,j]T (24)
{b}j＝[b-1,jb0,jb1,j…bN+1,j]T (26)
{c}j＝[c-1,jc0,jc1,j…cN+1,j]T (28)
其中，x方向函数φi(x)为符合悬臂梁边界条件的三次B样条基函数，具体为：
y方向函数ψj(y)为符合自由梁边界条件的三次B样条基函数，具体为：
其中，三次样条函数具体为：
步骤2.3、令
又令Φ1(x,y)＝Φ2(x,y)＝Φ3(x,y)＝Φ(x,y)(33)
则薄板的关于w1、w2、u3的离散化变形场函数为：
其中，Φ1(x,y)∈R1×(M+3)(N+1)、Φ2(x,y)∈R1×(M+3)(N+1)、Φ3(x,y)∈R1×(M+3)(N+1)分别为薄板纵向和横向变形的B样条插值函数的行矢量，q1(t)∈R(M+3)(N+1)、q2(t)∈R(M+3)(N+1)、q3(t)∈R(M+3)(N+1)分别为纵向和横向变形的时间变化量的列矢量。
步骤3、将薄板的关于w1、w2、u3的离散化变形场函数代回到薄板的动能和势能表达式，并根据拉格朗日方程建立薄板的刚柔耦合动力学方程，具体步骤为：
步骤3.1、将薄板的关于w1、w2、u3的离散化变形场函数代入到薄板任意一点纵向变形位移表达式u1、u2中，重新得到薄板变形量u1、u2为：
式(35)中，H1(x,y)、H2(x,y)为耦合形函数，具体为：
式(36)下标中“，”表示对坐标求偏导；
步骤3.2、将步骤3.1中的薄板变形量u1、u2代入到步骤1中的薄板动能的运动学表达式和薄板势能的运动学表达式，并取广义坐标将代入后的薄板动能的运动学表达式和薄板势能的运动学表达式带入拉格朗日方程建立系统刚柔耦合动力学方程，具体为：
式(37)中，M11＝W11，M22＝W22，M33＝W33 (38)
式中,a01、a02、a03为基点加速度在连体坐标系中的投影，相关常数阵为
Ci＝∫∫Aρh·HidA(i＝1,2) (55)
D1i＝∫∫Aρh·x·HidA(i＝1,2) (56)
D2i＝∫∫Aρh·y·HidA(i＝1,2) (57)
S1i＝∫∫Aρh·x·ΦidA(i＝1,2,3) (58)
S2i＝∫∫Aρh·y·ΦidA(i＝1,2,3) (59)
Yi＝∫∫Aρh·ΦidA(i＝1,2,3) (60)
实施例1
图2为绕y轴作定轴转动的悬臂薄板模型图。本实施例中薄板的长度a＝1.8288m，宽度b＝1.2192m，厚度h＝0.00254m，密度ρ＝2×103kg/m3，弹性模量E＝7×1010N/m2，泊松比μ＝0.3。基频为0.75Hz(4.71rad/s)。
图2中坐标原点O的加速度为0，即a01＝a02＝a03＝0,ω1＝ω3＝0,ω2＝ω,这里忽略柔性薄板的纵向变形，仅考虑横向变形，则式(37)中关于q1、q2项删去，因此得到系统的动力学方程为
式(65)中
K33＝Kf33-ω2W33+ω2D11 (66)
图2所示悬臂薄板转动规律为：
其中，T＝30s,总仿真时间t＝40s。
图3表示式(68)中的Ω取0.8Hz时薄板外侧中点z方向变形量(面外或横向变形)在B样条法、假设模态法以及有限元法三种离散方法下的结果。三种离散方法的模型均为考虑耦合变形量的一次近似模型。由图3可知，B样条法的仿真结果同假设模态法和有限元法的仿真结果吻合度很高，证明了本发明的正确性。
为了使薄板外侧角点具有较大横向变形，于是将薄板弹性模量E减小至3×1010N/m2，薄板的其他参数保持不变。薄板的旋转规律仍然采用式(68)的形式，仿真中取T＝30s,总仿真时间t＝40s。图4表示式(68)中的Ω取0.8Hz时薄板外侧中点z方向变形量(面外或横向变形)在B样条法、假设模态法以及有限元法三种离散方法下的结果。三种离散方法的模型均为考虑耦合变形量的一次近似模型。由图4可知，B样条法的仿真结果同有限元法的仿真结果吻合度很高，而假设模态法仿真结果偏小，误差较大。证明当薄板外侧角点横向变形量较大时B样条法比假设模态法的精度更有优势。
本发明基于MATLAB程序建立了图形用户界面，初始界面如图5所示。图中可输入的参数分别为：薄板长度a(m)，薄板宽度b(m)，薄板厚度h(m)，薄板弹性模量E(N/m2)和角速度Ω(rad/s)。通过输入以上参量然后点击“画图”按钮便可以在GUI界面中显示采用B样条法计算得到的薄板外侧中点z方向变形量(面外或横向变形)的图像。图6为输入一组参量后点击“画图”按钮得到的薄板外侧中点z方向变形量(面外或横向变形)随时间的变化曲线图。
由前述可知，本发明具有如下效果：基于B样条研究柔性矩形薄板刚柔耦合动力学响应的方法，创新地将B样条用来描述柔性矩形薄板的变形场，建立柔性矩形薄板的刚柔耦合动力学方程，通过算例进行仿真，将仿真结果与假设模态法及有限元法的结果进行对比表明B样条在计算柔性矩形薄板刚柔耦合动力学响应问题上的正确性，且当薄板横向变形较大时B样条法的计算精度高于假设模态法，一定程度上丰富和拓展了柔性体变形场离散方法。此外，本发明基于MATLAB程序建立图形用户界面，使得本领域研究人员在进行动力学分析时能够方便地通过修改参数进行仿真计算，同时计算的结果曲线图也能直观地展现在图形用户界面中。