[发明专利]基于五次埃尔米特形函数的柔性体结构高阶非线性有限元数值模拟方法

权利要求书

1.一种基于五次埃尔米特形函数的柔性体结构高阶非线性有限元数值模拟方法，其特征在于具体步骤为：

步骤一：基于五次埃尔米特(Quintic Hermite)高阶插值理论，推导五次埃尔米特插值形函数，其包括6个五次形函数具体为：

其中，x_a和x_b分别为空间有限元单元节点的起点坐标和终点坐标，总体坐标系x下的坐标区间范围为x∈[x_a,x_b]；局部坐标系下ξ的坐标区间范围为ξ∈[0,1]；

所述五次埃尔米特插值形函数的属性方程为：

步骤二：确定柔性体结构单元和该柔性体结构单元需要求解的参数，根据柔性体结构单元尺寸参数，确定该柔性体结构单元的物理模型，建立柔性体结构单元在空间和时间的数学运动偏微分方程；并分析得到该柔性体结构单元在轴向的变形特性，获取柔性体结构单元沿着轴向的拉伸限制条件；

其中柔性体结构单元在空间和时间的数学运动偏微分方程具体为：

其中，M_m为质量矩阵；

r、q均为随弧坐标系下的对弧长s和时间t变化的变量；

其中，r为柔性体结构单元在变形后的空间瞬时位置向量；

B为柔性体结构单元的抗弯刚度，B＝EI，E为柔性体结构单元材料的杨氏模量，I为柔性体结构单元对中性轴的惯性矩；

为柔性体结构单元轴向的有效张力；

q为外界横向载荷的沿着柔性体结构单元长度方向的载荷分布；

符号顶端的黑点表示对时间的导数，两个点表示对时间的二阶导数；符号右上角的'表示对空间弧长的一阶导数；符号右上角的”表示对空间弧长的二阶导数；符号右上角的”'表示对空间弧长的三阶导数；符号右上角的””表示对空间弧长的四阶导数；

步骤三：将步骤一得到的五次埃尔米特插值形函数(1)代入柔性体结构单元在空间和时间的数学运动偏微分方程(3)，采用分部积分对整个柔性体结构单元进行积分，并结合步骤一中所述五次埃尔米特插值形函数属性方程(2)，建立柔性体结构单元的空间单元上的非线性系统方程：

其中，L表示柔性体结构单元的长度；

a_i(s)代表五次埃尔米特插值形函数在总体弧坐标系中的表达式；

a_i(s)＝J_i·φ_i(ξ)(i＝1～6)，建立了总体坐标系和单元坐标系下形函数表达式的转换关系；

J_i表示在不同坐标系下的雅克比变换系数；

下标i与柔性体结构单元的长度L之间的关系式为：

步骤四：施加柔性体结构单元两端的边界条件，将方程(5)等号右侧定为柔性体结构单元的广义力F_i，该柔性体结构单元的广义力F_i表达式为：

步骤五：运用张量运算法则，离散化步骤三所述非线性系统方程(5)中的参数：

通过五次埃尔米特插值形函数模拟得到柔性体结构单元的空间瞬时位置分布r(s，t)及高阶导数r⁽ⁿ⁾(s，t)；

通过二次插值函数p_m(s)模拟得到柔性体结构单元内部有效张力分布柔性体结构单元内部质量分布M(s，t)、柔性体结构单元内部抗弯刚度分布B(s)；

具体模拟物理量方程为：

步骤六：将步骤五得到的数值模拟物理量方程(7)带入柔性体结构单元的空间单元上的非线性系统方程(5)，并结合步骤五中柔性体结构单元的广义力F_i表达式(6)和五次埃尔米特插值形函数的属性方程(2)，推导得出柔性体结构单元内的离散方程：

α_ikm为伽辽金法产生的弯曲刚度项相关的积分系数矩阵；

β_ikm为伽辽金法产生的拉力项相关的积分系数矩阵；

γ_ikm为伽辽金法产生的惯性力项相关的积分系数矩阵；

q_in为外部横向分布载荷作用到柔性体单元产生的向量；

F_in为单元节点处的广义力向量；

步骤七：对所述柔性体结构单元内的离散方程(8)进行坐标系转换：将五次埃尔米特插值形函数的坐标系从空间总体坐标系通过雅克比转换系数J_i转换到单元局部坐标系，得到局部坐标系下的柔性体结构单元内的离散方程，用于对方程(8)中的刚度系数矩阵进行求解；

其中，局部坐标系下的柔性体结构单元内的离散系统方程(9)为：

方程(8)中的所有刚度系数矩阵求解关系式为：

步骤八：联立步骤二建立柔性体结构单元沿着轴向的拉伸限制条件，重复步骤三到步骤五，将所述柔性体结构单元沿着轴向的拉伸限制条件进行离散化，得到所述柔性体结构单元沿着轴向的拉伸限制条件对应的离散方程(11)；

步骤九：根据所述局部坐标系下的柔性体结构单元内的离散系统方程(9)和方程柔性体结构单元沿着轴向的拉伸限制条件对应的离散方程(11)，用牛顿迭代法对柔性体结构单元的空间单元上的非线性系统耦合方程进行求解。