海冰動(dòng)力學(xué)數(shù)值方法研究進(jìn)展

1、第19卷第6期2004 年12月1001-8166 0963-08 文章編號(hào)：（2004）06-地球科學(xué)進(jìn)展A DVAN CE S I N E AR TH S C I E N C EV o l .19 N o.6D e c. ，2004海冰動(dòng)力學(xué)數(shù)值方法研究進(jìn)展季順迎，岳前進(jìn)，王瑞學(xué)）（大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧大連116023摘要：在海冰動(dòng)力學(xué)數(shù)值模擬和預(yù)測(cè)研究中，人們將海冰視為連續(xù)介質(zhì)分別建立了歐拉坐標(biāo)下S PH ）方法、的有限差分（FD ）方法、拉格朗日坐標(biāo)下的光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)（ 歐拉和拉格朗日坐標(biāo)P I C ）相結(jié)合的質(zhì)點(diǎn)網(wǎng)格法（ ，近年來(lái)又發(fā)展了基于非連續(xù)介質(zhì)

2、的顆粒流（G F ）方法。對(duì)以上幾種海FD 、P I C 和 S PH 方法可適用于中長(zhǎng)期海冰動(dòng)力學(xué)數(shù)值方法的特點(diǎn)和適用性進(jìn)行了討論，結(jié)果表明：S PH 方法的計(jì)算效率需進(jìn)一步提高；G F 方法在不同尺度下的海冰動(dòng)力學(xué)數(shù)冰動(dòng)力學(xué)數(shù)值模擬，但值模擬中的計(jì)算精度均有很強(qiáng)的適用性，但目前較適用于小尺度下海冰動(dòng)力學(xué)基本特性的數(shù)值試驗(yàn)研究，計(jì)算時(shí)效還不能滿(mǎn)足實(shí)際海冰數(shù)值模擬和預(yù)測(cè)的要求。為進(jìn)一步提高海冰動(dòng)力學(xué)模擬的精度和適用性，在不同時(shí)空尺度下分別發(fā)展與其相適應(yīng)的數(shù)值方法是必要的。關(guān)鍵詞：海冰動(dòng)力學(xué)；有限差分法；光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)；質(zhì)點(diǎn)網(wǎng)格法；顆粒流方法P 731.15文獻(xiàn)標(biāo)識(shí)碼：A 中圖分類(lèi)號(hào)：述54

3、，0前言為研究極區(qū)海冰對(duì)全球氣候的影響，也為解決中緯度冰區(qū)航運(yùn)、資源開(kāi)發(fā)中的海冰問(wèn)題，人們對(duì)不同尺度下海冰的動(dòng)力學(xué)特性進(jìn)行了大量的實(shí)驗(yàn)研究和數(shù)值模擬。在風(fēng)、浪、流和潮汐等海洋環(huán)境條件的作用下，海冰的動(dòng)力過(guò)程非常復(fù)雜，在不同尺度下的分布特性也有很大的差異。對(duì)不同尺度下的海冰動(dòng)力學(xué)過(guò)程進(jìn)行精確地?cái)?shù)值模擬和預(yù)測(cè)，建立和完善合理有效的數(shù)值方法一直是海冰動(dòng)力學(xué)研究的一個(gè)重要問(wèn)題。從20世紀(jì)70年代以來(lái)，人們?cè)趯?duì)大、中尺度下海冰動(dòng)力演化過(guò)程的數(shù)值模擬和預(yù)測(cè)中，大多將海冰視為二維流體或流變體，并采用歐拉坐標(biāo)下的有限差分法進(jìn)行計(jì)算，其模擬結(jié)果可較好地表征海冰31動(dòng)力演化的基本規(guī)律。但人們也注意到，在用有限差分

4、法計(jì)算海冰動(dòng)量方程中的平流項(xiàng)時(shí)有明顯的數(shù)值擴(kuò)散現(xiàn)象。這不僅限制了其準(zhǔn)確預(yù)測(cè)冰邊緣線的運(yùn)動(dòng)，對(duì)冰區(qū)的動(dòng)力特性也不能很好地描2003-06-16 ；2004-02-25收稿日期：修回日期：。人們?cè)诓捎脷W拉差分法于是，從90年代初期，進(jìn)行海冰數(shù)值模擬的同時(shí)，也在嘗試?yán)窭嗜辗椒ā?P r it c h a r d 等使用了一種自適應(yīng)的拉格朗日網(wǎng)如格單元進(jìn)行冰邊緣線的模擬；基于光滑質(zhì)點(diǎn)流體動(dòng)S h e n 力學(xué)方法，發(fā)展了一種純拉格朗日離散單D PM ）元模型（ ，較好地解決了數(shù)值擴(kuò)散問(wèn)題，并用N i aga r a 來(lái)模擬尼亞加拉（ ）河冰輸運(yùn)以及渤海海冰的動(dòng)力學(xué)過(guò)程。針對(duì)冰邊緣線復(fù)雜，且冰水兩相介

5、4F l a t o 將歐拉坐標(biāo)和拉格朗日坐標(biāo)質(zhì)共存的特點(diǎn)，P I C ）引入到海冰動(dòng)力學(xué)模擬耦合的質(zhì)點(diǎn)網(wǎng)格法（7，85，中，并在極區(qū)和中緯度海域的海冰數(shù)值模擬中得到了較為廣泛的應(yīng)用，其模擬精度明顯優(yōu)于有限差分9R h ee m 等發(fā)展了一種計(jì)算中尺度海冰動(dòng)力學(xué)法；D M D F ），的分布質(zhì)量/離散塊體模型（ 并用于鄂霍次O kh o t s k 克（ ）海的海冰動(dòng)力過(guò)程模擬。以上海冰動(dòng)力學(xué)模擬的數(shù)值方法大多是考慮大、中尺度下海冰的連續(xù)性分布而建立的。然而，在B a lti c 冰邊緣區(qū)、渤海和波羅的（ ）海等中緯度海域，*基金項(xiàng)目：40206004 國(guó)家自然科學(xué)基金項(xiàng)目“中小尺度海冰動(dòng)力學(xué)

6、本構(gòu)模型及數(shù)值方法研究”（編號(hào)：）資助.1972-E - m a i l j i s y d l ut . e du. c n 作者簡(jiǎn)介：季順迎（），男，河北省武邑人，講師，主要從事工程海冰數(shù)值模擬和顆粒流體動(dòng)力學(xué)研究. ： 649地球科學(xué)進(jìn)展第19卷海冰在強(qiáng)烈的風(fēng)、浪作用下表現(xiàn)出明顯的區(qū)域性和間斷性分布特點(diǎn)。為此，S h e n 等10，11將顆粒流理論用于模擬海冰與波浪的動(dòng)力作用過(guò)程，并建立了海冰碰撞流變學(xué)。隨后，H o pk i ns 12，13經(jīng)進(jìn)一步完善后，應(yīng)用于海冰重疊和堆積的數(shù)值模擬，并通過(guò)現(xiàn)場(chǎng)觀測(cè)和室內(nèi)實(shí)驗(yàn)對(duì)模擬結(jié)果進(jìn)行了驗(yàn)證。海冰動(dòng)力學(xué)數(shù)值方法是與海冰本構(gòu)模型、海冰熱力學(xué)等其

