CN101887474B - 基于有限元法与广义傅里叶级数法的结构振动分析方法 - Google Patents

Abstract

本发明的目的在于提供基于有限元法与广义傅里叶级数法的结构振动分析方法。分为以下步骤：将需要进行振动分析的结构区域划分，分别形成相应的有限元表述区域与广义傅里叶级数表述区域；对有限元表述区域进行有限元网格划分并形成相应的质量刚度矩阵，依据广义傅里叶级数表述区域特点选择相应的假设位移形式，形成质量刚度矩阵；之后在两个区域之间建立虚拟弹簧利用能量变分方法，将虚拟弹簧势能转化为总体耦合刚度矩阵；然后对形成的质量刚度矩阵根据位移进行排列，形成总体结构质量刚度矩阵；求解线性方程组得到相应的节点位移和级数展开中的未知系数。本发明对大型较复杂结构，既可以获得比有限元法更高的精度，又可以节省大量的计算成本。

Description

基于有限元法与广义傅里叶级数法的结构振动分析方法

技术领域

本发明涉及的是一种应用于工程力学和振动工程领域的结构振动分析方法。

背景技术

作为一种数值计算方法，有限元法在结构振动领域有着广泛的应用。在处理不规则复杂结构，该方法有着解析方法无法比拟的优势。但有限元法仍然存在一些缺点，如在计算大型复杂结构特别是高频问题甚至中频问题，由于现在计算技术限制，计算所需网格过多，导致求解非常的困难甚至无法求解。

在结构振动分析中，解析法具有结果精确，使用计算机求解具有计算结果的频带宽、精确可信以及占有计算机资源小，计算速度快等优点。但是，解析法往往只适用于简单规则结构。在解析方法中，傅里叶级数解法是近年来受到重视的方法之一，它可以适用于各种边界条件、运算方便，并由于级数间具有正交性，可以使计算量大为减小并保证很高精度。美国韦恩州立大学李文龙提出一种广义傅里叶级数方法成功解决了任意边界条件下横梁的弹性振动(W.L.Li，Vibrationanalysis of rectangular plates with general elastic boundary supports，Journal of Sound andVibration 273(2004)619-635.)，求解出的固有频率和振型的精度及级数收敛速度都达到了十分理想的效果，此解法还被应用于求解和分析多跨度桥梁接受运动负载时的振动问题(W.L.Li，M.Daniels，A Fourier series method for the vibrations of elasticallyrestrained plates arbitrarily loaded with springs and masses，Journal of Sound and Vibration 252(2002)768-781.)。文献Vibrations of rectangular plates with arbitrary non-uniform elastic edgerestraints(X.Zhang，Wen L.Li*，Journal of Sound and Vibration 326(2009)221-234)中，国内杜敬涛等人将傅里叶级数方法解决了任意边界弹性边界矩形板振动问题，不均匀边界问题乃至是板与板耦合振动分析问题。对于梁、板、圆柱壳等规则结构，这些结构的微分方程是四阶的，展开级数具有四阶(或更高阶)的逐项可导的性质，这些结构均可用广义傅里叶级数方法进行结构振动求解。对于不规则结构或者规则结构的布尔运算，如矩形板上有三角形孔，广义傅里叶级数方法还无法求解，原因是一般结构的高阶逐项可导的条件很难满足，所以，傅里叶级数解法解决问题的范围受到很大限制。

将有限元法解决复杂结构能力强的优点和傅里叶级数展开法计算精度高、计算资源消耗小且计算速度快的优点有机结合，可以解决大型较复杂结构振动问题。目前还没有将两种方法结合的技术出现以及相关文献报道。

发明内容

本发明的目的在于提供可用来解决大型较复杂结构振动特别是中高频振动计算难题的基于有限元法与广义傅里叶级数法的结构振动分析方法。

本发明的目的是这样实现的：

本发明基于有限元法与广义傅里叶级数法的结构振动分析方法，其特征是：

(1)将待分析的结构划分为规则结构与非规则结构两部分，两部分之间利用虚拟弹簧连接；

(2)规则结构部分设定位移函数，即设定傅里叶级数及其附加容许函数，依据能量原理将傅里叶级数及其容许函数中的系数作为未知量，规则结构部分转化为当量刚度矩阵与质量矩阵，非规则结构部分采用有限元法构造该区域的总体刚度矩阵与质量矩阵；

(3)规则结构部分与非规则结构部分之间的虚拟弹簧两端位移分别以有限元节点位移以及级数展开的系数表示，利用弹簧两端位移表示弹簧储存势能，利用能量方法中的变分方法得到有限元离散区域与级数展开区域之间的耦合质量矩阵与刚度矩阵；

