[asking]

Hi,

I am implementing a UDF of a cylinder subjected to a fluid flow in 2D. My UDF basically calculates the forces on the cylinder, uses a numerical scheme (runge-kutta) to determine the velocity and displacement of the cylinder, and updates the mesh accordingly.

The problem I have is that the simulation crashes after ~0.3 seconds. Looking at the force plots, it seems that after a few iterations, the lift force increases to unrealistic high values and then oscillates between positive and negative values at each time step. This produces huge displacement fluctuations and, thus, negative cell volumes. I've tried reducing the timestep, implicit updates every 5 iterations (it helps a little bit), different pressure-velocity solving schemes, but the problem persist.

The integration scheme used comes from: https://www.sciencedirect.com/scienc...41118717306454

They also used FLUENT, so it should work.

Here is the UDF being used:

#include "udf.h"
#include "mem.h"
#define PI 3.141592654
#define zoneID 8
FILE *outRK;
static real yRK = 0.0;
static real vyRK = 0.0;
static real current_time = -1;

DEFINE_CG_MOTION(cylinder_mesh_MC_3,dt,vel,omega,time,dtime)
{
real ctime = RP_Get_Real("flow-time");
real ctimestep = RP_Get_Integer("time-step");
real niter = N_ITER;

if (current_time < ctimestep)
{
current_time = ctimestep;

/*Define variables*/

/*Mesh variables*/
real cg[3],vcg[3];

/*Cylinder variables*/
real diameter = 0.01;
real fn = 14;
real density = 998.2;
real length = 1;
real water_depth = 1;

real mass = density*pow((0.5*diameter),2)*PI*length;
real ad_mass = density*pow((0.5*diameter),2)*PI*water_depth;
real total_mass = mass + ad_mass;
real k = 4*pow((PI*fn),2)*total_mass;
real c = 2 * 0.04 * sqrt(k*total_mass);

/*Force calculation. Force = F_pressure + F_viscous*/
real fy = 0.0;
real fvy = 0.0;
int i;

#if !RP_HOST
Thread *tc,*thread;
Domain *d = Get_Domain(1);
face_t f;
tc = Lookup_Thread(d,zoneID);
thread = DT_THREAD(dt);
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
real NV_VEC(A);

begin_f_loop(f,tc)
{
if (PRINCIPAL_FACE_P(f,tc))
{
fvy = F_STORAGE_R_N3V(f,tc,SV_WALL_SHEAR)[1]*-1;
F_AREA(A,f,tc);

/*Force calculation with a depth of 1m*/
fy += F_P(f,tc)*A[1] + fvy;
}
}
end_f_loop(f,tc)
#endif

#if RP_NODE
fy = PRF_GRSUM1(fy);
#endif

/*Dynamic mesh position*/
#if!RP_HOST
for (i=0;i<3;i++)
{
cg[i]=DT_CG(dt)[i];
vcg[i] = DT_VEL_CG(dt)[i];
}

Message("Position CG: %f \n",cg[1]);
#endif

node_to_host_real_2(fy,cg[1]);

/*Numerical methods*/
/*Runge-kutta 4th order*/
real K1 = (fy*water_depth - c*vyRK - k*cg[1]) / total_mass;
real K2 = (fy*water_depth - c*(vyRK+dtime*0.5*K1) - k*(cg[1]+dtime*0.5*vyRK)) / total_mass;
real K3 = (fy*water_depth - c*(vyRK+dtime*0.5*K2) - k*(cg[1]+dtime*0.5*vyRK+dtime*dtime*K1/8)) / total_mass;
real K4 = (fy*water_depth - c*(vyRK+dtime*K3) - k*(cg[1]+dtime*vyRK+dtime*dtime*K1/2)) / total_mass;

yRK = yRK + vRK*dtime + dtime*dtime*(K1 + K2 + K3 + K4)/6;
vyRK = vyRK + dtime*(K1 + K2 + K3 + K4)/6;

/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vyRK;

/*Save files*/
#if !RP_NODE
if(NULL == (outRK = fopen("dataRK.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outRK,"%.4f %.1f %.3f %.7f %.7f %.7f %.7f \n", ctime,niter, fy , cg[1], yRK, vyRK, vyRK2);
fclose(outRK);
#endif
}
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vyRK;
}

Any ideas/recommendation is greatly appreciated.

[obscureed]

Hi Asking,

I don't have much to suggest -- I've looked at the UDF, and nothing leaps out at me as being wrong. (Oh, one thing: "current_time = ctimestep;" is a strange thing to do, and mixes real and integer values. I doubt you do enough timesteps for this to go wrong. "current_time = ctime;" would make more sense.)

About the timestep: you need a timestep such that (velocity*timestep) is significantly less than mesh size. You can also look at the characteristic time of the spring constant.

You could try Euler integration (that is, the simplest possible scheme: v(t+dt) = v(t) + dt * F/m) as a comparison with the Runge-Kutta. If the timestep is reasonable, the answers should not be very different.

I suspect that the issue may be flow convergence, not the UDF. ("the lift force increases to unrealistic high values") One thing you could try would be to write a different UDF to apply an example of the motion that you are expecting -- for example, a pre-defined sinusoidal motion. See whether the lift force stays reasonable. (Integrate the effect of that lift force outside of Fluent.)

In desperation, you could turn off the damping and spring constant, just to see what happens.

Good luck!
Ed

[asking]

Originally Posted by obscureed

Hi Asking,

I don't have much to suggest -- I've looked at the UDF, and nothing leaps out at me as being wrong. (Oh, one thing: "current_time = ctimestep;" is a strange thing to do, and mixes real and integer values. I doubt you do enough timesteps for this to go wrong. "current_time = ctime;" would make more sense.)

About the timestep: you need a timestep such that (velocity*timestep) is significantly less than mesh size. You can also look at the characteristic time of the spring constant.

You could try Euler integration (that is, the simplest possible scheme: v(t+dt) = v(t) + dt * F/m) as a comparison with the Runge-Kutta. If the timestep is reasonable, the answers should not be very different.

I suspect that the issue may be flow convergence, not the UDF. ("the lift force increases to unrealistic high values") One thing you could try would be to write a different UDF to apply an example of the motion that you are expecting -- for example, a pre-defined sinusoidal motion. See whether the lift force stays reasonable. (Integrate the effect of that lift force outside of Fluent.)

In desperation, you could turn off the damping and spring constant, just to see what happens.

Good luck!
Ed

Thank you for your answer. You gave me ideas to work on (which is great since I have been stuck on this problem for quite some time).

Before trying the RK4 method, I implemented the Newmark-Beta (NB) method. Which is an implicit numerical scheme. I didn't get any numerical error with this method and the cylinder response seemed good. Nevertheless, as I reduced the timestep (everything else equal), the cylinder displacement didn't converge. I attached the displacements and forces I got at different timesteps.

The cylinder response at the highest tested timestep (0.005 s) is closely similar to the paper I am comparing my results with https://link.springer.com/content/pd...017-0010-9.pdf (Page 4, Figure 4). But this motion does not converge as the timestep decreases.

The only differences with this paper are that I scaled the problem changing the cylinder diameter (D=0.01m vs D=1m) and velocity (V=0.588 m/s vs V=58.8m/s). I also changed the mesh distribution (similar to the image attached).

From those results, I thought that maybe it was my implementation of NB. So, that's why I am trying to use RK4.

Some parameters:

D=0.01m
V=0.588 m/s
1st inflation layer = 0.001 mm (yplus = 1)
inner ring size = 0.04m

k-w SST turbulent model
Coupled Implicit second order solver
40 iterations and convergence at 1e-4

I am going to try a few things as you suggested and hopefully improve my results. Thanks!

[eagle_001]

Hi, South Yorkshire !
Thank you very much for your useful UDF, hope you have solved the problem completely. In your another post:
Runge kutta 4th order: Cylinder,
a gentleman named blackmask suggested you to use DEFINE_EXECUTE_AT_END to update the vy. I wonder whether you have tried this strategy. Did it real work? Thank you again and best wishes!

[asking]

Originally Posted by eagle_001

Hi, South Yorkshire !
Thank you very much for your useful UDF, hope you have solved the problem completely. In your another post:
Runge kutta 4th order: Cylinder,
a gentleman named blackmask suggested you to use DEFINE_EXECUTE_AT_END to update the vy. I wonder whether you have tried this strategy. Did it real work? Thank you again and best wishes!

Hi eagle,

I completely forgot about my thread! I had to suspend my numerical modelling since I needed to run some experiments. I forgot all the problems I had (I just remember they where many! haha) and how I solved them, unfortunately. Here is the UDF, it has the NB and RK4 schemes

#include "udf.h"
#include "mem.h"
#define PI 3.141592654
#define zoneID 11
FILE *outNB,*outRK;
static real y = 0.0;
static real yRK = 0.0;
static real dy = 0.0;
static real vy = 0.0;
static real vyRK = 0.0;
static real vyRK2 = 0.0;
static real ay = 0.0;
static real current_time = -1;

DEFINE_CG_MOTION(cylinder_mesh_MC_3,dt,vel,omega,time,dtime)
{
real ctime = RP_Get_Real("flow-time");
real ctimestep = RP_Get_Integer("time-step");
real niter = N_ITER;

if (current_time < ctimestep)
{
current_time = ctimestep;

/*Define variables*/

/*Mesh variables*/
real cg[3],vcg[3];

/*Cylinder variables*/
real diameter = 0.02;
real fn = 4.4;
real density = 998.2;
real length = 1;
real water_depth = 1;
real mass_ratio = 0.85;
real damping_ratio = 0.0032;

real mass = mass_ratio*density*pow((0.5*diameter),2)*PI*length;
real ad_mass = mass*(1/3); /*density*pow((0.5*diameter),2)*PI*water_depth;*/
real total_mass = mass + ad_mass;
real k = 4*pow((PI*fn),2)*total_mass;
real c = 2 * damping_ratio * sqrt(k*total_mass);

/*Force calculation. Force = F_pressure + F_viscous*/
real fy = 0.0;
real fvy = 0.0;
int i;

#if !RP_HOST
Thread *tc,*thread;
Domain *d = Get_Domain(1);
face_t f;
tc = Lookup_Thread(d,zoneID);
thread = DT_THREAD(dt);
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
real NV_VEC(A);

begin_f_loop(f,tc)
{
if (PRINCIPAL_FACE_P(f,tc))
{
fvy = F_STORAGE_R_N3V(f,tc,SV_WALL_SHEAR)[1]*-1;
F_AREA(A,f,tc);

/*Force calculation with a depth of 1m*/
fy += F_P(f,tc)*A[1] + fvy;
}
}
end_f_loop(f,tc)
#endif

#if RP_NODE
fy = PRF_GRSUM1(fy);
#endif

/*Dynamic mesh position*/
#if!RP_HOST
for (i=0;i<3;i++)
{
cg[i]=DT_CG(dt)[i];
vcg[i] = DT_VEL_CG(dt)[i];
}

Message("Position CG: %f \n",cg[1]);
#endif

node_to_host_real_2(fy,cg[1]);

/*Numerical methods*/
/*Numark-beta*/
real beta = 0.25;
real gamma = 0.5;
real term0 = (1/(beta*dtime*dtime))*(mass+ad_mass) + (gamma/(beta*dtime))*c;
real term1 = (1/(beta*dtime))*(mass+ad_mass) + ((gamma/beta)-1)*c;
real term2 = ((1/(2*beta))-1)*(mass+ad_mass) + dtime*((gamma/(2*beta))-1)*c;

real Keff = k + term0;
real Reff = fy*water_depth + term0*cg[1] + term1*vy + term2*ay;
Message("Velocity: %f \n",vy);
dy = Reff/Keff - cg[1];
y += dy;

real vprev = vy;
vy = (gamma/(beta*dtime))*dy + (1-(gamma/beta))*vy + dtime*(1-(gamma/(2*beta)))*ay;
ay = (1/(beta*dtime*dtime))*dy - (1/(beta*dtime))*vprev - ((1/(2*beta))-1)*ay;

/*Runge-kutta 4th order*/
real K1 = (fy*water_depth - c*vyRK - k*yRK) / total_mass;
real K2 = (fy*water_depth - c*(vyRK+dtime*0.5*K1) - k*(yRK+dtime*0.5*vyRK)) / total_mass;
real K3 = (fy*water_depth - c*(vyRK+dtime*0.5*K2) - k*(yRK+dtime*0.5*vyRK+dtime*dtime*K1/8)) / total_mass;
real K4 = (fy*water_depth - c*(vyRK+dtime*K3) - k*(yRK+dtime*vyRK+dtime*dtime*K1/2)) / total_mass;

yRK = yRK + vyRK*dtime + dtime*dtime*(K1 + K2 + K3 + K4)/6;
vyRK = vyRK + dtime*(K1 + K2 + K3 + K4)/6;
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vy;

/*Save files*/
#if !RP_NODE
/*Message ("Force = %f, pos = %f, vel = %f, acc = %f\n", fy, cg[1], y, vy);*/

if(NULL == (outNB = fopen("dataNB.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outNB,"%.4f %.1f %.3f %.7f %.7f %.7f \n", ctime,niter, fy , cg[1], y, vy);
fclose(outNB);

if(NULL == (outRK = fopen("dataRK.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outRK,"%.4f %.1f %.3f %.7f %.7f %.7f \n", ctime,niter, fy , cg[1], yRK, vyRK);
fclose(outRK);
#endif
}
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vy;
}

You need to change the /*Cylinder variables*/ part and the #define zoneID 11. Be careful about how the parameters are defined! (mass of the cylinder, for example). I need to change those in the future to avoid any confusion.

It works with the dynamic mesh diffusion method (before it only worked with the Spring method). Currently it is for 1DOF only (transverse motion) but it can be easily adapted for 2DOF (if you assume that the streamwise and crossflow motions are independent).

The last problem I had was related to convergence. The displacement/force didn't stabilise even after 30+ seconds of simulations. To speed this process, I usually run the first 0.1 s with a small transverse inlet velocity to speed up the vortex shedding formation. Then, I remove that transverse velocity and keep running the whole thing.

My goals now are:
- check with another paper that I am getting decent results.
- Adapting it to 2DOF.
- Implement the strip method. Maybe calculating 4 or 5 strips (2d simulations) along the span of a cylinder. It's just an idea though. I have no Idea how to do that.

Hope this helps and sorry for not updating my own thread!
Regards,

[eagle_001]

Hi, South Yorkshire!
Thanks a lot for your updated UDF and detailed explanations. I will try this immediately.
If you want to calculate the 3D case of cylinder, how about to use the WORKBENCH system coupling module? It will be much more easier to perform this in WORKBENCH, but maybe much more computation time consuming.
We also use the strip theory in aeronautics for the high-aspect-ratio wing, and get acceptable results. We often treat it as a low fidelity model, because it ignores the 3D effect of the flow field.
Best wishes!

[asking]

Originally Posted by eagle_001

Hi, South Yorkshire!
Thanks a lot for your updated UDF and detailed explanations. I will try this immediately.
If you want to calculate the 3D case of cylinder, how about to use the WORKBENCH system coupling module? It will be much more easier to perform this in WORKBENCH, but maybe much more computation time consuming.
We also use the strip theory in aeronautics for the high-aspect-ratio wing, and get acceptable results. We often treat it as a low fidelity model, because it ignores the 3D effect of the flow field.
Best wishes!

No problem! Glad to help.
What you said is true. A full 3d model would be ideal. One issue I have is that I want to try a pinned cylinder, where its displacement varies along the span of the cylinder. This means that the mesh will have to deform against the cylinder rotation along its fixed end and I have no idea how to consider that with a UDF.
That's why I want to move to the strip method. Do you know a tutorial about how to implement this method on ANSYS?

Regards,

[arunraj]

Hello guys,

Could you please elaborate in detail how to set up dynamic mesh method for free vibration problem with inflation layer. It would be a great help for many researchers in this forum.

[David070]

Hi Asking,

I am currently doing a similar problem as yours.

The problem I model is an object that is attached with a torsional spring, and the other end of the spring is under a prescribed motion. So the object will move under the interaction of the fluid and the elastic reaction force. I use the DEFINE_ZONE_MOTION to update the velocity of the object. (I use the algorithm of v(t+dt) = v(t) + dt * F/m). The problem is that this macro is called every time step, and the fluid force may not converge or let's say may be higher than reality, then the object will move fast because of the fluid force, then at the second timestep the fast speed will cause an even large fluid force, and then the calculation would soon diverge.

Do you have any suggestions for this?

I want to modify the program like this: the position of the object will update at every time step, but the velocity of the dynamic motion(fsi) is updated at every iteration. Once the dynamic motion is converged, then it enters the next time step. By doing this, the dynamic motion of the object will easily converge at every timestep. Then the calculation should be work.

So, do you know if there is a macro that can update the velocity every iteration if we use this way to model the problem?

Thank you very much!

Originally Posted by asking

Hi eagle,

I completely forgot about my thread! I had to suspend my numerical modelling since I needed to run some experiments. I forgot all the problems I had (I just remember they where many! haha) and how I solved them, unfortunately. Here is the UDF, it has the NB and RK4 schemes

#include "udf.h"
#include "mem.h"
#define PI 3.141592654
#define zoneID 11
FILE *outNB,*outRK;
static real y = 0.0;
static real yRK = 0.0;
static real dy = 0.0;
static real vy = 0.0;
static real vyRK = 0.0;
static real vyRK2 = 0.0;
static real ay = 0.0;
static real current_time = -1;

DEFINE_CG_MOTION(cylinder_mesh_MC_3,dt,vel,omega,time,dtime)
{
real ctime = RP_Get_Real("flow-time");
real ctimestep = RP_Get_Integer("time-step");
real niter = N_ITER;

if (current_time < ctimestep)
{
current_time = ctimestep;

/*Define variables*/

/*Mesh variables*/
real cg[3],vcg[3];

/*Cylinder variables*/
real diameter = 0.02;
real fn = 4.4;
real density = 998.2;
real length = 1;
real water_depth = 1;
real mass_ratio = 0.85;
real damping_ratio = 0.0032;

real mass = mass_ratio*density*pow((0.5*diameter),2)*PI*length;
real ad_mass = mass*(1/3); /*density*pow((0.5*diameter),2)*PI*water_depth;*/
real total_mass = mass + ad_mass;
real k = 4*pow((PI*fn),2)*total_mass;
real c = 2 * damping_ratio * sqrt(k*total_mass);

/*Force calculation. Force = F_pressure + F_viscous*/
real fy = 0.0;
real fvy = 0.0;
int i;

#if !RP_HOST
Thread *tc,*thread;
Domain *d = Get_Domain(1);
face_t f;
tc = Lookup_Thread(d,zoneID);
thread = DT_THREAD(dt);
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
real NV_VEC(A);

begin_f_loop(f,tc)
{
if (PRINCIPAL_FACE_P(f,tc))
{
fvy = F_STORAGE_R_N3V(f,tc,SV_WALL_SHEAR)[1]*-1;
F_AREA(A,f,tc);

/*Force calculation with a depth of 1m*/
fy += F_P(f,tc)*A[1] + fvy;
}
}
end_f_loop(f,tc)
#endif

#if RP_NODE
fy = PRF_GRSUM1(fy);
#endif

/*Dynamic mesh position*/
#if!RP_HOST
for (i=0;i<3;i++)
{
cg[i]=DT_CG(dt)[i];
vcg[i] = DT_VEL_CG(dt)[i];
}

Message("Position CG: %f \n",cg[1]);
#endif

node_to_host_real_2(fy,cg[1]);

/*Numerical methods*/
/*Numark-beta*/
real beta = 0.25;
real gamma = 0.5;
real term0 = (1/(beta*dtime*dtime))*(mass+ad_mass) + (gamma/(beta*dtime))*c;
real term1 = (1/(beta*dtime))*(mass+ad_mass) + ((gamma/beta)-1)*c;
real term2 = ((1/(2*beta))-1)*(mass+ad_mass) + dtime*((gamma/(2*beta))-1)*c;

real Keff = k + term0;
real Reff = fy*water_depth + term0*cg[1] + term1*vy + term2*ay;
Message("Velocity: %f \n",vy);
dy = Reff/Keff - cg[1];
y += dy;

real vprev = vy;
vy = (gamma/(beta*dtime))*dy + (1-(gamma/beta))*vy + dtime*(1-(gamma/(2*beta)))*ay;
ay = (1/(beta*dtime*dtime))*dy - (1/(beta*dtime))*vprev - ((1/(2*beta))-1)*ay;

/*Runge-kutta 4th order*/
real K1 = (fy*water_depth - c*vyRK - k*yRK) / total_mass;
real K2 = (fy*water_depth - c*(vyRK+dtime*0.5*K1) - k*(yRK+dtime*0.5*vyRK)) / total_mass;
real K3 = (fy*water_depth - c*(vyRK+dtime*0.5*K2) - k*(yRK+dtime*0.5*vyRK+dtime*dtime*K1/8)) / total_mass;
real K4 = (fy*water_depth - c*(vyRK+dtime*K3) - k*(yRK+dtime*vyRK+dtime*dtime*K1/2)) / total_mass;

yRK = yRK + vyRK*dtime + dtime*dtime*(K1 + K2 + K3 + K4)/6;
vyRK = vyRK + dtime*(K1 + K2 + K3 + K4)/6;
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vy;

/*Save files*/
#if !RP_NODE
/*Message ("Force = %f, pos = %f, vel = %f, acc = %f\n", fy, cg[1], y, vy);*/

if(NULL == (outNB = fopen("dataNB.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outNB,"%.4f %.1f %.3f %.7f %.7f %.7f \n", ctime,niter, fy , cg[1], y, vy);
fclose(outNB);

if(NULL == (outRK = fopen("dataRK.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outRK,"%.4f %.1f %.3f %.7f %.7f %.7f \n", ctime,niter, fy , cg[1], yRK, vyRK);
fclose(outRK);
#endif
}
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vy;
}

You need to change the /*Cylinder variables*/ part and the #define zoneID 11. Be careful about how the parameters are defined! (mass of the cylinder, for example). I need to change those in the future to avoid any confusion.

It works with the dynamic mesh diffusion method (before it only worked with the Spring method). Currently it is for 1DOF only (transverse motion) but it can be easily adapted for 2DOF (if you assume that the streamwise and crossflow motions are independent).

The last problem I had was related to convergence. The displacement/force didn't stabilise even after 30+ seconds of simulations. To speed this process, I usually run the first 0.1 s with a small transverse inlet velocity to speed up the vortex shedding formation. Then, I remove that transverse velocity and keep running the whole thing.

My goals now are:
- check with another paper that I am getting decent results.
- Adapting it to 2DOF.
- Implement the strip method. Maybe calculating 4 or 5 strips (2d simulations) along the span of a cylinder. It's just an idea though. I have no Idea how to do that.

Hope this helps and sorry for not updating my own thread!
Regards,

[AlexanderZ]

macro DEFINE_CG_MOTION is executed on each timestep
here

vel[0] = 0.0;
vel[1] = vy;

velocity was updated

isn't it what you are asking for? so what's your question?

[David070]

Hi Alexander,

Thank you for your reply. Actually, I want to update the velocity every iteration, but the location of the object updates every timestep.
This page describes my question:
Update every iteration!!!!!

Originally Posted by AlexanderZ

macro DEFINE_CG_MOTION is executed on each timestep
here

vel[0] = 0.0;
vel[1] = vy;

velocity was updated

isn't it what you are asking for? so what's your question?