[doc格式] 用于ALE有限元模拟的网格更新方法

[doc格式] 用于ALE有限元模拟的网格更新方法 用于ALE有限元模拟的网格更新方法 第4()卷第2期 2008年3月 力学 ChineseJournalofTheoreticalandAppliedMechanics 用于ALE有限元模拟的网格更新方法 周宏李俊峰王天舒. (清华大学航天航空学院,北京100084) 摘要任意拉格朗日欧拉法(ALE)可以通过定义参考网格的运动,实现自由液面跟踪,完成液体晃动的数值 计算.综合用于更新网格节点的3种基本计算方法,将多方向更新网格速度的技术应用于任意拉格朗日欧拉 网格节点的速度计算.给出了水平圆柱形贮箱和椭圆形贮箱内液体晃动算例,实现了多方向更新网格运动与晃 动流场计算的耦合,使ALE方法能胜任复杂几何边界下的自由液面流动的数值模拟. 关键词ALE方法,有限元方法,网格更新,液体晃动,水平圆柱形贮箱 中图分类号:O130.25文献标识码:A文章编号:0459—1879(2008)02—0267—06 引言 晃动(sloshing)是指有限体积容器内液体自由液 面的运动.采用任意拉格朗日欧拉方法(简称ALE 方法)描述带自由液面的晃动问题,能将自由面上 流体通量为零的运动学边界条件,转化为自由面上 ALE网格节点速度的法向投影等于质点速度的法向 投影.因此在用ALE有限元方法仿真晃动现象时, 计算区域内将存在动态的ALE网格与流体运动之间 的域耦合.首先,根据流场的质点速度计算自由面 上网格节点速度,并估算下一时刻自由界面的空间 位置和几何形状.然后从已获得的自由界面的边界 条件中计算ALE网格内部节点的速度,尤其是靠近 自由液面区域的网格节点速度.最终,在更新后的 计算网格上完成流场控制方程的空间离散和计算. 1981年,Hughes等【lJ提出了L—E(Lagrange— Euler)矩阵方法,通过对单元节点定义不同的L—E 系数与相应的流场质点速度,共同确定网格内部节 点的运动.1982年,Donea等L2.提出了平均法, 即用上一时刻相邻网格点的速度平均值确定内部网 格节点速度.1988年,Huerta等L3.提出了混合 法,即只在自由液面或运动边界上进行L—E网格更 新,内部网格节点速度用势流方程或代数方法计算 得到.用于ALE有限元计算的网格更新方法都是以 此3种方法为基础,其中Huerta提出的混合法应用 最广.另外,在混合法中,基于代数思想的平滑技术 多用于几何区域内的流场分析;而利用势流方 程等微分方法的平滑技术在非直壁面区域的流场计 2007-01—29收到第1稿,2007-12-21收到修改稿 1)国家自然科学基金资助项目(10572022,10302013). 2)E-mail:tswang@tsinghua.edu.cn 算中具有更强的适用性-4J. 微分方法根据更新对象的不同,可分为更新网 格的位移和速度两种:第1种,将网格看作假定机 构(Pseudo—structure),比如将计算网格理想为弹簧 振子系统L5J,计算网格节点的位移;或者假定为弹性 体,利用弹性体的本构关系和几何关系直接在力平 衡微分方程中计算节点位移.计算网格节点的形变 量L6’7J,获得网格节点位移.第2种方法是利用势流 原理直接获得网格节点的速度,只要网格节点速度 的梯度在计算域内为常数,那么网格将会保持原来 的尺寸比,即网格长宽比不发生大的变化,更新后 的单元接近原有单元作刚体运动后获得的新位置而 没有发生大变形. 本文将重点研究ALE网格内部节点速度更新 的数值方法.这一步骤关系到能否在自由面附近获 得性态良好的计算单元,能否降低大幅度液体晃动 计算的网格重构频率,同时网格节点速度的数值还 将影响到控制方程中对流项的量级.更重要的是作 为用于ALE方法中的网格更新方法,改进自由液面 网格节点速度的求解,能使其在更多样的几何 边界和运动边界问题中得到广泛应用. 1ALE方法描述下自由液面网格速度及边界条件 任意拉格朗日欧拉法(ALE)描述的有限元计算 方案,在1982年首次被Donea用于分析液体一结构 耦合的瞬态问题[2】.该文详细阐述了ALE方法的基 本理论,写出了在ALE参考坐标下的液体运动的质 268力学2008年第4O卷 量,动量和能量守恒式.根据文章中提到的ALE坐 标系下的质点导数理论,假设自由液面的曲面方程 为F(x,t)=0,那么对方程两边求导 dF : I+(—tI,).F:o(1)一=一l+IT,一TI,l?V=lIIlldtat1w’, 显然,ALE坐标下OF w = .,根据上式得到自由 面上流场应满足的运动学边界条件是 W)?F=0 (— 即网格节点速度在自由面法向上投影和自由液面流 体质点速度的法向投影相等.通过迭代能获得时刻 t的自由液面上节点的运动速度. 同样,根据流场在固壁面上的边界条件为可滑 移边界条件:(1)壁面法向上液体的网格速度投影 和刚体或固体的速度投影相等;(2)在切向上不对液 体的网格运动速度作限制.在非运动边界上或者远 离运动边界的区域上,还可将网格速度定义为零,这 就是ALE网格与Euler网格的分界线.用数学式可 以将ALE网格区域的边界条件表示成以下形式:在 液体的自由面上 tI,l,=l,(3) 其中,西l,的值通过公式(1)获得,或者直接取液 体自由液面的速度.在运动界面上,有 l=l(4)anIanI, 是液体网格节点与运动界面重合处的运动速度. 在距离运动界面一定距离后可定义网格速度为 tI,lb=0(5) 除了以上这些边界速度满足条件以外,在一些 特殊的区域,如自由液面与运动界面相交的接触线 区域,还需要对网格速度进行特殊的数值处理,以 保证ALE网格不在这些区域出现网格缠结等畸变. 2应用于ALE有限元的移动网格方法 如何在ALE区域内,将边界上的网格速度均匀 到区域内部,以平滑边界位置变化带来的网格变形, 是本节主要关注的问题.现有两种更新网格内部节 点的基本方法:一种是以更新网格位置为目标的节 点位移更新方法,另一种是以获得网格节点运动速 度为目标的节点速度更新方法. 节点位移更新方法是广泛用于有限元计算的网 格优化技术,它常与网格自动生成技术一起,成为 有限元计算前处理的核心步骤.这类方法在ALE有 限元方法和时空有限元方法出现后得到推广,但在 ALE有限元方法中使用节点位移更新方法,是期望 通过计算上一时刻网格节点的位移,获得ALE网格 的节点运动速度.另外,位移更新的网格优化技术, 会给数值仿真带来额外的计算量,一般都是在代数 方法不适用的情况下,才使用节点位移更新方法, 而且还需要尽量使网格位移更新计算代价较小.这 种将ALE计算格式与节点更新技术结合的计算方案 在带运动界面的外流场数值模拟问题中有着大量的 应用. 节点位移更新方法的通常做法是将ALE有限元 网格看作”假想弹性体”,将自由面,运动界面以及 固定面上的运动学边界条件作为”弹性体”变形的边 界条件.在不可压缩液体的晃动问题中,采用线弹 性模型能够保证液体体积守恒,而在液体不满足不 可压缩的条件时,液体体积守恒不是必要的约束, 可以用非线弹性模型.本文将以满足不可压缩条件 的液体网格为例说明该方法在ALE有限元分析中的 应用.网格节点位移的数学模型为 ? +,=0(6) 其中,是柯西应力张量,,是作用在网格上的外 力,界面上的纽曼边界条件 n?=Y(ii) 这里,是结构在运动界面的法向投影.要使这类 方法真正达到均匀网格变形的目的,需要给出假想 第2期周宏等:用于ALE有限元模拟的网格更新方法269 网格的相关物理参数.在非结构网格中,影响材料 参数的因素包括网格尺寸,长宽比和变形率.在变形 剧烈的区域里网格细密,且期望网格变形率为零; 在变形舒缓的区域尽量保持网格良好的长宽比.所 以要求尺寸小的网格单元有良好的刚性,而尺寸大 的网格单元能吸收变形.主要的处理技巧有两类, 控制单元雅克比值的方法Is],和针对单元定义不同 的本构关系矩阵的方法-9J. 有限元方法中,为了完成单元积分会将单元从 总体坐标向局部坐标转换,获得一个几何形状规则 的单元.这个转换过程中出现的转换矩阵叫雅克比 矩阵,它的行列式值能表征单元是不是发生了缠结 和不良的长宽比.因此在求解式(6)的边值问题时, 通过改变单元的雅克比行列式的值就能改变单元在 总体坐标下的形状.在文献【10】中首先采用了这种 方法,即在小尺寸的单元中,通过改变雅克比行列 式的值增加的方法,在Etein的 文章中)(的取值范围可以从0.0到2.0.这样能根据 具体问题,改变刚性参数)(的取值以适应不同的界 面耦合问题,而不需要对模型进行大规模的改变.同 时以单元尺寸为基础的改进方法,能在运动界面附 近更细密的单元,而不增加网格重构的频率, 文章的结论表明在运动界面发生大的位移后,靠近 运动界面的小尺寸单元依旧保持很好的单元性状. 2000年,Chiandussi[引总结了定义网格结构材 料参数的主要方法.文中将单元材料的杨氏模量分 别定义为单元重心到运动界面距离的线性关系,平 方关系和指数关系.根据单元不同的特性(单元应变 量,单元应变能和单元变形能密度),给出每个单元 的材料参数,然后计算出单元的应变场.首先假设 网格各向同性,杨氏模量已知,为常数E,获得一个 应变场e,应力为=Ee.要优化网格变形,就应在 相同的应力状态下,存在常应变,使’=E,满足 这个条件的弹性模量将是我们需要的,即 E:E E(14) 将E带入网格运动的模型获得新的应变场,作为网 格运动位移的依据.文章比较了以结构参数和以几 何参数为选取网格弹性参数的有效性,前者更 高效,尤其在泊松比取0.32时,用变形能密度计算 弹性模量效果最好. 节点速度更新方法能直接求出ALE参考系的运 动速度,在很多情况下,更新速度比更新位移更简 洁.在形状简单的计算区域内,引言中提到的代 数方法就能完成网格速度的向内均匀,即在某一平 行于网格坐标的直线上,以节点距离自由面或距离 运动界面的数值为变量计算内部节点的网格速度, 直到ALE网格的边界.但是这类方法的适用范围 有限,需要更一般的方法来计算内部网格的运动速 度.在1982年Donea就以平滑运动界面临近区域节 点速度的思想,给出了ALE内部网格节点的速度 求解公式 珐=?Jc+0.1L}J?击 (15) 这里?是与节点通过单元各边相邻的节点数目, LIj是节点与相邻节点J问的距离,是t时刻 节点总位移.式(15)的设计思想就是通过平均与 相邻节点上一时间步的节点速度来获得节点当前 的运动速度.并且在考虑运动精度和数值稳定性的 前提下,给出了一个经验不等式来确定网格节点速 度与液体实际流动速度的比值-2].该方法能自然的 形成ALE网格区域与Euler网格区域的界限. 另外一种方法是利用势流方程的微分技术,在 单元网格尺寸一致的前提下,如果速度的梯度在整 个区域内为常值,那么单元将能保持问题初始时刻 的长宽比.在数学上,很容易联想到拉普拉斯算子, 它表示了变量梯度的散度,这样在ALE网格区域内 平滑网格速度的问题就能归结为求解一个矢量的拉 普拉斯方程的边值问题. ? ()=?=0(16) 方程的边界条件为第2节中的式(3)一(5).文 献【11】在上式中加入了一个表征网格特征尺寸的参 数=h-p,P>0,实现了在小尺寸单元中,速度梯 度变化小,在大尺度单元中,速度梯度变化大.从而 270力学2008年第40卷 避免了均匀的网格速度梯度给尺寸大小不同的单元 造成不同程度变形率的影响.写成数学式子为 ? ()=0(17) 其中为对角阵,当i=1时,上式是关于网格速 度的一个拉普拉斯方程,当i?1时,需要通过求 解泊松方程来获得新的网格速度.本文将定义为 网格节点相对运动界面的距离的函数,并给出了 采用这种修正方法与标准的拉普拉斯速度平滑方法 的结果比较,表明此法在非结构网格中能起到均匀 网格速度,减少网格更新频率的目的. 3数值算例 更新规则容器内的运动网格使用的方法相对简 洁,用简单的代数平均思想就能将自由面的运动速 度向内部节点均匀,这种方法与自由面的平滑技术 相结合,就能很好地处理自由液面的晃动问题.但 是在复杂的贮箱内或者在壁面与参考坐标系不平行 的圆柱形或球形贮箱情况下,采用上面提到的网格 更新方法就显得非常必要. 本文将以一个水平放置的圆柱形边界的计算网 格为例,采用上节提到的节点速度更新方法将式(17) 不同时刻内的网格形状和节点速度表示出来.在图 1中是当晃动振幅与圆柱直径的比值大于15%后, 网格的形状和网格节点的速度场. 在图1中,可以看到自由液面上网格节点和内 部的网格节点速度都具有X,两个方向的自由度. 大部分同类文献只能完成矩形或直立圆柱形液体晃 动算例的数值模拟,就是因为网格速度只能完成单 方向更新.算例表明多方向更新ALE网格节点速度 的方法能用于此类边界形状的液体晃动计算. 本文第2个算例为长短轴比例为2的水平放置 椭圆腔内的数值晃动模拟.椭圆长轴a=800mm, 短轴b=400mm,充液深度h=800mm,在液深 为长轴半径的情况下,文献[12】提供的直立放置椭 圆形贮箱的一阶近似晃动频率的实验数值为u= 5.1972rad/s.算例外加激励为 = Agsin(wt+)(18) 其中,是激励的无量纲幅值,它是一个常数,该算 例中A=0.05;g=9.8m/s是重力加速度;=7r 是外加激励的相位;t是时间变量.此算例中,采用 了可滑移壁面边界条件,边界上的网格速度除了满 足公式(4)的要求外,对边界上网格节点的切向速 O O — O — O — O O O — O — O — O O O — O — O — O O O — O — O — O —— 0.6——0.20.20.6——0.6——0.20.20.6 T=137r/7T=27r 图1圆形贮箱内网格节点在T/2个周期内的速度场 Fig.1Nodalvelocityfieldsduringthehalfperiodictimein cylindricalcontainer 筹I)(等I)(19) 定义的分段线性函数,当=詈,:2a,有 c= 1 3 - 兰)<h,ic2. 第2期周宏等:用于ALE有限元模拟的网格更新方法271 图2给出了椭圆形贮箱内液体网格节点的速度 场.图3给出了对应时刻液体晃动的质量节点速度 场.可以看到在自由面上,网格节点的速度场与质 量节点的速度场都满足了自由液面的运动学边界条 件,通过测定计算网格的体积变化能检验不可压缩 条件是否得到满足,该算例液体体积保持在4-0.2% 的范围内.同时液面的运动特征反应了在晃动幅度 超过自由面特征尺寸的15%时具有的非线性特性. 不足的是,包含接触点的单元各节点都在边界上, 当流场速度值较大时,该网格的刚性位移和变形都 较周边单元大,只有局部加密单元或局部重构单元 才能获得更光滑的接触线界面区域. 0.5 O 0.5 0.5 O () 0.5 0 O O 图2椭圆形贮箱内液体晃动的速度场 Fig.2Fluidvelocityfieldsofsloshinginellipticalcontainer 0.5 0 — 0.5 0.5 0 O.5 0.5 O — O.5 0.5 0 — 0.5 0.5 0 0.5 0.5 O 图3椭圆形贮箱内液体网格节点的速度场 Fig.3Nodalvelocityfieldsintheellipticalcontainer 4结论 文章综述了近20年用于ALE有限元方法中的 网格更新技术,这些方法都在避免引进大规模的网 格计算上作了改进.本文通过修正ALE描述下的自 由面运动学边界条件,增加了自由面上网格节点运 动的自由度,同时了运动界面上的网格节点速 度;最后以求解Laplace方程的边值问题获得内部网 格的节点速度更新.文中给出了采用以上计算方案 完成的带弧面贮箱内的液体晃动的计算结果,可以 272力学2008年第4O卷 看到在晃动幅度超过自由液面特征尺寸15%之后, 自由液面依旧光滑,网格形状保持了很好的形态. 另外,在ALE有限元方法中使用更新网格速度还是 更新网格位移的方法,具有问题依赖性.在液一固耦 合问题中,界面上直接传递固体的变形量,因此位移 更新方法更适用.可见,无论使用哪一类网格更新方 法,都是使ALE描述方法获得更广泛应用的一项关 键技术.如果将该网格更新方法向三维情况推广, 将能够在弧形贮箱内完成如文献【13~15]中实现的 液体晃动仿真计算,有助于对液体晃动问题的进一 步认识. 参考文献 1HughesTJR.LiuWK.Lagrangian—Eulerianfiniteelement formulationforincompressibleviscousflows.Comput MethsApplMechEngrg,1981,29:329~349 2DoneaJ,GiulianiS.AnarbitraryLagrangian—Eulerianfi— niteelementmethodfortransientdynamicfluid—structure interaction.ComputMethsApplMechEngrg,1982,33: 689723 3HuertaALiuWK.Viscousflowwithlargefreesurfacefluid flow.ComputMethsApplMechEngrg,1988,69:277~324 4YamamotoK,KawaharaM.Structuraloscillationcontrol usingtunedliquiddamper.ComputStruct,1999,71: 435~446 5FarhatC.Parallelanddistributedsolutionofcouplednon— lineardynamicaeroelasticresponseproblems.InParal— lelSolutionMethodsinComputationalMechanics,Pa— padrakakisM,ed.NewYork:JohnWiley,1988.243~301 6JohnsonAA.TezduyarTE,Meshupdatesstrategiesin parallelfiniteelementcomputationsofflowproblemswith movingboundariesandinterfaces.ComputMethsAppl MechEngrg,1994,119:73~94 7WaIIWA.RammE.Fluid.structureinteractionbasedupon astabilized(ALE1finiteelementmethod.Proc.4thWl0rid CongressonComputationalMechanics—BuenosAires,Idel— sohnSR,OnateE,DvorkinEN,eds.Barcelona,Spain: CIMNE.1998 8SteinetaK.Automaticmeshupdateswiththesolid. extensionmeshmovingtechnique.ComputMethsAppl MechEngrg,2004,193:2019~2032 9ChiandussiG,BugedaG,OnateE.Asimplemethodfor automaticupdateoffiniteelementmeshes.Commun?. “erMethEngng.2000.16:119 10LohnerR.YangC.ImprovedALEmeshvelocitiesformov. ingbodies.CommunicationsinNumericalMethodsinEn— gineering,1996,12:599~608 11TezduyarTE,BehrM,MittalS.eta1.Computation ofunsteadyincompressibleflowswiththefiniteelement methods--space.timeformulations.iterativestrategiesand massivelyparallelimplementations,newmethodsintran. sientanalysis.SmolinskiPLiuWKHulbertG.eds.New York:ASME.1992AMD一143.724 12DodgeFT.Thenew”Dynamicbehaviorofliquidsinmov— ingcontainers”.SWRI.2000 13岳宝增,李俊峰.三维液体非线性晃动及其复杂现象.力学学 报,2002,34(6):949~955(YueBaozeng,LiJun~ng.The threedimensionalliquidnonlinearsloshinganditscomplex phenomena.ActaMechanicaSinica,2002,34(6):949~955 (inChinese)) 14郭正,刘君,瞿章华.非结构动网格在三维可动边界问题中的 应用.力学,2003,35(2):140~146(GuoZheng,Liu Jun,QuZhanghua,Dynamicunstructuredgridmethod withapplicationsto3Dunsteadyflowsinvolvingmoving boundaries.ActaMechanicaSinica,2003,35(2):14O146 (inChinese)) 15岳宝增.俯仰激励下三维液体大幅晃动问题研究.力学, 2005,37(3):199~203(YueBaozeng.Threedimensional largeamplitudeliquidsloshingunderpitchingexcitation. ActaMechanicaSinica,2005,37(3):199~203(inChinese)) MESHUPDTEALG0RITHMINALEFINITEMETH0DWITHINFREE SURFACEFL0W) ZhouHongLiJunfengWangTianshu.) (SchoololAerospace,TsinghuaUniversity,Beijing100084,China) AbstractThearbitraryLagrange—Eulermethodaffordstrackingthemotion offreeinterfaceinliquidslosh problemsthroughdefiningthenodalvelocitiesofreferenceframe.Basedonreviewingthethreeelementary algorithmsappliedinmeshupdating,themethodthatcomputationalnodescanmovemorethanonedirection isusedinspecialnumericalexample.ThisachievementdemonstratesthattheALEmethodisavailableformore complexgeometricalboundaries.Thenodalvelocitiesandmeshconfigurationsduringthesloshinginhorizontal cylinderareillustratedatthesectionofnuemericalexample. KeywordsALEmethod,finiteelementmethod,meshupdate,liquidslosh,horizontalcylindercontainer Received29January2007,revised21December2007. 1)TheprojectsupportedbytheNationalNaturalScienceFoundationofChina (10572022,10302013) 2)E-mail:tswang@tsinghua.edu.ca