(4)有限元质量刚度矩阵、广义傅里叶级数展开质量刚度矩阵和虚拟弹簧的耦合质量刚度矩阵依据位移排列组合得到结构总体质量刚度矩阵；

(5)由结构总体质量刚度矩阵得到线性方程组，求解线性方程组得到相应的节点位移和级数展开中的未知系数。

本发明的优势在于：与有限元法相比，傅里叶级数展开方法作为解析方法虽然只能用于一些规则结构但该方法无需任何网格且有着收敛快速，计算所需资源少的优点；而有限元方法作为一种成熟技术，在结构领域已经得到非常广泛的应用，适用于任意形状的结构。由于结构的大型化复杂化，加上计算机的计算能力限制，导致有限元法计算有一定误差以及计算效率低下，对于一些复杂结构的中高频甚至无法求解。通过本发明提出的将两者结合的连接方法，对大型较复杂结构，既可以获得比有限元法更高的精度，又可以节省大量的计算成本。

附图说明

图1为本发明的流程示意图；

图2是本发明实施方式1的结构示意图。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

结合图1，本发明可以分为以下步骤：

1、建立求解域，并将之分为规则结构如梁、矩形板等与非规则结构两部分，二者之间利用虚拟无穷大弹簧连接；

2、非规则部分采用有限元法。将非规则结构离散化成有限元单元，即将该区域分解成节点和单元；

3、假设代表单元物理行为的形函数，即假设代表单元解的近似连续函数，并对单元建立方程；

4、将单元组合成总体的问题，构造该区域的总体刚度矩阵与质量矩阵；

5、规则区域根据其控制方程特点，设定位移函数，即设定傅里叶级数及其附加容许函数；

6、依据能量原理将傅里叶级数及其容许函数中的系数作为未知量，规则区域转化为当量刚度矩阵与质量矩阵；

7、应用边界条件到不同区域的质量与刚度矩阵之中；

8、规则区域与非规则区域之间的虚拟弹簧两端位移分别以有限元节点位移以及级数展开的系数表示之；

9、利用弹簧两端位移，表示弹簧储存势能，同样的利用能量方法中的变分方法对之进行处理，得到有限元离散区域与级数展开区域之间的耦合质量矩阵与刚度矩阵。

10、分别利用形成的有限元质量刚度矩阵、傅里叶级数展开质量刚度矩阵以及连接二者的虚拟弹簧的耦合矩阵形成总体结构的总体质量与刚度矩阵。

11、求解得到的线性方程组，得到相应的节点位移和级数展开中的未知系数。

12、通过计算结果得到其他信息，如频率、各阶模态，亦可进一步做结构响应计算。

实施方式1：

结合图2，悬臂梁支撑左侧部分由有限元fa表示，支撑右侧部分由傅里叶级数展开方法表述。梁弯曲问题的基本方程可表示如下：

几何关系 $\mathrm{\κ}=-\frac{d^{2}w}{{\mathrm{dx}}^2}---\left(1\right)$

应力应变关系 $M=\mathrm{EI\κ}=-\mathrm{EI}\frac{d^{2}w}{d{x}^2}---\left(2\right)$

平衡方程 $Q=\frac{\mathrm{dM}}{\mathrm{dx}}=-\mathrm{EI}\frac{d^{2}w}{{\mathrm{dx}}^2}---\left(3\right)$

$-\frac{\mathrm{dQ}}{\mathrm{dx}}=\mathrm{EI}\frac{d^{4}w}{{\mathrm{dx}}^4}=q\left(x\right)---\left(4\right)$

$w=\stackrel{\‾}{w}$ $\frac{\mathrm{dw}}{\mathrm{dx}}=\stackrel{\‾}{\mathrm{\θ}}$

边界条件

或 $w=\stackrel{\‾}{w}$ $M=\stackrel{\‾}{M}---\left(5\right)$

或 $Q=\stackrel{\‾}{Q}$ $M=\stackrel{\‾}{M}$

以上各式中κ是梁中面变形后的曲率；M和Q分别是截面上的弯矩和横向剪力；I是截面上的惯性矩；

分别是边界上给定的挠度、转动、弯矩和剪力。

用有限元法分析梁弯曲问题时，采用Hermite多项式作为单元的插值函数。对于两个节点的一维单元，函数φ采用Hermite多项式的插值表达式可写成

