本公開涉及大氣科學(xué)領(lǐng)域，尤其涉及一種模擬復(fù)雜流體運動的高精度求解方法。

背景技術(shù)：

數(shù)值模式和數(shù)值模擬是定量研究復(fù)雜流體運動如天氣和氣候變化的主要工具，但是由于數(shù)學(xué)模型中各變量之間復(fù)雜的非線性作用關(guān)系，導(dǎo)致模擬結(jié)果中存在著不確定性。其中計算誤差對數(shù)值模擬結(jié)果的影響可以從復(fù)雜的大氣環(huán)流模式、耦合模式運行結(jié)果看出，也可以從簡單的混沌動力系統(tǒng)、準(zhǔn)地轉(zhuǎn)模式的數(shù)值試驗得到驗證。因此，如何采取有效的方法來控制計算誤差的增長，對天氣過程和復(fù)雜流體運動的長時間的數(shù)值計算和精確的數(shù)值模擬至關(guān)重要。

許多流動運動(從復(fù)雜的n-s方程到簡單的平流方程)都可以寫為算子形式：為了使用計算機求解此類方程，目前已有一些研究引進了高精度算法來對其進行數(shù)值模擬。但這些研究在空間差分上使用高精度格式的較多，而時間積分方案多為較低精度的計算格式(習(xí)慣上將階數(shù)3階或以下的算法稱為低精度或低階算法；高精度一般為4-9階算法，也稱高階算法；10階以上稱為超高階或超高精度算法)。這就導(dǎo)致了空間、時間計算精度有時不匹配，會影響模式長時間計算的準(zhǔn)確性。只有整體提高算法的精度(即同時提高時間和空間精度)，才可使數(shù)值模擬中的總誤差得到控制，即需要研究適用于大氣科學(xué)數(shù)值模式的時/空精度匹配的高階算法，以減小天氣和氣候數(shù)值模式中由于計算誤差導(dǎo)致的不確定性。

公開內(nèi)容

(一)要解決的技術(shù)問題

本公開提供了一種模擬復(fù)雜流體運動的高精度taylor-li求解方法，以至少部分解決以上所提出的技術(shù)問題。

(二)技術(shù)方案

根據(jù)本公開的一個方面，提供了一種模擬復(fù)雜流體運動的高精度求解方法，針對流體運動方程f為關(guān)于時間和空間的變量，l為包含f以及f關(guān)于空間變量導(dǎo)數(shù)的算子，先計算高精度的空間導(dǎo)數(shù)，進而計算出高精度的時間導(dǎo)數(shù)，再使用taylor格式求出下一步計算解；包括：

步驟a，給定初始條件，旋轉(zhuǎn)合適的空間步長h和時間步長τ，設(shè)定積分的終值時刻，最高時間積分精度為m階；

步驟b，在選定的m階時間精度下，采用n階高階精度空間差分格式進行空間差分計算，n與最高時間積分精度m的值相匹配；

步驟c，設(shè)置時間積分精度為k階，k＝1，2，…m，利用高精度的空間導(dǎo)數(shù)，進而計算出高精度的k階時間導(dǎo)數(shù)；

步驟d，增加時間積分階數(shù)k，重復(fù)步驟c，直到k＝m循環(huán)結(jié)束，采用taylor級數(shù)法進行時間積分計算一步時間積分；

步驟e，重復(fù)上述步驟b～d，積分得到指定目標(biāo)時刻表示流體速度方向的數(shù)值解。

在本公開的一些實施例中，所述步驟a中，最高時間積分精度m為大于等于3的整數(shù)；空間差分階數(shù)n為大于等于6的整數(shù)，例如最高時間積分精度m為5，10或20；空間差分階數(shù)n為30，50，100或500。當(dāng)最高時間積分精度m越大時，設(shè)置匹配的空間差分精度n階越高。

在本公開的一些實施例中，給定的所述初始條件為連續(xù)光滑的，周期性的初始條件。

