CFD中文网

大家好，请问大家在使用sixDoFRigidBodyMotion进行流固耦合的时候，在仿真完毕之后如果想要续算是怎么设置的？

我使用的是OpenFOAM V2306，使用的求解以下是我的dynamicMeshDict文件：

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1912                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dynamicFvMesh   dynamicOversetFvMesh;

motionSolverLibs (sixDoFRigidBodyMotion);

motionSolver    sixDoFRigidBodyMotion;

dynamicFvMesh       dynamicOversetFvMesh;


sixDoFRigidBodyMotionCoeffs
{
    patches         (cylinder);
    innerDistance   100;
    outerDistance   101;

    mass            10;
   centreOfMass    (0.0 0.0 0.0);
   momentOfInertia (1 1 1.5625);
   g               (0 0 0);
   rho             rhoInf;
   rhoInf          1;
   report          on;

   solver
   {
       type Newmark;
   }

   constraints
   {
    yLine
    {
       sixDoFRigidBodyMotionConstraint   line;
       direction ( 0 1 0);
    }
    rotation
    {
        sixDoFRigidBodyMotionConstraint axis;
        axis (0 0 1);
    }
   }

  restraints
  {
      spring1
      {
          sixDoFRigidBodyMotionRestraint linearSpring;

          anchor          (0.0 0.0 0.0);
          refAttachmentPt (0.0 0.0 0.0);
          stiffness       10.96622711232151;
          damping         0.6283185307179585;
          restLength      0;
      }
      spring2
      {
       sixDoFRigidBodyMotionRestraint        linearAxialAngularSpring;
       axis                                  (0 0 1);
       stiffness                             1.7134729863002358;
       damping                               0.09817477042468103;
      }
  }
}
// ************************************************************************* //

以下是正常仿真300s之后的log.overPimpleDyMFoam中最后一步和sixDoFRigidMotion和overset有关的输出：

Restraint spring1:  attachmentPt - anchor (0 -0.367945 0) spring length 0.367945 force (-0 4.0584 -0)
Restraint spring2:  angle 0.516284 moment (0 0 -0.883591)
6-DoF rigid body motion
    Centre of rotation: (0 -0.368716 0)
    Centre of mass: (0 -0.368716 0)
    Orientation: (0.869769 -0.493459 0 0.493459 0.869769 0 0 0 1)
    Linear velocity: (0 -0.372264 0)
    Angular velocity: (0 0 -0.107739)
inverseDistance : detected 2 mesh regions
    zone:0 nCells:50721  voxels:(40 40 1) bb:(16.9781 3.77358 -0.500006) (20.9906 8.00001 0.500006)
    zone:1 nCells:2280  voxels:(40 40 1) bb:(-72 -40 -2.50004) (-40 -24 -1.49996)
Overset analysis : nCells : 53001
    calculated   : 52701
    interpolated : 189 (from local:0  mixed local/remote:0  remote:189)
    hole         : 111

如果我想续算，直接修改controlDict中的endTime显然是会出现问题的，它会使用原来的dynamicMeshDict中的数据作为300s的初值，导致计算无法正常进行。

而如果我参照最后一步的输出修改dynamicMeshDict文件：

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1912                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dynamicFvMesh   dynamicOversetFvMesh;

motionSolverLibs (sixDoFRigidBodyMotion);

motionSolver    sixDoFRigidBodyMotion;


sixDoFRigidBodyMotionCoeffs
{
    patches         (cylinder);
    innerDistance   100;
    outerDistance   101;

    mass            10;
   // centreOfMass    (0.0 0.0 0.0);
   centreOfMass    (0 -0.368716 0);
   momentOfInertia (1 1 1.5625);
   orientation (0.869769 -0.493459 0 0.493459 0.869769 0 0 0 1);
   velocity   (0 -0.372264 0);
   angularMomentum  (0 0 -0.1683421875);
   // angularVelocity (0 0 -0.107739);

   g               (0 0 0);
   rho             rhoInf;
   rhoInf          1;
   report          on;

   solver
   {
       type Newmark;
   }

   constraints
   {
    yLine
    {
       sixDoFRigidBodyMotionConstraint   line;
       direction ( 0 1 0);
    }
    rotation
    {
        sixDoFRigidBodyMotionConstraint axis;
        axis (0 0 1);
    }
   }

  restraints
  {
      spring1
      {
          sixDoFRigidBodyMotionRestraint linearSpring;

          anchor          (0.0 0.0 0.0);
          // anchor          (0 -0.368716 0);
          // refAttachmentPt (0.0 0.0 0.0);
          refAttachmentPt    (0 -0.368716 0);
          stiffness       10.96622711232151;
          damping         0.06283185307179585;
          restLength      0;
      }
      spring2
      {
       sixDoFRigidBodyMotionRestraint        linearAxialAngularSpring;
       axis                                  (0 0 1);
       stiffness                             1.7134729863002358;
       damping                               0.009817477042468103;
       // referenceOrientation                  0;
      }
  }
}
// ************************************************************************* //

