Fluent二维单向流泥沙冲刷作用下河床变形代码

#include "udf.h"

#define Rho 1002  /*海水密度*/

#define g 9.81    /*重力加速度*/

#define s 2.65    /*泥沙相对密度*/

#define m 0.4     /*泥沙孔隙率*/

/*#define Nu 0.001148  海水动力粘性系数*/

#define d 0.0005    /*泥沙中值直径*/

/*#define mu 1.146e-6  海水运动粘性系数*/

#define C 2        /*输沙率计算经验系数*/

#define N 1000

/*

Theta—泥沙床面希尔兹数

D—无量纲泥沙粒径

Qb—推移质单宽体积输沙率

Qo—平面单宽体积输沙率

Theta_cr0—水平床面泥沙临界起动希尔兹数

Theta_cr—泥沙临界起动希尔兹数

h—海床面高程

Gamma—海床面倾角

Phi—泥沙休止角

wall_shear_force—壁面剪切力

wall-shear_stree—壁面剪切应力

Theta—网格面泥沙希尔兹数

Theta_v—网格节点泥沙希尔兹数

X1,X2,Y1,Y2—网格节点坐标

Delta_x—网格节点x方向坐标之差

Delta_y—网格节点y方向坐标之差

Delta_Qo—网格节点平面输沙率之差

Delta_h—网格改变量

*/


DEFINE_GRID_MOTION(wall_motion,domain,dt,time,dtime)

{

    Thread *tf=DT_THREAD(dt);

    face_t f;

    Node *v;

    real A[ND_ND];

    real area,wall_shear_force;

    real Theta_cr0,Phi,Gamma;

    real Theta[N]={0.0};

    real Theta_cr[N]={0.0};

    real X[N]={0.0};

    real Y[N]={0.0};

    real Qo[N]={0.0};

    real Delta_h[N]={0.0};

    real Delta_h_v[N]={0.0};

    real wall_shear_stress[N]={0.0};

    real X1,X2,Y1,Y2;

    real a,b,c,e,z1,z2,z3,z4,w1,w2,w3,w4,v1,v2,v3,v4,p1,p2,p3;

    real NV_VEC(Delta_y);

    int i,j,n;

    i=0;

    j=0;

    X2=0;

    Y2=0;

    NV_D(Delta_y,=,0.0,0.0,0.0);

    Theta_cr0=0.047;

    Phi=(32.5+1.27*d*1000)/180;

SET_DEFORMING_THREAD_FLAG(THREAD_T0(tf));

    begin_f_loop(f,tf)

    {

        i++;

        j++;

        F_AREA(A,f,tf);

        area=NV_MAG(A);

        wall_shear_force=F_STORAGE_R_N3V(f,tf,SV_WALL_SHEAR)[0];

        wall_shear_stress[i]=wall_shear_force/area;

        Theta[i]=fabs(wall_shear_stress[i])/(g*(s-1)*d*Rho);

    }

    end_f_loop(f,tf)

    i=0;

    begin_f_loop(f,tf)

    {

        i++;

        f_node_loop(f,tf,n)

        {

            X1=X2;

            Y1=Y2;

            v=F_NODE(f,tf,n);

            X2=NODE_X(v);

            Y2=NODE_Y(v);


            if(n!=0)

            {



                Gamma=atan(fabs(Y2-Y1)/fabs(X2-X1));

                a=1-0.85*0.85*tan(Phi)*tan(Phi);

                b=-sin(Gamma)+tan(Phi)*tan(Phi)*0.85*cos(Gamma);

                c=sin(Gamma)*sin(Gamma)-tan(Phi)*tan(Phi)*cos(Gamma)*cos(Gamma);

                e=tan(Phi)*(1-0.85*tan(Phi));

                Theta_cr[i]=Theta_cr0*((sqrt(b*b-a*c)-b)/e);

            }

        }

    }

    end_f_loop(f,tf)

    i=0;

    begin_f_loop(f,tf)

    {

        i++;

        if(Theta[i]>Theta_cr[i])

        Qo[i]=12*sqrt(g*(s-1)*pow(d,3.0)*Theta[i])*(Theta[i]-Theta_cr[i]);

        else

        Qo[i]=0;

    }

    end_f_loop(f,tf)

    i=1;

    begin_f_loop(f,tf)

    {

        i++;

        f_node_loop(f,tf,n)

        {


        v=F_NODE(f,tf,n);

            if(n!=0)

            {

                X[i]=NODE_X(v);

                Y[i]=NODE_Y(v);

            }


        }

    }

    end_f_loop(f,tf)

    i=1;

    begin_f_loop(f,tf)

    {

        i++;

        if(i<j-1)

        {

            p1=(Qo[i+1]-Qo[i-1])/(2*(X[i+1]-X[i]));

            p2=(Y[i+1]-Y[i-1])/(2*(X[i+1]-X[i]));

            p3=(Y[i+1]-2*Y[i]+Y[i-1])/((X[i+1]-X[i])*(X[i+1]-X[i]));



            if(wall_shear_stress[i]>0)

                Delta_h[i]=dtime*(-p1+C*p1*p2+C*Qo[i]*p3)/(m-1);

            else

                Delta_h[i]=dtime*(p1+C*p1*p2+C*Qo[i]*p3)/(m-1);

        }

    }

    end_f_loop(f,tf)

    i=4;

    begin_f_loop(f,tf)

    {

    i++;

        f_node_loop(f,tf,n)

        {

            if(i<j-4)

            {

                z1=(X[i-2]+X[i-1])/2;

                z2=(X[i-1]+X[i])/2;

                z3=(X[i]+X[i+1])/2;

                z4=(X[i+1]+X[i+2])/2;

                w1=(X[i]-z2)*(X[i]-z3)*(X[i]-z4);

                w2=(X[i]-z1)*(X[i]-z3)*(X[i]-z4);

                w3=(X[i]-z1)*(X[i]-z2)*(X[i]-z4);

                w4=(X[i]-z1)*(X[i]-z2)*(X[i]-z3);

                v1=(z1-z2)*(z1-z3)*(z1-z4);

                v2=(z2-z1)*(z2-z3)*(z2-z4);

                v3=(z3-z1)*(z3-z2)*(z3-z4);

                v4=(z4-z1)*(z4-z2)*(z4-z3);

                Delta_h_v[i]=w1*Delta_h[i-2]/v1+w2*Delta_h[i-1]/v2+w3*Delta_h[i]/v3+w4*Delta_h[i+1]/v4;

            }

        }


    }

    end_f_loop(f,tf)


    i=4;

    begin_f_loop(f,tf)

    {

        i++;

        f_node_loop(f,tf,n)

        {

            v=F_NODE(f,tf,n);

            if(i>4&&NODE_POS_NEED_UPDATE(v)&&i<j-5)

            {

                NODE_POS_UPDATED(v);

                Delta_y[1]=Delta_h_v[i];

                NV_V(NODE_COORD(v),+=,Delta_y);

            }

        }

    }

    end_f_loop(f,tf)

}