第四章 求解水动力弥散方程的有限单元法2

第四章 求解水动力弥散方程的有限单元法里兹法有限元和加权余量法有限元都能有效地用于求解地下水运动问题。

但里兹法仅适用于对称正定算子，因此用里兹法求解水动力弥散方程时，寻求相应的泛函比较困难，常要求对方程进行适当交换，去掉一阶导数数项，使变换后的方程为自伴的，用有限元求得的方程组对称正定，易于求解。

但是，这种变换多为指数变换，当Peclet 数较大时，变换会引入较大误差，有时由于指数过大，数值处理会产生溢出。

加权余量法可以不加限制地将一般的微分方程通过方程余量和权函数正交化途径，化为代数方程组而获得近似解。

由于它不要求寻找与微分方程相应的泛函，无论对自伴、非自伴或非线性问题都可以求解。

因此，加权余量法，尤其迦辽金法有限元在求解水动力弥散方程中应用较广。

有关有限单元法的一些基本概念，已在前两章节中做了介绍，本章将主要介绍迦辽金有限单无法在求解水动力弥散方程中的应用求解对水动力弥散方程的几种特殊解法。

第一节 用迦辽金有限元求解水质模型的基本方法如前所述，为求得计算区域溶质浓度的变化，应首先建立相应的数学模型，并给出相应的初始和边界条件。

