本发明专利技术公开了一种地下水溶质运移数值模拟中减少误差的方法，其是在水动力弥散系数上加上一个数值弥散估算值，得到一个修正弥散系数，用其替代方程中有明确物理意义的水动力弥散系数进行计算。并提出了一个新的参数－数值弥散因子（ＮＤＦ），可以根据研究需要进行参数分区并适当调节该因子的大小，从而达到控制数值振荡，减小数值弥散的目的。从一维到二维的多个数值算例的模拟计算结果表明，新方法能在消除数值振荡的基础上，较好地减少数值弥散，达到令人满意的精度。
【技术实现步骤摘要】
本专利技术涉及一种数值模拟中减少误差的方法，具体是一种在。
技术介绍
解决环境污染控制问题的一个重要技术过程就是将所研究的环境系统行为抽象为数学模型，这是进行定量研究工作的基础。在水环境污染模拟问题中，对流-弥散方程是溶质(污染物)迁移的基本方程，其求解结果的精度可能会影响对水环境系统的分析，其被广泛用于对溶质输运问题的刻画。近年来，国内外学者对对流-弥散方程的数值解法进行了广泛而深入的研究，提出了许多方法。这些方法各有利弊，当对流-弥散方程中的弥散项占主导地位时，使用各种数值方法均能得到较为满意的结果，反之，若对流项占主导地位，各种常见的数值解法都会遇到困难(孙讷正.地下水污染-数学模型和数值方法.北京地质出版社，1989)。因此，对流-弥散方程的数值解通常被认为是“令人困窘”的难题。 对流-弥散方程的数值法大体分三类Euler法、Lagrange法和Euler-Lagrange法(成建梅，胡进武.饱和水流溶质运移问题数值解法综述.水文地质工程地质，2003.99.)。各种数值方法都包含程度不同的数值弥散和振荡两类误差，这两类误差有共同的形成原因。Pinder指出，当-种差分格式使数值弥散最小时，就会出现解的振荡；而当控制了解的振荡，就要以增大数值弥散作为代价(Pinder GF，&Gray，W.G.Finiteelement simulation in surface and subsurface hydrology.Academic Press，1977.)。常用的Euler法主要包括有限差分法，有限单元法(Hossain M.A.，&Barber，M.E.Optimized Petrov-Galerkin model for advective-dispersive transport.Applied Mathematics and Computation，2000.115(1)，1-10.)等。这类方法应用比较方便，但易受数值振荡和数值弥散的影响。本专利技术主要是在Euler法中的上游加权法基础上进行改进。 为减少数值振荡与弥散，提高计算精度，通常采用两类方法。一类方法是在空间离散时加密网格，并减小时间步长，但这种方法有时在实践上会有困难，特别是对于一些区域性的污染问题，过小的空间步长和时间步长有时难以做到(薛禹群，谢春红.地下水数值模拟.北京科学出版社，2007.)；另外一类方法是使用更高精度的差分格式，如20世纪80年代发展起来的一类流体力学计算的差分格式TVD(Total VariationDiminishing)(水鸿寿.一维流体力学差分方法.北京国防工业出版社，1998.)，这类高阶方法常常减少了数值弥散，但是却引入了振荡。又有学者提出了通量限制器来消除振荡，保存陡的浓度锋面。如由Leonard等人开发的ULTIMATE(对流运移方程瞬时插值模拟的通用限制器)(Zheng C.，& Bennett G.D.Applied Contaminant TransportModeling.San FranciscoWiley Inter-science，2002.)，这种格式不会引入过多的数值弥散，且基本没有振荡。但在实际计算时耗时过多。 数值弥散一般认为是由时间和空间导数差分后的二阶截断误差引起的，其大小可以通过分析二阶截断误差近似得到。数值弥散使一个以弥散系数为D的对流-弥散问题变成为一个有新的弥散系数D*的对流-弥散问题(Noorishad J.，Tsang C.F.，Perrochet，P.，&Musy，A.A perspective on the numerical solution of convection-dominatedtransport problemsA price to pay for the easy way out.Water ResourcesResearch，1992.28(2)，551-561.)。
技术实现思路
本专利技术所要解决的技术问题是提供一种运用对流-弥散方程对地下水溶质运移进行数值模拟过程中，在保证无数值振荡的前提下，可以较大程度地减少数值弥散的方法。 本专利技术所述地下水中溶质运移数值模拟中减少误差的方法，包括以下步骤 1)从水文地质概念模型出发，用对流-弥散方程描述地下水流运动和溶质运移的基本特征，对流-弥散方程的一维溶质运移数学模型为 0≤x≤L，0≤t≤Tsum (1) C(x，0)＝Cx(0) C(0，t)＝C0，C(l，t)＝Cl 式中参数，C为溶质(污染物)浓度，t为时间，DL为弥散系数，u为实际平均流速，Cx(0)为初始时刻x处的浓度，C0为边界x＝0上的浓度，Cl为边界x＝l上的浓度； 2)为了得到(1)式左端浓度的时间导数项的差分格式，对j时刻和j+1时刻在i点的浓度进行泰勒展开，并引入时间权因子θ得 (2) 将(2)式右端第二项浓度对时间的二阶导数项可以转化成对空间的导数 (3) (3)式是个加权六点格式，当θ＝0时，为显式格式，当θ＝1/2时，为Crank-Nicolson格式，当θ＝1时，为隐式格式； 3)对(3)式的对流项使用泰勒展开得到第i和i+1节点的表达式，并引入空间加权因子α得 (4) 上式右端第三项为二阶截断误差； 4)将(3)、(4)式代回(1)式，且将弥散项用传统差分格式代替，得 (5) 其中截断误差εL为 (6) 其中Cr为科朗数，(6)式中二阶截断误差大小与时间加权因子θ，空间加权因子α，空间步长Δx，地下水实际平均流速u，以及时间步长t有关，且时间步长越长，空间步长以及地下水实际平均流速越大，则二阶截断误差越大，亦即数值弥散越大； 令D′＝uΔx，则 (6)式右端高于二阶的截断误差，以及空间与时间的离散对计算结果的影响可以通过在D’上乘以一个因子体现出来，定义这个因子为NDF，且令Dnum＝NDF×D′ (7) (8) 将(8)代回上游加权差分格式，令修正弥散系数D*＝D-Dnum，得到新计算格式如下 (9) 本专利技术提出了一个新的参数数值弥散因子(NDF)。其取值范围在0～1之间，根据不同的问题其取值不同。数值计算结果显示，使用修正弥散系数对传统的隐式格式有较大的改进。在保持无数值振荡的前提下，可以较大程度地减少数值弥散。本专利技术既不需要以加密网格为代价，又不需要引入复杂的差分格式，可以比较简单、方便地提高模拟计算精度。 附图说明 图1对隐式差分格式使用新方法计算结果对比，图1a的NDF＝1，图1b的NDF＝0.7～0.8 图2Peclet数为50时新旧方法计算结果对比，图2a是解析法、传统上游加权算法以及新方法(修正弥散系数取NDF＝1)计算结果对比，图2b是解析法、传统上游加权算法以及新方法(修正弥散系数取分别NDF＝1和NDF＝0.9)计算结果对比， 图3Peclet数为100与150时新旧方法对比，图3a是Peclet数为100，图3b是Peclet数为150， 图4本方法与旧方法计算结果对比图(lst year)(b本文档来自技高网...
【技术保护点】
一种地下水溶质运移数值模拟中减少误差的方法，其特征在于包括以下步骤： １）从水文地质概念模型出发，用对流－弥散方程描述地下水流运动和溶质运移的基本特征，对流－弥散方程的一维溶质运移数学模型为： ＊Ｃ／＊ｔ＝Ｄ↓［Ｌ］＊↑［２］Ｃ ／＊ｘ↑［２］－ｕ＊Ｃ／＊ｘ ０≤ｘ≤Ｌ，０≤ｔ≤Ｔ↓［ｓｕｍ］ （１） Ｃ（ｘ，０）＝Ｃ↓［ｘ］（０） Ｃ（０，ｔ）＝Ｃ↓［０］，Ｃ（ｌ，ｔ）＝Ｃ↓［ｌ］ 式中参数，Ｃ为溶质浓度，ｔ为时间，Ｄ↓［Ｌ］为弥散系数，ｕ为实 际平均流速，Ｃ↓［ｘ］（０）为初始时刻ｘ处的浓度，Ｃ↓［０］为边界ｘ＝０上的浓度，Ｃ↓［ｌ］为边界ｘ＝ｌ上的浓度； ２）为了得到（１）式左端浓度的时间导数项的差分格式，对ｊ时刻和ｊ＋１时刻在ｉ点的浓度进行泰勒展开，并引入时间权因子θ得 ： ＊＊＊ （２） 将（２）式右端第二项浓度对时间的二阶导数项转化成对空间的导数： ＊＊＊ （３） （３）式是个加权六点格式，当θ＝０时，为显式格式，当θ＝１／２时，为Ｃｒａｎｋ－Ｎｉｃｏｌｓｏｎ格式，当θ＝１时 ，为隐式格式； ３）对（３）式的对流项使用泰勒展开得到第ｉ和ｉ＋１节点的表达式，并引入空间加权因子α得： ＊Ｃ／＊ｘ＝α（Ｃ↓［ｉ］－Ｃ↓［ｉ－１］／Δｘ）＋（１－α）（Ｃ↓［ｉ＋１］－Ｃ↓［ｉ］／Δｘ）＋（２α－１）Δｘ／２＊ ↑［２］Ｃ／＊ｘ↑［２］＋… （４） 上式右端第三项为二阶截断误差； ４）将（３）、（４）式代回（１）式，且将弥散项用传统差分格式代替，得： （Ｃ↓［ｉ］↑［κ＋１］－Ｃ↓［ｉ］↑［κ］）／Δｔ＝θ｛Ｄ／（Δｘ）↑［２］ （Ｃ↓［ｉ＋１］↑［κ＋１］－２Ｃ↓［ｉ］↑［κ＋１］＋Ｃ↓［ｉ－１］↑［κ＋１］）－ｕ／Δｘ［αＣ↓［ｉ］↑［κ＋１］＋（１－α）Ｃ↓［ｉ＋１］↑［κ＋１］－αＣ↓［ｉ－１］↑［κ＋１］－（１－α）Ｃ↓［ｉ］↑［κ＋１］］｝＋（１－θ）｛Ｄ／（Δｘ）↑［２］（Ｃ↓［ｉ＋１］↑［κ］－２Ｃ↓［ｉ］↑［κ］＋Ｃ↓［ｉ－１］↑［κ］）－ｕ／Δｘ［αＣ↓［ｉ］↑［κ］＋（１－α）Ｃ↓［ｉ＋１］↑［κ］－αＣ↓［ｉ－１］↑［κ］－（１－α）Ｃ↓［ｉ］↑［κ］］｝－ε↓［Ｌ］ （５） 其中截断误差ε↓［Ｌ］为： ε↓［Ｌ］＝１／２［（２θ－１）ｕ↑［２］Δｔ＋（２α－１）ｕΔｘ］＊↑［２］Ｃ／＊ｘ↑［２］＋...
【技术特征摘要】
【专利技术属性】
技术研发人员：吴吉春，祝晓彬，梅一，
申请(专利权)人：南京大学，
类型：发明
国别省市：84[中国|南京]
还没有人留言评论。发表了对其他浏览者有用的留言会获得科技券。