# Getting a function/solution for an ODE over 2 dimensions

I am trying to calculate the laminar flow field within a tubular reactor. The relation between axial velocities, viscosity and pressure drop is given by this differential equation with the axial length x and the radius r:

For the pressure drop (d_p/d_x) I use a fixed value, thus transforming the equation into an ODE. The viscosity (eta) increases with the radius leading to an elongated velocity profile, as the solution of this ODE shows:

That is the code I used in Mathematica for this solution with a no slip boundary condition at the wall (vx==0) and the symmetry condition at the centre (d_vx/d_r==0):

eta0 = 1*10^-4
rmax = 0.0024
rmin = 1*10^-10
alpha = 5000
dpdx = -8
eta[r_] := eta0*E^(alpha*r)
lsg1d = NDSolve[{
dpdx == 1/r*D[r*eta[r]*D[vx[r], r], r],
DirichletCondition[vx[r] == 0, r == rmax],
Derivative[1][vx][rmin] == 0},
{vx}, {r, rmin, rmax}]
Plot[vx[r] /. lsg1d, {r, rmin, rmax}, AxesLabel -> Automatic]


So far so good, but for my actual problem the viscosity changes at every position x and the velocity profile has to be recalculated and becomes also a function of x. A naive attempt and extension of the previous code looks like this:

lsg2d = NDSolve[{
dpdx == 1/r*D[r*eta[r]*D[vx[x, r], r], r],
DirichletCondition[vx[x, r] == 0, r == rmax],
Derivative[0, 1][vx][x, rmin] == 0},
{vx}, {x, 0, 2}, {r, rmin, rmax}]


I would expect to get a surface along x and r with the same values as from the result above. Unfortunately this attempt is giving me errors in the boundary conditions. As there is no additional derivative present I think the boundary conditions should still be okay.

I would like to use the velocity profile later on in a system of partial differential equations describing convection-diffusion-reaction so having vx as a function of x and r would be great.

Any ideas to solve this problem are highly appreciated.

• You seem to try and solve a 2D Poisson equation. For this problem to be well-posed, you need boundary conditions along the entire boundary of your (now 2D) domain.
– Pirx
Commented Jul 18, 2016 at 17:49
• How does the viscosity change with x? Commented Jul 18, 2016 at 20:18
• The viscosity change is a quite complex function of the conversion. But an often used approximation is eta[x_,r_] := eta0*E^(alpha * r * x) Commented Jul 19, 2016 at 7:14
• Thanks for the answer, I am also thinking about additional boundary conditions but I think as long as I know the viscosity at a distinct position x the shape of vx is determined. So I should not need additional conditions. The error Mathematica is giving me seems to indicate that the Neumann condition (d_vx/d_r==0) is used as a Dirichlet condition which is of course nonsense. So I think it might be a also a problem of syntax. Commented Jul 19, 2016 at 7:19
• In principle I would like to span the result of the ODE over the axial dimension x in order to get a plane. For demonstration purposes it should be okay to keep the viscosity constant. Is there really no possibility? Commented Jul 19, 2016 at 13:27

I'm not sure what you're trying to achieve in your real problem, but according to your comment, in lsg2d, eta is a function of x and r and the function definition is
eta[x_,r_] := eta0*E^(alpha * r * x)

lsg2d = Table[NDSolve[{dpdx == 1/r*D[r*eta[x, r]*D[vx[r], r], r],