$\mathrm{\φ}\left(\mathrm{\ξ}\right)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{H}_i^{\left(0\right)}\left(\mathrm{\ξ}\right){\mathrm{\φ}}_i+\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{H}_i^{\left(1\right)}\left(\mathrm{\ξ}\right){\left(\frac{\mathrm{d\φ}}{\mathrm{d\ξ}}\right)}_i---\left(6\right)$

因此单元内挠度函数的插值表示如下

$w\left(\mathrm{\ξ}\right)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{H}_i^{\left(0\right)}\left(\mathrm{\ξ}\right)w_i+\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{H}_i^{\left(1\right)}\left(\mathrm{\ξ}\right){\left(\frac{\mathrm{dw}}{\mathrm{d\ξ}}\right)}_i$ (7)

$=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{N}_i\left(\mathrm{\ξ}\right)a_i={\mathrm{Na}}^e$

对泛函

${\mathrm{\Π}}_p\left(w\right)={\∫}_0^l\frac{1}{2}\mathrm{EI}{\left(\frac{d^{2}w}{{\mathrm{dx}}^2}\right)}^2\mathrm{dx}-{\∫}_0^l\frac{1}{2}\mathrm{EI}{\left(\frac{d^{2}w}{{\mathrm{dx}}^2}\right)}^2\mathrm{dx}-\underset{j}{\mathrm{\Σ}}{P}_j{w}_j+\underset{k}{\mathrm{\Σ}}{M}_k{\left(\frac{\mathrm{dw}}{\mathrm{dx}}\right)}_k---\left(8\right)$

取变分可以得到有限元的刚度矩阵

$K_e=\frac{\mathrm{EI}}{l^3}\left[\begin{array}{cccc}12& 6l& -12& 6l\\ 6l& {4l}^2& -6l& {2l}^2\\ -12& -6l& 12& -6l\\ 6l& 2l^2& -6l& {4l}^2\end{array}\right]---\left(9\right)$

同理得到单元的质量矩阵

$M_e=\frac{\mathrm{\ρAl}}{420}\left[\begin{array}{cccc}156& 22l& 54& -13l\\ 22l& 4l^2& 13l& -{3l}^2\\ 54& 13l& 156& -22l\\ -13l& -3l^2& -22l& {4l}^2\end{array}\right]---\left(10\right)$

由此，有限元形成结构自由振动方程如式(11)所示：

$(K_e-M_e{\mathrm{\ω}}^2)\stackrel{\→}{u_e}=0---\left(11\right)$

其中

代表有限单元中的每个单元节点上的值；ω为振动的圆频率。

解决弹性支承问题通常采用修改刚度矩阵的方法。以上述连续梁为例，可在相应节点上主刚度加相应值。

依据广义傅里叶级数方法及梁的特性假设位移条件形式如下

$w\left(x\right)=\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}{A}_m\cos {\mathrm{\λ}}_{\mathrm{am}}x+\underset{l=1}{\overset{4}{\mathrm{\Σ}}}\left(C_l{\mathrm{\ξ}}_a^l\left(x\right)\right)---\left(\text{12}\right)$

其中

为容许函数，吸收傅里叶级数方法边界不连续，加快傅里叶级数的收敛性。

${\mathrm{\ξ}}_a^1\left(x\right)=\frac{9a}{4\mathrm{\π}}\sin \left(\frac{\mathrm{\πx}}{2a}\right)-\frac{a}{12\mathrm{\π}}\sin \left(\frac{3\mathrm{\πx}}{2a}\right)---(13.a)$

${\mathrm{\ξ}}_a^2\left(x\right)=-\frac{9a}{4\mathrm{\π}}\cos \left(\frac{\mathrm{\πx}}{2a}\right)-\frac{a}{12\mathrm{\π}}\cos \left(\frac{3\mathrm{\πx}}{2a}\right)---(13.b)$

${\mathrm{\ξ}}_a^3\left(x\right)=\frac{a^3}{{\mathrm{\π}}^3}\sin \left(\frac{\mathrm{\πx}}{2a}\right)-\frac{a^3}{3{\mathrm{\π}}^3}\sin \left(\frac{3\mathrm{\πx}}{2a}\right)---(13.c)$

${\mathrm{\ξ}}_a^4\left(x\right)=-\frac{a^3}{{\mathrm{\π}}^3}\cos \left(\frac{\mathrm{\πx}}{2a}\right)-\frac{a^3}{3{\mathrm{\π}}^3}\cos \left(\frac{3\mathrm{\πx}}{2a}\right)---(13.d)$

因此梁的曲率可表达如下：

