海洋數(shù)值模型的理論及應(yīng)用.ppt

第一節(jié)海洋數(shù)值模型的發(fā)展概況 觀測 試驗理論分析數(shù)值模型 海洋的研究手段 海洋數(shù)值模型發(fā)展的客觀要求 海水運動方程推導(dǎo)過程中的動力學(xué)假設(shè) 比如海水靜力學(xué)假設(shè) 都是建立在大部分空間尺度大于10km 時間尺度大于慣性周期 2 f 的情況下的海水運動方程不存在精確解 1 偏微分方程本身是非線性的 2 將一個實際問題具體化時所需要的環(huán)境場如海底地形 海岸線幾何形狀 表面驅(qū)動力等都無法給出解析函數(shù)表達式 只能在離散化的空間和時間區(qū)域取得因此 求解海水運動方程只能通過近似求解的辦法來確定 而這些數(shù)值解又可能引入明顯的近似誤差 海洋數(shù)值模型發(fā)展的歷史 海洋模型到按其水平網(wǎng)格的離散方式以及所使用的垂向坐標系的不同大致經(jīng)歷了如下幾個發(fā)展階段最早出現(xiàn)并且還在使用的海洋模型是Bryan等人開發(fā)的基于原始方程的低階精度的有限差分模型 它在水深方向采用z坐標系目前常見的HOPS HarvardOceanPredictionSystem MOM GFDLModularOceanModel POP ParallelOceanProgram NCOM NCARCommunityOceanModel 模型在某種程度上都可以認為是該模型改進版 海洋數(shù)值模型發(fā)展的歷史 上個世紀70年代 sigma坐標系開始應(yīng)用于海洋模型在水深方向 比如目前被廣泛使用的POM PrincetonOceanModel ECOM EstuarineCoastalandOceanModel ROMS RegionalOceanModelingSystem 模型都屬于這種類型的模型與傳統(tǒng)的z坐標系相比 sigma坐標可以更好地貼合海底地形變化 海洋數(shù)值模型發(fā)展的歷史 自上世紀90年代以來 非結(jié)構(gòu)化網(wǎng)格技術(shù)開始在海洋模型中得到應(yīng)用 相應(yīng)的也出現(xiàn)了基于非結(jié)構(gòu)化網(wǎng)格離散方式的有限元模型 如SEOM模型 SpectralFiniteElementOceanModel 和有限體積模型 如FVCOM雖然與傳統(tǒng)的結(jié)構(gòu)化網(wǎng)格相比 非結(jié)構(gòu)化網(wǎng)格可以更好地擬合陸地邊界 但是代碼實現(xiàn)上的困難以及計算穩(wěn)定性的問題使其迄今還沒有得到非常廣泛的應(yīng)用 新一代的海洋數(shù)值模型 新一代的海洋模型廣泛采用隨地坐標系 terrain followingcoordinates 進而促進了有關(guān)時間步長 對流項和壓力梯度項等數(shù)值算法的改進進入新世紀以來 下一代的海洋數(shù)值動力模型正在緊鑼密鼓的研制中 代表性的是TOMS Terrain followingOceanModelingSystem 它融合了目前最先進的物理知識 數(shù)值方法和數(shù)據(jù)同化技術(shù)http www aos princeton edu WWWPUBLIC htdocs pom TOMS htm 海洋模型的商業(yè)軟件 就海流的仿真建模研究而言 現(xiàn)階段國外已形成了不少具有代表性的 具有強大前后處理功能及核心仿真建模技術(shù)的軟件或系統(tǒng) 商業(yè)軟件 如加拿大DalhousieUniversity研究的DALCOAST河口海岸預(yù)報系統(tǒng) 丹麥的DHIWater Environment機構(gòu)開發(fā)的的MIKE系列軟件系統(tǒng) 荷蘭的WLDelftHydraulics開發(fā)的的Delft 3D軟件 第二節(jié)基于POM模型的一般介紹 水動力模型基本特點 垂直方向的 坐標變換水平網(wǎng)格采用正交曲線坐標和ArakawaC差分格式水平時間差分采用顯格式 而垂直差分為隱格式自由表面可以模擬水位變化垂直和水平方向的混合擴散分別采用2 5階的Mellor Yamada湍流閉合模式和Smagorinski模式內(nèi)外模態(tài)分別處理速度較慢的內(nèi)重力波和速度較快的外重力波以提高整個模式計算效率包含了海水的熱動力過程 Sigma坐標變換 模型采用 a 垂向sigma坐標系和 b 水平ArakawaC網(wǎng)格 1 控制方程 連續(xù)方程 海水運動方程 2 3 4 溫度方程 湍流動能方程 鹽度方程 湍流混合長度方程 5 6 7 8 其中 是基于sigma坐標系的垂向流速矢量 它與笛卡兒坐標系下的垂向流速W之間的關(guān)系可表示為 另外 方程中分別表示海水的實際密度 Boussinesq近似密度以及密度擾動 其中 是實際水深 9 M Y湍流閉合模型 profq 在湍流動能及混合長度方程中分別表示湍流混合系數(shù) 熱擴散系數(shù)和湍流動能的垂直擴散系數(shù) 墻近似函數(shù) 其中 0 4是vonKarman常數(shù)湍流動能和混合長度方程組由下列等式關(guān)系閉合 其中 是穩(wěn)定函數(shù) 10 在2 5階M Y湍流閉合模型中 可表示為 其中 在不穩(wěn)定層次條件下 的上限值為0 023 對應(yīng)的值分別為2 0145和2 4401 在穩(wěn)定層次條件下 它下限值為 0 28 對應(yīng)的值分別為0 0470和0 0461 是2 5階M Y湍流閉合模式參數(shù) 由實驗測得 11 M Y湍流模型的適用條件 M Y湍流模型已被廣泛應(yīng)用于淺海潮汐 風(fēng)生表面混合層的模擬以及底邊界層的研究 該模型在混合較弱的層化流體中對湍流混合系數(shù)計算效果不佳 應(yīng)用于河口的可行性也有待于探討 雖然該模型無論是在物理上還是在數(shù)學(xué)上都存在著明顯的缺陷 但是目前依然得到廣泛應(yīng)用 只有深入了解海水的微細結(jié)構(gòu)以及海洋湍流結(jié)構(gòu)和性質(zhì)后 才有可能構(gòu)建出更為完善的湍流閉合模型 海水運動方程的水平擴散項 其中水平切應(yīng)力是水平渦粘性系數(shù) 由Smagorinsky模式計算 12 13 Smagorinsky水平擴散模式 其中參數(shù)C為無量綱值 其取值為0 1或0 2 如果計算網(wǎng)格劃分得足夠細的話 C可以取值為0由Smagorinsky模式計算的水平渦粘性系數(shù)隨著網(wǎng)格精度的改善和速度梯度的減小而減小 14 溫度鹽度的水平擴散項 其中是水平熱擴散系數(shù) 一般地 是一個小值 取0 1或0 2 甚至在某些情況下可以取為0 15 16 垂向邊界條件 1 海面 17 2 海底 其中 分別是湍流閉合模型參數(shù)和摩擦速度 分別是海面和海底的摩擦系數(shù) 18 模型的外模態(tài) 控制近岸環(huán)流的海水運動方程包含了運動速度較快的外重力波和速度較慢的內(nèi)重力波 因此進行模式分裂可以大大提高模型的計算效率模型的外模態(tài)就是計算自由水位和垂向平均的流速變化 需要較小的時間步長內(nèi)模態(tài)主要模擬流速矢量 溫度 鹽度等三維結(jié)構(gòu)變化 對時間步長的要求相對寬松采用模式分裂技巧可以有效地減少因外模態(tài)計算所消耗的計算時間 從而提高計算效率 垂向積分的海水運動方程 其中 連續(xù)方程 海水運動方程 19 20 21 水平擴散項 水平數(shù)值耗散項 22 23 模型代碼中的一些符號 im jm kb imm1 jmm1 kbm1 mode isplitdte dti days umolz zz dz dzzaam2d art aru arv dum dvm fsmdx dy h el d ua va ut vt et corswradwusurf wvsurf wubot wvbot wtsurf wssurf radu v t s rho l q2km kh kq aam aahrmean tclim sclim 數(shù)值積分流程圖 外模態(tài)和內(nèi)模態(tài)的相互作用 內(nèi)模態(tài)的計算方法 三維變量的計算 以T為例 可以分解為垂直擴散項的計算和水平對流擴散項的計算 以溫度方程為例方程寫為 求解過程分為兩步 第一步計算水平對流擴散項 采用中心差分格式 由advt子程序?qū)崿F(xiàn) 24 25 水平對流擴散項 垂直擴散項 第二步計算垂向擴散項 數(shù)值方法采用 蛙跳 格式 由proft子程序?qū)崿F(xiàn) 由于 蛙跳 格式在奇數(shù)時間步時的解與在偶數(shù)時間步時的解會發(fā)生偏離 因此每一時間步的計算結(jié)束后還需對上述解進行平滑處理 即 為一個小值 一般取0 05 處理后 26 27 計算網(wǎng)格的安排 二維外模態(tài)網(wǎng)格分布 三維內(nèi)模態(tài)網(wǎng)格分布 水平對流項的計算 雖然模型采用有限差分計算格式 但是對于每個網(wǎng)格對流項的計算都是按照有限容積的方法進行處理 即溫度對流過程可表示為 速度的對流過程與溫度相似 可寫為 其中 是由z坐標轉(zhuǎn)化為 坐標后產(chǎn)生的彎曲項 28 29 垂直擴散項的計算 垂直擴散項的計算公式 第k層 1 k kb 1 可以改寫為 30 31 展開后合并同類項得到其中 32 33 現(xiàn)在假定溫度解的形式為 由 34 得到的代入到 32 就可以得到其中的值由 33 求得 的值由上一層的值確定 34 35 當(dāng)k 1時 即表層海水 海水溫度主要取決于海表面溫度通量 公式 34 的解可以近似寫為 這樣 表層海水的溫度只能根據(jù)式 31 求解 即 上式可以進一步表示為 36 37 式 37 還可以進一步寫為 再與溫度的通解 34 比較 可以得到 34 38 39 當(dāng)短波輻射通量的計算公式 r ad1 ad2均為光輻射常數(shù) 根據(jù)不同水質(zhì)取值如下 40 ntp12345JerlovtypeIIaIbIIIIIr0 580 620 670 700 78ad1 m 0 350 601 01 51 4ad2 m 23 020 017 014 07 9 來源 Jerlov 1976 PaulsonandSimpson 1977 當(dāng)k kb 1時 即底層海水 假定海底的熱通量為0 根據(jù)前面的推導(dǎo)方式同樣可以得到底層海水溫度的計算公式 對于鹽度方程垂向擴散項的計算與溫度方程相同 只是不考慮太陽短波輻射這一項 41 時間步長的CFL限制條件 由于模型水平方向采用顯格式 因此時間步長的選取必須要滿足CFL穩(wěn)定條件 對于外模態(tài) 時間步長的限制條件為其中 可能預(yù)見到的最大流速 對于內(nèi)模態(tài) 時間步長的限制條件較外模態(tài)的情形寬松很多 主要是速度較快的外重力波已經(jīng)在外模態(tài)中考慮了 一般取值30 50即可 42 時間步長的其它限制條件 對于動量或標量還有其它一些時間限制條件其中或 以及其中分別為地球角速率和地理緯度 43 44 側(cè)開邊界條件 bcond 陸地及岸線是由dum dvm fsm控制的 在陸地上這些變量的值設(shè)為0 有水的地方設(shè)為1子程序bcond idx idx 1對應(yīng)的水位邊界條件 idx 2對應(yīng)垂向平均流速邊界條件 idx 3對應(yīng)三維水平流速邊界條件 idx 4對應(yīng)溫鹽邊界條件 idx 5對應(yīng)垂向流速邊界條件 idx 6對應(yīng)湍流動能和湍流混合長度邊界條件模型對開邊界條件的要求很高 而開邊界條件本身又具有很大的不確定性 因此有必要對模型的外模態(tài)和內(nèi)模態(tài)的開邊界分別進行處理 數(shù)值積分流程圖 外模態(tài)開邊界條件一 外模態(tài)開邊界條件二 流速水位 流速 外模態(tài)開邊界條件三 水位 內(nèi)模態(tài)開邊界條件 一 內(nèi)模態(tài)開邊界條件 二 第三節(jié)POM模型的實際應(yīng)用 渤海河流羽狀流的數(shù)值模擬 渤海及入海河流 河流徑流輸入模型 將河流視為點源 而點源的模擬采用的是一個水平網(wǎng)格的面積徑流動量對動量方程的貢獻可以忽略不計受河流徑流影響的二維連續(xù)方程 三維連續(xù)方程可表示為內(nèi)外模式相容 30天時純浮力羽流條件下渤海的鹽度分布 等值線間隔為0 5psu 最外層等值線為30psu a 海水表層 b 海水底層 30天時潮致羽流條件下渤海的鹽度分布 等值線間隔為0 5psu 最外層等值線為30psu a 海水表層 b 海水底層 東風(fēng) 南風(fēng) 風(fēng)場加入10