工程有限元方法动力响应分析与自由度缩减

1、有限元方法与应用有限元方法与应用结构结构动力有限元分析动力有限元分析 张有为张有为 工程力学工程力学系系内容回顾内容回顾 动力学两类基本问题动力学两类基本问题动载荷（又称动力分析）动载荷（又称动力分析）固有特性分析固有特性分析响应分析响应分析固固有有频频率率振振型型位位移移响响应应速速度度响响应应加加速速度度响响应应动动应应变变动动应应力力固有特性：固有特性：是一组模态参数构成，它由结构本身（质量与刚度分布）决定，而与外是一组模态参数构成，它由结构本身（质量与刚度分布）决定，而与外部载荷无关，但决定了结构对动载荷的响应；部载荷无关，但决定了结构对动载荷的响应；响应分析：响应分析：是计算结构对给

2、定动载荷的各种响应特性。是计算结构对给定动载荷的各种响应特性。模态分析模态分析时程分析时程分析自由振动问题受迫振动问题本节要点本节要点 动力响应分析动力响应分析 动力响应分析动力响应分析概述概述 直接积分方法直接积分方法 Newmark算法算法 中心差分中心差分 Wilson-法法 平均值法平均值法 基于基于ANSYS的动力响应分析的动力响应分析 自由度自由度缩减方法缩减方法 振型分解法振型分解法 位移凝聚法位移凝聚法动力响应分析概述动力响应分析概述 0000mxcxkxf txxxv 动力响应分析动力响应分析21,242ccmkim 12(cossin)txeCtCt24;22cmkcmm

3、通解部分200;1;2kcmmk012(cossin)txeCtCt00000(cossin)tvxxetxt比例阻尼系数 mxcxkxf t ckxxxf tmm 2002xxxf t动力响应分析概述动力响应分析概述 0tpxh tfd脉冲响应函数 cosf tt特解部分杜阿梅尔积分 01sin000tettmh tt 000002200coscos1sin1cost costsinttpttxh tdtdmeftdmm201直接积分法直接积分法首先首先给定待求时间长度给定待求时间长度T， 在在其中取一系列离散点其中取一系列离散点（i=0,1,n），我们不去），我们不去求求 ，只求，只求 即

4、即可，即给出可，即给出待求响应在各离散时刻的近似值。两个离散时刻间隔称为待求响应在各离散时刻的近似值。两个离散时刻间隔称为步步长长，步长可以相同（等步长），也可以不同（变步长）。，步长可以相同（等步长），也可以不同（变步长）。 0000MxCxKxF txxxv0, tT( )ix t 基本思想基本思想( )x t直接积分法直接积分法 时域中心差分法时域中心差分法原理原理1121122kkkkkkkXXXXtXXXtkkkkMXCXKXF基本假设：基本假设：112221121122kkkkMC XFKMXMC Xttttt1kX的等效刚度阵仅由的等效刚度阵仅由M、C组成，求解时组成，求解时K阵

5、不用求逆，阵不用求逆，显式算法显式算法联立联立1111kkkkkkkkkXXXtXXXtXXXt直接积分法直接积分法 时域中心差分法时域中心差分法d原理原理1121122kkkkkkkXXXXtXXXtkkkkMXCXKXF基本假设：基本假设：112221121122kkkkMC XttttFMtKXMC X联立联立流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，并，并计算下列积分常数计算下列积分常数 计算计算 计算等效刚度阵计算等效刚度阵2. 对每一时间步计算对每一时间步计算 tk时刻的等效力时刻的等效力 tk+1时

6、刻的位移时刻的位移 tk时刻的速度和加速度时刻的速度和加速度 0000MXFCXKX0X0120322111,2,2aaaa atat10030XXXta X 1X01Ka MaC2011kkkkFFKa M Xa MaC X1kkKXF0111112kkkkkkkXaXXXXaXX0a0a1010211022XXXXtXXXt1a1a2a直接积分法直接积分法 时域中心差分法稳定性讨论时域中心差分法稳定性讨论考虑无阻尼单自由度自由振动问题考虑无阻尼单自由度自由振动问题采用中心差分离散格式得到采用中心差分离散格式得到假定解的形式为假定解的形式为将（将（3）代入（）代入（2）式得到特征方程）式得到

