CN108875195A - 一种考虑接触的三维力学随机振动仿真模拟方法 - Google Patents

Abstract

本发明属于三维力学振动分析数值求解技术领域，涉及一种考虑接触的三维力学随机振动仿真模拟方法。本发明采用一种线性接触的分析方法把接触问题引入线性模态分析，获得包含接触的振型，因此在采用模态叠加法进行随机振动分析时接触问题也随之包含进去，最终做相应的接触处理则获得考虑接触的随机响应。通过上述方法能够高效的获得高精度的数值计算结果，因此考虑接触的模态叠加法是提高随机振动数值模拟方法精度的有效途径。

Description

一种考虑接触的三维力学随机振动仿真模拟方法

技术领域

本发明属于三维力学结构振动分析数值求解技术领域，具体涉及一种考虑接触的三维力学随机振动仿真模拟方法。

背景技术

随机振动是一种统计意义下描述的振动，其特点在于：在任何给定的时刻，振动的幅值和相位均是不可知的,即无法用确定的数学公式或图形表示的那种振动，随机振动通过其振动幅值的统计特性进行描述。在工程实践中经常遇到这种振动，如：大型桥梁结构、电视塔或高层建筑在风作用下的振动，船舶或海上石油平台在海浪作用下的振动，车辆在不平路面上行驶产生的颠簸，地震响应等，这些都属于随机振动，人们一般采用理论计算法、实验测量法和数值计算法三种方法对上述现象进行分析。理论计算法只能对简单的结构进行分析，实验测量法成本太高，随着计算机技术的发展，随机振动的数值模拟方法显得越来越重要，尤其是复杂结构，如水坝、地基、核电站和水库等系统的三维随机振动分析。随着科学技术的进步，对结构的力学振动特性提出了更高的要求，因此需要高精度的随机振动数值模拟方法。然而大部分的结构振动是一种多体运动，部件之间存在接触关系，因此研究考虑接触的随机振动数值模拟方法是提高计算精度的一种途径。

目前，对复杂结构进行随机振动分析，有直接时域求解结构运动方程的方法，也有频域的功率谱密度方法。时域求解方法由于随机信号在时域的不确定性，一般是很少使用，因此频域求解方法是目前应用最广泛的。在频域求解方法中，一般包括频率响应法和模态叠加法。频率响应法，需要对随机激励信号进行傅里叶展开，这无疑增加了计算量，同时一般的频率响应法，是针对系统函数的方法，无法考虑接触问题。模态叠加法是基于线性模态分析获得的振动振型来进行随机振动分析，是一种线性分析过程；接触分析是一个非线性分析过程，因此现有的线性模态分析方法无法考虑接触，从而在随机振动分析时也无法考虑接触。这一缺点导致了现有随机振动分析方法的精度不高。

发明内容

针对上述存在问题或不足，为提高现有随机振动分析方法的计算精度；本发明提供了一种考虑接触的三维力学随机振动仿真模拟方法。

本发明技术方案，包括以下步骤：

A.将目标结构进行建模，引入位移边界条件和随机激励建立对应的有限元几何结构模型；

B.根据A步骤建立的几何结构模型各部件之间的相互关系，生成各部件之间的接触关系；

C.对A步骤所建几何结构模型采用四面体网格进行剖分，将连续的几何结构空间转化为离散空间；

D.根据B步骤生成的接触关系和C步骤生成的网格生成接触面上的联合面网格；

E.利用有限元法，将平衡微分方程、几何方程和物理方程等效的结构力学边值问题在步骤C建立的网格空间进行离散，结合步骤D建立的联合面网格，建立考虑接触问题的结构的有限元方程；

F.对步骤E得到的有限元方程引入目标结构的惯性力和阻尼力，得到有限元运动方程，并简化有限元运动方程获得结构的自由振动有限元广义本征方程；

G.求F步骤所获得的有限元广义本征方程，获得前q阶特征值λ_j(j＝1,2,…q)和对应的特征向量即振幅向量；

