下载文档原格式
模态叠加法一.思想要点是在积分运动方程以前，利用系统自由振动的固有振型将方程组转换为n 个相互不耦合的方程，对这种方程可以解析或数值地进行积分。
对于每个方程可以采用各自不同的时间步长，即对于低阶振型可采用较大的时间步长。
当实际分析的时间历程较长，同时又只需要少数较低阶振型的结果时，采用振型叠加法将是十分有利的。
求解步骤：1.求解系统的固有频率和振型2.求解系统的动力响应二.求解固有频率与振型（求解不考虑阻尼影响的振动方程） ..()(){0}M a t Ka t += 解可假设为：0sin ()a t t φω=-φ是n 阶向量，ω是向量φ的振动频率，t 是时间变量，0t 是由初始条件确定的时间常数。
代入振动方程，得到一个广义特征值问题：20K M φωφ-=求解可得n 个特征解221122(,),(,),ωφωφ···2,(,)n n ωφ120ωω≤<<···n ω< 特征向量12,,φφ···,n φ代表系统的n 个固有振型，幅度可按以下要求规定T i i M φφ=1（i=1，2，···，n ），这样规定的固有振型又称正则振型。
将22(,)(,)i i j j ωφωφ代回特征方程，得：2i i i K M φωφ= 2j j j K M φωφ=前式两边前乘以j φT，后式两边前乘以i φT ，得：2j i i j i K M φφωφφTT = 2i j i i jK M φφωφφT T = 由()TTj i j i i j K K K φφφφφφT T==得：22i j i j i j M K ωφφωφφT T =，推出22()0i j j i M ωωφφT-=当i j ωω≠时，有0j i M φφT =这表明固有振型对于矩阵M 是正交的，可表示为：1 ()0 ()i j i j M i j φφT=⎧=⎨≠⎩得：2 ()0 ()i i j i j K i j ωφφT ⎧==⎨≠⎩如果定义123n [ ]φφφφΦ=K21222 0 0 n ωωω⎡⎤⎢⎥⎢⎥Ω=⎢⎥⎢⎥⎢⎥⎣⎦O则特征解的性质可表示成：M K T T ΦΦ=I ΦΦ=Ω原特征值问题可表示为：K M Φ=ΦΩ三.求解动力响应1.位移基向量的变换引入变换()()1ni i i a t x t x φ==Φ=∑其中()[]12 n x t x x x =L代入运动方程，并两边前乘以T Φ，可得：()()()()()...x t C x t x t Q t R t T T +ΦΦ+Ω=Φ= 初始条件相应地转换成：..0000 x x Ma M a T T =Φ=Φ 阻尼为振型阻尼，则：()()2 i=j 0 i j i i ij C ωξφφT ⎧⎪=⎨≠⎪⎩ 或11222 0 2 0 2n n C ωξωξωξT ⎡⎤⎢⎥⎢⎥ΦΦ=⎢⎥⎢⎥⎣⎦O 其中i ξ（i=1,2,···,n ）是第i 阶振型阻尼比，可得n 个相互不耦合的二阶常微分方程()()()()...22i i i i i i i x t x t x t r t ωξω++= （i=1，2，···，n ）若C 是Rayleigh 阻尼，即C M K αβ=+根据试验或相近似结构的资料已知两个振型的阻尼比i ξ和j ξ，可得22222()()2()()i j j i i j j i j j i i j i ξωξωαωωωωξωξωβωω-=--=-2.求解单自由度系统振动方程在振动分析中常常采用杜哈美（Duhamel ）积分，又称叠加积分，其基本思想是将任意激振力()i r t 分解为一系列微冲量的连续作用，分别求出系统对每个微冲量的响应，然后根据线性系统的叠加原理，将它们叠加起来，得到系统对任意激振的响应。
结构动力分析研究结构在动荷载作用的响应（如位移、应力、加速度等的时间历程），以确定结构的承载能力和动力特性等。
ANSYS动力分析方法有以下几种，现分别做简要介绍。
1.模态分析用模态分析可以确定设计中的结构或机器部件的振动特性（固有频率和振型）。
它也可以作为其他更详细的动力学分析的起点，例如瞬态动力学分析、谐响应分析、谱分析。
用模态分析可以确定一个结构的固有频率和振型。
固有频率和振型是承受动态荷载结构设计中的重要参数。
如果要进行谱分析或模态叠加法谐响应分析或瞬态动力学分析，固有频率和振型也是必要的。
ANSYS的模态分析是一线性分析，任何非线性特性（如塑性和接触单元）即使定义了也将忽略。
可进行有预应力模态分析、大变形静力分析后有预应力模态分析、循环对称结构的模态分析、有预应力的循环对称结构的模态分析、无阻尼和有阻尼结构的模态分析。
模态分析中模态的提取方法有七种，即分块兰索斯法、子空间迭代法、缩减法或凝聚法、PowerDynamics法、非对称法、阻尼法、QR阻尼法，缺省时采用分块兰索斯法。
2.谐响应分析任何持续的周期荷载将在结构中产生持续的周期响应（谐响应）。
谐响应分析使设计人员能预测结构的持续动力特性，从而使设计人员能够验证其设计能否成功地克服共振、疲劳及其他受迫振动引起的有害效果。
谐响应分析是用于确定线性结构在承受随时间按正弦（简谐）规律变化的荷载时的稳态响应的一种技术。
分析的目的是计算出结构在几种频率下的响应并得到一些响应值（通常是位移）对频率的曲线。
从这些曲线上可以找到“峰值”响应，并进一步观察频率对应的应力。
这种分析技术只计算结构的稳态受迫振动。
发生在激励开始时的瞬态振动不在谐响应分析中考虑。
谐响应分析是一种线性分析。
任何非线性特性，如塑性和接触（间隙）单元，即使被定义了也将被忽略，但在分析中可以包含非对称系统矩阵，如分析流体—结构相互作用问题。
谐响应分析同样也可以分析有预应力结构，如小提琴的弦（假定简谐应力比预加的拉伸应力小得多）。
高层建筑结构模态及固有频率计算与分析秦泗建【期刊名称】《《四川建材》》【年(卷),期】2019(045)011【总页数】2页(P60-61)【关键词】高层建筑; 固有频率; 三维模型; 应力【作者】秦泗建【作者单位】江苏省连云港市赣榆区柘汪镇建管所江苏连云港 222113; 江苏凯港集团有限公司江苏连云港 222113【正文语种】中文【中图分类】TU973.20 前言对于高层建筑结构而言，影响其质量的关键因素包括：强度性能、刚度性能以及稳定性等。
各设计要素决定着该建筑的使用寿命和安全性。
因此，在高层建筑的结构设计方面，如何提升建筑的安全系数，一直都是设计研发的关键。
一直以来，人们在如何确保建筑物的良好强度性能方面，做了大量的研究工作。
例如：文献[1]采用了三种混凝土强度检测方案，为建筑工程的材料选择和施工提供了可靠的质量保证；文献[2]从提升建筑结构稳定性方面，提出了关于建筑结构的力学分析法，达到了进一步提升建筑结构强度的目标；文献[3]针对钢结构建筑，在考虑震动因素前提下，提出了关于承载力分析的有限元法，结果的计算准确率和效率都有较大的提高。
上述研究表明，人们在建筑结构的强度方面，取得了丰富的成果，极大地提升了建筑的安全系数。
然而，当建筑物施工完毕并投入使用之后，实际承受的载荷情况是多样化的。
例如：当风作用在建筑物上时，建筑物承受的是风载荷；当地震来临时，建筑物承受对应的动载荷。
在建筑设计和选材的时候，由于材料的许用应力较大，通常情况下这些载荷不会超过建筑物本体的许用应力。
但是，这些载荷的动力源却存在共同的特点，即频率特征，无论是风力还是地震的作用，都会使高层建筑发生振动。
在这种条件下，如果建筑物自身的固有频率接近于外在动力源的固有频率值，则会产生共振的风险。
其结果是，建筑物本身产生破坏、断裂以及变形等，严重危害居住者的生命财产安全。
因此，就高层建筑结构而言，如何避免共振，是一个务必要研究的问题。
模态和固有频率关系关于模态和固有频率之间的关系，我们需要先了解模态和固有频率的概念以及它们在不同领域中的应用。
接下来，我们将逐步讨论模态和固有频率之间的关系，并详细解释它们在物理学、工程学和音乐中的应用。
下面是一步一步回答这个问题的具体步骤。
第一步：介绍模态和固有频率的定义模态是指物体振动、震动或挠曲时所处的特定状态或形态。
在物理学中，模态是指系统在特定激励下的运动方式。
在机械振动中，模态是指由自由度、质量和刚度决定的系统特征。
而固有频率是指物体在某一模态下振动的频率，也可以理解为系统在该模态下的固有振动频率。
第二步：阐述模态和固有频率之间的关系模态和固有频率之间存在着密切的关系。
一个系统可以具有多种不同的模态，每种模态都对应着不同的固有频率。
可以说，模态和固有频率是相互依存的关系。
不同的模态对应着不同的固有频率，而固有频率决定了系统的振动行为。
第三步：在物理学中的应用在物理学中，模态和固有频率的研究对于了解物体的振动特性、结构的稳定性和强度等方面至关重要。
例如，在工程结构设计中，通过分析结构的模态和固有频率，可以确定结构在不同频率下的振动模式，从而评估结构的稳定性和耐久性。
在微观领域中，模态和固有频率的研究也可以帮助研究人员了解物质分子的振动特性和结构稳定性。
第四步：在工程学中的应用在工程学中，模态和固有频率的研究被广泛应用于结构动力学、振动噪声控制、飞行器设计等领域。
通过模态分析和固有频率的计算，可以预测结构的动力特性、避免共振现象的发生、并优化设计方案。
在飞机和船舶等交通工具的设计中，模态和固有频率的分析可以帮助设计师避免共振现象，提高结构的稳定性和耐久性。
第五步：在音乐中的应用在音乐学中，模态和固有频率的概念与音乐的音高和音阶有关。
在西方音乐中，固有频率对应着音高的概念，不同音高的音符对应着不同的固有频率。
而在传统音乐中，模态则是指具有不同音阶形式的曲调。
例如，西方音乐中的大调和小调就是两种不同的模态。
模态叠加法计算地震加速度时程反应的几个问题潘旦光;李雪菊;芦盼【摘要】针对模态叠加法进行结构加速度时程反应分析,讨论了广义坐标积分步长、模态截断和离散位移时程求导三个问题对加速度时程反应误差的影响;并提出了振型加速度贡献系数的模态截断指标.以简谐荷载作用下单自由度体系的地震反应和三条地震波作用下的5层框架结构的地震反应为例,进行数值计算.算例分析结果表明:当广义单自由度计算的离散时间步长小于1/32的自振周期时,加速度反应的误差小于5% ;对于加速度时程反应而言,基于振型参与质量选取的模态数偏少,应基于振型加速度贡献系数作为模态截断的依据;由离散位移时程经中心差分法所得加速度的误差与直接由模态叠加法所得的基本相同,而快速傅里叶变换( FFT)方法所得加速度时程存在虚假的高频振动而误差较大,采用低通滤波可有效降低误差.【期刊名称】《科学技术与工程》【年(卷),期】2019(019)010【总页数】8页(P155-162)【关键词】模态叠加法;积分步长;模态截断;中心差分法;滤波【作者】潘旦光;李雪菊;芦盼【作者单位】北京科技大学土木工程系,北京100083;同济大学土木工程防灾国家重点实验室,上海200092;北京科技大学土木工程系,北京100083;北达科他州立大学物流运输系, Fargo 58108-6050【正文语种】中文【中图分类】TU311.3在地震等复杂荷载作用下结构时程反应计算方法很多[1—4]，其中，模态叠加法理论简单、计算效率高而成为最常用的方法。
对于线弹性体系，模态叠加法理论上严密而精确，但少数低阶模态动力反应的解作为结构动力反应解必然存在模态截断误差[5,6]。
为在进行动力反应分析前，预估模态截断所产生的误差，不同学者提出不同指标以选择合适的模态阶数。
振型参与质量[2,7]作为最常用的模态截断依据而被规范和教材所采用。
Wilson[2]在分析振型参与质量比的基础上进一步提出静荷载参与比和动荷载参与比指标，对于地震而言，动荷载参与比和振型参与质量本质上是相同的。
STAAD模态分析与固有频率求解方法概述模态是结构的固有振动特性，每一个模态具有特定的固有频率，阻尼比和模态振型，获取这些结构振动特性的过程称之为模态分析或频率振型分析。
结构分析中经常会用到结构的这些固有振动特性，比如底部剪力法求解地震作用时需要用到结构的基本自振周期，再比如说利用振型分解法求解多自由度体系的各种动力分析都需要用到结构的各阶周期和振型。
因此，模态分析不仅是求解结构振动特性的方法，也是动力分析的基础。
本文将模态分析的求解方法进行全面介绍。
STAAD提供了两种求解结构模态的方法，分别是瑞利法和特征向量法。
1.CALCULATE RAYLEIGH (FREQUENCY) 瑞利法2.. MODAL CALCULATION (REQUESTED) 特征向量法瑞利法一般来说，工程结构的基频或者前几阶固有频率比较重要，瑞利法就是一种计算结构基频的常用近似算法。
瑞利法又叫做能量法，其核心思想是基于边界条件假定一个基频振型函数，然后利用能量守恒原理（最大动能和最大势能相等），从而求出结构的第一阶固有频率。
工程中，常选用构件在重力荷载工况下的静力挠度曲线作为基频振型曲线，而这个挠度曲线越接近实际得出来的基频越准确，当不能判断第一振型样子的时候，需要设置多种工况（比方自重在三个方向的三种工况），在每个工况中使用该命令，频率低者更接近基频。
在STAAD.Pro中，命令CALCULATE RAYLEIGH (FREQUENCY) 用于计算此命令之前的荷载工况在相应于变形方向上的结构近似频率。
在一组荷载的作用下结构将产生位移，程序中将这一位移作为振型，并计算出与该振型所对应的结构自振频率。
因此，这一命令应紧跟于其所在荷载工况的后面给出。
示例：LOADING 1 DEAD LOADSELFWEIGHT Y 1CALCULATE RAYLEIGH FREQLOADING 1 DEAD LOADSELFWEIGHT Z 1CALCULATE RAYLEIGH FREQ在这个例子中，程序将会分别计算荷载工况1和2的瑞利频率，并输出该频率的数值（单位是周/秒）及沿着整体坐标方向的最大挠度数值和所在的节点号。
ansys模态叠加法
ANSYS模态叠加法是一种结构动力学分析方法，其基本原理是将结构的自由振动模态按照一定的比例相加，从而得到结构在外力作用下的响应。
该方法通常用于求解结构的自由振动响应、地震响应以及材料疲劳寿命等问题。
在ANSYS中，模态叠加法可通过建立有限元模型、求解结构的固有频率和振动模态、以及进行模态叠加计算等步骤实现。
具体而言，该方法包括以下步骤：
1. 建立有限元模型：将结构分割成若干个有限元，并对其进行网格剖分和材料属性定义。
2. 求解结构的固有频率和振动模态：在ANSYS中，利用求解器求解结构的固有频率和振动模态。
3. 进行模态叠加计算：将结构的不同振动模态按照一定的比例相加，得到结构在外力作用下的响应。
ANSYS模态叠加法具有计算精度高、计算速度快等优点，可以广泛应用于结构动力学分析和相关工程领域。
- 1 -。
ansys 模态叠加法频率位移曲线ANSYS是一款广泛应用的有限元分析软件，它可以用于结构动力学、流体力学、电磁场、声学等多种领域的仿真计算。
在结构动力学中，ANSYS可以进行模态分析、谐响应分析、随机响应分析、响应谱分析等多种类型的动力学分析。
本文主要介绍ANSYS中的谐响应分析，特别是模态叠加法的原理和方法，以及如何绘制频率位移曲线。
## 谐响应分析的原理谐响应分析是一种分析结构在正弦激励下的动态响应的方法，它可以用于评估结构的振动特性、应力分布、疲劳寿命等。
谐响应分析的基本假设是：- 结构的物理性质和几何形状不随时间变化；- 结构的外载荷是随时间正弦变化的，且具有相同的频率； - 结构的响应也是随时间正弦变化的，且具有相同的频率和相位角；- 结构的响应是稳态的，即不考虑初始条件和瞬态效应。
基于以上假设，结构的运动方程可以表示为：$$[M]\ddot{u}+[C]\dot{u}+[K]u=F\sin(\omega t+\phi)$$其中，$[M]$是结构的质量矩阵，$[C]$是结构的阻尼矩阵，$[K]$是结构的刚度矩阵，$u$是结构的位移向量，$F$是结构的外载荷幅值，$\omega$是结构的外载荷频率，$\phi$是结构的外载荷相位角。
由于结构的响应也是正弦变化的，可以假设：$$u=U\sin(\omega t+\theta)$$其中，$U$是结构的响应幅值，$\theta$是结构的响应相位角。
将上式代入运动方程，并利用三角函数的恒等式，可以得到： $$[-\omega^2[M]+i\omega[C]+[K]]U=F\cos(\phi-\theta)$$其中，$i$是虚数单位，满足$i^2=-1$。
上式是一个复数方程，可以分解为实部和虚部两个方程，分别表示结构的响应幅值和相位角与外载荷频率的关系，即：$$[-\omega^2[M]+[K]]U=F\cos(\phi-\theta)$$$$\omega[C]U=F\sin(\phi-\theta)$$上述两个方程可以用矩阵求解法或者迭代求解法求解，得到结构在不同频率下的响应幅值和相位角，进而可以得到结构的位移、速度、加速度、应力等响应。
探讨木结构古建筑结构模型固有频率的几种计算方法赵鸿铁;张风亮;薛建阳;隋龑;张锡成;马辉【期刊名称】《西安建筑科技大学学报（自然科学版）》【年(卷),期】2012(044)002【摘要】结构的固有频率是研究结构动力性能的最基本参数.采用激振锤人工敲击法和有限元软件模态分析法计算出木结构古建筑结构模型的固有频率,并将二者所得结果进行比较;同时,根据振动台试验,提出了一种计算木结构古建筑固有频率的新方法—振动台试验法,给出了计算简图并提出了此方法在实际试验中可能会遇到的几个问题.通过研究表明:人工敲击法测得振动台试验模型震前的固有频率为2.05 Hz,比地震最大加速度为100 gal震后的固有频率(2.1 Hz)略小；ansys有限元软件模态分析的结构一阶固有频率为2.226Hz,稍微高于激振锤人工敲击法测得结果；由于木结构古建筑构件之间的连续性和互相约束性能较差以及谱曲线放大以后,肉眼读取的数据结果存在一定的误差,可能会导致振动台试验法无法准确的计算出结构固有频率.为木结构古建筑的研究和修缮加固提供参考依据.【总页数】6页(P159-163,176)【作者】赵鸿铁;张风亮;薛建阳;隋龑;张锡成;马辉【作者单位】西安建筑科技大学土木工程学院,陕西西安710055;西部建筑科技国家重点实验室(筹),陕西西安710055;西安建筑科技大学土木工程学院,陕西西安710055;西安建筑科技大学土木工程学院,陕西西安710055;西部建筑科技国家重点实验室(筹),陕西西安710055;西安建筑科技大学土木工程学院,陕西西安710055;西安建筑科技大学土木工程学院,陕西西安710055;西安建筑科技大学土木工程学院,陕西西安710055【正文语种】中文【中图分类】TU366;TU311.1【相关文献】1.古建筑木结构柱础连接力学模型 [J], 王子昂;郭瑞;钱行健;杨睿昕;代健波;何佳勇2.木结构古建筑安全性评估方法的层次分析法模型 [J], 朱晓强3.CFRP材料在古建筑木结构保护领域的适用性探讨——以留园曲溪楼的加固修复为例 [J], 徐恺4.中国木结构古建筑抗震性能的探讨 [J], 卢蓬蓬5.古建筑木结构典型拼合构件抗弯性能计算方法研究 [J], 贾肖虎;淳庆;张承文因版权原因，仅展示原文概要，查看原文内容请购买。
固有频率的计算公式
1. 弹簧振子系统。
- 对于水平放置的弹簧振子（忽略摩擦力等阻力），其固有频率f的计算公式为f = (1)/(2π)√(frac{k){m}}，其中k是弹簧的劲度系数，单位是N/m；m是振子的质量，单位是kg。
- 推导过程：根据胡克定律F=-kx（F是弹簧的弹力，x是振子偏离平衡位置的位移），结合牛顿第二定律F = ma（a是加速度），可得ma=-kx，即a =-(k)/(m)x。
对于简谐振动，其加速度a = - ω^2x（ω是角频率），所以ω=√(frac{k){m}}。
又因为f=(ω)/(2π)，所以f = (1)/(2π)√(frac{k){m}}。
2. 单摆系统（小角度摆动情况，摆角θ<5^∘近似认为是简谐运动）
- 固有频率f=(1)/(2π)√(frac{g){l}}，其中g是重力加速度，g = 9.8m/s^2（在地球表面附近），l是单摆的摆长，单位是m。
- 推导过程：单摆的回复力F = -mgsinθ，当θ<5^∘时，sinθ≈θ（θ用弧度制表示），设摆长为l，x = lθ（x是偏离平衡位置的位移），则F=-(mg)/(l)x。
根据牛顿第二定律F = ma，可得a =-(g)/(l)x，对比简谐运动的加速度公式a = - ω^2x，可得
ω=√(frac{g){l}}，再由f=(ω)/(2π)，得到f=(1)/(2π)√(frac{g){l}}。
optistruct模态叠加法
OptiStruct是一种用于结构分析和优化的有限元分析软件。
在OptiStruct中，模态叠加法是一种用于计算结构的模态响应的方法。
模态叠加法基于模态分析的原理，通过计算结构的固有频率、振型和阻尼比来确定结构的模态响应。
在模态叠加法中，结构的响应可以表示为各个模态振型的线性组合。
模态叠加法的基本步骤包括：
1. 模态分析：使用OptiStruct进行模态分析，计算结构的固有频率、振型和阻尼比。
2. 模态叠加：将结构的响应表示为各个模态振型的线性组合。
通常，前几个固有频率较低的模态对结构的响应起主导作用，因此可以选择使用较少的模态进行叠加。
3. 叠加系数：确定每个模态振型的叠加系数。
可以使用初始条件、边界条件或外部激励来确定叠加系数。
4. 响应计算：根据叠加系数和模态振型计算结构的响应。
可以计算结构的位移、速度、加速度或其他感兴趣的响应。
5. 结果分析：分析计算得到的响应结果，评估结构的性能和安全性。
模态叠加法在结构分析和优化中广泛应用，可以用于研究结构的动
力响应、模态超验、共振、振动幅值等问题。
通过模态叠加法，可以更好地理解和预测结构在不同工况下的响应行为，从而优化结构设计和改进结构性能。
桥梁模态分析方法及应用1.桥梁结构建模：首先，需要将桥梁结构进行合理的简化和离散化处理，将其转化为一个由节点和单元组成的有限元模型。
节点表示结构中的关键位置，而单元则表示结构中的连接部分。
同时，还需考虑结构材料的物理性质和边界条件。
2.模态分析求解：基于桥梁结构的有限元模型，采用模态分析方法，计算出结构的振动特性。
常用的求解方法包括传统的模态超级位置法和模态叠加法，以及现代的模态综合法和模态相对位移法等。
3.模态参数提取：通过模态分析求解，可以得到各个模态的频率、振型和阻尼比等参数。
频率表示结构振动的频率，振型表示结构振动的模态形态，阻尼比表示结构振动的耗能能力。
4.模态分析结果分析：根据模态分析提取出的模态参数，可以对桥梁结构的振动特性进行分析和评估。
例如，可以分析结构的固有频率范围，评估结构的稳定性；可以分析结构的振型形态，对结构的设计进行优化。
1.结构稳定性评估：通过模态分析，可以对桥梁结构的固有频率进行计算和分析。
当结构的固有频率接近外界激励频率时，会产生共振现象，从而对结构的稳定性造成威胁。
通过模态分析，可以评估结构的固有频率范围，及时发现潜在的共振问题。
2.结构安全性评估：桥梁结构在交通载荷和地震荷载等外部力的作用下，会发生振动现象。
模态分析可以计算得到结构的固有频率和振型，进而对结构在不同工况下的振动反应进行分析和评估。
通过模态分析，可以确定结构的应力、挠度等响应情况，从而评估结构的安全性。
3.结构设计优化：模态分析可以提供结构的振动特性，进而可以对结构进行优化设计。
通过调整结构的材料、截面形状和布置，可以改变结构的固有频率和振型，从而达到减小振动响应和提高结构的动力性能的目的。
4.结构加固与改造：对于已有桥梁结构，在其设计寿命内，可能需要进行加固和改造工作。
模态分析可以帮助评估结构的性能和弱点，进而指导结构的加固和改造方案的确定。
通过改变结构的刚度和阻尼特性，可以减小结构的振动响应，提高结构的承载能力和耐久性。
模态叠加法算法理论及其编程实现模态叠加法（Modal Superposition Method）是一种广泛应用于结构动力学计算中的数值分析方法，用于求解结构物的自由振动和响应。
该方法基于弹性力学原理，将结构物的振动模态进行叠加求解，得到结构物的整体振动响应。
模态叠加法的理论基础是振动理论和线性时变系统的特性。
在模态叠加法中，首先需要进行模态分析，即求解结构物的固有振动模态。
固有振动模态是结构物在无外界扰动的情况下自发振动的模式，可以通过有限元方法等手段进行求解。
固有振动模态是结构物的基础振动形态，通过线性组合这些基础振动形态，可以得到任意时刻结构物的振动情况。
在模态叠加法中，结构物的振动可以表示为各个模态振动的叠加。
每个模态表示一个固有振动模态，由振形函数和频率确定。
假设有n个模态，则结构物的振动响应可以表示为：\[u(t)=\sum_{i=1}^{n} A_{i}\sin(\omega_{i}t+\phi_{i})\]其中，A_i为振幅，\omega_i为频率，\phi_i为初始相位。
模态叠加法的关键是确定各个模态的振幅和初始相位。
确定各个模态的振幅和初始相位可以通过结构物的初始条件和激励情况来确定。
当结构物受到初始条件的影响时，振动模态的振幅和初始相位可以由初始条件确定。
当结构物受到外界激励时，振动模态的振幅和初始相位可以由结构物的动态响应计算得到。
根据叠加原理，结构物的振动响应可以表示为各个模态响应的叠加。
通过求解每个模态的振动响应，再进行叠加，可以得到结构物的整体振动响应。
在进行模态叠加法的编程实现时，一般可以采用以下步骤：1.进行结构物的模态分析，求解固有振动模态。
2.根据激励情况和初始条件，确定各个模态的振幅和初始相位。
3.对每个模态进行振动响应分析，求解振动模态的振动响应。
4.将各个模态的振动响应进行叠加，得到结构物的整体振动响应。
在实际编程实现中，可以利用数值计算软件或编程语言来实现模态叠加法。
固有频率的计算公式固有频率是指系统在没有外界干扰下，根据其自身的特性和参数计算得到的频率。
它可以用于描述机械、电子、光学等不同领域中的系统振动频率。
在物理学中，一般使用弹性力学理论来计算固有频率。
弹性力学是研究物体在受力作用下发生形变和振动的力学学科，其中固有频率的计算是其中一个重要的问题。
以下是计算固有频率的一些常见公式：1.自由振动的单自由度系统：对于一个自由振动的单自由度系统，固有频率可以通过以下公式计算：f=1/(2π)*√(k/m)其中，f表示固有频率，k表示系统的弹簧常数，m表示系统的质量。
2.自由振动的多自由度系统：对于一个自由振动的多自由度系统，固有频率可以通过以下公式计算：ω^2*x=K*x其中，ω表示固有角频率，x表示系统的位移矢量，K表示系统的刚度矩阵。
上述方程可以通过对系统动力学方程进行求解，得到系统的固有角频率和振型。
3.机械振动系统：对于机械振动系统，可以使用以下公式计算固有频率：f=1/(2π)*√(K/m)其中，f表示固有频率，K表示系统的刚度，m表示系统的质量。
在机械工程中，刚度可以通过计算杆件的刚度、弹簧的刚度、轮毂刚度等来确定。
4.电磁振动系统：电磁振动系统的固有频率可以通过以下公式计算：f=1/(2π)*√(1/LC)其中，f表示固有频率，L表示电感，C表示电容。
该公式适用于电路中的振荡电路，如LC振荡电路和LCR振荡电路。
除了上述公式，还有许多其他的计算固有频率的公式，具体的计算方法取决于系统的特性和参数。
需要根据具体情况选择合适的公式进行计算。
此外，在实际应用中，还可以通过实验的方法来测量固有频率，以获得更准确的结果。
用模态叠加法求固有频率一、 模态分析法（振型叠加法）原理对于n 个自由度系统，其在广义坐标系下的运动微分方程为[]{}[]{}{}()M x k x F t += （1-1）设在t=0时，有初始条件：{}{}(0)0x x = 和 {}{}(0)0x x =通过求解特征值问题，可得系统的固有频率和振型向量{},(1,2,,)u i n ni iω= {}}(1,2,,)u i n i i ϕ=以正则振型矩阵[]ϕ作为变换矩阵，令{}[]{}x z ϕ= （a ）代入方程（1-1），并前乘以正则振型矩阵的转置T ϕ⎡⎤⎣⎦，得[][][]{}[][][]{}[]{}()T T T M z k z F t ϕϕϕϕϕ+= （b ）∵ [][][][]T M I ϕϕ=[][][][]21222n T k n nn ωωϕϕω⎡⎤⎢⎥⎢⎥=Λ=⎢⎥⎢⎥⎢⎥⎣⎦令 {}[]{}()()T P t F t ϕ= ---- 是正则坐标系下的激励。
则方程（b ）为{}[]{}{}()z z P t +Λ= （c ）展开后，得2()11112()22222()z z P t n z z P t n z z P t n nn n nωωω⎧+=⎪⎪+=⎨⎪⎪+=⎩ （1-2）式中 {}{}()()(1,2,,)T P t F t i n i i ϕ== ，为对应第i 个正则坐标的激励。
对于方程（1-2）是一组n 个独立的方程，每个方程和单自由度系统的强迫振动相同，因此可按单自由度系统中的方法独立地求解每个方程。
则由杜哈美积分得方程（1-2）的通解()0()cos sin 01()sin ()1,2,,0z i z t z t t i i ni ni nit P t d i n ni i ni ωωωτωττω=+⎰+-=式中0z i 和 0z i 是第i 个正则坐标的初始位移和初始速度。
∵ {}[]{}x z ϕ=∴ {}[]{}00x z ϕ= （d ）和 {}[]{}00x z ϕ= （e ）用 [][]T M ϕ 前乘以式（d ）两端，得[][]{}[][][]{}00T T M x M z ϕϕϕ=∴ {}[][]{}00T z M x ϕ=同理，有 {}[][]{}00T z M x ϕ=写成分量形式{}[]{},(1,2,,)00T z M x i n i i ϕ=={}[]{},(1,2,,)00T z M x i n i i ϕ==最后，由方程（a ），将正则坐标的解{}z 变换到原广义坐标{}x ，就得到方程（1-1）的解。
（2）在许多工程问题中系统的自由度很多，要想求出系统的所有固有频率和振型向量，计算成本很大，有时甚至是不可能的。
由于激励的高频成分很微弱，或者由于系统的高频振动没有激发出来，总之系统的响应中只有较低的几阶振型分量。
因此，使用振型叠加法可以使计算大大地简化。
例如，若系统为n 自由度，且只需考虑前p （p<<n ）阶固有特性，即只需计算出系统前p 阶固有频率和振型向量：{},(1,2,,)i p ni i ωϕ=则振型矩阵 []{}{}{}12p ϕϕϕϕ⎡⎤=⎢⎥⎣⎦是一个n ×p 的矩阵。
作如下的线性变换 {}[]{}11x z n p n p ϕ=⨯⨯⨯代入方程（1-1），并前乘以正则振型矩阵的转置[]T ϕ，得[][][]{}[][][]{}[]{}()T T T M z k z F t ϕϕϕϕϕ+= （b ）[][][][]T M I p n n n n p p p ϕϕ=⨯⨯⨯⨯[][][][]21222n T n k p n n n n p p p pp p pωωϕϕω⎡⎤⎢⎥⎢⎥==Λ⎢⎥⨯⨯⨯⨯⎢⎥⎢⎥⎣⎦⨯ {}[]{}()()11T P t F t p n p n ϕ=⨯⨯⨯∴ {}[]{}{}()111z z P t p p p p p +Λ=⨯⨯⨯⨯展开后，得2()11112()22222()z z P t n z z P t n z z P t ppp p p ωωω⎧+=⎪⎪+=⎨⎪⎪+=⎩ 则原广义坐标下的响应为{}[]{}{}{}{}121112x z z z z p n p p n p ϕϕϕϕ==+++⨯⨯⨯这种处理方法称为振型截断法。
二、振型叠加法计算步骤1. 建立广义坐标系下的运动微分方程[]{}[]{}{}()M x k x F t +=给出广义坐标系下初始条件：{}0x 和 {}0x2. 计算固有频率和振型向量{},(1,2,,)u i n ni i ω=计算主质量： {}[]{}(1,2,,)T M u M u i n a i i i== 得正则振型向量： {}}(1,2,,)u i n i iϕ= 3. 初始条件变换{}[]{}(1,2,,)00T z M x i n i i ϕ== {}[]{}(1,2,,)00T z M x i n i i ϕ== 4. 激励变换{}{}()()(1,2,,)T P t F t i n i i ϕ==5. 计算正则坐标系下的响应正则坐标系下运动微分方程2()(1,2,,)z z P t i n i ni i i ω+==利用杜哈美积分写出正则坐标系下的响应0()cos sin 01()sin ()(1,2,,)0z i z t z t t i i ni ni nit P t d i n ni i ni ωωωτωττω=+⎰+-=6. 写出原广义坐标系下的响应{}[]{}{}{}{}1212x z z z z n n ϕϕϕϕ==+++三、源程序******************************* Program to generate the response of an N-DOF structure * to any arbitrary applied loading using the* MODE SUPERPOSITION METHOD******************************* Program Author : Sanjoy Chakraborty* Auburn University******************************program mode_suse msimslimplicit real *8 (a-h,o-z)real *8 mr,krcharacter *12 flin,foutdimension mr(10,1),cr(10,1),kr(10,1),phi(10,10),tt(20,10),1 ft(20,10),fr(6000,10),q(6000,10),q1(6000,10),q2(6000,10),2 qt(10,1),q1t(10,1),q2t(10,1),nfor(10),xmin(10,1),vmin(10,1)data ndim/10/******************************** Assign I/O Devices and read in system parameters*******************************write(*,50)50format(///,' Program for Eigensolution and Mode Superposition'1 ,//,' Program Author : Sanjoy Chakraborty, Auburn University'2 ,//////,' Enter Input File Name : ')read(*,'(a)')flinopen(1,file=flin,status='old')write(*,'(//a)')' Enter Output File Name : 'read(*,'(a)')foutopen(5,file=fout,status='unknown')call input1(mr,cr,kr,phi,tt,ft,nd,nfor,ndim,delt,ttot,xmin,vmin)open(2,file='x.out',status='unknown')open(3,file='x1.out',status='unknown')open(4,file='x2.out',status='unknown')******************************** Generate Modal Force Vectors at discrete time intervals* fr(i,j) contains the modal force vector Fr at time step 'i'* j=1,nd corresponds to the DOF at which fr is applied * {fr} = [Phi]'.{ft}* where, {ft} = physical force vector at time step 'i'* [Phi] = modal transformation matrix for given system* [Phi]' = transpose of [Phi]*******************************call force1(tt,ft,delt,phi,nd,fr,nfor,ndim,ttot)******************************** Generate MODAL time history for all DOF's* q(i,j), q1(i,j), q2(i,j) contain the modal* displacements, velocities and accelerations resp.* i --> corresponds to time step number* j --> corresponds to a particular DOF******************************** Direct Integration of the uncoupled equations of motion* are carried out using the Newmark-B method*******************************do 10 i=1,ndcall newb(i,mr(i,1),cr(i,1),kr(i,1),fr,q,q1,q2,delt,ndim,ttot,1 xmin(i,1),vmin(i,1))10continue******************************** Transform MODAL time histories back to PHYSICAL Coords.* {x(t)} = [Phi].{q(t)} ..etc.* Print Output*******************************ic=1do 20 t=0.,ttot,deltdo 30 i=1,ndqt(i,1)=q(ic,i)q1t(i,1)=q1(ic,i)q2t(i,1)=q2(ic,i)30continuecall matmul(phi,qt,nd,mr,ndim)write(2,43)t,(mr(nn,1),nn=1,nd)43format(f8.4,10e14.6)do 40 i=1,ndq(ic,i)=mr(i,1)40continuecall matmul(phi,q1t,nd,mr,ndim)write(3,43)t,(mr(nn,1),nn=1,nd)do 41 i=1,ndq1(ic,i)=mr(i,1)41continuecall matmul(phi,q2t,nd,mr,ndim)write(4,43)t,(mr(nn,1),nn=1,nd)do 42 i=1,ndq2(ic,i)=mr(i,1)42continueic=ic+120continue*************************call res_spec(ic-1,q,q1,q2,delt,ndim,nd)************************** Generate values of MAXIMUM RESPONSE* @ each DOF*************************stopend*************************subroutine input1(m,c,k,phi,t,f,nd,nfor,ndim,delt,ttot,xmin,1 vmin)************************** [m1] is the system mass matrix (input)* [k1] is the system stiffness matrix (input)* [m] contains the modal mass parameters (generated)* [c] contains the modal damping parameters (generated)* [k] contains the modal stiffness parameters (generated)* alpha, beta = parameters to generate proportional damping matrix * [c] = alpha[m] + beta[k]* [phi] contains the modal transformation matrix* {t} contains the times at which the load vector* {f} is specified* nd = the number of degrees of freedom* nsol = 1 for EIGENSOLUTION ONLY* nsol = 2 for EIGENSOLUTION + MODE SUPERPOSITION ANALYSIS* delt = step size* ttot = total time for which response will be evaluated************************** This program uses the IMSL Math Subroutine DGVCSP to evaluate* the eigenvalues and eigenvectors of the dynamic system*************************implicit real *8 (a-h,o-z)real *8 m,k,m1,k1external dgvcspdimension m(ndim,1),c(ndim,1),k(ndim,1),phi(ndim,ndim),t(20,10)dimension nfor(10),f(20,10),m1(10,10),c1(10,10)dimension k1(10,10),omega(10)dimension phit(10,10),eva(10),evec(10,10)dimension temp1(10,10),temp2(10,10),temp3(10,10),xin(10,1),1 vin(10,1),xmin(10,1),vmin(10,1)read(1,*)nsol,nddo 90005 i=1,ndread(1,*)(m1(i,j),j=1,nd)90005continuedo 90006 i=1,ndread(1,*)(k1(i,j),j=1,nd)90006continue********************do 825 i=1,nddo 825 j=1,ndtemp1(i,j)=m1(i,j)temp2(i,j)=k1(i,j)825continuec Call the IMSL routine DGVCSP to evaluate the eigenvalues and c eigenvectors of the system.call dgvcsp(nd,temp2,ndim,temp1,ndim,eva,evec,ndim)do 903 i=1,ndomega(i)=sqrt(eva(nd-i+1))903continuec write(*,*)(omega(i),i=1,nd)c pausec stopdo 904 i=1,nddo 904 j=1,ndphi(i,j)=evec(i,nd-j+1)904continue******************************* Normalize Eigenvectors wrt 1st row values******************************do 826 j=1,ndtemp=phi(1,j)do 826 i=1,ndphi(i,j)=phi(i,j)/temp826continue******************************* Print EIGEN OUTPUT* Terminate if NSOL=1******************************write(5,905)905format(/,' E I G E N S O L U T I O N',//, 1' Mode # Omega(rad/sec)',/,2' ----------------------------')do 906 i=1,ndwrite(5,907)i,omega(i)907format(i6,e22.8)906continuewrite(5,912)912format(//,' Mode Shapes',/)do 908 i=1,ndwrite(5,'(10E14.6)')(phi(i,j),j=1,nd)908continueif(nsol.eq.1)stop****************************************read(1,*)delt,ttotread(1,*)alpha,betacall trans(phi,phit,nd,ndim)call matmul1(phit,m1,temp1,nd,ndim)call matmul1(temp1,phi,temp2,nd,ndim)do 909 i=1,ndm(i,1)=temp2(i,i)909continuecall matmul1(phit,k1,temp1,nd,ndim)call matmul1(temp1,phi,temp2,nd,ndim)do 910 i=1,ndk(i,1)=temp2(i,i)910continuedo 911 i=1,ndzeta=0.5*(alpha/omega(i)+beta*omega(i))c(i,1)=2.0*m(i,1)*omega(i)*zeta911continuedo 921 i=1,nddo 921 j=1,ndc1(i,j)=alpha*m1(i,j)+beta*k1(i,j)921continuewrite(5,922)nd922format(//,' S Y S T E M P R O P E R T I E S',1 //,' Degrees of Freedom =',i3,//,' Mass Matrix',/)do 913 i=1,ndwrite(5,914)(m1(i,j),j=1,nd)914format(10e15.6)913continuewrite(5,915)915format(//,' Stiffness Matrix',/)do 916 i=1,ndwrite(5,914)(k1(i,j),j=1,nd)916continuewrite(5,917)917format(//,' Damping Matrix',/)do 918 i=1,ndwrite(5,914)(c1(i,j),j=1,nd)918continuewrite(5,923)(m(i,1),i=1,nd)923format(//,' Modal Masses',//,10e15.6)write(5,924)(k(i,1),i=1,nd)924format(//,' Modal Stiffnesses',//,10e15.6)write(5,925)(c(i,1),i=1,nd)925format(//,' Modal Damping',//,10e15.6)read(1,*)(xin(i,1),i=1,nd)read(1,*)(vin(i,1),i=1,nd)do 851 i=1,nddo 851 j=1,ndtemp1(i,j)=0.0temp2(i,j)=0.0temp3(i,j)=0.0851continuedo 852 i=1,ndtemp1(i,i)=1.0/m(i,1)852continuecall matmul1(temp1,phit,temp2,nd,ndim)call matmul1(temp2,m1,temp3,nd,ndim)call matmul(temp3,xin,nd,xmin,ndim)call matmul(temp3,vin,nd,vmin,ndim)********************************** Input of forcing Function*********************************do 90011 i=1,ndnfor(i)=2t(1,i)=0.0t(2,i)=ttot+100.f(1,i)=0.0f(2,i)=0.090011continue90014read(1,*,err=90012)i,nfor(i)do 90013 j=1,nfor(i)read(1,*)t(j,i),f(j,i)90013continuegoto 9001490012continuereturnend**************************subroutine force1(tt,f,dt,p,nd,fr,nfor,ndim,ttot)**************************implicit real *8 (a-h,o-z)dimension tt(20,10),f(20,10),fr(6000,ndim),p(ndim,ndim),1 f1(10,1),pt(10,10),nfor(ndim),ft(10,1)do 90030 i=1,nddo 90030 j=1,ndpt(i,j)=p(j,i)90030continueic=1do 90020 t=0.,ttot,dt******************************** Interpolate to obtain value of forcing function Ft* at time = t*******************************do 90022 nf=1,nddo 90024 n=1,nfor(nf)-1if(t.ge.tt(n,nf).and.t.le.tt(n+1,nf))thenft(nf,1)=f(n,nf)+(f(n+1,nf)-f(n,nf))*1 (t-tt(n,nf))/(tt(n+1,nf)-tt(n,nf))goto 90022endif90024continue90022continue******************************** Generate MODAL force vector at time = t*******************************call matmul(pt,ft,nd,f1,ndim)do 90026 ii=1,ndfr(ic,ii)=f1(ii,1)90026continueic=ic+190020continuereturnend*******************************subroutine newb(ii,m,c,k,f,x,x1,x2,del,ndim,ttot,xin,vin) ******************************** Direct Integration of the uncoupled equations of motion * using the NEWMARK-b method******************************implicit real *8 (a-h,o-z)real *8 m,k,ksdimension f(6000,ndim),x(6000,ndim),x1(6000,ndim),x2(6000,ndim)******************************* Assign Initial Conditions (@ t=0.0)******************************x(1,ii)=xinx1(1,ii)=vinx2(1,ii)=(f(1,ii)-c*x1(1,ii)-k*x(1,ii))/mks=k+2.0*c/del+4.0*m/del**2******************************* Obtain time history at discrete time intervals of del******************************ic=1do 101 t=0.,ttot,deldps=f(ic+1,ii)+m*(4.0*x(ic,ii)/del**2+4.0*x1(ic,ii)1 /del+x2(ic,ii))+c*(2.0*x(ic,ii)/del+x1(ic,ii))x(ic+1,ii)=dps/ksx2(ic+1,ii)=(4.0/del**2)*(x(ic+1,ii)-x(ic,ii)-1 del*x1(ic,ii))-x2(ic,ii)x1(ic+1,ii)=(2.0/del)*(x(ic+1,ii)-x(ic,ii))-x1(ic,ii)c write(5,105)ic,t,f(ic+1,ii),f(ic,ii),dps,dxc105 format('**',i6,e17.7,/,5x,2f15.9,/,5x,2e20.10)ic=ic+1101continuereturnend*****************************subroutine matmul(a,b,nd,c,ndim)******************************* Used to multiply two matrices* [C] = [A].{B}*****************************implicit real *8 (a-h,o-z)dimension a(ndim,ndim),b(ndim,1),c(ndim,1)do 9030 i=1,nddo 9030 j=1,1c(i,j)=0.0do 9030 ii=1,ndc(i,j)=c(i,j)+a(i,ii)*b(ii,j)9030continuereturnend*****************************subroutine res_spec(n,x,x1,x2,del,ndim,nd)****************************** Obtain the maximum value of the response* (displacement, velocity or acceleration)* at each DOF for the time span considered****************************implicit real *8 (a-h,o-z)dimension x(6000,ndim),x1(6000,ndim),x2(6000,ndim)dimension xl(10),x1l(10),x2l(10),t(10),t1(10),t2(10)do 9070 i=1,ndxl(i)=x(1,i)x1l(i)=x1(1,i)x2l(i)=x2(1,i)t(i)=0.t1(i)=0.t2(i)=0.9070continuedo 9060 i=2,ndo 9060 j=1,ndif(abs(x(i,j)).gt.abs(xl(j)))thenxl(j)=x(i,j)t(j)=float(i-1)*delendifif(abs(x1(i,j)).gt.abs(x1l(j)))thenx1l(j)=x1(i,j)t1(j)=float(i-1)*delendifif(abs(x2(i,j)).gt.abs(x2l(j)))thenx2l(j)=x2(i,j)t2(j)=float(i-1)*delendif9060continuewrite(5,9068)9068format(///,' R E S P O N S E S P E C T R A',///, 1' DOF Maxm. @ t Maxm. @ t ', 2' Maxm. @ t',/,3' X Vel. ',4' Accln.',/,5' ----------------------------------------------------', 6'----------------')do 9065 i=1,ndwrite(5,9066)i,xl(i),t(i),x1l(i),t1(i),x2l(i),t2(i)9066format(i5,e15.7,f7.3,e15.7,f7.3,e15.7,f7.3)9065continuereturnend*************************************subroutine matmul1(a,b,c,nd,ndim)*************************************implicit real*8 (a-h,o-z)dimension a(10,10),b(10,10),c(10,10)do 80010 i=1,nddo 80010 j=1,ndc(i,j)=0.0do 80010 k=1,ndc(i,j)=c(i,j)+a(i,k)*b(k,j)80010continuereturnend***************************************subroutine trans(p,pt,nd,ndim)***************************************implicit real*8 (a-h,o-z)dimension p(10,10),pt(10,10)do 80001 i=1,nddo 80001 j=1,ndpt(i,j)=p(j,i)80001continuereturnend***********************************************。
