结构动力学大作业Word文档格式.docx

结构动力学大作业Word文档格式.docx

《结构动力学大作业Word文档格式.docx》由会员分享，可在线阅读，更多相关《结构动力学大作业Word文档格式.docx（9页珍藏版）》请在冰豆网上搜索。

其中τ=θt，由式（1-2）、（1-3）可以解出：

（1-4）

（1-5）

将式（1-4）、（1-5）带入运动方程：

（1-6）

（1-7）

注意到此时的式子为{

}和上一个时刻

、

以及t+θ△t时刻的荷载

相关，可以运用迭代的思想来求解，下图给出线弹性条件下Wilson-θ法的流程图：

图1-2Wilson-θ法流程图

2.Wilson-θ程序验算

对线弹性条件下的Wilson-θ法进行MATLAB编程，源代码见附录。

选取如下算例进行验证。

对于一个单自由度的无阻尼结构，当其受到一个周期荷载时，其结构响应分为稳态解和瞬态解，由于没有阻尼的影响，其瞬态解并不会衰减，其理论表达式为：

（2-1）

式中，

为位移响应，

为激励，

为刚度，

为荷载频率与固有振动频率之比，

为荷载频率，

为结构固有频率。

现令

为1，

为1，则

取为2/3。

程序求得的解与解析解对比如图2-1所示（由于理论解与程序基本重合，所以将理论解乘以-1，方便比较）：

a）位移

b）速度

c）加速度

图2-1Wilson-θ法结果验证

2.1△t的影响

上述算例验证时选择的△t非常小，因此看不出理论解与Wilson-θ法的求解区别，以下改变△t的取值，探讨△t对迭代的影响。

图2-2△t对位移曲线的影响

可以看出并不是△t太大时计算结果很不准确，偏小，反映不出周期特征；

当△t合适时正好基本和理论解重合，也不是△t越小越好越小时越能反映出一些细部特征，但这也不是很准确。

2.2θ的影响

当θ>

1.37时，该算法是无条件稳定的算法，以下探讨θ对算法的影响。

图2-3θ对位移曲线的影响

由上图可知随着θ值越大，位移的周期变大。

3.非线性问题求解

由于实际结构并不一定为线性，其刚度会随着位移的的变化而改变，下图为求解非线性问题时的Wilson-θ法流程。

此处要说明的是，刚度矩阵[K]y（t）是与位移相关的量，判断那时候速度的大小是为了确定其是否处于卸载段。

具体可能得根据实际情况求解。

图3-1Wilson-θ法解非线性问题

修改MATLAB程序，并用该程序来计算如下例题：

对该问题采用Wilson-θ法非线性方式计算，采用△t=0.1s和△t=0.05s两种方式，计算位移、速度和加速度曲线如下图所示：

a）位移

图3-2非线性分析结果

由上图可知，结构在0.6s时达到位移极值，在△t=0.1s和△t=0.05s算得的值分别为0.096m和0.108m，速度极值在0.9s取到分别为-0.468和-0.580，加速度极值在△t=0.1s时为0.7s时取到，为-2.127，在△t=0.05s极值在0.75s时取到，为-2.9。

4.附录

Wilson-θ法源程序

function[y_1,y_2,y_3]=wilson_theta（p,m,c,k,dt,v0,y0,a_0,theta）

%p代表输入的荷载，c为阻尼矩阵，dt为时间间隔,m为质量矩阵,k为刚度矩阵

%v0为初始的速度，y0为初始的位移.a_0为初始加速度

%输出的矩阵y_1代表位移,y_2代表速度，y_3代表加速度

ifnargin<

9

theta=1.4;

end

[L,r]=size（p）;

y_1=NaN（L,r）;

y_2=NaN（L,r）;

y_3=NaN（L,r）;

y_1（:

1）=y0;

y_2（:

1）=v0;

y_3（:

1）=a_0;

%计算积分常数

a0=6/（（theta*dt）^2）;

a1=3/theta/dt;

a2=2*a1;

a3=theta*dt/2;

a4=a0/theta;

a5=-a2/theta;

a6=1-3/theta;

a7=dt/2;

a8=dt^2/6;

%计算拟刚度矩阵

k0=k+a0*m+a1*c;

%计算拟荷载

fori=1:

r-1

R=p（:

i）+theta*（p（:

i+1）-p（:

i））+m*（a0*y_1（:

i）+a2*y_2（:

i）+2*y_3（:

i））+c*（a1*y_1（:

i）+2*y_2（:

i）+a3*y_3（:

i））;

y_theta=k0\R;

y_3（:

i+1）=a4*（y_theta-y_1（:

i））+a5*y_2（:

i）+a6*y_3（:

i）;

y_2（:

i+1）=y_2（:

i）+a7*（y_3（:

i）+y_3（:

i+1））;

y_1（:

i+1）=y_1（:

i）+y_2（:

i）*dt+a8*（y_3（:

i+1）+2*y_3（:

上述代码只适合分析线弹性结构，对于非线性结构，编写起来比较繁琐，针对不同的情况可能需要具体处理，所以本文只给出了针对本文例题的代码。

与上述代码不同的是以下代码增加了一个判断的语句。

function[y_1,y_2,y_3]=wilson_theta2（p,m,c,dt,v0,y0,a_0,theta）

%p代表输入的荷载，c为阻尼矩阵，dt为时间间隔,m为质量矩阵

ify_2（:

i）>

k=60*（y_1（:

i）<

=0.05）+3/y_1（:

i）*（y_1（:

0.05）;

else

k=60;

end

%计算拟刚度矩阵

k0=k+a0*m+a1*c;

%计算拟荷载