班级： 学号：
姓名： 目录。
1. wilson-θ法原理简介 2
2. wilson-θ程序验算 3
2.1 △t的影响 4
2.2 θ的影响 5
3. 非线性问题求解 5
4. 附录 8
wilson-θ法源程序 8
wilson-θ法是基于对加速度a的插值近似得到的，图1-1为wilson-θ法的原理示意图。
推导由t时刻的状态求t+△t时刻的状态的递推公式：
对τ积分可得速度与位移的表达式如下：
其中τ=θt，由式(1-2)、(1-3)可以解出：
将式(1-4)、(1-5)带入运动方程：
注意到此时的式子为{}和上一个时刻、、以及t+θ△t时刻的荷载相关，可以运用迭代的思想来求解，下图给出线弹性条件下wilson-θ法的流程图：
对线弹性条件下的wilson-θ法进行matlab编程，源**见附录。选取如下算例进行验证。
对于一个单自由度的无阻尼结构，当其受到一个周期荷载时，其结构响应分为稳态解和瞬态解，由于没有阻尼的影响，其瞬态解并不会衰减，其理论表达式为：
式中，为位移响应，为激励，为刚度，为荷载频率与固有振动频率之比，为荷载频率，为结构固有频率。
现令为1，为1，则为1，取为2/3。程序求得的解与解析解对比如图2-1所示（由于理论解与程序基本重合，所以将理论解乘以-1，方便比较）：
上述算例验证时选择的△t非常小，因此看不出理论解与wilson-θ法的求解区别，以下改变△t的取值，**△t对迭代的影响。
可以看出并不是△t太大时计算结果很不准确，偏小，反映不出周期特征；当△t合适时正好基本和理论解重合，也不是△t越小越好越小时越能反映出一些细部特征，但这也不是很准确。
当θ>1.37时，该算法是无条件稳定的算法，以下**θ对算法的影响。
由上图可知随着θ值越大，位移的周期变大。
由于实际结构并不一定为线性，其刚度会随着位移的的变化而改变，下图为求解非线性问题时的wilson-θ法流程。此处要说明的是，刚度矩阵[k]y(t)是与位移相关的量，判断那时候速度的大小是为了确定其是否处于卸载段。具体可能得根据实际情况求解。
修改matlab程序，并用该程序来计算如下例题：
对该问题采用wilson-θ法非线性方式计算，采用△t=0.1s和△t=0.05s两种方式，计算位移、速度和加速度曲线如下图所示：
由上图可知，结构在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。
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代表加速度。
if nargin<9
theta=1.4;
endl,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;
计算拟荷载。
for i=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(:,i));
end上述**只适合分析线弹性结构，对于非线性结构，编写起来比较繁琐，针对不同的情况可能需要具体处理，所以本文只给出了针对本文例题的**。与上述**不同的是以下**增加了一个判断的语句。
function [y_1,y_2,y_3]=wilson_theta2(p,m,c,dt,v0,y0,a_0,theta)
p代表输入的荷载，c为阻尼矩阵，dt为时间间隔，m为质量矩阵。
v0为初始的速度，y0为初始的位移。a_0为初始加速度。
输出的矩阵y_1代表位移，y_2代表速度，y_3代表加速度。
if nargin<9
theta=1.4;
endl,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;
for i=1:r-1
if y_2(:,i)>0
k=60*(y_1(:,i)<=0.05)+3/y_1(:,i)*(y_1(:,i)>0.05);
elsek=60;
end %计算拟刚度矩阵。
k0=k+a0*m+a1*c;
%计算拟荷载。
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(:,i));end
结构动力学大作业参考
题目信息 1 试设计一个3层框架，给出框架结构的一致质量矩阵 一致刚度矩阵，建立框架结构的运动微分方程，求出该框架结构的各阶频率和振型 并采用振型分解法，选定一个正弦动力荷载，求3层框架对于该荷载的位移反应 2 试设计一个3层框架，采用时程分析法，输入 波，求出所设计框架各层的非线性位移时程反应，要...
动力学作业
1 如图4 1所示，质量m1 40公斤的人用轻绳通过轻滑轮拉一个质量m2 20公斤的物体，在拉至图示位置 300时，因地面打滑无法前进，此时人对地面的压力f 人与地面间的静摩擦系数 2 质量m 3公斤的物体，t 0时，xo 0，vo 0，受到沿x轴的外力牛顿的作用，则t 3秒时质点的速度 质点的位置...
结构动力学考试提纲
结构动力学考试提纲 闭卷考试 说明 1.闭卷考试。2.考试时间 1月13日周一晚19 00 21 003.总成绩 80 的卷面成绩 20 的平时成绩 作业，出勤 4.考试内容与提纲完全一致。一，简答部分 50 1.动力学与静力学概念相关。2.动力自由度相关。3.rayleigh阻尼模型的确定。4.不...