$\frac{d^{2}w}{{\mathrm{dx}}^2}=-\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}{A}_m{{\mathrm{\λ}}_{\mathrm{am}}}^2\cos {\mathrm{\λ}}_{\mathrm{am}}x+\underset{l=1}{\overset{4}{\mathrm{\Σ}}}{C}_l{\left({\mathrm{\ξ}}_a^l\left(x\right)\right)}^{\′\′}---\left(14\right)$

梁的势能动能表达式如下

$U=\frac{1}{2}\underset{0}{\overset{L}{\∫}}\mathrm{EI}{\left(\frac{d^{2}w}{{\mathrm{ds}}^2}\right)}^2\mathrm{dx}$ (15)

$T=\frac{1}{2}{\mathrm{\ω}}^2\mathrm{\ρA}\underset{0}{\overset{L}{\∫}}\left(u^2\right)\mathrm{dx}$

一个边界上有两个弹簧：分别为线性弹簧k_w0，k_w1，扭簧K_w0，K_w1。这样的话，边界上的能量表达式可以写为如下形式：

$U_{b0}=\frac{1}{2}{k}_{w0}{w}^2+\frac{1}{2}{K}_{w0}{\left(\frac{\mathrm{dw}}{\mathrm{dx}}\right)}^2,(x=0)$ (16)

$U_{b1}=\frac{1}{2}{k}_{w1}{w}^2+\frac{1}{2}{K}_{w1}{\left(\frac{\mathrm{dw}}{\mathrm{dx}}\right)}^2,(x=a)$

因此，系统中储存的总势能为V＝U+U_b0+U_b1，直梁拉格朗日函数可以表示为如下形式

L＝V-T    (17)

将式(12)～(16)代入式17，采用Rayleigh-Ritz法使拉格朗日函数对每个未知Fourier系数取极值，我们可以得到8个线性方程组，进一步写为矩阵表达式形式：

$(K_s-M_s{\mathrm{\ω}}^2)\stackrel{\→}{A_s}=0---\left(18\right)$

其中 $\stackrel{\→}{A_s}={[A_0,A_1,...,A_M,C_1,C_2,C_3,C_4]}^{\′}$

通过以上推导，分别得到有限元与广义傅里叶级数展开方法的质量刚度矩阵。推导过程可以发现其质量刚度矩阵都由有变分导出，因此对于左右梁之间的弹簧k_c，K_c作一个相同的处理，从而得到其耦合刚度与质量矩阵。

耦合项能量，

$U_{c1}=\frac{1}{2}{k}_c{(w_n-(\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}{A}_m+\underset{l=1}{\overset{4}{\mathrm{\Σ}}}\left(C_l{\mathrm{\ξ}}_a^l\left(0\right)\right)))}^2$ (19)

$U_{c2}=\frac{1}{2}{K}_c{({\mathrm{\θ}}_n-\underset{l=1}{\overset{4}{\mathrm{\Σ}}}\left({C_l\left({\mathrm{\ξ}}_a^l\left(0\right)\right)}^{\′}\right))}^2$

在耦合能量项中，存在左端有限元梁的最后一个节点的两个位移以及右端的Fourier系数，同样的对每个未知系数取极值，可以得到5个线性方程组。通过这五个线性方程组，将左侧的有限元梁与右侧的广义傅里叶级数梁耦合在一起使之成为一体。这五个线性方程组与式(11)、(18)合在一起整理之，可以得到这样的总体质量刚度矩阵

$(K-M{\mathrm{\ω}}^2)\stackrel{\→}{A}=0---\left(20\right)$

其中 $\stackrel{\→}{A}={[w_1,{\mathrm{\θ}}_1,w_2,{\mathrm{\θ}}_2...w_n,{\mathrm{\θ}}_n,A_0,A_1,...,A_M,C_1,C_2,C_3,C_4]}^{\′}$

通过式(20)即可求解相应结构的振动特性。

根据以上推导，对一根直梁几种不同边界(S-S，C-C，F-F)进行了验算，其参数如下所示：

长度为L＝1m；

横截面积＝0.0003；

惯性矩I＝7.0e-009；

杨氏模量E＝2.1e11。

计算所得结果如下所示，耦合方法结果与有限元，广义傅里叶级数方法结果相差很小。

两边简支边界

两边固支边界

两边自由边界

双边自由边界条件第五十阶达到精确值0.001％精度时的计算时间比较对比结果如下：

Claims (1)

1.基于有限元法与广义傅里叶级数法的结构振动分析方法，其特征是：

