ABSTRACT

The Galerkin finite element method for the diffusion of water contamination is analyzed. The time variable is discretized by backward Euler method. Finally, we did numerical simulations with a concrete example.

Keywords:Water Contamination Diffusion Equation, Galerkin Finite Element Method, Backward Euler Method

水污染问题的向后Euler-Galerkin 方法及数值模拟

娜扎开提·阿迪力，伊马木·麦麦提，阿布都热西提·阿布都外力^*

新疆大学数学与系统科学学院，新疆 乌鲁木齐

收稿日期：2017年11月29日；录用日期：2017年12月15日；发布日期：2017年12月22日

摘 要

分析了浅水流动中的水污染扩散问题的Galerkin有限元方法，对于时间变量用向后Euler法进行离散。最后通过具体例子进行了数值模拟。

关键词 :水污染扩散方程，Galerkin有限元方法，向后Euler法

Copyright © 2017 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

浅水流中的污染物扩散问题有比较重要的实用价值。在相对较浅的自由表面流动中，水平速度为主要分量，问题可以适当的近似为二维问题。Navier-Stokes方程是流体力学中的基本方程，因此可当该方程为出发点，应用沿深度积分的方法可得到浅水流方程 [1] 。浅水流作为污染物运输的载体，可以写出一种污染物扩散问题的方程 [1] ：

$\frac{\partial \left(HC\right)}{\partial t}+\frac{\partial \left(H{U}_{1}C\right)}{\partial x}+\frac{\partial \left(H{U}_{2}C\right)}{\partial y}=\frac{\partial }{\partial x}\left(H{k}_1\frac{\partial C}{\partial x}\right)+\frac{\partial }{\partial y}\left(H{k}_2\frac{\partial C}{\partial y}\right)+Q,(x,y)\in D,t\in I=(0,T).$ (1)

其中C为污染物浓度；Q为源项； $k_1$ 、 $k_2$ 分别为x、y方向的扩散系数；H为水深； $U_1$ 、 $U_2$ 分别为流速在x、y方向的分量。水深与流速可由浅水流方程确定，本文的数值例子中将它们预先给出。

利用深度平均的连续守恒方程

$\frac{\partial H}{\partial t}+\frac{\partial \left(H{U}_1\right)}{\partial x}+\frac{\partial \left(H{U}_2\right)}{\partial y}=0,$

将(1)式进一步改写成

$H\frac{\partial C}{\partial t}+HU\cdot \nabla C-\nabla \cdot (HK\nabla C)+Q=0,(x,y)\in D,t\in I.$

其中 $U={\left(U_1,U_2\right)}^{\text{T}}$ ， $K=diag\left(k_1,k_2\right)$ 。

王焕焕和冯新龙 [2] 分析了水污染扩散问题的两种有限差分格式，证明了其相容性、稳定性与收敛性。王焕 [3] 研究了该问题的特征有限元方法，得到了 $H^1$ 误差和 $L_2$ 误差估计。

为简便，本文中以图1所示的河的示意图为例，假设河内没有污染源。令扩散系数 $k_1=k_2=k$ (常数)。假设 $C_2$ 与边界中的污染源充分远，可认为污染物浓度在研究的时间段内为零， $C_1$ 是河岸。因此问题可写成

