海冰動力過程的改進(jìn)離散元模型及在渤海的應(yīng)用

季順迎，王安良，米麗麗，劉煜，李寶輝

(1. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連 116023；2. 國家海洋環(huán)境預(yù)報中心 國家海洋局海洋災(zāi)害預(yù)報技術(shù)研究重點實驗室，北京 100081)

海冰動力過程的改進(jìn)離散元模型及在渤海的應(yīng)用

季順迎1，王安良2，米麗麗1，劉煜2，李寶輝2

(1. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連 116023；2. 國家海洋環(huán)境預(yù)報中心 國家海洋局海洋災(zāi)害預(yù)報技術(shù)研究重點實驗室，北京 100081)

海冰的斷裂、重疊和堆積等離散分布特性廣泛地存在于極區(qū)和副極區(qū)的不同海域，并對海冰的生消、運(yùn)移過程有著重要影響。針對海冰在不同尺度下的離散分布特點，發(fā)展海冰動力過程的離散元方法有助于完善海冰數(shù)值模式，提高海冰數(shù)值模擬的計算精度。為此，本文針對海冰生消運(yùn)移過程中的非連續(xù)分布和形變特性，發(fā)展了適用于海冰動力過程的改進(jìn)離散元模型(MDEM)。不同于傳統(tǒng)離散元方法，該模型將海冰離散為具有一定厚度、尺寸和密集度的圓盤單元。海冰單元設(shè)為諸多浮冰塊的集合體，其在運(yùn)移和相互接觸碰撞過程中，依照質(zhì)量守恒發(fā)生單元尺寸、密集度和厚度的相應(yīng)變化?；诤１x散性和流變性的特點，該模型采用黏彈性接觸本構(gòu)模型計算單元間的作用力，并依據(jù)Mohr-Coulomb準(zhǔn)則計算海冰法向作用下的塑性變形及切向摩擦力。為驗證該模型的可靠性，本文對海冰在規(guī)則水域內(nèi)的運(yùn)移和堆積過程進(jìn)行了分析，離散元計算結(jié)果與解析值相一致；此外，對旋轉(zhuǎn)風(fēng)場下海冰漂移規(guī)律的模擬進(jìn)一步驗證了本文方法的精確性。在此基礎(chǔ)上，對渤海遼東灣的海冰動力過程進(jìn)行了48 h數(shù)值分析，計算結(jié)果與衛(wèi)星遙感資料和油氣作業(yè)區(qū)的海冰現(xiàn)場監(jiān)測數(shù)據(jù)吻合良好。在下一步工作中將考慮海冰離散元模擬中的熱力因素影響，發(fā)展具有凍結(jié)、斷裂效應(yīng)的海冰離散元模型，更精確地模擬海冰動力-熱力耦合作用下的生消和運(yùn)移過程。

海冰動力學(xué)；離散元模型；堆積冰；渤海

1 引言

近年來，隨著極區(qū)海冰覆蓋面積和冰量的急劇變化，海洋資源開采向寒區(qū)的不斷延伸，北極夏季航道的初步開通，海冰這一自然現(xiàn)象及其引發(fā)的工程災(zāi)害引起了人們更加密切的關(guān)注。海冰在極區(qū)與大氣、海洋的相互作用是影響全球氣候演變趨勢的重要因素，在渤海、Baltic海和Okhotsk海等副極區(qū)海域又是影響油氣開發(fā)、冰區(qū)航運(yùn)和水產(chǎn)養(yǎng)殖的重要環(huán)境條件[1—3]。海冰在不同尺度下均表現(xiàn)出很強(qiáng)的離散分布特性，如大中尺度下的非連續(xù)流變行為，以及小尺度下的斷裂、重疊和堆積現(xiàn)象[4—8]。為對不同尺度下海冰的動力-熱力過程進(jìn)行精確、高效的數(shù)值模擬，需要針對海冰的離散分布規(guī)律及其動力學(xué)特性建立相適應(yīng)的數(shù)值方法和本構(gòu)模型。

自20世紀(jì)70年代以來，對極區(qū)及副極區(qū)的海冰演化過程進(jìn)行數(shù)值模擬，人們分別建立了歐拉坐標(biāo)下的有限差分法(FDM)[9—12]、拉格朗日坐標(biāo)下的光滑質(zhì)點流體動力學(xué)方法(SPH)[13—14]，以及兩種坐標(biāo)相結(jié)合的質(zhì)點網(wǎng)格法(PIC)[15—18]。在以上數(shù)值方法中，F(xiàn)DM方法是發(fā)展最早，也是目前應(yīng)用范圍最廣的海冰數(shù)值方法，并已應(yīng)用于北極、Baltic海和渤海等海域[9—10，12，19]。最近，有限體積法、有限元方法、質(zhì)量點法和混合拉格朗日-歐拉法等也被應(yīng)用于海冰的動力學(xué)數(shù)值計算[20—23]。在海冰動力學(xué)的數(shù)值方法中，一般采用黏塑性[9，12，24—25]、彈塑性[26]或黏彈塑性[14，27—28]本構(gòu)模型進(jìn)行海冰內(nèi)力的計算。特別是，光滑質(zhì)點流體動力學(xué)方法可有效克服因網(wǎng)格差分而產(chǎn)生的平流項擴(kuò)散的現(xiàn)象[29—30]，可提高對局地海域海冰動力現(xiàn)象模擬以及邊緣線預(yù)測的精度，而被引入到海冰動力學(xué)過程的數(shù)值模擬中。

