CN113656984A

CN113656984A - 一种基于蒙特卡洛法的梁结构固有频率和振型计算方法

Info

Publication number: CN113656984A
Application number: CN202111011060.XA
Authority: CN
Inventors: 祝磊; 张建勋; 孙海林
Original assignee: Beijing University of Civil Engineering and Architecture; China Architecture Design and Research Group Co Ltd
Current assignee: Beijing University of Civil Engineering and Architecture; China Architecture Design and Research Group Co Ltd
Priority date: 2021-08-31
Filing date: 2021-08-31
Publication date: 2021-11-16
Anticipated expiration: 2041-08-31
Also published as: CN113656984B

Abstract

本发明公开了一种基于蒙特卡洛法的梁结构固有频率和振型计算方法，该计算方法包括如下步骤：步骤一,设定初始参数；步骤二，选择合适的函数类型，随机生成初始模态函数；步骤三，求解体系的动能、应变能和固有频率；步骤四，校验所得结果是否满足收敛条件；步骤五，得到体系的固有频率和模态函数。该发明计算方法可根据使用者所期望的计算精度和时间成本自由控制计算误差限，以获得最合适的结果；利用MATLAB编程实现了蒙特卡洛法对梁结构的固有频率和振型曲线计算过程，可以快速准确得到计算结果，效果显著。

Description

一种基于蒙特卡洛法的梁结构固有频率和振型计算方法

技术领域

本发明涉及梁结构动力学分析领域，具体是一种基于蒙特卡洛法的梁结构固有频率和振型计算方法。

背景技术

梁是在土木和机械等工程领域应用最为广泛的构件之一，如建筑结构中的主梁次梁，单层厂房结构中的吊车梁，多跨桥梁结构，吊车和吊塔的悬臂，风力发电机中的叶片和钻探石油钻井的钻杆等，这些结构和构件在通常工况或极限工况下都不可避免地会产生振动，作为重要的受力构件，对其振动特性进行的研究具有重要的理论价值和应用价值。

能影响梁结构振动特性的因素较多，如横截面形状、截面面积、材料弹性模量、密度、梁长等因素。现有研究中已对一些经典约束条件和部分弹性约束条件特殊情况下的梁振动特性利用解析法获得了固有频率和振型函数，但此类方法对应的固有频率特征方程往往很复杂，且考虑到在理论上确定梁横向弯曲振动偏微分方程边界条件的困难，往往只能解决结构形状较简单的等截面梁的问题，对于变截面梁试图按经典的弹性力学方法获得解析解是十分困难的。

在理论研究中利用具体程序对其进行实现的方法相对较少，因此难以在工程实际计算中得到应用并指导技术问题的解决。对于较复杂的梁结构，一般可使用有限元法进行求解，如ANSYS和ABAQUS等。但此类大型通用有限元软件虽然性能强大但学习成本较高。基于此，发明一种适用于一般梁弯曲振动的基频和振型的计算方法和程序尤为必要。

发明内容

本发明的目的在于提供一种基于蒙特卡洛法的梁结构固有频率和振型计算方法，以解决背景技术中的问题。

为实现上述目的，本发明提供如下技术方案：

一种基于蒙特卡洛法的梁结构固有频率和振型计算方法，包括如下步骤：

步骤一,设定初始参数

(3)几何属性参数和材料属性参数；

(4)将梁沿长度方向划分成等长的n段，并在各节点处生成坐标元素，合成坐标向量x；

步骤二，选择合适的函数类型，随机生成初始模态函数

(3)随机生成初始模态函数系数向量；

(4)选择满足计算体系边界条件的坐标函数，与坐标函数合成模态函数；

步骤三，求解体系的动能、应变能和固有频率

(4)计算各节点的横向位移；

(5)利用位移计算体系的最大动能T_max和利用位移计算体系的最大应变能V_max；

(6)通过瑞利法求解体系固有频率ω；

步骤四，校验所得结果是否满足收敛条件

