Introduction to Numerical Methods in Chemical Engineering_P. Ahuja

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