上述數(shù)值方法均將海冰視為連續(xù)介質(zhì)，能夠理想地模擬海冰在宏觀尺度下的流變特性和平均分布狀態(tài)，具有計算效率高的特點。然而，大量的現(xiàn)場監(jiān)測和衛(wèi)星遙感資料表明，海冰在極區(qū)及其冰緣區(qū)、副極區(qū)等不同海域均表現(xiàn)出很強(qiáng)的離散分布特性，且平整冰、堆積冰、冰脊和水道相互交織出現(xiàn)[31—33]。冰塊尺寸在不同尺度下有很大的變化范圍，可從極區(qū)大尺度下的100 km以上到小尺度下的1 km以下[3，34]。在波浪作用下，海冰還可破碎為10 m甚至更小的冰塊。海冰的這種離散分布特性早在20世紀(jì)80年代初就得到了人們的重視，并明確指出海冰具有類似顆粒材料的力學(xué)行為[4，35]。Shen等最早將顆粒流理論應(yīng)用到海冰動力學(xué)中，建立了海冰碰撞流變學(xué)，并用于碎冰區(qū)的海冰動力學(xué)數(shù)值模擬[36]。針對海冰的離散分布特性，基于顆粒介質(zhì)理論的海冰動力學(xué)本構(gòu)模型也相繼發(fā)展起來[37-38]。但這些本構(gòu)模型并不能考慮細(xì)觀尺度下海冰的分布形態(tài)、碰撞接觸作用，其計算參數(shù)也沒有考慮海冰類型、冰塊尺寸等細(xì)觀因素的影響。最近，Hopkins的一系列研究表明，海冰離散元模型不但可用于小尺度下冰塊間的碰撞作用，中尺度下冰脊、冰隙的演化規(guī)律，甚至還可對極區(qū)大尺度下的海冰動力過程進(jìn)行數(shù)值模擬[39—40]。在以上海冰的離散元模型中，依據(jù)海冰的分布狀態(tài)將其離散成多個具有厚度、尺寸和速度的計算單元。當(dāng)考慮冰塊凍結(jié)時，該方法還能夠合理地模擬海冰動力學(xué)過程中的斷裂、重疊和堆積現(xiàn)象，具有物理意義明確、計算精度高的優(yōu)點。然而，由于海冰在不同尺度下分布形態(tài)的差異以及相互作用的復(fù)雜性，對非規(guī)則形狀海冰單元的構(gòu)造，以及單元間的接觸碰撞模型還需要進(jìn)一步研究。

近年來，離散元模型的發(fā)展和完善對其在海冰動力學(xué)數(shù)值模擬中的應(yīng)用有很大的促進(jìn)作用。但也注意到，離散元模型在海冰動力方程的計算中采用顯式計算，時間步長小，其計算效率明顯低于其他計算模型，目前還不能滿足海冰數(shù)值預(yù)報和長期海冰動力學(xué)模擬的時效要求。此外，采用傳統(tǒng)的離散元方法需要對每個冰塊建立相應(yīng)的計算模型。由于海冰單元數(shù)目巨大，這在當(dāng)前計算條件下難以實現(xiàn)；對每個冰塊的初始位置、尺寸、厚度等參數(shù)進(jìn)行精確確定也是目前海冰監(jiān)測技術(shù)不能達(dá)到的。為此，基于海冰動力學(xué)SPH數(shù)值模擬的思路，將每個海冰單元設(shè)為諸多浮冰塊的集合體，其剛度設(shè)為密集度的函數(shù)并考慮海冰單元的塑性變形。通過對傳統(tǒng)海冰離散元模型進(jìn)行改進(jìn)可對海冰動力學(xué)過程進(jìn)行相對精確、高效的數(shù)值模擬。

本文通過發(fā)展海冰動力過程的改進(jìn)離散元模型，對不同邊界條件下海冰運(yùn)移、堆積過程進(jìn)行數(shù)值分析以檢驗該模型的可靠性。通過對遼東灣海冰動力演化過程的離散元模擬，對油氣作業(yè)區(qū)的海冰參數(shù)進(jìn)行提取，并與現(xiàn)場海冰監(jiān)測資料進(jìn)行了對比分析，進(jìn)一步驗證改進(jìn)離散元模型的有效性和精確性。

2 海冰動力過程的改進(jìn)離散元模型

在海冰動力學(xué)的改進(jìn)離散元模型中，海冰單元設(shè)為諸多冰塊的集合體，相互接觸、碰撞作用受其平均厚度、密集度、單元尺寸等參數(shù)影響。當(dāng)單元間的法向應(yīng)力超過其壓縮強(qiáng)度時，海冰發(fā)生塑性變形，并由此導(dǎo)致單元尺寸的變化。海冰單元的運(yùn)動規(guī)律受控于風(fēng)和流的拖曳力、海面傾斜引起的海冰壓強(qiáng)梯度力、冰間作用力等因素。

2.1 海冰單元的運(yùn)動方程

海冰的運(yùn)動主要受控于單元間的作用力、風(fēng)和流的拖曳力、科氏力，以及海面傾斜引起的壓強(qiáng)梯度力，其動力方程由牛頓定律來描述，即：

(1)

式中，M為單位面積海冰質(zhì)量，即M=NρiHi，V為海冰速度矢量，f為科氏參數(shù)，且f=2ωesinφ，ωe為地球轉(zhuǎn)角速度，φ為地理緯度，K為垂直于海面的單位矢量，ξw為瞬時海面梯度，海冰在單位面積上風(fēng)和流的拖曳力為Vai和Vwi，其中Ca和Cw為風(fēng)和流的拖曳系數(shù)，Vai和Vwi為相對于海冰的風(fēng)速和流速矢量，·σ為海冰單元間相互作用引起的單位面積上的內(nèi)力矢量。