计算的收敛通过固有频率是否满足收敛要求来判断，各次循环过程中计算所得频率最小值为ω_min，最后一次循环中计算所得频率值为ω_j，若满足ω_min＝ω_j，则结果收敛，终止计算；

设定一个误差限ε得到近似平衡条件，通过限定平衡方程等式左侧计算结果中两个频率的相对误差限ε’来判断系统是否平衡，若ε’＜ε则进行步骤五，若ε’＞ε，则进行步骤二；

步骤五，得到体系的固有频率和模态函数

输出各节点坐标，根据坐标绘制梁结构第一阶振型，得到最终的振动模态，并得到计算过程中其他的过程变量的变化情况。

在上述技术方案的基础上，本发明还提供以下可选技术方案：

在一种可选方案中：所述步骤一中几何参数包括梁长L、宽b、高h、截面积A、截面惯性矩I和划分段数n；材料属性参数包括材料弹性模量E、密度ρ和梁每延米的质量m。

在一种可选方案中：所述步骤二中模态函数的维数等同于组合为模态函数的坐标函数个数。

在一种可选方案中：所述步骤五中在计算梁的第一个和最后一个节点应变能时，第1个节点的应变能V₁按照式V₁＝V₂计算,第n+1个节点的应变能V_n+1按照式V_n+1＝V_n计算。

在一种可选方案中：计算单元的刚度S时，第k个单元的刚度按照式S_k＝E_k·A_k/lo_k计算。

相较于现有技术，本发明的有益效果如下：

1.可根据使用者所期望的计算精度和时间成本自由控制计算误差限，以获得最合适的结果；

2.利用MATLAB编程实现了蒙特卡洛法对梁结构的固有频率和振型曲线计算过程，可以快速准确得到计算结果，效果显著。

附图说明

图1为本发明的一个实施例中的具体实施方式的流程图。

图2为本发明的实施例2中锥形截面悬臂梁示意图。

图3为本发明的锥形截面悬臂梁主视图。

图4为本发明的锥形截面悬臂梁俯视图。

图5为等截面简支梁工况下本发明所求振型函数与理论解的对比图。

图6为等截面悬臂梁工况下本发明所求振型函数的对比图。

图7为变截面悬臂梁工况下本发明所求振型函数的对比图。

具体实施方式

本发明所列举的各实施例仅用以说明本发明，并非用以限制本发明的范围。对本发明所作的任何显而易知的修饰或变更都不脱离本发明的精神与范围。

实施例1：

如图1所示，本发明提出一种基于蒙特卡洛法的梁结构固有频率和振型计算方法，包括如下步骤：

步骤一、设定初始参数：

设定梁的几何属性：在等截面梁中，记梁长为L，宽为b，高为h，截面积为：

A＝bh

截面惯性矩为：

I＝bh³/12；

考虑到变截面梁中变截面的形式较多，在此以梁高按线性变化，厚度不变为例：记梁长为L，宽为b；

梁高为：

截面积为：

截面惯性矩为：

设定梁的材料属性：记材料弹性模量为E，密度为ρ；

等截面梁每延米的质量为：

m＝ρA

变截面梁每延米的质量为：

考虑到梁结构本身属于分布参数体系的一种，为使固有频率收敛于精确解，从理论上需要取无限个坐标，因此在该方法中直接求解其精确解是不现实的，可以考虑求解其近似的数值解，故对连续体系进行离散化处理，将梁划分为等长的n段，则共有n+1个节点，均沿梁长度方向上分布，各节点坐标易得，将其储存在坐标向量x中。

步骤二、选择合适的函数类型，随机生成初始模态函数：

确定模态函数中所包含的坐标函数个数，随机生成介于-1到1之间的随机数，其个数等同于坐标函数个数。

根据待求解的构件边界条件确定所选择的模态函数，部分经典边界下的模态函数设定如下：

两端固定梁：

简支梁：

悬臂梁：

选择满足计算体系边界条件的坐标函数，与坐标函数合成模态函数；

步骤三、求解体系的动能、应变能和固有频率：

