Water gas shift equilibria via the NIST Webbook

December 13, 2011 at 01:49 AM | categories: miscellaneous, nonlinear algebra | View Comments

Water gas shift equilibria via the NIST Webbook

Water gas shift equilibria via the NIST Webbook

John Kitchin

The NIST webbook provides parameterized models of the enthalpy, entropy and heat capacity of many molecules. In this example, we will examine how to use these to compute the equilibrium constant for the water gas shift reaction $CO + H_2O \rightleftharpoons CO_2 + H_2$ in the temperature range of 500K to 1000K.

Parameters are provided for:

Contents

   Cp = heat capacity (J/mol*K)
   H° = standard enthalpy (kJ/mol)
   S° = standard entropy (J/mol*K)

with models in the form: $Cp^\circ = A + B*t + C*t^2 + D*t^3 + E/t^2$

$H^\circ - H^\circ_{298.15}= A*t + B*t^2/2 + C*t^3/3 + D*t^4/4 - E/t + F - H$

$S^\circ = A*ln(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G$

where $t=T/1000$, and $T$ is the temperature in Kelvin. We can use this data to calculate equilibrium constants in the following manner. First, we have heats of formation at standard state for each compound; for elements, these are zero by definition, and for non-elements, they have values available from the NIST webbook. There are also values for the absolute entropy at standard state. Then, we have an expression for the change in enthalpy from standard state as defined above, as well as the absolute entropy. From these we can derive the reaction enthalpy, free energy and entropy at standard state, as well as at other temperatures.

We will examine the water gas shift enthalpy, free energy and equilibrium constant from 500K to 1000K, and finally compute the equilibrium composition of a gas feed containing 5 atm of CO and H2 at 1000K.

close all; clear all; clc;

T = linspace(500,1000); % degrees K
t = T/1000;

Setup equations for each species

First we enter in the parameters and compute the enthalpy and entropy for each species.

Hydrogen

link

% T = 298-1000K valid temperature range
A =  33.066178;
B = -11.363417;
C =  11.432816;
D = -2.772874;
E = -0.158558;
F = -9.980797;
G =  172.707974;
H =  0.0;

Hf_29815_H2 = 0.0; % kJ/mol
S_29815_H2 = 130.68; % J/mol/K

dH_H2 = A*t + B*t.^2/2 + C*t.^3/3 + D*t.^4/4 - E./t + F - H;
S_H2 = (A*log(t) + B*t + C*t.^2/2 + D*t.^3/3 - E./(2*t.^2) + G);

H2O

link Note these parameters limit the temperature range we can examine, as these parameters are not valid below 500K. There is another set of parameters for lower temperatures, but we do not consider them here.

% 500-1700 K valid temperature range
A =   30.09200;
B =   6.832514;
C =   6.793435;
D =  -2.534480;
E =   0.082139;
F =  -250.8810;
G =   223.3967;
H =  -241.8264;

Hf_29815_H2O = -241.83; %this is Hf.
S_29815_H2O = 188.84;

dH_H2O = A*t + B*t.^2/2 + C*t.^3/3 + D*t.^4/4 - E./t + F - H;
S_H2O = (A*log(t) + B*t + C*t.^2/2 + D*t.^3/3 - E./(2*t.^2) + G);

CO

link

% 298. - 1300K valid temperature range
A =   25.56759;
B =   6.096130;
C =   4.054656;
D =  -2.671301;
E =   0.131021;
F =  -118.0089;
G =   227.3665;
H = -110.5271;

Hf_29815_CO = -110.53; %this is Hf kJ/mol.
S_29815_CO = 197.66;

dH_CO = A*t + B*t.^2/2 + C*t.^3/3 + D*t.^4/4 - E./t + F - H;
S_CO = (A*log(t) + B*t + C*t.^2/2 + D*t.^3/3 - E./(2*t.^2) + G);

CO2

link

% 298. - 1200.K valid temperature range
A =   24.99735;
B =   55.18696;
C =  -33.69137;
D =   7.948387;
E =  -0.136638;
F =  -403.6075;
G =   228.2431;
H =  -393.5224;

Hf_29815_CO2 = -393.51; %this is Hf.
S_29815_CO2 = 213.79;

dH_CO2 = A*t + B*t.^2/2 + C*t.^3/3 + D*t.^4/4 - E./t + F - H;
S_CO2 = (A*log(t) + B*t + C*t.^2/2 + D*t.^3/3 - E./(2*t.^2) + G);

Standard state heat of reaction

We compute the enthalpy and free energy of reaction at 298.15 K for the following reaction $CO + H2O \rightleftharpoons H2 + CO2$.

Hrxn_29815 = Hf_29815_CO2 + Hf_29815_H2 - Hf_29815_CO - Hf_29815_H2O;
Srxn_29815 = S_29815_CO2 + S_29815_H2 - S_29815_CO - S_29815_H2O;
Grxn_29815 = Hrxn_29815 - 298.15*(Srxn_29815)/1000;

sprintf('deltaH = %1.2f',Hrxn_29815)
sprintf('deltaG = %1.2f',Grxn_29815)

ans =

deltaH = -41.15


ans =

deltaG = -28.62

Correcting $\Delta H$ and $\Delta G$

we have to correct for temperature change away from standard state. We only correct the enthalpy for this temperature change. The correction looks like this:

$$ \Delta H_{rxn}(T) = \Delta H_{rxn}(T_{ref}) + \sum_i \nu_i (H_i(T)-H_i(T_{ref}))$$

Where $\nu_i$ are the stoichiometric coefficients of each species, with appropriate sign for reactants and products, and $(H_i(T)-H_i(T_{ref})$ is precisely what is calculated for each species with the equations

Hrxn = Hrxn_29815 + dH_CO2 + dH_H2 - dH_CO - dH_H2O;

The entropy is on an absolute scale, so we directly calculate entropy at each temperature. Recall that H is in kJ/mol and S is in J/mol/K, so we divide S by 1000 to make the units match.

Grxn = Hrxn - T.*(S_CO2 + S_H2 - S_CO - S_H2O)/1000;

Plot how the $\Delta G$ varies with temperature

over this temperature range the reaction is exothermic, although near 1000K it is just barely exothermic. At higher temperatures we expect the reaction to become endothermic.

figure; hold all
plot(T,Grxn)
plot(T,Hrxn)
xlabel('Temperature (K)')
ylabel('(kJ/mol)')
legend('\Delta G_{rxn}', '\Delta H_{rxn}', 'location','best')

Equilibrium constant calculation

Note the equilibrium constant starts out high, i.e. strongly favoring the formation of products, but drops very quicky with increasing temperature.

R = 8.314e-3; %kJ/mol/K
K = exp(-Grxn/R./T);

figure
plot(T,K)
xlim([500 1000])
xlabel('Temperature (K)')
ylabel('Equilibrium constant')

Equilibrium yield of WGS

Now let's suppose we have a reactor with a feed of H2O and CO at 10atm at 1000K. What is the equilibrium yield of H2? Let $\epsilon$ be the extent of reaction, so that $F_i = F_{i,0} + \nu_i \epsilon$. For reactants, $\nu_i$ is negative, and for products, $\nu_i$ is positive. We have to solve for the extent of reaction that satisfies the equilibrium condition.

A = CO
B = H2O
C = H2
D = CO2
Pa0 = 5; Pb0 = 5; Pc0 = 0; Pd0 = 0;  % pressure in atm
R = 0.082;
Temperature = 1000;

% we can estimate the equilibrium like this. We could also calculate it
% using the equations above, but we would have to evaluate each term. Above
% we simple computed a vector of enthalpies, entropies, etc...
K_Temperature = interp1(T,K,Temperature);

% If we let X be fractional conversion then we have $C_A = C_{A0}(1-X)$,
% $C_B = C_{B0}-C_{A0}X$, $C_C = C_{C0}+C_{A0}X$, and $C_D =
% C_{D0}+C_{A0}X$. We also have $K(T) = (C_C C_D)/(C_A C_B)$, which finally
% reduces to $0 = K(T) - Xeq^2/(1-Xeq)^2$ under these conditions.

f = @(X) K_Temperature - X^2/(1-X)^2;

Xeq = fzero(f,[1e-3 0.999]);
sprintf('The equilibrium conversion for these feed conditions is: %1.2f',Xeq)

ans =

The equilibrium conversion for these feed conditions is: 0.55

Compute gas phase pressures of each species

Since there is no change in moles for this reaction, we can directly calculation the pressures from the equilibrium conversion and the initial pressure of gases. you can see there is a slightly higher pressure of H2 and CO2 than the reactants, consistent with the equilibrium constant of about 1.44 at 1000K. At a lower temperature there would be a much higher yield of the products. For example, at 550K the equilibrium constant is about 58, and the pressure of H2 is 4.4 atm due to a much higher equilibrium conversion of 0.88.

P_CO = Pa0*(1-Xeq)
P_H2O = Pa0*(1-Xeq)
P_H2 = Pa0*Xeq
P_CO2 = Pa0*Xeq

P_CO =

    2.2748


P_H2O =

    2.2748


P_H2 =

    2.7252


P_CO2 =

    2.7252

Compare the equilibrium constants

We can compare the equilibrium constant from the Gibbs free energy and the one from the ratio of pressures. They should be the same!

K_Temperature
(P_CO2*P_H2)/(P_CO*P_H2O)

K_Temperature =

    1.4352


ans =

    1.4352

Summary

The NIST Webbook provides a plethora of data for computing thermodynamic properties. It is a little tedious to enter it all into Matlab, and a little tricky to use the data to estimate temperature dependent reaction energies. A limitation of the Webbook is that it does not tell you have the thermodynamic properties change with pressure. Luckily, those changes tend to be small.

% categories: Miscellaneous, nonlinear algebra
% tags: thermodynamics, reaction engineering


Published with MATLAB® 7.12

Read and Post Comments

Method of continuity for solving nonlinear equations - Part II

November 02, 2011 at 01:25 PM | categories: nonlinear algebra | View Comments

Method of continuity for solving nonlinear equations - Part II

Method of continuity for solving nonlinear equations - Part II

John Kitchin

clear all; close all; clc

Yesterday in Post 1324 we looked at a way to solve nonlinear equations that takes away some of the burden of initial guess generation. The idea was to reformulate the equations with a new variable $\lambda$, so that at $\lambda=0$ we have a simpler problem we know how to solve, and at $\lambda=1$ we have the original set of equations. Then, we derive a set of ODEs on how the solution changes with $\lambda$, and solve them.

Today we look at a simpler example and explain a little more about what is going on. Consider the equation: $f(x) = x^2 - 5x + 6 = 0$, which has two roots, $x=2$ and $x=3$. We will use the method of continuity to solve this equation to illustrate a few ideas. First, we introduce a new variable $\lambda$ as: $f(x; \lambda) = 0$. For example, we could write $f(x;\lambda) = \lambda x^2 - 5x + 6 = 0$. Now, when $\lambda=0$, we hve the simpler equation $- 5x + 6 = 0$, with the solution $x=6/5$. The question now is, how does $x$ change as $\lambda$ changes? We get that from the total derivative of how $f(x,\lambda)$ changes with $\lambda$. The total derivative is:

$$\frac{df}{d\lambda} = \frac{\partial f}{\partial \lambda} + \frac{\partial f}{\partial x}\frac{\partial x}{\partial \lambda}=0$$

We can calculate two of those quantities: $\frac{\partial f}{\partial \lambda}$ and $\frac{\partial f}{\partial x}$ analytically from our equation and solve for $\frac{\partial x}{\partial \lambda}$ as

$$ \frac{\partial x}{\partial \lambda} = -\frac{\partial f}{\partial \lambda}/\frac{\partial f}{\partial x}$$

That defines an ordinary differential equation that we can solve by integrating from $\lambda=0$ where we know the solution to $\lambda=1$ which is the solution to the real problem. For this problem: $\frac{\partial f}{\partial \lambda}=x^2$ and $\frac{\partial f}{\partial x}=-5 + 2\lambda x$.

dxdL = @(lambda,x) -x^2/(-5 + 2*lambda*x);

[L,x] = ode45(dxdL,[0 1], 6/5);

plot(L,x)
xlabel('\lambda')
ylabel('x')
fprintf('One solution is at x = %1.2f\n', x(end))

One solution is at x = 2.00

We found one solution at x=2. What about the other solution? To get that we have to introduce $\lambda$ into the equations in another way. We could try: $f(x;\lambda) = x^2 + \lambda(-5x + 6)$, but this leads to an ODE that is singular at the initial starting point. Another approach is $f(x;\lambda) = x^2 + 6 + \lambda(-5x)$, but now the solution at $\lambda=0$ is imaginary, and we don't have a way to integrate that! What we can do instead is add and subtract a number like this: $f(x;\lambda) = x^2 - 4 + \lambda(-5x + 6 + 4)$. Now at $\lambda=0$, we have a simple equation with roots at $\pm 2$, and we already know that $x=2$ is a solution. So, we create our ODE on $dx/d\lambda$ with initial condition $x(0) = -2$.

dxdL = @(lambda,x) (5*x-10)/(2*x-5*lambda);

[L,x] = ode45(dxdL,[0 1], -2);
plot(L,x)
xlabel('\lambda')
ylabel('x')
fprintf('One solution is at x = %1.2f\n', x(end))

One solution is at x = 3.00

Now we have the other solution. Note if you choose the other root, $x=2$, you find that 2 is a root, and learn nothing new. You could choose other values to add, e.g., if you chose to add and subtract 16, then you would find that one starting point leads to one root, and the other starting point leads to the other root. This method does not solve all problems associated with nonlinear root solving, namely, how many roots are there, and which one is "best"? But it does give a way to solve an equation where you have no idea what an initial guess should be. You can see, however, that just like you can get different answers from different initial guesses, here you can get different answers by setting up the equations differently.

% categories: Nonlinear algebra
% tags: math


Published with MATLAB® 7.13

Read and Post Comments

Method of continuity for nonlinear equation solving

November 01, 2011 at 07:53 PM | categories: nonlinear algebra | View Comments

method_of_continuity

Contents

Method of continuity for nonlinear equation solving

John Kitchin

Adapted from Perry's Chemical Engineers Handbook, 6th edition 2-63.

function main

We seek the solution to the following nonlinear equations:

$2 + x + y - x^2 + 8 x y + y^3 = 0$

$1 + 2x - 3y + x^2 + xy - y e^x = 0$

f = @(x,y) 2 + x + y - x^2 + 8*x*y + y^3;
g = @(x,y) 1 + 2*x - 3*y + x^2 + x*y - y*exp(x);

In principle this is easy, we simply need some initial guesses. The challenge here is what would you guess? The equations are implicit, so it is not easy to graph them, but lets give it a shot, starting on the x range -5 to 5. The idea is set a value for x, and then solve for y in each equation

x = linspace(-5,5,500);
% at x = 0, we have 2 + y + y^3 = 0
fy = arrayfun(@(x) fzero(@(y) f(x,y), 0.1),x);
gy = arrayfun(@(x) fzero(@(y) g(x,y), 0),x);
plot(x,fy,x,gy)
xlabel('x')
ylabel('y')

graphical solution

you can see there is a solution near x = -1, y = 0, because both functions equal zero there. We can even use that guess with fsolve. It is disappointly easy! But, keep in mind that in 3 or more dimensions, you cannot perform this visualization, and another method could be required.

f1 = @(X) [f(X(1),X(2)); g(X(1),X(2))];
fsolve(f1,[-1 0])

Equation solved at initial point.

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




ans =

    -1     0

We explore a method that bypasses this problem today. The principle is to introduce a new variable, $\lambda$, which will vary from 0 to 1. at $\lambda=0$ we will have a simpler equation, preferrably a linear one, which can be solved, which can be analytically solved. At $\lambda=1$, we have the original equations. Then, we create a system of differential equations that start at this solution, and integrate from $\lambda=0$ to $\lambda=1$, to recover the final solution.

We rewrite the equations as:

$f(x,y) = (2 + x + y) + \lambda(- x^2 + 8 x y + y^3) = 0$

$g(x,y) = (1 + 2x - 3y) + \lambda(x^2 + xy - y e^x) = 0$

Now, at $\lambda=0$ we have the simple linear equations:

$x + y = -2$

$2x + 3y = -1$

x0 = [[1 1]; [2 -3]]\[-2; -1]

x0 =

   -1.4000
   -0.6000

Construct the system of ODEs

We form the system of ODEs by differentiating the new equations with respect to $\lambda$. Why do we do that? The solution, (x,y) will be a function of $\lambda$. From calculus, you can show that:

$\frac{\partial f}{\partial x}\frac{\partial x}{\partial \lambda}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial \lambda}=-\frac{\partial f}{\partial \lambda}$

$\frac{\partial g}{\partial x}\frac{\partial x}{\partial \lambda}+\frac{\partial g}{\partial y}\frac{\partial y}{\partial \lambda}=-\frac{\partial g}{\partial \lambda}$

Now, solve this for $\frac{\partial x}{\partial \lambda}$ and $\frac{\partial y}{\partial \lambda}$. You can use Cramer's rule to solve for these to yield: $\begin{array}{c} \frac{\partial x}{\partial \lambda} = \frac{\partial f/\partial y \partial g/\partial \lambda - \partial f/\partial \lambda \partial g/\partial y}{\partial f/\partial x \partial g/\partial y - \partial f/\partial y \partial g/\partial x }\\ \frac{\partial y}{\partial \lambda} = \frac{\partial f/\partial \lambda \partial g/\partial x - \partial f/\partial x \partial g/\partial \lambda}{\partial f/\partial x \partial g/\partial y - \partial f/\partial y \partial g/\partial x } \end{array}$

For this set of equations: $\begin{array}{l} \partial f/\partial x = 1 - 2\lambda x + 8\lambda y \\ \partial f/\partial y = 1 + 8 \lambda x + 3 \lambda y^2 \\ \partial g/\partial x = 2 + 2 \lambda x + \lambda y - \lambda y e^x\\ \partial g/\partial y = -3 + \lambda x - \lambda e^x \end{array}$

Now, we simply set up those two differential equations on $frac{\partial x}{\partial \lambda}$ and $frac{\partial y}{\partial \lambda}$, with the initial conditions at $\lambda = 0$ which is the solution of the simpler linear equations, and integrate to $\lambda = 1$, which is the final solution of the original equations!

    function dXdlambda = ode(lambda, X)
        x = X(1);
        y = X(2);

        pfpx = 1 - 2*lambda*x + 8*lambda*y;
        pfpy = 1 + 8*lambda*x + 3*lambda*y^2;
        pfplambda = -x^2+8*x*y+y^3;

        pgpx = 2 + 2*lambda*x + lambda*y - lambda*y*exp(x);
        pgpy = -3 + lambda*x - lambda*exp(x);
        pgplambda = x^2 + x*y -y*exp(x);

        dxdlambda = (pfpy*pgplambda - pfplambda*pgpy)/(pfpx*pgpy-pfpy*pgpx);
        dydlambda = (pfplambda*pgpx-pfpx*pgplambda)/(pfpx*pgpy-pfpy*pgpx);

        dXdlambda = [dxdlambda; dydlambda];
    end

[lambda X] = ode45(@ode,[0 1],x0);

the solution

The last X value corresponds to $\lambda = 1$, which is the solution to the original equations.

Xsol = X(end,:);
x = Xsol(1)
y = Xsol(2)

x =

   -1.0004


y =

 -2.9244e-004

evaluate solution

You can see the solution is somewhat approximate; the true solution is x = -1, y = 0. The approximation could be improved by lowering the tolerance on the ODE solver.

f(x,y)
g(x,y)

ans =

  7.1778e-004


ans =

    0.0013

Summary

This is a fair amount of work to get a solution! The idea is to solve a simple problem, and then gradually turn on the hard part by the lambda parameter. What happens if there are multiple solutions? The answer you finally get will depend on your $\lambda=0$ starting point, so it is possible to miss solutions this way. For problems with lots of variables, this would be a good approach if you can identify the easy problem.

end

% categories: Nonlinear algebra


Published with MATLAB® 7.13

Read and Post Comments

Smooth transitions between discontinuous functions

October 30, 2011 at 05:12 PM | categories: miscellaneous, nonlinear algebra | View Comments

transition_discontinuous_functions

Contents

Smooth transitions between discontinuous functions

John Kitchin

In Post 1280 we used a correlation for the Fanning friction factor for turbulent flow in a pipe. For laminar flow (Re < 3000), there is another correlation that is commonly used: $f_F = 16/Re$. Unfortunately, the correlations for laminar flow and turbulent flow have different values at the transition that should occur at Re = 3000. This discontinuity can cause a lot of problems for numerical solvers that rely on derivatives.

Today we examine a strategy for smoothly joining these two functions.

function main

clear all; close all; clc

Friction factor correlations

    function f = fF_laminar(Re)
        f = 16./Re;
    end

    function f = fF_turbulent(Re)
        % Nikuradse correlation for turbulent flow
        function f = internalfunc(Re)
            f = fzero(@(f) 1/sqrt(f) - (4.0*log10(Re*sqrt(f))-0.4), 0.01);
        end
        % this makes the function act as if it were vectorized, so we can
        % get an array of friction factors
        f = arrayfun(@internalfunc,Re);
    end

Plot the correlations

Re1 = linspace(500,3000);
f1 = fF_laminar(Re1);

Re2 = linspace(3000,10000);
f2 = fF_turbulent(Re2);

figure(1)
plot(Re1,f1,Re2,f2)
xlabel('Re')
ylabel('f_F')

You can see the discontinuity at Re = 3000. What we need is a method to join these two functions smoothly. We can do that with a sigmoid function.

Sigmoid functions

A sigmoid function smoothly varies from 0 to 1 according to the equation: $\sigma(x) = \frac{1}{1 + e^{-(x-x0)/\alpha}}$. The transition is centered on $x0$, and $\alpha$ determines the width of the transition.

x = linspace(-4,4);
y = 1./(1+exp(-x/0.1));
figure(2)
plot(x,y)
xlabel('x'); ylabel('y'); title('\sigma(x)')

If we have two functions, $f_1(x)$ and $f_2(x)$ we want to smoothly join, we do it like this: $f(x) = (1-\sigma(x))f_1(x) + \sigma(x)f_2(x)$. There is no formal justification for this form of joining, it is simply a mathematical convenience to get a numerically smooth function. Other functions besides the sigmoid function could also be used, as long as they smoothly transition from 0 to 1, or from 1 to zero.

Smoothly transitioning function to estimate friction factor

    function f = fanning_friction_factor(Re)
        % combined, continuous correlation for the fanning friction factor.
        % the alpha parameter is chosen to provide the desired smoothness.
        % The transition region is about +- 4*alpha. The value 450 was
        % selected to reasonably match the shape of the correlation
        % function provided by Morrison (see last section of this file)
        sigma =  1./(1+exp(-(Re - 3000)/450));
        f = (1-sigma).*fF_laminar(Re) + sigma.*fF_turbulent(Re);
    end

Re = linspace(500,10000);
f = fanning_friction_factor(Re);

figure(1); hold all
plot(Re,f)
xlabel('Re')
ylabel('f_F')
legend 'laminar' 'turbulent' 'smooth transition'

Exiting fzero: aborting search for an interval containing a sign change
    because complex function value encountered during search.
(Function value at -0.0028 is -5.2902-21.627i.)
Check function or try again with a different starting value.

You can see that away from the transition the combined function is practically equivalent to the original two functions. That is because away from the transition the sigmoid function is 0 or 1. Near Re = 3000 is a smooth transition from one curve to the other curve.

Morrison derived a single function for the friction factor correlation over all Re: $f = \frac{0.0076\left(\frac{3170}{Re}\right)^{0.165}}{1 + \left(\frac{3171}{Re}\right)^{7.0}} + \frac{16}{Re}$. Here we show the comparison with the approach used above. The friction factor differs slightly at high Re, becauase Morrison's is based on the Prandlt correlation, while the work here is based on the Nikuradse correlation. They are similar, but not the same.

h = plot(Re,16./Re + (0.0076*(3170./Re).^0.165)./(1 + (3170./Re).^7))

% Append an entry to the legend
[LEGH,OBJH,OUTH,OUTM] = legend;
% Add object with new handle and new legend string to legend
legend([OUTH;h],OUTM{:},'Morrison')

h =

  415.0735

Summary

The approach demonstrated here allows one to smoothly join two discontinuous functions that describe physics in different regimes, and that must transition over some range of data. It should be emphasized that the method has no physical basis, it simply allows one to create a mathematically smooth function, which could be necessary for some optimizers or solvers to work.

end % main

% categories: Miscellaneous, nonlinear algebra
% tags: Fluids


Published with MATLAB® 7.12

Read and Post Comments

compute pipe diameter

October 27, 2011 at 07:26 PM | categories: nonlinear algebra | View Comments

compute pipe diameter

compute pipe diameter

A heat exchanger must handle 2.5 L/s of water through a smooth pipe with length of 100 m. The pressure drop cannot exceed 103 kPa at 25 degC. Compute the minimum pipe diameter required for this application.

Adapted from problem 8.8 in Problem solving in chemical and Biochemical Engineering with Polymath, Excel, and Matlab. page 303.

Contents

Solution

We need to estimate the Fanning friction factor for these conditions so we can estimate the frictional losses that result in a pressure drop for a uniform, circular pipe. The frictional forces are given by $F_f = 2f_F \frac{\Delta L v^2}{D}$, and the corresponding pressure drop is given by $\Delta P = \rho F_f$. In these equations, $\rho$ is the fluid density, $v$ is the fluid velocity, $D$ is the pipe diameter, and $f_F$ is the Fanning friction factor. The average fluid velocity is given by $v = \frac{q}{\pi D^2/4}$.

For laminar flow, we estimate $f_F = 16/Re$, which is a linear equation, and for turbulent flow ($Re > 2100$) we have the implicit equation $\frac{1}{\sqrt{f_F}}=4.0 \log(Re \sqrt{f_F})-0.4$. Of course, we define $Re = \frac{D v\rho}{\mu}$ where $\mu$ is the viscosity of the fluid.

It is known that $\rho(T) = 46.048 + 9.418 T -0.0329 T^2 +4.882\times10^{-5}-2.895\times10^{-8}T^4$ and $\mu = \exp\left({-10.547 + \frac{541.69}{T-144.53}}\right)$ where $\rho$ is in kg/m^3 and $\mu$ is in kg/(m*s).

function main

clear all; close all

u = cmu.units;

T = u.degC(25);
Q = 2.5*u.L/u.s;
deltaP = 103*u.kPa;
deltaL = 100*u.m;

% Note these correlations expect dimensionless T, where the magnitude
% of T is in K
rho = @(T) (46.048 + 9.418*T -0.0329*T^2 +4.882e-5*T^3-2.895e-8*T^4)*u.kg/u.m^3;
mu = @(T) exp(-10.547 + 541.69/(T-144.53))*u.Pa*u.s;

    function fF = fanning_friction_factor(Re)
        if Re < 2100
            error('Flow is probably not turbulent')
        end
        % solve the Nikuradse correlation to get the friction factor
        fz = @(f) 1/sqrt(f) - (4.0*log10(Re*sqrt(f))-0.4);
        fF = fzero(fz, 0.01);
    end

plot fanning friction factor

Re = linspace(2200,9000);
f = arrayfun(@fanning_friction_factor,Re);
plot(Re,f)
xlabel('Re')
ylabel('fanning friction factor')
% You can see why we use 0.01 as an initial guess for solving for the
% Fanning friction factor; it falls in the middle of ranges possible
% for these Re numbers.

Solve the problem

The aim is to find $D$ that solves: $ \Delta p = \rho 2 f_F \frac{\Delta L v^2}{D}$. This is a nonlinear equation in $D$, since D affects the fluid velocity, the Re, and the Fanning friction factor. Here is the function that we need to solve:

    function z = objective(D)
        v = Q/(pi*D^2/4);
        Re = D*v*rho(double(T))/mu(double(T));

        fF = fanning_friction_factor(Re);

        z = deltaP - 2*fF*rho(double(T))*deltaL*v^2/D;
    end

D = fzero(@objective,0.04*u.m)
fprintf('The minimum pipe diameter is %s\n', D.as(u.mm,'%1.2f'))

0.0389653*m
The minimum pipe diameter is 38.97*mm

Any pipe diameter smaller than that value will result in a larger pressure drop at the same volumetric flow rate, or a smaller volumetric flowrate at the same pressure drop. Either way, it will not meet the design specification.

Notes:

If D happened to be large enough to result in laminar flow, this solution would break down, because we didn't define a fanning friction factor correlation for laminar flow. We could do that in our fanning friction factor function, but some care should be taken to ensure it is continuous.

end % main

% categories: Nonlinear algebra
% tags: fluids


Published with MATLAB® 7.13

Read and Post Comments
Next Page »