本发明属于三维力学振动分析数值求解技术领域，涉及一种基于三层预处理子的三维力学模态仿真模拟方法。
背景技术：
飞行器在受到各种激励作用时会发生振颤，例如亚音速导致全机振动；低超音速时引起方向舵和垂尾振动；起落架放下时的整机振动等等。在上述激励的作用下飞行器会受力失衡、材料失效、甚至结构解体，这将不利于飞行器的气动一体化设计。解决或减小机身结构的各类振动问题，往往不是单纯增加结构刚度、强度所能凑效的，有时还可能适得其反，因为结构振动问题都与机身及其局部结构的固有频率和振型有关，所以进行飞行器的振动模态分析是非常必要的。模态分析可以获得飞行器的振动特性，是其结构设计的重要环节，主要校核结构承受各种载荷作用时，能否避免共振。因此快速有效的模态分析技术是飞行器设计的关键。一般实验测试方法由于其成本高昂，周期长难以大规模应用，而数值模拟方法由于其成本低、周期短、便于实现得到大规模应用。
目前，利用各种力学数值计算方法对飞行器结构进行模态仿真分析时，都是采用的有限元本征分析方法，而且大部分的结构分析商业软件也采用的是有限元方法，例如ansys、abaqus等。有限元分析一般包括，建模、网格划分、单元计算、矩阵集成、引入边界和激励、求解广义本征方程、后处理等几个步骤。在上述过程中，求解广义本征方程是最消耗计算时间和内存的步骤，因此寻求一种高效的求解技术是关键。
目前广泛采用的是一种隐式重启arnoldi迭代方法来求解，在这一过程中需要反复求解一个线性方程。预处理共轭梯度法广泛应用于线性方程的求解，而且对于本征问题，由于每次本征迭代时线性方程系数矩阵不变，因此只需要做一次预处理，这样就大大提高了计算效率。然而预处理共轭梯度法的思想是通过预处理子把原来的线性方程变换到一个容易求解的新线性方程，这说明一个好的预处理子对整个arnoldi过程的重要性。特别是当飞行器结构复杂时，网格划分较多，线性方程往往达到千万阶，一般的预处理子在计算时间和内存消耗上往往十分巨大，这将无法快速完成飞行器结构的模态分析。因此需要构造一个高效的预处理子来提高预处理共轭梯度法的求解效率，从而提高arnoldi本征迭代的速度，最终实现飞行器结构的快速模态分析。
技术实现要素：
针对上述存在问题或不足，为解决构造有限元模态分析中的高效预处理子，实现飞行器结构的快速模态分析；本发明提供了一种基于三层预处理子的三维力学模态仿真模拟方法。
该基于三层预处理子的三维力学模态仿真模拟方法，包括以下步骤：
a.将目标飞行器结构进行有限元建模，引入位移边界条件或者应力边界条件建立对应的几何结构模型；
b.采用四面体网格剖分求解域；
c.选择二阶标量叠层基函数，并运用标准的有限元本征分析方法得到目标飞行器结构的有限元本征方程；
d.采用arpack软件包中的隐式重启arnoldi迭代方法求解c步骤所获得的有限元本征方程，获得特征值λ和对应的特征向量α即振幅向量，arnoldi迭代方法每一步迭代所产生的线性方程采用基于三层预处理子的预处理共轭梯度法求解，三层预处理子的具体构造方法如下：
在隐式重启arnoldi迭代的过程中需要反复求解以下线性方程
px＝myj(1)
其中矩阵p为待预处理的矩阵，m为系统质量矩阵，x为待求中间向量，yj为迭代向量，采用预处理共轭梯度方法求解(1)式，在预处理共轭梯度法中需要求解如下预处理系统
nz＝r(2)
其中n为预处理子，z待求预处理向量，r为残差向量；
首先把矩阵p按照x,y,z三个方向分块排序，获得如下分块矩阵形式
其中pxx为x方向的矩阵，pyy为y方向的矩阵，pzz为z方向的矩阵，pxy和pyx为x方向和y方向的耦合矩阵，pxz和pzx为x方向和z方向的耦合矩阵，pyz和pzy为y方向和z方向的耦合矩阵，取第一层预处理子的形式为
对(4)式中的块对角矩阵pxx，pyy，pzz采用基于叠层基函数的性质来构造第二层预处理子：首先构造pxx的第二层预处理子，按照一阶基函数、二阶基函数的排列顺序把pxx写成如下分块形式
其中p11为一阶基函数形成的矩阵，p22为二阶基函数形成的矩阵，p12和p21为一阶基函数和二阶基函数之间的耦合矩阵，取第二层预处理子为如下形式
矩阵pyy，pzz所对应的第二层预处理子的构造方法和矩阵pxx的第二层预处理子的构造方法是一样的，分别给他们所对应的第二层预处理子取名为ny，nz。
对于(6)式中的块对角矩阵p11和p22采用多波前不完全cholesky分解来构造第三层预处理子，矩阵p11和p22所对应的第三层预处理子分别为nl和nh：首先对矩阵p11进行多波前不完全cholesky分解，获得如下第三层预处理子的形式
nl＝(d0p0)-1llt(p0td0t)-1.(7)
其中d0为对角矩阵，p0为重排序矩阵，l为矩阵p11的cholesky分解下三角矩阵，nh的构造方法和nl的构造方法一样，完成三层预处理子的构造；
e.对d步骤获得的特征值和对应特征向量进行后处理获得振动模态频率和对应振动振型。
综上所述，本发明构造的三层预处理子的三维力学模态仿真模拟方法能大幅提高arnoldi迭代方法的求解速度，实现快速的模态分析。
附图说明
图1是本发明基于三层预处理子的三维力学模态仿真模拟方法的流程图；
图2是实施例的有限元模型图；
图3是实施例与ansys软件的计算性能对比图。
具体实施方式
下面结合附图和实施例来详细说明本发明的技术方案。
参照图1，一种基于三层预处理子的三维力学模态仿真模拟方法，包括以下步骤：
a.将目标飞行器结构进行有限元建模，引入位移边界条件或者应力边界条件建立对应的几何结构模型。
根据飞行器的特性，引入位移边界条件建立对应的几何结构模型仿真整个结构的自由振动特性，如图2所示。
b.采用四面体网格剖分求解域。
剖分后的求解域被人为分割为多个三维四面体网格，从而将连续的几何结构空间转化为离散的网格空间。
c.选择二阶标量叠层基函数，并运用标准的有限元本征分析方法得到目标飞行器结构的有限元本征方程。
最终的有限元本征分析方程
(-ω2m+k)α＝0(8)
其中m为质量矩阵、k为刚度矩阵统称为系统的特征矩阵，α为待求位移幅值向量，ω为待求自然频率。上述的m和k矩阵是系统的特征矩阵，其矩阵元素和单元基函数以及单元离散有关，这里我们采用二阶叠层基函数，基函数形式为：
其中
d.采用arpack软件包中的隐式重启arnoldi迭代方法求解c步骤所获得的本征方程，获得特征值λ和对应的特征向量α即振幅向量，arnoldi迭代方法每一步迭代所产生的线性方程采用基于三层预处理子的预处理共轭梯度法求解。
对式(8)的本征问题，我们采用隐式重启arnoldi迭代方法求解以下变形方程
其中ω0为预估最低模态频率。在隐式重启arnoldi迭代的过程中需要反复求解以下方程
px＝myj(11)
其中矩阵
nz＝r(12)
其中n为预处理子，z待求预处理向量，r为残差向量。
方程(12)式的求解依赖于预处理子n的性质，这里我们提出一种三层预处理子的方法来构造：首先把矩阵p按照x,y,z三个方向分块排序，获得如下分块矩阵形式
其中pxx为x方向的矩阵，pyy为y方向的矩阵，pzz为z方向的矩阵，pxy和pyx为x方向和y方向的耦合矩阵，pxz和pzx为x方向和z方向的耦合矩阵，pyz和pzy为y方向和z方向的耦合矩阵。我们取第一层预处理子的形式为
这是第一层预处理子，对于(14)式中的块对角矩阵pxx，pyy，pzz我们采用基于叠层基函数的性质来构造第二层预处理子，矩阵pxx，pyy，pzz所对应的第二层预处理子构造方法是一样的，我们分别给他们所对应的第二层预处理子取名为nx，ny，nz，这里以矩阵pxx为例来说明。首先按照一阶基函数、二阶基函数的排列顺序把pxx写成如下分块形式
其中p11为一阶基函数形成的矩阵，p22为二阶基函数形成的矩阵，p12和p21为一阶基函数和二阶基函数之间的耦合矩阵。我们取第二层预处理子为如下形式
对于(16)式中的块对角矩阵p11和p22我们采用多波前不完全cholesky分解来构造第三层预处理子，矩阵p11和p22所对应的第三层预处理子分别为nl和nh，他们的构造方法是一样的，同样我们以p11为例来说明。我们矩阵p11对进行多波前不完全cholesky分解，然后可以获得如下第三层预处理子的形式
其中d0为对角矩阵，p0为重排序矩阵，l为矩阵p11的cholesky分解下三角矩阵。这样三层预处理子构造完毕。把预处理子n应用于(12)式的求解，这样我们就能快速的求解(11)式，从而加快隐式重启arnoldi迭代的收敛速度，实现快速模态分析。通过arpack软件包中的隐式重启arnoldi迭代方法我们最终求出一系列特征值λi(i＝1,2,…,n)和对应的特征向量αi＝(i＝1,2,…,n)即振幅向量，其中n为所关心的模态数。
e.对d步骤获得的特征值和对应特征向量进行后处理获得振动模态频率和对应振动振型。
对d步骤获得特征值λi进行处理，对应的振动模态频率为
根据对应的特征向量αi，结合基函数性质进行有限元后处理可以得到求解域内的位移分布，这就是对应振动模态频率的振动振型。
图2是为实施例的有限元模型图；图3示出了本发明方案的实施例和ansys软件计算性能的对比，结果可以看出具体实施例的计算速度比ansys软件要快2倍，而内存消耗仅仅是ansys软件的76％。