7、特征方程通过求解（通过求解（4）的特征方程得到）的特征方程得到为了使得解具有振荡特性且不发散，要求为了使得解具有振荡特性且不发散，要求又因为又因为 同时同时 所以所以因此，中心差分方法为保持稳定性，必须满足因此，中心差分方法为保持稳定性，必须满足其中其中 是临界步长，是临界步长， 是系统最小固有周期是系统最小固有周期20kkkXX22112kkkkXtXX 11,kkkkXXXX（1）（2）（3）2210kp 22kkpt （4）21,22242kkpp2240kp 22kkpt 2 /kkT/ktT （5）（6）（7）ncrTtt crtnT（8）21,0,kMCK上页上页直接积分法直接积分

8、法 时域中心差分法注意事项时域中心差分法注意事项1. 中心差分法是显式中心差分法是显式算法，刚度阵不出现在递推公式的左端算法，刚度阵不出现在递推公式的左端2. 中心差分法是条件稳定算法。时间步长中心差分法是条件稳定算法。时间步长t必须必须小于由小于由该问题求解方程性质所该问题求解方程性质所 决定的某个临界值决定的某个临界值tcr，否则不稳定。稳定条件：，否则不稳定。稳定条件： 其中，其中，n是系统最高阶固有频率，是系统最高阶固有频率，Tn是最小固有周期。是最小固有周期。3. 显式算法用于求解梁、板、壳等结构的动态响应时，如果对角化后的质量矩显式算法用于求解梁、板、壳等结构的动态响应时，如果对角

9、化后的质量矩 阵阵M中已略去了与转动自由度相关的项，则中已略去了与转动自由度相关的项，则M的阶数仅对于位移自由度阶数，的阶数仅对于位移自由度阶数， 这时为了使这时为了使显式算法显式算法能够进行，可以采用主从自由度方法将转动自由度作为能够进行，可以采用主从自由度方法将转动自由度作为 从自由度在单元层次凝聚掉。从自由度在单元层次凝聚掉。4. 中心差分法比较适用于由冲击、爆炸类型载荷引起的波传播问题。中心差分法比较适用于由冲击、爆炸类型载荷引起的波传播问题。5. 对于结构动力学问题，一般不采用中心差分法。因为结构动力响应中主要成对于结构动力学问题，一般不采用中心差分法。因为结构动力响应中主要成 分通

10、常是低频，从计算精度考虑，允许采用较大的时间步长。分通常是低频，从计算精度考虑，允许采用较大的时间步长。 2ncrnTtt 直接积分法直接积分法 时域时域Newmark法法原理原理1kX基本假设：基本假设：的等效刚度阵由的等效刚度阵由M、C、K组成，求解时组成，求解时K阵需要求逆，阵需要求逆，隐式算法隐式算法联立联立21122111121112kKkkkkKkMXXXttKMC XFttCXXtXt11211112kkkkkkkkkXXXXtXXXtXXt 1111kkkkMXCXKXF直接积分法直接积分法 时域时域Newmark法法原理原理流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K

11、，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，参数，参数、并计算下列积分常数并计算下列积分常数计算等效刚度计算等效刚度阵阵2. 对每一时间步计算对每一时间步计算 tk+1时刻的等效力时刻的等效力 tk+1时刻的位移时刻的位移 tk+1时刻的速度和加速度时刻的速度和加速度基本假设：基本假设：联立联立 0000MXFCXKX0X21122111121112kKkkkkKkMXXXttKMC XFttCXXtXt11211112kkkkkkkkkXXXXtXXXtXXt 1111kkkkMXCXKXF012324567111,121,2 ,1,2aaaattttaaatat

12、 02311145kKkkkkKkM a Xa Xa XFFC a Xa XatX11kkKXF101231671kkkkkkkkkXaXXa Xa XXXa Xa X01KKa MaC直接积分法直接积分法 时域时域Newmark法稳定性讨论法稳定性讨论考虑无阻尼单自由度自由振动问题考虑无阻尼单自由度自由振动问题采用采用Newmark离散格式得到离散格式得到根据根据Newmark的假设公式，进一步的假设公式，进一步假定解的形式为假定解的形式为将（将（4）代入（）代入（3）式得到特征方程）式得到特征方程通过求解（通过求解（5）的特征方程得到）的特征方程得到其中，其中，20kkkXX222111/

