【OpenFOAM】-olaFlow-算例6- waveFloatingObject

算例路径： olaFlow\tutorials\waveFloatingObject
算例描述： 波浪作用下的浮体的刚体运动，属于流固耦合（FSI）问题
学习目标： 动网格设置和使用，网格变形控制，浮体的物理参数设置，浮体做刚体运动的约束设置
算例快照：
在这里插入图片描述

图1 波浪与浮体的相互作用

在这里插入图片描述

图2 波浪与浮体的相互作用过程中的网格变形

文件结构：

.
├── 0.org
│   ├── U
│   ├── alpha.water
│   ├── epsilon
│   ├── k
│   ├── nut
│   ├── p_rgh
│   └── pointDisplacement
├── cleanCase
├── constant
│   ├── dynamicMeshDict
│   ├── dynamicMeshDict.sixDoF
│   ├── dynamicMeshDict.solid
│   ├── g
│   ├── transportProperties
│   ├── turbulenceProperties
│   └── waveDict
├── runCase
└── system
    ├── blockMeshDict
    ├── controlDict
    ├── decomposeParDict
    ├── fvSchemes
    ├── fvSolution
    ├── setFieldsDict
    └── topoSetDict

算例文件解析：

【runCase】

#!/bin/bash
mkdir 0

echo blockMesh meshing...
blockMesh > blockMesh.log

echo Creating floating object...
topoSet > topoSet.log    
subsetMesh -overwrite c0 -patch floatingObject > subsetMesh.log
# 用 toposet 工具选中Box c0
# 用 subsetMesh 工具扣除掉 c0 包含的单元，并将暴露出的内部面形成的 patch 重命名为 floatingObject

echo Preparing 0 folder...
rm -fr 0
cp -r 0.org 0

echo Setting the fields...
setFields > setFields.log

echo Running...
olaDyMFlow > olaDyMFlow.log

echo Simulation complete.

【0.org\U】

dimensions      [0 1 -1 0 0 0 0];
internalField   uniform (0 0 0);
boundaryField
{
    inlet
    {
        type            waveVelocity;
        waveDictName    waveDict;
        value           uniform (0 0 0);
    }
    outlet
    {
        type            waveAbsorption2DVelocity;
        absorptionDir   666.0;
        value           uniform (0 0 0);
    }
    stationaryWalls
    {
        type            noSlip;
    }
    atmosphere
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }
    floatingObject
    {
        type            movingWallVelocity;   //运动边界的速度设置 movingWallVelocity
        value           uniform (0 0 0);
    }
}

【0.org\pointDisplacement】

// 位移初、边值条件
dimensions      [0 1 0 0 0 0 0];
internalField   uniform (0 0 0);
boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
...   
    floatingObject
    {
        type            calculated;  // calculated
        value           uniform (0 0 0);
    }
}

【0.org\k】
【0.org\nut】
【0.org\epsilon】
对浮体结构表面使用了壁面函数边界。

【0.org\p_rgh】
【0.org\alpha.water】
参考【OpenFOAM】-olaFlow-算例1- baseWaveFlume

【constant\dynamicMeshDict】
动网格相关参数定义参考 Parameter Definitions - dynamicMotionSolverFvMesh

dynamicFvMesh       dynamicMotionSolverFvMesh;
motionSolverLibs    ("libsixDoFRigidBodyMotion.so"); // 引入运动求解器的库
motionSolver        sixDoFRigidBodyMotion;  // 指定运动求解器，六自由度刚体运动

