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

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

Info

Publication number
CN111859748A
CN111859748A CN202010681033.2A CN202010681033A CN111859748A CN 111859748 A CN111859748 A CN 111859748A CN 202010681033 A CN202010681033 A CN 202010681033A CN 111859748 A CN111859748 A CN 111859748A
Authority
CN
China
Prior art keywords
vertical
equation
formula
term
static
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.)
Granted
Application number
CN202010681033.2A
Other languages
English (en)
Other versions
CN111859748B (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

Images

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)所示：
Figure BDA0002585827680000021
Figure BDA0002585827680000022
式中，u为x方向的速度分量；v为y方向的速度分量；w为z坐标系中z方向的速度分量；ω为垂向混合坐标系中垂向速度分量，ω＝w-(sx+ηx)u-(sy+ηy)v-(st+ηt)；
Figure BDA0002585827680000023
Figure BDA0002585827680000024
η为水位；s为垂向混合坐标系中的垂向坐标；k为垂向网格的编号；x、y为水平坐标，z为z坐标系中的垂向坐标，t为时间；
(2)动量方程包括水平动量方程和垂向动量方程，所述水平动量方程如式(3)和式(4)所示，所述垂向动量方程如式(5)所示：
Figure BDA0002585827680000025
Figure BDA0002585827680000031
Figure BDA0002585827680000032
式中，f为科氏参数；g为重力加速度；ρ0为参考密度；ρ′为海水扰动密度；dk′为积分微元；Q为动压；KM为垂向湍粘性系数；Fx->s为u的水平粘性项；Fy->s为v的水平粘性项；
(3)温盐方程如式(6)和式(7)所示：
Figure BDA0002585827680000033
Figure BDA0002585827680000034
式中，T为海温；KH为垂向湍扩散系数；FT->s为温度的水平湍扩散项；R为短波辐射通量；S为盐度；FS->s为盐度的水平湍扩散项；
(4)湍方程如式(8)和式(9)所示：
Figure BDA0002585827680000035
Figure BDA0002585827680000036
式中，q2为湍动能的2倍；Kq为垂向湍流混合系数；
Figure BDA0002585827680000037
为绝热递减率修正后的垂向密度梯度；
Figure BDA0002585827680000038
l为湍流混合长度；Fq->s为湍动能的水平湍扩散项；B1、E1、E3均为经验参数；
Figure BDA0002585827680000039
为面壁近似函数；Fl->s为湍混合长的水平湍扩散项；
(5)状态方程如式(10)所示：
ρ＝ρ(T,S,P) (10)
式中，ρ为海水密度；P为压强。
其中，步骤2中，所述的对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理包括：
确定所要研究的区域，对该研究区域的地形进行模式网格化处理；
查找该研究区域的温度、盐度、密度分布资料，设置海洋模式的初始温度场、初始盐度场和初始密度场；查找该研究区域的风场、降水资料，设置海洋模式的强迫场；根据实际情况设置海洋模式的边界场；查找该研究区域的正压潮资料，设置海洋模式的初始速度场；根据内波模拟要求和实际情况，设置海洋模式的基本参数，所述海洋模式的基本参数包括底粗糙度、时间步长、积分步数、科氏参数和湍混合参数。
其中，步骤3进一步包括：
基于垂向混合坐标的非静力海洋模式利用分列算子将基于垂向混合坐标的非静力海洋模式控制方程组中不同的物理过程分开计算：
步骤3-1，显格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的平流项、科氏力项、斜压项和水平扩散项；
步骤3-2，隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项；
步骤3-3，隐格式求解水位；
步骤3-4，求解垂向混合坐标系中水位压强梯度力作用后的垂向速度
Figure BDA0002585827680000041
和z坐标系中水位压强梯度力作用后的垂向速度
Figure BDA0002585827680000042
步骤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)所示：
Figure BDA0002585827680000051
Figure BDA0002585827680000052
式中，(usk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(usk)的中间状态；(vsk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(vsk)的中间状态；n-1为第n-1时刻；Δt为时间步长；n为第n时刻；(usk)n-1为积分第n-1时刻的(usk)的值；
Figure BDA0002585827680000053
为积分第n-1时刻u的水平粘性项；(vsk)n-1为积分第n-1时刻的(vsk)的值；
Figure BDA0002585827680000054
为积分第n-1时刻v的水平粘性项；
步骤3-2中，所述的隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项如式(13)和式(14)所示：
Figure BDA0002585827680000055
Figure BDA0002585827680000061
式中，
Figure BDA0002585827680000062
为垂向扩散项作用后的u的中间状态；
Figure BDA0002585827680000063
为垂向扩散项作用后的v的中间状态；
Figure BDA0002585827680000064
为垂向扩散项作用后的sk的中间状态；
步骤3-3中，所述的隐格式求解水位包括：
首先，给出水位求解方程如式(15)和式(16)所示：
Figure BDA0002585827680000065
Figure BDA0002585827680000066
式中，
Figure BDA0002585827680000067
为水位压强梯度力作用后的u的中间状态；
Figure BDA0002585827680000068
为水位压强梯度力作用后的v的中间状态；
Figure BDA0002585827680000069
为水位压强梯度力作用后的sk的中间状态；
其次，先对式(15)和式(16)进行差分再垂向积分，可得到两个垂向积分式子；再将这两个垂向积分式子带入到式(2)垂向积分形式的差分格式中；整理该差分格式后化为五点隐格式，且关于水位的系数矩阵是严格对角占优的，用超松弛迭代方法得到水位的收敛解；
步骤3-4中，所述的求解垂向混合坐标系中水位压强梯度力作用后的垂向速度
Figure BDA00025858276800000610
和z坐标系中水位压强梯度力作用后的垂向速度
Figure BDA00025858276800000611
包括：
首先利用式(2)求解垂向混合坐标系中水位压强梯度力作用后的垂向速度
Figure BDA00025858276800000612
然后根据式子w＝ω+(sx+ηx)u+(sy+ηy)v+(st+ηt)进行差分，计算得到z坐标系中水位压强梯度力作用后的垂向速度
Figure BDA00025858276800000613
步骤3-5中，所述的隐格式求解动压可采用15对角线求解法或7对角线求解法求解动压；
步骤3-7中，所述的求解基于垂向混合坐标的非静力海洋模式控制方程组中的湍方程式(8)和式(9)为：湍方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解；
步骤3-8中，所述的求解基于垂向混合坐标的非静力海洋模式控制方程组中的温盐方程式(6)和式(7)为：温盐方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解。
其中，步骤3-5中，所述的15对角线求解法求解动压包括：
对连续性方程(1)左右同除sk得到变形后的连续性方程如式(17)所示：
Figure BDA0002585827680000071
以动压Q为中心，对变形后的连续性方程(17)进行差分，得到的结果如式(18)至式(22)所示：
Figure BDA0002585827680000072
Figure BDA0002585827680000073
Figure BDA0002585827680000074
Figure BDA0002585827680000081
Figure BDA0002585827680000082
Figure BDA0002585827680000083
式中，dx为x方向的网格间距；dy为y方向的网格间距；dz为垂向相邻两个网格面的间距；dzz为垂向相邻两个网格中心的间距；
Figure BDA0002585827680000084
Figure BDA0002585827680000085
i为x方向的网格编号，j为y方向的网格编号；
对式(18)至式(22)进行整理：式(18)+式(19)-式(20)+式(21)-式(22)，得：
Figure BDA0002585827680000091
其中，
Figure BDA0002585827680000092
Figure BDA0002585827680000093
Figure BDA0002585827680000094
Figure BDA0002585827680000095
Figure BDA0002585827680000096
Figure BDA0002585827680000097
Figure BDA0002585827680000098
Figure BDA0002585827680000099
Figure BDA0002585827680000101
Figure BDA0002585827680000102
Figure BDA0002585827680000103
Figure BDA0002585827680000104
Figure BDA0002585827680000105
Figure BDA0002585827680000106
Figure BDA0002585827680000107
Figure BDA0002585827680000108
式(23)为15点隐式方程，该15点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(23)时，得到动压的收敛解。
其中，步骤3-5中，所述的7对角线求解法求解动压包括：
Figure BDA0002585827680000109
Figure BDA00025858276800001010
对于w方程，忽略掉平流项：
Figure BDA00025858276800001011
将式(40)、式(41)和式(42)进行差分，然后带入式(1)半隐格式的差分形式中，经整理得：
Figure BDA0002585827680000111
其中，
Figure BDA0002585827680000112
Figure BDA0002585827680000113
Figure BDA0002585827680000114
Figure BDA0002585827680000115
Figure BDA0002585827680000116
Figure BDA0002585827680000117
Figure BDA0002585827680000118
Figure BDA0002585827680000119
式中，ZZ＝s+η；
式(43)为7点隐式方程，该7点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(43)时，得到动压的收敛解。
本发明的有益效果是：
本发明基于z坐标和σ坐标组成的垂向混合坐标系构建了一种非静力近似方法，应用于海洋数值模式，对海洋内波进行模拟计算。相比于传统的基于σ坐标和基于z坐标的非静力海洋模型对内波的模拟效果，本发明改善了海底边界层虚假混合、虚假流动以及陡峭陆坡地形处斜压梯度力误差的问题，提高了利用非静力海洋数值模式模拟预报内波的精度。
附图说明
图1：本发明一种基于垂向混合坐标的海洋内波模拟方法流程示意图。
具体实施方式
为能进一步了解本发明的发明内容、特点及功效，兹例举以下实施例，并配合附图详细说明如下：
为了解决非静力海洋模型中常用垂向坐标系存在的问题，提高海洋内波模拟的精度，可采用垂向混合坐标系，把z坐标和σ坐标统一到一个非静力数值模型中，在保持各自优点的同时，弥补各自的缺陷。因此，针对垂向混合坐标的优势和海洋内波模拟的现状，本发明提出了一种基于垂向混合坐标的海洋内波模拟方法。
如附图1所示，一种基于垂向混合坐标的海洋内波模拟方法，包括以下步骤，其中，垂向混合坐标系下考虑非静力的控制方程组求解算法是本发明的重要内容。
步骤1，基于任意因变量在z坐标系和垂向混合坐标系间的转换关系，得到基于垂向混合坐标的非静力海洋模式控制方程组，所述基于垂向混合坐标的非静力海洋模式控制方程组包括连续性方程、动量方程、温盐方程、湍方程和状态方程：
(1)连续性方程有两种形式，如式(1)和式(2)所示：
Figure BDA0002585827680000121
Figure BDA0002585827680000131
式中，u为x方向的速度分量；v为y方向的速度分量；w为z坐标系中z方向的速度分量；ω为垂向混合坐标系中垂向速度分量，ω＝w-(sx+ηx)u-(sy+ηy)v-(st+ηt)；
Figure BDA0002585827680000132
Figure BDA0002585827680000133
η为水位；s为垂向混合坐标系中的垂向坐标；k为垂向网格的编号；x、y为水平坐标，z为z坐标系中的垂向坐标，t为时间。
(2)动量方程包括水平动量方程和垂向动量方程，所述水平动量方程如式(3)和式(4)所示，所述垂向动量方程如式(5)所示：
Figure BDA0002585827680000134
Figure BDA0002585827680000135
Figure BDA0002585827680000136
式中，f为科氏参数；g为重力加速度；ρ0为参考密度；ρ′为海水扰动密度；dk′为积分微元；Q为动压；KM为垂向湍粘性系数；Fx->s为u的水平粘性项；Fy->s为v的水平粘性项。
(3)温盐方程如式(6)和式(7)所示：
Figure BDA0002585827680000137
Figure BDA0002585827680000141
式中，T为海温；KH为垂向湍扩散系数；FT->s为温度的水平湍扩散项；R为短波辐射通量；S为盐度；FS->s为盐度的水平湍扩散项。
(4)湍方程如式(8)和式(9)所示：
Figure BDA0002585827680000142
Figure BDA0002585827680000143
式中，q2为湍动能的2倍；Kq为垂向湍流混合系数；
Figure BDA0002585827680000144
为绝热递减率修正后的垂向密度梯度；
Figure BDA0002585827680000145
l为湍流混合长度；Fq->s为湍动能的水平湍扩散项；B1、E1、E3均为经验参数；W～为面壁近似函数；Fl->s为湍混合长的水平湍扩散项；垂向动量方程暂时仅保留了全导数项和非静力项。
(5)状态方程如式(10)所示：
ρ＝ρ(T,S,P) (10)
式中，ρ为海水密度；P为压强。
步骤2，根据步骤1构建的基于垂向混合坐标的非静力海洋模式控制方程组建立基于垂向混合坐标的非静力海洋模式，即，将基于垂向混合坐标的非静力海洋模式控制方程组离散化处理再用计算机语言编程。并且，获取内波相关数据，对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理，其中，所述的内波相关数据包括研究区域的地形、温度、盐度、密度、风场、降水和正压潮。
首先，确定所要研究的区域，对该研究区域的地形进行模式网格化处理。其次，查找该研究区域的温度、盐度、密度分布资料，设置海洋模式的初始温度场、初始盐度场和初始密度场；查找该研究区域的风场、降水资料，设置海洋模式的强迫场；根据实际情况设置海洋模式的边界场；查找该研究区域的正压潮资料，设置海洋模式的初始速度场；根据内波模拟要求和实际情况，设置海洋模式的基本参数，所述海洋模式的基本参数包括底粗糙度、时间步长、积分步数、科氏参数和湍混合参数。
步骤3，采用步骤2所建立的基于垂向混合坐标的非静力海洋模式模拟内波，获得内波数据，其中，所述的内波数据包括各方向流速、海温、密度、动压和水位、湍动能、湍流混合长度、垂向湍粘性系数和垂向湍扩散系数。
基于垂向混合坐标的非静力海洋模式利用分列算子将基于垂向混合坐标的非静力海洋模式控制方程组中不同的物理过程分开计算：
步骤3-1，显格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的平流项、科氏力项、斜压项和水平扩散项，如式(11)和式(12)所示：
Figure BDA0002585827680000151
Figure BDA0002585827680000152
式中，(usk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(usk)的中间状态；(vsk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(vsk)的中间状态；n-1为第n-1时刻；Δt为时间步长；n为第n时刻；(usk)n-1为积分第n-1时刻的(usk)的值；
Figure BDA0002585827680000153
为积分第n-1时刻u的水平粘性项；(vsk)n-1为积分第n-1时刻的(vsk)的值；
Figure BDA0002585827680000154
为积分第n-1时刻v的水平粘性项。
步骤3-2，隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项，如式(13)和式(14)所示：
Figure BDA0002585827680000161
Figure BDA0002585827680000162
式中，
Figure BDA0002585827680000163
为垂向扩散项作用后的u的中间状态；
Figure BDA0002585827680000164
为垂向扩散项作用后的v的中间状态；
Figure BDA0002585827680000165
为垂向扩散项作用后的sk的中间状态。
步骤3-3，隐格式求解水位：
首先，给出水位求解方程如式(15)和式(16)所示：
Figure BDA0002585827680000166
Figure BDA0002585827680000167
式中，
Figure BDA0002585827680000168
为水位压强梯度力作用后的u的中间状态；
Figure BDA0002585827680000169
为水位压强梯度力作用后的v的中间状态；
Figure BDA00025858276800001610
为水位压强梯度力作用后的sk的中间状态；
其次，先对式(15)和式(16)进行差分再垂向积分，可得到两个垂向积分式子；再将这两个垂向积分式子带入到式(2)垂向积分形式的差分格式中；整理该差分格式后化为五点隐格式，且关于水位的系数矩阵是严格对角占优的，用超松弛迭代方法得到水位的收敛解。其中，所谓的严格对角占优是指矩阵中每个主对角元素的模都大于与它同行的其他元素的模的总和；对列同样成立。
步骤3-4，求解垂向混合坐标系中水位压强梯度力作用后的垂向速度
Figure BDA00025858276800001611
和z坐标系中水位压强梯度力作用后的垂向速度
Figure BDA00025858276800001612
首先利用式(2)求解垂向混合坐标系中水位压强梯度力作用后的垂向速度
Figure BDA00025858276800001613
然后根据式子w＝ω+(sx+ηx)u+(sy+ηy)v+(st+ηt)进行差分，计算得到z坐标系中水位压强梯度力作用后的垂向速度
Figure BDA00025858276800001614
步骤3-5，隐格式求解动压：
方法一：15对角线求解法
对连续性方程(1)左右同除sk得到变形后的连续性方程如式(17)所示：
Figure BDA00025858276800001615
以动压Q为中心，对变形后的连续性方程(17)进行差分，得到的结果如式(18)至式(22)所示：
Figure BDA0002585827680000171
Figure BDA0002585827680000172
Figure BDA0002585827680000173
Figure BDA0002585827680000181
Figure BDA0002585827680000182
式中，dx为x方向的网格间距；dy为y方向的网格间距；dz为垂向相邻两个网格面的间距；dzz为垂向相邻两个网格中心的间距；
Figure BDA0002585827680000183
Figure BDA0002585827680000184
i为x方向的网格编号，j为y方向的网格编号；
对式(18)至式(22)进行整理：式(18)+式(19)-式(20)+式(21)-式(22)，得：
Figure BDA0002585827680000185
其中，
Figure BDA0002585827680000191
Figure BDA0002585827680000192
Figure BDA0002585827680000193
Figure BDA0002585827680000194
Figure BDA0002585827680000195
Figure BDA0002585827680000196
Figure BDA0002585827680000197
Figure BDA0002585827680000198
Figure BDA0002585827680000199
Figure BDA00025858276800001910
Figure BDA0002585827680000201
Figure BDA0002585827680000202
Figure BDA0002585827680000203
Figure BDA0002585827680000204
Figure BDA0002585827680000205
Figure BDA0002585827680000206
式(23)为15点隐式方程，该15点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(23)时，可得到动压的收敛解。
方法二：7对角线求解法
Figure BDA0002585827680000208
Figure BDA0002585827680000209
对于w方程，忽略掉平流项：
Figure BDA00025858276800002010
将式(40)、式(41)和式(42)进行差分，然后带入式(1)半隐格式的差分形式中，经整理得：
Figure BDA0002585827680000211
其中，
Figure BDA0002585827680000212
Figure BDA0002585827680000213
Figure BDA0002585827680000214
Figure BDA0002585827680000215
Figure BDA0002585827680000216
Figure BDA0002585827680000217
Figure BDA0002585827680000218
Figure BDA0002585827680000219
式中，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 (7)

1.一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，包括以下步骤：
步骤1，构建基于垂向混合坐标的非静力海洋模式控制方程组；
步骤2，根据步骤1构建的基于垂向混合坐标的非静力海洋模式控制方程组建立基于垂向混合坐标的非静力海洋模式；并且，获取内波相关数据，对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理，其中，所述的内波相关数据包括研究区域的地形、温度、盐度、密度、风场、降水和正压潮；
步骤3，采用步骤2所建立的基于垂向混合坐标的非静力海洋模式模拟内波，获得内波数据，其中，所述的内波数据包括各方向流速、海温、密度、动压和水位、湍动能、湍流混合长度、垂向湍粘性系数和垂向湍扩散系数；
步骤4，对步骤3获得的内波数据进行分析，用于对海洋内波进行预报。
2.根据权利要求1所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤1进一步包括：
基于任意因变量在z坐标系和垂向混合坐标系间的转换关系，得到基于垂向混合坐标的非静力海洋模式控制方程组，所述基于垂向混合坐标的非静力海洋模式控制方程组包括连续性方程、动量方程、温盐方程、湍方程和状态方程：
(1)连续性方程有两种形式，如式(1)和式(2)所示：
Figure FDA0002585827670000011
Figure FDA0002585827670000012
式中，u为x方向的速度分量；v为y方向的速度分量；w为z坐标系中z方向的速度分量；ω为垂向混合坐标系中垂向速度分量，ω＝w-(sx+ηx)u-(sy+ηy)v-(st+ηt)；
Figure FDA0002585827670000013
Figure FDA0002585827670000014
η为水位；s为垂向混合坐标系中的垂向坐标；k为垂向网格的编号；x、y为水平坐标，z为z坐标系中的垂向坐标，t为时间；
(2)动量方程包括水平动量方程和垂向动量方程，所述水平动量方程如式(3)和式(4)所示，所述垂向动量方程如式(5)所示：
Figure FDA0002585827670000021
Figure FDA0002585827670000022
Figure FDA0002585827670000023
式中，f为科氏参数；g为重力加速度；ρ0为参考密度；ρ′为海水扰动密度；dk′为积分微元；Q为动压；KM为垂向湍粘性系数；Fx->s为u的水平粘性项；Fy->s为v的水平粘性项；
(3)温盐方程如式(6)和式(7)所示：
Figure FDA0002585827670000024
Figure FDA0002585827670000025
式中，T为海温；KH为垂向湍扩散系数；FT->s为温度的水平湍扩散项；R为短波辐射通量；S为盐度；FS->s为盐度的水平湍扩散项；
(4)湍方程如式(8)和式(9)所示：
Figure FDA0002585827670000026
Figure FDA0002585827670000031
式中，q2为湍动能的2倍；Kq为垂向湍流混合系数；
Figure FDA0002585827670000032
为绝热递减率修正后的垂向密度梯度；q3为
Figure FDA0002585827670000033
l为湍流混合长度；Fq->s为湍动能的水平湍扩散项；B1、E1、E3均为经验参数；
Figure FDA0002585827670000034
为面壁近似函数；Fl->s为湍混合长的水平湍扩散项；
(5)状态方程如式(10)所示：
ρ＝ρ(T,S,P) (10)
式中，ρ为海水密度；P为压强。
3.根据权利要求1所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤2中，所述的对所建立的基于垂向混合坐标的非静力海洋模式进行初始化处理包括：
确定所要研究的区域，对该研究区域的地形进行模式网格化处理；
查找该研究区域的温度、盐度、密度分布资料，设置海洋模式的初始温度场、初始盐度场和初始密度场；查找该研究区域的风场、降水资料，设置海洋模式的强迫场；根据实际情况设置海洋模式的边界场；查找该研究区域的正压潮资料，设置海洋模式的初始速度场；根据内波模拟要求和实际情况，设置海洋模式的基本参数，所述海洋模式的基本参数包括底粗糙度、时间步长、积分步数、科氏参数和湍混合参数。
4.根据权利要求2所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤3进一步包括：
基于垂向混合坐标的非静力海洋模式利用分列算子将基于垂向混合坐标的非静力海洋模式控制方程组中不同的物理过程分开计算：
步骤3-1，显格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的平流项、科氏力项、斜压项和水平扩散项；
步骤3-2，隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项；
步骤3-3，隐格式求解水位；
步骤3-4，求解垂向混合坐标系中水位压强梯度力作用后的垂向速度
Figure FDA0002585827670000042
和z坐标系中水位压强梯度力作用后的垂向速度
Figure FDA0002585827670000043
步骤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直至满足内波模拟所需要的总的时间积分步数，输出内波数据。
5.根据权利要求4所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤3-1中，所述的显格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的平流项、科氏力项、斜压项和水平粘性项如式(11)和式(12)所示：
Figure FDA0002585827670000041
Figure FDA0002585827670000051
式中，(usk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(usk)的中间状态；(vsk)*为水平动量方程中平流项、科氏力项、斜压项和水平粘性项共同作用产生的(vsk)的中间状态；n-1为第n-1时刻；Δt为时间步长；n为第n时刻；(usk)n-1为积分第n-1时刻的(usk)的值；
Figure FDA0002585827670000052
为积分第n-1时刻u的水平粘性项；(vsk)n-1为积分第n-1时刻的(vsk)的值；
Figure FDA0002585827670000053
为积分第n-1时刻v的水平粘性项；
步骤3-2中，所述的隐格式求解基于垂向混合坐标的非静力海洋模式控制方程组中的水平动量方程的垂向扩散项如式(13)和式(14)所示：
Figure FDA0002585827670000054
Figure FDA0002585827670000055
式中，
Figure FDA0002585827670000056
为垂向扩散项作用后的u的中间状态；
Figure FDA0002585827670000057
为垂向扩散项作用后的v的中间状态；
Figure FDA0002585827670000058
为垂向扩散项作用后的sk的中间状态；
步骤3-3中，所述的隐格式求解水位包括：
首先，给出水位求解方程如式(15)和式(16)所示：
Figure FDA0002585827670000059
Figure FDA00025858276700000510
式中，
Figure FDA00025858276700000511
为水位压强梯度力作用后的u的中间状态；
Figure FDA00025858276700000512
为水位压强梯度力作用后的v的中间状态；
Figure FDA00025858276700000513
为水位压强梯度力作用后的sk的中间状态；
其次，先对式(15)和式(16)进行差分再垂向积分，可得到两个垂向积分式子；再将这两个垂向积分式子带入到式(2)垂向积分形式的差分格式中；整理该差分格式后化为五点隐格式，且关于水位的系数矩阵是严格对角占优的，用超松弛迭代方法得到水位的收敛解；
步骤3-4中，所述的求解垂向混合坐标系中水位压强梯度力作用后的垂向速度
Figure FDA0002585827670000061
和z坐标系中水位压强梯度力作用后的垂向速度
Figure FDA0002585827670000062
包括：
首先利用式(2)求解垂向混合坐标系中水位压强梯度力作用后的垂向速度
Figure FDA0002585827670000063
然后根据式子w＝ω+(sx+ηx)u+(sy+ηy)v+(st+ηt)进行差分，计算得到z坐标系中水位压强梯度力作用后的垂向速度
Figure FDA0002585827670000064
步骤3-5中，所述的隐格式求解动压可采用15对角线求解法或7对角线求解法求解动压；
步骤3-7中，所述的求解基于垂向混合坐标的非静力海洋模式控制方程组中的湍方程式(8)和式(9)为：湍方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解；
步骤3-8中，所述的求解基于垂向混合坐标的非静力海洋模式控制方程组中的温盐方程式(6)和式(7)为：温盐方程中的垂向扩散项采用隐格式求解，其它均采用显格式求解。
6.根据权利要求5所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤3-5中，所述的15对角线求解法求解动压包括：
对连续性方程(1)左右同除sk得到变形后的连续性方程如式(17)所示：
Figure FDA0002585827670000065
以动压Q为中心，对变形后的连续性方程(17)进行差分，得到的结果如式(18)至式(22)所示：
Figure FDA0002585827670000066
Figure FDA0002585827670000071
Figure FDA0002585827670000072
Figure FDA0002585827670000073
Figure FDA0002585827670000081
Figure FDA0002585827670000082
式中，dx为x方向的网格间距；dy为y方向的网格间距；dz为垂向相邻两个网格面的间距；dzz为垂向相邻两个网格中心的间距；
Figure FDA0002585827670000083
Figure FDA0002585827670000084
i为x方向的网格编号，j为y方向的网格编号；
对式(18)至式(22)进行整理：式(18)+式(19)-式(20)+式(21)-式(22)，得：
Figure FDA0002585827670000085
其中，
Figure FDA0002585827670000086
Figure FDA0002585827670000087
Figure FDA0002585827670000091
Figure FDA0002585827670000092
Figure FDA0002585827670000093
Figure FDA0002585827670000094
Figure FDA0002585827670000095
Figure FDA0002585827670000096
Figure FDA0002585827670000097
Figure FDA0002585827670000098
Figure FDA0002585827670000099
Figure FDA00025858276700000910
Figure FDA0002585827670000101
Figure FDA0002585827670000102
Figure FDA0002585827670000103
Figure FDA0002585827670000104
式(23)为15点隐式方程，该15点隐式方程的系数矩阵是严格对角占优的，用Gauss-Seidel迭代方法求解式(23)时，得到动压的收敛解。
7.根据权利要求5所述的一种基于垂向混合坐标的海洋内波模拟方法，其特征在于，步骤3-5中，所述的7对角线求解法求解动压包括：
Figure FDA0002585827670000105
Figure FDA0002585827670000106
对于w方程，忽略掉平流项：
Figure FDA0002585827670000107
将式(40)、式(41)和式(42)进行差分，然后带入式(1)半隐格式的差分形式中，经整理得：
Figure FDA0002585827670000111
其中，
Figure FDA0002585827670000112
Figure FDA0002585827670000113
Figure FDA0002585827670000114
Figure FDA0002585827670000115
Figure FDA0002585827670000116
Figure FDA0002585827670000117
Figure FDA0002585827670000118
Figure FDA0002585827670000119
式中，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 true CN111859748A (zh) 2020-10-30
CN111859748B 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)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112612027A (zh) * 2020-12-15 2021-04-06 中国科学院声学研究所 一种浅海环境下利用声能量起伏的海洋内波监测方法
CN113468679A (zh) * 2021-09-06 2021-10-01 中国空气动力研究与发展中心计算空气动力研究所 一种基于s-a模型的湍流长度尺度计算方法
CN114814276A (zh) * 2022-03-21 2022-07-29 汕头大学 海上风电设备运行导致周边海水垂向运动流速的计算方法
CN115391069A (zh) * 2022-10-27 2022-11-25 山东省计算中心（国家超级计算济南中心） 基于海洋模式roms的并行通讯方法及系统

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模式不同垂向坐标设置对比研究", 厦门大学学报. 自然科学版, vol. 53, no. 6 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112612027A (zh) * 2020-12-15 2021-04-06 中国科学院声学研究所 一种浅海环境下利用声能量起伏的海洋内波监测方法
CN113468679A (zh) * 2021-09-06 2021-10-01 中国空气动力研究与发展中心计算空气动力研究所 一种基于s-a模型的湍流长度尺度计算方法
CN114814276A (zh) * 2022-03-21 2022-07-29 汕头大学 海上风电设备运行导致周边海水垂向运动流速的计算方法
CN114814276B (zh) * 2022-03-21 2023-08-18 汕头大学 海上风电设备运行导致周边海水垂向运动流速的计算方法
CN115391069A (zh) * 2022-10-27 2022-11-25 山东省计算中心（国家超级计算济南中心） 基于海洋模式roms的并行通讯方法及系统
CN115391069B (zh) * 2022-10-27 2023-02-03 山东省计算中心（国家超级计算济南中心） 基于海洋模式roms的并行通讯方法及系统

Also Published As

Publication number Publication date
CN111859748B (zh) 2024-04-02

Similar Documents

Publication Publication Date Title
CN111859748A (zh) 一种基于垂向混合坐标的海洋内波模拟方法
Kunz et al. Simulation of the wind field in Athens using refined boundary conditions
McGillicuddy Jr et al. Coupled physical and biological modelling of the spring bloom in the North Atlantic (II): three dimensional bloom and post-bloom processes
Hecht et al. Ocean modeling in an eddying regime
Kordzadze et al. Operational forecast of hydrophysical fields in the Georgian Black Sea coastal zone within the ECOOP
Özgökmen et al. Dynamics of two-dimensional turbulent bottom gravity currents
Kanarska et al. Modelling of seasonal exchange flows through the Dardanelles Strait
Pidcock et al. The spatial variability of vertical velocity in an Iceland basin eddy dipole
Józsa On the internal boundary layer related wind stress curl and its role in generating shallow lake circulations
Michallet et al. Physical modeling of three‐dimensional intermediate beach morphodynamics
Jiang et al. Effects of topographic steering and ambient stratification on overflows on continental slopes: A model study
Song et al. An embedded bottom boundary layer formulation for z-coordinate ocean models
Bryan Potential vorticity in models of the ocean circulation
Liu et al. Bulk, spectral and deep water approximations for Stokes drift: Implications for coupled ocean circulation and surface wave models
Ferziger et al. Numerical simulation of geophysical turbulence
Higdon Numerical modelling of ocean circulation
Väli et al. Simulation of nutrient transport from different depths during an upwelling event in the Gulf of Finland
Dalu A parameterization of heat convection for a numerical sea breeze model
Allahdadi et al. Effect of stratification on current hydrodynamics over Louisiana shelf during Hurricane Katrina
Torres et al. Stratified rotating flow over complex terrain
Cervantes et al. A modeling study of Eulerian and Lagrangian aspects of shelf circulation off Duck, North Carolina
Stern et al. Modeling ice shelf cavities and tabular icebergs using Lagrangian elements
Choboter et al. On the baroclinic instability of axisymmetric rotating gravity currents with bottom slope
Styles et al. The sensitivity of an idealized Weddell Gyre to horizontal resolution
Roy et al. Large eddy simulation of wind flow over complex terrain: The Bolund Hill case

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