13、2kkkkktXXtXt X 11,kkkkXXXX（1）（2）（3）21 21.2211/20kkkppp 22kkpt （4）21,2224 12ggh（5）（6）111 21/2211/20kkkkkpXpXX 1/21kkpgp1/21kkphp（7）21,0,kMCK上页上页直接积分法直接积分法 时域时域Newmark法稳定性讨论法稳定性讨论1. 解在小阻尼情况下应当具有振荡性质，因此解在小阻尼情况下应当具有振荡性质，因此为为复数，即有复数，即有24(1)(2)hg241/24kp 当当 很大时候，即很大时候，即 不受限制不受限制要求，上式成立必须要求要求，上式成立必须要求kpt21

14、1/242. 稳定的解必须不是无限增长的，因此稳定的解必须不是无限增长的，因此|0.5时候，时候，|0.5这一人为因素而引起的一种这一人为因素而引起的一种“人工人工”阻阻尼，被称为尼，被称为“数值阻尼数值阻尼”。 这种数值阻尼在一定条件下是有用的。在实际积分时候选取的这种数值阻尼在一定条件下是有用的。在实际积分时候选取的t通常远大于系统最高通常远大于系统最高阶频率对应的周期，而此对应的响应是不可靠的，并产生数值干扰，通过引入数值阻阶频率对应的周期，而此对应的响应是不可靠的，并产生数值干扰，通过引入数值阻尼使得高频干扰迅速尼使得高频干扰迅速衰减，这在分析冲击响应时较为不利。衰减，这在分析冲击响应

15、时较为不利。直接积分法直接积分法 时域时域Wilson-法法2266263322kkkkkkkkMXXXttMCK XFtttCXXXt1=kkkkkkkkMXCXKXFFFFF其中，原理原理流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，参数，参数(一般取一般取1.4)，并计算下列积分，并计算下列积分常数常数计算等效刚度阵计算等效刚度阵2. 对每一时间步计算对每一时间步计算 tk+时刻的等效力时刻的等效力 tk+ 时刻的位移时刻的位移 tk+1时刻的速度和时刻的速度和加速度加速度基本假设：基本假设：2322126kk

16、kkkkkkkkkkkkkXXXXtXXXXXtXXXXXXt联立联立 0000MXFCXKX0X注意事项注意事项在在t,t+t区间内任一时刻区间内任一时刻t+, Wilson- 法要求加速度、速度、位移法要求加速度、速度、位移满足满足当当 =t 时时26kkkkkkkkkkkXXtXXXXtXXXtXX 推出推出 后，代入求出后，代入求出 ，把，把 代入，并令代入，并令=t，可，可求得求得 、 、 。kXkXkX1kX1kX1kXkX的等效刚度阵由的等效刚度阵由M、C、K组成，求解时组成，求解时K阵需要求逆，阵需要求逆，隐式算法隐式算法0121322024567863,2,23,1,26ta

17、aaa atttaataaaaa 01KKa MaC0212kKkkkkKkM a Xa XXFFC a XXatXkkKXkkkkkkkkkkkkkXaXXa Xa XXXaXXXXtXaXX直接积分法直接积分法 时域平均值法时域平均值法(Nastran)2111212213111123311123kkkkkkMK XtMCK XFFFttMCK Xtt原理原理流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，计算下列积分常数，计算下列积分常数 计算计算 、 、 计算等效刚度阵计算等效刚度阵2

18、. 对每一时间步计算对每一时间步计算 tk时刻的等效力时刻的等效力 tk+1时刻的位移时刻的位移 tk时刻的速度和加速度时刻的速度和加速度112111111221313kkkkkkkkkkkkkkkXXXXtXXXtXXXXFFFFkkkkMXCXKXF1kX基本假设：基本假设：的等效刚度阵仅由的等效刚度阵仅由M、C、K组成，求解时组成，求解时K阵需要求逆，阵需要求逆，隐式算法隐式算法联立联立 0000MXFCXKX0X012030211,2,2aaaa aatt 1X0113Ka MaCK211311131313kkkkkka MK XFFFFa MaCK X1kkKXF0111112kkk