$\{\begin{array}{l}H\frac{\partial C}{\partial t}+HU\cdot \nabla C-\nabla \cdot \left(Hk\nabla C\right)=0,\left(x,y\right)\in D,t\in I,(2)\\ C\left(x,y,0\right)=C_0\left(x,y\right),\left(x,y\right)\in D,(3)\\ {\frac{\partial }{\partial n}C\left(x,y,t\right)|}_{C_1}=0,t\in I,(4)\\ {C\left(x,y,t\right)|}_{C_2}=0，t\in I.(5)\end{array}$

2. 向后Euler-Galerkin有限元形式

本文使用下列空间，对于非负整数 $m$ 和实数 $k\in \left[1,\infty \right]$ ，定义

$W^{m,r}\left(D\right)=\left\{u\in L^r\left(D\right)|{\partial }^{\alpha }u\in L^r\left(D\right),\left|\alpha \right|\le m\right\}.$ [4]

图1. 河的形状示意图

在此空间上定义范数如下：

${\Vert u\Vert }_{m,r,D}={\left({\displaystyle \sum _{\left|\alpha \right|\le m}{\left|{\partial }^{\alpha }u\right|}^r\text{d}x}\right)}^{1/r},1\le r<\infty ,$

当 $r=2$ 时，记 $W^{m,r}\left(D\right)$ 为 $H^m$ ，相应的范数为 ${\Vert \cdot \Vert }_m$ 。

现在推导(2)~(5)的变分形式。记 $H_0^1=\left\{v\in H^1;{v|}_{C_2}=0\right\}$ ，用函数 $v\left(x,y\right)\in H_0^1$ 与方程(2)的两端做内积

$\left(H\frac{\partial C}{\partial t},v\right)+\left(HU\cdot \nabla C,v\right)-\left(\nabla \cdot \left(Hk\nabla C\right),v\right)=0,t\in I.$

因此(2)~(5)的变分形式为：求 $C\left(t\right):I\to H_0^1$ ，使得

$\{\begin{array}{l}\left(H\frac{\partial C}{\partial t},v\right)+A\left(C,v\right)=0,\forall v\in H_0^1,t\in I,\\ C\left(0\right)=C_0\left(x,y\right).\end{array}$ (6)

其中 $A\left(\cdot ,\cdot \right)$ 为下式所定义的双线性形式

$A\left(u,v\right)={\displaystyle {\int }_D\left[-v\nabla \cdot \left(Hk\nabla C\right)+HvU\cdot \nabla C\right]\text{d}x\text{d}y},\forall u,v\in H_0^1.$

$A\left(\cdot ,\cdot \right)$ 利用Green公式及(4)、(5)式可简化成

$A\left(C,v\right)={\displaystyle {\int }_D\left[Hk\nabla C\cdot \nabla v+HvU\cdot \nabla C\right]\text{d}x\text{d}y}.$

设 $T_h=\left\{\tau \right\}$ 为 $\overline{D}$ 上的拟均匀剖分，网格参数 $h\le h_0$ ，选取的有限元空间为

$V_h=\left\{v|v\in H_0^1\left(D\right),{v|}_{\tau }\in P_r\left(\tau \right),\forall \tau \in T_h\right\}.$

因此对于问题(2)~(5)的Galerkin有限元形式是：求 $C_h\left(t\right):I\to V_h$ ，使得

$\{\begin{array}{l}\left(H\frac{\partial C_h}{\partial t},v_h\right)+A\left(C_h,v_h\right)=0,\forall v_h\in V_h,t\in I,\\ C_h\left(0\right)=C_{0h}.\end{array}$ (7)

其中 $C_{0h}\in V_h$ 是 $C_0$ 的某个近似。

现设 $\left\{{\varphi }_1\left(x,y\right),{\varphi }_2\left(x,y\right),\cdots ,{\varphi }_N\left(x,y\right)\right\}$ 为空间 $V_h$ 的基底，将未知函数 $C_h$ 表示为

$C_h\left(x,y,t\right)={\displaystyle \sum _{i=1}^N{c}_i\left(t\right){\varphi }_i\left(x,y\right)},\left(x,y\right)\in D.$

并在(7)中取 $v_h={\varphi }_j\left(x,y\right),j=1,2,\cdots ,N$ ，则有

$\{\begin{array}{l}{\displaystyle \sum _{i=1}^N\frac{\text{d}{c}_i\left(t\right)}{\text{d}t}\left(H{\varphi }_i,{\varphi }_j\right)}+{\displaystyle \sum _{i=1}^N{c}_i\left(t\right)\left[\left(Hk\nabla {\varphi }_i,\nabla {\varphi }_j\right)+\left(HU\cdot \nabla {\varphi }_i,{\varphi }_j\right)\right]}=0,j=1,2,\cdots N,t\in I,\\ c_i\left(0\right)=c_i^0,\end{array}$ (8)

其中 $c_i^0$ 为 $C_{0h}$ 按基函数 $\left\{{\varphi }_i\left(x,y\right)\right\}$ 展开的系数。

引进矩阵和向量记号： $M={\left(m_{ij}\right)}_{N\times N}$ ， $A={\left(a_{ij}\right)}_{N\times N}$ ， $c={\left(c_1,c_2,\cdots ,c_N\right)}^{\text{T}}$ ，其中 $m_{ij}=\left(H{\varphi }_i,{\varphi }_j\right)={\displaystyle {\int }_{D}H{\varphi }_i{\varphi }_j\text{d}x\text{d}y}$ ， $a_{ij}=\left(Hk\nabla {\varphi }_i,\nabla {\varphi }_j\right)+\left(HU\cdot \nabla {\varphi }_i,{\varphi }_j\right)$ 。记M为质量矩阵，A为刚度矩阵。借助上述记号，可将(8)改写成

$\{\begin{array}{l}M\frac{\text{d}c}{\text{d}t}+Ac=0\begin{array}{c},\begin{array}{c}t\in I,\end{array}\end{array}\\ c\left(0\right)=c^0={\left(c_1^0,c_2^0,\cdots ,c_N^0\right)}^{\text{T}}.\end{array}$

如果再对时间变量用有限元离散会使计算量有成量级的增加。因此在时间 $t=t_n$ 上我们采用向后Euler差商，从而得到的全离散格式为

$\{\begin{array}{l}M\left(\frac{c^n-c^{n-1}}{\Delta t}\right)+A{c}^n=0,\\ c^0={\left(c_1^0,c_2^0,\cdots ,c_N^0\right)}^{\text{T}}.\end{array}$ (9)

其中 $\Delta t$ 为时间步长。

根据文献 [4] 可知，当有限元形式(7)中的初值 $C_{0h}$ 满足条件

$\Vert C_{0h}-C_0\Vert \le \alpha h^r{\Vert C_0\Vert }_r\left(\alpha 为常数\right)$

时，相应的近似解 $C_h\left(t\right)$ 与变分形式(6)的真解 $C\left(t\right)$ 之间的误差估计为

$\Vert C_h\left(t\right)-C\left(t\right)\Vert \le \alpha h^r\left\{\Vert C_0\Vert +{\displaystyle {\int }_0^t\Vert C_t\left(s\right)\Vert \text{d}s}\right\}.$

易知全离散格式(9)对时间t是一阶精度的。

3. 数值实验

针对图1所示的河，考虑初边值问题(2)~(5)的数值计算。设 $U\approx \left(0\text{m/s,}0\text{m/s}\right)$ 、k = 15 m²/s，初始条件为

$C_0\left(x,y\right)=\{\begin{array}{l}200\text{g}/{\text{m}}^{\text{3}},x=0,y=0,\\ 0,其他\end{array}$

为计算简便，将水深函数 $H\left(x,y\right)\left(0\le y\le 300,单位:\text{m}\right)$ 取为

$H\left(x,y\right)=\{\begin{array}{l}2,0\le x\le 200;\\ 4,200\le x\le 400;\\ 6,400\le x\le 600;\\ 8,600\le x\le 800;\end{array}$

x轴方向20等分，y轴方向15等分便得到矩形单元，这里有限元空间取为

$V_h=\left\{v|v\in H_0^1\left(\Omega \right),{v|}_{\tau }\in P_1\left(\tau \right),\forall \tau \in T_h\right\}.$

图2. 河中污染物浓度随时间变化的分布情况

取时间步长 $\Delta t=0.1\text{s}$ ，用全离散格式(9)进行计算，就可以得到数值解。图2(a)~图2(d)给出了当时间

t =100 s；t = 500 s；t = 1000 s；t = 2000 s时河中污染物浓度分布情况，图中的高度表示浓度大小。

4. 结论

本文中推导出了水污染扩散问题的向后Euler-Galerkin方法。通过数值模拟的结果看到效果是较好的。

文章引用

娜扎开提·阿迪力,伊马木·麦麦提,阿布都热西提·阿布都外力. 水污染问题的向后Euler-Galerkin方法及数值模拟
Numerical Analysis of Backward Euler-Galerkin Method for Water Contamination Problem[J]. 应用数学进展, 2017, 06(09): 1146-1150. http://dx.doi.org/10.12677/AAM.2017.69139

参考文献 (References)