2.2 海冰單元間的接觸模型

海冰單元之間的作用力主要包括單元間在法向和切向上因相互重疊而引起的彈性力和因相對速度引起的黏性力。此外，切線方向上還需考慮基于Mohr-Coulomb準(zhǔn)則的滑動摩擦。海冰單元MA和MB間的接觸模型如圖1所示，其中Kn和Kt分別為海冰單元間的法向和切向剛度系數(shù)，Cn和Ct分別為海冰單元間的法向和切向黏性系數(shù)。

圖1 海冰單元間的接觸模型Fig.1 Contact model of sea ice elements

在法線方向，海冰單元間的接觸力為：

(2)

在切線方向，基于Mohr-Coulomb摩擦定律，海冰單元的作用力為：

(3)

海冰單元的法向和切向剛度與其密集度密切相關(guān)。這里依據(jù)河冰和海冰力學(xué)性質(zhì)與其密集度的對應(yīng)關(guān)系，可取[14，41]：

(4)

式中，E為海冰彈性模量，Hi為海冰厚度，Ni為海冰單元的密集度，Nmax為其最大值，這里取Nmax=100%，j值通常取為15。此外，海冰力學(xué)參數(shù)還可設(shè)為密集度的負(fù)指數(shù)函數(shù)[9]。在離散元模型中，一般取單元的切向剛度Kt=0.5Kn。

海冰單元的黏性系數(shù)Cn可由其質(zhì)量、法向剛度和無量綱黏性系數(shù)確定，即[42]：

(5)

當(dāng)采用線黏彈性模型計算海冰單元的作用力時，兩個自由單元間的接觸時間被定義二元接觸時間，即[42]：

(6)

式中，Tbc為海冰單元的二元接觸時間，即兩個海冰單元從碰撞到分離的接觸時間。在線黏彈性模型中，它是一個與單元大小和材料性質(zhì)相關(guān)的常數(shù)，其可用于確定海冰動力學(xué)計算的時間步長Δt。一般取時間步長Δt為二元接觸時間Tbc的1/50。

2.3 海冰單元的塑性變形

當(dāng)海冰單元在外力作用下相互擠壓時，單元間的法向應(yīng)力可由法向力Fn和接觸面積S確定，即σn=Fn/S。在海冰動力學(xué)的改進(jìn)離散元模型中，海冰單元為諸多碎冰塊的集合體。其法向壓縮強(qiáng)度可根據(jù)單元內(nèi)浮冰塊在豎直方向上的靜水壓力，并由Mohr-Coulomb摩擦定律進(jìn)行確定[14，39]。該方法已成功地應(yīng)用于海冰動力學(xué)的SPH數(shù)值模擬中[13—14]。當(dāng)海冰單元的法應(yīng)力超過其壓縮強(qiáng)度時，將會發(fā)生塑性形變。相應(yīng)地，海冰單元將縮小，并引起單元密集度和厚度的改變。

依據(jù)顆粒材料的塑性失效準(zhǔn)則并考慮冰內(nèi)靜水壓力效應(yīng)，海冰單元的壓縮強(qiáng)度σc寫為[41]：

(7)

(8)

式中，ρi,ρw分別為海冰和海水的密度，g為重力加速度。

當(dāng)海冰單元發(fā)生塑性形變后，其面積會相應(yīng)減小，并由此導(dǎo)致密集度的增加，海冰發(fā)生輻合現(xiàn)象。當(dāng)海冰密集度達(dá)到最大值Nmax時，若海冰內(nèi)力進(jìn)一步增大，此時單元內(nèi)的海冰將發(fā)生堆積現(xiàn)象。相反，當(dāng)海冰單元間的應(yīng)力小于其壓縮強(qiáng)度時，海冰單元在靜水壓力作用下發(fā)生輻散現(xiàn)象，并由此導(dǎo)致海冰單元面積的增加和密集度的降低。該海冰在內(nèi)力作用下的輻合和輻散過程示意圖如圖2所示。在該輻合和輻散過程中，海冰單元的質(zhì)量保持恒定，而其面積、密集度和厚度發(fā)生相應(yīng)變化。

圖2 海冰輻合-輻散過程及海冰單元相應(yīng)面積、密集度的變化Fig.2 Divergent and convergent processes and the corresponding area，concentration of sea ice elements

3 規(guī)則區(qū)域內(nèi)海冰動力過程的數(shù)值模擬

為驗證MDEM模型的可靠性和準(zhǔn)確性，下面分別對規(guī)則區(qū)域內(nèi)兩個典型的海冰動力過程進(jìn)行數(shù)值分析。一是變寬度水道中海冰在風(fēng)和流拖曳下的漂移和堆積過程，另一個是矩形區(qū)域內(nèi)海冰在旋轉(zhuǎn)風(fēng)場作用下的動力過程。

3.1 變寬度水道內(nèi)海冰的漂移和堆積過程

在海冰動力學(xué)研究中，水道內(nèi)海冰漂移和堆積的動力過程分析是檢驗海冰數(shù)值方法的一個有效途徑，其平均冰厚的分布特征得到了現(xiàn)場觀測數(shù)據(jù)、室內(nèi)模型實驗和理論模型解析解的驗證[23，41，43—44]。在此基礎(chǔ)上，本文采用變寬度水道形式，并在水道內(nèi)放置一個障礙物。通過分析海冰的繞流特征進(jìn)一步檢驗改進(jìn)離散元模型對海冰動力過程模擬的適用性。

