拉格朗日、欧拉和任意拉格朗日-欧拉描述的有限元分析

暨两岸船舶与海洋工程水动力学研讨会文集拉格朗日、欧拉和任意拉格朗日-欧拉描述的有限元分析孙江龙1杨文玉2 杨侠3（1 华中科技大学船舶与海洋工程学院，武汉 430074；2 华中科技大学机械科学与工程学院，武汉 430074；3 武汉工程大学机电工程学院，武汉 430073）摘要：对拉格朗日、欧拉和任意拉格朗日–欧拉三种描述方法进行了分析，为了便于理解给出了三种描述的参考构形和参考坐标系，在参考坐标系下根据物质导数的定义分别得到相应的速度和加速度，并进行比较，将三种描述方法的区别列于表中，清晰地阐述了三种描述之间的相互关系，并进行了有限元分析。

关键词：拉格朗日;欧拉;任意拉格朗日–欧拉;有限元法1 引言自由液面大晃动引起的强非线性往往给问题的求解造成很大困难，对大晃动问题进行数值模拟，要先解决描述方法的选择问题。

过去通常采用欧拉法[1-3]和拉格朗日法[4-5]来描述非定常自由面流体流动，它们有着各自的优势和局限性。

采用固定网格的欧拉描述，整个计算过程中计算网格始终保持初始状态，从而可以描述流体质点运动的急剧变化，如碎波等现象。

欧拉描述虽然可以有效地分析整个流场内部的运动，但很难精确跟踪流体的自由液面，即很难给出准确的自由面形状和位置。

在拉格朗日描述中，网格结点与流体质点在整个运动过程中始终保持重合，流体质点与网格结点之间不存在相对运动，因此很容易跟踪自由液面，适用于线性小晃动问题。

这不仅大大地简化了控制方程地求解，而且还能有效地跟踪流体质点的运动轨迹，准确地描述波动的自由液面。

但是，在涉及求解带自由面流体大幅运动时，此时的晃动已经具有很强的非线性特征，如果还采用拉格朗日描述，由于流体质点运动的急剧变化，将导致计算网格的扭曲，会面临网格奇异问题，从而使计算无法继续进行。

拉格朗日描述和欧拉描述虽有各自的优点，但也存在较大的缺陷，如果将它们有机地结合在一起，充分利用各自的优点并克服其缺点，则可以解决各自都难于解决的问题，任意拉格朗日–欧拉描述[6－7]（ALE）方法就是基于该思路提出的。

在任意拉格朗日–欧拉描述中，网格结点的运动方式比较灵活，网格结点可以跟随流体质点一起运动，也可以固定不变，甚至可以采用网格结点在一个方向上固定而在其他方向上随流体质点一起运动等方式。

为了更加清晰地理解这三种描述方法，本研究从以下几个方面进行阐述和比较。

- 164 -暨两岸船舶与海洋工程水动力学研讨会文集- 165 -2 坐标系如图1所示，将连续介质在初始t 0时刻的形状称为初始构形，记为XΩv ，并引入拉格朗日坐标系（运动坐标系）123OX X X ，则在该参考系下物质点P 的位置矢量可表示为123(),,X X X X =v，坐标(123),,i X i =称为物质坐标，表示物质点P 在物质坐标系中的初始位置。

