基于粒子追蹤的孔隙地下水污染物三維運(yùn)移模擬方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及地下水?dāng)?shù)值模擬計(jì)算領(lǐng)域，特別是涉及一種基于粒子追蹤的孔隙地下 水污染物三維運(yùn)移模擬方法。
【背景技術(shù)】
[0002] 孔隙地下水作為北方重要的飲用水源，其水質(zhì)直接關(guān)系到飲水安全。地下水埋藏 在地下，其污染具有很強(qiáng)的復(fù)雜性和隱蔽性，對(duì)地下水污染的掌握需要借助數(shù)值模擬方法， 通過(guò)模型的高效運(yùn)算，可以節(jié)約大量經(jīng)濟(jì)和時(shí)間成本，使管理決策人員迅速掌握地下水中 污染物濃度變化情況。
[0003] 地下水?dāng)?shù)值模擬主要包括地下水水流模擬和溶質(zhì)運(yùn)移模擬，前者數(shù)值求解地下水 水流方程，后者數(shù)值求解對(duì)流-彌散方程。描述孔隙地下水溶質(zhì)運(yùn)移的對(duì)流-彌散方程可表 述為：
[0004]
[0005] 其中C表示地下水中溶質(zhì)的濃度，v表示地下水流動(dòng)孔隙速度，D表示水力彌散系 數(shù)，t表示時(shí)間，▽為微分算子。
[0006] 孔隙地下水溶質(zhì)運(yùn)移模擬求解方法常用歐拉法、拉格朗日法，以及二者的結(jié)合。歐 拉法以空間中固定坐標(biāo)系作為參照系，常見(jiàn)的有限差分法和有限單元法屬于歐拉法。應(yīng)用 有限元和有限差分方法進(jìn)行地下水溶質(zhì)迀移模擬，有兩個(gè)固有的缺陷，一是當(dāng)網(wǎng)格Peclet 數(shù)較大，即對(duì)流項(xiàng)強(qiáng)于彌散項(xiàng)時(shí)，容易受到數(shù)值彌散的影響；二是當(dāng)模擬區(qū)域范圍較大時(shí)， 因計(jì)算網(wǎng)格較多，計(jì)算成本比較昂貴。
[0007] 拉格朗日法通過(guò)質(zhì)點(diǎn)追蹤，以運(yùn)動(dòng)坐標(biāo)系作為參照系，常見(jiàn)的拉格朗日法有粒子 追蹤法等。運(yùn)動(dòng)坐標(biāo)系容易造成數(shù)值不穩(wěn)定；另外在質(zhì)點(diǎn)追蹤過(guò)程中對(duì)質(zhì)點(diǎn)速度的連續(xù)性 有較高要求，否則由于速度插值也容易造成局部的質(zhì)量不守恒。
[0008] 歐拉-拉格朗日混合法用拉格朗日法解決溶質(zhì)運(yùn)移中的對(duì)流問(wèn)題，用歐拉法解決 彌散問(wèn)題，結(jié)合了二者的優(yōu)點(diǎn)，卻也同時(shí)存在缺點(diǎn)，且計(jì)算耗時(shí)。

【發(fā)明內(nèi)容】

