该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
#include
#include
#define PI 3.1415926
void main()
{
int i,j,k,l,m,n,p,q,u;//p,q为i点和j点的读取指针
int t,step=1;//定义时间步迭代变量,及总的迭代时间步数,
double L,B,H;//定义长宽高的迭代数
double Length=0.1,Broad=0.05,High=0.005,density=2400,dt=0.0000001;//定义材料的长宽高及密度,dt为时间间隔
double inx=0.05,iny=0.025,inz=0.005;//定义入射点坐标
double Grid=0.0005,Hor=0.0015,R=0.005,r;//定义晶格长度和近场范围半径,以及冲击物的半径,r为距离冲击物中心的距离
double Fr,sumf[3];//定义规定的面力密度,及各点的初始冲击力
double x,y,z,xj,yj,zj;//定义i点及j点的实际坐标形式,。
static double a[300000][3],v[300000][3],s[300000][3];//a为加速度,v为速度,s为位移
double K=20*pow(10,9),c,s0=0.005,vr,ks=pow(10,17);//K为体积模量,c为弹性常量,ks为力常数
double rp[3],rd[3],ap,s1,ad;//定义相对位置及相对位移,以及变形后的距离,s1为伸长率,ad为相对位置
static double sc[300000][3];
c=18*K/(PI*pow(Hor,4));
L=Length/Grid;B=Broad/Grid;H=High/Grid;
for(t=0;t<=step;t++)
{
p=0;
for(i=0;i<=(int)(L-1);i++)
for(j=0;j<=(int)(B-1);j++)
for(k=0;k<=(int)(H-1);k++)
{
x=i*Grid+Grid/2.0;y=j*Grid+Grid/2.0;z=k*Grid+Grid/2.0;//计算i点的实际坐标
sumf[0]=0;sumf[1]=0;sumf[2]=0;
if(t==0)
{
r=sqrt(pow(x-inx,2)+pow(y-iny,2)+pow(z-inz,2));//计算计算点离入射点的距离
if(r
Fr=-ks*pow(r-R,2);
else
Fr=0;
a[p][0]=0;a[p][1]=0;a[p][2]=Fr/(density*pow(Grid,3));
v[p][0]=0;v[p][1]=0;v[p][2]=a[p][2]*dt;
s[p][0]=0;s[p][1]=0;s[p][2]=v[p][2]*dt/2.0;
p++;
}
else
{
q=0;
for(l=0;l<=L-1;l++)
for(m=0;m<=B-1;m++)
for(n=0;n<=H-1;n++)
{
xj=l*Grid+Grid/2.0;yj=m*Grid+Grid/2.0;zj=n*Grid+Grid/2.0;//计算j点的实际左边
rp[0]=xj-x;rp[1]=yj-y;rp[2]=zj-z;//计算相对位置
rd[0]=s[q][0]-s[p][0];rd[1]=s[q][1]-s[p][1];rd[2]=s[q][2]-s[p][2];//计算相对位移
ap=sqrt(pow(rp[0]+rd[0],2)+pow(rp[1]+rd[1],2)+pow(rp[2]+rd[2],2));
ad=sqrt(rp[0]*rp[0]+rp[1]*rp[1]+rp[2]*rp[2]);
s1=(ap-ad)/ad;
if(s1
u=1;
else
u=0;
if(ad<=Hor-0.5*Grid)
vr=1;
else
if(ad>(Hor-0.5*Grid)&&ad<=(Hor+0.5*Grid))
vr=(Hor+0.5*Grid-ad)/Grid;
else
vr=0;
sumf[0]=sumf[0]+c*s1*u*(rd[0]+rp[0])*vr*pow(Grid,3)/ap;
sumf[1]=sumf[1]+c*s1*u*(rd[1]+rp[1])*vr*pow(Grid,3)/ap;
sumf[2]=sumf[2]+c*s1*u*(rd[2]+rp[2])*vr*pow(Grid,3)/ap;
q++;
}
a[p][0]=sumf[0]/density;a[p][1]=sumf[1]/density;a[p][2]=sumf[2]/density;
s[p][0]=s[p][0]+v[p][0]*dt+a[p][0]*dt*dt/2.0;s[p][1]=s[p][1]+v[p][1]*dt+a[p][1]*dt*dt/2.0;s[p][2]=s[p][2]+v[p][2]*dt+a[p][2]*dt*dt/2.0;
v[p][0]=v[p][0]+a[p][0]*dt;v[p][1]=v[p][1]+a[p][1]*dt;v[p][2]=v[p][2]+a[p][2]*dt;
p++;
}
}
}
p=0;
for(i=0;i<=L-1;i++)
for(j=0;j<=B-1;j++)
for(k=0;k<=H-1;k++)
{
x=i*Grid+Grid/2.0;y=j*Grid+Grid/2.0;z=k*Grid+Grid/2.0;
sc[p][0]=x+s[p][0];sc[p][1]=y+s[p][1];sc[p][2]=z+s[p][2];
printf("%.16lf\t%.16lf\t%.16lf\n",sc[p][0],sc[p][1],sc[p][2]);
p++;
}
}