کد MATLAB محاسبه حداقل سرعت سیالی شدن

در مکانیک سیالات 2 موضوعی تحت عنوان بسترهای سیالی (Fluidized bed) در سر فصلی با عنوان جریان از روی اجسام غوطه ور وجود دارد که این برنامه مربوط به این قسمت می باشد.
بستر سیالی شده
در یک بستر سیال , با عبور سیال از میان مواد جامد گرانوله سیستم به صورت سوسپانسیون در می اید.و در نتیجه مواد جامد حالت مایع پیدا می کند و این روی کمیات ترمودینامیکی و مکانیکی ان تاثیر می گذارد.بستر های سیالی در صنعت موارد استفاده زیادی دارند مثل راکتورها,حمام های اختلاط با توزیع دما -پوشش های پودری ,خنک کن ها و کوره ها.
مشخصه بستر سیالی افت فشار سیال عبوری از بستر است .وقتی که یک سیال در میان یک توده به جریان در می اید .ابتدا فشار از پایین توده افزایش می یابد تا جایی که نیروی فشار برابر وزن و توده جامد شود و توده به صورت سوسپانسه در اید.افزایش بیشتر شدت جریان باعث می شود که لایه توده شروع به حرکت کند و به شرایط سیال شدن برسد در این شرایط افت فشار تقریبا ثابت می ماند و با افزایش بیشتر شدت جریان تغییر نمی کند.
این برنامه ابتدا ورودی های لازم مانند دانسیته,ویسکوزیته ی سیال که در اینجا اب می باشد و دانسیته و سایز ذرات ,درصد انبساط و کلیه ورودی های لازم را دریافت کرده و با فرمول های مورد استفاده حداقل سرعت سیالی شدن را محاسبه می کند. و همچنین در قسمت دوم برنامه سرعتی که لازم است تا  بستر به اندازه مورد نظر منبسط شود محاسبه شده است.
همچنین صورت مسئله و خود برنامه به صورت pdf جهت دانلود قرار داده شده است.

WWW.CHEMICAL-ENG.IR
%This script calculate Vmin
clear
clc
%dp=-1;Dp=-1;mu=-1;dw=-1;Em=-1;n=-1;le=-1
dp=input('Enter particle density(g/cubic cm):');while ~(~isempty(dp) && (dp > 0)) dp=input('Enter particle density(g/cubic cm):');
end
Dp=input('Enter particle size(mm):');while ~(~isempty(Dp) && (Dp > 0)) Dp=input('Enter particle size(mm):');
end
mu=input('Enter water viscosity(P):');while ~(~isempty(mu) && (mu > 0)) mu=input('Enter water viscosity(P):');
end
dw=input('Enter water density(g/cubic cm):');while ~(~isempty(dw) && (dw > 0)) dw=input('Enter water density(g/cubic cm):');
end
Em=input('Enter Epsilon:');while ~(~isempty(Em) && (Em > 0))Em=input('Enter Epsilon:');
end
n=input('Enter degree that polynomial fit(at least 3):');while ~(~isempty(n) && (n > 0))n=input('Enter degree that polynomial fit(at least 3):');
end
le=input('Enter percent expantion(1-100):');while ~(~isempty(le) && (le > 0))le=input('Enter percent expantion(1-100):');
end
delRho=dp-dw;
g=980;
%This data is obtained from the curve , x-Axis log scale ,y-Axis normal
%scale
x=[0.101024;0.115331;0.135751;0.159786;0.193915;0.240177;0.26866;0.297476;0.343082;0.403826;
۰.۵۰۰۱۶۷;۰.۵۸۸۷۲۴;۰.۶۹۲۹۶;۰.۸۳۲۴۴۲;۱;۱.۲۱۳۵۹;۱.۴۷۲۷۹;۱.۷۶۹۲۴;۲.۱۴۷۱۳;۲.۵۷۹۳۲;۳.۱۳۰۲۲;۳.۹۱۶۷;
۴.۷۵۳۲۶;۵.۷۶۸۴۹;۶.۹۲۹۶;۸.۸۴۹۱۷;۱۰.۶۳۰۴;۱۲.۷۷۰۱;۱۵.۳۴۰۵;۱۸.۸۰۷۷;۲۳.۰۵۸۵;۲۸.۸۵۲۱;۳۶.۴۷۰۹;۴۳.۸۱۲;
۵۳.۷۱۴۱;۶۵.۸۵۴۳;۸۱.۵۶۵۱;۱۰۰;۱۲۶.۴۰۷;۱۵۳.۴۰۵;۱۹۹.۹۳۳;۲۴۰.۱۷۷;۲۸۵.۵۹۶;۳۴۶.۵۹۶;۴۱۲.۱۴;۵۲۰.۹۷۲;
۶۰۰.۸۴۳;۷۳۶.۶۴۲;۹۸۹.۸۶۳];
y=[4.49408;4.51183;4.51183;4.49408;4.50296;4.45858;4.46746;4.4497;4.43195;4.39645;4.36982;4.32544;
۴.۲۸۹۹۴;۴.۲۶۳۳۱;۴.۲۱۸۹۳;۴.۱۶۵۶۸;۴.۱۱۲۴۳;۴.۰۵۰۳;۳.۹۹۷۰۴;۳.۹۴۳۷۹;۳.۸۸۱۶۶;۳.۸۰۱۷۸;۳.۷۳۹۶۴;۳.۶۸۶۳۹;
۳.۶۰۶۵۱;۳.۵۳۵۵;۳.۴۷۳۳۷;۳.۴۲۰۱۲;۳.۳۶۶۸۶;۳.۳۱۳۶۱;۳.۲۵۱۴۸;۳.۱۸۰۴۷;۳.۱۳۶۰۹;۳.۰۹۱۷۲;۳.۰۲۹۵۹;۲.۹۶۷۴۶;
۲.۹۲۳۰۸;۲.۸۶۹۸۲;۲.۸۳۴۳۲;۲.۷۸۹۹۴;۲.۷۴۵۵۶;۲.۷۱۰۰۶;۲.۶۸۳۴۳;۲.۶۵۶۸;۲.۶۳۰۱۸;۲.۵۷۶۹۲;۲.۵۶۸۰۵;۲.۵۳۲۵۴;۲.۵۰۵۹۲];
logx=log10(x);
C=[1.75*dw/(Dp/10*Em^3),150*mu*(1-Em)/(((Dp/10)^2)*Em^3),-g*delRho];
P=polyfit(logx,y,n);
z=-1:0.1:3;
p=polyval(P,z);
plot(z,p,logx,y,'r')
legend('polyfit curve','data curve')
v=roots(C);
figure
semilogx(x,y)
grid on
if v(1,1)>=0
Vm=v(1,1)
else
Vm=v(2,1)
end
NRe=Dp/10*Vm*dp/mu
m=polyval(P,log10(NRe))
E=1-(1-Em)/(1+le/100)
V0=Vm*10*(E/Em)^3.9;
sprintf('V0=%f mm/s',V0)