哈工大结构动力学作业-威尔逊-θ法
哈工大机械原理大作业——连杆机构运动分析16___2014

Harbin Institute of Technology机械原理大作业——连杆机构运动分析课程名称:机械原理院系:能源科学与工程学院班级:完成者:学号:题号: 16任课教师:丁刚完成内容:在完成题目计算要求的同时,扩展了内容,程序为该结构的通用程序,可解决机构在不同条件下的运动情况,文本最末为几种情况的分析哈尔滨工业大学16、如图所示机构,已知机构各构件的尺寸为,试求构件5的角位移、角速度和角加速度,并对计算结构进行分析。
(1)、结构分析从侧面看原机构为此机构分为级杆组(原动件1),级杆组RRP(2号套筒、3号杆),级杆组RRP(4号套筒、5号杆)(2)、建立坐标系(3)、各个杆组的运动分析采用逆推法,从RRP杆组(4号套筒、5号杆)开始分析已知,,,,现在假定已知,,其中,,,即消去,可得可求得,也可以通过书上3-23式求得通过正弦定理可求得再来看看角速度关系对于加速度,有如下关系其中到此4、5杆就分析完毕了,别忘记之前的假设,我假设了已知,,为求,,,现在来分析RRP杆组(2号套筒、3号杆)已知,,,已知,,,,其中,,,即消去,可得反解,即可求得,也可以通过书上3-23式求得通过正弦定理可求得继续,我们来看看角速度关系对于加速度,有如下关系其中现在,只需将所求得的,,和,,关联起来这是同一根杆,,,现在来看,,,由题目得,,和是未知的,但不影响整体,不然给一个初值,,当然,这是可以随意更改的。
基于以上的基本原理,matlab R2012b程序如下syms theta theta1 theta2 lamuda lamuda1 lamuda2 sigma sigma1 sigma2 beta beta1 beta2 l1 l11 l2 l21 t output itheta1=10;theta2=0;i=0;for theta3=60:420theta=theta3/180*pi;beta=asin((100/200)*sin(theta))+theta;l1=0.2*sin(beta)/sin(theta);beta1=(-theta1*(l1*sin(theta))*sin(theta)+theta1*(l1*cos(theta))*cos(theta))/(0.2*(sin(theta)*sin(b eta)+cos(theta)*cos(beta)));l11=-(theta1*(l1*sin(theta))*l1*cos(beta)+theta1*(l1*cos(theta))*l1*sin(beta))/(0.2*(sin(theta)*si n(beta)+cos(theta)*cos(beta)));C=(theta1^2)*0.2*cos(beta)-theta2*l1*sin(theta)-(theta1^2)*l1*cos(theta)-2*l11*theta1*sin(theta) ;D=(theta1^2)*0.2*cos(beta)+theta2*l1*sin(theta)-(theta1^2)*l1*cos(theta)+2*l11*theta1*sin(thet a);beta2=(-C*sin(theta)+D*cos(theta))/(0.2*(sin(theta)*sin(beta)+cos(theta)*cos(beta)));lamuda=beta-pi/2;lamuda1=beta1;lamuda2=beta2;sigma=asin((100/200)*sin(lamuda))+lamuda;l2=0.2*sin(sigma)/sin(lamuda);sigma1=(-lamuda1*(l2*sin(lamuda))*sin(lamuda)+lamuda1*(l2*cos(lamuda))*cos(lamuda))/(0.2 *(sin(lamuda)*sin(sigma)+cos(lamuda)*cos(sigma)));l21=-(lamuda1*(l2*sin(lamuda))*l2*cos(sigma)+lamuda1*(l2*cos(lamuda))*l2*sin(sigma))/(0.2* (sin(lamuda)*sin(sigma)+cos(lamuda)*cos(sigma)));A=(lamuda1^2)*0.2*cos(sigma)-lamuda2*l2*sin(lamuda)-(lamuda1^2)*l2*cos(lamuda)-2*l21*la muda1*sin(lamuda);B=(lamuda1^2)*0.2*cos(sigma)+lamuda2*l2*sin(lamuda)-(lamuda1^2)*l2*cos(lamuda)+2*l21*l amuda1*sin(lamuda);sigma2=(-A*sin(lamuda)+B*cos(lamuda))/(0.2*(sin(lamuda)*sin(sigma)+cos(lamuda)*cos(sigma )));i=i+1;output(i,1)=fix(theta/pi*180);output(i,2)=fix(sigma/pi*180);output(i,3)=fix(sigma1);output(i,4)=fix(sigma2);endoutputa=output(:,1);b=output(:,2);c=output(:,3);d=output(:,4);h1=plot(a,b);hold on;h2=plot(a,c);hold on;h3=plot(a,d);hold on;set(h1,'color',[1 0 0],'linewidth',2);set(h2,'color',[0 1 1],'linewidth',1);set(h3,'color',[0 0 1],'linewidth',2);m=legend('角位移','角速度','角加速度');x label('θ');title('平面连杆机构运动分析');figure;h1=plot(a,b);hold on;x label('θ');ylabel('角位移');title('平面连杆机构运动角度——角位移图');figure;h2=plot(a,c);hold on;x label('θ');ylabel('角速度');title('平面连杆机构运动角度——角速度图'); figure;h3=plot(a,d);hold on;x label('θ');ylabel('角加速度');title('平面连杆机构运动角度——角加速度图');汇总图各自的图像结果分析,上面的图形只是在一个初值,的条件下得出的,为了能解决所有问题,修改程序如下syms theta theta1 theta2 lamuda lamuda1 lamuda2 sigma sigma1 sigma2 beta beta1 beta2 l1 l11 l2 l21 t output iprompt={'输入:', '输入' ,'输入' };%设置提示字符串name='输入初值';%设置标题 numlines=1;%指定输入数据的行数 defAns={'60','10','0'};%设定默认值 Resize='on';%设定对话框尺寸可调节answer=inputdlg(prompt,name,numlines,defAns,'on');%创建输入对话框 h= str2num(answer{1}); theta1= str2num(answer{2}); theta2= str2num(answer{3}); i=0;for theta3=h:(360+h) theta=theta3/180*pi;beta=asin((100/200)*sin(theta))+theta; l1=0.2*sin(beta)/sin(theta);beta1=(-theta1*(l1*sin(theta))*sin(theta)+theta1*(l1*cos(theta))*cos(theta))/(0.2*(sin(theta)*sin(b eta)+cos(theta)*cos(beta)));l11=-(theta1*(l1*sin(theta))*l1*cos(beta)+theta1*(l1*cos(theta))*l1*sin(beta))/(0.2*(sin(theta)*si n(beta)+cos(theta)*cos(beta)));C=(theta1^2)*0.2*cos(beta)-theta2*l1*sin(theta)-(theta1^2)*l1*cos(theta)-2*l11*theta1*sin(theta) ;D=(theta1^2)*0.2*cos(beta)+theta2*l1*sin(theta)-(theta1^2)*l1*cos(theta)+2*l11*theta1*sin(thet a);beta2=(-C*sin(theta)+D*cos(theta))/(0.2*(sin(theta)*sin(beta)+cos(theta)*cos(beta)));lamuda=beta-pi/2;lamuda1=beta1;lamuda2=beta2;sigma=asin((100/200)*sin(lamuda))+lamuda;l2=0.2*sin(sigma)/sin(lamuda);sigma1=(-lamuda1*(l2*sin(lamuda))*sin(lamuda)+lamuda1*(l2*cos(lamuda))*cos(lamuda))/(0.2 *(sin(lamuda)*sin(sigma)+cos(lamuda)*cos(sigma)));l21=-(lamuda1*(l2*sin(lamuda))*l2*cos(sigma)+lamuda1*(l2*cos(lamuda))*l2*sin(sigma))/(0.2* (sin(lamuda)*sin(sigma)+cos(lamuda)*cos(sigma)));A=(lamuda1^2)*0.2*cos(sigma)-lamuda2*l2*sin(lamuda)-(lamuda1^2)*l2*cos(lamuda)-2*l21*la muda1*sin(lamuda);B=(lamuda1^2)*0.2*cos(sigma)+lamuda2*l2*sin(lamuda)-(lamuda1^2)*l2*cos(lamuda)+2*l21*l amuda1*sin(lamuda);sigma2=(-A*sin(lamuda)+B*cos(lamuda))/(0.2*(sin(lamuda)*sin(sigma)+cos(lamuda)*cos(sigma )));i=i+1;output(i,1)=fix(theta/pi*180);output(i,2)=fix(sigma/pi*180);output(i,3)=fix(sigma1);output(i,4)=fix(sigma2);endoutputa=output(:,1);b=output(:,2);c=output(:,3);d=output(:,4);h1=plot(a,b);hold on;h2=plot(a,c);hold on;h3=plot(a,d);hold on;set(h1,'color',[1 0 0],'linewidth',2);set(h2,'color',[0 1 1],'linewidth',1);set(h3,'color',[0 0 1],'linewidth',2);m=legend('角位移','角速度','角加速度');x label('θ');title('平面连杆机构运动分析');figure;h1=plot(a,b);hold on;xlabel('θ');y label('角位移');title('平面连杆机构运动角度——角位移图');figure;h2=plot(a,c);hold on;xlabel('θ');y label('角速度');title('平面连杆机构运动角度——角速度图');figure;h3=plot(a,d);hol d on;xlabel('θ');y label('角加速度');title('平面连杆机构运动角度——角加速度图');这样,在运行程序时就会弹出一个如下图所示的对话框,可以任意给定初值,解决不同问题。
哈工大结构动力学第一讲概述

