import math
from numpy import zeros
from scipy.optimize import fsolve
#INPUT
Lmf=2.4; #Length of bed at minimum fluidized condition in m
uo=[2.,4.,6.]; #Superficial gas velocity in m/s
GsII=100.; #Solid circulation rate in kg/m**2 s for Mode II
uoIII=4.; #Superficial gas velocity in m/s for Mode III
GsIII=[42.,50.,100.,200.,400.]; #Solid circulation rate in kg/m**2 s for Mode III
GsIV=[70.,100.,120.]; #Solid circulation rate in kg/m**2 s for Mode IV
dt=0.4; #Column diamter in m
Ht=10.; #Height of column in m
rhos=1000.; #Density of solid in kg/m**3
dpbar=55.; #Particle diameter in micrometer
ephsilonmf=0.5; #Void fraction at minimum fluidization condition
#CALCULATION
#Mode I
ephsilonstar=0.01; #Saturation carrying capacity of gas
ephsilonsd=[0.2,0.16,0.14]; #Solid holdup in lower dense region from Fig.8(b) for various uo
n=len(uo);
i=0;
Hfguess=2.; #Guess value of height
a = zeros(n)
Hf = zeros(n)
Hd = zeros(n)
ephsilonse = zeros(n)
GsI = zeros(n)
# endfunction
while i<n:
a[i]= 3./uo[i]; #Decay consmath.tant
def solver_func(Hf): #Function defined for solving the system
return Lmf*(1-ephsilonmf)-((ephsilonsd[i]-(ephsilonstar+(ephsilonsd[i]-ephsilonstar)*math.exp(-a[i]*Hf)))/a[i])-Ht*ephsilonsd[i]+Hf*(ephsilonsd[i]-ephsilonstar);
Hf[i] = fsolve(solver_func,1E-6)
Hd[i]=Ht-Hf[i];#Height of lower densce region
ephsilonse[i]=ephsilonstar+(ephsilonsd[i]-ephsilonstar)*math.exp(-a[i]*Hf[i]);#Solid holdup at exit
GsI[i]=rhos*uo[i]*ephsilonse[i];#Solid circulation rate from Eqn.(4)
i=i+1;
#Mode II
i=0;
Hfguess2=2;#Guess value of height
ephsilonseII = zeros(n)
HdII = zeros(n)
LmfII = zeros(n)
HfII = zeros(n)
while i<n:
ephsilonseII[i]=GsII/(rhos*uo[i]);#Solid holdup at exit
def solver_func1(Hf): #Function defined for solving the system
return ephsilonseII[i]-ephsilonstar-(ephsilonsd[i]-ephsilonstar)*math.exp(-a[i]*Hf);#From Eqn.(7)
HfII[i] = fsolve(solver_func1,1E-6)
HdII[i]=Ht-HfII[i];#Height of lower dense region
#Length of bed minimum fluidization condtion
LmfII[i]=(1-ephsilonmf)**-1*((ephsilonsd[i]-ephsilonseII[i])/a[i])+Ht*ephsilonsd[i]-HfII[i]*(ephsilonsd[i]-ephsilonstar);
i=i+1;
#Mode III
aIII = 3./uoIII; #Decay consmath.tant
ephsilonsdIII=0.16; #Solid holdup at lower dense region
i=0;
m=len(GsIII);
Hfguess3=2;#Guess value of height
ephsilonseIII = zeros(m)
HdIII = zeros(m)
ephsilonseIII = zeros(m)
LmfIII = zeros(m)
HfIII = zeros(m)
while i<m:
ephsilonseIII[i]=GsIII[i]/(rhos*uoIII);#Solid holdup at exit
def solver_func2(Hf): #Function defined for solving the system
return ephsilonseIII[i]-ephsilonstar-(ephsilonsdIII-ephsilonstar)*math.exp(-aIII*Hf);#From Eqn.(7)
HfIII[i] = fsolve(solver_func2,1E-6)
HdIII[i]=Ht-HfIII[i]; #Height of lower dense region
#Length of bed at minimum fluidization condition
LmfIII[i]=(1-ephsilonmf)**-1*(((ephsilonsdIII-ephsilonseIII[i])/aIII)+Ht*ephsilonsdIII-HfIII[i]*(ephsilonsdIII-ephsilonstar));
i=i+1;
#Mode IV
i=0;
Hfguess4=2;#Guess value of height
aIV = zeros(n)
ephsilonseIV = zeros(n)
HdIV = zeros(n)
HfIV = zeros(n)
LmfIV = zeros(n)
while i<n:
aIV[i]=3./uo[i]; #Decay consmath.tant
ephsilonseIV[i]=GsIV[i]/(rhos*uo[i]); #Solid holdup at exit
def solver_func3(Hf): #Function defined for solving the system
return ephsilonseIV[i]-ephsilonstar-(ephsilonsd[i]-ephsilonstar)*math.exp(-aIV[i]*Hf); #From Eqn.(7)
HfIV[i] = fsolve(solver_func3,1E-6)
HdIV[i]=Ht-HfIV[i];#Height of lower dense region
#Length of bed at minimum fluidization condition
LmfIV[i]=(1-ephsilonmf)**-1*(((ephsilonsd[i]-ephsilonseIV[i])/aIV[i])+Ht*ephsilonsd[i]-HfIV[i]*(ephsilonsd[i]-ephsilonstar));
i=i+1;
#OUTPUT
print 'Mode I';
print '\tuom/s\t\tephsilonse-\tHfm\t\tHdm\t\tGskg/m**2 s';
i=0;
while i<n:
print '\t%f\t%f\t%f\t%f\t%f'%(uo[i],ephsilonse[i],Hf[i],Hd[i],GsI[i]);
i=i+1;
print 'Mode II';
print '\tuom/s\t\tephsilonse-\tHfm\t\tHdm\t\tLmfm';
i=0;
while i<n:
print '\t%f\t%f\t%f\t%f\t%f'%(uo[i],ephsilonseII[i],HfII[i],HdII[i],LmfII[i]);
i=i+1;
print 'Mode III';
print '\tGskg/m** s\tephsilonse-\tHfm\t\tHdm\t\tLmfm';
i=0;
while i<m:
print '\t%f\t%f\t%f\t%f\t%f'%(GsIII[i],ephsilonseIII[i],HfIII[i],HdIII[i],LmfIII[i]);
i=i+1;
print 'Mode IV';
print '\tuom/s\t\tGskg/m**2 s\tephsilonse-\tHfm\t\tLmfm';
i=0;
while i<n:
print '\t%f\t%f\t%f\t%f\t%f'%(uo[i],GsIV[i],ephsilonseIV[i],HfIV[i],LmfIV[i]);
i=i+1;
#====================================END OF PROGRAM ======================================================