?? 2.txt

?? 基于CIP的海嘯波浪數(shù)值模擬
?? TXT
?? 第 1 頁 / 共 5 頁
字號:
}
/* RIGHT side */
/*... inlet&open boundary ...*/
for(j=2; j<=JMAX-3; j++){
pn[IMAX-1][j]=pn[IMAX-3][j];
pn[IMAX-2][j]=pn[IMAX-3][j];
}
/* BOTTM */
for(i=2; i<=IMAX-3; i++){
/* bottom */
pn[i][0]=pn[i][3];
pn[i][1]=pn[i][2];
}
/* TOP */
for(i=0; i<=IMAX-1; i++){
pn[i][JMAX-1]=Patm;
pn[i][JMAX-2]=Patm;
}
/* corner */
/* left & bottom */
pn[0][0]=pn[3][3];
pn[0][1]=pn[3][2];
pn[1][0]=pn[2][3];
pn[1][1]=pn[2][2];
/* right & bottom */
pn[IMAX-1][0]=pn[IMAX-4][3];
pn[IMAX-1][1]=pn[IMAX-4][2];
pn[IMAX-2][0]=pn[IMAX-3][3];
pn[IMAX-2][1]=pn[IMAX-3][2];
/********** B.C. for Pressure *****************************/
return(0);
}


int psor(double rho[IMAX][JMAX], double u[IMAX][JMAX], double v[IMAX][JMAX], double phi[IMAX][JMAX], double p[IMAX][JMAX], double pn[IMAX][JMAX], int rtt){
	int i, j;
	int ite=0;
	double oldp[IMAX][JMAX], op;
	double res=10., rhs=0., bnm=0.,npt;
	double fdx, fdy, fdx2, fdy2, DT2, fDT, fDT2;
	double b0, b1, b2, b3, b4, b;
	double div;
	double rhol, rhor, rhou, rhod;
	double dx2, dy2, cs2, rhoav;
	for(i=0; i<IMAX; i++){
	for(j=0; j<JMAX; j++){
	oldp[i][j]=p[i][j];
	}
	}
	dx2=DX*DX;
	dy2=DY*DY;
	fdx=1./DX;
	if(DX<=0.){
	printf("invalid fdx(=%16.8e). aborted! in psor()\n",fdx);
	exit(-1);
	}
	fdy=1./DY;
	if(DY<=0.){
	printf("invalid fdy(=%16.8e). aborted! in psor()\n",fdy);
	exit(-1);
	}
	fdx2=1./dx2;
	if(pow(DX,2)<=0.){
	printf("invalid fdx2(=%16.8e). aborted! in psor()\n",fdx2);
	exit(-1);
	}
	fdy2=1./dy2;
	if(pow(DY,2)<=0.){
	printf("invalid fdy2(=%16.8e). aborted! in psor()\n",fdy2);
	exit(-1);
	}
	DT2=pow(DT,2);
	if(DT2<=0.){
	printf("invalid DT2(=%16.8e). aborted! in psor()\n",DT2);
	exit(-1);
	}
	fDT=1./DT;
	if(DT<=0.){
	printf("invalid fDT(=%16.8e). aborted! in psor()\n",fDT);
	exit(-1);
	}
	fDT2=1./DT2;
	if(fDT2<=0.){
	printf("invalid DT2(=%16.8e). aborted! in psor()\n",DT2);
	exit(-1);
	}
	while ( (res>EPS) && (ite<MITE)){
	for(i=0; i<IMAX; i++){
	for(j=0; j<JMAX; j++){
	p[i][j]=ALPHA*pn[i][j]+(1.-ALPHA)*p[i][j];
	/*
	sample
	pp p
	p pn
	yp oldp
	*/
	}
	}
	ite++;
	res=0.;
	for(j=2; j<=JMAX-3; j++){
	for(i=2; i<=IMAX-3; i++){
	cs2=sound_velocity(i,j,rho,oldp,phi);
	rhol=(rho[i][j]+rho[i-1][j ])*0.5;
	rhor=(rho[i][j]+rho[i+1][j ])*0.5;
	rhod=(rho[i][j]+rho[i ][j-1])*0.5;
	rhou=(rho[i][j]+rho[i ][j+1])*0.5;
	rhoav=0.25*(rhol+rhor+rhou+rhod);
	if(rhol<=0.){
	printf("Negative Density at [%d][%d] rhol=%30.22e\n",i,j,rhol);
	exit(-1);
	}
	if(rhor<=0.){
	printf("Negative Density at [%d][%d] rhor=%30.22e\n",i,j,rhor);
	exit(-1);
	}
	if(rhod<=0.){
	printf("Negative Density at [%d][%d] rhod=%30.22e\n",i,j,rhod);
	exit(-1);
	}
	if(rhou<=0.){
	printf("Negative Density at [%d][%d] rhou=%30.22e rho=%30.22e rho(j+1)=%30.22e\n",i,j,rhou,rho[i][j],rho[i][j+1]);
	exit(-1);
	}
	/*
	sample
	pp p
	p pn
	yp oldp
	*/
	/* original sample SOR ( not included sound speed ) */
	// div=1.0+GAMMA*oldp[i][j]*DT*DT*((1.0/rhor+1.0/rhol)/dx2+(1.0/rhou+1.0/rhod)/dy2);
	// pn[i][j]=(oldp[i][j]-GAMMA*oldp[i][j]*DT*((u[i+1][j]-u[i][j])/DX - (p[i+1][j]/rhor+p[i-1][j]/rhol)*DT/dx2
	// +(v[i][j+1]-v[i][j])/DY - (p[i][j+1]/rhou + p[i][j-1]/rhod)*DT/dy2))/div;
	/* (improved) sample SOR, including sound speed */
	div=1.0+cs2*rho[i][j]*DT*DT*((1.0/rhor+1.0/rhol)/dx2+(1.0/rhou+1.0/rhod)/dy2);
	pn[i][j]=(oldp[i][j]-cs2*rho[i][j]*DT*((u[i+1][j]-u[i][j])/DX - (p[i+1][j]/rhor+p[i-1][j]/rhol)*DT/dx2
	+(v[i][j+1]-v[i][j])/DY - (p[i][j+1]/rhou + p[i][j-1]/rhod)*DT/dy2))/div;
	if(pn[i][j]<=0.){
	printf("Negative Pressure at [%d][%d] p=%30.22e\n",i,j,pn[i][j]);
	return(-1);
	}
	res+=pow(p[i][j]-pn[i][j],2);
	}
	}
	if(ite%200==0)printf("--(SOR)%4d st.ite., res.=%16.8e (%6.1f t.l.t. EPS)...(SOR)\n",ite,res,(res/EPS));
	/* B.C. */
	bc_p(pn);
	/* ... end of pressure iteration loop ... */
	}
	if(ite==MITE)printf("--not converged within %4d steps iterations, residual=%16.8e....(SOR)\n",MITE,res);
	else printf("--converged in %4d steps iterations. res=%16.8e...(SOR)\n",ite,res);
	for(i=0; i<IMAX; i++){
	for(j=0; j<JMAX; j++){
	p[i][j]=oldp[i][j];
	}
	}
	return(0);
	}