在本公開的一些實施例中，所述高精度求解方法使用了multipleprecision(mp)庫，采用1024二進制位精度。

在本公開的一些實施例中，所述步驟b中，空間差分計算采用遞歸微分的方式計算，空間差分方法采用以下公式：

其中u是一個有關(guān)x，t的函數(shù)u＝u(x,t)，“ux⁽¹⁾(xi)”即為在xi格點上的函數(shù)u對x的一階導(dǎo)數(shù)或者空間差分方法采用緊致空間差分方案。

在本公開的一些實施例中，所述步驟d中，時間積分計算采用以下公式：

其中k＝1,2,…m，

在本公開的一些實施例中，所述步驟c中，對于線性方程描述的復(fù)雜流體運動采用以下公式求得時間導(dǎo)數(shù)

其中k＝1,2,…m，為空間導(dǎo)數(shù)。

在本公開的一些實施例中，對于非線性方程描述的復(fù)雜流體運動其中a為非線性算子，采用以下公式求得時間導(dǎo)數(shù)：

其中k＝1,2,…m，為空間導(dǎo)數(shù)。

(三)有益效果

從上述技術(shù)方案可以看出，本公開的模擬復(fù)雜流體運動的高精度求解方法至少具有以下有益效果其中之一：

(1)通過靈活地調(diào)整計算精度階數(shù)，相比常規(guī)方法使用的固定低階時間積分，可以實現(xiàn)5階、10階甚至20階的時間積分算法，從而相應(yīng)的空間差分精度也能從6階拓展到30階、50階、100階甚至500階，因此該時/空精度匹配的高階算法使得算法的精度得到整體提高，數(shù)值模擬的總誤差得到控制；

(2)通過采用連續(xù)光滑、周期性較好的初值，可以獲得很好的計算效果，精度階數(shù)越高，誤差越小。對于線性和非線性問題的計算復(fù)雜流體模擬結(jié)果表明，空間差分的精度可以遠(yuǎn)超6階；

(3)使用了遞歸微分來提高計算速度，從而所實現(xiàn)的快速高精度差分格式比直接計算高階空間微分項的方案速度有顯著提升。這種計算方法不僅內(nèi)存占用小，計算速度快，而且超過20階的空間差分精度，對于常規(guī)方案已經(jīng)無法在常規(guī)時間內(nèi)完成的計算，本公開的方法仍然有效，當(dāng)精度提高到100階時仍能夠在常規(guī)時將內(nèi)完成計算；

(4)由于求解方法使用了multipleprecision(mp)庫，采用1024二進制位精度，相當(dāng)于200位以上的十進制有效數(shù)字，因此足以區(qū)分絕對誤差小到10^-200的數(shù)值解。

附圖說明

圖1為本公開第一實施例對流體運動方程的計算流程圖。

圖2為本公開第二實施例對線性方程的計算流程圖。

圖3為本公開第二實施例線性情形下計算誤差隨空間精度階數(shù)的變化。

圖4為本公開第三實施例對非線性方程的計算流程圖。

圖5為本公開第三實施例非線性情形下計算誤差隨空間精度階數(shù)的變化。

具體實施方式

為使本公開的目的、技術(shù)方案和優(yōu)點更加清楚明白，以下結(jié)合具體實施例，并參照附圖，對本公開進一步詳細(xì)說明。

本公開某些實施例于后方將參照所附附圖做更全面性地描述，其中一些但并非全部的實施例將被示出。實際上，本公開的各種實施例可以許多不同形式實現(xiàn)，而不應(yīng)被解釋為限于此數(shù)所闡述的實施例；相對地，提供這些實施例使得本公開滿足適用的法律要求。

在本公開的第一個示例性實施例中，提供了一種模擬復(fù)雜流體運動的高精度求解方法，圖1為本公開第一實施例流體方程計算方法流程圖。如圖1所示：針對流體運動方程先計算高精度的空間導(dǎo)數(shù)，進而計算出高精度的時間導(dǎo)數(shù)，再使用taylor格式求出下一步計算解。具體步驟如下：

步驟a.給定初始條件，旋轉(zhuǎn)合適的空間步長h和時間步長τ，設(shè)定積分的終值時刻，最高時間積分精度為m階；

其中：

最高時間積分精度m為大于等于3的整數(shù)，例如，m＝5,10或20；

給定的初始條件為連續(xù)光滑的，周期性好的初始條件；

步驟b.在選定的m階時間精度下，采用n階高階精度空間差分格式進行空間差分計算，n與最高時間積分精度m的值相匹配；

其中：n為大于等于6的整數(shù)，例如n＝30，50，100或500，當(dāng)最高時間積分精度m越大時，設(shè)置空間差分精度n階越高；

空間差分計算采用遞歸微分的方式計算，空間差分方法采用以下公式：

其中(i＝j(luò))，“ux⁽¹⁾(xi)”即為在xi格點上的函數(shù)u對x的一階導(dǎo)數(shù)

本公開在計算空間高階導(dǎo)數(shù)uy^(m)(yi)時盡量避免使用系數(shù)來直接計算，而是使用來重復(fù)計算m次，此即遞歸微分的方式。這樣做的好處如下：僅使用了兩個控制參數(shù)m和i，因此比使用4個參數(shù)n，m，i，j的循環(huán)要高效。當(dāng)m和n較大時，若在求解程序中計算是較費時間的，特別是當(dāng)m大于10階以上時更為明顯。而且所需的內(nèi)存為m×(n+1)×(n+1)×(n+1)，也較大，還會影響cpu的cache使用率，從而降低計算速度。此外，在實際計算時公式簡單，速度快，而且存儲一次即可反復(fù)使用，內(nèi)存占用比小1-2個量級。

步驟c.設(shè)置時間積分精度為k階，k＝1，2，…m，利用高精度的空間導(dǎo)數(shù)，進而計算出高精度的k階時間導(dǎo)數(shù)；

步驟d.增加時間積分階數(shù)k，重復(fù)步驟c，直到k＝m循環(huán)結(jié)束，采用taylor級數(shù)法進行時間積分計算一步時間積分；

步驟e.重復(fù)上述步驟b～d，積分得到指定目標(biāo)時刻的數(shù)值解。

為解決高精度求解方法在計算時還會碰到的舍入誤差問題，如果僅采用雙精度計算，會發(fā)現(xiàn)高精度算法得到的絕對誤差會徘徊在10^-15～10^-16，這恰好是雙精度計算時的相對誤差極限。所以為了得到超高精度方案的計算效果，采用多精度計算是必要的。本公開的求解方法使用了multipleprecision(mp)庫，采用1024二進制位精度，相當(dāng)于200位以上的十進制有效數(shù)字，足以區(qū)分絕對誤差小到10^-200的數(shù)值解。

至此，本公開第一實施例線性方程情況下的模擬復(fù)雜流體運動的高精度求解方法介紹完畢。

在本公開的第二個示例性實施例中，還提供了一種模擬復(fù)雜流體運動的高精度求解方法。

為了便于理解，現(xiàn)以一個復(fù)雜流體運動的典型方程(burgers方程)作為算例來進行說明。本實施例中，以線性為例進行說明。

圖2為本公開第二實施例線性方程計算方法流程圖。如圖2所示：

第一步，對線性方程(1)進行空間差分計算，

公式中u是一個有關(guān)x，t的函數(shù)u＝u(x,t)。

在本方案中，空間差分方法選用任意階精度的有限差分公式(2)。其優(yōu)點是可以實現(xiàn)高階精度(甚至超高階精度)的空間差分格式，從而方便尋找與高階時間積分匹配的精度更高的差分精度。

其中

第二步是計算時間積分，本求解方法采用的是taylor級數(shù)法。taylor級數(shù)法的核心是計算導(dǎo)數(shù)因此可以利用線性方程即對時間導(dǎo)數(shù)進行變換，得到公式(3)：

