Advanced Mathematics for Engineers and Scientists/Parallel Plate Flow: Easy IC

Parallel Plate Flow: Easy IC

Formulation

As with ODEs, separation of variables is easy to understand and works well whenever it works. For ODEs, we use the substitution rule to allow antidifferentiation, but for PDEs it's a very different process involving letting dependencies pass through the partial derivatives.

A fluid mechanics example will be used.

Consider two plates parallel to each other of huge extent, separated by a distance of 1. Fluid is smoothly flowing between these two plates in only one direction (call it x). This may be seen in the picture below.

Visualization of the parallel plate flow problem.

After some assumptions, the following PDE may be obtained to describe the fluid flow:

${\frac {\partial u}{\partial t}}=\nu {\frac {\partial ^{2}u}{\partial y^{2}}}-{\frac {P_{x}}{\rho }}\,$

This linear PDE is the result of simplifying the Navier-Stokes equations, a large nonlinear PDE system which describes fluid flow. u is the velocity of the fluid in the x direction, ρ is the density of the fluid, ν is the kinematic viscosity (metaphorically speaking, how hard the molecules of the fluid rub against each other), and P_x is the pressure gradient or pressure gradient vector field. Note that u = u(y, t), there is no dependence on x. In other words, the state of the fluid upstream is no different from the state downstream. Notice that we do not consider turbulences and that the state of the flow varies between the upper plate and the lower plate as the speed u is 0 close to the plates.

u(y, t) is a velocity profile. Fluid mechanics typically is concerned with velocity fields, contrary to rigid body mechanics in which the position of an object is what is important. In other words, with a rigid object all 'points' of the object move at the same speed. With a fluid we get a velocity field and each point is moving at its own speed.

The ratio P_x/ρ describes the driving force; it's a pressure change (gradient) along the x direction. If P_x is negative, then the pressure downstream (positive x) is smaller than the pressure upstream (negative x) and the fluid will flow left to right, i.e., u(y, t) will generally be positive.

Now on to create a specific problem: let's say that a constant negative pressure gradient was applied for a long time, until the velocity profile was steady (steady means "not changing with time"). Then the pressure gradient is suddenly removed, and without this driving force the fluid will slow down and stop. That is the assumption we are going to use for our example calculations. We assume the flow to be steady and then decaying as the pressure (expressed as P_x/ρ) is removed. Hence we can remove the term -P_x/ρ from our model as well, see (PDE) below.

Initial flow profile.

Let's say that before the pressure was removed, the velocity profile was u(y, t) = sin(π y). With velocity profile we mean, as the velocities are measured on a cross section of the flow, they form a hump. This would make sense: the friction dictates less motion near the plates (see next paragraph), so we could expect a maximum velocity near the centerline (y = 1/2). This assumed profile isn't really correct, but will serve as an example for now. It's graphed at right in the domain of interest.

Before getting into the math, one more thing is needed: boundary conditions. In this case, the BC is called the no slip condition, which states that the velocity of a fluid at a wall (boundary) is equal to the velocity of the wall. If we weren't making this assumption, the fluid would be moving like a rigid object. Since the velocities of the walls (or plates) in this problem are both zero, the velocity of the fluid must be zero at these two boundaries. The BCs are then u(0, t) = 0 (bottom plate) and u(1, t) = 0 (top plate).

The IBVP is:

${\frac {\partial u}{\partial t}}=\nu {\frac {\partial ^{2}u}{\partial y^{2}}}\qquad {\mbox{(PDE)}}\,$

$u(y,0)=\sin(\pi y)\qquad {\mbox{(IC)}}\,$

$u(0,t)=0\,$

$u(1,t)=0\qquad {\mbox{(BCs)}}\,$

Notice that the initial condition IC is the same hump from the first section of this book, just turned sideways. This IC implies that the flow is already flowing when we start our calculations. We are kind of calculating the decay of the flow speed.

Separation

