Plane poiseuelle flow solved by finite difference

September 30, 2011 at 01:05 PM | categories: odes | View Comments

poiseuille_flow_finite_difference

Contents

Plane poiseuelle flow solved by finite difference

John Kitchin

http://www.physics.arizona.edu/~restrepo/475B/Notes/source/node24.html

solve a linear boundary value problem of the form: y'' = p(x)y' + q(x)y + r(x) with boundary conditions y(x1) = alpha and y(x2) = beta.

For this example, we resolve the plane poiseuille flow problem we previously solved in Post 878 with the builtin solver bvp5c, and in Post 1036 by the shooting method. An advantage of the approach we use here is we do not have to rewrite the second order ODE as a set of coupled first order ODEs, nor do we have to provide guesses for the solution.

we want to solve u'' = 1/mu*DPDX with u(0)=0 and u(0.1)=0. for this problem we let the plate separation be d=0.1, the viscosity $\mu = 1$, and $\frac{\Delta P}{\Delta x} = -100$.

clc; close all; clear all;

we define function handles for each term in the BVP.

we use the notation for y'' = p(x)y' + q(x)y + r(x)

p = @(x) 0;
q = @(x) 0;
r = @(x) -100;

define the boundary conditions

we use the notation y(x1) = alpha and y(x2) = beta

x1 = 0; alpha = 0;
x2 = 0.1; beta = 0;

define discretization

here we define the number of intervals to discretize the x-range over. You may need to increase this number until you see your solution no longer depends on the number of discretized points.

npoints = 10;

No user modification required below here

Below here is just the algorithm for solving the finite difference problem. To solve another kind of linear BVP, just modify the variables above according to your problem.

The idea behind the finite difference method is to approximate the derivatives by finite differences on a grid. See here for details. By discretizing the ODE, we arrive at a set of linear algebra equations of the form $A y = b$, where $A$ and $b$ are defined as follows.

To solve our problem, we simply need to create these matrices and solve the linear algebra problem.

% compute interval width
h = (x2-x1)/npoints;

% preallocate and shape the b vector and A-matrix
b = zeros(npoints-1,1);
A = zeros(npoints-1, npoints-1);
X = zeros(npoints-1,1);

now we populate the A-matrix and b vector elements

for i=1:(npoints-1)
    X(i,1) = x1 + i*h;

    % get the value of the BVP Odes at this x
    pi = p(X(i));
    qi = q(X(i));
    ri = r(X(i));

    if i == 1
        b(i) = h^2*ri-(1-h/2*pi)*alpha; % first boundary condition
    elseif i == npoints - 1
        b(i) = h^2*ri-(1+h/2*pi)*beta; % second boundary condition
    else
        b(i) = h^2*ri; % intermediate points
    end

    for j = 1:(npoints-1)
        if j == i % the diagonal
            A(i,j) = -2+h^2*qi;
        elseif j == i-1 % left of the diagonal
            A(i,j) = 1-h/2*pi;
        elseif j == i+1 % right of the diagonal
            A(i,j) = 1 + h/2*pi;
        else
            A(i,j) = 0; % off the tri-diagonal
        end
    end
end

solve the equations A*y = b for Y

Y = A\b;

construct the solution

the equations were only solved for X between x1 and x2, so we add the boundary conditions back to the solution on each side

x = [x1; X; x2];
y = [alpha; Y; beta];

plot the numerical solution

plot(x,y)
xlabel('distance between plates')
ylabel('fluid velocity')
% the analytical solution to this problem is known, and we can plot it on
% top of the numerical solution
hold all
mu = 1;
d = 0.1
x = linspace(0,0.1);
Pdrop = -100; % this is DeltaP/Deltax
u = -(Pdrop)*d^2/2/mu*(x/d-(x/d).^2);
plot(x,u,'r--')

d =

    0.1000