19、kkkkXaXXXXaXX起步计算起步计算：100110000XXXtFKXCXFKXCX1F0F直接积分法直接积分法中心差分法中心差分法Newmark法法Wilson- 法法平均值法平均值法(Nastran)直接积分法直接积分法 数值算例数值算例频率分析频率分析s88.10211Ts44. 4222Ts63. 3233T1231,2,33001002100030142000102260:0,0XXtXX 直接积分法直接积分法 数值算例中心差分法数值算例中心差分法30.36310Tt 0217.589at111.3772at20215.178aa3210.066aa流程流程1.初始计算初始计算

20、 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，并，并计算下列积分常数计算下列积分常数 计算计算 计算等效刚度阵计算等效刚度阵2. 对每一时间步计算对每一时间步计算 tk时刻的等效力时刻的等效力 tk+1时刻的位移时刻的位移 tk时刻的速度和加速度时刻的速度和加速度 0000MXFCXKX0X0120322111,2,2aaaa atat10030XXXta X 1X01Ka MaC2011kkkkFFKa M Xa MaC X1kkKXF0111112kkkkkkkXaXXXXaXX0100003000016X 0006X 210142022K

21、 100030001M直接积分法直接积分法 数值算例中心差分法数值算例中心差分法流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，并，并计算下列积分常数计算下列积分常数 计算计算 计算等效刚度阵计算等效刚度阵2. 对每一时间步计算对每一时间步计算 tk时刻的等效力时刻的等效力 tk+1时刻的位移时刻的位移 tk时刻的速度和加速度时刻的速度和加速度 0000MXFCXKX0X0120322111,2,2aaaa atat10030XXXta X 1X01Ka MaC2011kkkkFFKa M Xa MaC X1kkKXF

22、0111112kkkkkkkXaXXXXaXX003000000.363 00.066 0006000.396tXXtXa X 011007.589 0300017.58900022.7670007.589Ka MaC直接积分法直接积分法 数值算例中心差分法数值算例中心差分法流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，并，并计算下列积分常数计算下列积分常数 计算计算 计算等效刚度阵计算等效刚度阵2. 对每一时间步计算对每一时间步计算 tk时刻的等效力时刻的等效力 tk+1时刻的位移时刻的位移 tk时刻的速度和加速度

23、时刻的速度和加速度 0000MXFCXKX0X0120322111,2,2aaaa atat10030XXXta X 1X01Ka MaC2011kkkkFFKa M Xa MaC X1kkKXF0111112kkkkkkkXaXXXXaXX2010210100100014215.178 0307.589 0306022001001013.178107.5890141.534260213.178tttttttttFFKa M Xa MaC XXXX 00022.7670007.589ttX7.58900022.7670007.589tttXF直接积分法直接积分法 数值算例中心差分法数值算例中心

24、差分法流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，并，并计算下列积分常数计算下列积分常数 计算计算 计算等效刚度阵计算等效刚度阵2. 对每一时间步计算对每一时间步计算 tk时刻的等效力时刻的等效力 tk+1时刻的位移时刻的位移 tk时刻的速度和加速度时刻的速度和加速度 0000MXFCXKX0X0120322111,2,2aaaa atat10030XXXta X 1X01Ka MaC2011kkkkFFKa M Xa MaC X1kkKXF0111112kkkkkkkXaXXXXaXXtX1X2X30000t00

25、0.3945142t00.0346521.4755853t0.0045650.1928092.9674784t0.0333310.577864.5187475t0.1294481.2595115.8220876t0.3573982.2364196.7136827t0.7858263.4248657.2160118t1.4584194.678497.5098699t2.3631135.831877.8481610t3.4134796.751338.445663直接积分法直接积分法 数值算例数值算例Newmark法法30.3630.250.510Tt 0100003000016X 0006X 21

