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

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

Plane Poiseuille flow - BVP

August 11, 2011 at 04:49 PM | categories: odes | View Comments

Plane Poiseuille flow - BVP

Plane Poiseuille flow - BVP

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):

$$\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

we start by defining an initial grid to get the solution on

x = linspace(0,d);

Next, we use the bvpinit command to provide an initial guess using the guess function defined below

solinit = bvpinit(x,@guess);

we also need to define the ode function (see odefun below).

and a boundary condition function (see bcfun below).

now we can get the solution.

sol = bvp5c(@odefun,@bcfun,solinit)

sol = 

    solver: 'bvp5c'
         x: [1x100 double]
         y: [2x100 double]
     idata: [1x1 struct]
     stats: [1x1 struct]

sol is a Matlab struct with fields. the x field contains the x-values of the solution as a row vector, the y field contains the solutions to the odes. Recall that U1 is the only solution we care about, as it is the original u(x) we were interested in.

y = sol.x;
u1 = sol.y(1,:);
u2 = sol.y(2,:);
plot(y,u1); figure(gcf)
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
mu = 1;
Pdrop = -100; % this is DeltaP/Deltax
u = -(Pdrop)*d^2/2/mu*(y/d-(y/d).^2);
hold on
plot(y,u,'r--')

As we would expect, 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 Uinit = guess(x)
% this is an initial guess of the solution. We need to provide a guess of
% what the function values of u1(y) and u2(y) are for a value of y. In some
% cases, this initial guess may affect the ability of the solver to get a
% solution. In this example, we guess that u1 and u2 are constants
U1 = 0.1;
U2 = 0;
Uinit = [U1; U2]; % we return a column vector

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 res = bcfun(U_0,U_d)
% this function defines the boundary conditions for each ode.
% we define the velocity at y=0 and y=d to be zero here. If the plates are
% moving, you would substitute the velocity of each plate for the zeros.
res = [0 - U_0(1);
       0 - U_d(1)];

% categories: ODEs
% tags: math, fluids


Published with MATLAB® 7.10

Read and Post Comments