本文采用的變寬度水道如圖3所示，總長度L=170 km，左側(cè)和右側(cè)寬度分別為W1=30 km和W1=14 km。海冰初始厚度Hi0=1.0 m，初始密集度Ni0=100%，冰區(qū)初始長度L1=30 km，其他主要計算參數(shù)列于表1中。數(shù)值模擬的第1 d、2 d、3 d和7 d的海冰單元速度矢量和平均厚度分布如圖4所示。在風(fēng)和流的作用下，海冰繞過三角形障礙物向右漂移，并在右側(cè)邊界發(fā)生堆積。堆積高度隨時間的推移而不斷增加，由此導(dǎo)致海冰內(nèi)力不斷增強(qiáng)。當(dāng)海冰內(nèi)力與風(fēng)、流的拖曳力達(dá)到平衡時，海冰單元運(yùn)動趨于靜止，堆積高度達(dá)到穩(wěn)定。在海冰漂移和堆積過程中，海冰單元的速度、尺寸和厚度均發(fā)生改變。

圖3 變寬度水道內(nèi)海冰的初始分布Fig.3 Initial distribution of ice cover in a various-width channel

表1 海冰動力過程的改進(jìn)離散元模擬中主要計算參數(shù)

Tab.1 Main computational parameters in MDEM simulation of sea ice dynamic processes

參數(shù)定義數(shù)值參數(shù)定義數(shù)值E彈性模量/MPa100μp海冰單元間摩擦系數(shù)0 5Ca風(fēng)的拖曳系數(shù)0 0015Cw流的拖曳系數(shù)0 0045Va風(fēng)速/m·s-120 0Vw流速/m·s-10 0

續(xù)表1

圖4 不同時刻的海冰速度和平均冰厚分布Fig.4 The distribution of ice velocity and average ice thickness at different time

以上結(jié)果表明，海冰單元在風(fēng)和流的驅(qū)動下向右側(cè)漂移過程中，海冰的速度、厚度不斷變化。特別在第1 d繞過三角形障礙物過程中，海冰在障礙物附近的速度明顯減小，冰厚也相對較厚。在障礙物后方形成一個明顯的水域，海冰呈現(xiàn)顯著的繞流特性。此外，受水道邊界摩擦的影響，邊界附近的冰速明顯減小，其平均冰厚也小于水道中部冰厚。

為進(jìn)一步驗證模型對冰厚模擬的準(zhǔn)確性，這里將海冰穩(wěn)定后的平均厚度與解析解進(jìn)行對比。本文海冰穩(wěn)定條件是指海冰顆粒單元的速度的最大值接近于零?？紤]邊界摩擦?xí)r，海冰的穩(wěn)態(tài)堆積高度可由經(jīng)典的冰壩理論計算[27]。它將海冰應(yīng)力在水道寬度方向進(jìn)行平均，并通過求解風(fēng)和流的拖曳力、海冰內(nèi)力以及邊界摩擦力的平衡方程確定海冰的堆積高度，即[14]：

(9)

式中，Hi0為海冰的初始厚度，B為水道的寬度，這里參照圖3知B=W2。

將以上水道內(nèi)海冰厚度的離散元模擬結(jié)果與式(9)的解析解列于圖5中，可以發(fā)現(xiàn)平均冰厚與解析解有較高的擬合度，由此驗證了改進(jìn)離散元模擬在海冰動力學(xué)模擬中的準(zhǔn)確性。

圖5 海冰平均厚度的離散元計算值與解析解的對比分析Fig.5 Comparison of average ice thickness simulated with DEM to analytical solution

3.2 旋轉(zhuǎn)風(fēng)場作下海冰動力過程

為檢驗海冰動力學(xué)數(shù)值方法的有效性，F(xiàn)lato建立了經(jīng)典的旋轉(zhuǎn)風(fēng)場作用下的海冰運(yùn)動算例[15]。在500 km×500 km的計算域內(nèi)，風(fēng)場中心位于(250 km，220 km)處，海冰初始分布于計算域的上半?yún)^(qū)域，如圖6所示。海冰初始厚度Hi0=1.0 m，初始密集度Ni0=80%，海冰單元的初始直徑D0=4.0 km。計算域中位置r處的風(fēng)速為：

(10)

式中，W為風(fēng)速矢量，K為垂直水面的單位向量，r和r分別為計算域中任一點的位置矢量與其到風(fēng)場中心的距離，風(fēng)速參數(shù)ω=0.5×10-3rad/s，λ=8×105m2/s。

圖6 旋轉(zhuǎn)風(fēng)場及海冰初始分布Fig.6 The initial distribution of sea ice in rotational wind field

采用改進(jìn)離散元模型對旋轉(zhuǎn)風(fēng)場作用下的海冰動力過程進(jìn)行了10 d數(shù)值計算，第2 d、5 d、10 d的冰速和平均冰厚分布如圖7所示?？梢园l(fā)現(xiàn)，海冰在風(fēng)和流的拖曳作用下在旋轉(zhuǎn)風(fēng)場內(nèi)圍繞風(fēng)場中心做渦旋運(yùn)動。在海冰運(yùn)動過程中，左側(cè)邊界附近的海冰單元與邊界發(fā)生擠壓和堆積使平均冰厚增加，而右側(cè)邊界附近的海冰單元因遠(yuǎn)離邊界漂移使平均冰厚減??；此外，從平均冰厚的等值線分布可以發(fā)現(xiàn)，在渦旋中心兩側(cè)會形成兩個明顯的渦，并在海冰漂移過程中發(fā)生移動和變形。以上計算結(jié)果表明，改進(jìn)的離散元模型可對旋轉(zhuǎn)風(fēng)場作用下海冰的動力演化過程進(jìn)行精確數(shù)值模擬。