26、0142022K 100030001M流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，参数，参数、并计算下列积分常数并计算下列积分常数计算等效刚度阵计算等效刚度阵2. 对每一时间步计算对每一时间步计算 tk+1时刻的等效力时刻的等效力 tk+1时刻的位移时刻的位移 tk+1时刻的速度和加速度时刻的速度和加速度 0000MXFCXKX0X012324567111,121,2 ,1,2aaaattttaaatat 02311145kKkkkkKkM a Xa Xa XFFC a Xa XatX11kkKXF10123167

27、1kkkkkkkkkXaXXa Xa XXXa Xa X01KKa MaC0123456730.3565.51011.0191.01.00.00.18150.1815a a aa a aa a直接积分法直接积分法 数值算例数值算例Newmark法法流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，参数，参数、并计算下列积分常数并计算下列积分常数计算等效刚度计算等效刚度阵阵2. 对每一时间步计算对每一时间步计算 tk+1时刻的等效力时刻的等效力 tk+1时刻的位移时刻的位移 tk+1时刻的速度和加速度时刻的速度和加速度 0

28、000MXFCXKX0X012324567111,121,2 ,1,2aaaattttaaatat 02311145kKkkkkKkM a Xa Xa XFFC a Xa XatX11kkKXF101231671kkkkkkkkkXaXXa Xa XXXa Xa X01KKa MaC21010014230.36 03002200132.3610195.0720232.36K 0100003030.3611.021.06001tttttFXXX 30.3611.021.00.180.18tttttttttttttXXXXXXXXXttttKXF直接积分法直接积分法 数值算例数值算例Newmark

29、法法流程流程1.初始计算初始计算 生成刚度阵生成刚度阵K，质量阵，质量阵M和阻尼阵和阻尼阵C 计算计算 选择合适的选择合适的t，参数，参数、并计算下列积分常数并计算下列积分常数计算等效刚度阵计算等效刚度阵2. 对每一时间步计算对每一时间步计算 tk+1时刻的等效力时刻的等效力 tk+1时刻的位移时刻的位移 tk+1时刻的速度和加速度时刻的速度和加速度 0000MXFCXKX0X012324567111,121,2 ,1,2aaaattttaaatat 02311145kKkkkkKkM a Xa Xa XFFC a Xa XatX11kkKXF101231671kkkkkkkkkXaXXa X

30、a XXXa Xa X01KKa MaCtX1X2X30000t0.0002410.0078140.3713112t0.0027290.059121.391393t0.015390.2259772.8265424t0.0582930.5947464.3708265t0.1681631.2290125.7400396t0.3961912.1385676.755037t0.7962133.2658647.3870288t1.4047114.4951957.7509359t2.2185635.6817658.05037110t3.1792926.690798.494545直接积分法直接积分法 数值算

31、例结果比较数值算例结果比较0123456789024681012位移时间 (t)Newmark U1Newmark U2Newmark U3中心差分 U1中心差分 U2中心差分 U3直接积分法直接积分法 工程工程算算例例1Fig. 1 信号塔简化模型信号塔简化模型Fig. 2 地震载荷加速度曲线地震载荷加速度曲线单 元 类单 元 类型型CROD网格网格20 节点节点，72 单元单元材 料 属材 料 属性性E = 200 Gpa, A = 0.01 m2, = 0.3, = 7800 kg/m3边 界 条边 界 条件件底端全约束底端全约束载荷载荷x方向受如图方向受如图2所示地震载荷所示地震载荷时

32、 间 步时 间 步timeSteps = 5000, t = 0.002sFig. 3 信号塔顶部位移时程曲线信号塔顶部位移时程曲线(x方向方向)Fig. 4 动态效果图动态效果图直接积分法直接积分法 工程算工程算例例2Fig. 5 小雁塔简化模型小雁塔简化模型Fig. 7 小雁塔顶部位移时程曲线小雁塔顶部位移时程曲线(x方向方向)-1.5-0.50.51.50.00.40.81.21.62.0Acceleration(m/s2)t(s)单 元 类 型单 元 类 型CTETRA网格网格3734 节点节点，15272 单元单元材 料 属 性材 料 属 性E = 1790MPa, = 0.1， =

