[asking]

Hi,

I wrote a UDF to calculate the forces of a 2D cylinder subjected to an incoming flow, determine the cylinder's response (acceleration, velocity, and displacement) and prescribe a rigid motion on the cylinder.

The response of the cylinder is determined by solving the dynamic equation:

F = my'' + cy' + ky

The Newmark-beta numerical method seems to work but I want to compare it with another method called runge-kutta 4th order. Considering the same modelling conditions, this numerical method crashes after ~40 iterations due to extremely high displacement values.

I can't figure it out where is the problem. Maybe something related to C, where I am only starting to use it.

Here is my implementation of the RK4 method:

/*Runge Kutta 4th order*/

#if!RP_HOST
func1 = vy;
func2 = (fy*water_depth - c*vy - k*y) / total_mass;

k0 = dtime * func1;
m0 = dtime * func2;

k1 = dtime * (func1 + m0*0.5);
m1 = dtime * (func2 - (c/total_mass)*0.5*m0 - (k/total_mass)*0.5*k0);

k2 = dtime * (func1 + m1*0.5);
m2 = dtime * (func2 - (c/total_mass)*0.5*m1 - (k/total_mass)*0.5*k1);

k3 = dtime * (func1 + m2);
m3 = dtime * (func2 - (c/total_mass)*m2 - (k/total_mass)*k2);

y = y + (k0 + 2*k1 + 2*k2 + k3)/6;
vy = vy + (m0 + 2*m1 + 2*m2 + m3)/6;

func1 = func2 = 0.0;
k0 = k1 = k2 = k3 = 0.0;
m0 = m1 = m2 = m3 = 0.0;
#endif

node_to_host_real_4(fy,cg[1],y,vy);

And here is the full code:

#include "udf.h"
#include "mem.h"
#define PI 3.141592654
#define zoneID 8
FILE *out;
static real y = 0.0;
static real dy = 0.0;
static real vy = 0.0;
static real ay = 0.0;
static real current_time = -1;

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

if (current_time < ctimestep)
{
current_time = ctimestep;

real k0,k1,k2,k3;
real m0,m1,m2,m3;
real func1,func2;

/*Define variables*/
real cg[3], vcg[3];
real fy,fvy;
real Keff,Reff,vprev;
real time1, niter;
real cg_pos;

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

real mass = density*pow((0.5*diameter),2)*PI*length;
real ad_mass = 1000*PI*pow((0.5*diameter),2)*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);

fy = 0.0;
fvy = 0.0;
int i;

real ctime = RP_Get_Real("flow-time");
real ctimestep = RP_Get_Integer("time-step");

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

/*Force calculation*/
#if !RP_HOST
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_HOST
for (i=0;i<3;i++)
{
cg[i]=DT_CG(dt)[i];
vcg[i] = DT_VEL_CG(dt)[i];
}
#endif

#if RP_NODE
fy = PRF_GRSUM1(fy);
/*cg_pos = PRF_GRSUM1(cg[1]);*/
#endif

/*Runge Kutta 4th order*/

#if!RP_HOST
func1 = vy;
func2 = (fy*water_depth - c*vy - k*y) / total_mass;

k0 = dtime * func1;
m0 = dtime * func2;

k1 = dtime * (func1 + m0*0.5);
m1 = dtime * (func2 - (c/total_mass)*0.5*m0 - (k/total_mass)*0.5*k0);

k2 = dtime * (func1 + m1*0.5);
m2 = dtime * (func2 - (c/total_mass)*0.5*m1 - (k/total_mass)*0.5*k1);

k3 = dtime * (func1 + m2);
m3 = dtime * (func2 - (c/total_mass)*m2 - (k/total_mass)*k2);

y = y + (k0 + 2*k1 + 2*k2 + k3)/6;
vy = vy + (m0 + 2*m1 + 2*m2 + m3)/6;

func1 = func2 = 0.0;
k0 = k1 = k2 = k3 = 0.0;
m0 = m1 = m2 = m3 = 0.0;
#endif

node_to_host_real_4(fy,cg[1],y,vy);
vel[0] = 0.0;
vel[1] = vy;
}

vel[1] = vy;
}

It is not a force calculation problem since the values are well calculated when the Newmark-beta method is used.

Can anyone help me with my code?

Thank you.

[blackmask]

Is there any particular reason that you exclude pressure contribution from the force calculation? Is time step the same for the Newmark-beta and the Runge-Kutta scheme? The latter might need a smaller time step.

[asking]

Originally Posted by blackmask

Is there any particular reason that you exclude pressure contribution from the force calculation? Is time step the same for the Newmark-beta and the Runge-Kutta scheme? The latter might need a smaller time step.

I believe I am considering the pressure contribution here

F_P(f,tc)*A[1]

Before this, I considered a fixed cylinder and obtained the same force values from my UDF and Fluent (report - force in Y direction). The pressure reference value is zero.

Yeah, unfortunately, I tried with the same and smaller time step values for the Runge-Kutta scheme with the same results. The weird thing is that I implemented this scheme in matlab at first and works fine.

[blackmask]

You are right. I overlooked that line completely.



I think your implementation is right. I suggest that you can reduce the amplitude of oscillatory motion by setting the initial value of position to the stationary point at the first time step, i.e., at first time step.



Also, you can monitor the values of (y,y',y'') with R-K scheme during the Newmark-beta simulation, and check out where the difference come from.

[asking]

Originally Posted by blackmask

You are right. I overlooked that line completely.



I think your implementation is right. I suggest that you can reduce the amplitude of oscillatory motion by setting the initial value of position to the stationary point at the first time step, i.e., at first time step.



Also, you can monitor the values of (y,y',y'') with R-K scheme during the Newmark-beta simulation, and check out where the difference come from.

Thank you Blackmask.

After implementing your suggestion I confirmed that both numerical schemes are correctly implemented. Both achieve similar displacements.

I think the main problem with my UDF is that it is called at the first iteration of a given time step. I confirm this adding N_ITER inside the if loop. This means that the forces are being calculated with "rough" estimations of the flow. This means that the cylinder will move as a function of a force that is not converged.

Is it possible to possible to call the UDF at the last iterations at each time step?

Regards,

[blackmask]

Yes, you can update "vy" in

DEFINE_EXECUTE_AT_END

[FDE]

Hi,
I need a UDF code to solve the ODE equation (for example: y'=y(y-1), y(0)=...) in FLUENT.
I know that the ADJUST macro must be used. But, I have some problems in writing this code.
Thanks.

[AlexanderZ]

show your code