计算科学 - 数值求解椭圆 (?) Navier-Stokes 方程的线性化系统（浅水导出） - 吾爱随笔录

在我的博士论文中，我的导师要求我构建一个求解器，其灵感来自“ HF Radar，JL Devonon，1989 年将最优控制理论应用于潮汐流映射的客观分析”一文。在本文中，在数值流体求解部分，作者通过去除对流项并使用线性摩擦系数来线性化浅水二维方程。然后他对方程使用傅里叶变换，最终得到系统：

\begin{aligned} - i σ \vec{U} + \vec{f} \times \vec{U} + g \vec{\nabla} h & = - r \vec{U} \\ i σ h + \vec{\nabla} \cdot (H \vec{U}) & = 0 \end{aligned}

$\begin{aligned} -i\sigma\vec{U}+\vec{f}\times\vec{U}+g\vec{\nabla}h&=-r\vec{U}\\\\ i\sigma h+\vec{\nabla}\cdot\left(H\vec{U}\right)&=0 \end{aligned}$

这里 $i$ 是虚数的基础， $\sigma$ 给定傅里叶分量的频率， $\vec{f}$ 科里奥利力矢量， $g$ 重力加速度， $h$ 当地海平面高度与“平均当地水平”相比， $r$ 线性摩擦系数， $H$ 局部深度和 $\vec{U}$ 局部速度。

然后，将该系统简化为形式为的单向量方程：

F (\vec{V}) = λ \vec{V} + \vec{f} \times \vec{V} + \vec{\nabla} (\vec{\nabla} \cdot (γ \vec{V})) = \vec{0}

$F(\vec{V})=\lambda\vec{V}+\vec{f}\times\vec{V}+\vec{\nabla}\left(\vec{\nabla}\cdot\left(\gamma\vec{V}\right)\right)=\vec{0}$

在哪里：

λ = r H - i σ, γ = \frac{- i g H}{σ}

$\lambda = rH-i\sigma,\quad \gamma=\frac{-igH}{\sigma}$

和 $\vec{V}$ 是传输，作为速度（假设在整个水柱上为常数）和水柱的局部高度（假设波动）的乘积给出 $h$ 与整体深度相比很小）。

然后，他解决了具有指定边界条件的狄利克雷问题。为此，他使用了一种迭代松弛方法（来自 Smith, GD, Numerical solution of partial different equations: Finite Difference methods, 1978）。并且只使用中心离散化方案。

然后我编写了自己的代码来解决它。我的第一次尝试是在具有固定边界条件的矩形盆上。但是，它似乎不适用于与 0 不同的任何条件。我还试图找出分析解决方案，但我发现的唯一（非常明显）解决方案是到处都是空解决方案。这是我使用的代码示例：

############################################
# First test of implementing Devenon's     #
# Equation system and solving them in a clo#
# -sed basin                               #
############################################

#*****Packages import**********************#
import numpy as np



#****Functions*****************************#


def  relaxation_scheme_Jean_Luc_overrelaxed(U,V,dy,dx,nx,ny,gamma,lamda,omega):
 """Function use to solve 2D linearized equation
 of Navier-Stokes for a rigid lid. From Devenon 1989
 . We're trying an overrelaxation technique with omega
  as the relaxation parameter.

 INPUT:
 -----
 U,V: Zonal and meridional speed, size nx*ny, also the output,complex valued
 nx: Number of grid cell in x direction
 ny: Number of grid cell in y direction
 dx: Size of a grid cell in x direction
 dy: Size of a grid cell in y direction
 b: Forcing source function, of size nx x ny
 gamma: Local complex phase speed
 lamda : Local complex wave vector
 omega: Relaxation factor, must belong to [1:2]

 OUTPUT:
 -------
 U,V: Already specified earlier """

 # First we compute the non-homogeneous terms

 #U_b=f*V[1:-1,1:-1]+ (gamma[:-2,:-2]*V[:-2,:-2] + 
 #       gamma[2:,2:]*V[2:,2:]-gamma[2:,:-2]*V[2:,:-2]-
 #       gamma[:-2,2:]*V[:-2,2:])/(4*dy*dx) 

 #V_b= -f*U[1:-1,1:-1]+ (gamma[:-2,:-2]*U[:-2,:-2] +
 #         gamma[2:,2:]*U[2:,2:]- gamma[2:,:-2]*U[2:,:-2]-
 #         gamma[:-2,2:]*U[:-2,2:])/(4*dx*dy)

 U_b=V[1:-1,1:-1]*0.
 V_b=V[1:-1,1:-1]*0.
 # Now we can use relaxation technique with an external term
 U[1:-1,1:-1] = 1/lamda[1:-1,1:-1] * (omega*gamma[2:,1:-1]*U[2:,1:-1]+
                2*(1-omega)*gamma[1:-1,1:-1]*U[1:-1,1:-1]+(omega-1)*gamma[:-2,1:-1]*U[:-2,1:-1]) \
                /dx**2  +1/lamda[1:-1,1:-1]*U_b

 V[1:-1,1:-1] = 1/lamda[1:-1,1:-1] * (omega*gamma[1:-1,2:]*V[1:-1,2:]+
                 2*(1-omega)*(gamma[1:-1,1:-1]*V[1:-1,1:-1])+(omega-1)*gamma[1:-1,:-2]*V[1:-1,:-2]) \
                 /dy**2 +1/lamda[1:-1,1:-1]*V_b

 return U,V

