main.cpp

用于快速的估计任意形状的物体的运动情况受力分子

#include "MESH.h"
#include "math.h"


int main()
{
    mesh data;	            	                                // 网格数据
    double *F = (double *)malloc(3*sizeof(double));		        // 力矢量
    double *M = (double *)malloc(3*sizeof(double));             // 力矩矢量
	double *dF = (double *)malloc(3*sizeof(double));            // 面元力
	double *dMcg = (double *)malloc(3*sizeof(double));          // 面元力矩
    int trans_calc_force(double,double,double,double,double,double,double,double,double,mesh,double *,double *); //气动计算程序

    // 建立椭球体表面网格
    data.NBLOCK = 1;

    data.IMAX = (int *)malloc(data.NBLOCK*sizeof(int));
    data.JMAX = (int *)malloc(data.NBLOCK*sizeof(int));
    data.KMAX = (int *)malloc(data.NBLOCK*sizeof(int));

    data.IMAX[0] = 81;
    data.JMAX[0] = 81;
    data.KMAX[0] = 1;

    data.gridnum = data.IMAX[0] * data.JMAX[0] * data.KMAX[0];    //节点数
    data.facenum = (data.IMAX[0]-1) * (data.JMAX[0]-1);          //面元数

    //------------------------------------------------------------------------
    //节点坐标

    data.gridx = (double *)malloc(data.gridnum*sizeof(double));
    data.gridy = (double *)malloc(data.gridnum*sizeof(double));
    data.gridz = (double *)malloc(data.gridnum*sizeof(double));

    double a = 0.5;
    double b = 0.5;

    int num = 0;
	int prev_num;
	int i,j,k;
	double theta,phi;

    for (k=0;k<data.KMAX[0];k++)
    {
        for (j=0;j<data.JMAX[0];j++)
        {
            phi = 2.0 * PI / float(data.JMAX[0]-1) * j;

            for (i=0;i<data.IMAX[0];i++)
            {
                theta = PI / float(data.IMAX[0]-1) * i;

                data.gridx[num] = a * cos(theta);
                data.gridy[num] = b * sin(theta) * cos(phi);
                data.gridz[num] = b * sin(theta) * sin(phi);

                num += 1;
            }
        }
    }


    //---------------------------------------------------------
    //      面元顶点编号

    data.node1 = (int *)malloc(data.facenum*sizeof(double));
    data.node2 = (int *)malloc(data.facenum*sizeof(double));
    data.node3 = (int *)malloc(data.facenum*sizeof(double));
    data.node4 = (int *)malloc(data.facenum*sizeof(double));

    num = 0;
    prev_num = 0;

	int iblock;

    for(iblock=0;iblock<data.NBLOCK;iblock++)
    {
        for(j=0;j<data.JMAX[iblock]-1;j++)
        {
            for(i=0;i<data.IMAX[iblock]-1;i++)
            {
                data.node1[num] = prev_num + i     + j*(data.IMAX[iblock]);
                data.node4[num] = prev_num + (i+1) + j*(data.IMAX[iblock]);
                data.node2[num] = prev_num + i     + (j+1)*(data.IMAX[iblock]);
                data.node3[num] = prev_num + (i+1) + (j+1)*(data.IMAX[iblock]);
                num = num + 1;
            }
        }

        prev_num += data.IMAX[iblock] * data.JMAX[iblock];
    }

    //---------------------------------------------------------
    //      最大值与最小值，中心坐标，半径

    data.xmax = 0.5;     data.xmin = -0.5;
    data.ymax = 0.5;     data.ymin = -0.5;
    data.zmax = 0.5;     data.zmin = -0.5;

    // 求解条件

    double h = 100;		// 高度
    double L = 1;		// 特征长度
    double Tw = 195.08;	// 壁面温度
    double u = 7500;	// 沿x轴速度
    double v = 0;		// 沿y轴速度
    double w = 0;		// 沿z轴速度
    double omega_x = 0;	        // 绕x轴角速度
    double omega_y = 0;	        // 绕y轴角速度
    double omega_z = 0;	        // 绕z轴角速度

    //赋初值
    F[0] = 0;
    F[1] = 0;
    F[2] = 0;
    M[0] = 0;
    M[1] = 0;
    M[2] = 0;

    // 调用计算子程序
    trans_calc_force(h, L, Tw, u, v, w, omega_x, omega_y, omega_z, data, F, M);

    return 0;

}

