張仁徽，陳宇里，耿釗

上海海事大學(xué)商船學(xué)院，上海 201306

網(wǎng)格尺寸對(duì)Delft3D有限元水流場(chǎng)仿真精度影響的分析

張仁徽，陳宇里，耿釗

上海海事大學(xué)商船學(xué)院，上海 201306

針對(duì)有限元水流場(chǎng)仿真中有效地選擇網(wǎng)格尺寸、平衡計(jì)算精度和計(jì)算成本的問(wèn)題，采用在網(wǎng)格形狀、邊界激勵(lì)類型、邊界激勵(lì)量大小等因素相同的條件下，輸入不同的網(wǎng)格尺寸分析了仿真的精度誤差的方法，建立了海南某港的Delft3D-FLOW模塊水流場(chǎng)仿真模型。研究結(jié)果表明：仿真計(jì)算的絕對(duì)誤差基本呈非線性關(guān)系，相同條件下靠近岸線的位置仿真誤差大于海域中間的位置仿真誤差。

有限元；水流場(chǎng)；網(wǎng)格尺寸；精度分析；Delft3D

隨著計(jì)算技術(shù)的發(fā)展，使用有限元方法對(duì)海洋流場(chǎng)仿真已經(jīng)成為一種實(shí)用的研究手段。對(duì)水流場(chǎng)的仿真經(jīng)歷了從二維到三維的發(fā)展階段，目前三維模式應(yīng)用已經(jīng)比較普遍，可以通過(guò)專業(yè)軟件的仿真更精確地反映出潮流的運(yùn)動(dòng)規(guī)律。水流場(chǎng)的精確仿真不僅是溢油模擬、海洋生態(tài)模擬等研究的基礎(chǔ)，更能為深水港、跨海大橋、航道治理等工程建設(shè)，貨運(yùn)船舶及軍艦的操作提供參考數(shù)據(jù)。

一般設(shè)想中，有限元仿真模擬使用的網(wǎng)格越精密其計(jì)算結(jié)果產(chǎn)生的誤差就越小，同時(shí)也帶來(lái)了更大規(guī)模的計(jì)算量，需要更多的計(jì)算時(shí)間。然而實(shí)際操作中，由于計(jì)算機(jī)數(shù)據(jù)字長(zhǎng)及對(duì)其舍入的作用，部分誤差反而會(huì)因?yàn)榫W(wǎng)格單元?jiǎng)澐诌^(guò)小而累積。因此在實(shí)際的工程運(yùn)用中，使用多大尺寸的網(wǎng)格能獲得計(jì)算精度、計(jì)算效率和技術(shù)成本的最好平衡，就成了技術(shù)人員面對(duì)的一個(gè)難題。針對(duì)這一問(wèn)題，本文對(duì)具體案例進(jìn)行水流場(chǎng)模擬，通過(guò)比對(duì)不同尺寸的網(wǎng)格得出的計(jì)算結(jié)果，對(duì)網(wǎng)格單元的劃分與計(jì)算精度的關(guān)系進(jìn)行探討。

1 基本原理

1.1 Delft3D模型

Delft3D水動(dòng)力學(xué)模型是由荷蘭水利研究所開發(fā)的一款水動(dòng)力學(xué)模型，可以在二維或三維層面通過(guò)有限元單元法對(duì)水流場(chǎng)進(jìn)行仿真。該仿真系統(tǒng)大致包含6個(gè)某塊：水流、水動(dòng)力、波浪、泥沙、水質(zhì)、生態(tài)等，在國(guó)內(nèi)外都有廣闊的運(yùn)用領(lǐng)域。我國(guó)從20世紀(jì)80年代中期就已經(jīng)開始將該系統(tǒng)運(yùn)用于相關(guān)項(xiàng)目，如長(zhǎng)江口、長(zhǎng)江重慶城區(qū)段、杭州灣、渤海灣、滇池、遼河、三江平原等，都相繼取得了不錯(cuò)的效果。

文中的模擬實(shí)驗(yàn)主要使用的是Delft3D中的FLOW模塊（水動(dòng)力模塊）。FLOW模塊的主要功能是模擬眾多不同的水流條件，如湖泊和大洋中的風(fēng)生流、河口及入?？诘某毕骱蛯?shí)驗(yàn)水槽中的紊流等等。FLOW模塊能夠建立不同規(guī)模的直線或曲線網(wǎng)格來(lái)計(jì)算非穩(wěn)定流，并且其在計(jì)算中提供了豐富的開邊界條件和初始條件，如不同來(lái)源的水流流速、流量、糙率或是天文潮等，可以模擬風(fēng)對(duì)水體表面的作用，甚至能在每一個(gè)網(wǎng)格點(diǎn)上設(shè)置獨(dú)立的水深數(shù)據(jù)?？梢哉f(shuō)，F(xiàn)LOW模塊是Delft3D系統(tǒng)的一個(gè)核心模塊。