#*****Parameters and variables*************#
nx = 40
ny = 40 # Number of grid points in each direction
nt  = 100 # Number of iteration of the relaxation method
xmin = 0.
xmax = 12000.
ymin = 0.
ymax = 12000. # Dimension (in meters) of the grid

dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)

# Initialization
# I understood that each speed would be complex-valued
# However, I wondered if U and V could only be the real and imaginary
# Part of the same variable in this case
U,V  = np.zeros((nx,ny),dtype=complex),np.zeros((nx,ny),dtype=complex)
x  = np.linspace(xmin, xmax, nx)
y  = np.linspace(xmin, xmax, ny)
X, Y = np.meshgrid(x, y)


H=10 # Setting constant depth for testing
sigma=2.31484e-5 # Approximately semidiurnal frequency,  for tide (in hertz)
g=9.81 # Gravity acceleration
K=2.5e-2 #Friction coefficient
f=2e-4
lamda=np.ones((nx,ny),dtype=complex)*(K-sigma*1j)# (Friction and wavenumber)
gamma=np.ones((nx,ny),dtype=complex)*-1j*g*H/sigma #(Gravity and wavenumber)

#Initial conditions
U[0,:]=1.
U[-1,:]=1.
U[:,0]=1.
U[:,-1]=1.

V[0,:]=1.
V[-1,:]=1.
V[:,0]=1.
V[:,-1]=1.

#********Solving****************************#

for i in range(nt):    
    U,V=relaxation_scheme_Jean_Luc_overrelaxed(U,V,dy,dx,nx,ny,gamma,lamda,2.)

对于任何情况，我都使用了除了常数之外的情况，我发现是非常稳定的行为（系统完全保持在其初始位置）还是非常不稳定的行为（我的不稳定性从传播的两侧开始增长，很快给我异常值）。这种行为似乎是由参数 lambda、gamma 以及单元格的大小驱动的。但是，我尝试过的任何组合都没有帮助我找到解决方案。这让我相信我在编码方面做得非常糟糕。

我尝试搜索该站点（以及一般的网络以寻找答案），并找到了许多不错的 EDP 甚至 Navier-Stokes 求解器的链接。虽然我没有得到任何可以链接到我的问题的东西，可能是由于对它的理解不足。话虽这么说，我还有另一个限制，我的顾问希望这成为一个教育和简单的问题。因此，我必须了解求解器的内部功能。

然后我有几个问题要问你：

1）第一个很明显：为什么它不起作用？我在编码方面做得不好吗？我应该纠正什么以获得正确的工作算法？

2）第二个不太明显：我可以用一种复杂的方式写出我的向量 V，使用传输的一个分量作为实部，第二个作为虚部吗？这对我的顾问来说似乎很明显，但对我来说一点也不明显。

3）该论文将其称为“椭圆”问题，接近亥姆霍兹问题。但是，我只看到两个带有耦合项的 ODE。在这种情况下我可以谈谈椭圆问题吗？