Practical tutorial and numerical implementation of multi-phase flow simulation based on phase-field model
Date:
This tutorial aims to introduce a simple and practical numerical framework for multi-phase incompressible and immiscible fluid simulation. The physical problem raised in this tutorial is N-component (N>=3) fluid mixing in a tilted channel. The mathematical model for capturing multi-phase fluid interfaces is the N-component conservative Allen-Cahn (CAC) equation, the fluid dynamics is governed by the incompressible Navier-Stokes equations.
The multi-component CAC equation is solved with a temporally first-order accurate semi-implicit operator splitting method. The incompressible Navier-Stokes equations with variable density and viscosity ratios are solved with a temporally first-order accurate pressure projection method. The governing equations are discretized in space with standard finite difference method. We describe the mathematical models and numerical implementations in details. The C (for computation) and Matlab (for post-treatment) codes of a four-component case are provided.
- The tutorial can be found in https://cfdyang521.github.io/files/paper5.pdf
/* The bnsch.c for computation is given as follows: */
#include
int nx, ny, n_level, p_relax, c_relax, count; double pi, xleft, xright, yleft, yright, h, h2, dt, gam, Cahn, rho1, rho2, rho3, rho4, vis1, vis2, vis3, vis4, Pe, Re, We, Fr, adv_u, **adv_v, **adv_c, **adv_c2, **adv_c3, **fx, **fy, **tu, **tv, **worku, **workv, **workp, **rho, **vis, **c4, **ct, **sc, **smu, **mu, title, Sum11, Sum12, Sum13, Sum21, Sum22, Sum23, SumF, **alp; char bufferu[100], bufferv[100],bufferp[100],bufferc[100], bufferc2[100], bufferc3[100], bufferc4[100]; int main() { extern int nx, ny, n_level, p_relax, c_relax, count; extern double pi, xleft, xright, yleft, yright, h, h2, dt, gam, Cahn, rho1, rho2, rho3, rho4, vis1, vis2, vis3, vis4, Pe, Re, We, Fr, **adv_u, **adv_v, **adv_c, **adv_c2, **adv_c3, **fx, **fy, **tu, **tv, **worku, **workv, **workp, **rho, **vis, **c4, **ct, **sc, **smu, **mu, title, Sum11, Sum12, Sum13, Sum21, Sum22, Sum23, SumF, **alp; extern char bufferu[100], bufferv[100],bufferp[100],bufferc[100], bufferc2[100], bufferc3[100], bufferc4[100]; int it, max_it, ns, i, j; double **u, **nu, **v, **nv, **oc, **nc, **oc2, **nc2, **oc3, **nc3,p;
FILE *fu, *fv, *fc, *fc2, *fc3, *fc4, *fp;
/*****/
nx = gnx;
ny = gny;
n_level = (int)(log(nx)/log(2)-0.9);
p_relax = 7;
c_relax = 5;
pi = 4.0*atan(1.0);
xleft = 0.0, xright = 8.0;
yleft = 0.0, yright = xright*(double)ny/(double)nx;
count = 0;
/***********************/
h = xright/(double)nx;
h2 = h*h;
dt = 0.02*h;
gam = 4.0*h/(4.0*sqrt(2.0)*atanh(0.9));
Cahn = pow(gam,2);
max_it = 16000;
ns = max_it/80;
rho1 = 4.0;
rho2 = 3.0;
rho3 = 2.0;
rho4 = 1.0;
vis1 = 1.0;
vis2 = 1.0;
vis3 = 1.0;
vis4 = 1.0;
Pe = 5.0/gam;
Re = 3000.0;
We = 50.0;
Fr = 1.0;
title = 60.0*pi/180.0;
/***********************/
printf("nx = %d, ny = %d\n", nx, ny);
printf("n_level = %d\n", n_level);
printf("dt = %f\n", dt);
printf("max_it = %d\n", max_it);
printf("ns = %d\n", ns);
printf("Reynolds number = %f\n", Re);
printf("Weber number = %f\n", We);
u = dmatrix(-1, nx+1, 0, ny+1);
nu = dmatrix(-1, nx+1, 0, ny+1);
v = dmatrix(0, nx+1, -1, ny+1);
nv = dmatrix(0, nx+1, -1, ny+1);
oc = dmatrix(0, nx+1, 0, ny+1);
nc = dmatrix(0, nx+1, 0, ny+1);
oc2 = dmatrix(0, nx+1, 0, ny+1);
nc2 = dmatrix(0, nx+1, 0, ny+1);
oc3 = dmatrix(0, nx+1, 0, ny+1);
nc3 = dmatrix(0, nx+1, 0, ny+1);
c4 = dmatrix(0, nx+1, 0, ny+1);
p = dmatrix(1, nx, 1, ny);
alp = dmatrix(0, nx+1, 0, ny+1);
adv_u = dmatrix(0, nx, 1, ny);
adv_v = dmatrix(1, nx, 0, ny);
adv_c = dmatrix(1, nx, 1, ny);
adv_c2 = dmatrix(1, nx, 1, ny);
adv_c3 = dmatrix(1, nx, 1, ny);
fx = dmatrix(0, nx, 1, ny);
fy = dmatrix(1, nx, 0, ny);
tu = dmatrix(0, nx, 1, ny);
tv = dmatrix(1, nx, 0, ny);
worku = dmatrix(0, nx, 1, ny);
workv = dmatrix(1, nx, 0, ny);
workp = dmatrix(1, nx, 1, ny);
rho = dmatrix(0, nx+1, 0, ny+1);
vis = dmatrix(0, nx+1, 0, ny+1);
ct = dmatrix(1, nx, 1, ny);
sc = dmatrix(1, nx, 1, ny);
smu = dmatrix(1, nx, 1, ny);
mu = dmatrix(1, nx, 1, ny);
zero_matrix(mu, 1, nx, 1, ny);
sprintf(bufferu,"datau.m");
sprintf(bufferv,"datav.m");
sprintf(bufferp,"datap.m");
sprintf(bufferc,"datac.m");
sprintf(bufferc2,"datac2.m");
sprintf(bufferc3,"datac3.m");
sprintf(bufferc4,"datac4.m");
fu = fopen(bufferu,"w");
fv = fopen(bufferv,"w");
fp = fopen(bufferp,"w");
fc = fopen(bufferc,"w");
fc2 = fopen(bufferc2,"w");
fc3 = fopen(bufferc3,"w");
fc4 = fopen(bufferc4,"w");
fclose(fu);
fclose(fv);
fclose(fp);
fclose(fc);
fclose(fc2);
fclose(fc3);
fclose(fc4);
initialization(u, v, oc, oc2, oc3, c4, p);
ijloop{
Sum11 = Sum11 + oc[i][j];
Sum12 = Sum12 + oc2[i][j];
Sum13 = Sum13 + oc3[i][j];
}
mat_copy(nu, u, 0, nx, 1, ny);
mat_copy(nv, v, 1, nx, 0, ny);
mat_copy(nc, oc, 1, nx, 1, ny);
mat_copy(nc2, oc2, 1, nx, 1, ny);
mat_copy(nc3, oc3, 1, nx, 1, ny);
print_data(nu, nv, nc, nc2, nc3, c4, p);
for (it=1; it<=max_it; it++) {
cal_den(oc, oc2, oc3, c4, rho);
augmenc(rho, nx, ny);
cal_vis(oc, oc2, oc3, c4, vis);
augmenc(vis, nx, ny);
full_step(u, v, oc, oc2, oc3, nu, nv, nc, nc2, nc3, p);
ijloop{ c4[i][j] = 1.0 - nc[i][j] - nc2[i][j] - nc3[i][j]; }
mat_copy(u, nu, 0, nx, 1, ny);
mat_copy(v, nv, 1, nx, 0, ny);
mat_copy(oc, nc, 1, nx, 1, ny);
mat_copy(oc2, nc2, 1, nx, 1, ny);
mat_copy(oc3, nc3, 1, nx, 1, ny);
printf("it = %d\n", it);
if (it % ns == 0) {
count++;
print_data(nu, nv, nc, nc2, nc3, c4, p);
printf("print out counts %d\n", count);
}
}
return 0; }
void initialization(double **u, double **v, double **c, double **c2, double **c3, double **c4, double **p) { extern double xright, yright, h, gam;
int i, j;
double x, y;
ijloop {
x = ((double)i-0.5)*h;
y = ((double)j-0.5)*h;
c[i][j] = 0.5 + 0.5*tanh((x-6.0)/(sqrt(2.0)*gam));
c2[i][j] = (1.0-c[i][j])*( 0.5 + 0.5*tanh((x-4.0)/(sqrt(2.0)*gam)) );
c3[i][j] = (1.0-c[i][j]-c2[i][j])*( 0.5 + 0.5*tanh((x-2.0)/(sqrt(2.0)*gam)) );
c4[i][j] = 1.0 - c[i][j] - c2[i][j] - c3[i][j];
p[i][j] = 0.0;
}
zero_matrix(u, 0, nx, 1, ny);
zero_matrix(v, 1, nx, 0, ny);
}
void cal_den(double **c, double **c2, double **c3, double **c4, double **rho) { extern double rho1, rho2, rho3, rho4;
int i, j;
ijloop {
rho[i][j] = ( rho1*c[i][j] + rho2*c2[i][j] + rho3*c3[i][j] + rho4*c4[i][j] )/rho1;
}
}
void cal_vis(double **c, double **c2, double **c3, double **c4, double **vis) { extern double vis1, vis2, vis3, vis4;
int i, j;
ijloop {
vis[i][j] = ( vis1*c[i][j] + vis2*c2[i][j] + vis3*c3[i][j] + vis4*c4[i][j] )/vis1;
}
}
void full_step(double **u, double **v, double **c, double **c2, double **c3, double **nu, double **nv, double **nc, double **nc2, double **nc3, double **p) { extern int nx, ny; extern double Sum11, Sum12, Sum13, dt, h2, h, **alp, **adv_u, **adv_v, **adv_c, **adv_c2, **adv_c3, **fx, **fy, **tu, **tv, **worku, **workv, **rho, title, Cahn; double x,y, beta1, beta2, beta3;
int i, j;
advection_step(u, v, c, c2, c3, adv_u, adv_v, adv_c, adv_c2, adv_c3);
// sf_force(c, fx, fy);
i0jloop {
fx[i][j] = -sin(title)/Fr;
}
ij0loop {
fy[i][j] = -cos(title)/Fr;
}
temp_uv(tu, tv, u, v, adv_u, adv_v, fx, fy);
Poisson(tu, tv, p);
grad_p(p, worku, workv, nx, ny);
i0jloop {
nu[i][j] = tu[i][j] - dt*worku[i][j]/(0.5*(rho[i+1][j]+rho[i][j]));
}
ij0loop {
nv[i][j] = tv[i][j] - dt*workv[i][j]/(0.5*(rho[i][j+1]+rho[i][j]));
}
iloop {
nv[i][0] = nv[i][ny] = 0.0;
}
//update CAC equations for three components
ijloop{
alp[i][j] = c[i][j](c[i][j]-0.5)(c[i][j]-1.0) + c2[i][j](c2[i][j]-0.5)(c2[i][j]-1.0)
+ c3[i][j](c3[i][j]-0.5)(c3[i][j]-1.0) + (1.0-c[i][j]-c2[i][j]-c3[i][j])( (1.0-c[i][j]-c2[i][j]-c3[i][j])-0.5)( (1.0-c[i][j]-c2[i][j]-c3[i][j])-1.0) ;
}
augmenc(c, nx, ny); augmenc(c2, nx, ny); augmenc(c3, nx, ny);
// Step 1:
// Semi-implicit scheme heat2(c, adv_c, nc); heat2(c2, adv_c2, nc2); heat2(c3, adv_c3, nc3);
//Step 2:
ijloop{
nc[i][j] = 0.5 + (nc[i][j]-0.5)/sqrt( exp(-dt/(2.0*Cahn*Pe)) + pow(2.0*nc[i][j]-1.0,2)*(1.0 - exp(-dt/(2.0*Cahn*Pe) )) );
nc2[i][j] = 0.5 + (nc2[i][j]-0.5)/sqrt( exp(-dt/(2.0*Cahn*Pe)) + pow(2.0*nc2[i][j]-1.0,2)*(1.0 - exp(-dt/(2.0*Cahn*Pe) )) );
nc3[i][j] = 0.5 + (nc3[i][j]-0.5)/sqrt( exp(-dt/(2.0*Cahn*Pe)) + pow(2.0*nc3[i][j]-1.0,2)*(1.0 - exp(-dt/(2.0*Cahn*Pe) )) );
}
// Step 3:
ijloop{
Sum21 = Sum21 + (c[i][j] - nc[i][j]);
Sum22 = Sum22 + (c2[i][j] - nc2[i][j]);
Sum23 = Sum23 + (c3[i][j] - nc3[i][j]);
SumF = SumF + sqrt( 0.25pow(nc[i][j],2)pow(1.0-nc[i][j],2) ) + sqrt( 0.25pow(nc2[i][j],2)pow(1.0-nc2[i][j],2) ) + sqrt( 0.25pow(nc3[i][j],2)pow(1.0-nc3[i][j],2) ) + sqrt( 0.25pow(1.0-nc[i][j]-nc2[i][j]-nc3[i][j],2)pow(1.0-(1.0-nc[i][j]-nc2[i][j]-nc3[i][j]),2) );
}
beta1 = Sum21/(dtSumF + 1.0e-16); beta2 = Sum22/(dtSumF + 1.0e-16); beta3 = Sum23/(dt*SumF + 1.0e-16);
ijloop{
c[i][j] = nc[i][j] + dtbeta1( sqrt( 0.25pow(nc[i][j],2)pow(1.0-nc[i][j],2) ) + sqrt( 0.25pow(nc2[i][j],2)pow(1.0-nc2[i][j],2) ) + sqrt( 0.25pow(nc3[i][j],2)pow(1.0-nc3[i][j],2) ) + sqrt( 0.25pow(1.0-nc[i][j]-nc2[i][j]-nc3[i][j],2)pow(1.0-(1.0-nc[i][j]-nc2[i][j]-nc3[i][j]),2) ) );
c2[i][j] = nc2[i][j] + dtbeta2( sqrt( 0.25pow(nc[i][j],2)pow(1.0-nc[i][j],2) ) + sqrt( 0.25pow(nc2[i][j],2)pow(1.0-nc2[i][j],2) ) + sqrt( 0.25pow(nc3[i][j],2)pow(1.0-nc3[i][j],2) ) + sqrt( 0.25pow(1.0-nc[i][j]-nc2[i][j]-nc3[i][j],2)pow(1.0-(1.0-nc[i][j]-nc2[i][j]-nc3[i][j]),2) ) );
c3[i][j] = nc3[i][j] + dtbeta3( sqrt( 0.25pow(nc[i][j],2)pow(1.0-nc[i][j],2) ) + sqrt( 0.25pow(nc2[i][j],2)pow(1.0-nc2[i][j],2) ) + sqrt( 0.25pow(nc3[i][j],2)pow(1.0-nc3[i][j],2) ) + sqrt( 0.25pow(1.0-nc[i][j]-nc2[i][j]-nc3[i][j],2)pow(1.0-(1.0-nc[i][j]-nc2[i][j]-nc3[i][j]),2) ) );
}
ijloop{
nc[i][j] = c[i][j];
nc2[i][j] = c2[i][j];
nc3[i][j] = c3[i][j];
}
}
void advection_step(double **u, double **v, double **c, double **c2, double **c3, double **adv_u, double **adv_v, double **adv_c, double **adv_c2, double **adv_c3) { extern int nx, ny; extern double h; double a, b, d, **ux, **uy, **vx, **vy;
int i, j, k;
ux = dmatrix(-1, nx, 1, ny);
uy = dmatrix(0, nx, 0, ny);
vx = dmatrix(0, nx, 0, ny);
vy = dmatrix(1, nx, -1, ny);
augmenuv(u, v);
augmenc(c, nx, ny);
augmenc(c2, nx, ny);
augmenc(c3, nx, ny);
ijloop { // conservative for phase field
adv_c[i][j] = (u[i][j]*(c[i+1][j]+c[i][j])-u[i-1][j]*(c[i][j]+c[i-1][j]))/(2.0*h)
+ (v[i][j]*(c[i][j+1]+c[i][j])-v[i][j-1]*(c[i][j]+c[i][j-1]))/(2.0*h);
adv_c2[i][j] = (u[i][j]*(c2[i+1][j]+c2[i][j])-u[i-1][j]*(c2[i][j]+c2[i-1][j]))/(2.0*h)
+ (v[i][j]*(c2[i][j+1]+c2[i][j])-v[i][j-1]*(c2[i][j]+c2[i][j-1]))/(2.0*h);
adv_c3[i][j] = (u[i][j]*(c3[i+1][j]+c3[i][j])-u[i-1][j]*(c3[i][j]+c3[i-1][j]))/(2.0*h)
+ (v[i][j]*(c3[i][j+1]+c3[i][j])-v[i][j-1]*(c3[i][j]+c3[i][j-1]))/(2.0*h);
}
//2nd_order ENO scheme
for (i=0; i<=nx-1; i++)
for (j=1; j<=ny; j++){
// u-part
if( (u[i][j] + u[i+1][j]) >= 0.0 ) k = i;
else k = i+1;
a = (u[k][j] - u[k-1][j])/h;
b = (u[k+1][j]-u[k][j])/h;
if( fabs(a) <= fabs(b) ) d = a;
else d = b;
ux[i][j] = u[k][j] + 0.5*h*d*( 1.0-2.0*(k-i) );
}
for (i=0; i<=nx; i++) for (j=1; j<=ny-1; j++){
if( (v[i][j] + v[i+1][j]) >= 0.0 ) k = j;
else k = j+1;
a = (u[i][k] - u[i][k-1])/h;
b = (u[i][k+1]-u[i][k])/h;
if( fabs(a) <= fabs(b) ) d = a;
else d = b;
uy[i][j] = u[i][k] + 0.5*h*d*( 1.0-2.0*(k-j) );
}
// v-part
for (i=1; i<=nx-1; i++) for (j=0; j<=ny; j++){
if( (u[i][j] + u[i][j+1]) >= 0.0 ) k = i;
else k = i+1;
a = (v[k][j] - v[k-1][j])/h;
b = (v[k+1][j]-v[k][j])/h;
if( fabs(a) <= fabs(b) ) d = a;
else d = b;
vx[i][j] = v[k][j] + 0.5*h*d*( 1.0-2.0*(k-i) );
}
for (i=1; i<=nx; i++)
for (j=0; j<=ny-1; j++){
if( (v[i][j] + v[i][j+1]) >= 0.0 ) k = j;
else k = j+1;
a = (v[i][k] - v[i][k-1])/h;
b = (v[i][k+1]-v[i][k])/h;
if( fabs(a) <= fabs(b) ) d = a;
else d = b;
vy[i][j] = v[i][k] + 0.5*h*d*( 1.0-2.0*(k-j) );
}
for (j=1; j<=ny; j++){ ux[-1][j] = ux[nx-1][j]; ux[nx][j] = ux[0][j]; }
for (i=0; i<=nx; i++){ uy[i][0] = 0.0; uy[i][ny] = 0.0; }
for (j=0; j<=ny; j++){ vx[0][j] = vx[1][j]; vx[nx][j] = vx[nx-1][j]; }
for (i=1; i<=nx; i++){ vy[i][-1] = 0.0; vy[i][ny] = 0.0; }
i0jloop { adv_u[i][j] = u[i][j](ux[i][j] - ux[i-1][j])/h + 0.25(v[i][j-1]+v[i+1][j-1]+v[i][j]+v[i+1][j])*(uy[i][j] - uy[i][j-1])/h; }
ij0loop{ adv_v[i][j] = 0.25(u[i-1][j]+u[i][j]+u[i-1][j+1]+u[i][j+1])(vx[i][j] - vx[i-1][j])/h + v[i][j]*(vy[i][j] - vy[i][j-1])/h; }
free_dmatrix(ux, -1, nx, 1, ny);
free_dmatrix(uy, 0, nx, 0, ny);
free_dmatrix(vx, 0, nx, 0, ny);
free_dmatrix(vy, 1, nx, -1, ny);
}
void augmenuv(double **u, double **v) { extern int nx, ny;
int i, j;
double aa;
for (j=1; j<=ny; j++) {
// aa = 0.5*(u[0][j]+u[nx][j]);
u[0][j] = u[nx][j] = 0.0;
u[-1][j] = -u[1][j];
u[nx+1][j] = u[nx-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+1][j] = -v[nx][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];
}
}
void augmenc(double **c, int nxt, int nyt) { int i, j;
for (j=1; j<=nyt; j++) {
c[0][j] = 2.0*c[1][j] - c[2][j];
c[nxt+1][j] = 2.0*c[nxt][j] - c[nxt-1][j];
}
for (i=0; i<=nxt; i++) {
c[i][0] = 2.0*c[i][1] - c[i][2];
c[i][nyt+1] = 2.0*c[i][nyt] - c[i][nyt-1];
}
}
void sf_force(double **c, double **fx, double **fy) { extern int nx, ny; extern double h, gam, We, Fr, **rho;
int i, j, zero_norm, m, n;
double fac, dx1[2][2], dy1[2][2], normd[2][2], adx, ady, dxx, dyy, **cx, **cy;
cx = dmatrix(0, nx+1, 0, ny+1);
cy = dmatrix(0, nx+1, 0, ny+1);
fac = 3.0*sqrt(2.0)*gam/(4.0*We);
augmenc(c, nx, ny);
ijloop {
dx1[0][0] = (c[i][j]+c[i][j-1]-c[i-1][j]-c[i-1][j-1])/(2.0*h);
dx1[0][1] = (c[i][j+1]+c[i][j]-c[i-1][j+1]-c[i-1][j])/(2.0*h);
dx1[1][0] = (c[i+1][j]+c[i+1][j-1]-c[i][j]-c[i][j-1])/(2.0*h);
dx1[1][1] = (c[i+1][j+1]+c[i+1][j]-c[i][j+1]-c[i][j])/(2.0*h);
dy1[0][0] = (c[i][j]-c[i][j-1]+c[i-1][j]-c[i-1][j-1])/(2.0*h);
dy1[0][1] = (c[i][j+1]-c[i][j]+c[i-1][j+1]-c[i-1][j])/(2.0*h);
dy1[1][0] = (c[i+1][j]-c[i+1][j-1]+c[i][j]-c[i][j-1])/(2.0*h);
dy1[1][1] = (c[i+1][j+1]-c[i+1][j]+c[i][j+1]-c[i][j])/(2.0*h);
zero_norm = 0;
for (m=0; m<2; m++) {
for (n=0; n<2; n++) {
normd[m][n] = sqrt(pow(dx1[m][n],2)+pow(dy1[m][n],2));
if (normd[m][n] <= 1.0e-9)
zero_norm = 1;
}
}
adx = (dx1[0][0]+dx1[0][1]+dx1[1][0]+dx1[1][1])/4.0;
ady = (dy1[0][0]+dy1[0][1]+dy1[1][0]+dy1[1][1])/4.0;
if (sqrt(pow(adx,2)+pow(ady,2)) <= 1.0e-9)
zero_norm = 1;
if ((!zero_norm) && (fabs(c[i][j]) < 0.98)) {
dxx = ((dx1[1][1]/normd[1][1]+dx1[1][0]/normd[1][0])-
(dx1[0][1]/normd[0][1]+dx1[0][0]/normd[0][0]))/(2.0*h);
dyy = ((dy1[1][1]/normd[1][1]+dy1[0][1]/normd[0][1])-
(dy1[1][0]/normd[1][0]+dy1[0][0]/normd[0][0]))/(2.0*h);
cx[i][j] = -fac*(dxx+dyy)*sqrt(pow(adx,2)+pow(ady,2))*adx;
cy[i][j] = -fac*(dxx+dyy)*sqrt(pow(adx,2)+pow(ady,2))*ady;
}
else
cx[i][j] = cy[i][j] = 0.0;
}
augmenc(cx, nx, ny);
augmenc(cy, nx, ny);
i0jloop {
fx[i][j] = 0.5*(cx[i+1][j]+cx[i][j])/(0.5*(rho[i+1][j]+rho[i][j]));
}
ij0loop {
fy[i][j] = 0.5*(cy[i][j+1]+cy[i][j])/(0.5*(rho[i][j+1]+rho[i][j]));
}
free_dmatrix(cx, 0, nx+1, 0, ny+1);
free_dmatrix(cy, 0, nx+1, 0, ny+1); }
void temp_uv(double **tu, double **tv, double **u, double **v, double **adv_u, double **adv_v, double **fx, double **fy) { extern double h, Re, **rho, **vis;
int i, j;
double vis_u, vis_d, vis_r, vis_l;
i0jloop {
vis_u = 0.25*(vis[i][j]+vis[i+1][j]+vis[i][j+1]+vis[i+1][j+1]);
vis_d = 0.25*(vis[i][j-1]+vis[i+1][j-1]+vis[i][j]+vis[i+1][j]);
tu[i][j] = u[i][j] + dt*( ( 2.0*(vis[i+1][j]*(u[i+1][j]-u[i][j])-vis[i][j]*(u[i][j]-u[i-1][j]))
+ vis_u*(u[i][j+1]-u[i][j]) - vis_d*(u[i][j]-u[i][j-1])
+ vis_u*(v[i+1][j]-v[i][j]) - vis_d*(v[i+1][j-1]-v[i][j-1]) )/(h*h*(0.5*(rho[i+1][j]+rho[i][j]))*Re)
- adv_u[i][j] + fx[i][j] );
}
ij0loop {
vis_r = 0.25*(vis[i][j]+vis[i+1][j]+vis[i][j+1]+vis[i+1][j+1]);
vis_l = 0.25*(vis[i-1][j]+vis[i][j]+vis[i-1][j+1]+vis[i][j+1]);
tv[i][j] = v[i][j] + dt*( ( 2.0*vis[i][j+1]*(v[i][j+1]-v[i][j])-2.0*vis[i][j]*(v[i][j]-v[i][j-1])
+ vis_r*(v[i+1][j]-v[i][j]) - vis_l*(v[i][j]-v[i-1][j])
+ vis_r*(u[i][j+1]-u[i][j]) - vis_l*(u[i-1][j+1]-u[i-1][j]) )/(h*h*(0.5*(rho[i][j+1]+rho[i][j]))*Re)
- adv_v[i][j] + fy[i][j]);
}
}
void Poisson(double **tu, double **tv, double **p) { extern int nx, ny; extern double **workp;
source_uv(tu, tv, workp, nx, ny);
MG_Poisson(p, workp);
}
void source_uv(double **tu, double **tv, double **divuv, int nxt, int nyt) { extern double dt;
int i, j;
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 MG_Poisson(double **p, double **f) { extern int nx, ny; extern double **rho;
int it_mg = 1, max_it = 100;
double resid = 1.0, resid2 = 10.0, tol = 1.0e-5, **sor;
sor = dmatrix(1, nx, 1, ny);
mat_copy(sor, p, 1, nx, 1, ny);
while (it_mg <= max_it && resid >= tol) {
vcycle_uv(p, f, rho, nx, ny, 1);
pressure_update(p);
mat_sub(sor, sor, p, 1, nx, 1, ny);
resid = mat_max(sor, 1, nx, 1, ny);
mat_copy(sor, p, 1, nx, 1, ny);
if (resid > resid2)
it_mg = max_it;
else
resid2 = resid;
it_mg++;
}
printf("Mac pressure iteration = %d residual = %16.14f\n", it_mg-1, resid);
free_dmatrix(sor, 1, nx, 1, ny);
return; }
void vcycle_uv(double **uf, double **ff, double **wf, int nxf, int nyf, int ilevel) { extern int n_level;
relax_p(uf, ff, wf, ilevel, nxf, nyf);
if (ilevel < n_level) {
int nxc, nyc;
double **rf, **fc, **uc, **wc;
nxc = nxf/2, nyc = nyf/2;
rf = dmatrix(1, nxf, 1, nyf);
fc = dmatrix(1, nxc, 1, nyc);
uc = dmatrix(1, nxc, 1, nyc);
wc = dmatrix(0, nxc+1, 0, nyc+1);
residual_den(rf, uf, ff, wf, nxf, nyf);
restrict1(rf, fc, nxc, nyc);
restrict1(wf, wc, nxc, nyc);
augmenc(wc, nxc, nyc);
zero_matrix(uc, 1, nxc, 1, nyc);
vcycle_uv(uc, fc, wc, nxc, nyc, ilevel+1);
prolong(uc, rf, nxc, nyc);
mat_add(uf, uf, rf, 1, nxf, 1, nyf);
relax_p(uf, ff, wf, ilevel, nxf, nyf);
free_dmatrix(rf, 1, nxf, 1, nyf);
free_dmatrix(fc, 1, nxc, 1, nyc);
free_dmatrix(uc, 1, nxc, 1, nyc);
free_dmatrix(wc, 0, nxc+1, 0, nyc+1);
}
}
void relax_p(double **p, double **f, double **w, int ilevel, int nxt, int nyt) { extern int ny, p_relax; extern double xright, Fr, title;
int i, j, iter;
double ht, ht2, a[4], sorc, coef;
ht = xright/(double)nxt;
ht2 = pow(ht,2);
for (iter=1; iter<=p_relax; iter++) {
ijloopt {
a[0] = 2.0/(w[i+1][j]+w[i][j]);
a[1] = 2.0/(w[i][j]+w[i-1][j]);
a[2] = 2.0/(w[i][j+1]+w[i][j]);
a[3] = 2.0/(w[i][j]+w[i][j-1]);
sorc = f[i][j];
coef = 0.0;
if (i==1) {
if(nxt == nx)
sorc -= (a[0]*p[i+1][j])/ht2 + sin(title)/(ht*Fr);
else
sorc -= (a[0]*p[i+1][j])/ht2;
coef -= a[0]/ht2;
}
else if (i==nxt) {
if(nxt == nx)
sorc -= -sin(title)/(ht*Fr) + (a[1]*p[i-1][j])/ht2;
else
sorc -= (a[1]*p[i-1][j])/ht2;
coef -= a[1]/ht2;
}
else {
sorc -= (a[0]*p[i+1][j] + a[1]*p[i-1][j])/ht2;
coef -= (a[0] + a[1])/ht2;
}
if (j==1) {
if (nyt==ny)
sorc -= a[2]*p[i][j+1]/ht2 + cos(title)/(ht*Fr);
else
sorc -= a[2]*p[i][j+1]/ht2;
coef -= a[2]/ht2;
}
else if (j==nyt) {
if (nyt==ny)
sorc -= -cos(title)/(ht*Fr) + a[3]*p[i][j-1]/ht2;
else
sorc -= a[3]*p[i][j-1]/ht2;
coef -= a[3]/ht2;
}
else {
sorc -= (a[2]*p[i][j+1] + a[3]*p[i][j-1])/ht2;
coef -= (a[2] + a[3])/ht2;
}
p[i][j] = sorc/coef;
}
}
}
void residual_den(double **r, double **u, double **f, double **den, int nxt, int nyt) { int i, j; double **dpdx, **dpdy;
dpdx = dmatrix(0, nxt, 1, nyt);
dpdy = dmatrix(1, nxt, 0, nyt);
grad_p(u, dpdx, dpdy, nxt, nyt);
i0jloopt {
dpdx[i][j] = dpdx[i][j]/(0.5*(den[i+1][j]+den[i][j]));
}
ij0loopt {
dpdy[i][j] = dpdy[i][j]/(0.5*(den[i][j+1]+den[i][j]));
}
div_uv(dpdx, dpdy, r, nxt, nyt);
mat_sub(r, f, r, 1, nxt, 1, nyt);
free_dmatrix(dpdx, 0, nxt, 1, nyt);
free_dmatrix(dpdy, 1, nxt, 0, nyt); }
void grad_p(double **p, double **dpdx, double **dpdy, int nxt, int nyt) { extern double xright, Fr, **rho, title;
int i, j;
double ht;
ht = xright/(double)nxt;
i0jloopt {
if (i==0)
if(nxt == nx)
dpdx[0][j] = -0.5*(rho[0][j]+rho[1][j])*sin(title)/Fr;
else
dpdx[0][j] = 0.0;
else if (i==nxt)
if(nxt == nx)
dpdx[nxt][j] = -0.5*(rho[nxt+1][j]+rho[nxt][j])*sin(title)/Fr;
else
dpdx[nxt][j] = 0.0;
else
dpdx[i][j] = (p[i+1][j] - p[i][j])/ht;
}
ij0loopt {
if (j==0) {
if (nyt==ny)
dpdy[i][0] = -0.5*(rho[i][1]+rho[i][0])*cos(title)/Fr;
else
dpdy[i][0] = 0.0;
}
else if (j==nyt) {
if (nyt==ny)
dpdy[i][nyt] = -0.5*(rho[i][nyt+1]+rho[i][nyt])*cos(title)/Fr;
else
dpdy[i][nyt] = 0.0;
}
else
dpdy[i][j] = (p[i][j+1] - p[i][j])/ht;
}
}
void restrict1(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 heat2(double **c_old, double **adv_c, double **c_new) { int max_it_CH = 50,it_mg = 1; double tol = 1.0e-6,resid = 1.0;
mat_copy(ct,c_old,1,nx,1,ny);
source(c_old, adv_c, sc);
while (it_mg <= max_it_CH && resid > tol){
vcycle(c_new, sc, nx ,ny,1);
resid = error(ct,c_new,nx,ny);
mat_copy(ct,c_new,1,nx,1,ny);
it_mg++;
}
printf("error %12.10f %d \n",resid,it_mg-1); }
double error(double **c_old, double **c_new, int nxt, int nyt) { double **r, res; r = dmatrix(1, nxt, 1, nyt); mat_sub(r, c_new, c_old, 1, nxt, 1, nyt); res = mat_max(r, 1, nxt, 1, nyt); free_dmatrix(r, 1, nxt, 1, nyt); return res; }
void source(double **c_old, double **adv_c, double **src_c) { extern int nx, ny; extern double dt, Cahn, Pe, **alp;
int i, j;
ijloop{
src_c[i][j] = c_old[i][j]/dt - adv_c[i][j] + alp[i][j]*c_old[i][j]/(Cahn*Pe);
}
}
void laplace_cac(double **a, double **lap_a, int nxt, int nyt) { extern double xright;
int i, j;
double ht2, dadx_L, dadx_R, dady_B, dady_T;
ht2 = pow(xright/(double)nxt,2);
ijloopt {
if (i > 1)
dadx_L = a[i][j] - a[i-1][j];
else
dadx_L = a[i+1][j] - a[i][j];
if (i < nxt)
dadx_R = a[i+1][j] - a[i][j];
else
dadx_R = a[i][j] - a[i-1][j];
if (j > 1)
dady_B = a[i][j] - a[i][j-1];
else
dady_B = a[i][j+1] - a[i][j];
if (j < nyt)
dady_T = a[i][j+1] - a[i][j];
else
dady_T = a[i][j] - a[i][j-1];
lap_a[i][j] = (dadx_R - dadx_L + dady_T - dady_B)/ht2;
}
}
void vcycle(double **cf_new, double **sc, int nxf,int nyf,int ilevel) { relax(cf_new, sc, 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, cf_new, sc, nxf, nyf);
restrict1(rf, fc, nxc, nyc);
zero_matrix(uc, 1, nxc, 1, nyc);
vcycle(uc, fc, nxc, nyc, ilevel + 1);
prolong(uc, rf, nxc, nyc);
mat_add(cf_new, cf_new, rf, 1, nxf, 1, nyf);
relax(cf_new, sc, nxf, nyf);
free_dmatrix(rf, 1, nxf, 1, nyf);
free_dmatrix(uc, 1, nxc, 1, nyc);
free_dmatrix(fc, 1, nxc, 1, nyc);
} }
void residual(double **r, double **u, double **f, int nxt, int nyt) { extern double dt, Pe; int i, j;
laplace_cac(u, r, nxt, nyt); ijloopt{ r[i][j] = f[i][j] - u[i][j]/dt + r[i][j]/Pe; } }
void relax(double **p,double **f,int nxt,int nyt) { extern double dt, Pe, Cahn;
int i,j,iter;
double a, ht2, coef, src;
ht2 = pow((xright-xleft) / (double) nxt,2);
a = 1.0/Pe;
for (iter=1; iter<=2; iter++) {
ijloopt{
src = f[i][j];
coef = 1.0/dt;
if (i > 1 && i < nxt) src += a*(p[i-1][j]+p[i+1][j])/ht2;
if (j > 1 && j < nyt) src += a*(p[i][j-1]+p[i][j+1])/ht2;
if (i>1 && i<nxt) coef += 2.0*a/ht2;
if (j>1 && j<nyt) coef += 2.0*a/ht2;
p[i][j] = src/coef;
}
} }
/***** util ******/ double **dmatrix(long nrl, long nrh, long ncl, long nch) { double **m; long i, nrow=nrh-nrl+1+NR_END, ncol=nch-ncl+1+NR_END;
m=(double **) malloc((nrow)*sizeof(double*));
m+=NR_END;
m-=nrl;
m[nrl]=(double *) malloc((nrow*ncol)*sizeof(double));
m[nrl]+=NR_END;
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-NR_END); free(m+nrl-NR_END);
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_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 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_mat(FILE *fptr, double **a, int nrl, int nrh, int ncl, int nch) { int i, j;
for (i=nrl; i<=nrh; i++) {
for (j=ncl; j<=nch; j++)
fprintf(fptr, " %16.14f", a[i][j]);
fprintf(fptr, "\n");
}
return; }
void print_data(double **u, double **v, double **c, double **c2, double **c3, double **c4, double **p) { extern char bufferu[100], bufferv[100],bufferp[100],bufferc[100], bufferc2[100], bufferc3[100], bufferc4[100]; int i, j; FILE *fu, *fv, *fp, *fc, *fc2, *fc3, *fc4; fu = fopen(bufferu,”a”); fv = fopen(bufferv,”a”); fp = fopen(bufferp,”a”); fc = fopen(bufferc,”a”); fc2 = fopen(bufferc2,”a”); fc3 = fopen(bufferc3,”a”); fc4 = fopen(bufferc4,”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(fc, " %16.14f", c[i][j]);
fprintf(fc2, " %16.14f", c2[i][j]);
fprintf(fc3, " %16.14f", c3[i][j]);
fprintf(fc4, " %16.14f", c4[i][j]);
}
fprintf(fv, "\n");
fprintf(fp, "\n");
fprintf(fc, "\n");
fprintf(fc2, "\n");
fprintf(fc3, "\n");
fprintf(fc4, "\n");
}
fclose(fu);
fclose(fv);
fclose(fp);
fclose(fc);
fclose(fc2);
fclose(fc3);
fclose(fc4);
return; }
void pressure_update(double **a) { extern int nx, ny;
int i, j;
double ave = 0.0;
ijloop {
ave = ave + a[i][j];
}
ave /= (nx+0.0)*(ny+0.0);
ijloop {
a[i][j] -= ave;
}
return; }
/* The bnsch.h is given as follows: */
void initialization(double **u, double **v, double **c, double **c2, double **c3, double **c4, double **p);
void cal_den(double **c, double **c2, double **c3, double **c4, double **rho);
void cal_vis(double **c, double **c2, double **c3, double **c4, double **vis);
void full_step(double **u, double **v, double **c, double **c2, double **c3, double **nu, double **nv, double **nc, double **nc2, double **nc3, double **p);
void advection_step(double **u, double **v, double **c, double **c2, double **c3, double **adv_u, double **adv_v, double **adv_c, double **adv_c2, double **adv_c3);
void augmenuv(double **u, double **v);
void augmenc(double **c, int nxt, int nyt);
void sf_force(double **c, double **fx, double **fy);
void temp_uv(double **tu, double **tv, double **u, double **v, double **adv_u, double **adv_v, double **fx, double **fy);
void Poisson(double **tu, double **tv, double **p);
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 MG_Poisson(double **p, double **f);
void vcycle_uv(double **uf, double **ff, double **wf, int nxf, int nyf, int ilevel);
void relax_p(double **p, double **f, double **w, int ilevel, int nxt, int nyt);
void residual_den(double **r, double **u, double **f, double **den, int nxt, int nyt);
void grad_p(double **p, double **dpdx, double **dpdy, int nxt, int nyt);
void grad_p2(double **p, double **dpdx, double **dpdy, int nxt, int nyt);
void restrict1(double **u_fine, double **u_coarse, int nxt, int nyt);
void prolong(double **u_coarse, double **u_fine, int nxt, int nyt);
void heat2(double **c_old, double **adv_c, double **c_new);
double error(double **c_old, double **c_new, int nxt, int nyt);
void source(double **c_old, double **adv_c, double **src_c);
void laplace_cac(double **a, double **lap_a, int nxt, int nyt);
void vcycle(double **cf_new, double **sc, int nxf,int nyf,int ilevel);
void residual(double **r, double **u, double **f, int nxt, int nyt);
void relax(double **p,double **f,int nxt,int nyt);
/* The util.h is given as follows: */
#define gnx 512 #define gny 64
#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
double **dmatrix(long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
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_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 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_mat(FILE *fptr, double **a, int nrl, int nrh, int ncl, int nch);
void print_data(double **u, double **v, double **c, double **c2, double **c3, double **c4, double **p);
void pressure_update(double **a);
/* The Matlab code for post-treatment is given as follows */
clear all; xright=8.0; yright=1.0; xleft = 0; yleft = 0; nx=2562; ny = 322;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(‘…/datac.m’); phi = load(ss); ss=sprintf(‘…/datac2.m’); phi2 = load(ss); ss=sprintf(‘…/datac3.m’); phi3 = load(ss); ss=sprintf(‘…/datac4.m’); phi4 = load(ss); n=size(phi,1)/nx;
ii=0; for i = 15
clf A = phi((i-1)nx+1:inx,:); B = phi2((i-1)nx+1:inx,:); C = phi3((i-1)nx+1:inx,:); D = phi4((i-1)nx+1:inx,:);
for i = 1:nx for j = 1:ny
if( A(i,j) > 1.0)
A(i,j) = 1.0;
elseif( A(i,j) < 0.0)
A(i,j) = 0.0;
else
A(i,j) = A(i,j);
end
if( B(i,j) > 1.0)
B(i,j) = 1.0;
elseif( A(i,j) < 0.0)
B(i,j) = 0.0;
else
B(i,j) = B(i,j);
end
if( C(i,j) > 1.0)
C(i,j) = 1.0;
elseif( A(i,j) < 0.0)
C(i,j) = 0.0;
else
C(i,j) = C(i,j);
end
if( D(i,j) > 1.0)
D(i,j) = 1.0;
elseif( A(i,j) < 0.0)
D(i,j) = 0.0;
else
D(i,j) = D(i,j);
end
end end
hh=contourf(xx,yy,A’,[0.5 0.5],’facecolor’,’r’,’edgecolor’,[0 0 0]); hold on; hh2=contourf(xx,yy,B’,[0.5 0.5],’facecolor’,[0.3 0.75 0.93],’edgecolor’,[0 0 0]); hold on; hh3=contourf(xx,yy,C’,[0.5 0.5],’facecolor’,’g’,’edgecolor’,[0 0 0]); hold on; hh4=contourf(xx,yy,D’,[0.5 0.5],’facecolor’,’none’,’edgecolor’,[0 0 0]); hold on;
axis image; axis image; axis([0 8 0 1]); %axis off box on; camroll(60) pause(0.6);
end