【JY】结构动力学初步-单质点结构的瞬态动力学分析

最新推荐文章于 2022-08-19 17:01:18 发布

建源之光

最新推荐文章于 2022-08-19 17:01:18 发布

阅读量1.6k

点赞数 1

文章标签：人工智能数据分析数学建模编程语言信号处理

本文链接：https://blog.csdn.net/Gavinson/article/details/116021283

版权

嘿朋友~记得先点蓝字关注我哦～

简介

单质点体系振动是最为简单的振动，通常在学习结构动力学中也是最开始学习这部分的知识和内容，这部分内容最为基础，也非常重要。它包括单自由度体系振动分析中涉及的物理量和基本概念，而且实际运动中，许多的问题也可按单自由度体系计算，比如普通的隔震结构、多自由度在正则化坐标系下的各个自由度均为解耦的单自由度体系。单质点的动力特性在隔震设计中起到指导性的作用，因此获取可靠准确的单质点结构分析结果十分重要，工程中最常用的3款软件为SAP2000、OpenSees、ANSYS，文章在对理论介绍后介绍了三种软件建模过程，起到对结构动力学学习的参考作用。

理论分析：

当结构为单质点自由振动体系时候，通过达朗贝尔原理或拉格朗日方程均可列得到下式：

如果设

可得到

补充：在隔震结构中，若上部结构和含隔震层时的结构周期差别较大时，上部结构可当成一个刚体，可视为单质点的剪切型模型分析。假设一个2自由度（带隔震层的结构和纯上部结构）隔震结构的等效地震反应分析。

隔震结构的周期：

上部结构的周期

图中△为：

当周期比值

很小的时候，上部结构可以近似为一个刚体，即为单自由度模型，具体的下回再具体解说。

接下来讲解下四种方法：杜哈梅积分（分段解析法）、纽马克β法、威尔逊θ法、Hilber-Hughes-Taylor a（HHT）方法。

GO!GO!GO!

杜哈梅积分：

在实际工程中，很多动力荷载不是简谐荷载，也不是周期荷载，而是随时间任意变化的荷载，此时可采用的计算方法是杜哈梅积分法。强迫运动方程是线性的，可以运用叠加原理。体系在随时间任意变化的动力荷载作用下的响应，可视作一系列独立瞬时冲量连续作用下响应的总和。因此，只需对瞬时冲量作用引起的微分响应进行积分，便可得到一般动力荷载作用下的响应。如下图所示：

通过冲量可以得到该式子：

即：

式中，dv为瞬时冲量引起的速度增量。此时质体的位移增量可由上式积分求得，它是时间的二阶微量，可以略去。因此，瞬时冲量作用所引起的微分响应dy（t）即为以dv为初始速度的微幅自由振动，可得：

如果冲量是在t=x时作用在体系上的。则在冲量结束以后的任一时刻，则位移表达式为：

于是，在任意荷载下的位移有：

当考虑阻尼影响时，杜哈梅积分位移计算式为：

进而，通过数值解法，如梯形公式、辛普森公式等，即可求解上述方程，可得到结果，也可通过对荷载进行分段求解分段解析解组合而成，即为分段解析法。

纽马克β法：

为了方便大家理解，小编此处仅罗列纽马克β法的步骤，详情可以翻阅《结构动力学》

（1）初始计算

①组成该体系的刚度矩阵[K]、质量矩阵[M]和阻尼矩阵[C]。

②选择时间步长Δt和积分参数r、β；（为了稳定收敛，通常取r=1/2，β=1/4）

③计算积分常数

④形成有效刚度矩阵

⑤对上矩阵进行三角分解

⑥指定初始条件：

（2）对每一时间步长

① 计算时刻t+At的有效荷载1

②求解t+Δt的位移

③计算t+Δt的速度和加速度

转步骤（2）的①，t=t+Δt，迭代循环，即可求得结果。

威尔逊θ法：