4 渤海海冰動力過程的數(shù)值模擬

為進(jìn)一步驗證改進(jìn)離散元模型在海冰動力學(xué)模擬中的適用性，下面將其用于渤海海冰動力演化過程的數(shù)值分析，并通過海冰衛(wèi)星遙感資料和油氣作業(yè)區(qū)海冰監(jiān)測數(shù)據(jù)對計算結(jié)果進(jìn)行檢驗。

4.1 遼東灣海冰分布的演化

在海冰動力學(xué)的MDEM模擬中，氣象條件采用遼東灣JZ20-2平臺上實測的時間間隔為10 min的風(fēng)速，水文條件采用Leendertse二維非線性長波模式在歐拉坐標(biāo)計算得到的潮流場。

海冰初始冰厚和密集度可由油氣平臺的海冰現(xiàn)場監(jiān)測圖像和衛(wèi)星遙感圖像提取[45—47]。這里對2010年1月22日10時50分的衛(wèi)星遙感圖像進(jìn)行海冰的參數(shù)提取。海冰衛(wèi)星遙感圖像如圖8a所示，由此提取的海冰單元分布如圖8b所示，其顏色定性地表示海冰單元的局部厚度。依據(jù)遼東灣JZ20-2油氣平臺上實測的海冰厚度，由海冰單元參數(shù)確定的海冰厚度和密集度分布如圖8b～d所示。在海冰動力學(xué)模擬中，海冰單元的初始直徑為2.0 km，共有5 618個單元，其他主要計算參數(shù)取表1中的數(shù)值。通過采用改進(jìn)離散單元模型對海冰動力過程的48 h數(shù)值模擬，海冰單元位置、海冰密集度、冰厚和冰速分布分別如圖9和10所示。

從模擬結(jié)果可以發(fā)現(xiàn)，海冰單元位置分布規(guī)律與衛(wèi)星遙感資料基本一致。在數(shù)值模擬的48 h內(nèi)，風(fēng)速在2～8 m/s之間，氣溫在-5～-8℃之間，氣象條件變化平緩，海冰運(yùn)動主要受控于潮流，且冰厚受熱力因素影響較小。對于海冰的動力過程，MDEM方法較好地模擬了遼東灣海冰的動力演化規(guī)律。它通過對海冰單元在拉格朗日坐標(biāo)下的運(yùn)移、碰撞和重疊特性的離散元分析，通過改變單元的位置、尺度、密集度和冰厚等海冰參數(shù)，較精確地描述海冰的運(yùn)移過程。特別是，它通過對海冰單元漂移速度和運(yùn)動軌跡的數(shù)值計算，能較準(zhǔn)確地確定冰緣線位置，而不存在歐拉坐標(biāo)下的數(shù)值擴(kuò)展現(xiàn)象。

但我們也發(fā)現(xiàn)模擬的冰緣線面積較實際冰緣線面積略大。這主要是因為由拉格朗日坐標(biāo)下海冰單元向歐拉坐標(biāo)下海冰參數(shù)的等值線插值過程中，有一個向外擴(kuò)散的趨勢。此外，本文采用MDEM方法對海冰演化過程中尚未考慮熱力因素的影響，還需在后續(xù)工作中進(jìn)一步完善以更合理地模擬海冰的生消和運(yùn)移過程。

圖7 旋轉(zhuǎn)風(fēng)場作用下海冰速度和平均冰厚分布Fig.7 Distributions of ice velocity and average ice thickness in the rotational wind field

圖9 海冰衛(wèi)星遙感圖像及MDEM模擬的海冰位置Fig.9 Satellite images of sea ice and the simulated location of sea ice elements with MEDM

圖10 MDEM模擬的海冰密集度、厚度和速度分布Fig.10 Distribution of ice concentration，thickness and velocity simulated with MDEM

4.2 遼東灣JZ20-2油氣海域的海冰參數(shù)演化

在渤海海冰動力學(xué)模擬中，油氣作業(yè)區(qū)海冰參數(shù)的確定對保障油氣安全作業(yè)具有很強(qiáng)的工程意義，也是冰期油氣作業(yè)中海冰管理工作的重要組成部分。在以上遼東灣海冰動力過程的改進(jìn)離散元模擬中，可以對JZ20-2油氣作業(yè)區(qū)(40°30′N，121°21′E)的海冰參數(shù)進(jìn)行提取。這里采用SPH方法中的Gauss函數(shù)由該油氣作業(yè)區(qū)的海冰單元插值得到該海域的海冰厚度、密集度、速度等參數(shù)。圖11為JZ20-2油氣作業(yè)區(qū)附近一個海冰單元的漂移軌跡，可以發(fā)現(xiàn)其受往復(fù)潮流和西北向風(fēng)的拖曳作用，在往復(fù)運(yùn)動過程中向西南方向漂移。圖12和圖13給出了該海域插值得到的海冰速度和厚度在48 h內(nèi)的變化過程，實測值由該油氣平臺上的海冰現(xiàn)場監(jiān)測及數(shù)字圖像處理系統(tǒng)獲得[45]?？梢园l(fā)現(xiàn)，海冰在半日潮作用下往復(fù)運(yùn)動，其厚度保持相對穩(wěn)定。結(jié)果表明，MDEM方法對該海域海冰參數(shù)的模擬結(jié)果與現(xiàn)場監(jiān)測數(shù)值基本一致。