7、它問(wèn)題的研究同步開(kāi)展的，并相互促進(jìn)。本文將主要對(duì)目前海冰動(dòng)力學(xué)所采用的數(shù)值方法進(jìn)行討論，進(jìn)而對(duì)海冰動(dòng)力學(xué)數(shù)值方法的進(jìn)一步發(fā)展進(jìn)行分析和展望。1海冰動(dòng)力學(xué)的有限差分法在海冰動(dòng)力學(xué)數(shù)值模擬中，發(fā)展最早、應(yīng)用范圍最廣的是歐拉坐標(biāo)下的有限差分法。有限差分計(jì)算模式一般包括海冰的動(dòng)量方程、連續(xù)方程轉(zhuǎn)換的冰厚和密集度守恒方程和海冰本構(gòu)模型。1.1傳統(tǒng)的 H i b l e r 有限差分方法為研究北極大尺度下海冰的演化過(guò)程，H i b l e r 1建立了粘塑性海冰動(dòng)力學(xué)本構(gòu)模型，并采用有限差分方法進(jìn)行數(shù)值模擬。該有限差分法中的海冰變量采用如圖1所示的Ar a k a w s - B 網(wǎng)格分布。圖中u 和v

8、分別為x 和y 方向冰速，m 為單位面積海冰質(zhì)量，h 和A 分別為海冰厚度和密集度，P 為冰內(nèi)壓力項(xiàng)，和分別為海冰塊體和剪切粘性系數(shù)，為海冰應(yīng)變率，D 1和D 2分別為計(jì)算中的一階和二階擴(kuò)散系數(shù)。圖1H i b l e r （1979 ）有限差分法中的網(wǎng)格模式F i g.1 G r i d sc h e m e i n t h e F D m e t h od o f H i b l e r （1979 ）計(jì)算中采用隱式的超松馳迭代方法求解動(dòng)力方程以確定冰速，用修正的半隱式La x- W e nd r o f f 方法求解連續(xù)方程以確定冰厚和密集度分布，進(jìn)而可采用較長(zhǎng)的時(shí)間步長(zhǎng)。在Hi b l

9、 e r 1對(duì)北極海冰的數(shù)值模擬中，x 和y 方向的網(wǎng)格尺度均為125 k m ，時(shí)間步長(zhǎng)為24h ，開(kāi)闊水和海冰的臨界冰厚h 00.5m。經(jīng)實(shí)測(cè)資料驗(yàn)證，該計(jì)算方法能夠很好地模擬北極海冰的動(dòng)力演化過(guò)程。H i b l e r 的有限差分模式還被進(jìn)一步應(yīng)用于冰緣區(qū)（M I Z ）、波羅的海和渤海的中尺度海冰動(dòng)力學(xué)模擬中2，3，1416 。在對(duì)冰邊緣區(qū)（MI Z ）和東格陵蘭（G r ee n l a nd ）海的數(shù)值模擬中，L e pp a r a n t a 等2，14縮小了計(jì)算的時(shí)空尺度，并將Hi b l e r 原有的差分格式引入到拉格朗日網(wǎng)格中，以降低平流項(xiàng)計(jì)算的擴(kuò)散現(xiàn)象。在波羅的海的

10、海冰模擬中，網(wǎng)格尺度被取為18.5 k m 和5 k m 。在渤海海冰數(shù)值模擬中，將海冰模式與數(shù)值天氣預(yù)報(bào)模式、大氣邊界層模式和潮流模式聯(lián)結(jié)，并采用密集度和冰厚的形變函數(shù)來(lái)處理計(jì)算區(qū)域內(nèi)開(kāi)闊水、平整冰和堆積冰的演化過(guò)程3，17，18。在每計(jì)算步對(duì)2個(gè)時(shí)間層進(jìn)行積分時(shí)，在第二個(gè)時(shí)間層加入了二階和四階擴(kuò)散項(xiàng)以增強(qiáng)計(jì)算的穩(wěn)定性。為縮短計(jì)算時(shí)間，在開(kāi)邊界采用了一種移動(dòng)邊界技術(shù)。計(jì)算中的網(wǎng)格距為8.64k m ×11.11 k m ，時(shí)間步長(zhǎng)為3h 或1h ，臨界冰厚h 02c m 。通過(guò)衛(wèi)星遙感和現(xiàn)場(chǎng)監(jiān)測(cè)資料驗(yàn)證，該模式已成功地用于渤海海冰業(yè)務(wù)數(shù)值預(yù)報(bào)17，18。另外，在渤海早期的海冰數(shù)值模

11、擬中，王仁樹(shù)等19也是采用有限差分法，并用交錯(cuò)網(wǎng)格模式和線粘性本構(gòu)模型，空間步長(zhǎng)為20k m ，時(shí)間步長(zhǎng)為24h 。為研究潮流對(duì)海冰漂移的影響，又將時(shí)間步長(zhǎng)取為15mi n 。經(jīng)與實(shí)測(cè)資料對(duì)比，該模式能夠反映海冰漂移過(guò)程的基本特征20。1.2 H i b l e r 有限差分方法的改進(jìn)在研究極區(qū)海冰對(duì)全球氣候的影響時(shí)，需要進(jìn)行長(zhǎng)期海冰數(shù)值模擬。當(dāng)采用Hi b l e r 的粘塑性本構(gòu)進(jìn)行海冰動(dòng)力學(xué)有限差分計(jì)算時(shí)，如果采用顯式計(jì)算，則要取用很小的時(shí)間步長(zhǎng)，進(jìn)而影響計(jì)算時(shí)效；如果采用隱式超松馳迭代法，可取較大的時(shí)間步長(zhǎng)，但需要大量的迭代計(jì)算，受網(wǎng)格尺度約束也較大，更不能適用于并行機(jī)上進(jìn)行計(jì)算。I

12、p 等21在研究基于Mo hr - C o u l o m b 屈服準(zhǔn)則的粘塑性海冰動(dòng)力學(xué)本構(gòu)模型時(shí)，將H i b l e r 的 A r a ka- w a- B 型網(wǎng)格用 A r a k a wa- C 型取代。它是將網(wǎng)格角點(diǎn)上的x 和y 向冰速分別置于網(wǎng)格邊界的左右和上下中心處。這不但簡(jiǎn)化了冰內(nèi)應(yīng)力的計(jì)算，同時(shí)也消除了原B 型格式在計(jì)算壓力場(chǎng)時(shí)因網(wǎng)格交替所引起的不穩(wěn)定性。由于計(jì)算中采用了顯式差分格式，所選用的時(shí)間步長(zhǎng)很小，這在很大程度上限制了Z hang Z H. O n m ode li ng i ce dynam i cs of s em i -e nc l o s ed s ea

