1.1 流固耦合动力学的研究进展
1.1.1 计算流体动力学的研究进展
流体动力学这门科学描述液体和气体的运动及其和固体的相互作用。它几乎涉及日常生活的方方面面,并在大部分科学和工程领域中处于核心地位,所以,这是一门宽广的多学科的学科领域[1]。对流体动力学深入理解的需求,促进了应用数学、计算物理和试验技术的不断进展。核心问题是它的基本方程(Navier-Stokes equations,NSE)没有一般的解析解,数值求解也具有挑战性。当前,流体动力学的成果丰硕,令人振奋,部分原因是新近发展的试验诊断技术和并行计算机在流动模拟与分析中得到了广泛的应用,所以,研究者可以获得流体动力学某些复杂、丰富、完整的细节[2]。
鉴于数值模拟具有很好的重复性,模拟时没有物理条件的限制等许多优点,计算流体动力学(computational fluid dynamics,CFD)在流体动力学研究方面凸显出更加重要的作用。CFD是20世纪60年代初伴随计算机技术发展而新兴起来的一门学科。计算流体力学的发展与计算机技术的发展直接相关。这是因为采用数值方法可能模拟物理问题的复杂程度、解决问题的广度和深度以及所能给出数值解的精度都与计算机的速度、内存和外围设备密切相关。计算流体动力学研究主要集中在数学物理模型、计算方法以及网格技术这几个方面。
1.数学物理模型
流体力学中的Navier-Stokes方程在CFD研究基础上分为5个阶段:
(1)求解线性无黏流方程,如小扰动位势流方程。这类方法较成熟,早期应用较多。其局限性较明显,计算精度较差。
(2)求解非线性无黏流方程,如全位势方程,欧拉方程。这类方法较上一类方法有较大的改进,通常是对主流的简化与对壁面附近的边界层的特殊处理相结合,要取得较好结果的关键在于对边界层的合理处理。
(3)求解黏性、时间平均即雷诺时均Navier-Stokes方程。对紊流的Navier-Stokes方程采用时间平均后,出现脉动应力项(雷诺应力),要使控制方程组封闭,必须进行假设,提出紊流模型。目前存在的紊流模型包括了从简单的代数模型(零方程模型)到标准k-ε模型(两方程模型)直至复杂的高阶矩封闭模型,而k-ε模型或其他形式的两方程模型是目前在工程领域中应用最为广泛的紊流模型。
(4)求解黏性、滤波后的Navier-Stokes方程,即大涡模拟(large eddy simulation,LES)。
(5)求非定常Navier-Stokes方程,即直接数值模拟(direct numerical simulation,DNS)。
除此以外,构造符合实际流动工况的湍流模型也是CFD数学模型的重要组成部分。目前,物理模型研究以考虑更多流动机制,如各向异性的非线性(应力应变关系)湍流研究为重点。
2.计算方法
为了实现上述模型方程的数学计算,还必须对这些方程作适当的离散,这就是CFD的计算方法。计算技术主要由两部分组成:方程的离散及离散方程组的求解。解的精度取决于前者,而求解的效率则取决于后者。在CFD中应用比较成熟和普遍的离散方法包含:有限差分法、有限体积法、有限元法。有限差分法是一种古老的也是最为简单的离散方法,它用有限差分来直接近似控制方程中的导数项,由其截断误差来评定精度。有限体积法是有限差分法与流动控制方程所表述的守恒律相结合的产物,因而具有牢固的数学基础,又有明确的物理背景。有限元法则是基于泛函分析理论之上的一种离散方法,在固体力学中有着广泛的应用。大约自20世纪60年代起,就有人开始着手将有限元法应用于求解Navier-Stokes方程的研究。有限元法的最大优点是能适应各种复杂的流动边界,其主要缺点是需直接求解大型线性方程组,需要大的存储。
流动控制方程对于不可压缩流体与可压缩流体的流动所表现的不同的性质导致解法上的差异。对于低速不可压缩流动,连续方程与动量方程便可构成封闭方程组。但速度场必须满足连续方程的约束,而连续方程却与压力没有直接关系从而导致求解的困难。针对这一问题出现了多种解法,如早期的流函数-涡量法,目前被采用的人工可压缩方法、压力校正法与时间分裂法等。而由Patankar和Spalding[3]在1972年提出的压力-速度校正法(semi-implicit method for pressure-linked equations,SIMPLE)得到最广泛的应用并统治了不可压缩流动数值模拟领域。这种方法又发展了几种新版本,SIMPLER(SIMPLE revised)算法、SIMPLEC(SIMPLE consistent)算法和PISO(pressure implicit with splitting of operators)算法。该方法采用的离散格式有迎风格式、混合格式、指数格式、QUICK(quadratic upwind interpolation of convective kinematics)格式、斜分格式等[4]。它可以采用有限体积法,也可以使用有限差分法离散。
对于可压缩流动,连续方程、动量方程与能量方程联立求解,称作耦合解法。耦合解法中可以将时间和空间混在一起离散,如20世纪70年代十分流行的Lax-Wondroff格式、MacCormack预测校正格式和Richtmyer两步格式等;也可将时间与空间分开处理,先单独对空间离散,离散后的方程成为对时间的常微分方程组,再用人们熟悉的常微分方程组的数值积分法对时间积分进行求解。由于后者更便于程序处理,现在被广泛采用。20世纪80年代后期,基于总变差减小(total variation diminishing,TVD)与矢通量分裂(flux vector splitting,FVS)、通量差分分裂(flux difference splitting,FDS)等方法的高精致离散格式(high resolution scheme,HRS)较好地解决了流体力学中跨、超音速的激波精确捕获问题。近年在研究非定常多尺度复杂流动时,出现了紧致格式。紧致格式有着精度高和网格基架少的优点,是目前研究的一种主要格式。除此之外,计算方法研究还涉及带限制器的高阶插值、谱方法、拉格朗日方法、时空守恒元方法等。值得一提的是将基因算法与传统计算流体力学结合在一起,在流体机械最优化设计许多方面显示良好的应用前景。
3.网格技术
在计算流体力学中,按照一定规律分布于流场中的离散点的集合叫网格(grid),分布这些网格结点的过程叫网格生成。网格生成对CFD至关重要,直接关系到CFD计算问题的成败。1971年Murman采用网格的计算首次表明可以取得正确的物理解。1974年Thompson等提出采用求解椭圆型方程方法生成贴体网格,在网格生成技术的发展中起到了开创作用。随后Steger等又提出采用求解双曲型方程方法生成贴体网格。在1976年,Reyhner率先采用非贴体的笛卡尔网格模拟绕进口的跨音速流动,他预见性地深入讨论了发展笛卡尔网格生成方法所要面临的许多问题,包括准确地确定物面边界条件、相临网格单元间光滑过渡和网格自适应等问题。20世纪90年代以来迅速发展的非结构网格和自适应笛卡尔网格等方法,使复杂外形的网格生成技术呈现出了更加繁荣发展的局面[5]。现在网格生成技术已经发展成为CFD的一个重要分支,它也是计算流体力学近20年来一个取得较大进展的领域。也正是网格生成技术的迅速发展,才实现了流场解的高质量,使工业界能够将CFD的研究成果——求解Euler或Navier-Stokes方程方法应用于产品型号设计中。
目前广泛采用的仍是结构型网格。对于较复杂的求解域,构造结构型网格时要根据其拓扑性质分成若干子域,各子域间采用分区对接或分区重叠技术来连接。非结构网格不受求解域的拓扑结构与边界形状限制,构造起来方便得多,而且便于生成自适应网格,能根据流场特征自动调整网格密度,对提高局部区域计算精度十分有利,然而,非结构网格所需内存量和计算量都比结构型网格大得多。
当前,网格技术方面重点突出网格与流动特征的相容性、分块网格以及混合网格技术,众多研究人员对复杂外形的网格生成技术从分块结构网格、非结构网格和笛卡尔网格3个不同的方向展开研究。分块网格主要用于处理复杂几何形式,也用于并行计算。混合网格技术包括结构网格和非结构网格的混合使用。非结构网格生成方法在其生成过程中采用一定的准则进行优化判断,因而能生成高质量的网格,很容易控制网格的大小和节点的密度,它采用随机的数据结构有利于进行网格自适应。近年来人们开始采用自适应的笛卡尔网格来计算复杂几何形状的流场,即在原始的均匀笛卡尔网格基础上根据几何外形特点或流场特点在局部区域内不断进行网格细化,得到精度符合要求、分布又是最理想的非均匀笛卡尔网格,达到准确模拟外形和捕捉激波等目的。
总之,计算流体力学主要向两个方面发展:一方面是研究流动非定常稳定特性、分叉解及湍流流动的机理,为流动控制提供理论依据,开展更为复杂的非定常、多尺度的流动特征和高精度、高分辨率的计算方法以及并行算法的研究;另一方面是将计算流体力学直接用于模拟各种实际流动,解决工业生产中提出来的各种问题。
1.1.2 计算结构动力学的研究进展
由于实际工程结构日益大型化、复杂化,动力学问题在国民经济和科学技术的发展中有着广泛的应用领域。这些结构的破裂、倾覆和垮塌等事故的发生,将给人民的生命财产造成巨大的损失。正确分析和设计这类结构,在理论上和实际上都是非常有意义的课题。
而其中非线性的动力学分析已经成为现代工程结构优化设计、响应预报及控制等问题中普遍存在的问题。非线性动力学分析,由于问题的复杂性,数值分析方法已经成为主要的手段之一。当前学者们所研究的结构,已不再是以前那种简单结构,而是越来越向复杂的方向发展,研究内容包括复杂的几何性质、复杂的材料行为以及复杂的载荷历程等。一般地,它们都具有非线性的性质。非线性问题与线性问题相比具有许多独特的性质,但是也使问题的复杂程度成倍增加[6]。在对结构进行非线性动力响应计算时,必须考虑到上述性质,针对不同的性质,在计算时应有不同的对策。
结构动力学分析中,经典的时间积分方法有中心差分法、Runge-Kuta法、Wilson-θ法、Newmark-β法、加权残值法等。近些年来,在对这些基本算法的原理取得深入的了解之后,科研工作者们针对各种问题提出了多种时间积分方法,其中许多方法通过对经典方法的修正和完善,在解决非线性、冲击等问题方面提高了计算性能,获得了更高的计算效率。如修正的梯形方法[7]、其他形式的中心差分方法等[8]。另外一些方法,开辟了新的求解思路,如时间域有限元方法[9]、精细时程积分方法[10]、Taylor级数方法[11]以及微分积分法[12,13]等。
对于非线性结构动力学系统,数值仿真算法常用的有以下三大类:
(1)计算数学领域研究的针对一阶常微分方程发展的方法。该类方法是将原运动微分方程转化为状态方程。求系统的响应就是解一阶矩阵微分方程的初值问题,其数值解法的实质就是将微分方程化成差分方程求解。常微分方程初值问题的数值解法有多种,最常用的如欧拉法、中点法和显式、隐式和半隐式Runge-kuta类方法,以及最近发展起来的精细时程积分方法。总体来说,一阶常微分方程组初值问题的数值解法可以处理时变、非线性系统。但对多自由度系统来说,要选取一个保证算法稳定的步长比较困难。
(2)力学界主要研究的逐步积分类方法。该类方法的基本原理是:设定位移、速度和加速度的近似表达式,如有限差分表达式、插值函数等;将这些表达式代入运动微分方程,导出积分格式;最后利用积分格式进行逐步求解,即由前一个或几个时间离散点上的位移、速度、加速度推算下一时间点上的位移、速度和加速度。这半个多世纪以来,尽管涌现出了大量的相关研究成果,目前在像ANSYS、ADINA等这样著名的大型商用有限元分析软件均采用逐步积分类方法,从这个角度上看,逐步积分类方法在非线性结构动力分析中是很有竞争力的,因此对这类方法的设计和算法分析理论的深入研究是非常有意义的。
(3)时间有限元法。将计算的时间区间[t0,t1]等分成n个子区间[tj-1,tj](j=1,2,3,…,n),每一个子区间称为时间有限元。将系统在[tj-1,tj]内任意时刻的坐标、速度和加速度通过插值函数用该时刻的坐标、速度和加速度表示,将该时刻的激振力、质量、阻尼和刚度矩阵也通过插值函数分别表示成各自在时元起、终点所取值的函数。将上述表达式代入运动微分方程,通过上面提到的一些方法可以推导出时间有限元法的递推计算格式。由这种递推格式就可以从前一离散时刻的位移、速度值计算出后一离散时刻的位移、速度值。
所有这些方法从分类上说有显式和隐式积分两大类。显式积分每一步效率较高,但时间步必须取得相当小才能保证其稳定性。隐式积分法则可以通过恰当的参数选择,保证积分的数值稳定性,因此积分步长大一些也可适用。熟知的Newmark-β法及Wilson-θ法都是隐式积分格式,然而这些积分法也有其弱点,每一步至少做一次系数矩阵的三角分解,这对规模较大的问题来说是很费时的。另外由于积分步长太大时,一些需要的高频振动分量不能正确地反映出来。如果用于保守系统,使用某些算法经过若干步的计算后,系统的能量不能保持守恒,这种情况被称之为“人工阻尼”或“算法阻尼”,这是数值积分方法本身带来的,是不可避免的。这些方法都可能因为逐步计算的特点而逐渐累积计算误差。因此,对不同的问题需要使用不同特性的算法,这就要求使用者对不同的算法的本质都要有好的了解,这样对不同的具体的问题就能正确地选择相应的数值算法。
求解时间离散后的非线性有限元方程组,通常又有3类方法,即增量法、迭代法和混合法。而基于牛顿-拉夫逊法的迭代方法在求解非线性有限元时使用较广。
1.1.3 计算流固耦合力学的研究进展
流固耦合问题是目前多物理场、多学科交叉领域研究的热点问题之一。自20世纪70年代流固耦合问题提出以来,研究势头不断上升。在80年代中,该领域的研究又得到了迅速发展。进入90年代后,随着计算技术的进步,研究范围不断扩大。从总体上看,流固耦合问题按耦合机理可分为两大类:第一类问题的特征是两相域部分或全部重叠在一起难以明显地分开,其耦合效应体现在其控制方程中;第二大类问题的特征是耦合作用仅仅发生在两相交界面上,耦合是通过两相交界面上的平衡与相容条件反映的。第二大类的问题又可分为3个子类:第一类是两相间有较大的相对运动的问题,包括水动弹性力学、气动弹性力学、输液管动力学等问题;第二类是流体有限位移的短期问题,如爆炸载荷作用下流固耦合系统的瞬态响应问题(冲击问题);第三类是流体有限位移的长期问题(响应问题)。本书主要研究第二大类问题中的第一子类问题。从数值方法上看,目前研究流固耦合问题主要采用两类方法:第一类是基于移动网格技术的耦合方法,如最具代表性的是任意拉格朗日-欧拉法(arbitrary Lagrangian Eulerian,ALE);第二类主要是基于固定网格技术的耦合方法,如浸入边界法[14-16](immersed boundary method,IBM)、浸入边界有限元法[17-19](immersed finite element method,IFEM)以及基于IBM的无网格法[20,21]等。本书主要研究采用第一类基于移动网格技术的耦合方法,以便准确捕捉流体和固体运动的边界。
历史上,人们对于流固耦合现象的早期认识源于航空工程中的气动弹性问题[22]。一般认为1933年Westergarrd[23]首次对作用在刚性垂直坝体上的动水压力的研究,标志着流固耦合理论研究的开始。之后,Kotsudo[24]和Chopra[25]改进了Westergarrd的数学模型,研究了在地震过程中,作用在坝体上的动水压力。在此期间,Blevin[26]对输液管道的振动问题做了大量的研究,Haskind[27]对船水耦合动力学问题做了开创性研究。同时我国学者也做了大量的研究,郑哲敏[28]对不可压缩的理想水体,应用速度势函数研究了半无限水域悬臂板的振动问题,居荣初[29]对水中圆柱体结构的流固耦合问题做了大量研究。直到20世纪60年代,有限元和边界元等数值方法出现以后,流固耦合问题才又为人们所关注。自70年代开始,由于数值计算方法和计算机技术突飞猛进的发展,给流固耦合理论研究带来了新的活力,众多学者采用有限元、边界元及其混合法等数值方法对流固耦合问题做了大量工作,使流固耦合数值方法研究得到较大的发展。如Zienkiewicz[30]等人在用有限元分析储液罐振动问题时,认为可以用一个附加质量阵来修改固体质量矩阵,液固相互作用通过流体的附加质量来反映,这种方法形式简单,计算量小,在工程中已被广泛接受,但附加质量的大小不好确定。刘正兴[31]采用位移-压力格式,对于流体介质中结构的动力特性及响应进行了计算分析,并提出了一般结构物附加水质量计算范围的确定方法,这对实际工程计算有一定的借鉴作用。Chwang[32]分析板在水中运动时,首次考虑了非线性表面波对动水压力的影响。Bath等[33]人采用位移-位移格式,建立了有限元方程,对浮动结构进行了流固耦合非线性分析。国内学者在20世纪70年代后,也进行了很多有价值的工作。如沈惠明、赵德有等人[34]采用流体边界元和结构有限元法求解了双面浸水板的流固耦合问题。傅作新[35]在挡水坝动水压力方面做了较多的研究工作,他采用流体边界元和结构有限元法分别探讨了刚性坝和弹性坝的动水压力,并考虑了地基与库水的相互作用及库底泥沙吸收的影响。张升明、吴士冲[36]以流体附加压力作为基本未知数,利用极值原理建立离散形式的流体运动方程,分析了弹性结构与流体相互作用的耦合振动问题,建立了流固耦合振动的一般方程式,并编制了相应的计算程序,对储液器进行了水弹性分析。黄玉盈[37]采用解析法分析了中厚浮板与水的耦联振动问题。杜修力、陈厚群[38]等人分析了拱坝系统三维非线性地震波动问题,但在水体的处理上采用了附加质量法。我国的李遇春[39]在研究渡槽内水体的非线性晃动问题时,也考虑了水体的非线性表面波的影响。张立翔[40]结合管道系统,在流体-结构互动的基础理论、计算方法、试验以及工程应用等方面做了大量的研究工作。刘晓旭[41]针对油气田开发,发展了以流体孔隙压力,温度和孔隙介质位移为基本变量的三场耦合的数学模型研究及有限元解法。杨向龙[42]研究了人体内分叉流动和人工心脏中的流固耦合问题,使用CFD求解器Fluent来求解黏性流体的流动,用自编制三维非线性有限元程序来求解薄壁结构的非线性变形,通过交替计算的方法来实现流体和结构的耦合。
从流固耦合模型发展历程来看,对气动弹性(或水动弹性)研究,发展了几种有效的计算模型。从流体动力学方面看,从最简单的势流模型发展到复杂的三维雷诺平均的湍流模型(Reynolds averaged Navier-Stokes,RANS),从结构动力学方面看,从简单的线性有限元梁模型发展到复杂三维实体有限元模型,同时,用以联系流体和固体求解器的有效的界面交换技术也得到不断地发展。流固耦合研究的发展历程可简单总结见图1.1。最初,Cunningham等[43]将势流流体控制方程与结构自振模态耦合起来,发展了一种计算跨音速气动弹性分析的耦合计算格式(以下简称为CAP-TSD模型)。该方法将升力面简化为薄板平面,简化了网格生成任务,当界面速度边界条件应用于薄板中性平面时,不需要额外的边界移动格式。但对于模拟存在强的撞击或黏性不可忽略的流动时,CAP-TSD模型显得无能为力。为了克服上述模型不足,Schuster等[44]将三维可压缩RANS湍流模型与线性化结构模态相耦合,来执行FSI计算(该模型简称ENS3DAE)。Beam-Warming隐式积分格式用于时间积分,代数剪切(Algebraic Shearing)技术用于流场中网格更新。同时,为了更好地预测流固耦合中流体边界层流动,Lee-Rausch和Batina[45,46]将三维可压缩边界层RANS湍流模型与线性正交化结构动力模态计算模型相耦合(该模型简称CFL3DAE)。采用具有二阶精度的时间差分格式,在每一个时间步内采用虚拟-时间(Pseudo-time)子迭代步以加快收敛速度。基于弹簧法则(spring analogy)的移动网格技术用于边界层网格更新。该模型能很好地预测机翼颤振下的边界层流动。Liu[47,48]采用并行、多块以及多重网格技术求解全三维的Navier-Stokes方程,并耦合结构动力学方程,采用双重时间步积分格式(dual time stepping scheme),保证流体解与固体解同步推进,采用TFI(transfinite interpolation)与弹簧法则相结合的动网格模式。此方法主要用分析飞机机翼的气动弹性问题。Farhat等[49,50]将固体和流体部分都采用三维非结构网格,将ALE方法用于变形的网格系统,分析F-16战斗机在飞行过程中的气动弹性问题。最近,Kamakoti等[51,52]发展了一种求解RANS方程的多块、结构化的CFD求解器,采用全隐或半隐的同步时间积分格式将固体动力学求解器与CFD求解器耦合在一起。采用基于有限单元的线性梁单元,对AGARD445.6机翼进行颤振分析。从以上分析可见,目前在流固耦合的CFD求解方法上,多采用基于RANS方程对湍流流场进行分析,而更高级别的流体数值模拟方法,如LES和DNS报道较少。
图1.1 流体和固体求解器同界面插值技术发展历程
总的来说,流固耦合研究成为多学科交叉领域的研究热点之一。纵观研究历程,早期的研究普遍采用解析法,对流体和结构运动规律进行大量简化,大部分研究都是将流体视为无旋、无黏的介质,而现在一些学者开始以黏性复杂的流体运动为研究对象;其次,早期的研究基本上都是针对平衡位置附近的小扰动、小变形进行简化分析,而现在主要针对大的相对运动、大变形进行研究;再者,以前的许多工作主要是线性频域的研究,而现在开始转向时间域的非线性问题。总之,简化越来越少,研究范围越来越广,研究的几何模型和物理模型也更趋复杂。
1.1.4 计算流固耦合力学的研究特色
下面将对流固耦合过程中遇到的特色问题如任意拉格朗日-欧拉描述、动态网格几何守恒定律和网格更新算法以及耦合界面不匹配网格间运动和载荷的传递等进行概述。
1.任意拉格朗日-欧拉描述
在流固耦合系统中,首先遇到的问题是坐标系描述问题。结构分析中习惯用拉格朗日(Lagrange)坐标系,着眼于物质点。流体力学中则大多使用欧拉(Euler)坐标系,着眼于空间点。这种描述方法的差异,对小变形问题,区别不大,但对大变形非线性问题就需要引起注意。为了统一耦合界面上的运动描述,Hirt等[53]、Hughes等[54]、Donea[55]、Ramaswamy and Kawahara[56,57]、Huerta and Liu[58]、Soulaimani等[59]提出一种任意拉格朗日-欧拉描述,充分利用了拉格朗日描述能精确描述边界运动和欧拉描述不会导致流场网格畸变的优点,能有效解决边界运动的流动问题。Tezduyar等[60,61]、Masud and Hughes[62]和Hansbo[63]将ALE描述引入到有限元方法中,来求解有限元时-空动边界流动问题。最近,Soulaimani等[64]和Braess等[65]采用有限单元法离散ALE控制方程,来求解自由表面流动问题,Belytschko等[66]和Dettmer等[67,68]将ALE描述下的流体控制方程用来求解流固耦合问题。
2.动态网格的更新与几何守恒定律
在流固耦合问题中,耦合界面上的运动导致另外一个问题,即流体动态网格更新。早期的耦合问题研究中多假设结构做微幅振动,网格只有速度没有位移,因此不存在单元畸变和网格移动问题。现在一般多采用动态网格方法,即网格既有速度描述又有位移描述,以提高对大变形问题的处理能力。为保持整个流体域的网格不因边界运动而导致的单元畸变,需要及时更新网格以追踪耦合界面并保持良好的单元形态。最初Batina[69]和Robinson等[70]分别采用弹簧模拟法来移动非结构网格和结构网格,假设流体网格的节点间是用虚拟的弹簧相连,弹簧的刚度与节点间的距离相关联。此种方法能处理大的网格变形,但计算时采用迭代方式并耦合椭圆网格生成器,网格数量太大时,计算比较费时。Eriksson[71]提出三步无限插值方法生成单块区域网格,由于该方法不能用于多块网格区域,Hartwich和Agrawal[72]组合弹簧模拟和无限插值两种方法,用于更新多块网格区域,区域边界线网格由弹簧模拟法完成,区域表面及其内部体积网格由无限插值方法完成。Farhat[73]采用拟结构方法更新网格,即把流体作为一个准结构来看待,给定虚拟的材料常数,通过求解每一时间步的该拟结构的静态平衡方程来获得更新后的流体网格。另一类较为经典的重新生成网格的方式是求解控制移动网格的偏微分方程[74-76]。通过建立移动网格区域内所用网格节点的控制方程并求解,并辅以控制参数来光滑求解后的网格。
对于基于动态网格的流动问题,除了需满足三大物理定律外,还需要满足流体网格运动与流体积分求解方法之间的一致性和相容性,即要满足几何守恒定律[77](geometric conservation law,GCL)。违反GCL条件可能导致额外的数值振荡从而导致计算失败。具有一阶时间精度并满足几何守恒定律的有限体积和有限单元数值格式详见[78,79],具有二阶时间精度并满足几何守恒定律的有限体积数值格式详见[80]。Farhat等[81]基于ALE法则总结了6种不同的时间积分格式(其中一些时间积分格式并不满足GCL条件)对求解精度的影响。以上这些讨论主要是基于密度的流体求解器,Kamakoti和Shyy[82]提出基于压力的流体求解方式下GCL条件的执行,并分别讨论了具有一阶时间精度隐式积分格式、一阶精度时间平均积分格式、二阶时间精度隐式积分格式和二阶精度时间平均积分格式。最近,文献[83]和文献[84]提出了满足GCL条件下基于压力的流体求解器,并被广泛用于复杂气弹性研究。
3.耦合界面插值算法
对于同步强耦合或交错弱耦合计算,得到3个主要的单独计算模型,即CFD模型、计算结构动力学(computational structure dynamics,CSD)模型和动网格模型,之后需要发展一种有效的界面耦合程序,以联系这些单独的计算模型。耦合界面不匹配网格间信息的准确交换,直接影响耦合计算精度,是获得真实物理解的关键环节之一。
目前,几种有效的界面插值程序已得到广泛的应用,例如无限平面插值(infinite plate splines,IPS)、有限平面插值(finite plate splines,FPS)、平面薄板插值(thin plate splines,TPS)、多重二次曲面插值(multiquadric biharmonics,MQ)、非一致B样条曲线\面(non-uniform B-splines,NUBS)插值、逆向等参插值(inverse isoparametric mapping,IIM)等,其优缺点见表1.1。但在实际耦合计算过程中,选择何种界面信息交换方法,往往是界面信息交换计算所需的时间和计算所需精度的折中。
表1.1 几种常用的流固耦合界面插值技术汇总
1.1.5 流固耦合力学模型及算法
流固耦合系统中,通过彼此相连的耦合界面,流体力作用在结构表面影响结构的运动,而固体的运动造成流体计算域及边界条件的改变从而影响流场。在耦合界面上的这种非定常流体力和边界运动状态事先都无从得知,只有通过耦合系统的求解才能获得,这正是耦合系统的非线性特征。不过,由于耦合问题特别是动态耦合问题的复杂性,在整个研究历程中,总是根据当时的计算能力作出必要的假设,以期在揭示物理本质与可忍受的计算量之间找到最佳的平衡。
流固耦合问题的数值研究方法,其模型经历了由简单到复杂(图1.1)、从弱(松)耦合发展到同步迭代强耦合,最后到整体积分紧耦合这样一个过程。从众多的文献资料来看,根据这一发展过程大致可把各种方法分为以下4类。
1.古典分析方法
古典分析方法严格来说不能称之为一种流固耦合方法。因为该方法通常是假定结构以给定的频率和幅值进行振动,通过计算非定常力来判断系统的稳定性。它把一个原来内部耦合的非线性问题分解成了两个独立的解耦问题,因此不能用来预测非线性极限环振荡现象。而且由于流体-结构间不存在封闭的反馈,因而不能准确反映两种连续介质间的能量传递。
常用的古典分析方法有能量法、特征值法等。这些方法通常以某一刚性特征型面进行分析,通过系统总阻尼来判断其动态稳定性。变形激盘法是其中的一个典型代表。该方法属于一种介于经验法和计算流体力学方法之间的半经验方法。该方法从流体力学的基本规律出发,有其严格而概念清晰的一面,可以用一些试验数据来处理复杂的流场现象。
2.交错积分弱耦合法
交错积分弱耦合法和下面的整体积分(或完全积分)耦合法都属于双向的耦合方法,因此可以考虑流体和结构的各种非线性因素,也可以预测到极限环振荡现象。两者的差别主要体现在时间积分推进是否同步上。
交错积分耦合法中,流体和结构用各自的求解器在时域积分,交错时间推进,其典型过程是:①在流体力作用下对结构响应进行积分,推进至下一时刻;②把结构的边界位移和运动传递给流体系统;③更新流体动态网格;④流体积分,计算新的流体压力和应力场;⑤把流体压力和应力转换成结构载荷,传递给结构。文献[85]称这一过程为顺序交错耦合(conventional serial staggered,CSS)。这种方法能充分利用现有的计算流体动力学和计算结构动力学的方法和程序,只需做少量的修改,从而保持程序的模块化。它的一个主要缺点是,由于积分是交错进行的,流体、结构的时间推进积分总是存在时间上滞后(time-lap),耦合界面上的能量不能保持守恒,无法满足动平衡。即使流体或结构解具有二阶时间精度,但耦合解也只具有一阶时间精度,因此该算法是一种弱耦合过程。由于流体和结构的特征时间尺度不同,为节省CPU时间及I/O文件传递,引入流体积分的子循环方法可使积分过程中采用不同的时间步长。文献[86]~文献[89]进一步详细描述了交错积分耦合算法、固体和流体之间能量的传递以及载荷与运动在非匹配网格间的转换方式及应用等问题。
3.同步迭代强耦合方法
为提高算法的精度和稳定性,在CSS耦合的每一个循环过程中,采用带预估和预估迭代的交错方法(predictor multi-corrector algorithm,PMA),在每一个时间步上进行多步迭代,直到满足计算精度要求后,流体和结构解再整体向前推进。此种强耦合方法克服了CSS耦合中流体、结构的时间推进积分总是存在时间上滞后的缺点,又没有整体积分紧耦合方法在数学模型等方面的困难,是目前流固耦合研究领域中最活跃的研究方向之一,并可以用于较为复杂结构和流动之间相互耦合计算。最初Brooks[90]将PMA方法用于求解有限元离散的不可压缩Navier-Stokes方程,后来Hughes[91]将PMA方法用于求解有限元离散的结构动力学方程,最近文献[92]~文献[94]将PMA方法用于求解多物理场的耦合系统。
4.整体积分紧耦合方法
整体积分紧耦合方法是把流体和结构看作通过耦合界面连接的单一连续介质,用单一的算子来描述控制方程。由于时间积分完全同步,且不存在任何时间滞后和能量不守恒现象,因此是一种具有相当吸引力的耦合方法。然而,整体积分紧耦合方法一方面要在欧拉坐标下处理流体运动方程,又要在拉格朗日坐标下处理结构运动方程,同时,在解单一的耦合控制方程时,导致结构部分的刚度矩阵远远大于流体部分的刚度矩阵,因此目前采用整体积分紧耦合方法来求解复杂结构与流场的耦合分析还不现实。最初Bathe等[95]采用基于有限单元的紧耦合算法,求解具有自由表面的可压和不可压缩流体与刚性固体之间的流固耦合计算。后来Rugonyi等[96]进一步发展整体积分的紧耦合有限单元法,求解考虑黏性流体与弹性固体相互作用的整体耦合系统。最近Heil[97]和Hubner等[98]采用紧耦合算法,执行黏性流体与大变形结构之间的流固耦合计算。但总的来说,在求解整体耦合方程时还存在许多求解上的困难,同时由于紧耦合计算相当耗时,目前主要用于简单的二维问题,关于这方面的研究国内外报道还比较少。