• 返回主页
  • shallowWaterFoam解析

    Reviewed by:
    求解器类型:不可压缩浅水方程求解器。

    “浅水方程”是一系列的偏微分方程组。用来描述计算域的高度远小于宽度的不可压缩流体流动。近岸前的海啸问题可以很好的用浅水方程来模拟。

    在下文中我们用“水”来代表需要模拟的流体。考虑一个二维、无旋的浅水域。水的密度,$\rho$,为常数。如图(1)所示$^1$:
    我们假定水是无粘的。浅水方程中的一个重要假定为:水是非常“浅”的,以至于在垂直的方向上流体的速度是相同的。进一步的,在水表处的压力即为大气的压力,$p_0$,其为一个常数且各处均一。

    在上图中,我们考虑$A$和$B$内的封闭区域。这个区域的宽为$\mathrm{d}x$,高为$\mathrm{d}h$。在界面$A$处的速度为$u$,界面$B$处的速度为$u+\mathrm{d}u$。那么在这个封闭区域内的水的质量则为:$\mathrm{d}m=\rho h\mathrm{d}x$。 区域内的质量并不能任意的产生或者消失。因此进入封闭区域的质量和离开封闭区域的质量应该是相等的。在单位时间内,进入$A$的水的质量为$\rho u(x) h(x)$,离开$B$的水的质量为$\rho u(x+\mathrm{d}x) h(x+\mathrm{d}x)$。那么,封闭区域内积累(消失)的质量可以表示为: \begin{equation} \frac{\partial 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$,将其代入到公式(1)有: \begin{equation} \frac{\partial h}{\partial t}=-\frac{\partial uh}{\partial x} \end{equation} 也即: \begin{equation} \frac{\partial h}{\partial t}+\nabla \cdot (uh)=0 \end{equation} 其中$h$为水面高度。公式(3)即为浅水方程中的水深方程。

    回到图(1),水内的压力可以表示为高度的函数: \begin{equation} p(x,z)=p_0+\rho g \left(h\left(x \right)-z\right) \end{equation} 其中$g$为重力加速度。依据牛顿第二定律:$m\frac{\mathrm{D}u}{\mathrm{D}t}=F$,其中$F$为作用在控制体上的力。因为我们假定是无粘的,因此不存在剪切力。故$F$只包含了压力。我们分析图(4)中的控制体:
    在这个控制体上,界面$A$所受到的压力为$F_1$,也即$F_1=\int_0^hp(x,z)\mathrm{d}z$。同理,界面$B$所受到的压力为$F_2=\int_0^hp(x+\mathrm{d}x,z)\mathrm{d}z$。下面我们来看$F_s$,其为一个作用在水表的向内的压力,其大小为$p_0$。首先,作用在弧$\mathrm{d}l$上的压力为:$F_l=p_0 \mathrm{d}l$。图(2)中$F_l$在$x$方向上的分量为$F_s$,且$F_s=p_0\mathrm{d}l \mathrm{sin} \alpha$。将$\mathrm{d}l=\frac{\mathrm{d}x}{\mathrm{cos} \alpha}$代入有:$F_s=p_0 \mathrm{d}x \mathrm{tan} \alpha=p_0 \frac{\partial h}{\partial x}\mathrm{d}x$。将所有的力加和有: \begin{equation} F=p_0 \frac{\partial h}{\partial x}\mathrm{d}x +\int_0^hp(x,z)\mathrm{d}z -\int_0^hp(x+\mathrm{d}x,z)\mathrm{d}z \end{equation} 分析方程(5)中右边的第二项,将公式(4)代入有: \begin{equation} \int_0^hp(x,z)\mathrm{d}z=\int_0^h p_0 \mathrm{d}z+\rho g \int_0^h \left(h\left(x \right)-z\right) \mathrm{d}z=p_0 h(x)+\frac{1}{2}\rho g h(x)^2 \end{equation} 参考公式(6),公式(5)右边的后面俩项可以化为: \begin{equation} \int_0^hp(x,z)\mathrm{d}z -\int_0^hp(x+\mathrm{d}x,z)\mathrm{d}z =p_0 h(x)+\frac{1}{2}\rho g h(x)^2-p_0 h(x+\mathrm{d}x)-\frac{1}{2}\rho g h(x+\mathrm{d}x)^2 \\ =-p_0\frac{\partial h}{\partial x}\mathrm{d}x - \rho g h \frac{\partial h}{\partial x} \mathrm{d}x \end{equation} 将公式(7)代入到公式(5)有: \begin{equation} F=- \rho g h \frac{\partial h}{\partial x} \mathrm{d}x \end{equation} 依据牛顿第二定律,也即: \begin{equation} m\frac{\mathrm{D}u}{\mathrm{D}t}=- \rho g h \frac{\partial h}{\partial x} \mathrm{d}x \end{equation} 将$m=\rho h \mathrm{d}x$代入到公式(9)有: \begin{equation} \frac{\mathrm{D}u}{\mathrm{D}t}=- g \frac{\partial h}{\partial x} \end{equation} 也即: \begin{equation} \frac{\partial u}{\partial t}+u \cdot \nabla u=- g \nabla h \end{equation} 目前公式(11)中我们只考虑了压力。对于有旋的情况,还需要添加科氏力源项。考虑科氏力源项,公式(11)演变为: \begin{equation} \frac{\partial u}{\partial t}+u \cdot \nabla u=- g \nabla h - f \times u \end{equation} 结合连续性方程,最终我们有: \begin{equation} \frac{\partial u}{\partial t}+\nabla \cdot (u u)=- g \nabla h - f \times u \end{equation} 其中: \begin{equation} f=2 \Omega \mathrm{sin} \phi \end{equation} 其中$\phi$为纬度。需要注意的是,在OpenFOAM中,$f=2\Omega \cdot \frac{g}{|g|} \frac{g}{|g|}$。目前并不清楚OpenFOAM中的定义从何而来。公式(3)以及公式(13)即为浅水方程的控制方程。
    下面我们进入代码分析,首先进入shallowWaterFoam.C文件,其首先通过readGravitationalAcceleration.H$^2$创建重力加速度。然后进入createFields.C

    Info<< "Reading field h\n" << endl;
    volScalarField h//水深(水高)
    (
        IOobject
        (
            "h",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
    
    Info<< "Reading field h0 if present\n" << endl;
    volScalarField h0//基准高度,例如床高
    (
        IOobject
        (
            "h0",
            runTime.findInstance("polyMesh", "points"),
            mesh,
            IOobject::READ_IF_PRESENT
        ),
        mesh,
        dimensionedScalar("h0", dimLength, 0.0)
    );
    
    Info<< "Creating field hU\n" << endl;
    volVectorField hU
    (
        IOobject
        (
            "hU",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
    
    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        hU/h
    );
    
    Info<< "Creating field hTotal for post processing\n" << endl;
    volScalarField hTotal//水深+床高,用于后处理
    (
        IOobject
        (
            "hTotal",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        h+h0
    );
    hTotal.write();
    
    #include "createPhi.H"//$(hu)_f \cdot S_f$
    
    Info<< "Creating Coriolis Force" << endl;
    const dimensionedVector F("F", ((2.0*Omega) & gHat)*gHat);//科氏力项,OpenFOAM中的定义
    
    mesh.setFluxRequired(h.name());
    

    回到shallowWaterFoam.C,在进入求解流程之前,OpenFOAM对求解的代码对应于公式(14)左右乘以了$h$,即: \begin{equation} \frac{\partial h u}{\partial t}+\nabla \cdot (h u u)=- g h \nabla h - f \times h u \end{equation} 下面进入求解流程:

        while (runTime.loop())
        {
            Info<< "\n Time = " << runTime.timeName() << nl << endl;
    
            #include "CourantNo.H"
    
            // --- Pressure-velocity PIMPLE corrector loop
            while (pimple.loop())
            {
                surfaceScalarField phiv("phiv", phi/fvc::interpolate(h));//\frac{(hu)_f \cdot S_f}{h_f}
    
                fvVectorMatrix hUEqn
                (
                    fvm::ddt(hU)//公式(15)左边第一项
                  + fvm::div(phiv, hU)//公式(15)左边第二项
                );
    
                hUEqn.relax();
    
                if (pimple.momentumPredictor())
                {
                    if (rotating)
                    {
                        solve(hUEqn + (F ^ hU) == -magg*h*fvc::grad(h + h0));//考虑科氏力
                    }
                    else
                    {
                        solve(hUEqn == -magg*h*fvc::grad(h + h0));//不考虑科氏力
                    }
    
                    // Constrain the momentum to be in the geometry if 3D geometry
                    if (mesh.nGeometricD() == 3)
                    {
                        hU -= (gHat & hU)*gHat;
                        hU.correctBoundaryConditions();
                    }//在3D网格中,将hU中和重力加速度矢量存在值的分量的值归0
                }
    

    上述方程即为通过公式(15)对速度,$hu$,进行动量预测。在有了预测的动量速度之后,需要通过水深方程也即公式(3)来求水深。遵循OpenFOAM组建压力泊松方程的步骤$^3$,对公式(3)重新整理有我们有最终的方程: \begin{equation} \frac{\partial h}{\partial t} + \nabla \cdot (\mathrm{HbyA}) = \nabla \cdot \left(\frac{1}{A_p} gh\nabla h \right) \end{equation} 回到shallowWaterFoam.C,其对应为:

                        fvScalarMatrix hEqn
                        (
                            fvm::ddt(h)
                          + fvc::div(phiHbyA)
                          - fvm::laplacian(ghrAUf, h)
                        );
    

    求解后我们有水深场。
    如何引用:
    李东岳. shallowWaterFoam解析[OL]. http://dyfluid.com/shallowWaterFoam.html, 2016.3.26
    [1]. The Shallow Water Equations[OL] http://eaps.mit.edu/~rap/courses/12333_notes/A2%20SWeqs.pdf
    [2]. 李东岳. twoLiquidMixingFoam解析[OL]. http://dyfluid.com/twoLiquidMixingFoam.html, 2016.3.10
    [3]. 李东岳. icoFoam解析[OL]. http://dyfluid.com/icoFoam.html, 2016.3.20


    东岳流体
    info@dyfluid.com