13、s ona ll y i ce-c ove r ed s eas R . Repo r t S e r i es i n G eophy s i cs ， H e l s i nki ，2000.第6期季順迎等：海冰動(dòng)力學(xué)數(shù)值方法研究進(jìn)展 659它的實(shí)際應(yīng)用。另外，H unke 等22通過(guò)2種途徑來(lái)提高海冰動(dòng)力學(xué)模擬的效率：一是在粘塑性模型的基礎(chǔ)上發(fā)展了彈粘塑性本構(gòu)模型，從而將海冰動(dòng)力學(xué)方程由原來(lái)的拋物線型轉(zhuǎn)變?yōu)殡p曲線方程，進(jìn)而采用顯式差分格式進(jìn)行并行計(jì)算；二是采用共軛梯度法來(lái)并行求解海冰動(dòng)力學(xué)方程。以上兩種處理方法在保證求解精度的條件下，計(jì)算效率有很大提高。在Hi b l e r 1的超松馳

14、迭代運(yùn)算中，對(duì)每個(gè)網(wǎng)格點(diǎn)冰速進(jìn)行逐個(gè)迭代求解，因此可稱(chēng)作點(diǎn)超松馳迭代（P S O R ）。計(jì)算中繁多的迭代過(guò)程在一定程度上限制了海冰動(dòng)力學(xué)的計(jì)算效率。近年來(lái)，Z h a ng 等23，24發(fā)展的線超松馳迭代法（LS O R ）和交替方向隱式法（A D I ）使海冰動(dòng)力學(xué)的模擬效率有了很大提高。在L S O R 方法中，每個(gè)時(shí)間步可引入幾個(gè)或幾十個(gè)擬時(shí)間步，在x 和y 方向上可分別通過(guò)一行或一列對(duì)冰速進(jìn)行迭代求解，因此計(jì)算效率明顯提高。當(dāng)網(wǎng)格尺度分別選用80km 和40 k m ，對(duì)應(yīng)的時(shí)間步長(zhǎng)分別為和24h 和12 h 時(shí)，分別對(duì)北極海冰進(jìn)行了30天和10天的數(shù)值模擬。計(jì)算結(jié)果表明：相對(duì)于PS

15、 O R 方法，L S O R 法的計(jì)算效率提高了3倍多，且能夠更合理地模擬海冰的塑性流變特性22。在海冰動(dòng)力學(xué)模擬的A D I 法中，由于采用x 和y 方向的隱顯交替差分格式，這樣動(dòng)力方程在上、下半時(shí)間步長(zhǎng)內(nèi)均被離散成三對(duì)角方程。它可通過(guò)Th o m a s 運(yùn)算直接求解冰速，而不需要迭代，進(jìn)而更加明顯地提高了計(jì)算效率24。Zh a ng 等23選用不同的網(wǎng)格尺度（160 k m ，40 k m ）、時(shí)間步長(zhǎng)（24h 、12 h 、6h 、3h 、1.5 h 、0.75 h ）和擬時(shí)間步（0、5、15），分別采用A D I 和 L S O R 方法對(duì)北極海冰動(dòng)力過(guò)程進(jìn)行了數(shù)值實(shí)驗(yàn)。計(jì)算結(jié)果均

16、表明：在不采用擬時(shí)間步和采用15個(gè)擬時(shí)間步時(shí)，A D I 方法的計(jì)算效率分別是L S O R 方法的18和8倍，并且計(jì)算精度也有明顯提高。另外，A D I 方法不但可以精確地模擬海冰漂移過(guò)程，且在采用各向異性海冰動(dòng)力學(xué)本構(gòu)模型時(shí)，也能夠大體模擬水道、冰脊和重疊冰的形成24，25。由此可見(jiàn)，歐拉坐標(biāo)系下的有限差分法在海冰動(dòng)力學(xué)模擬中有很高的計(jì)算效率，也能夠較理想地模擬海冰的動(dòng)力學(xué)特征。它是目前應(yīng)用范圍最廣，發(fā)展也相對(duì)最為成熟的海冰數(shù)值方法。有限差分法不但在大尺度長(zhǎng)期海冰數(shù)值模擬中占有主導(dǎo)地位，同時(shí)在中尺度海冰數(shù)值模擬和預(yù)測(cè)中也取得了良好的應(yīng)用效果。但還應(yīng)看到：由于該方法將海冰視為連續(xù)介質(zhì)，描述的

17、是網(wǎng)格內(nèi)海冰的平均狀態(tài)，因此不能精確地模擬海冰動(dòng)力破壞時(shí)的區(qū)域性特點(diǎn)；在計(jì)算海冰動(dòng)力方程中的平流項(xiàng)時(shí)存在數(shù)值擴(kuò)散現(xiàn)象，尤其是處理冰邊緣線時(shí)更有很大的誤差4，8。因此，人們?cè)趹?yīng)用有限差分法進(jìn)行不同尺度下海冰動(dòng)力學(xué)數(shù)值模擬的同時(shí)，一直在發(fā)展能夠精確模擬海冰動(dòng)力學(xué)特征的拉格朗日方法，或拉格朗日與歐拉坐標(biāo)相耦合的方法4，6，7。2海冰動(dòng)力學(xué)的光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)方法光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)（SPH ）方法是 G i n go l d 和M o n ag h a n （1977 ）等在研究天體物理流體動(dòng)力學(xué)中建立和發(fā)展起來(lái)的。海冰動(dòng)力學(xué)SPH 模擬的基本思想是將模擬的海冰區(qū)域分解成足夠多個(gè)海冰質(zhì)點(diǎn)2），每個(gè)質(zhì)點(diǎn)

18、都具有各自的質(zhì)量、動(dòng)量和能量。從圖2中可以看出，海冰質(zhì)點(diǎn)能夠較精確地表示出冰邊緣線位置。近年來(lái)，S PH 方法在河冰和海冰的數(shù)值模擬和預(yù)測(cè)方面得到了廣泛的應(yīng)用5，7，8，26。計(jì)算域內(nèi)任一位置處的海冰參量（冰厚、密集度、冰速以及冰內(nèi)應(yīng)力等）均可通過(guò)一個(gè)選定的權(quán)函數(shù)（或稱(chēng)為核函數(shù)）W 來(lái)由其鄰近海冰質(zhì)點(diǎn)的相應(yīng)參量積分確定。海冰參量的導(dǎo)數(shù)則可通過(guò)分布積分，由權(quán)函數(shù)的導(dǎo)數(shù)計(jì)算。權(quán)函數(shù)一般采用二維G a u s s 函數(shù)，因?yàn)槠渚哂泻芎玫倪B續(xù)性和插值精度，其空間分布如圖3所示。各海冰質(zhì)點(diǎn)的速度可由海冰動(dòng)量方程確定。由于海冰動(dòng)量方程是在拉格朗日坐標(biāo)下求解的，因此不存在F D M 模擬中的平流項(xiàng)，從而消除

