非线性动力学方程的求解方法.doc

非线性动力学方程的求解方法1、概述在工程实际问题中，我们常常面临这样的选择：我们所遇的问题究竟是静力的还是动力的。静力问题与动力问题，从力学的角度看就是是否考虑与加速度有关的力，而从数学求解方法看则是一个三维边值问题还是一个四维边值-初值问题。在这个问题的选择上没有固定的原则，一般取决于我们研究者、分析者对工程问题的判断。一般认为，实际工程大都是处于动力环境之中，因而属于动力问题。但是，由于时间、经费等方面的原因的限制，我们不可能把所有的问题都按照动力问题的方法来分析。对于许多具体的问题，与速度和加速度有关的力足够小，但是又影响结构分析结果的，将采用静力假定来模拟这些力。线性的动力有限元控制方程如式(1-1)所示。 (1-1)式中MDK分别为结构的质量、阻尼和刚度矩阵，R为荷载列矢量，、和分别是加速度、速度和位移列矢量。式(1-1)的解法大体上可以分为两类：直接积分法和模态叠加法。直接积分法在对控制方程进行数值积分之前不对方程做任何形式的变换，直接用数值积分的方法在时域上一步一步地对方程进行积分。模态叠加法是在求解之前对方程进行某种数学变换，使基底降低，或使矩阵的带宽减小，再进行求解。这两种方法在形式上不同，但是密切相关。上述每一类求解方法中又有许多具体的解法，每一种解法又有各自的特点。因此我们在选择一种方法求解一个问题时，要对该方法的收敛性、稳定性、效率、精度和费用等进行一些分析，讨论它对所求问题的有效性，从而使我们能够针对某一特定的问题，选择合适的方法。直接积分法基于以下两条：(1)不是在求解时间区间内任意时刻t都满足式(1-1)，而是在相隔t上的一些离散时刻满足式(1-1)。(2)对位移、速度和加速度在每一时间区间t内变化的形式进行假设，事实上若把式(1-1)看成一个常系数微分方程组，便可以用任何一种有限差分格式通过位移来近似表示速度和加速度，因此不同的差分格式就得到不同的方法。从差分格式上看，分显式与隐式两大类方法，所谓显式差分法就是不必对方程进行求解，而是由前一时刻t的平衡条件直接可求解t增量后t+t时刻的各参数解。而隐式差分法则必须对方程进行求解。显式差分法有中心差分法等。隐式差分法有Wilson-法和Newmark法。不同方法解的精度、稳定性、收敛性、效率不同。2、中心差分法 2.1中心差分法简述中心差分法的差分格式为： (2-1) (2-2) 认为和在任意时刻t上满足平衡方程(1-1)，即： (2-3)将式(2-1)和式(2-2)代入到式(2-3)中，可得：(2-4)2.2中心差分法的特点(1)的解是基于利用在时刻t的平衡条件，即是在假定t时刻式(2-3)成立的条件下来计算的，该积分过程称为显式积分法，因此中心差分法为显式差分法。但是，用式(2-4)计算时，不仅要用到而且还要用到。所以，计算在时刻t+t的解，必须存在一个具体的起始过程，即必须利用初始条件。这样，已知和，可由式(2-4)求得，由式(2-1)和式(2-3)可求得，如式(2-5)所示。 (2-5) 我们称式(2-5)为差分格式(2-4)的初始条件。(2)在求解时不需对刚度矩阵K进行三角分解。在应用此法时一般采用集束质量矩阵或称对角质量阵，而阻尼矩阵也通常为对角形式，这样，在运用式(2-4)时，就不需要对等号右端的系数矩阵)进行三角分解，从而可以节省计算时间。 如果不考虑系统的阻尼，则，这时式(2-4)可简化为： (2-6) 其中： (2-7)从式(2-6)和式(2-7)可以看出，若质量阵M为对角矩阵，即在系统中只考虑集束质量的情况，则方程(2-6)实际上是解耦的，即方程组中是相互独立的，因此只需进行矩阵相乘就可得到。把式(2-6)改写为： 写成分量形式，即： (2-8) 其中，和分别为和的第i个分量，而为质量矩阵的第i个对角线元素。对于，我们可以用静凝聚的方法进行处理。 若对总体刚度和质量矩阵不必进行三角分角，就意味着不必形成总体的K和M，因此，式(2-7)的计算只需在单元一级进行即可，即： (2-9) 式中Ke和Me分别为与待求节点自由度有关的单元刚度矩阵与质量矩阵。所以，使用式(2-8)和式(2-9)构成的中心差分具有明显的优点。它不需要计算总体刚度矩阵K和质量矩阵M，求解过程在单元一级进行，计算效率较高。而且，如果所有单元刚度和质量矩阵都相同时，或有大部分的单元都相同时，可以用一个或少数单元矩阵来求解整个系统。在工程中，常会遇到重复对称式结构，这时用这种方法来求解将具有很高的效率。(3)加速度和速度的差分格式都具有(t2)阶的误差。(4)可以证明，中心差分法所要求的步长为： (2-10)其中，Tn为有限无限无集合体的最小周期，n为单元系统的阶，为临界时间步长。 周期Tn可利用如下任何一种方法来计算：(1) 广义Jacobi法；(2) 逆迭代法；(3) 正迭代法；(4) Rayleigh商迭代法；(5) 多项式迭代法；(6) 斯图姆序列对分法；(7) 行列式搜索法；(8) 子空间迭代法。事实上，Tn为系统最大特征值的倒数。要求使用的时间步长t小于临界时间步长的差分格式，称为条件稳定的。若使用一个大于的时间步长，则积分将是不稳定的，此时由数值积分在计算机上导致的舍入误差会增大，从而使结果失去意义。中心差分法是一种条件稳定的方法，这一缺点导致了它在使用上具有很大的局限性。首选是时间步长选择困难，当矩阵阶数较高时，而且同时原系统的质量分布不均匀，将使得M矩阵的个别对角上的元素很小，从而使系统的高阶特征趋于无限大，Tn趋于零，这样便使积分过程成为不可能。所以，仅仅因为一个十分小的质量元素，就会导致计算效率的降低。当这种情况出现时，通常采用静凝聚的方法使矩阵的阶数降低，但这仅仅能在低频段模拟原系统。因此，中心差分法是条件稳定的非线性动力学分析方法，有其局限性。3、Wilson-法Wilson-法是线性加速度方法的扩展，基于以下两条假定：(1)当时间从t至变化到时间()时，加速度从变化为；(2)在时间区段内，结构的刚度、阻尼、地面运动加速度都无改变。则由上述两条假定可得： (3-1) 式中，为参数，可根据积分的稳定性和精度进行选取，一般要求大于1。 对式(3-1)两端的进行积分，可得： (3-2) 再对式(3-2)两端的进行积分，可得： (3-3) 在式(3-2)和式(3-3)中令，则可得在()时刻的速度和位移分别如式(3-4)和式(3-5)所示： (3-4) (3-5) 将加速度和速度写成由位移表达的形式，如式(3-6)和式(3-7)所示： (3-6) (3-7) Wilson-法认为在()时刻，系统满足动力平衡方程(1-1)，并由Wilson-法的第二条假定，我们把式(3-6)和式(3-7)代入到式(1-1)中，可得到式(3-8)。 (3-8) 其中， (3-9)求()时刻的位移、速度和加速度是我们所要得到的。求解的过程如下所示：(1) 将式(3-6)和式(3-7)代入到式(3-8)中，并加以整理可得：= + (3-10)由式(3-10)可以求得。(2) 将式(3-6)分别代入式(3-1)、式(3-2)和式(3-3)中，并令，可得： (3-11) (3-12) (3-13)Wilson-法的特点：(1) 该法是一种隐式积分法，因为刚度矩阵K是未知位移矢量的系数矩阵。(2) 该法不需要特别的初始过程，因为在时刻()的位移、速度和加速度只是利用在时刻t的相同的量来表示。(3) 可以证明，当时，该法为无条件稳定的差分格式。一般情况下，取。目前在许多大型通用程序中，其缺省值都为1.4。(4) 由于不必为选取的大小而担心，所以该法在使用上是比较方便的，即使有时可能选取的值较大，但是结果仍是稳定的，影响的只是精度。(5) 由于其是隐式差分格式，所以每一步都必须对方程求解，因而计算效率相对于中心差分法要低一些。4、Newmark法Newmark法也可以认为是线性加速度法的一种发展。其通常采用的平均常数加速度方案。Newmark法的差分格式为： (4-1) (4-2) 式(4-1)和式(4-2)中和是参数，根据积分的精度和稳定性的要求来确定这两个参数。当=1/2, =0时，加速度是常数，并在任一上等于；当=1/2, =1/8时，加速度是常数，且等于由变化到的时间间隔中点的值；当=1/2, =1/16时，加速度由变化到的变化是线性的；当=1/2, =1/6时，加速度是常数，且等于由变化到的时间间隔中点的值；当=1/2, =1/4时，相当于假定加速度保持平均值为(+)/2。我们考虑时刻()的平衡方程： (4-3)由式(4-2)通过求出，然后把代入式(4-1)中，就得到和的方程，它们仅仅通过位移来表示。把这两个关于和的关系式代到到式(4-3)中，求出，然后利用式(4-1)和式(4-3)就可以算出和，最终方程为：=+= (4-4) (4-5) (4-6) 将它们代入到运动方程(4-3)，得：=+ (4-7)由式(3-6)，解得后，代入式(4-4)和(4-5)，求得()时刻的速度和加速度。Newmark法的特点：(1)该法是一个隐式积分法，因为刚度矩阵K是未知位移矢量的系数矩阵，必须求解式(4-4)才可得到。(2)不需要特别的初始过程，因为在时刻()的位移、速度和加速度只是利用在t时刻的量来表示。(3)可以证明，纽马克法当和时，积分格式是无条件稳定的。一般取和=0.25即可。在大型通用程序中，、的缺省值也分别为0.5和0.25。但对于动力接触问题，、应分别取0.5。(4)由于该法是无条件稳定的，所以避免了选取上的麻烦。同样，由于每一步都必须对原方程进行求解，尽管对K和M的三角分解只需进行一次，但每一步都必须进行回代。所以，计算效率对有条件稳定的中心差分法要低。(5)该法与Wilson-法的求解过程基本相同，因而求解效率也基本相同。5、直接积分法的应用 对于非线性的情况，最通常的方法是采用纽马克法(梯形法)，并取=0.5和=0.25，把刚度矩阵K换成切线刚度矩阵，用代替，而且在每一时间步进行迭代。动力平衡方程改写为： (5-1)其中： (5-2)右上表K及K+1代表迭代次数。 由梯形法则有： (5-3) (5-4) 由式(5-2)、式(5-3)和式(5-4)可得： (5-5) 把式(5-2)代入式(5-1)可得： (5-6)其中： (5-7) 式(5-1)至式(5-7)的是时刻的切线刚度矩阵。式(5-6)与静力非线性迭代方程一致，可用各种Newton-Ramphson法求解。6、模态叠加法 对直接积分法的稳定性、精度和效率等方面的讨论可以看到，隐式积分法在进行每一步的积分时，都必须按照时间增量来求解动力平衡方程，也就每一时间步长都要象求解静力问题那样对方程左端系数矩阵进行三角分解，对于大型矩阵或带宽的矩阵来说，计算量是相当大的。即使对于显式积分法的中心差分法来说，由于为了保证计算的稳定性，必须选择得很小，所以一步步的计算量加起来的总计算量也很可观，而且还存在由于积累误差而造成的精度不高的问题。因此，如何找到大型矩阵或带宽较大的矩阵所构成的式(1-1)的有效的求解方法成为问题的关键。 通过数学变换的方法可以降低方程(1-1)的矩阵阶数，或减小矩阵的带宽。如果能变换成对角矩阵，即带宽为1，将是最理想的结果，这样方程组中每一个方程都是独立的，我们只需求解一维问题即可。 在选择数学变换时需要首先注意的问题是数学变换不应使原系统的固有特性有所改变，也就是说，数学变换改变的仅仅是系统在空间域的描述方式，而不是系统的物理性质。 为此，我们假设对自然节点位移列矢量做如下变换： (6-1)其中为自然节点位移n阶列矢量，P为一个非奇异的矩阵，称为广义位移。也就是说，经过变换矩阵P，把原来位于由自然节点位移列矢量构成的n维空间中的M、D和K矩阵，变换到由广义位移列矢量所构成的m维空间中去。 我们把式(6-1)代入到式(1-1)中，并前乘，则得： (6-2)其中： (6-3) 从数学上看，可以找到许多不同的变换矩阵使原系统带宽减小或阶数降低。但从力学上看，一个有效的变换来自无阻尼的自由振动平衡方程： (6-4) 方程(6-4)的解可假定为下面的形式： (6-5)其中是n阶向量，t是时间变量，是时间常数，是表示列矢量的振动频率的常量。将式(6-5)代入到式(6-4)中，我们得到广义特征问题： (6-6) 求解式(6-6)中得到m个特征解：()。其中特征矢量是M规格化正交的，即： (6-7) 并且， (6-8) 在力学上，我们称为第i阶模态或振型，是相应的振动频率。 对于式(6-6)的特征矩阵，我们有： (6-9)其中： 另外有： (6-10) 所以，很明显特征矩阵就是式(6-1)的一个合适的变换矩阵，它已经可以把K和M变换成带宽最小的对角矩阵形式。 因此，我们有： (6-11) 将式(6-11)代入到(1-1)中，并利用式(6-10)，可得： (6-12)其中： (6-13)的初始条件可用式(6-11)得到： (6-14)在实际使用时，式(6-11)中的不是取全部特征矢量，而是m个部分特征矢量，则式(6-12)变成m阶方程组，使方程组阶数降低。 如果考虑无阻尼的情况下，那以经过变换的方程组(6-12)中的各方程是独立的，因此有限元方程就成为非耦合的，每一个方程都视为一个单自由度系统进求解，求解将变得非常容易。 对于无阻尼的式(6-12)中的每一个方程有： (6-15) 求解式(6-15)的方法有两种，一种是用直接积分法，另一种是用Duhamal积分法求得。 对于完整的响应，必须算出式(6-12)中的所有解，并将每一模态的响应叠加起来就很到有限元的节点位移，即： (6-16)式(6-16)称为模态叠加法。当我们所取的模态数m小于n时，提高了计算的效率。但如何选取m值使效率提高，又同时可保证足够的精度是一个相当复杂的问题。模态叠加法计算系统响应的过程是： (1) 求解方程(6-6)的特征值和特征矢量。(2) 求解不耦合的平衡方程组(6-15)。(3) 将每一个模态以式(6-19)的形式进行叠加。 当考虑阻尼时，模态叠加法的过程将变得复杂一些，也可以利用梯形法则。 通常，利用模态叠加法解非线性动力学问题仅适用于系统处于局部非线性的情况并且仅考虑少数几个模态。比如，局部接触振动问题和撞击问题等。7、精细积分法 钏万勰院士提出了精细积分算法在国际上享有一定的盛誉。其提出的基于类算法的精细积分法，所得几乎是计算机上的精确解，是一种优于传统差分方法的积分算法。这类方法为解决非线动力学问题提供了有效的途径。 实现精细算法高精度、高效率的内在机理和根本原因如下所述： (1)，指数矩阵的级数展开式绝对收敛； (2)初始级数展开式中的有效展开项总数能够通过递推算法以指数方式扩展； (3)新境有效展开项的系数能够通过递推算法以指数或拟指数保留基数M或递推阶数N的增大，精细算法的逼近误差上界以指数方式减小。 在计算机有效精度范围内，对于常系数线性微分方程，它的解与精确解完全一致。精细积分法同时具有可以用大步长、绝对稳定的优点。对于线性常系数刚性方程，任何时刻可以一次求