哈工大结构动力学作业-威尔逊-θ法

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('平面连杆机构运动角度——角加速度图');这样，在运行程序时就会弹出一个如下图所示的对话框，可以任意给定初值，解决不同问题。

中心差分法、纽马克法和威尔逊-θ法与精确解的误差分析作者：于津津贾慧敏宋敏来源：《教育周报·教育论坛》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。

结构动力学大作业对于如下结构，是研究质量块的质量变化和在简支梁上位置的变化对整个系统模态的影响。

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 +++=阵。

第五章 结构动力学中的常用数值方法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 -θ法程序编写【摘要】本文主要介绍了通过使用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.具体应用算例如图所示，两自由度框架结构，其中初始静止，求各层位移。

第三章静定结构的位移计算Displacement of Statically Determinate Structures１．　弹性杆件的变形与变形能计算２．　变形体虚功原理３．　单位荷载法４．　图乘法５．　其他外因引起的位移计算６．　互等定理７．　结论与讨论1 结构位移计算概述一、结构的位移 (Displacement of Structures)x Δy ΔA A ′βF P 线位移，角位移，相对线位移、角位移等统称广义位移线位移角位移DC ΔΔ+相对线位移CD D ΔC ΔαF P γ相对角位移γα +制造误差 δ 等铁路工程技术规范规定: 二、 计算位移的目的引起结构位移的原因(1) 刚度要求如：荷载、温度改变 ΔT 、支座移动 Δc 、在工程上，吊车梁允许的挠度< 1/600 跨度；桥梁在竖向活载下，钢板桥梁和钢桁梁最大挠度 < 1/700 和1/900跨度高层建筑的最大位移< 1/1000 高度。

最大层间位移< 1/800 层高。

（3）理想联结 (Ideal Constraint)。

三、 本章位移计算的假定(2) 超静定、动力和稳定计算(3）施工要求叠加原理适用（principle of superposition)(1) 线弹性 (Linear Elastic),(2) 小变形 (Small Deformation),返首2 变形体虚功原理 (Principle of Virtual Work)一、功(Work)、实功(Real Work)和虚功(Virtual Work)两种状态力状态位移状态F P F P /2F P /2（虚）力状态（虚力状态）（虚位移状态）无关（虚）位移状态q注意：（1）属同一体系；（2）均为可能状态。

即位移应满足变形协调条件；力状态应满足平衡条件。

（3）位移状态与力状态完全无关；一些基本概念：实功：广义力在自身所产生的位移上所作的功功：力×力方向位移之总和广义力：功的表达式中，与广义位移对应的项功：广义力×广义位移之总和虚功：广义力与广义位移无关时所作的功W=F P ×Δ/2W=F P1×Δ11 /2W=F P2×Δ22 /2W=F P1×Δ12or W=F P2×Δ21变力功（1）质点系的虚位移原理具有理想约束的质点系，在某一位置处于平衡的必要和充分条件是：1P F 2N F 1N F 2P F 1m 2m 二、变形杆件的虚功原理Σf i δr i =0→→.对于任何可能的虚位移，作用于质点系的主动力所做虚功之和为零。

基于Wilson-θ和Newmark-β法的非线性动力学方程改进算法刘广;刘济科;陈衍茂【摘要】When using Wilson-θ or Newmark-β method to solve nonlinear dynamic equations,usually,we rewrite the equations in the form of incremental equilibrium equations.The coefficient matrix has to be updated in each integral step according to the state variables.In essence,this procedure is to linearize the considered nonlinear system in a single time step.It is usually difficult,however,to handle some strongly nonlinear problems with multiple-degrees-of-freedom.By using an incremental process,in the paper a new fast algorithm was proposed based on Wilson-θ method or Newmark-β method.According to the obtained solution at one time point,we present an initial guessed solution at the next time instant.Then the guessed solution can converge to the true solution via iterative corrections.Numerical examples show that,using the presented fast algorithm together with Wilson-θ method or Newmark-β method,one can get highly accurate solutions.Moreover,the presented algorithm can provide us with a simple way to adjust the convergence as necessary.As the presented methods avoid linearization of the considered nonlinear dynamic systems,they are not only more applicable but also more computationally efficient.%Wilson-θ法和Newmark-β法是非线性动力学方程求解的常用方法.它们的一个基本步骤是,将方程改写为增量平衡的形式,在每一个积分步长内用状态参量修正平衡方程的系数矩阵,其本质是在单个步长内对系统的非线性环节进行了线性化处理.本文基于增量思想分别改进了Wilson-θ法和Newmark-β法,根据即时解给出下一步的猜测解,然后对猜测解进行迭代校正,最终得到收敛的近似解.算例表明,改进算法的精度更高,且收敛准则简单.更为重要的是,本文方法无须对非线性项进行线性化处理,因而计算效率更高,适应范围更广.【期刊名称】《计算力学学报》【年(卷),期】2017(034)004【总页数】7页(P433-439)【关键词】Wilson-θ法;Newmark-β法;非线性动力学方程;猜测解;迭代校正【作者】刘广;刘济科;陈衍茂【作者单位】中山大学力学系,广州 510275;中山大学力学系,广州 510275;中山大学力学系,广州 510275【正文语种】中文【中图分类】O322Wilson-θ 法和 Newmark-β 法都是动力学方程数值计算常用的经典方法[1-5]。

