return 0;
wmake

shallowWaterFoam解析

如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次


引言

浅水方程用来描述计算域的宽度远大于高度的不可压缩流体流动现象。主要用于全球气候问题、台风预测、洋流问题、近岸问题等的数值模拟。OpenFOAM中的潜水方程求解器为shallowWaterFoam,附加考虑地球旋转引起的科氏力效应,后演变为单独的AtmosFoam程序,其为OpenFOAM创始人Henry Weller以及其夫人Hilary Weller在2008年共同开发的求解器(Weller and Weller, 2008)。本文从流动微元开始,推导浅水方程的控制方程,并介绍其在OpenFOAM中的离散求解。

如果打不开图像,请右键在新标签页打开图像后刷新几次
通过shallowWaterFoam模拟的波浪

控制方程

本文用水来代表需要模拟的流体。考虑一个二维、无旋的浅水域。引入如下假定

如果打不开图像,请右键在新标签页打开图像后刷新几次

图1. 虚线内的封闭区域宽为$\mathrm{d}x$。左侧高为$h(x)$,速度为$u(x)$。右侧高为$h(x+dx)$,速度为$u(x+\mathrm{d}x)$。
虚线区域内体积为$h(x)\mathrm{d}x$(横截面积为1),质量为:$\mathrm{d}m=\rho h(x)\mathrm{d}x$。

参考连续性方程的推导思路。连续性原理要求区域内的质量并不能任意的产生或者消失。因此进入封闭区域的质量和离开封闭区域的质量为区域内质量的变化。现考虑$x$方向,单位时间内,进入$A$的水的质量为$\rho u(x) h(x)$,离开$B$的水的质量为$\rho u(x+\mathrm{d}x) h(x+\mathrm{d}x)$。那么,封闭区域内质量的变化率可表示为:

\begin{equation}\label{cont1} \frac{\rd m}{\partial x}=\rho u(x) h(x)-\rho u(x+\mathrm{d}x) h(x+\mathrm{d}x) \\ =-\rho\frac{\partial (uh)}{\partial x} \mathrm{d}x \end{equation}

由于$\mathrm{d}m=\rho h\mathrm{d}x$,将其代入到方程\eqref{cont1}有:

\begin{equation} \frac{\partial h}{\partial t}=-\frac{\partial uh}{\partial x} \end{equation}

同理,$y$方向上有:

\begin{equation} \frac{\partial h}{\partial t}=-\frac{\partial vh}{\partial y} \end{equation}

对$x,y$方向加和有:

\begin{equation}\label{cont2} \frac{\partial h}{\partial t}+\nabla \cdot (\bfU h)=0 \end{equation}

方程\eqref{cont2}即为浅水方程中的水深方程。其物理意义很明确。如果把水深$h$理解为一个标量,那么如果进来较多的标量,必然引起标量的堆积,也即水深的增加。

现推导$x$方向的动量方程。水面下的压力可以表示为高度的函数:

\begin{equation}\label{bnl} p(x,z)=p_0+\rho g \left(h\left(x \right)-z\right) \end{equation}

其中$g$为重力加速度,$z$表示水面下某点的高度。$\rho$为水的密度。依据牛顿第二定律:

\begin{equation}\label{mom} m\frac{\mathrm{D}u}{\mathrm{D}t}=F_x\end{equation}

其中$F_x$为作用在区域上$x$方向的力。因为浅水方程假定是无粘的,因此不存在剪切力。故$F_x$只是压力的分量。现分析下面的控制体:


图2. 虚线内的封闭区域左侧受力为$F(x)$,右侧受力为$F(x+dx)$。
上方曲线$\rd l$受力为$F(f)$,其$x$反向受力为$F(f)_x$。

在这个控制体上,左侧所受到的总压力为:

\begin{equation}\label{x0} F(x)=\int_0^{h(x)}p(x,z)\mathrm{d}z \end{equation}

将方程\eqref{bnl}代入到方程\eqref{x0}有

\begin{equation}\label{x1} F(x)=\int_0^{h(x)}\left(p_0+\rho g \left(h\left(x \right)-z\right)\right)\mathrm{d}z \\ =p_0h(x)+\rho g \int_0^{h(x)}\left(h(x)-z)\right)\rd z \\ =p_0h(x)+\frac{1}{2}\rho g h(x)^2 \end{equation}

