基于四阶Runge-Kutta的2DOF翼型仿真阻尼异常问题排查

我正在开展两自由度（2DOF）翼型仿真项目，其运动由对应运动方程描述。为数值求解这组包含4个一阶初值问题（IVP）的常微分方程（ODE）系统，我选用了四阶Runge-Kutta方法，但仿真结果出现了非预期的阻尼现象，怀疑Runge-Kutta代码实现存在错误却无法定位。以下是我的Fluent UDF代码，同时附上预期结果图与实际结果图，恳请帮忙排查代码中的问题：

#include "udf.h"
#include "stdio.h"
#include "math.h"

#define PI 3.141592654

static real x1_0 = 0; //  y, 垂向位移
static real x2_0 = 0; //  y_dot, 垂向速度
static real x3_0 = 0;//-15*PI/180; //  theta, 俯仰位移
static real x4_0 = 30*PI/180; //  theta_dot, 俯仰角速度
static real x5_0 = 0; //  V, 电压

//static real current_time = -1;

DEFINE_CG_MOTION(flappingLinear1,dt,vel,omega,time,dtime)
{


//  real ctime = RP_Get_Real("flow-time");
//  real ctimestep = RP_Get_Integer("time-step");
//
//  if (current_time < ctimestep)
//    {
//      current_time = ctimestep;
      Thread *t;
      Domain *d;
      FILE *fp;

      real cg[3];
      real force[3]; 
      real moment[3];

      // 变量定义

      real S_star = -0.04;
      real Ured = 5.0;
      real eta_h = 0.1;
      real eta_theta = 0.1;
      real alpha = 0; // 垂向非线性阻尼系数，线性模型时为0
      real beta = 0;  // 俯仰非线性阻尼系数，线性模型时为0
      real freq_h = 1/Ured;
      real freq_dash = 1.0;
      real m_star = 2.0;
      real r_theta_sq = 0.2;
      
      real theta_L = 0.00155; // 1.55E-3 V/N
      real m_h = 12.387;
      real Cap = 0.00000012; // 120 nF
      real Resis = 100000; // 100 K-Ohm

      real sigma_1 = pow(theta_L,2)/(m_h*Cap*pow(freq_h,2));
      real sigma_2 = 1/(Resis*Cap*freq_h);

      real k_theta3 = 0; // 俯仰非线性刚度系数，线性模型时为0
      real k_h3 = 0; // 垂向非线性刚度系数，线性模型时为0

      real A11 = S_star;
      real A12 = pow(2*PI/Ured,2);
      real A13 = 4*PI*eta_h/Ured;
      real A14 = 4*PI*alpha*pow(freq_h,2)*Ured;
      real A15 = 2/(PI*m_star);

      real B11 = S_star/r_theta_sq;
      real B12 = pow(2*freq_dash*PI/Ured,2);
      real B13 = 4*eta_theta*freq_dash*PI/Ured;
      real B14 = 4*PI*beta*freq_dash*pow(freq_dash,2)*Ured;
      real B15 = 2/(PI*r_theta_sq*m_star);

      real D = 1 - A11*B11*pow(cos(x3_0),2);
      real D1 = 1 - A11*B11;

      real cRv = 1/pow(Ured,2);


      t = DT_THREAD(dt);
      d = Get_Domain(1);
      
      Compute_Force_And_Moment(d,t,cg,force,moment,TRUE);

      real Coeff_L = 2 * force[1];
      real Coeff_D = 2 * force[0];
      real Coeff_M = 2 * moment[2];

      #define a(x1,x2,x3,x4,x5) x2 //x1_dot
      #define b(x1,x2,x3,x4,x5) (1/D) * (A15*Coeff_L - A11*pow(x4_0,2)*sin(x3_0) \
                                + (A11*cos(x3_0) * (B15*Coeff_M - B12*(x3_0 + k_theta3*pow(x3_0,3)) \
                                - B13*x4_0 - B14*pow(x4_0,3))) - A12*(x1_0 + k_h3*pow(x1_0,3)) - A13*x2_0 \
                                - A14*pow(x2_0,3) + cRv*x5_0)
      #define c(x1,x2,x3,x4,x5) x4 //x3_dot
      #define d(x1,x2,x3,x4,x5) (1/D) * (B15*Coeff_M + (B11*cos(x3_0) * (A15*Coeff_L - A12*(x1_0 + k_h3*pow(x1_0,3)) \
                                - A13*x2_0 - A14*pow(x2_0,3))) - B12*(x3_0 + k_theta3*pow(x3_0,3)) - B13*x4_0 \
                                - B14*pow(x4_0,3) + cRv*x5_0)
      #define e(x1,x2,x3,x4,x5) -(sigma_1*x2_0) - (sigma_2*x5_0/Ured) //x5_dot

      // 声明所有用到的变量
      double k1, l1, m1, n1, o1,
             k2, l2, m2, n2, o2,
             k3, l3, m3, n3, o3,
             k4, l4, m4, n4, o4,
             k, l, m, n, o,
             x1, x2, x3, x4, x5;

        k1 = dtime * (a(x1_0, x2_0, x3_0, x4_0, x5_0));
        l1 = dtime * (b(x1_0, x2_0, x3_0, x4_0, x5_0));
        m1 = dtime * (c(x1_0, x2_0, x3_0, x4_0, x5_0));
        n1 = dtime * (d(x1_0, x2_0, x3_0, x4_0, x5_0));
        o1 = dtime * (e(x1_0, x2_0, x3_0, x4_0, x5_0));


        k2 = dtime * (a((x1_0 + k1 / 2), (x2_0 + l1 / 2),(x3_0 + m1 / 2),(x4_0 + n1 / 2),(x5_0 + o1 / 2)));
        l2 = dtime * (b((x1_0 + k1 / 2), (x2_0 + l1 / 2),(x3_0 + m1 / 2),(x4_0 + n1 / 2),(x5_0 + o1 / 2)));
        m2 = dtime * (c((x1_0 + k1 / 2), (x2_0 + l1 / 2),(x3_0 + m1 / 2),(x4_0 + n1 / 2),(x5_0 + o1 / 2)));
        n2 = dtime * (d((x1_0 + k1 / 2), (x2_0 + l1 / 2),(x3_0 + m1 / 2),(x4_0 + n1 / 2),(x5_0 + o1 / 2)));
        o2 = dtime * (e((x1_0 + k1 / 2), (x2_0 + l1 / 2),(x3_0 + m1 / 2),(x4_0 + n1 / 2),(x5_0 + o1 / 2)));

        k3 = dtime * (a((x1_0 + k2 / 2), (x2_0 + l2 / 2),(x3_0 + m2 / 2),(x4_0 + n2 / 2),(x5_0 + o2 / 2)));
        l3 = dtime * (b((x1_0 + k2 / 2), (x2_0 + l2 / 2),(x3_0 + m2 / 2),(x4_0 + n2 / 2),(x5_0 + o2 / 2)));
        m3 = dtime * (c((x1_0 + k2 / 2), (x2_0 + l2 / 2),(x3_0 + m2 / 2),(x4_0 + n2 / 2),(x5_0 + o2 / 2)));
        n3 = dtime * (d((x1_0 + k2 / 2), (x2_0 + l2 / 2),(x3_0 + m2 / 2),(x4_0 + n2 / 2),(x5_0 + o2 / 2)));
        o3 = dtime * (e((x1_0 + k2 / 2), (x2_0 + l2 / 2),(x3_0 + m2 / 2),(x4_0 + n2 / 2),(x5_0 + o2 / 2)));

        k4 = dtime * (a((x1_0 + k3), (x2_0 + l3),(x3_0 + m3),(x4_0 + n3),(x5_0 + o3)));
        l4 = dtime * (b((x1_0 + k3), (x2_0 + l3),(x3_0 + m3),(x4_0 + n3),(x5_0 + o3)));
        m4 = dtime * (c((x1_0 + k3), (x2_0 + l3),(x3_0 + m3),(x4_0 + n3),(x5_0 + o3)));
        n4 = dtime * (d((x1_0 + k3), (x2_0 + l3),(x3_0 + m3),(x4_0 + n3),(x5_0 + o3)));
        o4 = dtime * (d((x1_0 + k3), (x2_0 + l3),(x3_0 + m3),(x4_0 + n3),(x5_0 + o3)));

        k = (k1 + 2 * k2 + 2 * k3 + k4) / 6;
        l = (l1 + 2 * l2 + 2 * l3 + l4) / 6;
        m = (m1 + 2 * m2 + 2 * m3 + m4) / 6;
        n = (n1 + 2 * n2 + 2 * n3 + n4) / 6;
        o = (o1 + 2 * o2 + 2 * o3 + o4) / 6;

        x1 = x1_0 + k;
        x2 = x2_0 + l;
        x3 = x3_0 + m;
        x4 = x4_0 + n;
        x5 = x5_0 + o;

        // 更新状态变量
        x1_0 = x1;
        x2_0 = x2;
        x3_0 = x3;
        x4_0 = x4;
        x5_0 = x5;

        vel[1] = x2;
        omega[2] = x4;

        // 弧度转角度
        real pitch_angle = x3*57.3;


        #if RP_HOST
        fp = fopen("data.txt", "a");
        fprintf(fp, "%.5f %.10f %.10f %.10f %.10f %.10f %.10f %.10f %.10f %.10f\n", 
          time, x1, x2, x3, x4, x5, Coeff_L, Coeff_D, Coeff_M, pitch_angle);
        fclose(fp);
        #endif
//      }
//
//  vel[1] = x2_0;
//  omega[2] = x4_0;

}

DEFINE_CG_MOTION(FLBoundary1,dt,vel,omega,time,dtime)
{
  vel[1] = x2_0;
  omega[2] = x4_0;
}

预期结果图 #

实际结果图 #

内容的提问来源于stack exchange，提问作者Amir Khan