newtonfluidflow.m

matlab在化学工程中的应用实例程序,对大家研究学习有指导作用。

function NewtonFluidFlow
% 水平管中不可压缩牛顿流体的层流流动
% 求解动量守衡方程以获得剪切力和速度分布
% Laminar Flow in a Horizontal Pipe (Newtonian Fluid)
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/05/30 $

clear all
clc

global deltaP L mu R
deltaP = 500;
L = 10;
mu = 8.937e-4;
R = 0.009295;
a = 0;
b = R;

% 求解ODE-BVP问题
% initialize of solution with a guess of y1(r)=0,y2(r)=0,y3(r)=0  
solinit = bvpinit(linspace(a,b,100),[0 0 0]);
sol = bvp4c(@ODEs,@BCfun,solinit);

% 分析解
r = sol.x;
TauAnal = (deltaP/(2*L))*r;
uAnal = (deltaP*R^2/(4*mu*L))*(1-(r/R).^2);
umAnal = deltaP*R^2/(8*mu*L);

% 绘制图
% 剪应力分布
tau = sol.y(1,2:end)./sol.x(2:end);
tau = [0 tau];
plot(r,tau,'b-',r,TauAnal,'r-.')
xlabel('管径位置r，m')
ylabel('剪应力，kg/(m s^2)')
legend('数值解','分析解')
figure
% 剪切速率分布
plot(r,sol.y(2,:),'b-',r,uAnal,'r-.')
xlabel('管径位置r，m')
ylabel('剪切速率，m/s')
legend('数值解','分析解')

% 计算催化剂的有效因子
fprintf('\t平均剪切速率为: um = %.4f',sol.y(3,end))

% ------------------------------------------------------------------
function dydr = ODEs(r,y)
global deltaP L mu R
rTau = y(1);
u = y(2);
um = y(3);
drTaudr = deltaP*r/L; 
if r > 0             % 避免除数为零
    Tau = rTau/r;
else
    Tau = 0;
end
dudr = -Tau/mu;
dumdr = u*2*r/R^2;

dydr = [drTaudr; dudr; dumdr];

% ------------------------------------------------------------------
function bc = BCfun(ya,yb)
bc = [ya(1); yb(2); ya(3)];