The two solutions agree. The advantages of this method over the shooting method (:postid:1036`) and builtin solver ( Post 878 ) are that you do not need to break the 2nd order ODE into a system of 1st order ODEs, and there is no "guessing" of the solution to get the solver started. Some disadvantages are the need to have a linear BVP in the form described at the beginning. If your BVP is nonlinear, you can't solve it with linear algebra! Another disadvantage is the amount of code we had to write. If you need to solve this kind of problem often, you probably would want to encapsulate the algorithm code into a reusable function.

'done'

% categories: ODEs
% tags: BVP

% post_id = 1193; %delete this line to force new post;
% permaLink = http://matlab.cheme.cmu.edu/2011/09/30/plane-poiseuelle-flow-solved-by-finite-difference/;

ans =

done


Published with MATLAB® 7.12

Read and Post Comments

Plane Poiseuille flow - BVP solve by shooting method

September 08, 2011 at 06:24 PM | categories: odes | View Comments

Plane Poiseuille flow - BVP solve by shooting method

Plane Poiseuille flow - BVP solve by shooting method

In Post 878 learned how to use the BVP solver in Matlab to solve a boundary value problem. Another approach is to use the shooting method. The reason we can't use an initial value solver for a BVP is that there is not enough information at the initial value to start. In the shooting method, we take the function value at the initial point, and guess what the function derivatives are so that we can do an integration. If our guess was good, then the solution will go through the known second boundary point. If not, we guess again, until we get the answer we need. In this example we repeat the pressure driven flow example, but illustrate the shooting method.

In the pressure driven flow of a fluid with viscosity $\mu$ between two stationary plates separated by distance $d$ and driven by a pressure drop $\Delta P/\Delta x$, the governing equations on the velocity $u$ of the fluid are (assuming flow in the x-direction with the velocity varying only in the y-direction):

Contents

$$\frac{\Delta P}{\Delta x} = \mu \frac{d^2u}{dy^2}$$

with boundary conditions $u(y=0) = 0$ and $u(y=d) = 0$, i.e. the no-slip condition at the edges of the plate.

we convert this second order BVP to a system of ODEs by letting $u_1 = u$, $u_2 = u_1'$ and then $u_2' = u_1''$. This leads to:

$\frac{d u_1}{dy} = u_2$

$\frac{d u_2}{dy} = \frac{1}{\mu}\frac{\Delta P}{\Delta x}$

with boundary conditions $u_1(y=0) = 0$ and $u_1(y=d) = 0$.

for this problem we let the plate separation be d=0.1, the viscosity $\mu = 1$, and $\frac{\Delta P}{\Delta x} = -100$.

function main

clc; close all; clear all;

d = 0.1; % plate thickness

First guess

we need u_1(0) and u_2(0), but we only have u_1(0). We need to guess a value for u_2(0) and see if the solution goes through the u_2(d)=0 boundary value.

u1_0 = 0; % known
u2_0 = 1; % guessed

init = [u1_0 u2_0];
dspan = [0 d];
[y,U] = ode45(@odefun,dspan,init);
plot(y,U(:,1),[d],[0],'ro')
legend 'solution #1' 'u_2(d)=0'
xlabel('y')
ylabel('u_1')
% It appears we have undershot, so our initial slope was not large
% enough

second guess

u1_0 = 0; % known
u2_0 = 10; % guessed

init = [u1_0 u2_0];
dspan = [0 d];
[y,U] = ode45(@odefun,dspan,init);
plot(y,U(:,1),[d],[0],'ro')
legend 'solution #2' 'u_2(d)=0'
xlabel('y')
ylabel('u_1')
% Now we have overshot, the initial slope is too high. The right
% answer must be between 1 and 10 You could continue to iterate on
% this manually until the value of u_1(d) = 0. Instead, we use a
% solver to do that, now that we have a reasonable guess for u_2(0).

Using a solver

we have to make a function that is equal to zero when u_2(d) = 0.

u2_guess = 2;
u2_0_solved = fsolve(@myfunc,u2_guess)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.




u2_0_solved =

    5.0000

Plot final solution

init = [u1_0 u2_0_solved];
dspan = [0 d];
[y,U] = ode45(@odefun,dspan,init);

figure
plot(y,U(:,1),d,0,'ro')
xlabel('y')
ylabel('u_1')

Known analytical solution

the analytical solution to this problem is known, and we can plot it on top of the numerical solution

mu = 1;
Pdrop = -100; % this is DeltaP/Deltax
u = -(Pdrop)*d^2/2/mu*(y/d-(y/d).^2);
hold all
plot(y,u,'r--')
legend 'solution' 'u_2(d)=0' 'analytical solution'

'done'

ans =

done

Finally, the solutions agree, and show that the velocity profile is parabolic, with a maximum in the middle of the two plates. You can also see that the velocity at each boundary is zero, as required by the boundary conditions

function dUdy = odefun(y,U)
% here we define our differential equations
u1 = U(1);
u2 = U(2);

mu = 1;
Pdrop = -100; % this is DeltaP/Deltax
du1dy = u2;
du2dy = 1/mu*Pdrop;
dUdy = [du1dy; du2dy];


function Z = myfunc(u2_guess)
% this function is zero when u_1(d) = 0.
u1_0 = 0; % known
u2_0 = u2_guess; % guessed
d = 0.1;

init = [u1_0 u2_0];
dspan = [0 d];
sol = ode45(@odefun,dspan,init);

u = deval(sol,d); % u = [u_1(d) u_2(d)]
u1_d = u(1);
Z = u1_d;

% categories: ODEs
% tags: math, fluids, bvp


% post_id = 1036; %delete this line to force new post;


Published with MATLAB® 7.12

Read and Post Comments