1.2 FLOW模塊的計(jì)算原理

FLOW模塊建立在Navier-Stokes方程的基礎(chǔ)上，采用有限差分法中交替隱式法（ADI）對(duì)相應(yīng)坐標(biāo)系下的控制方程組進(jìn)行離散求解。忽略在垂直方向動(dòng)量方程中垂直加速度的影響，進(jìn)而推導(dǎo)出靜水壓強(qiáng)假定下的水流方程。

本文的仿真計(jì)算中，主要涉及兩組方程（連續(xù)方程與動(dòng)量守恒方程）與邊界條件表達(dá)式。

1.2.1 模型控制方程

波動(dòng)方程（Cartesian坐標(biāo)系）：

式中：η為平均海平面以上的水位，m；H為總水深，H＝h＋η，m；u、v為深度平均速度的東分量和北分量，m／s；x、y為空間坐標(biāo)，m；t為時(shí)間坐標(biāo)，s；h為凈水深，m；f為Coriolis參數(shù)，1／s；g為重力加速度，m／s2；為風(fēng)應(yīng)力分量；為底應(yīng)力分量。

式中：C為Chezy系數(shù)，C＝nH／6，其中n為Manning系數(shù)；ρ為海水密度，kg／m3。

1.2.2 界條件和初始條件

潮波計(jì)算采用下列邊界條件。沿閉邊界（水陸邊界）垂直海岸的流通量等于零：

Vn＝0

式中：V＝（u，v），n是指向邊界外的單位法向量。

關(guān)于對(duì)流項(xiàng)，在開邊界當(dāng)海水向計(jì)算區(qū)域流進(jìn)時(shí)，法向流速的導(dǎo)數(shù)等于零。

2 仿真過(guò)程及分析

2.1 案例介紹

本文所模擬的水流場(chǎng)為海南省西南部某港，東經(jīng)108°56′30″～109°48′28″，北緯18°09′34″～18°37′27″。該水域?qū)俨徽?guī)日潮混合潮型，以日潮為主，且有明顯的日潮不等現(xiàn)象。根據(jù)該港波浪觀測(cè)報(bào)告，本區(qū)潮流具有往復(fù)流特征，流向大都集中在290°～350°和120°～180°之間。大潮流速明顯大于小潮流速。大潮最大流速為0.76 m／s，小潮最大流速為0.42 m／s；大潮最大垂線平均流速為0.67 m／s，小潮最大垂線平均流速為0.32 m／s。受風(fēng)、浪及地形等影響，本區(qū)亦有余流存在，但流速較小，一般在8～10 cm／s。

2.2 物理參數(shù)

由于地理結(jié)構(gòu)的變化，不同算例中的物理參數(shù)和邊界條件都不盡相同。為有效控制仿真計(jì)算的穩(wěn)定性，本文通過(guò)查閱大量文獻(xiàn)，在仿真中將相關(guān)物理參數(shù)設(shè)定如表1所示。

表1 物理參數(shù)值的設(shè)定

本文選用的仿真模型在數(shù)值計(jì)算上使用流體表征粘性影響系數(shù)為chezy系數(shù)，在本案例中chezy值恒為65。水平渦流黏度為1㎡／s。

2.3 網(wǎng)格形狀的選擇

有限元仿真常用的網(wǎng)格有二維三角形，二維四邊形和三維四面體元、五面體元和六面體元。他們的邊界形狀主要有直線型、曲線型和曲面型。單元最佳形狀是正多邊形或正多面體，其具有良好相容性、逼近精確性和剖分過(guò)渡性和自適應(yīng)性，單元之間過(guò)渡相對(duì)平穩(wěn)。以二維三角形和二維四邊形為例，相關(guān)文獻(xiàn)表明，在步長(zhǎng)基本相同的情況下，四邊形網(wǎng)格的計(jì)算時(shí)間、解的精度等方面均優(yōu)于三角形網(wǎng)格，尤其是收斂速度方面，明顯優(yōu)于三角形網(wǎng)格。

本文所選的有限元仿真模型所使用的網(wǎng)格即為二維四邊形網(wǎng)格，能夠較好地滿足計(jì)算穩(wěn)定性和精度的要求。

2.4 網(wǎng)格規(guī)模的確定