19、了數(shù)值擴(kuò)散。計(jì)算中可考慮風(fēng)和流的拖曳力、海面傾斜力以及冰內(nèi)應(yīng)力對(duì)海冰質(zhì)點(diǎn)的作用，其中冰內(nèi)應(yīng)力可根據(jù)Hi b l e r 的粘塑性本構(gòu)模型或 M o hr - C o u l o m b 摩擦準(zhǔn)則進(jìn)行計(jì)算1，7。由于各離散質(zhì)點(diǎn)的疏密不同，且其位置也在隨時(shí)間不斷變化，因此，在計(jì)算過(guò)程中應(yīng)不斷調(diào)整各質(zhì)點(diǎn)的光滑長(zhǎng)度，進(jìn)而得到較精確的計(jì)算值。另外，在搜索海冰質(zhì)點(diǎn)的鄰近質(zhì)點(diǎn)時(shí)，一般搜索區(qū)域越大，計(jì)算結(jié)果越精確。但若搜索區(qū)域太大，則將花費(fèi)過(guò)多的計(jì)算時(shí)間，而當(dāng)搜索區(qū)域太小時(shí)，將降低計(jì)算結(jié)果的精度。一般當(dāng)搜索區(qū)域?yàn)?h 0×4h 0時(shí)（這里h 0為權(quán)函數(shù)W 中的光滑長(zhǎng)度），計(jì)算精度便可滿(mǎn)足要求。另外，

20、在進(jìn)行質(zhì)點(diǎn)搜索時(shí)，不同的搜索方式所占用的搜索時(shí)間是不同的?；谝陨蟂 PH 方法，S h e n 等（1993）發(fā)展了一種拉格朗日坐標(biāo)下的離散塊體模式（D P M ），并成功（圖 圖2光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)中的海冰質(zhì)點(diǎn)分布示意圖F i g.2 D i s t r i bu t i on o f se a i c e p a r t i c l e s i n S P H圖3二維高斯分布示意圖F i g.3 S k e t c h o f G a u ss i an d i s t r i bu t i on i n 2 D地用于河冰、湖冰和海冰動(dòng)力學(xué)特性的模擬5，7，8。隨后，G u tf r a

21、 i nd 等27采用SP H 方法對(duì)冰邊緣線的動(dòng)力作用進(jìn)行了一系列的數(shù)值試驗(yàn)。最近，S P H 方法又被應(yīng)用于渤海海冰動(dòng)力學(xué)的數(shù)值模擬中，其中初始海冰質(zhì)點(diǎn)尺寸為2k m ×2 k m ，計(jì)算的時(shí)間步長(zhǎng)為40s8。將 S P H 與F D 方法的計(jì)算結(jié)果對(duì)比發(fā)現(xiàn)：S P H 方法能夠更合理地模擬出海冰在風(fēng)和潮流作用下的輻合和輻散過(guò)程，且能夠精確地模擬冰邊緣線的位置。S P H 方法克服了有限差分法計(jì)算平流項(xiàng)時(shí)的擴(kuò)散現(xiàn)象，能夠精確地處理冰邊緣線的移動(dòng)位置，在海冰動(dòng)力學(xué)模擬中有較高的計(jì)算精度。但是，由于S P H 在進(jìn)行冰間相互作用的計(jì)算時(shí)，需要大量的循環(huán)計(jì)算，進(jìn)而限制了其在海冰數(shù)值預(yù)測(cè)

22、和長(zhǎng)期數(shù)值模擬中的適用性。因此，提高SP H 模擬的計(jì)算效率，特別是提高冰內(nèi)應(yīng)力計(jì)算效率，將會(huì)進(jìn)一步推進(jìn)S P H 方法在海冰動(dòng)力學(xué)中的應(yīng)用。3海冰動(dòng)力學(xué)的質(zhì)點(diǎn)網(wǎng)格法拉格朗日坐標(biāo)系適合于擾動(dòng)不是很大的多介質(zhì)流體運(yùn)動(dòng)，歐拉坐標(biāo)系則適用于大擾動(dòng)單介質(zhì)流體運(yùn)動(dòng)28。但是在海冰數(shù)值模擬中，一般擾動(dòng)較大，同時(shí)又是冰水兩相介質(zhì)共存狀態(tài)，所以將拉格朗日與歐拉坐標(biāo)系耦合是比較理想的方法。網(wǎng)格質(zhì)點(diǎn)方法（P I C ）就是一種很有效的耦合方法。在 P I C 方法中，流體具有兩重性，即一方面把流體看成連續(xù)介質(zhì)，從而在沒(méi)有物質(zhì)輸運(yùn)的情況下計(jì)算流場(chǎng)變化；另一方面又把流體看成若干個(gè)帶有一定質(zhì)量的質(zhì)點(diǎn)，然后在固定的歐拉矩

23、形網(wǎng)格上研究這些質(zhì)點(diǎn)的運(yùn)動(dòng)，以及質(zhì)量、動(dòng)量和能量的輸運(yùn)。因此，P I C 方法具有一般歐拉方法的優(yōu)點(diǎn)，能夠計(jì)算扭曲比較嚴(yán)重的二維流體力學(xué)模型；同時(shí)由于引入了相當(dāng)于拉格朗日質(zhì)量團(tuán)的質(zhì)點(diǎn)，具有計(jì)算多種物質(zhì)和處理自由面運(yùn)動(dòng)的能力。Fl a t o 4最早將PI C 方法引入海冰動(dòng)力學(xué)數(shù)值模擬中。該方法可以處理海冰動(dòng)力學(xué)中冰水介質(zhì)的兩相流問(wèn)題，同時(shí)又能夠精確模擬冰邊緣線的動(dòng)力過(guò)程。在海冰動(dòng)力學(xué)的P I C 模擬中，首先將歐拉網(wǎng)格內(nèi)的海冰離散成一個(gè)個(gè)獨(dú)立的質(zhì)點(diǎn)，然后將這些質(zhì)點(diǎn)插值到固定歐拉網(wǎng)格內(nèi)進(jìn)行動(dòng)量計(jì)算（圖4）。通過(guò)求解動(dòng)量方程得到固定網(wǎng)格的速度場(chǎng)，再將該速度場(chǎng)插值到各質(zhì)點(diǎn)所在位置得到各質(zhì)點(diǎn)的速度，

