圆柱结构在水流冲击下被动摆动的UDF，还请大家帮忙找下问题，谢谢！
-
我写了一个udf，想描述水流冲击下圆柱的摆动。但是现在稳定是稍微一刺激越摆越大，直接就发散了，我感觉应该是摆动方程有问题，还请大佬指点一下。这是我的udf：
#include "udf.h" #define zoneID1 2 //浮子的编号 #define IXX 0.003 //转动惯量 #define mg 0.98 //重力 #define r 0.25 //回转半径 #define sp 1 //弹簧回转力矩 real v_body1; //转速 DEFINE_EXECUTE_AT_END(exe_end) { int i; real tm,dtm; FILE *fd1; real f1[3],m1[3]; real x1[3]; Domain *domain=Get_Domain(1); Thread *tf1; real dx = 0.0;//x方向转角增量 real x = 0.0;//x方向转角 real ve = 0.0;//x方向角速度 real a = 0.0;//x方向角加速度 real ve_before; real a_before; tf1=Lookup_Thread(domain,zoneID1); if (!Data_Valid_P()) return; for(i=0;i<=2;i++) //三个方向（x y x方向） { f1[i]=0; m1[i]=0; x1[i]=0; } Compute_Force_And_Moment(domain,tf1,x1,f1,m1,TRUE);//每个柱上的受力、力矩 ve_before = ve, a_before = a; dtm=RP_Get_Real("physical-time-step");//获取当前时间步长 tm=RP_Get_Real("flow-time");//获取当前计算时间 a = (m1[0]+mg*sin(x)*r-sp*x ) / IXX;//摆动方程，求解转速加速度，第一项是输出转矩，第二项是重力，第三项是弹簧回转力 ve = ve + (a + a_before)*dtm / 2; //转速 dx = (ve + ve_before)*dtm / 2; //移动角度增量 x = x + dx; //移动角度 v_body1=ve; //输结果 fd1 = fopen("SwingResult.txt", "a"); fprintf(fd1, "%.5f %.5f %.6f %.6f %.6f %.6f\n", dtm, tm, a, ve, dx, x); fclose(fd1); } DEFINE_ZONE_MOTION(shake1234, omega, axis, origin, velocity, time, dtime) { origin[0] = 0; origin[1] = 0; origin[2] = 0; axis[0] = 1; axis[1] = 0; axis[2] = 0; *omega =v_body1; return; }本人刚开始做这个，我感觉是48到51行摆动方程有问题。
目前我考虑的是，圆柱摆动过程，受到1）水流冲击转矩，2）自身重力和3）弹簧回摆力，即48到51的代码。
请问是否应该将连续的方程，离散成模拟过程中的方程，这个应该怎么做，能否提供一些学习资料或者参考？感谢指点！
-
目前还在继续做，还是实现不了。我将运动方程中加入了浮力，并将方程通过4阶龙格库塔进行了展开，但是目前还是稳定不了。
目前我想的平衡方程是：
角加速度=（来流转矩+重力转矩-浮力转矩-弹簧转矩-介质阻尼）/转动惯量
见程序51行。我通过调整转动惯量进行测试，但是发现有一个问题：
1）转动惯量为0.015的时候，模型直接运动不起来，摆动的角位移一直为0，角速度倒是会在很小的范围内振荡；
2）转动惯量调整到0.01的时候，会摆动越来越大一直到发散。在这个过程中调整弹簧刚度和介质阻尼，都无法阻止这个发散，只是会减慢发散速度。这是目前的程序，还请大佬指点！
#include "udf.h" #define zoneID1 2 // 浮子的编号 #define IXX 0.003 // 转动惯量 #define mg -0.01 // 重力减去浮力=0.1-0.09 #define r 0.25 // 回转半径 #define sp 1 // 弹簧回转力矩 #define D 0.1 // 阻尼系数 real v_body1 = 0.0; // 转速 real theta = 0.0; // 初始角位移 real omega = 0.0; // 初始角速度 real a=0.1; // 角加速度 DEFINE_EXECUTE_AT_END(exe_end) { int i; real tm, dtm; FILE *fd1; real f1[3], m1[3]; real x1[3]; Domain *domain = Get_Domain(1); Thread *tf1; real omegak1, omegak2, omegak3, omegak4; real thetak1, thetak2, thetak3, thetak4; tf1 = Lookup_Thread(domain, zoneID1); if (!Data_Valid_P()) return; for (i = 0; i < 3; i++) // 三个方向（x y z方向） { f1[i] = 0; m1[i] = 0; x1[i] = 0; } Compute_Force_And_Moment(domain, tf1, x1, f1, m1, TRUE); // 每个柱上的受力、力矩 dtm = RP_Get_Real("physical-time-step"); // 获取当前时间步长 tm = RP_Get_Real("flow-time"); // 获取当前计算时间 // RK4方法，求解角加速度和更新角速度和角位移 // 第一步的角加速度 thetak1 = omega; omegak1 = (m1[0] +mg * r * sin(theta) - sp * theta - D * omega) / IXX; // 第二步的角加速度 thetak2 = omega + 0.5 * dtm * omegak1; omegak2 = (m1[0] + mg * r * sin(theta + 0.5 * dtm * thetak1) - sp * (theta + 0.5 * dtm * thetak1) - D * (omega + 0.5 * dtm * omegak1)) / IXX; // 第三步的角加速度 thetak3 = omega + 0.5 * dtm * omegak2; omegak3 = (m1[0] + mg * r * sin(theta + 0.5 * dtm * thetak2) - sp * (theta + 0.5 * dtm * thetak2) - D * (omega + 0.5 * dtm * omegak2)) / IXX; // 第四步的角加速度 thetak4 = omega + dtm * omegak3; omegak4 = (m1[0] + mg * r * sin(theta + dtm * thetak3) - sp * (theta + dtm * thetak3) - D * (omega + dtm * omegak3)) / IXX; // 使用RK4方法更新角速度和角位移 omega = omega + (dtm / 6) * (omegak1 + 2 * omegak2 + 2 * omegak3 + omegak4); theta = theta + (dtm / 6) * (thetak1 + 2 * thetak2 + 2 * thetak3 + thetak4); v_body1 = omega; // 输出结果 fd1 = fopen("SwingResult.txt", "a"); fprintf(fd1, "%.5f %.5f %.6f %.6f \n", dtm, tm, omega, theta); fclose(fd1); } DEFINE_ZONE_MOTION(shake1234, omega, axis, origin, velocity, time, dtime) { origin[0] = 0; origin[1] = 0; origin[2] = 0; axis[0] = 1; axis[1] = 0; axis[2] = 0; *omega = v_body1; return; }