溶质运移理论-（五）水动力弥散方程的数值解法

地下水流数值模拟方法,第六章 水动力弥散方程的数值解法,中国地质大学环境学院2014春,2,一、有限差分法,基本思想： 将空间划分成若干网格，把时间分成若干小时段，每一个网格中心点（格点）处的未知变量（H或C）视为该网格平均值。利用差商代替微商。基本步骤： （1）剖分（2）建立差分方程组（3）求解,3,一、有限差分法-导数的有限差分近似,图中，去x轴上任意一点i，其坐标为 在改点左右相聚为 处分别取（i-1）和（i+1），其坐标分别为 和以i为中心，泰勒展开C(x),4,一、有限差分法-导数的有限差分近似,整理并略去余项（6-1）-（6-2），再除以 略去余项（6-1）-（6-2），再除以 略去余项：,分别一阶导数的向前差分、向后差分及中心差分,二阶中心差分,5,一、有限差分法-一维水动力弥散的差分解法,一维水动力弥散方程,（6-7）,（1）显格式,式中仅有一个未知数,，解得,式（6-7）中的对流项取中心差分,6,可以证明，稳定性准则要求,即,即,。,（1）,（2）,。,由条件（2）,，格距要求很小；由（2）可知，鉴于,较小，导致,不能太大，将（2）代入（1）式中，得到,7,若对流项改为向后差分,解得,稳定性要求,不难看出，稳定性限制比对流项取中心差分有所改善。,8,（2）隐式格式,整理后得到,隐格式是无条件稳定的。,9,（2）Grank-Nicolson格式,整理后得到,取隐式和显示的平均，即Grank-Nicolson格式,10,一、有限差分法-二维水动力弥散的差分解法,二维水动力弥散方程,（4-56）,（1）显格式,式中仅有一个未知数,式（4-56）中的对流项取中心差分,11,一、有限差分法-二维水动力弥散的差分解法,化简后，有,涉及以（i，j）为中心的5个网格点在tn时刻的已知浓度,12,一、有限差分法-二维水动力弥散的差分解法,（2）隐格式,式（4-56）中右端的对流项取中心差分，右端个C的时阶均取n+1水平,13,一、有限差分法-二维水动力弥散的差分解法,（2）隐格式,整理后,收敛且无条件稳定涉及以（i，j）为中心的5个网格点在tn+1时刻的未知浓度,14,一、有限差分法-二维水动力弥散的差分解法,（3）Grank-Nicolson格式,将隐式格式的两式相加除以2,15,一、有限差分法-二维水动力弥散的差分解法,（3）Grank-Nicolson格式,整理得,16,一、有限差分法-二维水动力弥散的差分解法,（4）交替方向隐式法（ADI法）,分两次对三对角矩阵求逆，将一个t分成两个t/2计算 第一个半时间步，对x方向的偏导数采用隐式差分，对y方向的偏导数采用显示差分。,17,一、有限差分法-二维水动力弥散的差分解法,（4）交替方向隐式法（ADI法）,整理得,第二个半时间步，对y方向的偏导数采用隐式差分，对x方向的偏导数采用显示差分。,18,一、有限差分法-二维水动力弥散的差分解法,（4）交替方向隐式法（ADI法）,整理得,收敛且无条件稳定,19,一、有限差分法-求解差分方程的计算机程序举例,算例 见教材P81-85,20,二、有限单元法-伽辽金有限单元法,伽辽金法 对均匀多孔介质，一维稳定流场中二维水动力弥散方程 取x轴方向与地下水流向相同，记伽辽金法是寻找一个级数形式的试探函数作微分方程的近似解，并使其满足边界条件,（6-26）,21,二、有限单元法-伽辽金有限单元法,上式一般不会满足方程，因为 仅是近似解，得到一个剩余误差函数 R(x，y) ，在平均意义下迫使误差为0 我们加以权，令剩余的加权积分为0，W是一组权函数。加权剩余法,（6-29）,22,二、有限单元法-伽辽金有限单元法,在整个研究区D上，基函数NL(x,y)分段定义(1)函数区域连续（2）一阶导数不连续（3）二阶导不容易确定故采用分步积分的方法使二阶导数降阶,23,二、有限单元法-伽辽金有限单元法,对一维积分整理后分步积分推广到二维的情况,（6-32）,代入（6-29）,24,二、有限单元法-伽辽金有限单元法,有限单元剖分与基函数（1）三角形单元 见地下水流动问题数值模拟（2）矩形单元将区域记为D，边界记为B。要求各单元均质、等厚，即T、为常数。结点（内结点、边界结点）,25,二、有限单元法-伽辽金有限单元法,构造单元函数NL基函数取为 “双线性插值”,基函数NL（x，y）形状如同一顶高等于1、有4条直线斜边和4条下凹型曲边的尖顶斗笠,26,二、有限单元法-伽辽金有限单元法,子区DL以结点L为其共同结点的所有矩形单元（0）,27,二、有限单元法-伽辽金有限单元法,典型矩形单元 该单元中心位于坐标原点处，且4条边分别平行x，y轴 ,结点从左下角开始按逆时针方向编号i,j,m,n。沿x方向长2a，y方向宽2b。,28,二、有限单元法-伽辽金有限单元法,按上述要求所构造的NeL形式为对典型单元，令 则,29,二、有限单元法-伽辽金有限单元法,其中,即,（6-34）,30,二、有限单元法-伽辽金有限单元法,Nei(x,y)在结点i处为1，其它3处为0根据上式，有平行x或y方向上，Nei(x,y)线性变化故单元基函数为双线性插值函数若结点坐标用 表示结点浓度用 表示,31,二、有限单元法-伽辽金有限单元法,结点浓度CL(t)和单元基函数NeL(x,y)来确定单元内的近似解 ，即对于区域D，结点L的基函数在子区DL内各单元分块确定。,称为矩阵单元e上的基函数,（6-36）,32,二、有限单元法-伽辽金有限单元法,伽辽金有限单元法 左端积分范围为D，由于在 上NL=0,改写在子区DL的积分，对矩形单元逐个积分、求和。 设子区DL上有neL个单元 ，有,（6-39）,33,二、有限单元法-伽辽金有限单元法,由（6-36）得,由（6-34）得,再由（6-36）得,（6-40）,（6-42）,34,二、有限单元法-伽辽金有限单元法,将（6-40） （6-42） 代入（6-39）得,35,二、有限单元法-伽辽金有限单元法,（6-44）,36,二、有限单元法-伽辽金有限单元法,用矩阵形式表示一阶微分方程组,将一阶导数用差分形式表示,按矩阵符号有,若C在新时间水平t+t上取值，则,37,二、有限单元法-伽辽金有限单元法,（1）给定初始条件和边界条件，可求（6-46）各个系数矩阵，于是可解得第一个 末时刻各结点的浓度 ；（2）以此浓度为已知浓度 ，计算新的系数矩阵 ；（3）再次求解，得到第二个时间步长末时刻结点浓度,（6-46）,浓度 可在 之间任意位置近似表示,38,二、有限单元法-伽辽金有限单元法,按单元顺序计算系数矩阵每个矩阵单元涉及4个结点i，j，m，n，涉及4个方程单元e的作用，可用4行与4列的矩阵来表示，称为单元矩阵,39,二、有限单元法-伽辽金有限单元法,算例：单元弥散矩阵与整体弥散矩阵的关系设研究区划分为8个单元，16个结点，四周为隔水边界编号如下：,40,二、有限单元法-伽辽金有限单元法,算例：单元弥散矩阵与整体弥散矩阵的关系计算每个单元弥散矩阵，按双下标已知，放置单元弥散矩阵并叠加，得到总体弥散矩阵。,41,二、有限单元法-伽辽金有限单元法,从（6-44）知,将坐标原点移到矩形单元中心,有,42,二、有限单元法-伽辽金有限单元法,通过双下标一致的元素求和，得到总体弥散矩阵诸元素,43,二、有限单元法-伽辽金有限单元法,积分求解可用高斯求积方法，详见P97-105,44,二、有限单元法-伽辽金有限单元法,45,二、有限单元法-伽辽金有限单元法,算例：二维水动力弥散的伽辽金法计算机程序 数学模型写成,剖分成9个单元共20个结点，空间步长取x=y=10m，时间步长t=5d，渗透流速u=0.01m/d 具体代码见P106-111,46,三、过程与数值弥散,过量现象 如一维流动水动力弥散中单位浓度梯度函数注入的情况。 对比解析解与数值解C-x曲线，在浓度封面附近数值计算的浓度超过最大浓度值1和小于最小浓度0，失去物理意义。,数值弥散 若纯对流方程采用有限差分逼近，其结果呈现似有弥散的存在。,47,三、过程与数值弥散,48,三、过程与数值弥散,过量现象-解析 由于时间步长与空间尺度匹配不佳，因而含水层在数值上不能“吸收”住污染物所引起。解决办法： 将时间增量选择为几何级数的子项； 做试验校正时间步长和差分网格；,49,三、过程与数值弥散,数值弥散-解析 由截断误差引起。 对流项一阶空间导数采用向后差分逼近，浓度Ci-1泰勒展开 变形得 向后差分的截断误差主部分为,相当于弥散系数的弥散项,50,三、过程与数值弥散,数值弥散-解析 由故纯对流方程有限差分近似表示后，结果如同对流-弥散方程的结果，记用数值弥散系数与物理弥散系数的比值来评价数值弥散对物理弥散的干涉程度,此条件下的数值弥散系数,51,三、过程与数值弥散,数值弥散-解析 定义为Peclet数， 时，弥散占优； 时，对流占优。速度u与物理弥散系数DL比越大， 越大。对于小 数，数值解与精确解非常接近，对于大 数，数值解过渡带变宽，浓度锋面变平缓，在锋面的前后还出现解的振动过量现象。,控制x的大小从而减少数值弥散,52,三、过程与数值弥散,对纯对流方程,53,三、过程与数值弥散,对纯对流方程,结果也造成数值弥散,此时，总的数值弥散是两者之和,54,四、对流占主导地位时的数值方法,上游加权法,对流项的有限差分化过程中乘上一个适当的权因子，则波动和过量得到改善或消除,一维水动力弥散中对流项取中心隐式差分,若相邻格点i-1或（和）i+1的浓度越大，则格点i的浓度也越大，反映到差分方程上，系数Gi-1和Gi+1不得与Gi同号，上式要求 不能满足时，对流占优势,55,四、对流占主导地位时的数值方法,上游加权法,假定DL和u为常数，取等格距差分网格。一维对流-弥散方程在 积分，有对流-弥散通量用有限差分近似表示,56,四、对流占主导地位时的数值方法,上游加权法,为权因子；当=1/2时，对流项相当于中心差分；当=1时，对流项相当于取向后差分， 当=0时，对流项相当于向前差分。速度不同， 取值不同,57,四、对流占主导地位时的数值方法,上游加权法,上游浓度权因子大于1/2。物理意义：加强上游格点浓度值在对流项中的作用可通过浓度的线性（或二次）插值具体确定在t期间C在格边上随时间的变化，从而确定值。,58,四、对流占主导地位时的数值方法,仍用Gi-1、Gi和Gi+1分别表示上式左端第一、二、三项括号内系数。则采用上游加权法有限差分格式，可消除过量和波动，但截断误差增大，数值弥散反而增大。中心差分格式，截断误差为二阶上游加权有限差分截断误差为一阶,59,四、对流占主导地位时的数值方法,上游权因子的物理实质：对流项浓度C在某断面，例如 断面的中心差分意味着对流项的浓度在t时段内均以i与i-1两格点的平均浓度，即 通过 断面；实际为： t时段的初始时刻是该平均浓度值，随后在对流的作用下，该断面的浓度显然会向上游格点的浓度Ci-1 (u0时)的变化。,60,四、对流占主导地位时的数值方法,特征值方法 令源汇项为0，展开对流-弥散方程,61,四、对流占主导地位时的数值方法,特征值方法 假定流体密度值不变（不可压缩），含水层不可压缩 ，则 ，（6-75）简化为基本思路：物质导数表示一个沿轨迹（特征曲线） 运动的研究对象浓度变化率。由弥散作用引起。若存在运会想，还将受到沿轨迹源与汇、化学、生物化学反应、吸附与解析等。对流项用沿特征线的示踪质点的运动描述,（6-76）,62,四、对流占主导地位时的数值方法,特征值方法 特征点均匀分布研究区，赋一个初始浓度。 含水层中未受污染部分特征点的初始浓度为零，特征点浓度眼轨迹的瞬时浓度变化，借助矩形网格差分格式计算。,63,四、对流占主导地位时的数值方法,64,四、对流占主导地位时的数值方法,模型精确度 取决于网格距大小。可以在同一网格上计算速度场。如果采用较少的网格距，可以改善流场的描述。 坐标为 的特征点所在的单元，可用下标i,j确定 精确度受到每一网格内所选择的最初的特征点数目影响。,65,四、对流占主导地位时的数值方法,缺点： （1）实际程序冗长，需记录特征点的踪迹； （2）抽水点或出流边界的特征点需要消除，入渗边界须有特征点生成，不透水边界特征点须反射回来； （3）要通过消去特征点来避免滞留点附近特征点堆积； （4）特征点耗尽的单元，必须产生新的特征点。,66,四、对流占主导地位时的数值方法,动坐标系方法与网格变形方法 能消除过量现象，还保持较陡的浓度锋面。 对式 设坐标变换 得， 取空间步长 ，则,67,四、对流占主导地位时的数值方法,动坐标系方法与网格变形方法 在固定坐标系中取格距为 ，有 动坐标系中格点运动与固定坐标系格点运动重合。 设固定坐标系中格点xj在tn+1时刻与动坐标系中Xi重合，有,动坐标系下浓度结点的表达式,68,四、对流占主导地位时的数值方法,动坐标系方法与网格变形方法 式变成 对半无限长砂柱，当DL=0时,69,四、对流占主导地位时的数值方法,网格变形方法 设tn时刻有限元结点在空间中的分布是对应浓度分布为,70,四、对流占主导地位时的数值方法,网格变形方法 主要步骤：结点位置随时间变化，基函数及有限元方程的系数都依赖于时间，每推进一个 t都需要重新计算。二维问题还须进行畸形判断,71,四、对流占主导地位时的数值方法,结合动坐标的网格变形方法 一维对流-弥散在某个动点在t时刻位置x可写成 其中，x0是动点的初始位置，有取时间的一阶偏导数对流-弥散方程写成,若控制动点速度dx/dt 并让它接近流速u，则可转换成弥散为主的的方程式,72,四、对流占主导地位时的数值方法,若使用伽辽金法，基函数因依赖于结点位置而成为时间的函数。若用 表示随时间变化的基函数，则对任意t，有 故,73,四、对流占主导地位时的数值方法,随机步行法 采用示踪剂描述污染物运移。只有被污染的特征点在流场中运动。每个特征点都给定一个不变的污染物质量，此质量的和等于排进含水层的总污染物质量。 取许多单个特征点轨迹（随机步行）的平均值，可得到一个有弥散的特征点的分布。要得到一个浓度分布，就要先叠加一个网格并算每个网格单元所含有的特征点数。,74,四、对流占主导地位时的数值方法,随机步行法-以一维运移为例在x=0处，一个瞬时注入的、质量为 M的理想失踪及，在t时刻浓度分布对一个固定时间t，看做正态分布,75,四、对流占主导地位时的数值方法,随机步行法-以一维运移为例在时间t=0和位置x=0处，放置具有污染物质量为M/N的N个特征点，每一个特征点移动距离x，而在时刻t到达它们各自位置，即Z是一个平均值为0，标准差为1的正态分布随机变量 所得的频率分布，通过一个标准化因子，有对有限数量的颗粒，可近似去顶C(x,t)，将空间变量划分为长度x的若干区间，区间中点浓度,76,四、对流占主导地位时的数值方法,随机步行法-以一维运移为例若孔隙平均流速是时间和位置的函数，那么时间为0开始的颗粒P在时间间隔t内的随机轨迹xp(t)为,77,四、对流占主导地位时的数值方法,随机步行法-二维流场 首先假定流动方向平行于x轴，得到Z和Z是正态分布的随机变量的两个值。对任意方向的流动，需考虑弥散张量特性。,78,四、对流占主导地位时的数值方法,随机步行法-二维流场 首先假定流动方向平行于x轴，得到Z和Z是正态分布的随机变量的两个值。对任意方向的流动，需考虑弥散张量特性。,对流运动可确定在水流流动方向及其正交方向的弥散运动,79,四、对流占主导地位时的数值方法,随机步行法,80,四、对流占主导地位时的数值方法,随机步行法 数值结果的质量主要取决于特征点数目N，要考虑网格的选择。大网格间距会产生极平均的结果，小网格间距会产生十分粗糙的分布。 在稳定流动和一个源或数个源都有固定强度的情况下，可以节省很多特征点。,81,82,四、对流占主导地位时的数值