33、 1900 kg/m3边 界 条 件边 界 条 件底端全约束底端全约束载荷载荷x方向受如图方向受如图6所示的地震载荷所示的地震载荷时间步时间步timeSteps = 200, t = 0.01sFig. 6 地震载荷加速度曲线地震载荷加速度曲线Fig. 8 动态效果图动态效果图-0.08-0.0400.040.080.00.40.81.21.62.0Displacement(m)t(s)NastranNewmarkwilson-GenNewmarkBeta直接积分法直接积分法 工程算工程算例例3单 元 类 型单 元 类 型CHEXA网格网格4508 节点，4392 单元材 料 属 性材 料 属

34、 性E = 100GPa, = 0.3, = 3000 kg/m3边 界 条 件边 界 条 件底部以及前后两端全约束载荷载荷x方向施加如图10所示地震载荷时间步时间步timeSteps = 500, t = 0.002sFig. 9 某桥简化模型某桥简化模型 Fig. 11 大桥中间位移时程曲线（大桥中间位移时程曲线（x方向）方向）Fig. 10 地震载荷加速度曲线地震载荷加速度曲线Fig. 12 动态效果图动态效果图-4.0E-05-2.0E-050.0E+002.0E-054.0E-050.00.20.40.60.81.0Displacement(m)t(s)NastranNewmarkw

35、ilson-GenNewmarkBeta直接积分法直接积分法 工程算工程算例例4单元单元CBAR, CTRIAR( DKT Shell)网格网格2041 节点节点，1250 CBAR, 3530 DKT 材料材料C B A RE = 200GPa， = 0.3， = 7800 kg/m3, A = 0.3*0.3mCTRIARE = 2GPa， = 0.3， = 3900 kg/m3, t = 0.2m边界条件边界条件底部全约束底部全约束载荷载荷x方向受如图方向受如图14所示地震载荷所示地震载荷时间步时间步timeSteps = 1000, t = 0.002sFig. 13 房子简化模型房子

36、简化模型Fig. 15 房子顶部位移时程曲线房子顶部位移时程曲线(x方向方向)Fig. 14 地震载荷加速度曲线地震载荷加速度曲线Fig. 16 动态效果图动态效果图-1.0-0.50.00.51.00.00.51.01.52.0Acceleration(m/s2)t(s)-1.5E-03-1.0E-03-5.0E-040.0E+005.0E-041.0E-031.5E-030.00.40.81.21.62.0Displacement(m)t(s)NastranNewmarkWilson-GenNewmarkBeta基于基于ANSYS的动力响应分析的动力响应分析重力坝截面结构如图所示，坝高重力

37、坝截面结构如图所示，坝高120m，坝底，坝底宽为宽为76m，坝顶，坝顶10m，两侧的坡度如图所示。，两侧的坡度如图所示。大坝抗震性能分析的计算条件如下：大坝抗震性能分析的计算条件如下：（1）假设大坝的基础是嵌入到岩基中，地基）假设大坝的基础是嵌入到岩基中，地基是刚性的。是刚性的。（2）大坝采用的材料参数为：弹性模量）大坝采用的材料参数为：弹性模量35GPa，泊松比，泊松比0.2，容重，容重25KN/m3（3）计算分析大坝水位为）计算分析大坝水位为120m（4）水的质量密度）水的质量密度1000kg/m3（5）计算分析大坝在地震载荷下的响应）计算分析大坝在地震载荷下的响应基于基于ANSYS的动力

38、响应分析的动力响应分析自由度缩减方法自由度缩减方法 模态叠加法（振型分解法）模态叠加法（振型分解法） tMuCuKuf动力学方程TT2; MI K振型正交性uq令 tMqCqKqf代入 TTTtq C M Kqfq左乘T T2Ttq Cq q f整理自由度缩减方法自由度缩减方法 tMuCuKuf动力学方程 T2Ttq Cq q f整理后123n qq结构自由度n选取模态数qTq nn nn qq qMI作用1：降低计算自由度T2 K作用2：解耦T? C自由度缩减方法自由度缩减方法CMK瑞利阻尼T2 CI110miiiCM M K柯西阻尼T2 C比例阻尼12q自由度缩减方法自由度缩减方法 凝聚凝