多自由度体系wilson-θ法程序编写多自由度体系Wilson-θ法是一种广泛应用于多体动力学和结构动力学领域的数值计算方法。

本文将介绍如何使用Python编程语言编写多自由度体系Wilson-θ法的程序。

一、引言多自由度体系Wilson-θ法是一种基于有限元法的数值计算方法，适用于求解多体动力学和结构动力学中的问题。

该方法通过将体系分解为一系列有限元子系统，并采用θ矩阵方法进行求解，能够有效地处理大规模的多自由度体系。

二、程序编写1. 导入必要的库和模块在编写程序之前，需要导入必要的库和模块，包括numpy、scipy 和matplotlib等。

这些库提供了必要的数学运算、数值分析和图形绘制等功能。

2. 定义体系结构和有限元节点首先需要定义多自由度体系的结构和有限元节点的位置。

可以使用网格划分工具将体系划分为有限元网格，并定义每个节点的位置和编号。

3. 构建有限元矩阵和求解器使用Wilson-θ法进行数值计算，需要构建有限元矩阵和求解器。

该矩阵可以采用三角矩阵的形式进行表示，并使用θ矩阵方法进行求解。

在程序中，需要实现矩阵的构建、求解器的初始化等操作。

4. 迭代求解体系响应使用构建好的矩阵和求解器，可以进行迭代求解多自由度体系的响应。

在每次迭代中，需要输入当前时刻的体系响应作为初值，并输出下一时刻的响应结果。

5. 结果可视化最后，可以使用matplotlib等库将求解得到的响应结果进行可视化。

可以将时间历程、振型、频率响应等结果进行绘制，以便更好地分析体系的动态特性。

三、示例代码以下是一个简单的示例代码，用于演示如何使用Python编程语言编写多自由度体系Wilson-θ法的程序。

代码中假设体系由3个自由度的弹簧-质量系统组成，采用三角形矩阵进行求解。

```pythonimport numpy as npfrom scipy.sparse import csc_matrix, dia_matriximport matplotlib.pyplot as plt# 定义体系结构和有限元节点nodes = np.array([[0], [0.5], [1]]) # 节点位置数组degrees_of_freedom = 3 # 自由度数量system_size = len(nodes) # 体系大小node_indices = np.arange(system_size) # 节点编号数组# 构建有限元矩阵和求解器theta_matrix = csc_matrix(dia_matrix(system_size - degrees_of_freedom, 0)) # θ矩阵mass_matrix = csc_matrix(np.diag([0.5, 0.5, 1])) # 质量矩阵solution = np.zeros((system_size, degrees_of_freedom)) # 初始响应数组forces = np.zeros((system_size, degrees_of_freedom)) # 输出力数组forces[:degrees_of_freedom] = np.zeros((system_size, degrees_of_freedom)) # 初始输出力数组为零向量solver = theta_matrix.dot(solution) +theta_matrix.dot(forces) + mass_matrix # 初始化求解器theta_vector = np.zeros(system_size) # θ向量用于控制有限元矩阵的构造和更新# 进行迭代求解体系响应for iteration in range(100): # 迭代次数限制为100次response = solver.dot(theta_vector) # 输入当前时刻的响应作为初值进行迭代求解下一时刻的响应结果输出为力向量output_forces在每个节点上作用在体系的上结果可与theta向量用于控制有限元矩阵的构造和更新为了演示程序的基本结构和流程以上给出了一个简单的示例代码其中包含的主要内容有定义体系结构和有限元节点构建有限元矩阵和求解器以及进行迭代求解体系响应结果可视化等当然在实际应用中可能还需要考虑更多的因素例如如何处理边界条件如何处理体系的非线性特性等等因此在实际应用中需要根据具体问题对程序进行适当的修改和优化以下是一些可能需要的注意事项和技巧：1. 选择合适的有限元网格划分工具和算法，以确保计算的精度和效率。