% Poiseuille_startup.m

%

% This MATLAB script file plots the transient velocity

% profile for the start-up of pressure-driven flow in

% a pipe.  The series expansion is truncated after the

% seventh Bessel function.

%

% K. Beers. MIT 8/1/2002

% First, set the physical parameters.

G = -1;  % pressure gradient

rho = 1;  % density

mu = 1; % viscosity

nu = mu/rho;  % kinematic viscosity

Rp = 1;  % pipe radius

% Next, we need to find the first leading zeros of

% the Bessel function.

J_roots_guess = [2.4048, 5.5201, 8.6537, ...

                11.7915, 14.9309, 18.0711, 21.2116];

i_truncate = length(J_roots_guess);

J_roots = zeros(i_truncate,1);

for k=1:i_truncate

    J_roots(k) = fzero('besselj(0,x)',J_roots_guess(k));

end

% Now, we make a vector of the time values at which to

% plot the results.

time_end = 2;

num_plot = 20;

time = linspace(0,time_end,num_plot);

dt = time(2) - time(1);

% We also make a plot of the radial coordinate.

r = linspace(0,Rp,100);

% We calculate the steady-state velocity profile.

vz_inf = -G/4/mu*(Rp.^2 - r.^2);

% Calculate the phi factors of the Bessel function expansion.

Phi = zeros(i_truncate,1);

for k=1:i_truncate

    integrand = (Rp^2 - r.^2).*r.*besselj(0,J_roots(k)*r/Rp);

    integral = trapz(r,integrand);

    Phi(k) = 2/(Rp*Rp*besselj(1,J_roots(k)))^2 * integral;

end

% For each time, the transient velocity profile is calcualted

% and added to the figure;

figure;

hold on;

for itime=1:num_plot

    vz = vz_inf;

    % For each order in the expansion

    for k=1:i_truncate

        % Get time factor

        var_time = -nu*(J_roots(k)^2)/(Rp^2);

        t_factor = exp(var_time*time(itime));

        % add contribution of this order to solution

        vz = vz + (G*Rp*Rp/4/mu)*Phi(k)* ...

            t_factor*besselj(0,J_roots(k)*r/Rp);

    end

    plot(r,vz);

end

xlabel('r');

ylabel('v_z(r,t)');

phrase1 = ['Start-up of Poiseuille flow : ', ...

        '\Deltap/\Deltaz = ', num2str(G), ...

        ', \rho = ', num2str(rho), ...

        ', \mu = ', num2str(mu), ...

        ', R = ', num2str(Rp)];

title(phrase1);

phrase2 = ['\Delta t = ', num2str(dt), ...

        ', t_{end} = ', num2str(time_end), ...

        ', # of Bessel functions = ', int2str(i_truncate)];

gtext(phrase2);