ME19a.			Mathematica Code			Nov. 24, 2009.  Due Dec. 4

This notebook is used to solve the potential flow problem given as me19a-hw8:

PROBLEM 30 (Value 4 times the regular homework problems.)

Obtain a numerical solution to the planar, potential flow around a bend using the grid sketched below.

The upstream and downstream boundaries to be used are ϕ=0 and ϕ=1 respectively.  You are not required use any special procedures to treat the singular behaviour near the projecting corner.  You are advised to use an iterative method with over-relaxation and, by trial and error, to find an effective value of the over-relaxation factor which improves the convergence of the iterative process. Find:

(a) The velocity potential at all of the nodes and, in particular, the node A.

(b) If the volume averaged velocity, U, at inlet (B C) is unity find the velocity distributions along both the interior and exterior walls.  It is best to evaluate the velocity at points on the walls midway between the nodes.

(c) If the pressure coefficient is defined as (p-(Subscript[p, B]))/0.5ρ (U^2) where p is the pressure, (Subscript[p, B]) is the pressure at the point B and ρ is the fluid density, find the pressure coefficient along both the interior and exterior walls.  It is best to evaluate the pressures at the same points at which the velocities were evaluated.

The solution method is outlined in the file me19a-hw8-solutions.pdf and is implemented below.

(a) Find the velocity potential at all of the nodes and, in particular, the node A.

A sensitivity analysis will reveal that a relaxation factor of alpha=1.825 will lead to fastest convergence, without regard for the initialized value of the computational domain.  This value is used below and the computational domain is initialized to 0.6.

(* initialize domain *)
phi=Table[0.6,{15},{15}];

(* set inlet condition at nodes, i=15,j=1-5 *)
phi[[15,1;;5]]=0;

(* set outlet condition at nodes, i=1-9,j=15 *)
phi[[1;;9,15]]=1;

(* relaxation value *)
alpha=1.825;

(* counter to track number of iterations *)
count=0;

(* allowable error to check convergence *)
tolerance=10^-6;

(* iteration loop *)
phiNew=phi;While[True,{

totErrorSq=0;

For[j=14,j>0,j--,
For[i=14,i>0,i--,{

(* solve only for i,j in domain of interest *)
If[!(i>9&&j>5),{

If[i==1&&j==1(* top left corner *),
phiNew[[i,j]]=0.25*(2*phi[[i+1,j]]+2*phi[[i,j+1]]),

(* else1 *)If[j==1(* top wall *),
phiNew[[i,j]]=0.25*(phi[[i+1,j]]+phi[[i-1,j]]+2*phi[[i,j+1]]),

(* else2 *)If[i==1(* left wall *),
phiNew[[i,j]]=0.25*(phi[[i,j+1]]+phi[[i,j-1]]+2*phi[[i+1,j]]),

(* else3 *)If[i==9&&j>5(* right wall *),
phiNew[[i,j]]=0.25*(phi[[i,j+1]]+phi[[i,j-1]]+2*phi[[i-1,j]]),

(* else4 *)If[i>9&&j==5(* bottom wall *),
phiNew[[i,j]]=0.25*(phi[[i+1,j]]+phi[[i-1,j]]+2*phi[[i,j-1]]),

(* else the rest are interior nodes *)
phiNew[[i,j]]=0.25*(phi[[i+1,j]]+phi[[i-1,j]]+phi[[i,j+1]]+phi[[i,j-1]])

](* end4 *)
](* end3 *)
](* end2 *)
](* end1 *)
](* end *);

(* use relaxation constant to calculate new phi *)
phiNew[[i,j]]=(1-alpha)*phi[[i,j]]+alpha*phiNew[[i,j]];

(* sum the errors *)
errorSq=(phiNew[[i,j]]-phi[[i,j]])^2;
totErrorSq+=errorSq;

(* update solution for phi *)
phi[[i,j]]=phiNew[[i,j]];

}](* end if *);

}](* end for i *)
](* end for j *);

count++;

(* convergence check *)
If[Sqrt[totErrorSq]<tolerance,Break[]];

}];(* end while *)

Print["Domain initialized to ",phi[[15,15]]]
Print["Solution converged in ",count," iterations with relaxation constant ",alpha]
Print["Total error is ",Sqrt[totErrorSq]]
Print["The solution at node A is ",phi[[9,5]]]

solution=Table[phi[[i,j]],{j,15},{i,15}];
Style[TableForm[solution/.{phi[[15,15]]->""},TableHeadings->{Table[j,{j,15}],Table[i,{i,15}]}],6.5]

phi0=Table[0.6,{15},{15}];
phi0[[14,1;;5]]=0;
phi0[[1;;9,14]]=1;
pphi0=Show[
ListContourPlot[Flatten[Table[{i,-j,phi0[[i,j]]},{j,14},{i,14}],1],Contours->20,ContourStyle->None,ColorFunction->"Rainbow"],
Graphics[{White,Rectangle[{9,-15},{15,-5}]}],
FrameTicks->None,Frame->False,AspectRatio->1
];

pphi=Show[
ListContourPlot[Flatten[Table[{i,-j,phi[[i,j]]},{j,14},{i,14}],1],Contours->50,ColorFunction->"Rainbow"],
Graphics[{White,Rectangle[{9,-15},{15,-5}]}],
FrameTicks->None,Frame->False,AspectRatio->1
];

