?? 555.cpp.bak
字號:
gv[0][JMAX-1]=-gv[3][JMAX-3];
gv[1][JMAX-1]=-gv[2][JMAX-3];
/* right top */
gv[IMAX-1][JMAX-1]=-gv[IMAX-4][JMAX-3];
gv[IMAX-2][JMAX-1]=-gv[IMAX-3][JMAX-3];
return(0);
}
int bc_velo_sc(double sc[IMAX][JMAX], double dvalue){
/* boundary condition for the velocity at Scaler point */
int i,j;
/* LEFTSIDE */
/*... outlet ...*/
for(j=2; j<=OUTLET+1; j++){
sc[0][j]=sc[2][j];
sc[1][j]=sc[2][j];
}
/*... wall ...*/
for(j=OUTLET+2; j<=JMAX-3; j++){
sc[0 ][j]=-sc[3 ][j];
sc[1 ][j]=-sc[2 ][j];
}
/* RIGHTSIDE */
/*... inlet & open boundary ...*/
for(j=2; j<=JMAX-3; j++){
sc[IMAX-1][j]=sc[IMAX-3][j];
sc[IMAX-2][j]=sc[IMAX-3][j];
}
/* BOTTOM */
for(i=2; i<=IMAX-3; i++){
/* bottom */
sc[i][0]=-sc[i][3];
sc[i][1]=-sc[i][2];
}
/* UPPER */
for(i=2; i<=IMAX-3; i++){
/* upper */
sc[i][JMAX-1]=sc[i][JMAX-3];
sc[i][JMAX-2]=sc[i][JMAX-3];
}
/* corner */
/* left & bottom */
sc[0][0]=sc[3][3];
sc[0][1]=sc[3][2];
sc[1][0]=sc[2][3];
sc[1][1]=sc[2][2];
/* right & bottom */
sc[IMAX-1][0]=sc[IMAX-4][3];
sc[IMAX-1][1]=sc[IMAX-4][2];
sc[IMAX-2][0]=sc[IMAX-3][3];
sc[IMAX-2][1]=sc[IMAX-3][2];
/* left & top (wall condition) */
sc[0][JMAX-1]=sc[3][JMAX-4];
sc[0][JMAX-2]=sc[3][JMAX-3];
sc[1][JMAX-1]=sc[2][JMAX-4];
sc[1][JMAX-2]=sc[2][JMAX-3];
/* right & top (wall condition)*/
sc[IMAX-1][JMAX-1]=sc[IMAX-4][JMAX-4];
sc[IMAX-1][JMAX-2]=sc[IMAX-4][JMAX-3];
sc[IMAX-2][JMAX-1]=sc[IMAX-3][JMAX-4];
sc[IMAX-2][JMAX-2]=sc[IMAX-3][JMAX-3];
return(0);
}
/******** YOU SHOULD ALSO CONFIRM BOUNDARY CONDITION FOR PRESSURE!! ’int bc_p()’ **********/
/******** YOU SHOULD ALSO CONFIRM BOUNDARY CONDITION FOR PRESSURE!! ’int bc_p()’ **********/
/******** YOU SHOULD ALSO CONFIRM BOUNDARY CONDITION FOR PRESSURE!! ’int bc_p()’ **********/
//---------- Boundary Condition ----------
int correct_phi_mizutani(double phin[IMAX][JMAX], double un[IMAX][JMAX], double T){
/* mizutani2002 conservative type */
int i,j,m=0,negative=0,mix[IMAX][JMAX];
double phi_n, phi_s, phi_w, phi_e;
double phi_se, phi_sw, phi_ne, phi_nw;
double vinitial=0., vcurrent=0., vcorrect=0., vcalc=0.;
double correctv=0.;
double in_u=0.;
/* initial volume */
vinitial=((IMAX-4)*INLET)+T*INLET*(-IN_U/DX)-T*OUTLET*(-OUT_U/DX);
// vinitial=((IMAX-4)*INLET)+T/DX*in_u-T/DX*OUTLET*(-OUT_U);
for(i=2; i<=IMAX-3; i++){
for(j=2; j<=JMAX-3; j++){
vcalc+=iphi(phin[i][j]);
// printf("(%3d,%3d) phin:%8.4f iphi(phin):%8.4f\n",i,j,phin[i][j],iphi(phin[i][j]));
}
}
printf("Vtheory=%8.4f Vcalc=%8.4f Vc/Vt=%6.2f%%\n",vinitial,vcalc,fabs(vcalc/vinitial*100.));
for(i=2; i<=IMAX-3; i++){
for(j=2; j<=JMAX-3; j++){
/*
+------+------+------+
im,jp i,jp ip,jp
+------NW-----NE-----+
im,j | i,j | ip,j
+------SW-----SE-----+
im,jm i,jm ip,jm
+------+------+------+
*/
phi_se=0.25*( iphi(phin[i][j]) + iphi(phin[i+1][j]) + iphi(phin[i+1][j-1]) + iphi(phin[i][j-1]) );
phi_ne=0.25*( iphi(phin[i][j]) + iphi(phin[i+1][j]) + iphi(phin[i+1][j+1]) + iphi(phin[i][j+1]) );
phi_sw=0.25*( iphi(phin[i][j]) + iphi(phin[i-1][j]) + iphi(phin[i-1][j-1]) + iphi(phin[i][j-1]) );
phi_nw=0.25*( iphi(phin[i][j]) + iphi(phin[i-1][j]) + iphi(phin[i-1][j+1]) + iphi(phin[i][j+1]) );
/* definition of mix.(air and water) cell */
if(
iphi(phin[i][j])==0.5 ||
(phi_sw< 0.5 && phi_ne> 0.5) ||
(phi_sw> 0.5 && phi_ne< 0.5) ||
(phi_se< 0.5 && phi_nw> 0.5) ||
(phi_se> 0.5 && phi_nw< 0.5)
){
mix[i][j]=1;
m++; /* count of mix. cells */
}
else{
if (iphi(phin[i][j])>0.5) phin[i][j]=dphi(1);
else if(iphi(phin[i][j])<0.5) phin[i][j]=dphi(0);
mix[i][j]=0;
}
}
}
if(m==0){
printf("no correction..... really?? \n");
return(0);
}
for(i=2; i<=IMAX-3; i++){
for(j=2; j<=JMAX-3; j++){
/* current volume */
vcurrent+=iphi(phin[i][j]);
}
}
correctv=(vinitial-vcurrent)/(double)m;
for(i=2; i<=IMAX-3; i++){
for(j=2; j<=JMAX-3; j++){
if(mix[i][j]==1) phin[i][j]=phin[i][j]+dphi(correctv);
vcorrect+=iphi(phin[i][j]);
if( iphi(phin[i][j])<0. || iphi(phin[i][j])>1.){
negative=1;
}
}
}
if(negative==1){
printf("WARNING(correct_phi_mizutani).....correct_phi(phin) was executed....\n");
correct_phi(phin);
}
bc_scalar(phin,dphi(1.));
return(0);
}
int new_gradient(double pvn[IMAX][JMAX], double pv[IMAX][JMAX], double gx[IMAX][JMAX], double gy[IMAX][JMAX]){
int i, j;
for(i=1; i<=IMAX-2; i++){
for(j=1; j<=JMAX-2; j++){
gx[i][j]=gx[i][j]+(pvn[i+1][j ]-pvn[i-1][j ]-pv[i+1][j ]+pv[i-1][j ])*0.5/DX;
gy[i][j]=gy[i][j]+(pvn[i ][j+1]-pvn[i ][j-1]-pv[i ][j+1]+pv[i ][j-1])*0.5/DY;
}
}
return(0);
}
void bodyforce(double w[IMAX][JMAX], double wn[IMAX][JMAX]){
int i,j;
for(i=2; i<=IMAX-3; i++){
for(j=3; j<JMAX; j++){
wn[i][j]=w[i][j]+(-G)*DT;
}
}
}
int output(double p[IMAX][JMAX], double rho[IMAX][JMAX], double u[IMAX][JMAX], double v[IMAX][JMAX], double phi[IMAX][JMAX],
int count, double time){
// for microAVS
FILE *fpfld, *fpdat;
int i, j;
int stp;
char avsfld[25]= "./Data/avs????.fld";
char avsdat[25]= "./Data/avs????.dat";
char avsfldl[25]= "./avs????.fld";
char avsdatl[25]= "./avs????.dat";
double us[IMAX][JMAX],vs[IMAX][JMAX],uv[IMAX][JMAX],vu[IMAX][JMAX];
stp=count/OUTPUT;
avsfld[10] = '0' + stp/1000;
avsfld[11] = '0' + (stp%1000)/100;
avsfld[12] = '0' + ((stp%1000)%100)/10;
avsfld[13] = '0' + ((stp%1000)%100)%10;
avsdat[10] = '0' + stp/1000;
avsdat[11] = '0' + (stp%1000)/100;
avsdat[12] = '0' + ((stp%1000)%100)/10;
avsdat[13] = '0' + ((stp%1000)%100)%10;
avsdatl[5] = '0' + stp/1000;
avsdatl[6] = '0' + (stp%1000)/100;
avsdatl[7] = '0' + ((stp%1000)%100)/10;
avsdatl[8] = '0' + ((stp%1000)%100)%10;
avsfldl[5] = '0' + stp/1000;
avsfldl[6] = '0' + (stp%1000)/100;
avsfldl[7] = '0' + ((stp%1000)%100)/10;
avsfldl[8] = '0' + ((stp%1000)%100)%10;
if ((fpfld=fopen (avsfld,"w+")) == NULL){
printf("%s was not opened!\n",avsfld);
printf("’Maybe you should make ’Data’ directory....\n");
exit(1);
}
fprintf(fpfld,"# AVS field file\n");
fprintf(fpfld,"#\n");
fprintf(fpfld,"ndim = 2\n");
fprintf(fpfld,"dim1 = %d\n",JMAX-4); /* remove boundary region */
fprintf(fpfld,"dim2 = %d\n",IMAX-4); /* remove boundary region */
fprintf(fpfld,"nspace = 2\n");
fprintf(fpfld,"veclen = 5\n");
fprintf(fpfld,"data = double\n");
fprintf(fpfld,"field = irregular\n");
fprintf(fpfld,"label = pressure density u v phi\n");
fprintf(fpfld,"\n");
fprintf(fpfld,"coord 1 file=%s filetype=ascii skip=2 offset=0 stride=7\n",avsdatl);
fprintf(fpfld,"coord 2 file=%s filetype=ascii skip=2 offset=1 stride=7\n",avsdatl);
fprintf(fpfld,"variable 1 file=%s filetype=ascii skip=2 offset=2 stride=7\n",avsdatl);
fprintf(fpfld,"variable 2 file=%s filetype=ascii skip=2 offset=3 stride=7\n",avsdatl);
fprintf(fpfld,"variable 3 file=%s filetype=ascii skip=2 offset=4 stride=7\n",avsdatl);
fprintf(fpfld,"variable 4 file=%s filetype=ascii skip=2 offset=5 stride=7\n",avsdatl);
fprintf(fpfld,"variable 5 file=%s filetype=ascii skip=2 offset=6 stride=7\n",avsdatl);
fclose(fpfld);
velo_calc(u,v,us,vs,uv,vu);
bc_velo_sc(us,1);
bc_velo_sc(vs,0);
/* no-slip case
bc_velo_sc_ns(us);
bc_velo_sc_ns(vs);
*/
if ((fpdat=fopen (avsdat,"w+")) == NULL){
printf("%s was not opened!\n",avsdat);
printf("’Maybe you should make ’Data’ directory....\n");
exit(1);
}
fprintf(fpdat,"t=%16.8e\n",time);
fprintf(fpdat,"%16s %16s %16s %16s %16s %16s %16s\n","x","y","pressure","density","u","w","phi");
for(i=2; i<IMAX-2; i++){/* remove boundary region */
for(j=2; j<JMAX-2; j++){/* remove boundary region */
fprintf(fpdat,"%16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",(i-2)*DX,(j-2)*DY,p[i][j],rho[i][j],us[i][j],vs[i][j],iphi(phi[i][j]));
}
}
fclose(fpdat);
return(0);
}
int output_monitor(double p[IMAX][JMAX], double rho[IMAX][JMAX], double u[IMAX][JMAX], double v[IMAX][JMAX], double phi[IMAX][JMAX],
int count, double time, int ml){
FILE *fp;
int i, j, d, ymax;
double mx, depth,d0,va=0.044;
char fpfld[25]= "./Data/x????.dat";
if(ml>9999 || ml>=IMAX){
printf("(output_monitor) Monitoring line definition ERROR. Limit is 9999 or %d(=IMAX), but your definition is %d.\n",IMAX-1,ml);
exit(-1);
}
fpfld[8] = '0' + ml/1000;
fpfld[9] = '0' + (ml%1000)/100;
fpfld[10] = '0' + ((ml%1000)%100)/10;
fpfld[11] = '0' + ((ml%1000)%100)%10;
if ((fp=fopen(fpfld,"a")) == NULL){
printf("%s was not opened!\n",fpfld);
printf("’Maybe you should make ’Data’ directory....\n");
exit(1);
}
mx=(ml-1.5)*DX;// slx.add.explain: mx=((ml+1)-2-0.5)*DX=(ml-1.5)*DX; 2代表邊界區域(要除去),0.5代表格子中心
d=0;j=2;
while(iphi(phi[ml][j])>=0.5){
d++; j++;
}
depth=d*DY;
ymax=INLET*2.5;
d0=INLET*DY;
if(T<=0.){
fprintf(fp,"%16s","d0=");
fprintf(fp,"%16.8e\n",d0);
fprintf(fp,"%16s","g=");
fprintf(fp,"%16.8e\n",G);
fprintf(fp,"%16s","V*=");
fprintf(fp,"%16.8e\n",va);
fprintf(fp,"%16s","x=");
fprintf(fp,"%16.8e",mx);
fprintf(fp,"%32s"," ");
for(j=2; j<ymax; j++) fprintf(fp, "%16s%16s%16s%16s","u-at(y=)","u/V*-at(z/d0=)","v-at(y=)","v/V*-at(z/d0=)");
fprintf(fp,"\n");
fprintf(fp,"%16s"," Time");
fprintf(fp,"%16s"," t.sqrt(g/do)");
fprintf(fp,"%16s"," Depth");
fprintf(fp,"%16s"," d/d0");
for(j=2; j<ymax; j++) fprintf(fp, "%16.8f%16.8f%16.8f%16.8f",(j-1.5)*DY,(j-1.5)*DY/d0,(j-2)*DY,(j-2)*DY/d0);
fprintf(fp,"\n");
}
fprintf(fp,"%16.8e",T);
fprintf(fp,"%16.8e",T*sqrt(G/d0));
fprintf(fp,"%16.8e",depth);
fprintf(fp,"%16.8e",depth/d0);
for(j=2; j<ymax; j++) fprintf(fp, "%16.8e%16.8e%16.8e%16.8e",-u[ml][j],-u[ml][j]/va,v[ml][j],v[ml][j]/va);
fprintf(fp,"\n");
fclose(fp);
return(0);
}
int diffusion_velo_MIXINGLENGTH_NOTFIXED(double u[IMAX][JMAX], double v[IMAX][JMAX],
double rho[IMAX][JMAX], double phi[IMAX][JMAX],
double qu[IMAX][JMAX], double qv[IMAX][JMAX]){
/* WARNING NO_FIXED version : Mixing Length */
int i,j;
double u_s[IMAX][JMAX], v_s[IMAX][JMAX], u_v[IMAX][JMAX], v_u[IMAX][JMAX];
double lmbx, lmx[IMAX][JMAX];
double lmby, lmy[IMAX][JMAX];
double dudx_u, dudy_u, dvdx_u, dvdy_u;
double dudx_v, dudy_v, dvdx_v, dvdy_v;
double xDxx, xDxy, xDyx, xDyy;
double yDxx, yDxy, yDyx, yDyy;
double etax, etay;
double rsux[IMAX][JMAX], rsuy[IMAX][JMAX], rsvx[IMAX][JMAX], rsvy[IMAX][JMAX];
double ud, vd, mu;
printf("WARNING! WARNING! WARNING! NO-FIXED mixing length version for diffusion_velo()\n");
velo_calc(u,v,u_s,v_s,u_v,v_u);
bc_velo(v_u, u_v);
for(i=0; i<=IMAX-1; i++){
for(j=0; j<=JMAX-1; j++){
lmx[i][j]=VKC*(j-2+0.5)*DY;
lmy[i][j]=VKC*(j-2)*DY;
}
}
/* lmx,y 0 - i - IMAX-1 , 1 - j - JMAX-2 */
for(i=1; i<=IMAX-2; i++){
for(j=1; j<=JMAX-2; j++){
/* dudy and dvdx at cell-center */
dudx_u=(u[i+1][j]-u[i-1][j])/(2.*DX);
dudy_u=(u[i][j+1]-u[i][j-1])/(2.*DY);
dvdx_u=(v_u[i+1][j]-v_u[i-1][j])/(2.*DX);
dvdy_u=(v_u[i][j+1]-v_u[i][j-1])/(2.*DY);
dudx_v=(u_v[i+1][j]-u_v[i-1][j])/(2.*DX);
dudy_v=(u_v[i][j+1]-u_v[i][j-1])/(2.*DY);
dvdx_v=(v[i+1][j]-v[i-1][j])/(2.*DX);
dvdy_v=(v[i][j+1]-v[i][j-1])/(2.*DY);
/* Dij = 0.5*(@ui/@xj + @uj/@xi) */
/* for ud */
xDxx=0.5*(dudx_u+dudx_u);
xDxy=0.5*(dudy_u+dvdx_u);
xDyx=0.5*(dvdx_u+dudy_u);
xDyy=0.5*(dvdy_u+dvdy_u);
/* for vd */
yDxx=0.5*(dudx_v+dudx_v);
yDxy=0.5*(dudy_v+dvdx_v);
yDyx=0.5*(dvdx_v+dudy_v);
yDyy=0.5*(dvdy_v+dvdy_v);
etax = 0.5*(rho[i][j]+rho[i-1][j])*pow(lmx[i][j],2)
*sqrt( 2.*(xDxx*xDxx + xDxy*xDxy + xDyx*xDyx + xDyy*xDyy));
etay = 0.5*(rho[i][j]+rho[i][j-1])*pow(lmy[i][j],2)
*sqrt( 2.*(yDxx*yDxx + yDxy*yDxy + yDyx*yDyx + yDyy*yDyy));
rsux[i][j]=etax*dudx_u;
rsuy[i][j]=etax*dudy_u;
rsvx[i][j]=etay*dvdx_v;
rsvy[i][j]=etay*dvdy_v;
}
}
for(i=2; i<=IMAX-3; i++){
for(j=2; j<=JMAX-3; j++){
ud=(u[i-1][j]-2.*u[i][j]+u[i+1][j])/pow(DX,2)+(u[i][j-1]-2.*u[i][j]+u[i][j+1])/pow(DY,2)
+(rsux[i+1][j]-rsux[i-1][j])/(2.*DX) + (rsuy[i][j+1]-rsuy[i][j-1])/(2.*DY);
vd=(v[i-1][j]-2.*v[i][j]+v[i+1][j])/pow(DX,2)+(v[i][j-1]-2.*v[i][j]+v[i][j+1])/pow(DY,2)
+(rsvx[i+1][j]-rsvx[i-1][j])/(2.*DX) + (rsvy[i][j+1]-rsvy[i][j-1])/(2.*DY);
if(iphi(phi[i][j])<0. || iphi(phi[i][j])>1.){
printf("(diff_velo_ML) INVALID iphi value (=%16.8e) at (%d,%d)\n",iphi(phi[i][j]),i,j);
exit(-1);
}
mu=(1.-iphi(phi[i][j]))*MUair+iphi(phi[i][j])*MUwater;
if(mu<MUair || mu>MUwater){
printf("INVALID MU value(=%16.8e) at (%d,%d)\n",mu,i,j);
exit(-1);
}
qu[i][j]=mu/(0.5*(rho[i][j]+rho[i-1][j ]))*ud;
qv[i][j]=mu/(0.5*(rho[i][j]+rho[i ][j-1]))*vd;
}
}
/* for i: 0,1, IMAX-2, IMAX-1 <- not yet calculated.... */
/* for j: 0,1, JMAX-2, JMAX-1 <- not yet calculated.... */
return(0);
}
int diffusion_velo_MixingLength20071116(double u[IMAX][JMAX], double v[IMAX][JMAX],
double rho[IMAX][JMAX], double phi[IMAX][JMAX],
double qu[IMAX][JMAX], double qv[IMAX][JMAX]){
/* Mixing Length (20071116-:SF)*/
int i,j;
double u_s[IMAX][JMAX], v_s[IMAX][JMAX], u_v[IMAX][JMAX], v_u[IMAX][JMAX];
double dudx, dvdy, dudy, dvdx, ddudyy;
double J, lm;
double nut[IMAX][JMAX], nu[IMAX][JMAX];
double d11[IMAX][JMAX], d12[IMAX][JMAX], d21[IMAX][JMAX], d22[IMAX][JMAX];
double nuyp, nutyp, d12yp, nuym, nutym, d12ym;
double nuxp, nutxp, d21xp, nuxm, nutxm, d21xm;
/* redefine velocity points from cell-boundary to cell-venter */
velo_calc(u,v,u_s,v_s,u_v,v_u);
/* boundary condition for redefined velocity points */
bc_velo(v_u, u_v);
for(i=0; i<IMAX-1; i++){
for(j=1; j<JMAX-1; j++){
if(iphi(phi[i][j])>1 ||iphi(phi[i][j])<0 ){
printf("INVALID define of phi(=%8.4f, iphi=%8.4f) at (%d,%d) in diffusion_velo()\n",
phi[i][j],iphi(phi[i][j]),i,j);
exit(-1);
}
/* Kinematic Eddy Viscosity : NUT= lm^2 * J */
/* dudy and dvdx at cell-center */
dudx=(u[i+1][j]-u[i][j])/DX;
dvdy=(v[i][j+1]-v[i][j])/DY;
dudy=(u_v[i][j+1]-u_v[i][j])/DY;
dvdx=(v_u[i+1][j]-v_u[i][j])/DX;
/* for J */
J=sqrt( 2.*pow(dudx,2) +2.*pow(dvdy,2) +pow(dvdx+dudy,2));
/*............ for lm : mixing length .............................*/
// ddudyy=fabs((u[i][j+1]-2.*u[i][j]+u[i][j-1])/pow(DY,2));
// if(ddudyy<=0.)ddudyy=1e-10;
// lm=VKC*fabs(dudy/ddudyy);
// if j is 2, ycdn(cell-center) is DY/2..
lm=VKC*(j-2+0.5)*DY;
/*............ for lm : mixing length .............................*/
/*............ for NUT .........................*/
nut[i][j]=pow(lm,2)*J;
// nut[i][j]=iphi(phi[i][j])*pow(lm,2)*J;
// nut[i][j]=0.0;
/*............ for NUT .........................*/
/* for Dij */
d11[i][j]=0.5*(dudx+dudx);
d12[i][j]=0.5*(dvdx+dudy);
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -