//Libraries reading # include # include # include # include # include #include /* ********************************************************************************************************** * THIS PROGRAM SOLVES NUMERICALLY THE BURGERS EQUATION FOR TEMPERATURE-INHOMOGENEOUS FLUID AS HOMOGENEOUS* * ********************************************************************************************************** */ /* Declaration of global variables and functions */ int M; double st,alfa,beta,cA,cB,p0; double w,g,R,TA,rA,L,Q,G,N,A; double s,n_steps,ss; double complex Pn[201]; /* Field for complex parts of the pressure - set by the norm*/ double complex TPn[201]; /* Temporary field for complex parts of the pressure - set by the norm*/ double abso(double q); /* Absolute value function */ double T0(double x); /* Temperature function */ double dT0(double x); /* Derivation of temperature function */ double c0(double x); /* Sound velocity function */ double r0(double x); /* Density function */ /* Complex right-hand side of differential equation- set by the norm*/ double complex Fn(int n, double x); const float pi = 3.1415926537; int main() { /* begin */ /* Declaration of local variables */ FILE *cR, *cI, *Data; int i; double complex k1[201]; double complex k2[201]; double complex k3[201]; double complex k4[201]; /* -------------------------------------------------------------------------*/ /* --------------------Program body starts here--------------------------- */ /* -------------------------------------------------------------------------*/ /* Set ZEROs to the fields */ for (i = 0; i <= 200; i++) { Pn[i] = 0.+0.*I; TPn[i] = 0.+0.*I; } /* Setting of input parameters */ p0 = 101325; R = 287.058; g = 1.402; L = 1.; TA = 298.15;//357.78; cA = 346.4; rA = 1.193; alfa = 0.00004578; beta = 1.201; w = 10000; A = 0.5; /* Boundary condition */ Pn[1] = -(1000/(rA*cA*cA))*I; //Pn[1] = -0.5*I; s = 0.; /* Setting of dimensionless distance */ //Setting of input parameters printf("Set the whole distance ss="); scanf("%lf", &ss); printf("Set the whole number of Fourier series members (max. 200) M="); scanf("%d", &M); printf("Set the number of steps="); scanf("%lf", &n_steps); /* the distance of the numerical step */ st = ss/n_steps; //Linear dependence of temperature on the distance: T0=TA*(1+a*x) cB = c0(L); //Parameters calculation Q = w*L/(2*cA); N = beta*w*L/(cA); G = alfa*w*w*L/(2*rA*cA*cA*cA); /* Pn = cR/2 - j*cI/2 */ Data = fopen("Data.txt", "wt"); /* Information about resonator and the calculation is been writing here */ fprintf(Data, "The number of Fourier series members M= %d \n",M); fprintf(Data, "The number of steps= %lf \n",n_steps); fprintf(Data, "step= %lf \n",st); fprintf(Data, "The whole distance ss=%.8lf \n",ss); fprintf(Data, "Q=%.8lf \n",Q); fprintf(Data, "G=%.10lf \n",G); fprintf(Data, "N=%.8lf \n",N); fprintf(Data, "A=%.8lf \n",A); fprintf(Data, "cA=%.8lf \n",cA); fprintf(Data, "cB=%.8lf \n",cB); fclose(Data); //Eigenvalue calculation starts here for (i = 1; i<=M; i++) { TPn[i] = Pn[i]; } while (s<=ss) { /*begin while */ //RK4 method for (i = 1; i<=M; i++) { k1[i] = Fn(i,s); } for (i = 1; i<=M; i++) { TPn[i] = Pn[i]+k1[i]/2.; } for (i = 1; i<=M; i++) { k2[i] = Fn(i,s+st/2.); } for (i = 1; i<=M; i++) { TPn[i] = Pn[i]+k2[i]/2.; } for (i = 1; i<=M; i++) { k3[i] = Fn(i,s+st/2.); } for (i = 1; i<=M; i++) { TPn[i] = Pn[i]+k3[i]; } for (i = 1; i<=M; i++) { k4[i] = Fn(i,s+st); } for (i = 1; i<=M; i++) { TPn[i] = Pn[i]+(1./6.)*(k1[i]+2.*k2[i]+2.*k3[i]+k4[i])*sin(i/25.)/(i/25.); } for (i = 1; i<=M; i++) { Pn[i] = TPn[i]; } printf(" s= %lf \n",s);//writing the distance to the display s = s+st; //distance incrementation } /* end while */ cR = fopen("cR.txt", "wt"); /* Real parts of the pressure harm. components are been writing here*/ cI = fopen("cI.txt", "wt"); /* Imaginary parts of the pressure harm. components are been writing here */ for (i = 1; i<=M; i++) if (i==M) { fprintf(cR, "%.8e \n",creal(Pn[i])); fprintf(cI, "%.8e \n",cimag(Pn[i])); } else { fprintf(cR, "%.8e \t",creal(Pn[i])); fprintf(cI, "%.8e \t",cimag(Pn[i])); } fclose(cR);/* finish of the text file for the real Fourier series members */ fclose(cI);/* finish of the text file for the imaginary Fourier series members */ } /* end */ /* Functions definitions */ double T0(double x) { double TT; TT = TA*(1+A*x); return TT; } double dT0(double x) { double ddT; ddT = TA*A*x; return ddT; } double c0(double x) { double cc; cc = sqrt(g*R*T0(x)); return cc; } double r0(double x) { double rr; rr = p0/(R*T0(x)); return rr; } double complex Fn(int n, double x) { int m; double complex konvoluce,F; konvoluce = 0.+0.*I; for (m = 1; m < n ; m++) konvoluce += TPn[m]*TPn[n-m]; for (m = (n+1); m <= M; m++) konvoluce += 2*(TPn[m]*conj(TPn[m-n])); // konvoluce = 0.+0.*I; //Solve linear equation by comment out F = I*(n*N*TA/(2.*T0(x)))*konvoluce -(n*n*G*TA/T0(x))*TPn[n] -dT0(x)*TPn[n]/(2*T0(x)) -I*n*Q*((TA-T0(x))/T0(x))*TPn[n]; F = st*F; return F; } double abso(double q) { double ab; ab = q; if (ab < 0) ab = -ab; return ab; }