?? setup2.c
字號:
#include "global_var.h"
void diflow()
{ double temp;
acof=diff;
if (flow==0) return;
temp=diff-fabs(flow)*0.1;
acof=0.0;
if(temp<0) return;
temp=temp/diff;
acof=diff*pow(temp,5.0);
return;
}
void setup2()
{ double rel,fl,flm,gm,gmm;
int i,j,n;
//u方程的aw(aim) ae(aip) an(ajp) as(ajm) 系數設定
nf=1;
if (lsolve[nf])
{
ist=3;
jst=2;
user_gamsor();
rel=1.0-relax[nf];
//以下 i的循環設定控制體u(3,2)到u(6,2)的as系數。見圖9
for (i=3;i<=l2;i++)
{
fl=xcvi[i]*v(i,2)*rho(i,1);
flm=xcvip[i-1]*v(i-1,2)*rho(i-1,1);
//FLOW=R(1)*(FL+FLM);
flow=(fl+flm);
//DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2)
diff=(xcvi[i]*gam(i,1)+xcvip[i-1]*gam(i-1,1))/ydif[2]; //ydif[2]=(δy)s
diflow();
ajm[i][2]=acof+__max(0.0,flow);
}
for (j=2;j<=m2;j++)
{
//以下語句設定u(3,j)控制體的aw(aim)系數，見圖10。
flow=arx[j]*u(2,j)*rho(1,j);
//DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J)) diff=arx[j]*gam(1,j)/(xcv[2]); //xcv[2]=(δx)w
diflow();
aim[3][j]=acof+__max(0.,flow);
for (i=3;i<=l2;i++)
{
if(i==l2)
{
flow=arx[j]*u(l1-1,j)*rho(l1-1,j);
//DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J))
diff=arx[j]*gam(l1,j)/(xcv[l2]);
}
else
//u(i,j)控制體的aw(aip)與u(i+1,j)控制體的ae(aim)系數。見圖11。
{
fl=u(i,j)*(fx[i]*rho(i,j)+fxm[i]*rho(i-1,j));
flp=u(i+1,j)*(fx[i+1]*rho(i+1,j)+fxm[i+1]*rho(i,j));
flow=arx[j]*0.5*(fl+flp);
//DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J))
diff=arx[j]*gam(i,j)/(xcv[i]); //xcv[i]=2*(δx)w? -=diff減半=- ?
}
diflow();
aim[i+1][j]=acof+__max(0.,flow); //u(i+1,j)的aw
aip[i][j]=aim[i+1][j]-flow; //u(i,j)ae
//*******************************************************************
//下面要計算的u控制體上s n面上的系數 ajm(as) ajp(an)
// if (j==m2) 要計算的是最右邊一列的u控制體的y方向的flow與diff，進而計算出ajm(as) ajp(an)
if (j==m2)
{
fl=xcvi[i]*v(i,m1)*rho(i,m1);
flm=xcvip[i]*v(i-1,m1)*rho(i-1,m1);
//DIFF=R(M1)*(XCVI[i]*GAM(I,M1)+xcvip[i-1]*GAM(I-1,M1))/YDIF(M1)
diff=(xcvi[i]*gam(i,m1)+xcvip[i-1]*gam(i-1,m1))/ydif[m1];
}
else
//u(i,j)控制體的an(ajp)與u(i,j+1)控制體的as(ajm)系數。flow計算見圖12。diff計算見圖13。
{
fl=xcvi[i]*v(i,j+1)*(fy[j+1]*rho(i,j+1)+fym[j+1]*rho(i,j));
flm=xcvip[i-1]*v(i-1,j+1)*(fy[j+1]*rho(i-1,j+1)+fym[j+1]*rho(i-1,j));
gm=gam(i,j)*gam(i,j+1)/(ycv[j]*gam(i,j+1)+ycv[j+1]*gam(i,j)+1.0e-30)*xcvi[i];
gmm=gam(i-1,j)*gam(i-1,j+1)/(ycv[j]*gam(i-1,j+1)+ycv[j+1]*gam(i-1,j)+1.e-30)*xcvip[i-1];
//diff=rmn(j+1)*2.*(gm+gmm);
diff=2.*(gm+gmm);
}
//flow=rmn(j+1)*(fl+flm);
flow=(fl+flm);
diflow();
ajm[i][j+1]=acof+__max(0.,flow);
ajp[i][j]=ajm[i][j+1]-flow;
vol=ycvr[j]*xcvs[i]; //ycvr:y方向上的寬度,xcvs與y垂直的面的面積
apt=(rho(i,j)*xcvi[i]+rho(i-1,j)*xcvip[i-1])/(xcvs[i]*dt); // apt=ap0/(ΔxΔy)=(ρ0/Δt)/(ΔxΔy) 見5.57e 不包括面積。
ap[i][j]=ap[i][j]-apt; // Sp-ρ/Δt. ap在user_gamsor中被賦值為Sp
con[i][j]=con[i][j]+apt*u(i,j); //con=b/(ΔxΔy)=Sc+(ap0*Φ0) con在user_gamsor中被賦值為Sc 5.57f
ap[i][j]=(-ap[i][j]*vol+aip[i][j]+aim[i][j]+ajp[i][j]+ajm[i][j])/relax[nf]; //5.57g 不過加上了松弛
con[i][j]=con[i][j]*vol+rel*ap[i][j]*u(i,j);
//DU(I,J)=VOL/(XDIF(I)*SX(J))
du[i][j]=vol/(xdif[i]); //面積Ae
con[i][j]=con[i][j]+du[i][j]*(p(i-1,j)-p(i,j)); //6.6
du[i][j]=du[i][j]/ap[i][j]; //du=dw=de 用于壓力校正方程 6.22與6.23
}//for i
}//for j
//u方程的aw ae an as系數設定－－結束
solve();
}// if lslove
//---v差分方程的系數設定----
nf=2;
if(lsolve[nf])
{
ist=2;
jst=3;
user_gamsor();
rel=1.-relax[nf];
for (i=2;i<=l2;i++) //圖14
{
//AREA=R(1)*XCV(I)
area=xcv[i];
flow=area*v(i,2)*rho(i,1);
diff=area*gam(i,1)/ycv[2]; //ycv[2]=(δy)s
diflow();
ajm[i][3]=acof+__max(0.,flow);
}
for (j=3;j<=m2;j++) //圖15
{
fl=arxj[j]*u(2,j)*rho(1,j);
flm=arxjp[j-1]*u(2,j-1)*rho(1,j-1);
flow=fl+flm;
//DIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J))
diff=(arxj[j]*gam(1,j)+arxjp[j-1]*gam(1,j-1))/(xdif[2]); //(xdif[2])=(δx)w
diflow();
aim[2][j]=acof+__max(0.,flow);
for (i=2;i<=l2;i++)
{
if(i==l2)
{
fl=arxj[j]*u(l1,j)*rho(l1,j);
flm=arxjp[j-1]*u(l1,j-1)*rho(l1,j-1);
//DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J))
diff=(arxj[j]*gam(l1,j)+arxjp[j-1]*gam(l1,j-1))/(xdif[l1]); //(xdif[l1-1])=(δx)w
}
else //圖16
{
fl=arxj[j]*u(i+1,j)*(fx[i+1]*rho(i+1,j)+fxm[i+1]*rho(i,j));
flm=arxjp[j-1]*u(i+1,j-1)*(fx[i+1]*rho(i+1,j-1)+fxm[i+1]*rho(i,j-1));
gm=gam(i,j)*gam(i+1,j)/(xcv[i]*gam(i+1,j)+xcv[i+1]*gam(i,j)+1.e-30)*arxj[j];
gmm=gam(i,j-1)*gam(i+1,j-1)/(xcv[i]*gam(i+1,j-1)+xcv[i+1]*gam(i,j-1)+1.0e-30)*arxjp[j-1];
//DIFF=2.*(GM+GMM)/SXMN(J)
diff=2.*(gm+gmm);
}
flow=fl+flm;
diflow();
aim[i+1][j]=acof+__max(0.,flow);
aip[i][j]=aim[i+1][j]-flow;
//aim aip計算完畢
//
if(j==m2)
{
//AREA=R(M1)*XCV(I)
area=xcv[i];
flow=area*v(i,m1)*rho(i,m1);
diff=area*gam(i,m1)/ycv[m2]; //ycv(m2)=(δy)n
}
else //圖17
{
//area=r[j]*xcv[i];
area=xcv[i];
//fl=v(i,j)*(fy[j]*rho(i,j)+fym[j]*rho(i,j-1))*rmn[j];
//flp=v(i,j+1)*(fy[j+1]*rho(i,j+1)+fym[j+1]*rho(i,j))*rmn(j+1);
fl=v(i,j)*(fy[j]*rho(i,j)+fym[j]*rho(i,j-1));
flp=v(i,j+1)*(fy[j+1]*rho(i,j+1)+fym[j+1]*rho(i,j));
flow=(fv[j]*fl+fvp[j]*flp)*xcv[i];
diff=area*gam(i,j)/ycv[j]; // ycv[j]=2(δy)n? -=diff減半=-?
}
diflow();
ajm[i][j+1]=acof+__max(0.,flow);
ajp[i][j]=ajm[i][j+1]-flow;
vol=ycvrs[j]*xcv[i];
//SXT=SX(J)
//IF(J .EQ. M2) SXT=SX(M1)
//SXB=SX(J-1)
//IF(J .EQ. 3) SXB=SX(1)
//APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)*RHO(I,J-1)*0.5*(SXB+SXMN(J)))/(YCVRS(J)*DT) apt=(arxj[j]*rho(i,j)+arxjp[j-1]*rho(i,j-1))/(ycvrs[j]*dt);
ap[i][j]=ap[i][j]-apt;
con[i][j]=con[i][j]+apt*v(i,j);
ap[i][j]=(-ap[i][j]*vol+aip[i][j]+aim[i][j]+ajp[i][j]+ajm[i][j])/relax[nf];
con[i][j]=con[i][j]*vol+rel*ap[i][j]*v(i,j);
dv[i][j]=vol/ydif[j];
con[i][j]=con[i][j]+dv[i][j]*(p(i,j-1)-p(i,j));
dv[i][j]=dv[i][j]/ap[i][j]; // dv=ds=dn 用于壓力校正方程 6.22與6.23
}//for i
}//for j
solve();
} //if lsolve
//*---壓力校正方程系數設定----
nf=3;
if(lsolve[nf])
{
ist=2;
jst=2;
user_gamsor();
smax=0.;
ssum=0.;
for (j=2;j<=m2;j++)
for (i=2;i<=l2;i++)
{
vol=ycvr[j]*xcv[i];
con[i][j]=con[i][j]*vol;
}
for (i=2;i<=l2;i++)
{
//ARHO=R(1)*XCV(I)*RHO(I,1)
arho=xcv[i]*rho(i,1);
con[i][2]=con[i][2]+arho*v(i,2);
ajm[i][2]=0.;
}
for (j=2;j<=m2;j++)
{
arho=arx[j]*rho(1,j);
con[2][j]=con[2][j]+arho*u(2,j);
aim[2][j]=0.;
for(i=2;i<=l2;i++)
{
if(i==l2)
{
arho=arx[j]*rho(l1,j);
con[i][j]=con[i][j]-arho*u(l1,j);
aip[i][j]=0.;
}
else
{
arho=arx[j]*(fx[i+1]*rho(i+1,j)+fxm[i+1]*rho(i,j));
flow=arho*u(i+1,j);
con[i][j]=con[i][j]-flow;
con[i+1][j]=con[i+1][j]+flow;
aip[i][j]=arho*du[i+1][j];
aim[i+1][j]=aip[i][j];
}
if(j==m2)
{
//ARHO=RMN(M1)*XCV(I)*RHO(I,M1)
arho=xcv[i]*rho(i,m1);
con[i][j]=con[i][j]-arho*v(i,m1);
ajp[i][j]=0.;
}
else
{
//ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
arho=xcv[i]*(fy[j+1]*rho(i,j+1)+fym[j+1]*rho(i,j));
flow=arho*v(i,j+1);
con[i][j]=con[i][j]-flow;
con[i][j+1]=con[i][j+1]+flow;
ajp[i][j]=arho*dv[i][j+1];
ajm[i][j+1]=ajp[i][j];
}
ap[i][j]=aip[i][j]+aim[i][j]+ajp[i][j]+ajm[i][j];
pc(i,j)=0.;
smax=__max(smax,fabs(con[i][j]));
ssum=ssum+con[i][j];
} //for i
}//for j
solve();
//---對壓力場與速度場進行校正----
for (j=2;j<=m2;j++)
for (i=2;i<=l2;i++)
{
p(i,j)=p(i,j)+pc(i,j)*relax[np];
if(i!=2) u(i,j)=u(i,j)+du[i][j]*(pc(i-1,j)-pc(i,j));
if(j!=2) v(i,j)=v(i,j)+dv[i][j]*(pc(i,j-1)-pc(i,j));
}
} //lsolve
//---其它通用變量φ方程系數設定----
ist=2;
jst=2;
for (n=4;n<=nfmax;n++)
{
nf=n;
if(!lsolve[nf]) continue;
user_gamsor();
rel=1.-relax[nf];
for (i=2;i<=l2;i++) //圖18
{
//AREA=R(1)*XCV(I)
area=xcv[i];
flow=area*v(i,2)*rho(i,1);
diff=area*gam(i,1)/ydif[2];
diflow();
ajm[i][2]=acof+__max(0.,flow);
} //for i
for (j=2;j<=m2;j++) //圖19
{
flow=arx[j]*u(2,j)*rho(1,j);
//DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J))
diff=arx[j]*gam(1,j)/(xdif[2]);
diflow();
aim[2][j]=acof+__max(0.,flow);
for (i=2;i<=l2;i++)
{
if(i==l2)
{
flow=arx[j]*u(l1,j)*rho(l1,j);
//DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J)) diff=arx[j]*gam(l1,j)/(xdif[l1]);
}
else //圖20
{
flow=arx[j]*u(i+1,j)*(fx[i+1]*rho(i+1,j)+fxm[i+1]*rho(i,j));
//diff=arx[j]*2.*gam(i,j)*gam(i+1,j)/((xcv[i]*gam(i+1,j)+xcv(i+1)*gam(i,j)+1.0e-30)*sx[j]);
diff=arx[j]*2.*gam(i,j)*gam(i+1,j)/(xcv[i]*gam(i+1,j)+xcv[i+1]*gam(i,j)+1.0e-30);
}
diflow();
aim[i+1][j]=acof+__max(0.,flow);
aip[i][j]=aim[i+1][j]-flow;
//area=rmn(j+1)*xcv[i];
area=xcv[i];
if(j==m2)
{
flow=area*v(i,m1)*rho(i,m1);
diff=area*gam(i,m1)/ydif[m1];
}
else //圖21
{
flow=area*v(i,j+1)*(fy[j+1]*rho(i,j+1)+fym[j+1]*rho(i,j));
diff=area*2.*gam(i,j)*gam(i,j+1)/(ycv[j]*gam(i,j+1)+ycv[j+1]*gam(i,j)+1.0e-30);
}
diflow();
ajm[i][j+1]=acof+__max(0.,flow);
ajp[i][j]=ajm[i][j+1]-flow;
vol=ycvr[j]*xcv[i];
apt=rho(i,j)/dt;
ap[i][j]=ap[i][j]-apt;
con[i][j]=con[i][j]+apt*f[i][j][nf];
ap[i][j]=(-ap[i][j]*vol+aip[i][j]+aim[i][j]+ajp[i][j]+ajm[i][j])/relax[nf];
con[i][j]=con[i][j]*vol+rel*ap[i][j]*f[i][j][nf];
} //for i
} //for j
solve();
} //for n
return;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -