基于混合有限元法求解力学问题的方法

1.一种基于混合有限元法求解力学问题的方法，其特征在于，所述方法包括如下求解步骤：步骤S1，获取划分受力体的多个单元的单元数据、材料参数、边界条件和体力，其中，单元数据至少包括：单元形状、节点坐标、单元数量、单元类型、单元维度，且单元类型包括：内蕴混合有限元类型、拉格朗日元类型，材料参数包括：杨氏模量、泊松比、密度，边界条件包括：面力、位移；步骤S2，循环每个单元，获取当前单元上的多个积分点，并结合步骤S1所获取的单元数据、材料参数和体力确定体力在各个积分点处的值以及与单元类型对应的有限元基函数在各个积分点处的值、梯度值或散度值，包括：内蕴混合有限元基函数在各个积分点处的值、内蕴混合有限元基函数在各个积分点处的散度值、间断拉格朗日元基函数在各个积分点处的值、连续拉格朗日元基函数在各个积分点处的值和连续拉格朗日元基函数在各个积分点处的梯度值；步骤S3，依据单元类型确定各个单元的有限元基函数的整体自由度，并结合步骤S2所得的体力在单元内各个积分点处的值以及与单元类型对应的有限元基函数在单元内各个积分点处的值、梯度值或散度值构造系数矩阵A和载荷向量b；步骤S4，对每个边界单元循环，根据边界单元上的边界条件的类型对每个边界条件进行处理，调整载荷向量b，同时得到约束方程组Gx=s，其中，G表示约束矩阵，x表示待求解向量，s表示约束向量；步骤S5，对拉格朗日元类型单元和内蕴混合有限元类型单元的每个单元交界面循环，获取拉格朗日元类型单元的阶次k1、内蕴混合有限元类型单元的阶次k2以及当前单元交界面上由内蕴混合有限元类型单元指向拉格朗日元类型单元的单位法向向量；步骤S6，计算k1阶拉格朗日元基函数、k2阶内蕴混合有限元基函数在当前单元交界面的值；步骤S7，结合步骤S5至S6确定当前单元交界面的局部系数矩阵C，包括：遍历每个单元交界面，执行如下步骤：步骤S71，计算步骤S5所得的单位法向向量与步骤S6所得的k2阶内蕴混合有限元基函数在当前单元交界面内各个积分点处的值的乘积；步骤S72，计算所得的乘积结果与步骤S6所得的k1阶拉格朗日元基函数在当前单元交界面内各个积分点处的值的乘积，并将所得的乘积结果作为当前交界面的局部系数矩阵C；步骤S8，基于步骤S3所得的单元的有限元基函数的整体自由度以及步骤S6确定当前单元交界面的k1阶拉格朗日元基函数的第一自由度和k2阶内蕴混合有限元基函数的第二自由度；步骤S9，基于步骤S8所得的当前单元交界面的第一自由度和第二自由度将步骤S7所得的当前单元交界面的局部系数矩阵C组装至系数矩阵A；步骤S10，利用步骤S4所得的约束方程组Gx=s，对步骤S9所得的系数矩阵A和步骤S4所得的载荷向量b进行处理，得到线性代数方程组Ax=b；步骤S11，求解线性代数方程组Ax=b，得到应力结果与位移结果。 2.根据权利要求1所述的方法，其特征在于，所述步骤S2包括如下步骤：遍历每个单元，执行如下步骤：步骤S21，对于每个当前单元，获取当前单元的单元体积、单元阶次k、多个积分点以及每个积分点的积分权重；步骤S22，利用步骤S1所获取的体力，确定体力在当前单元上各个积分点处的值；步骤S23，确定当前单元的单元类型；步骤S24，在单元类型为拉格朗日元类型的情况下，利用步骤S21计算当前单元的k阶拉格朗日元基函数在各个积分点处的值和k阶拉格朗日元基函数在各个积分点处的梯度值；步骤S25，在单元类型为内蕴混合有限元类型的情况下，获取当前单元的三维对称矩阵空间的张量基，并结合步骤S21计算所得的k阶拉格朗日元基函数在各个积分点处的值和k阶拉格朗日元基函数在各个积分点处的梯度值，得到当前单元的内蕴混合有限元基函数在各个积分点处的值和内蕴混合有限元基函数在各个积分点处的散度值，同时利用步骤S21计算当前单元的k-1阶间断拉格朗日元基函数在各个积分点处的值。 3.根据权利要求2所述的方法，其特征在于，所述步骤S25包括如下步骤：步骤S251，利用步骤S21计算k阶拉格朗日元基函数在各个积分点处的值、k阶拉格朗日元基函数在各个积分点处的梯度值、k-1阶间断拉格朗日元基函数在各个积分点处的值；步骤S252，确定当前单元上各个插值点的类型，其中，插值点的类型包括：单元顶点、单元内部点、单元边上的点以及单元面上的点；步骤S253，根据插值点的类型确定三维对称矩阵空间的张量基；步骤S254，将各类插值点对应的三维对称矩阵空间的张量基与步骤S251所得的k阶拉格朗日元基函数在各个积分点处的值和k阶拉格朗日元基函数在各个积分点处的梯度值结合，得到内蕴混合有限元基函数在各个积分点处的值和内蕴混合有限元基函数在各个积分点处的散度值。 4.根据权利要求3所述的方法，其特征在于，所述步骤S253包括如下步骤：第一步：在插值点的类型为单元顶点和单元内部点时，获取三阶对称张量空间的标准基；第二步：在插值点的类型为单元边上的点时，获取当前单元的六条边各自的一个单位切向向量和两个单位法向向量，并利用每条边的一个单位切向向量和两个单位法向向量确定该单元边的边基；第三步：在插值点的类型为单元面上的点时，获取当前单元的四个面各自的两个单位切向向量和一个单位法向向量，并利用每个面的两个单位切向向量和一个单位法向向量确定该单元面的面基。 5.根据权利要求1所述的方法，其特征在于，所述步骤S3包括如下步骤：遍历每个单元，执行如下步骤：步骤S31，确定当前单元的单元类型；步骤S32，在单元类型为拉格朗日元类型的情况下，按照如下步骤构造系数矩阵A和载荷向量b，包括：步骤S321，按照连续自由度顺序对k阶拉格朗日元基函数对应的插值点进行排序，得到当前单元的k阶拉格朗日元基函数的整体自由度；步骤S322，利用步骤S24所得的k阶拉格朗日元基函数在各个积分点处的梯度值计算单元刚度矩阵，并结合步骤S321所得的k阶拉格朗日元基函数的整体自由度，将所计算的单元刚度矩阵装至系数矩阵A；步骤S323，利用步骤S24所得的k阶拉格朗日元基函数和体力在各个积分点处的值计算单元载荷向量，并结合步骤S321所得的k阶拉格朗日元基函数的整体自由度，将所计算的单元载荷向量组装至载荷向量b；步骤S33，在单元类型为内蕴混合有限元类型的情况下，按照如下步骤构造系数矩阵A和载荷向量b，包括：步骤S331，确定当前单元的内蕴混合有限元基函数的整体自由度；步骤S332，利用步骤S25所得的内蕴混合有限元基函数在各个积分点处的值计算单元系数矩阵，并结合步骤S331构造系数矩阵A；步骤S333，按照间断自由度顺序对步骤S25所得的k-1阶间断拉格朗日元基函数对应的插值点进行排序，得到当前单元的k-1阶间断拉格朗日元基函数的整体自由度；步骤S334，利用步骤S25所得的k-1阶间断拉格朗日元基函数在各个积分点处的值、内蕴混合有限元基函数在各个积分点处的散度值计算单元散度算子矩阵，并结合步骤S331所得的内蕴混合有限元基函数的整体自由度、步骤S333所得的k-1阶间断拉格朗日元基函数的整体自由度将单元散度算子矩阵组装至系数矩阵A；步骤S335，利用步骤S25所得的k-1阶间断拉格朗日元基函数和体力在各个积分点处的值计算单元载荷向量，并结合步骤S333所得的k-1阶间断拉格朗日元基函数的整体自由度将单元载荷向量组装至载荷向量b。 6.根据权利要求5所述的方法，其特征在于，所述步骤S331包括如下步骤：第一步：确定内蕴混合有限元基函数对应的插值点的类型；第二步：在插值点的类型为单元顶点和单元内部点时，按照连续自由度顺序确定内蕴混合有限元基函数的整体自由度；第三步：在插值点的类型为单元边上的点或单元面上的点时，判断内蕴混合有限元基函数是否由纯切向向量组成；第四步：在内蕴混合有限元基函数是由纯切向向量组成的情况下，按照间断自由度顺序确定内蕴混合有限元基函数的整体自由度；第五步：在内蕴混合有限元基函数不是由纯切向向量组成的情况下，按照连续自由度顺序确定内蕴混合有限元基函数的整体自由度。 7.根据权利要求1所述的方法，其特征在于，所述步骤S4包括如下步骤：遍历每个边界单元的各个边界条件，执行如下步骤：步骤S41，确定当前边界条件所属边界单元的单元类型；步骤S42，在单元类型为内蕴混合有限元类型的情况下，按照如下步骤对当前边界条件进行处理，包括：步骤S421，在当前边界条件为位移时，获取当前边界的单位外法向向量，并调用内蕴混合有限元基函数，并将内蕴混合有限元基函数、单位外法向向量与位移的积分结果加载到载荷向量b中；步骤S422，在当前边界条件为面力时，调用最小二乘法对面力进行处理，并将所得的处理结果添加到约束方程组Gx=s中；步骤S43，在单元类型为拉格朗日元类型的情况下，按照如下步骤对当前边界条件进行处理，包括：步骤S431，在当前边界条件为面力时，调用拉格朗日元的k阶拉格朗日元基函数，并将k阶拉格朗日元基函数与面力的积分结果加载到载荷向量b中；步骤S432，在当前边界条件为位移时，调用最小二乘法对位移进行处理，并将所得的处理结果添加到约束方程组Gx=s中。 8.一种计算机程序产品，其特征在于，包括：计算机程序，其中，所述计算机程序被处理器执行时实现权利要求1至7中任意一项所述的基于混合有限元法求解力学问题的方法。 9.一种电子设备，其特征在于，包括：存储器和处理器，所述处理器用于运行存储在所述存储器中的程序，其中，所述程序运行时执行权利要求1至7中任意一项所述的基于混合有限元法求解力学问题的方法。