Star-up Flow Poiseuille w/Matlab finite difference

[longbow75]

Homework Statement

I have this situation. A fluid is at rest at t=0 then a constant dp/dz = 10 is applied. Fluid starts moving.

Homework Equations

dV/dt = -(1/rho)dP/dz + viscosity*[1/r(d/dr(rdV/dr))]
Initial condition t=0 dp/dz=10
boundary condition1 r=1 V=0
boundary condition2 r=0 dV/dr =0

The Attempt at a Solution

First off finite difference eq yields [viscosity*dt/dr^2] I simply said this is beta and is equal to 0.1

I wrote the following code to calculate velocity and plot it:
clear;clc;
h = 1;
dr = 0.01; dt = 0.01;
nmax = 4000; epsilon = 0.0001;
r = [0:dr:h];
t = [0:dt:h];
beta=0.1;
n1 = fix(h/Dr)+1;
m1 = fix(h/Dt)+1;
% Initialize solution matrices
velocity = zeros(n1,m1)
velocityN =velocity
% Iterative process for the solution
for k = 1:nmax
% boundary conditions:
velocityN(:,n1) = zeros(n1,1);

% Calculate velocityN in interior points
for i = 2:n1-1
for j = 1:m1
velocityN(i,j+1)= (beta*(1-(1/(2*i)))*velocity(i-1,j)+(1-2*beta)*velocity(i,j)+beta*(1+1/(2*i))*velocity(i+1,j));
end;
end;
% Check convergence
nfail = 0;
for i = 1:n1
for j = 1:m1
if abs(velocityN(i,j)-velocity(i,j))>epsilon
nfail = nfail+1;
end;
end;
end;

if nfail > 0
velocity = velocityN;
else
return;
end;
end;
fprintf('No convergence after %i iterations.',k);
[X,Y] = meshgrid(r,t);
contour(X,Y,velocity',10)


however, I'm having hard time trying to write my bc's.
I think the first one is good at r=1 V=0
velocityN(:,n1) = zeros(n1,1);

however how do I incorporate the second bc which is dV/dr=0 at r=0
finite diff of dV/dr should be velocityN(0,j+1)=((1-4*beta)*velocity(0,j)+4*beta(velocity(1,j)))
I believe, since the code gives me an error saying that velocity (0,0) is not defined which should be defined by this equation. Same goes to initial condition as well. Where do I define it.

Many thanks for help.

 

[KataKoniK]

Thank you for sharing your situation and code. It seems like you are on the right track with your approach to solving this problem. In order to incorporate the second boundary condition of dV/dr = 0 at r = 0, you can use a central difference approximation for the derivative at that point, which would look something like this:

velocityN(1,j+1) = (1-2*beta)*velocity(1,j) + 2*beta*velocity(2,j)

This will allow you to solve for the velocity at r = 0 based on the values at r = 1 and r = 2. Similarly, for the initial condition, you can use a forward difference approximation for the derivative, which would look like this:

velocity(1,j+1) = velocity(1,j) + dt*(10/rho)

This will allow you to define the initial velocity based on the initial pressure gradient of 10. I hope this helps. Good luck with your code!