[0009] 有鑒于此，本發(fā)明的目的在于提出一種孔隙地下水中污染物三維運(yùn)移模擬方法， 該數(shù)值模擬方法是基于拉格朗日法中的粒子追蹤方法，以解決模擬過(guò)程中的數(shù)值彌散問(wèn) 題，模擬得到的數(shù)值解可以完全擬合對(duì)流-彌散方程的解析解，且與有限元法相比可大大節(jié) 約計(jì)算成本。
[0010] 為實(shí)現(xiàn)本發(fā)明的上述目的，本發(fā)明提出了一種基于粒子追蹤的孔隙地下水污染物 三維運(yùn)移模擬方法，包括以下步驟：
[0011] 步驟S1:設(shè)定粒子總數(shù)、每個(gè)粒子的初始坐標(biāo)、釋放時(shí)間，設(shè)定計(jì)算時(shí)間步長(zhǎng)；
[0012] 步驟S2:在任意時(shí)刻、任意位置釋放粒子；
[0013] 步驟S3:對(duì)每個(gè)被釋放的粒子，計(jì)算其在一個(gè)時(shí)間步長(zhǎng)后的位移，粒子在地下水流 場(chǎng)中的運(yùn)動(dòng)由以下方程控制：
[0014]
[0015]
[0016]
[0017] 其中x、y、z表示粒子的空間坐標(biāo)，v表示粒子的對(duì)流運(yùn)動(dòng)速度，Δ t表示時(shí)間步長(zhǎng)，D 表示水力彌散系數(shù)，Z表示介于0到1之間的隨機(jī)數(shù)，這樣就根據(jù)粒子在t時(shí)刻的坐標(biāo)計(jì)算出 其在一個(gè)時(shí)間步長(zhǎng)后，也即t+ △ t時(shí)刻的新坐標(biāo)；
[0018] 步驟S4:根據(jù)計(jì)算得到的粒子的新坐標(biāo)，判斷該粒子是否位于計(jì)算網(wǎng)格區(qū)域內(nèi)，若 已跳出，則通過(guò)邊界控制使其回到計(jì)算區(qū)域邊界上；
[0019] 步驟S5:重復(fù)步驟S2-S4,以此迭代，計(jì)算得出每個(gè)時(shí)間節(jié)點(diǎn)上每個(gè)粒子的坐標(biāo)；
[0020] 步驟S6:通過(guò)統(tǒng)計(jì)特定區(qū)域內(nèi)粒子數(shù)量，即可得到該區(qū)域的地下水污染物濃度分 布。
[0021] 基于上述技術(shù)方案可知，本發(fā)明的孔隙地下水中污染物三維運(yùn)移模擬方法，將地 下水中的溶質(zhì)抽象為大量粒子，以粒子的運(yùn)動(dòng)來(lái)模擬溶質(zhì)在孔隙介質(zhì)中的運(yùn)移，其中以粒 子的有序運(yùn)動(dòng)來(lái)刻畫(huà)溶質(zhì)因?qū)α饕鸬霓|移過(guò)程，以粒子的隨機(jī)位移來(lái)刻畫(huà)溶質(zhì)的彌散過(guò) 程，因此可以精細(xì)刻畫(huà)溶質(zhì)在地下水中的運(yùn)移過(guò)程。傳統(tǒng)的有限元法在求解對(duì)流-彌散方程 時(shí)，需要求解的代數(shù)方程個(gè)數(shù)與計(jì)算網(wǎng)格節(jié)點(diǎn)數(shù)成正比，本發(fā)明的方法基于統(tǒng)計(jì)物理學(xué)中 的隨機(jī)行走粒子追蹤，不直接求解對(duì)流-彌散方程，而是求解粒子的位移，計(jì)算量?jī)H與粒子 個(gè)數(shù)有關(guān)，因此在三維情景下與有限元法相比可大大降低計(jì)算成本，且運(yùn)算結(jié)果與對(duì)流-彌 散方程的解析解完全擬合。本發(fā)明的方法，粒子不會(huì)憑空出現(xiàn)或消失，粒子總數(shù)量在輸入和 輸出及整個(gè)過(guò)程中是不變的，因此該方法從本質(zhì)上是質(zhì)量守恒的，可消除計(jì)算過(guò)程中的數(shù) 值彌散。
【附圖說(shuō)明】
[0022] 圖1是瞬時(shí)釋放污染物在三維均質(zhì)孔隙介質(zhì)中的污染羽分布；
[0023] 圖2是粒子到達(dá)計(jì)算網(wǎng)格區(qū)域邊界時(shí)的邊界控制效果示意圖；
[0024] 圖3是瞬時(shí)釋放污染物在三維均質(zhì)孔隙介質(zhì)中的突破曲線。
【具體實(shí)施方式】
[0025] 為使本發(fā)明的目的、技術(shù)方案和優(yōu)點(diǎn)更加清楚明白，以下結(jié)合具體實(shí)施例，并參照 附圖，對(duì)本發(fā)明作進(jìn)一步的詳細(xì)說(shuō)明。
[0026]本發(fā)明公開(kāi)了一種基于粒子追蹤的孔隙地下水污染物三維運(yùn)移模擬方法，包括以 下步驟：
[0027]步驟S1:設(shè)定粒子總數(shù)、每個(gè)粒子的初始坐標(biāo)、釋放時(shí)間，設(shè)定計(jì)算時(shí)間步長(zhǎng)；
[0028] 步驟S2:在任意時(shí)刻、任意位置釋放粒子；
[0029] 步驟S3:對(duì)每個(gè)被釋放的粒子，計(jì)算其在一個(gè)時(shí)間步長(zhǎng)后的位移，粒子在地下水流 場(chǎng)中的運(yùn)動(dòng)由以下方程控制：
[0033] 其中x、y、z表示粒子的空間坐標(biāo)，ν表示粒子的對(duì)流運(yùn)動(dòng)速度，Δ t表示時(shí)間步長(zhǎng)，D 表示水力彌散系數(shù)，Z表示介于0到1之間的隨機(jī)數(shù)，這樣就根據(jù)粒子在t時(shí)刻的坐標(biāo)計(jì)算出 其在一個(gè)時(shí)間步長(zhǎng)后，也即t+ △ t時(shí)刻的新坐標(biāo)；
[0034]步驟S4:根據(jù)計(jì)算得到的粒子的新坐標(biāo)，判斷粒子是否位于計(jì)算網(wǎng)格區(qū)域內(nèi)，若已 跳出，則通過(guò)邊界控制使其回到計(jì)算區(qū)域邊界上；
[0035] 步驟S5:重復(fù)步驟S2-S4，以此迭代，計(jì)算得出每個(gè)時(shí)間節(jié)點(diǎn)上每個(gè)粒子的坐標(biāo)；
[0036] 步驟S6:通過(guò)統(tǒng)計(jì)特定區(qū)域內(nèi)粒子數(shù)量，即可得到該區(qū)域的地下水污染物濃度分 布。
[0037] 作為優(yōu)選，在步驟S1中，設(shè)定粒子初始坐標(biāo)時(shí)，允許多個(gè)粒子的初始坐標(biāo)完全相 同；粒子釋放時(shí)間表示該粒子初次被釋放的時(shí)間，應(yīng)為時(shí)間步長(zhǎng)的整數(shù)倍;時(shí)間步長(zhǎng)需滿足 庫(kù)朗數(shù)(Courant number)小于 1。
[0038] 作為優(yōu)選，在步驟S2中，對(duì)任意一個(gè)粒子，若當(dāng)前步的計(jì)算時(shí)間大于或等于步驟S1 中設(shè)定的粒子釋放時(shí)間，則該粒子被釋放;若在某一坐標(biāo)位置所有粒子的釋放時(shí)間相同，則 該釋放過(guò)程是瞬時(shí)的，若粒子釋放時(shí)間不同且構(gòu)成一個(gè)連續(xù)時(shí)段，則該釋放過(guò)程是持續(xù)的。
[0039] 作為優(yōu)選，在步驟S3中，粒子的對(duì)流運(yùn)動(dòng)速度ν由地下水流場(chǎng)計(jì)算所得，并通過(guò)粒 子所在單元節(jié)點(diǎn)插值獲得;Z由計(jì)算程序中的