CN111859748B - 一种基于垂向混合坐标的海洋内波模拟方法 - Google Patents

一种基于垂向混合坐标的海洋内波模拟方法 Download PDF

Info

Publication number
CN111859748B
CN111859748B CN202010681033.2A CN202010681033A CN111859748B CN 111859748 B CN111859748 B CN 111859748B CN 202010681033 A CN202010681033 A CN 202010681033A CN 111859748 B CN111859748 B CN 111859748B
Authority
CN
China
Prior art keywords
vertical
equation
formula
term
coordinates
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010681033.2A
Other languages
English (en)
Other versions
CN111859748A (zh
Inventor
胡旭辉
张学峰
王奭璠
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tianjin University
Original Assignee
Tianjin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tianjin University filed Critical Tianjin University
Priority to CN202010681033.2A priority Critical patent/CN111859748B/zh
Publication of CN111859748A publication Critical patent/CN111859748A/zh
Application granted granted Critical
Publication of CN111859748B publication Critical patent/CN111859748B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于垂向混合坐标的海洋内波模拟方法：(1)构建基于垂向混合坐标的非静力海洋模式控制方程组；(2)根据构建的基于垂向混合坐标的非静力海洋模式控制方程组建立基于垂向混合坐标的非静力海洋模式；并且，获取内波相关数据，对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理；(3)采用建立的基于垂向混合坐标的非静力海洋模式模拟内波，获得内波数据；(4)对获得的内波数据进行分析，用于对海洋内波进行预报。本发明把z坐标和σ坐标统一到一个非静力数值模型中，在保持各自优点的同时，弥补各自的缺陷，可提高利用非静力海洋数值模式模拟预报内波的精度。

Description

一种基于垂向混合坐标的海洋内波模拟方法
技术领域
本发明涉及海洋环境工程，特别涉及一种基于垂向混合坐标的海洋内波模拟方法。
背景技术
内波是发生在海水密度层结稳定的海洋中的一种波动，是海洋中普遍存在的一种动力学过程，其对海洋混合、全球的温盐环流起着重要作用，对海洋生态环境、海洋工程和军事等有重大影响。目前通常利用非静力海洋模型对海洋内波进行数值模拟，分析物理机制。但非静力近似计算目前只应用于σ坐标和z坐标，其使用范围受到限制。常用的z坐标，虽然对海洋上混合层和密度跃层刻画精度较高，但其在海底边界层存在阶梯效应，会造成虚假混合和虚假流动，使通量失衡。常用的σ坐标，虽然能很好地拟合复杂地形，但在陡峭的陆坡地形会造成斜压梯度力误差。
发明内容
本发明的目的在于避免使用单一垂向坐标系的非静力海洋模型模拟内波时难以模拟底部边界层或在陡峭地形和粗分辨率时造成斜压梯度力误差的问题，提供一种在z坐标和σ坐标组成的垂向混合坐标系下模拟海洋内波的方法，本发明方法把z坐标和σ坐标统一到一个非静力数值模型中，在保持各自优点的同时，弥补各自的缺陷，可提高利用非静力海洋数值模式模拟预报内波的精度。
本发明所采用的技术方案是：一种基于垂向混合坐标的海洋内波模拟方法，包括以下步骤：
步骤1，构建基于垂向混合坐标的非静力海洋模式控制方程组；
步骤2，根据步骤1构建的基于垂向混合坐标的非静力海洋模式控制方程组建立基于垂向混合坐标的非静力海洋模式；并且，获取内波相关数据，对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理，其中，所述的内波相关数据包括研究区域的地形、温度、盐度、密度、风场、降水和正压潮；
步骤3，采用步骤2所建立的基于垂向混合坐标的非静力海洋模式模拟内波，获得内波数据，其中，所述的内波数据包括各方向流速、海温、密度、动压和水位、湍动能、湍流混合长度、垂向湍粘性系数和垂向湍扩散系数；
步骤4，对步骤3获得的内波数据进行分析，用于对海洋内波进行预报。
其中，步骤1进一步包括：
基于任意因变量在z坐标系和垂向混合坐标系间的转换关系，得到基于垂向混合坐标的非静力海洋模式控制方程组，所述基于垂向混合坐标的非静力海洋模式控制方程组包括连续性方程、动量方程、温盐方程、湍方程和状态方程：
(1)连续性方程有两种形式，如式(1)和式(2)所示：
式中，u为x方向的速度分量；v为y方向的速度分量；w为z坐标系中z方向的速度分量；ω为垂向混合坐标系中垂向速度分量，ω＝w-(sx+ηx)u-(sy+ηy)v-(st+ηt)；
η为水位；s为垂向混合坐标系中的垂向坐标；k为垂向网格的编号；x、y为水平坐标，z为z坐标系中的垂向坐标，t为时间；
(2)动量方程包括水平动量方程和垂向动量方程，所述水平动量方程如式(3)和式(4)所示，所述垂向动量方程如式(5)所示：
式中，f为科氏参数；g为重力加速度；ρ0为参考密度；ρ′为海水扰动密度；dk′为积分微元；Q为动压；KM为垂向湍粘性系数；Fx->s为u的水平粘性项；Fy->s为v的水平粘性项；
(3)温盐方程如式(6)和式(7)所示：
式中，T为海温；KH为垂向湍扩散系数；FT->s为温度的水平湍扩散项；R为短波辐射通量；S为盐度；FS->s为盐度的水平湍扩散项；
(4)湍方程如式(8)和式(9)所示：
式中，q2为湍动能的2倍；Kq为垂向湍流混合系数；为绝热递减率修正后的垂向密度梯度；/>l为湍流混合长度；Fq->s为湍动能的水平湍扩散项；B1、E1、E3均为经验参数；/>为面壁近似函数；Fl->s为湍混合长的水平湍扩散项；
(5)状态方程如式(10)所示：
ρ＝ρ(T,S,P) (10)
式中，ρ为海水密度；P为压强。
其中，步骤2中，所述的对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理包括：
确定所要研究的区域，对该研究区域的地形进行模式网格化处理；
查找该研究区域的温度、盐度、密度分布资料，设置海洋模式的初始温度场、初始盐度场和初始密度场；查找该研究区域的风场、降水资料，设置海洋模式的强迫场；根据实际情况设置海洋模式的边界场；查找该研究区域的正压潮资料，设置海洋模式的初始速度场；根据内波模拟要求和实际情况，设置海洋模式的基本参数，所述海洋模式的基本参数包括底粗糙度、时间步长、积分步数、科氏参数和湍混合参数。
其中，步骤3进一步包括：
基于垂向混合坐标的非静力海洋模式利用分列算子将基于垂向混合坐标的非静力海洋模式控制方程组中不同的物理过程分开计算：
步骤3-1，显格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的平流项、科氏力项、斜压项和水平扩散项；
步骤3-2，隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项；
步骤3-3，隐格式求解水位；
步骤3-4，求解垂向混合坐标系中水位压强梯度力作用后的垂向速度和z坐标系中水位压强梯度力作用后的垂向速度/>
步骤3-5，隐格式求解动压；
步骤3-6，求解下一时刻的x方向的速度分量un+1、下一时刻的y方向的速度分量vn +1、下一时刻的z坐标系中z方向的速度分量wn+1以及下一时刻的垂向混合坐标系中垂向速度分量ωn+1；
步骤3-7，求解基于垂向混合坐标的非静力海洋模式控制方程组中的湍方程式(8)和式(9)；
步骤3-8，求解基于垂向混合坐标的非静力海洋模式控制方程组中的温盐方程式(6)和式(7)；
步骤3-9，求解基于垂向混合坐标的非静力海洋模式控制方程组中的状态方程式(10)；
步骤3-10，重复步骤3-1至步骤3-9直至满足内波模拟所需要的总的时间积分步数，输出内波数据。
其中，步骤3-1中，所述的显格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的平流项、科氏力项、斜压项和水平粘性项如式(11)和式(12)所示：
式中，(usk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(usk)的中间状态；(vsk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(vsk)的中间状态；n-1为第n-1时刻；Δt为时间步长；n为第n时刻；(usk)n-1为积分第n-1时刻的(usk)的值；为积分第n-1时刻u的水平粘性项；(vsk)n-1为积分第n-1时刻的(vsk)的值；/>为积分第n-1时刻v的水平粘性项；
步骤3-2中，所述的隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项如式(13)和式(14)所示：
式中，为垂向扩散项作用后的u的中间状态；/>为垂向扩散项作用后的v的中间状态；/>为垂向扩散项作用后的sk的中间状态；
步骤3-3中，所述的隐格式求解水位包括：
首先，给出水位求解方程如式(15)和式(16)所示：
式中，为水位压强梯度力作用后的u的中间状态；/>为水位压强梯度力作用后的v的中间状态；/>为水位压强梯度力作用后的sk的中间状态；
其次，先对式(15)和式(16)进行差分再垂向积分，可得到两个垂向积分式子；再将这两个垂向积分式子带入到式(2)垂向积分形式的差分格式中；整理该差分格式后化为五点隐格式，且关于水位的系数矩阵是严格对角占优的，用超松弛迭代方法得到水位的收敛解；
步骤3-4中，所述的求解垂向混合坐标系中水位压强梯度力作用后的垂向速度和z坐标系中水位压强梯度力作用后的垂向速度/>包括：
首先利用式(2)求解垂向混合坐标系中水位压强梯度力作用后的垂向速度然后根据式子w＝ω+(sx+ηx)u+(sy+ηy)v+(st+ηt)进行差分，计算得到z坐标系中水位压强梯度力作用后的垂向速度/>
步骤3-5中，所述的隐格式求解动压可采用15对角线求解法或7对角线求解法求解动压；
步骤3-7中，所述的求解基于垂向混合坐标的非静力海洋模式控制方程组中的湍方程式(8)和式(9)为：湍方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解；
步骤3-8中，所述的求解基于垂向混合坐标的非静力海洋模式控制方程组中的温盐方程式(6)和式(7)为：温盐方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解。
其中，步骤3-5中，所述的15对角线求解法求解动压包括：
对连续性方程(1)左右同除sk得到变形后的连续性方程如式(17)所示：
以动压Q为中心，对变形后的连续性方程(17)进行差分，得到的结果如式(18)至式(22)所示：
式中，dx为x方向的网格间距；dy为y方向的网格间距；dz为垂向相邻两个网格面的间距；dzz为垂向相邻两个网格中心的间距； i为x方向的网格编号，j为y方向的网格编号；
对式(18)至式(22)进行整理：式(18)+式(19)-式(20)+式(21)-式(22)，得：
其中，
式(23)为15点隐式方程，该15点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(23)时，得到动压的收敛解。
其中，步骤3-5中，所述的7对角线求解法求解动压包括：
对于w方程，忽略掉平流项：
将式(40)、式(41)和式(42)进行差分，然后带入式(1)半隐格式的差分形式中，经整理得：
其中，
式中，ZZ＝s+η；
式(43)为7点隐式方程，该7点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(43)时，得到动压的收敛解。
本发明的有益效果是：
本发明基于z坐标和σ坐标组成的垂向混合坐标系构建了一种非静力近似方法，应用于海洋数值模式，对海洋内波进行模拟计算。相比于传统的基于σ坐标和基于z坐标的非静力海洋模型对内波的模拟效果，本发明改善了海底边界层虚假混合、虚假流动以及陡峭陆坡地形处斜压梯度力误差的问题，提高了利用非静力海洋数值模式模拟预报内波的精度。
附图说明
图1：本发明一种基于垂向混合坐标的海洋内波模拟方法流程示意图。
具体实施方式
为能进一步了解本发明的发明内容、特点及功效，兹例举以下实施例，并配合附图详细说明如下：
为了解决非静力海洋模型中常用垂向坐标系存在的问题，提高海洋内波模拟的精度，可采用垂向混合坐标系，把z坐标和σ坐标统一到一个非静力数值模型中，在保持各自优点的同时，弥补各自的缺陷。因此，针对垂向混合坐标的优势和海洋内波模拟的现状，本发明提出了一种基于垂向混合坐标的海洋内波模拟方法。
如附图1所示，一种基于垂向混合坐标的海洋内波模拟方法，包括以下步骤，其中，垂向混合坐标系下考虑非静力的控制方程组求解算法是本发明的重要内容。
步骤1，基于任意因变量在z坐标系和垂向混合坐标系间的转换关系，得到基于垂向混合坐标的非静力海洋模式控制方程组，所述基于垂向混合坐标的非静力海洋模式控制方程组包括连续性方程、动量方程、温盐方程、湍方程和状态方程：
(1)连续性方程有两种形式，如式(1)和式(2)所示：
式中，u为x方向的速度分量；v为y方向的速度分量；w为z坐标系中z方向的速度分量；ω为垂向混合坐标系中垂向速度分量，ω＝w-(sx+ηx)u-(sy+ηy)v-(st+ηt)；
η为水位；s为垂向混合坐标系中的垂向坐标；k为垂向网格的编号；x、y为水平坐标，z为z坐标系中的垂向坐标，t为时间。
(2)动量方程包括水平动量方程和垂向动量方程，所述水平动量方程如式(3)和式(4)所示，所述垂向动量方程如式(5)所示：
/>
式中，f为科氏参数；g为重力加速度；ρ0为参考密度；ρ′为海水扰动密度；dk′为积分微元；Q为动压；KM为垂向湍粘性系数；Fx->s为u的水平粘性项；Fy->s为v的水平粘性项。
(3)温盐方程如式(6)和式(7)所示：
式中，T为海温；KH为垂向湍扩散系数；FT->s为温度的水平湍扩散项；R为短波辐射通量；S为盐度；FS->s为盐度的水平湍扩散项。
(4)湍方程如式(8)和式(9)所示：
式中，q2为湍动能的2倍；Kq为垂向湍流混合系数；为绝热递减率修正后的垂向密度梯度；/>l为湍流混合长度；Fq->s为湍动能的水平湍扩散项；B1、E1、E3均为经验参数；W～为面壁近似函数；Fl->s为湍混合长的水平湍扩散项；垂向动量方程暂时仅保留了全导数项和非静力项。
(5)状态方程如式(10)所示：
ρ＝ρ(T,S,P) (10)
式中，ρ为海水密度；P为压强。
步骤2，根据步骤1构建的基于垂向混合坐标的非静力海洋模式控制方程组建立基于垂向混合坐标的非静力海洋模式，即，将基于垂向混合坐标的非静力海洋模式控制方程组离散化处理再用计算机语言编程。并且，获取内波相关数据，对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理，其中，所述的内波相关数据包括研究区域的地形、温度、盐度、密度、风场、降水和正压潮。
首先，确定所要研究的区域，对该研究区域的地形进行模式网格化处理。其次，查找该研究区域的温度、盐度、密度分布资料，设置海洋模式的初始温度场、初始盐度场和初始密度场；查找该研究区域的风场、降水资料，设置海洋模式的强迫场；根据实际情况设置海洋模式的边界场；查找该研究区域的正压潮资料，设置海洋模式的初始速度场；根据内波模拟要求和实际情况，设置海洋模式的基本参数，所述海洋模式的基本参数包括底粗糙度、时间步长、积分步数、科氏参数和湍混合参数。
步骤3，采用步骤2所建立的基于垂向混合坐标的非静力海洋模式模拟内波，获得内波数据，其中，所述的内波数据包括各方向流速、海温、密度、动压和水位、湍动能、湍流混合长度、垂向湍粘性系数和垂向湍扩散系数。
基于垂向混合坐标的非静力海洋模式利用分列算子将基于垂向混合坐标的非静力海洋模式控制方程组中不同的物理过程分开计算：
步骤3-1，显格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的平流项、科氏力项、斜压项和水平扩散项，如式(11)和式(12)所示：
式中，(usk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(usk)的中间状态；(vsk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(vsk)的中间状态；n-1为第n-1时刻；Δt为时间步长；n为第n时刻；(usk)n-1为积分第n-1时刻的(usk)的值；为积分第n-1时刻u的水平粘性项；(vsk)n-1为积分第n-1时刻的(vsk)的值；/>为积分第n-1时刻v的水平粘性项。
步骤3-2，隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项，如式(13)和式(14)所示：
式中，为垂向扩散项作用后的u的中间状态；/>为垂向扩散项作用后的v的中间状态；/>为垂向扩散项作用后的sk的中间状态。
步骤3-3，隐格式求解水位：
首先，给出水位求解方程如式(15)和式(16)所示：
式中，为水位压强梯度力作用后的u的中间状态；/>为水位压强梯度力作用后的v的中间状态；/>为水位压强梯度力作用后的sk的中间状态；
其次，先对式(15)和式(16)进行差分再垂向积分，可得到两个垂向积分式子；再将这两个垂向积分式子带入到式(2)垂向积分形式的差分格式中；整理该差分格式后化为五点隐格式，且关于水位的系数矩阵是严格对角占优的，用超松弛迭代方法得到水位的收敛解。其中，所谓的严格对角占优是指矩阵中每个主对角元素的模都大于与它同行的其他元素的模的总和；对列同样成立。
步骤3-4，求解垂向混合坐标系中水位压强梯度力作用后的垂向速度和z坐标系中水位压强梯度力作用后的垂向速度/>
首先利用式(2)求解垂向混合坐标系中水位压强梯度力作用后的垂向速度然后根据式子w＝ω+(sx+ηx)u+(sy+ηy)v+(st+ηt)进行差分，计算得到z坐标系中水位压强梯度力作用后的垂向速度/>
步骤3-5，隐格式求解动压：
方法一：15对角线求解法
对连续性方程(1)左右同除sk得到变形后的连续性方程如式(17)所示：
以动压Q为中心，对变形后的连续性方程(17)进行差分，得到的结果如式(18)至式(22)所示：
/>
式中，dx为x方向的网格间距；dy为y方向的网格间距；dz为垂向相邻两个网格面的间距；dzz为垂向相邻两个网格中心的间距； i为x方向的网格编号，j为y方向的网格编号；
对式(18)至式(22)进行整理：式(18)+式(19)-式(20)+式(21)-式(22)，得：
其中，
/>
式(23)为15点隐式方程，该15点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(23)时，可得到动压的收敛解。
方法二：7对角线求解法
对于w方程，忽略掉平流项：
将式(40)、式(41)和式(42)进行差分，然后带入式(1)半隐格式的差分形式中，经整理得：
其中，
式中，ZZ＝s+η；
式(43)为7点隐式方程，该7点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(43)时，可得到动压的收敛解。
步骤3-6，求解下一时刻的x方向的速度分量un+1、下一时刻的y方向的速度分量vn +1、下一时刻的z坐标系中z方向的速度分量wn+1以及下一时刻的垂向混合坐标系中垂向速度分量ωn+1。
步骤3-7，求解基于垂向混合坐标的非静力海洋模式控制方程组中的湍方程式(8)和式(9)。
其中，湍方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解。
步骤3-8，求解基于垂向混合坐标的非静力海洋模式控制方程组中的温盐方程式(6)和式(7)。
其中，温盐方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解。
步骤3-9，求解基于垂向混合坐标的非静力海洋模式控制方程组中的状态方程式(10)。
步骤3-10，重复步骤3-1至步骤3-9直至满足内波模拟所需要的总的时间积分步数，输出内波数据。
步骤4，对步骤3获得的内波数据进行分析，用于对海洋内波进行预报。
针对基于垂向混合坐标的非静力海洋模式输出的各时刻变量数据进行处理，分析内波变化机制，对海洋内波进行预报。
尽管上面结合附图对本发明的优选实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，并不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可以做出很多形式，这些均属于本发明的保护范围之内。

Claims (5)

1.一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，包括以下步骤：
步骤1，构建基于垂向混合坐标的非静力海洋模式控制方程组，包括：
基于任意因变量在z坐标系和垂向混合坐标系间的转换关系，得到基于垂向混合坐标的非静力海洋模式控制方程组，所述基于垂向混合坐标的非静力海洋模式控制方程组包括连续性方程、动量方程、温盐方程、湍方程和状态方程：
(1)连续性方程有两种形式，如式(1)和式(2)所示：
式中，u为x方向的速度分量；v为y方向的速度分量；w为z坐标系中z方向的速度分量；ω为垂向混合坐标系中垂向速度分量，ω＝w-(sx+ηx)u-(sy+ηy)v-(st+ηt)；
η为水位；s为垂向混合坐标系中的垂向坐标；k为垂向网格的编号；x、y为水平坐标，z为z坐标系中的垂向坐标，t为时间；
(2)动量方程包括水平动量方程和垂向动量方程，所述水平动量方程如式(3)和式(4)所示，所述垂向动量方程如式(5)所示：
式中，f为科氏参数；g为重力加速度；ρ0为参考密度；ρ′为海水扰动密度；dk′为积分微元；Q为动压；KM为垂向湍粘性系数；Fx->s为u的水平粘性项；Fy->s为v的水平粘性项；
(3)温盐方程如式(6)和式(7)所示：
式中，T为海温；KH为垂向湍扩散系数；FT->s为温度的水平湍扩散项；R为短波辐射通量；S为盐度；FS->s为盐度的水平湍扩散项；
(4)湍方程如式(8)和式(9)所示：
式中，q2为湍动能的2倍；Kq为垂向湍流混合系数；为绝热递减率修正后的垂向密度梯度；q3为/>l为湍流混合长度；Fq->s为湍动能的水平湍扩散项；B1、E1、E3均为经验参数；/>为面壁近似函数；Fl->s为湍混合长的水平湍扩散项；
(5)状态方程如式(10)所示：
ρ＝ρ(T,S,P) (10)
式中，ρ为海水密度；P为压强；
步骤2，根据步骤1构建的基于垂向混合坐标的非静力海洋模式控制方程组建立基于垂向混合坐标的非静力海洋模式；并且，获取内波相关数据，对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理，其中，所述的内波相关数据包括研究区域的地形、温度、盐度、密度、风场、降水和正压潮；
步骤3，采用步骤2所建立的基于垂向混合坐标的非静力海洋模式模拟内波，获得内波数据，其中，所述的内波数据包括各方向流速、海温、密度、动压和水位、湍动能、湍流混合长度、垂向湍粘性系数和垂向湍扩散系数；
基于垂向混合坐标的非静力海洋模式利用分列算子将基于垂向混合坐标的非静力海洋模式控制方程组中不同的物理过程分开计算：
步骤3-1，显格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的平流项、科氏力项、斜压项和水平扩散项；
步骤3-2，隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项；
步骤3-3，隐格式求解水位；
步骤3-4，求解垂向混合坐标系中水位压强梯度力作用后的垂向速度和z坐标系中水位压强梯度力作用后的垂向速度/>
步骤3-5，隐格式求解动压；
步骤3-6，求解下一时刻的x方向的速度分量un+1、下一时刻的y方向的速度分量vn+1、下一时刻的z坐标系中z方向的速度分量wn+1以及下一时刻的垂向混合坐标系中垂向速度分量ωn+1；
步骤3-7，求解基于垂向混合坐标的非静力海洋模式控制方程组中的湍方程式(8)和式(9)；
步骤3-8，求解基于垂向混合坐标的非静力海洋模式控制方程组中的温盐方程式(6)和式(7)；
步骤3-9，求解基于垂向混合坐标的非静力海洋模式控制方程组中的状态方程式(10)；
步骤3-10，重复步骤3-1至步骤3-9直至满足内波模拟所需要的总的时间积分步数，输出内波数据；
步骤4，对步骤3获得的内波数据进行分析，用于对海洋内波进行预报。
2.根据权利要求1所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤2中，所述的对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理包括：
确定所要研究的区域，对该研究区域的地形进行模式网格化处理；
查找该研究区域的温度、盐度、密度分布资料，设置海洋模式的初始温度场、初始盐度场和初始密度场；查找该研究区域的风场、降水资料，设置海洋模式的强迫场；根据实际情况设置海洋模式的边界场；查找该研究区域的正压潮资料，设置海洋模式的初始速度场；根据内波模拟要求和实际情况，设置海洋模式的基本参数，所述海洋模式的基本参数包括底粗糙度、时间步长、积分步数、科氏参数和湍混合参数。
3.根据权利要求1所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤3-1中，所述的显格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的平流项、科氏力项、斜压项和水平粘性项如式(11)和式(12)所示：
式中，(usk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(usk)的中间状态；(vsk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(vsk)的中间状态；n-1为第n-1时刻；△t为时间步长；n为第n时刻；(usk)n-1为积分第n-1时刻的(usk)的值；为积分第n-1时刻u的水平粘性项；(vsk)n-1为积分第n-1时刻的(vsk)的值；/>为积分第n-1时刻v的水平粘性项；
步骤3-2中，所述的隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项如式(13)和式(14)所示：
式中，为垂向扩散项作用后的u的中间状态；/>为垂向扩散项作用后的v的中间状态；为垂向扩散项作用后的sk的中间状态；
步骤3-3中，所述的隐格式求解水位包括：
首先，给出水位求解方程如式(15)和式(16)所示：
式中，为水位压强梯度力作用后的u的中间状态；/>为水位压强梯度力作用后的v的中间状态；/>为水位压强梯度力作用后的sk的中间状态；
其次，先对式(15)和式(16)进行差分再垂向积分，可得到两个垂向积分式子；再将这两个垂向积分式子带入到式(2)垂向积分形式的差分格式中；整理该差分格式后化为五点隐格式，且关于水位的系数矩阵是严格对角占优的，用超松弛迭代方法得到水位的收敛解；
步骤3-4中，所述的求解垂向混合坐标系中水位压强梯度力作用后的垂向速度和z坐标系中水位压强梯度力作用后的垂向速度/>包括：
首先利用式(2)求解垂向混合坐标系中水位压强梯度力作用后的垂向速度然后根据式子w＝ω+(sx+ηx)u+(sy+ηy)v+(st+ηt)进行差分，计算得到z坐标系中水位压强梯度力作用后的垂向速度/>
步骤3-5中，所述的隐格式求解动压可采用15对角线求解法或7对角线求解法求解动压；
步骤3-7中，所述的求解基于垂向混合坐标的非静力海洋模式控制方程组中的湍方程式(8)和式(9)为：湍方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解；
步骤3-8中，所述的求解基于垂向混合坐标的非静力海洋模式控制方程组中的温盐方程式(6)和式(7)为：温盐方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解。
4.根据权利要求3所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤3-5中，所述的15对角线求解法求解动压包括：
对连续性方程(1)左右同除sk得到变形后的连续性方程如式(17)所示：
以动压Q为中心，对变形后的连续性方程(17)进行差分，得到的结果如式(18)至式(22)所示：
式中，dx为x方向的网格间距；dy为y方向的网格间距；dz为垂向相邻两个网格面的间距；dzz为垂向相邻两个网格中心的间距； i为x方向的网格编号，j为y方向的网格编号；
对式(18)至式(22)进行整理：式(18)+式(19)-式(20)+式(21)-式(22)，得：
其中，
式(23)为15点隐式方程，该15点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(23)时，得到动压的收敛解。
5.根据权利要求4所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤3-5中，所述的7对角线求解法求解动压包括：
对于w方程，忽略掉平流项：
将式(40)、式(41)和式(42)进行差分，然后带入式(1)半隐格式的差分形式中，经整理得：
其中，
式中，ZZ＝s+η；
式(43)为7点隐式方程，该7点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(43)时，得到动压的收敛解。
CN202010681033.2A 2020-07-15 2020-07-15 一种基于垂向混合坐标的海洋内波模拟方法 Active CN111859748B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010681033.2A CN111859748B (zh) 2020-07-15 2020-07-15 一种基于垂向混合坐标的海洋内波模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010681033.2A CN111859748B (zh) 2020-07-15 2020-07-15 一种基于垂向混合坐标的海洋内波模拟方法

Publications (2)

Publication Number Publication Date
CN111859748A CN111859748A (zh) 2020-10-30
CN111859748B true CN111859748B (zh) 2024-04-02

Family

ID=72984271

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010681033.2A Active CN111859748B (zh) 2020-07-15 2020-07-15 一种基于垂向混合坐标的海洋内波模拟方法

Country Status (1)

Country Link
CN (1) CN111859748B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112612027B (zh) * 2020-12-15 2021-09-10 中国科学院声学研究所 一种浅海环境下利用声能量起伏的海洋内波监测方法
CN113468679B (zh) * 2021-09-06 2021-11-12 中国空气动力研究与发展中心计算空气动力研究所 一种基于s-a模型的湍流长度尺度计算方法
CN114814276B (zh) * 2022-03-21 2023-08-18 汕头大学 海上风电设备运行导致周边海水垂向运动流速的计算方法
CN115391069B (zh) * 2022-10-27 2023-02-03 山东省计算中心（国家超级计算济南中心） 基于海洋模式roms的并行通讯方法及系统
CN116341408A (zh) * 2023-03-13 2023-06-27 中国人民解放军国防科技大学 一种距离相关水声学应用的海洋模式数据坐标转换方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20080052718A (ko) * 2006-12-08 2008-06-12 재단법인서울대학교산학협력재단 해양 방류시스템의 근역 및 원역 결합모의방법
KR20150117971A (ko) * 2014-04-11 2015-10-21 한국해양과학기술원 연안국지 해수순환 및 파랑 예측 방법 및 시스템
CN108562896A (zh) * 2018-04-23 2018-09-21 广西师范大学 一种基于三维正压浅海大陆架模型的深层海流反演方法
CN110096792A (zh) * 2019-04-28 2019-08-06 天津大学 一种计算非定常Langmuir环流的动态模拟方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20080052718A (ko) * 2006-12-08 2008-06-12 재단법인서울대학교산학협력재단 해양 방류시스템의 근역 및 원역 결합모의방법
KR20150117971A (ko) * 2014-04-11 2015-10-21 한국해양과학기술원 연안국지 해수순환 및 파랑 예측 방법 및 시스템
CN108562896A (zh) * 2018-04-23 2018-09-21 广西师范大学 一种基于三维正压浅海大陆架模型的深层海流反演方法
CN110096792A (zh) * 2019-04-28 2019-08-06 天津大学 一种计算非定常Langmuir环流的动态模拟方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
热带太平洋地区HYCOM模式不同垂向坐标设置对比研究;陈晓斌;周林;刘志宏;陈璇;厦门大学学报. 自然科学版;第53卷(第6期);全文 *

Also Published As

Publication number Publication date
CN111859748A (zh) 2020-10-30

Similar Documents

Publication Publication Date Title
CN111859748B (zh) 一种基于垂向混合坐标的海洋内波模拟方法
Parsons et al. Numerical modelling of flow structures over idealized transverse aeolian dunes of varying geometry
Chorley Models in geomorphology
Bristow et al. Turbulent flow structure associated with collision between laterally offset, fixed‐bed barchan dunes
CN106844859A (zh) 一种滑坡涌浪计算方法
Trowbridge A simple description of the deepening and structure of a stably stratified flow driven by a surface stress
CN115525978B (zh) 用于飞行器波浪情况水动力分析的速度入口数值造波方法
Zheng et al. A 3D unstructured grid nearshore hydrodynamic model based on the vortex force formalism
Simmons et al. Verification of a numerical ocean model of the Arabian Sea
Jin et al. Quasi-3D numerical modeling of shallow-water circulation
CN112765854B (zh) 一种路面内部裂缝数量预测方法
Wyngaard Large-eddy simulation: Guidelines for its application to planetary boundary layer research
Kanehira et al. The effects of smoothing length on the onset of wave breaking in smoothed particle hydrodynamics (SPH) simulations of highly directionally spread waves
Calantoni et al. Simulation of sediment motions using a discrete particle model in the inner surf and swash-zones
Schwab et al. Lagrangian comparison of objectively analyzed and dynamically modeled circulation patterns in Lake Erie
Guo et al. A new depth-integrated non-hydrostatic model for free surface flows
Price et al. Upper ocean dynamics
Jemison et al. A Prototype Incompressible Pressure-Based Solver for Free-Surface Flows in CREATETM-AV Kestrel
Kuebel Cervantes et al. Lagrangian characteristics of continental shelf flows forced by periodic wind stress
Roy et al. Large eddy simulation of wind flow over complex terrain: The Bolund Hill case
Skyllingstad et al. Nonlinear vertical mixing processes in the ocean: modeling and parameterization
Huang Enhancement of a Turbulence Sub-Model for More Accurate Predictions of Vertical Stratifications in 3D Coastal and Estuarine Modeling
XU et al. Three-dimensional modeling of hydrodynamics and dissolved oxygen transport in Tone River Estuary
Bao An improved immersed boundary method for atmospheric boundary layer simulations over complex terrain
Wen‐Bin et al. Application of Grid‐Nesting Technique on Sandwaves Migration Simulation II—Sandwaves Migration in Northern South China Sea

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant