第二题.4.c

来自「混凝土结构杆件的非线性计算」· C语言 代码 · 共 155 行

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define e 0.0022885
#define t 1150


double gauss(double(*f)( ), double a ,double b,int n)
{
	double *p,*q;
	double s,c,d;
	int i;
	p=(double*)malloc(n*sizeof(double));//存放权系数 
	q=(double*)malloc(n*sizeof(double));//存放高斯积分点 
	switch(n)//选择不同积分点数
	{
	case(1):
		p[0]=0.0;
		q[0]=2.0;
	break;
    case(2):
        p[0]=-0.5773502692;p[1]=0.5773502692;
		q[0]=1.0;q[1]=1.0;
    break;
    case(3):
        p[0]=-0.7745966692;p[1]=0.0;p[2]=0.7745966692;
        q[0]=5/9;q[1]=8/9;q[2]=5/9;
    break;
    case(4):
        p[0]=-0.8611363116;p[1]=-0.3399810436;
        p[2]=0.3399810436;p[3]=0.8611363116;
        q[0]=0.3478548451;q[1]=0.6521451549;
		q[2]=0.6521451549;q[3]=0.3478548451;
    break;
    case(5):
        p[0]=-0.9061798459;p[1]=-0.5384693101;p[2]=0.0;
        p[3]=0.5384693101;p[4]=0.9061798459;
        q[0]=0.2369268851;q[1]=0.4786286705;
		q[2]=0.5688888889;q[3]=0.4786286705;q[4]=0.2369268851;
    break;
    case(6):
        p[0]=-0.9324695142;p[1]=-0.6612093865;p[2]=-0.2366191861;
        p[3]=0.2366191861;p[4]=0.6612093865;p[5]=0.9324695142;
        q[0]=0.1713244924;q[1]=0.3607615730;q[2]=0.4679139346;
		q[3]=0.4679139346;q[4]=0.3607615730;q[5]=0.1713244924;
    break;
    default:
        p[0]=-0.9910791230;p[1]=-0.7415311856;p[2]=-0.4058451514;p[3]=0.0;
        p[4]=0.4058451514;p[5]=0.7415311856;p[6]=0.9910791230;
        q[0]=0.1294849662;q[1]=0.2797053915;q[2]=0.3818300505;q[3]=0.4179591337;
		q[4]=0.3818300505;q[5]=0.2797053915;q[6]=0.1294849662;
	}
	c=(b-a)/2; d=(b+a)/2;
	for(i=0,s=0.0; i<n; i++)//高斯公式
	{
		s+=q[i]*f(c*p[i]+d,b);
	}
	return(s*c);
}
double f(double y, double b)
{
	return((27088.30258*e*y/b-4132331.404*e*e*y*y/(b*b))*400/(1+445.324195*e*y/b));  //计算 积分σb dy
}
double f1(double y, double b)
{
	return((27088.30258*e*y/b-4132331.404*e*e*y*y/(b*b))*400*(565-b+y)/(1+445.324195*e*y/b)); //计算积分σb（h0-x+y) dy
}
  void main(void) 
  {
            
			FILE *fp,*fo;
			double a,b,b0,y,p,ff;
			double N,M,M2,M0,z,s1,s2,s;   //其中s1为受拉钢筋的合力;s2为受压钢筋的合力;s为钢筋屈服时的合力
	        double *x,*num;
	        double b1,*xh,NO[2];
	        int    n,i,j,k;
			
	        a=0.0;
	        n=5;
	        s=942.48*455;
            xh=(double*)malloc(601*sizeof(double));
	        for(i=0;i<=600;i++) 
				xh[i]=i;
			for(j=0; j<600;j++)
			{
				for(k=0; k<2; k++)
				{
					b1=xh[j];
			        y=gauss(f,a,b1,n);
			        z=gauss(f1,a,b1,n);
			        s1=200000*e*(565-b1)*942.48/b1;
			        s2=200000*e*(b1-35)*942.48/b1;
			        if(s2>=s) s2=s;
			        if(s1>=s) s1=s;
			        NO[k]=(y+s2-s1)/1000;          //计算轴力N
			        j++;
				}
				j=j-2;
				if(NO[0]<=t && NO[1]>t)
				{
					b0=xh[j];
					break;
				}
			}
			fp=fopen("data.txt","w");
	        num=(double*)malloc(40001*sizeof(double));
	        num[0]=b0;
	        for(k=0;k<=40000;k++)
			{ 
				num[k+1]=num[k]+0.000025;
				fprintf(fp,"%.10lf\n",num[k]);
			}
			fclose(fp);
	       fo=fopen("data.txt","r");
	       x=(double*)malloc(40001*sizeof(double));
	       for(i=0; i<=40000; i++ )
			   fscanf(fo,"%lf",x+i);
		   fclose(fo);
	       for(i=0; i<=40000;i++)
		   {
			   b=x[i];
		       y=gauss(f,a,b,n);   
		       z=gauss(f1,a,b,n); 
		       s2=200000*e*(b-35)*942.48/b;
		       s1=200000*e*(565-b)*942.48/b;
		       if(s2>=s) s2=s;
		       if(s1>=s) s1=s;
 	           N=(y+s2-s1)/1000;      //计算轴力N
     	       if(fabs(N-t)<=0.0001) 
			   {
				
				   p=e/b;
				   ff=0.1*p*7.2*7.2;
				   z=gauss(f1,a,b,n);
                   M=t*((z+s2*(565-35))/(N*1000)-265)/1000; //计算弯矩Ne0,p为曲率
				   M2=t*ff*1000;
				   M0=t*0.355+M2;
                   printf("e=%lf\n\nb=%lf\nN=%lf\nM=%lf\nM0=%lf\nM2=%lf\ns1=%lf\ns2=%lf\np=%.10lf\nff=%.10lf\n\n",e,b,N,M,M0,M2,s1/942.48,s2/942.48,p*1000000,ff*1000000);
				   
				   fo=fopen("data4.1.txt","w");
				   fprintf(fo,"e=%lf\n\nb=%lf\nN=%lf\nM=%lf\nM0=%lf\nM2=%lf\ns1=%lf\ns2=%lf\np=%.10lf\nff=%.10lf\n\n",e,b,N,M,M0,M2,s1/942.48,s2/942.48,p*1000000,ff*1000000);
				   fclose(fo);
			   }
		   }
  }