H.对G步骤获得的振幅向量进行质量矩阵归一化，获得q阶模态矩阵；

I.利用H步骤生成的模态矩阵采用模态叠加法对F步骤中获得的有限元运动方程进行解耦，获得q个互相独立的单自由度方程；

J.求解I步骤获得的q个互相独立的单自由度方程，获得相关函数矩阵，并根据维纳-辛钦关系获得随机位移响应的自功率谱密度矩阵S_αα(ω)；

K.对J步骤获得的自功率谱密度矩阵进行后处理获得目标结构的随机位移响应幅值的统计信息。

对于J步骤中的自功率谱密度矩阵S_αα(ω)其表现的是频域的信息，取其对角元素构成自功率谱密度自相关函数向量e_αα(ω)，在零均值的情况下获得随机位移响应插值系数的方差向量为

上式结果表示为，响应值小于一个标准偏差σ的概率为68.2％，响应值小于两个标准偏差2σ的概率为95.4％，响应值小于三个标准偏差3σ的概率为99.7％。上述计算值只是节点位移插值系数的方差，根据下式得到第m个单元的随机位移响应标准偏差分布函数

其中n_m为单元基函数的个数，为第m个单元中第i个节点上插值系数的标准偏差，N_i为体插值基函数。

对每一个单元执行上式，在包含接触的单元中，对应接触点上的插值系数为联合面网格中属于本单元的面网格上的插值系数。最终获得整体的位移标准偏差分布函数d(x,y,z)，代入结构内任意点的坐标值，即可获得该点的随机位移响应的幅值统计信息。

本发明采用一种线性接触的分析方法把接触问题引入线性模态分析，获得包含接触的振型，因此在采用模态叠加法进行随机振动分析时接触问题也随之包含进去，最终做相应的接触处理则获得考虑接触的随机响应。通过上述方法能够高效的获得高精度的数值计算结果，因此考虑接触的模态叠加法是提高随机振动数值模拟方法精度的有效途径。

附图说明

图1是本发明一种考虑接触的三维力学随机振动仿真模拟方法的流程图；

图2接触面联合网格示意图；

图3是实施例的有限元模型图；

图4是实施例与非接触三维力学随机振动仿真模拟方法计算值的对比图。

具体实施方式

下面结合附图和实施例来详细说明本发明的技术方案。

参照附图1，一种考虑接触的三维力学随机振动仿真模拟方法，包括以下步骤：

A.将目标结构进行建模，引入位移边界条件和随机激励建立对应的有限元几何结构模型；

建立目标结构的几何模型，根据结构的特性，引入位移边界条件和添加随机激励来仿真整个结构的随机振动特性，其结构如图1所示。具体的结构建模是结构力学数值计算中的一种公知过程，因此本步骤不再详细描述。

B.根据A步骤建立的几何结构模型各部件之间的相互关系，生成各部件之间的接触关系；

根据结构各部件的间的连接关系，生成相应的接触对。一般的前处理软件都能够实现该功能，这是一种公知的过程，因此本步骤不再详细阐述。

C.对A步骤所建几何结构模型采用四面体网格进行剖分，将连续的几何结构空间转化为离散空间；

采用四面体网格剖分仿真区域，剖分后的仿真区域被人为分割为多个三维四面体网格，从而将连续的几何结构空间转化为离散的网格空间。此时，接触面由一系列离散的面网格组成。由于四面体网格剖分是有限元方法中的一种公知过程，因此本步骤不再详细描述。

D.根据B步骤生成的接触关系和C步骤生成的网格生成接触面上的联合面网格；

在上述步骤C中生成的四面体网格在接触面上会存在一个公共面，由于接触面两边属于不同的部件，各部件具有不同的属性，因此接触面在不同的部件网格中具有不同的属性，需要建立一套联合的接触面网格。下面给出联合面网格的生成方法与步骤。

1)确定接触面网格的位置和所属四面体单元

根据步骤C生成的四面体网格信息和步骤B生成的接触对关系，确定每一个接触对上面网格以及其所属四面体单元。