首先根据最后一步的数据修改dynamicMeshDict文件我不知道会不会带来新的误差，然后这种方法也没有明显效果，以下是修改后续算的第一步输出：

Restraint spring1:  attachmentPt - anchor (0 -0.368716 0) spring length 0.368716 force (-0 4.06681 -0)
Restraint spring2:  angle 0.516062 moment (0 0 -0.883201)
6-DoF rigid body motion
    Centre of rotation: (0 -0.369639 0)
    Centre of mass: (0 -0.369639 0)
    Orientation: (0.869901 -0.493226 0 0.493226 0.869901 0 0 0 1)
    Linear velocity: (0 -0.371839 0)
    Angular velocity: (0 0 -0.10837)
inverseDistance : detected 2 mesh regions
    zone:0 nCells:50721  voxels:(230 230 1) bb:(-8.00004 -8.00004 -0.500036) (24 8.00004 0.500036)
    zone:1 nCells:2280  voxels:(230 230 1) bb:(-0.999904 -1.00093 -0.500003) (1.0001 0.99908 0.500003)
Overset analysis : nCells : 53001
    calculated   : 52683
    interpolated : 202 (from local:202  mixed local/remote:0  remote:0)
    hole         : 109

6-DoF相关的输出，通过修改dynamicMeshDict已经和前面比较接近了，但是zone 0， zone 1相关的输出有些差异，同时力系数突然变得很大：

forceCoeffs forceCoeffs write:
    Coefficient	Total	Pressure	Viscous	Internal
    Cd:	-734104	-733933	-171.164	0
    Cd(f):	-367052	-366966	-85.5821	0
    Cd(r):	-367052	-366966	-85.5821	0
    Cl:	-458471	-457934	-537.556	0
    Cl(f):	-295374	-294999	-375.408	0
    Cl(r):	-163097	-162935	-162.148	0
    CmPitch:	-66138.8	-66032.2	-106.63	0
    CmRoll:	-7.89329e-15	-7.89337e-15	7.8177e-20	0
    CmYaw:	8.40485e-15	8.40299e-15	1.86337e-18	0
    Cs:	2.00416e-14	2.00378e-14	3.83127e-18	0
    Cs(f):	1.84257e-14	1.84219e-14	3.779e-18	0
    Cs(r):	1.61596e-15	1.61591e-15	5.22695e-20	0

    writing force and moment coefficient files.

导致后续计算无法正常进行（时间步长变得很短）。

所以我想请教一下大家在使用sixDoFRigidMotion进行流固耦合仿真的时候，如果想要续算有什么好方法吗？因为似乎openfoam没有保存sixDoFRigidMotion的相关参数（如质心位置、速度、加速度等）。

@李东岳 老师，如果不使用refineMesh的话，结果会好很多，而且和论文中的Cd, Cl指数可以基本对上，但是在overset的边缘处还是会有一点问题。

我的oversetInterpolation的method使用的是inverseDistance.

大家好，我想仿真2D均匀流场中放置一个强迫旋转的椭圆，openfoam版本是v2306。

如果我使用AMI，那么可以得到和其他论文十分接近的结果，涡量图如下所示。

但是当我使用overset的时候，会出现比较大的偏差，在我对结果进行可视化之后，我发现在我使用局部加密的位置会出现很明显的间断，如下图所示，我觉得这可能是主要的误差来源。

以下是我的网格：

在我绘制涡量图时，会更加奇怪：

局部加密用的是openfoam自己的toposet和refineMesh，以下是对应的配置文件

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      refineMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Cells to refine; name of cell set
set refine3;

// Type of coordinate system:
// - global : coordinate system same for every cell. Usually aligned with
//   x,y,z axis. Specify in globalCoeffs section below.
// - patchLocal : coordinate system different for every cell. Specify in
//   patchLocalCoeffs section below.
// - fieldBased : uses the list of field names from the directions list for
//   selecting the directions to cut. Meant to be used with geometricCut, but
//   can also be used with useHexTopology.
coordinateSystem global;
//coordinateSystem patchLocal;
//coordinateSystem fieldBased;

