2次元キャビティ流れNO.2(C)プログラム
2003/09/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元キャビティ流れを速度-圧力法により差分法で計算する。
2. プログラム
//=================================================================C
// 2-Dimensional Cavity Flow by velocity and pressure C
// with FDM and SOR C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/09/10 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
#include <stdio.h>
#include <math.h>
#include <time.h>
// Global Define
#define NDX 301
#define NDY 301
double P[NDX][NDY], B[NDX][NDY] ;
double U[NDX][NDY], V[NDX][NDY] ;
FILE *FT1 ;
// Function Definition
int PSOR(int NX, int NY, double DT, double H, double ALP, double EPS) ;
void OUTUVP(int NX, int NY) ;
// Main Program
int main()
{ int i, j, k, NX, NY, NT, MT, ID, T1, T2 ;
int ITER, NTER, NI1, NI2, IER ;
double Re, DT, T, H, DH, H2, D2, D5, EPS, ALP, CPU, OVER ;
char FLOW2[13]={'F','L','O','W','2','-','?','?','.','t','x','t',} ;
char NTABL[10]={'0','1','2','3','4','5','6','7','8','9'} ;
// Initial Data
ID = 0 ;
IER = 0 ;
printf("Type In NX,Re \n") ;
scanf("%d %lf",&NX,&Re) ;
printf("DT,NT(Total NO.),MT(Output Interval) \n") ;
scanf("%lf %d %d",&DT,&NT,&MT) ;
if(NX >= NDX) { NX = NDX - 1 ; }
NY = NX ;
printf("# Cavity Flow Numerical Analysis \n") ;
printf("NX=%d NY=%d Re=%f DT=%f NT,MT=%d %d \n",NX,NY,Re,DT,NT,MT) ;
// Initial Constant
T = 0.0 ;
DH = NX ;
D5 = DH/2.0 ;
H = 1.0/DH ;
H2 = H*H ;
D2 = DH*DH ;
// Set Initial Condition
for (i=0; i<=NX; i++) {
for (j=0; j<=NY; j++) {
P[i][j] = 0.0 ;
U[i][j] = 0.0 ;
V[i][j] = 0.0 ; }
}
// U=1.0 on L1
for (i=0; i<=NX; i++) {
U[i][NY] = 1.0 ; }
// SOR Parameter (ALP)
ALP = 1.0 + log(NX*1.0)/log(NDX*1.2) ;
EPS = 1.0e-4 ;
// Main Loop
NTER = 0 ;
T1 = clock() ;
for (k=1; k<=NT; k++) {
T = T + DT ;
// Compute U,V
OVER = 0.0 ;
for (i=1; i<NX; i++) {
for (j=1; j<NY; j++) {
U[i][j] = U[i][j] + DT*((U[i-1][j]*U[i-1][j]-U[i+1][j]*U[i+1][j]
+ U[i][j-1]*V[i][j-1] - U[i][j+1]*V[i][j+1]
+ P[i-1][j] - P[i+1][j])*D5 - (4.0*U[i][j] - U[i][j-1]
- U[i-1][j] - U[i+1][j] - U[i][j+1])*D2/Re ) ;
V[i][j] = V[i][j] + DT*((V[i][j-1]*V[i][j-1]-V[i][j+1]*V[i][j+1]
+ U[i-1][j]*V[i-1][j] - U[i+1][j]*V[i+1][j]
+ P[i][j-1] - P[i][j+1])*D5 - (4.0*V[i][j] - V[i][j-1]
- V[i-1][j] - V[i+1][j] - V[i][j+1])*D2/Re ) ;
OVER = OVER + fabs(U[i][j]) + fabs(V[i][j]) ; }
}
// Check over-flow computation
if(OVER >= 1.0e10) {
printf("** Stop for over-flow computation ** \n") ;
printf(" You have to give smaller DT \n") ;
IER = 1 ;
goto M100 ; }
// Compute P
ITER = PSOR(NX, NY, DT, H, ALP, EPS) ;
NTER = NTER + ITER ;
// Output U,V,Phi
if( k%MT == 0) {
NI1 = ID/10 ;
NI2 = ID - NI1*10 ;
FLOW2[6] = NTABL[NI1] ;
FLOW2[7] = NTABL[NI2] ;
FT1 = fopen(FLOW2,"w") ;
fprintf(FT1,"# Cavity Flow Numerical Analysis \n") ;
fprintf(FT1,"NX=%d NY=%d Re=%f DT=%f NT,MT=%d %d \n",NX,NY,Re,DT,NT,MT) ;
fprintf(FT1,"#Time=%7.2f, Step=%5d SOR-loop=%6d \n",T,k,NTER) ;
printf("#Time=%7.2f, Step=%5d SOR-loop=%6d \n",T,k,NTER) ;
NTER = 0 ;
OUTUVP(NX, NY) ;
fclose(FT1) ;
if(ID < 99) { ID = ID + 1 ; }
}
}
T2 = clock() ;
CPU = (double)(T2 - T1)/CLOCKS_PER_SEC ;
printf(" NX=%4d, NT=%5d, Time(s)=%9.2f \n",NX,NT,CPU) ;
M100: return (IER) ;
}
//=================================================================C
int PSOR(int NX, int NY, double DT, double H, double ALP, double EPS)
//=================================================================C
// Solve Ax=b by SOR with 2 dimensional FDM C
// Given Omega ( Acceleration factor ) C
//-----------------------------------------------------------------C
// PHI,OMG ; Global Array C
// NX I*4, In, Grid Numbers on X-axis C
// NY I*4, In, Grid Numbers on Y-axis C
// DT R*8, In, Delta T C
// H R*8, In, H=1.0/NX C
// ALP R*8, In, SOR Acceleration factor C
// EPS R*8, In, if ||r||/||b|| <= EPS --> return C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/09/10 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
{ int i, j, k, IP, IM, JP, JM, ITER ;
double BN, RN, ERR, R ;
double UP0, UM0, U0P, U0M, VP0, VM0, V0P, V0M ;
// Set B and BN (B norm)
BN = 0.0 ;
for (i=0; i<=NX; i++) {
for (j=0; j<=NY; j++) {
// Boundary on L1
if(j != NY) { U0P=U[i][j+1] ; V0P=V[i][j+1] ; }
else { U0P=1.0 ; V0P=0.0 ; }
// Boundary on L2
if(i != 0) { UM0=U[i-1][j] ; VM0=V[i-1][j] ; }
else { UM0=U[1][j] ; VM0=-V[1][j] ; }
// Boundary on L3
if(j != 0) { U0M=U[i][j-1] ; V0M=V[i][j-1] ; }
else { U0M=-U[i][1] ; V0M=V[i][1] ; }
// Boundary on L4
if(i != NX) { UP0=U[i+1][j] ; VP0=V[i+1][j] ; }
else { UP0=U[NX-1][j] ; VP0=-V[NX-1][j] ; }
// Computation
B[i][j] = ((UP0-UM0)*(UP0-UM0) + (V0P-V0M)*(V0P-V0M))/4.0
+ (VP0 - VM0)*(U0P - U0M)/2.0
- (UP0 - UM0 + V0P - V0M)*H/(2.0*DT) ;
BN = BN + B[i][j]*B[i][j] ; }
}
// Main Loop
for (k=1; k<=NX*NY; k++) {
RN = 0.0 ;
for (i=0; i<=NX; i++) {
IP = i + 1 ;
IM = i - 1 ;
// Boundary on L1,L3
if(i == NX) { IP = NX - 1 ; }
if(i == 0) { IM = 1 ; }
for (j=0; j<=NY; j++) {
JP = j + 1 ;
JM = j - 1 ;
// Boundary on L2,L4
if(j == NY) { JP = NY - 1 ; }
if(j == 0) { JM = 1 ; }
// Computation
R = (B[i][j]+P[i][JM]+P[IM][j]+P[IP][j]+P[i][JP])/4.0-P[i][j] ;
P[i][j] = P[i][j] + ALP*R ;
// Reset on P1
if(i==0 & j==0) { P[i][j]=0.0 ; R=0.0 ; }
RN = RN + R*R*16.0 ;
} }
// if(ERR <= EPS) return
ERR = sqrt(RN/BN) ;
if(ERR <= EPS) goto M100 ;
}
M100: ITER = k ;
return (ITER) ;
}
//=================================================================C
void OUTUVP(int NX, int NY)
//=================================================================C
// Output U,V,P C
//-----------------------------------------------------------------C
// P,U,V ; Global Array C
// NX I*4, In, Grid Numbers on X-axis C
// NY I*4, In, Grid Numbers on Y-axis C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/09/10 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
{ int i, j ;
double X, Y, DX, DY ;
// Output
DX = 1.0/NX ;
DY = 1.0/NY ;
for (j=0; j<=NY; j++) {
Y = j*DY ;
for (i=0; i<=NX; i++) {
X = i*DX ;
fprintf(FT1,"%9.5f %9.5f %9.5f %9.5f %9.5f \n",
X,Y,P[i][j],U[i][j],V[i][j]) ; }
}
}