Wilson-θ法是常规的Newmark方法通过引入一个系数θ达到了无条件的稳定系数的提出是因为观察到一种现象，即不稳定的解趋向于在真实解附近振荡。因此，如果在时间增量内计算数值解，为使这种振荡减到最小，可以通过如下定义的时间步及荷载对Newmark方法作简单的修改来完成，即

（其中θ>=1.0）

对于Newmark方法进行修正后可得到下式子：

θ系数的使用有助于在体系的高阶振型中去除数值阻尼。如果θ等于1.0，就是Newmark方法。然而，对于高阶振型响应是很重要的问题，引入的误差可能比较大。除此之外，还不能在时间t处精度满足动力平衡方程。因此，若是存在高阶阵型的问题，则不再推荐使用系数θ。

Hilber-Hughes-Taylor a(HHT)方法：

a方法使用Newmark方法求解下列修正的运动方程：

其中a的取值参考范围在0~1/13。

当a等于零时，此方法还原为常量加速度方法。它在高阶振型中产生数值能量损耗。然而，它不能像用刚度比例阻尼那样的阻尼比对其进行预测。同时，它也不能在时间处求解基本平衡方程。然而，当前许多计算机程序都在使用它。该方法的效果似乎与使用刚度比例阻尼法的效果很相近。建议使用默认的HHT方法，除非对其他方法有特定的需要。.在非线性分析中，经常需要使用一个a负值来确保结果的收敏性a。

数值分析算例

【例题】为了对比各软件与理论计算结果的误差和展示建模过程，以一个简单的单质点体系作为算例，地震波选取三条人工波，单质点体系的自振周期为2.7s，阻尼比为0.05，提取质点的加速度及位移时程进行对比。

地震波时程

地震波反应谱

（为了对比简洁以下结果仅呈现RGB1的计算结果）

SAP2000

Sap2000模型参数设置：

SAP四种求解方法结果对比：

SAP2000加速度时程结果

SAP2000位移时程结果

ANSYS

ANSYS单质点模型

ANSYS命令流：

/prep7
!设置单位为N/m/kg
/units,si
!建立节点
n,1,0,0,0
n,2,0,1,0
!定义质量单元
et,1,mass21,
r,1,0,1846,0
!定义弹簧单元
et,2,combin14
r,2,10000
!建立质量单元
type,1
real,1
e,2
!建立弹簧单元
type,2
real,2
e,1,2
!约束自由度
d,1,all
d,2,all
ddel,2,uy
finish
!瞬态动力学求解
/SOLU   
ANTYPE,trans
TRNOPT,FULL
!阻尼比0.05
ALPHAD,0.08976
BETAD,0.0274
!时间步及间隔
t=6000
ti=0.005
!定义地震波数组
*DIM,wavey,ARRAY,t,1,1
!导入地震波数据
*create,readwave
*vread,wavey(1),'RGB1','txt',' ',
(e20.8)
*end
/input,readwave
!分步求解
*do,i,1,t,1   
time,ti*i
kbc,1
acel,,wavey(i),
nsel,all
OUTRES,all,all
solve   
*enddo
finish

ANSYS质量单元为MASS21单元，弹簧单元为COMBIN14单元，单质点模型周期为2.7s，阻尼为瑞丽阻尼。

ANSYS加速度时程结果

ANSYS位移时程结果

OpenSees

Opensees的代码如下：

####### ================ 数据包录入 ====================== ########
from openseespy.opensees import *
import numpy as np
import matplotlib.pyplot as plt
## =========== 参数 ============ ##
def Link(m,damp,k,wave,dt,Amp,g):
    wipe()
    model('basic', '-ndm', 2, '-ndf', 2)
    node(1, 0, 0)
    node(2, 0, 1)
    mass(2, m, m, 0)
    fix(1, 1, 1)
    uniaxialMaterial('Elastic',1,k)
    element("twoNodeLink", 1 ,1 ,2 ,'-mat' ,1 ,'-dir', 2)
    wavef = np.mat(np.loadtxt(wave)).T #地震波数据读入
    nwave = len(wavef)
    timeSeries('Path', 1, '-filePath', wave, '-dt', dt, '-factor', g)
    pattern('UniformExcitation' , 1, 1, '-accel' , 1, '-fact', Amp)