由于受海冰監(jiān)測資料的限制，本文尚未對大風(fēng)、寒流條件下的極端天氣條件進(jìn)行MDEM數(shù)值模擬。這在后續(xù)研究中將重點關(guān)注不同天氣條件下的海冰動力演化規(guī)律，對海冰的重疊、堆積過程進(jìn)行數(shù)值分析和實測數(shù)據(jù)驗證，進(jìn)一步檢驗MDEM方法對海冰動力過程的適用性。

5 結(jié)論

針對海冰的離散分布特點以及海冰漂移、重疊中的形變特征，本文發(fā)展了海冰動力過程的改進(jìn)離散單元模型。它在拉格朗日坐標(biāo)下將海冰離散為多個具有密集度、厚度的海冰單元，通過海冰單元間的摩擦、碰撞計算海冰內(nèi)力。海冰單元間采用具有Mohr-Coulomb摩擦準(zhǔn)則的黏彈性接觸模型，并考慮海冰單元在法向力作用下的塑性變形。在海冰單元質(zhì)量守恒條件下，計算單元相應(yīng)的尺寸、密集度和厚度，由此獲得海冰動力過程中的形變規(guī)律和重疊堆積特征。為驗證該模型的可靠性，對變寬度水道內(nèi)海冰的漂移和堆積過程進(jìn)行了數(shù)值模擬，冰厚分布與解析值相吻合；通過對旋轉(zhuǎn)風(fēng)場內(nèi)海冰漂移過程的數(shù)值分析，進(jìn)一步驗證了該模型對描述海冰動力學(xué)特性的精確性。最后，采用改進(jìn)的離散單元模型對渤海海冰的動力過程進(jìn)行了48 h數(shù)值模擬，計算結(jié)果與衛(wèi)星遙感資料和油氣作業(yè)區(qū)海冰現(xiàn)場監(jiān)測數(shù)據(jù)相吻合。

海冰動力過程的改進(jìn)離散單元模型具有物理意義明確、計算簡單、精度高等優(yōu)點，可對海冰在重疊堆積過程中產(chǎn)生的塑性變形進(jìn)行有效的計算。在本文工作基礎(chǔ)上將進(jìn)一步考慮海冰的動力輻散現(xiàn)象，并針對熱力因素的影響建立海冰單元的凍結(jié)和斷裂模型，更好地應(yīng)用于海冰的動力-熱力過程的數(shù)值分析。

圖11 遼東灣JZ20-2海域海冰單元漂移軌跡Fig.11 Drifting trajectory of one sea ice element at the JZ20-2 sea area of Liaodong Bay

圖12 遼東灣JZ20-2海域48 h模擬和實測冰速Fig.12 Simulated and measured ice velocity at the JZ20-2 sea area of Liaodong Bay for 48 hours

圖13 JZ202海域48 h內(nèi)的模擬和實測冰厚Fig.13 Simulated and measured ice thickness at the JZ20-2 sea area of Liaodong Bay for 48 hours

[1] 季順迎，岳前進(jìn). 工程海冰數(shù)值模型及應(yīng)用[M]. 北京：科學(xué)出版社，2011．

Ji Shunying，Yue Qianjin. Engineering Sea Ice Numerical Model and its Application[M]. Beijing: Science Press，2011．

[2] Holland M M，Stroeve J. Changing seasonal sea ice predictor relationships in a changing Arctic climate[J]. Geophysical Research Letters，2011，38: L18501．

[3] Dempsey J P. Research trends in ice model[J]. International Journal of Solids and Structures，2000，37(2): 131-153．

[4] Lepparanta M，Lensu M，Lu Q M. Shear flow of sea ice in the Marginal Ice Zone with collision theology[J]. Geophysica，1990，25(1/2):57-74．

[5] Tremblay L B，Mysak L A. Modeling sea ice as a granular material，including the dilatancy effect[J]. Journal of Physical Oceanography，1997，27(2): 2342-2360．

[6] OverLand J E，McNutt S L，Salo S，et al. Arctic sea ice sa a granular plastic[J]. Journal of Geophysical Research，1998，103(C10): 21845-21868．

[7] Hibler W D. Sea ice fracturing on the large scale[J]. Engineering Fracture Mechanics，2001，68(4): 2013-2043．

[8] Schulson E M. Compressive shear fault within arctic sea ice: Fracture on scale large and small[J]. Journal of Geophysical Research，2004，109(C07016): 1-23．

[9] Hibler W D. A Dynimic Thermodynamic sea ice model[J]. Journal of Geophysical Oceanography，1979，9(3):817-846．

[10] Zhang Z H，Lepparanta M. Modeling the influence of ice on sea level variations in the Baltic Sea[J]. Journal of Geophysical Research，1995，31(2): 31-45．

[11] 吳輝碇. 海冰的動力-熱力過程的數(shù)學(xué)處理[J]. 海洋與湖沼，1991，20(4): 321-327．

Wu Huiding. Mathematic representations of sea ice dynamic thermodynamic processes[J]. Oceanologia et Limnologia Sinica，1991，20(4): 321-327．

[12] 蘇潔，吳輝碇，白珊，等. 渤海冰-海洋耦合模式:Ⅱ.個例試驗[J]. 海洋學(xué)報，2005，27(2): 18-27．

Su Jie，Wu Huiding，Bai Shan，et al. A coupled ice-ocean model for the Bohai Sea:Ⅱ. Case study[J]. Haiyang Xuebao，2005，27(2): 18-27．

[13] Gutfraind R，Savage S B. Smoothed Particle Hydrodynamics for the simulation of broken-ice field: Mohr-Coulomb-Type rheology and frictional boundary conditions[J]. Journal of Computational Physics，1997，134(3): 203-215．