24、并根據(jù)網(wǎng)格內(nèi)各海冰質(zhì)點(diǎn)的重新分布確定各網(wǎng)格的平均冰厚和密集度。利用P I C 方法，F(xiàn) l a t o 4首先進(jìn)行了渦旋風(fēng)場(chǎng)作用下矩形區(qū)域內(nèi)海冰的動(dòng)力演化過(guò)程的數(shù)值試驗(yàn)，計(jì)算中采用了Hi b l e r 的粘塑性海冰本構(gòu)模型。在該數(shù)值試驗(yàn)中，計(jì)算區(qū)域?yàn)?00km 的正方形海域，海冰均勻覆蓋上半部計(jì)算域，初始冰厚為1.0 m ，時(shí)間步長(zhǎng)為2.4h。當(dāng)海冰質(zhì)點(diǎn)按5 k m 間隔初始均勻分布，有限差分的網(wǎng)格尺寸分別取為10km 和20 k m 時(shí)，這樣每個(gè)網(wǎng)格內(nèi)分別有4和16個(gè)海冰質(zhì)點(diǎn)。由此模擬的5天后的海冰厚度和密集度均非常理想，且當(dāng)取用20k m 網(wǎng)格尺寸時(shí)的冰厚等直線分布更為光滑。當(dāng)海冰質(zhì)點(diǎn)和

25、網(wǎng)格尺寸均取為10 k m 時(shí)，得不到光滑的冰厚等直線。由此得出：當(dāng)每個(gè)網(wǎng)格內(nèi)有10左右海冰質(zhì)點(diǎn)時(shí)，可得到較光滑的海冰厚度等直線，而網(wǎng)格內(nèi)的質(zhì)點(diǎn)數(shù)目對(duì)冰邊緣線的影響不大。隨后，F(xiàn) l a t o 4對(duì) B ea u f o r t 海的海冰動(dòng)力學(xué)過(guò)程進(jìn)行了9天的PI C 數(shù)值模擬。盡管計(jì)算中 沒(méi)有考慮海冰熱力作用過(guò)程，計(jì)算得到的冰邊緣線位置與實(shí)際結(jié)果有一定的差異，但其精度已明顯高于有限差分的計(jì)算結(jié)果。圖4海冰動(dòng)力學(xué)質(zhì)點(diǎn)網(wǎng)格法中的質(zhì)點(diǎn)和網(wǎng)格布局示意圖F i g.4 S c h e m e of p a r t i c l e a n d c e l l i n P I Cme t h od f

26、or s e a i c e d y n am c i s 為進(jìn)一步提高P I C 方法在海冰動(dòng)力學(xué)中的計(jì)算精度，H u a ng 等29采用了Zh a ng 和 H i b l e r （1997 ）建立的考慮擬時(shí)間步的LS O R 差分方法進(jìn)行海冰動(dòng)量方程的計(jì)算，并在每個(gè)網(wǎng)格內(nèi)放置更多的海冰質(zhì)點(diǎn)100個(gè)）。通過(guò)對(duì)冰邊緣區(qū)動(dòng)力學(xué)的F D 和P I C 方法的數(shù)值對(duì)比發(fā)現(xiàn)，這兩種方法的計(jì)算結(jié)果均能達(dá)到理想的結(jié)果，但PI C 方法由于避免了F D 中的平流擴(kuò)散影響，計(jì)算的冰邊緣線位置更為精確。近年來(lái)，在渤海海冰數(shù)值模擬中也在嘗試P I C 方法，以提高海冰動(dòng)力學(xué)的計(jì)算精度，滿(mǎn)足冰期油氣開(kāi)發(fā)的工程

27、要求。雖然P I C 方法在歐拉網(wǎng)格和拉格朗日質(zhì)點(diǎn)間來(lái)回插值的過(guò)程中仍受數(shù)值發(fā)散影響，但它具有一般歐拉方法的優(yōu)點(diǎn)，能夠計(jì)算扭曲比較嚴(yán)重的二維流體力學(xué)模型；同時(shí)由于引入了拉格朗日質(zhì)點(diǎn)，具有計(jì)算多種物質(zhì)和處理自由面運(yùn)動(dòng)的能力；另外，P I C 法的另一個(gè)重要優(yōu)點(diǎn)是對(duì)于重點(diǎn)分析區(qū)域可采用較多的質(zhì)點(diǎn)進(jìn)行精確計(jì)算，而在其它區(qū)域只保持較少的質(zhì)點(diǎn)以減小計(jì)算量。目前，P I C 方法越來(lái)越受到海冰數(shù)值模擬工作者的重視。在海冰動(dòng)力學(xué)的歐拉和拉格朗日耦合方法研究中，除了PI C 方法外，P r it c h a r d 等6建立了一種自適應(yīng)網(wǎng)格方法。該方法采用拉格朗日網(wǎng)格，網(wǎng)格隨冰邊緣線的運(yùn)動(dòng)而自動(dòng)調(diào)整形狀和位置

28、。盡管這種方法采用了運(yùn)動(dòng)的網(wǎng)格，為跟蹤冰邊緣線移動(dòng)提供了一個(gè)有效的途徑，但數(shù)值發(fā)散依然存在，而且這種運(yùn)動(dòng)網(wǎng)格法只適用于冰緣運(yùn)動(dòng)較光滑平緩的區(qū)域。當(dāng)冰邊緣線運(yùn)動(dòng)急劇和不規(guī)則時(shí)，運(yùn)動(dòng)網(wǎng)格將會(huì)嚴(yán)重扭曲，甚至可能導(dǎo)致網(wǎng)格點(diǎn)的重疊而使模擬失敗。另外，R h ee m 等9借鑒任意拉格朗日歐拉法（A L E ）的基本思想，并考慮海冰在中、小尺度下分別表現(xiàn)出的離散分布和流變連續(xù)特性，發(fā)展了一種計(jì)算中尺度海冰動(dòng)力學(xué)的分布質(zhì)量/離散塊體模型（DM D F ）。該方法采用拉格朗日法模擬離散性的大面積海冰，并對(duì)海冰的局部流變過(guò)程進(jìn)行歐拉有限差分計(jì)算。該模型已應(yīng)用于Okh o t s k 海的海冰動(dòng)力過(guò)程模擬，計(jì)算結(jié)

29、果與衛(wèi)星遙感資料相吻合。該模型能夠較好地反應(yīng)海冰在自然條件下的離散和連續(xù)狀態(tài)，有很強(qiáng)的靈活性，且計(jì)算方便。4海冰動(dòng)力學(xué)的顆粒流模型為研究海冰在波浪作用下的碰撞，S h e n 等10，11將顆粒流體動(dòng)力學(xué)理論應(yīng)用到海冰動(dòng)力學(xué)中，并建立了二維冰塊碰撞流變學(xué)模型。在最初的二維海冰碰撞流變學(xué)中，海冰被離散為大小相同的圓盤(pán)，各冰盤(pán)在運(yùn)動(dòng)中相互碰撞、摩擦，其相互作用采用H i b l e r 的粘塑性本構(gòu)模型，其中的粘性系數(shù)是冰塊大小、厚度和回彈系數(shù)的函數(shù)11。L u 等30和L e pp a r a n t a 等31采用大小不同的圓盤(pán)對(duì)海冰與波浪的相互作用以及碎冰區(qū)的動(dòng)力特征做了進(jìn)一步的研究?？紤]不