加分情况
• • • • • 独到的作业形式与新颖的解题方法; 对老师教学与学生学习的好建议; 高质量的CAI与自学报告; 好的问题。 加分幅度视情况而定;
几点说明
• • • • • 按时保质保量完成作业, 双周交作业. 答疑时间: 周二下午 13:30 – 17:00 联系电话: 86283199 电子信箱: guangchun_zhou@ 办公室: 土木工程学院509室.
• 纸面作业: 20%、CAI作业: 15%.
作业做法: (1) 对各章全部作业题进行分类分析, 自 主选择各个类别应做题数; (2) 不批改对错, 对错自 行检验, 答案可以与老师和助教探讨; (3) 可以改造 、自编作业题(酌情加分);(4)CAI作业可以自 行选择分析内容,但指定的须完成;(6)自学内 容交报告。
违纪严重者,视情况取消考试资格
结构动力学入门问答
请简要回答下面问题: 1. 结构动力问题的特性体现在哪里(对比 静力问题)? 2. 怎样解决结构动力问题? 3. 你认为学习结构动力学能带给你什么? 4. 你认为学好结构动力学最重要的什么?
Turkey_17-8-1999
India_26-1-2001
参考教材
• <结构动力学> 第二版(修订版), R.W.克拉 夫著 王光远译, 高等教育教育出版社, 2006. • <结构动力学> Anil K.Chopra, 谢礼立, 吕 大刚等译, 高等教育教育出版社, 2007. • <结构动力学> 邹经湘主编 • <建筑结构振动计算续编> 郭长城编著
成绩组成
结构:建筑物中承受荷载起骨架作用的部分.几何性质 ,材料性质,支撑连接方式. 动荷载:大小,方向,作用点随时间变化的荷载.其引 起的结构中的质量效应(惯性力)不可忽视时. 响应:位移,速度,加速度,内力,应力,应变等. 规律:响应与荷载及结构性质的定量关系.
中心差分法、纽马克法和威尔逊-θ法与精确解的误差分析

中心差分法、纽马克法和威尔逊-θ法与精确解的误差分析作者:于津津贾慧敏宋敏来源:《教育周报·教育论坛》2018年第05期摘要:在动荷载作用下的物体位移、速度和加速度的计算中,中心差分法、纽马克法和威尔逊-θ法三类方法都是可取的,为结构动力学的理论研究提供了参考。
但三类方法与精确值之间均存在一定的误差,本文基于这一问题进行研究和计算,通过图表展示这三类方法与精确值之间的关系。
关键词:结构动力学;中心差分法;纽马克法;威尔逊-θ法一、引言:结构动响应的数值计算问题,主要针对多自由或者连续体经过空间散离后建立的二阶常微分方程组形式的运动控制方程:[M] {¨x}+[ C] {﹒x}+[ K ]{x}=Q ;;(1)为了探究三种方法相较于精确解的误差,用如下具体问题进行具体分析。
如图1所示,该体系在冲击荷载 p(t)=[0 10]T 作用下,求该体系的位移反应表达式,质量单位Kg,弹簧k单位N/cm。
另:自由振动的周期T1=4.45,T2=2.8,使用中心差分法计算,取时间步长Δt=0.1,T2=0.28,并假定X0=0;V0=0试计算这个系统在前12个时间步长的反应。
取δ=0.25,γ=0.5,用纽马克法计算该系统的动力反应。
取θ=1.4,用威尔逊-θ法计算该系统的反应。
二、计算方法简介1、精确解计算根据精确解的计算公式可得:X1(t)=1-5/3×cos(2^0.5×t)+2/3×cos(5^0.5×t)X2(t);=3-5/3×cos(2^0.5×t)-4/3×cos(5^0.5×t)速度的计算公式为位移的导数,加速度的公式为速度的导数。
(下同)2、中心差分法用位移向前一步的差分表示的速度后一步的差分表示的速度的平均来确定当前时刻的速度,得到以当前时刻t为中心的前后时刻位移的差分表示的速度,即:若:x=x0-1/(2×a1)×d x0+1/(2×a0)×d2x0; ;x1(t)=x0(1);x2(t)=x0(2);3、纽马克法当在t时刻的响应{x}t,{﹒x}t,{¨x}t,已知时,要求下一时刻t+Δt的响应值{x} t+Δt,{﹒x} t+Δt,{¨x} t+Δt,令在待求时刻动力学方程成立,即:{﹒x} t+Δt={﹒x}t+Δt(1-γ){¨x} t +γΔt{¨x} t+Δt ;;(2){x} t+Δt={x}t+{﹒x}tΔt+(0.5-δ){¨x} tΔt^2+δ{¨x} t+ΔtΔt^2 ;(3)β,γ为按积分精度和稳定性要求而确定的参数,由式3可解得:1/{¨X}t+Δt;=βΔt 2({X}t+Δt -{x}t)-βΔt ×1/{﹒x}t-(2β-1) ×1/{¨x}t ;;(4)将(4)带入(2)得:{﹒x}t+Δt =γ/βΔt 2×({x}t+Δt -{x}t)+(1 –γ/β){﹒x}t +(1 -1/2β)t{¨x}t ;;(5){x}t +Δt 可由t +Δt 时刻的运动方程求得,即:[M]{¨X}t+Δt +[C]{¨X}t +Δt +[K]{X}t +Δt =[F] t +Δt ;;(6)將式(4)、式(5)代入式(6),可求得求得{X}t+Δt,后求{﹒X}t +Δt 和{¨X}t +Δt。
哈工大结构动力学大作业2012春