sixDoFRigidBodyMotionCoeffs  // 指定六自由度刚体运动系数
{
    // 指定浮体的 patch 名称，主要有以下几层含义：
    // (1) 围绕该 patch 网格发生变形；
    // (2) 6Dof 体的速度作为流体求解的边界条件；
    // (3) 沿该 patch 表面积分来计算流体作用力。
    patches         (floatingObject);  // 浮体 patch 名称，而非边界条件
    
    // innerDistance 和 outerDistance 用于设置浮体运动时的网格变形。 
    // innerDistance 距离内的网格节点全部做刚体运动；
    // 介于 innerDistance 和 outerDistance 之间的网格会发生变形；
    // outerDistance 之外的网格不会发生变形。效果如图2所示。
    // outerDistance 必须大于 innerDistance。
    innerDistance   0.05;    
    outerDistance   0.35;

    centreOfMass    (0.5 0.45 0.35);  // 全局坐标系下的质心坐标，为初始位置
    
    // Cuboid dimensions
    Lx              0.3;
    Ly              0.2;
    Lz              0.5;

    // Density of the solid
    rhoSolid        500;  // 浮体的密度

    // Cuboid mass
    mass            #calc "$rhoSolid*$Lx*$Ly*$Lz";  // 计算质量

    // Cuboid moment of inertia about the centre of mass
    momentOfInertia #codeStream   // 计算惯性矩
    {
        codeInclude
        #{
            #include "diagTensor.H"
        #};

        code
        #{
            scalar sqrLx = sqr($Lx);
            scalar sqrLy = sqr($Ly);
            scalar sqrLz = sqr($Lz);
            os  <<
                $mass
               *diagTensor(sqrLy + sqrLz, sqrLx + sqrLz, sqrLx + sqrLy)/12.0;
        #};
    };
    
    report          on; // 输出开关，将 6DoF 体的状态写入时间文件。在非 0 时间文件夹下的 uniform 文件夹中提供了 6DoF 体的运动状态。
    
   // 以下这两个参数用于确保 6Dof 求解器的稳定性。在高加速度的情况下，6Dof 求解器会变得不稳定。在边界上，高加速度会反馈给流体高的速度，反过来会产生大的流体作用力，导致迅速的减速，这种情形会使 6Dof 求解器发散。因此OpenFOAM 采用以下两个参数来壁面上述情形发生。
   // accelerationDamping 主要用于消除突然加速产生的发散，该参数降低了浮体上的计算加速度，但它与加速度的值成比例，类似于阻尼系数，因此，突然的加速会形成较大的松弛。正常的加速度是典型的时间历程的结果，形成的加速相对较小。此参数在 0 至 1.0 之间取值，推荐值为 0.9 至 1.0。
   // accelerationRelaxation 会直接降低加速度，计算过程中实际采用的加速度会在计算值的基础上直接乘该松弛系数，适用于新衍生出的边界条件和物体速度，须小心使用，太小的值会使物体不能正确响应流体作用力，若计算的流体作用力与产生的运动之间的差异太大会引起发散。此参数在 0 至 1.0 之间取值，推荐值为 0.9 至 1.0。
    accelerationRelaxation 0.7;
    //accelerationDamping 0;

    solver  // 指定 6Dof 运动的求解算法
    {
        type Newmark;  // options:  Newmark(二阶隐式算法), CrankNicolson(二阶隐式算法，默认参数时，与Newmark法相同), symplectic(二阶显示算法)
    }
    
   // constrains 直接控制浮体许可的运动，与 restraints 不同, 后者施加的是反作用力。constrains 不产生任何力，直接限制物体运动。物体运动类型关键字包括：
   // (1) axis：仅约束转动，不约束平动。仅允许绕所定义的轴转动，将物体运动限制为一个绕该轴的转动自由度加三个平动自由度
   // (2) line：与 axis 类似，不过约束的是平动，转动不受限制。仅允许沿所定义的向量方向的平动，将物体运动限制为一个沿定义方向的平动自由度加三个转动自由度
   // (3) plane：允许 2D 平动，限制了刚体在第3维上的运动。运动方向由平面方向(法向量)定义，矢量的原点取在刚体重心的初始位置，允许刚体在定义平面内任意平动。
   // (4) point：限制了刚体平动，允许刚体绕所定义的旋转中心的任意方向的转动。不要求旋转中心位于刚体重心处，可用于定义球形连接( spherical joint)，使刚体绕节点转动。
   // (5) orientation：限制了刚体转动，平动不受约束。若未定义旋转中心，则默认刚体质心。
   // 以上约束可以组合使用，如 “point + axis”的组合模拟了销钉连接（铰接）。
   
    constraints   // 浮体的运动约束条件
    {
        // fixedPoint
        // {
        //     sixDoFRigidBodyMotionConstraint point;
        //     centreOfRotation (0.5 0.45 0.1);
        // }

        fixedLine  // 自定义约束名称
        {
            sixDoFRigidBodyMotionConstraint line; 
            centreOfRotation (0.5 0.45 0.1);  // 旋转中心
            direction (0 1 0);  // 允许沿 y 轴平动
        }

        fixedAxis
        {
            sixDoFRigidBodyMotionConstraint axis;
            axis (0 1 0);  // 允许绕 y 轴转动和平动
        }
    }
}

【constant\dynamicMeshDict.sixDoF】

dynamicFvMesh       dynamicMotionSolverFvMesh;
motionSolverLibs    ("libsixDoFRigidBodyMotion.so");
motionSolver        sixDoFRigidBodyMotion;

sixDoFRigidBodyMotionCoeffs
{
    patches         (floatingObject);
    innerDistance   0.05;
    outerDistance   0.35;

    centreOfMass    (0.5 0.45 0.35);

    // Cuboid dimensions
    Lx              0.3;
    Ly              0.2;
    Lz              0.5;

    // Density of the solid
    rhoSolid        500;

    // Cuboid mass
    mass            #calc "$rhoSolid*$Lx*$Ly*$Lz";

    // Cuboid moment of inertia about the centre of mass
    momentOfInertia #codeStream
    {
        codeInclude
        #{
            #include "diagTensor.H"
        #};

        code
        #{
            scalar sqrLx = sqr($Lx);
            scalar sqrLy = sqr($Ly);
            scalar sqrLz = sqr($Lz);
            os  <<
                $mass
               *diagTensor(sqrLy + sqrLz, sqrLx + sqrLz, sqrLx + sqrLy)/12.0;
        #};
    };

    report          on;
    accelerationRelaxation 0.7;
    //accelerationDamping 0;

    solver
    {
        type Newmark;
    }

    constraints
    {
        // fixedPoint
        // {
        //     sixDoFRigidBodyMotionConstraint point;
        //     centreOfRotation (0.5 0.45 0.1);
        // }

        fixedLine
        {
            sixDoFRigidBodyMotionConstraint line;
            centreOfRotation (0.5 0.45 0.1);
            direction (0 1 0);
        }

        fixedAxis
        {
            sixDoFRigidBodyMotionConstraint axis;
            axis (0 1 0);
        }
    }
}

