本發(fā)明屬于流體動(dòng)力學(xué)數(shù)值模擬與計(jì)算領(lǐng)域，尤其涉及一種基于交替tvd差分算法的地球流體物質(zhì)平流輸運(yùn)模擬方法。

背景技術(shù)：

數(shù)值模式是研究和預(yù)報(bào)自然界中河流、海洋等運(yùn)動(dòng)以及其內(nèi)物質(zhì)輸運(yùn)的重要工具，隨著計(jì)算機(jī)技術(shù)和社會(huì)經(jīng)濟(jì)的發(fā)展，數(shù)值模式在科學(xué)研究、業(yè)務(wù)預(yù)報(bào)以及湖庫(kù)、海洋等環(huán)境管理與規(guī)劃等行業(yè)中發(fā)揮著日益重要的作用。

平流是地球流體中攜載物質(zhì)進(jìn)行空間輸送的主要過(guò)程，平流過(guò)程可表示為方程(1)。數(shù)值模式中，通常采用差分的方法將(1)式進(jìn)行離散，得到如式(2)形式的平流差分方程，然后對(duì)式(2)進(jìn)行數(shù)值迭代，得到物質(zhì)濃度c的空間分布及隨時(shí)間變化。

其中，c為某輸運(yùn)物質(zhì)濃度，u、v、w分別為x、y、z方向的流速。

其中，為半層網(wǎng)格處的物質(zhì)濃度，δt為離散的時(shí)間步長(zhǎng)，δx、δy、δz分別為x、y、z方向的網(wǎng)格間距，標(biāo)號(hào)上標(biāo)n表示時(shí)間步數(shù)，下標(biāo)i、j、k分別表示x、y、z方向的網(wǎng)格標(biāo)號(hào)。

在差分方程(2)中，如何計(jì)算半層網(wǎng)格處的物質(zhì)濃度是不同差分格式的主要差別所在。例如，迎風(fēng)格式使用的是上游網(wǎng)格濃度，中差格式使用的是上游和下游網(wǎng)格的平均濃度。這些算法雖然計(jì)算簡(jiǎn)便，但存在著精度低或強(qiáng)數(shù)值頻散的問(wèn)題，導(dǎo)致模擬結(jié)果失真較強(qiáng)且會(huì)影響模型的穩(wěn)定性。

總變差消減(totalvariationdiminishing，簡(jiǎn)稱(chēng)tvd)差分格式是一類(lèi)高精度無(wú)頻散的平流差分算法，其被廣泛應(yīng)用在海洋和大氣數(shù)值模式及相關(guān)研究中。tvd算法主要是在原有的低精度迎風(fēng)格式的基礎(chǔ)上加入反擴(kuò)散的高精度項(xiàng)，提高平流格式的精度，同時(shí)引入限制函數(shù)來(lái)控制加入的高精度項(xiàng)的強(qiáng)度，以此，保持整個(gè)向前積分過(guò)程中，數(shù)值格式的總變差保持不增，從而有效抑制高精度格式通常具有的數(shù)值頻散性。

tvd平流格式的算法如下(以x方向?yàn)槔?：

如式(3)所示，在計(jì)算時(shí)，在迎風(fēng)格式的基礎(chǔ)上加上高精度修正項(xiàng)，并在修正項(xiàng)前加入限制函數(shù)ψ(ri+1/2)，由此得到：

其中，

加入限制函數(shù)后保證計(jì)算結(jié)果的總變差保持不增，即：

tv(cⁿ⁺¹)≤tv(cⁿ)(4)

需使得限制函數(shù)滿(mǎn)足下述條件：

其中，r表示濃度梯度比，即或

雖然tvd限制函數(shù)有多種，但單一的使用一種限制函數(shù)不能很好的控制tvd格式中加入的高精度項(xiàng)的強(qiáng)度，由此造成了以上tvd數(shù)值格式仍具有一定的數(shù)值擴(kuò)散(如，minmod、vanleer、muscl、hsimt)和反擴(kuò)散(如，superbee)偽效應(yīng)。另外，之前tvd格式可以保證一維情況下的總變差不增，而在多維情況下則差強(qiáng)人意，仍然存在數(shù)值頻散。因此，利用現(xiàn)有tvd平流差分格式對(duì)地球流體中的物質(zhì)輸運(yùn)進(jìn)行長(zhǎng)時(shí)間模擬時(shí)，仍存在不能令人滿(mǎn)意的模擬結(jié)果。

技術(shù)實(shí)現(xiàn)要素：

為了解決現(xiàn)有技術(shù)的不足，本發(fā)明提供了一種基于交替tvd差分算法的地球流體物質(zhì)平流輸運(yùn)模擬方法，其提供一種有效減弱數(shù)值擴(kuò)散和反擴(kuò)散偽效應(yīng)并解決多維頻散平流格式的差分方法-tvdal(tvdwithalternatinglimiters)，將其應(yīng)用于地球流體中對(duì)物質(zhì)輸運(yùn)的模擬和預(yù)報(bào)。

本發(fā)明的一種基于交替tvd差分算法的地球流體物質(zhì)平流輸運(yùn)模擬方法，包括：

構(gòu)建地球流體物質(zhì)的平流輸送方程，進(jìn)而對(duì)所述平流輸送方程進(jìn)行數(shù)值離散，得到地球流體物質(zhì)輸送的離散方程；所述地球流體物質(zhì)輸送的離散方程的維數(shù)至少為二維；

將地球流體物質(zhì)輸送的離散方程拆解成與其維數(shù)相等的步數(shù)；

獲取地球流體物質(zhì)的流場(chǎng)數(shù)據(jù)并輸入至地球流體物質(zhì)輸送的離散方程中進(jìn)行數(shù)值迭代計(jì)算，預(yù)測(cè)出各個(gè)預(yù)設(shè)迭代時(shí)間段內(nèi)的地球流體物質(zhì)濃度；

在數(shù)值迭代計(jì)算的過(guò)程中，在不同的預(yù)設(shè)迭代時(shí)間段內(nèi)分別在每個(gè)維度上按照與其維數(shù)相等的步數(shù)交替使用不同類(lèi)型的總變差消減限制函數(shù)來(lái)計(jì)算下一個(gè)時(shí)間段內(nèi)地球流體物質(zhì)濃度，最終模擬和預(yù)測(cè)出地球流體物質(zhì)平流輸運(yùn)過(guò)程中流體物質(zhì)濃度的變化。

進(jìn)一步的，基于變差消減格式的構(gòu)造方法，對(duì)地球流體物質(zhì)的平流輸送方程進(jìn)行數(shù)值離散。

進(jìn)一步的，總變差消減限制函數(shù)包括數(shù)值擴(kuò)散性質(zhì)的限制函數(shù)和數(shù)值反擴(kuò)散性質(zhì)的限制函數(shù)。

進(jìn)一步的，地球流體物質(zhì)濃度是關(guān)于時(shí)間和空間的函數(shù)。

進(jìn)一步的，求取地球流體物質(zhì)的平流輸送方程的過(guò)程為：

地球流體物質(zhì)濃度分別對(duì)表征空間的方向參數(shù)的導(dǎo)函數(shù)，再與相應(yīng)表征空間的方向參數(shù)的流速相乘后，最后累加地球流體物質(zhì)濃度對(duì)時(shí)間的導(dǎo)函數(shù)，得到地球流體物質(zhì)的平流輸送方程。

其中，所述數(shù)值擴(kuò)散性質(zhì)的限制函數(shù)為minmod函數(shù)。

所述數(shù)值擴(kuò)散性質(zhì)的限制函數(shù)為vanleer函數(shù)。

所述數(shù)值擴(kuò)散性質(zhì)的限制函數(shù)為muscl函數(shù)。

所述數(shù)值擴(kuò)散性質(zhì)的限制函數(shù)為hsimt函數(shù)。

所述數(shù)值反擴(kuò)散性質(zhì)的限制函數(shù)包括superbee函數(shù)。

其中，在實(shí)際計(jì)算中，常采用的幾種經(jīng)典的tvd限制函數(shù)主要有：

superbee：

ψ(r)^superbee＝max[0,min(2r,1),min(r,2)](6)

vanleer(orharmonic)：

ψ(r)^vanleer＝(r+|r|)/(r+1)(7)

muscl：

minmod：

ψ(r)^minmod＝max[0,min(r,1)](9)

hsimt：

ψ(r)^hsimt＝max[0,min(2r,2,β)](10)

其中，

與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果是：

(1)本方法并沒(méi)有更改原有算式增加額外計(jì)算量，而是在不同時(shí)間步交替使用不同性質(zhì)的算法，很大程度地減弱了傳統(tǒng)tvd平流差分算法的數(shù)值擴(kuò)散和反擴(kuò)散偽效應(yīng)，使得模擬結(jié)果的誤差更小。

(2)本方法中將交替限制函數(shù)與離散方程拆分相結(jié)合，基于得到的多維無(wú)頻散平流算法，有效消除tvd平流格式在多維計(jì)算情況下物質(zhì)輸運(yùn)模擬的數(shù)值頻散問(wèn)題。

附圖說(shuō)明

構(gòu)成本申請(qǐng)的一部分的說(shuō)明書(shū)附圖用來(lái)提供對(duì)本申請(qǐng)的進(jìn)一步理解，本申請(qǐng)的示意性實(shí)施例及其說(shuō)明用于解釋本申請(qǐng)，并不構(gòu)成對(duì)本申請(qǐng)的不當(dāng)限定。

圖1是本發(fā)明的流程示意圖；

圖2(a)是初始時(shí)刻物質(zhì)濃度分布圖；

圖2(b)是基于傳統(tǒng)tvd算法的superbee模擬結(jié)果；

圖2(c)是基于傳統(tǒng)tvd算法的hsimt模擬結(jié)果；

圖2(d)是傳統(tǒng)tvd算法結(jié)合離散方程拆分算法的superbee模擬結(jié)果；

圖2(e)是傳統(tǒng)tvd算法結(jié)合離散方程拆分算法的hsimt模擬結(jié)果；

圖2(f)是本發(fā)明方法的模擬結(jié)果。

具體實(shí)施方式

應(yīng)該指出，以下詳細(xì)說(shuō)明都是例示性的，旨在對(duì)本申請(qǐng)?zhí)峁┻M(jìn)一步的說(shuō)明。除非另有指明，本文使用的所有技術(shù)和科學(xué)術(shù)語(yǔ)具有與本申請(qǐng)所屬技術(shù)領(lǐng)域的普通技術(shù)人員通常理解的相同含義。

需要注意的是，這里所使用的術(shù)語(yǔ)僅是為了描述具體實(shí)施方式，而非意圖限制根據(jù)本申請(qǐng)的示例性實(shí)施方式。如在這里所使用的，除非上下文另外明確指出，否則單數(shù)形式也意圖包括復(fù)數(shù)形式，此外，還應(yīng)當(dāng)理解的是，當(dāng)在本說(shuō)明書(shū)中使用術(shù)語(yǔ)“包含”和/或“包括”時(shí)，其指明存在特征、步驟、操作、器件、組件和/或它們的組合。

地球流體物質(zhì)輸送的離散方程的維數(shù)至少為二維，下面以地球流體物質(zhì)輸送的離散方程的維數(shù)是二維為例：

如圖1所示，本發(fā)明的基于交替tvd差分算法的地球流體物質(zhì)平流輸運(yùn)模擬方法，包括：

步驟1：構(gòu)建地球流體物質(zhì)的平流輸送方程，進(jìn)而對(duì)所述平流輸送方程進(jìn)行數(shù)值離散，得到地球流體物質(zhì)輸送的離散方程。

構(gòu)建地球流體物質(zhì)的平流輸送方程為公式(1)：

對(duì)所述平流輸送方程進(jìn)行數(shù)值離散，得到地球流體物質(zhì)輸送的離散方程為公式(2)：

其中，為半層網(wǎng)格處的物質(zhì)濃度，δt為離散的時(shí)間步長(zhǎng)，δx、δy分別為x、y方向的網(wǎng)格間距，標(biāo)號(hào)上標(biāo)n表示時(shí)間步數(shù)，下標(biāo)i、j分別表示x、y方向的網(wǎng)格標(biāo)號(hào)。

步驟2：將地球流體物質(zhì)輸送的離散方程拆解成與其維數(shù)相等的步數(shù)。

步驟3：獲取地球流體物質(zhì)的流場(chǎng)數(shù)據(jù)并輸入至地球流體物質(zhì)輸送的離散方程中進(jìn)行數(shù)值迭代計(jì)算，預(yù)測(cè)出各個(gè)預(yù)設(shè)迭代時(shí)間段內(nèi)的地球流體物質(zhì)濃度；

在數(shù)值迭代計(jì)算的過(guò)程中，在不同的預(yù)設(shè)迭代時(shí)間段內(nèi)分別在每個(gè)維度上按照與其維數(shù)相等的步數(shù)交替使用不同類(lèi)型的總變差消減限制函數(shù)來(lái)計(jì)算下一個(gè)時(shí)間段內(nèi)地球流體物質(zhì)濃度，最終模擬和預(yù)測(cè)出地球流體物質(zhì)平流輸運(yùn)過(guò)程中流體物質(zhì)濃度的變化。

具體地，輸入流速u(mài)、v信息，對(duì)離散方程進(jìn)行迭代并判斷迭代步數(shù)的奇偶性，若為奇數(shù)計(jì)算則采用下面步驟a)計(jì)算，若為偶數(shù)則采用步驟b)計(jì)算。

步驟a)：

將離散方程式(2)拆解為式(3-1)和式(3-2)兩步計(jì)算：

其中，的上標(biāo)com和dif分別表示利用數(shù)值反擴(kuò)散性質(zhì)(例如，superbee)和數(shù)值擴(kuò)散性質(zhì)(例如，minmod)的tvd限制函數(shù)來(lái)計(jì)算半層網(wǎng)格處物質(zhì)濃度