其中k為時間積分的精度，若我們的計算中最高為m階時間積分精度，那么k＝1,2,…m。計算方程(3)的關(guān)鍵是計算導(dǎo)數(shù)項這可使用公式(2)來計算，“ux⁽¹⁾(xi)”即為在xi格點上的函數(shù)u對x的一階導(dǎo)數(shù)

計算時k＝1,2,…m進行循環(huán)，當(dāng)k＝m時，循環(huán)結(jié)束。此時，可以利用公式(4)求得下一時刻的數(shù)值解。

其中

這就在格點位置為xi完成了時間步長為τ的m階的taylor格式的一步時間積分。重復(fù)上述過程，即可積分得到指定目標(biāo)時刻的數(shù)值解。

為驗證結(jié)果中誤差隨時間積分精度的變化可以采用以下步驟：

a.給定初始條件，旋轉(zhuǎn)合適的空間步長和時間步長，設(shè)定積分的終值時刻；

b.選擇時將積分精度為3階，進行計算。

c.在選定的時間精度下，改變不同的空間差分精度，由3階到30階，進行計算。

d.分析結(jié)果中誤差的隨空間精度的變化情況。

e.改變時間積分階數(shù)，重復(fù)步驟b～d，得到結(jié)果中誤差隨時間積分精度的變化。

圖3所示為第二實施例線性情形下計算誤差隨空間精度階數(shù)的變化，橫坐標(biāo)為空間精度階數(shù)，縱坐標(biāo)為誤差取對數(shù)(log10)，圖中，“□”、“●”、“*”連線分別代表時間精度為3、4、5、6階。

圖3所示求解平流方程(1)時，使用gauss波初始條件：計算時x方向的空間步長為h＝1/200，時間步長為τ＝1/400，計算區(qū)域為[0,1]，計算400步，從圖5可見，若m＝3階，n達到6階時，計算誤差與n大于6階時的計算誤差相差不多，這與feng07的結(jié)果一致；但m＝4階時，n達到8階左右計算誤差接近平緩變化；m＝5階時，n達到10階左右計算誤差接近平緩變化；m＝6階時，n達到12階左右計算誤差接近平緩變化；這些結(jié)果說明，試驗中總計算誤差受時間方向差分精度的影響較大，當(dāng)時間方向精度足夠高時，空間方向的格式精度可以遠(yuǎn)超過6階。

至此，本公開第二實施例線性方程情況下的模擬復(fù)雜流體運動的高精度求解方法介紹完畢。

在本公開的第三個示例性實施例中，提供了一種非線性方程情況下的模擬復(fù)雜流體運動的高精度求解方法。

為了便于理解，現(xiàn)以一個復(fù)雜流體運動的典型方程(burgers方程)作為算例來進行說明。本實施例中，以非線性為例進行說明。

圖4為本公開第三實施例線性方程計算方法流程圖。如圖4所示，本實施例中，對于一個描述復(fù)雜流體的非線性方程(5)，其空間差分方法也選為任意階精度的有限差分公式(2)。

其中a為非線性算子，計算時間積分時，非線性方程與線性方程有所不同，使用的是推導(dǎo)得到的公式(6)。

其中對非線性方程，仍然用空間差分的公式(2)來計算方程(6)中的使k＝1,2,…m循環(huán)計算，直到k＝m時，循環(huán)結(jié)束，即可完成對非線性方程中的一步計算。重復(fù)上述過程，即可積分得到指定目標(biāo)時刻的數(shù)值解。

為驗證結(jié)果中誤差隨時間積分精度的變化可以采用以下步驟：

a.給定初始條件，旋轉(zhuǎn)合適的空間步長和時間步長，設(shè)定積分的終值時刻；

b.選擇時將積分精度為3階，進行計算。

c.在選定的時間精度下，改變不同的空間差分精度，由3階到30階，進行計算。

d.分析結(jié)果中誤差的隨空間精度的變化情況。

e.改變時間積分階數(shù)，重復(fù)步驟b～d，得到結(jié)果中誤差隨時間積分精度的變化。

