算例路径: olaFlow\tutorials\waveFloatingObject
算例描述: 波浪作用下的浮体的刚体运动,属于流固耦合(FSI)问题
学习目标: 动网格设置和使用,网格变形控制,浮体的物理参数设置,浮体做刚体运动的约束设置
算例快照:
文件结构:
.
├── 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;
}
}