任意拉格朗日—欧拉法

任意拉格朗日—欧拉法 任意拉格朗日—欧拉法 第14卷第1期 计算力学学报 Vol.14No.1任意拉格朗日2欧拉描述法研究进展 张 雄 陆明万 王建军 (清华大学工程力学系,北京,100084)Ξ 摘 要 任意拉格朗日2欧拉(ALE)描述综合了纯拉格朗日和纯欧拉描述的优点,克服了各 自的缺点,成为非线性连续介质力学中大变形的非常有效的方法。本文论述了ALE法 的研究进展及其在流体动力学、流体2结构相互作用、加工成型、碰撞、接触等大变形问中的 应用。 关键词 非线性连续介质力学;任意拉格朗日2欧拉描述;非线性有限元 分类号 O33;O35;O242121 1 引 言 非线性连续介质力学的有限元分析方法已取得了很多令人瞩目的进展,其中在大变形问题中绝大多数研究工作都采用拉格朗日(主要用在固体力学中)或欧拉(主要用在流体力学中)描述方法。 在拉格朗日描述(即物质描述或L描述)中,计算网格固定在物体上——————————————————————————————————————————————— 随物体一起运动,即网格点与物质点在物体的变形过程中始终保持重合,因此物质点与网格点之间不存在相对运动(即迁移运动,也称对流运动)。这大大地简化了控制方程的求解过程,而且能准确描述物体的移动界面,并可跟踪质点的运动轨迹。但在涉及到特大变形的问题中,物质的扭曲将导致计算网格的畸形而使得计算失败。 在欧拉描述(即空间描述或E描述)中,网格固定在空间中,即计算网格在物体的变形过程中保持不变,因此可很容易处理物质的扭曲。但对运动界面需要引入非常复杂的数学映射,将可能导致较大的误差。另外当使用普通的伽辽金离散时,由于迁移项的影响,有限元方程中的系数矩阵是非对称的,而且还可能得到振荡解。 纯拉格朗日和纯欧拉描述都存在严重的缺陷,但也具有各自的优势。如果能将二者有机地结合,充分吸收各自的优势,克服各自的缺点,则可解决一大批只用纯拉格朗日和纯欧拉描述所解决不了的问题。任意拉格朗日2欧拉(ALE)方法就是基于此目的最早由Noh(1964)[1]以耦合欧拉2拉格朗日的术语提出的,并用有限差分法求解带有移动边界的二维流体动力学问题。在Noh的研究工作中,网格点可以随物质点一起运动,但也可以在空间中固定不动,甚至网格点可以在一个方向上固定,而在另一个方向上随物体一起运动,因此ALE描述也被称为耦合欧拉2拉格朗日描述。例如在液体表面波的传播问题中,网格点在垂向随物质点一起运动,而在Ξ核工业科学基金资助。 本文于1996年5月27日收到。 ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights ——————————————————————————————————————————————— reserved. 水平方向上固定不动。这样可很容易描述液体表面的运动,而且网格不会发生扭曲。 在ALE描述中计算网格可以在空间中以任意的形式运动,即可以独立于物质坐标系和空间坐标系运动。这样通过规定合适的网格运动形式可以准确地描述物体的移动界面,并维持单元的合理形状。类似于E描述,在ALE描述下的控制方程中也将出现对流项,因此也可能得到振荡解,需要进行相应的数值处理。纯拉格朗日和纯欧拉描述实际上是ALE描述的两个特例,即当网格的运动速度等于物体的运动速度时就退化为拉格朗日描述,而当网格固定于空间不动时就退化为欧拉描述。 [128]ALE最早是为了解决流体动力学问题而引入的,并且使用有限差分法。由于核反应堆 [13,14]结构安全分析的需要,Donea[12]、等人分别将ALE法引入有限元法中,用以求Belytschko [9]解流体与结构相互作用问题[12230]。Hughes等人建立了ALE描述的运动学理论,并使用有限 元法解决了粘性不可压缩流体流动和自由表面流动问题。 ALE描述中,参考构形是已知的,而初始构形和现时构形都是待求解的。因此ALE法尤其适合于在初始构形和现时构形都未知的问题中使用,如接触问题等。 目前,ALE法已被应用于固体力学领域中求解大变形问题,如碰撞——————————————————————————————————————————————— [36][55]、接触和弹性断裂力学[31,34]、路径相关材料(如弹塑性材料等)[35,39]、加工成型[40,47]等。另外对网格的运动以及控制方程的求解策略等问题也进行了多方面的研究[53,58]。 本文论述了ALE法的基本理论和研究进展及其在流体动力学、流体2结构相互作用、加工成型、碰撞、接触等大变形问题中的应用。 2 ALE运动学描述 将连续体在初始时刻t0的构形(称初始构形)记为8X,将t时刻的构形(称现时构形)记为8x。为了确定各物质点的位置,引入Lagrange或物质坐标系OX1X2X3(简记为L坐标系),它是和物体固结在一起并随物体一起运动的。质点在初始时刻to的位置由其在L坐标系中的位置矢量X完全确定,因此X也可用来识别物体中的不同质点。为了描述现时构形,引入空间坐标系或Euler坐标系ox1x2x3(简记为E坐标系),它是和空间固结在一起的。空间中各点的位置由其在E坐标系中的位置矢量x确定,因此可用x表示空间中的几何点。 拉格朗日描述以初始构形为参考构形研究物质点X在空间中的运动规律,即: x=x(X,t)(11a) 而欧拉描述则是以现时构形为参考构形来研究空间点x上物质点的运动规律,即研究: (11b)X=X(x,t) 式(11a)描述了同一质点X在不同时刻的空间位置,而式(11b)则描述了同一空间点x在不同时刻被物质点占据的情况。 ——————————————————————————————————————————————— 不同于拉格朗日和欧拉描述,ALE描述另外引入了一个可以独立于初始构形和现时构形运动的参考构形,记为8Ν。在物体的变形过程中,观察者始终跟随参考构形运动,因而对观察者而言参考构形是固定不动的,而初始构形和现时构形则都相对于参考构形运动。为了确定参考构形中各参考点的位置,引入参考坐标系oΝ1Ν2Ν3,参考构形中各点的位置由其在参考坐标系中的位置矢量Ν确定,因此 Ν=Ν(X,t) 描述了物质点X在参考坐标系中的运动规律,而 ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.(11c) x=x(Ν,t)(11d) 则描述参考点Ν在空间中的运动规律。 在物质描述方法中的有限单元剖分是对物体进行剖分的,网格点就是物质点,即网格是随物体一起运动的。在空间描述方法中的有限单元剖分是对空间进行的,网格点就是空间点,因此网格是固定在空间中不随物体运动的。而在ALE描述方法中的有限单元剖分是对参考构形进行的,网格点就是参考点,网格是独立于物体和空间运动的,可以根据需要自由选择。 以上映射应该是一对一映射,雅可比行列式J(描述从初始构形到现时构形的映射关系)和混合雅可比行列式Jδ(描述从参考构形到现时构形的映射关系)都应不等于0,即: J=?????0, Jδ=?????0Xjj(2) ——————————————————————————————————————————————— 根据定义t时刻某质点X在空间中的运动速度v等于质点X在空间的位置矢量x(X,t)对时间的导数,即: (3)v=??Xtδ则等于参考点(即网 参考构形中某点Ν在空间中的运动速度(也就是网格点的运动速度)v 格点)Ν在空间中的位置矢量x(Ν,t)对时间的导数,即 δ(4)v=??Νt 物质点X在参考坐标系中的位置矢量Ν(X,t)对时间的导数 (5)w=??Xt 则为物质点X在参考坐标系中的速度。 在ALE描述中,参考构形(也就是计算网格)的运动规律可以是任意给定的,指定特殊的网格运动规律可以将ALE描述退化为L描述和E描述。δ=0,即计算网格在空间中固定不动,退化为欧拉描述。?? vδ=v,即计算网格随同物体一起运动,退化为拉格朗日描述。?? vδ?v?0,即计算网格在空间中独立运动,对应于一般的ALE描述。?? v 在ALE描述中参考构形是已知的,各物理量用参考坐标Ν描述比较方便,即F=F(Ν,t)。因此各物理量的物质导数(即某固定质点X的物理量对时间的变化率)应通过该物理量的参考导数(即某固定参考点Ν处的物理量对时间的变化率,也称混合导数)来计算。为了表达方便,本文采用指标求和法则。对x求物质导数,得: δ(6)vi=vi+??wjj 将式(6)整理得: (7)ci=jjδ为物质点相对于网格点的运动速度,即迁移速度(也称对——————————————————————————————————————————————— 流速度)。式中:ci=vi-v从式(7)可i 以看出,迁移速度c实质上是将质点在参考坐标系中的运动速度w变换到现时坐标系中而得到的,也就是说c是在现时坐标系中物质点相对于参考构形的运动速度。式(7)也可化为: (8)ci=??Xt ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. 利用式(7)可将任意物理量F(Ν,t)的物质导数表示为: (9)??X=??Ν+cittxi 式(9)是ALE描述中很重要的关系式,利用它可将控制方程转化到参考坐标系下求解。δALE描述的重要特征是可以根据需要给定合适的网格运动速度v,以维持计算网格的合i理形状并准确描述物体的移动界面。所以网格运动算法在ALE描述中占很重要的地位。3 ALE描述的控制方程 为了建立在ALE描述下的控制方程,利用非线性连续介质力学的方法以及高斯定理可以推出物理量 G(t)=?g(Ν,t)dV 和 H(t)=?h(x,t)dV8ΝΝ8xx(10) 的物质导数为: ??X=t ??Xt 311 质量守恒方程(连续性方程)?(??=?t8Ν()dV??Ν+tixΝ(111a)(111b)8x+)dVxix ——————————————————————————————————————————————— 考察一连续体,用58X、58x和58Ν分别表示物质域8X、空间域8x和参考域8Ν的边界,而用δ分别表示连续体各构形的密度。ΘΘ和Θ在不同构形中该连续体的质量M可以写为:0、 M= δ(Νδ(x,t)式中: Θ,t)=JΘ(x,t), Θ0(X,t)=JΘ 质量守恒定律表明,质量的整体变化率(即物质导数)应等于零。利用式(12)和式(111a)可以得到在参考坐标系下的质量守恒方程为:δδ(131a)??Ν+=0 在8Ν域中ti 类似于修正拉格朗日(UL)方法,有时在空间域8x中求解控制方程是比较方便的。利用式(12)、式(111b)和式(9)可以把质量守恒方程写成:?8ΝδdV=ΘΝ?ΘdV=?ΘdV8xx8X0X(12) ??Ν+ci+Θ=0 在8x域中5t5xi5xi 312 动量守恒方程(平衡方程)(131b) 动量守恒定律表明,在t时刻占据参考域8Ν的物体的总动量的整体变化率等于施加在物体上的外力之和,即: ??Xt?8ΝδvdV=ΘiΝ?58ΝδtidSΝ+?8ΝδfdVΘiΝ(14) 式中δti是作用在参考域8Ν的边界58Ν的单位表面上的力,fi是作用于物体中单位质量的体力。利用散度定理以及式(111a)、式(131a)可将上式写成: δδδδf 在8域中(151a)Θ??Ν+Θwj=+ΘiΝtjj ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. ——————————————————————————————————————————————— 式中Tδji是定义在参考构形下的第一类皮奥拉2克希荷夫应力张量(即拉格朗日应力张量)。它与真实应力(即柯西应力)张量Ρij之间的关系为: δ=Jδ(16)Ρkjijxk 如果在参考域中直接求解式(151a),Tδji的不对称性将会给求解带来一定的困难。由于定T义在现时构形中的柯西应力是真实应力,因此在某些情况下将该式转化到空间域中求解是方便的。将式(151a)两边同乘Jδ-1后整理可得: Θ 313 能量守恒方程??Ν+Θcj=+Θfi 在8x域中txjxj(151b) 能量守恒定律指出,能量(包括内能和动能两部分)对时间的整体变化率应等于单位时间内外力对物体所作的功与单位时间内传入物体的热量之和(这里没有考虑物体从外界得到的其它能量,如电磁能等)。即: δEdV=δδfvdV-(17)??XΘtividSΝ+ΘqinidSΝΝiiΝ8Ν58Ν8Ν58Νt???? (131a)、(151a)和式(9)以及式中E是物体中单位质量上的能量,qi为热通量。利用式(111a)、 散度定理可把上式写为: δ??+Θδδ-Θwi=TjiΝtiji 式中e为物体单位质量中的内能。同样可将上式转化到现时构形中,即: ——————————————————————————————————————————————— Θ??Ν+Θci=Ρij-txixjxi(181a)(181b) 事实上能量方程只有当存在不同类型的能量转化时才起作用。 在上述控制方程中都出现了与迁移速度有关的迁移项,因此ALE方法同样具有类似于欧拉方法的缺点。但对某些特殊情况有可能通过选择合适的网格运动方式来降低迁移速度以克服数值求解困难,如文献[35]中的一维弹塑性波的传播问题。 由于ALE描述下的控制方程在形式上同欧拉描述中的控制方程很相似,因此有的研究者也将ALE法称为拟欧拉方法[13,14]。另外因为可以把从初始构形到现时构形的映射分解为从初始构形到参考构形的映射和从参考构形到现时构形的映射两部分,所以ALE法有时也被称为混合欧拉2拉格朗日法[34]。 314 边界条件 事实上,边界条件只与所求解的问题有关,而与所使用的描述无关。因而在欧拉描述和拉格朗日描述中使用的边界条件都可用于ALE描述中。动力学边界条件仍然是自然边界条件,它在求解时自动满足。这里主要讨论运动学边界条件。 在欧拉描述中,对移动边界的描述需引入非常复杂的数学映射,而在拉格朗日描述中则很简单。引入ALE描述的目的之一就是为了解决这一问题的,因此可在物质表面上使用拉格朗日描述。即:δ(191a)v=v 在物质表面上 上式是物质表面应满足的边界条件。但在涉及到大变形流动(如——————————————————————————————————————————————— 流体流动及金属加工成型等)问题时,为了既能准确描述物体的边界,又能维持网格的合理形状,一般只在物质表面的法向 ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. 使用物质描述,而允许网格沿切向运动。即: c??nx=0或w??nΝ=0(191b) 式中nx和nΝ分别是物体的现时构形和参考构形的表面 外法向矢量。上式需要在自由表面上求解以确定自由表 面的位置。 另一类常见的边界是物体间的交界面,如流体2结 构相互作用问题和摩擦接触问题。在界面上的每个点处 一般需设置两个节点,它们分别属于各自的物体,如图1 所示。在交界面上两物体可以沿切向滑动、相互粘接或 相互脱离。如果两物体完全粘接,则界面上节点对的物 质速度应相等。即: v??=v??图1 物体的界面(201a) (201b) 如果两物体沿界面的切向滑移,则界面上节点对的法向物质速度应相等,即:(v??-v??)??n=0 式中n为界面的法向矢量。如果两物体相互脱离,则界面上节点对的物质速度是相互独立的。 在ALE格式中,计算网格的选取可以独立于物体自身的运动。因此可令界面上节点对的网格运动速度相等,使得界面上的节点对在物——————————————————————————————————————————————— 体的变形过程中始终保持重合,以简化流体2结构相互作用和摩擦接触等问题中界面的处理过程。即:δδv=v????(21) 式(21)使得ALE描述很好地克服了拉格朗日法处理接触问题和流体2结构相互作用问题时存在的困难。在流体2结构相互作用问题中的流固界面上,固体界面一般采用物质边界条件(191a),而流体界面则采用条件(191b),即对流体域只在表面的法向使用物质描述,而在表面的切向上网格的运动与流体自身的运动无关,以维持流体域网格的合理形状。 在ALE描述的控制方程中存在的与网格和物体之间的迁移速度有关的迁移项将会给方程的求解带来很大的困难。如果使用普通的Galerkin有限元离散,则得到的有限元方程的系数矩阵将是非对称的,而且当网格与物体之间的迁移速度较大时,将可能得到空间振荡解。虽然可使用很精密的网格来消除解的空间振荡,但却大幅度地提高了计算量。 这个困难最初是在有限差分法中通过对迁移项使用迎风差分[62]的方法解决的。七十年代这一思想被Heinrich等人引入有限元法中,称为迎风有限元法[63]。另外,Petrov2Galerkin[64](将权函数和试函数取为不同的函数)和Taylor2Galerkin法[6]也是解决解的空间振荡问题的有效方法。 只有当引入物质的本构关系后,整个控制方程才能求解。根据物质客观性原理可知,本构方程是由物质决定的,而与所采用的运动学描述无关。因此在ALE描述中物质的本构方程与其它描述中的本构方程完全相同,可以直接引用。 ——————————————————————————————————————————————— 4 网格运动 在ALE描述中,网格可以自由运动这一特性给大变形问题的求解带来了很大的方便,但也给求解带来了困难,增加了计算量。因此如何选取合理的、有效的网格运动算法是ALE描 ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. 述中的一个很重要的问题。δ或控制质点在参考系下的运动速度w来实现网格的运动是通过控制网格的运动速度v 的,运动学边界条件对网格的运动方式有很大的影响。在某些情况下,边界的运动是已知的,这时可根据所求解问题的特殊性预先指定网格的运动速度或位移。如Liu等人[35]在求解弹塑性压力波的传播问题时采用了匀速网格运动以消除解的振荡。 [9]Hughes等人在研究粘性流动问题时提出了拉格朗日2欧拉矩阵法。他们把w表示为: (22)wi=(?ij-Αij)vj 式中Α。当Αij是对角矩阵,称为拉格朗日2欧拉矩阵ij=?ij时得到拉格朗日描述,当Αij=0时得到欧拉描述,而当Α。ij取其它值时,网格可以任意运动δ的微分方程组。将式(22)代入式(7)后可得到一组关于网格运动速度v在求解该微分方程组时需要将参数Α。解方程求得每一个节ij在求解域内进行有限元离散,并建立有限元方程 点处的参数后即可确定网格的运动规律。 在求解自由表面流动问题时,为了能准确描述自由表面的运动,一——————————————————————————————————————————————— 般在自由表面的法向(假设为Ν。同时为了避免网格的扭曲而在切向上(Ν3)上使用物质描述,即Α33=11和Ν2)使用了欧拉描述,即Α也就是说网格只能沿自由表面的法向运动,因此只需由式(7)和11=Α22=0。δ。(22)求解Ν方向的网格运动速度v拉格朗日2欧拉矩阵法需要在整个区域内求解网格运动δ,计算量比较大。速度v另外只通过给定拉格朗日2欧拉系数来维持网格的合理形状是比较困3 难的,且不直观。因此Huerta等人引入了混合法[11]。 不同于拉格朗日2欧拉矩阵法,混合法在自由表面的法向Ν3上使用物质描述(满足边界条δ和件(191b),但在切向Ν和Ν上需根据具体问题给定(或与其它节点耦合)网格的运动速度vv2。当由式(7)求得自由表面的法向网格运动速度后,区域内部的网格运动可用其它方法确定,33δ121如求解拉普拉斯方程[43,49]或直接用边界上的网格运动速度进行线性插值[23]等。 [55]Benson采用了一种简单但又很有效的网格运动算法。网格的扭曲一般有两种形式,即 压缩扭曲和剪切扭曲。当某节点周围最小单元的面积和最大单元的面积之比小于某一给定值或当该节点周围单元的顶角小于某一给定值时,则认为该节点周围的单元发生畸形,该节点需要移动以调整网格的形状。而节点的移动方法可采用等势松驰法[59,60]等方法。 Giuliani提出了另外一种度量单元压缩扭曲和剪切扭曲的指标[58],并将这两个指标的平方和作为评价网格形状好坏的目标函数。最优的网格应该使得这个目标函数取最小值,因此可以利用最优化方法——————————————————————————————————————————————— 寻求最优的网格点位移。该算法已在瞬态流体2结构相互作用计算程序EU2 [53,54]和Schreurs等人[40,44]也提出了相应的网格点运动算RDYN21M中实现。另外Ponthon 法。 在拉格朗日描述中,当网格出现扭曲后,通常需要重新对网格进行划分。因为方程是建立在网格点上的,因此当网格重新划分后,所有变量都要进行插值以求其在新网格点上的值。这种方法将带来额外误差。在ALE描述中的每一个时间步上网格都可以运动,但不需要插值运算来处理迁移运动,因为这一项内容已经包括在控制方程中。 5 流体动力学和流体2结构相互作用问题 流体动力学问题一般采用欧拉描述,对移动边界需要进行复杂的处理。1964年Noh[1]首? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. 次提出了ALE描述的思想,并引入有限差分程序CEL中,用来求解带移动边界的二维流体动力学问题。虽然Noh建立了在可以自由运动的参考构形中的控制方程,但在程序实现时实际上假定各计算区域只能使用拉格朗日描述和欧拉描述,而不能使网格自由运动。Trulio[2]研制 [3,4]了用ALE求解可压缩流体流动问题的有限差分程序AFTON。将ALEHirt、Amsden等人 与隐式积分格式结合解决了具有任意流动速度的二维流体流动——————————————————————————————————————————————— 问题,并研制了相应的有限差 [5][6]分程序YAQUI。Pracht、Stein等人将ALE推广到具有任意形状的移动边界的三维流体流 [7,8]动问题中。等人研制了二维、三维简化ALE格式(认为每个单元只有一种材料组成Amsden 且在边界使用物质描述)的有限差分程序SALE和SALE23D。 在这些研究工作中都使用有限差分法,并且只考虑了无粘性可压缩流体。Hughes、Ra2maswamy、Huerta[9,11]等人建立了ALE描述下的运动学理论,并解决了粘性不可压缩流体流动和自由表面流动问题。Huerta和Liu[11]用ALE描述的Petrov2Galerkin有限元很好地解决了海啸、坝的决口、流体的大幅度晃动等流体动力学问题。 [12,14]基于核反应堆结构安全分析的需要,Donea、将ALE描述法引入了Belytschko等人 有限元法中,用以求解瞬态流体2结构相互作用问题,较好地解决了容器中液面的大幅度晃动[19,24]、液体中高压气泡的扩展[20,23]、边界的大幅度运动[21,23]和容器中的液体泄露[23]等问题。 在流体2结构相互作用问题中,流体的运动幅度比较大。因此在固体域中一般使用拉格朗日描述,而在流体域内使用欧拉描述。在流体2结构交界面上的每个点处需设置两个节点,即流体节点和固体节点。固体节点应满足物质表面条件(191a),而流体节点应满足界面条件(191b)、(20)和(21)。也就是说在整个变形过程中,在流体2结构交界面上流体质点不能穿越网格,而且也不能穿越固体表面,同时应使固体网——————————————————————————————————————————————— 格点和流体网格点在界面上始终保持重合。6 固体力学问题 由于几何复杂性、大变形、有限滑移、材料非线性等因素的影响,固体间的接触问题是很复杂的。有限元法是分析接触问题的有效方法,主要有界面元法、变分不等式或参变量变分原理以及应力杂交元法等。用有限元求解接触问题时,在界面上需设置分别属于各自物体的两个节点。为了准确满足界面相容条件和计算接触应力,不论两个物体是粘接还是相互滑移,在变形过程中节点对必须保持重合。在拉格朗日描述中,物体的接触区域在求解前是未知的,因此很难使接触区的节点对在变形过程中始终保持重合。为了求解接触问题,一般需要使用非常精密的网格,而且对有限滑移问题还需求解单边位移约束的最优化问题。 在ALE描述中物体的初始构形和现时构形都是需要求解的部分,因此ALE法很好地克服了这个困难。Haber等人将位移矢量分解为拉格格朗日位移(实质上是网格点的位移)和欧拉位移(实质上是物质点与网格点之间的相对位移)之和[31]。在接触面上使节点对的拉格朗日位移相同,但各自可以具有独立的切向欧拉位移。这样在变形过程中接触面上的节点对始终保 (20)、(21)是完全等价持重合,而物质点可以沿切向相对于网格点运动。这种方法与式(19)、 的,很容易处理物体之间的摩擦接触[31,34]、加工成型[41,47]及碰撞[36,49,55]等大变形问题。 在弹性断裂力学中,应力强度因子可以用多种方法求得,如奇异单——————————————————————————————————————————————— 元法、J积分法、能量释放率法等。在拉格朗日位移法中,裂纹的长度不能连续变化,因此一般使用数值微分来求解由裂纹扩展引起的能量释放率。而在ALE法中,裂纹的长度是可以连续变化的。文献[32]导出了由裂纹引起的能量释放率的显式公式,将ALE应用于弹性断裂力学中求解应力强度因子,得 ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. 到了相当好的结果[34]。 对胡克类固体和牛顿流体等线性路径无关材料,其应力状态仅由各时刻的位移场和速度场唯一确定。与格拉朗日描述不同,ALE描述中的网格点在变形过程中并不代表同一物质点。因此当使用ALE描述求解记忆材料等路径相关材料时,单元的状态变量将受到物质点的迁移运动的影响[35,38]。这些物质点具有不同的应力和应变历史。这个问题同样存在于用欧拉描述分析塑性加工成型过程中。Liu等人建立了路径相关材料ALE描述的一般格式和非线性 [35]ALE有限元分析的显式计算过程。 7 结 语 ALE法吸收了拉格朗日描述和欧拉描述的优势,克服了各自的缺陷,因此成为非线性连续介质力学中很有前景的分析方法。由于ALE描述的计算网格可以独立运动,使得ALE描述成功地解决了一大批自由表面问题、摩擦接触问题、流体2结构相互作用问题。但ALE法使得计算较为复杂,而且也存在迁移的影响,因此需要研究简单的、有效——————————————————————————————————————————————— 的ALE计算格式(包括网格运动算法、方程的求解策略、迁移影响的 处理等)。文献[53,58]在这方面作了很多研究。 参 考 文 献 1 NohWF.CEL:Atime2dependenttwo2space2dimensionalcoupledEulerian2lagrangiancode,in:B. Alder,S.FernbachandM.Rotenberg,eds.MethodsinComputationalPhysics3,(Academicpress,NewYork,1964) 2 TrulioJG.TheoryandstructureoftheAFTONcodes,Rept.ASWL2TR266219,AirForceWeapons Laboratory,1966 3 AmsdenAAandHirtCW.YAQUI:AnarbitraryLagrangian2Euleriancomputerprogramforfluidflow atallspeeds,ReportLA25100,LosAlamosScientificLaboratory,LosAlamos,NM,1973 4 HirtCW,AmsdenAAandCookJL.AnarbitraryLagrangian2Euleriancomputingmethodforallflow ——————————————————————————————————————————————— ,253speed,J.Comput.Phys.1974,14:227 5 PrachtWE.Calculatingthree2dimensionalfluidflowsatallflowspeedswithanEulerian2Lagrangian ,159computionmesh1J.Comput.Phys.1975,17:132 6 SteinLR,GentryRAandHirtCW.Computationalsimulationoftransientblastloadingonthree2di2 , 74mensionalstructures,Comput.Meths.Appl.Mech.Engrg.,1977,11:57 7 AmsdenAA,RuppelHM,etal1SALE23D:AsimplifiedALEcomputerprogramforfluidflowatall flowspeed.LosAlmosScienticLaboratory,LosAlamos,NM,1980 8 AmsdenAA,RuppelHM,etal.SALE23D:AsimplifiedALEcomputerprogramforcalculatingthree2 dimensionalfluidflow,LosAlmosScienticLaboratory,LosAlamos,NM,1981 9 HughesTJR,LiuWKandZimmermanTK.Lagrangian2Eulerianfiniteelementformulationforincom2 ——————————————————————————————————————————————— , 349pressibleviscousflows,Comput,Meths.ApplMech.Engrg.,1981,29:329 10 RamaswamyB,KawaharaM.ArbitraryLagrangian2Eulerianfiniteelementmethodforunsteady,con2 vective,incompressibleviscousfreesurfacefluidflow,Int. 1053,1075 11 HuertaAandLiuWK.Viscousflowwithlargefreesurfacemotion.Comput.Meths.Appl.Mech.En2J.Num.MethodsFluids,1987,7: ,324grg.,1988,69:277 ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. 12 DoneaJ,Fasoli2Stellap,etal1LagrangianandEulerianfiniteelementtechniquesfortransientfluid2 2,USA,1977,15, 19structureinteractionproblems,in:Trans.SMiRT24,PaperB1?? 13 BelytschkoTandKennedyJM.Computermodelsforsubassemblysimulation1Nucl.Engrg.Design, 1978,49:17,38 ——————————————————————————————————————————————— 14 BelytschkoT,KennedyJM,etal.Quasi2Eulerianfiniteelementformulationforfluid2structureinterac2 ,69tion.J.PressureVesselTechnol.ASME,1980,102:62 15 ChangYW,ChuHY,etal.EvaluationofLagrangian,Eulerian,ArbitraryLagrangian2Eulerianmeth2 .SMiRT25,PaperB1??1,Ger2odsforfluid2structureinteractionproblemsinHCDAanalysis,in:Trans many,1979 16 DoneaJ.ArbitraryLagrangian2Eulerianfiniteelementmethods,in:T.BelytschkoandT.J.R.Hugh2 , 516es,eds.ComputationalMethodsforTransientAnalysis,North2Holland,Amsterdam,1983:473 17 LiuWK.Finiteelementproceduresforfluid2structureinteractionswithapplicationtoliquidstorage ,238tanks1Nucl.Engrg.Design,1981,65:221 18 LiuWKandMaDC.Computerimplementationaspectsforfluid2structureinter——————————————————————————————————————————————— actionproblems1Com2 .Meths.Appl.Mech.Engrg.,1982,31:129,148put 19 LiuWKandGvildysJ.Fluid2structureinteractionoftankswithaneccentriccorebarrel1Comput. ,57Meths.Appl.Mech.Engrg.,1986,58:51 20 BelytschkoT,FlanaganDPandKennedyJM.Finiteelementmethodswithuser2controlledmeshesfor , 688fluid2structureinteraction,Comput.Meths.Appl.Mech.Engrg.,1982,33:669 21 HuertaAandLiuWK.ALEformulationforlargeboundarymotion,In:Trans.SMiRT210,Vol.B, ,346USA,1989:335 22 NomuraT.ALEfiniteelementcomputaionsoffluid2structureinteractionproblems1Comput.Meths. ,308Appl.Mech.Engrg.,1994,112:291 23 NitikitpaiboonCandBatheKJ.AnarbitraryLagrangian2Eulerianvelocitypote——————————————————————————————————————————————— ntialformulationforflu2 5):871, 891id2structureinteraction1ComputersandStructures,1993,47(4?? 24 LiuWKandBelytschkoT.Fluidstructureinteractionwithsloshing,in:Trans.SMiRT27,Vol.B, ,18Chicago,IL,1983:11 25 DoneaJ,Fasoli2StellaP,etal.AnarbitraryLagrangianEulerianfiniteelementprocedurefortransient 3,Trans.SMiRT25,Germany,1979dynamicfluid2structureinteractionproblems,PaperB1?? 26 DoneaF,GiulianiS,etal.AnarbitraryLagrangian2Eulerianfiniteelementmethodfortransientdynam2 .Comput.Meths.Appl.Mech.Engrg,1982,33:689, 723icfluid2structureinteractions 27 NomuraT,HughesTJR.AnarbitraryLagrangian2Eulerianfiniteelementmethodforinteractionof , ——————————————————————————————————————————————— 138fluidandarigidbody.Comput.Meths.Appl.Mech.Engrg.1992,95:115 28 HuertaA,LiuWK.Viscousflowstructureinteraction.J.PressureVesselTechnol.ASME110,1988: 15,21 29 DoneaJ,Fasoli2StellaP,etal.ThecomputercodeEURDYN21Mfortransientdynamicfluid2structure interaction,EUR6751EN,commissionoftheEuropeancommunities,JointResearchCentIspraEstab2lishment,Italy,1980 30 LiuWK.Developmentoffiniteelementproceduresandcomputerimplementationaspectsinfluid2struc2 .eds.Newconceptsinfiniteelementanalysis,1981:247, 268tureinteractions,in:T.J.R.Hughesetal 31 HaberRBandAbelJ.Contactslipanalysisusingmixeddisplacements1J.Engrg.Mech.ASCE, 1983,109(2):411,429 32 HaberRBandKohHM.Explicitexpressionsforenergyreleaseratesusingvirtual——————————————————————————————————————————————— crackextensions1 ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. ,315Int.J.Num.Meth.Engrg.,1985,21:301 33 HaberRBandHariandjaBH.AnarbitraryLagrangian2Eulerianfiniteelementapproachtolargedefor2 , 201mationfrictionalcontact.Computers&Structures,1985,20:193 34 HaberRB.Amixedeulerian2Lagrangiandisplacementmodelforlarge2deformationanalysisinsolidme2 ,292chanics1Comput.Meths.Appl.Mech.Engrg,1984,43:277 35 LiuWK,BelytschkoT,etal.AnarbitraryLagrangian2Eulerianfiniteelementmethodforpath2depen2 , 245dentmaterials1Comput.Meths.Appl.Mech.Engrg1,1986,58:227 36 LiuWK,ChangH,etal.ArbitraryLagrangian2EulerianPetrov2Galerkinfiniteelementsfornonlinear ,310continua.Comput.Meths.Appl.Mech.Engrg.,1988,68:259 ——————————————————————————————————————————————— 37 ChenJS,LiuWK,etal.ArbitraryLagrangian2Eulerianmethodsformaterialswithmemoryandfric2 ,46tion.RecentDev.Comput.FluidDynam.AMD1989,95:41 38 GhoshS,KikuchiN.AnarbitraryLagrangian2Eulerianfiniteelementmethodforlargedeformation .Comput.Meths.Appl.Mech.Engrg.,1991,86:127, 188analysisofelastic2viscoplasticsolids 39 YamadaTandKikuchiF.AnarbitraryLagrangian2Eulerianfiniteelementmethodforincompressible , 177hyperelasticity.Comput.Meths.Appl.Mech.Engrg.,1993,102:149 40 SchreursPJGandVeidpausFE,etal.AnarbitraryEulerian2Lagrangianfiniteelementmodelforsim2 ulationofgeometricalnon2linearhyperelasticandelasto2plasticdeformationprocesses,in:J.F.T. .eds.NumericalMethodsinIndustrialFormingProcesses,PineridgePress,Swansea,U.Pittmanetal ——————————————————————————————————————————————— ,509K.1982:501 41 HuetinkJ.AnalysisofmetalformingprocessesbasedonacombinedEulerian2Lagrangianfiniteelement .eds.NumericalMethodsinIndustrialFormingProcesses,Piner2formulation,in:J.F.T.Pittmanetal ,509idgePress,Swansea,U.K.1982:501 42 HuetinkJ.VreedePT,etal.ProgressinMixedEulerian2Lagrangianfiniteelementsimulationofform2 .Int.J.Num.Meth.Engrg,1990,30:1441,1457ingprocesses 43 LiuWKChenJS,etal.AdaptiveALEfiniteelementswithparticularreferencetoexternalworkrate , 216onfrictionalinterface,Comput.Meths.Appl.Mech.Engrg.,1991,93:189 44 SchreursPJG,VeldpausFE,etal.Simulationofformingprocesses,usingthearbitraryEulerian2La2 , 36grangianformulation,Comput.Meths.Appl.Mech.Engrg,1986,58:19 45 ——————————————————————————————————————————————— HuetinkvanJ,derLugt,etal.AmixedEuler2Lagrangiancontactelementtodescribeboundaryandin2 .Proc.NUMETA,1987,1,Nijhof,TheNetherlandsterfacebehaviourinformingprocesses 46 HuetinkJ,VreedePT,etal.Thesimulationofcontactproblemsinformingprocessedusingamixed .eds.NumericalMethodsinIndus2Euler2Lagrangianfiniteelementmethod,in:E.G.Thompsonetal ,554trialFormingProcesses.Balkema,Roterdam,1989:549 47 LiuWK,ChangH,etal.ArbitraryLagrangian2Eulerianstressupdateforformingsimulations,Proc. Symp,RecentAdvancesinInelasticAnalysis,ASME,Boston,USA,1987 48 KennedymJM,BelytschkoTB.TheoryandapplicationofafiniteelementmethodforarbitratyLa2 ,146grangian2Eulerianandsolids,NuclEngrg.Design,1981,68:129 49 BensonDJ.AddinganALEcapabilitytoDYNA2d:ExperiencesandConclusions,Post2Conferenceon ——————————————————————————————————————————————— IMPACT9thInt.Conf.onStructuralMech.inReactorTechnology,Switzerland,1987 50 HuertaA.ProgressinarbitraryLagrangian2Euleriananalysisoffluidandsolidproblems,In:Trans. 1,Germany,1993:1,12SMiRT212,PaperB01?? 51 HuertaA,CabotGPandBodeL.ALEformulationinnonlocalstrainsofteningsolids,in:D.R.J. , 2268Owenetaleds.ComputationalPlasticity:FundamentalsandApplications.1992,2:2249 ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. 52 DoneaJT,BelytschkoT.Advancesincomputationalmechanics.Nucl.Engrg.Design,1992,134:1, 22 53 PonthonJP.AmethodtoreducecostofmeshdeformationinEulerian2Lagrangianformulation,in:J.L. , ——————————————————————————————————————————————— 74ChenotandE.Onate,eds.ModelingofMetalFormingProcesses.1988:65 54 PonthonJP.EfficientmeshmanagementinEulerian2Lagrangianmethodforlargedeformationanalysis, .eds.NumericalMethodsinIndustrialFormingProcesses.Balkema,1989:in:E.G.Thompsonetal 203,210 55 BensonDJ.Anefficient,accurate,si.Comput.mpleALEmethodfornonlinearfiniteelementprograms ,350Meths.Appl.Mech,Engrg,1989,72:305 56 BensonDJ,DavidJ.VectorizationtechniquesforexplicitarbitraryLagrangian2Euleriancalculations. , 328ComputerMethodsinAppliedMechanicsandEngineering,1992,96:303 57 ChengJH,KikuchiN.Ameshrezoningtechniqueforfiniteelementsimulationsofmetalformingpro2 .Int.J.Num.Meth.Engrg.,1989,23:219,228cesses 58 GiulianiS.Analgorithmforcontinuousrezoningofthehydrodynamicgridinarb——————————————————————————————————————————————— itraryLagrangian2Eule2 ,212riancomputercodes.Nucl.Engrg.Design.,1982,72:205 59 WinslowAM.“Equipotential” Zoningoftwo2dimensionalmeshes,UCRL27312,LawrenceRadiation Laboratory.1963 60 WinslowAMandBartonRT.Rescalingofequipotentialsmoothing,UCID219486,LawrenceLiver2 moreNationalLaboratory,Livermore,CA,1982 61 王天宏1接触问题的ALE方法研究:〔硕士学位论文〕.西安: 西安交通大学,199415 62 RitchmyerRDandMortonKW.Differencemethodsforinitial2valueproblems(Interscience,New York,2nd.ed.,1967) 63 HeinrichJC,HuyakornPS,etal.Anupwindfiniteelementschemefortwo2dimensionalconvective ,145transport.Int.J.Num.Meths.Engrg.1977,11:131 64 LohnerR,MorganK,etal.Thesolutionofnonlinearhyperbolicequationssystemsbythefiniteele2 ——————————————————————————————————————————————— ,1063mentmethod.Int.J.Num.Meths.Fluids,1984,4:1043 ResearchprogressinarbitraryLagrangian-Eulerianmethod ZhangXiong LuMingwan WangJianjun (DepartmentofEngineeringMechanics,TsinghuaUniversity,Beijing100084) Abstract ArbitraryLagrangian2Eulerian(ALE)methodincorporatedtheadvantagesofLa2grangianandEulerianmethodandovercametheirdisadvantages,soitbecameaveryefficient .Theresearchmethodforthelarge2deformationanalysisinnonlinearcontinuummechanics progressofALEanditsapplicatonsinhydrodynamics,fluid2structureinteraction,forming .processes,impactandcontactproblems Keywords:nonlinearcontinuummechanics;arbitrarylagrangian2euleriandescription;non2 linearfiniteelement ? 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. ——————————————————————————————————————————————— ———————————————————————————————————————————————