2.5 算例分析
2.5.1 任意载荷振动问题分析
1．线性载荷
线性载荷的常微分方程为
初始条件为x（0）=0，
（0）=0。其精确解为x =e′ -e-2t /（2-t-1）。
图2.6所示为用Chebyshev谱元法求解线性载荷振动问题。在区间［0，6］上，对不同单元和不同插值次数，采用不同方案计算线性载荷振动问题的位移最大绝对误差，如表2.1所示。
图2.6 Chebyshev谱元法求解线性载荷振动问题
表2.1 采用不同方案计算线性载荷振动问题的位移最大绝对误差
从表2.1可以看出，Chebyshev谱元法计算的线性载荷振动问题有很高的精度，单元数从10到100，插值次数从10到500，位移最大绝对误差均达到10-10数量级。从图2.7可以进一步看出，本章方法从h收敛和p收敛两个方面都获得了很高的精度，但是p收敛更稳定。而配点法，插值次数从10到5000，最高精度达到了10-8数量级，但随着插值次数增加，误差越来越大，甚至获得了错误的解，当插值次数为5000时，位移最大绝对误差为396.928。
图2.7 线性载荷作用振动问题Chebyshev谱元法两种收敛曲线
2．三角载荷
线性无阻尼系统在三角外力的作用下振动，其振动微分方程可以描述为
初始条件为x（0）=0、
，三角外力为
式中，F0=100，t0=0.4。
振动微分方程的精确解为
式中，k、tω分别表示系统的刚度和振动周期。
采用不同的单元数和插值次数分析0～1.2s内的位移，得到不同方案计算三角载荷振动问题的位移最大绝对误差，如表2.2所示。图2.8所示为采用Chebyshev谱元法求解三角载荷振动问题。
图2.8 采用Chebyshev谱元法求解三角载荷振动问题
表2.2 采用不同方案计算三角载荷振动问题的位移最大绝对误差
从表2.2中可以看出，用Chebyshev谱元法求解三角载荷作用下的振动问题可获得比求解线性载荷作用下的振动问题更高的精度。当单元数为12，插值次数为300时，位移最大绝对误差为4.146×10-12；最小的最大绝对误差是单元数为18、插值次数为12时，为7.549×10-15；最大的最大绝对误差是单元数为10、插值数为10时，为0.0036。而配点法求得的最小的位移最大绝对误差为4.267×10-6，最大的位移最大绝对误差为0.03118。因此，Chebyshev谱元法比配点法有更高的精度。从图2.9可以进一步说明Chebyshev谱元法不仅精度高，而且稳定。比较图2.9（a）和图2.9（b），可以看到h收敛更稳定。
图2.9 三角载荷作用振动问题Chebyshev谱元法两种收敛曲线
3．半正弦波脉冲载荷
控制方程为
初始条件为x（0）=0，
，ωn=1。半正弦波脉冲载荷为
其精确解为
当t1=10时，采用Chebyshev谱元法计算半正弦波脉冲载荷作用下的振动问题，可以得到如图2.10和表2.3的结果。
图2.10 Chebyshev谱元法求解半正弦波脉冲载荷作用下的振动问题
表2.3 采用不同方案计算半正弦波脉冲载荷作用下振动问题的位移最大绝对误差
从表2.3中可以看出，对Chebyshev谱元法单元数比较小，插值次数比较大时，所得结果精度低，如单元数为6，插值次数为100时，位移最大绝对误差为0.00018；而单元数比较大，插值次数比较小时，所得结果精度高，如单元数为100，插值次数为12时，位移最大绝对误差为4.961×10-15。对于配点法，当插值次数为2400时，可获得最高精度，其位移最大绝对误差为7.957×10-7；当插值次数为40时，精度最低，位移最大绝对误差为0.00771。总之，用Chebyshev谱元法，只要合理选择单元数和插值次数，就可以获得10-10数量级以上的精度。
4．悬臂梁
在空间上由传统的有限元法离散微分方程，形成关于节点变形二次微分方程［7］［7］：
式中，｛w｝为节点变形。每个节点包括6个变形、3个位移｛x，y，z｝和3个旋转｛θx，θy，θz｝，这里没有考虑阻尼。悬臂梁的几何尺寸和物理参数，如表2.4所示。
表2.4 悬臂梁的几何尺寸和物理参数
悬臂梁空间上可用Euler梁单元来离散，每个节点有6个自由度。这里把悬臂梁离散为10个单元，即11个节点，左端固定，共有60个自由度。悬臂梁端点受力为f（t）=bsin（ωt），方向竖直向上。
当载荷频率ω为10π、176.64π、184π时，用Chebyshev谱元法求解获得悬臂梁自由端竖直方向的位移响应，并且与ANSYS的分析结果进行比较，如图2.11所示。
图2.11 Chebyshev谱元法求解悬臂梁在正弦载荷作用下的振动问题
图2.11 Chebyshev谱元法求解悬臂梁在正弦载荷作用下的振动问题（续）
从图2.11可以看出，对于悬臂梁在正弦载荷作用下的振动问题，应用Chebyshev谱元法所得的结果与ANSYS的结果一致性较好。
2.5.2 聚集单元谱元法
为了验证聚集单元谱元法的优势，本章对4个不同实例进行分析，并与等距单元谱元法进行比较。
1．标准形式
对于突变载荷，常规的等距单元谱元法不能获得好的收敛效果。对于图2.12所示的冲击载荷，无论多么复杂的系统，其动态响应方程最终可以转化为式（2.57）所示的标准形式。其中，在0.5s时，产生很大的冲击力，冲击力的幅值可以由系数决定，冲击力的作用时间可以通过调整ε来改变，本例中ε=0.0001125。解析解可以通过积分因子法并借用MATLAB中的erf函数获得：
图2.12 冲击载荷
式中，
。
图2.13所示为等距单元谱元法的计算结果，很明显其误差很大，特别是在载荷突变处，当Nel=30时，误差极大；而当Nel=100时，误差减小了很多，但还是较大。图2.14所示为聚集单元谱元法的计算结果，在单元数和插值次数相等的前提下，该方法可消除误差。图2.15和图2.16更能说明聚集单元谱元法的优势。图2.15将等距单元谱元法与聚集单元谱元法进行比较，对单元数为50的情况，当单元插值次数为18时，前者误差为0.02506，后者误差为1.382×10-10，优势明显。图2.16将聚集单元谱元法与等距单元谱元法进行比较，发现对单元数为80的情况，当单元插值次数分别为5、10、18时，聚集单元谱元法误差分别为3×10-9、3.628×10-12、2.668×10-12，而当单元插值次数为18次时，误差为0.01124，可以看出聚集单元谱元法在求解冲击载荷响应有很大优势。
图2.13 等距单元谱元法计算结果
图2.14 聚集单元谱元法计算结果
图2.15 聚集单元谱元法与等距单元谱元法比较（单元数为50）
图2.16 等距单元谱元法与聚集单元谱元法比较（单元数为80）
2．线性单自由度系统
图2.17所示为线性单自由度系统。其中，固定质量m=1kg，弹簧刚度系数为k=0.9，阻尼器系数为c=0.9。当t=0时，系统以v=1m/s的速度撞击在一个固定的阻碍物上，同时所受冲击载荷
。系统的运动方程为
图2.17 线性单自由度系统
式（2.58）的自由振动解析解为
线性单自由度系统虽然既有初速度，又同时受冲击载荷的作用，但冲击载荷的中心在6s处，且其宽度由ε=0.01125来控制，很窄，因此在4s附近系统还属于自由振动。图2.18所示为线性单自由度系统的动态位移响应，从图中可以看出，当单元数为12，单元插值次数为12时，聚集单元谱元法明显比等距单元谱元法精度高。
图2.18 线性单自由度系统的动态位移响应
3．杆桁架结构
124杆桁架结构包含49个铰链、94个自由度（见图2.19）。其弹性模量E=207GPa，泊松比v=0.3，密度ρ =7850kg/m3，杆的截面积为0.645×10-4m2。在节点1、20、19、18、17、16、15的X正方向上作用相同的动态载荷，在节点1、2、3、4、5的y负方向上也作用相同的动态载荷。动态载荷
，其中，b=1000，ε=0.0001125。
图2.19 124杆桁架结构
从工程角度，将节点1、2、3、4、5、15、16、17、18、19、20作为重点考察位置（称为关键位置），同时也是冲击载荷作用的位置。在冲击载荷的作用下，124杆桁架结构关键位置节点X方向的动态位移响应如图2.20所示。将其中节点1X方向的位移单独考虑，并局部放大，如图2.21所示，可以明显看出，等距单元谱元法的计算结果在部分点偏离精确值，而聚集单元谱元法没有出现这种现象。实际上，关键位置的每个节点都出现了类似现象，限于篇幅，没有全部列出。
图2.20 124杆桁架结构关键位置节点X方向的动态位移响应
图2.21 124杆桁架结构节点1X方向的动态位移响应
4．连杆小头
连杆小头及其相连的杆身部分如图2.22所示。其采用PLANE42单元，右端固定，有111个节点、212个自由度。连杆小头内圆中间所受的水平向右的冲击载荷，
、ε=0.0001125。连杆材料的弹性模量E=207GPa，泊松比v=0.3，密度ρ =7850kg/m3。
图2.22 连杆小头及其相连的杆身部分
将连杆小头圆角处作为关键位置，考察节点15、42、49、93在X方向的动态位移响应。在图2.23中，若将模态叠加法获得的结果看作精确解，那么可以看出，在相同的单元数、尺寸、插值次数下，聚集单元谱元法的结果更精确，和模态叠加法保持一致，而等距单元谱元法的结果存在较大误差。由于位移响应本身很小，图2.23也进行了局部放大。
图2.23 连杆小头4个关键位置节点X方向的动态位移响应
2.5.3 非线性振动分析
1．Duffing型非线性振动方程
Duffing型非线性振动方程可以写为
式中，ε、F是给定的常数；ω是外载荷的频率，也是常数。
（1）当ε=-1/6、F=0、ω=0.7时，初始条件为
式（2.60）的近似解析解为
此时，Duffing型非线性振动问题的响应和Newton-Raphson迭代过程分别如图2.24和图2.25所示。
图2.24 Duffing型非线性振动问题的响应（第一种初始条件）
（2）当ε=-1/6、F=2、ω=1时，初始条件为
式（2.60）的近似解析解为
图2.25 Duffing型非线性振动问题的Newton-Raphson迭代过程（第一种初始条件）
此时，Duffing型非线性振动问题的响应和Newton-Raphson迭代过程分别如图2.26和图2.27所示。
图2.26 Duffing型非线性振动问题的响应（第二种初始条件）
图2.27 Duffing型非线性振动问题的Newton-Raphson迭代过程（第二种初始条件）
从图2.24和图2.26可以看出，本章方法获得的结果与近似精确解非常吻合。从图2.25和图2.27可以看出，Duffing型非线性振动在两组不同参数下都有很好的收敛。
2．单摆的非线性振动
单摆的非线性振动方程为
式中，g为重力加速度；l为摆长；θ为摆角。初始条件为θ（0）=θ0，
。
当单元数为10，插值次数为6时，通过Galerkin离散方案得到非线性方程组，然后利用Newton-Raphson法求解，可得到当初始摆角θ0=π/5时的摆角、角速度和角加速度，如图2.28所示，并且将其与ODE45求解器的计算结果比较，两者很吻合。图2.29为单摆非线性振动的Newton-Raphson迭代过程，可以看出其非线性有很好的收敛。
图2.28 单摆非线性振动的响应θ0=π/5
图2.29 单摆非线性振动的Newton-Raphson迭代过程
求出摆角响应θ（t）后，可以获得两个时间点ti、t j，它们满足θ（ti）＞0，θ（t j）＜0且ti＜tj。在区间［ti，tj］上进行重心Lagrange插值，然后采用二分法求解θ（t）=0，可获得tθ=0，那么具有某一初始摆角的单摆非线性振动的周期T=4tθ=0。将其和单摆线性振动的周期
相比如表2.5所示，并与精确解、二阶摄动解及DQ法进行比较。
根据表2.5的数据可得，当初始摆角θ0＜135°时，谱元法可以获得最大的绝对误差为0.0013，而二阶摄动解最大的绝对误差为0.1042，DQ法最大的绝对误差为0.0016；当θ0=150°时，谱元法获得最大的绝对误差为0.0016，而二阶摄动解最大的绝对误差为0.1585，DQ法最大的绝对误差为0.0125。
表2.5 单摆非线性振动的初始摆角和固有频率的比值
