基于 MATLAB 实现对结构动力响应的几种算法的验证
基于MATLAB实现对结构动力响应的几种算法的验证
基于 MATLAB 实现对结构动力响应的几种算法的验证1. 算例首先，本文给出一算例， 结构在外力谐振荷载 P (t ) = P 0 sin θt 作用下，分别利用理论解法，杜哈梅积分， Wilson-θ 法求出该结构的位移时程反应。
其中：m = 3.5×103 kg ， P 0 = 1.0×104 N ， k =1.3584515×107 ，ξ=0.05 ，θ=52.3s −1 ，ω=62.3s −1 ，⋅D ω= ω 1-ξ2=62.222 ，初始位移、速度v (0) = 0 ，v (0) = 0 ；2. 算法验证2.1 理论解法运动方程为： mv+cv+kv=0P sin t ϑ由线性代数解出其理论解为:]cos 2)sin -[(]4)-[(t]sin )-(2cos [2]4)-[(t]sin )0()0(cos )0([2222222202222222222t t m P t m P e v v t v e v D DD tD DD t θξωθθθωθωξθωωωθωωξωξωθωξθωθωωξωωξωξω-++-+⋅++++=--由于初始位移v(0) =0 ，v(0) =0 ；则：]cos 2)sin -[(]4)-[(t]sin )-(2cos [2]4)-[(2222222202222222222t t m P t m P e v D DD tθξωθθθωθωξθωωωθωωξωξωθωξθωθξω-++-+⋅+=-v(t ) =-3.115te⋅ 1.05269898×410-⋅[6.230cos62.222t −18.106sin 62.222t]+2.012808757 610-⨯⋅[1146sin 52.3t −325.829cos52.3t]可以用 MA TLAB 进行编程分析，画图位移时程图，详细程序见附录。
2.2Wilson-θ法Wilson-θ法是Wilson 于1966年基于线性加速度法的基础上提出一种无条件收敛的计算方法。
结构动力响应数值计算方法对比分析
结构动力响应数值计算方法对比分析作者：李涵来源：《青年生活》2019年第21期摘要：中心差分法、纽马克法、威尔逊-法是结构动力学中常用的三种方法，为了系统的比较其优缺性，本文针对一个双自由度的体系，首先根据已知条件计算出振动微分方程，运用Matlab计算出可求出12个步长内相应的位移值，即精确解。
然后分别运用中心差分法，纽马克法，威尔逊-法求出其近似解;最后通过三种方法的近似解与精确解相对比，进而分析出三种计算方法的优缺性，为结构动力计算提供依据。
关键词：动力计算、中心差分法、纽马克法、威尔逊-法1、动力体系概况2、精确解推导针对该双自由度体系，理论推导出系统的位移表达式，通过代入各时刻周期得出位移在各时刻的具体数值，即位移精确解。
对位移方程求一阶导数得出速度方程，求二阶导数求出加速度方程。
代入各时刻的周期值，通过Matlab计算得出位移、速度、加速度的数值如下：3、三种数值计算方法3.1、中心差分法中心差分法是基于用有限差分代替位移对时间的求导，对位移一阶求导得到速度，对位移二阶求导得加速度。
通过Matlab计算出前4个步长所对应的位移响应。
3.2、纽马克法纽马克-β法是一种将线性加速度方法普遍化的方法。
通过Matlab计算出前4个步长所对应的位移响应。
3.3威尔逊-法通过Matlab计算出前4个步长所对应的位移响应。
4、近似解与精确解对比分析从上述结构的位移、速度、加速度可以看出，三种方法都能大致表示该体系大体运动趋势，并且误差较小。
其中，在描述物体位移时，中心差分法较后两种方法更为精确。
然而在描述速度和加速度时，中心差分法表现出了较大的误差，而纽马克和威尔逊法则能更详尽的表征物体速度和加速度。
5、结论中心差分法、纽马克法和威尔逊-θ法均是结构动力计算中的常用方法。
本文针对具体的计算实例，分别计算出三种方法的动力响应结果，并与精确解进行对比。
经过分析，中心差分法能更精确的表示物体位移响应，而纽马克和威尔逊法在表征物体速度和加速度方面相较于中心差分法更为精确，三种方法，各有其优缺点，应视具体情况采用相应的计算方法。
基于 matlab 对结构动力特性分析的求解
山 西 建 筑
Ｓ ＨＡＮ ＸＩ ＡＲＣ ＨＩ Ｔ ＥＣ Ｔ ＵＲ Ｅ
Ｖｏ １ ． ４１ Ｎｏ． １ ５
Ｍａ ｙ ． ２ ０ １ ５
・ ４ ９・
文章编号 ： １ ０ ０ ９ － ６ ８ ２ ５ 【 ２ ０ １ ５ ） １ ５ － ０ ０ ４ ９ － ０ ２
点层 的串联刚片系模型 ， 每个节点层考虑 ２个方 向的 自由度 （ ２个 柱的试验研 究及 承载 力 分析 ［ Ｊ ］ ． 土 木 工程 学报 ， ２ ０ ０ ７ ， ４ ０
（ ３ ） ： ２ ５ — ３ １ ．
及 变形性 能研 究［ Ｊ ］ ． 混凝土 ， ２ ０ ０ ２ （ １ ０ ） ： １ ０ ５ — １ ０ ７ ．
计算 结果进行 比较 ， 说明 了此种方 法使 用的范围和缺 陷。 关键 词 ： 结构 动力 ， 柔度 法 ， 计算 ， 周期 中图分类号 ： Ｔ Ｕ ３ １ １ ． ３ 文献标识 码 ： Ａ
对 于高耸结构而 言 ， 结构 的周期 和振 型是进行结 构计 算 的重 ２ 计 算 算例
要步骤 。通常结构设 计 时所建 立 的三维 精细 有 限元模 型 自由度 数量极 为巨大 ， 在对结 构进 行 随机顺 风 向风振 计算 时极 为 不便 ， 从结构 动力分析 的角度 来看 ， 是 不 必要 的 ， 因此需 要对 模 型进 行
１８ １ ７
某地 一座 ９ ９ ｍ高发 射塔 ， 塔 总图如图 １ 所示 。
的关 系 ， 由此 可得 到 频 率 方 程 如 下 ：
ｌ ［ ６ ］ ［ Ｍ］一 Ａ［ ， ］ Ｉ ＝ ０
１ ∞ －
（ ２ ）
一
０ ． ０ ０ ０
动力学系统时域响应计算的四种方法Matlab源程序(接上)
Fu=varargin{3};
t=varargin{5};
[n,n]=size(M);[n,m]=size(Fu);
nstep=size(t');
[V,D]=eig(K,M);
[lambda,k]=sort(diag(D)); % Sort the eigenvalues and eigenvalues
% q0,dq0 - Initial conditions
% nummode - the number of extracted modes
% a, b - Parameters for proportional damping [C]=a[M]+b[K]
% Output parameters
error('Incorrect number of input arguments')
else
switch nargin
case 9
%------------------------------------------------------------------------
Modalresponse模态叠加法法Matlab源程序
function varargout=Modalresponse(varargin)
%-------------------------------------------------------------------------
% Solve the eigenvalue problem and normalized the eigenvectors
%-----------------------------------------------------------------------
基于MATLAB的四层框架结构动力响应与研究
暨南大学研究生课程论文课程：结构动力学姓名：许可悦学号：1634361002学院：力学与建筑工程学院专业：建筑与土木工程任课教师：李雪艳基于MATLAB的四层框架结构动力响应与研究许可悦（暨南大学理工学院力学与土木工程学院，广州 51063）摘要：本文用MATLAB语言对四层建筑结构进行编程，计算结构的自振频率、振型，分析该结构在自由振动和一般激励下的动力响应。
采用了Newmark-β法计算了在简谐正弦激励作用下结构的位移响应，并以此为初始条件结合瑞利阻尼矩阵计算了结构在简谐正弦荷载卸载后的结构自由振动的位移响应。
关键词：MATLAB、Newmark-β法、瑞利阻尼矩The four layers of frame structure dynamic responsebased on MATLAB and researchXu Keyue(Jinan university institute of mechanics and civil engineering department, Guangzhou)Abstract：This paper uses MATLAB language to program the the four layers of frame structure , calculates the self-vibration frequency and vibration mode of the structure, and analyzes the dynamic response of the structure under free vibration and general excitation. Adopted the Newmark - beta method to calculate the displacement of the structure under the action of a harmonic sine excitation response, and the initial conditions in combination with the Rayleigh damping matrix to calculate the structure in the structure of harmonic sine load after unloading free vibration displacement response.Key words：MATLAB; Newmark-βmethod;Rayleigh orthogonal damping1 引言在社会发展的今天，很多科技人员都会遇到数值分析计算机应用等问题，一些传统的高级程序语言如FORTRAN 等虽然能在一定程度上减轻计算量,但它们要求应用人员要具有较强的编程能力和对算法有深入的研究. 另外,在运用这些高级程序语言进行计算结果的可视化分析及图形处理方面,对非计算机专业的普通用户来说,存在着很大的难度. MATLAB 正是在这一应用要求背景下产生的数学类科技应用软件。
基于MATLAB的Newmark-β动力反应数值分析法的精度稳定性分析
５编程计算并绘制时程 曲线 分别编 写当 ｓ ０ １ ＝ ． Ｙ＝ ／ Ｂ＝ ／ 、／ １２ １４ １６时计算 ｕ ｖ ｗ的主程序 和 当 ｓＯ ３ Ｙ １２ ＝ ． ２ ＝ ／ Ｂ＝ ／ 、 ／ １４ ｉ６计算 Ｕ ｖＷ的主程序 。程序清单 （ 。 略）
Ｍ Ｔ Ａ 是矩阵实验 室（ａ ｒ ｘ Ｌ ｂｒ ｔｒ ） Ａ ＬＢ Ｍ ｔ ｉ ａｏ ａ ｏ ｙ 的简称 ， 是美 国 Ｍ ｔ Ｗ ｒ ｓ ａ ｈ ｏｋ 公司出品的商业数学软件 ， 于算法 开发、 用 数据可视 化、 数据分析 以及数 值
Ｔ＝ ＊ ｉ （／ ） ０ ５ ｎ２ ｐ ／ ｋ ｍ ‘．
Ｔ ：．３１ ｎ Ｏ ５ ４
当 Ｂ １４ ＝ ／ 时算法无条件稳定 ；当 Ｂ １ ６ ＝ ／ 时，算法稳定性条件 为 ｓ ≤
ｆ Ｏ５ Ｔ ＝ ． ５１ ｎ。
计算稳定性条件 ｆ ：
ｆ ０ ５ １Ｔ ＝ ． ５ ｎ ｆ０２ ３ ＝ ． ９４
然后分 别绘制 时程位移 曲线、 时程速度 曲线和绘制时程加速度 曲线 ， 并
计算 的高级技术计算语言和交互式环境。２ 世 纪 ７ 年代 ， 国新墨西哥大 Ｏ Ｏ 美
学计算机科学系主任 Ｃ ｅｅ Ｍ ｌｒ为 了减轻学生编程 的负担，用 ＦＲ Ｒ Ｎ ｌｖ ｏ ｅ Ｏ ＴＡ 编写 了最早的 Ｍ ＴＡ 。１ ８ ｉ ｆ ＬＢ ９ ４年由 Ｌ ｔ ｌ 、ｏ ｅ 、 ｔｖ ａ ｇ ｒ ｉ ｔｅ Ｍ ｌｒ Ｓ ｅ ｅ Ｂ ａ ｅ ｔ合作成立 了的 Ｍ ｔＷ ｒ ｓ公司正式把 Ｉ ＬＢ推 向市场 。到 ２ ａｈ ｏ ｋ ￣Ｔ Ａ Ｏ世纪 ９ Ｏ年代 ，Ａ ＬＢ Ｍ Ｔ Ａ
MATLAB中常见的结构动力学分析技巧
MATLAB中常见的结构动力学分析技巧引言：结构动力学是工程领域中重要的研究分支，它主要涉及结构物在外界作用下的运动和响应。
而在 MATLAB 软件中，我们可以通过各种功能强大的工具和函数来进行结构动力学的分析和模拟。
本文将介绍一些 MATLAB 中常见的结构动力学分析技巧，帮助读者更好地利用 MATLAB 进行结构工程相关研究。
一、模型建立与导入1. 建立结构模型在 MATLAB 中，我们可以使用各种方法建立结构模型，比如使用节点和单元建立有限元模型。
通过确定节点的坐标和连接关系，我们可以使用有限元方法对结构进行分析和计算。
除了有限元法，还有其他建模方法，如刚体或连续参数模型等，可以根据实际需要选择。
2. 导入外部模型如果我们已经在其他软件中建立好了结构模型，并希望在 MATLAB 中进行进一步的分析，可以通过导入外部模型来实现。
在 MATLAB 中，可以使用函数如“importgeometry”或“importFiniteElementModel”等，将已有的模型导入到 MATLAB 中进行后续处理。
二、动力学分析1. 自由振动分析自由振动是指结构在没有外力作用下的振动情况，通过这种分析可以得到结构的固有频率和模态形式。
在 MATLAB 中，我们可以使用函数“eig”或“eigs”来计算结构的特征值和特征向量。
进一步，通过可视化这些特征向量，我们可以观察到结构在不同固有频率下的振动模式。
2. 强迫振动响应分析强迫振动响应分析是指结构在外力作用下的振动情况，可以通过求解结构的动力学方程来获得。
在 MATLAB 中，我们可以使用函数如“ode45”、“ode23”等，通过数值解法求解二阶或高阶动力学方程。
这些函数可以根据结构的特定运动方程和边界条件，计算出结构的响应。
三、频域分析1. 傅里叶变换傅里叶变换是一种将信号从时域转换到频域的方法。
在结构动力学中，我们可以使用傅里叶变换来分析结构的振动特性。
结构动力学使用中心差分法计算单自由度体系动力反应的MATLAB程序
中心差分法计算单自由度体系动力反映的报告前言基于叠加原理的时域积分法与频域积分法一样，都假设结构在在全部反应过程中都是线性的。
而时域逐步积分法只是假设结构本构关系在一个微小的时间步距内是线性的，相当于分段直线来逼近实际的曲线。
时域逐步积分法是结构动力问题中研究并应用广泛的课题。
中心差分法是一种目前发展的一系列结构动力反应分析的时域逐步积分法的一种，时域逐步积分法还包括分段解析法、平均常加速度法、线性加速度法、Newmarket−β和Wilson−θ法等。
中心差分法(central difference method)原理中心差分法的基本思路将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。
中心差分法是一种显示的积分法，它基于用有限差分代替位移对时间的求导（即速度和加速度）。
如果采用等时间步长，∆t i=∆t（∆t为常数），则速度与加速度的中心差分近似为u i=u i+1+u i−12∆t(1)üi=u i+1−2u i+u i−1∆t2(2)用u表示位移，离散时间点的运动为：u i=u(t i)，u i=u̇(t i)，u i=ü(t i)(i=0,1,2…)体系的运动方程为mü(t)+cu̇(t)+ku(t)=P(t)(3)将速度和加速度的差分近似公式(1)和(2)代入(3)中得出在t i时刻的运动方程，将方程整理得到u i+1由u i 和u i−1表示的两步法的运动方程(4)：(m ∆t2+c2∆t)u i+1=P i−(k−2m∆t2)u i−(m∆t2−c2∆t)u i−1(4)这样就可以根据t i及以前的时刻的运动求得t i+1时刻的运动。
中心差分法属于两步法，用两步法计算时存在起步问题，必须要给出相邻的两个时刻的位移值，才能逐步计算。
对于地震作用下结构的反应问题和一般的零初始条件下的动力问题，可以用(4)直接计算，因为总可以假设初始的两个时间点（一般取i=0,−1）的位移等于零。
用matlab编程实现法计算多自由度体系的动力响应
用matlab 编程实现Newmark -β法计算多自由度体系的动力响应姓名:学号: 班级: 专业:用matlab 编程实现Newmark -β法 计算多自由度体系的动力响应一、Newmark -β法的基本原理Newmark-β法就是一种逐步积分的方法,避免了任何叠加的应用,能很好的适应非线性的反应分析。
Newmark-β法假定:t u u u u t t t t t t ∆ββ∆∆]}{}){1[(}{}{+++-+=&&&&&&(1-1)2]}{}){21[(}{}{}{t u u t u u u t t t t t t ∆γγ∆∆∆+++-++=&&&&& (1-2)式中,β与γ就是按积分的精度与稳定性要求进行调整的参数。
当β=0、5,γ=0、25时,为常平均加速度法,即假定从t 到t +∆t 时刻的速度不变,取为常数)}{}({21t t t u u ∆++&&&&。
研究表明,当β≥0、5, γ≥0、25(0、5+β)2时,Newmark-β法就是一种无条件稳定的格式。
由式(2-141)与式(2-142)可得到用t t u ∆+}{及t u }{,t u}{&,t u }{&&表示的t t u ∆+}{&,t t u ∆+}{&&表达式,即有t t t t t t t u u t u u tu }){121(}{1)}{}({1}{2&&&&&----=++γ∆γ∆γ∆∆ (1-3)t t t t t t t u t u u u t u }{)21(}){1()}{}({}{&&&&∆γβγβ∆γβ∆∆-+-+-=++ (1-4)考虑t +∆t 时刻的振动微分方程为:t t t t t t t t R u K u C u M ∆∆∆∆++++=++}{}]{[}]{[}]{[&&& (1-5)将式(2-143)、式(2-144) 代入(2-145),得到关于u t +∆t 的方程t t t t R u K ∆∆++=}{}]{[ (1-6)式中][][1][][2C t M t K K ∆γβ∆γ++=)}{)12(}){1(}{]([)}){121(}{1}{1]([}{}{2t t t t t t t t u t u u t C u u t u tM R R &&&&&&∆γβγβ∆γβγ∆γ∆γ∆-+-++-+++=+求解式(2-146)可得t t u ∆+}{,然后由式(2-143)与式(2-144)可解出t t u∆+}{&与t t u ∆+}{&&。
状态转移矩阵法求结构响应matlab
【结构响应的状态转移矩阵法求解及其在Matlab中的应用】1. 引言在工程结构领域中，研究结构在外部荷载作用下的响应是十分重要的。
状态转移矩阵法是一种常用的分析方法，能够有效地求解结构的响应。
本文将从状态转移矩阵法的基本原理入手，结合Matlab编程实例，探讨其在求解结构响应中的应用。
2. 状态转移矩阵法的基本原理在动力学分析中，结构的状态可以用位移和速度表示。
状态转移矩阵法通过建立结构的运动微分方程组，进而构建状态转移矩阵，从而描述结构在不同状态下的响应。
其基本原理是根据结构的自由振动特性和外部荷载作用下的影响，推导结构的运动方程，再通过求解状态转移矩阵，得到结构在任意时刻的响应。
3. 在Matlab中求解结构响应在Matlab中，可以利用状态转移矩阵法求解结构的响应。
需要建立结构的运动微分方程组，考虑其初始条件和外部荷载，然后通过Matlab编程计算状态转移矩阵，最后得到结构在不同时刻的位移、速度和加速度响应。
在具体编程实现过程中，需要考虑结构的参数、初始条件和外部荷载等因素，以及求解微分方程组的数值方法和稳定性。
4. 案例分析以某桥梁结构为例，假设结构参数已知，初始条件和外部荷载也已给定。
利用Matlab编程实现状态转移矩阵法，求解结构在不同时刻的位移和速度响应。
通过分析结果，可以得到结构的动态特性，包括振动频率、模态形态和响应时程等。
也可以对不同外部荷载情况下的结构响应进行比较分析，为结构设计和安全评估提供依据。
5. 总结与展望本文从状态转移矩阵法的基本原理入手，结合Matlab编程实例，探讨了其在求解结构响应中的应用。
通过对状态转移矩阵法的理论分析和具体案例的应用分析，可以更加全面、深刻地理解结构的动态特性和响应规律。
在未来的研究中，可以进一步探讨结构响应的优化控制方法，提高结构的动态性能和安全可靠性。
6. 个人观点与理解对于结构动力学分析和响应求解，状态转移矩阵法是一种有效的方法。
结合Matlab编程，可以实现结构响应的精确计算和可视化展示，为工程实践和科研研究提供了重要的技术支持。
Matlab在结构动力学与振动控制中的应用案例
Matlab在结构动力学与振动控制中的应用案例引言结构动力学研究了物体在受力作用下的运动规律，而振动控制则关注如何通过各种手段来控制结构的振动。
在过去的几十年里，Matlab作为一款功能强大的数值计算和数据可视化工具，被广泛地应用于结构动力学与振动控制领域。
本文将通过一些典型的案例，探讨Matlab在这些领域中的应用。
案例一：辛普森建筑物模型辛普森建筑物模型是用于研究地震对建筑物结构的影响的经典案例。
在这个模型中，建筑物底部通过弹簧与地面相连，顶部有一个质量为m的质点。
通过求解二阶常微分方程，在Matlab中可以得到建筑物的振动响应。
通过修改建筑物的初始参数和地震输入信号，我们可以得到不同条件下的振动响应，并进一步分析建筑物的安全性能。
Matlab提供了一系列用于求解常微分方程的函数，如ode45和ode15s等。
结合Matlab的图形界面，我们可以方便地可视化建筑物的振动响应。
通过修改建筑物模型的材料参数、形状和地震输入，我们可以直观地感受到这些因素对振动响应的影响，从而为结构的设计和改进提供参考。
案例二：滑模控制器的设计滑模控制是一种常用的控制方法，在结构振动控制中也被广泛应用。
在滑模控制中，通过引入一个滑模面，使得系统状态在滑模面上快速地滑动，从而实现对系统的控制。
在振动控制中，我们常用滑模控制器来实现结构的抑制和消除。
在Matlab中，我们可以借助控制系统工具箱，便捷地设计和分析滑模控制器。
通过建立结构的数学模型，并在Matlab中使用滑模控制器设计函数，我们可以得到系统的闭环响应，并评估控制器的性能指标，如响应时间、超调量和控制能力等。
在实际应用中，我们可以结合传感器和执行器等硬件装置，与Matlab相结合，实现闭环控制。
这为我们实现各种振动控制策略提供了一个方便而高效的平台。
案例三：有限元分析与模态分析有限元分析是结构工程中常用的一种分析方法，用于预测结构的应力、变形和振动等特性。
在Matlab中，有限元分析可以通过编写相应的程序实现。
应用MatLab软件探讨结构动力响应时域和频域数值模拟教学
应用MatLab软件探讨结构动力响应时域和频域数值模拟教学1. 引言结构动力学是工程力学中的重要分支，在现代工程设计中起着至关重要的作用。
传统的结构动力学教学方法主要基于理论推导和实验验证，然而这种教学方法存在时空限制和人力资源限制的问题。
随着计算机技术和数值模拟技术的不断发展，提高结构动力响应数值模拟技能已经成为了现代工程设计的重要课题，MatLab作为一种常用的科学计算软件，因其性能优良、易于使用和广泛应用而备受青睐。
本文借助MatLab软件，通过对结构动力响应的时域和频域数值模拟探究，讨论如何优化结构动力学的教学方法。
2. 时域数值模拟时域数值模拟是研究结构动力响应的常用方法，其原理是通过分析机械系统的运动学和动力学方程，在时域内求解结构响应的变化规律，这种方法对于研究结构响应的时间变化非常有帮助。
MatLab软件作为一种功能强大的科学计算软件，拥有丰富的函数库和数据分析工具，能够支持多种结构动力响应的时域数值模拟。
通过MatLab中的ode45函数，结构动力响应的时域数值模拟可以轻松完成，可以设置不同的初值条件和参数来模拟不同的结构响应情况。
3. 频域数值模拟频域数值模拟是研究结构动力响应的另一种常用方法，其原理是将结构响应在频域内进行分析，分析结构系统在不同频率下的动态响应特性，这种方法对于研究结构响应的频率特性非常有帮助。
MatLab软件作为一种功能强大的科学计算软件，同样能够支持多种结构动力响应的频域数值模拟。
通过MatLab中的FFT函数，结构动力响应的频域分析可以轻松完成，同时还可以使用MatLab的模拟仿真工具Simulink，进行结构动力响应的实时仿真。
4. 结构动力学教学优化通过MatLab软件进行结构动力响应的数值模拟，不仅可以提高学生对结构动力学教学的理解，同时也能够培养学生的计算机仿真技能和科学研究能力。
综合来看，基于MatLab软件进行结构动力学的教学方法能够更好地满足现代工程设计中对结构动力学技能的需求，是一种具有广阔前景的教学改革方法。
应用MatLab软件探讨结构动力响应时域和频域数值模拟教学
应用MatLab软件探讨结构动力响应时域和频域数值模拟教学【摘要】本文探讨了应用MatLab软件进行结构动力响应时域和频域数值模拟教学的相关内容。
在介绍了研究背景、研究意义和研究目的。
接着在详细讨论了MatLab软件在结构动力响应研究中的应用，时域和频域数值模拟教学的探讨，以及通过结构动力响应数值模拟案例分析来展示方法和应用。
最后对教学效果进行了评估。
结论部分总结了MatLab 软件在结构动力响应教学中的价值，并展望了未来的研究方向。
通过本文的研究，可以更好地掌握MatLab软件在结构动力响应领域的应用，提高教学效果和研究水平，促进结构动力学领域的发展。
【关键词】MatLab软件、结构动力响应、时域、频域、数值模拟、教学、案例分析、教学效果评估、教学价值、研究展望、结论总结。
1. 引言1.1 研究背景结构动力学是研究结构在外力作用下的振动特性和响应行为的一门学科，其研究对于评估结构工程的稳定性、安全性和耐久性具有重要意义。
随着科学技术的发展和应用需求的增加，人们对结构动力响应的研究需求也越来越迫切。
结构动力响应的研究不仅有助于优化结构设计和改进结构性能，还可以为减灾防灾和工程管理提供重要参考。
本研究旨在探讨利用MatLab软件进行结构动力响应的时域和频域数值模拟教学，旨在提高学生对结构动力学理论的理解和实践能力，促进结构动力响应研究的发展和应用。
1.2 研究意义结构动力响应是结构工程领域一个重要的研究方向，涉及到结构在外力作用下的振动响应。
结构动力响应的研究对于提高结构的安全性、可靠性及性能具有重要意义。
通过对结构动力响应进行研究，可以更好地了解结构在振动状态下的性能表现，从而为结构设计、改进和维护提供参考。
探讨MatLab软件在结构动力响应研究中的应用具有重要的研究意义。
通过研究MatLab软件在结构动力响应数值模拟教学中的应用，能够提高学生对结构振动与响应的理解，培养学生的实际应用能力，促进结构动力学教育的发展和完善。
使用MATLAB进行地震响应分析与结构优化
使用MATLAB进行地震响应分析与结构优化地震响应分析与结构优化是地震工程中的重要研究内容，通过使用MATLAB这一强大的工具，可以对结构在地震作用下的响应行为进行深入分析，并进一步进行结构的优化设计。
本文将探讨如何使用MATLAB进行地震响应分析与结构优化的步骤和方法，并通过实例进行说明。
一、地震响应分析地震响应分析是指在地震波作用下，结构的应力、应变和位移等响应情况的分析。
地震波是地震时由震源产生的地震能量在地球内部传播所产生的波动。
在进行地震响应分析时，需要先选取合适的地震波记录作为输入，然后进行模拟计算，并获取结构的响应结果。
在MATLAB中进行地震响应分析，可以利用其中的工具箱和函数。
首先，需要将地震波记录进行读取和处理，可以使用MATLAB中的文件读取函数和信号处理函数。
读取地震波记录后，可以进行滤波和降尺度处理，以使其适用于结构响应分析。
接着，需要进行指定结构的建模和定义，可以使用MATLAB中的结构类对象进行建模，并设置结构的材料和截面属性。
之后，可以利用地震波输入对结构进行地震动力分析，得到结构的响应结果。
二、结构优化设计结构优化设计是指通过调整结构的材料和几何参数，使其满足特定的性能要求，以达到最佳的设计效果。
在地震工程中，结构的优化设计常常与地震响应特性有关，可以通过优化设计来提高结构的地震抗震性能。
在MATLAB中进行结构优化设计，可以利用其中的优化工具箱和函数。
首先，需要确定设计目标和约束条件，例如最小化结构的质量或最大化结构的刚度，并设置约束条件，如保持结构的安全性能。
接着，可以利用MATLAB中的优化算法进行设计优化，例如遗传算法、模拟退火算法等。
在迭代优化过程中，可以通过定义适应度函数来评估设计的优劣，使优化算法能够搜索到最优解。
三、综合实例为了更好地说明使用MATLAB进行地震响应分析与结构优化的方法和步骤，下面将以一座钢筋混凝土框架结构为例进行说明。
首先，需要在MATLAB中导入地震波记录，并进行预处理，包括去除基线漂移、滤波以及降尺度处理。
MATLAB中常见的结构动力学分析技巧
MATLAB中常见的结构动力学分析技巧在MATLAB中进行结构动力学分析时，常见的技巧如下：1. 建立模型：使用MATLAB的建模工具箱，如Simulink或Simscape，可以方便地建立结构动力学模型。
这些工具箱提供了丰富的图形界面和预定义的组件，使得建立模型更加简单和直观。
2. 选择合适的数值方法：对于结构动力学分析，常用的数值方法包括有限元法、边界元法和拉格朗日法等。
在MATLAB中，可以利用有限元分析工具箱（FEATool）或边界元分析工具箱（BEMTool）来实现这些数值方法的应用。
此外，MATLAB还提供了一系列数值方法的函数库，如ode45、ode23等，可以用于解动力学方程。
3. 系统的模态分析：模态分析是结构动力学中的一项重要技术，用于确定系统的固有频率和振型。
MATLAB提供了现成的函数库，如eig、eigs等，用于求解系统的特征值和特征向量。
通过分析这些特征值和特征向量，可以确定系统的固有频率和振型，进而了解结构的动态特性。
4. 频谱分析：频谱分析是结构动力学中常用的分析方法之一，用于确定结构在不同频率下的响应。
在MATLAB中，可以利用FFT（快速傅里叶变换）函数和相关的函数库，如pwelch、spectrogram等，对信号进行频谱分析。
这些函数库提供了多种频谱分析方法，如功率谱密度、频谱图和甚至是瀑布图等。
5. 响应分析：响应分析是结构动力学中常用的分析方法之一，用于确定结构在特定外力作用下的响应情况。
在MATLAB中，可以使用ode函数和相关的函数库，如lsim、bode等，对结构进行响应分析。
这些函数库提供了多种分析方法，如自由振动、强迫振动和频率响应等。
6. 参数识别：参数识别是结构动力学中的关键技术之一，用于确定模型的未知参数。
在MATLAB中，可以使用系统辨识工具箱（System Identification Toolbox）进行参数识别。
该工具箱提供了多种参数识别方法，如最小二乘法、极大似然估计法和频率域法等。
基于MATLAB的结构动力特性分析
第１ 期
马志贵 ： 基 于 ＭＡ Ｔ Ｌ Ａ Ｂ的结构动力特性分析
＝
・１ ９・
２ 用 ＭＡ Ｔ Ｌ Ａ Ｂ建 立 平 面 梁 单 元 的 刚 度 矩 阵及 质 量 矩 阵
要求解式（ ６ ） ， 必须组集结构总刚度矩阵和结 构总质量矩阵 ， 用Ｍ Ａ Ｔ Ｌ Ａ Ｂ可 以方便的构造梁单
Ｆ ｅ ｂ ． ２ ０ １ ４
基 于 ＭＡ Ｔ Ｌ ＡＢ的 结构 动 力特 性 分 析
马 志贵 ， 刘世 忠
（ 兰州交通大学 土木工程学院 ， 甘肃 兰州 ７ ３ ０ ０ ７ ０ ）
摘要 ： 介 绍 了结构 动力特 性 分析 的基本 原理和 方 法 ， 利用 Ｍ Ａ Ｔ Ｌ Ａ Ｂ语 言编 写 了结构 自振 频 率 和振
Ｌ ０ ０ ； ０ １ Ｌ ０ Ｌ ＾ ２ Ｌ ３ ； ０ ０ １ ０ ２ Ｌ ３ Ｌ ２ ］ ； Ｈｕ ＝［ １ ０ ０ Ｘ ０ ０ ］ ； Ｈｖ ＝［ ０ １ Ｘ ０ ｘ ＾ ２ ｘ ３ ］ ；
０
，
式（ ２ ７ ） 中， Ｅ为 弹性模 量 ； Ａ为截 面面 积 ； ， 为 截 面对 主轴 的惯性 矩 ； ｚ 为梁单 元 长度 ． 刚度矩 阵 的程 序段 如下 ：
ｓ ｙ ｍｓ Ｅ Ｌ Ｘ Ｙ ｂ ｈ
～
一
，
一
ｏ 一 铆 丁
巨 一 ， 一
型向量的求解程序． 通过一个算例 ， 表 明编写的程序计算效率 高， 计算精度 高于用 Ａ Ｎ Ｓ Ｙ Ｓ所得结
果， 可为求 解结构 自振 频率 和振 型提 供 新 的思路 ．
关 键词 ： 有 限元 ； 动力特 ｌ ＃ － ； ＭＡ Ｔ Ｌ Ａ Ｂ；
中图分 类号 ： Ｔ Ｐ ３ ９ １ ． ７ 文 献标 志码 ： Ａ
基于MATLAB的结构动力响应分析方法
求 Ｄ特征值程序为 ：
１．］ ２２ ； ｌ ］ｃ ］Ｓ ］ 一［ ；一［ ；ｚ［ ；
选取旋转矩阵 Ｐ 一Ｒ ２ １ｅ）使 Ａ ’ ｋ 第 （，， ， ｎ 一Ｐ Ａ 一 列次对角元 ａ ｎ 一０ 再选取旋 转矩阵 Ｐ 一Ｒ（ ， ， 。’ ； 。 ３２ ０）使 Ａ 一ＰＡ（ ） ２， ’ ２ １ 的对 角 元 ａ ’ …… 如此 继 续 ３ 一０ ２
１ 运动方 程
பைடு நூலகம்
若使式 （） ５有非零解 ， 其系数矩 阵行列式值必须为 零。由此可求得 ∞的两个实根 ， 即体系的两个 自振圆 频率。然而对于 自由度数 目 多的结构物， 很 要解 ∞ 的 ｎ 次方程是非常 困难 的。我们 可以把这个问题归结为 求特征值问题。
２ ２ 主振 型 ．
维普资讯
１８ ８
西 部探 矿工 程
２０ ０ ７年第 ４期
基 于 ＭＡＴ Ａ Ｌ Ｂ的结 构 动 力响 应 分 析 方 法
马军庆 ， 姜鹏 飞
（ 放军理 工 大学 工程 兵工程 学 院， 苏 南京 ２００ ） 解 江 １０ ７
摘 要 ： 对结构动 力响应 分析 时 自振频 率的求解 问题 ， 针 简要介 绍 了 Ｑ Ｒ方法 的基本原理 ， 借助 ＭＡ Ｉ Ｂ编程 实现 了结构 自振频率的求解过程。通过算例介绍 了程序 的实现方法， Ｔ Ａ 并与 ＭＡ Ｉ Ｂ Ｔ Ａ
式 中 ： 、 Ｋ — 多 自由 度 体 系 的 质 量 矩 阵 、 尼 矩 Ｍ Ｃ、 — 阻 阵 、 度矩 阵 ； 刚 ｌ ｎ 维单 位列 向量 ； —— ×ｌ
维普资讯
２０ ０７年第 ４期
西部 探矿 工程
１９ ８
Ａ－ 一 ＱＩ１ 一 ， Ａ — Ｒ一 一 ， 中 一 为 正 交 ｋ１ ‘ １令 ｋ ｋ１ １其 一Ｒｋ １ 矩 阵 ， ｋ Ｒ一为上 三 角矩 阵 。这 样 做 成 的 矩 阵 序 列 Ａ 都 相似于 Ａ ， 从而与 Ａ有相同的特征值 。 Ｑ Ｒ方 法分解 Ａ 一 一 ＱＩ Ｒ 一 常 用旋 转 矩 阵 来 完 ‘ 一
用matlab编程实现法计算多自由度体系的动力响应
用matlab编程实现Newmark-β法计算多自由度体系的动力响应姓名：学号：班级：专业：用matlab 编程实现Newmark -β法 计算多自由度体系的动力响应一、Newmark -β法的基本原理Newmark-法是一种逐步积分的方法，避免了任何叠加的应用，能很好的适应非线性的反应分析。
Newmark-法假定：t u u u ut t t t t t ∆ββ∆∆]}{}){1[(}{}{+++-+= (1-1)2]}{}){21[(}{}{}{t u u t uu u t t t t t t ∆γγ∆∆∆+++-++= (1-2) 式中，和是按积分的精度和稳定性要求进行调整的参数。
当=0.5，=0.25时，为常平均加速度法，即假定从t 到t +t 时刻的速度不变，取为常数)}{}({21t t t u u ∆++ 。
研究表明，当≥0.5，≥0.25(0.5+)2时，Newmark-法是一种无条件稳定的格式。
由式(2-141)和式(2-142)可得到用t t u ∆+}{及t u }{，t u }{ ，t u }{ 表示的t t u ∆+}{ ，t t u ∆+}{ 表达式，即有t tt t t t t u u t u u t u}){121(}{1)}{}({1}{2----=++γ∆γ∆γ∆∆ (1-3) t t t t t t t u t uu u t u}{)21(}){1()}{}({}{ ∆γβγβ∆γβ∆∆-+-+-=++ (1-4) 考虑t +t 时刻的振动微分方程为：t t t t t t t t R u K u C uM ∆∆∆∆++++=++}{}]{[}]{[}]{[ (1-5) 将式(2-143)、式(2-144) 代入(2-145)，得到关于u t +t 的方程t t t t R u K ∆∆++=}{}]{[ (1-6)式中][][1][][2C t M tK K ∆γβ∆γ++= )}{)12(}){1(}{]([)}){121(}{1}{1]([}{}{2t t t t t t t t u t uu t C u u t u tM R R ∆γβγβ∆γβγ∆γ∆γ∆-+-++-+++=+求解式(2-146)可得t t u ∆+}{，然后由式(2-143)和式(2-144)可解出t t u ∆+}{ 和t t u∆+}{ 。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
基于 MATLAB 实现对结构动力响应的几种算法的验证1. 算例首先，本文给出一算例， 结构在外力谐振荷载 P (t ) = P 0 sin θt 作用下，分别利用理论解法，杜哈梅积分， Wilson-θ 法求出该结构的位移时程反应。
其中：m = 3.5×103 kg ， P 0 = 1.0×104 N ， k =1.3584515×107 ，ξ=0.05 ，θ=52.3s −1 ，ω=62.3s −1 ，⋅D ω= ω 1-ξ2=62.222 ，初始位移、速度v (0) = 0 ，v (0) = 0 ；2. 算法验证2.1 理论解法运动方程为： mv+cv+kv=0P sin t ϑ由线性代数解出其理论解为:]cos 2)sin -[(]4)-[(t]sin )-(2cos [2]4)-[(t]sin )0()0(cos )0([2222222202222222222t t m P t m P ev v t v e v D DD tD DD t θξωθθθωθωξθωωωθωωξωξωθωξθωθωωξωωξωξω-++-+⋅++++=--由于初始位移v(0) =0 ，v(0) =0 ；则：]cos 2)sin -[(]4)-[(t]sin )-(2cos [2]4)-[(22222222022222222220t t m P t m P e v D DD tθξωθθθωθωξθωωωθωωξωξωθωξθωθξω-++-+⋅+=-v(t ) =-3.115te⋅ 1.05269898×410-⋅[6.230cos62.222t −18.106sin 62.222t]+2.012808757 610-⨯⋅[1146sin 52.3t −325.829cos52.3t]可以用 MA TLAB 进行编程分析，画图位移时程图，详细程序见附录。
2.2Wilson-θ法Wilson-θ法是Wilson 于1966年基于线性加速度法的基础上提出一种无条件收敛的计算方法。
该方法假定在θt ∆时程步长内，体系的加速度反应按线性变化。
对于地震持续时间内的每一个微小时 段t ∆ ，从第一时段开始到最后一个时段，逐一的重复以下计算步骤，即得到结构地震反应的全过程。
下面以第i+1时段（时刻时刻到1+i i t t ）为例：{}{}[]{}[][][][]{}[]{}[]{}{}[]{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}[][]{}[][]{}附录。
移时程图，详细程序见进行编程分析，画图位可以用；时刻的相对加速度反应，求得第代入结构运动微分方程和最后将和速度增量移反应时刻的各质点的相对位得到第，和速度增量，分别加上位移增量时刻的地震反应为基础再以各质点在第；和速度增量内位移增量正常时段再求解第；内的加速度增量正常时段计算出第先利用位移增量的加长时段内各质点的计算MATLAB x xC xx x x x x x xx x x xx x x t xt x t x t x xxt x x x i xxx xx x xx C x x xP C K P K x x i i g i i i i i i i i i i i i i i i i i i i i i g iii )(-t )5(t )4(631)(2t 1)3()2(6t 1i )2()23()36(136;t )1(1i 11i 11i 1i 11i 1i 1i 11i 1i 1122122211,+-+-+++++++++++++***-*K M +M +=∆+=∆+=∆∆∆+∆+∆=∆+∆=∆∆∆∆+--∆=∆∆∆+∆+++M +∆M -=∆+M +K =∆=∆∆∆=+ τττττθτττττθτττττττ2.3 杜哈梅积分杜哈梅积分在考虑阻尼的情况是：ττωτωωωξωωτξωξωd t p e m v v t v e v D tt DD DD t )(sin )(1t]sin )0()0(cos )0([0)(-+++=⎰---可以用 MA TLAB 进行编程分析，画图位移时程图，详细程序见附录。
3. 位移时程反应对比分析利用MATLAB将理论解法，杜哈梅积分，Wilson-θ法求解出来的位移时程反应画在同一张图中，进行比较分析。
从图中可以看出，以上三种方法得出来的位移时程曲线基本吻合，误差基本保持在5%以内，所以以上几种方法在求解相关问题上都具有一定的作用效力。
4.结论本文通过一个简单的单自由度系统动力分析算例（仅作位移分析，其它分析雷同），基于MATLAB，将理论解法，杜哈梅积分法，逐步积分法（本文采用Wilson-θ法）进行相互验证，从最后的位移分析图对比上，可以很好的看出三种方法均能很好的彼此验证，从而说明了三种方法在相关问题上的作用效力。
附录：MA TLAB 源程序%理论解，杜哈梅积分，Wilson-θ法程序clc;clearh1=figure(8);set( h1, 'color','w')%理论解法t=0:0.01:1;v=1.05269898*10^(-4)*exp(-3.115*t).*(6.230*cos (62.222*t)-18.106*sin(62.222*t))+2.012808757*1 0^(-6)*(1146*sin(52.3*t)-325.829*cos(52.3*t)); plot(t,v,'k')hold on%杜哈梅积分法aa=1;%输入时间长度bb=0.01;%输入精度t=bb:bb:aa;t1=t;theta=52.3;%输入荷载频率w=62.3;%输入自振频率m=3500;%输入质量p0=10000;%输入荷载幅值p0=p0*ones(1,aa/bb);p=p0.*sin(theta*t);%荷载函数for i=1:(aa/bb)for j=1:icanshu1(j)=p(j)/(m*w)*bb*sin(w*(t(i)-t1(j)));%杜哈梅积分中的被积函数endy(i)=sum(canshu1);%%位移值endfor i=1:aa/bb-1v1(i)=(y(i+1)-y(i))/bb;%计算速度endfor i=1:(aa/bb-2)a(i)=(v1(i+1)-v1(i))/bb;%计算加速度endhold onplot(t,y,'r')%画位移图hold on%Wilson-θ法dt=0.01;ct=1.4;ndzh=100;k=13584515;c=21805;t=0:dt:ndzh*dt;ag=10000*sin(52.3*t);ag1=ag(1:ndzh);ag2=ag(2:ndzh+1);agtao=ct*(ag2-ag1);wyi1=0;sdu1=0;jsdu1=0;wyimt=0;sdumt=0;jsdumt=0;for i=1:ndzhkxin=k+(3/(ct*dt))*c+(6/(ct*ct*dt*dt))*m;%kxin为新的刚度dpxin=-m*agtao(i)+m*(6/(ct*dt)*sdu1+3*jsdu1)+c*(3*sdu1+ct*dt/2*jsdu1); %新的力增量dxtao=kxin\dpxin;dtjsdu=6*dxtao/(ct*(ct^2*dt^2))-6*sdu1/(ct*ct*dt)-(3/ct)*jsdu1;jsdu=jsdu1+dtjsdu;dtsdu=(dt/2)*(jsdu+jsdu1);sdu=sdu1+dtsdu;dtwyi=dt*sdu1+(1/3)*dt^2*jsdu1+(dt^2/6)*jsdu;wyi=wyi1+dtwyi;jsdu=-m\(m*ag2(i)+c*sdu+k*wyi); % 调整加速度wyi1=wyi;sdu1=sdu;jsdu1=jsdu;wyimt=[wyimt wyi];sdumt=[sdumt sdu];jsdumt=[jsdumt jsdu];endplot(t,-wyimt/3000,'b');hold offlegend('\fontsize{9}\fontname{ 黑体} 理论解','\fontsize{9}\fontname{黑体} 杜哈梅积分法','\fontsize{9}\fontname{黑体} wilson-{\theta}法')%基于双线性滞回模型的单自由度体系的地震能量分析程序%质量57041kg，阻尼36612 N·s/m,初始刚度2350000N/m,刚度折减系数0.2，屈服位移0.01m，采用ELCENTRO波%参数替换直接在下面修改，然后运行clcformat long;m=57041;%质量ug=importdata('ELCENTRO.txt');%地震波txt 文件ug=ug/100;P=-m*ug;num=size(P,1);c=36612;%阻尼k1=2350000;%初始刚度k2=k1*0.2;a=zeros(1,num);v=zeros(1,num);x=zeros(1,num);%加速度速度位移Ei=zeros(1,num);Ek=zeros(1,num);Ed=zeros(1,num);Eh=zeros(1,num);%输入能动能阻尼耗能滞回耗能EI=zeros(1,num);EK=zeros(1,num);ED=zeros(1,num);EH=zeros(1,num);%累积的各种能量time=zeros(1,num);a(1)=P(1)/m;hfl=zeros(1,num);h=0.02;%地震波采样间隔Xy=0.01;%屈服位移pxmax=0;nxmax=0;pd=0;%双线型滞回模型折线的标识，0表示弹性，1表示正向弹塑性，2表示反向弹性，-1表示反向弹塑性，-2表示正向弹性for ii=1:numif pd==0 %弹性阶段k=k1;if x(ii)>Xypd=1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-Xy];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp+k1*Xy;kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=Xy+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=k1*Xy+dtf;a(ii)=(P(ii)-hfl(ii+1)-c*v(ii))/m;elseif x(ii)<-Xypd=-1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)+Xy];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp-k1*Xy;kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=-Xy+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=-k1*Xy+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;endendi f pd==1 %正向弹塑性k=k2;if v(ii)<0pd=2;b=[(a(ii)-a(ii-1))/2/ha(ii-1) v(ii-1)];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;xp=x(ii-1)+v(ii-1)*dth+a(ii-1)*dth^2/2+(a(ii)-a(ii-1))*dth^3/h/6;ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*0+k1*Xy+k2*(xp-Xy);kd=k1+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*0/hp+3*ap)+c*(3*0+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*0-hp*ap/2;x(ii)=xp+dtx;v(ii)=0+dtv;dtf=k1*dtx;hfl(ii)=k1*Xy+k2*(xp-Xy)+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;pxmax=xp;endendif pd==2 %反向弹性k=k1;if x(ii)>pxmaxpd=1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-pxmax];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp+k1*Xy+k2*(pxmax-Xy);kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=pxmax+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=k1*Xy+k2*(pxmax-Xy)+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;elseif x(ii)<(pxmax-2*Xy)pd=-1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-pxmax+2*Xy];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp+(pxmax*k2-k1*Xy-k2*Xy);kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=pxmax-2*Xy+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=pxmax*k2-k1*Xy-k2*Xy+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;endendif pd==-1 %反向弹塑性k=k2;if v(ii)>0pd=-2;b=[(a(ii)-a(ii-1))/2/h a(ii-1)v(ii-1)];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;xp=x(ii-1)+v(ii-1)*dth+a(ii-1)*dth^2/2+(a(ii)-a(ii-1))*dth^3/h/6;ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*0-k1*Xy+k2*(xp+Xy);kd=k1+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*0/hp+3*ap)+c*(3*0+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*0-hp*ap/2;x(ii)=xp+dtx;v(ii)=0+dtv;dtf=k1*dtx;hfl(ii)=-k1*Xy+k2*(xp+Xy)+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;nxmax=x(ii);endendif pd==-2 %正向弹性k=k1;if x(ii)<nxmaxpd=-1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-nxmax];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp-k1*Xy+k2*(nxmax+Xy);kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=nxmax+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=-k1*Xy+k2*(nxmax+Xy)+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;elseif x(ii)>(nxmax+2*Xy)pd=1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-nxmax-2*Xy];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp+(k1*Xy+nxmax*k2+Xy*k2);kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=nxmax+2*Xy+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=k1*Xy+nxmax*k2+Xy*k2+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;endendi f ii>=numbreak;endkd=k+3*c/h+6*m/h/h;dtpd=P(ii+1)-P(ii)+m*(6*v(ii)/h+3*a(ii))+c*(3*v(ii)+h*a(ii)/2);dtx=dtpd/kd;dtv=3*dtx/h-3*v(ii)-h*a(ii)/2;x(ii+1)=x(ii)+dtx;v(ii+1)=v(ii)+dtv;dtf=k*dtx;hfl(ii+1)=hfl(ii)+dtf;a(ii+1)=(P(ii+1)-hfl(ii+1)-c*v(ii+1))/m;Ei(ii+1)=-m*ug(ii+1)*v(ii+1)*h;Ed(ii+1)=c*v(ii+1)^2*h;Ek(ii+1)=m*a(ii+1)*v(ii+1)*h;Eh(ii+1)=hfl(ii+1)*v(ii+1)*h;for jj=1:num-1EI(jj+1)=EI(jj)+Ei(jj+1);ED(jj+1)=ED(jj)+Ed(jj+1);EK(jj+1)=EK(jj)+Ek(jj+1);EH(jj+1)=EH(jj)+Eh(jj+1);endendfor j=1:numtime(j)=h*(j-1);endfid=fopen('EL波弹塑性各种能.txt','wt'); %输出结果到EL波弹塑性各种能.txt fprintf(fid,'%g %g %g %g%g\n',[time;EI;ED;EH;EK]);fclose(fid);plot(time,EI);hold onplot(time,ED);hold onplot(time,EK);hold onplot(time,EH);。