非线性体系随机响应计算的确定性积分方法
Classified Index: TU393.4
U.D.C.:681.014
Dissertation for the Master Degree in Engineering
DETERMINISTIC INTEGRATION
ALGORITHMS FOR STOCHASTIC
RESPONSE COMPUTATIONS OF
NON-LINEAR SYSTEMS
Candidate: Yue Juan
Supervisor: Prof. Wu Zhifeng
Academic Degree Applied for: Master of Engineering Speciality: Solid Mechanics
School of Civil Engineering Affiliation:
June,2006 Date of Defence:
Harbin Institute of TechnologyDegree-Conferring-Institution:
摘 要
摘 要
在随机振动领域中，随机响应主要由均值和协方差来描述。随着有限元 理论的发展，直接积分计算方法在随机响应
分析
定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析
中，特别是对于大型多自由 度体系和自由度数较大的有限元体系，已经有很多应用。如中心差分法、 Houbolt 法、 Newmark 法等。然而，这些方法在处理多自由度的有限元模 型，特别是自由度数 n 很大时，计算量大，而且需要很大的存储空间。
本课
题
快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题
研究的是一种结合 Karhunen-Loéve(K-L)展开的求解协方差矩阵 的确定性积分方法。引入的 K-L 展开向量可以用任何一种确定性的逐步积 分方法进行积分。对于大型有限元系统，协方差矩阵的维数 , 2n, 其中 n 为 系统自由度数，用远小于 n 的 m 个 K-L 向量就可以准确地表达协方差矩 阵。当 n 很大而 m 很小时，计算优势就会很明显。这样，这种方法就特别 适合求解大型有限元系统的协方差矩阵。
结合协方差矩阵 K-L 分解的确定性积分方法主要是利用 K-L 展开向量 的一些优势，将协方差矩阵用其特征值的 K-L 展开替代，然后用任意的确 定性逐步积分方法进行积分。这种方法可以解决有上千自由度的大型有限元 模型的随机响应问题，使有限元模型的计算在计算效率上得到很大的提高， 并且准确适用。
本文首先在第二章中介绍 K-L 展开在随机过程中的应用，将随机过程， 随机激励，随机响应用 K-L 展开表示，然后在第三章中介绍利用 K-L 展开 进行随机响应计算的方法，在第四章中，主要讨论 K-L 展开结合等效线性 化求解非线性系统二阶统计矩的计算过程，最后在第五章中，通过算例将这 种方法应用到实际结构中，并将计算结果进行系统分析。
关键词 随机响应;非线性体系;协方差矩阵;K-L 展开;确定性积分方法
-I-
哈尔滨工业大学工学硕士学位论文
Abstract
In stochastic analysis, the response is described by the mean and covariance function. Direct integration algorithms for the stochastic response computation, published since the development of the finite element methods, are applicable for MDOF-systems of FE-methods with few hundred degrees of freedom (DOF). First the stochastic central difference method has been suggested, followed by the stochastic Houbolt method and the stochastic Newmark algorithm. However, all these algorithms cannot be applied efficiently for FE-models with thousands DOFs, particularly for cases in which n is considerably larger than a few hundred, because of the huge number of required operations and also due to storage requirements.
The method introduced in this paper employs deterministic integration algorithms together with the Karhunen-Lo é ve (K-L) representation for the covariance matrix. The deterministic vectors for the K-L representation can be integrated by any available deterministic step-by-step procedure. A relatively small number m<
标准
excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载
化 ............................................................20
2.7 本章小结 ..................................................................................................21 第 3 章 李雅普诺夫微分方程的确定性积分求解 ............................................22
3.1 引言 ..........................................................................................................22
-IV-
目录
3.2 方差矩阵的李雅普诺夫微分方程 ...........................................................22
3.3 李雅普诺夫微分方程的求解 ...................................................................24
3.3.1 逐步积分法计算过程........................................................................24
3.3.2 基于 K-L 展开的逐步积分法求解李雅普诺夫微分方程 ................27
3.4 子空间分解 ..............................................................................................29 3.5 本章小结 ..................................................................................................31 第 4 章 非线性体系随机响应计算的求解方法 ................................................32
4.1 引言 ..........................................................................................................32 4.2 非线性系统运动方程 ..............................................................................32 4.3 等效线性化方法在非线性系统中的应用 ...............................................34
4.4 滞迟系统在平稳随机激励下的随机响应分析 .......................................36
4.5 本章小结 ..................................................................................................37 第 5 章 算例分析 ...............................................................................................39 5.1 引言 ..........................................................................................................39 5.2 单自由度体系的随机响应分析 ...............................................................39
5.3 多自由度体系的随机响应分析 ...............................................................44
5.4 本章小结 ..................................................................................................49 结 论 .................................................................................................................50 参考文献: .........................................................................................................51 致谢 .....................................................................................................................55
-V-
第 1 章 绪论
第1章 绪论
1.1 课题背景
1.1.1 随机振动理论的发展
随着科学技术的发展，随机振动理论已经越来越成熟，传统的随机振动 理论一般是研究确定性结构在随机激励下的响应计算问题;在自然界和
工程
路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理
中，存在一大类诱发结构系统振动的外因，如大气洄流，地面强风，地震动 等等，他们的一个共同特点就是随机性，即不能用确定性的时间与空间坐标 的函数来描述，这些都称为随机振源，由这些因素引起的结构系统的振动称 为随机振动[1]。
在结构工程领域主要研究的是地震引起的随机振动[2]。作用在结构体系 上的动荷载的发生机制可能是各种各样的，其中包括风或波的作用以及车辆 运动等，但对于结构物来说，地震无疑是最重要的动力输入形式。地震问题 之所以重要，部分原因是由于大地震在人口稠密地区所造成的可怕后果[3]， 因此，结构要有能抵抗强烈地面运动所产生的动荷载的能力。关于地震动的 随机振动理论在土木工程领域得到了很大的发展。
结构地震反应分析方法，可以分为两大类[4]:确定性分析和非确定性分 析。确定性分析是指结构在确定性地震作用下的反应分析，它的运动时间历 程可以用确定性函数来描述。而非确定性分析在地震工程领域里一般是指结 构在随机振动作用下的反应分析。众所周知，地震动是非常复杂的现象，具 有很强的随机性，即使在同一地点上，遭受相同震级且震源也相同的两次地 震，结果也不会完全一样[5]，也就是说具有非重复性，因此，随机振动理论 的分析方法可以较好的反映这些特性。
结构随机振动是上世纪 60 年代以后发展起来的一门新的学科[6]，线性 时不变系统随机振动的一般理论已经比较成熟，但非线性系统的随机振动问 题还没有得到满意的解决[7]，仍然没有一种成熟的一般解法，虽然研究人员 提出和发展了很多方法，但这些方法都各有利弊，非线性随机振动在某种意 义上说仍然有很大的研究空间。
-1-
哈尔滨工业大学工学硕士学位论文
1.1.2 非线性随机振动理论的研究
非线性随机振动理论是随机振动的一门重要的分支，在很多工程问题 中，如足以引起严重破坏的地震动作用下的建筑物反应，近海结构物在波浪 力作用下的动力反应等，由于体系在动力反应过程中，动力特性是变化的， 因此，适用于线性体系的叠加原理就不能应用了，为此需要发展适用于非线 性体系的分析方法。
严格的说，在土木、机械、航空、航天及海洋工程等几乎所有的工程领 域中，非线性现象普遍存在，结构体系的振动总不同程度地满足某种非线 性，只是在小变形的微幅振动下大多数体系的非线性特征不明显，我们可以 用较好的线性模型来模拟，当实际结构的非线性比较弱，或非线性对所研究 的问题不产生显著的影响时，线性系统作为非线性系统的近似模型而得到广 泛的应用，解决了大量工程实际问题。按照线性理论来模拟所得的结果，往 往是实际系统随机响应量的一次近似，这种近似在许多情况下是满意的，但 也有很多场合，问题的线性处理并不能给出满意的结果。当实际结构的非线 性比较强或存在跳跃、分叉、参激等本质非线性时，用线性系统分析难以对 现象做出正确的解释，必须依赖于非线性分析而使问题得到合理的解决。
在十九世纪末，Poincare 首先从几何和拓扑观点对天体力学问题进行定 性研究[8]，他的工作和思想对非线性科学的发展有深远的影响。随着近代科 学技术的迅速发展，许多非线性力学问题急需解决，大型高速电子计算机和 有效的大规模算法的出现，使得人们能够对非线性问题进行大量的数值计算 和仿真，揭示了极其丰富的非线性动力学现象。各种现代数学理论的发展， 为非线性的研究提供了强有力的工具。于是，非线性动力学无论从广度到深 度都以空前的速度向前发展，成为当前非常活跃的力学分支。非线性动力学 研究非线性力学系统各种运动状态的定性和定量变化规律，当前非线性动力 学的中心课题包括分叉，混沌，分形，孤立子，拟序结构，图像生成和演化 等[9]。概括地说，非线性问题涉及十分广泛的范围，有很大的研究空间。我 们知道随机振动是普遍存在的振动形式，非线性系统和结构又是普遍存在的 结构形式，非线性系统的随机响应计算就有着非常广泛的研究领域。
1.1.3 非线性随机振动研究过程中的主要难点
对于大变形振动或本身含有非线性元件的体系振动，我们必须用非线性
-2-
第 1 章 绪论
的模型来描述。非线性振动的反应分析比线性的情况要复杂得多，根本原因 是叠加原理不适用。特别是在计算非线性多自由度体系以及有几千个自由度 的有限元模型时，在计算上存在很大的困难。首先，随机激励的幅值往往没 有物理上限，大幅的响应总是可能的，虽然这种大幅响应出现的概率很小， 但它关系到结构的破坏，而非线性效应正是这些大幅响应中起着重要甚至决 定作用的。其次，非线性系统中，存在着许多线性系统中不会出现的现象， 如跳跃，自由振动，混沌等，因此，非线性随机响应计算在实际中会有很大 的困难出现。
非线性系统的数学模型由非线性微分方程来表示，对这些方程的求解存 在着许多困难，由于叠加原理不再适合非线性系统，在线性系统中经常使用 的十分有效的方法，如振型分解法等都不能应用在非线性系统中了。另外， 在 Gauss 随机激励下，非线性系统的响应一般不再是 Gauss 过程，由相应的 前二阶矩不能得到相应的全部统计特性。因而要研究相应的高阶矩，或直接 寻求相应的概率密度函数，这往往是非常困难的。
事实上，正是由于非线性系统在计算上存在这些困难，近二十年来，在 众多学者的努力下，已经发展了许多计算方法，其中一类是扩散过程理论方 法，即 FPK 方程法，体系的状态反应是马尔科夫过程矢量时，由于反应的 转移概率密度满足 FPK 方程，因此在一定初始条件下，我们可以通过求解 FPK 方程来得到反应的概率密度。另一类是从确定性非线性振动理论方法 推广来的，包括摄动法，等效线性化法，分段等效线性化方法，各种截断 法，拟静态法，及数值模拟法等。这些方法或是因为求解范围小，或是因其 计算繁琐而大多难以满足工程需要。
在非线性系统中，位移和速度反应的均值和协方差矩阵的计算是难点， 由于协方差矩阵的形状庞大，在存储和计算上都存在很大的问题 [10]。尤其 在处理大型多自由度系统，以及有上千自由度的有限元系统时，就会出现更 多问题。因为要计算系统的协方差矩阵 , ，而协方差矩阵 , 的维数与体系的 自由度数有关，一般来说协方差矩阵的维数不小于自由度数的二倍。当自由 度数很大时，计算很难得到有效的实施，将耗费很大的存储空间和计算机 时。因此需要一种合理有效的解决办法。如何解决这个问题，在学术上引起 很多学者的关注。
-3-
哈尔滨工业大学工学硕士学位论文
1.2 国内外在非线性随机振动方向的研究现状及发展综述
1.2.1 求解非线性问题的一般方法
非线性确定性振动的反应分析，实质上就是具体的求非线性微分方程的 解，由于非线性微分方程的求解很困难，目前非线性振动的反应分析还没有 一种能够普遍适用的方法，只能采用某些特殊的方法来尽可能揭示体系的运 动状态。早在四十几年前，非线性随机振动问题的研究就受到了重视，时至 今日，已发展了很多方法，例如:摄动法，统计矩截断法，随机因子法， FPK 方程法等，下面介绍将这些国内外经常使用的解决非线性问题的方法 做以简单的介绍。
1.2.1.1 摄动法
当结构的运动方程中非线性部分较弱时，用摄动法(erturbation Method)
解随机输入下非线性系统的响应可以得到很有成效的结果 [11]。摄动法认为 弱非线性系统的解接近于相应的线性解，因而将解过程展开成小参数 , 的幂 级数，并且有系统的运动微分方程中 , 的同次幂系数相等的条件，通过求解 一系列线性微分方程来确定幂级数中各项的系数，从而得到非线性微分方程 的近似解。因此，摄动法就是把一个求解非线性随机振动的问题转化为求解 一组线性随机矩阵微分方程的问题。
应当指出的是，从计算的观点来说，摄动法中涉及的高阶项的统计性质 是很难求得的，另外，对于幂级数的展开式，其收敛性无论在均方意义下还 是在其它意义下均无满意的证明，所以，要使摄动解有效，需要假设解过程 X (t) 的每个样本函数都能用 , 的收敛幂级数表示。通常，摄动法不但要求 , 充分小，而且要求干扰的振幅也要充分小。因此，摄动法虽然减少了计算工 作量，但对于结构中各种参数的随机性均是以一个小参数(即摄动量)笼统 的予以描述，由于其方法本身的局限性，
决定
郑伟家庭教育讲座全集个人独资股东决定成立安全领导小组关于成立临时党支部关于注销分公司决定
了它将无法反映出结构中某一 参数的随机性对结构动力参数的影响。
1.2.1.2 统计矩截断法
统计矩截断法是一类求解非线性系统反应矩的方法 [12]，比较常应用的 方法有四种:高斯截断方法，累积量截断法，中心矩截断法，和线性表出 法。在工程实际中，最关心的是反应的一阶矩和二阶矩，因此，我们可以用 统计矩截断法将反应的二阶矩截断，将三阶距和三阶以上的矩用前两阶矩表 示出来。同样的，在第 N 阶矩截断，也不能认为 N 阶以上的矩等于零，必
-4-
第 1 章 绪论
须通过某种方法，将 N 阶矩及 N 阶以上的矩用前 N 阶矩近似的表示出来， 然后利用方程求解。这就是统计矩截断法的基本思想。统计矩截断法都基于 响应为高斯分布的假设，可适用于单个或多个自由度非线性系统受平稳或非 平稳随机激励的情况。在 Gauss 随机激励下，非线性系统的响应一般不再是 Gauss 过程，由相应的前二阶矩不能得到相应的全部统计特性。矩函数截断 法的计算过程较为复杂，计算量大，很难在工程应用中推广。 1.2.1.3 随机因子法
对于摄动法中存在的问题，随机因子法可以很好的解决，该方法的基本 思想是通过将随机变量分解成一个随机因子与一个确定性量的乘积，并设法 将结构的动力特性分析结果表达成结构参数的随机因子与动力特性确定性量 乘积的函数 [13]，再利用求解随机变量函数的数字特征的矩法或代数综合 法，便可以求出结构动力特性随机变量的均值和方差等统计量，然而该方法 主要针对线性结构能提出随机因子的问题，在非线性问题上的解法需进一步 完善。
1.2.1.4 FPK 方程法
在众多求解非线性系统随机响应的方法中，FPK 方程法吸引了大批学 者的关注。从理论上说，它是非线性随机振动分析最严密的方法，可以得到 特殊问题的解析解。对它的研究也大大推动了非线性随机振动科学的发展， 以及对这类振动现象本质的认识[14]。
对于受白噪声激励，而响应为 Markov 过程的系统，用 FPK 方程法可以 解决平稳与非平稳、弱非线性系统。由于在均方意义下白噪声的积分无意 义，所以只有将振动微分方程化为伊藤型方程时，其响应随机过程才具有良 好的 Markov 随机过程，响应就可以完全由它的初始条件和转移概率密度函 数决定，然而，当外干扰力是平稳的正态的白噪声过程时，响应的转移概率 密度函数可以从求解 FPK 方程的方法中得到。一经响应的条件概率密度函 数确定之后，那么响应的时域及频域的各种信息就可以完全解出了。
可以看出， FPK 方程法能够解决的问题不多，而且附加了苛刻的条 件，其研究基本上局限在单自由度体系 [15]。非线性体系的随机振动的 FPK 方程法就是要通过求解体系的状态方程或扩充的状态矢量的 FPK 方程来求 得状态矢量，或扩充了的状态矢量的联合概率密度。因此，从理论上来说， FPK 方程法是非线性随机振动分析的最严密最完美的方法。然而，遗憾的 是，FPK 方程的精确解对于非线性体系的非平稳反应只是在少数一阶体系 的情况下才能得到，即使对于平稳反应的情况，也只有少数几类特殊的单自
-5-
哈尔滨工业大学工学硕士学位论文
由度和多自由度体系才能得到。因此，对于比较一般的非线性随机振动问 题，还只能用近似的分析方法。
1.2.2 大型多自由度系统随机响应计算的有效方法
1.2.2.1 Monte-Carlo 法
Monte-Carlo 方法是一种随机抽样的方法又称为概率模拟法，统计实验 法，它是求解数学、物理、及工程技术问题近似解的数值方法，属于实验数 学的一个分支[16]，它的应用非常广泛。Monte-Carlo 方法在非线性随机计算 中的应用越来越被重视。在应用时，首先需要建立一个概率模型，把需要求 解的问题与概率模型联系起来，然后通过随机抽样试验，得到某些统计特征 量作为待求问题的近似解，这就是 Monte-Carlo 方法的基本应用框架。 Monte-Carlo 法研究随机振动问题的基本思想就是概率论中的基于大数法则 的逼近原理，即系统响应的统计特性可以从大量响应样本中近似获取，其近 似程度随响应样本的增加而提高。在 Monte-Carlo 法中，采用简单抽样方式 进行随机变量的数值模拟，它所抽取的子样满足同分布性质(例如正态分 布)的独立随机变量，当抽取样本个数足够多时，样本均值将以概率 1 收敛 于分布均值，而事件 A 出现的频率 n/N 则以概率 1 收敛于事件 A 出现的概 率。这样就保证了 Monte-Carlo 方法的概率收敛性。利用 Monte-Carlo 方法 来进行随机响应分析，主要包括以下三个步骤[17]:
(1) 根据激励过程或场的统计特性产生激励样本;
(2) 对每个激励样本，求解运动微分方程，产生相应的数值样本;
(3) 在大量的响应样本中获取所需要的响应统计量。
目前，Monte-Carlo 方法已经成为研究随机振动，尤其是非线性随机振 动的重要方法之一。原因有两个:首先，Monte-Carlo 方法原则上适用任何 系统与任何激励，尤其是对于一些用现有方法很难处理的复杂问题，如强非 线性，强参数，严重的材料不均匀等。对于这些问题它往往是最可行的解决 办法。其次，该方法是检验各种近似方法的适用性与精度的有力工具。
在足够的样本情况下，该方法能够获得响应的任意高阶矩，对随机有限 元模型建立的新方法，一般是利用 Monte-Carlo 随机抽样检验其有效性[18]。 现今成熟的商用有限元分析软件为 Monte-Carlo 的实现提供了十分方便的途 径。Monte-Carlo 法具有普遍性，适用于各种问题，但是，Monte-Carlo 方法 的缺点如同他的优点一样明显，那就是计算的耗时性，往往需要海量计算， 但存在数据收敛判据方面的困难，对于具有一定规模的非线性结构计算问
-6-
第 1 章 绪论
题，单个样本的计算已经十分耗时，若要进行足够样本量的计算，在计算时 间上是绝不可行的。用该方法即使研究一个很简单的问题，也往往要耗费大 量机时，要使结果增加一位有效数字，计算时间要增加上百倍。这些都是在 工程上应用的主要障碍。尽管如此，对于复杂非线性问题，Monte-Carlo 方 法在计算上仍然有着不可替代的地位。Monte-Carlo 方法在非线性体系计算 中的核心问题就是如何减小计算时间，降低计算成本。其研究主要在两个方 面进行，一是基于误差分析的抽样法，二是并行算法。
在国外，随机动力响应计算在非线性体系中的应用，主要使用 Monte Carlo 模拟，或等效线性化的方法[19]。当自由度数较小，一般只有几十个自 由度的多自由度体系可以用简化的线性模型来模拟，但当有限元模型含有几 千自由度时就不能用线性化的结构进行计算了。在求解大型非线性结构时， 等效线性化法和 Monte-Carlo 法比较常用，与其他方法相比，等效线性化法 和 Monte-Carlo 法概念清晰，应用简单。其中，等效线性化法是对响应分布 的整体坐标下上的计算，而 Monte-Carlo 法则是在局部坐标下的计算。等效 线性化法计算的是随机响应的前两阶矩，Monte-Carlo 法可以得到随机动力 响应的更多信息，如动力可靠度等。
1.2.2.2 等效线性化法
等效线性化方法，又称为统计线性化方法或随机线性化方法，是工程中 应用的最广泛的一种计算非线性系统随机响应的近似解析法，其基本思想 是，用一个具有精确解的线性系统代替给定的非线性系统，使两方程之差在 统计意义上最小[20-25]。该方法最早由 Kazakov，Booton，和 Caughey 各自独 立提出的，当时的主要应用是在控制理论和自动化领域中。Booton 的理论 影响较小，Kazakov 的线性化理论被苏联的许多学者引用和发展。两位作者 都将方法应用在正态白噪声激励下的非线性系统中，在这种情况下，有一些 困难出现，为了计算线性化参数，要解非线性系统的可靠度密度函数，一般 来说，这个函数是未知的，要用线性系统的可靠度密度函数来代替非线性系 统的密度函数。这些计算要在一些假设下完成，即在进行线性化系统的可靠 度计算时，非线性参数的期望不变，这种替代和假设在计算非线性系统时非 常适用。Caughey 在总结 Krylov 和 Bogoliubov 的等效线性化理论时，将这 个理论进行了改进和发展。在线性化参数随时间变化的非稳态问题中，通过 在时间轴上的离散，等效线性化方法同样得到实现。这时先要解决期望问 题，高斯激励下的期望可以通过解均值和协方差矩阵的微分方程得到。对于 非高斯激励，同样可以计算，只是计算过程比较复杂。
-7-
哈尔滨工业大学工学硕士学位论文
另外，朱位秋等和赵雷等对等效线性化作了推广，发展了等效非线性系 统方法，该法的主要思想是以某个具有精确解的非线性系统代替给定的非线 性系统，使得两方程之差在统计意义上或平均意义上最小，它在一定程度上 能克服等效线性化法的局限性，但由于具有精确解的非线性系统的有限性， 使该法的应用也就非常有限。
1.2.2.3 复模态方法
目前，用等效线性化方法进行非线性结构体系的随机响应分析，还是有 一定限制的，一般限制在单自由度体系，或含有较少自由度的多自由度系 统，在这种情况下，有一种复模态分析方法比较可行 [26]，下面将其予以介 绍。
对于复模态分析方法，在计算精度的要求很高，积分工作量也可能是很 大时，是较为合适的积分方法。对线性系统和线性子空间系统的模型分析就 是这样一种适合的方法，它应用的范围比直接积分法要广泛一些，因为它的 积分是在转换的模态空间中进行的，计算与自由度数无关。对运动方程进行 某种数学变换，使基底降低，或使其矩阵的带宽变小，再进行求解计算。在 非线性体系中，由于刚度矩阵是变化的，所以在每一个构形下的系统的特征 值和特征向量都是不相同的，因此对于每一个构形都有不同的特征区间，计 算起来很复杂。另外，这种方法首先需要解决复杂特征值问题，如果不需要 频繁计算特征值，那么这种方法可以在大型有限元系统中应用 [27]。对于非 线性系统，在零均值问题中，模型特性不会发生改变，这种方法应用起来比 较方便，但在非零均值问题中，模型特征随着时间发生改变，这样就需要频 繁解决特征值问题，计算量相当大，对于大型有限元模型的计算不适用。复 模态方法使用的是状态空间的展开，利用模型转换使矩阵规模降低，以其传 统的结构特征值进行计算。对于大型有限元系统，因为特征值矩阵形状庞 大，复模态方法不太适用。当特征值问题的解法被改进了，复模态方法在非 线性体系中的应用才可行。但是，对于非零均值问题，模型特征在每一个时 间步长上都要改变，其应用也是很难进行的。
1.2.3 数值积分方法在随机响应计算中的应用
对于非线性体系的动力分析，最有效的方法是数值的逐步积分法，这种 方法把反应的时程划分为短的，相等的(也可以不相等)时段，对每一个时 段按照线性体系来计算其反应 [28-30]。这个线性体系的特性是时段开始时刻
-8-
第 1 章 绪论
限定的特性，时段结束时的特性按照那时体系的变形和应力状态来修正。这 样，非线性分析就近似为一系列依次变化的线性体系的分析。
逐步积分法特别适用于线性体系，此时因为不需要每步都需要修正结构 特性，使计算过程简化，当要计算受短时脉冲荷载作用的结构反应，或其它 可能激发起多个振型而只需计算很短时间的反应时，这种不需要求解振型和 频率的直接积分法就非常方便了 [31]。对于逐步积分法，需要计算协方差矩 阵，协方差矩阵的维数大于等于自由度数的二倍，计算起来很费时，而且存 储量大。
在直接积分法中有两种基本方法，一是显式的直接积分法，二是隐式的 直接积分法。在显式直接积分法中，动力响应的诸量直接用已求得的位移、 速度、与加速度的数值表示出;而在隐式直接积分方法中，瞬时的差分方程 与运动方程相耦合，需要求解这些方程组才能计算出各瞬时的位移量。
在动力分析中，有三种常用的显式方法，它们是:中心差分法，具有梯 形法则的二循环迭代法和四阶 Runge-Kutta 法。也有三种常用的隐式方法，
, 法，和 Houbolt 法等。 即:Newmark 法，Wilson-
直接积分法的主要特点就是不对系统的动力学方程进行坐标变换，而直 接求其在离散时间域上的数值解，使用直接积分法时，须对时间域进行离 散，并在每一段时间域内给出系统的坐标，速度与加速度的两个近似关系式 通常是一种差分方程，将这个差分方程代入系统动力学方程中，由此可得到 由系统在过去时刻的诸运动量求得下一时刻的诸运动量的递推关系式。直接 积分法是目前在工程上经常使用的方法，同时，各种各样的改进的或新的直 接积分法还在不断出现。
1.2.4 求解大型多自由度非线性系统随机响应时存在的困难
应该指出的是，虽然目前已有许多预测非线性系统随机响应的方法，但 没有一个是令人十分满意的，尤其是对于多自由度的非线性系统，各种方法 的计算都有很多问题，因此非线性系统随机振动还有待于进一步研究。用逐 步积分法计算响应的协方差矩阵时，协方差矩阵的维数有时会相当大，计算 起来很费时，而且需要很大的存储量。
随着有限元理论的发展，直接积分计算方法在随机反应分析中，特别是 对于大型多自由度体系和自由度数较大的有限元体系，应用得越来越多。如 中心差分法，Newmark 法，Houbolt 法，Runge-Kutta 法等等[32]，直接积分
-9-
哈尔滨工业大学工学硕士学位论文
法的主要特点就是直接求出系统运动方程在离散时间域上的数值解，在随机 响应计算中则主要计算系统运动方程的二阶矩，即系统的协方差矩阵。但是 在处理大型多自由度系统，以及有上千自由度的有限元系统时，就会出现很 多问题。因为要计算系统的协方差矩阵 , ，而协方差矩阵 , 的维数与体系的 自由度数有关，一般来说协方差矩阵的维数不小于自由度数的二倍 [33]。当 自由度数很大时，计算很难得到有效的实施，将耗费很大的存储空间和计算 机时。大型有限元系统有其确定性的动力分析过程，对于有限元模型，要在 频域内获得这样的稳态解是很困难的，需要很大的工作量。因此需要一种合 理有效的解决办法。
1.3 本文的主要研究工作
本课题研究的是一种基于 K-L 展开的求解非线性体系随机响应的确定性 积分方法。这种方法可以避免大量存储，节省计算时间。在这个方法里面主 要利用 K-L 展开，将协方差矩阵用其特征值的 K-L 展开向量替代，然后用 比较常见的确定性逐步积分方法进行积分。因为协方差矩阵的维数大于等于 系统自由度数的二倍，所以多自由度体系的维数较大，而替代时所用的 K-L 向量个数远小于协方差矩阵的维数，所以这种方法主要是解决多自由度和有 限元系统的随机响应计算问题。当系统自由度数较大，而 K-L 向量数相对 较小时，会使其计算优势更加明显，使有限元模型的计算在效率上得到很大 的提高，并且准确适用。
本文首先将 K-L 展开应用到连续随机过程的离散化上。用 K-L 向量和 的形式表示随机过程，随机激励，和随机响应，从而将连续的随机过程离散 化;通过截断较小数量级的 K-L 向量，使实际参与计算的 K-L 向量数减 小，达到使计算简便的目的;然后将振幅和频率随时间而发生变化的非平稳 随机激励，通过子空间分解转化到状态向量中，使新状态空间下的随机响应 通过基于 K-L 展开的确定性积分方法进行求解;将随机响应协方差矩阵所 满足的李雅普诺夫微分方程用基于 K-L 展开的确定性积分方法求解，对求 解过程进行详细论述;通过基于 K-L 展开的确定性积分方法结合等效线性 化方法来求解非线性系统随机响应。
最后，通过算例将上述方法在实际结构中进行计算分析。通过单自由度 体系随机响应计算的具体算例，将用基于 K-L 展开的确定性积分方法计算 随机响应协方差的计算结果与通过杜哈美积分求得的结果进行对比，得出结
- 10 -
第 1 章 绪论
论。进一步通过五自由度的非线性体系的具体算例，结合等效线性化方法， 对非线性体系的随机响应进行计算，将计算结果与相应的线性体系的随机响 应进行对比，得出结论。
为提高计算效率，本文使用了确定性的参数矩阵，非线性存储力都在运 动方程的荷载部分里考虑，非线性元素被等效线性化，其在时间上的改变是 通过在每一个时间步长上的独立积分得到的，方程之间的相互关系由迭代方 法解决。
- 11 -
哈尔滨工业大学工学硕士学位论文
第2章 K-L 展开在随机响应计算中的应用
2.1 引言
近年来， Karhunen-Loéve 展开在随机振动领域有很多应用。本章介绍 这种方法在随机过程中的应用。将随机激励，随机响应利用 K-L 展开进行 离散，特别是将随机响应的数字特征用 K-L 展开表示，可以使计算更加简 便。这种方法的最大优势在于提高了大型有限元体系随机响应协方差计算的 计算效率，使一些不能实现的计算成为可能。
2.2 随机过程及其数字特征
2.2.1 随机过程
如果一个随机变量 X 是一个变量 t 的函数，则称 X (t) 为随机函数。若自 变量 t 为时间坐标，常称 X (t) 为随机过程[34]。
随机过程 X (t) 可看成是一簇函数 x1 (t), x2 (t ), x3 (t ),L xn (t ) 组成的集合，其 中每一个样本函数都是一条确定性的时间历程曲线。
随机过程可以根据本身的取值情况分为连续随机过程和离散随机过程。 如果一个过程 X (t) 对任意的时刻 t i T ， X (t i ) 都是连续型随机变量，则称此 过程为连续型过程;如果 X (t i ) 都是离散随机变量，则称此过程为离散型过 程。
2.2.2 随机过程的数字特征
数字特征可以描述随机过程的概率特征，它们能够集中反映过程的主要 统计特性。在地震工程中，只有少数问题可以从解析的角度给出概率分布函 数解答，大多数实际问题仍然只能给出数字特征解答。
2.2.2.1 随机过程的均值
随机过程的均值(数学期望) , X (t ) 定义为:
， , X (t ) , E[ X (t )] , x , p X ( x, t ) , dx (2-1) ，
式中: p X ( x, t ) 为 X (t) 的一维概率分布密度函数[35]。
- 12 -
第 2 章 K-L 展开在随机响应计算中的应用
2.2.2.2 随机过程的方差
自协方差函数 ,XX (t1 , t 2 ) 定义为:
,XX (t1 , t 2 ) , E[( X (t1 ) , X (t1 )) , ( X (t 2 ) , X (t 2 ))]
， ， (2-2) ( x1 , X (t1 ))( x2 , X (t 2 )) p X ( x1 , t1 ; x2 , t 2 )dx1dx 2 , ， ， 方差函数是随机过程概率分布密度的二阶中心矩。
如果 t1 , t 2 , t ，则:
,XX (t, t ) , E[( X (t ) , X (t )) 2 ] , ,XX (t ) (2-3)
称 ,XX (t ) 为随机过程 X (t) 的方差函数，它表示随机过程的离散范围。
另外，随机过程的标准差定义为:
, X (t ) , ,XX (t ) (2-4)
对于两个随机过程，定义互协方差函数为:
,xy (t1,t2 ) , E[(X (t) , x (t1 ))(Y (t2 ) , y (t2 )]
， ， (2-5)
, ， ， (x , x (t1 ))(y , y (t2 )) pxy (x,t1; y,2t )dxdy
式中 p xy ( x, t1 ; y, t 2 ) 为随机过程 X (t) 和 Y (t) 联合概率分布密度函数。
2.2.2.3 随机过程的相关函数
相关函数定义为:对于随机过程，截取任意两个不同时刻 t1 和 t 2 ，得到 两个随机变量 X (t1 ) 和 X (t 2 ) ，这两个随机变量的相关程度可以用自相关函数
Rx (t1 , t 2 ) 表示，其定义如下:
， ， R X (t 1 , t 2 ) , E[ X ( 1t ) X ( 2t )] , x 2 p X (x 1 , t 1 ; x 2 , t 2 )dx 1 dx 2 (2-6) 1 ， ， x
对于两个不同的随机过程，可定义互相关函数为:
RXY (t1 , t 2 ) , E[ X (t1 ) X (t 2 )]
， ， (2-7) , p XY ( x1 , t1; y2 , t 2 )dx1dy2 ， ， x1 y2
- 13 -
哈尔滨工业大学工学硕士学位论文
2.3 随机过程及数字特征的 K-L 展开表示
2.3.1 随机过程的 K-L 展开表示
K-L 展开是一种类似于级数展开的计算过程，近年来在数值模拟中被广 泛应用，下面以一个连续随机过程 x(t ) 为例，用 K-L 展开可以将这个连续的 随机过程转换成一个离散的随机过程:
N x N x (0) x , , x ， , j ， j , j , x ， , j x ( j ) (2-8)
j ,1 j ,1
式中: , x ——x 的均值向量;
N x ——x 所包含的元素个数;
x ( j ) ——K-L 向量.
x ( j ) , ， j , , j (2-9) 式中: ， j ,, j 分别为协方差矩阵 ,xx 的特征值和特征向量，可以通过求解特征 值问题求解。
,xx , , j , ， j , , j (2-10)
其中特征向量 , j 有如下正交关系:
i , j ,1 (2-11) i , j ;0
方程(2-8)中的正态随机变量 , 定义如下:
1 (2-12) , j , , x T , , j ， j
i , j ,1 并且有 E,, i , , ,j , , (2-13) E{, j } , 0 i , j ;0
N
随机变量 ,, j ,j ,1可以被假设是统计独立的。
这样，就将连续分布的随机过程 x(t ) 用 K-L 展开独立定义了，下面介绍
- 14 -
第 2 章 K-L 展开在随机响应计算中的应用
随机过程数学特征的 K-L 展开表示。
2.3.2 数字特征的 K-L 展开表示
2.3.2.1 均值和协方差的 K-L 展开表示
均值可以用如下的 K-L 向量表示:
,(t) , x (0) (t) (2-14) 协方差矩阵的 K-L 向量表示:
N T (2-15) ,xx (t ) , x (i )(t) , x (i )(t ) i ,1 K-L 向量 x (i ) (t ) 直接与特征向量 , i (t ) 和特征值 ， (t ) 有关:
,xx (t) , , i (t) , , i (t) , ，i (t) (2-16)
因此，协方差矩阵可以用如下和的形式表示: ns
(2-17) ,xx (t ) , ，i (t ) ,,i (t ),iT (t ) i ,1
, i T , , i , 1， n x 表示状态向量 x(t)所含元素 其中特征向量有正交关系
个数。事实上，上式最根本的特点在于，当 ns 较大时，所获得的特征值 ，i 在大小上差距较大，一般来说，最大的特征值要比较小的大很多个数量级， 这时，可以将特别小的省略，取 N << ns 作为求和的数量，使计算得到化 简。
取一个较小的正数 , ，假设是 10 10 , , , 10 5 ，当 ，i , N , , , ，1 时，其中 ，1 , ，max ，可以得到下面的大小排列:
，1 , ，2 , L， N , ,，1 , ， N ，1 (2-18)
N T 用 N 代替 ns 可得到: (2-19) ,xx (t ) , ，i (t ), i (t ), i (t ) i ,1 对于(2-15)中的 K-L 向量为: (2-20)x (i ) (t) , , i (t ) ，i (t )
- 15 -
哈尔滨工业大学工学硕士学位论文
特征向量: (2-21) , i (t) , x (i ) (t ) / ，i (t ) ns 2 特征值: (2-22) j ,1
2.3.2.2 线性关系的特征值的 K-L 展开表示
对于线性关系 z , ci xi ， ci 是确定的参数，则可以通过下式计算 z
i
的均值和方差:
均值: , z , ci xi( 0) (2-23) i N s 方差: ( (2-24) j ,1 i
如果有两个线性关系: (2-25) , ci xi ， z 2 , ck xk z1i k
则对于确定的参数 ci 和 ck ，有下式:
N x ( (2-26) 1 2 j ,1 i k
2.4 随机响应在状态空间下的 K-L 展开表示
T , , 对于状态向量 x(t ) , z T (t ) z& T (t ) v T (t ) q T (t) ，我们得到它的响应均
值和协方差如下:
均值: (2-27) , x (t ) , E{x(t )} , {, zT , ,& zT , , vT , , qT }T
方差: ,xx (t ) , E{~x (t )~x T (t )} (2-28)
其中: x z z& v q (2-29) ~(t ) , x(t ) , x (t ) , {~ T , ~ T , ~ T , ~ T }T
对于状态向量之间的协方差有下面关系:
- 16 -
第 2 章 K-L 展开在随机响应计算中的应用
E{u i (t )} , ui(0) (t ) 位移的均值: (2-30)
速度的均值: (2-31) E{u&i (t)} , u&i(0) (t ) N 位移的协方差: (2-32) E{ui ~(t )u j ~(t)} , u i( k )(t ), u (jk )(t ) k ,1
N 速度的协方差: (2-33) E{ui ~( t&,) u~ j&(t ,)} ui&( k ()t ,) u &(jk )(t ) k ,1
N
位移和速度的协方差: (2-34) E{ui ~(t ), u~& j(t )} , ui( k ()t ), u& (jk ()t ) k ,1
N 速度和位移的协方差: (2-35) E{u~&i (t ), q~ j(t)} , u&i( k, ) q (jk ( )t ) k ,1
2.5 随机激励的 K-L 展开表示
2.5.1 地震动随机模型
2.5.1.1 一维及多维地震动随机模型
(1) 一维地震动就是一个方向的地震动，可以是水平方向，也可以是与 水平垂直的方向。大多数结构的地震反应分析和抗震设计都只考虑一维地震 动的作用[36]。
(2) 在任一固定点的实际地震动通常具有三个分量，其中两个水平分量 分别记为 a x (t ) 和 a y (t ) ，竖向分量记为 a z (t ) 。三个地震动分量组成的三维地 震动向量记为: a g (t ) , ,a x (t ) y ( at ) z ( at, )。三维地震动分量可以分别用一 维的地震动随机模型表示。
2.5.1.2 平稳及非平稳随机模型
(1) 平稳随机模型
大量统计表明，地震动的均值是零，因此地震动的平稳随机模型主要是 相关函数，或谱密度的具体表达形式。比较常见的是白噪声模型和过滤白噪 声模型。
(2) 非平稳随机模型
从一些地震记录上可以看到非平稳随机过程。这种非平稳性分为三个阶 段:开始阶段，地震动速度从小变大;接着是平稳阶段，地震动大体保持其 平均强度不变;然后是衰减阶段，地震动比较缓慢的逐渐减小。因此，研究
- 17 -
哈尔滨工业大学工学硕士学位论文
者们通常采用随时间变化的强度函数和平稳过程的乘积来反映这种非平稳 性，其一般形式为:
Y (t) , f (t) X (t) (2-36) 式中 X (t) 为平稳随机过程， f (t ) 为确定性的强度包络函数。 2.5.1.3 白噪声及过滤白噪声模型
白噪声是地震地面加速度过程的一个近似模型，这种模型极其简单，用 此模型来分析结构的随机地震反应非常方便。白噪声模型表示的地震地面加
速度 a g s (t) 的谱密度为:
(2-37) S a g s (， ) , S 0
它在 ( ,， ) 的频域内是常数，与频率无关，称为无频率特性的平稳地震动 模型。
a g s (, ) , 2 S 0, (, ) Ra g s (t) 的相关函数为: (2-38)
过滤白噪声模型中有一种由日本学者田治见宏和金井清提出的模型能较 好的反映地震动频率特性。这种模型是基于这样的事实提出来的:地震引起 基岩的运动，然后通过地震土层传到结构，先假定地震基岩加速度过程为白 噪声，地表土层相当于一个滤波器，处理为单自由度线性体系，其运动方程 可以表示为:
(2-39) X&& g ， 2, g ， g X& g ， ， g2 X g , Z&&(t )
式中 X&& g (t ) 是地震地面相对于基岩的加速度; , g 和 ， g 是作为单自由度体 系的地表土层阻尼比和固有圆频率; Z&&(t ) 是地震基岩加速度，其均值为零， 谱密度为 S 0 ，于是绝对加速度可以表示成:
a g (t ) , Z&&(t ) ， X&& g (t ) (2-40)
即地震通过基岩将震动传给上层建筑，结构的绝对加速度等于地震加速度和 基岩加速度的和，这就是过滤白噪声模型。
2.5.2 过滤白噪声随机激励的 K-L 展开
认为地震动的随机模型是三维的，地震加速度包括水平方向的 a x (t ) 和 a y (t ) ，以及垂直方向的 a z (t ) ，并假设地震动在 x -， y -， z -三个方向上的
- 18 -
第 2 章 K-L 展开在随机响应计算中的应用
作用是各自独立的，这些加速度都作为统计独立的随机过程，被模拟成过滤 白噪声:
a(t) , ,12g, f 1 ， 2, 1g ,1g, f 2 , 22 g, f 3 2, 2 g , 2 g, f 4 (2-41)
滤波器满足下面的微分方程:
1 0 0 , ,, f 1 , , 0 , ,, f 1 , , 0 , 2 , ,, , , 0 2, 1g ,1g d ,, f 2 , , ,1g 0 , , ,, f 2 ,, ， ,,e(t,,)(，t) (2-42) , , , ,, 0 0 1 0 dt ,, f 3 , , ,, f 3 , , 0 , , , 2 2, 2 g , 2 g , 22 g ,; f 4,,, 2, 2 g , 2 g ,, ,; f 4 ,,, ,; 0 ,, ,, ,1g
式中: ， (t ) ——白噪声;
e(t ) ——强度包络函数。
将(2-41)式简写成: (2-43) a(t ) , Q f, (t )
其中, (t ) 满足下面的微分方程:
(2-44) ,&(t ) , A f, (t ) ， b f (t ) , ， (t ) 式中:, ——状态向量; A f ——系统矩阵;
b f ——滤波器分布矩阵。 ，(t ) 表示白噪声，且有:
E{，(t)， T (t ， , )} , I (t), (, ), , I (t) , 2 S 0 (t) , (2-45)
式中: I (t ) ——白噪声强度，
, (, ) ——Dirac’s 函数，
S 0 (t) , S 0 ——白噪声强度的固定谱密度。
对于每一个时间点 t，滤波器参数的相关性可以用协方差 ,,, (t) 来计算， ,,, (t ) 满足 Lyapunov 微分方程:
,&,, (t) , A f ,,, (t) ， ,,, (t) ATf ， b f (t)Ib Tf (t) (2-46) 对于 ,,, (t1 , t 2 ) 则表示在两个时间间隔 t1 ， t 2 上的滤波器参数的相关特性。它
- 19 -
哈尔滨工业大学工学硕士学位论文
满足如下微分方程:
,,, (t1 , t 2 ) , A f ,,, (t1 , t 2 ) (2-47) 1 t
协方差矩阵 ,aa (t1 , t 2 ) 可通过方程(2-43)表示的关系来求解。 ,aa (t1 , t 2 ) 满足:
,aa (t1, t 2 ) , E{Q f, (t1 ), T (t 2 )QTf }
, Q f E{Q f, (t1 ), T (t 2 )}QTf (2-48)
, Q f ,,, (t1, t 2 )QTf
这样，当 ,aa (t1 , t 2 ) 解出后，继续求解它的特征值及特征向量，代入 (2- 43)过滤白噪声 a(t ) 就可以用确定的 K-L 向量来表示了，就完成了对地震加 速度 a(t)的 K-L 分解。
2.6 K-L 展开在计算中的量纲标准化
方程(2-8)和方程(2-19)中的截断后的数量是 N，它代替了 ns ，就是将对 协方差矩阵影响不大的量去除掉，这样，在应用 K-L 展开时，状态向量元 素 {xk }kN,1 的各量就在计算中起主要作用。为了实现一阶状态空间的降阶，必 须将状态空间中的各量进行统一标准化，这样当进行 K-L 展开时，计算的 结果就会更加准确。因为状态向量里包含位移、速度、加速度、滤波器、响 应、含有阻尼参数的非线性关系等，这些不同的量要用不同的单位来表示， 如秒、米、小时等等。因为这些单位的使用，会使状态向量里的各个量之间 相差很多个数量级，这就在进行 K-L 展开时影响了计算的准确性，这是我 们不希望有的，必须加以改进，办法就是引进一种标准化的方法，使状态向 量在数量级上得到统一。我们用一个参照向量将状态向量 x(t ) 转化成标准化 的向量 y(t) :
x (t ) (2-49) y k (t ) , k 或者 xk (t ) , s k y k (t ) s k
将上式应用到 K-L 展开中，就可以得到下面的标准化的状态向量:
N
(2-50) y(t ) , y (0) (t ) ， , j , y ( j )(t ) j ,1
由上文所提到的原因，其中很重要的一点是，K-L 向量的数量 N 要由 标准化的状态向量 y 来求得，不是由 x 来求得的。这样，对于标准化的状
- 20 -
第 2 章 K-L 展开在随机响应计算中的应用
态向量 y ( j ) (t) 来说，它本身是正交的，而对于原始的状态向量 x ( j ) (t ) 就不再 是相互正交的了。
根据有限的选择合适参数向量 s 的累积经验，我们应用以下的选取方式:
(2-51) si , ci ， , i
这里 ci 表示一个很小的正数，而 , i 是响应 xi 的标准差。 ci 可以是一个 比标准差 , i 小几个数量级的正数，这样就可以避免在方程(2-49)中的除数为 零。同时保证滤波器状态向量 v j 标准化:
s j nv ) s j , (2-52) 4 n x
x 表示状态向量 x 的数量， nv 表示滤波器 v 的状态向量个数。 这里 n
2.7 本章小结
本章引入 K-L 展开的表示方法，将随机过程，随机激励，随机响应都 用 K-L 展开表示，在计算上起到了简化的作用。通过 K-L 展开可以使连续 的随机过程离散化;将大型的协方差矩阵在省略一些对计算影响较小的量 后，使矩阵规模减小，从而使计算简便;对于连续的过滤白噪声激励，通过 K-L 展开使其离散，在进行大型有限元系统随机响应计算时，会使计算得到 简化。本章最后介绍了量纲标准化的方法，在进行 K-L 展开计算过程中， 由于各个物理量单位不同，可能会产生计算误差，将量纲标准化可以使计算 误差减小，这在进行大型结构的计算时，保证结果收敛是非常有用的。
- 21 -
哈尔滨工业大学工学硕士学位论文
第3章 李雅普诺夫微分方程的确定性积分求解
3.1 引言
对于很多工程问题而言，结构随机响应的方差分析十分重要。但是除了 少量简单问题可以直接计算出响应方差以外，大多数情况下的方差求解是很 困难的。在本章中，主要讨论结构随机响应的协方差矩阵所满足的李雅普诺 夫方程的求解过程。这里引入一种基于 K-L 展开的逐步积分方法，来求解 李雅普诺夫方程，即计算随机响应二阶统计矩。下面将计算过程作以详细的 介绍。
3.2 方差矩阵的李雅普诺夫微分方程
这里首先介绍由随机运动微分方程建立李雅普诺夫微分方程。
离散化 n 自由度结构受均匀调制演变随机激励 { f (t)} 时的运动微分方程 可表示为:
(3-1) [M ]{&y&} ， [C ]{y&} ， [K ]{y} , { f (t)} , g (t){x(t )} 其中 {x(t)}为平稳高斯白噪声随机过程向量， g (t) 为调制函数。该方程可表 示为:
y, , 0 , I ,, y, , 0 , & ,x(t ,) (3-2) , , ， , , , , 1 1 , M K M C ,,; y& , ,, g (t ),,I ;&,y&
其中 ,I ,为 n 阶单位矩阵，记:
y 0 I , 0 , ,z, , ,, ,, , ,A, , ,, ,G, , ,, (3-3) 1 1 , , , ; y& , , M K M C , ,g (t ),I
则式(3-2)可以表示为:
(3-4) ,z&, , ,A, ,z,， ,G, ,x,
当式(3-1)的右端为零时，方程(3-4)的解为:
,,(t,0 t ), , ,,(t t0 ), , exp([ A ]( t0t )) (3-5)
而当式(3-1)的右端不为零时，式(3-4)的解为:
t
(3-6) ,z(t,) , ,,(t, 0t ), ,z 0 , ， ，t ,,(t, s, ,) f (s,) ds 0
- 22 -
第 3 章 李雅普诺夫微分方程的确定性积分求解
其中 ,z 0 ,为 ,z(t ),在初始时刻 t 0 的值。通常可认为初始条件与外激励无关，从
而由式(3-6)可得响应的平均向量:
t
(3-7) E[{z(t )}] , [,(t, t0 )]E[{ z0 }] ， ， [,(t, s)]E[{ f (s)}]ds t0 以及相关矩阵:
[Rzz (t1 , t 2 )] , E[{z(t1 )}{z(t 2 )}T ] t t 2 , [,(t1 , t 0 )][Rzz (t 0 , t 0 )][,(t 2 , t 0 )]T ， ， ds1 ， [,(t1 , s1 )] t0 t0 T T (3-8) , E[{ f (s1 )}{ f (s 2 )} ][,(t 2 , s 2 )] ds 2 t1 t 2 , [,(t1 , t 0 )][Rzz (t 0 , t 0 )][,(t 2 , t 0 )]T ， ， ds1 ， [,(t1 , s1 )][G(s1 )] t0 t0
T T T , E[{x(s1 )}{x(s 2 )} ][G(s 2 )] [,(t 2 , s 2 )] ds 2
因 {x(t )] 为平稳白噪声过程，设其强度函数为[D]，则其相关函数为:
[Rxx (s1 (3-9) , s 2 )] , E[{x(s1 )}{x(s 2 )}T ] , [D], (s 2 s1 ) 从而(3-8)为: t1 t 2 [Rzz (t1 , t 2 )] , [,(t, t0 )][Rzz (t 0 , t 0 )][,(t 2 , t 0 )]T ， ， ds1 ， [,(t1 , s1 )] t 0 t0 (3-10) T T , [G(s1 )][D,] (s 2 s1 )[G(s 2 )] [,(t 2 , s 2 )] ds 2 若再令 t1 , t 2 , t ，则式(3-10)成为:
[,] , E[{z(t )}{z(t )}T ] , [,(t, t0 )][Rzz (t0 , t0 )][,(t, t0 )]T
t (3-11) ， ， [,(t, s)][G(s)][D][G(T s[)],(t, s)]T ds t0 对 t 求导得:
[ ,& ] , E [{ z& ( t )}{ z ( t )} T ] ， E [{ z ( t )}{ z& ( t )} T ] , [ ,& ( t , t 0 )][ R zz ( t 0 , t 0 )][ , ( t , t 0 )] T
， [ , ( t , t 0 )][ R zz ( t 0 , t 0 )][ , ( t , t 0 )] T & (3-12) t
， [ ,& ( t , s )][ G ( s )][ D ][ G ( s )] T [ , ( t , s )] T ds ， t 0
t ， [ , ( t , s )][ G ( s )] [ D ][ G ( s )] T [ ,& ( t , s )] T ds ， t 0
， [ , ( t , t )][ G ( t )][ D ][ G ( t )] T [ , ( t , t )] T
- 23 -
哈尔滨工业大学工学硕士学位论文
因为基本解 [,]满足[37]:
[,& (t, t 0 )] , [ A(t )][,(t, t 0 )] ，
[,(t, t)] , [I ] (3-13)
故得方差矩阵 [,] 满足 Lyapunov 微分方程: [,& ] , [ A][,] ， [,][ A]T ， [G(t)][D][G(t)]T (3-14)
进一步等于:
[,& ] , [ A][,]T ， [,][ A]T ， [B] (3-15)
,0 0, 0 0 ,0 0 0 0,, 式中: B(t ) , , (3-16) 0 GI (t )G T ,0 0,
, , 0 0 ,0 0,
E{w(t )wT (t ， , )} , I (t ), (t ) (3-17) I (t) , 2 S 0 (t)
这里 I(t)表示白噪声强度， , (t ) 是 Dirac’s 函数， S 0 (t ) 是白噪声强度的固定
谱密度。
3.3 李雅普诺夫微分方程的求解
在本小节中，介绍基于 K-L 展开的确定性积分方法求解李雅普诺夫方 程，在这个求解方法中，主要用到 K-L 向量和逐步积分法，K-L 向量在上 一章中已经有了介绍，下面先介绍一下逐步积分法的计算过程。
3.3.1 逐步积分法计算过程
这里先介绍方法中涉及的逐步积分法。目前，在结构振动中微分方程的 求 解 需 要 逐 步 积 分 方 法 来 实 现 [38-39] ， 其 中 常 见 的 有 : 线 性 加 速 度 法 ， Wilson , 法，平均加速度法， Newmark , 法，龙格-库特法等。对于多自 由度体系的动力分析，从稳定性与计算精度的角度考虑，多采用平均加速度 法 ， 即 Newmark , 法 中 取 , , 0.25 。 为 了 说 明 问 题 ， 将 本 文 所 使 用 的 Newmark , 法进行详细说明， Newmark , 法是一种将线性加速度方法普 遍化的方法，首先介绍线性加速度法。
- 24 -
第 3 章 李雅普诺夫微分方程的确定性积分求解
3.3.1.1 线性加速度法
这种方法的基本假定是，质点的加速度反应在任意微小时段 t 内的变
化呈线性关系，假设在时段 t , t j ，1 t j 内，结构加速度反应是关于时间 , 的
线性函数，即取:
{U&&} j ，1 {U&&} j { U&&} j , , 常量 (3-18) {U&&&} j , t t
将位移 {U } 按泰勒级数在 t j 附近展开:
{U j & }{U j &&}{U j &&&}, ， , 2 ， {U ( jt ， , )} , {U } j ， (3-19) , 3 ， L 1 ! 2 ! 3 !
{U j &&&}对 , 求导，得: (3-20) {U& ( j ，t , )} , {U j& } ， {U&&} j , ， , 2 ， L 2
当 , , t 时，显然， {U (t j ， , } , {U } j ，1而式(3-19)和式(3-20)分别演化为:
(3-21) { U } j , {U& } j t ， 1 / 2{U&&} j t 2 ， 1 / 6{U&&} j t 2
和: (3-22) { U& } j , {U&&} j t ， 1/ 2{ U&&} j t 2
{ U } j {U& } j 6 (3-23) { U&&} j , 6 3{U&&} j t 2 t
U } j {{ U& } j , 3 j 1 / 2{U&&} j t 3{U& }进而得到: (3-24) t
6 3 [M ] ， [C ] ， [K ] 令: (3-25) [K ] j , 2 t t
(3-26) [K ] j X (t ， t ) , P(t ， t )
t ( t ) 2 ([M ] ， [C ] ， [K ]),U&&,j ，1 2 6 (3-27)
t {U } j ) [K ]({U j } ， t{U } j ， , , P, j [C ]({U j& } ， && & ( t ) 2 {U&& j } ) 2 3
从而通过等效刚度 [K ] 和等效荷载 P(t ， t ) 就可以将 t ， t 时刻的位移，
- 25 -
哈尔滨工业大学工学硕士学位论文
速度，及加速度求出来。
3.3.1.2 Newmark , 法
Newmark , 法是一种将线性加速度方法普遍化的方法，可以认为是概 括了平均加速度和线性加速度算法的广义算法。在 Newmark , 法中，有下 面的式子:
(3-28) {U } j ，1 , {U } j ， {U& } j t ， (1 / 2 , ){U&&} j ( t ) 2 ， , {U&&} j ，1 ( t ) 2
(3-29) {U& } j ，1 , {U& } j ， (1 ， ){U&&} j t ， ，{U&&} j ，1 t
Newmark , 法当取 ， , 1/ 2 和 , , 1 / 4 时(平均加速度法)，为无条件稳定 的，当取 ， , 1 / 2 和 , , 1 / 6 时为有条件稳定的(线性加速度法)，一般来说 , 取值在 1 / 4 ~ 1 / 6 之间比较好。
计算步骤如下:
(1)形成刚度矩阵[K]，质量矩阵[M]和阻尼矩阵[C];
(2)计算初始值 {U 0 },{U& 0 },{U&&0 };
(3)选择时间步长 t ，参数 ， 和 , ，以及计算积分常数:
1 , 1 1 , 1 a0 , ， a1 , , a2 , , a3 , 1, a 4 , 2 2， ， t ， t ， ， t
t , ( 2), a6 , t (1 , ), a7 , , t ; a5 , 2 ,
) ) (4)形成有效刚度矩阵 [K :] [K ] , [K ，] a0 [M ] ， a1[C ] ; ) (5)对 [K ]进行三角分解: [Kˆ ] , [L][D][L]T
(6)计算在时刻 t ， t 有效荷载:
,P ， t ,, ,P ， t ,， [M0 ]( {a tU } ， a 2 { tU } ， a3{U t }) ， [C ](1{ tUa } ， a 4 { tU } ， a5{ tU }) ˆt ˆ & && & &&t
(7)求解在时刻 t ， t 的位移: [L][D][L]T {U t ， t } , {Pˆt ， t }
(8)计算在时刻 t ， t 的速度和加速度:
{U&&t ， t } , a6 ({U t ， t } {U t }) a 2 {U& t } a3{U&&t }
} {U& & && && } , {U t } ， a6 {U t } ， a7 {U t ， t t ， t
- 26 -
第 3 章 李雅普诺夫微分方程的确定性积分求解
3.3.2 基于 K-L 展开的逐步积分法求解李雅普诺夫微分方程
首先将运动方程转化为一阶状态方程形式:
(3-30) z&(t ) , A(t ) z(t) ， h(t )
其中: z(t) , x(t ) , (t ) (3-31)
T T T T (G， (t)) h(t , ) {0 0 (3-32) 0 }
而 ,(t) , E{z(t ) z T (t )} 并假设在 t 时刻已知， A(t ) 因为在计算中不涉及，所以
在这里不作进一步分析，根据上面的状态空间展开，我们可以得到下面的关
于协方差 ,(t) 的李雅普诺夫方程:
,& (t) , A(t),(t) ， ,(t) AT (t) ， B(t ) (3-33)
0 0 ,0 0,
0 0 ,0 0,, 其中: B(t ) , , (3-34) ,0 0 GI (t )G T 0,
, , 0 0 ,0 0,
因为李雅普诺夫微分方程是线性的，其解可以由下面的两个部分和的形式表
示:
(3-35) ,(t k ， , ) , ,1 (t k ， , ) ， ,2 (t k ， , )
初始条件为: (3-36) ,1 (t k ) , ,(t k ); ,2 (t k ) , 0
微分方程可以化为下面两个方程:
,&1 (t k ， , ) , A(t k ， , ),1 (t k ， , ) ， ,1 (t k ， , ) AT (t k ， , ) (3-37)
,& 2 (tk ， , ) , A(kt ， , ),2 (tk ， , ) ， ,2 (tk ， , ) AT (tk ， , ) ， B(tk ， , ) (3-38)
首先考虑方程(3-37)的解，它表示在零荷载时对于一个时间间隔 , 的协
方差矩阵 ,1 (t k ， , ) 的计算。对于每一个时间点 t k 都有下式:
N
(3-39) z(t) , , j x ( j )(t k ) j ,1 很明显，对于零外荷载在时间点 (t k ， , ) 上有下式: N
(3-40)z (t k ， , ) , , j , x ( j )(t k ， , ) j ,1
- 27 -
哈尔滨工业大学工学硕士学位论文
有初始条件: x ( j ) (t k ) , x ( j ) (t k ) ，对于零随机荷载 G， (t) , 0 ，且 1 , j , N ，
每一个固定的向量 x ( j ) (t k ， , ) 都满足线性运动方程:
(3-41) Mu&&( j ) (t ) ， Cu& ( j ) (t ) ， Ku ( j ) (t ) , Qv ( j ) (t ) ， Rq ( j ) (t )
v& ( j ) (t) , Fv ( j ) (t) (3-42)
(3-43) q& ( j ) (t ) , Uu ( j ) (t ) ， Vu ( j ) (t ) ， Wq ( j ) (t )
有: (3-44) x ( j )T , {u ( j )T , u& ( j )T , v ( j )T , q ( j )T } 或者: (3-45) x ( j ) ,T {u ( j ) ,T u& ( j ) ,T u&&( j ) ,T v ( j ) ,T q ( j ) } T
(t k ， , ) 通过确定的逐步积分法来解，在算例中用上小节提到的 K-L 向量 x ( j )T
Newmark 方法来进行计算，这样对于每个时间步长 t 对向量 x ( j ) (t k ， , ) 进行
积分得到 ,1 (t k ， t ) 的解:
N
(3-46) ,1 (t k ， t ) , x ( j )(t k ， t ) x ( j ) T(t k ， t) j ,1
下面，考虑第二部分 ,2 (t k ， t) ，因为 ,2 (t k ) , 0 ，所以对于很小的时间步长
t 和连续白噪声激励有:
(3-47) ,2 (t k ， t) , tB(t k ， t )
对于离散的白噪声，独立的随机激励只作用在离散的时间点 t k , k t 上，那
么，(3-47)就可以就有确切的解
即: (3-48) ,2 (t k ， t) , tB(t k ， t )
考虑更为特殊的情况，在每一个时间点 t k , k t 上协方差矩阵都不连续，
- 28 -
第 3 章 李雅普诺夫微分方程的确定性积分求解
就有: (3-49) ,2 (t k ， t ) , 0 ;
(3-50) ,2 (t k ， t ， ) , tB(t k ， t ) ;
(3-51) ,(t k ， t ， ) , ,1 (t k ， t ) ， tB(t k ， t ) 这里:“-”表示在时间点之前，“+”表示在时间点之后，为了使连续白噪声 的计算更加精确，我们采用脉冲前后的平均值来计算:
,(t k ， t ) , ,1 (t k ， t ) ， 0.5 tB(t k ， t ) (3-52)
3.4 子空间分解
通过逐步积分法对方程(3-40)中的 K-L 向量 x ( j ) (t k ， , ) 进行积分时，用 到的是最原始的状态空间。但是，当要解决较大系统时，要对 K-L 向量进 行进一步更新，必须对状态向量 x(t k ) 作标准化处理。因此，要用到标准化 的状态向量的协方差矩阵。
N1
(3-53) ,(t k ) , tB ( tk ) ， y ( j )(t k ) y ( j ) T(t k ) j ,1
这里用到方程(3-46)和(2-49)的标准化通过 x ( j ) (t k ) 得到向量 y ( j ) (t k ) ，系数 B
定义如下: (3-54) Bij , Bij /(s i s j )
下面主要计算在展开中要用到的 K-L 向量 y ( j ) (t k ) :
N 2
(3-55) ,(t k ) , y ( j )(t k ) y ( j )T (t k ) j ,1
要通过解方程 (2-16)标准特征值问题得到向量 y ( j ) (t k ) ，从方程 (3-53)可以通 过 K-L 向量得到子空间由两部分构成，一部分是由 K-L 向量 y ( j ) (t k ) 组成 的，一部分是由矩阵 B(t k ) 或等效的矩阵 B (t k ) 组成的。方程 (3-34)表明矩阵 B(t k ) 是离散的，矩阵维度 nw 等效于白噪声向量 w(t ) 中的非零元素数量，一 般是 1，或 2，很少能大于 10。这样，子空间的维度是 N1 ， nw ，这样，特征 值问题可以在减小的空间下解决，同时不会引起误差。矩阵 Y 包含所有的向 量 y ( j ) ，矩阵 Yˆ 包含能够代表空间 B 的向量。
- 29 -
哈尔滨工业大学工学硕士学位论文
(3-56) Y , [ y (1) , y ( 2) L y ( N1 ) ]
Yˆ , [ y ( N1 ，1) , y ( N1 ， 2) ,L y ( N1 ， nw ) ] (3-57)
Yˆ , [0T ,0T , G T ,0T ] 由方程(3-42)得 (3-58)
这两部分的向量合起来构成矩阵 Y ，进一步代替协方差矩阵 C
Y , [Y ,Yˆ ] (3-59)
因为 Y 中的向量完全与协方差矩阵 , 相关，原始的 n x 维状态空间的特征值问 题就可以转换成 N1 ， nw 维的特征值问题:
O ,, , , , , (3-60)
D , Y T ,Y , PP T ， tHIH T , , P , Y T Y , (3-61) T , H , Y G , , Gij , Gij / si , 原始状态空间的特征向量就得到如下的转换: , , Y, (3-62)
其中 K-L 向量: (3-63) y ( j ) , , j ， j
方程(3-55)中所需要的向量数 N 2 通过应用(2-18)的截断法则进行调整。如果 取 N 2 , N1 ，n w ，则没有截断误差，但是为了避免 N 成线性增加，最小的特 征值 ， j 和特征向量, j 要被省略.
N 2 ， nw 即: , , , ，1 ， j (3-64)
j , N 2 ，1
这种截断将引入对方差在每一个时间步上的低估，但是误差很小忽略不计。 另外，值得注意的是，非线性系统在等效线性化中引入的误差将会很大，这 在计算中是不可避免的。最后，向量 y ( j ) 在进行下一步计算前转换回原始的 K-L 向量 x ( j ) 。因为标准特征值问题可以在一个减小的空间中进行计算，所 以任何一种方法都能适用。
特征值问题最后通过 Jacobi 转换进行解决，这个过程可以使计算更加有效 率，因为下三角矩阵 D 很小，解可以很快收敛。另外，在迭代中要计算线
- 30 -
第 3 章 李雅普诺夫微分方程的确定性积分求解
性化参数，这些参数与协方差矩阵直接相关，采用前面的解作为初始特征 值，可以使结果很快收敛。
3.5 本章小结
本章介绍了基于 K-L 展开的确定性积分方法来求解随机响应协方差所 满足的李雅普诺夫微分方程。这种方法利用 K-L 展开，将响应统计量用 K- L 向量表示，利用逐步积分法进行积分。在计算中，可以将一些特别小的 K-L 向量省略不计，这样实际计算中涉及的 K-L 向量与原协方差矩阵的维 数相比就很小了，这种求解方法可以应用到大型有限元体系的随机响应计算 中。同时介绍了子空间分解的方法，这都是针对为减小大型有限元体系的随 机响应计算难度提出来的。李雅普诺夫微分方程在实际工程中有很大的应用 空间，新的求解方法可以提高方程的求解效率，有很好的推广价值。
- 31 -
哈尔滨工业大学工学硕士学位论文
第4章 非线性体系随机响应计算的求解方法
4.1 引言
本章对非线性体系随机响应计算进行论述，主要结合等效线性化方法， 将非线性体系线性化，然后利用上一章提出的基于 K-L 展开的确定性积分 方法，计算出结构随机响应的二阶统计矩。
4.2 非线性系统运动方程
对于运动方程: (4-1) Mu&& ， Du& ， Ku , f (t )
非线性项在方程右边的 f (t ) 中考虑，并假设结构的刚度矩阵 K ，阻尼矩阵 D ，质量矩阵 M 是固定的，向量 u 和 f 分别表示系统的广义位移和广义力 向量，以上的线性运动方程与下面的一个线性滤波器结合:
, , F, ， Gw (4-2) &
, 表示滤波器的状态向量， ， 是独立的白噪声激励向量， F 表示滤波 式中
器， G 是与白噪声激励有关的矩阵。为计算方便， F 和 G 都假设是随时间 变化的固定矩阵，所有的非线性因素都在非线性向量函数 g 中给与考虑:
q& , g (u, u&, q) (4-3)
运动方程(4-1)中的 f (t ) 包括外部固定的荷载 f d (t ) ，零均值过滤白噪声 Q, (t )
和非线性存储力 Rq(t ) :
f (t ) , f d (t ) ， Q, (t ) ， Rq(t ) (4-4)
矩阵 Q, R 都假设是固定的，用过滤白噪声模拟有色噪声的随机激励，可以 在实际工程中得以更广泛的应用，非线性系统的响应用一个简单的状态向量 来描述: x , {u T , u& T , v T , q T }T 由于激励是一个随机的过滤白噪声，所以状态 向量 X (t) 也是一个随机的向量。随机响应由可靠度密度函数来计算，但对 于大型有限元系统或多自由度系统这种方法就不适宜了，应采用统计的方法
- 32 -
第 4 章 非线性体系随机响应计算的求解方法
来计算。用随机响应的前二阶矩即:均值、方差、和相关性来近似表示体系 的随机响应特征，下面以协方差矩阵 , 的计算为例来说明。
x(t)的均值: (4-5) , (t ) , E{( x(t)} , {, uT (t ), ,& uT (t ),0, , qT (t )}T
,(t) , E{( x(t ) , (t )) , ( x(t ) , (t ))} x(t)的方差 (4-6) , u (t ) , E{u(t )}, , 其中: (4-7)
E{, (t )} , 0 ,,
将均值代入运动方程中，得到:
(4-8) M,&&u (t ) ， D,& u (t ) ， K, u (t ) , f (t ) ， R, q (t )
其中: ,& q (t ) , g (, u (t ), ,& u (t), , q (t )) ， b(t ) (4-9)
(t ) 是一个很小的量，(4-8),(4-9)这两个方程通过方程(4-1),(4-3),(4-4)的 向量 b
期望来解，其中方程(4-2)定义的是一个零均值的过程，即 E{v(t )} , 0 。因为 方程 (4-9)是非线性的，所以这两个关系要通过迭代来求解。任何一种可进 行迭代求解的方法都适合，如 Newton-Raphson 法等都可以求解。当非线性 不是很强的时候，下面的简化过程可以使解得到收敛，否则需要更复杂的方 法进行计算。
如 Newmark 法， 首先，方程 (4-8)需要一个无条件收敛的逐步积分法 ,
Wilson-Theta 法等，由 t 时刻的解来求出 t ， t 时刻的解:
即: [, u (t ), ,& u (t ), ,&&u (t )] [, u (t ， t ), ,& u (t ， t ), ,&&u (t ， t )]
在迭代计算开始时，假设 q(t ， t) , q(t) 。根据前一时刻的解，非线性
方程 (4-9)可以通过假设 , u (, ) 和 ,& u (, ) 已知来进行积分求解，这里的 , 满足
t , , , t ， t ，即是关于 , q (, ) 的微分方程的积分。非线性方程通常是解耦的
局部非线性方程的集合，这样就不需要在解方程 (4-8)时需要的无条件稳定 的方法，方程 (4-9)可以用任何一种正确的，方便的方法进行积分求解，如 Runge-Kutta 法等。这个积分进行完将得到一个对于 q(t ， t ) 的新的估计， 同时方程 (4-8)右手边的 Rq(t ， t ) 得到新的解。迭代开始于解方程 (4-8)，直
- 33 -
哈尔滨工业大学工学硕士学位论文
到向量 q(t ， t ) 的范数在两个子迭代中完全拟合。方程(4-9)中的向量 b(t)取 决于反应的统计量
(4-10) b(t ) , E{g (u(t ), u&(t ), q(t ))} g (, u (t ), ,& u (t ), , q (t )) 当协方差 ,(t) 确定下来后， b(t)可求。有时假设 b(t , , , t ， t ) , b(t ) 也可以 得到准确的计算，因此，要计算均值 , (t ， t ) 的迭代可以通过确定协方差 ,(t) 的迭代来求解。
到目前为止，对于多自由度非线性体系，或者大型非线性有限元体系， 还没有一种精确的解协方差矩阵的方法，但可以通过等效线性化的方法来计 算非线性多自由度系统的二阶矩，显然，通过这样的计算可以使非线性系统 由统计的等效线性化系统来模拟，同时需要引入参数 U (t) , V (t ) , W (t )
g (u(t ), u&(t ), q(t ) g (, u (t ), ,& u (t ), , q (t ))
(4-11) , b(t ) ， U (t )(u(t ) , u (t )) ， V (t)(u&( t) ,& u (t)) ， W (t )(q(t ) , q (t ))
式(4-11)结合方程(4-1),(4-2)就产生一个线性的运动方程系统，要注意的是， (4-11)两边都取期望就可以得到 (4-10)式。线性化参数的计算也需要进行迭 代，因为这些参数取决于以协方差为代表的响应的统计，这种统计线性化的 方法可以在文献 [40-43]中找到详细的说明。线性系统的协方差我们可以用 第 3 章提出的方法进行计算。
结合均值向量，运动方程可以有下面的引申:
M (u&&(t) ,&&u (t)) ， D(u&(t) ,& u (t)) ， K (u(t) , u (t)) , Qv(t) ， R(q(t) , q (t)) (4-12)
(4-13) v&(t) , Fv(t) ， Gw(t)
q& (t ) ,& q (t ) , U (u( t ), u (t )) ， V (u&( t , )& u (t )) ， W (q( t ) , q (t )) (4-14)
线性化的系统协方差矩阵的具体计算，已在第 3 章中介绍。
4.3 等效线性化方法在非线性系统中的应用
工程结构的非线性性质主要分为两类:几何非线性和材料非线性 [44]。 几何非线性主要发生于几何大变形问题。材料非线性主要体现在阻尼项和刚 度项中。当结构的非线性项表现在阻尼项，刚度项或两者的耦合项中时，分
- 34 -
第 4 章 非线性体系随机响应计算的求解方法
别称为阻尼非线性，刚度非线性，或滞变非线性 [45-48]。从非线性随机振动 工程分析的角度来看，如果非线性程度不是太强，不考虑分叉，跳跃等本质 非线性现象，则等效线性化已取得了较大的进展，是一种相对有效，且简便 可行的方法，是工程实践中最具应用潜质的预测非线性系统随机响应的近似 解法。其基本思想是用一个有精确解的线性系统代替给定的非线性系统，并 使两方程之差在某种意义上最小。
若 n 自由度非线性系统方程可表示为:
(4-15) [M ]{&y&} ， [C ]{y&} ， [K ]{y} ， {Q} , {F (t)}
其中 {Q} , {q1 , q 2 ,L, q n }T 为非线性项向量， qi 是 { y }{ y& }及 { &y& }的函数; ,F (t, )为零均值平稳矢量随机过程。设式 (4-15)的平稳响应存在，用以下的
等效线性系统近似替代(4-15)
(4-16) ([M ] ， [M e ]){Y&&} ， ([C ] ， [Ce ]){Y&} ， ([K ] ， [K e ]){Y } , {F (t)}
],[Ce ],[K e ]，分别是 n 维等效质量，阻尼及刚度矩阵。 式中 [M e
为使上式成为(4-15)在最小二乘意义上的最优等效线性系统，须令两方程之
差: {e} , {Q} [M e ]{Y&&} [Ce ]{Y&} [K e ]{Y } (4-17)
的均方值 E[{e}T {e}] 为最小，其必要条件是:
E[{eT {}e}] , E[{e}T {e}]
[M e ]ij [Ce ]ij
, E[{e}T {e}] (4-18) [K e ]ij
, 0(i, j , 1,2,L, n)
[K e ]ij , E[ qi / Y j ]
求解上式可得 : (4-19) [Ce ]ij , E[ qi / Y&j ]
[M e ]ij , E[qi / Y& ] & j
由于等效参数表达式 (4-19)中含有等效线性系统 (4-16)的响应方差矩阵的元 素，因此用等效线性化法求解需要迭代进行。
- 35 -
哈尔滨工业大学工学硕士学位论文
4.4 滞迟系统在平稳随机激励下的随机响应分析
许多工程结构在严重荷载的作用下(如强地震引起的地面运动)，由于
材料的屈服，内摩擦或者构件结合面间的摩擦等原因，往往会呈现出滞迟特
性 [47]，滞迟系统的恢复力是一个多值函数 [48]，不仅取决于系统的瞬时位
移，而且取决于系统响应的历史。滞迟效应主要表现为刚度的减小，与能量
耗散能力的增加。
光滑滞迟系统的等效线性化分析过程如下:
滞迟系统的恢复力可表示为:
(4-20) q( y, y&, t ) , g ( y, y& ) ， cz(t )
式中 g 为非滞迟分量，它是瞬时位移与速度的函数; c , z(t ) 为滞迟分量，满
足如下非线性微分方程
1 n 1 n (4-21) z , &[ Ay ,& (, y& z z , y& z )] ,
其中参数 , 与 , 控制滞迟回线的形状; A,, 与 , 支配滞迟的强度与刚度; n
决定从弹性区到塑性区过渡的光滑性，适当选取这些参数，可给出具有各种
能量耗散能力的滞迟系统。
, , , , , 1时，近于理想弹塑性系统的单自由度滞迟系统的运 当上式中的 n
动微分方程可写为:
, &y& ， 2,， 0 y& ， ，， 0 y ， (1 ， )， 02 z , x(t ) / m , (4-22) z& , Ay& , y& z ,y& z ,
这里为屈服后与屈服前的刚度比，并设为零均值平稳或非平稳随机过程。上
式的等效线性系统为:
, &y& ，,， 2 0 y& ， ，， 02 y ， (1 ， )，02 z , x(t ) / m (4-23) , z& ， C1 y& ， C2 z , 0 ,
等效参数为:
, 2 , E[ y&z] ,, C1 , ， ,, z , A (4-24) ,,
- 36 -
第 4 章 非线性体系随机响应计算的求解方法
2 , E[ y&z], (4-25) C 2 , ,,, y& ， , , , z ,,
T 引入状态向量 q , ,y, y&, ,z ，式(4-23)可以重写为: d q , Gq ， f (4-26) dt
, 0 0 , 1 2 其中: 2,， 0 (4-27) (1 ， )， ,,0 G , ,, ，， 02
, 0 C1 C 2 , , , T f , ,0 x(t) / ,m 0(4-28)
设 q 的协方差矩阵为 , ,即:
yz, , yy yy&
(4-29) y&y& ,(t) , E ,, y&y y&z,,
zy& ,, zy zz ,,
若白噪声强度为 I (t) , 2 S (t ) ，则 , 满足如下的微分方程:
(4-30) ,& ， G, ， ,G T , B 其中: ,0 0 0,
(4-31) B , ,,0 I (t ) 0,,
,0 0 0, , ,
接下来的求解过程就等同于线性的求解过程，但因为方程(4-23)的参数 C1, 和
C 2 取决于协方差矩阵，所以要通过迭代求解。具体计算过程及计算结果在
算例中给出，这里就不加以叙述。
4.5 本章小结
本章主要介绍了求解非线性系统随机响应的一种方法，即结合等效线性
- 37 -
哈尔滨工业大学工学硕士学位论文
化和第三章提出的基于 K-L 展开的确定性积分方法进行计算。首先利用等
效线性化的方法将非线性系统线性化，然后将线性化的系统用基于 K-L 展
开的确定性积分方法进行计算。对于非线性系统，线性化参数直接与协方差
矩阵有关，要用迭代的方法进行求解，这样，非线性系统的随机响应就得到
了很好的计算。
- 38 -
第 5 章 算例分析
第5章 算例分析
5.1 引言
以线性单自由度体系和非线性多自由度体系的随机振动为例，用前面提 出的方法利用 Matlab 软件 [49-50]对其随机响应进行计算，并将所求得的单自 由度体系在白噪声激励下的响应与用解析解求得的响应进行对比。将非线性 多自由度体系的响应与相应的线性体系的响应进行对比，给出计算框图，及 计算结果的图表。
5.2 单自由度体系的随机响应分析
如图所示的一理想化受水平地震作用的单层结构，体系的抗侧移刚度由 无重钢架提供，体系的质量 M 集中作用在横梁中部。假定杆件不发生轴向 变形，对于静力分析问题来说，该结构具有三个静力自由度:一个水平线位 移，两个节点转角;而对于动力分析问题来说，在水平地面运动作用下，一 个独立的水平位移即可以确定质点在水平惯性力的作用下的位移。因此，它 只有一个水平侧移动力自由度。
u
a(t )
图 5-1 单自由度体系的计算模型
运动方程为: (5-1) Mu&& ， Cu& ， Ku , f (t )
各参数假定为: M =1kg, C =0.6 kg , s 1 , K =100 kg , s 2
- 39 -
哈尔滨工业大学工学硕士学位论文
外荷载假定为: f (t ) , Ma(t )
加速度为: a(t ) , w(t )U (t )
U (t) ——单位阶跃函数;
w(t ) ——连续的白噪声激励，其强度为 I 0 , 2 S 0
S 0 ——固定的谱密度，
且有: (5-2) E{w(t )w(t ， , )} , I 0, (, )
其中 , (,) 为 Dirac’ 函数
基于 K-L 展开的确定性积分求解方法如下:
将运动方程转化成状态空间方程:
(5-3) x& , As x ， bs f (t )
,u, , 0 , 其中， x , , , bs (t ) , , , 1;u&, ;M ,
, 0 1 , As , , 1 1 M C ,, , M K
得到关于系统响应方差的李雅普诺夫微分方程:
T T (5-4) ,& xx (t ) , As ,xx (t) ， ,xx (t) As ， bs (t)I 0bs (t )
通过第三章介绍的求解李雅普诺夫微分方程的方法，可以求得上述方程的数
值解，其计算过程流程图由图 5-6 表示出来。
为了对计算结果进行分析，下面利用由杜哈美积分求得的该体系的位移
和速度的方差的解析解进行计算 。
该方程的解析解可以用如下的方程表示:
2 ， ,,, I 0 ,, 2,，t ,， 2 2 (5-5) ,1 e , 2 (1 , cos 2， d t) ， u ，,， d d ,,, 4,， 3 ,;
， ,,, I 0, 2 ,， (5-6) , , u2& (t ) , 1 e 2,，t , 2 (1 , 2 cos 2， d t ) ， d ,， d ,,, 4,， ,;
其中: ， d , ， 1 , 2 ; , 为体系的阻尼比
- 40 -
第 5 章 算例分析
得到计算结果如下
0.030
0.020 0.025
0.020 0.015
0.015 解析解 解析解 0.010 数值解 数值解
0.010 位移标准差(m) 0.005 0.005
0.000 0.000 0 5 10 15 20 0 5 10 15 20
时间(s) 时间(s)
图 5-2 位移标准差( , , 0.03) 图 5-3 位移标准差( , , 0.05)
表 5-1 位移标准差 (, , 0.03)
(m)
t=0s t=5s t=10s t=15s t=20s 位移标准差(m)
解析解 0 0.024982 0.025722 0.025767 0.025769
积分解 0 0.025180 0.025789 0.025818 0.02582
表 5-2 位移标准差 (, , 0.05)
(m)
t=0s t=5s t=10s t=15s t=20s
解析解 0 0.019919 0.02 0.02 0.02
积分解 0 0.01986 0.01996 0.019961 0.019961
- 41 -
速度标准差(m/s)
哈尔滨工业大学工学硕士学位论文
0.30 解析解 0.20 数值解 0.25 解析解 0.15 数值解 0.20 0.15 0.10
0.10 速度标准差(m/s) 0.05
0.05
0.00 0.00 0 5 10 15 20 0 5 10 15 20
时间(s) 时间(s)
图 5-4 速度标准差( , , 0.03 ) 图 5-5 速度标准差( , , 0.05 )
-3 速度标准差( , , 0.03 ) 表 5
(m/s)
0 5s 10s 15s 20s 时 间
解析解 0 0.25158 0.25787 0.25818 0.25820
积分解 0 0.24927 0.25650 0.25691 0.25694
表 5-4 速度标准差( , , 0.05 )
(m/s)
0s 5s 10s 15s 20s 时 间
解析解 0 0.19998 0.20087 0.20087 0.20088
积分解 0 0.199250. 0.2 0.2 0.2
- 42 -
第 5 章 算例分析
, (t)
求 , (t)的特征值 V 特征向量 D
0 , ,V1 , D11 D12 ,
, 0 V2 ,, ,D21 D22 ,, , ,
u(1) (t) , V1 D11 u& (1) (t) , V1 D21
u(2) (t) , V2 D21 u& (2) (t) , V2 D22
( j ) (t ) ， ku ( j ) (t ) , 0 代入方程 mu&&( j ) (t ) ， cu&
用 Newmark 解 u ( j ) (t ， t ), u& ( j ) (t ， t )
z ( j )T (t ， t) , [u ( j ) (t ， t ) u ( j ) (t ， t )]
2
,1 (t ， t ) , z ( j )(t ， t ) z ( j ) T(t ， t) j ,1
,(t ， t ) , ,1 (t ， t) ， 0.5 tB(t ， t )
图 5-6 线性单自由度体系随机响应计算框图
- 43 -
哈尔滨工业大学工学硕士学位论文
5.3 多自由度体系的随机响应分析
如图 5-7 所示, 一个 5 个自由度的滞迟系统,
受地面过滤白噪声激励, 系统各参数取无量纲
形式，具体取值如下:
m1 , m2 , m3 , m4 , m5 , 1.0
c1 , c2 , c3 , c4 , c5 , 0.6
k1 , k 2 , k3 , k 4 , k5 , 200
, i , 0.2 , i , , i , 0.5
Ai , 1.0
i , 1,2,L,5
, g , 0.5
应用基于 K-L 展开的确定性积分方法,结合等
效线性化的方法，解这个多自由度非线性体
系，求解过程及计算结果分析如下:
图 5-7 多自由度体系计算模型
1 系统的运动方程为:
[M ]{Y&&} ， [C ]{Y&} ， [K ]{Y } , [M ]{E}&x&g ， R{Z}
2,通过等效线性化化成线性运动方程组
[M ]{Y&&} ， [C ]{Y&} ， [K e ]{Y } ， [K h ]{Z} , [M ]{E}&x&g
{Z&} ， [P]{Y&} ， [H ]{Z} , 0
式中 [K e ]——弹性刚度阵，
[K h ] ——滞迟刚度阵，
{E}——惯性力指示向量
质量矩阵:
- 44 -
第 5 章 算例分析
0 0 0 0 , ,m1
, 0 00 m2 , 0 ,,
0 0 ] , , 0 [M m3 , 0
, , 0 0 m4 , 0 0 ,
0 0 0 m5 ,, 0 ,,
阻尼矩阵:
0 0 0 , ,c1 ， c2 c2
, c 0 c3 c2 ， c3 2 0 ,, ,
[c] , , 0 0 , 4 c c3 c3 ， c4 , , 0 0 c4 ， c4 c5 c5 , ,
0 , 0 0 c5 ,, c5 ,
弹性刚度阵;
0 0 , ,，1k1 ，，2k2 0 ，2k2 , 0 ，2k2 ，，3k3 ，3k3 0 ,, ,
[Ke ] , , ，3k3 ，，4k4 ，4k4 0 ,
, , Sym ,
，5k5 ,, ,,
滞迟刚度阵:
,(1 ，1)k1 ， (1 ，2 )k2 (1 ，2 )k2 0 0 0 , , , 0 (1 ，2 )k2 ， (1 ，3 )k3 O 0 , ,
0 , [K h ] , , O O , ,
O (1 ，5 ))k5 , ,
, Sym (1 ，5 )k5 ,, ,
T 参数矩阵: {E} , {1 1 1 1 1}
T {Y } , ,y1 y5 , y 2 y3 y 4 T ,Z, , ,z1 z 5 , z 2 z3 z 4
- 45 -
哈尔滨工业大学工学硕士学位论文
0 , 0 , , p1 0 0 0 ,h1 , , 0 0 , p p 2 h2 2 0 ,, , , , 0 , p3 p3 0 , [P] , , 0 [H ] , , h3 , , , , 0 p 4 p 4 h4 , 0 0 , , ,
0 0 p5 ,, 0 p5 ,, ,, 0 h5 ,,
地面加速度由统计独立的随机过程来描述，定义成过滤白噪声:
&x&g (t ) , ,12g v f 1 ， 2, 1g ,1g v f 2 , 22 g v f 3 2, 2 g , 2 g v f 4
滤波器满足下面的微分方程:
1 0 0 , ,v ,v f 1 , , 0 f 1 , , 0 , , , , 2 , 0 , ,v , , 2, 0 1g ,1g d ,v f 2 , , ,1g , , ,, f 2 ,, ， ,,e(t ),,w(t ) , , , 0 1 0 dt ,v f 3 , , 0 , ,v f 3 , , 0 , , ,1g , 22 g 2, 1g,;v f 4 ,, ,,, ,12g 2, 2 g, 2 g,, ,;v f 4 ,, ,; 0 ,,
其中 : ,1g , 15.6rad / s, , 1g , 0.6, , 2 g , 1.0rad / s, , 2 g , 0.995
I , 0.08m 2 / s 2
地震动在本质上是非平稳的随机激励，这种非平稳性质可以通过一个包
络函数体现出来，为计算简单，我们取下面的包络函数:
e at e bt e(t ) , max(e at e bt )
取 a , 0.2 ， b , 0.4 其包络图像见图 5-8
1 .0
0 .8
0 .6 e(t) 0 .4
0 .2
0 .0 0 5 10 15
时 间 ( s)
图 5-8 白噪声激励的包络函数 e(t)
- 46 -
第 5 章 算例分析
,(t) , 0.5 tB
求 , (t)的 特 征 值 [V]特
征向量[D]
z(t) ,[u(t),u&(t),q(t),v(t)]T , V1 D1L Vj Dj ,L Vn Dn
假设 q (j) (t ， t) , q ( j ) (t )代入
m u&&( j ) (t ) ， cu& ( j ) (t ) ， k e u ( j ) (t ) , k h q ( j ) (t ) ， Qv ( j ) (t ) j) j) v& ( j ) (t ) , Fv ( j ) (t ) 求出 u(( t ， t)， u&(( t ， t)
z ( j )* (Tt ， t) , {u ( j )* (t ， t),u& ( j )(t ， t), v ( j )* (t ， t), q ( j )(t ， t)} N T ,1 (t ， t) , z ( j )(t ， t) , z ( j )(t ， t) j,1
,* (t ， t) , ,1* (t ， t) ， 0.5 tB
2 * E[ yi zi ] &[, i ， ,i, zi ] Ai pi (t ， t) , , ui &
2 E[iu zi ] &* (, i, i ， ,i ) hi (t ， t) , , z
q& ， pu& ， hq , 0 求出 q (t ， t)
* ( j ) * 若 q ( ， t t) q(t ， t ) , , 则 z ， (t t ) , z ( j )(t ， t ) n T ,1 (t ， t ) , z ( j ) (t ， t ) z ( j ) (t ， t ) j ,1
,(t ， t ) , ,1 (t ， t ) ， 0.5 tB
图 5-9 非线性多自由度体系计算框图
- 47 -
哈尔滨工业大学工学硕士学位论文
计算流程图如图 5-9 所示
为了分析非线性系统在非平稳随机激励作用下的随机响应的特点，取结构
体系的第一层，第三层，第五层的位移标准差与线性体系随机响应进行对
比，对比结果见图 5-10
各层线性 一层非线性 0.03 三层非线性 五层非线性
0.02
0.01 位移标准差(m)
0.00 0 5 10 15
时间 (s)
图 5-10 非线性体系位移标准差与线性体系位移标准差
表 5-5 非线性体系位移标准差与线性体系位移标准差
(m)
时间 一层 三层 五层
线性 3s 0.0066 0.011 0.0145
非线性 0.00528 0.00935 0.01305 线性 0.011 0.0184 0.0243 6s
非线性 0.0088 0.01564 0.02187
线性 0.012 0.0201 0.0265 9s
非线性 0.0096 0.01708 0.02385 线性 0.0111 0.0185 0.0244 12s 非线性 0.00888 0.01572 0.02196 线性 0.0093 0.0155 0.0204 15s
非线性 0.00744 0.01317 0.01836
由图 5-10 各层绝对位移的标准差可以得出:基于 K-L 展开的确定性积分方
- 48 -
第 5 章 算例分析
法结合等效线性化方法可以求解多自由度非线性系统。由于阻尼增大，迟滞 系统的最大响应时间比相应的线性系统的最大响应时间要早，在滞迟系统 中，最大绝对位移标准差大约出现在 t , 6.5s 时，而线性系统的最大位移大 约出现在 t , 7.8s 时。
5.4 本章小结
通过单自由度体系的算例，通过基于 K-L 展开的确定性积分方法计算 单自由度线性体系的位移和速度标准差，并将其计算结果与解析解进行对 比。可以看到两种方法的计算结果十分吻合，即这种新的求解随机响应的方 法在求解线性系统时非常适用。通过第二个多自由度非线性体系的算例，可 以看出，新方法可以用来求解非线性系统在非平稳随机激励下的随机响应。 通过非线性系统的计算结果与相应线性系统的计算结果的对比，可以看出非 线性系统的响应总体趋势与线性系统一致，由于非线性存储力对阻尼的影 响，使最大响应出现时间提前，最大响应降低。其原因是非线性存储力增加 了阻尼影响，使能量分散，因此降低了响应。这样，可以得出:基于 K-L 展开的确定性积分方法可以应用于线性和非线性体系中，可以计算平稳和非 平稳随机激励，这种方法在实际工程中有很好的应用价值。
- 49 -
哈尔滨工业大学工学硕士学位论文
结 论
本文的主要工作是将一种基于 K-L 展开的确定性积分方法应用到结构 随机响应计算当中，得出以下结论:
(1)通过 K-L 展开，可以将连续的随机过程用其自身的 K-L 向量表示， 即可以将连续的随机过程离散化;可以用 K-L 向量将地震动荷载白噪声或 过滤白噪声进行离散;结构随机响应的一阶、二阶统计矩也可以通过 K-L 向量表示。在应用 K-L 展开时，只采用对计算结果有影响的向量，对于数 量级较小的忽略不计，这样进行实际计算的 K-L 向量在数量上得到一次简 化，对于大型有限元系统，该简化可以节省很多计算机时和存储量。
(2)随机响应的协方差矩阵可以通过它的 K-L 向量表示，而 K-L 向量又 取决于协方差矩阵本身的特征值，特征向量，这样引出一种巧妙的求解协方 差矩阵所满足的李雅普诺夫微分方程的新方法。这种方法的基本原理是将协 方差矩阵用其 K-L 向量表示，将每一个 K-L 向量带入运动微分方程中，结 合逐步积分方法求出下一时间步长上的 K-L 向量，即求出下一时间步上的 协方差矩阵。这种方法的优势在于它适于求解自由度数较大的多自由度体系 和有限元体系。求解过程比较简便，可以提高计算效率。
(3)对于非平稳的随机激励，通过子空间分解在状态方程中考虑，计算 过程与平稳的随机激励几乎相同，使非平稳随机激励下结构响应计算得到简 化，尤其在进行大型有限元体系的非稳态随机响应计算时会很好地体现出计 算优势。
(4)非线性体系的计算，要先通过等效线性化将非线性系统线性化，然 后对于线性化的体系进行基于 K-L 展开的逐步积分方法计算。其中线性化 参数与协方差矩阵相关，需要迭代求解。从算例可以看出，因为滞迟系统使 阻尼增加，能量分散，减小了响应，滞迟系统的随机响应比相应的线性系统 的随机响应要小。
(5)基于 K-L 展开的确定性积分方法在计算随机响应时的限制条件很 少，可以采用很多常用的积分方法进行计算。此方法可以求解线性系统和非 线性系统的随机响应，可以求解平稳和非平稳随机激励下的结构响应，适于 应用在实际工程中。
- 50 -
参考文献:
参考文献:
陈英俊, 甘幼深, 于希哲. 结构随机振动. 人民交通出版社, 1993: 349 1
2 黄玉盈. 结构振动分析基础. 华中工学院出版社，1988: 56
3 胡少伟 , 苗同臣 . 结构振动理论及应用 . 中国建筑工业出版社 , 2005:
165~170 方同. 工程随机振动. 国防出版社, 1995: 26 4
5 R.W.克 拉 夫 , J.彭 津 . 结 构 动 力 学 . 王 光 远 等 译 . 科 学 出 版 社 , 1983:
405~406 朱位秋. 非线性随机动力学与控制—Hamilton 理论体系框架. 科学出版 6 社, 2003: 186
L. D. Lutes, S. Sarkani. Stochastic Analysis of Structural and Mechanical 7 Vibrations. Prentice-Hall. 1997: 4
H. J. Pradwarter, G. I. Schueller. Stochastic Structural Dynamics—A Primer 8 for Practical Applications Uncertainty Modeling in Vibration. Control and
Fuzzy Analysis of Structural Systems. World Scientific. 1997, vol.10: 128 C. A. Schenk, H. J. Pradlwarter, G. I. Schueller. On the Dynamic Stochastic 9 Response of FE-Models. Probabilistic Engineering Mechanics. 2004: 161
李贵青. 抗震结构计算理论和方法. 地震出版社, 1985: 89~140 10
11 庄表中 , 陈乃立 , 高瞻 . 非线性结构振动理论及应用 . 浙江大学出版社 ,
1986: 15
欧进萍, 王光远. 结构随机振动. 高等教育出版社, 1998: 180~199 12
13 徐日旭等编著. 结构振动与振动测试. 东北林业大学出版社, 1993: 89
14 李国豪. 工程结构抗震动力学. 上海科学技术出版社, 1980: 11~45
15 帕兹 M. 结构动力学. 高等教育出版社, 1993: 258~260
16 Ray. W. Clough, Joseph. Penzien. Dynamics of Structures. Computers&
Structures, 1995: 2~5
H. J. Pradwarter, G. I. Schueller, C. A. Schenk. A Computational Procedure to 17
Estimate the Stochastic Dynamic Response of Large Non-linear FE-Models. Comput Methods Appl. Mech. Engrg, 192(2003):777-801 C. Prope, H. J. Pradlwarter, G. I. Schueller. Equivalent Linearization and 18
Monte Carlo Simulation in Stochastic Dynamics. Probabilistic Engineering
- 51 -
哈尔滨工业大学工学硕士学位论文
Mechanics, 2003: 1-15
19 H. J. Pradlwarter. Deterministic Integration Algorithms for Stochastic
Response Computations of FE-systems. Computers and Structures. 80 2002, 80: 1489-1502 Ray .W .Clough, Joseph. Penzien. Dynamics. of. Structures. Computers & 20 Structures, 1995: 65~69 Y. Ohtori, B. F. Spencer. Semi-Implicit Integration Algorithm for Solution of 21 Non-linear Stochastic Vibration Problems. In 8th ASCE Specialty Conference on Probabilistic Mechanics and Structural Reliability. 2000, Vol. CD-ROM,
PMC: 323 H. J. Pradlwarter, G. I. Schueller, U. Dorka. Reliability of MDOF-Systems 22 ith Hysteretic Devices. Eng Struct. 1998, 20(8): 685-91 w G. I. Scheller, H. J. Pradlwarter. On the Stochastic Response of Non-linear 23 FE-Models. Arch Appl Mech. 1999, 69: 765-84
M. Vasta, G. I. Schueller. Phase Space Reduction in Stochastic Dynamics. J 24 Eng Mech. 2000, 126(6): 626
H. H. Emam, H. P. Pradwarter, G. I. Schueller. A Computational Procedure 25
for the Implementation of Equivalent Linearization in Finite Element
Analysis. Earthq Engng Struct Dyn. 2000, 29: 1-7
C. A. Schenk, H. J. Pradlwarter, G. I. Schueller. Nonstationary Response of 26
Large Non-Linear FE Systems Under Stochastic Loading. 2003: submitted
A. W. Smyth, S. F. Masri. Nonstationary Response of Nonlinear Systems 27 Using Equivalent Linearization with a Compact Analytical Form of the
Excitation Process. Probab Engng Mech. 2002, 17: 97-108 H. H. Emam, H. P. Pradlwarter, G. I. Schueller. A Computational Procedure 28
for the Implementation of Equivalent Linearization in Finite Element
Analysis. J .Earthquake Eng. Struct. Dyn. 2000, 29: 1-17 Y. Ohtori, B. F. Spencer. Semi-Implicit Integration Algorithm for Solution of 29
Nonlinear Stochastic Vibration Problems. 8th ASCE Specialty Conference on Probobilistic Mechanics and Structural Reliability. 2000，volume PMC: 223 L. Zhang, J. W. Zu, Z. Zheng. The stochastic Newmark Algorithm for 30
Random Analysis of Multi-Degree-of-Freedom Nonlinear System. Comput
Struct. 1999, 70: 557
- 52 -
参考文献:
31 To CWS. A Stochastic Version of the Newmark Family of Algorithm for
Discretized Dynamics Systems. Comput Struct 1992, 44(3): 667 32 吴勃英主编. 数值分析原理. 科学出版社, 2003: 73~75
33 H. J. Pradlwarter, G. I. Schueller. Stochastic Response Evaluation of Large
FE-Models with Nonlinear Hysteretic Elements by Complex Mode Synthesis.
Spanos PD, editor 3rd International Conference on Computational Stochastic
Mechasnics, Santorini, Greeese. The Netherlands, A. A. Balkema Publ, 1993:
551
34 H. J. Pradlwarter, Li W-L. On the Computation of the Stochastic Response of
Highly Nonlinear Large MDOF-Systems. Probabilist Eng Mech. 1991, 6(2):
109
35 To CWS. Recursive Expression for Random Response of Nonlinear Systems.
Comput Struct. 1988, 29(3): 451
36 李宏南, 李忠献等. 结构振动与控制. 中国建筑工业出版社, 2005: 149
37 徐仲, 张凯院, 陆全, 冷国伟. 矩阵论简明教程. 科学出版社，2005: 25
38 王忠, 王雅琳, 王芳. 非平稳激励下结构动力响应的状态方程精细积分法.
计算力学学报. 2002 年 11 月，Vol.19，No.4
39 薛笑荣 . 非线性动力学方程的高精度计算研究 . 西北工业大学硕士毕业
论文. 2001: 2
40 R. Ghanem, P. D. Spanos. Stochastic Finite Elements: a Spectral Approach.
Heidelberg, Springer, 1991: 26~30
41 M. A. Crisfield. Non-linear Finite Element Analysis of Solids and Structures.
Essentials New York; John Wiley and Sons, 1991: 15
42 K. L. Bathe. Finite Element Procedures. Prentice-Hall; 1999: 3 43 J. B. Roberts, P. D. Spanos. Random Vibration and Statistical Linearization.
John Wiley and Sons, 1990: 4~5
44 殷学纲, 陈淮, 蹇开林编著. 结构振动分析的子结构方法. 中国铁道出版
社，1991: 69~72
45 林家浩, 张亚辉. 随机振动的虚拟激励法. 科学出版社, 2004: 366~385
46 张春宇 . 大型结构局部非线性随机地震响应分析 . 大连理工大学硕士毕
业论文. 2002: 16~19
47 K. Subbaraj, M. A. Dokainish. A Survey of Direct Time-Integration Methods
in Computational Structural Dynamics-?Implicit Method. Comput Struct
- 53 -
哈尔滨工业大学工学硕士学位论文
1989, 32(6): 1387
48 B. Miao. Direct Integration Variance Prediction of Random Response of
Nonlinear Systems. Comput Struct. 1993, 46(6): 979 49 薛定宇 , 陈阳泉 . 高等应用数学问题的 Matlab 求解 . 清华大学出版社 ,
2004: 198
50 王沫然. Matlab 与科学计算. 第二版. 电子工业出版社, 2005: 19~185
- 54 -
致谢
致谢
本文的工作是在导师吴知丰教授的悉心指导下顺利完成的，这里首先要 感谢导师吴知丰教授。两年的研究生学习期间，我很庆幸自己能得到吴老师 的耐心指导，使自己在学习方法和思维方式上都得到新的启迪，思路和视野 得到开阔。吴老师丰厚渊博的知识，敏捷开阔的思维，锐意开拓的创新精 神，让学生特别倾佩，老师踏实求是的治学态度和宽厚待人的人格魅力，更 让学生受益匪浅，是学生终生学习的榜样。在此向导师致以最衷心的感谢。
同时感谢同教研室的师兄、师妹、师弟的支持和帮助。
感谢每一位在课题进展过程中给予我帮助和鼓励的同学和朋友，特别是 孙亚民，王健，何长江等同学的鼎力相助，使很多理论和程序上的问题都得 到了顺利的解决。
感谢学校的培养，感谢每一位给予我指导和帮助的老师。
感谢默默支持我的父母和家人。
岳 娟
二〇〇六年六月二十日
- 55 -
浙公网安备 33021202002483