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 Finite Difference method to find

the flow rate if the pressure difference is 820 psi over the length of

60 mm

SOLUTION:

This method is to convert the differential equations into algebraic equations in terms of solution at each grid point. Then, the functions are expanded in Taylor series. These algebraic equations form a tridiagonal matrix, which is solved by using the "tridiag" function as shown below.

The code for solving using Matlab:

--------------THE MAIN FILE

global afd bfd cfd dfd

H=0.001/2; %mm

delta_y=H/20;

delta_P=820*101325/14.696;

L=0.06;

visc_o=8818;

D=(-1)*(delta_y^2)*(delta_P)/(L*visc_o);

for i=1:21

y(i)=delta_y*i;

afd(i)=1;

bfd(i)=-2;

cfd(i)=1;

dfd(i)=D;

end

afd(n)=0;

cfd(1)=2;

dfd(n)=0;

tridiag(n)

dfd'

-----THE "function tridiag" FILE

function tridiag(n)

global afd bfd cfd dfd

% this subroutine solves a tridiagonal system of

% equations of the type that often occur with the

% finite difference method.

% input

% a(n), b(n), c(n), d(n) in the following equaiton:

% a(i)*y(i-1) + b(i)*y(i) + c(i)*y(i+1) = d(i)

% n - the size of the system being solved.

% output

% afd(n), bfd(n), cfd(n) this is the lu decomposition.

% the original afd, bfd, cfd are destroyed.

% dfd(n) is the solution.

for l=2:n

s = afd(l)/bfd(l-1);

bfd(l) = bfd(l)-s*cfd(l-1);

afd(l) = s;

end

% this is the lower decomposition. the routine could

% be separated here if you have several right-hand sides

% with the same matrix.

% forward sweep tridiagonal matrix

for l=2:n

dfd(l) = dfd(l)-afd(l)*dfd(l-1);

end

% back substitution

dfd(n) = dfd(n)/bfd(n);

for l=2:n

k = n-l+1;

dfd(k) = (dfd(k)-cfd(k)*dfd(k+1))/bfd(k);