采用overset重叠网格和UDF函数模拟串列三圆柱的涡激振动，其中UDF分别尝试了Newmark-Beta方法和4阶Runge Kutta法获取圆柱的振动响应，通过DEFINE_CG_MOTION宏赋予三个圆柱及component cells的运动速度，算例的雷诺数约200，采用k-omega sst模型，考虑水作为来流介质，计算过程中尝试了时间步长从1.0e-3缩小到1.0e-5等多个量级，但是求解过程总是出现升力和升力矩突然骤增，继而导致圆柱运动速度过大，最终计算发散。
请各位大佬帮忙看看是哪里出问题？
具体的UDF和部分设置如下：
#include "udf.h"
#include "sg_mem.h"
#include "dynamesh_tools.h"
#define PI 3.141592654
#define zoneID_1 4
#define zoneID_2 16
#define zoneID_3 20
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 = 5;
static real y2 = 0.0;
static real y2RK = 0.0;
static real dy2 = 0.0;
static real vy2 = 0.0;
static real vy2RK = 0.0;
static real vy2RK2 = 0.0;
static real ay2 = 0.0;
static real current_time2 = 5;
static real y3 = 0.0;
static real y3RK = 0.0;
static real dy3 = 0.0;
static real vy3 = 0.0;
static real vy3RK = 0.0;
static real vy3RK2 = 0.0;
static real ay3 = 0.0;
static real current_time3 = 5;
DEFINE_CG_MOTION(cylinder_1,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.063;
real fn = 1.0892;
real density = 998.2;
real length = 1;
real water_depth = 1;
real mass_ratio = 0.3937;
real damping_ratio = 0.01;
real mass = mass_ratio*density*pow((0.5*diameter),2)*PI*length;
real ad_mass = mass*(0); /*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_1);
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; /*“*-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/4)) / 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] = vyRK;
/*Save files*/
#if !RP_NODE
/*Message ("Force = %f, pos = %f, vel = %f, acc = %f\n", fy, cg[1], y, vy);*/
if(NULL == (outNB = fopen("dataNB1.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outNB,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y, vy);
fclose(outNB);
if(NULL == (outRK = fopen("dataRK1.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outRK,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], yRK, vyRK);
fclose(outRK);
#endif
}
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vyRK;
}
DEFINE_CG_MOTION(cylinder_1_frontgrid_1,dt,vel,omega,time,dtime)
{
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
vel[0]=0.0;
vel[1]=vyRK;
}
DEFINE_CG_MOTION(cylinder_1_overset_2,dt,vel,omega,time,dtime)
{
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
vel[0]=0.0;
vel[1]=vyRK;
}
DEFINE_ZONE_MOTION(cylinder_1_zone,omega,axis,origin,velocity,time,dtime)
{
N3V_D(velocity, =, 0, 0, 0);
N3V_S(origin, =, -0.32);
N3V_D(axis, =, 0.0, 0.0, 1.0);
velocity[1]=vyRK;
}
DEFINE_CG_MOTION(cylinder_2,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_time2 < ctimestep)
{
current_time2 = ctimestep;
/*Define variables*/
/*Mesh variables*/
real cg[3],vcg[3];
/*Cylinder variables*/
real diameter = 0.063;
real fn = 1.0892;
real density = 998.2;
real length = 1;
real water_depth = 1;
real mass_ratio = 0.3937;
real damping_ratio = 0.01;
real mass = mass_ratio*density*pow((0.5*diameter),2)*PI*length;
real ad_mass = mass*(0); /*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_2);
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; /*“*-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*vy2 + term2*ay2;
Message("Velocity: %f \n",vy2);
dy2 = Reff/Keff - cg[1];
y2 += dy2;
real vprev = vy2;
vy2 = (gamma/(beta*dtime))*dy2 + (1-(gamma/beta))*vy2 + dtime*(1-(gamma/(2*beta)))*ay2;
ay2 = (1/(beta*dtime*dtime))*dy2 - (1/(beta*dtime))*vprev - ((1/(2*beta))-1)*ay2;
/*Runge-kutta 4th order*/
real K1 = (fy*water_depth - c*vy2RK - k*yRK) / total_mass;
real K2 = (fy*water_depth - c*(vy2RK+dtime*0.5*K1) - k*(y2RK+dtime*0.5*vy2RK)) / total_mass;
real K3 = (fy*water_depth - c*(vy2RK+dtime*0.5*K2) - k*(y2RK+dtime*0.5*vy2RK+dtime*dtime*K1/4)) / total_mass;
real K4 = (fy*water_depth - c*(vy2RK+dtime*K3) - k*(y2RK+dtime*vy2RK+dtime*dtime*K1/2)) / total_mass;
y2RK = y2RK + vy2RK*dtime + dtime*dtime*(K1 + K2 + K3 + K4)/6;
vy2RK = vy2RK + dtime*(K1 + K2 + K3 + K4)/6;
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vy2RK;
/*Save files*/
#if !RP_NODE
/*Message ("Force = %f, pos = %f, vel = %f, acc = %f\n", fy, cg[1], y, vy);*/
if(NULL == (outNB = fopen("dataNB2.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outNB,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y2, vy2);
fclose(outNB);
if(NULL == (outRK = fopen("dataRK2.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outRK,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y2RK, vy2RK);
fclose(outRK);
#endif
}
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vy2RK;
}
DEFINE_CG_MOTION(cylinder_2_frontgrid_1,dt,vel,omega,time,dtime)
{
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
vel[0]=0.0;
vel[1]=vy2RK;
}
DEFINE_CG_MOTION(cylinder_2_overset_2,dt,vel,omega,time,dtime)
{
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
vel[0]=0.0;
vel[1]=vy2RK;
}
DEFINE_ZONE_MOTION(cylinder_2_zone,omega,axis,origin,velocity,time,dtime)
{
N3V_D(velocity, =, 0, 0, 0);
N3V_S(origin, =, 0.0);
N3V_D(axis, =, 0.0, 0.0, 1.0);
velocity[1]=vy2RK;
}
DEFINE_CG_MOTION(cylinder_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_time3 < ctimestep)
{
current_time3 = ctimestep;
/*Define variables*/
/*Mesh variables*/
real cg[3],vcg[3];
/*Cylinder variables*/
real diameter = 0.063;
real fn = 1.0892;
real density = 998.2;
real length = 1;
real water_depth = 1;
real mass_ratio = 0.3937;
real damping_ratio = 0.01;
real mass = mass_ratio*density*pow((0.5*diameter),2)*PI*length;
real ad_mass = mass*(0); /*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_3);
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; /*“*-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*vy3 + term2*ay3;
Message("Velocity: %f \n",vy3);
dy3 = Reff/Keff - cg[1];
y3 += dy3;
real vprev = vy3;
vy3 = (gamma/(beta*dtime))*dy3 + (1-(gamma/beta))*vy3 + dtime*(1-(gamma/(2*beta)))*ay3;
ay3 = (1/(beta*dtime*dtime))*dy3 - (1/(beta*dtime))*vprev - ((1/(2*beta))-1)*ay3;
/*Runge-kutta 4th order*/
real K1 = (fy*water_depth - c*vy3RK - k*y3RK) / total_mass;
real K2 = (fy*water_depth - c*(vy3RK+dtime*0.5*K1) - k*(y3RK+dtime*0.5*vy3RK)) / total_mass;
real K3 = (fy*water_depth - c*(vy3RK+dtime*0.5*K2) - k*(y3RK+dtime*0.5*vy3RK+dtime*dtime*K1/4)) / total_mass;
real K4 = (fy*water_depth - c*(vy3RK+dtime*K3) - k*(y3RK+dtime*vy3RK+dtime*dtime*K1/2)) / total_mass;
y3RK = y3RK + vy3RK*dtime + dtime*dtime*(K1 + K2 + K3 + K4)/6;
vy3RK = vy3RK + dtime*(K1 + K2 + K3 + K4)/6;
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vy3RK;
/*Save files*/
#if !RP_NODE
/*Message ("Force = %f, pos = %f, vel = %f, acc = %f\n", fy, cg[1], y, vy);*/
if(NULL == (outNB = fopen("dataNB3.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outNB,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y3, vy3);
fclose(outNB);
if(NULL == (outRK = fopen("dataRK3.txt","a")))
{
Error("Could not open file for append!\n");
}
fprintf(outRK,"%16.4e %12.1f %16.3e %16.7f %16.7f %16.7f \n", ctime,niter, fy , cg[1], y3RK, vy3RK);
fclose(outRK);
#endif
}
/*Transfer result to the dynamic mesh*/
vel[0] = 0.0;
vel[1] = vy3RK;
}
DEFINE_CG_MOTION(cylinder_3_frontgrid_1,dt,vel,omega,time,dtime)
{
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
vel[0]=0.0;
vel[1]=vy3RK;
}
DEFINE_CG_MOTION(cylinder_3_overset_2,dt,vel,omega,time,dtime)
{
NV_S(vel, =, 0.0);
NV_S(omega, =, 0.0);
vel[0]=0.0;
vel[1]=vy3RK;
}
DEFINE_ZONE_MOTION(cylinder_3_zone,omega,axis,origin,velocity,time,dtime)
{
N3V_D(velocity, =, 0, 0, 0);
N3V_S(origin, =, 0.32);
N3V_D(axis, =, 0.0, 0.0, 1.0);
velocity[1]=vy3RK;
}
运动速度
运动位移
压力系数
动网格设置