3
$\begingroup$

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:

enter image description here

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:

velocity profile from ODE

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.

Improve this question
$\endgroup$
8
  • $\begingroup$ 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. $\endgroup$
    – Pirx
    Commented Jul 18, 2016 at 17:49
  • $\begingroup$ How does the viscosity change with x? $\endgroup$
    – Young
    Commented Jul 18, 2016 at 20:18
  • $\begingroup$ 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) $\endgroup$
    – Sebbo
    Commented Jul 19, 2016 at 7:14
  • $\begingroup$ 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. $\endgroup$
    – Sebbo
    Commented Jul 19, 2016 at 7:19
  • $\begingroup$ 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? $\endgroup$
    – Sebbo
    Commented Jul 19, 2016 at 13:27

1 Answer 1

Reset to default
1
$\begingroup$

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)

Then why not a simple 

lsg2d = Table[NDSolve[{dpdx == 1/r*D[r*eta[x, r]*D[vx[r], r], r], 
  vx[rmax] == 0, Derivative[1][vx][rmin] == 0}, {vx}, {r, rmin, rmax}], {x, 0, 2, 1/25}]
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for looking into this after some time passed. :-) As the question is formulated I think your approach is answering it. My actual problem is now a bit different and I am already thinking if "Table" can help there. I may come back soon with a new question and giving more details. :-) $\endgroup$
    – Sebbo
    Commented Oct 25, 2016 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.