第二题.2.c

混凝土结构杆件的非线性计算

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define e 0.0026
#define t 3500 


	/*double gauss(double (*f)(),double a,double b,int n)//高斯公式积分
	{
	double c,d,s;
	double	y[5],h[5];
	int i;
	c=(b-a)/2;
	d=(b+a)/2;
	switch(n)
	{
	case(1):
		y[0]=0;  	h[0]=2;
		break;
	case(2): 
		y[0]=-1/sqrt(3);    h[0]=1;
		y[1]=1/sqrt(3); 	h[1]=1;
		break;
	case(3):
		y[0]=-sqrt(15)/5;   h[0]=5.0/9;
		y[1]=0;             h[1]=8.0/9;
		y[2]=sqrt(15)/5;	h[2]=5.0/9;
		break;
	case(4):
		y[0]=0.8611363;   h[0]=0.3478548;
		y[1]=0.3398810;	  h[1]=0.6521452;
		y[2]=-0.8611363;  h[2]=0.3478548;
		y[3]=-0.3398810;  h[3]=0.6521452;
		break;
	default:
		y[0]=0.9061793;	  h[0]=0.2369269;
		y[1]=0.5384693;   h[1]=0.4786287;
		y[2]=0;           h[2]=0.5688889;
		y[3]=-0.9061793;  h[3]=0.2369269;
		y[4]=-0.5384693;  h[4]=0.4786287;
    	break;
	}
	for(i=0,s=0;  i<n;  i++)
		s+=h[i]*f(c*y[i]+d,b);
//	s*=c;
	return(s*c);
}*/
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,N,M,M1,M2,z,s1,s2,s,p;   //其中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*1000/b;
				   z=gauss(f1,a,b,n);
                   M=(z+s2*(565-35)-t*1000*265)/1000000; //计算弯矩Ne
				   M2=0.1*N*p*7.2*7.2;
				   M1=M-M2;
				   printf("e=%lf\n\nb=%lf\nN=%lf\nM=%lf\nM1=%lf\nM2=%lf\ns1=%lf\ns2=%lf\np=%.8lf\n",e,b,N,M,M1,M2,s1/942.48,s2/942.48,p);
				   fo=fopen("data2.6.txt","w");
				   fprintf(fo,"e=%lf\n\nb=%lf\nN=%lf\nM=%lf\nM1=%lf\nM2=%lf\ns1=%lf\ns2=%lf\np=%.8lf\n",e,b,N,M,M1,M2,s1/942.48,s2/942.48,p);
				   fclose(fo);
			   }
		   }
  }