LDPE is flowing in a die that has a height of 1mm (H=1mm/2), a width of 30mm(W), and is 60mm long (L). Use Shooting method to find the flow rate if the pressure difference is 820 psi over the length of 60 mm.

SOLUTION: Since alpha is unknown, it is initially guessed. The equations are then solved. The Newton_Ralphson method is then used to check if the velocity in the x direction at y=H is 0, (a sensivity equation with respect to alpha is derived for this Newton-Ralphson method). If not, alpha is adjusted and the integration is repeated.

The code for solving using Matlab

------THE MAIN FILE

%Use Newton-Ralphson method to check for V1=0 at y=H

%We want V1(H)=f(alpha)=0

%Iterate on alpha

%set guess of alpha before calling

global alpha W

alpha = 1;

diff = 1;

while (diff>1.e-8)

problem_22

nn = size(V);

n = nn(1) ; %first value of V

fn = V(n,1) ; %last value of V1

fnprime = V(n,3) ; %last value of V3

diff = fn/fnprime;

alpha = alpha - diff;

diff = abs(diff);

end

Velocity_1 = V(:,1); %values of V1= u = velocity in x

direction, dV1_dy=shear rate

Velocity_5 = V(:,5);

flowrate=W.*Velocity_5;

s=size(V,1); %return the number of rows in V

Flow_rate=flowrate(s,1)

figure(1)

plot(y,Velocity_1,'r')

grid

xlabel('POSITION, y (m)');

ylabel('VELOCITY IN X DIRECTION, V1=u (m/s)');

-----THE "PROBLEM" FILE

global alpha

%set initial conditions

V0(1) = alpha;

V0(2) = 0;

V0(3) = 1;

V0(4) = 0;

V0(5)=0;

H=1.e-3; %m

%integrate

[y,V]=ode45('ode_22',0,H/2,V0, 1.e-15);

----THE "ode_22" FILE

function dV_dy=ode_22(y,V)

global alpha W

no=8818;

%initial viscosity, in Pa.s

deltaP=820*101325/14.696; %convert from psi to Pa

L=60.e-3;

W=30.e-3;

%-------------------------------------

dV_dy(1)=V(2);

dV_dy(2)=-deltaP/(L*no);

dV_dy(3)=V(4);

dV_dy(4)=0;