將離散方程拆解成(3-1)和(3-2)兩步，依次計(jì)算x和y兩個(gè)方向上平流過(guò)程：第一步計(jì)算x方向上平流過(guò)程時(shí)，先使用一種數(shù)值反擴(kuò)散性質(zhì)的tvd限制函數(shù)，計(jì)算得到中間變量c*；第二步將中間變量c*代入y方向差分方程，計(jì)算y方向上平流過(guò)程，使用一種數(shù)值擴(kuò)散性質(zhì)的tvd限制函數(shù)計(jì)算得到新時(shí)間步的濃度cⁿ⁺¹。

步驟b)：

同樣分(4-1)和(4-2)兩步計(jì)算，但第一步先計(jì)算y方向平流過(guò)程，使用一種數(shù)值反擴(kuò)散性質(zhì)的tvd限制函數(shù)，計(jì)算得到中間變量c*；第二步將中間變量c*代入x方向差分方程，計(jì)算x方程平流過(guò)程，使用一種數(shù)值擴(kuò)散性質(zhì)的tvd限制函數(shù)，計(jì)算得到新時(shí)間步的物質(zhì)濃度。

其中，的上標(biāo)com和dif分別表示利用數(shù)值反擴(kuò)散性質(zhì)(例如，superbee)和數(shù)值擴(kuò)散性質(zhì)(例如，minmod)的tvd限制函數(shù)來(lái)計(jì)算半層網(wǎng)格處物質(zhì)濃度

最后判斷積分時(shí)長(zhǎng)是否達(dá)到模式的設(shè)定值，沒(méi)有則返回步驟2繼續(xù)，若達(dá)到則終止計(jì)算。

本發(fā)明的內(nèi)容可以通過(guò)以下的模擬結(jié)果進(jìn)一步進(jìn)行說(shuō)明。

1、模擬實(shí)驗(yàn)：應(yīng)用本發(fā)明方法與其他方法，對(duì)給定某一濃度分布的物質(zhì)平流過(guò)程進(jìn)行模擬。

2、模擬結(jié)果

在一個(gè)水平二維海域，其內(nèi)存在以12小時(shí)為周期沿45°方向的往復(fù)流動(dòng)，初始時(shí)刻物質(zhì)濃度分布給定如圖2(a)所示分布，使用不同平流差分算法得到的100個(gè)周期后模擬結(jié)果分別如圖2(b)-圖2(f)所示。其中，圖2(b)和2(c)為基于傳統(tǒng)tvd算法superbee和hsimt的模擬結(jié)果，圖2(d)與圖2(e)為傳統(tǒng)tvd算法結(jié)合離散方程拆分算法的superbee和hsimt模擬結(jié)果；圖2(f)為本發(fā)明方法的模擬結(jié)果。

從圖2(a)-圖2(f)可以看出，基于傳統(tǒng)的tvd算法的模擬結(jié)果存在明顯的數(shù)值頻散，使用離散方程拆分算法后數(shù)值頻散得到改善，但仍存在數(shù)值反擴(kuò)散和擴(kuò)散的偽效應(yīng)，使用本發(fā)明方法則基本保持了真解的分布，減弱了傳統(tǒng)tvd算法的數(shù)值反擴(kuò)散和擴(kuò)散偽效應(yīng)，同時(shí)計(jì)算結(jié)果無(wú)數(shù)值頻散，取得了令人滿(mǎn)意的模擬結(jié)果。

綜上，本發(fā)明提出的該方法有效地減弱數(shù)值擴(kuò)散和反擴(kuò)散偽效應(yīng)并解決多維頻散平流格式的差分算法-tvdal的物質(zhì)平流模擬方法。本發(fā)明的方法通過(guò)交替使用性質(zhì)相反偽效應(yīng)的傳統(tǒng)tvd限制函數(shù)，有效的減弱了數(shù)值擴(kuò)散和反擴(kuò)散，減小了平流模擬結(jié)果的誤差，同時(shí)，將多維離散方程進(jìn)行拆分計(jì)算，保證了多維平流模擬是無(wú)數(shù)值頻散產(chǎn)生，以上改進(jìn)通過(guò)物質(zhì)平流實(shí)驗(yàn)的模擬結(jié)果可以看出。

上述雖然結(jié)合附圖對(duì)本發(fā)明的具體實(shí)施方式進(jìn)行了描述，但并非對(duì)本發(fā)明保護(hù)范圍的限制，所屬領(lǐng)域技術(shù)人員應(yīng)該明白，在本發(fā)明的技術(shù)方案的基礎(chǔ)上，本領(lǐng)域技術(shù)人員不需要付出創(chuàng)造性勞動(dòng)即可做出的各種修改或變形仍在本發(fā)明的保護(hù)范圍以?xún)?nèi)。