在仿真中，不合適的網(wǎng)格規(guī)模將影響網(wǎng)格的正交性，甚至引起邊界網(wǎng)格發(fā)生畸變，直接影響仿真的計(jì)算精度和計(jì)算用時(shí)。對(duì)此，本文所選用的仿真模型對(duì)網(wǎng)格正交性的最低要求為：cos?＜0.02，?為網(wǎng)格橫縱方向的夾角。

本文案例中仿真的水域面積約為50 000平方公里，通過(guò)實(shí)驗(yàn)表明，當(dāng)網(wǎng)格尺寸小于2 000，即網(wǎng)格長(zhǎng)度小于0.6 km時(shí)，網(wǎng)格的正交性才能滿足cos?＜0.02。因此本文選用的網(wǎng)格尺寸，最大為2 000。

2.5 仿真過(guò)程

本文在對(duì)選定水域進(jìn)行仿真使用的開邊界初始條件中的主要激勵(lì)，為對(duì)該水域產(chǎn)生影響最大的4個(gè)主要天文分潮：O1、K1、S2、M2。4個(gè)天文分潮的振幅、遲角的數(shù)據(jù)由潮汐表上潮汐調(diào)和常數(shù)求得。具體數(shù)值為M2分潮振幅0.2，遲角345°；S2分潮振幅0.07，遲角70°；O1分潮振幅0.3，遲角295°；K1分潮振幅0.3，遲角345°。

為便于比較不同網(wǎng)格計(jì)算精度，本文以精度輸出和預(yù)估計(jì)算量為考量標(biāo)準(zhǔn)，采用了若干種不同疏密程度網(wǎng)格模型。全部尺寸的網(wǎng)格采用的都是均勻網(wǎng)格。為研究精度較大時(shí)網(wǎng)格尺寸改變帶來(lái)的誤差變化坡度，在500以前網(wǎng)格尺寸以100為單位減小，500以后網(wǎng)格尺寸以倍數(shù)減小。網(wǎng)格尺寸的具體參數(shù)見(jiàn)表2。

表2 網(wǎng)格尺寸參數(shù)

為確保網(wǎng)格邊界的計(jì)算畸變不影響計(jì)算結(jié)果，在建立網(wǎng)格時(shí)盡可能包含附近較長(zhǎng)的海岸線，并且在保證靠近邊界的網(wǎng)格正交性的前提下使其盡可能與岸線貼合。

在模擬中，本文共選取了3個(gè)測(cè)點(diǎn)：測(cè)點(diǎn)1靠近海岸，測(cè)點(diǎn)2在海域中間，測(cè)點(diǎn)3在遠(yuǎn)離輸入激勵(lì)的邊界。如圖1所示。

圖1 仿真中3個(gè)觀測(cè)點(diǎn)的位置

建立網(wǎng)格模型和輸入邊界激勵(lì)條件后，在工程機(jī)上運(yùn)行Delft3D的FLOW模塊進(jìn)行有限元仿真計(jì)算，設(shè)定的仿真時(shí)間段為30 d。在經(jīng)過(guò)一段時(shí)間的激勵(lì)輸入后，可得較為穩(wěn)定的仿真狀態(tài)。同時(shí)可以輸出海域中每個(gè)網(wǎng)絡(luò)節(jié)點(diǎn)的即時(shí)潮流全息數(shù)據(jù)（流速、流向等）。圖2為觀測(cè)點(diǎn)1在狀態(tài)穩(wěn)定后某一天24 h之內(nèi)的潮流信息。

圖2 觀測(cè)點(diǎn)1在1 d內(nèi)的潮流信息

在輸出的觀測(cè)點(diǎn)數(shù)據(jù)中選取13 d，每天24 h的模擬穩(wěn)定潮流信息。可得到3個(gè)觀測(cè)點(diǎn)在激勵(lì)條件和仿真時(shí)間相同的情況下，使用不同尺寸的網(wǎng)格仿真計(jì)算出的節(jié)點(diǎn)流速和流向。

本文以觀測(cè)站的實(shí)際觀測(cè)資料作為標(biāo)準(zhǔn)數(shù)據(jù)，將不同網(wǎng)格尺寸的數(shù)據(jù)與實(shí)際測(cè)量值對(duì)比，分別計(jì)算它們的相對(duì)誤差，并對(duì)所有誤差進(jìn)行橫向分析。

圖3 觀測(cè)點(diǎn)2在100×100網(wǎng)格條件下的潮流信息對(duì)比

驗(yàn)證資料取用2010年6月10日～6月11日（大潮）的水文測(cè)量資料，對(duì)大潮的潮位過(guò)程進(jìn)行潮位驗(yàn)證。流速驗(yàn)證取用3個(gè)測(cè)點(diǎn)的垂線平均流速值。分別將3個(gè)觀測(cè)點(diǎn)各種網(wǎng)格精度下的仿真計(jì)算結(jié)果與真實(shí)觀測(cè)值進(jìn)行比對(duì)。圖3為觀測(cè)點(diǎn)2在網(wǎng)格精度為100×100條件下的矢量對(duì)比圖。

依次將3個(gè)觀測(cè)點(diǎn)的仿真流向和流速與觀測(cè)值進(jìn)行對(duì)比，得到仿真計(jì)算精度誤差變化。圖4由上至下依次為不同尺寸的網(wǎng)格模型在3個(gè)觀測(cè)點(diǎn)的流向和流速計(jì)算精度誤差變化。

圖4 觀測(cè)點(diǎn)計(jì)算精度誤差變化趨勢(shì)

2.6 仿真結(jié)果分析

由圖4可以看出，在FLOW模塊的仿真中，總體趨勢(shì)上尺寸越精細(xì)，仿真計(jì)算值與觀測(cè)值接近度越高，網(wǎng)格尺寸越大，絕對(duì)誤差也呈增大趨勢(shì)?？梢?jiàn)在仿真中由計(jì)算機(jī)數(shù)據(jù)步長(zhǎng)和取舍造成的累積誤差并不明顯。

在網(wǎng)格尺寸由100到500增大的過(guò)程中，誤差曲線較為陡峭，誤差的上升速度相對(duì)較快。而在網(wǎng)格尺寸到達(dá)500以后，雖然尺寸增大的速度增加了，誤差曲線卻趨于平緩，誤差的增速趨于穩(wěn)定。

對(duì)比3個(gè)不同的觀測(cè)點(diǎn)，在網(wǎng)格尺寸相同的條件下，由于所處位置的不同，輸入的仿真數(shù)據(jù)誤差量也有著較大區(qū)別。以水流速度誤差為例，可以看到靠近岸線的觀測(cè)點(diǎn)1的最大誤差超過(guò)了300%，遠(yuǎn)離初始邊界的觀測(cè)點(diǎn)3最大誤差也達(dá)到120%，而處于仿真海域中間位置的觀測(cè)點(diǎn)2即使在網(wǎng)格尺寸相當(dāng)大的情況下，仿真計(jì)算誤差仍然沒(méi)有超過(guò)20%，并且增長(zhǎng)曲線非常平緩。這是由于靠近岸線的觀測(cè)點(diǎn)所處的網(wǎng)格需要貼合岸線的曲線，當(dāng)網(wǎng)格尺寸變大的時(shí)候，會(huì)導(dǎo)致邊界網(wǎng)格的正交性變差、不平滑，甚至出現(xiàn)畸變的情況，而遠(yuǎn)離激勵(lì)輸入的初始邊界的觀測(cè)點(diǎn)所處位置受到激勵(lì)驅(qū)動(dòng)所達(dá)到的狀態(tài)不夠穩(wěn)定，這都是導(dǎo)致仿真誤差增大的原因。

圖5 同一觀測(cè)點(diǎn)不同網(wǎng)格尺寸下的流速仿真值

由圖5可以看到，對(duì)于同一個(gè)觀測(cè)點(diǎn)，隨著網(wǎng)格尺寸的增大，仿真模擬的穩(wěn)定性也大大降低。在實(shí)際仿真輸出中，常取某個(gè)時(shí)刻的數(shù)值作為參考，因?yàn)榉€(wěn)定性的下降直接導(dǎo)致仿真誤差的加劇。而在網(wǎng)格尺寸等仿真條件相同時(shí)，某點(diǎn)仿真流速和流向的輸出誤差值，基本處于同一數(shù)量級(jí)，沒(méi)有明顯差異。

3 結(jié)束語(yǔ)

通過(guò)對(duì)仿真結(jié)果的精度進(jìn)行比較，工程技術(shù)人員在使用FLOW進(jìn)行水流場(chǎng)模擬的時(shí)候，若工程對(duì)象為碼頭、溢油模擬等靠近岸線的內(nèi)容，則應(yīng)選用較為精細(xì)的網(wǎng)格；若工程對(duì)象為在海域當(dāng)中，如為模擬船舶操縱而進(jìn)行水流場(chǎng)仿真等，可以適當(dāng)考慮使用較粗糙的網(wǎng)格，以在可接受的誤差范圍內(nèi)節(jié)省仿真時(shí)間和仿真成本。