2)复制接触单元面网格

对上述D步骤中的1)步骤所确定的接触单元面网格进行共形复制，如附图2所示面网格A'B'C'是接触单元面网格ABC的一个复制，附图2中面网格A'B'C'和面网格ABC在空间几何位置上是重叠的，之间是没有缝隙的，图中只是为了示意方便，描述成那样。这样在接触面上就形成了一套共形网格，其中面网格ABC属于四面体ABCD，面网格A'B'C'属于四面体A'B'C'D'，它们可以分别携带不同的材料属性。

E.利用有限元法，将平衡微分方程、几何方程和物理方程等效的结构力学边值问题在步骤C建立的网格空间进行离散，结合步骤D建立的联合面网格，建立考虑接触问题的结构的有限元方程；

对于空间边值问题，在结构(弹性体)内部我们要考虑静力学、几何学、物理学三方面条件，分别建立三套方程；并给定约束和激励，建立位移边界条件或应力边界条件。具体如下：

平衡微分方程

几何方程

物理方程

位移边界条件

应力边界条件

上述(1)(2)(3)(4)(5)式中，σ_x,σ_y,σ_z,τ_xy＝τ_yx,τ_yz＝τ_zy,τ_zx＝τ_xz表示求解区域中6个应力分量，ε_x,ε_y,ε_z,γ_xy,γ_yz,γ_zx表示求解区域中6个形变应力分量，u,v,w表示求解区域中3个位移分量。E是求解区域中结构的杨氏弹性模量，μ是求解区域中结构的泊松比，S_u表示位移边界面，S_σ表示应力边界面。l＝cos(n′,x),m＝cos(n′,y),n＝cos(n′,z)，表示应力边界面S_σ上的方向余弦，其中n′为应力边界面S_σ的外法线，x,y,z为应力边界面S_σ上三个方向的坐标值。为位移边界面S_u上的位移值，f_x,f_y,f_z为求解区域内结构受到的各个方向的体力，为应力边界面S_σ上受到的各个方向的面力，具体推导过程为一种公知过程，这里不再阐述。

经过有限元推导我们得到如下每一个四面体单元中的有限元方程，有限元法是一种众所周知的近似求解数理边值问题的数值技术，在结构力学中的应用已经很成熟，这里不再具体赘述。

K^mα^m＝F^m (6)

其中m为第几个单元，K^m为第m个单元的刚度矩阵，α^m为第m个单元的位移向量，F^m为第m个单元的外载荷量。具体表达式如下

K^m＝∫∫∫_ΩN^TL^TDLNdV (7)

式(7)、(8)、(9)中T是矩阵转置符号，为第m个单元的体力向量，为第m个单元的面力向量；α^m中位移的下标表示第几个插值点，n为体插值基函数(插值点)的个数，u,v,w表示三个位移分量。

N为体插值基函数的矩阵形式

N＝[N₁ N₂ … N_i … N_n] (12)

(13)式中N_i为体插值基函数，下标表示第几个插值点。对所有的四面体单元进行编号，同时对四面体单元内的插值点进行编号，最后去除重复的插值点，得到一组全局编号，该编号的个数即为整体系统的自由度，然后通过有限元系统装配我们可以得到如下整体结构有限元方程，具体装配过程是一种公知的过程，这里不再描述。

Kα＝F (14)

其中K为弹性体的刚度矩阵，α为结构位移向量，F为外载荷量。

其中n_f为系统总自由度。

通常我们将联合面网格h上的两个接触的点P和Q构成接触面上接触点对如附图2所示，他们的位移分别是和其位移可由所在接触网格面上的节点位移插值得到，则有

其中L_i是面插值基函数，为接触点所在单元的节点的位移矢量，n_S为面基函数(插值点)的个数，下标表示第几个插值点。这样一来，对于接触点P和Q间的相对位移可以表示为

其中

上面式(18)至(21)是在总体坐标系中定义，为方便引入接触条件，需要将其转换到局部坐标系中，即