Variables are separated the following way: we assume that $u(y,t)=Y(y)T(t)$, where Y and T are (unknown) functions respectively of y and t. This form is substituted into the PDE:

${\frac {\partial }{\partial t}}(u(y,t))=\nu {\frac {\partial ^{2}}{\partial y^{2}}}(u(y,t))\,$

Using $u(y,t)=Y(y)T(t)$ yields:

${\frac {\partial }{\partial t}}(Y(y)T(t))=\nu {\frac {\partial ^{2}}{\partial y^{2}}}(Y(y)T(t))\,$

$Y(y){\frac {\partial }{\partial t}}(T(t))=\nu T(t){\frac {\partial ^{2}}{\partial y^{2}}}(Y(y))\,$

$Y(y)T'(t)=\nu T(t)Y''(y)\,$

${\frac {T'(t)}{\nu T(t)}}={\frac {Y''(y)}{Y(y)}}\,$

Look carefully at the last equation: the left side of the equation depends strictly on t, and the right side strictly on y, and they are equal. t may be varied independently of y and they'd still be equal, and y may be varied independently of t and they'd still be equal. This can only happen if both sides are constant. This may be shown as follows:

${\frac {d}{dt}}\left({\frac {T'(t)}{\nu T(t)}}\right)={\frac {d}{dt}}\left({\frac {Y''(y)}{Y(y)}}\right)\,$

${\frac {d}{dt}}\left({\frac {T'(t)}{\nu T(t)}}\right)=0\,$

Taking the derivative of both sides, turns the right-hand-side into a constant, 0. The left-hand-side is left as is. Then taking the integrals of both sides yields:

${\frac {T'(t)}{\nu T(t)}}={\mbox{a constant}}\,$

Integration of the ordinary derivative recovers the left side but leaves the right side a constant. It follows by similarity that Y''/Y is a constant as well.

The constant in question is called the separation constant. We could simply give it a letter, such as A, but a good choice of the constant will make work easier later. In this case the best choice is -k². This will be justified later (but it should be reemphasized that it may be notated any way you want, assuming it can span the domain).

${\frac {T'(t)}{\nu T(t)}}={\frac {Y''(y)}{Y(y)}}=-k^{2}\,$

${\Big \Downarrow }$

${\frac {T'(t)}{\nu T(t)}}=-k^{2}\,$

${\frac {Y''(y)}{Y(y)}}=-k^{2}\,$

The variables are now separated. The last two equations are two ODEs which may be solved independently (in fact, the Y equation is an eigenvalue problem), though they both contain an unknown constant. Note that ν was kept for the T equation. This choice makes the solution slightly easier, but is again completely arbitrary.

${\frac {T'(t)}{\nu T(t)}}=-k^{2}\,$

Rearrange and note that the notion buried in the expression below is that a function is equal to its own derivative, which kind of strikes the Euler number bell.

$T'(t)=-k^{2}\nu T(t)\,$

Then solve.

$T(t)=C_{1}e^{-k^{2}\nu t}\,$

Here is a solution plucked from Ted Woollett's 'Maxima by Example', Chapter 3, '3.2.3 Exact Solution Using desolve'. Maxima's output is not shown but should be easily reconcilable with what has been written.

(%i1) de:(-%k^2*v*T(t))-('diff(T(t),t)) = 0;

(%i2) gsoln:desolve(de,T(t));

The same goes for our right-hand-side.

${\frac {Y''(y)}{Y(y)}}=-k^{2}\,$

$Y''(y)=-k^{2}Y(y)\,$

And solve.

$Y(y)=C_{2}\cos(ky)+C_{3}\sin(ky)\,$

In Maxima the solution goes like this. When Maxima is asking for 'zero' or 'nonzero' then type 'nonzero;':

(%i1)de:(-%k^2*Y(y))-('diff(Y(y),y,2));

(%i2)atvalue('diff(Y(y),y), y=0, Y(0)*%k )$

(%i3)gsoln:desolve(de,Y(y));

The overall solution will be $u(y,t)$, still with unknown constants, and as of now the product of Y and T. So we are plugging the partial solutions for $Y(y)$ and $T(t)$ back in:

$u(y,t)=Y(y)T(t)=C_{1}e^{-k^{2}\nu t}(C_{2}\cos(ky)+C_{3}\sin(ky))$

$u(y,t)=e^{-k^{2}\nu t}(A\cos(ky)+B\sin(ky))$

Note that C₁ has been multiplied into C₂ and C₃, reducing the number of arbitrary constants. Because an unknown constant multiplied by an unknown constant yields still an unknown constant.

The IC or BCs should now be applied. If the IC was applied first, coefficients would be equated and all of the constants would be determined. However, the BCs may or may not have been fulfilled (in this case they would, but you're not generally so lucky). So to be safe, the BCs will be applied first:

$u(0,t)=e^{-k^{2}\nu t}(A\cos(k\cdot 0)+B\sin(k\cdot 0))=0$

$A\cdot 1+B\cdot 0=0\Rightarrow A=0\qquad {\mbox{(no slip along lower boundary)}}\,$

So A being zero eliminates the $\cos(k\cdot 0)$ term.

$u(1,t)=e^{-k^{2}\nu t}B\sin(k\cdot 1)=0\,$

$B\sin(k\cdot 1)=0\Rightarrow B=0\quad {\mbox{or}}\quad k=n\pi \qquad {\mbox{(no slip along upper boundary)}}\,$

If we took B = 0, the solution would have just been u(y, t) = 0 (often called the trivial solution), which would satisfy the BCs and the PDE but couldn't possibly satisfy the IC. So, we take k = nπ instead, where n is any integer. After applying the BCs we have:

Decaying flow.

$u(y,t)=e^{-(n\pi )^{2}\nu t}B\sin(n\pi y)\,$

Then we need to apply the IC to it. Per the IC from above $u(y,0)$ is:

$u(y,0)=\sin(\pi y)\,$

Per the BC $u(y,0){\mbox{with t = 0}}$ is, see above. Setting those two definitions of $u(y,0)$ equal is:

$e^{0}B\sin(n\pi y)=\sin(\pi y)\,$

Since t = 0 as per the IC assumption, $e^{-(n\pi )^{2}\nu t}$ is becoming $e^{0}$. The equality can only hold if B = 1 and n = 1, allowing us to simply remove those two constants from the function in its BC incarnation. That's it! The complete solution is:

$u(y,t)=e^{-\pi ^{2}\nu t}\sin(\pi y)\,$

It's worth verifying that the IC, BCs, and PDE are all satisfied by this. Also notice that the solution is a product of a function of t and a function of y. The graph at the right is illustrating this. Observe that the profile is plotted for different values of νt, rather than specifying some ν and graphing different values for t. Remember that v is the kinematic viscosity. Hence the decay of flow speed also depends on it. Looking at the solution, t and ν appear only once and they're multiplying, so it's natural to do this. A dimensionless time could have been introduced from the beginning.

So what happens? The fluid starts with its initial profile and slows down exponentially. Note that with x replaced with y and t replaced with νt, this is exactly the same as the result of heat flow in a rod as shown in the introduction. This is not a coincidence: the PDE for the rod describes diffusion of heat, the PDE for the parallel plates describes diffusion of momentum.

Take a second look at the separation constant, -k². The square is convenient, without it the solution for Y(y) would have involved square roots. Without the negative sign, the solution would have involved exponentials instead of sinusoids, so the constant would have come out imaginary.

The assumption that u(y, t) = Y(y)T(t) is justified by the physics of the problem: it would make sense that the profile would keep its general shape (due to Y(y)), only it'd get flattened over time as the fluid slows down (due to T(t)).

Quickly recap what we did. First we took a model PDE from Navier-Stokes and simplified it by making an assumption about the pressure being removed. Then we sort of first-stage solved it by separating Y and T. Then we applied the BC (boundary conditions) and the IC (initial condition) yielding our final solution.