system('BandSPD')
numberer("RCM")
    constraints("Plain")
    integrator("LoadControl", 1.0)
    test('NormDispIncr', 1.0e-15,  10 )
    algorithm("Linear")
    #integrator('CentralDifference' )
    #algorithm('NewtonLineSearch')
    #algorithm('SecantNewton')
    integrator('Newmark',  0.5,  0.25 )
    #integrator('HHT',  1.0,  0.5, 0.25 )
    analysis('Transient')
    C1 = 2*damp*(k/m)**0.5
    rayleigh(C1, 0, 0, 0)
    tFinal = (nwave-1) * dt
    tCurrent = getTime()
    ok = 0
    time = []
    Acce = []
    Dise = []
    Foree = []
    while ok == 0 and tCurrent < tFinal:
        ok = analyze(1, 0.005)
        if ok != 0:
            test('NormDispIncr', 1.0e-12,  100, 0)
            algorithm('ModifiedNewton', '-initial')
            ok =analyze( 1, 0.005)
            if ok == 0:
                print("ok")
            test('NormDispIncr', 1.0e-12,  10 )
            algorithm('Newton')  
        tCurrent = getTime()
        time.append(tCurrent)
        Acce.append(nodeAccel(2,1))
        Dise.append(nodeDisp(2,1))
        Foree.append(eleDynamicalForce(1,1))
    Acc = np.mat([Acce]).T
    Acc = Acc + wavef*Amp*g
    Dis = np.mat(Dise).T*(-1)
    Fore = np.mat(Foree)
    return time,Acc,Dis,Fore

为了方便看清对比，将计算结果以1.5m/s2进行相上偏移对比（均以Newmark-β法计算）

OpenSees加速度时程结果

为了方便看清对比，将计算结果以100mm进行相上偏移对比（均以Newmark-β法计算）

OpenSees位移时程结果

Matlab

理论计算均利用专业软件MATLAB完成，计算方法包括分段解析法、杜哈梅积分法、威尔逊-q法、纽马克法。

对比结果：

数值软件与理论计算加速度时程对比结果

数值软件与理论计算速度时程对比结果

数值软件与理论计算位移时程对比结果

数值软件与理论计算位移时程偏移对比结果

结论

各算法下，求解结果高度一致，且误差几乎可以忽略。其中可以发现，加速度/速度/位移在算法中若是通过每一步迭代出来时，即在每个时刻均可通过相应算法得到加速度/速度/位移时，则加速度结果未呈现由于后期求导而产生的背景噪声。若是求解方程的位移结果，则通过位移结果进行求导得到的加速度/速度，会出现较严重的背景噪声（可能由于信号求导产生的高阶误差），且该特性在长周期单质点更为明显。

其中SAP2000为步步算得加速度/速度/位移，而Opensees和ANSYS则采用求导形式，在接下来的文章中会介绍由位移时程获取准确加速度时程的方法，赶紧点个关注吧，这样就不会错过后续的文章啦！

【JY】从一根悬臂梁说起

【JY】2B青年欢乐多之Matlab篇

【JY】Abaqus6.14-4如何关联fortran？

【JY】位移角还是有害位移角？

【JY】如何利用python来编写GUI？

【JY】施工的竖向变形分析在分析什么？

【JY|土木】失稳你过来，我们谈谈吧。

【JY】OpenSEES简介及分析流程

【JY】如何解决MATLAB GUI编程软件移植运行问题？

~求关注转赞评~

建源之光

关注

1
点赞
踩
19

收藏

觉得还不错? 一键收藏
0
评论
【JY】结构动力学初步-单质点结构的瞬态动力学分析

嘿朋友~记得先点蓝字关注我哦～简介单质点体系振动是最为简单的振动，通常在学习结构动力学中也是最开始学习这部分的知识和内容，这部分内容最为基础，也非常重要。它包括单自由度体系...
复制链接

扫一扫