结构动力学多自由度线性体系Wilson-θ法程序编写
结构动力学大作业2
结构动力学大作业班级：学号：姓名：目录1. Wilson-θ法原理简介 (2)2. Wilson-θ程序验算 (3)2.1△t的影响 (4)2.2 θ的影响 (5)3. 非线性问题求解 (5)4. 附录 (8)Wilson-θ法源程序 (8)1. Wilson -θ法原理简介图1-1Wilson-θ法示意图Wilson-θ法是基于对加速度a 的插值近似得到的，图1-1为Wilson-θ法的原理示意图。
推导由t 时刻的状态求t +△t 时刻的状态的递推公式：{}{}{}{}()t tt t t y y y y tτθτθ++∆=+-∆ (1-1)对τ积分可得速度与位移的表达式如下：{}{}{}{}{}2()2t t t t t t yy y y ytτθττθ++∆=++-∆ (1-2){}{}{}{}{}{}23()26t t t t t t t y y y y y ytτθτττθ++∆=+++-∆ (1-3)其中τ=θt ，由式(1-2)、(1-3)可以解出：{}{}{}{}{}266()2()t t t tt t t y y y y y t tθθθθ+∆+∆=---∆∆(1-4){}{}{}{}{}3()22t t t t t t t tyy y y y t θθθθ+∆+∆∆=---∆(1-5)将式(1-4)、(1-5)带入运动方程：[]{}[]{}[]{}{}m y C y k y P ++=(1-6)[]{}[]{}[]{}{}t t t t t t t tm y C y k y P θθθθ+∆+∆+∆+∆++= (1-7)注意到此时的式子为{{}t t y θ+∆}和上一个时刻{}t y 、{}t y、{}t y 以及t +θ△t 时刻的荷载{}t t P θ+∆相关，可以运用迭代的思想来求解，下图给出线弹性条件下Wilson -θ法的流程图：图1-2Wilson-θ法流程图2.Wilson-θ程序验算对线弹性条件下的Wilson-θ法进行MATLAB编程，源代码见附录。
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 方法等。
哈工大结构动力学作业-威尔逊-θ法
结构动力学大作业(威尔逊- 法）姓名：学号:班级：专业：威尔逊—θ法原理及应用【摘要】在求解单自由度体系振动方程时我们用了常加速度法及线加速度法等数值分析方法。
在多自由度体系中，也有类似求解方法,即中心差分法及威尔逊—θ法。
实际上后两种方法也能求解单自由度体系振动方程。
对于数值方法，有三个重要要求：收敛性、稳定性及精度。
本文推导了威尔逊-θ法的公式，并利用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
中心差分法、纽马克法和威尔逊-θ法与精确解的误差分析
中心差分法、纽马克法和威尔逊-θ法与精确解的误差分析作者：于津津贾慧敏宋敏来源：《教育周报·教育论坛》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。
第十四章结构动力学多自由度
第 j 阶振型的广义质量 第 j 阶振型的广义刚度
j
K
* j
M
* j
32 / 67
§14-6 多自由度结构的自由振动
5、多自由度体系振型正交性应用举例
例题1 检验框架结构振型的正确性
m2
A1
1 1.618
A2
1 0.618
h2
m1
M
m
0
0 m
K
12EI l3
2 1
1
1
h1
A1TM A2 0
1 2 n
K 2M A 0
M
1
2
I
A
0
A1,A2, ,An
25 / 67
§14-6 多自由度结构的自由振动
注意 自振频率与主振型一一对应 振型只表明振动的形状，不能唯一确定其幅值 振型是多自由度特有的概念
3、振型的标准化 补充条件，使主振型用确定的幅值表示 方法1：规定振型中某元素为1，其它元素就是 相对于它的比值。通常选第一个元素或最大一 个元素，令其等于1。
k21 L kn1
k22 L kn2
L L L
k1n y1 FP1
k2n L knn
yM2 yn
FMP2 FPn
简写为 M Y&& Y& K Y FP
其中[M]为质量矩阵，在集中质点的结构中它是对角矩阵；[β] 为阻尼矩阵；[K]为刚度矩阵，根据反力互等定理，它是对称矩 阵。
一、刚度法
取任意质量mi列动力平衡方程：
FIi FRi Fei FPi 0 n
根据叠加原理有
Fei kij y j
jn1
FRi ij yj
j 1
结构动力学多自由度
▪ 振型方程：
(K i2M)ji 0 (i 1, 2, 3, n)
▪∵
K 2i M 0
▪ ∴ 第i 个振型方程中的n 个方程中只有n-1个是独立的！ ▪ ——无法得到j1i、 j2i、 … 、 jni 的确定值， ▪ 但可以确定各质点振幅之间的相对比值： ▪ —— 振型的幅值是任意的，但形状是惟一的。
一致质量矩阵：
L
pava m13v1 0 fI ( x)v( x)dx
L
0
m( x) 3( x)v3
L
1( x)v1dx
mij 0 EI ( x)i ( x) j ( x)dx
L
cij 0 c( x) i ( x) j ( x)dx
其中，c(x)表示分布的粘滞阻尼特性。
一致节点荷载
L
vˆ 表示体系的形状，不随时间变化。
v 2vˆ sin(t ) 2v 2mvˆ sin(t ) kvˆ sin(t ) 0
k 2m vˆ 0
k 2m vˆ 0
即： k 2m 0
上式的N个根，表述体系可能存在的N个振型的频率。
1
2
3
N
2)
2
)
y32
(t
)
jˆ
32
s
in
(
2
t
2
)
1
jˆ
2i
yi
(t
)
jˆ3i
s
in(i
t
i
)
jˆ ni
1
jˆ 21
jˆ 31
jˆ 32
1
jˆ 22
将N个振型中的每一振型形式，用F表示N个振型所组成的方阵。
11 12 13 1N
结构动力响应数值计算方法对比分析
结构动力响应数值计算方法对比分析作者：李涵来源：《青年生活》2019年第21期摘要：中心差分法、纽马克法、威尔逊-法是结构动力学中常用的三种方法，为了系统的比较其优缺性，本文针对一个双自由度的体系，首先根据已知条件计算出振动微分方程，运用Matlab计算出可求出12个步长内相应的位移值，即精确解。
然后分别运用中心差分法，纽马克法，威尔逊-法求出其近似解;最后通过三种方法的近似解与精确解相对比，进而分析出三种计算方法的优缺性，为结构动力计算提供依据。
关键词：动力计算、中心差分法、纽马克法、威尔逊-法1、动力体系概况2、精确解推导针对该双自由度体系，理论推导出系统的位移表达式，通过代入各时刻周期得出位移在各时刻的具体数值，即位移精确解。
对位移方程求一阶导数得出速度方程，求二阶导数求出加速度方程。
代入各时刻的周期值，通过Matlab计算得出位移、速度、加速度的数值如下：3、三种数值计算方法3.1、中心差分法中心差分法是基于用有限差分代替位移对时间的求导，对位移一阶求导得到速度，对位移二阶求导得加速度。
通过Matlab计算出前4个步长所对应的位移响应。
3.2、纽马克法纽马克-β法是一种将线性加速度方法普遍化的方法。
通过Matlab计算出前4个步长所对应的位移响应。
3.3威尔逊-法通过Matlab计算出前4个步长所对应的位移响应。
4、近似解与精确解对比分析从上述结构的位移、速度、加速度可以看出，三种方法都能大致表示该体系大体运动趋势，并且误差较小。
其中，在描述物体位移时，中心差分法较后两种方法更为精确。
然而在描述速度和加速度时，中心差分法表现出了较大的误差，而纽马克和威尔逊法则能更详尽的表征物体速度和加速度。
5、结论中心差分法、纽马克法和威尔逊-θ法均是结构动力计算中的常用方法。
本文针对具体的计算实例，分别计算出三种方法的动力响应结果，并与精确解进行对比。
经过分析，中心差分法能更精确的表示物体位移响应，而纽马克和威尔逊法在表征物体速度和加速度方面相较于中心差分法更为精确，三种方法，各有其优缺点，应视具体情况采用相应的计算方法。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
多自由度线性体系Wilson -θ法程序编写
【摘要】本文主要介绍了通过使用Matlab 软件，Wilson-θ法编写多自由度线性
体系的程序的原理、流程图、具体算例以及使用注意事项。
通过该程序可以得到剪切型结构在任意函数荷载作用下各质点的位移函数。
【关键词】Matlab ；多自由度；Wilson-θ法
1.wilson-θ法原理
wilson-θ法中最主要的步骤就是推导由t 时刻的状态求t t ∆+时刻的状态的递推公
式，现推导如下：
对τ积分
解出
代入
整理，得
其中
本程序的核心就是对以上公式的循环使用。
{}{}{}{})(t t t t t y y t
y
y -∆+=∆++θτθτ
t ∆=θτ{}{}{}{}{})(22
t t t t t t y
y t y y y
-∆++=∆++θτθττ{}{}{}{}{}{})(623
2t t t t t t t y
y t
y y y y -∆+++=∆++θτ
θτττ{}{}{}{}{})(21
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 2
2)(3∆---∆=∆+∆+θθθθ[]{}[]{}[]{}{}t t t t t t t t P y k y C y
m ∆+∆+∆+∆+=++θθθθ []{}[]{}[]{}{}P y k y C y
m =++ []{}[]
R y k t t =∆+θ[]
[][][]
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 ∆++∆++∆+∆+-+=∆+θθθθθ{}{}{}{})
(t t t t t t P P P P -+=∆+∆+θθ
2.程序流程图
求出各常数值
For I=1 to n
[][][][]
c a m a k k 1
++=
3.具体应用算例
如图所示，两自由度框架结构，其中
初始静止，求各层位移。
将相应的数据输入到程序中，得出各层位移关于时间的图像。
图1为第一层，图2为第二层。
将所得数值解与精确解相比较，图中实线为数值解，虚线为精确解。
由两张图，我们可以看出数值解大致是与精确相近的，但是仍然有些许的不同，这可能是算法中仍然有缺陷，说明程序仍然有待改善。
kg
1020021==m m kN/m 30001=k kN/m 20002=k kN 217.0=P s
/145.11=
θ图1
4.程序使用注意事项
（1）本程序针对于剪切型刚架结构，对于其他结构无法使用。
（2）本程序中各质点的荷载必须是函数的形式（包括常数），即对于只有某些点的荷载无法使用，且荷载函数输入时，必须采用inline语句。
例如荷载为常数10，则输入inline（’10’）；如荷载函数为sin（at），则输入inline（’sin(a*t)’,’t’）。
（3）本程序主要针对无阻尼情况，若有阻尼，只需输入阻尼矩阵即可。
（4）θ的值应大于1.37，通常取1.4，优化值为1.420815。
（5）从第一层开始为m
1，m
2
……
【参考文献】
[1] 王焕定. 结构力学（第3版）[M]. 北京：高等教育出版社，2010.
[2] Anil K.Chopra. 结构动力学理论及其在地震工程中的应用（第2版）[M]. 北京：高等教育出版社，2007.
图2。