将坐标向量x带入所选取的模态函数中，计算各节点的横向位移；利用位移计算体系的最大动能，公式为：

利用位移计算体系的最大应变能，仅考虑弯曲变形，公式为：

通过瑞利法求解体系固有频率，可得：

在离散化的体系中，以上各项的计算公式分别为：

最大动能：

最大应变能：

固有频率：

以上式中φ”在离散体系中无法直接通过运算给出，利用中心差分法求近似解，对第2至第n个节点有：

令第1个节点φ₁”＝φ₂”，第n+1个节点φ_n”＝φ_n+1”，当划分段数足够多时这种近似带来的误差微乎其微。

步骤四、校验所得结果是否满足收敛条件：

计算是否收敛通过固有频率是否满足收敛条件来判断，在整个计算过程中所得频率最小值记为ω_min，最后一次计算所得频率值为ω_j，设定可接受的最大相对误差限为ε；

通过平衡方程式：

ε′＝|ω_min-ω_j|/ω_min

计算其相对误差限ε’，若满足ε’≤ε则进行步骤五，若ε’＞ε则进行步骤二；

该判断方法有极小概率造成结果失真，从理论上讲固有频率最小值也是最接近结构真实自然频率的，但在获得较为接近结构真实频率前若蒙特卡洛方法在随机选取模态函数的过程中有两次选取的极为接近满足收敛条件且为计算过程中的最小值时，将直接输出错误的计算结果。对此，进行多次运算或是提高收敛条件的严苛性均能极大程度地减小失真的可能性。一般来讲由此产生的失真均发生在计算的初期，所得结果误差很大，也可根据经验直接判断出来。

步骤五、得到体系的固有频率和模态函数：

输出各节点坐标，根据坐标绘制梁结构第一阶振型并进行归一化处理，至此即可得到最终的固有频率和振动模态等参数。

实施例2

本实施方式提供了一种基于蒙特卡洛法的梁结构固有频率和振型计算方法的MATLAB编程实现，包括如下步骤：

步骤一、设定初始参数，以图2所示锥形悬臂梁为例：

1.1、锥形梁的几何属性设定，限于篇幅只给出部分参数：

梁长、梁固定端高/2和划分段数:

L＝100；h＝10；b＝5；n＝1000；

1.2、材料属性设定：

梁材料密度、梁质量微元和梁材料弹性模量

rho＝1；m＝rho*2*h*b*(1-x/L)；E＝3；

1.3、梁惯性矩和截面抗弯刚度:

I＝b*(2*h*(1-x/L)).^3/12；EI＝E*I；

1.4、后续计算中初始参数设定：

固有频率初值和第j次循环中的固有频率:omega＝20；bound＝20；

要求固有频率初值大于第一阶固有频率，可取足够大

步骤二、选择合适的函数类型，随机生成初始模态函数：

2.1、利用MATLAB中的随机数生成器生成系数

aa(1)＝2*rand(1)-1；aa(2)＝2*rand(1)-1；aa(3)＝2*rand(1)-1；aa(4)＝2*rand(1)-1；％随机生成介于-1到1之间的随机数作为模态函数中的系数；

2.2、选择合适的满足边界条件的振型函数，并将系数带入。

phi1＝a(1)*(x/L).^2+a(2)*(x/L)^2.*(1-x/L)+a(3)*(x/L)^2*(1-x/L).^2+a(4)*(x/L)^2.*(1-x/L).^3；％悬臂梁振型函数，在此以多项式形式来表达；

步骤三、求解体系的动能、应变能和固有频率：

3.1、先计算离散化后各点的弯矩。

for i＝2:n

aa(i)＝(ph i 1(i-1)+ph i 1(i+1)-2*ph i 1(i))/dx^2；

c(i)＝E I(i)*(aa(i))^2*dx；

End％第2至第n个点的弯矩；

c(1)＝c(2)；c(n+1)＝c(n)；％第1和第n+1个点的弯矩

3.2、求解体系的最大应变能。

d＝sum(c)；％求解应变能，已完成约分化简