研究結(jié)果表明，相同條件下靠近岸線的位置仿真誤差誤差曲線較為陡峭，誤差的上升速度相對(duì)較快；海域中間的位置仿真誤差增長(zhǎng)曲線非常平緩。

本次研究在相同的水文條件下，重點(diǎn)關(guān)注不同網(wǎng)格尺寸對(duì)于精度誤差的影響，在今后的研究中可以更多分析不同形狀的網(wǎng)格對(duì)于仿真精度的影響，以及不同形狀的網(wǎng)格適用于何種工程仿真之中。

［1］徐騰飛，趙人達(dá).計(jì)算機(jī)精度對(duì)隨機(jī)場(chǎng)網(wǎng)格尺寸的影響［J］.計(jì)算機(jī)工程與應(yīng)用，2011，47（17）：210-211.

［2］陳榮昌，趙前，鄧健，等.Delft3D和Oilmap在內(nèi)河溢油模擬中的聯(lián)合應(yīng)用研究［J］.中國(guó)水運(yùn)，2011，11（4）：65-67.

［3］楊益洪.西湖混合流三維數(shù)值模擬研究［D］.杭州：浙江大學(xué)，2007：24-36.

［4］修榮榮，徐海明，黃善波.網(wǎng)格單元形狀對(duì)數(shù)值計(jì)算的影響［J］.工程熱物理學(xué)報(bào)，2004，25（2）：317-319.

［5］郭軍朝，谷正氣，孟慶超.基于湍流模型在轎車底部流場(chǎng)的仿真分析［J］.計(jì)算機(jī)仿真，2007，24（4）：258-261.

［6］龔曙光，謝桂蘭，邱愛(ài)紅，等.CAE仿真分析中計(jì)算精度與網(wǎng)格劃分關(guān)系的探討［J］.設(shè)計(jì)與研究，2008，21（12）：35-38.

［7］秦權(quán)，林道錦，梅剛.結(jié)構(gòu)可靠度隨機(jī)有限理論及工程應(yīng)用［M］.北京：清華大學(xué)出版社，2006：78-92.

［8］史恒通，王成華.土坡有限元穩(wěn)定分析若干問(wèn)題探討［J］.巖土力學(xué)，2000，21（2）：152-155.

［9］裴利劍，蘇蓮萍，張瑋，等.網(wǎng)格疏密和模型邊界對(duì)強(qiáng)度折減法計(jì)算精度的影響［J］.昆明冶金高等專科學(xué)校學(xué)報(bào)，2010，26（5）：64-68.

［10］裴利劍，姚永仲.邊坡穩(wěn)定有限元中泊松比的影響［J］.昆明冶金高等?？茖W(xué)校學(xué)報(bào)，2009，25（3）：56-59.

［11］ZHANG J，ELLINGWOOD B.Orthogonal series expansions of random fields in reliability analysis［J］.Journal of Engi-neering Mechanics，1993，120（12）：2660-2677.

［12］MCCAY D F.Development and application of damage as-sessmentmodeling：examp le assessment for the North Cape oil spill［J］.Marine Pollution Bulletin，2003，47（9／10／11／12）：341-359.

Influence analysis of themesh size on the precision of flow field simulation of the the Delft3D finite element

ZHANG Renhui，CHEN Yuli，GENG Zhao
Merchant Marine College，Shanghai Maritime University，Shanghai201306，China

：In order to select themesh size efficiently，and balance the calculation precision and computational ex-pense in the simulation of finite element flow field，a Delft3D-FLOW field simulation model at a portof Hainan was established.The precision error of the simulation was analyzed by inputting variousmesh sizes with constantmesh shapes，types and values of boundary incentive.The analysis result shows that，the absolute simulation errors are non-linear，and the simulation error near the shoreline is greater than that in themiddle of the ocean area under the same simulation conditions.

finite element；flow field；mesh size；precision analysis；Delft3D

TB115

：A

：1009-671X（2015）01-057-05

10.3969／j.issn.1009-671X.201402003

http：//www.cnki.net／kcms／detail／23.1191.U.20150112.1530.005.htm l

2014-02-24.

日期：2015-01-12.

交通運(yùn)輸部應(yīng)用基礎(chǔ)研究基金資助項(xiàng)目（2013329810300）；上海市教委科研創(chuàng)新基金資助重點(diǎn)項(xiàng)目（13ZZ124）.

張仁徽（1989-），男，碩士研究生；陳宇里（1979-）男，副教授，博士.

張仁徽，E-mail：306852381＠qq.com.