有效
位移特征值的后处理方法
胡俊、陈伟、徐士雷
北京大学
摘要
本申请公开了一种位移特征值的后处理方法。包括:利用结构体的多个单元的单元数据、材料参数和边界条件构造矩阵特征值问题,求解该问题得到有限元类型对应的特征值和特征向量,并利用有限元类型对应的有限元基函数确定对应的离散特征函数;基于离散特征函数的k阶导函数在各个单元内各个积分点处的值确定重构离散k阶导函数;利用有限元类型对应的特征值、特征函数及重构离散k阶导函数确定有限元类型对应的后验误差估计子;按照后处理方法和有限元类型对应的后验误差估计子对特征值进行后处理,得到最终的特征值。本申请解决了相关技术在利用非协调有限元处理结构体的振动问题时,难以得到高精度特征值的技术问题。
1.一种位移特征值的后处理方法,其特征在于,所述方法包括如下步骤:步骤S1,获取划分结构体所得的多个单元的单元数据、材料参数、边界条件以及后处理方法,其中,结构体包括:薄膜或薄板,单元数据包括:单元形状、节点坐标、单元数量、单元维度,边界条件包括:力约束、位移约束,材料参数包括:与薄膜相关的第一材料参数、与薄板相关的第二材料参数,且第一材料参数包括:表面张力系数、单位质量密度,第二材料参数包括:抗弯刚度、单位质量密度,后处理方法包括:重构特征值或组合特征值;步骤S2,依据材料参数确定求解结构体对应的特征值问题所使用的有限元类型和求导阶次k,其中,有限元类型的类型包括:与薄膜对应的Crouzeix-Raviart元或增广Crouzeix-Raviart元、协调线性元、与薄板对应的Morley元,求导阶次k包括:薄膜对应的求导阶次k1或薄板对应的求导阶次k2;步骤S3,对结构体的每个单元循环,获取当前单元的单元面积、各个积分点以及每个积分点的积分权重,并结合步骤S1所得的单元数据与步骤S2所得的有限元类型和求导阶次k确定与有限元类型对应的有限元基函数及有限元基函数的k阶导数在当前单元内各个积分点处的值;步骤S4,利用步骤S3所得的与有限元类型对应的有限元基函数及有限元基函数的k阶导数在当前单元内各个积分点处的值确定当前单元的单元刚度矩阵和单元质量矩阵;步骤S5,确定有限元的自由度,并结合由步骤S3至S4所得的各个单元的单元刚度矩阵和单元质量矩阵,得到整体刚度矩阵A和整体质量矩阵M;步骤S6,根据边界条件的类型对步骤S5所得的整体刚度矩阵A和整体质量矩阵M进行处理,得到对应的矩阵特征值问题 ;步骤S7,求解矩阵特征值问题 ,得到与有限元类型对应的离散特征值 和离散特征向量x,并结合步骤S3所得的与有限元类型对应的有限元基函数确定与有限元类型对应的离散特征函数;步骤S8,确定与有限元类型对应的离散特征函数的k阶导函数,并利用与有限元类型对应的离散特征函数的k阶导函数在各个单元内各个积分点处的值,确定属于Crouzeix-Raviart元离散空间的重构离散k阶导函数;步骤S9,利用与有限元类型对应的离散特征值 、离散特征函数和重构离散k阶导函数确定与有限元类型对应的后验误差估计子;步骤S10,按照步骤S1所确定的后处理方法,并利用步骤S9所得的与有限元类型对应的后验误差估计子对相应有限元类型对应的离散特征值 进行后处理,得到最终的特征值。
2.根据权利要求1所述的方法,其特征在于,所述步骤S4包括如下步骤:遍历每个单元,执行如下步骤:步骤S41,利用步骤S3所得的与有限元类型对应的有限元基函数在当前单元内各个积分点处的值和步骤S1所得的材料系数确定单元质量矩阵;步骤S42,利用步骤S3所得的与有限元类型对应的有限元基函数的k阶导数在当前单元内各个积分点处的值和步骤S1所得的材料系数确定单元刚度矩阵。
3.根据权利要求1所述的方法,其特征在于,所述步骤S5包括如下步骤:步骤S51,确定有限元类型;步骤S52,在有限元类型为Crouzeix-Raviart元的情况下,将函数在单元内每条单元边上的积分平均值作为Crouzeix-Raviart元的自由度;步骤S53,在有限元类型为增广Crouzeix-Raviart元的情况下,将函数在单元内每条单元边上的积分平均值以及函数在单元上的积分平均值作为增广Crouzeix-Raviart元的自由度;步骤S54,在有限元类型为协调线性元的情况下,将函数在单元内每个单元顶点处的函数值作为协调线性元的自由度;步骤S55,在有限元类型为Morley元的情况下,将函数在单元内每个单元顶点处的函数值及函数在单元内每条单元边的法向导数的积分平均值作为Morley元的自由度;步骤S56,依据上述步骤S52至S55所得的各类有限元的自由度,对各个单元的单元刚度矩阵和单元质量矩阵进行组装,得到整体刚度矩阵A和整体质量矩阵M。
4.根据权利要求1所述的方法,其特征在于,所述步骤S6包括如下步骤:遍历每个边界条件,执行如下步骤:步骤S61,确定当前边界条件的类型;步骤S62,在所述边界条件的类型为位移约束时,按照如下规则对整体刚度矩阵A和整体质量矩阵M进行处理:步骤S621,确定当前边界条件所属单元对应的有限元类型;步骤S622,在有限元类型为Crouzeix-Raviart元或增广Crouzeix-Raviart元的情况下,将对应于边界边的所有自由度作为边界自由度,删除步骤S5所得的整体刚度矩阵A和整体质量矩阵M中对应于边界自由度的行和列;步骤S623,在有限元类型为Morley元的情况下,将对应于边界边和边界边上的顶点的所有自由度作为边界自由度,删除整体刚度矩阵A和整体质量矩阵M中对应于边界自由度的行和列;步骤S63,在所述边界条件的类型为力约束的情况下,保持整体刚度矩阵A和整体质量矩阵M不变。
5.根据权利要求1所述的方法,其特征在于,所述步骤S7包括如下步骤:步骤S71,求解矩阵特征值问题 ,得到与有限元类型对应的离散特征值 和离散特征向量x;步骤S72,循环有限元类型对应的每个自由度,确定自由度的类型;步骤S73,在自由度的类型为边界自由度的情况下,利用对应的边界条件确定离散特征函数在该自由度上的值;步骤S74,在自由度的类型为非边界自由度的情况下,利用离散特征向量x确定离散特征函数在该自由度上的值;步骤S75,利用步骤S73和步骤S74所得的离散特征函数在自由度上的值以及有限元类型对应的有限元基函数确定与有限元类型对应的离散特征函数。
6.根据权利要求1所述的方法,其特征在于,所述步骤S8包括如下步骤:步骤S81,确定与有限元类型对应的离散特征函数的k阶导函数;步骤S82,遍历单元的每条单元边,确定单元边的类型,其中,单元边的类型包括:公共边、边界边;步骤S83,在单元边的类型为公共边的情况下,在相邻两个单元中求解离散特征函数在单元边的中点处的k阶导数,并将两个k阶导数的平均值作为重构离散k阶导函数在单元边对应的Crouzeix-Raviart元自由度的值;步骤S84,在单元边的类型为边界边的情况下,从与单元边共顶点的其他单元边的多个中点内确定与单元边中点共线的两个中点,并利用这两个中点的k阶导数的平均值进行外插,得到重构离散k阶导函数在单元边对应的Crouzeix-Raviart元自由度的值;步骤S85,结合步骤S83和步骤S84所得的重构离散k阶导函数在各个单元边对应的Crouzeix-Raviart元自由度的值,得到重构离散k阶导函数。
7.根据权利要求1所述的方法,其特征在于,所述步骤S9包括如下步骤:步骤S91,计算非协调有限元对应的重构离散k阶导函数在当前单元内各个积分点处的值与非协调有限元对应的离散特征函数的k阶导数在当前单元内各个积分点处的值的差值,并利用该差值的平方以及各个积分点对应的积分权重确定在结构体所在区域上的第一积分值,其中,非协调有限元包括:Crouzeix-Raviart元、增广Crouzeix-Raviart元、Morley元;步骤S92,在采用Crouzeix-Raviart元或者增广Crouzeix-Raviart元求解薄膜对应的特征值问题的情况下,遍历每个单元,并按照如下步骤确定Crouzeix-Raviart元或者增广Crouzeix-Raviart元对应的第一类后验误差估计子:步骤S921,根据有限元类型确定仅依赖离散特征函数的k+1阶导数的第一类多项式;步骤S922,利用Crouzeix-Raviart元或者增广Crouzeix-Raviart元对应的重构离散k阶导函数计算当前单元上的重构离散k+1阶导数值,并利用所得的重构离散k+1阶导数值计算第一类多项式在当前单元上的值;步骤S923,利用第一类多项式在当前单元上的值、Crouzeix-Raviart元或者增广Crouzeix-Raviart元对应的离散特征函数在当前单元上各个积分点处的值以及各个积分点对应的积分权重确定第二积分值;步骤S924,计算Crouzeix-Raviart元或者增广Crouzeix-Raviart元对应的离散特征值与所有单元的第二积分值之和的乘积,并利用第一积分值减去乘积的二倍,得到Crouzeix-Raviart元或者增广Crouzeix-Raviart元对应的第一类后验误差估计子;步骤S93,在采用Crouzeix-Raviart元或增广Crouzeix-Raviart元求解薄膜对应的特征值问题的情况下,遍历每个单元,并按照如下方法确定Crouzeix-Raviart元或增广Crouzeix-Raviart元对应的第二类后验误差估计子:步骤S931,根据协调线性元确定仅依赖离散特征函数的k+1阶导数的第二类多项式;步骤S932,利用Crouzeix-Raviart元或增广Crouzeix-Raviart元对应的重构离散k阶导函数计算当前单元上重构离散k+1阶导数值,并利用所得的重构离散k+1阶导数值计算第二类多项式在当前单元的值、第二类多项式在当前单元的每条单元边上的积分点处的值;步骤S933,遍历当前单元上的每条单元边,确定第二类多项式在当前单元边所属的相邻两个单元的平均值,并计算该平均值与Crouzeix-Raviart元或增广Crouzeix-Raviart元对应的离散特征函数在当前单元边上的法向导数的跳跃的乘积,利用所得乘积结果、各个积分点对应的积分权重计算当前单元边的第三积分值;步骤S934,计算Crouzeix-Raviart元或增广Crouzeix-Raviart元对应的离散特征函数与离散特征值的乘积,并将该乘积与对应的离散特征函数的拉普拉斯相加,计算所得相加结果与第二类多项式的乘积在当前单元的值,结合当前单元内各个积分点对应的积分权重,得到第四积分值;步骤S935,将第一积分值加上所有单元边的第三积分值之和的二倍,再减去所有单元的第四积分值之和的二倍,得到Crouzeix-Raviart元或增广Crouzeix-Raviart元对应的第二类后验误差估计子;步骤S94,在采用Morley元求解薄板对应的特征值问题的情况下,基于步骤S91中的第一积分值确定Morley元对应的第三类后验误差估计子。
8.根据权利要求1所述的方法,其特征在于,所述步骤S10包括如下步骤:步骤S101,在步骤S1所确定的后处理方法为重构特征值的情况下,将非协调有限元类型对应的离散特征值和后验误差估计子相加,得到最终的特征值;步骤S102,在步骤S1所确定的后处理方法为组合特征值的情况下,按照预设的组合系数将任意两种有限元类型对应的离散特征值进行线性组合,得到最终的特征值,其中,组合系数由有限元对应的后验误差估计子确定。
9.一种非易失性存储介质,其特征在于,所述非易失性存储介质中存储有计算机程序,其中,所述非易失性存储介质所在设备通过运行所述计算机程序执行权利要求1至8中任意一项所述的位移特征值的后处理方法。
10.一种计算机程序产品,其特征在于,包括:计算机程序,其中,所述计算机程序被处理器执行时实现权利要求1至8中任意一项所述的位移特征值的后处理方法。




