本代码转载自：http://www.codeforge.cn/read/260028/keti_2d_b.c__html
#include
"udf.h"
#define Rho 1002
#define g 9.81
#define s 2.65
#define m 0.4
#define d 0.0005
#define C 2
#define N
1000
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)
}
