2.124程序.3问.c

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

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define e 0.002
#define e0 355.0  

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);
}

double f(double y, double b)
{
 	return((0.131*e*y/b-20*e*e*y*y/(b*b))*400/(0.00000484+0.002155*e*y/b));  //计算 积分σb dy
 }
double f1(double y, double b)
{
 	return((0.131*e*y/b-20*e*e*y*y/(b*b))*400*(560-b+y)/(0.00000484+0.002155*e*y/b)); //计算积分σb（h0-x+y) dy
 }

void main(void)
{
	FILE *fp,*fo;
	double a,b,b0,y,N,M,z,s1,s2,s,ec;   //其中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+1];
			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]=(z+s2*(565-35))/(y+s2-s1);        //计算轴力ec
			j++;
		}
		j=j-2;
		if(NO[0]<=e0 && NO[1]>e0)
		{
			b0=xh[j];
			break;
		}
	}
	fp=fopen("datas.text","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("datas.text","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;
 	    ec=(z+s2*(565-35))/(y+s2-s1);        //计算轴力ec
     	if(fabs(ec-e0)<=0.02) 
		{
			    N=(y+s2-s1)/1000;      //计算轴力N
				z=gauss(f1,a,b,n);
                M=N*e0/1000000; //计算弯矩Ne
				printf("e=%.10lf\n\nb=%.10lf\nN=%.10lf\nM=%.10lf\ns1=%lf\ns2=%lf\nec=%lf\n",e,b,N,M,s1/942.48,s2/942.48,ec);
		    break;
		}
		
	}
    	
		
}