Solution of a Nonlinear Boundary Value Problem:
Analysis of Steady Inviscid Axisymmetric
Fluid Flow With Swirl

Mathematica in Education
Vol.2 No.4
(October 1993)

by
Rob ert Spall
Department of Mechanical Engineering
University of South Alabama
spall@cmb00.larc.nasa.gov
or
spall@me1.me.usouthal.edu

SpallFig.1.eps

Needs["Graphics`Arrow`"]

?Arrow

a1 = Arrow[{0,0}, {8,0}];
a2 = Arrow[{0,0}, {0,3.5}];
a3 = Arrow[{.6,0}, {.6,2.5}];
a4 = Arrow[{1.2,0}, {1.2,1.2}];
a5 = Arrow[{6.5,0}, {6.5,1.2}];
a6 = Arrow[{7.1,0}, {7.1,2.7}];
a7 = Arrow[{1.7,1.7}, {5.6,1.7}];
t1 = Text[FontForm["Flow direction", {"Helvetica", 8}],
		{3.65,1.9}];
t2 = Text[FontForm["r", {"Times-Italic", 7}],
		{.25,3.6}];
t3 = Text[FontForm["a1", {"Times-Italic", 7}],
		{.9,2.1}];
t4 = Text[FontForm["a2", {"Times-Italic", 7}],
		{1.4,.8}];
t5 = Text[FontForm["b2", {"Times-Italic", 7}],
		{6.2,.8}];
t6 = Text[FontForm["b1", {"Times-Italic", 7}],
		{7.5,2.4}];
t7 = Text[FontForm["x", {"Times-Italic", 7}],
		{7.9,.2}];
l1 = Line[{{.9,1.2},{6.8,1.2}}];
l2 = Line[{{.9,2.5},{3.2,2.5},{4.4,2.7},{6.8,2.7}}];
Show[Graphics[{a1,a2,a3,a4,a5,a6,a7,
	t1,t2,t3,t4,t5,t6,t7,
	l1,l2}]];

Wall geometry

a1 = 4.0;
b1 = 4.1;
a2 = 10.^-10;
b2 = 10.^-10;

Vortex Distribution

V0 = 1.5;
U = 1;
psi[r_] := 0.5 r^2 U;

Equations and Boundary Conditions

bc1 = f[b1] == 0.5 a1^2 U;
d = 0.5(a2^2 U);

eqn = f[r](f''[r] - 1/r f'[r] + 
         2 V0^2/U Exp[-2/U f[r]] * 
           (1.-Exp[-2/U f[r]])) + 
             V0^2 r^2( Exp[-4/U f[r]] - 
               Exp[-2/U f[r]]) == 0;

Nonlinear ODE Solver Code

errmax = 10^-8;
guess1 = b1 U;
guess2 = 1.05 b1 U;
bc2 = f'[b1] == guess1;

soln1 = NDSolve[{eqn, bc1, bc2}, f, {r, b1, b2}];

result1 = First[f[b2] /. soln1];
bc2 = f'[b1] == guess2;
soln2 = NDSolve[{eqn,bc1,bc2}, f, {r,b1,b2}];

result2 = First[f[b2] /. soln2];
iter = 0;
error = Abs[result2 - d]

While[error>errmax && iter < 10,
	guess = guess1 + (guess2 - guess1)/(result2-result1)(d-result1);
	bc2 = f'[b1] == guess;
	soln = NDSolve[{eqn,bc1,bc2}, f, {r,b1,b2}];
	guess1 = guess2;
	guess2 = guess;
	result1 = result2;
	result2 = First[f[b2] /. soln];
	error = Abs[result2 - d];
	iter = iter+1;
	Print["iter=", iter, " error=", error, " guess=", guess]]

Plotting

u[r_] := 1/r (First[f'[r] /. soln])

SpallFig.2.eps

V0 = 0.0;
plot1 = Plot[Evaluate[u[r] /. soln], {r,b2,b1}, 
			Frame -> True,
			GridLines -> Automatic,
			FrameLabel -> {"r", "U"},
			PlotRange -> {-1.25,1.25},
			DisplayFunction->Identity];

V0 = 0.49;
plot2 = Plot[Evaluate[u[r] /. soln], {r,b2,b1}, 
			Frame -> True,
			GridLines -> Automatic,
			FrameLabel -> {"r", "U"},
			PlotRange -> {-1.25,1.25},
			DisplayFunction->Identity];

V0 = 1.0;
plot3 = Plot[Evaluate[u[r] /. soln], {r,b2,b1}, 
			Frame -> True,
			GridLines -> Automatic,
			FrameLabel -> {"r", "U"},
			PlotRange -> {-1.25,1.25},
			DisplayFunction->Identity];

V0 = 1.2;
plot4 = Plot[Evaluate[u[r] /. soln], {r,b2,b1}, 
			Frame -> True,
			GridLines -> Automatic,
			FrameLabel -> {"r", "U"},
			PlotRange -> {-1.25,1.25},
			DisplayFunction->Identity];

V0 = 1.5;
plot5 = Plot[Evaluate[u[r] /. soln], {r,b2,b1}, 
			Frame -> True,
			GridLines -> Automatic,
			FrameLabel -> {"r", "U"},
			PlotRange -> {-1.25,1.25},
			DisplayFunction->Identity];

Show[plot1,plot2,plot3,plot4,plot5,
	DisplayFunction->$DisplayFunction];