2次元角柱周りの流れNO.3(C)プログラム
2003/09/17 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元角柱周りの流れを速度-圧力法により差分法で計算する。
2. プログラム
//=================================================================C
// 2-Dimensional Flow around Square a Pillar by velocity C
// and pressure with FDM and SOR C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/09/14 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
#include <stdio.h>
#include <math.h>
#include <time.h>
// Global Define
#define NDX 302
#define NDY 302
#define P(i,j) P[j+1][i+1]
#define B(i,j) B[j+1][i+1]
#define U(i,j) U[j+1][i+1]
#define V(i,j) V[j+1][i+1]
#define PHI(i,j) PHI[j+1][i+1]
#define A(i,j,k) A[k-1][j+1][i+1]
#define OMG(i,j) PHI[j+1][i+1]
double P[NDY][NDX], A[5][NDY][NDX], B[NDY][NDX] ;
double U[NDY][NDX], V[NDY][NDX], PHI[NDY][NDX] ;
int NX, NY, NX1, NX2, NY1, NY2 ;
FILE *FT1 ;
// Function Definition
int PSOL(double hx, double hy, double dt, double Re, double OMG,
double EPS, int ITER) ;
int PHISOL(double hx, double hy, double OMG, double EPS, int ITER) ;
void OUTPUT(double hx, double hy) ;
// Main Program
int main()
{ int i, j, k, NT, MT, ID, T1, T2 ;
int ITER, NTER, ILP, NI1, NI2, IER ;
double Re, dt, T, hx, hy, dx, dy, ddx, ddy, EPS, OMG ;
double XM ,YM, OVER, CPU ;
char FLOW3[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 ;
XM = 2.0 ;
YM = 1.0 ;
printf("Type In NX,NY \n") ;
scanf("%d %d",&NX,&NY) ;
printf("Type in Re,DT,NT(Total NO.),MT(Output Interval) \n") ;
scanf("%lf %lf %d %d",&Re,&dt,&NT,&MT) ;
// Initial Data
if(NX >= NDX-1) { NX = NDX - 2 ; }
if(NY >= NDY-1) { NY = NDY - 2 ; }
NX1 = NX*2/10 ;
NX2 = NX*3/10 ;
NY1 = NY*2/5 ;
NY2 = NY*3/5 ;
printf("Flow around Square a Pillar Numerical Analysis \n") ;
printf("XM=%f NX1,NX2,NX=%d %d %d \n",XM,NX1,NX2,NX) ;
printf("YM=%f NY1,NY2,NY=%d %d %d \n",YM,NY1,NY2,NY) ;
printf("Re=%f dt=%f NT,MT=%d %d \n",Re,dt,NT,MT) ;
// Initial Constant
IER = 0 ;
T = 0.0 ;
hx = XM/NX ;
hy = YM/NY ;
dx = 1.0/(2.0*hx) ;
dy = 1.0/(2.0*hy) ;
ddx = 1.0/(hx*hx) ;
ddy = 1.0/(hy*hy) ;
EPS = 1.0e-4 ;
OMG = 1.0 + log((NX+NY)*1.0)/log((NDX+NDY)*1.0) ;
// Set Initial U,V
ILP = NX*NY ;
ITER = PHISOL(hx, hy, OMG, EPS, ILP) ;
// Set Initial P
for (j=-1; j<=NY+1; j++) {
for (i=-1; i<=NX+1; i++) {
P(i,j) = 0.0 ; }
}
// Main Loop
T1 = clock() ;
NTER = ITER ;
for (k=1; k<=NT; k++) {
T = T + dt ;
// Compute P
ILP = NX*NY ;
ITER = PSOL(hx, hy, dt, Re, OMG, EPS, ILP) ;
NTER = NTER + ITER ;
// Compute U,V
OVER = 0.0 ;
// Lower-Y Area
for (j=1; j<NY1; j++) {
for (i=1; i<=NX; i++) {
#include "UV.h"
} }
// Middle-Y Area
for (j=NY1; j<=NY2; j++) {
for (i=1; i<NX1; i++) {
#include "UV.h"
}
for (i=NX2+1; i<=NX; i++) {
#include "UV.h"
} }
// Upper-Y Area
for (j=NY2+1; j<NY; j++) {
for (i=1; i<=NX; i++) {
#include "UV.h"
} }
// Check Divergent
if(OVER > 1.0e10) {
printf("** Stop for over flow computation ** \n") ;
printf(" You have to give smaller dt \n") ;
IER = 1 ;
goto M100 ; }
// Output U,V,Phi
if( k%MT == 0) {
NI1 = ID/10 ;
NI2 = ID - NI1*10 ;
FLOW3[6] = NTABL[NI1] ;
FLOW3[7] = NTABL[NI2] ;
FT1 = fopen(FLOW3,"w") ;
fprintf(FT1,"# Flow around Square a Pillar Numerical Analysis \n") ;
fprintf(FT1,"XM=%f NX1,NX2,NX=%d %d %d \n",NX1,NX2,NX) ;
fprintf(FT1,"YM=%f NY1,NY2,NY=%d %d %d \n",NY1,NY2,NY) ;
fprintf(FT1,"Re=%f dt=%f NT,MT=%d %d \n",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 ;
OUTPUT(hx, hy) ;
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 PSOL(double hx, double hy, double dt, double Re, double OMG,
double EPS, int ITER)
//=================================================================C
// Solve Ax=b by SOR with 2 dimensional FDM C
// Solve P by U,V with -div(P)=f(U,V) C
//-----------------------------------------------------------------C
// P,A,B,U,V ; Global array C
// hx R*8, In. Delta x C
// hy R*8, In, Delta y C
// dt R*8, In, Delta t C
// Re R*8, In, Re Number C
// OMG R*8, In, Acceleration parameter for SOR C
// EPS R*8, In, if ||r||/||b|| <= EPS --> return C
// ITER I*4, In, Maximum Number of Iteration C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/09/14 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
{ int i, j, k, ILP ;
double DYX, DXY, DHX, DHY, BN, RN, ERR, A1, A2, A3, A4, A5 ;
double R, UP0, UM0, U0P, U0M, VP0, VM0, V0P, V0M ;
// Set Value
DYX = hy/hx ;
DXY = hx/hy ;
DHX = 1.0/(2.0*hx) ;
DHY = 1.0/(2.0*hy) ;
// Set Outer Boundary U,V,P
for (i=-1; i<=NX+1; i++) {
P(i,-1) = P(i,1) ;
P(i,NY+1) = P(i,NY-1) ;
U(i,-1) = U(i,0) ;
V(i,-1) = V(i,0) ;
U(i,NY+1) = U(i,NY) ;
V(i,NY+1) = V(i,NY) ; }
for (j=-1; j<=NY+1; j++) {
// P=0 on L3
P(NX,j) = 0.0 ;
P(-1,j) = P(1,j) ;
P(NX+1,j) = P(NX-1,j) ;
U(-1,j) = U(0,j) ;
V(-1,j) = V(0,j) ;
U(NX+1,j) = U(NX,j) ;
V(NX+1,j) = V(NX,j) ; }
// Set A and B
BN = 0.0 ;
for (j=0; j<=NY; j++) {
for (i=0; i<NX; i++) {
// Inner
A1 = -DXY ;
A2 = -DYX ;
A3 = 2.0*(DXY + DYX) ;
A4 = -DYX ;
A5 = -DXY ;
// on B1
if(i == NX1 & j >= NY1 & j <= NY2) {
A3 = A3 - DYX ;
A4 = 0.0 ;
UP0 = 0.0 ;
VP0 = 0.0 ; }
else {
UP0 = U(i+1,j) ;
VP0 = V(i+1,j) ; }
// on B2
if(j == NY1 & i >= NX1 & i <= NX2) {
A3 = A3 - DXY ;
A5 = 0.0 ;
U0P = 0.0 ;
V0P = 0.0 ; }
else {
U0P = U(i,j+1) ;
V0P = V(i,j+1) ; }
// on B3
if(i == NX2 & j >= NY1 & j <= NY2) {
A2 = 0.0 ;
A3 = A3 - DYX ;
UM0 = 0.0 ;
VM0 = 0.0 ; }
else {
UM0 = U(i-1,j) ;
VM0 = V(i-1,j) ; }
// on B4
if(j == NY2 & i >= NX1 & i <= NX2) {
A1 = 0.0 ;
A3 = A3 - DXY ;
U0M = 0.0 ;
V0M = 0.0 ; }
else {
U0M = U(i,j-1) ;
V0M = V(i,j-1) ; }
// Set A,B
A(i,j,1) = A1 ;
A(i,j,2) = A2 ;
A(i,j,3) = A3 ;
A(i,j,4) = A4 ;
A(i,j,5) = A5 ;
B(i,j) = ((UP0-UM0)*(UP0-UM0)*DYX + (V0P-V0M)*(V0P-V0M)*DXY )/4.0
+ (VP0-VM0)*(U0P-U0M)/2.0
- ((UP0-UM0)*hy + (V0P-V0M)*hx)/(2.0*dt) ;
// in Square a Pillar
if( (i > NX1 & i < NX2) & (j > NY1 & j < NY2) ) {
A(i,j,1) = 0.0 ;
A(i,j,2) = 0.0 ;
A(i,j,3) = 1.0 ;
A(i,j,4) = 0.0 ;
A(i,j,5) = 0.0 ;
B(i,j) = P(NX1,NY1) ; }
BN = BN + B(i,j)*B(i,j)/(A(i,j,3)*A(i,j,3)) ;
} }
// Solve Ax=b
for (k=1; k<=ITER; k++) {
RN = 0.0 ;
for (j=0; j<=NY; j++) {
for (i=0; i<NX; i++) {
R = ( B(i,j) - A(i,j,1)*P(i,j-1) - A(i,j,2)*P(i-1,j)
- A(i,j,4)*P(i+1,j) - A(i,j,5)*P(i,j+1) )/A(i,j,3) - P(i,j) ;
P(i,j) = P(i,j) + OMG*R ;
RN = RN + R*R ;
} }
// Check Conversion
ERR = sqrt(RN/BN) ;
// printf("k=%d ERR=%e RN=%e BN=%e \n",k,ERR,RN,BN) ;
if(ERR <= EPS) goto M100 ;
}
M100: ILP = k ;
return (ILP) ;
}
//=================================================================C
int PHISOL(double hx, double hy, double OMG, double EPS, int ITER)
//=================================================================C
// Solve Ax=b by SOR with 2 dimensional FDM C
// Solve PHI and Compute U,V C
// -div(PHI)=0 and d(PHI)/dx=U, d(PHI)/dy=V C
//-----------------------------------------------------------------C
// PHI,A,B,U,V ; Global array C
// hx R*8, In. Delta x C
// hy R*8, In, Delta y C
// OMG R*8, In, Acceleration parameter for SOR C
// EPS R*8, In, if ||r||/||b|| <= EPS --> return C
// ITER I*4, In, Number of Iteration C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/09/14 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
{ int i, j, k, ILP ;
double DYX, DXY, DHX, DHY, A1, A2, A3, A4, A5, BV ;
double BN, RN, R, ERR ;
// Set Value
DYX = hy/hx ;
DXY = hx/hy ;
DHX = 1.0/(2.0*hx) ;
DHY = 1.0/(2.0*hy) ;
// Set Initial PHI=X
for (j=0; j<=NY; j++) {
for (i=0; i<=NX+1; i++) {
PHI(i,j) = i*hx ;
} }
// Set A and B
BN = 0.0 ;
for (j=1; j<NY; j++) {
for (i=1; i<=NX; i++) {
// Inner
A1 = -DXY ;
A2 = -DYX ;
A3 = 2.0*(DXY + DYX) ;
A4 = -DYX ;
A5 = -DXY ;
BV = 0.0 ;
// on L3
if(i == NX) {
A3 = A3 - DYX ;
A4 = 0.0 ;
BV = BV + hy ; }
// on B1
if(i == NX1 & j >= NY1 & j <= NY2) {
A3 = A3 - DYX ;
A4 = 0.0 ; }
// on B2
if(j == NY1 & i >= NX1 & i <= NX2) {
A3 = A3 - DXY ;
A5 = 0.0 ; }
// on B3
if(i == NX2 & j >= NY1 & j <= NY2) {
A2 = 0.0 ;
A3 = A3 - DYX ; }
// on B4
if(j == NY2 & i >= NX1 & i <= NX2) {
A1 = 0.0 ;
A3 = A3 - DXY ; }
// in Square a Pillar
if( (i > NX1 & i < NX2) & (j > NY1 & j < NY2) ) {
A1 = 0.0 ;
A2 = 0.0 ;
A3 = 1.0 ;
A4 = 0.0 ;
A5 = 0.0 ;
BV = 0.0 ; }
// Set A,B
A(i,j,1) = A1 ;
A(i,j,2) = A2 ;
A(i,j,3) = A3 ;
A(i,j,4) = A4 ;
A(i,j,5) = A5 ;
B(i,j) = BV ;
BN = BN + BV*BV/(A3*A3) ;
} }
// Solve Ax=b
for (k=1; k<=ITER; k++) {
RN = 0.0 ;
for (j=1; j<NY; j++) {
for (i=1; i<=NX; i++) {
R = ( B(i,j) - A(i,j,1)*PHI(i,j-1) - A(i,j,2)*PHI(i-1,j)
- A(i,j,4)*PHI(i+1,j)-A(i,j,5)*PHI(i,j+1) )/A(i,j,3) - PHI(i,j) ;
PHI(i,j) = PHI(i,j) + OMG*R ;
RN = RN + R*R ;
} }
// Check Conversion
ERR = sqrt(RN/BN) ;
if(ERR <= EPS) goto M100 ;
// printf("kk=%d ERR=%e RN=%e BN=%e \n",k,ERR,RN,BN) ;
}
M100: ILP = k ;
// Compute U,V
for (j=1; jNY; j++) {
for(i=1; i<=NX; i++) {
U(i,j) = (PHI(i+1,j) - PHI(i-1,j))*DHX ;
V(i,j) = (PHI(i,j+1) - PHI(i,j-1))*DHY ;
} }
// in Square a Pillar
for (j=NY1; j<=NY2; j++) {
for (i=NX1; i<=NX2; i++) {
U(i,j) = 0.0 ;
V(i,j) = 0.0 ;
} }
// on L1 Boundary
for (j=1; j<NY; j++) {
U(0,j) = 1.0 ;
V(0,j) = 0.0 ; }
// on L2,L4 Boundary
for (i=0; i<=NX; i++) {
U(i,0) = 1.0 ;
V(i,0) = 0.0 ;
U(i,NY) = 1.0 ;
V(i,NY) = 0.0 ; }
return (ILP) ;
}
//=================================================================C
void OUTPUT(double hx, double hy)
//=================================================================C
// Compute PHI and Output PHI,U,V C
// OMG(i,j) = dV/dx - dU/dy C
//-----------------------------------------------------------------C
// hx R*8, In, Delta x C
// hy R*8, In, Delta y C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/09/14 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
{ int i, j ;
double X, Y ;
// Set Omega
for (j=0; j<=NY; j++) {
for (i=0; i<=NX; i++) {
OMG(i,j) = (V(i+1,j)-V(i-1,j))/hx - (U(i,j+1)-U(i,j-1))/hy ;
} }
// Output
for (j=0; j<=NY; j++) {
Y = j*hy ;
for (i=0; i<=NX; i++) {
X = i*hx ;
fprintf(FT1,"%10.4f %10.4f %10.4f %10.4f %10.4f \n",
X,Y,OMG(i,j),U(i,j),V(i,j)) ; }
}
}
3. include file（ファイル名；UV.h)
U(i,j) = U(i,j) - dt*((U(i+1,j)*U(i+1,j) - U(i-1,j)*U(i-1,j))*dx
+ (U(i,j+1)*V(i,j+1) - U(i,j-1)*V(i,j-1))*dy
+ (P(i+1,j) - P(i-1,j))*dx
+ ( (2.0*U(i,j) - U(i-1,j) - U(i+1,j))*ddx
+ (2.0*U(i,j) - U(i,j-1) - U(i,j+1))*ddy )/Re ) ;
V(i,j) = V(i,j) - dt*((V(i,j+1)*V(i,j+1) - V(i,j-1)*V(i,j-1))*dy
+ (U(i+1,j)*V(i+1,j) - U(i-1,j)*V(i-1,j))*dx
+ (P(i,j+1) - P(i,j-1))*dy
+ ( (2.0*V(i,j) - V(i-1,j) - V(i+1,j))*ddx
+ (2.0*V(i,j) - V(i,j-1) - V(i,j+1))*ddy )/Re ) ;
OVER = OVER + fabs(U(i,j)) + fabs(V(i,j)) ;