将连续介质在现时t 时刻的形状称为现时构形，记为x Ωv ，并引入欧拉坐标系（空间坐标系）ox x x 123，则在该参考系下，物质点P 位于p 处，其位置矢量可表示为),,(321x x x x =ϖ，坐标)3,2,1(=i x i 称为空间坐标，表示物质点P 在空间坐标系中的现时位置。

最后任意选择一个独立于初始构形和现时构形的任意构形，记为χΩv ，并引入任意拉格朗日-欧拉坐标系（任意坐标系）o χχχ123，则在该参考系下，物质点P 位于q 处，其位置矢量可表示为),,(321χχχχ=ϖ，坐标123(),,i i χ=称为任意坐标，表示物质点P 在任意坐标系中的现时位置。

图1 参考坐标系拉格朗日描述（简称为L 描述）是在初始构形XΩv 的基础上来研究物质点X v在空间坐标系中的运动规律，其数学表达为(),X t x x =v v v ，描述了同一物质点X v 在空间坐标系中的位置矢量随时间的变化情况。

欧拉描述（简称为E 描述）是在现时构形x Ωv 的基础上来研究空间点x ϖ处物质点的运动规律，其数学表达为(),X X x t =v v v ，描述了同一空间点x ϖ处物质点的各物理量随时间的变化情况。

任意拉格朗日-欧拉描述（简称为ALE 描述）则是在任意构形χΩv 的基础上来研究任意点χϖ在暨两岸船舶与海洋工程水动力学研讨会文集- 166 -空间坐标系中的运动规律，其数学表达为()t x x ,χϖϖϖ=，描述了同一任意点χϖ在空间坐标系中的位置矢量随时间的变化情况。

而另一数学表达式为(),X t χχ=v v v ，描述了同一物质点X v 在任意坐标系中的位置矢量随时间的变化情况。

3 物质导数：速度和加速度某一物质点X ϖ在空间坐标系中的运动速度u ϖ等于该物质点在空间坐标系中的位置矢量),(t X x x ϖϖϖ=对时间t 的导数X t t X x u ϖϖϖϖ∂∂=/),(，式中的X ϖ表示物质点坐标X ϖ固定时的值。

某一物质点X ϖ在任意坐标系中的运动速度v ϖ等于该物质点在任意坐标系中的位置矢量),(t X ϖϖϖχχ=对时间t 的导数X t t X v ϖϖϖϖ∂∂=/),(χ，式中的Xϖ表示物质点坐标X ϖ固定时的值。

某一任意点χϖ在空间坐标系中的运动速度w ϖ等于该任意点在空间坐标系中的位置矢量),(t x x χϖϖϖ=对时间t 的导数χχϖϖϖϖt t x w ∂∂=/),(，式中的χϖ表示任意点坐标χϖ固定时的值。

在物质点描述（L 描述）方法中，有限单元剖分是对整个运动流域进行剖分的，其网格点是随同流体质点一起运动的，网格点就是物质点。

在空间点描述（E 描述）方法中，有限单元剖分是对整个固定流域进行剖分的，其网格点是固定在空间中不动的，网格点就是空间点。

在任意点描述（ALE 描述）方法中，有限单元剖分是对整个任意流域进行剖分的，其网格点是独立于流体质点和空间点任意运动的，可以根据需要自由选择运动方式，网格点就是任意点。

在ALE 描述中，由于网格点的运动规律可以是任意给定的，那么指定网格点特殊的运动规律就可以将ALE 描述退化为L 描述和E 描述。

u w ϖϖ=，即网格点是随同物质点一起运动，此时退化为L 描述。

0=w ϖ，即网格点是固定在空间中不动，此时退化为E 描述。

0≠≠u w ϖϖ，即网格点是独立于流体质点和空间点随同任意点一起运动，此时对应于一般的ALE 描述。

各物理量的实质导数等于物理量的网格导数与物理量的迁移导数之和，实质导数为某一固定物质点X ϖ的物理量对时间的变化率；网格导数为某一固定物质点X ϖ在网格点处的物理量对时间的变化率；迁移导数为某一固定物质点X ϖ相对于网格点的速度（对流速度）和物理量对空间的变化率（迁移加速度）的乘积。

在ALE 描述中任意构形是已知的，所以各物理量用物质点在任意坐标系中的位置矢量),(t X ϖϖϖχχ=来表示会比较方便，即用)),,((t t X f f ϖϖχ=表示各物理量。

则f 对t 的全导数（物暨两岸船舶与海洋工程水动力学研讨会文集- 167 -质导数）为/////X Df Dt f tf t f t χ∂∂∂∂∂∂χ∂χ∂==+v v v vg ，X ϖ表示物质点坐标X ϖ固定时的值，χϖ表示任意点坐标χϖ固定时的值。

为了研究方便，下面引入求和法则，各物理量可用)),,((t t X f f j ϖχ=来表示，则可得/////j j j X Df Dt f t f t f t χ∂∂∂∂∂∂χ∂χ∂==+v g 。

那么对空间坐标系中的位置矢量)),,((t t X x j i ϖχ求全导数后(即为物质点X ϖ在空间坐标系中的运动速度i u )可得()(,),////ji i j i i j j X u x X t t t x t x t χχχχ=∂∂=∂∂+∂∂∂∂v v g ，整理后得//i i i j j u w x t χχ−=∂∂∂∂g 。

令i i ic w u =−，则/i i j j c x v χ=∂∂g ，i c 为物质点相对于网格点的运动速度，称为对流速度。

由式/i i j j c x v χ=∂∂g 可得//j i j i t c x χχ∂∂=∂∂g ，通过整理后可得物质点Xϖ的任意物理量的物质导数///j i i X f t f t c f x χ∂∂=∂∂+∂∂v g 。

那么物质点X ϖ在空间坐标系中运动速度)),,((t t X u j i ϖχ的全导数（即为物质点X ϖ在空间坐标系中的加速度）为//j i i k i k a u t c u x χ=∂∂+∂∂g ，其中j t u i χ∂∂/为网格导数，等于空间坐标系下物质点的运动速度i u 在网格点坐标j χ固定时对时间的偏导数；k c 为对流速度，等于空间坐标系下物质点的运动速度k u 与网格点的运动速度k w 之差。

在L 描述下（以运动坐标系为参考的，网格点就是物质点），网格导数j t u i χ∂∂/就变为X i t u ϖ∂∂/，等于空间坐标系下物质点的运动速度i u 对时间的全导数；对流速度kc 就为0，网格点的运动速度k k u w =，有Xi i t u a ϖ∂∂=/。

在E 描述下（以空间坐标系为参考的，网格点就是空间点），网格导数j t u i χ∂∂/就变为j x i tu ∂∂/，等于空间坐标系下物质点的运动速度i u 对时间的偏导数；对流速度k c 就为k u ，网格点的运动速度0=k w ，有//j i i k i k x a u t u u x =∂∂+∂∂g 。

4 有限元分析在进行流体的有限元分析时，如何对流体的不可压条件进行处理是非常关键的问题。

处理的方法有很多种，目前主要有泊松方程法[8]和分步法[9-10] 。

前者的缺点是只能采用速度和压力的不等阶插值，大大增加了计算量。

而后者则利用不可压条件对NS 方程采用分步求解，将压力和速度变量分开后再进行方程的有限元离散，对速度和压力采用等阶插值，这样不仅实现方便，而且计算精度高。

采用有限元法来建立粘性流体运动问题的有限元方程组的步骤如下：1）选择描述方法，导出在该描述下的流体运动控制方程，包括连续性方程和运动方程。

2）采用直接法（直接从运动方程获得压力-泊松方程）或中间速度法[11]（其中中间速度法又暨两岸船舶与海洋工程水动力学研讨会文集可分为两种：在NS方程中改变压力项来获得压力-泊松方程或在NS方程中略去压力项来获得压力-泊松方程）对运动方程进行时域上的差分离散和空间域上的有限元离散，建立不可压黏性流体运动问题的有限元方程组。

3）针对建立的有限元方程组，对边界条件（包括本质边界条件和自然边界条件）进行适当处理。

用ALE有限元法求解带自由面流动问题时，无需直接求解自由面运动方程，只要通过自由面上网格结点的运动来描述自由面变化。