39、聚法法静凝聚法静凝聚法：又称为盖恩法：又称为盖恩法(Guyan)子结构界面位移凝聚子结构界面位移凝聚：是一种动力子结构方法，将各子结构的刚：是一种动力子结构方法，将各子结构的刚度质量特征反映到其边界点或度质量特征反映到其边界点或“出口点出口点”。将所有子结构拼装成。将所有子结构拼装成整体之后，总刚度矩阵与总质量矩阵的阶数就是所有这些出口自整体之后，总刚度矩阵与总质量矩阵的阶数就是所有这些出口自由度的总和，达到降阶目的。由度的总和，达到降阶目的。动凝聚法动凝聚法：又称：又称库哈法库哈法(Kuhar)将结构中起主要作用的自由度称为主自由度将结构中起主要作用的自由度称为主自由度(xm)，把不，把不起

40、主要作用的自由度称为副自由度起主要作用的自由度称为副自由度(xs)，并通过合成凝，并通过合成凝聚矩阵将聚矩阵将副副自由度消除掉，以此降低结构计算自由度。自由度消除掉，以此降低结构计算自由度。自由度缩减方法自由度缩减方法 静凝聚静凝聚法法假定作用在副自由度上的外力为假定作用在副自由度上的外力为0，结构的静力方程写成分块形式，结构的静力方程写成分块形式sssmsmsmmmmKKx0KKxF将上述方程写成如下两个方程将上述方程写成如下两个方程ssssmmK xKx0mssmmmmK xKxF（1）（2）（3）将（将（2）式改写为）式改写为smxTx1sssm TK K其中其中（4）（5） KxF自由

41、度缩减方法自由度缩减方法将（将（4）式代入（）式代入（3）式，得到只包含主自由度的简化静力方程）式，得到只包含主自由度的简化静力方程其中其中 是凝聚刚度是凝聚刚度1mmmssssmKKK K KmmKxF（6）（7）KmssmmmmK xKxFsmxTx（4）（3）联立以下方程联立以下方程smxTxsmxxxmmTxxTxI（8）自由度缩减方法自由度缩减方法mxTxsssmsmsmmmmKKx0KKxFm0KxFmm0KTxFTTmm0T KTxTIFTTTmmT KTxFmmKxFTKT KT凝聚刚度矩阵凝聚刚度矩阵自由度缩减方法自由度缩减方法 静凝聚实例静凝聚实例1122233442000

42、00200000200000000 xxkkmxxkkkmxxkkkmxxkkm四自由度弹簧质量系统的广义特征值问题四自由度弹簧质量系统的广义特征值问题选取第二个和第四个质量为主自由度，将广义特征值问题重写为选取第二个和第四个质量为主自由度，将广义特征值问题重写为自由度缩减方法自由度缩减方法11332224420101000021101001120001001010001xxxxkmxxxx将上式写成分块形式将上式写成分块形式2sssmssssmsmsmmmmsmmmKKxMMxKKxMMx所以所以11102002021122sssmkkkkk TK K自由度缩减方法自由度缩减方法T10201

43、0211110102111122222111201100110222010101kkKT KTT1010002113110010011222422100101500110244000101mmMT MT自由度缩减方法自由度缩减方法31241544mM1121122k K根据缩减以后的质量阵和刚度阵，求解广义特征值问题根据缩减以后的质量阵和刚度阵，求解广义特征值问题缩减后的频率缩减后的频率原系统的频率原系统的频率210.1233722km221.1180071km210.1206148km22km相对误差相对误差2.3%11.8%自由度缩减方法自由度缩减方法 动凝聚动凝聚法法动动凝聚法与静凝聚法的区别凝聚法与静凝聚法的区别：静凝聚法只对系统的刚度矩阵：静凝聚法只对系统的刚度矩阵 应用应用凝聚技术，凝聚技术，动凝聚法对系统的动力矩阵动凝聚法对系统的动力矩阵 应用应用凝聚技术凝聚技术K2KM2sssmssssmsmsmmmmsmmmKKxMMxKKxMMx根据如下广义特征值问题根据如下广义特征值问题将上述方程第一行展开将上述方程第一行展开2ssssmmssssmmK xK xM xM x自由度缩减