C computation code and MATLAB posttreatment code: Immersed Boundary Method for fluid flows passing through a cylinder
Date:
//Introduction: The first-order accurate pressure projection method is used to update solutions. The intermediate velocities are obtained by an explicit Euler scheme. The pressure Poisson equation is iteratively solved with a linear Multigrid algorithm. The advection term in NS equations is approximated with a first-order upwind scheme and the viscosity term is approximated with a five-point central difference scheme. The IBM force model in “ An Immersed Boundary Method with Formal Second-Order Accuracy and Reduced Numerical Viscosity “ is adopted.
The C codes are as follows:
/***************** Immersed Boundary Method for the incompressible Navier-Stokes Equations: flow passing a cylinder *********************/
#include
#include
#include
#include
#include “ns.h”
int nx, ny, nxy, n_level, p_relax, nt;
double **fxx, **fyy, **bd_posi, **ebd_posi, **bd_velo, **bd_force, **fx, **fy, pi, **ux, **uy, **sor,
h, h2, **tu, **tv, **workp, **worku, **workv, **adv_u, **adv_v, dt, xleft, xright, yleft, yright, Re, We;
int main()
{
extern int nx, ny, nxy, n_level, p_relax, nt;
extern double **fxx, **fyy, **bd_posi, **ebd_posi, **bd_velo, **bd_force, **fx, **fy, pi, **ux, **uy, **sor,
h, h2, **tu, **tv, **workp, **worku, **workv, **adv_u, **adv_v, dt, xleft, xright, yleft, yright, Re, We;
int ia, ja, k, i, j, it, max_it, ns, count = 1;
double **u, **v, **nu, **nv, **p, **c, x, y;
FILE *fu, *fv, *fp, *fuv, *fvv, *fxy, *my, *my2;
pi = 4.0*atan(1.0);
p_relax = 3;
nx = gnx;
ny = gny;
nxy = nx * ny;
n_level = (int)(log(nx)/log(2)-0.9);
xleft=0.0, xright= 32.0;
yleft=0.0, yright= 8.0;
/********/
h = xright / (double)nx;
h2 = pow(h, 2);
dt = 2.0hh;
/********/
c = dmatrix(0, nx+1, 0, ny+1);
fx = dmatrix(1, nx, 1, ny);
fy = dmatrix(1, nx, 1, ny);
fxx = dmatrix(0, nx, 1, ny);
fyy = dmatrix(1, nx, 0, ny);
ux = dmatrix(0, nx+1, 0, ny+1);
uy = dmatrix(0, nx+1, 0, ny+1);
workp = dmatrix(0, nx+1, 0, ny+1);
worku = dmatrix(0, nx+1, 0, ny+1);
workv = dmatrix(0, nx+1, 0, ny+1);
u = dmatrix(-1, nx+1, 0, ny+1);
v = dmatrix(0, nx+1, -1, ny+1);
tu = dmatrix(0, nx, 1, ny);
tv = dmatrix(1, nx, 0, ny);
p = dmatrix(0, nx, 0, ny);
nu = dmatrix(-1, nx+1, 0, ny+1);
nv = dmatrix(0, nx+1, -1, ny+1);
adv_u = dmatrix(0, nx, 1, ny);
adv_v = dmatrix(1, nx, 0, ny);
sor = dmatrix(1, nx, 1, ny);
zero_matrix(fx, 1, nx, 1, ny);
zero_matrix(fy, 1, nx, 1, ny);
zero_matrix(fxx, 0, nx, 1, ny);
zero_matrix(fyy, 1, nx, 0, ny);
max_it = 30720;
ns = 1280; //max_it;
nt = 100;
Re = 100.0;
We = 1.0;
printf(“nx = %d , ny = %d\n”, nx, ny);
printf(“dt = %f\n”,dt);
printf(“max_it = %d\n”, max_it);
printf(“ns = %d\n”, ns);
printf(“n_level = %d\n\n”, n_level);
printf(“Reynolds number = %f \n”,Re);
bd_posi = dmatrix(0, nt+1, 1, 2);
ebd_posi = dmatrix(0, nt+1, 1, 2);
bd_velo = dmatrix(1, nt, 1, 2);
bd_force = dmatrix(1, nt, 1, 2);
fu = fopen(“u.m”,”w”);
fv = fopen(“v.m”,”w”);
fuv = fopen(“uv.m”,”w”);
fvv = fopen(“vv.m”,”w”);
fp = fopen(“p.m”,”w”);
fxy = fopen(“bd_posi.m”,”w”);
fprintf(fu,”\n”);
fprintf(fv,”\n”);
fprintf(fuv,”\n”);
fprintf(fvv,”\n”);
fprintf(fp,”\n”);
fprintf(fxy,”\n”);
fclose(fu);
fclose(fv);
fclose(fuv);
fclose(fvv);
fclose(fp);
fclose(fxy);
initialization(p, u, v, bd_posi);
print_data1(u, v, p);
print_data2(bd_posi);
mat_copy(nu, u, 0, nx, 1, ny);
mat_copy(nv, v, 1, nx, 0, ny);
for (it=1; it<=max_it; it++) {
boundary_force(bd_posi, bd_force);
fluid_force(bd_posi, bd_force);
full_step(u, v, nu, nv, p);
mat_copy(u, nu, 0, nx, 1, ny);
mat_copy(v, nv, 1, nx, 0, ny);
move_boundary(bd_posi, bd_velo, u, v);
if (it % ns==0) {
count++;
print_data1(nu, nv, p);
printf("print out counts %d \n",count);
}
if (it % ns==0) {
print_data2(bd_posi); }
printf(" %d \n",it);
}
fu = fopen(“u.m”,”a”);
fv = fopen(“v.m”,”a”);
fuv = fopen(“uv.m”,”a”);
fvv = fopen(“vv.m”,”a”);
fp = fopen(“p.m”,”a”);
fxy = fopen(“bd_posi.m”,”a”);
fprintf(fu,” \n”);
fprintf(fv,”\n”);
fprintf(fuv,” \n”);
fprintf(fvv,”\n”);
fprintf(fp,”\n”);
fprintf(fxy,”“,nt,nt);
fclose(fu);
fclose(fv);
fclose(fuv);
fclose(fvv);
fclose(fp);
fclose(fxy);
return 0;
}
void initialization(double **p, double **u, double **v, double **bd_posi)
{
extern int nx, ny, nt;
extern double ds, pi, h, xright, yright, **ebd_posi;
int i, j, k;
double x, y, the, aa, r, a1=1.0, ep1=0.2*a1;
FILE *my;
for (k=1; k<=nt; k++) {
the = 2.0*pi*k/nt;
r = a1; // + ep1*cos(2.0*the);
bd_posi[k][1] = r*cos(the) + 7.0;
bd_posi[k][2] = r*sin(the) + 4.0;
}
bd_posi[0][1] = bd_posi[nt][1];
bd_posi[nt+1][1] = bd_posi[1][1];
bd_posi[0][2] = bd_posi[nt][2];
bd_posi[nt+1][2] = bd_posi[1][2];
//define fixed solid interface points
for (k = 0; k <= nt+1; k++) {
ebd_posi[k][1] = bd_posi[k][1];
ebd_posi[k][2] = bd_posi[k][2];
}
ijloop {
p[i][j] = 0.0;
}
i0jloop {
y = ((double)j - 0.5)*h;
u[i][j] = 1.0;
}
ij0loop {
v[i][j] = 0.0;
}
jloop{
u[0][j] = 1.0;
}
ijloop{
x = ((double)i - 0.5)*h;
y = ((double)j - 0.5)*h;
if( sqrt( pow(x-7.0,2) + pow(y-4.0,2) ) <= r)
u[i][j] = 0.0;
}
}
void boundary_force(double **bd_posi, double **bd_force)
{
extern int nt;
extern double We, **ebd_posi;
int k;
double fac, kappa, ds1, ds2, x1, y1, x2, y2;
fac = 1.0;
kappa = 9600.0;
for (k=1; k<=nt; k++) {
x1 = bd_posi[k+1][1] - bd_posi[k][1];
y1 = bd_posi[k+1][2] - bd_posi[k][2];
x2 = bd_posi[k][1] - bd_posi[k-1][1];
y2 = bd_posi[k][2] - bd_posi[k-1][2];
ds1 = sqrt(x1*x1 + y1*y1);
ds2 = sqrt(x2*x2 + y2*y2);
//solid boundary force
bd_force[k][1] = kappa*(ebd_posi[k][1] - bd_posi[k][1]);
bd_force[k][2] = kappa*(ebd_posi[k][2] - bd_posi[k][2]);
}
}
void fluid_force(double **bd_posi, double **bd_force)
{
extern int nx, ny, nt;
extern double **fx, **fy, **fxx, **fyy, h, h2;
int i, j, k, ia, ja;
double x1, x2, y1, y2, ds1, ds2, d, d1, a;
FILE *my;
zero_matrix(fx, 1, nx, 1, ny);
zero_matrix(fy, 1, nx, 1, ny);
for (k=1; k<=nt; k++) {
x1 = bd_posi[k+1][1] - bd_posi[k][1];
y1 = bd_posi[k+1][2] - bd_posi[k][2];
x2 = bd_posi[k][1] - bd_posi[k-1][1];
y2 = bd_posi[k][2] - bd_posi[k-1][2];
ds1 = sqrt(x1*x1 + y1*y1);
ds2 = sqrt(x2*x2 + y2*y2);
ia = (int)(bd_posi[k][1]/h)-1;
ja = (int)(bd_posi[k][2]/h)-1;
for (i=ia; i<=ia+3; i++) {
d1 = delta((((double)i-0.5)*h -bd_posi[k][1])/h)/h2;
for (j=ja; j<=ja+3; j++) {
d = d1*delta((((double)j-0.5)*h-bd_posi[k][2])/h);
fx[i][j] += d*bd_force[k][1]*0.5*(ds1+ds2);
fy[i][j] += d*bd_force[k][2]*0.5*(ds1+ds2);
}
}
}
ijloop {
if (i==nx)
fxx[0][j] = fxx[nx][j] = 0.5*(fx[nx][j]+fx[1][j]);
else
fxx[i][j] = 0.5*(fx[i][j]+fx[i+1][j]);
if (j==ny)
fyy[i][0] = fyy[i][ny] = 0.5*(fy[i][ny]+fy[i][1]);
else
fyy[i][j] = 0.5*(fy[i][j]+fy[i][j+1]);
}
}
void move_boundary(double **bd_posi, double **bd_velo, double **u, double **v)
{
extern int nx, ny, nt;
extern double h, dt;
int i, j, k, ia, ja;
double d, d1;
for (k=1; k<=nt; k++) {
bd_velo[k][1] = 0.0;
bd_velo[k][2] = 0.0;
ia = (int)(bd_posi[k][1]/h)-1;
ja = (int)(bd_posi[k][2]/h)-1;
for (i=ia; i<=ia+3; i++) {
d = delta((((double)i-0.5)*h -bd_posi[k][1])/h);
for (j=ja; j<=ja+3; j++) {
d1 = d*delta((((double)j-0.5)*h-bd_posi[k][2])/h);
bd_velo[k][1] = bd_velo[k][1] + 0.5*(u[i][j]+u[i-1][j])*d1;
bd_velo[k][2] = bd_velo[k][2] + 0.5*(v[i][j]+v[i][j-1])*d1;
}
}
bd_posi[k][1] = bd_posi[k][1] + dt*bd_velo[k][1];
bd_posi[k][2] = bd_posi[k][2] + dt*bd_velo[k][2];
}
bd_posi[0][1] = bd_posi[nt][1];
bd_posi[nt+1][1] = bd_posi[1][1];
bd_posi[0][2] = bd_posi[nt][2];
bd_posi[nt+1][2] = bd_posi[1][2];
}
double delta(double x)
{
double value;
if (fabs(x) > 2.0)
value=0.0;
else if (fabs(x) < 1.0)
value=(3.0-2.0*fabs(x)+sqrt(1.0+4.0*fabs(x)-4.0*x*x))/8.0;
else
value=0.5-(3.0-2.0*(2.0-fabs(x))+sqrt(1.0+4.0*(2.0-fabs(x))-4.0*pow(2.0-fabs(x),2)))/8.0;
return value;
}
void full_step(double **u, double **v, double **nu, double **nv, double **p)
{
extern int nx, ny;
extern double dt, **worku, **workv, **adv_u, **adv_v, **tu, **tv;
int i, j;
FILE *my;
advection_step(u, v, adv_u, adv_v);
temp_uv(tu, tv, u, v, adv_u, adv_v);
Poisson(tu, tv, p);
grad_p(p, worku, workv, nx, ny);
i0jloop {
nu[i][j] = tu[i][j] - dt*worku[i][j];
}
ij0loop {
nv[i][j] = tv[i][j] - dt*workv[i][j];
}
}
void advection_step(double **u, double **v, double **adv_u, double **adv_v)
{
extern int nx, ny;
extern double h;
int i, j;
augmenuv(u, v, nx, ny);
i0jloop {
if (u[i][j]>0.0)
adv_u[i][j] = u[i][j]*(u[i][j]-u[i-1][j])/h;
else
adv_u[i][j] = u[i][j]*(u[i+1][j]-u[i][j])/h;
if (v[i][j-1]+v[i+1][j-1]+v[i][j]+v[i+1][j]>0.0)
adv_u[i][j] += 0.25*(v[i][j-1]+v[i+1][j-1]+v[i][j]+v[i+1][j])*(u[i][j]-u[i][j-1])/h;
else
adv_u[i][j] += 0.25*(v[i][j-1]+v[i+1][j-1]+v[i][j]+v[i+1][j])*(u[i][j+1]-u[i][j])/h;
}
ij0loop {
if (u[i-1][j]+u[i][j]+u[i-1][j+1]+u[i][j+1]>0.0)
adv_v[i][j] = 0.25*(u[i-1][j]+u[i][j]+u[i-1][j+1]+u[i][j+1])*(v[i][j]-v[i-1][j])/h;
else
adv_v[i][j] = 0.25*(u[i-1][j]+u[i][j]+u[i-1][j+1]+u[i][j+1])*(v[i+1][j]-v[i][j])/h;
if (v[i][j]>0.0)
adv_v[i][j] += v[i][j]*(v[i][j] - v[i][j-1])/h;
else
adv_v[i][j] += v[i][j]*(v[i][j+1] - v[i][j])/h;
}
}
void temp_uv(double **tu, double **tv, double **u, double **v, double **adv_u, double **adv_v)
{
extern int nx, ny;
extern double h, Re, **fxx, **fyy, We;
int i, j;
FILE *my;
augmenuv(u, v, nx, ny);
i0jloop {
tu[i][j] = u[i][j] + dt*( (u[i+1][j]+u[i-1][j]-4.0*u[i][j]+u[i][j+1]+u[i][j-1])/(Re*h*h)
-adv_u[i][j] + fxx[i][j]/We);
}
ij0loop {
tv[i][j] = v[i][j] + dt*( (v[i+1][j]+v[i-1][j]-4.0*v[i][j]+v[i][j+1]+v[i][j-1])/(Re*h*h)
-adv_v[i][j] + fyy[i][j]/We);
}
}
void Poisson(double **tu, double **tv, double **p)
{
extern int nx, ny;
extern double **workp, h;
int i, j;
FILE *my;
source_uv(tu, tv, workp, nx, ny);
MG_Poisson(p, workp);
}
void MG_Poisson(double **p, double **f)
{
extern int nx, ny;
extern double **sor, **workv;
int i, j, max_it = 2000, it_mg = 1;
double tol = 1.0e-5, resid = 1.0;
FILE *my;
mat_copy(workv, p, 1, nx, 1, ny);
while (it_mg <= max_it && resid >= tol) {
vcycle_uv(p, f, nx, ny, 1);
pressure_update(p);
ijloop
sor[i][j] = workv[i][j] - p[i][j];
resid = mat_max(sor, 1, nx, 1, ny);
mat_copy(workv, p, 1, nx, 1, ny);
it_mg++;
}
printf("Mac iteration = %d residual = %16.15f \n",it_mg, resid);
return;
}
void vcycle_uv(double **uf, double **ff, int nxf, int nyf, int ilevel)
{
extern int n_level;
relax(uf, ff, nxf, nyf);
if (ilevel < n_level) {
int nxc, nyc;
double **rf, **uc, **fc;
nxc=nxf / 2;
nyc=nyf / 2;
rf=dmatrix(1, nxf, 1, nyf);
uc=dmatrix(1, nxc, 1, nyc);
fc=dmatrix(1, nxc, 1, nyc);
residual(rf, uf, ff, nxf, nyf);
restrict(rf, fc, nxc, nyc);
zero_matrix(uc, 1, nxc, 1, nyc);
vcycle_uv(uc, fc, nxc, nyc, ilevel + 1);
prolong(uc, rf, nxc, nyc);
mat_add(uf, uf, rf, 1, nxf, 1, nyf);
relax(uf, ff, nxf, nyf);
free_dmatrix(rf, 1, nxf, 1, nyf);
free_dmatrix(uc, 1, nxc, 1, nyc);
free_dmatrix(fc, 1, nxc, 1, nyc);
}
}
void relax(double **p, double **f, int nxt, int nyt)
{
extern int p_relax, nx;
extern double xright;
int i, j, iter;
double ht, ht2, coef, src;
FILE *my;
ht = xright / (double) nxt;
ht2 = pow(ht, 2);
// coef = -0.0/ht2;
for (iter=1; iter<=p_relax; iter++) {
ijloopt {
coef = 0.0;
src = f[i][j];
if (i==1) {
src -= (p[2][j] /*+ p[nxt][j]*/)/ht2;
coef += -1.0/ht2;
}
else if (i==nxt) {
src -= (p[nxt-1][j] /*+ p[1][j]*/)/ht2;
coef += -1.0/ht2;
}
else {
src -= (p[i+1][j] + p[i-1][j])/ht2;
coef += -2.0/ht2;
}
if (j==1){
src -= (p[i][2] )/ht2;
coef += -1.0/ht2; }
else if (j==nyt) {
src -= (p[i][nyt-1] )/ht2;
coef += -1.0/ht2; }
else {
src -= (p[i][j+1] + p[i][j-1])/ht2;
coef += -2.0/ht2; }
p[i][j] = src / coef;
}
}
}
void restrict(double **u_fine, double **u_coarse, int nxt, int nyt)
{
int i, j;
ijloopt {
u_coarse[i][j] = 0.25*(u_fine[2*i-1][2*j-1]
+u_fine[2*i-1][2*j]
+u_fine[2*i][2*j-1]
+u_fine[2*i][2*j]);
}
}
void prolong(double **u_coarse, double **u_fine, int nxt, int nyt)
{
int i, j;
ijloopt {
u_fine[2*i-1][2*j-1] = u_fine[2*i-1][2*j] =
u_fine[2*i][2*j-1] = u_fine[2*i][2*j] = u_coarse[i][j];
}
}
void grad_p(double **p, double **dpdx, double **dpdy, int nxt, int nyt)
{
extern int nx;
extern double xright, Re, **pbx;
int i, j;
double ht;
ht = xright / (double) nxt;
i0jloopt {
if (i==0)
dpdx[0][j] = 0.0; //(p[1][j] - p[nxt][j])/ht;
else if (i==nxt)
dpdx[nxt][j] = 0.0; //(p[1][j] - p[nxt][j])/ht;
else
dpdx[i][j] = (p[i+1][j] - p[i][j])/ht;
}
ij0loopt {
if (j==0)
dpdy[i][0] = 0.0;
dpdy[i][nyt] = 0.0;
else
dpdy[i][j] = (p[i][j+1] - p[i][j])/ht;
}
}
void source_uv(double **tu, double **tv, double **divuv, int nxt, int nyt)
{
extern double dt;
int i, j;
FILE *my;
div_uv(tu, tv, divuv, nxt, nyt);
ijloopt {
divuv[i][j] /= dt;
}
}
void div_uv(double **tu, double **tv, double **divuv, int nxt, int nyt)
{
extern double xright;
int i, j;
double ht;
ht = xright / (double) nxt;
ijloopt {
divuv[i][j] = (tu[i][j] - tu[i-1][j] + tv[i][j] - tv[i][j-1])/ht;
}
}
void residual(double **r, double **u, double **f, int nxt, int nyt)
{
laplace(u, r, nxt, nyt);
mat_sub(r, f, r, 1, nxt, 1, nyt);
}
void laplace(double **p, double **lap_p, int nxt, int nyt)
{
double **dpdx, **dpdy;
dpdx = dmatrix(0, nxt, 1, nyt);
dpdy = dmatrix(1, nxt, 0, nyt);
grad_p(p, dpdx, dpdy, nxt, nyt);
div_uv(dpdx, dpdy, lap_p, nxt, nyt);
free_dmatrix(dpdx, 0, nxt, 1, nyt);
free_dmatrix(dpdy, 1, nxt, 0, nyt);
}
double *dvector (long nl, long nh)
{
double *v;
v=(double ) malloc((nh-nl+1+1)sizeof(double));
return v-nl+1;
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
{
double **m;
long i, nrow=nrh-nrl+1+1, ncol=nch-ncl+1+1;
m=(double *) malloc((nrow)sizeof(double*));
m+=1;
m-=nrl;
m[nrl]=(double ) malloc((nrowncol)*sizeof(double));
m[nrl]+=1;
m[nrl]-=ncl;
for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;
return m;
}
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
{
free(m[nrl]+ncl-1);
free(m+nrl-1);
return;
}
void mat_add(double **a, double **b, double **c,
int xl, int xr, int yl, int yr)
{
int i, j;
for (i=xl; i<=xr; i++)
for (j=yl; j<=yr; j++){
a[i][j] = b[i][j]+c[i][j];
}
return;
}
void mat_add2(double **a, double **b, double **c,
double **a2, double **b2, double **c2,
int xl, int xr, int yl, int yr)
{
int i, j;
for (i=xl; i<=xr; i++)
for (j=yl; j<=yr; j++){
a[i][j] = b[i][j]+c[i][j];
a2[i][j] = b2[i][j]+c2[i][j];
}
return;
}
void zero_vector(double *a, int xl, int xr)
{
int i, j;
for (i=xl; i<=xr; i++)
a[i] = 0.0;
return;
}
void zero_matrix(double **a, int xl, int xr, int yl, int yr)
{
int i, j;
for (i=xl; i<=xr; i++)
for (j=yl; j<=yr; j++){
a[i][j] = 0.0;
}
return;
}
void mat_copy(double **a, double **b,
int xl, int xr, int yl, int yr)
{
int i, j;
for (i=xl; i<=xr; i++)
for (j=yl; j<=yr; j++)
a[i][j]=b[i][j];
return;
}
void mat_copy2(double **a, double **b,
double **a2, double **b2,
int xl, int xr, int yl, int yr)
{
int i, j;
for (i=xl; i<=xr; i++)
for (j=yl; j<=yr; j++){
a[i][j]=b[i][j];
a2[i][j]=b2[i][j];
}
return;
}
void mat_sub(double **a, double **b, double **c,
int nrl, int nrh, int ncl, int nch)
{
int i,j;
for (i=nrl; i<=nrh; i++)
for (j=ncl; j<=nch; j++)
a[i][j]=b[i][j]-c[i][j];
return;
}
void mat_sub2(double **a, double **b, double **c,
double **a2, double **b2, double **c2,
int nrl, int nrh, int ncl, int nch)
{
int i,j;
for (i=nrl; i<=nrh; i++)
for (j=ncl; j<=nch; j++) {
a[i][j]=b[i][j]-c[i][j];
a2[i][j]=b2[i][j]-c2[i][j];
}
return;
}
double mat_max(double **a,
int nrl, int nrh, int ncl, int nch)
{
int i, j;
double x = 0.0;
for(i = nrl; i <= nrh; i++)
for(j = ncl; j <= nch; j++){
if (fabs(a[i][j]) > x)
x = fabs(a[i][j]);
}
return x;
}
void print_data1(double **u, double **v, double **p)
{
extern int nx, ny, nt;
int i, j, k;
FILE *fu, *fv, *fp, *fuv, *fvv;
fu = fopen("u.m", "a");
fv = fopen("v.m", "a");
fuv = fopen("uv.m", "a");
fvv = fopen("vv.m", "a");
fp = fopen("p.m", "a");
iloop {
jloop {
fprintf(fu, " %16.14f", 0.5*(u[i][j]+u[i-1][j]));
fprintf(fv, " %16.14f", 0.5*(v[i][j]+v[i][j-1]));
fprintf(fp, " %16.14f", p[i][j]);
}
fprintf(fu, "\n");
fprintf(fv, "\n");
fprintf(fp, "\n");
}
i0loop {
jloop {
fprintf(fuv, " %16.14f", u[i][j]);
}
fprintf(fuv, "\n");
}
iloop {
j0loop {
fprintf(fvv, " %16.14f", v[i][j]);
}
fprintf(fvv, "\n");
}
fclose(fu);
fclose(fv);
fclose(fp);
fclose(fuv);
fclose(fvv);
return;
}
void print_data2(double **bd_posi)
{
extern int nt;
int k;
FILE *fxy;
fxy = fopen("bd_posi.m", "a");
for (k=1; k<=nt; k++) {
fprintf(fxy, " %16.14f %16.14f \n", bd_posi[k][1], bd_posi[k][2]);
fprintf(fxy, "\n");}
fclose(fxy);
return;
}
void pressure_update(double **a)
{
extern int nx, ny;
int i, j;
double ave = 0.0;
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++){
ave = ave + a[i][j];
}
ave /= (nx+0.0)*(ny+0.0);
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++){
a[i][j] -= ave;
}
return;
}
void augmenuv(double **u, double **v, int nx, int ny)
{
int i, j;
double max_val = 0.0;
for (j=1; j<=ny; j++) {
u[-1][j] = 2.0 - u[1][j]; //u[nx-1][j];
u[0][j] = 1.0;
u[nx+1][j] = u[nx][j]; //u[1][j];
}
for (i=-1; i<=nx+1; i++) {
u[i][0] = u[i][1];
u[i][ny+1] = u[i][ny];
}
for (j=0; j<=ny; j++) {
v[0][j] = -v[1][j]; //v[nx][j];
v[nx+1][j] = v[nx][j]; //v[1][j];
}
for (i=0; i<=nx+1; i++) {
v[i][0] = v[i][ny] = 0.0;
v[i][-1] = -v[i][1];
v[i][ny+1] = -v[i][ny-1];
}
}
This is an associated h file of C code:
#define gnx 512
#define gny 128
#define iloop for(i=1;i<=gnx;i++)
#define i0loop for(i=0;i<=gnx;i++)
#define jloop for(j=1;j<=gny;j++)
#define j0loop for(j=0;j<=gny;j++)
#define ijloop iloop jloop
#define i0jloop i0loop jloop
#define ij0loop iloop j0loop
#define iloopt for(i=1;i<=nxt;i++)
#define i0loopt for(i=0;i<=nxt;i++)
#define jloopt for(j=1;j<=nyt;j++)
#define j0loopt for(j=0;j<=nyt;j++)
#define ijloopt iloopt jloopt
#define i0jloopt i0loopt jloopt
#define ij0loopt iloopt j0loopt
void initialization(double **p, double **u, double **v, double **bd_posi);
void boundary_force(double **bd_posi, double **bd_force) ;
void fluid_force(double **bd_posi, double **bd_force);
void move_boundary(double **bd_posi, double **bd_velo, double **u, double **v);
double delta(double x);
void temp_uv(double **tu, double **tv, double **u, double **v, double **adv_u, double **adv_v);
void Poisson(double **tu, double **tv, double **p);
void MG_Poisson(double **p, double **f);
void source_uv(double **tu, double **tv, double **divuv, int nxt, int nyt);
void div_uv(double **tu, double **tv, double **divuv, int nxt, int nyt);
void grad_p(double **p, double **dpdx, double **dpdy, int nxt, int nyt);
void vcycle_uv(double **uf, double **ff, int nxf, int nyf, int ilevel);
void residual(double **r, double **u, double **f, int nxt, int nyt);
void laplace(double **p, double **lap_p, int nxt, int nyt);
void relax(double **p, double **f, int nxt, int nyt);
void restrict(double **u_fine, double **u_coarse, int nxt, int nyt);
void full_step(double **u, double **v, double **nu, double **nv, double **p);
void advection_step(double **u, double **v, double **adv_u, double **adv_v);
void prolong(double **u_coarse, double **u_fine, int nxt, int nyt);
double *dvector (long nl, long nh);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void mat_add(double **a, double **b, double **c,
int xl, int xr, int yl, int yr);
void mat_add2(double **a, double **b, double **c,
double **a2, double **b2, double **c2,
int xl, int xr, int yl, int yr);
void zero_vector(double *a, int xl, int xr);
void zero_matrix(double **a, int xl, int xr, int yl, int yr);
void mat_copy(double **a, double **b,
int xl, int xr, int yl, int yr);
void mat_copy2(double **a, double **b,
double **a2, double **b2,
int xl, int xr, int yl, int yr);
void mat_sub(double **a, double **b, double **c,
int nrl, int nrh, int ncl, int nch);
void mat_sub2(double **a, double **b, double **c,
double **a2, double **b2, double **c2,
int nrl, int nrh, int ncl, int nch);
double mat_max(double **a,
int nrl, int nrh, int ncl, int nch);
void print_data1(double **u, double **v, double **p);
void print_data2(double **bd_posi);
void pressure_update(double **a);
void augmenuv(double **u, double **v, int nx, int ny);
The MATLAB codes for posttreatment are as follows:
clear all;
xright=32.0; yright=8.0; xleft = 0; yleft = 0; nx=512; ny = 128;h= (xright-xleft)/nx;
x=linspace(xleft+0.5h,xright-0.5h,nx); y=linspace(yleft+0.5h,yright-0.5h,ny);
[xx,yy]=meshgrid(x,y);
ss=sprintf(‘…/bd_posi.m’);
C1=load(ss);
ss=sprintf(‘…/u.m’);
uuu=load(ss);
ss=sprintf(‘…/v.m’);
vvv=load(ss);
nt = 100;
for i= 1:11
clf;
xp1 = C1(1+nt*(1-1):nt*1,1);
yp1 = C1(1+nt*(1-1):nt*1,2);
plot(xp1,yp1,'k.-','linewidth',2,'markersize',2.0);hold on;
U = uuu((i-1)*nx+1:i*nx,:);
V = vvv((i-1)*nx+1:i*nx,:);
U = U';
V = V';
% Show velocity vector
kk=6; kk2 =6; h = quiver(xx(2:kk:end,2:kk2:end),yy(2:kk:end,2:kk2:end),...
0.6*U(2:kk:end,2:kk2:end),0.6*V(2:kk:end,2:kk2:end),0,'b');hold on
%Calculate vorticity
% wow(1:512,1:128) = 0;
% for i = 2:nx-1
% for j = 2:ny-1
% wow(i,j) = (V(i+1,j)-V(i-1,j))/(2h) - (U(i,j+1)-U(i,j-1))/(2h);
% if( sqrt((x(i)-7)^2 + (y(j)-4)^2) < 0.9)
% wow(i,j) = 0;
% end
% end
% end
%Show vorticity
% for mm = -2.7:0.4:2.7
% hold on;contour(xx,yy,wow’,[mm mm],’edgecolor’,’b’);hold on;
% end
for ii = 1:100
ttx(ii) = xp1(ii);
tty(ii) = yp1(ii);
end
patch(ttx,tty,’k’);hold on;
view(0,90);
axis image;
axis([0 32 0 8]);
set(gca,'fontsize',11);
pause(1);
end
This is an associated h file of C code:
The MATLAB codes for posttreatment are as follows: