Difference between revisions of "Inline oscillating cylinder by Michael Alletto"

From OpenFOAM Wiki
Jump to navigation Jump to search
 
(14 intermediate revisions by 2 users not shown)
Line 3: Line 3:
 
* '''affiliation''': WRD MS
 
* '''affiliation''': WRD MS
 
* '''contact''': <mail address='michael.alletto@gmx.de' description='author'>click here for email address</mail>
 
* '''contact''': <mail address='michael.alletto@gmx.de' description='author'>click here for email address</mail>
* '''OpenFOAM versions''': v1912
+
* '''OpenFOAM versions''': v2212
 
* '''Published under''': CC BY-NC-ND license ([https://creativecommons.org/licenses creative commons licenses])
 
* '''Published under''': CC BY-NC-ND license ([https://creativecommons.org/licenses creative commons licenses])
  
Line 10: Line 10:
 
Go back to [https://wiki.openfoam.com/Meshing Meshing].
 
Go back to [https://wiki.openfoam.com/Meshing Meshing].
  
Go back to [https://wiki.openfoam.com/Turbulence Turbulence].
+
=[https://github.com/jnmlujnmlu/OpenFOAMTeaching/blob/master/MichaelAlletto Inline oscillating cylinder]=
 
 
 
 
=[https://github.com/jnmlujnmlu/OpenFOAMTeaching/blob/master/MichaelAlletto 2D flapping wing]=
 
  
 
You can download the case file [https://github.com/jnmlujnmlu/OpenFOAMTeaching/blob/master/MichaelAlletto/OscillatingCylinderSend.tar.gz OscillatingCylinderSend.tar.gz] here.
 
You can download the case file [https://github.com/jnmlujnmlu/OpenFOAMTeaching/blob/master/MichaelAlletto/OscillatingCylinderSend.tar.gz OscillatingCylinderSend.tar.gz] here.
Line 19: Line 16:
 
==Introduction==
 
==Introduction==
  
In this tutorial we will examine the inline oscillatory motion of a cylinder in a uniform ow. Inline means that the cylinder is oscillating parallel to the incoming flow. This kind of configuration is often used to validate solvers
+
In this tutorial we will examine the inline oscillatory motion of a cylinder in a uniform flow. Inline means that the cylinder is oscillating parallel to the incoming flow. This kind of configuration is often used to validate solvers
 
aiming to simulate large body motions (see e.g. [1, 4, 5]). The reason is that that despite its simplicity this configuration is associated with a variety of complex physical phenomena.
 
aiming to simulate large body motions (see e.g. [1, 4, 5]). The reason is that that despite its simplicity this configuration is associated with a variety of complex physical phenomena.
  
Line 25: Line 22:
  
 
Another interesting phenomena is the increase of the lift coefficient C<sub>l</sub> when the cylinder is oscillating at 2f<sub>st</sub> [6]. For frequencies which are either lower or higher than 2f<sub>st</sub> C<sub>l</sub> is much lower [6]. [3] made simulations with an in line oscillating cylinder at a Reynolds number of Re = 175. For this Re number the ow around a fixed cylinder is two dimensional. For an oscillation frequency f<sub>c</sub> of 2f<sub>st</sub> and an amplitude of 65% of the cylinder diameter the frequency spectrum of the lift force increases considerably. For lower amplitudes of the oscillation only a few peaks are present in the frequency spectrum.
 
Another interesting phenomena is the increase of the lift coefficient C<sub>l</sub> when the cylinder is oscillating at 2f<sub>st</sub> [6]. For frequencies which are either lower or higher than 2f<sub>st</sub> C<sub>l</sub> is much lower [6]. [3] made simulations with an in line oscillating cylinder at a Reynolds number of Re = 175. For this Re number the ow around a fixed cylinder is two dimensional. For an oscillation frequency f<sub>c</sub> of 2f<sub>st</sub> and an amplitude of 65% of the cylinder diameter the frequency spectrum of the lift force increases considerably. For lower amplitudes of the oscillation only a few peaks are present in the frequency spectrum.
In order to simulate this kind of ow congurations, OpenFOAM offers two different approaches: the overset methodology and the morphing mesh method. The latter can be used if the cylinder motion is not to large to deteriorate the mesh quality. Regarding the overset method we have one or more meshes which displace relative to a background mesh. The coupling between the moving meshes and the background mesh is achieved via interpolation of the  
+
In order to simulate this kind of flow configurations, OpenFOAM offers two different approaches: the overset methodology and the morphing mesh method. The latter can be used if the cylinder motion is not too large to deteriorate the mesh quality. Regarding the overset method we have one or more meshes which displace relative to a background mesh. The coupling between the moving meshes and the background mesh is achieved via interpolation of the  
 
ow quantities. Regarding the morphing mesh, the connectivity between the cells remain unchanged.However the position of the points building a cell can move relative to each other.
 
ow quantities. Regarding the morphing mesh, the connectivity between the cells remain unchanged.However the position of the points building a cell can move relative to each other.
  
Line 32: Line 29:
 
==Set up==
 
==Set up==
  
The Simulation setup is inspired by the experiments of [3]. The experiments consisted in a flat plate which rotates  about the stroke axis with an angle of 0 <= &Phi; <= 90°  and pitches about the rotational axis with an angle of &beta; = +- 50°. The details of the kinematics are described in  [3]. The wing span was R = 107mm and the thickness was s = 3.4mm (see [4] since in [3] the thickness of the plate was not described). The distance from the rotational axis to the wing root was &Delta;R = 71.5mm (private communications). In order to perform a two--dimensional simulation, the rotatory motion of the point located at mid-span is projected on a horizontal plane. That means that the displacement of the point located at the center of the wing in the experiment is the same as the horizontal displacement of the wing in the simulation. The pitching motion is the same in the experiment and in the simulation.  
+
The figure below shows the setup of the simulation: At the inflow a uniform velocity of U<sub>inf</sub> = 1 m/s in x-direction is applied. The pressure is set to zero gradient at this boundary. At the top and bottom walls a slip boundary condition is applied. At the outlet the pressure is to be set equal to zero and for the velocity a zero gradient boundary condition is applied. The kinematic viscosity &nu; = 0.01 m^2/s is set in order that the Reynolds number Re = U<sub>inf</sub> D / &nu; based on the cylinder diameter D = 1 m is equal to Re = 100. For this Reynolds number the flow over a stationary cylinder is two dimensional. For this reason the flow is assumed to be two dimensional. This is imposed by using only one cell in span wise direction and setting both boundary conditions in this direction to the type empty. At the cylinder wall a no slip boundary condition is applied for the velocity and zero gradient for the pressure. The cylinder oscillates in x-direction with a prescribed sinosoidal motion: x<sub>c</sub> = A  sin(2&pi; f<sub>c</sub> t) = A sin (&omega; t). Four different frequencies are computed: f<sub>c</sub> = 0, 0.16, 0.33, 0.66 s^(-1). The amplitude is set to A = 0.14m. The purpose of this tutorial was to check whether the increase of the lift at f<sub>c</sub> = 2f<sub>st</sub> can be captured and if both techniques to simulate moving bodies deliver comparable results.
  
===Mesh generation===
+
[[File:Setup_alettoinline.png|450px|center|setup of the simulation]]
  
Since the interpolation errors in overset grids is proportional to the cell size difference between moving and background mesh (of course also to the cell size itself), the purpose of the mesh generation was to minimize this errors. In order to achieve this goal, the size of the cells of the background mesh and the moving grid have to be similar in the region where the interpolation is performed. For this reason the moving mesh was generated with snappyHexMesh in combination with different refinement regions. Since snappyHexMesh can perform only three--dimensional refinements (a hexahedron is split in eight cells with no option for the user to impose only splitting in two directions) the utility extrudeMesh is used. With this utility the 2D mesh at the symmetry plane can be extruded having a mesh with only one cell in span wise direction. The mesh is shown in the figure below.
+
== Mesh generation ==
  
[[File:MovingMesh_alettoflap.png|900px|center|mesh of the flapping wing]]
+
The mesh used for the overset mesh simulation is shown in the figure below. We can see that it is formed by the mesh surrounding the cylinder and the background mesh. The mesh surrounding the cylinder is generated by simply stretching the circular form of the cylinder outward. The background mesh consist of a rectangular block mesh with different refinement zones. In the outer part of the rectangle where the gradients of the flow field are smaller, a coarser mesh is used. Closer to the region surrounding the cylinder the mesh is finer. Note that when generating an overset mesh, it has to be paid attention that the cells at the intersection of the moving and the background mesh should have roughly the same size to minimize the interpolation errors.
  
===Simulation set up===
+
[[File:meshOverset2p_alettoinline.png|900px|center|mesh used for the overset method]]
  
As for the other tutorials presented, the different zones (background mesh and the moving meshes) should be differentiated by the scalar field zoneID. This field should have a different number (starting with zero) for each mesh (see figure below). The overset method implemented in OpenFOAM detects itself the different cell types. The cells adjacent to the patch denominated as overset are chosen as interpolated. For this cells the matrix coefficients are not derived from the finite volume discretization but using one of overset interpolation method chosen by the user. The cells which are inside the wing are flagged as hole cells and exclude from the calculation. The other cells are treated as calculated (the matrix coefficients are derived from a finite volume discretization). The cell types are visualized on the right of the figure below.
+
The figure below shows the morphing mesh. It is generated using the utility blockMesh. Note that the dictionary for the generation of the mesh was taken from [https://curiosityfluids.com/2016/07/19/oscillating-cylinder-in-crossflow-pimpledymfoam/ here] and slightly modified in order to adept it to the present test case. It is evident that the constitution of the mesh is much more complex compared to the overset mesh. The number of cells of both meshes are similar: roughly 10 000 for the overset mesh and roughly 12 000 for the morphing mesh.
  
[[File:ZoneIDflap.png|900px|center|zoneID and cell types (blue calculated, gray interpolated and red hole cells)]]
+
[[File:meshMorphing2_alettoinline.png|900px|center|mesh used for the with the morphing mesh method]]
  
The kinematics of the wing are imposed by choosing the adequat motion solver. Since the wing executes a sinusoidal movement in x-direction and pitches at the and of the half stroke with a linear change of the angle &beta; (see figure below on wing kinematics where the evolution of the x-position of and the angle &beta; during the first cycle is shown), the tabulated6DoFMotion function is chosen. This function allows to impose the temporal evolution of the displacement in x-,y- and z-direction and also the rotation about the same axis. The displacements and rotation are relative to the initial position.
+
== Simulation set up ==
The rotation is described about the point denoted as CofG. Since the rotation in the experiment was about an axis located a forth of the chord length from the leading edge, CofG was set to at this position in the initial state. The dynamic mesh dictionary is shown below.
 
  
 +
Since the procedure how to set up a overset mesh is described in other tutorials, here we will concentrate on the difference on how to describe the body motion between the overset mesh and the morphing mesh.
  
 +
Regarding the overset mesh the solid body motion has to be prescribed in the constant/dynamicMeshDict:
  
  <TT>dynamicFvMesh      dynamicOversetFvMesh;  
+
  <TT>motionSolverLibs ( "libfvMotionSolvers.so" );
 +
solver          displacementLaplacian;
 +
displacementLaplacianCoeffs
 +
{
 +
    diffusivity    uniform 1;
 +
}
 +
dynamicFvMesh      dynamicOversetFvMesh;
 +
dynamicOversetFvMeshCoeffs
 +
{
 +
//    layerRelax 0.3;
 +
}
 
  solver          multiSolidBodyMotionSolver;
 
  solver          multiSolidBodyMotionSolver;
 
  multiSolidBodyMotionSolverCoeffs
 
  multiSolidBodyMotionSolverCoeffs
 
  {
 
  {
     movingZone
+
    movingzone
 +
    {
 +
      solidBodyMotionFunction oscillatingLinearMotion;
 +
      oscillatingLinearMotionCoeffs
 +
      {
 +
        amplitude ( -0.14 0 0);
 +
        omega 1.04;
 +
      }
 +
    }
 +
}
 +
</TT>
 +
 
 +
The type of motion and the mesh zone which has to perform this prescribed motion has to be specified. In our case the mesh zone named movingzone has to describe an oscillatory motion in x-direction with an amplitude of A = 0.14 m and a frequency of &omega; = 1.04 rad/s.
 +
 
 +
Regarding the morphing mesh, the displacement is described via the boundary condition of the field pointDisplacement:
 +
 
 +
<TT>
 +
dimensions      [0 1 0 0 0 0 0];
 +
internalField  uniform (0 0 0);
 +
boundaryField
 +
{
 +
     topAndBottom
 
     {
 
     {
         solidBodyMotionFunction tabulated6DoFMotion;
+
         type            uniformFixedValue;
         CofG            (0.00546 0.00651 0);
+
         uniformValue    (0 0 0);
        timeDataFileName "<constant>/6DoF.dat";
 
 
     }
 
     }
 +
    inlet
 +
    {
 +
        type            uniformFixedValue;
 +
        uniformValue    (0 0 0);
 +
    }
 +
    outlet
 +
    {
 +
        type            uniformFixedValue;
 +
        uniformValue    (0 0 0);
 +
    }
 +
    cylinder
 +
    {
 +
        type            oscillatingDisplacement;
 +
        amplitude      ( -0.14 0 0 );
 +
        omega          1.04;
 +
        value          uniform (0 0 0);
 +
    }
 +
    frontAndBack
 +
    {
 +
        type            empty;
 +
    }
 +
}
 +
</TT>
 +
 +
We see that the in this case the displacement of the cylinder is imposed by specifying the type of motion the boundary patch named cylinder has to perform. All other boundary patches do not move and hence the value of the field pointDisplacement is set to uniform (0 0 0).
 +
 +
The method used to relax the cylinder movement towards the fixed outer boundary is specified in the constant/dynamicMeshDict:
 +
 +
<TT>
 +
dynamicFvMesh  dynamicMotionSolverFvMesh;
 +
motionSolverLibs ( "libfvMotionSolvers.so" );
 +
solver          displacementLaplacian;
 +
displacementLaplacianCoeffs
 +
{
 +
    diffusivity  inverseDistance 1(cylinder);
 
  }
 
  }
 
  </TT>
 
  </TT>
  
[[File:Kinematics_alettoflap.png|400px|center|kinematics of the wing]]
+
== Results ==
 +
 
 +
The table below shows the comparison between the present simulations and reference simulation and experiments for a fixed cylinder. The variables which are compared are the Strouhal number, the mean drag coefficient C<sub>d,mean</sub> and the maximum lift coefficient C<sub>l,max</sub>. The Strouhal number is computed by analyzing the signal of the lift (force in z-direction). Since the purpose of this tutorial is to evaluate the capability of OpenFOAM to simulate an oscillating cylinder, the first necessary step is to check the quality of the simulations for a fixed cylinder. It is evident that the agreement is reasonable. The results of the overset mesh and the mesh generated with blockMesh are very close. This means that the overset interpolation does not deteriorate the quality of the results for this case. The results shown are computed with a python script which analyzes the forces computed around the cylinder. Of course the python script is also provided within the downloadable files.
 +
 
 +
[[File:table1_alettoinline.png|500px|center|Comparison of the Strouhal number, the mean drag coefficient C<sub>d,mean</sub> and the maximum lift coefficient C<sub>l,max</sub> with experiments (exp) and other simulation (sim) ]]
 +
 
 +
The next table below shows the comparison of the mean drag coefficient C<sub>d,mean</sub> and the maximum lift coefficient C<sub>l,max</sub> with other simulation (sim) for different oscillation frequencies f<sub>c</sub>. For the present simulations also the Strouhal number computed with the frequency of the lift signal with the maximum power is shown. It can see seen that both mean drag coefficient C<sub>d,mean</sub> and the maximum lift coefficient C<sub>l,max</sub> have a maximum at a frequency twice the shedding frequency of a fixed cylinder f<sub>c</sub> = 2 f<sub>st</sub>. The difference between the morphing mesh and the overset mesh are marginal and the agreement with the reference simulations are reasonable. All simulation (the one of this tutorials and also the reference simulations) report an almost tripling of the lift coefficient with respect to the stationary cylinder (see table above). The Strouhal number computed with the frequency with the maximum power of the lift signal remains unchanged when the oscillation frequency f<sub>c</sub> is increased.
  
===Case specific postprocessing===
+
[[File:table2_alettoinline.png|350px|center|Comparison of the mean drag coefficient  C<sub>d,mean</sub> and the maximum lift coefficient C<sub>l,max</sub> with other simulation (sim) for different oscillation frequencies f<sub>c</sub>]]
  
In the reference experiments of [3] the phase averaged lift coefficient C<sub>l</sub>  = F<sub>y</sub>/(0.5&rho;U*<sup>2</sup>) and the phase averaged drag coefficient C<sub>d</sub> = F<sub>x</sub>/(0.5&rho;U*<sup>2</sup>) were measured. U* is defined as U* = πfR where f is the flapping frequency and R is the wing span. In the following simulations the force coefficients are computed by scaling the forces by the mean flapping velocity and the density. Since the forces scale with the square of the velocity, the measured force coefficients have to be rescaled by the mean flapping velocity at the second moment of area radius R<sub>2</sub> (see e.g. [5]) to be compared with the present simulation. A small python script is included in the tutorial which computes the phase averaged force coefficients of the simulation and compares it with the measurements.  
+
The figure below shows the comparison of the temporal evolution of drag coefficient in the last 20 s of the simulation. It is evident that for increasing frequency f<sub>c</sub> the amplitude of the drag coefficient increases since the velocity at which the cylinder is displacing through the fluid increases. The results of the morphing mesh and the overset mesh are very close.
  
Furthermore phase averaged PIV measurements of the normalized vorticity &omega; = &omega;* c/U* where shown in the reference experiment. c represents the chord length. In order to compare this quantity with the experiments, a small coded function object is included in the controlDict which does the phase averaging at phase times of t = 0.2, 0.64, 1.24, 1.6 and 1.8 s. This times correspond to the non--dimensional times of t/T = 0.05, 0.16, 0.31, 0.4 and 0.45 of the reference experiments of [3].
+
[[File:cdall_alettoinline.png|900px|center|comparison of the drag coeffcient of the morphing mesh (left) and overset mesh (right)]]
  
==Results==
+
The figure below  shows the comparison of the temporal evolution of lift coefficient. It is obvious that for a frequency of f<sub>c</sub> = 2 f<sub>st</sub> the amplitude of the lift coefficient is highest. Again the results of the morphing mesh and the overset mesh are very similar.
  
The figure below shows the temporal evolution of the force coefficient (left) and the comparison of the phase averaged force coefficients with the experiments (right). Looking at the evolution of the force coefficient it is evident that the signal is not completely periodic as expected for this low Reynolds number. The interaction of the wing with the residual vorticities from the previous strokes seems to lead to a not completely periodic evolution of the forces on the wing. The agreement of the phase averaged coefficients with the experiments is satisfactory especially if we keep in mind that a 2D simulation is compared with a 3D experiment. In the simulation the influence of the flow from the root of the wing to the tip which is induced by a mean spanwise pressure gradient is completely disregarded. The peak of the force coefficients at the end of the cycle probably results from an too high influence of the wing rotation at the end of the cycle. 
+
[[File:clall_alettoinline.png|900px|center|comparison of the liftcoeffcient of the morphing mesh (left) and overset mesh (right)]]
  
[[File:ClCd_allettoflap.png|900px|center|left: temporal evolution of C<sub>d</sub> (blue) and C<sub>l</sub> (green); right: phase averaged C<sub>d</sub> and C<sub>l</sub>]]
+
The figure below  shows the power spectrum of the lift for the four oscillation frequencies f<sub>c</sub> analyzed and both mesh types. The results of both mesh types are again very similar. The amount of peaks increase with increasing f<sub>c</sub>.
  
According to [3] the cycle of a flapping wing can be divided into five stages: the emergence of the leading edge vortex (LEV), the growth of the LEV, the lift--off of the LEV, the breakdown and decay of the LEV. For each of this stages the phase averaged scaled vorticity &omega; (left column of the figure below) and the pressure scaled by the mean flapping velocity squared and the density (right column of the figure below) is shown. The interpretation of the results and its comparison with the reference experiment is left to the interested reader.
+
[[File:frall_alettoinline.png|900px|center|comparison of the power spectrum of lift for the morphing mesh (left) and the overset mesh (right)]]
  
[[File:Bigpic_allettoflap.png|600px|center|phase averaged normalized vorticity &omega; (left column) and pressure (right column)]]
+
A final remark should be made regarding the computational time. For the fixed cylinder and the case with a oscillation frequency f<sub>c</sub> = 0.16, the simulations with the morphing mesh are much faster. For increasing frequency however, the maximum Courant number Co<sub>max</sub> had to be decreased from 0.75 for f<sub>c</sub> = 0.16 to Co<sub>max</sub> = 0.5 for f<sub>c</sub> = 0.3 and Co<sub>max</sub>  = 0.25 for f<sub>c</sub> = 0.66 for the morphing mesh. The purpose was to avoid nonphysical pressure peaks which occurred for higher Co<sub>max</sub>. For this reason the simulation time for the cases with the highest frequency was comparable between both mesh types. For the overset mesh the maximum Courant number Co<sub>max</sub> = 0.75 was the same for all frequencies.
  
 
==References==
 
==References==
[1] Dominic D Chandar, Bharathi Boppana, and Vasanth Kumar. A comparative study of di�erent overset grid solvers between openfoam, star-ccm+ and ansys-fluent. In 2018 AIAA Aerospace Sciences Meeting, page 1564, 2018.
+
[1] Dominic D Chandar, Bharathi Boppana, and Vasanth Kumar. A comparative study of different overset grid solvers between openfoam, star-ccm+ and ansys-fluent. In 2018 AIAA Aerospace Sciences Meeting, page 1564, 2018.
  
[2] Owen M Gri�n and MS Hall. Vortex shedding lock-on and flow control in bluff body wakes. 1991.
+
[2] Owen M Griffn and MS Hall. Vortex shedding lock-on and flow control in bluff body wakes. 1991.
  
 
[3] Justin S Leontini, David Lo Jacono, and Mark C Thompson. A numerical study of an inline oscillating cylinder in a free stream. Journal of Fluid Mechanics, 688:551-568, 2011.
 
[3] Justin S Leontini, David Lo Jacono, and Mark C Thompson. A numerical study of an inline oscillating cylinder in a free stream. Journal of Fluid Mechanics, 688:551-568, 2011.

Latest revision as of 19:59, 22 January 2023

Go back to Collection by authors.

Go back to Meshing.

Inline oscillating cylinder

You can download the case file OscillatingCylinderSend.tar.gz here.

Introduction

In this tutorial we will examine the inline oscillatory motion of a cylinder in a uniform flow. Inline means that the cylinder is oscillating parallel to the incoming flow. This kind of configuration is often used to validate solvers aiming to simulate large body motions (see e.g. [1, 4, 5]). The reason is that that despite its simplicity this configuration is associated with a variety of complex physical phenomena.

One of this phenomena are frequency lock in: When the frequency of the oscillating cylinder is close to twice the shedding frequency of the stationary cylinder fst, the frequency of the vortices shed became equal to the one of the cylinder. The amplitude of the cylinder motion which is required to maintain the lock in regime is increasing with departure from twice the shedding frequency of a fixed cylinder [2].

Another interesting phenomena is the increase of the lift coefficient Cl when the cylinder is oscillating at 2fst [6]. For frequencies which are either lower or higher than 2fst Cl is much lower [6]. [3] made simulations with an in line oscillating cylinder at a Reynolds number of Re = 175. For this Re number the ow around a fixed cylinder is two dimensional. For an oscillation frequency fc of 2fst and an amplitude of 65% of the cylinder diameter the frequency spectrum of the lift force increases considerably. For lower amplitudes of the oscillation only a few peaks are present in the frequency spectrum. In order to simulate this kind of flow configurations, OpenFOAM offers two different approaches: the overset methodology and the morphing mesh method. The latter can be used if the cylinder motion is not too large to deteriorate the mesh quality. Regarding the overset method we have one or more meshes which displace relative to a background mesh. The coupling between the moving meshes and the background mesh is achieved via interpolation of the ow quantities. Regarding the morphing mesh, the connectivity between the cells remain unchanged.However the position of the points building a cell can move relative to each other.

Since both methods described above can be used to simulate the flow of an inline oscillating cylinder, both methods can be compared against each other besides the validation against experiments and other simulations. The comparison of this two methods and the validation against experiments is the purpose of this tutorial.

Set up

The figure below shows the setup of the simulation: At the inflow a uniform velocity of Uinf = 1 m/s in x-direction is applied. The pressure is set to zero gradient at this boundary. At the top and bottom walls a slip boundary condition is applied. At the outlet the pressure is to be set equal to zero and for the velocity a zero gradient boundary condition is applied. The kinematic viscosity ν = 0.01 m^2/s is set in order that the Reynolds number Re = Uinf D / ν based on the cylinder diameter D = 1 m is equal to Re = 100. For this Reynolds number the flow over a stationary cylinder is two dimensional. For this reason the flow is assumed to be two dimensional. This is imposed by using only one cell in span wise direction and setting both boundary conditions in this direction to the type empty. At the cylinder wall a no slip boundary condition is applied for the velocity and zero gradient for the pressure. The cylinder oscillates in x-direction with a prescribed sinosoidal motion: xc = A sin(2π fc t) = A sin (ω t). Four different frequencies are computed: fc = 0, 0.16, 0.33, 0.66 s^(-1). The amplitude is set to A = 0.14m. The purpose of this tutorial was to check whether the increase of the lift at fc = 2fst can be captured and if both techniques to simulate moving bodies deliver comparable results.

setup of the simulation

Mesh generation

The mesh used for the overset mesh simulation is shown in the figure below. We can see that it is formed by the mesh surrounding the cylinder and the background mesh. The mesh surrounding the cylinder is generated by simply stretching the circular form of the cylinder outward. The background mesh consist of a rectangular block mesh with different refinement zones. In the outer part of the rectangle where the gradients of the flow field are smaller, a coarser mesh is used. Closer to the region surrounding the cylinder the mesh is finer. Note that when generating an overset mesh, it has to be paid attention that the cells at the intersection of the moving and the background mesh should have roughly the same size to minimize the interpolation errors.

mesh used for the overset method

The figure below shows the morphing mesh. It is generated using the utility blockMesh. Note that the dictionary for the generation of the mesh was taken from here and slightly modified in order to adept it to the present test case. It is evident that the constitution of the mesh is much more complex compared to the overset mesh. The number of cells of both meshes are similar: roughly 10 000 for the overset mesh and roughly 12 000 for the morphing mesh.

mesh used for the with the morphing mesh method

Simulation set up

Since the procedure how to set up a overset mesh is described in other tutorials, here we will concentrate on the difference on how to describe the body motion between the overset mesh and the morphing mesh.

Regarding the overset mesh the solid body motion has to be prescribed in the constant/dynamicMeshDict:

motionSolverLibs ( "libfvMotionSolvers.so" );
solver          displacementLaplacian;
displacementLaplacianCoeffs
{
    diffusivity     uniform 1;
}
dynamicFvMesh       dynamicOversetFvMesh;
dynamicOversetFvMeshCoeffs
{
//    layerRelax 0.3;
}
solver          multiSolidBodyMotionSolver;
multiSolidBodyMotionSolverCoeffs
{
    movingzone
    { 
     solidBodyMotionFunction oscillatingLinearMotion;
     oscillatingLinearMotionCoeffs
     {
        amplitude ( -0.14 0 0);
        omega 1.04;
     }
    }
}

The type of motion and the mesh zone which has to perform this prescribed motion has to be specified. In our case the mesh zone named movingzone has to describe an oscillatory motion in x-direction with an amplitude of A = 0.14 m and a frequency of ω = 1.04 rad/s.

Regarding the morphing mesh, the displacement is described via the boundary condition of the field pointDisplacement:


dimensions      [0 1 0 0 0 0 0];
internalField   uniform (0 0 0);
boundaryField
{
   topAndBottom
   {
       type            uniformFixedValue;
       uniformValue    (0 0 0);
   }
   inlet
   {
       type            uniformFixedValue;
       uniformValue    (0 0 0);
   }
   outlet
   {
       type            uniformFixedValue;
       uniformValue    (0 0 0);
   }
   cylinder
   {
       type            oscillatingDisplacement;
       amplitude       ( -0.14 0 0 );
       omega           1.04; 
       value           uniform (0 0 0);
   }
   frontAndBack
   {
       type            empty;
   }
}

We see that the in this case the displacement of the cylinder is imposed by specifying the type of motion the boundary patch named cylinder has to perform. All other boundary patches do not move and hence the value of the field pointDisplacement is set to uniform (0 0 0).

The method used to relax the cylinder movement towards the fixed outer boundary is specified in the constant/dynamicMeshDict:


dynamicFvMesh   dynamicMotionSolverFvMesh;
motionSolverLibs ( "libfvMotionSolvers.so" );
solver          displacementLaplacian;
displacementLaplacianCoeffs
{
    diffusivity  inverseDistance 1(cylinder);
}

Results

The table below shows the comparison between the present simulations and reference simulation and experiments for a fixed cylinder. The variables which are compared are the Strouhal number, the mean drag coefficient Cd,mean and the maximum lift coefficient Cl,max. The Strouhal number is computed by analyzing the signal of the lift (force in z-direction). Since the purpose of this tutorial is to evaluate the capability of OpenFOAM to simulate an oscillating cylinder, the first necessary step is to check the quality of the simulations for a fixed cylinder. It is evident that the agreement is reasonable. The results of the overset mesh and the mesh generated with blockMesh are very close. This means that the overset interpolation does not deteriorate the quality of the results for this case. The results shown are computed with a python script which analyzes the forces computed around the cylinder. Of course the python script is also provided within the downloadable files.

Comparison of the Strouhal number, the mean drag coefficient Cd,mean and the maximum lift coefficient Cl,max with experiments (exp) and other simulation (sim)

The next table below shows the comparison of the mean drag coefficient Cd,mean and the maximum lift coefficient Cl,max with other simulation (sim) for different oscillation frequencies fc. For the present simulations also the Strouhal number computed with the frequency of the lift signal with the maximum power is shown. It can see seen that both mean drag coefficient Cd,mean and the maximum lift coefficient Cl,max have a maximum at a frequency twice the shedding frequency of a fixed cylinder fc = 2 fst. The difference between the morphing mesh and the overset mesh are marginal and the agreement with the reference simulations are reasonable. All simulation (the one of this tutorials and also the reference simulations) report an almost tripling of the lift coefficient with respect to the stationary cylinder (see table above). The Strouhal number computed with the frequency with the maximum power of the lift signal remains unchanged when the oscillation frequency fc is increased.

Comparison of the mean drag coefficient Cd,mean and the maximum lift coefficient Cl,max with other simulation (sim) for different oscillation frequencies fc

The figure below shows the comparison of the temporal evolution of drag coefficient in the last 20 s of the simulation. It is evident that for increasing frequency fc the amplitude of the drag coefficient increases since the velocity at which the cylinder is displacing through the fluid increases. The results of the morphing mesh and the overset mesh are very close.

comparison of the drag coeffcient of the morphing mesh (left) and overset mesh (right)

The figure below shows the comparison of the temporal evolution of lift coefficient. It is obvious that for a frequency of fc = 2 fst the amplitude of the lift coefficient is highest. Again the results of the morphing mesh and the overset mesh are very similar.

comparison of the liftcoeffcient of the morphing mesh (left) and overset mesh (right)

The figure below shows the power spectrum of the lift for the four oscillation frequencies fc analyzed and both mesh types. The results of both mesh types are again very similar. The amount of peaks increase with increasing fc.

comparison of the power spectrum of lift for the morphing mesh (left) and the overset mesh (right)

A final remark should be made regarding the computational time. For the fixed cylinder and the case with a oscillation frequency fc = 0.16, the simulations with the morphing mesh are much faster. For increasing frequency however, the maximum Courant number Comax had to be decreased from 0.75 for fc = 0.16 to Comax = 0.5 for fc = 0.3 and Comax = 0.25 for fc = 0.66 for the morphing mesh. The purpose was to avoid nonphysical pressure peaks which occurred for higher Comax. For this reason the simulation time for the cases with the highest frequency was comparable between both mesh types. For the overset mesh the maximum Courant number Comax = 0.75 was the same for all frequencies.

References

[1] Dominic D Chandar, Bharathi Boppana, and Vasanth Kumar. A comparative study of different overset grid solvers between openfoam, star-ccm+ and ansys-fluent. In 2018 AIAA Aerospace Sciences Meeting, page 1564, 2018.

[2] Owen M Griffn and MS Hall. Vortex shedding lock-on and flow control in bluff body wakes. 1991.

[3] Justin S Leontini, David Lo Jacono, and Mark C Thompson. A numerical study of an inline oscillating cylinder in a free stream. Journal of Fluid Mechanics, 688:551-568, 2011.

[4] Chuan-Chieh Liao, Yu-Wei Chang, Chao-An Lin, and JM McDonough. Simulating flows with moving rigid boundary using immersed-boundary method. Computers & Fluids, 39(1):152-167, 2010. 11

[5] Shen-Wei Su, Ming-Chih Lai, and Chao-An Lin. An immersed boundary technique for simulating complex flows with rigid boundary. Computers & Fluids, 36(2):313-324, 2007.

[6] Y Tanida, A Okajima, and Y Watanabe. Stability of a circular cylinder oscillating in uniform ow or in a wake. Journal of Fluid Mechanics, 61(4):769-784, 1973.

[7] Charles HK Williamson. Vortex dynamics in the cylinder wake. Annual review of Fluid mechanics, 28(1):477-539, 1996.