int isoentro_rho(double rho[IMAX][JMAX], double rhon[IMAX][JMAX],
double p[IMAX][JMAX], double pn[IMAX][JMAX],
double phi[IMAX][JMAX],
double un[IMAX][JMAX], double vn[IMAX][JMAX]){
int i,j;
double cs2;
double pnl, pnr, pnd, pnu, pna;
double pl, pr, pd, pu, pa;
// --------------------------------------------------------------------------------
printf("(isoentro_rho) this function is NOT CONSIDERED to use....\n");
exit(-1);
// --------------------------------------------------------------------------------
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
// ----- shimizu2001 ---------------------------------------------------------------------
// cs2=sound_velocity(i, j, rho, p, phi);
// rhon[i][j]=rho[i][j]+1./cs2*(pn[i][j]-p[i][j]);
// ---------------------------------------------------------------------------------------
// ----- mizutani2002 --------------------------------------------------------------------
rhon[i][j]=rho[i][j]-rho[i][j]*((un[i+1][j]-un[i][j])/DX+(vn[i][j+1]-vn[i][j])/DY)*DT;
if(rhon[i][j]<=RHOair)rhon[i][j]=RHOair;
if(rhon[i][j]>=RHOwater)rhon[i][j]=RHOwater;
// ---------------------------------------------------------------------------------------
}
}
return(0);
}
int rho_from_phi(double rhon[IMAX][JMAX], double phin[IMAX][JMAX]){
int i,j;
double tmpphi[IMAX][JMAX];
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
rhon[i][j]=iphi(phin[i][j])*RHOwater + (1.-iphi(phin[i][j]))*RHOair;
if(rhon[i][j]<=0.){
printf("(rho_from_phi) Negative Density(=%8.4f) at (%d,%d), phin=%8.4f\n",
rhon[i][j], i, j, phin[i][j]);
exit(-1);
}
}
}
return(0);
}
int rho_from_phi_filtered(double rhon[IMAX][JMAX], double phin[IMAX][JMAX]){
int i,j, mgc=0, agc=0, wgc=0;
double tmpphi[IMAX][JMAX];
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
tmpphi[i][j]=0.2*(iphi(phin[i][j])+iphi(phin[i][j])+iphi(phin[i][j+1])+iphi(phin[i+1][j])+iphi(phin[i+1][j+1]));
}
}
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
if(tmpphi[i][j]<=0.2){
phin[i][j]=dphi(0.);
rhon[i][j]=RHOair;
agc++;
}
else if(tmpphi[i][j]>=0.8){
phin[i][j]=dphi(1.);
rhon[i][j]=RHOwater;
wgc++;
}
else{
rhon[i][j]=iphi(phin[i][j])*RHOwater + (1.-iphi(phin[i][j]))*RHOair;
phin[i][j]=dphi(0.5);
mgc++;
}
if(rhon[i][j]<=0.){
printf("(rho_from_phi) Negative Density(=%8.4f) at (%d,%d), phin=%8.4f\n",
rhon[i][j], i, j, phin[i][j]);
exit(-1);
}
}
}
printf("(rho_from_phi) Air=%d Water=%d Mid=%d\n",agc,wgc,mgc);
return(0);
}
int correct_phi(double phin[IMAX][JMAX]){
/* original */
int i,j,nom=0;
double mixture[IMAX][JMAX], vlm=0., tvlm=0., difvlm, bi;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
if(iphi(phin[i][j])<0.){/* air cell */
phin[i][j]=dphi(0.);
}
else if(iphi(phin[i][j])>1.){/* water cell */
phin[i][j]=dphi(1.);
}
if(i>=2 && i<IMAX-1 && j>=2 && j<JMAX-1) vlm+=iphi(phin[i][j]); /* total volume */
}
}
return(0);
}
int correct_phi_isoentro(double phin[IMAX][JMAX]){
int i,j;
int na=0, nf=0, nm=0;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
if(iphi(phin[i][j])<=0.){
phin[i][j]=dphi(0.);
na++;
}
else if(iphi(phin[i][j])>=1.){
phin[i][j]=dphi(1.);
nf++;
}
else{
phin[i][j]=dphi(0.5);
nm++;
}
}
}
printf("Air:%3d(%8.4f%%) Water:%3d(%8.4f%%) Mix:%3d(%8.4f%%)\n",
na,(double)na/(double)(IMAX*JMAX)*100., nf,(double)nf/(double)(IMAX*JMAX)*100, nm,(double)nm/(double)(IMAX*JMAX)*100.);
return(0);
}
int correct_phi_adv(double phin[IMAX][JMAX]){
int i,j;
int na=0, nf=0, nm=0;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
if(iphi(phin[i][j])<0.5){
phin[i][j]=dphi(0.);
na++;
}
else if(iphi(phin[i][j])>0.5){
phin[i][j]=dphi(1.);
nf++;
}
else{
phin[i][j]=dphi(0.5);
nm++;
}
}
}
printf("Air:%3d(%8.4f%%) Water:%3d(%8.4f%%) Mix:%3d(%8.4f%%)\n",
na,(double)na/(double)(IMAX*JMAX)*100., nf,(double)nf/(double)(IMAX*JMAX)*100, nm,(double)nm/(double)(IMAX*JMAX)*100.);
return(0);
}
int isoentro_phi(double phi[IMAX][JMAX], double phin[IMAX][JMAX], double rho[IMAX][JMAX],
double p[IMAX][JMAX], double pn[IMAX][JMAX],
double un[IMAX][JMAX], double vn[IMAX][JMAX]){
int i,j;
double cs2;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
// ----- shimizu2001 ---------------------------------------------------------------------
// cs2=sound_velocity(i, j, rho, p, phi);
// phin[i][j]=phi[i][j]+phi[i][j]/(rho[i][j]*cs2)*(pn[i][j]-p[i][j]);
// ----- shimizu2001 ---------------------------------------------------------------------
// ----- mizutani2002 ---------------------------------------------------------------------
phin[i][j]=phi[i][j]-DT* phi[i][j]*((un[i+1][j]-un[i][j])/DX+ (vn[i][j+1]-vn[i][j])/DY);
// ----- mizutani2002 ---------------------------------------------------------------------
}
}
correct_phi(phin);
return(0);
}
int isoentro_update(double u[IMAX][JMAX], double un[IMAX][JMAX],
double v[IMAX][JMAX], double vn[IMAX][JMAX],
double rho[IMAX][JMAX], double rhon[IMAX][JMAX],
double phi[IMAX][JMAX], double phin[IMAX][JMAX]){
int i,j;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
u[i][j]= un[i][j];
v[i][j]= vn[i][j];
rho[i][j]=rhon[i][j];
phi[i][j]=phin[i][j];
}
}
return(0);
}
int diffusion_update(double u[IMAX][JMAX], double un[IMAX][JMAX],
double v[IMAX][JMAX], double vn[IMAX][JMAX]){
int i,j;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
u[i][j]=un[i][j];
v[i][j]=vn[i][j];
}
}
return(0);
}
int sc_update(double scn[IMAX][JMAX], double sc[IMAX][JMAX]){
int i,j;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
sc[i][j]=scn[i][j];
}
}
return(0);
}
//---------- Boundary Condition ----------
int bc_scalar(double sc[IMAX][JMAX], double dvalue){
int i,j;
/***** LEFT SIDE *****/
/*... outlet ...*/
for(j=2; j<=OUTLET+1; j++){
/* left */
sc[0 ][j]=sc[2 ][j];
sc[1 ][j]=sc[2 ][j];
}
/*... wall ...*/
for(j=OUTLET+2; j<=JMAX-1; j++){
/* left */
sc[0 ][j]=sc[3 ][j];
sc[1 ][j]=sc[2 ][j];
}
/***** RIGHT SIDE *****/
/*... inlet ...*/
for(j=2; j<=INLET+1; j++){
sc[IMAX-1][j]=dvalue;
sc[IMAX-2][j]=dvalue;
}
/*... openboundary on inlet ...*/
for(j=INLET+2; j<=JMAX-1; 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];
}
/***** TOP (open boundary) *****/
for(i=2; i<=IMAX-3; i++){
/* top.. wall b.c. */
sc[i][JMAX-1]=sc[i][JMAX-3];
sc[i][JMAX-2]=sc[i][JMAX-3];
}
/* corner */
/* left & bottom
| | | | |
+---+---+---+---+--
3 | | | B | A |
+---+---+---+---+--
2 | | | D | C |
+---+---+---+---+--
1 | C | D | | |
+---+---+---+---+--
0 | A | B | | |
+---+---+---+---+--
0 1 2 3
*/
sc[0][0]=sc[3][3]; /* A */
sc[0][1]=sc[3][2]; /* C */
sc[1][0]=sc[2][3]; /* B */
sc[1][1]=sc[2][2]; /* D */
/* right & bottom
| | | |
-+---+---+---+---
3 | B | A | |
-+---+---+---+---
2 | D | C | |
-+---+---+---+---
1 | | | C | D
-+---+---+---+---
0 | | | A | B
-+---+---+---+---
-4 -3 -2 -1
*/
sc[IMAX-1][0]=sc[IMAX-4][3]; /* B */
sc[IMAX-1][1]=sc[IMAX-4][2]; /* D */
sc[IMAX-2][0]=sc[IMAX-3][3]; /* A */
sc[IMAX-2][1]=sc[IMAX-3][2]; /* C */
/* left & top (wall condition)
-1 | A | B | |
+---+---+---+---
?? 文件大小 58 K
?? 上傳用戶 csytml
?? 所屬分類其他
??? 相關(guān)標(biāo)簽

#CIP #海 #數(shù)值模擬
?? 快捷鍵說明

復(fù)制代碼 Ctrl + C
搜索代碼 Ctrl + F
全屏模式 F11
切換主題 Ctrl + Shift + D
顯示快捷鍵 ?
增大字號 Ctrl + =
減小字號 Ctrl + -
亚洲欧美第一页_禁久久精品乱码_粉嫩av一区二区三区免费野_久草精品视频

?? 2.txt

?? 快捷鍵說明