同理,右侧所受到的总压力为(注意反向的负号)

\begin{equation}\label{x2} F(x+\rd x)=-p_0h(x+\rd x)-\frac{1}{2}\rho g h(x+\rd x)^2 \end{equation}

现分析$F(f)$,其为一个作用在水表的向内的压力,大小为

\begin{equation} F(f)=p_0\rd l \end{equation}

那么其在$x$方向上的分量为$F(f)_x$为

\begin{equation}\label{x3} F(f)_x=F(f)\sin\alpha=p_0\rd l\frac{\rd h}{\rd l}=p_0\rd h=p_0\frac{\partial h}{\partial x}\rd x \end{equation}

将方程\eqref{x1},\eqref{x2}和\eqref{x3}加和有$x$方向所有的受力:

\begin{equation}\label{x} F_x=p_0\left(h(x)-h(x+\rd x)\right)+\frac{1}{2}\rho g \left(h(x)^2-h(x+\rd x)^2\right)+p_0\frac{\partial h}{\partial x}\rd x \\ =-p_0\frac{\partial h}{\partial x}\rd x-\rho g h\frac{\partial h}{\partial x}\rd x+p_0\frac{\partial h}{\partial x}\rd x \\ =-\rho g h\frac{\partial h}{\partial x}\rd x \end{equation}

将方程\eqref{x}代入到方程\eqref{mom}有

\begin{equation}\label{} m\frac{\mathrm{D}u}{\mathrm{D}t}=-\rho g h\frac{\partial h}{\partial x}\rd x \end{equation}

\begin{equation}\label{momx} h\frac{\mathrm{D}u}{\mathrm{D}t}=- gh \frac{\partial h}{\partial x} \end{equation}

同理,$y$方向有

\begin{equation}\label{momy} h\frac{\mathrm{D}v}{\mathrm{D}t}=- gh \frac{\partial h}{\partial y} \end{equation}

对方程\eqref{momx},\eqref{momy}加和有

\begin{equation}\label{momm} \frac{\partial h\bfU}{\partial t}+h\bfU \cdot \nabla \bfU=- gh \nabla h \end{equation}

目前方程\eqref{momm}中只考虑了压力。对于毛细波,可能还需要考虑表面张力。对于有旋的情况,还需要添加科氏力。若考虑科氏力源项,方程\eqref{momm}演变为:

\begin{equation}\label{momm2} \frac{\partial h\bfU}{\partial t}+h\bfU \cdot \nabla \bfU=- gh \nabla h - 2h\Omega\sin\phi \times \bfU \end{equation}

其中$\Omega$为地球旋转角速度,$\phi$为纬度。结合连续性方程有最终的动量方程

\begin{equation}\label{momm3} \frac{\partial h\bfU}{\partial t}+\nabla\cdot(h\bfU \bfU)=- gh \nabla h - 2h\Omega\sin\phi \times \bfU \end{equation}

同时,若水底本身存在高度$h_0$,则应在动量方程上附加相应的压力变化:

\begin{equation}\label{momm4} \frac{\partial h\bfU}{\partial t}+\nabla\cdot(h\bfU \bfU)=- gh \nabla (h+h_0) - 2h\Omega\sin\phi \times \bfU \end{equation}

方程离散以及求解

浅水方程的离散和普通N-S方程的离散大同小异。只不过动量方程的未知量变为$h\bfU$,压力方程的变量为$h$(通过连续性方程组建)。首先对方程\eqref{momm3}进行半离散化,求得预测的$h\bfU$场($\mathrm{HbyA}$)。然后通过方程\eqref{cont2}构建水深方程,其最终形式为

\begin{equation}\label{lap} \frac{\partial h}{\partial t}+\nabla\cdot\left(\mathrm{HbyA}-gh\nabla h-gh\nabla h_0\right)=0 \end{equation}

参考文献

Weller, H., Weller, H., 2008. A high-order arbitrarily unstructured finite‐volume model of the global atmosphere: Tests solving the shallow‐water equations. International journal for numerical methods in fluids, 56, 1589-1596.

更新历史
2018.10.29大修

东岳流体®版权所有
勘误、讨论、补充内容请前往CFD中文网

');