Simulation of Vortex-Induced Vibration on a 2D Cylinder with a Degree of Freedom

[Manu-CFD]

Hello, I am trying to implement vortex-induced vibration on a cylinder for my Engineering thesis, and I have some questions about what codes need to be modified. In principle, what I did was modify the geocomplex_cyl subroutine of the BC-Cylinder code and the ana_x_cyl and ana_y_cyl subroutines of the ibm code, directly solving the mass-spring-damper differential equation with Runge-Kutta 4 to obtain the vertical position of the cylinder, from the following way:

L1 = (1/M)(LIFT(t/dt)-Cubcy-Kceyy)
K1 = ubcy
L2 = (1/M)(LIFT(t/dt+1)-C*(ubcy+0.5*(2dt)L1)-K(ceyy+0.5(2dt)K1))
K2 = (ubcy+0.5(2dt)L1)
L3 = (1/M)(LIFT(t/dt+1)-C*(ubcy+0.5*(2dt)L2)-KR(ceyy+0.5(2dt)K2))
K3 = (ubcy+0.5(2dt)L2)
L4 = (1/M)(LIFT(t/dt+2)-C*(ubcy+L3*(2dt))-K(ceyy+K3*(2dt))
K4 = (ubcy+L3(2dt))
ceyy = ceyy + (2dt)0.16667(K1+2K2+2K3+K4)
ubcy = ubcy + (2dt)0.16667(L1+2L2+2*L3+L4)
cexx = cex

Where M, C and K are the mass, damping and spring constants respectively, for which I adopted the following values: M = 1, C = 0.02 and K = 0.05. Using the LIFT function I save the lift coefficient values ​​to be able to solve numerically.
When running the program, in some cases the solution diverges, and in others I get a segmentation error that does not allow the simulation to be carried out. I think I am incorrectly modifying the codes for the simulation to function correctly, maybe someone can help me.
On the other hand, I think that the way I define the speed could be wrong, which I do by assigning said value to the ubcy parameter, and I would like to know if that value should represent the initial speed in the direction of the "y" axis or if it should be the instantaneous speed.
I also leave here the values ​​of the relevant parameters that I am using to run the simulation:

! Mesh
nx=361 ! X-direction nodes
ny=325 ! Y-direction nodes
nz=1 ! Z-direction nodes

! Domain
xlx = 20. ! Lx (Size of the box in x-direction)
yly = 18. ! Ly (Size of the box in y-direction)
zlz = 1. ! Lz (Size of the box in z-direction)

! Boundary conditions
nclx1 = 2
nclxn = 2
ncly1 = 0
nclyn = 0
nclz1 = 0
nclzn = 0

! Flow parameters
iin = 1 ! Inflow conditions (1: classic, 2: turbinit)
re = 185. ! nu=1/re (Kinematic Viscosity)
u1 = 1.0 ! u1 (max velocity) (for inflow condition)
u2 = 1.0 ! u2 (min velocity) (for inflow condition)
init_noise = 0.0 ! Turbulence intensity (1=100%) !! Initial condition
inflow_noise = 0.0 ! Turbulence intensity (1=100%) !! Inflow condition

! Time stepping
dt = 0.0005 ! Time step
ifirst = 1 ! First iteration
ilast = 100000 ! Last iteration

iibm=3

!================
&ibmstuff
!================
imove = 1 ! (0: Stationary Objects, 1: Moving Objects)
ubcx = 0 ! Prescribed Velocity (ux)
ubcz = 0 ! Prescribed Velocity (uz)

cex=8. ! when simulating a cylinder x coordinate of the center
cey=9. ! when simulating a cylinder y coordinate of the center
ra=0.5 ! when simulating a cylinder, radius

iforces=1
nvol=1 !Number of volumes for computing force balance

&ForceCVs

xld(1) = 1.0 !X left for volume control
xrd(1) = 15.0 !X right for volume control
yld(1) = 1.0 !Y bottom for volume control
yud(1) = 17.0 !Y top for volume control

Thanks in advance, regards.

[mathrack]

Hello,

To avoid potential issues connected with the coupling of the spring-mass solver, have you tried to perform a simulation with an imposed periodic motion for the cylinder ? Using the natural frequency associated with the spring and the mass ?

We do not provide support for large modifications in the code, but changing the imposed motion is a minor modification and it is easier for us to help.

If we could reproduce the segmentation fault with this simpler case it would be interesting.

EDIT : ubcx and ubcy are the instantaneous velocities, I think this is correct

Regards,
Mathrack

[Manu-CFD]

Hello, thanks for answering, yes, I previously simulated the forced oscillation for a range of frequencies including the natural frequency, and apart from some differences in the lift and drag values ​​compared to those obtained from the literature, I did not have any problems with segmentation and the simulation seemed to be correct.

[mathrack]

Without a simple reproducible example, I am not sure we can help.

As a side remark, please note that the number of nodes in y is 325 (5x5x13), this is probably not the best choice for a periodic direction.

To fix segmentation faults, you can try to compile the code in debug mode and run it on a single CPU with gdb.

[mathrack]

Hello,

I think for some specific cases, the variables ipol and jpol can be used without being defined, see
https://github.com/xcompact3d/Incompact3d/blob/master/src/ibm.f90#L524
and
https://github.com/xcompact3d/Incompact3d/blob/master/src/ibm.f90#L691

Changing ipol to ipolf and jpol to jpolf in the lines above might solve your segmentation fault

EDIT : this will not fix your issue if ipolf and jpolf are incorrect...

[Manu-CFD]

Thank you very much, I'm going to try those modifications, greetings

[mathrack]

Here is a sketch for an imposed periodic motion : https://github.com/mathrack/Incompact3d-1/tree/moving_cylinder

I have not validated it but the code is running on my laptop without segmentation fault.

[Manu-CFD]

Hello, thanks for all the contribution, I was trying to run the simulation with the modifications you proposed, but I get an error related to reading the input file, both when I try to run my VIV code and when I test your forced oscillation code . The error that appears with my code is the following:
Error termination. Backtrace:
At line 221 of file /home/manukr98/NEW/Incompact3d/src/parameters.f90 (unit = 10, file = 'input.i3d')
Fortran runtime error: End of file

And when I run your code it is similar:
Termination error. Backtrace:
At line 113 of file /home/manukr98/NEW/Incompact3d/src/parameters.f90 (unit = 10, file = 'input.i3d')
Fortran runtime error: Cannot match namelist object name cyl_oscil

I checked the definition of the new parameters and the input file but I can't find where the error comes from.

[mathrack]

I am not sure to understand the error you have obtained. Your version of the code is not available on github so I can not look at the content of the provided line.

In the modified version of the code I have suggested, line 113 of parameters.f90 is a comment. This should not trigger some error. If you are running the modified version I have suggested, please make sure you are running the source available on github and provide

• input file (if different from the input files included in examples/)
• the output of the code, including a few lines before the error
• the number of CPUs

Cheers

[Manu-CFD]

Hello, the problem has finally been resolved, it is now working perfectly having downloaded the latest version of the code and using your modifications. Apparently the error I was getting was also related to the position and velocity values ​​of the cylinder diverging due to a poor choice of the parameters I was using when numerically solving the mass-spring-damper equation. Thank you very much for your contributions, greetings

[mathrack]

You will also need to fix the computation of ypraf in genepsi3d. You can use      do j=1,nyraf         ypraf(j) = (j-1)*dyraf       enddo Sep 10, 2024 21:16:53 Manuel Krapf ***@***.***>:

…

Thank you very much, I'm going to try those modifications, greetings — Reply to this email directly, view it on GitHub[#292 (comment)], or unsubscribe[https://github.com/notifications/unsubscribe-auth/ADHFJ623C2TBG6VZ7Y5NTJTZV5ASDAVCNFSM6AAAAABNW43Y4KVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTANRQGY2DMNI]. You are receiving this because you commented. [Tracking image][https://github.com/notifications/beacon/ADHFJ657DB3P65OL4QWVHU3ZV5ASDA5CNFSM6AAAAABNW43Y4KWGG33NNVSW45C7OR4XAZNRIRUXGY3VONZWS33OINXW23LFNZ2KUY3PNVWWK3TUL5UWJTQAUHLYC.gif]