结构动力学大作业对于如下结构,是研究质量块的质量变化和在简支梁上位置的变化对整个系统模态的影响。
1以上为一个简支梁结构。
集中质量块放于梁上,质量块距简支梁的左端点距离为L.将该简支梁简化为欧拉伯努利梁,并离散为N 个单元。
每个单元有两个节点,四个自由度。
单元的节点位移可表示为:]1122,,,e v v δθθ⎡=⎣则单元内一点的挠度可计作:带入边界条件:1332210)(x a x a x a a x v +++=01)0(a v x v ===3322102)(L a L a L a a v L x v +++===110d d a xv x ===θ2321232d d L a L a a xv Lx ++===θ10v a =[]1234N N N N N =建立了单元位移模式后,其动能势能均可用节点位移表示。
单元的动能为:00111()222l l T TT ke e e e e y E dx q N Ndxq q mq t ρρ∂===∂⎰⎰ 其中m 为单元质量阵,并有:lT m N Ndx ρ=⎰带入公式后积分可得:2222156225413224133541315622420133224l l l l l l l m l l ll ll ρ-⎡⎤⎢⎥-⎢⎥=⎢⎥-⎢⎥---⎣⎦单元势能可表示为2220011()()222T l lT Te pe e e e q y E EI dx EI N N dxq q Kq x ∂''''===∂⎰⎰其中K 为单元刚度矩阵,并有()lT K EI N N dx ''''=⎰2232212612664621261266264l l l l l l EI k l l l l lll -⎡⎤⎢⎥-⎢⎥=⎢⎥---⎢⎥-⎣⎦以上为单元类型矩阵,通过定义全局位移矩阵,可以得到系统刚度矩阵和系统质量矩11θ=a )2(1)(3211222θθ+--=Lv v L a )(1)(22122133θθ++-=Lv v L a 1232133222231)(θ⎪⎪⎭⎫ ⎝⎛+-+⎪⎪⎭⎫ ⎝⎛+-=L x L x x v L x L x x v 22232332223θ⎪⎪⎭⎫ ⎝⎛-+⎪⎪⎭⎫ ⎝⎛-+L x L x v L x L x 24231211)()()()()(θθx N v x N x N v x N x v +++=阵。
哈尔滨工业大学 结构力学II 第二套张金生 结构动力学-9

X 2
1 1.78 2.21 1 1.8 2.24
X DX
3
2
X 3
2.算例: 用迭代法计算图示体系的基频和基本振型.
m m m
解:
m m m m
1 1 1 1 2 2 1 k 1 2 3
X X a
~ X
0 0
T X 1 mX 0 0 X X 1 *
1
4.667 m 8.334 归一化 k 10.334 4.99 m 8.98 归一化 k 11.19
X 2
X DX
3
2
X 3
2.算例: 用迭代法计算图示体系的基频和基本振型.
m m m
解:
m m m m
y(t ) X i i cos( i t i )
动能为
y2 (t )
速度为
m1
y1 (t )
1 1 1 2 2 2 Ti (t ) m1 y1 (t ) m2 y2 (t ) mN y N (t ) 2 2 2 1 T y (t ) m y (t ) 2 1 T X i mX i i2 cos2 ( i t i ) 2 势能为 1 T U i (t ) X i k X i sin 2 ( i t i ) 2
a 0.0328 k / m b 0.0591 m / k 1 2 3 (a bቤተ መጻሕፍቲ ባይዱ 3 ) 0.0624 2 3
m
k
m m m m
2 1 0 k 1 2 1 k 0 1 1 0 0.151 0.0591 c am bk 0.0591 0.151 0.0591 mk 0 0.0591 0.0919
结构动力学中的常用数值方法

第五章 结构动力学中的常用数值方法5.1.结构动力响应的数值算法....0()(0)(0)M x c x kx F t x a x v ⎧++=⎪⎪=⎨⎪=⎪⎩当c 为比例阻尼、线性问题→模态叠加最常用。
但当C 无法解耦,有非线性存在,有冲击作用(激起高阶模态,此时模态叠加法中的高阶模态不可以忽略)。
此时就要借助数值积分方法,在结构动力学问题中,有一类方法称为直接积分方法最为常用。
所识直接是为模态叠加法相对照来说,模态叠加法在求解之前,需要对原方程进行解耦处理,而本节的方法不用作解耦的处理,直接求解。
(由以力学,工程中的力学问题为主要研究对象的学者发展出来的)中心差分法的解题步骤1. 初始值计算(1) 形成刚度矩阵K ,质量矩阵M 和阻尼矩阵C 。
(2) 定初始值0x ,.0x ,..0x 。
(3) 选择时间步长t ∆,使它满足cr t t ∆<∆,并计算 021()a t =∆,112a t=∆,202a a =(4) 计算...0011122t x x x x a a -∆=-+(5) 形成等效质量阵01M a M a C -=+ (6) 对M -阵进行三角分解T M LDL -= 2.对每一时间步长(1) 计算时刻t 的等效载荷21()()t t t tt Q Q K a M x a Ma C x --∆=---- (2) 求解t t +∆时刻的位移 ()Tt t t L D L x Q -+∆=(3) 如需要计算时刻t 的速度和加速度值,则.1()t t t t t x a x x +∆-∆=-..0(2)t t t t t t x a x x x +∆-∆=-+若系统的质量矩阵和阻尼矩阵为对角阵时,则计算可进一步简化。
纽马克法的解题步骤1.初始值计算(1)形成系统刚度矩阵K ,质量矩阵M 和阻尼矩阵C(2)定初始值0x ,.0x ,..0x 。
(3)选择时间步长t ∆,参数γ、σ。
结构动力学多自由度线性体系Wilson-θ法程序编写

多自由度线性体系Wilson -θ法程序编写【摘要】本文主要介绍了通过使用Matlab 软件,Wilson-θ法编写多自由度线性体系的程序的原理、流程图、具体算例以及使用注意事项。
通过该程序可以得到剪切型结构在任意函数荷载作用下各质点的位移函数。
【关键词】Matlab ;多自由度;Wilson-θ法1.wilson-θ法原理wilson-θ法中最主要的步骤就是推导由t 时刻的状态求t t ∆+时刻的状态的递推公式,现推导如下:对τ积分解出代入整理,得其中本程序的核心就是对以上公式的循环使用。
{}{}{}{})(t t t t t y y tyy -∆+=∆++θτθτt ∆=θτ{}{}{}{}{})(22t t t t t t yy t y y y-∆++=∆++θτθττ{}{}{}{}{}{})(6232t t t t t t t yy ty y y y -∆+++=∆++θτθτττ{}{}{}{}{})(21t t t t t t t yy t y t y y -∆+∆+=∆+∆+θθθθ{}{}{}{}{})2(6)(2t t t t t t t yy t y t y y +∆+∆+=∆+∆+θθθθ{}{}{}{}{}t t t t t t t y y t y y t y 26)()(62-∆--∆=∆+∆+θθθθ{}{}{}{}{}t t t t t t t yty y y t y 22)(3∆---∆=∆+∆+θθθθ[]{}[]{}[]{}{}t t t t t t t t P y k y C ym ∆+∆+∆+∆+=++θθθθ []{}[]{}[]{}{}P y k y C ym =++ []{}[]R y k t t =∆+θ[][][][]c tm t k k ∆+∆+=θθ3)(62[]{}{}{}[]{}{}{}[]{}{}{})223()26)(6()(2t t t t t t t t t t yt y y t c y y t y t m P P P R ∆++∆++∆+∆+-+=∆+θθθθθ{}{}{}{})(t t t t t t P P P P -+=∆+∆+θθ2.程序流程图求出各常数值For I=1 to n[][][][]c a m a k k 1++=3.具体应用算例如图所示,两自由度框架结构,其中初始静止,求各层位移。
哈尔滨工业大学结构动力学PPT课件

x0 x0 , x0 x0 xt c1n cosnt c2n sinnt
c1 x0 n , c2 x0
第36页/共42页
x
t
x0
n
sin nt
x0
cos nt
令
x0 cos n
, x0 sin
则可化为
其中:
xt sinnt
2
x02
x0
n
tg x0n arctg x0n
T1
1 2
l 0
d
l
2
x2
1 2
(1 3
l)x2
1 m1 23
x2TΒιβλιοθήκη T1Tm1 2
m1 3
m
x2
1 2
meq x2
又因为: 弹簧的势能与弹簧质量无关, 则
V 1 kx2 2
由能量法,可得
meq x kx 0 弹性元件质量不能忽略时,利用等
效质量,将质量折算到质量块上, 弹性元件仍看作无质量的。
• 18世纪线性振动理论成熟期。
第11页/共42页
• 19世纪非线性振动理论,各种工程实际结构振动的近似 求解方法。
• 20世纪50年代初由于航空航天工程的发展,原本确定性 理论无法解释包含随机变化的工程问题,发展了随机振 动理论。
• 20世纪后期计算机技术的飞速发展,数值计算方法和理 论成为主要研究方法之一。
第7页/共42页
三、结构动力学研究的内容
结构动力学就是研究结构系统在激励力作用下产生的响 应规律的科学,研究激励力、结构和响应三者关系的科 学。
现代结构动力学主要研究以下三个方面的内容 第一类问题:响应分析(结构动力计算)
输入 (动力荷载)
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
结构动力学大作业(威尔逊- 法)
姓名:
学号:
班级:
专业:
威尔逊—θ法原理及应用
【摘要】在求解单自由度体系振动方程时我们用了常加速度法及线加速度法等数值分析方法。
在多自由度体系中,也有类似求解方法,即中心差分法及威尔逊—θ法。
实际上后两种方法也能求解单自由度体系振动方程。
对于数值方法,有三个重要要求:收敛性、稳定性及精度。
本文推导了威尔逊-θ法的公式,并利用MATLAB 编程来研究单自由度体系的动力特性。
【关键词】威尔逊—θ法 冲击荷载 阻尼比
【正文】威尔逊-θ法可以很方便的求解任意荷载作用下单自由度体系振动问题。
实际上,当 1.37θ>时,威尔逊—θ法是无条件收敛的. 一、威尔逊—θ法的原理
威尔逊-θ法是线性加速度法的一种拓展(当1θ=时,两者相同),其基本思路和实现方法是求出在时间段[],t t t θ+∆时刻的运动,其中1θ≥,然后通过内插得到
i t t +∆时刻的运动(见图 1。
1)。
图 1。
1
1、公式推导
推导由t 时刻的状态求t t θ+∆时刻的状态的递推公式:
{}{}{}{})(t t t t t y
y t y y -∆+=∆++θτθτ
对τ积分
{}{}{}{}{})(22
t t t t t t y
y t y y y
-∆++=∆++θτθττ
{}{}{}{}{}{})(623
2
t t t t t t t y
y t y y y y -∆+++=∆++θτ
θτττ
t ∆=θτ
{}{}{}{}{})(2
1
t t t t t t t y
y t y t y y -∆+∆+=∆+∆+θθθθ
{}{}{}{}{})2(6)(2
t t t t t t
t y
y t y t y y +∆+∆+=∆+∆+θθθθ {}{}{}{}{}t t t t t t t y y t y y t y
26
)()(62
-∆--∆=∆+∆+θθθθ
{}{}{}{}{}t t t t t t t y
t
y y y t y
22)(3∆---∆=∆+∆+θθθθ
[]{}[]{}[]{}{}P y k y C y
m =++ []{}[]{}[]{}{}t t t t t t t t P y k y C y m ∆+∆+∆+∆+=++θθθθ
[]{}
{}t t t
t R y k ∆+∆+=θθ
[][][][]
c t
m t k k ∆+∆+=θθ3)(6
2
[]{}{}
{}[]{}{}{}[]{}{}{})223()26)(6(
)(2t t
t t t t t t
t t
y t
y y t c y y t y t m P P P R ∆++∆++∆+∆+-+=∆+θθθθθ
2、MA TLAB 源程序: clc;clear;
K=input (’请输入结构刚度k (N/m )'); M=input ('请输入质量(kg )');
C=input (’请输入阻尼(N *s/m )'); t=sym (’t ’);%产生符号对象t Pt=input(’请输入荷载);
Tp=input (’请输入荷载加载时长(s)'); Tu=input ('请输入需要计算的时间长度(s ) ’); dt=input ('请输入积分步长(s)'); Sita=input('请输入θ’);
uds=0:dt:Tu;%确定各积分步时刻
pds=0:dt:Tp;
Lu=length(uds);
Lp=length(pds);
if isa(Pt,'sym')%荷载为函数
P=subs(Pt,t,uds); %将荷载在各时间步离散
if Lu〉Lp
P(Lp+1:Lu)=0;
end
elseif isnumeric(Pt)%荷载为散点
if Lu〈=Lp
P=Pt(1:Lu);
else
P(1:Lp)=Pt;
P(Lp+1:Lu)=0;
end
end
y=zeros(1,Lu);%给位移矩阵分配空间
y1=zeros(1,Lu);%给速度矩阵分配空间
y2=zeros(1,Lu);%给加速度矩阵分配空间
pp=zeros(1,Lu-1);%给广义力矩阵分配空间
yy=zeros(1,Lu-1);%给y(t+theta*t)矩阵分配
FF=zeros(1,Lu);%给内力矩阵分配空间
y(1)=input('请输入初位移(m)’);
y1(1)=input(’请输入初速度(m/s)');
%——-—-——-———--———--初始计算-—-—------———————--——--——
y2(1)=(P(1)—C*y1(1)-K*y(1))/M;%初始加速度
FF(1)=P(1)-M*y2(1);
l=6/(Sita*dt)^2;
q=3/(Sita*dt);
r=6/(Sita*dt);
s=Sita*dt/2;
for z=1:Lu—1
kk=K+l*M+q*C;
pp(z)=P(z)+Sita*(P(z+1)—P(z))+(l*y(z)+r*y1(z)+2*y2(z))*M+(q*y(z)+2*y1(z)+s*y2(z))*C;
yy(z)=pp(z)/kk;
y2(z+1)=l/Sita*(yy(z)—y(z))-l*dt*y1(z)+(1-3/Sita)*y2(z);
y1(z+1)=y1(z)+dt/2*(y2(z+1)+y2(zp));
y(z+1)=y(z)+y1(z)*dt+dt*dt/6*(y2(z+1)+2*y2(z));
FF(z+1)=P(z+1)—M*y2(z+1);
end
plot (uds ,y ,’r ’),xlabel('时间 t ’),ylabel('位移 y ’),title ('位移图形’) 二、利用威尔逊-θ法求冲击荷载下的结构反应
1、矩形脉冲
研究不同时长脉冲作用下,体系振动位移。
取单自由度刚度为1N/m ,质量为1/(4*pi^2)kg ,频率为2*pi 1s -,周期为1s ,阻尼c=0,荷载为1N ,积分步长为0。
1,θ=1。
42,初位移为0,初速度为0时的质点位移时间图如下:
图2.1 1/41/4d n t s T ==
图2。
2 1/21/2d n t s T ==
图2.3 1d n t s T ==
图2。
4 1.25 1.25d n t s T ==
图2.5 1.5 1.5d n t s T ==
由图形可看出:当/2d n t T >时,最大位移发生在荷载离开前; 当/2d n t T <时,最大位移发生在荷载离开后; 当/2d n t T =时,最大位移发生在荷载离开时。
特别的,当d n t T =时,d t t =后没有位移。
2、其他脉冲
图2。
6 负斜率直线
图2。
7 正斜率直线图2。
8 二次抛物线
图2.9 5次抛物线
三、利用威尔逊—θ法求不同阻尼下结构自振反应
本体系度刚度为1N/m,质量为1/(4*pi^2)kg ,频率为2*pi 1s -。
故其临界阻尼
21/r c m ω==pi.分别取结构阻尼c 为:0。
05/pi,0。
1/pi ,1/pi ,1。
5/pi,进行计算。
计算结果见下图:
图3.1 0.05ξ=
ξ=
图3。
2 0.1
ξ=图3。
3 1
图3。
4 1.5ξ=
由图形对比可知,当1ξ<时,阻尼越大,结构运动衰减越快,此时结构处于小阻尼状态;当1ξ=体系不能振动,此时的c 为临界阻尼;当1ξ>时,为超阻尼系统,不发生自由振动.
四、利用威尔逊-θ法求实例
已知:如下图,W=438。
18kN;k=40181.1kN/m
求:体系位移反应.
图4。
1 结构 图4.2 荷载
由MATLAB 计算出的结果见图 4.3。
图4。
3 位移曲线
手算结果:
kN 1088.433⨯==g W m
s /119.301088.4310437=⨯⨯=
ω s
208.02==ωπ
T
T t <<=05.01
3310.02543010 5.375102
S =⨯⨯⨯=⨯ 3
335.37510 4.0571043.881030.19
S A m m ω-⨯===⨯⨯⨯ 手算结果与电算结果近似相等,说明当冲击荷载作用时间远小于结构自振周期是,可近似认为最大位移S A m ω
=(S 为荷载与坐标轴所围成的面积).。