Attempting to write a code that numerically solves for Poiseuille Flow, which is this the right methodology to take?

Question

I am looking at writing a C++ code to solve for a Poiseuille Flow between two parallel plates. The way I am thinking about doing this is to break up the second order DE into a system of 2 DE's and then solving the first DE by using Euler's method and then the second equation using RK4 or to run the DE twice through a solver (say Euler for example). So the system is for velocity is: $$ \frac{d^2 U_x}{dy^2} = \frac{1}{\mu} \frac{d P}{dx} $$ Where $\frac{d P}{dx}$ is a constant and $\mu$ is the viscosity. So since both of these values are essentially known and constant the system is really reduced to: $$ \frac{d^2 U_x}{dy^2} = C $$ where C is a constant.

Now there are two methods that I was looking at doing.

The first was to run the second order ode through a solver (say Euler for convince) get values, and then run it through the solver a second time.

Or, I was looking at breaking it into a system of first order ODE's but the system the occurs doesn't seem to help with the numerical process.

I know this has an analytical solution but what I am doing is for understanding and trying to see how different boundary conditions (i.e velocity at the boundaries, or different pressure gradients) affect the outcome. My main goal would be able to adapt the code to a general case so I can see the difference between Poiseuille flow, Couette Flow, etc. and see how the velocity profiles change. I also know there is a program like OpenFOAM that can do this but it seems overkill to me to use OpenFOAM for this simple of a problem, because it is essentially solving an ODE. Any help would be greatly appreciated.

Answer 1

The system you want to solve can be summarized as $$\frac{d^2 u}{dx^2} = c,~~~u(0) = u(L) = 0$$ which has the analytical solution $u(x) = \frac{c}{2}x(x-L)$.

This is a boundary value problem apposed to the more common initial value problem where all the conditions are specified at a single point.

The ODE is second order so to evolve it from $x=0$ we need to know both $u(0)$ and $u'(0)$. The first one we know, but the second one is determined by the condition that it gives $u(L) = 0$. A common method to solve such an ODE is to use a shooting method. In this method we simply guess a value for $u'(0)$, evolve the ODE to $x = L$ and check if $u(L) = 0$. If not we go back and guess another value for $u'(0)$ and repeat the process till we have convergence. The "guess" step usually means using a root finding method like bisection or the secant method to determine the new value of $u'(0)$ based on the old value and the difference $\epsilon = u(L) - 0$ (this is a root finding problem since $\epsilon = 0$ when we have the correct $u'(0)$).

As for solving the ODE you can use any method you like. For this purpose writing it as a system of first order equations

$$\frac{d}{dx}\pmatrix{u(x)\\z(x)} = \pmatrix{z(x)\\c},~~~\pmatrix{u(0)\\z(0)} = \pmatrix{u(0)\\u'(0)}$$ might be useful.

Here is an example of how to perform a shooting method using bisection in Mathematica

(* Define the system *)
L = 1;
c = 2;
odesystem[udot_] := NDSolve[{u''[x] == c, u[0] == 0, u'[0] == udot}, u, {x, 0, 1}];
uofL[udot_] := u[L] /. odesystem[udot][[1]];

(* We are to use bisection so check that u[L] < 0 for u'(0) = ulow and u[L] > 0 for u'(0) = uhigh *)
ulow = -2 c L; (* Pick some large negative value *)
uhigh = 2 c L; (* Pick some large positive value *)
Print[uofL[ulow]];
Print[uofL[uhigh]];

(* Do the shooting using bisection *)
epsilon = 0.0001; (* Precision in the u'(0) value we want: for simplicity we use this as convergence criterion *)
While[Abs[uhigh - ulow] > epsilon,

  (* The new value of u'(0) *)
  udotcur = (ulow + uhigh)/2.0;

  (* Solve the ODE and find the new value of u[L] *)
  uofLcur = uofL[udotcur];

  (* Bisection step: update ulow or uhigh depending on if uofLcur is positive or negative *)
  If[uofLcur < 0.0, ulow = udotcur, uhigh = udotcur ];
  ];

(* Output the values we find in the end *)
Print["The u'(0) value we need: ", ulow];
Print["The analytical u'(0) value: ", -c L/2];

which after $16$ iterations outputs:

The u'(0) value we need: -1.00006
The analytical u'(0) value: -1