euler中点隐式差分格式及其保辛算法

euler中点隐式差分格式及其保辛算法

数值优化算法在动态方程的研究中，r-k方法、中心差异分法、绅士模型方法和新格式法被用来分解系统的r-k方法、中央差异分法、精确误差法和新格式法。这些强大的算法本身消耗系统的能量，导致动态响应的相位延迟，因此不知道该怎么办的长期跟踪能力。为了提高计算的精度，可以通过降低时间间隔来实现。降低分散相位需要降低计算量和工程应用价值的成本。冯康等建立了哈密顿系统的辛几何算法,从理论上清楚地阐明了传统数值方法导致能量耗散的根本原因,即其截断项是耗散项.而辛算法对应的差分格式严格保持哈密顿系统的辛结构,有限阶辛算法的截断部分不会导致系统能量发生线性变化.在长时间数值稳定性方面,辛算法具有独特的优越性.关于辛算法近几年的发展和应用在文献中已有综述和回顾.对于多数结构动力学分析问题都要借助于有限元方法.因此,研究用辛差分格式求解有限元动力学平衡方程的有关问题非常重要.本文用Euler中点隐式差分格式来求解动力学方程有两个目的,一是研究用对保守系统保辛的差分格式来求解非保守系统时算法的特性及其求解精度,包括算法的Jacobi矩阵性质,算法是否保能量等;二是研究辛算法及其在计算稳态响应时的精度问题.Kane等证明了隐式中点格式和Newmark格式都是变分积分格式,对某些参数Newmark格式为辛形式.本文在用Euler中点隐式格式求解有限元动力学平衡方程的公式推导过程中,直接证明了(δ=0.5和α=0.25的Newmark算法就是Euler中点隐式差分格式.1算法的特性分析考虑如下线性动力学方程本文只考虑结构刚度矩阵K,质量矩阵M和阻尼矩阵C皆为对称正定的情况.引入广义动量方程(1)变为令状态向量为zT={xTyT},方程(2)可以写为式中方程(3)可以用Euler中点隐式差分格式来求解,即式中:h=tk+1-tk(k≥0)为差分步长,差分格式(5)的精度为O(h2).把式(4)代入方程(5)中得到式中:Jacobi矩阵,矩阵L和广义载荷列向量γk+1分别为式中:得到了广义位移xk+1和广义动量yk+1之后,再根据式(2)的第2式就可以得到广义加速度式(6)和式(7)就是用Euler中点隐式差分格式求解动力学方程(1)的递推格式.需要的已知条件为初始速度和初始位移,不需要初始加速度.下面分两种情况对算法的特性进行讨论.情况1C=0(1)对任意差分步长h,不难证明矩阵A满足下面辛矩阵条件即矩阵A为辛矩阵.(2)由于A为辛矩阵,即ATJA=J,其中J为标准辛矩阵,因此可以得到即辛矩阵A与其逆矩阵A-1具有相同的本征值λ,因此λ2=1.若λ为实数,则其只能等于1或-1;若λ为复数,则必成对出现为eiθ和e-iθ.由此可以得到下面的结论,辛矩阵A所有本征值的模为1,其谱半径为1.容易证明,当h→0或h→∞时,辛矩阵A是正交辛矩阵.当h、0时,A的所有本征值为1,当h→∞时,其所有本征值为-1.由于辛矩阵的行列式为1或-1,不存在零本征值,因此辛矩阵A的本征值在从-1到1的变化过程中,必须经过复数域.(3)不难证明在积分过程中线性保守系统能量是严格守恒的,即Ek+1=Ek,其中(4)由于矩阵A与R无关,因此即使R≠0,矩阵A一定也是辛矩阵.当载荷列向量R=0时,式(2)就是Hamilton系统正则方程.情况2C≠0(1)当阻尼C为比例阻尼时,对任意差分步长h,容易证明矩阵A满足式(())给出的关系.(2)根据平衡方程(1)对式(6)的第1式进行变换可以得到不难看出,式(11)就是当δ=0.5,α=(0.25时的Newmark算法由此可见,δ=0.5和α=0.25的Newmark算法就是Euler中点隐式格式,对保守系统是一种辛算法.Newmark算法的经典计算格式需要初始加速度,而格式(6)不需要初始加速度.2时间步长的确定考虑如下动力学方程设系统的初始条件为,简谐激励为RT={010}sinωt.考虑如下正定比例阻尼C1和负定比例阻尼C2:系统的固有频率和固有模态向量分别为和.系统(12)固有振动的最小周期T2为2π/ω2=2.8,因此下面选用的时间步长若是0.28,则相当于是T2/10.2.1jacobi矩阵a的谱半径无条件稳定算法的Jacobi矩阵的谱半径ρ对任何时间步长都要小于等于1.2.1.1-110-n对于任意比例阻尼矩阵,式(9)是自动满足的.若步长h=10-n(n为正整数),则|ρ-1|≈10-n,即随着步长的减小,式(8)逐渐被满足,A逐渐逼近于辛矩阵,其所有本征值的模皆趋近1,ρ也逐渐趋近于1,参见图1和表1.2.1.2=0.5的rema算法的特征引入广义动量y=Mx,同样可以把Newmark差分格式变成为与式(6)相同的递推格式,图2给出了不考虑阻尼时其Jacobi矩阵的谱半径ρ和参数δ的关系.值得强调的是,即使不考虑阻尼,式(9)也得不到满足,因此δ≠0.5的Newmark算法与Euler中点格式也就是与δ=0.5的Newmark算法具有本质的区别,当步长h逐渐减小时,A也逼近于辛矩阵,其所有本征值的模趋近1.从图2可以看出,谱半径ρ也逐渐趋近于1.这意味者,随着步长h的减小,算法固有的“人工阻尼”也会愈来愈小,但永远无法从根本上消除.因此对于δ≠0.5的Newmark算法,只有通过减小步长h才能减小人工阻尼,以提高计算精度,但代价是计算量的迅速增大,这种做法的实用价值是有局限性的.不过对于合适的差分步长,算法的阻尼对滤掉响应的高频成分和暂态响应是有作用的.2.2包辛和人工电阻为了清楚起见,下面分几种情况进行讨论.2.2.1阻尼系统的附加(1)保守系统的自由振动对于保守系统,Euler中点格式能够保证线性系统的能量严格守恒,由此可以得到一个重要的结论,那就是Euler中点辛格式本身不存在任何“阻尼”,否则它不会使保守系统的能量严格守恒,这不同于一般的R-K方法、Wilson-θ方法和δ≠0.5的Newmark方法.(2)非保守系统的自由振动把Euler中点格式用于阻尼系统,矩阵A不再是辛矩阵.虽然Euler中点格式本身不存在阻尼,但把它用于阻尼系统时,相位滞后也会给系统带来附加的阻尼.对于正阻尼系统,相位滞后使阻尼耗散的能量滞后,这相当于相位滞后给系统带来附加的负阻尼,参见图3,其中E为系统的总机械能;对负阻尼系统,相位滞后使阻尼向系统输入的能量滞后,这相当于给系统提供了附加的正阻尼,参见图4.因此,对于非保守系统,Euler中点隐式格式也具有广义的保“结构”特征,也就是该格式总是有要保持系统能量不变的特性.而δ≠0.5的Newmark方法本身是存在阻尼的,它耗散系统能量的特性与物理系统本身阻尼特性无关.值得指出的是,虽然保辛算法能够使保守系统的能量守恒,但这并不意味着保证位移或速度振幅不变,因为位移和速度响应振幅同时变化也会使系统的能量守恒,当步长h比较大时就会出现这种情况.2.2.2动态响应的特性在下面的讨论中,不考虑初始条件.(1)无阻尼系统的受迫振动Euler中点格式:对于无阻尼的受迫振动系统,矩阵A是辛矩阵,不论载荷性质如何,其误差主要是存在相位滞后.对于简谐激励作用,用该格式可以求得包括暂态响应在内的一般解.若只求解稳态响应,需要施加人工阻尼,但幅值和相位的准确性得不到保证.δ≠0.5的Newmark方法:算法阻尼使暂态响应逐渐衰减至忽略不计,因此只能求得稳态解.值得指出的是,δ的变化对稳态响应基本没有影响.(2)阻尼系统的受迫振动对于一般非交变激励作用下用不同方法求解的动态响应的特性,与阻尼系统自由振动响应类似,这里不再赘述.对于简谐激励作用下系统的响应,由于系统的暂态响应逐渐衰减,因此用不同格式得到的响应都是稳态响应,只是精度不同.伴随自由振动中有频率比因子,ωi为系统的固有频率,的大小直接决定在振动初始阶段稳态响应的幅值和暂态响应幅值的相对大小,也决定着不同差分方法的求解精度.当时,暂态响应的幅值通常不大于稳态的幅值,Euler中点方法和δ≠0.5的Newmark方法求解的稳态响应振幅精度都比较高,但相位滞后的程度不同;当时,暂态响应的幅值通常不小于稳态的幅值,算法人工阻尼使伴随自由振动不正常衰减,致使δ≠0.5的Newmark方法的相位超前,参见图5;当ω1≤ω≤ωn时,δ≠0.5的Newmark方法振幅衰减严重,相位可能超前或滞后,参见图6.值得指出的是,在的不同取值区间内,Euler中点方法的振幅和相位的求解精度都要高于δ≠0.5的Newmark方法,两种方法的相位误差都不累积,但前者对频率比不敏感.当时,为了保证数值积分精度,积分步长要由激励周期来确定.另外还有个现象,在求解简谐激励响应时,保持步长不变,对于Euler中点格式阻尼系统的求解精度大于无阻尼系统的求解精度,这是因为系统阻尼可以滤去暂态响应;对于δ≠0.5的Newnark方法,无阻尼系统的求解精度比阻尼系统的求解精度高,这是因为不论δ大小,算法阻尼也都能滤去系统的暂态响应.3算法的相位误差本文给出了用Euler中点隐式差分格式求解动力学方程的递推格式,得到了如下结论:(1)辛矩阵所有本征值的模等于1;Euler中点辛格式使线性保守系统能量严格守恒;对于阻尼系统自由振动问题,Euler算法相位误差给系统带来了附加能量,但其总是具有试图保持系统总能量不变的特性;对于自由振动、非简谐受迫振动,Euler算法和Newmark算法都存在相位误差累积;其他结论参见表2.(2)证明了δ=0.5和α=0.25的Newmark算法就是Euler中点辛差分格式,δ≠0.5的Newmark算法相位误差存在超前或滞后两种可能性,但