圖5所示為非線性情形下計算誤差隨空間精度階數(shù)的變化，橫坐標(biāo)為空間精度階數(shù)，縱坐標(biāo)為誤差取對數(shù)(log10)，圖中，“□”、“●”、“*”連線分別代表時間精度為3、4、5、6階。

圖5所示求解無粘性burgers方程(7)時，初始條件為：u(x,0)＝-sin(x)。試驗時，計算區(qū)域選為[-π,π]，網(wǎng)格數(shù)為n＝800，時間步長為τ＝0.001，側(cè)邊界條件為周期邊界條件。計算時進行檢驗計算誤差的時刻選為：t＝0.8，共計算800步。從圖5可以看出，當(dāng)時間精度為3階時，計算誤差在空間精度達到6階后就基本保持不變；當(dāng)時間精度為4、5、6階時，有效的空間差分精度階數(shù)分別可達11、18、33。這個試驗說明了對非線性burgers方程，只要時間積分精度階數(shù)足夠高，空間方向的差分精度階數(shù)可以超過6階。

為了達到簡要說明的目的，上述第一實施例中任何可作相同應(yīng)用的技術(shù)特征敘述皆并于此，無需再重復(fù)相同敘述。

至此，本公開第三實施例非線性方程情況下的模擬復(fù)雜流體運動的高精度求解方法介紹完畢。

在本公開的第四個示例性實施例中，提供了一種模擬復(fù)雜流體運動的高精度求解方法，先計算高精度的空間導(dǎo)數(shù)，進而計算出高精度的時間導(dǎo)數(shù)，再使用taylor格式求出下一步計算解。具體步驟如下：

步驟a.給定初始條件，旋轉(zhuǎn)合適的空間步長和時間步長，設(shè)定積分的終值時刻，最高時間積分精度為m階；

步驟b.在選定的m階時間精度下，采用n階高階精度致空間差分格式，通過遞歸微分進行空間差分計算，n與最高時間積分精度m的值相匹配；

步驟c.設(shè)置時間積分精度為k階，k＝1，2，…m，利用計算的高精度的空間導(dǎo)數(shù)，進而計算出高精度的時間導(dǎo)數(shù)；

步驟d.增加時間積分階數(shù)k，重復(fù)步驟c，直到k＝m循環(huán)結(jié)束，采用taylor級數(shù)法進行時間積分計算一步時間積分；

步驟e.重復(fù)上述步驟b～d，積分得到指定目標(biāo)時刻表示流體速度方向的數(shù)值解。

為了達到簡要說明的目的，上述第一實施例中任何可作相同應(yīng)用的技術(shù)特征敘述皆并于此，無需再重復(fù)相同敘述。

至此，本公開第四實施例的模擬復(fù)雜流體運動的高精度求解方法介紹完畢。

表1.采用不同算法實現(xiàn)的求解程序(線性情形)，計算時所需要的墻鐘時間，其中dp含義為doubleprecision，mp為multipleprecision，試驗平臺為使用intele5-26402.6ghzcpu的linux系統(tǒng)，時間單位為：秒。

除非有所知名為相反之意，本說明書及所附權(quán)利要求中的數(shù)值參數(shù)是近似值，能夠根據(jù)通過本公開的內(nèi)容所得的所需特性改變。具體而言，所有使用于說明書及權(quán)利要求中表示組成的含量、反應(yīng)條件等等的數(shù)字，應(yīng)理解為在所有情況中是受到「約」的用語所修飾。一般情況下，其表達的含義是指包含由特定數(shù)量在一些實施例中±10％的變化、在一些實施例中±5％的變化、在一些實施例中±1％的變化、在一些實施例中±0.5％的變化。

再者，單詞“包含”不排除存在未列在權(quán)利要求中的元件或步驟。位于元件之前的單詞“一”或“一個”不排除存在多個這樣的元件。

此外，除非特別描述或必須依序發(fā)生的步驟，上述步驟的順序并無限制于以上所列，且可根據(jù)所需設(shè)計而變化或重新安排。并且上述實施例可基于設(shè)計及可靠度的考慮，彼此混合搭配使用或與其他實施例混合搭配使用，即不同實施例中的技術(shù)特征可以自由組合形成更多的實施例。

在此提供的算法和顯示不與任何特定計算機、虛擬系統(tǒng)或者其它設(shè)備固有相關(guān)。各種通用系統(tǒng)也可以與基于在此的啟示一起使用。根據(jù)上面的描述，構(gòu)造這類系統(tǒng)所要求的結(jié)構(gòu)是顯而易見的。此外，本公開也不針對任何特定編程語言。應(yīng)當(dāng)明白，可以利用各種編程語言實現(xiàn)在此描述的本公開的內(nèi)容，并且上面對特定語言所做的描述是為了披露本公開的最佳實施方式。

本公開可以借助于包括有若干不同元件的硬件以及借助于適當(dāng)編程的計算機來實現(xiàn)。本公開的各個部件實施例可以以硬件實現(xiàn)，或者以在一個或者多個處理器上運行的軟件模塊實現(xiàn)，或者以它們的組合實現(xiàn)。本領(lǐng)域的技術(shù)人員應(yīng)當(dāng)理解，可以在實踐中使用微處理器或者數(shù)字信號處理器(dsp)來實現(xiàn)根據(jù)本公開實施例的相關(guān)設(shè)備中的一些或者全部部件的一些或者全部功能。本公開還可以實現(xiàn)為用于執(zhí)行這里所描述的方法的一部分或者全部的設(shè)備或者裝置程序(例如，計算機程序和計算機程序產(chǎn)品)。這樣的實現(xiàn)本公開的程序可以存儲在計算機可讀介質(zhì)上，或者可以具有一個或者多個信號的形式。這樣的信號可以從因特網(wǎng)網(wǎng)站上下載得到，或者在載體信號上提供，或者以任何其他形式提供。

本領(lǐng)域那些技術(shù)人員可以理解，可以對實施例中的設(shè)備中的模塊進行自適應(yīng)性地改變并且把它們設(shè)置在與該實施例不同的一個或多個設(shè)備中?？梢园褜嵤├械哪K或單元或組件組合成一個模塊或單元或組件，以及此外可以把它們分成多個子模塊或子單元或子組件。除了這樣的特征和/或過程或者單元中的至少一些是相互排斥之外，可以采用任何組合對本說明書(包括伴隨的權(quán)利要求、摘要和附圖)中公開的所有特征以及如此公開的任何方法或者設(shè)備的所有過程或單元進行組合。除非另外明確陳述，本說明書(包括伴隨的權(quán)利要求、摘要和附圖)中公開的每個特征可以由提供相同、等同或相似目的的替代特征來代替。并且，在列舉了若干裝置的單元權(quán)利要求中，這些裝置中的若干個可以是通過同一個硬件項來具體體現(xiàn)。

類似地，應(yīng)當(dāng)理解，為了精簡本公開并幫助理解各個公開方面中的一個或多個，在上面對本公開的示例性實施例的描述中，本公開的各個特征有時被一起分組到單個實施例、圖、或者對其的描述中。然而，并不應(yīng)將該公開的方法解釋成反映如下意圖：即所要求保護的本公開要求比在每個權(quán)利要求中所明確記載的特征更多的特征。更確切地說，如下面的權(quán)利要求書所反映的那樣，公開方面在于少于前面公開的單個實施例的所有特征。因此，遵循具體實施方式的權(quán)利要求書由此明確地并入該具體實施方式，其中每個權(quán)利要求本身都作為本公開的單獨實施例。

以上所述的具體實施例，對本公開的目的、技術(shù)方案和有益效果進行了進一步詳細(xì)說明，所應(yīng)理解的是，以上所述僅為本公開的具體實施例而已，并不用于限制本公開，凡在本公開的精神和原則之內(nèi)，所做的任何修改、等同替換、改進等，均應(yīng)包含在本公開的保護范圍之內(nèi)。