您所在位置：主题：[原创]一个N个自称C高手的哥们一周了还没解决的问题!!寻达人指点!!
yamalu
[专家分：0] 发布于 2007-05-24 16:44:00
寻找达人指点迷津!!这个程序的原程序在"i"不分奇偶的情况下可以运行,但是我现在想改变算法,"i"分奇偶计算,结果输出确实乱码?请教过N位高手,皆沉!!!故来次寻找高高手!!谢过
#define M 20
#define N 250000
#define eps 1e-6
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include <time.h>
main()
{
int i,j;
FILE *fp;
double D,h,t,A;
double R,P1,P2,Cv,miu,Pr,c,rou0;
double *theta,*rou,*u,*p,*Re,*lamta,*Nu,*H,*k,*u_ever;
double *theta1,*rou1,*u1,*p1,*Re1,*lamta1,*Nu1,*H1,*k1,*u_ever1;
double *theta2,*rou2,*u2,*p2,*Re2,*lamta2,*Nu2,*H2,*k2,*u_ever2;
double start,finish; /* 开始时间,结束时间 */
start=(double)clock(); /*精确到ms */
// _t tnow1,tnow2;
theta=(double *)malloc((M+2)*sizeof(double));
rou=(double *)malloc((M+2)*sizeof(double));
u=(double *)malloc((M+2)*sizeof(double));
p=(double *)malloc((M+2)*sizeof(double));
Re=(double *)malloc((M+2)*sizeof(double));
lamta=(double *)malloc((M+2)*sizeof(double));
Nu=(double *)malloc((M+2)*sizeof(double));
H=(double *)malloc((M+2)*sizeof(double));
k=(double *)malloc((M+2)*sizeof(double));
u_ever=(double *)malloc((M+2)*sizeof(double));
theta1=(double *)malloc((M+2)*sizeof(double));
rou1=(double *)malloc((M+2)*sizeof(double));
u1=(double *)malloc((M+2)*sizeof(double));
p1=(double *)malloc((M+2)*sizeof(double));
Re1=(double *)malloc((M+2)*sizeof(double));
lamta1=(double *)malloc((M+2)*sizeof(double));
Nu1=(double *)malloc((M+2)*sizeof(double));
H1=(double *)malloc((M+2)*sizeof(double));
k1=(double *)malloc((M+2)*sizeof(double));
u_ever1=(double *)malloc((M+2)*sizeof(double));
/*
依次定义：
变量 管道长度，截面面积，管道直径，仿真总时间，空间步长，时间步长；
气体常数，入口压强，出口压强，定容比热容，动力黏度系数，雷诺数，摩擦系数，普朗特数，努塞尔系数、空气热传导率，传热率
数组 温度，密度，流速*/
D=0.004;
A=3.1415926/4*D*D;
/*
h=L/M; /*空间步长h=0.1* /
t=T/N; /*时间步长t=1e-5,t/h=1e-4*/
h=1;
t=1e-5;
R=287.1; P1=200000; P2=101325; Cv=718; Pr=0.72;
c=0.2;rou0=1.204;
for(i=0;i<=M+1;i++) /*所有格子零时刻的初始值*/
{
theta[i]=293.15;
rou[i]=P1/(R*theta[i]);
p[i]=200000;
u[i]=0;
}
/*计算数值解主体程序*/
for(j=0;j<=N;j++)
{
theta1[0]=theta[1]; /*第0号格子所有时刻对应的值*/
rou1[0]=rou[1];
p1[0]=p[1];
u1[0]=0;
for(i=1;i<=M;i++)
{
if(i%2==0)
{
rou1[i]=(rou1[i-1]+rou1[i+1])/2 ;
u1[i]=(u1[i-1]+u1[i+1])/2 ;
lamta1[i]=(lamta1[i-1]+lamta1[i+1])/2;
p1[i]=(p1[i-1]+p1[i+1])/2 ;
}
if(i%2==1||i==20)
{
miu=(4.8e-8)*(theta[0]-273.15)+1.71e-5; /*动力黏度系数*/
//if(i%2==1){
Re1[i]=0.5*(rou[i]+rou[i-1])*fabs(u[i])*D/miu; /*雷诺数*/
if(Re1[i]<eps) lamta1[i]=0;
else lamta1[i]=0.3164*pow(Re1[i],-0.25); /*摩擦系数*/
/*普朗特数*/
Nu1[i]=0.023*pow(Re1[i],0.8)*pow(Pr,0.4); /*努塞尔系数*/
k1[i]=7.95e-5*theta[i]+2.0465e-3; /*空气热传导率*/
H1[i]=2*Nu1[i]*k1[i]/D; /* 传热率*/
if(u[i]>0)
{
u1[i]=u[i]-t/h*u[i]*(u[i]-u[i-1])-2/(rou[i]+rou[i-1])*t/h*(p[i]-p[i-1])-lamta1[i]*t/(2*D)*pow(u[i],2);
}
else
{
u1[i]=u[i]-t/h*u[i]*(u[i+1]-u[i])-2/(rou[i]+rou[i-1])*t/h*(p[i]-p[i-1])+lamta[i]*t/(2*D)*pow(u[i],2);
}
u_ever[i]=0.5*(u[i]+u[i+1]);
if(u_ever[i]>0)
{
rou1[i]=rou[i]-t/h*rou[i]*(u[i+1]-u[i])-t/h*u_ever[i]*(rou[i]-rou[i-1]);
theta1[i]=theta[i]-t*4*H1[i]*(theta[i]-293.15)/(rou[i]*Cv*D)-t/h*u_ever[i]*(theta[i]-theta[i-1])-t/h*R*theta[i]*(u[i+1]-u[i])/Cv+0.5*lamta1[i]*t/(Cv*D)*pow(u_ever[i],3);
}
else
{
rou1[i]=rou[i]-t/h*rou[i]*(u[i+1]-u[i])-t/h*u_ever[i]*(rou[i+1]-rou[i]);
theta1[i]=theta[i]-t*4*H1[i]*(theta[i]-273.15)/(rou[i]*Cv*D)-t/h*u_ever[i]*(theta[i+1]-theta[i])-t/h*R*theta[i]*(u[i+1]-u[i])/Cv+0.5*lamta1[i]*t/(Cv*D)*fabs(u_ever[i])*pow(u_ever[i],2);
}
p1[i]=rou1[i]*R*theta1[i];
}
}
/*for(i=1;i<=M;i++)
{
if(i%2==0)
{
rou1[i]=(rou1[i-1]+rou1[i+1])/2 ;
u1[i]=(u1[i-1]+u1[i+1])/2 ;
lamta1[i]=(lamta1[i-1]+lamta1[i+1])/2;
p1[i]=(p1[i-1]+p1[i+1])/2 ;
}
}*/
if(P2/p1[M]<=0.5)
u1[M+1]=1e-8*c*rou0*287.1*theta1[M]/A*sqrt(293.15/theta1[M]);
else if
(P2/p1[M]<=1)
u1[M+1]=1e-8*c*rou0*287.1*theta1[M]/A*sqrt(293.15/theta1[M])*sqrt(1-4*pow((P2/p1[M]-0.5),2));
else
u1[M+1]=-P2/p1[M]*1e-8*c*rou0*287.1*theta1[M]/A*sqrt(293.15/theta1[M])*sqrt(1-4*pow((p1[M]/P2-0.5),2));
theta1[M+1]=theta[M]; /*第M+1号格子其它参数的值*/
rou1[M+1]=rou[M];
p1[M+1]=p[M];
if(j%1000==0) printf("%d\n", j);
if(j==0)
{
fp=fopen("D:\\theta1.txt","w+");
fprintf(fp,"%f %f %f %f %f %f %f %f %f \n",p1[1],p1[M],theta1[1],theta1[M],u1[1],u1[M],rou1[1],rou1[M],j*t);
}
else if(j%100==0)
{
fp=fopen("D:\\theta1.txt","a+");
fprintf(fp,"%f %f %f %f %f %f %f %f %f \n",p1[1],p1[M],theta1[1],theta1[M],u1[1],u1[M],rou1[1],rou1[M],j*t);
}
fclose(fp);
theta2=theta;
rou2=rou;
u2=u;
p2=p;
Re2=Re;
lamta2=lamta;
Nu2=Nu;
H2=H;
k2=k;
u_ever2=u_ever;
theta=theta1;
rou=rou1;
u=u1;
p=p1;
Re=Re1;
lamta=lamta1;
Nu=Nu1;
H=H1;
k=k1;
u_ever=u_ever1;
theta1=theta2;
rou1=rou2;
u1=u2;
p1=p2;
Re1=Re2;
lamta1=lamta2;
Nu1=Nu2;
H1=H2;
k1=k2;
u_ever1=u_ever2;
}
finish=(double)clock();
printf("%.4fms",(finish-start));
//(&tnow2);
//printf("程序运行时间为 ： %d\n",tnow2-tnow1);
free(theta);
free(rou);
free(u);
free(p);
free(Re);
free(lamta);
free(Nu);
free(H);
free(u_ever);
free(k);
free(rou1);
free(u1);
free(p1);
free(Re1);
free(lamta1);
free(Nu1);
free(H1);
free(u_ever1);
free(k1);
free(theta1);
}
#define M 20
#define N 250000
#define eps 1e-6
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include <time.h>
main()
{
int i,j;
FILE *fp;
double D,h,t,A;
double R,P1,P2,Cv,miu,Pr,c,rou0;
double *theta,*rou,*u,*p,*Re,*lamta,*Nu,*H,*k,*u_ever;
double *theta1,*rou1,*u1,*p1,*Re1,*lamta1,*Nu1,*H1,*k1,*u_ever1;
double *theta2,*rou2,*u2,*p2,*Re2,*lamta2,*Nu2,*H2,*k2,*u_ever2;
double start,finish; /* 开始时间,结束时间 */
start=(double)clock(); /*精确到ms */
// _t tnow1,tnow2;
theta=(double *)malloc((M+2)*sizeof(double));
rou=(double *)malloc((M+2)*sizeof(double));
u=(double *)malloc((M+2)*sizeof(double));
p=(double *)malloc((M+2)*sizeof(double));
Re=(double *)malloc((M+2)*sizeof(double));
lamta=(double *)malloc((M+2)*sizeof(double));
Nu=(double *)malloc((M+2)*sizeof(double));
H=(double *)malloc((M+2)*sizeof(double));
k=(double *)malloc((M+2)*sizeof(double));
u_ever=(double *)malloc((M+2)*sizeof(double));
theta1=(double *)malloc((M+2)*sizeof(double));
rou1=(double *)malloc((M+2)*sizeof(double));
u1=(double *)malloc((M+2)*sizeof(double));
p1=(double *)malloc((M+2)*sizeof(double));
Re1=(double *)malloc((M+2)*sizeof(double));
lamta1=(double *)malloc((M+2)*sizeof(double));
Nu1=(double *)malloc((M+2)*sizeof(double));
H1=(double *)malloc((M+2)*sizeof(double));
k1=(double *)malloc((M+2)*sizeof(double));
u_ever1=(double *)malloc((M+2)*sizeof(double));
/*
依次定义：
变量 管道长度，截面面积，管道直径，仿真总时间，空间步长，时间步长；
气体常数，入口压强，出口压强，定容比热容，动力黏度系数，雷诺数，摩擦系数，普朗特数，努塞尔系数、空气热传导率，传热率
数组 温度，密度，流速*/
D=0.004;
A=3.1415926/4*D*D;
/*
h=L/M; /*空间步长h=0.1* /
t=T/N; /*时间步长t=1e-5,t/h=1e-4*/
h=1;
t=1e-5;
R=287.1; P1=200000; P2=101325; Cv=718; Pr=0.72;
c=0.2;rou0=1.204;
for(i=0;i<=M+1;i++) /*所有格子零时刻的初始值*/
{
theta[i]=293.15;
rou[i]=P1/(R*theta[i]);
p[i]=200000;
u[i]=0;
}
/*计算数值解主体程序*/
for(j=0;j<=N;j++)
{
theta1[0]=theta[1]; /*第0号格子所有时刻对应的值*/
rou1[0]=rou[1];
p1[0]=p[1];
u1[0]=0;
for(i=1;i<=M;i++)
{
if(i%2==0)
{
rou1[i]=(rou1[i-1]+rou1[i+1])/2 ;
u1[i]=(u1[i-1]+u1[i+1])/2 ;
lamta1[i]=(lamta1[i-1]+lamta1[i+1])/2;
p1[i]=(p1[i-1]+p1[i+1])/2 ;
}
if(i%2==1||i==20)
{
miu=(4.8e-8)*(theta[0]-273.15)+1.71e-5; /*动力黏度系数*/
//if(i%2==1){
Re1[i]=0.5*(rou[i]+rou[i-1])*fabs(u[i])*D/miu; /*雷诺数*/
if(Re1[i]<eps) lamta1[i]=0;
else lamta1[i]=0.3164*pow(Re1[i],-0.25); /*摩擦系数*/
/*普朗特数*/
Nu1[i]=0.023*pow(Re1[i],0.8)*pow(Pr,0.4); /*努塞尔系数*/
k1[i]=7.95e-5*theta[i]+2.0465e-3; /*空气热传导率*/
H1[i]=2*Nu1[i]*k1[i]/D; /* 传热率*/
if(u[i]>0)
{
u1[i]=u[i]-t/h*u[i]*(u[i]-u[i-1])-2/(rou[i]+rou[i-1])*t/h*(p[i]-p[i-1])-lamta1[i]*t/(2*D)*pow(u[i],2);
}
else
{
u1[i]=u[i]-t/h*u[i]*(u[i+1]-u[i])-2/(rou[i]+rou[i-1])*t/h*(p[i]-p[i-1])+lamta[i]*t/(2*D)*pow(u[i],2);
}
u_ever[i]=0.5*(u[i]+u[i+1]);
if(u_ever[i]>0)
{
rou1[i]=rou[i]-t/h*rou[i]*(u[i+1]-u[i])-t/h*u_ever[i]*(rou[i]-rou[i-1]);
theta1[i]=theta[i]-t*4*H1[i]*(theta[i]-293.15)/(rou[i]*Cv*D)-t/h*u_ever[i]*(theta[i]-theta[i-1])-t/h*R*theta[i]*(u[i+1]-u[i])/Cv+0.5*lamta1[i]*t/(Cv*D)*pow(u_ever[i],3);
}
else
{
rou1[i]=rou[i]-t/h*rou[i]*(u[i+1]-u[i])-t/h*u_ever[i]*(rou[i+1]-rou[i]);
theta1[i]=theta[i]-t*4*H1[i]*(theta[i]-273.15)/(rou[i]*Cv*D)-t/h*u_ever[i]*(theta[i+1]-theta[i])-t/h*R*theta[i]*(u[i+1]-u[i])/Cv+0.5*lamta1[i]*t/(Cv*D)*fabs(u_ever[i])*pow(u_ever[i],2);
}
p1[i]=rou1[i]*R*theta1[i];
}
}
/*for(i=1;i<=M;i++)
{
if(i%2==0)
{
rou1[i]=(rou1[i-1]+rou1[i+1])/2 ;
u1[i]=(u1[i-1]+u1[i+1])/2 ;
lamta1[i]=(lamta1[i-1]+lamta1[i+1])/2;
p1[i]=(p1[i-1]+p1[i+1])/2 ;
}
}*/
if(P2/p1[M]<=0.5)
u1[M+1]=1e-8*c*rou0*287.1*theta1[M]/A*sqrt(293.15/theta1[M]);
else if
(P2/p1[M]<=1)
u1[M+1]=1e-8*c*rou0*287.1*theta1[M]/A*sqrt(293.15/theta1[M])*sqrt(1-4*pow((P2/p1[M]-0.5),2));
else
u1[M+1]=-P2/p1[M]*1e-8*c*rou0*287.1*theta1[M]/A*sqrt(293.15/theta1[M])*sqrt(1-4*pow((p1[M]/P2-0.5),2));
theta1[M+1]=theta[M]; /*第M+1号格子其它参数的值*/
rou1[M+1]=rou[M];
p1[M+1]=p[M];
if(j%1000==0) printf("%d\n", j);
if(j==0)
{
fp=fopen("D:\\theta1.txt","w+");
fprintf(fp,"%f %f %f %f %f %f %f %f %f \n",p1[1],p1[M],theta1[1],theta1[M],u1[1],u1[M],rou1[1],rou1[M],j*t);
}
else if(j%100==0)
{
fp=fopen("D:\\theta1.txt","a+");
fprintf(fp,"%f %f %f %f %f %f %f %f %f \n",p1[1],p1[M],theta1[1],theta1[M],u1[1],u1[M],rou1[1],rou1[M],j*t);
}
fclose(fp);
theta2=theta;
rou2=rou;
u2=u;
p2=p;
Re2=Re;
lamta2=lamta;
Nu2=Nu;
H2=H;
k2=k;
u_ever2=u_ever;
theta=theta1;
rou=rou1;
u=u1;
p=p1;
Re=Re1;
lamta=lamta1;
Nu=Nu1;
H=H1;
k=k1;
u_ever=u_ever1;
theta1=theta2;
rou1=rou2;
u1=u2;
p1=p2;
Re1=Re2;
lamta1=lamta2;
Nu1=Nu2;
H1=H2;
k1=k2;
u_ever1=u_ever2;
}
finish=(double)clock();
printf("%.4fms",(finish-start));
//(&tnow2);
//printf("程序运行时间为 ： %d\n",tnow2-tnow1);
free(theta);
free(rou);
free(u);
free(p);
free(Re);
free(lamta);
free(Nu);
free(H);
free(u_ever);
free(k);
free(rou1);
free(u1);
free(p1);
free(Re1);
free(lamta1);
free(Nu1);
free(H1);
free(u_ever1);
free(k1);
free(theta1);
}