（固体力学专业论文）一种时域有限元法求解动力学和传热学问题.pdf

大连理工大学硕士学位论文 萄矍 工程中存在大量与时间有关的问题例如热传导或者弹性动力学这类的瞬态场和动态 问题。时域问题的研究，具有重要的理论探讨价值和实际工程意义。 本文给出一种在运用适当的空间离散技术后建立一般的数值正演模型的解决方法， 在时域上利用加权余量原理，给出基于时域的递推格式，推导出一种时域有限元算法。 该方法具有简便、实用、精度较高的特点。 本文所作的主要工作如下： 一、基于加权余量原理，针对在空间域进行有限元离散后的一组常微分方程组， 在时间域上利用拉格朗日线性插值离散和有限元单元技术，推导了基于一阶和 二阶系统的时域有限元的递推格式。 二、 在上面得到的理论基础上，推导出了求解一维动力问题的递推格式的时域有 限元算法，并对其进行了稳定性分析。 三、 基于相同原理，推导出了求解二维抛物型瞬态传热问题和二维双曲型传热问 题的递推格式的时域有限元算法。 四、 对动力学问题和传热问题分别进行了数值验算，数值结果表明，本文的计算 方法能较好的反映出问题的特性，是一种适合用于解决经过空间域离散后的常 微分方程组的有效方法，具有可靠的计算精度和工程应用价值。 关键词：加权余量，时域有限元，弹性动力学，传热 基于加权余量原理的时域有限元分析 a b s t r a c t t i m e d e p e n d e n tp r o b l e m sp l a ya l li m p o r t a n tr o l ei nm a n ye n g i n e e r i n ga s p e c t s s u c ha s e l a s t i cd y n a m i c s ，h e a tc o n d u c t i o ne t c s t u d yo nt i m e - d e p e n d e n tp r o b l e m si s d r a m a t i c a l l y s i g n i f i c a n tt h e o r e t i c a l l y a n d p r a c t i c a l l y am e t h o do f s e t t i n gu pn u m e r i c a lm o d e l si sp r e s e n t e di nt h i sp a p e rb ya p p l y i n ga d a p t i v e s p a c e d i s c r e f i z a t i o n a s e r i e so fr e c u r s i v ef o r m u l a sa r ed e d u c e dt h r o u g hac o n v e n t i o n a l w e i g h t e dr e s i d u a lt e c h n i q u ei nt i m ed o m a i n c o n v e n i e n c e 、p r a c t i c a l i t ya n dh i g hp r e c i s i o no f t h i sm e t h o da r ee x h i b i t e db yn u m e r i c a l e x a m p l e s m a i nw o r ko f t h i s p a p e r i n c l u d e s 1 、a p p l y i n g ac o n v e n t i o n a lw e i g h t e dr e s i d u a lt e c h n i q u e ，t oas e r i e so f o r d i n a r yd i f f e r e n t i a l e q u a t i o nt h r o u g h t h es p a t i a ld i s e r e t i z a t i o n ，u t i l i z i n gt h ed i s c r e t i z a t i o no f l a g r a n g ei n t e r p o l a t i o n a n dt h et e c h n i q u eo ff i n i t ee l e m e n tm e t h o d ，d e d u c e dt h er e c u r r e n c ef o r mo ft h et i m ef i n i t e e l e m e n tm e t h o da b o v et h ef i r s t - o r d e ra n ds e c o n d o r d e rs y s t e m 2 、b a s e do nt h et h e o r i e sm e n t i o n e da b o v e ，d e d u c e dt h et i m ef i n i t ee l e m e n tm e t h o da b o u t t h es o l u t i o no f t h e o n e - d i m e n s i o n d y n a m i c s p r o b l e m a n d d o t h e r e s e a r c h o f s t a b i l i t y 3 、b a s e do nt h es a m et h e o r i e s ，d e d u c e dt h et i m ef i n i t ee l e m e n tm e t h o da b o u tt h es o l u t i o n o f t h e p a r a b o l i ca n dh y p e r b o l i c h e a tt r a n s f e r p r o b l e m s 4 、ac o u p l eo f n u m e r i c a le x a m p l e so i ld y n a m i cp r o b l e m sa n dh e a tt r a n s f e rp r o b l e m sa r e g i v e na n dt h en u m e r i c a lr e s u l t ss h o wt h a t t h i sm e t h o dc a r lb eu s e dt o s o l v et h eo r d i n a r y d i f f e r e n t i a le q u a t i o n sw i t hs p a t i a ld i s c r e t i z a t i o n ，n u m e r i c a ls o l u t i o n ss h o wi t sr e l i a b l ep r e c i s i o n a n d g r e a ta p p l i c a t i o n s1 1 2e n g i n e e r i n g - k e y w o r d s ：w e i g h t e d r e s i d u a l t e c h n q u e ；t i m e f i n i t e e l e m e n t m e t h o d ；e l a s t i c d y n a m i c s ；h e a t h - a n s f e r n 独创性说明 作者郑重声明：本硕士学位论文是我个人在导师指导下进行的研究工作 及取得研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文 中不包含其他人已经发表或撰写的研究成果，也不包含为获得大连理工大学 或者其他单位的学位或证书所使用过的材料。与我一同工作的同志对本研究 所做的贡献均己在论文中做了明确的说明并表示了谢意。 作者签名：蕴：塑日期：竺互：i ：7 大连理工大学硕士学位论文 1 绪论 1 1 时域有限元的发展状况 1 1 1 问题的提出 工程中存在大量与时间有关的问题，例如热传导或者弹性动力学这类的瞬态场和动 态问题，这类问题也可以采用与静态分析中的有限元法相类似的方法，是结构离散化， 即把整体结构离散为有限个连续单元结构的组合体，这样，无限自由度系统的问题就变 成了有限个自由度的问题，从数学角度看，就是把偏微分方程的求解问题转化为代数方 程组的求解问题，最后用矩阵分析理论和计算机工具进行各种结构的动态分析。 动态问题和静态问题的主要差别在于动态问题具有随时间而变化的性质，即出现了 时间变量，由于时间变量的影响，显然动态问题不能像静态问题那样具有单一的解，而 必须建立相应于反应时程中感兴趣的全部时间的一系列解答，因雨增加了分析问题的难 度，并且出现了静态问题所没有的特殊问题，例如：特征值分析、阻尼计算等问题。目 前动态问题的分析方法，对空间域的离散与静力学中的空间离散没有多大差别，而对时 间域的离散是动力分析中的新问题。 对于动态问题，在空间上的处理大致上可以分为四类【lj 。即直接法，差分法，加权 余量法和变分法。直接法一是一种物理近似，是一种古典的近似方法，只对结构空间进 行离散，主要用于杆系结构。差分法一是一种数学上的离散近似方法，它将微分化为差 分，然后直接对运动微分方程进行处理。如弹性薄板的弯曲问题，由于弹性薄板的边界 常常是比较规则的，因此，常用差分格式处理该问题。加权余量法一也是一种数学上的 离散近似方法，它可以直接从微分方程中求得近似解。变分法一是研究泛函取极值的方 法。它从能量观点出发，运用变分方法得到能量变分原理。能量变分原理有势能原理， 余能原理和广义变分原理等。 利用有限元法仅对空间域进行离散而不对时域离散，可以得到目前常见的动态问题 的运动方程： 阻如 + c 】f 毋卜医肛 ： p tjlj 基于加权余量原理的时域有限元分析 其中i m i 、i c 、i k l 分别是结构体系经过空间有限元离散后的己知矩阵，妙 r 、r 、 、 u 、 u 分别是有限元集合体的变量、速度和加速度向量， p 则是等效右端项向量。 ljlj 方程( 1 1 ) 是通过考虑在时刻t 的静态平衡而推导出来的，则方程( 1 1 ) 原则上可 以认为是作用体系在t 时刻的静态平衡【2 】表达式。 从数学的角度看，方程( 1 1 ) 是一个二阶线性微分方程组，原则上可以用求解常系 数微分方程组的标准过程求得方程组的解。但是，如果矩阵的阶数很高，则采用一般微 分方程组的求解过程要付出许多艰辛，除非质量矩阵、阻尼矩阵和刚度矩阵有特殊的性 质可以利用。目前求解方程( 1 1 ) 的方法可以分成两大类：振型迭加法和直接积分法。 振型迭加法即先求出结构的一部分低阶特征对，然后将求出的振型作为广义坐标的基 底进行坐标变换，实际上就是一种r i t z 变换，由原来的物理空间转换到模态空间去求解， 这样将体系转化为若干个单自由度体系的组合问题。显然，该方法的计算精度取决于所 假设的模态，并且振型迭加法并非总是有效的。直接积分法一是指在进行数值积分前没 有像振型迭加法那样把方程进行另一种形式的变换，思想是：( 1 ) 只是在相隔f 的一 些离散的时间区间上而不是试图在任意时刻t 上满足方程( 1 1 ) 。即静态平衡是在求解 区间上一些离散时刻点上获得的，因此，似乎在静态分析中使用过的所有求解方法，在 直接积分法中也能有效的应用。( 2 ) 假定位移( 温度) 、速度、加速度在每一时间区间 出内都是变化的。该方法的特点是对时域进行离散，即将方程离散为各离散时刻的方程。 这样，方程( 11 ) 就转化为一个由位移组成的在该离散时刻的线性代数方程组，然后利 用初始条件，可逐步求出在任一系列离散时刻上的动力响应值。而直接积分法有中心差 分法、h o u b o # 法、w i l s d n0 法和n e w m a r k 法。 时间有限元就其概念上来讲与常用的空间有限元一样，就是在我们感兴趣的时间域 内用有限元来进行离散，时域有限元与空间域有限元一样都属于有限元法的范畴。紧随 着在空间域划分单元，在时域内应用有限单元的思想也就开始了，但它不依赖于空间有 限元独立地被研究发展，时间有限元法能得出有高阶精度数值稳定的算法，一开始在实 际应用中受到一定限制主要是因为计算机技术的落后，而现在随着高速高存储计算机技 术的发展，时间有限元法己被开始应用在热问题、弹性动力学及流体力学等许多领域中， 时域有限元算法形成步骤与空间有限元类似，最后得到的是所有变量同时得到求解的算 法，其中有的由于系数矩阵构成的特殊可形成递推算法。 有限元法发展至今，已经成为工程数值分析的有力工具。特别是在固体力学和结构 分析的领域内，利用它已经成功地解决了一大批有重大意义的问题，而根据结构力学与 大连理工大学硕士学位论文 最优控制模拟理论【3 1 ，空间坐标有限元的方法也可用于时间坐标。工程中大量与时间有 关的现象例如瞬态传热问题或弹性动力学这类的瞬态场和动态问题，在空间用有限元法 离散后导致一组关于时间的常微分方程组从而问题转化为常微分方程组初值问题的求 解。 近些年来，人们又开始探讨和研究一种求解方程( 1 1 ) 的新途径，即从平衡方程的 积分形式入手，应用变分原理，在由空间和时间张成的四维广义空间上，然后类似建立 常规的空间有限元模型一样，对时间域进行离散，建立起时间有限元模型。 1 1 2 时域有限元的发展 当对空间和时间域同时离散，即采用同时包含空间坐标和时间坐标的形函数，得出 所谓的时间一空间有限元模型，对结构动力学问题，时空有限元可看作通常的空间有限 元法在时域的一种扩展，只是单元又增加了一个时间维，如梁单元变成矩形单元、板单 元就变成三维棱柱体等。位移插值函数的一般形式为： u ( x ，”乙t ) = n ( x ，y 而f 过 ( 1 2 ) 但是由于这种方法需要针对不同的单元模型，不具备通用性。因此，通常不采用根据时 间坐标和空间坐标同时确定的形函数来描述整个空间一时问域的有限元离散化模型，其 中主要的原因还在于：在一般情况下，同时对空间、时间域进行离散所涉及到的变量太 多，仅就几何性质而言，时间域相对比较简单，使得采用此方法的可行性很小。而 且z i e n k i e w i c z 证实了这种对空间一时间域同时进行离散所得到的结果与分别在空间域 和时间域进行离散所得到的结果基本一致1 4 】，所以将空间域和时间域相继进行单独离散 的方法更比较符合实际发展趋势。 有限元法一般都是建立在极值原理、变分方法或加权余量法基础上的，建立在极值 原理上的方法在某种意义上对真实解有良好的逼近，然而对初值问题没有极值原理可以 应用；有限单元法最初是为了近似求解边值问题而发展起来的，因此试图用有限元法这 种近似的方法解决初值问题也十分自然。z i e n k i e w i c 2 4 】和o d e n 5 均对此进行过总结。早 先的有限元法大多基于变分原理，用于边值问题的变分原理有很多，面用于初值问题的 就很少。 目前关于对空间域及时间域相继单独进行离散的变分方法中，所采用的变分原理主 要是h a m 2 t o n 原理和h a m i l t o n 定律。肌，和b a r u c # 6 1 则将时间域ht ，l 划分成若干个 基于加权余量原理的时域有限元分析 时间区段，在时间区段k 。，t ，l 上用三次h e r m i 抬插值函数逼近广义坐标。并以三次 h e r m i t e 插值函数的二阶导数遥近广义坐标的变分，然后对h e r m i t e 定律的变分式进行 变分计算，得到f 。和f i 时刻的位移，速度及其变分为参变量的表达式，然后按照有限元 方法集成整体方程。并类似在有限元法中引入支撑条件那样，引入了初始条件，同样建 立了一种时间有限元法。但h e r m i t e 原理考虑的时时间域的边值问题，然而实际的动力 学问题常常只能预先给出初位移分布，而不能给出末位移的分布。 我们知道的有限元法一般都是建立在极值原理基础上的，建立在极值原理上的方法 在某种意义上对真实解有良好的逼近，然而对初值问题就没有极值原理可以应用，从目 前情况看，针对边值问题的变分原理很多，而针对初值问题的变分原理却是寥寥无几。 g u r ；i n 型 7 】i s 变分原理的泛函形式及证明较为复杂，但它被认为是目前唯一能表征有关 初值问题的全部特征的变分原理，它通过卷积理论将运动平衡方程和初始条件合并成一 个场方程从而将混合初边值问题转化为等价的边值问题，不过用其解决问题的工作还比 较少。由于泛函中有双重卷积给时间离散带来很大困难，所以近几年才逐渐有所发展， 最厶w i l s o n 和胧c k e l l 首先将g u r t i n 型变分原理用于瞬态热传导问题的有限元分析中 ( 9 l 。之后g h a b o u s s i 和w i l s o 一：o i 又以g u r t i n 型交分原理为基础建立了多孔液体饱和弹性 体动力分析的有限元格式。在以上的分析中，采用的方法与通常的瞬态阀题的有限元法 相同，即首先用g u r t i n 型变分原理建立离散体的平衡微分方程，然后用逐步积分法进行 求解。正如g h a b o u s s i 等人”叫讨论的那样，采用修正的h a m i l t o n 原理可以建立起相同有 限元格式。这种应用g u r t i n 型变分原理的方法由于以新的泛函o ( 0 通过q = t + 0 来 定义) 代替原有命题的泛函q ，致使在最终的微分方程中不再含有卷积积分。也因此消 除了g u r ：i n 型变分原理初值条件这一突出的特点，马立明等【i l 】通过在空间构造l a p l a c e 泛函，得到了几个具有单重卷积的与最小势能原理及修正h e l l i g e r - r e i s s n e r 广义变分 原理相对应的g u r t i n 型和有位移、应力两类变量的g u r t i n 混合型变分原理泛函。而在 文献 12 】中基于位移型g u r t i n 型变分原理，对时域进行离散时，采用具有非时间步参数的 插值函数逼近广义节点位移，给出了一种新的计算弹性动力学初值问题的非时间步参数 的无条件稳定的时间有限元法。 我们已经知道，能量变分求解和微分方程求解两者在一定条件下是等价的，它们是 同一个问题的两种不同表现形式。但是微分方程的形式更广泛些，与其相应的泛函有些 现在还未找到，有些则可能不存在。所以能量变分出现后不久，就有人从微分方程出发 来寻求级数形式的近似解析解，这就是加权余量法( n 】，它在空间域上是求解线性、非线 性微分方程的有效方法，是基于等效积分格式的近似方法也是通用的数值计算方法，它 大连理工大学硕士学位论文 可以用于广泛的方程类型，通过采用等效积分的弱形式，可以降低对近似函数的连 续性要求，从而可以在更广泛的范围内选择试探函数，当近似函数满足连续性和完备性 要求、试探函数的项数不断增加时，近似解可以趋近于精确解n “，它是着眼于区域的整 体，所以对于矩形区域及简单边界条件的二维偏微分方程问题都做出了成功的范例，而 当边界条件复杂或者区域的结构复杂时，我们把区域划分成很多小块( 有限单元) ，然 后对每- 4 , 块做计算，最后总体合成，那么任何复杂问题在小块中都会显得很简单，甚 至不必再去精心选取试探函数，而用最简单的线性插值函数就可得到足够的精度 1 5 】。在 固体力学中的加权余量法，具有原理统一、方法简便、灵活多样、适应广泛、余量可知、 工作量少、计算精度高等特点，特别是将该方法和计算机相结合，程序设计简短，占用 内存不大，数据准备少，计算速度快。因而在国内有了广泛深入的发展，形成了计算力 学的一个新分支，引起国际计算力学界的极大兴趣和关注，但是利用加权余量原理基础 上的时域有限元法进行初值问题的研究，还是一个比较新型的领域。 而运用加权余量法构造时间有 艮元的优点是可用于一般的动态方程，概念明确，应 用灵活，一般的都是假定所研究的全域相当于一个时间单元，在单元内对位移或者惯性 力或者阻尼力以简单的线性、高阶拉格朗目插值、埃米特插值或者三次样条插值作为试 函数，然后由权函数与运动方程相乘并在时间区间上积分得到逐步或者递归的计算格式。 z y e n k i e w i c z 已经提出过基于加权余量原理的时域有限元 1 q ；他首先给出了基于加 权余量形式的时间坐标的有限元，对位移使用拉格朗臼线性及高阶插值给出了单步和多 部递推格式，而后又以t a y l o r 公式逼近响应得出一族单步方法 17 】 1 8 】。我国的蔡承文等 1 9 i 2 0 在更广泛的意义上对位移以任意k 2 次的样条函数逼近，在k = 3 时得出与 e l e n k i e w i c z 等价的多步和单步格式，并表明这两种格式在本质上时互相联系的；这类 在位移空间用一个权函数建立的递推格式不会导致算法矩阵规模的扩大，但无条件稳定 的算法精度不会超过二阶，而在状态空间中用两个或更多的权函数推导公式会导致系数 矩阵规模的扩大却能得到无条件稳定的高阶精度的算法，研究表明时间有限元法有以下 特点： 1 、在一个时间段内引入时域有限元，可以对所有时空变量同时求解进行并行算法或 运用子结构算法。 2 、理论上高次多项式作为基函数可获得任意阶精度，实际可用的算法与基于差分方 法的时间积分方法相比有相对较高的精度。 3 、时间域的“几何”性质简单，所以无须处理复杂的边界。 4 、当问题是激励和材料特性在时间上变化很快时，方法尤其精确。 基于加权余量原理的时域有限元分析 5 、可用于非线性问题，而且能导致更多的可能性( 2 l 】。 1 2 本文求解的问题 对于形如方程( 1 ，1 ) 的问题，可以是结构动力学问题，也可以是双曲型瞬态传热问 题，而当i ml _ 0 的时候也可以是对抛物型瞬态传热问题的求解。 在瞬态传热问题的求解中，一般需要求解与时间有关的偏微分方程。由于实际问题 的复杂性，即使在空间用有限元等方法离散后，在时域也还是要做离散求解。采用高等 数学方法分析热传导过程，一般要假定所研究的对象是个连续体，即认为所研究的对 象内各点上的温度、密度、速度等都是空间坐标的连续函数。实际上，只要被研究物体 的几何尺度远远大于分子的平均自由行程，连续体的假定即可成立口2 】。就本文涉及的内 容而言，连续体的假定都是成立的。 当稳态导热方程加上a 丁a 项以后，就变成了抛物型方程。此时除了边界条件必须 已知外，初始条件也必须是已知的。通常称它为初边值问题。在本文算法中，把要研究 的时间域划分为有限个子域，而把每个子域当作一个时间单元来求解，求解就从初始温 度场开始：每隔一个时间单元，求解下一个单元的温度场，这种递推格式算法的优点是 可以节省计算机内存；此类问题的一般求解方法的特点是在空间域内用有限单元网格划 分，而在时间域内则用有限差分网格划分。实质上就是用有限单元和有限差分的混合解 法。因为它充分利用了有限单元在空间域内划分中的优点和有限差分法在时间推进中的 优点，所以它算是一种成功的结合；但是这种方法的缺点是稳定性较差，容易产生振荡， 甚至发散；在文献 2 3 】中提出将小波理论应用与求解热传导方程，作者将初始波动状态称 为“热小波”( h e a t l e t ) ，在2 空间，热小波可以看作是原子，这样的话任何的热传导 方程都可以用热小波通过一些简单规则“组装”来表示，而本文提出的算法在解决抛物 型传热方程时，在空间域用高精度的八节点等参单元进行离散，而在时间域用加权余量 原理构造时间有限元，当时间步取到较小时，即便用线性拉格朗日插值也可得到精确度 较高的结果。 而当传热问题方程中加上a 2 t & 2 项后，方程就变成双曲型传热问题，与抛物型传 热问题相比，双曲传热问题由于其存在温度场对时间的二阶导数而相对要困难的多。双 啦传热问题是一个有着广泛实际应用前景的交又研究领域，是一个涉及到传热学、物理、 数学、计算机、实验技术等学科的交叉领域。它与许多实际工程相联系，例如：激光在 生物医学中的应用，激光在生物医学中的应用已经有三十多年的历史，而由于这些应用 大部分与热效应有关，因此吸引了不少工程热物理领域的专家加入到这个研究行列中来， 大连理工大学硕士学位论文 一般认为激光照射到生物组织表面后，只有约百分之几的能量被表面反射回来，剩下的 大部分能量进入生物组织内部进行传播。在传播过程中，不断发生散射和吸收。而被吸 收的那部分能量就转变为热能，用于加热生物组织，所以激光照射生物组织的时候所生 成的热能实际上遍布于生物组织的各个角落，是一个容积吸收过程，所以可以将激光照 射处理为传热学中的第二类边界条件 ”。激光制作半导体材料的淬火，低温工程，稀薄 材料的高速传热，以及高强度的电磁辐射等问题。 在这个领域中活跃着许多数学、机械、力学方面的学者，特别是近二十年来，随着 计算科学的迅猛发展，计算速度的大大提高都为它的进一步发展带来了契机，对于这类 问题的求解，许多学者已经做了大量的工作p s i 2 6 1 p t 。对于求解双曲传热问题，许多学者 进行了大量的探讨和研究，文献口8 】是用拉普拉斯变换和控制元相结合的方法对一维双曲 传热问题进行近似求解，为了得至更精确的估计，作者将要研究的整个时间域分成许多 子域来分析；文献【2 9 】发现双曲传热问题中初始温度分布舻和源项厂对于时间域r 的变化 的影响跟温度的初始变化率有关，这个发现揭示了温度域的结构并且显著简化了双曲传 热方程解决方法的发展；而文献口o 中提出用时域精细算法来对二阶双曲传热问题进行求 解；这些文献都为本文算法的形成提供了思路和重要的理论依据，而在这些求解双曲传 热问题的方法中有限单元法是被广泛采用的数值方法之一，在空间域经过有限元离散后 的双曲传热方程，是一组关于时间的常微分方程组，从而问题转化成为常微分方程组初 值问题的求解。 在本文中利用加权余量原理将初边值问题转化为系列具有递推形式的边值问题， 从而可以方便地利用各种求解边值问题的技巧进行求解，这是可以对所有变量同时得到 求解的算法，在本文算法中利用加权余量原理对时间域进行离散，并以拉格朗目插值函 数作为试探函数，然后由权函数与运动方程相乘并在时间区间上积分得到逐步或者递归 的计算格式。而本文则是将整个时间域划分为有限个类同子域，然后将子域作为一个时 间单元，通过子域来描写整个时间域。 1 3 本文所作的主要工作 本文采用高精度的八节点等参单元，建立一般的动力学、瞬态传热问题的数值正演 模型。针对动力学问题和瞬态传热问题中存在的关于时间场的导数的问题，本文主要做 了以下工作： 一、借助恰当的空间离散方法，对空间域进行离散，即在空间将模型离散成为有限个单 元体；在本文中对空间域的离散采用的是有限单元法，利用高精度单元( 如八节点 基于加权余量原理的时域有限元分析 等参单元) 对模型在空间域进行离散，并考虑到各类边界条件和初始条件。建立便 于求解的正演模型。 二、空间坐标有限元的方法也可用于时间坐标，离散时间坐标后，在典型单元内即可通 过加权余量原理得出平衡方程，并在空间上将全域化分为有限个类同子域，而将单 个子域作为一个时间单元，通过得到的逐步计算格式描述子域，最后得出基于时域 的递推格式。给出一种较新型的时域有限元法，并对本文方法进行了稳定性方面的 研究分析，本文方法具有计算理论简单可靠，计算方法简便实用、计算结果精度较 高等特点，是适用于动力学问题和瞬态传热问题的有效方法。本文的工作，为加权 余量原理的应用方面开拓了一条新的途径。 三、在上面获得的理论基础上，针对以上工作，编写相应的动力分析和传热分析程序。 对一些动力学问题、抛物型传热和双曲型传热问题分别进行了数值验算，验算结果 与其他算法进行对比，对比结果表明，本文的计算方法能较好的反映出问题的特性， 是一种适合用于解决经过空间域离散后的常微分方程组的有效方法，具有可靠的计 算精度和工程应用价值。 大连理工大学硕士学位论文 2 动力学分析以及时域有限元的方法介绍 2 1 引言 运用加权余量法构造时间有限元的优点是可用于一般的动态方程，概念明确，应用 灵活，原理简单等优点，首先对问题进行空间离散，本文使用的高精度有限单元( 如八 节点等参单元) ，针对得出的常微分方程，运用加权余量原理，将初边值问题转化成一 系列具有递推形式的边值问题，而在本文算法中介绍的时域有限元中使用的加权余量法 同传统的加权余量原理有所不同，传统的方法一般都是假定所研究的全域相当于个时 间单元，而在本文算法中则是将要研究的整个时间域划分为有限个相等子域，而将子域 作为一个时间单元，在单元内运用加权余量原理建立平衡方程，给出一种时域有限元方 法，并且研究了该时域有限元方法的稳定性问题。而对于时间域，由于时间域的“几何” 性质很简单，所以根本无需处理复杂的边界，而且当时间步比较小的时候，即便对位移、 阻尼力或者惯性力以简单的线性拉格朗曰插值函数作为试探函数也可获得较高的精度， 然后由权函数与运动方程相乘并在时间区间上积分得到逐步或者递归的计算格式。 2 2 动力学问题 22 1 弹性动力学基本方程 设有弹性体n 其边界为量v 的闭包矿= v u s 。又s 可分为两部分，即给定位移 的边界s 。和给定面力的边界s ，且有： s = s 。u s 。；s 。n s 。= 0 ( 2 1 1 ) 弹性体所占的时间一空间域用v x 0 ，c 。) 表示。各场变量均为空间坐标j 以及关于 时间t 的函数，其中x = g ，x ：，_ ) ，且应满足： 动力平衡方程 盯+ 工= p 矿0 2 u i ( 在y 。，o 。) 上) 应变一位移关系 ( 2 1 2 ) 基于加权余量原理的时域有限元分析 白：委+ “) ( 在y x o ，。) 上) 应力一应变关系 。b2 e 一“ 或= c f k t o “ ( 在矿 0 ，。) 上) 式中，e 州和c 州分别是体系的刚度系数和柔度系数。 力的边界条件 ( 2 1 3 ) ( 2 1 4 ) l = ”= 霉( 在s ，x o ，o 。) 上) ( 2 1 5 ) 式中，为边界面外法线方向单位矢量，霉是应力边界s ，上给定的面力。 位移边界条件 暖是位移边界s 。上给定的位移。 初始条件 ( 在& 0 ，m ) 上) g ，o ) = 玩，g ) 一g ，o ) = 瓦，g )( 在v x o ，m ) 上) 以上各式构成了弹性动力学的混合初值一边值问题。 2 22 算法步骤 本节就一个简单的杆问题进行讨论； 问题的基本方程( 考虑阻尼影响) ： ( 2 1 6 ) ( 2 1 7 ) 大连理工大学硕士学位论文 瓦o c t + 矗= p 矿0 2 u + c 詈c d r 甜 其中p 、c 为已知参数 对方程( 2 2 1 ) 使用加权余量法，则有 j o i o v 瓦o 仃+ 占一p 矿0 2 u c 如= 。 其中：v 为权函数向量，盯= e s ，s = 罢。 o x 分部积分后式( 2 2 2 ) 有： 吼。呻卜始卜 p e v 7 鲁也一c e v 7d u + p b d u = 。 有限元离散，即 “= 坚v = ! 其中 】= i1 一手，引 “= “o + 亭 d u = 蟛， f 为空间杆单元长度， “。为杆单元初始位移。 则式( 2 2 3 ) 成为： ( 2 2 1 ) ( 2 2 2 ) ( 2 2 3 ) ( 2 2 4 ) 基于加权余量原理的时域有限元分析 ( 1 7 i ) ，：k 一也7 i ) 。b 一羔7f n r e n d 考u 一引n t c 嘴鼍_ v r l n to n d l 等十蛐n t b d = q 即得到在空间已经过有限元离散的离散方程组 m i i + c u + k u + f = 0 ( 2 2 5 ) ( 2 2 6 ) 其中“为待求变量，m ，c ，足为已知的经过有限元空间离散后得出的质量、阻尼和刚 度矩阵，厂为考虑初始条件的已知右端项向量。 对式( 2 2 6 ) 应用加权余量法有： w 7 ( m i i + c u + k u + f ) d r = 0 其中w 为权函数向量。 分部积分后式( 2 2 7 ) 有 ( 2 2 7 ) c w 7 埘牲一c w 7 谢牝。一j ! ( 害 7 m 鲁d r + c w 7 ，k c w 7 c “，k 一! ：| 警丁函d r + j ：1w 7 觑d r + r w 7 仁。 c z z s ， 现将向量“进步在时域上做有限元离散，即 “= n uw = n w 则式( 2 2 8 ) 成为 ( 2 2 9 ) 大连理工大学硕士学位论文 ( 型7 m n ) 1一( 型7 m ) l一旦f 7 m n d r u + ( 型7 c h ) 一( 型7 c u ) l l f = “i f = t o o l 扫t n jr - f o 一型7p 7 c n d r u + w 7 7 k d c d c u + w 7 厂= 0 注意w 是将w 也按与“相同的时间有限元离散后的权向量。 于是式( 2 2 1 0 ) 可写为： w 7 ( 丝+ 一整) 堑一w 7 ( 胁+ c ) l= 丝7 f 一型7 ( m u + c “) i l 扣f 1 f = t o 将丝重新排列，t f 是按照时间排序，重新排列即是要化为坚，转化成按照空间排序， 即：原来的“的排列顺序是第一时刻一节点、二节点至第t 2 节点，转化以后变成丝的顺 序是第一节点的一时刻、二时刻至第n 时刻，重新排列“的同时也将型7 重新排列后，原 方程就化为； 一w “g u 一1 4 2 “f 7 曩刘 = 一w l ； i 其中丛= p 7 m n d f c = p 7 c n d f 丝= 7 k n d f f ? = 妒f d r 婪7 是位移在时间域的排列展开： 型7 = ( 2 2 1 2 ) 基于加权余量原理的时域有限元分析 消去w “，于是式( 2 2 1 2 ) 成为 r 埘o + c u o 、 啦7 ：，一 j l l 一( m f i “+ c u “) j 其中上标的0 、 分别表示一个时间单元内的初始节点和最后节点。 将式( 2 2 1 5 ) 分块为： 雕鳓制一 一篡， 凝聚“” “= g ：c f “一g ：。“。一g ：。“) 则由( 2 2 1 6 ) 可得 g i 。“。+ g i ：g 丢p “一g ：。“。一g 2 3 un ) + q 。“= f 。一 叛。一函。 ( g 1 ，一g l ：g ；：g 2 ，l “= f 。一g l ：，“+ ( g l ：唰g 2 。一c g 3 , 。一蚴。 继而得出 ( 2 2 1 5 ) ( 2 2 1 6 ) ( 2 ，2 。1 7 ) “8 = ( g ，一g 1 2 g 丢g ：，厂1 p 。一g 1 2 g ：- ：1 f ”+ ( g 。：g 丢g 。一c g 。- 。一慨。) = 曩瓦+ 2 “o + ，i 。( 2 ，2 1 83 进而得出 1 4 大连理工大! 学硕士学位论文 m u ”= 一( f ”一g 3 l “一g 3 2 “一g 3 3 “4 + c u “) = 一( f ”一g 3 ，“。一g 3 ：g ：p “一g 2 ，“。一g 2 3 u ) 一g l “”+ c u “) = 一( f ”一g 3 。“。一g 3 ：f “+ g 3 ：g 2 。“。+ ( g 3 ：g ：g 2 。一g 3 ，+ c k ”) = 一( f ”一6 ：g 丢f ”+ ( g 3 ：g 墓g 2 ，一g 3 ，- 。+ ( g 3 ：6 葛g ：，一g 3 ，+ c b “) = 一f ”+ g 3 ：吲f “一( g 3 ：g 丢g 2 一g 3 。- 。 一( g 3 ：q - ：1 g 2 ，一g 3 3 + c ( j i z l 局+ h 2 u 。+ h 3 n 。 = 一( f “一g 3 ：呸- ：t ，i n + ( g 3 ：呸- ：1 g 2 ，一g 3 3 + c h r ) 一( g 3 ：g 2 1 + 帜：g 2 3 一q ，+ c ) i z ：一g j 。- 。一心：g ：，一g 3 3 + c ) ( 2 2 1 9 ) 于是解出： i “= e + h 4 u o + 5 i o 其中： 只：m 。f ，一一g 3 ：g 丢f m + ( g 3 ：吲，一g 3 3 + g 2c h 民1只= m 。i ，”一g 3 ：g 丢f “+ 【g 3 ：吲，一c 虹只l h 。= 埘1 ( g 3 ：g _ ：g 2 ，+ ( g 3 ：g ；g g 2 ，一g 3 ，+ c k ：一g 3 。) h ，= 划1 ( g 3 ：g 丢g ：，一g ，+ c k ； 联立( 2 2 1 8 ) 与( 2 2 ，2 0 ) 式，得： ( 州鸶h i 上述推导过程对二阶、一阶的方程均适用。 1 5 ( 2 2 2 0 ) ( 2 2 2 1 ) ( 2 2 2 2 ) ， 踟砒 v 几 基于加权余量原理的时域有限元分析 23 算法的稳定性分析 积分格式的稳定性是指对很大步长值告的方程，一定不要人为的放大其初始条件， 以致于使低阶振型的积分毫无价值。稳定生也意味着：在时刻t 的位移、速度由于计算 机上的舍入而产生的误差不会因为积分而增大。 积分方法可分为无条件稳定和条件稳定两种；如果在任何初始条件下，对于任意步 a + 长a t ，特别是詈很大的时候，解都不会无限的增大或者振荡，则这种积分方法世无条 件稳定的。否则，如果仅当等小于一个通常成为稳定极限的确定值的时候，才不会出现 解的无限增大或者振荡的情况，则这种积分方法是条件稳定的。 种算法的积分格式是否稳定，通常是研究单自由度、无阻尼、自由振动的情况。 看数值解法的传递矩阵的谱半径是否小于l ，若对任意时间步长f ，传递矩阵的谱半径 都小于1 ，则这种积分方法格式是无条件稳定的。若只有当等小于某一确定值的时候， a f 传递矩阵的谱半径才小于1 ，而等大于该值的时候，传递矩阵的谱半径就不小于1 ，则 这种积分格式就是条件稳定的。 下面就是本文方法的稳定性分析； 运动方程组经变换为n 个不相耦合的微分方程后，其性质不变，因此可以方便的对 非耦合的微分方程讨论解的稳定性。又因为月个方程是相似的，故仅需分析其中一个作 为代表，将它写成 i i 十2 j j i j + 磷x j = r 或鼍+ c ：i 。+ 砰一= l ( 2 3 1 ) 讨论解的稳定性实质上是讨论误差引起的响应，所以在上式中可令r = 0 。 另一方面由于在正阻尼的情况下，阻尼对解的稳定性是有利的，所以在讨论解的稳定性 时，总可令e = 0 ，基于上述两点，要讨论的方程是 j ，+ 岔x i = 0 ( 2 3 2 ) 大连理工大学硕士学位论文 利用时域有限元法对( 2 3 2 ) 进行求解，可以写出 ( _ l 一华a 。+ ( 2 一半a + ( _ 1 一华a 矿。 ( 23 。3 ) 假定解的形式为 ( a i ) ，+ r = 咒( 以f ) ，( 口f ) ，= 兄( q ) 卜出 将( 2 3 4 ) 代入( 2 3 2 ) ，则可得到特征方程 ( 一6 一p ，) 2 + ( 1 2 8 p i ) 2 + ( 一6 一a ) = o 其中p 。= 纽2 砰 解出上式的解 ，一1 2 + 8 p f + ( p 。一2 ) 2 4 2 可i i 厂 ，一1 2 + 8 p f 一0 。一2 ) 2 4 如2 可i i 厂 ( 2 3 4 ) ( 2 3 5 ) ( 2 3 6 ) a 的根关系到解的性质，首先为使在小阻尼情况下的解具有震荡特性，a 必须是复 数，这就要求： ( a 一2 ) 2 4 0 因为p 。= a t 2 砰，所以 p 。 4 同时= 2 7 r t ，所以上式的到 ( 2 3 7 ) 基于加权余量原理的时域有限元分析 r 至 万 其次为了使问题的解不会无限地增长，还应要求 州1 则( 2 3 6 ) 表示的 和如在只 4 情况下，自动满足川1 要求。 因此为保持解的稳定性，时间步长必须服从以下条件 出0 ：玉 石 ( 2 3 8 ) ( 2 3 9 ) ( 2 3 1 0 ) 其中，是临界时间步长，z 是系统的最小固有周期。实际上并不需要从求解整个 系统的固有值问题得到。关于l 的估计，i r o n s 已证明系统的最小固有周期总是大于 等于最小尺寸单元的最小固有振动周期，因此可以简便地将后者用于( 2 3 1 ) 式，其结 果总是偏于安全的。 2 4 数值算例 1 、一维算例 首先考虑一维( 不含空间) 问题： 运动方程为： 馨+ 6 呈：，t 矿+ 6 百”- ， d = b 2 4 c ，分别考虑下面几种情况 例1 、d ob = 4c = 3 1 3 初始条件为： o = 1 0 u 0 = 2 0 基于加权余量原理的时域有限元分析 该问题精确解为： “一e 一3 。+ e 一+ 1 用本文方法跟精确解进行比较，取时间步为0 5 秒，共计算5 0 步。计算结果与精确 解曲线比较结果见图( 2 3 ) 图2 2 本文数值解和精确解的比较 f i g u r e2 2t h ec o n t r a d i s t i n c t i v er e s e to f n u m e r i c a ls o l u t i o na n da n a l y s i ss o l u