30、同尺寸和形狀的海冰動(dòng)力學(xué)顆粒流模型如圖5所示。圖5不同尺寸和形狀海冰的顆粒流模型示意圖F i g.5 S k e t c h o f se a i c e g r a nu l ar f l ow d y n am i c s w i t hd i ff ere n t s h a p e a n d s i z e p a r t i c l e s海冰重疊和堆積是海冰動(dòng)力學(xué)中的一個(gè)重要現(xiàn)象。由于它們形成過(guò)程復(fù)雜，且對(duì)冰區(qū)結(jié)構(gòu)物有很大威脅，因此人們一直致力于它們的力學(xué)成因、物理（968 地球科學(xué)進(jìn)展第 19 卷 力學(xué)性質(zhì)的現(xiàn)場(chǎng)監(jiān)測(cè)及其形成過(guò)程的數(shù)值模擬工 33 32， 34 Hopkins

31、等 。最近， 建立了三維海冰碰撞 作 流變學(xué)， 它不但可用于研究小尺度下海冰與波浪的 相互作用， 同時(shí)還應(yīng)用于中尺度下冰脊和重疊冰形 成過(guò)程的數(shù)值模擬 34 顆粒流模型能夠合理地模擬海冰動(dòng)力學(xué)過(guò)程中的斷 裂、 重疊和堆積現(xiàn)象， 具有物理意義明確、 計(jì)算精度高 的優(yōu)點(diǎn)。但顆粒流模型在計(jì)算效率上要明顯低于其 它模型， 目前還不能用于海冰數(shù)值預(yù)報(bào)和長(zhǎng)期海冰動(dòng) 力學(xué)模擬的時(shí)效要求； 另外， 由于它需要輸入精確的 海冰初始信息， 這在現(xiàn)時(shí)海冰監(jiān)測(cè)工作中也有很大困 難。因此， 在保證顆粒流模型計(jì)算精度的前提下， 提 高計(jì)算時(shí)效和簡(jiǎn)化輸入?yún)?shù)是必要的。 。 在海冰動(dòng)力學(xué)的顆粒流模型中， 海冰重疊和堆 積的力

32、學(xué)模型主要包括小變形時(shí)的粘彈性和發(fā)生冰 脊或擠壓破壞時(shí)的塑性?xún)刹糠?。采用以上彈?Hopkins 塑性模型， 并考慮海冰在堆積過(guò)程中的能量守恒， 35 100 km）下進(jìn)行了一系列 在中尺度 10 （ 5 結(jié)語(yǔ)與展望 在對(duì)不同尺度海冰數(shù)值模擬和預(yù)測(cè)中， 針對(duì)海 冰動(dòng)力學(xué)的不同特點(diǎn)和海冰數(shù)值模擬的應(yīng)用要求， 人們建立了一系列的海冰動(dòng)力學(xué)數(shù)值方法。這些數(shù) 拉 值方法主要包括歐拉坐標(biāo)下的有限差分法 FD）、 （ 格朗日坐標(biāo)下的光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)方法 SPH）、 （ 歐拉與拉格朗日坐標(biāo)耦合的質(zhì)點(diǎn)網(wǎng)格方法 PIC）， （ 以及將海冰視離散介質(zhì)的顆粒流方法 GF （ ）。 有限差分法是發(fā)展最早、 相對(duì)比較

33、完善的計(jì)算 方法。它計(jì)算效率高、 穩(wěn)定性好， 已廣泛地應(yīng)用于 大、 中尺度的長(zhǎng)、 中期海冰動(dòng)力學(xué)模擬和預(yù)測(cè)中。但 由于它存在動(dòng)量方程平流項(xiàng)計(jì)算的擴(kuò)散現(xiàn)象和冰邊 緣線確定的不精確等缺點(diǎn)， 為此人們也在發(fā)展新的 計(jì)算方法以提高海冰動(dòng)力學(xué)模擬的精度。近來(lái)年， SPH 方法在河冰和海冰動(dòng)力學(xué)模擬中得到了廣泛的 應(yīng)用。它能夠精確地模擬冰邊緣線的位置， 且克服 了 FD 方法的數(shù)值擴(kuò)散現(xiàn)象。但在海冰參數(shù)， 特別 是冰內(nèi)應(yīng)力的插值迭代過(guò)程中， 大量的繁瑣計(jì)算降 低了它的計(jì)算效率。PIC 方法由于同時(shí)具有歐拉和 拉格朗日方法的優(yōu)點(diǎn)， 能確較精確地確定冰邊緣線 位置， 也極大地降低了網(wǎng)格間的擴(kuò)散現(xiàn)象。雖然它 在

34、網(wǎng)格和質(zhì)點(diǎn)的相互插值過(guò)程中仍存在一定的擴(kuò)散 現(xiàn)象， 但在大中尺度海冰動(dòng)力學(xué)模擬中有很好的應(yīng) 用前景。 以上 3 種計(jì)算方法均是將海冰視為連續(xù)介質(zhì)條 件下進(jìn)行的。但自然條件下， 海冰在大中尺度下表 現(xiàn)出流變特性的同時(shí)， 在風(fēng)、 流等海洋條件作用 浪、 下還不斷地?cái)嗔选?重疊和堆積， 在小尺度下又顯示出 很強(qiáng)的不連續(xù)特性。顆粒流體方法就是針對(duì)這種現(xiàn) 象建立， 有非常明確的物理意義， 并已成功地進(jìn)行了 海冰動(dòng)力過(guò)程的數(shù)值模擬。但由于它將海冰離散為 諸多均勻或不同尺寸形狀的顆粒體， 顆粒間的相互 作用計(jì)算需要取用很小的時(shí)間步長(zhǎng)和大量計(jì)算來(lái)完 成， 這在三維顆粒流模型中表現(xiàn)的更為突出。因此， 目前該方法