(b) Find the velocity distributions along both the interior and exterior walls.

Since velocities are calculated between nodes, u[[i,j]] corresponds to the velocity between nodes i and i+1, and v[[i,j]] corresponds to the velocity between nodes j+1 and j.

dx=dy=1;

(* initialize u *)
u=Table[0,{14},{14}];

(* initialize v *)
v=Table[0,{14},{14}];

For[j=14,j>0,j--,
For[i=14,i>0,i--,{

(* solve only for i,j in domain of interest *)
If[!(i>9&&j>5),{

u[[i,j]]=(phi[[i,j]]-phi[[i+1,j]])/dx;

v[[i,j]]=(phi[[i,j+1]]-phi[[i,j]])/dy

}](* end if *);

}](* end for i *)
](* end for j *);

UInlet=0.25*(0.5*u[[14,1]]+u[[14,2]]+u[[14,3]]+u[[14,4]]+0.5*u[[14,5]]);

u=u/UInlet;
u[[9,6;;14]]=0;

v=v/UInlet;
v[[10;;14,5]]=0;

usolution=Table[u[[i,j]],{j,14},{i,14}];
Style[TableForm[usolution/.{0->""},TableHeadings->{Table[j,{j,15}],Table[ToString[i]<>"-"<>ToString[i+1],{i,14}]}],6.5]

Note that u[[i,15]]=0 for all i because phi[[i,15]]=phi[[i+1,15]]=1 at the outlet.  Positive u is defined to be to the left.

vsolution=Table[v[[i,j]],{j,14},{i,14}];
Style[TableForm[vsolution/.{0->""},TableHeadings->{Table[ToString[j]<>"-"<>ToString[j+1],{j,14}],Table[i,{i,14}]}],6.5]

Note that v[[15,j]]=0 for all j because phi[[15,j]]=phi[[15,j+1]]=0 at the inlet.  Positive v direction is down.

pu=Show[
ListContourPlot[Flatten[Table[{i,-j,u[[i,j]]},{j,14},{i,14}],1],Contours->100,ContourStyle->None,ColorFunction->"Rainbow"],
Graphics[{White,Rectangle[{9,-14},{14,-5}]}],
FrameTicks->None,Frame->False,AspectRatio->1
];

pv=Show[
ListContourPlot[Flatten[Table[{i,-j,v[[i,j]]},{j,14},{i,14}],1],Contours->100,ContourStyle->None,ColorFunction->"Rainbow"],
Graphics[{White,Rectangle[{9,-14},{14,-5}]}],
FrameTicks->None,Frame->False,AspectRatio->1
];

pvel=Show[
ListContourPlot[Flatten[Table[{i,-j,Sqrt[u[[i,j]]^2+v[[i,j]]^2]},{j,14},{i,14}],1],Contours->100,ContourStyle->None,ColorFunction->"Rainbow"],
Graphics[{White,Rectangle[{9,-14},{14,-5}]}],
FrameTicks->None,Frame->False,AspectRatio->1
];

(c) Find the pressure coefficient, (p-(Subscript[p, B]))/0.5`ρ (U^2), along both the interior and exterior walls.

Along a streamline from B we use Bernoulli's to obtain (Subscript[p, B])+1/2ρ ((Subscript[u, B])^2)=p+1/2ρ ((u^2)+(v^2)) or rearranging (Subscript[c, p])=(p-(Subscript[p, B]))/0.5ρ (U^2)=(((Subscript[u, B])^2)-((u^2)+(v^2)))/(U^2).  Note that we have already normalized our calculated velocities, so cp[[i,j]]=u[[14,1]]^2-(u[[i,j]]^2+v[[i,j]]^2).

(* initialize cp *)
cp=Table[0,{14},{14}];

For[j=14,j>0,j--,
For[i=14,i>0,i--,{

(* solve only for i,j in domain of interest *)
If[!(i>9&&j>5),{

cp[[i,j]]=(u[[14,1]]^2+v[[14,1]]^2)-(u[[i,j]]^2+v[[i,j]]^2)

}](* end if *);

}](* end for i *)
](* end for j *);

cpsolution=Table[cp[[i,j]],{j,14},{i,14}];
Style[TableForm[cpsolution/.{0->""},TableHeadings->{Table[ToString[j]<>"-"<>ToString[j+1],{j,14}],Table[ToString[i]<>"-"<>ToString[i+1],{i,14}]}],6.5]

pcp=Show[
ListContourPlot[Flatten[Table[{i,-j,cp[[i,j]]},{j,14},{i,14}],1],Contours->100,ContourStyle->None,ColorFunction->"Rainbow"],
Graphics[{White,Rectangle[{9,-14},{14,-5}]}],
FrameTicks->None,Frame->False,AspectRatio->1
];

Colorful contour plots

GraphicsGrid[{{pphi0,pphi,pcp},{pu,pv,pvel}},Spacings->0]

The top row shows (left) the initialized computational domain with the given boundary conditions, (center) the solution for the velocity potential where spacing of the contour lines indicates the relative flow velocity, and (right) the pressure coefficient --- there is a large pressure drop around the corner at node A.  The bottom row of images shows (left) the u velocity component, (center) the v velocity component, and (right) the speed or magnitude of the velocity.