[14] 季順迎，沈洪道，王志聯(lián)，等. 基于Mohr-Coulomb準(zhǔn)則的黏彈塑性海冰動力學(xué)本構(gòu)模型[J]. 海洋學(xué)報，2005，27(4): 19-30．

Ji Shunying，Shen Hongdao，Wang Zhilian，et al. A viscoelastic-plastic constitutive model with Mohr-Coulomb yielding criterion for sea ice dynamics[J]. Haiyang Xuebao，2005，27(4): 19-30．

[15] Flato G M. A particle-in-cell sea-ice model[J]. Atmosphere and Oceanography，1993，31(3): 339-358．

[16] Huang Z J，Savage S B.Particle-in-cell and finite difference approaches for the study of marginal ice zone problems[J]. Cold Regions Science and Technology，1998，28(1):1-28．

[17] 劉煜，吳輝碇，張占海，等. 基于質(zhì)點-網(wǎng)格模式的海冰厚度變化過程數(shù)值模擬[J]. 海洋學(xué)報，2006，28(2): 14-21．

Liu Yu，Wu Huiding，Zhang Zhanhai，et al. Modeling for the dynamic process of ice thickness variation using a particle-in-cell ice model[J]. Haiyang Xuebao，2006，28(2): 14-21．

[18] Kubat I，Sayed M，Savage S，et al. Numerical simulations of ice thickness redistribution in the Gulf of St. Lawrence[J]. Cold Regions Science and Technology，2010，60(2): 15-28．

[19] Zhang J，Rothrock D A. Effect of sea ice rheology in numerical investigations of climate[J]. Journal of Geophysical Research，2005，110: C08014．

[20] Hutchings J K，Jasak H，Laxon S W. A strength implicit correction scheme for the viscous-plastic sea ice model[J]. Ocean Modelling，2004，7(2): 111-133．

[21] Wang L R，Ikeda M. A Lagrangian description of sea ice dynamics using the finite element method[J]. Ocean Modelling，2004，7(3): 21-38．

[22] Sulsky D，Schreyer H，Peterson K，et al. Using the material-point method to model sea ice dynamics[J]. Journal of Geophysical Research，2007，112: C02S90．

[23] 李海，季順迎，沈洪道，等. 海冰動力學(xué)的混合拉格朗日-歐拉數(shù)值方法[J]. 海洋學(xué)報，2008，30(2):1-11．

Li Hai，Ji Shunying，Shen Hongdao，et al. A hybrid Lagrangian-Eulerian numerical model for sea ice dynamics[J]. Haiyang Xuebao，2008，30(2):1-11．

[24] Lemieux J F，Tremblay B，Sedlcek J，et al. Improving the numerical convergence of viscous-plastic sea ice models with the Jacobian-free Newton Krylov method[J]. Journal of Computational Physics，2010，229: 2840-2852．

[25] Tsamados M，F(xiàn)eltham D L，Wilchinsky A V. Impact of a new anisotropic rheology on simulations of Arctic sea ice[J]. Journal of Geophysical Research，2013，doi:10.1029/2012JC007990．

[26] Coon M D，Knoke G S，Echert D C，et al. The architecture of an anisotropic elastic-plastic sea ice mechanics constitutive law[J]. Journal of Geophysical research，1998，103(C10): 21915-21925．

[27] Hunk E C，Dukowicz J K. An elastic-viscous-plastic model for sea ice dynamics[J]. Journal of Physical Oceanography，1997，27(2): 1849-1867．

[28] Losch M，Danilov S. On solving the momentum equations of dynamic sea ice models with implicit solvers and the elastic-viscous-plastic technique[J]. Ocean Modelling，2012，41(4): 42-52．

[29] 王剛，岳前進(jìn)，李海，等. 基于SPH方法的渤海海冰動力學(xué)數(shù)值模擬[J]. 大連理工大學(xué)學(xué)報，2007，47(3): 323-328．

Wang Gang，Yue Qianjin，Li Hai，et al. Numerical simulation of sea ice dynamics with Sph approach in Bohai Sea[J]. Journal of Dalian University of Technology，2007，47(3): 323-328．

[30] Liu M B，Liu G R. Smoothed Particle Hydrodynamics(SPH): an overview and recent developments[J]. Archives of Computational Methods in Engineering，2010，17(3):25-76．

[31] Dai M，Shen H H，Hopkins M A，et al. Wave rafting and the equilibrium pancake ice cover thickness[J]. Journal of Geophysical Research，2004，109: C07023．

[32] Herman A. Influence of ice concentration and floe-size distribution on cluster formation in sea-ice floes[J]. Central European Journal of Physics，2012，10(3): 715-722．

[33] Wilchinsky A V，F(xiàn)eltham D L. Rheology of discrete failure regimes of anisotropic sea ice[J]. Journal of Physical Oceanography，2012，42(4): 1065-1082．

[34] 季順迎，李春花，劉煜. 海冰離散元模型的研究回顧及展望[J]. 極地研究，2012，24(4): 315-330．

Ji Shunying，Li Chunhua，Liu Yi. A review of advances in sea-ice discrete element models[J]. Chinese Journal of Polar Research，2012，24(4): 315-330．

[35] Rothrock D，Thomdike A S. Measuring the sea ice grian size distribution[J]. Journal of Geophysical Research，1984，89(C4):6477-6486．

[36] Shen H H，Hibler W D，Lepparanta M. On applying granular flow theory to a deforming broken ice field[J]. Acta Mechanics，1986，63(3): 143-160．