3.3、求解体系的最大动能。

e＝sum(dx*m.*(phi1).^2)；％求解动能，已完成约分化简；

3.4、求解固有频率

w＝sqrt(d/e)；％瑞利法求频率

步骤四、校验所得结果是否满足收敛条件：

4.1、确定最小频率和待验证收敛性的频率；

4.2、判断固有频率是否满足收敛条件

在此相对误差限取为0.001；

if abs((omega-bound)/omega)<＝0.001

break

end

步骤五、得到体系的固有频率和模态函数：

5.1给出梁的固有频率并绘制梁的模态曲线。

omega

figure

plot(x,mod)。

实施例3

本实施方式以高度线性变化的锥形悬臂梁和等截面简支梁为算例，给出了基于蒙特卡洛法的梁结构固有频率和振型计算的过程和结果，包括如下步骤：

步骤一、设定初始参数

锥形梁的参数设定同实施例2，简支梁参数设定：梁长150；每延米质量10；截面抗弯刚度EI＝20000，其余参数同上。

步骤二、选择合适的函数类型，随机生成初始模态函数

简支梁的初始模态函数为：

步骤三、求解体系的动能、应变能和固有频率

按离散化公式进行求解，与前文相同

步骤四、校验所得结果是否满足收敛条件

步骤五、给出体系的固有频率和模态函数

简支梁算得的最终结果如图3(第一次计算所得振型)和下表所示，简支梁振型函数理论解易得，图中给出对比。部分振型函数难以求得理论解的梁给出程序多次计算得到的结果，以验证结果的准确性和稳定性。

表1简支梁工况下的固有频率

第十次计算结果即为前文提到的失真情况，直接忽略即可。

表2悬臂梁工况下的固有频率，几何参数与简支梁相同，振型函数如图4所示：

表3锥形悬臂梁工况下的固有频率，几何参数与具体实施方式二相同，振型函数如图5所示：

以上所述，仅为本公开的具体实施方式，但本公开的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本公开揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本公开的保护范围之内。因此，本公开的保护范围应以权利要求的保护范围为准。

Claims

1.一种基于蒙特卡洛法的梁结构固有频率和振型计算方法，其特征在于，包括如下步骤：

步骤一,设定初始参数

(1)几何属性参数和材料属性参数；

(2)将梁沿长度方向划分成等长的n段，并在各节点处生成坐标元素，合成坐标向量x；

步骤二，选择合适的函数类型，随机生成初始模态函数

(1)随机生成初始模态函数系数向量；

(2)选择满足计算体系边界条件的坐标函数，与坐标函数合成模态函数；

步骤三，求解体系的动能、应变能和固有频率

(1)计算各节点的横向位移；

(2)利用位移计算体系的最大动能T_max和利用位移计算体系的最大应变能V_max；

(3)通过瑞利法求解体系固有频率ω；

步骤四，校验所得结果是否满足收敛条件

步骤五，得到体系的固有频率和模态函数

2.根据权利要求1所述的基于蒙特卡洛法的梁结构固有频率和振型计算方法，其特征在于，所述步骤一中几何参数包括梁长L、宽b、高h、截面积A、截面惯性矩I和划分段数n；材料属性参数包括材料弹性模量E、密度ρ和梁每延米的质量m。

3.根据权利要求1所述的基于蒙特卡洛法的梁结构固有频率和振型计算方法，其特征在于，所述步骤二中模态函数的维数等同于组合为模态函数的坐标函数个数。

4.根据权利要求1所述的基于蒙特卡洛法的梁结构固有频率和振型计算方法，其特征在于，所述步骤五中在计算梁的第一个和最后一个节点应变能时，第1个节点的应变能V₁按照式V₁＝V₂计算,第n+1个节点的应变能V_n+1按照式V_n+1＝V_n计算。

5.根据权利要求4所述的基于蒙特卡洛法的梁结构固有频率和振型计算方法，其特征在于，计算单元的刚度S时，第k个单元的刚度按照式S_k＝E_k·A_k/lo_k计算。