computer programma's ( build 2016,update 07/07/2016)
1e) Compartmental distillation model.
info :staat tussen de code zoals weergegeven op de startpagina.
numerical method : Euler Integration (oplossing van system of differential equations)
dowload code:
programma compartmoddist
output view : (video,*.png)
view compartmoddist
sub: liquid enthalpy
info:staat tussen de code zoals weergegeven op de startpagina.
polynoom: enthalpy_liquid (T) = a(1)+a(2)*T+a(3)*T^2
download code:
liquidenthalpy
output view:
view liquidenthalpy
sub: vapour enthalpy
info:staat tussen de code zoals weergegeven op de startpagina.
polynoom:enthalpy_vap(T)=a(1)+a(2)*T+a(3)*T^2+a(4)*T^3+a(5)*T^4
download code:
vapourenthalpy
output view:
view vapourenthalpy
plot : polynoom op de canvas (= deel van het computer scherm)
download code : polynoom canvas
output view : view polynoom canvas
2e) Molecular weight and density of tray liquid (binary)
info:staat tussen de code zoals weergegeven op de startpagina.
download code:
molecular weight and density
output view:
view molecular weight and density
3e) Tray hydraulics (Francis weir formula)
info:staat tussen de code zoals weergegeven op de startpagina.
download code:
tray hydraulics
output view:
view tray hydraulic
4e) Equilibrium vapour composition and tray temperature
info: Dit programma kan eveneens gebruikt worden voor het
kookpunt van niet ideale vloeistoffen mits aanpassing van
de code. (binary system of liquids)
download code: bubblepointnonideal
output view : 1e) view optie1
2e) view optie2
5e)Stripping column : version 1.0 without graphics.
info:staat tussen de code zoals weergegeven op de startpagina.
Doel: Het programma vindt het aantal schotels in de stripping kolom.
download code: stripping column no
output view : view stripping column no
6e)Stripping column : version 1.0 with graphics.
info:staat tussen de code zoals weergegeven op de startpagina.
Doel: Het programma vindt het aantal schotels in de stripping kolom.
download code : stripping column yes
output view : view stripping column yes
7e) Thiele-Geddes method (n-butane,i-pentane,n-pentane,n-hexane).
info:staat tussen de code zoals weergegeven op de startpagina.
info:Thiele-Geddes op de extrapagina.
ref :Design of equilibrium stage processes,Budford D.Smith,1963,McGraw-Hill.
download code: thiele-geddes method
output view : 1e) optie view n-butane
2e) optie view i-pentane
3e) optie view n-pentane
4e) optie view n-hexane
8e) McCabe-Thiele method (binary solution ,full graphics ).
info:staat tussen de code zoals weergegeven op de startpagina.
download code: McCabe-Thiele method
output view : 1e) first input
2e) graphics output
9e) Draw tielines from x-y diagram in the enthalpy diagram (example).
info:staat tussen de code zoals weergegeven op de startpagina.
downloadcode: tielinesnew14
output view: graphic tielines
10e) Algo: TDMA =TriDiagonal Matrix Algorithm.
info:staat tussen de code zoals weergegeven op de startpagina.
downloadcode: TDMA code
downloadcode: TDMA sub code
output view: view TDMA
11e) One-Dimensional Steady Heat Conduction (use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
problem: 1-dimensional steady state conduction without heat generation
in a rectangular slab. (k=thermal conductivity)
equation : k*d/dx(dT/dx)=0
boundary condition (=bc) : left side temperature of the slab = 100 °C
right side temperature of the slab = 200 °C
downloadcode: 1Dheatsteady
output view : result1Dheatsteady
12e) One-Dimensional Convection-Diffusion Steady (use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
problem:1-dimensional steady state convection-diffusion of component C is
transported through convection and diffusion. ( u =speed,D=diffusion coef.)
equation : -u*dC/dx+D*d/dx(dC/dx) = 0 (FD,central differences both)
boundary condition (=bc): x=0 concentration at this point c0 = 0
x=L concentration at this point cL = 1
downloadcode: 1Dconvectiondiffusion
output view: result 1Dconvectiondiffusion
13e ) Tubular reactor with axial dispersion 'first order kinetics' (use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
problem:1-dimensional steady state tubular reactor (u=mean axial velocity,
D=axial dispersion coefficient,k=reaction rate )
equation : D*d/dx(dCA/dx)-u*dCA/dx-k*CA=0 ( FD,central,upwind)
boundary condition (=bc): proposed by Danckwerts ( BVP problem)
x=0 inlet u*CA_in=u*CA-D*dCA/dx ( Mixed bc )
x=L outlet dCA/dx=0
downloadcode: 1DTubularreactorfirstorderkinetics
output view : profiletubularreactorfirstorder
14e ) Tubular reactor with axial dispersion 'second order kinetics' (use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
problem:1-dimensional steady state tubular reactor (u=mean axial velocity,
D=axial dispersion coefficient,k=reaction rate )
equation : D*d/dx(dCA/dx)-u*dCA/dx-k*CA^2=0 ( FD,central,upwind)
=> main diagonal is not constant ( do loop for ' iteration ' )
result : output view
boundary condition (=bc): proposed by Danckwerts ( BVP problem)
x=0 inlet u*CA_in=u*CA-D*dCA/dx (Mixed bc)
x=L outlet dCA/dx=0
downloadcode: 1DTubularreactorsecondorderkinetics
output view : profiletubularreactorsecondorder
15e ) 1D :Transient ( dT/dt <> 0 ) heat Conduction in a reactangular slab : T=T(t,x) (use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
problem:1-dimensional transient heat conduction in a reactangular slab .
alpha = k/(rho* Cp). ( t=time,x=space,T=temperature => T=T(t,x))
M° = M = alpha*dt/dx^2
equation : dT/dt = alpha*d/dx(dT/dx) ( FD,forward time,central in space)
initial condition (= ic): T=T(t=0, x ; given Temperature)
boundary condition (=bc): ( BVP problem)
x=0 left T0= temperature (fixed condition)
x=L right TL= temperature
solve : implicit discretization method ( unconditionally stable )
T(n+1,i) -T(n,i) = M*(T(n+1,i+1)-2*T(n+1,i)+T(n+1,i-1))
info : explicit discretization is stable for the below condition
dt and dx is not indepent in the above condition M° =< 1/2.
T(n+1,i) = T(n,i)+M*(T(n,i+1)-2*T(n,i)+T(n,i-1))
single time : dt
downloadcode: heatconductionimplicittimeinonestepfinalcontrolfinal
output view : viewheatconductionimplicittimeinonestepfinalcontrolfinal
time interval : number of time steps (r-1)*dt
downloadcode: heatconductionimplicittimeintervalfinalfinalfinalfinal
output view: viewheatconductionimplicittimeintervalfinalfinalfinalfinal
programming can be used for a brick wall with some thickness.
16e ) 1D :Transient ( dT/dt <> 0 ) heat Conduction in a cylinder : T=T(t,r) (use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
problem:1-dimensional transient heat conduction in a cylinder .
alpha = k/(rho* Cp). ( t=time,r=space,T=temperature => T=T(t,r))
M = alpha*dt/dr^2
equation : dT/dt = alpha*(d/dr(dT/dr)+1/r*dT/dr) ( FD,forward time,central in space(both))
h=1/r*dT/dr (problem when r=0 in the center of the cylinder)
from calculus solve this with the L'Hopital's rule (0/0) at node 1.
result : for node 1 is dT/dt=alpha*(d/dr(dT/dr)+d/dr(dT/dr) = alpha*2*(d/dr(dT/dr))
initial condition (= ic): T=T(t=0, r ; given Temperature inside the cylinder) = ( T(n,i) see below why )
boundary condition (=bc): ( BVP problem)
r=0 center T0= temperature ? (see below)
r=m*dr rand TL= temperature = ( T_surface ,see below) ,( m ? see below :example)
solve: implicit discretization method
T(n+1,i) -T(n,i) = M*[(T(n+1,i+1)-2*T(n+1,i)+T(n+1,i-1)) +1/(2*m)*(T(n+1,i+1)-T(n+1,i-1)) (***)
node 1 :1e) T(n+1,i) -T(n,i) = 2*M*(T(n+1,i+1)-2*T(n+1,i)+T(n+1,i-1)) (*)
We put i=1 in the term T(n+1,i-1) => T(n+1,0) 'hypothetical node 0'
2e) r=0 , no flux => dT/dr=0 (central in space) [ ref node : i=1 ( use node = 0 and node = 2) ]
discretization version of dT/dr=0 :( T(n+1,2) - T(n+1,0) )/2*dr =0
result 2e : T(n+1,0) = T(n+1,2) (**)
result : Put in i=1 (*) and then use (**) to eliminate T(n+1,0) we get the first equation from the matrix.
(-4*M-1)*T(n+1,1)+4*M*T(n+1,2) = - T(n,1) *
node j : ( j <> 1 and j <> node last) => m=j-1 ( tridigonal system to use TDMA a,b,c) ' M= stability <> m '
a*T(n+1,i-1)+b*T(n+1,i)+c*T(n+1,i+1) = d use (***) to find a,b,c,d ( TDMA general equation = TDMA_G)
a = (M-M/(2*m))
b = -(2*M+1)
c = (M+M/(2*m))
d = - T(n,i) ( initial values ; given temperature)
node last : How to find last node,value for m ? radius cylinder = r and dr ( m_last = integer = r/dr - 1)
example : given cylinder with radius =0.4 cm ,dr=0.1 cm => m_last = 0.4/0.1-1=3
node 1 : no m in (*)
node j : m=1,2,..m_last -1
node last : m=m_last
number of equations = no = m_last +1
How to find a,b,c,d for last equation (a*T(n+1,i-1)+b*T(n+1,i)+c*T(n+1,i+1)= d) ?
example : i_last = i=4 = r/dr , m_last =3
a=(M-M/(2*m_last)) => a=(M-M/6)=5*M/6
b=-(2*M+1)=-2*M-1=> b=-2*M-1
c= to the other side of the equation in d
T(n+1,i+1) = Temperature on the surface of the cylinder = T_surface
d = -T(n,i) - (M+M/(2*m_last))*T_surface => d= - T(n,i)-(M+M/6)*T_surface
d= -T(n,i)-7M/6*T_surface
The values of tridiagonal matrix depent on M and m.
single time : dt
downloadcode : heatconductionimplicitsingletimecylinder
output view : viewheatconductionimplicitsingletimecylinder
time interval : t=tbegin,t=t+dt..tend
downloadcode : heatconductionimplicitintervaltimecylinder
output view : viewheatconductionimplicitintervaltimecylinder
17e ) 1D :Transient ( dC/dt <> 0 ) diffusion in sphere : C=C(t,r) (use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
problem:1-dimensional transient diffusion in sphere .
D = diffusion coefficient ( t=time,r=space,C=concentration => C=C(t,r))
equation : dC/dt = D*(d/dr(dC/dr+2/r*dC/dr) * ( FD,forward time = FT,central in space =CS)
initial condition (= ic): C=C(t=0, r ; given concentration) = Co , 0 =< r =<R
boundary condition (=bc): ( BVP problem)
r=0 inside sphere dC/dr = 0
discretization , c(n,-1) and c(n+1,-1) are ghost nodes
r=0 => C=concentration ( c(n,1) - c(n,-1))/2*dr =0 => c(n,-1) = c(n,+1) or c(n+1,-1) = c(n+1,+1) **
r=R on (s)urface of the sphere C=Cs = 0
M = D*dt/dr^2
solve : implicit discretization method
node i : i=1,2,.......
discretization version
C(n+1,i) -C(n,i) = M*(C(n+1,i+1)-2*C(n+1,i)+C(n+1,i-1) + 1/m*(C(n+1,i+1) - C(n+1,i-1))) ***
q1=1+1/m
q2=1-1/m
use tdma ?
-q2*C(n+1,i-1)+(1/M+2)*C(n+1,i) - q1*C(n+1,i+1) = 1/M*C(n,i)
initial concentration : mass/volume = C(n,i) =Co
-q2*C(n+1,i-1)+(1/M+2)*C(n+1,i) - q1*C(n+1,i+1) =1/M*Co
a = -q2
b = (1+M/2)
c = -q1
d= 1/M*Co
a*C(n+1,i-1)+b*C(n+1,i)+c*C(n+1,i+1) = d
node i=0 : use L'Hospital's rule (0/0) " first equation of matrix "
dC/dt= D*(d/dr(dC/dr)+2*d/dr(dC/dr)) = 3*D*d/dr(dC/dr)
discretization version (use **)
C(n+1,0) - C(n,0) =3*M *( C(n+1,1) - 2*C(n+1,0) + C(n+1,-1))
C(n+1,0) - C(n,0) =3*M *( C(n+1,1) - 2*C(n+1,0) + C(n+1,1)) = 6*M*(C(n+1,1) - C(n+1,0))
1/(6*M)*(C(n+1,0)-C(n,0))=C(n+1,1) - C(n+1,0)
(1+1/(6*M)) *C(n+1,0) - C(n+1,1) = 1/ (6*M)*C(n,0)
intial concentration : mass/volume = C(n,0) =Co
(1+1/(6*M)) *C(n+1,0) - C(n+1,1) = 1/(6*M) *Co
b = (1+1/(6*M))
c = -1
d = 1/(6*M)*Co
example : Determination concentration of drug in centre of a bead (=parel ,Co = mass/volume (sphere) )
1e) parel = gel matrix + ( drug content = Co )
2e) D of drug in the gel matrix
result : When the drug is on surface then concentration is zero
( takes infinite time for all drug content )
The values of tridiagonal matrix depent on M and m.
single time : dt
downloadcode: TRANSIENTDIFFUSIONIMPLICITSINGLETIMESPHEREFINAL
output view : viewTRANSIENTDIFFUSIONIMPLICITSINGLETIMESPHEREFINAL
18e) 1D-Parabolic equation ( Heat conduction,diffusion equation...) ,non-dimensional.
info:staat tussen de code zoals weergegeven op de startpagina.
method explicit : the parabolic equation ( see ' 15e)' )
time interval : t=0, dt.....t_end
downloadcode : parabolicequationexplicitnondimensional
output view :viewparabolicequationexplicitnondimensional
problem : general computational mathematics for pde = (p)artial (d)ifferential (e)quation .
u=u(x,t)
discretization : ( use of gridpoints),general
uij=u(i,j) : C(Center) 0 -4
x = i*dx ( i = 0,1,2.....)
t = j*dt ( j = 0,1,2.....)
neighborhood : N (North) : u(i,j+1) dy 1
S (South) : u(i,j-1) -dy 1
W(West) : u(i-1,j) -dx 1
E(East ) : u(i+1,j) dx 1
pde : du/dt=d/dx(du/dx) ' 1D Parabolic equation '
discretization : (use of gridpoints,general (f)inite-(d)ifference (= FD) approximation)
(u(i,j+1)-u(i,j))/dt =1/dx^2*[(p*(u(i+1,j+1)-2*u(i,j+1)+u(i-1,j+1))+(1-p)*(u(i+1,j)-2*u(i,j)+u(i-1,j))]
domain : 0 = < p =<1
p = 0 explicit method
p = 1/2 Crank-Nicolson
p = 1 fully implicit backward time-difference method.
19e) 2D-Laplace equation (steady = dT/dt = 0).
info:staat tussen de code zoals weergegeven op de startpagina.
example : How to make the calculation for the 'temperature distribution'
on a square plate ( 2D). (T = temperature )
method : iteration method(without source)=2D Laplace equation
equation : pde d/dx(dT/dx)+d/dy(dT/dy) = 0
initial condition (=ic) : temperature inside = 0 °C
boundary condition (=bc) : temperature:top tp
bottom tb
left tl
right tr
downloadcode : laplacechem2dfinal
output view : viewsteadylaplace2Diteration
method : iteration method (q = source term) = 2D Poisson equation
equation : pde d/dx(dT/dx)+d/dy(dT/dy) = -q
initial condition (=ic) : temperature inside = 0 °C
boundary condition (=bc) : temperature:top tp
bottom tb
left tl
right tr
downloadcode : laplacechem2dfinalsource
output view : viewsteadylaplace2Diterationsource
method : iteration methods (sor : relaxation method)
sor = 0 'Jacobi method'
sor = 1 'Gauss siedel method'
sor < 1 'underrelaxation method'
sor > 1 'overrelaxation method'
equation : pde d/dx(dT/dx)+d/dy(dT/dy) = 0
initial condition (=ic) : temperature inside = 0 °C
boundary condition (=bc) : temperature:top tp
bottom tb
left tl
right tr
downloadcode : laplacechem2dfinaldifferentmethod
output view : viewsteadylaplace2Diteration
new: 08/11/2017: Laplace equation 2D (numerical scheme),sor factor
Dirichlet Boundary condition ( = fixed value on the boundary),corrected 'version 1.1'
print "---------------------------------------------------------------------------"
print " sor = 0 , jacobi; sor = 1 Gauss Seidel "
print " sor < 1 ;underrelaxation,sor > 1 ;overrelaxation "
print "----------------------------------------------------------------------------"
download (TrueBasic,code) :laplacechem2dfinaldifferentmethodsfinalversion11
download (TrueBasic,outp) :laplacechem2dfinaldifferentmethodsfinalversion11outp
20e) 2D-Laplace equation ,steady (dT/dt =0),( ADI method use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
example : How to make the calculation for the 'temperature distribution'
on a square plate ( 2D). (T = temperature )
method : ADI method(without source)=2D Laplace equation
ADI method : sweeping method : j sweep (horizontal)
T(i-1,j) -4*T(i,j)+T(i+1,j) = -T(i,j-1) -T(i,j+1)
j -> i=1..div-1 => (TDMA) ' right' (i-1,i,i+1)
i sweep (vertical )
T(i,j-1) -4*T(i,j)+T(i,j+1) = -T(i-1,j) -T(i+1,j)
i -> j=1..div-1 => (TDMA) 'right'(j-1,j,j+1)
aa *x =d (aa ,tridiagonal system 'a,b,c')
'a,b,c' => (i-1,i,i+1)
'a,b,c' => (j-1,j,j+1)
ref : seventh edition (Wiley)
advanced engineering mathematics
Erwin Kreyszig
equation : pde d/dx(dT/dx)+d/dy(dT/dy) = 0
initial condition (=ic) : temperature inside =to
boundary condition (=bc) : temperature:top tt
bottom bt
left lt
right rt
downloadcode : laplace2DsteadyADIfinal
output view :viewlaplace2DsteadyADIfinal
21e) 2D-Laplace equation ,unsteady (dT/dt <>0),( ADI method one step to cvg use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
example : How to make the calculation for the 'temperature distribution'
on a plate ( 30 by 40).
equation : pde dT/dt=d/dx(dT/dx)+d/dy(dT/dy)
initial condition (=ic) : temperature inside =to
boundary condition (=bc) : given on the boundary
method : ADI method(without source)=2D Laplace equation
ADI method : sweeping method :x sweep (horizontal)
( calculation : TDMA)
y sweep (vertical ) use values from xsweep
( calculation : TDMA)
(without use ADI method : in general we get ' pentadiagonal matrix ' )
pentadiagonal matrix structure : ' s1 0 TDMA 0 s2 ' , s1 and s2 are off diagonal from matrix.
, 0 = diagonal with only zero elements.
, TDMA = ' a b c '
opmerking : general solving 2D problems use 'Finite Element Method'
one step = ' xsweep and ysweep '
downloadcode : LAPLACE2DUNSTEADYADIonestep
output view :viewLAPLACE2DUNSTEADYADIonestep
update ?.................... (one step to cvg ) 'finite differences'
downloadcode : LAPLACE2DUNSTEADYADICVGFINAL
output view :viewLAPLACE2DUNSTEADYADICVGFINAL
22e) 2D-Heat equation ,unsteady (dT/dt <>0),( ADI simulation use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
example : How to make the calculation for the 'temperature distribution'
on a square plate. (given BC,IC)
equation : pde dT/dt=alpha*(d/dx(dT/dx)+d/dy(dT/dy) )
initial condition (=ic) : temperature inside the square =to
boundary condition (=bc) : given on the boundary
top : fixed temperature
bottom : dT/dy =0
left : dT/dx =0
right : fixed temperature
method : ADI method( heat 2D equation )
uniform grid : dx=dy
CFL = Courant Fredric Levi number = M =m= alpha*dt/dx^2
m1 = 2+2/m ,m2 = 2/m - 2
T = T(i,j,time)
ADI method : first step : timestep ( n => n+1/2 ='*' ) " x sweep "
-T(i-1,j,*)+m1*T(i,j,*)-T(i+1,j,*) = T(i,j+1,n)+T(i,j-1,n)+m2*T(i,j,n) (1)
TDMA = (a,b,c) ='(i-1,i,i+1)'
a = -1
b = m1
c = -1
( for boundary nodes use , ghost nodes => go to 17e) ' see code = downloadcode ' )
second step :timestep ( n+1/2 => n +1='+ ' ) " y sweep "
-T(i,j-1,+)+m1*T(i,j,+)-T(i,j+1,+) = T(i+1,j,*)+T(i-1,j,*)+m2*T(i,j,*) (2)
TDMA = (a,b,c)='(j-1,j,j+1)'
a = -1
b = m1
c = -1
( for boundary nodes use , ghost nodes=>go to 17e) ' see code = downloadcode ')
" time for simulation " = 0,dt..t_max
downloadcode : LAPLACE2DUNSTEADYADITIMEINTERVALfinal
output view : viewLAPLACE2DUNSTEADYADITIMEINTERVALfinal
23e) intro example of finite elements,( use TDMA)
info:staat tussen de code zoals weergegeven op de startpagina.
example : How to make the calculation for the 'temperature distribution'
in a composite wall. ' chemical plant construction'
equation : q = -k*A*dT/dx
q = heat flow rate
k = conduction coef. ,k(1),k(2),....
T = temperature
x = thickness , dx(1),dx(2),....
boundary condition (=bc) : given on the boundary
outside composite wall = T0 (°C) ( x=0)
first layer : k(1),dx(1)
second layer : 'insulation',k(2),dx(2)
thridth layer : k(3),dx(3) =dx(1)
outside composite wall = TL (°C) (x=L)
L= dx(1)+dx(2)+dx(3)
insulation principles : k(2) << k(1) = k(3)
principles : 'simple finite elements'
(a) Pre-processing : discretization of selected elements.
(b) Calculation of element matrices,vectors.
(c) Assembly of the element matrices (example = 3) and vectors (global stiffness matrix).
(d) Incorporation of the boundary conditions into the global equations.
(e) Solutions of the equations (use TDMA,other),unknown nodal values (T2,T3)
(f) Post-processing : heat flows. ( = > temperature distribution ).
downloadcode : introfiniteelementfinal
output view : viewintrofiniteelementfinal
24e) finite elements 'shape functions "
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Fundamentals of the finite element method for heat and fluid flow
Roland W.Lewis,Perumal Nithiarasu,Kankanhalli N.Seetharamu
Wiley .
example : How to make the calculation on a square plate (a*b=a) and
find the temperature at some location (xo,yo) = global.
T1 = T1(0,0)
T2 = T2(a,0)
T3 = T3(a,b)
T4 = T4(0,b)
T=a1+a2*x+a3*y+a4*x*y
unknown : a1,a2,a3,a4
to find ' T ' use T1,T2,T3,T4.
N1,N2,N3,N4 ' shape function ' (inside the code )'
convert (xo,yo) to (xi,yi) =local ' (inside the code )'
T=N1*T1+N2*T2+N3*T3+N4*T4
dT/dx = dN1/dx*T1+dN2/dx*T2+dN3/dx*T3+dN4/dx*T4
dT/dy = dN1/dy*T1+dN2/dy*T2+dN3/dy*T3+dN4/dy*T4
! Calculation of the heat fluxes.
qx = -kx*dT/dx
qy = -ky*dT/dy
equation : q = -k*dT/dx | (xi,yi)
q = heat flow rate,qx,qy
k = conduction coef,kx,ky
T = temperature
x ,y = direction (space)
downloadcode : heatsquareheatflux_x_y
output view : viewheatsquareheatflux_x_y
25e) (1D Laplace equation without source )finite element formulation using galerkin method (use TDMA).
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Petrov–Galerkin method
solution : a) Pre-processing :Discretize the domain (between x=0 and x=L).
b) Weak form ( PI) 'Galerkin's weighted residual method ' (use='Ni,Nj')
PI = Integration by parts ( result *)
weighted function = 'shape function,see above programming'
c) T (element ) = NiTi + NjTj (see also: programming 23e))
example : 5 elements ( N1N2,N2N3,N3N4,N5N6)
i=1,2,..5 ,j=1,2,..6
d) The characteristics for each element (e)= 'NiNj'
example : [k(e)] { phi(e)} ={f(e)} ,{f(e)} = column vector (2,1) = Load vector
,[k(e)] = matrix (2,2)
, {phi(e)} = column vector (2,1)
m = k*A/L ( 1D : A = 'm^2'=1)
matrix ( 2,2) =
ke(1,1) = 1*m,ke(1,2) = -1*m,ke(2,1) = -1*m,ke(2,2)=1*m
f(e):vector (2,1) =
f(1,1) = -k*dT/dx | xi ( begin of element =xi :coordinate,i=1,2..)
f(1,2) = k*dT/dx | xj ( end of element =xj:coordinate,j=1,2..)
phi(e):vector(2,1) =
phi(1,1) = Ti
phj(1,2) = Tj
e) We assemble all elements and recast the matrix in a global system of equations =
Stiffness matrix.
f) We incorporate the boundary condition (=bc) in the global system.
g) We choose a solver (= TDMA,other inv(A)*b=x,when Ax=b)
h) Post-processing : heat flows slab ( = > temperature distribution ) and heat transfer.
heat transfer : qx| x=0 = -k*dT/dx| x=0 ?
qx | x=L= -K*dT/dx | x=L ?
example : 1D Steady State conduction in plane slab.
L = 10 cm,k = 1.2 W/cm °C,A=1 'm^2',elements = 5
equation : d/dx(k*dT/dx) +S=0
S = 0 ( source term)
boundary condition (=bc) : left (x=0) side temperature of the slab = T0 °C
right (x=L) side temperature of the slab = TL °C
k= inside 1D laplace equation
find : * dT/dx |(x=0) and dT/dx | (x=L) 'unknown'
downloadcode: GALERKIN1DSTEADYSTATECONDUCTIONS0FINAL (version 1.0)
downloadcode: GALERKIN1DSTEADYSTATECONDUCTIONS0FINALsubfinal (version 1.1,sub)
output view : viewGALERKIN1DSTEADYSTATECONDUCTIONS0FINAL
26e) (1D Laplace equation with source )finite element formulation using galerkin method (use TDMA).
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Fundamentals of the finite element method for heat and fluid flow
Roland W.Lewis,Perumal Nithiarasu,Kankanhalli N.Seetharamu
Wiley .
( nuclear reactors,chemical and biological reactors)
example : 1D Steady State conduction in a plane wall with internal heat source (s).
L *=30 mm,k = 21 W/m °C,A=1 'm^2',elements = 4,dx=7.5mm
equation : k* d/dx(dT/dx) +S=0
S=0.3 MW/m^3 "internal heat source "
principal of the loadvector ( 'see inside the code' )
boundary condition (=bc) : * symmetric = -L,L (center x=0 to the wall x=L 'TL=Twall')
downloadcode: GALERKIN1Dconductionsourcetfinal
test :matrixinv
output view :viewGALERKIN1Dconductionsourcetfinal
test view :viewmatrixinv
27e) (1D heat generation ,solid cylinder)finite element formulation using galerkin method (use TDMA).
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Fundamentals of the finite element method for heat and fluid flow
Roland W.Lewis,Perumal Nithiarasu,Kankanhalli N.Seetharamu
Wiley .
( nuclear reactors,chemical and biological reactors)
example : 1D Steady State conduction in a cylinder with internal heat source (G).
equation : k* d/dr(dT/dr+1/r*dT/dr) +G=0
G "internal heat source "
boundary condition (=bc): r=ro ,T=Twall cylinder=Tw
r=ro ,-k*dT/dr=ho*(Tw-Ta)
principal of the loadvector ( 'see inside the code' )
principal :stiffness matrix ('see inside the code')
downloadcode: GALERKIN1Dconductioncylindersourcefinalsub
output view : viewGALERKIN1Dconductioncylindersourcefinalsub
28e) (2D laplace equation ,galerkin method:global stiffness matrix ) one element.
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Fundamentals of the finite element method .
Hartley Gradin,Jr. (pg:333-350)
Macmillan Publishing Company
(Heat Transfer and fluid flow in 2 Dimensions ,'Steady-State Heat flow)
principal of local nodes ( 'see inside the code')
(three-node triangle ) '1,2,3'
principal of global nodes ('see inside the code')
(global nodes ) '1,3,4' (example side : local 1,2 =inside => no convection)
example : triangle (b) = ' 29e) '
pt 1: (0,0),pt 2 (2.5,2),pt 3 (0,4)
side 1 to 2 : local = inside the structure
side 2 to 3 : local = convection Ta,h
side 3 to 1: local = insulation
result : globalstiffness matrix = part of conduction +part of convection
downloadcode:globalstiffnessgalerkinlaplaceoneelement
output view :viewglobalstiffnessgalerkinlaplaceoneelement
29e) (2D laplace equation ,galerkin method:global stiffness matrix ) two element.
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Fundamentals of the finite element method .
Hartley Gradin,Jr.(pg 333--350)
Macmillan Publishing Company
principal of two-element assemblage equation ('see inside the code')
(three-node triangle) 'triangle a' + (three-node triangle)'triangle b'
example:triangle(a)
pt 1:(0,0),pt 2:(5,0),pt 3(2.5,2)
side 1 to 2 : local = insulation
side 2 to 3 : local = convection Ta,h
side 3 to 1 : local = inside the structure
result : global stiffness (a,b) = part of conduction (a,b)+ part of convection(a,b)
a = triangle ( local 1,2,3 => global 1,2,3)
b = triangle ( local 1,2,3 => global 1,3,4) 'see 28e)'
(assemblage equation ,number of points in the structure ,matrix = (4,4))
downloadcode : globalstiffnessgalerkinlaplacetwoelementfinal
output view :viewglobalstiffnessgalerkinlaplacetwoelementfinal
30e) Generation of Mesh Data " A square mesh of right-angled triangles (finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also :Finite Element Methods for Engineers
Roger T Fenner ( Chapter 4 :Finite Element Meshes,Fortran*)
Imperial College Press,1996
* A square mesh of right-angled triangles.
* A square mesh of mainly isosceles triangular elements.
* A triangular mesh of equilateral elements.
* A circular mesh.
* Mesh Modification.
downloadcode : meshgridfinal
output view : viewmeshgridfinal
31e) 1D "polynoom " Simpson rule 1/3 integration (finite element)
info:staat tussen de code zoals weergegeven op de startpagina.
downloadcode :TESTINTEGRATIONSIMPSONFINAL
output view :viewTESTINTEGRATIONSIMPSONFINAL
32e) 2D Trapezoidal rule for integration (finite element)
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Numerical methods for computer science,engineering,and mathematics.
John H.Mathews
Prentice Hall
info: 2D Trapezoidal rule (ref : algo.)
downloadcode:TRAPZOIDAL2DRULEINTEGRATIONfinal
output view :viewTRAPZOIDAL2DRULEINTEGRATIONfinal control:check
33e) 2D galerkin method:temperature distribution in one element (finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Finite Element Analysis Theory and practice .
M J FAGAN
Longman
see also : 28e),29e)
= One element =
0d : i => point : same value
1D : i,j =>line : stiffness matrix(2,2):load vector(1,2),(Ni,Nj)
2D : i,j,k => triangle : stiffness matrix(3,3):load vector(1,3),(Ni,Nj,Nk)
3D : i,j,k,l => tetrahedron : stiffness matrix(4,4):load vector(1,4),(Ni,Nj,Nk,Nl)
= More than one element = structure =
structure = with 'n points'
global stiffness matrix (n,n)
2D : The connection of a element ('matrix (3,3)') in to the (' matrix (n,n)') in function of the connection
of the elements in the structure. (use three-node triangles). = assembly
example : structure of four points = ' square ' n=4 divided in ' two triangles '
direction = counter clockwise = ccw (for each triangle)
global stiffness matrix :
* * * *
* * * *
* * * *
* * * *
stiffness matrix ' triangle 1'
first triangle : global nodes 1,2,3 => local nodes 1,2,3 (*)
* * * *
* * * *
* * * *
* * * *
stiffness matrix ' triangle 2'
second triangle : global nodes 1,3,4 => local nodes 1,2,3 (*)
* * * *
* * * *
* * * *
* * * *
global stiffness matrix = stiffness matrix ' triangle 1' + stiffness matrix ' triangle 2'
example : computation temperature on a triangle (i,j,k)
side 1 : ij ,side 2 : jk ,side 3 : ki
convection : edge ij
convection : surface (ijk)
Ta = Temperature ambient,h,k given
downloadcode : 2Dfiniteelementoneelement
output view : view2Dfiniteelementoneelement
34e) recast = smaller matrix in a bigger matrix (finite element)
info:staat tussen de code zoals weergegeven op de startpagina.
example : We put mat(3,3) in mat(4,4) . (see: output view)
downloadcode : recast43finalsubfinal
output view : viewrecast43finalsubfinal
35e) recast = smaller matrix (m,m) in a bigger matrix (n,n) (finite element)
option nolet
! n > m
dim a(8,8) ! structure
n=UBOUND(a,1) ! max value for index = i=j
mat print a
! put in a matrix b (4*4)
dim b(4,4)
m = Ubound(b,1)
for i=1 to m
for j=1 to m
b(i,j)= i*j
next j
next i
dim dx(1),dy(1)
mat redim dx(m),dy(m)
! dx(i) ,i=1,2,..m : x position
! dy(i) ,i=1,2,..m : y position
dx(1)=3
dx(2)=6
dx(3)=7
dx(4)=8
dy(1)=4
dy(2)=6
dy(3)=7
dy(4)=8
for i=1 to m
for j=1 to m
a(dx(i),dy(j))=b(i,j)
next j
next i
mat print a
end
output view : recastmxinxn
36e) recast = Example 5.2.1 * (finite element)
info:staat tussen de code zoals weergegeven op de startpagina.
see also : * Fundamentals of the finite element method for heat and fluid flow
Roland W.Lewis,Perumal Nithiarasu,Kankanhalli N.Seetharamu
Wiley .
( nuclear reactors,chemical and biological reactors)
example : *How to find equation (5.16) from pg 131 from example 5.2.1.
download code :recastexample521
output view :viewrecastexample521
37e) Gaussian Elimination and pivoting (finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Numerical methods for computer science,engineering,and mathematics.
John H.Mathews
Prentice Hall .
example : Find the parabola y=A+B*x+C*x^2 that passes through the three
points (1,1),(2,-1),(3,1).
download code : gausseliminationparabolafinal
output view : viewgausseliminationparabolafinal
38e) Determinant (use LU factorization) (finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Numerical methods for computer science,engineering,and mathematics.
John H.Mathews
Prentice Hall .
example : Use above example for explination of the principal of the determinant.
concept : Determinant can be used for Cramer's rule (solve linear systems).
download code :DETERMINANTparabolafinal.
output view :viewDETERMINANTparabolafinal.
39e) LU for (linear systems) (finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Numerical methods for computer science,engineering,and mathematics.
John H.Mathews
Prentice Hall .
example : Solve 37e) with LU factorization ,the same result.
AA*X=B=LU*X=B
a) L*Y= B => Y
b) U*X=Y => X
download code :LUfactorizationparabolafinal
output view :viewLUfactorizationparabolafinal
40e) LU Find L,U ' Crout's method' (linear systems) (finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also:Computational methods for heat and mass transfer.
Pradip Majumdar ( pg 83 to 87)
Taylor & Francis
algo : 'Crout's method ' see inside the code ' for LU Decomposition .
download code :LUdecompositionCROUTfinal
output view :viewLUdecompositionCROUTfinal
41e) LU Find L,U ' Crout's method ' solve linear system A*x=b,(finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also:Computational methods for heat and mass transfer.
Pradip Majumdar ( pg 83 to 87)
Taylor & Francis
algo 1e : 'Crout's method ' see inside the code ' for LU Decomposition.
algo 2e : step lin 1,step lin 2 'see inside the code ' solve the system inside the code.
download code :LUdecompositionCROUTfinaluselinsystemfinal
output view :viewLUdecompositionCROUTfinaluselinsystemfinal
42e) LU Find 'inverse matrix A^(-1) , A*A^(-1) = I",(finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also:Analytical ,Numerical and computational methods for science and engineering
Gene H.Hostetter ................... ( pg 281 to 283 )
Prentice Hall
algo 1e : 'Crout's method ' see inside the code ' for LU Decomposition.
algo 2e: ' see inside the code ' A=LU => A^(-1) = (L*U)^(-1) = U^(-1)*L^(-1)
download code : LUmatrixinversefinal
output view : viewLUmatrixinversefinal
43e) LL^T Find 'A= L*L^T",(finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also:Digital computation for chemical engineers.
Leon Lapidus ................... ( pg 249 )
McGraw-Hill Series in Chemical Engineering
algo 1e : 'Cholesky Decomposition ' A=L*L^T' .
use: x^T*A*x>0 and a(i,j)=a(j,i) , a(i,j) >0
download code:Cholesky decompositionfinal
output view :viewCholesky decompositionfinal
44e) LDL^T Find 'A= L*D*L^T",(finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also: Numerical algorithms with fortran ( c )
Gisela Engeln-Mullges , 'Frank Uhlig'
Springer ( pg 81)
algo 1e : 'Cholesky Decomposition 'A=L*D*L^T' .
download code:Cholesky decompositiondiagonalfinal
output view :viewCholesky decompositiondiagonalfinal
45e) LL^T Find 'A= L*L^T : banded cholesky for Ax=b",(finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also: algo: find on the internet from a university website .
algo 1e) ' LB = Lowerband Storage = see inside the code 'use to solve Ax=b .
download code:CHOLESKY DECOMPOSITIONbandedfinal
output view :viewCHOLESKY DECOMPOSITIONbandedfinal
46e) LR Find 'five diagonal system A= L*R solve for Ax=b",(finite element).
info:staat tussen de code zoals weergegeven op de startpagina.
see also: Numerical algorithms with fortran ( c )
Gisela Engeln-Mullges , 'Frank Uhlig'
Springer ( pg 98-101)
algo 1e) ' five-diagonal systems = see inside the code 'use to solve Ax=b .
use: strongly nonsingular matrix A
download code :fivediagonalsystemfinal
output view :viewfivediagonalsystemfinal
47e) Find 'secant algorithm solution of two non-linear equations".
info:staat tussen de code zoals weergegeven op de startpagina.
example : eq1: f1(x,y) = x*x+y*y-50 = 0
eq2: f2(x,y) = x*y-25 = 0
find : xi,yi
result : f1(xi,yi) = 0
f2(xi,yi) = 0
algo : ' secant algorithm '
use: Jacobian approximation = ' Finite Differences '
download code : nonlinearusejacobifinitedifferencefinal
output view : viewnonlinearusejacobifinitedifferencefinal
48e) Find 'jacobian (symbolic *) algorithm solution of two non-linear equations".
info:staat tussen de code zoals weergegeven op de startpagina.
example : eq1: f1=f1(x,y) = x*x+y*y-50 = 0
eq2: f2=f2(x,y) = x*y-25 = 0
need symbolic : differentation to build the jacobian matrix.
*************** software symbolic computation ****************
use : symbolic * algebra = " WMAXIMA = use 'see inside the code' "
( Calculation : Partial Derivatives =>)
calculation : ((df1/dx,df1/dy),(df2/dx,df2/dy)) = J = 'Jacobian'**
**********************************************************
find : xi,yi
result : f1(xi,yi) = 0
f2(xi,yi) = 0
** every step of the computation : ' new J and J^(-1) '
algo: ' Jacobian algorithm ' = ' Newton-Raphson algorithm in multidimensional space'
Newton method : x(i+1) = x(i) -f(x(i))/df(x(i)) = x(i) - (df(x(i))^(-1) )*f(x(i))
Jacobian method : X(i+1) = X(i) -J^(-1) *F(X(i))
x(i) = number => X(i) = vector (2,1)
f(x(i)) = number => F(X(i)) = vector(2,1)
df(x(i)) =number => J = matrix(2,2) ' (df(x(i))^(-1)) = number'
info : J*J^(-1) = unity (2,2) -> number = 1 => matrix = ' a(i,i) = 1.0 and i<>j zero'
unity (2,2) = ((1,0),(0,1)) = eye(2)
download code :nonlinearusejacobisymbolicfinal
output view :viewnonlinearusejacobisymbolicfinal
49e) Find 'broyden (symbolic *) algorithm for the solution of three non-linear equations".
info:staat tussen de code zoals weergegeven op de startpagina.
example : eq1: f1=f1(x,y,z) = 3*x-cos(y*z)-1/2=0
eq2: f2=f2(x,y,z) = x^2-81*(y+0.1)^2+sin(z)+1.06=0
eq3: f3=f3(x,y,z) = exp(-x*y)+20*z+(10*pi-3)/3 = 0
need symbolic : differentation to build the jacobian matrix.
*************** software symbolic computation ****************
use : symbolic * algebra = " WMAXIMA = use 'see inside the code' "
( Calculation : Partial Derivatives =>)
calculation : ((df1/dx,df1/dy,df1/dz),(df2/dx,df2/dy,df2/dz),(df3/dx,df3/dy,df3/dz)) = J = 'Jacobian'**
**********************************************************
find : xi,yi,zi (example : 0.1,0.1,-0.1) ' start value '
result : f1(xi,yi,zi) = 0
f2(xi,yi,zi) = 0
f3(xi,yi,zi) = 0
** one time computation : ' from J with values above we get J^(-1) '
use : Inversion formula from 'Sherman and Morrison'
use : general formula '3 by 3 matrix inversion' = 'use wMaxima'
algo : simple Broyden ' see inside the code '
cvg : check with 2-norm
and
Numerical Analysis (sec edition)
Richard L.Burden,J.Douglas Faires
Prindle ,Weber&Schmidt (' Boston ,Massachusetts)
download code:simplebroydenfinal
output view :viewsimplebroydenfinal
50e) Binary distillation " find number of stages or trays = Smoker equation"
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Design equilibrium stage process : pg 155-161
Buford D.Smith
McGraw -Hill
example : Find number of theoretical stages or trays through the use of the Smoker
equation.
algo : Smoker equation
input : R,q,xd,zf,xb,alpha
output : Ns = Number of stages from stripping section
Ne = Number of stages from enriching section
Nt = Ne + Ns
info : ' "click" on the word Smoker'
download code :smokerdistbinfinal
output view :viewsmokerdistbinfinal
51e) Binary distillation "Theory Txy diagram (pressure = cte)".
info:staat tussen de code zoals weergegeven op de startpagina.
download code :TXYDIAGRAMtheoryfinal
output view :viewTXYDIAGRAMtheoryfinal
52e) Binary distillation "database to draw Txy diagram(pressure = cte)".
info:staat tussen de code zoals weergegeven op de startpagina.
download code :databaseTXYDIAGRAMTHEORYfinal
output view : viewdatabaseTXYDIAGRAMTHEORYfinal
53e) Binary distillation "database to draw Pxy diagram (temperature =cte)".
info:staat tussen de code zoals weergegeven op de startpagina.
download code :databasePXYDIAGRAMTHEORYfinal
output view : viewdatabasePXYDIAGRAMTHEORYfinal
54e) find "bubblepoint of Azeotrope,pressure given (Benzene-Ethanol)",
info:staat tussen de code zoals weergegeven op de startpagina.
see also: Introductory chemical engineering thermodynamics.(sec edition)
J.Richard Elliott.Carl T .Lira
Prentice Hall.
example : Find the bubblepoint "Benzene,Ethanol" by given pressure
Thermodynamic model : Margules 2 for binary systems and
for the vapor pressure ' Antoine eq.'
download code : azeotropemargules2final
output view : viewazeotropemargules2final
55e) find "isothermal flash (ideal),pressure given (feed:Benzene,toluene,o-xylene)",
info:staat tussen de code zoals weergegeven op de startpagina.
see also: Introductory chemical engineering thermodynamics.(sec edition)
J.Richard Elliott.Carl T .Lira
Prentice Hall.
example : Isothermal flash calculation of three products ( vaporization of the feed).
pressure : dew point < p < bubblepoint : Two phase region. *
method:* for cvg use interhalving method to find the sum x(i).
after this y(i) = x(i)*ps(i)/p => y(1),..y(nc)
download code : ISOTHERMALFLASHPRESSUREfinal
output view : viewISOTHERMALFLASHPRESSUREfinal
56e) find "isothermal flash"newton raphson" (ideal),pressure given ,
info:staat tussen de code zoals weergegeven op de startpagina.
see also: Introductory chemical engineering thermodynamics.(sec edition)
J.Richard Elliott.Carl T .Lira
Prentice Hall.
example : Isothermal flash calculation of three products ( vaporization of the feed).
pressure : dew point < p < bubblepoint : Two phase region. *
method:* for cvg use newton raphson method to find the sum x(i).
after this y(i) = x(i)*k(i) => y(1),..y(nc)
download code : isothermalflashnewtonfinal
output view : viewisothermalflashnewtonfinal
57e) Binary distillation "database (use linspace) to draw Txy diagram(pressure = cte)".
info:staat tussen de code zoals weergegeven op de startpagina.
! How to set dimension from the matrix 'see inside the code'
nv = number of values
tb(i) , i = 1..2
temp(i) ,i = 1..nv
sub linspace(nv,tb(),temp())
ymax = max(tb(1),tb(2))
ymin = min(tb(1),tb(2))
for i=1 to nv
m = (i-1)/(nv-1)
temp(i)= ymin + m *(ymax-ymin)
next i
end sub
see also : 52e)
download code :drawtxytableslinspacefinal
output view :viewdrawtxytableslinspacefinal
58e) Binary distillation rating for the (system)"interhalving method" ,pressure given ,
info:staat tussen de code zoals weergegeven op de startpagina.
see also: Introductory chemical engineering thermodynamics.(sec edition)
J.Richard Elliott.Carl T .Lira
Prentice Hall.
example : We want to find r = reflux ratio for the system of ' methanol /aq.'
Given : xb ( = .001) ,xd ( = .9990)
Given column : a) number of plates ( = ' bottom ' + nt '=0 + 40 =40)
b) number of feed plate ( = ' nf ' =15)
system : The van Laar equation for binair system ( methanol/aq)
method : rating => interhalving method ( => dr=dr/2)
* r , r-dr,r+dr 'see code below'
equilibrium => newton raphson method ( => t=t-f/df )
solution problem : Find the reflux ratio to achieve product purities of
xd,xb for given column starting from initial guess of
r .If the correct value has been guessed then vapor leaving
tray nt is equal to xd,otherwise a other r value guessed
based on difference between two compostions ( =.0001).
download code* : bindistratingfinaloutputfinalfinal
output view (nice output) : bindistratingniceoutputversion1
output view : bindistratingoutputversion1
59e) Binary distillation design ( find : number of plates,number of the feed plate)
info:staat tussen de code zoals weergegeven op de startpagina.
see also: Introductory chemical engineering thermodynamics.(sec edition)
J.Richard Elliott.Carl T .Lira
Prentice Hall.
example : We want to find number of plates and number of the feed plate.
Given : xb ( = .02) ,xd ( = .98),alpha( = 2.00),z ( = .5)
Given column : q ( = 1.0), rr = reflux ratio ( = 2.560)
system : equil => y=alpha*x/(1+(alpha-1)*x)
solution problem : find intersection of 'sol line and q line'
sol line = stripping section operating line
q line = quality of the feed operating line
We get the part below the feed ( = stripping section).
first part : we get (nf =feed plate)
We get the part above the feed (= rectifying section).
second part: we get (nt = number of plates )
download code: bindistdesignfinal
output view : bindistdesignfinalout
60e) Multicomponent distillation design Thiele/Geddes ( find :estimate w, new k values )
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Distillation Dynamics and Control
Pradeep B.Deshpande
Instrument Society oF America
example : We estimate molal flows of the bottoms (=w) of the column and find so new
values for k ( = equilibrium equation).We find this g(i) and then sum(g(i),i=1..nc) <> ? 1.
( check with do loop)
source code : ( info see code below)
nc = number of components.
pt = total pressure (atm)
t = t(°C)+273 ( Kelvin)
i = 1...nc
k (i) = (exp(a(i) +b(i)/t))/pt
? what : g(i)
set : g(i) = xw , e= w/f, f = w+d
a) f*xf = d*xd + w*xw
b) f*xf/xw = d*xd/xw +w*xw/xw
set : db = xd/xw
c) f*xf/xw = d*db +w
set: d=f-w
d) f*xf/xw = ( f-w)*db+w
division by f ,e =w/f
e) xf/xw =(1-e)*db +e
f) xf/xw =( (1-db)*e+db)
m = ((1-db)*e+db)
g) xw = xf/m = g(i)
download code: simplethielegeddesfinal
output view (1) : simplethielegeddesfinalniceoutput
output view (2) : simplethielegeddesfinalnormaloutput
61e) Diagonal Dominant Matrix = A, 'Linear system :Ax=b' (Sub )
info:staat tussen de code zoals weergegeven op de startpagina.
see also : Computational Mathematics
B.P.Demidovich
Translated from Russian By Geoorge Yankovsky
MIR PUBLISHERS .MOSCOW
example: A = a(i,j) , i=1..n,j=1..n ( n = 3) 'Read - Data Statement'
AlGO : abs(a(i,i)) > sum(abs(a(i,j)),j=1..n,j<>i)
Check : Assume 'n' Linear equations in 'n' unknowns
can be solved with 'Iterative Method'
see also : COMPUTER PROGRAMMA'S
download code: diagonaldominantversion1final
output view : diagonaldominantoutput
62e) Dynamic Model Cross-Flow Heat Exchanger (Finite Differences)
info:staat tussen de code zoals weergegeven op de startpagina.
example : Heat Exchanger (Cross-Flow :type)
( liquid flow in the pipe cool down along the length with 20°C'Air')
Why Cross-Flow : liquid stream perpendicular to cooling
PDE : Heat Balance
Heat Balance : (1)=(2)+(3)
(1) accumulation(heat) 'volume element'
(2) Heat loss :due enegy transfer to air
(3) Heat gain: entering fluid (100°C)
(1) : pi*D^2/4*dx*rho*cp*dt
(2) : F*cp*(T-(T+DT))*dt
(3) : U*pi*dx*dt*(Ta-T(x,t))
PDE : computation (discrete version)
( first order approximation)
FTBS = Forward in time and Back wards in Space
( 'see code to understand this')
i = space
j = time
dtemp_dx=(temp(j,i)-temp(j,i-1))/dx !
temp(j+1,i)=temp(j,i)+(dt/pipe)*(fcp*dtemp_dx+upid*(ta-temp(i,j)))
download code: DynamicmodelCrossflowheatexfinalfinalfinal.tru
output view : DynamicmodelCrossflowheatexout
63e) IVP RK4 'first Order' RK4 (Finite Differences)
info:staat tussen de code zoals weergegeven op de startpagina.
example : IVP =Initial Value Problem
ODE : dy/dt = f(t,y) =(t-y)/2
ic : y(t=0)
ic = initial condition
Method : RK4 method Integration
download code: IVPRK4FINALFINALFINAL.tru
output view : ivprk4finalfinalfinal
64e) BVP Shooting' Second Order Equation ' Euler(Finite Differences)
info:staat tussen de code zoals weergegeven op de startpagina.
info :youtube :numericalmethodsguy
Shooting Method: Example: Part 1 of 4
https://www.youtube.com/watch?v=ZMgikZ-lcS8
Shooting Method: Example: Part 2 of 4
https://www.youtube.com/watch?v=R7I8FrlB_KM
Shooting Method: Example: Part 3 of 4
https://www.youtube.com/watch?v=RNBe0qOOIls
Shooting Method: Example: Part 4 of 4
https://www.youtube.com/watch?v=Sc0Z_qIndz0
print " What is the shooting method => "
example : BVP => convert ' 2 IVP =Initial Value Problem '
BVP : Boundary Value Problem
ODE : d(2)y/dx(2) =2*y+8*x*(9-x)
bc: y(0)=0 and y(9)=0
bc = boundary conditions
Method : Euler method integration ( system 2 first order IVP)
youtube : version 1.0
download code: SHOOTINGMETHODfinalnotupdatefinalversion1.tru
optimal : version 1.0a 'code'
download code: SHOOTINGMETHODfinalnotupdatefinalversion1a.tru
output view : SHOOTINGMETHODupdatefinalversion1
65e) Dynamic Model Co-Current Heat Exchanger,Euler (Finite Differences)
info:staat tussen de code zoals weergegeven op de startpagina.
DYNAMIC MODEL : PDE'S 'see code'
( Co-current double pipe heat exchanger )
Co –current : (entry segment same side,numbering 1 =< i =< N)
S =(S)hell T = (T)ube , Ts = temperature shell ,Tt=temperature tube
Tsi = temperature shell input
Tti = temperature tube input
Entry segment :
dTs(1)/dt = bs*(Tsi - Ts(1)) –cs*(Ts(1) –Tt(1))
dTt(1)/dt = bt*(Tti - Ts(1)) +ct*(Ts(1) –Tt(1))
I th segment : i = 2 to N-2
dTs(i)/dt = bs*(Ts(i -1) - Ts(i)) –cs*(Ts(i) –Tt(i))
dTt(i)/dt = bt*(Tt(i -1) - Ts(i)) +ct*(Ts(i) –Tt(i))
exit segment : ( above i = N )
dTs(N)/dt = bs*(Ts(N-1) - Ts(N)) –cs*(Ts(N) –Tt(N))
dTt(N)/dt = bt*(Tt(N-1) - Ts(N)) +ct*(Ts(N) –Tt(N))
Method : Euler Integration : 'see code '
download code: COCURRENTHEATEXCHANGERFINALFINALFINAL.tru
output view : cocurrentheatexchangerversion1
66e) Dynamic Model Counter-Current Heat Exchanger,Euler (Finite Differences)
info:staat tussen de code zoals weergegeven op de startpagina.
DYNAMIC MODEL : PDE'S 'see code'
(Counter-current double pipe heat exchanger)
Counter–Current :(entry segment different side,numbering 1 =< i =< N)
S =(S)hell T = (T)ube , Ts = temperature shell ,Tt=temperature tube
Tsi ,Tso = temperature shell input,output
Tti ,Tto = temperature tube input,output
Entry segment:
dTs(N)/dt = bs*(Tsi - Ts(N)) –cs*(Ts(N) –Tt(N))
dTt(1)/dt = bt*(Tti - Ts(1)) +ct*(Ts(1) –Tt(1))
I th segment :
dTs(i)/dt = bs*(Ts(i +1) - Ts(i)) –cs*(Ts(i) –Tt(i))
dTt(i)/dt = bt*(Tt(i -1) - Ts(i)) +ct*(Ts(i) –Tt(i))
Exit segment :
dTs(1)/dt = bs*(Ts(2) - Ts(1)) –cs*(Ts(1) –Tt(1))
dTt(N)/dt = bt*(Tt(N-1) - Ts(N)) +ct*(Ts(N) –Tt(N))
Method : Euler Integration : 'see code '
download code:COUNTERCURRENTHEATEXCHANGERFINALFINALFINAL
output view : countercurrentheatexchangerversion1 matrixvectorwrite.png
67e) Complex Domain Geometry (data matrix ) 'finite element method'
info:staat tussen de code zoals weergegeven op de startpagina.
Use CFD = Computational Fluid Dynamics (complex boundary)
Application : ref : The finite element method in engineering (5ed)
Singiresu S.Rao,Elsevier (Pg 580)
download code: FORMATDATAFLOW2DFINAL3VERSION1
datafile :femflow2d.dat
output view1 :datamatrixchoose1
output view2:datamatrixchoose2
write : (matrix and vector) to file :
download code: MATRIXVECTORWRITEDATAFINALFINAL
datafile (matrix a) :MAT_A
datafile (vector b) :VEC_B
output view:matrixvectorwrite
68e) Potential flow around a circular cylinder 'finite elemeFnt method,CD'
info:staat tussen de code zoals weergegeven op de startpagina.
Singiresu S.Rao,Elsevier (Pg 571-584)
download code:FLOW2DOBJECTFINALVERSION1FINAL
output view: potentiaflowover2dobject
69e) the viscosity (mu) and distance between two parallel plates(laminar flow,heated)
info:staat tussen de code zoals weergegeven op de startpagina.
Application : ref : Basics of the finite element method
Solid Mechanics,Heat Transfer and Fluid Mechanics
Paul E.Allaire (pg 198)
download code:VISCOSITYTEMPERATUREDISTANCEFROMWALLFINAL
output view: viscositydistancefinalwithoutgraphics