[37] Wichinsky A V，F(xiàn)eltham D L. Anisotropic model for granulated sea ice dynamics[J]. Journal of the Mechanics and Physics of Solids，2006，54(4):1147-1185．

[38] Sedlacek J，Lemieux J F，Mysak L A，et al. The granular sea ice model in spherical coordinates and its application to a global climate model[J]. Journal of Climate，2007，20(1):5946-5961．

[39] Hopkins M A，F(xiàn)rankenstein S，Thorndike A S. Formation of an aggregate scale in Arctic sea ice[J]. Journal of Geophysical Research，2004，109: C01032．

[40] Hopkins M A，Thorndike A S. Floe formation in Arctic sea ice[J]. Journal of Geophysical Research，2006，111: C11S23．

[41] Shen H T，Shen H H，Tsai S M. Dynamic transport of river ice[J]. Journal of Hydraulic research，1990，28(6): 659-671．

[42] Babic M，Shen H H，Shen H T. The stress tensor in granular shear flows of uniform，deformable disks at high solids concentrations[J]. Journal of Fluid Mechanics，1990，219(4): 81-118．

[43] Pariset E，Hausser R，Gagnon A. Formation of ice cover and ice jams in rivers[J]. Journal of Hydraulics Division，ASCE，1966，92(HY6):1-24．

[44] 沈洪道. 冰動力學(xué)的拉格朗日離散元模式[J]. 海洋預(yù)報，1999，16(3): 71-84．

Shen Hongdao. Lagrangian discrete-parcel model for ice dynamics[J].Marine Forecasts，1999，16(3): 71-84．

[45] 季順迎，王安良，王宇新，等. 渤海海冰現(xiàn)場監(jiān)測的數(shù)字圖像技術(shù)及其應(yīng)用[J]. 海洋學(xué)報，2011，33(4): 79-87．

Ji Shunying，Wang Anliang，Wang Yuxin，et al. A digital image technology and its application for the sea ice field observation in the Bohai Sea[J]. Haiyang Xuebao，2011，33(4): 79-87．

[46] Kwok R，Cunningham G F. ICES at over Arctic sea ice: Estimation of snow depth and ice thickness[J]. Journal of Geophysical Research，2008，113:C08010．

[47] Sun B，Wen J，He M，et al. Sea ice thickness measurement and its underside morphology analysis using radar penetration in the Arctic Ocean[J]. Science China，2003，46(11): 1151-1160．

Modified discrete element model for sea ice dynamics and its applications in the Bohai Sea

Ji Shunying1，Wang Anliang2，Mi Lili1，Liu Yu2，Li Baohui2

(1.StateKeyLaboratoryofStructuralAnalysisforIndustrialEquipment，DalianUniversityofTechnology，Dalian116023，China; 2.KeyLaboratoryofResearchonMarineHazardsForecasting，StateOceanicAdministration,NationalMarineEnvironmentalForecastingCenter，Beijing100081，China)

Breakup，rafting and ridging of ice cover exists widely in the polar and sub-polar regions. These processes affact the growth，vanishing and drifting of sea ice significantly. Considering the discrete distribution of sea ice on various scales，a discrete element model (DEM) should be developed to improve the sea ice numerical model and its computational precission. Thus，a modified discrete element model (MDEM) is established in this study to simulate the sea ice dynamics. Different with the traditional DEM，the ice cover is subdivided into a series of disks with their own characteristics including thickness，velocity，size and concentration，adopting the concept of smoothed particle hydrodynamics (SPH) of sea ice dynamics. Each sea ice element which is an assembly of ice floes changes in its size，concentration and thickness，according to the mass conservation law during drifting and inter-element collisions. According to the non-continuous distribution and rheology characteristics of ice cover，the viscous-elastic constitutive model is adopted. And the Mohr-Coulomb friction law is considered to determine the plastic deformation and tangential friction. To assess the reliability of this MDEM for sea ice dynamics，the drifting and ridging of ice cover in a various-width channel is simulated，and the simulated distribution of ice thickness is validated by the analytical solution. The drifting of sea ice in a rotational wind field is also simulated efficiently with high precision. Moreover，the sea ice dynamics in the Bohai Sea is simulated for 48 h. The simulated results match well with the satellite remote images and field observed data. In the future study，the MDEM will be improved by coupling dynamics and thermodynamics of sea ice. The growth，vanishing and drifting of sea ice will be simulated more accurately by consideringthe refrozen effect and breakage feature of ice cover.

sea ice dynamics; discrete element model; ice ridge; Bohai Sea

10.3969/j.issn.0253-4193.2015.05.006

2013-04-28；

2015-01-11。

國家自然科學(xué)基金項目(41176012);國家海洋公益性行業(yè)科研專項經(jīng)費(fèi)項目(201105016, 201205007);高等學(xué)校博士學(xué)科點專項科研基金項目(20130041110010)。

季順迎(1972年—)，男，河北省武邑縣人，博士，教授，從事工程海冰數(shù)值模型研究。E-mail:jisy@dlut.edu.cn

P731.15

A

0253-4193(2015)05-0054-14

季順迎，王安良，米麗麗，等. 海冰動力過程的改進(jìn)離散元模型及在渤海的應(yīng)用[J]. 海洋學(xué)報，2015，37(5)：54-67，

Ji Shunying，Wang Anliang，Mi Lili，et al.Modified discrete element model for sea ice dynamics and its applications in the Bohai Sea[J]. Haiyang Xuebao，2015，37(5)：54-67，doi：10.3969/j.issn.0253-4193.2015.05.006