35、在計(jì)算時(shí)效上還不能滿(mǎn)足海冰數(shù)值預(yù)測(cè) 海冰堆積數(shù)值模擬實(shí)驗(yàn)， 并對(duì)計(jì)算區(qū)域內(nèi)海冰應(yīng)力 進(jìn)行了統(tǒng)計(jì)分析， 并得到了海冰在離散狀態(tài)下發(fā)生 25 重疊和堆積的屈服曲線。Hibler 發(fā)現(xiàn)該屈服曲線 同其改進(jìn)的粘塑性模型比較接近， 能夠很好地反映 海冰在堆積過(guò)程中的塑性應(yīng)力狀態(tài)。 Hop利用 海 冰 重 疊 和 堆 積 的 顆 粒 流 方 法， 16 kins 對(duì)冰池內(nèi)海冰堆積過(guò)程的物理實(shí)驗(yàn)進(jìn)行了數(shù) 值模擬。該實(shí)驗(yàn)是將一規(guī)則區(qū)域內(nèi)的碎冰塊以一定 的速率從冰池一端向另一端推動(dòng)， 在推進(jìn)過(guò)程中對(duì) 表征冰內(nèi)應(yīng)力的平均推力和厚度分布進(jìn)行了測(cè)量。 Hopkins 的模擬結(jié)果與實(shí)驗(yàn)相一致， 即在開(kāi)始推進(jìn)階 段， 海

36、冰彈性應(yīng)力迅速增加， 當(dāng)海冰發(fā)生堆積后， 應(yīng) 16 力增長(zhǎng)率趨近于零， 海冰進(jìn)入塑性階段 。最近， 海冰數(shù)值模擬的顆粒流方法已越來(lái)越引起人們的興 趣， 并嘗試擴(kuò)展到大尺度海冰 動(dòng)力學(xué)的數(shù)值 模擬 36 38 。 中 39 40 Loset ， 建 類(lèi)似于海冰動(dòng)力學(xué)的顆粒流模型， 立了碎冰區(qū)海冰動(dòng)力學(xué)的離散元模型 DEM ） （ 。該模 型將散狀分布的冰塊視為離散材料， 也采用不同直徑 的圓盤(pán)來(lái)模擬冰塊的相互碰撞。冰塊間碰撞時(shí)的作 Loset 用力采用粘彈塑性本構(gòu)模型計(jì)算。采用該模型， 還進(jìn)行了冰塊對(duì)攔冰柵作用的數(shù)值模擬研究， 并分析 了密集度、 冰速和冰塊大小對(duì)冰力的影響。 海冰動(dòng)力學(xué)的顆粒流

37、模型在計(jì)算方法方面異于 建立在連續(xù)介質(zhì)力學(xué)基礎(chǔ)上的有限差分法、 光滑質(zhì) 點(diǎn)流體動(dòng)力學(xué)方法和質(zhì)點(diǎn)網(wǎng)格方法， 同時(shí)在本構(gòu)模 型方面也不同于傳統(tǒng)的粘塑性和彈塑性本構(gòu)模型。 它將海冰離散成多個(gè)質(zhì)點(diǎn)， 該質(zhì)點(diǎn)由最初的二維均 一圓盤(pán)發(fā)展到考慮冰塊尺寸和形狀的三維模型， 計(jì) 算精度也在不斷提高。顆粒流模型也由最初用于處 理波浪和海冰的相互作用， 發(fā)展到研究海冰重疊和 Hopkins 又將其進(jìn)一步發(fā)展， 堆積過(guò)程。近年來(lái)， 并 試圖應(yīng)用于大中尺度下的海冰動(dòng)力學(xué)的數(shù)值模擬。 大量的數(shù)值試驗(yàn)、 模型實(shí)驗(yàn)和現(xiàn)場(chǎng)監(jiān)測(cè)資料均表明： 第 6 期季順迎等： 海冰動(dòng)力學(xué)數(shù)值方法研究進(jìn)展 969 和中長(zhǎng)期海冰動(dòng)力學(xué)模擬的要求。

38、 通過(guò)對(duì)以上海冰動(dòng)力學(xué)數(shù)值方法的分析和討 論， 并結(jié)合海冰動(dòng)力學(xué)數(shù)值模擬中的工作， 建議在 中、 長(zhǎng)期海冰動(dòng)力學(xué)模擬和預(yù)測(cè)中， 采用連續(xù)介質(zhì)條 PIC 以保證計(jì)算效率， 其中 件下的 FD 、 或 SPH 方法， SPH 方法計(jì)算效率相對(duì)較低， 但其計(jì)算冰邊緣線精 rine Science and Technology， 1997 ， （2 2 ）： 101 121 . 10 Shen H H ，H ibler W D ，Leppara nta M . On applying granular flow theory to a deform in g broken ice field J .

39、Acta Mechanics ， 1 986 ， ： 63 143 160 . 11 Shen H H ，Hibler W D ，Leppar anta M . The role of floe collisions in sea ice r heology J . Journa l of Geophysical R esearch， 198 7， 2 C10 9（ 7 ）：7 085 - 096 . 12 Hopkins M A ，Hibler W D ，F(xiàn)lato G M . On the num erical sim ulation of the sea ice ridging pr o

40、cess J . Journa l of Geophysical Resear ch， 1991 ， （C3）： 809 - 82 0. 96 4 4 13 Hopkins M A ，Tuh kuri J，Lensu M . Rafting and riding of th in ice sheets J . Jou rnal of Geophysical Res ea rch， 1999 ， 104 （C6）：13 605 -13 613 . 14 Lepparanta M ，H ib ler W D . Mesoscale sea ice defor m ation in m argin

41、al ice zone dynam ics J . Jour nal of Geophys ical Re sear ch ， 1987 ， ： 06 0- 070 . 92 7 7 15 Zhang Z H ，Lepparanta M . Mod eling the influence of ice on sea levelvariations in the Baltic Sea J . G eoph ysica， 199 5， （2）： 31 31 - . 45 度要明顯高于其它兩種方法， 且不存在平流項(xiàng)的擴(kuò) GF 散； 模型可用于研究不同尺度， 特別是小尺度下 的海冰重疊動(dòng)力學(xué)特性， 但

42、其計(jì)算時(shí)效還需要進(jìn)一 步改進(jìn)以擴(kuò)展它的適用范圍。 在海冰動(dòng)力學(xué)模擬中， 除了采用現(xiàn)有方法外， 還 應(yīng)針對(duì)海冰的動(dòng)力作用強(qiáng)烈、 冰水移動(dòng)邊界復(fù)雜的 特征， 深入研究適用于中小尺度下海冰動(dòng)力學(xué)模擬 的其它數(shù)值方法， 如任意拉格朗日歐拉 ALE）數(shù)值 （ 方法。ALE 方法可根據(jù)海冰在不同尺度下分布的 連續(xù)和非連續(xù)性特征， 在不同的網(wǎng)格尺寸下靈活地 采用歐拉坐標(biāo)系和拉格朗日坐標(biāo)系進(jìn)行海冰動(dòng)力學(xué) 精確描述。 海冰動(dòng)力學(xué)數(shù)值方法是與海冰本構(gòu)模型， 甚至 海冰熱力過(guò)程等密切相關(guān)的。因此， 在開(kāi)展海冰動(dòng) 力學(xué)數(shù)值方法研究時(shí)， 除考慮不同時(shí)空尺度下海冰 的動(dòng)力特性、 海冰模擬的精度要求等因素外， 還應(yīng)與 其它

