校核压力.c

这个源码是从其它网站上转载来的

#include "math.h"
#include "stdio.h"
double fr(k,d,Re,Rr)
double k,d,Re,Rr;
{
	double result;
	result=1/sqrt(Rr)+2.01*log10(k/(3.7065*d)+2.52/(Re*sqrt(Rr)));
	return(result);
}
void main()
{
	double u,x,y,cu,z,Pz,P[3];
	double Pi,T,R0;
	double R,R1,R2,R3,d;
	double de1,c,c0,k;
	double l[4]={35.6*1e3,62.5*1e3,50.8*1e3,59.1*1e3};
	double q[4]={2.5*1e8,2.8*1e8,3.3*1e8,6.39*1e8};
	double Re[4];
	int i;
	c0=0.03848;  de1=0.5624;  T=300.15;  k=0.05*1e-3;
	Pz=2.0*1e6;
	Pi=3.1415926;  R0=0.7272;
	z=9.192528e-001;       d=0.365;
	x=2.57+0.2781*de1+1063.6/T;
	y=1.11+0.04*x;
	cu=2.415*(7.77+0.1844*de1)*pow(T,1.5)/(122.4+377.58*de1+1.8*T)*1e-4;
	u=cu*exp(x*pow((R0/1000),y));
	for(i=0;i<4;i++)
	{
	q[i]=q[i]/(350.0*24*3600);
	Re[i]=1.534*q[i]*de1*1000/(d*u);
	R1=0.001;
	R2=1.0;
	while(1)
	{
		R3=(R1+R2)/2;
		if(fr(k,d,Re[i],R3)>1e-6)   R1=R3;
        if(fr(k,d,Re[i],R3)<-1e-6)  R2=R3;
        if(fabs(fr(k,d,Re[i],R3))<1e-6)  break;
	}
	R=R3;
	c=R*z*de1*T/(c0*c0*pow(d,5));
	P[i]=sqrt(Pz*Pz+c*l[i]*q[i]*q[i]);
	printf("P[%d]=%ePa\n\n\n",i,P[i]);
	Pz=P[i];
	}
}