// .. and its coefficients. x,y in this case. (normal direction is calculated
// as tan1^tan2)
globalCoeffs
{
    tan1 (1 0 0);
    tan2 (0 1 0);
}


// List of directions to refine, if global or patchLocal
directions
(
    tan1
    tan2
    normal
);

// List of directions to refine, if "fieldBased". Keep in mind that these
// fields must be of type "vectorField", not "volVectorField".
//directions
//(
//    radialDirectionFieldName
//    angularDirectionFieldName
//    heightDirectionFieldName
//);

// Whether to use hex topology. This will
// - if patchLocal: all cells on selected patch should be hex
// - split all hexes in 2x2x2 through the middle of edges.
useHexTopology  true;

// Cut purely geometric (will cut hexes through vertices) or take topology
// into account. Incompatible with useHexTopology
geometricCut    false;

// Write meshes from intermediate steps
writeMesh       false;

// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      topoSetDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

actions
(
    {
        name    refine3;
        type    cellSet;
        action  new;
        source  boxToCell;
        min     (-2 -2 -0.5);
        max     (2 2 0.5);
    }
);

// ************************************************************************* //

我想请问一下如果我使用局部加密网格，那么我有哪些需要注意的地方吗？

以上，希望各位大佬可以拨冗指点迷津。

各位老师好，我想请问一下在openfoam当中有比较容易的方法结合主动控制和被动控制吗？

例如目前我想用openfoamV2306的overPimpleDyMFoam进行这样一个仿真：

在均匀来流的流场中放置一个椭圆，只放开它垂直来流方向的自由度，但是给它一个恒定的角速度，让它在流场的作用下进行VIV(vortex induced vibration，涡激振动)，如下图所示：

我一开始的想法是直接在dynamicDict文件中缝合solidbody进行主动控制和sixDoFRigidBodyMotion进行被动运动，因为这两种方法单独使用都是可行的，但是结合起来无法进行运算，具体表现为计算发散，Co突增，以下是我尝试的dynamicMeshDict

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dynamicFvMesh       dynamicOversetFvMesh;

solver          multiSolidBodyMotionSolver;

solvers
{
    viv
    {
        motionSolver sixDoFRigidBodyMotion;
        patches         (cylinder);
        innerDistance   100;
        outerDistance   101;
        report          on;

        mass            2;
       centreOfMass    (0.0 0.0 0.0);
       momentOfInertia (1 1 1);

       g               (0 0 0);
       rho             rhoInf;
       rhoInf          1;

       solver
       {
           type Newmark;
       }

     constraints
     {

      yLine
      {
         sixDoFRigidBodyMotionConstraint   line;
         direction ( 0 1 0);
      }

      rotation
      {
          // sixDoFRigidBodyMotionConstraint orientation;
          sixDoFRigidBodyMotionConstraint axis;
          // centreOfRotation (0 0 0);
          axis (0 0 1);
      }
     }
    }
    overset_internal
    {
        motionSolver solidBody;
        solidBodyMotionFunction rotatingMotion;
        rotatingMotionCoeffs
        {
            origin      (0. 0. 0.);
            axis        (0 0 1);
            omega       1.570;
        }
    }
}

我也在网上搜索了很多相关的问题，但是都没有得到一个满意的答复。
如本站的sixDoFRigidBodyMotion问题求助这篇帖子中，李东岳老师提到，在openfoam中主动控制和被动移动的结合起来不难，但是没有给出具体的方法。

在github上有个feedback-flow-control-openfoam仓库，他的算例和我的需求比较相似。他提到在openfoam中并不支持同时使用主动和被动控制，而他是通过自定义边界条件来实现主动控制（也就是强迫旋转），这样可以只使用sixDoFRigidBodyMotion这一个solver，但是他的方法对我的问题并没有作用，自定义边界条件并没有改变我的物体的位置，也就是说在我看来，他的方法是一种等效的方法，他没有实际让物体运动起来，而是改变物体表面的速度分布来等效物体运动。我觉得这种方法对于轴对称的圆柱来说可能有效果，但是我的物体是一个椭圆，同时我也不想使用那些非惯性系之类的方法。

在cfd-online中其实也有一些相似的问题，但是那些问题都没有得到回答，比如：
https://www.cfd-online.com/Forums/openfoam-solving/238053-prescribed-rotation-6dof-possible.html
https://www.cfd-online.com/Forums/openfoam/253447-integrating-sixdofrigidbodymotion-rotatingmotion-openfoam.html

以上，希望各位老师可以不吝赐教，万分感激！