(1)将待分析的结构划分为规则结构与非规则结构两部分，两部分之间利用虚拟弹簧连接；

(2)规则结构部分设定位移函数，即设定傅里叶级数及其附加容许函数，依据能量原理将傅里叶级数及其容许函数中的系数作为未知量，规则结构部分转化为当量刚度矩阵与质量矩阵，非规则结构部分采用有限元法构造该区域的总体刚度矩阵与质量矩阵；

(3)规则结构部分与非规则结构部分之间的虚拟弹簧两端位移分别以有限元节点位移以及级数展开的系数表示，利用弹簧两端位移表示弹簧储存势能，利用能量方法中的变分方法得到有限元离散区域与级数展开区域之间的耦合质量矩阵与刚度矩阵；

(4)有限元质量刚度矩阵、广义傅里叶级数展开质量刚度矩阵和虚拟弹簧的耦合质量刚度矩阵依据位移排列组合得到结构总体质量刚度矩阵：

一个边界上有两个弹簧：分别为线性弹簧k_w0，k_w1、扭簧K_w0，K_w1，边界上的能量表达式为： $\begin{array}{cc}{U}_{b0}=\frac{1}{2}{k}_{w0}{w}^2+\frac{1}{2}{K}_{w0}{\left(\frac{\mathrm{dw}}{\mathrm{dx}}\right)}^2& (x=0)\\ U_{b1}=\frac{1}{2}{k}_{w1}{w}^2+\frac{1}{2}{K}_{w1}{\left(\frac{\mathrm{dw}}{\mathrm{dx}}\right)}^2& (x=a)\end{array},$ w为梁的横向位移，k_w0为分布在梁结构x＝0处的线性弹簧，k_w1为分布在梁结构x＝a处的线性弹簧，K_w0为分布在梁结构x＝0处的扭转弹簧，K_w1为分布在梁结构x＝a处的扭转弹簧，x为坐标，a为梁结构的长度，

系统中储存的总势能为V＝U+U_b0+U_b1，直梁拉格朗日函数表示为如下形式L＝V-T，其中T为系统的动能，V为系统势能，U为梁结构本身的势能，U_b0+U_b1为储存在边界弹簧上的势能，

采用Rayleigh-Ritz法使拉格朗日函数对每个未知Fourier系数取极值，得到8个线性方程组，进一步写为矩阵表达式形式：

其中

K_s为形成的系统刚度矩阵，M_s为形成的系统质量矩阵，ω为系统圆频率，

为未知的傅里叶级数相关系数，A₀，A₁，…，A_M为原傅里叶级数的系数，C₁，C₂，C₃，C₄为附加项的系数，

分别得到有限元与广义傅里叶级数展开方法的质量刚度矩阵，耦合项能量：

$\begin{array}{c}{U}_{c1}=\frac{1}{2}{k}_c{(w_n-(\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}{A}_m+\underset{l=1}{\overset{4}{\mathrm{\Σ}}}\left(C_l{\mathrm{\ξ}}_a^l\left(0\right)\right)))}^2\\ U_{c2}=\frac{1}{2}{K}_c{({\mathrm{\θ}}_n-\underset{l=1}{\overset{4}{\mathrm{\Σ}}}\left(C_l{\left({\mathrm{\ξ}}_a^l\left(0\right)\right)}^{\′}\right))}^2\end{array},$ k_c为假设的连接两种结构的线弹簧，K_c为假设的连接两种结构的扭转弹簧，U_c1、U_c2分别为储存在线弹簧k_c与扭转弹簧K_c上的势能，w_n、θ_n为有限元结构n节点上的横向位移与转角位移，A_m、C_l为傅里叶系数，

为用于提高傅里叶级数收敛性的附加函数，

在耦合能量项中，存在左端有限元梁的最后一个节点的两个位移以及右端的Fourier系数，同样的对每个未知系数取极值，得到5个线性方程组，通过这五个线性方程组，将左侧的有限元梁与右侧的广义傅里叶级数梁耦合在一起使之成为一体，最后得到总体质量刚度矩阵：

其中

K为形成的总体系统刚度矩阵，M为形成的总体系统质量矩阵，ω为系统圆频率，

为表述结构位移的未知参数，w₁，θ₁，w₂，θ₂…w_n，θ_n为有限元结构中的n个节点的横向位移与转角位移，A₀，A₁，…，A_M为原傅里叶级数的系数，C₁，C₂，C₃，C₄为附加项的系数；

(5)由结构总体质量刚度矩阵得到线性方程组，求解线性方程组得到相应的节点位移和级数展开中的未知系数。

Images