（固体力学专业论文）输流管道在粘性流动下耦合振动问题的数值方法研究.pdf

y8 3 2 7 0 7输流管道在粘性流动下耦合振动问题的数值方法研究固体力学专业研究生邱清水指导教师于建华输流管道耦合振动问题是近五十年柬管道流体力学、非线性动力学等领域研究的热点。由于该问题的复杂性，目前对它的处理主要集中在数值研究方面，但现有的研究方法对该问题采用太多简化，尚未能真实地反映耦合性质。本文深入臻究了鲥方法。构造了基于p e t r o v 一( ；a l e r k i n 格式畿迎风有限元列式。在计算上，流体区域内采用n e w t o nr a p h s o n 迭代求解非线性的n a v i e r s t o k e s 方程组。根据网格扭曲的两种情形一剪切扭曲和压缩扭曲，提出网格运动的局部线性插值与重构方法。所提方法设计简便、易于程序化，同时具有较高的汁算效率。本文还对a l e 有限元在湍流计算中的应用进行了初步研究，分析表明a l e 方法能够方便地被推广到湍流的计算中去。最后，采用面向对象的方法，在c + + b u i l d e r 编译环境下编制了相应的程序。算铡模拟结果表觋，本文方法具有撮好能效率和精度。关键词：流固耦合a l e 有限元网格运动输流管道t h en u m e r i c a lm e t h o d sr e s e a r c h o fc o u p l e dv ib r a t i o no fv i s c o u sf l o wf l u i d c o n v e y i n gp i p e ss o l i dm e c h a n i c sg r a d u a t eq i uq i n g s h u is u p e r v i s o ry uj i a n h u at h eq u e s t i o no fc o u p l e dv i b r a t i o no ff l u i d c o n v e y i n gp i p e si sh o t s p o ti nt h er e s e a r c ha r e a so fp i p e l i n ef l u i dm e c h a n i c sa n dn o n l i n e a rd y n a m i c sf o rl a t ef i f t yy e a r s ，b e c a u s eo ft h ec o m p l e x i t y ，t h er e s e a r c hi sf o c u s e do nn u m e r i c a lp r o c e s sb u tt h ee x i s t i n gp r o c e s sm e t h m sa p p l yt o om u c hs i m p l i f i c a t i o n8 0t h ec o u p l i n gc a n r o tb er e f l e c t e dt r u t h f u l l y b a s e do np e t r o v g a t e r k i nu p w i n df e m ( f i n i t ee l e m e n tm e t h o d ) f o r m a t ，t h ep a p e rm a d e o m er e s e a r c ho na l em e t h o dw h e nc a l c u l a t i o n n e w t o n r a p h s o ni t e r a t i o nw a sa p p l i e dt os o l v et h es e to fn o n l i n e a rn a v i e r - s t o k e se q u a t i o n si nt i l ef l o wa r e a 、a c c o r d i n gt ot w oc o n d i t i o n s ，s h e a r i n gd e f o r ma n dc o m p r e s sd e f o r m ，t h ep a p e rp u tf o r w a r dt h em e t h o do fs e c t i o n a ln u m e r i c a lp l u ga n dr e c o n s t r u c t i o no fm e s hm o t i o n t h i sm e t h o di ss i m p l et od e s i g na n dp r o g r a ma n dh a v eg o o de f f i c i e n c yo fc a l c u l a t i o n ，t h ep a p e ra l s om a d es o m ep r e l i m i n a r yr e s e a r c ho nt h ea p p l i c a t i o no fa l ef e mi nt h ec a l c u l a t i o no ft u r b u l e n tf l o w , w h i c hs h o w e dt h a tt h em e t h o do fa l ec o a lb ec o n v e n i e n t l ya p p l i e dt ot h ec a l c u l a t i o no ft u r b u l e n tf l o wf i n a l l y , a p p l y i n go o p ( o r i e n t e d - o b j e c tp r o g r a m i n g ) ，im a d et h ec o r r e s p o n d i n gp r o g r a mi nt h ec o m p i l a t i o ne n v i r o n m e n to fc + + b u i l d e rt h es i m u l a t i o nr e s u l t so fe x a m p l e ss h o w e dt h a tt h i sm e t h o dh a dv e wg o o de f f i c i e n c ya n dp r e c i s i o n k e yw o r d s ：f l u i d s t r u c t u r ei n t e r a c t i o n ；a r b i t r a r yl a g r a n g i a n e u l e r i a nf i n i t ee l e m e n tm e t h o d ；m e s hm o t i o n ；f l u i d c o n v e y i n gp i p e s四川大学硕士学位论文1 概述输流管道在众多工业领域，如：船舶动力、石油化工、水利工程、航空航天等都有着十分广泛的应用，发挥着极为重要的作用；在生物医学工程研究中，如人体血管系统受力特性和病理学，也涉及到变截面输流管道的力学分析。但输流管道在工作过程中由于流体的流动状态的变化而引起流体、管道问的耦合振动，其动力学行为十分复杂”1 。因此有关输流管道的耦合振动研究一直是理论界和工程界广泛关注的课题。本文主要研究输流管道耦合振动问题的数值计算方法。1 1 本文研究的目的和意义与流体相接触的结构在受到水流波浪和冲击力等流体动力的激励或其他非流体的激励而发生振动时，会使其周围的流场发生变化。这种流场的变化反过来使结构受到的流体动力发生变化，形成反馈的流体一结构相互作用( f l u i d s t r u c t u r ei n t e r a c t i o n ，简称f s i ) 的耦合问题。输流管道的祸合振动是典型的f s i 问题。在这一类问题中，水力暂态。1 诱发管道振动，振动又引起新的水力暂态，这种相互间的作用构成一个具有强烈耦合效应的非线性动力学系统，是近5 0 年来管道流体力学、非线性动力学等领域研究的热点。f s i 是输流管道固有的力学特性，研究输流管道中f s i 现象对了解输流管道系统动力学性态，稳定系统运行、提高运行可靠度、避免管道出现超常压力、应力、振动噪音等具有重要学术及应用价值。变截面输流管道内的流体流动在工业工程和生物医学中具有许多应用背景。在工程上，采用变截面管制作热交换器可以改善其热交换性能。在生物医学中，动脉粥样硬化的形成会使动脉出现局部狭窄，局部狭窄的出现会改变心血管系统的流动与生理特性“l 。介于实验条件的局限，当前对输流管道中f s i 现象的研究主要集中在数值研究方面a 由于研究对象的复杂性，现有的研究仍然有相当多的按照计算流体动力学( c f d ) 的方法处理，将管壁刚性化而不计管壁的变形嘲：有的研究考虑了管壁的变形，采用线性假设，但其处理方法仍然简单化，乃将流体域与固体域分开并单独求解。这二种方法不仅忽略了f s i 的重要特性，而且在计算上将可能由于管壁的位移和较强的非线性使得网格的修改及防止流动畸变出现不少输流管道枉粘性流动下耦合振动问题的数值方法研究的问题，造成较大误差。非线性的数值建模方法在流体和同体各自的领域中取得了可喜的成绩，但对于f s i 问题，遇到的最大困难就在于采用统一的坐标系及两相界面的协调处理问题。为了解决不同坐标系的困难，a l e ( a r b i t r a r yl a g r a n g i a n e u l e r i a n )被引入到流固耦合分析中用以处理界面协调及自由面问题1 。与其他耦合一样，由于非线性使叠加原理失效，使得非线性输流管道耦合振动问题的求解需要探求全场求解方案。综上所述，本文的研究具有。阔的工程背景和重要的理论参考价值。1 2 流固耦合定义和特征流固耦合力学是流体力学与圆体力学交叉形成的一门力学分支。在学科上它还涉及到计算力学、振动力学等学科知识：在工程技术上它与许多工程领域，如十木工程、生物工程、核动力等均有不同程度的联系。流固耦合力学的一个最重要特征是流体与固体两相介质之问的相互作用，即可变形的同体在流体荷载的作用下会产生变形或运动，而变形或运动又反过来影响流场的状态，从而改变流体荷载的分布和大小“1 。正是这种两相介质之间的相互作用将在不同条件下产生形式各异的流固耦合现象。一般而言，流吲耦合问题可由其耦合方程来定义，这组耦合方程的定义域同时具有流体域与固体域。由于在描述流体的控制方程中有描述固体运动的未知变量，在描述固体的控制方程中亦有描述流体运动的未知变量，因此其耦合方程通常具有以下两个特钲：i ) 这组藕合方程在流体域或固体域中均不能单独地进行求解；2 ) 无法显式地从耦合方程中消去描述流体运动的独立变量或描述固体运动的独立变最使方程解耦。按猫合机理的不同，流固稿合问题大体上可以分为两大类。第一类问题的特征是流固两相介质部分或者全部地重叠在一起，例如，渗流问题”1 、多孔介质中的波的传播问题”1 等。在这类问题中两相介质难以明显地分开，使得描述其物理现象的方程，尤其是本构方程需要针对具体的物理现象来建立。第二类问题的特征是耦合作用仅仅发生在两相介质的交界面上，例如，输流管道的耦合振动问题”1 。1 等。这类问题在运动学和动力学上，流固两相介质均有各自的拙述方法和力学特性，其联系仅在于耦合面上的力学平衡及运动协调关系。四川大学硕士学位论文1 3 流固耦合问题的研究方法在流固耦合过程中，流体与固体两相介质交互作用，其非线性特性是非常明显的。不仅流体运动方程是非线性的，而且耦合运动的特性将随着结构振动的幅值不同而变化，这就使得耦合过程的力学行为呈现出十分复杂的特征，其数学描述可能是高度非线性的。由于流固耦合系统的复杂性，当前对这类问题的处理主要集中在数值分析方面，其方法主要有：有限元方法“”1 、边界元法“、边界元与有限元相结合的方法“”“、流体元的有限差分法以及流体元与有限元的联合运用等。与其它方法相比，有限元方法具有对复杂边界条件和几何形状良好适应性和更容易处理复杂本构关系和非线性问题的优点，因此它能够灵活地模拟和逼近复杂的求解区域”。有限元方法首先在飞机结构分析中被提出来，尔后在固体力学中得到广泛的应用，并且逐步渗透到其它领域，如传热学、流体力学、电磁学等计算当中。在过去的五十年中，随着计算机的出现和发展，有限元方法已取得了很多令人瞩目的进展。在国外已经研制了不少大型结构有限元分析软件，著名的如：n a s t r a n 、a n s y s 、a d i n a 、s a p 等。在国内，对有限元法的研究也已取得了丰硕的成果，但在研制大型通用结构分析软件方面，尚还有很多工作要做“”。在有限元计算当中，根据描述物体运动方法的不同，可以归为两大类，即拉格朗日描述方法和欧拉描述方法。其中，在固体力学大变形问题中绝大多数研究工作都使用拉格朗日描述方法，而在流体力学有限元中主要采用欧拉描述方法1 8 3 1 。拉格朗日描述方法又称物质描述或l 描述，它是以初始构形为参考构形来建立物体的运动学和动力学方程。拉格朗日描述的特点是计算网格固定在物体上随物体一起运动，即网格点与物质点在物体的变形过程中始终保持重合，因此物质点与网格点之间不存在相对运动( 或称迁移运动、对流运动) 。这就大大地简化了描述物体运动的控制方程的求解过程，并且能够准确描述物体的移动界面，和跟踪质点的运动轨迹。但在涉及到大变形的问题中，物体的扭曲将可能导致计算网格的畸形而使得计算失败。欧拉描述方法又称空间描述或e 描述，它是以现时构形为参考构形来建立物体的运动学和动力学方程。欧拉描述的计算网格是固定在空间的，即计算网格在物体的变形过程中始终保持不变，因此可以很容易地处理物体的扭曲。但输流管道在粘性流动p 耦台振动问题的数值方法研究对移动界面的处理需要引入非常复杂的数学映射，这将可能导致较大的误差。另外，在欧拉描述中使用普通的g a l e r k i n 离散时，由于迁移项的影响，会导致有限元方程中的系数矩阵是非对称的，增加了求解的困难，而且还可能得到振荡数值解。以上两种运动描述方法卜的差异，使得流固耦合界面上点的协调比起单纯的固体或流体非线性问题更加困难。但同时两种运动描述方法又都有各自的优势，如果能够将二者有机地结合起来，充分利用各自处理问题的优势，克服各自的缺陷，则可解决许多只用纯拉格朗日和纯欧拉描述所解决不了的问题。a l e描述f 是基于此目的而被提出。自其被提出以来，a l e 描述获得了快速的发展，并与有限差分、边界元、有限元等方法相结合，在许多流固耦合问题中得到了应用。目前a l e 描述有限元是研究流固祸合问题有效方法之一。14a l e 描述研究进展141a l e 描述法拉格朗日描述方法和欧拉描述方法是连续介质力学中两种经典的描述连续体运动的方法。前者以连续体的初始构形为参考构形来研究物质点在空间中的运动规律，它描述了同一质点在不同时刻的空问位置；后者则是以连续体的现时构形为参考构形来研究某空间点上物质点的运动规律，它描述了同一空间点在不同时刻被物质点所占据的情况。彳i 同于拉格朗h 描述和欧拉描述，a l e 描述是在连续体现时构形和初始构形之外引入一个可以独立于现时构形和初始构形运动的参考构形。在物体的变形和运动过程中，参考构形是始终固定不动的，而现时构形和初始构形都相对于参考构形运动。a i j e 描述法综合了纯拉格朗f 1 描述法和纯欧拉描述法的优点，克服了各自的缺陷，它在固体区域采用拉格朗日描述，在流体区域采用欧拉描述方法，而在流固交界面上使用混合欧拉一拉格朗日描述法。正是这一在不同的区域使用不同的描述的特点，使得a i 。f ：描述能够准确地描述流固耦合现象。a l e 描述为拉格朗日描述和欧拉描述的联合使用提供了一条有效的途径。它的一个重要特征是计算冽格是基于参考坐标而划分的，也就是说计算网格可以独立于物质坐标系和空间坐标系运动。这样通过指定合适的网格运动形式就四川大学颈士掌位论文一一可以准确地描述物体的移动界面，并使单元在运动过程中保持合理形状，克服了纯拉格朗日描述和纯欧拉描述的缺陷。实际上，纯粹的拉格朗日描述和欧拉描述是a l e 描述的两个特例，即当网格的运动速度等于物体的运动速度时，a l e描述就退化为拉格朗日描述：而当网格固定在空间不动时，a l e 描述就退化为欧拉描述，因此a l e 描述提供了两种描述方法的一种统一的描述。1 4 2a l e 描述研究进展a l e 描述方法最初出现在求解流体动力学问题的有限差分方法中，它是由n o hwf c e l ( 1 9 6 4 年) 以耦合欧拉一拉格朗日( c o u p l e de u l e r i a n l a n g r a n g i a n )的术语提出的，并且与有限差分法相结合，求解一个带有移动边界的二维流体动力学问题“。在n o h 的研究工作中，计算网格点可以随物质点起运动，但也可以在空间中固定不动，甚至计算网格可以在一个方向上固定，而在另一个方向上随物体一起运动，因此a l e 描述也被称为耦合欧拉一拉格朗日描述。虽然n o h 建立了基于可以自由运动的参考构形上的控制方程，但在程序实现时实际上是假定整个计算区域只能使用拉格朗日描述和欧拉描述中的一种，而不能使网格自由运动。1 9 7 4 年h i r tgw ，a m a d e n 从将a l e 描述与求解方程的隐式积分格式结合，解决了具有任意流动速度的二维流体流动问题1 ，并研制了相应的有限差分程序y a q u i 。1 9 7 7 年s t e i nl e l a n dr ，g e n t r yrha n dh i r tcw 将a l e 描述推广到具有任意形状活动边界的三维时间依赖性流动问题中，用以处理三维边界条件，并通过个六面体受水平波作用问题的实例，得到了与实验相符合的有关结论。“。h m s d e nha ，r u p p e lhm 研制了二维、三维简化a l e 格式有限差分程序s a l e 和s a l e 一3 d ，他们在设计程序时认为每个单元只由一种材料组成”m 1 。以上这些研究工作都是将a l e 描述方法嵌入有限差分方法中，并且只考虑了无粘性可压缩流体，解决了部分流固耦合问题。a l e 描述方法被引入到有限单元中，最初是为了满足核反应堆结构安全分析中的非线性数值模拟技术的需要，并用它求解流体与结构相互作用瞬态问题。l i uw k ”铡用a l e 有限元对贮液箱在地震荷载作用下产生的液面晃动、动态弯曲、贮液箱与液体交界面处结构的动力学行为进行了研究，给出了相应的a l e 处理方法以及数值计算格式。l i uw k ，c h e nj s 汹3 基于a l e 有限元和自适应输流管道在粘性流动下耦合振动问题的数值方法研究网格技术对材料加工成型过程进行数值模拟，在他们的研究当中，对流固交界面的处理采用了a l e 描述和外力虚功率联合运用方法，同时给出了一个改进的牛顿迭代公式，讨论了a l e 有限元在非线性问题中应用的策略问题，并给出了金属加工成型的示例。t h o m a sh u g h e sj r 2 7 研究了不可压缩粘性流体流动的问题和带有自由液面液体流动问题，并初步建立了a l e 描述法的一般运动学理论框架。h u e r t aa ，l i uw ”发展了a l e 描述的一般理论框架，相应地推导了在自由液面大幅度晃动下非线性粘性流动的a l e p g 有限元计算格式，还在程序实现时使用预测一多步校正的方法，最后对海啸、坝库相互作用等大r e y n o l d s数大晃动问题进行了研究。a l e 描述的一般理论框架的建立为a l e 的发展与广泛应用提供坚实的理论基础。在l i uw k ，h e r m a nc h a n g ，e t a l “的研究中，对a l e 描述的运动学方程进行了线性化处理，给出了切线刚度矩阵的积分形式，并对算法的稳定性作了初步的研究。l i uw k ，b e l y t s c h k ot ，c h a n gh 3 0 将a l e 描述方法应用于路径相关性材料问题。他们建立了路径相关材料在a l e 描述下有限元计算的一般格式，导出了a l e 有限元显示时间积分算法，并将其理论在一维弹塑性波动动力学问题中加以应用，得到了满意的结果。t a k a s h in o m u r a 、t h o m a sjrh u g h e s o “。t a k a s h in o m u r a ”1 研究了受弹簧支撑的剐体在粘性不可压流体中的耦合运动，在研究中采用a l e 描述方法处理耦合的边界条件、使用s u p g ( s t r e a m l i n eu p w i n d p e t r o v g a l e r k i n ) 方法在流体域上进行空间离散、并运用预测一多步校正方法进行时间步积分。k j b a t h e 等人。”选取速度势和密度作为流体区域的变量，导出了相应的a l e 有限元计算格式。j c k o o 、e p t a h r e n t h o l d 。”从能量观点出发，应用虚功原理建立了可压缩粘性流动的a l e 动力学方程，并直接对哈密顿正则方程进行有限元离散。国内学者在这方面的研究始于上世纪九十年代中期。刘志宏、黄玉盈。m 将a l e 描述法与边界元方法相结合，利用a l e 网格可任意移动的特性，对自由液面的大晃动问题进行研究，实现了控制点的自由跟踪、避免了单元的畸变。曾江红、王照林”“在研究具有自由液面大幅晃动的n a v i e r s t o k e s ( 简称n s )方程的数值求解时，对流体区域采用a l e 运动学描述，使网格自由移动，并采用流线迎风( s u p g ) 加权余量法处理强对流项。温德超等人m 3 用a l e 方法研究三维自由液面问题。岳宝曾、王照林、匡金炉3 将a l e 边界元法应用于非线四川i 大学硕士学位论文性晃动问题，并应用g m e r k i n 加权余量方法对动力学边界条件进行有限元离散，在时间上采用n e w m a r k 方法进行离散，推导了预测一多步校正的计算格式。郑群、刘顺隆。”“0 1 推导了旋转叶轮内流湍流n s 方程计算的a l e 格式，并以亚网格( s g s ) 模式模拟了湍流运动。刘向军、徐旭常”“选用与流动方向较一致的网格体系，应用a l e 算法对四角切向燃烧煤粉锅炉炉膛内冷态流场进行了模拟，得到稳定的数值计算结果。文 4 2 将a l e 描述方法与有限体积法相结合对流体的控制方程进行空间离散，采用预测一多步校正的方法求解流体和结构的运动方程。王建军等人“”采用g a l e r k i n 格式和n e w m a r k 法对a l e 描述的n s 方程进行空间域和时间域离散，得到相应的n e w m a r k 一压力修正求解公式，通过算例的分析得到较好的结果。岳宝曾等人“”对a l e 描述法在非线性流固耦合和有限元中的运用作了进一步的研究。他们在计算中采用时间分步格式，在空间域上利用g a l e r k i n 加权余量法对系统方程进行数值离散，并分别对圆形贮箱中三维液体非线性晃动和结构与t l d ( t u n e dl i q u i dd a m p e r ) 装置之间的耦合问题进行了数值模拟。孙江龙、叶恒奎“”直接从n s 方程和连续性方程出发，采用四边形单元对求解区域进行网格划分，并借助g a l e r k i n 加权余量法导出了相应的a l e 有限元控制方程。他们在计算过程中速度项和压力项均采用线性插值方法进行插值。143a l e 描述的网格运动a l e 描述中网格的自由移动给流固耦合问题的处理带来了很大的方便，但同时也给方程的求解增加了难度，增加了计算的工作量“”，且存在迁移颈的影响。因此选取合理的、有效的网格运动算法是a l e 描述中的一个重要的问题。a l e 的网格运动是通过控制网格的运动速度或者控制质点在参考坐标系下的运动速度来实现的，运动学边界条件对网格的运动形式有着重要的影响。在边界运动已知的情况下，可以根据所求解的问题的特殊性预先指定网格的运动形式( 网格的速度或位移) 。如l i u 等人。”在求解弹塑性压力波的传播问题时采用了匀速网格运动以消除数值解的振荡。在h u g h e s 等人“7 1 的研究中，对网格运动，首次提出了拉格朗日一欧拉矩阵法，其基本思想是在整个求解区域把网格的速度表示成：= ( ，- a ) v( 1 1 )7输流管道在粘性流动1 _ _ 耦合振动问题的数值方法研究式中v 为物质速度，i 为单位阵，彳为对角阵称为拉格朗日一欧拉矩阵a 式( 卜1 )表明：当a ：1 时得到拉格朗日描述，当| = d 时得到欧拉描述方法，而当d 取其它的任意值时，计算网格可以任意的运动。拉格朗日一欧拉矩阵法需要在整个求解区域内求解网格的速度。h u e r t a 等人”在求解自由表面大幅晃动问题时，根据问题的特殊性，引入了混合方法计算网格速度。混合法的特点是在液体自由表面的法向上使用物质描述，但在切向方向上要根据具体的情况给定网格的运动速度。在求得自由表面的法向网格速度后，区域内部的网格运动则可以用其它的方法加以确定，如直接应用边界上的网格速度进行线性插值。”等。网格的扭曲一般有两种形式：压缩扭曲和剪切扭曲的原理。b a y i dj ，b e n s o n “”采用了一秭基于以上两种网格扭曲形式的网格运动算法。当某节点周围最小单元的面积和最大单元的面积之比小于某一给定值或者当该节点周围单元的顶角小于某一给定值时，则认为该节点周围的单元发生畸形变形，该节点需要移动以调整网格的形状来消除单元的畸形。节点的移动方法可以采用等势松弛法等方法。在g i u l i a n is 4 9 的研究中，提出了另外一种度量单元压缩扭曲和剪切扭曲的指标，并将这两种指标的平方作为评价网格形状好坏的目标函数，利用最优化的方法，馒这个目标函数取最小值以便网格形状最优。另外s c h r e u r s 等人6 ”也提出相应的网格运动算法。国内的王跃先等人“基于三角形畸变能，设计网格运动算法，该算法将求解节点位移( 速度) 的问题转化为一个无约束最优化问题。目前，a l e 描述与有限元相结合的方法越来越多地被从事流固耦合力学研究的工作者所接受，a l e 描述也大量地被应用于流固耦合问题的研究中“”。虽然a l e 描述提供了条有效的途径将拉氏系与欧氏系联系起来，但是要真正地将两者各自有效的方法褶结合求解非线性流固耦合问题，其任务仍然十分艰巨。当前a l e 描述方法的研究动态是寻求简单而有效的a l e 计算格式( 这包括网格运动算法、方程的求解策略、迁移项影响的处理以及交界面上相容条件的处理等) 来进一步完善a l e 有限元方法在流固耦合动力学领域中豹应用。其中，网格的运动算法是目前a l e 描述有限元研究的比较活跃的分支，而设计合理、有效的网格运动算法更是a l e 有限元法研究的重要课题。四川大学硕士学位论文1 5 本文主要工作本文深入探讨了a l e 有限元在流固耦合问题中的应用，并详细推导了基于s u p g 方法的a l e 有限元列式、离散和计算。着重研究了a l e 网格运动的插值、重构以及非线性对流项处理。本文主要工作如下：1 ) 、在阅读大量的国内外有关文献的基础上简单地介绍了流固耦合问题的定义、特征及其研究方法；着重地阐述了a l e 描述法的基本概念、研究进展和目前研究动态：详细地研究了a l e 网格运动的插值、算法与重构。2 ) 、基于网格自由运动这一特点，提出在a l e 有限元的计算过程中采用网格单元的局部线性插值方法，对计算网格进行插值和重构，3 ) 、采用以流线迎风p e t r o v g a l e r k i n 格式为基础的s u p g 方法，将权函数和试函数选取为不同的函数，解决由于流体控制方程中非线性对流项而带来的数值振荡。4 ) 、在借鉴前入成果的基础上，编制了基于c + + b u i l d e r 5 0 运行平台的面向对象a l e 有限元程序。5 ) 、对a l e 有限元方法在湍流计算中的应用推广做了初步的探讨。最后，通过对典型数值算例的模拟，证踢了本文s u p g a l e 有限元及网格运动的单元局部线性插值的正确性和有效性。输流管道在粘性流动下耦合振动问题的数值方法研究2 任意拉格朗日一欧拉( 矗l e ) 描述法众所周知，l a g r a n g e 描述法着眼于质点，e u l e r 描述法着眼于空间点，它们在求解非线性流固耦合问题时都存在严重的缺陷。a l e 描述综合了它们的优点屏弃了它们缺点，为两者的联合应用提供了一条有效的途径。由于计算网格可以以任意形式运动，使得a l e 描述能够在不同的求解区域使用不同的描述，并能很好地处理流固交界面的耦合问题以及网格的扭曲问题。本章详细研究了a l e 运动学描述以及在参考坐标系下的运动学、动力学方程的建立，深入探讨了基于流线迎风p e t r o v g a l e r k i n 格式的a l e 有限元列式的形成。重点研究了计算网格的插值与重构、非线性对流项的处理，提出单元局部线性插值方法对网格进行插值与重构。2 ，1a l e 运动学描述m 1在连续介质力学“”中，将连续体在初始时刻b 的构形( 称为初始构形) 记为口。，将连续体在r 时刻的构形( 称为现时构形) 记为d 。，为了确定连续体各物质点的位置。；f 入l a g r a n g e 坐标系或物质坐标系o x 名，丘( 简记为l 坐标系) 。l a g r a n g e 坐标系是固定在物体上并随物体一起运动的，它可以跟踪物质点在运动过程中的整个运动轨迹。质点在初始时刻“的位置由其在l 坐标系中的位置矢量x 完全确定，因i = x 也可被用来识别物体中的不同质点。先了籀述连续体的现时构形，引入空间坐标系或e u l e r 坐标系o x 。x 2 ( 简称e 坐标系) 。e u l e r坐标系设法在空间中的每点上描述出物质运动随时间的变化状况，是固定在空间中不动的坐标系。空间中各点的位置可以由其在e 坐标系中的位置矢量，完全确定，因此可用x 表示空间中的几何点。拉格朗臼描述是以初始构形为参考构形研究物质点x 在空间中的运动规律，因此物质点x 在空间中的位最可以定义为：工= 工x ，力c 2 一i )而欧拉描述则是以现时构形为参考构形来研究空间点x 上物质点的运动规律，即x = x b ，t )( 2 2 )式( 2 一1 ) 描述了同一物质点x 在不同时刻的空间位置，而式( 2 2 ) 则描述了同一空间点x 在不同时刻被物质点所占据的情况。四川大学硕士学位论文与拉格朗日描述和欧拉描述不同，a l e 描述方法另外引入了一个可以独立于初始构形和现时构形运动的参考构形，记作啦。在物体的变形过程中，观察者始终跟随参考构形运动，因而对观察者而言参考构形是固定不动的，而初始构形和现时构形都是相对于参考构形运动。为了确定参考构形中各参考点的位置，引入参考坐标系d 茧f ：6 ，参考构形中各点的位置由其在参考坐标系中的位置矢量喜确定，因此f = 孝( x ，f )( 2 - 3 )描述了物质点z 在参考坐标系中的运动规律，而x = x ( f ，f )( 2 4 )描述了参考点f 在空间中的运动规律。建立在拉格朗日描述方法基础上的有限元法的单元剖分是针对物体进行的，网格点就是物质点，即网格点和物质点是固结在一起的，随物体一起运动。建立在空间描述方法基础上的有限元剖分是对空间进行的，网格点就是空间点，因此网格点是固定在空间中不随物体运动的。与以上两者不同，建立于a l e 描述方法基础上的有限元的单元剖分是对参考构形进行的，网格点就是参考点，网格是独立于物体和空间运动的，可以根据需要自由选择。以上这些映射( 式( 2 1 ) 一式( 2 4 ) ) 都是一一映射。描述从初始构形到现时构形映射关系的雅可比行列式，和描述从参考构形到现时构形映射关系的混合雅可比行列式了都应不等于零，即：扛引剐0 2 引则( 2 - 5 )根据定义，r 时刻某质点x 在空间中的运动速度v 等于质点x 在空间的位置矢量x ( x ，) 对时间的导数，即：垒延! 判a fj 。( 2 6 )参考构形中某点毒在空间中的运动速度，即网格点的运动速度矿等于参考点( 即网格点) 孝在空间中的位置矢量x ( 品f ) 对时间的导数，即旷：掣7 )西i f”“物质点x 在参考坐标系中的位置矢量毒( x ，f ) 对时间的导数输流管道在粘性流动下耦合振动问题的数值方法研究w =( 2 8 )则为物质点在参考坐标系中的速度。在a l e 描述中，计算网格( 即参考构形) 的运动规律可以是任意给定的。通过指定特殊的网格运动规律可以将a l e 描述转化为拉格朗日描述和欧拉描述。旷= 0 ，即计算网格在空间中固定不动，转化为欧拉描述，矿= y ，即计算网格随物体一起运动，转化为拉格朗日描述。矿v 0 ，即计算网格在空间中独立运动，则对应于一般的a l e 描述。在a l e 描述中参考构形是已知的，各个物理量用参考坐标孝描述比较方便，即f = f ( 孝，) ，因此各个物理量的物质导数应通过该物理量的参考导数来计算。因此对工求物质导数，有：嘲+ 掣整理得到：( 2 9 )铲掣掣( 2 - l o )j上式中q = k k 为物质点相对于网格点的运动速度，即迁移速度或称为对流速度。从( 2 - 9 ) 式可以看出，迁移速度c 实质上是将质点在参考坐标系中的运动速度w 变换到现时坐标系中而得到的，也就是说c 是现时坐标系中物质点相对于参考构形的运动速度。因此式( 2 - 9 ) 可以化为：。：鲥挚删( 2 - i i )u i x利用式( 2 - 9 ) 可将物理量f = f ( f ，) 的物质导数表示为：i o f i ：掣i 托i o f ( 2 - 1 2 )a ；y夙 ，敏，式( 2 1 1 ) 是a l e 描述中一个很重要的关系式，利用它可以将控制方程转化为参考坐标系下求解。a l e 描述的重要特征是可以根据需要给定合适的网格运动速度旷，以维持计算网格的合理形状并准确描述物体的移动暴厦，所以网格运动算法在a l e 描述中占有很重的地位。r，：、一似百管一四川大学硕士学位论文2 2a l e 描述的控制方程2 2 1a l e 描述的一般方程物质的本构关系是由物质本身的特性决定的，它与所采用的运动学描述方法无关。因此在a l e 描述中物质的本构方程与其它描述中的本构方程相同，可以直接引用。为了建立在a l e 运动学描述下的连续性方程、平衡方程、能量守恒方程，这里首先利用非线性连续介质力学的方法和高斯定理给出以下物理量：g ( r ) = 口，g 皓，r ) d 和( f ) = 口，h ( x ，灿夏( 2 - 1 3 )的物质导数。“：剖x = j 吼+ 彬毒卜剖。= 双o h 。+ v , 吼a h ，h g2 2 1 1 连续性方程( 质量守恒方程)取一任意的连续体作为研究对象，分别用厂。、以表示描述连续体的物质域g 毛、空间域试和参考域哦的边界，用p 。、p 和多分别表示连续体各构形的密度。根据连续介质力学知识可以得到在不同构形中连续体的质量m 可以表示为：m2j 。p 2j 址p d g 。j 。p 。x( 2 1 6 )式中：p ( o ，) = 7 p ( x ，) ，p 。c 彤，f ) = ，p ( x ，f ) 。质量守恒定律表明，物质体在运动过程中物质质量不生不灭，即物质质量对时间的整体变化率( 即物质导数) 等于零。因此，利用式( 2 1 6 ) 和式( 2 一1 4 ) 可以得到在参考坐标系下的质量守恒方程为：望o tj + 彬害= o在g 域中j 。p ( 2 一1 7 )与修正拉格朗日方法相似，有时在空间域g 中求解控制方程是比较方便输流管道在粘性流动下耦合振动问题的数值方法研究的。利用式( 2 f 曰、式( 2 i5 ) 和式( 2 一i 2 ) 可以把质量守恒方程写成：静c ，毒+ p 瓦a k = 。姐域中( 2 - 1 8 )2 ，2 1 2 平衡万程( 动量守恒方程)动量守恒定律表明。在t 时刻占据着参考域g 的物体的总动量的整体变化率等于施加在物体上的外力之和，即：割，lp 杉d 壤2 峨+ l 西d 壤( 2 圳)式中t 是作用在参考域g 的边界上单位表面上的力，z 式作用于物体中单位质量的体力，由散度定理及式( 2 - 4 ) 和式( 2 - i 7 ) 可以姆上式改写为：乃矧_ a g 一= 警+ 西砚域中( 2 _ 2 0 )上式中一，是定义在参考构形下的第一类克希荷夫应力张量( 即拉格朗日应力张量) ，它与真实应力 邵橱西应力j 张量a 具有如下的关系：，= 7 警( 2 - 2 1 )在式( 2 一) 中东约不对稼性将会给方程的求解带来一定的困难，茵茈为了便于求解，在某些情况下要将该式转化到空间域中。在式( 2 2 0 ) 两边同时乘以了一t 并整理可得：蕾毒+ 兢在q 域中( 2 - 2 2 )2 2 1 3 能量守恒方程能量守恒定律指出，能量( 包括内能和动能两部分) 对时间的整体变化率等于单位时间内外力对物体所作的功与单位时间内外界传入物体的热量之和( 这里不考虑物体从外界获得的其它能量，如电磁能等) 即：甍x l 西e ；= t i ? y d s ：+ l 。j m v t ；一lr j i 蕊：心一2 蟛一篓西p四川大学硕士学位论文式中e 是物体中单位质量上的能量，q 。为热通量。利用式( 2 1 4 ) 、式( 2 - 2 0 ) 和式( 2 - 1 2 ) 以及散度定理可以把上式写为：铡；彤毒2 署一虿o q j 在g 域中式e 中为物体单位质量中的内能。将上式转化到空间坐标系中为：f 2 2 4 )p 安j + 腭警：芋一i o q ,在g域中(2-25)p 瓦i + 腭瓦2 d i 一瓦在g 域中事实上，能量方程只有当存在不同类型的能量转换时才起作用。2 ，2 2 不可压缩粘性流体流动控制方程的a l e 描述设粘性流体流动占据着空间区域q ，其边界为f ，则由a l e 描述法的基本理论可以推得，在a l e 描述下粘性流动的运动方程的分量形式为：尸警+ p 眈一哆尝= 挈+ ，在q 上2 s )不可压缩连续性方程的分量形式为：i d v i = 0r 在9 上( 2 2 7 )其中，k ，分别为流体的流动速度分量和计算网格运动速度分量，m s 一；口为流体的密度，k g m 一；盯。为c o u c h y 应力张量的分量，k g 。m s ；f 为体力向量，ms ，c o u c h y 应力张量的分量盯。可以表示为：6 ”。一p 气+ 2 k ，)( 2 2 8 )其中，p 表示压力，k g m r s ：为流体的动力粘性系数，k g l m 一, s - 1 ；所1为变形速度张量的分量，即，) = 三( 誓+ 鲁】，显然二阶张量d 为对称张量。方程( 2 2 6 ) 、( 2 2 7 ) 和( 2 2 8 ) 所组成的流体力学基本方程组俗称n s 方程。相应的流体边界上有o i r i c h l e t 边界条件y 2 y在厂l 上( 2 2 9 )输流管道在粘性流动下耦合振动问题的数值方法研究口| i 尘：t在j r ，上( 2 3 0 )其中，厂，和r ：是r 的互补子集。初始条件为：在臼上( 2 - 3 1 )2 2 3 固体区域的基本方程由于在固体部分的运动学描述仍然采用l a g r a n g e 描述方法因此固体区域基本方程可以遵照弹性力学的方法来建立，对于小变形问题，在不考虑阻尼影响的情况下，由弹性力学基本理论可以得到固体区域的基本方程组为：平衡方程挈+ 谚：竹( 2 3 2 )x i几何方程铲糖+ 等j。，物理方程= a 岛f + 2 岛( 2 3 4 )边界条件“t l 。2h l ，”止，= ( 2 3 5 )其中盯扩s 。分别为g a u c h y 应力和应变张量，晚为体力分量。忱为流体旖加于结构上的动压力荷载分量。丑= 硐o e、= 焘分别为材料的l a m e 常数。e 为材料弹性模量，u 为泊松比。2 3a l e 控制方程的有限元离散2 3 1a l e 描述的n - s 方程的有限元离散纵观有关文献不难发现，固体力学有限单元列式的建立几乎都是以变分原理为基础。而在流体力学中则由于流体控制方程的变分原理在与其对应的e u l e r 坐标系中难以建立，因而大多数流体有限元都以g a l e r k i n 法等加权余量法为基础1 。不可压n s 方程的有限元分析，根据所选取的基本未知量的不同，可以分为以流函数作为未知函数的流函数公式、以流函数和涡度作为未知函数的公式6四川大学硕士学位论文以及以流场速度和压力作为基本未知量的速度压力公式。这三种有限元列式各有优缺点。由于使用了流函数的概念，前两种列式能够自动满足不可压缩条件，但是在三维的情况下。流函数的构造非常困难。速度一压力列式直接从n s 方程出发，基本概念清晰，并可以用虚功率原理来解释。另外，边界条件的处理也较前二者简单明了。本文引入s u p g 迎风修正，建立在a l e 描述体系下，具有迎风修正的g a l e r k i n 格式速度一压力有限元列式。s u p g 方法是8 0 年代初期由h u g h e s 和b r o o k s 共同提出的其目的是使n s 方程的数值求解格式具有人工耗散能力以解决速度场的数值振荡问题“。其基本思想是，采用试函数的导数进行权函数摄动修正。便可在完全不触动原控制方程的情况下，将所需要的流线耗散项，很自然地融入g a l e r k i n 方法中，从而得到一种精度和稳定性都较好的求解方法，即s u p g 方法。有关s u p g 更加详细的阐述，请参阅文献 5 3 。下面对a l e 描述的n s 方的有限元列式作一简要阐述，其详细的推导过程参阅文献 5 4 。a l e 有限元的网格划分是在参考坐标系下进行的，其计算网格的运动是相对于参考坐标系的，即独立于空问和物质坐标系。将分量形式的运动方程( 2 2 6 ) 转换为矢量形式后，在q 。中积分，可得与其等价的积分形式缈 彤- 幻一) v 矿j 十彤一p 詈一厂一v 盯) p = o( z 嗡)式中，吒是计算网格的运动速度，其具体形式可任意指定；西乞是对应于对流项的权函数而彤则为对应于体力和扩散项的权函数，在经典的g a l e r k i n 格式中，a , o 和巧k 具有相同的形式；而若这两个权函数独立选取，则将导出迎风修正的o a l e r k i n 格式。式( 2 3