非线性动力学系统高性能计算方法

第1章绪论 1.1非线性动力学系统计算方法概述 工程中几乎所有的问题本质上都是非线性的，线性系统是在一定假设条件下对非线性系统的理想化近似。对于非线性动力学系统而言，一般由以时间为自变量的非线性微分方程描述。根据研究对象的自由度划分，非线性动力学系统可分为有限自由度系统（或集中质量系统），由一个或者一组非线性常微分方程所描述；无限自由度系统（或分布质量系统），由一个或一组非线性偏微分方程所描述。有限自由度系统只含有一个时间自变量，因此动力学方程表现为常微分方程的形式，研究对象一般用质点、质点系或刚体代替。例如，在轨道上运动的卫星与地球组成的系统，在不考虑其他星体引力作用时可以看成两个质点构成的质点系，此时卫星位置可用三个坐标完全确定，如果考虑卫星姿态，此时卫星看成具有六个自由度的刚体。但是不论将卫星看成质点还是刚体，本质上都是由非线性常微分方程所描述的有限自由度系统。无限自由度系统除了时间自变量之外还有其他自变量，对于力学系统来说，通常以空间坐标为其他自变量，因此动力学方程表现为偏微分方程的形式，研究对象表现为弹性板、壳、梁，或者流体等连续介质。例如，研究柔性航天器、飞机气动弹性变形等问题时，所面临的就是这类无限自由度系统。通常情况下，无限自由度系统可以通过空间离散的方法（如假设模态法）退化为有限自由度系统，得到退化后的只有一个时间自变量的微分方程，因此，只需要处理非线性微分方程即可。 本书主要讨论非线性常微分方程动力学系统的高性能解算方法。在实际问题中，很多动力学系统用以时间为自变量的常微分方程组x＝f（x，t）描述。其中，x是N维状态矢量，它是关于时间t的函数。根据系统中是否显含时间t，可分为自治系统和非自治系统。工程中绝大多数问题都是非线性的，即f是x的非线性函数。由于非线性的存在破坏了线性系统叠加性，方程组x＝f（x，t）精确解析解一般不可得。所以非线性动力学系统的研究主要依赖近似解法或数值解法。 对非线性动力学系统的求解需求由来已久，*早可以追溯到300多年前牛顿时代。1687年，牛顿发表了其*重要的著作《自然哲学的数学原理》，总结了近代天体力学和地面力学的成就，提出了力学的三大定律和万有引力定律，使经典力学成为一个完整的理论体系。牛顿在这本著作中提出了一个基础的天体力学问题： 如果两个物体以连心力互相作用，力的大小与它们的距离平方呈反比，那么两个物体会如何运动？这就是天体力学领域的“二体问题”。17世纪下半叶，牛顿和莱布尼茨发明了微积分，为动力学的数学化描述奠定了基础。二体问题的动力学方程可用经典牛顿定律得出，如图1－1所示，其中G为万有引力常数，r为两质点的距离。 图1－1二体问题及其数学描述 1710年，瑞士数学家约翰 伯努利首先解出了二体问题的解析解，得到了问题的答案“如果在一个质点上观察另一个质点，那么另一个质点的轨道一定是一条直线，或一条抛物线，或一个椭圆，或双曲线的一支”。二体问题也是极少数可以精确解析求解的非线性动力学方程的经典案例。 从历史发展的脉络来看，非线性动力学系统的求解方法经历了解析法、早期数值法、近似解析法、计算机辅助数值法四个时期。17世纪末，以严格的微分方程形式对动力学系统进行描述后，如何对微分方程进行求解是直接面临的难题，以牛顿、伯努利为代表的科学家利用变换与积分等手段对特殊类别的方程进行了解析求解。然而，绝大多数的非线性方程无法求解甚至不存在闭合解析解，这种情况下，18世纪上半叶，泰勒、欧拉等开展了早期的数值求解方法研究，如欧拉积分法就是建立在差分原理上逐步递推的数值方法，沿着这条思路发展出了RungeKutta法等耳熟能详的经典方法。受限于手推计算，早期数值计算的精度不高，解算过程极其烦琐。近似解析法出现在19世纪。众所周知，二体问题早在牛顿时代已被解决，根据牛顿的万有引力定律，学过高中物理的学生都不难列出三体问题的动力学方程，它是含有九个方程的微分方程组。但是，求解这个方程则是难上加难。在著名的三体问题牵引下，Poincare提出了摄动法等近似解析法。摄动法又称小参数展开法，利用小参数展开求解非线性方程的渐近近似解析解。近似解析法使用过程中需要大量的公式推导，往往只能进行前几阶计算，且对强非线性问题求解效果不佳。20世纪中叶，现代计算机技术的出现为非线性动力学系统的求解打开了一扇窗户，依托计算机辅助计算，数值计算方法得到了前所未有的蓬勃发展。 在飞机、导弹等国防工业需求牵引下，计算机辅助数值计算的研究热潮首先出现在固体、流体等连续介质力学领域，求解对象在数学本质上为偏微分方程。求解思路是通过有限差分法、有限元法等将空间坐标离散，得到相应的代数方程组（静力学问题）或者常微分方程组（动力学问题）。有限差分法是*早发展的空间离散算法，基于广义加权残余法的思想，通过离散化空间域中残差为零的不同处理形式，发展了一系列数值算法，其中主要包括： ① Galerkin法（Galerkin method），试函数与权函数可以相同也可以不同，但二者均为任意阶连续的全局函数，当试函数与权函数不同时，即为Petrov－Galerkin法（Petrov－Galerkin method）；② 配点法（collocation method），试函数采用任意阶连续的局部或全局函数，权函数采用狄拉克δ函数；③ 有限体积法（finite volume method），试函数采用任意阶连续的局部或全局函数，权函数采用局部Heaviside阶跃函数；④ 初等Galerkin有限元法（primal Galerkin finite element method），试函数与权函数均为相同的全局任意阶连续函数；⑤ 无网格局部Petrov－Galerkin方法（meshless local Petrov Galerkin method），Satya N. Atluri与合作者自1998年以来发展了一系列此类方法，其中试函数采用单位划分、移动*小二乘或径向基函数等无网格函数，权函数则采用狄拉克δ函数、Heaviside函数、径向基函数等。将偏微分方程通过上述任何一种方法进行空间坐标离散，可以得到非线性常微分方程，即非线性动力学系统，即本书主要研究内容。 本书中前5章所涉及的时域配点法、谐波平衡法和高维谐波平衡法均为全局法，即求解动力学方程时，以整个时间域为对象进行求解，不进行子域分割。全局法适用于求解具有周期特性的非线性动力学系统。当非线性动力学系统的解具有周期性时，传统求解方法是摄动法和谐波平衡法。摄动法受限于复杂符号推导和弱非线性，不利于高精度计算。谐波平衡法本质上是时域的Galerkin法。谐波平衡法是求解非线性动力学系统周期解的强有力的工具。谐波平衡法的原理是，首先假设系统的近似解为Fourier级数形式，然后将其代入动力学方程令各次谐波自相平衡得到以Fourier系数为变量的非线性代数方程组。该方法不受非线性强弱的限制，且对于不太复杂的周期响应一般只需几个谐波即可得到高精度的解。但是，谐波平衡法在对非线性项进行Fourier级数展开过程中涉及烦冗的公式推导，并且随着谐波个数的增多会变得越来越复杂。为了克服这个缺点，2002年美国杜克大学的Hall等[1]提出了高维谐波平衡法，该方法随后被广泛应用到各种动力学问题的求解中。高维谐波平衡法避免了传统谐波平衡法的复杂公式推导，但是产生了多余的非物理解。之后，Liu等[2]指出，高维谐波平衡法的非物理“假解”源于高次谐波对低次谐波的混淆（即著名的混淆现象）。但是，高次谐波到底以何种规律混入低次谐波这一问题没有得到解答。本书将深入研究高维谐波平衡法、谐波平衡法等非线性动力学系统的半解析求解方法，在2.3节对高维谐波平衡法与谐波平衡法的关系进行证明，澄清了困扰学术界十余年的“高维谐波平衡法是谐波平衡法变种”的误解。在3.4节给出高维谐波平衡法产生混淆现象的机理，解释“假解”产生的本质原因（图12）。 图1－2高维谐波平衡法的混淆机理示意图 本书第6～8章介绍的变分迭代积分法、改进RungeKutta法等数值积分法，均为局部法。局部法即求解问题时其时间域先被离散为有限个子域，然后求解；每个子域上又可以采用各种全局方法进行求解，例如，本书第7章提出的变分迭代配点法就利用变分迭代和时域配点的思想进行快速计算。可见，全局法和局部法是不可分割的，从其发展历程来看，二者是相互促进相互依存的关系（图13）。借助现代化计算机，数值积分法成为研究非线性动力系统响应*为直接的方法。在文献中存在着数以百计的数值积分法，其中大部分是有限差分法，如RungeKutta法、HilbertHughesTaylor法及Newmark法。这些方法直接使用了微分的定义以及Taylor级数理论对原始的常微分方程离散化，并在时间间隔较小的离散点上给出数值估计。基于每个时间间隔中的估计方法是前向估计还是后向估计，一般可以将其分为显式方法和隐式方法。一般来说，隐式方法比显式方法更加稳定，但是前者一般需要求解相关的非线性代数方程组。 图1－3非线性系统数值积分方法发展历程[38] 1.2航空航天领域的典型问题 航空航天领域存在着大量非线性动力学问题，结合作者十余年的研究经历，本书主要围绕四个典型问题进行计算方法的理论研究和应用分析。这四个问题分别为经典Duffing方程、非线性气动弹性问题、轨道递推与转移问题、近地航天器的相对运动。其中轨道转移问题是两点边值问题，其他为初值问题。 Duffing方程是工程中*常见的非线性方程之一。很多航空航天领域中的问题都可以直接用Duffing方程描述或者处理后退化为Duffing方程。例如，使用Galerkin法将超声速流中板的颤振模型进行空间离散，得到的模态方程就是耦合Duffing方程组。对于二元机翼模型，它经过积分变换后可以转化为含有三次非线性项的微分方程组，如果考虑二元机翼受简谐作用外力的话，该系统中的非线性方程将表现为Duffing方程形式，如果不考虑二元机翼系统的气动力，那么该模型振动问题将退化为两自由度耦合Duffing振子模型，鉴于Duffing方程在非线性动力学领域的广泛代表性，本书将其作为计算方法的验证原型反复使用。此外，非线性气动弹性问题也是航空航天领域重点关注的问题之一，由于机翼气动弹性引起的颤振、抖振等效应是飞行器设计过程中必须考虑的问题之一，对气动弹性问题的研究及高精度的气动弹性仿真计算结果将对飞行器设计产生重要影响。本书将围绕亚声速二元机翼动力学模型开展相关的计算方法研究，并对机翼颤振问题进行探究。卫星轨道递推与轨道转移问题是航天动力学与控制领域的基础问题。航天器正常工作的基础在于其始终保持在目标工作轨道位置，这需要通过轨道递推方程不断计算航天器在当前以及未来时刻所应到达的目标位置，并通过计算当前实际位置与目标位置的距离偏差进行调整。轨道转移计算可以看成一类受到两点边值约束轨道递推问题，在转移轨道过程中受到与轨道递推问题相同的二体运动方程约束，在诸如火星探测等远距离深空探测任务中，航天器轨道运算结果的精度对任务成败起到决定性作用。而星际探测航天器往往面临着星载计算机计算能力低下的硬件限制，因此，研究可在计算能力低下的星载计算机中实现的高精度轨道递推计算的数值算法，是星际探测任务中亟须解决的重要问题。此外，近年来兴起的空间卫星编队技术以及空间在轨服务技术对近地航天器相对运动问题的高效高精度求解提出了迫切需要，适用于星载计算机快速求解复杂近地航天器相对运动方程的高效高精度数值计算方法研究成为相关技术发展的核心问题之一。 1.2.1经典Duffing方程 Duffing方程是一种描述强迫振动的微分方程模型。该模型是非线性理论中常用的一种典型非线性振动系统模型，这一系统虽然形式简单，但却具有高度的代表性。事实上，工程实际中诸多非线性振动问题的数学模型都可以归结为