43、海冰研究?jī)?nèi)容相結(jié)合， 進(jìn)而全面地提高海冰模 擬和預(yù)測(cè)的精度。 ） 參考文獻(xiàn)（References ： 1 Hibler W D . A dynam ic ther m odynam ic sea ice m ode l J . Jour na l of G eoph ysical Ocean ogr aphy， 197 9，9： 17- 6. 8 84 2Leppar anta M ，H ib ler W D . The r ole of plastic ice interaction in Margin al Ice Zone dynam ics J . Jour nalof Geophys

44、ical R es earch， 1985 ， （C6 90 1 ）： 1 899 -11 909 . 3 W u H uiding （吳輝碇 ），Bai Shan （白珊），Zhang Zhanh ai （張占 海 . Num erical sim ulation for dynam icalpr ocesses ofsea ice J . ） Acta Oceanologica Sinica 1988 ， （3）： - （in 16 303 325 （海洋學(xué)報(bào)）， Chinese ） 4 Flato G M . A particle-in -cell sea-ic e m od el J

45、. Atm osphere and Oceanography，993 ， （3）：339 - 58. 1 31 3 C 5 Shen H T ， hen Y C ，W ake A ，et al. Lagrangian discrete parcel sim u lation of river ice dynam ics J . Inter national Jour nal of Off s hor e an d P olar E ngineering， 93 ，（4 19 3 ）： 832 332 . 6 Pritchard R ，Muelle r A ，Hanzlick D ，et al.

46、 Fore castin g Bering Sea ice edge behavior J . Journalo f Geophys ical Research ， 0， 199 95：775 - 88. 7 7Sh en H T ，Su Ju nshan ，Liu Lian wu. S PH sim ulation of river ice dynam ics J . Jour nal of Com putational Physics 2000 ，16 5： ， 752 770 . Sh 8 W ang Z ， en H T ，W u H . A Lagr angian sea ice m

47、 odelw ith discr ete parcel m ethod A. I n：Shen ， in Su rface W ater C . ed Ice USA ： Rotterdam ， 1998 .313 - 20 . 3 16 Zhang Zha nhai （張占 海），B ai Sh an （白珊），W u H uid in g （吳輝 et 碇）， a l. Num er ical forecasting system for sea ice op eration in the Bohai Sea J . Mar in e For eca sts （海 洋 預(yù) 報(bào)），1994

48、，11 11 （2）： -18（in Chinese ） 17 W u Hu idin g （吳 輝 碇），Mathem atical representations of sea ice therm odynam ic processes J . Oceanologia et lim nodynam ic logia sinca 海 洋 與 湖 沼），1991 ，2 2 4）： - 27（in Chi 321 3 （ （ nese . ） 18 Bai Shan（白 珊），W u Huiding（吳 輝 碇）. Num erical sea ice forecast for the Boha

49、i Se a J . Acta Meteor olog ica Sinica （氣 象學(xué) 1998 ， （2 13 9- 1 in C hinese . 56 ）： 15 （ 報(bào)）， ） 19 W ang Renshu （王仁 樹(shù)）. Num ericaltests ofsea ice fo r the B ohai Sea J . Acta Oceanologica Sinica （海 洋 學(xué) 報(bào)），1984 ，6（4）： 572 - （in Ch in ese . 580 ） Liu 20 W ang Renshu （王仁 樹(shù)）， Qinzhen g （劉欽 政），Chen W eib i

50、n （陳 偉斌），et al. Num erical sim u lation an d tests of sea ice drift process in the Bohai Sea J . Ocean olog ia et lim nologia Sinica 199 （海 洋與 湖沼）， 4，25 3）：40 8-415（in C hin ese . （ ） 21Ip C F ，Hibler W D ，F(xiàn)lato G M . On the effectofthe rheology on seasonal sea ice sim ulations J . Annu al Glaciolog

51、y，1991 ，15： 17 - . 25 22 Hunk e E C ，Dukowicz J K . An elastic-viscous-plastic m odel for sea ice dynam ics J . Journal of Phys ical Oceanography， 97， 19 27 ： 849 - 86 7. 1 1 23 Zh ang J，H ib ler W D . On an efficientnum ericalm ethod for m odeling sea ice dynam ics J . Journal of Geophysical R esea

52、rch， 199 7， 02 C4 8 691 - 702 . 1 （ ）： 8 24 Zh ang J，Roth rock D . Mod eling Arctic sea ice w ith an efficient plastic solution J . Jour nal of Geophys ical R es earch， 000 ， 2 105 3 （C2）： 3 25-3 338 . 25 Hibler W D . Sea ice fractur ing on the lar ge scale J . E ngineer ing Fracture Mechanics 2001

53、， ： 013 - 043 . 68 2 2 ， Yue Qianjin 26 Ji Shun ying （季 順迎）， （岳 前進(jìn)）. Num erical sim ulation of localdr ifting sea ice in Liaodong Bay by Sm ooth ed Parti cle H ydr odynam ics Method J . Hydro-science a nd Engineering 2001 ， 4）：8- （in Ch inese . 15 （水 利水 運(yùn)工 程學(xué) 報(bào)）， （ ） 27 Gutfr aind R ，Savage S B . Ma

54、rginal ice zone rheology：Com parison of resu lts from con tinuu m -plastic m odels and discr eteparticel sim ulation s J . Jour nal of G eophys ical Res ea rch，1 997 ，102 12 12 （C6）： 647 - 661 . 9 Rheem C K ，Yam aguchi H ， Kato H . Distr ibuted m ass /discrete floe m odel for pack ice r heology com

55、pu tation J . Journal of Ma 28 Li Deyuan （李 德元），Xu Guorong （徐 國(guó)榮），et al. Num erical 970 地球科學(xué)進(jìn)展第 19 卷 Method for 2 D Un steady Fluid Mec han ics M . Bei jing：Scien ce Pr ess， 199 8.45 - （in Chinese 68 ） ead ice an d floes 35 Hopkins M A . On th e m e sosacle inter action of l 19 J . Jour nalof Geophy

56、s ica l Research ， 96，101（C8）：18 31 5 18 326 . 29 Huang Z J，Savage SB . Par ticle-in-cell and finite d ifferen ce app roach es for the stu dy of m argin al ice zone problem s J . Cold Regions S cience and Techn olog y， 998 ， ： - . 1 28 1 28 Tryde P . On the r ole of ice interaction due to 30 Lu Q M ，Larsen J， floe collisions in m ar ginalice zone dynam ics J . Jou rnalof Geo physical Research ， 1989 ， ： 5 25- 537 . 94 14 14 31 Lepparan ta M ，Lensu M ，Lu Q M . Shear flow of sea ice in the