Plate startup

This notebook has been written in Mathematica by 

Mark J. McCready
Professor of Chemical and Biomolecular Engineering
University of Notre Dame
Notre Dame IN 46556
USA

Mark.J.McCready.1@nd.edu
http://www.nd.edu/~mjm/


It is copyrighted to the extent allowed by whatever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy, that this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version:  12/2/17

Startup flow between two parallel plates that have a spacing of b.

Suppose that we want to solve the parallel plate flow problem where the fluid is initially at rest and then at t=0, the bottom plate is accelerated to u = u0.

You can simplify the Navier-Stokes equations to obtain the “Heat Equation”.  (So named because this problem is analogous to heat conduction in a solid where the temperature at one boundary is suddenly changed and held constant.)

DEigensystem[{-Laplacian[u[y],{y}],DirichletCondition[u[y]==0,True]},u[y],{y,0,b},4]

heqn=D[u[y,t],t]==ν  D[u[y,t],{y,2}];

The initial condition is, ux= 0

ic=u[y,0]==0;

The boundary conditions for a channel of height b are:

bc={u[0,t]==u0,u[b,t]==0};

For these boundary conditions, the problem is easily solved

sol=DSolve[{heqn,ic,bc},u[y,t],{y,t}]

We need to pick the number of terms in the sum.  You can change this to see what happens

asolp1=(u[y,t]/.sol[[1]]/.{∞->100}//Activate)/.{ν->1,b->1,u0->1};

Here is a plot of the velocity profile at three times.

Plot[{Evaluate[asolp1/.t->.01],Evaluate[asolp1/.t->.05],Evaluate[asolp1/.t->.1]},{y,0,1},PlotLegends->{"t=.01","t=.05","t=.1"},AxesLabel->{"y","ux"}]

Note the the “y” axis is the moving plate, which is moving upward.  The disturbance extends further into the fluid as the time increases, but because of the no slip condition, the velocity of the plate is always 1. 

By t =1 the disturbance is felt essentially to top wall.

Let’s use a less viscous fluid.  The efficiency of “momentum transport” should be less.

asolp1=(u[y,t]/.sol[[1]]/.{∞->100}//Activate)/.{ν->.01,b->1,u0->1};

Plot[{Evaluate[asolp1/.t->.1],Evaluate[asolp1/.t->.5],Evaluate[asolp1/.t->1]},{y,0,.5},PlotLegends->{"t=.1","t=.5","t=1"},AxesLabel->{"y","ux"},PlotRange->{0,1}]

With ν=.01, the disturbance is getting only about .35 in 1 sec.

What do these two results tell you about the time scale for this flow?

Plot3D[asolp1//Evaluate,{y,0,1},{t,0,20},Exclusions->None,PlotRange->All,AxesLabel->{"y","t","ux"}]

Let’s get the similarity solution for the startup flow

Here is the “Heat Equation”, the pde that governs the problem.  The fluid is motionless for t<0.  At t = 0, the plate moves at U.

heateq =  D[u[t,y],t] - ν   D[u[t,y],{y,2}]

We are going to hypothesize that in the absence of a geometric length scale, there will be a relation between distance and time that generates self-similar profiles in a new variable, η = y/(Sqrt[4 ν  t])

We call this eta and use as needed.  The new ode will be in the function, “f(η)”, which equals u[t,y]/U -- a nondimensional velocity.

eta = y/Sqrt[4 ν t]

Now use the chain rule on the derivatives converting these to derivatives of η.  Mathematica can use the chain rule with the new η[t,y], without evaluating the result

ueq1=heateq/.{Derivative[a1_,a2_][u][t,y]:>  (∂_({t,a1},{y,a2}))f[η[t,y]]}

Now we substitute the explicit functional form of η which simplifies the situation greatly.

ueq2 =ueq1/.Derivative[a1_,a2_][η][t,y]:>  (∂_({t,a1},{y,a2}))eta

It is easier to use a single variable for η[t,y], let’s call it ζ.

ueq3 = ueq2/.η[t,y]-> ζ

Now clear the equation.

ueq4 = ueq3/.y-> ζ 2 Sqrt[t ν]

The final ode is.

ueq5 = Expand[ueq4  4 t]

The boundary condition for vx[0,t] = U, is now just f[0] = 1,  because y = 0 is the same as η =0.  However, note that y--> Infinity and t =0 both make η--> Infinity.  This is is good because you had 2 boundary conditions and now you can use only 1 boundary condition.  Thus it is also fortunate that the other boundary conditions are the same: vx[y, t=0] ==0 and vx[infinity,t] ==0.  So the other boundary condition for f is f[∞] = 0.

This ode is easily solved with these boundary conditions.

DSolve[{ueq5==0, f[0] ==1,f[∞] ==0},f[ζ],ζ]

ans = f[ζ]/.%[[1]]

ans = 1-Erf[ζ]

Plot[ans,{ζ,0,3}, AxesLabel->{"y/\!\(\*SqrtBox[\(4\\ ν\\ t\)]\)" ,"ux/U"}]

So we see that for any combination of y and t that give η = 2, the velocity will be its undisturbed original value.  This makes it easy to get a relation for the boundary layer thickness as a function of time!