其中T是两种坐标系之间的转换矩阵,上标T是矩阵转置符号。

为局部坐标系的三个单位基矢量，和式局部坐标系下P点和Q点的位移。

在黏结接触状态下局部坐标中一个接触点对的接触力引起的等效节点力向量为

其中和为第h个联合面网格单元的罚系数。

进一步可以得到第h个联合面网格单元整体坐标系下接触力等效节点力向量

或者写为

其中为第h个联合面网格单元的接触刚度矩阵。

对所有联合面网格单元的接触节点计算，并按照联合面网格所属体单元单元的编号，以及面插值点在体单元的位置组装到式(14)，则可得到系统的运动方程，即

(K+K_c)α＝F (27)

F.对步骤E得到的有限元方程引入目标结构的惯性力和阻尼力，得到有限元运动方程，并简化有限元运动方程获得结构的自由振动有限元广义本征方程；

当研究结构振动问题时，上述E步骤的α位移向量为时间的函数，我们重新定义

则根据E步骤讨论得到的有限元方程(27)，引入物体的惯性力和阻尼力得到有限元运动方程如下

其中

M＝∫∫∫_ΩρN^TNdΩ (30)

C＝∫∫∫_ΩνN^TNdΩ (31)

M为质量矩阵，ρ为求解区域物体的密度，ν为阻尼比，为α(t)对时间的二阶导数，为α(t)对时间的一阶导数。为时间相关的载荷向量，E为外力方向指示性向量，为运动加速度。

当物体自由振动时，此时F(t)＝0且无阻尼，方程(29)退化为

其振动形式叫做自由振动，该方程有解的形式

这是简谐振动的形式，其中ω为常数，为振幅向量，将其代入式(32)中，有

消去e^jωt后，有

该方程有非零解的条件是

|(K+K_c)-ω²M|＝0 (36)

这就是包含接触问的有限元广义本征方程。

G.求F步骤所获得的有限元广义本征方程，获得前q阶特征值λ_j(j＝1,2,…q)和对应的特征向量即振幅向量；

求解F步骤得到的有限元广义本征方程(36)，得到一系列的特征值λ_j(j＝1,2,…q)和对应的特征向量具体的求解方法有很多且是一种公知的过程，这里不再详细阐述。由于这里得到的特征向量包含了接触单元节点上的值，因此特征向量有如下表达形式

H.对G步骤获得的振幅向量进行质量矩阵归一化，获得q阶模态矩阵；

对于所有前q阶特征向量进行比例缩放获得质量归一化特征向量使其满足下述质量归一化方程

取前q阶质量归一化特征向量组成如下形式的模态矩阵

I.利用H步骤生成的模态矩阵采用模态叠加法对F步骤中获得的有限元运动方程进行解耦，获得q个互相独立的单自由度方程；

令

ξ(t)＝[ξ(t)₁ ξ(t)₂ … ξ(t)_j … ξ(t)_q]^T (41)

其中ξ_j(t)为第j阶模态的模态坐标。将Φ^T左乘式(29)各项，并将式(41)代入式(29)得

其中

为ξ(t)对时间的二阶导数，为ξ(t)对时间的一阶导数，ν_j为第j阶模态对应的阻尼比。

方程(42)可以被分解成q个相互独立的单自由度方程

其中γ_j为第阶模态的参与系数

J.求解I步骤获得的q个互相独立的单自由度方程，获得相关函数矩阵，并根据维纳-辛钦关系获得随机位移响应的自功率谱密度矩阵S_αα(ω)；

式(46)的解在时域内表示为

其中h_j(t)为第j阶模态的脉冲响应函数。将式(48)代入式(40)得到

于是相关函数矩阵为

其中为加速度的相关函数。应用维纳-辛钦关系，对上式进行傅里叶变换得到随机位移响应的自功率谱密度矩阵

其中上标“*”代表取复共轭，H_j(ω)为第j阶模态的频响函数，为加速度的功率谱密度函数。

