下载文档原格式
《高等结构动力学》课程Case Study斜坡缓冲车辆的运动学模型与缓冲距离影响因素分析姓名： 班号： 学号：摘要： 为了计算无动力车辆在斜坡上的最小缓冲距离，本文建立了斜坡行驶车辆的半车模型的运动学方程，采用龙格库达法对微分方程求解，得到初始速度与斜坡角度对缓冲距离的影响规律。
1.引言为防止制动失灵的车辆冲下山谷, 盘山公路的下行方向每隔一定距离需要设置一个缓冲区，如图1所示。
缓冲区一般由一段具有上升坡度的渣土路面形成。
制动失灵的车辆驶入缓冲区后，其动能一部分转换成势能，一部分由车轮与路面的摩擦耗散。
图2所示为一车辆简化模型，车体高h=1.8m, 长b=5m 。
已知前轮刚度K 1=5.5*105N/m, 前轮阻尼系数C 1=8*104N •s/m, 后轮刚度K 2=8.5*105N/m, 后轮阻尼系数C 2=C 1; 车体按匀质记，总重10吨，质心距地面高度H=1.5m 。
摩擦力按下式计算：()()i i f t N t μ=⋅ i=1, 2 μ—摩擦系数，μ=0.3N i -- 车轮所受地面的正压力。
图1 盘山公路缓冲区示意图 图2 车辆简化模型假设: ① 车辆行驶过程中的车体变形很小,可忽略不计。
② 车轮质量与车身质量相比很小,可忽略不计。
分别给出缓冲区坡度为300和450时的车辆驶入速度与缓冲区长度的关系曲线以及车速为70Km/小时时缓冲区的最小长度。
2 斜坡行驶车辆的动力学模型斜坡行驶车辆的物理模型与力学模型分别如图3和图4所示。
图3 斜坡行驶车辆物理模型图4 斜坡行驶车辆力学模型2）.图5 斜坡行驶车辆模型受力分析建立如图5所示的斜坡行驶车辆的力学模型，以质心C 点垂直方向坐标cy 和转角c θ为广义坐标，1y 和2y 分别为弹簧位置垂直方向坐标（均取在弹簧原长的位置处），采用达朗贝尔原理建立车辆运动的微分方程如下。
以C 点垂直斜面方向的力平衡方程：111222()()cos 0c my k y cy k y cy mg α++++-=(1)以C 点沿斜面方向的力平衡方程：111222()sin 0c mx k y cy k y cy mg μα+++++=(2)以质心C 点取矩的力矩平衡方程：111222111222()()()()()()022cc c L LJ cy k y cy k y H y k y cy H y k y cy θμμ-+-++-++-+=+(3)式中，车辆转动惯量22()12m J b h =+；A 点坐标12c c Ly y θ=-，B 点坐标12c c Ly y θ=+，坐标几何关系如图6所示。
结构动力学大作业班级：学号：姓名：目录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编程，源代码见附录。
高等结构动力学大作业引言：高等结构动力学是土木工程中的重要学科，涉及到结构的振动和响应分析。
为了加深学生对该学科的理解和运用能力，设计一份详细具体的大作业是非常有益的。
本文将介绍一个高等结构动力学大作业的设计，包括作业目标、内容和评价方式。
一、作业目标1.理论掌握：通过大作业，学生需要巩固和应用所学的高等结构动力学理论，提高对结构振动和响应分析方法的理解和运用能力。
2.实践能力培养：作业要求学生进行实际案例的分析和计算，培养他们的实践能力和问题解决能力。
3.创新思考：作业鼓励学生从不同的角度进行创新性思考，提出改进或优化现有结构的方法或方案。
4.报告撰写能力：作业要求学生以报告形式呈现研究成果，培养他们的科学写作能力和沟通表达能力。
二、作业内容1.理论分析：作业可以要求学生选择一个特定的结构，如悬索桥、高层建筑等，进行结构振动和响应分析。
学生需要运用所学的高等结构动力学理论，计算结构的固有频率、模态形态等。
2.实验模拟：作业可以设计实验模拟任务，要求学生使用相关软件或设备进行结构的振动试验，获取结构的模态参数和响应曲线数据。
3.结构优化：作业可以要求学生对给定的结构进行优化设计，以降低结构的振动响应或改善结构的抗震性能。
学生需要提出具体的优化方案，并进行相应的计算和分析。
4.报告撰写：作业最终要求学生将研究成果整理成报告。
报告应包括问题陈述、理论分析或实验过程、计算方法和结果分析等内容，以及对结论和进一步研究的讨论。
三、评价方式1.报告评估：根据学生的报告内容、结构分析和计算准确性、结果分析等方面，评估学生对高等结构动力学的理解和应用能力。
可以采用定量评价指标和评分标准进行评估。
2.讨论与答辩：在评价阶段，可以组织学生进行讨论和答辩，让学生互相交流和分享研究成果，进一步加深对问题的理解和探讨。
3.同伴评价：可以引入同伴评价的方式，让学生互相评价和给出建议，促进学生之间的交流和学习。
4.教师评价：教师对学生的报告进行评价，包括对报告内容、分析思路和计算方法的评估，提供及时的反馈和指导。
高等结构动力学大作业
高安槽钢是一种优质的建筑钢材，广泛应用于建筑、桥梁、输电塔等领域。
它以其优良的性能和合理的价格赢得了广大用户的好评。
一、高安槽钢简介
高安槽钢，全称高安热轧槽钢，是一种热轧成型的槽形钢材。
它通常由碳素结构钢、优质碳素结构钢、低合金结构钢等材料制成，具有较好的力学性能和耐腐蚀性能。
二、高安槽钢的参数
高安槽钢的参数主要包括：材质、规格、形状、尺寸等。
其中，材质决定了槽钢的力学性能和耐腐蚀性能；规格和形状则决定了槽钢在使用过程中的具体用途；尺寸则影响了槽钢的承载能力和使用寿命。
三、高安槽钢的应用领域
高安槽钢广泛应用于建筑、桥梁、输电塔、石油、化工、船舶、机车等领
域。
例如，在建筑领域，高安槽钢可用于搭建建筑框架、支撑结构等；在桥梁领域，高安槽钢可用于桥梁的支撑结构、加固结构等。
四、高安槽钢的优势与特点
高安槽钢具有以下优势和特点：
1.良好的力学性能：高安槽钢具有较高的抗拉强度、屈服强度和耐压强度，能够满足各种工程结构的使用要求。
2.耐腐蚀性能好：高安槽钢采用优质钢材制成，具有良好的耐腐蚀性能，可适用于各种环境。
3.尺寸精度高：高安槽钢采用先进的生产工艺，保证了产品的尺寸精度，便于施工安装。
4.质量稳定：高安槽钢的生产过程严格控制，保证了产品质量的稳定。
五、高安槽钢的生产厂家及联系方式
高安槽钢的生产厂家众多，其中以我国大型钢铁企业为主。
结构动力学试题及答案（本文按试题和答案格式进行编写）试题一：1. 请问什么是结构动力学？2. 简述结构动力学的研究对象和主要内容。
3. 结构动力学分析常用的方法有哪些？4. 结构动力学分析中常用的数学模型有哪些？5. 结构动力学的应用领域有哪些？答案一：1. 结构动力学是研究结构在外力作用下的动态响应及其稳定性的学科。
2. 结构动力学的研究对象是各种工程结构，主要内容包括结构的振动、冲击响应、瞬态响应和稳态响应等。
3. 结构动力学分析常用的方法有模态分析法、频率响应分析法、时程分析法等。
4. 结构动力学分析中常用的数学模型有单自由度体系、多自由度体系、连续体系等。
5. 结构动力学的应用领域广泛，包括建筑结构工程、桥梁工程、风力发电机组、地震工程等。
试题二：1. 结构动力学分析中，模态分析的基本原理是什么？2. 简述模态分析的步骤和计算方法。
3. 常用的模态分析软件有哪些？4. 请问什么是结构的固有频率和阻尼比？5. 结构的模态振型对结构动力响应有什么影响？答案二：1. 模态分析是基于结构的振动特性，通过求解结构的固有频率、模态振型和阻尼比等参数，来研究结构的动力响应。
2. 模态分析的步骤包括建立结构有限元模型、求解结构的固有频率和模态振型、计算结构的阻尼比等。
常用的计算方法有有限元法、拉普拉斯变换法等。
3. 常用的模态分析软件有ANSYS、ABAQUS、MSC.NASTRAN等。
4. 结构的固有频率是结构在无外力作用下自由振动的频率，阻尼比是结构振动过程中能量耗散的程度。
5. 结构的模态振型对结构动力响应有很大影响，不同的模态振型会导致不同的振动特性和反应。
试题三：1. 结构动力学分析中，频率响应分析的基本原理是什么？2. 简述频率响应分析的步骤和计算方法。
3. 频率响应分析和模态分析有什么区别？4. 结构的频率响应函数和传递函数有什么区别？5. 频率响应分析在结构设计中的应用有哪些？答案三：1. 频率响应分析是研究结构在单频激励下的响应特性，通过求解结构的频率响应函数，来获得结构的响应。
结构动力学大作业------------------------------------------作者xxxx------------------------------------------日期xxxx结构动力学大作业班级土木卓越１２0１班学号U20１２10323姓名陈祥磊指导老师叶昆20１４。
１2.30 结构动力学大作业-—SDO Ｆ体系在任意荷载作用下的动力响应 一、结构参数计算结构为右图所示的 １、kg m 3101000⨯=m N k /1020006⨯= ２、m m m m N =⋅⋅⋅⋅⋅⋅==21 k k k k N λ==⋅⋅⋅⋅⋅⋅==213、结构参数中5=N ;0.1=λ。
二、确定各阶频率和振型多自由度体系自由振动时的运动方程为012121111=+⋅⋅⋅+++n n y k y k y k y m 022221212=+⋅⋅⋅+++n n y k y k y k ym .。
．.。
.12jN-1N02211=+⋅⋅⋅+++n nn n n n y k y k y k y m 写成矩阵形式即为[]{}[]{}{}0=+y K yM 假设此方程的解答为{}{}()αω+=t Y y sin ,带入到运动方程中得到振动方程[][](){}{}02=-Y M K ω此方程要有非零解必须满足频率方程[][]02=-M K ω,可解得各阶主频率i ω再根据 [][](){}(){}02=-i i Y M K ω可求出结构的主振型。
在主振型中,通常将最后一个位移值设定为1,只要在程序中加入下列语句：MDOF ．YMa ｔrix(：，ｉ)=ＭDO Ｆ.YMat ｒix(：,i ）/ＭDOF 。
YMatr ｉx(MD ＯF 。
ND,i)运行程序之后得到如下结果： 1、各阶频率i ω和周期i TW1 12.7２9０2６1 T1 ０。
49３６10８43W ２ 37.1558483２T ２ 0。
Advanced Structural Dynamics ProjectThe dynamic response and stability analysis of the beam under vertical excitationInstructor：Dr. Li WeiName：Student ID：1.Problem description and the purpose of the project1.1 calculation modelAn Eular beam subjected to an axial force. Please build the differential equation of motion and use a proper difference method to solve this differential equation. Study the dynamic stability of the beam related to the frequency and amplitude of the force. As shown in the Fig 1.1.Fig1.11.2 purpose and process arrangementa.learning how to create mathematical model of the continuoussystem and select proper calculation method to solve it.b.learning how to build beam vibration equation and solve Mathieuequation.ing Floquet theory to judge vibration system’s stability andanalyze the relationship among the frequency and amplitude of the force and dynamic response.This project will introduce the establishment of the mathematical model of the continuous system in section 2, the movement equation and the numerical solution of using MATLAB in section 3, Applying Floquent theory to study the dynamic stability of the beam related to the frequency and amplitude of the force in section 4. In the last of the project, we get some conclusions in section 5.2. the mathematical model of the systemThe geometric model of the beam and force-simplified diagram is shown in Fig.2.1.We assume that its stiffness(EI) is constant and the deflection of the beam is small, and the boundary conditions is simply support. Now the beam subjected to an axial force. We assume the force is equal to 0cos P t ω. F=f 0coswt yxFig.2.1We select the length of x ∆ in any position of the beam, the free-body diagram is shown in Fig.2.2.Fig.2.2Using equations of movement equilibrium, that is to say:+↑()y y F m a ∑=∆ （1）0G M +∑= （2） From equation (1), we will get:22),(),(ty x A t x x S t x S ∂∂∆=∆+-ρ （3）Divide equation (3) by x ∆ and take the limit:22xy A x S ∂∂=∂∂-ρ （4） Then synthesize equation (2),we can get:0),()],(),([),(),(=∆∆+--∆++-∆+x t x x S t x y t x x y F t x M t x x M （5）Divide equation (5) by x ∆ and take the limit:S xy F x M =∂∂+∂∂ （6） Combine equation (4) with (6):0222222=∂∂+∂∂+∂∂ty A x y F x M ρ （7） And 22),(xy EI t x M ∂∂= （8） Combine equation (7) with (8):0)(22222222=∂∂+∂∂+∂∂∂∂ty A x y F x y EI x ρ （9） We know EI is a constant, so0)(222244=∂∂+∂∂+∂∂ty m x y t F x y EI （10） In equation (10), m is the mass of unit length. Now we will use assumed-modes method. Named lx n t T t x y n πsin )(),(=,so: 0sin )(22244422=⎥⎦⎤⎢⎣⎡-+l x n T l n t F T l n EI dt T d m n n n πππ （11） 0))(1(222=-+n non n T F t F P dt T d n=1,2,...... （12） In the equation (12)222222,l EI n F m EI l n P n on ππ==And t F t F ϖcos )(= ,so0)cos 1(222=-+n non n T t F F P dt T d ϖ n=1,2,...... （13） 0)cos (22=-+n n T t dtT d ϖεδ （14） In the equation (12)4)(L n A EI πρδ= 22)(Ln F πε= Equation (14) is the Mathieu equation. it is difficult to solve the analytical solution directly, thus, we use the approximate derivative namely an average acceleration method to get the numerical solution from the reference.3. Numerical solution3.1 using MATLAB to solve equationWe will use the Newmark-β method [1] to solve equation (14). We can use the initial condition 00u u 和to integrate the move equation:m 0u cu ku ++= （15）Fig.3.1As shown in Fig.3.1))(2(11+++∆+=i i i i i u u t u u （16） )(4121+++∆+∆+=i i i i i i i u u t t u u u （17） ()0cos =-+i i u wt u εδ （18）From equation (16), (17) and (18), we will get:()i i u t wn u ∆∆--=∆cos εδ （19）i i ii u u t u 2)2(-∆∆=∆ （20） ()()22]cos [44]cos [)(2tt wn t u t wn u t u i i i ∆∆-+∆+∆-∆-=∆εδεδ （21） When applying the MATLAB, we need discrete the processing time t, get time step 02.0=∆t .When solving the vibration stability interval, there are three variables to participate in the discussion, namely w c ,,δ. So take a particular w first and discuss the remaining two parameters.From Floquent theory [2]，we can use parameter A to judge stability.Equation 0)()(22=++y t dt dy t dt y d ξξ（22） Take two sets of special solution:1)0(,0)0(0)0(,1)0(2211====yy yy （20） Parameter [2] )]()([2121T y T y A +=（21） If abs (A) is less than 1, the system is stability. And if abs (A) is greater than 1, the system is instability. When abs(A) is equal to 1, the system is critical state.We use MATLAB Codes to solve equations. We use ω=2 Math ieu Equation to judge the validity of the codes. From Fig.3.2 and Fig.3.3, wecan consider the codes are correct.In these follow figures, ω=2, the horizontal axis is δ, vertical axis is ε.Fig.3.2. stable domain in reference [3]and[2]Fig.3.3. stable domain in MATLAB solutionCompared Fig.3.2 with Fig.3.3, we can see that the stability domain of numerical solutions applying average acceleration method are consistent with the standard solutions. it can concluded that when the system have solution whose cycle is equal to π or 2π, )3,2,1(2 ==n n δThis chapter discusses the accuracy of the vibration stability determination with Floquent theory. The next chapter will discuss the numerical solution and stable domain and two parameters ’ influences on the stability for this question.4. Parameters influenceIn this part, we only consider two parameters, namely the frequency and amplitude of the force.4.1the influence of the force’s frequency4.1.1the stability of the systemWhen we discuss the stability of the system related to the frequency of the force, we should select some different frequencies, so we choose ω=1,2,4,6,8 and10. Using MATLAB codes, we can obtain the figs of the stability. We can know the stable region is bigger with the increase of the frequency in Fig.4.1.ω=1 ω=2ω=4 ω=6ω=8 ω=10Fig.4.1 stable domain with different ω4.1.2the response of the systemWhen we discuss the response of the system, the system should be stable. So we choose 7δ=，1ε=-，ω=2,4 and 6. In Fig.4.2, the cycle of the response increase and the range of the reactive amplitude is smaller with the increase of the frequency.Fig.4.2 responses of the system with different ω4.2the influence of the force’s amplitudeThe εis related to the force’s amplitude P. The cycle of the response a little increase and the range of the reactive amplitude is bigger with the increase of the force’s amplitude, in Fig.4.3.and Fig4.4.Fig.4.3 vibration response curve with different δThe red curve is w=2, δ=12,ε=1; the blue curve is w=2, δ=14, ε=1.This figure state that the vibration cycle is smaller and the amplitude have a little change with the increase of the δ.Fig.4.3 vibration response curve with different εThe red curve is w=2, δ=12,ε=1; the blue curve is w=2, δ=12, ε=5. This figure state that the amplitude is smaller and the vibration cycle have a little change with the increase of the ε.5. Conclusion(1) With the increase of the frequency, the stable region and the cycle of the response are bigger, but the range of the reactive amplitude is smaller.(2)With the incr ease of the force’s amplitude, the cycle of the response a little increase and the range of the reactive amplitude is much bigger. (3)Vibration in the stable region, the vibration cycle is smaller with the increase of the δ; the amplitude is smaller with the increase of the ε.AcknowledgementsThe author is grateful for upperclassman Li Yong, he give me much assistance. And the author is also grateful for Doctor Li Wei. In his classes, I felt very happy and can understand his class effectively. At last, the author is also grateful for classmates in the same laboratory, they give me much guidance and encourage. I gain a lot of knowledge through this study and I will work harder in the future.References1、ROY R.CRSIG, Jr. STRUCTURAL DYNAMICS. New York: John Wiley & Sons.2、王海期. 非线性振动. 北京: 高等教育出版社, 19923、顾志平. 非线性振动. 北京: 中国电力出版社, 2012。
高等结构动力学大作业在高等结构动力学课程的学习过程中，我们将接触到许多有关结构动力学的理论和方法。
本文将围绕高等结构动力学的内容，探讨其在工程实践中的应用和未来的发展趋势。
一、结构动力学简介结构动力学是研究结构在受到外界力作用下的响应和振动特性的学科。
它广泛应用于桥梁、建筑物、飞机、船舶等工程结构的设计和分析过程中。
在实际工程中，结构动力学的研究对于保证结构的安全性、提高结构的抗震性能至关重要。
二、结构动力学的应用领域1. 桥梁工程：结构动力学在桥梁工程中有着广泛的应用。
通过结构动力学分析，可以评估桥梁的振动响应，预测桥梁的疲劳寿命，并优化桥梁的设计参数，提高桥梁的安全性和使用寿命。
2. 建筑物工程：结构动力学在建筑物工程中也起到关键的作用。
通过结构动力学分析，可以评估建筑物在风荷载和地震荷载下的响应，为建筑物的设计提供科学依据，确保建筑物具备足够的抗震性能和舒适性。
3. 航空航天工程：在航空航天工程中，结构的振动特性和动态响应对于飞行安全至关重要。
结构动力学可以用于评估飞行器的疲劳寿命、优化飞行器的设计，提高飞行器的结构强度和稳定性。
三、结构动力学的方法和技术1. 动力学数学模型：结构动力学利用数学模型描述结构在受力作用下的运动规律。
常见的数学模型包括单自由度振动系统、多自由度振动系统以及连续体振动系统等。
2. 振动试验技术：振动试验技术是结构动力学研究中常用的方法之一。
通过振动试验可以获取结构的振动特性和模态参数，为结构分析和设计提供实验数据支持。
3. 数值计算方法：结构动力学的研究也离不开数值计算方法的支持。
常用的数值计算方法包括有限元法、边界元法、模态超级元法等。
这些方法可以用于求解结构的静力响应和动力响应，预测结构的疲劳寿命和抗震性能等。
四、结构动力学的挑战与前景1. 疲劳寿命与保养：在长期使用过程中，结构的疲劳寿命是一个需要关注的问题。
结构动力学可以通过疲劳寿命评估和振动监测技术帮助我们预测结构的损伤情况，以及制定合理的结构维修和保养策略。
结构动力学作业姓名：学号：目录1.力插值法 (1)1.1分段常数插值法 (1)1.2分段线性插值法 (4)2.加速度插值法 (7)2.1常加速度法 (7)2.2线加速度法 (9)附录 (12)分段常数插值法源程序 (12)分段线性插值法源程序 (12)常加速度法源程序 (13)线加速度法源程序 (13)1.力插值法力插值法对结构的外荷载进行插值，分为分段常数插值法和分段线性插值法，这两种方法均适用于线性结构的动力反应计算。
1.1分段常数插值法图1-1为一个单自由度无阻尼系统，结构的刚度为k ，质量为m ，位移为y (t )，施加的外力为P (t )。
图1-2为矩形脉冲荷载的示意图，图中t d 表示作用的时间，P 0表示脉冲荷载的大小。
图1-1 单自由度无阻尼系统示意图图1-2 矩形脉冲荷载示意图对于一个满足静止初始条件的无阻尼单自由度体系来说，当施加一个t d 时间的矩形脉冲荷载，此时结构在t d 时间内的位移反应可以用杜哈梅积分得到：0()sin ()2 (1cos )(1cos ) (0)tst st d P y t t d m ty t y t t Tωττωπω=-=-=-≤≤⎰(1-1)如果结构本身有初始的位移和速度，那么叠加上结构自由振动的部分，结构的位移反应为：02()cos sin (1cos) (0)st d y ty t y t t y t t Tπωωω=++-≤≤ (1-2)图1-3 分段常数插值法微段示意图对于施加于结构任意大小的力，将其划分为Δt 的微段，每一段的荷载都为一个常数(每段相当于一个矩形的脉冲荷载)，如图1-3所示，则将每一段的位移和速度写成增量的形式为：1cos t sin t (1cos t)iii i y P y y kωωωω+=∆+∆+-∆ (1-3)i+1/sin t cos t sin t iii y P y y kωωωωω=-∆+∆+∆ (1-4)程序流程图如下i+1cos t sin t (1cos t)iii y P y y kωωωω=∆+∆+-∆i+1/sin t cos t sin ti i i y Py y kωωωωω=-∆+∆+∆图1-4 分段常数插值法流程图根据流程图可以编写相应的算法，利用MATLAB 进行编程，程序源代码见附录。
为了验证程序的正确性，本文选取的以下的例题进行验证。
对于一个单自由度的无阻尼结构，当其受到一个周期荷载时，其结构响应分为稳态解和瞬态解，由于没有阻尼的影响，其瞬态解并不会衰减，其理论表达式为：021()()(sin sin )1p x t t t k ωβωβ=-- (1-5)式中，()x t 为位移响应，0p 为激励，k 为刚度，β为荷载频率与固有振动频率之比，ω为荷载频率，ω为结构固有频率。
现令0p 为1，k 为1，则ω为1，ω取为2/3。
程序求得的解与解析解对比如图1-5所示（由于理论解与程序基本重合，所以将理论解乘以-1，方便比较）：位移y时间ta )位移速度v时间tb )速度图1-5 分段常数插值法结果验证由图1-5可知理论解与程序算得的解基本重合，可以验证程序的准确性。
1.2分段线性插值法与分段常数插值法不同，分段线性插值法将每一微段的力当成一个线性的直线，对于每一个微段，可看成一个矩形和一个三角形脉冲的叠加。
图1-6为分段线性插值微段示意图。
图1-6 分段线性插值法微段示意图对于无阻尼的体系，后一个时间步的位移和速度可由前一个时间步相应的值求得：11cos sin (1cos )(1sin )ii i i i y P P y y t t t t k k tωωωωωω+∆=∆+∆+-∆+-∆∆(1-6) 11/sin cos sin (1cos )i i i i i y P P y y t t t t k k tωωωωωωω+∆=-∆+∆+∆+-∆∆(1-7) 分段线性插值法的流程图如图1-7所示，与分段常数插值法仅仅是迭代的方式有所不一样。
11cos sin (1cos )(1sin )ii i i i y P P y y t t t t k k tωωωωωω+∆=∆+∆+-∆+-∆∆11/sin cos sin (1cos )i i i i i y P P y y t t t t k k tωωωωωωω+∆=-∆+∆+∆+-∆∆图1-7 分段线性插值法流程图程序源代码见附录，同样利用1.1节的算例进行验证，所得的结果如图1-8所示。
位移y时间ta )位移速度v时间tb )速度图1-8 分段线性插值法结果验证由上图可知程序的正确性。
2.加速度插值法加速度插值法也叫逐步积分法，其对加速度进行插值，可分为常加速度法和线加速度法。
2.1常加速度法图2-1常加速度法微段示意图对于一个单自由度结构，其运动方程为：my()y()()()t c t ky t p t ++=(2-1)将式（1-1）转变为增量方程：m y c y k y p ∆+∆+∆=∆(2-2)在通过逐步积分，将时间转化为一系列微小的时间段t ∆ ，如图2-1所示，现令111()()2i i i i y t y y t t t ++=+<≤，则t 时间的速度可表示为：111()(()())21()()(()())2i i t ti i t t i i i y t dt y t y t dty t y t y t y t t++=+-=+∆⎰⎰令t =t i +1,则i +1时刻速度可以表示为：111()2ii i i y y y y t ++=++∆(2-3)同理，位移可以表示为：2111()4i i i i i y y y t y y t ++=+∆++∆(2-4)将式(2-3)、(2-4)代入式(2-1)，即：1111my i i i i cy ky p ++++++=(2-5)此时，式(2-5)中只有1i y +为未知变量，可直接求出1i y +，之后再利用式(2-3)、(2-4)，可求出t i +1时刻的速度与位移。
算法的流程图如下所示：将荷载作用的时间划分为0243y t =∆t t∆∆*4(2)2i i i i mP P c y my t∆=∆+++∆22i i iy y y t∆=∆-∆24()2i i i i y y y t y t ∆=∆-∆-∆图2-2 常加速度法流程图算例验证的结果如下图所示，说明了该程序的正确性。
由于需要对加速度进行插值，此处增加了加速度验证。
位移y时间ta )位移速度v时间tb )速度加速度a时间tc )加速度图2-3 常加速法结果验证2.2线加速度法线加速度法与常加速度法原理类似，其速度与位移的增量方程与常加速度法对应的方程略有不同，图2-4为线加速度法微段示意图。
a )位移b )速度图2-4线加速度法微段示意图关于具体原理不在赘述，下图为线加速度法的流程图。
将荷载作用的时间划分为n 个分段，并计算263()()()()k t k t m c t t t=++∆∆6()()[()3()]()[3()()]2tP t P t m y t y t c t y t y t t ∆∆=∆++++∆/y P k∆=∆266()()()3()()y t y t y t y t tt∆=∆--∆∆3()()3()()2ty t y t y t y t t ∆∆=∆--∆0002332()2y y y t t=-∆∆图2-5 线加速度法流程图算例验证的结果如图2-6所示。
位移y时间ta )位移 速度v时间tb )速度加速度a时间tc )加速度图2-6 线加速法结果验证附录分段常数插值法源程序function x=inter_force_constant(p,w,dt,k,v0,y0)%分段常系数插值法%p代表输入的力，w为结构基本频率，dt为时间间隔,k为结构刚度%v0为初始的速度，y0为初始的位移%输出的矩阵第一列代表时间，第二列代表位移，第三列代表速度[r,~]=size(p);x=NaN(r,3);x(:,1)=p(:,1);x(1,2)=y0;x(1,3)=v0;for i=1:r-1x(i+1,2)=x(i,2)*cos(w*dt)+x(i,3)/w*sin(w*dt)+p(i,2)/k*(1-cos(w*dt));x(i+1,3)=(-x(i,2)*sin(w*dt)+x(i,3)/w*cos(w*dt)+p(i,2)/k*sin(w*dt))*w;end分段线性插值法源程序function x=inter_force_line(p,w,dt,k,v0,y0)%分段线性插值法%p代表输入的力，w为结构基本频率，dt为时间间隔,k为结构刚度%v0为初始的速度，y0为初始的位移%输出的矩阵第一列代表时间，第二列代表位移，第三列代表速度[r,~]=size(p);x=NaN(r,3);x(:,1)=p(:,1);x(1,2)=y0;x(1,3)=v0;for i=1:r-1x(i+1,2)=x(i,2)*cos(w*dt)+x(i,3)/w*sin(w*dt)+p(i,2)/k*(1-cos(w*dt))+(p(i+1,2)-p(i,2))/k/w/dt*(w*dt-sin(w*dt));x(i+1,3)=(-x(i,2)*sin(w*dt)+x(i,3)/w*cos(w*dt)+p(i,2)/k*sin(w*dt)+(p(i+1,2)-p(i,2))/k/w/dt*(1-cos(w*dt)))*w;end常加速度法源程序function x=inter_a_constant(p,w,m,keci,dt,v0,y0,k)%p代表输入的荷载，w为结构基本频率，keci为阻尼比，dt为时间间隔,m为结构质量,k 为结构刚度%v0为初始的速度，y0为初始的位移%输出的矩阵第一列代表时间，第二列代表位移,第三列代表速度，第四列代表加速度[r,~]=size(p);x=NaN(r,4);x(:,1)=p(:,1);x(1,2)=y0;x(1,3)=v0;x(1,4)=4/3/dt/dt*y0;c=2*keci*w;for i=1:r-1K=k+2*c/dt+4*m/dt/dt;dP=p(i+1,2)-p(i,2)+(4*m/dt+2*c)*x(i,3)+2*m*x(i,4);dY=dP/K;dV=2/dt*dY-2*x(i,3);dA=4/dt/dt*(dY-x(i,3)*dt)-2*x(i,4);x(i+1,2)=x(i,2)+dY;x(i+1,3)=x(i,3)+dV;x(i+1,4)=x(i,4)+dA;end线加速度法源程序function x=inter_a_line(p,w,m,keci,dt,v0,y0,k)%p代表输入的荷载，w为结构基本频率，keci为阻尼比，dt为时间间隔,m为结构质量,k 为结构刚度%v0为初始的速度，y0为初始的位移%输出的矩阵第一列代表时间，第二列代表位移,第三列代表速度，第四列代表加速度[r,~]=size(p);x=NaN(r,4);x(:,1)=p(:,1);x(1,2)=y0;x(1,3)=v0;x(1,4)=3/2/dt/dt*y0-3/2/dt*v0;c=2*keci*w;for i=1:r-1K=k+3*c/dt+6*m/dt/dt;dP=p(i+1,2)-p(i,2)+m*(6/dt*x(i,3)+3*x(i,4))+c*(3*x(i,3)+dt/2*x(i,4));dY=dP/K;dV=3/dt*dY-3*x(i,3)-dt/2*x(i,4);dA=6/dt/dt*dY-6/dt*x(i,3)-3*x(i,4);x(i+1,2)=x(i,2)+dY;x(i+1,3)=x(i,3)+dV;x(i+1,4)=x(i,4)+dA;end。
