跳转至

最小风应力海洋惯性流

背景

风通过海面向海洋输入水平动量。一次风事件过后，近似均匀的上层海水会在科里奥利力作用下， 以接近当地惯性频率的速度旋转。研究者通常先用单层模型，也称平板海洋模型，理解这一响应， 再进一步加入垂向结构和湍流混合。

本例研究一个逆向强迫问题：位于北纬 \(45^\circ\)、深 50 m 的混合层从静止开始，怎样选择空间 均匀的二维海面应力，才能在 24 小时后形成 \(0.20\ \mathrm{m/s}\) 的正东向流，同时使应力平方 积分最小？把完整的应力时间历程作为决策变量后，可以得到一个紧凑且有闭式解的最优控制 基准问题。

这里的最优应力是理想化的科学强迫实验，并不意味着真实风场能够按需操纵或被精确预测。

平板海洋模型

令 \(u\) 和 \(v\) 分别为混合层的东向、北向流速，\(\tau_x\) 和 \(\tau_y\) 为对应的海面应力分量。 时间以小时计，带阻尼的 \(f\) 平面动量方程为

\[ \begin{aligned} \dot u &= f v-r u+b\tau_x,\\ \dot v &=-f u-r v+b\tau_y, \end{aligned} \qquad b=\frac{3600}{\rho H}. \]

其中，\(f=2\Omega\sin\phi\) 从 \(\mathrm{s^{-1}}\) 换算为 \(\mathrm{h^{-1}}\)，\(r\) 是线性 Rayleigh 动量阻尼率，\(\rho\) 是海水密度，\(H\) 是混合层深度。符号约定对应北半球的东-北 坐标：无强迫水流会顺时针旋转。

数值参数为

\[ \begin{gathered} \phi=45^\circ,\qquad \Omega=7.2921159\times10^{-5}\ \mathrm{rad/s},\\ f=0.3712539314\ \mathrm{h^{-1}},\qquad \frac{2\pi}{f}=16.9242\ \mathrm h,\\ \rho=1025\ \mathrm{kg/m^3},\qquad H=50\ \mathrm m,\qquad r=(48\ \mathrm h)^{-1},\\ b=0.07024390244\, \frac{\mathrm{m/s}}{\mathrm h\,(\mathrm{N/m^2})}. \end{gathered} \]

目标函数、端点与路径边界

固定时域为 \(T=24\ \mathrm h\)。归一化的平均强迫代价为

\[ \min_{\tau_x,\tau_y} J =\frac{1}{T}\int_0^T \left[ \left(\frac{\tau_x}{\tau_s}\right)^2 +\left(\frac{\tau_y}{\tau_s}\right)^2 \right]\,\mathrm dt, \qquad \tau_s=0.25\ \mathrm{N/m^2}. \]

尺度 \(\tau_s\) 使 \(J\) 成为无量纲量，但不改变最优应力历程。端点条件为

\[ \begin{bmatrix}u(0)\\v(0)\end{bmatrix} =\begin{bmatrix}0\\0\end{bmatrix}\ \mathrm{m/s}, \qquad \begin{bmatrix}u(T)\\v(T)\end{bmatrix} =\begin{bmatrix}0.20\\0\end{bmatrix}\ \mathrm{m/s}. \]

逐分量路径边界把配点问题限制在明确的物理范围内：

\[ |u(t)|,|v(t)|\le0.35\ \mathrm{m/s}, \qquad |\tau_x(t)|,|\tau_y(t)|\le0.30\ \mathrm{N/m^2}. \]

这些是分量限制，而不是欧氏模长限制。无约束精确最优解严格位于全部边界以内，因此它同时也是 这个约束问题的最优解。

可控性 Gramian 精确解

记 \(\mathbf{x}=[u,v]^\mathsf{T}\)， \(\boldsymbol{\tau}=[\tau_x,\tau_y]^\mathsf{T}\)，则

\[ \dot{\mathbf{x}}=A\mathbf{x}+B\boldsymbol{\tau}, \qquad A=\begin{bmatrix}-r&f\\-f&-r\end{bmatrix}, \qquad B=bI_2. \]

状态转移矩阵为

\[ \Phi(t)=e^{At} =e^{-rt} \begin{bmatrix} \cos(ft)&\sin(ft)\\ -\sin(ft)&\cos(ft) \end{bmatrix}. \]

由于阻尼和强迫都是各向同性的，有限时域可控性 Gramian 可以化为一个标量乘单位矩阵：

\[ W(T)=\int_0^T \Phi(T-s)BB^\mathsf{T}\Phi(T-s)^\mathsf{T}\,\mathrm ds =wI_2, \qquad w=b^2\frac{1-e^{-2rT}}{2r}. \]

当 \(\mathbf{x}(0)=0\) 且终点 \(\mathbf{x}_T\) 固定时，唯一的最小能量应力为

\[ \boldsymbol{\tau}^*(t) =B^\mathsf{T}\Phi(T-t)^\mathsf{T}W(T)^{-1}\mathbf{x}_T, \]

其未归一化能量为

\[ E^*=\int_0^T \boldsymbol{\tau}^{*\mathsf{T}}\boldsymbol{\tau}^*\,\mathrm dt =\mathbf{x}_T^\mathsf{T}W(T)^{-1}\mathbf{x}_T. \]

对应流速也有闭式表达式：

\[ \mathbf{x}^*(t)= \frac{b^2e^{-r(T-t)}(1-e^{-2rt})}{2rw} \begin{bmatrix} \cos(f(T-t))&-\sin(f(T-t))\\ \sin(f(T-t))&\cos(f(T-t)) \end{bmatrix} \mathbf{x}_T. \]

这组解析解既提供可行初值，更重要的是为数值最优解提供独立于配点离散的参考。

变量与单位

符号 含义 单位
\(t,T\) 时间和强迫时域 h
\(u,v\) 混合层东向、北向流速 m/s
\(\tau_x,\tau_y\) 东向、北向海面应力 N/m²
\(\tau_s\) 目标函数归一化尺度 N/m²
\(f\) 小时时间坐标下的科里奥利参数 h⁻¹
\(r\) 线性动量阻尼率 h⁻¹
\(\rho\) 海水密度 kg/m³
\(H\) 混合层深度 m
\(b\) 应力到流速变化率的增益 \((\mathrm{m/s})/[\mathrm h(\mathrm{N/m^2})]\)

离散与独立检查

默认求解使用 48 个均匀 Lobatto 网格区间，每个区间含四个点。脚本在 4,001 个时刻重构 连续解，并在非线性规划收敛判断之外执行四类检查：

  1. 逐点比较数值流速、应力与可控性 Gramian 精确解。
  2. 对优化应力作线性插值，再用独立的高精度 DOP853 求解器重新积分两条物理动量方程。
  3. 比较配点目标值、密集梯形积分应力能量和精确最小能量。
  4. 在密集网格上检查端点残差与全部分量边界。

--quick 模式使用 32 个区间和相同的四点多项式阶次。它仍会与解析解比较，而不只是检查 求解器是否报告收敛。

验证结果

指标 默认网格结果
归一化目标函数 \(J\) \(0.3562380575\)
精确归一化目标函数 \(0.3562380568\)
应力能量 \(E\) \(0.5343570891\ (\mathrm{N/m^2})^2\,\mathrm h\)
流速模长峰值 \(0.20000000\ \mathrm{m/s}\)
应力模长峰值 \(0.18767921\ \mathrm{N/m^2}\)
相对精确解的最大流速误差 \(1.279\times10^{-7}\ \mathrm{m/s}\)
相对精确解的最大应力误差 \(9.652\times10^{-6}\ \mathrm{N/m^2}\)
独立前向积分最大误差 \(1.705\times10^{-7}\ \mathrm{m/s}\)
密集路径边界最大违反量 \(0\)

流速轨迹图直观显示了科里奥利旋转。最优应力本身也在旋转，其包络向终端时刻增大，因为较晚 施加的强迫所产生的流速响应在到达目标前受到的阻尼衰减更小。

最优海洋混合层流速轨迹、风应力分量与累计强迫代价

建模范围与局限

平板模型假定水层水平一致、垂向均匀，并采用恒定深度、密度、科里奥利参数和线性阻尼。 它省略压强梯度强迫、非线性平流、层结、卷吸、湍流垂向结构、海岸线、地形、波浪和海气 反馈。Rayleigh 阻尼只是耗散的降阶表示，不是湍流闭合模型。

计算还假定优化器能够使用完整且确定的应力时间历程。真实海洋预报必须处理大气强迫不确定性、 数据同化，通常还需要反馈或滚动时域估计。本例适合验证最优控制离散并理解旋转线性动力学， 不能作为业务化风场或海洋预报。

运行示例

运行经过验证的默认网格：

python -m examples.ocean_inertial_current

运行较快的解析烟雾测试：

python -m examples.ocean_inertial_current --quick --no-show

保存图片且不打开窗口：

python -m examples.ocean_inertial_current --save ocean-inertial-current.png --no-show

源代码

完整可运行示例见： examples/ocean_inertial_current.py。