非靜壓模型垂向網(wǎng)格分離計算方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及一種非靜壓模型垂向網(wǎng)格分離計算方法，屬于技術(shù)領(lǐng)域。
【背景技術(shù)】
[0002] 目前，靜壓假設(shè)的數(shù)學(xué)模型已在河口海岸水動力研究中得到廣泛應(yīng)用。但在地形 劇烈變化、密度梯度較大、水流急劇變化等具有垂向流速變化較大的區(qū)域，靜壓模型難以準(zhǔn) 確模擬水動力的變化。故此，如何發(fā)展非靜壓模型，已成為近年河口海岸數(shù)值模擬研究的熱 點(diǎn)方向。大多數(shù)非靜壓模型壓力值的計算采用分步法，即單獨(dú)計算壓力的動壓部分與靜壓 部分，而求解動壓時需解泊松方程，此部分在非靜壓模型中耗時嚴(yán)重，影響了非靜壓模型的 計算效率及其應(yīng)用范圍。

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

[0003] 為解決現(xiàn)有技術(shù)的不足，本發(fā)明的目的在于提供一種非靜壓模型垂向網(wǎng)格分離計 算方法，通過分離求解泊松方程和動量方程的計算網(wǎng)格，解決非靜壓模型計算過程中耗時 過多的問題。
[0004] 為了實(shí)現(xiàn)上述目標(biāo)，本發(fā)明采用如下的技術(shù)方案：
[0005] -種非靜壓模型垂向網(wǎng)格分離計算方法，其特征是，在垂直方向采用不同網(wǎng)格求 解壓力泊松方程和動量方程，在求解泊松方程時采用粗網(wǎng)格，求解動量方程時采用細(xì)網(wǎng)格， 具體包括如下步驟：
[0006] 1)采用細(xì)網(wǎng)格求解不含動壓項(xiàng)的動量方程得到速度中間值；
[0007] 2)采用粗網(wǎng)格求解壓力泊松方程，通過分段平均得到粗網(wǎng)格速度中間值，得到粗 網(wǎng)格上的動壓值，并將此動壓值插值到細(xì)網(wǎng)格中得到細(xì)網(wǎng)格的動壓值；
[0008] 3)通過動量方程的動壓項(xiàng)修正速度場，得到最終的速度值。
[0009] 前述的非靜壓模型垂向網(wǎng)格分離計算方法，其特征是，所述步驟2)具體包括如下 步驟：
[0010]2. 1)通過分段平均法求得粗網(wǎng)格速度中間值if，分段平均法表達(dá)式如下：
[0011] (1)
[0012] 式中i表示粗網(wǎng)格上的網(wǎng)格編號，隊(duì)表示在粗網(wǎng)格第i個網(wǎng)格內(nèi)細(xì)網(wǎng)格的個數(shù)，Ji 表示第i個粗網(wǎng)格中細(xì)網(wǎng)格編號的最小值；
[0013]2. 2)通過泊松方程（3)求解粗網(wǎng)格上的動壓值p'，
[0014]
(2)
[0015] 式中D表示水深，P表示密度，x，y，0分別表示水平方向坐標(biāo)和垂向0坐 標(biāo)，X#，/表示笛卡爾坐標(biāo)系水平方向坐標(biāo)，u$，/，V表示速度中間值；
[0016] 2. 3)由粗網(wǎng)格的動壓值通過插值得到細(xì)網(wǎng)格上的動壓值p。
[0017] 前述的非靜壓模型垂向網(wǎng)格分離計算方法，其特征是，所述步驟2)中插值的求取 采用三次樣條插值法，插值函數(shù)需滿足如下條件：
[0018] r3)
[0019] (4)
[0020] (5)
[0021] 式中Z= -h表示位于水底處，h表示靜水面以下水深，n表示水體表面水位；p(1) 表示粗網(wǎng)格第i個網(wǎng)格中插值函數(shù)，p(1)'表示插值函數(shù)p(1)的一階導(dǎo)數(shù)，P(1)"表示插值函 數(shù)P(1)的二階導(dǎo)數(shù)，0i表示粗網(wǎng)格第i網(wǎng)格底面。坐標(biāo)，Pl表示粗網(wǎng)格第i網(wǎng)格動壓值。
[0022] 本發(fā)明所達(dá)到的有益效果：通過減少計算泊松方程的計算網(wǎng)格，在計算精度沒有 降低的條件下，顯著減少非靜壓模型的計算時間，測試算例顯示，在計算泊松方程的垂向網(wǎng) 格數(shù)最少可以僅為動量方程網(wǎng)格數(shù)的十分之一，采用網(wǎng)格分離方法后，計算時間最多可減 少約80%。
【附圖說明】
[0023] 圖1是非靜壓模型網(wǎng)格分離方法的網(wǎng)格布置圖；
[0024]圖2是非靜壓模型網(wǎng)格分離方法計算步驟圖；
[0025] 圖3是Lock-Exchange算例初始狀態(tài)圖；
[0026] 圖4是不同網(wǎng)格設(shè)置密度計算結(jié)果。
【具體實(shí)施方式】
[0027] 下面結(jié)合附圖對本發(fā)明作進(jìn)一步描述。以下實(shí)施例僅用于更加清楚地說明本發(fā)明 的技術(shù)方案，而不能以此來限制本發(fā)明的保護(hù)范圍。
[0028]本發(fā)明涉及的一種非靜壓模型垂向網(wǎng)格分離計算方法，在垂直方向采用不同網(wǎng)格 求解壓力泊松方程和動量方程，減少非靜壓模型求解壓力泊松方程的時間，從而提高計算 效率。
[0029] 求解泊松方程的網(wǎng)格數(shù)較少，稱為粗網(wǎng)格；求解動量方程的網(wǎng)格數(shù)較多，稱為細(xì)網(wǎng) 格，網(wǎng)格布置圖如圖1所示，實(shí)線表示粗網(wǎng)格，虛線表示細(xì)網(wǎng)格，P'，u'，w'分別表示粗網(wǎng)格 上動壓值和水平、垂向速度值，P，u,w表示細(xì)網(wǎng)格上動壓值和水平、垂向速度值。
[0030] 非靜壓模型目前主要的計算方法為分步法，即將壓力值分為動壓部分和靜壓部 分，第一步通過求解不含動壓項(xiàng)的動量方程得到速度中間值，然后求解壓力泊松方程，得到 動壓值，最后通過動量方程的動壓項(xiàng)修正速度場，得到最終的速度值。
[0031] 而本發(fā)明采用的非靜壓模型垂向網(wǎng)格分離計算方法是在非靜壓模型計算分步法 的基礎(chǔ)上進(jìn)行改進(jìn)，分步法的第一步不變，第二步由于求解泊松方程采用粗網(wǎng)格，需要求解 粗網(wǎng)格的速度中間值，因此增加通過分段平均得到粗網(wǎng)格速度中間值的步驟，在粗網(wǎng)格上 求解壓力泊松方程得到動壓值后需要增加動壓值插值到細(xì)網(wǎng)格的步驟，而后的第三步修正 速度場不變。
[0032] 第一步和第三步由于都不變，就不進(jìn)行闡述了，進(jìn)行改進(jìn)的第二步中包括如下步 驟：
[0033] 1)在細(xì)網(wǎng)格中求解不含動壓項(xiàng)的動量方程（1)，得到速度中間值U%
[0034]
[0035] 式中，分別為動量方程中底摩阻項(xiàng)、斜壓項(xiàng)和紊流擴(kuò)散項(xiàng)，上標(biāo)n表示模 型計算的第n時間步，F(xiàn),G,H表示通量項(xiàng)；
[0036]2)通過分段平均法求得粗網(wǎng)格速度中間值if，分段平均法表達(dá)式如下：
[0037] (2)
[0038] 式中i表示粗網(wǎng)格上的網(wǎng)格編號，隊(duì)表示在粗網(wǎng)格第i個網(wǎng)格內(nèi)細(xì)網(wǎng)格的個數(shù)，Ji 表示第i個粗網(wǎng)格中細(xì)網(wǎng)格編號的最小值；
[0039]3)通過泊松方程（3)求解粗網(wǎng)格上的動壓值p'，
[0040] (3)；
[0041] 4)由粗網(wǎng)格的動壓值通過插值得到細(xì)網(wǎng)格上的動壓值p，插值方法采用三次樣條 插值，插值函數(shù)需滿足如下條件：
[0042] C4)
[0043] (5)
[0044] 邊界條件：
（6)
[0045] 5)通過上一步得到的動壓值p，利用方程（7)修正細(xì)網(wǎng)格上的速度場，并利用方程 (8)更新自由表面水位：
[0046]
C7) (8)
[0048] 公式（7)中U(1)龍格庫塔方法第一階段的速度值，Sf:為動壓源項(xiàng)，At為時間步 長；公式（8)中D表示水深，u，v為水平流速；
[0049]6)此步驟根據(jù)具體算例確定，考慮紊流時，可求解紊流模型，包含k-G，大渦模擬 等；在包含鹽度、泥沙的模擬算例中，可通過鹽度、泥沙對流擴(kuò)散方程（9)求解鹽度、泥沙等 的分布和輸移。
[0050]
C9)
[0051] 其中可表示為溫度、鹽度、泥沙濃度等，i，j，k為水