📄 exam10.c
字号:
#include "engine.h"
#include "matrix.h"
#include "global_var.h"
void user_could_reset_some_setting_here()
{
int i;
char *string1=" temprature ";
char *string2=" vel u ";
char *string3=" str fn ";
char *string4=" vel v ";
char *string5=" pressure ";
char *string6=" density ";
char *string7=" k.energy ";
char *string8=" disipat ";
char *string9=" turbulent visc ";
title[1]=strdup(string2);
title[2]=strdup(string4);
title[3]=strdup(string3);
title[4]=strdup(string1);
title[5]=strdup(string7);
title[6]=strdup(string8);
title[7]=strdup(string9);
title[11]=strdup(string5);
title[12]=strdup(string6);
relax[1]=0.5;
relax[2]=0.5;
relax[11]=0.8;
relax[5]=0.4;
relax[6]=0.4;
relax[13]=0.5;
for (i=1;i<=6;i++)
{
lsolve[i]=1;
lprint[i]=1;
}
lprint[11]=1;
lprint[7]=1;
lprintdebuginfo=0;
mode=1;
last=1;
dt=1e30;
xl=1.0;
yl=4.0;
l1=59;
m1=59;
amu=1.e-6;
cmu=.09;
c1=1.44;
c2=1.92;
prt=.9;
prk=1.;
prd=1.3;
pr=.7;
prprt=pr/prt;
pfn=9.*(prprt-1.)/pow(prprt,0.25);
cmu4=pow(cmu,0.25);
for (i=1;i<=nfmax;i++) ntimes[i]=30;
}
void user_generate_velocity_grid() //速度u v的控制网格
{ int i,j;
double dx,dy;
xu[2]=0; // 第二个主控制容积在x方向的w界面位置，即速度u(i,j)所在位置的x坐标＝0
// 难道说，第一个控制容积留给了虚拟的边壁？ ???
dx=xl/(double)(l1-2); // x方向上的网格宽度
for (i=3;i<=l1;i++) // l1-1就是最后一个控制体（虚拟边界控制体）的i坐标了。
//计数的习惯与C的数组的序列号保持一致，从0开始。
{
xu[i]=xu[i-1]+dx; //网格在x方向是均匀，所以第三列～第七列控制容积的w面简单以dx递增
}
// NOTE : xu[1] 没有被赋值
yv[2]=0;
dy=yl/(double)(m1-2);
for (j=3;j<=m1;j++)
{
yv[j]=yv[j-1]+dy;
}
// NOTE :xu[1] yv[1] 没有被赋值
}
void user_initial_grid_value()
{
int i,j;
for (j=1;j<=m1;j++)
for (i=1;i<=l1;i++)
{
u(i,j)=0;
v(i,j)=5;
v(1,j)=0;
v(i,2)=5;
if (i>(l1/2)) v(i,2)=10;
t(i,j)=100;
t(1,j)=0;
if (i>(l1/2)) t(i,1)=400;
ake(i,j)=0.005*v(i,2)*v(i,2);
dis(i,j)=0.1*ake(i,j)*ake(i,j);
}
}
void user_correct_density()
{
}
void user_bound_condition_setup()
{
static double flowin;
double factor,fl,afl,vmin;
int i,j;
if (iter1==0) //每一个时间片开始后，重新计算入口流量。因为入口条件可能会变
{
flowin=0;
for (i=2;i<=l2;i++) flowin+=rho(i,1)*v(i,2)*xcv[i];
}
fl=0;
afl=0.;
vmin=0.;
for (i=2;i<=l2;i++)
{
if(v(i,m2)<0) vmin=__max(vmin,-v(i,m2));
afl=afl+rho(i,m1)*xcv[i];
fl=fl+rho(i,m1)*v(i,m2)*xcv[i];
}
factor=flowin/fl;
for (i=2;i<=l2;i++) v(i,m1)=(v(i,m2)+vmin)*factor;
for (j=2;j<=m2;j++) dis(l1,j)=dis(l2,j);
}
void user_gamsor()
{
double rel,amt,factor,xplus[MAXGRIDSIZE],gamp,gamm,dudx,dvdy,dudy,dvdx,diss;
int i,j;
if(nf==3) return;
if(nf==1)
{
rel=1.-relax[ngam]; //relax[13]=0.5
for (j=1;j<=m1;j++)
for (i=1;i<=l1;i++)
{
amt=cmu*rho(i,j)*ake(i,j)*ake(i,j)/(dis(i,j)+1.e-30);
if(iter1==0&&iter==0) amut(i,j)=amt;
amut(i,j)=relax[ngam]*amt+rel*amut(i,j);
}
}
factor=1.;
if(nf==4) factor=1./prt;
if(nf==5) factor=1./prk;
if(nf==6) factor=1./prd;
for (j=1;j<=m1;j++)
for (i=1;i<=l1;i++)
{
gam(i,j)=amut(i,j)*factor; //amut在nf==1时计算后一直保留
if (nf!=1) gam(l1,j)=0.;
gam(i,m1)=0;
}
for (j=2;j<=m2;j++)
{
switch(nf)
{
case 1: //u方程
gam(1,j)=0;
break;
case 2: //v方程
xplus[j]=rho(2,j)*sqrt(ake(2,j))*cmu4*xdif[2]/amu;
gam(1,j)=amu;
if(xplus[j]>11.5) gam(1,j)=amu*xplus[j]/(log(9.*xplus[j])*2.5);
break;
case 3: //压力校正
gam(1,j)=0;
break;
case 4: //温度
gam(1,j)=amu/pr;
if(xplus[j]>11.5) gam(1,j)=amu/prt*xplus[j]/(2.5*log(9.*xplus[j])+pfn);
break;
case 5: //k
gam(1,j)=0;
break;
case 6: //diss
gam(1,j)=0;
break;
}
}
switch(nf)
{
case 1:
{
for (j=2;j<=m2;j++)
for (i=3;i<=l2;i++)
{
con[i][j]=(gam(i,j)*(u(i+1,j)-u(i,j))/xcv[i]-gam(i-1,j)*(u(i,j)-u(i-1,j))/xcv[i-1])/xdif[i];
gamp=gam(i,j+1)*gam(i-1,j+1)/(gam(i,j+1)+gam(i-1,j+1)+1.e-30);
gamp=gamp+gam(i,j)*gam(i-1,j)/(gam(i,j)+gam(i-1,j)+1.e-30);
gamm=gam(i,j-1)*gam(i-1,j-1)/(gam(i,j-1)+gam(i-1,j-1)+1.e-30);
gamm=gamm+gam(i,j)*gam(i-1,j)/(gam(i,j)+gam(i-1,j)+1.e-30);
con[i][j]=con[i][j]+(gamp*(v(i,j+1)-v(i-1,j+1))-gamm*(v(i,j)-v(i-1,j)))/(ycv[j]*xdif[i]);
ap[i][j]=0;
}
break;
}
case 2:
{
for (j=3;j<=m2;j++)
for (i=2;i<=l2;i++)
{
con[i][j]=(gam(i,j)*(v(i,j+1)-v(i,j))/ycv[j]-gam(i,j-1)*(v(i,j)-v(i,j-1))/ycv[j-1])/ydif[j];
gamp=gam(i+1,j)*gam(i+1,j-1)/(gam(i+1,j)+gam(i+1,j-1)+1.e-30);
gamp=gamp+gam(i,j)*gam(i,j-1)/(gam(i,j)+gam(i,j-1)+1.e-30);
gamm=gam(i-1,j)*gam(i-1,j-1)/(gam(i-1,j)+gam(i-1,j-1)+1.e-30);
gamm=gamm+gam(i,j)*gam(i,j-1)/(gam(i,j)+gam(i,j-1)+1.e-30);
con[i][j]=con[i][j]+(gamp*(u(i+1,j)-u(i+1,j-1))-gamm*(u(i,j)-u(i,j-1)))/(xcv[i]*ydif[j]);
ap[i][j]=0.;
}
break;
}
case 4:
{
for (j=2;j<=m2;j++)
for (i=2;i<=l2;i++)
{
con[i][j]=0.;
ap[i][j]=0.;
}
break;
}
case 5:
{
for (j=2;j<=m2;j++)
for (i=2;i<=l2;i++)
{
dudx=(u(i+1,j)-u(i,j))/xcv[i];
dvdy=(v(i,j+1)-v(i,j))/ycv[j];
if(j==2) dudy=(0.5*(u(i,j+1)-u(i,j))+0.5*(u(i+1,j+1)-u(i+1,j)))/(ydif[j+1]);
else if(j==m2) dudy=(0.5*(u(i,j)-u(i,j-1))+0.5*(u(i+1,j)-u(i+1,j-1)))/ydif[j];
else dudy=(0.5*(u(i,j+1)-u(i,j-1))+0.5*(u(i+1,j+1)-u(i+1,j-1)))/(ydif[j]+ydif[j+1]);
if(i==2) dvdx=(0.5*(v(i+1,j)-v(i,j))+0.5*(v(i+1,j+1)-v(i,j+1)))/xdif[i+1];
else if(i==l2) dvdx=(0.5*(v(i,j)-v(i-1,j))+0.5*(v(i,j+1)-v(i-1,j+1)))/xdif[i];
else dvdx=(0.5*(v(i+1,j)-v(i-1,j))+0.5*(v(i+1,j+1)-v(i-1,j+1)))/(xdif[i]+xdif[i+1]);
gen(i,j)=2.*(dudx*dudx+dvdy*dvdy)+(dudy+dvdx)*(dudy+dvdx);
con[i][j]=gen(i,j)*amut(i,j);
ap[i][j]=-rho(i,j)*dis(i,j)/(ake(i,j)+1.e-30);
}
break;
}
case 6:
{
for (j=2;j<=m2;j++)
for (i=2;i<=l2;i++)
{
con[i][j]=c1*gen(i,j)*cmu*rho(i,j)*ake(i,j);
ap[i][j]=-c2*rho(i,j)*dis(i,j)/(ake(i,j)+1.e-30);
}
for (j=2;j<m2;j++)
{
diss=cmu*pow(ake(2,j),1.5)/(0.4*cmu4*xdif[2]);
con[2][j]=1.e+30*diss;
ap[2][j]=0.-1.e+30;
}
break;
}
} //switch nf
}
void user_output()
{
//fprintf(filehandle8,"\niter=%3d SMAX=%7.3f SSUM=%9.4f v(4,7)=%9.4f T(4,7)=%9.4f F(4,7,3)=%9.4f\n",iter,smax,ssum,v(4,7),t(4,7),f[4][7][3]);
fprintf(filehandle8,"\niter==%3d time=%7.3f\n",iter,time);
printinfo();
user_plot();
}
void give_matlab_a_double_matrix(Engine* ep, char * matrix_name, int matrix_x_dimension, int matrix_y_dimension, double * databuffer)
{ mxArray *T;
T=mxCreateDoubleMatrix(matrix_x_dimension,matrix_y_dimension,mxREAL);
mxSetName(T, matrix_name);
memcpy((void*)mxGetPr(T),(void*)databuffer,matrix_x_dimension*matrix_y_dimension*sizeof(double));
engPutArray(ep, T);
mxDestroyArray(T);
}
void give_matlab_a_interger(Engine* ep, char * matrix_name, int data)
{ char matlabsentence[255];
sprintf(matlabsentence,"%s=%d;",matrix_name,data);
engEvalString(ep,matlabsentence);
}
void user_plot() //MATLAB画图
{
int i,j,matrixcount,velocv_ibeg,velocv_iend,velocv_jbeg,velocv_jend,maincv_ibeg,maincv_iend,maincv_jbeg,maincv_jend;
double matrixX[MAXGRIDSIZE*MAXGRIDSIZE],matrixY[MAXGRIDSIZE*MAXGRIDSIZE],matrixU[MAXGRIDSIZE*MAXGRIDSIZE];
double matrixV[MAXGRIDSIZE*MAXGRIDSIZE],matrixP[MAXGRIDSIZE*MAXGRIDSIZE],matrixT[MAXGRIDSIZE*MAXGRIDSIZE],matrixGAM[MAXGRIDSIZE*MAXGRIDSIZE];
static Engine *ep;
velocv_ibeg=2;
velocv_iend=l1-1;
velocv_jbeg=2;
velocv_jend=m1-1;
maincv_ibeg=1;
maincv_iend=l1-1;
maincv_jbeg=1;
maincv_jend=m1-1;
if (iter==0)
{
ep = engOpen(NULL);
give_matlab_a_interger(ep,"last",last);
give_matlab_a_interger(ep,"iter",iter);
engEvalString(ep,"M=moviein(last+1);iter=1");
}
matrixcount=0;
for (j=velocv_jbeg;j<=velocv_jend;j++)
for (i=velocv_ibeg;i<=velocv_iend;i++)
{
matrixX[matrixcount]=xu[i];
matrixY[matrixcount]=yv[j];
matrixU[matrixcount]=u(i,j);
matrixV[matrixcount]=v(i,j);
matrixcount++;
}
give_matlab_a_double_matrix(ep,"VELOX",velocv_iend-velocv_ibeg+1,velocv_jend-velocv_jbeg+1,matrixX);
give_matlab_a_double_matrix(ep,"VELOY",velocv_iend-velocv_ibeg+1,velocv_jend-velocv_jbeg+1,matrixY);
give_matlab_a_double_matrix(ep,"U",velocv_iend-velocv_ibeg+1,velocv_jend-velocv_jbeg+1,matrixU);
give_matlab_a_double_matrix(ep,"V",velocv_iend-velocv_ibeg+1,velocv_jend-velocv_jbeg+1,matrixV);
matrixcount=0;
for (j=maincv_jbeg;j<=maincv_jend;j++)
for (i=maincv_ibeg;i<=maincv_iend;i++)
{
matrixX[matrixcount]=x[i];
matrixY[matrixcount]=y[j];
matrixT[matrixcount]=t(i,j);
matrixP[matrixcount]=p(i,j);
matrixGAM[matrixcount]=gam(i,j);
matrixcount++;
}
give_matlab_a_double_matrix(ep,"X",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixX);
give_matlab_a_double_matrix(ep,"Y",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixY);
give_matlab_a_double_matrix(ep,"T",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixT);
give_matlab_a_double_matrix(ep,"P",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixP);
give_matlab_a_double_matrix(ep,"GAM",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixGAM);
if (iter==last-1)
{
engEvalString(ep,"hold off;pcolor(X,Y,T);shading interp;caxis([100 400]);colorbar;axis equal;hold on;quiver(VELOX,VELOY,U,V);");
}
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -