Wilson-θ法计算步骤
哈工大结构动力学作业-威尔逊-θ法
结构动力学大作业(威尔逊- 法）姓名：学号:班级：专业：威尔逊—θ法原理及应用【摘要】在求解单自由度体系振动方程时我们用了常加速度法及线加速度法等数值分析方法。
在多自由度体系中，也有类似求解方法,即中心差分法及威尔逊—θ法。
实际上后两种方法也能求解单自由度体系振动方程。
对于数值方法，有三个重要要求：收敛性、稳定性及精度。
本文推导了威尔逊-θ法的公式，并利用MATLAB 编程来研究单自由度体系的动力特性。
【关键词】威尔逊—θ法 冲击荷载 阻尼比【正文】威尔逊-θ法可以很方便的求解任意荷载作用下单自由度体系振动问题。
实际上,当 1.37θ>时，威尔逊—θ法是无条件收敛的. 一、威尔逊—θ法的原理威尔逊-θ法是线性加速度法的一种拓展（当1θ=时,两者相同），其基本思路和实现方法是求出在时间段[],t t t θ+∆时刻的运动，其中1θ≥,然后通过内插得到i t t +∆时刻的运动（见图 1。
1)。
图 1。
11、公式推导推导由t 时刻的状态求t t θ+∆时刻的状态的递推公式:{}{}{}{})(t t t t t yy t y y -∆+=∆++θτθτ对τ积分{}{}{}{}{})(22t t t t t t yy t y y y-∆++=∆++θτθττ{}{}{}{}{}{})(6232t t t t t t t yy t y y y y -∆+++=∆++θτθτττt ∆=θτ{}{}{}{}{})(21t t t t t t t yy t y t y y -∆+∆+=∆+∆+θθθθ{}{}{}{}{})2(6)(2t t t t t tt yy t y t y y +∆+∆+=∆+∆+θθθθ {}{}{}{}{}t t t t t t t y y t y y t y26)()(62-∆--∆=∆+∆+θθθθ{}{}{}{}{}t t t t t t t yty y y t y22)(3∆---∆=∆+∆+θθθθ[]{}[]{}[]{}{}P y k y C ym =++ []{}[]{}[]{}{}t t t t t t t t P y k y C y m ∆+∆+∆+∆+=++θθθθ[]{}{}t t tt R y k ∆+∆+=θθ[][][][]c tm t k k ∆+∆+=θθ3)(62[]{}{}{}[]{}{}{}[]{}{}{})223()26)(6()(2t tt t t t t tt ty ty 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〉LpP（Lp+1:Lu)=0;endelseif isnumeric（Pt）%荷载为散点if Lu〈=LpP=Pt(1:Lu)；elseP(1:Lp)=Pt;P(Lp+1：Lu）=0；endendy=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—1kk=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）;endplot （uds ，y ，’r ’）,xlabel('时间 t ’）,ylabel('位移 y ’），title （'位移图形’） 二、利用威尔逊-θ法求冲击荷载下的结构反应1、矩形脉冲研究不同时长脉冲作用下，体系振动位移。
威尔逊—θ法在matlab中的实现
^
^
k u(ti +θ∆t) = P(ti +θ∆t)
其中
^
k=k+
6
m+ 3 c
(θ∆t)2 θ∆t
（10） （11）
•
^
P(ti
+ θ∆t )
=
Pi
+θ
( Pi +1
−
Pi
)
+
6 [ (θ∆t ) 2
ui
+
6 θ∆t
•
ui
+
••
2 ui
]m
+
(3 θ∆t
ui
+
2 ui
+
θ∆t 2
••
ui
)c
（12）
（2）
•
τ 2 ••
τ 3 ••
••
u(ti +τ ) = u(ti ) +τ u(ti ) + 2 u(ti ) + 6θ∆t [u(ti +θ∆t) − u(ti )]
（3）
当 τ = θ∆t 时，由式（2）和式（3）得到
•
u(ti
+ θ∆t )
=
•
u(ti
)
+ θ∆t
••
u(ti
)
+
θ∆t 2
加速度 0.00E+00 -7.41E-04 -1.21E-03 -1.42E-03 -1.58E-03 -1.54E-03 -1.44E-03 -1.18E-03 -8.64E-04 -4.67E-04 -2.70E-05 3.80E-04 8.26E-04 1.10E-03 2.48E-03 7.78E-03 1.49E-02 8.72E-03 -1.65E-03 -6.50E-04 -1.54E-03 -2.86E-03 -4.51E-03 -6.70E-03 -8.43E-03 -3.74E-03 1.78E-03
【国家自然科学基金】_wilson-θ法_基金支持热词逐年推荐_【万方软件创新助手】_20140801
2011年 序号 1 2 3 4 5 6 7 8 9
2011年 科研热词 立式储罐 地震动响应 土与结构相互作用 非粘滞阻尼 竖向隔震 时程积分 指数型阻尼 微分求积法 基础隔震 推荐指数 2 2 2 1 1 1 1 1 1
2012年 序号 1 2 3 4 5 6 7 8 9
科研热词 推荐指数 高层框架 1 预估校正 1 非线性动力极限承载能力 1 直接积分法 1 流固耦合 1 波浪动力响应 1 样条加权残数法 1 qr法 1 morison公式 1
2008年 序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
科研热词 推荐指数 高速列车 1 静动力分析 1 铁道工程 1 邓肯-张模型 1 逐步数值积分法 1 轨段单元 1 谱半径 1 空间振动 1 稳定性 1 状态传递算子 1 液化场地 1 液化判别 1 桩-土-桥梁结构地震相互作用 1 板式轨道 1 有限元法 1 数值模拟 1 大型振动台模型试验 1 wilson-θ 法 1
推荐指数 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1
2014年 序号 1 2 3 4 5 6 7 8 9 10
2014年 科研热词 面向对象 非线性 轴向槽径向气体轴承 算法构架 混沌 有限元 微分变换法 动力时程分析 分岔 sipesc.fems 推荐指数 1 1 1 1 1 1 1 1 1 1
推荐指数 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
2010年 序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14
科研热词 推荐指数 非传统hamilton型变分原理 3 辛算法 3 相空间 3 动力响应 3 初值-边值问题 3 非线性响应 1 附加运动 1 陀螺效应 1 约束条件 1 稳定性 1 平衡方程 1 分岔 1 修正计算 1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
子结构拟动力试验的CD-Wilson θ法研究
Ｒｅ ｓ ｅ ａ ｒ ｃ ｈ ｏ ｎ ＣＤ － Ｗ ｉ ｌ ｓ ｏ ｎ ０ ｍｅ ｔ ｈ ｏ ｄ ｉ ｎ ｓ ｕ ｂ ｓ ｔ ｒ ｕ ｃ ｔ ｕｒ ｅ ｐ ｓ ｅ ｕ ｄｏ ｄ ｙ ｎ ａ ｍｉ ｃ ｔ ｅ ｓ ｔ ＷＡ Ｎ Ｇ Ｘ ｉ ａ ｏ ｆ ｅ ｎ ｇ ， Ｃ Ａ Ｉ Ｘ ｉ ｎ ｊ ｉ ａ ｎ ｇ ， Ｔ Ｉ Ａ Ｎ Ｓ ｈ ｉ ｚ ｈ ｕ ’
Ａｂ ｓ ｔ ｒ ａ ｃ ｔ ： Ｔ ｈ ｅ ａ ｄ ｖ ａ ｎ ｔ ａ ｇ ｅ ｓ ｏ ｆ ｃ ｏ ｍｂ ｉ ｎ ａ ｔ ｉ ｏ ｎ ｏ ｆ ｎ ｕ ｍｅ ｉ ｒ ｃ ａ ｌ ｉ ｎ ｔ ｅ ｇ ｒ ａ ｔ ｉ ｏ ｎ ｍｅ ｔ ｈ ｏ ｄ ｃ ａ ｌ ｌ ｂ ｅ ｃ ｏ ｎ ｓ ｉ ｄ ｅ ｒ ｅ ｄ ａ ｓ ｔ ｈ ｅ ｅ ｘ ｐ ｌ ｉ ｃ ｉ ｔ ｎ ｕ ｍｅ ｉｃ ｒ ｌ ａ ｉ ｎ ｔ ｅ ｒａ ｇ ｔ ｉ ｏ ｎ ｍｅ ｔ ｈ ｄ ｏ ｗ ｉ ｈｏ ｔ ｕ ｔ ｉ ｔ ｅ ｒ ａ ｔ ｉ ｏ ｎ ａ ｎ ｄ ｔ ｈ ｅ ｉ ｍｐ ｉ ｌ ｃ ｉ ｔ ｎ ｕ ｍｅ ｒ ｉ ｃ ｌ ａ ｉ ｎ ｔ ｅ ｒａ ｇ ｔ ｉ ｏ ｎ ｍｅ ｔ ｈ ｄ ｏ ｗ ｉ ｔ ｈ ｕ ｎ ｃ ｏ ｎ ｄ ｉ ｔ ｉ ｏ ｎ ａ ｌ ｓ ｔ ａ ｂ ｉ ｉ ｌ ｔ ｙ ｉ ｎ ｔ ｈ ｅ ｓ ｕ ｂ ｓ ｔ ｒ ｕ ｃ ｔ ｕ ｒ ｅ ｐ ｓ ｅ ｕ ｄ ｄｙ ｏ ｎ ａ ｍｉ ｃ ｔ ｅ ｓ ｔ ． Ａｃ ｃ ｏ ｒ ｄ ｉ ｎ ｇ ｔ ｏ ｈｅ ｔ Ｃ Ｄ・ Ｎｅ ｗ ｍａ ｒｋ ｍｅ ｔ ｈ ｄ ｏ ｉｔ ｗ ｈ ｌ ｒｇ ａ ｅ ｕ ｐ ｐ ｅ ｒ ｆ ｅｑ ｒ ｕ ｅ ｎ ｃ ｅ ｔ ｈ ａ ｔ ｌ ｅ ａ ｄ ｓ ｔ ｏ ｔ ｈ ｅ ｔ ｉ ｍｅ ｓ ｔ ｅ ｐ ｂ ｅ ｃ ｏ ｍｅ ｓ ｍａ ｌ ｌ ｅ ｒ ， ｔ ｈ ｅ ＣＤ－ Ｗｉ ｌ ｓ ｏ ｎ ０ ｍｅ ｔ ｈ ｄ ｏ ｉ ｓ ｃ ｏ ｍｂ ｉ ｎ ｅ ｄ ｂ ｙ ｅ ｘ ｐ ｌ ｉ ｃ ｉ ｔ ｃ ｅ ｎ ｔ ｅ ｒ ｄ ｉ ｆ ｅ ｒ ｅ ｎ ｃ ｅ ｍｅ ｔ ｈ ｄ ｏ ｎｄ ａ ｉ ｍｐ ｌ ｉ ｃ ｉ ｔ Ｗｉ ｌ ｓ ｏ ｎ － ０ ｍｅ ｔ ｈ ｄ． ｏ Ａｆ ｔ ｅ ｒ ｃ ｌｃ ａ ｕ ｌ ａ ｔ ｉ ｏ ｎ ａ ｎ ｄ ｍｏ ｄｉ ｉ ｆ ｃ ａ ｔ ｉ ｏ ｎ ｔ ｈ ｅ Ｗｉ ｌ ｓ ｏ ｎ － ０ ｍｅ ｔ ｈ ｄ， ｏ
4.2 时域分析法--直接积分法
3. 对每一时间步计算等效荷载增量
4. 解方程计算位移 t τ
u
K Ft τ
-1
5.计算时刻t+Δ t的响应
用增量形式表示的Wilsonθ法 （推导省略）
增量方程
• 计算公式
K Δ ut τ Δ Ft τ
t τ
• 得到
时刻的响应，
• 再转变成t+Δ t的响应：
直接积分法的补充说明
显式积分、隐式积分
• 显式积分是在第i步计算中状态ti满足运动方程式的
计算方法。
ui1 ui Δ t f ti , ui
•当前时刻的位移只与前一时刻的加速度和位移有关，这就意味着当前 时刻的位移求解无需迭代过程。 另外，只要将运动方程中的质量矩阵 和阻尼矩阵对角化（线性无关），前一时刻的加速度求解无需解联立 方程组，从而使问题大大简化，这就是所谓的显式求解法。
x
x (t )
f s
结构在t时刻的刚度矩阵 由t时刻结构各构件的切线刚度确定
x(t t )
x (t )
(t ) [ K (t )]x(t ) P(t ) [M ] x(t ) [C(t )]x
----增量方程(3)
方程左边的力增量表达式是近似的！
常用的隐式积分法
例：
g (t )kN P(t ) mx
2.5
求位移时程曲线， 恢复力时程曲线， 最大位移，最大恢复力， 开始时静止。 x
g
W=15kN
x(t)
计算步骤： 1.确定积分步长t 2.确定当前积分步长内结构的 质量，刚度和阻尼矩阵以及 阻尼力和恢复力 3.计算初始加速度 4.确定等效刚度K*和等效荷载 矩阵P*
wilson法和newmark法的理论过程
第三章离散化结构动力方程的解法(2013.4.24）§3.1 绪言对于一个实际结构，由有限元法离散化处理后，应用瞬时最小势能原理可导出动力方程[]{}[]{}[]{}{}++=（3.1）M u C u K u F(t)这里，{}u、{}u、{}u及{}F t分别表示加速度、速度、位移及所()作用的外力矢量，他们都是与时间有关的。
从数学的角度来看，式（3.1）是一个常系数的二阶线性常微分方程组，对于它的求解原则上并无困难。
但是，由于[]M、[]C 和[]K的阶数非常高，使得式（3.1）的求解必须花费很大的代价，便促使人们去寻求一些效率高的近似计算方法。
目前，用于求解式（3.1）的方法，大致可分为两大类。
一是坐标变换法，它是对结构动力方程式（3.1），在求解之前，进行模态坐标变换，实际上就是一种Ritz变换，即把原物理空间的动力方程变换到模态空间中去求解。
现在，普遍使用的方法是模态（振型）迭加法，即用结构的前q阶实际主模态集（主振型阵）构成坐标变换阵进行变换。
通过这一变换，实现降阶，求较好的近似解，而且，还用解除耦合的办法，简化方程的计算。
还有一种所谓假设模态法，即是用一组假设模态，构成模态坐标变换阵进行变换，获得一组降阶的而不解耦的模态基坐标方程。
显然，这种方法的计算精度，取决于所假设的模态。
用Ritz矢量法求解的近似模态作为假设模态，可得到满足要求的精度。
二是直接积分法，它是对式（3.1）在求解之前，不进行坐标变换，直接进行数值积分计算。
这种方法的特点是对时域进行离散，将式（3.1）分为各离散时刻的方程，然后，将该时刻的加速度和速度用相邻时刻的各位移线性组合而成，于是，式（3.1）就化为一个由位移组成的该离散时刻上的响应值，通常又称为逐步积分法。
线性代数方程组的解法与静力时刻的位移来线性组合，就导致了各种不同的方法。
主要有中央差分法，Houbolt 方法，Wilson -θ法和Newmark 方法等。
wilson法和newmark法的理论过程要点
第三章离散化结构动力方程的解法(2013.4.24）§3.1 绪言对于一个实际结构，由有限元法离散化处理后，应用瞬时最小势能原理可导出动力方程[]{}[]{}[]{}{}（3.1）++=M u C u K u F(t)这里，{}u 、{}u 、{}u及{}F t分别表示加速度、速度、位移及所()作用的外力矢量，他们都是与时间有关的。
从数学的角度来看，式（3.1）是一个常系数的二阶线性常微分方程组，对于它的求解原则上并无困难。
但是，由于[]M、[]C 和[]K的阶数非常高，使得式（3.1）的求解必须花费很大的代价，便促使人们去寻求一些效率高的近似计算方法。
目前，用于求解式（3.1）的方法，大致可分为两大类。
一是坐标变换法，它是对结构动力方程式（3.1），在求解之前，进行模态坐标变换，实际上就是一种Ritz变换，即把原物理空间的动力方程变换到模态空间中去求解。
现在，普遍使用的方法是模态（振型）迭加法，即用结构的前q阶实际主模态集（主振型阵）构成坐标变换阵进行变换。
通过这一变换，实现降阶，求较好的近似解，而且，还用解除耦合的办法，简化方程的计算。
还有一种所谓假设模态法，即是用一组假设模态，构成模态坐标变换阵进行变换，获得一组降阶的而不解耦的模态基坐标方程。
显然，这种方法的计算精度，取决于所假设的模态。
用Ritz矢量法求解的近似模态作为假设模态，可得到满足要求的精度。
二是直接积分法，它是对式（3.1）在求解之前，不进行坐标变换，直接进行数值积分计算。
这种方法的特点是对时域进行离散，将式（3.1）分为各离散时刻的方程，然后，将该时刻的加速度和速度用相邻时刻的各位移线性组合而成，于是，式（3.1）就化为一个由位移组成的该离散时刻上的响应值，通常又称为逐步积分法。
线性代数方程组的解法与静力时刻的位移来线性组合，就导致了各种不同的方法。
主要有中央差分法，Houbolt 方法，Wilson -θ法和Newmark 方法等。
结构动力学-第五章 数值分析方法 (Part 2)
§5.5 Wilson-θ 法
ti +1 时刻的解
}i +1 = {u }i +1 {u 6 θ 3 Δt 2
({u}
i +θ
{u}i +1
Δt }i + ({u }i +1 + {u }i ) = {u 2 Δt 2 }i + }i +1 + 2 {u }i ) = {u}i + Δt {u {u ( 6
结构动力学
第五章 动力反应数值分析方法
11 of 23
华南理工大学
土木与交通学院
土木工程系
§5.5 Wilson-θ 法
不同数值积分法计算精度的比较
(0) = 0 考虑无阻尼自由振动问题： mu + ku = 0 u (0) = 1, u
步长： Δt = 0.1 × Tn
结构动力学 第五章 动力反应数值分析方法 12 of 23 华南理工大学 土木与交通学院 土木工程系
i +θ i
i +1
− {P}i ) +
⎡ 6 ⎤ 6 }i + 2 {u }i ⎥ + u + {u [M ] ⎢ 2 { }i θ Δt ⎢ ⎥ ⎣ (θ Δt ) ⎦ ⎛ 3 ⎞ θ Δt }i + }i ⎟ {u [C ] ⎜ {u}i + 2 {u 2 ⎝ θ Δt ⎠
结构动力学 第五章 动力反应数值分析方法 9 of 23 华南理工大学 土木与交通学院 土木工程系
4 of 23
华南理工大学
土木与交通学院
土木工程系
§5.5 Wilson-θ 法
加速度变化规律
( ti ) + ατ a ( ti + τ ) = u (0 ≤ τ ≤ θ Δt )