//===============================================================================================//
//      计算物体三个方向的受力与力矩                                                             //
//-----------------------------------------------------------------------------------------------//
//      输入：                                                                                   //
//              h			    飞行高度（km）		                                             //
//				L				物体特征长度（m）												 //
//				Tw				物体表面温度（K）												 //
//				u				x方向速度（m/s）												 //
//				v				y方向速度（m/s）												 //
//				w				z方向速度（m/s）												 //
//				b				绕x轴角速度（rad/s）											 //
//				q				绕y轴角速度（rad/s）											 //
//				r				绕z轴角速度（rad/s）											 //
//				data			网格数据												         //
//																								 //		
//      输出：                                                                                   //
//              F               力矢量												             //
//				M				力矩矢量                                                         //
//              Cd              阻力系数														 //
//              Cl              升力系数                                                         //
//              Cz              侧向力系数                                                       //
//              CMx             X方向力矩系数                                                    //
//              CMy             y方向力矩系数                                                    //
//              CMz             z方向力矩系数                                                    //
//===============================================================================================//
int trans_calc_force(double h, double L, double Tw, double u, double v, double w, double b, double q, double r, mesh data, double *F, double *M)
{
		// 各面元法向量
		double *normal_x = (double *)malloc(data.facenum*sizeof(double));
		double *normal_y = (double *)malloc(data.facenum*sizeof(double));
		double *normal_z = (double *)malloc(data.facenum*sizeof(double));
        
		//各面元切向
	    double *parallel_x = (double *)malloc(data.facenum*sizeof(double));
		double *parallel_y = (double *)malloc(data.facenum*sizeof(double));
		double *parallel_z = (double *)malloc(data.facenum*sizeof(double));
        
		//各面元面积
		double *A_face = (double *)malloc(data.facenum*sizeof(double));
        
		//来流方向和各面元切向夹角的正、余弦值
        double *sin_face = (double *)malloc(data.facenum*sizeof(double));
        double *cos_face = (double *)malloc(data.facenum*sizeof(double));

		//各个面元所受力与力矩
		double *dF_x = (double *)malloc(data.facenum*sizeof(double));
		double *dF_y = (double *)malloc(data.facenum*sizeof(double));
		double *dF_z = (double *)malloc(data.facenum*sizeof(double));
		double *dM_x = (double *)malloc(data.facenum*sizeof(double));
		double *dM_y = (double *)malloc(data.facenum*sizeof(double));
		double *dM_z = (double *)malloc(data.facenum*sizeof(double));
        
		//牛顿法计算各面元压力系数
         double *cp_face = (double *)malloc(data.facenum*sizeof(double));
        
		//自由分子流各面元压力
        double *p_face = (double *)malloc(data.facenum*sizeof(double));
        
		//自由分子流各面元剪切应力
        double *tao_face = (double *)malloc(data.facenum*sizeof(double));

		//桥函数中各面元的C_cont
		double *Cp_facecont = (double *)malloc(data.facenum*sizeof(double));
        double *Ctao_facecont = (double *)malloc(data.facenum*sizeof(double));
		
        //桥函数中各面元的C_fm
		double *Cp_facefm = (double *)malloc(data.facenum*sizeof(double));
		double *Ctao_facefm = (double *)malloc(data.facenum*sizeof(double));
		
   	   //最终气动系数
		 //面元
		double *Cp_face = (double *)malloc(data.facenum*sizeof(double));
        double *Ctao_face = (double *)malloc(data.facenum*sizeof(double));
		//总体
		double Cl_cont, Cd_cont, Cz_cont, CMx_cont, CMy_cont, CMz_cont;
		double Cl_fm, Cd_fm, Cz_fm, CMx_fm, CMy_fm, CMz_fm;
		double Cl, Cd, Cz, CMx, CMy, CMz;
	

		// 当地大气状态参数：密度、压力、温度、分子平均自由程、相对分子质量、相对分子常数
        double *rho, *p, *T, *lamda, *Me, R;

		//高度
		double *h_cont, *h_fm;
		
		// 流场特征参数
        double *Ma_face = (double *)malloc(data.facenum*sizeof(double)); 
		double Ma,Kn;
		
		//当地速度比
        double *S_face = (double *)malloc(data.facenum*sizeof(double));

        double gama = 1.4;
        
		int env_calc_atmos(double, double *, double *, double *, double *, double *);  //环境函数
        double erf(double);                                                            //误差函数
		int env_calc_h(double, double *);                                              //反算高度函数

		rho = (double *)malloc(sizeof(double));                  //密度
		p = (double *)malloc(sizeof(double));                    //压力
		T = (double *)malloc(sizeof(double));                    //温度
		lamda = (double *)malloc(sizeof(double));                //分子平均自由程
		Me = (double *)malloc(sizeof(double));                   //相对分子质量

		h_cont = (double *)malloc(sizeof(double));               //下界高度
	    h_fm = (double *)malloc(sizeof(double));                 //上界高度

		double D = 0.5; //半径
	    double ref_A = PI * D * D;  //参考面积

        // 计算大气状态参数
		env_calc_atmos(h, rho, p, T, lamda, Me);

		// 努森数
		Kn = *lamda / L;

        // 来流速度大小
		double vmo;
        vmo = sqrt(u*u + v*v +w*w);

		// 来流马赫数
        R = 8.31 * 1000 / *Me;
        Ma = vmo/sqrt(gama * R * *T);

		// 自由分子流理论相关参数
		//-----------------------------------------------------------------------------------
        const double k = 1.38e-23;       //波尔兹曼常数
        double S , xigema, xigema1;    
        S = vmo/sqrt(2 * R * *T);        //来流速度比
        xigema = 1;                      //法向动量适应系数
        xigema1 = 1;                     //切向动量适应系
   

		int index,num;
        double x1, y1, z1;
        double x2, y2, z2;
	    double x3, y3, z3;
		double x4, y4, z4;
	    double x0, y0, z0;
	    double r1x, r1y, r1z;
	    double r2x, r2y, r2z;
	    double r0x, r0y, r0z;
	    double nx, ny, nz, nr;
	    double tx, ty, tz, tr;
        double dot;
		double test1, test2;
	    double u1,v1,w1;
	    double V1mo;
		double a1,a2,b1,b2,c,S1,S2;


		double xc, yc, zc;

        //坐标轴原点
        xc = 0.5 * ( data.xmax + data.xmin );
        yc = 0.5 * ( data.ymax + data.ymin );
        zc = 0.5 * ( data.zmax + data.zmin );

    for(num=0;num<data.facenum;num++)
    {
        // 面元四个角点的坐标
        index = data.node1[num];
        x1 = data.gridx[index];
        y1 = data.gridy[index];
        z1 = data.gridz[index];

        index = data.node2[num];
        x2 = data.gridx[index];
        y2 = data.gridy[index];
        z2 = data.gridz[index];

        index = data.node3[num];
        x3 = data.gridx[index];
        y3 = data.gridy[index];
        z3 = data.gridz[index];
        
        index = data.node4[num];
        x4 = data.gridx[index];
        y4 = data.gridy[index];
        z4 = data.gridz[index];

        // 面元对角线向量
        r1x = x3 - x1;
        r1y = y3 - y1;
        r1z = z3 - z1;

        r2x = x4 - x2;
        r2y = y4 - y2;
        r2z = z4 - z2;

        //  叉乘: (r3-r1) x (r4-r2)求面元单位外法向
        nx = r1y * r2z - r2y * r1z;
        ny = r1z * r2x - r2z * r1x;
        nz = r1x * r2y - r2x * r1y;
        nr = sqrt(nx*nx + ny*ny + nz*nz);

        if (fabs(nr)<1e-6)
        {
            printf("计算有误");
            return -1;
        }
        
		//面元中心点
        x0 = 0.25 * ( x1 + x2 + x3 + x4 );
        y0 = 0.25 * ( y1 + y2 + y3 + y4 );
        z0 = 0.25 * ( z1 + z2 + z3 + z4 );
        
		//面元中心点矢径
        r0x = x0 - xc;
        r0y = y0 - yc;
        r0z = z0 - zc;

        dot = nx * r0x + ny * r0y + nz * r0z;

        if (dot>0)