由第一章可知，对非均质流体，该数学模型由以下方程组成：()5713))((--∂∂+∂∂-=j j j s iji x z g x p k v ρμθ重力加速度g j 是矢量，所以g j ；只在x j =z 对取值。

()5813),(),,(--==p c p c μμρρ对低浓度不可压缩流体，可用其一阶近似式；()()6013)()(5913)(),(0000--∂∂-∂∂∂∂=∂∂---+=-+=i ii ij i cv x x c D x t c c c c c βμμαρρ或()()6213)])([(6113)())(()(1--∂∂+∂∂∂∂=∂∂--∂∂-∂∂∂∂=∂∂jj j ij i i i j c ij i c x z g x p k x t p C v x x D x t ρμρμθθθ对均质流体该模型由以下方程组成：()()()63136513))((6413))((--∂∂-=--∂∂∂∂=∂∂--∂∂∂∂=∂∂ii ii di i s x h k v x h kh x t h x h KM x t h θμμ潜水承压水对流弥散方程同式（3－l －60）或式（3－l －61）。

式中各符号意义同第一章说明。

在均质流体中，由于ρ和μ为常数，浓度分布不影响流速分布，水流问题和浓度问题可以分别独立求解。

即首先由水流方程求得流速分布，再将流速作为已知量，代人对流一弥散方程求得浓度分布。

而在非均质流体情况下，由于ρ和μ与浓度c 有关，浓度变化亦影响流速，因此水流问题和溶质运移问题必须联立求解，其中最常用的是迭代求解法。

在这种求解过程中，首先给出密度ρ和粘度μ，由运动方程求出流速公布，再求出浓度分布，然后代入状态方程求出ρ、μ，与初值比较，进一步修正ρ、μ值。

在一个计算时段内重复上述计算过程，直到收敛为止。

由上述介绍可见，无论对均质还是非均质流体，都须分别计算流速分布和浓度分布。

以下具体介绍迦辽金有限元解法。

一、压力p 的计算根据迦辽金有限元法的概念，将压力p 的试函数p ~通过基函数的线性组合来表示，即 ()()143,...,2,1--=Φ=≈N J p p p j j式中，J Φ为终点J 的基函数。

将水流连续方程（3－l －62）写成迦辽金积分表达式，即()2430}))(({1--=ΩΦ∂∂-⎥⎥⎦⎤⎢⎢⎣⎡+∂∂∂∂⎰Ωd t pg x p k x L i j ij jμρμρ上式二阶导数应用分部积分和格林公式，并将积分转化到各个单元上进行，可得到下式：()()343,....,3,2,10~~~21--==Φ⎪⎪⎭⎫ ⎝⎛+∂∂⎪⎪⎭⎫ ⎝⎛+Ω⎥⎦⎤Φ∂∂-∂Φ∂⎢⎢⎣⎡⎪⎭⎫ ⎝⎛+∂∂⎪⎪⎭⎫ ⎝⎛-∑∑⎰ΓΩN I dL n g x p k d t p x g x p k L i i j ij i L L i j ij ρμρμρμρ式中，n i 为边界民的外法向分量；L 为单元e 中二类边界长度。

将压力p 的试函数J J p p Φ=~代入立式得：()443~1--Φ⎪⎪⎭⎫ ⎝⎛+∂∂-Ω∂Φ∂+Ω∂Φ∂∂Φ∂+ΩΦΦ∂∂∑⎰∑⎰∑⎰∑⎰ΓΩΓΩΩdL n g x p k d x g k d p x x k d t p I i j j ij Ij ij J i J j J ij I J Jρμρρμρμρμ在相应定解条件下求解上式，可得到区域Ω内各节点的压力分布P J ：二、流速分布的计算流速计算有两种方法。

第一种方法是将求得的压力p j 代入式（3－4－1），再代人运动方程（3－l －57）即得：()543--⎪⎪⎭⎫ ⎝⎛+∂Φ∂⎪⎪⎭⎫ ⎝⎛-=j j jJ iiji g x p k v ρμθ由上式求得的流速在单元内部是常数，但在单元与单元之间交界处不连续，这对用某些方法计算浓度会有一定影响。

第二种方法是对运动方程（3－1－57）进行迦辽金有限元求解。

首先将流速用基函数的()643~--Φ==Jij i i v v v线性组合来表示，即然后将运动方程写成迦辽金积分表达式：⎰Ω==ΩΦ+∂∂+Φ) ..., ,2 ,1( 0)]~)(([n i g g x pk v I i j i ijJ ij ρμθ （3-4-7）把式（3－4－l ）代入式（3－4－7），并在各单元上进行积分即得到计算各结点流速的有限元方程：0=ΩΦ⎪⎪⎭⎫ ⎝⎛+∂Φ∂+ΩΦΦ∑⎰∑⎰ΓΩΓΩd g x p k d v I j j JJ i ij I J ij ρμθ 三、水动力弥散方程的求解为进行有限元计算：：同样将浓度C 写成基函数的线性组合，即()n J C C C J J ,...,2,1~=Φ=≈将对流一弥散方程写成迦辽金积分表达式；()()()1043,...,2,10~~~--==ΩΦ⎥⎥⎦⎤⎢⎢⎣⎡∂∂-∂∂-⎪⎪⎭⎫ ⎝⎛∂∂∂∂⎰Ωn I d t c c v x x C D x I i i j ij i对上式二阶导数应用分部积分和格林公式，则有()1143~~~--Φ∂∂+Ω∂Φ∂∂∂-=ΩΦ⎪⎪⎭⎫ ⎝⎛∂∂∂∂⎰⎰⎰ΓΓΩdL n x c D d x x c D d x c D x I i j ij i I j ij I j ij i式中， Γ为区域Ω的二类或三类边界段； n i 为边界Γ的外法向矢量。

将武（3－4－－11）代入式（3－4～10），并将积分转为在单元上进行，则得到；()()12430~~~~--=Φ∂∂+Ω⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧∂Φ∂∂∂-Φ⎥⎦⎤⎢⎣⎡∂∂-∂∂-∑⎰∑⎰ΓΓΓΩdL n x c D d x x c D t c C v x I i j ij i I j ij I i i上式左侧最后一项面积积分表示弥散通量的边界条件。

Γ表示单元e 与二类边界Γ相交的部分，若两者不相交，则该线积分为零。

对于包含有三类边界，即已知对流一弥散通量的边界，为了边界条件处理方便，可以对弥散方程中的二阶和一阶导教体应用格林公式，即有：()()1343~~~~--Φ⎪⎪⎭⎫ ⎝⎛-∂∂+Ω∂Φ∂-=ΩΦ⎥⎥⎦⎤⎢⎢⎣⎡∂∂-⎪⎪⎭⎫ ⎝⎛∂∂∂∂⎰⎰⎰ΓΩΩdL n C v x c D d x D d C v x x c D x I i i jij j I ij I i i j ij i此式代入式（3－4－10），并在单元中进行积分，则得：()14430~~~2--=Φ⎪⎪⎭⎫ ⎝⎛-∂∂+Ω⎥⎥⎦⎤⎢⎢⎣⎡∂Φ∂⎪⎪⎭⎫ ⎝⎛-∂∂-Φ∂∂-∑⎰∑⎰ΓΩdL n C v x c D d x C v x c D t c I i i j ij i I i j ij I 将近似函数J J c C Φ=~代人式（3－4－12）和式（3－4－14），分别得到二类和三类边界的迦辽金有限元积分或为二类边界：()()15430~--=∂∂-Ω∂∂+Ω⎥⎥⎦⎤⎢⎢⎣⎡∂∂∂∂+∂∂∑⎰∑⎰∑⎰ΩΩΓe e e I i j j J J I J i I j J ij I J i i e e dL n x cDi d t C d C x x D v x φφφφφφφ三类边界：()()16430~--=∂∂-Ω∂∂+Ω⎥⎥⎦⎤⎢⎢⎣⎡∂∂∂∂+∂∂∑⎰∑⎰∑⎰ΓΩΩdL n x cD d t c d C x x D v x I i j ij j J I j i I j J ij i j i i φφφφφφφ式（3－4－16）所表示的处理方法对流速的求解比较敏感。

三类边界问题还可以用以下处理办法。

如三类边界条件表示为()1743)(''3--=-∂∂ΓC n v n C v x cD i i i i jij其中C’’表示自边界外侧区域进入边界处的液体浓度。

由此可得：''C n v n C v n x cD i i J i j i i jij+=∂∂φ 此式代人式（3－4－16），整理后得到三类边界的另一种计算式：()()18430''3--=-Ω∂∂+-Ω⎥⎥⎦⎤⎢⎢⎣⎡∂∂∂∂+∂∂∑⎰∑⎰∑⎰∑⎰ΩΓΓΩεεεφφφφφφφφφdL C n v d t C dL C n v d C x x D v x I i i JJ I J J I i i J j J i I ij I J i i 上述式（3－4－4）、式（3－4－5）［或式（3－4－8）］、式（3－4－15）及式（3－4－16）或式（3－4－18）即是对水质模型进行迦辽金有限元计算的基本公式。

在具体应用中，由于剖分单元的形状和单元基函数的类型不同，由上述基本算式可以得到相应的有限元方程。

关于水流方程的有限元计算问题，已在第一篇第五章作了介绍，以下仅就二维对流一弥散方程的有限元计算方法作一说明。

设剖分单元为三角形，单元编号用e 表示，单元中三个结点编号分别以i 、j 、m 表示。

因是二维问题，j i x x ,的下标i 、j=1,2，相应x 1＝x ，x 2=y ，由此式（3-4-12）变为()()()()1943,....,2,10~~~~~~~~~~--==Γ⎥⎥⎦⎤⎢⎢⎣⎡⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+⎭⎬⎫⎩⎨⎧⎪⎪⎭⎫ ⎝⎛∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂∂∂-⎥⎦⎤⎢⎣⎡∂∂-∂∂-∂∂-∑⎰∑⎰Γ∆n i d n y cD x c D n y c D x c D y y c D y x c D x y c D x x c D t c C v y C v x i y yy yx x xy xx t iyy i yx i xy i xx i y x φφφφφφε为了进行有限单元的计算，需构造三角形单元的基函数，其中最简单的是线性插值。