You can regard the trz
separetly.
First I rationalize all parameters and rewrite equation for simple testing later.
ode = {0 == -τrz[
r] + η/(1 + (ϵ λ/η) τzz) u'[r],
1/r D[r τrz[r], r] == dpdz, u[Ri] == 0, u[Ro] == 0}
Since trz
does not depend on u
, treat it separatly.
dsol = DSolve[τrz[r]/r + Derivative[1][τrz][r] ==
1000, τrz, r]
{* {{τrz -> Function[{r}, 500 r + C[1]/r]} *}
Insert the result and apply a starting condition instead of boundary condition, since NDSolve
can not find an explicit equation for u'[r]
and therefore switches to differential-algbraic equation.
ode2 = Join[Most[ode] /. First@dsol /. C[1] -> c, {u'[Ri] == aa}]
Omitt the trz
part
ode3 = ode2[[{1, 3, 4}]]
Find a solution for the c
parameter in dependance of aa
with a little trick, since Solve
does not find a solution for undefined aa
.
cc2[aa_] =
c /. First@
Solve[ode3[[1]] /. r -> Ri /. u'[Ri] -> -Pi, c, Reals] /.
Pi -> -aa;
Since shooting method did not work with Version 8.0, I guessed a value for aa
. Find the exact aa
with FindRoot
.
ndsol = NDSolve[ode3 /. c -> cc2[aa] /. aa -> -88.1, u, {r, Ri, Ro}]
The solution satisfies ode.
u[Ri] /. ndsol
u[Ro] /. ndsol
Plot[u[r] /. ndsol, {r, Ri, Ro},
Epilog -> {Red, Point[{Ro, u[Ro] /. First@ndsol}]}]

and the error (would be lower with aa found by FindRoot
)
Plot[Evaluate[
ode3[[1, 2]] /. c -> cc2[aa] /. aa -> -88.1 /. First@ndsol], {r, Ri,
Ro}, PlotRange -> All]
