目前还在继续做,还是实现不了。我将运动方程中加入了浮力,并将方程通过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;
}