多孔课程设计-206310111-李睿.doc

多孔介质污染物迁移动力学课程设计报告姓 名：李 睿学 号：2006310111班 级：环 研 三指导教师：王 洪 涛2007年1月28目 录一 数学模型11.1 假设条件11.2 地下水运动数学模型11.3 污染物迁移数学模型2二 有限元法32.1 剖分研究域32.2 构造基函数32.3 形成有限元方程6三算法实现103.1 剖分研究域103.2 存储原始数据103.3 生成系数矩阵103.4 处理特殊节点113.5 程序设计流程113.6 源程序11四 程序验证124.1 地下水运动数学模型程序验证124.2 污染物迁移数学模型程序验证14五 自拟问题165.1 问题描述165.2 问题求解17六 建议小结19附录20一 数学模型1.1 假设条件1) 含水介质与参数：含水层是非均质各向异性承压含水层，水运动和污染物迁移参数已知；2) 水流和污染物迁移条件：地下水的运动和污染物的迁移是平面二维的，是非稳定的，含水层中仅有一种可溶污染物；3) 源汇条件：研究域中存在源汇的作用，包括点源(抽水井和注水井)、面源等作用，污染物迁移考虑吸附作用；4) 初始条件：在初始时刻，全研究域上的水头分布和污染物浓度分布已知；5) 边界条件：研究域的边界可以是第一、二、三类（水流问题仅有一、二类）边界或其组合，也可以有内边界，对于给定边界类型，边界条件已知。1.2 地下水运动数学模型在上述假设条件下，地下水二维渗流可以用下面数学模型描述：微分方程 （1.1）初始条件 （1.2）第一类边界条件 （1.3）第二类边界条件 （1.4）式中，：待求水头，L；S：贮水系数，无量纲；Txx，Tyy：导水系数分量，L2T-1；WM：源汇项，LT-1，为单位时间单位面积含水层得到的源汇水体积；H0(x,y)：给定研究域上的初始水头，L；H1(x,y,t)：第一类边界G1上给定的水头，L；q(x,y,t)：在第二类边界G2上给定的水通量，L2T-1，为单位时间垂直通过单位边界长度进入研究域的水体积；G：二维平面研究域；G：研究域的边界，；cos(n,x)，cos(n,y)：边界外法线向量与坐标轴正向之间夹角的余弦。1.3 污染物迁移数学模型非稳定二维对流弥散可以用下面数学模型描述：微分方程 （1.5）式中，L为算子 （1.6）为了书写简便起见，记，则有 （1.7）式中， （1.8）初始条件 （1.9）第一类边界条件 （1.10）第二类边界条件 （1.11）第三类边界条件 （1.12）式中，：浓度，ML-3；Dij：水动力弥散系数张量的坐标分量，L2T-1；IM：源汇项，ML-2T-1，为单位时间、单位面积含水层得到的污染物质量，污染物进入研究域时为正，流出研究域时为负；We：单位时间单位面积含水层得到的水量，LT-1；Ce：We所含污染物的浓度，ML-3；Wo：单位时间单位面积含水层的开采水量（开采为负），LT-1，开采水的浓度为C；Rd：滞留因子，无量纲；l：一级化学反应常数，T-1；n：空隙度；ux、uy：x和y方向上的水流实际速度，LT-1；C0(x,y)：给定研究域上的初始浓度，ML3；C1(x,y,t)：第一类边界G1上给定的浓度，ML-3；f(x,y,t)：在第二类边界G2上给定的弥散通量，ML-1T-1，为由于弥散作用在单位时间垂直通过单位边界长度进入研究域的污染物质量，进入研究域为正，离开为负；g(x,y,t)：在第三类边界G3上给定的污染物通量，ML-1T-1，为单位时间垂直通过单位边界长度进入研究域的污染物质量，流入研究域为正，流出为负；b：含水层厚度，L；对于承压水问题，；对于潜水（无压水）问题，h为潜水含水层厚度，H为水头，z为含水层底板高度；G：研究域的边界，；其余符号同前。由于污染物迁移模型通过简化就可以得到渗流模型，所以我们讨论污染物迁移模型的有限差分方程，然后通过简化得到渗流模型的有限差分解。二 有限元法拟用有限单元法求解上述地下水运动和污染物迁移数学模型。不难看出，上述两个数学模型是同时存在的。为了简化问题，通常的方法是先求地下水运动数学模型的解H，由此，与水流有关的参量如ux、uy、b就是已知的。比较可知，水运动方程（1.1）是对流弥散方程（1.5）在忽略对流项条件下的一个特例。因此，我们将从污染物迁移模型出发，讨论其有限元数值解。同时，给出水流模型的解，但并不重复推导过程。2.1 剖分研究域剖分研究域G成三角形网，共剖分为Ne个三角形单元，表示为：，第b三角形的面积用表示。共剖分为Np个节点，表示为：。其中已知浓度的节点为m个。 如果能求出节点上的浓度，那么研究域上任一点的浓度都可以用节点浓度的线性组合表示出来，因此试探解为 （2.1）式中，为内插值函数，称为基函数；Ci(t)为节点i上的浓度。用节点浓度的线性组合来表示问题的解，这种有限元法称为线性元法。2.2 构造基函数线性元基本思想：是用一个平面来近似地代替这个曲面，且在单元的三个节点上两者重合，即节点浓度等于真实浓度。每个三角形都用一个小平面代替浓度曲面，则浓度曲面就可以用连在一起的小平面近似表示。 任选三角形单元，其三个节点（）的坐标是已知的，分别为、和；相应的浓度为、和，是待求的未知量。注意，单元的三个节点（）必须按逆时针排列，依次为p1、p2、p3。空间上通过三点的平面方程为对于内的任一给定点p(x,y)，其浓度C可由此方程计算出来。为了确定常数、，把（）三节点的坐标值和浓度值代入上面方程中，有联立求解方程组可以求得常数、，从而 （2.2）式中，为（p,p2,p3）三点构成的三角形的面积，见图2.1。 （2.3）为（p,p3,p1）三角形的面积；为（p,p1,p2）三角形的面积 （2.4） （2.5）为三角形的面积图2.1 单元面积构成 （2.6）其中，；。 显然， 与前述的方程（2.1）比较可知，若，则试探解应与上式相同。就是说，当时，p点的浓度仅与三角形三个节点的浓度有关，与其它节点无关。所以，（2.2）在上应展开为 （2.7）比较可知；而其余的（）。这样就得到了基函数在单元上的定义。图2.2 以i节点为定点的单元下面考虑任节点i的基函数，它是针对i节点定义的。假设以i为顶点的单元有个，如图2.2所示。若且计算点p位于三角形中，则。若p点位于以i为顶点的另一个三角形中，则，依此类推，显然的定义为 （2.8）由此得到了（）在研究域上的全部定义。基函数构造完毕。2.3 形成有限元方程用加权余量法形成有限元方程。由于试探解仅仅是近似解，所以将（2.7）代入微分方程中将产生误差，称为残差，用函数R(x, y,t)表示，即有 （2.9）式中，为试探解。残差的加权积分为零可以表示为 （2.10）式中，W(x,y)：权函数。Galerkin法：权函数与基函数相同，即 （2.11）将权函数代入上式可以得到m个方程即也即（为了书写方便期间，略去C上面的一杠） （2.12）引进记号逐个节点建立有限元方程，。若j =内部节点或第三类边界点，则对弥散项和对流项实施分部积分；若j第二类边界节点，则仅对弥散项实施分部积分。下面以j =内部节点或第三类边界点的情况为例加以讨论。得运用格林公式将二维面积分化为沿边界的线积分，得到考虑沿边界的积分，当j内部节点时，；当j三类边界节点时，且在第三类边界上，进入研究域的对流弥散通量已知，为g(x,y,t)。第三类边界条件为代入得微分方程的等效积分弱形式（2.13）把试探解代入，并交换积分与求和的顺序，注意Ni(x,y)与t无关；Ci(t)与（x,y）无关，得 （2.14）命（当j=第二类边界节点时，其余不变。）得 （2.15）初始条件变为 （2.16）方程（6.3.18a）中关于时间的导数用差分近似。将0t时间分成Nt个时间步，用k来计数，其中，tk1和tk为k+1和k时间步所对应的时间。在方程（2.15）中简单地用差分代替导数，并使用隐式差分格式，得整理得 （2.17）这就是二维问题的有限元方程。式中，：k+1时间步的待求浓度；：k时间步的已知浓度。当k=0时，用初始条件代入，其中为初始条件在i节点处的值。上面方程中，Mi,j、Ai,j和Fj等都是已知的系数，只有Ci是待求的节点浓度。对于给定k而言，有限元方程是线性方程组。有m个未知数（虽然方程左端有Np个浓度值，但只有m个是未知的），m个方程。求解之得节点上的浓度，问题得解。有限元方程可简写为 （2.18）式中，：质量矩阵；：弥散矩阵；：对流矩阵；：与浓度有关的源汇矩阵；：浓度列向量；：右端项列向量。下面讨论有限元方程组（2.17）系数的计算方法，这里直接给出结果，不作推导。 （2.19）式中， Ni,j为i和j共用三角形个数，为三角形的面积；、分别为三角形上含水层的孔隙度和滞留因子。 （2.20） （2.21）式中，、为单元上的水流实际速度，由水流场确定，是已知的，为 （2.22） （2.23）其中，N3j为以j为顶点的第三类（）边界边的个数，为第l个边界边的长度；为对应上的边界通量；为以j为顶点的三角形个数；为第单元的源汇强度。到此为止，我们得到了有限元方程组中各系数的计算公式。以上是关于污染物迁移的有限元方程，对于水运动的有限元方程，其形式上与之类似，具体变化为：质量矩阵M中的滞留因子Rd换成贮水系数S；弥散矩阵D中的水动力弥散系数D换成导水系数T；对流矩阵U和与浓度有关的源汇矩阵W为0；浓度列向量C换成水头列向量H；右端列向量F中的第三类边界换成第二类边界。三 算法实现3.1 剖分研究域采用Matlab自带的软件包pdetool进行研究域的剖分，并输出剖分后生成的Np2维节点坐标矩阵，Ne3维单元信息矩阵以及Nm7维边界信息矩阵。3.2 存储原始数据1) 存储单元信息用Np2维矩阵存储各个节点的坐标值，用Ne3维矩阵存储各单元三个节点的编号。2) 存储区域参数区域参数包括导水系数T，弹性贮水系数S，水动力弥散系数D，空隙率n，含水层厚度b，反应常数和滞留因子Rd，均采用Np行的矩阵（或向量）存储。其中T和D具有两个分量，因此用2列的矩阵存储，其他各参数均用一列向量存储。行维以单元编号为序。3) 存储边界条件用Matlab剖分生成的edges存储边界条件。4) 存储源汇项源汇项包括面源和点源。存在源汇的地方其值为源汇强度，不存在源汇的地方其值为零。3.3 生成系数矩阵由(2.18)可知，有限元法的关键在于各系数矩阵的生成，而各系数矩阵的表达形式是类似的。在此以质量矩阵M为例探讨一种简单有效的计算方法。由(2.19)可知，计算需要首先判断i点和j点是否在同一个单元中，如果是，则找出共用这两个节点的单元数进行计算。一种常用的思路是先查找整个单元节点信息矩阵，找出以i为顶点的单元，再在这些单元中查找含有另一个顶点j的单元。这种做法比较容易理解，但当剖分单元数较多时会消耗大量的时间。一种简单有效的计算方法是首先针对单个单元计算其三个顶点两两之间的质量矩阵（共9个），再将各个单元的质量矩阵按节点编号对应到矩阵中进行累加。例如，某单元三个节点的编号分别是2、5、6，则可以按如图3.1所示的方式进行对应。这种方法绕开了遍历查找操作，代之以按序号查询，大大缩短了计算。因此按这种算法编写程序。图3.1 Mij矩阵的形成3.4 处理特殊节点特殊节点指的是已知污染物浓度或水头值的节点。为方便起见，程序采用混合节点编号。当生成NpNp维的系数矩阵后，应进行行循环，判断是否为已知点，如果是已知点，则删去这一行，得到mNp维矩阵。再进行列循环，如果为已知点，则删去这一列，并将该列的系数向量值累计到右端项中，得到mm维矩阵。这样形成的系数矩阵有m行m列，方程个数也是m，可以通过直接法或迭代法求解m个未知值。3.5 程序设计流程如前所述，地下水运动模型和污染物迁移模型在形式上是一致的，因此两者的程序设计流程也是相同的，如图3.2所示。3.6 源程序 详见附录。剖分研究域读入单元信息设定区域参数、边界条件、源汇项划分时间段、设定时间步生成系数矩阵(NpNp)、右端列向量(Np1)时间段循环转化系数矩阵(mm)、右端列向量(m1)时间步循环求解目标向量(m1)扩充目标向量(Np1)输出结果图3.2 程序设计流程图四 程序验证4.1 地下水运动数学模型程序验证采用Theis公式对水流部分进行验证。其中Theis公式用Flow2D程序实现。检验问题如下：研究域为8080m的正方形区域；K=1m/d；M=1m；m*=0.01；H0=50m；井位于研究域的中央；开采量Q=80 m3/d；分成6个时段，时段长度为Dt=0.1，0.2，0.4，0.6，0.8，1.0 d，每个时段等分成10个时间步。下面仅以累积时间分别为0.3d，1.3d和2.1d时的数据为例，将有限元程序与Theis公式模拟结果进行比较，如图4.14.6所示：图4.1 有限元法模拟结果(0.3d) 图4.2 Theis公式计算结果(0.3d) 图4.3 有限元法模拟结果(1.3d) 图4.4 Theis公式计算结果(1.3d) 图4.5 有限元法模拟结果(2.1d) 图4.6 Theis公式计算结果(2.1d)可以看出，地下水运动数学模型的有限元法数值模拟结果与Theis公式解析计算结果符合得很好，证实了源程序的正确性。4.2 污染物迁移数学模型程序验证使用下面的等速一维流场中一端注入污染物的一维对流弥散问题解析解对污染物迁移模型进行检验。解析解公式式中，C为浓度，ML-3；C0为注入污染物的浓度，ML-3；DL为水动力弥散系数，L2T-1；u为水流实际速度，LT-1；x为距注入点的距离，L；t为计算时间；erfc(y)为误差余函数，Matlab软件中实现了此函数，可以使用。检验问题如下：研究域为x=0，20m；初始时刻全域无污染物。DL=1.0 m2/d；u=0.2m/d；n=0.2；C0=10.0g/m3；Rd=1；l=0.0。t=10d，分成10个时段，时段长度为Dt=0.1，0.15，0.25，0.40，0.60，0.85，1.20，1.60，2.10，2.75。每个时段等分成5个时间步（要保证Peclet数=1.0e-6 x0=y; y=lw*x0+f; n=n+1;end污染物迁移数学模型主程序function pollutant %主程序data_water;data_pollutant;%原始数据M=zeros(Np); %初始化M矩阵D=zeros(Np); %初始化D矩阵U=zeros(Np); %初始化U矩阵W=zeros(Np); %初始化W矩阵F=zeros(Np,1); %初始化F矩阵area=zeros(Ne,1); %初始化面积矩阵for Tt=1:Nt %时间段循环 dt=TT(Tt,1)/N; %时间步长for n=1:Ne t=beta(n,:); i=t(1); j=t(2); k=t(3); x1=xy(beta(n,1),1); y1=xy(beta(n,1),2); x2=xy(beta(n,2),1); y2=xy(beta(n,2),2); x3=xy(beta(n,3),1); y3=xy(beta(n,3),2); area(n)=abs(det(1,x1,y1;1,x2,y2;1,x3,y3)/2; %求解单元面积 M(i,i)=M(i,i)+Rd(n)*area(n)/6; M(j,j)=M(j,j)+Rd(n)*area(n)/6; M(k,k)=M(k,k)+Rd(n)*area(n)/6; M(i,j)=M(i,j)+Rd(n)*area(n)/12; M(j,i)=M(j,i)+Rd(n)*area(n)/12; M(j,k)=M(j,k)+Rd(n)*area(n)/12; M(k,j)=M(k,j)+Rd(n)*area(n)/12; M(k,i)=M(k,i)+Rd(n)*area(n)/12; M(i,k)=M(i,k)+Rd(n)*area(n)/12; %迭代求M矩阵 D(i,i)=D(i,i)+(Dxy(n,1)*(xy(j,2)-xy(k,2)*(xy(j,2)-xy(k,2)+Dxy(n,2)*(xy(k,1)-xy(j,1)*(xy(k,1)-xy(j,1)/(4*area(n); D(j,j)=D(j,j)+(Dxy(n,1)*(xy(k,2)-xy(i,2)*(xy(k,2)-xy(i,2)+Dxy(n,2)*(xy(i,1)-xy(k,1)*(xy(i,1)-xy(k,1)/(4*area(n); D(k,k)=D(k,k)+(Dxy(n,1)*(xy(i,2)-xy(j,2)*(xy(i,2)-xy(j,2)+Dxy(n,2)*(xy(j,1)-xy(i,1)*(xy(j,1)-xy(i,1)/(4*area(n); D(i,j)=D(i,j)+(Dxy(n,1)*(xy(j,2)-xy(k,2)*(xy(k,2)-xy(i,2)+Dxy(n,2)*(xy(k,1)-xy(j,1)*(xy(i,1)-xy(k,1)/(4*area(n); D(j,i)=D(j,i)+(Dxy(n,1)*(xy(k,2)-xy(i,2)*(xy(j,2)-xy(k,2)+Dxy(n,2)*(xy(i,1)-xy(k,1)*(xy(k,1)-xy(j,1)/(4*area(n); D(j,k)=D(j,k)+(Dxy(n,1)*(xy(k,2)-xy(i,2)*(xy(i,2)-xy(j,2)+Dxy(n,2)*(xy(i,1)-xy(k,1)*(xy(j,1)-xy(i,1)/(4*area(n); D(k,j)=D(k,j)+(Dxy(n,1)*(xy(i,2)-xy(j,2)*(xy(k,2)-xy(i,2)+Dxy(n,2)*(xy(j,1)-xy(i,1)*(xy(i,1)-xy(k,1)/(4*area(n); D(k,i)=D(k,i)+(Dxy(n,1)*(xy(i,2)-xy(j,2)*(xy(j,2)-xy(k,2)+Dxy(n,2)*(xy(j,1)-xy(i,1)*(xy(k,1)-xy(j,1)/(4*area(n); D(i,k)=D(i,k)+(Dxy(n,1)*(xy(j,2)-xy(k,2)*(xy(i,2)-xy(j,2)+Dxy(n,2)*(xy(k,1)-xy(j,1)*(xy(j,1)-xy(i,1)/(4*area(n); %迭代求D矩阵 ux=Kxy(n,1)*(HH(i,Tt)*(xy(j,2)-xy(k,2)+xy(j)*(xy(k,2)-xy(i,2)+HH(k,Tt)*(xy(i,2)-xy(j,2)/(2*n(n)*area(n); uy=Kxy(n,2)*(HH(i,Tt)*(xy(k,1)-xy(j,1)+xy(j)*(xy(i,1)-xy(k,1)+HH(k,Tt)*(xy(j,1)-xy(i,1)/(2*n(n)*area(n); %迭代求水流实际流速 U(i,i)=U(i,i)+(ux*(xy(j,2)-xy(k,2)+uy*(xy(k,1)-xy(j,1)/6; U(j,j)=U(j,j)+(ux*(xy(k,2)-xy(i,2)+uy*(xy(i,1)-xy(k,1)/6; U(k,k)=U(k,k)+(ux*(xy(i,2)-xy(j,2)+uy*(xy(j,1)-xy(i,1)/6; U(i,j)=U(i,j)+(ux*(xy(k,2)-xy(i,2)+uy*(xy(i,1)-xy(k,1)/6; U(j,i)=U(j,i)+(ux*(xy(j,2)-xy(k,2)+uy*(xy(k,1)-xy(j,1)/6; U(j,k)=U(j,k)+(ux*(xy(i,2)-xy(j,2)+uy*(xy(j,1)-xy(i,1)/6; U(k,j)=U(k,j)+(ux*(xy(k,2)-xy(i,2)+uy*(xy(i,1)-xy(k,1)/6; U(k,i)=U(k,i)+(ux*(xy(j,2)-xy(k,2)+uy*(xy(k,1)-xy(j,1)/6; U(i,k)=U(i,k)+(ux*(xy(i,2)-xy(j,2)+uy*(xy(j,1)-xy(i,1)/6; %迭代求U矩阵 W(i,i)=W(i,i)+r(n)*area(n)/6; W(j,j)=W(j,j)+r(n)*area(n)/6; W(k,k)=W(k,k)+r(n)*area(n)/6; W(i,j)=W(i,j)+r(n)*area(n)/12; W(j,i)=W(j,i)+r(n)*area(n)/12; W(j,k)=W(j,k)+r(n)*area(n)/12; W(k,j)=W(k,j)+r(n)*area(n)/12; W(k,i)=W(k,i)+r(n)*area(n)/12; W(i,k)=W(i,k)+r(n)*area(n)/12; %迭代求W矩阵 F(i)=F(i)+I1(n)*area(n)/3+I2(n)/3; F(j)=F(j)+I1(n)*area(n)/3+I2(n)/3; F(k)=F(k)+I1(n)*area(n)/3+I2(n)/3; %迭代求F矩阵endF1=F;M1=M;A=D+U+W+M./dt;A=sparse(A);a=0;b=0;for i=1:Np a=a+1; b=b+1; if ID(a)=0; A(b,:)=; M(b,:)=; F(b)=; F=F-C1(a)*A(:,b)+C1(a)*M(:,b)/dt; A(:,b)=; M(:,b)=; b=b-1; endend %M、A转变为mm矩阵，F转为m维向量for t=1:N %