基于 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）的位移等于零。
- 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);。