K.对J步骤获得的自功率谱密度矩阵进行后处理获得结构的随机位移响应幅值的统计信息。

对于J步骤中的自功率谱密度矩阵S_αα(ω)其表现的是频域的信息，取其对角元素构成自功率谱密度自相关函数向量e_αα(ω)，在零均值的情况下可以获得随机位移响应插值系数的方差向量为

上式结果表示为，响应值小于一个标准偏差σ的概率为68.2％，响应值小于两个标准偏差2σ的概率为95.4％，响应值小于三个标准偏差3σ的概率为99.7％。上述计算值只是节点位移插值系数的方差，根据下式得到第m个单元的随机位移响应方差分布函数

其中n_m为单元基函数的个数，为第m个单元中第i个节点上插值系数的方差。对每一个单元执行(53)式，在包含接触的单元中，对应接触点上的插值系数为联合面网格中属于本单元的面网格上的插值系数。最终获得整体的位移方差分布函数d(x,y,z)，代入结构内任意点的坐标值，即可获得该点的随机位移响应的幅值统计信息。

图3是为实施例的有限元模型图；图4示出了实施例与非接触三维力学随机振动仿真模拟方法计算值的对比图。从图4的结果中可以看出具体实施例的计算值比非接触三维力学随机振动仿真模拟方法的计算值要大，这更加接近真实值。

Claims (1)

1.一种考虑接触的三维力学随机振动仿真模拟方法，包括以下步骤：

A.将目标结构进行建模，引入位移边界条件和随机激励建立对应的有限元几何结构模型；

B.根据A步骤建立的几何结构模型各部件之间的相互关系，生成各部件之间的接触关系；

C.对A步骤所建几何结构模型采用四面体网格进行剖分，将连续的几何结构空间转化为离散空间；

D.根据B步骤生成的接触关系和C步骤生成的网格生成接触面上的联合面网格；

E.利用有限元法，将平衡微分方程、几何方程和物理方程等效的结构力学边值问题在步骤C建立的网格空间进行离散，结合步骤D建立的联合面网格，建立考虑接触问题的结构的有限元方程；

F.对步骤E得到的有限元方程引入目标结构的惯性力和阻尼力，得到有限元运动方程，并简化有限元运动方程获得结构的自由振动有限元广义本征方程；

G.求F步骤所获得的有限元广义本征方程，获得前q阶特征值λ_j(j＝1,2,…q)和对应的特征向量即振幅向量；

H.对G步骤获得的振幅向量进行质量矩阵归一化，获得q阶模态矩阵；

I.利用H步骤生成的模态矩阵采用模态叠加法对F步骤中获得的有限元运动方程进行解耦，获得q个互相独立的单自由度方程；

J.求解I步骤获得的q个互相独立的单自由度方程，获得相关函数矩阵，并根据维纳-辛钦关系获得随机位移响应的自功率谱密度矩阵S_αα(ω)；

K.对J步骤获得的自功率谱密度矩阵进行后处理获得目标结构的随机位移响应幅值的统计信息；

所述后处理具体如下：

对于J步骤中的自功率谱密度矩阵S_αα(ω)其表现的是频域的信息，取其对角元素构成自功率谱密度自相关函数向量e_αα(ω)，在零均值的情况下获得随机位移响应插值系数的方差向量为

上式结果表示为，响应值小于一个标准偏差σ的概率为68.2％，响应值小于两个标准偏差2σ的概率为95.4％，响应值小于三个标准偏差3σ的概率为99.7％；上述计算值只是节点位移插值系数的方差，根据下式得到第m个单元的随机位移响应标准偏差分布函数

其中n_m为单元基函数的个数，为第m个单元中第i个节点上插值系数的标准偏差，N_i为体插值基函数。

对每一个单元执行上式，在包含接触的单元中，对应接触点上的插值系数为联合面网格中属于本单元的面网格上的插值系数；最终获得整体的位移标准偏差分布函数d(x,y,z)，代入结构内任意点的坐标值，即可获得该点的随机位移响应的幅值统计信息。