% polymer_flow_1D.m
% This MATLAB program uses finite differences
% to model the flow rate of a poly(styrene)
% melt between two parallel flat plates,
% using the Carreau-Yasuda model of a shear
% thinning fluid.
% This version uses one-sided approximations
% for the first-derivatives to avoid oscillations
% and provides a sparsity pattern to fsolve to
% speed up Jacobian calculations.
% K.J. Beers. MIT ChE. 2/17/2005

function iflag_main = polymer_flow_1D();
iflag_main = 0;

% set the fluid parameters of poly(styrene)
% melt at 453K and plot predicted shear-rate
% dependent viscosity
Fluid.eta_0 = 1.48e4;  % in Pa*s
Fluid.eta_inf = 0;
Fluid.lambda = 1.04; % in s
Fluid.a = 2;
Fluid.n = 0.398;

% Ask the user to input the problem
% specifications
Sys.B = input(...
    'Enter distance in m between plates : ');

% Ask the user for a target mid-point velocity.
% Then, for reference, the pressure gradient
% necessary to achieve this velocity, neglecting
% shear-thinning effects, is printed to the screen.
v_tar = input(...
    'Enter target mid-point velocity in m/s : ');
dP_dx_tar = -8*Fluid.eta_0*v_tar/Sys.B/Sys.B;
disp(' ');
disp(...
    'Neglecting shear-thinning effects, the dP_dx');
disp(['necessary to achieve this velocity is : ', ...
    num2str(dP_dx_tar)]);
disp(' ');
% Then, compute approximate Reynolds' # to check
% for whether these flow conditions are at all
% consistent with laminar flow
Re_tar = 1000*v_tar*Sys.B/Fluid.eta_0;
disp(['approx. target Reynolds number : ', ...
    num2str(Re_tar)]);
disp(' ');
Sys.dP_dx = input(...
    'Enter dynamic pressure gradient in Pa/m : ');

% ask user to input grid parameters
Grid.N = input(...
    'Enter number of grid points : ');
% generate the grid data
Grid.dy = Sys.B/(Grid.N + 1);
Grid.y = linspace(Grid.dy,Sys.B-Grid.dy,Grid.N);

% Make plot of shear-dependent viscosity
gammaDotVals = logspace(-1,2,250);
viscVals = calc_viscosity(gammaDotVals,Fluid);
figure;  loglog(gammaDotVals,viscVals);
xlabel('Shear rate in s^{-1}');
ylabel('Viscosity in Pa*s');
title('Viscosity of poly(styrene) melt at 453K');

% make initial guess based on velocity profile
% of Newtonian fluid with the zero shear-rate
% viscosity.
v_0 = 1/2/Fluid.eta_0*Sys.dP_dx.* ...
    (Grid.y.*(Grid.y - Sys.B));

% generate matrix showing sparsity
% pattern of the Jacobian
Jac_S = Jac_pattern(Grid);

% set solver parameters
Options = optimset('LargeScale','off', ...
    'JacobPattern', Jac_S);
% call fsolve
[v, f] = fsolve(@calc_f, v_0, Options, ...
    Sys, Fluid, Grid);

% check infinity norm of solution vector
f_norm = norm(f,inf);
disp(['function infinity norm = ', ...
    num2str(f_norm)]);

% compute for the solution the shear rate
% and viscosity at each grid point
gamma_dot = calc_gamma_dot(v, Grid);
eta = calc_viscosity(gamma_dot, Fluid);
shear_stress = Sys.dP_dx.*(Grid.y - Sys.B/2);

% compute local Reynolds' numbers and return
% max. value as check on validity of
% laminar flow calculations
Re = 1000.*(v.^2)./eta./gamma_dot;
Re_max = max(Re);
disp(['max. local Reynolds number = ', ...
    num2str(Re_max)]);

% plot the results
figure;  subplot(2,2,1);  plot(v,Grid.y);
hold on;  plot(v_0, Grid.y, '-.');
xlabel('v_x(y) (m/s)');   ylabel('y (m)');
title('1-D laminar flow of PS melt at 453K');
subplot(2,2,2);  plot(shear_stress, Grid.y);
xlabel('\tau_{yx} (Pa)');  ylabel('y (m)');
subplot(2,2,3);  plot(gamma_dot, Grid.y);
xlabel('shear rate (s^{-1})');  ylabel('y (m)');
subplot(2,2,4);  plot(eta, Grid.y);
xlabel('viscosity (Pa*s)');  ylabel('y (m)');

% save results in .mat file
save polymer_flow_1D_results.mat;

iflag_main = 1;
return;


% ====================
% This routine estimates the shear-rate dependent
% viscosity of a fluid using the Carreau-Yasuda
% model.
function visc = calc_viscosity(gamma_dot,Fluid);

var1 = 1 + (Fluid.lambda.*gamma_dot).^Fluid.a;
var2 = var1.^((Fluid.n - 1)/Fluid.a);

visc = Fluid.eta_inf + ...
    (Fluid.eta_0 - Fluid.eta_inf).*var2;
return;

% ====================
% This routine returns a sparse matrix showing
% the locations of the non-zero elements in
% the Jacobian
function Jac_S = Jac_pattern(Grid);

Jac_S = spalloc(Grid.N,Grid.N,2*Grid.N);
for k=1:Grid.N
    Jac_S(k,k) = 1;
    if(k > 1)
        Jac_S(k,k-1) = 1;
    end
end

return;

% ====================
% This routine returns the vector of function
% values for the set of nonlinear algebraic equations.
function f = calc_f(v, Sys, Fluid, Grid);

f = zeros(Grid.N,1);
Jac = zeros(Grid.N,Grid.N);

% compute shear rates at each grid point
gamma_dot = calc_gamma_dot(v, Grid);

% compute viscosities at each grid point
eta = calc_viscosity(gamma_dot, Fluid);

var1 = 1/Grid.dy;
% compute function values
f(1) = eta(1)*v(1)*var1 ...
    - Sys.dP_dx*(Grid.y(1) - Sys.B/2);
for k=2:Grid.N
    f(k) = eta(k)*(v(k)-v(k-1))*var1 ...
        - Sys.dP_dx*(Grid.y(k) - Sys.B/2);
end

return;


% ====================
% This routine estimates the shear rates
% at each grid point.
function gamma_dot = calc_gamma_dot(v, Grid);

gamma_dot = zeros(size(v));
gamma_dot(1) = v(1);
for k=2:Grid.N
    gamma_dot(k) = v(k) - v(k-1);
end
gamma_dot = gamma_dot.*(1/Grid.dy);
gamma_dot = abs(gamma_dot);

return;