【constant\dynamicMeshDict.solid】

dynamicFvMesh       dynamicMotionSolverFvMesh;
motionSolverLibs   ("librigidBodyMeshMotion.so");
motionSolver       rigidBodyMotion;

rigidBodyMotionCoeffs
{
    report          on;

    solver
    {
        type Newmark;
    }

    accelerationRelaxation 0.7;

    bodies
    {
        floatingObject
        {
            type            cuboid;
            parent          root;

            // Cuboid dimensions
            Lx              0.3;
            Ly              0.2;
            Lz              0.5;

            // Density of the cuboid
            rho             500;

            // Cuboid mass
            mass            #calc "$rho*$Lx*$Ly*$Lz";
            L               ($Lx $Ly $Lz);
            centreOfMass    (0 0 0.25);
            transform       (1 0 0 0 1 0 0 0 1) (0.5 0.45 0.1);

            joint
            {
                type            composite;
                joints
                (
                    {
                        type Py;
                    }
                    {
                        type Ry;
                    }
                );
            }

            patches         (floatingObject);
            innerDistance   0.05;
            outerDistance   0.35;
        }
    }
}

【constant\turbulenceProperties】

simulationType  laminar; // RAS;
RAS
{
    RASModel        kEpsilon;
    turbulence      on;
    printCoeffs     on;
}

【constant\waveDict】

waveType        regular;
waveTheory      StokesIII;
genAbs          1;
absDir          0.0;
nPaddles        1;
wavePeriod      1.0;
waveHeight      0.1;
waveDir         0.0;
wavePhase       1.57079633;
tSmooth         0.0;

【constant\g】
【constant\transportProperties】
以上均与一般水槽设置相同，参考【OpenFOAM】-olaFlow-算例1- baseWaveFlume

【system\blockMeshDict】
【system\controlDict】
【system\decomposeParDict】
【system\setFieldsDict】
参考【OpenFOAM】-olaFlow-算例1- baseWaveFlume

【system\topoSetDict】

actions
(
    {
        name    c0;
        type    cellSet;
        action  new;
        source  boxToCell; // 选中如下box框选中的cell，并将单元集命名为c0
        sourceInfo
        {
            box (0.35 0.35 0.1) (0.65 0.55 0.6);
        }
    }

    {
        name    c0;
        type    cellSet;
        action  invert;  // 反向选择，即求补集
    }
);

【system\fvSchemes】

ddtSchemes
{
    default         CrankNicolson 0.5;
}
gradSchemes
{
    default         Gauss linear;
}
divSchemes
{
    div(rhoPhi,U)  Gauss vanLeerV;
    div((rhoPhi|interpolate(porosity)),U)  Gauss vanLeerV;
    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss linear;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
    default         Gauss linear corrected;
}
interpolationSchemes
{
    default         linear;
}
snGradSchemes
{
    default         corrected;
}

【system\fvSolution】

solvers
{
    "alpha.water.*"
    {
        nAlphaCorr      2;
        nAlphaSubCycles 1;
        cAlpha          1;

        MULESCorr       yes;
        nLimiterIter    5;
        alphaApplyPrevCorr  yes;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0;
    }

    "pcorr.*"
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-5;
            relTol          0;
            smoother        DICGaussSeidel;
            cacheAgglomeration no;
        }

        tolerance       1e-05;
        relTol          0;
        maxIter         100;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-8;
        relTol          0.01;
        smoother        DIC;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-8;
            relTol          0;
            nVcycles        2;
            smoother        DICGaussSeidel;
            nPreSweeps      2;
        }

        tolerance       1e-8;
        relTol          0;
        maxIter         20;
    }

    "(U|k|epsilon)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
        nSweeps         1;
    }

    "(U|k|epsilon)Final"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-6;
        relTol          0;
        nSweeps         1;
    }
}

PIMPLE
{
    momentumPredictor   no;
    nOuterCorrectors    3;
    nCorrectors         1;
    nNonOrthogonalCorrectors 0;
    correctPhi          yes;
    moveMeshOuterCorrectors yes;
}

relaxationFactors
{
    equations
    {
        ".*" 1;
    }
}