Scilab Textbook Companion for Introduction To Numerical Methods In Chemical Engineering by P. Ahuja1 Created by Chandra Prakash Sipani B.TECH Part-2 Chemical Engineering IIT-BHU College Teacher Ahuja, Pradeep Cross-Checked by Pradeep Ahuja August 10, 2013
1
Funded by a grant from the National Mission on Education through ICT, http://spoken-tutorial.org/NMEICT-Intro. This Textbook Companion and Scilab codes written in it can be downloaded from the ”Textbook Companion Project” section at the website http://scilab.in
Book Description Title: Introduction To Numerical Methods In Chemical Engineering Author: P. Ahuja Publisher: PHI Learning, New Delhi Edition: 1 Year: 2010 ISBN: 8120340183
1
Scilab numbering policy used in this document and the relation to the above book. Exa Example (Solved example) Eqn Equation (Particular equation of the above book) AP Appendix to Example(Scilab Code that is an Appednix to a particular
Example of the above book) For example, Exa 3.51 means solved example 3.51 of this book. Sec 2.3 means a scilab code whose theory is explained in Section 2.3 of the book.
2
Contents Listt of Scilab Codes Lis
4
1 lin linear ear alge algebrai braic c equati equations ons
7
2 NONLI NONLINEAR NEAR ALGEBR ALGEBRAIC AIC EQUA EQUATIONS TIONS
14
3 CHEMI CHEMICAL CAL ENGIN ENGINEERIN EERING G THERMOD THERMODYNAMI YNAMICS CS
19
4 INI INITIA TIAL L VALU VALUE E PROBL PROBLEMS EMS
37
5 BOUND BOUNDAR ARY Y VALUE VALUE PROBL PROBLEMS EMS
58
6 CON CONVEC VECTIO TION N DIFFUSI DIFFUSION ON PROBLE PROBLEMS MS
68
7 TUB TUBULA ULAR R REACTO REACTOR R WITH AXIAL AXIAL DISPERS DISPERSION ION
74
8 CHEMIC CHEMICAL AL REACTI REACTION ON AND DIFFUS DIFFUSION ION IN SPHERI SPHERI-CAL CATALYST PELLET
78
9 ONE DIMENSI DIMENSIONAL ONAL TRANSIE TRANSIENT NT HEA HEAT CONDUCTION CONDUCTION
83
10 TWO DIMENSIONAL STEADY AND TRANSIENT HEAT CONDUCTION 90
3
List of Scilab Codes Exa 1. Exa 1.11 Exaa 1. Ex 1.22 Exaa 1. Ex 1.33 Exaa 1. Ex 1.44 Exaa 1. Ex 1.55 Exaa 1. Ex 1.66 Exaa 2. Ex 2.11 Exaa 2. Ex 2.22 Exaa 2. Ex 2.33 Exaa 2. Ex 2.44 Exaa 2. Ex 2.55 Exaa 2. Ex 2.66 Exaa 3. Ex 3.11 Exaa 3. Ex 3.22 Exa 3.3 Exa 3.4 Exa 3.5 Exa 3.6 3.6 Exa 3.7 Exa 3.9 3.9 Exa 3.1 3.100 Exa 4.1 Exa 4.2 Exaa 4. Ex 4.33 Exaa 4. Ex 4.44 Exaa 4. Ex 4.55 Exaa 4. Ex 4.66 Exa 4.7
TDM DMA A me meth thod od . . . . . . . . . . . . . . . . . . . . . . gaus ga usss el elim imin inat atio ion n me meth thod od . . . . . . . . . . . . . . . . gaus ga usss el elim imin inat atio ion n me meth thod od . . . . . . . . . . . . . . . . gaus ga usss el elim imin inat atio ion n me meth thod od . . . . . . . . . . . . . . . . gaus ga usss se seid idel el me meth thod od . . . . . . . . . . . . . . . . . . . gaus ga usss se seid idel el me meth thod od . . . . . . . . . . . . . . . . . . . alge al gebr brai aicc eq equa uati tion onss . . . . . . . . . . . . . . . . . . . . alge al gebr brai aicc eq equa uati tion onss . . . . . . . . . . . . . . . . . . . . alge al gebr brai aicc eq equa uati tion onss . . . . . . . . . . . . . . . . . . . . alge al gebr brai aicc eq equa uati tion onss . . . . . . . . . . . . . . . . . . . . alge al gebr brai aicc eq equa uati tion onss . . . . . . . . . . . . . . . . . . . . alge al gebr brai aicc eq equa uati tion onss . . . . . . . . . . . . . . . . . . . . ther th ermo mody dyna nami mics cs . . . . . . . . . . . . . . . . . . . . . ther th ermo mody dyna nami mics cs . . . . . . . . . . . . . . . . . . . . . flash flas h cal calcul culati ations ons usi using ng Rao Raoult ult la law w . . . . . . . . . . . BPT and DPT cal calcul culati ation on usi using ng modifi modified ed rao raoult ult la law w . flash flas h calc calcula ulatio tions ns usi using ng modi modified fied Rao Raoult ult la law w . . . . . . vapou apourr press pressure ure cal calcul culati ation on usin usingg cubic cubic equ equati ation on of sta state te pressu pre ssure re x y dia diagra gram m usin usingg gamma gamma phi app approa roach ch . . . chemi ch emical cal reac reactio tion n engine engineeri ering ng 2 simu simulta ltaneo neous us react reaction ionss adiaba adi abatic tic flam flamee tempe temperat rature ure . . . . . . . . . . . . . . . soluti sol ution on of ord ordina inary ry diff differe erent ntial ial eq equat uation ion . . . . . . . . soluti sol ution on of ord ordina inary ry diff differe erent ntial ial eq equat uation ion . . . . . . . . doub do uble le pi pipe pe he heat at ex exccha hang nger er . . . . . . . . . . . . . . . stir st irre red d ta tank nk wi with th co coil il he heat ater er . . . . . . . . . . . . . . . stir st irre red d ta tank nk wi with th co coil il he heat ater er . . . . . . . . . . . . . . . pneu pn euma mati ticc co conv nvey eyin ingg . . . . . . . . . . . . . . . . . . . simult sim ultane aneous ous ord ordina inary ry diff differe erent ntial ial equ equati ations ons . . . . . . 4
7 8 9 10 12 12 14 15 15 16 17 18 19 21 23 24 27 29 31 33 35 37 38 39 40 41 42 43
Exa 4.8 Exa 4.9 Exa 4.10 4.10 Exa 4.1 4.111 Exa 4.1 4.122 Exa 4.1 4.133 Exa 4.1 4.144 Exaa 4.1 Ex 4.155 Exaa 4.1 Ex 4.166 Exaa 4.1 Ex 4.177 Exa 4.1 4.188 Exaa 5. Ex 5.11 Exaa 5. Ex 5.22 Exaa 5. Ex 5.33 Exaa 5. Ex 5.44 Exaa 5. Ex 5.55 Exaa 5.6 Ex 5.6 Exaa 5.7 Ex 5.7 Exa 5.8 Exaa 6. Ex 6.11 Exaa 6. Ex 6.22 Exa 7.1 7.1 Exa 7.2 Exa 8. Exa 8.11 Exaa 8. Ex 8.22 Exaa 8. Ex 8.33 Exaa 8. Ex 8.44 Exa 9.1 Exa 9.2 Exa 9.3 Exa 9.4 Exaa 9. Ex 9.55 Exa 10.1 10.1 Exa 10.2 10.2 Exa 10. 10.33 Exa 10. 10.88 Exa 10. 10.99
simultane simult aneous ous ord ordina inary ry diff differe erent ntial ial equ equati ations ons . . . . . . simult sim ultane aneous ous ord ordina inary ry diff differe erent ntial ial equ equati ations ons . . . . . . series ser ies of of stirre stirred d tank tank with with coil coil heate heaterr . . . . . . . . . . batch bat ch and stir stirred red tan tank k react reactor or . . . . . . . . . . . . . batch bat ch and stir stirred red tan tank k react reactor or . . . . . . . . . . . . . batch bat ch and stir stirred red tan tank k react reactor or . . . . . . . . . . . . . batch bat ch and stir stirred red tan tank k react reactor or . . . . . . . . . . . . . plug pl ug flo flow w rea react ctor or . . . . . . . . . . . . . . . . . . . . . plug pl ug flo flow w rea react ctor or . . . . . . . . . . . . . . . . . . . . . plug pl ug flo flow w rea react ctor or . . . . . . . . . . . . . . . . . . . . . non isot isother hermal mal plug plug flow flow rea react ctor or . . . . . . . . . . . . disc di scre reti tiza zati tion on in 1 D sp spac acee . . . . . . . . . . . . . . . . disc di scre reti tiza zati tion on in 1 D sp spac acee . . . . . . . . . . . . . . . . disc di scre reti tiza zati tion on in 1 D sp spac acee . . . . . . . . . . . . . . . . disc di scre reti tiza zati tion on in 1 D sp spac acee . . . . . . . . . . . . . . . . disc di scre reti tiza zati tion on in 1 D sp spac acee . . . . . . . . . . . . . . . . 1 D st stea eady dy st stat atee hea heatt con condu duct ctio ion n . . . . . . . . . . . . 1 D st stea eady dy st stat atee hea heatt con condu duct ctio ion n . . . . . . . . . . . . chemi ch emical cal rea reacti ction on and diff diffusi usion on in pore . . . . . . . . . upwi up wind nd sc sche heme me . . . . . . . . . . . . . . . . . . . . . . upwi up wind nd sc sche heme me . . . . . . . . . . . . . . . . . . . . . . boundary value proble problem m in che chemical mical react reaction ion engin engineeri eering ng boundar boun dary y valu valuee proble problem m in chem chemica icall reacti reaction on engin enginee eerring second order . . . . . . . . . . . . . . . . . . . . . first fir st or orde derr re reac acti tion on . . . . . . . . . . . . . . . . . . . . seco se cond nd or orde derr re reac acti tion on . . . . . . . . . . . . . . . . . . non no n is isot othe herm rmal al co cond ndit itio ion n . . . . . . . . . . . . . . . . non no n is isot othe herm rmal al co cond ndit itio ion n . . . . . . . . . . . . . . . . transi tra nsien entt con conduc ductio tion n in rec rectan tangul gular ar sla slab b . . . . . . . . transi tra nsien entt con conduc ductio tion n in rec rectan tangul gular ar sla slab b . . . . . . . . transi tra nsien entt con conduc ductio tion n in cyl cylind inder er . . . . . . . . . . . . . transi tra nsien entt con conduc ductio tion n in sph sphere ere . . . . . . . . . . . . . tran tr ansi sien entt di diffu ffusi sion on in sp sphe here re . . . . . . . . . . . . . . . discre dis cretiz tizati ation on in 2 D space space gauss gauss seid seidel el method method . . . . discre dis cretiz tizati ation on in 2 D space space gauss gauss seid seidel el method method . . . . discre dis cretiz tizati ation on in in 2 D space space . . . . . . . . . . . . . . . . discre dis cretiz tizati ation on in in 2 D space space . . . . . . . . . . . . . . . . discre dis cretiz tizati ation on in in 2 D space space . . . . . . . . . . . . . . . . 5
44 45 46 47 48 49 50 52 52 54 55 58 59 59 60 62 63 64 66 68 71 74 76 78 79 80 81 83 84 86 87 88 90 91 92 93 94
Exa 10. 10.10 10 ADI meth method od . . . . . . . . . . . . . . . . . . . . . . . Exa 10.11 ADI method for transien transientt heat conduction conduction . . . . . . .
6
95 96
Chapter 1 linear algebraic equations
Scilab code Exa 1.1 TDMA method
1 / / ch ch 1 ex 1.1 − s o l v i n g
s e t o f a l g e b r a i c e qu q u a t io i o n s by T r i d i a g o n a l M at a t ri r i x A l g or o r i th t h m Method .
2 clc 3 disp ( ” t he h e s o l n o f e g 1. 1 .1−− > ” ) ; 4 f o r i = 2: / / s ub ub d i a g o n a l 2: 7 , a ( i ) = 1 ;
assignment 5 end 6 f o r j = 1: 1: 7 , b ( j ) = - 2 ;
/ / ma m a in in d i a g o n a l
assignment 7 end 8 f o r k = 1: 1: 6 , c ( k ) = 1 ;
/ / s u p er er d i a g o n a l
assignment 9 end 10 d(1)=-240
//given values
assignment 11 12 13 14 15 16
f o r l = 2: 2: 6 , d ( l ) = - 4 0; 0; end d(7)=-60 i=1; n=7;
7
// i n i t i a l b is equal
17 beta1(i)=b(i);
t o b et e t a s i n c e a 1 =0 / / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); 26 e n d 27 28 disp ( ” t h e s o l u t i o n o f e x 1 . 1 b y TD TDMA me me t h o d i s ” ) ; 29 f o r i=1:7, disp (x(i)); 30 e n d 18 19 20 21 22 23 24 25
Scilab code Exa 1.2 gauss elimination method
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
/ / ch ch 1 e x 1 . 2
clc h e s o l n o f e g 1. 1 .2−− > ” ) ; disp ( ” t he a 1 = 10 10 , a 2 =1 =1 , a 3 =2 =2 , b 1 = 2 , b 2 =1 =1 0 , b 3 =1 =1 , c 1 = 1 , c 2 =2 =2 , c 3 =1 =1 0 , d 1 = 44 4 4 , d 2 =5 =5 1 , d 3 =6 =6 1 ,
b3=b3-(b1/a1)*a3 b2=b2-(b1/a1)*a2 d2=d2-(b1/a1)*d1 b1=b1-(b1/a1)*a1
/ / 1 s t ro w / / 2 n d r ow ow / / 3 r d ro w // g i v e n v a l u e s / / f o r m ak a k in in g b 1 =0
/ / f o r m ak a k i ng n g c 1= 1 =0
c3=c3-(c1/a1)*a3 c2=c2-(c1/a1)*a2 d3=d3-(c1/a1)*d1
8
17 18 19 20 21 22 23 24 25 26
c1=c1-(c1/a1)*a1
/ / f o r m ak ak in in g c 2= 2 =0
c3=c3-(c2/b2)*b3 d3=d3-(c2/b2)*d2 c2=c2-(c2/b2)*b2
/ / f i n a l v a l ue ue s o f x x3=d3/c3; x2=(d2-(b3*x3))/b2; x1=(d1-(x3*a3)-(x2*a2))/a1; h e s o l u t i o n u s i n g g au a u ss ss e l i m i n a t i o n disp (x3,x2,x1, ” t he m et et ho ho d i s ” ) ;
Scilab code Exa 1.3 gauss elimination method
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
/ / ch ch 1 e x 1 . 3
clc h e s o l n o f e g 1. 1 .3−− > ” ) ; disp ( ” t he / / 1 s t ro w a 1 = 3 , a 2 =1 =1 , a 3 = -2 -2 , / / 2 n d r ow ow b 1 = - 1, 1 , b 2 =4 =4 , b 3 = -3 -3 , / / 3 r d ro w c 1 = 1 , c 2 = -1 -1 , c 3 =4 =4 , //given values d 1 = 9 , d 2 = -8 -8 , d 3 =1 =1 ,
b3=b3-(b1/a1)*a3 b2=b2-(b1/a1)*a2 d2=d2-(b1/a1)*d1 b1=b1-(b1/a1)*a1
c3=c3-(c1/a1)*a3 c2=c2-(c1/a1)*a2 d3=d3-(c1/a1)*d1 c1=c1-(c1/a1)*a1
/ / f o r m ak a k in in g b 1 =0
/ / f o r m ak a k i ng n g c 1= 1 =0
/ / f o r m ak ak in in g c 2= 2 =0
c3=c3-(c2/b2)*b3 d3=d3-(c2/b2)*d2 c2=c2-(c2/b2)*b2
9
22 23 24 25 26
/ / f i n a l v a l ue ue s o f x x3=d3/c3; x2=(d2-(b3*x3))/b2; x1=(d1-(x3*a3)-(x2*a2))/a1; h e s o l u t i o n u s i n g g au a u ss ss e l i m i n a t i o n disp (x3,x2,x1, ” t he m et et ho ho d i s ” ) ;
Scilab code Exa 1.4 gauss elimination method
1 2 3 4
/ / ch 1 ex 1 .4 clc h e s o l u t i o n o f e g 11..4−− > ” ) ; disp ( ” t he a 1 =. = . 35 35 , a 2 =. = . 16 16 , a 3 =. = . 21 21 , a 4 = .0 .0 1
//1 st
row b 1 =. = . 54 54 , b 2 =. = . 42 42 , b 3 =. = . 54 54 , b 4 =. =. 1 c 1 =. = . 04 04 , c 2 =. = . 24 24 , c 3 =. =. 1 , c 4 =. = . 65 65 d 1 =. = . 07 07 , d 2 =. = . 18 18 , d 3 =. = . 15 15 , d 4 = .2 .2 4 r 1 =1 = 1 4 , r 2 =2 =2 8 , r 3 =1 = 1 7. 7. 5 , r 4 = 10 10 .5 .5
5 6 7 8
/ / 2 n d r ow ow / / 3 r d row / / 4 t h row // giv en
values 9 10 b4=b4-(b1/a1)*a4
// f o r
making b1= b1=00 11 12 13 14 15 16
b3=b3-(b1/a1)*a3 b2=b2-(b1/a1)*a2 r2=r2-(b1/a1)*r1 b1=b1-(b1/a1)*a1
// f o r
c4=c4-(c1/a1)*a4
m a ki k i n g c 1 =0 =0 17 18 19 20 21 22
c3=c3-(c1/a1)*a3 c2=c2-(c1/a1)*a2 r3=r3-(c1/a1)*r1 c1=c1-(c1/a1)*a1
// f o r
d4=d4-(d1/a1)*a4
10
making d1= d1=00 23 24 25 26 27 28
d3=d3-(d1/a1)*a3 d2=d2-(d1/a1)*a2 r4=r4-(d1/a1)*r1 d1=d1-(d1/a1)*a1
// f o r
c4=c4-(c2/b2)*b4
m a ki k i n g c 2 =0 =0 29 30 31 32 33
c3=c3-(c2/b2)*b3 r3=r3-(c2/b2)*r2 c2=c2-(c2/b2)*b2
/ / f o r m a ki ki n g
d4=d4-(d2/b2)*b4
d2=0 34 35 36 37 38
d3=d3-(d2/b2)*b3 r4=r4-(d2/b2)*r2 d2=d2-(d2/b2)*b2
/ / f o r m a ki ki ng ng
d4=d4-(d3/c3)*c4
d3=0 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
r4=r4-(d3/c3)*r3 d3=d3-(d3/c3)*c3
B2=r4/d4; D2=(r3-(c4*B2))/c3; B1=(r2-(D2*b3)-(B2*b4))/b2; D1=(r1-(B2*a4)-(D2*a3)-(B1*a2))/a1; RATES o f disp (B2,D2,B1,D1, ” th e v a l u e s o f MOLAR FLOW RA D1 , B1 , D2 , B2 r e s p e c t i v e l y a re re ” );
B=D2+B2; x 1B 1B = ( . 21 2 1 * D 2 + . 01 0 1 * B2 B2 ) / B ; x 2B 2B = ( . 54 5 4 * D 2 + . 1* 1* B 2 ) /B /B ; x 3B 3B = ( . 1* 1* D 2 + . 65 6 5 * B 2 )/ )/ B ; x 4B 4B = ( . 15 1 5 * D 2 + . 24 2 4 * B2 B2 ) / B ; h e c o m p os o s i ti t i o n o f s tr t r e aam m B disp ( x4 x 4 B , x 3 B , x 2 B , x 1 B , ” t he i s ” );
54
11
55 56 57 58 59 60
D=D1+B1; x 1D 1D = ( . 35 3 5 * D 1 + . 16 1 6 * B1 B1 ) / D ; x 2D 2D = ( . 54 5 4 * D 1 + . 42 4 2 * B1 B1 ) / D ; x 3D 3D = ( . 04 0 4 * D 1 + . 24 2 4 * B1 B1 ) / D ; x 4D 4D = ( . 07 0 7 * D 1 + . 18 1 8 * B1 B1 ) / D ; h e c o m p os o s i ti t i o n o f s tr t r e aam m D disp ( x4 x 4 D , x 2 D , x 2 D , x 1 D , ” t he i s ” );
Scilab code Exa 1.5 gauss seidel method
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
clc h e s o l n o f e g 1. 1 .5−− > Ga G a us us s S e i d e l Meth od ” ) ; disp ( ” t he f o r i=1:3,xnew(i)=2,e(i)=1 end x=1e-6 while e ( 1 ) > x & e ( 2) 2) > x & e ( 3 ) > x d o f o r i = 1: 1 : 3 , x o ld ld ( i ) = x n e w ( i ) , e n d xnew(1)=(44-xold(2)-2*xold(3))/10 xnew(2)=(-2*xnew(1)+51-xold(3))/10 xnew(3)=(-2*xnew(2)-xnew(1)+61)/10 f o r i=1:3,e(i)= a b s (xnew(i)-xold(i)) end end disp ( ” t h e v a l u e s o f x 1 , x 2 , x 3 r e s p e c t i v e l y i s ” ) ; f o r i=1:3, disp (xnew(i)); end
Scilab code Exa 1.6 gauss seidel method
1 2 3 4
clc h e s o l n o f e g 1. 1 .6−− > Ga G a us us s S e i d e l Meth od ” ) ; disp ( ” t he f o r i=1:3,xnew(i)=2,e(i)=1 end
12
5 6 7 8 9 10 11 12 13 14 15 16
x=1e-6 while e ( 1 ) > x & e ( 2) 2) > x & e ( 3 ) > x d o f o r i = 1: 1 : 3 , x o ld ld ( i ) = x n e w ( i ) , e n d xnew(1)=(9-xold(2)+2*xold(3))/3 xnew(2)=(xnew(1) -8+3*xold(3))/4 -8+3*xold(3))/4 xnew(3)=(xnew(2)-xnew(1)+1)/4 f o r i=1:3,e(i)= a b s (xnew(i)-xold(i)) end end disp ( ” t h e v a l u e s o f x 1 , x 2 , x 3 r e s p e c t i v e l y f o r i=1:3, disp (xnew(i)); end
13
i s ” );
Chapter 2 NONLINEAR ALGEBRAIC EQUATIONS
Scilab code Exa 2.1 algebraic equations
1 2 3 4 5 6 7 8
/ / ch ch 2 e x 2 . 1 s o l v i n g u s i n g n e w t o n ’ s m et h o d .
clc h e s o l n o f e q n 2. 2 .1−− > Newton Newton Method” Method” ) ; disp ( ” t he // i n i t i a l va lu e x=.5 xnew=0 e =1 while e > 10 1 0 ^ - 4 d o x = xn x n ew ew , function y = F a ( x ) , // defi nin g y=x^3-5*x+1;
fn 9 endfunction 10 der= derivative ( F a , x ) ,
//
d i f f e r e n ti a t i n g the fn 11 12 13 14
xnew=x-Fa(x)/der, e = abs ( x n e w - x ) , end he r o o t o f t h e eq n i s ” ) ; disp (xnew, ” t he
14
Scilab code Exa 2.2 algebraic equations
1 clc 2 disp ( ” t h e s o l u t i o n o f e x 2 . 2 −−> P r es e s s u r e Drop i n P i p e ”); 3 meu=1.79*10^-5 4 rough=.0000015 // roug hne ss 5 dia=.004 6 e_by_D=rough/dia 7 rho=1.23 8 v=50 // v e l o c i t y o f a i r 9 l =1 10 Re=(rho*v*dia)/meu // Reynold ’ s num numbe berr 11 ffnew=0.01 12 e =1 13 t1=e_by_D/3.7 / / t er er m 1 o f e q n . 14 t2=2.51/Re / / t er er m 2 o f e q n . 15 disp (Re, ” t h e R ey e y no n o l ds ds no . i s ” ) ; 16 funcprot ( 0 ) 17 while e > 1e 1 e - 6 d o f f = ff ff ne ne w , function y = F h ( f f ) , 18 t3 = sqrt ( f f ) , 19 y=1/t3+2* l o g (t1+t2/t3)/2.3,
/ / d i v i d e by 2 . 3 s i n c e l o g i s b a s e ’ e ’ i n s t e a d o f 10 20 21 22 23 24 25 26
endfunction ; fdash= derivative (Fh,ff); ffnew=ff-Fh(ff)/fdash; e = a b s (ff-ffnew) end disp (ff, ” t h e f a n n i n g f r i c t i o n f a c t o r delta_p=(ff*l*v^2*rho)/(2*dia)
// f ’ ( f f )
is ”)
p r e s s u r e d ro ro p 27 disp (del ta_p , ” t he h e p r e s s u re re drop i n p a s ca l s i s ”);
Scilab code Exa 2.3 algebraic equations
15
//
1 clc 2 disp ( ” t h e s o l u t i o n
fluidization 3 4 5 6 7 8 9 10 11 12 13
o f e g 2 .3 .3 v e l o ci t y ” ); / / g i v e n d a ta ta
−−>
minimum
P=2*101325 T=298.15 M=28.97*10^-3 R=8.314 rho=(P*M)/(R*T) rho_p=1000 dia=1.2*10^-4 / / v o id id f r a c t i o n ep=.42 sph=.88 meu=1.845*10^-5 t1=1.75*rho*(1-ep)/(sph*dia*ep^3)
/ / t h e s e a r e t he h e t er e r ms m s o f t he he function . 14 15 16 17 18 19
t2=150*meu*(1-ep)^2/(sph^2*dia^2*ep^3) t3=(1-ep)*(rho_p-rho)*9.8 vnew=0.1 e1=1 while e1 e 1 > 1 e - 6 d o v = vn v n ew ew , function y=Fb(v); // defi nin g y=t1*v^2+t2*v-t3,
fn 20 endfunction , endfunction , 21 vdash= derivative ( F b , v ) ,
//
d i f f e r e n ti a t i n g the fn 22 23 24 25
vnew=v-Fb(v)/vdash, e1 = abs ( v n e w - v ) , end h e m in im um f l u i d i s a t i o n disp ( v , ” t he );
Scilab code Exa 2.4 algebraic equations
1 clc
16
v e l o c i t y i n m/ m/ s i s ”
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
h e s o l n o f e g 2. 2 .4−− > T er e r m in in a l V e l o c i t y ” ) ; disp ( ” t he dia=2*10^-3 / / g i v e n d a ta ta P=101325 T=298.15 M=28.89*10^-3 R=8.314 rho=(P*M)/(R*T) rho_oil=900 meu=1.85*10^-5 Re_oil_by_v=rho*dia/meu v n e w = 0. 0. 1 , e = 1 while e > 1 e - 6 d o v = vn v n ew ew , R e_ e _ oi o i l = R e _o _ o i l_ l _ b y_ y_ v * v , Cd=24*(1+.15*Re_oil^.687)/Re_oil, vnew= sqrt (4*(rho_oil (4*(rho_oil -rho) -rho) *9.81*dia/(3*Cd*rho)) , e = abs ( v n e w - v ) , end h e t e r m i n a l v e l o c i t y i n m/ m/ s i s ” ) ; disp ( v , ” t he
Scilab code Exa 2.5 algebraic equations
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
clc h e s o l n o f e g 2. 2 .5−− > no n o n l i n e a r e q u a t io io n s ”); disp ( ” t he x n e w = 0. 0. 1 , y n e w = 0 .5 .5 , e 1 =1 =1 , e 2 = 1 while e1 e 1 > 1 e - 6 & e 2 > 1 e - 6 d o x = xn x n ew ew , y = yn yn ew ew , y 1 = e x p (x)+x*y-1, d_fx= e x p ( x ) + y d_fy=x y 2 = s i n (x*y)+x+y-1, d_gx=y* c o s (x*y)+1 d_gy=x* c o s (x*y)+1 t 1 = ( y 2 * d _ fy f y ) , t 2 = ( y1 y1 * d _ g y ) , D1=d_fx*d_gy-d_fy*d_gx delta_x=(t1-t2)/D1 t3=(y1*d_gx),t4=(y2*d_fx) delta_y=(t3-t4)/D1
17
16 17 18 19 20 21 22
xnew=x+delta_x ynew=y+delta_y e 1 = a b s (x-xnew) e 2 = a b s (y-ynew) end a n d y r e s p e c t i v e l y a re re ” ) ; disp (y,x, ” t h e v a l u e s o f x an al ue o f x can be c o ns i d e r e d a s disp ( ” s u c h s m a l l v al ze ro . ”)
Scilab code Exa 2.6 algebraic equations
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
clc h e s o l n o f e g 2. 2 .6−− > ” ) ; disp ( ” t he x n e w = 0. 0. 1 , y n e w = 0 .5 .5 , e 1 =1 =1 , e 2 = 1 while e 1 > 10 10 ^ -6 - 6 & e 2 > 10 1 0 ^ -6 - 6 d o x = xn xn ew ew , y = yn yn e w , y1=3*x^3+4*y^2-145, d_fx=9*x^2 d_fy=8*y y2=4*x^2-y^3+28, d_gx=8*x d_gy=-3*y^2 D2=d_fx*d_gy-d_gx*d_fy delta_x=(y2*d_fy-y1*d_gy)/D2 delta_y=(y1*d_gx-y2*d_fx)/D2 xnew=x+delta_x ynew=y+delta_y e 1 = a b s (xnew-x) e 2 = a b s (ynew-y) end h e v a lu l u e s o f x a nd nd y a re re r e s p e c t i v e l y ” ); disp (y,x, ” t he
18
Chapter 3 CHEMICAL ENGINEERING THERMODYNAMICS
Scilab code Exa 3.1 thermodynamics
clc h e s o l n o f e g 3. 3 .1−− > C ub u b ic ic e q n . o f s t a t e ” ) ; disp ( ” t he mo l ”) ; disp ( ” a l l v a l u e s i n m3 / mo T = 37 3 7 3. 3 . 15 1 5 , P = 10 1 0 13 1 3 25 2 5 , T c = 64 64 7. 7 . 1 , P c = 22 2 2 0. 0 . 55 5 5 * 10 10 ^5 ^5 , w = .3 . 3 4 5 , R = 8 . 3 14 14 5 T r = T / Tc // redu ced Tc , P r = P / Pc Pc , 1 2 3 4
p r e s s u r e & r e du d u c ed e d t e m pe p e r a tu tu r e 6 7 8 9 10 11
12 13 14 15 16
b0=.083-.422*Tr^-1.6 b1=.139-.172*Tr^-4.2 B=(b0+w*b1)*R*Tc/Pc vnew=1 e 1 = 1 , v o ld ld = R * T / P + B o l um um e i n disp (vold, ” t h e s o l n b y v i r i a l g a s e q n . o f v ol m3 / m ol ol b y Z ( T , P) P) i s ” ) ; while e 1 > 1e 1 e - 6 d o v ol o l d = vn v n ew ew , function y = F h ( v o l d ) , y=P*vold/(R*T)-1-B/vold endfunction ; ydash= derivative (Fh,vold); vnew=vold-Fh(vold)/ydash;
19
17 e 1 = a b s (vold-vnew) 18 e n d 19 disp (vold, ” t h e s o l n b y v i r i a l g a s e q n . o f v ol o l um um e i n m3 / m ol ol b y Z ( T ,V ,V) i s ” ) ; 20 / / by b y p en en g r o b i n s o n m e t h o d 21 k=.37464+1.54226*w-.26992*w^2 22 s=1+k*(1-Tr^.5) 23 lpha=s^2 24 a=.45724*R^2*Tc^2*lpha/Pc 25 b=.0778*R*Tc/Pc 26 vnew=b,e2=1, 27 v o l = b , f e =0 =0 , f d =0 =0 28 disp ( ” t he h e v o l u m e o f s a t ur u r a t ed e d l i q . a n d s a t ur u r a t ed ed v a po p o u r b y p en en g− r o b in i n s o n method r e s p e c t i v e l y i s ” ) 29 f o r i = 0: 0: 2 d o 30 vol=vnew 31 y2=vol^3*P+vol^2*(P*b-R*T)-vol*(3*P*b^2+2*b* R*T-a)+P*b^3+b^2*R*T-a*b 32 ydash2=3*P*vol^2+(P*b-R*T)*2*vol-(3*P*b^2+2*b*R*T-a) 33 vnew=vol-y2/ydash2, 34 e 2 = a b s (vol-vnew) 35 i f i = =1 =1 & e 2 < 1 e -6 -6 then fd=vnew,vnew=R*T/P, 36 e n d 37 e n d 38 disp (vol,fd); 39 funcprot ( 0 ) 40 / / by by r e d l i c h −kw kwong ong method 41 i =0 42 a=.42748*R^2*Tc^2.5/Pc 43 b=.08664*R*Tc/Pc 44 Vnew=b,e3=1 45 disp ( ” t he h e v o l o f s a t ur u r a t ed e d l i q . a n d v a po po u r by r e d l i c h k w o n g m et et h o d r e s p e c t i v e l y a r e ” ) ; 46 f o r i = 0: 0: 2 d o V = Vn V n ew ew , function y3=gh(V), 47 y3=V^3*P-R*T*V^2-V*(P*b^2+b*R*T-a/ sqrt ( T ))-a*b/ sqrt ( T ) 48 endfunction , endfunction , 49 deriv= derivative (gh,V);
20
50 51 52
Vnew=V-gh(V)/deriv; e3 = a b s (Vnew-V) if i = =1 =1 & e 3 < 1 e -6 -6 then de=Vnew,Vnew=R*T/P
/ / f o r s a t ur u r a t ed ed l i q . 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
end end disp ( V , d e ) ;
/ / v a n d e r w a a l s m et et ho ho d
71
i =0 a=27*R^2*Tc^2/(64*Pc) b=R*Tc/(8*Pc) vnew=b, v=b,e=1 f o r i = 0: 0: 2 d o v = vn v n ew ew , function v3=bh(v), v3=v^3*P-v^2*(P*b+R*T)+a*v-a*b, endfunction deriva= derivative ( b h , v ) , vnew=v-bh(v)/deriva e4 = a b s ( v - v n e w ) , if i = =1 = 1 & e 4 < 1 0^ 0^ - 6 then sol=vnew,vnew=R*T/P end end h e v o l o f s a t ur u r a t ed e d l i q . a n d v a po p o u r b y v a nd nd e r disp ( ” t he w a al al m e th o d r e s p e c t i v e l y a r e ” ) ; disp (vnew,sol);
Scilab code Exa 3.2 thermodynamics
1 clc 2 disp ( ” t he h e s o l n o f e g 3. 3 .2−− > BPT & DPT Ca lc . ” ) ; 3 z1=.5,P=101.325
/ / g i v e n d a ta ta 4 5 6 7
a 1 = 1 4. 4 . 3 91 91 6 , b 1 = 2 79 7 9 5 .3 .3 2 , c 1 = 2 30 30 . 00 00 2 , a2=16.262,b2=3799.887,c2=226.346 x 1 = z1 z1 , x 2 =1 =1 - z 1 T1sat=b1/(a1- l o g (P))-c1
21
8 9 10 11 12 13 14 15 16 17 18 19 20
T2sat=b2/(a2- l o g (P))-c2 Told=273,e=1 Tnew=x1*T1sat+x2*T2sat while e > 1 e -6 -6 d o T ol o l d = Tn T n ew ew , function y1=fw(Told), P1sat= e x p (a1-b1/(Told+c1)), P2sat= e x p (a2-b2/(Told+c2)), y1=P-x1*P1sat-x2*P2sat endfunction ydashd= derivative (fw,Told) Tnew=Told-fw(Told)/ydashd e = a b s (Told-Tnew) end h e b ub u b bl b l e p o in i n t temp . i n C e l s i u s i s ” disp (Tnew, ” t he );
//
21
calc o f dew de w point 22 y 1 = z 1 , y 2 =1 =1 - z 1 , e = 1 23 Tnew=y1*T1sat+y2*T2sat 24 while e > 1 e -6 -6 d o T ol o l d = Tn T n ew ew , function y11=fw1(Told), 25 P1sat= e x p ( a 1 - b 1 / ( T o l d + c 1 ) ) , 26 P2sat= e x p ( a 2 - b 2 / ( T o l d + c 2 ) ) , 27 y11=1/P-y1/P1sat-y2/P2sat 28 endfunction 29 ydashd= derivative (fw1,Told) 30 Tnew=Told-fw1(Told)/ydashd 31 e = a b s (Told-Tnew) 32 end 33 disp (Tnew, ” t h e d e w p o i n t t e m p . i n C e l s i u s i s ” ) ;
22
Scilab code Exa 3.3 flash calculations using Raoult law
1 clc 2 disp ( ” t he h e s o l n o f e g 3. 3 .3−− > F la l a sh s h C al a l c . u s i ng n g R ao a o u lt lt s Law” ) 3 psat(1)=195.75,psat(2)=97.84,psat(3)=50.32
/ / g i v e n d a ta ta 4 z(1)=.3,z(2)=.3,z(3)=.4 5 bpp=z(1)*psat(1)+z(2)*psat(2)+z(3)*psat(3)
/ / B ub u b b le l e p o i nt nt p r e s s u r e 6 dpp=1/(z(1)/psat(1)+z(2)/psat(2)+z(3)/psat(3))
/ / de de w p o i nt nt p r e s s u r e 7 disp (dpp,bpp, ” t h e b ub u b bl b l e p o i nt n t p r e s s u r e a n d de d ew p o i n t p r e s s u re r e r e s p e c t i v e l y a re re ” ); 8 P = 10 10 0 , v =1 =1 , v ne n e w =1 = 1 , e =1 = 1 , y 2 =0 =0 , d er er = 0
//
g i v e n p r e s s u r e i s b et e t we w e en e n BPP & DPP 9 f o r i = 1 :3 :3 k ( i ) = p s at at ( i ) / P 10 e n d 11 while e > 1e 1 e - 6 d o v = vn v n ew ew , f o r i=1:3, 12 t1=(1-v+v*k(i)),y2=y2+z(i)*(k(i)-1)/t1, 13 der=der-z(i)*(k(i)-1)^2/t1^2, 14 end 15 vnew=v-y2/der,e= a b s ( v n e w - v ) , 16 end 17 f o r i = 1: 1 : 3 , x ( i ) = z ( i ) /( /( 1 - v + v * k ( i ) ) , y ( i ) = x ( i )* )* k ( i ) 18 end 19 mo l f r c t n o f a ce c e to t o ne ne i n l i q . pha se i s disp ( x ( 1 ) , ” mo ” ); 20 moo l f r c t n o f a c e to t o n e i n v a po p o u r p ha h a se se disp ( y ( 1 ) , ” m i s ” ); 21 disp ( x ( 2 ) , ” m o l f r c t n o f a c e t o n i t r i l e i n l i q . p ha h a se se i s ”); 22 mo l f r c t n o f a c e t o n i t r i l e i n v a p o u r . disp ( y ( 2 ) , ” mo
23
23
24
p ha h a se se i s ”); moo l f r c t n o f n it i t ro ro m mee th t h an an e i n l i q . disp ( x ( 3 ) , ” m p ha h a se se i s ”); m o l f r c t n o f n i tr t r o me m e th t h a ne n e i n v ap a p ou ou r disp ( y ( 3 ) , ” mo p ha h a se se i s ”);
Scilab code Exa 3.4 BPT and DPT calculation using modified raoult law
1 a 1 2 = 4 3 7 . 98 98 * 4 .1 .1 8 6 , a 2 1 = 1 23 2 3 8 *4 * 4 . 1 86 86 , v 1 = 7 6. 6. 9 2 , v 2 = 1 8 .0 .0 7 2
// ca o f BP
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
clc h e s o l n o f e g 3. 3 .4−− > ” ) ; disp ( ” t he t=100 x 1 = .5 . 5 , R = 8 . 3 14 14 a1=16.678,b1=3640.2,c1=219.61 a2=16.2887,b2=3816.44,c2=227.02 x2=1-x1 p1sat= e x p (a1-b1/(c1+t)) p2sat= e x p (a2-b2/(c2+t)) h12=v2* e x p (-a12/(R*(t+273.15)))/v1 h21=v1* e x p (-a21/(R*(t+273.15)))/v2 m=h12/(x1+x2*h12)-h21/(x2+x1*h21) g 1 = e x p ( - l o g (x1+x2*h12)+x2*m) g 2 = e x p ( - l o g (x2+x1*h21)-x1*m) p=x1*g1*p1sat+x2*g2*p2sat disp ( p , ” b o i l i n g p o i n t p r e s s u r e i n k P a i s ” ) ;
// ca
24
o f BP 20 21 22 23 24 25 26
27 28 29 30 31
p = 1 0 1 .3 .3 2 5 , x 1 = . 5 , e = 1 x2=1-x1 t1sat=b1/(a1- l o g (p))-c1 t2sat=b2/(a2- l o g (p))-c2 tnew=x1*t1sat+x2*t2sat while e > 1 0^ 0^ - 4 d o t ol o l d = tn t n ew ew , p1sat= e x p (a1-b1/(c1+told)),p2sat= e x p (a2-b2/(c2+ told)), p1sat=p/(g1*x1+g2*x2*(p2sat/p1sat)) tnew=b1/(a1- l o g (p1sat))-c1, e = a b s (tnew-told) end i n t t em e m pe p e ra r a tu t u re re i n C e l s i u s i s ”) disp (tnew, ” b o i l i n g p o in ;
//
32
cal o f dp p 33 34 35 36 37 38 39 40 41 42 43 44
e 1 =1 = 1 , e 2 =1 =1 , e 3 =1 =1 , p ol ol d = 1 t=100,y1=.5 y2=1-y1 p1sat= e x p (a1-b1/(c1+t)) p2sat= e x p (a2-b2/(c2+t)) g 1 =1 = 1 , g 2 =1 =1 , g 11 11 = 1 , g 22 22 = 1 pnew=1/(y1/(g1*p1sat)+y2/(g2*p2sat)) while e 1 > . 00 0 0 01 0 1 d o p ol o l d = pn p n ew ew , while e2 e 2 > . 0 00 00 1 & e 3 > . 00 0 0 0 1 d o g 1 = g1 g 1 1 , g 2 = g 22 22 , x1=y1*pold/(g1*p1sat) x2=y2*pold/(g2*p2sat) x1=x1/(x1+x2) x2=1-x1
25
45
46
47
48
49
h12=v2* e x p (-a12/(R*(t +273.15)))/v1 h21=v1* e x p (-a21/(R*(t +273.15)))/v2 m=h12/(x1+x2*h12)-h21 /(x2+x1*h21) g11= e x p ( - l o g (x1+x2*h12 )+x2*m) g22= e x p ( - l o g (x2+x1*h21 )-x1*m) e2 = a b s ( g 11 11 - g 1 ) , e 3 = a b s (g22-g2)
50 51 52 53 54 55 56
end pnew=1/(y1/(g1*p1sat)+y2/(g2*p2sat)) e 1 = a b s (pnew-pold) end disp (pnew, ” d e w p o i n t p r e s s u r e
i n k Pa i s ” ) ; // calc dpt
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
p = 1 0 1 .3 .3 2 5 , y 1 = . 5 , e 4 =1 = 1 , e 5 =1 =1 , e 6 = 1 y2=1-y1 t1sat=b1/(a1- l o g (p))-c1 t2sat=b2/(a2- l o g (p))-c2 tnew=y1*t1sat+y2*t2sat g 1 1 = 1 , g 22 22 = 1 while e 4 > . 00 0 0 01 0 1 d o t ol o l d = tn t n ew ew , p1sat= e x p (a1-b1/(c1+told)) p2sat= e x p ( a 2 - b 2 / ( c 2 + t o l d ) ) , while e5 e 5 > . 00 0 0 01 01 & e 6 > .0 . 0 00 0 0 1 d o g 1 =g = g 11 11 , g 2 = g2 g2 2 , x1=y1*p/(g1*p1sat) x2=y2*p/(g2*p2sat) x1=x1/(x1+x2) x2=1-x1 h12=v2* e x p (-a12/(R*(t+273.15)))/v1 h21=v1* e x p (-a21/(R*(t+273.15)))/v2
26
72
73
74
m=h12/(x1+x2*h12)-h21 /(x2+x1*h21) g11= e x p ( - l o g (x1+x2*h12 )+x2*m) g22= e x p ( - l o g (x2+x1*h21 )-x1*m) e5 = a b s ( g 11 11 - g 1 ) , e 6 = a b s (g22-g2)
75 76 77 78 79 80 81 82
end
p1sat=p*(y1/g1+y2*p1sat/(g2*p2sat)) tnew=b1/(a1- l o g (p1sat))-c1 e 4 = a b s (tnew-told) end de w p o in i n t t em e m pe p e ra r a tu t u re re disp (tnew, ” de ;
i n C e l s i u s i s ”)
Scilab code Exa 3.5 flash calculations using modified Raoult law
1 clc 2 disp ( ” t he he s o l n
o f e g 3. 3 .5−− > F la l a sh s h c a l c . u s in in g M o d if i f i e d R ao a o ul u l t l aw aw ” ) ;
3 a 12 1 2 = 2 92 9 2 . 6 6* 6* 4 .1 . 1 8 , a 21 2 1 = 1 4 4 5 . 26 26 * 4 .1 .1 8 , v 1 = 7 4. 4. 0 5* 5* 1 0^ 0^ - 6 , v 2 = 1 8. 8 . 0 7* 7* 1 0^ 0^ - 6 , R = 8 . 3 14 14 4 t=100,z1=.3, 5 z2=1-z1 6 a 1 = 1 4. 4 . 3 91 9 1 5 5 , a 2 = 1 6. 6. 2 62 62 , b 1 = 2 79 7 9 5 .8 .8 2 , b 2 = 3 79 79 9 .8 .8 9 , c 1 = 2 30 30 . 00 00 2 , c 2 = 2 2 6. 6. 3 5 7 e1=1,e2=1,e3=1,e4=1,e5=1,e6=1,vnew=0 8 // ca lc
o f BPP BP P 9 x 1 =z = z 1 , x 2 = z2 z2 10 p1sat= e x p (a1-(b1/(t+c1))) 11 p2sat= e x p (a2-(b2/(t+c2)))
27
12 13 14 15 16 17 18 19
h12=v2* e x p (-a12/(R*(t+273.15)))/v1 h21=v1* e x p (-a21/(R*(t+273.15)))/v2 m=(h12/(x1+x2*h12))-(h21/(x2+x1*h21)) g 1 = e x p ( - l o g (x1+x2*h12)+x2*m) g 2 = e x p ( - l o g (x2+x1*h21)-x1*m) p=x1*g1*p1sat+x2*g2*p2sat h e b ub u b bl b l e p o in in t p r e s s u r e i s ” ); disp ( p , ” t he b pp pp = p , g b1 b1 = g1 g 1 , g b2 b2 = g 2
//g1
& g 2 a re re a c t i v i t y co− e f f i c i e n t s //
20
calc o f DPP DP P 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
y 1 =z = z 1 , y 2 = z2 z2 g 1 = 1 , g 2 =1 =1 pnew=1/(y1/(g1*p1sat)+y2/(g2*p2sat)) g 1 n = g1 g1 , g 2 n = g2 g2 while e 1 > . 00 0 0 01 0 1 d o p ol o l d = pn p n ew ew , while e 2 > . 0 0 01 0 1 & e 3 > . 0 00 00 1 d o g 1 =g = g 1n 1n , g 2 = g2 g2 n , x1=y1*pold/(g1*p1sat) x2=y2*pold/(g2*p2sat) x1=x1/(x1+x2) x2=1-x1 g1n= e x p ( - l o g (x1+x2*h12)+x2*m) g2n= e x p ( - l o g (x2+x1*h21)-x1*m) e2 = a b s ( g1 g 1 n - g 1 ) , e 3 = a b s (g2n-g2) end pnew=1/(y1/(g1n*p1sat)+y2/(g2n*p2sat)) e 1 = a b s (pnew-pold) end h e d e w p o in in t p r e s s u r e i s ” ); disp (pnew, ” t he d p p = p ne ne w , g d 1 = g1 g1 n , g d2 d2 = g 2 n p=200 v=(bpp-p)/(bpp-dpp) g1=((p-dpp)*(gb1-gd1))/(bpp-dpp)+gd1
28
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
g2=((p-dpp)*(gb2-gd2))/(bpp-dpp)+gd2
/ / c a l c o f d i s t r i b u t i o n co− e f f i c i e n t s while e4 e 4 > . 00 00 01 0 1 & e 5 > .0 . 0 00 0 0 1 d o g 1n 1n = g1 g 1 , g 2 n = g2 g2 , k1=g1n*p1sat/p k2=g2n*p2sat/p while e6 e 6 > . 00 00 0 01 01 d o v = vn v n ew ew , function f = f v ( v ) , y1=(k1*z1)/(1-v+v*k1) y2=(k2*z2)/(1-v+v*k2) x1=y1/k1 x2=y2/k2 f=y1-x1+y2-x2 endfunction derv= derivative (fv,v) vnew=v-fv(v)/derv e6 = a b s (vnew-v) end h12=v2* e x p (-a12/(R*(t+273.15)))/v1 h21=v1* e x p (-a21/(R*(t+273.15)))/v2 m=(h12/(x1+x2*h12))-(h21/(x2+x1*h21)) g1 = e x p ( - l o g (x1+x2*h12)+x2*m) g2 = e x p ( - l o g (x2+x1*h21)-x1*m) e4 = a b s ( g1 g1 - g 1 n ) , e 5 = a b s (g2-g2n) end h e n o . o f m ol o l es e s i n v a p ou o u r p ha h a se se i s ” ) ; disp ( v , ” t he x 1 a n d y 1 r e s p e c t i v e l y a re re ” ); disp (y1,x1, ” x1
Scilab code Exa 3.6 vapour pressure calculation using cubic equation of
state 1 clc 2 disp ( ” t he he s o l n
o f e g 3. 3 .6−− > V a p o u r P r e s s u r e u s i n g C u b ic i c Eq E qn . o f S t a te te ” )
3 t = 3 7 3 .1 .1 5 , t c = 64 6 4 7 .1 .1 , p c = 2 2 0. 0 . 5 5* 5 * 1 0 ^5 ^5 , w = . 34 34 5 , R = 8 . 3 1 4
29
// gi ven //
4 f 1 = 1 , e 1 =1 =1 , e 2 =1 =1 , v ne n e w =1 = 1 , p ne ne w = 1
a ss s s um u m ed ed v a l u e s 5 6 7 8 9 10 11
k=.37464+1.54226*w-.26992*w*2 s=(1+k*(1-(t/tc)^.5))^2, a=.45724*R*R*tc*tc*s/pc b=.0778*R*tc/pc
// c a l c o f v ol . o f l i q . while f 1 > 1 0^ 0^ - 4 d o p = pn p n ew ew , v ne n e w =b =b , t1=(p*b-R*t)
/ / c o− e f f i c i e n t s
o f v o l . i n t he he eq n . 12 13 14 15 16 17 18 19 20 21 22 23
t2=3*p*b^2+2*b*R*t-a t3=p*b^3+b^2*R*t-a*b while e 1 >1 > 1 e -6 - 6 & v n e w > 0 d o v ol ol d = vn vn e w , y1=vold^3*p+vold^2*t1-vold*t2+t3, der=3*vold^2*p+2*vold*t1-t2 vnew=vold-y1/der e 1 = a b s (vnew-vold) end vliq=vold,
// f u g a c i t y o f l i q .
zliq=p*vliq/(R*t) c=(a/(2*1.414*b*R*t))*( l o g ((vliq+(1+ sqrt (2))*b)/( vliq+(1- sqrt (2))*b))), 24 t5=zliq-p*b/(R*t) 25 fl2=p* e x p (zliq-1- l o g ( t 5 ) - c ) , 26 vvnew=R*t/p, / / a ss s s um u m ed ed v a l u e c l o s e
t o t h e i d e a l v a lu lu e 27 / / c a l c o f v o l . o f v a p ou ou r 28 29 30 31 32 33 34 35 36
while e 2 > 1e 1 e - 6 d o v vo v o ld l d = v vn v n ew ew , x2=vvold^3*p+vvold^2*t1-vvold*t2+t3, der1=3*vvold^2*p+2*vvold*t1-t2, vvnew=vvold-x2/der1 e 2 = a b s (vvnew-vvold) end
/ / f u g a c i t y o f v ap a p ou ou r vvap=vvold,zvap=p*vvap/(R*t), d=(a/(2* sqrt (2)*b*R*t))*( l o g ((zvap+(1+ sqrt (2))*b*p/(
30
R*t))/(zvap+(1- sqrt (2))*p*b/(R*t)))) 37 t6=zvap-p*b/(R*t) 38 fv2=p* e x p (zvap-1- l o g (t6)-d) 39 pnew=p*fl2/fv2 // updating the
v a l ue ue o f P 40 f 1 = a b s ((fl2-fv2)/fv2) 41 e n d 42 disp ( p , ” t he h e v a po p o ur u r p r e s s u r e o f w at a t er e r i n Pa i s ” ) ;
Scilab code Exa 3.7 pressure x y diagram using gamma phi approach
1 clc 2 disp ( ” t he he s o l n
o f e g 3. 3 .7−− >P−x−y c a l c . u s i n g Gamma− P hi h i a p pr p r o ac ac h ” ) ; 3 e t = 1 , e r = 1 , s o ld //assumed ld = 0 , s n ew ew = 0 values 4 R = 8 . 31 31 4 , t = 1 0 0 , x 1 = .9 . 9 5 8 , a 12 1 2 = 1 0 7 . 38 3 8 * 4 .1 .1 8 , a 21 21 =469.55*4.18, 5 tc1=512.6,tc2=647.1,pc1=80.97*10^5,pc2=220.55*10^5, w1=.564,w2=.345,zc1=.224,zc2=.229,v1=40.73*10^-6, / / g i v e n d a ta ta v2=18.07*10^-6 6 x2=1-x1, 7 a 1 = 1 6. 6 . 5 93 93 8 , a 2 = 1 6. 6. 2 62 62 , b 1 = 3 64 64 4 .3 .3 , b 2 = 3 79 7 9 9 .8 .8 9 , c 1 = 2 39 39 . 76 76 , c 2 = 2 2 6. 6. 3 5 8 p1sat= e x p (a1-b1/(c1+t))*1000
//Saturation Pressure 9 p2sat= e x p (a2-b2/(c2+t))*1000 10 t=t+273.15
//Temp i n K el e l vi v i n r eq eq . 11 12 13 14
h12=(v2/v1)* e x p (-a12/(R*t)) h21=(v1/v2)* e x p (-a21/(R*t)) z=h12/(x1+x2*h12)-h21/(x2+x1*h21) g 1 = e x p ( - l o g (x1+x2*h12)+x2*z)
/ / A c t i v i t y Co− e f f i c i e n t s 31
15 g 2 = e x p ( - l o g (x2+x1*h21)-x1*z) 16 tr1=t/tc1
//Reduced
Temp. 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
b0=.083-.422*tr1^-1.6 b1=.139-.172*tr1^-4.2 b11=(R*tc1/pc1)*(b0+b1*w1) tr2=t/tc2 b0=.083-.422*tr2^-1.6 b1=.139-.172*tr2^-4.2 b22=(R*tc2/pc2)*(b0+b1*w2) w12=(w1+w2)^.5 tc12=(tc1*tc2)^.5 zc12=(zc1+zc2)^.5 v c1 c1 = z c 1 * R * t c1 c 1 / p c1 c 1 , v c2 c2 = z c 2 * R * t c2 c2 / p c 2 vc12=((vc1^.33+vc2^.33)/2)^3 pc12=zc12*R*tc12/vc12 tr12=t/tc12 b0=.083-.422*tr12^-1.6 b1=.139-.172*tr12^-4.2 b12=(R*tc12/pc12)*(b0+b1*w12) d12=2*b12-b11-b22 p=x1*g1*p1sat+x2*p2sat*g2 y 1 = x 1 * g1 g1 * p 1 s a t / p , y 2 = x2 x 2 * g 2 * p 2 s at at / p pnew=p,
/ / c a l c o f P r es es s ur e while et e t > 1 e - 6 d o p = pn p n ew ew , f1=p1sat*( e x p (b11*p1sat/(R*t)))*( e x p ((v1*(pp1sat)/(R*t)))) f2=p2sat*( e x p (b22*p2sat/(R*t)))*( e x p ((v2*(p-p2sat)/( R*t)))) while e r > 1e 1 e - 6 d o s ol o l d = sn s n ew ew , fc1= e x p ((p/(R*t))*(b11+y2^2*d12)) fc2= e x p ((p/(R*t))*(b22+y1^2*d12)) k1=g1*f1/(fc1*p) k2=g2*f2/(fc2*p) snew=x1*k1+x2*k2 y1=x1*k1/snew y2=x2*k2/snew
32
50 51 52 53 54 55 56 57
e r = a b s (snew-sold) end pnew=(x1*g1*f1/fc1)+(x2*g2*f2/fc2) y1=x1*g1*f1/(fc1*pnew) y2=x2*g2*f2/(fc2*pnew) e t = a b s (pnew-p) end a m t . o f m et e t ha h a no n o l i n v ap a p ou o u r p ha h a se se disp (p,y1, ” t h e am s y s t e m p r e s s u r e i n Pa P a r e s p e c t i v e l y a re re ” )
and
Scilab code Exa 3.9 chemical chemical reaction reaction engineering engineering 2 simultane simultaneous ous reacreac-
tions 1 clc 2 disp ( ” t he he s o l n 3 4 5 6 7 8 9
o f e g 3. 3 .9−− > C h em e m i ca ca l R e a c t i o n E q u i l i b r i u m −2 S i m u l ta t a n e ou o u s R e a c t i on on s ” ) / / l e t x 1 a n d x 2 b e t he h e r e a c t i o n c o − o r d i na na t e f o r 1 s t a n d 2 nd nd r e a c t i o n s //assumed x 1 ne n e w = .9 . 9 , x 2 ne n e w = .6 .6 , r 1 = 1 , r 2 = 1 values / / s i n c e P=1 a t m Kp=1 // giv en K 1 = . 57 5 7 4 , K 2 = 2 .2 .2 1 / / a t eqm . K y e 1 = K1 K1 , K y e2 e2 = K 2
while r 1 > 1 e -6 -6 & r 2 > 1 e -6 - 6 , x 1 = x1 x1 ne ne w , x 2 = x2 x 2 ne ne w , m_CH4=1-x1,m_H2O=5-x1-x2,m_CO=x1-x2,m_H2=3*x1+x2, / / m o le l e s o f r e a c t a n t s a n d p r od o d u ct ct s m_ m _ C O 2 = x2
a t e qm . 10 total=m_CO2+m_H2+m_CO+m_H2O+m_CH4 11 Ky1=m_CO*m_H2^3/(m_CH4*m_H2O*total^2) 12 Ky2=m_CO2*m_H2/(m_CO*m_H2O) 13 f1= Ky1 -.574
//1 st
f u n c t i o n i n x1 and x2
f2= Ky2 -2.21
14
//2nd
f u n c t i o n i n x1 and x2 15
d3=((3*x1+x2)^2*(12*x1-8*x2))/((1-x1)*(5-x1-x2)
33
16
17 18
*(6+2*x1)^2) d4=(3*x1+x2)^3*(x1-x2)*(8*x1^2+6*x1*x2-24*x1+2* x2-16) d5=((1-x1)^2)*((5-x1-x2)^2)*((6+2*x1)^3) // df 1 df1_dx1=d3-(d4/d5)
/dx1 − p a r t i a l d e r i v a t i v e o f f 1 w r t t o x 1 19 20 21
22
d6=3*(x1-x2) *((3*x1+x2)^2) -(3*x1+x2)^3 -(3*x1+x2)^3 d7=(1-x1)*(5-x1-x2)*((6+x1*2)^2) d8=((x1-x2)*(3*x1+x2)^3)/((1-x1)*((5-x1-x2)^2) *(6+2*x1)^2) // df 1 df1_dx2=(d6/d7)+d8
/dx2 − p a r t i a l d e r i v a t i v e o f f 1 w r t t o x 2 23 24
d9=(x1-x2)*(5-x1-x2) df2_dx1=3*x2/d9-(x2*(3*x1+x2)*(5-2*x1))/(d9^2)
/ / d f 1 / d x2 x 2− p a r t i a l d e r i v a t i v e o f f 1 wrt t o x2 25 26 27
d10=(3*x1+2*x2)/d9 d11=x2*(3*x1+x2)*(2*x2-5)/(d9^2) df2_dx2=d10-d11
// d f 1 / d x 2− p a r t i a l d e r i v a t i v e o f f 1 w r t t o x 2
28 dm=df1_dx1*df2_dx2 -df1_dx2* df2_dx1 29 delta_x1=(f2*df1_dx2 -f1*df2_dx2)/ dm 30 delta_x2=(f1*df2_dx1 -f2*df1_dx1)/ dm 31 x1new=x1+delta_x1
// updati ng
t h e v a lu l u e s o f x1 & x2 x2new=x2+delta_x2 r 1 = a b s ( x1 x 1 - x 1 n ew ew ) , r 2 = a b s (x2new-x2) end h e v a l u e o f X1 a n d X2 X2 r e s p e c t i v e l y i s ” ) disp (x2,x1, ” t he ; 36 m_CH4=1-x1,m_H2O=5-x1-x2,m_CO=x1-x2,m_H2=3*x1+x2, / / m o le l e s o f r e a c t a n t s a n d p r od o d u ct ct s m_ m _ C O 2 = x2 32 33 34 35
a t e qm . 37 total=m_CO2+m_H2+m_CO+m_H2O+m_CH4 38 disp (m_CO2,m_H2,m_CO,m_H2O,m_CH4, ” t h e m ol o l es e s a t eq e qm o f CH4, H2O,CO,H2, O,CO,H2, CO2 ar e ” ) 39 disp (total, ” t o t a l n u m b e r o f m ol o l es e s a t eq m . i s ” )
34
Scilab code Exa 3.10 adiabatic flame temperature
1 clc 2 disp ( ” t he h e s o l n o f e g 3. 3 .1 0−− > A d i a b a t i c F l a me me Temperature” ); 3 u1=1,u2=3.5,u3=2,u4=3 //moles
g i v e n 1−C2H6 , 2−O2 , 3−CO2, 4−H2 H2O O 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19
a1=1.648,a2=6.085,a3=5.316,a4=7.7 b1=4.124e-2,b2=.3631e-2,b3=1.4285e-2,b4=.04594e-2 c1=-1.53e-5,c2=-.1709e-5,c3=-.8362e-5,c4=.2521e-5 d1=1.74e-9,d2=.3133e-9,d3=1.784e-9,d4=-.8587e-9 n1=1,n2=4,n3=10,n4=0,t0=298.15,t1=25,e1=1 t1=t1+273.15
/ / c a l c . o f s um o f c o − e f f i c i e n t s o f h e a t c a p a c i t y o f t h e r xn xn . sa=n1*a1+n2*a2+n3*a3+n4*a4 sb=n1*b1+n2*b2+n3*b3+n4*b4 sc=n1*c1+n2*c2+n3*c3+n4*c4 sd=n1*d1+n2*d2+n3*d3+n4*d4 da=u4*a4+u3*a3-u2*a2-u1*a1 db=u4*b4+u3*b3-u2*b2-u1*b1 dc=u4*c4+u3*c3-u2*c2-u1*c1 dd=u4*d4+u3*d3-u2*d2-u1*d1 h0=(u4*(-57.798)+u3*(-94.05)-u2*0-u1*(-20.236))*1000
/ / e n t h al a l p y o f t h e r xn xn . 20 tnew=1000 21 while e1 e 1 > 1 e - 6 d o t = tn t n ew ew , 22 function f = f t ( t ) , 23 f=sa*(t-t1)+(sb/2)*(t^2-t1^2)+(sc/3)*(t^3-t1^3) +(sd/4)*(t^4-t1^4)+h0+da*(t-t0)+(db/2)*(t^2t0^2)+(dc/3)*(t^3-t0^3)+(dd/4)*(t^4-t0^4) 24 endfunction 25 d r = derivative ( f t , t ) , 26 tnew=t-ft(t)/dr,
35
27 e 1 = a b s ((tnew-t)/tnew) 28 e n d 29 disp (tnew, ” t he h e a d i a b a t i c f la la m mee t e m p i n K i s ” ) ;
36
Chapter 4 INITIAL VALUE PROBLEMS
Scilab code Exa 4.1 solution of ordinary differential equation
1 clc 2 disp ( ” t he he s o l u t i o n
o f e . g . 4 . 1 −−> I n t e g r a t i o n o f O r di d i na n a ry r y D i f f e r e n t i a l E qu q u at a t io io n ” ) 3 / / i n t h i s p r o bl b l e m d y / dx d x=x+y 4 x_0=0 // i n i t i a l va lu es giv en 5 6 7 8 9 10 11 12
y_0=0 function ydash=fs(x,y), ydash=x+y, endfunction f o r x_0=0:0.1:0.2, h=0.1
// st ep
i n c re r e me m e nt nt o f 0 . 1 13 14 15 16 17 18 19 e n d
f_0=fs(x_0,y_0) k1=h*f_0 k2=h*fs(x_0+h/2,y_0+k1/2) k3=h*fs(x_0+h/2,y_0+k2/2) k4=h*fs(x_0+h,y_0+k3) y_0=y_0+(k1+2*k2+2*k3+k4)/6
37
20 y_0=y_0-(k1+2*k2+2*k3+k4)/6
//to
g e t v a lu l u e a t x = 0. 0. 2 21 disp (y_0, ” t h e v a lu l u e o f y a t x =. = . 2 u s i ng n g Runge K u t t a m et et ho ho d i s ” ) ; 22 a e = e x p (x_0)-x_0-1
/ / a n a l y t i c a l e q n g i ve ve n 23 disp (ae, ” t he h e v al a l ue u e o f y a t x =0 =0 . 2 f r o m a n a l y t i c a l eqn i s ” );
Scilab code Exa 4.2 solution of ordinary differential equation
1 clc 2 disp ( ” t he he s o l u t i o n
o f e . g . 4 . 2 −−>O r d i n a r y D i f f e r e n t i a l E qn . − R un un ge ge K u t ta t a m et et ho ho d ” ) 3 / / i n t h i s p ro r o bl b l em e m d y / dx dx=−y / ( 1 + x ) 4 x_0=0 // i n i t i a l va lu es giv en 5 6 7 8 9 10
y_0=2 function ydash=fr(x,y), ydash=-y/(1+x), endfunction f o r x_0=0:0.01:2.5, h=0.01
// st ep
i nc n c re r e me m e nt nt o f 0 . 01 11 f_0=fr(x_0,y_0) 12 k1=h*f_0 13 k2=h*fr(x_0+h/2,y_0+k1/2) 14 k3=h*fr(x_0+h/2,y_0+k2/2) 15 k4=h*fr(x_0+h,y_0+k3) 16 y_0=y_0+(k1+2*k2+2*k3+k4)/6 17 e n d 18 y_0=y_0-(k1+2*k2+2*k3+k4)/6
//
f i n a l v al a l ue u e a t x =2 =2 . 5 19 disp (y_0, ” t h e v a lu l u e o f y a t x =2 = 2 . 5 u s i ng n g Ru nge K u t t a m et et ho ho d i s ” ) ;
38
Scilab code Exa 4.3 double pipe heat exchanger
1 clc 2 disp ( ” t he h e s o l u t i o n o f e . g . 4 . 3 −−>D o ub u b l e P i p e H ea ea t E x c h a n g e r ”); 3 r h o = 1 00 // 00 0 , v = 1 , d i a = 2. 2 . 4 *1 *1 0 ^ - 2 , C p = 4 18 18 4
g i v e n d at at a 4 5 6 7 8 9 10 11 12 13 14 15
mdot=rho*v*%pi*dia^2/4 t1=mdot*Cp U=200 Ts=250 // i n i t i a l z =0
va lu es giv en / / dT dT// dz= dz=U U∗ p i ∗ d i a ∗ ( Ts−T) /( mdot ∗ Cp ) function Tgrad=fr(z,T), Tgrad=U*%pi*dia*(Ts-T)/(mdot*Cp), endfunction T=20 f o r z=0:0.01:10, h=0.01
// st ep
i nc n c re r e me m e nt nt o f 0 . 01 16 17 18 19 20 21 22
k1=h*fr(z,T) k2=h*fr(z+h/2,T+k1/2) k3=h*fr(z+h/2,T+k2/2) k4=h*fr(z+h,T+k3) T=T+(k1+2*k2+2*k3+k4)/6 if z==5 then T=T-(k1+2*k2+2*k3+k4)/6, h e v al a l ue u e o f T i n d e g C e l s i u s a t z =5 =5 disp ( T , ” t he m u s i n g R u n g e K ut u t ta ta m e t h o d i s ” ) ; end
23 24 e n d 25 T=T-(k1+2*k2+2*k3+k4)/6
// f i n a l
v a l u e a t z =1 =1 0 m 26 disp ( T , ” t he h e v al a l u e o f T i n d e g C e l s i u s a t z =1 0 m u s i n g R u n g e K ut u t ta t a m et et ho ho d i s ” ) ; 39
Scilab code Exa 4.4 stirred tank with coil heater
1 clc 2 disp ( ” t he he s o l u t i o n o f e . g . C o i l H ea e a te te r ” ) 3 vol=.5*.5*2 4 rho=1000 5 m=vol*rho 6 vol_rate=.001 7 mdot=vol_rate*rho 8 out_A=1 9 U=200 10 Cp=4184 11 T 1 = 20 20 , T s = 25 25 0 , T _ e xi xi t = 2 8
4 .4
−−>
S t i r r e d T a n k w it it h / / g i v e n d a ta ta
/ / te tem p i n C e l s i u s 12 t1=(mdot*Cp*T1+U*out_A*Ts)/(m*Cp)
/ / t er e r ms m s o f t he h e f u n c ti ti o n 13 14 15 16 17 18
t2=(mdot*Cp+U*out_A)/(m*Cp)
/ / d t / d t = ( m d o t ∗ Cp(T1−T)+ T)+Q Q dot ) /m∗ Cp function tgrad=fv(tm,T), tgrad=t1-t2*T endfunction T=20
//
i n i t i a l va lu e 19 funcprot ( 0 ) 20 f o r tm=0:1:5000, 21 h =1
// st ep
i n c re r e me m e nt nt o f 1 s e c 22 23 24 25 26 27
k1=h*fv(tm,T) k2=h*fv(tm+h/2,T+k1/2) k3=h*fv(tm+h/2,T+k2/2) k4=h*fv(tm+h,T+k3) e1 = a b s (T-T_exit) h e t i me me a t w h i c h e x i t if e1<1e-3 t he h e n d i sp sp (tm, ” t he
40
temp . i n s e c . i s 2 8 C i s ” ) 28 end 29 T=T+(k1+2*k2+2*k3+k4)/6 30 end 31 T=T-(k1+2*k2+2*k3+k4)/6
// f i n a l
s t e a dy d y s t a t e temp . 32 disp ( T , ” t he h e v a l u e o f s te t e a d y Te T emp i n C e l s i u s i s ” ) ;
Scilab code Exa 4.5 stirred tank with coil heater
1 clc 2 disp ( ” t he he s o l u t i o n
o f e . g . 4 .5 w iitt h C o i l H ea e a te te r ” );
−−>
S t i r r e d V e s s el el //
3 m=760
g i v e n d at at a 4 5 6 7 8 9 10
mdot=12/60 U_into_A=11.5/60 Cp=2.3 T 1 = 25 2 5 , T s = 1 50 50 t1=(mdot*Cp*T1+U_into_A*Ts)/(m*Cp) t2=(mdot*Cp+U_into_A)/(m*Cp)
12 13 14 15 16 17
function tgrade=fg(t,T); tgrade=(t1-t2*T), endfunction T=25 f o r t=0:1:3000, h =1
/ / u s i n g e n e r g y b a l a n c e we k n o w a c c u m u l a ti t i o n= n = i n pu pu t − output 11 / /T /T i s t h e t e m p . o f f l u i d i n s t i r r e d t a n k
// st ep
i n c re r e me m e nt nt o f 1 s e c 18 19 20 21
k1=h*fg(t,T) k2=h*fg(t+h/2,T+k1/2) k3=h*fg(t+h/2,T+k2/2) k4=h*fg(t+h,T+k3)
41
22 T=T+(k1+2*k2+2*k3+k4)/6 23 end 24 T=T-(k1+2*k2+2*k3+k4)/6
/ / t o g et et
v a l ue u e a t x =0 =0 . 2 25 disp ( T , ” t he h e v al al u e o f T i n d eg C a f t e r 50 mins i s ” ); 26 T_steady=(mdot*Cp*T1+U_into_A*Ts)/(mdot*Cp+U_into_A) 27 disp (T_ steady , ” t he h e s t e aad d y s t a t e temp i n d e g C i s ” ) ;
Scilab code Exa 4.6 pneumatic conveying
1 2 3 4 5 6 7 8 9 10 11 12 13 14
clc h e s o l n o f e g 4. 4 .6−− > P n e um u m a t ic i c C o n v ey ey i n g ” ) disp ( ” t he / / g i v e n d a ta ta dia=3*10^-4 v_sprfcl=12 rho_p=900 meu=1.8*10^-5 P=101325 T=298.15 R=8.314 M=28.84*10^-3 rho_air=P*M/(R*T) proj_A=%pi*(dia^2)/4 volm=%pi*(dia^3)/6 t1=rho_air*proj_A/(volm*rho_p)
/ / t er e r m s o f t he he f u n c t i o n 15 16 17 18 19 20 21 22 23 24
t2=((rho_air/rho_p)-1)*9.81*2 y =0 f o r z=.01:.01:10, h=.01 vn1= sqrt ( y ) Re=rho_air*(12-vn1)*dia/meu Cd=24*(1+.15*Re^.687)/Re function dy_by_dz=fy(z,y), dy_by_dz=t1*Cd*(12- sqrt (y))^2+t2, endfunction
42
25 kk1=h*fy(z,y) 26 kk2=h*fy(z+h/2,y+kk1/2) 27 kk3=h*fy(z+h/2,y+kk2/2) 28 kk4=h*fy(z+h,y+kk3) 29 y=y+(kk1+2*kk2+2*kk3+kk4)/6 30 end 31 v = sqrt ( y ) // f i n a l
v al a l ue ue o f
velocity 32 disp ( v , ” t he h e v a l u e o f v a t t he h e e n d o f t he h e p n eu e u m at at i c c o nv n v e yo yo r i s ” ) ;
Scilab code Exa 4.7 simultaneous ordinary differential equations
1 2 3 4 5 6 7 8 9 10 11 12
clc h e s o l n o f e g 4. 4 .7−− > S i m u l t a n e o u s O . D. D. E . ” ) disp ( ” t he function dx_dt=fw(t,x,y); dx_dt=x+2*y, endfunction function dy_dt=fq(t,x,y); dy_dt=3*x+2*y endfunction // i n i t i a l va lu es y=4,x=6
/ / s o l v i n g b y R u n g e−K u t ta t a m et et ho ho d f o r t=0:.1:.2, h=.1
// st ep
i n c re r e me m e nt nt o f 0 . 1 13 14 15 16 17 18 19 20 21
k1=h*fw(t,x,y) l1=h*fq(t,x,y) k2=h*fw(t+h/2,x+k1/2,y+l1/2) l2=h*fq(t+h/2,x+k1/2,y+l1/2) k3=h*fw(t+h/2,x+k2/2,y+l2/2) l3=h*fq(t+h/2,x+k2/2,y+l2/2) k4=h*fw(t+h,x+k3,y+l3) l4=h*fq(t+h,x+k3,y+l3) x=x+(k1+2*k2+2*k3+k4)/6
43
22 23 24 25 26 27 28 29 30
y=y+(l1+2*l2+2*l3+l4)/6 end x=x-(k1+2*k2+2*k3+k4)/6 y=y-(l1+2*l2+2*l3+l4)/6 re ” ); disp (y,x, ” t h e v a l u e s o f x a n d y r e p e c t i v e l y a re t_an=.2 x_an=4* e x p (4*t)+2* e x p ( - t ) y_an=6* e x p (4*t)-2* e x p ( - t ) re disp (y_an,x_an, ” t h e a n a l y t i c a l v a l u e s o f x a n d y a re r e s p e c t i v e l y ” );
Scilab code Exa 4.8 simultaneous ordinary differential equations
1 2 3 4 5 6 7 8 9 10 11
clc h e s o l n o f e g 4. 4 .8−− > S i m u l t a n e o u s O . D. D. E . ” ) disp ( ” t he function dy_dx=fw(x,y,z); dy_dx=z, endfunction function dz_dx=fq(x,y,z); dz_dx=-y endfunction // i n i t i a l va lu es y=2,z=1 f o r x=0:.1:3, // st ep h=.1
i n c re r e me m e nt nt o f 0 . 1 12 13 14 15 16 17 18 19 20 21
k1=h*fw(x,y,z) l1=h*fq(x,y,z) k2=h*fw(x+h/2,y+k1/2,z+l1/2) l2=h*fq(x+h/2,y+k1/2,z+l1/2) k3=h*fw(x+h/2,y+k2/2,z+l2/2) l3=h*fq(x+h/2,y+k2/2,z+l2/2) k4=h*fw(x+h,y+k3,z+l3) l4=h*fq(x+h,y+k3,z+l3) y=y+(k1+2*k2+2*k3+k4)/6 z=z+(l1+2*l2+2*l3+l4)/6
44
22 23 24 25 26
end y=y-(k1+2*k2+2*k3+k4)/6 z=z-(l1+2*l2+2*l3+l4)/6 re ” ) ; disp (z,y, ” t h e v a l u e s o f y a n d z r e s p e c t i v e l y a re
27 28 29 30 31
/ / f o r t he h e g iv i v en e n a n a l y t i c a l e q n s t h e v a l u es es o f A and a l ph p h a c a n b e d e ttee r mi m i n ed ed u s i n g i n i t i a l v a l u e s o f y and z alpha= atan ( 2 ) A=2/ s i n (alpha) y_an=A* s i n (alpha+x) z_an=A* c o s (alpha+x) re disp (z_an,y_an, ” t h e a n a l y t i c a l v a l u e s o f y a n d z a re ”);
Scilab code Exa 4.9 simultaneous ordinary differential equations
1 clc 2 disp ( ” t he h e s o l n o f e g 4. 4 .9−− > S i m u l t a n e o u s O . D. D. E . ” ) 3 function dy_dx=fw(x,y,z); // l e t u s
h a ve v e d y / dx d x=z , t h e r e f o r e d 2y 2y / d x2 x 2=d z / dx dx 4 5 6 7 8 9 10 11
dy_dx=z, endfunction function dz_dx=fq(x,y,z); dz_dx=-y*x, endfunction y=2,z=1 f o r x=0:.1:3, h=.1
// st ep
i n c re r e me m e nt nt o f 0 . 1 12 13 14 15 16 17
k1=h*fw(x,y,z) l1=h*fq(x,y,z) k2=h*fw(x+h/2,y+k1/2,z+l1/2) l2=h*fq(x+h/2,y+k1/2,z+l1/2) k3=h*fw(x+h/2,y+k2/2,z+l2/2) l3=h*fq(x+h/2,y+k2/2,z+l2/2)
45
18 k4=h*fw(x+h,y+k3,z+l3) 19 l4=h*fq(x+h,y+k3,z+l3) 20 y=y+(k1+2*k2+2*k3+k4)/6 21 z=z+(l1+2*l2+2*l3+l4)/6 22 end 23 y=y-(k1+2*k2+2*k3+k4)/6 24 z=z-(l1+2*l2+2*l3+l4)/6 25 disp (z,y, ” t h e v a l u e s o f y an a n d z r e p e c t i v e l y a re re ” );
Scilab code Exa 4.10 series of stirred tank with coil heater
1 clc 2 disp ( ” t h e s o l u t i o n
o f eg 4 .1 0 T a n k s w it i t h C o i l H e a te te r s ” )
−−>
S e r ie s o f S t ir r ed
3 4 Cp=2000,A=1,U=200,m=1000,mdot=2,Ts=250
/ / g i v e n d a ta ta 5 T 0 = 20 20 , T 1 =0 =0 , T 2 =0 =0 , T 3 =0 =0 6 7 / / fr fr o m e n nee rg r g y b a l a n c es e s f o r t he h e t a nk n k s we w e h av av e
a c c u m u l a t i o n = i n l e t −o u t l e t 8 T1_steady=(mdot*Cp*(T0)+U*A*(Ts))/(mdot*Cp+U*A) 9 disp (T1_steady , ” t he h e s t e aad d y s t a t e t em e m pe p e ra r a tu t u re r e o f t a nk nk 1 i s ”); 10 T2_steady=(mdot*Cp*(T1_steady)+U*A*(Ts))/(mdot*Cp+U* A) 11 disp (T2_steady , ” t he h e s t e aad d y s t a t e t em e m pe p e ra r a tu t u re r e o f t a nk nk 2 i s ”); 12 T3_steady=(mdot*Cp*(T2_steady)+U*A*(Ts))/(mdot*Cp+U* A) 13 disp (T3_steady , ” t he h e s t e aad d y s t a t e t em e m pe p e ra r a tu t u re r e o f t a nk nk 3 i s ”); 14 final_T3=.99*T3_steady 15 function dT1_by_dt=f1(t,T1,T2,T3), 16 dT1_by_dt=(mdot*Cp*(T0-T1)+U*A*(Ts-T1))/(m*Cp),
46
17 18 19 20 21 22 23 24 25 26 27
endfunction function dT2_by_dt=f2(t,T1,T2,T3), dT2_by_dt=(mdot*Cp*(T1-T2)+U*A*(Ts-T2))/(m*Cp), endfunction function dT3_by_dt=f3(t,T1,T2,T3), dT3_by_dt=(mdot*Cp*(T2-T3)+U*A*(Ts-T3))/(m*Cp), endfunction T1=20,T2=20,T3=20
/ / s o l v i n g b y N e w t o n ’ s M e t ho ho d f o r t=0:1:10000, h =1
// st ep
i n c re r e me m e nt nt o f 1 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
k1=h*f1(t,T1,T2,T3) l1=h*f2(t,T1,T2,T3) m1=h*f3(t,T1,T2,T3) k2=h*f1(t+h/2,T1+k1/2,T2+l1/2,T3+m1/2) l2=h*f2(t+h/2,T1+k1/2,T2+l1/2,T3+m1/2) m2=h*f3(t+h/2,T1+k1/2,T2+l1/2,T3+m1/2) k3=h*f1(t+h/2,T1+k2/2,T2+l2/2,T3+m2/2) l3=h*f2(t+h/2,T1+k2/2,T2+l2/2,T3+m2/2) m3=h*f3(t+h/2,T1+k2/2,T2+l2/2,T3+m2/2) k4=h*f1(t+h,T1+k3,T2+l3,T3+m3) l4=h*f2(t+h,T1+k3,T2+l3,T3+m3) m4=h*f3(t+h,T1+k3,T2+l3,T3+m3) T1=T1+(k1+2*k2+2*k3+k4)/6 T2=T2+(l1+2*l2+2*l3+l4)/6 e1 = a b s (T3-final_T3) p p ro r o x . t i me m e wh e n if e1<1e-3 t he h e n d i sp sp ( t , ” t h e a pp
T e mp m p er e r at a t ur u r e i n 3 r d t an a n k i s 9 9% o f s t ea ea d y v a lu l u e i s ” ) ; break 44 end 45 T3=T3+(m1+2*m2+2*m3+m4)/6 46 e n d
Scilab code Exa 4.11 batch and stirred tank reactor
47
1 clc 2 // / / b a tc tc h r e a c t o r s 3 disp ( ” t he he s o l u t i o n
o f e . g . 4 . 11
−−>
B at a t c h a nd nd S t i r r e d
Ta n k R e a c t o r s ” ) 4 / / r xn x n g i ve v e n A−−> B 5 6 7 8 9 10 11
rate_const_k=1 function dCa_by_dt=fs1(t,Ca), dCa_by_dt=-rate_const_k*Ca, endfunction Ca=1 f o r t=0.1:0.1:3, h=0.1
// st ep
i n c re r e me m e nt nt o f 0 . 1 12 k1=h*fs1(t,Ca) 13 k2=h*fs1(t+h/2,Ca+k1/2) 14 k3=h*fs1(t+h/2,Ca+k2/2) 15 k4=h*fs1(t+h,Ca+k3) 16 Ca=Ca+(k1+2*k2+2*k3+k4)/6 17 e n d 18 disp (Ca, ” t h e v a lu l u e o f c on o n c . a t t =3 =3 u s i ng n g Runge K u t t a m et et ho ho d i s ” ) ; 19 Ca_anl= e x p (-t) // a n a l y t i c a l
solution 20 disp (Ca_anl, ” t he he a n a l y t i c a l
s o ln . i s ”)
Scilab code Exa 4.12 batch and stirred tank reactor
1 2 3 4
clc
/ / rx r x n A−−>B // i np u t= t=FC FCa0 a0 , ou tp ut= ut=F FCa / / a p p l y i n g m a s s b a l a n c e o f c o mp m p o ne n e n t A we we g e t d (V ( V∗ Ca )/dt=F ∗ Ca0 −F ∗ Ca−k ∗ Ca ∗V 5 disp ( ” t he h e s o l u t i o n o f e . g . 4 . 1 2 −−>B at a t c h a nd nd S t i r r e d Ta n k R e a c t o r s ” ) 6 rate_const_k=1
48
7 8 9 10 11 12 13 14
C a 0 = 1 , F =1 =1 , V = 1 0 function dVCa_by_dt=fr(t,Ca1), dVCa_by_dt=F*Ca0-F*Ca1-rate_const_k*Ca1*V, endfunction Ca1=1 f o r t=0.1:0.1:10, // st ep h=0.1
i n c re r e me m e nt nt o f 0 . 1 15 k1=h*fr(t,Ca1) 16 k2=h*fr(t+h/2,Ca1+k1/2) 17 k3=h*fr(t+h/2,Ca1+k2/2) 18 k4=h*fr(t+h,Ca1+k3) 19 Ca1=Ca1+(k1+2*k2+2*k3+k4)/6 20 e n d // f i n a l v al ue 21 disp (Ca1, ” t h e v a lu l u e o f Ca a t t =1 =10 s u s i ng n g Rung e K u tt t t a m et et ho ho d i s ” ) ; 22 Ca_steady=F*Ca0/(F+rate_const_k*V) 23 disp (Ca_steady , ” t he h e s te t e ad a d y s t a t e v al a l ue ue o f conc . i s ” ) ;
Scilab code Exa 4.13 batch and stirred tank reactor
1 2 3 4 5
clc
/ / g i v e n r xn x n A−−>B−−>C / / g i v e n d a ta ta k 1 = 1 , k 2 =1 =1 a t ch ch R e a c t o r s ” ) disp ( ” t h e s o l u t i o n o f e g 4 . 1 3 −−>B at // function dA_by_dt=f1a(t,A,B,C), f u n c t i o n s d e f in e d
6 7 8 9 10 11
dA_by_dt=-A, endfunction function dB_by_dt=f2a(t,A,B,C), dB_by_dt=A-B, endfunction function dC_by_dt=f3a(t,A,B,C),
49
12 dC_by_dt=B, 13 endfunction 14 A=1,B=0,C=0 15 f o r t=0.1:.1:10, 16 h=.1
// i n i t i a l va lu es // st ep
i n c re r e me m e nt nt o f 0 . 1 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
k1=h*f1a(t,A,B,C) l1=h*f2a(t,A,B,C) m1=h*f3a(t,A,B,C) k2=h*f1a(t+h/2,A+k1/2,B+l1/2,C+m1/2) l2=h*f2a(t+h/2,A+k1/2,B+l1/2,C+m1/2) m2=h*f3a(t+h/2,A+k1/2,B+l1/2,C+m1/2) k3=h*f1a(t+h/2,A+k2/2,B+l2/2,C+m2/2) l3=h*f2a(t+h/2,A+k2/2,B+l2/2,C+m2/2) m3=h*f3a(t+h/2,A+k2/2,B+l2/2,C+m2/2) k4=h*f1a(t+h,A+k3,B+l3,C+m3) l4=h*f2a(t+h,A+k3,B+l3,C+m3) m4=h*f3a(t+h,A+k3,B+l3,C+m3) A=A+(k1+2*k2+2*k3+k4)/6 B=B+(l1+2*l2+2*l3+l4)/6 C=C+(m1+2*m2+2*m3+m4)/6 if t = = . 5 | t = = 1| 1 | t = = 2 | t = = 5 t h en en d i sp sp ( C , B , A , ” s e c s i s ” , t , ” t h e c on o n c . o f A , B ,C ,C a f t e r ” ); end
33 34 e n d 35 disp (C,B,A, ” t he he c o n c . r e s p e c t i v e l y i s ” );
o f A , B ,C ,C a f t e r 10 s e c s
Scilab code Exa 4.14 batch and stirred tank reactor
1 2 3 4
clc
/ / g i v e n r x n A+B−−k1 −−>C // B+C−−k2 −−>D //given rate
k 1 = 1 , k 2 =1 =1
constants 50
5 6 7 8 9 10 11 12 13 14 15 16 17 18
a t ch ch R e a c t o r s ” ) disp ( ” t h e s o l u t i o n o f e g 4 . 1 4 −−>B at function dA_by_dt=f1a(t,A,B,C,D), dA_by_dt=-A*B, endfunction function dB_by_dt=f2a(t,A,B,C,D), dB_by_dt=-A*B-B*C, endfunction function dC_by_dt=f3a(t,A,B,C,D), dC_by_dt=A*B-B*C, endfunction function dD_by_dt=f4a(t,A,B,C,D), dD_by_dt=B*C, endfunction // i n i t i a l A=1,B=1,C=0,D=0 values
19 f o r t=.1:.1:10, 20 h=.1
// st ep
i n c re r e me m e nt nt o f 0 . 1 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
k1=h*f1a(t,A,B,C,D) l1=h*f2a(t,A,B,C,D) m1=h*f3a(t,A,B,C,D) n1=h*f4a(t,A,B,C,D) k2=h*f1a(t+h/2,A+k1/2,B+l1/2,C+m1/2,D+n1/2) l2=h*f2a(t+h/2,A+k1/2,B+l1/2,C+m1/2,D+n1/2) m2=h*f3a(t+h/2,A+k1/2,B+l1/2,C+m1/2,D+n1/2) n2=h*f4a(t+h/2,A+k1/2,B+l1/2,C+m1/2,D+n1/2) k3=h*f1a(t+h/2,A+k2/2,B+l2/2,C+m2/2,D+n2/2) l3=h*f2a(t+h/2,A+k2/2,B+l2/2,C+m2/2,D+n2/2) m3=h*f3a(t+h/2,A+k2/2,B+l2/2,C+m2/2,D+n2/2) n3=h*f4a(t+h/2,A+k2/2,B+l2/2,C+m2/2,D+n2/2) k4=h*f1a(t+h,A+k3,B+l3,C+m3,D+n3) l4=h*f2a(t+h,A+k3,B+l3,C+m3,D+n3) m4=h*f3a(t+h,A+k3,B+l3,C+m3,D+n3) n4=h*f4a(t+h,A+k3,B+l3,C+m3,D+n3) A=A+(k1+2*k2+2*k3+k4)/6 B=B+(l1+2*l2+2*l3+l4)/6 C=C+(m1+2*m2+2*m3+m4)/6 D=D+(n1+2*n2+2*n3+n4)/6
51
if t = = . 5 | t = = 1| 1 | t = = 2 | t = = 5 t h en en d i sp sp (D,C,B,A, ” s e c s i s ” ,t , ” t h e c on o n c . o f A , B ,C , C ,D ,D a f t e r ” ) ; end
41
42 43 e n d 44 h e c oon nc. disp (D,C,B,A, ” t he r e s p e c t i v e l y i s ” );
o f A ,B ,C, D a f t e r 10 s e c s
Scilab code Exa 4.15 plug flow reactor
1 2 3 4 5 6 7 8 9 10
clc l u g F l o w R e a ct ct o r ” ) disp ( ” t h e s o l u t i o n o f e g 4 . 1 5 −−>P lu / / g i ve v e n r a t e c o n st st a n t rc_k1=1 / / me me a n a x i a l v e l o c i t y u =1 function dCa_by_dx=fm(x,Ca), dCa_by_dx=-Ca, endfunction Ca=1 f o r x=.1:.1:10, // st ep h=0.1
i n c re r e me m e nt nt o f 0 . 1 11 12 13 14 15 16
k1=h*fm(x,Ca) k2=h*fm(x+h/2,Ca+k1/2) k3=h*fm(x+h/2,Ca+k2/2) k4=h*fm(x+h,Ca+k3) Ca=Ca+(k1+2*k2+2*k3+k4)/6 if x==.5|x==1|x==2|x==5 t he h e n d i sp sp (Ca, ” l e n g t h h e v a l u e o f c o nc nc . a t ” ) ; ,x , ” t he end
is ”
17 18 e n d 19 disp (Ca, ” t h e v a l u e o f Ca a t x =1 0 u s i n g R u n g e K ut ut t a m e t h o d i n p lu lu g f lo w r e a c t o r i s ” );
Scilab code Exa 4.16 plug flow reactor
52
1 clc 2 / / g i v e n r xn x n A−−>B−−>C 3 r c _ k 1 = 1 , r c _k _k 2 = 1
// gi ven
r a t e c o n s t a nt nt s / /m / m ea n a x i a l
4 u =1
velocity 5 disp ( ” t h e s o l u t i o n o f e g 4 . 1 6 6 7 8 9 10 11 12 13 14 15 16 17
−−>
function dA_by_dx=f1e(x,A,B,C), dA_by_dx=-A, endfunction function dB_by_dx=f2e(x,A,B,C), dB_by_dx=A-B, endfunction function dC_by_dx=f3e(x,A,B,C), dC_by_dx=B, endfunction A=1,B=0,C=0 f o r x=.1:.1:10, h=.1
P lu l u g F l o w R e a ct ct o r ” )
// st ep
i n c re r e me m e nt nt o f 0 . 1 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
k1=h*f1e(x,A,B,C) l1=h*f2e(x,A,B,C) m1=h*f3e(x,A,B,C) k2=h*f1e(x+h/2,A+k1/2,B+l1/2,C+m1/2) l2=h*f2e(x+h/2,A+k1/2,B+l1/2,C+m1/2) m2=h*f3e(x+h/2,A+k1/2,B+l1/2,C+m1/2) k3=h*f1e(x+h/2,A+k2/2,B+l2/2,C+m2/2) l3=h*f2e(x+h/2,A+k2/2,B+l2/2,C+m2/2) m3=h*f3e(x+h/2,A+k2/2,B+l2/2,C+m2/2) k4=h*f1e(x+h,A+k3,B+l3,C+m3) l4=h*f2e(x+h,A+k3,B+l3,C+m3) m4=h*f3e(x+h,A+k3,B+l3,C+m3) A=A+(k1+2*k2+2*k3+k4)/6 B=B+(l1+2*l2+2*l3+l4)/6 C=C+(m1+2*m2+2*m3+m4)/6 mt r i s if x = = . 5 | x = = 1| 1 | x = = 2 | x = = 5 t h en en d i sp sp ( C , B , A , ” mt ” , x , ” t he h e c oon n c . o f A , B , C a t a d i s t a n ce ce o f ” ); end
53
35 e n d 36 he c o n c . disp ( C , B , A , ” t he 1 0 mtr i s ” ) ;
o f A , B ,C , C a t a d i s t a n ce ce o f
Scilab code Exa 4.17 plug flow reactor
1 2 3 4
clc
/ / g i v e n r x n A+B−−k1 −−>C // B+C−−k2 −−>D // ra te
rc_k1=1,rc_k2=1
constants 5 disp ( ” t h e s o l u t i o n o f e g 4 . 1 7 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
−−>
function dA_by_dx=f1a(x,A,B,C,D), dA_by_dx=-A*B, endfunction function dB_by_dx=f2a(x,A,B,C,D), dB_by_dx=-A*B-B*C, endfunction function dC_by_dx=f3a(x,A,B,C,D), dC_by_dx=A*B-B*C, endfunction function dD_by_dx=f4a(x,A,B,C,D), dD_by_dx=B*C, endfunction A=1,B=1,C=0,D=0 f o r x=.1:.1:10, h=.1
P lu l u g F l o w R e a ct ct o r ” )
// st ep
i n c re r e me m e nt nt o f 0 . 1 21 22 23 24 25 26 27
k1=h*f1a(x,A,B,C,D) l1=h*f2a(x,A,B,C,D) m1=h*f3a(x,A,B,C,D) n1=h*f4a(x,A,B,C,D) k2=h*f1a(x+h/2,A+k1/2,B+l1/2,C+m1/2,D+n1/2) l2=h*f2a(x+h/2,A+k1/2,B+l1/2,C+m1/2,D+n1/2) m2=h*f3a(x+h/2,A+k1/2,B+l1/2,C+m1/2,D+n1/2)
54
28 29 30 31 32 33 34 35 36 37 38 39 40 41
n2=h*f4a(x+h/2,A+k1/2,B+l1/2,C+m1/2,D+n1/2) k3=h*f1a(x+h/2,A+k2/2,B+l2/2,C+m2/2,D+n2/2) l3=h*f2a(x+h/2,A+k2/2,B+l2/2,C+m2/2,D+n2/2) m3=h*f3a(x+h/2,A+k2/2,B+l2/2,C+m2/2,D+n2/2) n3=h*f4a(x+h/2,A+k2/2,B+l2/2,C+m2/2,D+n2/2) k4=h*f1a(x+h,A+k3,B+l3,C+m3,D+n3) l4=h*f2a(x+h,A+k3,B+l3,C+m3,D+n3) m4=h*f3a(x+h,A+k3,B+l3,C+m3,D+n3) n4=h*f4a(x+h,A+k3,B+l3,C+m3,D+n3) A=A+(k1+2*k2+2*k3+k4)/6 B=B+(l1+2*l2+2*l3+l4)/6 C=C+(m1+2*m2+2*m3+m4)/6 D=D+(n1+2*n2+2*n3+n4)/6 if x = = . 5 | x = = 1| 1 | x = = 2 | x = = 5 t h en en d i sp sp (D,C,B,A, ” s e c s i s ” ,x , ” t h e c on o n c . o f A , B ,C , C ,D ,D a f t e r ” ) ; end
42 43 e n d 44 h e c oon nc. disp (D,C,B,A, ” t he r e s p e c t i v e l y i s ” );
o f A ,B ,C, D a f t e r 10 s e c s
Scilab code Exa 4.18 non isothermal plug flow reactor
1 clc 2 disp ( ” t he h e s o l u t i o n o f e g 4. 4 .18 −− >Non− I s o t h e rm r m a l P lu lu g F lo lo w R e a c t o r ” ) 3 T=294.15 4 / / rx r x n A−−>B 5 R = 8. 8 . 31 31 4 , r ho ho = 9 80 80 .9 . 9 , M W =2 = 2 00 00 , U = 19 1 9 00 00 , C p =1 = 1 5. 5. 7 , H _r _ r xn xn = 92 9 2 90 9 0 0 , T 1 = 38 38 8. 8 . 71 7 1 , m do d o t = 1. 1 . 26 26 , d ia ia = 2 .5 . 5 4* 4 * 10 10 ^ - 2 , E / / g i v e n d a ta ta =108847 6 b=E/(R*T1),k1=5.64*10^13* e x p ( - b) b ) , A = %p % p i * di di a ^2 ^2 /4 /4 , n a0 a0 = m d o t * 1 0 0 0/ 0/ M W , T s = 3 8 8. 8. 7 1 7 k = k 1 * e x p (b*(1-T1/T)) 8 9 //dX by dV=ra/na0
55
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
//dX by dV=k ∗ (1 −X)/F / / f ro r o m e ne n e rg rg X b a l a n c e //mdot ∗ Cp ∗ d T b y d z + r a ∗A∗ H RXN−q=0 //q=U∗ %p i ∗ d i a ∗ ( Ts−T ) / /−mdot ∗ Cp ∗ dT by dV+4 ∗U / d i a ∗ ( Ts−T )− r a ∗ H rxn=0 F=mdot/rho t1=A*k1/F s1=mdot*Cp/A s2=4*U/dia s3=H_rxn*t1 function dX_by_dz=fg1(z,X,T), dX_by_dz=t1*(1-X)* e x p (b*(1-T1/T)) endfunction function dT_by_dz=fd1(z,X,T), ra=na0/A*(t1*(1-X)* e x p (b*(1-T1/T))) dT_by_dz=(ra*H_rxn-s2*(Ts-T))/-s1 endfunction X=0,T=294.15 f o r z=0:.1:350, h=.1
// sz ep
i nc n c re r e me m e nz nz o f 0 . 1 34 35 36 37 38 39 40 41 42 43 44 45
k1=h*fg1(z,X,T) l1=h*fd1(z,X,T) k2=h*fg1(z+h/2,X+k1/2,T+l1/2) l2=h*fd1(z+h/2,X+k1/2,T+l1/2) k3=h*fg1(z+h/2,X+k2/2,T+l2/2) l3=h*fd1(z+h/2,X+k2/2,T+l2/2) k4=h*fg1(z+h,X+k3,T+l3) l4=h*fd1(z+h,X+k3,T+l3) X=X+(k1+2*k2+2*k3+k4)/6 T=T+(l1+2*l2+2*l3+l4)/6
/ / c o n d i t i o n f o r h e ig i g h t c a l c . f o r 9 0% c o n ve v e r s io io n he h e i g h t o f t h e if X > . 9 & X < . 9 00 0 0 5 t he h e n d i sp sp ( z , ” t he t ow o w er e r r e q u i r e d f o r 9 0% c o n ve v e r s io io n i n mtrs i s ” 56
46 47
) ; break end end
57
Chapter 5 BOUNDARY VALUE PROBLEMS
Scilab code Exa 5.1 discretization in 1 D space
1 clc 2 / / f i n i t e d i f f e r e n c e m et et ho ho d 3 disp ( ” t he h e s o l u t i o n o f e g 55..1−− > D i s c r e t i z a t i o n i n 1−D s p a c e ”); 4 / / g i v e n d 2 y b y d x2 x 2 − 2=0 h e n ce ce i t i s d i r i c h l e t ’ s
problem 5 6 x _ 1 = 0 , y _1 _1 = 0
// i n i t i a l boun dary
conditions 7 x _ 3 = 1 , y _3 _3 = 0 8 9 delta_x=.5
/ / s i n c e we h a ve ve t o f i n d
y 2 a t x =. =. 5 s o x 2 = . 5 10 / / i n t he he c e n t r a l d i f f e r e n c e meth od s u b s t i t u t i n g we h a v e 11 / / ( y 3 −2∗ y 2 +y +y 1 ) / ( d e l t a x ) ˆ 2= 2= 2 12 / / s i n c e y 1 & y 3 =0 =0 a s p er e r B .C .C . 13 y_2=(y_3+y_1-2*delta_x^2)/2 14 disp (y_2, ” t h e v a l u e o f y a t x =. = . 5 f ro ro m f i n i t e
58
i =2 =2
d i f f e r e n c e meth od i s ” ) ; 15 x=.5 16 y_exact=x^2-x 17 disp (y_e xact , ” t he h e v al a l ue u e o f y f r o m e xa x a ct c t s o l n a t x =. =. 5 i s ” );
Scilab code Exa 5.2 discretization in 1 D space
1 clc 2 disp ( ” t h e s o l u t i o n o f e g 5 .2 . 2 −−> D i s c r e t i z a t i o n i n 1− D s p a c e ”); 3 / / bo b o un u n da d a ry r y c o n d i t i o n s a r e : x =0 a t y =0 = 0 ; d y / dx dx =1 a t x
=1 4 disp ( ” t o s o l v e
t h i s p r o b l e m we w e w i l l t ak a k e d e l t a x =. =. 5 s i n c e we h av a v e t o f i n d t he h e v a l u e a t x =. =. 5 ” ) ;
5 6 7 8 9
delta_x=.5 y_1=0
// u si ng c e n t r a l d i f f e r e n c e eqn / / a t x =1 =1 , i =3 =3 dy_by_dx=1 // y 4= 4=dy/dx dy/dx ∗ 2 ∗ d e l t a x + y 2 sin cefro m B.C. a t node 3
10 11 / / y 2 = d e l t a x ˆ 2+ y 3 − d e l t a x
on
s u b s t i t u t i n g t h e v al a l ue ue o f y 4 12 y_3=-(2*delta_x^2+2*(delta_x^2-delta_x)-y_1) / / o n
s u b s t i t u t i n g t h e v al a l ue ue o f y 2 13 y_2=delta_x^2+y_3-delta_x 14 disp (y_2, ” t he h e v al a l u e o f y a t x =. =. 5 i s ” ) ;
Scilab code Exa 5.3 discretization in 1 D space
1 clc
59
2 disp ( ” t h e s o l u t i o n o f e g 5 .3 . 3 −−> D i s c r e t i z a t i o n i n 1− D s p a c e ”); 3 / / g i v e n t h e s o u r c e t er e r m f ( x ) =4 =4 x ˆ2 ˆ2 − 2 x −4 4 / / g i v e n e qn q n d 2y 2y / d x2 x2 −2y 2y=f =f ( x ) 5 y_1=0 6 y_4=-1 7 delta_x=1/3 / / s i n c e g i v en e n 3 p a r ts ts and
length=1 8 9 10 11 12 13 14 15 16 17 18
f o r i=0:3,j=0:delta_x:1; end
/ / g iv i v en en t o d i v i d e t h e l i n e i n 3 p a r t s / / a t n od od e 2 / / y 3 − 2∗ y 2
function d = f x 3 ( x ) , d=(4*x^2-2*x-4) endfunction f2=fx3(j(2)) f3=fx3(j(3)) y_3=((f2)*delta_x^2+(2+2*delta_x^2)*((f3)*delta_x^2y_4)-y_1)/(1-(2+2*delta_x^2)^2) 19 y_2=(f3+2*y_3)*delta_x^2+2*y_3-y_4 20 disp (y_3,y_2, ” i s r e s p e c t i v e l y ” ,j(3),j(2), ” t h e v a lu lu e o f y a t x=” ) ;
Scilab code Exa 5.4 discretization in 1 D space
1 2 3 4 5 6
// f ( x )= )=4x2 4x2 −2x −4
clc disp ( ” t h e s o l u t i o n y_1=0 dy_by_dx=-3 delta_x=1/3
o f e x 5 . 4 b y TD TDMA me me t h o d i s ” ) ; / / a t x =1 / / s i n c e g i v en e n 3 p a r ts ts and
length=1 7 f o r i=0:3,j=0:delta_x:1; 8 end
60
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
/ / g iv i v en en t o d i v i d e t h e l i n e i n 3 p a r t s
function d = f x 3 ( x ) , d=(4*x^2-2*x-4) endfunction f2=fx3(j(2)) f3=fx3(j(3)) f4=fx3(j(4)) a c t a n a l y t i c a l s o l n a re re ” ); disp ( ” t h e e x ac f o r i=1:3,y(i)=-2*j(i+1)^2+j(i+1), disp ( y ( i ) ) ; end
/ / u s i n g B . C. C. − 2 a t n o d e 4 we w e g et et / / ( y 5 − y 3 ) /(2 ∗ d e l t a x ) =−3 / / eq eq n a t n o d e 2 / / − 20 ∗ y 2 + 9∗ y 3 = f 2 / / a t n od od e 3 / / 9 ∗ y 2 − 20 ∗ y 3 + 9∗ y 4 = f 3 / / a t n od od e 4 / / 1 8 ∗ y 3 − 20 ∗ y 4 = 1 6 h e v al a l u e o f f ( x ) a t n o d e 2 3 and 4 disp (f4,f3,f2, ” t he a r e ”); / / no no w s o l v i n g t h e e q u a t i o n s u s i n g TD TDMA m et et ho ho d
28 29 30 a(2)=9,a(3)=18 31 32 f o r j = 1: 1: 3 , b ( j ) = - 2 0; 0;
/ / s ub u b d i a g o n a l a s si s i g nm n m e nt nt / / ma ma i n d i a g o n a l
assignment 33 e n d 34 f o r k = 1: 1: 2 , c ( k ) = 9 ;
/ / s u p er er d i a g o n a l
assignment 35 e n d 36 d(1)=f2
//given values
assignment 37 38 39 40 41 42
d(2)=f3 d(3)=16 i=1; n=3; beta1(i)=b(i);
// i n i t i a l b is equal 61
t o b et e t a s i n c e a 1 =0 43 44 45 46 47 48 49 50 51 52 53 54
/ / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); end TDMA me method a r e ” ) ; disp ( ” t h e v a l u e s o f y 2 , y 3 , y 4 b y TD f o r i=1:3, disp (x(i)); end
Scilab code Exa 5.5 discretization in 1 D space
1 2 3 4 5 6 7 8
clc h e s o l n o f e q n 5. 5 .5−− > ” ) ; disp ( ” t he delta_x=.1 y_11=1 / / d y / dx dx a t x =0 dy_by_dx_1=0
/ / g i v e n e qn q n d 2y 2 y / d x2 x 2=y re ”); disp ( ” t h e a n a l y t i c a l s o l n a re
f o r i = 1: 1: 1 0 , y _ an an ( i ) = cosh ((i-1)/10)/ cosh (1), disp ( y_an(i)); 9 end 10 / / u s i ng n g c e n t r a l d i f f e r e n c e m e t h o d we h a ve ve 11 / / y ( i − 1) − ( 2+ 2+ d e l a t x ˆ 2 ) y ( i ) + y ( i + 1) 1) =0 =0 12 / / t h e r e f o r e t h e e q ns n s c an a n b e s o l v e d u s i n g TD TDMA m e t h o d 13 f o r i = 2: / / s ub ub d i a g o n a l 2: 1 0 , a ( i ) = 1
assignment 14 e n d 15 f o r j = 1 :1 :1 0 , b ( j ) = - 2 . 01 01 ;
/ / ma ma i n d i a g o n a l
assignment 62
16 e n d 17 c(1)=2 18 f o r k = 2: 2: 9 , c ( k ) = 1 ;
/ / s u p er er d i a g o n a l
assignment 19 e n d 20 f o r l = 1: 1: 9 , d ( l ) = 0 ; 21 end
//given values
assignment 22 23 24 25
d(10)=-1 i=1; n=10; beta1(i)=b(i);
// i n i t i a l b is equal
t o b et e t a s i n c e a 1 =0 26 27 28 29 30 31 32 33 34 35 36 37
/ / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); end r o m y 1 t o y 10 1 0 b y TDMA me me t h o d disp ( ” t h e v a l u e s o f y f ro a r e ”); f o r i=1:10, disp ( x ( i ) ) ; end
Scilab code Exa 5.6 1 D steady state heat conduction
1 clc 2 disp ( ” t he h e s o l n o f e q n 5. 5 .6−− >1−D S t ea e a d y H ea ea t Conduction” ); 3 / / i n t h e g i v e n p ro r o bl b l em em 4 T _ 1 = 1 00 00 , T _ 10 10 = 2 0 0
63
// s i n c e 9
5 delta_x=(T_10-T_1)/9
d i v i s i o n s a r e t o b e ma made 6 / / u s in i n g c e n t r a l d i f f e r e n c e m e t h o d we g e t 7 / / f o r n od o d e 2−−> 2 ∗ T 2 −T 3=100 8 f o r i = 2 :8 / / s ub ub d i a g o n a l :8 , a ( i ) = -1 -1 assignment 9 end 10 f o r j = 1: 1: 8 , b ( j ) = 2 ;
/ / ma ma i n d i a g o n a l
assignment 11 e n d 12 f o r k = 1: 1: 7 , c ( k ) = - 1 ;
/ / s u p er er d i a g o n a l
assignment 13 14 15 16 17 18 19
end d ( 1 ) = 1 00 00 , d ( 8 ) = 2 00 00 f o r l = 2: 2: 7 , d ( l ) = 0 ; end i=1; n=8; beta1(i)=b(i);
/ / g i v e n v a l u e s a s s ig i g n m en en t
// i n i t i a l b is equal
t o b et e t a s i n c e a 1 =0 20 21 22 23 24 25 26 27 28 29 30 31
/ / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); end r o m T2 t o T9 b y TD TDMA me m e th th od od disp ( ” t h e v a l u e s o f T f ro a r e ”); f o r i=1:8, disp (x(i)); end
64
Scilab code Exa 5.7 1 D steady state heat conduction
1 clc 2 disp ( ” t he h e s o l n o f e q n 5. 5 .7−− >1−D S t ea e a d y H ea ea t Conduction” ); 3 dia=.02 4 l=.05 5 T_0=320 6 delta_x=l/4 7 k=50 8 h=100 9 T_surr=20 10 / / B . C−−> d ( t h e t a ) / d x+ x+h ( t h e t a ) / k= k =0 a t x = 0 . 0 5 11 // l e t m=sq rt (hP (hP/k /kA) A) 12 P=%pi*dia 13 A=%pi*dia^2/4 14 m = sqrt (h*P/(k*A)); 15 / / u s in i n g c e n t r a l d i f f e r e n c e m e t h o d we g e t 16 / /− t h e t a ( i − 1)+(2+(m∗ d e l t a x ˆ 2 ) ∗ th et a ( i )+th )+th et a ( i +1 +1)) )
=0 17 theta_0=T_0-T_surr 18 / / u s i ng n g B . C. C . a t n o d e 4 we g et e t −−> t h e t a ( 5 ) =t =t h e t a ( 3 ) −2
h ∗ t h e t a ( 4 ) ∗ d e l t a x / k 19 / / no n o w t h e e q n s c an a n b e s o l v e d u s i n g TDMA me me t ho ho d 20 f o r i = 2 :3 / / s ub ub d i a g o n a l :3 , a ( i ) = -1 -1 assignment 21 e n d 22 a(4)=-2 23 f o r j = 1 :3 :3 , b ( j ) = 2 . 0 62 62 5 ;
/ / ma ma i n d i a g o n a l
assignment 24 e n d 25 b(4)=2.1125 26 f o r k = 1: 1: 3 , c ( k ) = - 1 ;
/ / s u p er er d i a g o n a l
assignment 27 e n d 28 f o r l = 2: 2: 4 , d ( l ) = 0 ; 29 end
//given values
assignment 65
30 31 32 33
d(1)=300 i=1; n=4; beta1(i)=b(i);
// i n i t i a l b is equal
t o b et e t a s i n c e a 1 =0 34 35 36 37 38 39 40 41 42 43 44 45
/ / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); end h e v a l u es e s o f T f ro r o m T1 t o T4 i n C e l s i u s b y disp ( ” t he TDMA method a r e ” ) ; f o r i=1:4, disp (x(i)-T_surr); end
Scilab code Exa 5.8 chemical reaction and diffusion in pore
1 clc 2 disp ( ” t he he s o l n
o f e g 5. 5 .8−− > C he h e mi m i ca ca l R e a c t i o n and D i f f u s i o n i n Pore ”);
3 4 5 6 7 8 9
lnght=.001 k_const=.001 D=10^-9 delta_x=l/100 C =1
// in mol/m mol/m33 / /B / B . C . a r e C=1 a t x =0 // dC / d x =0 a t x =1 0ˆ −3 s i n c e a t t h e end p o in i n t c o nc nc . i s c o n s t . 10 / / u s in i n g c e n t r a l d i f f e r e n c e m e t h o d we g e t t he he f o l l o w i n g e qn q n s w hi h i ch c h c an a n b e s o l v e d u s i n g TD TDMA 66
method / / s ub ub d i a g o n a l
11 f o r i = 2: 2: 9 9 , a ( i ) = 1
assignment 12 e n d 13 a(100)=2
/ / s i n c e C 9 9=C 1 00 00 u s i n g B . C
. 14 f o r j = 1 : 10 10 0 , b ( j ) = - 2 . 00 00 01 01 ,
/ / ma m a in in d i a g o n a l
assignment 15 e n d 16 f o r k = 1: 1 : 9 9 , c ( k ) = 1; 1;
/ / s u p er er d i a g o n a l
assignment 17 18 19 20 21 22 23
end d(1)=-1 f o r l = 2 :1 :1 00 0 0 , d ( l ) = 0; 0; end i=1; n=100; beta1(i)=b(i);
/ / g i v e n v a l u e s a s s ig i g n m en en t
// i n i t i a l b is equal
t o b et e t a s i n c e a 1 =0 / / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); 32 e n d 33 disp ( x ( 5 0 ) , ” t h e v a l u e s o f c on o n c . a t x =. = . 5mm or o r a t t he he 5 0 t h n o de de i s ” ) ; 24 25 26 27 28 29 30 31
67
Chapter 6 CONVECTION DIFFUSION PROBLEMS
Scilab code Exa 6.1 upwind scheme
1 / / g i ve v e n c o n v e c t iv iv e
d i f f u s i v e e q n −−>
u ∗ ( dc /dx )+ )+D D∗ (
−
d2C/dx2)=0 2 3 4 5 6 7
/ / a t x =0 C_ini=0 / / a t x =1 C_end=1 clc disp ( ” t h e s o l n o f e g 6 . 1 ” ) ;
/ / u s i n g c e n t r a l d i f f e r e n c e metho d f o r b o t h d i f f u s i o n a n d c o n v e c t i v e t e rm rm 8 / /−u ∗ (C( i +1)−C ( i − 1)) /(2 ∗ d e l t a x ) + D ∗ (C ( i +1 )+ )+C( C( i − 1) − 2 ∗C ( i ) ) / d e l t a x ˆ2 = 0 9 delta_x=1/50 10 / /o / o n s o l v i n g t he h e g i ve v e n e qn q n s a nd nd b y u s in i n g t he h e g i v en en
b o un u n d ar a r y e q n s we h a ve ve 11 Pe=50 12 Pe_local=50*delta_x 13 alpha= Pe_local -2
// gi ven / / u /D /D=50 a s l = 1 / / c o− e f f o f C( i
+1 ) 14 Beta=Pe_local+2
/ / c o− e f f o f C( i 68
1) / / m u l t i p l i n g w it i t h −2∗ d e l t a x ˆ 2/ 2 /D w e g e t / / −( P e l o c a l + 2 ) ∗C ( i − 1 ) + 4 ∗ C ( i ) + ( P e l o c a l − 2) ∗C ( i +1)=0 / / s o l v i n g e q n s u s i n g TDMA m et et ho ho d / / s ub ub d i a g o n a l f o r i = 2 :4 :4 9 , a ( i ) = - B et et a assignment −
15 16 17 18
19 e n d 20 f o r j = 1: 1: 4 9 , b ( j ) =4 =4 ,
/ / ma m a in in d i a g o n a l
assignment 21 e n d 22 f o r k = 1: 1 : 4 8 , c ( k ) = a l ph ph a ;
/ / s u p er er d i a g o n a l
assignment 23 24 25 26 27 28 29
end d ( 1 ) = B e ta ta * C _ in i n i , d ( 4 9 ) = - a lp l p h a * C _ e nd nd f o r l = 2: 2 : 4 8 , d ( l ) = 0; 0; / / g i v e n v a l u e s a s s ig i g n m en en t end i=1; n=49; // i n i t i a l b is equal beta1(i)=b(i);
t o b et e t a s i n c e a 1 =0 30 31 32 33 34 35 36 37 38 39
/ / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); end o n c . u s i n g CD CDS me m e t h o d f o r x =. = . 84 84 disp ( ” t h e v a l u e s o f c on t o 1 a r e ”) f o r i=42:49, disp ( x ( i ) ) ; end
40 41 42 / / p a r t ( i i ) u s i n g CDS a nd n d UDS m et et ho ho d 43 / / m u l t i p l i n g w it i t h − d e l t a x ˆ 2/ 2 / D we g e t 44 / / −( P e l o c a l + 1 ) ∗C ( i − 1 ) + ( P e l o c a l + 2 ) ∗C( i ) −C( i +1) =0
69
45 BEta=Pe_local+2 46 Gamma=Pe_local+1 47 f o r i = 2 :4 :4 9 , a ( i ) = - G a mm mm a
/ / s ub ub d i a g o n a l
assignment 48 e n d 49 f o r j = 1: 1 : 4 9 , b ( j ) = B Et Et a ,
/ / ma ma i n d i a g o n a l
assignment 50 e n d 51 f o r k = 1: 1 : 4 8 , c ( k ) = - 1; 1;
/ / s u p er er d i a g o n a l
assignment 52 53 54 55 56 57 58
end d ( 1 ) = G a mm m m a * C _ in in i , d ( 4 9 ) = C _ en en d f o r l = 2: 2 : 4 8 , d ( l ) = 0; 0; / / g i v e n v a l u e s a s s ig i g n m en en t end i=1; n=49; // i n i t i a l b is equal beta1(i)=b(i);
t o b et e t a s i n c e a 1 =0 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
/ / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); end n c . u s i n g CDS/UDS m et et h o d f o r x disp ( ” t h e v a l u e s o f c o nc = . 8 4 t o 1 a re re ” ) f o r i=42:49, disp ( x ( i ) ) ; end disp ( ” t h e a n a l y t i c a l s o l n i s ” ) ; f o r i=42:49; C_an(i)=C_ini+(( e x p (Pe*.02*i)-1)*(C_end-C_ini) / ( e x p (Pe)-1)); disp (C_an(i)); end 70
Scilab code Exa 6.2 upwind scheme
1 / / g i ve v e n c o n v e c t iv iv e
d i f f u s i v e e q n −−>
u ∗ ( dc /dx )+ )+D D∗ (
−
d2C/dx2)=0 2 3 4 5 6 7
/ / a t x =0 C_ini=0 / / a t x =1 C_end=1 clc h e s o l n o f e g 6. 6 .2−− > ” ) ; disp ( ” t he
/ / u s i n g c e n t r a l d i f f e r e n c e metho d f o r b o t h d i f f u s i o n a n d c o n v e c t i v e t e rm rm 8 / /−u ∗ (C( i +1)−C ( i − 1)) /(2 ∗ d e l t a x ) + D ∗ (C ( i +1 )+ )+C( C( i − 1) − 2 ∗C ( i ) ) / d e l t a x ˆ2 = 0 9 delta_x=1/50 10 / /o / o n s o l v i n g t he h e g i ve v e n e qn q n s a nd nd b y u s in i n g t he h e g i v en en
b o un u n d ar a r y e q n s we h a ve ve 11 Pe=500 12 Pe_local=500*delta_x 13 alpha= Pe_local -2
// giv en / / u /D /D=5 0 a s l =1 / / c o− e f f o f C( i
+1 ) 14 Beta=Pe_local+2
/ / c o− e f f o f C( i
1) / / m u l t i p l i n g w it i t h −2∗ d e l t a x ˆ 2/ 2 /D w e g e t / / −( P e l o c a l + 2 ) ∗C ( i − 1 ) + 4 ∗ C ( i ) + ( P e l o c a l − 2) ∗C ( i +1)=0 / / s o l v i n g e q n s u s i n g TDMA m et et ho ho d / / s ub ub d i a g o n a l f o r i = 2 :4 :4 9 , a ( i ) = - B et et a assignment −
15 16 17 18
19 e n d 20 f o r j = 1: 1: 4 9 , b ( j ) =4 =4 ,
/ / ma m a in in d i a g o n a l
assignment 21 e n d 22 f o r k = 1: 1 : 4 8 , c ( k ) = a l ph ph a ;
/ / s u p er er d i a g o n a l
assignment 71
23 24 25 26 27 28 29
end d ( 1 ) = B e ta ta * C _ in i n i , d ( 4 9 ) = - a lp l p h a * C _ e nd nd f o r l = 2: 2 : 4 8 , d ( l ) = 0; 0; / / g i v e n v a l u e s a s s ig i g n m en en t end i=1; n=49; // i n i t i a l b is equal beta1(i)=b(i);
t o b et e t a s i n c e a 1 =0 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
/ / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); end o n c . u s i n g CD CDS me m e t h o d f o r x =. = . 84 84 disp ( ” t h e v a l u e s o f c on t o 1 a r e ”) f o r i=42:49, disp ( x ( i ) ) ; end / / p a r t ( i i ) u s i n g CDS a nd n d UDS m et et ho ho d / / m u l t i p l i n g w it i t h − d e l t a x ˆ 2/ 2 / D we g e t / / −( P e l o c a l + 1 ) ∗C ( i − 1 ) + ( P e l o c a l + 2 ) ∗C( i )−C( i +1) =0 BEta=Pe_local+2 Gamma=Pe_local+1 f o r i = 2 :4 :4 9 , a ( i ) = - G a mm mm a
/ / s ub ub d i a g o n a l
assignment 48 e n d 49 f o r j = 1: 1 : 4 9 , b ( j ) = B Et Et a ,
/ / ma ma i n d i a g o n a l
assignment 50 e n d 51 f o r k = 1: 1 : 4 8 , c ( k ) = - 1; 1;
/ / s u p er er d i a g o n a l
assignment 52 e n d 53 d ( 1 ) = G a mm m m a * C _ in in i , d ( 4 9 ) = C _ en en d 54 f o r l = 2: 2 : 4 8 , d ( l ) = 0; 0;
72
55 56 57 58
end i=1; n=49; beta1(i)=b(i);
/ / g i v e n v a l u e s a s s ig i g n m en en t
// i n i t i a l b is equal
t o b et e t a s i n c e a 1 =0 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
/ / s i n c e c 7= 7 =0 gamma1(i)=d(i)/beta1(i); m=i+1; f o r j = m : n , b e ta ta 1 ( j ) = b ( j ) - a( a ( j ) * c ( j - 1) 1) / b e t a 1 ( j - 1) 1) ; gamma1(j)=(d(j)-a(j)*gamma1(j-1))/beta1(j); end / / s i n c e c 7= 7 =0 x(n)=gamma1(n); n1=n-i; f o r k = 1 : n1 n1 , j = n - k ; x ( j ) = g a mm m m a 1 ( j ) - c ( j )* ) * x ( j + 1) 1) / b e t a 1 ( j); end n c . u s i n g CDS/UDS m et et h o d f o r x disp ( ” t h e v a l u e s o f c o nc = . 8 4 t o 1 a re re ” ) f o r i=42:49, disp ( x ( i ) ) ; end disp ( ” t h e a n a l y t i c a l s o l n i s ” ) ; f o r i=42:49; C_an(i)=C_ini+(( e x p (Pe*.02*i)-1)*(C_end-C_ini) / ( e x p (Pe)-1)); disp (C_an(i)); end
73
Chapter 7 TUBULAR REACTOR WITH AXIAL DISPERSION
Scilab code Exa 7.1 boundary value problem in chemical reaction engi-
neering 1 clc 2 disp ( ” s o l n o f e g 7.1−− > F i r s t O r de d e r R ea e a c ti t i o n −5 0 p a r t s ”) 3 e1=1,e2=1
/ / a ss s s um u m ed ed v a l u e s 4 u=1,D=10^-4,k=1,C_a_in=1,delta_x=10/50
/ / g i v e n d a ta ta 5
cf_ca1_n1=-2*D/delta_x ^2-3*u/delta_x ^2-3*u/delta_x -k-2*u^2/D
/ / c o − e f f i c i e n t o f C−A1 a t n o d e 1 6 cf_ca2_n1=2*D/delta_x^2+u/delta_x 7 cf_da1_n1=-(2*u^2/D+2*u/delta_x)*C_a_in
/ / r i g h t hand s i d e co− efficient 8 cf_ca1_n2=D/delta_x^2+u/delta_x 9 cf_ca2_n2=-2*D/delta_x ^2-u/delta_x ^2-u/delta_x -k 10 cf_ca3_n2=D/delta_x^2 11 cf_da1_n2=0
74
12 13 14 15 16
cf_ca2_n3=cf_ca1_n2 cf_ca3_n3=cf_ca2_n2 cf_ca4_n3=cf_ca3_n2 cf_da1_n3=0 cf_ca50_n51=2*D/delta_x^2+u/delta_x
/ / c o − e f f i c i e n t o f C−A5 A500 a t node 51 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
cf_ca51_n51=-2*D/ delta_x^2-u/delta_ delta_x^2-u/delta_x x -k cf_da51_n51=0 f o r i=2:50, a(i)=cf_ca1_n2 , end a(51)=cf_ca2_n1 ,c(1)= cf_ca2_n1 f o r i=2:50,c(i)= cf_ca3_n2 cf_ca3_n2 , end d(1)=cf_da1_n1 f o r i=2:51,d(i)=cf_da1_n2 end f o r i=1:51,x(i)=0, end b(1) =cf_ca1_n1 , f o r i=2:51,b(i)= cf_ca2_n2 cf_ca2_n2 , e n d while e1 e 1 > 1 e -6 -6 & e 2 > 1 e -6 -6 d o f o r i = 1 : 5 1 , x 1 ( i ) = x ( i ) , en d , i =1 = 1 , n = 51 51 , B et e t a ( i )= ) = b ( i) i) , Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) end e 1 = a b s ( x ( 4 2 ) - x1 x 1 ( 4 2 ) ) , e 2 = a b s (x(18)-x1(18)) end f o r i=1:51, disp (x(i)) end
75
Scilab code Exa 7.2 boundary value problem in chemical reaction engi-
neering second order 1 clc 2 disp ( ” s o l n o f e g 7.2−− > S e co c o n d O rd r d er e r R e a c t io i o n − 20 par ts ”) 3 e1=1,e2=1 4 u=1,D=10^-4,k=1,C_a_in=1,delta_x=10/20 5 cff_ca2_n1=2*D/delta_x^2+u/delta_x
/ / c o − e f f i c i e n t o f C−A2 a t n od od e 1 6 cff_da1_n1=-(2*u^2/D+2*u/delta_x)*C_a_in
/ / r i g h t hand s i d e co− e f f i c i e n t 7 cff_ca1_n2=D/delta_x^2+u/delta_x 8 cff_ca3_n2=D/delta_x^2
//co− e f f i c i e n t o f C−A3 a t n o d e 2 9 10 11 12 13
cff_da1_n2=0 cff_ca2_n3=cf_ca1_n2 cff_ca4_n3=cf_ca3_n2 cff_da1_n3=0 cff_ca20_n21=2*D/delta_x^2+u/delta_x
/ / c o − e f f i c i e n t o f C−A2 0 a t n od od e 2 1 14 15 16 17 18 19 20 21 22 23
cff_da21_n21=0 f o r i=2:20, a(i)=cff_ca1_n2 , end a(21)=cff_ca2_n1 ,c(1)= cff_ca2_n1 f o r i=2:20,c(i)= cff_ca3_n2 cff_ca3_n2 , end d(1)=cff_da1_n1 f o r i=2:21,d(i)=cff_da1_n2 end f o r i=1:21,x(i)=0,
76
24 e n d 25 while e1 e 1 > 1 e -6 -6 & e 2 > 1 e -6 -6 d o f o r i = 1 : 2 1 , x 1 ( i ) = x ( i ) , en d , 26 cff_ca1_n1=-2*D/ delta_x^2-3*u/delta_x -x1(1) -2*u -2*u / / ma ma i n d i a g o n a l e l e m en en t s ^2/D
d e p e nd nd e n c e o n c o n c . 27 b(1) =cff_ca1_n1 , 28 f o r i=2:21,b(i) =-2*D/delta_x^2-u/ delta_x delta_x -x(i), e n d 29 30 // s o l v i n g by TD TDMA met method hod 31 i =1 = 1 , n = 21 21 , B et e t a ( i) i ) = b (i (i ) , 32 Gamma(i)=d(i)/Beta(i) 33 i1=i+1, 34 f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , 35 Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), 36 end 37 x(n)=Gamma(n) 38 n1=n-i, 39 f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) 40 end 41 e 1 = a b s ( x ( 1 ) - x 1 ( 1) 1) ) , e 2 = a b s (x(21)-x1(21)) 42 e n d 43 f o r i=1:21, disp (x(i)) 44 e n d
77
Chapter 8 CHEMICAL REACTION AND DIFFUSION IN SPHERICAL CATALYST PELLET
Scilab code Exa 8.1 first order reaction
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
clc h e s o l n o f e g 8. 8 .1−− > ” ) ; disp ( ” t he f o r i = 1: 1 : 1 00 00 , x ( i ) = 0 end i te t e r =0 = 0 , e 1 =1 =1 , f = 1 while e1 e1 > 1 e - 6 & f > 1 e -6 -6 d o i te t e r = i te te r +1 + 1 , f o r i=1:100, x1(i)=x(i), en d , f o r i = 2 :1 : 1 0 0 , a ( i ) = 1 - ( 1/ 1/ ( i - 1 ) ) en d , b ( 1) 1 ) = - 6 .0 .0 1 , f o r i = 2 :1 : 1 0 0 , b ( i ) = - 2 . 01 01 en d , c ( 1) 1 ) = 6 , f o r i=2:99,c(i)=1+(1/(i-1)) en d , f o r i=1:99,d(i)=0, en d , d ( 1 0 0) 0 ) = - 10 1 0 0/ 0/ 9 9 , i =1 = 1 , n = 10 10 0 , B et e t a ( i )= ) = b ( i) i) , Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end
78
17 18 19
x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) end e1 = a b s ( x ( 1 ) - x 1 ( 1 ) ) , f = abs ( x ( 1 0 0 ) - x 1 ( 1 0 0 ) ) ,
20 21 22 23 e n d 24 disp (iter, ” t h e s o l u t i o n b y TDMA o f n od od e 7 7 t o 9 9 b y 1 s t o r d er e r r xn xn . i s ” ) ; 25 f o r i=78:100, 26 disp ( x ( i ) ) ; 27 e n d
Scilab code Exa 8.2 second order reaction
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
clc h e s o l n o f e g 8. 8 .2−− > ” ) ; disp ( ” t he f o r i = 1: 1 : 1 00 00 , x ( i ) = 0 end i te t e r =0 = 0 , e 1 =1 =1 , f = 1 k=.1,D=10^-9,r=.01,delta_r=r/10,t1=k*delta_r^2/D while e1 e1 > 1 e - 6 & f > 1 e -6 -6 d o i te t e r = i te te r +1 + 1 , f o r i=1:100, x1(i)=x(i), en d , f o r i = 2 :1 : 1 0 0 , a ( i ) = 1 - ( 1/ 1/ ( i - 1 ) ) en d , b ( 1 ) = - 6 -. - . t 1 * x 1 ( 1) 1 ) , f o r i = 2 :1 :1 0 0 , b ( i ) = -2 -2 - t 1 * x1(i) en d , c ( 1) 1 ) = 6 , f o r i=2:99,c(i)=1+(1/(i-1)) en d , f o r i=1:99,d(i)=0, en d , d ( 1 0 0) 0 ) = - 10 1 0 0/ 0/ 9 9 , i =1 = 1 , n = 10 10 0 , B et e t a ( i )= ) = b ( i) i) , Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end
79
18 19 20
x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) end e1 = a b s ( x ( 1 ) - x 1 ( 1 ) ) , f = abs ( x ( 1 0 0 ) - x 1 ( 1 0 0 ) ) ,
21 22 23 24 e n d 25 disp ( ” t h e s o l u t i o n b y TDMA o f n od o d e 7 7 t o 9 9 b y 2 nd nd o r d er er rx n . i s ” ); 26 f o r i=77:100, 27 disp ( x ( i ) ) ; 28 e n d
Scilab code Exa 8.3 non isothermal condition
1 clc 2 disp ( ” t he h e s o l n o f e g 8. 8 .3−− > ” ) ; 3 f o r i = 1: 1 : 1 00 00 , x ( i ) = 0
// i n i t i a l
values 4 end 5 e 2 = 1 , f 1 =1 =1 , i t e r = 0
//assumed
values 6 k = . 1* 1 * 10 1 0 ^ - 2 , D = 10 1 0 ^ - 9 , r = . 01 01 , d e l ta ta _ r = r / 10 10 0 , t 1 = k * d e l t a_ a_ r / / g i v e n d a ta ta ^2/D 7 / / no no w s o l v i n g t he h e e qn q n s f o r a l l t he h e n od o d es e s a n d t he he n
s i m p l i f y i n g we g et et t h e f o l l o w i n g r e l a t i o n s 8 while e2 e 2 > 1 e - 6 & f 1 > 1 e - 6 d o i te t e r = i te te r +1 + 1 , f o r i=1:100, x1(i)=x(i), 9 en d , f o r i = 2 :1 : 1 0 0 , a ( i ) = 1 - ( 1/ 1/ ( i - 1 ) ) 10 en d , b ( 1) 1 ) = - 6 - t1 t1 * e x p ((1-x1(1))/(2-x1(1))), f o r i = 2 :1 :1 0 0 , b ( i ) = -2 - 2 - t 1 * e x p ((1-x(i))/(2-x(i))) 11 en d , c ( 1) 1 ) = 6 , f o r i=2:99,c(i)=1+(1/(i-1)) 12 en d , f o r i=1:99,d(i)=0, en d , d ( 1 0 0) 0 ) = - 10 1 0 0/ 0/ 9 9 , 13 i =1 = 1 , n = 10 10 0 , B et e t a ( i )= ) = b ( i) i) ,
80
14 15 16 17 18 19 20 21
Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) end e2 = a b s ( x ( 1 ) - x 1 ( 1 ) ) , f1 = a b s ( x ( 1 0 0 ) - x 1 ( 1 0 0 ) ) ,
22 23 24 25 e n d 26 disp ( ” t h e s o l u t i o n b y TDMA o f n od o d e 7 7 t o 1 00 00 by 1 s t o r d er er rx n . i s ” ); 27 f o r i=76:100, 28 disp ( x ( i ) ) ; 29 e n d
Scilab code Exa 8.4 non isothermal condition
1 2 3 4 5 6 7 8
9 10 11
clc h e s o l n o f e g 8. 8 .4−− > ” ) ; disp ( ” t he f o r i = 1: 1 : 1 00 00 , x ( i ) = 0 end i te t e r =0 = 0 , e 1 =1 =1 , f 1 =1 =1 while e1 e 1 > 1 e - 6 & f 1 > 1 e - 6 d o i te t e r = i te te r +1 + 1 , f o r i=1:100, x1(i)=x(i), en d , f o r i = 2 :1 : 1 0 0 , a ( i ) = 1 - ( 1/ 1/ ( i - 1 ) ) en d , b ( 1 ) = - 6 - . 01 0 1 * e x p ((10-10*x1(1))/(11-10*x1(1)) ) , f o r i = 2 :1 : 1 0 0 , b ( i ) = - 2 - . 01 0 1 * e x p ((10-10*x1(i)) /(11-10*x1(i))) en d , c ( 1) 1 ) = 6 , f o r i=2:99,c(i)=1+(1/(i-1)) en d , f o r i=1:99,d(i)=0, en d , d ( 1 0 0) 0 ) = - 10 1 0 0/ 0/ 9 9 , i =1 = 1 , n = 10 10 0 , B et e t a ( i )= ) = b ( i) i) ,
81
12 13 14 15 16 17 18 19
Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) end e1 = a b s ( x ( 1 ) - x 1 ( 1 ) ) , f1 = a b s ( x ( 1 0 0 ) - x 1 ( 1 0 0 ) ) ,
20 21 22 23 e n d 24 disp ( ” t h e s o l u t i o n b y TDMA o f n od od e 7 7 t o 9 9 b y 1 s t o r d er er rx n . i s ” ); 25 f o r i=76:100, 26 disp ( x ( i ) ) ; 27 e n d
82
Chapter 9 ONE DIMENSIONAL TRANSIENT HEAT CONDUCTION
Scilab code Exa 9.1 transient conduction in rectangular slab
1 clc 2 disp ( ” t he he s o l n
o f e g 9. 9 .1−− > T r a n s i en e n t h e at a t c o n d u ct ct i o n i n a R ec e c t an a n g ul ul a r s l a b ”) 3 disp ( ” t h e v a l u e s o f Temp m ea e a su s u re r e d f r o m c e n t r e a t t he he g a p o f . 2 cm a re re ” ) 4 a l p h a = 10 // 1 0 ^ - 5 , d e lt l t a _ t = .1 .1 , d e l ta ta _ x = 1 0^ 0^ - 3 g i v e n d at at a 5 m=alpha*delta_t/(delta_x^2) 6 f o r i = 2: 2 : 4 , a ( i )= )= m
//
s ub−d i a g o n a l 7 end 8 b(1)=(-2*m-1)/2 9 f o r i=2:4,b(i)=-2*m-1
// dia gon al 10 e n d 11 f o r i = 1: 1 : 3 , c ( i )= )= m
//
s u p e r −d i a g o n a l 83
12 e n d 13 f o r i = 1: 1: 4 , x ( i ) = 2 0
//
i n i t i a l tem per atu re 14 e n d 15 f o r t=0.1:.1:3.1, f o r i = 1: 1: 4 , y ( i ) = x ( i ) , e n d
/ /T /TD DMA me th od d(1)=-.5*y(1), d(2)=-y(2) d(3)=-y(3) d(4)= -y(4) -300 i =1 =1 , n = 4 Beta(i)=b(i), Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) 30 end 31 f o r i=1:4, disp (x(i)); 16 17 18 19 20 21 22 23 24 25 26 27 28 29
// s o l u t i o n o f temperature 32 33 34 35
end disp ( ”−−−−−−−−−−−−−−−−−” ) end disp ( ”−−−−−−−−−−END−−−−−−−−− ” ) ;
Scilab code Exa 9.2 transient conduction in rectangular slab
1 clc 2 disp ( ” t he he s o l n
o f e g 9. 9 .2−− > T r a n s i en e n t C on o n du d u ct c t io io n i n R e c ta t a n g u la l a r S l ab ab ” ); 84
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
d e l t a _t _ t = 1 , d e l ta ta _ x = . 05 05 , a l ph ph a = 1 0^ 0^ - 5 t1=alpha*delta_t/delta_x^2 f o r i = 2: 2: 9 , a ( i ) = - t 1 end f o r i=1:9,b(i)=1+2*t1 end f o r i = 1: 1: 8 , c ( i ) = - t 1 end t=1,tf=3000 f o r i = 1: 1: 9 , x ( i ) = 3 00 00 end e1=425, i m e w h e n c e n tr t r e temp i s 4 2 5 K i n s e c s . i s ” ) disp ( ” t im f o r t=1:1:tf, f o r i = 1: 1: 9 , y ( i ) = x ( i ) , e n d d(1)=y(1)+1.7, d(9)=y(9)+2.4, f o r i=2:8,d(i)=y(i) end i =1 =1 , n = 9 Beta(i)=b(i), Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) end i f a bs bs ( x ( 5) 5) - e 1 ) >0 > 0 & a b s (x(5)-e1)<.1 t he h e n d i sp sp ( t ) ; break ; end end h e v a l u e s o f t e m p . a t r eq e q . t im im e i s ” ) ; disp ( ” t he f o r i=1:9, disp (x(i)); e n d
85
Scilab code Exa 9.3 transient conduction in cylinder
1 clc 2 disp ( ” t he h e s o l n o f e g 9. 9 .3−− > T r a n s i en e n t C on o n du d u ct c t io io n i n c y l i n d e r ”); 3 d e l t a _ t = .1 .1 , d e l ta ta _ r = . 00 00 1 , a l ph p h a = 1 0^ 0^ - 5
/ / g i v e n d a ta ta 4 t1=alpha*delta_t/delta_r^2 5 a(2)=.5*t1,a(3)=.75*t1,a(4)=.833*t1
/ / s u b− d i a g o n a l 6 b(1)=-4*t1-1 7 f o r i=2:4,b(i)=-2*t1-1
/ / ma m a in in d i a g o n a l 8 end 9 c(1)=4*t1,c(2)=1.5*t1,c(3)=1.25*t1
/ / s u p e r −d i a g o n a l 10 11 12 13
f o r i = 1: 1: 4 , x ( i ) = 2 0 end T 2 , T2 T 2 & T4 a t t i me me i n t e r v a l o f . 1 s e c i s ” ) disp ( ” T1 , T2 f o r t=.1:.1:2.1, f o r i = 1: 1: 4 , y ( i ) = x ( i ) , e n d
/ /T /TD DMA Me Meth thod od 14 15 16 17 18 19 20 21 22 23 24 25 26
d(4)=-y(4)-7*t1*300/6, f o r i=1:3,d(i)=-y(i) end i =1 =1 , n = 4 Beta(i)=b(i), Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a
86
27 28 29 30
(j) end disp (x(4),x(3),x(2),x(1)); disp ( ”−−−−−−−−−−−−−−−−−−−− ” ) ; end
Scilab code Exa 9.4 transient conduction in sphere
1 clc 2 disp ( ” t he h e s o l n o f e g 9. 9 .4−− > T r an a n s i en e n t c o nd n d u ct c t i on on i n S p h e r e ” ); 3 d e l t a _ t = .1 .1 , d e l ta ta _ r = . 00 00 1 , a l ph p h a = 1 0^ 0^ - 5 4 t1=alpha*delta_t/delta_r^2 5 a(2)=0,a(3)=.5,a(4)=.667 6 b(1)=-7*t1 7 f o r i=2:4,b(i)=-3 8 end 9 c(1)=6,c(2)=2,c(3)=1.5 10 f o r i = 1: 1: 4 , x ( i ) = 2 0 11 e n d 12 disp ( ” T1 , T2 T 2 , T2 T 2 & T4 a t t i me me i n t e r v a l o f . 1 s e c i s ” ) 13 f o r t=.1:.1:1.4, f o r i = 1: 1: 4 , y ( i ) = x ( i ) , e n d 14 d(4)=-y(4)-400, 15 f o r i=1:3,d(i)=-y(i) 16 end 17 i =1 =1 , n = 4 18 Beta(i)=b(i), 19 Gamma(i)=d(i)/Beta(i) 20 i1=i+1, 21 f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , 22 Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), 23 end 24 x(n)=Gamma(n) 25 n1=n-i, 26 f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a
87
27 28 29 30
(j) end disp (x(4),x(3),x(2),x(1)); disp ( ”−−−−−−−−−−−−−−−−−−−” ) ; end
Scilab code Exa 9.5 transient diffusion in sphere
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
clc h e s o l n o f e g 9. 9 .5−− > ” ) ; disp ( ” t he R = . 3 26 26 , D = 3 * 10 10 ^ - 7 d e l ta ta _ t = 1 , d e l ta ta _ r = . 03 0 3 2 6 , c o n c_ c _ i n i = 1 0 / (1 ( 1 . 3 3 * % pi pi * R ^ 3 ) t1=D*delta_t/delta_r^2 n c . o f d ru ru g i s ” ) ; disp (co nc_ini , ” t h e i n i t i a l c o nc f o r i = 2 :1 : 1 0 , a ( i ) = - ( 1 - 1/ 1/ ( i - 1 ) ) end b(1)=591.4 f o r i=2:10,b(i)=3544.5 end c(1)=-1, f o r i = 2 :9 :9 , c ( i ) = - ( 1 + 1/ 1/ ( i - 1 ) ) end f o r i = 1 :1 : 1 0 , x ( i ) = c o n c_ c_ i n i end Co n c . a t c e n t r e a t t =3 =3 h r , 1 2 h r , 2 4 h r , 4 8 h r disp ( ” Co i s ” ); f o r t=1:delta_t:172800, f o r i = 1: 1 : 10 10 , y ( i ) = x ( i ) , e n d d(1)=y(1)*590.4, f o r i=2:10,d(i)=3542.5*y(i) end i =1 = 1 , n = 10 10 Beta(i)=b(i), Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j),
88
27 28 29 30
31 32 33 34 e n d
end x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) end if t = = 1 0 8 00 00 | t = = 4 3 2 00 00 | t = = 8 6 4 00 00 | t = = 1 7 2 8 00 00 then disp (x(6)); end
89
Chapter 10 TWO DIMENSIONAL STEADY AND TRANSIENT HEAT CONDUCTION
Scilab code Exa 10.1 discretization in 2 D space gauss seidel method
1 clc 2 disp ( ” t he h e s o l n o f e g 10 1 0.1−− >2 −D s t e ad a d y h ea ea t c o n d u c t i o n −G a us us s S e i d e l m e t h o d ” ) ; 3 f o r i=1:9,tnew(i)=101,e(i)=1
/ / a ss s s um u m ed ed v a l u e s 4 end 5 t=1e-6 6 / / s i nc e
a l l t h e n o d e s a re re i n t e r i o r n o d es s o d i s c r e t i z e d eqn us e d i s eqn 1 0 . 10
7 while e ( 1 ) > t & e ( 2) 2) > t & e ( 3 ) > t & e ( 4) 4) > t & e ( 5 ) > t & e ( 6 ) > t & e ( 7) 7) > t & e (8 ( 8 ) >t > t & e ( 9) 9) > t d o 8 f o r i = 1: 1 : 9 , t o ld ld ( i ) = t n e w ( i ) , e n d 9 tnew(1)=(told(2)+40+told(4))/4
/ / on on s o l v i n g e q n s f o r v a r i o u s n od o d es e s we g et et , 10 11
tnew(2)=(tnew(1)+told(3)+told(5)+20)/4 tnew(3)=(tnew(2)+told(6)+420)/4
90
12 13 14 15 16 17 18 19 20 21 22 23
tnew(4)=(told(5)+tnew(1)+told(7)+20)/4 tnew(5)=(tnew(2)+told(8)+told(6)+tnew(4))/4 tnew(6)=(tnew(3)+tnew(5)+told(9)+400)/4 tnew(7)=(tnew(4)+told(8)+40)/4 tnew(8)=(tnew(5)+tnew(7)+told(9)+20)/4 tnew(9)=(tnew(6)+420+tnew(8))/4 f o r i=1:9,e(i)= a b s (tnew(i)-told(i)) end
end h e v a l u es e s o f T f ro r o m 1 s t e l em e m en en t t o l a s t disp ( ” t he f o r i=1:9, disp (tnew(i)); end
i s ” );
Scilab code Exa 10.2 discretization in 2 D space gauss seidel method
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
clc h e s o l n o f e g 10 1 0. 2−− > ” ) ; disp ( ” t he f o r i=1:5,tnew(i)=101,e(i)=1 end t=1e-6
/ / s i n c e t h e r e i s n o s o ur u r ce c e t e r m s o we g et e t t he he f o l l o w i n g e qn qn s . while e ( 1) 1) > t & e ( 2 ) > t & e ( 3) 3) > t & e ( 4 ) > t & e (5 (5 ) > t d o f o r i = 1: 1 : 5 , t o ld ld ( i ) = t n e w ( i ) , e n d tnew(1)=(told(2)*2+300)/4 tnew(2)=(tnew(1)+told(3)+300)/4 tnew(3)=(tnew(2)+told(4)+200)/4 tnew(4)=(told(5)+tnew(3)+300)/4 tnew(5)=(2*tnew(4)+300)/4 f o r i=1:5,e(i)= a b s (tnew(1)-told(i)) end end h e v a l u es e s o f T f ro r o m 1 s t e l em e m en en t t o l a s t disp ( ” t he f o r i=1:5, disp (tnew(i)); end
91
i s ” );
Scilab code Exa 10.3 discretization in 2 D space
1 clc 2 disp ( ” t he h e s o l n o f e g 10 1 0. 3−− > Red −B l a c k G au a u ss ss − S e i d e l Method” ) ; 3 f o r j=1:5,tn(j,1)=400, e n d //
d e f i ni n g c o n di t i on s 4 5 6 7 8 9 10
f o r j=2:4,tn(1,j)=20,tn(5,j)=20,tn(j,5)=20, e n d f o r i=1:9,e(i)=1 end f o r i=2:4,j=2:4,tn(i,j)=60 end t1=1e-6 while e ( 1 ) > t 1 & e (2 ( 2 ) > t 1 & e ( 3) 3) > t 1 & e ( 4 ) > t1 t1 & e ( 5 ) > t 1 & e ( 6 ) > t 1& 1 & e ( 7) 7) > t1 t 1 & e ( 8) 8) > t1 t 1 & e ( 9) 9) > t1 t1 do f o r i = 2 : 4 , j =2:4,t(i,j)=tn(i,j), e n d
11 / / u s i n g r ed e d − b l a c k g au a u ss ss − s e i d e l method 12 f o r i=2:4,j=2:4,tn(i,j)=(tn(i+1,j)+tn(i-1,j)+tn(i,j +1)+tn(i,j-1))/4, e n d 13 // / / no w g e t t i n g t h e a b s o l u t e d i f f e r e n c e o f t h e 9
variables 14 e ( 1 ) = a b s (t(2,2)-tn(2,2)),e(2)= a b s (t(2,3)-tn(2,3)),e (3)= a b s (t(2,4)-tn(2,4)),e(4)= a b s ( t ( 3 , 2 ) - t n ( 3 , 2 ) ) , e(5)= a b s (t(3,3)-tn(3,3)),e(6)= a b s (t(3,4)-tn(3,4)) ,e(7)= a b s (t(4,2)-tn(4,2)),e(8)= a b s (t(4,3)-tn(4,3) ),e(9)= a b s ( t ( 4 , 4 ) - t n ( 4 , 4 ) ) , 15 e n d 16 / /n / no w d e f i n i n g p o s i t i o n s o f t h e v a ri ri o us nodes
a c c o rd r d i n g t o r ed e d b l a ck c k c om o m bi b i n at a t io io n 17 R1=t(2,4),R2=t(4,4),R3=t(3,3),R4=t(2,2),R5=t(2,4),B1 =t(4,3),B2=t(3,2),B3=t(3,4),B4=t(2,3) 18 disp (R5,R4,R3,R2,R1, ” t h e v a l u e o f RE RED p o i n t s r e s p e c t i v e l y i s ” ); 19 disp (B4,B3,B2,B1, ” t h e v a l u e o f BLUE p o i n t s
92
r e s p e c t i v e l y i s ” );
Scilab code Exa 10.8 discretization in 2 D space
1 clc 2 disp ( ” t he h e s o l n o f e g 10 1 0. 8−− > G au a u ss ss S e i d e l M e t h o d ” ) ; 3 f o r i=1:9,tnew(i)=101,e(i)=1 //
a ss s s um u m ed ed v a l u e s 4 end 5 t=1e-6 6 while e ( 1 ) > t & e ( 2) 2) > t & e ( 3 ) > t & e ( 4) 4) > t & e ( 5 ) > t & e ( 6 ) > t & e ( 7) 7) > t & e (8 ( 8 ) >t > t & e ( 9) 9) > t d o 7 f o r i = 1: 1 : 9 , t o ld ld ( i ) = t n e w ( i ) , e n d 8 / / u s i n g e qn q n 1 0 . 1 0 f o r t h e i n t e r i o r n o d es e s an and
c o n ve v e c t i v e b o u n d a r y c o n d i t i o n s f o r c o rn rn e r nodes 9 10 11 12
13 14 15 16
17 18 19 20 21 22 23
tnew(1)=(told(2)+50+.5*told(4)+100/3)*3/7 tnew(2)=(tnew(1)+told(3)+told(5)+100)/4 tnew(3)=(tnew(2)+told(6)+600)/4 tnew(4)=(told(5)+.5*tnew(1)+.5*told(7)+100/3) *3/7 tnew(5)=(tnew(2)+told(8)+told(6)+tnew(4))/4 tnew(6)=(tnew(3)+tnew(5)+told(9)+500)/4 tnew(7)=(.5*tnew(4)+.5*told(8)+100/3)*3/4 tnew(8)=(tnew(5)+.5*tnew(7)+.5*told(9)+100/3) *3/7 tnew(9)=(tnew(6)+100/3+.5*tnew(8)+250)*3/7 f o r i=1:9,e(i)= a b s (tnew(i)-told(i)) end
end h e v a l u es e s o f T f ro r o m 1 s t e l em e m en en t t o l a s t disp ( ” t he f o r i=1:9, disp (tnew(i)); end
93
i s ” );
Scilab code Exa 10.9 discretization in 2 D space
1 clc 2 disp ( ” t he h e s o l n o f e g 10 1 0. 9−− > S te t e ad a d y s t a t e h ea ea t c o n du d u c t io i o n i n L−s h a p ed e d b od od y ” ) ; 3 f o r i=1:9,tnew(i)=101,e(i)=1 //assumed
values 4 end 5 t=1e-6 6 while e ( 1 ) > t & e ( 2) 2) > t & e ( 3 ) > t & e ( 4) 4) > t & e ( 5 ) > t & e ( 6 ) > t & e ( 7) 7) > t & e (8 ( 8 ) >t > t & e ( 9) 9) > t d o 7 f o r i = 1: 1 : 9 , t o ld ld ( i ) = t n e w ( i ) , e n d 8 / / u si s i ng ng t h e b a s i c d i s c r e t i z a t i o n eqn . f o r a l l
t he h e n od o d es e s s i n c e t he h e b o un u n d ar a r y c o n d i t i o n s v ar ar y f o r e aacc h p o in in t 9 10
11 12 13 14
15 16 17 18 19 20 21 22 23
tnew(1)=(told(2)+1.25+told(4))/2.05 tnew(2)=(.5*tnew(1)+.5*told(3)+told(5)+1.25) /2.05 tnew(3)=(.5*tnew(2)+.5*told(6)+1.25)/1.05 tnew(4)=(told(5)+.5*tnew(1)+45)/2 tnew(5)=(tnew(2)+told(8)+90+tnew(4))/4 tnew(6)=(.5*tnew(3)+tnew(5)+.5*told(7)+91.25) /3.05 tnew(7)=(.5*tnew(6)+.5*told(8)+91.25)/2.05 tnew(8)=(91.25+.5*tnew(7)+.5*told(9))/2.05 tnew(9)=(47.125+.5*tnew(8))/1.025 f o r i=1:9,e(i)= a b s (tnew(i)-told(i)) end
end h e v a l u es e s o f T f ro r o m 1 s t e l em e m en en t t o l a s t disp ( ” t he f o r i=1:9, disp (told(i)); end
94
i s ” );
Scilab code Exa 10.10 ADI method
1 clc 2 disp ( ” t h e s o l n o f e g 1 0 . 1 0 −− > A l t er er n a ti n g D i r e ct i o n I m p l i c i t Method ” ) ; 3 nmax=25 4 a ( 2 ) = 1 , a ( 3) // 3) = 1 ,
d e f i ni n g v a r i a b l e s 5 b(1)=-4,b(2)=-4,b(3)=-4 6 c(1)=1,c(2)=1 7 t(1,2)=20,t(1,3)=20,t(1,4)=20,t(2,1)=20,t(3,1)=20,t (4,1)=20,t(5,2)=20,t(5,3)=20,t(5,4)=20,t(2,5) =400,t(3,5)=400,t(4,5)=400 8 tstar(1,2)=20,tstar(1,3)=20,tstar(1,4)=20,tstar(2,1) =20,tstar(3,1)=20,tstar(4,1)=20,tstar(5,2)=20, tstar(5,3)=20,tstar(5,4)=20,tstar(2,5)=400,tstar (3,5)=400,tstar(4,5)=400 9 f o r i=2:4, f o r j = 2 :4 : 4 t ( i , j ) = 20 20 10 end 11 e n d 12 // s o l v i n g by TD TDMA Me Meth thod od 13 f o r nn=1:nmax, f o r jj=2:4,d(1)=-t(1,jj)-t(2,jj+1)tstar(2,jj-1), 14 d(2)=-t(3,jj+1)-tstar(3,jj-1), 15 d(3)=-t(5,jj)-t(4,jj+1)-tstar(4,jj-1) 16 i=1,n=3 17 Beta(i)=b(i), 18 Gamma(i)=d(i)/Beta(i) 19 i1=i+1, 20 f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , 21 Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), 22 end 23 x(n)=Gamma(n) 24 n1=n-i,
95
25
f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) end tstar(2,jj)=x(1) tstar(3,jj)=x(2) tstar(4,jj)=x(3)
26 27 28 29 30 e n d 31 f o r ii=2:4,d(1)=-t(ii,1)-tstar(ii+1,2)-tstar(ii-1,2) , 32 d(2)=-tstar(ii+1,3)-tstar(ii-1,3) 33 d(3)=-t(ii,5)-tstar(ii+1,4)-tstar(ii-1,4) 34 i=1,n=3 35 Beta(i)=b(i), 36 Gamma(i)=d(i)/Beta(i) 37 i1=i+1, 38 f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , 39 Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), 40 end 41 x(n)=Gamma(n) 42 n1=n-i, 43 f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) 44 end 45 t(ii,2)=x(1) 46 t(ii,3)=x(2) 47 t(ii,4)=x(3) 48 e n d 49 e n d 50 disp ( ” t h e s o l n b y ADI m e t h o d i s ” ) ; 51 disp (t(2,4),t(2,3),t(2,2)); 52 disp ( ”−−−−−−−−−−−−−− ” ) ; 53 disp (t(3,4),t(3,3),t(3,2)); 54 disp ( ”−−−−−−−−−−−−−− ” ) ; 55 disp (t(4,4),t(4,3),t(4,2));
96
Scilab code Exa 10.11 ADI method for transient heat conduction
1 2 3 4 5 6 7 8 9 10 11 12 13 14
15 16 17 18 19 20 21 22 23 24 25 26 27 28
clc disp ( ” t h e s o l n o f e g 1 0 . 1 1 −− > ” ) ; f o r k = 2 :1 :1 0 , a ( k )= ) = -1 - 1 , b ( k ) =2 =2 .4 .4 , c ( k ) = -1 -1 end alpha=1,delta_t=.05,delta_x=.1 m=alpha*delta_t/delta_x^2 b(1)=2.4 c(1)=-2 f o r k = 1: 1: 1 1 , t ( 11 11 , k ) = 4 0 0 , t s ta t a r ( 11 1 1 , k ) = 4 0 0 , t ( k , 1 1) 1) =400,tstar(k,11)=400 end f o r i=1:10, f o r j = 1: 1: 1 0 , t ( i , j ) =0 =0 , t s ta ta r ( i , j ) = 0 end end f o r tm=.05:.05:.5, f o r jj=1:10, if jj==1 t h en en f or or i i =1:10, =1:10, d(ii)=2*t(ii ,2) -1.6 -1.6*t(ii *t(ii ,1), en d ,d(10)=d(10) +400, e ls l s e f or or i i = 1: 1 : 10 1 0 , d ( i i ) = t ( ii ii , j j + 1 ) + t ( ii ii , j j -1) -1.6*t(ii, jj), en d ,d(10)=d(10)+400, en d , i =1 = 1 , n = 10 10 Beta(i)=b(i), Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) en d , f o r c o u n t = 1 :1 :1 0 , t s ta t a r ( c o un un t , j j ) = x ( c o u n t ) , end end f o r ii=1:10, if ii==1 t h en en f or or jj=1:10,d(jj)=2*tstar(ii +1,1) -1.6*tstar(ii -1.6*tstar(ii ,1), en d , d ( 1 0 ) = d ( 10 10 ) +400, e ls l s e f or o r jj=1:10,d(jj)=tstar(ii-1,jj
97
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
)+tstar(ii+1,jj)-1.6*tstar(ii,jj), en d , d (10)=d(10)+400, e n d i =1 = 1 , n = 10 10 Beta(i)=b(i), Gamma(i)=d(i)/Beta(i) i1=i+1, f o r j = i 1 : n , B e ta ta ( j ) = b ( j ) - a ( j )* ) * c ( j - 1) 1 ) / B e t a ( j - 1) 1) , Gamma(j)=(d(j)-a(j)*Gamma(j-1))/Beta(j), end x(n)=Gamma(n) n1=n-i, f o r k = 1 : n1 n1 , j = n - k , x ( j )= ) = G a m ma ma ( j ) - c ( j ) * x ( j + 1) 1) / B e t a (j) en d , f o r c o u n t = 1 :1 :1 0 , t ( ii i i , c o un u n t ) = x ( c o un un t ) , e n d end end disp ( ” t h e s o l n b y ADI i s ” ) ; f o r i=1:10,j=1:10, disp (t(i,j)); end h e t a b l e i s a c t u a l l y i n t er e r c h a n g e d row & disp ( ” t he
c o l u m n . . h o r i z o n t a l l y a r e t he h e e l e me m e n t s o f t he he column” ) ;
98