CN105512417A - 基于粒子追踪的孔隙地下水污染物三维运移模拟方法 - Google Patents
基于粒子追踪的孔隙地下水污染物三维运移模拟方法 Download PDFInfo
- Publication number
- CN105512417A CN105512417A CN201510947333.XA CN201510947333A CN105512417A CN 105512417 A CN105512417 A CN 105512417A CN 201510947333 A CN201510947333 A CN 201510947333A CN 105512417 A CN105512417 A CN 105512417A
- Authority
- CN
- China
- Prior art keywords
- particle
- delta
- time
- coordinate
- dimensional
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/30—Circuit design
- G06F30/36—Circuit design at the analogue level
- G06F30/367—Design verification, e.g. using simulation, simulation program with integrated circuit emphasis [SPICE], direct methods or relaxation methods
Landscapes
- Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Microelectronics & Electronic Packaging (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
一种基于粒子追踪的孔隙地下水污染物三维运移模拟方法,包括:设定粒子总数、每个粒子的初始坐标、释放时间,设定计算时间步长;在任意时刻、任意位置释放粒子;对每个被释放的粒子,计算其在一个时间步长后的位移;根据计算得到的粒子的新坐标,判断该粒子是否位于计算网格区域内,若已跳出,则通过边界控制使其回到计算区域边界上;重复上述步骤,以此迭代,计算得出每个时间节点上每个粒子的坐标;通过统计特定区域内粒子数量,即可得到该区域的地下水污染物浓度分布。与传统方法相比,本发明的方法可大大降低计算成本,且利用该方法得到的模拟结果与对流弥散方程的解析解有很好的匹配。
Description
技术领域
本发明涉及地下水数值模拟计算领域,特别是涉及一种基于粒子追踪的孔隙地下水污染物三维运移模拟方法。
背景技术
孔隙地下水作为北方重要的饮用水源,其水质直接关系到饮水安全。地下水埋藏在地下,其污染具有很强的复杂性和隐蔽性,对地下水污染的掌握需要借助数值模拟方法,通过模型的高效运算,可以节约大量经济和时间成本,使管理决策人员迅速掌握地下水中污染物浓度变化情况。
地下水数值模拟主要包括地下水水流模拟和溶质运移模拟,前者数值求解地下水水流方程,后者数值求解对流-弥散方程。描述孔隙地下水溶质运移的对流-弥散方程可表述为:
其中C表示地下水中溶质的浓度,v表示地下水流动孔隙速度,D表示水力弥散系数,t表示时间,为微分算子。
孔隙地下水溶质运移模拟求解方法常用欧拉法、拉格朗日法,以及二者的结合。欧拉法以空间中固定坐标系作为参照系,常见的有限差分法和有限单元法属于欧拉法。应用有限元和有限差分方法进行地下水溶质迁移模拟,有两个固有的缺陷,一是当网格Peclet数较大,即对流项强于弥散项时,容易受到数值弥散的影响;二是当模拟区域范围较大时,因计算网格较多,计算成本比较昂贵。
拉格朗日法通过质点追踪,以运动坐标系作为参照系,常见的拉格朗日法有粒子追踪法等。运动坐标系容易造成数值不稳定;另外在质点追踪过程中对质点速度的连续性有较高要求,否则由于速度插值也容易造成局部的质量不守恒。
欧拉-拉格朗日混合法用拉格朗日法解决溶质运移中的对流问题,用欧拉法解决弥散问题,结合了二者的优点,却也同时存在缺点,且计算耗时。
发明内容
有鉴于此,本发明的目的在于提出一种孔隙地下水中污染物三维运移模拟方法,该数值模拟方法是基于拉格朗日法中的粒子追踪方法,以解决模拟过程中的数值弥散问题,模拟得到的数值解可以完全拟合对流-弥散方程的解析解,且与有限元法相比可大大节约计算成本。
为实现本发明的上述目的,本发明提出了一种基于粒子追踪的孔隙地下水污染物三维运移模拟方法,包括以下步骤:
步骤S1:设定粒子总数、每个粒子的初始坐标、释放时间,设定计算时间步长;
步骤S2:在任意时刻、任意位置释放粒子;
步骤S3:对每个被释放的粒子,计算其在一个时间步长后的位移,粒子在地下水流场中的运动由以下方程控制:
其中x、y、z表示粒子的空间坐标,v表示粒子的对流运动速度,Δt表示时间步长,D表示水力弥散系数,Z表示介于0到1之间的随机数,这样就根据粒子在t时刻的坐标计算出其在一个时间步长后,也即t+Δt时刻的新坐标;
步骤S4:根据计算得到的粒子的新坐标,判断该粒子是否位于计算网格区域内,若已跳出,则通过边界控制使其回到计算区域边界上;
步骤S5:重复步骤S2-S4,以此迭代,计算得出每个时间节点上每个粒子的坐标;
步骤S6:通过统计特定区域内粒子数量,即可得到该区域的地下水污染物浓度分布。
基于上述技术方案可知,本发明的孔隙地下水中污染物三维运移模拟方法,将地下水中的溶质抽象为大量粒子,以粒子的运动来模拟溶质在孔隙介质中的运移,其中以粒子的有序运动来刻画溶质因对流引起的迁移过程,以粒子的随机位移来刻画溶质的弥散过程,因此可以精细刻画溶质在地下水中的运移过程。传统的有限元法在求解对流-弥散方程时,需要求解的代数方程个数与计算网格节点数成正比,本发明的方法基于统计物理学中的随机行走粒子追踪,不直接求解对流-弥散方程,而是求解粒子的位移,计算量仅与粒子个数有关,因此在三维情景下与有限元法相比可大大降低计算成本,且运算结果与对流-弥散方程的解析解完全拟合。本发明的方法,粒子不会凭空出现或消失,粒子总数量在输入和输出及整个过程中是不变的,因此该方法从本质上是质量守恒的,可消除计算过程中的数值弥散。
附图说明
图1是瞬时释放污染物在三维均质孔隙介质中的污染羽分布;
图2是粒子到达计算网格区域边界时的边界控制效果示意图;
图3是瞬时释放污染物在三维均质孔隙介质中的突破曲线。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白,以下结合具体实施例,并参照附图,对本发明作进一步的详细说明。
本发明公开了一种基于粒子追踪的孔隙地下水污染物三维运移模拟方法,包括以下步骤:
步骤S1:设定粒子总数、每个粒子的初始坐标、释放时间,设定计算时间步长;
步骤S2:在任意时刻、任意位置释放粒子;
步骤S3:对每个被释放的粒子,计算其在一个时间步长后的位移,粒子在地下水流场中的运动由以下方程控制:
其中x、y、z表示粒子的空间坐标,v表示粒子的对流运动速度,Δt表示时间步长,D表示水力弥散系数,Z表示介于0到1之间的随机数,这样就根据粒子在t时刻的坐标计算出其在一个时间步长后,也即t+Δt时刻的新坐标;
步骤S4:根据计算得到的粒子的新坐标,判断粒子是否位于计算网格区域内,若已跳出,则通过边界控制使其回到计算区域边界上;
步骤S5:重复步骤S2-S4,以此迭代,计算得出每个时间节点上每个粒子的坐标;
步骤S6:通过统计特定区域内粒子数量,即可得到该区域的地下水污染物浓度分布。
作为优选,在步骤S1中,设定粒子初始坐标时,允许多个粒子的初始坐标完全相同;粒子释放时间表示该粒子初次被释放的时间,应为时间步长的整数倍;时间步长需满足库朗数(Courantnumber)小于1。
作为优选,在步骤S2中,对任意一个粒子,若当前步的计算时间大于或等于步骤S1中设定的粒子释放时间,则该粒子被释放;若在某一坐标位置所有粒子的释放时间相同,则该释放过程是瞬时的,若粒子释放时间不同且构成一个连续时段,则该释放过程是持续的。
作为优选,在步骤S3中,粒子的对流运动速度v由地下水流场计算所得,并通过粒子所在单元节点插值获得;Z由计算程序中的随机函数产生;对流弥散系数D由下式计算获得:
其中δij为克罗内克符号,αL为纵向弥散系数,αT为横向弥散系数,Dd为分子扩散系数,vi为i方向平均孔隙速度。
作为优选,在步骤S4中,通过粒子位置坐标与计算网格坐标进行比较,确定粒子是否位于计算网格之内;边界控制的实现是通过计算粒子在该时间步长内的位移轨迹与计算网格边界的交点,将交点坐标作为粒子的新坐标。
作为优选,在步骤S6中,特定区域内粒子数量与被释放粒子总量的比值,即为该区域内地下水污染物的相对浓度。
作为优选,以一定数量的粒子代表地下水中污染物(溶质)的一定浓度,以粒子的位置代表地下水中污染物(溶质)的空间分布。其中,该粒子没有体积和质量;不同的粒子可以按组分类,属于不同组别的粒子可以表示相同类型,也可以表示不同类型的地下水污染物。
下面通过本发明的一个优选实施例对本发明做进一步说明。
如图1所示,选取正方体结构的计算区域,该区域为均质孔隙介质,地下水流向为对角线方向。在该区域上游即正方体一个顶点附近瞬时释放粒子,模拟地下水中污染物在该区域的分布。具体步骤如下:
步骤S1-S2:设定粒子总数为500,每个粒子的初始坐标为(0.01,0.01,0.01),在初始时刻瞬时释放;设定时间步长为一个粒子通过一个计算网格长度所需时间的1/10;
步骤S3:计算一个时间步长后每个粒子的位移,根据粒子在t时刻的坐标计算出其在一个时间步长后,也即t+Δt时刻的新坐标;
步骤S4:根据计算得到的粒子的新坐标,判断粒子是否位于计算网格区域内,若已跳出,则通过边界控制使其回到计算区域边界,如图2所示;
步骤S5:重复步骤S2-S4,以此迭代,计算得出每个时间节点上每个粒子的坐标;
步骤S6:通过统计特定区域内粒子数量,该数量与500的比值即可认为是任意位置的溶质浓度。
由该区域下游即正方体对角顶点处的污染物浓度变化,获取突破曲线,模拟结果与解析解对比如图3所示,其中红色曲线为解析解,蓝色曲线为模拟结果。可见该方法获得的模拟结果与解析解完全拟合。
以上所述的具体实施例,对本发明的目的、技术方案和有益效果进行了进一步详细说明,应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限制本发明,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (9)
1.一种基于粒子追踪的孔隙地下水污染物三维运移模拟方法,包括以下步骤:
步骤S1:设定粒子总数、每个粒子的初始坐标、释放时间,设定计算时间步长;
步骤S2:在任意时刻、任意位置释放粒子;
步骤S3:对每个被释放的粒子,计算其在一个时间步长后的位移,粒子在地下水流场中的运动由以下方程控制:
其中x、y、z表示粒子的空间坐标,v表示粒子的对流运动速度,Δt表示时间步长,D表示水力弥散系数,Z表示介于0到1之间的随机数,这样就根据粒子在t时刻的坐标计算出其在一个时间步长后,也即t+Δt时刻的新坐标;
步骤S4:根据计算得到的粒子的新坐标,判断该粒子是否位于计算网格区域内,若已跳出,则通过边界控制使其回到计算区域边界上;
步骤S5:重复步骤S2-S4,以此迭代,计算得出每个时间节点上每个粒子的坐标;
步骤S6:通过统计特定区域内粒子数量,即可得到该区域的地下水污染物浓度分布。
2.如权利要求1所述的孔隙地下水污染物三维运移模拟方法,其特征在于,所述步骤S1中,设定粒子初始坐标时,允许多个粒子的初始坐标完全相同;粒子释放时间表示该粒子初次被释放的时间,应为时间步长的整数倍;时间步长需满足库朗数小于1。
3.如权利要求1所述的孔隙地下水污染物三维运移模拟方法,其特征在于,在步骤S2中,对任意一个粒子,若当前步的计算时间大于或等于步骤S1中设定的粒子释放时间,则该粒子被释放;若在某一坐标位置所有粒子的释放时间相同,则该释放过程是瞬时的,若粒子释放时间不同且构成一个连续时段,则该释放过程是持续的。
4.如权利要求1所述的孔隙地下水污染物三维运移模拟方法,其特征在于,在步骤S3中,粒子的对流运动速度v由地下水流场计算得到,并通过粒子所在单元节点插值获得;Z由计算程序中的随机函数产生;对流弥散系数D由下式计算获得:
其中,δij为克罗内克符号,αL为纵向弥散系数,αT为横向弥散系数,Dd为分子扩散系数,vi为i方向平均孔隙速度。
5.如权利要求1所述的孔隙地下水污染物三维运移模拟方法,其特征在于,在步骤S4中,通过粒子位置坐标与计算网格坐标进行比较,确定粒子是否位于计算网格之内;边界控制的实现是通过计算粒子在该时间步长内的位移轨迹与计算网格边界的交点,将交点坐标作为粒子的新坐标。
6.如权利要求1所述的孔隙地下水污染物三维运移模拟方法,其特征在于,在步骤S6中,特定区域内粒子数量与被释放粒子总量的比值,即为该区域内地下水污染物的相对浓度。
7.如权利要求1所述的孔隙地下水污染物三维运移模拟方法,其特征在于,以一定数量的所述粒子代表地下水中污染物的一定浓度,以所述粒子的位置代表地下水中污染物的空间分布。
8.如权利要求1所述的孔隙地下水污染物三维运移模拟方法,其特征在于,所述粒子没有体积和质量。
9.如权利要求1所述的孔隙地下水污染物三维运移模拟方法,其特征在于,所述粒子能够按组分类,属于不同组别的所述粒子表示不同类型的地下水污染物。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510947333.XA CN105512417B (zh) | 2015-12-17 | 2015-12-17 | 基于粒子追踪的孔隙地下水污染物三维运移模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510947333.XA CN105512417B (zh) | 2015-12-17 | 2015-12-17 | 基于粒子追踪的孔隙地下水污染物三维运移模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105512417A true CN105512417A (zh) | 2016-04-20 |
CN105512417B CN105512417B (zh) | 2018-09-18 |
Family
ID=55720397
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510947333.XA Active CN105512417B (zh) | 2015-12-17 | 2015-12-17 | 基于粒子追踪的孔隙地下水污染物三维运移模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105512417B (zh) |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105928836A (zh) * | 2016-04-26 | 2016-09-07 | 中山大学 | 基于3d打印及spt技术的岩层液体扩散系数测量方法及装置 |
CN106529078A (zh) * | 2016-11-28 | 2017-03-22 | 西安天圆光电科技有限公司 | 一种适用于实时仿真的快速面源红外诱饵的粒子建模方法 |
CN106886682A (zh) * | 2017-01-04 | 2017-06-23 | 中国环境科学研究院 | 用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法 |
CN108595842A (zh) * | 2018-04-25 | 2018-09-28 | 中国科学院合肥物质科学研究院 | 模拟粒子辐照损伤的时间步长优化方法 |
CN111351731A (zh) * | 2020-04-16 | 2020-06-30 | 中煤科工集团重庆研究院有限公司 | 矿山工作面粉尘危害在线监测系统及方法 |
CN113298419A (zh) * | 2021-06-15 | 2021-08-24 | 辽宁工程技术大学 | 典型条件地下水污染羽空间分布插值技术的优选方法 |
CN115983160A (zh) * | 2023-02-09 | 2023-04-18 | 溧阳气动创新研究院有限公司 | 一种基于拉格朗日法的水滴轨迹模拟计算方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102073796A (zh) * | 2011-02-21 | 2011-05-25 | 南京大学 | 一种模拟溶质三维运移过程的格子行走方法 |
US20110275878A1 (en) * | 2008-09-01 | 2011-11-10 | Rainer Meckenstock | Method for the Degradation of Pollutants in Water and/or Soil |
-
2015
- 2015-12-17 CN CN201510947333.XA patent/CN105512417B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20110275878A1 (en) * | 2008-09-01 | 2011-11-10 | Rainer Meckenstock | Method for the Degradation of Pollutants in Water and/or Soil |
CN102073796A (zh) * | 2011-02-21 | 2011-05-25 | 南京大学 | 一种模拟溶质三维运移过程的格子行走方法 |
Non-Patent Citations (1)
Title |
---|
魏恒: "地下水溶质迁移模拟研究进展", 《冰川冻土》 * |
Cited By (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105928836A (zh) * | 2016-04-26 | 2016-09-07 | 中山大学 | 基于3d打印及spt技术的岩层液体扩散系数测量方法及装置 |
CN106529078A (zh) * | 2016-11-28 | 2017-03-22 | 西安天圆光电科技有限公司 | 一种适用于实时仿真的快速面源红外诱饵的粒子建模方法 |
CN106529078B (zh) * | 2016-11-28 | 2019-04-09 | 西安天圆光电科技有限公司 | 一种适用于实时仿真的快速面源红外诱饵的粒子建模方法 |
CN106886682A (zh) * | 2017-01-04 | 2017-06-23 | 中国环境科学研究院 | 用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法 |
CN106886682B (zh) * | 2017-01-04 | 2020-01-07 | 中国环境科学研究院 | 用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法 |
CN108595842A (zh) * | 2018-04-25 | 2018-09-28 | 中国科学院合肥物质科学研究院 | 模拟粒子辐照损伤的时间步长优化方法 |
CN108595842B (zh) * | 2018-04-25 | 2021-12-21 | 中国科学院合肥物质科学研究院 | 模拟粒子辐照损伤的时间步长优化方法 |
CN111351731A (zh) * | 2020-04-16 | 2020-06-30 | 中煤科工集团重庆研究院有限公司 | 矿山工作面粉尘危害在线监测系统及方法 |
CN113298419A (zh) * | 2021-06-15 | 2021-08-24 | 辽宁工程技术大学 | 典型条件地下水污染羽空间分布插值技术的优选方法 |
CN115983160A (zh) * | 2023-02-09 | 2023-04-18 | 溧阳气动创新研究院有限公司 | 一种基于拉格朗日法的水滴轨迹模拟计算方法 |
Also Published As
Publication number | Publication date |
---|---|
CN105512417B (zh) | 2018-09-18 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105512417A (zh) | 基于粒子追踪的孔隙地下水污染物三维运移模拟方法 | |
Juez et al. | A 2D weakly-coupled and efficient numerical model for transient shallow flow and movable bed | |
Beyers et al. | Modeling transient snowdrift development around complex three-dimensional structures | |
Aliabadi et al. | Effects of roof-edge roughness on air temperature and pollutant concentration in urban canyons | |
Hanna et al. | Comparisons of JU2003 observations with four diagnostic urban wind flow and Lagrangian particle dispersion models | |
CN108021780B (zh) | 一种基于无规则非结构网格模型的山洪动态仿真方法 | |
CN109614638B (zh) | 一种非直接建模的城市风环境cfd模拟方法 | |
CN103234883A (zh) | 一种基于道路交通流量实时估算中心城区pm2.5浓度的方法 | |
CN104008229A (zh) | 一种街区污染物浓度分布模型建立方法 | |
Katsuki et al. | Cellular model for sand dunes with saltation, avalanche and strong erosion: collisional simulation of barchans | |
Kenjereš et al. | Capturing transient effects in turbulent flows over complex urban areas with passive pollutants | |
CN106886682A (zh) | 用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法 | |
Assimakopoulos et al. | Experimental validation of a computational fluid dynamics code to predict the wind speed in street canyons for passive cooling purposes | |
Liu et al. | Multi-block lattice Boltzmann simulations of solute transport in shallow water flows | |
Hong et al. | Measurement and prediction of soil erosion in dry field using portable wind erosion tunnel | |
CN110990926B (zh) | 一种基于面积修正率的城市地表建筑水动力学仿真方法 | |
Liu et al. | Lattice Boltzmann method for the age concentration equation in shallow water | |
CN102073796B (zh) | 一种模拟溶质三维运移过程的格子行走方法 | |
CN113298419A (zh) | 典型条件地下水污染羽空间分布插值技术的优选方法 | |
CN105808812A (zh) | 一种地表水水龄二维介观数值模拟方法 | |
CN105160088A (zh) | 地下水流量的计算方法 | |
CN105740588B (zh) | 用于矿井下多巷道任意角度耦合的水害漫延方法 | |
CN101866392B (zh) | 一种模拟溶质一维运移过程的格子行走方法 | |
Li | Tracer mixing at fracture intersections | |
CN104156572B (zh) | 保守性污染物对流扩散数值模拟方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |