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

CFD中文网

liuxiaodong