fluent重叠网格和UDF模拟串列三圆柱涡激振动的计算发散

[欣晴]

采用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;
}

运动速度

运动位移

压力系数

动网格设置